Kernel Warping

High Point Density left a wide strain distribution unexplained. Its 11024 Gauss points ranged from about -16400 to 106100 microstrain. VIC-2D's own measurements of the same stretch stayed between about 17300 and 23100. That page blamed element size: a fixed tracking error, divided by a small element size $L$ , becomes a large strain error.

That explanation leaves out the tracking error's shape. This page measures the error against each point's known true position. The $x$ error follows a wave, with an amplitude of about 0.13 px, tied to where each point's true position falls within its pixel. That wave explains High Point Density's vertical stripes.

This page introduces locate_warp, a subpixel method different from locate_subpixel. It cuts the $x$ error's standard deviation from 0.109 px to 0.008 px. Relative to the locate_subpixel baseline, locate_warp shrinks the strain spread by a factor of 14, from a standard deviation of 16531 to 1161 microstrain.

Pixel Locking

This page's 2% stretch sends every reference point to a new known $x$ location, $x = 1.02 X$ ( $y = Y$ remains constant for all points). The High Point Density page's own locate_subpixel call already tracked all 2862 of those points. Subtracting each point's known true position from its tracked one gives the tracking error:

x error: mean -0.0110 px, std 0.1091 px
y error: mean +0.0005 px, std 0.0385 px

Each row's x error vs. the all-row mean: correlation 0.82 to 0.94; row 27 (y = 151 px): 0.82

Saved: kernel_warping_pixel_locking.png

two panels. Left: x error along one grid row against true x, oscillating between about -0.18 and +0.17 px, with dotted vertical lines every 51 px marking where each cycle begins. Right: x error of every tracked point plotted against the fractional part of its true x position, gray dots with black binned means tracing one full sine-like cycle, about -0.14 px near 0.3 and +0.12 px near 0.7, crossing zero near 0.5 — **Figure 9.3-1:** Left: each point's $x$ error along row 27 of the grid, at $y = 151$ px. Right: the same error for all points, against the fractional part of the true $x$ (gray), with the mean of each tenth of a pixel (black). Row 27 is the middle of the grid's 54 rows, chosen to stay clear of the image's top and bottom edges. The wave depends only on $x$ , so every row shows it. Each row's error correlates 0.82 to 0.94 with the average over all rows. Row 27 sits at the low end, 0.82. So the left panel shows the least regular of the 54 rows. Dotted lines mark every 51 px, the distance over which the true position's fractional part completes one cycle.

The right panel shows the error depends on where the true position sits within its pixel. Points with a fractional part between 0 and 0.5 read low, by as much as 0.14 px near 0.3. Points between 0.5 and 1 read high, by as much as 0.12 px near 0.7. Points at an exact half pixel read about right. This is pixel locking: a subpixel estimator's answer drifts toward whole-pixel values. locate_subpixel refines the peak of an upsampled phase correlation surface. That surface comes from a zero-padded kernel, and its peak shifts in a way that depends on the true fractional position. About half of this wave's size comes from how the stretched image itself was made: image.stretch uses bilinear interpolation. The Residual Wave separates the two causes.

The left panel shows why this matters for strain. Every reference position $X$ on the grid is a whole pixel. So the true position, $x = 1.02 X = X + 0.02 X$ , has the same fractional part as $0.02 X$ . That fractional part climbs from 0 to 1, then wraps back to 0, each time $0.02 X$ grows by 1. That takes $Δ X = 1/0.02 = 50$ reference pixels. The stretch scales those 50 px by 1.02, so one cycle spans $1.02 \times 50 = 51$ pixels in the current image. So the bias itself forms a wave along $x$ , with an amplitude of about 0.13 px and a period of 51 px. Strain measures the slope of displacement. The slope of that wave swings to about $\pm 0.13 \cdot 2 π /51 \approx \pm 0.016$ , or 16000 microstrain. That is the same order as High Point Density's measured standard deviation, 16531 microstrain. The same wave explains the vertical stripes in that page's field figure, one stripe pair per cycle.

It also explains why Simple Stretch Revisited never saw this. That page spaced its points 50 px apart, one full bias cycle. Every true position landed on a whole pixel, so every node sat at the same point on the wave. High Point Density's 5 px spacing spread the nodes across the whole wave, and the error showed up in every element.

Inverse Compositional Gauss-Newton

Commercial DIC codes such as VIC-2D warp the kernel itself, instead of refining a correlation peak (as done by locate_subpixel). dictk.grid.locate_warp warps the kernel too, as VIC-2D does. It works in two stages:

First, dictk.translation.locate finds the nearest whole pixel.
Second, inverse compositional Gauss-Newton (IC-GN) refines it.

IC-GN lets the reference kernel deform by an affine warp, $W$ , with six parameters: $u$ , $u_{x}$ , $u_{y}$ , $v$ , $v_{x}$ , and $v_{y}$ . $u$ and $v$ are the displacements along $x$ and $y$ . The subscripts mark their gradients across the kernel, so $u_{x} = \partial u / \partial x$ and $v_{y} = \partial v / \partial y$ :

$W (Δ x, Δ y; p) = 1 + u_{x} v_{x} 0 u_{y} 1 + v_{y} 0 u v 1 Δ x Δ y 1, p = (u, u_{x}, u_{y}, v, v_{x}, v_{y}) .$

Here $(Δ x, Δ y)$ labels each pixel of the reference kernel by its fixed offset from the kernel's center, from $- 13$ to $+ 13$ px here. The warp maps each such pixel to its position in the current image. IC-GN searches for the $p$ that minimizes the zero-normalized sum of squared differences (ZNSSD) between the reference kernel $f$ and the current image sampled through the warp, $g = g (W)$ :

$C (p) = Δ x, Δ y \sum [\frac{f - f ˉ}{∥ f - f ˉ ∥} - \frac{g - g ˉ}{∥ g - g ˉ ∥}]^{2} .$

Here $\overset{ˉ}{f}$ and $\overset{g}{ˉ}$ are the mean intensities over the 729 pixels of $f$ and $g$ (the kernel spans 13 px on each side of its center pixel, so $(2 \cdot 13 + 1)^{2} = 729$ ). The double bars $∥ \cdot ∥$ denote the Euclidean norm of those pixel values.

Each iteration solves for a small update $Δ p$ to the reference kernel. It then composes that update's inverse into the current warp, $W \leftarrow W \circ W (Δ p)^{- 1}$ . The inverse compositional form puts all gradient work on $f$ . It computes the Hessian once per point, not once per iteration. Cubic B-splines sample $g$ at the warped, fractional positions. The tracked position is the warp's own translation, $(u, v)$ .

Deriving the Update

Four steps lead from $C (p)$ to the update IC-GN applies. dictk.warp's refinement follows them in order. Below, $ξ = (Δ x, Δ y)$ names one pixel's offset in the kernel, and every sum runs over all $(2 \cdot 13 + 1)^{2} = 729$ pixels.

1. Reduce ZNSSD to least squares. Write $f = f - \overset{ˉ}{f}$ and $g = g - \overset{g}{ˉ}$ for the zero-mean kernels. Then write $a = f / ∥ f ∥$ and $b = g / ∥ g ∥$ for the same kernels scaled to unit length, so $\sum a^{2} = 1$ and $\sum b^{2} = 1$ . The bracket in $C$ is $a - b$ . Its square expands to $(a - b)^{2} = a^{2} - 2 ab + b^{2}$ . Summing over the pixels gives

$C (p) = 1 \sum a^{2} - 2 \sum a b + 1 \sum b^{2} = 2 - 2 \frac{\sum f g}{∥ f ∥ ∥ g ∥} = 2 [1 - ZNCC (p)] .$

So minimizing $C$ maximizes ZNCC, the criterion from Correlation Criteria. ZNCC runs from $- 1$ to $+ 1$ , because $\sum a b$ is the dot product of two unit-length vectors. So $C = 2 [1 - ZNCC]$ runs from 0 to 4. $C = 0$ is a perfect match, where $ZNCC = 1$ . $C = 4$ is a perfectly inverted pattern, where $ZNCC = - 1$ .

Next, let $∥ f ∥^{2} = \sum_{ξ} f (ξ)^{2}$ , the sum of the squared zero-mean reference intensities. It depends only on the reference kernel, not on $p$ , so it is a constant. Multiplying $C$ by it leaves the minimizer unchanged:

