/normxcorr/trunk

To get this branch, use:
bzr branch http://suren.me/webbzr/normxcorr/trunk

« back to all changes in this revision

Viewing changes to grid_generator.m

  • Committer: Suren A. Chilingaryan
  • Date: 2009-01-15 13:50:29 UTC
  • Revision ID: csa@dside.dyndns.org-20090115135029-wleapicg9a4593tp
Initial import

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
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