function [IVDAngle,InfrRot,Supr,xgrid,ygrid,zgrid,DAnterior, DCenter, DPosterior] =...
    getIVDAngles(j,x_coord, y_coord, z_coord, displayflag, InfCenter, SupCenter, InfPlaneN, SupPlaneN,...
    SuprXgrid,SuprYgrid,SuprZgrid,InfrXgrid,InfrYgrid,InfrZgrid)

% FitPlane(x_coord(j,3:4,:),y_coord(j,3:4,:),z_coord(j,3:4,:),displayflag,'r.',displayflag,0.1);
% FitPlane(x_coord(j+1,1:2,:),y_coord(j+1,1:2,:),z_coord(j+1,1:2,:),displayflag,'g.',displayflag,0.5);
a = 20;
if displayflag
    
    vectarrow([SupCenter(1,j+1) SupCenter(2,j+1) SupCenter(3,j+1)],...    % L(j+1) Superior
        [a*SupPlaneN(1,j+1)+SupCenter(1,j)...
        a*SupPlaneN(2,j+1)+SupCenter(2,j)...
        a*SupPlaneN(3,j+1)+SupCenter(3,j)],'c');
    hold on
    vectarrow([InfCenter(1,j) InfCenter(2,j) InfCenter(3,j)],...    % Lj Inferior
        [a*InfPlaneN(1,j)+InfCenter(1,j)...
        a*InfPlaneN(2,j)+InfCenter(2,j)...
        a*InfPlaneN(3,j)+InfCenter(3,j)],'r');
    hold on
end

Supr = reshape(SuprXgrid,prod(size(SuprXgrid)),1);
Supr = [Supr reshape(SuprZgrid,prod(size(SuprZgrid)),1)];
Supr = [Supr reshape(SuprYgrid,prod(size(SuprYgrid)),1)];

Infr = reshape(InfrXgrid,prod(size(InfrXgrid)),1);
Infr = [Infr reshape(InfrZgrid,prod(size(InfrZgrid)),1)];
Infr = [Infr reshape(InfrYgrid,prod(size(InfrYgrid)),1)];


if sum(sum(Infr,1),2)==0 | sum(sum(Supr,1),2)==0
    
    IVDAngle=0;
    InfrRot=[];
    Supr=[];
    xgrid=[];
    ygrid=[];
    zgrid=[];
    DAnterior=0;
    DCenter=0;
    DPosterior=0;
    
else
    
%     [psi1, theta1, phi1] = getAngle(Infr,Supr);
%     [psi1, theta1, phi1]*180/pi
%     
%     %%% Get Phi
%     % Use to remove sagittal
%     [InfrRot,NewInfPlaneN(:,j),InfNewCenter(:,j),xgrid,ygrid,zgrid] = rotate(Infr, [0,-theta1,0],displayflag, 'b.',displayflag);
%     [psi1, theta1, phi1] = getAngle(InfrRot,Supr);
%     [psi1, theta1, phi1]*180/pi
%     
%     
%     % Remove Theta
%     [InfrRot,NewInfPlaneN(:,j),InfNewCenter(:,j),xgrid,ygrid,zgrid] = rotate(InfrRot, [0,-theta1,0],displayflag, 'b.',displayflag);
%     [psi1, theta1, phi1] = getAngle(InfrRot,Supr);
%     [psi1, theta1, phi1]*180/pi
%     
%     CorAngle=-phi1*180/pi;
%     
%     
%     
%     %%% Get Theta
%     % Use to remove sagittal
%     [InfrRot,NewInfPlaneN(:,j),InfNewCenter(:,j),xgrid,ygrid,zgrid] = rotate(Infr, [-phi1,0,-psi],displayflag, 'b.',displayflag);
%     [psi1, theta1, phi1] = getAngle(InfrRot,Supr);
%     [psi1, theta1, phi1]*180/pi
%     
%     
%     % Remove Theta
%     [InfrRot,NewInfPlaneN(:,j),InfNewCenter(:,j),xgrid,ygrid,zgrid] = rotate(InfrRot, [-phi1,0,-psi1],displayflag, 'b.',displayflag);
%     [psi1, theta1, phi1] = getAngle(InfrRot,Supr);
%     [psi1, theta1, phi1]*180/pi
%     
%     AxAngle=-theta1*180/pi;
    
    %%% Get Sagittal (Psi)
    % Use to remove coronal and axial rotations
    [psi1, theta1, phi1] = getAngle(Infr,Supr);
%     [psi1, theta1, phi1]*180/pi
    
    [InfrRot,NewInfPlaneN(:,j),InfNewCenter(:,j),xgrid,ygrid,zgrid] = rotate(Infr, [-phi1,-theta1,0],displayflag, 'b.',displayflag);
    [psi1, theta1, phi1] = getAngle(InfrRot,Supr);
    IVDAngle = psi1*180/pi;
    
    
    if displayflag
        vectarrow([InfNewCenter(1,j) InfNewCenter(2,j) InfNewCenter(3,j)],...
            [a*NewInfPlaneN(1,j)+InfNewCenter(1,j)...
            a*NewInfPlaneN(2,j)+InfNewCenter(2,j)...
            a*NewInfPlaneN(3,j)+InfNewCenter(3,j)],'k');
    end
    
    
    
    DAnterior = sqrt(sum((Supr(3,:)-Infr(3,:)).^2));
    DCenter = sqrt(sum((Supr(13,:)-Infr(13,:)).^2));
    DPosterior = sqrt(sum((Supr(23,:)-Infr(23,:)).^2));
    
end