Page MenuHomec4science
Files F61098826

PearsonsCorrelation.java
No OneTemporary
Actions

File Metadata

Created
Sat, May 4, 13:00

PearsonsCorrelation.java
View Options

/*-
* #%L
* Fiji's plugin for colocalization analysis.
* %%
* Copyright (C) 2009 - 2017 Fiji developers.
* %%
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3 of the
* License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public
* License along with this program. If not, see
* <http://www.gnu.org/licenses/gpl-3.0.html>.
* #L%
*/
package sc.fiji.coloc.algorithms;
import net.imglib2.RandomAccessibleInterval;
import net.imglib2.TwinCursor;
import net.imglib2.algorithm.math.ImageStatistics;
import net.imglib2.type.logic.BitType;
import net.imglib2.type.numeric.RealType;
import net.imglib2.view.Views;
import sc.fiji.coloc.gadgets.DataContainer;
import sc.fiji.coloc.gadgets.MaskFactory;
import sc.fiji.coloc.gadgets.ThresholdMode;
import sc.fiji.coloc.results.ResultHandler;
/**
* A class that represents the mean calculation of the two source
* images in the data container.
*
* @param <T>
*/
public class PearsonsCorrelation<T extends RealType< T >> extends Algorithm<T> {
// Identifiers for choosing which implementation to use
public enum Implementation {Classic, Fast};
// The member variable to store the implementation of the Pearson's Coefficient calculation used.
Implementation theImplementation = Implementation.Fast;
// resulting Pearsons value without thresholds
double pearsonsCorrelationValue;
// resulting Pearsons value below threshold
double pearsonsCorrelationValueBelowThr;
// resulting Pearsons value above threshold
double pearsonsCorrelationValueAboveThr;
/**
* Creates a new Pearson's Correlation and allows us to define
* which implementation of the calculation to use.
* @param implementation The implementation of Pearson's Coefficient calculation to use.
*/
public PearsonsCorrelation(Implementation implementation) {
super("Pearson correlation");
this.theImplementation = implementation;
}
/**
* Creates a new Pearson's Correlation with default (fast) implementation parameter.
*/
public PearsonsCorrelation() {
this(Implementation.Fast);
}
@Override
public void execute(DataContainer<T> container) throws MissingPreconditionException {
// get the 2 images for the calculation of Pearson's
RandomAccessibleInterval<T> img1 = container.getSourceImage1();
RandomAccessibleInterval<T> img2 = container.getSourceImage2();
RandomAccessibleInterval<BitType> mask = container.getMask();
// get the thresholds of the images
AutoThresholdRegression<T> autoThreshold = container.getAutoThreshold();
if (autoThreshold == null ) {
throw new MissingPreconditionException("Pearsons calculation need AutoThresholdRegression to be run before it.");
}
T threshold1 = autoThreshold.getCh1MaxThreshold();
T threshold2 = autoThreshold.getCh2MaxThreshold();
if (threshold1 == null || threshold2 == null ) {
throw new MissingPreconditionException("Pearsons calculation needs valid (not null) thresholds.");
}
/* Create cursors to walk over the images. First go over the
* images without a mask. */
TwinCursor<T> cursor = new TwinCursor<T>(
img1.randomAccess(), img2.randomAccess(),
Views.iterable(mask).localizingCursor());
MissingPreconditionException error = null;
if (theImplementation == Implementation.Classic) {
// get the means from the DataContainer
double ch1Mean = container.getMeanCh1();
double ch2Mean = container.getMeanCh2();
try {
cursor.reset();
pearsonsCorrelationValue = classicPearsons(cursor,
ch1Mean, ch2Mean);
} catch (MissingPreconditionException e) {
// probably a numerical error occurred
pearsonsCorrelationValue = Double.NaN;
error = e;
}
try {
cursor.reset();
pearsonsCorrelationValueBelowThr = classicPearsons(cursor,
ch1Mean, ch2Mean, threshold1, threshold2, ThresholdMode.Below);
} catch (MissingPreconditionException e) {
// probably a numerical error occurred
pearsonsCorrelationValueBelowThr = Double.NaN;
error = e;
}
try {
cursor.reset();
pearsonsCorrelationValueAboveThr = classicPearsons(cursor,
ch1Mean, ch2Mean, threshold1, threshold2, ThresholdMode.Above);
} catch (MissingPreconditionException e) {
// probably a numerical error occurred
pearsonsCorrelationValueAboveThr = Double.NaN;
error = e;
}
}
else if (theImplementation == Implementation.Fast) {
try {
cursor.reset();
pearsonsCorrelationValue = fastPearsons(cursor);
} catch (MissingPreconditionException e) {
// probably a numerical error occurred
pearsonsCorrelationValue = Double.NaN;
error = e;
}
try {
cursor.reset();
pearsonsCorrelationValueBelowThr = fastPearsons(cursor,
threshold1, threshold2, ThresholdMode.Below);
} catch (MissingPreconditionException e) {
// probably a numerical error occurred
pearsonsCorrelationValueBelowThr = Double.NaN;
error = e;
}
try {
cursor.reset();
pearsonsCorrelationValueAboveThr = fastPearsons(cursor,
threshold1, threshold2, ThresholdMode.Above);
} catch (MissingPreconditionException e) {
// probably a numerical error occurred
pearsonsCorrelationValueAboveThr = Double.NaN;
error = e;
}
}
// if an error occurred, throw it one level up
if (error != null)
throw error;
}
/**
* Calculates Pearson's R value without any constraint in values, thus it uses no thresholds.
* If additional data like the images mean is needed, it is calculated.
*
* @param <S> The images base type.
* @param img1 The first image to walk over.
* @param img2 The second image to walk over.
* @return Pearson's R value.
* @throws MissingPreconditionException
*/
public <S extends RealType<S>> double calculatePearsons(
RandomAccessibleInterval<S> img1, RandomAccessibleInterval<S> img2)
throws MissingPreconditionException {
// create an "always true" mask to walk over the images
final long[] dims = new long[img1.numDimensions()];
img1.dimensions(dims);
RandomAccessibleInterval<BitType> alwaysTrueMask = MaskFactory.createMask(dims, true);
return calculatePearsons(img1, img2, alwaysTrueMask);
}
/**
* Calculates Pearson's R value without any constraint in values, thus it uses no
* thresholds. A mask is required to mark which data points should be visited. If
* additional data like the images mean is needed, it is calculated.
*
* @param <S> The images base type.
* @param img1 The first image to walk over.
* @param img2 The second image to walk over.
* @param mask A mask for the images.
* @return Pearson's R value.
* @throws MissingPreconditionException
*/
public <S extends RealType<S>> double calculatePearsons(
RandomAccessibleInterval<S> img1, RandomAccessibleInterval<S> img2,
RandomAccessibleInterval<BitType> mask) throws MissingPreconditionException {
TwinCursor<S> cursor = new TwinCursor<S>(
img1.randomAccess(), img2.randomAccess(),
Views.iterable(mask).localizingCursor());
double r;
if (theImplementation == Implementation.Classic) {
/* since we need the means and apparently don't have them,
* calculate them.
*/
double mean1 = ImageStatistics.getImageMean(img1);
double mean2 = ImageStatistics.getImageMean(img2);
// do the actual calculation
r = classicPearsons(cursor, mean1, mean2);
} else {
r = fastPearsons(cursor);
}
return r;
}
/**
* Calculates Pearson's R value with the possibility to constraint in values.
* This could be useful of one wants to apply thresholds. You need to provide
* the images means, albeit not used by all implementations.
*
* @param <S> The images base type.
* @param cursor The cursor to walk over both images.
* @return Pearson's R value.
* @throws MissingPreconditionException
*/
public <S extends RealType<S>> double calculatePearsons(TwinCursor<S> cursor,
double mean1, double mean2, S thresholdCh1, S thresholdCh2,
ThresholdMode tMode) throws MissingPreconditionException {
if (theImplementation == Implementation.Classic) {
// do the actual calculation
return classicPearsons(cursor, mean1, mean2,
thresholdCh1, thresholdCh2, tMode);
} else {
return fastPearsons(cursor, thresholdCh1,
thresholdCh2, tMode);
}
}
/**
* Calculates Person's R value by using a Classic implementation of the
* algorithm. This method allows the specification of a TwinValueRangeCursor.
* With such a cursor one for instance can combine different thresholding
* conditions for each channel. The cursor is not closed in here.
*
* @param <T> The image base type
* @param cursor The cursor that defines the walk over both images.
* @param meanCh1 Mean of channel 1.
* @param meanCh2 Mean of channel 2.
* @return Person's R value
*/
public static <T extends RealType<T>> double classicPearsons(TwinCursor<T> cursor,
double meanCh1, double meanCh2) throws MissingPreconditionException {
return classicPearsons(cursor, meanCh1, meanCh2, null, null, ThresholdMode.None);
}
public static <T extends RealType<T>> double classicPearsons(TwinCursor<T> cursor,
double meanCh1, double meanCh2, final T thresholdCh1, final T thresholdCh2,
ThresholdMode tMode) throws MissingPreconditionException {
// the actual accumulation of the image values is done in a separate object
Accumulator<T> acc;
if (tMode == ThresholdMode.None) {
acc = new Accumulator<T>(cursor, meanCh1, meanCh2) {
@Override
final public boolean accept(T type1, T type2) {
return true;
}
};
} else if (tMode == ThresholdMode.Below) {
acc = new Accumulator<T>(cursor, meanCh1, meanCh2) {
@Override
final public boolean accept(T type1, T type2) {
return type1.compareTo(thresholdCh1) < 0 ||
type2.compareTo(thresholdCh2) < 0;
}
};
} else if (tMode == ThresholdMode.Above) {
acc = new Accumulator<T>(cursor, meanCh1, meanCh2) {
@Override
final public boolean accept(T type1, T type2) {
return type1.compareTo(thresholdCh1) > 0 ||
type2.compareTo(thresholdCh2) > 0;
}
};
} else {
throw new UnsupportedOperationException();
}
double pearsonsR = acc.xy / Math.sqrt(acc.xx * acc.yy);
checkForSanity(pearsonsR, acc.count);
return pearsonsR;
}
/**
* Calculates Person's R value by using a fast implementation of the
* algorithm. This method allows the specification of a TwinValueRangeCursor.
* With such a cursor one for instance can combine different thresholding
* conditions for each channel. The cursor is not closed in here.
*
* @param <T> The image base type
* @param cursor The cursor that defines the walk over both images.
* @return Person's R value
*/
public static <T extends RealType<T>> double fastPearsons(TwinCursor<T> cursor)
throws MissingPreconditionException {
return fastPearsons(cursor, null, null, ThresholdMode.None);
}
public static <T extends RealType<T>> double fastPearsons(TwinCursor<T> cursor,
final T thresholdCh1, final T thresholdCh2, ThresholdMode tMode)
throws MissingPreconditionException {
// the actual accumulation of the image values is done in a separate object
Accumulator<T> acc;
if (tMode == ThresholdMode.None) {
acc = new Accumulator<T>(cursor) {
@Override
final public boolean accept(T type1, T type2) {
return true;
}
};
} else if (tMode == ThresholdMode.Below) {
acc = new Accumulator<T>(cursor) {
@Override
final public boolean accept(T type1, T type2) {
return type1.compareTo(thresholdCh1) < 0 ||
type2.compareTo(thresholdCh2) < 0;
}
};
} else if (tMode == ThresholdMode.Above) {
acc = new Accumulator<T>(cursor) {
@Override
final public boolean accept(T type1, T type2) {
return type1.compareTo(thresholdCh1) > 0 ||
type2.compareTo(thresholdCh2) > 0;
}
};
} else {
throw new UnsupportedOperationException();
}
// for faster computation, have the inverse of N available
double invCount = 1.0 / acc.count;
double pearsons1 = acc.xy - (acc.x * acc.y * invCount);
double pearsons2 = acc.xx - (acc.x * acc.x * invCount);
double pearsons3 = acc.yy - (acc.y * acc.y * invCount);
double pearsonsR = pearsons1 / (Math.sqrt(pearsons2 * pearsons3));
checkForSanity(pearsonsR, acc.count);
return pearsonsR;
}
/**
* Does a sanity check for calculated Pearsons values. Wrong
* values can happen for fast and classic implementation.
*
* @param val The value to check.
*/
private static void checkForSanity(double value, int iterations) throws MissingPreconditionException {
if ( Double.isNaN(value) || Double.isInfinite(value)) {
/* For the _fast_ implementation this could happen:
* Infinity could happen if only the numerator is 0, i.e.:
* sum1squared == sum1 * sum1 * invN
* and
* sum2squared == sum2 * sum2 * invN
* If the denominator is also zero, one will get NaN, i.e:
* sumProduct1_2 == sum1 * sum2 * invN
*
* For the classic implementation it could happen, too:
* Infinity happens if one channels sum of value-mean-differences
* is zero. If it is negative for one image you will get NaN.
* Additionally, if is zero for both channels at once you
* could get NaN. NaN
*/
throw new MissingPreconditionException("A numerical problem occured: the input data is unsuitable for this algorithm. Possibly too few pixels (in range were: " + iterations + ").");
}
}
@Override
public void processResults(ResultHandler<T> handler) {
super.processResults(handler);
handler.handleValue("Pearson's R value (no threshold)", pearsonsCorrelationValue, 2);
handler.handleValue("Pearson's R value (below threshold)", pearsonsCorrelationValueBelowThr, 2);
handler.handleValue("Pearson's R value (above threshold)", pearsonsCorrelationValueAboveThr, 2);
}
public double getPearsonsCorrelationValue() {
return pearsonsCorrelationValue;
}
public double getPearsonsCorrelationBelowThreshold() {
return pearsonsCorrelationValueBelowThr;
}
public double getPearsonsCorrelationAboveThreshold() {
return pearsonsCorrelationValueAboveThr;
}
}

Event Timeline

c4science · Help