$∥ f ∥^{2} C (p) = ξ \sum [f (ξ) - λ g (ξ; p)]^{2}, λ = \frac{∥ f ~ ∥}{∥ g ∥} .$

2. Put the increment on the reference kernel. A forward, additive Gauss-Newton step would perturb the current side, $g (W (ξ; p + Δ p))$ . Its linearization needs the gradient of $g$ at the warped positions. Those positions move every iteration, so the 6x6 Hessian would need rebuilding every iteration too. The inverse compositional form perturbs $f$ instead. It seeks the incremental warp $W (ξ; Δ p)$ that makes the reference kernel match the current one:

$Δ p min ξ \sum [f (W (ξ; Δ p)) - λ g (ξ; p)]^{2} .$

If $f \circ W (Δ p)$ matches $g \circ W (p)$ , then $f$ matches $g \circ W (p) \circ W (Δ p)^{- 1}$ . That gives the update, $W (p) \leftarrow W (p) \circ W (Δ p)^{- 1}$ . In the $3 \times 3$ matrix form above, composition is matrix multiplication. At convergence, $Δ p = 0$ , and the two forms share the same minimum.

3. Linearize about $Δ p = 0$ . $W (ξ; 0)$ is the identity, so a first-order Taylor expansion gives

$f (W (ξ; Δ p)) \approx f (ξ) + \nabla f (ξ) \frac{\partial W}{\partial p} Δ p, \frac{\partial W}{\partial p} = [10 Δ x 0 Δ y 0 01 0 Δ x 0 Δ y],$

with $\nabla f = (f_{x}, f_{y})$ , the reference image's gradient. The product is one row of six numbers per pixel, its steepest-descent row:

$s (ξ) = [f_{x}, f_{x} Δ x, f_{x} Δ y, f_{y}, f_{y} Δ x, f_{y} Δ y] .$

The expansion holds $\overset{ˉ}{f}$ and $∥ \tilde{f} ∥$ fixed. That approximation vanishes as $Δ p \to 0$ .

Figure 9.3-2 shows the six components of $s$ as images, for the reference kernel at $(150, 150)$ . Rigid Versus Warped Kernel below tracks that same kernel. The $u$ and $v$ images are the image gradient itself. The other four multiply that gradient by $Δ x$ or $Δ y$ . So they vanish along the kernel's center column or center row. A pixel 13 px from the center weighs 13 times more than a pixel 1 px from it. The kernel's edges carry the information about stretch and shear.

six square heat maps, 27 by 27 pixels each, in two rows of three, colored blue for negative and red for positive with white at zero. Top row: f_x, f_x times delta x, f_x times delta y. Bottom row: f_y, f_y times delta x, f_y times delta y. The f_x and f_y panels show speckle-scale structure everywhere. The four weighted panels fade to white along the center row or column and grow strongest toward the kernel edges — **Figure 9.3-2:** The six steepest-descent images for the reference kernel at $(150, 150)$ , one per component of $p = (u, u_{x}, u_{y}, v, v_{x}, v_{y})$ . Each panel uses its own symmetric color scale. The Hessian $H$ sums products of these six images, pixel by pixel.

4. Solve the normal equations. Define the residual $r (ξ) = λ g (ξ; p) - f (ξ)$ . Substituting step 3 into step 2 leaves a linear least-squares problem,

$Δ p min ξ \sum [s (ξ) Δ p - r (ξ)]^{2} .$

Setting its gradient with respect to $Δ p$ to zero gives

$H Δ p = ξ \sum s (ξ)^{⊤} r (ξ), H = ξ \sum s (ξ)^{⊤} s (ξ) .$

$H$ is $6 \times 6$ . It depends only on $f$ and $ξ$ , never on $p$ . So IC-GN builds and inverts it once per point, before the first iteration. That is the payoff of step 2.

Iterating. The first $p$ comes from locate's whole-pixel answer: $u$ and $v$ take its displacement, and all four gradient terms start at zero. Each iteration then does five things:

Sample $g$ at the 729 warped positions with cubic B-splines.
Form $\tilde{g}$ , $λ$ , and the residual $r$ .
Solve for $Δ p = H^{- 1} \sum s^{⊤} r$ .
Build $W (Δ p)$ , invert that $3 \times 3$ matrix, and compose it into $W (p)$ .
Stop once every component of $Δ p$ falls below $1 0^{- 6}$ , or after 50 iterations.

One example shows those iterations in action: the point at $(150, 150)$ , under a 7% stretch, a 0.05 shear, and a $(6.43, 4.92)$ px shift. Rigid Versus Warped Kernel below explains how that example was chosen. dictk.warp.fit with max_iterations set to $k$ returns the warp after $k$ iterations. Iteration 0 is locate's whole-pixel start. Figure 9.3-3 shows the residual after each iteration, and Figure 9.3-4 the error:

IC-GN reaches its final warp after 9 iterations:

Iteration $k$	ZNSSD $C$	Translation error (px)	Gradient error
0	1.77e-01	4.37e-01	7.00e-02
1	2.20e-03	1.76e-02	4.08e-03
2	1.21e-03	2.69e-03	8.24e-04
3	1.19e-03	3.19e-03	2.88e-04
4	1.19e-03	2.79e-03	1.37e-04
5	1.19e-03	2.87e-03	1.24e-04
6	1.19e-03	2.85e-03	1.17e-04
7	1.19e-03	2.86e-03	1.16e-04
8	1.19e-03	2.85e-03	1.16e-04
9	1.19e-03	2.86e-03	1.16e-04

four square heat maps of the residual r, 27 by 27 pixels each, on one shared blue-white-red color scale from about -130 to 130. At k = 0, with C = 0.18, strong red and blue speckles cover the kernel. At k = 1, 2, and 3, with C = 0.0022, 0.0012, and 0.0012, the maps are almost entirely white — **Figure 9.3-3:** The residual $r = λ \tilde{g} - \tilde{f}$ after $k$ iterations, on one shared color scale, with each panel's ZNSSD $C$ in its title. At $k = 0$ , the whole-pixel start, speckle structure fills the kernel. One iteration removes it.

two line plots against iteration k from 0 to 9, both with logarithmic vertical axes. Left, translation error: 0.44 px at k = 0, 0.018 px at k = 1, then flat near 0.003 px from k = 2 on. Right, gradient error: 0.07 at k = 0, falling steadily to about 0.00012 by k = 6, then flat — **Figure 9.3-4:** Error against iteration $k$ . Left: the distance from the fitted translation to the true shift $(6.43, 4.92)$ . Right: the largest difference between the fitted gradient terms and the true deformation, $1.07$ stretch and $0.05$ shear.

One iteration cuts $C$ from 0.18 to 0.0022, and the translation error from 0.44 px to 0.018 px. By $k = 2$ , $C$ levels off at 0.0012, and the translation error settles at about 0.003 px. The gradient terms take until $k = 6$ to settle, at about $1.2 \times 1 0^{- 4}$ . IC-GN stops after 9 iterations, the first one where every parameter's update falls below the default tolerance of $1 0^{- 6}$ . The error stops improving at $k = 6$ , but the updates take three more iterations to shrink below that tolerance. The 50-iteration cap never applies.

Neither error reaches zero. This example's current image comes from a quintic spline, while IC-GN samples it with a cubic one. That mismatch sets the 0.003 px floor. Rigid Versus Warped Kernel below finds the same 0.003 px as locate_warp's median error across 30 shifts.

Rigid Versus Warped Kernel

locate and locate_subpixel never warp their kernel. Each cuts an axis-aligned square from reference_image. Each slides that square, unchanged, across the search area. The only thing either one finds is a translation. In contrast, locate_warp finds how the kernel translates and deforms.

One point, at $(150, 150)$ , under an exaggerated deformation shows the difference. The current image stretches 7% along $x$ and shears by 0.05. The 7% stretch is 3.5 times this page's 2%, large enough to see the warp.

A shift then moves the point by a fractional number of pixels. A whole-pixel shift would hand locate an exact answer by construction, since locate only returns whole pixels. So the first table tracks the point under 30 random fractional shifts, each between 4 and 8 px along $x$ and between 2 and 6 px along $y$ . Figure 9.3-5 uses the one shift whose three errors sit closest to the three median errors in that table. The second table tracks that shift. Then come the warp matrices, $W$ , from locate_subpixel and locate_warp, followed by the true one:

