Skip to main content

LatticeEnhancedFilterV2


This is a band-enhanced filter.
ref: Ondrej L. Krivanek et. al., Nature 464, 571-574

============================================================

// Lattice enhanced filter
// Reference : Ondrej L. Krivanek et. al., Nature 464, 571-574
//
// Version 2
// The filter enhance the lattice feature without removing background.
// Set the weight factor in the background and high frequency region  = 1
// And the weight factor of the enhanced feature  = 3.5
//
// The first window will show the radial average of the FFT of the image
// enhanced feature (look for the peak) then key the value into the window
//
// 2021/03/15
// Renfong
// windless@gmail.com



// sub-function1 : radial average
// This script is retrived and  modified from 
// http://www.gatan.com/sites/default/files/Scripts/Rotational%20Average.s
image Rotational_Average(image img)
{
// Define neccessary parameters and constants
number xscale, yscale, xsize, ysize
number centerx, centery, halfMinor
number scale = img.ImageGetDimensionScale(0)
string unit = img.ImageGetDimensionUnitString(0)

// Likewise, declare intermediate images
image rotational_average, dst, line_projection
// If the source image is complex, take the modulus
if ( img.ImageIsDataTypeComplex( ))
img := modulus(img)

// Get the dimension sizes, and determine half the smallest dimension 
img.Get2dSize( xsize, ysize )
halfMinor = min( xsize, ysize )/2

// Find the centre of the image
centerx = xsize / 2
centery = ysize / 2

// Convert the image to polar co-ordinates...
dst := RealImage( "dst", 4, halfMinor, xsize )
dst = warp( img, icol*sin(irow*2*pi()/xsize) + \
centerx, icol*cos(irow*2*pi()/xsize) + centery )

// and create a line projection using the icol intrinsic variable, 
// normalising with the sampling density
line_projection := RealImage( "line projection", 4, halfMinor, 1 )
line_projection = 0
line_projection[icol,0] += dst
line_projection /= xsize
line_projection.ImageCopyCalibrationFrom(img)
return line_projection
}

// sub-function2 : Gaussian blur
// http://www.dmscripting.com/files/Gaussian_Blur.s
image GaussianBlur(image sourceimg, number standarddev)
{
if(standarddev<=0) return sourceimg
number xsize, ysize
getsize(sourceimg, xsize, ysize)

// Create the gaussian kernel using the same dimensions as the expanded image
image kernelimg:=realimage("",4,xsize,ysize)
kernelimg=1/(2*pi()*standarddev**2)*exp(-1*(iradius**2/(2*standarddev**2)))


// Carry out the convolution in Fourier space
compleximage fftkernelimg:=realFFT(kernelimg)
compleximage FFTSource:=realfft(sourceimg)
compleximage FFTProduct:=FFTSource*fftkernelimg.modulus().sqrt()
realimage invFFT:=realIFFT(FFTProduct)
invFFT=(invFFT-min(invFFT))/(max(invFFT)-min(invFFT))
invFFT=(max(sourceimg)-min(sourceimg))*invFFT+min(sourceimg)
return invFFT
}


// Main function
image LatticeEnhancedFilter(image img, number mu, number f0)
// LatticeEnhancedFilter(image img, number mu, number f0)
// img : input image
//  mu : lattice component (unit: pixel)
//  f0 : background region (unit: pixel)
{
// Set parameters
number sig=sqrt((mu-f0)**2/2.50553)
number sx, sy
getsize(img,sx,sy)
// Creat mask
///////////////////////////////////////////////////////////////
image filimg := img*0
filimg=(1/(sqrt(2*pi()*sig**2)))*exp(-0.5*((iradius-mu)/sig)**2)
// Set background region weight=1
number v0=(1/(sqrt(2*pi()*sig**2)))*exp(-0.5*((f0-mu)/sig)**2)
filimg=tert(iradius<f0,v0,filimg)
filimg=tert(iradius>(2*mu-f0),v0,filimg)
filimg=filimg/v0
// smooth by Gaussian blur
image filimg2=GaussianBlur(filimg,f0/2)
filimg2=filimg2/filimg2.GetPixel(sx/2,sy/2)
filimg2.SetName("Mask")
image LP_filimg2=Rotational_Average(filimg2)
LP_filimg2.SetName("Mask Rotational Profile")
// LP_filimg2.showimage()

// Apply mask
image fft_mask_img=RealFFT(img)*filimg2
image mask_img=RealIFFT(fft_mask_img)
return mask_img
}

