gtReadBoundaryProperties.m

function gtReadBoundaryProperties(phaseid, boundaries_structure, vol_out, r_vectors)
% gtReadBoundaryProperties(phaseid, boundaries_structure, vol_out, r_vectors)

disp('have changed this to axis convention matlab1 // instX , matlab2 // instY, matlab3 // instZ')
disp('AK, 11/8/2010')
pause(3)

%continue to add to boundaries_structure.mat - orientations,
%disorientation, r-vectors etc.
parameters=[];
load('parameters.mat');
spacegroup = parameters.cryst(phaseid).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

% Get symmetry operators
symm = gtCrystGetSymmetryOperators([], spacegroup);
symm = {symm.g3};

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

    if boundaries_structure(ii).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
        grain1=boundaries_structure(ii).grain1;
        grain2=boundaries_structure(ii).grain2;

        if grain1==0
            g1=[];
            g1equiv=[];
        elseif grain1==-1
            continue;
        else
            R1= r_vectors(grain1,2:4);
            boundaries_structure(ii).grain1_R_vector=R1;
            g1 = gtMathsRod2RotMat(R1);
        end

        if grain2==0
            g2=[];
        elseif grain2==-1
            continue
        else
            R2= r_vectors(grain2,2:4);
            boundaries_structure(ii).grain2_R_vector=R2;
            g2 = gtMathsRod2RotMat(R2);
        end

        % calculate the misorientation angle and axis between these grains
        % TODO: replace by a MISORIENTATION function
        if ~isempty(g1) && ~isempty(g2)
            % If we have both grains
            % Need to search symmetry equivalents for the minimum misorientation
            rot_offset=[];
            % warning('g*symm maybe should be symm*g')
            for jj=1:length(symm)
                g1equiv = g1*symm(jj);
                netg = g1equiv\g2;
                rot_offset(jj) = acos((trace(netg)-1)/2);
            end

            dummy = find(rot_offset ==min(rot_offset));
            g1equiv = g1*symm(dummy);
            netg = g1equiv\g2;

            misorientation_angle = acosd((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(ii).eigVal=eigVal(p,p);
                end
            end
            boundaries_structure(ii).misorientation_axis=misorientation_axis';
            boundaries_structure(ii).misorientation_angle=misorientation_angle;
        else
            boundaries_structure(ii).misorientation_axis=[];
            boundaries_structure(ii).misorientation_angle=[];
            boundaries_structure(ii).eigVal=[];
        end

        %even if we have only one grain, can deal with the boundary plane
        test = vol_out==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(ii).gb_norm = gb_normal';
        boundaries_structure(ii).gb_centre = origin';

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

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

        if spacegroup==225 || spacegroup==229
            % if cubic then this is the index plane
            boundaries_structure(ii).gb_norm_grain1=plane1';
            boundaries_structure(ii).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.cryst.latticepar(1)));
                hkil(4)   = hkil(4) / (1/parameters.cryst.latticepar(3));

                boundaries_structure(ii).gb_norm_hex_grain1=hkil;
            else
                boundaries_structure(ii).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.cryst.latticepar(1)));
                hkil(4)   = hkil(4) / (1/parameters.cryst.latticepar(3));
                boundaries_structure(ii).gb_norm_hex_grain2=hkil;
            else
                boundaries_structure(ii).gb_norm_hex_grain2=[];
            end

            if ~isempty(boundaries_structure(ii).misorientation_axis)
                hkl=boundaries_structure(ii).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.cryst.latticepar(1)));
                hkil(4)   = hkil(4) / (1/parameters.cryst.latticepar(3));
                boundaries_structure(ii).misorientation_axis_hex=hkil;
            end
        else
            disp('crystallography not supported!')
        end
    end % skip small boundaries
end

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

end % end of function