isotropic non-linear diffusion smoothing (Perona-Malik)

edit

it is my answer :

#include <opencv2/opencv.hpp> 
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp> 
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include "opencv2/core/ocl.hpp"

using namespace cv; 
using namespace std;


int main (int argc,char **argv)
{

// Refernce http://image.diku.dk/imagecanon/material/PeronaMalik1990.pdf (IEEE PAMI v12 n 7 1990)
Mat x = imread("f:/lib/opencv/samples/data/lena.jpg",CV_LOAD_IMAGE_GRAYSCALE);
Mat x0;
x.convertTo(x0, CV_32FC1);


double t=0;
double lambda=0.25; // Defined in equation (7)
double K=10; // defined after equation(13) in text
imshow("Original",x);

Mat x1,xc;

 while (t<20)
 {
Mat D; // defined just before equation (5) in text
Mat gradxX,gradyX; // Image Gradient t time 
Sobel(x0,gradxX,CV_32F,1,0,3);
Sobel(x0,gradyX,CV_32F,0,1,3);
D = Mat::zeros(x0.size(),CV_32F);
for (int i=0;i<x0.rows;i++)
    for (int j = 0; j < x0.cols; j++)
    {
        float gx = gradxX.at<float>(i, j), gy = gradyX.at<float>(i,j);
        float d;
        if (i==0 || i== x0.rows-1 || j==0 || j==x0.cols-1) // conduction coefficient set to 1 p633 after equation 13
            d=1;
        else
            d =1.0/(1.0+abs((gx*gx + gy*gy))/(K*K)); // expression of g(gradient(I))
            //d =-exp(-(gx*gx + gy*gy)/(K*K)); // expression of g(gradient(I))
        D.at<float>(i, j) = d;
   }
x1 = Mat::zeros(x0.size(),CV_32F);
for (int i = 1; i < x0.rows-1; i++)
{
    float *u1 = (float*)x1.ptr(i);
    u1++;
    for (int j = 1; j < x0.cols-1; j++,u1++)
    {
        // Value of I at (i+1,j),(i,j+1)...(i,j)
        float ip10=x0.at<float>(i+1, j),i0p1=x0.at<float>(i, j+1);
        float im10=x0.at<float>(i-1, j),i0m1=x0.at<float>(i, j-1),i00=x0.at<float>(i, j);
        // Value of D at at (i+1,j),(i,j+1)...(i,j)
        float cp10=D.at<float>(i+1, j),c0p1=D.at<float>(i, j+1);
        float cm10=D.at<float>(i-1, j),c0m1=D.at<float>(i, j-1),c00=D.at<float>(i, j);
        // Equation (7) p632
        *u1 = i00 + lambda/4*( (cp10+c00)*(ip10-i00) + (c0p1+c00)*(i0p1-i00) + (cm10+c00)*(im10-i00)+ (c0m1+c00)*(i0m1-i00));
        // equation (9)
   }
}
x1.copyTo(x0);
x0.convertTo(xc,CV_8U);
imshow("Perrony x0",xc);
cout << "*";
waitKey(10);
t=t+lambda;
}

waitKey();
return 0;
}