Skip to content

Commit b01bb1a

Browse files
elevanselevans
committed
Convert SACA ImageJ Op to SciJava Op
This commit converts the imagej op copied from the saca branch to a SciJava op. This Op works, however it is not multithreaded and it produces NaN results at the edge of the heatmap.
1 parent a1b67de commit b01bb1a

File tree

3 files changed

+197
-134
lines changed

3 files changed

+197
-134
lines changed

‎scijava-ops-image/src/main/java/org/scijava/ops/image/coloc/saca/AdaptiveSmoothedKendallTau.java‎

Lines changed: 95 additions & 82 deletions
Original file line numberDiff line numberDiff line change
@@ -1,114 +1,125 @@
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+
*/
129

2-
package net.imagej.ops.coloc.saca;
30+
package org.scijava.ops.image.coloc.saca;
331

32+
import java.lang.Math;
433
import java.util.ArrayList;
534
import java.util.List;
635
import java.util.Random;
7-
import java.util.function.Function;
836

9-
import net.imglib2.Cursor;
10-
import net.imglib2.Localizable;
1137
import net.imglib2.RandomAccess;
12-
import net.imglib2.RandomAccessibleInterval;
38+
import net.imglib2.Localizable;
1339
import net.imglib2.loops.LoopBuilder;
1440
import net.imglib2.type.numeric.RealType;
15-
import net.imglib2.type.numeric.real.DoubleType;
16-
import net.imglib2.util.Intervals;
41+
import net.imglib2.RandomAccessibleInterval;
42+
import net.imglib2.img.ImgFactory;
1743
import net.imglib2.util.Localizables;
1844
import net.imglib2.util.Util;
1945
import net.imglib2.view.IntervalView;
2046
import net.imglib2.view.Views;
2147
import net.imglib2.view.composite.CompositeIntervalView;
2248
import net.imglib2.view.composite.GenericComposite;
49+
import net.imglib2.type.numeric.real.DoubleType;
2350

