Search in sources :

Example 6 with Nullable

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

the class ResultsMatchCalculator method getRoiCoordinates.

@Nullable
private static double[] getRoiCoordinates(String line) {
    // Extract the startT and x,y coordinates from the first available point in the pair
    final int[] index = { 1, 4 };
    final String[] fields = line.split("\t");
    final int startT = Integer.parseInt(fields[0]);
    for (final int i : index) {
        if (i < fields.length) {
            if (fields[i].equals("-")) {
                continue;
            }
            final double x = Double.parseDouble(fields[i]);
            final double y = Double.parseDouble(fields[i + 1]);
            return new double[] { startT, x, y };
        }
    }
    return null;
}
Also used : PeakResultPoint(uk.ac.sussex.gdsc.smlm.results.PeakResultPoint) Point(java.awt.Point) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 7 with Nullable

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

the class JumpDistanceAnalysis method doFitJumpDistancesMle.

/**
 * Fit the jump distances using a maximum likelihood estimation with the given number of species.
 *
 * <p>Results are sorted by the diffusion coefficient ascending.
 *
 * @param jumpDistances The jump distances (in um^2)
 * @param estimatedD The estimated diffusion coefficient
 * @param n The number of species in the mixed population
 * @return Array containing: { D (um^2), Fractions }. Can be null if no fit was made.
 */
