Skip to content

Commit

Permalink
Merge pull request #5 from tomka/add-costes-auto-threshold-alternative
Browse files Browse the repository at this point in the history
Add and default to regression pattern used by Costes
  • Loading branch information
tomka committed Jul 29, 2015
2 parents 62fa69d + 646d38d commit 5665600
Show file tree
Hide file tree
Showing 14 changed files with 262 additions and 70 deletions.
27 changes: 24 additions & 3 deletions src/main/java/Coloc_2.java
@@ -1,5 +1,6 @@
import algorithms.Algorithm;
import algorithms.AutoThresholdRegression;
import algorithms.AutoThresholdRegression.Implementation;
import algorithms.CostesSignificanceTest;
import algorithms.Histogram2D;
import algorithms.InputCheck;
Expand Down Expand Up @@ -45,7 +46,6 @@
import net.imglib2.type.NativeType;
import net.imglib2.type.logic.BitType;
import net.imglib2.type.numeric.RealType;
import net.imglib2.util.Util;
import net.imglib2.view.Views;
import results.PDFWriter;
import results.ResultHandler;
Expand Down Expand Up @@ -119,11 +119,16 @@ protected enum RoiConfiguration {
// A list of all ROIs/masks found
protected ArrayList<MaskInfo> masks = new ArrayList<MaskInfo>();

// default indices of image, mask and ROI choices
// A list of auto threshold implementations
protected String[] regressions = new String[
AutoThresholdRegression.Implementation.values().length];

// default indices of image, mask, ROI and regression choices
protected static int index1 = 0;
protected static int index2 = 1;
protected static int indexMask = 0;
protected static int indexRoi = 0;
protected static int indexRegr = 0;

// the images to work on
protected Img<T> img1, img2;
Expand Down Expand Up @@ -203,6 +208,13 @@ public boolean showDialog() {
}
}

// find all available regression strategies
Implementation[] regressionImplementations =
AutoThresholdRegression.Implementation.values();
for (int i=0; i<regressionImplementations.length; ++i) {
regressions[i] = regressionImplementations[i].toString();
}

