use of 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
*
* @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 };
}
use of 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.
* @return
*/
private double calculateDrift(int[] imageT, float scale, double[] dx, double[] dy, double[] originalDriftTimePoints, double smoothing, int iterations, final ImageProcessor[] images, FloatProcessor reference, boolean includeCurrentDrift) {
// Align
tracker.status("Aligning images");
final AlignImagesFFT aligner = new AlignImagesFFT();
aligner.init(reference, WindowMethod.NONE, false);
final Rectangle alignBounds = AlignImagesFFT.createHalfMaxBounds(reference.getWidth(), reference.getHeight(), reference.getWidth(), reference.getHeight());
List<double[]> alignments = Collections.synchronizedList(new ArrayList<double[]>(images.length));
List<Future<?>> futures = new LinkedList<Future<?>>();
int imagesPerThread = getImagesPerThread(images);
for (int i = 0; i < images.length; i += imagesPerThread) {
futures.add(threadPool.submit(new ImageAligner(aligner, images, imageT, alignBounds, alignments, i, i + imagesPerThread)));
}
Utils.waitForCompletion(futures);
tracker.progress(1);
// Used to flag when an alignment has failed
originalDriftTimePoints = Arrays.copyOf(originalDriftTimePoints, originalDriftTimePoints.length);
double[] newDx = new double[dx.length];
double[] newDy = new double[dy.length];
int ok = 0;
for (double[] result : alignments) {
int t = (int) result[2];
if (Double.isNaN(result[0])) {
// TODO: 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) {
double d1 = dx[t] - newDx[t];
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;
}
Aggregations