Scripting Resources for DigitalMicrograph™

banner

Dave Mitchell's DigitalMicrograph™ Scripting Website

Home | Scripts | Examples | Functions | Recent Updates | Tutorials | Resources | Publications | Consulting | Projects | Contact & Bio | Private | Search

 

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)