Speed up LinearIntensityMap8BitFilter

pom.xml

@@ -5,7 +5,7 @@
     <parent>
         <groupId>org.scijava</groupId>
         <artifactId>pom-scijava</artifactId>
-        <version>38.0.1</version> <!-- 38.0.1 was released on 2024-07-08, see https://mvnrepository.com/artifact/org.scijava/pom-scijava -->
+        <version>39.0.0</version> <!-- 38.0.1 was released on 2024-07-08, see https://mvnrepository.com/artifact/org.scijava/pom-scijava -->
         <relativePath />
     </parent>
 
@@ -141,6 +141,9 @@
         -->
         <n5-version>3.2.0</n5-version>
 
+        <imglib2.version>7.1.4</imglib2.version>
+        <imglib2-algorithm.version>0.17.2</imglib2-algorithm.version>
+
         <jackson-version>2.14.3</jackson-version> <!-- NOTE: 2.15.3 causes NullPointerException in maven enforcer plugin -->
 
         <swagger-version>1.6.2</swagger-version>

render-app/pom.xml

@@ -103,6 +103,11 @@
             <artifactId>imglib2</artifactId>
         </dependency>
 
+        <dependency>
+            <groupId>net.imglib2</groupId>
+            <artifactId>imglib2-algorithm</artifactId>
+        </dependency>
+
         <dependency>
             <groupId>net.imglib2</groupId>
             <artifactId>imglib2-realtransform</artifactId>

render-app/src/main/java/org/janelia/alignment/filter/LinearIntensityMap8BitFilter.java

@@ -2,16 +2,22 @@
 
 import ij.ImagePlus;
 import ij.ImageStack;
+import ij.process.ByteProcessor;
 import ij.process.FloatProcessor;
 import ij.process.ImageProcessor;
 
+import net.imglib2.algorithm.blocks.BlockSupplier;
+import net.imglib2.type.numeric.integer.UnsignedByteType;
+import org.janelia.alignment.intensity.Coefficients;
 import org.janelia.alignment.intensity.LinearIntensityMap;
 
 import net.imglib2.img.Img;
 import net.imglib2.img.array.ArrayImgs;
 import net.imglib2.img.imageplus.ImagePlusImgs;
 import net.imglib2.type.numeric.real.FloatType;
 