Across 30 random fractional shifts:

Function	Smallest error (px)	Median error (px)	Largest error (px)
`locate`	0.081	0.462	0.658
`locate_subpixel`	0.054	0.105	0.150
`locate_warp`	0.001	0.003	0.006

The most typical shift, (6.43, 4.92) px, puts the true position at (156.43, 154.92):

Function	Found $x$	Found $y$	Error (px)
`locate`	156.000	155.000	0.437
`locate_subpixel`	156.340	154.970	0.103
`locate_warp`	156.432	154.918	0.003

$W_{locate_subpixel} = 100010 6.3400 4.9700 1 (u_{x} = u_{y} = v_{x} = v_{y} = 0)$

$W_{locate_warp} = 1.0701 - 0.0001 0 0.0499 1.0000 0 6.4318 4.9178 1$

$W_{true} = 1.07 00 0.05 10 6.43 4.92 1$

two zoomed panels of the same deformed speckle image around one point. In both, a dashed white square marks the reference kernel at (150, 150), an orange arrow points to the tracked position, and a red cross marks the true position at (156.43, 154.92). Left: an orange axis-aligned square with a 10x10 grid of orange dots, the same shape as the reference square, moved by the arrow. Right: an orange parallelogram, slightly wider than the square and leaning right toward the bottom, with its 10x10 grid of dots stretched and sheared to match — **Figure 9.3-5:** Both panels show the deformed image around one point. The dashed white square marks the reference kernel at its original position. The orange arrow points to each function's tracked position, and the red cross marks the true one. Left: `locate` and `locate_subpixel` move the square without changing its shape (drawn at `locate_subpixel`'s position). Right: `locate_warp` also stretches and shears the kernel. Its outline and sample points use the warp IC-GN fitted, $W_{locate_warp}$ .

Across all 30 shifts, locate misses by a median 0.462 px, since it rounds to a whole pixel. locate_subpixel misses by a median 0.105 px, because its rigid square no longer matches the stretched and sheared content. locate_warp misses by a median 0.003 px, about 35 times smaller than locate_subpixel's. Even locate_warp's largest miss over all 30 shifts, 0.006 px, is about 9 times smaller than locate_subpixel's smallest, 0.054 px.

$W_{locate_subpixel}$ keeps an identity block in its upper left. A rigid kernel can only translate, so its four gradient terms stay at zero. $W_{locate_warp}$ is the warp IC-GN fitted. dictk.warp.fit returns it directly. Its gradient block matches $W_{true}$ to about 0.0001. Its translation matches to within 0.003 px.

Two things reduce the bias. The spline samples the image itself at fractional positions, so no correlation peak needs interpolating. The affine warp also lets the kernel stretch along with the material, instead of forcing a rigid match. The call mirrors locate_subpixel, with the same grid and margins:

Show the locate_warp call

from dictk.grid import locate_warp

found = locate_warp(
    reference_image=reference_image,
    current_image=current_image,
    reference_points=points,
    kernel_margin_width=13,
    kernel_margin_height=13,
    search_margin_width=25,
    search_margin_height=25,
)

The same fractional-position error curve, now for both trackers:

Tracker	Time (s)	Mean $x$ error (px)	Std $x$ error (px)	Max \| $x$ error\| (px)
`locate_subpixel`	1.6	-0.0110	0.1091	0.3200
`locate_warp`	1.7	-0.0000	0.0080	0.0234

mean x error against the fractional part of true x for two trackers. locate_subpixel, circles, swings from about -0.14 to +0.12 px. locate_warp, squares, stays within about 0.01 px of zero across the whole pixel — **Figure 9.3-6:** Mean $x$ error per tenth of a pixel, for `locate_subpixel` (circles) and `locate_warp` (squares), across all 2862 points.

The error standard deviation falls from 0.109 px to 0.008 px. The table lists both run times. For all 2862 points they differ by less than a second. A residual wave, about 0.01 px tall, remains. The Residual Wave traces it to the test image.

Figure 9.3-7 shows the 2862 locate_warp errors behind that 0.008 px. They come from the same points as Figure 9.2-1's locate_subpixel errors:

2862 points: mean x error -0.00004 px, std 0.0080 px, range [-0.023, +0.021] px

histogram of dictk.grid.locate_warp's 2862 x errors in pixels, a narrow distribution with two broad lobes either side of zero, spanning roughly -0.023 to +0.021 px, with a black vertical line at the mean of 0.0000 px, a pale blue band one standard deviation of 0.0080 px either side of it, and a dashed red vertical line at exact tracking, 0 px — **Figure 9.3-7:** Distribution of `dictk.grid.locate_warp`'s $x$ error across all 2862 points (gray, 60 bins). The dashed red line marks exact tracking, 0 px. The solid black line marks the mean of the 2862 errors, 0.0000 px. The blue band spans one standard deviation, 0.0080 px, on each side of that mean. Range is [-0.023, +0.021] px.

The mean error is 0.0000 px to four decimals. The standard deviation is 0.0080 px, 13.6 times smaller than locate_subpixel's 0.1090 px. The errors span [-0.023, +0.021] px, against locate_subpixel's [-0.29, +0.32] px.

Strain `locate_warp`

High Point Density's strain pipeline stays the same: dictk.grid.elements for connectivity, then gauss_point_log_strains at each Gauss point. Only the tracked positions change:

Source	Mean (µε)	Std (µε)	Min (µε)	Max (µε)	Inside [17560, 22360]
`locate_subpixel` + Q4	20465	16531	-16446	106134	9.5%
VIC-2D	19876	1385	17330	23137	98.2%
`locate_warp` + Q4	19851	1161	16562	23160	97.4%
Analytical, $ln (1.02)$	19803	0

The spread falls from 16531 microstrain for locate_subpixel to 1161 microstrain for locate_warp. The range now sits where VIC-2D's does. Of the 11024 Gauss points, 97.4% land inside VIC-2D's own colorbar range, up from 9.5%. The standard deviation, 1161 microstrain, comes in under VIC-2D's 1385. The element code is unchanged from High Point Density, so the tracker accounts for the whole drop from 16531 to 1161 microstrain.

VIC-2D's measured exx field, colorbar 17560 to 22360 microstrain, a striped pattern in VIC-2D's own 16-band magenta-blue-cyan-green-yellow-orange-red color scale — **Figure 9.3-8:** VIC-2D's $e_{xx}$ (left) and `dictk`'s $E_{11}$ from `locate_warp` (right). Both use 17560 to 22360 microstrain and VIC-2D's own 16-band color scale. In High Point Density's version of the right panel, only 9.5% of the Gauss points fell inside this range. The rest clipped to solid magenta or solid red. This one uses the whole scale.

dictk's own E11 field from locate_warp, on the same 17560-22360 microstrain range and VIC-2D's 16-band colormap, showing the full range of colors in narrow vertical stripes, with almost no clipped solid magenta or solid red — **Figure 9.3-8:** VIC-2D's $e_{xx}$ (left) and `dictk`'s $E_{11}$ from `locate_warp` (right). Both use 17560 to 22360 microstrain and VIC-2D's own 16-band color scale. In High Point Density's version of the right panel, only 9.5% of the Gauss points fell inside this range. The rest clipped to solid magenta or solid red. This one uses the whole scale.

two overlaid normalized histograms of log strain in microstrain. dictk, gray filled, one broad peak from about 17000 to 23000 centered near 20000. VIC-2D, black outline, five separated clusters over the same 17300 to 23100 range. A dashed red line at 19803 microstrain — **Figure 9.3-9:** Distribution of `dictk`'s $E_{11}$ from `locate_warp` (gray, 11024 Gauss points) and VIC-2D's $e_{xx}$ (black outline, 2682 points), each normalized to unit area. The dashed red line marks $ln (1.02) \approx 19803$ microstrain.

The two distributions now cover the same range. Their shapes still differ. VIC-2D's splits into five clusters. dictk's forms one broad peak. Both fields keep vertical stripes. Both tools measured the same bilinear images, and The Residual Wave below shows those images carry a distortion of their own, tied to each point's fractional position. Each tracker responds to it differently.

