1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
| #include <stdio.h>
#include <iostream>
#include <string>
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
using namespace cv;
using namespace std;
double computeDiffusion(double v, double lambda);
double flux(const Mat& img, int x, int y, double lambda);
void peronaMalik(Mat& out, const Mat& in, double dt, int nbrIter, double lambda);
int main(int argc , char** argv)
{
//Initialization of the data
Mat img = imread("../noisi.jpg", CV_LOAD_IMAGE_GRAYSCALE);
Mat out;
double dt = 0.1, lambda = 5.;
int nbrIter = 20;
img.convertTo(img, CV_64F, 1./255); //Convert the image source in CV_64F (i.e float with a range of [0;1])
peronaMalik(out, img, dt, nbrIter, lambda);
imshow("Non linear anysotropic diffusion", out);
waitKey(0);
return 0;
}
double computeDiffusion(double v, double lambda)
{
return 1. /( 1 + ( pow(v/lambda, 2.) ) );
}
double flux(const Mat& img, int i, int j, double lambda)
{
double current = img.at<double>(i, j);
int width = img.cols;
int height = img.rows;
int prevI = i - 1;
int nextI = i + 1;
int prevJ = j - 1;
int nextJ = j + 1;
//Condition at the border of the image
if(prevI < 0)
prevI = 0;
if(nextI >= height)
nextI = height - 1;
if(prevJ < 0)
prevJ = 0;
if(nextJ >= width)
nextJ = width - 1;
//Get the pixel in the neighborhoor of (i,j)
double pixelN = img.at<double>(prevI, j);
double pixelS = img.at<double>(nextI, j);
double pixelE = img.at<double>(i, nextJ);
double pixelW = img.at<double>(i, prevJ);
double pixelNW = img.at<double>(prevI, prevJ);
double pixelSW = img.at<double>(nextI, prevJ);
double pixelNE = img.at<double>(prevI, nextJ);
double pixelSE = img.at<double>(nextI, nextJ);
//Compute the conduction coefficients
double diffN = computeDiffusion(abs(pixelN - current), lambda);
double diffS = computeDiffusion(abs(pixelS - current), lambda);
double diffE = computeDiffusion(abs(pixelE - current), lambda);
double diffW = computeDiffusion(abs(pixelW - current), lambda);
double diffNW = computeDiffusion(abs(pixelNW - current), lambda);
double diffSW = computeDiffusion(abs(pixelSW - current), lambda);
double diffNE = computeDiffusion(abs(pixelNE - current), lambda);
double diffSE = computeDiffusion(abs(pixelSE - current), lambda);
//Compute the flux
double delta = (diffN * (pixelN - current))
+ (diffS * (pixelS - current))
+ (diffE * (pixelE - current))
+ (diffW * (pixelW - current))
+ 0.5 * ( (diffNW *(pixelNW - current))
+ (diffSW *(pixelSW - current))
+ (diffNE *(pixelNE - current))
+ (diffSE *(pixelSE - current)) );
return delta;
}
/* function which compute the anysotropic diffusion filter*/
void peronaMalik(Mat& out, const Mat& in, double dt, int nbrIter, double lambda)
{
Mat imgTemp;
in.copyTo(imgTemp);
in.copyTo(out);
for(int n = 0; n < nbrIter; ++n)
{
for(int i = 0; i < in.rows; ++i)
for(int j = 0; j < in.cols; ++j)
{
double value = (dt * flux(out, i, j, lambda)) + out.at<double>(i, j);
imgTemp.at<double>(i, j) = value;
}
imgTemp.copyTo(out);
}
return;
} |
Partager