+import static org.janelia.alignment.intensity.FastLinearIntensityMap.linearIntensityMap;
+
 public class LinearIntensityMap8BitFilter
         extends IntensityMap8BitFilter {
 
@@ -59,40 +65,61 @@ public void after(final LinearIntensityMap8BitFilter otherFilter) {
     public void process(final ImageProcessor ip,
                         final double scale) {
 
-        final FloatProcessor as = new FloatProcessor(numberOfRegionColumns, numberOfRegionRows);
-        final FloatProcessor bs = new FloatProcessor(numberOfRegionColumns, numberOfRegionRows);
-
-        // this filter always produces 8-bit corrected output
-        final FloatProcessor fp = ip.convertToFloatProcessor();
-        fp.resetMinAndMax();
-        final double min = 0;
-        final double max = 255;
-
-        for (int i = 0; i < coefficients.length; ++i) {
-            final double[] ab = coefficients[i];
+        if (ip instanceof ByteProcessor) {
+            final byte[] pixels = (byte[]) ip.getPixels();
+            final Img<UnsignedByteType> img = ArrayImgs.unsignedBytes(pixels, ip.getWidth(), ip.getHeight());
 
             /* coefficients mapping into existing [min, max] */
-            as.setf(i, (float) ab[0]);
-            bs.setf(i, (float) ((max - min) * ab[1] + min - ab[0] * min));
+            final double min = 0;
+            final double max = 255;
+            final Coefficients.CoefficientFunction fn = (coefficientIndex, flattenedFieldIndex) -> {
+                final double[] ab = coefficients[flattenedFieldIndex];
+                return coefficientIndex == 0
+                        ? ab[0]
+                        : ((max - min) * ab[1] + min - ab[0] * min);
+            };
+
+            final Coefficients c = new Coefficients(fn, 2, numberOfRegionColumns, numberOfRegionRows );
+            BlockSupplier.of(img)
+                    .andThen(linearIntensityMap(c, img))
+                    .tile(256)
+                    .copy(img, pixels);
+        } else {
+            final FloatProcessor as = new FloatProcessor(numberOfRegionColumns, numberOfRegionRows);
+            final FloatProcessor bs = new FloatProcessor(numberOfRegionColumns, numberOfRegionRows);
+
+            // this filter always produces 8-bit corrected output
+            final FloatProcessor fp = ip.convertToFloatProcessor();
+            fp.resetMinAndMax();
+            final double min = 0;
+            final double max = 255;
+
+            for (int i = 0; i < coefficients.length; ++i) {
+                final double[] ab = coefficients[i];
+
+                /* coefficients mapping into existing [min, max] */
+                as.setf(i, (float) ab[0]);
+                bs.setf(i, (float) ((max - min) * ab[1] + min - ab[0] * min));
+            }
+            final ImageStack coefficientsStack = new ImageStack(numberOfRegionColumns, numberOfRegionRows);
+            coefficientsStack.addSlice(as);
+            coefficientsStack.addSlice(bs);
+
+            final LinearIntensityMap<FloatType> map =
+                    new LinearIntensityMap<FloatType>(
+                            ImagePlusImgs.from(
+                                    new ImagePlus("", coefficientsStack)
+                            )
+                    );
+
+            final long[] dims = new long[]{ip.getWidth(), ip.getHeight()};
+            final Img<FloatType> img = ArrayImgs.floats((float[]) fp.getPixels(), dims);
+
+            map.run(img);
+
+            // Need to reset intensity range back to full 8-bit before converting to byte processor!
+            fp.setMinAndMax(0, 255);
+            ip.setPixels(0, fp);
         }
-        final ImageStack coefficientsStack = new ImageStack(numberOfRegionColumns, numberOfRegionRows);
-        coefficientsStack.addSlice(as);
-        coefficientsStack.addSlice(bs);
-
-        final LinearIntensityMap<FloatType> map =
-                new LinearIntensityMap<FloatType>(
-                        ImagePlusImgs.from(
-                                new ImagePlus("", coefficientsStack)
-                        )
-                );
-
-        final long[] dims = new long[]{ip.getWidth(), ip.getHeight()};
-        final Img<FloatType> img = ArrayImgs.floats((float[]) fp.getPixels(), dims);
-
-        map.run(img);
-
-        // Need to reset intensity range back to full 8-bit before converting to byte processor!
-        fp.setMinAndMax(0, 255);
-        ip.setPixels(0, fp);
     }
 }

render-app/src/main/java/org/janelia/alignment/filter/QuadraticIntensityMap8BitFilter.java

@@ -2,15 +2,22 @@
 
 import ij.ImagePlus;
 import ij.ImageStack;
+import ij.process.ByteProcessor;
 import ij.process.FloatProcessor;
 import ij.process.ImageProcessor;
 
-import org.janelia.alignment.intensity.QuadraticIntensityMap;
+import net.imglib2.algorithm.blocks.BlockSupplier;
+import net.imglib2.type.numeric.integer.UnsignedByteType;
+import org.janelia.alignment.intensity.Coefficients;
+import org.janelia.alignment.intensity.Coefficients.CoefficientFunction;
 
 import net.imglib2.img.Img;
 import net.imglib2.img.array.ArrayImgs;
 import net.imglib2.img.imageplus.ImagePlusImgs;
 import net.imglib2.type.numeric.real.FloatType;
+import org.janelia.alignment.intensity.QuadraticIntensityMap;
+
+import static org.janelia.alignment.intensity.FastQuadraticIntensityMap.quadraticIntensityMap;
 
 public class QuadraticIntensityMap8BitFilter extends IntensityMap8BitFilter {
 
@@ -38,39 +45,67 @@ public QuadraticIntensityMap8BitFilter(final int numberOfRegionRows,
     public void process(final ImageProcessor ip,
                         final double scale) {
 
-        final FloatProcessor as = new FloatProcessor(numberOfRegionColumns, numberOfRegionRows);
-        final FloatProcessor bs = new FloatProcessor(numberOfRegionColumns, numberOfRegionRows);
-        final FloatProcessor cs = new FloatProcessor(numberOfRegionColumns, numberOfRegionRows);
-
-        final FloatProcessor fp = ip.convertToFloatProcessor();
-        fp.resetMinAndMax();
-        final double min = 0;
-        final double max = 255;
-        final double delta = max - min;
-
-        for (int i = 0; i < coefficients.length; ++i) {
-            final double[] abc = coefficients[i];
+        if (ip instanceof ByteProcessor) {
+            final byte[] pixels = (byte[]) ip.getPixels();
+            final Img<UnsignedByteType> img = ArrayImgs.unsignedBytes(pixels, ip.getWidth(), ip.getHeight());
 
             // mapping coefficients of polynomial on [0, 1] x [0, 1]
             // to coefficients of polynomial of the same shape on [min, max] x [min, max]
-            final double anew = abc[0] / delta;
-            as.setf(i, (float) anew);
-            bs.setf(i, (float) (min * anew * (min / delta - 2) + abc[1]));
-            cs.setf(i, (float) (min * (min * anew - abc[1]) + delta * abc[2] + min));
-        }
-        final ImageStack coefficientsStack = new ImageStack(numberOfRegionColumns, numberOfRegionRows);
-        coefficientsStack.addSlice(as);
-        coefficientsStack.addSlice(bs);
-        coefficientsStack.addSlice(cs);
+            final double min = 0;
+            final double max = 255;
+            final double delta = max - min;
+            final CoefficientFunction fn = (coefficientIndex, flattenedFieldIndex) -> {
+                final double[] abc = coefficients[flattenedFieldIndex];
+                final double anew = abc[0] / delta;
+                if (coefficientIndex == 0) {
+                    return anew;
+                } else if (coefficientIndex == 1) {
+                    return (min * anew * (min / delta - 2) + abc[1]);
+                } else {
+                    return (min * (min * anew - abc[1]) + delta * abc[2] + min);
+                }
+            };
 
-        final QuadraticIntensityMap<FloatType> map =
-                new QuadraticIntensityMap<FloatType>(ImagePlusImgs.from(new ImagePlus("", coefficientsStack)));
+            final Coefficients c = new Coefficients(fn, 3, numberOfRegionColumns, numberOfRegionRows );
+            BlockSupplier.of(img)
+                    .andThen(quadraticIntensityMap(c, img))
+                    .tile(256)
+                    .copy(img, pixels);
+        } else {
+            final FloatProcessor as = new FloatProcessor(numberOfRegionColumns, numberOfRegionRows);
+            final FloatProcessor bs = new FloatProcessor(numberOfRegionColumns, numberOfRegionRows);
+            final FloatProcessor cs = new FloatProcessor(numberOfRegionColumns, numberOfRegionRows);
 
-        final long[] dims = new long[]{ip.getWidth(), ip.getHeight()};
-        final Img<FloatType> img = ArrayImgs.floats((float[])fp.getPixels(), dims);
+            final FloatProcessor fp = ip.convertToFloatProcessor();
+            fp.resetMinAndMax();
+            final double min = 0;
+            final double max = 255;
+            final double delta = max - min;
 
-        map.run(img);
-        fp.setMinAndMax(0, 255);
-        ip.setPixels(0, fp);
+            for (int i = 0; i < coefficients.length; ++i) {
+                final double[] abc = coefficients[i];
+
+                // mapping coefficients of polynomial on [0, 1] x [0, 1]
+                // to coefficients of polynomial of the same shape on [min, max] x [min, max]
+                final double anew = abc[0] / delta;
+                as.setf(i, (float) anew);
+                bs.setf(i, (float) (min * anew * (min / delta - 2) + abc[1]));
+                cs.setf(i, (float) (min * (min * anew - abc[1]) + delta * abc[2] + min));
+            }
+            final ImageStack coefficientsStack = new ImageStack(numberOfRegionColumns, numberOfRegionRows);
+            coefficientsStack.addSlice(as);
+            coefficientsStack.addSlice(bs);
+            coefficientsStack.addSlice(cs);
+
+            final QuadraticIntensityMap<FloatType> map =
+                    new QuadraticIntensityMap<FloatType>(ImagePlusImgs.from(new ImagePlus("", coefficientsStack)));
+
+            final long[] dims = new long[]{ip.getWidth(), ip.getHeight()};
+            final Img<FloatType> img = ArrayImgs.floats((float[]) fp.getPixels(), dims);
+
+            map.run(img);
+            fp.setMinAndMax(0, 255);
+            ip.setPixels(0, fp);
+        }
     }
 }

render-app/src/main/java/org/janelia/alignment/intensity/Coefficients.java

@@ -0,0 +1,60 @@
+package org.janelia.alignment.intensity;
+
+import net.imglib2.util.IntervalIndexer;
+import net.imglib2.util.Intervals;
+
+import static net.imglib2.util.Util.safeInt;
+
+/**
+ * Holds flattened coefficient array and dimensions.
+ */
+public class Coefficients {
+
+    private final int[] size;
+    private final int[] strides;
+
+    /**
+     * {@code flattenedCoefficients[i]} holds the flattened array of the i-the coefficient.
+     * That is, for linear map {@code y=a*x+b}, {@code flattenedCoefficients[0]} holds all the {@code a}s and
+     * {@code flattenedCoefficients[1]} holds all the {@code b}s.
+     */
+    final float[][] flattenedCoefficients;
+
+    public Coefficients(
+            final double[][] coefficients,
+            final int... fieldDimensions) {
+        this((c, f) -> coefficients[f][c], coefficients[0].length, fieldDimensions);
+    }
+
+    @FunctionalInterface
+    public interface CoefficientFunction {
+        double apply(int coefficientIndex, int flattenedFieldIndex);
+    }
+
+    public Coefficients(
+            final CoefficientFunction coefficients,
+            final int numCoefficients,
+            final int... fieldDimensions) {
+        final int numElements = safeInt(Intervals.numElements(fieldDimensions));
+        size = fieldDimensions.clone();
+        strides = IntervalIndexer.createAllocationSteps(size);
+        flattenedCoefficients = new float[numCoefficients][numElements];
+        for (int j = 0; j < numCoefficients; ++j) {
+            for (int i = 0; i < numElements; ++i) {
+                flattenedCoefficients[j][i] = (float) coefficients.apply(j, i);
+            }
+        }
+    }
+
+    int size(final int d) {
+        return size[d];
+    }
+
+    int stride(final int d) {
+        return strides[d];
+    }
+
+    int numDimensions() {
+        return size.length;
+    }
+}