Search in sources :

Example 1 with FitResult

use of gdsc.smlm.fitting.FitResult in project GDSC-SMLM by aherbert.

the class GaussianFit method fit.

/**
	 * Fits a single 2D Gaussian to the data. The fit is initialised at the highest
	 * value and then optimised.
	 * <p>
	 * Data must be arranged in yx block order, i.e. height rows of width.
	 * <p>
	 * The angle parameter is only returned if using elliptical Gaussian fitting.
	 * <p>
	 * Note: The fit coordinates should be offset by 0.5 if the input data represents pixels
	 * 
	 * @return Array containing the fitted curve data: Background, Amplitude, PosX, PosY, StdDevX, StdDevY, Angle. Null
	 *         if no fit is possible.
	 */
public double[] fit(float[] data, int width, int height) {
    if (data == null || data.length != width * height)
        return null;
    // Get the limits
    float max = Float.MIN_VALUE;
    int maxIndex = -1;
    for (int i = data.length; i-- > 0; ) {
        float f = data[i];
        if (max < f) {
            max = f;
            maxIndex = i;
        }
    }
    if (maxIndex < 0) {
        return null;
    }
    Gaussian2DFitter gf = createGaussianFitter(false);
    FitResult fitResult = gf.fit(Utils.toDouble(data), width, height, new int[] { maxIndex });
    if (fitResult.getStatus() == FitStatus.OK) {
        chiSquared = fitResult.getError();
        double[] params = fitResult.getParameters();
        // Check bounds
        if (params[3] < 0 || params[3] >= width || params[4] < 0 || params[4] >= height)
            return null;
        // Re-arrange order for backwards compatibility with old code.
        return new double[] { params[0], params[1], params[3], params[4], Gaussian2DFitter.fwhm2sd(params[5]), Gaussian2DFitter.fwhm2sd(params[6]), params[2] };
    }
    return null;
}
Also used : Gaussian2DFitter(gdsc.smlm.fitting.Gaussian2DFitter) FitResult(gdsc.smlm.fitting.FitResult)

Example 2 with FitResult

use of gdsc.smlm.fitting.FitResult in project GDSC-SMLM by aherbert.

the class SpotAnalysis method updateCurrentSlice.