The dictk result is not smoothed. Commercial codes usually also apply a strain window, fitting displacement over a neighborhood of points before differentiating. The fit averages the tracking errors of every point in the window. It also blurs any real strain change inside the window into one value. Strain Window derives a strain window and measures it on these tracked positions. Fixing the bias first matters. A window spanning a whole bias cycle would also flatten High Point Density's stripes. That averaging would hide the error while leaving it in the tracked positions.

The Residual Wave

Until now, this page treated the tracker as the source of the error. It compared locate_subpixel with locate_warp, and locate_warp left a 0.01 px wave. This section isolates the other error contributor, the generator: the bilinear interpolation image.stretch uses to build the deformed image. Pixel Locking flagged it earlier. The measurement here separates the two causes.

dictk.image.stretch builds the current image with bilinear interpolation. Bilinear interpolation moves content by a fractional pixel, and it also blurs that content. The blur depends on the fraction. A whole-pixel position gets none, and a half-pixel position gets the most. So the stretched image differs from a true stretch, by an amount that follows each point's fractional position. Any tracker reads that difference as a position-dependent error. The table below measures how large.

A quintic spline has no fraction-dependent blur like bilinear interpolation's. It builds a second version of the same 2% stretch. Tracking both versions shows how much of each tracker's error came from the image:

Current image	Tracker	Std $x$ error (px)	Std $E_{11}$ (µε)
bilinear (`image.stretch`)	`locate_subpixel`	0.1091	16531
bilinear (`image.stretch`)	`locate_warp`	0.0080	1161
quintic spline	`locate_subpixel`	0.0568	10274
quintic spline	`locate_warp`	0.0015	260

With the quintic-spline image, locate_warp's $x$ error falls from 0.008 px to 0.0015 px. Its strain spread falls from 1161 to 260 microstrain. Removing the generator's distortion cuts IC-GN's $x$ error by a factor of 5. locate_subpixel improves too, from 0.109 px to 0.057 px, but its strain spread stays above 10000 microstrain. About half of phase correlation's error came from the generator. The other half belongs to the tracker.

Every other figure on this page keeps image.stretch. Every synthetic image in this book uses bilinear interpolation, not a quintic spline, and VIC-2D measured those same bilinear images. Setting locate_warp on the quintic image (260 microstrain) against VIC-2D on a bilinear image (1385 microstrain) would be an unfair comparison. The two tools would track different images. So the drop from 1161 to 260 microstrain shows how much of locate_warp's spread the generator causes. It does not belong next to VIC-2D's result. The fair comparison stays locate_warp on the bilinear image, 1161 against VIC-2D's 1385 microstrain.

dictk's own tests of fractional translation avoid the generator entirely. They shift a smoothed speckle image by an exact fraction of a pixel, in the Fourier domain. There, locate_warp recovers every tested shift to within 0.005 px, and beats locate_subpixel every time.

Summary

The table collects this page's results, one column per approach. Each column has a smaller $x$ error standard deviation and a smaller $E_{11}$ standard deviation than the column before it.

The table leaves out whole-pixel locate, the baseline. locate rounds every position, so its $x$ error has a standard deviation of 0.30 px on these 2862 points, with the same kernel and search margins as the other columns. Rounding to the nearest integer alone gives $1/ 12 \approx 0.29$ px.

Subpixel Accuracy compares locate with locate_subpixel. VIC-2D is the reference tool. locate_subpixel and locate_warp track the same bilinear image.

The shaded column tracks a different image, a quintic-spline version of the same 2% stretch. It is not comparable with the other columns. VIC-2D, locate_subpixel, and locate_warp all tracked the bilinear images that image.stretch builds. The shaded column shows what locate_warp would reach if the generator did not distort the image: 0.0015 px and 260 microstrain, down from 0.0080 px and 1161. VIC-2D was not run on the quintic image, so the table makes no comparison.

	`locate_subpixel`	VIC-2D	`locate_warp`	`locate_warp`, quintic image
$x$ error Figure	Figure 9.2-1	Figure 8.1-6	Figure 9.3-7
Mean $x$ error (px)	-0.0110	-0.0003	0.0000	+0.0001
Std $x$ error (px)	0.1091	0.0154	0.0080	0.0015
$x$ error range (px)	[-0.29, +0.32]	[-0.039, +0.035]	[-0.023, +0.021]	[-0.004, +0.004]
$E_{11}$ Figure	Figure 9.2-5	Figure 8.1-8	Figure 9.3-9
Mean $E_{11}$ (microstrain)	20465	19876	19851	19796
Std $E_{11}$ (microstrain)	16531	1385	1161	260
$E_{11}$ range (microstrain)	[-16446, 106134]	[17330, 23137]	[16562, 23160]	[18717, 20618]

Each $x$ error is the tracked $x$ minus the true $x = 1.02 X$ . The closed-form strain for this 2% stretch is $E_{11} = ln (1.02) \approx 19803$ microstrain. Each mean $E_{11}$ above can be compared with that value. The locate_subpixel column uses upsample_factor=100, which refines the correlation peak to within 0.01 px. High Point Density uses the same value. The VIC-2D column uses its 2682 valid points for every row. The dictk columns use all 2862 points for $x$ and 11024 Gauss points for $E_{11}$ .

`kernel_warping_pixel_locking.py`

Show the code

"""Pixel locking: `grid.locate_subpixel`'s own x-error against each
point's known true position, at High Point Density's full 53x54 grid.
Plots the error along one row of the grid, and against the true
position's fractional part.
"""

import matplotlib.pyplot as plt
import numpy as np

from dictk.grid import generate, locate_subpixel
from dictk.image import PixelCoordinate, read, stretch

plt.rcParams.update({"font.family": "serif", "mathtext.fontset": "cm"})

FACTOR_X = 1.02

reference_image = read(path="astronaut0.png")
current_image = stretch(arr=reference_image, factor_x=FACTOR_X)
points = generate(
    origin=PixelCoordinate(x=18, y=16),
    count_x=53,
    count_y=54,
    spacing_x=5,
    spacing_y=5,
)
found = locate_subpixel(
    reference_image=reference_image,
    current_image=current_image,
    reference_points=points,
    kernel_margin_width=13,
    kernel_margin_height=13,
    search_margin_width=25,
    search_margin_height=25,
    upsample_factor=100,
)

true_x = np.array([p.x * FACTOR_X for p in points])
error_x = np.array([f.x for f in found]) - true_x
error_y = np.array([f.y - p.y for f, p in zip(found, points)])
fraction = true_x % 1

print(f"* x error: mean {error_x.mean():+.4f} px, std {error_x.std():.4f} px")
print(f"* y error: mean {error_y.mean():+.4f} px, std {error_y.std():.4f} px")

row = 27
rows = error_x.reshape(54, 53)
row_mean = rows.mean(axis=0)
correlations = [np.corrcoef(r, row_mean)[0, 1] for r in rows]
print()
print(
    f"Each row's x error vs. the all-row mean: correlation "
    f"{min(correlations):.2f} to {max(correlations):.2f}; "
    f"row {row} (y = {points[row * 53].y} px): {correlations[row]:.2f}"
)

bins = np.linspace(0, 1, 11)
centers = (bins[:-1] + bins[1:]) / 2
means = [
    error_x[(fraction >= a) & (fraction < b)].mean() for a, b in zip(bins, bins[1:])
]

fig, (row_axes, fraction_axes) = plt.subplots(
    1, 2, figsize=(10, 4), constrained_layout=True
)
row_x = true_x.reshape(54, 53)[row]
row_axes.plot(row_x, error_x.reshape(54, 53)[row], color="black", marker=".")
row_axes.axhline(0, color="red", linestyle="--", linewidth=1)
for boundary in np.arange(np.ceil(row_x[0]), row_x[-1], 50 * FACTOR_X):
    row_axes.axvline(boundary, color="gray", linestyle=":", linewidth=1)
row_axes.set_xlabel(r"true $x$ (px)")
row_axes.set_ylabel(r"$x$ error (px)")

fraction_axes.scatter(fraction, error_x, s=2, color="gray", alpha=0.4)
fraction_axes.plot(centers, means, color="black", marker="o")
fraction_axes.axhline(0, color="red", linestyle="--", linewidth=1)
fraction_axes.set_xlabel(r"fractional part of true $x$ (px)")
fraction_axes.set_ylabel(r"$x$ error (px)")


