David Joy, Yeong-Uk Koa and Justin Hwu’s algorithm for determining the resolution of scanning electron beam instruments. The algorithm was written for NIH or SCION image.
This document describes its design and use.
{SMART - Scanning Microscope Analysis and Resolution Testing
(c) v6.0c DCJ December 2000
Thanks to Brendan Griffin, Travis Barroni, Kimeo Kanda and Gian Lorusso for
enhancements and corrections to the code.
Revision v50.b fixes the 'copy' bug in the two image mode
Revision v6.0b adds the J. Frank and L. Al-Ali S/N code
Nature vol 256,p376 (1975) modified for SEM use
Revision v6.0c December 2000 modifies the autocorrelation routine
}
var {global }
i,j,k,m,n,MyPicID,PixCount:integer;
x,y,ROIW,ROIH,HM:integer;
centerx,centery,radius,w1,h1:integer;
Pixperum,width,CutOff,resolution,ratio,angle,height:real;
procedure roi(size: integer); {sets up square ROI}
var
width, height: integer;
begin
GetPicSize(width, height);
if size > width then
size := width;
if width = 0 then begin
PutMessage('No image window open.');
exit;
end;
MakeRoi(width/2 - size/2, height/2 - size/2, size, size);
end;
macro '*** SMART V6.0 gamma ***';
macro '(---'; begin end;
macro 'Make square ROI...'
macro ' 128 x 128 pixel'; begin roi(128) end;
macro ' 256 x 256 pixel'; begin roi(256) end;
macro ' 512 x 512 pixel'; begin roi(512) end;
macro '(---'; begin end;
macro 'Manual Resolution Test';
{use FFT method to measure information cut-off}
var
radX,radY,pic3:integer;
pixsize,rfactor,theta:real;
left,top,m,n:integer;
begin
SetScale(1,'pixel'); {set size unit to pixels}
MyPicID:=PidNumber; {identify front image}
SelectPic(MyPicID);
GetPicSize(w1,h1); {get width of image in pixels}
beep;
{get calibration data on the image}
width:=GetNumber('Image width in um',1.0);
PixCount:=w1;
Pixperum:=PixCount/width;
GetROI(x,y,ROIW,ROIH); {get size of ROI}
centerx:=ROIW/2;
centery:=ROIH/2;
{now do the FFT of this area}
begin
fft('foreward');
end;
{clean up the screen}
SelectPic(MyPicID);
Dispose;
SetThreshold(128); {high contrast image}
{get the calibration factor for the FFT display}
rfactor:=(1000*ROIW)/pixperum;
PutMessage('Adjust Threshold');
{pauses program to allow for adjustment}
{crude hack to create a pseudo-table for index lookup}
rUser2[1]:=0.1;rUser2[2]:=0.2;rUser2[3]:=0.3;
rUser2[4]:=0.5;rUser2[5]:=0.7;rUser2[6]:=1;
rUser2[7]:=1.5;rUser2[8]:=2;rUser2[9]:=2.5;
rUser2[10]:=3;rUser2[11]:=4;rUser2[12]:=5;
rUser2[13]:=7;rUser2[14]:=10;rUser2[15]:=20;
rUser2[16]:=30;rUser2[17]:=50;rUser2[18]:=70;
rUser2[19]:=100;rUser2[20]:=150;
{draw the calibration rings}
for m:=1 to 20 do begin
n:=21-m;
radius:=round((1000*ROIW)/(pixperum*rUser2[n]));
if ((radius<ROIW/2) and (radius>ROIW/10)) then
begin
MakeOvalROI(ROIW/2-radius,ROIH/2-radius, 2*radius,2*radius);
DrawBoundary; {shows the circle}
{now annotate the rings}
theta:=1.57*(1-random); {arbitrary angle}
radX:=round(radius*cos(theta)-1);
radY:=round(radius*sin(theta)-1);
MoveTo(ROIW/2+radX,ROIH/2-radY);
Write(rUser2[n]:1:1); {label the ring}
end;
end;
end;
macro 'AutoMeasure Resolution';
{uses FFT method to measure high frequency information cut-off}
begin
SetScale(1,'pixel'); {set size unit to pixels}
MyPicID:=PidNumber; {identify front image}
SelectPic(MyPicID);
GetPicSize(w1,h1); {get width of image in pixels}
beep;
{get calibration data on the image}
width:=GetNumber('Image width in um',1.0);
PixCount:=w1;
Pixperum:=PixCount/width;
GetROI(x,y,ROIW,ROIH); {get size of ROI}
centerx:=ROIW/2;
centery:=ROIH/2;
{now do the FFT of this area}
begin
fft('foreward');
end;
{clean up the screen}
SelectPic(MyPicID);
Dispose;
SetThreshold(128); {high contrast image}
PutMessage('Adjust Threshold,then click');
{pauses program to allow for adjustment}
SetBinaryCount(4); {maximum number of open/close loops}
MakeBinary; {binary the FFT prior to fitting}
{now clean up the FFT}
i:=1;
while i<=5 do begin {cleaning up the FFT}
Erode;
Dilate;
i:=i+1;
end;
{now perform the analysis}
SetOptions('Major;Minor;Angle'); {set up options in results window}
SetParticleSize(500,999999); {exclude small outlying areas}
AnalyzeParticles; {stereology routine for FFT}
{now compute the resolution - the information limit is....}
CutOff:=(rMajor[1]+rMinor[1])/2; {major axis of ellipse}
ratio:=(rMajor[1]-rMinor[1])/rMajor[1]; {measure of astigmatism}
resolution:=(1000*ROIW)/(Pixperum*CutOff); {in nm}
ShowMessage(' Resolution=',resolution:1:2,' nm','\ Eccentricity =',ratio:1:2);
end;
macro '(---'; begin end;
macro 'AutoCorrelation';
{measure resolution by autocorrelation of the image FFT}
var
pic1,pic2,fft1,fft2,out:integer;
x,y,ROIW,ROIH,wide:integer;
mid,count,ppv,min,max,start:integer;
pkpos,upper,lower,HM:integer;
pixsize,SN:real;
begin
SetScale(1,'pixels'); {set size unit to pixels}
Pic1:=PidNumber; {identify working image}
SelectPic(Pic1);
GetPicSize(w1,h1); {get width of image in pixels}
beep;
{calibrate the image size}
width:=GetNumber('Image width in um',1.0);
PixCount:=w1;
PixPerum:=PixCount/width;
{now get the position of the ROI}
GetROI(x,y,ROIW,ROIH);
mid:=ROIH/2; {y-coordinate for line profile}
start:=1; {start x-coordinate for profile}
wide:=2; {width of line profile}
{do the first FFT}
fft('foreward'); {transform of the image into Fourier space}
fft1:=pidNumber;
ImageMath('cmul',fft1,fft1,1,0,'FFT2'); {now do the autocorrelation}
out:=pidNumber;
{clean up the screen}
SelectPic(fft1);
Dispose;
SelectPic(out);
fft('Inverse');
fft('SwapQuadrants');
SetPicName('AutoCorrelation');
{now get the line profile values across this}
MakeLineRoi(start,mid,ROIW,mid);
GetPlotData(count,ppv,min,max);
{error trap.....}
if min<=0 then min:=1;
{get the cut-off value for the Rayleigh criterion}
HM:=(max-min)/5;
i:=ROIW/2;
{find the width of the peak = upper side first}
while ((PlotData[i]>HM) and (i<ROIW)) do begin
i:=i+1;
end;
upper:=i;
{now the lower side}
j:=(ROIW/2);
while ((PlotData[j]>HM) and (0<j))do begin
j:=j-1;
end;
lower:=j;
{now compute the resolution}
resolution:=(1000*128)/((upper -lower)*pixperum);
ShowMessage('\ resolution=',resolution:1:2,' nm');
end;
macro '(---'; begin end;
macro 'Resolution from Two Images';
{computes resolution from composite of two sequential images}
var
pic1,pic2,pic3,pic4,fft1,fft2,out:integer;
centerx,centery,newx,newy:integer;
mid,count,ppv,min,max,start:integer;
upper,lower,HM:integer;
pixsize,rfactor,theta:real;
x,y,v:integer;
radX,radY,pic3:integer;
left,top,m,n:integer;
begin
if nPics<>2 then begin
PutMessage('Exactly two images required');
exit;
end;
SetScale(1,'pixels');
pic1:=pidNumber; {identify the first image}
GetROI(x,y,ROIW,ROIH);
Copy; {image in the ROI area}
centerx:=ROIW/2;
centery:=ROIH/2;
SetNewSize(ROIW,ROIH);
MakeNewWindow('First'); {copy it to a new image}
Paste;
pic4:=pidNumber; {get id for this new image}
{get calibration data on the image}
SelectPic(1);
GetPicSize(w1,h1); {get width of image in pixels}
width:=GetNumber('Image width in um',1.0);
PixCount:=w1;
Pixperum:=PixCount/width;
pic2:=pidNumber; {and then the second image}
{shift the ROI by 16 pixels}
MakeROI(x+16,y,ROIW,ROIH);
Copy; {to clipboard}
SetNewSize(ROIW,ROIH);
MakeNewWindow('Second');
Paste;
pic3:=pidNumber;
ImageMath('add',pic3,pic4,0.5,0,'Overlay');
out:=pidNumber;
{get the calibration factor for the FFT display}
rfactor:=(1000*ROIW)/pixperum;
{now clean up the screen}
SelectPic(pic1);
Dispose;
SelectPic(pic2);
Dispose;
SelectPic(pic3);
Dispose;
SelectPic(pic4);
Dispose;
{now do the FFT}
begin
fft('foreward');
end;
{clean up the screen}
SelectPic(out);
Dispose;
SetThreshold(128);
PutMessage('Adjust threshold to enhance fringes - then click');
{pauses program to allow for adjustment}
{crude hack to create a pseudo-table for index lookup}
rUser2[1]:=0.1;rUser2[2]:=0.2;rUser2[3]:=0.3;
rUser2[4]:=0.5;rUser2[5]:=0.7;rUser2[6]:=1;
rUser2[7]:=1.5;rUser2[8]:=2;rUser2[9]:=2.5;
rUser2[10]:=3;rUser2[11]:=4;rUser2[12]:=5;
rUser2[13]:=7;rUser2[14]:=10;rUser2[15]:=20;
rUser2[16]:=30;rUser2[17]:=50;rUser2[18]:=70;
rUser2[19]:=100;rUser2[20]:=150;
{draw the calibration rings}
for m:=1 to 20 do begin
n:=21-m;
radius:=round((1000*ROIW)/(pixperum*rUser2[n]));
if ((radius<ROIW/2) and (radius>ROIW/10)) then
begin
MakeOvalROI(ROIW/2-radius,ROIH/2-radius, 2*radius,2*radius);
DrawBoundary; {shows the circle}
{now annotate the rings}
theta:=1.57*(1-random); {arbitrary angle}
radX:=round(radius*cos(theta)-1);
radY:=round(radius*sin(theta)-1);
MoveTo(ROIW/2+radX,ROIH/2-radY);
Write(rUser2[n]:1:1); {label the ring}
end;
end;
end;
macro '(---'; begin end;
macro 'S/N ratio ';
{computes S/N from covariance of sequential line scans
this is a beta version and may need more work}
var
varT1,T1,dummy:real;
varT2,T2,covT1T2:real;
pic1,pic2:integer;
x,y,ROIW,ROIH,wide:integer;
startx,starty,nmax:integer;
SNR,SNRAv:real;
begin
{calibrate in pixels for convenience)
SetScale (1,'pixel');
pic1:=pidNumber; {identify the first image}
SelectPic(pic1);
GetPicSize(w1,h1); {get width of image in pixels}
varT1:=0;
T1:=0;
nmax:=100; {test vakue}
{randomly choose an area to analyze}
startx:=round(w1*random);
starty:=round(h1*random);
{ensure these values are inside image}
if (startx+nmax)>=w1 then startx:=w1-nmax-1; {error trap}
if (starty+nmax)>=h1 then starty:=h1-nmax-1;
for i:=startx to (startx+nmax) do begin {1st line}
T1:=T1+GetPixel(i,starty);
end;
T1:=T1/nmax; {mean value}
{now compute variance}
k:=1;
for i:=startx to (startx+nmax) do begin
rUser1[k]:=GetPixel(i,starty);
varT1:=varT1+(rUser1[k]-T1)*(rUser1[k]-T1);
k:=k+1;
end;
varT1:=varT1/nmax; {variance 1}
{now do the next line}
varT2:=0;
T2:=0;
for i:=startx to (startx+nmax) do begin
T2:=T2+GetPixel(i,starty+1);
end;
T2:=T2/nmax; {mean value}
{now compute variance}
k:=1;
for i:=startx to (startx+nmax) do begin
rUser2[k]:=GetPixel(i,starty+1);
varT2:=varT2+(rUser2[k]-T2)*(rUser2[k]-T2);
k:=k+1;
end;
varT2:=varT2/nmax; {variance 2}
{now compute the covariance}
k:=1;
covT1T2:=0;
for k:=1 to nmax do begin
covT1T2:=covT1T2+(rUser1[k]-T1)*(rUser2[k]-T2);
k:=k+1;
end;
covT1T2:=covT1T2/nmax; {covariance T1 and T2}
{now get the SNR value}
dummy:=sqrt(varT1*varT2);
SNR:=covT1T2/(dummy-covT1T2);
if SNR<0 then ShowMessage('\ Inappropriate Conditions'); {error trap}
SNR:=sqrt(SNR);
{now repeat this process for n lines }
SNRAv:=SNR;
m:=2;
for n:=1 to 50 do begin {the loop of pairs of lines}
T1:=0;
for i:=startx to (startx+nmax) do begin {1st line}
T1:=T1+GetPixel(i,starty+m);
end;
T1:=T1/nmax; {mean value}
{now compute variance}
k:=1;
for i:=startx to (startx+nmax) do begin
rUser1[k]:=GetPixel(i,starty+m);
varT1:=varT1+(rUser1[k]-T1)*(rUser1[k]-T1);
k:=k+1;
end;
varT1:=varT1/nmax; {variance 1}
{now compute the covariance}
k:=1;
covT1T2:=0;
for k:=1 to nmax do begin
covT1T2:=covT1T2+(rUser1[k]-T1)*(rUser2[k]-T2);
k:=k+1;
end;
covT1T2:=covT1T2/nmax; {covariance T1 and T2}
{now get this SNR value}
dummy:=sqrt(varT1*varT2);
SNR:=covT1T2/(dummy-covT1T2);
SNR:=sqrt(SNR);
SNRAv:=SNRAv+SNR;
m:=m+1; {advance line by one pixel}
{now do the second line}
varT2:=0;
T2:=0;
for i:=startx to (startx+nmax) do begin
T2:=T2+GetPixel(i,starty+m);
end;
T2:=T2/nmax; {mean value}
{now compute variance}
k:=1;
for i:=startx to (startx+nmax) do begin
rUser2[k]:=GetPixel(i,starty+m);
varT2:=varT2+(rUser2[k]-T2)*(rUser2[k]-T2);
k:=k+1;
end;
varT2:=varT2/nmax; {variance 2}
{now compute the covariance}
k:=1;
covT1T2:=0;
for k:=1 to nmax do begin
covT1T2:=covT1T2+(rUser1[k]-T1)*(rUser2[k]-T2);
k:=k+1;
end;
covT1T2:=covT1T2/nmax; {covariance T1 and T2}
{now get the SNR value}
dummy:=sqrt(varT1*varT2);
SNR:=covT1T2/(dummy-covT1T2);
SNR:=sqrt(SNR);
SNRAv:=SNRAv+SNR;
m:=m+1; {advance down one line}
ShowMessage('\ still calculating');
end; {of line loop}
SNRAv:=SNRAv/n; {average all the data lines}
ShowMessage('\\ SNR=',SNRAv:1:2);
end; {of this macro}
end;