Humboldt-Universität zu Berlin - Faculty of Mathematics and Natural Sciences - Strukturforschung / Elektronenmikroskopie


// Test Script for testing fftw3 routines in Transforms.dll
// Author: Christoph T. Koch
// Max Planck Institut für Metallforschung
// Stuttgart, Germany
// Email:
// This dll provides DM wrapper routines for FFTW3 (single and double precision)
// See for details on these very fast fft routines
// FFTW's main feature of interest here is that
// input/output arrays may be of any size (not just powers of 2)
// Only a tiny subset of the FFTW functionality is provided here,
// keeping the usefulness/effort-ratio at maximum.

void DoShiftImageCenter(complexImage &img, number shift)
complexImage imgBuffer
number width, height, w2s, h2s, w2l, h2l

getSize(img, width, height)

w2s = floor(width/2)
h2s = floor(height/2)
w2l = ceil(width/2)
h2l = ceil(height/2)

imgBuffer = img

img[h2s,w2s,height,width] = imgBuffer[ 0, 0, h2l, w2l] // tl to lr
img[h2s, 0,height, w2s] = imgBuffer[ 0,w2l, h2l,width] // tr to ll
img[0, w2s, h2s,width] = imgBuffer[h2l, 0,height, w2l] // ll to tr
img[0, 0, h2s, w2s] = imgBuffer[h2l,w2l,height,width] // lr to tl
img[h2l,w2l,height,width] = imgBuffer[ 0, 0, h2s, w2s] // tl to lr
img[h2l, 0,height, w2l] = imgBuffer[ 0,w2s, h2s,width] // tr to ll
img[0, w2l, h2l,width] = imgBuffer[h2s, 0,height, w2s] // ll to tr
img[0, 0, h2l, w2l] = imgBuffer[h2s,w2s,height,width] // lr to tl


number sizeX=400;
number sizeY=300;
number shiftX = 103;
number shiftY = 84;
number noise = 0.5;
number mval,mx,my;

// toggle the verbose mode, default os OFF (i.e. no messages).
// T_toggleVerbose();
// T_toggleVerbose();

