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
Login/Signup to Answer

Question Tools



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

Seen: 911 times

Last updated: Sep 16 '15