Examples with ImageExtractor uk.ac.sussex.gdsc.core.utils.ImageExtractor used on opensource projects

Example 6 with ImageExtractor

use of uk.ac.sussex.gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class GaussianFit method runFinal.

/**
 * Perform fitting using the chosen maxima. Update the overlay if successful.
 *
 * @param ip The input image
 */
private void runFinal(ImageProcessor ip) {
    ip.reset();
    final Rectangle bounds = ip.getRoi();
    // Crop to the ROI
    final float[] data = ImageJImageConverter.getData(ip);
    final int width = bounds.width;
    final int height = bounds.height;
    // Sort the maxima
    float[] smoothData = data;
    if (getSmooth() > 0) {
        // Smoothing destructively modifies the data so create a copy
        smoothData = Arrays.copyOf(data, width * height);
        final BlockMeanFilter filter = new BlockMeanFilter();
        if (settings.smooth <= settings.border) {
            filter.stripedBlockFilterInternal(smoothData, width, height, (float) settings.smooth);
        } else {
            filter.stripedBlockFilter(smoothData, width, height, (float) settings.smooth);
        }
    }
    SortUtils.sortIndices(maxIndices, smoothData, true);
    // Show the candidate peaks
    if (maxIndices.length > 0) {
        final String message = String.format("Identified %d peaks", maxIndices.length);
        if (isLogProgress()) {
            IJ.log(message);
            for (final int index : maxIndices) {
                IJ.log(String.format("  %.2f @ [%d,%d]", data[index], bounds.x + index % width, bounds.y + index / width));
            }
        }
        // Check whether to run if the number of peaks is large
        if (maxIndices.length > 10) {
            final GenericDialog gd = new GenericDialog("Warning");
            gd.addMessage(message + "\nDo you want to fit?");
            gd.showDialog();
            if (gd.wasCanceled()) {
                return;
            }
        }
    } else {
        IJ.log("No maxima identified");
        return;
    }
    results = new ImageJTablePeakResults(settings.showDeviations, imp.getTitle() + " [" + imp.getCurrentSlice() + "]");
    final CalibrationWriter cw = new CalibrationWriter();
    cw.setIntensityUnit(IntensityUnit.COUNT);
    cw.setDistanceUnit(DistanceUnit.PIXEL);
    cw.setAngleUnit(AngleUnit.RADIAN);
    results.setCalibration(cw.getCalibration());
    results.setPsf(PsfProtosHelper.getDefaultPsf(getPsfType()));
    results.setShowFittingData(true);
    results.setAngleUnit(AngleUnit.DEGREE);
    results.begin();
    // Perform the Gaussian fit
    long ellapsed = 0;
    final FloatProcessor renderedImage = settings.showFit ? new FloatProcessor(ip.getWidth(), ip.getHeight()) : null;
    if (!settings.singleFit) {
        if (isLogProgress()) {
            IJ.log("Combined fit");
        }
        // Estimate height from smoothed data
        final double[] estimatedHeights = new double[maxIndices.length];
        for (int i = 0; i < estimatedHeights.length; i++) {
            estimatedHeights[i] = smoothData[maxIndices[i]];
        }
        final FitConfiguration config = new FitConfiguration();
        setupPeakFiltering(config);
        final long time = System.nanoTime();
        final double[] params = fitMultiple(data, width, height, maxIndices, estimatedHeights);
        ellapsed = System.nanoTime() - time;
        if (params != null) {
            // Copy all the valid parameters into a new array
            final double[] validParams = new double[params.length];
            int count = 0;
            int validPeaks = 0;
            validParams[count++] = params[0];
            final double[] initialParams = convertParameters(fitResult.getInitialParameters());
            final double[] paramsDev = convertParameters(fitResult.getParameterDeviations());
            final Rectangle regionBounds = new Rectangle();
            final float[] xpoints = new float[maxIndices.length];
            final float[] ypoints = new float[maxIndices.length];
            int npoints = 0;
            for (int i = 1, n = 0; i < params.length; i += Gaussian2DFunction.PARAMETERS_PER_PEAK, n++) {
                final int y = maxIndices[n] / width;
                final int x = maxIndices[n] % width;
                // Check the peak is a good fit
                if (settings.filterResults && config.validatePeak(n, initialParams, params, paramsDev) != FitStatus.OK) {
                    continue;
                }
                if (settings.showFit) {
                    // Copy the valid parameters before there are adjusted to global bounds
                    validPeaks++;
                    for (int ii = i, j = 0; j < Gaussian2DFunction.PARAMETERS_PER_PEAK; ii++, j++) {
                        validParams[count++] = params[ii];
                    }
                }
                final double[] peakParams = extractParams(params, i);
                final double[] peakParamsDev = extractParams(paramsDev, i);
                addResult(bounds, regionBounds, peakParams, peakParamsDev, npoints, x, y, data[maxIndices[n]]);
                // Add fit result to the overlay - Coords are updated with the region offsets in addResult
                final double xf = peakParams[Gaussian2DFunction.X_POSITION];
                final double yf = peakParams[Gaussian2DFunction.Y_POSITION];
                xpoints[npoints] = (float) xf;
                ypoints[npoints] = (float) yf;
                npoints++;
            }
            setOverlay(npoints, xpoints, ypoints);
            // Draw the fit
            if (validPeaks != 0) {
                addToImage(bounds.x, bounds.y, renderedImage, validParams, validPeaks, width, height);
            }
        } else {
            if (isLogProgress()) {
                IJ.log("Failed to fit " + TextUtils.pleural(maxIndices.length, "peak") + ": " + getReason(fitResult));
            }
            imp.setOverlay(null);
        }
    } else {
        if (isLogProgress()) {
            IJ.log("Individual fit");
        }
        int npoints = 0;
        final float[] xpoints = new float[maxIndices.length];
        final float[] ypoints = new float[maxIndices.length];
        // Extract each peak and fit individually
        final ImageExtractor ie = ImageExtractor.wrap(data, width, height);
        float[] region = null;
        final Gaussian2DFitter gf = createGaussianFitter(settings.filterResults);
        double[] validParams = null;
        final ShortProcessor renderedImageCount = settings.showFit ? new ShortProcessor(ip.getWidth(), ip.getHeight()) : null;
        for (int n = 0; n < maxIndices.length; n++) {
            final int y = maxIndices[n] / width;
            final int x = maxIndices[n] % width;
            final long time = System.nanoTime();
            final Rectangle regionBounds = ie.getBoxRegionBounds(x, y, settings.singleRegionSize);
            region = ie.crop(regionBounds, region);
            final int newIndex = (y - regionBounds.y) * regionBounds.width + x - regionBounds.x;
            if (isLogProgress()) {
                IJ.log("Fitting peak " + (n + 1));
            }
            final double[] peakParams = fitSingle(gf, region, regionBounds.width, regionBounds.height, newIndex, smoothData[maxIndices[n]]);
            ellapsed += System.nanoTime() - time;
            // Output fit result
            if (peakParams != null) {
                if (settings.showFit) {
                    // Copy the valid parameters before there are adjusted to global bounds
                    validParams = peakParams.clone();
                }
                double[] peakParamsDev = null;
                if (settings.showDeviations) {
                    peakParamsDev = convertParameters(fitResult.getParameterDeviations());
                }
                addResult(bounds, regionBounds, peakParams, peakParamsDev, n, x, y, data[maxIndices[n]]);
                // Add fit result to the overlay - Coords are updated with the region offsets in addResult
                final double xf = peakParams[Gaussian2DFunction.X_POSITION];
                final double yf = peakParams[Gaussian2DFunction.Y_POSITION];
                xpoints[npoints] = (float) xf;
                ypoints[npoints] = (float) yf;
                npoints++;
                // Draw the fit
                if (settings.showDeviations) {
                    final int ox = bounds.x + regionBounds.x;
                    final int oy = bounds.y + regionBounds.y;
                    addToImage(ox, oy, renderedImage, validParams, 1, regionBounds.width, regionBounds.height);
                    addCount(ox, oy, renderedImageCount, regionBounds.width, regionBounds.height);
                }
            } else if (isLogProgress()) {
                IJ.log("Failed to fit peak " + (n + 1) + ": " + getReason(fitResult));
            }
        }
        // Update the overlay
        if (npoints > 0) {
            setOverlay(npoints, xpoints, ypoints);
        } else {
            imp.setOverlay(null);
        }
        // Create the mean
        if (settings.showFit) {
            for (int i = renderedImageCount.getPixelCount(); i-- > 0; ) {
                final int count = renderedImageCount.get(i);
                if (count > 1) {
                    renderedImage.setf(i, renderedImage.getf(i) / count);
                }
            }
        }
    }
    results.end();
    if (renderedImage != null) {
        ImageJUtils.display(TITLE, renderedImage);
    }
    if (isLogProgress()) {
        IJ.log("Time = " + (ellapsed / 1000000.0) + "ms");
    }
}