// main script
image img:=getfrontimage()
image fft_img=RealFFT(img)
number dp_sx, dp_sy
fft_img.GetSize(dp_sx,dp_sy)
image LP1

//  Display the DP radial average
if(dp_sx>1024)
{
LP1=Rotational_Average(fft_img)[0,0,1,round(dp_sx/4)]
}
else
{
LP1=Rotational_Average(fft_img)
}

LP1[0,0,1,8]=0
LP1.setname("radial average of DP")
LP1.showimage()

// parameter input
number f0, mu
GetNumber("Select enhanced center",30,mu)
f0=10
result("mu="+mu+", f0="+f0+"\n")
LP1.DeleteImage()
image LE_img=LatticeEnhancedFilter(img,mu,f0)


// copy calibration info
string str=img.GetName()
LE_img.setname(str+"Lattice-Enhanced Mask2")
LE_img.ImageCopyCalibrationFrom(img)
LE_img.showimage()

Comments

Post a Comment

Popular posts from this blog

Top hat filter

The top_hat filter can be used to detect the relatively small edges/peaks superimposed on large background signals. The concept came from the EELS workshop during IMC19. Thanks to Prof. Nestor J. Zaluzec. -- // Using Top_hat digital filter to detect the  relatively small edges  //    superimposed on large background signals. // // ref: Ultramicroscopy 18 (1985) 185-190  //      Digital Filters  for Application to Data Analysis in EELS //      by Nestor J. ZALUZEC // Parameters: // win_s: signal window (default:3) // win_b: background window (default:3) //  a_s : amplitude of signal (fixed value) //  a_b : amplitude of background  (fixed value) // Renfong 2018/10/11 // Main function image Top_Hat_Filter(image img, number win_s, number win_b) { // read image string fname=img.GetName() number sx,sy img.getsize(sx,sy) // filter image img2 := imageclone(img)*0 //the area between...
Post a Comment
Read more

HyperSpy - read the calibration information in a dm3/dm4 file

Some example of dm3 file reading by using Python HyperSpy package, which can read the detail information of the dm file. -- # import packages import numpy as np import hyperspy.api as hs # load file sp=hs.load('sp.dm3') # Read the axis information      # Print all the calibration detail print(sp.axes_manager) ''' <Axes manager, axes: (272|2042)>             Name |   size |  index |  offset |   scale |  units  ================ | ======= | ====== | ======= | ======= | ======                     x |    272 |      0 |       -0 |  0.0025 |     µm   --------------- |  ------ | ----- |  ------ | ------- | ------    Energy loss |  2042 |         | 3.2e+02 |       1 |     eV...
Post a Comment
Read more

MLLS in matlab

MLLS stands for  multiple linear least squares fitting, which is the common strategy for the solving EELS edge overlapping and which is also built-in the GMS software. The target spectrum Y and the reference spectrum X Y = A * X Assuming Y is 1*256 matrix and we have three reference spectrums, ie, X is 3*256 matrix. So A is 1*3 matrix. The target is to solve A. If Y and X are n*n matrices, we can use the simple formula Y * inv(X) = A * X * inv(X), ie., A = Y * inv(X). However, Y and X are not n*n  matrices, it is necessary to have some trick to solve it. We can multiply the transpose matrix to produce n*n matrix. Y * X' = A * X * X'  (ps X' means the transpose matrix of X) so A = Y * X' * inv(X * X') Here is the Matlab code: =========  % create target spectrum x=0:256; c=[90,120,155]; sig=[5,10,8]; int=[5,10,8]; xn=zeros(size(x)); ref=zeros(length(c),length(x)); factor=rand(size(c))'; for i=1:length(c)     xn=xn+int(i)*ex...
Post a Comment
Read more