From 88eb5896fa54fa3a8e5bfcdab56a06ed36a2a751 Mon Sep 17 00:00:00 2001
From: Austin Sanders <austinsanders1993@gmail.com>
Date: Tue, 25 Feb 2025 13:00:49 -0700
Subject: [PATCH] Cahvor fix (#490)
* Removed tolerance check from optical shift addition
* Updated changelog
---
CHANGELOG.md | 5 ++-
src/Distortion.cpp | 96 ++++++++++++++++++++++------------------------
2 files changed, 49 insertions(+), 52 deletions(-)
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2c170f8..82891ec 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -35,6 +35,9 @@ release.
## [Unreleased]
+### Fixed
+- Fixed CAHVOR model optical shifts by removing tolerance check [#488](https://github.com/DOI-USGS/usgscsm/issues/488)
+
## [2.0.1] - 2024-01-23
### Changed
@@ -51,4 +54,4 @@ release.
- Updated ALE submodule [#470](https://github.com/DOI-USGS/usgscsm/pull/470)
### Fixed
-- Fixed issue with radial distortion computation [#464](https://github.com/DOI-USGS/usgscsm/pull/464)
\ No newline at end of file
+- Fixed issue with radial distortion computation [#464](https://github.com/DOI-USGS/usgscsm/pull/464)
diff --git a/src/Distortion.cpp b/src/Distortion.cpp
index e69a5cd..2f9a382 100644
--- a/src/Distortion.cpp
+++ b/src/Distortion.cpp
@@ -328,16 +328,13 @@ void removeDistortion(double dx, double dy, double &ux, double &uy,
double shiftedDy = dy - opticalDistCoeffs[4];
double rr = shiftedDx * shiftedDx + shiftedDy * shiftedDy;
- if (rr > tolerance)
- {
- double dr = opticalDistCoeffs[0] +
- (rr * (opticalDistCoeffs[1] + rr * opticalDistCoeffs[2]));
+ double dr = opticalDistCoeffs[0] +
+ (rr * (opticalDistCoeffs[1] + rr * opticalDistCoeffs[2]));
- ux = shiftedDx * (1.0 - dr);
- uy = shiftedDy * (1.0 - dr);
- ux += opticalDistCoeffs[3];
- uy += opticalDistCoeffs[4];
- }
+ ux = shiftedDx * (1.0 - dr);
+ uy = shiftedDy * (1.0 - dr);
+ ux += opticalDistCoeffs[3];
+ uy += opticalDistCoeffs[4];
}
break;
@@ -587,51 +584,48 @@ void applyDistortion(double ux, double uy, double &dx, double &dy,
double shiftedUx = ux - opticalDistCoeffs[3];
double shiftedUy = uy - opticalDistCoeffs[4];
double rp2 = (ux * ux) + (uy * uy);
+ double rp = sqrt(rp2);
+ // Compute first fractional distortion using rp
+ double drOverR =
+ opticalDistCoeffs[0] +
+ (rp2 * (opticalDistCoeffs[1] + (rp2 * opticalDistCoeffs[2])));
- if (rp2 > tolerance)
+ // Compute first distorted point estimate, r
+ double r = rp + (drOverR * rp);
+ double r_prev, r2_prev;
+ int iteration = 0;
+
+ do
{
- double rp = sqrt(rp2);
- // Compute first fractional distortion using rp
- double drOverR =
- opticalDistCoeffs[0] +
- (rp2 * (opticalDistCoeffs[1] + (rp2 * opticalDistCoeffs[2])));
-
- // Compute first distorted point estimate, r
- double r = rp + (drOverR * rp);
- double r_prev, r2_prev;
- int iteration = 0;
-
- do
+ // Don't get in an end-less loop. This algorithm should
+ // converge quickly. If not then we are probably way outside
+ // of the focal plane. Just set the distorted position to the
+ // undistorted position. Also, make sure the focal plane is less
+ // than 1km, it is unreasonable for it to grow larger than that.
+ if (iteration >= 20 || r > 1E9)
{
- // Don't get in an end-less loop. This algorithm should
- // converge quickly. If not then we are probably way outside
- // of the focal plane. Just set the distorted position to the
- // undistorted position. Also, make sure the focal plane is less
- // than 1km, it is unreasonable for it to grow larger than that.
- if (iteration >= 20 || r > 1E9)
- {
- drOverR = 0.0;
- break;
- }
-
- r_prev = r;
- r2_prev = r * r;
-
- // Compute new fractional distortion:
- drOverR = opticalDistCoeffs[0] +
- (r2_prev *
- (opticalDistCoeffs[1] + (r2_prev * opticalDistCoeffs[2])));
-
- // Compute new estimate of r
- r = rp + (drOverR * r_prev);
- iteration++;
- } while (fabs(r - r_prev) > desiredPrecision);
-
- dx = shiftedUx / (1.0 - drOverR);
- dy = shiftedUy / (1.0 - drOverR);
- dx += opticalDistCoeffs[3];
- dy += opticalDistCoeffs[4];
- }
+ drOverR = 0.0;
+ break;
+ }
+
+ r_prev = r;
+ r2_prev = r * r;
+
+ // Compute new fractional distortion:
+ drOverR = opticalDistCoeffs[0] +
+ (r2_prev *
+ (opticalDistCoeffs[1] + (r2_prev * opticalDistCoeffs[2])));
+
+ // Compute new estimate of r
+ r = rp + (drOverR * r_prev);
+ iteration++;
+ } while (fabs(r - r_prev) > desiredPrecision);
+
+ dx = shiftedUx / (1.0 - drOverR);
+ dy = shiftedUy / (1.0 - drOverR);
+ dx += opticalDistCoeffs[3];
+ dy += opticalDistCoeffs[4];
+
}
break;
--
GitLab