Skip to content

Commit 9fe89b8

Browse files
ctruedenctrueden
committed
Merge branch 'pvalue'
This branch adds a coloc.pValue op that computes the p-value of a colocalization (or other) measure that produces a Double value. It does this by shuffling the input image in blocks and repeatedly applying the measure, then comparing against the original.
2 parents 2d52043 + e7b3dd3 commit 9fe89b8

File tree

7 files changed

+615
-2
lines changed

7 files changed

+615
-2
lines changed

‎src/main/java/net/imagej/ops/OpEnvironment.java‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@
3737
import java.util.Map;
3838

3939
import net.imagej.ops.cached.CachedOpEnvironment;
40+
import net.imagej.ops.coloc.ColocNamespace;
4041
import net.imagej.ops.convert.ConvertNamespace;
4142
import net.imagej.ops.copy.CopyNamespace;
4243
import net.imagej.ops.create.CreateNamespace;
@@ -829,6 +830,11 @@ default <I, O> RandomAccessibleInterval<O> slice(
829830

830831
// -- Operation shortcuts - other namespaces --
831832

833+
/** Gateway into ops of the "coloc" namespace. */
834+
default ColocNamespace coloc() {
835+
return namespace(ColocNamespace.class);
836+
}
837+
832838
/** Gateway into ops of the "copy" namespace. */
833839
default CopyNamespace copy() {
834840
return namespace(CopyNamespace.class);

‎src/main/java/net/imagej/ops/coloc/ColocNamespace.java‎

Lines changed: 56 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -29,14 +29,17 @@
2929

3030
package net.imagej.ops.coloc;
3131

32-
import org.scijava.plugin.Plugin;
33-
3432
import net.imagej.ops.AbstractNamespace;
3533
import net.imagej.ops.Namespace;
3634
import net.imagej.ops.OpMethod;
35+
import net.imagej.ops.special.function.BinaryFunctionOp;
36+
import net.imglib2.Dimensions;
37+
import net.imglib2.RandomAccessibleInterval;
3738
import net.imglib2.type.numeric.RealType;
3839
import net.imglib2.type.numeric.real.DoubleType;
3940

41+
import org.scijava.plugin.Plugin;
42+
4043
/**
4144
* The coloc namespace contains ops that facilitate colocalization analysis. b
4245
*
@@ -76,6 +79,57 @@ public <T extends RealType<T>, U extends RealType<U>> Double kendallTau(final It
7679
return result;
7780
}
7881

82+
// -- pValue --
83+
84+
@OpMethod(op = net.imagej.ops.coloc.pValue.PValue.class)
85+
public <T extends RealType<T>, U extends RealType<U>> Double pValue(
86+
final RandomAccessibleInterval<T> in1,
87+
final RandomAccessibleInterval<U> in2,
88+
final BinaryFunctionOp<Iterable<T>, Iterable<U>, Double> op)
89+
{
90+
final Double result = (Double) ops().run(
91+
net.imagej.ops.coloc.pValue.PValue.class, in1, in2, op);
92+
return result;
93+
}
94+
95+
@OpMethod(op = net.imagej.ops.coloc.pValue.PValue.class)
96+
public <T extends RealType<T>, U extends RealType<U>> Double pValue(
97+
final RandomAccessibleInterval<T> in1,
98+
final RandomAccessibleInterval<U> in2,
99+
final BinaryFunctionOp<Iterable<T>, Iterable<U>, Double> op,
100+
final int nrRandomizations)
101+
{
102+
final Double result = (Double) ops().run(
103+
net.imagej.ops.coloc.pValue.PValue.class, in1, in2, op, nrRandomizations);
104+
return result;
105+
}
106+
107+
@OpMethod(op = net.imagej.ops.coloc.pValue.PValue.class)
108+
public <T extends RealType<T>, U extends RealType<U>> Double pValue(
109+
final RandomAccessibleInterval<T> in1,
110+
final RandomAccessibleInterval<U> in2,
111+
final BinaryFunctionOp<Iterable<T>, Iterable<U>, Double> op,
112+
final int nrRandomizations, final Dimensions psfSize)
113+
{
114+
final Double result = (Double) ops().run(
115+
net.imagej.ops.coloc.pValue.PValue.class, in1, in2, op, nrRandomizations,
116+
psfSize);
117+
return result;
118+
}
119+
120+
@OpMethod(op = net.imagej.ops.coloc.pValue.PValue.class)
121+
public <T extends RealType<T>, U extends RealType<U>> Double pValue(
122+
final RandomAccessibleInterval<T> in1,
123+
final RandomAccessibleInterval<U> in2,
124+
final BinaryFunctionOp<Iterable<T>, Iterable<U>, Double> op,
125+
final int nrRandomizations, final Dimensions psfSize, final long seed)
126+
{
127+
final Double result = (Double) ops().run(
128+
net.imagej.ops.coloc.pValue.PValue.class, in1, in2, op, nrRandomizations,
129+
psfSize, seed);
130+
return result;
131+
}
132+
79133
// -- Namespace methods --
80134

81135
@Override

‎src/main/java/net/imagej/ops/coloc/ShuffledView.java‎

Lines changed: 162 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,162 @@
1+
2+
package net.imagej.ops.coloc;
3+
4+
import java.util.Collections;
5+
import java.util.List;
6+
import java.util.Random;
7+
8+
import net.imglib2.AbstractInterval;
9+
import net.imglib2.Interval;
10+
import net.imglib2.Point;
11+
import net.imglib2.RandomAccess;
12+
import net.imglib2.RandomAccessibleInterval;
13+
import net.imglib2.Sampler;
14+
import net.imglib2.View;
15+
import net.imglib2.util.IntervalIndexer;
16+
17+
import org.scijava.util.IntArray;
18+
19+
/**
20+
* Randomly shuffles an image blockwise.
21+
*
22+
* @author Curtis Rueden
23+
* @author Ellen T Arena
24+
*
25+
* @param <T> Type of image to be shuffled.
26+
*/
27+
public class ShuffledView<T> extends AbstractInterval implements
28+
RandomAccessibleInterval<T>, View
29+
{
30+
31+
private static Random rng;
32+
private final RandomAccessibleInterval<T> image;
33+
private List<Integer> blockIndices;
34+
private int[] blockSize;
35+
private int[] blockDims;
36+
37+
public ShuffledView(final RandomAccessibleInterval<T> image,
38+
final int[] blockSize, final long seed)
39+
{
40+
this(image, blockSize, null, seed);
41+
}
42+
43+
public ShuffledView(final RandomAccessibleInterval<T> image,
44+
final int[] blockSize, final List<Integer> blockIndices)
45+
{
46+
this(image, blockSize, blockIndices, 0);
47+
}
48+
49+
private ShuffledView(final RandomAccessibleInterval<T> image,
50+
final int[] blockSize, final List<Integer> blockIndices, final long seed)
51+
{
52+
super(image); // uses same bounds as the input image
53+
this.image = image;
54+
this.blockSize = blockSize;
55+
56+
// compute some info about our block sizes
57+
final int numDims = image.numDimensions();
58+
blockDims = new int[numDims];
59+
long totalBlocks = 1;
60+
for (int d = 0; d < numDims; d++) {
61+
final long blockDim = image.dimension(d) / blockSize[d];
62+
if (blockDim * blockSize[d] != image.dimension(d)) {
63+
throw new IllegalArgumentException("Image dimension #" + d +
64+
" is not evenly divisible by block size:" + blockSize[d]);
65+
}
66+
if (blockDim > Integer.MAX_VALUE) {
67+
throw new UnsupportedOperationException("Block dimension #" + d +
68+
" is too large: " + blockDim);
69+
}
70+
blockDims[d] = (int) blockDim;
71+
totalBlocks *= blockDims[d];
72+
}
73+
if (totalBlocks > Integer.MAX_VALUE) {
74+
throw new UnsupportedOperationException("Too many blocks: " +
75+
totalBlocks);
76+
}
77+
if (blockIndices == null) {
78+
this.blockIndices = createBlocks((int) totalBlocks);
79+
rng = new Random(seed);
80+
shuffleBlocks();
81+
} else {
82+
this.blockIndices = blockIndices;
83+
}
84+
}
85+
86+
private static List<Integer> createBlocks(final int blockCount)
87+
{
88+
// generate the identity mapping of indices
89+
final IntArray blocks = new IntArray();
90+
blocks.ensureCapacity(blockCount);
91+
for (int b = 0; b < blockCount; b++)
92+
blocks.addValue(b);
93+
return blocks;
94+
}
95+
96+
public void shuffleBlocks() {
97+
if (rng == null) {
98+
throw new IllegalStateException("No seed provided. Cannot shuffle.");
99+
}
100+
Collections.shuffle(blockIndices, rng);
101+
}
102+
103+
@Override
104+
public RandomAccess<T> randomAccess() {
105+
return new ShuffledRandomAccess();
106+
}
107+
108+
@Override
109+
public RandomAccess<T> randomAccess(final Interval interval) {
110+
return randomAccess(); // FIXME
111+
}
112+
113+
private class ShuffledRandomAccess extends Point implements RandomAccess<T> {
114+
115+
public ShuffledRandomAccess() {
116+
super(image.numDimensions());
117+
}
118+
119+
@Override
120+
public T get() {
121+
// Convert from image coordinates to block coordinates.
122+
final long[] blockPos = new long[position.length];
123+
final long[] blockOffset = new long[position.length];
124+
for (int d = 0; d < position.length; d++) {
125+
blockPos[d] = position[d] / blockSize[d];
126+
blockOffset[d] = position[d] % blockSize[d];
127+
}
128+
129+
// Convert N-D block coordinates to 1D block index.
130+
final int blockIndex = IntervalIndexer.positionToIndex(blockPos,
131+
blockDims);
132+
133+
// Map block index to shuffled block index.
134+
final int shuffledBlockIndex = blockIndices.get(blockIndex);
135+
136+
// Now convert our 1D shuffled block index back to N-D block
137+
// coordinates.
138+
final long[] shuffledBlockPos = new long[position.length];
139+
IntervalIndexer.indexToPosition(shuffledBlockIndex, blockDims,
140+
shuffledBlockPos);
141+
142+
// Finally, position the original image according to our shuffled
143+
// position.
144+
final RandomAccess<T> imageRA = image.randomAccess();
145+
for (int d = 0; d < position.length; d++) {
146+
final long pd = shuffledBlockPos[d] * blockSize[d] + blockOffset[d];
147+
imageRA.setPosition(pd, d);
148+
}
149+
return imageRA.get();
150+
}
151+
152+
@Override
153+
public Sampler<T> copy() {
154+
throw new UnsupportedOperationException();
155+
}
156+
157+
@Override
158+
public RandomAccess<T> copyRandomAccess() {
159+
throw new UnsupportedOperationException();
160+
}
161+
}
162+
}

‎src/main/java/net/imagej/ops/coloc/pValue/PValue.java‎

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
2+
package net.imagej.ops.coloc.pValue;
3+
4+
import net.imagej.ops.Ops;
5+
import net.imagej.ops.coloc.ShuffledView;
6+
import net.imagej.ops.special.function.AbstractBinaryFunctionOp;
7+
import net.imagej.ops.special.function.BinaryFunctionOp;
8+
import net.imglib2.Dimensions;
9+
import net.imglib2.IterableInterval;
10+
import net.imglib2.RandomAccessibleInterval;
11+
import net.imglib2.type.numeric.RealType;
12+
import net.imglib2.util.Intervals;
13+
import net.imglib2.view.Views;
14+
15+
import org.scijava.plugin.Parameter;
16+
import org.scijava.plugin.Plugin;
17+
18+
/**
19+
* This algorithm repeatedly executes a colocalization algorithm, computing a
20+
* p-value. It is based on a new statistical framework published by Wang et al
21+
* (2017) IEEE Signal Processing "Automated and Robust Quantification of
22+
* Colocalization in Dual-Color Fluorescence Microscopy: A Nonparametric
23+
* Statistical Approach".
24+
*/
25+
@Plugin(type = Ops.Coloc.PValue.class)
26+
public class PValue<T extends RealType<T>, U extends RealType<U>> extends
27+
AbstractBinaryFunctionOp<RandomAccessibleInterval<T>, RandomAccessibleInterval<U>, Double>
28+
implements Ops.Coloc.PValue
29+
{
30+
31+
@Parameter
32+
private BinaryFunctionOp<Iterable<T>, Iterable<U>, Double> op;
33+
34+
@Parameter(required = false)
35+
private int nrRandomizations = 1000;
36+
37+
@Parameter(required = false)
38+
private Dimensions psfSize;
39+
40+
@Parameter(required = false)
41+
private long seed = 0x27372034;
42+
43+
@Override
44+
public Double calculate(final RandomAccessibleInterval<T> image1,
45+
final RandomAccessibleInterval<U> image2)
46+
{
47+
final int[] blockSize = blockSize(image1, psfSize);
48+
final RandomAccessibleInterval<T> trimmedImage1 = trim(image1, blockSize);
49+
final RandomAccessibleInterval<U> trimmedImage2 = trim(image2, blockSize);
50+
51+
final ShuffledView<T> shuffled = new ShuffledView<>(image1, blockSize,
52+
seed);
53+
final IterableInterval<T> shuffledIterable = Views.iterable(shuffled);
54+
final double[] sampleDistribution = new double[nrRandomizations];
55+
56+
final IterableInterval<T> iterableImage1 = Views.iterable(trimmedImage1);
57+
final IterableInterval<U> iterableImage2 = Views.iterable(trimmedImage2);
58+
final double value = op.calculate(iterableImage1, iterableImage2);
59+
60+
for (int i = 0; i < nrRandomizations; i++) {
61+
shuffled.shuffleBlocks();
62+
sampleDistribution[i] = op.calculate(shuffledIterable, iterableImage2);
63+
}
64+
return calculatePvalue(value, sampleDistribution);
65+
}
66+
67+
private double calculatePvalue(final double input,
68+
final double[] distribution)
69+
{
70+
double count = 0;
71+
for (int i = 0; i < distribution.length; i++) {
72+
if (distribution[i] > input) {
73+
count++;
74+
}
75+
}
76+
final double pvalue = count / distribution.length;
77+
return pvalue;
78+
}
79+
80+
private static int[] blockSize(final Dimensions image,
81+
final Dimensions psfSize)
82+
{
83+
if (psfSize != null) return Intervals.dimensionsAsIntArray(psfSize);
84+
85+
final int[] blockSize = new int[image.numDimensions()];
86+
for (int d = 0; d < blockSize.length; d++) {
87+
final long size = (long) Math.floor(Math.sqrt(image.dimension(d)));
88+
if (size > Integer.MAX_VALUE) {
89+
throw new IllegalArgumentException("Image dimension #" + d +
90+
" is too large: " + image.dimension(d));
91+
}
92+
blockSize[d] = (int) size;
93+
}
94+
return blockSize;
95+
}
96+
97+
private static <V> RandomAccessibleInterval<V> trim(
98+
final RandomAccessibleInterval<V> image, final int[] blockSize)
99+
{
100+
final long[] min = Intervals.minAsLongArray(image);
101+
final long[] max = Intervals.maxAsLongArray(image);
102+
for (int d = 0; d < blockSize.length; d++) {
103+
final long trimSize = image.dimension(d) % blockSize[d];
104+
final long half = trimSize / 2;
105+
min[d] += half;
106+
max[d] -= trimSize - half;
107+
}
108+
return Views.interval(image, min, max);
109+
}
110+
}

‎src/main/templates/net/imagej/ops/Ops.list‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ namespaces = ```
2929
[name: "kendallTau", iface: "KendallTau"],
3030
[name: "manders", iface: "Manders"],
3131
[name: "pearsons", iface: "Pearsons"],
32+
[name: "pValue", iface: "PValue"],
3233
[name: "spearman", iface: "Spearman"],
3334
]],
3435
[name: "convert", iface: "Convert", ops: [

0 commit comments

Comments
 (0)