/normxcorr/trunk

To get this branch, use:
bzr branch http://suren.me/webbzr/normxcorr/trunk
1 by Suren A. Chilingaryan
Initial import
1
function [grid_x,grid_y]=grid_generator(FileNameBase,PathNameBase);
2
3
% Code to generate the DIC analysis grid
4
% Completely rewritten by Chris
5
% Programmed first by Dan and Rob 
6
% 
7
% Last revision: 12/27/06
8
9
% The grid_generator function will help you create grids of markers. The
10
% dialog has different options allowing you to create a marker grid which is rectangular,
11
% circular, a line or two rectangels of a shape or contains only of two
12
% markers. After choosing one of the shapes you will be asked for the base
13
% image which is typically your first image. After opening that image you
14
% will be asked to click at the sites of interest and the markers will be
15
% plotted on top of your image. You can choose if you want to keep these
16
% markers or if you want to try again.
17
% It has to be noted that you can
18
% always generate your own marker positions. Therefore the marker position
19
% in pixel has to be saved as a text based format where the x-position is
20
% saved as grid_x.dat and the y-position saved as grid_y.dat.
21
%
22
23
24
25
% Prompt user for base image
26
if exist('FileNameBase')==0
27
[FileNameBase,PathNameBase] = uigetfile( ...
28
    {'*.bmp;*.tif;*.jpg;*.TIF;*.BMP;*.JPG','Image files (*.bmp,*.tif,*.jpg)';'*.*',  'All Files (*.*)'}, ...
29
    'Open base image for grid creation');
30
31
end
32
cd(PathNameBase)
33
im_grid = imread(FileNameBase);
34
35
[grid_x,grid_y,FileNameBase,PathNameBase] = gridtypeselection(FileNameBase, PathNameBase, im_grid);
36
37
close all
38
39
%-------------------------------
40
%
41
% Decide which type of grid you want to create
42
43
function [grid_x,grid_y,FileNameBase,PathNameBase] = gridtypeselection(FileNameBase, PathNameBase, im_grid);
44
45
hold off
46
imshow(im_grid,'truesize');
47
48
gridselection = menu(sprintf('Which type of grid do you want to use'),...
49
    'Rectangular','Circular','Two Markers','Line','Two Rectangles of Markers','Cancel');
50
51
52
if gridselection==1
53
    [grid_x,grid_y,FileNameBase,PathNameBase] = rect_grid(FileNameBase, PathNameBase, im_grid);
54
    return
55
end
56
57
if gridselection==2
58
    [grid_x,grid_y,FileNameBase,PathNameBase] = circ_grid(FileNameBase, PathNameBase, im_grid);
59
    return
60
end
61
62
if gridselection==3
63
    [grid_x,grid_y,FileNameBase,PathNameBase] = twop_grid(FileNameBase, PathNameBase, im_grid);
64
    return
65
end
66
67
if gridselection==4
68
    [grid_x,grid_y,FileNameBase,PathNameBase] = line_grid(FileNameBase, PathNameBase, im_grid);
69
    return
70
end
71
72
if gridselection==5
73
    [grid_x,grid_y,FileNameBase,PathNameBase] = tworect_grid(FileNameBase, PathNameBase, im_grid);
74
    return
75
end
76
77
if gridselection==6
78
    return;
79
end
80
81
82
83
%-------------------------------
84
%
85
% Define two rectangles and add them to one marker array
86
87
function [grid_x,grid_y,FileNameBase,PathNameBase] = tworect_grid(FileNameBase, PathNameBase, im_grid);
88
89
[grid_x1,grid_y1,FileNameBase,PathNameBase] = rect_grid(FileNameBase, PathNameBase, im_grid);
90
imshow(im_grid,'truesize');
91
[grid_x2,grid_y2,FileNameBase,PathNameBase] = rect_grid(FileNameBase, PathNameBase, im_grid);
92
93
grid_x1=reshape(grid_x1,[],1);
94
grid_x2=reshape(grid_x2,[],1);
95
grid_y1=reshape(grid_y1,[],1);
96
grid_y2=reshape(grid_y2,[],1);
97
98
grid_x=[grid_x1; grid_x2];
99
grid_y=[grid_y1; grid_y2];
100
101
imshow(im_grid,'truesize');
102
hold on
103
plot(grid_x,grid_y,'.')
104
title(['Selected grid has ',num2str(length(grid_x)), ' rasterpoints'])    % plot a title onto the image
105
106
% Accept the chosen markers, try again or give up 
107
108
confirmcircselection = menu(sprintf('Do you want to use these markers?'),...
109
    'Yes','No, try again','Go back to grid-type selection');
