function nablaI = Gradient(I); 
    [N1,N2]=size(I); nablaI = zeros(N1,N2,2);  
    S1 = [N1 1:N1-1]; S2 = [N2 1:N2-1]; 
    nablaI(:,:,1) = cat(2,I-I(:,S2)); nablaI(:,:,2) = cat(2, I-I(S1,:)); 
endfunction
