bzr branch
http://suren.me/webbzr/normxcorr/trunk
1
by Suren A. Chilingaryan
Initial import |
1 |
function [validx,validy]=automate_image(grid_x,grid_y,filenamelist,validx,validy); |
2 |
||
3 |
% Code to start actual image correlation
|
|
4 |
% Programmed by Chris and Rob
|
|
5 |
% Last revision: 09/10/08
|
|
6 |
||
7 |
% The automation function is the central function and processes all markers and
|
|
8 |
% images by the use of the matlab function cpcorr.m.
|
|
9 |
% Therefore the Current directory in matlab has to be the folder where
|
|
10 |
% automate_image.m finds the filenamelist.mat, grid_x.dat and grid_y.dat as well
|
|
11 |
% as the images specified in filenamelist.mat. Just type automate_image; and
|
|
12 |
% press ENTER at the command line of matlab.
|
|
13 |
% At first, automate_image.m will open the first image in the filenamelist.mat and
|
|
14 |
% plot the grid as green crosses on top. The next step will need some time since
|
|
15 |
% all markers in that image have to be processed for the first image. After correlating
|
|
16 |
% image one and two the new raster positions will be plotted as red crosses. On top
|
|
17 |
% of the image and the green crosses. The next dialog will ask you if you want to
|
|
18 |
% continue with this correlation or cancel. If you press continue, automate_image.m
|
|
19 |
% will process all images in the filenamelist.mat. The time it will take to process
|
|
20 |
% all images will be plotted on the figure but can easily be estimated by knowing the
|
|
21 |
% raster point processing speed (see processing speed).
|
|
22 |
% Depending on the number of images and markers you are tracking, this process
|
|
23 |
% can take between seconds and days. For 100 images and 200 markers a decent
|
|
24 |
% computer should need 200 seconds. To get a better resolution you can always
|
|
25 |
% run jobs overnight (e.g. 6000 markers in 1000 images) with higher resolutions.
|
|
26 |
% Keep in mind that CORRSIZE which you changed in cpcorr.m will limit your
|
|
27 |
% resolution. If you chose to use the 15 pixel as suggested a marker distance of
|
|
28 |
% 30 pixel will lead to a full cover of the strain field. Choosing smaller marker
|
|
29 |
% distances will lead to an interpolation since two neighboring markers share
|
|
30 |
% pixels. Nevertheless a higher marker density can reduce the noise of the strain field.
|
|
31 |
% When all images are processed, automate_image will write the files validx.mat,
|
|
32 |
% validy.mat, validx.txt and validy.txt. The text files are meant to store the result in a
|
|
33 |
% format which can be accessed by other programs also in the future.
|
|
34 |
||
2
by Suren A. Chilingaryan
Support for different optimization modes |
35 |
% Changed by Suren A. Chilingaryan <csa@dside.dyndns.org> to optimize performance
|
36 |
% by tighter integration with Matlab images toolkit (few sources from the toolkit
|
|
37 |
% are moved to the current source tree and adjusted to benefit from knowledge of
|
|
38 |
% tasks we are solving here. As well CUDA framework, if available, is used to
|
|
39 |
% speed-up core computations using parallel processors of latest NVidia video
|
|
40 |
% cards.
|
|
41 |
% OPTIMIZE parameter is controlling which optimizations should be used.
|
|
42 |
% 0 - The original version, no optimizations
|
|
43 |
% 1 - The modified sources from images toolkit are used
|
|
44 |
% 2 - The CUDA matlab plugin is used to compute FFT's
|
|
45 |
% 3 - Most of computations are performed on NVidia card
|
|
16
by Suren A. Chilingaryan
Optimize FFT size |
46 |
% 4 - Optimize FFT sizes for better performance (affects precision)
|
24
by Suren A. Chilingaryan
Optimization mode 5 for matlab: image loading in C code |
47 |
% 5 - Load images in C code
|
2
by Suren A. Chilingaryan
Support for different optimization modes |
48 |
% You can safely set an optimization level to 3. If NVidia card is not available
|
49 |
% the level will be automatically reduced to 1.
|
|
6
by Suren A. Chilingaryan
A little more computations are moved to CUDA |
50 |
%
|
51 |
% SAVEDATA controls if data is saved per iteration or at the end
|
|
52 |
% 0 - in the end
|
|
53 |
% 1 - per iteration
|
|
2
by Suren A. Chilingaryan
Support for different optimization modes |
54 |
|
29
by Suren A. Chilingaryan
Implementation of image and fragment modes, support for non-cacheable grids |
55 |
OPTIMIZE = 4; |
2
by Suren A. Chilingaryan
Support for different optimization modes |
56 |
CORRSIZE = 15; |
1
by Suren A. Chilingaryan
Initial import |
57 |
PRECISION = 1000; |
4
by Suren A. Chilingaryan
Instead of transfer compute local sums and denormals on board |
58 |
SILENT = 1; |
6
by Suren A. Chilingaryan
A little more computations are moved to CUDA |
59 |
SAVEDATA = 0; |
2
by Suren A. Chilingaryan
Support for different optimization modes |
60 |
|
31
by Suren A. Chilingaryan
CUDAfication of real-time module |
61 |
warning off Images:initSize:adjustingMag |
62 |
||
1
by Suren A. Chilingaryan
Initial import |
63 |
% exist('grid_x')
|
64 |
% exist('grid_y')
|
|
65 |
% exist('filenamelist')
|
|
66 |
% exist('validx')
|
|
67 |
% exist('validy')
|
|
68 |
||
69 |
if exist('images', 'dir') |
|
70 |
imagedir = 'images/'; |
|
71 |
else
|
|
72 |
imagedir = ''; |
|
73 |
end
|
|
74 |
||
31
by Suren A. Chilingaryan
CUDAfication of real-time module |
75 |
if exist('data', 'dir') |
76 |
datadir = 'data/'; |
|
77 |
else
|
|
78 |
if exist([imagedir, 'grid_x.dat'], 'file') |
|
79 |
datadir = imagedir; |
|
80 |
else
|
|
81 |
datadir = ''; |
|
82 |
end
|
|
83 |
end
|
|
84 |
||
85 |
||
1
by Suren A. Chilingaryan
Initial import |
86 |
% Load necessary files
|
87 |
if exist('grid_x')==0 |
|
88 |
load([datadir, 'grid_x.dat']) % file with x position, created by grid_generator.m |
|
89 |
end
|
|
90 |
if exist('grid_y')==0 |
|
91 |
load([datadir, 'grid_y.dat']) % file with y position, created by grid_generator.m |
|
92 |
end
|
|
93 |
if exist('filenamelist')==0 |
|
94 |
load([datadir, 'filenamelist']) % file with the list of filenames to be processed |
|
95 |
end
|
|
96 |
||
97 |
||
98 |
resume=0; |
|
99 |
if exist('validx')==1 |
|
100 |
if exist('validy')==1 |
|
101 |
resume=1; |
|
102 |
[Rasternum Imagenum]=size(validx); |
|
103 |
end
|
|
104 |
end
|
|
105 |
||
106 |
[r,c]=size(filenamelist); % this will determine the number of images we have to loop through |
|
107 |
||
108 |
% Open new figure so previous ones (if open) are not overwritten
|
|
109 |
if (~SILENT) |
|
24
by Suren A. Chilingaryan
Optimization mode 5 for matlab: image loading in C code |
110 |
if OPTIMIZE > 4 |
111 |
OPTIMIZE = 4 |
|
112 |
end
|
|
113 |
||
1
by Suren A. Chilingaryan
Initial import |
114 |
h=figure; |
115 |
imshow([imagedir, filenamelist(1,:)]) % show the first image |
|
116 |
title('Initial Grid For Image Correlation (Note green crosses)') % put a title |
|
117 |
hold on |
|
118 |
plot(grid_x,grid_y,'g+') % plot the grid onto the image |
|
119 |
hold off |
|
120 |
||
121 |
% Start image correlation using cpcorr.m
|
|
122 |
g = waitbar(0,sprintf('Processing images')); % initialize the waitbar |
|
123 |
set(g,'Position',[275,50,275,50]) % set the position of the waitbar [left bottom width height] |
|
124 |
end
|
|
125 |
||
126 |
firstimage=1; |
|
127 |
||
128 |
if resume==1 |
|
129 |
firstimage=Imagenum+1 |
|
130 |
end
|
|
131 |
||
132 |
%Initializing hardware code
|
|
133 |
||
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
134 |
ncp = prod(size(grid_x)); |
1
by Suren A. Chilingaryan
Initial import |
135 |
|
2
by Suren A. Chilingaryan
Support for different optimization modes |
136 |
if OPTIMIZE > 2 |
137 |
hwid = normxcorr_hw(); |
|
4
by Suren A. Chilingaryan
Instead of transfer compute local sums and denormals on board |
138 |
|
2
by Suren A. Chilingaryan
Support for different optimization modes |
139 |
if hwid > 0 |
16
by Suren A. Chilingaryan
Optimize FFT size |
140 |
err = normxcorr_hw(hwid, 1, ncp, CORRSIZE, PRECISION, OPTIMIZE); |
3
by Suren A. Chilingaryan
Add error checking to initialization process |
141 |
if (err ~= 0) |
142 |
normxcorr_hw(1); |
|
143 |
OPTIMIZE = 2; |
|
144 |
end
|
|
2
by Suren A. Chilingaryan
Support for different optimization modes |
145 |
else
|
146 |
if hwid < 0 |
|
147 |
OPTIMIZE = 1; |
|
148 |
else
|
|
149 |
OPTIMIZE = 2; |
|
150 |
end
|
|
151 |
end
|
|
152 |
else
|
|
153 |
hwid = 0; |
|
154 |
end
|
|
1
by Suren A. Chilingaryan
Initial import |
155 |
|
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
156 |
% Initialize variables |
157 |
if OPTIMIZE>2 |
|
158 |
base_points_x=single(grid_x); |
|
159 |
base_points_y=single(grid_y); |
|
160 |
else
|
|
161 |
base_points_x=grid_x; |
|
162 |
base_points_y=grid_y; |
|
163 |
end
|
|
164 |
||
165 |
if resume==1 |
|
166 |
if OPTIMIZE>2 |
|
167 |
input_points_x=single(validx(:,Imagenum)); |
|
168 |
input_points_y=single(validy(:,Imagenum)); |
|
169 |
else
|
|
170 |
input_points_x=validx(:,Imagenum); |
|
171 |
input_points_y=validy(:,Imagenum); |
|
172 |
end
|
|
173 |
inputpoints=1; |
|
174 |
else
|
|
175 |
if OPTIMIZE>2 |
|
176 |
input_points_x=single(grid_x); |
|
177 |
input_points_y=single(grid_y); |
|
178 |
else
|
|
179 |
input_points_x=grid_x; |
|
180 |
input_points_y=grid_y; |
|
181 |
end
|
|
182 |
end
|
|
183 |
||
184 |
||
185 |
if OPTIMIZE > 2 |
|
186 |
normxcorr_hw(hwid, 3, base_points_x, base_points_y); |
|
24
by Suren A. Chilingaryan
Optimization mode 5 for matlab: image loading in C code |
187 |
if OPTIMIZE > 4 |
188 |
normxcorr_hw(hwid, 4, strcat(imagedir, filenamelist(1,:))); |
|
189 |
else
|
|
190 |
base = uint8(mean(double(imread([imagedir, filenamelist(1,:)])),3)); % read in the base image ( which is always image number one. You might want to change that to improve correlation results in case the light conditions are changing during the experiment |
|
29
by Suren A. Chilingaryan
Implementation of image and fragment modes, support for non-cacheable grids |
191 |
|
192 |
% Verification of GPU local sum and denom computations
|
|
193 |
% dic_basecorr3(CORRSIZE, PRECISION, OPTIMIZE, hwid, base_points_x, base_points_y, base);
|
|
194 |
||
24
by Suren A. Chilingaryan
Optimization mode 5 for matlab: image loading in C code |
195 |
normxcorr_hw(hwid, 4, base); |
196 |
end
|
|
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
197 |
else
|
24
by Suren A. Chilingaryan
Optimization mode 5 for matlab: image loading in C code |
198 |
base = uint8(mean(double(imread([imagedir, filenamelist(1,:)])),3)); % read in the base image ( which is always image number one. You might want to change that to improve correlation results in case the light conditions are changing during the experiment |
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
199 |
base_points_for(:,1)=reshape(base_points_x,[],1); |
200 |
base_points_for(:,2)=reshape(base_points_y,[],1); |
|
201 |
||
202 |
%some checks (moved here from cpcorr)
|
|
203 |
if any(base_points_for(:)<0.5) || any(base_points_for(:,1)>size(base,2)+0.5) || any(base_points_for(:,2)>size(base,1)+0.5) |
|
204 |
msg = sprintf('In function %s, Control Points must be in pixel coordinates.',mfilename); |
|
205 |
eid = sprintf('Images:%s:cpPointsMustBeInPixCoord',mfilename); |
|
206 |
error(eid, msg); |
|
207 |
end
|
|
208 |
||
209 |
if OPTIMIZE > 0 |
|
210 |
data_base = struct; |
|
1
by Suren A. Chilingaryan
Initial import |
211 |
|
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
212 |
%moved from normxcorr2
|
213 |
% We can precalculate it here, since near-edge images are dropped and, therefore, A,T always have the same size
|
|
214 |
% n is width and height of T, but following computations are done for A
|
|
215 |
% outsize is sum if widths and heights of T and A - 1 (outsize = size(A) + size(T) - 1)
|
|
216 |
||
217 |
outsize = 6 * CORRSIZE + 1; |
|
218 |
n = 2*CORRSIZE + 1; |
|
219 |
mn = n*n; |
|
220 |
||
221 |
rects_base = dic_calc_rects(base_points_for,2*CORRSIZE,base); |
|
222 |
||
223 |
for icp = 1:ncp |
|
224 |
if isequal(rects_base(icp,3:4),[0 0]) |
|
225 |
%near edge
|
|
226 |
data_base(icp).skip = 1; |
|
227 |
continue
|
|
4
by Suren A. Chilingaryan
Instead of transfer compute local sums and denormals on board |
228 |
end
|
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
229 |
|
230 |
sub_base = imcrop(base,rects_base(icp,:)); |
|
231 |
||
232 |
data_base(icp).skip = 0; |
|
233 |
||
4
by Suren A. Chilingaryan
Instead of transfer compute local sums and denormals on board |
234 |
double_sub_base = double(sub_base); |
235 |
||
236 |
if (numel(sub_base) < 2) || (std(double_sub_base(:)) == 0) |
|
237 |
% This check is moved here for normxcorr2
|
|
238 |
eid = sprintf('Images:%s:sameElementsInTemplate',mfilename); |
|
239 |
msg = 'The values of TEMPLATE cannot all be the same.'; |
|
240 |
error(eid,'%s',msg); |
|
241 |
end
|
|
242 |
||
243 |
local_sum_A = dic_local_sum(double_sub_base,n,n); |
|
244 |
local_sum_A2 = dic_local_sum(double_sub_base.*double_sub_base,n,n); |
|
245 |
||
246 |
% Note: diff_local_sums should be nonnegative, but may have negative
|
|
247 |
% values due to round off errors. Below, we use max to ensure the
|
|
248 |
% radicand is nonnegative.
|
|
249 |
denom_A = sqrt( max(( local_sum_A2 - (local_sum_A.^2)/mn ) / (mn-1), 0) ); |
|
250 |
i_nonzero = find(denom_A~=0); |
|
251 |
||
2
by Suren A. Chilingaryan
Support for different optimization modes |
252 |
data_base(icp).denom_A = denom_A; |
253 |
data_base(icp).i_nonzero = i_nonzero; |
|
254 |
||
255 |
data_base(icp).local_sum_A = local_sum_A; |
|
256 |
data_base(icp).sub_base = sub_base; |
|
257 |
||
258 |
if (OPTIMIZE > 1) |
|
259 |
data_base(icp).sub_base_fft = fft2_cuda(double_sub_base, outsize, outsize); |
|
260 |
else
|
|
261 |
data_base(icp).sub_base_fft = fft2(double_sub_base, outsize, outsize); |
|
262 |
end
|
|
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
263 |
|
264 |
data_base(icp).base_fractional_offset = base_points_for(icp,:) - round(base_points_for(icp,:)*PRECISION)/PRECISION; |
|
2
by Suren A. Chilingaryan
Support for different optimization modes |
265 |
end
|
266 |
end
|
|
1
by Suren A. Chilingaryan
Initial import |
267 |
end
|
268 |
||
269 |
%normxcorr_hw(hwid);
|
|
270 |
%return;
|
|
271 |
||
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
272 |
if OPTIMIZE > 2 |
273 |
normxcorr_hw(hwid, 12, input_points_x, input_points_y); |
|
274 |
else
|
|
6
by Suren A. Chilingaryan
A little more computations are moved to CUDA |
275 |
input_correl(:,1)=reshape(input_points_x,[],1); % we reshape the input points to one row of values since this is the shape cpcorr will accept |
276 |
input_correl(:,2)=reshape(input_points_y,[],1); |
|
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
277 |
end
|
6
by Suren A. Chilingaryan
A little more computations are moved to CUDA |
278 |
|
1
by Suren A. Chilingaryan
Initial import |
279 |
for i=firstimage:(r-1) % run through all images |
280 |
tic % start the timer |
|
281 |
||
282 |
||
24
by Suren A. Chilingaryan
Optimization mode 5 for matlab: image loading in C code |
283 |
if OPTIMIZE > 4 |
284 |
normxcorr_hw(hwid, 13, strcat(imagedir, filenamelist((i+1),:))); |
|
285 |
input_correl = normxcorr_hw(hwid, 14); |
|
286 |
elseif OPTIMIZE > 2 |
|
29
by Suren A. Chilingaryan
Implementation of image and fragment modes, support for non-cacheable grids |
287 |
input = uint8(mean(double(imread([imagedir, filenamelist((i+1),:)])),3)); % read in the image which has to be correlated |
288 |
||
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
289 |
%Validate findpeak
|
290 |
%dic_cpcorr3(CORRSIZE, PRECISION, OPTIMIZE, hwid, ncp, input);
|
|
24
by Suren A. Chilingaryan
Optimization mode 5 for matlab: image loading in C code |
291 |
|
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
292 |
normxcorr_hw(hwid, 13, input) |
293 |
input_correl = normxcorr_hw(hwid, 14); |
|
6
by Suren A. Chilingaryan
A little more computations are moved to CUDA |
294 |
elseif OPTIMIZE > 0 |
24
by Suren A. Chilingaryan
Optimization mode 5 for matlab: image loading in C code |
295 |
input = uint8(mean(double(imread([imagedir, filenamelist((i+1),:)])),3)); % read in the image which has to be correlated |
6
by Suren A. Chilingaryan
A little more computations are moved to CUDA |
296 |
input_correl(:,:)=dic_cpcorr(CORRSIZE, PRECISION, OPTIMIZE, hwid, data_base, input_correl, input); |
2
by Suren A. Chilingaryan
Support for different optimization modes |
297 |
else
|
24
by Suren A. Chilingaryan
Optimization mode 5 for matlab: image loading in C code |
298 |
input = uint8(mean(double(imread([imagedir, filenamelist((i+1),:)])),3)); % read in the image which has to be correlated |
6
by Suren A. Chilingaryan
A little more computations are moved to CUDA |
299 |
input_correl(:,:)=cpcorr(input_correl, base_points_for, input, base); |
300 |
end
|
|
19
by Suren A. Chilingaryan
Provide stand-alone library |
301 |
|
8
by Suren A. Chilingaryan
Complete elimination of cpcorr |
302 |
%We need to copy |
6
by Suren A. Chilingaryan
A little more computations are moved to CUDA |
303 |
validx(:,i)=double(input_correl(:,1)); % the results we get from cpcorr for the x-direction |
304 |
validy(:,i)=double(input_correl(:,2)); % the results we get from cpcorr for the y-direction |
|
305 |
||
306 |
if (SAVEDATA) |
|
307 |
savelinex=validx(:,i)'; |
|
308 |
dlmwrite([datadir, 'resultsimcorrx.txt'], savelinex , '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 ;-) |
|
309 |
||
310 |
saveliney=validy(:,i)'; |
|
311 |
dlmwrite([datadir, 'resultsimcorry.txt'], saveliney , 'delimiter', '\t', '-append'); |
|
312 |
end
|
|
1
by Suren A. Chilingaryan
Initial import |
313 |
|
314 |
if (~SILENT) |
|
315 |
waitbar(i/(r-1)) % update the waitbar |
|
316 |
||
317 |
imshow([imagedir, filenamelist(i+1,:)]) % update image |
|
318 |
hold on |
|
319 |
plot(grid_x,grid_y,'g+') % plot start position of raster |
|
6
by Suren A. Chilingaryan
A little more computations are moved to CUDA |
320 |
plot(input_correl(:,1),input_correl(:,2),'r+') % plot actual postition of raster |
1
by Suren A. Chilingaryan
Initial import |
321 |
hold off |
322 |
drawnow
|
|
323 |
time(i)=toc; % take time |
|
324 |
estimatedtime=sum(time)/i*(r-1); % estimate time to process |
|
325 |
title(['# Im.: ', num2str((r-1)),'; Proc. Im. #: ', num2str((i)),'; # Rasterp.:',num2str(row*col), '; Est. Time [s] ', num2str(round(estimatedtime)), '; Elapsed Time [s] ', num2str(round(sum(time)))]); % plot a title onto the image |
|
326 |
drawnow
|
|
327 |
end
|
|
328 |
end
|
|
329 |
||
2
by Suren A. Chilingaryan
Support for different optimization modes |
330 |
if OPTIMIZE > 2 |
331 |
normxcorr_hw(hwid); |
|
332 |
end
|
|
1
by Suren A. Chilingaryan
Initial import |
333 |
|
334 |
if (~SILENT) |
|
335 |
close(g) |
|
336 |
end
|
|
337 |
||
338 |
close all |
|
339 |
||
340 |
% save
|
|
341 |
||
342 |
if (~SILENT) |
|
343 |
save ([datadir,'time.dat'], 'time', '-ascii', '-tabs'); |
|
344 |
end
|
|
345 |
||
346 |
save ([datadir,'validx.dat'], 'validx', '-ascii', '-tabs'); |
|
347 |
save ([datadir,'validy.dat'], 'validy', '-ascii', '-tabs'); |