Coverage for src/dictk/image.py: 97%

1"""Image I/O, grayscale conversion, combination, and inspection utilities.""" 

2 

3import base64 

4import importlib.resources 

5from collections.abc import Sequence 

6from pathlib import Path 

7from typing import NamedTuple 

8 

9import imageio.v3 as iio 

10import numpy as np 

11from matplotlib import patches 

12from matplotlib import pyplot as plt 

13from scipy.interpolate import RegularGridInterpolator 

14from scipy.ndimage import zoom 

15 

16 

17class PixelCoordinate(NamedTuple): 

18 """A pixel location in an image's reference frame. 

19 

20 Attributes: 

21 x: Horizontal pixel coordinate, increasing left to right. 

22 y: Vertical pixel coordinate, increasing top to bottom. 

23 """ 

24 

25 x: int 

26 y: int 

27 

28 

29# Shared figure scale (pixels of image data per inch) for 

30# subimage_bounds_plot() and subimage_plot(), so each saved figure's size 

31# is proportional to its actual pixel content rather than a fixed default 

32# figure size -- letting the two be visually compared for relative size 

33# instead of both rendering at roughly the same size regardless of how 

34# many pixels each actually covers. 

35_FIGURE_PIXELS_PER_INCH = 100 

36 

37 

38def subimage( 

39 *, image: np.ndarray, origin: PixelCoordinate, width: int, height: int 

40) -> np.ndarray: 

41 """Extract a width x height crop of `image` with its top-left corner at `origin`. 

42 

43 `origin` is expressed in `image`'s own pixel reference frame (`image`'s 

44 top-left corner is `(0, 0)`) and may lie partially or entirely outside 

45 `image`'s bounds. Wherever the requested region doesn't overlap 

46 `image`, the corresponding output pixels are zero (black) rather than 

47 raising an error — so a subimage straddling an edge, or lying 

48 completely outside `image`, is always a valid, well-shaped result. 

49 

50 Args: 

51 image: Source 2D grayscale image array. 

52 origin: Top-left corner of the region to extract, in `image`'s 

53 pixel reference frame. 

54 width: Width of the extracted region in pixels. Must be >= 1. 

55 height: Height of the extracted region in pixels. Must be >= 1. 

56 

57 Returns: 

58 A 2D array of shape (height, width), same dtype as `image`. 

59 

60 Raises: 

61 ValueError: If `image` is not 2D, or width or height is less than 1. 

62 """ 

63 if image.ndim != 2: 

64 raise ValueError(f"image must be 2D, got shape {image.shape}") 

65 if width < 1: 

66 raise ValueError(f"width {width} must be >= 1") 

67 if height < 1: 

68 raise ValueError(f"height {height} must be >= 1") 

69 

70 image_height, image_width = image.shape 

71 result = np.zeros((height, width), dtype=image.dtype) 

72 

73 src_x_start = max(origin.x, 0) 

74 src_y_start = max(origin.y, 0) 

75 src_x_end = min(origin.x + width, image_width) 

76 src_y_end = min(origin.y + height, image_height) 

77 

78 if src_x_start >= src_x_end or src_y_start >= src_y_end: 

79 return result # requested region doesn't overlap image at all 

80 

81 dst_x_start = src_x_start - origin.x 

82 dst_y_start = src_y_start - origin.y 

83 dst_x_end = dst_x_start + (src_x_end - src_x_start) 

84 dst_y_end = dst_y_start + (src_y_end - src_y_start) 

85 

86 result[dst_y_start:dst_y_end, dst_x_start:dst_x_end] = image[ 

87 src_y_start:src_y_end, src_x_start:src_x_end 

88 ] 

89 return result 

90 

91 

92def subimage_bounds_plot( 

93 *, 

94 image: np.ndarray, 

95 origin: PixelCoordinate, 

96 width: int, 

97 height: int, 

98 path: Path, 

99 dpi: int = 300, 

100) -> None: 

101 """Save a figure overlaying a source image's bounds and a subimage region. 

102 

103 Draws `image`'s own bounds in blue and the requested `origin`/`width`/ 

104 `height` region in red, on top of `image` itself — including cases 

105 where the red region extends partially or completely outside the blue 

106 one, to visualize what `subimage()` would crop from. An `'o'` marker 

107 is also drawn at each rectangle's own origin: blue at `(0, 0)` for 

108 `image`, red at `origin` for the subimage — both in `image`'s pixel 

109 reference frame. 

110 

111 Args: 

112 image: Source 2D grayscale image array. 

113 origin: Top-left corner of the region, in `image`'s pixel 

114 reference frame; see `subimage()`. 

115 width: Width of the region in pixels. Must be >= 1. 

116 height: Height of the region in pixels. Must be >= 1. 

117 path: Output file path for the figure; format is inferred from 

118 the extension by matplotlib's savefig (e.g. .png), not 

119 dictk's own write/write_svg. 

120 dpi: Resolution of the saved figure. 

121 

122 Raises: 

123 ValueError: If width or height is less than 1. 

124 """ 

