Examples with AlignImagesFft uk.ac.sussex.gdsc.core.ij.AlignImagesFft used on opensource projects

Example 1 with AlignImagesFft

use of uk.ac.sussex.gdsc.core.ij.AlignImagesFft 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 };
}

Example 2 with AlignImagesFft

use of uk.ac.sussex.gdsc.core.ij.AlignImagesFft in project GDSC-SMLM by aherbert.

the class DriftCalculator method calculateDrift.

/**
 * Calculate the drift of images to the reference image. Update the current drift parameters.
 *
 * @param imageT The frame number for each image
 * @param scale The image scale (used to adjust the drift to the correct size)
 * @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 LOESS smoothing parameter
 * @param iterations LOESS iterations parameter
 * @param images The images to align
 * @param reference The reference image
 * @param includeCurrentDrift Set to true if the input images already have the current drift
 *        applied. The new drift will be added to the current drift.
 * @param ticker the ticker
 * @return the change in the drift
 */
private double calculateDrift(int[] imageT, float scale, double[] dx, double[] dy, double[] originalDriftTimePoints, double smoothing, int iterations, final ImageProcessor[] images, FloatProcessor reference, boolean includeCurrentDrift, Ticker ticker) {
    // Align
    tracker.status("Aligning images");
    final AlignImagesFft aligner = new AlignImagesFft();
    aligner.initialiseReference(reference, WindowMethod.NONE, false);
    final Rectangle alignBounds = AlignImagesFft.createHalfMaxBounds(reference.getWidth(), reference.getHeight(), reference.getWidth(), reference.getHeight());
    final List<double[]> alignments = Collections.synchronizedList(new ArrayList<double[]>(images.length));
    final List<Future<?>> futures = new LinkedList<>();
    final int imagesPerThread = getImagesPerThread(images);
    for (int i = 0; i < images.length; i += imagesPerThread) {
        futures.add(executor.submit(new ImageAligner(aligner, images, imageT, alignBounds, alignments, i, i + imagesPerThread, ticker, settings.subPixelMethod)));
    }
    ConcurrencyUtils.waitForCompletionUnchecked(futures);
    tracker.progress(1);
    // Used to flag when an alignment has failed
    originalDriftTimePoints = Arrays.copyOf(originalDriftTimePoints, originalDriftTimePoints.length);
    final double[] newDx = new double[dx.length];
    final double[] newDy = new double[dy.length];
    int ok = 0;
    for (final double[] result : alignments) {
        final int t = (int) result[2];
        if (Double.isNaN(result[0])) {
            // Q: How to ignore bad alignments?
            // Only do smoothing where there was an alignment?
            originalDriftTimePoints[t] = 0;
            tracker.log("WARNING : Unable to align image for time %d to the overall projection", t);
        } else {
            ok++;
            newDx[t] = result[0] / scale;
            newDy[t] = result[1] / scale;
            if (includeCurrentDrift) {
                // New drift = update + previous drift
                newDx[t] += dx[t];
                newDy[t] += dy[t];
            }
        }
    }
    if (ok < 2) {
        tracker.log("ERROR : Unable to align more than 1 image to the overall projection");
        return Double.NaN;
    }
    // Store the pure drift values for plotting
    calculatedTimepoints = Arrays.copyOf(originalDriftTimePoints, originalDriftTimePoints.length);
    lastdx = Arrays.copyOf(newDx, newDx.length);
    lastdy = Arrays.copyOf(newDy, newDy.length);
    // Perform smoothing
    if (smoothing > 0) {
        tracker.status("Smoothing drift");
        if (!smooth(newDx, newDy, originalDriftTimePoints, smoothing, iterations)) {
            return Double.NaN;
        }
    }
    // Interpolate values for all time limits
    tracker.status("Interpolating drift");
    interpolate(newDx, newDy, originalDriftTimePoints);
    // Average drift correction for the calculated points should be zero to allow change comparison
    normalise(newDx, originalDriftTimePoints);
    normalise(newDy, originalDriftTimePoints);
    // Calculate change and update the input drift parameters
    double change = 0;
    for (int t = 0; t < dx.length; t++) {
        if (originalDriftTimePoints[t] != 0) {
            final double d1 = dx[t] - newDx[t];
            final double d2 = dy[t] - newDy[t];
            change += Math.sqrt(d1 * d1 + d2 * d2);
        }
        // Update all points since interpolation has already been done
        dx[t] = newDx[t];
        dy[t] = newDy[t];
    }
    tracker.status("");
    return change;
}

Example 3 with AlignImagesFft

use of uk.ac.sussex.gdsc.core.ij.AlignImagesFft in project GDSC-SMLM by aherbert.

the class PsfCreator method align2.

/**
 * Align the PSFs with the combined PSF using the uk.ac.sussex.gdsc.core.ij.AlignImagesFFT class.
 *
 * @param combined the combined
 * @param psfs the psfs
 * @return The XYZ translations for each PSF
 */
@SuppressWarnings("unused")
private float[][] align2(ExtractedPsf combined, final ExtractedPsf[] psfs) {
    final int n = psfs.length * 3;
    final List<Future<?>> futures = new LocalList<>(n);
    final AlignImagesFft[] align = new AlignImagesFft[3];
    final Rectangle[] bounds = new Rectangle[3];
    for (int i = 0; i < 3; i++) {
        align[i] = new AlignImagesFft();
        final FloatProcessor fp1 = combined.getProjection(i, true);
        final FloatProcessor fp2 = psfs[0].getProjection(i, true);
        align[i].initialiseReference(fp1, WindowMethod.TUKEY, true);
        bounds[i] = AlignImagesFft.createHalfMaxBounds(fp1.getWidth(), fp1.getHeight(), fp2.getWidth(), fp2.getHeight());
    }
    final float[][] results = new float[psfs.length][3];
    for (int j = 0; j < psfs.length; j++) {
        final int jj = j;
        for (int i = 0; i < 3; i++) {
            final int ii = i;
            futures.add(threadPool.submit(() -> {
                final ExtractedPsf psf = psfs[jj];
                final double[] result = align[ii].align(psf.getProjection(ii, true), WindowMethod.TUKEY, bounds[ii], SubPixelMethod.CUBIC);
                // We just average the shift from each projection. There should be
                // two shifts for each dimension
                results[jj][Projection.getXDimension(ii)] -= result[0] / 2;
                results[jj][Projection.getYDimension(ii)] -= result[1] / 2;
            // psfs[index].show(TITLE + index);
            }));
        }
    }
    ConcurrencyUtils.waitForCompletionUnchecked(futures);
    return results;
}