fig.savefig("kernel_warping_pixel_locking.png", dpi=300)
print()
print("Saved: kernel_warping_pixel_locking.png")

`kernel_warping_panels.py`

Show the code

"""Side-by-side illustration of the kernel `locate`/`locate_subpixel`
match (left) against the kernel `locate_warp` warps (right), on one
point under an exaggerated affine deformation.

First sweeps 30 random fractional shifts under that deformation and
reports each function's median error. The figure then uses the shift
whose three errors sit closest to those medians, so the example is
typical, not hand-picked.
"""

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import Polygon
from scipy.ndimage import affine_transform

from dictk import translation, warp
from dictk.image import PixelCoordinate, read

plt.rcParams.update({"font.family": "serif", "mathtext.fontset": "cm"})

# Deformation gradient in (x, y) order: 7% stretch along x plus 0.05
# shear, exaggerated beyond this page's 2% so the warp is visible, with
# distinct values so each one is easy to spot in the printed matrices.
# 15% with 0.15 shear defeats the whole-pixel starting match for about
# a third of shifts; 7% with 0.05 shear defeats none of the 30 below.
F = np.array([[1.07, 0.05], [0.0, 1.0]])
POINT = PixelCoordinate(x=150, y=150)
MARGIN = 13
SEARCH_MARGIN = 25
SHIFT_COUNT = 30
DOT_COUNT = 10
ZOOM = 34

reference_image = read(path="astronaut0.png")
pivot = np.array([POINT.x, POINT.y], dtype=np.float64)
# affine_transform maps each output (row, column) back to its input
# position, so it takes the inverse deformation, in (y, x) order.
inverse = np.linalg.inv(F)[::-1, ::-1]


def deformed(shift):
    offset = pivot[::-1] - inverse @ (pivot[::-1] + shift[::-1])
    return affine_transform(
        reference_image.astype(np.float64), inverse, offset=offset, order=5
    )


def tracked(current_image):
    kwargs = dict(
        reference_image=reference_image,
        current_image=current_image,
        reference_point=POINT,
        search_center=POINT,
        kernel_margin_width=MARGIN,
        kernel_margin_height=MARGIN,
        search_margin_width=SEARCH_MARGIN,
        search_margin_height=SEARCH_MARGIN,
    )
    return (
        translation.locate(**kwargs),
        translation.locate_subpixel(**kwargs),
        warp.locate(**kwargs),
        warp.fit(**kwargs),
    )


names = ["`locate`", "`locate_subpixel`", "`locate_warp`"]
# Rounded to 0.01 px, so each shift prints exactly as used.
shifts = np.round(
    np.random.default_rng(0).uniform([4, 2], [8, 6], size=(SHIFT_COUNT, 2)), 2
)
errors = np.array(
    [
        [
            np.hypot(*(np.array([f.x, f.y]) - pivot - shift))
            for f in tracked(deformed(shift))[:3]
        ]
        for shift in shifts
    ]
)
medians = np.median(errors, axis=0)
print(f"Across {SHIFT_COUNT} random fractional shifts:")
print()
print("| Function | Smallest error (px) | Median error (px) | Largest error (px) |")
print("|---|---|---|---|")
for name, smallest, median, largest in zip(
    names, errors.min(axis=0), medians, errors.max(axis=0)
):
    print(f"| {name} | {smallest:.3f} | {median:.3f} | {largest:.3f} |")

typical = int(np.argmin((np.abs(errors - medians) / medians).sum(axis=1)))
TRANSLATION = shifts[typical]
current_image = deformed(TRANSLATION)
true_position = pivot + TRANSLATION
found_locate, found_subpixel, found_warp, fitted = tracked(current_image)

print()
print(
    f"The most typical shift, ({TRANSLATION[0]:.2f}, {TRANSLATION[1]:.2f}) px, "
    f"puts the true position at ({true_position[0]:.2f}, {true_position[1]:.2f}):"
)
print()
print("| Function | Found $x$ | Found $y$ | Error (px) |")
print("|---|---|---|---|")
for name, found in zip(names, [found_locate, found_subpixel, found_warp]):
    error = np.hypot(found.x - true_position[0], found.y - true_position[1])
    print(f"| {name} | {found.x:.3f} | {found.y:.3f} | {error:.3f} |")

rigid = np.array(
    [
        [1.0, 0.0, found_subpixel.x - POINT.x],
        [0.0, 1.0, found_subpixel.y - POINT.y],
        [0.0, 0.0, 1.0],
    ]
)
true = np.eye(3)
true[:2, :2] = F
true[:2, 2] = TRANSLATION


def latex(matrix, exact):
    """`matrix` as a LaTeX bmatrix. Entries where `exact` is True are
    fixed by construction and print as given (`1`, `0`, `1.15`), not
    padded to four decimals."""
    rows = [
        " & ".join(
            f"{value:g}" if fixed else f"{value:.4f}"
            for value, fixed in zip(row, fixed_row)
        )
        for row, fixed_row in zip(matrix, exact)
    ]
    return "\\begin{bmatrix}" + " \\\\ ".join(rows) + "\\end{bmatrix}"


bottom_row = np.zeros((3, 3), dtype=bool)
bottom_row[2] = True
rigid_exact = bottom_row.copy()
rigid_exact[:2, :2] = True
for name, matrix, exact, note in [
    (r"locate\_subpixel", rigid, rigid_exact, r" \quad (u_x = u_y = v_x = v_y = 0)"),
    (r"locate\_warp", fitted, bottom_row, ""),
    ("true", true, np.ones((3, 3), dtype=bool), ""),
]:
    print()
    print("$$")
    print(f"\\boldsymbol{{W}}_\\text{{{name}}} = {latex(matrix, exact)}{note}")
    print("$$")

offsets = np.linspace(-MARGIN, MARGIN, DOT_COUNT)
grid = np.array([(dx, dy) for dy in offsets for dx in offsets])
square = np.array(
    [[-MARGIN, -MARGIN], [MARGIN, -MARGIN], [MARGIN, MARGIN], [-MARGIN, MARGIN]]
)

fig, (left, right) = plt.subplots(1, 2, figsize=(10, 5.2), constrained_layout=True)
for ax in (left, right):
    ax.imshow(current_image, cmap="gray", vmin=0, vmax=255)
    ax.set_xlim(pivot[0] - ZOOM, pivot[0] + ZOOM)
    ax.set_ylim(pivot[1] + ZOOM, pivot[1] - ZOOM)
    ax.set_xlabel("x (pixels)")
    ax.set_ylabel("y (pixels)")
    ax.add_patch(
        Polygon(
            pivot + square,
            fill=False,
            edgecolor="white",
            linestyle="--",
            linewidth=1.5,
        )
    )
    ax.plot(*true_position, marker="+", color="red", markersize=14, mew=2)

for ax, found, corners, dots, title in [
    (
        left,
        found_subpixel,
        square,
        grid,
        "locate, locate_subpixel: rigid square kernel",
    ),
    (
        right,
        found_warp,
        square @ fitted[:2, :2].T,
        grid @ fitted[:2, :2].T,
        "locate_warp: affine-warped kernel",
    ),
]:
    center = np.array([found.x, found.y])
    ax.annotate(
        "",
        xy=center,
        xytext=pivot,
        arrowprops=dict(arrowstyle="-|>", color="orange", linewidth=2),
    )
    ax.add_patch(Polygon(center + corners, fill=False, edgecolor="orange", linewidth=2))
    ax.scatter(*(center + dots).T, s=8, color="orange")
    ax.set_title(title)

fig.savefig("kernel_warping_panels.png", dpi=300)

`kernel_warping_steepest.py`

Show the code

"""The six steepest-descent images of IC-GN's update, for the reference
kernel at the point Kernel Warping's panel figure tracks. Each image is
one component of s(xi), in the same order as p = (u, u_x, u_y, v, v_x,
v_y). The Hessian H is built from these six images alone.
"""

import matplotlib.pyplot as plt
import numpy as np

from dictk.image import read

plt.rcParams.update({"font.family": "serif", "mathtext.fontset": "cm"})

X0, Y0 = 150, 150
MARGIN = 13

