gtAnalyseBoundaryIndexPlane.m

function [voldev, volhkl, poles]=gtAnalyseBoundaryIndexPlane(bounds_vol, grain_vol, boundaries_structure, boxsize, list)

%think about analysing the local variations in character of boundaries, in
%terms of their index plane.
%consider a local neighbourhood of <boxsize> aound each point on the
%boundary.  Calculate how far each section of boundary is from a low index
%plane, and which boundary it is.

%this is information which can also be usefully displayed on a pole figure,
%giving a pole figure weighted by area (amount of boundary) rather than by
%frequency of boundaries.  Collect this info at the same time.

%considering every voxel if very slow and unnecessary.  Need to do some
%kind of coarser grid.  How to do this, knowing that the boundaries are not
%necessarily flat?

%looking at a (very) few examples - it seems that an area of around 10x10
%voxels is needed to avoid a lot of noise in the data.  However, it can be
%nice to sample more points than just one every 10x10 voxels.  Really want
%a grid spacing - ie a list of points to look at - and then a sample area - around each point.
%exporting a long list of pole data makes gtMakePoleFigure very slow to
%run.  Would be better to collect the data on the fly in the form of the
%phi/psi/density or x/y/density variable inside gtMPF, and export this.


if isempty(list)
    list=1:max(bounds_vol(:));
    disp('doing all boundaries...');
end

voldev=zeros(size(bounds_vol));%output vol
volhkl=zeros(size(bounds_vol));%output vol
poles=[];

% Low index hkls to consider
hkls=[1 0 0; 1 1 0; 1 1 1; 2 1 0; 2 1 1; 2 2 1; 3 1 0; 3 1 1; 3 2 0; 3 2 1; 3 2 2; 1 1 4];
% Expand to all variants
all_hkls=[];
hkl_type=[];
for ii = 1:size(hkls,1)
    all_hkls = [all_hkls ; gtGetReflections(hkls(ii, :))];
    hkl_type = [hkl_type; repmat(ii, size( gtGetReflections(hkls(ii,:)), 1),1)]; %list of hkl type
end

% Normalise hkl vectors
tmp = sqrt(sum((all_hkls.*all_hkls),2));
normalised_hkls = all_hkls./(repmat(tmp,1,3));

% Compute all orientation matrices g
all_g = gtMathsRod2RotMat(R);

for ii = 1:length(list) % list(i) is the id of the boundary in the volume
  
    fprintf('doing boundary %d... %d of %d', list(ii), ii, length(list));

    % Get boundaries_structure data
    grainids=[boundaries_structure(list(ii)).grain1 boundaries_structure(list(ii)).grain2];
    R=[boundaries_structure(list(ii)).grain1_R_vector; boundaries_structure(list(ii)).grain2_R_vector];
  
    % Skip external boundaries / small boundaries
    if any(grainids==0) || boundaries_structure(list(ii)).count<4
        continue;
    end

    warning('doing one half only');
    for jj = 2; %do one half of the boundary at a time
  
        % Pick out the boundary of interest
        [y,x,z]=ind2sub(size(bounds_vol),find(bounds_vol==list(ii) & grain_vol==grainids(jj))); %voxels of boundary

        % Could reduce points to be considered here by using something like:
        % find(mod(x,2)==0 & mod(y,2)==0 & mod(z,2)==0)
        % or by using a grid of the boxsize (/2?) to create subregions

        g = all_g(:, :, ii);

        % Loop through this volume, calculating local character
        % may want to do this in larger steps

        for k = 1:length(x)
            xx=x(k); yy=y(k); zz=z(k);   
        
            % from gtReadBoundaryProperties - reference
            %[x,y,z]=ind2sub(size(test), find(test));
            %[origin, a]=lsplane([x y z]);
            %gb_normal=a;
            %gb_normal(3)=-gb_normal(3);

            dummy=find(x>=xx-boxsize & x<=xx+boxsize & y>=yy-boxsize & y<=yy+boxsize & z>=zz-boxsize & z<=zz+boxsize);

            if length(dummy)>3
                [origin, a]=lsplane([y(dummy) x(dummy) z(dummy)]);
                a(3)=-a(3);
                local_normal=a; % Should be equivalent to previous

                % Express normal in crystallographic axes
                plane1 = gtVectorLab2Cryst(local_normal, g);

                %% Find closest low index plane
                %dev=acosd(dot(repmat(plane1',size(normalised_hkls,1),1), normalised_hkls,2));
                %[ang, ind]=min(dev);
                %hkl=hkl_type(ind);

                %% Collect pole data
                poles=[poles; plane1'];

                %% Write output - if there is no low index plane close,
                %% assign 15 - non low index
                %if ang<5
                %    voldev(yy,xx,zz)=ang;
                %    volhkl(yy,xx,zz)=hkl;
                %else
                %    voldev(yy,xx,zz)=ang;
                %    volhkl(yy,xx,zz)=15;
                %end
            end
        end
    end
end

% Apply the symmetry operators to pole data
symm = gtCrystGetSymmetryOperators('cubic');
poles2 = [];
for ii = 1:length(symm)
    poles2 = [poles2; poles*symm(ii).g3];
end

% Add this to correspond to line in gtMakePoleFigure
poles2=[poles2; -poles2];
%poles=unique(poles2, 'rows');
poles=poles2;