Basics of Edge Preserving Smoothing with Non-linear Diffusion

To remove the noise from the images, common filters (such as Gaussian Smoothing) may end up losing useful information by smoothing across the edges.

To prevent this phenomenon, edge-preserving smoothing methods developed by researchers. Today we will focus on a simple non-linear diffusion, also known as the Perona-Malik approach; one of the first non-linear diffusion techniques.

First, let's start from scratch.

What is Diffusion?

From Wikipedia we can get an intuitive definition;

Diffusion is the net movement of anything (for example, atoms, ions, molecules, energy) generally from a region of higher concentration to a region of lower concentration.

You may say that the definition is irrelevant to images since they are static and don’t contain any atom or molecule. But, images also have higher concentration and lower concentration “pixel” regions, and by diffusion, we cant help local brightness to diffuse neighboring pixels which should result in a smoother image.

Let's took at the equations we will use to describe diffusion.

Fick’s Law

Here, D will determine the speed of diffusion (i. e. diffusivity) and it can be either a constant number, a space-depending function g(x), an image-dependent function g(u), or a matrix. In the 3-rd case, the equation hence the diffusion is non-linear.

We will combine Fick’s Law with the continuity equation;

Continuity Equation

Note that for our initial purpose, we need to choose D so that we will have less diffusion when there is an edge!

So what we really want is D to be inversely proportional to the gradient of our image.

Perona-Malik suggested the following diffusion coefficient;

λ in here is the ‘contrast parameter’ because areas where |∇u| ≫ λ will not be affected that much during diffusion.

With the chosen diffusion coefficient wecan rewrite the flux(or j) as,

Also the continuity equation becomes;

Now since we obtained all we need, we can use Forward Euler Method to apply diffusion to our image!

I will be providing a simple python function below along with the results,I must note that the code written for gray scale images. Also, code uses another function to calculate sobel filters of arbitrary sizes; which will be a topic of another story.

Results

The full code with GUI is in below link

https://github.com/bertaye/Edge-Preserving-Nonlinear-Diffusion

def nonlinear_diffusion_filter(img, filter_size, n_iter, time_step=0.01):
    img = img.astype(np.float32)
    sobel_x = get_sobel_x(filter_size)
    sobel_y = get_sobel_y(filter_size)
    for i in range(n_iter):
        # Compute gradients in x and y directions
        grad_x = cv2.filter2D(src=img, ddepth=-1, kernel=sobel_x)
        grad_y = cv2.filter2D(src=img, ddepth=-1, kernel=sobel_y)
        kappa = 50; #or lambda
        grad_magnitude = np.sqrt(grad_x**2 + grad_y**2)
        diffusivity = 1 / (1 + (grad_magnitude/kappa)**2)
        
        img = img + time_step * (cv2.filter2D(grad_x * diffusivity, ddepth=-1, kernel=sobel_x) + cv2.filter2D(grad_y * diffusivity, ddepth=-1, kernel=sobel_y))
        img = np.clip(img, 0, 255)
        
    return img.astype(np.uint8)

def get_sobel_x(size_n):
    custom_sobel = np.ndarray((size_n,size_n))
    for i in range(size_n):
        for j in range(size_n):
            if j != size_n//2:
                custom_sobel[i,j] = (j-size_n//2)/((i-size_n//2)**2 + (j-size_n//2)**2)
            else:
                custom_sobel[i,j] = 0
    return custom_sobel

def get_sobel_y(size_n):
    custom_sobel = np.ndarray((size_n,size_n))
    for i in range(size_n):
        for j in range(size_n):
            if i != size_n//2:
                custom_sobel[i,j] = (i-size_n//2)/((i-size_n//2)**2 + (j-size_n//2)**2)
            else:
                custom_sobel[i,j] = 0
    return -1*custom_sobel