Skip to content

Commit 1aeff37

Browse files
ctruedenctrueden
authored
Merge pull request #145 from scijava/scijava-ops-image/decon-fixes-and-qol
Fix deconvolution namespace and quality test
2 parents 179a9dc + 4ee41c0 commit 1aeff37

File tree

9 files changed

+393
-138
lines changed

9 files changed

+393
-138
lines changed

ScriptingInFiji.md renamed to docs/ops/doc/ScriptingInFiji.md

File renamed without changes.

‎docs/ops/doc/examples/deconvolution.rst‎

Lines changed: 149 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,149 @@
1+
=============================================
2+
Richardson-Lucy Total Variation Deconvolution
3+
=============================================
4+
5+
In this example we will use SciJava Ops to perform Richardson-Lucy (RL) deconvolution on a 3D dataset (X, Y, Z) of
6+
a HeLa cell nucleus stained with DAPI (4′,6-diamidino-2-phenylindole) and imaged on an epifluorescent microscope at 100x.
7+
The SciJava Ops framework currently supports the standard RL algorithm as well as the Richardson-Lucy Total Variation (RLTV)
8+
algorithm, which utilizes a regularization factor to limit the noise amplified by the RL algorithm :sup:`1`. Typically,
9+
the RLTV algorithm returns improved axial and lateral resolution when compared to RL.
10+
11+
You can download the 3D HeLa cell nuclus dataset `here`_.
12+
13+
.. figure:: https://media.imagej.net/scijava-ops/1.0.0/rltv_example_1.gif
14+
15+
Results of RLTV deconvolution on the sample data.
16+
17+
RLTV parameter descriptions
18+
===========================
19+
20+
The table below contains the necessary parameter values needed for the ``kernelDiffraction`` Op to create the synthetic
21+
point spread function (PSF) using the Gibson-Lanni model :sup:`2` for the sample HeLa cell nucleus dataset.
22+
23+
+--------------------------------------+-------+
24+
| Parameter | Value |
25+
+======================================+=======+
26+
| Iterations | 15 |
27+
+--------------------------------------+-------+
28+
| Numerical aperature | 1.45 |
29+
+--------------------------------------+-------+
30+
| Emission wavelength (nm) | 457 |
31+
+--------------------------------------+-------+
32+
| Refractive index (Immersion) | 1.5 |
33+
+--------------------------------------+-------+
34+
| Refractive index (Sample) | 1.4 |
35+
+--------------------------------------+-------+
36+
| Lateral resolution (μm/pixel) | 0.065 |
37+
+--------------------------------------+-------+
38+
| Axial resolution (μm/pixel) | 0.1 |
39+
+--------------------------------------+-------+
40+
| Particle/sample position (μm/pixel) | 0.0 |
41+
+--------------------------------------+-------+
42+
| Regularization factor | 0.002 |
43+
+--------------------------------------+-------+
44+
45+
SciJava Ops via Fiji's scripting engine with `script parameters`_:
46+
47+
.. tabs::
48+
49+
.. code-tab:: scijava-groovy
50+
51+
#@ OpEnvironment ops
52+
#@ ImgPlus img
53+
#@ Integer iterations(label="Iterations", value=30)
54+
#@ Float numericalAperture(label="Numerical Aperture", style="format:0.00", min=0.00, value=1.45)
55+
#@ Integer wavelength(label="Emission Wavelength (nm)", value=550)
56+
#@ Float riImmersion(label="Refractive Index (immersion)", style="format:0.00", min=0.00, value=1.5)
57+
#@ Float riSample(label="Refractive Index (sample)", style="format:0.00", min=0.00, value=1.4)
58+
#@ Float lateral_res(label="Lateral resolution (μm/pixel)", style="format:0.0000", min=0.0000, value=0.065)
59+
#@ Float axial_res(label="Axial resolution (μm/pixel)", style="format:0.0000", min=0.0000, value=0.1)
60+
#@ Float pZ(label="Particle/sample Position (μm)", style="format:0.0000", min=0.0000, value=0)
61+
#@ Float regularizationFactor(label="Regularization factor", style="format:0.00000", min=0.00000, value=0.002)
62+
#@output ImgPlus psf
63+
#@output ImgPlus result
64+
65+
import net.imglib2.FinalDimensions
66+
import net.imglib2.type.numeric.real.FloatType
67+
import net.imglib2.type.numeric.complex.ComplexFloatType
68+
69+
// convert input image to float
70+
img_float = ops.op("create.img").arity2().input(img, new FloatType()).apply()
71+
ops.op("convert.float32").arity1().input(img).output(img_float).compute()
72+
73+
// use image dimensions for PSF size
74+
psf_size = new FinalDimensions(img.dimensionsAsLongArray())
75+
76+
// convert the input parameters to meters (m)
77+
wavelength = wavelength.toFloat() * 1E-9
78+
lateral_res = lateral_res * 1E-6
79+
axial_res = axial_res * 1E-6
80+
pZ = pZ * 1E-6
81+
82+
// create the synthetic PSF
83+
psf = ops.op("create.kernelDiffraction").arity9().input(psf_size,
84+
numericalAperture,
85+
wavelength,
86+
riSample,
87+
riImmersion,
88+
lateral_res,
89+
axial_res,
90+
pZ,
91+
new FloatType()).apply()
92+
93+
// deconvolve image
94+
result = ops.op("deconvolve.richardsonLucyTV").arity8().input(img_float, psf, new FloatType(), new ComplexFloatType(), iterations, false, false, regularizationFactor).apply()
95+
96+
.. code-tab:: python
97+
98+
#@ OpEnvironment ops
99+
#@ ImgPlus img
100+
#@ Integer iterations(label="Iterations", value=30)
101+
#@ Float numericalAperture(label="Numerical Aperture", style="format:0.00", min=0.00, value=1.45)
102+
#@ Integer wavelength(label="Emission Wavelength (nm)", value=550)
103+
#@ Float riImmersion(label="Refractive Index (immersion)", style="format:0.00", min=0.00, value=1.5)
104+
#@ Float riSample(label="Refractive Index (sample)", style="format:0.00", min=0.00, value=1.4)
105+
#@ Float lateral_res(label="Lateral resolution (μm/pixel)", style="format:0.0000", min=0.0000, value=0.065)
106+
#@ Float axial_res(label="Axial resolution (μm/pixel)", style="format:0.0000", min=0.0000, value=0.1)
107+
#@ Float pZ(label="Particle/sample Position (μm)", style="format:0.0000", min=0.0000, value=0)
108+
#@ Float regularizationFactor(label="Regularization factor", style="format:0.00000", min=0.00000, value=0.002)
109+
#@output ImgPlus psf
110+
#@output ImgPlus result
111+
112+
from net.imglib2 import FinalDimensions
113+
from net.imglib2.type.numeric.real import FloatType
114+
from net.imglib2.type.numeric.complex import ComplexFloatType
115+
116+
# convert input image to float
117+
img_float = ops.op("create.img").arity2().input(img, FloatType()).apply()
118+
ops.op("convert.float32").arity1().input(img).output(img_float).compute()
119+
120+
# use image dimensions for PSF size
121+
psf_size = FinalDimensions(img.dimensionsAsLongArray())
122+
123+
# convert the input parameters to meters (m)
124+
wavelength = float(wavelength) * 1E-9
125+
lateral_res = lateral_res * 1E-6
126+
axial_res = axial_res * 1E-6
127+
pZ = pZ * 1E-6
128+
129+
# create the synthetic PSF
130+
psf = ops.op("create.kernelDiffraction").arity9().input(psf_size,
131+
numericalAperture,
132+
wavelength,
133+
riSample,
134+
riImmersion,
135+
lateral_res,
136+
axial_res,
137+
pZ,
138+
FloatType()).apply()
139+
140+
# deconvolve image
141+
result = ops.op("deconvolve.richardsonLucyTV").arity8().input(img_float, psf, FloatType(), ComplexFloatType(), iterations, False, False, regularizationFactor).apply()
142+
143+
| :sup:`1`: `Dey et. al, Micros Res Tech 2006`_
144+
| :sup:`2`: `Gibson & Lanni, JOSA 1992`_
145+
146+
.. _`Dey et. al, Micros Res Tech 2006`: https://pubmed.ncbi.nlm.nih.gov/16586486/
147+
.. _`Gibson & Lanni, JOSA 1992`: https://pubmed.ncbi.nlm.nih.gov/1738047/
148+
.. _`here`: https://media.imagej.net/sample_data/3d/hela_nucleus.tif
149+
.. _`script parameters`: https://imagej.net/scripting/parameters

