OpenCV实践（2）- 矩阵的掩码操作

1 描述

在对数字图像进行处理时，我们一般都会在空间域（spatial domain）或者频域（frequency domain）中进行。所谓“空间域”，实际上指的是图像本身，在空间域上的操作常常是改变像素点的值，也就是经过一个映射（我们所做的变换，如滤波等），将原来的f(x,y)变换为新的g(x,y)。而“频域”，它的数学基础是法国学者傅里叶提出的傅里叶级数和随后发展起来的傅里叶变换。在这其中起到重要作用的，就是电子计算机的不断完善和快速傅里叶变换（FFT）算法的提出。这些使得傅里叶变换成为了一种有力的分析和变换工具。就像一列波，我们在时间上观察，每个时刻的幅值是一个时间的函数。而当我们变换角度，从频率域上去看，又会发现它是一系列正弦波的叠加，而这些正弦波的频率都会是某个基波频率的整数倍。可谓“横看成岭侧成峰”！

在空间域的操作主要可以分为两类：第一类是所谓的“图像强度变换”（Intensity Transform），另一类是所谓的“空间域图像滤波”（Spatial Filtering）。这两者的区别主要是处理方法的不同。前者对单个像素点进行操作，例如通过阈值函数实现图形的二值化，实现灰度平均等。而后者建立在邻域（neighborhood）的概念上，讲究的是利用一个矩阵核（Kernel）对一个小区域进行操作。今天这篇文章主要介绍的是后者，以及如何用OpenCV中的函数去实现。

2 测试demo描述

我们来看一下图像对比度增强方法的问题。主要地就是为图像中的每一个像素应用下面这个公式：



从公式中看到，把掩码矩阵的中心点放到你要计算的点上，然后加上用重叠的矩阵值相乘得到的像素值。用矩阵的表达形式，更容易查看。

现在让我们看看，分别用像素的基本访问方法和filter2D函数是怎样实现上面的运算的。

3 基本方法

下面是我们自己实现上面公式的代码：

void Sharpen(const Mat& myImage, Mat& Result)
{
CV_Assert(myImage.depth() == CV_8U);  // accept only uchar images

Result.create(myImage.size(), myImage.type());
const int nChannels = myImage.channels();

for(int j = 1; j < myImage.rows - 1; ++j)
{
const uchar* previous = myImage.ptr<uchar>(j - 1);
const uchar* current  = myImage.ptr<uchar>(j);
const uchar* next     = myImage.ptr<uchar>(j + 1);

uchar* output = Result.ptr<uchar>(j);

for(int i = nChannels; i < nChannels * (myImage.cols - 1); ++i)
{
*output++ = saturate_cast<uchar>(5 * current[i]
-current[i - nChannels] - current[i + nChannels] - previous[i] - next[i]);
}
}

Result.row(0).setTo(Scalar(0));
Result.row(Result.rows - 1).setTo(Scalar(0));
Result.col(0).setTo(Scalar(0));
Result.col(Result.cols - 1).setTo(Scalar(0));
}

首先，我们需要确定输入图像数据是unsigned char 格式。对此，我们使用CV_Assert函数，当其表达式是false时，会抛出错误。

CV_Assert(myImage.depth() == CV_8U);  // accept only uchar images

我们创建与数据矩阵相同大小和类型的矩阵，作为输出矩阵。当然了，图像在内存中的存储与通道数也是有关的。

Result.create(myImage.size(), myImage.type());
const int nChannels = myImage.channels();

我们使用[]操作符访问像素。因为我们需要同时访问多行，所以需要获取每一行的指针（前一行，后一行，和当前行）。我们还需要另一个指针保存我们的计算结果。然后使用[]操作符简化右值操作。

for(int j = 1; j < myImage.rows - 1; ++j)
{
const uchar *previous = myImage.ptr<uchar>(j - 1);
const uchar *current  = myImage.ptr<uchar>(j);
const uchar *next     = myImage.ptr<uchar>(j + 1);

uchar *output = Result.ptr<uchar>(j);

for(int i = nChannels; i < nChannels * (myImage.cols - 1); ++i)
{
*output++ = saturate_cast<uchar>(5 * current[i] - current[i - nChannels] - current[i + nChannels] - previous[i] - next[i]);
}
}

