Skip to main content

Posts

Showing posts from 2018

k-means clustering

In order to clustering data to reduce the noise, we can use the simple k-means clustering algorithm. The idea of k-means is quite simple. Here is the step of the k-means algorithm. 1. Randomly pick samples (depending on how many groups we want to cluster) as the references. 2. Compute the distance between each data point and references. 3. Comparing the distance to each reference, grouping the data points from the shortest distance. 4. Compute the centroid of each group as the new reference. 5. Repeat 2-4, until the centroids are the same with the previous result. Here is the Matlab code: ======================================= % An example for k-means clustering % % Renfong 2018/12/19 % % create test data % There are 3 types of sigal % 1-20, 36-70, 91-100 are the group 1 % 21-35 are group 2 % 71-90 are group 3 xx=0:1:1024; cen=[120, 360, 780]; amp=[100, 60, 55]; sig=[50, 10, 30]; % peak 1+3 for i=1:20     sp(i,:)=amp(1)*exp((-...

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...

Image binng in matlab

Here are my own image binning codes. The first one was built on 2015/3/2. And the second one was written on 2018/11/5. I think I have a big improvement. Ha~ ==== function [ox,oy]=binning(x,y,nbins) % This function is used to bin data to average %  % x: x-axis data % y: y-axis data % nbins: binning factor % % 2015/03/02  % Renfong m=max(size(x)); n0=fix(m/nbins); n1=mod(m,nbins); if n1==0 ox=zeros(n0,1); oy=zeros(n0,1); for ii=1:n0 for jj=1:nbins ox(ii)=ox(ii)+x((ii-1)*nbins+jj); oy(ii)=oy(ii)+y((ii-1)*nbins+jj); end ox(ii)=ox(ii); oy(ii)=oy(ii); end else ox=zeros(n0+1,1); oy=zeros(n0+1,1); for ii=1:n0 for jj=1:nbins ox(ii)=ox(ii)+x((ii-1)*nbins+jj); oy(ii)=oy(ii)+y((ii-1)*nbins+jj); end ox(ii)=ox(ii); oy(ii)=oy(ii); end for ii=1:n1 ox(n0+1)=ox(n0+1)+x(n0*nbins+ii); oy(n0+1)=oy(n0+1)+y(n0*nbins+ii); end ox(n0+1)=ox(n0+1)*nbins/n1; oy(n0+1)=oy(n0+1)*nbins/n1; end ==== function out=xy_bins(img,nbins) ...

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(:))...

Non-local means filter --> Image denoise

There are several denoise filters such as Gaussian, bilateral and so on[1]. However, they may deteriorate the image quality. Non-local means filter is a powerful (but time-consuming >"<) denoise filter, which can preserve the detail of the image by averaging neighbor pixels with similar neighborhoods! The Matlab code can be found from the website listed in ref 2. For more detail refer to ref 3. Example: [1] https://web.cs.hacettepe.edu.tr/~erkut/bil717.s12/w09-bilateral-nlmeans.pdf [2] https://pastebin.com/28GsvWkW [3] https://www.youtube.com/watch?v=9tUns4HYtcw

Some test for Top-hat filter

1. Peak standard deviation:     Set the peak std for 1, 2.5, 5, 10 and 20. After std=10, the signals are difficult to identify. 2. Noise:     Set different noise for 0%, 25%, 50%, 75% and 100%. The peaks are hard to identify from 75% noise. 3. Different THF signal window:     Set the noise to 100% and use the different signal window to identify the signal. The result shows the wider window can improve the peak identification. Summary: The THF can help us to identify the weak signal from the high exponential background. It is noted that the S/N ratio would affect the THF result. If we want to identify the rare element in the material, we may use the high S/N ratio spectrum for applying THF or using the large signal window for noisy data.

The difficulty of electron tomography

The difficulties limit the direct application to electron tomography The tilt range has to be from -90 degree to +90 degree The number of projections in a tilt series needs to be 2N for an N*N object The grid points past the resolution circle cannot be experimentally determined These limitations are overcome by combining the PPFFT with an iterative process. Ref: http://www.physics.ucla.edu/research/imaging/EST/index.html

Mask creation

// A practice to create mask by built-in function // iradius, itheta, tert image img:=RealImage("Mask",4,256,256) // for STEM equivalent detector image BF, ABF, LAADF, HAADF BF=img*0 BF=(iradius<10) BF.SetName("BF") ABF=img*0 ABF=(iradius>=10 & iradius<20) ABF.SetName("ABF") LAADF=img*0 LAADF=(iradius>=60 & iradius <80) LAADF.SetName("LAADF") HAADF=img*0 HAADF=(iradius>=100 & iradius<135) HAADF.SetName("HAADF") BF.ShowImage() ABF.ShowImage() LAADF.ShowImage() HAADF.ShowImage() // for DPC detector, 4 segment // set rotate angle number ang=15 if(ang>=90) { ang=0 } ang=pi()*ang/180 // Note that the itheta is varied counterclockwise from 0 to -pi() when y>0  //                                                        and clockwise from 0 to +pi() when y<0  ...

Creating new orthogonal cell via VESTA

The VESTA can be downloaded via http://jp-minerals.org/vesta/en/download.html Original cell Si, space group:Fd-3m, a=5.43. The cell would be looked like... Transform For new cell new a axis: original [1 -1 0] direction new b axis: original [1 1 -2] direction new c axis: original [1 1 1] direction Result new cell: a* projection: [1-10] zone axis b* projection: [11-2] zone axis c* projection: [111] zone axis

Ptychography (1) Basic Concept

Coherent Diffractive Imaging [1] Also known as lensless imaging , is based on the retrieval of a complex sample structure (phase) from a measurement of scattered radiation intensity from a coherently illuminated sample out of the imaging plane. [2] Coherent diffractive imaging (CDI) is an alternative technique that promises, in theory, to allow imaging of specimens at a resolution limited only by the wavelength of the illuminating radiation. CDI replaces the imaging optics and their associated imperfections, with a solution to an inverse problem (phase problem) . The inversion step involves determining the complex exit-surface wave function from the recorded intensity, which in the case of diffractive imaging is in the far field.  Reconstruction of the full complex exit wave recovers the amplitude and phase of the specimen transmission function, related to physical properties of the specimen. CDI involves a non-linear inverse problem and the consequent issue of the...

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...

Practice of dm-script build-in function: icol & calibration information

Using icol function to create peak. Gaussian peak I=A*exp((x-mu)^2/(2*sigma^2) where mu is the peak center            sigma is the divergence of the peak -- image img=RealImage("peak",4,1024,1) //  peak 1 parameters // N_K edge number int1=30 number mu1=300 number sig1=5 // peak 2 parameters // O_K edge number int2=150 number mu2=432 number sig2=0.9 number int3=50 number mu3=438 number sig3=2.5 // add noise number noise=0.5 // EELS signal creation img=int1*exp(-1*(icol-mu1)**2/(2*sig1**2)) img+=int2*exp(-1*(icol-mu2)**2/(2*sig2**2)) img+=int3*exp(-1*(icol-mu3)**2/(2*sig3**2)) img+=0.05*exp(10-icol/100)+noise*int1*random() // Write calibration img.ImageSetDimensionOrigin(0,100) img.ImageSetDimensionScale(0,1) img.ImageSetDimensionUnitString(0, "eV" ) img.ImageSetIntensityUnitString( "e-" ) img.setname("Test Spec + " + noise*100 + "% noise") img.showimage() -- Result