reference_image = read(path="astronaut0.png").astype(np.float64)
gradient_y, gradient_x = np.gradient(reference_image)
rows = slice(Y0 - MARGIN, Y0 + MARGIN + 1)
columns = slice(X0 - MARGIN, X0 + MARGIN + 1)
fx = gradient_x[rows, columns]
fy = gradient_y[rows, columns]
dy, dx = np.mgrid[-MARGIN : MARGIN + 1, -MARGIN : MARGIN + 1].astype(np.float64)

components = [
    (r"$u$:  $f_x$", fx),
    (r"$u_x$:  $f_x\,\Delta x$", fx * dx),
    (r"$u_y$:  $f_x\,\Delta y$", fx * dy),
    (r"$v$:  $f_y$", fy),
    (r"$v_x$:  $f_y\,\Delta x$", fy * dx),
    (r"$v_y$:  $f_y\,\Delta y$", fy * dy),
]
extent = (-MARGIN - 0.5, MARGIN + 0.5, MARGIN + 0.5, -MARGIN - 0.5)
fig, axes = plt.subplots(2, 3, figsize=(10, 6.4), constrained_layout=True)
for ax, (title, values) in zip(axes.flat, components):
    limit = np.abs(values).max()
    shown = ax.imshow(values, cmap="RdBu_r", vmin=-limit, vmax=limit, extent=extent)
    ax.set_title(title)
    ax.set_xlabel(r"$\Delta x$ (px)")
    ax.set_ylabel(r"$\Delta y$ (px)")
    fig.colorbar(shown, ax=ax, shrink=0.8)
fig.savefig("kernel_warping_steepest.png", dpi=200)

`kernel_warping_iterations.py`

Show the code

"""IC-GN iteration by iteration, on Kernel Warping's panel-figure example:
the point at (150, 150) under a 7% stretch, 0.05 shear, and the typical
shift (6.43, 4.92) that kernel_warping_panels.py selects. warp.fit with
max_iterations=k returns the warp after k iterations. Iteration 0 is the
whole-pixel start from translation.locate.
"""

import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import affine_transform, map_coordinates

from dictk import translation, warp
from dictk.image import PixelCoordinate, read

plt.rcParams.update({"font.family": "serif", "mathtext.fontset": "cm"})

F = np.array([[1.07, 0.05], [0.0, 1.0]])
SHIFT = np.array([6.43, 4.92])
POINT = PixelCoordinate(x=150, y=150)
MARGIN = 13
SEARCH_MARGIN = 25
RESIDUAL_ITERATIONS = [0, 1, 2, 3]

reference_image = read(path="astronaut0.png").astype(np.float64)
pivot = np.array([POINT.x, POINT.y], dtype=np.float64)
inverse = np.linalg.inv(F)[::-1, ::-1]
offset = pivot[::-1] - inverse @ (pivot[::-1] + SHIFT[::-1])
current_image = affine_transform(reference_image, inverse, offset=offset, order=5)
kwargs = dict(
    reference_image=reference_image,
    current_image=current_image,
    reference_point=POINT,
    search_center=POINT,
    kernel_margin_width=MARGIN,
    kernel_margin_height=MARGIN,
    search_margin_width=SEARCH_MARGIN,
    search_margin_height=SEARCH_MARGIN,
)

dy, dx = np.mgrid[-MARGIN : MARGIN + 1, -MARGIN : MARGIN + 1].astype(np.float64)
f = reference_image[
    POINT.y - MARGIN : POINT.y + MARGIN + 1, POINT.x - MARGIN : POINT.x + MARGIN + 1
]
f_tilde = f - f.mean()


def warp_at(iterations):
    if iterations == 0:
        start = translation.locate(**kwargs)
        return np.array(
            [[1.0, 0.0, start.x - POINT.x], [0.0, 1.0, start.y - POINT.y], [0, 0, 1.0]]
        )
    return warp.fit(**kwargs, max_iterations=iterations)


def residual(matrix):
    xs = POINT.x + matrix[0, 0] * dx + matrix[0, 1] * dy + matrix[0, 2]
    ys = POINT.y + matrix[1, 0] * dx + matrix[1, 1] * dy + matrix[1, 2]
    g = map_coordinates(current_image, [ys, xs], order=3)
    g_tilde = g - g.mean()
    znssd = np.sum(
        (f_tilde / np.linalg.norm(f_tilde) - g_tilde / np.linalg.norm(g_tilde)) ** 2
    )
    scale = np.linalg.norm(f_tilde) / np.linalg.norm(g_tilde)
    return scale * g_tilde - f_tilde, znssd


final = warp.fit(**kwargs)
iterations = [0]
while not np.allclose(warp_at(iterations[-1]), final, rtol=0, atol=1e-12):
    iterations.append(iterations[-1] + 1)
records = []
for k in iterations:
    matrix = warp_at(k)
    r, znssd = residual(matrix)
    translation_error = np.hypot(*(matrix[:2, 2] - SHIFT))
    gradient_error = np.abs(matrix[:2, :2] - F).max()
    records.append((k, znssd, translation_error, gradient_error, r))

print(f"IC-GN reaches its final warp after {iterations[-1]} iterations:")
print()
print("| Iteration $k$ | ZNSSD $C$ | Translation error (px) | Gradient error |")
print("|---|---|---|---|")
for k, znssd, translation_error, gradient_error, _ in records:
    print(f"| {k} | {znssd:.2e} | {translation_error:.2e} | {gradient_error:.2e} |")

shown = [record for record in records if record[0] in RESIDUAL_ITERATIONS]
limit = np.abs(shown[0][4]).max()
fig, axes = plt.subplots(
    1, len(shown), figsize=(10, 3.1), constrained_layout=True, squeeze=False
)
extent = (-MARGIN - 0.5, MARGIN + 0.5, MARGIN + 0.5, -MARGIN - 0.5)
for ax, (k, znssd, _, _, r) in zip(axes.flat, shown):
    image = ax.imshow(r, cmap="RdBu_r", vmin=-limit, vmax=limit, extent=extent)
    ax.set_title(f"$k = {k}$,  $C = {znssd:.2g}$")
    ax.set_xlabel(r"$\Delta x$ (px)")
axes.flat[0].set_ylabel(r"$\Delta y$ (px)")
fig.colorbar(image, ax=axes.flat[-1], shrink=0.85, label="residual $r$")
fig.savefig("kernel_warping_residuals.png", dpi=200)

ks = [record[0] for record in records]
fig, (left, right) = plt.subplots(1, 2, figsize=(10, 3.6), constrained_layout=True)
left.semilogy(ks, [record[2] for record in records], "o-", color="black")
left.set_title("translation error")
left.set_ylabel("error (px)")
right.semilogy(ks, [record[3] for record in records], "o-", color="black")
right.set_title("gradient error")
right.set_ylabel(r"largest $|W_{ij} - F_{ij}|$")
for ax in (left, right):
    ax.set_xlabel("iteration $k$")
    ax.set_xticks(ks)
    ax.grid(True, which="major", color="0.85", linewidth=0.8)
fig.savefig("kernel_warping_convergence.png", dpi=200)

`kernel_warping_bias.py`

Show the code

"""The same fractional-position error curve as
kernel_warping_pixel_locking.py, for both `grid.locate_subpixel` and
`grid.locate_warp`, on the same axes.
"""

import time

import matplotlib.pyplot as plt
import numpy as np

from dictk.grid import generate, locate_subpixel, locate_warp
from dictk.image import PixelCoordinate, read, stretch

plt.rcParams.update({"font.family": "serif", "mathtext.fontset": "cm"})

FACTOR_X = 1.02

reference_image = read(path="astronaut0.png")
current_image = stretch(arr=reference_image, factor_x=FACTOR_X)
points = generate(
    origin=PixelCoordinate(x=18, y=16),
    count_x=53,
    count_y=54,
    spacing_x=5,
    spacing_y=5,
)
kwargs = dict(
    reference_image=reference_image,
    current_image=current_image,
    reference_points=points,
    kernel_margin_width=13,
    kernel_margin_height=13,
    search_margin_width=25,
    search_margin_height=25,
)

true_x = np.array([p.x * FACTOR_X for p in points])
fraction = true_x % 1
bins = np.linspace(0, 1, 11)
centers = (bins[:-1] + bins[1:]) / 2

