diff --git a/CMakeLists.txt b/CMakeLists.txt
index 904d731..08b93b3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,17 +1,15 @@
 cmake_minimum_required(VERSION 3.10)
 project(${SKBUILD_PROJECT_NAME} VERSION ${SKBUILD_PROJECT_VERSION})
 
-if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES)
-    set(CMAKE_CUDA_ARCHITECTURES "native")
+find_package(Ceres REQUIRED)
+if(NOT TARGET Ceres::ceres)
+    # Older Ceres versions don't come with an imported interface target.
+    add_library(Ceres::ceres INTERFACE IMPORTED)
+    target_include_directories(
+        Ceres::ceres INTERFACE ${CERES_INCLUDE_DIRS})
+    target_link_libraries(
+        Ceres::ceres INTERFACE ${CERES_LIBRARIES})
 endif()
-find_package(COLMAP REQUIRED)
-if(${COLMAP_VERSION} VERSION_LESS "3.9.0")
-    message( SEND_ERROR "COLMAP version >= 3.9.0 required, found ${COLMAP_VERSION}." )
-endif()
-if(${CERES_VERSION} VERSION_LESS "2.1.0")
-    message( SEND_ERROR "Ceres version >= 2.1 required, found ${CERES_VERSION}." )
-endif()
-
 
 find_package(Python REQUIRED COMPONENTS Interpreter Development.Module)
 
@@ -19,6 +17,6 @@ find_package(pybind11 2.11.1 REQUIRED)
 
 pybind11_add_module(pyceres _pyceres/bindings.cc)
 target_include_directories(pyceres PRIVATE ${PROJECT_SOURCE_DIR})
-target_link_libraries(pyceres PRIVATE colmap::colmap glog::glog Ceres::ceres)
+target_link_libraries(pyceres PRIVATE glog::glog Ceres::ceres)
 target_compile_definitions(pyceres PRIVATE VERSION_INFO="${PROJECT_VERSION}")
 install(TARGETS pyceres LIBRARY DESTINATION .)
diff --git a/Dockerfile b/Dockerfile
index d57da3d..5f56853 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -2,9 +2,6 @@ ARG UBUNTU_VERSION=22.04
 ARG NVIDIA_CUDA_VERSION=12.3.1
 FROM nvidia/cuda:${NVIDIA_CUDA_VERSION}-devel-ubuntu${UBUNTU_VERSION} as builder
 
-ARG COLMAP_VERSION=3.9.1
-ARG CUDA_ARCHITECTURES=70
-ENV CUDA_ARCHITECTURES=${CUDA_ARCHITECTURES}
 ENV QT_XCB_GL_INTEGRATION=xcb_egl
 
 # Prevent stop building ubuntu at time zone selection.
@@ -17,21 +14,12 @@ RUN apt-get update && \
         cmake \
         ninja-build \
         build-essential \
-        libboost-program-options-dev \
-        libboost-filesystem-dev \
-        libboost-graph-dev \
-        libboost-system-dev \
         libeigen3-dev \
-        libflann-dev \
-        libfreeimage-dev \
-        libmetis-dev \
         libgoogle-glog-dev \
+        libgflags-dev \
         libgtest-dev \
-        libsqlite3-dev \
-        libglew-dev \
-        qtbase5-dev \
-        libqt5opengl5-dev \
-        libcgal-dev \
+        libatlas-base-dev \
+        libsuitesparse-dev \
         python-is-python3 \
         python3-minimal \
         python3-pip \
@@ -47,16 +35,7 @@ RUN apt-get install -y --no-install-recommends --no-install-suggests wget && \
     cmake ../ceres-solver-2.1.0 -GNinja && \
     ninja install
 
-# Install Colmap.
-RUN wget "https://github.com/colmap/colmap/archive/refs/tags/${COLMAP_VERSION}.tar.gz" -O colmap-${COLMAP_VERSION}.tar.gz && \
-    tar zxvf colmap-${COLMAP_VERSION}.tar.gz && \
-    mkdir colmap-build && \
-    cd colmap-build && \
-    cmake ../colmap-${COLMAP_VERSION} -GNinja -DCMAKE_CUDA_ARCHITECTURES=${CUDA_ARCHITECTURES} && \
-    ninja install
-
-
 # Build pyceres.
 COPY . /pyceres
 WORKDIR /pyceres
