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

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

Drift correction in Matlab

In order to improve S/N ratio, microscopist uses several short acquisition time images, and then sum them up. So the drift correction is very important. Here is the demo of how to use Matlab do drift correction. In the first, load an image, and using circshift function to shift the object in the image. Then use fft cross correlation to compute the moving distance. Finally, shift the object to the origial position. Here is the testing code: == % Demo of drift correction % 2018/11/15  by Renfong im1= imread('cameraman.tif');    % reference image [sy,sx]=size(im1); % shift the object im2=circshift(im1,[20,10]);    % the object moved down 20 pixels and moved right 10 pixels. figure(1); subplot(121);imshow(im1); subplot(122);imshow(im2); % Using fft cross correlations to detect the moving distance[1] fftim1=fft2(im1); fftim2=fft2(im2); cc=fftshift(ifft2(fftim1.*conj(fftim2))); [shiftY,shiftX]=find(cc==max(cc(:))...
Post a Comment
Read more