Ask Your Question

isotropic linear diffusion smoothing

asked 2015-09-10 06:08:09 -0500

theodore gravatar image

updated 2015-10-02 10:19:00 -0500

I want to apply the denoising filter I named in the title which is based on the following equations:

image description

where d is a scalar constant diffusivity parameter, I(x, y) is the initial noisy image, and u(x, y, t) is the image obtained after a diffusion time t lets say 5, 10 and 30. However, I am quite confused about which function to use and how, in order to achieve this in OpenCV. I have the feeling that it is quite simple but for some reason I am confused. Does anyone have an idea?

Here is a sample image:

image description


this is the result images after following @LBerger 's code, they are in time t = 0/t = 4 and d = 1:

image description image description

, is it expected to be like that?

I think that something is wrong, because I am trying also to compare it with the gaussian smoothing. And according to the following formula:

image description

where G√2t (x, y) is the Gaussian kernel. This proves that performing isotropic linear diffusion for a time t with d = 1 is exactly equivalent to performing Gaussian smoothing with a σ = √(2t)

and applying the gaussian filter with the following code:

void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, const float sigma, const int ksize_x = 0, const int ksize_y = 0)
    int ksize_x_ = ksize_x, ksize_y_ = ksize_y;

    // Compute an appropriate kernel size according to the specified sigma
    if (sigma > ksize_x || sigma > ksize_y || ksize_x == 0 || ksize_y == 0)
        ksize_x_ = (int)ceil(2.0f*(1.0f + (sigma - 0.8f) / (0.3f)));
        ksize_y_ = ksize_x_;

    // The kernel size must be and odd number
    if ((ksize_x_ % 2) == 0)
        ksize_x_ += 1;

    if ((ksize_y_ % 2) == 0)
        ksize_y_ += 1;

    // Perform the Gaussian Smoothing
    GaussianBlur(src, dst, Size(ksize_x_, ksize_y_), sigma, sigma, BORDER_DEFAULT);

    // show result
    std::ostringstream out;
    out << std::setprecision(1) << std::fixed << sigma;
    String title = "sigma: " + out.str();
    imshow(title, dst);
    imwrite("gaussian/" + title + ".png", dst);


and calling it with gaussian_2D_convolution(img, smoothed, sqrt(2*5)); the two results of gaussian smoothing and isotropic linear smoothing in time t = 5 are respectively:

image description image description

which are for sure not similar :-(.

edit retag flag offensive close merge delete



Something like this :

int main (int argc,char **argv)
Mat x = imread("C:/Users/Laurent.PC-LAURENT-VISI/Downloads/1441883251877551.png",CV_LOAD_IMAGE_GRAYSCALE);
Mat x0;

double t=0,dt=1;
double diffusivity=0.1;
Mat x1,xc;

while (t < 5)
    Mat l;
    imshow("Image t=" + to_string(t),xc);
waitKey(0) ;
return 0;
LBerger gravatar imageLBerger ( 2015-09-10 12:40:29 -0500 )edit

@LBerger thanks for the response ;-), would be able though to provide some comments on the code.

theodore gravatar imagetheodore ( 2015-09-10 14:01:05 -0500 )edit

Vectorial computation is very far but I remember well your equation is


You can use finite difference equation and it becomes :

u(t+1)-u(t) = k' laplacian(u(t)) and hence u(t+1) =u(t) + k' Laplacian(u(t))

k' is k by time unit. As you are not doing a real simulation of diffusivity you don't care.

Of course that's a brute force method I think you have read paper on subject. If you want to write something you can have a look here p 847

LBerger gravatar imageLBerger ( 2015-09-10 14:11:00 -0500 )edit

@LBerger thank you again for your help, running though your code I getting the result images that I have added in the update part of my initial question. They seem much more noisy than the original image, is that expected?

theodore gravatar imagetheodore ( 2015-09-11 12:47:18 -0500 )edit

Yes with your first example it's working. Is it? But noisiy image is not very far from answer. But in this example if you use t=50 instead of 50 it's wrong. My algorithm is not good because it does not converge. Criteria given page 847 is not respected at all

LBerger gravatar imageLBerger ( 2015-09-11 13:35:01 -0500 )edit