110
111
if confirmcircselection==2
112
    close all
113
    hold off
114
    imshow(im_grid,'truesize');
115
    tworect_grid(FileNameBase, PathNameBase, im_grid);
116
end
117
118
if confirmcircselection==3
119
    close all
120
    gridtypeselection(FileNameBase, PathNameBase, im_grid);
121
end
122
123
if confirmcircselection==1
124
    close all
125
    save grid_x.dat grid_x -ascii -tabs
126
    save grid_y.dat grid_y -ascii -tabs
127
end
128
129
%-------------------------------
130
%
131
% Define line and create markers
132
133
function [grid_x,grid_y,FileNameBase,PathNameBase] = line_grid(FileNameBase, PathNameBase, im_grid);
134
135
title(sprintf('Pick two points on the sample.') )
136
137
[x(1,1),y(1,1)]=ginput(1);
138
hold on
139
plot(x(1,1),y(1,1),'+g')
140
141
[x(2,1),y(2,1)]=ginput(1);
142
plot(x(2,1),y(2,1),'+g')
143
144
145
linelength=sqrt((x(2,1)-x(1,1))*(x(2,1)-x(1,1))+(y(2,1)-y(1,1))*(y(2,1)-y(1,1)));
146
lineslope=(y(2,1)-y(1,1))/(x(2,1)-x(1,1));
147
intersecty=y(1,1)-lineslope*x(1,1);
148
ycalc=zeros(2,1);
149
ycalc=lineslope*x+intersecty;
150
plot(x(:,1),ycalc(:,1),'-b')
151
152
153
prompt = {'Enter the number of intersections between markers on the line:'};
154
dlg_title = 'Input for grid creation';
155
num_lines= 1;
156
def     = {'30'};
157
answer = inputdlg(prompt,dlg_title,num_lines,def);
158
linediv = str2num(cell2mat(answer(1,1)));
159
linestep=((max(x)-min(x))/linediv);
160
grid_x(1:linediv+1)=min(x)+linestep*(1:linediv+1)-linestep;
161
grid_y=lineslope*grid_x+intersecty;
162
163
plot(grid_x,grid_y,'ob')
164
title(['Selected grid has ',num2str(linediv), ' rasterpoints'])    % plot a title onto the image
165
166
% Accept the chosen markers, try again or give up 
167
168
confirmcircselection = menu(sprintf('Do you want to use these markers?'),...
169
    'Yes','No, try again','Go back to grid-type selection');
170
171
if confirmcircselection==2
172
    close all
173
    hold off
174
    imshow(im_grid,'truesize');
175
    twop_grid(FileNameBase, PathNameBase, im_grid);
176
end
177
178
if confirmcircselection==3
179
    close all
180
    gridtypeselection(FileNameBase, PathNameBase, im_grid);
181
end
182
183
if confirmcircselection==1
184
    save grid_x.dat grid_x -ascii -tabs
185
    save grid_y.dat grid_y -ascii -tabs
186
end
187
188
%-------------------------------
189
%
190
% Select two markers
191
192
function [grid_x,grid_y,FileNameBase,PathNameBase] = twop_grid(FileNameBase, PathNameBase, im_grid);
193
194
title(sprintf('Pick two points on the sample.') )
195
196
[x(1,1),y(1,1)]=ginput(1);
197
hold on
198
plot(x(1,1),y(1,1),'+g')
199
200
[x(2,1),y(2,1)]=ginput(1);
201
plot(x(2,1),y(2,1),'+g')
202
203
% Accept the chosen markers, try again or give up 
204
205
confirmcircselection = menu(sprintf('Do you want to use these two markers?'),...
206
    'Yes','No, try again','Go back to grid-type selection');
207
208
if confirmcircselection==2
209
    close all
210
    hold off
211
    imshow(im_grid,'truesize');
212
    twop_grid(FileNameBase, PathNameBase, im_grid);
213
end
214
215
if confirmcircselection==3
216
    close all
