function [normalout,coeffout,out,Center] = FitPlaneDB(x_coord,y_coord,z_coord,supinf,displayflag)%, displayflag,clor,displayflag2,transp)
%FitPlaneDB This funcion will use a procrustes analysis to fit a plane to
%the inferior and superior endplate at each vertebral layer.
%The output normalout defines the 3 vectors that describe the plane in
%space

%This function is flexible for both superior and inferior vertebral
%endplate.  The superior or inferior endplate will be evaluated as
%described by the variable 'supinf'
if supinf=='S'
    x_coord_=x_coord(:,1:2,:);
    y_coord_=y_coord(:,1:2,:);
    z_coord_=z_coord(:,1:2,:);
elseif supinf=='I'
    x_coord_=x_coord(:,3:4,:);
    y_coord_=y_coord(:,3:4,:);
    z_coord_=z_coord(:,3:4,:);
end

%For each level
for i=1:size(x_coord,1)
    %Have temporary variables to for x,y,z coordinates
    x=x_coord_(i,:,:); y=y_coord_(i,:,:); z=z_coord_(i,:,:);
    
    %Get rid of slices that dont have ROIs in them 
    a=find(sum(x)==0); x(:,:,a)=[]; a=find(sum(y)==0); y(:,:,a)=[];
    a=find(sum(z)==0); z(:,:,a)=[];
    x=reshape(x,size(x,3)*2,1); y=reshape(y,size(y,3)*2,1);
    z=reshape(z,size(z,3)*2,1);
    X=[x y z];  Center(i,:)=mean(X)';%Centroid of endplate
    
    %Principle coefficient analysis
    [coeff,score,roots] = princomp(X);
    basis = coeff(:,1:2);%2 vectors that are "square" to the plane
    normal(:,1) = coeff(:,3);%normal vector
    [n,p] = size(X);
    meanX = mean(X,1);
    Xfit = repmat(meanX,n,1) + score(:,1:2)*coeff(:,1:2)';
    
    %Set up the surfaces for the figure if displayflag is on
    [xgrid,ygrid] = meshgrid(linspace(min(Xfit(:,1)),max(Xfit(:,1)),5), ...
        linspace(min(Xfit(:,2)),max(Xfit(:,2)),5));
    zgrid=(1/normal(3)).*(meanX*normal-(xgrid.*normal(1)+ygrid.*normal(2)));
    r = 1.2;%radius of spheres in plot
    [spx,spy,spz]=sphere;
    spx = spx*r;
    spy = spy*r;
    spz = spz*r;
    
    %Output of function, 
    out{i}=[xgrid ygrid zgrid];
    normalout{i}=normal;
    coeffout{i}=coeff';
    
    %plotting
    X=[x y -z];%for asthetics
    if displayflag
        
        for j = 1:length(X)
            surf(spx+X(j,1),spy+X(j,2),spz+X(j,3),'Facecolor','Green','Edgecolor','none')
            hold on
            %             camlight(100,45);
            camlight 'headlight';
            lighting phong
        end
        surf(xgrid,ygrid,-zgrid,'EdgeColor',[0 0 0]);
        %quiver3(Center(1),Center(2),-Center(3),normal(1), normal(2), -normal(3),20)
        xlabel('x'), ylabel('y'), zlabel('z')
    end
end



