Examples with SelectedResultStore uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter.SelectedResultStore used on opensource projects

Example 1 with SelectedResultStore

use of uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter.SelectedResultStore in project GDSC-SMLM by aherbert.

the class FitWorker method run.

/**
 * Locate all the peaks in the image specified by the fit job.
 *
 * <p>WARNING: The FitWorker fits a sub-region of the data for each maxima. It then updates the
 * FitResult parameters with an offset reflecting the position. The initialParameters are not
 * updated with this offset unless configured.
 *
 * @param job The fit job
 */
public void run(FitJob job) {
    final long start = System.nanoTime();
    job.start();
    this.job = job;
    benchmarking = false;
    this.slice = job.slice;
    // Used for debugging
    // if (logger == null) logger = new gdsc.fitting.logging.ConsoleLogger();
    // Crop to the ROI
    cc = new CoordinateConverter(job.bounds);
    // Note if the bounds change for efficient caching.
    newBounds = !cc.dataBounds.equals(lastBounds);
    if (newBounds) {
        lastBounds = cc.dataBounds;
    }
    final int width = cc.dataBounds.width;
    final int height = cc.dataBounds.height;
    borderLimitX = width - border;
    borderLimitY = height - border;
    data = job.data;
    // This is tied to the input data
    dataEstimator = null;
    // relative to the global origin.
    if (isFitCameraCounts) {
        cameraModel.removeBias(cc.dataBounds, data);
    } else {
        cameraModel.removeBiasAndGain(cc.dataBounds, data);
    }
    final FitParameters params = job.getFitParameters();
    this.endT = (params != null) ? params.endT : -1;
    candidates = indentifySpots(job, width, height, params);
    if (candidates.getSize() == 0) {
        finishJob(job, start);
        return;
    }
    fittedBackground = new Statistics();
    // Always get the noise and store it with the results.
    if (params != null && !Float.isNaN(params.noise)) {
        noise = params.noise;
        fitConfig.setNoise(noise);
    } else if (calculateNoise) {
        noise = estimateNoise();
        fitConfig.setNoise(noise);
    }
    // System.out.printf("Slice %d : Noise = %g\n", slice, noise);
    if (logger != null) {
        LoggerUtils.log(logger, Level.INFO, "Slice %d: Noise = %f", slice, noise);
    }
    final ImageExtractor ie = ImageExtractor.wrap(data, width, height);
    double[] region = null;
    final float offsetx = cc.dataBounds.x;
    final float offsety = cc.dataBounds.y;
    if (params != null && params.fitTask == FitTask.MAXIMA_IDENITIFICATION) {
        final float sd0 = (float) xsd;
        final float sd1 = (float) ysd;
        for (int n = 0; n < candidates.getSize(); n++) {
            // Find the background using the perimeter of the data.
            // TODO - Perhaps the Gaussian Fitter should be used to produce the initial estimates but no
            // actual fit done.
            // This would produce coords using the centre-of-mass.
            final Candidate candidate = candidates.get(n);
            int x = candidate.x;
            int y = candidate.y;
            final Rectangle regionBounds = ie.getBoxRegionBounds(x, y, fitting);
            region = ie.crop(regionBounds, region);
            final float b = (float) Gaussian2DFitter.getBackground(region, regionBounds.width, regionBounds.height, 1);
            // Offset the coords to the centre of the pixel. Note the bounds will be added later.
            // Subtract the background to get the amplitude estimate then convert to signal.
            final float amplitude = candidate.intensity - ((relativeIntensity) ? 0 : b);
            final float signal = (float) (amplitude * 2.0 * Math.PI * sd0 * sd1);
            final int index = y * width + x;
            x += offsetx;
            y += offsety;
            final float[] peakParams = new float[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
            peakParams[Gaussian2DFunction.BACKGROUND] = b;
            peakParams[Gaussian2DFunction.SIGNAL] = signal;
            peakParams[Gaussian2DFunction.X_POSITION] = x + 0.5f;
            peakParams[Gaussian2DFunction.Y_POSITION] = y + 0.5f;
            // peakParams[Gaussian2DFunction.Z_POSITION] = 0;
            peakParams[Gaussian2DFunction.X_SD] = sd0;
            peakParams[Gaussian2DFunction.Y_SD] = sd1;
            // peakParams[Gaussian2DFunction.ANGLE] = 0;
            final float u = (float) Gaussian2DPeakResultHelper.getMeanSignalUsingP05(signal, sd0, sd1);
            sliceResults.add(createResult(x, y, data[index], 0, noise, u, peakParams, null, n, 0));
        }
    } else {
        initialiseFitting();
        // Smooth the data to provide initial background estimates
        final float[] smoothedData = backgroundSmoothing.process(data, width, height);
        final ImageExtractor ie2 = ImageExtractor.wrap(smoothedData, width, height);
        // Perform the Gaussian fit
        // The SpotFitter is used to create a dynamic MultiPathFitResult object.
        // This is then passed to a multi-path filter. Thus the same fitting decision process
        // is used when benchmarking and when running on actual data.
        // Note: The SpotFitter labels each PreprocessedFitResult using the offset in the FitResult
        // object.
        // The initial params and deviations can then be extracted for the results that pass the
        // filter.
        MultiPathFilter filter;
        final IMultiPathFitResults multiPathResults = this;
        final SelectedResultStore store = this;
        coordinateStore = coordinateStore.resize(cc.dataBounds.x, cc.dataBounds.y, width, height);
        if (params != null && params.fitTask == FitTask.BENCHMARKING) {
            // Run filtering as normal. However in the event that a candidate is missed or some
            // results are not generated we must generate them. This is done in the complete(int)
            // method if we set the benchmarking flag.
            benchmarking = true;
            // Filter using the benchmark filter
            filter = params.benchmarkFilter;
            if (filter == null) {
                // Create a default filter using the standard FitConfiguration to ensure sensible fits
                // are stored as the current slice results.
                // Note the current fit configuration for benchmarking may have minimal filtering settings
                // so we do not use that object.
                final FitConfiguration tmp = new FitConfiguration();
                final double residualsThreshold = 0.4;
                filter = new MultiPathFilter(tmp, createMinimalFilter(PrecisionMethod.POISSON_CRLB), residualsThreshold);
            }
        } else {
            // Filter using the configuration.
            if (this.filter == null) {
                // This can be cached. Q. Clone the config?
                this.filter = new MultiPathFilter(fitConfig, createMinimalFilter(fitConfig.getPrecisionMethod()), config.getResidualsThreshold());
            }
            filter = this.filter;
        }
        // If we are benchmarking then do not generate results dynamically since we will store all
        // results in the fit job.
        dynamicMultiPathFitResult = new DynamicMultiPathFitResult(ie, ie2, !benchmarking);
        // dynamicMultiPathFitResult = new DynamicMultiPathFitResult(ie, false);
        // The local background computation is only required for the precision method.
        // Also compute it when benchmarking.
        localBackground = benchmarking || fitConfig.getPrecisionMethodValue() == PrecisionMethod.MORTENSEN_LOCAL_BACKGROUND_VALUE;
        // Debug where the fit config may be different between benchmarking and fitting
        if (slice == -1) {
            fitConfig.initialise(1, 1, 1);
            final String newLine = System.lineSeparator();
            final String tmpdir = System.getProperty("java.io.tmpdir");
            try (BufferedWriter writer = Files.newBufferedWriter(Paths.get(tmpdir, String.format("config.%d.txt", slice)))) {
                JsonFormat.printer().appendTo(config.getFitEngineSettings(), writer);
            } catch (final IOException ex) {
                logger.log(Level.SEVERE, "Unable to write message", ex);
            }
            FileUtils.save(Paths.get(tmpdir, String.format("filter.%d.xml", slice)).toString(), filter.toXml());
            // filter.setDebugFile(String.format("/tmp/fitWorker.%b.txt", benchmarking));
            final StringBuilder sb = new StringBuilder();
            sb.append((benchmarking) ? ((uk.ac.sussex.gdsc.smlm.results.filter.Filter) filter.getFilter()).toXml() : fitConfig.getSmartFilterString()).append(newLine);
            sb.append(((uk.ac.sussex.gdsc.smlm.results.filter.Filter) filter.getMinimalFilter()).toXml()).append(newLine);
            sb.append(filter.residualsThreshold).append(newLine);
            sb.append(config.getFailuresLimit()).append(newLine);
            sb.append(config.getDuplicateDistance()).append(":");
            sb.append(config.getDuplicateDistanceAbsolute()).append(newLine);
            if (spotFilter != null) {
                sb.append(spotFilter.getDescription()).append(newLine);
            }
            sb.append("MaxCandidate = ").append(candidates.getSize()).append(newLine);
            for (int i = 0, len = candidates.getLength(); i < len; i++) {
                TextUtils.formatTo(sb, "Fit %d [%d,%d = %.1f]%n", i, candidates.get(i).x, candidates.get(i).y, candidates.get(i).intensity);
            }
            FileUtils.save(Paths.get(tmpdir, String.format("candidates.%d.xml", slice)).toString(), sb.toString());
        }
        FailCounter failCounter = config.getFailCounter();
        if (!benchmarking && params != null && params.pass != null) {
            // We want to store the pass/fail for consecutive candidates
            params.pass = new boolean[candidates.getLength()];
            failCounter = new RecordingFailCounter(params.pass, failCounter);
            filter.select(multiPathResults, failCounter, true, store, coordinateStore);
        } else {
            filter.select(multiPathResults, failCounter, true, store, coordinateStore);
        }
        // Note: We go deeper into the candidate list than max candidate
        // for any candidate where we have a good fit result as an estimate.
        // Q. Should this only be for benchmarking?
        // if (benchmarking)
        // System.out.printf("Slice %d: %d + %d\n", slice, dynamicMultiPathFitResult.extra,
        // candidates.getSize());
        // Create the slice results
        final CandidateList fitted = gridManager.getFittedCandidates();
        sliceResults.ensureCapacity(fitted.getSize());
        for (int i = 0; i < fitted.getSize(); i++) {
            if (fitted.get(i).fit) {
                sliceResults.push(createResult(offsetx, offsety, fitted.get(i)));
            }
        }
        if (logger != null) {
            LoggerUtils.log(logger, Level.INFO, "Slice %d: %d / %d = %s", slice, success, candidates.getSize(), TextUtils.pleural(fitted.getSize(), "result"));
        }
    }
    this.results.addAll(sliceResults);
    finishJob(job, start);
}