Coverage for src/dictk/imaging.py: 94%
153 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-24 22:28 +0000
1"""Image I/O, grayscale conversion, combination, and inspection utilities."""
3import base64
4import importlib.resources
5from pathlib import Path
7import imageio.v3 as iio
8import numpy as np
9from matplotlib import pyplot as plt
10from scipy.interpolate import RegularGridInterpolator
11from scipy.ndimage import zoom
14def checkerboard(
15 width: int, height: int, count_x: int = 8, count_y: int = 8
16) -> np.ndarray:
17 """Generate a black-and-white checkerboard test image.
19 Same parameters and pixel values as the `dictk checkerboard` CLI
20 command, minus the file write — see the `dictk` package docstring for
21 why the API layer stops at the array. count_x and count_y are named for
22 the rectangles they divide the image into, not "squares": unequal
23 counts (or a width/height ratio that doesn't match count_x/count_y)
24 produce rectangular cells, not square ones.
26 Args:
27 width: Image width in pixels.
28 height: Image height in pixels.
29 count_x: Number of rectangles along the width. Must be >= 1.
30 count_y: Number of rectangles along the height. Must be >= 1.
32 Returns:
33 A 2D uint8 array of shape (height, width) with values 0 or 255.
35 Raises:
36 ValueError: If count_x or count_y is less than 1.
37 """
38 if count_x < 1:
39 raise ValueError(f"count_x {count_x} must be >= 1")
40 if count_y < 1:
41 raise ValueError(f"count_y {count_y} must be >= 1")
43 rows = (np.arange(height) * count_y // height) % 2
44 cols = (np.arange(width) * count_x // width) % 2
45 pattern = np.logical_xor(rows[:, None], cols[None, :])
46 return (pattern * 255).astype(np.uint8)
49def astronaut(width: int = 512, height: int = 512) -> np.ndarray:
50 """Load a bundled real-world grayscale reference image.
52 Same parameters and pixel values as the `dictk astronaut` CLI
53 command, minus the file write — see the `dictk` package docstring for
54 why the API layer stops at the array.
56 The source is a NASA portrait of astronaut Eileen Collins, from the
57 NASA Great Images database ("No known copyright restrictions, released
58 into the public domain."). It's bundled as a color image and converted
59 with `rgba_to_gray`; unlike `rosta` and `checkerboard`, it isn't
60 procedurally generated, so `width`/`height` other than the native
61 512x512 resize the source image (via `scipy.ndimage.zoom`) rather than
62 computing a fresh pattern at that size.
64 Args:
65 width: Image width in pixels. Must be >= 1.
66 height: Image height in pixels. Must be >= 1.
68 Returns:
69 A 2D uint8 array of shape (height, width).
71 Raises:
72 ValueError: If width or height is less than 1.
73 """
74 if width < 1:
75 raise ValueError(f"width {width} must be >= 1")
76 if height < 1:
77 raise ValueError(f"height {height} must be >= 1")
79 asset_path = importlib.resources.files("dictk") / "data" / "astronaut.png"
80 with importlib.resources.as_file(asset_path) as path:
81 color = read_image(path)
82 gray = rgba_to_gray(color)
84 native_height, native_width = gray.shape
85 if (width, height) == (native_width, native_height):
86 return gray
88 zoom_factors = (height / native_height, width / native_width)
89 resized = zoom(gray.astype(np.float64), zoom_factors, order=3)
90 return np.clip(resized, 0, 255).astype(np.uint8)
93def is_rgba(arr: np.ndarray) -> bool:
94 """Check whether an image array is in RGB or RGBA format.
96 Args:
97 arr: Input image array.
99 Returns:
100 True if the array is 3D with 3 or 4 channels, False otherwise.
101 """
102 return arr.ndim == 3 and arr.shape[2] in (3, 4)
105def rgba_to_gray(arr: np.ndarray) -> np.ndarray:
106 """Convert an RGB(A) image to grayscale by averaging the RGB channels.
108 Args:
109 arr: Input image array, either 2D (grayscale) or 3D (color).
111 Returns:
112 A 2D grayscale image array. If the input is already 2D, it is
113 returned unchanged.
115 Raises:
116 ValueError: If the array is neither 2D nor a 3-or-4-channel 3D array.
117 """
118 if arr.ndim == 2:
119 return arr
121 if is_rgba(arr):
122 return np.mean(arr[:, :, :3], axis=2).astype(arr.dtype)
124 raise ValueError(
125 "Input array must be either 2D (grayscale) or 3D with 3 or 4 channels (color)."
126 )
129def combine_images(a: np.ndarray, b: np.ndarray) -> np.ndarray:
130 """Combine two images by averaging their pixel values.
132 Args:
133 a: First image, 2D grayscale or 3D color.
134 b: Second image, same shape as `a` once converted to grayscale.
136 Returns:
137 A 2D uint8 array, normalized to the range [0, 255].
139 Raises:
140 ValueError: If the grayscale-converted images don't share a shape.
141 """
142 gray_a = rgba_to_gray(a)
143 gray_b = rgba_to_gray(b)
145 if gray_a.shape != gray_b.shape:
146 raise ValueError(
147 f"shape mismatch: a.shape={gray_a.shape}, b.shape={gray_b.shape}"
148 )
150 combined = gray_a.astype(np.float64) + gray_b.astype(np.float64)
151 return (combined / combined.max() * 255).astype(np.uint8)
154def brightness(arr: np.ndarray, factor: float) -> np.ndarray:
155 """Adjust image brightness by an additive shift, clipped to [0, 255].
157 Brightness *translates* the pixel-intensity histogram: every pixel is
158 shifted by the same amount, so dark areas lighten right along with
159 bright ones (unlike a multiplicative scale, where a black pixel would
160 stay exactly black). factor=1.0 leaves the image unchanged; factor=1.5
161 shifts every pixel by +127.5*0.5 = +63.75, factor=2.0 by +127.5.
163 Args:
164 arr: A 2D grayscale image array.
165 factor: Brightness factor. 1.0 is unchanged; > 1.0 brightens
166 (shifts toward white); < 1.0 darkens (shifts toward black).
168 Returns:
169 A 2D uint8 array, same shape as `arr`.
170 """
171 max_pixel_value = 255.0
172 offset = (factor - 1.0) * (max_pixel_value / 2.0)
173 shifted = arr.astype(np.float64) + offset
174 return np.clip(shifted, 0, max_pixel_value).astype(np.uint8)
177def contrast(arr: np.ndarray, factor: float) -> np.ndarray:
178 """Adjust image contrast by scaling around the mean, clipped to [0, 255].
180 Contrast *stretches* the pixel-intensity histogram outward from its own
181 mean, rather than shifting it: the mean stays roughly the same, while
182 values spread further from it. factor=1.0 leaves the image unchanged;
183 factor=0.0 collapses every pixel to the mean (a flat gray image);
184 factor > 1.0 pushes values further toward 0 and 255.
186 Args:
187 arr: A 2D grayscale image array.
188 factor: Contrast factor. 1.0 is unchanged; > 1.0 increases
189 contrast; between 0.0 and 1.0 decreases it.
191 Returns:
192 A 2D uint8 array, same shape as `arr`.
193 """
194 mean = arr.astype(np.float64).mean()
195 stretched = (arr.astype(np.float64) - mean) * factor + mean
196 return np.clip(stretched, 0, 255).astype(np.uint8)
199def _backward_map(
200 arr: np.ndarray, xs_source: np.ndarray, ys_source: np.ndarray
201) -> np.ndarray:
202 """Sample `arr` via bilinear interpolation at (xs_source, ys_source).
204 Shared by geometric-transform functions (`stretch`, `translate`, ...):
205 each computes where every output pixel's source coordinate falls under
206 its own inverse transform, then hands the resulting coordinate grids
207 here to do the actual sampling. Coordinates outside `arr`'s bounds are
208 filled with black (0).
210 Args:
211 arr: A 2D grayscale image array.
212 xs_source: Source x-coordinate for each output pixel, shape (height, width).
213 ys_source: Source y-coordinate for each output pixel, shape (height, width).
215 Returns:
216 A 2D uint8 array, same shape as `arr`.
217 """
218 height, width = arr.shape
219 interpolator = RegularGridInterpolator(
220 points=(np.arange(height), np.arange(width)),
221 values=arr.astype(np.float64),
222 method="linear",
223 bounds_error=False,
224 fill_value=0.0,
225 )
226 points = np.stack((ys_source.ravel(), xs_source.ravel()), axis=1)
227 deformed = interpolator(points).reshape(arr.shape)
228 return np.clip(deformed, 0, 255).astype(np.uint8)
231def stretch(
232 arr: np.ndarray, factor_x: float = 1.0, factor_y: float = 1.0
233) -> np.ndarray:
234 """Apply a uniaxial or biaxial stretch, pivoting on the image origin.
236 Mimics a continuum-mechanics stretch deformation gradient
237 diag(factor_x, factor_y), anchored at the image's top-left corner
238 (x=0, y=0): that corner stays fixed, and content grows (factor > 1.0)
239 or shrinks (factor < 1.0) away from it along each axis. Uses backward
240 mapping — for each output pixel, the inverse of the stretch locates
241 its source coordinate in `arr`, with bilinear interpolation for
242 non-integer source coordinates — so the result has no gaps, unlike
243 naively moving each source pixel forward. A factor < 1.0 shrinks
244 content toward the origin, leaving black (fill value 0) margins along
245 the far (bottom/right) edges.
247 Args:
248 arr: A 2D grayscale image array.
249 factor_x: Stretch factor along the x-axis. Must be > 0.
250 factor_y: Stretch factor along the y-axis. Must be > 0.
252 Returns:
253 A 2D uint8 array, same shape as `arr`.
255 Raises:
256 ValueError: If factor_x or factor_y is <= 0.
257 """
258 if factor_x <= 0:
259 raise ValueError(f"factor_x {factor_x} must be > 0")
260 if factor_y <= 0:
261 raise ValueError(f"factor_y {factor_y} must be > 0")
263 height, width = arr.shape
264 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
266 # Backward mapping, pivoting on the origin (0, 0): for each output
267 # pixel, the inverse of the stretch gives the coordinate to sample
268 # from in the original image.
269 xs_source = xs / factor_x
270 ys_source = ys / factor_y
272 return _backward_map(arr, xs_source, ys_source)
275def translate(arr: np.ndarray, dx: float = 0.0, dy: float = 0.0) -> np.ndarray:
276 """Apply a rigid-body translation: every pixel shifts by (dx, dy).
278 A pure displacement, with no change in shape or size — the simplest
279 transformation category. Uses the same backward-mapping approach as
280 `stretch`, so non-integer displacements are handled with bilinear
281 interpolation rather than rounding. Content shifted in from outside
282 the original bounds is filled with black (fill value 0).
284 Args:
285 arr: A 2D grayscale image array.
286 dx: Displacement in pixels along the x-axis; positive moves
287 content right.
288 dy: Displacement in pixels along the y-axis; positive moves
289 content down.
291 Returns:
292 A 2D uint8 array, same shape as `arr`.
293 """
294 height, width = arr.shape
295 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
297 # Backward mapping: the source of output pixel (x, y) is (x - dx, y - dy).
298 xs_source = xs - dx
299 ys_source = ys - dy
301 return _backward_map(arr, xs_source, ys_source)
304def rotate(arr: np.ndarray, angle: float) -> np.ndarray:
305 """Apply a rigid-body rotation, pivoting on the image origin.
307 Rotates content by `angle` degrees, positive counterclockwise,
308 pivoting on the image's top-left corner (0, 0) rather than its
309 center — consistent with `stretch` and `translate`'s pivot choice in
310 this codebase, but unlike a typical "object spins in place" rotation
311 example: most content swings away from that fixed corner, similar to
312 a door on a hinge. Uses the same backward-mapping approach as
313 `stretch` and `translate`, so non-integer source coordinates are
314 bilinearly interpolated, and any pixel with no corresponding source
315 coordinate within `arr`'s bounds is filled with black (fill value 0).
317 Args:
318 arr: A 2D grayscale image array.
319 angle: Rotation angle in degrees; positive is counterclockwise.
321 Returns:
322 A 2D uint8 array, same shape as `arr`.
323 """
324 height, width = arr.shape
325 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
327 theta = np.deg2rad(angle)
328 cos_theta = np.cos(theta)
329 sin_theta = np.sin(theta)
331 # Backward mapping: applying the inverse (-angle) rotation to each
332 # output coordinate gives the coordinate to sample from.
333 xs_source = cos_theta * xs + sin_theta * ys
334 ys_source = -sin_theta * xs + cos_theta * ys
336 return _backward_map(arr, xs_source, ys_source)
339def shear(arr: np.ndarray, shear_x: float = 0.0, shear_y: float = 0.0) -> np.ndarray:
340 """Apply a simple shear, pivoting on the image origin.
342 Mimics a continuum-mechanics shear deformation gradient
343 [[1, shear_x], [shear_y, 1]], anchored at the image's top-left corner
344 (x=0, y=0), consistent with `stretch`, `translate`, and `rotate`'s
345 pivot choice in this codebase: horizontal planes slide relative to
346 each other by an amount proportional to their y-coordinate
347 (shear_x), and/or vertical planes slide by an amount proportional to
348 their x-coordinate (shear_y). Uses the same backward-mapping approach
349 as the other transform functions, so non-integer source coordinates
350 are bilinearly interpolated, and any pixel with no corresponding
351 source coordinate within `arr`'s bounds is filled with black (fill
352 value 0).
354 Args:
355 arr: A 2D grayscale image array.
356 shear_x: Horizontal shear factor.
357 shear_y: Vertical shear factor.
359 Returns:
360 A 2D uint8 array, same shape as `arr`.
362 Raises:
363 ValueError: If shear_x * shear_y == 1, which makes the
364 deformation gradient singular (non-invertible).
365 """
366 determinant = 1.0 - shear_x * shear_y
367 if determinant == 0:
368 raise ValueError(
369 f"shear_x {shear_x} and shear_y {shear_y} produce a singular "
370 "deformation gradient (shear_x * shear_y == 1)"
371 )
373 height, width = arr.shape
374 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
376 # Backward mapping: invert the shear deformation gradient
377 # [[1, shear_x], [shear_y, 1]] to find each output pixel's source.
378 xs_source = (xs - shear_x * ys) / determinant
379 ys_source = (-shear_y * xs + ys) / determinant
381 return _backward_map(arr, xs_source, ys_source)
384def complex_deform(
385 arr: np.ndarray, factor_x: float = 1.0, factor_y: float = 1.0, angle: float = 0.0
386) -> np.ndarray:
387 """Apply an anisotropic stretch composed with a rotation, in one pass.
389 Composes a stretch (deformation gradient diag(factor_x, factor_y))
390 with a rotation (`angle` degrees, counterclockwise): the combined
391 deformation gradient is F = R(angle) @ diag(factor_x, factor_y), i.e.
392 the stretch is applied first and the rotation second. Applying both
393 in a single backward-mapping pass (rather than calling `stretch` and
394 then `rotate` separately) avoids the extra blur of interpolating
395 twice. Represents realistic loading scenarios where materials
396 experience multiple simultaneous deformation modes — typically the
397 hardest case for correlation algorithms. Pivots on the image's
398 top-left corner (0, 0), consistent with `stretch`, `translate`,
399 `rotate`, and `shear`'s pivot choice in this codebase.
401 Args:
402 arr: A 2D grayscale image array.
403 factor_x: Stretch factor along the x-axis, applied before the
404 rotation. Must be > 0.
405 factor_y: Stretch factor along the y-axis, applied before the
406 rotation. Must be > 0.
407 angle: Rotation angle in degrees, applied after the stretch;
408 positive is counterclockwise.
410 Returns:
411 A 2D uint8 array, same shape as `arr`.
413 Raises:
414 ValueError: If factor_x or factor_y is <= 0.
415 """
416 if factor_x <= 0:
417 raise ValueError(f"factor_x {factor_x} must be > 0")
418 if factor_y <= 0:
419 raise ValueError(f"factor_y {factor_y} must be > 0")
421 height, width = arr.shape
422 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
424 theta = np.deg2rad(angle)
425 cos_theta = np.cos(theta)
426 sin_theta = np.sin(theta)
428 # Backward mapping for F = R(angle) @ diag(factor_x, factor_y):
429 # F_inv = diag(1/factor_x, 1/factor_y) @ R(-angle), applied to each
430 # output coordinate to find its source. Un-rotate first, then unscale.
431 xs_rotated = cos_theta * xs + sin_theta * ys
432 ys_rotated = -sin_theta * xs + cos_theta * ys
433 xs_source = xs_rotated / factor_x
434 ys_source = ys_rotated / factor_y
436 return _backward_map(arr, xs_source, ys_source)
439def crack_dislocation(arr: np.ndarray, offset: float = 8.0) -> np.ndarray:
440 """Apply a discontinuous vertical-crack displacement field.
442 Splits the image with a vertical crack at x = width / 2: the left
443 half shifts down by `offset` pixels and the right half shifts up by
444 `offset` pixels, producing a displacement field that jumps
445 discontinuously across the crack line — unlike every other transform
446 in this module, which varies smoothly. Standard DIC assumes smooth
447 displacements and cannot capture this jump; cases like this motivate
448 the Heaviside finite-element formulation. Uses the same backward
449 mapping as the other transform functions, just with a piecewise
450 (rather than single-matrix) displacement field, so non-integer
451 source coordinates are bilinearly interpolated, and any pixel with
452 no corresponding source coordinate within `arr`'s bounds is filled
453 with black (fill value 0).
455 Args:
456 arr: A 2D grayscale image array.
457 offset: Displacement in pixels; the left half (x < width / 2)
458 shifts down by this amount, the right half shifts up by it.
460 Returns:
461 A 2D uint8 array, same shape as `arr`.
462 """
463 height, width = arr.shape
464 xs, ys = np.meshgrid(np.arange(width), np.arange(height))
466 xs_source = xs.astype(np.float64)
467 ys_source = ys.astype(np.float64)
469 left_half = xs < (width / 2)
470 ys_source[left_half] -= offset
471 ys_source[~left_half] += offset
473 return _backward_map(arr, xs_source, ys_source)
476def read_image(path: Path) -> np.ndarray:
477 """Read an image file into a NumPy array.
479 Args:
480 path: Path to the image file.
482 Returns:
483 The image as a NumPy array.
484 """
485 return iio.imread(path)
488def write_svg(arr: np.ndarray, path: Path) -> None:
489 """Write a NumPy array to an SVG file.
491 SVG is a vector format with no native pixel-grid concept, so the array
492 is PNG-encoded and embedded as a base64 data URI inside a minimal SVG
493 wrapper (the standard way to carry raster data in SVG) rather than
494 traced into vector shapes.
496 Args:
497 arr: The image array to save.
498 path: The output file path.
499 """
500 height, width = arr.shape[:2]
501 png_bytes = iio.imwrite("<bytes>", arr, extension=".png")
502 encoded = base64.b64encode(png_bytes).decode("ascii")
504 svg = (
505 '<?xml version="1.0" encoding="UTF-8"?>\n'
506 f'<svg xmlns="http://www.w3.org/2000/svg" '
507 f'xmlns:xlink="http://www.w3.org/1999/xlink" '
508 f'width="{width}" height="{height}" viewBox="0 0 {width} {height}">\n'
509 f' <image width="{width}" height="{height}" '
510 f'xlink:href="data:image/png;base64,{encoded}"/>\n'
511 "</svg>\n"
512 )
513 Path(path).write_text(svg, encoding="ascii")
516def write_image(arr: np.ndarray, path: Path) -> None:
517 """Write a NumPy array to an image file.
519 Dispatches on the file extension: `.svg` is handled by `write_svg`
520 (embedding a raster PNG in an SVG wrapper); every other extension
521 (`.tiff`, `.png`, `.jpg`, ...) is handled by imageio directly.
523 Args:
524 arr: The image array to save.
525 path: The output file path.
526 """
527 if Path(path).suffix.lower() == ".svg":
528 write_svg(arr, path)
529 return
531 iio.imwrite(path, arr)
534def describe_image(arr: np.ndarray) -> str:
535 """Format a description of an image array's type, shape, and color format.
537 Args:
538 arr: The image array to describe.
540 Returns:
541 A multi-line description string.
542 """
543 lines = [
544 f"Type: {type(arr)}",
545 f"Shape: {arr.shape}",
546 f"Dtype: {arr.dtype}",
547 ]
549 match arr.shape:
550 case (_height, _width):
551 lines.append("The image is grayscale.")
552 case (_height, _width, 3):
553 lines.append("The image is color (RGB).")
554 case (_height, _width, 4):
555 lines.append("The image is color (RGBA, includes alpha channel).")
556 case _:
557 lines.append("The image has an unsupported format.")
559 return "\n".join(lines)
562def save_histogram(arr: np.ndarray, path: Path, dpi: int = 300) -> None:
563 """Save a histogram of pixel intensities [0, 255] for a grayscale image.
565 Args:
566 arr: The 2D grayscale image array, expected type uint8, range [0, 255].
567 path: The output file path for the histogram image.
568 dpi: Resolution of the saved figure.
569 """
570 plt.figure()
571 plt.hist(arr.ravel(), bins=256, range=(0, 255), color="black", alpha=0.7)
572 plt.title(f"pixel histogram intensity\n{path}", fontsize=8)
573 plt.xlabel("pixel intensity (0-255)")
574 plt.ylabel("frequency")
575 plt.savefig(path, dpi=dpi)
576 plt.close()