Skip to content

Commit c462a93

Browse files
gselzerwiedenm
gselzer
authored and
wiedenm
committed
Port DefaultBilateral filter as proof of concept
1 parent 56fc47a commit c462a93

2 files changed

Lines changed: 319 additions & 0 deletions

File tree

‎src/main/java/net/imagej/ops/filter/bilateral/DefaultBilateral.java‎

Lines changed: 139 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,139 @@
1+
/*
2+
* #%L
3+
* ImageJ software for multidimensional image processing and analysis.
4+
* %%
5+
* Copyright (C) 2014 - 2018 ImageJ developers.
6+
* %%
7+
* Redistribution and use in source and binary forms, with or without
8+
* modification, are permitted provided that the following conditions are met:
9+
*
10+
* 1. Redistributions of source code must retain the above copyright notice,
11+
* this list of conditions and the following disclaimer.
12+
* 2. Redistributions in binary form must reproduce the above copyright notice,
13+
* this list of conditions and the following disclaimer in the documentation
14+
* and/or other materials provided with the distribution.
15+
*
16+
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17+
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18+
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19+
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
20+
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21+
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
22+
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
23+
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24+
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25+
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
26+
* POSSIBILITY OF SUCH DAMAGE.
27+
* #L%
28+
*/
29+
30+
package net.imagej.ops.filter.bilateral;
31+
32+
import net.imglib2.Cursor;
33+
import net.imglib2.FinalInterval;
34+
import net.imglib2.Interval;
35+
import net.imglib2.RandomAccess;
36+
import net.imglib2.RandomAccessibleInterval;
37+
import net.imglib2.algorithm.neighborhood.Neighborhood;
38+
import net.imglib2.algorithm.neighborhood.RectangleNeighborhood;
39+
import net.imglib2.algorithm.neighborhood.RectangleNeighborhoodFactory;
40+
import net.imglib2.type.numeric.RealType;
41+
import net.imglib2.view.Views;
42+
43+
import org.scijava.ops.core.Op;
44+
import org.scijava.ops.core.computer.Computer4;
45+
import org.scijava.param.Parameter;
46+
import org.scijava.plugin.Plugin;
47+
import org.scijava.struct.ItemIO;
48+
49+
/**
50+
* Performs a bilateral filter on an image.
51+
*
52+
* @author Gabe Selzer
53+
* @param <I>
54+
* @param <O>
55+
*/
56+
57+
@Plugin(type = Op.class, name = "filter.bilateral")
58+
@Parameter(key = "inputRAI", description = "the input data")
59+
@Parameter(key = "sigmaR", description = "range smoothing param, larger sigma means larger effect of intensity differences.")
60+
@Parameter(key = "sigmaS", description = "spatial smoothing param, larger sigma means smoother image.")
61+
@Parameter(key = "radius", description = "defines size of the square of pixels considered at each iteration.")
62+
@Parameter(key = "outputRAI", type = ItemIO.BOTH)
63+
public class DefaultBilateral<I extends RealType<I>, O extends RealType<O>>
64+
implements Computer4<RandomAccessibleInterval<I>, Double, Double, Integer, RandomAccessibleInterval<O>> {
65+
66+
public final static int MIN_DIMS = 2;
67+
68+
public final static int MAX_DIMS = 2;
69+
70+
private static double gauss(final double x, final double sigma) {
71+
final double mu = 0.0;
72+
return 1 / (sigma * Math.sqrt(2 * Math.PI)) * Math.exp(-0.5 * (x - mu) * (x - mu) / (sigma * sigma));
73+
}
74+
75+
private double getDistance(long[] x, long[] y) {
76+
double distance = 0;
77+
78+
for (int i = 0; i < x.length; i++) {
79+
double separation = x[i] - y[i];
80+
if (separation != 0) {
81+
distance += separation * separation;
82+
}
83+
84+
}
85+
86+
return Math.sqrt(distance);
87+
}
88+
89+
@Override
90+
public void compute(final RandomAccessibleInterval<I> input, final Double sigmaR, final Double sigmaS,
91+
final Integer radius, final RandomAccessibleInterval<O> output) {
92+
93+
final long[] size = new long[input.numDimensions()];
94+
input.dimensions(size);
95+
96+
final RandomAccess<O> outputRA = output.randomAccess();
97+
final Cursor<I> inputCursor = Views.iterable(input).localizingCursor();
98+
final long[] currentPos = new long[input.numDimensions()];
99+
final long[] neighborhoodPos = new long[input.numDimensions()];
100+
final long[] neighborhoodMin = new long[input.numDimensions()];
101+
final long[] neighborhoodMax = new long[input.numDimensions()];
102+
Neighborhood<I> neighborhood;
103+
Cursor<I> neighborhoodCursor;
104+
final RectangleNeighborhoodFactory<I> fac = RectangleNeighborhood.factory();
105+
while (inputCursor.hasNext()) {
106+
inputCursor.fwd();
107+
inputCursor.localize(currentPos);
108+
double distance;
109+
inputCursor.localize(neighborhoodMin);
110+
inputCursor.localize(neighborhoodMax);
111+
neighborhoodMin[0] = Math.max(0, neighborhoodMin[0] - radius);
112+
neighborhoodMin[1] = Math.max(0, neighborhoodMin[1] - radius);
113+
neighborhoodMax[0] = Math.min(input.max(0), neighborhoodMax[0] + radius);
114+
neighborhoodMax[1] = Math.min(input.max(1), neighborhoodMax[1] + radius);
115+
final Interval interval = new FinalInterval(neighborhoodMin, neighborhoodMax);
116+
neighborhood = fac.create(currentPos, neighborhoodMin, neighborhoodMax, interval, input.randomAccess());
117+
neighborhoodCursor = neighborhood.localizingCursor();
118+
double weight, v = 0.0;
119+
double w = 0.0;
120+
do {
121+
neighborhoodCursor.fwd();
122+
neighborhoodCursor.localize(neighborhoodPos);
123+
distance = getDistance(currentPos, neighborhoodPos);
124+
weight = gauss(distance, sigmaS);// spatial kernel
125+
126+
distance = Math.abs(inputCursor.get().getRealDouble() - neighborhoodCursor.get().getRealDouble());// intensity
127+
// difference
128+
weight *= gauss(distance, sigmaR);// range kernel, then exponent addition
129+
130+
v += weight * neighborhoodCursor.get().getRealDouble();
131+
w += weight;
132+
} while (neighborhoodCursor.hasNext());
133+
outputRA.setPosition(currentPos);
134+
outputRA.get().setReal(v / w);
135+
}
136+
137+
}
138+
139+
}

