Coverage for src/dictk/strain_window.py: 100%
45 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-10-30 13:51 +0000
1r"""Strain from a least-squares displacement fit over a window of grid points.
3Each grid point gets the affine fit of $u$ and $v$ over the
4$n \times n$ points around it, with $n = 2m + 1$ and $m$ the window's
5`margin`. The fit's slopes are the displacement gradient
6$\boldsymbol{\nabla}_0\,\boldsymbol{u}$ at that point. Near the grid
7boundary the window truncates to the points that exist, and the fit
8solves the full $3 \times 3$ normal equations.
10See [Strain Window](../getting_started/strain_window.html) for the
11derivation this module implements.
12"""
14from collections.abc import Sequence
16import numpy as np
18from dictk.element import deformation_gradient, log_strain
19from dictk.image import PixelCoordinate, SubpixelCoordinate
22def _check(*, count_x: int, count_y: int, margin: int, point_count: int) -> None:
23 if margin < 1:
24 raise ValueError(f"margin {margin} must be >= 1")
25 size = 2 * margin + 1
26 if size > count_x or size > count_y:
27 raise ValueError(
28 f"window of {size}x{size} points (margin {margin}) must fit inside "
29 f"the {count_x}x{count_y} grid"
30 )
31 if point_count != count_x * count_y:
32 raise ValueError(
33 f"got {point_count} points, must be count_x * count_y = {count_x * count_y}"
34 )
37def _windows(*, count_x: int, count_y: int, margin: int) -> list[np.ndarray]:
38 windows = []
39 for j in range(count_y):
40 rows = np.arange(max(j - margin, 0), min(j + margin + 1, count_y))
41 for i in range(count_x):
42 columns = np.arange(max(i - margin, 0), min(i + margin + 1, count_x))
43 windows.append((rows[:, None] * count_x + columns[None, :]).ravel())
44 return windows
47def displacement_gradients(
48 *,
49 reference_points: Sequence[PixelCoordinate | SubpixelCoordinate],
50 current_points: Sequence[PixelCoordinate | SubpixelCoordinate],
51 count_x: int,
52 count_y: int,
53 margin: int,
54) -> np.ndarray:
55 r"""The displacement gradient $\boldsymbol{\nabla}_0\,\boldsymbol{u}$ at every grid point.
57 See [General
58 Window](../getting_started/strain_window.html#general-window) and
59 [Window Edges](../getting_started/strain_window.html#window-edges).
61 Args:
62 reference_points: The grid's reference positions, in
63 [`dictk.grid.generate`](./grid.html#generate)'s row-major
64 order.
65 current_points: The same points' current positions, same order.
66 count_x: Number of grid points along x.
67 count_y: Number of grid points along y.
68 margin: The window's half-width $m$, in points, not counting its
69 center. The window has $2m + 1$ points per side. Must be >= 1.
71 Returns:
72 Shape `(count_x * count_y, 2, 2)`. Entry `k` is
73 `[[du/dX, du/dY], [dv/dX, dv/dY]]` at point `k`, the same layout
74 as [`dictk.element.displacement_gradient`](./element.html#displacement_gradient).
76 Raises:
77 ValueError: If `margin` is less than 1, the $(2m + 1) \times
78 (2m + 1)$ window does not fit inside the grid,
79 `reference_points` does not hold `count_x * count_y` points,
80 or `current_points` does not match its length.
81 """
82 _check(
83 count_x=count_x,
84 count_y=count_y,
85 margin=margin,
86 point_count=len(reference_points),
87 )
88 if len(current_points) != len(reference_points):
89 raise ValueError(
90 f"current_points length ({len(current_points)}) must match "
91 f"reference_points length ({len(reference_points)})"
92 )
94 reference = np.array([[p.x, p.y] for p in reference_points], dtype=float)
95 displacements = (
96 np.array([[p.x, p.y] for p in current_points], dtype=float) - reference
97 )
99 gradients = np.empty((len(reference), 2, 2))
100 for k, window in enumerate(
101 _windows(count_x=count_x, count_y=count_y, margin=margin)
102 ):
103 offsets = reference[window] - reference[k]
104 design = np.column_stack([np.ones(len(window)), offsets])
105 coefficients = np.linalg.lstsq(design, displacements[window], rcond=None)[0]
106 gradients[k] = coefficients[1:].T
107 return gradients
110def log_strains(
111 *,
112 reference_points: Sequence[PixelCoordinate | SubpixelCoordinate],
113 current_points: Sequence[PixelCoordinate | SubpixelCoordinate],
114 count_x: int,
115 count_y: int,
116 margin: int,
117) -> np.ndarray:
118 r"""The logarithmic strain $\boldsymbol{E}^{(0)} = \ln\boldsymbol{U}$ at every grid point.
120 Passes each point's `displacement_gradients` entry through
121 [`dictk.element.deformation_gradient`](./element.html#deformation_gradient)
122 and [`dictk.element.log_strain`](./element.html#log_strain). See
123 [From Gradient to
124 Strain](../getting_started/strain_window.html#from-gradient-to-strain).
126 Args:
127 reference_points: The grid's reference positions, in
128 [`dictk.grid.generate`](./grid.html#generate)'s row-major
129 order.
130 current_points: The same points' current positions, same order.
131 count_x: Number of grid points along x.
132 count_y: Number of grid points along y.
133 margin: The window's half-width $m$, in points, not counting its
134 center. Must be >= 1.
136 Returns:
137 Shape `(count_x * count_y, 2, 2)`, one strain tensor per point.
139 Raises:
140 ValueError: Same conditions as `displacement_gradients`.
141 """
142 return np.array(
143 [
144 log_strain(
145 deformation_gradient=deformation_gradient(displacement_gradient=g)
146 )
147 for g in displacement_gradients(
148 reference_points=reference_points,
149 current_points=current_points,
150 count_x=count_x,
151 count_y=count_y,
152 margin=margin,
153 )
154 ]
155 )
158def std_factors(*, count_x: int, count_y: int, margin: int) -> np.ndarray:
159 r"""Each point's gradient standard deviation, relative to a full interior window.
161 Assumes independent tracking errors with equal variance. A full
162 window gives `1.0`. A window truncated by the grid boundary gives
163 more. On a uniformly spaced grid, an edge point at $m = 7$ gives
164 2.58 for the gradient across the edge and a corner point gives
165 3.54. See [Window
166 Edges](../getting_started/strain_window.html#window-edges).
168 Correlated tracking error can exceed these factors. A full window
169 cancels most of a periodic error, and a truncated window does not.
170 On `dictk.grid.locate_warp` positions under a 2% stretch, whose
171 error follows a 51 px wave, the measured $E_{11}$ error at an edge
172 point is 4.03 times the interior's, against a predicted 2.58. Three
173 columns in, it is 6.04 times, against a predicted 1.60. See [Edges
174 and Corners](../getting_started/strain_window.html#edges-and-corners).
176 Args:
177 count_x: Number of grid points along x.
178 count_y: Number of grid points along y.
179 margin: The window's half-width $m$, in points. Must be >= 1.
181 Returns:
182 Shape `(count_x * count_y, 2)`, in row-major order. Column 0 is
183 the factor for $\partial / \partial X$, column 1 for
184 $\partial / \partial Y$.
186 Raises:
187 ValueError: If `margin` is less than 1 or the $(2m + 1) \times
188 (2m + 1)$ window does not fit inside the grid.
189 """
190 _check(
191 count_x=count_x,
192 count_y=count_y,
193 margin=margin,
194 point_count=count_x * count_y,
195 )
197 indices = np.array(
198 [[i, j] for j in range(count_y) for i in range(count_x)], dtype=float
199 )
200 variances = np.empty((len(indices), 2))
201 for k, window in enumerate(
202 _windows(count_x=count_x, count_y=count_y, margin=margin)
203 ):
204 offsets = indices[window] - indices[k]
205 design = np.column_stack([np.ones(len(window)), offsets])
206 variances[k] = np.diag(np.linalg.inv(design.T @ design))[1:]
208 interior = margin * count_x + margin
209 return np.sqrt(variances / variances[interior])