图像锐化

1. 图像锐化

图像锐化也称边缘增强。锐化技术用于加强图像中的边界和细节信息。由于边界和细节信息对应频域中的高频部分，所以在频域中通常对图像进行高通滤波，在空间域则进行微分处理。

2. 空间域方法

在进行频率域处理前，先简单介绍两种空间域的图像锐化算子。

2.1. 梯度模算子

由数学知识，我们知道方向导数取得最大值的方向是梯度方向，这个最大值便是梯度的模。在边界和细节部分进行梯度运算，得到的值较大，而在较平缓的地方，会得到较小的梯度。为了运算简便，梯度模$G[f(i,j)]$可以进行如下三种近似运算：

其中$(m,n)$为$(i,j)$领域中的点，若$(m,n)$在角上要加权$1/\sqrt{2}$。
这里实现另外两种近似运算，分别是Roberts梯度算子和Sobel梯度算子。

1. Roberts梯度算子：
   $$G[f(i,j)]=max[|f(i,j)-f(i+1,j+1)|,|f(i+1,j)-f(i,j+1)|]$$
   实现代码：

#include<opencv2/opencv.hpp>
#include <iostream>
using namespace cv;
int main() {
	Mat src = imread("1.jpg",0);
	if (src.empty()) {
		std::cout << "加载图片失败" << std::endl;
		exit(0);
	}
	namedWindow("原图", WINDOW_NORMAL);
	imshow("原图", src);
	Mat dst(src.size(),src.type(),Scalar(0));
	int rows = src.rows;
	int cols = src.cols;
	for (int i = 0; i < rows-1; i++) {
		for (int j = 0; j < cols-1; j++) {
			uchar f1 = abs(src.at<uchar>(i, j) - src.at<uchar>(i + 1, j + 1));
			uchar f2 = abs(src.at<uchar>(i+1, j) - src.at<uchar>(i , j + 1));
			dst.at<uchar>(i, j) =  f1 > f2 ? f1 : f2;
		}
	}
	namedWindow("Roberts梯度增强", WINDOW_NORMAL);
	imshow("Roberts梯度增强", dst);
	waitKey(0);
	return 1;
}

Roberts梯度增强效果：

Roberts梯度增强对斜边有较好的增强性能，对水平方向和垂直方向的边界增强效果就不明显了。

2. Sobel梯度算子：