@Nullable
private double[][] doFitJumpDistancesMle(double[] jumpDistances, double estimatedD, int n) {
    final MaxEval maxEval = new MaxEval(20000);
    final CustomPowellOptimizer powellOptimizer = createCustomPowellOptimizer();
    calibrated = isCalibrated();
    if (n == 1) {
        try {
            final JumpDistanceFunction function = new JumpDistanceFunction(jumpDistances, estimatedD);
            // The Powell algorithm can use more general bounds: 0 - Infinity
            final SimpleBounds bounds = new SimpleBounds(function.getLowerBounds(), function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY));
            final PointValuePair solution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(function.guess()), bounds, new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MAXIMIZE);
            final double[] fitParams = solution.getPointRef();
            ll = solution.getValue();
            lastFitValue = fitValue = MathUtils.getAkaikeInformationCriterion(ll, 1);
            final double[] coefficients = fitParams;
            final double[] fractions = new double[] { 1 };
            LoggerUtils.log(logger, Level.INFO, "Fit Jump distance (N=1) : %s, MLE = %s, Akaike IC = %s (%d evaluations)", formatD(fitParams[0]), MathUtils.rounded(ll, 4), MathUtils.rounded(fitValue, 4), powellOptimizer.getEvaluations());
            return new double[][] { coefficients, fractions };
        } catch (final TooManyEvaluationsException ex) {
            LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : Too many evaluation (%d)", powellOptimizer.getEvaluations());
        } catch (final TooManyIterationsException ex) {
            LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : Too many iterations (%d)", powellOptimizer.getIterations());
        } catch (final ConvergenceException ex) {
            LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : %s", ex.getMessage());
        }
        return null;
    }
    final MixedJumpDistanceFunction function = new MixedJumpDistanceFunction(jumpDistances, estimatedD, n);
    final double[] lB = function.getLowerBounds();
    int evaluations = 0;
    PointValuePair constrainedSolution = null;
    try {
        // The Powell algorithm can use more general bounds: 0 - Infinity
        constrainedSolution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(function.guess()), new SimpleBounds(lB, function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY)), new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MAXIMIZE);
        evaluations = powellOptimizer.getEvaluations();
        LoggerUtils.log(logger, Level.FINE, "Powell optimiser fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
    } catch (final TooManyEvaluationsException ex) {
        LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : Too many evaluation (%d)", n, powellOptimizer.getEvaluations());
    } catch (final TooManyIterationsException ex) {
        LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : Too many iterations (%d)", n, powellOptimizer.getIterations());
    } catch (final ConvergenceException ex) {
        LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : %s", n, ex.getMessage());
    }
    if (constrainedSolution == null) {
        LoggerUtils.log(logger, Level.INFO, "Trying CMAES optimiser with restarts ...");
        final double[] uB = function.getUpperBounds();
        final SimpleBounds bounds = new SimpleBounds(lB, uB);
        // Try a bounded CMAES optimiser since the Powell optimiser appears to be
        // sensitive to the order of the parameters. It is not good when the fast particle
        // is the minority fraction. Could this be due to too low an upper bound?
        // The sigma determines the search range for the variables. It should be 1/3 of the initial
        // search region.
        final double[] s = new double[lB.length];
        for (int i = 0; i < s.length; i++) {
            s[i] = (uB[i] - lB[i]) / 3;
        }
        final OptimizationData sigma = new CMAESOptimizer.Sigma(s);
        final OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(function.x.length))));
        // Iterate this for stability in the initial guess
        final CMAESOptimizer cmaesOptimizer = createCmaesOptimizer();
        for (int i = 0; i <= fitRestarts; i++) {
            // Try from the initial guess
            try {
                final PointValuePair solution = cmaesOptimizer.optimize(new InitialGuess(function.guess()), new ObjectiveFunction(function), GoalType.MAXIMIZE, bounds, sigma, popSize, maxEval);
                if (constrainedSolution == null || solution.getValue() > constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    LoggerUtils.log(logger, Level.FINE, "CMAES optimiser [%da] fit (N=%d) : MLE = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
                }
            } catch (final TooManyEvaluationsException ex) {
            // No solution
            }
            if (constrainedSolution == null) {
                continue;
            }
            // Try from the current optimum
            try {
                final PointValuePair solution = cmaesOptimizer.optimize(new InitialGuess(constrainedSolution.getPointRef()), new ObjectiveFunction(function), GoalType.MAXIMIZE, bounds, sigma, popSize, maxEval);
                if (solution.getValue() > constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    LoggerUtils.log(logger, Level.FINE, "CMAES optimiser [%db] fit (N=%d) : MLE = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
                }
            } catch (final TooManyEvaluationsException ex) {
            // No solution
            }
        }
        if (constrainedSolution != null) {
            try {
                // Re-optimise with Powell?
                final PointValuePair solution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(constrainedSolution.getPointRef()), new SimpleBounds(lB, function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY)), new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MAXIMIZE);
                if (solution.getValue() > constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    LoggerUtils.log(logger, Level.INFO, "Powell optimiser re-fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
                }
            } catch (final TooManyEvaluationsException ex) {
            // No solution
            } catch (final TooManyIterationsException ex) {
            // No solution
            } catch (final ConvergenceException ex) {
            // No solution
            }
        }
    }
    if (constrainedSolution == null) {
        LoggerUtils.log(logger, Level.INFO, "Failed to fit N=%d", n);
        return null;
    }
    final double[] fitParams = constrainedSolution.getPointRef();
    ll = constrainedSolution.getValue();
    // Since the fractions must sum to one we subtract 1 degree of freedom from the number of
    // parameters
    fitValue = MathUtils.getAkaikeInformationCriterion(ll, fitParams.length - 1);
    final double[] d = new double[n];
    final double[] f = new double[n];
    double sum = 0;
    for (int i = 0; i < d.length; i++) {
        f[i] = fitParams[i * 2];
        sum += f[i];
        d[i] = fitParams[i * 2 + 1];
    }
    for (int i = 0; i < f.length; i++) {
        f[i] /= sum;
    }
    // Sort by coefficient size
    sort(d, f);
    final double[] coefficients = d;
    final double[] fractions = f;
    LoggerUtils.log(logger, Level.INFO, "Fit Jump distance (N=%d) : %s (%s), MLE = %s, Akaike IC = %s (%d evaluations)", n, formatD(d), format(f), MathUtils.rounded(ll, 4), MathUtils.rounded(fitValue, 4), evaluations);
    if (isValid(d, f)) {
        lastFitValue = fitValue;
        return new double[][] { coefficients, fractions };
    }
    return null;
}
Also used : MaxEval(org.apache.commons.math3.optim.MaxEval) InitialGuess(org.apache.commons.math3.optim.InitialGuess) SimpleBounds(org.apache.commons.math3.optim.SimpleBounds) CMAESOptimizer(org.apache.commons.math3.optim.nonlinear.scalar.noderiv.CMAESOptimizer) ObjectiveFunction(org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction) PointValuePair(org.apache.commons.math3.optim.PointValuePair) TooManyEvaluationsException(org.apache.commons.math3.exception.TooManyEvaluationsException) ConvergenceException(org.apache.commons.math3.exception.ConvergenceException) TooManyIterationsException(org.apache.commons.math3.exception.TooManyIterationsException) OptimizationData(org.apache.commons.math3.optim.OptimizationData) CustomPowellOptimizer(uk.ac.sussex.gdsc.smlm.math3.optim.nonlinear.scalar.noderiv.CustomPowellOptimizer) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 8 with Nullable

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

the class GaussianFit method fitSingle.

/**
 * Fits a 2D Gaussian to the given data. Fits all the specified peaks.
 *
 * <p>Data must be arranged in yx block order, i.e. height rows of width.
 *
 * @param gf the gf
 * @param data the data
 * @param width the width
 * @param height the height
 * @param index Index of the data to fit
 * @param estimatedHeight Estimated height for the peak (input from smoothed data)
 * @return Array containing the fitted curve data: The first value is the Background. The
 *         remaining values are Amplitude, PosX, PosY, StdDevX, StdDevY for each fitted peak. Null
 *         if no fit is possible.
 */
