function [ df ] = Derivatives_NGF_v( f,g, resolutions )

eps = 0.005;

[rows,columns] = size(f);
t_y = -(rows-1)/2:(rows-1)/2; %defining vectors for mesh grid
t_x = -(columns-1)/2:(columns-1)/2; %defining vectors for mesh grid
MN = columns*rows;
[x_1,x_2] = meshgrid(t_x*resolutions(1),t_y*resolutions(2));

[Rx1, Rx2] = gradient(g, resolutions(1), resolutions(2));
[Tx1, Tx2] = gradient(f, resolutions(1), resolutions(2));
[Tx11, Tx12] = gradient(Tx1, resolutions(1), resolutions(2));
[Tx21, Tx22] = gradient(Tx2, resolutions(1), resolutions(2));

Tx1 = Tx1(:);
Tx2 = Tx2(:);
Rx1 = Rx1(:);
Rx2 = Rx2(:);

Tx11 = Tx11(:);
Tx12 = Tx12(:);
Tx21 = Tx21(:);
Tx22 = Tx22(:);

x_1 = x_1(:);
x_2 = x_2(:);

N0 = (Tx1.*Rx1 + Tx2.*Rx2);
N0 = [N0,N0]';

N1 = [-Rx1.*Tx11, -Rx1.*Tx12]';    
N2 = [-Rx2.*Tx21, -Rx2.*Tx22]'; 

Di = sqrt(Tx1.^2 + Tx2.^2 + eps^2).*sqrt(Rx1.^2 + Rx2.^2 + eps^2);
Di = [Di,Di]';

C1 = [Tx11.*Tx1, Tx12.*Tx1]';
C2 = [Tx21.*Tx2, Tx22.*Tx2]';

C = C1 + C2;

norm_grad_T = sqrt(Tx1.^2 + Tx2.^2 + eps^2);
norm_grad_T = [norm_grad_T,norm_grad_T]';
norm_grad_R = sqrt(Rx1.^2 + Rx2.^2 + eps^2);
norm_grad_R = [norm_grad_R,norm_grad_R]';

A = (N1 + N2)./Di;

B = C.*N0./(norm_grad_T.^3.*norm_grad_R);

d = ngf_D(f,g,resolutions);

d = [d(:),d(:)]'; 
% d = 1;
df = (1/(rows*columns))*d.*(A+B);

% df =(A+B);