Example 7 with ImageExtractor

use of uk.ac.sussex.gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class PerPixelCameraModelTest method canCropAndConvertDataWithCropBounds.

@SeededTest
void canCropAndConvertDataWithCropBounds(RandomSeed seed) {
    final PerPixelCameraModelTestData data = (PerPixelCameraModelTestData) dataCache.computeIfAbsent(seed, PerPixelCameraModelTest::createData);
    final PerPixelCameraModel model = new PerPixelCameraModel(w, h, data.bias, data.gain, data.variance);
    final UniformRandomProvider rand = RngUtils.create(seed.getSeed());
    final ImageExtractor ie = ImageExtractor.wrap(data.bias, w, h);
    for (int j = 0; j < 10; j++) {
        final Rectangle bounds = getBounds(rand, ie);
        checkConversion(data, bounds, model.crop(bounds, false));
    }
}

Example 8 with ImageExtractor

use of uk.ac.sussex.gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class PerPixelCameraModelTest method canGetMeanVarianceData.

private static void canGetMeanVarianceData(RandomSeed seed, boolean initialise, boolean normalised) {
    final PerPixelCameraModelTestData data = (PerPixelCameraModelTestData) dataCache.computeIfAbsent(seed, PerPixelCameraModelTest::createData);
    final PerPixelCameraModel model = createModel(data, initialise);
    final UniformRandomProvider rand = RngUtils.create(seed.getSeed());
    final ImageExtractor ie = ImageExtractor.wrap(data.bias, w, h);
    for (int i = 0; i < 10; i++) {
        final Rectangle bounds = getBounds(rand, ie);
        final float[] v = (normalised) ? model.getNormalisedVariance(bounds) : model.getVariance(bounds);
        final double e = MathUtils.sum(v) / v.length;
        final double o = (normalised) ? model.getMeanNormalisedVariance(bounds) : model.getMeanVariance(bounds);
        Assertions.assertEquals(e, o);
    }
}

