|
9 | 9 | # ---------------------------------------------------------------------------- |
10 | 10 | from .S1_gmtsar import S1_gmtsar |
11 | 11 | from .PRM import PRM |
| 12 | +import numpy as np |
| 13 | + |
| 14 | + |
| 15 | +def _read_slc_patch(src, cy: int, cx: int, half: int) -> np.ndarray: |
| 16 | + """ |
| 17 | + Read a complex SLC patch from an open rasterio dataset. |
| 18 | +
|
| 19 | + Handles both formats: |
| 20 | + - 1 band complex (S1 geotiffs: complex_int16 → complex64) |
| 21 | + - 2 bands real/imag (NISAR: int16 pairs) |
| 22 | +
|
| 23 | + Parameters |
| 24 | + ---------- |
| 25 | + src : rasterio.DatasetReader |
| 26 | + Open rasterio dataset. |
| 27 | + cy, cx : int |
| 28 | + Center coordinates of patch. |
| 29 | + half : int |
| 30 | + Half patch size. |
| 31 | +
|
| 32 | + Returns |
| 33 | + ------- |
| 34 | + np.ndarray |
| 35 | + Complex64 patch of shape (2*half, 2*half). |
| 36 | + """ |
| 37 | + from rasterio.windows import Window |
| 38 | + window = Window(cx - half, cy - half, 2 * half, 2 * half) |
| 39 | + data = src.read(window=window) |
| 40 | + if src.count == 1: |
| 41 | + # Single band complex (S1) |
| 42 | + return data[0].astype(np.complex64) |
| 43 | + else: |
| 44 | + # Two bands: real, imag (NISAR) |
| 45 | + return (data[0] + 1j * data[1]).astype(np.complex64) |
| 46 | + |
| 47 | + |
| 48 | +def _xcorr_refine_slc(ref_path: str, rep_path: str, |
| 49 | + ashift: float, rshift: float, |
| 50 | + stretch_a: float = 0.0, stretch_r: float = 0.0, |
| 51 | + a_stretch_a: float = 0.0, a_stretch_r: float = 0.0, |
| 52 | + patch_size: int = 256, |
| 53 | + min_response: float = 0.1, debug: bool = False) -> dict: |
| 54 | + """ |
| 55 | + Run xcorr refinement by reading patches directly from geotiff files. |
| 56 | +
|
| 57 | + Reads only the required patches from disk, avoiding full image load. |
| 58 | + Handles both S1 (1 band complex) and NISAR (2 bands real/imag) formats. |
| 59 | +
|
| 60 | + Parameters |
| 61 | + ---------- |
| 62 | + ref_path : str |
| 63 | + Path to reference SLC geotiff. |
| 64 | + rep_path : str |
| 65 | + Path to repeat SLC geotiff. |
| 66 | + ashift, rshift : float |
| 67 | + Geometry-based azimuth and range shifts. |
| 68 | + stretch_a, stretch_r, a_stretch_a, a_stretch_r : float |
| 69 | + Geometry-based stretch parameters. |
| 70 | + patch_size : int |
| 71 | + Xcorr patch size. Default 256. |
| 72 | + min_response : float |
| 73 | + Minimum correlation response. |
| 74 | + debug : bool |
| 75 | + Print debug info. |
| 76 | +
|
| 77 | + Returns |
| 78 | + ------- |
| 79 | + dict |
| 80 | + Correction coefficients (same as xcorr_fitoffset output). |
| 81 | +
|
| 82 | + Raises |
| 83 | + ------ |
| 84 | + RuntimeError |
| 85 | + If xcorr failed (insufficient valid patches). |
| 86 | + """ |
| 87 | + import rasterio |
| 88 | + from .utils_satellite import xcorr_patch, xcorr_fitoffset |
| 89 | + |
| 90 | + half = patch_size // 2 |
| 91 | + hann = np.outer(np.hanning(patch_size), np.hanning(patch_size)).astype(np.float32) |
| 92 | + results = [] |
| 93 | + |
| 94 | + with rasterio.open(ref_path) as src_ref, rasterio.open(rep_path) as src_rep: |
| 95 | + ny_ref, nx_ref = src_ref.height, src_ref.width |
| 96 | + ny_rep, nx_rep = src_rep.height, src_rep.width |
| 97 | + |
| 98 | + # Auto-compute grid: ~2x patch spacing, minimum 4 patches per dimension |
| 99 | + n_rows = max(4, (ny_ref - patch_size) // (2 * patch_size) + 1) |
| 100 | + n_cols = max(4, (nx_ref - patch_size) // (2 * patch_size) + 1) |
| 101 | + grid = (n_rows, n_cols) |
| 102 | + |
| 103 | + if debug: |
| 104 | + print(f"Xcorr refinement: {grid[0]}×{grid[1]} = {grid[0]*grid[1]} patches, size {patch_size}") |
| 105 | + print(f"Image sizes: ref={ny_ref}×{nx_ref}, rep={ny_rep}×{nx_rep}") |
| 106 | + print(f"Geometry params: ashift={ashift:.2f}, rshift={rshift:.2f}") |
| 107 | + |
| 108 | + n_rows, n_cols = grid |
| 109 | + for row in range(n_rows): |
| 110 | + cy1 = int((row + 0.5) * ny_ref / n_rows) |
| 111 | + for col in range(n_cols): |
| 112 | + cx1 = int((col + 0.5) * nx_ref / n_cols) |
| 113 | + |
| 114 | + # Apply geometry offset - compute float position first |
| 115 | + cy2_float = cy1 + ashift + stretch_a * cx1 + a_stretch_a * cy1 |
| 116 | + cx2_float = cx1 + rshift + stretch_r * cx1 + a_stretch_r * cy1 |
| 117 | + |
| 118 | + # Truncate to integer for patch reading |
| 119 | + cy2 = int(cy2_float) |
| 120 | + cx2 = int(cx2_float) |
| 121 | + |
| 122 | + # Track truncation artifact - phaseCorrelate will "find" this sub-pixel |
| 123 | + # and we need to subtract it to get the TRUE residual |
| 124 | + frac_a = cy2_float - cy2 |
| 125 | + frac_r = cx2_float - cx2 |
| 126 | + |
| 127 | + # Bounds check |
| 128 | + if cy1 < half or cy1 > ny_ref - half: |
| 129 | + continue |
| 130 | + if cy2 < half or cy2 > ny_rep - half: |
| 131 | + continue |
| 132 | + if cx1 < half or cx1 > nx_ref - half: |
| 133 | + continue |
| 134 | + if cx2 < half or cx2 > nx_rep - half: |
| 135 | + continue |
| 136 | + |
| 137 | + # Read and correlate patches |
| 138 | + patch1 = _read_slc_patch(src_ref, cy1, cx1, half) |
| 139 | + patch2 = _read_slc_patch(src_rep, cy2, cx2, half) |
| 140 | + |
| 141 | + result = xcorr_patch(patch1, patch2, hann, min_response=min_response) |
| 142 | + if result is not None: |
| 143 | + # Compensate for int() truncation artifact |
| 144 | + # phaseCorrelate "finds" the sub-pixel that was lost by truncation |
| 145 | + # Subtract it to get the TRUE residual beyond the geometry |
| 146 | + result['dy'] -= frac_a |
| 147 | + result['dx'] -= frac_r |
| 148 | + result['cy1'] = cy1 |
| 149 | + result['cx1'] = cx1 |
| 150 | + results.append(result) |
| 151 | + |
| 152 | + if debug: |
| 153 | + print(f"Xcorr results: {len(results)} with response > {min_response}") |
| 154 | + |
| 155 | + # Fit bilinear using full radar extent for normalization |
| 156 | + corrections = xcorr_fitoffset(results, nx=nx_ref, ny=ny_ref, debug=debug) |
| 157 | + |
| 158 | + if corrections is None: |
| 159 | + raise RuntimeError("Xcorr fitoffset failed - insufficient valid patches. Scene alignment cannot continue.") |
| 160 | + |
| 161 | + if debug: |
| 162 | + print(f"Xcorr fitoffset result:") |
| 163 | + print(f" ashift={corrections['ashift']:.4f}, stretch_a={corrections['stretch_a']:.8f}, a_stretch_a={corrections['a_stretch_a']:.8f}") |
| 164 | + print(f" rshift={corrections['rshift']:.4f}, stretch_r={corrections['stretch_r']:.8f}, a_stretch_r={corrections['a_stretch_r']:.8f}") |
| 165 | + |
| 166 | + return corrections |
12 | 167 |
|
13 | 168 |
|
14 | 169 | class S1_align(S1_gmtsar): |
@@ -250,7 +405,6 @@ def align_rep(self, burst_rep: str, burst_ref: str, prm_ref: "PRM", |
250 | 405 |
|
251 | 406 | # Xcorr refinement: measure actual offsets and correct geometry alignment |
252 | 407 | if xcorr is not None: |
253 | | - from .utils_satellite import xcorr_refine_slc |
254 | 408 | import os |
255 | 409 | import math |
256 | 410 |
|
@@ -305,7 +459,7 @@ def align_rep(self, burst_rep: str, burst_ref: str, prm_ref: "PRM", |
305 | 459 | print(f" Geometry (TIFF): ashift={geom_ashift_tiff:.2f}") |
306 | 460 |
|
307 | 461 | # Run xcorr WITH geometry pre-applied (in TIFF space) to find small residuals |
308 | | - xcorr_params = xcorr_refine_slc( |
| 462 | + xcorr_params = _xcorr_refine_slc( |
309 | 463 | ref_tiff, rep_tiff, |
310 | 464 | ashift=geom_ashift_tiff, rshift=geom_rshift, |
311 | 465 | stretch_a=geom_stretch_a, stretch_r=geom_stretch_r, |
|
0 commit comments