Search in sources :

Example 16 with Nullable

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

the class DriftCalculator method calculateUsingMarkers.

/**
 * Calculates drift using the feducial markers within ROI.
 *
 * <p>Adapted from the drift calculation method in QuickPALM.
 *
 * @param results the results
 * @param limits the limits
 * @param rois the rois
 * @return the drift { dx[], dy[] }
 */
@Nullable
private double[][] calculateUsingMarkers(MemoryPeakResults results, int[] limits, Roi[] rois) {
    final Spot[][] roiSpots = findSpots(results, rois, limits);
    // Check we have enough data
    if (roiSpots.length == 0) {
        IJ.error("No peak fit results in the selected ROIs");
        return null;
    }
    final double[] dx = new double[limits[1] + 1];
    final double[] dy = new double[dx.length];
    final double[] sum = new double[roiSpots.length];
    final double[] weights = calculateWeights(roiSpots, dx.length, sum);
    final double smoothing = updateSmoothingParameter(weights);
    lastdx = null;
    double change = calculateDriftUsingMarkers(roiSpots, weights, sum, dx, dy, smoothing, settings.iterations);
    if (Double.isNaN(change) || tracker.isEnded()) {
        return null;
    }
    ImageJUtils.log("Drift Calculator : Initial drift " + MathUtils.rounded(change));
    for (int i = 1; i <= settings.maxIterations; i++) {
        change = calculateDriftUsingMarkers(roiSpots, weights, sum, dx, dy, smoothing, settings.iterations);
        if (Double.isNaN(change)) {
            return null;
        }
        if (converged(i, change, getTotalDrift(dx, dy, weights))) {
            break;
        }
    }
    if (tracker.isEnded()) {
        return null;
    }
    interpolate(dx, dy, weights);
    plotDrift(limits, dx, dy);
    saveDrift(weights, dx, dy);
    return new double[][] { dx, dy };
}
Also used : Point(java.awt.Point) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 17 with Nullable

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

the class MultiPathFilter method acceptAnyDoublet.

/**
 * Check any new and all existing results within the multi-doublet fit results are valid. Returns
 * the new results. Coordinate shift filter is disabled for the doublet results.
 *
 * <p>New results and validated candidates that fail the primary filter can be filtered using the
 * minimal filter and sent to the store. The store can be used to determine if a fit for a
 * different candidate has been performed already.
 *
 * @param multiPathResult the multi path result
 * @param validateCandidates Set to true to validate the candidates
 * @param store the store
 * @param candidateId the candidate id
 * @return The new results that pass the filter
 */
@Nullable
private PreprocessedPeakResult[] acceptAnyDoublet(final MultiPathFitResult multiPathResult, boolean validateCandidates, SelectedResultStore store, final int candidateId) {
    final FitResult multiDoubletFitResult = multiPathResult.getMultiDoubletFitResult();
    if (multiDoubletFitResult == null || multiDoubletFitResult.getResults() == null) {
        return null;
    }
    // Doublets may drift further than single spot candidates.
    // We must validate the doublet spot without shift filtering.
    // Note: Only disable shift for the doublet results.
    // doublets = len(multi-doublet) - len(multi) + 1
    final PreprocessedPeakResult[] results = multiDoubletFitResult.getResults();
    final int nDoublets = results.length - multiPathResult.getMultiFitResult().getResults().length + 1;
    setupFilter(FilterValidationOption.NO_SHIFT);
    validationResults = new int[results.length];
    for (int i = 0; i < nDoublets; i++) {
        validationResults[i] = filter.validate(results[i]);
    }
    restoreFilterState();
    for (int i = nDoublets; i < results.length; i++) {
        validationResults[i] = filter.validate(results[i]);
    }
    return acceptAnyInternal(candidateId, multiDoubletFitResult, validateCandidates, store);
}
Also used : FitResult(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFitResult.FitResult) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 18 with Nullable

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

the class MultiPathFilter method score.

/**
 * Score the assignments (TP/FP) and then clear the list.
 *
 * @param assignments The assignments
 * @param score Scores array to accumulate TP/FP scores
 * @param predicted The number of predictions
 * @param save Set to true to save the scored assignments
 * @param actual The number of actual results in the frame
 * @return the fractional assignments
 */
@Nullable
private static FractionalAssignment[] score(final ArrayList<FractionalAssignment> assignments, final double[] score, final int predicted, boolean save, int actual) {
    if (assignments.isEmpty()) {
        return null;
    }
    final FractionalAssignment[] tmp = assignments.toArray(new FractionalAssignment[0]);
    final RankedScoreCalculator calc = RankedScoreCalculator.create(tmp, actual, predicted);
    final double[] result = calc.score(predicted, false, save);
    score[0] += result[0];
    score[1] += result[1];
    score[2] += result[2];
    score[3] += result[3];
    assignments.clear();
    return calc.getScoredAssignments();
}
Also used : FractionalAssignment(uk.ac.sussex.gdsc.core.match.FractionalAssignment) RankedScoreCalculator(uk.ac.sussex.gdsc.core.match.RankedScoreCalculator) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 19 with Nullable

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

