Revision history [back]

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