Examples with FloatProcessor ij.process.FloatProcessor used on opensource projects

Example 76 with FloatProcessor

use of ij.process.FloatProcessor 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 77 with FloatProcessor

use of ij.process.FloatProcessor in project GDSC-SMLM by aherbert.

the class ImageKernelFilter method setNPasses.

@Override
public void setNPasses(int passes) {
    // Create the kernel from the image
    boolean build = kernelImp.getID() != lastId || settings.method != lastMethod || settings.filter != lastFilter;
    build = build || (settings.method == Settings.METHOD_SPATIAL && kf == null);
    build = build || (settings.method == Settings.METHOD_FHT && ff == null);
    if (build) {
        final Operation operation = Operation.forOrdinal(settings.filter);
        FloatProcessor fp = kernelImp.getProcessor().toFloat(0, null);
        if (settings.method == Settings.METHOD_SPATIAL) {
            if (kf == null || kernelImp.getID() != lastId || settings.zero != lastZero) {
                fp = pad(fp);
                final int kw = fp.getWidth();
                final int kh = fp.getHeight();
                final float[] kernel = (float[]) fp.getPixels();
                kf = (settings.zero) ? new ZeroKernelFilter(kernel, kw, kh) : new KernelFilter(kernel, kw, kh);
            }
            switch(operation) {
                case CONVOLUTION:
                    kf.setConvolution(true);
                    break;
                case CORRELATION:
                    kf.setConvolution(false);
                    break;
                case DECONVOLUTION:
                default:
                    // Spatial filtering does not support anything other than convolution or correlation.
                    ImageJUtils.log("Unsupported operation (%s), default to correlation", operation.getName());
                    kf.setConvolution(false);
                    break;
            }
        } else {
            if (ff == null || kernelImp.getID() != lastId) {
                final int kw = fp.getWidth();
                final int kh = fp.getHeight();
                final float[] kernel = (float[]) fp.getPixels();
                ff = new FhtFilter(kernel, kw, kh);
                ff.initialiseKernel(dataImp.getWidth(), dataImp.getHeight());
            }
            ff.setOperation(operation);
        }
        lastId = kernelImp.getID();
        lastMethod = settings.method;
        lastFilter = settings.filter;
        lastZero = settings.zero;
    }
    ticker = ImageJUtils.createTicker(passes, passes);
}

Example 78 with FloatProcessor

use of ij.process.FloatProcessor in project GDSC-SMLM by aherbert.

the class TraceMolecules method createBilinearPlot.

private FloatProcessor createBilinearPlot(List<double[]> results, int width, int height) {
    final FloatProcessor ip = new FloatProcessor(width, height);
    // Create lookup table that map the tested threshold values to a position in the image
    final int[] xLookup = createLookup(timeThresholds, settings.getMinTimeThreshold(), width);
    final int[] yLookup = createLookup(ddistanceThresholds, settings.getMinDistanceThreshold(), height);
    origX = (settings.getMinTimeThreshold() != 0) ? xLookup[1] : 0;
    origY = (settings.getMinDistanceThreshold() != 0) ? yLookup[1] : 0;
    final int gridWidth = timeThresholds.length;
    final int gridHeight = ddistanceThresholds.length;
    for (int y = 0, prevY = 0; y < gridHeight; y++) {
        for (int x = 0, prevX = 0; x < gridWidth; x++) {
            // Get the 4 flanking values
            final double x1y1 = results.get(prevY * gridWidth + prevX)[2];
            final double x1y2 = results.get(y * gridWidth + prevX)[2];
            final double x2y1 = results.get(prevY * gridWidth + x)[2];
            final double x2y2 = results.get(y * gridWidth + x)[2];
            // Pixel range
            final int x1 = xLookup[x];
            final int x2 = xLookup[x + 1];
            final int y1 = yLookup[y];
            final int y2 = yLookup[y + 1];
            final double xRange = x2 - x1;
            final double yRange = y2 - y1;
            for (int yy = y1; yy < y2; yy++) {
                final double yFraction = (yy - y1) / yRange;
                for (int xx = x1; xx < x2; xx++) {
                    // Interpolate
                    final double xFraction = (xx - x1) / xRange;
                    final double v1 = x1y1 * (1 - xFraction) + x2y1 * xFraction;
                    final double v2 = x1y2 * (1 - xFraction) + x2y2 * xFraction;
                    final double value = v1 * (1 - yFraction) + v2 * yFraction;
                    ip.setf(xx, yy, (float) value);
                }
            }
            prevX = x;
        }
        prevY = y;
    }
    // Convert to absolute for easier visualisation
    final float[] data = (float[]) ip.getPixels();
    for (int i = 0; i < data.length; i++) {
        data[i] = Math.abs(data[i]);
    }
    return ip;
}

Example 79 with FloatProcessor

use of ij.process.FloatProcessor in project GDSC-SMLM by aherbert.

the class DriftCalculator method calculateDriftUsingImageStack.

/**
 * Calculate the drift of images to the reference image. If no reference is provided then produce
 * a combined z-projection. Update the current drift parameters.
 *
 * @param reference the reference
 * @param images The images to align
 * @param fhtImages The images to align (pre-transformed to a FHT)
 * @param blockT The frame number for each image
 * @param dx The X drift
 * @param dy The Y drift
 * @param originalDriftTimePoints Non-zero when the frame number refers to an aligned image frame
 * @param smoothing the smoothing
 * @param iterations the iterations
 * @return The change in the drift (NaN is an error occurred)
 */