the class FindFociBaseProcessor method findCentreGaussianFit.

/**
 * Finds the centre of the image using a 2D Gaussian fit to projection along the Z-axis.
 *
 * @param subImage the sub image
 * @param dimensions the dimensions
 * @param projectionMethod (0) Average value; (1) Maximum value
 * @return the centre of the image
 */
@Nullable
private int[] findCentreGaussianFit(float[] subImage, int[] dimensions, int projectionMethod) {
    if (FindFoci_PlugIn.IS_GAUSSIAN_FIT_ENABLED < 1) {
        return null;
    }
    final int blockSize = dimensions[0] * dimensions[1];
    final float[] projection = new float[blockSize];
    if (projectionMethod == 1) {
        // Maximum value
        for (int z = dimensions[2]; z-- > 0; ) {
            int index = blockSize * z;
            for (int i = 0; i < blockSize; i++, index++) {
                if (projection[i] < subImage[index]) {
                    projection[i] = subImage[index];
                }
            }
        }
    } else {
        // Average value
        for (int z = dimensions[2]; z-- > 0; ) {
            int index = blockSize * z;
            for (int i = 0; i < blockSize; i++, index++) {
                projection[i] += subImage[index];
            }
        }
        for (int i = blockSize; i-- > 0; ) {
            projection[i] /= dimensions[2];
        }
    }
    final GaussianFit_PlugIn gf = new GaussianFit_PlugIn();
    final double[] fitParams = gf.fit(projection, dimensions[0], dimensions[1]);
    int[] centre = null;
    if (fitParams != null) {
        // Find the centre of mass along the z-axis
        centre = convertCentre(new double[] { fitParams[2], fitParams[3] });
        // Use the centre of mass along the projection axis
        double com = 0;
        long sum = 0;
        for (int z = dimensions[2]; z-- > 0; ) {
            final int index = blockSize * z;
            final float value = subImage[index];
            if (value > 0) {
                com += z * value;
                sum += value;
            }
        }
        if (sum == 0) {
            centre[2] = dimensions[2] / 2;
        } else {
            // Avoid clipping
            centre[2] = MathUtils.clip(0, dimensions[2] - 1, round(com / sum));
        }
    }
    return centre;
}
Also used : GaussianFit_PlugIn(uk.ac.sussex.gdsc.ij.utils.GaussianFit_PlugIn) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 20 with Nullable

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

the class FindFociBaseProcessor method getSortedMaxpoints.

/**
 * Find all local maxima (irrespective whether they finally qualify as maxima or not).
 *
 * @param pixels The image to be analyzed
 * @param maxima the maxima
 * @param types A byte image, same size as ip, where the maximum points are marked as MAXIMUM
 * @param globalMin The image global minimum
 * @param threshold The threshold below which no pixels are processed.
 * @return Maxima sorted by value.
 */
