From ecc1455dff2814780440bc379cb56bb8992e1f1f Mon Sep 17 00:00:00 2001 From: s-sasaki-earthsea-wizard Date: Thu, 13 Aug 2026 14:10:28 +0900 Subject: [PATCH 1/2] Vectorize the per-pixel loop in generate_insar_mask generate_insar_mask built the InSAR mask layer with a pure-Python double loop over every output pixel, calling the scalar SubSwaths.get_sample_sub_swath twice per pixel. On a NISAR L-band frame (RIFG interferogram grid 6840x10581) this costs ~4.2 us/px, ~300 s per product, and makes prepare_insar_hdf5 the largest non-GPU-addressable stage of the InSAR workflow (~560 s, CPU/GPU parity since the cost is interpreter-bound). Replace the inner loop with numpy array operations: - Fetch each sub-swath's per-line [start, end) valid-sample interval array once via the existing SubSwaths.get_valid_samples_array API and evaluate membership as vectorized interval tests (_subswath_numbers), preserving get_sample_sub_swath semantics: out-of-bounds -> 0, first-match-wins ordering, empty-array short-circuit, and no-sub-swath-information -> 1. - Preserve both rounding rules exactly: int(x + 0.5) (truncation toward zero) for the sub-swath lookup via np.trunc, and Python round() (round-half-even) for the exception-mask lookup via np.rint. - Widen the exception-mask bytes to uint32 before the << 16 / << 8 shifts. Under NumPy >= 2.0 (NEP 50) the previous uint8-scalar shifts overflow to 0 and silently drop those bits (#335); the explicit widening produces the intended packing under both promotion regimes. The GDAL row reads and the produced mask values are unchanged; output is bitwise-identical to the previous implementation (verified against a frozen copy of the scalar loop on fixtures covering the rounding edges, empty/missing sub-swath layouts, out-of-swath rows/columns, out-of-bounds secondary indices, and high exception-mask bits, under NumPy 1.26). --- python/packages/nisar/products/insar/utils.py | 175 ++++++++++++------ 1 file changed, 123 insertions(+), 52 deletions(-) diff --git a/python/packages/nisar/products/insar/utils.py b/python/packages/nisar/products/insar/utils.py index 0643bde0c..2d61b79b7 100644 --- a/python/packages/nisar/products/insar/utils.py +++ b/python/packages/nisar/products/insar/utils.py @@ -494,6 +494,62 @@ def generate_dem_rdr(radar_grid_obj, dem_src = None +def _subswath_numbers(subswaths, + intervals, + azi_idx_arr, + rg_idx_arr): + """ + Vectorized equivalent of SubSwaths.get_sample_sub_swath over index + arrays. + + Returns 0 for out-of-swath samples, otherwise the 1-based number of + the first sub-swath whose per-line valid-sample interval + [start, end) contains the sample. An empty interval array claims + every in-bounds sample (matching the scalar API's short-circuit), + and a dataset without sub-swath information assigns 1 everywhere in + bounds. + + Parameters + ---------- + subswaths : isce3.product.SubSwaths + The subswath object of the RSLC + intervals : list of numpy.ndarray + Per-sub-swath [start, end) valid-sample interval arrays, i.e. + [subswaths.get_valid_samples_array(s) for s = 1..num_sub_swaths] + azi_idx_arr : numpy.ndarray + Integer azimuth indices + rg_idx_arr : numpy.ndarray + Integer slant range indices + + Returns + ---------- + numpy.ndarray + int64 sub-swath numbers, same shape as the index arrays + """ + in_bounds = ((azi_idx_arr >= 0) & (azi_idx_arr < subswaths.length) & + (rg_idx_arr >= 0) & (rg_idx_arr < subswaths.width)) + numbers = np.zeros(azi_idx_arr.shape, dtype=np.int64) + if not intervals: + return np.where(in_bounds, np.int64(1), numbers) + + # Clipped so the per-line gather stays legal; out-of-bounds samples + # are excluded through in_bounds + azi_gather = np.clip(azi_idx_arr, 0, subswaths.length - 1) + for number, interval in enumerate(intervals, start=1): + if interval.size == 0: + claimed = in_bounds + else: + claimed = (in_bounds & + (rg_idx_arr >= interval[azi_gather, 0]) & + (rg_idx_arr < interval[azi_gather, 1])) + unassigned = numbers == 0 + numbers[unassigned & claimed] = number + if not unassigned.any(): + break + + return numbers + + def generate_insar_mask(ref_rslc_obj, sec_rslc_obj, ref_rslc_h5_obj, @@ -558,59 +614,74 @@ def _load_exception_mask(h5_obj, rslc_obj, swath): sec_rslc_obj, sec_swath) - mask = [] - for i in azi_idx_arr: - # Check if the azimuth index is within the radar grid - if i >= 0 and i < ref_swath.lines: - range_off = \ - range_offset_band.ReadAsArray(0, - int(i), - ref_swath.samples, - 1) - azimuth_off = \ - azimuth_offset_band.ReadAsArray(0, - int(i), - ref_swath.samples, - 1) - for j in rg_idx_arr: - - # Initialize the all mask ids to be 0 - mask_id = 0 - subswath_mask_id = 0 - ref_input_exception_mask_id = 0 - sec_input_exception_mask_id = 0 - - # Check if the range index is within the swath - if j >= 0 and j < ref_swath.samples: - subswath_mask_id = _compute_subswath_mask_id(int(i),int(j), - azimuth_off[0,int(j)], - range_off[0,int(j)], - ref_subswaths, - sec_subswaths) - - # reference RSLC input exception mask id - ref_input_exception_mask_id = ref_input_exception_mask[int(i),int(j)] << 16 - - # secondary RSLC input exception mask id - sec_i = round(i + azimuth_off[0,int(j)]) - sec_j = round(j + range_off[0,int(j)]) - if ((sec_i >=0 and sec_i < sec_swath.lines) and - (sec_j >=0 and sec_j < sec_swath.samples)): - sec_input_exception_mask_id = sec_input_exception_mask[sec_i,sec_j] << 8 - - # mask id - mask_id = subswath_mask_id | ref_input_exception_mask_id | sec_input_exception_mask_id - - # append the mask id - mask.append(mask_id) - - # The azimuth index is not in the radar grid meaning no subswath mask - else: - mask += [0] * len(rg_idx_arr) + # Fetch each sub-swath's per-line valid-sample interval array once + # (1-based API); the per-sample sub-swath tests below then run as + # numpy array operations instead of two scalar + # SubSwaths.get_sample_sub_swath calls per output pixel + ref_intervals = [ref_subswaths.get_valid_samples_array(s) + for s in range(1, ref_subswaths.num_sub_swaths + 1)] + sec_intervals = [sec_subswaths.get_valid_samples_array(s) + for s in range(1, sec_subswaths.num_sub_swaths + 1)] + + azi_idx_arr = np.asarray(azi_idx_arr, dtype=np.float64) + rg_idx_arr = np.asarray(rg_idx_arr, dtype=np.float64) + + # int() truncates toward zero, as does astype on non-negative and + # negative values alike + rg_idx_int = rg_idx_arr.astype(np.int64) + col_in_swath = (rg_idx_arr >= 0) & (rg_idx_arr < ref_swath.samples) + # Clipped copy so the per-row gathers stay legal; out-of-swath + # columns are zeroed through col_in_swath at the end + rg_gather = np.clip(rg_idx_int, 0, ref_swath.samples - 1) + + mask = np.zeros((len(azi_idx_arr), len(rg_idx_arr)), dtype=np.uint32) + for row, i in enumerate(azi_idx_arr): + # The azimuth index is not in the radar grid meaning no + # subswath mask + if not (0 <= i < ref_swath.lines): + continue + + i_int = int(i) + range_off = range_offset_band.ReadAsArray( + 0, i_int, ref_swath.samples, 1)[0] + azimuth_off = azimuth_offset_band.ReadAsArray( + 0, i_int, ref_swath.samples, 1)[0] + rg_off = range_off[rg_gather] + az_off = azimuth_off[rg_gather] + + # Sub-swath numbers of the reference RSLC and, through the + # nearest neighbor of the geometric coregistration offsets, of + # the secondary RSLC (int(x + 0.5) of the scalar code = + # truncation toward zero) + ref_num = _subswath_numbers( + ref_subswaths, ref_intervals, + np.full(rg_gather.shape, i_int, dtype=np.int64), rg_idx_int) + sec_num = _subswath_numbers( + sec_subswaths, sec_intervals, + np.trunc(i_int + az_off + 0.5).astype(np.int64), + np.trunc(rg_idx_int + rg_off + 0.5).astype(np.int64)) + mask_row = (10 * ref_num + sec_num).astype(np.uint32) + + # Reference RSLC input exception mask bits; widened to uint32 + # before the shift so the packing is safe under NEP 50 scalar + # promotion as well + mask_row |= (ref_input_exception_mask[i_int, rg_gather] + .astype(np.uint32) << 16) + + # Secondary RSLC input exception mask bits; round() of the + # scalar code is round-half-even, as is np.rint + sec_i = np.rint(i + az_off).astype(np.int64) + sec_j = np.rint(rg_idx_arr + rg_off).astype(np.int64) + sec_in_swath = ((sec_i >= 0) & (sec_i < sec_swath.lines) & + (sec_j >= 0) & (sec_j < sec_swath.samples)) + sec_exception = sec_input_exception_mask[ + np.clip(sec_i, 0, sec_swath.lines - 1), + np.clip(sec_j, 0, sec_swath.samples - 1)].astype(np.uint32) << 8 + mask_row |= np.where(sec_in_swath, sec_exception, np.uint32(0)) + + mask[row] = np.where(col_in_swath, mask_row, np.uint32(0)) del ref_input_exception_mask del sec_input_exception_mask - return np.array(mask).reshape( - (len(azi_idx_arr), - len(rg_idx_arr))).astype(np.uint32) \ No newline at end of file + return mask \ No newline at end of file From dbb4e8dbc0b056ca1fd6c7a2f55018b9c80fd3ac Mon Sep 17 00:00:00 2001 From: s-sasaki-earthsea-wizard Date: Wed, 19 Aug 2026 11:10:00 +0900 Subject: [PATCH 2/2] Add regression tests for generate_insar_mask Cover the semantics the vectorized implementation must preserve: - _subswath_numbers against the scalar SubSwaths.get_sample_sub_swath API as oracle: out-of-bounds -> 0, first-match-wins ordering, the empty valid-samples-array short-circuit, and the no-sub-swath-information -> 1 path; - generate_insar_mask against a per-pixel scalar reference (_compute_subswath_mask_id plus Python-int bit packing) on synthetic fixtures with adversarial offsets (exact k + 0.5 half-integers, large out-of-swath pushes), differing reference/secondary dimensions, empty and absent sub-swath layouts, and the missing inputDataExceptionMask dataset path; - exact uint32 packing of exception-mask bytes with the MSB set, which the previous uint8-scalar << 16 / << 8 shifts silently drop under NumPy >= 2 scalar promotion (#335); - the rounding asymmetry between the sub-swath lookup (int(x + 0.5), truncation toward zero) and the exception-mask lookup (round-half-even), including a negative secondary index where truncation toward zero and floor diverge. The seven behavior tests pass unchanged against the pre-vectorization scalar implementation under NumPy 1.26, confirming they encode the existing semantics rather than the new implementation's. --- tests/python/packages/CMakeLists.txt | 1 + .../packages/nisar/products/insar/utils.py | 349 ++++++++++++++++++ 2 files changed, 350 insertions(+) create mode 100644 tests/python/packages/nisar/products/insar/utils.py diff --git a/tests/python/packages/CMakeLists.txt b/tests/python/packages/CMakeLists.txt index ddf0e7402..7ef3ce0c1 100644 --- a/tests/python/packages/CMakeLists.txt +++ b/tests/python/packages/CMakeLists.txt @@ -45,6 +45,7 @@ nisar/noise/noise_estimation_from_raw.py nisar/pointing/doppler_lut_from_raw.py #nisar/pointing/el_null_range_from_raw_ant.py nisar/products/granule_id.py +nisar/products/insar/utils.py nisar/products/readers/antenna_parser.py nisar/products/readers/attitude.py nisar/products/readers/GSLC.py diff --git a/tests/python/packages/nisar/products/insar/utils.py b/tests/python/packages/nisar/products/insar/utils.py new file mode 100644 index 000000000..a47de8ad6 --- /dev/null +++ b/tests/python/packages/nisar/products/insar/utils.py @@ -0,0 +1,349 @@ +"""Tests for the mask-generation helpers in nisar.products.insar.utils. + +The vectorized generate_insar_mask() and _subswath_numbers() must +reproduce the semantics of the original scalar per-pixel loop exactly: + +- SubSwaths.get_sample_sub_swath equivalence: out-of-bounds -> 0, + first-match-wins sub-swath ordering, an empty valid-samples array + claims every in-bounds sample, and a dataset without sub-swath + information assigns 1 everywhere in bounds; +- the two distinct rounding rules: the sub-swath lookup uses + int(x + 0.5) (truncation toward zero) while the exception-mask + lookup uses round() (round-half-even); +- rows/columns outside the reference swath produce whole-pixel 0, and + secondary indices outside the secondary swath drop the secondary + exception bits while keeping the sub-swath contribution; +- the inputDataExceptionMask bytes are packed into the uint32 mask as + (ref << 16) | (sec << 8), which must hold for mask bytes >= 2**(8-k) + under NumPy 2 scalar promotion (NEP 50) as well. +""" +import numpy as np +import h5py +import pytest +from osgeo import gdal + +import isce3 +from nisar.products.insar.utils import (_compute_subswath_mask_id, + _subswath_numbers, + generate_insar_mask) + +SWATH_PATH = "/science/LSAR/RSLC/swaths" + +# Unique names for the in-memory h5py files +_h5_counter = 0 + + +class FakeSwath: + """Duck-typed stand-in for the Swath metadata generate_insar_mask + reads (lines, samples, sub_swaths).""" + + def __init__(self, lines, samples, subswaths): + self.lines = lines + self.samples = samples + self._subswaths = subswaths + + def sub_swaths(self): + return self._subswaths + + +class FakeSLC: + SwathPath = SWATH_PATH + + def __init__(self, swath): + self._swath = swath + + def getSwathMetadata(self, freq): + return self._swath + + +def make_h5(freq, exception_mask): + """In-memory h5py file; exception_mask=None omits the dataset.""" + global _h5_counter + _h5_counter += 1 + f = h5py.File(f"insar-mask-test-{_h5_counter}", "w", + driver="core", backing_store=False) + if exception_mask is not None: + f.create_dataset( + f"{SWATH_PATH}/frequency{freq}/inputDataExceptionMask", + data=exception_mask) + return f + + +def make_offset_raster(path, data): + drv = gdal.GetDriverByName("ENVI") + ds = drv.Create(str(path), data.shape[1], data.shape[0], 1, + gdal.GDT_Float64) + ds.GetRasterBand(1).WriteArray(data) + ds.FlushCache() + ds = None + + +def make_subswaths(rng, lines, samples, n_sub, empty_at=(), no_info=False): + """Random sub-swath layout with a few fully-invalid lines; empty_at + sub-swath numbers get an empty valid-samples array, no_info builds a + SubSwaths without any sub-swath information.""" + if no_info: + return isce3.product.SubSwaths(lines, samples, []) + arrays = [] + for s in range(1, n_sub + 1): + if s in empty_at: + arrays.append(np.empty((0, 0), dtype=np.int32)) + continue + start = rng.integers(0, samples, size=lines) + width = rng.integers(0, samples // 2 + 1, size=lines) + end = np.minimum(start + width, samples) + invalid = rng.random(lines) < 0.1 + end[invalid] = start[invalid] + arrays.append(np.stack([start, end], axis=1).astype(np.int32)) + return isce3.product.SubSwaths(lines, samples, arrays) + + +def build_offsets(rng, lines, samples, scale): + """Offset field with adversarial values: smooth random, offsets that + land (index + offset) exactly on k + 0.5 half-integers where the two + rounding rules diverge, and large pushes outside the secondary + swath.""" + off = rng.normal(0.0, scale, size=(lines, samples)) + jj = np.arange(samples, dtype=np.float64) + half_rows = rng.choice(lines, size=max(1, lines // 5), replace=False) + for r in half_rows: + targets = rng.integers(-3, samples + 3, + size=samples).astype(np.float64) + 0.5 + sel = rng.random(samples) < 0.3 + off[r, sel] = (targets - jj)[sel] + blow = rng.random((lines, samples)) < 0.02 + off[blow] = rng.choice([-1.0, 1.0], size=blow.sum()) * (samples + lines) + return off + + +def scalar_reference_mask(ref_swath, sec_swath, ref_exception_mask, + sec_exception_mask, range_off, azimuth_off, + azi_idx_arr, rg_idx_arr): + """Per-pixel reference implementation of the intended mask + semantics: _compute_subswath_mask_id for the sub-swath digits and + Python-int bit packing of the exception-mask bytes (immune to fixed + width scalar promotion).""" + ref_subswaths = ref_swath.sub_swaths() + sec_subswaths = sec_swath.sub_swaths() + mask = np.zeros((len(azi_idx_arr), len(rg_idx_arr)), dtype=np.uint32) + for row, i in enumerate(azi_idx_arr): + if not (0 <= i < ref_swath.lines): + continue + for col, j in enumerate(rg_idx_arr): + if not (0 <= j < ref_swath.samples): + continue + az_off = azimuth_off[int(i), int(j)] + rg_off = range_off[int(i), int(j)] + mask_id = _compute_subswath_mask_id( + int(i), int(j), az_off, rg_off, + ref_subswaths, sec_subswaths) + mask_id |= int(ref_exception_mask[int(i), int(j)]) << 16 + sec_i = round(i + az_off) + sec_j = round(j + rg_off) + if (0 <= sec_i < sec_swath.lines and + 0 <= sec_j < sec_swath.samples): + mask_id |= int(sec_exception_mask[sec_i, sec_j]) << 8 + mask[row, col] = mask_id + return mask + + +class TestSubswathNumbers: + """_subswath_numbers against the scalar pybind oracle + SubSwaths.get_sample_sub_swath.""" + + @pytest.mark.parametrize("n_sub,empty_at,no_info", [ + (3, (), False), + (1, (), False), + (3, (2,), False), + (0, (), True), + ]) + def test_matches_scalar_api(self, n_sub, empty_at, no_info): + rng = np.random.default_rng(12345 + n_sub + 100 * no_info) + lines, samples = 40, 56 + subswaths = make_subswaths(rng, lines, samples, n_sub, + empty_at=empty_at, no_info=no_info) + intervals = [subswaths.get_valid_samples_array(s) + for s in range(1, subswaths.num_sub_swaths + 1)] + + # every index pair from 4 outside the swath on either side, + # plus random scattered pairs + azi = np.arange(-4, lines + 4, dtype=np.int64) + rg = np.arange(-4, samples + 4, dtype=np.int64) + azi_grid, rg_grid = np.meshgrid(azi, rg, indexing="ij") + + actual = _subswath_numbers(subswaths, intervals, + azi_grid, rg_grid) + + expected = np.array( + [[subswaths.get_sample_sub_swath(int(a), int(r)) for r in rg] + for a in azi], dtype=np.int64) + np.testing.assert_array_equal(actual, expected) + + +class TestGenerateInsarMask: + """generate_insar_mask against the per-pixel scalar reference.""" + + @pytest.mark.parametrize( + "name,seed,ref_dims,sec_dims,n_sub,empty_at,no_info,no_masks," + "off_scale", + [ + ("random_3sub", 1, (40, 56), (37, 59), 3, (), False, False, 2.5), + ("empty_mid_subswath", 2, (32, 48), (32, 48), 3, (2,), False, + False, 2.5), + ("no_subswath_info", 3, (32, 48), (30, 44), 0, (), True, False, + 2.5), + ("no_exception_masks", 4, (32, 48), (32, 48), 2, (), False, True, + 2.5), + ("large_offsets", 5, (36, 42), (22, 30), 2, (), False, False, + 25.0), + ]) + def test_matches_scalar_reference(self, tmp_path, name, seed, ref_dims, + sec_dims, n_sub, empty_at, no_info, + no_masks, off_scale): + rng = np.random.default_rng(seed) + ref_lines, ref_samples = ref_dims + sec_lines, sec_samples = sec_dims + + ref_swath = FakeSwath( + ref_lines, ref_samples, + make_subswaths(rng, ref_lines, ref_samples, n_sub, + empty_at=empty_at, no_info=no_info)) + sec_swath = FakeSwath( + sec_lines, sec_samples, + make_subswaths(rng, sec_lines, sec_samples, n_sub, + no_info=no_info)) + + if no_masks: + ref_exc = np.zeros((ref_lines, ref_samples), dtype=np.uint8) + sec_exc = np.zeros((sec_lines, sec_samples), dtype=np.uint8) + ref_h5 = make_h5("A", None) + sec_h5 = make_h5("A", None) + else: + ref_exc = rng.integers(0, 256, size=(ref_lines, ref_samples), + dtype=np.uint8) + sec_exc = rng.integers(0, 256, size=(sec_lines, sec_samples), + dtype=np.uint8) + ref_h5 = make_h5("A", ref_exc) + sec_h5 = make_h5("A", sec_exc) + + range_off = build_offsets(rng, ref_lines, ref_samples, off_scale) + azimuth_off = build_offsets(rng, ref_lines, ref_samples, off_scale) + rg_off_path = tmp_path / f"{name}_range.off" + az_off_path = tmp_path / f"{name}_azimuth.off" + make_offset_raster(rg_off_path, range_off) + make_offset_raster(az_off_path, azimuth_off) + + # integral-float index arrays extending past the swath on both + # sides, matching the np.round(...) arrays the callers build + azi_idx = np.round(np.linspace(-3, ref_lines + 3, ref_lines + 8)) + rg_idx = np.round(np.linspace(-3, ref_samples + 3, ref_samples + 8)) + + actual = generate_insar_mask( + FakeSLC(ref_swath), FakeSLC(sec_swath), ref_h5, sec_h5, + str(rg_off_path), str(az_off_path), "A", azi_idx, rg_idx) + expected = scalar_reference_mask( + ref_swath, sec_swath, ref_exc, sec_exc, range_off, + azimuth_off, azi_idx, rg_idx) + + ref_h5.close() + sec_h5.close() + assert actual.dtype == np.uint32 + np.testing.assert_array_equal(actual, expected) + + def test_exception_mask_bit_packing(self, tmp_path): + """Exact uint32 packing, including bytes with the MSB set whose + << 16 / << 8 shifts overflow a fixed-width uint8 scalar + (regression test for silently dropped exception bits).""" + lines, samples = 4, 6 + # single sub-swath covering every sample -> sub-swath digits 11 + full = np.tile(np.array([[0, samples]], dtype=np.int32), + (lines, 1)) + swath = FakeSwath( + lines, samples, + isce3.product.SubSwaths(lines, samples, [full])) + + ref_exc = np.zeros((lines, samples), dtype=np.uint8) + sec_exc = np.zeros((lines, samples), dtype=np.uint8) + ref_exc[1, 2] = 0xAB + sec_exc[1, 2] = 0xCD + ref_exc[2, 3] = 0x80 + sec_exc[2, 3] = 0xFF + ref_h5 = make_h5("A", ref_exc) + sec_h5 = make_h5("A", sec_exc) + + zeros = np.zeros((lines, samples), dtype=np.float64) + rg_off_path = tmp_path / "packing_range.off" + az_off_path = tmp_path / "packing_azimuth.off" + make_offset_raster(rg_off_path, zeros) + make_offset_raster(az_off_path, zeros) + + idx_azi = np.arange(lines, dtype=np.float64) + idx_rg = np.arange(samples, dtype=np.float64) + mask = generate_insar_mask( + FakeSLC(swath), FakeSLC(swath), ref_h5, sec_h5, + str(rg_off_path), str(az_off_path), "A", idx_azi, idx_rg) + ref_h5.close() + sec_h5.close() + + assert mask[1, 2] == (0xAB << 16) | (0xCD << 8) | 11 + assert mask[2, 3] == (0x80 << 16) | (0xFF << 8) | 11 + assert mask[0, 0] == 11 + + def test_rounding_rules_diverge(self, tmp_path): + """The sub-swath lookup truncates int(x + 0.5) while the + exception-mask lookup rounds half to even; an exact +0.5 range + offset at an even column exercises both: the sub-swath lookup + reads column j + 1 while the exception mask reads column j.""" + lines, samples = 3, 8 + # sub-swath 1 = columns [0, 4), sub-swath 2 = columns [4, 8) + s1 = np.tile(np.array([[0, 4]], dtype=np.int32), (lines, 1)) + s2 = np.tile(np.array([[4, 8]], dtype=np.int32), (lines, 1)) + swath = FakeSwath( + lines, samples, + isce3.product.SubSwaths(lines, samples, [s1, s2])) + + sec_exc = np.zeros((lines, samples), dtype=np.uint8) + sec_exc[1, 2] = 0x11 # read by round-half-even (2.5 -> 2) + sec_exc[1, 3] = 0x22 # NOT read: int(2.5 + 0.5) = 3 is the + # sub-swath lookup only + ref_h5 = make_h5("A", np.zeros((lines, samples), dtype=np.uint8)) + sec_h5 = make_h5("A", sec_exc) + + rg_off = np.zeros((lines, samples), dtype=np.float64) + rg_off[1, 2] = 3.5 # column 2 -> secondary range 5.5: + # sub-swath int(6.0) = 6 -> sub-swath 2, + # exception round(5.5) = 6 (half-even) + rg_off[1, 4] = -1.5 # column 4 -> secondary range 2.5: + # sub-swath int(3.0) = 3 -> sub-swath 1, + # exception round(2.5) = 2 + rg_off[1, 6] = -6.7 # column 6 -> secondary range -0.7: + # int(-0.2) truncates toward zero to 0 + # (in bounds, sub-swath 1) while + # exception round(-0.7) = -1 is out of + # bounds and drops the secondary bits + az_off = np.zeros((lines, samples), dtype=np.float64) + rg_off_path = tmp_path / "rounding_range.off" + az_off_path = tmp_path / "rounding_azimuth.off" + make_offset_raster(rg_off_path, rg_off) + make_offset_raster(az_off_path, az_off) + + idx_azi = np.arange(lines, dtype=np.float64) + idx_rg = np.arange(samples, dtype=np.float64) + mask = generate_insar_mask( + FakeSLC(swath), FakeSLC(swath), ref_h5, sec_h5, + str(rg_off_path), str(az_off_path), "A", idx_azi, idx_rg) + ref_h5.close() + sec_h5.close() + + # column 2: ref sub-swath 1, sec sub-swath 2, sec exception + # byte from round-half-even column 6 (zero) + assert mask[1, 2] == 12 + # column 4: ref sub-swath 2, sec sub-swath 1, sec exception + # byte 0x11 from round-half-even column 2 + assert mask[1, 4] == (0x11 << 8) | 21 + # column 6: ref sub-swath 2, sec column int(-0.2) = 0 -> + # sub-swath 1 (a floor would give -1 -> 0), no exception bits + assert mask[1, 6] == 21 + # zero-offset columns keep matching digits and no exception bits + assert mask[1, 0] == 11 + assert mask[1, 5] == 22