private void updateCurrentSlice(int slice) {
    if (slice != currentSlice) {
        currentSlice = slice;
        double signal = getSignal(slice);
        double noise = smoothSd[slice - 1];
        currentLabel.setText(String.format("Frame %d: Signal = %s, SNR = %s", slice, Utils.rounded(signal, 4), Utils.rounded(signal / noise, 3)));
        drawProfiles();
        // Fit the PSF using a Gaussian
        float[] data2 = (float[]) rawImp.getImageStack().getProcessor(slice).getPixels();
        double[] data = Utils.toDouble(data2);
        FitConfiguration fitConfiguration = new FitConfiguration();
        fitConfiguration.setFitFunction(FitFunction.FIXED);
        fitConfiguration.setBackgroundFitting(true);
        fitConfiguration.setSignalStrength(0);
        fitConfiguration.setCoordinateShift(rawImp.getWidth() / 4.0f);
        fitConfiguration.setComputeResiduals(false);
        fitConfiguration.setComputeDeviations(false);
        Gaussian2DFitter gf = new Gaussian2DFitter(fitConfiguration);
        double[] params = new double[7];
        double psfWidth = Double.parseDouble(widthTextField.getText());
        params[Gaussian2DFunction.BACKGROUND] = smoothMean[slice - 1];
        params[Gaussian2DFunction.SIGNAL] = (gain * signal);
        params[Gaussian2DFunction.X_POSITION] = rawImp.getWidth() / 2.0f;
        params[Gaussian2DFunction.Y_POSITION] = rawImp.getHeight() / 2.0f;
        params[Gaussian2DFunction.X_SD] = params[Gaussian2DFunction.Y_SD] = psfWidth;
        FitResult fitResult = gf.fit(data, rawImp.getWidth(), rawImp.getHeight(), 1, params, new boolean[1]);
        if (fitResult.getStatus() == FitStatus.OK) {
            params = fitResult.getParameters();
            final double spotSignal = params[Gaussian2DFunction.SIGNAL] / gain;
            rawFittedLabel.setText(String.format("Raw fit: Signal = %s, SNR = %s", Utils.rounded(spotSignal, 4), Utils.rounded(spotSignal / noise, 3)));
            ImageROIPainter.addRoi(rawImp, slice, new PointRoi(params[Gaussian2DFunction.X_POSITION], params[Gaussian2DFunction.Y_POSITION]));
        } else {
            rawFittedLabel.setText("");
            rawImp.setOverlay(null);
        }
        // Fit the PSF using a Gaussian
        if (blurImp == null)
            return;
        data2 = (float[]) blurImp.getImageStack().getProcessor(slice).getPixels();
        data = Utils.toDouble(data2);
        params = new double[7];
        //float psfWidth = Float.parseFloat(widthTextField.getText());
        params[Gaussian2DFunction.BACKGROUND] = (float) smoothMean[slice - 1];
        params[Gaussian2DFunction.SIGNAL] = (float) (gain * signal);
        params[Gaussian2DFunction.X_POSITION] = rawImp.getWidth() / 2.0f;
        params[Gaussian2DFunction.Y_POSITION] = rawImp.getHeight() / 2.0f;
        params[Gaussian2DFunction.X_SD] = params[Gaussian2DFunction.Y_SD] = psfWidth;
        fitResult = gf.fit(data, rawImp.getWidth(), rawImp.getHeight(), 1, params, new boolean[1]);
        if (fitResult.getStatus() == FitStatus.OK) {
            params = fitResult.getParameters();
            final double spotSignal = params[Gaussian2DFunction.SIGNAL] / gain;
            blurFittedLabel.setText(String.format("Blur fit: Signal = %s, SNR = %s", Utils.rounded(spotSignal, 4), Utils.rounded(spotSignal / noise, 3)));
            ImageROIPainter.addRoi(blurImp, slice, new PointRoi(params[Gaussian2DFunction.X_POSITION], params[Gaussian2DFunction.Y_POSITION]));
        } else {
            blurFittedLabel.setText("");
            blurImp.setOverlay(null);
        }
    }
}
Also used : FitConfiguration(gdsc.smlm.fitting.FitConfiguration) Gaussian2DFitter(gdsc.smlm.fitting.Gaussian2DFitter) FitResult(gdsc.smlm.fitting.FitResult) PointRoi(ij.gui.PointRoi)

Example 3 with FitResult

use of gdsc.smlm.fitting.FitResult in project GDSC-SMLM by aherbert.

the class TraceMolecules method fitTraces.

