1212namespace aliceVision {
1313namespace camera {
1414
15- Vec2 Equidistant::project (const Eigen::Matrix4d& pose, const Vec4& pt, bool applyDistortion) const
15+ Vec2 Equidistant::transformProject (const Eigen::Matrix4d& pose, const Vec4& pt, bool applyDistortion) const
1616{
1717 const double rsensor = std::min (sensorWidth (), sensorHeight ());
1818 const double rscale = sensorWidth () / std::max (w (), h ());
@@ -38,7 +38,31 @@ Vec2 Equidistant::project(const Eigen::Matrix4d& pose, const Vec4& pt, bool appl
3838 return pt_ima;
3939}
4040
41- Eigen::Matrix<double , 2 , 9 > Equidistant::getDerivativeProjectWrtRotation (const Eigen::Matrix4d& pose, const Vec4& pt) const
41+ Vec2 Equidistant::project (const Vec4& pt, bool applyDistortion) const
42+ {
43+ const double rsensor = std::min (sensorWidth (), sensorHeight ());
44+ const double rscale = sensorWidth () / std::max (w (), h ());
45+ const double fmm = _scale (0 ) * rscale;
46+ const double fov = rsensor / fmm;
47+
48+ // Compute angle with optical center
49+ const double angle_Z = std::atan2 (sqrt (pt (0 ) * pt (0 ) + pt (1 ) * pt (1 )), pt (2 ));
50+
51+ // Ignore depth component and compute radial angle
52+ const double angle_radial = std::atan2 (pt (1 ), pt (0 ));
53+
54+ const double radius = angle_Z / (0.5 * fov);
55+
56+ // radius = focal * angle_Z
57+ const Vec2 P{cos (angle_radial) * radius, sin (angle_radial) * radius};
58+
59+ const Vec2 pt_disto = applyDistortion ? this ->addDistortion (P) : P;
60+ const Vec2 pt_ima = this ->cam2ima (pt_disto);
61+
62+ return pt_ima;
63+ }
64+
65+ Eigen::Matrix<double , 2 , 9 > Equidistant::getDerivativeTransformProjectWrtRotation (const Eigen::Matrix4d& pose, const Vec4& pt) const
4266{
4367 Eigen::Matrix4d T = pose;
4468 const Vec4 X = T * pt; // apply pose
@@ -87,7 +111,7 @@ Eigen::Matrix<double, 2, 9> Equidistant::getDerivativeProjectWrtRotation(const E
87111 return getDerivativeCam2ImaWrtPoint () * getDerivativeAddDistoWrtPt (P) * d_P_d_angles * d_angles_d_X * d_X_d_R;
88112}
89113
90- Eigen::Matrix<double , 2 , 16 > Equidistant::getDerivativeProjectWrtPose (const Eigen::Matrix4d& pose, const Vec4& pt) const
114+ Eigen::Matrix<double , 2 , 16 > Equidistant::getDerivativeTransformProjectWrtPose (const Eigen::Matrix4d& pose, const Vec4& pt) const
91115{
92116 Eigen::Matrix4d T = pose;
93117 const Vec4 X = T * pt; // apply pose
@@ -136,7 +160,7 @@ Eigen::Matrix<double, 2, 16> Equidistant::getDerivativeProjectWrtPose(const Eige
136160 return getDerivativeCam2ImaWrtPoint () * getDerivativeAddDistoWrtPt (P) * d_P_d_angles * d_angles_d_X * d_X_d_T.block <3 , 16 >(0 , 0 );
137161}
138162
139- Eigen::Matrix<double , 2 , 16 > Equidistant::getDerivativeProjectWrtPoseLeft (const Eigen::Matrix4d& pose, const Vec4& pt) const
163+ Eigen::Matrix<double , 2 , 16 > Equidistant::getDerivativeTransformProjectWrtPoseLeft (const Eigen::Matrix4d& pose, const Vec4& pt) const
140164{
141165 Eigen::Matrix4d T = pose;
142166 const Vec4 X = T * pt; // apply pose
@@ -185,7 +209,7 @@ Eigen::Matrix<double, 2, 16> Equidistant::getDerivativeProjectWrtPoseLeft(const
185209 return getDerivativeCam2ImaWrtPoint () * getDerivativeAddDistoWrtPt (P) * d_P_d_angles * d_angles_d_X * d_X_d_T.block <3 , 16 >(0 , 0 );
186210}
187211
188- Eigen::Matrix<double , 2 , 4 > Equidistant::getDerivativeProjectWrtPoint (const Eigen::Matrix4d& pose, const Vec4& pt) const
212+ Eigen::Matrix<double , 2 , 4 > Equidistant::getDerivativeTransformProjectWrtPoint (const Eigen::Matrix4d& pose, const Vec4& pt) const
189213{
190214 Eigen::Matrix4d T = pose;
191215 const Vec4 X = T * pt; // apply pose
@@ -235,7 +259,7 @@ Eigen::Matrix<double, 2, 4> Equidistant::getDerivativeProjectWrtPoint(const Eige
235259 return getDerivativeCam2ImaWrtPoint () * getDerivativeAddDistoWrtPt (P) * d_P_d_angles * d_angles_d_X * d_X_d_pt;
236260}
237261
238- Eigen::Matrix<double , 2 , 3 > Equidistant::getDerivativeProjectWrtPoint3 (const Eigen::Matrix4d& T, const Vec4& pt) const
262+ Eigen::Matrix<double , 2 , 3 > Equidistant::getDerivativeTransformProjectWrtPoint3 (const Eigen::Matrix4d& T, const Vec4& pt) const
239263{
240264 const Vec4 X = T * pt; // apply pose
241265
@@ -284,7 +308,7 @@ Eigen::Matrix<double, 2, 3> Equidistant::getDerivativeProjectWrtPoint3(const Eig
284308 return getDerivativeCam2ImaWrtPoint () * getDerivativeAddDistoWrtPt (P) * d_P_d_angles * d_angles_d_X * d_X_d_pt;
285309}
286310
287- Eigen::Matrix<double , 2 , 3 > Equidistant::getDerivativeProjectWrtDisto (const Eigen::Matrix4d& pose, const Vec4& pt) const
311+ Eigen::Matrix<double , 2 , 3 > Equidistant::getDerivativeTransformProjectWrtDisto (const Eigen::Matrix4d& pose, const Vec4& pt) const
288312{
289313 Eigen::Matrix4d T = pose;
290314 const Vec4 X = T * pt; // apply pose
@@ -308,7 +332,7 @@ Eigen::Matrix<double, 2, 3> Equidistant::getDerivativeProjectWrtDisto(const Eige
308332 return getDerivativeCam2ImaWrtPoint () * getDerivativeAddDistoWrtDisto (P);
309333}
310334
311- Eigen::Matrix<double , 2 , 2 > Equidistant::getDerivativeProjectWrtScale (const Eigen::Matrix4d& pose, const Vec4& pt) const
335+ Eigen::Matrix<double , 2 , 2 > Equidistant::getDerivativeTransformProjectWrtScale (const Eigen::Matrix4d& pose, const Vec4& pt) const
312336{
313337 Eigen::Matrix4d T = pose;
314338 const Vec4 X = T * pt; // apply pose
@@ -343,22 +367,22 @@ Eigen::Matrix<double, 2, 2> Equidistant::getDerivativeProjectWrtScale(const Eige
343367 return getDerivativeCam2ImaWrtPoint () * getDerivativeAddDistoWrtPt (P) * d_P_d_radius * d_radius_d_fov * d_fov_d_scale;
344368}
345369
346- Eigen::Matrix<double , 2 , 2 > Equidistant::getDerivativeProjectWrtPrincipalPoint (const Eigen::Matrix4d& pose, const Vec4& pt) const
370+ Eigen::Matrix<double , 2 , 2 > Equidistant::getDerivativeTransformProjectWrtPrincipalPoint (const Eigen::Matrix4d& pose, const Vec4& pt) const
347371{
348372 return getDerivativeCam2ImaWrtPrincipalPoint ();
349373}
350374
351- Eigen::Matrix<double , 2 , Eigen::Dynamic> Equidistant::getDerivativeProjectWrtParams (const Eigen::Matrix4d& pose, const Vec4& pt3D) const
375+ Eigen::Matrix<double , 2 , Eigen::Dynamic> Equidistant::getDerivativeTransformProjectWrtParams (const Eigen::Matrix4d& pose, const Vec4& pt3D) const
352376{
353377 Eigen::Matrix<double , 2 , Eigen::Dynamic> ret (2 , getParams ().size ());
354378
355- ret.block <2 , 2 >(0 , 0 ) = getDerivativeProjectWrtScale (pose, pt3D);
356- ret.block <2 , 2 >(0 , 2 ) = getDerivativeProjectWrtPrincipalPoint (pose, pt3D);
379+ ret.block <2 , 2 >(0 , 0 ) = getDerivativeTransformProjectWrtScale (pose, pt3D);
380+ ret.block <2 , 2 >(0 , 2 ) = getDerivativeTransformProjectWrtPrincipalPoint (pose, pt3D);
357381
358382 if (_pDistortion != nullptr )
359383 {
360384 const size_t distortionSize = _pDistortion->getDistortionParametersCount ();
361- ret.block (0 , 4 , 2 , distortionSize) = getDerivativeProjectWrtDisto (pose, pt3D);
385+ ret.block (0 , 4 , 2 , distortionSize) = getDerivativeTransformProjectWrtDisto (pose, pt3D);
362386 }
363387
364388 return ret;
@@ -436,6 +460,27 @@ Eigen::Matrix<double, 3, 2> Equidistant::getDerivativetoUnitSphereWrtScale(const
436460 return d_ret_d_angles * d_angles_d_fov * d_fov_d_scale;
437461}
438462
463+ Eigen::Matrix<double , 3 , Eigen::Dynamic> Equidistant::getDerivativeBackProjectUnitWrtParams (const Vec2& pt2D) const
464+ {
465+ size_t disto_size = getDistortionParamsSize ();
466+
467+ const Vec2 ptMeters = ima2cam (pt2D);
468+ const Vec2 ptUndist = removeDistortion (ptMeters);
469+ const Vec3 ptSphere = toUnitSphere (ptUndist);
470+
471+ Eigen::Matrix<double , Eigen::Dynamic, Eigen::Dynamic> J (3 , getParams ().size ());
472+
473+ J.block <3 , 2 >(0 , 0 ) = getDerivativetoUnitSphereWrtScale (ptUndist);
474+ J.block <3 , 2 >(0 , 2 ) = getDerivativetoUnitSphereWrtPoint (ptUndist) * getDerivativeRemoveDistoWrtPt (ptMeters) * getDerivativeIma2CamWrtPrincipalPoint ();
475+
476+ if (disto_size > 0 )
477+ {
478+ J.block (0 , 4 , 3 , disto_size) = getDerivativetoUnitSphereWrtPoint (ptUndist) * getDerivativeRemoveDistoWrtDisto (ptMeters);
479+ }
480+
481+ return J;
482+ }
483+
439484double Equidistant::imagePlaneToCameraPlaneError (double value) const { return value / _scale (0 ); }
440485
441486Vec2 Equidistant::cam2ima (const Vec2& p) const { return _circleRadius * p + getPrincipalPoint (); }
@@ -459,7 +504,7 @@ bool Equidistant::isVisibleRay(const Vec3& ray) const
459504 if (std::abs (angle) > 1.2 * (0.5 * fov))
460505 return false ;
461506
462- const Vec2 proj = project (Eigen::Matrix4d::Identity (), ray.homogeneous (), true );
507+ const Vec2 proj = transformProject (Eigen::Matrix4d::Identity (), ray.homogeneous (), true );
463508 const Vec2 centered = proj - Vec2 (_circleCenter (0 ), _circleCenter (1 ));
464509 return centered.norm () <= _circleRadius;
465510}
0 commit comments