Dave Mitchell's DigitalMicrograph™ Scripting Website

// Script to create an edge direction map, coloured to show the direction of the edges.

// D. R. G. Mitchell, 5th Nov. 2011

// adminnospam@dmscripting.com (remove the nospam to make this address work)

// version:2011105, v1.0

// www.dmscripting.com

// function to create and apply a Gaussian blurring function to an imaage.

// This based on Andrey Chuvlin's similar function. Slight blurring prior to edge

// finding reduces the sensitivity to noise.

// Sourceimg is the image to be filtered; standard deviation is that of the Gaussian

// kernel: 1=minimal blurring, 3=mild blurring, 10=severe blurring.

// Blurring the image before edge detection reduces the sensitivity to noise.

image GaussianConvolution(image sourceimg, number standarddev)

{

// Trap for a standard deviation value of zero entered

if(standarddev==0) return sourceimg

// get the size of the source image. If it is not a power of 2 in dimension

// warp it so that it is

number xsize, ysize, div2size, expandx, expandy, logsize

getsize(sourceimg, xsize, ysize)

expandx=xsize

expandy=ysize

// Check the x axis for power of 2 dimension - if it is not, round up to the next size

// eg if it is 257 pixels round it up to 512.

logsize=log2(xsize)

if(mod(logsize,1)!=0) logsize=logsize-mod(logsize,1)+1

expandx=2**logsize

// Check the y axis for power of 2 dimension - if it is not, round up to the next size

// eg if it is 257 pixels round it up to 512.

logsize=log2(ysize)

if(mod(logsize,1)!=0) logsize=logsize-mod(logsize,1)+1

expandy=2**logsize

// Use the Warp function to stretch the image to fit into the revised dimensions

image warpimg=realimage("",4,expandx, expandy)

warpimg=warp(sourceimg, icol*xsize/expandx, irow*ysize/expandy)

// Create the gaussian kernel using the same dimensions as the expanded image

image kernelimg:=realimage("",4,expandx,expandy)

number xmidpoint=xsize/2

number ymidpoint=ysize/2

kernelimg=1/(2*pi()*standarddev**2)*exp(-1*(((icol-xmidpoint)**2+(irow-ymidpoint)**2)/(2*standarddev**2)))

// Carry out the convolution in Fourier space

compleximage fftkernelimg:=realFFT(kernelimg)

compleximage FFTSource:=realfft(warpimg)

compleximage FFTProduct:=FFTSource*fftkernelimg.modulus().sqrt()

realimage invFFT:=realIFFT(FFTProduct)

// Warp the convoluted image back to the original size

image filter=realimage("",4,xsize, ysize)

filter=warp(invFFT,icol/xsize*expandx,irow/ysize*expandy)

return filter

}

// function to create a colour map showing the direction of the edges found by Sobel filtering.

// The horizontal and vertical derivatives are computed and the gradient and direction are obtained as

// G=sqrt(dx**2+dy**2) and theta=atan2(dy,dx) respectively.

// The edge direction map is then rounded up/down to the nearest angle of 0, 45, 90 and 135 degs.

// The gradient image shows the edges (Sobel filtering). This is converted into a binary image by thresholding. The

// edges are then colour coded according to their orientation in one of the four angles (above).

// Sourceimg is the image to be filtered; threshold is the numerical value of the threshold.

// If set to zero, the mean intensity value of the edgegradient image will be used (usually about right).

image createedgedirectionmap(image sourceimg, number threshold)