217
    gridtypeselection(FileNameBase, PathNameBase, im_grid);
218
end
219
220
if confirmcircselection==1
221
    grid_x=x;
222
    grid_y=y;
223
    save grid_x.dat grid_x -ascii -tabs
224
    save grid_y.dat grid_y -ascii -tabs
225
end
226
%-------------------------------
227
%
228
% Select a circular area
229
230
function [grid_x,grid_y,FileNameBase,PathNameBase] = circ_grid(FileNameBase, PathNameBase, im_grid);
231
232
title(sprintf('Pick three points on the circle in clockwise order at the upper boundary of the sample.') )
233
234
[x(1,1),y(1,1)]=ginput(1);
235
hold on
236
plot(x(1,1),y(1,1),'+g')
237
238
[x(2,1),y(2,1)]=ginput(1);
239
plot(x(2,1),y(2,1),'+g')
240
241
[x(3,1),y(3,1)]=ginput(1);
242
plot(x(3,1),y(3,1),'+g')
243
244
xnew=x;
245
ynew=y;
246
247
% Calculate center between the 3 sorted points and the normal slope of the vectors
248
slope12=-1/((ynew(2,1)-ynew(1,1))/(xnew(2,1)-xnew(1,1)));
249
slope23=-1/((ynew(3,1)-ynew(2,1))/(xnew(3,1)-xnew(2,1)));
250
center12(1,1)=(xnew(2,1)-xnew(1,1))/2+xnew(1,1);
251
center12(1,2)=(ynew(2,1)-ynew(1,1))/2+ynew(1,1);
252
center23(1,1)=(xnew(3,1)-xnew(2,1))/2+xnew(2,1);
253
center23(1,2)=(ynew(3,1)-ynew(2,1))/2+ynew(2,1);
254
% plot(center12(1,1),center12(1,2),'+b')
255
% plot(center23(1,1),center23(1,2),'+b')
256
257
if slope12==slope23
258
    return
259
end
260
261
% Calculate the crossing point of the two vectors
262
achsenabschnitt1=center12(1,2)-center12(1,1)*slope12;
263
achsenabschnitt2=center23(1,2)-center23(1,1)*slope23;
264
xdata=min(x):max(x);
265
ydata1=achsenabschnitt1+slope12*xdata;
266
ydata2=achsenabschnitt2+slope23*xdata;
267
% plot(xdata,ydata1,'-b')
268
% plot(xdata,ydata2,'-b')
269
xcross=(achsenabschnitt2-achsenabschnitt1)/(slope12-slope23);
270
ycross=slope12*xcross+achsenabschnitt1;
271
plot(xcross,ycross,'or')
272
273
% Calculate radius and plot circle
274
R=sqrt((xcross-xnew(1,1))*(xcross-xnew(1,1))+(ycross-ynew(1,1))*(ycross-ynew(1,1)));
275
% ydata=ycross-sqrt(R*R-(xdata-xcross).*(xdata-xcross));
276
% plot(xdata,ydata,'-b')
277
278
% Calculate angle between vectors
279
xvector=[1;0];
280
x1vec(1,1)=xnew(1,1)-xcross;x1vec(2,1)=ynew(1,1)-ycross
281
x3vec(1,1)=xnew(3,1)-xcross;x3vec(2,1)=ynew(3,1)-ycross
282
alpha13=acos((dot(x1vec,x3vec))/(sqrt(x1vec'*x1vec)*sqrt(x3vec'*x3vec)))*180/pi;
283
alpha01=acos((dot(xvector,x1vec))/(sqrt(x1vec'*x1vec)*sqrt(xvector'*xvector)))*180/pi;
284
alpha03=acos((dot(xvector,x3vec))/(sqrt(xvector'*xvector)*sqrt(x3vec'*x3vec)))*180/pi;
285
totalangle=alpha13;
286
minangle=alpha01;
287
maxangle=alpha03;
288
angldiv=abs(round(totalangle))*10;
289
anglstep=(totalangle/angldiv);
290
anglall(1:angldiv+1)=maxangle+anglstep*(1:angldiv+1)-anglstep;
291
xcircle(1:angldiv+1)=xcross+R*cos(-anglall(1:angldiv+1)/180*pi);
292
ycircle(1:angldiv+1)=ycross+R*sin(-anglall(1:angldiv+1)/180*pi);
293
plot(xcircle,ycircle,'-b')
294
drawnow
295
296
title(['Segment of circle spreads over ',num2str(totalangle),'°'])
297
298
299
% Accept the chosen circle, try again or give up 
300
301
confirmcircselection = menu(sprintf('Do you want to use this circle as basis?'),...
302
    'Yes','No, try again','Go back to grid-type selection');
