Search in sources :

Example 86 with Nullable

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

the class DriftCalculator method calculateUsingFrames.

/**
 * Calculates drift using images from N consecutive frames aligned to the overall image.
 *
 * @param results the results
 * @param limits the limits
 * @param reconstructionSize the reconstruction size
 * @return the drift { dx[], dy[] }
 */
@Nullable
private double[][] calculateUsingFrames(MemoryPeakResults results, int[] limits, int reconstructionSize) {
    // Extract the localisations into blocks of N consecutive frames
    final BlockPeakResultProcedure p = new BlockPeakResultProcedure(settings);
    results.sort();
    results.forEach(p);
    final List<List<Localisation>> blocks = p.blocks;
    if (blocks.size() <= 1) {
        tracker.log("ERROR : Require at least 2 images for drift calculation");
        return null;
    }
    // Check the final block has enough localisations
    final List<Localisation> nextBlock = p.nextBlock;
    if (nextBlock.size() < settings.minimimLocalisations) {
        blocks.remove(blocks.size() - 1);
        if (blocks.size() <= 1) {
            tracker.log("ERROR : Require at least 2 images for drift calculation");
            return null;
        }
        final List<Localisation> combinedBlock = blocks.get(blocks.size() - 1);
        combinedBlock.addAll(nextBlock);
    }
    // Find the average time point for each block
    final int[] blockT = new int[blocks.size()];
    int time = 0;
    for (final List<Localisation> block : blocks) {
        long sum = 0;
        for (final Localisation r : block) {
            sum += r.time;
        }
        blockT[time++] = (int) (sum / block.size());
    }
    // Calculate a scale to use when constructing the images for alignment
    final Rectangle bounds = results.getBounds(true);
    final ResultsImageSettings.Builder builder = ResultsImageSettings.newBuilder().setImageSizeMode(ResultsImageSizeMode.IMAGE_SIZE).setImageSize(reconstructionSize);
    final float scale = ImagePeakResultsFactory.getScale(builder, bounds, 1);
    executor = Executors.newFixedThreadPool(Prefs.getThreads());
    final double[] dx = new double[limits[1] + 1];
    final double[] dy = new double[dx.length];
    final double[] originalDriftTimePoints = getOriginalDriftTimePoints(dx, blockT);
    lastdx = null;
    final double smoothing = updateSmoothingParameter(originalDriftTimePoints);
    double change = calculateDriftUsingFrames(blocks, blockT, bounds, scale, dx, dy, originalDriftTimePoints, smoothing, settings.iterations);
    if (Double.isNaN(change) || tracker.isEnded()) {
        return null;
    }
    plotDrift(limits, dx, dy);
    ImageJUtils.log("Drift Calculator : Initial drift " + MathUtils.rounded(change));
    for (int i = 1; i <= settings.maxIterations; i++) {
        change = calculateDriftUsingFrames(blocks, blockT, bounds, scale, dx, dy, originalDriftTimePoints, smoothing, settings.iterations);
        if (Double.isNaN(change)) {
            return null;
        }
        plotDrift(limits, dx, dy);
        if (converged(i, change, getTotalDrift(dx, dy, originalDriftTimePoints))) {
            break;
        }
    }
    if (tracker.isEnded()) {
        return null;
    }
    plotDrift(limits, dx, dy);
    return new double[][] { dx, dy };
}
Also used : Rectangle(java.awt.Rectangle) ResultsImageSettings(uk.ac.sussex.gdsc.smlm.data.config.ResultsProtos.ResultsImageSettings) Point(java.awt.Point) List(java.util.List) ArrayList(java.util.ArrayList) LinkedList(java.util.LinkedList) LocalList(uk.ac.sussex.gdsc.core.utils.LocalList) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 87 with Nullable

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

the class BlinkEstimator method fit.

/**
 * Fit the dark time to counts of molecules curve. Only use the first n fitted points.
 *
 * <p>Calculates:<br> N = The number of photoblinking molecules in the sample<br> nBlink = The
 * average number of blinks per flourophore<br> tOff = The off-time
 *
 * @param td The dark time
 * @param ntd The counts of molecules
 * @param numberOfFittedPointsSetting the number of fitted points
 * @param log Write the fitting results to the ImageJ log window
 * @return The fitted parameters [N, nBlink, tOff], or null if no fit was possible
 */
@Nullable
public double[] fit(double[] td, double[] ntd, int numberOfFittedPointsSetting, boolean log) {
    blinkingModel = new BlinkingFunction();
    blinkingModel.setLogging(true);
    for (int i = 0; i < numberOfFittedPointsSetting; i++) {
        blinkingModel.addPoint(td[i], ntd[i]);
    }
    final LevenbergMarquardtOptimizer optimiser = new LevenbergMarquardtOptimizer(INITIAL_STEP_BOUND_FACTOR, COST_RELATIVE_TOLERANCE, PAR_RELATIVE_TOLERANCE, ORTHO_TOLERANCE, THRESHOLD);
    try {
        final double[] obs = blinkingModel.getY();
        // @formatter:off
        final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(1000).start(new double[] { ntd[0], 0.1, td[1] }).target(obs).weight(new DiagonalMatrix(blinkingModel.getWeights())).model(blinkingModel, blinkingModel::jacobian).build();
        // @formatter:on
        blinkingModel.setLogging(false);
        final Optimum optimum = optimiser.optimize(problem);
        final double[] parameters = optimum.getPoint().toArray();
        double mean = 0;
        for (final double d : obs) {
            mean += d;
        }
        mean /= obs.length;
        double ssTotal = 0;
        for (final double o : obs) {
            ssTotal += (o - mean) * (o - mean);
        }
        // This is true if the weights are 1
        final double ssResiduals = optimum.getResiduals().dotProduct(optimum.getResiduals());
        r2 = 1 - ssResiduals / ssTotal;
        adjustedR2 = getAdjustedCoefficientOfDetermination(ssResiduals, ssTotal, obs.length, parameters.length);
        if (log) {
            ImageJUtils.log("  Fit %d points. R^2 = %s. Adjusted R^2 = %s", obs.length, MathUtils.rounded(r2, 4), MathUtils.rounded(adjustedR2, 4));
            ImageJUtils.log("  N=%s, nBlink=%s, tOff=%s (%s frames)", MathUtils.rounded(parameters[0], 4), MathUtils.rounded(parameters[1], 4), MathUtils.rounded(parameters[2], 4), MathUtils.rounded(parameters[2] / msPerFrame, 4));
        }
        return parameters;
    } catch (final TooManyIterationsException ex) {
        if (log) {
            ImageJUtils.log("  Failed to fit %d points: Too many iterations: (%s)", blinkingModel.size(), ex.getMessage());
        }
    } catch (final ConvergenceException ex) {
        if (log) {
            ImageJUtils.log("  Failed to fit %d points", blinkingModel.size());
        }
    }
    return null;
}
Also used : 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) LeastSquaresBuilder(org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Aggregations

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