{

// use the Sobel filter to compute the horizontal (dy) and vertical (dx) derivatives

number scalex, scaley, xsize, ysize

string unitstring

getscale(sourceimg,scalex, scaley)

getunitstring(sourceimg, unitstring)

getsize(sourceimg, xsize, ysize)

realimage edgegradient=exprsize(xsize,ysize,1)

realimage edgedirection=exprsize(xsize,ysize,1)

realimage dx=exprsize(xsize,ysize,1)

realimage dy=exprsize(xsize,ysize,1)

dx= offset(sourceimg,-1,-1) - offset(sourceimg,1,-1) + \

2*(offset(sourceimg,-1,0) - offset(sourceimg,1,0)) + \

offset(sourceimg,-1,1) - offset(sourceimg,1,1)

dy= offset(sourceimg,-1,-1) - offset(sourceimg,-1,1) + \

2*(offset(sourceimg,0,-1) - offset(sourceimg,0,1)) + \

offset(sourceimg,1,-1) - offset(sourceimg,1,1)

edgegradient=sqrt(dx**2+dy**2)

edgedirection=atan2(dy, dx)

// Convert the -3.142 to +3.142 direction image from radians to degrees, running 0 to 360 (easier to work with)

edgedirection=(edgedirection+pi())*(180/pi())

// round the direction image up/down to 0, 45, 90, 135 degree angles.

// ie any point lying between 337.5 and 22.5 degs is rounded to 0 (angles

// run from 0 to 360 degrees

edgedirection=tert(edgedirection>337.5 || edgedirection<=22.5, 0, edgedirection)

edgedirection=tert(edgedirection>22.5 && edgedirection<=67.5, 45, edgedirection)

edgedirection=tert(edgedirection>67.5 && edgedirection<=112.5, 90, edgedirection)

edgedirection=tert(edgedirection>112.5 && edgedirection<=157.5, 135, edgedirection)

edgedirection=tert(edgedirection>157.5 && edgedirection<=202.5, 0, edgedirection)

edgedirection=tert(edgedirection>202.5 && edgedirection<=247.5, 45, edgedirection)

edgedirection=tert(edgedirection>247.5 && edgedirection<=292.5, 90, edgedirection)

edgedirection=tert(edgedirection>292.5 && edgedirection<=337.5, 135, edgedirection)

// Create a thresholded image on the basis of the gradient image - after thresholding and

// then colour coding according to the edge's direction

// define the colours for the various angles

rgbnumber red=rgb(255,0,0) // 0 degs

rgbnumber green=rgb(0,255,0) // 45 degs

rgbnumber blue=rgb(0,0,255) // 90 degs

rgbnumber yellow=rgb(255,255,0) // 135 degs

// Threshold the gradient image on the supplied threshold value (if 0 then the image's mean value is used)

// and then colour the thresholded image according to the direction

if(threshold==0) threshold=mean(edgegradient)

rgbimage thresholdimg:=rgbimage("",4,xsize, ysize)

thresholdimg=tert(edgegradient>=threshold && edgedirection==0, red, thresholdimg) // 0 degs are red

thresholdimg=tert(edgegradient>=threshold && edgedirection==45, green, thresholdimg) // 45 degs are green

thresholdimg=tert(edgegradient>=threshold && edgedirection==90, blue, thresholdimg) // 90 degs are blue

thresholdimg=tert(edgegradient>=threshold && edgedirection==135, yellow, thresholdimg) // 135 degs are yellow

return(thresholdimg)

}

// Main program

// source the front-most image and the standard deviation for Gaussian blurring

number xsize, ysize, standarddev, threshold

image front:=getfrontimage()

getsize(front, xsize, ysize)

string imgname=getname(front)

if(!getnumber("Select the standard deviation of the Gaussian kernel:",3,standarddev)) exit(0)

// Apply the Gaussina blur (reduce edge finding sensitivity to noise) and display the result

image blur:=GaussianConvolution(front, standarddev)

// Source the threshold value for converting the Sobel filtered image into an edge map.

// Note: when the threshold is set to zero the above function will use the edgegradient image's mean

// as the threshold value. The thresholded edges will then be coloured to show the direction thereof.

if(!getnumber("Select the intensity threshold (0 uses image mean):",0,threshold)) exit(0)

image edgegradientmap=createedgedirectionmap(blur, threshold)

setname(edgegradientmap, imgname+" - Edge Gradient Map")

showimage(edgegradientmap)

setwindowposition(edgegradientmap, 172,54)