Dave Mitchell's DigitalMicrograph™ Scripting Website

// Cross correlate ROI with image

// Takes a region of interest in an image and uses it as a motif to look for similar regions in the rest

// of the image, using cross correlation. For versions of GMS prior to 2, the image will need to be square

// and a power of 2 in size (eg 1024 x 1024). For later versions of GMS the image can be any shape.

// A small (10-20% of the image size) region of interest must be present on the front-most image.

// This can be any size shape, but the final motif image will be a circular region from within that ROI.

// The radius of the circular region is the minimum of the half the smallest side length of the ROI.

// The motif image excises the ROI, offsets it to be centred and applies a Butterworth filter to

// it so that the edges of the ROI are feathered down to zero to eliminate sharp edge artefacts.

// This feathered motif image is then cross correlated with the original

// image. Regions which appear bright are where the structure in the motif and the image correlate

// strongly. This may be useful for finding defects or distorted regions of lattice images etc.

// D. R. G. Mitchell, adminnospam@dmscripting.com (remove the nospam to make this address work)

// version:20160611, v1.1, June 2016, www.dmscripting.com

// Function to create a Temperature colour lookup table (CLUT)

image CreateTemperatureCLut()

{

image customclut=rgbimage("", 4,256,1)

customclut=tert(icol<51, rgb(0,0,((icol/51)*255)),customclut) // black to blue rgb (0,0,0) -> rgb(0,0,1)

customclut=tert(icol>=51 && icol<102, rgb(0, (((icol-51)/51)*255), 255-(((icol-51)/51)*255)), customclut) // blue to green

customclut=tert(icol>=102 && icol<153, rgb(((icol-102)/51)*255, (255-(((icol-102)/51)*255)), 0), customclut) // green to red

customclut=tert(icol>=153 && icol<204, rgb(255, (((icol-153)/51)*255),0), customclut) // red to yellow

customclut=tert(icol>=204 && icol<=255, rgb(255,255,(((icol-204)/51)*255)), customclut) // yellow to white

return customclut

}

// Function to create a butterworth filter. Imgxsize and imgysize are the sizes of the filter image

// bworthorder is a numerical value (1-6 is good), which defines the rate at which the edge of the filter

// decays to zero. Low values give shallow slopes. zeroradius specifies the radius of the filter.

image butterworthfilter(number imgxsize, number imgysize, number bworthorder, number zeroradius)

{

// See John Russ's Image Processing Handbook, 2nd Edn, p 31

image butterworthimg=realimage("",4,imgxsize, imgysize)

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

}

// Main program starts here

// These two flags set the behaviour of the script

number clipnegativevalues =1 // Set this to 0 (no) or 1 (yes) to clip out all negative values from the final image

number applyclut=1 // set this to 0 (no) or 1 (yes) to apply a colour lookup table to the final image

// Source the frontimage

number nodocs=countdocumentwindowsoftype(5)

if(nodocs<1)

{

showalert("Ensure an image with a small region of interest (ROI) is displayed.",2)

exit(0)

}

image front:=getfrontimage()

// Check that an ROI is present on the image

imagedisplay frontdisp=front.imagegetimagedisplay(0)

number norois=frontdisp.imagedisplaycountrois()

if(norois<1)

{

showalert("Ensure a square ROI defines a region of typical structure.",2)

exit(0)

}

image roi:=getfrontimage()[]

// Get some information on the ROI

roi theroi=frontdisp.imagedisplaygetroi(0)

number roit, roil, roib, roir

theroi.roigetrectangle(roit, roil, roib, roir)

number roicntrx=((roir-roil)/2)+roil

number roicntry=((roib-roit)/2)+roit

number imgx, imgy, roix, roiy

getsize(front, imgx, imgy)

getsize(roi, roix, roiy)

number imgcntrx=imgx/2

number imgcntry=imgy/2

// Move the ROI to the centre of an otherwise empty image

number offsetx=roicntrx-imgcntrx

number offsety=roicntry-imgcntry

image markup=imageclone(front)*0

markup=offset(front, offsetx, offsety)

number minoffset, maxoffset

minmax(markup, minoffset, maxoffset)

markup=(markup-minoffset)/(maxoffset-minoffset)

// Create a butterworth mask for the above image

number butterworthorder=3 // values between 1 and 6 are sensible. 1 produces a very gentle slope to the

// feathering at the edge of the butterworth filter. 3 is a fairly steep roll off and 6 is a very abrupt edge

number roixradius=(roir-roil)/2

number roiyradius=(roib-roit)/2

number zeroradius=min(roixradius, roiyradius)

image bwfilter=butterworthfilter(imgx, imgy, butterworthorder, zeroradius)

// Create a motif image and cross correlate it with the original image

image motif=bwfilter*markup

image frontscale=imageclone(front)

number frontmin, frontmax

minmax(frontscale, frontmin, frontmax)

frontscale=(frontscale-frontmin)/(frontmax-frontmin)

image crosscorr=crosscorrelate(motif, frontscale)

string imgname=getname(front)

setname(crosscorr, "Cross Correlation of "+imgname)

imagecopycalibrationfrom(crosscorr, front)

taggroup fronttags=front.imagegettaggroup()

taggroup cctags=crosscorr.imagegettaggroup()

taggroupcopytagsfrom(cctags, fronttags)

// Rotate the image and display in temperature colour

crosscorr=rotate(crosscorr, pi()) // the cross correlation is rotated with respect to the original

if(clipnegativevalues!=0) crosscorr=tert(crosscorr<0,0,crosscorr) // clip out negatives if the flag is set to 1

// set at the start of the main script

showimage(crosscorr)

imagedisplay ccdisp=crosscorr.imagegetimagedisplay(0)

image clut=CreateTemperatureCLut()

if(applyclut!=0) ccdisp.imagedisplaysetinputcolortable(clut) // colour the image if the flag has been set to 1

// at the start of the program