function [curvatures] = calculateCurvatures(dt)

sz = size(dt);

for x = 1:sz(2)
  for y = 1:sz(1)
     for z = 1:sz(3)

       if(x==1)
         Dy(x,y,z) = (dt(x+1,y,z)-dt(x,y,z));
       elseif(x==sz(2))
         Dy(x,y,z) = (dt(x,y,z)-dt(x-1,y,z));
       else
         Dy(x,y,z) = (dt(x+1,y,z)-dt(x-1,y,z))/2;
       end

       if(y==1)
         Dx(x,y,z) = (dt(x,y+1,z)-dt(x,y,z));
       elseif(y==sz(1))
         Dx(x,y,z) = (dt(x,y,z)-dt(x,y-1,z));
       else
         Dx(x,y,z) = (dt(x,y+1,z)-dt(x,y-1,z))/2;
       end

       if(z==1)
         Dz(x,y,z) = (dt(x,y,z+1)-dt(x,y,z));
       elseif(z==sz(3))
         Dz(x,y,z) = (dt(x,y,z)-dt(x,y,z-1));
       else
         Dz(x,y,z) = (dt(x,y,z+1)-dt(x,y,z-1))/2; 
       end

     end
  end
end

for x = 1:sz(2)
  for y = 1:sz(1)
     for z = 1:sz(3)

       if(x==1)
         Dyy(x,y,z) = (Dy(x+1,y,z)-Dy(x,y,z));
	 Dxy(x,y,z) = (Dx(x+1,y,z)-Dx(x,y,z));
	 Dzy(x,y,z) = (Dz(x+1,y,z)-Dz(x,y,z));
       elseif(x==sz(2))
         Dyy(x,y,z) = (Dy(x,y,z)-Dy(x-1,y,z));
	 Dxy(x,y,z) = (Dx(x,y,z)-Dx(x-1,y,z));
	 Dzy(x,y,z) = (Dz(x,y,z)-Dz(x-1,y,z));
       else
         Dyy(x,y,z) = (Dy(x+1,y,z)-Dy(x-1,y,z))/2;
         Dxy(x,y,z) = (Dx(x+1,y,z)-Dx(x-1,y,z))/2;
         Dzy(x,y,z) = (Dz(x+1,y,z)-Dz(x-1,y,z))/2;
       end

       if(y==1)
         Dxx(x,y,z) = (Dx(x,y+1,z)-Dx(x,y,z));
	 Dyx(x,y,z) = (Dy(x,y+1,z)-Dy(x,y,z));
	 Dzx(x,y,z) = (Dz(x,y+1,z)-Dz(x,y,z));
       elseif(y==sz(1))
         Dxx(x,y,z) = (Dx(x,y,z)-Dx(x,y-1,z));
	 Dyx(x,y,z) = (Dy(x,y,z)-Dy(x,y-1,z));
	 Dzx(x,y,z) = (Dz(x,y,z)-Dz(x,y-1,z));
       else
         Dxx(x,y,z) = (Dx(x,y+1,z)-Dx(x,y-1,z))/2;
         Dyx(x,y,z) = (Dy(x,y+1,z)-Dy(x,y-1,z))/2;
         Dzx(x,y,z) = (Dz(x,y+1,z)-Dz(x,y-1,z))/2;
       end

       if(z==1)
         Dzz(x,y,z) = (Dz(x,y,z+1)-Dz(x,y,z));
	 Dxz(x,y,z) = (Dx(x,y,z+1)-Dx(x,y,z));
	 Dyz(x,y,z) = (Dy(x,y,z+1)-Dy(x,y,z));
       elseif(z==sz(3))
         Dzz(x,y,z) = (Dz(x,y,z)-Dz(x,y,z-1));
	 Dxz(x,y,z) = (Dx(x,y,z)-Dx(x,y,z-1));
	 Dyz(x,y,z) = (Dy(x,y,z)-Dy(x,y,z-1));
       else
         Dzz(x,y,z) = (Dz(x,y,z+1)-Dz(x,y,z-1))/2;
         Dxz(x,y,z) = (Dx(x,y,z+1)-Dx(x,y,z-1))/2;
         Dyz(x,y,z) = (Dy(x,y,z+1)-Dy(x,y,z-1))/2;
       end

     end
  end
end

