isotropic linear diffusion smoothing

edit

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;
x.convertTo(x0,CV_32F);

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

while (t < 5)
{
    Mat l;
    Laplacian(x0,l,CV_32F);
    x1=x0+diffusivity*l;
    x1.convertTo(xc,CV_8U);
    x0=x1;
    imshow("Image t=" + to_string(t),xc);
    waitKey(10);
    t=t+dt;
}
waitKey(0) ;
return 0;
 }

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

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

du/dt=k.Laplacian(u)

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 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?

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

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.

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)

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 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.

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);
        waitKey();
    }

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