private void fitTraces(MemoryPeakResults results, Trace[] traces) {
    // Check if the original image is open and the fit configuration can be extracted
    ImageSource source = results.getSource();
    if (source == null)
        return;
    if (!source.open())
        return;
    FitEngineConfiguration config = (FitEngineConfiguration) XmlUtils.fromXML(results.getConfiguration());
    if (config == null)
        return;
    // Show a dialog asking if the traces should be refit
    ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
    gd.addMessage("Do you want to fit the traces as a single peak using a combined image?");
    gd.addCheckbox("Fit_closest_to_centroid", !fitOnlyCentroid);
    gd.addSlider("Distance_threshold", 0.01, 3, distanceThreshold);
    gd.addSlider("Expansion_factor", 1, 4.5, expansionFactor);
    // Allow fitting settings to be adjusted
    FitConfiguration fitConfig = config.getFitConfiguration();
    gd.addMessage("--- Gaussian fitting ---");
    String[] filterTypes = SettingsManager.getNames((Object[]) DataFilterType.values());
    gd.addChoice("Spot_filter_type", filterTypes, filterTypes[config.getDataFilterType().ordinal()]);
    String[] filterNames = SettingsManager.getNames((Object[]) DataFilter.values());
    gd.addChoice("Spot_filter", filterNames, filterNames[config.getDataFilter(0).ordinal()]);
    gd.addSlider("Smoothing", 0, 2.5, config.getSmooth(0));
    gd.addSlider("Search_width", 0.5, 2.5, config.getSearch());
    gd.addSlider("Border", 0.5, 2.5, config.getBorder());
    gd.addSlider("Fitting_width", 2, 4.5, config.getFitting());
    String[] solverNames = SettingsManager.getNames((Object[]) FitSolver.values());
    gd.addChoice("Fit_solver", solverNames, solverNames[fitConfig.getFitSolver().ordinal()]);
    String[] functionNames = SettingsManager.getNames((Object[]) FitFunction.values());
    gd.addChoice("Fit_function", functionNames, functionNames[fitConfig.getFitFunction().ordinal()]);
    String[] criteriaNames = SettingsManager.getNames((Object[]) FitCriteria.values());
    gd.addChoice("Fit_criteria", criteriaNames, criteriaNames[fitConfig.getFitCriteria().ordinal()]);
    gd.addNumericField("Significant_digits", fitConfig.getSignificantDigits(), 0);
    gd.addNumericField("Coord_delta", fitConfig.getDelta(), 4);
    gd.addNumericField("Lambda", fitConfig.getLambda(), 4);
    gd.addNumericField("Max_iterations", fitConfig.getMaxIterations(), 0);
    gd.addNumericField("Fail_limit", config.getFailuresLimit(), 0);
    gd.addCheckbox("Include_neighbours", config.isIncludeNeighbours());
    gd.addSlider("Neighbour_height", 0.01, 1, config.getNeighbourHeightThreshold());
    gd.addSlider("Residuals_threshold", 0.01, 1, config.getResidualsThreshold());
    //gd.addSlider("Duplicate_distance", 0, 1.5, fitConfig.getDuplicateDistance());
    gd.addMessage("--- Peak filtering ---\nDiscard fits that shift; are too low; or expand/contract");
    gd.addCheckbox("Smart_filter", fitConfig.isSmartFilter());
    gd.addCheckbox("Disable_simple_filter", fitConfig.isDisableSimpleFilter());
    gd.addSlider("Shift_factor", 0.01, 2, fitConfig.getCoordinateShiftFactor());
    gd.addNumericField("Signal_strength", fitConfig.getSignalStrength(), 2);
    gd.addNumericField("Min_photons", fitConfig.getMinPhotons(), 0);
    gd.addSlider("Min_width_factor", 0, 0.99, fitConfig.getMinWidthFactor());
    gd.addSlider("Width_factor", 1.01, 5, fitConfig.getWidthFactor());
    gd.addNumericField("Precision", fitConfig.getPrecisionThreshold(), 2);
    gd.addCheckbox("Debug_failures", debugFailures);
    gd.showDialog();
    if (!gd.wasOKed()) {
        source.close();
        return;
    }
    // Get parameters for the fit
    fitOnlyCentroid = !gd.getNextBoolean();
    distanceThreshold = (float) gd.getNextNumber();
    expansionFactor = (float) gd.getNextNumber();
    config.setDataFilterType(gd.getNextChoiceIndex());
    config.setDataFilter(gd.getNextChoiceIndex(), Math.abs(gd.getNextNumber()), 0);
    config.setSearch(gd.getNextNumber());
    config.setBorder(gd.getNextNumber());
    config.setFitting(gd.getNextNumber());
    fitConfig.setFitSolver(gd.getNextChoiceIndex());
    fitConfig.setFitFunction(gd.getNextChoiceIndex());
    fitConfig.setFitCriteria(gd.getNextChoiceIndex());
    fitConfig.setSignificantDigits((int) gd.getNextNumber());
    fitConfig.setDelta(gd.getNextNumber());
    fitConfig.setLambda(gd.getNextNumber());
    fitConfig.setMaxIterations((int) gd.getNextNumber());
    config.setFailuresLimit((int) gd.getNextNumber());
    config.setIncludeNeighbours(gd.getNextBoolean());
    config.setNeighbourHeightThreshold(gd.getNextNumber());
    config.setResidualsThreshold(gd.getNextNumber());
    fitConfig.setSmartFilter(gd.getNextBoolean());
    fitConfig.setDisableSimpleFilter(gd.getNextBoolean());
    fitConfig.setCoordinateShiftFactor(gd.getNextNumber());
    fitConfig.setSignalStrength(gd.getNextNumber());
    fitConfig.setMinPhotons(gd.getNextNumber());
    fitConfig.setMinWidthFactor(gd.getNextNumber());
    fitConfig.setWidthFactor(gd.getNextNumber());
    fitConfig.setPrecisionThreshold(gd.getNextNumber());
    // Check arguments
    try {
        Parameters.isAboveZero("Distance threshold", distanceThreshold);
        Parameters.isAbove("Expansion factor", expansionFactor, 1);
        Parameters.isAboveZero("Search_width", config.getSearch());
        Parameters.isAboveZero("Fitting_width", config.getFitting());
        Parameters.isAboveZero("Significant digits", fitConfig.getSignificantDigits());
        Parameters.isAboveZero("Delta", fitConfig.getDelta());
        Parameters.isAboveZero("Lambda", fitConfig.getLambda());
        Parameters.isAboveZero("Max iterations", fitConfig.getMaxIterations());
        Parameters.isPositive("Failures limit", config.getFailuresLimit());
        Parameters.isPositive("Neighbour height threshold", config.getNeighbourHeightThreshold());
        Parameters.isPositive("Residuals threshold", config.getResidualsThreshold());
        Parameters.isPositive("Coordinate Shift factor", fitConfig.getCoordinateShiftFactor());
        Parameters.isPositive("Signal strength", fitConfig.getSignalStrength());
        Parameters.isPositive("Min photons", fitConfig.getMinPhotons());
        Parameters.isPositive("Min width factor", fitConfig.getMinWidthFactor());
        Parameters.isPositive("Width factor", fitConfig.getWidthFactor());
        Parameters.isPositive("Precision threshold", fitConfig.getPrecisionThreshold());
    } catch (IllegalArgumentException e) {
        IJ.error(TITLE, e.getMessage());
        source.close();
        return;
    }
    debugFailures = gd.getNextBoolean();
    if (!PeakFit.configureSmartFilter(globalSettings, filename))
        return;
    if (!PeakFit.configureDataFilter(globalSettings, filename, false))
        return;
    if (!PeakFit.configureFitSolver(globalSettings, filename, false))
        return;
    // Adjust settings for a single maxima
    config.setIncludeNeighbours(false);
    fitConfig.setDuplicateDistance(0);
    // Create a fit engine
    MemoryPeakResults refitResults = new MemoryPeakResults();
    refitResults.copySettings(results);
    refitResults.setName(results.getName() + " Trace Fit");
    refitResults.setSortAfterEnd(true);
    refitResults.begin();
    // No border since we know where the peaks are and we must not miss them due to truncated searching 
    FitEngine engine = new FitEngine(config, refitResults, Prefs.getThreads(), FitQueue.BLOCKING);
    // Either : Only fit the centroid
    // or     : Extract a bigger region, allowing all fits to run as normal and then 
    //          find the correct spot using Euclidian distance.
    // Set up the limits
    final double stdDev = FastMath.max(fitConfig.getInitialPeakStdDev0(), fitConfig.getInitialPeakStdDev1());
    float fitWidth = (float) (stdDev * config.getFitting() * ((fitOnlyCentroid) ? 1 : expansionFactor));
    IJ.showStatus("Refitting traces ...");
    List<JobItem> jobItems = new ArrayList<JobItem>(traces.length);
    int singles = 0;
    int fitted = 0;
    for (int n = 0; n < traces.length; n++) {
        Trace trace = traces[n];
        if (n % 32 == 0)
            IJ.showProgress(n, traces.length);
        // Skip traces with one peak
        if (trace.size() == 1) {
            singles++;
            // Use the synchronized method to avoid thread clashes with the FitEngine
            refitResults.addSync(trace.getHead());
            continue;
        }
        Rectangle bounds = new Rectangle();
        double[] combinedNoise = new double[1];
        float[] data = buildCombinedImage(source, trace, fitWidth, bounds, combinedNoise, false);
        if (data == null)
            continue;
        // Fit the combined image
        FitParameters params = new FitParameters();
        params.noise = (float) combinedNoise[0];
        float[] centre = trace.getCentroid();
        if (fitOnlyCentroid) {
            int newX = (int) Math.round(centre[0]) - bounds.x;
            int newY = (int) Math.round(centre[1]) - bounds.y;
            params.maxIndices = new int[] { newY * bounds.width + newX };
        } else {
            params.filter = new ArrayList<float[]>();
            params.filter.add(new float[] { centre[0] - bounds.x, centre[1] - bounds.y });
            params.distanceThreshold = distanceThreshold;
        }
        // This is not needed since the bounds are passed using the FitJob
        //params.setOffset(new float[] { bounds.x, bounds.y });
        int startT = trace.getHead().getFrame();
        params.endT = trace.getTail().getFrame();
        ParameterisedFitJob job = new ParameterisedFitJob(n, params, startT, data, bounds);
        jobItems.add(new JobItem(job, trace, centre));
        engine.run(job);
        fitted++;
    }
    engine.end(false);
    IJ.showStatus("");
    IJ.showProgress(1);
    // Check the success ...
    FitStatus[] values = FitStatus.values();
    int[] statusCount = new int[values.length + 1];
    ArrayList<String> names = new ArrayList<String>(Arrays.asList(SettingsManager.getNames((Object[]) values)));
    names.add(String.format("No maxima within %.2f of centroid", distanceThreshold));
    int separated = 0;
    int success = 0;
    final int debugLimit = 3;
    for (JobItem jobItem : jobItems) {
        int id = jobItem.getId();
        ParameterisedFitJob job = jobItem.job;
        Trace trace = jobItem.trace;
        int[] indices = job.getIndices();
        FitResult fitResult = null;
        int status;
        if (indices.length < 1) {
            status = values.length;
        } else if (indices.length > 1) {
            // Choose the first OK result. This is all that matters for the success reporting
            for (int n = 0; n < indices.length; n++) {
                if (job.getFitResult(n).getStatus() == FitStatus.OK) {
                    fitResult = job.getFitResult(n);
                    break;
                }
            }
            // Otherwise use the closest failure. 
            if (fitResult == null) {
                final float[] centre = traces[id].getCentroid();
                double minD = Double.POSITIVE_INFINITY;
                for (int n = 0; n < indices.length; n++) {
                    // Since the fit has failed we use the initial parameters
                    final double[] params = job.getFitResult(n).getInitialParameters();
                    final double dx = params[Gaussian2DFunction.X_POSITION] - centre[0];
                    final double dy = params[Gaussian2DFunction.Y_POSITION] - centre[1];
                    final double d = dx * dx + dy * dy;
                    if (minD > d) {
                        minD = d;
                        fitResult = job.getFitResult(n);
                    }
                }
            }
            status = fitResult.getStatus().ordinal();
        } else {
            fitResult = job.getFitResult(0);
            status = fitResult.getStatus().ordinal();
        }
        // All jobs have only one peak
        statusCount[status]++;
        // Debug why any fits failed
        if (fitResult == null || fitResult.getStatus() != FitStatus.OK) {
            refitResults.addAll(trace.getPoints());
            separated += trace.size();
            if (debugFailures) {
                FitStatus s = (fitResult == null) ? FitStatus.UNKNOWN : fitResult.getStatus();
                // Only display the first n per category to limit the number of images
                double[] noise = new double[1];
                if (statusCount[status] <= debugLimit) {
                    Rectangle bounds = new Rectangle();
                    buildCombinedImage(source, trace, fitWidth, bounds, noise, true);
                    float[] centre = trace.getCentroid();
                    Utils.display(String.format("Trace %d (n=%d) : x=%f,y=%f", id, trace.size(), centre[0], centre[1]), slices);
                    switch(s) {
                        case INSUFFICIENT_PRECISION:
                            float precision = (Float) fitResult.getStatusData();
                            IJ.log(String.format("Trace %d (n=%d) : %s = %f", id, trace.size(), names.get(status), precision));
                            break;
                        case INSUFFICIENT_SIGNAL:
                            if (noise[0] == 0)
                                noise[0] = getCombinedNoise(trace);
                            float snr = (Float) fitResult.getStatusData();
                            IJ.log(String.format("Trace %d (n=%d) : %s = %f (noise=%.2f)", id, trace.size(), names.get(status), snr, noise[0]));
                            break;
                        case COORDINATES_MOVED:
                        case OUTSIDE_FIT_REGION:
                        case WIDTH_DIVERGED:
                            float[] shift = (float[]) fitResult.getStatusData();
                            IJ.log(String.format("Trace %d (n=%d) : %s = %.3f,%.3f", id, trace.size(), names.get(status), shift[0], shift[1]));
                            break;
                        default:
                            IJ.log(String.format("Trace %d (n=%d) : %s", id, trace.size(), names.get(status)));
                            break;
                    }
                }
            }
        } else {
            success++;
            if (debugFailures) {
                // Only display the first n per category to limit the number of images
                double[] noise = new double[1];
                if (statusCount[status] <= debugLimit) {
                    Rectangle bounds = new Rectangle();
                    buildCombinedImage(source, trace, fitWidth, bounds, noise, true);
                    float[] centre = trace.getCentroid();
                    Utils.display(String.format("Trace %d (n=%d) : x=%f,y=%f", id, trace.size(), centre[0], centre[1]), slices);
                }
            }
        }
    }
    IJ.log(String.format("Trace fitting : %d singles : %d / %d fitted : %d separated", singles, success, fitted, separated));
    if (separated > 0) {
        IJ.log("Reasons for fit failure :");
        // Start at i=1 to skip FitStatus.OK
        for (int i = 1; i < statusCount.length; i++) {
            if (statusCount[i] != 0)
                IJ.log("  " + names.get(i) + " = " + statusCount[i]);
        }
    }
    refitResults.end();
    MemoryPeakResults.addResults(refitResults);
    source.close();
}
Also used : FitParameters(gdsc.smlm.engine.FitParameters) ArrayList(java.util.ArrayList) Rectangle(java.awt.Rectangle) FitStatus(gdsc.smlm.fitting.FitStatus) MemoryPeakResults(gdsc.smlm.results.MemoryPeakResults) ParameterisedFitJob(gdsc.smlm.engine.ParameterisedFitJob) FitEngineConfiguration(gdsc.smlm.engine.FitEngineConfiguration) ExtendedGenericDialog(ij.gui.ExtendedGenericDialog) ClusterPoint(gdsc.core.clustering.ClusterPoint) Trace(gdsc.smlm.results.Trace) FitEngine(gdsc.smlm.engine.FitEngine) FitConfiguration(gdsc.smlm.fitting.FitConfiguration) FitResult(gdsc.smlm.fitting.FitResult) ImageSource(gdsc.smlm.results.ImageSource)

