Scripting Resources for DigitalMicrograph™ 
Example: FFT Manipulation 

Function 
Demonstrates various FFT processes, such as Butterworth filtering and computing the inverse FFT, phase, power spectrum etc. 
Version 
version:20070216, v1.0 
Author 
D. R. G. Mitchell 
Acknowledgements 
Thanks to Bernhard Schaffer and Kazuo Ishizuka for tips on using these commands. 
Comments 
To play with this script have an image shown which contains periodic information  such as a HRTEM. 
System Requirements 
Should be compatible with all recent versions of DigitalMicrograph. 
Known Issues 
The image must be 2^n x 2^n in dimension eg 256 x 256, 512 x 512 etc. 
Supported 
Yes 
Included Files 
Main script file. 
Source Code 
// Some example functions which show how to do FFT/IFFT manipulations // Please use/develop/expand and correct these examples at will.
// D. R. G. Mitchell, adminnospam@dmscripting.com (remove the nospam to make this email address work) // v1.0, February 2007 // version:20070216 // Thanks to Bernhard Schaffer and Kazuo Ishizuka for advice and assistance
// To use this script, have an image with periodic information in it foremost. Something like a HRTEM image will // do fine. The image must be square and have dimensions an integral power of 2 eg 512 x 512, 1024 x 1024 etc.
// The script shows how to :
// Calculate the FFT of an image // Calculate the inverse FFT of an FFT // Create a Butterworth filter to use for high frequency filtering // Apply the Butterworth filter to a FFT // Carry out the inverse FFT of the above to show the effect of filtering // Extract the real component of a complex image // Extract the imaginary component of a complex image // Compute the magnitude component of the FFT // Compute the Phase component of the FFT  don't stare at it too long // Compute the Power Spectrum of an image.
// Some functions:
// This function tests the passed in image to make sure it is the right data type, is square and is of dimension 2n x 2n
void testfftcompatibility(image frontimage) { number xsize, ysize, imagetype, modlog2value getsize(frontimage, xsize, ysize)
// Get the datatype of the image imagetype=imagegetdatatype(frontimage)
// Trap for complex or RGB images  these are not compatible with FFTing if(imagetype==3 imagetype==13  imagetype==23) // 3=packed complex, 13=complex 16, 23=RGB { showalert("The foremost image must be of type Integer or Real!",0) exit(0) }
// Checks to make sure that the image dimensions are an integral power of 2 modlog2value=mod(log2(xsize), 1)
if (xsize!=ysize  modlog2value!=0) { showalert("Image dimensions must be an integral power of 2 for FFT!",0) exit(0) } }
// This function carried out the forward FFT. The function used (realFFT()) requires a real image // so a clone of the passed in image is created and converted to a real image
image forwardfft(realimage frontimage) {
// Get some info on the passed in image number xsize, ysize, imagetype string imgname=getname(frontimage) getsize(frontimage, xsize, ysize)
// create a complex image of the correct size to store the result of the FFT compleximage fftimage=compleximage("",8,xsize, ysize)
// Clone the passed in image and convert it to type real (required for realFFT()) image tempimage=imageclone(frontimage) converttofloat(tempimage) fftimage=realfft(tempimage) deleteimage(tempimage)
return fftimage }
// The Butterworth Filter Function  this creates a filter which can be applied to a FFT // to exclude the high frequency component  ie remove noise.
image butterworthfilter(number imgsize, number bworthorder, number zeroradius) { // See John Russ's Image Processing Handbook, 2nd Edn, p 316 image butterworthimg=realimage("",4,imgsize, imgsize) butterworthimg=0
// note the halfpointconst value sets the value of the filter at the halfway point // ie where the radius = zeroradius. A value of 0.414 sets this value to 0.5 // a value of 1 sets this point to root(2)
number halfpointconst=0.414 butterworthimg=1/(1+halfpointconst*(iradius/zeroradius)**(2*bworthorder)) return butterworthimg }
// Function to carry out image processing on the FFT and return the result // The passed in images are the original HRTEM image and a Butterworth filter to remove the high frequency component // Note if the Butterworth image is inverted then the low frequency component is filtered and the high frequencies are retained. // Be aware if the central region of the FFT is removed by a mask, then weird things will happed to // the resulting inverse image. It is better to leave a pinhole in the mask to allow the very lowest // frequencies through. This pinhole should have a gradual edge to avoid ringing. An example of this is // shown in my HRTEM Filter script on this database.
image FFTfiltering(image frontimage, image butterworthimg) { number xsize, ysize getsize(frontimage, xsize, ysize)
// Compute the FFT of passed in image, then mulitply it by the Butterworth filter image compleximage fftimage=forwardfft(frontimage) compleximage maskedfft=fftimage*butterworthimg return maskedfft }
// Main program starts here
// Check to make sure an image is shown
number nodocs nodocs=countdocumentwindowsoftype(5) if(nodocs==0) { showalert("There are no images displayed!",0) exit(0) }
// Get the foremost image and some data from it
image front:=getfrontimage() number xsize, ysize getsize(front, xsize, ysize)
// Position the foremost image
setwindowposition(front, 142, 24) updateimage(front) string imgname=getname(front)
// Check to make sure the image is compatible with FFT
testfftcompatibility(front)
// Carry out the forward FFT
image fftimage=forwardfft(front)
// packed complex images are Hermitian  see DM help for details // Creating packed complex FFTs enables use of the packedIFFT() function  see below // which will return a real image from a FFT. Use of the IFFT() function returns a complex image // and if used, the converttopackedcomplex() line below should be omitted.
converttopackedcomplex(fftimage)
setname(fftimage, "FFT of "+imgname) showimage(fftimage) setwindowposition(fftimage, 172,54)
// Do the inverse FFT  there is an alternative function  IFFT(). However, // this requires a complex, rather than packed complex, image as its argument // and returns the original image as a complex image rather than a real image.
image inversefftimg=packedIFFT(fftimage) showimage(inversefftimg) setname(inversefftimg, "Inverse FFT of "+imgname) setwindowposition(inversefftimg, 202, 84)
// Carry out filtering by masking out parts of the FFT then doing the inverse // In this case a Butterworth filter is used. This selects the low frequency // part of the FFT and filters out the high frequency stuff  removing noise. // The Butterworth filter, has a value of 1 in the central region, then rolls off to a value // of zero near the edge of the FFT. Graduated filters like this are essential in FFT filtering // If sharp cut offs are used, 'ringing' artefacts may appear in the inverse FFT. Changing the // values of Butterworth order (currently=3) and zero radius (currently = xsize/5) will change // the slope of the roll off and the radius at which the filter's value drops to half respectively.
image butterworthimage=butterworthfilter(xsize, 3, xsize/5) showimage(butterworthimage) setname(butterworthimage, "Butterworth Filter") setwindowposition(butterworthimage, 232, 114)
// Apply the Butterworth filter image to the orignal front image // and show the resulting masked FFT of the frontimage
compleximage maskimage=fftfiltering(front, butterworthimage) converttopackedcomplex(maskimage) showimage(maskimage) setname(maskimage, "Butterworth Masked FFT of "+imgname) setwindowposition(maskimage, 262, 144)
// Do the inverse FFT on the masked FFT to illustrate the effect of removing the high frequency component
image invfilteredimg=packedIFFT(maskimage) showimage(invfilteredimg) setname(invfilteredimg, "Butterworth Filtered "+imgname) setwindowposition(invfilteredimg, 292, 174)
// FFTs are complex images containing a real and an imaginary part. FFTs contain two components  magnitude and phase. // The magnitude component is usually displayed as a power spectrum. This shows the relative magnitudes of the frequencies // contributing to the image. The FFT you see in Digital Micrograph is actually the absolute value of the magnitude. The code below // shows how to access the real and imaginary parts of a complex image and compute the magnitude, phase and power spectrum // Note the phase image does not impart any intuitive information  unless you are into psychedelic drugs.
// The magnitude of a FFT = sqt (real component**2 + imaginary component**2) // The phase of a FFT = atan2(real component / imaginary component) // The Power Spectrum of a FFT = magnitude**2
// Extract and display the real part of the complex image (FFT) image realimg=real(fftimage) showimage(realimg) setwindowposition(realimg, 322,204) setname(realimg, "Real part of FFT of "+imgname)
// Extract and display the real part of the complex image (FFT) image imagimg=imaginary(fftimage) showimage(imagimg) setwindowposition(imagimg, 352,234) setname(imagimg, "Imaginary part of FFT of "+imgname)
// Compute the magnitude image image magimg=sqrt(real(fftimage)**2+imaginary(fftimage)**2) showimage(magimg) setwindowposition(magimg, 382, 264) setname(magimg, "Magnitude of FFT of "+imgname)
// Compute the phase image image phaseimg=atan(real(fftimage)/imaginary(fftimage)) showimage (phaseimg) setwindowposition(phaseimg, 412, 294) setname(phaseimg, "Phase of FFT of "+imgname)
// Compute the Power Spectrum image image psimg=magimg**2 showimage(psimg) setwindowposition(psimg, 442, 324) setname(psimg, "Power Spectrum of "+imgname) 