Search in sources :

Example 1 with IJLogger

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());
}
Also used : FitConfiguration(gdsc.smlm.fitting.FitConfiguration) IJLogger(gdsc.core.ij.IJLogger)

Example 2 with IJLogger

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());
}
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 3 with IJLogger

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;
}
Also used : JumpDistanceAnalysis(gdsc.smlm.fitting.JumpDistanceAnalysis) IJLogger(gdsc.core.ij.IJLogger)

Example 4 with IJLogger

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;
}
Also used : Color(java.awt.Color) BinomialFitter(gdsc.smlm.fitting.BinomialFitter) Plot2(ij.gui.Plot2) ClusterPoint(gdsc.core.clustering.ClusterPoint) IJLogger(gdsc.core.ij.IJLogger) PointValuePair(org.apache.commons.math3.optim.PointValuePair)

Example 5 with IJLogger

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);
}
Also used : FitConfiguration(gdsc.smlm.fitting.FitConfiguration) Gaussian2DFitter(gdsc.smlm.fitting.Gaussian2DFitter) IJLogger(gdsc.core.ij.IJLogger)

Aggregations

IJLogger (gdsc.core.ij.IJLogger)5 FitConfiguration (gdsc.smlm.fitting.FitConfiguration)3 ClusterPoint (gdsc.core.clustering.ClusterPoint)1 Logger (gdsc.core.logging.Logger)1 BasePoint (gdsc.core.match.BasePoint)1 BinomialFitter (gdsc.smlm.fitting.BinomialFitter)1 FitStatus (gdsc.smlm.fitting.FitStatus)1 Gaussian2DFitter (gdsc.smlm.fitting.Gaussian2DFitter)1 JumpDistanceAnalysis (gdsc.smlm.fitting.JumpDistanceAnalysis)1 PeakResultPoint (gdsc.smlm.ij.plugins.ResultsMatchCalculator.PeakResultPoint)1 Plot2 (ij.gui.Plot2)1 Color (java.awt.Color)1 Rectangle (java.awt.Rectangle)1 ArrayList (java.util.ArrayList)1 PointValuePair (org.apache.commons.math3.optim.PointValuePair)1