Example 4 with FitResult

use of gdsc.smlm.fitting.FitResult 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 data
	 * @param width
	 * @param 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.
	 *         <p>
	 *         Null if no fit is possible.
	 */
private double[] fitSingle(Gaussian2DFitter gf, float[] data, int width, int height, int index, double estimatedHeight) {
    this.fitResult = gf.fit(Utils.toDouble(data), width, height, new int[] { index }, new double[] { estimatedHeight });
    if (fitResult.getStatus() == FitStatus.OK) {
        chiSquared = fitResult.getError();
        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.getParameterStdDev(), fitResult.getNumberOfPeaks(), fitResult.getNumberOfFittedParameters(), fitResult.getStatusData(), fitResult.getIterations(), fitResult.getEvaluations());
            return null;
        }
        return params;
    }
    return null;
}
Also used : FitResult(gdsc.smlm.fitting.FitResult)

Example 5 with FitResult

use of gdsc.smlm.fitting.FitResult in project GDSC-SMLM by aherbert.

the class FitWorker method add.

/*
	 * (non-Javadoc)
	 * 
	 * @see gdsc.smlm.results.filter.MultiPathFilter.SelectedResultStore#add(gdsc.smlm.results.filter.MultiPathFilter.
	 * SelectedResult)
	 */