@Nullable
protected Coordinate[] getSortedMaxpoints(Object pixels, int[] maxima, byte[] types, float globalMin, float threshold) {
    final ArrayList<Coordinate> maxpoints = new ArrayList<>(500);
    // working list for expanding local plateaus
    final IntArrayList pointList = new IntArrayList();
    int id = 0;
    final int[] xyz = new int[3];
    setPixels(pixels);
    if (is2D()) {
        for (int i = maxxByMaxyByMaxz; i-- > 0; ) {
            if ((types[i] & (EXCLUDED | MAX_AREA | PLATEAU | NOT_MAXIMUM)) != 0) {
                continue;
            }
            final float v = getf(i);
            if (v < threshold || v == globalMin) {
                continue;
            }
            getXy(i, xyz);
            final int x = xyz[0];
            final int y = xyz[1];
            // Check whether we have a local maximum.
            final boolean isInnerXy = (y != 0 && y != ylimit) && (x != 0 && x != xlimit);
            boolean isMax = true;
            boolean equalNeighbour = false;
            for (int d = 8; d-- > 0; ) {
                if (isInnerXy || isWithinXy(x, y, d)) {
                    final float vNeighbor = getf(i + offset[d]);
                    if (vNeighbor > v) {
                        isMax = false;
                        break;
                    } else if (vNeighbor == v) {
                        // Neighbour is equal, this is a potential plateau maximum
                        equalNeighbour = true;
                    } else {
                        // This is lower so cannot be a maxima
                        types[i + offset[d]] |= NOT_MAXIMUM;
                    }
                }
            }
            if (isMax) {
                id++;
                if (id >= searchCapacity) {
                    log(() -> "The number of potential maxima exceeds the search capacity: " + searchCapacity + ". Try using a denoising/smoothing filter or increase the capacity.");
                    return null;
                }
                if (equalNeighbour) {
                    // Search the local area marking all equal neighbour points as maximum
                    if (!expandMaximum(maxima, types, globalMin, threshold, i, v, id, maxpoints, pointList)) {
                        // Not a true maximum, ignore this
                        id--;
                    }
                } else {
                    types[i] |= MAXIMUM | MAX_AREA;
                    maxima[i] = id;
                    maxpoints.add(new Coordinate(getIndex(x, y), id, v));
                }
            }
        }
    } else {
        for (int i = maxxByMaxyByMaxz; i-- > 0; ) {
            if ((types[i] & (EXCLUDED | MAX_AREA | PLATEAU | NOT_MAXIMUM)) != 0) {
                continue;
            }
            final float v = getf(i);
            if (v < threshold || v == globalMin) {
                continue;
            }
            getXyz(i, xyz);
            final int x = xyz[0];
            final int y = xyz[1];
            final int z = xyz[2];
            // Check whether we have a local maximum.
            final boolean isInnerXy = (y != 0 && y != ylimit) && (x != 0 && x != xlimit);
            final boolean isInnerXyz = (zlimit == 0) ? isInnerXy : isInnerXy && (z != 0 && z != zlimit);
            boolean isMax = true;
            boolean equalNeighbour = false;
            for (int d = 26; d-- > 0; ) {
                if (isInnerXyz || (isInnerXy && isWithinZ(z, d)) || isWithinXyz(x, y, z, d)) {
                    final float vNeighbor = getf(i + offset[d]);
                    if (vNeighbor > v) {
                        isMax = false;
                        break;
                    } else if (vNeighbor == v) {
                        // Neighbour is equal, this is a potential plateau maximum
                        equalNeighbour = true;
                    } else {
                        // This is lower so cannot be a maxima
                        types[i + offset[d]] |= NOT_MAXIMUM;
                    }
                }
            }
            if (isMax) {
                id++;
                if (id >= searchCapacity) {
                    log(() -> "The number of potential maxima exceeds the search capacity: " + searchCapacity + ". Try using a denoising/smoothing filter or increase the capacity.");
                    return null;
                }
                if (equalNeighbour) {
                    // Search the local area marking all equal neighbour points as maximum
                    if (!expandMaximum(maxima, types, globalMin, threshold, i, v, id, maxpoints, pointList)) {
                        // Not a true maximum, ignore this
                        id--;
                    }
                } else {
                    types[i] |= MAXIMUM | MAX_AREA;
                    maxima[i] = id;
                    maxpoints.add(new Coordinate(getIndex(x, y, z), id, v));
                }
            }
        }
    }
    if (ImageJUtils.isInterrupted()) {
        return null;
    }
    for (int i = maxxByMaxyByMaxz; i-- > 0; ) {
        // reset attributes no longer needed
        types[i] &= ~NOT_MAXIMUM;
    }
    Collections.sort(maxpoints, Coordinate::compare);
    // Build a map between the original id and the new id following the sort
    final int[] idMap = new int[maxpoints.size() + 1];
    // Label the points
    for (int i = 0; i < maxpoints.size(); i++) {
        final int newId = (i + 1);
        final Coordinate c = maxpoints.get(i);
        final int oldId = c.id;
        idMap[oldId] = newId;
        c.setId(newId);
    }
    reassignMaxima(maxima, idMap);
    return maxpoints.toArray(new Coordinate[0]);
}
Also used : ArrayList(java.util.ArrayList) IntArrayList(it.unimi.dsi.fastutil.ints.IntArrayList) IntArrayList(it.unimi.dsi.fastutil.ints.IntArrayList) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Aggregations

Nullable (uk.ac.sussex.gdsc.core.annotation.Nullable)49 Point (java.awt.Point)10 ArrayList (java.util.ArrayList)8 LinkedList (java.util.LinkedList)6 ImagePlus (ij.ImagePlus)5 Rectangle (java.awt.Rectangle)5 FloatProcessor (ij.process.FloatProcessor)4 ImageProcessor (ij.process.ImageProcessor)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 File (java.io.File)3 ConcurrentRuntimeException (org.apache.commons.lang3.concurrent.ConcurrentRuntimeException)3 LeastSquaresBuilder (org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder)3 Optimum (org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum)3 LeastSquaresProblem (org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem)3 LevenbergMarquardtOptimizer (org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer)3