bzr branch
http://suren.me/webbzr/normxcorr/trunk
1
by Suren A. Chilingaryan
Initial import |
1 |
% written by Chris
|
2 |
||
3 |
function [validx,validy]=peak_labelling; |
|
4 |
||
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.
|
|
16 |
%
|
|
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.
|
|
20 |
||
21 |
||
22 |
[Image,PathImage] = uigetfile('*.tif','Open Image'); |
|
23 |
cd(PathImage) |
|
24 |
load('filenamelist') |
|
25 |
filenumber=length(filenamelist); |
|
26 |
% filelist_generator
|
|
27 |
||
28 |
I=imread(Image); % read Image |
|
29 |
Itemp=mean(double(I),3); |
|
30 |
I=Itemp; |
|
31 |
figure, image(I); %show Image |
|
32 |
axis('equal'); |
|
33 |
drawnow
|
|
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 |
|
36 |
xmin = xprof(1,1); |
|
37 |
xmax = xprof(2,1); |
|
38 |
ymin = yprof(2,1); |
|
39 |
ymax= yprof(1,1); |
|
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 |
|
43 |
close(msgboxwicon) |
|
44 |
image(I2); %show with subtracted background |
|
45 |
axis('equal'); |
|
46 |
t(1,1)=toc; |
|
47 |
tic; |
|
48 |
drawnow
|
|
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 |
|
55 |
counter=0; |
|
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) |
|
59 |
counter=counter+1; % |
|
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 |
|
68 |
end
|
|
69 |
||
70 |
||
71 |
% cropping in x direction to reduce to the area of interest
|
|
72 |
% crop away peaks which are too small or too big
|
|
73 |
||
74 |
Amin=10; %minimum peaksize 10 pixel --> only testwise |
|
75 |
Amax=1000; %maximum peaksize --> only testwise |
|
76 |
counter=0 |
|
77 |
i=0; |
|
78 |
||
79 |
% throw away useless peaks (defined by position and size)
|
|
80 |
for i=1:countermax; % loop through all points |
|
81 |
||
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 |
|
88 |
counter=counter+1; |
|
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 |
|
97 |
end
|
|
98 |
end
|
|
99 |
end
|
|
100 |
end
|
|
101 |
end
|
|
102 |
end
|
|
103 |
end
|
|
104 |
||
105 |
% clear variables
|
|
106 |
clear powderxy |
|
107 |
clear powderdata |
|
108 |
clear powderarea |
|
109 |
clear powdercentroid |
|
110 |
clear powderboundingbox |
|
111 |
clear counter |
|
112 |
clear countermax |
|
113 |
clear Amin |
|
114 |
clear Amax |
|
115 |
clear roi |
|
116 |
||
117 |
close all |
|
118 |
||
119 |
t(1,2)=toc; % stop timer |
|
120 |
||
121 |
||
122 |
% Start fitting process of the peaks, labeled by bwlabel
|
|
123 |
||
124 |
tic; % start timer |
|
125 |
counter=0 |
|
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) |
|
132 |
||
133 |
% get line profile in x direction for the fitting routine
|
|
134 |
||
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 |
|
137 |
||
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
|
|
141 |
||
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 |
|
146 |
||
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 |
|
149 |
||
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
|
|
155 |
% hold on
|
|
156 |
% plot(xtest,ytest,'r') % plot the fitted function
|
|
157 |
% plot(xtest,yguess,'b') % plot the guessed function
|
|
158 |
% drawnow
|
|
159 |
% hold off
|
|
160 |
||
161 |
% fitting in y-direction
|
|
162 |
% guess parameters for the fitting routine --> bad guesses lead to
|
|
163 |
% an error message which stops the fitting
|
|
164 |
||
165 |
% get line profile in x direction for the fitting routine
|
|
166 |
||
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 |
|
169 |
||
170 |
% fitting in y-direction
|
|
171 |
% guess parameters for the fitting routine --> bad guesses lead to
|
|
172 |
% an error message which stops the fitting
|
|
173 |
||
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 |
|
178 |
||
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 |
|
181 |
||
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
|
|
188 |
% hold on
|
|
189 |
% plot(xtest,ytest,'g') % plot the fitted function
|
|
190 |
% plot(xguess,yguess,'b') % plot the guessed function
|
|
191 |
% drawnow
|
|
192 |
% hold off
|
|
193 |
||
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
|
|
196 |
||
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)
|
|
203 |
counter=counter+1;
|
|
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
|
|
215 |
end
|
|
216 |
end
|
|
217 |
end
|
|
218 |
end
|
|
219 |
end
|
|
220 |
end
|
|
221 |
|
|
222 |
end
|
|
223 |
||
224 |
||
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
|
|
229 |
sum(t);
|
|
230 |
total_time_h=sum(t)
|
|
231 |
close all
|
|
232 |
||
233 |
% plot image with peaks labeled by bwlabel (crosses) and the chosen points
|
|
234 |
% which are easy to fit with a gaussian distribution (circles)
|
|
235 |
||
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']) |
|
238 |
axis('equal');
|
|
239 |
hold on;
|
|
240 |
plot(cropxy(:,2),cropxy(:,3),'+','Color','white') % peaks from bwlabel |
|
241 |
plot(fitxy(:,3),fitxy(:,7),'o','Color','white'); % "good" points |
|
242 |
drawnow
|
|
243 |
||
244 |
total_progress=1/filenumber; |
|
245 |
||
246 |
pause
|
|
247 |
||
248 |
close all |
|
249 |
fitlength=size(fitxy); |
|
250 |
fitcounter=fitlength(1,1) |
|
251 |
% again for all images
|
|
252 |
for m=1:(filenumber-1) % loop through all images |
|
253 |
tic; %start timer |
|
254 |
counter=0; |
|
255 |
I = imread(filenamelist(m,:)); %read image |
|
256 |
Itemp=mean(double(I),3); |
|
257 |
I=Itemp; |
|
258 |
f = waitbar(0,'Working on Image'); |
|
259 |
||
260 |
||
261 |
% loop number
|
|
262 |
for c=1:fitcounter %loop trough all points |
|
263 |
waitbar(c/(fitcounter-1)); %progress bar |
|
264 |
||
265 |
% load variables
|
|
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); |
|
277 |
||
278 |
% crop the area around the point to fit
|
|
279 |
||
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
|
|
283 |
% imshow(cropI)
|
|
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')
|
|
293 |
% hold on
|
|
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))])
|
|
297 |
% drawnow
|
|
298 |
% hold off
|
|
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')
|
|
307 |
% hold on
|
|
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))])
|
|
311 |
% drawnow
|
|
312 |
% hold off
|
|
313 |
||
314 |
||
315 |
if exitflagx>0 |
|
316 |
if exitflagy>0 |
|
317 |
counter=counter+1; |
|
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; |
|
330 |
||
331 |
end
|
|
332 |
end
|
|
333 |
||
334 |
||
335 |
end
|
|
336 |
||
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))]) |
|
339 |
fitcounter=counter; |
|
340 |
drawnow
|
|
341 |
time(m)=toc; |
|
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 |
|
348 |
close(f); |
|
349 |
||
350 |
end
|
|
351 |
||
352 |
% save the stuff |
|
353 |
save fitxy.dat fitxy -ascii -tabs |
|
354 |
||
355 |
[validx,validy]=sortvalidpoints(fitxy); |
|
356 |
title(['Processing Images finished!']) |
|
357 |
||
358 |
save('validx'); |
|
359 |
save('validy'); |
|
360 |
||
361 |
save validx.dat validx -ascii -tabs; |
|
362 |
save validy.dat validy -ascii -tabs; |