Examples with FitStatus uk.ac.sussex.gdsc.smlm.fitting.FitStatus used on opensource projects

Example 1 with FitStatus

use of uk.ac.sussex.gdsc.smlm.fitting.FitStatus in project GDSC-SMLM by aherbert.

the class DoubletAnalysis method runAnalysis.

/**
 * Run analysis.
 */
private void runAnalysis() {
    final ReferenceResults results = referenceResults.get();
    if (results == null) {
        IJ.error(TITLE, "No doublet results in memory");
        return;
    }
    // Ask the user to set filters
    if (!showAnalysisDialog()) {
        return;
    }
    showResults(results.doubletResults, settings.analysisShowResults);
    // Store the effect of fitting as a doublet
    final ArrayList<DoubletBonus> data = new ArrayList<>(results.doubletResults.size());
    // True positive and False positives at residuals = 0
    double tp = 0;
    double fp = 0;
    final Logger logger = (settings.analysisLogging) ? ImageJPluginLoggerHelper.getLogger(getClass()) : 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
    final FitConfiguration filterFitConfig2 = filterFitConfig.createCopy();
    filterFitConfig2.setCoordinateShift(Integer.MAX_VALUE);
    final int size = 2 * config.getFittingWidth() + 1;
    final Rectangle regionBounds = new Rectangle(0, 0, size, size);
    final double otherDriftAngle = 180 - settings.analysisDriftAngle;
    final FitConfiguration fitConfig = config.getFitConfiguration();
    filterFitConfig.setup();
    // Process all the results
    for (final DoubletResult result : results.doubletResults) {
        // Filter the singles that would be accepted
        if (result.good1) {
            filterFitConfig.setNoise(result.noise);
            final FitStatus fitStatus0 = filterFitConfig.validatePeak(0, result.fitResult1.getInitialParameters(), result.fitResult1.getParameters(), result.fitResult1.getParameterDeviations());
            double tp1 = 0;
            double fp1 = 0;
            if (fitStatus0 == FitStatus.OK) {
                tp1 = result.tp1;
                fp1 = result.fp1;
            } else if (settings.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;
                double fp2 = 0;
                // Basic spot criteria (SNR, Photons, width)
                filterFitConfig2.setNoise(result.noise);
                final FitStatus fitStatus1 = filterFitConfig2.validatePeak(0, result.fitResult2.getInitialParameters(), result.fitResult2.getParameters(), result.fitResult2.getParameterDeviations());
                final FitStatus fitStatus2 = filterFitConfig2.validatePeak(1, result.fitResult2.getInitialParameters(), result.fitResult2.getParameters(), result.fitResult2.getParameterDeviations());
                // Log basic failures
                final boolean[] accept = new boolean[2];
                if (fitStatus1 == FitStatus.OK) {
                    accept[0] = true;
                } else if (settings.analysisLogging) {
                    logFailure(logger, 1, result, fitStatus1);
                }
                if (fitStatus2 == FitStatus.OK) {
                    accept[1] = true;
                } else if (settings.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]) && (result.gap < settings.minGap)) {
                    accept[0] = accept[1] = false;
                    if (logger != null) {
                        LoggerUtils.log(logger, Level.INFO, "Reject Doublet (%.2f): Fitted coordinates below min gap (%g<%g)", result.getMaxScore(), result.gap, settings.minGap);
                    }
                }
                if (accept[0] || accept[1]) {
                    // The logic in here will be copied to the FitWorker.quadrantAnalysis routine.
                    final double[] params = result.fitResult1.getParameters();
                    final 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 * Math.min(regionBounds.width, regionBounds.height);
                    }
                    // Set an upper limit on the shift that is not too far outside the fit window
                    final double maxShiftX;
                    final double maxShiftY;
                    final double factor = Gaussian2DFunction.SD_TO_HWHM_FACTOR;
                    if (fitConfig.isXSdFitting()) {
                        // 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.getInitialXSd();
                        maxShiftY = regionBounds.height * 0.5 + factor * fitConfig.getInitialYSd();
                    }
                    for (int n = 0; n < 2; n++) {
                        if (!accept[n]) {
                            continue;
                        }
                        // Reset
                        accept[n] = false;
                        final double xShift = newParams[Gaussian2DFunction.X_POSITION + n * Gaussian2DFunction.PARAMETERS_PER_PEAK] - params[Gaussian2DFunction.X_POSITION];
                        final double yShift = newParams[Gaussian2DFunction.Y_POSITION + n * Gaussian2DFunction.PARAMETERS_PER_PEAK] - params[Gaussian2DFunction.Y_POSITION];
                        if (Math.abs(xShift) > maxShiftX || Math.abs(yShift) > maxShiftY) {
                            if (logger != null) {
                                LoggerUtils.log(logger, Level.INFO, "Reject P%d (%.2f): Fitted coordinates moved outside fit region (x=%g,y=%g)", 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
                        (result.angle[n] > settings.analysisDriftAngle && result.angle[n] < otherDriftAngle)) {
                            if (logger != null) {
                                LoggerUtils.log(logger, Level.INFO, "Reject P%d (%.2f): Fitted coordinates moved into wrong quadrant" + " (x=%g,y=%g,a=%f)", n + 1, result.getMaxScore(), xShift, yShift, result.angle[n]);
                            }
                            continue;
                        }
                        // 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, results.numberOfMolecules, settings.useMaxResiduals);
    if (settings.showJaccardPlot) {
        plotJaccard(residualsScore, (settings.useMaxResiduals) ? results.residualsScoreMax : results.residualsScoreAv);
    }
    final StringBuilder sb = new StringBuilder(results.analysisPrefix);
    sb.append(settings.analysisTitle).append('\t');
    sb.append((settings.useMaxResiduals) ? "Max" : "Average").append('\t');
    sb.append(Settings.SELECTION_CRITERIAS[settings.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.getMaxWidthFactor()).append('\t');
        sb.append(filterFitConfig.getPrecisionThreshold()).append('\t');
        sb.append(filterFitConfig.getPrecisionMethod()).append('\t');
    }
    sb.append(settings.analysisDriftAngle).append('\t');
    sb.append(settings.minGap).append('\t');
    addJaccardScores(sb);
    createAnalysisTable().append(sb.toString());
    saveTemplate(sb.toString());
}

Example 2 with FitStatus

use of uk.ac.sussex.gdsc.smlm.fitting.FitStatus in project GDSC-SMLM by aherbert.

the class BaseFunctionSolverTest method fitGaussian.

double[] fitGaussian(FunctionSolver solver, double[] data, double[] params, double[] expected) {
    // logger.info(FunctionUtils.getSupplier("%s : Expected %s", solver.getClass().getSimpleName(),
    // Arrays.toString(expected)));
    params = params.clone();
    final FitStatus status = solver.fit(data, null, params, null);
    if (status != FitStatus.OK) {
        Assertions.fail(FunctionUtils.getSupplier("Fit Failed: %s i=%d: %s != %s", status.toString(), solver.getIterations(), Arrays.toString(params), Arrays.toString(expected)));
    }
    return params;
}

Example 3 with FitStatus

use of uk.ac.sussex.gdsc.smlm.fitting.FitStatus in project GDSC-SMLM by aherbert.

the class BaseFunctionSolver method fit.

@Override
public FitStatus fit(double[] y, double[] fx, double[] a, double[] parametersVariance) {
    // Reset the results
    numberOfFittedPoints = y.length;
    iterations = 0;
    evaluations = 0;
    value = 0;
    lastY = null;
    lastA = null;
    preProcess();
    final FitStatus status = computeFit(y, fx, a, parametersVariance);
    if (status == FitStatus.OK) {
        if (lastY == null) {
            lastY = y;
        }
        if (lastA == null) {
            lastA = a;
        }
        postProcess();
    }
    return status;
}

Example 4 with FitStatus

use of uk.ac.sussex.gdsc.smlm.fitting.FitStatus in project GDSC-SMLM by aherbert.

the class NonLinearFit method computeFit.

/**
 * Uses Levenberg-Marquardt method to fit a nonlinear model with coefficients (a) for a set of
 * data points (x, y).
 *
 * <p>It is assumed that the data points x[i] corresponding to y[i] are consecutive integers from
 * zero.
 *
 * @param y Set of n data points to fit (input)
 * @param fx Fitted data points (output)
 * @param a Set of m coefficients (input/output)
 * @param parametersVariance Standard deviation of the set of m coefficients (output)
 * @return The fit status
 */
@Override
public FitStatus computeFit(double[] y, double[] fx, final double[] a, final double[] parametersVariance) {
    final int n = y.length;
    final int nparams = function.gradientIndices().length;
    // Create dynamically for the parameter sizes
    calculator = GradientCalculatorUtils.newCalculator(nparams, isMle());
    // Initialise storage.
    // Note that covar and da are passed to EJMLLinerSolver and so must be the correct size.
    beta = new double[nparams];
    da = new double[nparams];
    covar = new double[nparams][nparams];
    alpha = new double[nparams][nparams];
    ap = new double[a.length];
    // Store the { best, previous, new } sum-of-squares values
    sumOfSquaresWorking = new double[3];
    boolean copyYfit = true;
    if (isMle()) {
        // We must have positive data
        y = ensurePositive(n, y);
        // Store the function values for use in computing the log likelihood
        lastY = y;
        if (fx == null) {
            // Re-use space
            lastFx = SimpleArrayUtils.ensureSize(lastFx, y.length);
            fx = lastFx;
            // We will not need to copy fx later since lastFx is used direct
            copyYfit = false;
        }
    }
    final FitStatus result = doFit(n, y, fx, a, parametersVariance, sc);
    this.evaluations = this.iterations = sc.getIteration();
    // Ensure we have a private copy of fx since the any calling code may modify it
    if (isMle() && copyYfit) {
        lastFx = SimpleArrayUtils.ensureSize(lastFx, y.length);
        System.arraycopy(fx, 0, lastFx, 0, y.length);
    }
    return result;
}