Coverage for paintmaskanimator/mask_transform.py: 93%
270 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 04:50 +0000
1"""Binary-preserving mask transformation (inverse mapping + area coverage).
3The TP-quality transform used to project every color mask *forward* through the
4deformation (``QPainter.setTransform`` / per-mask triangle rasterization) and
5then rebuild a binary image from blurred grayscale masks. Forward mapping skips
6or double-writes output pixels depending on the rotation angle, and the
7blur-then-``> 0`` threshold inflated every shape by the blur radius, so rotated
8edges came out as irregular noise instead of a clean staircase.
10This module replaces that with the standard approach for binary raster
11resampling:
131. **Inverse mapping** - every *output* pixel is projected back into source
14 space, so no output pixel is ever skipped or written twice.
152. **Area coverage** - each output pixel is probed at ``n x n`` subsample
16 positions and each color counts how many probes it won.
173. **Majority resolution** - the color covering most of the pixel takes it
18 whole. For a two-region boundary that is exactly a 50% area threshold, and
19 the result is strictly binary: no intermediate colors, no partial alpha.
21Everything here is pure NumPy: no Qt, no Pillow. The functions operate on a
22*label image* (``uint16``, 0 = outside/background, 1..N = color ids) instead of
23one grayscale mask per color, which also makes the cost independent of the
24palette size.
25"""
26from __future__ import annotations
28import math
30import numpy as np
32from . import native
34#: Subsample probes per axis. 4 -> 16 probes per output pixel, which resolves
35#: coverage in 1/16 steps: enough to place a 50% boundary within a pixel.
36DEFAULT_SUBSAMPLES = 4
38#: Output rows processed per band. Bounds peak memory: a band allocates
39#: ``band_rows * width * subsamples^2`` entries rather than the whole image.
40DEFAULT_BAND_ROWS = 64
43def homography_matrix(source_points, target_points):
44 """Return the 3x3 matrix mapping ``source_points`` onto ``target_points``.
46 Points are ``(x, y)`` pairs. Four correspondences are required. Raises
47 ``numpy.linalg.LinAlgError`` for degenerate quads.
48 """
49 rows = []
50 values = []
51 for (x, y), (u, v) in zip(source_points, target_points):
52 x, y, u, v = float(x), float(y), float(u), float(v)
53 rows.append([x, y, 1, 0, 0, 0, -u * x, -u * y])
54 values.append(u)
55 rows.append([0, 0, 0, x, y, 1, -v * x, -v * y])
56 values.append(v)
57 solved = np.linalg.solve(
58 np.asarray(rows, dtype=np.float64),
59 np.asarray(values, dtype=np.float64),
60 )
61 matrix = np.empty((3, 3), dtype=np.float64)
62 matrix.flat[:8] = solved
63 matrix[2, 2] = 1.0
64 return matrix
67def _probe_coordinates(start, stop, subsamples):
68 """Probe centers, in output-pixel units, for the range ``[start, stop)``."""
69 step = 1.0 / float(subsamples)
70 offsets = (np.arange(subsamples, dtype=np.float64) + 0.5) * step
71 return (
72 np.arange(start, stop, dtype=np.float64)[:, None] + offsets[None, :]
73 ).reshape(-1)
76class QuadMapper:
77 """Inverse map for the free / scale / rotate transform (a planar quad).
79 ``target_quad`` are the four destination corners in output-image
80 coordinates, in the same order as the source rectangle corners
81 (top-left, top-right, bottom-right, bottom-left).
82 """
84 def __init__(self, target_quad, source_width, source_height):
85 self.source_width = int(source_width)
86 self.source_height = int(source_height)
87 source_quad = [
88 (0.0, 0.0),
89 (float(source_width), 0.0),
90 (float(source_width), float(source_height)),
91 (0.0, float(source_height)),
92 ]
93 target = [(float(p[0]), float(p[1])) for p in target_quad[:4]]
94 # Solve destination -> source directly: the inverse map is what the
95 # renderer needs, and solving it straight avoids inverting a matrix
96 # that may be near-singular in one direction only.
97 self.matrix = homography_matrix(target, source_quad)
99 def map_band(self, width, y0, y1, subsamples):
100 # Each coordinate is separable - ``m0*x + m1*y + m2`` - so the per-axis
101 # terms are computed once as vectors and only combined when the 2-D
102 # probe grid is finally needed. That keeps one full-size temporary per
103 # coordinate instead of one per term.
104 probe_x = _probe_coordinates(0, width, subsamples)
105 probe_y = _probe_coordinates(y0, y1, subsamples)
106 matrix = self.matrix
108 def combine(row):
109 # Kept in float64: this path is the reference the Rust renderer is
110 # checked against, so probe positions must not round differently.
111 along_x = matrix[row, 0] * probe_x
112 along_y = matrix[row, 1] * probe_y + matrix[row, 2]
113 return along_x[None, :] + along_y[:, None]
115 source_x = combine(0)
116 source_y = combine(1)
117 if abs(matrix[2, 0]) < 1e-12 and abs(matrix[2, 1]) < 1e-12:
118 # Affine (scale / rotate / shear): the perspective divisor is
119 # constant, so skip the division across the whole grid.
120 scale = float(matrix[2, 2])
121 if abs(scale) < 1e-12: 121 ↛ 122line 121 didn't jump to line 122 because the condition on line 121 was never true
122 return source_x, source_y, np.zeros(source_x.shape, dtype=bool)
123 if abs(scale - 1.0) > 1e-12: 123 ↛ 124line 123 didn't jump to line 124 because the condition on line 123 was never true
124 source_x /= scale
125 source_y /= scale
126 return source_x, source_y, True
127 denominator = combine(2)
128 # A zero denominator is the projective horizon; mark those probes
129 # invalid rather than letting the division produce infinities.
130 safe = np.abs(denominator) > 1e-12
131 np.copyto(denominator, 1.0, where=~safe)
132 source_x /= denominator
133 source_y /= denominator
134 return source_x, source_y, safe
136 def forward(self, source_x, source_y):
137 """Map source positions to output positions (the inverse direction)."""
138 try:
139 matrix = np.linalg.inv(self.matrix)
140 except np.linalg.LinAlgError:
141 return None, None
142 denominator = (
143 matrix[2, 0] * source_x + matrix[2, 1] * source_y + matrix[2, 2]
144 )
145 denominator = np.where(np.abs(denominator) > 1e-12, denominator, np.nan)
146 return (
147 (matrix[0, 0] * source_x + matrix[0, 1] * source_y + matrix[0, 2])
148 / denominator,
149 (matrix[1, 0] * source_x + matrix[1, 1] * source_y + matrix[1, 2])
150 / denominator,
151 )
154class MeshMapper:
155 """Inverse map for the mesh transform, from the dense curved grid.
157 The dense grid is rasterized triangle by triangle in *output* space; the
158 barycentric weights of each covered output probe give its source position.
159 That is already a backward map, so unlike the quad case there is nothing to
160 invert - but it must be rasterized rather than evaluated analytically.
162 ``dense_points`` are ``(x, y)`` pairs in output-image coordinates, laid out
163 row-major over ``dense_cols x dense_rows``.
164 """
166 def __init__(
167 self, dense_points, dense_cols, dense_rows, source_width, source_height
168 ):
169 self.points = np.asarray(dense_points, dtype=np.float64).reshape(-1, 2)
170 self.dense_cols = int(dense_cols)
171 self.dense_rows = int(dense_rows)
172 self.source_width = int(source_width)
173 self.source_height = int(source_height)
175 def map_band(self, width, y0, y1, subsamples):
176 out_width = int(width) * int(subsamples)
177 out_height = (int(y1) - int(y0)) * int(subsamples)
178 source_x = np.zeros((out_height, out_width), dtype=np.float64)
179 source_y = np.zeros((out_height, out_width), dtype=np.float64)
180 valid = np.zeros((out_height, out_width), dtype=bool)
182 scale = float(subsamples)
183 # Output-space triangle vertices scaled into probe coordinates, with the
184 # band origin removed so each band indexes from 0.
185 points = self.points * scale
186 origin_y = float(y0) * scale
187 cols, rows = self.dense_cols, self.dense_rows
188 span_x = max(1, cols - 1)
189 span_y = max(1, rows - 1)
190 last_x = float(self.source_width - 1)
191 last_y = float(self.source_height - 1)
193 def raster(target_indices, source_triangle):
194 triangle = points[list(target_indices)]
195 target_x = triangle[:, 0]
196 target_y = triangle[:, 1] - origin_y
197 min_x = max(0, int(math.floor(float(target_x.min()))))
198 max_x = min(out_width - 1, int(math.ceil(float(target_x.max()))))
199 min_y = max(0, int(math.floor(float(target_y.min()))))
200 max_y = min(out_height - 1, int(math.ceil(float(target_y.max()))))
201 if max_x < min_x or max_y < min_y:
202 return
203 denominator = (
204 (target_y[1] - target_y[2]) * (target_x[0] - target_x[2])
205 + (target_x[2] - target_x[1]) * (target_y[0] - target_y[2])
206 )
207 if abs(denominator) < 1e-8: 207 ↛ 208line 207 didn't jump to line 208 because the condition on line 207 was never true
208 return
209 yy, xx = np.mgrid[min_y:max_y + 1, min_x:max_x + 1]
210 # Probe centers sit half a subsample step inside the cell.
211 xx = xx + 0.5
212 yy = yy + 0.5
213 weight0 = (
214 (target_y[1] - target_y[2]) * (xx - target_x[2])
215 + (target_x[2] - target_x[1]) * (yy - target_y[2])
216 ) / denominator
217 weight1 = (
218 (target_y[2] - target_y[0]) * (xx - target_x[2])
219 + (target_x[0] - target_x[2]) * (yy - target_y[2])
220 ) / denominator
221 weight2 = 1.0 - weight0 - weight1
222 inside = (
223 (weight0 >= -1e-6) & (weight1 >= -1e-6) & (weight2 >= -1e-6)
224 )
225 if not np.any(inside):
226 return
227 mapped_x = (
228 weight0 * source_triangle[0][0]
229 + weight1 * source_triangle[1][0]
230 + weight2 * source_triangle[2][0]
231 )
232 mapped_y = (
233 weight0 * source_triangle[0][1]
234 + weight1 * source_triangle[1][1]
235 + weight2 * source_triangle[2][1]
236 )
237 region_x = source_x[min_y:max_y + 1, min_x:max_x + 1]
238 region_y = source_y[min_y:max_y + 1, min_x:max_x + 1]
239 region_valid = valid[min_y:max_y + 1, min_x:max_x + 1]
240 # The grid is anchored on source pixel *centers*; shift to the
241 # area convention the sampler uses (pixel n spans [n, n+1)).
242 region_x[inside] = mapped_x[inside] + 0.5
243 region_y[inside] = mapped_y[inside] + 0.5
244 region_valid[inside] = True
246 for grid_y in range(rows - 1):
247 top = last_y * grid_y / span_y
248 bottom = last_y * (grid_y + 1) / span_y
249 row_base = grid_y * cols
250 for grid_x in range(cols - 1):
251 left = last_x * grid_x / span_x
252 right = last_x * (grid_x + 1) / span_x
253 index = row_base + grid_x
254 raster(
255 (index, index + 1, index + cols + 1),
256 ((left, top), (right, top), (right, bottom)),
257 )
258 raster(
259 (index, index + cols + 1, index + cols),
260 ((left, top), (right, bottom), (left, bottom)),
261 )
262 return source_x, source_y, valid
264 def forward(self, source_x, source_y):
265 """Map source positions to output positions.
267 The dense grid samples the source on a regular lattice, so the forward
268 direction is a bilinear lookup in that grid - no inversion needed.
269 """
270 cols, rows = self.dense_cols, self.dense_rows
271 span_x = max(1, cols - 1)
272 span_y = max(1, rows - 1)
273 last_x = max(1e-9, float(self.source_width - 1))
274 last_y = max(1e-9, float(self.source_height - 1))
275 grid_x = np.clip(
276 np.asarray(source_x, dtype=np.float64) / last_x * span_x,
277 0.0, span_x,
278 )
279 grid_y = np.clip(
280 np.asarray(source_y, dtype=np.float64) / last_y * span_y,
281 0.0, span_y,
282 )
283 x0 = np.minimum(grid_x.astype(np.int64), cols - 2)
284 y0 = np.minimum(grid_y.astype(np.int64), rows - 2)
285 fraction_x = (grid_x - x0)[..., None]
286 fraction_y = (grid_y - y0)[..., None]
287 points = self.points.reshape(rows, cols, 2)
288 top = (
289 points[y0, x0] * (1.0 - fraction_x) + points[y0, x0 + 1] * fraction_x
290 )
291 bottom = (
292 points[y0 + 1, x0] * (1.0 - fraction_x)
293 + points[y0 + 1, x0 + 1] * fraction_x
294 )
295 mapped = top * (1.0 - fraction_y) + bottom * fraction_y
296 return mapped[..., 0], mapped[..., 1]
299def _sample_labels(label_image, source_x, source_y, valid):
300 """Nearest-neighbour label lookup; probes outside the source read 0.
302 ``valid`` may be a scalar ``True`` when the mapping cannot fail.
303 """
304 height, width = label_image.shape
305 x = np.floor(source_x, out=source_x).astype(np.int32)
306 y = np.floor(source_y, out=source_y).astype(np.int32)
307 inside = (x >= 0) & (x < width) & (y >= 0) & (y < height)
308 if valid is not True:
309 inside &= valid
310 np.clip(x, 0, width - 1, out=x)
311 np.clip(y, 0, height - 1, out=y)
312 # One flat gather beats 2-D fancy indexing: the index arithmetic stays in
313 # int32 and NumPy takes a single pass over the label image.
314 y *= width
315 y += x
316 labels = np.take(label_image.reshape(-1), y)
317 labels[~inside] = 0
318 return labels
321def _probe_axis(labels, subsamples):
322 """Reshape ``(H*n, W*n)`` probes into ``(H, W, n*n)`` per-pixel probes."""
323 high_rows, high_cols = labels.shape
324 rows = high_rows // subsamples
325 cols = high_cols // subsamples
326 return (
327 labels.reshape(rows, subsamples, cols, subsamples)
328 .transpose(0, 2, 1, 3)
329 .reshape(rows, cols, subsamples * subsamples)
330 )
333def _majority(values):
334 """Per-row winner of ``values`` (shape ``(M, S)``), lowest value wins ties.
336 Returns ``(winner, count)``. Implemented by sorting each row and measuring
337 run lengths, which is O(S log S) with S = 16 and fully vectorized.
338 """
339 count = values.shape[1]
340 order = np.sort(values, axis=1)
341 positions = np.arange(count)
342 starts_run = np.empty(order.shape, dtype=bool)
343 starts_run[:, 0] = True
344 np.not_equal(order[:, 1:], order[:, :-1], out=starts_run[:, 1:])
345 ends_run = np.empty(order.shape, dtype=bool)
346 ends_run[:, :-1] = starts_run[:, 1:]
347 ends_run[:, -1] = True
348 run_start = np.maximum.accumulate(
349 np.where(starts_run, positions, 0), axis=1
350 )
351 run_end = np.minimum.accumulate(
352 np.where(ends_run, positions, count - 1)[:, ::-1], axis=1
353 )[:, ::-1]
354 run_length = run_end - run_start + 1
355 # argmax returns the *first* maximum, and rows are sorted ascending, so ties
356 # resolve to the smallest value - i.e. the highest-priority rank.
357 best = np.argmax(run_length, axis=1)
358 rows = np.arange(order.shape[0])
359 return order[rows, best], run_length[rows, best]
362def _area_scale(mapper, source_width, source_height):
363 """Rough output/source length ratio, measured at the middle of the source."""
364 center_x = source_width / 2.0
365 center_y = source_height / 2.0
366 span = max(1.0, min(source_width, source_height) / 8.0)
367 xs = np.array([center_x, center_x + span, center_x])
368 ys = np.array([center_y, center_y, center_y + span])
369 mapped_x, mapped_y = mapper.forward(xs, ys)
370 if mapped_x is None or not np.all(np.isfinite(mapped_x)): 370 ↛ 371line 370 didn't jump to line 371 because the condition on line 370 was never true
371 return 1.0
372 edge_x = (mapped_x[1] - mapped_x[0], mapped_y[1] - mapped_y[0])
373 edge_y = (mapped_x[2] - mapped_x[0], mapped_y[2] - mapped_y[0])
374 area = abs(edge_x[0] * edge_y[1] - edge_x[1] * edge_y[0])
375 return math.sqrt(area) / span
378def _splat_label(label_image, mapper, label, out_width, out_height):
379 """Output pixels hit by forward-mapping every source pixel of ``label``.
381 Coverage alone cannot keep a line that shrinks below one pixel: its area
382 falls under every threshold and the line dissolves into dashes. Mapping the
383 source pixels *forward* fixes that, because a connected run of source pixels
384 lands on output pixels at most one step apart, so the line stays connected
385 however far it is shrunk. Forward mapping leaves holes when enlarging,
386 which is exactly the case where coverage already does the right thing - so
387 this is only used on a shrink, and only as a union with the coverage result.
388 """
389 source_y, source_x = np.nonzero(label_image == label)
390 if source_y.size == 0: 390 ↛ 391line 390 didn't jump to line 391 because the condition on line 390 was never true
391 return None
392 mapped_x, mapped_y = mapper.forward(
393 source_x.astype(np.float64) + 0.5, source_y.astype(np.float64) + 0.5
394 )
395 if mapped_x is None: 395 ↛ 396line 395 didn't jump to line 396 because the condition on line 395 was never true
396 return None
397 finite = np.isfinite(mapped_x) & np.isfinite(mapped_y)
398 x = np.floor(mapped_x[finite]).astype(np.int64)
399 y = np.floor(mapped_y[finite]).astype(np.int64)
400 inside = (x >= 0) & (x < out_width) & (y >= 0) & (y < out_height)
401 hit = np.zeros((out_height, out_width), dtype=bool)
402 hit[y[inside], x[inside]] = True
403 return hit
406def render_labels_numpy(
407 label_image,
408 mapper,
409 out_width,
410 out_height,
411 *,
412 subsamples=DEFAULT_SUBSAMPLES,
413 coverage_labels=(),
414 protect_labels=(),
415 band_rows=DEFAULT_BAND_ROWS,
416 progress=None,
417):
418 """Transform ``label_image`` through ``mapper`` onto an output raster.
420 ``label_image`` is ``uint16`` with 0 meaning "nothing here"; label ids must
421 already be ordered by priority (a lower id wins a tie). Returns
422 ``(labels, coverage)`` where ``labels`` is ``(out_height, out_width)``
423 ``uint16`` and ``coverage`` maps each id in ``coverage_labels`` to a
424 ``uint8`` array holding that label's true area coverage (0-255) - not a
425 blurred mask, so thresholding it stays meaningful.
427 Ids listed in ``protect_labels`` additionally keep their coverage ridge, so
428 a line thinned below one pixel by a shrink stays connected instead of
429 breaking into dashes. Protected ids must also appear in ``coverage_labels``.
430 """
431 label_image = np.ascontiguousarray(label_image, dtype=np.uint16)
432 out_width = max(1, int(out_width))
433 out_height = max(1, int(out_height))
434 subsamples = max(1, int(subsamples))
435 probes = subsamples * subsamples
437 result = np.zeros((out_height, out_width), dtype=np.uint16)
438 coverage = {
439 int(label): np.zeros((out_height, out_width), dtype=np.uint8)
440 for label in coverage_labels
441 }
443 band_rows = max(1, int(band_rows))
444 band_count = max(1, math.ceil(out_height / band_rows))
445 for band_index, y0 in enumerate(range(0, out_height, band_rows)):
446 y1 = min(out_height, y0 + band_rows)
447 if progress is not None: 447 ↛ 448line 447 didn't jump to line 448 because the condition on line 447 was never true
448 progress(band_index, band_count)
449 source_x, source_y, valid = mapper.map_band(
450 out_width, y0, y1, subsamples
451 )
452 sampled = _sample_labels(label_image, source_x, source_y, valid)
453 per_pixel = _probe_axis(sampled, subsamples)
454 rows = y1 - y0
455 flat = per_pixel.reshape(rows * out_width, probes)
457 # Only pixels straddling a boundary need their probes counted; every
458 # interior pixel has all probes on one color. Boundaries are a small
459 # fraction of any real drawing, so resolving just those is much cheaper
460 # than sorting the whole band.
461 low = flat.min(axis=1)
462 mixed = low != flat.max(axis=1)
463 winner = low
464 if probes > 1 and np.any(mixed):
465 winner = winner.copy()
466 winner[mixed], _count = _majority(flat[mixed])
467 result[y0:y1] = winner.reshape(rows, out_width)
469 for label, buffer in coverage.items():
470 hits = np.where(low == label, probes, 0).astype(np.uint16)
471 if np.any(mixed):
472 hits[mixed] = np.count_nonzero(
473 flat[mixed] == label, axis=1
474 ).astype(np.uint16)
475 buffer[y0:y1] = ((hits * 255) // probes).astype(
476 np.uint8
477 ).reshape(rows, out_width)
478 if progress is not None: 478 ↛ 479line 478 didn't jump to line 479 because the condition on line 478 was never true
479 progress(band_count, band_count)
481 _protect_lines(
482 result, coverage, label_image, mapper, protect_labels,
483 out_width, out_height,
484 )
485 return result, coverage
488def _protect_lines(
489 result, coverage, label_image, mapper, protect_labels, out_width, out_height
490):
491 """Union the forward-mapped source pixels of ``protect_labels`` into ``result``.
493 Only a shrink can thin a line below one pixel. At 1:1 or larger the
494 coverage result is already complete, and splatting there would only pad the
495 edge of lines that are perfectly fine.
496 """
497 shrinking = protect_labels and _area_scale(
498 mapper, label_image.shape[1], label_image.shape[0]
499 ) < 0.999
500 for label in protect_labels if shrinking else ():
501 label = int(label)
502 hit = _splat_label(label_image, mapper, label, out_width, out_height)
503 if hit is None: 503 ↛ 504line 503 didn't jump to line 504 because the condition on line 503 was never true
504 continue
505 result[hit] = np.uint16(label)
506 buffer = coverage.get(label)
507 if buffer is not None: 507 ↛ 500line 507 didn't jump to line 500 because the condition on line 507 was always true
508 # Report full coverage there so the caller's line-width threshold
509 # cannot drop the pixels that keep the line connected.
510 buffer[hit] = 255
513def render_labels(
514 label_image,
515 mapper,
516 out_width,
517 out_height,
518 *,
519 subsamples=DEFAULT_SUBSAMPLES,
520 coverage_labels=(),
521 protect_labels=(),
522 band_rows=DEFAULT_BAND_ROWS,
523 progress=None,
524):
525 """Transform ``label_image`` through ``mapper``; see :func:`render_labels_numpy`.
527 Both mappers are resolved by the Rust renderer. The planar quad (scale /
528 rotate / free) maps probes straight from the 3x3 inverse matrix, so no band
529 of probe coordinates is ever held in memory; the mesh is rasterized
530 triangle by triangle into one band of probe *labels*, which is a quarter of
531 what the coordinates alone would take. :func:`render_labels_numpy` stays as
532 the reference the two are checked against.
533 """
534 label_image = np.ascontiguousarray(label_image, dtype=np.uint16)
535 out_width = max(1, int(out_width))
536 out_height = max(1, int(out_height))
537 subsamples = max(1, int(subsamples))
538 # The renderer resolves the whole raster in one call, so progress can only
539 # report the step as a whole.
540 if progress is not None:
541 progress(0, 1)
542 if isinstance(mapper, QuadMapper):
543 result, coverage = native.render_quad_labels(
544 label_image, mapper.matrix, out_width, out_height,
545 subsamples, coverage_labels,
546 )
547 else:
548 result, coverage = native.render_mesh_labels(
549 label_image, mapper.points, mapper.dense_cols, mapper.dense_rows,
550 out_width, out_height, subsamples, max(1, int(band_rows)),
551 coverage_labels,
552 )
553 if progress is not None:
554 progress(1, 1)
555 _protect_lines(
556 result, coverage, label_image, mapper, protect_labels,
557 out_width, out_height,
558 )
559 return result, coverage
562def labels_to_rgba(labels, colors):
563 """Expand a label image into binary RGBA using ``colors`` indexed by label.
565 ``colors`` is an ``(N, 4) uint8`` array whose row 0 is the background.
566 Output alpha is only ever 0 or 255.
567 """
568 palette = np.asarray(colors, dtype=np.uint8)
569 labels = np.asarray(labels)
570 clipped = np.clip(labels, 0, len(palette) - 1)
571 rgba = palette[clipped]
572 return np.ascontiguousarray(rgba)