Performance tune for prepare_insar_hdf5, resample_slc_v2, and rubbersheet

[lijun99]

This PR aims to speed up three modules in insar workflows, prepare_insar_hdf5, resample_slc_v2, and rubbersheet. I put them together for easier testing. Please let me know if separate PRs are preferred.

1. prepare_insar_hdf5
   Replace the per-pixel Python loop (and its per-pixel calls into SubSwaths.get_sample_sub_swath) with vectorized numpy operations, reading offset raster rows only as needed to keep memory bounded. Tests show ~3.5-10x speedup on the prepare_insar_hdf5 step, and output verified bit-for-bit identical to the original.

2. resample_slc_v2
   Replace the legacy HDF5 reader with HDF5OptimizedReader used in unwrap.py, crossmul.py, and others. This cuts the resample step's SLC I/O read time by ~5-10x. Output is verified bit-for-bit identical to the original.

3. rubbersheet
   Update rubbersheet polyfit critical_value from 0.1 to Baarda's convention 3.29.
   0.1 (keeping only 8%) is far tighter than the standard w-test threshold, causing excessive outlier-removal iterations without a corresponding accuracy benefit. 3.29 (alpha=0.001) is the conventional geodetic data-snooping critical value. Recommended range 1.69-3.29.
   Here is a test on a scene from Afar, Ethiopia,

