Skip to content

Camera model fixes and unit test overhaul #2202

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
servantftransperfect merged 11 commits into from
Oct 16, 2026
Merged

Camera model fixes and unit test overhaul #2202

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
ef2d3c4
[camera] Use the largest image side when setting the focal length in mm
Oct 12, 2026
a1ce2c7
[camera] Rescale the undistortion with the intrinsic
Oct 12, 2026
b291a66
[camera] Equidistant: compare and rescale the image circle
Oct 12, 2026
e73b9e9
[camera] Deep copy the distortion when assigning intrinsics
Oct 12, 2026
9ddd761
[camera] Fix the Equirectangular Jacobians
Oct 12, 2026
44f2f27
[camera] Exact Jacobians of removeDistortion for all distortion models
Oct 12, 2026
bc5b3a6
[camera] Compose Distortion and Undistortion, differentiate the Disto…
Oct 12, 2026
ec7dc0f
[sfm] Intrinsics cost functions: compute residuals in undistorted space
Oct 12, 2026
658fb9b
[camera] Rework the camera unit tests
Oct 12, 2026
5e82364
[sfm] Add unit tests for the intrinsics cost functions with Undistortion
Oct 12, 2026
b34db54
[camera] Remove the Pinhole(w, h, K) constructor
Oct 14, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 7 additions & 7 deletions src/aliceVision/camera/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,13 +52,13 @@ alicevision_add_library(aliceVision_camera
)

# Unit tests
alicevision_add_test(distortion_test.cpp NAME "camera_distortionRadial" LINKS aliceVision_camera)
alicevision_add_test(pinholeBrown_test.cpp NAME "camera_pinholeBrown" LINKS aliceVision_camera)
alicevision_add_test(pinholeFisheye_test.cpp NAME "camera_pinholeFisheye" LINKS aliceVision_camera)
alicevision_add_test(pinholeFisheye1_test.cpp NAME "camera_pinholeFisheye1" LINKS aliceVision_camera)
alicevision_add_test(pinholeRadial_test.cpp NAME "camera_pinholeRadial" LINKS aliceVision_camera)
alicevision_add_test(equidistant_test.cpp NAME "camera_equidistant" LINKS aliceVision_camera)
alicevision_add_test(undistortion_test.cpp NAME "camera_undistortion" LINKS aliceVision_camera)
alicevision_add_test(distortion_test.cpp NAME "camera_distortionRadial" LINKS aliceVision_camera)
alicevision_add_test(undistortion_test.cpp NAME "camera_undistortion" LINKS aliceVision_camera)
alicevision_add_test(intrinsicProjection_test.cpp NAME "camera_intrinsicProjection" LINKS aliceVision_camera)
alicevision_add_test(intrinsicDerivatives_test.cpp NAME "camera_intrinsicDerivatives" LINKS aliceVision_camera)
alicevision_add_test(intrinsicState_test.cpp NAME "camera_intrinsicState" LINKS aliceVision_camera)
alicevision_add_test(cameraCommon_test.cpp NAME "camera_cameraCommon" LINKS aliceVision_camera)
alicevision_add_test(cameraUndistortImage_test.cpp NAME "camera_undistortImage" LINKS aliceVision_camera aliceVision_image)


# SWIG Binding
Expand Down
15 changes: 8 additions & 7 deletions src/aliceVision/camera/DistortionBrown.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
// You can obtain one at https://mozilla.org/MPL/2.0/.

#include "DistortionBrown.hpp"
#include <aliceVision/system/Logger.hpp>
#include <vector>

namespace aliceVision {
Expand Down Expand Up @@ -111,16 +110,18 @@ Vec2 DistortionBrown::removeDistortion(const Vec2& p) const
return undistorted_value;
}