$$G_x=\left[\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right],G_y=\left[\begin{array}{ccc}1& 0& -1\\ 2& 0& -2\\ 1& 0& -1\end{array}\right]$$
则Sobel梯度为：
$$G[f(i,j)]=[G_x[f(i,j{)}^2]+G_y[f(i,j{)}^2{]}^{\frac{1}{2}}$$
或
$$G[f(i,j)]=|G_x[f(i,j)]|+|G_y[f(i,j)]|$$
这个算法的特点是对称一阶差分，中心加权，同时对图像具有平滑作用。
Sobel梯度算子在opencv中可通过Sobel函数实现。实现代码如下：

#include<opencv2/opencv.hpp>
#include <iostream>
using namespace cv;
int main() {
	Mat src = imread("1.jpg",0);
	if (src.empty()) {
		std::cout << "加载图片失败" << std::endl;
		exit(0);
	}
	namedWindow("原图", WINDOW_NORMAL);
	imshow("原图", src);
	Mat dst;
	Sobel(src, dst, src.depth(), 1,1, 3);
	namedWindow("Sobel梯度增强", WINDOW_NORMAL);
	imshow("Sobel梯度增强", dst);
	waitKey(0);
	return 1;
}

Sobel梯度增强效果：

3. 梯度模算子的其他情况
   对图像使用梯度模算子后，输出图像与输入图像关系为：
   $$g(i,j)=G[f(i,j)]$$
   式中，$G$为梯度模算子。
   有时候为了不破坏原图像较平缓的区域而又能有效加强边界和图像细节，可以适当选取门限$T$使得：
   
   其中，$L(i,j)$可根据实际需要选为某规定值或则为$G[f(i,j)]$。
   当只关心边界的灰度变化时可令：
   
   其中，$L_B$根据实际需要确定。
   若只对边界的位置感兴趣，则令：
   
   其中，$L_G,L_B$根据需要确定其大小。
   以最后一种情况为例，即：
   

这里取$L_B=0，L_G=255,T=20$。
代码如下：

#include<opencv2/opencv.hpp>
#include <iostream>
using namespace cv;
int main() {
	Mat src = imread("a.jpg",0);
	if (src.empty()) {
		std::cout << "加载图片失败" << std::endl;
		exit(0);
	}
	namedWindow("原图", WINDOW_NORMAL);
	imshow("原图", src);
	Mat dst;
	uchar T = 20;
	uchar LG = 255;
	uchar LB = 0;
	Sobel(src, dst, src.depth(), 1,1, 3);
	for (int i = 0; i < dst.rows; i++) {
		for (int j = 0; j < dst.cols; j++) {
			if (dst.at<uchar>(i, j) > T) {
				dst.at<uchar>(i, j) = LG;
			}
			else
			{
				dst.at<uchar>(i, j) = LB;
			}
		}
	}
	namedWindow("第一种情况效果", WINDOW_NORMAL);
	imshow("第一种情况效果", dst);
	waitKey(0);
	return 1;
}

效果：

2.2. Laplace算子

拉普拉斯(Laplace)是著名物理学家和数学家，laplace定义如下：

对于离散函数$f(i,j)$,laplace算子定义为：

因为

同理可得：

所以有

由此可得Laplace算子矩阵为
$$\left[\begin{array}{ccc}0& 1& 0\\ 1& -4& 1\\ 0& 1& 0\end{array}\right],$$
Laplace算子的实现可以通过opencv提供的Laplacian函数进行处理，也可使用opencv提供的filter2D进行掩码操作，两种方法的实现效果以及两种方法生成的图片的差值如下：

代码如下：

#include<opencv2/opencv.hpp>
#include <iostream>
using namespace cv;
int main() {
	Mat src = imread("a.jpg",0);
	if (src.empty()) {
		std::cout << "加载图片失败" << std::endl;
		exit(0);
	}
	namedWindow("原图", WINDOW_NORMAL);
	imshow("原图", src);
	Mat dst;
	Laplacian(src, dst, src.depth());
	namedWindow("laplace效果", WINDOW_NORMAL);
	imshow("laplace效果", dst);
	Mat kernel = (Mat_<char>(3, 3) << 0,1,0,1,-4,1,0,1,0);
	Mat dst1;
	filter2D(src, dst1, src.depth(), kernel);
	namedWindow("掩码实现效果", WINDOW_NORMAL);
	imshow("掩码实现效果", dst1);
	for (int i = 0; i < dst.rows; i++) {
		for (int j = 0; j < dst.cols;j++) {
			if (dst.at<uchar>(i, j) >= dst1.at<uchar>(i, j)) {
				dst.at<uchar>(i, j) -= dst1.at<uchar>(i, j);
			}
			else {
				dst.at<uchar>(i, j) = dst1.at<uchar>(i, j)- dst.at<uchar>(i, j);
			}
		}
	}
	namedWindow("两种方法生成的图片的差值", WINDOW_NORMAL);
	imshow("两种方法生成的图片的 差值", dst);
	waitKey(0);
	return 1;
}

由以下公式
还可以得到另一种算子矩阵：

$$\left[\begin{array}{ccc}0& -1& 0\\ -1& 5& -1\\ 0& -1& 0\end{array}\right],$$

若$f$表示边界模糊图像，由数学分析可知：当像素点$(i,j)$在边界处低灰度侧时，由于邻点灰度大于或等于$f(i,j)$，所以有：

产生下冲；
当像素点$(i,j)$在边界处高灰度侧时，由于邻点灰度小于或等于$f(i,j)$，所以有：

产生上冲
而在其他交均匀区域有：

不改变其内容。
由以上的分析可知，上冲和下冲会使图像边界变得更加陡峭，从而达到增强边界的效果，且不破化图像较平缓的区域。这个算子矩阵的实现可参考opencv+c++实现掩膜操作

3. 频域中图像锐化的方法

在频域中进行图像锐化主要是对图像进行高通滤波或高频提升滤波来实现。下面是几种常见滤波器。

3.1. 理想高通滤波器

理想高通滤波器(ideal high pass filter,IHPF)的传递函数为：

其中$D_0$为截止频率，$D(u，v)=(u^2+v^2{)}^{\frac{1}{2}}$。$H(u,v)$的三维示意图及二维剖面图如下：

算法的关键程序代码如下：

Mat ImgSharpening::idealHPF(Mat input, double D0)
{
	Mat src = input.clone();
	if (input.channels() != 1) {
		cvtColor(input, src, COLOR_BGR2GRAY);//转为灰度图
	}
	Mat srcDft = DFT::dft2(src);
	srcDft = DFT::shift(srcDft);//将零频率分量移到频谱中心
	int M = srcDft.rows;
	int N = srcDft.cols;
	int m1 = M / 2;
	int n1 = N / 2;
	int channels = srcDft.channels();
	for (int m = 0; m < M; m++) {
		double* current = srcDft.ptr<double>(m);
		for (int n = 0; n < N; n++) {
			if (sqrt((m - m1) * (m - m1) + (n - n1) * (n - n1)) <= D0) {
				current[n * channels] = 0.0;
				current[n * channels + 1] = 0.0;
			}
		}
	}
	srcDft = DFT::shift(srcDft);//还原
	Mat dst;
	idft(srcDft, dst, DFT_SCALE | DFT_REAL_OUTPUT);
	dst.convertTo(dst, CV_8U);
	return dst;
}

其中DFT::dft2和 DFT::shift定义如下：

Mat DFT::dft2(Mat input)
{
	if (input.channels() != 1) {
		std::cout << "只处理单通道的图像" << std::endl;
		exit(0);
	}
	Mat src = input.clone();
	//扩充图像尺寸为DFT的最佳尺寸
	Mat padded;
	int newRows = getOptimalDFTSize(src.rows);
	int newCols = getOptimalDFTSize(src.cols);
	copyMakeBorder(src, padded, newRows - src.rows, 0, newCols - src.cols, 0, BORDER_CONSTANT, Scalar(0));
	//将planes融合合并成一个多通道数组
	Mat plans[] = { Mat_<double>(padded),Mat::zeros(padded.size(),CV_64F) };
	Mat mergeArray;
	merge(plans, 2, mergeArray);
	//傅里叶变换
	dft(mergeArray, mergeArray);
	return mergeArray;
}


Mat DFT::shift(Mat input)
{
	Mat src = input.clone();
	src = src(Range(0, src.rows & -2), Range(0, src.cols & -2));
	int y0 = src.rows / 2;
	int x0 = src.cols / 2;
	Mat q0(src, Rect(0, 0, x0, y0));       //左上角
	Mat q1(src, Rect(x0, 0, x0, y0));      //右上角
	Mat q2(src, Rect(0, y0, x0, y0));      //左下角
	Mat q3(src, Rect(x0, y0, x0, y0));     //右下角
	swapArea(q0, q3);
	swapArea(q1, q2);
	return src;
}					

void DFT::swapArea(Mat& area1, Mat& area2)
{
	Mat temp;
	area1.copyTo(temp);
	area2.copyTo(area1);
	temp.copyTo(area2);
}

理想高通滤波效果：

3.2. 巴特沃斯高通滤波器

巴特沃斯高通滤波器(butterworth high pass filter,BHPF)的传递函数为：

其中$D_0$为截止频率，$D(u，v)=(u^2+v^2{)}^{\frac{1}{2}}$。$H(u,v)$的三维示意图及二维剖面图如下：

算法的关键程序代码如下：

Mat ImgSharpening::ButterworthHPF(Mat input, double D0, int n0)
{
	Mat src = input.clone();
	if (input.channels() != 1) {
		cvtColor(input, src, COLOR_BGR2GRAY);//转为灰度图
	}
	Mat srcDft = DFT::dft2(src);
	srcDft = DFT::shift(srcDft);//将零频率分量移到频谱中心
	int M = srcDft.rows;
	int N = srcDft.cols;
	int m1 = M / 2;
	int n1 = N / 2;
	int channels = srcDft.channels();
	for (int m = 0; m < M; m++) {
		double* current = srcDft.ptr<double>(m);
		for (int n = 0; n < N; n++) {
			double D = sqrt((m - m1) * (m - m1) + (n - n1) * (n - n1));
			double H = 1 / (1 + (sqrt(2) - 1) * pow(D0 / D, 2 * n0));
			current[n * channels] *= H;
			current[n * channels + 1] *= H;
		}
	}
	srcDft = DFT::shift(srcDft);//还原
	Mat dst;
	idft(srcDft, dst, DFT_SCALE | DFT_REAL_OUTPUT);
	dst.convertTo(dst, CV_8U);
	return dst;
}

巴特沃斯高通滤波效果：（取$D_0=20,n=2$）

3.3. 指数高通滤波器

指数高通滤波器的传递函数为：

其中$D_0$为截止频率，$D(u，v)=(u^2+v^2{)}^{\frac{1}{2}}$。$H(u,v)$的三维示意图及二维剖面图如下：（$D_0=10,n=2$）

算法的关键程序代码如下：

Mat ImgSharpening::expHPF(Mat input, double D0, int n0)
{
	Mat src = input.clone();
	if (input.channels() != 1) {
		cvtColor(input, src, COLOR_BGR2GRAY);//转为灰度图
	}
	Mat srcDft = DFT::dft2(src);
	srcDft = DFT::shift(srcDft);//将零频率分量移到频谱中心
	int M = srcDft.rows;
	int N = srcDft.cols;
	int m1 = M / 2;
	int n1 = N / 2;
	int channels = srcDft.channels();
	for (int m = 0; m < M; m++) {
		double* current = srcDft.ptr<double>(m);
		for (int n = 0; n < N; n++) {
			double D = sqrt((m - m1) * (m - m1) + (n - n1) * (n - n1));
			double H = exp((log(1 / sqrt(2))) * pow(-D0 / D, n0));
			current[n * channels] *= H;
			current[n * channels + 1] *= H;
		}
	}
	srcDft = DFT::shift(srcDft);//还原
	Mat dst;
	idft(srcDft, dst, DFT_SCALE | DFT_REAL_OUTPUT);
	dst.convertTo(dst, CV_8U);
	return dst;
}

指数高通滤波效果：（取$D_0=20,n=2$）

3.4. 梯形滤波器

梯形滤波器的传递函数为：

其中$D_0$为截止频率，$D_1<D_0$,$D_1$根据要求选择，$D(u，v)=(u^2+v^2{)}^{\frac{1}{2}}$。$H(u,v)$的三维示意图及二维剖面图如下：

算法的关键程序代码如下：

Mat ImgSharpening::trapezoidHPF(Mat input, double D0, double D1)
{//
	Mat src = input.clone();
	if (input.channels() != 1) {
		cvtColor(input, src, COLOR_BGR2GRAY);//转为灰度图
	}
	if (D0 <= D1) {
		std::cout << "输入参数有误，参数必须是D1<D0" << std::endl;
		exit(0);
	}
	Mat srcDft = DFT::dft2(src);
	srcDft = DFT::shift(srcDft);//将零频率分量移到频谱中心
	int M = srcDft.rows;
	int N = srcDft.cols;
	int m1 = M / 2;
	int n1 = N / 2;
	int channels = srcDft.channels();
	for (int m = 0; m < M; m++) {
		double* current = srcDft.ptr<double>(m);
		for (int n = 0; n < N; n++) {
			double D = sqrt((m - m1) * (m - m1) + (n - n1) * (n - n1));
			double H;
			if (D > D0) {
				H = 1;
			}
			else if (D < D1) {
				H = 0;
			}
			else {
				H = (D - D1) / (D0 - D1);
			}
			current[n * channels] *= H;
			current[n * channels + 1] *= H;
		}
	}
	srcDft = DFT::shift(srcDft);//还原
	Mat dst;
	idft(srcDft, dst, DFT_SCALE | DFT_REAL_OUTPUT);
	dst.convertTo(dst, CV_8U);
	return dst;
}

梯形高通滤波效果：（取$D_0=20,D_1=10$）

3.5. 高频加强滤波器

因为高通滤波器对噪声没有抑制作用，所以若简单的使用高通滤波器，图片可能受到噪声污染而变差。为了既加强图片细节又抑制噪声，可使用高频加强滤波器，这种滤波器实际上是由一个低通滤波器$H_l(u,v)$,一个高通滤波器$H_h(u,v)$和一个全通滤波器$H_a(u,v)$构成，传递函数为：
$$H(u,v)=H_l(u,v)H_h(u,v)+H_a(u,v)$$
其中
$$H_a(u,v)=a,(0<a<1)$$
低通滤波器$H_l(u,v)$的增益为$\sqrt{1-a}$，截止频率为$D_l$,高通滤波器$H_h(u,v)$的增益为$\sqrt{1-a}$，截止频率为$D_h$,
$H(u,v)$的径向截面图如下：

算法的关键程序代码如下：

Mat ImgSharpening::HFEF(Mat input, double Dh, double DL, double a)
{//高频加强滤波器（High frequency enhanced filter）
	//DL和Dh分别为低通，高通对应的截止频率，a为全通滤波器增益
	Mat src = input.clone();
	if (input.channels() != 1) {
		cvtColor(input, src, COLOR_BGR2GRAY);//转为灰度图
	}
	if (DL <= Dh) {
		std::cout << "输入参数有误，你的输入参数是DL <= Dh,退化成了全通滤波器" << std::endl;
		exit(0);
	}
	Mat srcDft = DFT::dft2(src);
	srcDft = DFT::shift(srcDft);//将零频率分量移到频谱中心
	int M = srcDft.rows;
	int N = srcDft.cols;
	int m1 = M / 2;
	int n1 = N / 2;
	int channels = srcDft.channels();
	for (int m = 0; m < M; m++) {
		double* current = srcDft.ptr<double>(m);
		for (int n = 0; n < N; n++) {
			double D = sqrt((m - m1) * (m - m1) + (n - n1) * (n - n1));
			double H;
			if (D > DL||D<Dh) {
				H = a;
			}
			else {
				H = 1;
			}
			current[n * channels] *= H;
			current[n * channels + 1] *= H;
		}
	}
	srcDft = DFT::shift(srcDft);//还原
	Mat dst;
	idft(srcDft, dst, DFT_SCALE | DFT_REAL_OUTPUT);
	dst.convertTo(dst, CV_8U);
	return dst;
}

图像增强的方法往往不是单一的，经常在一种方法使用完后再用另一种方法作后续处理，以达到最终目的。下面是一副图片经过高频加强滤波器处理后再进行直方图均化的效果。关于直方图均化可参考修正直方图增强对比度。