Skip to content

Revert "Compressed slc forward" #749

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mirzaees merged 1 commit into from
Oct 27, 2026
Merged

Revert "Compressed slc forward" #749

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 8 additions & 15 deletions src/dolphin/phase_link/_compress.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
def compress(
slc_stack: ArrayLike,
pl_cpx_phase: ArrayLike,
is_real_mask: ArrayLike | None = None,
first_real_slc_idx: int = 0,
slc_mean: ArrayLike | None = None,
reference_idx: int | None = None,
):
Expand All @@ -24,12 +24,10 @@ def compress(
pl_cpx_phase : ArrayLike
The estimated complex phase from phase linking.
shape = (nslc, rows // strides.y, cols // strides.x)
is_real_mask : ArrayLike, optional
Boolean mask of shape (nslc,), where `True` marks a "real" (non-compressed)
SLC within `slc_stack`. Entries marked `False` (already-Compressed SLCs)
are excluded from the dot product during compression, regardless of their
position in the stack.
If None, all entries are treated as real.
first_real_slc_idx : int
Index within `slc_stack` where the "real" SLCs start.
Indexes before this (which are assumed to be already Compressed SLCs) are
excluded from the dot product during compression
slc_mean : ArrayLike, optional
The mean SLC magnitude, shape (rows, cols), to use as output pixel magnitudes.
If None, the mean is computed from the input SLC stack.
Expand All @@ -49,14 +47,9 @@ def compress(
else:
pl_referenced = pl_cpx_phase

if is_real_mask is None:
is_real_mask = np.ones(np.asarray(slc_stack).shape[0], dtype=bool)
else:
is_real_mask = np.asarray(is_real_mask, dtype=bool)

# Exclude the already-compressed layers *after* doing the reference.
pl_referenced = pl_referenced[is_real_mask, :, :]
slcs = slc_stack[is_real_mask]
# Slice away the compressed layers *after* doing the reference.
pl_referenced = pl_referenced[first_real_slc_idx:, :, :]
slcs = slc_stack[first_real_slc_idx:]
# If the output is downsampled, we need to make `pl_cpx_phase` the same shape
# as the output
pl_estimate_upsampled = upsample_nearest(pl_referenced, slcs.shape[1:])
Expand Down
45 changes: 19 additions & 26 deletions src/dolphin/phase_link/_core.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@ def run_phase_linking(
avg_mag: ArrayLike | None = None,
use_slc_amp: bool = False,
baseline_lag: Optional[int] = None,
last_compressed_slc_idx: Optional[int] = None,
first_real_slc_idx: int = 0,
compute_crlb: bool = True,
) -> PhaseLinkOutput:
"""Estimate the linked phase for a stack of SLCs.
Expand Down Expand Up @@ -144,12 +144,10 @@ def run_phase_linking(
or to set the SLC amplitude to 1.0. By default False.
baseline_lag : int, optional, default=None
lag for temporal baseline to do short temporal baseline inversion (STBAS)
last_compressed_slc_idx : int, optional, default = None
The index of the most recent (last) compressed SLC in the stack.
This is only used for the CRLB computation, to pick which date's phase
to treat as the fixed/known reference (0 variance) when inverting the
Fisher Information Matrix. If None (no compressed SLC in the stack),
the first date (index 0) is used.
first_real_slc_idx : int, optional, default = 0
The index of the first real SLC in the stack.
This is only used for the CRLB computation.
By default 0.
compute_crlb : bool, optional
Whether to compute the CRLB, by default True

Expand Down Expand Up @@ -207,7 +205,7 @@ def run_phase_linking(
reference_idx=reference_idx,
neighbor_arrays=neighbor_arrays,
baseline_lag=baseline_lag,
last_compressed_slc_idx=last_compressed_slc_idx,
first_real_slc_idx=first_real_slc_idx,
compute_crlb=compute_crlb,
)

Expand Down Expand Up @@ -266,7 +264,7 @@ def run_cpl(
reference_idx: int = 0,
neighbor_arrays: Optional[np.ndarray] = None,
baseline_lag: Optional[int] = None,
last_compressed_slc_idx: Optional[int] = None,
first_real_slc_idx: int = 0,
compute_crlb: bool = True,
) -> PhaseLinkOutput:
"""Run the Combined Phase Linking (CPL) algorithm.
Expand Down Expand Up @@ -304,12 +302,10 @@ def run_cpl(
StBAS parameter to include only nearest-N interferograms for phase linking.
A `baseline_lag` of `n` will only include the closest `n` interferograms.
`baseline_line` must be positive.
last_compressed_slc_idx : int, optional, default = None
The index of the most recent (last) compressed SLC in the stack.
This is only used for the CRLB computation, to pick which date's phase
to treat as the fixed/known reference (0 variance) when inverting the
Fisher Information Matrix. If None (no compressed SLC in the stack),
the first date (index 0) is used.
first_real_slc_idx : int, optional, default = 0
The index of the first real SLC in the stack.
This is only used for the CRLB computation.
By default 0.
compute_crlb : bool, optional
Whether to compute the CRLB, by default True

Expand Down Expand Up @@ -358,7 +354,7 @@ def run_cpl(
zero_correlation_threshold=zero_correlation_threshold,
reference_idx=reference_idx,
num_looks=num_looks,
last_compressed_slc_idx=last_compressed_slc_idx,
first_real_slc_idx=first_real_slc_idx,
compute_crlb=compute_crlb,
)
# Get the temporal coherence
Expand Down Expand Up @@ -398,7 +394,7 @@ def run_cpl(
"beta",
"reference_idx",
"num_looks",
"last_compressed_slc_idx",
"first_real_slc_idx",
"compute_crlb",
),
)
Expand All @@ -409,7 +405,7 @@ def process_coherence_matrices(
zero_correlation_threshold: float = 0.0,
reference_idx: int = 0,
num_looks: int = 1,
last_compressed_slc_idx: Optional[int] = None,
first_real_slc_idx: int = 0,
compute_crlb: bool = True,
) -> tuple[Array, Array, Array, Array]:
"""Estimate the linked phase for a stack of coherence matrices.
Expand Down Expand Up @@ -438,12 +434,10 @@ def process_coherence_matrices(
num_looks : int, optional
The number of looks used to form the input correlation data, used
during CRLB computation.
last_compressed_slc_idx : int, optional, default = None
The index of the most recent (last) compressed SLC in the stack.
This is only used for the CRLB computation, to pick which date's phase
to treat as the fixed/known reference (0 variance) when inverting the
Fisher Information Matrix. If None (no compressed SLC in the stack),
the first date (index 0) is used.
first_real_slc_idx : int, optional, default = 0
The index of the first real SLC in the stack.
This is only used for the CRLB computation.
By default 0.
compute_crlb : bool, optional
Whether to compute the CRLB
Default is True.
Expand Down Expand Up @@ -535,8 +529,7 @@ def process_coherence_matrices(
if compute_crlb:
# Build X once and do the inverse-free CRLB from X
X = crlb._build_fisher_from_abs_gamma(Gamma, Gamma_inv, num_looks)
crlb_ref_idx = 0 if last_compressed_slc_idx is None else last_compressed_slc_idx
crlb_std_dev = crlb._crlb_from_x(X, crlb_ref_idx, 0, 1e-6)
crlb_std_dev = crlb._crlb_from_x(X, max(first_real_slc_idx - 1, 0), 0, 1e-6)

else:
crlb_std_dev = jnp.zeros(C_arrays.shape[:-1], dtype=jnp.float32)
Expand Down
100 changes: 25 additions & 75 deletions src/dolphin/stack.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,29 +136,10 @@ def first_real_slc_idx(self) -> int:
msg = "No real SLCs in ministack"
raise ValueError(msg) from e

@property
def is_real_mask(self) -> np.ndarray:
"""Boolean mask marking the real (non-compressed) SLCs.

Unlike `first_real_slc_idx`, this is correct even when compressed SLCs
are not a contiguous prefix of `file_list`.
"""
return ~np.array(self.is_compressed)

@property
def last_compressed_slc_idx(self) -> Optional[int]:
"""Index of the most recent (last) compressed SLC.

or None if there isn't one.
"""
idxs = np.where(np.array(self.is_compressed))[0]
return int(idxs[-1]) if idxs.size else None

@property
def real_slc_date_range(self) -> tuple[DateOrDatetime, DateOrDatetime]:
"""Date range of the real SLCs in the ministack."""
real_idxs = np.where(self.is_real_mask)[0]
return (self.dates[real_idxs[0]][0], self.dates[real_idxs[-1]][-1])
return (self.dates[self.first_real_slc_idx][0], self.dates[-1][-1])

@property
def real_slc_date_range_str(self) -> str:
Expand Down Expand Up @@ -506,56 +487,31 @@ def plan(
# TODO: will we ever actually need to read the old metadata here?
compressed_slc_infos.append(CompressedSlcInfo.from_filename(f))

# Solve each ministack using current chunk (and the previous compressed SLCs).
# Chunk over the *real*-SLC-flagged positions directly (gathering by index,
# not a contiguous `range()`/`slice()`) since a compressed SLC can be
# interleaved among real SLCs elsewhere in `self.file_list` (e.g.
# forward-mode/incremental reprocessing inputs sorted by
# `opera_utils.sort_files_by_date`, which has no compressed-first
# guarantee). A raw positional slice could otherwise pick up an
# interleaved compressed SLC a second time, or come up short on real SLCs.
real_idxs = np.where(~np.array(self.is_compressed))[0]
real_idx_chunks = [
real_idxs[i : i + ministack_size]
for i in range(0, len(real_idxs), ministack_size)
]
# Solve each ministack using current chunk (and the previous compressed SLCs)
ministack_starts = range(
self.first_real_slc_idx, len(self.file_list), ministack_size
)

for cur_real_idxs in real_idx_chunks:
cur_files = [self.file_list[i] for i in cur_real_idxs]
cur_dates = [self.dates[i] for i in cur_real_idxs]
cur_is_compressed = [False] * len(cur_real_idxs)
for full_stack_idx in ministack_starts:
cur_slice = slice(full_stack_idx, full_stack_idx + ministack_size)
cur_files = list(self.file_list[cur_slice]).copy()
cur_dates = list(self.dates[cur_slice]).copy()

# Read compressed*.tif files and if they do not exist use the compressed*.h5
# Limit the num comp slcs to `max_num_compressed`
cur_comp_slc_infos = compressed_slc_infos[-self.max_num_compressed :]
# num_ccslc = len(cur_comp_slc_infos)

# Chronologically merge the existing compressed SLCs (keyed by their
# `reference_date`) with this chunk's real SLCs, instead of blindly
# putting all compressed SLCs first. A compressed SLC's reference date
# can be later than some real SLC dates also present here, so a plain
# prefix-concat would silently mislabel which entries are "real".
# `sort()` is stable, so ties keep compressed-before-real (today's
# behavior), making this a no-op reordering in the fully-contiguous case.
entries = [
(c.dates[0], c.path, c.dates, True) for c in cur_comp_slc_infos
] + [
(d[0], f, d, is_comp)
for f, d, is_comp in zip(
cur_files, cur_dates, cur_is_compressed, strict=False
)
]
entries.sort(key=lambda e: e[0])

combined_files = [e[1] for e in entries]
combined_dates = [e[2] for e in entries]
combined_is_compressed = [e[3] for e in entries]

# Position of the most recent compressed SLC in the merged list.
# Equals `num_ccslc - 1` in the fully-contiguous case (preserving
# today's behavior exactly); may differ once interleaved.
_comp_idxs = [i for i, c in enumerate(combined_is_compressed) if c]
last_ccslc_idx = _comp_idxs[-1] if _comp_idxs else None
comp_slc_files = [c.path for c in compressed_slc_infos]
# Add the existing compressed SLC files to the start, but
# limit the num comp slcs `max_num_compressed`
cur_comp_slc_files = comp_slc_files[-self.max_num_compressed :]
combined_files = cur_comp_slc_files + cur_files

combined_dates = [
c.dates for c in compressed_slc_infos[-self.max_num_compressed :]
] + cur_dates

num_ccslc = len(cur_comp_slc_files)
combined_is_compressed = num_ccslc * [True] + list(
self.is_compressed[cur_slice]
)

# Make the current ministack output folder using the start/end dates
new_date_str = format_dates(
Expand All @@ -567,13 +523,9 @@ def plan(
compressed_reference_idx = compressed_idx
elif self.compressed_slc_plan == CompressedSlcPlan.ALWAYS_FIRST:
# Here, CompSLCs have same base phase, but different "residual" added on
compressed_reference_idx = (
last_ccslc_idx if last_ccslc_idx is not None else 0
)
compressed_reference_idx = max(0, num_ccslc - 1)
elif self.compressed_slc_plan == CompressedSlcPlan.LAST_PER_MINISTACK:
# Here, CompSLCs have same base phase, but different "residual" added on
# `-1` means "chronologically last entry", which stays correct
# after sorting by date.
compressed_reference_idx = -1

# Set the `output_reference_idx`, used for making interferograms
Expand All @@ -586,9 +538,7 @@ def plan(
# looking like single-reference, relative to day 1
# For `LAST_PER_MINISTACK`, the interferograms are formed
# which are the shortest possible temporal baseline for the given inputs
output_reference_idx = (
last_ccslc_idx if last_ccslc_idx is not None else 0
)
output_reference_idx = max(0, num_ccslc - 1)

cur_ministack = MiniStackInfo(
file_list=combined_files,
Expand Down
46 changes: 18 additions & 28 deletions src/dolphin/workflows/single.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,16 +96,13 @@ def run_wrapped_phase_single(
xhalf, yhalf = half_window["x"], half_window["y"]

# If we were passed any compressed SLCs in `input_slc_files`,
# then we want a mask for when we create new compressed SLCs.
# We skip the old compressed SLCs to create new ones.
# (A boolean mask, rather than a single boundary index, since a compressed
# SLC is not guaranteed to be a contiguous prefix of `input_slc_files`.)
is_real_mask = ministack.is_real_mask
num_real = int(is_real_mask.sum())
# then we want that index for when we create new compressed SLCs.
# We skip the old compressed SLCs to create new ones
first_real_slc_idx = ministack.first_real_slc_idx

msg = (
f"Processing {num_real} SLCs +"
f" {len(input_slc_files) - num_real} compressed SLCs. "
f"Processing {len(input_slc_files) - first_real_slc_idx} SLCs +"
f" {first_real_slc_idx} compressed SLCs. "
)
logger.info(msg)

Expand Down Expand Up @@ -262,7 +259,7 @@ def _process_block(
neighbor_arrays=neighbor_arrays,
baseline_lag=baseline_lag,
avg_mag=amp_mean[in_rows, in_cols] if amp_mean is not None else None,
last_compressed_slc_idx=ministack.last_compressed_slc_idx,
first_real_slc_idx=ministack.first_real_slc_idx,
compute_crlb=write_crlb,
)
except PhaseLinkRuntimeError as e:
Expand All @@ -285,34 +282,36 @@ def _process_block(

# Compress the ministack using only the non-compressed SLCs
# Get the mean to set as pixel magnitudes
abs_stack = np.abs(cur_data[is_real_mask][:, in_trim_rows, in_trim_cols])
abs_stack = np.abs(cur_data[first_real_slc_idx:, in_trim_rows, in_trim_cols])
cur_data_mean, cur_amp_dispersion, _ = calc_ps_block(abs_stack)
cur_comp_slc = compress(
# Get the inner portion of the full-res SLC data
cur_data[:, in_trim_rows, in_trim_cols],
pl_output.cpx_phase[:, out_trim_rows, out_trim_cols],
is_real_mask=is_real_mask,
first_real_slc_idx=first_real_slc_idx,
slc_mean=cur_data_mean,
reference_idx=ministack.compressed_reference_idx,
)

# Save each of the MLE estimates (ignoring those corresponding to
# compressed SLCs indexes)
assert num_real == len(phase_linked_slc_files)
assert len(pl_output.cpx_phase[first_real_slc_idx:]) == len(
phase_linked_slc_files
)
# ### Save results ###
with write_lock:
# ### Save results ###
for img, f in zip(
pl_output.cpx_phase[is_real_mask][:, out_trim_rows, out_trim_cols],
pl_output.cpx_phase[first_real_slc_idx:, out_trim_rows, out_trim_cols],
phase_linked_slc_files,
strict=True,
):
writer.queue_write(img, f, out_rows.start, out_cols.start)

if write_crlb:
for img, f in zip(
pl_output.crlb_std_dev[is_real_mask][
:, out_trim_rows, out_trim_cols
pl_output.crlb_std_dev[
first_real_slc_idx:, out_trim_rows, out_trim_cols
],
phase_linked_crlb_files,
strict=True,
Expand Down Expand Up @@ -492,26 +491,17 @@ def _get_amp_mean_variance(
return None, None


def _real_date_str_list(ministack: MiniStackInfo) -> list[str]:
"""Get the date strings for only the real (non-compressed) SLCs."""
return [
d
for d, is_real in zip(
ministack.get_date_str_list(), ministack.is_real_mask, strict=True
)
if is_real
]


def _name_slcs(ministack: MiniStackInfo) -> list[str]:
"""Generate SLC filenames for the ministack."""
date_strs = _real_date_str_list(ministack)
start_idx = ministack.first_real_slc_idx
date_strs = ministack.get_date_str_list()[start_idx:]
return [f"{Path(d).stem}.slc.tif" for d in date_strs]


def _name_crlbs(ministack: MiniStackInfo) -> list[str]:
"""Generate CRLB filenames for the ministack."""
date_strs = _real_date_str_list(ministack)
start_idx = ministack.first_real_slc_idx
date_strs = ministack.get_date_str_list()[start_idx:]
return [f"crlb_{Path(d).stem}.tif" for d in date_strs]


Expand Down
Loading
Loading