From 3f2980c7064a7af368272565ce2f91c1d2c8f7da Mon Sep 17 00:00:00 2001 From: Brett Johnston Date: Mon, 10 Aug 2026 12:50:06 -0700 Subject: [PATCH 1/4] Initial Image Alias Mitigation --- sarkit_processing/sicd_alias_mitigation.py | 410 +++++++++++++++++++++ tests/conftest.py | 34 ++ tests/test_sicd_alias_mitigation.py | 237 ++++++++++++ 3 files changed, 681 insertions(+) create mode 100644 sarkit_processing/sicd_alias_mitigation.py create mode 100644 tests/test_sicd_alias_mitigation.py diff --git a/sarkit_processing/sicd_alias_mitigation.py b/sarkit_processing/sicd_alias_mitigation.py new file mode 100644 index 0000000..edcfcf4 --- /dev/null +++ b/sarkit_processing/sicd_alias_mitigation.py @@ -0,0 +1,410 @@ +import argparse +import math +from typing import Sequence + +import lxml.etree +import numba +import numpy as np +import numpy.polynomial.polynomial as npp +import numpy.typing as npt +import sarkit.sicd as sksicd +import scipy.fft as spfft +import scipy.ndimage + +from sarkit_processing import sicd_pixel_type +from sarkit_processing.sicd_deskew import apply_phase_poly, get_deskew_phase_poly + +try: + from smart_open import open +except ImportError: + pass + + +def shift_poly_axis(coeffs: npt.ArrayLike, shift: float, axis: int) -> npt.NDArray: + """ + Return new coefficients so that: + + OrigPoly(..., x_axis + shift, ...) == NewPoly(..., x_axis, ...) + + coeffs[k0, k1, ...] is the coefficient for: + x0**k0 * x1**k1 * ... + + Parameters + ---------- + coeffs : array-like + N-dimensional polynomial coefficient tensor. + shift : float + Shift applied to the selected input axis. + axis : int + Axis whose variable is shifted. + + Returns + ------- + ndarray + Shifted coefficient tensor. + """ + coeffs = np.asarray(coeffs) + out = np.zeros_like(coeffs, dtype=np.result_type(coeffs, shift, float)) + + degree = coeffs.shape[axis] - 1 + + # Move target polynomial axis to front: shape (degree+1, ...) + c = np.moveaxis(coeffs, axis, 0) + o = np.moveaxis(out, axis, 0) + + powers = shift ** np.arange(degree + 1) + + for old_power in range(degree + 1): + # (x + shift)^old_power = sum new_power terms + for new_power in range(old_power + 1): + o[new_power] += ( + c[old_power] + * math.comb(old_power, new_power) + * powers[old_power - new_power] + ) + + return out + + +def polyscale2d(coeffs: npt.NDArray, scale_x: float, scale_y: float) -> npt.NDArray: + """ + Returns new polynomial with scaled coordinate axes so that: + + NewPoly(x, y) == OrigPoly(scale_x * x, scale_y * y). + + Parameters + ---------- + coeffs : ndarray + 2-D polynomial coefficients, where coeffs[i, j] is the coefficient + of x**i * y**j. + scale_x : float + Scale factor applied to the x coordinate. + scale_y : float + Scale factor applied to the y coordinate. + + Returns + ------- + ndarray + Scaled coefficient tensor. + """ + ix = np.arange(coeffs.shape[0])[:, np.newaxis] + iy = np.arange(coeffs.shape[1])[np.newaxis, :] + scale = (scale_x**ix) * (scale_y**iy) + return coeffs * scale + + +@numba.njit(parallel=True) +def _apply_row_phase_poly(array, phase_poly, col_0, col_ss): + """Apply a separate 1D phase polynomial to each row.""" + out = np.empty_like(array) + for rowidx in numba.prange(array.shape[0]): + coeffs = phase_poly[rowidx] + for colidx in range(array.shape[1]): + col_val = col_0 + colidx * col_ss + phase_val = coeffs[-1] + for ndx in range(coeffs.shape[0] - 1, 0, -1): + phase_val = phase_val * col_val + coeffs[ndx - 1] + out[rowidx, colidx] = array[rowidx, colidx] * np.exp( + 1j * 2 * np.pi * phase_val + ) + + return out + + +def _power(arr): + return np.abs(arr) ** 2 + + +def _powersum(arr): + return np.sum(np.abs(arr) ** 2) + + +def _normalize(arr): + out = np.zeros_like(arr) + np.divide(arr, np.abs(arr), out=out, where=arr != 0) + return out + + +def _unit(arr): + return arr / np.linalg.norm(arr, axis=-1) + + +def _apply_spectral_phase( + arr, phase_polys, fft_sizes, output_shape, col_0=0, col_ss=1.0 +): + arr = spfft.fftshift(spfft.fft(arr, axis=1, n=fft_sizes[1]), axes=1) + arr = spfft.fftshift(spfft.fft(arr.T, axis=1, n=fft_sizes[0]), axes=1) + arr = _apply_row_phase_poly(arr, phase_polys, col_0, col_ss) + arr = spfft.ifft(spfft.ifftshift(arr, axes=1), axis=1, n=fft_sizes[0])[ + :, : output_shape[0] + ] + arr = spfft.ifft(spfft.ifftshift(arr.T, axes=1), axis=1, n=fft_sizes[1])[ + :, : output_shape[1] + ] + + return arr + + +def _focus_alias(arr, sicd_xmltree, zone, prf_override): + sicdew = sksicd.ElementWrapper(sicd_xmltree.getroot()) + if prf_override is not None: + ipp_poly = [0, prf_override] + else: + ipp_poly = sicdew["Timeline"]["IPP"]["Set"][0]["IPPPoly"] # TODO which poly? + + coa_nom = npp.polyval2d(0, 0, sicdew["Grid"]["TimeCOAPoly"]) + prf_nom = npp.polyval(coa_nom, npp.polyder(ipp_poly)) + krg_ctr_nom = sicdew["Grid"]["Row"]["KCtr"] + alias_freq_nom = zone * prf_nom + range_rate_nom = alias_freq_nom / krg_ctr_nom + arp_poly = sicdew["Position"]["ARPPoly"] + pos_nom = npp.polyval(coa_nom, arp_poly) + vel_nom = npp.polyval(coa_nom, npp.polyder(arp_poly)) + scp = sicdew["GeoData"]["SCP"]["ECF"] + tan_pa_rate = np.linalg.norm( + np.cross(vel_nom, _unit(scp - pos_nom)) + ) / np.linalg.norm(scp - pos_nom) + phase_slope_poly = np.array( + [-alias_freq_nom / tan_pa_rate, range_rate_nom / tan_pa_rate] + ) + fft_size_0 = spfft.next_fast_len(arr.shape[0]) + fft_size_1 = spfft.next_fast_len(arr.shape[1]) + kaz_ss = 1.0 / (sicdew["Grid"]["Col"]["SS"] * fft_size_1) + kaz = (np.arange(fft_size_1) - fft_size_1 // 2) * kaz_ss + krg_ss = 1.0 / (sicdew["Grid"]["Row"]["SS"] * fft_size_0) + phase_vs_krg_poly = phase_slope_poly * kaz[:, np.newaxis] + phase_vs_samp_poly = polyscale2d( + shift_poly_axis(phase_vs_krg_poly, krg_ctr_nom, 1), 1.0, krg_ss / krg_ctr_nom + ) + fft_sizes = (fft_size_0, fft_size_1) + focus = _apply_spectral_phase(arr, phase_vs_samp_poly, fft_sizes, fft_sizes) + + return focus, fft_sizes, phase_vs_samp_poly + + +def _remove_alias(arr, sicd_xmltree, zone, thresh, dilation, prf_override): + # TODO remove intermediate variables + focus_alias, fft_sizes, phase_polys = _focus_alias( + arr, sicd_xmltree, zone, prf_override + ) + normalize = _normalize(arr) + focus_alias_norm = _apply_spectral_phase( + normalize, phase_polys, fft_sizes, fft_sizes + ) + power = _power(focus_alias_norm) + + # Gives array of bool + over_thresh = power > thresh # TODO numba + if dilation > 0: # TODO verify + dilated_mask = scipy.ndimage.binary_dilation( + over_thresh, + structure=np.ones((2 * dilation + 1, 2 * dilation + 1), dtype=bool), + ) + else: + dilated_mask = over_thresh + + removed_power = np.sum(np.abs(focus_alias) ** 2, where=dilated_mask) + under_thresh_vals = np.where(~dilated_mask, focus_alias, 0) + focus_img = _apply_spectral_phase( + under_thresh_vals, -phase_polys, fft_sizes, arr.shape + ) + kept_data_frac = np.mean(~dilated_mask) + + return focus_img, removed_power, kept_data_frac + + +def _mitigate(arr, sicd_xmltree, zones, thresh, dilation, prf_override=None): + removed_power = np.float32(0) + kept_data_frac = np.float64(1) + for zone in zones: + arr, zone_removed_power, zone_kept_data_frac = _remove_alias( + arr, sicd_xmltree, zone, thresh, dilation, prf_override + ) + removed_power += zone_removed_power + kept_data_frac *= zone_kept_data_frac + + return arr, removed_power, kept_data_frac + + +def _minimumpower(left, right): + assert left.shape == right.shape + assert left.dtype == right.dtype + return np.where(np.abs(left) <= np.abs(right), left, right) + + +def _trim_spectrum(arr, sicd_xmltree, axis): + sicdew = sksicd.ElementWrapper(sicd_xmltree.getroot()) + axis_idx = {"Row": 0, "Col": 1}[axis] + fft_size = spfft.next_fast_len(arr.shape[axis_idx]) + num_keep = int( + np.ceil( + sicdew["Grid"][axis]["ImpRespBW"] * sicdew["Grid"][axis]["SS"] * fft_size + ) + ) + phase_poly = get_deskew_phase_poly(sicd_xmltree, axis=axis) + arr, sicd_xmltree = apply_phase_poly(arr, phase_poly, sicd_xmltree) + if axis == "Row": + arr = arr.T + fftf0 = spfft.fft(arr, axis=1, n=fft_size) + fftf0[:, num_keep // 2 : fft_size - num_keep // 2] = 0 + ffti0 = spfft.ifft(fftf0, axis=1, n=fft_size)[:, : arr.shape[-1]] + if axis == "Row": + ffti0 = ffti0.T + arr, sicd_xmltree = apply_phase_poly(ffti0, -phase_poly, sicd_xmltree) + return arr, sicd_xmltree + + +def _trim_spectrum_2d(arr, sicd_xmltree): + arr, sicd_xmltree = _trim_spectrum(arr, sicd_xmltree, "Col") + arr, sicd_xmltree = _trim_spectrum(arr, sicd_xmltree, "Row") + + return arr + + +def prf_alias_removal( + image: npt.NDArray, + sicd_xmltree: lxml.etree.ElementTree, + num_iters: int, + zones: Sequence[float], + *, + threshold: float = 9.0, + dilate: int = 1, + prf: float | None = None, +) -> tuple[npt.NDArray, float, float]: + """ + Remove PRF alias energy + + Parameters + ---------- + image : ndarray + The complex array to apply the algorithm. + sicd_xmltree : lxml.etree.ElementTree + SICD XML ElementTree + num_iters : int + Number of iterations of cleaning. + zones : sequence of float + Sequence of alias zones that should be cleaned. The zone values need not be integers. + threshold : float, optional + Threshold to use when creating NIFT mask (units of stddev) (Default: 9.0). + dilate : int, optional + Number of pixels to dilate the masks (Default: 1). + prf : float or None, optional + Pulse repetition frequency (Hz). (Default: None). + If ``None``, the SICD's Timeline/IPP parameters are used. + + Returns + ------- + out : ndarray + The array with alias energy removed. + removed_pwr_frac : float + The fraction of the power that was removed. + removed_data_frac: float + The fraction of the bins that were removed. + """ + out = before = image.copy() + ps = _powersum(image) + removed_power = np.float32(0) + kept_data_frac = np.float64(1) + for _ in range(num_iters): + out, iter_removed_power, iter_kept_data_frac = _mitigate( + out, sicd_xmltree, zones, threshold, dilate, prf + ) + removed_power += iter_removed_power + kept_data_frac *= iter_kept_data_frac + # TODO: Verify trim + out = _minimumpower(out, before) + out = _trim_spectrum_2d(out, sicd_xmltree) + + removed_pwr_frac = removed_power / ps + removed_data_frac = 1.0 - kept_data_frac + + return out, removed_pwr_frac, removed_data_frac + + +def main(args=None): + parser = argparse.ArgumentParser( + description="Clean up PRF alias energy from a SICD." + ) + parser.add_argument("input_sicd_filename") + parser.add_argument("output_sicd_filename") + parser.add_argument( + "--threshold", + type=float, + default=9.0, + help=( + "Threshold to use when creating NIFT mask (units of stddev) " + "(default: %(default)g)" + ), + ) + parser.add_argument( + "--num-iters", + type=int, + default=2, + help="Number of iterations of cleaning to do (default: %(default)d)", + ) + parser.add_argument( + "--zones", + type=float, + default=[-2, -1, 1, 2], + nargs="+", + help="PRF multipliers to search for alias energy (default: %(default)s)", + ) + parser.add_argument( + "--symmetric", + action="store_true", + help="Search positive and negative of zones (default: listed zones only, not their inverse)", + ) + parser.add_argument( + "--dilate", + type=int, + default=1, + help="Number of pixels to dilate the masks (default: %(default)d)", + ) + parser.add_argument( + "--prf", type=float, default=None, help="Override the PRF with a fixed value" + ) + config = parser.parse_args(args) + + zones = config.zones.copy() + if config.symmetric: + zones += [-zone for zone in config.zones] + zones = sorted(set(zones)) + + with ( + open(config.input_sicd_filename, "rb") as file, + sksicd.NitfReader(file) as reader, + ): + xmltree = reader.metadata.xmltree + image = reader.read_image() + + image, xmltree = sicd_pixel_type.sicd_as_re32f_im32f(image, xmltree) + sicdew = sksicd.ElementWrapper(xmltree.getroot()) + assert sicdew["Grid"]["Row"]["Sgn"] == sicdew["Grid"]["Col"]["Sgn"] + if sicdew["Grid"]["Row"]["Sgn"] == 1: + image = np.conjugate(image) + sicdew["Grid"]["Row"]["Sgn"] = sicdew["Grid"]["Col"]["Sgn"] = -1 + + image, removed_pwr_frac, removed_data_frac = prf_alias_removal( + image.astype("complex64"), + xmltree, + config.num_iters, + zones, + threshold=config.threshold, + dilate=config.dilate, + prf=config.prf, + ) + metadata = reader.metadata + metadata.xmltree = xmltree + + with ( + open(config.output_sicd_filename, "wb") as file, + sksicd.NitfWriter(file, reader.metadata) as writer, + ): + writer.write_image(image) + + print("removed_pwr_frac = ", removed_pwr_frac) + print("removed_data_frac = ", removed_data_frac) + + +if __name__ == "__main__": + main() diff --git a/tests/conftest.py b/tests/conftest.py index 1e110f4..b6cf3e8 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -47,3 +47,37 @@ def example_sicd(tmp_path_factory): with open(tmp_sicd, "wb") as f, sksicd.NitfWriter(f, sicd_meta) as w: w.write_image(_random_array((nrows, ncols), dtype)) yield tmp_sicd + + +@pytest.fixture(scope="session") +def example_sicd_alias(tmp_path_factory): + rng = np.random.default_rng(seed=20260810) + sicd_etree = lxml.etree.parse(good_sicd_xml_path) + tmp_sicd = ( + tmp_path_factory.mktemp("data") / good_sicd_xml_path.with_suffix(".sicd").name + ) + shape = ( + int(sicd_etree.findtext("{*}ImageData/{*}NumRows")), + int(sicd_etree.findtext("{*}ImageData/{*}NumCols")), + ) + + image = ( + rng.uniform(-1.0, 1.0, size=shape) + 1j * rng.uniform(-1.0, 1.0, size=shape) + ).astype(np.complex64) + + fft1 = np.fft.fftshift(np.fft.fft2(image)) + fft1[900:1060, 400:610] += 1000 + image = np.fft.ifft2(np.fft.ifftshift(fft1)) + + sec = {"security": {"clas": "U"}} + sicd_meta = sksicd.NitfMetadata( + xmltree=sicd_etree, + file_header_part={"ostaid": "nowhere"} | sec, + im_subheader_part={"isorce": "this sensor"} | sec, + de_subheader_part=sec, + ) + + with open(tmp_sicd, "wb") as f, sksicd.NitfWriter(f, sicd_meta) as w: + w.write_image(image) + + yield tmp_sicd diff --git a/tests/test_sicd_alias_mitigation.py b/tests/test_sicd_alias_mitigation.py new file mode 100644 index 0000000..e967907 --- /dev/null +++ b/tests/test_sicd_alias_mitigation.py @@ -0,0 +1,237 @@ +import shutil +import subprocess +import sys + +import numpy as np +import numpy.polynomial.polynomial as npp +import pytest +import sarkit.sicd as sksicd + +import sarkit_processing.sicd_alias_mitigation as spsam +import tests.utils + + +@pytest.mark.parametrize("shift", [-2.5, -1.0, 0.0, 0.75, 3.0]) +def test_shift_poly_axis(shift): + rng = np.random.default_rng() + coeffs = rng.normal(size=6) + shifted = spsam.shift_poly_axis(coeffs, shift, 0) + + xvals = rng.normal(size=100) + + np.testing.assert_allclose( + npp.polyval(xvals + shift, coeffs), + npp.polyval(xvals, shifted), + rtol=1e-11, + atol=1e-11, + ) + + +@pytest.mark.parametrize("axis", [0, 1]) +def test_2d_shift_poly_axis(axis): + rng = np.random.default_rng() + coeffs = rng.normal(size=(5, 4)) + shift = 1.3 + + shifted = spsam.shift_poly_axis(coeffs, shift, axis) + + xvals = rng.normal(size=50) + yvals = rng.normal(size=50) + + if axis == 0: + expected = npp.polyval2d(xvals + shift, yvals, coeffs) + else: + expected = npp.polyval2d(xvals, yvals + shift, coeffs) + + actual = npp.polyval2d(xvals, yvals, shifted) + + np.testing.assert_allclose(actual, expected, rtol=1e-11, atol=1e-11) + + +@pytest.mark.parametrize( + "scale_x, scale_y", + [ + (1.0, 1.0), + (2.0, 3.0), + (0.5, 4.0), + (-2.0, 3.0), + (-1.5, -0.5), + ], +) +def test_polyscale2d(scale_x, scale_y): + rng = np.random.default_rng() + coeffs = rng.normal(size=(5, 4)) + scaled = spsam.polyscale2d(coeffs, scale_x, scale_y) + + xvals = rng.normal(size=100) + yvals = rng.normal(size=100) + + expected = npp.polyval2d(scale_x * xvals, scale_y * yvals, coeffs) + actual = npp.polyval2d(xvals, yvals, scaled) + + np.testing.assert_allclose(actual, expected, rtol=1e-11, atol=1e-11) + + +def test_main(tmp_path, example_sicd_alias): + output_file = tmp_path / "cleaned_example.sicd" + + subprocess.check_call( + [ + sys.executable, + "-m", + "sarkit_processing.sicd_alias_mitigation", + str(example_sicd_alias), + str(output_file), + "--threshold", + "7.0", + "--num-iters", + "3", + "--zones", + "-1.0", + "2.0", + "--symmetric", + "--dilate", + "3", + ], + cwd=tmp_path, + ) + + assert output_file.is_file() + + +def test_sicd_alias_mitigation_smart_open(tmp_path, example_sicd_alias): + output_file = tmp_path / "out.sicd" + + shutil.copyfile(example_sicd_alias, tmp_path / example_sicd_alias.name) + with tests.utils.static_http_server(tmp_path) as server_url: + subprocess.check_call( + [ + sys.executable, + "-m", + "sarkit_processing.sicd_alias_mitigation", + f"{server_url}/{example_sicd_alias.name}", + output_file, + ], + ) + + assert output_file.exists() + + +def test_sicd_alias_mitigation_thresh(example_sicd_alias): + with ( + open(example_sicd_alias, "rb") as file, + sksicd.NitfReader(file) as reader, + ): + xmltree = reader.metadata.xmltree + image = reader.read_image() + + num_iters = 2 + zones = [-1.0, 1.0] + _, removed_pwr_frac, removed_data_frac = spsam.prf_alias_removal( + image.astype("complex64"), xmltree, num_iters, zones, threshold=8.0 + ) + + _, removed_pwr_frac_1, removed_data_frac_1 = spsam.prf_alias_removal( + image.astype("complex64"), xmltree, num_iters, zones, threshold=7.0 + ) + + assert removed_pwr_frac < removed_pwr_frac_1 + assert removed_data_frac < removed_data_frac_1 + + +def test_sicd_alias_mitigation_dilate(example_sicd_alias): + with ( + open(example_sicd_alias, "rb") as file, + sksicd.NitfReader(file) as reader, + ): + xmltree = reader.metadata.xmltree + image = reader.read_image() + + num_iters = 2 + zones = [-1.0, 1.0] + _, removed_pwr_frac, removed_data_frac = spsam.prf_alias_removal( + image.astype("complex64"), + xmltree, + num_iters, + zones, + threshold=7.0, + ) + + _, removed_pwr_frac_1, removed_data_frac_1 = spsam.prf_alias_removal( + image.astype("complex64"), + xmltree, + num_iters, + zones, + threshold=7.0, + dilate=3, + ) + + assert removed_pwr_frac < removed_pwr_frac_1 + assert removed_data_frac < removed_data_frac_1 + + +def test_sicd_alias_mitigation_iters(example_sicd_alias): + with ( + open(example_sicd_alias, "rb") as file, + sksicd.NitfReader(file) as reader, + ): + xmltree = reader.metadata.xmltree + image = reader.read_image() + + zones = [-1.0, 1.0] + _, removed_pwr_frac, removed_data_frac = spsam.prf_alias_removal( + image.astype("complex64"), + xmltree, + 2, + zones, + threshold=4.0, + ) + + _, removed_pwr_frac_1, removed_data_frac_1 = spsam.prf_alias_removal( + image.astype("complex64"), + xmltree, + 3, + zones, + threshold=4.0, + ) + + assert removed_pwr_frac_1 > removed_pwr_frac + assert removed_data_frac_1 > removed_data_frac + + +def test_sicd_alias_mitigation_prf(example_sicd_alias): + with ( + open(example_sicd_alias, "rb") as file, + sksicd.NitfReader(file) as reader, + ): + xmltree = reader.metadata.xmltree + image = reader.read_image() + + sicdew = sksicd.ElementWrapper(xmltree.getroot()) + + num_iters = 2 + zones = [-1.0, 1.0] + _, removed_pwr_frac, removed_data_frac = spsam.prf_alias_removal( + image.astype("complex64"), + xmltree, + num_iters, + zones, + threshold=4.0, + ) + + time_poly = sicdew["Grid"]["TimeCOAPoly"] + ipp_poly = sicdew["Timeline"]["IPP"]["Set"][0]["IPPPoly"] + coa = npp.polyval2d(0, 0, time_poly) + prf = npp.polyval(coa, npp.polyder(ipp_poly)) + del sicdew["Timeline"]["IPP"] + _, removed_pwr_frac_1, removed_data_frac_1 = spsam.prf_alias_removal( + image.astype("complex64"), + xmltree, + num_iters, + zones, + threshold=4.0, + prf=prf, + ) + + assert removed_pwr_frac == removed_pwr_frac_1 + assert removed_data_frac == removed_data_frac_1 From 20d6127e99f399b1e8fae00881f29621343b4b65 Mon Sep 17 00:00:00 2001 From: Brett Johnston Date: Tue, 11 Aug 2026 16:52:21 -0700 Subject: [PATCH 2/4] Add alias mitigation to the single entry point --- CHANGELOG.md | 2 + sarkit_processing/__main__.py | 5 + sarkit_processing/sicd_alias_mitigation.py | 181 +++++++++++---------- tests/test_sicd_alias_mitigation.py | 25 +-- 4 files changed, 108 insertions(+), 105 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f9250af..14d181c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- `sicd_alias` CLI subcommand ## [0.1.0] - 2026-08-03 diff --git a/sarkit_processing/__main__.py b/sarkit_processing/__main__.py index e09933d..ea1d820 100644 --- a/sarkit_processing/__main__.py +++ b/sarkit_processing/__main__.py @@ -3,6 +3,7 @@ import sarkit_processing._coords import sarkit_processing._sicd_chip +import sarkit_processing.sicd_alias_mitigation from sarkit_processing import _cli @@ -23,6 +24,10 @@ def add_subcommand(name, sc: _cli.Subcommand) -> None: add_subcommand("coords", sarkit_processing._coords.CoordsSubcommand()) add_subcommand("sicd_chip", sarkit_processing._sicd_chip.SicdChipSubcommand()) + add_subcommand( + "sicd_alias", + sarkit_processing.sicd_alias_mitigation._SicdAliasSubcommand(), + ) config = parser.parse_args(args) return config.command_handler(config) diff --git a/sarkit_processing/sicd_alias_mitigation.py b/sarkit_processing/sicd_alias_mitigation.py index edcfcf4..2281e12 100644 --- a/sarkit_processing/sicd_alias_mitigation.py +++ b/sarkit_processing/sicd_alias_mitigation.py @@ -1,4 +1,3 @@ -import argparse import math from typing import Sequence @@ -11,7 +10,7 @@ import scipy.fft as spfft import scipy.ndimage -from sarkit_processing import sicd_pixel_type +from sarkit_processing import _cli, sicd_pixel_type from sarkit_processing.sicd_deskew import apply_phase_poly, get_deskew_phase_poly try: @@ -20,7 +19,7 @@ pass -def shift_poly_axis(coeffs: npt.ArrayLike, shift: float, axis: int) -> npt.NDArray: +def _shift_poly_axis(coeffs: npt.ArrayLike, shift: float, axis: int) -> npt.NDArray: """ Return new coefficients so that: @@ -66,7 +65,7 @@ def shift_poly_axis(coeffs: npt.ArrayLike, shift: float, axis: int) -> npt.NDArr return out -def polyscale2d(coeffs: npt.NDArray, scale_x: float, scale_y: float) -> npt.NDArray: +def _polyscale2d(coeffs: npt.NDArray, scale_x: float, scale_y: float) -> npt.NDArray: """ Returns new polynomial with scaled coordinate axes so that: @@ -173,8 +172,8 @@ def _focus_alias(arr, sicd_xmltree, zone, prf_override): kaz = (np.arange(fft_size_1) - fft_size_1 // 2) * kaz_ss krg_ss = 1.0 / (sicdew["Grid"]["Row"]["SS"] * fft_size_0) phase_vs_krg_poly = phase_slope_poly * kaz[:, np.newaxis] - phase_vs_samp_poly = polyscale2d( - shift_poly_axis(phase_vs_krg_poly, krg_ctr_nom, 1), 1.0, krg_ss / krg_ctr_nom + phase_vs_samp_poly = _polyscale2d( + _shift_poly_axis(phase_vs_krg_poly, krg_ctr_nom, 1), 1.0, krg_ss / krg_ctr_nom ) fft_sizes = (fft_size_0, fft_size_1) focus = _apply_spectral_phase(arr, phase_vs_samp_poly, fft_sizes, fft_sizes) @@ -321,90 +320,94 @@ def prf_alias_removal( return out, removed_pwr_frac, removed_data_frac -def main(args=None): - parser = argparse.ArgumentParser( - description="Clean up PRF alias energy from a SICD." - ) - parser.add_argument("input_sicd_filename") - parser.add_argument("output_sicd_filename") - parser.add_argument( - "--threshold", - type=float, - default=9.0, - help=( - "Threshold to use when creating NIFT mask (units of stddev) " - "(default: %(default)g)" - ), - ) - parser.add_argument( - "--num-iters", - type=int, - default=2, - help="Number of iterations of cleaning to do (default: %(default)d)", - ) - parser.add_argument( - "--zones", - type=float, - default=[-2, -1, 1, 2], - nargs="+", - help="PRF multipliers to search for alias energy (default: %(default)s)", - ) - parser.add_argument( - "--symmetric", - action="store_true", - help="Search positive and negative of zones (default: listed zones only, not their inverse)", - ) - parser.add_argument( - "--dilate", - type=int, - default=1, - help="Number of pixels to dilate the masks (default: %(default)d)", - ) - parser.add_argument( - "--prf", type=float, default=None, help="Override the PRF with a fixed value" - ) - config = parser.parse_args(args) - - zones = config.zones.copy() - if config.symmetric: - zones += [-zone for zone in config.zones] - zones = sorted(set(zones)) - - with ( - open(config.input_sicd_filename, "rb") as file, - sksicd.NitfReader(file) as reader, - ): - xmltree = reader.metadata.xmltree - image = reader.read_image() - - image, xmltree = sicd_pixel_type.sicd_as_re32f_im32f(image, xmltree) - sicdew = sksicd.ElementWrapper(xmltree.getroot()) - assert sicdew["Grid"]["Row"]["Sgn"] == sicdew["Grid"]["Col"]["Sgn"] - if sicdew["Grid"]["Row"]["Sgn"] == 1: - image = np.conjugate(image) - sicdew["Grid"]["Row"]["Sgn"] = sicdew["Grid"]["Col"]["Sgn"] = -1 - - image, removed_pwr_frac, removed_data_frac = prf_alias_removal( - image.astype("complex64"), - xmltree, - config.num_iters, - zones, - threshold=config.threshold, - dilate=config.dilate, - prf=config.prf, - ) - metadata = reader.metadata - metadata.xmltree = xmltree +class _SicdAliasSubcommand(_cli.Subcommand): + def get_argument_parser_kwargs(self): + return dict( + description="Clean up PRF alias energy from a SICD", + ) - with ( - open(config.output_sicd_filename, "wb") as file, - sksicd.NitfWriter(file, reader.metadata) as writer, - ): - writer.write_image(image) + def add_arguments(self, parser): + parser.add_argument("input_sicd_filename") + parser.add_argument("output_sicd_filename") + parser.add_argument( + "--threshold", + type=float, + default=9.0, + help=( + "Threshold to use when creating NIFT mask (units of stddev) " + "(default: %(default)g)" + ), + ) + parser.add_argument( + "--num-iters", + type=int, + default=2, + help="Number of iterations of cleaning to do (default: %(default)d)", + ) + parser.add_argument( + "--zones", + type=float, + default=[-2, -1, 1, 2], + nargs="+", + help="PRF multipliers to search for alias energy (default: %(default)s)", + ) + parser.add_argument( + "--symmetric", + action="store_true", + help="Search positive and negative of zones (default: listed zones only, not their inverse)", + ) + parser.add_argument( + "--dilate", + type=int, + default=1, + help="Number of pixels to dilate the masks (default: %(default)d)", + ) + parser.add_argument( + "--prf", + type=float, + default=None, + help="Override the PRF with a fixed value", + ) + + def run_command(self, config): + zones = config.zones.copy() + if config.symmetric: + zones += [-zone for zone in config.zones] + zones = sorted(set(zones)) + + with ( + open(config.input_sicd_filename, "rb") as file, + sksicd.NitfReader(file) as reader, + ): + xmltree = reader.metadata.xmltree + image = reader.read_image() + + image, xmltree = sicd_pixel_type.sicd_as_re32f_im32f(image, xmltree) + sicdew = sksicd.ElementWrapper(xmltree.getroot()) + assert sicdew["Grid"]["Row"]["Sgn"] == sicdew["Grid"]["Col"]["Sgn"] + if sicdew["Grid"]["Row"]["Sgn"] == 1: + image = np.conjugate(image) + sicdew["Grid"]["Row"]["Sgn"] = sicdew["Grid"]["Col"]["Sgn"] = -1 + + image, removed_pwr_frac, removed_data_frac = prf_alias_removal( + image.astype("complex64"), + xmltree, + config.num_iters, + zones, + threshold=config.threshold, + dilate=config.dilate, + prf=config.prf, + ) + metadata = reader.metadata + metadata.xmltree = xmltree - print("removed_pwr_frac = ", removed_pwr_frac) - print("removed_data_frac = ", removed_data_frac) + with ( + open(config.output_sicd_filename, "wb") as file, + sksicd.NitfWriter(file, reader.metadata) as writer, + ): + writer.write_image(image) + print("removed_pwr_frac = ", removed_pwr_frac) + print("removed_data_frac = ", removed_data_frac) -if __name__ == "__main__": - main() + return 0 diff --git a/tests/test_sicd_alias_mitigation.py b/tests/test_sicd_alias_mitigation.py index e967907..f491c67 100644 --- a/tests/test_sicd_alias_mitigation.py +++ b/tests/test_sicd_alias_mitigation.py @@ -1,12 +1,11 @@ import shutil -import subprocess -import sys import numpy as np import numpy.polynomial.polynomial as npp import pytest import sarkit.sicd as sksicd +import sarkit_processing.__main__ import sarkit_processing.sicd_alias_mitigation as spsam import tests.utils @@ -15,7 +14,7 @@ def test_shift_poly_axis(shift): rng = np.random.default_rng() coeffs = rng.normal(size=6) - shifted = spsam.shift_poly_axis(coeffs, shift, 0) + shifted = spsam._shift_poly_axis(coeffs, shift, 0) xvals = rng.normal(size=100) @@ -33,7 +32,7 @@ def test_2d_shift_poly_axis(axis): coeffs = rng.normal(size=(5, 4)) shift = 1.3 - shifted = spsam.shift_poly_axis(coeffs, shift, axis) + shifted = spsam._shift_poly_axis(coeffs, shift, axis) xvals = rng.normal(size=50) yvals = rng.normal(size=50) @@ -61,7 +60,7 @@ def test_2d_shift_poly_axis(axis): def test_polyscale2d(scale_x, scale_y): rng = np.random.default_rng() coeffs = rng.normal(size=(5, 4)) - scaled = spsam.polyscale2d(coeffs, scale_x, scale_y) + scaled = spsam._polyscale2d(coeffs, scale_x, scale_y) xvals = rng.normal(size=100) yvals = rng.normal(size=100) @@ -74,12 +73,9 @@ def test_polyscale2d(scale_x, scale_y): def test_main(tmp_path, example_sicd_alias): output_file = tmp_path / "cleaned_example.sicd" - - subprocess.check_call( + sarkit_processing.__main__.main( [ - sys.executable, - "-m", - "sarkit_processing.sicd_alias_mitigation", + "sicd_alias", str(example_sicd_alias), str(output_file), "--threshold", @@ -93,7 +89,6 @@ def test_main(tmp_path, example_sicd_alias): "--dilate", "3", ], - cwd=tmp_path, ) assert output_file.is_file() @@ -104,13 +99,11 @@ def test_sicd_alias_mitigation_smart_open(tmp_path, example_sicd_alias): shutil.copyfile(example_sicd_alias, tmp_path / example_sicd_alias.name) with tests.utils.static_http_server(tmp_path) as server_url: - subprocess.check_call( + sarkit_processing.__main__.main( [ - sys.executable, - "-m", - "sarkit_processing.sicd_alias_mitigation", + "sicd_alias", f"{server_url}/{example_sicd_alias.name}", - output_file, + str(output_file), ], ) From e5faab3550b5f0ebffec2b1c2f9316b04080f316 Mon Sep 17 00:00:00 2001 From: Brett Johnston Date: Thu, 13 Aug 2026 14:25:25 -0700 Subject: [PATCH 3/4] Add new processing node to the output sicd xml --- sarkit_processing/__main__.py | 4 +- ...sicd_alias_mitigation.py => sicd_alias.py} | 39 ++++++++++++++----- sarkit_processing/sicd_area_plane.py | 3 +- ...alias_mitigation.py => test_sicd_alias.py} | 38 +++++++++--------- 4 files changed, 53 insertions(+), 31 deletions(-) rename sarkit_processing/{sicd_alias_mitigation.py => sicd_alias.py} (91%) rename tests/{test_sicd_alias_mitigation.py => test_sicd_alias.py} (80%) diff --git a/sarkit_processing/__main__.py b/sarkit_processing/__main__.py index ea1d820..ca340d2 100644 --- a/sarkit_processing/__main__.py +++ b/sarkit_processing/__main__.py @@ -3,7 +3,7 @@ import sarkit_processing._coords import sarkit_processing._sicd_chip -import sarkit_processing.sicd_alias_mitigation +import sarkit_processing.sicd_alias from sarkit_processing import _cli @@ -26,7 +26,7 @@ def add_subcommand(name, sc: _cli.Subcommand) -> None: add_subcommand("sicd_chip", sarkit_processing._sicd_chip.SicdChipSubcommand()) add_subcommand( "sicd_alias", - sarkit_processing.sicd_alias_mitigation._SicdAliasSubcommand(), + sarkit_processing.sicd_alias._SicdAliasSubcommand(), ) config = parser.parse_args(args) diff --git a/sarkit_processing/sicd_alias_mitigation.py b/sarkit_processing/sicd_alias.py similarity index 91% rename from sarkit_processing/sicd_alias_mitigation.py rename to sarkit_processing/sicd_alias.py index 2281e12..9e60311 100644 --- a/sarkit_processing/sicd_alias_mitigation.py +++ b/sarkit_processing/sicd_alias.py @@ -1,3 +1,5 @@ +import datetime +import importlib.metadata import math from typing import Sequence @@ -194,7 +196,7 @@ def _remove_alias(arr, sicd_xmltree, zone, thresh, dilation, prf_override): # Gives array of bool over_thresh = power > thresh # TODO numba - if dilation > 0: # TODO verify + if dilation > 0: dilated_mask = scipy.ndimage.binary_dilation( over_thresh, structure=np.ones((2 * dilation + 1, 2 * dilation + 1), dtype=bool), @@ -268,7 +270,7 @@ def prf_alias_removal( *, threshold: float = 9.0, dilate: int = 1, - prf: float | None = None, + prf_override: float | None = None, ) -> tuple[npt.NDArray, float, float]: """ Remove PRF alias energy @@ -287,7 +289,7 @@ def prf_alias_removal( Threshold to use when creating NIFT mask (units of stddev) (Default: 9.0). dilate : int, optional Number of pixels to dilate the masks (Default: 1). - prf : float or None, optional + prf_override : float or None, optional Pulse repetition frequency (Hz). (Default: None). If ``None``, the SICD's Timeline/IPP parameters are used. @@ -306,7 +308,7 @@ def prf_alias_removal( kept_data_frac = np.float64(1) for _ in range(num_iters): out, iter_removed_power, iter_kept_data_frac = _mitigate( - out, sicd_xmltree, zones, threshold, dilate, prf + out, sicd_xmltree, zones, threshold, dilate, prf_override ) removed_power += iter_removed_power kept_data_frac *= iter_kept_data_frac @@ -363,7 +365,7 @@ def add_arguments(self, parser): help="Number of pixels to dilate the masks (default: %(default)d)", ) parser.add_argument( - "--prf", + "--prf-override", type=float, default=None, help="Override the PRF with a fixed value", @@ -382,13 +384,13 @@ def run_command(self, config): xmltree = reader.metadata.xmltree image = reader.read_image() - image, xmltree = sicd_pixel_type.sicd_as_re32f_im32f(image, xmltree) sicdew = sksicd.ElementWrapper(xmltree.getroot()) assert sicdew["Grid"]["Row"]["Sgn"] == sicdew["Grid"]["Col"]["Sgn"] if sicdew["Grid"]["Row"]["Sgn"] == 1: - image = np.conjugate(image) - sicdew["Grid"]["Row"]["Sgn"] = sicdew["Grid"]["Col"]["Sgn"] = -1 + # TODO: Handle SICDs with Sgn=+1 + raise NotImplementedError("SICDs with Sgn=+1 are not supported") + image, xmltree = sicd_pixel_type.sicd_as_re32f_im32f(image, xmltree) image, removed_pwr_frac, removed_data_frac = prf_alias_removal( image.astype("complex64"), xmltree, @@ -396,11 +398,30 @@ def run_command(self, config): zones, threshold=config.threshold, dilate=config.dilate, - prf=config.prf, + prf_override=config.prf_override, ) metadata = reader.metadata metadata.xmltree = xmltree + # Add processing node to describe alias mitigation results + now = datetime.datetime.now(datetime.timezone.utc).isoformat() + sicdew["ImageFormation"].add( + "Processing", + { + "Type": f"{__package__} {importlib.metadata.version(__package__)} | sicd_alias @ {now}", + "Applied": "true", + "Parameter": ( + ("num_iters", f"{config.num_iters}"), + ("zones", " ".join(map(str, zones))), + ("threshold", f"{config.threshold}"), + ("dilate", f"{config.dilate}"), + ("prf_override", f"{config.prf_override}"), + ("removed_pwr_frac", f"{removed_pwr_frac}"), + ("removed_data_frac", f"{removed_data_frac}"), + ), + }, + ) + with ( open(config.output_sicd_filename, "wb") as file, sksicd.NitfWriter(file, reader.metadata) as writer, diff --git a/sarkit_processing/sicd_area_plane.py b/sarkit_processing/sicd_area_plane.py index e3a7a63..ee44d44 100644 --- a/sarkit_processing/sicd_area_plane.py +++ b/sarkit_processing/sicd_area_plane.py @@ -171,8 +171,9 @@ def compute_suitable_rc_plane(sicd_xmltree): eigen_val, eigen_vec = np.linalg.eig( ground_resolution_proj @ ground_resolution_proj.T ) + eigen_val = eigen_val.real + eigen_vec = eigen_vec.real max_val = np.argmax(eigen_val) - scene_coord_angle = np.arctan2(eigen_vec[max_val][1], eigen_vec[max_val][0]) semi_major = np.sqrt(np.abs(eigen_val[max_val])) semi_minor = np.sqrt(np.abs(eigen_val[(max_val + 1) % 2])) diff --git a/tests/test_sicd_alias_mitigation.py b/tests/test_sicd_alias.py similarity index 80% rename from tests/test_sicd_alias_mitigation.py rename to tests/test_sicd_alias.py index f491c67..ca45437 100644 --- a/tests/test_sicd_alias_mitigation.py +++ b/tests/test_sicd_alias.py @@ -6,7 +6,7 @@ import sarkit.sicd as sksicd import sarkit_processing.__main__ -import sarkit_processing.sicd_alias_mitigation as spsam +import sarkit_processing.sicd_alias as spsa import tests.utils @@ -14,7 +14,7 @@ def test_shift_poly_axis(shift): rng = np.random.default_rng() coeffs = rng.normal(size=6) - shifted = spsam._shift_poly_axis(coeffs, shift, 0) + shifted = spsa._shift_poly_axis(coeffs, shift, 0) xvals = rng.normal(size=100) @@ -32,7 +32,7 @@ def test_2d_shift_poly_axis(axis): coeffs = rng.normal(size=(5, 4)) shift = 1.3 - shifted = spsam._shift_poly_axis(coeffs, shift, axis) + shifted = spsa._shift_poly_axis(coeffs, shift, axis) xvals = rng.normal(size=50) yvals = rng.normal(size=50) @@ -60,7 +60,7 @@ def test_2d_shift_poly_axis(axis): def test_polyscale2d(scale_x, scale_y): rng = np.random.default_rng() coeffs = rng.normal(size=(5, 4)) - scaled = spsam._polyscale2d(coeffs, scale_x, scale_y) + scaled = spsa._polyscale2d(coeffs, scale_x, scale_y) xvals = rng.normal(size=100) yvals = rng.normal(size=100) @@ -94,7 +94,7 @@ def test_main(tmp_path, example_sicd_alias): assert output_file.is_file() -def test_sicd_alias_mitigation_smart_open(tmp_path, example_sicd_alias): +def test_sicd_alias_smart_open(tmp_path, example_sicd_alias): output_file = tmp_path / "out.sicd" shutil.copyfile(example_sicd_alias, tmp_path / example_sicd_alias.name) @@ -110,7 +110,7 @@ def test_sicd_alias_mitigation_smart_open(tmp_path, example_sicd_alias): assert output_file.exists() -def test_sicd_alias_mitigation_thresh(example_sicd_alias): +def test_sicd_alias_thresh(example_sicd_alias): with ( open(example_sicd_alias, "rb") as file, sksicd.NitfReader(file) as reader, @@ -120,11 +120,11 @@ def test_sicd_alias_mitigation_thresh(example_sicd_alias): num_iters = 2 zones = [-1.0, 1.0] - _, removed_pwr_frac, removed_data_frac = spsam.prf_alias_removal( + _, removed_pwr_frac, removed_data_frac = spsa.prf_alias_removal( image.astype("complex64"), xmltree, num_iters, zones, threshold=8.0 ) - _, removed_pwr_frac_1, removed_data_frac_1 = spsam.prf_alias_removal( + _, removed_pwr_frac_1, removed_data_frac_1 = spsa.prf_alias_removal( image.astype("complex64"), xmltree, num_iters, zones, threshold=7.0 ) @@ -132,7 +132,7 @@ def test_sicd_alias_mitigation_thresh(example_sicd_alias): assert removed_data_frac < removed_data_frac_1 -def test_sicd_alias_mitigation_dilate(example_sicd_alias): +def test_sicd_alias_dilate(example_sicd_alias): with ( open(example_sicd_alias, "rb") as file, sksicd.NitfReader(file) as reader, @@ -142,7 +142,7 @@ def test_sicd_alias_mitigation_dilate(example_sicd_alias): num_iters = 2 zones = [-1.0, 1.0] - _, removed_pwr_frac, removed_data_frac = spsam.prf_alias_removal( + _, removed_pwr_frac, removed_data_frac = spsa.prf_alias_removal( image.astype("complex64"), xmltree, num_iters, @@ -150,7 +150,7 @@ def test_sicd_alias_mitigation_dilate(example_sicd_alias): threshold=7.0, ) - _, removed_pwr_frac_1, removed_data_frac_1 = spsam.prf_alias_removal( + _, removed_pwr_frac_1, removed_data_frac_1 = spsa.prf_alias_removal( image.astype("complex64"), xmltree, num_iters, @@ -163,7 +163,7 @@ def test_sicd_alias_mitigation_dilate(example_sicd_alias): assert removed_data_frac < removed_data_frac_1 -def test_sicd_alias_mitigation_iters(example_sicd_alias): +def test_sicd_alias_iters(example_sicd_alias): with ( open(example_sicd_alias, "rb") as file, sksicd.NitfReader(file) as reader, @@ -172,7 +172,7 @@ def test_sicd_alias_mitigation_iters(example_sicd_alias): image = reader.read_image() zones = [-1.0, 1.0] - _, removed_pwr_frac, removed_data_frac = spsam.prf_alias_removal( + _, removed_pwr_frac, removed_data_frac = spsa.prf_alias_removal( image.astype("complex64"), xmltree, 2, @@ -180,7 +180,7 @@ def test_sicd_alias_mitigation_iters(example_sicd_alias): threshold=4.0, ) - _, removed_pwr_frac_1, removed_data_frac_1 = spsam.prf_alias_removal( + _, removed_pwr_frac_1, removed_data_frac_1 = spsa.prf_alias_removal( image.astype("complex64"), xmltree, 3, @@ -192,7 +192,7 @@ def test_sicd_alias_mitigation_iters(example_sicd_alias): assert removed_data_frac_1 > removed_data_frac -def test_sicd_alias_mitigation_prf(example_sicd_alias): +def test_sicd_alias_prf(example_sicd_alias): with ( open(example_sicd_alias, "rb") as file, sksicd.NitfReader(file) as reader, @@ -204,7 +204,7 @@ def test_sicd_alias_mitigation_prf(example_sicd_alias): num_iters = 2 zones = [-1.0, 1.0] - _, removed_pwr_frac, removed_data_frac = spsam.prf_alias_removal( + _, removed_pwr_frac, removed_data_frac = spsa.prf_alias_removal( image.astype("complex64"), xmltree, num_iters, @@ -215,15 +215,15 @@ def test_sicd_alias_mitigation_prf(example_sicd_alias): time_poly = sicdew["Grid"]["TimeCOAPoly"] ipp_poly = sicdew["Timeline"]["IPP"]["Set"][0]["IPPPoly"] coa = npp.polyval2d(0, 0, time_poly) - prf = npp.polyval(coa, npp.polyder(ipp_poly)) + prf_override = npp.polyval(coa, npp.polyder(ipp_poly)) del sicdew["Timeline"]["IPP"] - _, removed_pwr_frac_1, removed_data_frac_1 = spsam.prf_alias_removal( + _, removed_pwr_frac_1, removed_data_frac_1 = spsa.prf_alias_removal( image.astype("complex64"), xmltree, num_iters, zones, threshold=4.0, - prf=prf, + prf_override=prf_override, ) assert removed_pwr_frac == removed_pwr_frac_1 From e92cc4a0184da3a79ff41aeabc47afc0962c4bfc Mon Sep 17 00:00:00 2001 From: Brett Johnston Date: Wed, 26 Aug 2026 09:37:43 -0700 Subject: [PATCH 4/4] Apply 1 suggestion(s) to 1 file(s) Co-authored-by: Daniel Pressler --- sarkit_processing/sicd_area_plane.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/sarkit_processing/sicd_area_plane.py b/sarkit_processing/sicd_area_plane.py index 3adbe33..2d49401 100644 --- a/sarkit_processing/sicd_area_plane.py +++ b/sarkit_processing/sicd_area_plane.py @@ -171,8 +171,6 @@ def compute_suitable_rc_plane(sicd_xmltree): eigen_val, eigen_vec = np.linalg.eigh( ground_resolution_proj @ ground_resolution_proj.T ) - eigen_val = eigen_val.real - eigen_vec = eigen_vec.real max_val = np.argmax(eigen_val) scene_coord_angle = np.arctan2(eigen_vec[max_val][1], eigen_vec[max_val][0]) semi_major = np.sqrt(np.abs(eigen_val[max_val]))