Eigen::Matrix2d DistortionBrown::getDerivativeRemoveDistoWrtPt([[maybe_unused]] const Vec2& p) const
Eigen::Matrix2d DistortionBrown::getDerivativeRemoveDistoWrtPt(const Vec2& p) const
{
ALICEVISION_THROW_ERROR("Brown inverse jacobian are not implemented");
return Eigen::Matrix2d::Identity();
// Inverse of the Jacobian of addDistortion, at the undistorted point.
return getDerivativeAddDistoWrtPt(removeDistortion(p)).inverse();
}

Eigen::MatrixXd DistortionBrown::getDerivativeRemoveDistoWrtDisto([[maybe_unused]] const Vec2& p) const
Eigen::MatrixXd DistortionBrown::getDerivativeRemoveDistoWrtDisto(const Vec2& p) const
{
ALICEVISION_THROW_ERROR("Brown inverse jacobian are not implemented");
return Eigen::Matrix<double, 2, 5>::Zero();
// Exact derivative of the inverse of addDistortion: with p = add(q, theta) and p fixed,
// dq/dtheta = -(dadd/dq)^-1 * dadd/dtheta, at the undistorted point q.
const Vec2 undistorted = removeDistortion(p);
return -getDerivativeAddDistoWrtPt(undistorted).inverse() * getDerivativeAddDistoWrtDisto(undistorted);
}

} // namespace camera
Expand Down
51 changes: 4 additions & 47 deletions src/aliceVision/camera/DistortionFisheye.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,53 +151,10 @@ Eigen::Matrix2d DistortionFisheye::getDerivativeRemoveDistoWrtPt(const Vec2& p)

Eigen::MatrixXd DistortionFisheye::getDerivativeRemoveDistoWrtDisto(const Vec2& p) const
{
double r_dist = sqrt(p(0) * p(0) + p(1) * p(1));
const double eps = 1e-8;
if (r_dist < eps)
{
return Eigen::Matrix<double, 2, 4>::Zero();
}

const Vec2 p_undist = removeDistortion(p);
const double r = sqrt(p_undist(0) * p_undist(0) + p_undist(1) * p_undist(1));

const double& k1 = _distortionParams.at(0);
const double& k2 = _distortionParams.at(1);
const double& k3 = _distortionParams.at(2);
const double& k4 = _distortionParams.at(3);

const double theta = std::atan(r);
const double theta2 = theta * theta;
const double theta3 = theta2 * theta;
const double theta4 = theta2 * theta2;
const double theta5 = theta4 * theta;
const double theta6 = theta3 * theta3;
const double theta7 = theta6 * theta;
const double theta8 = theta4 * theta4;
const double theta9 = theta8 * theta;
const double theta_dist = theta + k1 * theta3 + k2 * theta5 + k3 * theta7 + k4 * theta9;
const double r_coeff = theta_dist / r;
double d_r_coeff_d_theta_dist = 1.0 / r;

Eigen::Matrix<double, 1, 4> d_r_theta_dist_d_params;
d_r_theta_dist_d_params(0, 0) = theta3;
d_r_theta_dist_d_params(0, 1) = theta5;
d_r_theta_dist_d_params(0, 2) = theta7;
d_r_theta_dist_d_params(0, 3) = theta9;

Eigen::Matrix<double, 1, 4> d_rcoeff_d_params = d_r_coeff_d_theta_dist * d_r_theta_dist_d_params;

Eigen::Matrix<double, 2, 4> ret;
ret(0, 0) = -(p(0) * d_rcoeff_d_params(0, 0)) / (r_coeff * r_coeff);
ret(0, 1) = -(p(0) * d_rcoeff_d_params(0, 1)) / (r_coeff * r_coeff);
ret(0, 2) = -(p(0) * d_rcoeff_d_params(0, 2)) / (r_coeff * r_coeff);
ret(0, 3) = -(p(0) * d_rcoeff_d_params(0, 3)) / (r_coeff * r_coeff);
ret(1, 0) = -(p(1) * d_rcoeff_d_params(0, 0)) / (r_coeff * r_coeff);
ret(1, 1) = -(p(1) * d_rcoeff_d_params(0, 1)) / (r_coeff * r_coeff);
ret(1, 2) = -(p(1) * d_rcoeff_d_params(0, 2)) / (r_coeff * r_coeff);
ret(1, 3) = -(p(1) * d_rcoeff_d_params(0, 3)) / (r_coeff * r_coeff);

return ret;
// Exact derivative of the inverse of addDistortion: with p = add(q, theta) and p fixed,
// dq/dtheta = -(dadd/dq)^-1 * dadd/dtheta, at the undistorted point q.
const Vec2 undistorted = removeDistortion(p);
return -getDerivativeAddDistoWrtPt(undistorted).inverse() * getDerivativeAddDistoWrtDisto(undistorted);
}

} // namespace camera
Expand Down
34 changes: 4 additions & 30 deletions src/aliceVision/camera/DistortionFisheye1.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,36 +127,10 @@ Eigen::Matrix2d DistortionFisheye1::getDerivativeRemoveDistoWrtPt(const Vec2& p)

