Merge pull request #199 from tpietzsch/newexport-rendermesh

Speed up RenderTransformMeshMappingWithMasks

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,75 +105,72 @@ 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<AffineModel2D, double[][]> 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];
         RenderTransformMesh.calculateTargetBoundingBox(pq, min, max);
 
         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);
+            /* 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() + 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 lastIgnoredException;
+        return null;
     }
 
     private static final Logger LOG = LoggerFactory.getLogger(RenderTransformMeshMappingWithMasks.class);

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).
+ * <p>
+ * 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);
+        }
+    }
+}

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);
 
-}
+    /**
+     * 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}.
+         * <p>
+         * 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)}.
+         * <p>
+         * 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}.
+     * <p>
+     * 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;
+                }
+            };
+        }
+    }
+}