Implementing Edge Preserving Diffusion (Anisotropic Diffusion)

I am trying to implement edge preserving diffusion. Recall the general diffusion equation to be:

$$\DeclareMathOperator{\Div}{div}\delta_t u = \Div(g \nabla u)$$

Where $$g$$ is the speed of the diffusion process and $$u$$ is the image. Now if we want an edge preserving diffusion, we can use a non linear diffusion speed as in Scale Space and Edge Detection Using Anisotropic Diffusion:

$$g(|\nabla u|) = \frac{1}{\sqrt{1+\frac{|\nabla u|}{\lambda^2}}}$$

To implement this smoothing effect, we can simply rewrite the diffusion equation as:

$$\delta_t u = \delta_x(g(|\nabla u|)\delta_xu) + \delta_y(g(|\nabla u|)\delta_yu)$$

And then perform gradient descent: $$u_{t+1} = u_t + \tau\delta_t u$$

This is implemented as follows:

lam = 0.1
im = smisc.imread('lena.png', mode='L')
for i in range(100):
ux = cv2.Sobel(im, cv2.CV_64F, 1, 0, ksize=5)  # x
uy = cv2.Sobel(im, cv2.CV_64F, 0, 1, ksize=5)  # y

g = 1 / (1 + (ux ** 2 + uy ** 2) / lam ** 2) ** 0.5
utx = cv2.Sobel(g * ux, cv2.CV_64F, 1, 0, ksize=5)  # x
uty = cv2.Sobel(g * uy, cv2.CV_64F, 0, 1, ksize=5)  # y

im = im + (utx + uty) * 0.01  # gradient descent
plt.imshow(im, cmap='gray', vmin=0, vmax=255)
plt.show()


This works well and smooths the image while preserving the edges. I decided that I would also try to implement this a different way, in particular I would solve the equations until the end and implement the smoothing. In particular I did: \begin{align} \delta_t u &= \Div(g \nabla u) \\ &= \delta_x(g(|\nabla u|)\delta_xu) + \delta_y(g(|\nabla u|)\delta_yu) \\ &= u_{xx}\big(1 + \frac{u_x^2 + u_y^2}{\lambda^2}\big)^{-\frac{1}{2}} - \frac{1}{\lambda^2}u_x^2u_{xx}\big(1 + \frac{u_x^2 + u_y^2}{\lambda^2}\big)^{-\frac{3}{2}} \\ &+ u_{yy}\big(1 + \frac{u_x^2 + u_y^2}{\lambda^2}\big)^{-\frac{1}{2}} - \frac{1}{\lambda^2}u_y^2u_{yy}\big(1 + \frac{u_x^2 + u_y^2}{\lambda^2}\big)^{-\frac{3}{2}} \\ \end{align}

Using gradient descent this can be implemented as:

lam = 0.1
im = smisc.imread('lena.png', mode='L')
for i in range(100):
ux = cv2.Sobel(im, cv2.CV_64F, 1, 0, ksize=5)  # x
uxx = cv2.Sobel(ux, cv2.CV_64F, 1, 0, ksize=5)  # x
uy = cv2.Sobel(im, cv2.CV_64F, 0, 1, ksize=5)  # y
uyy = cv2.Sobel(uy, cv2.CV_64F, 0, 1, ksize=5)  # y

utx = uxx * (1 + (ux**2 + uy**2) / lam**2) ** -0.5 - 1 / lam**2 * ux**2 * uxx * (1 + (ux**2 + uy**2) / lam**2) ** (-1.5)
uty = uyy * (1 + (ux**2 + uy**2) / lam**2) ** -0.5 - 1 / lam**2 * uy**2 * uyy * (1 + (ux**2 + uy**2) / lam**2) ** (-1.5)

im = im + (utx + uty) * 0.01 # gradient descent
plt.imshow(im, cmap='gray', vmin=0, vmax=255)
plt.show()


Interestingly this smooths as well, however, there appears to be some weird dot artifacts appearing in some locations, does anyone know why this happens and what is wrong with the second implementations? I show the outputs below:

Original image:

Implementation 1:

Implementation 2:

• What do you mean by "I would solve the equations until the end and implement the smoothing"? – Royi May 11 '19 at 7:34
• What I mean is that I plugged in what $g$ is and differentiated this with respect to $x$ and $y$. Doing so gave me $\delta_t u$ with respect to the first ($u_x, u_y$) and second derivatives ($u_xx, u_yy$) of the image $u$. – jlcv May 12 '19 at 3:03
• Have you validated your results using numeric finite differences? – Royi May 12 '19 at 15:33
• What do you mean by that @Royi? – jlcv May 14 '19 at 19:33
• Are you sure you derived the equations correctly? Have you validated them using Finite Differences (Numeric Derivative)? – Royi May 15 '19 at 3:28