gtPlotHexagon.m

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:
%       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}
%       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
%                    {false}
%       section    = flag to draw a section on the yz-plane to follow
%                    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
%
%     Deals 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.
%
%     Reference hexagon coordinates see gtHexagonalUnitCell
%
%
%     Version 003 16-01-2013 by LNervo

app.grainids   = [];
app.cmap       = [];
app.angle      = 0;
app.tomo_setup = false;
app.reference  = false;
app.section    = false;
app.overlap    = false;
app.view       = [0 90]; % xz plane

app = parse_pv_pairs(app,varargin);


% keep info on grainID in the first column
if size(r_vectors,2) == 3
    r_vectors(:,2:4) = 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

hf = figure();

% find grid size
if length(app.grainids) > 5
    xdim = 5;
    ydim = ceil(length(app.grainids)/xdim);
else
    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));

% Compute all orientation matrices g
all_g = gtMathsRod2RotMat(r_vectors(:, 2:4));

% Express corners3 of each grain in cartesian SAMPLE CS
all_corners3a = gtVectorCryst2Lab(corners3, all_g);

for ii = 1:length(app.grainids)
    grainID = app.grainids(ii);
    jj = find(r_vectors(:, 1) == grainID);

    corners3a = all_corners3a(:, :, jj);

    % allow for the two coordinate systems.  Transform from instrument to
    % reconstructed tomo coordinates: x-->y y-->x z-->-z
    if app.tomo_setup
        corners3a=[corners3a(:,2) corners3a(:,1) -corners3a(:,3)]; 
    end
    % do the rotate to follow the section plane
    corners3b = ([ cosd(app.angle) sind(app.angle) 0; ...
                  -sind(app.angle) cosd(app.angle) 0; ...
                                 0               0 1]*corners3a')';
    % store info on vertices
    vertices(:,:,grainID) = corners3b;
    
    %subplot(ceil(ii,xdim), mod(ii,xdim), ii);
    subplot(ydim, xdim, ii);
    
    % reference hex
    if app.reference && ~app.overlap
        patch('Vertices', corners3, 'Faces', hex_faces, 'FaceColor', 'b', 'FaceAlpha', 0.4);
        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 && app.grainids(1) ~= grainID
        overlapID = app.grainids(1);
        kk = find(r_vectors(:, 1) == overlapID);
        corners3a = all_corners3(:, :, kk);
        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
        pos=min(mean(corners3b(hex_faces(1,:), 1)), mean(corners3b(hex_faces(2,:), 1)));
        section_plane_vertices=[pos 1 1; pos 1 -1; pos -1 1; pos -1 -1];
        section_plane_face=[1 2 4 3];
        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,:), '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');
    % view
    set(get(gcf, 'CurrentAxes'), 'view', app.view);

    axis tight;
    axis equal;
    axis square;
    axis vis3d;

    % activate rotating option
    h = rotate3d(gcf);
    set(h,'Enable','on');
    axis equal;
    
end % end for grains

set(gcf,'UserData',vertices);

hf=gcf;

end % end of function