Search in sources :

Example 1 with FailCounter

use of uk.ac.sussex.gdsc.smlm.results.count.FailCounter in project GDSC-SMLM by aherbert.

the class FailCountManager method analyseData.

private void analyseData() {
    final LocalList<FailCountData> failCountData = failCountDataRef.get();
    if (failCountData.isEmpty()) {
        IJ.error(TITLE, "No fail count data in memory");
        return;
    }
    if (!showAnalysisDialog()) {
        return;
    }
    final int maxCons = getMaxConsecutiveFailCount(failCountData);
    final int maxFail = getMaxFailCount(failCountData);
    // Create a set of fail counters
    final LocalList<FailCounter> counters = new LocalList<>();
    final TByteArrayList type = new TByteArrayList();
    for (int i = 0; i <= maxCons; i++) {
        counters.add(ConsecutiveFailCounter.create(i));
    }
    fill(type, counters, 0);
    // The other counters are user configured.
    // Ideally this would be a search to optimise the best parameters
    // for each counter as any enumeration may be way off the mark.
    // Note that 0 failures in a window can be scored using the consecutive fail counter.
    int max = Math.min(maxFail, settings.getRollingCounterMaxAllowedFailures());
    for (int fail = MathUtils.clip(1, maxFail, settings.getRollingCounterMinAllowedFailures()); fail <= max; fail++) {
        // Note that n-1 failures in window n can be scored using the consecutive fail counter.
        for (int window = Math.max(fail + 2, settings.getRollingCounterMinWindow()); window <= settings.getRollingCounterMaxWindow(); window++) {
            counters.add(RollingWindowFailCounter.create(fail, window));
        }
        switch(checkCounters(counters)) {
            case ANALYSE:
                break;
            case CONTINUE:
                break;
            case RETURN:
                return;
            default:
                throw new IllegalStateException();
        }
    }
    fill(type, counters, 1);
    max = Math.min(maxFail, settings.getWeightedCounterMaxAllowedFailures());
    for (int fail = MathUtils.min(maxFail, settings.getWeightedCounterMinAllowedFailures()); fail <= max; fail++) {
        for (int w = settings.getWeightedCounterMinPassDecrement(); w <= settings.getWeightedCounterMaxPassDecrement(); w++) {
            counters.add(WeightedFailCounter.create(fail, 1, w));
        }
        switch(checkCounters(counters)) {
            case ANALYSE:
                break;
            case CONTINUE:
                break;
            case RETURN:
                return;
            default:
                throw new IllegalStateException();
        }
    }
    fill(type, counters, 2);
    max = Math.min(maxFail, settings.getResettingCounterMaxAllowedFailures());
    for (int fail = MathUtils.min(maxFail, settings.getResettingCounterMinAllowedFailures()); fail <= max; fail++) {
        for (double f = settings.getResettingCounterMinResetFraction(); f <= settings.getResettingCounterMaxResetFraction(); f += settings.getResettingCounterIncResetFraction()) {
            counters.add(ResettingFailCounter.create(fail, f));
        }
        switch(checkCounters(counters)) {
            case ANALYSE:
                break;
            case CONTINUE:
                break;
            case RETURN:
                return;
            default:
                throw new IllegalStateException();
        }
    }
    fill(type, counters, 3);
    for (int count = settings.getPassRateCounterMinAllowedCounts(); count <= settings.getPassRateCounterMaxAllowedCounts(); count++) {
        for (double f = settings.getPassRateCounterMinPassRate(); f <= settings.getPassRateCounterMaxPassRate(); f += settings.getPassRateCounterIncPassRate()) {
            counters.add(PassRateFailCounter.create(count, f));
        }
        switch(checkCounters(counters)) {
            case ANALYSE:
                break;
            case CONTINUE:
                break;
            case RETURN:
                return;
            default:
                throw new IllegalStateException();
        }
    }
    fill(type, counters, 4);
    counters.trimToSize();
    // Score each of a set of standard fail counters against each frame using how
    // close they are to the target.
    final double[] score = new double[counters.size()];
    final double targetPassFraction = settings.getTargetPassFraction();
    final int nThreads = Prefs.getThreads();
    final ExecutorService executor = Executors.newFixedThreadPool(nThreads);
    final LocalList<Future<?>> futures = new LocalList<>(nThreads);
    final Ticker ticker = ImageJUtils.createTicker(failCountData.size(), nThreads);
    IJ.showStatus("Analysing " + TextUtils.pleural(counters.size(), "counter"));
    for (int i = 0; i < failCountData.size(); i++) {
        final FailCountData data = failCountData.unsafeGet(i);
        futures.add(executor.submit(() -> {
            if (IJ.escapePressed()) {
                return;
            }
            // TODO - Ideally this plugin should be run on benchmark data with ground truth.
            // The target could be to ensure all the correct results are fit
            // and false positives are excluded from incrementing the pass counter.
            // This could be done by saving the results from a benchmarking scoring
            // plugin to memory as the current dataset.
            data.initialiseAnalysis(targetPassFraction);
            // Score in blocks and then do a synchronized write to the combined score
            final Thread t = Thread.currentThread();
            final double[] s = new double[8192];
            int index = 0;
            while (index < counters.size()) {
                if (t.isInterrupted()) {
                    break;
                }
                final int block = Math.min(8192, counters.size() - index);
                for (int j = 0; j < block; j++) {
                    final FailCounter counter = counters.unsafeGet(index + j).newCounter();
                    s[j] = data.score(counter);
                }
                // Write to the combined score
                synchronized (score) {
                    for (int j = 0; j < block; j++) {
                        score[index + j] += s[j];
                    }
                }
                index += block;
            }
            ticker.tick();
        }));
    }
    executor.shutdown();
    ConcurrencyUtils.waitForCompletionUnchecked(futures);
    IJ.showProgress(1);
    if (IJ.escapePressed()) {
        IJ.showStatus("");
        IJ.error(TITLE, "Cancelled analysis");
        return;
    }
    IJ.showStatus("Summarising results ...");
    // TODO - check if the top filter is at the bounds of the range
    final int minIndex = SimpleArrayUtils.findMinIndex(score);
    ImageJUtils.log(TITLE + " Analysis : Best counter = %s (Score = %f)", counters.unsafeGet(minIndex).getDescription(), score[minIndex]);
    // Show a table of results for the top N for each type
    final int topN = Math.min(settings.getTableTopN(), score.length);
    if (topN > 0) {
        final byte[] types = type.toArray();
        final byte maxType = types[types.length - 1];
        final TextWindow resultsWindow = createTable();
        for (byte b = 0; b <= maxType; b++) {
            int[] indices;
            // Use a heap to avoid a full sort
            final IntDoubleMinHeap heap = new IntDoubleMinHeap(topN);
            for (int i = 0; i < score.length; i++) {
                if (types[i] == b) {
                    heap.offer(score[i], i);
                }
            }
            if (heap.getSize() == 0) {
                continue;
            }
            indices = heap.getItems();
            // Ensure sorted
            SortUtils.sortIndices(indices, score, false);
            final StringBuilder sb = new StringBuilder();
            try (final BufferedTextWindow tw = new BufferedTextWindow(resultsWindow)) {
                for (int i = 0; i < topN; i++) {
                    sb.setLength(0);
                    final int j = indices[i];
                    sb.append(i + 1).append('\t');
                    sb.append(counters.unsafeGet(j).getDescription()).append('\t');
                    sb.append(score[j]);
                    tw.append(sb.toString());
                }
            }
        }
    }
    // TODO - Save the best fail counter to the current fit configuration.
    IJ.showStatus("");
}
Also used : RollingWindowFailCounter(uk.ac.sussex.gdsc.smlm.results.count.RollingWindowFailCounter) FailCounter(uk.ac.sussex.gdsc.smlm.results.count.FailCounter) ResettingFailCounter(uk.ac.sussex.gdsc.smlm.results.count.ResettingFailCounter) ConsecutiveFailCounter(uk.ac.sussex.gdsc.smlm.results.count.ConsecutiveFailCounter) WeightedFailCounter(uk.ac.sussex.gdsc.smlm.results.count.WeightedFailCounter) PassRateFailCounter(uk.ac.sussex.gdsc.smlm.results.count.PassRateFailCounter) TByteArrayList(gnu.trove.list.array.TByteArrayList) BufferedTextWindow(uk.ac.sussex.gdsc.core.ij.BufferedTextWindow) Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) LocalList(uk.ac.sussex.gdsc.core.utils.LocalList) TextWindow(ij.text.TextWindow) BufferedTextWindow(uk.ac.sussex.gdsc.core.ij.BufferedTextWindow) IntDoubleMinHeap(uk.ac.sussex.gdsc.core.trees.heaps.IntDoubleMinHeap) ExecutorService(java.util.concurrent.ExecutorService) Future(java.util.concurrent.Future)

