3
function [validx, validy]=pickpeak;
7
[Image,PathImage] = uigetfile('*.tif','Open Image');
10
filenumber=size(filenamelist);
11
filenumber=filenumber(1);
15
I=imread(Image); % read Image
16
Itemp=mean(double(I),3);
18
[Isizey, Isizex]=size(I);
20
prompt = {'How many peaks do you want to follow?'};
21
dlg_title = 'Manual peak picking';
24
answer = inputdlg(prompt,dlg_title,num_lines,def);
25
numberofpeaks= str2num(cell2mat(answer(1,1)));
27
figure, imshow(uint8(I)); %show Image
30
title(sprintf('Mark the region of interest: Click on the on the lower left corner and and then on the upper right corner'));
34
for i=1:numberofpeaks;
36
[xprof, yprof]=ginput(2); % Get the Area of Interest
38
cropxy(i,2) = (round(xprof(2,1)-xprof(1,1))/2)+xprof(1,1);
39
cropxy(i,3) = (round(yprof(1,1)-yprof(2,1))/2)+yprof(2,1);
40
cropxy(i,4) = xprof(1,1);
41
cropxy(i,5) = yprof(2,1);
42
cropxy(i,6) = round((xprof(2,1)-xprof(1,1))/2);
43
cropxy(i,7) = round((yprof(1,1)-yprof(2,1))/2);
49
plot(cropxy(i,2),cropxy(i,3),'o');
54
tic; % start timer for time estimation
55
% I2 = imsubtract (I, imopen(I,strel('disk',15))); % subtract background
57
% [Isizey, Isizex]=size(I2);
58
% image(I2); %show with subtracted background5
66
% Start fitting process of the peaks, labeled by bwlabel
70
g = waitbar(0,'Processing image'); % nucleating the progress bar
71
fitcountertemp=size(cropxy); % number off peaks to cycle through
72
fitcounter=fitcountertemp(1,1); % number off peaks to cycle through
73
for c=1:fitcounter %start the loop to process all points
74
waitbar(c/(fitcounter-1)); % growth of the progress bar
75
cropI=imcrop(I2,[round(cropxy(c,4)) round(cropxy(c,5)) round(cropxy(c,6))*2 round(cropxy(c,7))*2]); % crop the region around the detected peak (bwlabel)
77
% get line profile in x direction for the fitting routine
79
xdata = [(round(cropxy(c,2))-round(cropxy(c,6))):1:(round(cropxy(c,2))+round(cropxy(c,6)))]; % x-coordinate for the fitting which is equivalent to the x coordinate in the image
80
ydata=sum(cropI)/(2*cropxy(c,7)); % y-coordinate for the fitting which is equivalent to the greyvalues in the image - integrated in y direction of the image
82
% fitting in x-direction
83
% guess some parameters for the fitting routine --> bad guesses lead to
84
% an error message which stops the fitting
86
back_guess=(ydata(1)+ydata(round(cropxy(c,6))*2))/2; % guess for the background level - average of the first and last greyvalue
87
sig1_guess=(cropxy(c,6)*2)/5; % guess for the peak width - take a fith of the cropping width
88
amp_guess1=ydata(round((length(ydata))/2))-back_guess; % guess for the amplitude - take the greyvalue at the peak position
89
mu1_guess=cropxy(c,2); % guess for the position of the peak - take the position from bwlabel
91
% start fitting routine
92
[x,resnormx,residual,exitflagx,output] = lsqcurvefit(@gauss_onepk, [amp_guess1 mu1_guess sig1_guess back_guess], xdata, ydata); %give the initial guesses and data to the gauss fitting routine
94
% show the fitting results
95
xtest = [(round(cropxy(c,2))-round(cropxy(c,6))):0.1:(round(cropxy(c,2))+round(cropxy(c,6)))]; % x values for the plot of the fitting result
96
ytest = (x(1)*exp((-(xtest-x(2)).^2)./(2.*x(3).^2))) + x(4); % y values of the fitting result
97
yguess=(amp_guess1*exp((-(xtest-mu1_guess).^2)./(2.*sig1_guess.^2))) + back_guess; %y values for the guess plot
98
plot(xdata,ydata,'o') % plot the experimental data
100
plot(xtest,ytest,'r') % plot the fitted function
101
plot(xtest,yguess,'b') % plot the guessed function
105
% fitting in y-direction
106
% guess parameters for the fitting routine --> bad guesses lead to
107
% an error message which stops the fitting
109
% get line profile in x direction for the fitting routine
111
xdata = [(round(cropxy(c,3))-round(cropxy(c,7))):1:(round(cropxy(c,3))+round(cropxy(c,7)))]; % x data in y direction of the image
112
ydata=sum(cropI')/(2*cropxy(c,6)); % integrate greyvalues in x direction and normalize it to the number of integrated lines
114
% fitting in y-direction
115
% guess parameters for the fitting routine --> bad guesses lead to
116
% an error message which stops the fitting
118
back_guess=(ydata(1)+ydata(round(cropxy(c,7))*2))/2; % guess for the background level - average of the first and last greyvalue
119
sig1_guess=(cropxy(c,6)*2)/5; % guess for the peak width - take a fith of the cropping width
120
amp_guess1=ydata(round((length(ydata))/2))-back_guess; % guess for the amplitude - take the greyvalue at the peak position
121
mu1_guess=cropxy(c,3); % guess for the position of the peak - take the position from bwlabel
123
% start fitting routine
124
[y,resnormy,residual,exitflagy,output] = lsqcurvefit(@gauss_onepk, [amp_guess1 mu1_guess sig1_guess back_guess], xdata, ydata); %give the initial guesses and data to the gauss fitting routine
126
% show the fitting results
127
xtest = [(round(cropxy(c,3))-round(cropxy(c,7))):0.1:(round(cropxy(c,3))+round(cropxy(c,7)))]; % x values for the plot of the fitting result
128
ytest = (y(1)*exp((-(xtest-y(2)).^2)./(2.*y(3).^2))) + y(4); % y values of the fitting result
129
yguess=(amp_guess1*exp((-(xtest-mu1_guess).^2)./(2.*sig1_guess.^2))) + back_guess; %y values for the guess plot
130
plot(xdata,ydata,'o') % plot the experimental data
132
plot(xtest,ytest,'g') % plot the fitted function
133
plot(xtest,yguess,'b') % plot the guessed function
137
% sort out the bad points and save the good ones in fitxy
138
% this matrix contains the to be used points from the first image
140
if exitflagx>0 % if the fitting routine didn't find end before the 4000th iteration (check that in lsqcurvefit.m) then exitflag will be equal or smaller then 0
142
if exitflagy>0 % the same for the y direction fitting
144
if x(3)>0.05 % the width of the peak should be wider than 1 pixel - this is negotiable: different powder particle or cameras can give back results with very narrow peaks
146
if y(3)>0.05 % the same for y direction fitting
148
if resnormx<5000 % A measure of the "goodness" of fit is the residual, the difference between the observed and predicted data. (in the help file: Mathematics: Data Analysis and Statistics: Analyzing Residuals)
150
if resnormy<5000 % the same for the y- direction - - - a good value is as far as I know until now between 30 and 50. The good fits stay well beyond that (between 0 and 10)
152
if (round(x(2))-round(cropxy(c,6)))>0
154
if (round(x(2))+round(cropxy(c,6)))<Isizex
156
if (round(y(2))-round(cropxy(c,7)))>0
158
if (round(y(2))+round(cropxy(c,7)))<Isizey
161
fitxy(counter,1)=c; % points get their final number
162
fitxy(counter,2)=x(1); % fitted amplitude x-direction
163
fitxy(counter,3)=abs(x(2)); % fitted position of the peak x-direction
164
fitxy(counter,4)=abs(x(3)); % fitted peak width in x-direction
165
fitxy(counter,5)=(x(4)); % fitted background in x-direction
166
fitxy(counter,6)=y(1); % fitted amplitude y-direction
167
fitxy(counter,7)=abs(y(2)); % fitted position of the peak y-direction
168
fitxy(counter,8)=abs(y(3)); % fitted peak width in y-direction
169
fitxy(counter,9)=abs(y(4)); % fitted background in y-direction
170
fitxy(counter,10)=cropxy(c,6); % cropping width in x-direction
171
fitxy(counter,11)=cropxy(c,7); % cropping width in ydirection
186
close(g) % close progress bar window
187
t(1,3)=toc; % stop timer
188
image_time_s=t(1,3); % take time per image
189
estimated_totaltime_h=image_time_s*filenumber/3600; % calculate estimated time
194
% plot image with peaks labeled by bwlabel (crosses) and the chosen points
195
% which are easy to fit with a gaussian distribution (circles)
197
figure, image(I2); %show Image
198
title(['Number of selected Images: ', num2str(filenumber), '; Estimated time [h] ', num2str((round(estimated_totaltime_h*10)/10)), ' Crosses are determined peaks, circles are chosen for the analysis. If you want to run the analysis hit ENTER'])
201
plot(cropxy(:,2),cropxy(:,3),'+','Color','white') % peaks from bwlabel
202
plot(fitxy(:,3),fitxy(:,7),'o','Color','white'); % "good" points
205
total_progress=1/filenumber;
210
fitlength=size(fitxy);
211
fitcounter=fitlength(1,1)
212
% again for all images
213
for m=1:(filenumber-1) % loop through all images
216
f = waitbar(0,'Working on Image');
217
I = imread(filenamelist(m,:)); %read image
218
Itemp=mean(double(I),3);
223
for c=1:fitcounter %loop trough all points
224
waitbar(c/(fitcounter-1)); %progress bar
227
pointnumber=fitxy(c,(m-1)*12+1);
228
amp_guess_x=fitxy(c,(m-1)*12+2);
229
mu_guess_x=fitxy(c,(m-1)*12+3);
230
sig_guess_x=fitxy(c,(m-1)*12+4);
231
back_guess_x=fitxy(c,(m-1)*12+5);
232
amp_guess_y=fitxy(c,(m-1)*12+6);
233
mu_guess_y=fitxy(c,(m-1)*12+7);
234
sig_guess_y=fitxy(c,(m-1)*12+8);
235
back_guess_y=fitxy(c,(m-1)*12+9);
236
crop_x=fitxy(c,(m-1)*12+10);
237
crop_y=fitxy(c,(m-1)*12+11);
239
% crop the area around the point to fit
241
cropedI=imcrop(I,[(round(mu_guess_x)-round(crop_x)) (round(mu_guess_y)-round(crop_y)) 2*round(crop_x) 2*round(crop_y)]);
242
% cropI=imsubtract (cropedI, imopen(cropedI,strel('disk',15))); % subtract background
243
cropI=cropedI;% imshow(cropI)
244
% get line profile in x direction
245
xdatax = [(round(mu_guess_x)-round(crop_x)):1:(round(mu_guess_x)+round(crop_x))];
246
ydatax=sum(cropI)/(2*(crop_y));
247
xguessx = [(round(mu_guess_x)-round(crop_x)):0.1:(round(mu_guess_x)+round(crop_x))];
248
yguessx = (amp_guess_x*exp((-(xguessx-mu_guess_x).^2)./(2.*sig_guess_x.^2))) + back_guess_x;
249
[x,resnormx,residualx,exitflagx,output] = lsqcurvefit(@gauss_onepk, [amp_guess_x mu_guess_x sig_guess_x back_guess_x], xdatax, ydatax);
250
xtestx = [(round(mu_guess_x)-round(crop_x)):0.1:(round(mu_guess_x)+round(crop_x))];
251
ytestx = (x(1)*exp((-(xtestx-x(2)).^2)./(2.*x(3).^2))) + x(4);
252
% plot(xdatax,ydatax,'o')
254
% plot(xtestx,ytestx,'r')
255
% plot(xguessx,yguessx,'b')
256
% title(['Filename: ',filenamelist(m,:), '; Progress [%]: ',num2str((round(total_progress*10))/10), '; Tot. t [h] ', num2str((round(total_time_h*10)/10)), '; Est. t [h] ', num2str((round(estimated_totaltime_h*10)/10))])
260
xdatay = [(round(mu_guess_y)-round(crop_y)):1:(round(mu_guess_y)+round(crop_y))];
261
ydatay=sum(cropI')/(2*(crop_y));
262
xguessy = [(round(mu_guess_y)-round(crop_y)):0.1:(round(mu_guess_y)+round(crop_y))];
263
yguessy = (amp_guess_y*exp((-(xguessy-mu_guess_y).^2)./(2.*sig_guess_y.^2))) + back_guess_y;
264
[y,resnormy,residualy,exitflagy,output] = lsqcurvefit(@gauss_onepk, [amp_guess_y mu_guess_y sig_guess_y back_guess_y], xdatay, ydatay);
265
xtesty = [(round(mu_guess_y)-round(crop_y)):0.1:(round(mu_guess_y)+round(crop_y))];
266
ytesty= (y(1)*exp((-(xtesty-y(2)).^2)./(2.*y(3).^2))) + y(4);
267
% plot(xdatay,ydatay,'o')
269
% plot(xtesty,ytesty,'g')
270
% plot(xguessy,yguessy,'b')
271
% title(['Filename: ',filenamelist(m,:), '; Progress [%]: ',num2str((round(total_progress*10))/10), '; Tot. t [h] ', num2str((round(total_time_h*10)/10)), '; Est. t [h] ', num2str((round(estimated_totaltime_h*10)/10))])
279
fitxy(counter,m*12+1)=pointnumber;
280
fitxy(counter,m*12+2)=abs(x(1));
281
fitxy(counter,m*12+3)=abs(x(2));
282
fitxy(counter,m*12+4)=abs(x(3));
283
fitxy(counter,m*12+5)=abs(x(4));
284
fitxy(counter,m*12+6)=abs(y(1));
285
fitxy(counter,m*12+7)=abs(y(2));
286
fitxy(counter,m*12+8)=abs(y(3));
287
fitxy(counter,m*12+9)=abs(y(4));
288
fitxy(counter,m*12+10)=crop_x;
289
fitxy(counter,m*12+11)=crop_y;
290
fitxy(counter,m*12+12)=resnormx;
297
imshow(uint8((I))); %show Image
298
title(['Filename: ',filenamelist(m,:), '; Progress [%]: ',num2str((round(total_progress*10))/10), '; Tot. t [h] ', num2str((round(total_time_h*10)/10)), '; Est. t [h] ', num2str((round(estimated_totaltime_h*10)/10))])
301
plot(fitxy(:,m*12+3),fitxy(:,m*12+7),'o','Color','white'); % "good" points
304
total_progress=1/filenumber;
309
% plot(fitxy(:,m*12+1),fitxy(:,m*12+12),'+');
310
% title(['Filename: ',filenamelist(m,:), '; Progress [%]: ',num2str((round(total_progress*10))/10), '; Tot. t [h] ', num2str((round(total_time_h*10)/10)), '; Est. t [h] ', num2str((round(estimated_totaltime_h*10)/10))])
314
total_time_s=sum(time);
315
total_time_h=sum(time)/3600;
316
image_time_s=total_time_s/m;
317
estimated_totaltime_h=image_time_s*(filenumber)/3600;
318
progress_percent=total_time_h/estimated_totaltime_h*100;
319
total_progress=(m+1)/(filenumber)*100;
324
save fitxy.dat fitxy -ascii -tabs
326
[validx,validy]=sortvalidpoints(fitxy);
327
title(['Processing Images finished!'])
332
save validx.dat validx -ascii -tabs;
333
save validy.dat validy -ascii -tabs;