gt_find_grain_boundaries.m

function [vol_out, boundaries, boundaries_count] = gt_find_grain_boundaries(vol)
% GT_FIND_GRAIN_BOUNDARIES  Starts writing the boundaries_structure.mat file
%     [vol_out, boundaries, boundaries_count] = gt_find_grain_boundaries(vol)
%     -----------------------------------------------------------------------
%     now need to read the grain boundaries from the stack with dilated grains
%     Saves boundaries_structure as tmp.mat file
%
%     INPUT:
%       vol = 'grains' volume contained in 5_reconstruction/volume.mat
%
%     OUTPUT:
%       vol_out          = boundaries volume
%       boundaries       = 
%       boundaries_count = 
%
%     this actually shares functionality with gtMakeBoundariesStructure, which
%     also suppresses single voxel "boundaries".  Should merge...
% 
%     Version 002 by LNervo

boundaries=[];
boundaries_count=0;
boundaries_structure=[];
boundary_number=1;
[size_x, size_y, size_z]=size(vol);
vol_out=zeros(size(vol));

% find the voxels that need to be studied
% dilates colours
vol_dilated = imdilate(vol, ones(3,3,3));
vol_eroded  = imerode(vol,ones(3,3,3));
dif         = vol_dilated-vol_eroded;
todo_voxels = find(dif~=0);


for i=1:length(todo_voxels)
    
    [x,y,z]=ind2sub(size(vol), todo_voxels(i));
    
    % look at voxel neighbourhood
    xstart=x-1;
    xend=x+1;
    ystart=y-1;
    yend=y+1;
    zstart=z-1;
    zend=z+1;
    if xstart==0
        xstart=1;
    end
    if ystart==0
        ystart=1;
    end
    if zstart==0
        zstart=1;
    end
    if xend>size_x
        xend=size_x;
    end
    if yend>size_y
        yend=size_y;
    end
    if zend>size_z
        zend=size_z;
    end
    
    neighbours=vol(xstart:xend, ystart:yend, zstart:zend);
    neighbours=neighbours(:);
    % neighbours(find(neighbours==0))=[]; % don't consider boundaries to zeros
    
    voxel_colour=vol(x,y,z);
    % find any different coloured neighbours
    neighbours=neighbours(find(neighbours~=voxel_colour));
    if ~isempty(neighbours)
        neighbour_colour=min(neighbours); % lowest non-zero neighbour colour
        
        % has this combination of voxel_colour and neighbour_colour already been
        % treated? look for rows in boundaries
        dummy1=find(boundaries(:)==voxel_colour);
        [dummy1,a]=ind2sub(size(boundaries),dummy1);
        dummy2=find(boundaries(:)==neighbour_colour);
        [dummy2,a]=ind2sub(size(boundaries),dummy2);
        
        % any common rows?
        boundary_number=dummy1(find(ismember(dummy1,dummy2)));
        
        if ~isempty(boundary_number) % a previously encountered boundary
            vol_out(x,y,z)=boundary_number;
            boundaries_count(boundary_number)=boundaries_count(boundary_number)+1;
            % add to boundaries_structure
            boundaries_structure(boundary_number).count=boundaries_structure(boundary_number).count+1;
            boundary_number=[];
        else % a new boundary
            boundary_number=size(boundaries,1)+1;
            disp(['now on boundary: ' num2str(boundary_number)])
            vol_out(x,y,z)=boundary_number;
            boundaries(boundary_number,:)=[voxel_colour, neighbour_colour];
            boundaries_count(boundary_number)=1;
            % add to boundaries_structure
            boundaries_structure(boundary_number).grain1=voxel_colour;
            boundaries_structure(boundary_number).grain2=neighbour_colour;
            boundaries_structure(boundary_number).count=1;
            boundary_number=[];
        end
        
    end % if a boundary voxel
    

end % loop over todo voxels

save('vol_out.mat','vol_out');
save('tmp.mat','boundaries_structure');

end % end of function