Newer
Older
Laura Nervo
committed
function [pof, good_plnorms, dsqr]=plot_rodrigues_consistancy_sortgrains(plane_normal,hkl,spacegroup,latticepar,plot_flag)
% [pof, good_plnorms, dsqr]=plot_rodrigues_consistancy_sortgrains(plane_normal,hkl,spacegroup,latticepar,plot_flag)
% -----------------------------------------------------------------------------------------------------------------
% Peter's modification of plot_rodrigues_consistancy_test for sorting grains;
%
% INPUT:
% plane_normal = list of plane normals in sample coords
% hkl = a vector of reflection type eg [1 1 1; 0 0 2; 2 2 0 etc...]
% for hexagonal this this be a 4 index notation
% spacegroup = parameters.cryst(phaseid).spacegroup
% latticepar = lattice parameters
% plot_flag = true to plot figure, false to not plot anything
%
% Version 001 19-10-2007 by AKing
% add consistancy check
% Version 002 19-10-2012 by LNervo
% add spacegroup argument
save_plot=[];
if plot_flag
figure;
hold on;
end
lines.cor=[];
lines.plnorm=[];
consistancy_test=0.2; %select those lines within distance consistancy_test*inc of the intersection (was 0.2)
switch spacegroup
case {225, 229} % ie cubic lattice
inc=0.1; % scan space in these steps (size of the sampling element)
inc_ext=1.1*inc; % extended space to determine precise point of intersection; (should be larger than 0.87*inc)
lim=sqrt(2)-1;
limz=sqrt(2)-1;
for j=1:size(plane_normal,1) % loop of plane_normals
%________calculate h vector___________
%from eq (3.7) in Poulsen
%B is identity matrix in cubic material
%gtnewGetReflections returns matrix of all reflections related by symmetry
%to hkl, using Riso GetCubicSymOp.m
Yoann Guilhem
committed
hkl_reflections = gtGetReflections(hkl(j,:), spacegroup);
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
if size(hkl_reflections,1)~=1%skip if hkl = [0 0 0]
for i=1:size(hkl_reflections,1) % loop of hkl reflections
h = hkl_reflections(i,:);
h = h./(sqrt(h*h'));
%________calculate y vector___________
% y is the normalised scattering vector in sample system
% == plane normal
y = plane_normal(j,:);
%_____calculate rodrigues space vector_____
%formulae (3.16), (3.17) in Poulsen
r_zero = cross(h,y)./(1+(dot(h,y)));
r_direction = (h + y)./(1+(dot(h,y)));
lines.cor=[lines.cor; r_zero r_direction];
lines.plnorm=[lines.plnorm; j];
%max line length needed ~1.5, say +-1.5 for saftey
t=(-1.5:0.01:1.5)';
r_plot = repmat(r_zero,length(t),1) + (t*r_direction);
%remove everything outside of a fundamental region (assume cube)
lim=sqrt(2)-1;
for i=1:3
dummy=find(r_plot(:,i)<-lim | r_plot(:,i)>lim);
r_plot(dummy,:)=[];
end
%remove the corners of this region also
lim=1/sqrt(3);
%calculate distance of each point along each <111> directions
vector=[1 1 1]./sqrt(3);%normalise
vector=repmat(vector,size(r_plot,1),1);
tmp = sum( ((abs(r_plot)).*vector),2); %projections along <111>
dummy=find(tmp>lim);
r_plot(dummy,:)=[];
%collect all plot data for analysis
save_plot=[save_plot; r_plot];
%____plot______
if plot_flag
if all(hkl(j,:)==[1 1 1]) %(111) - red
plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'r');
elseif all(hkl(j,:)==[0 0 2]) %(002) - blue
plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b');
elseif all(hkl(j,:)==[2 2 0]) %(220) - green
plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'g');
elseif all(hkl(j,:)==[3 1 1]) %(311) - black
plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'k');
else %all others in pink
plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m');
end
end
end%loop through all hkl variants
end%if hkl is okay (skip hkl = [0 0 0])
%uncomment this to see line by line in rodrigues space.
%k=waitforbuttonpress %can watch as each is plotted in turn
end%loop through all normals
%%%%%%%end case 225/ 229
case 663
inc=0.08;% scan space in these steps (size of the sampling element)
inc_ext=1.1*inc; % extended space to determine precise point of intersection; (should be larger than 0.87*inc)
lim=1;
limz=2-sqrt(3);
plane_normal_cart=plane_normal;
for j=1:size(plane_normal_cart,1) % loop of plane_normals
%________calculate h vector___________
%from eq (3.7) in Poulsen
%B is identity matrix in cubic material
%gtnewGetReflections returns matrix of all reflections related by symmetry
%to hkl, using Riso GetCubicSymOp.m
Yoann Guilhem
committed
hkil_reflections_hex = gtGetReflections(hkl(j,:), spacegroup);
if size(hkil_reflections_hex,1)~=1%skip if hkl = [0 0 0]
for i=1:size(hkil_reflections_hex,1) % loop of hkil reflections
% passage de l espace direct hexagonal en coordonnees cartesiennes
h(1)= hkil_reflections_hex(i,1) + 0.5 * hkil_reflections_hex(i,2);
h(2)= 3^0.5/2 *hkil_reflections_hex(i,2);
h(3)= hkil_reflections_hex(i,4);
% normalisation of the cartesian coordinates space
% allow for c/a ratio in hexagonal unit cell
Laura Nervo
committed
h(1)=h(1)*2/(sqrt(3)*latticepar(1));
h(2)=h(2)*2/(sqrt(3)*latticepar(1));
h(3)=h(3)/latticepar(3);
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
% normalise to unit vector
h = h./(sqrt(h*h'));
%________calculate y vector___________
% y is the normalised scattering vector in cartesian system
% == plane normal
y = plane_normal_cart(j,:);
y = y./(sqrt(y*y'));
%_____calculate rodrigues space vector_____
%formulae (3.16), (3.17) in Poulsen
r_zero = cross(h,y)./(1+(dot(h,y)));
r_direction = (h + y)./(1+(dot(h,y)));
lines.cor=[lines.cor; r_zero r_direction];
lines.plnorm=[lines.plnorm; j];
%max line length needed ~1.5, say +-1.5 for saftey
t=(-3:0.01:3)';
r_plot = repmat(r_zero,length(t),1) + (t*r_direction);
%remove everything outside of a fundamental region (assume cf Heinz_Neumann acta. crist.)
% base and top of the prism
dummy=find(abs(r_plot(:,1))>=1 | abs(r_plot(:,2))>=1 | abs(r_plot(:,3))>=(2-sqrt(3)) );
r_plot(dummy,:)=[];
% remove facette of the faces
dummy=find((r_plot(:,1)+sqrt(3)*r_plot(:,2)>=2) | (r_plot(:,1)-sqrt(3)*r_plot(:,2))>=2);
r_plot(dummy,:)=[];
dummy=find((-r_plot(:,1)+sqrt(3)*r_plot(:,2)>=2) | (-r_plot(:,1)-sqrt(3)*r_plot(:,2))>=2);
r_plot(dummy,:)=[];
dummy=find((r_plot(:,2)+sqrt(3)*r_plot(:,1))>=2 | (r_plot(:,2)-sqrt(3)*r_plot(:,1))>=2);
r_plot(dummy,:)=[];
dummy=find((-r_plot(:,2)+sqrt(3)*r_plot(:,1))>=2 | (-r_plot(:,2)-sqrt(3)*r_plot(:,1))>=2);
r_plot(dummy,:)=[];
%collect all plot data for analysis
save_plot=[save_plot; r_plot];
%____plot______ % adapted to the snow case
if plot_flag
% keyboard
%plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m');
hkil_hex=hkl;
if all(hkil_hex(j,:)==[0 0 0 2])
plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'r');
%figure(2)
%quiver3(0,0,0,h(1),h(2),h(3),'r')
elseif all(hkil_hex(j,:)==[1 1 -2 0])
plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b');
%figure(2)
%quiver3(0,0,0,h(1),h(2),h(3),'b')
elseif all(hkil_hex(j,:)==[1 -1 0 0])
plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'g');
%figure(2)
%quiver3(0,0,0,h(1),h(2),h(3),'g')
elseif all(hkil_hex(j,:)==[1 -1 0 1])
plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'k');
%figure(2)
%quiver3(0,0,0,h(1),h(2),h(3),'k')
elseif all(hkil_hex(j,:)==[1 1 -2 2])
plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'c');
%figure(2)
%quiver3(0,0,0,h(1),h(2),h(3),'c')
% elseif all(hkil_hex(j,:)==[1 1 -2 1])
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'k');
%
% elseif all(hkil_hex(j,:)==[1 -1 0 2])
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m');
%
% elseif all(hkil_hex(j,:)==[1 1 -2 3])
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'c');
%
% elseif all(hkil_hex(j,:)==[1 -1 0 3])
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'y');
%
% elseif all(hkil_hex(j,:)==[1 1 -2 4])
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'r:');
%
% elseif all(hkil_hex(j,:)==[1 -1 0 4])
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b:');
%
% elseif all(hkil_hex(j,:)==[0 0 0 4])
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'g:');
%
% elseif all(hkil_hex(j,:)==[-2 0 2 1])
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m:');
%
% elseif all(hkil_hex(j,:)==[-2 0 2 3])
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'c:');
%
% elseif all(hkil_hex(j,:)==[-3 1 2 1])
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'y:');
%
% elseif all(hkil_hex(j,:)==[-3 1 2 0])
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'b-');
% else
% plot3(r_plot(:,1),r_plot(:,2),r_plot(:,3),'m');
%
end
% legend(' 0 0 0 -2', '1 1 -2 0', '1 -1 0 0', '1 1 -2 1', '1 -1 0 2','1 1 -2 3', '1 1 -2 4','1 -1 0 4', ' 0 0 0 4','-2 0 2 1','-2 0 2 3','-3 1 2 3','-3 1 2 0');
end
end %loop through all hkl variants
end %if hkl is okay (skip hkl = [0 0 0])
% uncomment this to see line by line in rodrigues space.
% k=waitforbuttonpress %can watch as each is plotted in turn
end%loop through all normals
end % end switch space group
if plot_flag
axis equal
end
lines.cor;
lines.plnorm;
size(lines.cor)
%find intersections
tic
max_score=0;
max_coords=[0 0 0];
score=[];
count=1;
for i=-lim:inc:lim
for j=-lim:inc:lim
for k=-limz:inc:limz
dummy = find(save_plot(:,1)>i & save_plot(:,1)<(i+inc) &save_plot(:,2)>j & save_plot(:,2)<(j+inc) &save_plot(:,3)>k & save_plot(:,3)<(k+inc));
score(count)=length(dummy);
if (score(count)>max_score)
max_score=score;
max_coords=[i j k];
end
count=count+1;
end
end
end
max_score;
max_coords;
element = inc*[0 0 0; 0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1];
max_element = repmat(max_coords,8,1)+element;
max_coords = max_coords+0.5*[inc inc inc]; %this is the approximate solution
if plot_flag
plot3(max_element(:,1),max_element(:,2),max_element(:,3),'go')
plot3(max_coords(:,1),max_coords(:,2),max_coords(:,3),'co')
end
% Determine precise point of intersection inside the sampling element
lines.cor=normlinedir(lines.cor);
line_cand.cor=[];%line equations which pass close to max_coords
line_cand.plnorm=[];%to which plane normal (struct_id entry) does line belong - to ensure only one line from each normal
%for j=1:size(lines.cor,1) % picks lines which pass close to the sampling element: candidate lines
% if (pointtolinedist(max_coords, lines.cor(j,:)) < inc_ext)
% line_cand.cor=[line_cand.cor; lines.cor(j,:)];
% line_cand.plnorm=[line_cand.plnorm; lines.plnorm(j)];
% end
%end
distances=pointtolinedist(max_coords,lines.cor);
dummy=find(distances< inc_ext);
line_cand.cor = lines.cor(dummy,:);
line_cand.plnorm = lines.plnorm(dummy);
% If more than one line from a given plane is present in the set, then
% choose the one that fits the best.
pof_prel=pofintersectexp(line_cand.cor)';
% for k=1:size(plane_normal,1) % delete redundant lines from line_cand
% dis=1e+38;
% m=1;
% while (m <= size(line_cand.cor,1))
% if (line_cand.plnorm(m) == k)
% if pointtolinedist(pof_prel, line_cand.cor(m,:)) < dis
% dis=pointtolinedist(pof_prel, line_cand.cor(m,:));
% dis_line.cor=line_cand.cor(m,:);
% end
% line_cand.cor(m,:)=[];
% line_cand.plnorm(m)=[];
% m=m-1;
% end
% m=m+1;
% end
% if dis ~= 1e+38
% line_cand.cor=[line_cand.cor; dis_line.cor(1:6)];
% line_cand.plnorm=[line_cand.plnorm; k];
% end
% end
% Alternative (cleaner!) loop to find the closest line from a given plane
% to the point of intersection, if more than one is present
distances=pointtolinedist(pof_prel, line_cand.cor);%distance from each candidate line to point
dummy2=[];%list of lines to delete
for k=1:size(plane_normal,1)
dummy=find(plane_normal==k);
if length(dummy)>2%if more than one line, delete all but the closest
%- make a list of those to delete
dummy2 = [dummy2 find(distances~=min(distances(dummy)) & line_cand.plnorm ==k)];
end
end
line_cand.cor(dummy2,:)=[];%do the deletion
line_cand.plnorm(dummy2,:)=[];%do the deletion
%line_cand now contains information only about lines which pass near the
%approx intersection, and only the best line arising from each plane
%normal
%repeat loop, to try and choose more precisely
pof_prel=pofintersectexp(line_cand.cor)'; %least squares best intersection
%find distances of remaining lines to point
distances=pointtolinedist(pof_prel, line_cand.cor);%distance from each candidate line to point
test_distance = consistancy_test*inc;
dummy2=find(distances>test_distance);
%%display the distances and which plane normals they come from
%distances
%line_cand.plnorm
line_cand.cor(dummy2,:)=[];%do the deletion
line_cand.plnorm(dummy2,:)=[];%do the deletion
%%[line_cand.cor,line_cand.plnorm]
%final least squares fit with good lines only, display results
max_coords;
pof_prel;
pof=pofintersectexp(line_cand.cor)'
dsqr=pointtolinedistsqrsum(max_coords, line_cand.cor)
%consistancy check part
good_plnorms = 1:size(plane_normal,1);%index of plane normals
good_plnorms = ismember(good_plnorms,line_cand.plnorm);
if plot_flag
plot3(pof(1),pof(2),pof(3),'mo','markersize',25);
plot3(pof(1),pof(2),pof(3),'m*','markersize',25);
end
toc