Search in sources :

Example 1 with Nullable

use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.

the class PcPalmFitting method fitClusteredModel.

/**
 * Fits the correlation curve with r>0 to the clustered model using the estimated density and
 * precision. Parameters must be fit within a tolerance of the starting values.
 *
 * @param gr the correlation curve
 * @param sigmaS The estimated precision
 * @param proteinDensity The estimated protein density
 * @param resultColour the result colour
 * @return The fitted parameters [precision, density, clusterRadius, clusterDensity]
 */
@Nullable
private double[] fitClusteredModel(double[][] gr, double sigmaS, double proteinDensity, String resultColour) {
    final ClusteredModelFunctionGradient function = new ClusteredModelFunctionGradient();
    clusteredModel = function;
    ImageJUtils.log("Fitting %s: Estimated precision = %f nm, estimated protein density = %g um^-2", clusteredModel.getName(), sigmaS, proteinDensity * 1e6);
    clusteredModel.setLogging(true);
    for (int i = offset; i < gr[0].length; i++) {
        // error)
        if (gr[0][i] > sigmaS * settings.fitAboveEstimatedPrecision) {
            clusteredModel.addPoint(gr[0][i], gr[1][i]);
        }
    }
    double[] parameters;
    // The model is: sigma, density, range, amplitude
    final double[] initialSolution = new double[] { sigmaS, proteinDensity, sigmaS * 5, 1 };
    // Constrain the fitting to be close to the estimated precision (sigmaS) and protein density.
    // LVM fitting does not support constrained fitting so use a bounded optimiser.
    final SumOfSquaresModelFunction clusteredModelMulti = new SumOfSquaresModelFunction(clusteredModel);
    final double[] x = clusteredModelMulti.x;
    // Put some bounds around the initial guess. Use the fitting tolerance (in %) if provided.
    final double limit = (settings.fittingTolerance > 0) ? 1 + settings.fittingTolerance / 100 : 2;
    final double[] lB = new double[] { initialSolution[0] / limit, initialSolution[1] / limit, 0, 0 };
    // The amplitude and range should not extend beyond the limits of the g(r) curve.
    final double[] uB = new double[] { initialSolution[0] * limit, initialSolution[1] * limit, MathUtils.max(x), MathUtils.max(gr[1]) };
    ImageJUtils.log("Fitting %s using a bounded search: %s < precision < %s & %s < density < %s", clusteredModel.getName(), MathUtils.rounded(lB[0], 4), MathUtils.rounded(uB[0], 4), MathUtils.rounded(lB[1] * 1e6, 4), MathUtils.rounded(uB[1] * 1e6, 4));
    final PointValuePair constrainedSolution = runBoundedOptimiser(initialSolution, lB, uB, clusteredModelMulti);
    if (constrainedSolution == null) {
        return null;
    }
    parameters = constrainedSolution.getPointRef();
    int evaluations = boundedEvaluations;
    // Refit using a LVM
    if (settings.refitWithGradients) {
        ImageJUtils.log("Re-fitting %s using a gradient optimisation", clusteredModel.getName());
        final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
        Optimum lvmSolution;
        try {
            // @formatter:off
            final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(parameters).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, new MultivariateMatrixFunction() {

                @Override
                public double[][] value(double[] point) {
                    return function.jacobian(point);
                }
            }).build();
            // @formatter:on
            lvmSolution = optimizer.optimize(problem);
            evaluations += lvmSolution.getEvaluations();
            final double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
            if (ss < constrainedSolution.getValue()) {
                ImageJUtils.log("Re-fitting %s improved the SS from %s to %s (-%s%%)", clusteredModel.getName(), MathUtils.rounded(constrainedSolution.getValue(), 4), MathUtils.rounded(ss, 4), MathUtils.rounded(100 * (constrainedSolution.getValue() - ss) / constrainedSolution.getValue(), 4));
                parameters = lvmSolution.getPoint().toArray();
            }
        } catch (final TooManyIterationsException ex) {
            ImageJUtils.log("Failed to re-fit %s: Too many iterations (%s)", clusteredModel.getName(), ex.getMessage());
        } catch (final ConvergenceException ex) {
            ImageJUtils.log("Failed to re-fit %s: %s", clusteredModel.getName(), ex.getMessage());
        }
    }
    clusteredModel.setLogging(false);
    // Ensure the width is positive
    parameters[0] = Math.abs(parameters[0]);
    double ss = 0;
    final double[] obs = clusteredModel.getY();
    final double[] exp = clusteredModel.value(parameters);
    for (int i = 0; i < obs.length; i++) {
        ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
    }
    final double totalSumSquares = MathUtils.getTotalSumOfSquares(clusteredModel.getY());
    final double adjustedR2 = MathUtils.getAdjustedCoefficientOfDetermination(ss, totalSumSquares, clusteredModel.size(), parameters.length);
    final double fitSigmaS = parameters[0];
    final double fitProteinDensity = parameters[1];
    // The radius of the cluster domain
    final double domainRadius = parameters[2];
    // The density of the cluster domain
    final double domainDensity = parameters[3];
    // This is from the PC-PALM paper. However that paper fits the g(r)protein exponential convolved
    // in 2D with the g(r)PSF. In this method we have just fit the exponential
    final double nCluster = 2 * domainDensity * Math.PI * domainRadius * domainRadius * fitProteinDensity;
    final double e1 = parameterDrift(sigmaS, fitSigmaS);
    final double e2 = parameterDrift(proteinDensity, fitProteinDensity);
    ImageJUtils.log("  %s fit: SS = %f. Adj.R^2 = %f. %d evaluations", clusteredModel.getName(), ss, adjustedR2, evaluations);
    ImageJUtils.log("  %s parameters:", clusteredModel.getName());
    ImageJUtils.log("    Average precision = %s nm (%s%%)", MathUtils.rounded(fitSigmaS, 4), MathUtils.rounded(e1, 4));
    ImageJUtils.log("    Average protein density = %s um^-2 (%s%%)", MathUtils.rounded(fitProteinDensity * 1e6, 4), MathUtils.rounded(e2, 4));
    ImageJUtils.log("    Domain radius = %s nm", MathUtils.rounded(domainRadius, 4));
    ImageJUtils.log("    Domain density = %s", MathUtils.rounded(domainDensity, 4));
    ImageJUtils.log("    nCluster = %s", MathUtils.rounded(nCluster, 4));
    // Check the fitted parameters are within tolerance of the initial estimates
    valid2 = true;
    if (settings.fittingTolerance > 0 && (Math.abs(e1) > settings.fittingTolerance || Math.abs(e2) > settings.fittingTolerance)) {
        ImageJUtils.log("  Failed to fit %s within tolerance (%s%%): Average precision = %f nm (%s%%)," + " average protein density = %g um^-2 (%s%%)", clusteredModel.getName(), MathUtils.rounded(settings.fittingTolerance, 4), fitSigmaS, MathUtils.rounded(e1, 4), fitProteinDensity * 1e6, MathUtils.rounded(e2, 4));
        valid2 = false;
    }
    // positive
    if (domainRadius < fitSigmaS) {
        ImageJUtils.log("  Failed to fit %s: Domain radius is smaller than the average precision (%s < %s)", clusteredModel.getName(), MathUtils.rounded(domainRadius, 4), MathUtils.rounded(fitSigmaS, 4));
        valid2 = false;
    }
    if (domainDensity < 0) {
        ImageJUtils.log("  Failed to fit %s: Domain density is negative (%s)", clusteredModel.getName(), MathUtils.rounded(domainDensity, 4));
        valid2 = false;
    }
    if (adjustedR2 <= randomModelAdjustedR2) {
        ImageJUtils.log("  Failed to fit %s - Adjusted r^2 has decreased %s%%", clusteredModel.getName(), MathUtils.rounded((100 * (randomModelAdjustedR2 - adjustedR2) / randomModelAdjustedR2), 4));
        valid2 = false;
    }
    addResult(clusteredModel.getName(), resultColour, valid2, fitSigmaS, fitProteinDensity, domainRadius, domainDensity, nCluster, -1, adjustedR2);
    return parameters;
}
Also used : PointValuePair(org.apache.commons.math3.optim.PointValuePair) LeastSquaresBuilder(org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder) Optimum(org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum) LevenbergMarquardtOptimizer(org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer) DiagonalMatrix(org.apache.commons.math3.linear.DiagonalMatrix) ConvergenceException(org.apache.commons.math3.exception.ConvergenceException) TooManyIterationsException(org.apache.commons.math3.exception.TooManyIterationsException) LeastSquaresProblem(org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem) MultivariateMatrixFunction(org.apache.commons.math3.analysis.MultivariateMatrixFunction) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 2 with Nullable

