Coverage for src/dictk/discontinuity.py: 100%
62 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-09-09 23:57 +0000
1"""Locating a displacement-field discontinuity from correlation-surface peak structure."""
3from typing import NamedTuple
5import numpy as np
6from scipy.signal import find_peaks
8from dictk.correlation import WindowingMethod, _window, zncc
9from dictk.image import PixelCoordinate, SubpixelCoordinate, subimage
12def peak_ratio(
13 *, surface: np.ndarray, height_threshold: float = 0.2, distance: int = 2
14) -> float:
15 """Second-tallest / tallest resolvable peak height along `surface`'s own argmax column.
17 A correlation surface computed from a window straddling a
18 discontinuity shows two comparably-tall peaks, not one -- see
19 [Synthetic Dislocation](../getting_started/synthetic_dislocation.html)
20 and [Experimental
21 Dislocation](../getting_started/experimental_dislocation.html). This
22 turns that qualitative signature into a single number, generalizing
23 `two_peak_separation` (Synthetic Dislocation's own sweep script,
24 which reports the two peaks' separation, not their heights, and
25 exists only to check a *known* offset against a *known* separation).
26 No ground truth is needed here: a value near `1.0` means two
27 comparably-tall peaks -- an ambiguous window, straddling a
28 discontinuity. A value near `0.0` (or exactly `0.0`, one peak only)
29 means a single, unambiguous match -- no discontinuity in this
30 window.
32 Restricted to one column, not the full 2D surface, matching
33 `two_peak_separation`'s own approach: a straddling window's two
34 candidate matches differ in `y` (the two halves' opposite vertical
35 shifts, per
36 [`crack_dislocation`](./image.html#crack_dislocation)) but land at
37 the same `x` (neither half moves horizontally), so both peaks always
38 fall in the same column, the one through the surface's own global
39 maximum.
41 Args:
42 surface: A 2D correlation surface, e.g. from
43 [`dictk.correlation.zncc`](./correlation.html#zncc).
44 height_threshold: Minimum peak height to resolve, as a fraction
45 of the column's own maximum. Passed to
46 `scipy.signal.find_peaks`'s `height`. Default `0.2`, matching
47 `two_peak_separation`'s own threshold.
48 distance: Minimum separation, in pixels, between resolvable
49 peaks. Passed to `scipy.signal.find_peaks`'s `distance`.
50 Default `2`, matching `two_peak_separation`'s own value.
52 Returns:
53 The second-tallest peak's height divided by the tallest's, or
54 `0.0` if fewer than two peaks are resolvable along the column,
55 or if the column's own maximum is not positive (nothing to
56 divide by).
58 Raises:
59 ValueError: If `surface` is not 2D.
60 """
61 if surface.ndim != 2:
62 raise ValueError(f"surface must be 2D, got shape {surface.shape}")
64 x_max = int(np.argmax(surface.max(axis=0)))
65 column = surface[:, x_max]
66 column_max = column.max()
67 if column_max <= 0:
68 return 0.0
70 peaks, props = find_peaks(
71 column, height=height_threshold * column_max, distance=distance
72 )
73 if len(peaks) < 2:
74 return 0.0
76 order = np.argsort(props["peak_heights"])[::-1][:2]
77 tallest, second_tallest = props["peak_heights"][order]
78 return float(second_tallest / tallest)
81class DiscontinuitySweep(NamedTuple):
82 """One [`sweep`](#sweep) result: window-center positions paired with
83 each one's own [`peak_ratio`](#peak_ratio).
85 Attributes:
86 positions: Window-center positions, in `current_image`'s pixel
87 reference frame, in sweep order.
88 peak_ratios: Each position's own `peak_ratio`, index-aligned
89 with `positions`.
90 """
92 positions: list[PixelCoordinate]
93 peak_ratios: list[float]
96def _validate_margins(
97 *,
98 kernel_margin_width: int,
99 kernel_margin_height: int,
100 search_margin_width: int,
101 search_margin_height: int,
102) -> None:
103 """Shared margin validation, matching `translation.locate`'s own checks."""
104 if kernel_margin_width < 1:
105 raise ValueError(f"kernel_margin_width {kernel_margin_width} must be >= 1")
106 if kernel_margin_height < 1:
107 raise ValueError(f"kernel_margin_height {kernel_margin_height} must be >= 1")
108 if search_margin_width <= kernel_margin_width:
109 raise ValueError(
110 f"search_margin_width {search_margin_width} must be greater than "
111 f"kernel_margin_width {kernel_margin_width}"
112 )
113 if search_margin_height <= kernel_margin_height:
114 raise ValueError(
115 f"search_margin_height {search_margin_height} must be greater than "
116 f"kernel_margin_height {kernel_margin_height}"
117 )
120def sweep(
121 *,
122 reference_image: np.ndarray,
123 current_image: np.ndarray,
124 start: PixelCoordinate,
125 end: PixelCoordinate,
126 samples: int,
127 kernel_margin_width: int,
128 kernel_margin_height: int,
129 search_margin_width: int,
130 search_margin_height: int,
131 windowing: WindowingMethod | None = None,
132 height_threshold: float = 0.2,
133 distance: int = 2,
134) -> DiscontinuitySweep:
135 """Evaluate `peak_ratio` at `samples` window-center positions evenly spaced from `start` to `end`.
137 At each position, extracts a kernel from `reference_image` and a
138 search area from `current_image`, both centered on that position
139 (the same "center a window here" step [Synthetic
140 Dislocation](../getting_started/synthetic_dislocation.html#moving-the-window-off-the-crack)
141 does by hand at one fixed row), computes their ZNCC surface, and
142 reads its `peak_ratio`. A discontinuity crossing the `start`-`end`
143 line shows up as a rise in `peak_ratios` near the crossing, high on
144 an ambiguous, straddling window and low everywhere else -- see
145 [`locate`](#locate) to turn that rise into a single position.
147 Args:
148 reference_image: The reference (undeformed) 2D grayscale image.
149 current_image: The current (deformed) 2D grayscale image.
150 start: The first window-center position, in both images' shared
151 pixel reference frame.
152 end: The last window-center position.
153 samples: Number of window-center positions, evenly spaced from
154 `start` to `end` inclusive. Must be >= 2.
155 kernel_margin_width: Half each kernel's width, in pixels. Must be
156 >= 1.
157 kernel_margin_height: Half each kernel's height, in pixels. Must
158 be >= 1.
159 search_margin_width: Half each search area's width, in pixels.
160 Must be greater than `kernel_margin_width`.
161 search_margin_height: Half each search area's height, in pixels.
162 Must be greater than `kernel_margin_height`.
163 windowing: Passed straight through to
164 [`dictk.correlation.window`](./correlation.html#window) for
165 both the kernel and search area at every position. Default
166 `None` applies no windowing.
167 height_threshold: Passed straight through to each position's own
168 [`peak_ratio`](#peak_ratio) call.
169 distance: Passed straight through to each position's own
170 `peak_ratio` call.
172 Returns:
173 A `DiscontinuitySweep` pairing each swept position with its own
174 `peak_ratio`, in sweep order.
176 Raises:
177 ValueError: If `samples` is less than 2, or the margin arguments
178 are invalid (see `kernel_margin_width`/etc. above).
179 """
180 if samples < 2:
181 raise ValueError(f"samples {samples} must be >= 2")
182 _validate_margins(
183 kernel_margin_width=kernel_margin_width,
184 kernel_margin_height=kernel_margin_height,
185 search_margin_width=search_margin_width,
186 search_margin_height=search_margin_height,
187 )
189 positions = []
190 peak_ratios = []
191 for i in range(samples):
192 t = i / (samples - 1)
193 p0 = PixelCoordinate(
194 x=round(start.x + t * (end.x - start.x)),
195 y=round(start.y + t * (end.y - start.y)),
196 )
197 kernel = subimage(
198 image=reference_image,
199 origin=PixelCoordinate(
200 x=p0.x - kernel_margin_width, y=p0.y - kernel_margin_height
201 ),
202 width=2 * kernel_margin_width,
203 height=2 * kernel_margin_height,
204 )
205 search = subimage(
206 image=current_image,
207 origin=PixelCoordinate(
208 x=p0.x - search_margin_width, y=p0.y - search_margin_height
209 ),
210 width=2 * search_margin_width,
211 height=2 * search_margin_height,
212 )
213 kernel, search = _window(kernel=kernel, search=search, windowing=windowing)
214 surface = zncc(kernel=kernel, search=search)
215 positions.append(p0)
216 peak_ratios.append(
217 peak_ratio(
218 surface=surface, height_threshold=height_threshold, distance=distance
219 )
220 )
222 return DiscontinuitySweep(positions=positions, peak_ratios=peak_ratios)
225def _parabolic_vertex(*, y0: float, y1: float, y2: float) -> float:
226 """Fractional offset (in samples, from index 1) of the vertex of the parabola through three equally-spaced points.
228 Pure math, no DIC-specific meaning -- fits $y = a(n - n_0)^2 + c$
229 through `(0, y0)`, `(1, y1)`, `(2, y2)` and returns $n_0 - 1$, the
230 vertex's offset from the center sample. The same 3-point parabolic
231 refinement idea [Parallelism with PyTorch](../getting_started/parallelism_pytorch.html#subpixel-from-a-correlation-surface)
232 already found measurably more accurate than upsampling-based
233 subpixel refinement, applied here along a 1D `peak_ratio` sweep
234 instead of a 2D correlation surface.
236 Args:
237 y0: Value one sample before the center.
238 y1: Value at the center.
239 y2: Value one sample after the center.
241 Returns:
242 The vertex's fractional offset from index 1, in samples.
243 Positive means the true vertex sits past index 1, toward index
244 2. `0.0` if `y0`, `y1`, `y2` are collinear (no curvature to fit
245 -- a division by zero avoided, not a meaningful answer).
246 """
247 denominator = y0 - 2 * y1 + y2
248 if denominator == 0:
249 return 0.0
250 return 0.5 * (y0 - y2) / denominator
253def locate(
254 *,
255 reference_image: np.ndarray,
256 current_image: np.ndarray,
257 start: PixelCoordinate,
258 end: PixelCoordinate,
259 samples: int,
260 kernel_margin_width: int,
261 kernel_margin_height: int,
262 search_margin_width: int,
263 search_margin_height: int,
264 windowing: WindowingMethod | None = None,
265 height_threshold: float = 0.2,
266 distance: int = 2,
267) -> SubpixelCoordinate:
268 """Locate a discontinuity crossing the `start`-`end` line, to subpixel precision.
270 Calls [`sweep`](#sweep) with every argument passed straight through,
271 takes its `peak_ratios`' argmax -- the most ambiguous, most likely
272 straddling window swept -- and refines that sample to a fractional
273 position along the line with a 3-point parabolic fit
274 ([`_parabolic_vertex`](#_parabolic_vertex)) through it and its two
275 neighbors.
277 Args:
278 reference_image: The reference (undeformed) 2D grayscale image.
279 current_image: The current (deformed) 2D grayscale image.
280 start: The first window-center position swept.
281 end: The last window-center position swept.
282 samples: Number of window-center positions swept. Must be >= 2
283 (from `sweep`); in practice needs an argmax with an interior
284 neighbor on each side, so effectively >= 3.
285 kernel_margin_width: Half each kernel's width, in pixels. Must be
286 >= 1.
287 kernel_margin_height: Half each kernel's height, in pixels. Must
288 be >= 1.
289 search_margin_width: Half each search area's width, in pixels.
290 Must be greater than `kernel_margin_width`.
291 search_margin_height: Half each search area's height, in pixels.
292 Must be greater than `kernel_margin_height`.
293 windowing: Passed straight through to `sweep`.
294 height_threshold: Passed straight through to `sweep`.
295 distance: Passed straight through to `sweep`.
297 Returns:
298 The discontinuity's estimated crossing position along the
299 `start`-`end` line, generally fractional.
301 Raises:
302 ValueError: If the `peak_ratios` argmax sits at either sweep
303 endpoint (no interior neighbor to fit a parabola against --
304 the sweep never found a clear internal peak, so refining one
305 would be fabricating precision, not measuring it), or (from
306 `sweep`) if `samples` is less than 2 or the margin arguments
307 are invalid.
308 """
309 result = sweep(
310 reference_image=reference_image,
311 current_image=current_image,
312 start=start,
313 end=end,
314 samples=samples,
315 kernel_margin_width=kernel_margin_width,
316 kernel_margin_height=kernel_margin_height,
317 search_margin_width=search_margin_width,
318 search_margin_height=search_margin_height,
319 windowing=windowing,
320 height_threshold=height_threshold,
321 distance=distance,
322 )
323 ratios = result.peak_ratios
324 i = int(np.argmax(ratios))
325 if i == 0 or i == len(ratios) - 1:
326 raise ValueError(
327 f"peak_ratios argmax sits at sweep endpoint (index {i} of "
328 f"{len(ratios)}); no interior neighbor to refine against"
329 )
331 offset = _parabolic_vertex(y0=ratios[i - 1], y1=ratios[i], y2=ratios[i + 1])
332 fractional_index = i + offset
333 t = fractional_index / (samples - 1)
334 return SubpixelCoordinate(
335 x=start.x + t * (end.x - start.x),
336 y=start.y + t * (end.y - start.y),
337 )