// set up the users preferences
displayImages = Prefs.get(PREF_KEY+"displayImages", false);
autoSavePdf = Prefs.get(PREF_KEY+"autoSavePdf", true);
Expand All @@ -217,6 +229,7 @@ public boolean showDialog() {
boolean useCostes = Prefs.get(PREF_KEY+"useCostes", true);
int psf = (int) Prefs.get(PREF_KEY+"psf", 3);
int nrCostesRandomisations = (int) Prefs.get(PREF_KEY+"nrCostesRandomisations", 10);
indexRegr = (int) Prefs.get(PREF_KEY+"regressionImplementation", 0);

/* make sure the default indices are no bigger
* than the amount of images we have
Expand All @@ -230,6 +243,9 @@ public boolean showDialog() {
gd.addChoice("ROI_or_mask", roisAndMasks, roisAndMasks[indexMask]);
//gd.addChoice("Use ROI", roiLabels, roiLabels[indexRoi]);

gd.addChoice("Threshold regression", regressions,
regressions[indexRegr]);

gd.addCheckbox("Show_\"Save_PDF\"_Dialog", autoSavePdf);
gd.addCheckbox("Display_Images_in_Result", displayImages);
gd.addCheckbox("Display_Shuffled_Images", displayShuffledCostes);
Expand Down Expand Up @@ -324,6 +340,9 @@ else if (indexMask == 3)
masks.add(new MaskInfo(null, null));
}

// get information about the mask/ROI to use
indexRegr = gd.getNextChoiceIndex();

// read out GUI data
autoSavePdf = gd.getNextBoolean();
displayImages = gd.getNextBoolean();
Expand All @@ -340,6 +359,7 @@ else if (indexMask == 3)
nrCostesRandomisations = (int) gd.getNextNumber();

// save user preferences
Prefs.set(PREF_KEY+"regressionImplementation", indexRegr);
Prefs.set(PREF_KEY+"autoSavePdf", autoSavePdf);
Prefs.set(PREF_KEY+"displayImages", displayImages);
Prefs.set(PREF_KEY+"displayShuffledCostes", displayShuffledCostes);
Expand Down Expand Up @@ -424,7 +444,8 @@ public void colocalise(Img<T> img1, Img<T> img2, BoundingBox roi,
userSelectedJobs.add( container.setInputCheck(
new InputCheck<T>()) );
userSelectedJobs.add( container.setAutoThreshold(
new AutoThresholdRegression<T>(pearsonsCorrelation)) );
new AutoThresholdRegression<T>(pearsonsCorrelation,
AutoThresholdRegression.Implementation.values()[indexRegr])));

// add user selected algorithms
addIfValid(pearsonsCorrelation, userSelectedJobs);
Expand Down
97 changes: 43 additions & 54 deletions src/main/java/algorithms/AutoThresholdRegression.java
Expand Up @@ -6,7 +6,6 @@
import net.imglib2.TwinCursor;
import net.imglib2.type.logic.BitType;
import net.imglib2.type.numeric.RealType;
import net.imglib2.util.Util;
import net.imglib2.view.Views;
import results.ResultHandler;

Expand All @@ -15,6 +14,9 @@
* used for Person colocalisation calculation.
*/
public class AutoThresholdRegression<T extends RealType< T >> extends Algorithm<T> {
// Identifiers for choosing which implementation to use
public enum Implementation {Costes, Bisection};
Implementation implementation = Implementation.Bisection;
/* the threshold for y-intercept to y-max to
* raise a warning about it being to high.
*/
Expand All @@ -33,8 +35,13 @@ public class AutoThresholdRegression<T extends RealType< T >> extends Algorithm<
PearsonsCorrelation<T> pearsonsCorrellation;

public AutoThresholdRegression(PearsonsCorrelation<T> pc) {
this(pc, Implementation.Costes);
}

public AutoThresholdRegression(PearsonsCorrelation<T> pc, Implementation impl) {
super("auto threshold regression");
pearsonsCorrellation = pc;
implementation = impl;
}

@Override
Expand Down Expand Up @@ -107,27 +114,15 @@ public void execute(DataContainer<T> container)
final double m = num/denom;
final double b = ch2Mean - m*ch1Mean ;

// initialize some variables relevant for regression
// indicates whether the threshold has been found or not
boolean thresholdFound = false;
// the maximum number of iterations to look for the threshold
final int maxIterations = 100;
// the current iteration
int iteration = 0;
// current and last working threshold
double threshold1, threshold2;
// A stepper that walks thresholds
Stepper stepper;
// to map working thresholds to channels
ChannelMapper mapper;

// let working threshold walk on channel one if the regression line
// leans more towards the abscissa (which represents channel one) for
// positive and negative correlation.
if (m > -1 && m < 1.0) {
// Start at the midpoint of channel one
threshold1 = Math.abs(container.getMaxCh1() +
container.getMinCh1()) * 0.5;
threshold2 = container.getMaxCh1();

// Map working threshold to channel one (because channel one has a
// larger maximum value.
mapper = new ChannelMapper() {
Expand All @@ -142,12 +137,17 @@ public double getCh2Threshold(double t) {
return (t * m) + b;
}
};
} else {
// Start at the midpoint of channel two
threshold1 = Math.abs(container.getMaxCh2() +
container.getMinCh2()) * 0.5;
threshold2 = container.getMaxCh2();
// Select a stepper
if (implementation == Implementation.Bisection) {
// Start at the midpoint of channel one
stepper = new BisectionStepper(
Math.abs(container.getMaxCh1() + container.getMinCh1()) * 0.5,
container.getMaxCh1());
} else {
stepper = new SimpleStepper(container.getMaxCh1());
}

} else {
// Map working threshold to channel two (because channel two has a
// larger maximum value.
mapper = new ChannelMapper() {
Expand All @@ -162,6 +162,15 @@ public double getCh2Threshold(double t) {
return t;
}
};
// Select a stepper
if (implementation == Implementation.Bisection) {
// Start at the midpoint of channel two
stepper = new BisectionStepper(
Math.abs(container.getMaxCh2() + container.getMinCh2()) * 0.5,
container.getMaxCh2());
} else {
stepper = new SimpleStepper(container.getMaxCh2());
}
}

// Min threshold not yet implemented
Expand All @@ -182,53 +191,31 @@ public double getCh2Threshold(double t) {
final double maxVal = dummyT.getMaxValue();

// do regression
while (!thresholdFound && iteration<maxIterations) {
while (!stepper.isFinished()) {
// round ch1 threshold and compute ch2 threshold
ch1ThreshMax = Math.round(mapper.getCh1Threshold(threshold1));
ch2ThreshMax = Math.round(mapper.getCh2Threshold(threshold1));
ch1ThreshMax = Math.round(mapper.getCh1Threshold(stepper.getValue()));
ch2ThreshMax = Math.round(mapper.getCh2Threshold(stepper.getValue()));
/* Make sure we don't get overflow the image type specific threshold variables
* if the image data type doesn't support this value.
*/
thresholdCh1.setReal(clamp(ch1ThreshMax, minVal, maxVal));
thresholdCh2.setReal(clamp(ch2ThreshMax, minVal, maxVal));

// Person's R value
double currentPersonsR = Double.MAX_VALUE;
// indicates if we have actually found a real number
boolean badResult = false;
try {
// do persons calculation within the limits
currentPersonsR = pearsonsCorrellation.calculatePearsons(cursor,
ch1Mean, ch2Mean, thresholdCh1, thresholdCh2, ThresholdMode.Below);
final double currentPersonsR =
pearsonsCorrellation.calculatePearsons(cursor,
ch1Mean, ch2Mean, thresholdCh1, thresholdCh2,
ThresholdMode.Below);
stepper.update(currentPersonsR);
} catch (MissingPreconditionException e) {
/* the exception that could occur is due to numerical
* problems within the pearsons calculation.
*/
badResult = true;
}

/* If the difference between both thresholds is < 1, we consider
* that as reasonable close to abort the regression.
*/
final double thrDiff = Math.abs(threshold1 - threshold2);
if (thrDiff < 1.0)
thresholdFound = true;

// update working thresholds for next iteration
threshold2 = threshold1;
if (badResult || currentPersonsR < 0) {
// we went too far, increase by the absolute half
threshold1 = threshold1 + thrDiff * 0.5;
} else if (currentPersonsR > 0) {
// as long as r > 0 we go half the way down
threshold1 = threshold1 - thrDiff * 0.5;
* problems within the Pearsons calculation. */
stepper.update(Double.NaN);
}

// reset the cursor to reuse it
cursor.reset();

// increment iteration counter
iteration++;
}

/* Store the new results. The lower thresholds are the types
Expand Down Expand Up @@ -269,8 +256,9 @@ public double getCh2Threshold(double t) {

// add warnings if values are below lowest pixel value of images
if ( ch1ThreshMax < container.getMinCh1() || ch2ThreshMax < container.getMinCh2() ) {
addWarning("thresholds too low",
"The auto threshold method could not find a positive threshold.");
String msg ="The auto threshold method could not find a positive threshold.";
msg += implementation == Implementation.Costes ? "" : " Maybe you should try classic thresholding.";
addWarning("thresholds too low", msg);
}
}

Expand All @@ -291,6 +279,7 @@ public void processResults(ResultHandler<T> handler) {
handler.handleValue( "b to y-max ratio", bToYMaxRatio, 2 );
handler.handleValue( "Ch1 Max Threshold", ch1MaxThreshold.getRealDouble(), 2);
handler.handleValue( "Ch2 Max Threshold", ch2MaxThreshold.getRealDouble(), 2);
handler.handleValue( "Threshold regression", implementation.toString());
}

public double getBToYMaxRatio() {
Expand Down
70 changes: 70 additions & 0 deletions src/main/java/algorithms/BisectionStepper.java
@@ -0,0 +1,70 @@
package algorithms;

/**
* Try to converge a threshold based on an update value condition. If the
* update value is larger zero, the threshold is lowered by half the distance
* between the last thresholds. If the update value falls below zero or is not
* a number, the threshold is increased by such a half step.
*
* @author Tom Kazimiers
*
*/
public class BisectionStepper extends Stepper {
protected double threshold1;
protected double threshold2;
protected double thrDiff = Double.NaN;
protected int iterations = 0;
protected int maxIterations = 100;

/**
* Initialize the bisection stepper with a start threshold and its
* last threshold
*
* @param threshold The current threshold
* @param lastThreshold The last threshold
*/
public BisectionStepper(double threshold, double lastThreshold) {
threshold1 = threshold;
threshold2 = lastThreshold;
thrDiff = Math.abs(threshold1 - threshold2);
}

/**
* Update threshold by a bisection step. If {@value} is below zero or not
* a number, the step is made upwards. If it is above zero, the stoep is
* downwards.
*/
@Override
public void update(double value) {
// update working thresholds for next iteration
threshold2 = threshold1;
if (Double.NaN == value || value < 0) {
// we went too far, increase by the absolute half
threshold1 = threshold1 + thrDiff * 0.5;
} else if (value > 0) {
// as long as r > 0 we go half the way down
threshold1 = threshold1 - thrDiff * 0.5;
}
// Update difference to last threshold
thrDiff = Math.abs(threshold1 - threshold2);
// Update update counter
iterations++;
}

/**
* Get current threshold.
*/
@Override
public double getValue() {
return threshold1;
}

/**
* If the difference between both thresholds is < 1, we consider
* that as reasonable close to abort the regression.
* */
@Override
public boolean isFinished() {
return iterations > maxIterations || thrDiff < 1.0;
}
}
1 change: 0 additions & 1 deletion src/main/java/algorithms/CostesSignificanceTest.java
Expand Up @@ -22,7 +22,6 @@
import net.imglib2.type.NativeType;
import net.imglib2.type.logic.BitType;
import net.imglib2.type.numeric.RealType;
import net.imglib2.util.Util;
import net.imglib2.view.Views;
import results.ResultHandler;

Expand Down

0 comments on commit 5665600

Please sign in to comment.