Skip to content
Snippets Groups Projects
Commit a60b4b19 authored by Laura Nervo's avatar Laura Nervo Committed by Nicola Vigano
Browse files

Improved a bit the drawing of hexagonal unit cells 


gtPlotHexagon : inputs are r_vectors, varargin (added 'grainids','view' to options)
gtShowSampleSurface : added input 'lookID' to look for a specific value into the volume
gtHexagonalUnitCell : added some drawing options to varargin

Signed-off-by: default avatarLaura Nervo <laura.nervo@esrf.fr>
parent 6e38402d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
6_rendering/gtHexagonalUnitCell.m
+ 110
− 91
View file @ a60b4b19
function data = gtHexagonalUnitCell(varargin)
% GTHEXAGONALUNITCELL Draws the hexagonal unit cell with crystallographic axes
% for Euler angles (0,0,0)
%
% data = gtHexagonalUnitCell(varargin)
% ------------------------------------
% Draws the unit cell using this convention:
% Xs // to <2-1-10> or <100> axis (a-axis)
% Ys // to
%
% Crystallographic axes are Xc, Yc, Zc which corresponds to a,b,c.
% For hexagonal unit cells, there are the equivalents axes a1,a2,a3,c in
% which a1//a, a2 is a1 rotated 120 deg counterclockwise, a3 is a2 rotated
% 120 degrees counterclockwise.
% Sample axes are Xs,Ys,Zs. Please refer to the
% DCTdoc_Crystal_orientation_descriptions.pdf in the DCT_Documentation
% folder for the graintracking code.
%
app.centered = true;
app.draw = true;
app.ratio = 1.5857;
app.ind = 0;
app.vertices = false;
app.numbers = false;
app.patch = false;
app.axes = false;
app.centered = false;
app.draw = true;
app.ratio = 1.5857;
app.ind = 0;
app.vertices = false;
app.numbers = false;
app.patch = false;
app.sampleaxes = false;
app.crystaxes = true;
app.view = [103 20];
app.add_vertices = false;
app = parse_pv_pairs(app, varargin);
% lattice axes' length; c is normalized to 1
c=1;
a=1/app.ratio;
% apothem = c2 = r*sqrt(3)/2
cx=a*sqrt(3);
c2=cx/2;
% calculation of vertices
% first 6 are for 000-1 plane
% first is at X=0 and then counterclockwise
v=0:60:360;
cv=a*cosd(v);
sv=a*sind(v);
x0 = 0*cx+cv;
y0 = 0*cx+sv;
x0 = a*cosd(v);
y0 = a*sind(v);
if app.centered
z0 = repmat(-c/2,1,length(v));
z1 = repmat(+c/2,1,length(v));
......@@ -33,100 +53,99 @@ vertices = [[x0;y0;z0],[x0;y0;z1]]';
vertices(7,:)=[];
vertices(end,:)=[];
if app.add_vertices
vertices(end+1,:) = [0 0 z0(1)];
vertices(end+1,:) = [0 0 z1(1)];
vertices(end+1,:) = [0 0 (z1(1)-z0(1))/3];
vertices(end+1,:) = [0 0 (z1(1)-z0(1))/3*2];
vertices(end+1,:) = [0 0 (z1(1)-z0(1))/4];
vertices(end+1,:) = [0 0 (z1(1)-z0(1))/4*2];
vertices(end+1,:) = [0 0 (z1(1)-z0(1))/4*3];
end
% hexagon face directions
% vertices for patch
% hkil of plane normals
faces_sides = [1 2 8 7; 4 5 11 10; 5 6 12 11; 2 3 9 8; 3 4 10 9; 6 1 7 12];
faces_end = [1 2 3 4 5 6; 7 8 9 10 11 12]; %
hkil = [1 0 -1 0; -1 0 1 0; 0 -1 1 0; 0 1 -1 0; -1 1 0 0; 1 -1 0 0; ...
0 0 0 -1; 0 0 0 1]'; % column vectors
% for crystallographic axes and labels
orig=repmat([0 0 z0(1)],4,1);
xcrys=[1;-0.5;-0.5;0];
ycrys=[0;0.866;-0.866;0];
zcrys=[z0(1);z0(1);z0(1);z1(1)+0.4];
lcrys={'a1 [2-1-10]';...
'a2 [-12-10]';...
'a3 [-1-120]';...
'c [0001]'};
if app.draw
figure;
% basal planes
line(x0,y0,z0);
line(x0,y0,z1);
hold on;
if app.vertices
plot3(x0,y0,z0,'.k')
plot3(x0(1),y0(1),z0(1),'.b')
plot3(x0,y0,z1,'.k')
plot3(x0(1),y0(1),z1(1),'.b')
end
% prismatic planes
for ii=1:length(x0)
plot3([x0(ii) x0(ii)],[y0(ii) y0(ii)],[z0(ii) z1(ii)], 'b-', 'LineWidth', 1)
end
quiver3(0,0,z0(1),1,0,0,0,'k')
quiver3(0,0,z0(1),-0.5,0.866,0,0,'k')
quiver3(0,0,z0(1),-0.5,-0.866,0,0,'k')
quiver3(0,0,z0(1),0,0,1+0.4,0,'k')
text(1+0.1,0,z0(1),'a1','Color',[0 0 0])
text(-0.5-0.1,0.866+0.1,z0(1),'a2','Color',[0 0 0])
text(-0.5-0.1,-0.866-0.1,z0(1),'a3','Color',[0 0 0])
text(0,0,z1(1)+0.5,'c','Color',[0 0 0])
if app.numbers
for ii=1:12
text(vertices(ii,1),vertices(ii,2),vertices(ii,3)+0.1,num2str(ii),'Color',[0 0 0])
end
% cryst axes
if app.crystaxes
quiver3(orig(:,1),orig(:,2),orig(:,3),xcrys,ycrys,zcrys,0,'k')
text(xcrys-0.1,ycrys+0.1,zcrys,lcrys,'Color',[0 0 0])
end
view(103,20)
if app.axes
% lab system : x y z vectors
quiver3(-1,0,0,2,0,0,0,'r');
quiver3(0,-1,0,0,2,0,0,'k');
quiver3(0,0,-1,0,0,2,0,'b');
xlabel('x')
ylabel('y')
zlabel('z')
view(app.view)
if app.sampleaxes
% sample reference system : Xs, Ys, Zs vectors
quiver3(-1.5,0,0,3,0,0,0,'r');
quiver3(0,-1.5,0,0,3,0,0,'g');
quiver3(0,0,z0(1)-0.5,0,0,(z1(1)-z0(1))+1.5,0,'b');
text(1.5+0.1,0,0,'Xs','Color',[1 0 0])
text(0,1.5+0.1,0,'Ys','Color',[0 1 0])
text(0,0,z1(1)+1.1,'Zs','Color',[0 0 1])
end
end
% hexagon face directions
% +yz -yz -xz +xz +xy -xy
% +x -x -y +y +z -z
faces_sides = [1 2 8 7; 4 5 11 10; 5 6 12 11; 2 3 9 8; 3 4 10 9; 6 1 7 12];
faces_end = [1 2 3 4 5 6; 7 8 9 10 11 12]; % -
hkil = [1 0 -1 0; -1 0 1 0; 0 -1 1 0; 0 1 -1 0; -1 1 0 0; 1 -1 0 0; ...
0 0 0 -1; 0 0 0 1]'; % column vectors
if app.draw && app.patch
if app.ind >=1 && app.ind <=6
patch('Vertices', vertices, 'Faces', faces_sides(app.ind,:), ...
'FaceVertexCData', repmat([0 0.8 0],length(app.ind),1), 'FaceColor', 'flat');
elseif app.ind==7 || app.ind==8
patch('Vertices', vertices, 'Faces', faces_end(app.ind-6,:), ...
'FaceVertexCData', repmat([0 0.8 0],length(app.ind-6),1), 'FaceColor', 'flat');
elseif app.ind==0
patch('Vertices', vertices, 'Faces', faces_sides, ...
'FaceVertexCData', repmat([0 0.8 0],6,1), 'FaceColor', 'flat');
patch('Vertices', vertices, 'Faces', faces_end, ...
'FaceVertexCData', repmat([0 0.8 0],2,1), 'FaceColor', 'flat');
if app.patch
if app.ind >=1 && app.ind <=6
patch('Vertices', vertices, 'Faces', faces_sides(app.ind,:), ...
'FaceVertexCData', repmat([0 0.8 0],length(app.ind),1), 'FaceColor', 'flat');
elseif app.ind==7 || app.ind==8
patch('Vertices', vertices, 'Faces', faces_end(app.ind-6,:), ...
'FaceVertexCData', repmat([0 0.8 0],length(app.ind-6),1), 'FaceColor', 'flat');
elseif app.ind==0
patch('Vertices', vertices, 'Faces', faces_sides, ...
'FaceVertexCData', repmat([0 0.8 0],6,1), 'FaceColor', 'flat');
patch('Vertices', vertices, 'Faces', faces_end, ...
'FaceVertexCData', repmat([0 0.8 0],2,1), 'FaceColor', 'flat');
end
end
end
vertices(end+1,:) = [0 0 z0(1)];
vertices(end+1,:) = [0 0 z1(1)];
axis equal
set(gca,'Visible','off')
vertices(end+1,:) = [0 0 (z1(1)-z0(1))/3];
vertices(end+1,:) = [0 0 (z1(1)-z0(1))/3*2];
vertices(end+1,:) = [0 0 (z1(1)-z0(1))/4];
vertices(end+1,:) = [0 0 (z1(1)-z0(1))/4*2];
vertices(end+1,:) = [0 0 (z1(1)-z0(1))/4*3];
if app.draw
axis equal
set(gca,'Visible','off')
if app.numbers
for ii=13:14
text(vertices(ii,1),vertices(ii,2),vertices(ii,3)+0.1,num2str(ii),'Color',[0 0 0])
end
text(vertices(:,1),vertices(:,2),vertices(:,3)+0.1,num2str([1:length(vertices)]'),'Color',[0 0 0])
end
if app.vertices
for ii=13:length(vertices)
plot3(vertices(ii,1),vertices(ii,2),vertices(ii,3),'.k')
end
plot3(vertices(:,1),vertices(:,2),vertices(:,3),'.k')
end
end
end % end draw
if nargout == 1
data.vertices = vertices;
data.vertices = vertices;
data.faces_sides = faces_sides;
data.faces_end = faces_end;
data.hkil = hkil;
data.hf = gcf;
data.faces_end = faces_end;
data.hkil = hkil;
data.orig = orig;
data.xcrys = xcrys;
data.ycrys = ycrys;
data.zcrys = zcrys;
data.lcrys = lcrys;
if app.draw
data.hf = gcf;
end
end
end % end of function
6_rendering/gtPlotHexagon.m
+ 54
− 45
View file @ a60b4b19
function hf = gtPlotHexagon(grains, r_vectors, varargin)
% hf = gtPlotHexagon(grains, r_vectors, varargin)
% -----------------------------------------------
function hf = gtPlotHexagon(r_vectors, varargin)
% hf = gtPlotHexagon(r_vectors, varargin)
% ---------------------------------------
% Makes a figure containing subplots of the unit cells of the selected grains
%
% INPUT:
% grains = IDs of grains of interest <double 1xM>
% r_vectors = orientations for all the grains <double Nx4>
%
% OPTIONAL INPUT (varargin):
% grainids = IDs of grains of interest <double 1xM>
% cmap = colour map for all the grains (rgb) <double N+1x3> {[1 0 0]}
% angle = rotation around z-axis to visualize orientations from the
% section plane as defined in gtShowSampleSurface {0}
% hf = figure handle if already existing []
% tomo_setup = flag to match tomo reference system: x-->y y-->x z-->-z
% {false}
% reference = flag to plot the reference unit cell for each grain orientation
......@@ -20,9 +19,10 @@ function hf = gtPlotHexagon(grains, r_vectors, varargin)
% gtShowSampleSurface {false}
% overlap = flag to overlap the first grainID in the list to all the
% other grains {true}
% view = azimuth and elevation for viewing in 3D {[0 180]} xz plane
%
% OUTPUT:
% hf = figure handle
% hf = figure handle
%
% Deals with coordinate systems. R-vector orientation data in the instrument
% system: x//beam, z upwards, y towards door.
......@@ -34,58 +34,65 @@ function hf = gtPlotHexagon(grains, r_vectors, varargin)
%
% Version 003 16-01-2013 by LNervo
app.grainids = [];
app.cmap = [];
app.angle = 0;
app.hf = [];
app.tomo_setup = false;
app.reference = false;
app.section = false;
app.overlap = true;
app.overlap = false;
app.view = [0 90]; % xz plane
app = parse_pv_pairs(app,varargin);
if isempty(app.cmap)
%if no colourmap passed in, use red
app.cmap = zeros(size(r_vectors, 1), 3);
app.cmap(:,1) = 1;
%if overlapping,use different colour
if app.overlap
app.cmap(1,:) = [0 1 0];
end
end
if isempty(app.hf)
app.hf = figure();
end
% keep info on grainID in the first column
if size(r_vectors,2) == 3
r_vectors(:,2:4) = r_vectors;
r_vectors(:,1) = (1:length(r_vectors))';
r_vectors(:,1) = (1:size(r_vectors,1))';
end
if isempty(app.grainids)
app.grainids = 1:size(r_vectors,1);
end
if isempty(app.cmap)
app.cmap = zeros(size(r_vectors,1), 3);
app.cmap(:,1) = 1;
app.showbar = false;
else
app.showbar = false;
if size(app.cmap,1)==1
app.cmap = repmat(app.cmap,size(r_vectors,1),1);
end
% remove bkg color
if all(app.cmap(1,:) == 0) && size(app.cmap,1) == size(r_vectors,1)+1
app.cmap(1,:) = [];
end
end
% patch object hexagonal prism
data = gtHexagonalUnitCell('ratio',2,'draw',false);
corners3 = data.vertices(1:12,:);
rec_faces = data.faces_sides;
hex_faces = data.faces_end;
vertices = zeros(12,3,length(r_vectors));
hf = figure();
% find grid size
if length(grains) > 5
if length(app.grainids) > 5
xdim = 5;
ydim = ceil(length(grains)/xdim);
ydim = ceil(length(app.grainids)/xdim);
else
xdim = length(grains);
xdim = length(app.grainids);
ydim = 1;
end
set(hf, 'Position', [1 250*ydim 200*xdim 250*ydim])
% patch object hexagonal prism
data = gtHexagonalUnitCell('ratio',2,'draw',false,'centered',true);
corners3 = data.vertices(1:12,:);
rec_faces = data.faces_sides;
hex_faces = data.faces_end;
vertices = zeros(12,3,size(r_vectors,1));
set(app.hf, 'Position', [1 250*ydim 200*xdim 250*ydim])
% apply the R-vector
for ii=1:length(grains)
for ii=1:length(app.grainids)
grainID = grains(ii);
grainID = app.grainids(ii);
% takes the R_vector
R = r_vectors(find(r_vectors(:,1)==grainID),2:4);
% crystal to sample rotation matrix == g^-1
......@@ -113,11 +120,11 @@ for ii=1:length(grains)
patch('Vertices', corners3, 'Faces', rec_faces, 'FaceColor', 'b', 'FaceAlpha', 0.4)
end
% if overlap the first grain is used as reference for the others
if app.overlap && grains(1) ~= grainID
overlapID = grains(1);
if app.overlap && app.grainids(1) ~= grainID
overlapID = app.grainids(1);
corners3a = (Rod2g( r_vectors(find(r_vectors(:,1)==overlapID),2:4) )*corners3')';
patch('Vertices', corners3a, 'Faces', hex_faces, 'FaceColor', app.cmap(overlapID+1,:), 'FaceAlpha', 0.4)
patch('Vertices', corners3a, 'Faces', rec_faces, 'FaceColor', app.cmap(overlapID+1,:), 'FaceAlpha', 0.4)
patch('Vertices', corners3a, 'Faces', hex_faces, 'FaceColor', app.cmap(overlapID,:), 'FaceAlpha', 0.4)
patch('Vertices', corners3a, 'Faces', rec_faces, 'FaceColor', app.cmap(overlapID,:), 'FaceAlpha', 0.4)
end
% reference background like the section plane
if app.section
......@@ -127,25 +134,27 @@ for ii=1:length(grains)
patch('Vertices', section_plane_vertices, 'Faces', section_plane_face, 'FaceColor', 'k', 'FaceAlpha', 0.4)
end
% including the rotation for sectioning
patch('Vertices', corners3b, 'Faces', hex_faces, 'FaceColor', app.cmap(grainID+1,:), 'FaceAlpha', 1)
patch('Vertices', corners3b, 'Faces', rec_faces, 'FaceColor', app.cmap(grainID+1,:), 'FaceAlpha', 1)
patch('Vertices', corners3b, 'Faces', hex_faces, 'FaceColor', app.cmap(grainID,:), 'FaceAlpha', 1)
patch('Vertices', corners3b, 'Faces', rec_faces, 'FaceColor', app.cmap(grainID,:), 'FaceAlpha', 1)
% a bit of labeling
set(get(get(gcf, 'CurrentAxes'), 'Title'), 'String', sprintf('grain %d unit cell', grainID))
set(get(get(gcf, 'CurrentAxes'), 'xLabel'),'String', 'x axis')
set(get(get(gcf, 'CurrentAxes'), 'yLabel'),'String', 'y axis')
set(get(get(gcf, 'CurrentAxes'), 'zLabel'),'String', 'z axis')
% yz-plane view
set(get(gcf, 'CurrentAxes'), 'view', [90 0])
% view
set(get(gcf, 'CurrentAxes'), 'view', app.view)
axis tight;
axis equal;
axis square;
axis vis3d;
axis equal;
% activate rotating option
h = rotate3d(gcf);
set(h,'Enable','on');
axis equal;
end % end for grains
set(gcf,'UserData',vertices)
......
6_rendering/gtShowSampleSurface.m
+ 4
− 4
View file @ a60b4b19
function gtShowSampleSurface(vol, angle, depth_orig, cmap, r_vectors, grains)
% gtShowSampleSurface(vol, angle, depth_orig, cmap, r_vectors, grains)
function gtShowSampleSurface(vol, angle, depth_orig, cmap, r_vectors, grains, lookID)
% gtShowSampleSurface(vol, angle, depth_orig, cmap, r_vectors, grains, lookID)
% --------------------------------------------------------------------
% rotate volume by angle, remove layer depth pixels thick from surface
% use volume with surface added (surface label is 800, grains 1-700)
......@@ -22,7 +22,7 @@ vol=imrotate(vol, angle, 'nearest', 'crop');
%find edge of sample
outline_rot=sum(vol==800, 3);
outline_rot=sum(vol==lookID, 3);
%calculate depth of section
depth=depth_orig+find(sum(outline_rot, 1),1, 'first');
if isempty(depth)
......@@ -168,6 +168,6 @@ 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)
gtPlotHexagon(r_vectors, 'grainids', grains, 'angle', angle, 'cmap', cmap)
end % end of function
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment