Search in sources :

Example 66 with Nullable

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

the class PcPalmFitting method fitRandomModel.

/**
 * Fits the correlation curve with r>0 to the random 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 estimate protein density
 * @param resultColour the result colour
 * @return The fitted parameters [precision, density]
 */
@Nullable
private double[] fitRandomModel(double[][] gr, double sigmaS, double proteinDensity, String resultColour) {
    final RandomModelFunction function = new RandomModelFunction();
    randomModel = function;
    ImageJUtils.log("Fitting %s: Estimated precision = %f nm, estimated protein density = %g um^-2", randomModel.getName(), sigmaS, proteinDensity * 1e6);
    randomModel.setLogging(true);
    for (int i = offset; i < gr[0].length; i++) {
        // error)
        if (gr[0][i] > sigmaS * settings.fitAboveEstimatedPrecision) {
            randomModel.addPoint(gr[0][i], gr[1][i]);
        }
    }
    final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
    Optimum optimum;
    try {
        // @formatter:off
        final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(new double[] { sigmaS, proteinDensity }).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, function::jacobian).build();
        // @formatter:on
        optimum = optimizer.optimize(problem);
    } catch (final TooManyIterationsException ex) {
        ImageJUtils.log("Failed to fit %s: Too many iterations (%s)", randomModel.getName(), ex.getMessage());
        return null;
    } catch (final ConvergenceException ex) {
        ImageJUtils.log("Failed to fit %s: %s", randomModel.getName(), ex.getMessage());
        return null;
    }
    randomModel.setLogging(false);
    final double[] parameters = optimum.getPoint().toArray();
    // Ensure the width is positive
    parameters[0] = Math.abs(parameters[0]);
    final double ss = optimum.getResiduals().dotProduct(optimum.getResiduals());
    final double totalSumSquares = MathUtils.getTotalSumOfSquares(randomModel.getY());
    randomModelAdjustedR2 = MathUtils.getAdjustedCoefficientOfDetermination(ss, totalSumSquares, randomModel.size(), parameters.length);
    final double fitSigmaS = parameters[0];
    final double fitProteinDensity = parameters[1];
    // Check the fitted parameters are within tolerance of the initial estimates
    final double e1 = parameterDrift(sigmaS, fitSigmaS);
    final double e2 = parameterDrift(proteinDensity, fitProteinDensity);
    ImageJUtils.log("  %s fit: SS = %f. Adj.R^2 = %f. %d evaluations", randomModel.getName(), ss, randomModelAdjustedR2, optimum.getEvaluations());
    ImageJUtils.log("  %s parameters:", randomModel.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));
    valid1 = 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%%)", randomModel.getName(), MathUtils.rounded(settings.fittingTolerance, 4), fitSigmaS, MathUtils.rounded(e1, 4), fitProteinDensity * 1e6, MathUtils.rounded(e2, 4));
        valid1 = false;
    }
    if (valid1) {
        // ---------
        // TODO - My data does not comply with this criteria.
        // This could be due to the PC-PALM Molecule code limiting the nmPerPixel to fit the images in
        // memory thus removing correlations at small r.
        // It could also be due to the nature of the random simulations being 3D not 2D membranes
        // as per the PC-PALM paper.
        // ---------
        // Evaluate g(r)protein where:
        // g(r)peaks = g(r)protein + g(r)stoch
        // g(r)peaks ~ 1 + g(r)stoch
        // Verify g(r)protein should be <1.5 for all r>0
        final double[] grStoch = randomModel.value(parameters);
        final double[] grPeaks = randomModel.getY();
        final double[] radius = randomModel.getX();
        for (int i = 0; i < grPeaks.length; i++) {
            // Only evaluate above the fitted average precision
            if (radius[i] < fitSigmaS) {
                continue;
            }
            // Note the RandomModelFunction evaluates g(r)stoch + 1
            final double grProtein = grPeaks[i] - (grStoch[i] - 1);
            if (grProtein > settings.grProteinThreshold) {
                // Failed fit
                ImageJUtils.log("  Failed to fit %s: g(r)protein %s > %s @ r=%s", randomModel.getName(), MathUtils.rounded(grProtein, 4), MathUtils.rounded(settings.grProteinThreshold, 4), MathUtils.rounded(radius[i], 4));
                valid1 = false;
            }
        }
    }
    addResult(randomModel.getName(), resultColour, valid1, fitSigmaS, fitProteinDensity, 0, 0, 0, 0, randomModelAdjustedR2);
    return parameters;
}
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)

Example 67 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 68 with Nullable

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

the class PeakResultsReader method createNStormResult.

/**
 * Creates the NStorm result.
 *
 * @param line the line
 * @param readPhotons Set to {@code true} if there is a Photons field
 * @return the peak result
 */
