gtReadBoundaryProperties.m

function gtReadBoundaryProperties(boundaries_structure, vol_boundaries)%, r_vectors)


%
%continue to add to boundaries_structure.mat - orientations,
%disorientation, r-vectors etc.
load parameters;
spacegroup=parameters.acq.spacegroup;

%find the max grainID - ugly but ca va
maxgrain=0;
for i=1:length(boundaries_structure)
    maxgrain=max(maxgrain, boundaries_structure(i).grain1);
    maxgrain=max(maxgrain, boundaries_structure(i).grain2);
end

r_vectors=[];
if exist('r_vectors.mat', 'file')
    load r_vectors
else
    r_vectors=zeros(maxgrain,4);%keep a record to save loading .mats
end


%get sym operators
if spacegroup==225 || spacegroup==229
    % cubic material
    sym = gtGetCubicSymOp;
elseif spacegroup==167 || spacegroup==663 || spacegroup==194
    % hexagonal material - 4x4 element sym operators
    sym=gtGetHexagonalSymOp_sab;
else
    disp('unrecognised spacegroup! quitting...')
    return
end


% loop through the boundaries
for i=1:size(boundaries_structure, 2)

    if boundaries_structure(i).count<4 %can still be useful to calc index planes where we have only one of the grains
        continue
    else

        %Read in the r-vector data for the two grains
        i
        grain1=boundaries_structure(i).grain1;
        grain2=boundaries_structure(i).grain2;

        if grain1==0
            g1=[];
            g1equiv=[];
        else
            if all(r_vectors(grain1,:)==0)
                tmp1=load(sprintf('4_grains/grain%d_/grain%d_.mat',grain1,grain1), 'R_vector');
                R1=tmp1.R_vector;
                r_vectors(grain1,:)=[grain1 R1];
            else
                R1= r_vectors(grain1,2:4);
            end
            boundaries_structure(i).grain1_R_vector=R1;
            g1=Rod2g(R1);
        end

        if grain2==0
            g2=[];
        else
            if all(r_vectors(grain2,:)==0)
                tmp2=load(sprintf('4_grains/grain%d_/grain%d_.mat',grain2,grain2), 'R_vector');
                R2=tmp2.R_vector;
                r_vectors(grain2,:)=[grain2 R2];
            else
                R2= r_vectors(grain2,2:4);
            end
            boundaries_structure(i).grain2_R_vector=R2;
            g2=Rod2g(R2);
        end

% this only works for cubic systems until someone works out how to pdo te
% permutations of the g1 for the hexagonal system
if spacegroup==225 || spacegroup==229

        % calculate the misorientation angle and axis between these grains
        if ~isempty(g1) & ~isempty(g2)
            %if we have both grains
            %need to search symmetry equivelents for the minimum misorientation
            rot_offset=[];
            % warning('g*sym.g maybe should be sym.g*g')
            for j=1:24
                g1equiv = g1*sym(j).g;
                netg = inv(g1equiv)*g2;
                rot_offset(j) = acos((trace(netg)-1)/2);
            end

            dummy = find(rot_offset ==min(rot_offset));
            g1equiv = g1*sym(dummy).g;
            netg = inv(g1equiv)*g2;

            misorientation_angle = (180/pi)*acos((trace(netg)-1)/2);

            [eigVec, eigVal]=eig(netg);
            for p=1:3
                if isreal(eigVec(:,p))  %find the real eigenvector - should be just one
                    misorientation_axis=eigVec(:,p);
                end
            end
            
            boundaries_structure(i).misorientation_axis=misorientation_axis';
            boundaries_structure(i).misorientation_angle=misorientation_angle;
        else
            boundaries_structure(i).misorientation_axis=[];
            boundaries_structure(i).misorientation_angle=[];
        end
        
