forked from subangstrom/aci
-
Notifications
You must be signed in to change notification settings - Fork 0
/
find_atomic_columns.m
executable file
·130 lines (115 loc) · 4.73 KB
/
find_atomic_columns.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
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)
%sign=1 for adf
%sign=-1 for abf
%style=1 to use the correlated image for peak finding
%style=2 to use the original image for peak finding
warning('off','optimlib:lsqncommon:SwitchToLargeScale');
[ImageX, ImageY]=size(raw_image);
%find the best template used for searching peaks (atomic columns)
if(sigma==0)
template_search_range=10;
maxC=zeros(template_search_range,1);
for i=1:1:template_search_range
template=sign*fspecial('gaussian',20,i);
C=normxcorr2(template,raw_image);
maxC(i,1)=max(max(C));
end
sigma=find(maxC==max(maxC));
end
fprintf('gaussian template sigma=%d\n',sigma);
% use the default threshold to detect peaks
if(threshold==0)
threshold=0.75;
end
for i=1:1:1
template=sign*fspecial('gaussian',19,sigma);
C=normxcorr2(template,raw_image);%ser.data{1,i});
maxC=max(max(C));
D=C(10:10+ImageX-1,10:10+ImageY-1); % this might be a different way to get the peak position
%D=ImageSum;
C(C<threshold*maxC)=0;
C(C>threshold*maxC)=1;
image_peak=raw_image.*C(10:10+ImageX-1,10:10+ImageY-1);
if verbose ==1
figure;
colormap(gray);
imagesc(image_peak);
axis image;
end
CC=bwconncomp(image_peak);
areas_in_pixels = cellfun(@length, CC.PixelIdxList);
fprintf('total peaks found: %d areas range from %d to %d\n',CC.NumObjects,min(areas_in_pixels),max(areas_in_pixels));
fprintf('peaks with areas larger than 10 pixels: %d\n',length(areas_in_pixels(areas_in_pixels>10)));
fprintf('peaks with areas larger than 20 pixels: %d\n',length(areas_in_pixels(areas_in_pixels>20)));
area_size=[1:1:1000];
h_area=hist(areas_in_pixels,area_size);
if verbose ==1
figure;
plot(h_area);
end
% use the centroid first, see if we can improve the result later using
% fitting
centroid = regionprops(CC, 'centroid');
%area_threshold=250;
expand=0; %expand the area to include more pixels
% xpos=zeros(CC.NumObjects,1);
% ypos=zeros(CC.NumObjects,1);
% intensity=zeros(CC.NumObjects,1);
image1=zeros(ImageX,ImageY);
image2=zeros(ImageX,ImageY);
indexj=1;
for j=1:1:min(CC.NumObjects,max_peak_num)
if(areas_in_pixels(j)<area_threshold(1) || areas_in_pixels(j)>area_threshold(2) || areas_in_pixels(j)>1000 || areas_in_pixels(j)<10)
continue;
end
peak_index=CC.PixelIdxList{j};
peak_x=mod(peak_index,ImageX);
peak_x(peak_x==0)=ImageX;
peak_y=floor((peak_index-1)/ImageX)+1;
% peak_x_min=max(min(peak_x)-expand,1);
% peak_x_max=min(max(peak_x)+expand,ImageX);
% peak_y_min=max(min(peak_y)-expand,1);
% peak_y_max=min(max(peak_y)+expand,ImageY);
% fy=peak_x_min:1:peak_x_max;
% fx=peak_y_min:1:peak_y_max;
%peak_data=ser.data{1,1}(peak_x_min:peak_x_max,peak_y_min:peak_y_max);
%[peak_x_min,peak_x_max,peak_y_min,peak_y_max]
fx=peak_x;
fy=peak_y;
fz=peak_x;
for pixels=1:1:length(peak_index)
if style==1
fz(pixels)=sign*D(fx(pixels),fy(pixels));
end
if style==2
fz(pixels)=sign*raw_image(fx(pixels),fy(pixels));
end
end
% if style==1
% peak_data=sign*D(peak_x_min:peak_x_max,peak_y_min:peak_y_max);
% end
% if style==2
% peak_data=sign*raw_image(peak_x_min:peak_x_max,peak_y_min:peak_y_max);
% end
% fz=peak_data;
if ~iscell(initial_values)
[fitresult{indexj},zfit{indexj}, fiterr(indexj,1:7), zerr{indexj}, resnorm(indexj), rr(indexj),StartPoint{indexj}] = fmgaussfit_improved(fy,fx,fz,0,fit_shift);
else
[fitresult{indexj},zfit{indexj}, fiterr(indexj,1:7), zerr{indexj}, resnorm(indexj), rr(indexj),StartPoint{indexj}] = fmgaussfit_improved(fy,fx,fz,initial_values{indexj},fit_shift);
end
mass_center(indexj,1:2)=centroid(j).Centroid;
object_index(indexj)=j;
indexj=indexj+1;
for pixels=1:1:length(peak_index)
image1(fx(pixels),fy(pixels))=fz(pixels);
image2(fx(pixels),fy(pixels))=zfit{indexj-1}(pixels);
end
oimage{indexj}=fz;
%fprintf('find peak number %d in frame %d\n',j,i);
end
% figure;
% imagesc(image1);
% figure;
% imagesc(image2);
end
end