Fix SNAPHU IntegratePhase() to prevent accumulating large numerical error (#967)

* Update IntegratePhase() to resolve numerical precision issues

SNAPHU's IntegratePhase() function determines the unwrapped phase by
accumulating the discrete unwrapped phase gradients along an integration
path that spans the interferogram (tile). It first scans across the
first row, and then down each column.

In practice, this can lead to unacceptably large numerical precision
errors by summing thousands of floating point numbers.

This update modifies the algorithm to instead integrate an integer
number of cycles of 2pi in order to avoid the accumulation of floating
point rounding errors. The residual phase error is significantly
improved and, as a result, the absolute error tolerances used in the
SNAPHU unit tests could be decreased by an order of magnitude.

In addition, the order in which we integrate over the interferogram was
changed from column-by-column to row-by-row, in order to improve memory
locality for row-major data structures like those currently used in
SNAPHU.

* Update SNAPHU tests to use both abs. & rel. error tolerances

Author: gmgunter

cxx/isce3/unwrap/snaphu/snaphu_util.cpp

@@ -27,6 +27,8 @@ static
 int IsFalse(char *str);
 static
 double ModDiff(double f1, double f2);
+static
+int DiffNCycle(double f1, double f2);
 
 
 /* function: IsTrue()
@@ -105,6 +107,30 @@ double ModDiff(double f1, double f2){
 }
 
 
+/* function: DiffNCycle()
+ * -------------------
+ * Computes the number of cycles of 2pi that must be added to the
+ * difference between f1 and f2 in order to "wrap" it to the interval
+ * (-pi,pi].  f1 and f2 should be between [0,2pi).  Assumes that PI has
+ * been defined.
+ *
+ * See also: ModDiff()
+ */
+static
+int DiffNCycle(double f1, double f2){
+
+  const auto f3=f1-f2;
+  if(f3>PI){
+    return -1;
+  }
+  if(f3<=-PI){
+    return 1;
+  }
+
+  return(0);
+}
+
+
 /* function: WrapPhase()
  * ---------------------
  * Makes sure the passed float array is properly wrapped into the [0,2pi)
@@ -302,30 +328,27 @@ int CalcFlow(Array2D<float>& phase, Array2D<short>* flowsptr, long nrow, long nc
  * wrapped phase to create an unwrapped phase field.  The unwrapped
  * phase field will be the same size as the wrapped field.  The array
  * rowflow should have size N-1xM and colflow size NxM-1 where the
- * phase fields are NxM.  Output is saved to a file.
+ * phase fields are NxM.
  */
 int IntegratePhase(Array2D<float>& psi, Array2D<float>& phi, Array2D<short>& flows,
                    long nrow, long ncol){
 
-  long row, col;
-
-  auto rowflow = flows.block(0,0,nrow-1,ncol);
-  auto colflow = flows.block(nrow-1,0,nrow,ncol-1);
-
-  /* set first element as seed */
-  phi(0,0)=psi(0,0);
+  auto rowflow=flows.block(0,0,nrow-1,ncol);
+  auto colflow=flows.block(nrow-1,0,nrow,ncol-1);
 
-  /* integrate over first row */
-  for(col=1;col<ncol;col++){
-    phi(0,col)=phi(0,col-1)+(ModDiff(psi(0,col),psi(0,col-1))
-      +colflow(0,col-1)*TWOPI);
-  }
-
-  /* integrate over columns */
-  for(row=1;row<nrow;row++){
-    for(col=0;col<ncol;col++){
-      phi(row,col)=phi(row-1,col)+(ModDiff(psi(row,col),psi(row-1,col))
-        -rowflow(row-1,col)*TWOPI);
+  /* compute unwrapped phase by accumulating the total number of cycles to be
+     added along an integration path that spans each row */
+  long initcycles=0;
+  for(long row=0;row<nrow;row++){
+    if(row>0){
+      initcycles+=DiffNCycle(psi(row,0),psi(row-1,0))-rowflow(row-1,0);
+    }
+    long cycles=initcycles;
+    for(long col=0;col<ncol;col++){
+      if(col>0){
+        cycles+=DiffNCycle(psi(row,col),psi(row,col-1))+colflow(row,col-1);
+      }
+      phi(row,col)=psi(row,col)+TWOPI*cycles;
     }
   }
 

tests/python/packages/isce3/unwrap/snaphu.py

@@ -261,7 +261,7 @@ def test_smooth_cost(self, init_method):
             mphase = phase[mask]
             munw = unw_raster.data[mask]
             offset = mphase[0] - munw[0]
-            assert np.allclose(munw + offset, mphase, rtol=0.0, atol=1e-5)
+            assert np.allclose(mphase - offset, munw, rtol=1e-6, atol=1e-6)
 
         # Check the set of unique labels (masked-out pixels are labeled 0).
         unique_cc_labels = set(np.unique(ccl_raster.data))
@@ -371,7 +371,7 @@ def frac_nonzero(a):
         mphase = phase[mask]
         munw = unw_raster.data[mask]
         offset = mphase[0] - munw[0]
-        good_pixels = np.isclose(munw + offset, mphase, rtol=0.0, atol=1e-4)
+        good_pixels = np.isclose(mphase - offset, munw, rtol=1e-6, atol=1e-6)
         assert frac_nonzero(~good_pixels) < 1e-3
 
     def test_tile_mode(self):
@@ -437,4 +437,4 @@ def test_tile_mode(self):
             mphase = phase[mask]
             munw = unw_raster.data[mask]
             offset = mphase[0] - munw[0]
-            assert np.allclose(munw + offset, mphase, rtol=0.0, atol=1e-5)
+            assert np.allclose(mphase - offset, munw, rtol=1e-6, atol=1e-6)