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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
|
%% written by Chris
% updated last 9/8/2010
% this function allows to track multiple peaks along one axis. The
% resolution can be much higher than the cross correlation methods of
% digital image correlation.
% function goes here
function [validx,validy]=multipeak_tracking;
%% Load image and ask for orientation and pos/neg
if exist('filenamelist')==0
load('filenamelist') % file with the list of filenames to be processed
end
filenumber=size(filenamelist);
filenumber=filenumber(1);
h=figure;
imshow(filenamelist(1,:)) % show the first image
title('First Image') % put a title
ImageA=(uint8(mean(double(imread(filenamelist(1,:))),3)));
sel_pos_neg = menu(sprintf('Do you want use the negative image for analysis?'),...
'As is', 'Negativ', 'Cancel')
if sel_pos_neg==3
close all;
return
end
if sel_pos_neg==2
ImageA=255-ImageA;
imshow(ImageA) % show the first image
negativ=1;
end
if sel_pos_neg==1
negativ=0;
end
sel_horz_vert = menu(sprintf('Do you want to measure horizontal or vertical?'),...
'Horizontal', 'Vertical', 'Cancel')
if sel_horz_vert==3
close all;
return
end
if sel_horz_vert==2
ImageA=ImageA';
imshow(ImageA) % show the first image
orientation=90;
end
if sel_horz_vert==1
orientation=0;
end
[imr,imc]=size(ImageA);
%% select path for displacement extraction along markers
%Select point in middle of marker
xlabel('Location On Image [Pixels]')
ylabel('Location On Image [Pixels]')
title(['Click on the center of the sample. Width: ', num2str(imc),...
'; Height: ',num2str(imr)])
marker_pt=round(ginput(1));
x_mark = marker_pt(1);
y_mark = marker_pt(2);
hold on
line([1 imc],[y_mark y_mark],'Color','r')
%Prompt for integration width
prompt = {'Enter integration width [Pixels]:'};
dlg_title = 'Input integration width for the analysis';
num_lines= 1;
def = {'20'};
answer = inputdlg(prompt,dlg_title,num_lines,def);
int_width = str2num(cell2mat(answer(1)));
line([1 imc],[y_mark-int_width/2 y_mark-int_width/2],'Color','g')
line([1 imc],[y_mark+int_width/2 y_mark+int_width/2],'Color','g')
%% Calculate line profile data and select peaks
xdata = [1:1:imc];
ydata= sum(ImageA((y_mark-int_width/2):(y_mark+int_width/2),:),1)/int_width;
g=figure;
plot(xdata,ydata)
xlabel('Location On Image [Pixels]')
ylabel('Pixel Intensity Value')
sel_peakselect = menu(sprintf('Do you want to select single or all peaks?')...
,'Single Select', 'All Select', 'Cancel')
if sel_peakselect==3
close all;
return
end
% If single peaks should be chosen
if sel_peakselect==1
sel_peakselect_2=1
i=1
hold on
while sel_peakselect_2==1
title('Click on the left of the chosen peak.')
marker_pt=round(ginput(1));
x_peak_left = marker_pt(1);
y_peak_left = marker_pt(2);
line([x_peak_left x_peak_left], [1 max(ydata)],'Color','r')
title('Click on the right of the chosen peak.')
marker_pt=round(ginput(1));
x_peak_right = marker_pt(1);
y_peak_right = marker_pt(2);
line([x_peak_right x_peak_right], [1 max(ydata)],'Color','r')
plot(xdata(x_peak_left:x_peak_right),ydata(x_peak_left:x_peak_right),'k')
peak_pos(i,:)=[x_peak_left x_peak_right];
i=i+1;
sel_peakselect_2 = menu(sprintf('Do you want to select another peak?')...
,'One more', 'No, lets process')
end
end
% If all peaks should be found automatically
if sel_peakselect==2
end
%% Start first fitting
fitcountertemp=size(peak_pos); % number off peaks to cycle through
fitcounter=fitcountertemp(1,1); % number off peaks to cycle through
for c=1:fitcounter %start the loop to process all points
fitdatax=xdata(peak_pos(c,1):peak_pos(c,2));
fitdatay=ydata(peak_pos(c,1):peak_pos(c,2));
% fitting in x-direction
% guess some parameters for the fitting routine --> bad guesses lead to
% an error message which stops the fitting
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
sig1_guess=(peak_pos(c,2)-peak_pos(c,1))/2; % guess for the peak width - take half of the cropping width
amp_guess1=max(fitdatay); % guess for the amplitude - take the greyvalue at the peak position
mu1_guess=(peak_pos(c,2)+peak_pos(c,1))/2; % guess for the position of the peak - take the position from bwlabel
% start fitting routine
[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
% show the fitting results
xtest = fitdatax; % x values for the plot of the fitting result
ytest = (x(1)*exp((-(xtest-x(2)).^2)./(2.*x(3).^2))) + x(4); % y values of the fitting result
yguess=(amp_guess1*exp((-(xtest-mu1_guess).^2)./(2.*sig1_guess.^2))) + back_guess; %y values for the guess plot
plot(xtest,ytest,'r') % plot the fitted function
plot(xtest,yguess,'b') % plot the guessed function
drawnow
peaks(c,1,:)=x;
end
sel_go = menu(sprintf('Do you want proceed?')...
,'Do it!', 'No, stop!')
if sel_peakselect==2
close all;
return
end
cutout=peak_pos(:,2)-peak_pos(:,1)
for m=1:(filenumber-1) % loop through all images
ImageA=((mean(double(imread(filenamelist(m,:))),3)));
if sel_horz_vert==2
ImageA=ImageA';
imshow(ImageA) % show the first image
orientation=90;
end
if sel_pos_neg==2
ImageA=255-ImageA;
imshow(ImageA) % show the first image
negativ=1;
end
image=m
xdata = [1:1:imc];
ydata= sum(ImageA((y_mark-int_width/2):(y_mark+int_width/2),:),1)/int_width;
plot(xdata,ydata)
title(['Image Number: ', num2str(m)])
hold on
for c=1:fitcounter %start the loop to process all points
if (peaks(c,m,2)+cutout(c)/2)>imc
fitdatax=xdata((imc-cutout(c)):imc);
fitdatay=ydata((imc-cutout(c)):imc);
else
fitdatax=xdata((peaks(c,m,2)-cutout(c)/2):(peaks(c,m,2)+cutout(c)/2));
fitdatay=ydata((peaks(c,m,2)-cutout(c)/2):(peaks(c,m,2)+cutout(c)/2));
end
% fitting in x-direction
% guess some parameters for the fitting routine --> bad guesses lead to
% an error message which stops the fitting
back_guess=peaks(c,m,4); % guess for the background level - average of the first and last greyvalue
sig1_guess=peaks(c,m,3); % guess for the peak width - take half of the cropping width
amp_guess1=peaks(c,m,1); % guess for the amplitude - take the greyvalue at the peak position
mu1_guess=peaks(c,m,2); % guess for the position of the peak - take the position from bwlabel
% start fitting routine
[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
% show the fitting results
xtest = fitdatax; % x values for the plot of the fitting result
ytest = (x(1)*exp((-(xtest-x(2)).^2)./(2.*x(3).^2))) + x(4); % y values of the fitting result
yguess=(amp_guess1*exp((-(xtest-mu1_guess).^2)./(2.*sig1_guess.^2))) + back_guess; %y values for the guess plot
plot(xtest,ytest,'r') % plot the fitted function
plot(xtest,yguess,'g') % plot the guessed function
drawnow
peaks(c,m+1,:)=x;
end
hold off
Amplitude=peaks(:,m,1)';
Peakposition=peaks(:,m,2)';
Peakwidth=peaks(:,m,3)';
Background=peaks(:,m,4)';
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 ;-)
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 ;-)
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 ;-)
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 ;-)
end
%% Save data
validx=peaks(:,:,2);
validy=peaks(:,:,3);
Amplitude=peaks(:,:,1)';
Peakposition=peaks(:,:,2)';
Peakwidth=peaks(:,:,3)';
Background=peaks(:,:,4)';
save peaks;
save Amplitude.dat Amplitude -ascii -tabs;
save Peakposition.dat Peakposition -ascii -tabs;
save Peakwidth.dat Peakwidth -ascii -tabs;
save Background.dat Background -ascii -tabs;
save validx_mp.dat validx -ascii -tabs;
save validy_mp.dat validy -ascii -tabs;
|