Ask Your Question
4

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


UPDATE

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

    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:

image description image description

which are for sure not similar :-(.

edit retag flag offensive close merge delete

Comments

1

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

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

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

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

1 answer

Sort by » oldest newest most voted
2

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)
    {
        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();
edit flag offensive delete link more
Login/Signup to Answer

Question Tools

3 followers

Stats

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

Seen: 867 times

Last updated: Sep 16 '15