对于图像的边界值，上面的公式是无法实现计算的。所以，一个简单的方法就是直接对边界值进行赋值：

Result.row(0).setTo(Scalar(0));     // 顶行
Result.row(Result.rows - 1).setTo(Scalar(0)); // 低行
Result.col(0).setTo(Scalar(0));     // 最左列
Result.col(Result.cols - 1).setTo(Scalar(0)); // 最右列

4 filter2D函数

应用这样的过滤器和应用图像处理的其它模块是共同的。首先你需要定义一个Mat对象，持有一个这样的掩码：

Mat kern = (Mat_<char>(3,3) <<  0, -1,  0,
-1,  5, -1,
0, -1,  0);

然后，调用filter2D函数，指定输入和输出图像，和要使用的掩码：

filter2D(I, K, I.depth(), kern);

这个函数甚至有第五个参数，是可选的，用来指针kernel的中心点，第六个参数用来决定操作没有被定义到的区域做什么操作（这个区域就是边界）。用这个函数的优点就是：代码更短，也不复杂，由于内部许多优化技术的应用，导致比我们自己手写实现的代码拥有更快的执行速度。从我们的例子中可以看出，使用filter2D函数比我们自己实现的代码快6～7ms左右。

如图所示：



上图中，左图为处理前，右图为处理后。

运行效率比较：

（1） 自己写的代码: 0.00929756秒

（2）OpenCV的filter2D函数:0.00232559秒

5 源码

完整代码如下：

#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <iostream>

using namespace std;
using namespace cv;

static void help(char* progName)
{
cout << endl
<<  "This program shows how to filter images with mask: the write it yourself and the"
<< "filter2d way. " << endl
<<  "Usage:"                                                                        << endl
<< progName << " [image_name -- default lena.jpg] [G -- grayscale] "        << endl << endl;
}

void Sharpen(const Mat& myImage,Mat& Result);

int main( int argc, char* argv[])
{
help(argv[0]);
const char* filename = argc >=2 ? argv[1] : "lena.jpg";

Mat I, J, K;

if (argc >= 3 && !strcmp("G", argv[2]))
I = imread( filename, CV_LOAD_IMAGE_GRAYSCALE);
else
I = imread( filename, CV_LOAD_IMAGE_COLOR);

namedWindow("Input", WINDOW_AUTOSIZE);
namedWindow("Output", WINDOW_AUTOSIZE);

imshow("Input", I);
double t = (double)getTickCount();

Sharpen(I, J);

t = ((double)getTickCount() - t)/getTickFrequency();
cout << "Hand written function times passed in seconds: " << t << endl;

imshow("Output", J);
waitKey(0);

Mat kern = (Mat_<char>(3,3) <<  0, -1,  0,
-1,  5, -1,
0, -1,  0);
t = (double)getTickCount();
filter2D(I, K, I.depth(), kern );
t = ((double)getTickCount() - t)/getTickFrequency();
cout << "Built-in filter2D time passed in seconds:      " << t << endl;

imshow("Output", K);

waitKey(0);
return 0;
}
void Sharpen(const Mat& myImage,Mat& Result)
{
CV_Assert(myImage.depth() == CV_8U);  // accept only uchar images

const int nChannels = myImage.channels();
Result.create(myImage.size(),myImage.type());

for(int j = 1 ; j < myImage.rows-1; ++j)
{
const uchar* previous = myImage.ptr<uchar>(j - 1);
const uchar* current  = myImage.ptr<uchar>(j    );
const uchar* next     = myImage.ptr<uchar>(j + 1);

uchar* output = Result.ptr<uchar>(j);

for(int i= nChannels;i < nChannels*(myImage.cols-1); ++i)
{
*output++ = saturate_cast<uchar>(5*current[i]
-current[i-nChannels] - current[i+nChannels] - previous[i] - next[i]);
}
}

Result.row(0).setTo(Scalar(0));
Result.row(Result.rows-1).setTo(Scalar(0));
Result.col(0).setTo(Scalar(0));
Result.col(Result.cols-1).setTo(Scalar(0));
}