private double calculateDriftUsingImageStack(FloatProcessor reference, ImageProcessor[] images, Fht[] fhtImages, int[] blockT, double[] dx, double[] dy, double[] originalDriftTimePoints, double smoothing, int iterations) {
    Ticker ticker = Ticker.createStarted(tracker, images.length, true);
    if (reference == null) {
        // Construct images using the current drift
        tracker.status("Constructing reference image");
        // Build an image using the current drift
        final List<Future<?>> futures = new LinkedList<>();
        ticker = Ticker.createStarted(tracker, images.length * 2L, true);
        final ImageProcessor[] blockIp = new ImageProcessor[images.length];
        final double[] threadDx = new double[images.length];
        final double[] threadDy = new double[images.length];
        for (int i = 0; i < images.length; i++) {
            threadDx[i] = dx[blockT[i]];
            threadDy[i] = dy[blockT[i]];
        }
        final int imagesPerThread = getImagesPerThread(images);
        for (int i = 0; i < images.length; i += imagesPerThread) {
            futures.add(executor.submit(new ImageTranslator(images, blockIp, threadDx, threadDy, i, i + imagesPerThread, ticker, settings.interpolationMethod)));
        }
        ConcurrencyUtils.waitForCompletionUnchecked(futures);
        // Build an image with all results.
        reference = new FloatProcessor(blockIp[0].getWidth(), blockIp[0].getHeight());
        for (final ImageProcessor ip : blockIp) {
            reference.copyBits(ip, 0, 0, Blitter.ADD);
        }
    }
    // Ensure the reference is windowed
    AlignImagesFft.applyWindowSeparable(reference, WindowMethod.TUKEY);
    return calculateDrift(blockT, 1f, dx, dy, originalDriftTimePoints, smoothing, iterations, fhtImages, reference, false, ticker);
}

Example 80 with FloatProcessor

use of ij.process.FloatProcessor in project GDSC-SMLM by aherbert.

the class DriftCalculator method calculateUsingImageStack.

/**
 * Calculates drift using images from a reference stack aligned to the overall z-projection.
 *
 * @param stack the stack
 * @param limits the limits
 * @return the drift { dx[], dy[] }
 */
@Nullable
private double[][] calculateUsingImageStack(ImageStack stack, int[] limits) {
    // Update the limits using the stack size
    final int upperT = settings.startFrame + settings.frameSpacing * (stack.getSize() - 1);
    limits[1] = Math.max(limits[1], upperT);
    // TODO - Truncate the stack if there are far too many frames for the localisation limits
    tracker.status("Constructing images");
    executor = Executors.newFixedThreadPool(Prefs.getThreads());
    // Built an image and FHT image for each slice
    final ImageProcessor[] images = new ImageProcessor[stack.getSize()];
    final Fht[] fhtImages = new Fht[stack.getSize()];
    final List<Future<?>> futures = new LinkedList<>();
    final Ticker ticker = Ticker.createStarted(tracker, images.length, true);
    final int imagesPerThread = getImagesPerThread(images);
    final AlignImagesFft aligner = new AlignImagesFft();
    final FloatProcessor referenceIp = stack.getProcessor(1).toFloat(0, null);
    // We do not care about the window method because this processor will not
    // actually be used for alignment, it is a reference for the FHT size
    aligner.initialiseReference(referenceIp, WindowMethod.NONE, false);
    for (int i = 0; i < images.length; i += imagesPerThread) {
        futures.add(executor.submit(new ImageFhtInitialiser(stack, images, aligner, fhtImages, i, i + imagesPerThread, ticker)));
    }
    ConcurrencyUtils.waitForCompletionUnchecked(futures);
    tracker.progress(1);
    if (tracker.isEnded()) {
        return null;
    }
    final double[] dx = new double[limits[1] + 1];
    final double[] dy = new double[dx.length];
    final double[] originalDriftTimePoints = new double[dx.length];
    final int[] blockT = new int[stack.getSize()];
    for (int i = 0, t = settings.startFrame; i < stack.getSize(); i++, t += settings.frameSpacing) {
        originalDriftTimePoints[t] = 1;
        blockT[i] = t;
    }
    final double smoothing = updateSmoothingParameter(originalDriftTimePoints);
    lastdx = null;
    // For the first iteration calculate drift to the first image in the stack
    // (since the average projection may have a large drift blurring the image)
    double change = calculateDriftUsingImageStack(referenceIp, images, fhtImages, blockT, dx, dy, originalDriftTimePoints, smoothing, settings.iterations);
    if (Double.isNaN(change) || tracker.isEnded()) {
        return null;
    }
    plotDrift(limits, dx, dy);
    ImageJUtils.log("Drift Calculator : Initial drift " + MathUtils.rounded(change));
    for (int i = 1; i <= settings.maxIterations; i++) {
        change = calculateDriftUsingImageStack(null, images, fhtImages, blockT, dx, dy, originalDriftTimePoints, smoothing, settings.iterations);
        if (Double.isNaN(change)) {
            return null;
        }
        plotDrift(limits, dx, dy);
        if (converged(i, change, getTotalDrift(dx, dy, originalDriftTimePoints))) {
            break;
        }
    }
    if (tracker.isEnded()) {
        return null;
    }
    plotDrift(limits, dx, dy);
    return new double[][] { dx, dy };
}