bzr branch
http://suren.me/webbzr/normxcorr/trunk
41
by Suren A. Chilingaryan
Synchronize with Chris version, set precision to 1, initial implementation of labview wrapper |
1 |
%% written by Chris
|
2 |
||
3 |
% updated last 9/8/2010
|
|
4 |
% this function allows to track multiple peaks along one axis. The
|
|
5 |
% resolution can be much higher than the cross correlation methods of
|
|
6 |
% digital image correlation.
|
|
7 |
||
8 |
% function goes here
|
|
9 |
function [validx,validy]=multipeak_tracking; |
|
10 |
||
11 |
%% Load image and ask for orientation and pos/neg
|
|
12 |
||
13 |
if exist('filenamelist')==0 |
|
14 |
load('filenamelist') % file with the list of filenames to be processed |
|
15 |
end
|
|
16 |
||
17 |
filenumber=size(filenamelist); |
|
18 |
filenumber=filenumber(1); |
|
19 |
||
20 |
h=figure; |
|
21 |
imshow(filenamelist(1,:)) % show the first image |
|
22 |
title('First Image') % put a title |
|
23 |
||
24 |
ImageA=(uint8(mean(double(imread(filenamelist(1,:))),3))); |
|
25 |
||
26 |
sel_pos_neg = menu(sprintf('Do you want use the negative image for analysis?'),... |
|
27 |
'As is', 'Negativ', 'Cancel') |
|
28 |
if sel_pos_neg==3 |
|
29 |
close all; |
|
30 |
return
|
|
31 |
end
|
|
32 |
if sel_pos_neg==2 |
|
33 |
ImageA=255-ImageA; |
|
34 |
imshow(ImageA) % show the first image |
|
35 |
negativ=1; |
|
36 |
end
|
|
37 |
if sel_pos_neg==1 |
|
38 |
negativ=0; |
|
39 |
end
|
|
40 |
||
41 |
sel_horz_vert = menu(sprintf('Do you want to measure horizontal or vertical?'),... |
|
42 |
'Horizontal', 'Vertical', 'Cancel') |
|
43 |
if sel_horz_vert==3 |
|
44 |
close all; |
|
45 |
return
|
|
46 |
end
|
|
47 |
if sel_horz_vert==2 |
|
48 |
ImageA=ImageA'; |
|
49 |
imshow(ImageA) % show the first image |
|
50 |
orientation=90; |
|
51 |
end
|
|
52 |
if sel_horz_vert==1 |
|
53 |
orientation=0; |
|
54 |
end
|
|
55 |
||
56 |
[imr,imc]=size(ImageA); |
|
57 |
||
58 |
||
59 |
%% select path for displacement extraction along markers
|
|
60 |
%Select point in middle of marker
|
|
61 |
xlabel('Location On Image [Pixels]') |
|
62 |
ylabel('Location On Image [Pixels]') |
|
63 |
title(['Click on the center of the sample. Width: ', num2str(imc),... |
|
64 |
'; Height: ',num2str(imr)]) |
|
65 |
marker_pt=round(ginput(1)); |
|
66 |
x_mark = marker_pt(1); |
|
67 |
y_mark = marker_pt(2); |
|
68 |
hold on |
|
69 |
line([1 imc],[y_mark y_mark],'Color','r') |
|
70 |
||
71 |
||
72 |
%Prompt for integration width
|
|
73 |
prompt = {'Enter integration width [Pixels]:'}; |
|
74 |
dlg_title = 'Input integration width for the analysis'; |
|
75 |
num_lines= 1; |
|
76 |
def = {'20'}; |
|
77 |
answer = inputdlg(prompt,dlg_title,num_lines,def); |
|
78 |
int_width = str2num(cell2mat(answer(1))); |
|
79 |
||
80 |
line([1 imc],[y_mark-int_width/2 y_mark-int_width/2],'Color','g') |
|
81 |
line([1 imc],[y_mark+int_width/2 y_mark+int_width/2],'Color','g') |
|
82 |
||
83 |
||
84 |
%% Calculate line profile data and select peaks
|
|
85 |
xdata = [1:1:imc]; |
|
86 |
ydata= sum(ImageA((y_mark-int_width/2):(y_mark+int_width/2),:),1)/int_width; |
|
87 |
||
88 |
g=figure; |
|
89 |
plot(xdata,ydata) |
|
90 |
xlabel('Location On Image [Pixels]') |
|
91 |
ylabel('Pixel Intensity Value') |
|
92 |
||
93 |
sel_peakselect = menu(sprintf('Do you want to select single or all peaks?')... |
|
94 |
,'Single Select', 'All Select', 'Cancel') |
|
95 |
if sel_peakselect==3 |
|
96 |
close all; |
|
97 |
return
|
|
98 |
end
|
|
99 |
||
100 |
% If single peaks should be chosen
|
|
101 |
if sel_peakselect==1 |
|
102 |
sel_peakselect_2=1 |
|
103 |
i=1 |
|
104 |
hold on |
|
105 |
while sel_peakselect_2==1 |
|
106 |
title('Click on the left of the chosen peak.') |
|
107 |
marker_pt=round(ginput(1)); |
|
108 |
x_peak_left = marker_pt(1); |
|
109 |
y_peak_left = marker_pt(2); |
|
110 |
line([x_peak_left x_peak_left], [1 max(ydata)],'Color','r') |
|
111 |
||
112 |
title('Click on the right of the chosen peak.') |
|
113 |
marker_pt=round(ginput(1)); |
|
114 |
x_peak_right = marker_pt(1); |
|
115 |
y_peak_right = marker_pt(2); |
|
116 |
line([x_peak_right x_peak_right], [1 max(ydata)],'Color','r') |
|
117 |
||
118 |
plot(xdata(x_peak_left:x_peak_right),ydata(x_peak_left:x_peak_right),'k') |
|
119 |
peak_pos(i,:)=[x_peak_left x_peak_right]; |
|
120 |
||
121 |
i=i+1; |
|
122 |
sel_peakselect_2 = menu(sprintf('Do you want to select another peak?')... |
|
123 |
,'One more', 'No, lets process') |
|
124 |
end
|
|
125 |
end
|
|
126 |
||
127 |
% If all peaks should be found automatically
|
|
128 |
if sel_peakselect==2 |
|
129 |
||
130 |
end
|
|
131 |
||
132 |
%% Start first fitting |
|
133 |
fitcountertemp=size(peak_pos); % number off peaks to cycle through |
|
134 |
fitcounter=fitcountertemp(1,1); % number off peaks to cycle through |
|
135 |
||
136 |
for c=1:fitcounter %start the loop to process all points |
|
137 |
fitdatax=xdata(peak_pos(c,1):peak_pos(c,2)); |
|
138 |
fitdatay=ydata(peak_pos(c,1):peak_pos(c,2)); |
|
139 |
||
140 |
% fitting in x-direction
|
|
141 |
% guess some parameters for the fitting routine --> bad guesses lead to
|
|
142 |
% an error message which stops the fitting
|
|
143 |
||
144 |
back_guess=(ydata(peak_pos(c,2))+ydata(peak_pos(c,2)))/2; % guess for the background level - average of the first and last greyvalue |
|
145 |
sig1_guess=(peak_pos(c,2)-peak_pos(c,1))/2; % guess for the peak width - take half of the cropping width |
|
146 |
amp_guess1=max(fitdatay); % guess for the amplitude - take the greyvalue at the peak position |
|
147 |
mu1_guess=(peak_pos(c,2)+peak_pos(c,1))/2; % guess for the position of the peak - take the position from bwlabel |
|
148 |
||
149 |
% start fitting routine
|
|
150 |
[x,resnormx,residual,exitflagx,output] = lsqcurvefit(@gauss_onepk, [amp_guess1 mu1_guess sig1_guess back_guess], fitdatax, fitdatay); %give the initial guesses and data to the gauss fitting routine |
|
151 |
||
152 |
% show the fitting results
|
|
153 |
xtest = fitdatax; % x values for the plot of the fitting result |
|
154 |
ytest = (x(1)*exp((-(xtest-x(2)).^2)./(2.*x(3).^2))) + x(4); % y values of the fitting result |
|
155 |
yguess=(amp_guess1*exp((-(xtest-mu1_guess).^2)./(2.*sig1_guess.^2))) + back_guess; %y values for the guess plot |
|
156 |
plot(xtest,ytest,'r') % plot the fitted function |
|
157 |
plot(xtest,yguess,'b') % plot the guessed function |
|
158 |
drawnow
|
|
159 |
||
160 |
peaks(c,1,:)=x; |
|
161 |
end
|
|
162 |
||
163 |
sel_go = menu(sprintf('Do you want proceed?')... |
|
164 |
,'Do it!', 'No, stop!') |
|
165 |
||
166 |
if sel_peakselect==2 |
|
167 |
close all; |
|
168 |
return
|
|
169 |
end
|
|
170 |
||
171 |
cutout=peak_pos(:,2)-peak_pos(:,1) |
|
172 |
||
173 |
for m=1:(filenumber-1) % loop through all images |
|
174 |
ImageA=((mean(double(imread(filenamelist(m,:))),3))); |
|
175 |
||
176 |
if sel_horz_vert==2 |
|
177 |
ImageA=ImageA'; |
|
178 |
imshow(ImageA) % show the first image |
|
179 |
orientation=90; |
|
180 |
end
|
|
181 |
if sel_pos_neg==2 |
|
182 |
ImageA=255-ImageA; |
|
183 |
imshow(ImageA) % show the first image |
|
184 |
negativ=1; |
|
185 |
end
|
|
186 |
||
187 |
image=m
|
|
188 |
||
189 |
xdata = [1:1:imc]; |
|
190 |
ydata= sum(ImageA((y_mark-int_width/2):(y_mark+int_width/2),:),1)/int_width; |
|
191 |
||
192 |
plot(xdata,ydata) |
|
193 |
title(['Image Number: ', num2str(m)]) |
|
194 |
hold on |
|
195 |
||
196 |
||
197 |
for c=1:fitcounter %start the loop to process all points |
|
198 |
||
199 |
if (peaks(c,m,2)+cutout(c)/2)>imc |
|
200 |
fitdatax=xdata((imc-cutout(c)):imc); |
|
201 |
fitdatay=ydata((imc-cutout(c)):imc); |
|
202 |
else
|
|
203 |
fitdatax=xdata((peaks(c,m,2)-cutout(c)/2):(peaks(c,m,2)+cutout(c)/2)); |
|
204 |
fitdatay=ydata((peaks(c,m,2)-cutout(c)/2):(peaks(c,m,2)+cutout(c)/2)); |
|
205 |
end
|
|
206 |
||
207 |
% fitting in x-direction |
|
208 |
% guess some parameters for the fitting routine --> bad guesses lead to
|
|
209 |
% an error message which stops the fitting
|
|
210 |
||
211 |
back_guess=peaks(c,m,4); % guess for the background level - average of the first and last greyvalue |
|
212 |
sig1_guess=peaks(c,m,3); % guess for the peak width - take half of the cropping width |
|
213 |
amp_guess1=peaks(c,m,1); % guess for the amplitude - take the greyvalue at the peak position |
|
214 |
mu1_guess=peaks(c,m,2); % guess for the position of the peak - take the position from bwlabel |
|
215 |
||
216 |
% start fitting routine
|
|
217 |
[x,resnormx,residual,exitflagx,output] = lsqcurvefit(@gauss_onepk, [amp_guess1 mu1_guess sig1_guess back_guess], fitdatax, fitdatay); %give the initial guesses and data to the gauss fitting routine |
|
218 |
% show the fitting results
|
|
219 |
xtest = fitdatax; % x values for the plot of the fitting result |
|
220 |
ytest = (x(1)*exp((-(xtest-x(2)).^2)./(2.*x(3).^2))) + x(4); % y values of the fitting result |
|
221 |
yguess=(amp_guess1*exp((-(xtest-mu1_guess).^2)./(2.*sig1_guess.^2))) + back_guess; %y values for the guess plot |
|
222 |
plot(xtest,ytest,'r') % plot the fitted function |
|
223 |
plot(xtest,yguess,'g') % plot the guessed function |
|
224 |
drawnow
|
|
225 |
peaks(c,m+1,:)=x; |
|
226 |
end
|
|
227 |
hold off |
|
228 |
Amplitude=peaks(:,m,1)'; |
|
229 |
Peakposition=peaks(:,m,2)'; |
|
230 |
Peakwidth=peaks(:,m,3)'; |
|
231 |
Background=peaks(:,m,4)'; |
|
232 |
dlmwrite('backupPeakposition.txt', Peakposition' , 'delimiter', '\t', '-append'); % Here we save the result from each image; if you are desperately want to run this function with e.g. matlab 6.5 then you should comment this line out. If you do that the data will be saved at the end of the correlation step - good luck ;-) |
|
233 |
dlmwrite('backupAmplitude.txt', Amplitude' , 'delimiter', '\t', '-append'); % Here we save the result from each image; if you are desperately want to run this function with e.g. matlab 6.5 then you should comment this line out. If you do that the data will be saved at the end of the correlation step - good luck ;-) |
|
234 |
dlmwrite('backupPeakwidth.txt', Peakwidth' , 'delimiter', '\t', '-append'); % Here we save the result from each image; if you are desperately want to run this function with e.g. matlab 6.5 then you should comment this line out. If you do that the data will be saved at the end of the correlation step - good luck ;-) |
|
235 |
dlmwrite('backupBackground.txt', Background' , 'delimiter', '\t', '-append'); % Here we save the result from each image; if you are desperately want to run this function with e.g. matlab 6.5 then you should comment this line out. If you do that the data will be saved at the end of the correlation step - good luck ;-) |
|
236 |
end
|
|
237 |
||
238 |
%% Save data
|
|
239 |
validx=peaks(:,:,2); |
|
240 |
validy=peaks(:,:,3); |
|
241 |
||
242 |
Amplitude=peaks(:,:,1)'; |
|
243 |
Peakposition=peaks(:,:,2)'; |
|
244 |
Peakwidth=peaks(:,:,3)'; |
|
245 |
Background=peaks(:,:,4)'; |
|
246 |
||
247 |
save peaks; |
|
248 |
save Amplitude.dat Amplitude -ascii -tabs; |
|
249 |
save Peakposition.dat Peakposition -ascii -tabs; |
|
250 |
save Peakwidth.dat Peakwidth -ascii -tabs; |
|
251 |
save Background.dat Background -ascii -tabs; |
|
252 |
save validx_mp.dat validx -ascii -tabs; |
|
253 |
save validy_mp.dat validy -ascii -tabs; |