Newer
Older
function gtShowSampleSurface(vol, angle, depth_orig, cmap, r_vectors, grains, lookID)
% gtShowSampleSurface(vol, angle, depth_orig, cmap, r_vectors, grains, lookID)
Laura Nervo
committed
% --------------------------------------------------------------------
% 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==lookID, 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));
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
for ii = 1:length(grains)
if isempty(find(section==grains(ii), 1))
dum(ii) = 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('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');
% Compute all orientation matrices g
all_g = gtMathsRod2RotMat(r_vectors(:, 2:4));
% Express every c axis in cartesian SAMPLE CS
c_axes = gtVectorCryst2Lab([0 0 1], all_g);
for ii = 1:length(grains)
grainID = grains(ii);
jj = find(r_vectors(:, 1) == grainID);
c_axis = c_axes(jj, :);
% 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==grainID));
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', grainID);
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', grainID, x_tilt);
fprintf('y_tilt for grain %d should be %0.2f degrees\n', grainID, y_tilt);
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');
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(r_vectors, 'grainids', grains, 'angle', angle, 'cmap', cmap)