Skip to content

Commit bab64cf

Browse files
gselzergselzer
committed
WIP: port convolve, correlate, deconvolve Ops
TODO: test
1 parent 4dea941 commit bab64cf

26 files changed

Lines changed: 2688 additions & 432 deletions

‎src/main/java/net/imagej/ops/deconvolve/NonCirculantFirstGuess.java‎

Lines changed: 113 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,113 @@
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.deconvolve;
31+
32+
import java.util.function.BiFunction;
33+
34+
import net.imglib2.Dimensions;
35+
import net.imglib2.RandomAccessibleInterval;
36+
import net.imglib2.img.Img;
37+
import net.imglib2.type.numeric.ComplexType;
38+
import net.imglib2.type.numeric.RealType;
39+
import net.imglib2.view.Views;
40+
41+
import org.scijava.Priority;
42+
import org.scijava.ops.OpDependency;
43+
import org.scijava.ops.core.Op;
44+
import org.scijava.ops.core.computer.Computer;
45+
import org.scijava.ops.core.function.Function3;
46+
import org.scijava.param.Parameter;
47+
import org.scijava.plugin.Plugin;
48+
import org.scijava.struct.ItemIO;
49+
50+
/**
51+
* Calculate non-circulant first guess. This is used as part of the Boundary
52+
* condition handling scheme described here
53+
* http://bigwww.epfl.ch/deconvolution/challenge2013/index.html?p=doc_math_rl)
54+
*
55+
* @author Brian Northan
56+
* @param <I>
57+
* @param <O>
58+
* @param <K>
59+
* @param <C>
60+
*/
61+
62+
@Plugin(type = Op.class, name = "deconvolve.firstGuess", priority = Priority.LOW)
63+
@Parameter(key = "input")
64+
@Parameter(key = "outType")
65+
@Parameter(key = "k")
66+
@Parameter(key = "output", type = ItemIO.OUTPUT)
67+
public class NonCirculantFirstGuess<I extends RealType<I>, O extends RealType<O>, K extends RealType<K>, C extends ComplexType<C>>
68+
implements
69+
Function3<RandomAccessibleInterval<I>, O, Dimensions, RandomAccessibleInterval<O>>
70+
{
71+
72+
73+
@OpDependency(name = "create.img")
74+
BiFunction<Dimensions, O, Img<O>> create;
75+
76+
@OpDependency(name = "stats.sum")
77+
Computer<Iterable<I>, O> sum;
78+
79+
/**
80+
* k is the size of the measurement window. That is the size of the acquired
81+
* image before extension, k is required to calculate the non-circulant
82+
* normalization factor
83+
*/
84+
@Override
85+
public RandomAccessibleInterval<O> apply(RandomAccessibleInterval<I> in, final O outType, final Dimensions k) {
86+
87+
final Img<O> firstGuess = create.apply(in, outType);
88+
89+
// set first guess to be a constant = to the average value
90+
91+
// so first compute the sum...
92+
final O s = outType.createVariable();
93+
sum.compute(Views.iterable(in), s);
94+
95+
// then the number of pixels
96+
long numPixels = 1;
97+
98+
for (int d = 0; d < k.numDimensions(); d++) {
99+
numPixels = numPixels * k.dimension(d);
100+
}
101+
102+
// then the average value...
103+
final double average = s.getRealDouble() / (numPixels);
104+
105+
// set first guess as the average value computed above (TODO: use fill op)
106+
for (final O type : firstGuess) {
107+
type.setReal(average);
108+
}
109+
110+
return firstGuess;
111+
}
112+
113+
}

‎src/main/java/net/imagej/ops/deconvolve/NonCirculantNormalizationFactor.java‎

Lines changed: 239 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,239 @@
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.deconvolve;
31+
32+
import java.util.Iterator;
33+
import java.util.concurrent.ExecutorService;
34+
import java.util.function.BiFunction;
35+
import java.util.function.Function;
36+
37+
import net.imglib2.Cursor;
38+
import net.imglib2.Dimensions;
39+
import net.imglib2.FinalDimensions;
40+
import net.imglib2.FinalInterval;
41+
import net.imglib2.Interval;
42+
import net.imglib2.Point;
43+
import net.imglib2.RandomAccessibleInterval;
44+
import net.imglib2.img.Img;
45+
import net.imglib2.type.numeric.ComplexType;
46+
import net.imglib2.type.numeric.RealType;
47+
import net.imglib2.util.Util;
48+
import net.imglib2.view.Views;
49+
50+
import org.scijava.Priority;
51+
import org.scijava.ops.OpDependency;
52+
import org.scijava.ops.core.Op;
53+
import org.scijava.ops.core.computer.Computer3;
54+
import org.scijava.ops.core.computer.Computer7;
55+
import org.scijava.ops.core.inplace.Inplace6First;
56+
import org.scijava.param.Parameter;
57+
import org.scijava.plugin.Plugin;
58+
import org.scijava.struct.ItemIO;
59+
60+
/**
61+
* Calculate non-circulant normalization factor. This is used as part of the
62+
* Boundary condition handling scheme described here
63+
* http://bigwww.epfl.ch/deconvolution/challenge2013/index.html?p=doc_math_rl)
64+
*
65+
* @author Brian Northan
66+
* @param <I>
67+
* @param <O>
68+
* @param <K>
69+
* @param <C>
70+
*/
71+
72+
@Plugin(type = Op.class, name = "deconvolve.normalizationFactor",
73+
priority = Priority.LOW)
74+
@Parameter(key = "io", type = ItemIO.BOTH)
75+
@Parameter(key = "k")
76+
@Parameter(key = "l")
77+
@Parameter(key = "fftInput")
78+
@Parameter(key = "fftKernel")
79+
@Parameter(key = "executorService")
80+
public class NonCirculantNormalizationFactor<I extends RealType<I>, O extends RealType<O>, K extends RealType<K>, C extends ComplexType<C>>
81+
implements Inplace6First<RandomAccessibleInterval<O>, Dimensions, Dimensions, RandomAccessibleInterval<C>, RandomAccessibleInterval<C>, ExecutorService>
82+
{
83+
84+
/**
85+
* k is the size of the measurement window. That is the size of the acquired
86+
* image before extension, k is required to calculate the non-circulant
87+
* normalization factor
88+
*/
89+
private Dimensions k;
90+
91+
/**
92+
* l is the size of the psf, l is required to calculate the non-circulant
93+
* normalization factor
94+
*/
95+
private Dimensions l;
96+
97+
RandomAccessibleInterval<C> fftInput;
98+
99+
RandomAccessibleInterval<C> fftKernel;
100+
101+
// Normalization factor for edge handling (see
102+
// http://bigwww.epfl.ch/deconvolution/challenge2013/index.html?p=doc_math_rl)
103+
private Img<O> normalization = null;
104+
105+
@OpDependency(name = "create.img")
106+
private BiFunction<Dimensions, O, Img<O>> create;
107+
108+
@OpDependency(name = "copy.rai")
109+
private Function<RandomAccessibleInterval<O>, RandomAccessibleInterval<O>> copy;
110+
111+
@OpDependency(name = "filter.correlate")
112+
private Computer7<RandomAccessibleInterval<O>, RandomAccessibleInterval<K>, RandomAccessibleInterval<C>, RandomAccessibleInterval<C>, Boolean, Boolean, ExecutorService, RandomAccessibleInterval<O>> correlater;
113+
114+
// @OpDependency(name = "math.divide") TODO: allow the matcher to fix this
115+
private Computer3<Iterable<O>, Iterable<O>, Double, Iterable<O>> divide = (in1, in2, in3, out) -> {
116+
Iterator<O> itr1 = in1.iterator();
117+
Iterator<O> itr2 = in2.iterator();
118+
Iterator<O> itrout = out.iterator();
119+
120+
while(itr1.hasNext() && itr2.hasNext() && itrout.hasNext()) {
121+
Double val2 = itr2.next().getRealDouble();
122+
itrout.next().setReal(val2 == 0 ? in3 : itr1.next().getRealDouble() / val2);
123+
}
124+
};
125+
126+
/**
127+
* apply the normalization image needed for semi noncirculant model see
128+
* http://bigwww.epfl.ch/deconvolution/challenge2013/index.html?p=doc_math_rl
129+
*/
130+
@Override
131+
public void mutate(RandomAccessibleInterval<O> arg, final Dimensions k, final Dimensions l, final RandomAccessibleInterval<C> fftInput, final RandomAccessibleInterval<C> fftKernel, final ExecutorService es) {
132+
this.k = k;
133+
this.l = l;
134+
this.fftInput = fftInput;
135+
this.fftKernel = fftKernel;
136+
137+
// if the normalization image hasn't been computed yet, then compute it
138+
if (normalization == null) {
139+
this.createNormalizationImageSemiNonCirculant(arg, Util.getTypeFromInterval(arg), es);
140+
}
141+
142+
RandomAccessibleInterval<O> copyArg = copy.apply(arg);
143+
144+
// normalize for non-circulant deconvolution
145+
divide.accept(Views.iterable(copyArg), normalization, 0., Views.iterable(arg));
146+
}
147+
148+
protected void createNormalizationImageSemiNonCirculant(Interval fastFFTInterval, O type, ExecutorService es) {
149+
150+
// k is the window size (valid image region)
151+
final int length = k.numDimensions();
152+
153+
final long[] n = new long[length];
154+
final long[] nFFT = new long[length];
155+
156+
// n is the valid image size plus the extended region
157+
// also referred to as object space size
158+
for (int d = 0; d < length; d++) {
159+
n[d] = k.dimension(d) + l.dimension(d) - 1;
160+
}
161+
162+
// nFFT is the size of n after (potentially) extending further
163+
// to a fast FFT size
164+
for (int d = 0; d < length; d++) {
165+
nFFT[d] = fastFFTInterval.dimension(d);
166+
}
167+
168+
FinalDimensions fd = new FinalDimensions(nFFT);
169+
170+
// create the normalization image
171+
normalization = create.apply(fd, type);
172+
173+
// size of the measurement window
174+
final Point size = new Point(length);
175+
final long[] sizel = new long[length];
176+
177+
for (int d = 0; d < length; d++) {
178+
size.setPosition(k.dimension(d), d);
179+
sizel[d] = k.dimension(d);
180+
}
181+
182+
// starting point of the measurement window when it is centered in fft space
183+
final Point start = new Point(length);
184+
final long[] startl = new long[length];
185+
final long[] endl = new long[length];
186+
187+
for (int d = 0; d < length; d++) {
188+
start.setPosition((nFFT[d] - k.dimension(d)) / 2, d);
189+
startl[d] = (nFFT[d] - k.dimension(d)) / 2;
190+
endl[d] = startl[d] + sizel[d] - 1;
191+
}
192+
193+
// size of the object space
194+
final Point maskSize = new Point(length);
195+
final long[] maskSizel = new long[length];
196+
197+
for (int d = 0; d < length; d++) {
198+
maskSize.setPosition(Math.min(n[d], nFFT[d]), d);
199+
maskSizel[d] = Math.min(n[d], nFFT[d]);
200+
}
201+
202+
// starting point of the object space within the fft space
203+
final Point maskStart = new Point(length);
204+
final long[] maskStartl = new long[length];
205+
206+
for (int d = 0; d < length; d++) {
207+
maskStart.setPosition((Math.max(0, nFFT[d] - n[d]) / 2), d);
208+
maskStartl[d] = (Math.max(0, nFFT[d] - n[d]) / 2);
209+
}
210+
211+
final RandomAccessibleInterval<O> temp = Views.interval(normalization,
212+
new FinalInterval(startl, endl));
213+
final Cursor<O> normCursor = Views.iterable(temp).cursor();
214+
215+
// draw a cube the size of the measurement space
216+
while (normCursor.hasNext()) {
217+
normCursor.fwd();
218+
normCursor.get().setReal(1.0);
219+
}
220+
221+
final Img<O> tempImg = create.apply(fd, type);
222+
223+
// 3. correlate psf with the output of step 2.
224+
correlater.compute(normalization, null, fftInput, fftKernel, true, false, es, tempImg);
225+
226+
normalization = tempImg;
227+
228+
final Cursor<O> cursorN = normalization.cursor();
229+
230+
while (cursorN.hasNext()) {
231+
cursorN.fwd();
232+
233+
if (cursorN.get().getRealFloat() <= 1e-3f) {
234+
cursorN.get().setReal(1.0f);
235+
236+
}
237+
}
238+
}
239+
}

0 commit comments

Comments
 (0)