fig, ax = plt.subplots(figsize=(7, 4), constrained_layout=True)
print(
    "| Tracker | Time (s) | Mean $x$ error (px) | Std $x$ error (px) | Max \\|$x$ error\\| (px) |"
)
print("|---|---|---|---|---|")
for name, locate, extra, marker in [
    ("`locate_subpixel`", locate_subpixel, dict(upsample_factor=100), "o"),
    ("`locate_warp`", locate_warp, {}, "s"),
]:
    start = time.perf_counter()
    found = locate(**kwargs, **extra)
    elapsed = time.perf_counter() - start
    error_x = np.array([f.x for f in found]) - true_x
    print(
        f"| {name} | {elapsed:.1f} | {error_x.mean():+.4f} | {error_x.std():.4f} "
        f"| {np.abs(error_x).max():.4f} |"
    )
    means = [
        error_x[(fraction >= a) & (fraction < b)].mean() for a, b in zip(bins, bins[1:])
    ]
    ax.plot(centers, means, color="black", marker=marker, label=name.strip("`"))
ax.axhline(0, color="red", linestyle="--", linewidth=1)
ax.set_xlabel(r"fractional part of true $x$ (px)")
ax.set_ylabel(r"mean $x$ error (px)")
ax.legend(frameon=False)
fig.savefig("kernel_warping_bias.png", dpi=300)

`kernel_warping_strain.py`

Show the code

"""High Point Density's own strain pipeline, rerun on `grid.locate_warp`
positions: a summary table against `grid.locate_subpixel` and VIC-2D,
the field on VIC-2D's own colorbar, and the distribution against
VIC-2D's own.
"""

import csv

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import ListedColormap

from dictk.element import gauss_point_coordinates, gauss_point_log_strains
from dictk.grid import elements, generate, locate_subpixel, locate_warp
from dictk.image import PixelCoordinate, read, stretch
from dictk.plot import element_strain_plot

plt.rcParams.update({"font.family": "serif", "mathtext.fontset": "cm"})

FACTOR_X = 1.02
VMIN, VMAX = 17560, 22360

reference_image = read(path="astronaut0.png")
current_image = stretch(arr=reference_image, factor_x=FACTOR_X)
points = generate(
    origin=PixelCoordinate(x=18, y=16),
    count_x=53,
    count_y=54,
    spacing_x=5,
    spacing_y=5,
)
kwargs = dict(
    reference_image=reference_image,
    current_image=current_image,
    reference_points=points,
    kernel_margin_width=13,
    kernel_margin_height=13,
    search_margin_width=25,
    search_margin_height=25,
)
element_indices = elements(count_x=53, count_y=54)


def strains(found):
    values, coordinates = [], []
    for element in element_indices:
        current_corners = [found[i] for i in element]
        values.extend(
            strain[0, 0]
            for strain in gauss_point_log_strains(
                reference_points=[points[i] for i in element],
                current_points=current_corners,
            )
        )
        coordinates.extend(gauss_point_coordinates(points=current_corners))
    return np.array(values) * 1e6, coordinates


phase, _ = strains(locate_subpixel(**kwargs, upsample_factor=100))
found = locate_warp(**kwargs)
warp, coordinates = strains(found)

with open("../verification/simple_stretch_vic_out.csv") as f:
    rows = [{k.strip(' "'): v for k, v in row.items()} for row in csv.DictReader(f)]
vic = np.array([float(r["exx"]) * 1e6 for r in rows if float(r["sigma"]) != -1])
analytical = np.log(FACTOR_X) * 1e6

print(
    f"| Source | Mean (µε) | Std (µε) | Min (µε) | Max (µε) | Inside [{VMIN}, {VMAX}] |"
)
print("|---|---:|---:|---:|---:|---:|")
for name, values in [
    ("`locate_subpixel` + Q4", phase),
    ("VIC-2D", vic),
    ("`locate_warp` + Q4", warp),
]:
    inside = np.mean((values >= VMIN) & (values <= VMAX)) * 100
    print(
        f"| {name} | {values.mean():.0f} | {values.std():.0f} | {values.min():.0f} "
        f"| {values.max():.0f} | {inside:.1f}% |"
    )
print(f"| Analytical, $\\ln(1.02)$ | {analytical:.0f} | 0 | | | |")

vic2d_cmap = ListedColormap(
    [
        (0.8314, 0.0000, 1.0000),
        (0.5176, 0.0000, 1.0000),
        (0.1843, 0.0000, 1.0000),
        (0.0000, 0.1333, 1.0000),
        (0.0000, 0.4510, 1.0000),
        (0.0000, 0.7843, 1.0000),
        (0.0000, 1.0000, 0.8980),
        (0.0000, 1.0000, 0.5843),
        (0.0000, 1.0000, 0.2510),
        (0.0667, 1.0000, 0.0000),
        (0.3843, 1.0000, 0.0000),
        (0.7176, 1.0000, 0.0000),
        (1.0000, 0.9686, 0.0000),
        (1.0000, 0.6510, 0.0000),
        (1.0000, 0.3176, 0.0000),
        (1.0000, 0.0000, 0.0000),
    ]
)
element_strain_plot(
    points=found,
    elements=element_indices,
    coordinates=coordinates,
    values=warp,
    label=r"Log Strain, $E_{11}$ (microstrain)",
    image=current_image,
    dot_size=6,
    marker="s",
    show_mesh_lines=False,
    cmap=vic2d_cmap,
    vmin=VMIN,
    vmax=VMAX,
    figsize=(6.9, 6.0),
    path="kernel_warping_strain_vic_colorbar.png",
)

fig, ax = plt.subplots(figsize=(7, 4), constrained_layout=True)
bins = np.linspace(min(warp.min(), vic.min()), max(warp.max(), vic.max()), 61)
ax.hist(
    warp,
    bins=bins,
    density=True,
    color="gray",
    alpha=0.8,
    label="dictk.grid.locate_warp",
)
ax.hist(vic, bins=bins, density=True, histtype="step", color="black", label="VIC-2D")
ax.axvline(analytical, color="red", linestyle="--", linewidth=1.5)
ax.set_xlabel(r"Log strain $E_{11}$ (microstrain)")
ax.set_ylabel("density")
ax.legend(loc="lower center", bbox_to_anchor=(0.5, 1.0), ncol=2, frameon=False)
fig.savefig("kernel_warping_strain_histogram.png", dpi=300)

`kernel_warping_generator.py`

Show the code

"""The same 2% stretch, generated two ways: `dictk.image.stretch`'s own
bilinear interpolation, and a quintic spline. Tracks both with
`grid.locate_subpixel` and `grid.locate_warp`, then reports each
combination's x error and E11 spread.
"""

import numpy as np
from scipy.ndimage import map_coordinates

from dictk.element import gauss_point_log_strains
from dictk.grid import elements, generate, locate_subpixel, locate_warp
from dictk.image import PixelCoordinate, read, stretch

FACTOR_X = 1.02

reference_image = read(path="astronaut0.png")
height, width = reference_image.shape
ys, xs = np.mgrid[0:height, 0:width].astype(np.float64)
generators = {
    "bilinear (`image.stretch`)": stretch(arr=reference_image, factor_x=FACTOR_X),
    "quintic spline": map_coordinates(
        reference_image.astype(np.float64), [ys, xs / FACTOR_X], order=5
    ),
}
points = generate(
    origin=PixelCoordinate(x=18, y=16),
    count_x=53,
    count_y=54,
    spacing_x=5,
    spacing_y=5,
)
element_indices = elements(count_x=53, count_y=54)
true_x = np.array([p.x * FACTOR_X for p in points])

print("| Current image | Tracker | Std $x$ error (px) | Std $E_{11}$ (µε) |")
print("|---|---|---|---|")
for generator, current_image in generators.items():
    for name, locate, extra in [
        ("`locate_subpixel`", locate_subpixel, dict(upsample_factor=100)),
        ("`locate_warp`", locate_warp, {}),
    ]:
        found = locate(
            reference_image=reference_image,
            current_image=current_image,
            reference_points=points,
            kernel_margin_width=13,
            kernel_margin_height=13,
            search_margin_width=25,
            search_margin_height=25,
            **extra,
        )
        error_x = np.array([f.x for f in found]) - true_x
        strains = np.array(
            [
                strain[0, 0]
                for element in element_indices
                for strain in gauss_point_log_strains(
                    reference_points=[points[i] for i in element],
                    current_points=[found[i] for i in element],
                )
            ]
        )
        print(
            f"| {generator} | {name} | {error_x.std():.4f} "
            f"| {strains.std() * 1e6:.0f} |"
        )

