gtAnalyseCrackPlane.m

function [data, vol, map, surface_view] = gtAnalyseCrackPlane(crack)
%returns the crack plane data as a list [x y z nx ny nz]
%also returns a labelled volume and tle corresponding colourmap
%surface_view is a top down image

%read the x,y,z points of the crack
[y,x,z]=ind2sub(size(crack),find(crack)); %voxels of crack

boxsize=5;
%data out =[x y z nx ny nz]
data=zeros(length(find(crack)), 6);

for i=1:length(x)
    xx=x(i); yy=y(i); zz=z(i);
    dummy=find(x>xx-boxsize & x<xx+boxsize & y>yy-boxsize & y<yy+boxsize & z>zz-boxsize & z<zz+boxsize);

    if length(dummy)>3
        %this is right in tomography coordinates (ie as in volview, x l->r, y
        %top->bot, z up->down)
        [origin, a]=lsplane([x(dummy) y(dummy) z(dummy)]);

        %convert to instrument coordinates to agree with gtCaxisCmap conventions
        a=a([2 1 3]);
        a(3)=-a(3);

        data(i, :)=[xx yy zz a(1) a(2) a(3)];

        if mod(i, 1000)==0
            100*i/length(x)
        end
    end
end

%cut the empty bits of data
data(find(data(:,1)==0), :)=[];

%now produce a labelled volume, and the corresponding colormap
%volume
vol=zeros(size(crack));
for i=1:length(data)
    vol(data(i,2), data(i,1), data(i,3))=i;
end

% Colormap
rgb = gtVectorOrientationColor(data(:, 4:6));

%add black background
map=zeros(size(data,1)+1, 3);
map(2:end, :)=rgb; %ie rgb plus zero background

%read out "bird's eye" view of coloured crack;
surface_view=zeros(size(crack, 1), size(crack, 2));
for y=1:size(crack, 1)
    for x=1:size(crack, 2)
        a=squeeze(vol(y,x,:));
        a=a(find(a, 1, 'first'));
        if ~isempty(a)
            surface_view(y,x)=a;
        end
    end
end

end % end of function