From d9ff43cd7b3e2aa86284417c09b41da217468999 Mon Sep 17 00:00:00 2001 From: Gabriel Selzer Date: Tue, 11 Feb 2025 10:41:13 -0600 Subject: [PATCH 1/2] Make DefaultMedian faster --- .../scijava/ops/image/stats/DefaultMedian.java | 15 ++++++++++----- .../scijava/ops/image/stats/DefaultQuantile.java | 3 ++- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/scijava-ops-image/src/main/java/org/scijava/ops/image/stats/DefaultMedian.java b/scijava-ops-image/src/main/java/org/scijava/ops/image/stats/DefaultMedian.java index 1fce33822..e40eacef3 100644 --- a/scijava-ops-image/src/main/java/org/scijava/ops/image/stats/DefaultMedian.java +++ b/scijava-ops-image/src/main/java/org/scijava/ops/image/stats/DefaultMedian.java @@ -29,12 +29,10 @@ package org.scijava.ops.image.stats; -import java.util.ArrayList; - import net.imglib2.type.numeric.RealType; +import org.scijava.collections.DoubleArray; import org.scijava.function.Computers; -import org.scijava.ops.spi.Op; import org.scijava.ops.spi.OpDependency; /** @@ -55,6 +53,11 @@ public class DefaultMedian, O extends RealType> @OpDependency(name = "stats.quantile") private Computers.Arity2, Double, O> quantileOp; + /** + * TODO: One/many of these lists could get really big over time. Maybe that will be a problem? + */ + private static final ThreadLocal BUFFER = ThreadLocal.withInitial(DoubleArray::new); + /** * TODO * @@ -63,9 +66,11 @@ public class DefaultMedian, O extends RealType> */ @Override public void compute(final Iterable input, final O median) { - final var statistics = new ArrayList(); + final var statistics = BUFFER.get(); + statistics.clear(); +// final var statistics = new ArrayList(); - input.forEach(i -> statistics.add(i.getRealDouble())); + input.forEach(i -> statistics.addValue(i.getRealDouble())); final var k = statistics.size() / 2; var result = DefaultQuantile.select(statistics, k); diff --git a/scijava-ops-image/src/main/java/org/scijava/ops/image/stats/DefaultQuantile.java b/scijava-ops-image/src/main/java/org/scijava/ops/image/stats/DefaultQuantile.java index 5259b482e..257841580 100644 --- a/scijava-ops-image/src/main/java/org/scijava/ops/image/stats/DefaultQuantile.java +++ b/scijava-ops-image/src/main/java/org/scijava/ops/image/stats/DefaultQuantile.java @@ -33,6 +33,7 @@ import org.scijava.function.Computers; import java.util.ArrayList; +import java.util.List; import static java.util.Collections.swap; @@ -83,7 +84,7 @@ public void compute(final Iterable input, final N quantile, * This an all-in-one method version of your basic quick select algorithm. *

*/ - static double select(final ArrayList array, final int k) { + static double select(final List array, final int k) { var left = 0; var right = array.size() - 1; From ccd6ee6e4b90e65b24f5c93ee1310492e19ed1ab Mon Sep 17 00:00:00 2001 From: Edward Evans Date: Fri, 21 Feb 2025 08:00:58 -0600 Subject: [PATCH 2/2] WIP add sliding window median filter and test Working Op, however the output seems to differ from the original ImageJ implementation. Further testing is required. Additionally, it needs to be bench marked. --- .../ops/image/stats/SlidingWindowMedian.java | 175 ++++++++++++++++++ .../image/filter/NonLinearFiltersTest.java | 21 +++ 2 files changed, 196 insertions(+) create mode 100644 scijava-ops-image/src/main/java/org/scijava/ops/image/stats/SlidingWindowMedian.java diff --git a/scijava-ops-image/src/main/java/org/scijava/ops/image/stats/SlidingWindowMedian.java b/scijava-ops-image/src/main/java/org/scijava/ops/image/stats/SlidingWindowMedian.java new file mode 100644 index 000000000..c6150408f --- /dev/null +++ b/scijava-ops-image/src/main/java/org/scijava/ops/image/stats/SlidingWindowMedian.java @@ -0,0 +1,175 @@ +/*- + * #%L + * Image processing operations for SciJava Ops. + * %% + * Copyright (C) 2014 - 2025 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.image.stats; + +import net.imglib2.RandomAccessible; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.algorithm.neighborhood.RectangleShape; +import net.imglib2.outofbounds.OutOfBoundsMirrorFactory; +import net.imglib2.outofbounds.OutOfBoundsMirrorFactory.Boundary; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.view.Views; + +import org.scijava.function.Computers; +import org.scijava.function.Container; + +/** + * @author Edward Evans + * @author Gabriel Selzer + * @implNote op names='filter.median', priority='100.' + */ + +public class SlidingWindowMedian + implements Computers.Arity2, Integer, RandomAccessibleInterval> { + + /** + * Sliding window median filter + * + * @param input The input image. + * @param span The span value used to create a square kernel. The kernel + * shape is (size, size) where size is (span * 2 + 1). For example, + * a span value of 1 produces a (3, 3) kernel. + * @param output The output image/container. + */ + @Override + public void compute(RandomAccessible input, + final Integer span, @Container RandomAccessibleInterval output) { + // prep input image + var oobf = new OutOfBoundsMirrorFactory>(Boundary.SINGLE); + var extended = Views.extend((RandomAccessibleInterval) input, oobf); + + // create sliding window buffer + final int size = span * 2 + 1; + float[][] window = new float[size][size]; + float[] buffer = new float[size * size]; + int bufLen = 0; + int medianIndex = 0; + + // get random access + var raIn = extended.randomAccess(); + var raOut = output.randomAccess(); + + // slide window across the image + // here 'i' and 'j' are 'x' and 'y' on the output image + int oY = 0; + int wX = 0; + int wY = 0; + long[] dims = output.dimensionsAsLongArray(); + for (int i = 0; i < dims[0]; i++) { + // read a full window at the top of the image + // here 'k' and 'l' are 'x' and 'y' of the window + for (int l = -span + oY; l < size + (-span + oY); l++) { + // set window y-axis (row) + raIn.setPosition(l, 1); + for (int k = -span + i; k < size + (-span + i); k++) { + // set window x-axis (col) + raIn.setPosition(k, 0); + var v = raIn.get().get(); + window[wY][wX] = v; + insertSorted(buffer, v, bufLen++); + wX++; + } + wX = 0; + wY++; + } + wY = 0; + // compute median and store the output + medianIndex = buffer.length / 2; + raOut.setPositionAndGet(i, oY).set(buffer[medianIndex]); + // advance the window column wise (down, +y) remove old row, fill new row + oY++; + for (int j = size; oY < dims[1]; j++) { + var row = j % size; + // remove old data from the row from window and buffer + for (int w = 0; w < window[row].length; w++) { + removeSorted(buffer, window[row][w]); + bufLen--; + } + // set random access to new row y-axis pos + raIn.setPosition(j, 1); + for (int e = -span + i; e < size + (-span + i); e++) { + // set random acces to new row x-axis pos + raIn.setPosition(e, 0); + var v = raIn.get().get(); + window[row][wX] = v; + insertSorted(buffer, v, bufLen++); + wX++; + } + wX = 0; + // compute new median and store in output + medianIndex = buffer.length / 2; + raOut.setPositionAndGet(i, oY).set(buffer[medianIndex]); + oY++; + } + oY = 0; + bufLen = 0; + } + } + + public static void insertSorted(float[] arr, float value, int len) { + // find the sorted position the value should be inserted + int insertPos = 0; + while (insertPos < len && arr[insertPos] < value && arr[insertPos] != 0.0f) { + insertPos++; + } + + // If we're not inserting at the end, shift elements to make space + if (insertPos < len) { + // Find the last non-zero element + //int lastIndex = len - 1; + //while (lastIndex > 0 && arr[lastIndex] == 0.0f) { + // lastIndex--; + //} + + // Shift elements to the right + for (int i = len; i >= insertPos + 1; i--) { + arr[i] = arr[i - 1]; + } + } + + // Insert the value + arr[insertPos] = value; + } + + public static void removeSorted(float[] arr, float value) { + // Find the position where value should be inserted + int insertPos = 0; + while (insertPos < arr.length && arr[insertPos] != value) { + insertPos++; + } + + // Shift elements to the right + for (int i = insertPos; i < arr.length - 1; i++) { + arr[i] = arr[i + 1]; + } + } + +} diff --git a/scijava-ops-image/src/test/java/org/scijava/ops/image/filter/NonLinearFiltersTest.java b/scijava-ops-image/src/test/java/org/scijava/ops/image/filter/NonLinearFiltersTest.java index 164dcdf55..e5bcbc450 100644 --- a/scijava-ops-image/src/test/java/org/scijava/ops/image/filter/NonLinearFiltersTest.java +++ b/scijava-ops-image/src/test/java/org/scijava/ops/image/filter/NonLinearFiltersTest.java @@ -34,8 +34,11 @@ import java.util.ArrayList; import java.util.Collections; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.type.numeric.real.FloatType; import org.scijava.ops.image.AbstractOpTest; import org.scijava.ops.image.filter.sigma.DefaultSigmaFilter; +import org.scijava.ops.image.stats.SlidingWindowMedian; import org.scijava.ops.image.util.TestImgGeneration; import net.imglib2.algorithm.neighborhood.Neighborhood; import net.imglib2.algorithm.neighborhood.RectangleShape; @@ -131,6 +134,24 @@ public void testMeanFilterNullable() { ops.op("filter.mean").input(in, shape).output(out).compute(); } + /** + * @see NeighborhoodFilters#defaultMedian(Computers.Arity1, Neighborhood, + * Object) + */ + @Test + public void testFoo() { + var foo = new SlidingWindowMedian(); + var image = ops.op("convert.float32").input(in).apply(); + var oobf = new OutOfBoundsMirrorFactory>(Boundary.SINGLE); + var extend = Views.extend((RandomAccessibleInterval) image, oobf); + var output = ops.op("create.img").input(image).apply(); + foo.compute(extend, 1, (RandomAccessibleInterval) output); + + var expected = ops.op("filter.median").input(image, shape, oobf).apply(); + var raOut = ((RandomAccessibleInterval) output).randomAccess(); + var raExp = ((RandomAccessibleInterval) expected).randomAccess(); + } + /** * @see NeighborhoodFilters#defaultMedian(Computers.Arity1, Neighborhood, * Object)