3
function [validx,validy]=peak_labelling;
5
% The peak labelling function is an alternative methode to the
6
% automate_image function. The difference between the two functions is,
7
% that the automate_image function the correlation coefficient uses to
8
% track a fixed grid of markers, while the peak_labelling function is
9
% searching for peaks in a base image and tryes to fit a gauss function to
10
% it. Therefore the image should have a very low background and bright
11
% round peaks. If this is the case peak_labelling will find these peaks fit
12
% the gauss function in x and y direction to each of the peaks and then
13
% tracks all peaks in all images. The output files are fitxy.dat,
14
% validx.dat and validy.dat, which will end up in the Current directory of
15
% matlab. Attention with each run these files will be overwritten.
17
% The peak_labelling function is a bit less sensible to image noise, is
18
% only tracking markers in the image which are actually there and can under
19
% certain circumstances provide a higher accuracy for larger markers.
22
[Image,PathImage] = uigetfile('*.tif','Open Image');
25
filenumber=length(filenamelist);
28
I=imread(Image); % read Image
29
Itemp=mean(double(I),3);
31
figure, image(I); %show Image
34
title(sprintf('Mark the region of interest: Click on the on the lower left corner and and then on the upper right corner'))
35
[xprof, yprof]=ginput(2); % Get the Area of Interest
40
tic; % start timer for time estimation
41
msgboxwicon=msgbox('Subtracting image background, please wait.','Processing...')
42
I2 = imsubtract (I, imopen(I,strel('disk',15))); % subtract background
44
image(I2); %show with subtracted background
49
roi=(I2>10); %subtract greyvalues to work only with real peaks
50
[labeled,numObjects] = bwlabel(roi,8); %label all peaks - very important function, crucial for the whole process; see matlab manual
51
powderdata=regionprops(labeled,'basic') % get peak properties from bwlabel
52
powderarea=[powderdata.Area]; %define area variable
53
powdercentroid=[powderdata.Centroid]; %define position variable
54
powderboundingbox=[powderdata.BoundingBox]; %define bounding box variable
56
countermax=length(powdercentroid)/2;
57
powderxy=zeros(countermax,8);
58
for i=1:countermax; % get all data from the bwlabel (position, bounding box and area of peaks)
60
powderxy(i,1)=i; % number of the detected particle
61
powderxy(i,2)=powdercentroid(1, (i*2-1)); % x coordinate of particle position
62
powderxy(i,3)=powdercentroid(1, (i*2)); % y coordinate of particle position
63
powderxy(i,4)=powderboundingbox(1, (i*4)-3); % x coordinate of bounding box
64
powderxy(i,5)=powderboundingbox(1, (i*4)-2); % y coordinate of bounding box
65
powderxy(i,6)=powderboundingbox(1, (i*4)-1); % width (x) of bounding box
66
powderxy(i,7)=powderboundingbox(1, (i*4)); % height (y) of bounding box
67
powderxy(i,8)=powderarea(1, i); % area of bounding box
71
% cropping in x direction to reduce to the area of interest
72
% crop away peaks which are too small or too big
74
Amin=10; %minimum peaksize 10 pixel --> only testwise
75
Amax=1000; %maximum peaksize --> only testwise
79
% throw away useless peaks (defined by position and size)
80
for i=1:countermax; % loop through all points
82
if xmin<powderxy(i,2) % crop all points left from Region Of Interest (ROI)
83
if powderxy(i,2)<xmax % crop all points right from Region Of Interest (ROI)
84
if ymin<powderxy(i,3) % crop all points below the Region Of Interest (ROI)
85
if powderxy(i,3)<ymax % crop all points above the Region Of Interest (ROI)
86
if Amin<powderxy(i,8) % crop all points with a small peak area
87
if powderxy(i,8)<Amax % crop all points with a too big area
89
cropxy(counter,1)=counter; % peaks get a new number
90
cropxy(counter,2)=powderxy(i,2); % x
91
cropxy(counter,3)=powderxy(i,3); % y
92
cropxy(counter,4)=powderxy(i,4); % x bounding box
93
cropxy(counter,5)=powderxy(i,5); % y bounding box
94
cropxy(counter,6)=powderxy(i,6); % width (x) bounding box
95
cropxy(counter,7)=powderxy(i,7); % height (y) bounding box
96
cropxy(counter,8)=powderxy(i,8); % area bounding box
110
clear powderboundingbox
119
t(1,2)=toc; % stop timer
122
% Start fitting process of the peaks, labeled by bwlabel
126
g = waitbar(0,'Processing image'); % nucleating the progress bar
127
fitcountertemp=size(cropxy); % number off peaks to cycle through
128
fitcounter=fitcountertemp(1,1); % number off peaks to cycle through
129
for c=1:fitcounter %start the loop to process all points
130
waitbar(c/(fitcounter-1)); % growth of the progress bar
131
cropI=imcrop(I,[(round(cropxy(c,2))-round(cropxy(c,6))) (round(cropxy(c,3))-round(cropxy(c,7))) round(cropxy(c,6))*2 round(cropxy(c,7))*2]); % crop the region around the detected peak (bwlabel)
133
% get line profile in x direction for the fitting routine
135
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
136
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
138
% fitting in x-direction
139
% guess some parameters for the fitting routine --> bad guesses lead to
140
% an error message which stops the fitting
142
back_guess=(ydata(1)+ydata(round(cropxy(c,6))*2))/2; % guess for the background level - averadge of the first and last greyvalue
143
sig1_guess=(cropxy(c,6))/5; % guess for the peak width - take a fith of the cropping width
144
amp_guess1=ydata(round(cropxy(c,6))); % guess for the amplitude - take the greyvalue at the peak position
145
mu1_guess=cropxy(c,2); % guess for the position of the peak - take the position from bwlabel
147
% start fitting routine
148
[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
150
% show the fitting results
151
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
152
ytest = (x(1)*exp((-(xtest-x(2)).^2)./(2.*x(3).^2))) + x(4); % y values of the fitting result
153
yguess=(amp_guess1*exp((-(xtest-mu1_guess).^2)./(2.*sig1_guess.^2))) + back_guess; %y values for the guess plot
154
% plot(xdata,ydata,'o') % plot the experimental data
156
% plot(xtest,ytest,'r') % plot the fitted function
157
% plot(xtest,yguess,'b') % plot the guessed function
161
% fitting in y-direction
162
% guess parameters for the fitting routine --> bad guesses lead to
163
% an error message which stops the fitting
165
% get line profile in x direction for the fitting routine
167
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
168
ydata=sum(cropI')/(2*cropxy(c,6)); % integrate greyvalues in x direction and normalize it to the number of integrated lines
170
% fitting in y-direction
171
% guess parameters for the fitting routine --> bad guesses lead to
172
% an error message which stops the fitting
174
back_guess=(ydata(1)+ydata(round(cropxy(c,7))*2))/2; % guess for the background level - averadge of the first and last greyvalue
175
sig1_guess=(cropxy(c,6))/5; % guess for the peak width - take a fith of the cropping width
176
amp_guess1=ydata(round(cropxy(c,7))); % guess for the amplitude - take the greyvalue at the peak position
177
mu1_guess=cropxy(c,3); % guess for the position of the peak - take the position from bwlabel
179
% start fitting routine
180
[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
182
% show the fitting results
183
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
184
ytest = (y(1)*exp((-(xtest-y(2)).^2)./(2.*y(3).^2))) + y(4); % y values of the fitting result
185
xguess = [(round(cropxy(c,3))-round(cropxy(c,7))):0.1:(round(cropxy(c,3))+round(cropxy(c,7)))]; % x values for the guess
186
yguess = (amp_guess1*exp((-(xguess-mu1_guess.^2)./(2.*sig1_guess.^2))) + back_guess); % y values for the guess plpot
187
% plot(xdata,ydata,'o') % plot the experimental data
189
% plot(xtest,ytest,'g') % plot the fitted function
190
% plot(xguess,yguess,'b') % plot the guessed function
194
% sort out the bad points and save the good ones in fitxy
195
% this matrix contains the to be used points from the first image
197
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
198
if exitflagy>0 % the same for the y direction fitting
199
if x(3)>1 % 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
200
if y(3)>1 % the same for y direction fitting
201
if resnormx/cropxy(c,6)<10 % A measure of the "happyness" of a fit is the residual, the difference between the observed and predicted data. (in the help file: Mathematics: Data Analysis and Statistics: Analyzing Residuals)
202
if resnormy/cropxy(c,7)<10 % 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)
204
fitxy(counter,1)=c; % points get their final number
205
fitxy(counter,2)=abs(x(1)); % fitted amplitude x-direction
206
fitxy(counter,3)=abs(x(2)); % fitted position of the peak x-direction
207
fitxy(counter,4)=abs(x(3)); % fitted peak width in x-direction
208
fitxy(counter,5)=(x(4)); % fitted background in x-direction
209
fitxy(counter,6)=abs(y(1)); % fitted amplitude y-direction
210
fitxy(counter,7)=abs(y(2)); % fitted position of the peak y-direction
211
fitxy(counter,8)=abs(y(3)); % fitted peak width in y-direction
212
fitxy(counter,9)=abs(y(4)); % fitted background in y-direction
213
fitxy(counter,10)=cropxy(c,6); % cropping width in x-direction
214
fitxy(counter,11)=cropxy(c,7); % cropping width in ydirection
225
close(g) % close progress bar window
226
t(1,3)=toc; % stop timer
227
image_time_s=t(1,3); % take time per image
228
estimated_totaltime_h=image_time_s*filenumber/3600 % calculate estimated time
233
% plot image with peaks labeled by bwlabel (crosses) and the chosen points
234
% which are easy to fit with a gaussian distribution (circles)
236
figure, image(I2); %show Image
237
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'])
240
plot(cropxy(:,2),cropxy(:,3),'+','Color','white') % peaks from bwlabel
241
plot(fitxy(:,3),fitxy(:,7),'o','Color','white'); % "good" points
244
total_progress=1/filenumber;
249
fitlength=size(fitxy);
250
fitcounter=fitlength(1,1)
251
% again for all images
252
for m=1:(filenumber-1) % loop through all images
255
I = imread(filenamelist(m,:)); %read image
256
Itemp=mean(double(I),3);
258
f = waitbar(0,'Working on Image');
262
for c=1:fitcounter %loop trough all points
263
waitbar(c/(fitcounter-1)); %progress bar
266
pointnumber=fitxy(c,(m-1)*12+1);
267
amp_guess_x=fitxy(c,(m-1)*12+2);
268
mu_guess_x=fitxy(c,(m-1)*12+3);
269
sig_guess_x=fitxy(c,(m-1)*12+4);
270
back_guess_x=fitxy(c,(m-1)*12+5);
271
amp_guess_y=fitxy(c,(m-1)*12+6);
272
mu_guess_y=fitxy(c,(m-1)*12+7);
273
sig_guess_y=fitxy(c,(m-1)*12+8);
274
back_guess_y=fitxy(c,(m-1)*12+9);
275
crop_x=fitxy(c,(m-1)*12+10);
276
crop_y=fitxy(c,(m-1)*12+11);
278
% crop the area around the point to fit
280
cropI=imcrop(I,[(round(mu_guess_x)-round(crop_x)) (round(mu_guess_y)-round(crop_y)) 2*round(crop_x) 2*round(crop_y)]);
281
% 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)]);
282
% cropI=imsubtract (cropedI, imopen(cropedI,strel('disk',15))); % subtract background
284
% get line profile in x direction
285
xdatax = [(round(mu_guess_x)-round(crop_x)):1:(round(mu_guess_x)+round(crop_x))];
286
ydatax=sum(cropI)/(2*(crop_y));
287
xguessx = [(round(mu_guess_x)-round(crop_x)):0.1:(round(mu_guess_x)+round(crop_x))];
288
yguessx = (amp_guess_x*exp((-(xguessx-mu_guess_x).^2)./(2.*sig_guess_x.^2))) + back_guess_x;
289
[x,resnormx,residualx,exitflagx,output] = lsqcurvefit(@gauss_onepk, [amp_guess_x mu_guess_x sig_guess_x back_guess_x], xdatax, ydatax);
290
xtestx = [(round(mu_guess_x)-round(crop_x)):0.1:(round(mu_guess_x)+round(crop_x))];
291
ytestx = (x(1)*exp((-(xtestx-x(2)).^2)./(2.*x(3).^2))) + x(4);
292
% plot(xdatax,ydatax,'o')
294
% plot(xtestx,ytestx,'r')
295
% plot(xguessx,yguessx,'b')
296
% 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))])
299
xdatay = [(round(mu_guess_y)-round(crop_y)):1:(round(mu_guess_y)+round(crop_y))];
300
ydatay=sum(cropI')/(2*(crop_y));
301
xguessy = [(round(mu_guess_y)-round(crop_y)):0.1:(round(mu_guess_y)+round(crop_y))];
302
yguessy = (amp_guess_y*exp((-(xguessy-mu_guess_y).^2)./(2.*sig_guess_y.^2))) + back_guess_y;
303
[y,resnormy,residualy,exitflagy,output] = lsqcurvefit(@gauss_onepk, [amp_guess_y mu_guess_y sig_guess_y back_guess_y], xdatay, ydatay);
304
xtesty = [(round(mu_guess_y)-round(crop_y)):0.1:(round(mu_guess_y)+round(crop_y))];
305
ytesty= (y(1)*exp((-(xtesty-y(2)).^2)./(2.*y(3).^2))) + y(4);
306
% plot(xdatay,ydatay,'o')
308
% plot(xtesty,ytesty,'g')
309
% plot(xguessy,yguessy,'b')
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))])
318
fitxy(counter,m*12+1)=pointnumber;
319
fitxy(counter,m*12+2)=abs(x(1));
320
fitxy(counter,m*12+3)=abs(x(2));
321
fitxy(counter,m*12+4)=abs(x(3));
322
fitxy(counter,m*12+5)=abs(x(4));
323
fitxy(counter,m*12+6)=abs(y(1));
324
fitxy(counter,m*12+7)=abs(y(2));
325
fitxy(counter,m*12+8)=abs(y(3));
326
fitxy(counter,m*12+9)=abs(y(4));
327
fitxy(counter,m*12+10)=crop_x;
328
fitxy(counter,m*12+11)=crop_y;
329
fitxy(counter,m*12+12)=resnormx;
337
plot(fitxy(:,m*12+1),fitxy(:,m*12+12),'+');
338
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))])
342
total_time_s=sum(time);
343
total_time_h=sum(time)/3600;
344
image_time_s=total_time_s/m;
345
estimated_totaltime_h=image_time_s*(filenumber)/3600
346
progress_percent=total_time_h/estimated_totaltime_h*100;
347
total_progress=(m+1)/(filenumber)*100
353
save fitxy.dat fitxy -ascii -tabs
355
[validx,validy]=sortvalidpoints(fitxy);
356
title(['Processing Images finished!'])
361
save validx.dat validx -ascii -tabs;
362
save validy.dat validy -ascii -tabs;
b'\\ No newline at end of file'