Number version = T_GetVersion();
Result("Version of Transform: "+version+"\n");
if (version<1.00) {
OkDialog( "This script does not run with the current Version of Transform.dll." )

// first: 3 complex images and a few variables
ComplexImage inImg = CreateComplexImage("real space Image",sizeX,sizeY);
ComplexImage ifftImg = CreateComplexImage("inv. fft'd Image",sizeX,sizeY);
ComplexImage outImg = CreateComplexImage("reciprocal space Image",sizeX,sizeY);

// Do the same for double precision:
ComplexImage d_inImg := ComplexImage("real space Image (double prec.)",16,sizeX,sizeY);
ComplexImage d_ifftImg := ComplexImage("inv. fft'd Image (double prec.)",16,sizeX,sizeY);
ComplexImage d_outImg := ComplexImage("reciprocal space Image (double prec.)",16,sizeX,sizeY);

// create a test input image, something that produces diffraction spots, but
// also some diffuse signal in the FFT:
inImg = (cos(2*Pi()*25*icol/sizeX)*sin(2*Pi()*15*irow/sizeY)+PoissonRandom(0.1))*exp(-0.00005*((icol-sizeX/2)*(icol-sizeX/2)+(irow-sizeY/2)*(irow-sizeY/2)));
inImg = inImg*inImg;
setName(inImg,"original Image");

// The code below demonstrates the use of the complex <-> complex fft routines
// Within the Transforms DLL a flag indicates whether FFTW3 should first try to
// find the fastest possible algorithm ("fftw3SetMeasure"), or just guess
// ("fftw3SetEstimate") which one that might be.

// call the FFT routine. Note that the DC component (q=(0,0)) is at pixel (0,0)
// We therefore must also call ShiftCenter to shift the center of the FFT to the center
// of the output image, if we want to show the fft in the usual way.
setName(outImg,"FFT of Image");

// The fftw routines work for single as well as double precision arrays:
d_inImg = inImg[];
setName(d_outImg,"FFT of Image (double prec.)");
setName(d_inImg,"Image (double prec.)");

// We may make a faster FFTW3-plan by letting FFTW3 test different
// versions of its algorithm. Note that using a different array size
// will also require a new fftw3_plan, estimating, or measuring the performance,
// depending on the current setting.
// Using fftw3SetMeasure will make the first fft to be computed rather slowly, since
// fftw3 will first try which method is optimal for this array size.
// however, the next fft with the same array size will be computed in the fastest
// way possible.
// The actual performance optimization is not done until the fftw3 routine is
// called again. It will be done independently for each fftw-verson
// (e.g. forward- and backward fft will be optimized independently).
// T_fft_SetMeasure();
// However, I have seen a few numerical errors, when using T_fft_SetMeasure()
// I therefore discourage to use this option, suggetsing to
// ALWAYS USE T_fftSetEstimate()!!!

// do the inverse FFT, but this time optimize the performance first
// IMPORTANT: The inverse fft routine does not normalize, it is therefore necessary
// to normalize (i.e. divide by number of pixels) the image afterwards.
ifftImg = ifftImg/(sizeX*sizeY);
setName(ifftImg,"ifft of fft of Image");
ifftImg = ifftImg/(sizeX*sizeY);
d_inImg = d_inImg/(sizeX*sizeY);
setName(d_inImg,"ifft of fft of Image (double precision)");

// showimage(ifftImg);
DoShiftImageCenter(outImg, 1);
DoShiftImageCenter(d_outImg, 1)
// ShiftCenter(outImg); // only works for single precision

// Test cross correlation
// General correlation funtion:
// allows computation of cross-, mutual-, and phase correlation
// see: R.R. Meyer et al, Ultramicrocopy 92 (2002), 89-109
// Input: image1, image2, corr.image, exponent, damping coefficient
// image1, image2, and the pre-allocated correlation image must be real images of equal dimensions
// The exponent (p) is applied in reciprocal space (see paper by Meyer et al)
// Here a few cases with special names:
// exponent=1: cross correlation
// exponent=0: phase correlation
// exponent=0.5: mutual correlation
// The damping factor (c) defines a Gaussian in reciprocal space with which the product of img1* and img2*
// is multiplied to lower the effect of high-frequency noise.
// F'(q) = |F(q)|^p*exp(i*phi(q))*exp(-((c*q_x)^2+(c*q_y)^2)), q_x,q_y=-0.5..0.5
// Here F(q)=FFT(img1)*conj(FFT(img2))
// p=exponent, c=damping coefficient, q=reciprocalspace coordinate
// F'(r) = ifft(F'(q))
// or: XCorrelateGeneral(img2,img1,F',p,c);

result("Script for Testing different image alignment schemes\n");
result("Input shift: "+shiftX+", "+shiftY+"\n");

// create a few images
RealImage img1 = CreateFloatImage("Image 1",sizeX,sizeY);
RealImage xcorr = CreateFloatImage("Image 1",sizeX,sizeY);
img1 = (cos(2*Pi()*25*icol/sizeX)*sin(2*Pi()*15*irow/sizeY));
img1 = img1*sin(2*Pi()*2.73*sqrt((icol/sizeX-0.5)*(icol/sizeX-0.5)+(irow/sizeY-0.5)*(irow/sizeY-0.5)));
img1 = img1*img1;

// RealImage img2 = T_shiftImage(img1,shiftX,shiftY); // shifts with wrapping edges
RealImage img2 = T_shiftImageNoWrap(img1,shiftX,shiftY); // shifts without wrapping edges

// add some noise to images:
img1 = img1+PoissonRandom(noise*img1);
img2 = img2+PoissonRandom(noise*img2);
setName(img2,"image2 = image1 (aside from noise), but shifted");

result("Regular (damped) Cross-correlation has maximum value of "+mval+" at ("+mx+", "+my+")\n");
setName(xcorr,"Cross correlation of image1 and image2");

// Phase correlataion function:
if (1) {
setName(xcorr,"Phase correlation of image1 and image2");
result("Phase correlation has maximum value of "+mval+" at ("+mx+", "+my+")\n");

// Mutual correlation function:
if (1) {
result("Mutual correlation has maximum value of "+mval+" at ("+mx+", "+my+")\n");
setName(xcorr,"Mutual correlation of image1 and image2");


// Note that this DLL also includes the HoughTransform of square images
// by Vincent Hou (need square image to work properly)
// Image newImg = T_HoughTransform(inImg);
// showimage(newImg);