〔两行哥〕OpenCV4Android入门教程之算法系列(二):
卷积在信号处理领域有极其广泛的应用,也有严格的物理和数学定义。
OpenCV中对图像进行模糊操作,其背后的原理就是卷积运算,可是究竟卷积运算是什么,模糊的卷积算法又是如何实现的呢?本文将进行讨论。考虑到大部分读者的非专业性,本人将尽量不使用专业术语,而使用通俗易懂的“白话”进行讲述,若有任何疑问或建议,望读者斧正。
注:需要读者有一定数学基础,能有一定矩阵基础或矩阵的概念,同时对统计学正态分布有概念,这对理解卷积算法很有帮助。
一、概述
(一)一维卷积运算
卷积运算需要两个对象,一个是被卷积的对象(操作数矩阵),一个是决定卷积效果的对象(卷积数矩阵),两着共同决定卷积运算的结果。
假设P为一个一维操作数矩阵,其值为[1,2,3,4,5,6,7,8,9];
假设R为一个一维卷积数矩阵,其值为[-1,0,1]。
对操作数矩阵中的每个数值,计算它左右两个数值和卷积数矩阵中对应位置的数值的乘积,然后把结果相加,最终得到的值就作为操作数矩阵中每个数值的新值。
如图1所示,每行灰色方格为操作数矩阵,橙色方格为卷积数矩阵,橙色方格依次向右移动:
1.如图1第1行所示,卷积数矩阵[-1,0,1]覆盖了操作数矩阵中的[1,2,3],那么操作数矩阵中的第2个数,原本为2,经过卷积运算后的新值为-1 * 1 + 0 * 2 + 1 * 3 = 2。
2.如图1第2行所示,卷积数矩阵[-1,0,1]覆盖了操作数矩阵中的[2,3,4],那么操作数矩阵中的第3个数,原本为3,经过卷积运算后的新值为-1 * 2 + 0 * 3 + 1 * 4 = 2。
3.如图1第3行所示,卷积数矩阵[-1,0,1]覆盖了操作数矩阵中的[3,4,5],那么操作数矩阵中的第4个数,原本为4,经过卷积运算后的新值为-1 * 3 + 0 * 4 + 1 * 5 = 2。
4.以此类推,第7次运算,卷积数矩阵[-1,0,1]覆盖了操作数矩阵中的[7,8,9],那么操作数矩阵中的第8个数,原本为8,经过卷积运算后的新值为-1 * 7 + 0 * 8 + 1 * 8 = 2。
图1 一维卷积运算
综上,除去第1位数值1和第9位数值9没有运算,P经过R卷积后的矩阵为[1,2,2,2,2,2,2,2,9]。
1.我们可以抛弃首位1与末位9作为运算的结果,即[2,2,2,2,2,2,2]。
2.也可以保留首位1与末位9,同时对第1位数值1和第9位数值9补充运算,我们在P的首尾各加一个数值,首位增加1(与原P第1位数值相同),尾位增加9(与原P第9位数值相同),得到的新矩阵P'为[1,1,2,3,4,5,6,7,8,9,9],其经过上述卷积数矩阵运算后的结果变为[1,2,2,2,2,2,2,2,1]。
(二)二维(平面图像)卷积运算
以数字图像处理中基本的处理方法线性滤波为例,待处理图像可被看做一个大二维矩阵,图像的每个像素对应大矩阵的每个元素。假设图像的分辨率是1920 * 1080,那么对应的大矩阵的行数为1080,列数为1920。
用于滤波的是一个滤波器小矩阵(也叫卷积核)。滤波器小矩阵一般是个奇数方阵(行数和列数为相同奇数的矩阵),比如用于边缘检测的Sobel算子就是两个3 * 3的小矩阵。
进行滤波就是对大矩阵中的每个像素,计算它周围像素和滤波器矩阵对应位置元素的乘积,然后把结果相加,最终得到的值就作为该像素的新值,这样就完成了一个像素点的滤波,如图2所示。
举个例子,假设图像为3 * 3矩阵(仅有9个像素,如图3),滤波器矩阵如图4,那么图3中像素P5,经过一次滤波运算后的新值P5' = P1 * R1 + P2 * R2 + P3 * R3 + P4 * R4 + P5 * R5 + P6 * R6 + P7 * R7 + P8 * R8 + P9 * R9。需要留意的是,这里的运算与传统的矩阵乘法运算并不相同,仅仅是将两个矩阵对应位置元素的乘积相加。
图2 卷积运算
图3 待处理图像
图4 滤波器矩阵
那么如何对一个大矩阵进行卷积运算呢?与一维卷积运算类似,从大矩阵的左上角开始,依次右移至右上角,然后下移一个单位,从左边再次依次移至右端,以此类推,直至卷积矩阵移至右下角,如图5和图6所示。
图5 图像卷积1
图6 图像卷积2
对于边缘像素,采用与一维卷积运算类似的方法,补充与边缘像素一样的像素进行卷积运算即可,或直接丢弃边缘像素。
如上文所述,对图像大矩阵和滤波小矩阵对应位置元素相乘再求和的操作就叫卷积运算或协相关运算。两者很类似,两者唯一的差别就是卷积在计算前需要翻转卷积核,而协相关则不需要翻转。
注:这里的翻转,指的是绕卷积核中心旋转180度后得到的新矩阵,与矩阵运算中的转置不同。
二维卷积运算需要4个嵌套循环,所以它并不快,除非我们使用很小的卷积核。我们一般使用3x3或者5x5,而且也有其他的规则要求:
(1)滤波器的大小应该是奇数,这样它才有一个中心,例如3x3,5x5或者7x7。有中心,也就有了半径,例如5x5大小的核半径为2。
(2)滤波器矩阵所有的元素之和应该要等于1,这是为了保证滤波前后图像的亮度保持不变。当然这不是硬性要求。
(3)如果滤波器矩阵所有元素之和大于1,那么滤波后的图像就会比原图像更亮,反之,如果小于1,那么得到的图像就会变暗。如果和为0,图像不会变黑,但也会非常暗。
(4)若每个像素取值0~255,对于滤波后的结构,可能会出现负数或者大于255的数值。对这种情况,我们将其直接截断到0和255之间即可。对于负数,也可以取绝对值。
二、简单卷积操作
图像卷积操作最常见的效果就是图像的模糊与锐化。如下图7至图8的变换,即为图像的模糊。
图7 原图
图8 模糊图像
那么图像模糊的卷积算法是如何实现的呢?下面讲解几种图像卷积操作的算法。
注:为了保证处理后图片亮度一致,以下示例卷积核每个元素之和为1。
(一)均值模糊
均值模糊的卷积核矩阵如图9所示。
图9 均值模糊3*3卷积核矩阵
如上图所示,这是一个3 * 3的均值模糊矩阵,每个元素的值都为1/9。我们可以很明显地看出这个卷积核的数学特征,其本质就是将图像中某个像素点及周围8个像素点的值进行求和,然后求这9个像素点和的平均值,并将中心像素点的新值更改为平均值。同理,如果是5 * 5的卷积核,其每个元素应该为1/25。
让我们来想想,为什么图9这样的卷积核矩阵可以让图像看起来模糊了呢?可以看出,经过卷积运算后图像的各个像素点之间的“差异值”减小,即图像所有像素的方差减小,那么图像的观感就显得更加模糊。那如果我们使用 3 * 3 和 5 * 5两种大小的卷积核,哪种更加模糊呢?应该是第二种,各位读者可以算算两种卷积核处理后的差异。
其实,数值上,图像中每个像素点表现地更加平滑,图形上,图像中每个像素点表现地更加模糊,其本质就是“中间点”失去了细节。
(二)中值模糊
中值模糊与均值模糊不同,同样以3 * 3卷积核大小为例,其卷积算法可以表述为:将图像中某个像素点即周围8个像素点的值进行排序,然后求这9个像素点的中值,并将中心像素点的新值更改为中值。其效果与均值效果类似,均使得图像各个像素点之间的“差异值”减小,表现出模糊的视觉效果。
(三)边缘锐化
让我们看一个特殊的卷积核,如图10所示。
图10 示例卷积核
如果我们对一幅图像使用上述的卷积核,有什么变化呢?我们以图11的模糊图片为源图,图12展示了卷积处理后的结果,可以看到图12与图11相比,边缘轮廓更加清晰。
图11 模糊图像
图12 卷积处理后的图像
其实图10就是一个基本的3 * 3锐化滤波器矩阵,面积越大的卷积核锐化效果越明显。同样的,我们也可以通过调整卷积核来获得不同的锐化效果,如图13和图14所示。图10的锐化滤波器实际上是计算当前像素值和周围像素值的差别,然后将这个差别加到原来的位置上。另外,中间点的权值要比所有的权值和大于1,意味着这个像素是在保持原来值的基础上增加了这些差异值,同时因为卷积核的所有元素之和为1,也保证了其亮度不变。可能读者不是很理解,我们以图3的图像矩阵为例,使用图10的卷积核进行卷积运算:
P5 = (-1) * P1 + (-1) * P2 + (-1) * P3 + (-1) * P4 + 9 * P5 + (-1) * P6 + (-1) * P7 + (-1) * P8 + (-1) * P9 = P5 + (P5 - P1) + (P5 - P2) + (P5 - P3) + (P5 - P4) + (P5 - P6) + (P5 - P7) + (P5 - P8) + (P5 - P9)
(四)高斯模糊
与均值模糊不同的是,高斯模糊增加了对卷积核的权重考虑。均值模糊的卷积运算中,卷积核每个元素的权重相等,而高斯模糊,卷积核每个元素的权重参照了高斯分布。
我们先来看一下一维高斯分布(也叫做正态分布),如图13所示:
图13 一维高斯分布
正态分布中,越接近中心点,取值越大,越远离中心,取值越小。计算平均值的时候,我们只需要将“中心点”作为原点,其他点按照其在正态曲线上的位置来分配权重,就可以得到一个加权平均值。正态分布的密度函数也称为高斯函数,其一维表达式如图14所示:
图14 正态分布表达式
其中,μ是x的均值,σ是x的标准差。因为我们上文在计算平均值的时候,中心点就是原点,所以μ等于0,得到图15的函数:
图15 简化一维高斯函数
我们根据图15的表达式,将其推导至二维空间,便可用于二维图像的处理,如图16所示:
图16 简化二维高斯函数
在二维空间中,正态分布的图像如图17所示:
图17 二维正态分布
有了上述知识铺垫,让我们来计算一个3 * 3高斯模糊的卷积核。假设卷积核中心点坐标为(0,0),即原点为中心点,其他点坐标如图18所示:
图18 卷积核点坐标
再对图16中的σ的值取1,结合图18的坐标矩阵,得到权重矩阵,如图19所示:
图19 权重矩阵
接着对图19的权重矩阵进行归一化。这9个点的权重总和等于0.7794,如果只计算这9个点的加权平均,还必须让它们的权重之和等于1,因此上面9个值还要分别除以0.7794,得到最终的权重矩阵。之所以让权重矩阵所有元素的权重总值等于1,是为了保证图像亮度不变。否则的话,权重总值大于1会让图像偏亮,小于1会让图像偏暗。归一化后的权重矩阵如图20所示。
图20 归一化的权重矩阵
有了归一化后权重矩阵,就可以计算高斯模糊的值了。假设现有9个像素点,灰度值(0~255)如图21所示:
图21 待处理图像
每个点乘以图20权重矩阵对应位置的权重值,结果如图22所示:
图22 每个点乘以权重
将图22中9个元素值加起来,就是中心点的高斯模糊的值。
对所有点重复这个过程,就得到了高斯模糊后的图像。对于彩色图片来说,则需要对RGB三个通道分别做高斯模糊。对于边缘像素点,采用上文同样的处理方法。
三、OpenCV中高斯模糊算法优化
(一)参数指定
先来看一下OpenCV中高斯模糊Api的三个重载:
图23 高斯模糊Api
这里先以第一个参数最少的方法为例进行讲解:
-src:源图像;
-dst:输出图像;
-ksize:卷积核大小;
-sigmaX:高斯分布标准差,即图16中的σ。
这里需要注意的是,如果sigmaX的值为0或小于0,而ksize为有效值,OpenCV会根据ksize值自动计算sigmaX。如果sigmaX值大于0,OpenCV会忽略ksize的值,而根据sigmaX值自动计算ksize。两者的计算公式为:sigmaX = 0.3 * ((ksize - 1) * 0.5 - 1) + 0.8。
例如,如下三行代码执行效果是完全一样的:
Imgproc.GaussianBlur(mSrcMat, mTargetMat, new Size(5, 5), 0);
Imgproc.GaussianBlur(mSrcMat, mTargetMat, new Size(), 1.1);
Imgproc.GaussianBlur(mSrcMat, mTargetMat, new Size(55,55), 1.1);
(二)算法的代码实现
/**
* 构建权重矩阵,参照图16简化二维高斯函数
*/
public static double[][] getMatrix(int radius){
//根据半径(radius)来创建权重矩阵,size即为矩阵的边长
int size = 2 * radius + 1;
double[][] matrix = new double[size][size];
//通过图13可知,距离原点3σ以外的点,权重已经可以忽略不计。
//反推过来,以3σ长度作为radius,即σ取radius的1/3,取值比较适当,也简化了运算
double sigama = (double) radius / 3;
double sigamaDouble = 2 * sigama * sigama;
double sigamaPi = Math.PI * sigamaDouble;
int row = 0;
double sum = 0;
for(int i = -radius ; i <= radius ; i++){
int line = 0;
for(int j = -radius ; j <= radius ; j++){
double x = i * i;
double y = j * j;
matrix[row][line] = Math.exp(-(x + y)/sigamaDouble)/sigamaPi;
sum += matrix[row][line];
line++;
}
row++;
}
//归一处理,让权重总值等于1
//使用总值大于1的卷积核会让图像偏亮,小于1的卷积核会让图像偏暗
for(int i = 0 ; i < size ; i++){
for(int j = 0 ; j < size ; j++){
matrix[i][j] /= sum;
}
}
return matrix;
}
/**
* 对图像进行高斯模糊计算
*/
private void calcBlur() {
//获取权重矩阵
double[][] matrix = getMatrix(radius);
int width = srcBitmap.getWidth();
int height = srcBitmap.getHeight();
int[] currentPixels = new int[width * height];
srcBitmap.getPixels(currentPixels, 0, width, 0, 0, width, height)
for (int i = 0; i < width; i++) {
for (int j = 0; j < height; j++) {
int red = 0;
int green = 0;
int blue = 0;
int x = i - radius;
int y = j - radius;
//忽略边界值
if (x > 0 && y > 0 && (i + radius < width && j + radius < height)) {
for (int tempI = -radius; tempI <= radius; tempI++) {
for (int tempJ = -radius; tempJ <= radius; tempJ++) {
int color = currentPixels[(j + tempJ) * width + i + tempI];
red += (int) (Color.red(color) * matrix[tempI + radius][tempJ + radius]);
green += (int) (Color.green(color) * matrix[tempI + radius][tempJ + radius]);
blue += (int) (Color.blue(color) * matrix[tempI + radius][tempJ + radius]);
}
}
//将处理后的像素点放回原图像矩阵
currentPixels[j * width + i] = Color.rgb(red, green, blue);
}
}
}
}
这里的red、green、blue是通过Color.red()、Color.green()、Color.blue()方法转换而来,也可以直接通过二进制位运算得来:
red = (p & 0xff0000) >> 16;
green = (p & 0x00ff00) >> 8;
blue = (p & 0x0000ff);
如果不明白的话,请参阅:〔两行哥〕OpenCV4Android入门教程之API系列(二)
(三)算法优化:二维卷积转化为一维卷积
如果运行上文我们自己实现的高斯模糊逻辑,会发现效率非常低,执行时间很长。让我们分析一下calcBlur()方法的算法复杂度:
1.第一层循环for (int i = 0; i < width; i++) {}
2.第二层循环for (int j = 0; j < height; j++) {}
3.第三层循环for (int tempI = -radius; tempI <= radius; tempI++) {}
4.第四层循环for (int tempJ = -radius; tempJ <= radius; tempJ++) {}
那么上文高斯模糊逻辑的算法复杂度近似为O(width * height * (2 * radius) ^ 2)。上文的的处理逻辑是建立在二维的情况下进行的,而实际上高斯模糊也可以在二维图像上对两个独立的一维空间分别进行计算,即线性可分。也就是说,使用二维矩阵变换得到的效果可以通过在水平方向进行一维高斯矩阵变换叠加竖直方向进行的一维高斯矩阵变换得到。
/**
* 构建权重矩阵,参照图15简化一维高斯函数
*/
public static double[][] getMatrix(int radius){
//根据radius创建权重矩阵.
int size = 2 * radius + 1;
double[] matrix = new double[size];
double sigama = (double) radius / 3;
double sigamaDouble = 2 * sigama * sigama;
double sqlPi = Math.sqrt(2 * Math.PI);
double sigamaPi = sigama * sqlPi;
int row = 0;
double sum = 0;
for(int i = -radius ; i <= radius ; i++){
double x = i * i;
matrix[row] = Math.exp(-x/sigamaDouble)/sigamaPi;
sum += matrix[row];
row++;
}
//归一处理,让权重总值等于1
//使用总值大于1的卷积核会让图像偏亮,小于1的卷积核会让图像偏暗
for(int i = 0 ; i < size ; i++){
matrix[i] /= sum;
}
}
/**
* 对图像进行高斯模糊计算
*/
private void calcBlur() {
//获取权重矩阵
double[] matrix = getMatrix(radius);
int width = scaleBitmap.getWidth();
int height = scaleBitmap.getHeight();
int[] currentPixels = new int[width * height];
int red[] = new int[width * height];
int green[] = new int[width * height];
int blue[] = new int[width * height];
scaleBitmap.getPixels(currentPixels, 0, width, 0, 0, width, height);
//使用水平方向的一维高斯权重矩形进行滤波
for(int j = 0 ; j < height ; j++){
for(int i = 0 ; i < width ; i++){
int n = 0;
int x = i - radius;
int y = j - radius;
if(x >=0 && y >= 0 && (i+radius < width && j+radius < height)) {
for (int temp = -radius; temp <= radius; temp++) {
int point = temp + i;
int colorPoint = j * width + point;
int color = currentPixels[colorPoint];
red[colorPoint] += Color.red(color) * matrix[n];
green[colorPoint] += Color.green(color) * matrix[n];
blue[colorPoint] += Color.blue(color) * matrix[n];
n++;
}
}
}
}
//使用垂直方向的一维高斯权重矩阵进行滤波
for(int i = 0 ; i < width ; i++){
for(int j = 0 ; j < height ; j++){
int n = 0;
int r = 0 , b = 0 , g = 0;
int x = i - radius;
int y = j - radius;
if(x >=0 && y >= 0 && (i+radius < width && j+radius < height)) {
for (int temp = -radius; temp <= radius; temp++) {
int currentPoint = (j + temp) * width + i;
Log.e(TAG, "temp = " + temp + " i = " + i + " j : " + j + " currentPoint = " + currentPoint
);
r += red[currentPoint] * matrix[n];
g += green[currentPoint] * matrix[n];
b += blue[currentPoint] * matrix[n];
n++;
}
currentPixels[j*width + i] = Color.rgb(r, g, b);
}
}
}
}
优化后calcBlur()方法的算法复杂度:
1.使用水平方向的一维高斯权重矩形进行滤波
(1)第一层循环:for(int j = 0 ; j < height ; j++){}
(2)第二层循环:for(int i = 0 ; i < width ; i++){}
(3)第三层循环:for (int temp = -radius; temp <= radius; temp++) {}
2.使用垂直方向的一维高斯权重矩阵进行滤波
(1)第一层循环:for(int i = 0 ; i < width ; i++){}
(2)第二层循环:for(int j = 0 ; j < height ; j++){}
(3)第三层循环: for (int temp = -radius; temp <= radius; temp++) {}
其逻辑的算法复杂度近似为O(width * height * (2 * radius) * 2)。这里的radius总是大于1,即 (2 * radius) * 2总是小于(2 * radius) ^ 2,因此优化后的逻辑的算法复杂度比优化前的算法复杂度减少了radius倍。
实际上,OpenCV的源码也采用了一维的高斯权重矩阵进行了运算,而非二维权重矩阵。为了加快运算速度,OpenCV中将常用的3 * 3,5 * 5,7 * 7的权重矩阵作为常量写入了代码。
此时你也应该明白了,图23中为什么后两个OpenCV高斯模糊Api中有sigmaX和sigmaY两个参数了,因为OpenCV中已经将高斯卷积分解为横向一维高斯卷积和纵向一维卷积,那么用户就可以分别定义横向一维卷积时的sigmaX和纵向一维卷积时的sigmaY。默认情况下,如果没有指定sigmaY,则sigmaY=sigmaX。
那么如果sigmaX和sigmaY不一样会有什么样的图形效果呢?
图24 sigmaX大于sigmaY情况下高斯模糊
这是sigmaX远大于sigmaY时的高斯模糊图像,像不像目标在水平方向上高速运动时抓拍的图像?
本次关于图像卷积与模糊操作的讲解就到这里结束,其实图像卷积还有很多其他应用,我们下次再聊。
本人参考了以下文章,向原作者致谢:
Python下尝试实现图片的高斯模糊化
Opencv学习(1):高斯滤波
图像卷积与滤波的一些知识点
高斯模糊与图像卷积滤波一些知识点