Use time instead of range to test geo2rdr convergence (#886)

* Control geo2rdr convergence  with time.

* Clarify cuda Geo2rdr test (still meets old metric).

* Fix up pybind cuda Geo2rdr test

* Test using max rather than RMS error.

* Update documentation of geo2rdr threshold.

* Use consistent geo2rdr thresholds in bindings.

* Apply time step after convergence check so returned range and time are consistent.

* Make CPU, GPU, and Python geo2rdr tests consistent and allow for error slightly larger than tolerance.

* Update comment to better describe new test metric.

* Ensure (time, range) are consistent even when geo2rdr fails to converge.

* Repeat (t,r) consistency fix in other two geo2rdr implementations.

Author: bhawkins

cxx/isce3/cuda/geometry/gpuGeometry.cu

@@ -74,10 +74,14 @@ int geo2rdr(const Vec3& inputLLH, const isce3::core::Ellipsoid& ellipsoid,
     // Use mid-orbit epoch as initial guess
     aztime = orbit.midTime();
 
+    // Newton step, init to zero.
+    double dt = 0.0;
+
     // Begin iterations
     int converged = 0;
-    double slantRange_old = 0.0;
     for (int i = 0; i < maxIter; ++i) {
+        // Apply Newton step from previous iteration.
+        aztime -= dt;
 
         // Interpolate the orbit to current estimate of azimuth time
         Vec3 pos, vel;
@@ -95,16 +99,6 @@ int geo2rdr(const Vec3& inputLLH, const isce3::core::Ellipsoid& ellipsoid,
             return converged;
         }
 
-        // Check convergence
-        if (std::abs(slantRange - slantRange_old) < threshold) {
-            converged = 1;
-            *slantRange_result = slantRange;
-            *aztime_result = aztime;
-            return converged;
-        } else {
-            slantRange_old = slantRange;
-        }
-
         // Compute doppler
         const double dopfact = dr.dot(vel);
         const double fdop = doppler.eval(slantRange) * dopscale;
@@ -119,8 +113,17 @@ int geo2rdr(const Vec3& inputLLH, const isce3::core::Ellipsoid& ellipsoid,
         const double c2 = (fdop / slantRange) + fdopder;
         const double fnprime = c1 + c2 * dopfact;
 
-        // Update guess for azimuth time
-        aztime -= fn / fnprime;
+        // Compute Newton step.  Don't apply until start of next iteration so
+        // that returned (slantRange, aztime) are always consistent.
+        dt = fn / fnprime;
+
+        // Check convergence
+        if (std::abs(dt) < threshold) {
+            converged = 1;
+            *slantRange_result = slantRange;
+            *aztime_result = aztime;
+            return converged;
+        }
     }
 
     // If we reach this point, no convergence for specified threshold

cxx/isce3/geocode/GeocodePolygon.h

@@ -33,7 +33,7 @@ class GeocodePolygon {
      * @param[in]  orbit               Orbit
      * @param[in]  input_dop           Doppler LUT associated with the radar grid
      * @param[in]  dem_raster          Input DEM raster
-     * @param[in]  threshold           Distance threshold for convergence
+     * @param[in]  threshold           Azimuth time threshold for convergence (s)
      * @param[in]  num_iter            Maximum number of Newton-Raphson iterations
      * @param[in]  delta_range         Step size used for computing Doppler derivative
      */

cxx/isce3/geometry/RTC.h

@@ -155,7 +155,7 @@ void computeRtc(isce3::product::Product& product, isce3::io::Raster& dem_raster,
  * looks associated with the geogrid will be saved
  * @param[in]  rtc_memory_mode     Select memory mode
  * @param[in]  interp_method       Interpolation Method
- * @param[in]  threshold           Distance threshold for convergence
+ * @param[in]  threshold           Azimuth time threshold for convergence (s)
  * @param[in]  num_iter            Maximum number of Newton-Raphson iterations
  * @param[in]  delta_range         Step size used for computing derivative of
  * doppler
@@ -214,7 +214,7 @@ void computeRtc(const isce3::product::RadarGridParameters& radarGrid,
  * looks associated with the geogrid will be saved
  * @param[in]  rtc_memory_mode     Select memory mode
  * @param[in]  interp_method       Interpolation Method
- * @param[in]  threshold           Distance threshold for convergence
+ * @param[in]  threshold           Azimuth time threshold for convergence (s)
  * @param[in]  num_iter            Maximum number of Newton-Raphson iterations
  * @param[in]  delta_range         Step size used for computing derivative of
  * doppler
@@ -301,7 +301,7 @@ void computeRtcBilinearDistribution(isce3::io::Raster& dem_raster,
  * looks associated with the geogrid will be saved
  * @param[in] rtc_memory_mode      Select memory mode
  * @param[in] interp_method        Interpolation Method
- * @param[in] threshold            Distance threshold for convergence
+ * @param[in] threshold            Azimuth time threshold for convergence (s)
  * @param[in] num_iter             Maximum number of Newton-Raphson iterations
  * @param[in] delta_range          Step size used for computing derivative of
  * doppler

cxx/isce3/geometry/detail/Geo2Rdr.h

@@ -10,7 +10,7 @@ namespace isce3 { namespace geometry { namespace detail {
 
 /** \internal Root-finding configuration parameters for geo2rdr */
 struct Geo2RdrParams {
-    /** \internal Absolute slant range convergence tolerance (m) */
+    /** \internal Absolute azimuth time convergence tolerance (s) */
     double threshold = 1e-8;
 
     /** \internal Maximum number of Newton-Raphson iterations */

cxx/isce3/geometry/detail/Geo2Rdr.icc

@@ -136,9 +136,13 @@ geo2rdr(double* t, double* r, const isce3::core::Vec3& llh,
         }
     }
 
+    // Newton step, initialized to zero.
+    double dt = 0.0;
+
     // begin iterations
-    double r_old = 0.;
     for (int i = 0; i < params.maxiter; ++i) {
+        // apply Newton step here so that (r,t) are always consistent on return.
+        *t -= dt;
 
         // interpolate orbit
         Vec3 pos, vel;
@@ -156,21 +160,17 @@ geo2rdr(double* t, double* r, const isce3::core::Vec3& llh,
             }
         }
 
-        // check for convergence
-        if (std::abs(*r - r_old) < params.threshold) {
-            return ErrorCode::Success;
-        }
-
         // update guess for azimuth time
-        double dt;
         const auto status = computeDopplerAztimeDiff(&dt, *t, *r, rvec, vel,
                                                      doppler, wvl, params.delta_range);
         if (status != ErrorCode::Success) {
             return status;
         }
 
-        *t -= dt;
-        r_old = *r;
+        // check for convergence
+        if (std::abs(dt) < params.threshold) {
+            return ErrorCode::Success;
+        }
     }
 
     // if we reach this point, no convergence for specified threshold

cxx/isce3/geometry/geometry.cpp

@@ -198,9 +198,14 @@ int isce3::geometry::geo2rdr(const Vec3& inputLLH, const Ellipsoid& ellipsoid,
     if (error)
         return !converged;
 
-    double slantRange_old = slantRange;
+    // Newton step, initialized to zero.
+    double aztime_diff = 0.0;
+
     // Begin iterations
     for (int i = 0; i < maxIter; ++i) {
+        // Apply Newton step computed in previous iteration here so that
+        // (aztime, slantRange) are always consistent on return.
+        aztime -= aztime_diff;
 
         // Interpolate the orbit to current estimate of azimuth time
         orbit.interpolate(
@@ -219,17 +224,15 @@ int isce3::geometry::geo2rdr(const Vec3& inputLLH, const Ellipsoid& ellipsoid,
         }
 
         slantRange = dr.norm();
-        // Check convergence
-        if (std::abs(slantRange - slantRange_old) < threshold)
-            return converged;
-        else
-            slantRange_old = slantRange;
 
         // Update guess for azimuth time
-        double aztime_diff = _compute_doppler_aztime_diff(dr, satvel, doppler,
+        aztime_diff = _compute_doppler_aztime_diff(dr, satvel, doppler,
                 wavelength, aztime, slantRange, deltaRange);
 
-        aztime -= aztime_diff;
+        // Check convergence
+        if (std::abs(aztime_diff) < threshold) {
+            return converged;
+        }
     }
     // If we reach this point, no convergence for specified threshold
     return !converged;

cxx/isce3/geometry/metadataCubes.h

@@ -112,7 +112,7 @@ inline void writeVectorDerivedCubes(const int array_pos_i,
  * @param[out] elevation_angle_raster      Elevation angle (in degrees wrt
  * geodedic nadir) cube raster
  * @param[out] ground_track_velocity_raster Ground-track velocity raster
- * @param[in]  threshold_geo2rdr           Range threshold for geo2rdr
+ * @param[in]  threshold_geo2rdr           Azimuth time threshold for geo2rdr
  * @param[in]  numiter_geo2rdr             Geo2rdr maximum number of iterations
  * @param[in]  delta_range                 Step size used for computing
  * derivative of doppler
@@ -183,7 +183,7 @@ void makeRadarGridCubes(const isce3::product::RadarGridParameters& radar_grid,
  * @param[out] elevation_angle_raster    Elevation angle (in degrees wrt
  * geodedic nadir) cube raster
  * @param[out] ground_track_velocity_raster Ground-track velocity raster
- * @param[in]  threshold_geo2rdr         Range threshold for geo2rdr
+ * @param[in]  threshold_geo2rdr         Azimuth time threshold for geo2rdr
  * @param[in]  numiter_geo2rdr           Geo2rdr maximum number of iterations
  * @param[in]  delta_range               Step size used for computing
  * derivative of doppler

python/extensions/pybind_isce3/cuda/geocode/cuGeocode.cpp

@@ -5,6 +5,7 @@
 
 #include <isce3/container/RadarGeometry.h>
 #include <isce3/core/Constants.h>
+#include <isce3/geometry/detail/Geo2Rdr.h>
 #include <isce3/io/Raster.h>
 #include <isce3/product/GeoGridParameters.h>
 
@@ -14,6 +15,7 @@ using isce3::cuda::geocode::Geocode;
 
 void addbinding(pybind11::class_<Geocode>& pyGeocode)
 {
+    const isce3::geometry::detail::Geo2RdrParams defaults;
     pyGeocode
             .def(py::init<const isce3::product::GeoGridParameters&,
                          const isce3::container::RadarGeometry&,
@@ -28,8 +30,10 @@ void addbinding(pybind11::class_<Geocode>& pyGeocode)
                             isce3::core::BILINEAR_METHOD,
                     py::arg("dem_interp_method") =
                             isce3::core::BIQUINTIC_METHOD,
-                    py::arg("threshold") = 1e-8, py::arg("maxiter") = 50,
-                    py::arg("delta_range") = 10, py::arg("invalid_value") = 0.0,
+                    py::arg("threshold") = defaults.threshold,
+                    py::arg("maxiter") = defaults.maxiter,
+                    py::arg("delta_range") = defaults.delta_range,
+                    py::arg("invalid_value") = 0.0,
                     R"(
             Create CUDA geocode object.
 

python/extensions/pybind_isce3/cuda/geometry/geo2rdr.cpp

@@ -6,6 +6,7 @@
 #include <isce3/core/Ellipsoid.h>
 #include <isce3/core/LUT2d.h>
 #include <isce3/core/Orbit.h>
+#include <isce3/geometry/detail/Geo2Rdr.h>
 #include <isce3/geometry/geometry.h>
 #include <isce3/io/Raster.h>
 #include <isce3/product/RadarGridParameters.h>
@@ -16,6 +17,7 @@ namespace py = pybind11;
 
 void addbinding(py::class_<Geo2rdr> & pyGeo2Rdr)
 {
+    const isce3::geometry::detail::Geo2RdrParams defaults;
     pyGeo2Rdr
         .def(py::init([](const isce3::product::RadarGridParameters & radar_grid,
              const isce3::core::Orbit & orbit,
@@ -35,8 +37,8 @@ void addbinding(py::class_<Geo2rdr> & pyGeo2Rdr)
             py::arg("orbit"),
             py::arg("ellipsoid"),
             py::arg("doppler") = isce3::core::LUT2d<double>(),
-            py::arg("threshold") = 1e-8,
-            py::arg("numiter") = 25,
+            py::arg("threshold") = defaults.threshold,
+            py::arg("numiter") = defaults.maxiter,
             py::arg("lines_per_block") = 1000)
         .def("geo2rdr", py::overload_cast<isce3::io::Raster &, const std::string &,
                 double, double>

python/extensions/pybind_isce3/geocode/GeocodePolygon.cpp

@@ -47,7 +47,7 @@ void addbinding(py::class_<GeocodePolygon<T>> &pyGeocodePolygon)
         orbit               Orbit
         input_dop           Doppler LUT associated with the radar grid
         dem_raster          Input DEM raster
-        threshold           Distance threshold for convergence
+        threshold           Azimuth time threshold for convergence (s)
         num_iter            Maximum number of Newton-Raphson iterations
         delta_range         Step size used for computing Doppler derivative
         )")