elseif spacegroup==167 || spacegroup==663 || spacegroup==194
% same calculation but using hexagonal sym operations
     % calculate the misorientation angle and axis between these grains
        if ~isempty(g1) & ~isempty(g2)
            %if we have both grains
            %need to search symmetry equivelents for the minimum misorientation
            rot_offset=[];
            % warning('g*sym.g maybe should be sym.g*g')
            for j=1:length(sym)
                g1equiv = g1*sym(j).g3;
                netg = inv(g1equiv)*g2;
                rot_offset(j) = acos((trace(netg)-1)/2);
            end

            dummy = find(rot_offset ==min(rot_offset));
            g1equiv = g1*sym(dummy).g3;
            netg = inv(g1equiv)*g2;

            misorientation_angle = (180/pi)*acos((trace(netg)-1)/2);

            [eigVec, eigVal]=eig(netg);
            for p=1:3
                if isreal(eigVec(:,p))  %find the real eigenvector - should be just one
                    misorientation_axis=eigVec(:,p);
                    boundaries_structure(i).eigVal=eigVal(p,p);
                end
            end
            
            boundaries_structure(i).misorientation_axis=misorientation_axis';
            boundaries_structure(i).misorientation_angle=misorientation_angle;
            
        else
            boundaries_structure(i).misorientation_axis=[];
            boundaries_structure(i).misorientation_angle=[];
            boundaries_structure(i).eigVal=[];
        end
    
    
end






        %even if we have only one grain, can deal with the boundary plane
        test = vol_boundaries==i;
        [x,y,z]=ind2sub(size(test), find(test));
        [origin, a]=lsplane([x y z]);

        gb_normal=a;
        gb_normal(3)=-gb_normal(3);

        boundaries_structure(i).gb_norm = gb_normal';
        boundaries_structure(i).gb_centre = origin';

        if ~isempty(g1equiv)
            plane1=inv(g1equiv)*gb_normal;
        elseif ~isempty(g1)
            plane1=inv(g1)*gb_normal;
        else
            plane1=[];
        end

        if ~isempty(g2)
            plane2=inv(g2)*gb_normal;
        else
            plane2=[];
        end

        if spacegroup==225 || spacegroup==229
            % if cubic then this is the index plane
            boundaries_structure(i).gb_norm_grain1=plane1';
            boundaries_structure(i).gb_norm_grain2=plane2';

        elseif spacegroup==167 || spacegroup==663 || spacegroup==194
            % if hexagonal, keep cartesian vector, but also
            % convert the cartesian plane normal to hexagonal four index notation
            % ditto for misorientation axis

            boundaries_structure(i).gb_norm_grain1=plane1';
            boundaries_structure(i).gb_norm_grain2=plane2';

            if ~isempty(plane1)
                hkl=plane1;
                hkil=[];
                hkil(1) = hkl(1) - (1/sqrt(3))*hkl(2);
                hkil(2) = (2/sqrt(3))*hkl(2);
                hkil(3) = -(hkil(1) + hkil(2));
                hkil(4) = hkl(3);
                % account for lattice vector lengths
                hkil(1:3) = hkil(1:3) / (2/(sqrt(3)*parameters.acq.latticepar(1)));
                hkil(4)   = hkil(4) / (1/parameters.acq.latticepar(3));

                boundaries_structure(i).gb_norm_hex_grain1=hkil;
            else
                boundaries_structure(i).gb_norm_hex_grain1=[];
            end

            if ~isempty(plane2)
                hkl=plane2;

                hkil=[];
                hkil(1) = hkl(1) - (1/sqrt(3))*hkl(2);
                hkil(2) = (2/sqrt(3))*hkl(2);
                hkil(3) = -(hkil(1) + hkil(2));
                hkil(4) = hkl(3);
                % account for lattice vector lengths
                hkil(1:3) = hkil(1:3) / (2/(sqrt(3)*parameters.acq.latticepar(1)));
                hkil(4)   = hkil(4) / (1/parameters.acq.latticepar(3));
                boundaries_structure(i).gb_norm_hex_grain2=hkil;
            else
                boundaries_structure(i).gb_norm_hex_grain2=[];
            end
% 
%             if ~isempty(boundaries_structure(i).misorientation_axis)
%                 hkl=misorientation_axis';
% 
%                 hkil=[];
%                 hkil(1) = hkl(1) - (1/sqrt(3))*hkl(2);
%                 hkil(2) = (2/sqrt(3))*hkl(2);
%                 hkil(3) = -(hkil(1) + hkil(2));
%                 hkil(4) = hkl(3);
%                 % account for lattice vector lengths
%                 hkil(1:3) = hkil(1:3) / (2/(sqrt(3)*parameters.acq.latticepar(1)));
%                 hkil(4)   = hkil(4) / (1/parameters.acq.latticepar(3));
%                 boundaries_structure(i).misorientation_axis_hex=hkil;
%             end

        else
            disp('crystallography not supported!')
        end

    end % skip small boundaries

end


save('boundaries_structure.mat', 'boundaries_structure')
save('r_vectors.mat', 'r_vectors')