125 if width < 1: 

126 raise ValueError(f"width {width} must be >= 1") 

127 if height < 1: 

128 raise ValueError(f"height {height} must be >= 1") 

129 

130 image_height, image_width = image.shape 

131 

132 margin = max(width, height, image_width, image_height) * 0.05 

133 x_min = min(0, origin.x) - margin 

134 x_max = max(image_width, origin.x + width) + margin 

135 y_min = min(0, origin.y) - margin 

136 y_max = max(image_height, origin.y + height) + margin 

137 

138 figsize = ( 

139 (x_max - x_min) / _FIGURE_PIXELS_PER_INCH, 

140 (y_max - y_min) / _FIGURE_PIXELS_PER_INCH, 

141 ) 

142 fig, ax = plt.subplots(figsize=figsize) 

143 ax.imshow( 

144 image, cmap="gray", origin="upper", extent=(0, image_width, image_height, 0) 

145 ) 

146 

147 ax.add_patch( 

148 patches.Rectangle( 

149 (0, 0), 

150 image_width, 

151 image_height, 

152 edgecolor="blue", 

153 facecolor="none", 

154 linewidth=1.5, 

155 label="source image", 

156 ) 

157 ) 

158 ax.add_patch( 

159 patches.Rectangle( 

160 (origin.x, origin.y), 

161 width, 

162 height, 

163 edgecolor="red", 

164 facecolor="none", 

165 linewidth=1.5, 

166 label="subimage", 

167 ) 

168 ) 

169 ax.plot(0, 0, marker="o", color="blue", markersize=6) 

170 ax.plot(origin.x, origin.y, marker="o", color="red", markersize=6) 

171 

172 ax.set_xlim(x_min, x_max) 

173 ax.set_ylim(y_max, y_min) # inverted: image y increases downward 

174 

175 ax.set_xlabel("x (pixels)") 

176 ax.set_ylabel("y (pixels)") 

177 # Legend placed fully outside the axes (not just an "upper right"-style 

178 # corner) so it never overlaps the source or subimage rectangles, 

179 # regardless of how much of the frame they fill; bbox_inches="tight" 

180 # on save keeps it from being clipped off the saved figure. 

181 ax.legend(loc="center left", bbox_to_anchor=(1.02, 0.5), fontsize=8) 

182 ax.set_title(f"subimage bounds: origin=({origin.x}, {origin.y}), {width}x{height}") 

183 

184 plt.savefig(path, dpi=dpi, bbox_inches="tight") 

185 plt.close(fig) 

186 

187 

188def subimage_plot( 

189 *, 

190 image: np.ndarray, 

191 origin: PixelCoordinate, 

192 width: int, 

193 height: int, 

194 path: Path, 

195 dpi: int = 300, 

196) -> None: 

197 """Save a figure of the subimage itself, in its own local reference frame. 

198 

199 Extracts the `width` x `height` region at `origin` (via `subimage()`) 

200 and plots just that result, labeled with its own local pixel 

201 coordinates — `(0, 0)` at its own top-left corner — rather than 

202 `image`'s coordinates. An `'o'` marker is drawn at that local origin 

203 `(0, 0)`, matching the red origin marker `subimage_bounds_plot()` 

204 draws at the same subimage in `image`'s reference frame. Unlike 

205 `subimage_bounds_plot()`, which shows where the region falls relative 

206 to `image`, this shows what the extracted result actually looks like, 

207 including any zero (black) padding from parts of the region that fell 

208 outside `image`. 

209 

210 Args: 

211 image: Source 2D grayscale image array. 

212 origin: Top-left corner of the region, in `image`'s pixel 

213 reference frame; see `subimage()`. 

214 width: Width of the region in pixels. Must be >= 1. 

215 height: Height of the region in pixels. Must be >= 1. 

216 path: Output file path for the figure; format is inferred from 

217 the extension by matplotlib's savefig (e.g. .png), not 

218 dictk's own write/write_svg. 

219 dpi: Resolution of the saved figure. 

220 

221 Raises: 

222 ValueError: If width or height is less than 1. 

223 """ 

224 region = subimage(image=image, origin=origin, width=width, height=height) 

225 

226 # Same margin approach as subimage_bounds_plot(), so the red border 

227 # gets the same small breathing room from the figure edge instead of 

228 # sitting flush against it. 

229 margin = max(width, height) * 0.05 

230 x_min, x_max = -margin, width + margin 

231 y_min, y_max = -margin, height + margin 

232 

233 figsize = ( 

234 (x_max - x_min) / _FIGURE_PIXELS_PER_INCH, 

235 (y_max - y_min) / _FIGURE_PIXELS_PER_INCH, 

236 ) 

237 fig, ax = plt.subplots(figsize=figsize) 

