Coverage for src/dictk/translation.py: 100%
48 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-08-25 21:06 +0000
1"""Point translation tracking between a reference and current image."""
3import numpy as np
4from skimage.registration import phase_cross_correlation
6from dictk.correlation import WindowingMethod, _kernel_pad, _window
7from dictk.image import PixelCoordinate, SubpixelCoordinate, subimage
10def locate(
11 *,
12 reference_image: np.ndarray,
13 current_image: np.ndarray,
14 reference_point: PixelCoordinate,
15 search_center: PixelCoordinate,
16 kernel_margin_width: int,
17 kernel_margin_height: int,
18 search_margin_width: int,
19 search_margin_height: int,
20 windowing: WindowingMethod | None = None,
21) -> PixelCoordinate:
22 """Given a `reference_point` expressed in the `reference_image`
23 frame, find its position expressed in the `current_image` frame.
25 `reference_point` is a fixed, already-known point in `reference_image`.
26 Where that same physical point ends up in `current_image` is exactly
27 what this function finds — it is *not* an input, it's the return
28 value. `search_center` is a different thing: just where to center the
29 search area in `current_image`, i.e., a guess of roughly where to
30 look. It doesn't need to be exact — only close enough that the true
31 (unknown) position falls within the search area — and if no better
32 guess is available, passing `reference_point` again is a reasonable
33 default.
35 Concretely: extracts a rectangular kernel (also called a subset,
36 filter, or convolution matrix; `2 * kernel_margin_width` wide,
37 `2 * kernel_margin_height` tall) from `reference_image` centered at
38 `reference_point`, and a rectangular search area (also called a
39 search window, scanning zone, or area of interest (AOI);
40 `2 * search_margin_width` wide, `2 * search_margin_height` tall)
41 from `current_image` centered at `search_center`, then locates the
42 kernel within the search area via FFT-based phase cross-correlation
43 (`skimage.registration.phase_cross_correlation`). The kernel always
44 comes from `reference_image`; the search area always comes from
45 `current_image`. Integer-pixel precision only; subpixel refinement is
46 out of scope for now.
48 The true position can be anywhere within `search_margin_width`/
49 `search_margin_height` of `search_center` in every direction and
50 still be found correctly -- the recoverable range is symmetric,
51 bounded only by the search margins themselves, not by
52 `kernel_margin_width`/`kernel_margin_height`. See [Recoverable
53 Displacement Range](../getting_started/recoverable_displacement_range.html)
54 for why that's worth stating explicitly: an earlier version of this
55 function had a real, silent bug here -- it recovered a displacement
56 in the negative direction up to the full search margin, but capped
57 at exactly the kernel margin in the positive direction, past which it
58 returned a confidently wrong position instead of failing visibly.
60 Args:
61 reference_image: The reference (undeformed) 2D grayscale image.
62 current_image: The current (deformed) 2D grayscale image.
63 reference_point: The point's fixed, known position, in
64 `reference_image`'s pixel reference frame.
65 search_center: Where to center the search area, in
66 `current_image`'s pixel reference frame — a guess of roughly
67 where `reference_point` ended up, not the answer itself.
68 kernel_margin_width: Half the kernel's width, in pixels. Must
69 be >= 1.
70 kernel_margin_height: Half the kernel's height, in pixels. Must
71 be >= 1.
72 search_margin_width: Half the search area's width, in pixels.
73 Must be greater than `kernel_margin_width`.
74 search_margin_height: Half the search area's height, in pixels.
75 Must be greater than `kernel_margin_height`.
76 windowing: If given, both the kernel and search area are passed
77 through [`dictk.correlation.window`](./correlation.html#window)
78 with this method -- tapering their edges toward zero to
79 reduce spectral leakage -- before padding/FFT, the same
80 `windowing` parameter
81 [`dictk.correlation.phase_correlation`](./correlation.html#phase_correlation)
82 exposes for the surface this function doesn't return. Default
83 `None` applies no windowing, matching this function's
84 original behavior exactly.
86 Returns:
87 The point's location, in `current_image`'s pixel reference frame.
89 Raises:
90 ValueError: If `kernel_margin_width` or `kernel_margin_height` is
91 less than 1, or either `search_margin_width`/
92 `search_margin_height` is not greater than its kernel
93 counterpart.
94 """
95 if kernel_margin_width < 1:
96 raise ValueError(f"kernel_margin_width {kernel_margin_width} must be >= 1")
97 if kernel_margin_height < 1:
98 raise ValueError(f"kernel_margin_height {kernel_margin_height} must be >= 1")
99 if search_margin_width <= kernel_margin_width:
100 raise ValueError(
101 f"search_margin_width {search_margin_width} must be greater than "
102 f"kernel_margin_width {kernel_margin_width}"
103 )
104 if search_margin_height <= kernel_margin_height:
105 raise ValueError(
106 f"search_margin_height {search_margin_height} must be greater than "
107 f"kernel_margin_height {kernel_margin_height}"
108 )
110 kernel_width = 2 * kernel_margin_width
111 kernel_height = 2 * kernel_margin_height
112 search_width = 2 * search_margin_width
113 search_height = 2 * search_margin_height
115 kernel_origin = PixelCoordinate(
116 x=reference_point.x - kernel_margin_width,
117 y=reference_point.y - kernel_margin_height,
118 )
119 kernel = subimage(
120 image=reference_image,
121 origin=kernel_origin,
122 width=kernel_width,
123 height=kernel_height,
124 )
126 search_origin = PixelCoordinate(
127 x=search_center.x - search_margin_width,
128 y=search_center.y - search_margin_height,
129 )
130 search = subimage(
131 image=current_image,
132 origin=search_origin,
133 width=search_width,
134 height=search_height,
135 )
137 # phase_cross_correlation requires both images the same shape; only
138 # the kernel needs padding, since the search area is always larger.
139 # _window optionally tapers both first -- the same shared step
140 # `phase_correlation` uses -- then _kernel_pad grows kernel up to
141 # search's own shape.
142 #
143 # centered=True matters: phase_cross_correlation's FFT-based shift is
144 # only correct up to half the padded array's own size in each
145 # direction (it's circular/periodic) -- past that it wraps around to
146 # a confidently wrong answer instead of failing visibly. With the
147 # default centered=False (kernel's content anchored at the padded
148 # array's top-left corner), that safe range is asymmetric: unbounded
149 # in the negative direction, but capped at exactly kernel_margin_width/
150 # kernel_margin_height in the positive direction, no matter how large
151 # search_margin is set. Centering kernel's content in the padded array
152 # instead makes the safe range symmetric in both directions -- see
153 # [Recoverable Displacement
154 # Range](../getting_started/recoverable_displacement_range.html) for the
155 # derivation and how this was found.
156 kernel, search = _window(kernel=kernel, search=search, windowing=windowing)
157 kernel_padded, pad_before_height, pad_before_width = _kernel_pad(
158 kernel=kernel, shape=search.shape, centered=True
159 )
161 # search=reference_image, kernel_padded=moving_image, not the other
162 # way around: skimage docs define `shift` as "the shift required to
163 # register moving_image with reference_image", and what we actually
164 # want is where the kernel's content sits within the (larger) search
165 # area, so search must be reference_image and kernel_padded must be
166 # moving_image. Swapping them negates the result and returns the
167 # wrong location -- verified empirically, not just by this reasoning.
168 #
169 # normalization="phase" (skimage's own default) is used here
170 # deliberately -- we do NOT follow hdic's registration.py, which uses
171 # normalization=None. hdic's None came from its own pipeline running
172 # a separate ZNCC-style preprocess() step before calling
173 # phase_cross_correlation, specifically to avoid double-normalizing
174 # on top of that. locate() has no such preprocessing step, so that
175 # reasoning doesn't carry over here -- and empirically, sweeping
176 # kernel/search margin ratios against a known displacement,
177 # normalization="phase" matched or outperformed normalization=None,
178 # not the other way around.
179 #
180 # What normalization actually does: the core computation is a
181 # cross-power spectrum, image_product = FFT(reference) *
182 # conj(FFT(moving)). normalization="phase" divides this by its own
183 # magnitude at every frequency, discarding contrast/energy
184 # information and keeping only phase -- the classic Kuglin-Hines
185 # phase correlation, giving a sharp peak and robustness to
186 # illumination differences between the two images.
187 # normalization=None skips that division (plain unnormalized
188 # cross-correlation), which skimage's own docs describe as less
189 # robust to noise but sometimes preferable in high-noise scenarios.
190 # Which is better is genuinely content-dependent, not a settled
191 # default-is-always-better situation.
192 shift, _error, _phasediff = phase_cross_correlation(
193 reference_image=search, moving_image=kernel_padded, normalization="phase"
194 )
196 return PixelCoordinate(
197 x=search_origin.x + int(shift[1]) + kernel_margin_width + pad_before_width,
198 y=search_origin.y + int(shift[0]) + kernel_margin_height + pad_before_height,
199 )
202def locate_subpixel(
203 *,
204 reference_image: np.ndarray,
205 current_image: np.ndarray,
206 reference_point: PixelCoordinate,
207 search_center: PixelCoordinate,
208 kernel_margin_width: int,
209 kernel_margin_height: int,
210 search_margin_width: int,
211 search_margin_height: int,
212 windowing: WindowingMethod | None = None,
213 upsample_factor: int = 100,
214) -> SubpixelCoordinate:
215 """The subpixel-accurate sibling of [`locate`](#locate): recovers a
216 fractional position instead of rounding it away.
218 Same kernel/search extraction, windowing, and centered padding as
219 `locate` -- deliberately a separate function, not a parameter added
220 to it, so `locate`'s own return type (`PixelCoordinate`, always
221 `int`) never changes shape based on an argument. `locate` truncates
222 `phase_cross_correlation`'s shift to the nearest whole pixel with
223 `int()`; this function instead requests a refined, generally
224 fractional shift from `phase_cross_correlation` itself
225 (`upsample_factor`) and returns it directly, undiscarded.
227 `upsample_factor` does not "fix" `locate`'s own exact-integer
228 matching in the way it might sound like it should: if the true
229 displacement genuinely isn't an integer (a stretch's own resampling
230 can easily make this true even for an otherwise-integer-pixel
231 target -- see [Subpixel
232 Accuracy](../getting_started/subpixel_accuracy.html)), no amount of
233 refinement makes `locate`'s truncated answer correct, because
234 `locate`'s question ("which integer pixel?") isn't the right one to
235 ask anymore. What subpixel refinement *does* reliably improve is how
236 close the returned position lands to the true, generally fractional,
237 target -- a different, and for a real (non-integer) displacement, a
238 more honest question.
240 Args:
241 reference_image: The reference (undeformed) 2D grayscale image.
242 current_image: The current (deformed) 2D grayscale image.
243 reference_point: The point's fixed, known position, in
244 `reference_image`'s pixel reference frame.
245 search_center: Where to center the search area, in
246 `current_image`'s pixel reference frame -- see `locate`'s
247 own docstring for the full explanation of this parameter.
248 kernel_margin_width: Half the kernel's width, in pixels. Must
249 be >= 1.
250 kernel_margin_height: Half the kernel's height, in pixels. Must
251 be >= 1.
252 search_margin_width: Half the search area's width, in pixels.
253 Must be greater than `kernel_margin_width`.
254 search_margin_height: Half the search area's height, in pixels.
255 Must be greater than `kernel_margin_height`.
256 windowing: Passed straight through to
257 [`dictk.correlation.window`](./correlation.html#window),
258 same as `locate`'s own `windowing` parameter.
259 upsample_factor: Passed straight through to
260 `skimage.registration.phase_cross_correlation`'s own
261 `upsample_factor` -- images are registered to within
262 `1 / upsample_factor` of a pixel. Default `100`, matching
263 the point past which [Subpixel
264 Accuracy](../getting_started/subpixel_accuracy.html) found
265 diminishing returns (`10` was already close to `100`'s own
266 accuracy in that measurement). Must be >= 1.
268 Returns:
269 The point's location, in `current_image`'s pixel reference
270 frame, generally fractional.
272 Raises:
273 ValueError: If `kernel_margin_width` or `kernel_margin_height` is
274 less than 1, either `search_margin_width`/
275 `search_margin_height` is not greater than its kernel
276 counterpart, or `upsample_factor` is less than 1.
277 """
278 if kernel_margin_width < 1:
279 raise ValueError(f"kernel_margin_width {kernel_margin_width} must be >= 1")
280 if kernel_margin_height < 1:
281 raise ValueError(f"kernel_margin_height {kernel_margin_height} must be >= 1")
282 if search_margin_width <= kernel_margin_width:
283 raise ValueError(
284 f"search_margin_width {search_margin_width} must be greater than "
285 f"kernel_margin_width {kernel_margin_width}"
286 )
287 if search_margin_height <= kernel_margin_height:
288 raise ValueError(
289 f"search_margin_height {search_margin_height} must be greater than "
290 f"kernel_margin_height {kernel_margin_height}"
291 )
292 if upsample_factor < 1:
293 raise ValueError(f"upsample_factor {upsample_factor} must be >= 1")
295 kernel_width = 2 * kernel_margin_width
296 kernel_height = 2 * kernel_margin_height
297 search_width = 2 * search_margin_width
298 search_height = 2 * search_margin_height
300 kernel_origin = PixelCoordinate(
301 x=reference_point.x - kernel_margin_width,
302 y=reference_point.y - kernel_margin_height,
303 )
304 kernel = subimage(
305 image=reference_image,
306 origin=kernel_origin,
307 width=kernel_width,
308 height=kernel_height,
309 )
311 search_origin = PixelCoordinate(
312 x=search_center.x - search_margin_width,
313 y=search_center.y - search_margin_height,
314 )
315 search = subimage(
316 image=current_image,
317 origin=search_origin,
318 width=search_width,
319 height=search_height,
320 )
322 kernel, search = _window(kernel=kernel, search=search, windowing=windowing)
323 kernel_padded, pad_before_height, pad_before_width = _kernel_pad(
324 kernel=kernel, shape=search.shape, centered=True
325 )
327 shift, _error, _phasediff = phase_cross_correlation(
328 reference_image=search,
329 moving_image=kernel_padded,
330 normalization="phase",
331 upsample_factor=upsample_factor,
332 )
334 return SubpixelCoordinate(
335 x=search_origin.x + shift[1] + kernel_margin_width + pad_before_width,
336 y=search_origin.y + shift[0] + kernel_margin_height + pad_before_height,
337 )