Eigen::MatrixXd DistortionFisheye1::getDerivativeRemoveDistoWrtDisto(const Vec2& p) const
{
const double& k1 = _distortionParams.at(0);
double r_dist = sqrt(p(0) * p(0) + p(1) * p(1));
const double eps = 1e-8;
if (k1 * r_dist < eps)
{
return Eigen::Matrix<double, 2, 1>::Zero();
}

const Vec2 p_undist = removeDistortion(p);
const double r = sqrt(p_undist(0) * p_undist(0) + p_undist(1) * p_undist(1));

const double part1 = 2.0 * r * std::tan(0.5 * k1);
const double part2 = std::atan(part1);
const double part3 = k1 * r;
const double coef = part2 / part3;

double ca = cos(0.5 * k1);
double d_part1_d_params = r / (ca * ca);
double d_part3_d_params = r;
double d_part2_d_part1 = 1.0 / (part1 * part1 + 1.0);
double d_coef_d_part2 = 1.0 / part3;
double d_coef_d_part3 = -part2 / (part3 * part3);
double d_coef_d_params = d_coef_d_part3 * d_part3_d_params + d_coef_d_part2 * d_part2_d_part1 * d_part1_d_params;

// p'/coef
Eigen::Matrix<double, 2, 1> ret;
ret(0, 0) = -p(0) * d_coef_d_params / (coef * coef);
ret(1, 0) = -p(1) * d_coef_d_params / (coef * coef);

return ret;
// Exact derivative of the inverse of addDistortion: with p = add(q, theta) and p fixed,
// dq/dtheta = -(dadd/dq)^-1 * dadd/dtheta, at the undistorted point q.
const Vec2 undistorted = removeDistortion(p);
return -getDerivativeAddDistoWrtPt(undistorted).inverse() * getDerivativeAddDistoWrtDisto(undistorted);
}

} // namespace camera
Expand Down
114 changes: 12 additions & 102 deletions src/aliceVision/camera/DistortionRadial.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,30 +69,10 @@ Eigen::Matrix2d DistortionRadialK1::getDerivativeRemoveDistoWrtPt(const Vec2& p)

