From 7764e63239609b675dc61654b51956fdc5b962cc Mon Sep 17 00:00:00 2001 From: Gabriel Selzer Date: Mon, 11 Jun 2018 13:58:55 -0500 Subject: [PATCH 1/5] Add DefaultSmooth and DefaultSharpen Ops Adapted from the implementations created by Barry Dezonia --- .../imagej/ops/filter/FilterNamespace.java | 43 ++++- .../ops/filter/sharpen/DefaultSharpen.java | 164 ++++++++++++++++++ .../ops/filter/smooth/DefaultSmooth.java | 163 +++++++++++++++++ src/main/templates/net/imagej/ops/Ops.list | 2 + 4 files changed, 367 insertions(+), 5 deletions(-) create mode 100644 src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java create mode 100644 src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java diff --git a/src/main/java/net/imagej/ops/filter/FilterNamespace.java b/src/main/java/net/imagej/ops/filter/FilterNamespace.java index 7b49291dd1..0219bf92aa 100644 --- a/src/main/java/net/imagej/ops/filter/FilterNamespace.java +++ b/src/main/java/net/imagej/ops/filter/FilterNamespace.java @@ -6,13 +6,13 @@ * %% * 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 @@ -131,9 +131,16 @@ public , O extends RealType> O addPoissonNoise( return result; } - // -- bilateral -- - - /** Executes the "bilateral" filter on the given arguments. */ + /** + * Executes a bilateral filter on the given arguments. + * + * @param in + * @param out + * @param sigmaR + * @param sigmaS + * @param radius + * @return + */ @OpMethod(op = net.imagej.ops.filter.bilateral.DefaultBilateral.class) public , O extends RealType> RandomAccessibleInterval bilateral(final RandomAccessibleInterval out, @@ -1299,6 +1306,19 @@ public > RandomAccessibleInterval partialDerivative( return result; } + // -- Sharpen + + /** Executes the "sharpen" filter operation on the given arguments. */ + @OpMethod(op = net.imagej.ops.filter.sharpen.DefaultSharpen.class) + public > RandomAccessibleInterval sharpen( + final RandomAccessibleInterval in) + { + @SuppressWarnings("unchecked") + final RandomAccessibleInterval result = + (RandomAccessibleInterval) ops().run(Ops.Filter.Sharpen.class, in); + return result; + } + /** Executes the "sigma" filter operation on the given arguments. */ @OpMethod(op = net.imagej.ops.filter.sigma.DefaultSigmaFilter.class) public > IterableInterval sigma( @@ -1325,6 +1345,19 @@ public > IterableInterval sigma( return result; } + // -- Smooth + + /** Executes the "sharpen" filter operation on the given arguments. */ + @OpMethod(op = net.imagej.ops.filter.smooth.DefaultSmooth.class) + public > RandomAccessibleInterval smooth( + final RandomAccessibleInterval in) + { + @SuppressWarnings("unchecked") + final RandomAccessibleInterval result = + (RandomAccessibleInterval) ops().run(Ops.Filter.Smooth.class, in); + return result; + } + // -- Sobel @OpMethod(op = net.imagej.ops.filter.sobel.SobelRAI.class) diff --git a/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java b/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java new file mode 100644 index 0000000000..8b09e0837d --- /dev/null +++ b/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java @@ -0,0 +1,164 @@ +/* + * #%L + * ImageJ software for multidimensional image processing and analysis. + * %% + * Copyright (C) 2014 - 2018 ImageJ 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 net.imagej.ops.filter.sharpen; + +import net.imagej.Extents; +import net.imagej.Position; +import net.imagej.ops.Ops; +import net.imagej.ops.special.function.AbstractUnaryFunctionOp; +import net.imglib2.Cursor; +import net.imglib2.FinalInterval; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessible; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.algorithm.neighborhood.Neighborhood; +import net.imglib2.algorithm.neighborhood.RectangleNeighborhood; +import net.imglib2.algorithm.neighborhood.RectangleNeighborhoodFactory; +import net.imglib2.algorithm.neighborhood.RectangleShape.NeighborhoodsAccessible; +import net.imglib2.type.numeric.RealType; +import net.imglib2.util.Util; +import net.imglib2.view.Views; + +import org.scijava.plugin.Plugin; + +/** + * Op rendition of the SharpenDataValues plugin written by Barry Dezonia + * + * @author Gabe Selzer + */ +@Plugin(type = Ops.Filter.Sharpen.class) +public class DefaultSharpen> extends + AbstractUnaryFunctionOp, RandomAccessibleInterval> + implements Ops.Filter.Sharpen +{ + + final double[] kernel = { -1, -1, -1, -1, 12, -1, -1, -1, -1 }; + double scale; + + @Override + public RandomAccessibleInterval calculate( + final RandomAccessibleInterval input) + { + final RandomAccessibleInterval output = ops().copy().rai(input); + + if(input.numDimensions() > 2) { + throw new IllegalArgumentException("Input has too few dimensions! Only 2+ dimensional images allowed!"); + } + + final long[] planeDims = new long[input.numDimensions() - 2]; + for (int i = 0; i < planeDims.length; i++) + planeDims[i] = input.dimension(i + 2); + final Extents extents = new Extents(planeDims); + final Position planePos = extents.createPosition(); + if (planeDims.length == 0) { + computePlanar(planePos, input, output); + } + else { + while (planePos.hasNext()) { + planePos.fwd(); + computePlanar(planePos, input, output); + } + + } + return output; + } + + private void computePlanar(final Position planePos, + final RandomAccessibleInterval input, + final RandomAccessibleInterval output) + { + // TODO can we just set scale to 4? + scale = 0; + for (final double d : kernel) + scale += d; + + final T type = Util.getTypeFromInterval(input); + + final long[] imageDims = new long[input.numDimensions()]; + input.dimensions(imageDims); + + // create all objects needed for NeighborhoodsAccessible + RandomAccessibleInterval slicedInput = ops().copy().rai(input); + for (int i = planePos.numDimensions() - 1; i >= 0; i--) { + slicedInput = Views.hyperSlice(slicedInput, input.numDimensions() - 1 - i, + planePos.getLongPosition(i)); + } + + final RandomAccessible refactoredInput = Views.extendMirrorSingle( + slicedInput); + final RectangleNeighborhoodFactory factory = RectangleNeighborhood + .factory(); + final FinalInterval neighborhoodSpan = new FinalInterval(new long[] { -1, + -1 }, new long[] { 1, 1 }); + + final NeighborhoodsAccessible neighborhoods = + new NeighborhoodsAccessible<>(refactoredInput, neighborhoodSpan, factory); + + // create cursors and random accesses for loop. + final Cursor cursor = Views.iterable(input).localizingCursor(); + final RandomAccess outputRA = output.randomAccess(); + for (int i = 0; i < planePos.numDimensions(); i++) { + outputRA.setPosition(planePos.getLongPosition(i), i + 2); + } + final RandomAccess> neighborhoodsRA = neighborhoods + .randomAccess(); + + int algorithmIndex = 0; + double sum; + final double[] n = new double[9]; + while (cursor.hasNext()) { + cursor.fwd(); + if (cursor.getLongPosition(0) == 14 && cursor.getLongPosition(1) == 0) + System.out.println("Hit 14"); + neighborhoodsRA.setPosition(cursor); + final Neighborhood current = neighborhoodsRA.get(); + final Cursor neighborhoodCursor = current.cursor(); + + algorithmIndex = 0; + sum = 0; + while (algorithmIndex < n.length) { + neighborhoodCursor.fwd(); + n[algorithmIndex++] = neighborhoodCursor.get().getRealDouble(); + } + + for (int i = 0; i < kernel.length; i++) { + sum += kernel[i] * n[i]; + } + + double value = sum / scale; + + outputRA.setPosition(cursor.getLongPosition(0), 0); + outputRA.setPosition(cursor.getLongPosition(1), 1); + if (value > type.getMaxValue()) value = type.getMaxValue(); + if (value < type.getMinValue()) value = type.getMinValue(); + outputRA.get().setReal(value); + } + } +} diff --git a/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java b/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java new file mode 100644 index 0000000000..6bba85b36e --- /dev/null +++ b/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java @@ -0,0 +1,163 @@ +/* + * #%L + * ImageJ software for multidimensional image processing and analysis. + * %% + * Copyright (C) 2014 - 2018 ImageJ 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 net.imagej.ops.filter.smooth; + +import net.imagej.Extents; +import net.imagej.Position; +import net.imagej.ops.Ops; +import net.imagej.ops.special.function.AbstractUnaryFunctionOp; +import net.imglib2.Cursor; +import net.imglib2.FinalInterval; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessible; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.algorithm.neighborhood.Neighborhood; +import net.imglib2.algorithm.neighborhood.RectangleNeighborhood; +import net.imglib2.algorithm.neighborhood.RectangleNeighborhoodFactory; +import net.imglib2.algorithm.neighborhood.RectangleShape.NeighborhoodsAccessible; +import net.imglib2.type.numeric.RealType; +import net.imglib2.util.Util; +import net.imglib2.view.Views; + +import org.scijava.plugin.Plugin; + +/** + * Op rendition of the SmoothDataValues plugin written by Barry Dezonia + * + * @author Gabe Selzer + */ +@Plugin(type = Ops.Filter.Smooth.class) +public class DefaultSmooth> extends + AbstractUnaryFunctionOp, RandomAccessibleInterval> + implements Ops.Filter.Smooth +{ + + final double[] kernel = { 1, 1, 1, 1, 1, 1, 1, 1, 1 }; + double scale; + + @Override + public RandomAccessibleInterval calculate( + final RandomAccessibleInterval input) + { + final RandomAccessibleInterval output = ops().copy().rai(input); + + if(input.numDimensions() > 2) { + throw new IllegalArgumentException("Input has too few dimensions! Only 2+ dimensional images allowed!"); + } + + final long[] planeDims = new long[input.numDimensions() - 2]; + for (int i = 0; i < planeDims.length; i++) + planeDims[i] = input.dimension(i + 2); + final Extents extents = new Extents(planeDims); + final Position planePos = extents.createPosition(); + if (planeDims.length == 0) { + computePlanar(planePos, input, output); + } + else { + while (planePos.hasNext()) { + planePos.fwd(); + computePlanar(planePos, input, output); + } + + } + return output; + } + + private void computePlanar(final Position planePos, + final RandomAccessibleInterval input, + final RandomAccessibleInterval output) + { + // TODO can we just set scale to 4? + scale = 0; + for (final double d : kernel) + scale += d; + + final T type = Util.getTypeFromInterval(input); + + final long[] imageDims = new long[input.numDimensions()]; + input.dimensions(imageDims); + + // create all objects needed for NeighborhoodsAccessible + RandomAccessibleInterval slicedInput = ops().copy().rai(input); + for (int i = planePos.numDimensions() - 1; i >= 0; i--) { + slicedInput = Views.hyperSlice(slicedInput, input.numDimensions() - 1 - i, + planePos.getLongPosition(i)); + } + + final RandomAccessible refactoredInput = Views.extendMirrorSingle( + slicedInput); + final RectangleNeighborhoodFactory factory = RectangleNeighborhood + .factory(); + final FinalInterval neighborhoodSpan = new FinalInterval(new long[] { -1, + -1 }, new long[] { 1, 1 }); + + final NeighborhoodsAccessible neighborhoods = + new NeighborhoodsAccessible<>(refactoredInput, neighborhoodSpan, factory); + + // create cursors and random accesses for loop. + final Cursor cursor = Views.iterable(input).localizingCursor(); + final RandomAccess outputRA = output.randomAccess(); + for (int i = 0; i < planePos.numDimensions(); i++) { + outputRA.setPosition(planePos.getLongPosition(i), i + 2); + } + final RandomAccess> neighborhoodsRA = neighborhoods + .randomAccess(); + + int algorithmIndex = 0; + double sum; + final double[] n = new double[9]; + while (cursor.hasNext()) { + + cursor.fwd(); + neighborhoodsRA.setPosition(cursor); + final Neighborhood current = neighborhoodsRA.get(); + final Cursor neighborhoodCursor = current.cursor(); + + algorithmIndex = 0; + sum = 0; + while (algorithmIndex < n.length) { + neighborhoodCursor.fwd(); + n[algorithmIndex++] = neighborhoodCursor.get().getRealDouble(); + } + + for (int i = 0; i < kernel.length; i++) { + sum += kernel[i] * n[i]; + } + + double value = sum / scale; + + outputRA.setPosition(cursor.getLongPosition(0), 0); + outputRA.setPosition(cursor.getLongPosition(1), 1); + if (value > type.getMaxValue()) value = type.getMaxValue(); + if (value < type.getMinValue()) value = type.getMinValue(); + outputRA.get().setReal(value); + } + } +} diff --git a/src/main/templates/net/imagej/ops/Ops.list b/src/main/templates/net/imagej/ops/Ops.list index 15a6095fc4..8a7cb12b1f 100644 --- a/src/main/templates/net/imagej/ops/Ops.list +++ b/src/main/templates/net/imagej/ops/Ops.list @@ -123,7 +123,9 @@ namespaces = ``` [name: "paddingIntervalCentered", iface: "PaddingIntervalCentered"], [name: "paddingIntervalOrigin", iface: "PaddingIntervalOrigin"], [name: "tubeness", iface: "Tubeness"], + [name: "sharpen", iface: "Sharpen"], [name: "sigma", iface: "Sigma", aliases: ["sigmaFilter", "filterSigma"]], + [name: "smooth", iface: "Smooth"], [name: "sobel", iface: "Sobel"], [name: "variance", iface: "Variance", aliases: ["varianceFilter", "filterVariance", "var", "varFilter", "filterVar"]], [name: "frangiVesselness", iface: "FrangiVesselness"], From f5cff8c6d3938e3c8e4238ef5d216fdef2721e85 Mon Sep 17 00:00:00 2001 From: Gabriel Selzer Date: Mon, 14 Jan 2019 16:26:40 -0600 Subject: [PATCH 2/5] Preserve differentiation in summation from plugin --- .../imagej/ops/filter/sharpen/DefaultSharpen.java | 12 ++++++++++-- .../net/imagej/ops/filter/smooth/DefaultSmooth.java | 10 +++++++++- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java b/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java index 8b09e0837d..8d28e1c998 100644 --- a/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java +++ b/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java @@ -42,6 +42,7 @@ import net.imglib2.algorithm.neighborhood.RectangleNeighborhood; import net.imglib2.algorithm.neighborhood.RectangleNeighborhoodFactory; import net.imglib2.algorithm.neighborhood.RectangleShape.NeighborhoodsAccessible; +import net.imglib2.type.numeric.IntegerType; import net.imglib2.type.numeric.RealType; import net.imglib2.util.Util; import net.imglib2.view.Views; @@ -151,8 +152,15 @@ private void computePlanar(final Position planePos, for (int i = 0; i < kernel.length; i++) { sum += kernel[i] * n[i]; } - - double value = sum / scale; + + //find the value for the output + double value; + if(type instanceof IntegerType) { + value = (sum + scale / 2) / scale; + } + else { + value = sum / scale; + } outputRA.setPosition(cursor.getLongPosition(0), 0); outputRA.setPosition(cursor.getLongPosition(1), 1); diff --git a/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java b/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java index 6bba85b36e..d5169cae2a 100644 --- a/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java +++ b/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java @@ -42,6 +42,7 @@ import net.imglib2.algorithm.neighborhood.RectangleNeighborhood; import net.imglib2.algorithm.neighborhood.RectangleNeighborhoodFactory; import net.imglib2.algorithm.neighborhood.RectangleShape.NeighborhoodsAccessible; +import net.imglib2.type.numeric.IntegerType; import net.imglib2.type.numeric.RealType; import net.imglib2.util.Util; import net.imglib2.view.Views; @@ -151,7 +152,14 @@ private void computePlanar(final Position planePos, sum += kernel[i] * n[i]; } - double value = sum / scale; + //find the value for the output + double value; + if(type instanceof IntegerType) { + value = (sum + scale / 2) / scale; + } + else { + value = sum / scale; + } outputRA.setPosition(cursor.getLongPosition(0), 0); outputRA.setPosition(cursor.getLongPosition(1), 1); From b5b2822da3120724b4c8a0f8f622077065de5609 Mon Sep 17 00:00:00 2001 From: Gabriel Selzer Date: Fri, 18 Jan 2019 09:46:52 -0600 Subject: [PATCH 3/5] Rely on helper Ops to decrease amount of code --- .../ops/filter/sharpen/DefaultSharpen.java | 157 +++++++++-------- .../ops/filter/smooth/DefaultSmooth.java | 158 +++++++++--------- 2 files changed, 153 insertions(+), 162 deletions(-) diff --git a/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java b/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java index 8d28e1c998..0c34a54622 100644 --- a/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java +++ b/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java @@ -32,18 +32,19 @@ import net.imagej.Extents; import net.imagej.Position; import net.imagej.ops.Ops; +import net.imagej.ops.special.computer.Computers; +import net.imagej.ops.special.computer.UnaryComputerOp; import net.imagej.ops.special.function.AbstractUnaryFunctionOp; -import net.imglib2.Cursor; -import net.imglib2.FinalInterval; -import net.imglib2.RandomAccess; +import net.imagej.ops.special.function.UnaryFunctionOp; +import net.imagej.ops.special.inplace.Inplaces; +import net.imagej.ops.special.inplace.UnaryInplaceOp; +import net.imglib2.IterableInterval; import net.imglib2.RandomAccessible; import net.imglib2.RandomAccessibleInterval; -import net.imglib2.algorithm.neighborhood.Neighborhood; -import net.imglib2.algorithm.neighborhood.RectangleNeighborhood; -import net.imglib2.algorithm.neighborhood.RectangleNeighborhoodFactory; -import net.imglib2.algorithm.neighborhood.RectangleShape.NeighborhoodsAccessible; import net.imglib2.type.numeric.IntegerType; import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.integer.ByteType; +import net.imglib2.type.numeric.real.DoubleType; import net.imglib2.util.Util; import net.imglib2.view.Views; @@ -60,113 +61,107 @@ public class DefaultSharpen> extends implements Ops.Filter.Sharpen { - final double[] kernel = { -1, -1, -1, -1, 12, -1, -1, -1, -1 }; - double scale; + final double[][] kernel = { { -1, -1, -1 }, { -1, 12, -1 }, { -1, -1, -1 } }; + RandomAccessibleInterval kernelRAI; + + /** + * The sum of all of the kernel values + */ + double scale = 4; + + UnaryFunctionOp> kernelCreator; + UnaryComputerOp, RandomAccessibleInterval> convolveOp; + UnaryInplaceOp, IterableInterval> addConstOp; + UnaryInplaceOp, IterableInterval> divConstOp; + UnaryComputerOp clipTypesOp; + UnaryComputerOp, IterableInterval> convertOp; + + @SuppressWarnings({ "unchecked", "rawtypes" }) + @Override + public void initialize() { + T inType = Util.getTypeFromInterval(in()); + + // convolution kernel + kernelRAI = ops().create().kernel(kernel, new ByteType()); + + IterableInterval dummyDoubleType = ops().create().img(in(), new DoubleType()); + + convolveOp = (UnaryComputerOp) Computers.unary(ops(), + Ops.Filter.Convolve.class, RandomAccessibleInterval.class, + RandomAccessible.class, kernelRAI); + addConstOp = (UnaryInplaceOp) Inplaces.unary(ops(), Ops.Math.Add.class, + dummyDoubleType, new DoubleType(scale / 2)); + divConstOp = (UnaryInplaceOp) Inplaces.unary(ops(), Ops.Math.Divide.class, + dummyDoubleType, new DoubleType(scale)); + clipTypesOp = (UnaryComputerOp) Computers.unary(ops(), + Ops.Convert.Clip.class, new DoubleType(), inType); + convertOp = Computers.unary(ops(), + Ops.Convert.ImageType.class, Views.iterable(in()), dummyDoubleType, clipTypesOp); + + } + + @SuppressWarnings("unchecked") @Override public RandomAccessibleInterval calculate( final RandomAccessibleInterval input) { - final RandomAccessibleInterval output = ops().copy().rai(input); + // intermediate image to hold the convolution data. Depending on the input + // type we have to be able create an intermediate of wide enough type. + final RandomAccessibleInterval intermediate = ops().create().img(in(), new DoubleType()); + // output image + final RandomAccessibleInterval output = (RandomAccessibleInterval) ops().create().img(in()); if(input.numDimensions() > 2) { throw new IllegalArgumentException("Input has too few dimensions! Only 2+ dimensional images allowed!"); } + // compute the sharpening on 2D slices of the image final long[] planeDims = new long[input.numDimensions() - 2]; for (int i = 0; i < planeDims.length; i++) planeDims[i] = input.dimension(i + 2); final Extents extents = new Extents(planeDims); final Position planePos = extents.createPosition(); if (planeDims.length == 0) { - computePlanar(planePos, input, output); + computePlanar(planePos, input, intermediate); } else { while (planePos.hasNext()) { planePos.fwd(); - computePlanar(planePos, input, output); + computePlanar(planePos, input, intermediate); } } + + T inputType = Util.getTypeFromInterval(input); + IterableInterval iterableIntermediate = Views.iterable( + intermediate); + + // divide by the scale, if integerType input also add scale / 2. + if (inputType instanceof IntegerType) addConstOp.mutate( + iterableIntermediate); + divConstOp.mutate(iterableIntermediate); + + //convert the result back to the input type. + convertOp.compute(iterableIntermediate, Views.iterable(output)); + return output; } private void computePlanar(final Position planePos, final RandomAccessibleInterval input, - final RandomAccessibleInterval output) + final RandomAccessibleInterval intermediate) { - // TODO can we just set scale to 4? - scale = 0; - for (final double d : kernel) - scale += d; - - final T type = Util.getTypeFromInterval(input); - - final long[] imageDims = new long[input.numDimensions()]; - input.dimensions(imageDims); - - // create all objects needed for NeighborhoodsAccessible - RandomAccessibleInterval slicedInput = ops().copy().rai(input); + // create 2D slice in the case of a (N>2)-dimensional image. + RandomAccessible slicedInput = Views.extendMirrorSingle(input); + RandomAccessibleInterval slicedIntermediate = intermediate; for (int i = planePos.numDimensions() - 1; i >= 0; i--) { slicedInput = Views.hyperSlice(slicedInput, input.numDimensions() - 1 - i, planePos.getLongPosition(i)); + slicedIntermediate = Views.hyperSlice(slicedIntermediate, intermediate + .numDimensions() - 1 - i, planePos.getLongPosition(i)); } - final RandomAccessible refactoredInput = Views.extendMirrorSingle( - slicedInput); - final RectangleNeighborhoodFactory factory = RectangleNeighborhood - .factory(); - final FinalInterval neighborhoodSpan = new FinalInterval(new long[] { -1, - -1 }, new long[] { 1, 1 }); - - final NeighborhoodsAccessible neighborhoods = - new NeighborhoodsAccessible<>(refactoredInput, neighborhoodSpan, factory); - - // create cursors and random accesses for loop. - final Cursor cursor = Views.iterable(input).localizingCursor(); - final RandomAccess outputRA = output.randomAccess(); - for (int i = 0; i < planePos.numDimensions(); i++) { - outputRA.setPosition(planePos.getLongPosition(i), i + 2); - } - final RandomAccess> neighborhoodsRA = neighborhoods - .randomAccess(); - - int algorithmIndex = 0; - double sum; - final double[] n = new double[9]; - while (cursor.hasNext()) { - cursor.fwd(); - if (cursor.getLongPosition(0) == 14 && cursor.getLongPosition(1) == 0) - System.out.println("Hit 14"); - neighborhoodsRA.setPosition(cursor); - final Neighborhood current = neighborhoodsRA.get(); - final Cursor neighborhoodCursor = current.cursor(); - - algorithmIndex = 0; - sum = 0; - while (algorithmIndex < n.length) { - neighborhoodCursor.fwd(); - n[algorithmIndex++] = neighborhoodCursor.get().getRealDouble(); - } - - for (int i = 0; i < kernel.length; i++) { - sum += kernel[i] * n[i]; - } - - //find the value for the output - double value; - if(type instanceof IntegerType) { - value = (sum + scale / 2) / scale; - } - else { - value = sum / scale; - } - - outputRA.setPosition(cursor.getLongPosition(0), 0); - outputRA.setPosition(cursor.getLongPosition(1), 1); - if (value > type.getMaxValue()) value = type.getMaxValue(); - if (value < type.getMinValue()) value = type.getMinValue(); - outputRA.get().setReal(value); - } + convolveOp.compute(slicedInput, slicedIntermediate); } } diff --git a/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java b/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java index d5169cae2a..ad1c6f13ba 100644 --- a/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java +++ b/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java @@ -32,18 +32,19 @@ import net.imagej.Extents; import net.imagej.Position; import net.imagej.ops.Ops; +import net.imagej.ops.special.computer.Computers; +import net.imagej.ops.special.computer.UnaryComputerOp; import net.imagej.ops.special.function.AbstractUnaryFunctionOp; -import net.imglib2.Cursor; -import net.imglib2.FinalInterval; -import net.imglib2.RandomAccess; +import net.imagej.ops.special.function.UnaryFunctionOp; +import net.imagej.ops.special.inplace.Inplaces; +import net.imagej.ops.special.inplace.UnaryInplaceOp; +import net.imglib2.IterableInterval; import net.imglib2.RandomAccessible; import net.imglib2.RandomAccessibleInterval; -import net.imglib2.algorithm.neighborhood.Neighborhood; -import net.imglib2.algorithm.neighborhood.RectangleNeighborhood; -import net.imglib2.algorithm.neighborhood.RectangleNeighborhoodFactory; -import net.imglib2.algorithm.neighborhood.RectangleShape.NeighborhoodsAccessible; import net.imglib2.type.numeric.IntegerType; import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.integer.ByteType; +import net.imglib2.type.numeric.real.DoubleType; import net.imglib2.util.Util; import net.imglib2.view.Views; @@ -60,112 +61,107 @@ public class DefaultSmooth> extends implements Ops.Filter.Smooth { - final double[] kernel = { 1, 1, 1, 1, 1, 1, 1, 1, 1 }; - double scale; + final double[][] kernel = { { 1, 1, 1 }, { 1, 1, 1 }, { 1, 1, 1 } }; + RandomAccessibleInterval kernelRAI; + + /** + * The sum of all of the kernel values + */ + double scale = 9; + + UnaryFunctionOp> kernelCreator; + UnaryComputerOp, RandomAccessibleInterval> convolveOp; + UnaryInplaceOp, IterableInterval> addConstOp; + UnaryInplaceOp, IterableInterval> divConstOp; + UnaryComputerOp clipTypesOp; + UnaryComputerOp, IterableInterval> convertOp; + + @SuppressWarnings({ "unchecked", "rawtypes" }) + @Override + public void initialize() { + T inType = Util.getTypeFromInterval(in()); + + // convolution kernel + kernelRAI = ops().create().kernel(kernel, new ByteType()); + + IterableInterval dummyDoubleType = ops().create().img(in(), new DoubleType()); + + convolveOp = (UnaryComputerOp) Computers.unary(ops(), + Ops.Filter.Convolve.class, RandomAccessibleInterval.class, + RandomAccessible.class, kernelRAI); + addConstOp = (UnaryInplaceOp) Inplaces.unary(ops(), Ops.Math.Add.class, + dummyDoubleType, new DoubleType(scale / 2)); + divConstOp = (UnaryInplaceOp) Inplaces.unary(ops(), Ops.Math.Divide.class, + dummyDoubleType, new DoubleType(scale)); + clipTypesOp = (UnaryComputerOp) Computers.unary(ops(), + Ops.Convert.Clip.class, new DoubleType(), inType); + convertOp = Computers.unary(ops(), + Ops.Convert.ImageType.class, Views.iterable(in()), dummyDoubleType, clipTypesOp); + + } + + @SuppressWarnings("unchecked") @Override public RandomAccessibleInterval calculate( final RandomAccessibleInterval input) { - final RandomAccessibleInterval output = ops().copy().rai(input); + // intermediate image to hold the convolution data. Depending on the input + // type we have to be able create an intermediate of wide enough type. + final RandomAccessibleInterval intermediate = ops().create().img(in(), new DoubleType()); + // output image + final RandomAccessibleInterval output = (RandomAccessibleInterval) ops().create().img(in()); if(input.numDimensions() > 2) { throw new IllegalArgumentException("Input has too few dimensions! Only 2+ dimensional images allowed!"); } + // compute the smoothing on 2D slices of the image final long[] planeDims = new long[input.numDimensions() - 2]; for (int i = 0; i < planeDims.length; i++) planeDims[i] = input.dimension(i + 2); final Extents extents = new Extents(planeDims); final Position planePos = extents.createPosition(); if (planeDims.length == 0) { - computePlanar(planePos, input, output); + computePlanar(planePos, input, intermediate); } else { while (planePos.hasNext()) { planePos.fwd(); - computePlanar(planePos, input, output); + computePlanar(planePos, input, intermediate); } } + + T inputType = Util.getTypeFromInterval(input); + IterableInterval iterableIntermediate = Views.iterable( + intermediate); + + // divide by the scale, if integerType input also add scale / 2. + if (inputType instanceof IntegerType) addConstOp.mutate( + iterableIntermediate); + divConstOp.mutate(iterableIntermediate); + + //convert the result back to the input type. + convertOp.compute(iterableIntermediate, Views.iterable(output)); + return output; } private void computePlanar(final Position planePos, final RandomAccessibleInterval input, - final RandomAccessibleInterval output) + final RandomAccessibleInterval intermediate) { - // TODO can we just set scale to 4? - scale = 0; - for (final double d : kernel) - scale += d; - - final T type = Util.getTypeFromInterval(input); - - final long[] imageDims = new long[input.numDimensions()]; - input.dimensions(imageDims); - - // create all objects needed for NeighborhoodsAccessible - RandomAccessibleInterval slicedInput = ops().copy().rai(input); + // create 2D slice in the case of a (N>2)-dimensional image. + RandomAccessible slicedInput = Views.extendMirrorSingle(input); + RandomAccessibleInterval slicedIntermediate = intermediate; for (int i = planePos.numDimensions() - 1; i >= 0; i--) { slicedInput = Views.hyperSlice(slicedInput, input.numDimensions() - 1 - i, planePos.getLongPosition(i)); + slicedIntermediate = Views.hyperSlice(slicedIntermediate, intermediate + .numDimensions() - 1 - i, planePos.getLongPosition(i)); } - final RandomAccessible refactoredInput = Views.extendMirrorSingle( - slicedInput); - final RectangleNeighborhoodFactory factory = RectangleNeighborhood - .factory(); - final FinalInterval neighborhoodSpan = new FinalInterval(new long[] { -1, - -1 }, new long[] { 1, 1 }); - - final NeighborhoodsAccessible neighborhoods = - new NeighborhoodsAccessible<>(refactoredInput, neighborhoodSpan, factory); - - // create cursors and random accesses for loop. - final Cursor cursor = Views.iterable(input).localizingCursor(); - final RandomAccess outputRA = output.randomAccess(); - for (int i = 0; i < planePos.numDimensions(); i++) { - outputRA.setPosition(planePos.getLongPosition(i), i + 2); - } - final RandomAccess> neighborhoodsRA = neighborhoods - .randomAccess(); - - int algorithmIndex = 0; - double sum; - final double[] n = new double[9]; - while (cursor.hasNext()) { - - cursor.fwd(); - neighborhoodsRA.setPosition(cursor); - final Neighborhood current = neighborhoodsRA.get(); - final Cursor neighborhoodCursor = current.cursor(); - - algorithmIndex = 0; - sum = 0; - while (algorithmIndex < n.length) { - neighborhoodCursor.fwd(); - n[algorithmIndex++] = neighborhoodCursor.get().getRealDouble(); - } - - for (int i = 0; i < kernel.length; i++) { - sum += kernel[i] * n[i]; - } - - //find the value for the output - double value; - if(type instanceof IntegerType) { - value = (sum + scale / 2) / scale; - } - else { - value = sum / scale; - } - - outputRA.setPosition(cursor.getLongPosition(0), 0); - outputRA.setPosition(cursor.getLongPosition(1), 1); - if (value > type.getMaxValue()) value = type.getMaxValue(); - if (value < type.getMinValue()) value = type.getMinValue(); - outputRA.get().setReal(value); - } + convolveOp.compute(slicedInput, slicedIntermediate); } -} +} \ No newline at end of file From 65b9bb5e290e55c674dfbf73c2353584065fbfef Mon Sep 17 00:00:00 2001 From: Gabriel Selzer Date: Fri, 18 Jan 2019 13:47:18 -0600 Subject: [PATCH 4/5] Create regression tests --- .../ops/filter/sharpen/SharpenTest.java | 68 +++++++++++++++++++ .../imagej/ops/filter/smooth/SmoothTest.java | 68 +++++++++++++++++++ 2 files changed, 136 insertions(+) create mode 100644 src/test/java/net/imagej/ops/filter/sharpen/SharpenTest.java create mode 100644 src/test/java/net/imagej/ops/filter/smooth/SmoothTest.java diff --git a/src/test/java/net/imagej/ops/filter/sharpen/SharpenTest.java b/src/test/java/net/imagej/ops/filter/sharpen/SharpenTest.java new file mode 100644 index 0000000000..c9bf6a8cd3 --- /dev/null +++ b/src/test/java/net/imagej/ops/filter/sharpen/SharpenTest.java @@ -0,0 +1,68 @@ +/* + * #%L + * ImageJ software for multidimensional image processing and analysis. + * %% + * Copyright (C) 2014 - 2018 ImageJ 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 net.imagej.ops.filter.sharpen; + +import java.util.stream.IntStream; + +import net.imagej.ops.AbstractOpTest; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.type.numeric.integer.IntType; + +import org.junit.Test; + +/** + * Tests {@link DefaultSharpen}. + * + * @author Gabe Selzer + */ +public class SharpenTest extends AbstractOpTest { + + @Test + public void testRegression() { + int[] range = IntStream.rangeClosed(0, 24).toArray(); + int[] regressionData = { -9, -6, -5, -4, -2, 4, 7, 8, 9, 11, 9, 12, 13, 14, + 16, 14, 17, 18, 19, 21, 27, 29, 30, 31, 34 }; + Img img = ArrayImgs.ints(range, 5, 5); + + Img actual = (Img) ops.filter().sharpen(img); + Img expected = ArrayImgs.ints(regressionData, 5l, 5l); + + assertIterationsEqual(expected, actual); + } + + @Test (expected = NegativeArraySizeException.class) + public void testTooFewDimensions() { + int[] range = IntStream.rangeClosed(1, 7).toArray(); + Img img = ArrayImgs.ints(range, 7); + + Img output = (Img) ops.filter().sharpen(img); + } +} diff --git a/src/test/java/net/imagej/ops/filter/smooth/SmoothTest.java b/src/test/java/net/imagej/ops/filter/smooth/SmoothTest.java new file mode 100644 index 0000000000..c871ca4921 --- /dev/null +++ b/src/test/java/net/imagej/ops/filter/smooth/SmoothTest.java @@ -0,0 +1,68 @@ +/* + * #%L + * ImageJ software for multidimensional image processing and analysis. + * %% + * Copyright (C) 2014 - 2018 ImageJ 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 net.imagej.ops.filter.smooth; + +import java.util.stream.IntStream; + +import net.imagej.ops.AbstractOpTest; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.type.numeric.integer.IntType; + +import org.junit.Test; + +/** + * Tests {@link DefaultSmooth}. + * + * @author Gabe Selzer + */ +public class SmoothTest extends AbstractOpTest { + + @Test + public void testRegression() { + int[] range = IntStream.rangeClosed(0, 24).toArray(); + int[] regressionData = { 5, 5, 6, 7, 7, 6, 7, 8, 9, 9, 11, 12, 13, 14, 14, + 16, 17, 18, 19, 19, 18, 18, 19, 20, 21 }; + Img img = ArrayImgs.ints(range, 5, 5); + + Img actual = (Img) ops.filter().smooth(img); + Img expected = ArrayImgs.ints(regressionData, 5l, 5l); + + assertIterationsEqual(expected, actual); + } + + @Test(expected = NegativeArraySizeException.class) + public void testTooFewDimensions() { + int[] range = IntStream.rangeClosed(1, 7).toArray(); + Img img = ArrayImgs.ints(range, 7); + + Img output = (Img) ops.filter().smooth(img); + } +} From dd53936a078e131dbe54256e697acb44c34d3371 Mon Sep 17 00:00:00 2001 From: Gabriel Selzer Date: Tue, 19 May 2020 14:33:08 -0500 Subject: [PATCH 5/5] Fix dimensionality check The dimensionality check seemed to fail the Op when the input had MORE than two dimensions. It should fail the Op when the input has LESS than two dimensions --- src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java | 2 +- src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java b/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java index 0c34a54622..5fda418b4f 100644 --- a/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java +++ b/src/main/java/net/imagej/ops/filter/sharpen/DefaultSharpen.java @@ -112,7 +112,7 @@ public RandomAccessibleInterval calculate( // output image final RandomAccessibleInterval output = (RandomAccessibleInterval) ops().create().img(in()); - if(input.numDimensions() > 2) { + if(input.numDimensions() < 2) { throw new IllegalArgumentException("Input has too few dimensions! Only 2+ dimensional images allowed!"); } diff --git a/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java b/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java index ad1c6f13ba..40b3a2acbf 100644 --- a/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java +++ b/src/main/java/net/imagej/ops/filter/smooth/DefaultSmooth.java @@ -112,7 +112,7 @@ public RandomAccessibleInterval calculate( // output image final RandomAccessibleInterval output = (RandomAccessibleInterval) ops().create().img(in()); - if(input.numDimensions() > 2) { + if(input.numDimensions() < 2) { throw new IllegalArgumentException("Input has too few dimensions! Only 2+ dimensional images allowed!"); }