2451
/**
25-
* Adapted from Shulei's original Java code for AdaptiveSmoothedKendallTau from
26-
* his RKColocal R package.
27-
* (https://github.com/lakerwsl/RKColocal/blob/master/RKColocal_0.0.1.0000.tar.gz)
52+
* Helper class for Spatially Adaptive Colocalization Analysis (SACA) op.
2853
*
2954
* @author Shulei Wang
30-
* @author Curtis Rueden
31-
* @author Ellen T Arena
55+
* @author Ellen TA Dobson
3256
*/
57+
3358
public final class AdaptiveSmoothedKendallTau {
3459

3560
private AdaptiveSmoothedKendallTau() {}
3661

37-
public static <I extends RealType<I>, O extends RealType<O>> void execute(
38-
final RandomAccessibleInterval<I> image1,
39-
final RandomAccessibleInterval<I> image2, final I thres1, final I thres2,
40-
final RandomAccessibleInterval<O> result, final long seed)
41-
{
42-
execute(image1, image2, thres1, thres2, new DoubleType(), result, seed);
43-
}
44-
45-
public static <I extends RealType<I>, T extends RealType<T>, O extends RealType<O>> void execute(
46-
final RandomAccessibleInterval<I> image1,
47-
final RandomAccessibleInterval<I> image2, final I thres1, final I thres2,
48-
final T intermediate, final RandomAccessibleInterval<O> result,
49-
final long seed)
50-
{
51-
final Function<RandomAccessibleInterval<I>, RandomAccessibleInterval<T>> factory =
52-
img -> Util.getSuitableImgFactory(img, intermediate).create(img);
53-
execute(image1, image2, thres1, thres2, factory, result, seed);
54-
}
55-
56-
public static <I extends RealType<I>, T extends RealType<T>, O extends RealType<O>> void execute(
57-
final RandomAccessibleInterval<I> image1,
58-
final RandomAccessibleInterval<I> image2, final I thres1, final I thres2,
59-
Function<RandomAccessibleInterval<I>, RandomAccessibleInterval<T>> factory,
60-
final RandomAccessibleInterval<O> result, final long seed)
62+
// TODO: check that output float type is actually what we want here.
63+
public static <I extends RealType<I>> RandomAccessibleInterval<DoubleType>
64+
execute(final RandomAccessibleInterval<I> image1,
65+
final RandomAccessibleInterval<I> image2, final I thres1, final I thres2,
66+
final long seed)
6167
{
6268
final long nr = image1.dimension(1);
6369
final long nc = image1.dimension(0);
64-
final RandomAccessibleInterval<T> oldtau = factory.apply(image1);
65-
final RandomAccessibleInterval<T> newtau = factory.apply(image1);
66-
final RandomAccessibleInterval<T> oldsqrtN = factory.apply(image1);
67-
final RandomAccessibleInterval<T> newsqrtN = factory.apply(image1);
68-
final List<RandomAccessibleInterval<T>> stop = new ArrayList<>();
69-
for (int s = 0; s < 3; s++)
70-
stop.add(factory.apply(image1));
70+
final ImgFactory<DoubleType> factory = Util.getSuitableImgFactory(image1,
71+
new DoubleType());
72+
final RandomAccessibleInterval<DoubleType> oldtau = factory.create(image1);
73+
final RandomAccessibleInterval<DoubleType> newtau = factory.create(image1);
74+
final RandomAccessibleInterval<DoubleType> oldsqrtN = factory.create(
75+
image1);
76+
final RandomAccessibleInterval<DoubleType> newsqrtN = factory.create(
77+
image1);
78+
final RandomAccessibleInterval<DoubleType> result = factory.create(image1);
79+
final List<RandomAccessibleInterval<DoubleType>> stop = new ArrayList<>();
7180
final double Dn = Math.sqrt(Math.log(nr * nc)) * 2;
7281
final int TU = 15;
7382
final int TL = 8;
7483
final double Lambda = Dn;
75-
76-
LoopBuilder.setImages(oldsqrtN).forEachPixel(t -> t.setOne());
77-
84+
final double stepsize = 1.15;
85+
final Random rng = new Random(seed);
86+
boolean isCheck = false;
7887
double size = 1;
79-
final double stepsize = 1.15; // empirically the best, but could have users
80-
// set
8188
int intSize;
82-
boolean IsCheck = false;
8389

84-
final Random rng = new Random(seed);
90+
for (int s = 0; s < 3; s++)
91+
stop.add(factory.create(image1));
92+
93+
LoopBuilder.setImages(oldsqrtN).forEachPixel(t -> t.setOne());
8594

8695
for (int s = 0; s < TU; s++) {
8796
intSize = (int) Math.floor(size);
8897
singleiteration(image1, image2, thres1, thres2, stop, oldtau, oldsqrtN,
89-
newtau, newsqrtN, result, Lambda, Dn, intSize, IsCheck, rng);
98+
newtau, newsqrtN, result, Lambda, Dn, intSize, isCheck, rng);
9099
size *= stepsize;
91100
if (s == TL) {
92-
IsCheck = true;
101+
isCheck = true;
93102
LoopBuilder.setImages(stop.get(1), stop.get(2), newtau, newsqrtN)
94103
.forEachPixel((ts1, ts2, tTau, tSqrtN) -> {
95104
ts1.set(tTau);
96105
ts2.set(tSqrtN);
97106
});
98107
}
99108
}
109+
// get factory and create imgs
110+
return result;
100111
}
101112