‎src/test/java/net/imagej/ops/filter/bilateral/DefaultBilateralTest.java‎

Lines changed: 180 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,180 @@
1+
/*-
2+
* #%L
3+
* ImageJ software for multidimensional image processing and analysis.
4+
* %%
5+
* Copyright (C) 2014 - 2018 ImageJ developers.
6+
* %%
7+
* Redistribution and use in source and binary forms, with or without
8+
* modification, are permitted provided that the following conditions are met:
9+
*
10+
* 1. Redistributions of source code must retain the above copyright notice,
11+
* this list of conditions and the following disclaimer.
12+
* 2. Redistributions in binary form must reproduce the above copyright notice,
13+
* this list of conditions and the following disclaimer in the documentation
14+
* and/or other materials provided with the distribution.
15+
*
16+
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17+
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18+
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19+
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
20+
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21+
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
22+
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
23+
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24+
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25+
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
26+
* POSSIBILITY OF SUCH DAMAGE.
27+
* #L%
28+
*/
29+
package net.imagej.ops.filter.bilateral;
30+
31+
import static org.junit.Assert.assertEquals;
32+
33+
import net.imagej.ops.AbstractOpTest;
34+
import net.imglib2.Cursor;
35+
import net.imglib2.img.Img;
36+
import net.imglib2.img.array.ArrayImgs;
37+
import net.imglib2.type.numeric.integer.ByteType;
38+
39+
import org.junit.Test;
40+
41+
public class DefaultBilateralTest extends AbstractOpTest {
42+
43+
@Test
44+
public void testBigImage() {
45+
final byte[] data = { 7, 8, 9, 1, 2, 3, 7, 9, 8, 1, 3, 2, 8, 7, 9, 2, 1, 3, 8, 9, 7, 2, 3, 1, 9, 7, 8, 3, 1, 2,
46+
9, 8, 7, 3, 2, 1 };
47+
final Img<ByteType> in = ArrayImgs.bytes(data, 6, 6);
48+
final Img<ByteType> out = generateByteArrayTestImg(false, 6, 6);
49+
50+
ops.run("filter.bilateral", in, 15.0, 5.0, 2, out);
51+
52+
final byte[] expected = { 8, 7, 6, 4, 3, 2, 8, 7, 6, 4, 3, 2, 8, 7, 6, 4, 3, 2, 8, 7, 6, 4, 3, 2, 8, 7, 6, 4, 3,
53+
2, 8, 7, 6, 4, 3, 2 };
54+
55+
Cursor<ByteType> cout = out.cursor();
56+
for (int i = 0; i < expected.length; i++) {
57+
assertEquals(cout.next().get(), expected[i]);
58+
}
59+
}
60+
61+
@Test
62+
public void testMath() {
63+
final byte[] data = { 7, 4, 9, 1 };
64+
final Img<ByteType> in = ArrayImgs.bytes(data, 2, 2);
65+
final Img<ByteType> out = generateByteArrayTestImg(false, 2, 2);
66+
67+
ops.run("filter.bilateral", in, 15.0, 5.0, 1, out);
68+
69+
Cursor<ByteType> cout = out.cursor();
70+
final byte[] expected = { 5, 5, 5, 5 };
71+
int counter = 0;
72+
while (cout.hasNext()) {
73+
byte actual = cout.next().get();
74+
assertEquals(expected[counter++], actual);
75+
}
76+
}
77+
78+
@Test
79+
public void testArrayToCellImg() {
80+
81+
final byte[] data = { 7, 8, 9, 1, 2, 3, 7, 9, 8, 1, 3, 2, 8, 7, 9, 2, 1, 3, 8, 9, 7, 2, 3, 1, 9, 7, 8, 3, 1, 2,
82+
9, 8, 7, 3, 2, 1 };
83+
final Img<ByteType> in = ArrayImgs.bytes(data, 6, 6);
84+
final Img<ByteType> out = generateByteArrayTestImg(false, 6, 6);
85+
final Img<ByteType> cellOut = generateByteTestCellImg(false, 6, 6);
86+
87+
ops.run("filter.bilateral", in, 15.0, 5.0, 2, out);
88+
ops.run("filter.bilateral", in, 15.0, 5.0, 2, cellOut);
89+
90+
Cursor<ByteType> cout = out.cursor();
91+
Cursor<ByteType> cCellOut = cellOut.cursor();
92+
while (cout.hasNext()) {
93+
byte expected = cout.next().get();
94+
byte actual = cCellOut.next().get();
95+
assertEquals(expected, actual);
96+
}
97+
}
98+
99+
// @Test
100+
// public void testGaussianVsBilateral() {
101+
// final byte[] data = { 7, 8, 9, 1, 2, 3, 7, 9, 8, 1, 3, 2, 8, 7, 9, 2, 1, 3, 8, 9, 7, 2, 3, 1, 9, 7, 8, 3, 1, 2,
102+
// 9, 8, 7, 3, 2, 1 };
103+
// final Img<ByteType> in = ArrayImgs.bytes(data, 6, 6);
104+
// final Img<ByteType> gaussOut = generateByteArrayTestImg(false, 6, 6);
105+
// final Img<ByteType> bilateralOut = generateByteTestCellImg(false, 6, 6);
106+
//
107+
// ops.run("filter.bilateral", bilateralOut, in, 15, 5, 2);
108+
// final double sigma = 5;
109+
// ops.run(GaussRAISingleSigma.class, gaussOut, in, sigma);
110+
// assertEquals(areCongruent(gaussOut, bilateralOut, 0), false);
111+
// }
112+
113+
@Test
114+
public void testZeroes() {
115+
final byte[] data = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
116+
0, 0, 0, 0, 0, 0 };
117+
final Img<ByteType> in = ArrayImgs.bytes(data, 6, 6);
118+
final Img<ByteType> out = generateByteArrayTestImg(false, 6, 6);
119+
120+
ops.run("filter.bilateral", in, 15.0, 5.0, 2, out);
121+
122+
Cursor<ByteType> cout = out.cursor();
123+
while (cout.hasNext()) {
124+
byte expected = cout.next().get();
125+
assertEquals(expected, 0);
126+
}
127+
}
128+
129+
@Test
130+
public void testNegatives() {
131+
final byte[] data = { -7, -8, -9, -1, -2, -3, -7, -9, -8, -1, -3, -2, -8, -7, -9, -2, -1, -3, -8, -9, -7, -2,
132+
-3, -1, -9, -7, -8, -3, -1, -2, -9, -8, -7, -3, -2, -1 };
133+
final Img<ByteType> in = ArrayImgs.bytes(data, 6, 6);
134+
final Img<ByteType> out = generateByteArrayTestImg(false, 6, 6);
135+
136+
ops.run("filter.bilateral", in, 15.0, 5.0, 2, out);
137+
138+
final byte[] expected = { -8, -7, -6, -4, -3, -2, -8, -7, -6, -4, -3, -2, -8, -7, -6, -4, -3, -2, -8, -7, -6,
139+
-4, -3, -2, -8, -7, -6, -4, -3, -2, -8, -7, -6, -4, -3, -2 };
140+
141+
Cursor<ByteType> cout = out.cursor();
142+
for (int i = 0; i < expected.length; i++) {
143+
assertEquals(cout.next().get(), expected[i]);
144+
}
145+
}
146+
147+
//TODO the next few tests do not run because we do not have any form of Contingent.
148+
149+
// @Test(expected = IllegalArgumentException.class)
150+
// public void testTooManyDimensions() {
151+
// final byte[] data = { 2, 2, 2, 2, 2, 2, 2, 2 };
152+
// final Img<ByteType> in = ArrayImgs.bytes(data, 2, 2);
153+
// final Img<ByteType> out = generateByteArrayTestImg(false, 2, 2, 2);
154+
//
155+
// ops.run("filter.bilateral", in, 15.0, 5.0, 2, out);
156+
//
157+
// final byte[] expected = { 2, 2, 2, 2, 2, 2, 2, 2 };
158+
//
159+
// Cursor<ByteType> cout = out.cursor();
160+
// for (int i = 0; i < expected.length; i++) {
161+
// assertEquals(cout.next().get(), expected[i]);
162+
// }
163+
// }
164+
//
165+
// @Test(expected = IllegalArgumentException.class)
166+
// public void testMismatchedDimensions() {
167+
// final byte[] data = { 1, 1, 1, 1, 1, 1 };
168+
// final Img<ByteType> in = ArrayImgs.bytes(data, 2, 3);
169+
// final Img<ByteType> out = generateByteArrayTestImg(false, 3, 2);
170+
//
171+
// ops.run("filter.bilateral", in, 15.0, 5.0, 2, out);
172+
//
173+
// final byte[] expected = { 1, 1, 1, 1, 1, 1 };
174+
// Cursor<ByteType> cout = out.cursor();
175+
// for (int i = 0; i < expected.length; i++) {
176+
// assertEquals(cout.next().get(), expected[i]);
177+
// }
178+
// }
179+
180+
}

0 commit comments

Comments
 (0)