function [d_z] = delta_z(I, voisins, distances)

    I_max = max(I, [], "all");
    abs_grad_z = sqrt( (I_max ./ I).^2 - 1 );

    p1 = voisins(:, 1);
    p2 = voisins(:, 2);

    d_z = distances .* ( abs_grad_z(p1) + abs_grad_z(p2) ) / 2;

end