Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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
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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
function 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.
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
for x=[1:size(section, 2)] %along the rows
%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
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
disp(sprintf('c_axis of grain %d - in 3D coordinates, after sample rotation, is:', grain))
disp(sprintf('%0.2f %0.2f %0.2f', c_axis(1), c_axis(2), c_axis(3)))
disp(sprintf('x_tilt for grain %d should be %0.2f degrees', grain, x_tilt))
disp(sprintf('y_tilt for grain %d should be %0.2f degrees', 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)