crit_value Iterations Time % of data removed Removed points with corr_peak > 0.7
0.1 (current) 23,187 72.35s 58.0% 68.3% (mostly good matches!)
1.0 3,197 16.30s 8.0% 3.4%
2.0 3,111 15.76s 7.8% 2.8%
3.29 (Baarda's standard) 3,008 14.66s 7.5% 2.8%

It turns out that 0.1 removes 68.3% of good pixels for polyfit, not only makes the fitting less representative, but also leads to long iterations to finish.
In tests, choosing 3.29 leads to ~5-10x speedup for the polyfit. (Note that the rubbersheet procedure has heavy I/Os, the actual I/O performance depends on types of storage systems.) There is a small difference, ~0.2 rad across the scene, from a ramp, since the fitting parameters are slightly different.

BTW, I am also working on a gpu phase unwrapping. It shows some early success; but needs more test. Please give it a try when you have time.

[lijun99]

The test failures seem to be pre-existing. Please take a look at this PR, @hfattahi @bhawkins @xhuang-jpl.

[Tyler-g-hudson]

@lijun99 Could you rebase to the develop branch? That should resolve these test failures.

[Tyler-g-hudson]

It looks like Linux Latest Debug and (particularly) MacOS Latest Debug have large numbers of test failures outside of the ones that @bhawkins fixed

[Tyler-g-hudson]

We're getting a weird mix of test failures here.

Minimum Python Supported:

• ❌ Failed test.cxx.isce3.geometry.geometry.geometry
• ✅ No SEGFAULT in test.python.pybind.geometry.pntintersect

Linux Latest Release:

• ❌ Failed test.cxx.isce3.geometry.geometry.geometry
• ❌ SEGFAULT in test.python.pybind.geometry.pntintersect

Linux Latest Debug:

• ✅ No fail in test.cxx.isce3.geometry.geometry.geometry
• ✅ No SEGFAULT in test.python.pybind.geometry.pntintersect
• ❌ Lots of aborted subprocesses in python plus one in cxx
• ❌ Strange errors with files and headers in some python workflows

macOS Latest Release:

• ✅ No fail in test.cxx.isce3.geometry.geometry.geometry
• ❌ SEGFAULT in test.python.pybind.geometry.pntintersect

macOS Latest Debug:

• ❌ All Python tests fail with errors importing the pybinds:

../../../../../../install/packages/isce3/__init__.py:2: in <module>
    from .ext import extisce3
../../../../../../install/packages/isce3/ext/__init__.py:3: in <module>
    from . import isce3 as extisce3
E   ImportError: Internal error while parsing type signature (2)

EDIT: These appear to be universal to recent PR's. I will open an issue about it.

[Tyler-g-hudson]

This pattern of errors is recorded now in #368. I don't see anything in the content of the errors that would be special to these changes, but the resample_slc and resample_slc_v2 Python workflows are effected in Linux Latest Debug so there's a small chance that an error in this PR for Linux Latest Debug would be hidden by those errors.

@hfattahi how do you feel about proceeding with review with this in mind? I think such a hidden error is very improbable, so my position is that we can proceed.

[Tyler-g-hudson]

The CI error issue is fixed by #366 - please rebase when that PR merges and rebase again.

[lijun99]

Great, all test passed! Thank you, @bhawkins and @Tyler-g-hudson !

[Tyler-g-hudson]

With regards to the Resample update, this looks good!

[hfattahi]

@lijun99 thanks for this PR. I am ok with Resample and h5 prep changes but couple of us think that the polyfitting changes requires intense testing. To make progress faster, I suggest to decouple this PR to three separate PRs. Resample can be merged quickly. we will review h5 prep and we will need to discuss and test polyfitting further.

[lijun99]

I have spilt the PR into three separate parts.

2. resample_slc_v2, resubmitted as Use optimized HDF5 reader for resample_slc_v2 #372.

3. rubbersheet, per discussion with @hfattahi, changing the default critical_value should be cautious. Let's do more tests before changing it. But from my research into the topic, 0.1 doesn't make sense to me.

4. prepare_insar_hdf5, I noticed Add regression tests for generate_insar_mask #359 from @s-sasaki-earthsea-wizard, which is on the same line. I will check his approach before deciding whether I will need to submit a new one.

Thanks!

[s-sasaki-earthsea-wizard]

@lijun99 Thanks for the ping — and sorry for the overlap: I should have spotted #358 before filing #359 (I had been tracking this loop via #354 / #335 since Aug 13, but missed your PR when I filed mine two days after yours).

Since the two utils.py changes are independent implementations of the same fix, I ran yours through the benchmark harness I built for #354 to make the comparison concrete. Setup: utils.py from the #358 head (bdb5974), from the #359 head (dbb4e8d), and pristine develop (blob 0643bde0c) as baseline, all against the same compiled isce3 (develop, numpy 1.26.4, py3.12).

Timing — synthetic reproducer from #354 (6840 x 10581 = 72.4 Mpx RIFG-like grid, three sub-swaths, no NISAR granules needed):

utils.py variant	wall time (median of 3)	vs scalar	peak RSS
develop (scalar per-pixel loop)	572.6 s (single run)	1x	—
#358 (_get_sample_subswath_grid)	15.8 s	~36x	9.9 GB
#359 (_subswath_numbers)	10.0 s	~57x	1.6 GB

The output mask CRC is identical (0x295ba19b) across all three variants and all runs on this fixture, so on correctness of the hot path and on the size of the win we are in full agreement — both rewrites eliminate the per-pixel pybind calls, which is where ~all the time goes. The residual differences come from temporaries: #358 materializes several full-grid float64/int64 arrays via meshgrid (hence the peak RSS), while #359 processes row-wise. For reference, on a quieter host I previously measured the scalar loop at 383.3 s with the same ~57x ratio for #359, so the ratios are stable across host-load conditions.

A minor contract-fidelity note on the #358 version. I also ran both implementations against a self-contained scalar oracle that calls SubSwaths.get_sample_sub_swath directly (adversarial fixtures: exact k + 0.5 offsets where the two rounding rules diverge, empty sub-swath arrays, out-of-swath pushes). Both match the oracle on every file-loadable input — including an RSLC without sub-swath datasets, which Serialization.h deserializes as num_sub_swaths == 1 with an empty valid-samples array, a state both PRs handle correctly via the empty-array short-circuit. The one divergence I found is num_sub_swaths == 0: there the scalar API returns 1 for in-bounds samples (SubSwaths::getSampleSubSwath: "If the dataset does not have sub-swaths information, consider samples valid and belonging to the first sub-swath"), while _get_sample_subswath_grid's per-sub-swath loop never executes and leaves 0, dropping the sub-swath digits. To be clear about reachability: this state cannot come from RSLC deserialization (the numSubSwaths setter even throws for n <= 0), only from programmatic construction such as the pybind list constructor SubSwaths(length, width, []) — so it is a fidelity gap against the documented scalar API rather than a production-path bug, and it matters mostly for tests and library consumers that build SubSwaths objects directly. If you keep your implementation it is a one-line guard to make the helper a full drop-in.

Tests. #359 also carries a pytest suite (~350 lines, tests/python/packages/nisar/products/insar/utils.py) that validates the vectorized mask against the scalar pybind API as oracle, including the two distinct rounding rules (int(x + 0.5) truncation for the sub-swath lookup vs round() half-even for the exception-mask lookup), the empty-interval short-circuit, the num_sub_swaths == 0 fallback above, and the uint32 bit packing that #335 is about.

Happy to resolve the overlap either way: if you and the maintainers prefer the #358 implementation for the prepare part, I can rebase #359 down to a tests-only PR — the suite validates either implementation (modulo the fallback above). Otherwise #359 already carries implementation + tests. Whichever way, the cross-check script and logs are public in my benchmark repo (s-sasaki-earthsea-wizard/isce3-benchmark, scripts/crosscheck_insar_mask_oracle.py + artifacts/pr358-crosscheck-20260827/).

[s-sasaki-earthsea-wizard]

For completeness: the real-data leg of the comparison above (NISAR frame; bitwise + peak-RSS results for both implementations) is now posted in #359 (comment) — keeping the data in one place since that is where the prepare-part decision lands.