# isotropic linear diffusion smoothing

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

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:

## UPDATE

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

, 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:

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);

waitKey(260);
}


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:

which are for sure not similar :-(.

edit retag close merge delete

1

Something like this :

int main (int argc,char **argv)
{
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;
}

( 2015-09-10 12:40:29 -0500 )edit

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

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

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

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

( 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

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

( 2015-09-11 15:28:44 -0500 )edit
1

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)

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

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

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

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

Sort by » oldest newest most voted

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)
{
min=d;
tmin=t;
iAfter=0;
x0.copyTo(xSol);
}
else
iAfter++;
cout << " t= "<<t<<" diffusivity "<<diffusivity<<" Error "<<d<<"\n";
Mat l;
Laplacian(x0,l,CV_32F,1);
x1 =x0 +(diffusivity*dt)*l;
x1.copyTo(x0);
t=t+dt;
}
xSol.convertTo(xc,CV_8U);
cout << " tMin= "<<tmin<<" diffusivity "<<diffusivity<<" Error "<<min<<" with dt ="<<dt<<" \n";
imshow("x1",xc);
xDst.convertTo(xc,CV_8U);
imshow("Blur x0",xc);
waitKey();

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

( 2019-10-18 19:12:35 -0500 )edit

Official site

GitHub

Wiki

Documentation