-RUN pip install . -vv --config-settings=cmake.define.CMAKE_CUDA_ARCHITECTURES=${CUDA_ARCHITECTURES}
+RUN pip install . -vv
diff --git a/README.md b/README.md
index dfaa800..d004264 100644
--- a/README.md
+++ b/README.md
@@ -4,9 +4,7 @@ This repository provides minimal Python bindings for the [Ceres Solver](http://c
 
 ## Installation
 
-1. Install [COLMAP 3.9.1](https://colmap.github.io/)
-
-2. Clone the repository and build the package:
+Clone the repository and build the package:
 
 ```sh
 git clone https://github.com/cvg/pyceres.git
@@ -14,36 +12,22 @@ cd pyceres
 python -m pip install .
 ```
 
-### Docker image
-
 Alternatively, you can build the Docker image:
 
 ```sh
-export COLMAP_VERSION=3.9.1
-export CUDA_ARCHITECTURES=70
-docker build -t pyceres \
-    --build-arg COLMAP_VERSION=${COLMAP_VERSION}  \
-    --build-arg CUDA_ARCHITECTURES=${CUDA_ARCHITECTURES}  \
-    -f Dockerfile .
+docker build -t pyceres -f Dockerfile .
 ```
 
 ## Factor graph optimization
 
-For now we support the following cost functions, defined in `_pyceres/factors/`:
-- camera reprojection error (with fixed or variable pose)
-- rig reprojection error (with fixed or variable rig extrinsics)
-- relative pose prior
-- absolute pose prior
-
-All factors support basic observation covariances. Reprojection error costs rely on camera models defined in COLMAP. Absolute poses are represented as quaternions and are expressed in the sensor frame, so are pose residuals, which use the right-hand convention as in the [GTSAM library](https://github.com/borglab/gtsam).
-
-## Examples
-See the Jupyter notebooks in `examples/`.
+Factors may be defined in Python (see [`examples/test_python_cost.py`](./examples/test_python_cost.py)) or in C++ with associated Python bindings.
+[PyCOLMAP](https://github.com/colmap/colmap/tree/main/pycolmap) provides the following cost functions in `pycolmap.cost_functions`:
+- reprojection error for different camera models, with fixed or variable pose and 3D points
+- reprojection error for multi-camera rigs, with fixed or variable rig extrinsics
+- error of absolute and relative poses
+- Sampson error for epipolar geometry
 
-## TODO
-- [ ] Define a clean interface for covariances, like in GTSAM
-- [ ] Add bindings for Ceres covariance estimation
-- [ ] Proper benchmark against GTSAM
+See [`examples/`](./examples/) to use these factors.
 
 ## Credits
 Pyceres was inspired by the work of Nikolaus Mitchell for [ceres_python_bindings](https://github.com/Edwinem/ceres_python_bindings) and is maintained by [Philipp Lindenberger](https://github.com/Phil26AT) and [Paul-Edouard Sarlin](https://psarlin.com/).
diff --git a/_pyceres/factors/bindings.h b/_pyceres/factors/bindings.h
index 3f2809f..31e25b7 100644
--- a/_pyceres/factors/bindings.h
+++ b/_pyceres/factors/bindings.h
@@ -1,85 +1,24 @@
 #pragma once
 
-#include "_pyceres/factors/bundle.h"
-#include "_pyceres/factors/common.h"
-#include "_pyceres/factors/pose_graph.h"
 #include "_pyceres/helpers.h"
 
+#include <ceres/ceres.h>
+#include <ceres/normal_prior.h>
 #include <pybind11/eigen.h>
 #include <pybind11/pybind11.h>
-#include <pybind11/stl.h>
 
 namespace py = pybind11;
-using namespace colmap;
 
 void bind_factors(py::module& m) {
-  m.def("ReprojErrorCost",
-        &CreateCostFunction<ReprojErrorCostFunction, const Eigen::Vector2d&>,
-        py::arg("camera_model_id"),
-        py::arg("point2D"));
-  m.def("ReprojErrorCost",
-        &CreateCostFunction<ReprojErrorCostFunctionWithNoise,
-                            const Eigen::Vector2d&,
-                            const double>,
-        py::arg("camera_model_id"),
-        py::arg("point2D"),
-        py::arg("stddev"));
-  m.def("ReprojErrorCost",
-        &CreateCostFunction<ReprojErrorConstantPoseCostFunction,
-                            const Rigid3d&,
-                            const Eigen::Vector2d&>,
-        py::arg("camera_model_id"),
-        py::arg("cam_from_world"),
-        py::arg("point2D"));
-  m.def("ReprojErrorCost",
-        &CreateCostFunction<ReprojErrorConstantPoseCostFunctionWithNoise,
-                            const Rigid3d&,
-                            const Eigen::Vector2d&,
-                            const double>,
-        py::arg("camera_model_id"),
-        py::arg("cam_from_world"),
-        py::arg("point2D"),
-        py::arg("stddev"));
-
-  m.def("RigReprojErrorCost",
-        &CreateCostFunction<RigReprojErrorCostFunction, const Eigen::Vector2d&>,
-        py::arg("camera_model_id"),
-        py::arg("point2D"));
-  m.def("RigReprojErrorCost",
-        &CreateCostFunction<RigReprojErrorCostFunctionWithNoise,
-                            const Eigen::Vector2d&,
-                            const double>,
-        py::arg("camera_model_id"),
-        py::arg("point2D"),
-        py::arg("stddev"));
-  m.def("RigReprojErrorCost",
-        &CreateCostFunction<RigReprojErrorConstantRigCostFunction,
-                            const Rigid3d&,
-                            const Eigen::Vector2d&>,
-        py::arg("camera_model_id"),
-        py::arg("cam_from_rig"),
-        py::arg("point2D"));
-  m.def("RigReprojErrorCost",
-        &CreateCostFunction<RigReprojErrorConstantRigCostFunctionWithNoise,
-                            const Rigid3d&,
-                            const Eigen::Vector2d&,
-                            const double>,
-        py::arg("camera_model_id"),
-        py::arg("cam_from_rig"),
-        py::arg("point2D"),
-        py::arg("stddev"));
-
-  m.def("PoseGraphRelativeCost",
-        &PoseGraphRelativeCost::Create,
-        py::arg("j_from_i"),
-        py::arg("covariance"));
-  m.def("PoseGraphAbsoluteCost",
-        &PoseGraphAbsoluteCost::Create,
-        py::arg("cam_from_world"),
-        py::arg("covariance"));
-
-  m.def("NormalPrior",
-        &CreateNormalPrior,
-        py::arg("mean"),
-        py::arg("covariance"));
+  m.def(
+      "NormalPrior",
+      [](const Eigen::VectorXd& mean,
+         const Eigen::Matrix<double, -1, -1>& covariance) {
+        THROW_CHECK_EQ(covariance.cols(), mean.size());
+        THROW_CHECK_EQ(covariance.cols(), covariance.rows());
+        return new ceres::NormalPrior(covariance.inverse().llt().matrixL(),
+                                      mean);
+      },
+      py::arg("mean"),
+      py::arg("covariance"));
 }
diff --git a/_pyceres/factors/bundle.h b/_pyceres/factors/bundle.h
deleted file mode 100644
index 0cee120..0000000
--- a/_pyceres/factors/bundle.h
+++ /dev/null
@@ -1,128 +0,0 @@
-#pragma once
-
-#include "_pyceres/log_exceptions.h"
-
-#include <ceres/ceres.h>
-#include <ceres/rotation.h>
-#include <colmap/estimators/cost_functions.h>
-#include <colmap/geometry/rigid3.h>
-#include <colmap/sensor/models.h>
-
-using namespace colmap;
-
-class LinearCostFunction : public ceres::CostFunction {
- public:
-  LinearCostFunction(const double s) : s_(s) {
-    set_num_residuals(1);
-    mutable_parameter_block_sizes()->push_back(1);
-  }
-
-  bool Evaluate(double const* const* parameters,
-                double* residuals,
-                double** jacobians) const final {
-    *residuals = **parameters * s_;
-    if (jacobians && *jacobians) {
-      **jacobians = s_;
-    }
-    return true;
-  }
-
- private:
-  const double s_;
-};
-
-template <class CostFunction, typename... Args>
-class CostFunctionIsotropicNoise {
- public:
-  static ceres::CostFunction* Create(Args&&... args, const double stddev) {
-    THROW_CHECK_GE(stddev, 0.0);
-    ceres::CostFunction* cost_function =
-        CostFunction::Create(std::forward<Args>(args)...);
-    std::vector<ceres::CostFunction*> conditioners(
-        cost_function->num_residuals(), new LinearCostFunction(1.0 / stddev));
-    return new ceres::ConditionedCostFunction(
-        cost_function, conditioners, ceres::TAKE_OWNERSHIP);
-  }
-};
-
-template <typename CameraModel>
-class RigReprojErrorConstantRigCostFunction
-    : public ReprojErrorCostFunction<CameraModel> {
-  using Parent = ReprojErrorCostFunction<CameraModel>;
-
- public:
-  explicit RigReprojErrorConstantRigCostFunction(const Rigid3d& cam_from_rig,
-                                                 const Eigen::Vector2d& point2D)
-      : Parent(point2D), cam_from_rig_(cam_from_rig) {}
-
-  static ceres::CostFunction* Create(const Rigid3d& cam_from_rig,
-                                     const Eigen::Vector2d& point2D) {
-    return (new ceres::AutoDiffCostFunction<
-            RigReprojErrorConstantRigCostFunction<CameraModel>,
-            2,
-            4,
-            3,
-            3,
-            CameraModel::num_params>(
-        new RigReprojErrorConstantRigCostFunction(cam_from_rig, point2D)));
-  }
-
-  template <typename T>
-  bool operator()(const T* const rig_from_world_rotation,
-                  const T* const rig_from_world_translation,
-                  const T* const point3D,
-                  const T* const camera_params,
-                  T* residuals) const {
-    const Eigen::Quaternion<T> cam_from_world_rotation =
-        cam_from_rig_.rotation.cast<T>() *
-        EigenQuaternionMap<T>(rig_from_world_rotation);
-    const Eigen::Matrix<T, 3, 1> cam_from_world_translation =
-        cam_from_rig_.rotation.cast<T>() *
-            EigenVector3Map<T>(rig_from_world_translation) +
-        cam_from_rig_.translation.cast<T>();
-    return Parent::operator()(cam_from_world_rotation.coeffs().data(),
-                              cam_from_world_translation.data(),
-                              point3D,
-                              camera_params,
-                              residuals);
-  }
-
- private:
-  const Rigid3d& cam_from_rig_;
-};
-
-template <typename CameraModel>
-using ReprojErrorCostFunctionWithNoise =
-    CostFunctionIsotropicNoise<ReprojErrorCostFunction<CameraModel>,
-                               const Eigen::Vector2d&>;
-
-template <typename CameraModel>
-using ReprojErrorConstantPoseCostFunctionWithNoise =
-    CostFunctionIsotropicNoise<ReprojErrorConstantPoseCostFunction<CameraModel>,
-                               const Rigid3d&,
-                               const Eigen::Vector2d&>;
-
-template <typename CameraModel>
-using RigReprojErrorCostFunctionWithNoise =
-    CostFunctionIsotropicNoise<RigReprojErrorCostFunction<CameraModel>,
-                               const Eigen::Vector2d&>;
-
-template <typename CameraModel>
-using RigReprojErrorConstantRigCostFunctionWithNoise =
-    CostFunctionIsotropicNoise<
-        RigReprojErrorConstantRigCostFunction<CameraModel>,
-        const Rigid3d&,
-        const Eigen::Vector2d&>;
-
-template <template <typename> class CostFunction, typename... Args>
-ceres::CostFunction* CreateCostFunction(const CameraModelId camera_model_id,
-                                        Args&&... args) {
-  switch (camera_model_id) {
-#define CAMERA_MODEL_CASE(CameraModel)                                     \
-  case CameraModel::model_id:                                              \
-    return CostFunction<CameraModel>::Create(std::forward<Args>(args)...); \
-    break;
-    CAMERA_MODEL_SWITCH_CASES
-#undef CAMERA_MODEL_CASE
-  }
-}
diff --git a/_pyceres/factors/common.h b/_pyceres/factors/common.h
deleted file mode 100644
index fbef0e0..0000000
--- a/_pyceres/factors/common.h
+++ /dev/null
@@ -1,14 +0,0 @@
-#pragma once
-
-#include "_pyceres/helpers.h"
-
-#include <ceres/ceres.h>
-#include <ceres/normal_prior.h>
-
-ceres::CostFunction* CreateNormalPrior(
-    const Eigen::VectorXd& mean,
-    const Eigen::Matrix<double, -1, -1>& covariance) {
-  THROW_CHECK_EQ(covariance.cols(), mean.size());
-  THROW_CHECK_EQ(covariance.cols(), covariance.rows());
-  return new ceres::NormalPrior(covariance.inverse().llt().matrixL(), mean);
-}
diff --git a/_pyceres/factors/pose_graph.h b/_pyceres/factors/pose_graph.h
deleted file mode 100644
index 8d2707a..0000000
--- a/_pyceres/factors/pose_graph.h
+++ /dev/null
@@ -1,107 +0,0 @@
-#pragma once
-
-#include <Eigen/Core>
-#include <ceres/ceres.h>
-#include <ceres/rotation.h>
-#include <colmap/geometry/rigid3.h>
-
-using namespace colmap;
-
-using Matrix6d = Eigen::Matrix<double, 6, 6>;
-
-template <typename T>
-inline void EigenQuaternionToAngleAxis(const T* eigen_quaternion,
-                                       T* angle_axis) {
-  const T quaternion[4] = {eigen_quaternion[3],
-                           eigen_quaternion[0],
-                           eigen_quaternion[1],
-                           eigen_quaternion[2]};
-  ceres::QuaternionToAngleAxis(quaternion, angle_axis);
-}
-
-/* Adapted from:
- * https://github.com/ceres-solver/ceres-solver/blob/master/examples/slam/pose_graph_3d/pose_graph_3d_error_term.h
- * https://github.com/ceres-solver/ceres-solver/blob/master/internal/ceres/autodiff_benchmarks/relative_pose_error.h
- * Use right-hand multiplication like in GTSAM.
- * The covariance is expressed in the frame j.
- */
-struct PoseGraphRelativeCost {
- public:
-  explicit PoseGraphRelativeCost(const Rigid3d& j_from_i,
-                                 const Matrix6d covariance)
-      : j_from_i_(j_from_i),
-        sqrt_information_(covariance.inverse().llt().matrixL()) {}
-
-  static ceres::CostFunction* Create(const Rigid3d& j_from_i,
-                                     const Matrix6d covariance) {
-    return (
-        new ceres::AutoDiffCostFunction<PoseGraphRelativeCost, 6, 4, 3, 4, 3>(
-            new PoseGraphRelativeCost(j_from_i, covariance)));
-  }
-
-  template <typename T>
-  bool operator()(const T* const i_from_world_q,
-                  const T* const i_from_world_t,
-                  const T* const j_from_world_q,
-                  const T* const j_from_world_t,
-                  T* residuals_ptr) const {
-    const Eigen::Quaternion<T> i_from_j_q =
-        EigenQuaternionMap<T>(i_from_world_q) *
-        EigenQuaternionMap<T>(j_from_world_q).inverse();
-    Eigen::Matrix<T, 3, 1> i_from_j_t =
-        EigenVector3Map<T>(i_from_world_t) -
-        i_from_j_q * EigenVector3Map<T>(j_from_world_t);
-
-    const Eigen::Quaternion<T> mes_from_est_q =
-        j_from_i_.rotation.cast<T>() * i_from_j_q;
-    EigenQuaternionToAngleAxis(mes_from_est_q.coeffs().data(), residuals_ptr);
-
-    Eigen::Map<Eigen::Matrix<T, 3, 1>> mes_from_est_t(residuals_ptr + 3);
-    mes_from_est_t = j_from_i_.translation.cast<T>() +
-                     j_from_i_.rotation.cast<T>() * i_from_j_t;
-
-    Eigen::Map<Eigen::Matrix<T, 6, 1>> residuals(residuals_ptr);
-    residuals.applyOnTheLeft(sqrt_information_.template cast<T>());
-    return true;
-  }
-
- private:
-  const Rigid3d& j_from_i_;
-  const Matrix6d sqrt_information_;
-};
-
-struct PoseGraphAbsoluteCost {
- public:
-  explicit PoseGraphAbsoluteCost(const Rigid3d& cam_from_world,
-                                 const Matrix6d& covariance)
-      : cam_from_world_(cam_from_world),
-        sqrt_information_(covariance.inverse().llt().matrixL()) {}
-
-  static ceres::CostFunction* Create(const Rigid3d& cam_from_world,
-                                     const Matrix6d covariance) {
-    return (new ceres::AutoDiffCostFunction<PoseGraphAbsoluteCost, 6, 4, 3>(
-        new PoseGraphAbsoluteCost(cam_from_world, covariance)));
-  }
-
-  template <typename T>
-  bool operator()(const T* const cam_from_world_q,
-                  const T* const cam_from_world_t,
-                  T* residuals_ptr) const {
-    const Eigen::Quaternion<T> mes_from_est_q =
-        cam_from_world_.rotation.cast<T>() *
-        EigenQuaternionMap<T>(cam_from_world_q).inverse();
-    EigenQuaternionToAngleAxis(mes_from_est_q.coeffs().data(), residuals_ptr);
-
-    Eigen::Map<Eigen::Matrix<T, 3, 1>> mes_from_est_t(residuals_ptr + 3);
-    mes_from_est_t = cam_from_world_.translation.cast<T>() -
-                     mes_from_est_q * EigenVector3Map<T>(cam_from_world_t);
-
-    Eigen::Map<Eigen::Matrix<T, 6, 1>> residuals(residuals_ptr);
-    residuals.applyOnTheLeft(sqrt_information_.template cast<T>());
-    return true;
-  }
-
- private:
-  const Rigid3d& cam_from_world_;
-  const Matrix6d sqrt_information_;
-};
diff --git a/_pyceres/log_exceptions.h b/_pyceres/log_exceptions.h
index 454df1e..eb1f75d 100644
--- a/_pyceres/log_exceptions.h
+++ b/_pyceres/log_exceptions.h
@@ -4,7 +4,6 @@
 #include <sstream>
 #include <string>
 
-#include <colmap/util/misc.h>
 #include <pybind11/pybind11.h>
 
 namespace py = pybind11;
@@ -124,11 +123,3 @@ void __ThrowCheckOpImpl(const char* file,
 #define THROW_CHECK_LT(val1, val2) THROW_CHECK_OP(_LT, <, val1, val2)
 #define THROW_CHECK_GE(val1, val2) THROW_CHECK_OP(_GE, >=, val1, val2)
 #define THROW_CHECK_GT(val1, val2) THROW_CHECK_OP(_GT, >, val1, val2)
-
-#define THROW_CHECK_FILE_EXISTS(path) \
-  THROW_CHECK_MSG(ExistsFile(path),   \
-                  std::string("File ") + path + " does not exist.");
-
-#define THROW_CHECK_DIR_EXISTS(path) \
-  THROW_CHECK_MSG(ExistsDir(path),   \
-                  std::string("Directory ") + path + " does not exist.");