Grow conncomp

src/dolphin/_cli_unwrap.py

@@ -77,6 +77,16 @@ def get_parser(subparser=None, subcommand_name="unwrap") -> argparse.ArgumentPar
         default=UnwrapMethod.SNAPHU.value,
         help="Choice of unwrapping algorithm to use.",
     )
+    algorithm_opts.add_argument(
+        "--run-goldstein",
+        action="store_true",
+        help="Run Goldstein filter before unwrapping.",
+    )
+    algorithm_opts.add_argument(
+        "--run-interpolation",
+        action="store_true",
+        help="Run interpolation before unwrapping.",
+    )
 
     tophu_opts = parser.add_argument_group("Tophu options")
     # Add ability for downsampling/tiling with tophu

src/dolphin/timeseries.py

@@ -35,6 +35,26 @@ class ReferencePointError(ValueError):
     pass
 
 
+def argmin_index(arr: ArrayLike) -> tuple[int, ...]:
+    """Get the index tuple of the minimum value of the array.
+
+    If multiple occurrences of the minimum value exist, returns
+    the index of the first such occurrence in the flattened array.
+
+    Parameters
+    ----------
+    arr : array_like
+        The input array.
+
+    Returns
+    -------
+    tuple of int
+        The index of the minimum value.
+
+    """
+    return np.unravel_index(np.argmin(arr), np.shape(arr))
+
+
 @jit
 def weighted_lstsq_single(
     A: ArrayLike,
@@ -607,28 +627,33 @@ def correlation_to_variance(correlation: ArrayLike, nlooks: int) -> Array:
 
 def select_reference_point(
     ccl_file_list: Sequence[PathOrStr],
-    amp_dispersion_file: PathOrStr,
+    condition_file: PathOrStr,
     output_dir: Path,
+    condition_func: Callable[[ArrayLike], tuple[int, ...]] = argmin_index,
     block_shape: tuple[int, int] = (512, 512),
     num_threads: int = 5,
 ) -> ReferencePoint:
     """Automatically select a reference point for a stack of unwrapped interferograms.
 
-    Uses the amplitude dispersion and connected component labels, the point is selected
+    Uses the condition file and connected component labels, the point is selected
     which
 
     1. is within intersection of all nonzero connected component labels (always valid)
-    2. has the lowest amplitude dispersion
+    2. has the condition applied to condition file. for example: has the lowest
+       amplitude dispersion
 
     Parameters
     ----------
     ccl_file_list : Sequence[PathOrStr]
         List of connected component label phase files.
-    amp_dispersion_file: PathOrStr
-        Amplitude dispersion, stitched and multilooked to the same size each raster
-        in `ccl_file_list`
+    condition_file: PathOrStr
+        A file with the same size as each raster, like amplitude dispersion or
+        temporal coherence in `ccl_file_list`
     output_dir: Path
         Path to store the computed "conncomp_intersection.tif" raster
+    condition_func: Callable[[ArrayLike, ]]
+        The function to apply to the condition file,
+        for example numpy.argmin which finds the pixel with lowest value
     block_shape: tuple[int, int]
         Size of blocks to read from while processing `ccl_file_list`
         Default = (512, 512)
@@ -693,12 +718,12 @@ def intersect_conncomp(arr: np.ma.MaskedArray, axis: int) -> np.ndarray:
     # Create a mask of pixels with this label
     isin_largest_conncomp = label == largest_idx
 
-    amp_dispersion = io.load_gdal(amp_dispersion_file, masked=True)
+    condition_file_values = io.load_gdal(condition_file, masked=True)
     # Mask out where the conncomps aren't equal to the largest
-    amp_dispersion.mask = amp_dispersion.mask | (~isin_largest_conncomp)
+    condition_file_values.mask = condition_file_values.mask | (~isin_largest_conncomp)
 
-    # Pick the (unmasked) point with the lowest amplitude dispersion
-    ref_row, ref_col = np.unravel_index(np.argmin(amp_dispersion), amp_dispersion.shape)
+    # Pick the (unmasked) point with the condition applied to condition file
+    ref_row, ref_col = condition_func(condition_file_values)
 
     # Cast to `int` to avoid having `np.int64` types
     return ReferencePoint(int(ref_row), int(ref_col))

src/dolphin/unwrap/_snaphu_py.py

@@ -1,6 +1,9 @@
 from __future__ import annotations
 
 from pathlib import Path
+from typing import Optional
+
+from numpy.typing import ArrayLike
 
 from dolphin._log import get_log
 from dolphin._types import Filename
@@ -21,13 +24,13 @@ def unwrap_snaphu_py(
     ntiles: tuple[int, int] = (1, 1),
     tile_overlap: tuple[int, int] = (0, 0),
     nproc: int = 1,
-    mask_file: Filename | None = None,
+    mask_file: Optional[Filename] = None,
     zero_where_masked: bool = False,
-    unw_nodata: float | None = DEFAULT_UNW_NODATA,
-    ccl_nodata: int | None = DEFAULT_CCL_NODATA,
+    unw_nodata: Optional[float] = DEFAULT_UNW_NODATA,
+    ccl_nodata: Optional[int] = DEFAULT_CCL_NODATA,
     init_method: str = "mst",
     cost: str = "smooth",
-    scratchdir: Filename | None = None,
+    scratchdir: Optional[Filename] = None,
 ) -> tuple[Path, Path]:
     """Unwrap an interferogram using at multiple scales using `tophu`.
 
@@ -151,3 +154,74 @@ def unwrap_snaphu_py(
         return _zero_from_mask(unw_filename, cc_filename, mask_file)
 
     return Path(unw_filename), Path(cc_filename)
+
+
+def grow_conncomp_snaphu(
+    unw_filename: Filename,
+    corr_filename: Filename,
+    nlooks: float,
+    mask: Optional[ArrayLike] = None,
+    ccl_nodata: Optional[int] = DEFAULT_CCL_NODATA,
+    cost: str = "smooth",
+    min_conncomp_frac: float = 0.0001,
+    scratchdir: Optional[Filename] = None,
+) -> Path:
+    """Compute connected component labels using SNAPHU.
+
+    Parameters
+    ----------
+    unw_filename : Filename
+        Path to output unwrapped phase file.
+    corr_filename : Filename
+        Path to input correlation file.
+    nlooks : float
+        Effective number of looks used to form the input correlation data.
+    mask : Array, optional
+        binary byte mask array, by default None.
+        Assumes that 1s are valid pixels and 0s are invalid.
+    ccl_nodata : float, optional
+        Nodata value for the connected component labels.
+    cost : str
+        Statistical cost mode.
+        Default = "smooth"
+    min_conncomp_frac : float, optional
+        Minimum size of a single connected component, as a fraction of the total number
+        of pixels in the array. Defaults to 0.0001.
+    scratchdir : Filename, optional
+        If provided, uses a scratch directory to save the intermediate files
+        during unwrapping.
+
+    Returns
+    -------
+    conncomp_path : Filename
+        Path to output connected component label file.
+
+    """
+    import snaphu
+
+    unw_suffix = full_suffix(unw_filename)
+    cc_filename = str(unw_filename).replace(unw_suffix, CONNCOMP_SUFFIX)
+
+    with (
+        snaphu.io.Raster(unw_filename) as unw,
+        snaphu.io.Raster(corr_filename) as corr,
+        snaphu.io.Raster.create(
+            cc_filename,
+            like=unw,
+            nodata=ccl_nodata,
+            dtype="u2",
+            **DEFAULT_TIFF_OPTIONS_RIO,
+        ) as conncomp,
+    ):
+        snaphu.grow_conncomps(
+            unw=unw,
+            corr=corr,
+            nlooks=nlooks,
+            mask=mask,
+            cost=cost,
+            min_conncomp_frac=min_conncomp_frac,
+            scratchdir=scratchdir,
+            conncomp=conncomp,
+        )
+
+    return Path(cc_filename)

src/dolphin/unwrap/_unwrap.py

@@ -21,7 +21,7 @@
     UNW_SUFFIX,
     UNW_SUFFIX_ZEROED,
 )
-from ._snaphu_py import unwrap_snaphu_py
+from ._snaphu_py import grow_conncomp_snaphu, unwrap_snaphu_py
 from ._tophu import multiscale_unwrap
 from ._utils import create_combined_mask, set_nodata_values
 
@@ -37,7 +37,7 @@ def run(
     output_path: Filename,
     *,
     nlooks: float = 5,
-    mask_file: Optional[Filename] = None,
+    mask_filename: Optional[Filename] = None,
     zero_where_masked: bool = False,
     unwrap_method: UnwrapMethod = UnwrapMethod.SNAPHU,
     init_method: str = "mst",
@@ -69,7 +69,7 @@ def run(
         Path to output directory.
     nlooks : int, optional
         Effective number of looks used to form the input correlation data.
-    mask_file : Filename, optional
+    mask_filename : Filename, optional
         Path to binary byte mask file, by default None.
         Assumes that 1s are valid pixels and 0s are invalid.
     zero_where_masked : bool, optional
@@ -160,8 +160,8 @@ def run(
         out_files.append(outf)
     logger.info(f"{len(out_files)} left to unwrap")
 
-    if mask_file:
-        mask_file = Path(mask_file).resolve()
+    if mask_filename:
+        mask_filename = Path(mask_filename).resolve()
 
     # This keeps it from spawning a new process for a single job.
     Executor = ThreadPoolExecutor if max_jobs > 1 else DummyProcessPoolExecutor
@@ -176,7 +176,7 @@ def run(
                 init_method=init_method,
                 cost=cost,
                 unwrap_method=unwrap_method,
-                mask_file=mask_file,
+                mask_filename=mask_filename,
                 zero_where_masked=zero_where_masked,
                 downsample_factor=downsample_factor,
                 ntiles=ntiles,
@@ -197,7 +197,7 @@ def run(
             # We're not passing all the unw files in, so we need to tally up below
             _unw_path, _cc_path = fut.result()
 
-    if zero_where_masked and mask_file is not None:
+    if zero_where_masked and mask_filename is not None:
         all_out_files = [
             Path(str(outf).replace(UNW_SUFFIX, UNW_SUFFIX_ZEROED))
             for outf in all_out_files
@@ -219,7 +219,7 @@ def unwrap(
     corr_filename: Filename,
     unw_filename: Filename,
     nlooks: float,
-    mask_file: Optional[Filename] = None,
+    mask_filename: Optional[Filename] = None,
     zero_where_masked: bool = False,
     ntiles: Union[int, tuple[int, int]] = 1,
     tile_overlap: tuple[int, int] = (0, 0),
@@ -250,7 +250,7 @@ def unwrap(
         Path to output unwrapped phase file.
     nlooks : float
         Effective number of looks used to form the input correlation data.
-    mask_file : Filename, optional
+    mask_filename : Filename, optional
         Path to binary byte mask file, by default None.
         Assumes that 1s are valid pixels and 0s are invalid.
     zero_where_masked : bool, optional
@@ -306,9 +306,9 @@ def unwrap(
 
     Returns
     -------
-    unw_path : Path
+    unw_path : Filename
         Path to output unwrapped phase file.
-    conncomp_path : Path
+    conncomp_path : Filename
         Path to output connected component label file.
 
     """
@@ -320,13 +320,13 @@ def unwrap(
     unwrap_method = UnwrapMethod(unwrap_method)
 
     # Check for a nodata mask
-    if io.get_raster_nodata(ifg_filename) is None or mask_file is None:
+    if io.get_raster_nodata(ifg_filename) is None or mask_filename is None:
         # With no marked `nodata`, just use the passed in mask
-        combined_mask_file = mask_file
+        combined_mask_file = mask_filename
     else:
         combined_mask_file = Path(ifg_filename).with_suffix(".mask.tif")
         create_combined_mask(
-            mask_filename=mask_file,
+            mask_filename=mask_filename,
             image_filename=ifg_filename,
             output_filename=combined_mask_file,
         )
@@ -395,13 +395,15 @@ def unwrap(
             f"Masking pixels with correlation below {interpolation_cor_threshold}"
         )
         coherent_pixel_mask = corr[:] >= interpolation_cor_threshold
+
         logger.info(f"Interpolating {pre_interp_ifg_filename} -> {interp_ifg_filename}")
         modified_ifg = interpolate(
             ifg=ifg,
             weights=coherent_pixel_mask,
             weight_cutoff=interpolation_cor_threshold,
             max_radius=max_radius,
         )
+
         logger.info(f"Writing interpolated output to {interp_ifg_filename}")
         io.write_arr(
             arr=modified_ifg,
@@ -481,4 +483,19 @@ def unwrap(
             options=opts,
         )
 
+        # Regrow connected components after phase modification
+        corr = io.load_gdal(corr_filename)
+        mask = corr[:] > 0
+        # TODO decide whether we want to have the
+        # 'min_conncomp_frac' option in the config
+        conncomp_path = grow_conncomp_snaphu(
+            unw_filename=unw_filename,
+            corr_filename=corr_filename,
+            nlooks=nlooks,
+            mask=mask,
+            ccl_nodata=ccl_nodata,
+            cost=cost,
+            scratchdir=scratchdir,
+        )
+
     return unw_path, conncomp_path

src/dolphin/workflows/_cli_config.py

@@ -39,6 +39,8 @@ def create_config(
     no_unwrap: bool = False,
     no_inversion: bool = False,
     n_parallel_unwrap: int = 1,
+    run_goldstein: bool = False,
+    run_interpolation: bool = False,
     unwrap_method: UnwrapMethod = UnwrapMethod.SNAPHU,
     troposphere_files: Optional[list[str]] = None,
     tropo_date_fmt: str = "%Y%m%d",
@@ -104,6 +106,8 @@ def create_config(
             "n_parallel_jobs": n_parallel_unwrap,
             "run_unwrap": not no_unwrap,
             "zero_where_masked": zero_where_masked,
+            "run_goldstein": run_goldstein,
+            "run_interpolation": run_interpolation,
         },
         timeseries_options={
             "run_inversion": not no_inversion,
@@ -274,6 +278,16 @@ def get_parser(subparser=None, subcommand_name="run"):
         default=1,
         help="Number of interferograms to unwrap in parallel.",
     )
+    unwrap_group.add_argument(
+        "--run-goldstein",
+        action="store_true",
+        help="Run Goldstein filter before unwrapping.",
+    )
+    unwrap_group.add_argument(
+        "--run-interpolation",
+        action="store_true",
+        help="Run interpolation before unwrapping.",
+    )
 
     # Correction options
     correction_group = parser.add_argument_group("Correction options")

src/dolphin/workflows/unwrapping.py

@@ -78,7 +78,7 @@ def run(
         cor_filenames=cor_file_list,
         output_path=output_path,
         nlooks=nlooks,
-        mask_file=output_mask,
+        mask_filename=output_mask,
         zero_where_masked=unwrap_options.zero_where_masked,
         max_jobs=unwrap_options.n_parallel_jobs,
         ntiles=unwrap_options.ntiles,