for x = 1:sz(2)
  for y = 1:sz(1)
     for z = 1:sz(3)
     
     H(x,y,z) = (Dxx(x,y,z)*(Dy(x,y,z)^2+Dz(x,y,z)^2)+Dyy(x,y,z)*(Dx(x,y,z)^2+Dz(x,y,z)^2)+Dzz(x,y,z)*(Dx(x,y,z)^2+Dy(x,y,z)^2)-2*(Dx(x,y,z)*Dy(x,y,z)*Dxy(x,y,z)+Dy(x,y,z)*Dz(x,y,z)*Dyz(x,y,z)+Dx(x,y,z)*Dz(x,y,z)*Dxz(x,y,z)))/(2*(Dx(x,y,z)^2+Dy(x,y,z)^2+Dz(x,y,z)^2)^1.5);
     
     K(x,y,z)=(Dx(x,y,z)^2*(Dyy(x,y,z)*Dzz(x,y,z)-Dyz(x,y,z)^2)+Dy(x,y,z)^2*(Dxx(x,y,z)*Dzz(x,y,z)-Dxz(x,y,z)^2)+Dz(x,y,z)^2*(Dxx(x,y,z)*Dyy(x,y,z)-Dxy(x,y,z)^2))/((Dx(x,y,z)^2+Dy(x,y,z)^2+Dz(x,y,z)^2)^2)-(2*(Dx(x,y,z)*Dy(x,y,z)*(Dxz(x,y,z)*Dyz(x,y,z)-Dxy(x,y,z)*Dzz(x,y,z))+Dy(x,y,z)*Dz(x,y,z)*(Dxy(x,y,z)*Dxz(x,y,z)-Dyz(x,y,z)*Dxx(x,y,z))+Dx(x,y,z)*Dz(x,y,z)*(Dxy(x,y,z)*Dyz(x,y,z)-Dxz(x,y,z)*Dyy(x,y,z))))/((Dx(x,y,z)^2+Dy(x,y,z)^2+Dz(x,y,z)^2)^2);
     
     N1 =-(Dy(x,y,z)*Dx(x,y,z)*Dyy(x,y,z)*Dz(x,y,z)-Dy(x,y,z)^2*Dx(x,y,z)*Dyz(x,y,z)-Dz(x,y,z)^3*Dxy(x,y,z)+Dz(x,y,z)^2*Dxz(x,y,z)*Dy(x,y,z)+Dz(x,y,z)^2*Dyz(x,y,z)*Dx(x,y,z)-Dz(x,y,z)*Dzz(x,y,z)*Dx(x,y,z)*Dy(x,y,z)-Dy(x,y,z)^2*Dxy(x,y,z)*Dz(x,y,z)+Dy(x,y,z)^3*Dxz(x,y,z));
     
     N2 =-(Dy(x,y,z)*Dx(x,y,z)*Dxy(x,y,z)*Dz(x,y,z)+Dy(x,y,z)^2*Dxz(x,y,z)*Dx(x,y,z)-Dx(x,y,z)^2*Dyz(x,y,z)*Dy(x,y,z)-Dz(x,y,z)^3*Dxx(x,y,z)+2*Dz(x,y,z)^2*Dxz(x,y,z)*Dx(x,y,z)-Dz(x,y,z)*Dzz(x,y,z)*Dx(x,y,z)^2-Dy(x,y,z)^2*Dxx(x,y,z)*Dz(x,y,z));
     
     D=Dz(x,y,z)*(Dx(x,y,z)^2+Dy(x,y,z)^2+Dz(x,y,z)^2)^1.5;

     pk1(x,y,z) = H(x,y,z) + sqrt(H(x,y,z)^2-K(x,y,z));
     pk2(x,y,z) = H(x,y,z) - sqrt(H(x,y,z)^2-K(x,y,z));
     
     pd1x(x,y,z) = (-Dz(x,y,z)*N1); 
     pd1y(x,y,z) = (-Dz(x,y,z)*N2+Dz(x,y,z)*pk1(x,y,z)*D); 
     pd1z(x,y,z) = (-Dx(x,y,z)*N1+Dy(x,y,z)*(N2-pk1(x,y,z)*D));
     
     pd2x(x,y,z) = (-Dz(x,y,z)*N1); 
     pd2y(x,y,z) = (-Dz(x,y,z)*N2+Dz(x,y,z)*pk2(x,y,z)*D);
     pd2z(x,y,z) = (-Dx(x,y,z)*N1+Dy(x,y,z)*(N2-pk2(x,y,z)*D));
     
     end
  end
end

curvatures.Dx = Dx;
curvatures.Dy = Dy;
curvatures.Dz = Dz;
curvatures.Dxy = Dxy;
curvatures.Dyz = Dyz;
curvatures.Dxz = Dxz;
curvatures.H = H;
curvatures.K = K;
curvatures.principal1=pk1;
curvatures.principal2=pk2;
curvatures.direction1x=pd1x;
curvatures.direction1y=pd1y;
curvatures.direction1z=pd1z;
curvatures.direction2x=pd2x;
curvatures.direction2y=pd2y;
curvatures.direction2z=pd2z;