use of gdsc.core.ij.IJLogger in project GDSC-SMLM by aherbert.
the class PeakFit method updateFitConfiguration.
/**
* Updates the configuration for peak fitting. Configures the calculation of residuals, logging and peak validation.
*
* @return
*/
private void updateFitConfiguration(FitEngineConfiguration config) {
FitConfiguration fitConfig = config.getFitConfiguration();
// Adjust the settings that are relevant within the fitting configuration.
fitConfig.setComputeResiduals(config.getResidualsThreshold() < 1);
logger = (resultsSettings.logProgress) ? new IJLogger() : null;
fitConfig.setLog(logger);
fitConfig.setComputeDeviations(resultsSettings.showDeviations);
// Add the calibration for precision filtering
fitConfig.setNmPerPixel(calibration.getNmPerPixel());
fitConfig.setGain(calibration.getGain());
fitConfig.setBias(calibration.getBias());
fitConfig.setEmCCD(calibration.isEmCCD());
}
use of gdsc.core.ij.IJLogger in project GDSC-SMLM by aherbert.
the class DoubletAnalysis method runAnalysis.
/**
* Run analysis.
*/
private void runAnalysis() {
if (doubletResults == null) {
IJ.error(TITLE, "No doublet results in memory");
return;
}
// Ask the user to set filters
if (!showAnalysisDialog())
return;
showResults(doubletResults, analysisShowResults);
// Store the effect of fitting as a doublet
ArrayList<DoubletBonus> data = new ArrayList<DoubletBonus>(doubletResults.size());
// True positive and False positives at residuals = 0
double tp = 0;
double fp = 0;
Logger logger = (analysisLogging) ? new IJLogger() : null;
// Get filters for the single and double fits
// No coordinate shift for the doublet. We have already done simple checking of the
// coordinates to get the good=2 flag
FitConfiguration filterFitConfig2 = filterFitConfig.clone();
filterFitConfig2.setCoordinateShift(Integer.MAX_VALUE);
final int size = 2 * config.getRelativeFitting() + 1;
Rectangle regionBounds = new Rectangle(0, 0, size, size);
final double otherDriftAngle = 180 - analysisDriftAngle;
// Process all the results
for (DoubletResult result : doubletResults) {
// Filter the singles that would be accepted
if (result.good1) {
filterFitConfig.setNoise(result.noise);
FitStatus fitStatus0 = filterFitConfig.validatePeak(0, result.fitResult1.getInitialParameters(), result.fitResult1.getParameters());
double tp1 = 0, fp1 = 0;
if (fitStatus0 == FitStatus.OK) {
tp1 = result.tp1;
fp1 = result.fp1;
} else if (analysisLogging)
logFailure(logger, 0, result, fitStatus0);
// width diverged spots as OK for a doublet fit
if ((fitStatus0 == FitStatus.OK || fitStatus0 == FitStatus.WIDTH_DIVERGED) && selectFit(result) && result.good2) {
double tp2 = 0, fp2 = 0;
// Basic spot criteria (SNR, Photons, width)
filterFitConfig2.setNoise(result.noise);
FitStatus fitStatus1 = filterFitConfig2.validatePeak(0, result.fitResult2.getInitialParameters(), result.fitResult2.getParameters());
FitStatus fitStatus2 = filterFitConfig2.validatePeak(1, result.fitResult2.getInitialParameters(), result.fitResult2.getParameters());
// Log basic failures
boolean[] accept = new boolean[2];
if (fitStatus1 == FitStatus.OK) {
accept[0] = true;
} else if (analysisLogging)
logFailure(logger, 1, result, fitStatus1);
if (fitStatus2 == FitStatus.OK) {
accept[1] = true;
} else if (analysisLogging)
logFailure(logger, 2, result, fitStatus2);
// We can filter each spot with criteria such as shift and the angle to the quadrant.
if (accept[0] || accept[1]) {
if (result.gap < minGap) {
accept[0] = accept[1] = false;
if (analysisLogging)
logger.info("Reject Doublet (%.2f): Fitted coordinates below min gap (%g<%g)\n", result.getMaxScore(), result.gap, minGap);
}
}
if (accept[0] || accept[1]) {
// The logic in here will be copied to the FitWorker.quadrantAnalysis routine.
double[] params = result.fitResult1.getParameters();
double[] newParams = result.fitResult2.getParameters();
// Set up for shift filtering
double shift = filterFitConfig.getCoordinateShift();
if (shift == 0 || shift == Double.POSITIVE_INFINITY) {
// Allow the shift to span half of the fitted window.
shift = 0.5 * FastMath.min(regionBounds.width, regionBounds.height);
}
// Set an upper limit on the shift that is not too far outside the fit window
final double maxShiftX, maxShiftY;
final double factor = Gaussian2DFunction.SD_TO_HWHM_FACTOR;
if (fitConfig.isWidth0Fitting()) {
// Add the fitted standard deviation to the allowed shift
maxShiftX = regionBounds.width * 0.5 + factor * params[Gaussian2DFunction.X_SD];
maxShiftY = regionBounds.height * 0.5 + factor * params[Gaussian2DFunction.Y_SD];
} else {
// Add the configured standard deviation to the allowed shift
maxShiftX = regionBounds.width * 0.5 + factor * fitConfig.getInitialPeakStdDev0();
maxShiftY = regionBounds.height * 0.5 + factor * fitConfig.getInitialPeakStdDev1();
}
for (int n = 0; n < 2; n++) {
if (!accept[n])
continue;
// Reset
accept[n] = false;
final double xShift = newParams[Gaussian2DFunction.X_POSITION + n * 6] - params[Gaussian2DFunction.X_POSITION];
final double yShift = newParams[Gaussian2DFunction.Y_POSITION + n * 6] - params[Gaussian2DFunction.Y_POSITION];
if (Math.abs(xShift) > maxShiftX || Math.abs(yShift) > maxShiftY) {
if (analysisLogging)
logger.info("Reject P%d (%.2f): Fitted coordinates moved outside fit region (x=%g,y=%g)\n", n + 1, result.getMaxScore(), xShift, yShift);
continue;
}
if (Math.abs(xShift) > shift || Math.abs(yShift) > shift) {
// Allow up to a 45 degree difference to show the shift is along the vector
if (result.a[n] > analysisDriftAngle && result.a[n] < otherDriftAngle) {
if (analysisLogging)
logger.info("Reject P%d (%.2f): Fitted coordinates moved into wrong quadrant (x=%g,y=%g,a=%f)", n + 1, result.getMaxScore(), xShift, yShift, result.a[n]);
continue;
}
// Note: The FitWorker also checks for drift to another candidate.
}
// This is OK
accept[n] = true;
}
}
if (accept[0]) {
tp2 += result.tp2a;
fp2 += result.fp2a;
}
if (accept[1]) {
tp2 += result.tp2b;
fp2 += result.fp2b;
}
if (accept[0] || accept[1]) {
tp += tp2;
fp += fp2;
// Store this as a doublet bonus
data.add(new DoubletBonus(result.getMaxScore(), result.getAvScore(), tp2 - tp1, fp2 - fp1));
} else {
// No doublet fit so this will always be the single fit result
tp += tp1;
fp += fp1;
}
} else {
// No doublet fit so this will always be the single fit result
tp += tp1;
fp += fp1;
}
}
}
// Compute the max Jaccard
computeScores(data, tp, fp, numberOfMolecules, useMaxResiduals);
if (showJaccardPlot)
plotJaccard(residualsScore, (useMaxResiduals) ? _residualsScoreMax : _residualsScoreAv);
createAnalysisTable();
StringBuilder sb = new StringBuilder(analysisPrefix);
sb.append(analysisTitle).append('\t');
sb.append((useMaxResiduals) ? "Max" : "Average").append('\t');
sb.append(SELECTION_CRITERIA[selectionCriteria]).append('\t');
if (filterFitConfig.isSmartFilter()) {
sb.append(filterFitConfig.getSmartFilterName()).append("\t\t\t\t\t\t\t\t");
} else {
sb.append('\t');
sb.append(filterFitConfig.getCoordinateShiftFactor()).append('\t');
sb.append(filterFitConfig.getSignalStrength()).append('\t');
sb.append(filterFitConfig.getMinPhotons()).append('\t');
sb.append(filterFitConfig.getMinWidthFactor()).append('\t');
sb.append(filterFitConfig.getWidthFactor()).append('\t');
sb.append(filterFitConfig.getPrecisionThreshold()).append('\t');
sb.append(filterFitConfig.isPrecisionUsingBackground()).append('\t');
}
sb.append(analysisDriftAngle).append('\t');
sb.append(minGap).append('\t');
addJaccardScores(sb);
analysisTable.append(sb.toString());
saveTemplate(sb.toString());
}
use of gdsc.core.ij.IJLogger in project GDSC-SMLM by aherbert.
the class TraceDiffusion method fitJumpDistance.
/**
* Fit the jump distance histogram.
* <p>
* Update the plot by adding the fit line(s).
*
* @param jumpDistances
* (in um^2)
* @param jdHistogram
* @return The fitted coefficients and fractions
*/
private double[][] fitJumpDistance(StoredDataStatistics jumpDistances, double[][] jdHistogram) {
final double msd = jumpDistances.getMean();
final double meanDistance = Math.sqrt(msd) * 1e3;
// TODO:
// Q. Should the beta be expressed using the mean-distance or MSD?
// Q. Should it be normalised to the frame length. If not then the beta will be invariant on
// jump distance length
beta = meanDistance / precision;
Utils.log("Jump Distance analysis : N = %d, Time = %d frames (%s seconds). MSD = %s um^2/jump, Mean Distance = %s nm/jump, Precision = %s nm, Beta = %s", jumpDistances.getN(), settings.jumpDistance, Utils.rounded(settings.jumpDistance * exposureTime, 4), Utils.rounded(msd, 4), Utils.rounded(meanDistance, 4), Utils.rounded(precision, 4), Utils.rounded(beta, 4));
IJLogger logger = new IJLogger(debugFitting, debugFitting);
JumpDistanceAnalysis jd = new JumpDistanceAnalysis(logger);
jd.setFitRestarts(settings.fitRestarts);
jd.setMinFraction(minFraction);
jd.setMinDifference(minDifference);
jd.setMinN(myMinN);
jd.setMaxN(maxN);
// Update the plot with the fit
jd.setCurveLogger(this);
// Set the calibration
jd.setN(settings.jumpDistance);
jd.setDeltaT(exposureTime);
if (settings.precisionCorrection)
jd.setError(precision, true);
jd.setMsdCorrection(settings.msdCorrection);
double[][] fit;
if (settings.mle)
fit = jd.fitJumpDistancesMLE(jumpDistances.getValues(), jdHistogram);
else
fit = jd.fitJumpDistanceHistogram(jumpDistances.getMean(), jdHistogram);
// Get the raw fitted D and convert it to a calibrated D*
if (fit != null) {
fit[0] = jd.calculateApparentDiffusionCoefficient(fit[0]);
// Check the largest D
checkTraceDistance(fit[0][0]);
ic = jd.getInformationCriterion();
}
return fit;
}
use of gdsc.core.ij.IJLogger in project GDSC-SMLM by aherbert.
the class PCPALMClusters method fitBinomial.
/**
* Fit a zero-truncated Binomial to the cumulative histogram
*
* @param histogramData
* @return
*/
private double[] fitBinomial(HistogramData histogramData) {
// Get the mean and sum of the input histogram
double mean;
double sum = 0;
count = 0;
for (int i = 0; i < histogramData.histogram[1].length; i++) {
count += histogramData.histogram[1][i];
sum += histogramData.histogram[1][i] * i;
}
mean = sum / count;
String name = "Zero-truncated Binomial distribution";
Utils.log("Mean cluster size = %s", Utils.rounded(mean));
Utils.log("Fitting cumulative " + name);
// Convert to a normalised double array for the binomial fitter
double[] histogram = new double[histogramData.histogram[1].length];
for (int i = 0; i < histogramData.histogram[1].length; i++) histogram[i] = histogramData.histogram[1][i] / count;
// Plot the cumulative histogram
String title = TITLE + " Cumulative Distribution";
Plot2 plot = null;
if (showCumulativeHistogram) {
// Create a cumulative histogram for fitting
double[] cumulativeHistogram = new double[histogram.length];
sum = 0;
for (int i = 0; i < histogram.length; i++) {
sum += histogram[i];
cumulativeHistogram[i] = sum;
}
double[] values = Utils.newArray(histogram.length, 0.0, 1.0);
plot = new Plot2(title, "N", "Cumulative Probability", values, cumulativeHistogram);
plot.setLimits(0, histogram.length - 1, 0, 1.05);
plot.addPoints(values, cumulativeHistogram, Plot2.CIRCLE);
Utils.display(title, plot);
}
// Do fitting for different N
double bestSS = Double.POSITIVE_INFINITY;
double[] parameters = null;
int worse = 0;
int N = histogram.length - 1;
int min = minN;
final boolean customRange = (minN > 1) || (maxN > 0);
if (min > N)
min = N;
if (maxN > 0 && N > maxN)
N = maxN;
Utils.log("Fitting N from %d to %d%s", min, N, (customRange) ? " (custom-range)" : "");
// Since varying the N should be done in integer steps do this
// for n=1,2,3,... until the SS peaks then falls off (is worse then the best
// score several times in succession)
BinomialFitter bf = new BinomialFitter(new IJLogger());
bf.setMaximumLikelihood(maximumLikelihood);
for (int n = min; n <= N; n++) {
PointValuePair solution = bf.fitBinomial(histogram, mean, n, true);
if (solution == null)
continue;
double p = solution.getPointRef()[0];
Utils.log("Fitted %s : N=%d, p=%s. SS=%g", name, n, Utils.rounded(p), solution.getValue());
if (bestSS > solution.getValue()) {
bestSS = solution.getValue();
parameters = new double[] { n, p };
worse = 0;
} else if (bestSS < Double.POSITIVE_INFINITY) {
if (++worse >= 3)
break;
}
if (showCumulativeHistogram)
addToPlot(n, p, title, plot, new Color((float) n / N, 0, 1f - (float) n / N));
}
// Add best it in magenta
if (showCumulativeHistogram && parameters != null)
addToPlot((int) parameters[0], parameters[1], title, plot, Color.magenta);
return parameters;
}
use of gdsc.core.ij.IJLogger in project GDSC-SMLM by aherbert.
the class GaussianFit method createGaussianFitter.
private Gaussian2DFitter createGaussianFitter(boolean simpleFiltering) {
FitConfiguration config = new FitConfiguration();
config.setMaxIterations(getMaxIterations());
config.setSignificantDigits(getSignificantDigits());
config.setDelta(getDelta());
config.setInitialPeakStdDev(getInitialPeakStdDev());
config.setComputeDeviations(showDeviations);
config.setDuplicateDistance(0);
// Set-up peak filtering only for single fitting
config.setDisableSimpleFilter(!simpleFiltering);
setupPeakFiltering(config);
if (isLogProgress()) {
config.setLog(new IJLogger());
}
if (getFitCriteria() >= 0 && getFitCriteria() < FitCriteria.values().length) {
config.setFitCriteria(FitCriteria.values()[getFitCriteria()]);
} else {
config.setFitCriteria(FitCriteria.LEAST_SQUARED_ERROR);
}
if (getFitFunction() >= 0 && getFitFunction() < FitFunction.values().length) {
config.setFitFunction(FitFunction.values()[getFitFunction()]);
} else {
config.setFitFunction(FitFunction.CIRCULAR);
}
config.setBackgroundFitting(fitBackground);
return new Gaussian2DFitter(config);
}
Aggregations