102-
private static <I extends RealType<I>, T extends RealType<T>, O extends RealType<O>>
103-
void singleiteration(final RandomAccessibleInterval<I> image1,
104-
final RandomAccessibleInterval<I> image2, final I thres1, final I thres2,
105-
final List<RandomAccessibleInterval<T>> stop,
106-
final RandomAccessibleInterval<T> oldtau,
107-
final RandomAccessibleInterval<T> oldsqrtN,
108-
final RandomAccessibleInterval<T> newtau,
109-
final RandomAccessibleInterval<T> newsqrtN,
110-
final RandomAccessibleInterval<O> result, final double Lambda,
111-
final double Dn, final int Bsize, final boolean isCheck, final Random rng)
113+
private static <I extends RealType<I>> void singleiteration(
114+
final RandomAccessibleInterval<I> image1,
115+
final RandomAccessibleInterval<I> image2, final I thres1, final I thres2,
116+
final List<RandomAccessibleInterval<DoubleType>> stop,
117+
final RandomAccessibleInterval<DoubleType> oldtau,
118+
final RandomAccessibleInterval<DoubleType> oldsqrtN,
119+
final RandomAccessibleInterval<DoubleType> newtau,
120+
final RandomAccessibleInterval<DoubleType> newsqrtN,
121+
final RandomAccessibleInterval<DoubleType> result, final double Lambda,
122+
final double Dn, final int Bsize, final boolean isCheck, final Random rng)
112123
{
113124
final double[][] kernel = kernelGenerate(Bsize);
114125

@@ -127,25 +138,28 @@ void singleiteration(final RandomAccessibleInterval<I> image1,
127138
final double[] w2 = new double[totnum];
128139
final double[] cumw = new double[totnum];
129140

130-
RandomAccessibleInterval<T> workingImageStack = Views.stack(oldtau, newtau, oldsqrtN, newsqrtN, stop.get(0), stop.get(1), stop.get(2));
131-
CompositeIntervalView<T, ? extends GenericComposite<T>> workingImage =
141+
RandomAccessibleInterval<DoubleType> workingImageStack = Views.stack(oldtau,
142+
newtau, oldsqrtN, newsqrtN, stop.get(0), stop.get(1), stop.get(2));
143+
CompositeIntervalView<DoubleType, ? extends GenericComposite<DoubleType>> workingImage =
132144
Views.collapse(workingImageStack);
133145

134-
IntervalView<Localizable> positions = Views.interval( Localizables.randomAccessible(result.numDimensions() ), result );
146+
IntervalView<Localizable> positions = Views.interval(Localizables
147+
.randomAccessible(result.numDimensions()), result);
135148
final long nr = result.dimension(1);
136149
final long nc = result.dimension(0);
137150
final RandomAccess<I> gdImage1 = image1.randomAccess();
138151
final RandomAccess<I> gdImage2 = image2.randomAccess();
139-
final RandomAccess<T> gdTau = oldtau.randomAccess();
140-
final RandomAccess<T> gdSqrtN = oldsqrtN.randomAccess();
141-
LoopBuilder.setImages(positions, result, workingImage).forEachPixel((pos, resPixel, workingPixel) -> {
142-
T oldtauPix = workingPixel.get(0);
143-
T newtauPix = workingPixel.get(1);
144-
T oldsqrtNPix = workingPixel.get(2);
145-
T newsqrtNPix = workingPixel.get(3);
146-
T stop0Pix = workingPixel.get(4);
147-
T stop1Pix = workingPixel.get(5);
148-
T stop2Pix = workingPixel.get(6);
152+
final RandomAccess<DoubleType> gdTau = oldtau.randomAccess();
153+
final RandomAccess<DoubleType> gdSqrtN = oldsqrtN.randomAccess();
154+
LoopBuilder.setImages(positions, result, workingImage).forEachPixel((pos,
155+
resPixel, workingPixel) -> {
156+
DoubleType oldtauPix = workingPixel.get(0);
157+
DoubleType newtauPix = workingPixel.get(1);
158+
DoubleType oldsqrtNPix = workingPixel.get(2);
159+
DoubleType newsqrtNPix = workingPixel.get(3);
160+
DoubleType stop0Pix = workingPixel.get(4);
161+
DoubleType stop1Pix = workingPixel.get(5);
162+
DoubleType stop2Pix = workingPixel.get(6);
149163
final long row = pos.getLongPosition(1);
150164
updateRange(row, Bsize, nr, rowrange);
151165
if (isCheck) {
@@ -157,15 +171,15 @@ void singleiteration(final RandomAccessibleInterval<I> image1,
157171
updateRange(col, Bsize, nc, colrange);
158172
getData(Dn, kernel, gdImage1, gdImage2, gdTau, gdSqrtN, LocX, LocY, LocW,
159173
rowrange, colrange, totnum);
160-
newsqrtNPix.setReal(Math.sqrt(NTau(thres1, thres2, LocW, LocX,
161-
LocY)));
174+
newsqrtNPix.setReal(Math.sqrt(NTau(thres1, thres2, LocW, LocX, LocY)));
162175
if (newsqrtNPix.getRealDouble() <= 0) {
163176
newtauPix.setZero();
164177
resPixel.setZero();
165178
}
166179
else {
167-
final double tau = WtKendallTau.calculate(LocX, LocY, LocW, combinedData,
168-
rankedindex, rankedw, index1, index2, w1, w2, cumw, rng);
180+
final double tau = WtKendallTau.calculate(LocX, LocY, LocW,
181+
combinedData, rankedindex, rankedw, index1, index2, w1, w2, cumw,
182+
rng);
169183
newtauPix.setReal(tau);
170184
resPixel.setReal(tau * newsqrtNPix.getRealDouble() * 1.5);
171185
}
@@ -179,16 +193,15 @@ void singleiteration(final RandomAccessibleInterval<I> image1,
179193
newsqrtNPix.set(oldsqrtNPix);
180194
}
181195
}
182-
});
183-
196+
});
197+
184198
// TODO: instead of copying pixels here, swap oldTau and newTau every time.
185199
// :-)
186200
LoopBuilder.setImages(oldtau, newtau, oldsqrtN, newsqrtN).forEachPixel((
187201
tOldTau, tNewTau, tOldSqrtN, tNewSqrtN) -> {
188202
tOldTau.set(tNewTau);
189203
tOldSqrtN.set(tNewSqrtN);
190204
});
191-
192205
}
193206

194207
private static <I extends RealType<I>, T extends RealType<T>> void getData(

‎scijava-ops-image/src/main/java/org/scijava/ops/image/coloc/saca/SACA.java‎

Lines changed: 54 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
* #%L
33
* ImageJ software for multidimensional image processing and analysis.
44
* %%
5-
* Copyright (C) 2014 - 2018 ImageJ developers.
5+
* Copyright (C) 2014 - 2024 ImageJ developers.
66
* %%
77
* Redistribution and use in source and binary forms, with or without
88
* modification, are permitted provided that the following conditions are met:
@@ -27,62 +27,79 @@
2727
* #L%
2828
*/
2929

30-
package net.imagej.ops.coloc.saca;
30+
package org.scijava.ops.image.coloc.saca;
3131

32-
import net.imagej.ops.Ops;
33-
import net.imagej.ops.special.computer.AbstractBinaryComputerOp;
32+
import net.imglib2.FinalDimensions;
3433
import net.imglib2.RandomAccessibleInterval;
3534
import net.imglib2.histogram.Histogram1d;
3635
import net.imglib2.type.numeric.RealType;
36+
import net.imglib2.type.numeric.real.DoubleType;
3737
import net.imglib2.util.Intervals;
3838
import net.imglib2.view.Views;
3939

40-
import org.scijava.plugin.Parameter;
41-
import org.scijava.plugin.Plugin;
40+
import java.util.function.Function;
41+
42+
import org.scijava.function.Functions;
43+
import org.scijava.ops.spi.Nullable;
44+
import org.scijava.ops.spi.OpDependency;
4245

4346
/**
44-
* This algorithm is adapted from Spatially Adaptive Colocalization Analysis
45-
* (SACA) by Wang et al (2019); computes thresholds using Otsu method.
47+
* Spatially Adaptive Colocalization Analysis (SACA) Adapted from Shulei's
48+
* original Java code for AdaptiveSmoothedKendallTau from his RKColocal R
49+
* package.
50+
* (https://github.com/lakerwsl/RKColocal/blob/master/RKColocal_0.0.1.0000.tar.gz)
4651
*
47-
* @param <I> Type of the input images
48-
* @param <O> Type of the output image
52+
* @author Shulei Wang
53+
* @author Curtis Rueden
54+
* @author Ellen TA Dobson
55+
* @author Edward Evans
56+
* @param <I> input type
57+
* @implNote op names='coloc.saca', priority='100.'
4958
*/
50-
@Plugin(type = Ops.Coloc.SACA.class)
51-
public class SACA<I extends RealType<I>, O extends RealType<O>> extends
52-
AbstractBinaryComputerOp<RandomAccessibleInterval<I>, RandomAccessibleInterval<I>, RandomAccessibleInterval<O>>
53-
implements Ops.Coloc.SACA
59+
60+
public class SACA<I extends RealType<I>> implements
61+
Functions.Arity5<RandomAccessibleInterval<I>, RandomAccessibleInterval<I>, I, I, Long, RandomAccessibleInterval<DoubleType>>
5462
{
5563

56-
@Parameter(required = false)
57-
private I thres1;
64+
@OpDependency(name = "image.histogram")
65+
private Function<Iterable<I>, Histogram1d<I>> histOp;
5866

59-
@Parameter(required = false)
60-
private I thres2;
67+
@OpDependency(name = "threshold.otsu")
68+
private Function<Histogram1d<I>, I> otsuOp;
6169

62-
@Parameter(required = false)
63-
private long seed = 0xdeadbeef;
70+
/**
71+
* Spatially Adaptive Colocalization Analysis (SACA)
72+
*
73+
* @param image1 input image 1
74+
* @param image2 input image 2
75+
* @param thres1 threshold 1 value; otsu threshold applied if null
76+
* @param thres2 threshold 2 value; otsu threshold applied if null
77+
* @param seed seed to use; default 0xdeadbeefL
78+
* @return the output
79+
*/
6480

6581
@Override
66-
public void compute(final RandomAccessibleInterval<I> image1,
67-
final RandomAccessibleInterval<I> image2,
68-
final RandomAccessibleInterval<O> result)
82+
public RandomAccessibleInterval<DoubleType> apply(
83+
final RandomAccessibleInterval<I> image1,
84+
final RandomAccessibleInterval<I> image2, @Nullable I thres1,
85+
@Nullable I thres2, @Nullable Long seed)
6986
{
70-
71-
// check image sizes
72-
if (!(Intervals.equalDimensions(image1, image2))) {
73-
throw new IllegalArgumentException("Image dimensions do not match");
87+
// ensure images have the same dimensions
88+
FinalDimensions dims1 = new FinalDimensions(image1.dimensionsAsLongArray());
89+
FinalDimensions dims2 = new FinalDimensions(image2.dimensionsAsLongArray());
90+
if (!(Intervals.equalDimensions(dims1, dims2))) {
91+
throw new IllegalArgumentException(
92+
"Input image dimensions do not match.");
7493
}
7594

76-
// compute thresholds if necessary
77-
if (thres1 == null) thres1 = threshold(image1);
78-
if (thres2 == null) thres2 = threshold(image2);
79-
80-
AdaptiveSmoothedKendallTau.execute(image1, image2, thres1, thres2, result, seed);
81-
}
95+
// get seed and thresholds if necessary
96+
if (seed == null) seed = 0xdeadbeefL;
97+
if (thres1 == null) thres1 = otsuOp.apply(histOp.apply(Views.iterable(
98+
image1)));
99+
if (thres2 == null) thres2 = otsuOp.apply(histOp.apply(Views.iterable(
100+
image2)));
82101

83-
<V extends RealType<V>> V threshold(final RandomAccessibleInterval<V> image) {
84-
final Histogram1d<V> histogram = ops().image().histogram(Views.iterable(
85-
image));
86-
return ops().threshold().otsu(histogram);
102+
return AdaptiveSmoothedKendallTau.execute(image1, image2, thres1, thres2,
103+
seed);
87104
}
88105
}

0 commit comments

Comments
 (0)