function [Dxc,Dyc,Dzc] = calculateNormals(dt)


sz = size(dt);

Dxc(sz(1),sz(2),sz(3))=0;
Dyc(sz(1),sz(2),sz(3))=0;
Dzc(sz(1),sz(2),sz(3))=0;

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

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

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

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

         %Dyc(x,y,z) = (dt(x,y+1,z)-dt(x,y-1,z))/2;
         %Dzc(x,y,z) = (dt(x,y,z+1)-dt(x,y,z-1))/2;
     end
  end
end