@@ -328,16 +328,13 @@ void removeDistortion(double dx, double dy, double &ux, double &uy,
328328 double shiftedDy = dy - opticalDistCoeffs[4 ];
329329 double rr = shiftedDx * shiftedDx + shiftedDy * shiftedDy;
330330
331- if (rr > tolerance)
332- {
333- double dr = opticalDistCoeffs[0 ] +
334- (rr * (opticalDistCoeffs[1 ] + rr * opticalDistCoeffs[2 ]));
331+ double dr = opticalDistCoeffs[0 ] +
332+ (rr * (opticalDistCoeffs[1 ] + rr * opticalDistCoeffs[2 ]));
335333
336- ux = shiftedDx * (1.0 - dr);
337- uy = shiftedDy * (1.0 - dr);
338- ux += opticalDistCoeffs[3 ];
339- uy += opticalDistCoeffs[4 ];
340- }
334+ ux = shiftedDx * (1.0 - dr);
335+ uy = shiftedDy * (1.0 - dr);
336+ ux += opticalDistCoeffs[3 ];
337+ uy += opticalDistCoeffs[4 ];
341338 }
342339 break ;
343340
@@ -587,51 +584,48 @@ void applyDistortion(double ux, double uy, double &dx, double &dy,
587584 double shiftedUx = ux - opticalDistCoeffs[3 ];
588585 double shiftedUy = uy - opticalDistCoeffs[4 ];
589586 double rp2 = (ux * ux) + (uy * uy);
587+ double rp = sqrt (rp2);
588+ // Compute first fractional distortion using rp
589+ double drOverR =
590+ opticalDistCoeffs[0 ] +
591+ (rp2 * (opticalDistCoeffs[1 ] + (rp2 * opticalDistCoeffs[2 ])));
590592
591- if (rp2 > tolerance)
593+ // Compute first distorted point estimate, r
594+ double r = rp + (drOverR * rp);
595+ double r_prev, r2_prev;
596+ int iteration = 0 ;
597+
598+ do
592599 {
593- double rp = sqrt (rp2);
594- // Compute first fractional distortion using rp
595- double drOverR =
596- opticalDistCoeffs[0 ] +
597- (rp2 * (opticalDistCoeffs[1 ] + (rp2 * opticalDistCoeffs[2 ])));
598-
599- // Compute first distorted point estimate, r
600- double r = rp + (drOverR * rp);
601- double r_prev, r2_prev;
602- int iteration = 0 ;
603-
604- do
600+ // Don't get in an end-less loop. This algorithm should
601+ // converge quickly. If not then we are probably way outside
602+ // of the focal plane. Just set the distorted position to the
603+ // undistorted position. Also, make sure the focal plane is less
604+ // than 1km, it is unreasonable for it to grow larger than that.
605+ if (iteration >= 20 || r > 1E9 )
605606 {
606- // Don't get in an end-less loop. This algorithm should
607- // converge quickly. If not then we are probably way outside
608- // of the focal plane. Just set the distorted position to the
609- // undistorted position. Also, make sure the focal plane is less
610- // than 1km, it is unreasonable for it to grow larger than that.
611- if (iteration >= 20 || r > 1E9 )
612- {
613- drOverR = 0.0 ;
614- break ;
615- }
616-
617- r_prev = r;
618- r2_prev = r * r;
619-
620- // Compute new fractional distortion:
621- drOverR = opticalDistCoeffs[0 ] +
622- (r2_prev *
623- (opticalDistCoeffs[1 ] + (r2_prev * opticalDistCoeffs[2 ])));
624-
625- // Compute new estimate of r
626- r = rp + (drOverR * r_prev);
627- iteration++;
628- } while (fabs (r - r_prev) > desiredPrecision);
629-
630- dx = shiftedUx / (1.0 - drOverR);
631- dy = shiftedUy / (1.0 - drOverR);
632- dx += opticalDistCoeffs[3 ];
633- dy += opticalDistCoeffs[4 ];
634- }
607+ drOverR = 0.0 ;
608+ break ;
609+ }
610+
611+ r_prev = r;
612+ r2_prev = r * r;
613+
614+ // Compute new fractional distortion:
615+ drOverR = opticalDistCoeffs[0 ] +
616+ (r2_prev *
617+ (opticalDistCoeffs[1 ] + (r2_prev * opticalDistCoeffs[2 ])));
618+
619+ // Compute new estimate of r
620+ r = rp + (drOverR * r_prev);
621+ iteration++;
622+ } while (fabs (r - r_prev) > desiredPrecision);
623+
624+ dx = shiftedUx / (1.0 - drOverR);
625+ dy = shiftedUy / (1.0 - drOverR);
626+ dx += opticalDistCoeffs[3 ];
627+ dy += opticalDistCoeffs[4 ];
628+
635629 }
636630 break ;
637631
0 commit comments