@Nullable
private static PeakResult createNStormResult(String line, boolean readPhotons) {
    // Ywc
    try (Scanner scanner = new Scanner(line)) {
        scanner.useDelimiter(tabPattern);
        scanner.useLocale(Locale.US);
        // channelName
        scanner.next();
        // x
        scanner.nextFloat();
        // y
        scanner.nextFloat();
        final float xc = scanner.nextFloat();
        final float yc = scanner.nextFloat();
        final float height = scanner.nextFloat();
        final float area = scanner.nextFloat();
        final float width = scanner.nextFloat();
        // phi
        scanner.nextFloat();
        final float ax = scanner.nextFloat();
        final float bg = scanner.nextFloat();
        // intensity
        scanner.nextFloat();
        final int frame = scanner.nextInt();
        final int length = scanner.nextInt();
        double photons = 0;
        if (readPhotons) {
            // These are not needed
            // Link
            scanner.next();
            // Valid
            scanner.next();
            // Z
            scanner.next();
            // Zc
            scanner.next();
            photons = scanner.nextDouble();
        }
        // The coordinates are in nm
        // The values are in ADUs. The area value is the signal.
        // The following relationship holds when length == 1:
        // Area = Height * 2 * pi * (Width / (pixel_pitch*2) )^2
        // => Pixel_pitch = 0.5 * Width / sqrt(Area / (Height * 2 * pi))
        final float[] params = new float[SIZE_TWO_AXIS];
        params[PeakResult.BACKGROUND] = bg;
        params[PeakResult.INTENSITY] = area;
        params[PeakResult.X] = xc;
        params[PeakResult.Y] = yc;
        // Convert width (2*SD) to SD
        final float sd = width / 2f;
        // Convert to separate XY widths using the axial ratio
        if (ax == 1) {
            params[INDEX_SX] = sd;
            params[INDEX_SY] = sd;
        } else {
            // Ensure the axial ratio is long/short
            final double a = Math.sqrt((ax < 1) ? 1.0 / ax : ax);
            params[INDEX_SX] = (float) (sd * a);
            params[INDEX_SY] = (float) (sd / a);
        }
        // Store the photons in the error value
        return new ExtendedPeakResult(frame, (int) xc, (int) yc, height, photons, 0.0f, 0, params, null, frame + length - 1, 0);
    } catch (final NoSuchElementException ignored) {
    // Ignore
    }
    return null;
}
Also used : Scanner(java.util.Scanner) NoSuchElementException(java.util.NoSuchElementException) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 69 with Nullable

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

the class BenchmarkFilterAnalysis method score.

@Nullable
@Override
public SearchResult<FilterScore>[] score(double[][] points) {
    gaIteration++;
    SimpleFilterScore max = filterScoreOptimum;
    final FilterScoreResult[] scoreResults = scoreFilters(setStrength(new FilterSet(searchSpaceToFilters(points))), false);
    if (scoreResults == null) {
        return null;
    }
    @SuppressWarnings("unchecked") final SearchResult<FilterScore>[] scores = new SearchResult[scoreResults.length];
    for (int index = 0; index < scoreResults.length; index++) {
        final FilterScoreResult scoreResult = scoreResults[index];
        final SimpleFilterScore result = new SimpleFilterScore(scoreResult, true, scoreResult.criteria >= minCriteria);
        if (result.compareTo(max) < 0) {
            max = result;
        }
        scores[index] = new SearchResult<>(result.result.filter.getParameters(), result);
    }
    filterScoreOptimum = max;
    // Add the best filter to the table
    // This filter may not have been part of the scored subset so use the entire results set for
    // reporting
    final DirectFilter filter = max.result.filter;
    final FractionClassificationResult r = scoreFilter(filter, DEFAULT_MINIMUM_FILTER, gaResultsList, coordinateStore);
    final StringBuilder text = createResult(filter, r);
    add(text, gaIteration);
    gaWindow.accept(text.toString());
    return scores;
}
Also used : FilterSet(uk.ac.sussex.gdsc.smlm.results.filter.FilterSet) IDirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.IDirectFilter) DirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.DirectFilter) FractionClassificationResult(uk.ac.sussex.gdsc.core.match.FractionClassificationResult) SearchResult(uk.ac.sussex.gdsc.smlm.search.SearchResult) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 70 with Nullable

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

the class BenchmarkFilterAnalysis method readFilterSets.

