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; 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)