function gtShowSampleSurface(vol, angle, depth_orig, cmap, r_vectors, grains)
% 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
axis square
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)
end % end of function