Example 9 with ImageExtractor

use of uk.ac.sussex.gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class PsfCreator method fitSpot.

private MemoryPeakResults fitSpot(ImageStack stack, final int width, final int height, final int x, final int y) {
    Rectangle regionBounds = null;
    // Create a fit engine
    final MemoryPeakResults results = new MemoryPeakResults();
    final FitConfiguration fitConfig = config.getFitConfiguration();
    results.setCalibration(fitConfig.getCalibration());
    results.setPsf(fitConfig.getPsf());
    results.setSortAfterEnd(true);
    results.begin();
    final int threadCount = Prefs.getThreads();
    final FitEngine engine = FitEngine.create(config, SynchronizedPeakResults.create(results, threadCount), threadCount, FitQueue.BLOCKING);
    for (int slice = 1; slice <= stack.getSize(); slice++) {
        // Extract the region from each frame
        final ImageExtractor ie = ImageExtractor.wrap((float[]) stack.getPixels(slice), width, height);
        if (regionBounds == null) {
            regionBounds = ie.getBoxRegionBounds(x, y, boxRadius);
        }
        final float[] region = ie.crop(regionBounds);
        // Fit only a spot in the centre
        final FitParameters params = new FitParameters();
        params.maxIndices = new int[] { boxRadius * regionBounds.width + boxRadius };
        final ParameterisedFitJob job = new ParameterisedFitJob(slice, params, slice, region, regionBounds);
        // jobItems.add(job);
        engine.run(job);
    }
    engine.end(false);
    results.end();
    return results;
}

Example 10 with ImageExtractor

use of uk.ac.sussex.gdsc.core.utils.ImageExtractor 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);
}