Search in sources :

Example 1 with Logger

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());
}
Also used : ArrayList(java.util.ArrayList) Rectangle(java.awt.Rectangle) Logger(gdsc.core.logging.Logger) IJLogger(gdsc.core.ij.IJLogger) IJLogger(gdsc.core.ij.IJLogger) PeakResultPoint(gdsc.smlm.ij.plugins.ResultsMatchCalculator.PeakResultPoint) BasePoint(gdsc.core.match.BasePoint) FitStatus(gdsc.smlm.fitting.FitStatus) FitConfiguration(gdsc.smlm.fitting.FitConfiguration)

Example 2 with Logger

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;
    }
}
Also used : Logger(gdsc.core.logging.Logger)

Aggregations

Logger (gdsc.core.logging.Logger)2 IJLogger (gdsc.core.ij.IJLogger)1 BasePoint (gdsc.core.match.BasePoint)1 FitConfiguration (gdsc.smlm.fitting.FitConfiguration)1 FitStatus (gdsc.smlm.fitting.FitStatus)1 PeakResultPoint (gdsc.smlm.ij.plugins.ResultsMatchCalculator.PeakResultPoint)1 Rectangle (java.awt.Rectangle)1 ArrayList (java.util.ArrayList)1