public void add(SelectedResult selectedResult) {
    // TODO - Print the current state of the dynamicMultiPathFitResult to file.
    // This will allow debugging what is different between the benchmark fit and the PeakFit.
    // Output:
    // slice
    // candidate Id
    // Initial and final params for each fit result.
    // Details of the selected result.
    // Then try to figure out why the benchmark fit deviates from PeakFit.
    // Add to the slice results.
    final PreprocessedPeakResult[] results = selectedResult.results;
    if (results == null)
        return;
    final int currrentSize = sliceResults.size();
    final int candidateId = dynamicMultiPathFitResult.candidateId;
    final FitResult fitResult = (FitResult) selectedResult.fitResult.data;
    // The background for each result was the local background. We want the fitted global background
    final float background = (float) fitResult.getParameters()[0];
    final double[] dev = fitResult.getParameterStdDev();
    if (queueSize != 0)
        throw new RuntimeException("There are results queued already!");
    for (int i = 0; i < results.length; i++) {
        if (results[i].isExistingResult())
            continue;
        if (results[i].isNewResult()) {
            final double[] p = results[i].toGaussian2DParameters();
            // Store slice results relative to the data frame (not the global bounds)
            // Convert back so that 0,0 is the top left of the data bounds
            p[Gaussian2DFunction.X_POSITION] -= cc.dataBounds.x;
            p[Gaussian2DFunction.Y_POSITION] -= cc.dataBounds.y;
            final float[] params = new float[7];
            params[Gaussian2DFunction.BACKGROUND] = background;
            for (int j = 1; j < 7; j++) params[j] = (float) p[j];
            final float[] paramsDev;
            if (dev == null) {
                paramsDev = null;
            } else {
                paramsDev = new float[7];
                paramsDev[Gaussian2DFunction.BACKGROUND] = (float) dev[Gaussian2DFunction.BACKGROUND];
                final int offset = results[i].getId() * 6;
                for (int j = 1; j < 7; j++) paramsDev[j] = (float) dev[offset + j];
            }
            addSingleResult(results[i].getCandidateId(), params, paramsDev, fitResult.getError(), results[i].getNoise());
            if (logger != null) {
                // Show the shift, signal and width spread
                PreprocessedPeakResult peak = results[i];
                logger.info("Fit OK %d (%.1f,%.1f) [%d]: Shift = %.3f,%.3f : SNR = %.2f : Width = %.2f,%.2f", peak.getCandidateId(), peak.getX(), peak.getY(), peak.getId(), Math.sqrt(peak.getXRelativeShift2()), Math.sqrt(peak.getYRelativeShift2()), peak.getSNR(), peak.getXSDFactor(), peak.getYSDFactor());
            }
        } else {
        // This is a candidate that passed validation. Store the estimate as passing the primary filter.
        // We now do this is the pass() method.
        //storeEstimate(results[i].getCandidateId(), results[i], FILTER_RANK_PRIMARY);
        }
    }
    job.setFitResult(candidateId, fitResult);
    // Reporting
    if (this.counter != null) {
        FitType fitType = dynamicMultiPathFitResult.fitType;
        if (selectedResult.fitResult.getStatus() == 0) {
            fitType.setOK(true);
            if (dynamicMultiPathFitResult.getSuperMultiFitResult() == selectedResult.fitResult)
                fitType.setMultiOK(true);
            else if (dynamicMultiPathFitResult.getSuperMultiDoubletFitResult() == selectedResult.fitResult)
                fitType.setMultiDoubletOK(true);
            else if (dynamicMultiPathFitResult.getSuperDoubletFitResult() == selectedResult.fitResult)
                fitType.setDoubletOK(true);
        }
        add(fitType);
    }
    if (logger != null) {
        switch(fitResult.getStatus()) {
            case OK:
                // We log good results in the loop above. 
                break;
            case BAD_PARAMETERS:
            case FAILED_TO_ESTIMATE_WIDTH:
                logger.info("Bad parameters: %s", Arrays.toString(fitResult.getInitialParameters()));
                break;
            default:
                logger.info(fitResult.getStatus().toString());
                break;
        }
    }
    // Debugging
    if (logger2 != null) {
        double[] peakParams = fitResult.getParameters();
        if (peakParams != null) {
            // Parameters are the raw values from fitting the region. Convert for logging.
            peakParams = Arrays.copyOf(peakParams, peakParams.length);
            int npeaks = peakParams.length / 6;
            for (int i = 0; i < npeaks; i++) {
                peakParams[i * 6 + Gaussian2DFunction.X_POSITION] += cc.fromFitRegionToGlobalX();
                peakParams[i * 6 + Gaussian2DFunction.Y_POSITION] += cc.fromFitRegionToGlobalY();
                peakParams[i * 6 + Gaussian2DFunction.SHAPE] *= 180.0 / Math.PI;
            }
        }
        final int x = candidates.get(candidateId).x;
        final int y = candidates.get(candidateId).y;
        logger2.debug("%d:%d [%d,%d] %s (%s) = %s\n", slice, candidateId, cc.fromDataToGlobalX(x), cc.fromDataToGlobalY(y), fitResult.getStatus(), fitResult.getStatusData(), Arrays.toString(peakParams));
    }
    // Check if there were any new results
    int npeaks = sliceResults.size() - currrentSize;
    if (npeaks != 0) {
        success++;
        // Support for post-processing filter 
        if (resultFilter != null) {
            // Check all result peaks for the distance to the filter positions
            PeakResult[] peakResults = new PeakResult[npeaks];
            for (int i = sliceResults.size(); npeaks-- > 0; ) {
                peakResults[npeaks] = sliceResults.get(--i);
            }
            resultFilter.filter(fitResult, candidateId, peakResults);
        }
    }
}
Also used : PreprocessedPeakResult(gdsc.smlm.results.filter.PreprocessedPeakResult) FitResult(gdsc.smlm.fitting.FitResult) MultiPathFitResult(gdsc.smlm.results.filter.MultiPathFitResult) PreprocessedPeakResult(gdsc.smlm.results.filter.PreprocessedPeakResult) PeakResult(gdsc.smlm.results.PeakResult) IdPeakResult(gdsc.smlm.results.IdPeakResult) ExtendedPeakResult(gdsc.smlm.results.ExtendedPeakResult)

