-
Notifications
You must be signed in to change notification settings - Fork 8
/
ewpc_mapSpots.m
152 lines (126 loc) · 5.25 KB
/
ewpc_mapSpots.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
147
148
149
150
151
function [ spotMaps ] = ewpc_mapSpots( data4d,spotList,tol,valid)
%ewpc_mapSpots Produces maps of peak fits in CBEDFFT
% input:
% data4d -- the 4D dataset, ordered k1,k2,x1,x2
% spotList -- a struct array identifying the spots to be mapped.
% Fields must be:
% 'id' -- a name identifying the spot.
% 'spotRangeX1' -- index range around the spot in X1.
% 'spotRangeX2' -- index range around the spot in X2.
% tol -- (optional) convergence tolerance. Default is 1e-3.
% valid -- (optional) a mask to identify a subregion in which to
% perform the fitting procedure
% output:
% spotMaps -- struct array containing maps of maximum spot index Q1,Q2
% for each spot in spotList. Fields are:
% 'id' -- a name identifying the spot.
% 'Q1map' -- map of Q1 index value at peak maximum.
% 'Q2map' -- map of Q2 index value at peak maximum.
% 'x1map' -- map of x1 vector component for the spot.
% 'x2map' -- map of x2 vector component for the spot.
% 'spotlength' -- map of the spot vector length.
% 'spotangle' -- map of the spot vector angle in degrees.
%
%This function is part of the PC-STEM Package by Elliot Padgett in the
%Muller Group at Cornell University. Last updated June 26, 2019.
%initialize data info
[N_k1,N_k2,N_x1,N_x2]=size(data4d);
win=window2(N_k1,N_k2,@hann);
minval=min(data4d(:));
if nargin<3
tol=1e-3;
valid = ones(N_x1,N_x2);
elseif nargin<4
valid = ones(N_x1,N_x2);
end
%Set up Q indexing
q1range = 1:N_k1; q2range = 1:N_k2;
[Q2,Q1]=meshgrid(q2range,q1range);
%Prealocate Q maps
spotMaps = struct('id',{spotList.id},'Q1map',zeros(N_x1,N_x2),'Q2map',zeros(N_x1,N_x2));
fprintf('Starting CBEDFFT spot-tracking map.\n'); tic
totalspots = length(spotMaps)*sum(valid(:)); spotsfit = 0; % For tracking progress
%iterate through spatial positions
for j=1:N_x1
for k=1:N_x2
if valid(j,k)
%select local CBED
CBED = data4d(:,:,j,k);
EWPC = ewpc(CBED);
%define continuous Fourier transform
PeakFun = @(x) -abs(cft2(win.*log(CBED-minval+0.1),x(1),x(2),1)); %x is [q1,q2]
%iterate through spots of interest
for s = 1:length(spotList)
%Get spot locations from input struct
spot_ROI_q1 = spotList(s).spotRangeQ1;
spot_ROI_q2 = spotList(s).spotRangeQ2;
spotNeighborhood = EWPC(spot_ROI_q1,spot_ROI_q2);
%Find rough maximum of peak
[~,maxidx] = max(spotNeighborhood(:));
Q1_roi = Q1(spot_ROI_q1,spot_ROI_q2);
Q2_roi = Q2(spot_ROI_q1,spot_ROI_q2);
Q1max = Q1_roi(maxidx); Q2max = Q2_roi(maxidx);
%Search for spot peak in continuous Fourier transform
constrainedPeakFun = @(x) ConstrainedFun(x,PeakFun,...
[spot_ROI_q1(1),spot_ROI_q1(end)],[spot_ROI_q2(1),spot_ROI_q2(end)]);
[peakQ] = fminsearch(constrainedPeakFun,[Q1max,Q2max],optimset('TolX',tol,'TolFun',tol));
%Assign in maps
spotMaps(s).Q1map(j,k) = peakQ(1);
spotMaps(s).Q2map(j,k) = peakQ(2);
spotsfit = spotsfit+1;
end
else
for s = 1:length(spotList)
spotMaps(s).Q1map(j,k) = nan;
spotMaps(s).Q2map(j,k) = nan;
end
end
end
fprintf('Completed %.1f percent of map. About %d seconds remain.\n',...
100*spotsfit/totalspots,round(toc/spotsfit*(totalspots-spotsfit)));
end
t=toc;
fprintf('\nDone. Total time: %.1f s. Time per peak fit: %.3f s.\n',t,t/totalspots)
% Add vector-form maps to spotMaps struct
spotMaps = calculateSpotMapVectors( spotMaps);
end
%%
function y = ConstrainedFun(x,fun,win1,win2)
%adds constraint to objective function fun, which is assumed to be always
%negative, by adding a positive "cone of shame" outside the specified
%window
if x(1) < win1(1) || x(1) > win1(2) || x(2) < win2(1) || x(2) > win2(2)
cent = [mean(win1),mean(win2)];
y = norm(x-cent);
else
y = fun(x);
end
end
%%
function w = window2(N,M,w_func)
%Makes a 2D window function for FFT
wc=window(w_func,N);
wr=window(w_func,M);
[maskr,maskc]=meshgrid(wr,wc);
w=maskr.*maskc;
end
%%
function [ spotMaps_updated ] = calculateSpotMapVectors( spotMaps)
%calculateSpotMapVectors Calculates vector components, length, and angle
% assuming detector has 124 pixels, i.e. zero=(63,63). These are added to
% the spotMaps struct
numSpots = length(spotMaps);
spotMaps_updated=spotMaps;
for i=1:numSpots
x1map=spotMaps(i).Q1map;
x1map=x1map-63; x1map(x1map>62) = x1map(x1map>62)- 124;
x2map=spotMaps(i).Q2map;
x2map=x2map-63; x2map(x2map>62) = x2map(x2map>62)- 124;
spotlength=sqrt(x1map.^2+x2map.^2);
spotangle=atan2d(x1map,x2map);
spotMaps_updated(i).VectorX1 = x1map;
spotMaps_updated(i).VectorX2 = x2map;
spotMaps_updated(i).VectorLength = spotlength;
spotMaps_updated(i).VectorAngle = spotangle;
end
end