function [ d ] = ngf_D( f, g, resolutions )

[row, column] = size(f);
n = row*column;

eps = 0.005;

[fx, fy] = gradient(f,resolutions(1), resolutions(2));
[gx, gy] = gradient(g,resolutions(1), resolutions(2));

norm_g_f = sqrt(fx.^2 + fy.^2 + eps.^2); 
ngf_fx = fx./norm_g_f;
ngf_fy = fy./norm_g_f;
 
norm_g_g = sqrt(gx.^2 + gy.^2 + eps.^2); 
ngf_gx = gx./norm_g_g;
ngf_gy = gy./norm_g_g;

d = ngf_fx.*ngf_gx + ngf_fy.*ngf_gy;

D = (1/(2*n))*d(:)'*d(:);

end

