From c7cd55a77a9542ae0e54ec403120840f3b255bbd Mon Sep 17 00:00:00 2001 From: tpietzsch Date: Tue, 3 Dec 2024 11:02:30 +0100 Subject: [PATCH 1/2] Speed up RenderTransformMeshMappingWithMasks --- .../RenderTransformMeshMappingWithMasks.java | 103 +++++++-------- .../java/org/janelia/alignment/Triangle.java | 117 ++++++++++++++++++ .../janelia/alignment/mapper/PixelMapper.java | 63 +++++++++- 3 files changed, 227 insertions(+), 56 deletions(-) create mode 100644 render-app/src/main/java/org/janelia/alignment/Triangle.java diff --git a/render-app/src/main/java/org/janelia/alignment/RenderTransformMeshMappingWithMasks.java b/render-app/src/main/java/org/janelia/alignment/RenderTransformMeshMappingWithMasks.java index 49b9ddb55..10153fb26 100644 --- a/render-app/src/main/java/org/janelia/alignment/RenderTransformMeshMappingWithMasks.java +++ b/render-app/src/main/java/org/janelia/alignment/RenderTransformMeshMappingWithMasks.java @@ -1,4 +1,4 @@ -/** +/* * License: GPL * * This program is free software; you can redistribute it and/or @@ -21,10 +21,14 @@ import java.util.concurrent.atomic.AtomicInteger; import mpicbg.models.AffineModel2D; +import mpicbg.models.NoninvertibleModelException; import mpicbg.trakem2.util.Pair; import mpicbg.util.Util; +import net.imglib2.realtransform.AffineTransform2D; +import org.janelia.alignment.Triangle.Range; import org.janelia.alignment.mapper.PixelMapper; +import org.janelia.alignment.mapper.PixelMapper.LineMapper; import org.slf4j.Logger; import org.slf4j.LoggerFactory; @@ -101,13 +105,34 @@ final public void run() { } } + /** + * Extract the inverse {@code model} transform + * (mapping target to source coordinates). + * + * @param model + * @return the inverse model transform + * @throws NoninvertibleModelException + */ + private static AffineTransform2D getTransformToSource(AffineModel2D model) throws NoninvertibleModelException { + try { + final AffineTransform2D transform = new AffineTransform2D(); + final double[] m = new double[6]; + model.toArray(m); + transform.set(m[0], m[2], m[4], m[1], m[3], m[5]); + return transform.inverse(); + } catch (final Exception e) { + throw new NoninvertibleModelException("Model not invertible."); + } + } + + private static Exception mapTriangle(final Pair ai, final PixelMapper pixelMapper) { - final int w = pixelMapper.getTargetWidth() - 1; - final int h = pixelMapper.getTargetHeight() - 1; - final double[][] pq = ai.b; + final AffineModel2D model = ai.a; + + final Triangle triangle = new Triangle(pq[2][0], pq[3][0], pq[2][1], pq[3][1], pq[2][2], pq[3][2]); final double[] min = new double[2]; final double[] max = new double[2]; @@ -115,61 +140,29 @@ private static Exception mapTriangle(final Pair ai, final int minX = Math.max(0, Util.roundPos(min[0])); final int minY = Math.max(0, Util.roundPos(min[1])); - final int maxX = Math.min(w, Util.roundPos(max[0])); - final int maxY = Math.min(h, Util.roundPos(max[1])); + final int maxX = Math.min(pixelMapper.getTargetWidth() - 1, Util.roundPos(max[0])); + final int maxY = Math.min(pixelMapper.getTargetHeight() - 1, Util.roundPos(max[1])); + + final AffineTransform2D affineTransform2D; + try { + affineTransform2D = getTransformToSource(model); + } catch (final NoninvertibleModelException e) { + return e; + } + final LineMapper lineMapper = pixelMapper.createLineMapper(); + final double dx = affineTransform2D.d(0).getDoublePosition(0); + final double dy = affineTransform2D.d(0).getDoublePosition(1); final double[] source = new double[2]; - - // TODO: ask Saalfeld why we can't just throw the exception - are there common cases where ignoring makes sense? - Exception lastIgnoredException = null; - - if (pixelMapper.isMappingInterpolated()) { - - for (int targetY = minY; targetY <= maxY; ++targetY) { - for (int targetX = minX; targetX <= maxX; ++targetX) { - - if (RenderTransformMesh.isInTargetTriangle(pq, targetX, targetY)) { - - source[0] = targetX; - source[1] = targetY; - - try { - ai.a.applyInverseInPlace(source); - } catch (final Exception e) { - lastIgnoredException = e; - continue; - } - - pixelMapper.mapInterpolated(source[0], source[1], targetX, targetY); - } - } - } - - } else { - - for (int targetY = minY; targetY <= maxY; ++targetY) { - for (int targetX = minX; targetX <= maxX; ++targetX) { - - if (RenderTransformMesh.isInTargetTriangle(pq, targetX, targetY)) { - - source[0] = targetX; - source[1] = targetY; - - try { - ai.a.applyInverseInPlace(source); - } catch (final Exception e) { - lastIgnoredException = e; - continue; - } - - pixelMapper.map(source[0], source[1], targetX, targetY); - } - } - } - + for (int targetY = minY; targetY <= maxY; ++targetY) { + final Range xRange = triangle.intersect(targetY, minX, maxX + 1); + source[0] = xRange.from(); + source[1] = targetY; + affineTransform2D.apply(source, source); + lineMapper.map(source[0], source[1], dx, dy, xRange.from(), targetY, xRange.length()); } - return lastIgnoredException; + return null; } private static final Logger LOG = LoggerFactory.getLogger(RenderTransformMeshMappingWithMasks.class); diff --git a/render-app/src/main/java/org/janelia/alignment/Triangle.java b/render-app/src/main/java/org/janelia/alignment/Triangle.java new file mode 100644 index 000000000..17e59aaf9 --- /dev/null +++ b/render-app/src/main/java/org/janelia/alignment/Triangle.java @@ -0,0 +1,117 @@ +package org.janelia.alignment; + +import static org.janelia.alignment.Triangle.WindingOrder.CCW; + +/** + * A triangle in 2D, comprising points A, B, C specified by {@code ax, ay, + * bx, by, cx, cy}. The points are stored in counter-clockwise order (the + * constructor re-orders points if necessary). + *

+ * When intersecting horizontal scan-lines with a triangle edge, the points + * {@code (a, b)} of the edge are always re-ordered such that {@code ay <= + * by}. This is to prevent potential gaps between triangles introduced by + * rounding errors. When two neighboring triangles contain the same edge + * (same points in the same order), an X position on a horizontal line will + * never be between triangles. + */ +class Triangle { + + private final double ax; + private final double ay; + private final double bx; + private final double by; + private final double cx; + private final double cy; + + enum WindingOrder { + CW, CCW; + + static WindingOrder of(final double ax, final double ay, final double bx, final double by, final double cx, final double cy) { + final double P = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax); + return P > 0 ? CW : CCW; + } + } + + Triangle(final double ax, final double ay, final double bx, final double by, final double cx, final double cy) { + if (WindingOrder.of(ax, ay, bx, by, cx, cy) == CCW) { + this.ax = ax; + this.ay = ay; + this.bx = bx; + this.by = by; + } else { + this.ax = bx; + this.ay = by; + this.bx = ax; + this.by = ay; + } + this.cx = cx; + this.cy = cy; + } + + static class Range { + private final int from; + private final int length; + + private Range(final int min, final int max) { + this.from = min; + this.length = max - min; + } + + /** + * min (inclusive) of this Range. + */ + int from() { + return from; + } + + /** + * number of elements in this Range. + */ + int length() { + return length; + } + } + + /** + * Compute the intersection of horizontal scan-line at {@code y} with this triangle. + * The resulting discrete X range can is [{@code intersectionMinX(), intersectionMaxX()}) + * + * @param y Y coordinate of scanline to intersect with this triangle + * @param rangeMinX inclusive lower bound (minimum discrete X coordinate to consider) + * @param rangeMaxX exclusive upper bound (maximum discrete X coordinate to consider + 1) + * @return the discrete X range of the intersected scanline. + */ + Range intersect(final double y, final int rangeMinX, final int rangeMaxX) { + final double[] bounds = new double[]{rangeMinX, rangeMaxX}; + intersect(y, bounds); + return new Range((int) Math.ceil(bounds[0]), (int) Math.ceil(bounds[1])); + } + + void intersect(final double y, final double[] bounds) { + updateRange(ax, ay, bx, by, y, bounds); + updateRange(bx, by, cx, cy, y, bounds); + updateRange(cx, cy, ax, ay, y, bounds); + } + + private static void updateRange(final double ax, final double ay, final double bx, final double by, final double y, final double[] bounds) { + final boolean flip = ay > by; + final double px = flip ? ax : bx; + final double qx = flip ? bx : ax; + final double py = flip ? ay : by; + final double qy = flip ? by : ay; + + final double s = (qx - px) / (qy - py); + if (Double.isInfinite(s)) { + return; + } else if (Double.isNaN(s)) { + throw new IllegalArgumentException("zero length triangle edge"); + } + + final double x = px + s * (y - py); + if (flip) { + bounds[1] = Math.min(bounds[1], x); + } else { + bounds[0] = Math.max(bounds[0], x); + } + } +} diff --git a/render-app/src/main/java/org/janelia/alignment/mapper/PixelMapper.java b/render-app/src/main/java/org/janelia/alignment/mapper/PixelMapper.java index 94fd20c44..9a03a1363 100644 --- a/render-app/src/main/java/org/janelia/alignment/mapper/PixelMapper.java +++ b/render-app/src/main/java/org/janelia/alignment/mapper/PixelMapper.java @@ -50,4 +50,65 @@ void mapInterpolated(final double sourceX, final int targetX, final int targetY); -} \ No newline at end of file + /** + * Maps values for horizontal ranges of target pixels. + */ + interface LineMapper { + /** + * Maps value for a line of pixels in the target {@code (targetX + i, targetY)}, where {@code 0 <= i < length}. + *

+ * Stepping one pixel in X in the target, means stepping {@code (sourceStepX, sourceStepY)} in the source. + * That is, target pixel {@code (targetX + i, targetY)} corresponds to source pixel {@code (sourceX + i * sourceStepX, sourceY + i * sourceStepY)}. + *

+ * The interpolation method is determined when constructing the {@code LineMapper} instance is constructed (see {@link #createLineMapper(boolean)}). + * + * @param sourceX source X coordinate. + * @param sourceY source Y coordinate. + * @param sourceStepX source X offset corresponding to stepping 1 pixel in X in the target. + * @param sourceStepY source Y offset corresponding to stepping 1 pixel in X in the target. + * @param targetX local target X coordinate. + * @param targetY local target Y coordinate. + * @param length number of pixels to map. + */ + void map(double sourceX, double sourceY, double sourceStepX, double sourceStepY, int targetX, int targetY, int length); + } + + /** + * Create a {@code LineMapper}. + *

+ * If {@link #isMappingInterpolated()}{@code ==true} the {@code LineMapper} will use linear interpolation, + * otherwise nearest-neighbor interpolation. + * + * @return a new {@code LineMapper} + */ + default LineMapper createLineMapper() { + return createLineMapper(isMappingInterpolated()); + } + + /** + * Create a {@code LineMapper}. + * + * @param isMappingInterpolated if {@code true} the {@code LineMapper} will use linear interpolation, + * if {@code false} the {@code LineMapper} will use nearest-neighbor interpolation. + * @return a new {@code LineMapper} + */ + default LineMapper createLineMapper(final boolean isMappingInterpolated) { + if (isMappingInterpolated) { + return (sx, sy, sdx, sdy, tx, ty, length) -> { + for (int x = tx; x < (tx + length); ++x) { + PixelMapper.this.mapInterpolated(sx, sy, x, ty); + sx += sdx; + sy += sdy; + } + }; + } else { + return (sx, sy, sdx, sdy, tx, ty, length) -> { + for (int x = tx; x < (tx + length); ++x) { + PixelMapper.this.map(sx, sy, x, ty); + sx += sdx; + sy += sdy; + } + }; + } + } +} From 881edd0dc2c2e4b6412d87b4fd6e73f84197f47c Mon Sep 17 00:00:00 2001 From: tpietzsch Date: Mon, 6 Jan 2025 13:48:22 +0100 Subject: [PATCH 2/2] Make the tests pass --- .../RenderTransformMeshMappingWithMasks.java | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/render-app/src/main/java/org/janelia/alignment/RenderTransformMeshMappingWithMasks.java b/render-app/src/main/java/org/janelia/alignment/RenderTransformMeshMappingWithMasks.java index 10153fb26..2cd3eba2f 100644 --- a/render-app/src/main/java/org/janelia/alignment/RenderTransformMeshMappingWithMasks.java +++ b/render-app/src/main/java/org/janelia/alignment/RenderTransformMeshMappingWithMasks.java @@ -155,11 +155,19 @@ private static Exception mapTriangle(final Pair ai, final double dy = affineTransform2D.d(0).getDoublePosition(1); final double[] source = new double[2]; for (int targetY = minY; targetY <= maxY; ++targetY) { - final Range xRange = triangle.intersect(targetY, minX, maxX + 1); + final Range xRange = triangle.intersect(targetY, minX, maxX); + /* TODO: Actually, "maxX" should be "maxX + 1" here, but we compensate for the additional "+1" below. + "maxX + 1" would be correct because intersect range upper bound is exclusive. + */ source[0] = xRange.from(); source[1] = targetY; affineTransform2D.apply(source, source); - lineMapper.map(source[0], source[1], dx, dy, xRange.from(), targetY, xRange.length()); + lineMapper.map(source[0], source[1], dx, dy, xRange.from(), targetY, + xRange.length() + 1 + /* TODO: This "+1" is to fix the tests and be compatible to the old code. + However, it leads to pixels on the edge of neighboring triangles to be drawn twice. + */ + ); } return null;