‎docs/ops/doc/index.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@ The combination of these libraries allows declarative image analysis workflows,
2828
:maxdepth: 2
2929
:caption: ⚙️ Examples
3030

31+
examples/deconvolution
3132
examples/example_gaussian_subtraction
3233
examples/opencv_denoise
3334

‎scijava-ops-image/src/main/java/org/scijava/ops/image/deconvolve/AccelerationState.java‎

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
/*
2+
* #%L
3+
* ImageJ2 software for multidimensional image processing and analysis.
4+
* %%
5+
* Copyright (C) 2014 - 2023 ImageJ2 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 org.scijava.ops.image.deconvolve;
31+
32+
import net.imglib2.RandomAccessibleInterval;
33+
import net.imglib2.img.Img;
34+
import net.imglib2.type.numeric.RealType;
35+
36+
/**
37+
* Helper class for the Vector Accelerator Op. The Acceleration State stores the
38+
* Vector Accelerator's data while iterative deconvolution is applied.
39+
*
40+
* @author Curtis Rueden
41+
* @author Edward Evans
42+
* @param <T>
43+
*/
44+
45+
public class AccelerationState<T extends RealType<T>> {
46+
47+
private final RandomAccessibleInterval<T> ykIterated;
48+
private Img<T> xkm1Previous = null;
49+
private Img<T> ykPrediction = null;
50+
private Img<T> hkVector = null;
51+
private Img<T> gk;
52+
private Img<T> gkm1;
53+
private double accelerationFactor = 0.0f;
54+
55+
public AccelerationState(RandomAccessibleInterval<T> ykIterated) {
56+
this.ykIterated = ykIterated;
57+
}
58+
59+
public RandomAccessibleInterval<T> ykIterated() {
60+
return ykIterated;
61+
}
62+
63+
public Img<T> xkm1Previous() {
64+
return xkm1Previous;
65+
}
66+
67+
public void xkm1Previous(Img<T> xkm1Previous) {
68+
this.xkm1Previous = xkm1Previous;
69+
}
70+
71+
public Img<T> ykPrediction() {
72+
return ykPrediction;
73+
}
74+
75+
public void ykPrediction(Img<T> ykPrediction) {
76+
this.ykPrediction = ykPrediction;
77+
}
78+
79+
public Img<T> hkVector() {
80+
return hkVector;
81+
}
82+
83+
public void hkVector(Img<T> hkVector) {
84+
this.hkVector = hkVector;
85+
}
86+
87+
public Img<T> gk() {
88+
return gk;
89+
}
90+
91+
public void gk(Img<T> gk) {
92+
this.gk = gk;
93+
}
94+
95+
public Img<T> gkm1() {
96+
return gkm1;
97+
}
98+
99+
public void gkm1(Img<T> gkm1) {
100+
this.gkm1 = gkm1;
101+
}
102+
103+
public double accelerationFactor() {
104+
return accelerationFactor;
105+
}
106+
107+
public void accelerationFactor(double accelerationFactor) {
108+
this.accelerationFactor = accelerationFactor;
109+
}
110+
}