303
304
if confirmcircselection==2
305
    close all
306
    imshow(im_grid,'truesize');
307
    circ_grid(FileNameBase, PathNameBase, im_grid);
308
end
309
310
if confirmcircselection==3
311
    close all
312
    gridtypeselection(FileNameBase, PathNameBase, im_grid);
313
end
314
315
if confirmcircselection==1
316
    
317
    prompt = {'Enter the number of intersections between markers on the circle:'};
318
    dlg_title = 'Input for grid creation';
319
    num_lines= 1;
320
    def     = {'30'};
321
    answer = inputdlg(prompt,dlg_title,num_lines,def);
322
    angldiv = str2num(cell2mat(answer(1,1)));
323
    
324
    anglstep=(totalangle/angldiv);
325
    anglall(1:angldiv+1)=maxangle+anglstep*(1:angldiv+1)-anglstep;
326
    
327
    markerxpos(1:angldiv+1)=xcross+R*cos(-anglall(1:angldiv+1)/180*pi);
328
    markerypos(1:angldiv+1)=ycross+R*sin(-anglall(1:angldiv+1)/180*pi);
329
    
330
    plot(markerxpos,markerypos,'ob');
331
    
332
    % Pick the lower bound in the image
333
    title(sprintf('Pick three points lying on the circle in clockwise order. The first and last one define the width of the raster') )
334
    
335
    [x(4,1),y(4,1)]=ginput(1);
336
    hold on
337
    plot(x(1,1),y(1,1),'+r')
338
    
339
    lowboundx=x(4,1);
340
    lowboundy=y(4,1);
341
    
342
    R2=sqrt((xcross-lowboundx(1,1))*(xcross-lowboundx(1,1))+(ycross-lowboundy(1,1))*(ycross-lowboundy(1,1)));
343
    markerxposlb(1:angldiv+1)=xcross+R2*cos(-anglall(1:angldiv+1)/180*pi);
344
    markeryposlb(1:angldiv+1)=ycross+R2*sin(-anglall(1:angldiv+1)/180*pi);
345
    
346
    plot(markerxposlb,markeryposlb,'ob');
347
    
348
    prompt = {'Enter the number of intersections between the upper and lower bound:'};
349
    dlg_title = 'Input for grid creation';
350
    num_lines= 1;
351
    def     = {'5'};
352
    answer = inputdlg(prompt,dlg_title,num_lines,def);
353
    Rdiv = str2num(cell2mat(answer(1,1)));
354
    
355
    Rstep=((R-R2)/Rdiv);
356
    Rall(1:Rdiv+1)=R2+Rstep*(1:Rdiv+1)-Rstep;
357
    
358
    grid_x=ones(Rdiv+1,angldiv+1)*xcross;
359
    grid_y=ones(Rdiv+1,angldiv+1)*ycross;
360
    A=Rall;
361
    B=cos(-anglall(1:angldiv+1)/180*pi);
362
    C=A'*B;
363
    grid_x=grid_x+Rall'*cos(-anglall(1:angldiv+1)/180*pi);
364
    grid_y=grid_y+Rall'*sin(-anglall(1:angldiv+1)/180*pi);
365
    
366
    close all
367
    imshow(im_grid,'truesize');
368
    hold on
369
    plot(grid_x,grid_y,'.b')    
370
    
371
    title(['Selected grid has ',num2str(angldiv*Rdiv), ' rasterpoints'])    % plot a title onto the image
372
373
    
374
    % Do you want to keep the grid?
375
    confirmselection = menu(sprintf('Do you want to use this grid?'),...
376
        'Yes','No, try again','Go back to grid-type selection');
377
    
378
    if confirmselection==1
379
        % Save settings and grid files in the image directory for visualization/plotting later
380
        %         save settings.dat xspacing yspacing xmin_new xmax_new ymin_new ymax_new -ascii -tabs
381
        save grid_x.dat grid_x -ascii -tabs
