where N v = p × N n and N n is the number of neighbours for each source j, ∇D|v is the vth element of the gradient vector and Φv(u)=u2/[1+(uKv)2]. K v = α v × β v where α v depends on the distance between a source and its current neighbour and β v depends on the discrepancy regarding orientations of two sources considered. For small gradients the local cost is quadratic, thus producing areas with smooth spatial changes in intensity, whereas for higher gradients, the associated cost is finite: Φ v (u) ≈ Kv2, thus allowing the preservation of discontinuities. The estimator at the ith iteration is of the form