‎scijava-ops-image/src/main/java/org/scijava/ops/image/deconvolve/PadAndRichardsonLucy.java‎

Lines changed: 27 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@
3737
import net.imglib2.FinalDimensions;
3838
import net.imglib2.RandomAccessibleInterval;
3939
import net.imglib2.outofbounds.OutOfBoundsConstantValueFactory;
40+
import net.imglib2.outofbounds.OutOfBoundsMirrorFactory;
4041
import net.imglib2.outofbounds.OutOfBoundsFactory;
4142
import net.imglib2.type.NativeType;
4243
import net.imglib2.type.numeric.ComplexType;
@@ -63,7 +64,7 @@
6364
*/
6465
public class PadAndRichardsonLucy<I extends RealType<I> & NativeType<I>, O extends RealType<O> & NativeType<O>, K extends RealType<K> & NativeType<K>, C extends ComplexType<C> & NativeType<C>>
6566
implements
66-
Functions.Arity10<RandomAccessibleInterval<I>, RandomAccessibleInterval<K>, O, C, Integer, Boolean, Boolean, long[], OutOfBoundsFactory<I, RandomAccessibleInterval<I>>, OutOfBoundsFactory<K, RandomAccessibleInterval<K>>, RandomAccessibleInterval<O>>
67+
Functions.Arity9<RandomAccessibleInterval<I>, RandomAccessibleInterval<K>, O, C, Integer, Boolean, Boolean, long[], OutOfBoundsFactory<I, RandomAccessibleInterval<I>>, RandomAccessibleInterval<O>>
6768
{
6869

6970
private Computers.Arity1<RandomAccessibleInterval<O>, RandomAccessibleInterval<O>> computeEstimateOp =
@@ -78,9 +79,6 @@ public class PadAndRichardsonLucy<I extends RealType<I> & NativeType<I>, O exten
7879
@OpDependency(name = "math.multiply")
7980
private Computers.Arity2<RandomAccessibleInterval<O>, RandomAccessibleInterval<O>, RandomAccessibleInterval<O>> multiplyOp;
8081

81-
@OpDependency(name = "deconvolve.accelerate")
82-
private Inplaces.Arity1<RandomAccessibleInterval<O>> accelerator;
83-
8482
// TODO: can this go in AbstractFFTFilterF?
8583
@OpDependency(name = "create.img")
8684
private BiFunction<Dimensions, O, RandomAccessibleInterval<O>> outputCreator;
@@ -97,7 +95,7 @@ public class PadAndRichardsonLucy<I extends RealType<I> & NativeType<I>, O exten
9795
@OpDependency(name = "deconvolve.richardsonLucy")
9896
private Computers.Arity12<RandomAccessibleInterval<I>, RandomAccessibleInterval<K>, //
9997
RandomAccessibleInterval<C>, RandomAccessibleInterval<C>, Boolean, //
100-
Boolean, C, Integer, Inplaces.Arity1<RandomAccessibleInterval<O>>, //
98+
Boolean, C, Integer, Boolean, //
10199
Computers.Arity1<RandomAccessibleInterval<O>, RandomAccessibleInterval<O>>, //
102100
List<Inplaces.Arity1<RandomAccessibleInterval<O>>>, //
103101
RandomAccessibleInterval<O>, RandomAccessibleInterval<O>> richardsonLucyOp;
@@ -135,17 +133,17 @@ public class PadAndRichardsonLucy<I extends RealType<I> & NativeType<I>, O exten
135133

136134
return (input, kernel, out) -> {
137135
richardsonLucyOp.compute(input, kernel, fftImg, fftKernel, true, true,
138-
complexType, maxIterations, accelerate ? accelerator : null,
139-
computeEstimateOp, list, firstGuess.apply(raiExtendedInput, Util
140-
.getTypeFromInterval(out), out), out);
136+
complexType, maxIterations, accelerate, computeEstimateOp, list,
137+
firstGuess.apply(raiExtendedInput, Util.getTypeFromInterval(out),
138+
out), out);
141139
};
142140
}
143141

