Coverage for src/dictk/translation.py: 100%
25 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-27 23:53 +0000
1"""Point translation tracking between a reference and current image."""
3import numpy as np
4from skimage.registration import phase_cross_correlation
6from dictk.image import PixelCoordinate, subimage
9def locate(
10 *,
11 reference_image: np.ndarray,
12 current_image: np.ndarray,
13 reference_point: PixelCoordinate,
14 search_center: PixelCoordinate,
15 kernel_margin_width: int,
16 kernel_margin_height: int,
17 search_margin_width: int,
18 search_margin_height: int,
19) -> PixelCoordinate:
20 """Given a `reference_point` expressed in the `reference_image`
21 frame, find its position expressed in the `current_image` frame.
23 `reference_point` is a fixed, already-known point in `reference_image`.
24 Where that same physical point ends up in `current_image` is exactly
25 what this function finds — it is *not* an input, it's the return
26 value. `search_center` is a different thing: just where to center the
27 search area in `current_image`, i.e., a guess of roughly where to
28 look. It doesn't need to be exact — only close enough that the true
29 (unknown) position falls within the search area — and if no better
30 guess is available, passing `reference_point` again is a reasonable
31 default.
33 Concretely: extracts a rectangular kernel (also called a subset,
34 filter, or convolution matrix; `2 * kernel_margin_width` wide,
35 `2 * kernel_margin_height` tall) from `reference_image` centered at
36 `reference_point`, and a rectangular search area (also called a
37 search window, scanning zone, or area of interest (AOI);
38 `2 * search_margin_width` wide, `2 * search_margin_height` tall)
39 from `current_image` centered at `search_center`, then locates the
40 kernel within the search area via FFT-based phase cross-correlation
41 (`skimage.registration.phase_cross_correlation`). The kernel always
42 comes from `reference_image`; the search area always comes from
43 `current_image`. Integer-pixel precision only; subpixel refinement is
44 out of scope for now.
46 Args:
47 reference_image: The reference (undeformed) 2D grayscale image.
48 current_image: The current (deformed) 2D grayscale image.
49 reference_point: The point's fixed, known position, in
50 `reference_image`'s pixel reference frame.
51 search_center: Where to center the search area, in
52 `current_image`'s pixel reference frame — a guess of roughly
53 where `reference_point` ended up, not the answer itself.
54 kernel_margin_width: Half the kernel's width, in pixels. Must
55 be >= 1.
56 kernel_margin_height: Half the kernel's height, in pixels. Must
57 be >= 1.
58 search_margin_width: Half the search area's width, in pixels.
59 Must be greater than `kernel_margin_width`.
60 search_margin_height: Half the search area's height, in pixels.
61 Must be greater than `kernel_margin_height`.
63 Returns:
64 The point's location, in `current_image`'s pixel reference frame.
66 Raises:
67 ValueError: If `kernel_margin_width` or `kernel_margin_height` is
68 less than 1, or either `search_margin_width`/
69 `search_margin_height` is not greater than its kernel
70 counterpart.
71 """
72 if kernel_margin_width < 1:
73 raise ValueError(f"kernel_margin_width {kernel_margin_width} must be >= 1")
74 if kernel_margin_height < 1:
75 raise ValueError(f"kernel_margin_height {kernel_margin_height} must be >= 1")
76 if search_margin_width <= kernel_margin_width:
77 raise ValueError(
78 f"search_margin_width {search_margin_width} must be greater than "
79 f"kernel_margin_width {kernel_margin_width}"
80 )
81 if search_margin_height <= kernel_margin_height:
82 raise ValueError(
83 f"search_margin_height {search_margin_height} must be greater than "
84 f"kernel_margin_height {kernel_margin_height}"
85 )
87 kernel_width = 2 * kernel_margin_width
88 kernel_height = 2 * kernel_margin_height
89 search_width = 2 * search_margin_width
90 search_height = 2 * search_margin_height
92 kernel_origin = PixelCoordinate(
93 x=reference_point.x - kernel_margin_width,
94 y=reference_point.y - kernel_margin_height,
95 )
96 kernel = subimage(
97 image=reference_image,
98 origin=kernel_origin,
99 width=kernel_width,
100 height=kernel_height,
101 )
103 search_origin = PixelCoordinate(
104 x=search_center.x - search_margin_width,
105 y=search_center.y - search_margin_height,
106 )
107 search = subimage(
108 image=current_image,
109 origin=search_origin,
110 width=search_width,
111 height=search_height,
112 )
114 # phase_cross_correlation requires both images the same shape; only
115 # the kernel needs padding, since the search area is always larger.
116 pad_width = search_width - kernel_width
117 pad_height = search_height - kernel_height
118 kernel_padded = np.pad(kernel, ((0, pad_height), (0, pad_width)))
120 # search=reference_image, kernel_padded=moving_image, not the other
121 # way around: skimage docs define `shift` as "the shift required to
122 # register moving_image with reference_image", and what we actually
123 # want is where the kernel's content sits within the (larger) search
124 # area, so search must be reference_image and kernel_padded must be
125 # moving_image. Swapping them negates the result and returns the
126 # wrong location -- verified empirically, not just by this reasoning.
127 #
128 # normalization="phase" (skimage's own default) is used here
129 # deliberately -- we do NOT follow hdic's registration.py, which uses
130 # normalization=None. hdic's None came from its own pipeline running
131 # a separate ZNCC-style preprocess() step before calling
132 # phase_cross_correlation, specifically to avoid double-normalizing
133 # on top of that. locate() has no such preprocessing step, so that
134 # reasoning doesn't carry over here -- and empirically, sweeping
135 # kernel/search margin ratios against a known displacement,
136 # normalization="phase" matched or outperformed normalization=None,
137 # not the other way around.
138 #
139 # What normalization actually does: the core computation is a
140 # cross-power spectrum, image_product = FFT(reference) *
141 # conj(FFT(moving)). normalization="phase" divides this by its own
142 # magnitude at every frequency, discarding contrast/energy
143 # information and keeping only phase -- the classic Kuglin-Hines
144 # phase correlation, giving a sharp peak and robustness to
145 # illumination differences between the two images.
146 # normalization=None skips that division (plain unnormalized
147 # cross-correlation), which skimage's own docs describe as less
148 # robust to noise but sometimes preferable in high-noise scenarios.
149 # Which is better is genuinely content-dependent, not a settled
150 # default-is-always-better situation.
151 shift, _error, _phasediff = phase_cross_correlation(
152 reference_image=search, moving_image=kernel_padded, normalization="phase"
153 )
155 return PixelCoordinate(
156 x=search_origin.x + int(shift[1]) + kernel_margin_width,
157 y=search_origin.y + int(shift[0]) + kernel_margin_height,
158 )