use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.

the class PcPalmMolecules method plot.

@Nullable
private static double[][] plot(DoubleData stats, String label, boolean integerBins) {
    if (integerBins) {
        // The histogram is not need for the return statement
        new HistogramPlotBuilder(TITLE, stats, label).setMinBinWidth(1).show();
        return null;
    }
    // Show a cumulative histogram so that the bin size is not relevant
    final double[][] hist = MathUtils.cumulativeHistogram(stats.values(), false);
    // Create the axes
    final double[] xValues = hist[0];
    final double[] yValues = hist[1];
    // Plot
    final String title = TITLE + " " + label;
    final Plot plot = new Plot(title, label, "Frequency");
    plot.addPoints(xValues, yValues, Plot.LINE);
    ImageJUtils.display(title, plot);
    return hist;
}
Also used : Plot(ij.gui.Plot) HistogramPlot(uk.ac.sussex.gdsc.core.ij.HistogramPlot) HistogramPlotBuilder(uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 3 with Nullable

use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.

the class PcPalmMolecules method fitGaussian.

@Nullable
private static double[] fitGaussian(float[] x, float[] y) {
    final WeightedObservedPoints obs = new WeightedObservedPoints();
    for (int i = 0; i < x.length; i++) {
        obs.add(x[i], y[i]);
    }
    final Collection<WeightedObservedPoint> observations = obs.toList();
    final GaussianCurveFitter fitter = GaussianCurveFitter.create().withMaxIterations(2000);
    final GaussianCurveFitter.ParameterGuesser guess = new GaussianCurveFitter.ParameterGuesser(observations);
    double[] initialGuess = null;
    try {
        initialGuess = guess.guess();
        return fitter.withStartPoint(initialGuess).fit(observations);
    } catch (final TooManyEvaluationsException ex) {
        // Use the initial estimate
        return initialGuess;
    } catch (final Exception ex) {
        // Just in case there is another exception type, or the initial estimate failed
        return null;
    }
}
Also used : WeightedObservedPoints(org.apache.commons.math3.fitting.WeightedObservedPoints) TooManyEvaluationsException(org.apache.commons.math3.exception.TooManyEvaluationsException) WeightedObservedPoint(org.apache.commons.math3.fitting.WeightedObservedPoint) WeightedObservedPoint(org.apache.commons.math3.fitting.WeightedObservedPoint) ClusterPoint(uk.ac.sussex.gdsc.core.clustering.ClusterPoint) GaussianCurveFitter(org.apache.commons.math3.fitting.GaussianCurveFitter) DataException(uk.ac.sussex.gdsc.core.data.DataException) TooManyEvaluationsException(org.apache.commons.math3.exception.TooManyEvaluationsException) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 4 with Nullable

use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.

the class PcPalmAnalysis method computeAutoCorrelationCurveFft.

/**
 * Compute the auto-correlation curve using FFT. Computes the correlation image and then samples
 * the image at radii up to the specified length to get the average correlation at a given radius.
 *
 * @param im the image
 * @param wp the weighted image
 * @param maxRadius the max radius
 * @param nmPerPixel the nm per pixel
 * @param density the density
 * @return { distances[], gr[], gr_se[] }
 */
@Nullable
private double[][] computeAutoCorrelationCurveFft(ImageProcessor im, ImageProcessor wp, int maxRadius, double nmPerPixel, double density) {
    log("Performing FFT correlation");
    final FloatProcessor corrIm = computeAutoCorrelationFft(im);
    final FloatProcessor corrW = computeAutoCorrelationFft(wp);
    if (corrIm == null || corrW == null) {
        error("Unable to perform Fourier transform");
        return null;
    }
    final int centre = corrIm.getHeight() / 2;
    final Rectangle crop = new Rectangle(centre - maxRadius, centre - maxRadius, maxRadius * 2, maxRadius * 2);
    if (settings.showCorrelationImages) {
        displayCorrelation(corrIm, "Image correlation FFT", crop);
        displayCorrelation(corrW, "Window correlation FFT", crop);
    }
    log("  Normalising correlation");
    final FloatProcessor correlation = normaliseCorrelation(corrIm, corrW, density);
    if (settings.showCorrelationImages) {
        displayCorrelation(correlation, "Normalised correlation FFT", crop);
    }
    return computeRadialAverage(maxRadius, nmPerPixel, correlation);
}
Also used : FloatProcessor(ij.process.FloatProcessor) Rectangle(java.awt.Rectangle) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 5 with Nullable

use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.

the class TraceMatchCalculator method extractPulses.

@Nullable
private static Pulse[] extractPulses(MemoryPeakResults results) {
    if (results == null) {
        return null;
    }
    final Pulse[] pulses = new Pulse[results.size()];
    final Counter i = new Counter();
    results.forEach(DistanceUnit.PIXEL, (XyrResultProcedure) (x, y, result) -> pulses[i.getAndIncrement()] = new Pulse(x, y, result.getFrame(), result.getEndFrame()));
    return pulses;
}
Also used : CoordinateProvider(uk.ac.sussex.gdsc.smlm.utils.CoordinateProvider) TextWindow(ij.text.TextWindow) AtomicBoolean(java.util.concurrent.atomic.AtomicBoolean) WindowManager(ij.WindowManager) Point(java.awt.Point) HashMap(java.util.HashMap) ImageRoiPainter(uk.ac.sussex.gdsc.smlm.ij.utils.ImageRoiPainter) AtomicReference(java.util.concurrent.atomic.AtomicReference) MatchResult(uk.ac.sussex.gdsc.core.match.MatchResult) ArrayList(java.util.ArrayList) MemoryPeakResults(uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults) MathUtils(uk.ac.sussex.gdsc.core.utils.MathUtils) LinkedList(java.util.LinkedList) XyrResultProcedure(uk.ac.sussex.gdsc.smlm.results.procedures.XyrResultProcedure) ExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog) InputSource(uk.ac.sussex.gdsc.smlm.ij.plugins.ResultsManager.InputSource) Pulse(uk.ac.sussex.gdsc.core.match.Pulse) DistanceUnit(uk.ac.sussex.gdsc.smlm.data.config.UnitProtos.DistanceUnit) Coordinate(uk.ac.sussex.gdsc.core.match.Coordinate) WindowAdapter(java.awt.event.WindowAdapter) ConcurrencyUtils(uk.ac.sussex.gdsc.core.utils.concurrent.ConcurrencyUtils) WindowEvent(java.awt.event.WindowEvent) Consumer(java.util.function.Consumer) List(java.util.List) Counter(uk.ac.sussex.gdsc.smlm.results.count.Counter) PointPair(uk.ac.sussex.gdsc.core.match.PointPair) ImageJUtils(uk.ac.sussex.gdsc.core.ij.ImageJUtils) IJ(ij.IJ) MatchCalculator(uk.ac.sussex.gdsc.core.match.MatchCalculator) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable) PlugIn(ij.plugin.PlugIn) Collections(java.util.Collections) Counter(uk.ac.sussex.gdsc.smlm.results.count.Counter) Pulse(uk.ac.sussex.gdsc.core.match.Pulse) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Aggregations

