-
Notifications
You must be signed in to change notification settings - Fork 5
/
ACI_flow.m
146 lines (127 loc) · 4.9 KB
/
ACI_flow.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
%% load the RevSTEM image. The gray-scale image is loaded to the variable 'ImageSum_R'
% the RevSTEM image has not been filtered
% note that our code can also work on regular STEM images
% the 'serReader.m' script can be used to read a converntiaonl STEM image
% acquired by TIA
load example.mat
%% display the RevSTEM image
figure;
imagesc(ImageSum_R);
axis image;
colormap(gray);
%% find the atom column locations using normalized cross correlation data
% function [fitresult,oimage,zfit, fiterr, zerr, resnorm,...
% rr,image1,image2,object_index,mass_center,C,D,StartPoint,h_area]=...
% find_atomic_columns(raw_image,sigma,threshold,max_peak_num,sign,style,...
% area_threshold,initial_values,fit_shift,verbose);
% the main output 'fitresult_N' has a cell structure,
% each cell is a 1x7 array [amp, ang, sx, sy, xo, yo, zo] containing the peak fitting result
% amp: amplitude
% ang: rotation angle of the two main axes
% sx: sigma along the first main axis
% sy: sigma along the second main axis
% xo: x coordinate of the peak center
% yo: y coordinate of the peak center
% zo: background intensity
% Meaning of inputs of 'find_atom_clolumns'
% raw_image: The input STEM image
% sigma: sigma of the gaussian distribution for normalized cross-correlation,
% when the number is 0, the program decides the best sigma for ncc
% threshold: threshold to separate the atom columns
% max_peak_num: the limit of peak numbers
% sign: 1: find the peaks; 2: find the valleys (for ABF)
% style: 1: use ncc data to fit the peak; 2: use experimental data
% area_threshold: only areas within the area_threshold range are used for fitting
% initial_values: starting values for peak fitting, can be set as 0
% fit_shift: for future use, 0
% verbose: show the ncc threshold map and area size histogram
[fitresult_N,oimage,zfit, fiterr, zerr, resnorm, rr,image1,image2,object_index,mass_center_N]=...
find_atomic_columns(ImageSum_R,0,0.1,6000,1,1,[100 400],0,0,1);
%% find the atom column locations using experimental RevSTEM data
% style is set to 2 to use experimental data for fitting
% starting points are set to be fitresult_N
[fitresult_E,oimage_E,zfit_E, fiterr_E, zerr_E, resnorm_E, rr_E,image1_E,image2_E,startpoint_E,mass_center_E]=...
find_atomic_columns(ImageSum_R,0,0.1,6000,1,2,[100 400],fitresult_N,0,0);
%% distance histogram calculated from the fitting result
figure;
[xydist_E,h_E]=position_analysis(fitresult_E,ImageSum_R,300,30);
plot(h_E);
%% calculate PSD for the example image
figure;
[Iproj,Iavg,Istd]=project_image_RD(ImageSum_R,100,-90:1:90);
plot(Istd);
%% we can then use the point-PSD to find exactly locations of peaks at roughly 85 degrees and -5 degrees (175 degrees in the plot)
figure;
[d_E,proj_acc_E,projx_E,peak_index_E,row_map_E,col_map_E,mini_E,row_stat_E,col_stat_E,coord_angle_E]...
=assign_xy_to_peaks(ImageSum_R,fitresult_E,200,-5,85,1,1);
plot(projx_E);
%% note the output from matlab says 'find image aligned a long xxx degree at index yyy', use the index yyy to plot the projected profile
figure;
subplot(2,1,1); plot(proj_acc_E(40,:));
subplot(2,1,2); plot(proj_acc_E(143,:));
%% the matrix representation is stored in 'mini_E'
% mini_E is a 2D matrix with each node containing the index of the peak
% fitting result in fitresult_R
figure;
imagesc(mini_E>0);
daspect([1 sqrt(2) 1]);
colormap(gray);
%% calculate the intensity map easily from the matrix representation
col_int_E=get_col_int( ImageSum_R,fitresult_E,0,0);
% the circle around each atom coloumn shows the intensity
figure;
imagesc(ImageSum_R);
axis image
axis off
colormap(gray);
hold on;
[sx,sy]=size(mini_E);
hold all;
jets=jet(256);
upper=max(col_int_E);
lower=min(col_int_E);
for i=1:1:sx
for j=1:1:sy
if mini_E(i,j) == 0
continue;
end
p1=mini_E(i,j);
if col_int_E(p1)==0
continue;
end
color_temp=round((col_int_E(p1)-lower)/(upper-lower)*256);
if(color_temp<1) color_temp=1; end
if(color_temp>256) color_temp=256; end
plot(fitresult_E{p1}(5), fitresult_E{p1}(6),'or','color',jets(color_temp,:),'markersize',9,'linewidth',2);
end
end
title('intensity map','fontsize',20);
%% calculate the ratio map from the matrix representation
[ratio_E,ratio_E_matrix,col_int_matrix]=get_ratio(fitresult_E,mini_E,col_int_E,1,1);
figure;
imagesc(ImageSum_R);
axis image
axis off
colormap(gray);
hold on;
[sx,sy]=size(mini_E);
hold all;
jets=jet(256);
upper=max(ratio_E);
lower=min(ratio_E(ratio_E>0));
for i=1:1:sx
for j=1:1:sy
if mini_E(i,j) == 0
continue;
end
p1=mini_E(i,j);
if ratio_E(p1)==0
continue;
end
color_temp=round((ratio_E(p1)-lower)/(upper-lower)*256);
if(color_temp<1) color_temp=1; end
if(color_temp>256) color_temp=256; end
plot(fitresult_E{p1}(5), fitresult_E{p1}(6),'or','color',jets(color_temp,:),'markersize',9,'linewidth',2);
end
end
title('ratio map','fontsize',20);