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

Example 91 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();
    Rectangle bounds = ip.getRoi();
    // Crop to the ROI
    float[] data = ImageConverter.getData(ip);
    int width = bounds.width;
    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);
        AverageFilter filter = new AverageFilter();
        //filter.blockAverage(smoothData, width, height, smooth);
        if (smooth <= border)
            filter.stripedBlockAverageInternal(smoothData, width, height, (float) smooth);
        else
            filter.stripedBlockAverage(smoothData, width, height, (float) smooth);
    }
    Sort.sort(maxIndices, smoothData);
    // Show the candidate peaks
    if (maxIndices.length > 0) {
        String message = String.format("Identified %d peaks", maxIndices.length);
        if (isLogProgress()) {
            IJ.log(message);
            for (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) {
            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 IJTablePeakResults(showDeviations, imp.getTitle() + " [" + imp.getCurrentSlice() + "]");
    results.begin();
    // Perform the Gaussian fit
    long ellapsed = 0;
    if (!singleFit) {
        if (isLogProgress())
            IJ.log("Combined fit");
        // Estimate height from smoothed data
        double[] estimatedHeights = new double[maxIndices.length];
        for (int i = 0; i < estimatedHeights.length; i++) estimatedHeights[i] = smoothData[maxIndices[i]];
        FitConfiguration config = new FitConfiguration();
        setupPeakFiltering(config);
        long time = System.nanoTime();
        double[] params = fitMultiple(data, width, height, maxIndices, estimatedHeights);
        ellapsed = System.nanoTime() - time;
        if (params != null) {
            // Copy all the valid parameters into a new array
            double[] validParams = new double[params.length];
            int c = 0;
            int validPeaks = 0;
            validParams[c++] = params[0];
            double[] initialParams = convertParameters(fitResult.getInitialParameters());
            double[] paramsDev = convertParameters(fitResult.getParameterStdDev());
            Rectangle regionBounds = new Rectangle();
            int[] xpoints = new int[maxIndices.length];
            int[] ypoints = new int[maxIndices.length];
            int nMaxima = 0;
            for (int i = 1, n = 0; i < params.length; i += 6, n++) {
                int y = maxIndices[n] / width;
                int x = maxIndices[n] % width;
                // Check the peak is a good fit
                if (filterResults && config.validatePeak(n, initialParams, params) != FitStatus.OK)
                    continue;
                if (showFit) {
                    // Copy the valid parameters
                    validPeaks++;
                    for (int ii = i, j = 0; j < 6; ii++, j++) validParams[c++] = params[ii];
                }
                double[] peakParams = extractParams(params, i);
                double[] peakParamsDev = extractParams(paramsDev, i);
                addResult(bounds, regionBounds, data, peakParams, peakParamsDev, nMaxima, x, y, data[maxIndices[n]]);
                // Add fit result to the overlay - Coords are updated with the region offsets in addResult
                double xf = peakParams[3];
                double yf = peakParams[4];
                xpoints[nMaxima] = (int) (xf + 0.5);
                ypoints[nMaxima] = (int) (yf + 0.5);
                nMaxima++;
            }
            setOverlay(nMaxima, xpoints, ypoints);
            // Draw the fit
            if (showFit && validPeaks != 0) {
                double[] pixels = new double[data.length];
                EllipticalGaussian2DFunction f = new EllipticalGaussian2DFunction(validPeaks, width, height);
                invertParameters(validParams);
                f.initialise(validParams);
                for (int x = 0; x < pixels.length; x++) pixels[x] = f.eval(x);
                FloatProcessor fp = new FloatProcessor(width, height, pixels);
                // Insert into a full size image
                FloatProcessor fp2 = new FloatProcessor(ip.getWidth(), ip.getHeight());
                fp2.insert(fp, bounds.x, bounds.y);
                Utils.display(TITLE, fp2);
            }
        } else {
            if (isLogProgress()) {
                IJ.log("Failed to fit " + Utils.pleural(maxIndices.length, "peak") + getReason(fitResult));
            }
            imp.setOverlay(null);
        }
    } else {
        if (isLogProgress())
            IJ.log("Individual fit");
        int nMaxima = 0;
        int[] xpoints = new int[maxIndices.length];
        int[] ypoints = new int[maxIndices.length];
        // Extract each peak and fit individually
        ImageExtractor ie = new ImageExtractor(data, width, height);
        float[] region = null;
        Gaussian2DFitter gf = createGaussianFitter(filterResults);
        for (int n = 0; n < maxIndices.length; n++) {
            int y = maxIndices[n] / width;
            int x = maxIndices[n] % width;
            long time = System.nanoTime();
            Rectangle regionBounds = ie.getBoxRegionBounds(x, y, singleRegionSize);
            region = ie.crop(regionBounds, region);
            int newIndex = (y - regionBounds.y) * regionBounds.width + x - regionBounds.x;
            if (isLogProgress()) {
                IJ.log("Fitting peak " + (n + 1));
            }
            double[] peakParams = fitSingle(gf, region, regionBounds.width, regionBounds.height, newIndex, smoothData[maxIndices[n]]);
            ellapsed += System.nanoTime() - time;
            // Output fit result
            if (peakParams != null) {
                double[] peakParamsDev = null;
                if (showDeviations) {
                    peakParamsDev = convertParameters(fitResult.getParameterStdDev());
                }
                addResult(bounds, regionBounds, data, peakParams, peakParamsDev, n, x, y, data[maxIndices[n]]);
                // Add fit result to the overlay - Coords are updated with the region offsets in addResult
                double xf = peakParams[3];
                double yf = peakParams[4];
                xpoints[nMaxima] = (int) (xf + 0.5);
                ypoints[nMaxima] = (int) (yf + 0.5);
                nMaxima++;
            } else {
                if (isLogProgress()) {
                    IJ.log("Failed to fit peak " + (n + 1) + getReason(fitResult));
                }
            }
        }
        // Update the overlay
        if (nMaxima > 0)
            setOverlay(nMaxima, xpoints, ypoints);
        else
            imp.setOverlay(null);
    }
    results.end();
    if (isLogProgress())
        IJ.log("Time = " + (ellapsed / 1000000.0) + "ms");
}

Example 92 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
	 * 
	 * @param limits
	 * @return the drift { dx[], dy[] }
	 */
private double[][] calculateUsingImageStack(ImageStack stack, int[] limits) {
    // Update the limits using the stack size
    int upperT = startFrame + frameSpacing * (stack.getSize() - 1);
    limits[1] = FastMath.max(limits[1], upperT);
    // TODO - Truncate the stack if there are far too many frames for the localisation limits
    tracker.status("Constructing images");
    threadPool = 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()];
    List<Future<?>> futures = new LinkedList<Future<?>>();
    progressCounter = 0;
    totalCounter = images.length;
    int imagesPerThread = getImagesPerThread(images);
    final AlignImagesFFT aligner = new AlignImagesFFT();
    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.init(referenceIp, WindowMethod.NONE, false);
    for (int i = 0; i < images.length; i += imagesPerThread) {
        futures.add(threadPool.submit(new ImageFHTInitialiser(stack, images, aligner, fhtImages, i, i + imagesPerThread)));
    }
    Utils.waitForCompletion(futures);
    tracker.progress(1);
    if (tracker.isEnded())
        return null;
    double[] dx = new double[limits[1] + 1];
    double[] dy = new double[dx.length];
    double[] originalDriftTimePoints = new double[dx.length];
    int[] blockT = new int[stack.getSize()];
    for (int i = 0, t = startFrame; i < stack.getSize(); i++, t += frameSpacing) {
        originalDriftTimePoints[t] = 1;
        blockT[i] = t;
    }
    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, iterations);
    if (Double.isNaN(change) || tracker.isEnded())
        return null;
    plotDrift(limits, dx, dy);
    Utils.log("Drift Calculator : Initial drift " + Utils.rounded(change));
    for (int i = 1; i <= maxIterations; i++) {
        change = calculateDriftUsingImageStack(null, images, fhtImages, blockT, dx, dy, originalDriftTimePoints, smoothing, 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 };
}

Example 93 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
	 * @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
	 * @param 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) {
    progressCounter = 0;
    totalCounter = images.length;
    if (reference == null) {
        // Construct images using the current drift
        tracker.status("Constructing reference image");
        // Built an image using the current drift
        List<Future<?>> futures = new LinkedList<Future<?>>();
        totalCounter = images.length * 2;
        final ImageProcessor[] blockIp = new ImageProcessor[images.length];
        double[] threadDx = new double[images.length];
        double[] threadDy = new double[images.length];
        for (int i = 0; i < images.length; i++) {
            threadDx[i] = dx[blockT[i]];
            threadDy[i] = dy[blockT[i]];
        }
        int imagesPerThread = getImagesPerThread(images);
        for (int i = 0; i < images.length; i += imagesPerThread) {
            futures.add(threadPool.submit(new ImageTranslator(images, blockIp, threadDx, threadDy, i, i + imagesPerThread)));
        }
        Utils.waitForCompletion(futures);
        // Build an image with all results.
        reference = new FloatProcessor(blockIp[0].getWidth(), blockIp[0].getHeight());
        for (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);
}

Example 94 with FloatProcessor

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

the class PSFCombiner method combineImages.

private void combineImages() {
    double nmPerPixel = getNmPerPixel();
    if (nmPerPixel <= 0)
        return;
    double nmPerSlice = getNmPerSlice();
    if (nmPerPixel <= 0)
        return;
    // Find the lowest start point
    int min = 0;
    for (PSF psf : input) {
        if (min > psf.start)
            min = psf.start;
    }
    // Shift all stacks and find the dimensions
    final int shift = -min;
    int max = 0, size = 0;
    int totalImages = 0;
    for (PSF psf : input) {
        psf.start += shift;
        totalImages += psf.psfSettings.nImages;
        if (max < psf.getEnd())
            max = psf.getEnd();
        if (size < psf.getSize())
            size = psf.getSize();
    }
    // Create a stack to hold all the images
    ImageStack stack = new ImageStack(size, size, max);
    for (int n = 1; n <= max; n++) stack.setPixels(new float[size * size], n);
    // Insert all the PSFs
    IJ.showStatus("Creating combined image ...");
    int imageNo = 0;
    double fraction = 1.0 / input.size();
    for (PSF psf : input) {
        double progress = imageNo * fraction;
        ImageStack psfStack = psf.psfStack;
        final int offsetXY = (psf.getSize() - size) / 2;
        final int offsetZ = psf.start;
        final int w = psf.getSize();
        final double weight = (1.0 * psf.psfSettings.nImages) / totalImages;
        final double increment = fraction / psfStack.getSize();
        for (int n = 1; n <= psfStack.getSize(); n++) {
            IJ.showProgress(progress += increment);
            // Get the data and adjust using the weight
            float[] psfData = ImageConverter.getData(psfStack.getProcessor(n));
            for (int i = 0; i < psfData.length; i++) psfData[i] *= weight;
            // Insert into the combined PSF
            ImageProcessor ip = stack.getProcessor(n + offsetZ);
            ip.copyBits(new FloatProcessor(w, w, psfData, null), offsetXY, offsetXY, Blitter.ADD);
        }
        imageNo++;
    }
    // IJ.showStatus("Normalising ...");
    // PSFCreator.normalise(stack, 1 + shift);
    IJ.showProgress(1);
    IJ.showStatus("");
    ImagePlus imp = Utils.display("Combined PSF", stack);
    imp.setSlice(1 + shift);
    imp.resetDisplayRange();
    imp.updateAndDraw();
    final double fwhm = getFWHM();
    imp.setProperty("Info", XmlUtils.toXML(new PSFSettings(imp.getSlice(), nmPerPixel, nmPerSlice, totalImages, fwhm)));
    Utils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images, FWHM = %s\n", imp.getTitle(), imp.getSlice(), Utils.rounded(nmPerPixel), Utils.rounded(nmPerSlice), totalImages, Utils.rounded(fwhm));
}

Example 95 with FloatProcessor

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

the class PSFCreator method addToPSF.

private boolean addToPSF(int maxz, final int magnification, ImageStack psf, final double[] centre, final float[][] spot, final Rectangle regionBounds, double progress, final double increment, final boolean centreEachSlice) {
    // Calculate insert point in enlargement
    final int z = (int) centre[2];
    int insertZ = maxz - z + 1;
    // Enlargement size
    final int w = regionBounds.width, h = regionBounds.height;
    final int dstWidth = w * magnification;
    final int dstHeight = h * magnification;
    // Multi-thread for speed
    if (threadPool == null)
        threadPool = Executors.newFixedThreadPool(Prefs.getThreads());
    List<Future<?>> futures = new LinkedList<Future<?>>();
    for (int i = 0; i < spot.length; i++) {
        //final int slice = i + 1;
        final ImageProcessor ip = psf.getProcessor(insertZ++);
        final float[] originalSpotData = spot[i];
        futures.add(threadPool.submit(new Runnable() {

            public void run() {
                if (Utils.isInterrupted())
                    return;
                incrementProgress(increment);
                double insertX, insertY;
                // Enlarge
                FloatProcessor fp = new FloatProcessor(w, h, originalSpotData, null);
                fp.setInterpolationMethod(interpolationMethod);
                fp = (FloatProcessor) fp.resize(dstWidth, dstHeight);
                // In the case of Bicubic interpolation check for negative values
                if (interpolationMethod == ImageProcessor.BICUBIC) {
                    float[] pixels = (float[]) fp.getPixels();
                    for (int i = 0; i < pixels.length; i++) if (pixels[i] < 0)
                        pixels[i] = 0;
                }
                // when resizing and the CoM will move.
                if (centreEachSlice) {
                    final double[] com = calculateCenterOfMass(fp);
                    //System.out.printf("CoM %d : %f %f vs %f %f\n", slice, com[0], com[1],
                    //		centre[0] * magnification, centre[1] * magnification);
                    // Get the insert position by subtracting the centre-of-mass of the enlarged image from the 
                    // image centre + allow for a border of 1 pixel * magnification
                    insertX = magnification + dstWidth * 0.5 - com[0];
                    insertY = magnification + dstHeight * 0.5 - com[1];
                //Utils.log("Insert point = %.2f,%.2f => %.2f,%.2f\n", dstWidth * 0.5 - cx, dstHeight * 0.5 - cy,
                //		insertX, insertY);
                } else {
                    // Get the insert position from the stack centre using enlargement
                    insertX = getInsert(centre[0], (int) centre[0], magnification);
                    insertY = getInsert(centre[1], (int) centre[1], magnification);
                //Utils.log("Insert point = %.2f,%.2f => %.2f,%.2f\n", centre[0] - (int) centre[0], centre[1] - (int) centre[1], insertX, insertY);
                }
                // Copy the processor using a weighted image
                final int lowerX = (int) insertX;
                final int lowerY = (int) insertY;
                final double wx2 = insertX - lowerX;
                final double wx1 = 1 - wx2;
                final double wy2 = insertY - lowerY;
                final double wy1 = 1 - wy2;
                // Add to the combined PSF using the correct offset and the weighting
                copyBits(ip, fp, lowerX, lowerY, wx1 * wy1);
                copyBits(ip, fp, lowerX + 1, lowerY, wx2 * wy1);
                copyBits(ip, fp, lowerX, lowerY + 1, wx1 * wy2);
                copyBits(ip, fp, lowerX + 1, lowerY + 1, wx2 * wy2);
            //// Check CoM is correct. This is never perfect since the bilinear weighting 
            //// interpolates the data and shifts the CoM.
            //ImageProcessor ip2 = ip.createProcessor(ip.getWidth(), ip.getHeight());
            //copyBits(ip2, fp, lowerX, lowerY, wx1 * wy1);
            //copyBits(ip2, fp, lowerX + 1, lowerY, wx2 * wy1);
            //copyBits(ip2, fp, lowerX, lowerY + 1, wx1 * wy2);
            //copyBits(ip2, fp, lowerX + 1, lowerY + 1, wx2 * wy2);
            //
            //double[] com = getCoM((FloatProcessor) ip2);
            //System.out.printf("Inserted CoM %d : %f %f\n", slice, com[0], com[1]);
            }
        }));
        if (Utils.isInterrupted())
            break;
    }
    Utils.waitForCompletion(futures);
    return !Utils.isInterrupted();
}