I do not think that the it works even in the first example. Moreover, I do not think that the image should be noisy as it appears, rather than smoother like the result from the gaussian filter.

theodore gravatar imagetheodore ( 2015-09-11 15:28:44 -0500 )edit

I don't understand your last comment. I have changed my program. Now I'm using diffusitvity coefficient time dependant. You can download here.

If you run program you can see that initial parameter are important (dt and diffusivity). I don't know method to reduce noise without this choice ( frequency, size of window...)

Again program given is not a reference I hope it's only a test to understand diffusion equation in image processing.

A way to improve this program is to consider D as function (x,y). Then laplacian is replace by your orignal expression and filter will become an adaptive filter (i have read some paper about this)

LBerger gravatar imageLBerger ( 2015-09-12 03:07:05 -0500 )edit

Just as a side note: If you are interested also in non-linear diffusion you can have a look at the KAZE-detector / desciptor, it is based on a non-linear diffusion scale-space if I remember correctly.

Guanta gravatar imageGuanta ( 2015-09-13 07:39:06 -0500 )edit

@Guanta thanks for the hint, that's actually my next step after I manage to accomplish this task :-). @LBerger thanks again for your hints and your help. Actually, I have opened a thread in the stackoverflow and there with the help of the user Miki who suggested me a solution based on your approach I managed to get somehow a good result.

theodore gravatar imagetheodore ( 2015-09-13 14:46:58 -0500 )edit

However using the following code:

 void isotropicLinearDiffusionSmoothing(const cv::Mat& src, cv::Mat& dst, const double diffusivity = 1)
        Mat1f img;
        src.convertTo(img, CV_32F);
        float lambda = 0.25;
    //    lambda = fmax(0.001, std::fmin(lambda, 0.25)); // something in [0, 0.25]
        int iter = 5;
        while (iter--)
            Mat1f lap;
            Laplacian(img, lap, CV_32F);
            img += lambda * diffusivity * lap;

        img.convertTo(dst, CV_8U);

        imshow("ILDS", dst);

I cannot get the same result with the gaussian smoothing. An diffusion time of 5 corresponds to a σ = √(2t) = 3.2 for the gaussian

theodore gravatar imagetheodore ( 2015-09-13 14:50:32 -0500 )edit

1 answer

Sort by » oldest newest most voted

answered 2015-09-16 08:54:16 -0500

LBerger gravatar image

After many comments _it could be an answer:

double t=0,dx2=1;
double diffusivity=1.49,dt=dx2/(4*diffusivity)*0.9,sigma=diffusivity;

int ksize_x_ = 0, ksize_y_ = 0;
ksize_x_ = (int)ceil(2.0f*(1.0f + (sigma - 0.8f) / (0.3f)));
ksize_y_ = ksize_x_;
if ((ksize_x_  % 2) == 0)
    ksize_x_ += 1;

if ((ksize_y_ % 2) == 0)
    ksize_y_ += 1;
Mat xDst;
Mat x1,xc;
GaussianBlur(x0, xDst, Size(ksize_x_, ksize_y_), sigma, sigma, BORDER_DEFAULT);
double min=1e40,tmin;
int iAfter=0;
Mat xSol;
while (iAfter<3)
    double d=norm(x0-xDst);
    if (d < min)
    cout << " t= "<<t<<" diffusivity "<<diffusivity<<" Error "<<d<<"\n";
    Mat l;
    x1 =x0 +(diffusivity*dt)*l;
cout << " tMin= "<<tmin<<" diffusivity "<<diffusivity<<" Error "<<min<<" with dt ="<<dt<<" \n";
imshow("Blur x0",xc);
edit flag offensive delete link more


Hi, I have been thinking about this problem too. and I did not understand the time steps in the isotropic smoothing iterations. I mean how could you compare the time in the isotropic smoothing way and Gaussian filter? (every iteration in your loops corresponds to which sigma in Gaussian filter?)

Dr. Magas gravatar imageDr. Magas ( 2019-10-18 19:12:35 -0500 )edit
Login/Signup to Answer

Question Tools



Asked: 2015-09-10 06:08:09 -0500

Seen: 968 times

Last updated: Sep 16 '15