144142
// return a richardson lucy computer
145143
return (input, kernel, out) -> {
146144
richardsonLucyOp.compute(input, kernel, fftImg, fftKernel, true, true,
147-
complexType, maxIterations, accelerate ? accelerator : null,
148-
computeEstimateOp, null, null, out);
145+
complexType, maxIterations, accelerate, computeEstimateOp, null, null,
146+
out);
149147
};
150148
}
151149

@@ -207,30 +205,40 @@ private void computeFilter(final RandomAccessibleInterval<I> input,
207205
* @param accelerate indicates whether or not to use acceleration
208206
* @param borderSize
209207
* @param obfInput
210-
* @param obfKernel
211208
* @return the output
212209
*/
213210
@Override
214211
public RandomAccessibleInterval<O> apply(RandomAccessibleInterval<I> input,
215212
RandomAccessibleInterval<K> kernel, O outType, C complexType,
216213
Integer maxIterations, @Nullable Boolean nonCirculant,
217214
@Nullable Boolean accelerate, @Nullable long[] borderSize,
218-
@Nullable OutOfBoundsFactory<I, RandomAccessibleInterval<I>> obfInput,
219-
@Nullable OutOfBoundsFactory<K, RandomAccessibleInterval<K>> obfKernel)
215+
@Nullable OutOfBoundsFactory<I, RandomAccessibleInterval<I>> obfInput)
220216
{
221-
if (obfInput == null) obfInput = new OutOfBoundsConstantValueFactory<>(Util
222-
.getTypeFromInterval(input).createVariable());
223-
217+
// default to circulant
224218
if (nonCirculant == null) {
225-
this.nonCirculant = false;
219+
nonCirculant = false;
220+
this.nonCirculant = nonCirculant;
226221
}
227222
else {
228223
this.nonCirculant = nonCirculant;
229224
}
230225

231-
if (accelerate == null) {
232-
accelerate = false;
226+
// out of bounds factory will be different depending on if circulant or
227+
// non-circulant is used
228+
if (obfInput == null) {
229+
if (nonCirculant) {
230+
obfInput = new OutOfBoundsConstantValueFactory<>(Util
231+
.getTypeFromInterval(input).createVariable());
232+
}
233+
else {
234+
obfInput = new OutOfBoundsMirrorFactory<>(
235+
OutOfBoundsMirrorFactory.Boundary.SINGLE);
236+
}
233237
}
238+
239+
// default to no acceleration
240+
if (accelerate == null) accelerate = false;
241+
234242
this.maxIterations = maxIterations;
235243

236244
RandomAccessibleInterval<O> output = outputCreator.apply(input, outType);

0 commit comments

Comments
 (0)