382
        save grid_y.dat grid_y -ascii -tabs
383
    end
384
    
385
    if confirmselection==2
386
        close all
387
        hold off
388
        imshow(im_grid,'truesize');
389
        circ_grid(FileNameBase, PathNameBase, im_grid);
390
    end
391
    
392
    if confirmselection==3
393
        gridtypeselection(FileNameBase, PathNameBase, im_grid);
394
    end
395
    
396
end
397
398
399
return
400
401
402
403
%-------------------------------
404
%
405
406
function [grid_x,grid_y,FileNameBase,PathNameBase] = rect_grid(FileNameBase, PathNameBase, im_grid);
407
408
title(sprintf('Define the region of interest.  Pick (single click) a point in the LOWER LEFT region of the gage section.\n  Do the same for a point in the UPPER RIGHT portion of the gage section.'))
409
410
[x(1,1),y(1,1)]=ginput(1);
411
hold on
412
plot(x(1,1),y(1,1),'+b')
413
414
[x(2,1),y(2,1)]=ginput(1);
415
hold on
416
plot(x(2,1),y(2,1),'+b')
417
418
drawnow
419
420
xmin = min(x);
421
xmax = max(x);
422
ymin = min(y);
423
ymax = max(y);
424
425
lowerline=[xmin ymin; xmax ymin];
426
upperline=[xmin ymax; xmax ymax];
427
leftline=[xmin ymin; xmin ymax];
428
rightline=[xmax ymin; xmax ymax];
429
430
plot(lowerline(:,1),lowerline(:,2),'-b')
431
plot(upperline(:,1),upperline(:,2),'-b')
432
plot(leftline(:,1),leftline(:,2),'-b')
433
plot(rightline(:,1),rightline(:,2),'-b')
434
435
% closereq
436
437
cd(PathNameBase)
438
439
% Prompt user for grid spacing/resolution
440
prompt = {'Enter horizontal (x) resolution for image analysis [pixels]:', ...
441
        'Enter vertical (y) resolution for image analysis [pixels]:'};
442
dlg_title = 'Input for grid creation';
443
num_lines= 1;
444
def     = {'50','50'};
445
answer = inputdlg(prompt,dlg_title,num_lines,def);
446
xspacing = str2num(cell2mat(answer(1,1)));
447
yspacing = str2num(cell2mat(answer(2,1)));
448
449
% Round xmin,xmax and ymin,ymax "up" based on selected spacing
450
numXelem = ceil((xmax-xmin)/xspacing)-1;
451
numYelem = ceil((ymax-ymin)/yspacing)-1;
452
453
xmin_new = (xmax+xmin)/2-((numXelem/2)*xspacing);
454
xmax_new = (xmax+xmin)/2+((numXelem/2)*xspacing);
455
ymin_new = (ymax+ymin)/2-((numYelem/2)*yspacing);
456
ymax_new = (ymax+ymin)/2+((numYelem/2)*yspacing);
457
458
% Create the analysis grid and show user
459
[x,y] = meshgrid(xmin_new:xspacing:xmax_new,ymin_new:yspacing:ymax_new);
460
[rows columns] = size(x);
461
zdummy = 200.*ones(rows,columns);
462
imshow(FileNameBase)
463
title(['Selected grid has ',num2str(rows*columns), ' rasterpoints'])    % plot a title onto the image
464
hold on;
465
plot(x,y,'+b')
466
467
grid_x=x;
468
grid_y=y;
469
470
% Do you want to keep the grid?
471
confirmselection = menu(sprintf('Do you want to use this grid?'),...
472
    'Yes','No, try again','Go back to grid-type selection');
473
474
if confirmselection==1
475
    % Save settings and grid files in the image directory for visualization/plotting later
476
    save settings.dat xspacing yspacing xmin_new xmax_new ymin_new ymax_new -ascii -tabs
477
    save grid_x.dat x -ascii -tabs
478
    save grid_y.dat y -ascii -tabs
479
    close all
480
    hold off
481
end
482
483
if confirmselection==2
484
    close all
485
    hold off
486
    imshow(im_grid,'truesize');
487
    rect_grid(FileNameBase, PathNameBase, im_grid);
488
end
489
490
if confirmselection==3
491
    close all
492
    hold off
493
    gridtypeselection(FileNameBase, PathNameBase, im_grid);
494
end