Example 2 with FailCounter

use of uk.ac.sussex.gdsc.smlm.results.count.FailCounter 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);
}
Also used : Rectangle(java.awt.Rectangle) BufferedWriter(java.io.BufferedWriter) ImageExtractor(uk.ac.sussex.gdsc.core.utils.ImageExtractor) FailCounter(uk.ac.sussex.gdsc.smlm.results.count.FailCounter) SelectedResultStore(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter.SelectedResultStore) IOException(java.io.IOException) AreaStatistics(uk.ac.sussex.gdsc.core.filters.AreaStatistics) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) IDirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.IDirectFilter) MultiPathFilter(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter) MultiFilter(uk.ac.sussex.gdsc.smlm.results.filter.MultiFilter) MaximaSpotFilter(uk.ac.sussex.gdsc.smlm.filters.MaximaSpotFilter) MultiPathFilter(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter) IMultiPathFitResults(uk.ac.sussex.gdsc.smlm.results.filter.IMultiPathFitResults)

Example 3 with FailCounter

use of uk.ac.sussex.gdsc.smlm.results.count.FailCounter in project GDSC-SMLM by aherbert.

the class FitEngineConfiguration method getFailCounter.

/**
 * Gets the fail counter.
 *
 * @return the fail counter
 */