@Nullable
private double[] fitSingle(Gaussian2DFitter gf, float[] data, int width, int height, int index, double estimatedHeight) {
    this.fitResult = gf.fit(SimpleArrayUtils.toDouble(data), width, height, new int[] { index }, new double[] { estimatedHeight });
    if (fitResult.getStatus() == FitStatus.OK) {
        chiSquared = fitResult.getError();
        final double[] params = fitResult.getParameters();
        convertParameters(params);
        // Check the fit is within the data
        if (params[Gaussian2DFunction.X_POSITION] < 0 || params[Gaussian2DFunction.X_POSITION] > width || params[Gaussian2DFunction.Y_POSITION] < 0 || params[Gaussian2DFunction.Y_POSITION] > height) {
            fitResult = new FitResult(FitStatus.OUTSIDE_FIT_REGION, fitResult.getDegreesOfFreedom(), fitResult.getError(), fitResult.getInitialParameters(), fitResult.getParameters(), fitResult.getParameterDeviations(), fitResult.getNumberOfPeaks(), fitResult.getNumberOfFittedParameters(), fitResult.getStatusData(), fitResult.getIterations(), fitResult.getEvaluations());
            return null;
        }
        return params;
    }
    return null;
}
Also used : FitResult(uk.ac.sussex.gdsc.smlm.fitting.FitResult) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 9 with Nullable

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

the class ImageJ3DResultsViewer method createSphereSizeFromPrecision.

@Nullable
private static Point3f[] createSphereSizeFromPrecision(MemoryPeakResults results) {
    final PrecisionResultProcedure p = new PrecisionResultProcedure(results);
    try {
        final PrecisionMethod m = p.getPrecision();
        IJ.log("Using precision method " + FitProtosHelper.getName(m));
        final Point3f[] size = new Point3f[results.size()];
        for (int i = 0, j = 0; i < p.precisions.length; i++) {
            // Precision is in NM which matches the rendering
            final float v = (float) p.precisions[i];
            size[j++] = new Point3f(v, v, v);
        }
        return size;
    } catch (final DataException ex) {
        IJ.error(TITLE, "The results have no precision: " + ex.getMessage());
        return null;
    }
}
Also used : DataException(uk.ac.sussex.gdsc.core.data.DataException) Point3f(org.scijava.vecmath.Point3f) PrecisionResultProcedure(uk.ac.sussex.gdsc.smlm.results.procedures.PrecisionResultProcedure) PrecisionMethod(uk.ac.sussex.gdsc.smlm.data.config.FitProtos.PrecisionMethod) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 10 with Nullable

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

the class ImageJ3DResultsViewer method createAlphaFromSize.

@Nullable
private static float[] createAlphaFromSize(double minA, double maxA, Point3f[] sphereSize) {
    if (sphereSize == null) {
        return null;
    }
    final double[] d = new double[sphereSize.length];
    for (int i = 0; i < d.length; i++) {
        final Point3f p = sphereSize[i];
        if (p.x == p.y && p.y == p.z) {
            d[i] = 3.0 * p.x * p.x;
        } else {
            // Use the squared distance. This is the equivalent of the area of the shape projected to
            // 2D.
            d[i] = (double) p.x * p.x + p.y * p.y + p.z * p.z;
        }
    // Use the average radius. This is the equivalent of the mean radius of an enclosing
    // ellipsoid.
    // d[i] = Math.sqrt(d[i] / 3);
    }
    final double[] limits = MathUtils.limits(d);
    final double min = limits[0];
    final double max = limits[1];
    if (min == max) {
        ImageJUtils.log("No per-item transparency as size is fixed");
        return null;
    }
    final double range = (maxA - minA) / (max - min);
    final float[] alpha = new float[d.length];
    for (int i = 0; i < alpha.length; i++) {
        // Largest distance has lowest alpha (more transparent)
        alpha[i] = (float) (minA + range * (max - d[i]));
    }
    return alpha;
}
Also used : Point3f(org.scijava.vecmath.Point3f) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Aggregations

Nullable (uk.ac.sussex.gdsc.core.annotation.Nullable)83 Point (java.awt.Point)16 LinkedList (java.util.LinkedList)12 ArrayList (java.util.ArrayList)11 Rectangle (java.awt.Rectangle)9 ConvergenceException (org.apache.commons.math3.exception.ConvergenceException)8 TooManyIterationsException (org.apache.commons.math3.exception.TooManyIterationsException)8 FractionalAssignment (uk.ac.sussex.gdsc.core.match.FractionalAssignment)8 DirectFilter (uk.ac.sussex.gdsc.smlm.results.filter.DirectFilter)8 IDirectFilter (uk.ac.sussex.gdsc.smlm.results.filter.IDirectFilter)8 ImagePlus (ij.ImagePlus)7 FloatProcessor (ij.process.FloatProcessor)7 Plot (ij.gui.Plot)6 ConcurrentRuntimeException (org.apache.commons.lang3.concurrent.ConcurrentRuntimeException)6 HistogramPlot (uk.ac.sussex.gdsc.core.ij.HistogramPlot)6 ExtendedGenericDialog (uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog)6 MultiPathFilter (uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter)6 PeakFractionalAssignment (uk.ac.sussex.gdsc.smlm.results.filter.PeakFractionalAssignment)6 ImageProcessor (ij.process.ImageProcessor)5 LeastSquaresBuilder (org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder)5