@Nullable
@SuppressWarnings("unchecked")
private List<FilterSet> readFilterSets() {
    if (extraOptions) {
        final MultiPathFilter multiFilter = BenchmarkSpotFit.getMultiFilter();
        if (multiFilter != null) {
            final IDirectFilter f = multiFilter.getFilter();
            if (f instanceof DirectFilter) {
                final GenericDialog gd = new GenericDialog(TITLE);
                gd.addMessage("Use an identical filter to " + BenchmarkSpotFit.TITLE);
                gd.enableYesNoCancel();
                gd.hideCancelButton();
                gd.showDialog();
                if (gd.wasOKed()) {
                    final List<FilterSet> filterSets = new ArrayList<>(1);
                    final List<Filter> filters = new ArrayList<>(1);
                    filters.add((DirectFilter) f);
                    final FilterSet filterSet = new FilterSet(filters);
                    filterSets.add(filterSet);
                    resetParametersFromFitting();
                    createResultsPrefix2();
                    return filterSets;
                }
            }
        }
    }
    GUIFilterSettings filterSettings = SettingsManager.readGuiFilterSettings(0);
    final String filename = ImageJUtils.getFilename("Filter_File", filterSettings.getFilterSetFilename());
    if (filename != null) {
        IJ.showStatus("Reading filters ...");
        filterSettings = filterSettings.toBuilder().setFilterSetFilename(filename).build();
        // Allow the filters to be cached
        final Triple<String, Long, List<FilterSet>> filterCache = LAST_FILTER_LIST.get();
        if (isSameFile(filename, filterCache)) {
            final GenericDialog gd = new GenericDialog(TITLE);
            gd.hideCancelButton();
            gd.addMessage("The same filter file was selected.");
            gd.addCheckbox("Re-use_filters", settings.reUseFilters);
            gd.showDialog();
            if (!gd.wasCanceled()) {
                settings.reUseFilters = gd.getNextBoolean();
                if (settings.reUseFilters) {
                    SettingsManager.writeSettings(filterSettings);
                    return filterCache.getRight();
                }
            }
        }
        final Path path = Paths.get(filename);
        try (BufferedReader input = new BufferedReader(new UnicodeReader(Files.newInputStream(path), null))) {
            // Use the instance so we can catch the exception
            final Object o = FilterXStreamUtils.getXStreamInstance().fromXML(input);
            if (!(o instanceof List<?>)) {
                IJ.log("No filter sets defined in the specified file: " + filename);
                return null;
            }
            SettingsManager.writeSettings(filterSettings);
            List<FilterSet> filterSets = (List<FilterSet>) o;
            if (containsStandardFilters(filterSets)) {
                IJ.log("Filter sets must contain 'Direct' filters");
                return null;
            }
            // Check they are not empty lists
            final List<FilterSet> filterSets2 = new LinkedList<>();
            for (final FilterSet filterSet : filterSets) {
                if (filterSet.size() != 0) {
                    filterSets2.add(filterSet);
                } else {
                    IJ.log("Filter set empty: " + filterSet.getName());
                }
            }
            if (filterSets2.isEmpty()) {
                IJ.log("All Filter sets are empty");
                return null;
            }
            // Maintain the same list type
            filterSets.clear();
            filterSets.addAll(filterSets2);
            // Option to enumerate filters
            filterSets = expandFilters(filterSets);
            // Save for re-use
            LAST_FILTER_LIST.set(Triple.of(filename, getLastModified(path), filterSets));
            return filterSets;
        } catch (final Exception ex) {
            IJ.log("Unable to load the filter sets from file: " + ex.getMessage());
        } finally {
            IJ.showStatus("");
        }
    }
    return null;
}
Also used : Path(java.nio.file.Path) FilterSet(uk.ac.sussex.gdsc.smlm.results.filter.FilterSet) IDirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.IDirectFilter) DirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.DirectFilter) ArrayList(java.util.ArrayList) IDirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.IDirectFilter) UnicodeReader(uk.ac.sussex.gdsc.core.utils.UnicodeReader) LinkedList(java.util.LinkedList) XStreamException(com.thoughtworks.xstream.XStreamException) IOException(java.io.IOException) ConcurrentRuntimeException(org.apache.commons.lang3.concurrent.ConcurrentRuntimeException) GUIFilterSettings(uk.ac.sussex.gdsc.smlm.ij.settings.GUIProtos.GUIFilterSettings) Filter(uk.ac.sussex.gdsc.smlm.results.filter.Filter) IDirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.IDirectFilter) MultiPathFilter(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter) DirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.DirectFilter) MaximaSpotFilter(uk.ac.sussex.gdsc.smlm.filters.MaximaSpotFilter) GenericDialog(ij.gui.GenericDialog) NonBlockingGenericDialog(ij.gui.NonBlockingGenericDialog) ExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog) MultiPathFilter(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter) BufferedReader(java.io.BufferedReader) ArrayList(java.util.ArrayList) SettingsList(uk.ac.sussex.gdsc.core.utils.SettingsList) List(java.util.List) LinkedList(java.util.LinkedList) LocalList(uk.ac.sussex.gdsc.core.utils.LocalList) 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