public FailCounter getFailCounter() {
    // manually passed in. Either would allow more complex fail counters to be used.
    if (failCounter == null) {
        final int failuresLimit = getFailuresLimit();
        final FailCounter f1 = // Negative will disable
        (failuresLimit >= 0) ? ConsecutiveFailCounter.create(failuresLimit) : null;
        final double passRate = getPassRate();
        // TODO - the allowed counts could be an input
        final FailCounter f2 = (passRate > 0) ? PassRateFailCounter.create(5, passRate) : null;
        // All fail counters must pass to continue fitting
        failCounter = CombinedAndFailCounter.join(f1, f2);
        if (failCounter == null) {
            failCounter = NullFailCounter.INSTANCE;
        }
    }
    return failCounter.newCounter();
}
Also used : ConsecutiveFailCounter(uk.ac.sussex.gdsc.smlm.results.count.ConsecutiveFailCounter) CombinedAndFailCounter(uk.ac.sussex.gdsc.smlm.results.count.CombinedAndFailCounter) FailCounter(uk.ac.sussex.gdsc.smlm.results.count.FailCounter) PassRateFailCounter(uk.ac.sussex.gdsc.smlm.results.count.PassRateFailCounter) NullFailCounter(uk.ac.sussex.gdsc.smlm.results.count.NullFailCounter)

Aggregations

FailCounter (uk.ac.sussex.gdsc.smlm.results.count.FailCounter)3 ConsecutiveFailCounter (uk.ac.sussex.gdsc.smlm.results.count.ConsecutiveFailCounter)2 PassRateFailCounter (uk.ac.sussex.gdsc.smlm.results.count.PassRateFailCounter)2 TByteArrayList (gnu.trove.list.array.TByteArrayList)1 TextWindow (ij.text.TextWindow)1 Rectangle (java.awt.Rectangle)1 BufferedWriter (java.io.BufferedWriter)1 IOException (java.io.IOException)1 ExecutorService (java.util.concurrent.ExecutorService)1 Future (java.util.concurrent.Future)1 AreaStatistics (uk.ac.sussex.gdsc.core.filters.AreaStatistics)1 BufferedTextWindow (uk.ac.sussex.gdsc.core.ij.BufferedTextWindow)1 Ticker (uk.ac.sussex.gdsc.core.logging.Ticker)1 IntDoubleMinHeap (uk.ac.sussex.gdsc.core.trees.heaps.IntDoubleMinHeap)1 ImageExtractor (uk.ac.sussex.gdsc.core.utils.ImageExtractor)1 LocalList (uk.ac.sussex.gdsc.core.utils.LocalList)1 Statistics (uk.ac.sussex.gdsc.core.utils.Statistics)1 MaximaSpotFilter (uk.ac.sussex.gdsc.smlm.filters.MaximaSpotFilter)1 CombinedAndFailCounter (uk.ac.sussex.gdsc.smlm.results.count.CombinedAndFailCounter)1 NullFailCounter (uk.ac.sussex.gdsc.smlm.results.count.NullFailCounter)1