238 ax.imshow(region, cmap="gray", origin="upper", extent=(0, width, height, 0)) 

239 ax.add_patch( 

240 patches.Rectangle( 

241 (0, 0), 

242 width, 

243 height, 

244 edgecolor="red", 

245 facecolor="none", 

246 linewidth=1.5, 

247 ) 

248 ) 

249 ax.plot(0, 0, marker="o", color="red", markersize=6) 

250 

251 ax.set_xlim(x_min, x_max) 

252 ax.set_ylim(y_max, y_min) # inverted: image y increases downward 

253 

254 ax.set_xlabel("x (pixels)") 

255 ax.set_ylabel("y (pixels)") 

256 

257 plt.savefig(path, dpi=dpi, bbox_inches="tight") 

258 plt.close(fig) 

259 

260 

261def subimage_comparison_plot( 

262 *, 

263 image: np.ndarray, 

264 origin: PixelCoordinate, 

265 width: int, 

266 height: int, 

267 path: Path, 

268 dpi: int = 300, 

269) -> None: 

270 """Save a side-by-side comparison of a subimage's placement and its extraction. 

271 

272 The left panel matches `subimage_bounds_plot()`: `image`'s bounds in 

273 blue, the requested region in red, with an `'o'` marker at each 

274 rectangle's own origin. The right panel shows the extracted subimage 

275 (via `subimage()`) with a red border and origin marker, in its own 

276 local reference frame — but drawn using the *same* axis limits as the 

277 left panel, rather than being cropped or zoomed to the subimage's own 

278 size the way `subimage_plot()` is. That shared scale is what makes 

279 the two red boxes render at identical size by construction — same 

280 data units per pixel in both panels — rather than approximating it by 

281 sizing each panel's figure independently around its own content. 

282 

283 Args: 

284 image: Source 2D grayscale image array. 

285 origin: Top-left corner of the region, in `image`'s pixel 

286 reference frame; see `subimage()`. 

287 width: Width of the region in pixels. Must be >= 1. 

288 height: Height of the region in pixels. Must be >= 1. 

289 path: Output file path for the figure; format is inferred from 

290 the extension by matplotlib's savefig (e.g. .png), not 

291 dictk's own write/write_svg. 

292 dpi: Resolution of the saved figure. 

293 

294 Raises: 

295 ValueError: If width or height is less than 1. 

296 """ 

297 if width < 1: 

298 raise ValueError(f"width {width} must be >= 1") 

299 if height < 1: 

300 raise ValueError(f"height {height} must be >= 1") 

301 

302 region = subimage(image=image, origin=origin, width=width, height=height) 

303 image_height, image_width = image.shape 

304 

305 margin = max(width, height, image_width, image_height) * 0.05 

306 x_min = min(0, origin.x) - margin 

307 x_max = max(image_width, origin.x + width) + margin 

308 y_min = min(0, origin.y) - margin 

309 y_max = max(image_height, origin.y + height) + margin 

310 

311 panel_width = (x_max - x_min) / _FIGURE_PIXELS_PER_INCH 

312 panel_height = (y_max - y_min) / _FIGURE_PIXELS_PER_INCH 

313 fig, (ax_left, ax_right) = plt.subplots( 

314 1, 2, figsize=(panel_width * 2, panel_height), constrained_layout=True 

315 ) 

316 

317 ax_left.imshow( 

318 image, cmap="gray", origin="upper", extent=(0, image_width, image_height, 0) 

319 ) 

320 ax_left.add_patch( 

321 patches.Rectangle( 

322 (0, 0), 

323 image_width, 

324 image_height, 

325 edgecolor="blue", 

326 facecolor="none", 

327 linewidth=1.5, 

328 ) 

329 ) 

330 ax_left.add_patch( 

331 patches.Rectangle( 

332 (origin.x, origin.y), 

333 width, 

334 height, 

335 edgecolor="red", 

336 facecolor="none", 

337 linewidth=1.5, 

338 ) 

339 ) 

340 ax_left.plot(0, 0, marker="o", color="blue", markersize=6) 

341 ax_left.plot(origin.x, origin.y, marker="o", color="red", markersize=6) 

342 ax_left.set_xlim(x_min, x_max) 

343 ax_left.set_ylim(y_max, y_min) # inverted: image y increases downward 

344 ax_left.set_xlabel("x (pixels)") 

345 ax_left.set_ylabel("y (pixels)") 

346 ax_left.set_title("source image + subimage") 

347 

348 ax_right.imshow(region, cmap="gray", origin="upper", extent=(0, width, height, 0)) 

349 ax_right.add_patch( 

350 patches.Rectangle( 

351 (0, 0), 

352 width, 

353 height, 

354 edgecolor="red", 

355 facecolor="none", 

356 linewidth=1.5, 

357 ) 

358 ) 

359 ax_right.plot(0, 0, marker="o", color="red", markersize=6) 

