基本原理
各向異性擴散濾波主要是用來平滑圖像的,克服了高斯模糊的缺陷,各向異性擴散在平滑圖像時是保留圖像邊緣的,和雙邊濾波很像。各向異性擴散也叫P-M擴散,各向異性擴散(Anisotropic diffusion)的算法可以詳見論文:《Scale-space and edge detection using anisotropic diffusion》
通常我們有将圖像看作矩陣的,看作圖的,看作随機過程的。各向異性擴散濾波将圖像看作熱量場了。每個像素看作熱流,根據目前像素和周圍像素的關系,來确定是否要向周圍擴散。比如某個鄰域像素和目前像素差别較大,則代表這個鄰域像素很可能是個邊界,那麼目前像素就不向這個方向擴散了,這個邊界也就得到保留了。
論文有具體的推導公式,都是熱學上的,這裡隻介紹一下最終結論用到的公式。主要疊代方程如下:
I就是圖像了,因為是個疊代公式,是以有疊代次數t。
四個散度公式是在四個方向上對目前像素求偏導,news就是東南西北嘛,公式如下:
而cN/cS/cE/cW則代表四個方向上的導熱系數,邊界的導熱系數都是小的。公式如下:
其中g函數是
最後整個公式需要先前設定的參數主要有三個,疊代次數t,根據情況設定;導熱系數相關的k,取值越大越平滑,越不易保留邊緣;lambda同樣也是取值越大越平滑。
示例示範
各向異性濾波可以用于做人像美顔磨皮算法。工程代碼連結
#include <opencv2/opencv.hpp>
#include <iostream>
using namespace cv;
using namespace std;
float k = 15;
float lambda = 0.25;
int N = 20;
void BilateralFilter(Mat &image, Mat &result);
int main(int argc, char** argv)
{
Mat src = imread("D:\\TestData\\lena.jpg");
if (src.empty()) {
printf("could not load image...\n");
return -1;
}
namedWindow("input image", CV_WINDOW_AUTOSIZE);
imshow("input image", src);
vector<Mat> mv;
vector<Mat> results;
split(src, mv);
for (int n = 0; n < mv.size(); n++) {
Mat m = Mat::zeros(src.size(), CV_32FC1);
mv[n].convertTo(m, CV_32FC1);
results.push_back(m);
}
int w = src.cols;
int h = src.rows;
Mat copy = Mat::zeros(src.size(), CV_32FC1);
for (int i = 0; i < N; i++) {
BilateralFilter(results[0], copy);
copy.copyTo(results[0]);
BilateralFilter(results[1], copy);
copy.copyTo(results[1]);
BilateralFilter(results[2], copy);
copy.copyTo(results[2]);
}
Mat output;
normalize(results[0], results[0], 0, 255, NORM_MINMAX);
normalize(results[1], results[1], 0, 255, NORM_MINMAX);
normalize(results[2], results[2], 0, 255, NORM_MINMAX);
results[0].convertTo(mv[0], CV_8UC1);
results[1].convertTo(mv[1], CV_8UC1);
results[2].convertTo(mv[2], CV_8UC1);
Mat dst;
merge(mv, dst);
imshow("result", dst);
imwrite("./result.png", dst);
waitKey(0);
return 0;
}
void BilateralFilter(Mat &image, Mat &result) {
int width = image.cols;
int height = image.rows;
// 四鄰域梯度
float n = 0, s = 0, e = 0, w = 0;
// 四鄰域系數
float nc = 0, sc = 0, ec = 0, wc = 0;
float k2 = k * k;
for (int row = 1; row < height - 1; row++) {
for (int col = 1; col < width - 1; col++) {
// gradient
n = image.at<float>(row - 1, col) - image.at<float>(row, col);
s = image.at<float>(row + 1, col) - image.at<float>(row, col);
e = image.at<float>(row, col - 1) - image.at<float>(row, col);
w = image.at<float>(row, col + 1) - image.at<float>(row, col);
nc = exp(-n * n / k2);
sc = exp(-s * s / k2);
ec = exp(-e * e / k2);
wc = exp(-w * w / k2);
result.at<float>(row, col) = image.at<float>(row, col) + lambda * (n*nc + s * sc + e * ec + w * wc);
}
}
}