Aggregations

FitResult (gdsc.smlm.fitting.FitResult)6 FitConfiguration (gdsc.smlm.fitting.FitConfiguration)2 Gaussian2DFitter (gdsc.smlm.fitting.Gaussian2DFitter)2 MultiPathFitResult (gdsc.smlm.results.filter.MultiPathFitResult)2 ClusterPoint (gdsc.core.clustering.ClusterPoint)1 BasePoint (gdsc.core.match.BasePoint)1 FitEngine (gdsc.smlm.engine.FitEngine)1 FitEngineConfiguration (gdsc.smlm.engine.FitEngineConfiguration)1 FitParameters (gdsc.smlm.engine.FitParameters)1 ParameterisedFitJob (gdsc.smlm.engine.ParameterisedFitJob)1 FitStatus (gdsc.smlm.fitting.FitStatus)1 PeakResultPoint (gdsc.smlm.ij.plugins.ResultsMatchCalculator.PeakResultPoint)1 ExtendedPeakResult (gdsc.smlm.results.ExtendedPeakResult)1 IdPeakResult (gdsc.smlm.results.IdPeakResult)1 ImageSource (gdsc.smlm.results.ImageSource)1 MemoryPeakResults (gdsc.smlm.results.MemoryPeakResults)1 PeakResult (gdsc.smlm.results.PeakResult)1 Trace (gdsc.smlm.results.Trace)1 BasePreprocessedPeakResult (gdsc.smlm.results.filter.BasePreprocessedPeakResult)1 PreprocessedPeakResult (gdsc.smlm.results.filter.PreprocessedPeakResult)1