图像增强的目的主要包括:①消除噪声,改善图像的视觉效果;②突出边缘,有利于识别和处理。空间域指的是图像平面本身,空间域中的图像处理方法直接对图像中的像素进行处理。表达式为:
g(x, y)=T[f(x, y)],
式中,f(x,y)是输入图像,g(x,y)是输出图像,T是在点(x,y)的一个邻域上定义的针对f的算子,这个算子可应用于单幅图像的像素或一组图像的像素。
假设,邻域为3×3大小的方形,算子T为“计算邻域内像素的平均灰度”。考虑图像中的一个任意位置,例如(100,150)。在输出图像中,该位置的结果g(100,150)是f(100,150)与其8个邻点的和除以9.然后将这个领域中心移到下一个相邻的位置,并重复上面的过程以生成输出图像g的下一个值。通常,这个过程始于输入图像的左上角,并且水平(垂直)扫描、一次一行(列)地逐个像素进行处理。当最小邻域的大小为1×1时,也就是g只依赖于点(x,y)处的f值,此时的T称为灰度变换函数。
反转变换函数得到灰度级在区间[0, L-1]内的反转图像的形式为
s=L - 1 - r
采用这种方式反转图像的灰度级,会得到类似于照片底片的结果。
对数变换的通式为
s = clog(1 + r)
可以拓展图像中的暗像素值,压缩高灰度级值。比如在显示傅里叶谱时,范围0-10^6或更高的频谱值很常见。计算机可以处理这么大的数字,但不能忠实再现范围如此宽的数值,会损失灰度细节,用对数变换可以很好地解决。
变换形式为
s = c rγ
式中,c和γ是正常数
将图像变换到频率域,在频率域中处理,然后对结果进行反变换,把结果带回空间域。伽马变换常用来校正显示器,人眼对于较暗(接近0)的亮度值比较敏感,对于较亮(接近1)的亮度值则不太敏感,伽马变换可以提高暗像素的存储精度,但会损失亮像素的存储精度。
设原图像灰度值在[0,M]之间,感兴趣目标的灰度范围为[a,b],目标灰度范围为[c,d],变换形式为
如果感兴趣的灰度范围是原图的灰度范围,目标范围是0-255,可以直接用
normalize(src, dst, iLow, iHigh, NORM_MINMAX);
或者自己实现:
void grayStrech(Mat& srcImage, Mat& dstImage) { double minValue, maxValue; minMaxLoc(srcImage, &minValue, &maxValue); dstImage = srcImage.clone(); uchar* p = srcImage.data; uchar* pd = dstImage.data; for (int i = 0; i < srcImage.rows * srcImage.cols; i++) { *pd = 255 / (maxValue - minValue) * (*p - minValue); p++; pd++; } }
效果:
灰度直方图反映了数字图像中每一灰度级与其出现频率间的关系,它能描述该图像的概貌。通过修改直方图的方法增强图像是一种实用而有效的处理技术。令r表示一幅图像f(x, y)的灰度,f的非归一化直方图定义为h(r) = n,n就是灰度为r的像素的数量。类似的,归一化的直方图定义为:
p(r) = h(r)/M/N
M和N是图像的行数和列数。
直方图均衡化是将原图像通过某种变换,得到一幅灰度直方图为均匀分布的新图像的方法。
假设r为原图像的灰度,s为修正后的图像灰度,值域为[0,L-1]。灰度映射变换为
s = T(r), r∈[0, L-1]
对输入图像中的给定灰度值r,它将产生一个输出灰度值s。
T满足下列条件:
① 在 0 ≤ r ≤ L-1内为单调递增函数,保证灰度级从黑到白的次序不变;
② 在 0 ≤ r ≤ L-1内,有 0 ≤ T ( r ) ≤ L-1 ,确保映射后的像素灰度在允许的范围内。
反变换关系为
r = T-1(s)
当变换函数为r的累积直方图函数时,可以满足上面的条件,使直方图均衡化:
s = T(r) = ∫pr(r)dr
对于离散的数字图像,用频率来代替概率,则变换函数为
s = T(r) = Σpr(r)
例1:假设有一幅图大小为64×64像素,灰度级为[0, 7]。其灰度分布如下表所示:
rk | nk | pr(rk) | S(rk) | 四舍五入到0-7之间 | sk | nsk | pk(s) |
r0=0 | 790 | 0.19 | 0.19 | 0.19*7=1.33 => 1 | s0=1 | 790 | 0.19 |
r1=1 | 1023 | 0.25 | 0.44 | 0.44*7=3.08 => 3 | s1=3 | 1023 | 0.25 |
r2=2 | 850 | 0.21 | 0.65 | 0.65*7=4.55 => 5 | s2=5 | 850 | 0.21 |
r3=3 | 656 | 0.16 | 0.81 | 0.81*7=5.67 => 6 | |||
r4=4 | 329 | 0.08 | 0.89 | 0.89*7=6.23 => 6 | s3=6 | 985 | 0.24 |
r5=5 | 245 | 0.06 | 0.95 | 0.95*7=6.65 => 7 | |||
r6=6 | 122 | 0.03 | 0.98 | 0.98*7=6.86 => 7 | |||
r7=7 | 81 | 0.02 | 1.00 | 1.00*7=7.00 => 7 | s4=7 | 448 | 0.11 |
void equalizeHist( InputArray src, OutputArray dst );
结果:
在某些情况下,并不一定需要具有均匀直方图的图像,有时需要具有特定的直方图的图像,以便能够增强图像中某些灰度级。用于生成具有规定直方图的图像的方法,称为直方图匹配或直方图规定化。直方图均衡化处理是直方图规定化的一个特例。
使用如下步骤可以得到一幅灰度级具有规定概率密度函数的图像:
1) 由输入图像做直方图均衡化后得到pr(r)
2) 由规定的概率密度函数得到函数G(z)
3) 计算反变换z=G-1(s)。这是从s到z的映射,s是均衡化后的灰度值,z是规定概率密度函数下的灰度值
例:假设有一幅图大小为64×64像素,灰度级为[0, 7]。规定的概率密度如下所示,
过程如下表所示,做完均衡化后如下表前三列所示,中间三列为规定的概率密度得到的。查找最接近sk的值,比如0.15*(L-1)=1.05,最接近1,所以对应性s0=> z3。直方图均衡化后的图像中,值为1的每个像素将映射到直方图规定化后的图像中值为3的图像。
sk | nsk | pk(s) | zk | pz(z) | G(z) | 找映射 | nk | pz(z) |
s0=1 | 790 | 0.19 | z0=0 | 0.00 | 0.00 | z0 | 0 | 0.00 |
s1=3 | 1023 | 0.25 | z1=1 | 0.00 | 0.00 | z1 | 0 | 0.00 |
s2=5 | 850 | 0.21 | z2=2 | 0.00 | 0.00 | z2 | 0 | 0.00 |
z3=3 | 0.15 | 1.05 | z3=> s0=1 | 790 | 0.19 | |||
s3=6 | 985 | 0.24 | z4=4 | 0.20 | 2.45 | z4=> s1=3 | 1023 | 0.25 |
z5=5 | 0.30 | 4.55 | z5=> s2=5 | 850 | 0.21 | |||
z6=6 | 0.20 | 5.95 | z6=> s3=6 | 985 | 0.24 | |||
s4=7 | 448 | 0.11 | z7=7 | 0.15 | 7.00 | z7=> s4=7 | 448 | 0.11 |
空间滤波通过把每个像素的值替换为该像素及其邻域的的函数值来修改图像。其中,平滑滤波器用于降低灰度的急剧过渡。由于随机噪声通常由灰度的急剧过渡组成,因此平滑的一个明显应用就是降噪。另一个应用是平滑因灰度级数量不足导致的图像中的伪轮廓。
盒式滤波器系数的值相同,因为所有的行和列都相同,所以秩为1,可分离。
opencv函数:
void boxFilter( InputArray src, OutputArray dst, int ddepth, Size ksize, Point anchor = Point(-1,-1), bool normalize = true, int borderType = BORDER_DEFAULT );
这里的α取值规则如下:
核取不同大小时的效果如下所示:
均值滤波是一种线性滤波器,将一个窗口区域中的像素计算平均值,然后将窗口中计算得到的均值设置为锚点上的像素值。
opencv函数:
blur( InputArray src, OutputArray dst, Size ksize, Point anchor = Point(-1,-1), int borderType = BORDER_DEFAULT );
blur和boxFilter都是方框型滤波器,简称方波,唯一不同的是,blur是boxFilter的归一化形式。前者是归一化方框型滤波器,后者是非归一化方框型滤波器。在OpenCV的官方文档中有这样的说明:The call blur(src, dst, ksize, anchor, borderType) is equivalent to boxFilter(src, dst, src.type(), ksize, anchor, true, borderType)。也就是说当boxFilter中参数normalize被设置为true时,这两个算子是完全相同的。
盒式滤波器过于简单,适合快速试验,产生视觉上能够接受的平滑结果,但是它在很多应用上面有局限性。应用中的所选的核一般是圆对称的,而高斯核已被证明是唯一可分离的圆对称核。因为可分离,所以高斯滤波器的计算优势可以和盒式滤波器媲美,而且具有很多适合于图像处理的性质。
中值滤波器是统计排序(非线性)滤波器,其响应基于滤波器所包含区域内的像素的排序。低通滤波器模糊了图像,降噪性能较差,在这种情况下,中值滤波明显远远优于低通滤波。
opencv函数:
void medianBlur( InputArray src, OutputArray dst, int ksize );
效果(左:原图(有椒盐噪声);中:高斯;右:中值)
锐化的作用是突出灰度中的过渡。空间模糊在空间域可以通过平均(平滑)一个邻域内的像素来实现。由于平均运算类似于积分运算,因此我们可以合理地认为锐化能够通过空间微分来实现。平滑通常称为低通滤波,类似的,锐化通常被称为高通滤波。
边缘在灰度上通常类似于斜坡过渡,这时一阶导数会产生较宽的边缘,因为斜坡的导数非零。另一方面,二阶导数会产生宽度为1像素并由0分隔的双边缘。由此我们得出的结论是,与一阶导数相比,二阶导数可增强更精细的细节,因此是锐化图像的一个理想特性。另外,与实现一阶导数相比,实现二阶导数所需的运算量更少。
图像中的一个3×3区域的灰度值如下图所示
z1 | z2 | z3 |
z4 | z5 | z6 |
z7 | z8 | z9 |
一阶导数的最简近似是gx=z8-z5和gy=z6-z5。Roberts在图像处理早期的开发中,使用交叉值提出了其他两个定义:gx=z9-z5和gy=z8-z6,因此罗伯特交叉梯度算子就是
如前所述,我们更喜欢使用奇数大小的核,因为它们关于唯一的(整数)中心空间对称。我们感兴趣的最小核是大小为3×3的核。使用中心为z5的3×3邻域时,gx和gy近似如下:
gx = (z7+2z8+z9) - (z1+2z2+z3) gy = (z3+2z6+z9) - (z1+2z4+z7)在图像中任意坐标的值为两个核与f(x,y)的卷积的平方相加后再开方:
M(x, y) = [gx2+gy2]1/2
下面的核称为sobel算子。在中心系数中使用权值2的原因是,强调中心的重要程度来实现某种平滑。所有核的系数都是0,因此他们对恒定灰度区域给出零响应。gx和gy的计算是线性运算,用卷积实现的,但是M(x, y)的计算涉及平方和平方根或绝对值,这是非线性运算。
例:
int main() { Mat src = imread("D:/Backup/桌面/3.png", 0); Mat gau,dst,mask; medianBlur(src, src, 3); Mat resultX, resultY, resultXY; //-------------Sobel边缘检测-------------- //X方向一阶边缘 Sobel(src, resultX, CV_16S, 1, 0, 3); convertScaleAbs(resultX, resultX); //Y方向一阶边缘 Sobel(src, resultY, CV_16S, 0, 1, 3); convertScaleAbs(resultY, resultY); //整幅图像的一阶边缘 resultXY = resultX + resultY; imshow("src", src); imshow("result", resultXY); waitKey(0); return 0; }Sobel算子
效果:
对于两个变量的函数,拉普拉斯定义为
在离散情况下:
利用卷积来实现拉普拉斯核:
0 | 1 | 0 |
1 | -4 | 1 |
0 | 1 | 0 |
void Laplacian( InputArray src, OutputArray dst, int ddepth, int ksize = 1, double scale = 1, double delta = 0, int borderType = BORDER_DEFAULT );
拉普拉斯增强核
0 | -1 | 0 |
-1 | 5 | -1 |
0 | -1 | 0 |
opencv实现:
void LaplacianEnhance(Mat& srcImage, Mat& dstImage) { for (size_t i = 1; i < srcImage.rows - 1; i++) { for (size_t j = 1; j < srcImage.cols - 1; j++){ dstImage.at<uchar>(i, j) = saturate_cast<uchar>(5 * srcImage.at<uchar>(i, j) - srcImage.at<uchar>(i - 1, j) - srcImage.at<uchar>(i, j - 1) - srcImage.at<uchar>(i + 1, j) - srcImage.at<uchar>(i, j + 1)); } } }
结果:(左:原图;中:拉普拉斯;右:拉普拉斯增强)
用原图减去拉普拉斯图像,也可以得到锐化后的图像,但是效果不太好。
还可以用cv::filter2D()进行卷积:
Mat kernel = (cv::Mat_<double>(3, 3) << 0,-1,0,-1,5,-1,0,-1,0); filter2D(src, dst, CV_8UC1, kernel);
从原图中减去一幅钝化(平滑后的)图像,是20世纪30年代以来印刷和出版社一直用来锐化图像的过程。这个过程称为钝化掩蔽,它由如下步骤组成:
1)模糊原图像
2)从原图减去模糊后的图像(产生的差称为模板)
3)将模板和原图相加
例:
int main() { Mat src = imread("D:/Backup/桌面/2.png", 0); Mat gau,dst,mask; GaussianBlur(src, gau, Size(3, 3), 5, 5); subtract(src, gau, mask); dst = mask * 3 + src; imshow("src", src); imshow("dst", dst); waitKey(0); return 0; }
效果:
假定原图像为f(x,y),经傅立叶变换为F(u,v)。频率域增强就是选择合适的滤波器H(u,v)对F(u,v)的频谱成分进行处理,然后经逆傅立叶变换得到增强的图像g(x,y)。
由于噪声主要集中在高频部分,为去除噪声改善图像质量,滤波器采用低通滤波器H(u,v) 来抑制高频成分,通过低频成分,然后再进行逆傅立叶变换获得滤波图像,就可达到平滑图像的目的。常用的频率域低通滤波器H(u,v) 有以下几种:
代码实现:
Mat idealLowFilter(Mat& img, int radius) { Mat result = Mat::zeros(img.size(), CV_32F); for (int i = 0; i < img.cols; i++) { for (int j = 0; j < img.rows; j++) { if (sqrt(pow(i - img.cols / 2, 2) + pow(j - img.rows / 2, 2)) < radius) result.at<float>(j, i) = 1; } } return result; }
代码实现:
Mat butterworthLowFilter(Mat& img, int radius, int n) { Mat result = Mat::zeros(img.size(), CV_32F); for (int i = 0; i < img.cols; i++) { for (int j = 0; j < img.rows; j++) { double d = sqrt(pow(i - img.cols / 2, 2) + pow(j - img.rows / 2, 2)); result.at<float>(j, i) = 1 / pow(1 + d / radius, 2 * n); } } return result; }
代码实现:
Mat gaussianLowFilter(Mat& img, int radius) { Mat result = Mat::zeros(img.size(), CV_32F); for (int i = 0; i < img.cols; i++) { for (int j = 0; j < img.rows; j++) { double d2 = pow(i - img.cols / 2, 2) + pow(j - img.rows / 2, 2); result.at<float>(j, i) = exp(-d2 / (2 * pow(radius, 2.0))); } } return result; }
例:
#include <iostream> #include <opencv.hpp> #include <opencv2\opencv.hpp> #include <opencv2/highgui/highgui.hpp> #include <opencv2/core/core.hpp> #include <opencv2/highgui/highgui_c.h> using namespace cv; using namespace std; //1. 理想低通滤波 Mat idealLowFilter(Mat& img, int radius) { Mat result = Mat::zeros(img.size(), CV_32F); for (int i = 0; i < img.cols; i++) { for (int j = 0; j < img.rows; j++) { if (sqrt(pow(i - img.cols / 2, 2) + pow(j - img.rows / 2, 2)) < radius) result.at<float>(j, i) = 1; } } return result; } //2. 巴特沃斯低通滤波器 Mat butterworthLowFilter(Mat& img, int radius, int n) { Mat result = Mat::zeros(img.size(), CV_32F); for (int i = 0; i < img.cols; i++) { for (int j = 0; j < img.rows; j++) { double d = sqrt(pow(i - img.cols / 2, 2) + pow(j - img.rows / 2, 2)); result.at<float>(j, i) = 1 / pow(1 + d / radius, 2 * n); } } return result; } // 3. 高斯低通滤波器 Mat gaussianLowFilter(Mat& img, int radius) { Mat result = Mat::zeros(img.size(), CV_32F); for (int i = 0; i < img.cols; i++) { for (int j = 0; j < img.rows; j++) { double d2 = pow(i - img.cols / 2, 2) + pow(j - img.rows / 2, 2); result.at<float>(j, i) = exp(-d2 / (2 * pow(radius, 2.0))); } } return result; } void FourierTransform(cv::Mat& image) { image.convertTo(image, CV_32F); //选取最适合做fft的宽和高 int m1 = getOptimalDFTSize(image.rows); int n1 = getOptimalDFTSize(image.cols); Mat padded; //填充 copyMakeBorder(image, padded, 0, m1 - image.rows, 0, n1 - image.cols, BORDER_CONSTANT, Scalar::all(0)); Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) }; Mat complexI; merge(planes, 2, complexI); //planes[0], planes[1]是实部和虚部 dft(complexI, complexI, DFT_SCALE | DFT_COMPLEX_OUTPUT); split(complexI, planes); //定义幅度谱和相位谱 Mat ph, mag, idft; phase(planes[0], planes[1], ph); magnitude(planes[0], planes[1], mag); //由实部planes[0]和虚部planes[1]得到幅度谱mag和相位谱ph //重新排列傅里叶图像中的象限,使得原点位于图像中心 int cx = mag.cols / 2; int cy = mag.rows / 2; Mat q0(mag, Rect(0, 0, cx, cy)); //左上角图像划定ROI区域 Mat q1(mag, Rect(cx, 0, cx, cy)); //右上角图像 Mat q2(mag, Rect(0, cy, cx, cy)); //左下角图像 Mat q3(mag, Rect(cx, cy, cx, cy)); //右下角图像 //变换左上角和右下角象限 Mat tmp; q0.copyTo(tmp); q3.copyTo(q0); tmp.copyTo(q3); //变换右上角和左下角象限 q1.copyTo(tmp); q2.copyTo(q1); tmp.copyTo(q2); imshow("mag", mag); int radius = 80; int n = 1; Mat filter = idealLowFilter(mag, radius); multiply(filter, mag, mag); imshow("mag2", mag); q0.copyTo(tmp); q3.copyTo(q0); tmp.copyTo(q3); q1.copyTo(tmp); q2.copyTo(q1); tmp.copyTo(q2); //傅里叶逆变换 polarToCart(mag, ph, planes[0], planes[1]); //由幅度谱mag和相位谱ph恢复实部planes[0]和虚部planes[1] merge(planes, 2, idft); dft(idft, idft, DFT_INVERSE | DFT_REAL_OUTPUT); image = idft(Rect(0, 0, image.cols & -2, image.rows & -2)); image.convertTo(image, CV_8U); } void main() { Mat img = imread("D:/Backup/桌面/4.png",0); imshow("src", img); FourierTransform(img); imshow("DFT img", img); waitKey(); system("pause"); return; }低通滤波
效果:
一般来说,自然图片的光照一般是均匀渐变的,所以i应该是低频分量,而不同物体对光的反射是具有突变的,所以r是高频分量。现在我们对两边取对数,并做Fourier变换,得到线性组合的频率域。
lnf(x,y)=lni(x,y)+lnr(x,y)
FFT(lnf(x,y))=FFT(lni(x,y))+FFT(lnr(x,y))
我们希望对低频能量进行压制,这样就降低了动态范围,而要对高频进行提高,这样就增强了图像的对比度:
opencv实现:
// 1. 高斯高通滤波器 Mat gaussianHighFilter(Mat& img, int radius, double gammaH, double gammaL, double c) { Mat result = Mat::zeros(img.size(), CV_32F); for (int i = 0; i < img.cols; i++) { for (int j = 0; j < img.rows; j++) { double d2 = pow(i - img.cols / 2, 2) + pow(j - img.rows / 2, 2); result.at<float>(j, i) = (gammaH - gammaL) * (1 - exp(-d2 / (2 * pow(radius, 2.0)))) + gammaL; } } return result; } Mat homoMorphicFilter(const Mat& img, double D0, double gammaH, double gammaL, double c) { Mat image; img.convertTo(image, CV_32F); log(image + 1, image); //选取最适合做fft的宽和高 int m1 = getOptimalDFTSize(image.rows); int n1 = getOptimalDFTSize(image.cols); Mat padded; //填充 copyMakeBorder(image, padded, 0, m1 - image.rows, 0, n1 - image.cols, BORDER_CONSTANT, Scalar::all(0)); Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) }; Mat complexI; merge(planes, 2, complexI); //planes[0], planes[1]是实部和虚部 dft(complexI, complexI, DFT_SCALE | DFT_COMPLEX_OUTPUT); split(complexI, planes); //定义幅度谱和相位谱 Mat ph, mag, idft; phase(planes[0], planes[1], ph); magnitude(planes[0], planes[1], mag); //由实部planes[0]和虚部planes[1]得到幅度谱mag和相位谱ph //重新排列傅里叶图像中的象限,使得原点位于图像中心 int cx = mag.cols / 2; int cy = mag.rows / 2; Mat q0(mag, Rect(0, 0, cx, cy)); //左上角图像划定ROI区域 Mat q1(mag, Rect(cx, 0, cx, cy)); //右上角图像 Mat q2(mag, Rect(0, cy, cx, cy)); //左下角图像 Mat q3(mag, Rect(cx, cy, cx, cy)); //右下角图像 //变换左上角和右下角象限 Mat tmp; q0.copyTo(tmp); q3.copyTo(q0); tmp.copyTo(q3); //变换右上角和左下角象限 q1.copyTo(tmp); q2.copyTo(q1); tmp.copyTo(q2); imshow("mag", mag); int radius = 30; Mat result = gaussianHighFilter(mag, radius, gammaH, gammaL, c); multiply(result, mag, mag); imshow("mag2", mag); q0.copyTo(tmp); q3.copyTo(q0); tmp.copyTo(q3); q1.copyTo(tmp); q2.copyTo(q1); tmp.copyTo(q2); //傅里叶逆变换 polarToCart(mag, ph, planes[0], planes[1]); //由幅度谱mag和相位谱ph恢复实部planes[0]和虚部planes[1] merge(planes, 2, idft); dft(idft, idft, DFT_INVERSE | DFT_REAL_OUTPUT); image = idft(Rect(0, 0, img.cols, img.rows)); exp(image, image); normalize(image, image, 0, 1, NORM_MINMAX); return image; }高斯同态滤波
效果:
参考:
1. 【计算机视觉】数字图像处理(四)—— 图像增强
2. 傅里叶变换及低通滤波再反变换
3. OpenCV C++(十一)----频率域滤波