Coverage for src/dictk/image.py: 99%
172 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-08-25 21:06 +0000
1"""Image I/O, grayscale conversion, combination, and inspection utilities."""
3import base64
4import importlib.resources
5from pathlib import Path
6from typing import NamedTuple
8import imageio.v3 as iio
9import numpy as np
10from scipy.interpolate import RegularGridInterpolator
11from scipy.ndimage import zoom
14class PixelCoordinate(NamedTuple):
15 """A pixel location in an image's reference frame.
17 Attributes:
18 x: Horizontal pixel coordinate, increasing left to right.
19 y: Vertical pixel coordinate, increasing top to bottom.
20 """
22 x: int
23 y: int
26class SubpixelCoordinate(NamedTuple):
27 """A sub-pixel location in an image's reference frame.
29 Same reference frame as [`PixelCoordinate`](#PixelCoordinate), but
30 `x`/`y` are `float`, not `int` -- returned by
31 [`dictk.translation.locate_subpixel`](./translation.html#locate_subpixel)
32 and [`dictk.grid.locate_subpixel`](./grid.html#locate_subpixel),
33 which recover a fractional position via subpixel refinement.
34 `PixelCoordinate`'s own `int` fields would silently truncate exactly
35 the value those functions exist to preserve.
37 Attributes:
38 x: Horizontal pixel coordinate, increasing left to right.
39 y: Vertical pixel coordinate, increasing top to bottom.
40 """
42 x: float
43 y: float
46def subimage(
47 *, image: np.ndarray, origin: PixelCoordinate, width: int, height: int
48) -> np.ndarray:
49 """Extract a width x height crop of `image` with its top-left corner at `origin`.
51 `origin` is expressed in `image`'s own pixel reference frame (`image`'s
52 top-left corner is `(0, 0)`) and may lie partially or entirely outside
53 `image`'s bounds. Wherever the requested region doesn't overlap
54 `image`, the corresponding output pixels are zero (black) rather than
55 raising an error — so a subimage straddling an edge, or lying
56 completely outside `image`, is always a valid, well-shaped result.
58 Args:
59 image: Source 2D grayscale image array.
60 origin: Top-left corner of the region to extract, in `image`'s
61 pixel reference frame.
62 width: Width of the extracted region in pixels. Must be >= 1.
63 height: Height of the extracted region in pixels. Must be >= 1.
65 Returns:
66 A 2D array of shape (height, width), same dtype as `image`.
68 Raises:
69 ValueError: If `image` is not 2D, or width or height is less than 1.
70 """
71 if image.ndim != 2:
72 raise ValueError(f"image must be 2D, got shape {image.shape}")
73 if width < 1:
74 raise ValueError(f"width {width} must be >= 1")
75 if height < 1:
76 raise ValueError(f"height {height} must be >= 1")
78 image_height, image_width = image.shape
79 result = np.zeros((height, width), dtype=image.dtype)
81 src_x_start = max(origin.x, 0)
82 src_y_start = max(origin.y, 0)
83 src_x_end = min(origin.x + width, image_width)
84 src_y_end = min(origin.y + height, image_height)
86 if src_x_start >= src_x_end or src_y_start >= src_y_end:
87 return result # requested region doesn't overlap image at all
89 dst_x_start = src_x_start - origin.x
90 dst_y_start = src_y_start - origin.y
91 dst_x_end = dst_x_start + (src_x_end - src_x_start)
92 dst_y_end = dst_y_start + (src_y_end - src_y_start)
94 result[dst_y_start:dst_y_end, dst_x_start:dst_x_end] = image[
95 src_y_start:src_y_end, src_x_start:src_x_end
96 ]
97 return result
100def checkerboard(
101 *, width: int, height: int, count_x: int = 8, count_y: int = 8
102) -> np.ndarray:
103 """Generate a black-and-white checkerboard test image.
105 Same parameters and pixel values as the `dictk checkerboard` CLI
106 command, minus the file write — see the `dictk` package docstring for
107 why the API layer stops at the array. count_x and count_y are named for
108 the rectangles they divide the image into, not "squares": unequal
109 counts (or a width/height ratio that doesn't match count_x/count_y)
110 produce rectangular cells, not square ones.
112 Args:
113 width: Image width in pixels.
114 height: Image height in pixels.
115 count_x: Number of rectangles along the width. Must be >= 1.
116 count_y: Number of rectangles along the height. Must be >= 1.
118 Returns:
119 A 2D uint8 array of shape (height, width) with values 0 or 255.
121 Raises:
122 ValueError: If count_x or count_y is less than 1.
123 """
124 if count_x < 1:
125 raise ValueError(f"count_x {count_x} must be >= 1")
126 if count_y < 1:
127 raise ValueError(f"count_y {count_y} must be >= 1")
129 rows = (np.arange(height) * count_y // height) % 2
130 cols = (np.arange(width) * count_x // width) % 2
131 pattern = np.logical_xor(rows[:, None], cols[None, :])
132 return (pattern * 255).astype(np.uint8)
135def astronaut(*, width: int = 512, height: int = 512) -> np.ndarray:
136 """Load a bundled real-world grayscale reference image.
138 Same parameters and pixel values as the `dictk astronaut` CLI
139 command, minus the file write — see the `dictk` package docstring for
140 why the API layer stops at the array.
142 The source is a NASA portrait of astronaut Eileen Collins, from the
143 NASA Great Images database ("No known copyright restrictions, released
144 into the public domain."). It's bundled as a color image and converted
145 with `rgba_to_gray`; unlike `rosta` and `checkerboard`, it isn't
146 procedurally generated, so `width`/`height` other than the native
147 512x512 resize the source image (via `scipy.ndimage.zoom`) rather than
148 computing a fresh pattern at that size.
150 Args:
151 width: Image width in pixels. Must be >= 1.
152 height: Image height in pixels. Must be >= 1.
154 Returns:
155 A 2D uint8 array of shape (height, width).
157 Raises:
158 ValueError: If width or height is less than 1.
159 """
160 if width < 1:
161 raise ValueError(f"width {width} must be >= 1")
162 if height < 1:
163 raise ValueError(f"height {height} must be >= 1")
165 asset_path = importlib.resources.files("dictk") / "data" / "astronaut.png"
166 with importlib.resources.as_file(asset_path) as path:
167 color = read(path=path)
168 gray = rgba_to_gray(color)
170 native_height, native_width = gray.shape
171 if (width, height) == (native_width, native_height):
172 return gray
174 zoom_factors = (height / native_height, width / native_width)
175 resized = zoom(gray.astype(np.float64), zoom_factors, order=3)
176 return np.clip(resized, 0, 255).astype(np.uint8)
179def is_rgba(arr: np.ndarray) -> bool:
180 """Check whether an image array is in RGB or RGBA format.
182 Args:
183 arr: Input image array.
185 Returns:
186 True if the array is 3D with 3 or 4 channels, False otherwise.
187 """
188 return arr.ndim == 3 and arr.shape[2] in (3, 4)
191def rgba_to_gray(arr: np.ndarray) -> np.ndarray:
192 """Convert an RGB(A) image to grayscale by averaging the RGB channels.
194 Args:
195 arr: Input image array, either 2D (grayscale) or 3D (color).
197 Returns:
198 A 2D grayscale image array. If the input is already 2D, it is
199 returned unchanged.
201 Raises:
202 ValueError: If the array is neither 2D nor a 3-or-4-channel 3D array.
203 """
204 if arr.ndim == 2:
205 return arr
207 if is_rgba(arr):
208 return np.mean(arr[:, :, :3], axis=2).astype(arr.dtype)
210 raise ValueError(
211 "Input array must be either 2D (grayscale) or 3D with 3 or 4 channels (color)."
212 )
215def combine(*, a: np.ndarray, b: np.ndarray) -> np.ndarray:
216 """Combine two images by averaging their pixel values.
218 Args:
219 a: First image, 2D grayscale or 3D color.
220 b: Second image, same shape as `a` once converted to grayscale.
222 Returns:
223 A 2D uint8 array, normalized to the range [0, 255].
225 Raises:
226 ValueError: If the grayscale-converted images don't share a shape.
227 """
228 gray_a = rgba_to_gray(a)
229 gray_b = rgba_to_gray(b)
231 if gray_a.shape != gray_b.shape:
232 raise ValueError(
233 f"shape mismatch: a.shape={gray_a.shape}, b.shape={gray_b.shape}"
234 )
236 combined = gray_a.astype(np.float64) + gray_b.astype(np.float64)
237 return (combined / combined.max() * 255).astype(np.uint8)
240def brightness(*, arr: np.ndarray, factor: float) -> np.ndarray:
241 """Adjust image brightness by an additive shift, clipped to [0, 255].
243 Brightness *translates* the pixel-intensity histogram: every pixel is
244 shifted by the same amount, so dark areas lighten right along with
245 bright ones (unlike a multiplicative scale, where a black pixel would
246 stay exactly black). factor=1.0 leaves the image unchanged; factor=1.5
247 shifts every pixel by +127.5*0.5 = +63.75, factor=2.0 by +127.5.
249 Args:
250 arr: A 2D grayscale image array.
251 factor: Brightness factor. 1.0 is unchanged; > 1.0 brightens
252 (shifts toward white); < 1.0 darkens (shifts toward black).
254 Returns:
255 A 2D uint8 array, same shape as `arr`.
256 """
257 max_pixel_value = 255.0
258 offset = (factor - 1.0) * (max_pixel_value / 2.0)
259 shifted = arr.astype(np.float64) + offset
260 return np.clip(shifted, 0, max_pixel_value).astype(np.uint8)
263def contrast(*, arr: np.ndarray, factor: float) -> np.ndarray:
264 """Adjust image contrast by scaling around the mean, clipped to [0, 255].
266 Contrast *stretches* the pixel-intensity histogram outward from its own
267 mean, rather than shifting it: the mean stays roughly the same, while
268 values spread further from it. factor=1.0 leaves the image unchanged;
269 factor=0.0 collapses every pixel to the mean (a flat gray image);
270 factor > 1.0 pushes values further toward 0 and 255.
272 Args:
273 arr: A 2D grayscale image array.
274 factor: Contrast factor. 1.0 is unchanged; > 1.0 increases
275 contrast; between 0.0 and 1.0 decreases it.
277 Returns:
278 A 2D uint8 array, same shape as `arr`.
279 """
280 mean = arr.astype(np.float64).mean()
281 stretched = (arr.astype(np.float64) - mean) * factor + mean
282 return np.clip(stretched, 0, 255).astype(np.uint8)
285def _backward_map(
286 arr: np.ndarray, xs_source: np.ndarray, ys_source: np.ndarray
287) -> np.ndarray:
288 """Sample `arr` via bilinear interpolation at (xs_source, ys_source).
290 Shared by geometric-transform functions (`stretch`, `translate`, ...):
291 each computes where every output pixel's source coordinate falls under
292 its own inverse transform, then hands the resulting coordinate grids
293 here to do the actual sampling. Coordinates outside `arr`'s bounds are
294 filled with black (0).
296 Args:
297 arr: A 2D grayscale image array.
298 xs_source: Source x-coordinate for each output pixel, shape (height, width).
299 ys_source: Source y-coordinate for each output pixel, shape (height, width).
301 Returns:
302 A 2D uint8 array, same shape as `arr`.
303 """
304 height, width = arr.shape
305 interpolator = RegularGridInterpolator(
306 points=(np.arange(height), np.arange(width)),
307 values=arr.astype(np.float64),
308 method="linear",
309 bounds_error=False,
310 fill_value=0.0,
311 )
312 points = np.stack((ys_source.ravel(), xs_source.ravel()), axis=1)
313 deformed = interpolator(points).reshape(arr.shape)
314 return np.clip(deformed, 0, 255).astype(np.uint8)
317def stretch(
318 *, arr: np.ndarray, factor_x: float = 1.0, factor_y: float = 1.0
319) -> np.ndarray:
320 """Apply a uniaxial or biaxial stretch, pivoting on the image origin.
322 Mimics a continuum-mechanics stretch deformation gradient
323 diag(factor_x, factor_y), anchored at the image's top-left corner
324 (x=0, y=0): that corner stays fixed, and content grows (factor > 1.0)
325 or shrinks (factor < 1.0) away from it along each axis. Uses backward
326 mapping — for each output pixel, the inverse of the stretch locates
327 its source coordinate in `arr`, with bilinear interpolation for
328 non-integer source coordinates — so the result has no gaps, unlike
329 naively moving each source pixel forward. A factor < 1.0 shrinks
330 content toward the origin, leaving black (fill value 0) margins along
331 the far (bottom/right) edges.
333 Args:
334 arr: A 2D grayscale image array.
335 factor_x: Stretch factor along the x-axis. Must be > 0.
336 factor_y: Stretch factor along the y-axis. Must be > 0.
338 Returns:
339 A 2D uint8 array, same shape as `arr`.
341 Raises:
342 ValueError: If factor_x or factor_y is <= 0.
343 """
344 if factor_x <= 0:
345 raise ValueError(f"factor_x {factor_x} must be > 0")
346 if factor_y <= 0:
347 raise ValueError(f"factor_y {factor_y} must be > 0")
349 height, width = arr.shape
350 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
352 # Backward mapping, pivoting on the origin (0, 0): for each output
353 # pixel, the inverse of the stretch gives the coordinate to sample
354 # from in the original image.
355 xs_source = xs / factor_x
356 ys_source = ys / factor_y
358 return _backward_map(arr, xs_source, ys_source)
361def translate(*, arr: np.ndarray, dx: float = 0.0, dy: float = 0.0) -> np.ndarray:
362 """Apply a rigid-body translation: every pixel shifts by (dx, dy).
364 A pure displacement, with no change in shape or size — the simplest
365 transformation category. Uses the same backward-mapping approach as
366 `stretch`, so non-integer displacements are handled with bilinear
367 interpolation rather than rounding. Content shifted in from outside
368 the original bounds is filled with black (fill value 0).
370 Args:
371 arr: A 2D grayscale image array.
372 dx: Displacement in pixels along the x-axis; positive moves
373 content right.
374 dy: Displacement in pixels along the y-axis; positive moves
375 content down.
377 Returns:
378 A 2D uint8 array, same shape as `arr`.
379 """
380 height, width = arr.shape
381 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
383 # Backward mapping: the source of output pixel (x, y) is (x - dx, y - dy).
384 xs_source = xs - dx
385 ys_source = ys - dy
387 return _backward_map(arr, xs_source, ys_source)
390def rotate(*, arr: np.ndarray, angle: float) -> np.ndarray:
391 """Apply a rigid-body rotation, pivoting on the image origin.
393 Rotates content by `angle` degrees, positive counterclockwise,
394 pivoting on the image's top-left corner (0, 0) rather than its
395 center — consistent with `stretch` and `translate`'s pivot choice in
396 this codebase, but unlike a typical "object spins in place" rotation
397 example: most content swings away from that fixed corner, similar to
398 a door on a hinge. Uses the same backward-mapping approach as
399 `stretch` and `translate`, so non-integer source coordinates are
400 bilinearly interpolated, and any pixel with no corresponding source
401 coordinate within `arr`'s bounds is filled with black (fill value 0).
403 Args:
404 arr: A 2D grayscale image array.
405 angle: Rotation angle in degrees; positive is counterclockwise.
407 Returns:
408 A 2D uint8 array, same shape as `arr`.
409 """
410 height, width = arr.shape
411 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
413 theta = np.deg2rad(angle)
414 cos_theta = np.cos(theta)
415 sin_theta = np.sin(theta)
417 # Backward mapping: applying the inverse (-angle) rotation to each
418 # output coordinate gives the coordinate to sample from.
419 xs_source = cos_theta * xs + sin_theta * ys
420 ys_source = -sin_theta * xs + cos_theta * ys
422 return _backward_map(arr, xs_source, ys_source)
425def shear(*, arr: np.ndarray, shear_x: float = 0.0, shear_y: float = 0.0) -> np.ndarray:
426 """Apply a simple shear, pivoting on the image origin.
428 Mimics a continuum-mechanics shear deformation gradient
429 [[1, shear_x], [shear_y, 1]], anchored at the image's top-left corner
430 (x=0, y=0), consistent with `stretch`, `translate`, and `rotate`'s
431 pivot choice in this codebase: horizontal planes slide relative to
432 each other by an amount proportional to their y-coordinate
433 (shear_x), and/or vertical planes slide by an amount proportional to
434 their x-coordinate (shear_y). Uses the same backward-mapping approach
435 as the other transform functions, so non-integer source coordinates
436 are bilinearly interpolated, and any pixel with no corresponding
437 source coordinate within `arr`'s bounds is filled with black (fill
438 value 0).
440 Args:
441 arr: A 2D grayscale image array.
442 shear_x: Horizontal shear factor.
443 shear_y: Vertical shear factor.
445 Returns:
446 A 2D uint8 array, same shape as `arr`.
448 Raises:
449 ValueError: If shear_x * shear_y == 1, which makes the
450 deformation gradient singular (non-invertible).
451 """
452 determinant = 1.0 - shear_x * shear_y
453 if determinant == 0:
454 raise ValueError(
455 f"shear_x {shear_x} and shear_y {shear_y} produce a singular "
456 "deformation gradient (shear_x * shear_y == 1)"
457 )
459 height, width = arr.shape
460 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
462 # Backward mapping: invert the shear deformation gradient
463 # [[1, shear_x], [shear_y, 1]] to find each output pixel's source.
464 xs_source = (xs - shear_x * ys) / determinant
465 ys_source = (-shear_y * xs + ys) / determinant
467 return _backward_map(arr, xs_source, ys_source)
470def complex_deform(
471 *,
472 arr: np.ndarray,
473 factor_x: float = 1.0,
474 factor_y: float = 1.0,
475 angle: float = 0.0,
476) -> np.ndarray:
477 """Apply an anisotropic stretch composed with a rotation, in one pass.
479 Composes a stretch (deformation gradient diag(factor_x, factor_y))
480 with a rotation (`angle` degrees, counterclockwise): the combined
481 deformation gradient is F = R(angle) @ diag(factor_x, factor_y), i.e.
482 the stretch is applied first and the rotation second. Applying both
483 in a single backward-mapping pass (rather than calling `stretch` and
484 then `rotate` separately) avoids the extra blur of interpolating
485 twice. Represents realistic loading scenarios where materials
486 experience multiple simultaneous deformation modes — typically the
487 hardest case for correlation algorithms. Pivots on the image's
488 top-left corner (0, 0), consistent with `stretch`, `translate`,
489 `rotate`, and `shear`'s pivot choice in this codebase.
491 Args:
492 arr: A 2D grayscale image array.
493 factor_x: Stretch factor along the x-axis, applied before the
494 rotation. Must be > 0.
495 factor_y: Stretch factor along the y-axis, applied before the
496 rotation. Must be > 0.
497 angle: Rotation angle in degrees, applied after the stretch;
498 positive is counterclockwise.
500 Returns:
501 A 2D uint8 array, same shape as `arr`.
503 Raises:
504 ValueError: If factor_x or factor_y is <= 0.
505 """
506 if factor_x <= 0:
507 raise ValueError(f"factor_x {factor_x} must be > 0")
508 if factor_y <= 0:
509 raise ValueError(f"factor_y {factor_y} must be > 0")
511 height, width = arr.shape
512 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
514 theta = np.deg2rad(angle)
515 cos_theta = np.cos(theta)
516 sin_theta = np.sin(theta)
518 # Backward mapping for F = R(angle) @ diag(factor_x, factor_y):
519 # F_inv = diag(1/factor_x, 1/factor_y) @ R(-angle), applied to each
520 # output coordinate to find its source. Un-rotate first, then unscale.
521 xs_rotated = cos_theta * xs + sin_theta * ys
522 ys_rotated = -sin_theta * xs + cos_theta * ys
523 xs_source = xs_rotated / factor_x
524 ys_source = ys_rotated / factor_y
526 return _backward_map(arr, xs_source, ys_source)
529def crack_dislocation(*, arr: np.ndarray, offset: float = 8.0) -> np.ndarray:
530 """Apply a discontinuous vertical-crack displacement field.
532 Splits the image with a vertical crack at x = width / 2: the left
533 half shifts down by `offset` pixels and the right half shifts up by
534 `offset` pixels, producing a displacement field that jumps
535 discontinuously across the crack line — unlike every other transform
536 in this module, which varies smoothly. Standard DIC assumes smooth
537 displacements and cannot capture this jump; cases like this motivate
538 the Heaviside finite-element formulation. Uses the same backward
539 mapping as the other transform functions, just with a piecewise
540 (rather than single-matrix) displacement field, so non-integer
541 source coordinates are bilinearly interpolated, and any pixel with
542 no corresponding source coordinate within `arr`'s bounds is filled
543 with black (fill value 0).
545 Args:
546 arr: A 2D grayscale image array.
547 offset: Displacement in pixels; the left half (x < width / 2)
548 shifts down by this amount, the right half shifts up by it.
550 Returns:
551 A 2D uint8 array, same shape as `arr`.
552 """
553 height, width = arr.shape
554 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
556 xs_source = xs.astype(np.float64)
557 ys_source = ys.astype(np.float64)
559 left_half = xs < (width / 2)
560 ys_source[left_half] -= offset
561 ys_source[~left_half] += offset
563 return _backward_map(arr, xs_source, ys_source)
566def read(*, path: Path) -> np.ndarray:
567 """Read an image file into a NumPy array.
569 Args:
570 path: Path to the image file.
572 Returns:
573 The image as a NumPy array.
574 """
575 return iio.imread(path)
578def write_svg(*, arr: np.ndarray, path: Path) -> None:
579 """Write a NumPy array to an SVG file.
581 SVG is a vector format with no native pixel-grid concept, so the array
582 is PNG-encoded and embedded as a base64 data URI inside a minimal SVG
583 wrapper (the standard way to carry raster data in SVG) rather than
584 traced into vector shapes.
586 Args:
587 arr: The image array to save.
588 path: The output file path.
589 """
590 height, width = arr.shape[:2]
591 png_bytes = iio.imwrite("<bytes>", arr, extension=".png")
592 encoded = base64.b64encode(png_bytes).decode("ascii")
594 svg = (
595 '<?xml version="1.0" encoding="UTF-8"?>\n'
596 f'<svg xmlns="http://www.w3.org/2000/svg" '
597 f'xmlns:xlink="http://www.w3.org/1999/xlink" '
598 f'width="{width}" height="{height}" viewBox="0 0 {width} {height}">\n'
599 f' <image width="{width}" height="{height}" '
600 f'xlink:href="data:image/png;base64,{encoded}"/>\n'
601 "</svg>\n"
602 )
603 Path(path).write_text(svg, encoding="ascii")
606def write(*, arr: np.ndarray, path: Path) -> None:
607 """Write a NumPy array to an image file.
609 Dispatches on the file extension: `.svg` is handled by `write_svg`
610 (embedding a raster PNG in an SVG wrapper); every other extension
611 (`.tiff`, `.png`, `.jpg`, ...) is handled by imageio directly.
613 Args:
614 arr: The image array to save.
615 path: The output file path.
616 """
617 if Path(path).suffix.lower() == ".svg":
618 write_svg(arr=arr, path=path)
619 return
621 iio.imwrite(path, arr)
624def describe(arr: np.ndarray) -> str:
625 """Format a description of an image array's type, shape, and color format.
627 Args:
628 arr: The image array to describe.
630 Returns:
631 A multi-line description string.
632 """
633 lines = [
634 f"Type: {type(arr)}",
635 f"Shape: {arr.shape}",
636 f"Dtype: {arr.dtype}",
637 ]
639 match arr.shape:
640 case (_height, _width):
641 lines.append("The image is grayscale.")
642 case (_height, _width, 3):
643 lines.append("The image is color (RGB).")
644 case (_height, _width, 4):
645 lines.append("The image is color (RGBA, includes alpha channel).")
646 case _:
647 lines.append("The image has an unsupported format.")
649 return "\n".join(lines)