360 # Same limits as ax_left, not the subimage's own tight size -- this is 

361 # the whole point: identical data-units-per-pixel in both panels. 

362 ax_right.set_xlim(x_min, x_max) 

363 ax_right.set_ylim(y_max, y_min) 

364 ax_right.set_xlabel("x (pixels)") 

365 ax_right.set_ylabel("y (pixels)") 

366 ax_right.set_title("subimage") 

367 

368 plt.savefig(path, dpi=dpi, bbox_inches="tight") 

369 plt.close(fig) 

370 

371 

372class ArrowAnnotation(NamedTuple): 

373 """A labeled, colored arrow to overlay on an image, from `tail` to `head`. 

374 

375 Attributes: 

376 tail: Arrow's starting point, in the image's pixel reference frame. 

377 head: Arrow's ending point, in the image's pixel reference frame. 

378 color: Matplotlib color name for the arrow and its legend entry. 

379 label: Legend label for the arrow. 

380 """ 

381 

382 tail: PixelCoordinate 

383 head: PixelCoordinate 

384 color: str 

385 label: str 

386 

387 

388def point_plot( 

389 *, image: np.ndarray, arrows: Sequence[ArrowAnnotation], path: Path, dpi: int = 300 

390) -> None: 

391 """Save a figure overlaying one or more labeled arrows on `image`. 

392 

393 Each `ArrowAnnotation` draws a straight arrow from `tail` to `head`, in 

394 `image`'s own pixel reference frame — e.g. from the origin `(0, 0)` to a 

395 point of interest, or between two points to show a displacement. 

396 

397 Args: 

398 image: Source 2D grayscale image array. 

399 arrows: One or more arrows to overlay, each with its own color and 

400 legend label. 

401 path: Output file path for the figure; format is inferred from the 

402 extension by matplotlib's savefig (e.g. .png), not dictk's own 

403 write/write_svg. 

404 dpi: Resolution of the saved figure. 

405 

406 Raises: 

407 ValueError: If `arrows` is empty. 

408 """ 

409 if not arrows: 

410 raise ValueError("arrows must not be empty") 

411 

412 image_height, image_width = image.shape 

413 

414 endpoints_x = [pt.x for arrow in arrows for pt in (arrow.tail, arrow.head)] 

415 endpoints_y = [pt.y for arrow in arrows for pt in (arrow.tail, arrow.head)] 

416 margin = max(image_width, image_height) * 0.05 

417 x_min = min(0, *endpoints_x) - margin 

418 x_max = max(image_width, *endpoints_x) + margin 

419 y_min = min(0, *endpoints_y) - margin 

420 y_max = max(image_height, *endpoints_y) + margin 

421 

422 figsize = ( 

423 (x_max - x_min) / _FIGURE_PIXELS_PER_INCH, 

424 (y_max - y_min) / _FIGURE_PIXELS_PER_INCH, 

425 ) 

426 fig, ax = plt.subplots(figsize=figsize) 

427 ax.imshow( 

428 image, cmap="gray", origin="upper", extent=(0, image_width, image_height, 0) 

429 ) 

430 

431 for arrow in arrows: 

432 ax.annotate( 

433 "", 

434 xy=(arrow.head.x, arrow.head.y), 

435 xytext=(arrow.tail.x, arrow.tail.y), 

436 arrowprops={ 

437 "color": arrow.color, 

438 "width": 1.5, 

439 "headwidth": 8, 

440 "shrink": 0, 

441 }, 

442 ) 

443 # ax.annotate's arrows don't register with the legend on their own, so 

444 # proxy line handles stand in for each arrow's color/label. 

445 handles = [ 

446 plt.Line2D([0], [0], color=arrow.color, lw=1.5, label=arrow.label) 

447 for arrow in arrows 

448 ] 

449 

450 ax.set_xlim(x_min, x_max) 

451 ax.set_ylim(y_max, y_min) # inverted: image y increases downward 

452 ax.set_xlabel("x (pixels)") 

453 ax.set_ylabel("y (pixels)") 

454 # Same fully-outside legend placement as subimage_bounds_plot(), so it 

455 # never overlaps the image regardless of arrow placement. 

456 ax.legend( 

457 handles=handles, loc="center left", bbox_to_anchor=(1.02, 0.5), fontsize=8 

458 ) 

459 

460 plt.savefig(path, dpi=dpi, bbox_inches="tight") 

461 plt.close(fig) 

462 

463 

464def checkerboard( 

465 *, width: int, height: int, count_x: int = 8, count_y: int = 8 

466) -> np.ndarray: 

