cophus/scanning-drift-corr

Correction of nonlinear and linear scanning probe artifacts using orthogonal scan pairs.

scanning-drift-corr

Correction of nonlinear and linear scanning probe artifacts using orthogonal scan pairs.

This method will also work on most time series, though it does not yet implement the fast scan direction correction algorithm.

Running the code in Matlab

Running these command lines should reproduce the included example data in matlab.

First, load the data:

load('data_examples/nonlinear_drift_correction_synthetic_dataset_for_testing.mat')

Next, verify the input data has horizontal fast scan directions (i.e. most drift artifacts occur in the vertical / row direction). This can also be easily done using windowed FFT images.

figure('Name','Measured Data 0deg');imagesc(image00deg); axis equal off; colormap(gray(256)); figure('Name','Measured Data 90deg');imagesc(image90deg); axis equal off; colormap(gray(256)); Zoom in on these images, note the row-by-row jumps and overall drift.

The first SPmerge step is to initialize the alignment structure, and perform a search over linear drift vectors (TODO: make this much faster via FFT shear).

sMerge = SPmerge01linear([0 90],image00deg,image90deg);

Note that the alignment is overall quite poor and shows strong Moire artifacts. The yellow cross shows the reference position for the "wrinkle smoothing" initial alignment step. We should move it to a region of the image with good alignment before running the next step. For example:

sMerge.ref = [560 450];

UPDATE - This new version of the code does not perform the "wrinkle smoothing" alignment steps by default! In order to run some of these iterations, you must specify how many iterations you wish to perform using the third input argument. For example, try:

sMerge = SPmerge02(sMerge,0,8);

Much better! In the previous step, we ran 0 "final alignment" steps and 8 "wrinkle smoothing" steps. We can now do the main alignment:

sMerge = SPmerge02(sMerge);

The above command is equivalent to using sMerge = SPmerge02(sMerge,32); or sMerge = SPmerge02(sMerge,32,0); Assuming the alignment has converged correctly, we can now generate a final output image:

imageFinal = SPmerge03(sMerge);

Other notes - I've found that the new linear drift search finds a good initial alignment for most experimental datasets. This is the motivation behind defaulting to 0 iterations for the "wrinkle smoothing" step. However some long dwell time experiments will have slowly varying drift vectors, and so I've left this algorithm in this code.

The vast majority of bad alignments happen when the scan lines are not horizontal, or the images have a large amound of fast scan error (requires a different algorithm to correct). Most orthogonal image pairs can be well-aligned by playing around with this setting in SPmerge02.m: originWindowAverage = 4; Setting this value too large will render the algorithm unable to fit local origin shifts, and setting it too small could overfit to noise. The related value for the "wrinkle smoothing" algorithm is: originInitialAverage = size(sMerge.scanLines,1) / 16;

The other important parameters to check are those that change the initial guess for the alignment. sMerge.ref for example, or the linear search steps in SPmerge01.m: sMerge.linearSearch = linspace(-0.08,0.08,1+2*4);

It's also worth pointing out that the algoithm requires enough padding around the images to fit all of the image translations - make sure this value in SPmerge01.m is not too small: paddingScale = (1+1/2);

If you have image pairs that you are having difficulty with, please feel free to email them to me at [email protected] and I would be happy to have a look!

Authors

Colin Ophus

1. Original author of the Matlab code.
2. Wrote the drift correction paper.

Publication

A publication describing this method can be found here:

https://doi.org/10.1016/j.ultramic.2015.12.002

scanning-drift-corr's Issues

Reason for using "symmetric" in ifft?

Hi Colin,

Do you have a specific reason for using the "symmetric" argument in the ifft2 code?

scanning-drift-corr/matlab/SPmerge01linear.m

Line 106 in 2ae74ef

Icorr = ifft2(sqrt(abs(m)).*exp(1i*angle(m)),'symmetric');

Hybrid correlation?

Hiya @cophus ,

Just out of curiosity, where did you discover this algorithm? It looks like some sort of cross correlation, but I haven't been able to dig up any info on it. What is it?

m = fft2(w2.*sMerge.imageTransform(:,:,1)) ...
            .* conj(fft2(w2.*sMerge.imageTransform(:,:,2)));
Icorr  = ifft2(sqrt(abs(m)).*exp(1i*angle(m)),'symmetric')

scanning-drift-corr/matlab/SPmerge01linear.m

Line 106 in 2ae74ef

Icorr = ifft2(sqrt(abs(m)).*exp(1i*angle(m)),'symmetric');

Help with KDE

I'm not used to working with kernel density estimation, and I think it'd save time if @cophus explains what happens in the matlab code.

There are a few relevant functions:

• scipy.stats.gaussian_kde
• scipy.ndimage.gaussian_filter

I think we should be going straight for the kde function, but I'm unsure how to implement it.

@HamishGBrown I see you've made progress as well! I haven't stuck my code online yet, but I've got a slightly more literal python translation (less efficient) than you have up to the first link in this post. I suggest we both keep going independently and then compare drafts :p

Benchmarking

For (at least) the SPMerge01 function, the bottleneck lies in the image filtering.
In the matlab code, this happens here, and the bottleneck is the convolution step.

% Apply KDE
r = max(ceil(sMerge.KDEsigma*3),5);
sm = fspecial('gaussian',2*r+1,sMerge.KDEsigma);
sm = sm / sum(sm(:));
sig = conv2(sig,sm,'same');
count = conv2(count,sm,'same');
sub = count > 0;
sig(sub) = sig(sub) ./ count(sub);
sMerge.imageTransform(:,:,indImage) = sig;

The convolution with a gaussian kernel of size 2*r+1 and std KDESigma is implemented quite a bit faster in scipy and opencv.

I've attached a jupyter notebook with my findings when tried on a 4096x4096 image, with a kernel of size 11x11 and std=0.5.

In short, opencv's is fastest at 449ms. Half a second for such an image is still a very long time, since this function is repeated a lot.

Benchmark_Gaussians.zip