Eigen::MatrixXd DistortionRadialK1::getDerivativeRemoveDistoWrtDisto(const Vec2& p) const
{
double r_dist = sqrt(p(0) * p(0) + p(1) * p(1));
const double eps = 1e-8;
if (r_dist < eps)
{
return Eigen::Matrix<double, 2, 1>::Zero();
}

const Vec2 p_undist = removeDistortion(p);

const double& k1 = _distortionParams[0];

const double r = sqrt(p_undist(0) * p_undist(0) + p_undist(1) * p_undist(1));
const double r2 = r * r;

const double r_coeff = 1.0 + k1 * r2;

Eigen::Matrix<double, 1, 1> d_rcoeff_d_params;
d_rcoeff_d_params(0, 0) = r2;

Eigen::Matrix<double, 2, 1> ret;
ret(0, 0) = -(p(0) * d_rcoeff_d_params(0, 0)) / (r_coeff * r_coeff);
ret(1, 0) = -(p(1) * d_rcoeff_d_params(0, 0)) / (r_coeff * r_coeff);

return ret;
// Exact derivative of the inverse of addDistortion: with p = add(q, theta) and p fixed,
// dq/dtheta = -(dadd/dq)^-1 * dadd/dtheta, at the undistorted point q.
const Vec2 undistorted = removeDistortion(p);
return -getDerivativeAddDistoWrtPt(undistorted).inverse() * getDerivativeAddDistoWrtDisto(undistorted);
}

double DistortionRadialK1::getUndistortedRadius(double r) const
Expand Down Expand Up @@ -203,40 +183,10 @@ Eigen::Matrix2d DistortionRadialK3::getDerivativeRemoveDistoWrtPt(const Vec2& p)

Eigen::MatrixXd DistortionRadialK3::getDerivativeRemoveDistoWrtDisto(const Vec2& p) const
{
double r_dist = sqrt(p(0) * p(0) + p(1) * p(1));
const double eps = 1e-8;
if (r_dist < eps)
{
return Eigen::Matrix<double, 2, 3>::Zero();
}

Vec2 p_undist = removeDistortion(p);

const double& k1 = _distortionParams[0];
const double& k2 = _distortionParams[1];
const double& k3 = _distortionParams[2];

const double r = sqrt(p_undist(0) * p_undist(0) + p_undist(1) * p_undist(1));
const double r2 = r * r;
const double r4 = r2 * r2;
const double r6 = r4 * r2;

const double r_coeff = 1.0 + k1 * r2 + k2 * r4 + k3 * r6;

Eigen::Matrix<double, 1, 3> d_rcoeff_d_params;
d_rcoeff_d_params(0, 0) = r2;
d_rcoeff_d_params(0, 1) = r4;
d_rcoeff_d_params(0, 2) = r6;

Eigen::Matrix<double, 2, 3> ret;
ret(0, 0) = -(p(0) * d_rcoeff_d_params(0, 0)) / (r_coeff * r_coeff);
ret(0, 1) = -(p(0) * d_rcoeff_d_params(0, 1)) / (r_coeff * r_coeff);
ret(0, 2) = -(p(0) * d_rcoeff_d_params(0, 2)) / (r_coeff * r_coeff);
ret(1, 0) = -(p(1) * d_rcoeff_d_params(0, 0)) / (r_coeff * r_coeff);
ret(1, 1) = -(p(1) * d_rcoeff_d_params(0, 1)) / (r_coeff * r_coeff);
ret(1, 2) = -(p(1) * d_rcoeff_d_params(0, 2)) / (r_coeff * r_coeff);

return ret;
// Exact derivative of the inverse of addDistortion: with p = add(q, theta) and p fixed,
// dq/dtheta = -(dadd/dq)^-1 * dadd/dtheta, at the undistorted point q.
const Vec2 undistorted = removeDistortion(p);
return -getDerivativeAddDistoWrtPt(undistorted).inverse() * getDerivativeAddDistoWrtDisto(undistorted);
}