Nullable (uk.ac.sussex.gdsc.core.annotation.Nullable)34 Point (java.awt.Point)6 LinkedList (java.util.LinkedList)5 Rectangle (java.awt.Rectangle)4 ConvergenceException (org.apache.commons.math3.exception.ConvergenceException)4 TooManyIterationsException (org.apache.commons.math3.exception.TooManyIterationsException)4 FractionalAssignment (uk.ac.sussex.gdsc.core.match.FractionalAssignment)4 DirectFilter (uk.ac.sussex.gdsc.smlm.results.filter.DirectFilter)4 IDirectFilter (uk.ac.sussex.gdsc.smlm.results.filter.IDirectFilter)4 Plot (ij.gui.Plot)3 FloatProcessor (ij.process.FloatProcessor)3 ArrayList (java.util.ArrayList)3 ConcurrentRuntimeException (org.apache.commons.lang3.concurrent.ConcurrentRuntimeException)3 Point3f (org.scijava.vecmath.Point3f)3 ExtendedGenericDialog (uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog)3 Ticker (uk.ac.sussex.gdsc.core.logging.Ticker)3 ImagePlus (ij.ImagePlus)2 List (java.util.List)2 ArrayBlockingQueue (java.util.concurrent.ArrayBlockingQueue)2 TooManyEvaluationsException (org.apache.commons.math3.exception.TooManyEvaluationsException)2