use of gdsc.core.logging.Logger 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.logging.Logger in project GDSC-SMLM by aherbert.
the class JumpDistanceAnalysisTest method fit.
private void fit(String title, int samples, int n, double[] d, double[] f, boolean mle) {
// Used for testing
// @formatter:off
//if (!mle) return;
//if (mle) return;
// For easy mixed populations
//if (f.length == 2 && Math.min(f[0],f[1])/(f[0]+f[1]) <= 0.2) return;
//n = 2;
// @formatter:on
JumpDistanceAnalysis.sort(d, f);
double[] jumpDistances = createData(samples, d, f);
Logger logger = null;
//logger = new gdsc.smlm.utils.logging.ConsoleLogger();
JumpDistanceAnalysis jd = new JumpDistanceAnalysis(logger);
jd.setFitRestarts(3);
double[][] fit;
if (n == 0) {
jd.setMinFraction(0.05);
jd.setMinDifference(2);
jd.setMaxN(10);
fit = (mle) ? jd.fitJumpDistancesMLE(jumpDistances) : jd.fitJumpDistances(jumpDistances);
} else {
// No validation
jd.setMinFraction(0);
jd.setMinDifference(0);
fit = (mle) ? jd.fitJumpDistancesMLE(jumpDistances, n) : jd.fitJumpDistances(jumpDistances, n);
}
double[] fitD = (fit == null) ? new double[0] : fit[0];
double[] fitF = (fit == null) ? new double[0] : fit[1];
// Record results to file
if (out != null)
writeResult(title, sample.d, sample.f, samples, n, d, f, mle, fitD, fitF);
fitN = fitD.length;
AssertionError error = null;
try {
Assert.assertEquals("Failed to fit n", d.length, fitD.length);
for (int i = 0; i < d.length; i++) {
Assert.assertEquals("Failed to fit d", d[i], fitD[i], deltaD * d[i]);
Assert.assertEquals("Failed to fit f", f[i], fitF[i], deltaF * f[i]);
}
} catch (AssertionError e) {
error = e;
} finally {
double[] e1 = getPercentError(d, fitD);
double[] e2 = getPercentError(f, fitF);
log("%s %s N=%d sample=%d, n=%d : %s = %s [%s] : %s = %s [%s]\n", (error == null) ? "+++ Pass" : "--- Fail", title, d.length, samples, n, toString(d), toString(fitD), toString(e1), toString(f), toString(fitF), toString(e2));
if (error != null)
throw error;
}
}
Aggregations