پیاده سازی فیلتر Anisotropic diffusion

Question

سلام می خواستم بدونم کسی راجع به پیاده سازی فیلتر (انتشار نا همسانگرد) یا Anisotropic diffusion می تونه کمکم کنه فرمولش خیلی پیچیده است

توضیحات در اینجا

 

Answer 1

سلام.

اگر با OpenCV آشنا هستید کد زیر اینکارو براتون انجم بده.

#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;
}