double DistortionRadialK3::getUndistortedRadius(double r) const
Expand Down Expand Up @@ -361,50 +311,10 @@ Eigen::Matrix2d DistortionRadialK3PT::getDerivativeRemoveDistoWrtPt(const Vec2&

Eigen::MatrixXd DistortionRadialK3PT::getDerivativeRemoveDistoWrtDisto(const Vec2& p) const
{
const double& k1 = _distortionParams[0];
const double& k2 = _distortionParams[1];
const double& k3 = _distortionParams[2];

const double r_dist = sqrt(p(0) * p(0) + p(1) * p(1));
const double eps = 1e-8;
if (r_dist < eps)
{
return Eigen::Matrix<double, 2, 3>::Zero();
}

const Vec2 p_undist = removeDistortion(p);
const double r = sqrt(p_undist(0) * p_undist(0) + p_undist(1) * p_undist(1));

const double r2 = r * r;
const double r4 = r2 * r2;
const double r6 = r4 * r2;

const double denum = (1.0 + k1 + k2 + k3);
const double num = (1.0 + k1 * r2 + k2 * r4 + k3 * r6);
const double r_coeff = num / denum;

double denum2 = denum * denum;
Eigen::Matrix<double, 1, 3> d_num_d_params;
d_num_d_params(0, 0) = r2;
d_num_d_params(0, 1) = r4;
d_num_d_params(0, 2) = r6;

Eigen::Matrix<double, 1, 3> d_denum_d_params;
d_denum_d_params(0, 0) = 1;
d_denum_d_params(0, 1) = 1;
d_denum_d_params(0, 2) = 1;

Eigen::Matrix<double, 1, 3> d_rcoeff_d_params = (denum * d_num_d_params - num * d_denum_d_params) / denum2;

Eigen::Matrix<double, 2, 3> ret;
ret(0, 0) = -(p(0) * d_rcoeff_d_params(0, 0)) / (r_coeff * r_coeff);
ret(0, 1) = -(p(0) * d_rcoeff_d_params(0, 1)) / (r_coeff * r_coeff);
ret(0, 2) = -(p(0) * d_rcoeff_d_params(0, 2)) / (r_coeff * r_coeff);
ret(1, 0) = -(p(1) * d_rcoeff_d_params(0, 0)) / (r_coeff * r_coeff);
ret(1, 1) = -(p(1) * d_rcoeff_d_params(0, 1)) / (r_coeff * r_coeff);
ret(1, 2) = -(p(1) * d_rcoeff_d_params(0, 2)) / (r_coeff * r_coeff);

return ret;
// Exact derivative of the inverse of addDistortion: with p = add(q, theta) and p fixed,
// dq/dtheta = -(dadd/dq)^-1 * dadd/dtheta, at the undistorted point q.
const Vec2 undistorted = removeDistortion(p);
return -getDerivativeAddDistoWrtPt(undistorted).inverse() * getDerivativeAddDistoWrtDisto(undistorted);
}

Vec2 DistortionRadialK3PT::removeDistortion(const Vec2& p) const
Expand Down
64 changes: 36 additions & 28 deletions src/aliceVision/camera/Equidistant.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,15 +14,13 @@ namespace camera {

std::shared_ptr<Equidistant> Equidistant::cast(std::shared_ptr<IntrinsicBase> sptr) { return std::dynamic_pointer_cast<Equidistant>(sptr); }

Vec2 Equidistant::transformProject(const Eigen::Matrix4d& pose, const Vec4& pt, bool applyDistortion) const
Vec2 Equidistant::cameraPlanePoint(const Vec4& X) const
{
const double rsensor = std::min(sensorWidth(), sensorHeight());
const double rscale = sensorWidth() / std::max(w(), h());
const double fmm = _scale(0) * rscale;
const double fov = rsensor / fmm;

const Vec4 X = pose * pt;

// Compute angle with optical center
const double angle_Z = std::atan2(sqrt(X(0) * X(0) + X(1) * X(1)), X(2));

Expand All @@ -32,36 +30,24 @@ Vec2 Equidistant::transformProject(const Eigen::Matrix4d& pose, const Vec4& pt,
const double radius = angle_Z / (0.5 * fov);

// radius = focal * angle_Z
const Vec2 P{cos(angle_radial) * radius, sin(angle_radial) * radius};

const Vec2 pt_disto = applyDistortion ? this->addDistortion(P) : P;
const Vec2 pt_ima = this->cam2ima(pt_disto);
return Vec2(cos(angle_radial) * radius, sin(angle_radial) * radius);
}

return pt_ima;
Vec2 Equidistant::transformProject(const Eigen::Matrix4d& pose, const Vec4& pt, bool applyDistortion) const
{
return project(pose * pt, applyDistortion);
}

Vec2 Equidistant::project(const Vec4& pt, bool applyDistortion) const
{
const double rsensor = std::min(sensorWidth(), sensorHeight());
const double rscale = sensorWidth() / std::max(w(), h());
const double fmm = _scale(0) * rscale;
const double fov = rsensor / fmm;

// Compute angle with optical center
const double angle_Z = std::atan2(sqrt(pt(0) * pt(0) + pt(1) * pt(1)), pt(2));

// Ignore depth component and compute radial angle
const double angle_radial = std::atan2(pt(1), pt(0));

const double radius = angle_Z / (0.5 * fov);

// radius = focal * angle_Z
const Vec2 P{cos(angle_radial) * radius, sin(angle_radial) * radius};

const Vec2 P = cameraPlanePoint(pt);
const Vec2 pt_disto = applyDistortion ? this->addDistortion(P) : P;
const Vec2 pt_ima = this->cam2ima(pt_disto);
return this->cam2ima(pt_disto);
}

return pt_ima;
Vec2 Equidistant::transformProjectDistortionOnly(const Eigen::Matrix4d& pose, const Vec4& pt) const
{
return this->cam2ima(this->addDistortionOnly(cameraPlanePoint(pose * pt)));
}

Eigen::Matrix<double, 2, 3> Equidistant::getDerivativeTransformProjectWrtPoint3(const Eigen::Matrix4d& T, const Vec4& pt) const
Expand Down Expand Up @@ -263,7 +249,7 @@ Eigen::Matrix<double, 3, 2> Equidistant::getDerivativetoUnitSphereWrtScale(const
Eigen::Matrix<double, 3, Eigen::Dynamic> Equidistant::getDerivativeBackProjectUnitWrtParams(const Vec2& pt2D) const
{
const Vec2 ptMeters = ima2cam(pt2D);
const Vec2 ptUndist = removeDistortion(ptMeters);
const Vec2 ptUndist = removeDistortionOnly(ptMeters);

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> J(3, getParameters().size());

Expand All @@ -279,7 +265,7 @@ Eigen::Matrix<double, 3, Eigen::Dynamic> Equidistant::getDerivativeBackProjectUn
size_t disto_size = getDistortionParamsSize();

const Vec2 ptMeters = ima2cam(pt2D);
const Vec2 ptUndist = removeDistortion(ptMeters);
const Vec2 ptUndist = removeDistortionOnly(ptMeters);

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> J(3, disto_size);
J = getDerivativetoUnitSphereWrtPoint(ptUndist) * getDerivativeRemoveDistoWrtDisto(ptMeters);
Expand Down Expand Up @@ -317,6 +303,28 @@ bool Equidistant::isVisibleRay(const Vec3& ray) const

EINTRINSIC Equidistant::getType() const { return EINTRINSIC::EQUIDISTANT_CAMERA; }

bool Equidistant::operator==(const IntrinsicBase& otherBase) const
{
if (!IntrinsicScaleOffsetDisto::operator==(otherBase))
{
return false;
}

// Same type, checked by the parent class.
const Equidistant& other = static_cast<const Equidistant&>(otherBase);

return _circleRadius == other._circleRadius && _circleCenter == other._circleCenter;
}

void Equidistant::rescale(float factorW, float factorH)
{
IntrinsicScaleOffsetDisto::rescale(factorW, factorH);

_circleRadius *= factorW;
_circleCenter(0) *= factorW;
_circleCenter(1) *= factorH;
}

double Equidistant::getHorizontalFov() const
{
const double rsensor = std::min(sensorWidth(), sensorHeight());
Expand Down
Loading
Loading