Newer
Older
function gtShowSampleSurface(vol, angle, depth_orig, cmap, r_vectors, grains)
Laura Nervo
committed
% gtShowSampleSurface(vol, angle, depth_orig, cmap, r_vectors, grains)
% --------------------------------------------------------------------
% rotate volume by angle, remove layer depth pixels thick from surface
% use volume with surface added (surface label is 800, grains 1-700)
% show what's left
% highlight grains of interest
% show trace of basal planes on section, display inclination, using
% gtPlotHexagon to display crystallographic orientation.
%
% rotation convention for section: starting from 12-6 o'clock on grain map // y, angle
% is measured clockwise. Hence, first cut of interest is around 13.1
% degrees. Depth is removed from the left of the sample, after rotation
%
% Deal with coordinate systems. R-vector orientation data in the instrument
% system: x//beam, z upwards, y towards door.
% In the reconstructed tomography/grain map data, z is downwards, y//x
% instruments, x//y instrument.
% Finally present all information in the frame of the tomo data.
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
vol=imrotate(vol, angle, 'nearest', 'crop');
%find edge of sample
outline_rot=sum(vol==800, 3);
%calculate depth of section
depth=depth_orig+find(sum(outline_rot, 1),1, 'first');
if isempty(depth)
depth=depth_orig;
disp('no outline found - designed for outline==800 - guessing')
end
%get section, transpose for correct orientation and view direction
section=squeeze(vol(:,depth,:))';
%get surface grains not cut by section
section_outline=outline_rot(:,depth)';
%section is vertical as we look at plane through the sample. Thus x, y in
%the section correspond to y z in 3D
if 0 %simple style
section_min=find(section_outline, 1, 'first');
section_max=find(section_outline, 1, 'last');
for x=[1:section_min section_max:size(section, 2)] %along the rows
for y=1:size(section, 1) %down the columns
tmp=vol(x,:,y); %get the line of voxels
tmp(find(tmp==800))=0;%remove outline
if any(tmp~=0)
[v,v,v]=find(tmp, 1, 'first');
section(y,x)=v;
end
end
end
else %more difficult...
sample_dia=find(sum(outline_rot,1), 1, 'last')-find(sum(outline_rot,1), 1, 'first');
if isempty(sample_dia)
sample_dia=max(size(outline_rot));
disp('again, no outline (==800) found')
end
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
%estimate how thick a layer we should look at to get near surface
%grains
a=asin(x-(size(section,2)/2)/(sample_dia/2));
if ~isreal(a)
a=pi/2;
end
layer_thickness=5/cos(a);
layer_thickness=min(round(layer_thickness), 40);
for y=1:size(section, 1) %down the columns
%don't remove the edges of the section
if section(y,x)~=800
tmp=vol(x,:,y); %get the line of voxels
%at this point, is the surface the sample surface, or the section?
sam_edge=max(depth, find(tmp==800, 1, 'first'));
tmp=tmp(sam_edge:sam_edge+layer_thickness); %consider a 21 pixel zone
tmp(find(tmp==800))=0;%remove outline
if any(tmp~=0)
[v,v,v]=find(tmp, 1, 'first');
section(y,x)=v;
end
end
end
end
end
section_id=figure;
imshow(section,[0 max(vol(:))])
%in this image, y tomo is from left to right, z tomo is from top to bottom
%(note this is after the rotation of the sample to display the right
%section)
colormap(cmap)
hold on
set(get(get(section_id, 'CurrentAxes'), 'Title'), 'String',sprintf('Section at angle %0.1f degrees, %d pixels depth',angle,depth_orig'))
%set axis titles - note that because we are looking at a section, x->y,
%y->z
set(get(get(section_id, 'CurrentAxes'), 'yLabel'),'String', 'z axis')
set(get(get(section_id, 'CurrentAxes'), 'xLabel'),'String', 'y axis')
axis equal
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
if 0
dum=[];
for i=1:length(grains)
if isempty(find(section==grains(i), 1))
dum(i)=1;
end
end
grains(find(dum))=[];
end
dum=[];
%plot basal plane traces on the grains of interest
disp('note that the section plane that we look at is y (--> left to right) z (--> top to bottom')
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
for i=1:length(grains)
grain=grains(i);
R=r_vectors(find(r_vectors(:,1)==grain), 2:4);
g=Rod2g(R);
c_axis=(g*[0 0 1]')';
%this is in instrument coordinates... convert to tomo
c_axis=[c_axis(2) c_axis(1) -c_axis(3)];
%test inclination: suitable for fatigue?
a=abs(90-acosd(c_axis(3)));
if 1%a>30 & & a<60 % if the inclination is good
%rotate c_axis by angle (same sense as imrotate
c_axis = [cosd(angle) sind(angle) 0; -sind(angle) cosd(angle) 0;0 0 1]*c_axis';
%section plane is yz, normal // x
section_norm=[1 0 0];
%cross product gives the intesection of these two planes
basal_trace=cross(c_axis, section_norm);
%basal trace: [x y z] in 3d --> x_section=y, y_section=z
[y,x]=ind2sub(size(section), find(section==grain));
point=[mean(x), mean(y)];
point=point+50*basal_trace(2:3);
point(2,:)=point-100*basal_trace(2:3);
plot(point(:,1), point(:,2), 'y')
%other component of this basal plane
%if we can tilt around the 3D y-axis?
%define positive rotation of sample to be anticlockwise about positive
%y-tilt - around the left-right axis of the image as we look at the
%section
y_tilt=atand(-c_axis(1)/c_axis(3));
%x_tilt - this is the one we can see plotted on the image - hence from
%basal trace
x_tilt=atand(-basal_trace(3)/basal_trace(2));
%display for user
fprintf('c_axis of grain %d - in 3D coordinates, after sample rotation, is:\n', grain)
fprintf('%0.2f %0.2f %0.2f\n', c_axis(1), c_axis(2), c_axis(3))
fprintf('x_tilt for grain %d should be %0.2f degrees\n', grain, x_tilt)
fprintf('y_tilt for grain %d should be %0.2f degrees\n', grain, y_tilt)
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
else
dum(i)=1;
end
end
%again, remove grains with non-interesting inclinations
grains(find(dum))=[];
%plot the unit cell representations of the visible grains in second figure
%change gtPlotHexagon to take the two coordinate systems into account
gtPlotHexagon(grains, r_vectors, angle, cmap)