467 """Generate a black-and-white checkerboard test image. 

468 

469 Same parameters and pixel values as the `dictk checkerboard` CLI 

470 command, minus the file write — see the `dictk` package docstring for 

471 why the API layer stops at the array. count_x and count_y are named for 

472 the rectangles they divide the image into, not "squares": unequal 

473 counts (or a width/height ratio that doesn't match count_x/count_y) 

474 produce rectangular cells, not square ones. 

475 

476 Args: 

477 width: Image width in pixels. 

478 height: Image height in pixels. 

479 count_x: Number of rectangles along the width. Must be >= 1. 

480 count_y: Number of rectangles along the height. Must be >= 1. 

481 

482 Returns: 

483 A 2D uint8 array of shape (height, width) with values 0 or 255. 

484 

485 Raises: 

486 ValueError: If count_x or count_y is less than 1. 

487 """ 

488 if count_x < 1: 

489 raise ValueError(f"count_x {count_x} must be >= 1") 

490 if count_y < 1: 

491 raise ValueError(f"count_y {count_y} must be >= 1") 

492 

493 rows = (np.arange(height) * count_y // height) % 2 

494 cols = (np.arange(width) * count_x // width) % 2 

495 pattern = np.logical_xor(rows[:, None], cols[None, :]) 

496 return (pattern * 255).astype(np.uint8) 

497 

498 

499def astronaut(*, width: int = 512, height: int = 512) -> np.ndarray: 

500 """Load a bundled real-world grayscale reference image. 

501 

502 Same parameters and pixel values as the `dictk astronaut` CLI 

503 command, minus the file write — see the `dictk` package docstring for 

504 why the API layer stops at the array. 

505 

506 The source is a NASA portrait of astronaut Eileen Collins, from the 

507 NASA Great Images database ("No known copyright restrictions, released 

508 into the public domain."). It's bundled as a color image and converted 

509 with `rgba_to_gray`; unlike `rosta` and `checkerboard`, it isn't 

510 procedurally generated, so `width`/`height` other than the native 

511 512x512 resize the source image (via `scipy.ndimage.zoom`) rather than 

512 computing a fresh pattern at that size. 

513 

514 Args: 

515 width: Image width in pixels. Must be >= 1. 

516 height: Image height in pixels. Must be >= 1. 

517 

518 Returns: 

519 A 2D uint8 array of shape (height, width). 

520 

521 Raises: 

522 ValueError: If width or height is less than 1. 

523 """ 

524 if width < 1: 

525 raise ValueError(f"width {width} must be >= 1") 

526 if height < 1: 

527 raise ValueError(f"height {height} must be >= 1") 

528 

529 asset_path = importlib.resources.files("dictk") / "data" / "astronaut.png" 

530 with importlib.resources.as_file(asset_path) as path: 

531 color = read(path=path) 

532 gray = rgba_to_gray(color) 

533 

534 native_height, native_width = gray.shape 

535 if (width, height) == (native_width, native_height): 

536 return gray 

537 

538 zoom_factors = (height / native_height, width / native_width) 

539 resized = zoom(gray.astype(np.float64), zoom_factors, order=3) 

540 return np.clip(resized, 0, 255).astype(np.uint8) 

541 

542 

543def is_rgba(arr: np.ndarray) -> bool: 

544 """Check whether an image array is in RGB or RGBA format. 

545 

546 Args: 

547 arr: Input image array. 

548 

549 Returns: 

550 True if the array is 3D with 3 or 4 channels, False otherwise. 

551 """ 

552 return arr.ndim == 3 and arr.shape[2] in (3, 4) 

553 

554 

555def rgba_to_gray(arr: np.ndarray) -> np.ndarray: 

556 """Convert an RGB(A) image to grayscale by averaging the RGB channels. 

557 

558 Args: 

559 arr: Input image array, either 2D (grayscale) or 3D (color). 

560 

561 Returns: 

562 A 2D grayscale image array. If the input is already 2D, it is 

563 returned unchanged. 

564 

565 Raises: 

566 ValueError: If the array is neither 2D nor a 3-or-4-channel 3D array. 

567 """ 

568 if arr.ndim == 2: 

569 return arr 

570 

571 if is_rgba(arr): 

572 return np.mean(arr[:, :, :3], axis=2).astype(arr.dtype) 

573 

574 raise ValueError( 

575 "Input array must be either 2D (grayscale) or 3D with 3 or 4 channels (color)." 

576 ) 

577 

578 

579def combine(*, a: np.ndarray, b: np.ndarray) -> np.ndarray: 

580 """Combine two images by averaging their pixel values. 

581 

582 Args: 

583 a: First image, 2D grayscale or 3D color. 

584 b: Second image, same shape as `a` once converted to grayscale. 

585 

586 Returns: 

587 A 2D uint8 array, normalized to the range [0, 255]. 

588 

589 Raises: 

590 ValueError: If the grayscale-converted images don't share a shape. 

591 """ 

592 gray_a = rgba_to_gray(a) 

593 gray_b = rgba_to_gray(b) 

594 

595 if gray_a.shape != gray_b.shape: 

596 raise ValueError( 

597 f"shape mismatch: a.shape={gray_a.shape}, b.shape={gray_b.shape}" 

598 ) 

599 

600 combined = gray_a.astype(np.float64) + gray_b.astype(np.float64) 

601 return (combined / combined.max() * 255).astype(np.uint8) 

602 

603 

604def brightness(*, arr: np.ndarray, factor: float) -> np.ndarray: 

605 """Adjust image brightness by an additive shift, clipped to [0, 255]. 

606 

607 Brightness *translates* the pixel-intensity histogram: every pixel is 

608 shifted by the same amount, so dark areas lighten right along with 

609 bright ones (unlike a multiplicative scale, where a black pixel would 

610 stay exactly black). factor=1.0 leaves the image unchanged; factor=1.5 

611 shifts every pixel by +127.5*0.5 = +63.75, factor=2.0 by +127.5. 

612 

613 Args: 

614 arr: A 2D grayscale image array. 

615 factor: Brightness factor. 1.0 is unchanged; > 1.0 brightens 

616 (shifts toward white); < 1.0 darkens (shifts toward black). 

617 

618 Returns: 

619 A 2D uint8 array, same shape as `arr`. 

620 """ 

621 max_pixel_value = 255.0 

622 offset = (factor - 1.0) * (max_pixel_value / 2.0) 

623 shifted = arr.astype(np.float64) + offset 

624 return np.clip(shifted, 0, max_pixel_value).astype(np.uint8) 

625 

626 

627def contrast(*, arr: np.ndarray, factor: float) -> np.ndarray: 

628 """Adjust image contrast by scaling around the mean, clipped to [0, 255]. 

629 

630 Contrast *stretches* the pixel-intensity histogram outward from its own 

631 mean, rather than shifting it: the mean stays roughly the same, while 

632 values spread further from it. factor=1.0 leaves the image unchanged; 

633 factor=0.0 collapses every pixel to the mean (a flat gray image); 

634 factor > 1.0 pushes values further toward 0 and 255. 

635 

636 Args: 

637 arr: A 2D grayscale image array. 

638 factor: Contrast factor. 1.0 is unchanged; > 1.0 increases 

639 contrast; between 0.0 and 1.0 decreases it. 

640 

641 Returns: 

642 A 2D uint8 array, same shape as `arr`. 

643 """ 

644 mean = arr.astype(np.float64).mean() 

645 stretched = (arr.astype(np.float64) - mean) * factor + mean 

646 return np.clip(stretched, 0, 255).astype(np.uint8) 

647 

648 

649def _backward_map( 

650 arr: np.ndarray, xs_source: np.ndarray, ys_source: np.ndarray 

651) -> np.ndarray: 

652 """Sample `arr` via bilinear interpolation at (xs_source, ys_source). 

653 

654 Shared by geometric-transform functions (`stretch`, `translate`, ...): 

655 each computes where every output pixel's source coordinate falls under 

656 its own inverse transform, then hands the resulting coordinate grids 

657 here to do the actual sampling. Coordinates outside `arr`'s bounds are 

658 filled with black (0). 

659 

660 Args: 

661 arr: A 2D grayscale image array. 

662 xs_source: Source x-coordinate for each output pixel, shape (height, width). 

663 ys_source: Source y-coordinate for each output pixel, shape (height, width). 

664 

665 Returns: 

666 A 2D uint8 array, same shape as `arr`. 

667 """ 

668 height, width = arr.shape 

669 interpolator = RegularGridInterpolator( 

670 points=(np.arange(height), np.arange(width)), 

671 values=arr.astype(np.float64), 

672 method="linear", 

673 bounds_error=False, 

674 fill_value=0.0, 

675 ) 

676 points = np.stack((ys_source.ravel(), xs_source.ravel()), axis=1) 

677 deformed = interpolator(points).reshape(arr.shape) 

678 return np.clip(deformed, 0, 255).astype(np.uint8) 

679 

680 

681def stretch( 

682 *, arr: np.ndarray, factor_x: float = 1.0, factor_y: float = 1.0 

683) -> np.ndarray: 

684 """Apply a uniaxial or biaxial stretch, pivoting on the image origin. 

685 

686 Mimics a continuum-mechanics stretch deformation gradient 

687 diag(factor_x, factor_y), anchored at the image's top-left corner 

688 (x=0, y=0): that corner stays fixed, and content grows (factor > 1.0) 

689 or shrinks (factor < 1.0) away from it along each axis. Uses backward 

690 mapping — for each output pixel, the inverse of the stretch locates 

691 its source coordinate in `arr`, with bilinear interpolation for 

692 non-integer source coordinates — so the result has no gaps, unlike 

693 naively moving each source pixel forward. A factor < 1.0 shrinks 

694 content toward the origin, leaving black (fill value 0) margins along 

695 the far (bottom/right) edges. 

696 

697 Args: 

698 arr: A 2D grayscale image array. 

699 factor_x: Stretch factor along the x-axis. Must be > 0. 

700 factor_y: Stretch factor along the y-axis. Must be > 0. 

701 

702 Returns: 

703 A 2D uint8 array, same shape as `arr`. 

704 

705 Raises: 

706 ValueError: If factor_x or factor_y is <= 0. 

707 """ 

708 if factor_x <= 0: 

709 raise ValueError(f"factor_x {factor_x} must be > 0") 

710 if factor_y <= 0: 

711 raise ValueError(f"factor_y {factor_y} must be > 0") 

712 

713 height, width = arr.shape 

714 xs, ys = np.meshgrid(np.arange(width), np.arange(height)) 

715 

716 # Backward mapping, pivoting on the origin (0, 0): for each output 

717 # pixel, the inverse of the stretch gives the coordinate to sample 

718 # from in the original image. 

719 xs_source = xs / factor_x 

720 ys_source = ys / factor_y 

721 

722 return _backward_map(arr, xs_source, ys_source) 

723 

724 

725def translate(*, arr: np.ndarray, dx: float = 0.0, dy: float = 0.0) -> np.ndarray: 

726 """Apply a rigid-body translation: every pixel shifts by (dx, dy). 

727 

728 A pure displacement, with no change in shape or size — the simplest 

729 transformation category. Uses the same backward-mapping approach as 

730 `stretch`, so non-integer displacements are handled with bilinear 

731 interpolation rather than rounding. Content shifted in from outside 

732 the original bounds is filled with black (fill value 0). 

733 

734 Args: 

735 arr: A 2D grayscale image array. 

736 dx: Displacement in pixels along the x-axis; positive moves 

737 content right. 

738 dy: Displacement in pixels along the y-axis; positive moves 

739 content down. 

740 

741 Returns: 

742 A 2D uint8 array, same shape as `arr`. 

743 """ 

744 height, width = arr.shape 

745 xs, ys = np.meshgrid(np.arange(width), np.arange(height)) 

746 

747 # Backward mapping: the source of output pixel (x, y) is (x - dx, y - dy). 

748 xs_source = xs - dx 

749 ys_source = ys - dy 

750 

751 return _backward_map(arr, xs_source, ys_source) 

752 

753 

754def rotate(*, arr: np.ndarray, angle: float) -> np.ndarray: 

755 """Apply a rigid-body rotation, pivoting on the image origin. 

756 

757 Rotates content by `angle` degrees, positive counterclockwise, 

758 pivoting on the image's top-left corner (0, 0) rather than its 

759 center — consistent with `stretch` and `translate`'s pivot choice in 

760 this codebase, but unlike a typical "object spins in place" rotation 

761 example: most content swings away from that fixed corner, similar to 

762 a door on a hinge. Uses the same backward-mapping approach as 

763 `stretch` and `translate`, so non-integer source coordinates are 

764 bilinearly interpolated, and any pixel with no corresponding source 

765 coordinate within `arr`'s bounds is filled with black (fill value 0). 

766 

767 Args: 

768 arr: A 2D grayscale image array. 

769 angle: Rotation angle in degrees; positive is counterclockwise. 

770 

771 Returns: 

772 A 2D uint8 array, same shape as `arr`. 

773 """ 

774 height, width = arr.shape 

775 xs, ys = np.meshgrid(np.arange(width), np.arange(height)) 

776 

777 theta = np.deg2rad(angle) 

778 cos_theta = np.cos(theta) 

779 sin_theta = np.sin(theta) 

780 

781 # Backward mapping: applying the inverse (-angle) rotation to each 

782 # output coordinate gives the coordinate to sample from. 

783 xs_source = cos_theta * xs + sin_theta * ys 

784 ys_source = -sin_theta * xs + cos_theta * ys 

785 

786 return _backward_map(arr, xs_source, ys_source) 

787 

788 

789def shear(*, arr: np.ndarray, shear_x: float = 0.0, shear_y: float = 0.0) -> np.ndarray: 

790 """Apply a simple shear, pivoting on the image origin. 

791 

792 Mimics a continuum-mechanics shear deformation gradient 

793 [[1, shear_x], [shear_y, 1]], anchored at the image's top-left corner 

794 (x=0, y=0), consistent with `stretch`, `translate`, and `rotate`'s 

795 pivot choice in this codebase: horizontal planes slide relative to 

796 each other by an amount proportional to their y-coordinate 

797 (shear_x), and/or vertical planes slide by an amount proportional to 

798 their x-coordinate (shear_y). Uses the same backward-mapping approach 

799 as the other transform functions, so non-integer source coordinates 

800 are bilinearly interpolated, and any pixel with no corresponding 

801 source coordinate within `arr`'s bounds is filled with black (fill 

802 value 0). 

803 

804 Args: 

805 arr: A 2D grayscale image array. 

806 shear_x: Horizontal shear factor. 

807 shear_y: Vertical shear factor. 

808 

809 Returns: 

810 A 2D uint8 array, same shape as `arr`. 

811 

812 Raises: 

813 ValueError: If shear_x * shear_y == 1, which makes the 

814 deformation gradient singular (non-invertible). 

815 """ 

816 determinant = 1.0 - shear_x * shear_y 

817 if determinant == 0: 

818 raise ValueError( 

819 f"shear_x {shear_x} and shear_y {shear_y} produce a singular " 

820 "deformation gradient (shear_x * shear_y == 1)" 

821 ) 

822 

823 height, width = arr.shape 

824 xs, ys = np.meshgrid(np.arange(width), np.arange(height)) 

825 

826 # Backward mapping: invert the shear deformation gradient 

827 # [[1, shear_x], [shear_y, 1]] to find each output pixel's source. 

828 xs_source = (xs - shear_x * ys) / determinant 

829 ys_source = (-shear_y * xs + ys) / determinant 

830 

831 return _backward_map(arr, xs_source, ys_source) 

832 

833 

834def complex_deform( 

835 *, 

836 arr: np.ndarray, 

837 factor_x: float = 1.0, 

838 factor_y: float = 1.0, 

839 angle: float = 0.0, 

840) -> np.ndarray: 

841 """Apply an anisotropic stretch composed with a rotation, in one pass. 

842 

843 Composes a stretch (deformation gradient diag(factor_x, factor_y)) 

844 with a rotation (`angle` degrees, counterclockwise): the combined 

845 deformation gradient is F = R(angle) @ diag(factor_x, factor_y), i.e. 

846 the stretch is applied first and the rotation second. Applying both 

847 in a single backward-mapping pass (rather than calling `stretch` and 

848 then `rotate` separately) avoids the extra blur of interpolating 

849 twice. Represents realistic loading scenarios where materials 

850 experience multiple simultaneous deformation modes — typically the 

851 hardest case for correlation algorithms. Pivots on the image's 

852 top-left corner (0, 0), consistent with `stretch`, `translate`, 

853 `rotate`, and `shear`'s pivot choice in this codebase. 

854 

855 Args: 

856 arr: A 2D grayscale image array. 

857 factor_x: Stretch factor along the x-axis, applied before the 

858 rotation. Must be > 0. 

859 factor_y: Stretch factor along the y-axis, applied before the 

860 rotation. Must be > 0. 

861 angle: Rotation angle in degrees, applied after the stretch; 

862 positive is counterclockwise. 

863 

864 Returns: 

865 A 2D uint8 array, same shape as `arr`. 

866 

867 Raises: 

868 ValueError: If factor_x or factor_y is <= 0. 

869 """ 

870 if factor_x <= 0: 

871 raise ValueError(f"factor_x {factor_x} must be > 0") 

872 if factor_y <= 0: 

873 raise ValueError(f"factor_y {factor_y} must be > 0") 

874 

875 height, width = arr.shape 

876 xs, ys = np.meshgrid(np.arange(width), np.arange(height)) 

877 

878 theta = np.deg2rad(angle) 

879 cos_theta = np.cos(theta) 

880 sin_theta = np.sin(theta) 

881 

882 # Backward mapping for F = R(angle) @ diag(factor_x, factor_y): 

883 # F_inv = diag(1/factor_x, 1/factor_y) @ R(-angle), applied to each 

884 # output coordinate to find its source. Un-rotate first, then unscale. 

885 xs_rotated = cos_theta * xs + sin_theta * ys 

886 ys_rotated = -sin_theta * xs + cos_theta * ys 

887 xs_source = xs_rotated / factor_x 

888 ys_source = ys_rotated / factor_y 

889 

890 return _backward_map(arr, xs_source, ys_source) 

891 

892 

893def crack_dislocation(*, arr: np.ndarray, offset: float = 8.0) -> np.ndarray: 

894 """Apply a discontinuous vertical-crack displacement field. 

895 

896 Splits the image with a vertical crack at x = width / 2: the left 

897 half shifts down by `offset` pixels and the right half shifts up by 

898 `offset` pixels, producing a displacement field that jumps 

899 discontinuously across the crack line — unlike every other transform 

900 in this module, which varies smoothly. Standard DIC assumes smooth 

901 displacements and cannot capture this jump; cases like this motivate 

902 the Heaviside finite-element formulation. Uses the same backward 

903 mapping as the other transform functions, just with a piecewise 

904 (rather than single-matrix) displacement field, so non-integer 

905 source coordinates are bilinearly interpolated, and any pixel with 

906 no corresponding source coordinate within `arr`'s bounds is filled 

907 with black (fill value 0). 

908 

909 Args: 

910 arr: A 2D grayscale image array. 

911 offset: Displacement in pixels; the left half (x < width / 2) 

912 shifts down by this amount, the right half shifts up by it. 

913 

914 Returns: 

915 A 2D uint8 array, same shape as `arr`. 

916 """ 

917 height, width = arr.shape