`kernel_warping_x_error_histogram.py`

Show the code

"""Distribution of `grid.locate_warp`'s x error across the full 53x54
grid on `image.stretch`'s bilinear image: the tracked x minus the true
x = 1.02 X, with the mean and one standard deviation marked.
"""

import matplotlib.pyplot as plt
import numpy as np

from dictk.grid import generate, locate_warp
from dictk.image import PixelCoordinate, read, stretch

FACTOR_X = 1.02

reference_image = read(path="astronaut0.png")
current_image = stretch(arr=reference_image, factor_x=FACTOR_X)
points = generate(
    origin=PixelCoordinate(x=18, y=16),
    count_x=53,
    count_y=54,
    spacing_x=5,
    spacing_y=5,
)
found = locate_warp(
    reference_image=reference_image,
    current_image=current_image,
    reference_points=points,
    kernel_margin_width=13,
    kernel_margin_height=13,
    search_margin_width=25,
    search_margin_height=25,
)

error_x = np.array([f.x for f in found]) - np.array([p.x * FACTOR_X for p in points])
mean = error_x.mean()
std = error_x.std()
print(
    f"{len(points)} points: mean x error {mean:+.5f} px, std {std:.4f} px, "
    f"range [{error_x.min():+.3f}, {error_x.max():+.3f}] px"
)

plt.rcParams.update({"font.family": "serif", "mathtext.fontset": "cm"})
fig, ax = plt.subplots(figsize=(7, 4), constrained_layout=True)
ax.hist(error_x, bins=60, color="gray", alpha=0.8, label=r"$x$ error at 2862 points")
ax.axvspan(
    mean - std,
    mean + std,
    color="royalblue",
    alpha=0.15,
    label=rf"$\pm$ 1 standard deviation, {std:.4f} px",
)
mean_text = f"{mean:+.4f}" if round(mean, 4) else "0.0000"
ax.axvline(mean, color="black", linewidth=1.5, label=f"mean, {mean_text} px")
ax.axvline(0, color="red", linestyle="--", linewidth=1.5, label="exact tracking, 0 px")
ax.set_xlabel(r"$x$ error (px)")
ax.set_ylabel("frequency")
ax.set_ylim(top=ax.get_ylim()[1] * 1.45)
legend = ax.legend(loc="upper right", framealpha=1.0)
legend.set_zorder(10)
fig.savefig("kernel_warping_x_error_histogram.png", dpi=300)

`kernel_warping_summary.py`

Show the code

"""Summary of Kernel Warping's progression as one HTML table:
`locate_subpixel`, VIC-2D, and `locate_warp` on `image.stretch`'s
bilinear image, then `locate_warp` on a quintic-spline image. Reports the
mean and standard deviation of each tracker's x error and of its E11.
"""

import csv

import numpy as np
from scipy.ndimage import map_coordinates

from dictk.element import gauss_point_log_strains
from dictk.grid import elements, generate, locate_subpixel, locate_warp
from dictk.image import PixelCoordinate, read, stretch

FACTOR_X = 1.02

reference_image = read(path="astronaut0.png")
height, width = reference_image.shape
ys, xs = np.mgrid[0:height, 0:width].astype(np.float64)
bilinear_image = stretch(arr=reference_image, factor_x=FACTOR_X)
quintic_image = map_coordinates(
    reference_image.astype(np.float64), [ys, xs / FACTOR_X], order=5
)
points = generate(
    origin=PixelCoordinate(x=18, y=16),
    count_x=53,
    count_y=54,
    spacing_x=5,
    spacing_y=5,
)
element_indices = elements(count_x=53, count_y=54)
true_x = np.array([p.x * FACTOR_X for p in points])


def track(tracker, image, **extra):
    return tracker(
        reference_image=reference_image,
        current_image=image,
        reference_points=points,
        kernel_margin_width=13,
        kernel_margin_height=13,
        search_margin_width=25,
        search_margin_height=25,
        **extra,
    )


def statistics(found):
    error_x = np.array([f.x for f in found]) - true_x
    strains = (
        np.array(
            [
                strain[0, 0]
                for element in element_indices
                for strain in gauss_point_log_strains(
                    reference_points=[points[i] for i in element],
                    current_points=[found[i] for i in element],
                )
            ]
        )
        * 1e6
    )
    return (
        error_x.mean(),
        error_x.std(),
        (error_x.min(), error_x.max()),
        strains.mean(),
        strains.std(),
        (strains.min(), strains.max()),
    )


with open("../verification/simple_stretch_vic_out.csv") as f:
    rows = [{k.strip(' "'): v for k, v in row.items()} for row in csv.DictReader(f)]
valid = [r for r in rows if float(r["sigma"]) != -1]
vic_error_x = np.array(
    [float(r["x"]) + float(r["u"]) - FACTOR_X * float(r["x"]) for r in valid]
)
vic_strain = np.array([float(r["exx"]) * 1e6 for r in valid])

columns = [
    (
        "<code>locate_subpixel</code>",
        statistics(track(locate_subpixel, bilinear_image, upsample_factor=100)),
        "",
    ),
    (
        "VIC-2D",
        (
            vic_error_x.mean(),
            vic_error_x.std(),
            (vic_error_x.min(), vic_error_x.max()),
            vic_strain.mean(),
            vic_strain.std(),
            (vic_strain.min(), vic_strain.max()),
        ),
        "",
    ),
    ("<code>locate_warp</code>", statistics(track(locate_warp, bilinear_image)), ""),
    (
        "<code>locate_warp</code>, quintic image",
        statistics(track(locate_warp, quintic_image)),
        ' class="unreachable"',
    ),
]


def range_format(bounds):
    decimals = 2 if max(abs(b) for b in bounds) >= 0.1 else 3
    low, high = bounds
    return f"[{low:+.{decimals}f}, {high:+.{decimals}f}]"


labels = [
    ("Mean $x$ error (px)", lambda v: f"{v:+.4f}" if round(v, 4) else "0.0000"),
    ("Std $x$ error (px)", lambda v: f"{v:.4f}"),
    ("$x$ error range (px)", range_format),
    ("Mean $E_{11}$ (microstrain)", lambda v: f"{v:.0f}"),
    ("Std $E_{11}$ (microstrain)", lambda v: f"{v:.0f}"),
    ("$E_{11}$ range (microstrain)", lambda v: f"[{round(v[0])}, {round(v[1])}]"),
]
figure_rows = {
    0: (
        "$x$ error Figure",
        {
            "VIC-2D": "./simple_stretch.html#fig-vic-x-error-histogram",
            "<code>locate_subpixel</code>": "./high_point_density.html#fig-hpd-x-error-histogram",
            "<code>locate_warp</code>": "#fig-kw-x-error-histogram",
        },
    ),
    3: (
        "$E_{11}$ Figure",
        {
            "VIC-2D": "./simple_stretch.html#fig-vic-exx-histogram",
            "<code>locate_subpixel</code>": "./high_point_density.html#fig-hpd-strain-histogram",
            "<code>locate_warp</code>": "#fig-kw-strain-histogram",
        },
    ),
}

print('<table class="progression">')
print("<thead><tr><th></th>", end="")
for name, _, css in columns:
    print(f"<th{css}>{name}</th>", end="")
print("</tr></thead>")
print("<tbody>")
for row, (label, fmt) in enumerate(labels):
    if row in figure_rows:
        figure_label, targets = figure_rows[row]
        print(f"<tr><td>{figure_label}</td>", end="")
        for name, _, css in columns:
            link = f'<a href="{targets[name]}">Figure</a>' if name in targets else ""
            print(f"<td{css}>{link}</td>", end="")
        print("</tr>")
    print(f"<tr><td>{label}</td>", end="")
    for _, values, css in columns:
        print(f"<td{css}>{fmt(values[row])}</td>", end="")
    print("</tr>")
print("</tbody>")
print("</table>")

Keyboard shortcuts

dictk User Guide