OpenCV 源码分析

OpenCV 版本:4.5.1

源码架构

modules 目录下是主要的功能代码, 在此目录下又根据不同的功能分为几个模块。如:

core 定义了基本的数据结构,包括 Mat 和被其它模块使用的的基础函数。

imgproc 图像处理模块,包括线性和非线性滤波、图像的几何变换、色彩空间转换、直方图等等。

video 视频分析模块,包括运动估计、背景差分、对象追踪算法。

objdetect 检测预定义好的类的对象和实例,如面部、眼睛、人物、汽车等。

highgui 一个易于使用的简单的 UI 能力接口。

videoio 一个易于使用的视频获取和解码的接口。

每个模块内部又有大致如下结构,以 imgproc 模块为例:

doc 模块相关文档

include 模块的头文件

misc 一些杂项,一般是其它语言如 Java 、Python、Objective-C 的一些相关文件。

perf 包含性能测试的一些文件。

src 基本是功能实现代码。

test 测试文件。

Screenshot_20210114_220830

core 模块

core 定义了基本的数据结构,包括 Mat 和被其它模块使用的的基础函数。

在 OpevCV 中一个很重要的类就是 Mat。Mat 类替代了先前用来表示图像的 C 语言中 IplImage 和 CvMat 数据结构,不需要再手动管理内存。在 OpenCV 中被用来表示一副图像,同时也可用作普通矩阵。

Mat 类的定义在 modules/core/include/opencv2/core/mat.hpp 头文件中,实现在modules/core/src/matrix.cpp 中。

Mat 类的定义:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
class CV_EXPORTS Mat
{
public:
// 一系列函数
...
/*! 包括一些位域
- Mat 的标识
- 数据是否连续
- 深度
- 通道数目
*/
int flags;
//! 矩阵的维度,>= 2
int dims;
//! 矩阵的行数和列数,超过 2 维取值为 -1
int rows, cols;
//! 指向数据的指针
uchar* data;

//! 定位 ROI(Region of interest) 和调整 ROI 的帮助域
const uchar* datastart;
const uchar* dataend;
const uchar* datalimit;

//! 自定义分配器
MatAllocator* allocator;
//! 标准分配器
static MatAllocator* getStdAllocator();
static MatAllocator* getDefaultAllocator();
static void setDefaultAllocator(MatAllocator* allocator);

//! 与 UMat 的交互
UMatData* u;

MatSize size;
MatStep step;
...
};

Mat 的构造函数之一:

1
2
3
4
5
6
Mat::Mat(int _rows, int _cols, int _type)
: flags(MAGIC_VAL), dims(0), rows(0), cols(0), data(0), datastart(0), dataend(0),
datalimit(0), allocator(0), u(0), size(&rows), step(0)
{
create(_rows, _cols, _type);
}

创建行数为 _rows,列数为 _cols,类型为 _type 的 Mat。

其函数内部通过 create() 函数来创建对象:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
void Mat::create(int d, const int* _sizes, int _type)
{
int i;
// 断言检查
CV_Assert(0 <= d && d <= CV_MAX_DIM && _sizes);
// 通过宏,对 _type 进行运算获得 Mat 类型
_type = CV_MAT_TYPE(_type);

if( data && (d == dims || (d == 1 && dims <= 2)) && _type == type() )
{
if( d == 2 && rows == _sizes[0] && cols == _sizes[1] )
return;
for( i = 0; i < d; i++ )
if( size[i] != _sizes[i] )
break;
if( i == d && (d > 1 || size[1] == 1))
return;
}

int _sizes_backup[CV_MAX_DIM]; // #5991
if (_sizes == (this->size.p))
{
for(i = 0; i < d; i++ )
_sizes_backup[i] = _sizes[i];
_sizes = _sizes_backup;
}

release();
if( d == 0 )
return;
// 计算 flags
flags = (_type & CV_MAT_TYPE_MASK) | MAGIC_VAL;
// 设置 size
setSize(*this, d, _sizes, 0, true);

if( total() > 0 )
{
// 设置分配器
MatAllocator *a = allocator, *a0 = getDefaultAllocator();
#ifdef HAVE_TGPU
if( !a || a == tegra::getAllocator() )
a = tegra::getAllocator(d, _sizes, _type);
#endif
if(!a)
a = a0;
try
{
// 分配空间, data 初始化为 0
u = a->allocate(dims, size, _type, 0, step.p, ACCESS_RW /* ignored */, USAGE_DEFAULT);
CV_Assert(u != 0);
}
catch (...)
{
if (a == a0)
throw;
u = a0->allocate(dims, size, _type, 0, step.p, ACCESS_RW /* ignored */, USAGE_DEFAULT);
CV_Assert(u != 0);
}
CV_Assert( step[dims-1] == (size_t)CV_ELEM_SIZE(flags) );
}

// 增加引用计数
addref();
finalizeHdr(*this);
}

读取、写入、显示图像

从文件中读取图像:

1
Mat img = imread(filename);

灰度图像:

1
Mat img = imread(filename, IMREAD_GRAYSCALE);

将图像写入文件:

1
imwrite(filename, img);

展示图像:

1
imshow( "Title", img );

imread() 方法在 modules/imgcodecs/src/loadsave.cpp 中:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
Mat imread( const String& filename, int flags )
{
CV_TRACE_FUNCTION();

/// 创建基本的图像容器
Mat img;

/// 加载数据
imread_( filename, flags, img );

/// 如果 EXIF 文件设置了方向 flag 就应用
if( !img.empty() && (flags & IMREAD_IGNORE_ORIENTATION) == 0 && flags != IMREAD_UNCHANGED )
{
ApplyExifOrientation(filename, img);
}

/// 返回图像数据的引用
return img;
}

imread_() 方法实现:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
static bool
imread_( const String& filename, int flags, Mat& mat )
{
/// 找到对应的解码器处理图像
ImageDecoder decoder;

#ifdef HAVE_GDAL
if(flags != IMREAD_UNCHANGED && (flags & IMREAD_LOAD_GDAL) == IMREAD_LOAD_GDAL ){
decoder = GdalDecoder().newDecoder();
}else{
#endif
// 根据文件信息找到解码器
decoder = findDecoder( filename );
#ifdef HAVE_GDAL
}
#endif

/// 如果没找到对应的解码器,返回 nothing
if( !decoder ){
return 0;
}

// 缩放分母
int scale_denom = 1;
if( flags > IMREAD_LOAD_GDAL )
{
if( flags & IMREAD_REDUCED_GRAYSCALE_2 )
scale_denom = 2;
else if( flags & IMREAD_REDUCED_GRAYSCALE_4 )
scale_denom = 4;
else if( flags & IMREAD_REDUCED_GRAYSCALE_8 )
scale_denom = 8;
}

/// 设置缩放
decoder->setScale( scale_denom );

/// 设置文件名
decoder->setSource( filename );

try
{
// 读取文件头确保成功
if( !decoder->readHeader() )
return 0;
}
catch (const cv::Exception& e)
{
std::cerr << "imread_('" << filename << "'): can't read header: " << e.what() << std::endl << std::flush;
return 0;
}
catch (...)
{
std::cerr << "imread_('" << filename << "'): can't read header: unknown exception" << std::endl << std::flush;
return 0;
}


// 建立所需图像文件的大小
Size size = validateInputImageSize(Size(decoder->width(), decoder->height()));

// 获得解码类型
int type = decoder->type();
if( (flags & IMREAD_LOAD_GDAL) != IMREAD_LOAD_GDAL && flags != IMREAD_UNCHANGED )
{
if( (flags & IMREAD_ANYDEPTH) == 0 )
type = CV_MAKETYPE(CV_8U, CV_MAT_CN(type));

if( (flags & IMREAD_COLOR) != 0 ||
((flags & IMREAD_ANYCOLOR) != 0 && CV_MAT_CN(type) > 1) )
type = CV_MAKETYPE(CV_MAT_DEPTH(type), 3);
else
type = CV_MAKETYPE(CV_MAT_DEPTH(type), 1);
}

// 根据上面得到的信息创建存放图像的 Mat
mat.create( size.height, size.width, type );

// 开始读取图像数据
bool success = false;
try
{
//将数据放入容器
if (decoder->readData(mat))
success = true;
}
catch (const cv::Exception& e)
{
std::cerr << "imread_('" << filename << "'): can't read data: " << e.what() << std::endl << std::flush;
}
catch (...)
{
std::cerr << "imread_('" << filename << "'): can't read data: unknown exception" << std::endl << std::flush;
}
if (!success)
{
mat.release();
return false;
}

// 如果设置了缩放就进行调整
if( decoder->setScale( scale_denom ) > 1 ) // if decoder is JpegDecoder then decoder->setScale always returns 1
{
resize( mat, mat, Size( size.width / scale_denom, size.height / scale_denom ), 0, 0, INTER_LINEAR_EXACT);
}

return true;
}

写操作类似,找到对应的编码器,设置好文件名,具体的写入交给编码器来做。

imshow() 位于modules/highgui/src/window.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
void cv::imshow( const String& winname, InputArray _img )
{
CV_TRACE_FUNCTION();
const Size size = _img.size();
#ifndef HAVE_OPENGL
CV_Assert(size.width>0 && size.height>0);
{
Mat img = _img.getMat();
CvMat c_img = cvMat(img);
cvShowImage(winname.c_str(), &c_img);
}
#else
const double useGl = getWindowProperty(winname, WND_PROP_OPENGL);
CV_Assert(size.width>0 && size.height>0);

if (useGl <= 0)
{
Mat img = _img.getMat();
CvMat c_img = cvMat(img);
cvShowImage(winname.c_str(), &c_img);
}
else
{
const double autoSize = getWindowProperty(winname, WND_PROP_AUTOSIZE);

if (autoSize > 0)
{
resizeWindow(winname, size.width, size.height);
}

setOpenGlContext(winname);

cv::ogl::Texture2D& tex = ownWndTexs[winname];

if (_img.kind() == _InputArray::CUDA_GPU_MAT)
{
cv::ogl::Buffer& buf = ownWndBufs[winname];
buf.copyFrom(_img);
buf.setAutoRelease(false);

tex.copyFrom(buf);
tex.setAutoRelease(false);
}
else
{
tex.copyFrom(_img);
}

tex.setAutoRelease(false);

setOpenGlDrawCallback(winname, glDrawTextureCallback, &tex);

updateWindow(winname);
}
#endif
}

下面还有部分代码,总的逻辑是如果定义了 OpenGL,就使用 OpenGL 进行绘制。如果没有定义 OpenGL,下面还可以使用 Qt、Win32UI、GTK 等方式绘图。

视频

画图

图像的代数运算

OpenCV 提供了 add() subtract() multiply() divide() 等图像的加减乘除等运算。在modules/core/include/opencv2/core.hpp 中定义,实现在modules/core/src/arithm.cpp中。

比如 add() 源码:

1
2
3
4
5
6
7
void cv::add( InputArray src1, InputArray src2, OutputArray dst,
InputArray mask, int dtype )
{
CV_INSTRUMENT_REGION();

arithm_op(src1, src2, dst, mask, dtype, getAddTab(), false, 0, OCL_OP_ADD );
}

根据传入的图像参数和运算类型,调用arithm_op方法进行运算。

几何变换

直方图

滤波处理

filter2D

可以使用 2 维卷积来对 2D 图像进行平滑和锐化操作。OpenCV 中使用 filter2D() 方法对一幅 2D 图像进行卷积操作。

源码位于modules/imgproc/src/filter.dispatch.cpp,函数定义如下:

1
2
3
void filter2D(InputArray _src, OutputArray _dst, int ddepth,
InputArray _kernel, Point anchor0,
double delta, int borderType)

内部调用了:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
void filter2D(int stype, int dtype, int kernel_type,
uchar * src_data, size_t src_step,
uchar * dst_data, size_t dst_step,
int width, int height,
int full_width, int full_height,
int offset_x, int offset_y,
uchar * kernel_data, size_t kernel_step,
int kernel_width, int kernel_height,
int anchor_x, int anchor_y,
double delta, int borderType,
bool isSubmatrix)
{
bool res;
res = replacementFilter2D(stype, dtype, kernel_type,
src_data, src_step,
dst_data, dst_step,
width, height,
full_width, full_height,
offset_x, offset_y,
kernel_data, kernel_step,
kernel_width, kernel_height,
anchor_x, anchor_y,
delta, borderType, isSubmatrix);
if (res)
return;

/*CV_IPP_RUN_FAST(ippFilter2D(stype, dtype, kernel_type,
src_data, src_step,
dst_data, dst_step,
width, height,
full_width, full_height,
offset_x, offset_y,
kernel_data, kernel_step,
kernel_width, kernel_height,
anchor_x, anchor_y,
delta, borderType, isSubmatrix))*/

res = dftFilter2D(stype, dtype, kernel_type,
src_data, src_step,
dst_data, dst_step,
width, height,
full_width, full_height,
offset_x, offset_y,
kernel_data, kernel_step,
kernel_width, kernel_height,
anchor_x, anchor_y,
delta, borderType);
if (res)
return;
ocvFilter2D(stype, dtype, kernel_type,
src_data, src_step,
dst_data, dst_step,
width, height,
full_width, full_height,
offset_x, offset_y,
kernel_data, kernel_step,
kernel_width, kernel_height,
anchor_x, anchor_y,
delta, borderType);
}

先尝试使用replacementFilter2D进行滤波处理。如果没有计算出结果,则使用dftFilter2D,基于 DFT (离散傅立叶变换)的滤波方式。如果还是没有计算出结果,则采用原始的方式进行计算。

blur

可以是使用blur()进行均值滤波。

比如 C++ 中:

1
blur( src, dst, Size( 3, 3 ), Point(-1,-1) );

源码位于modules/imgproc/src/box_filter.dispatch.cpp

1
2
3
4
5
6
7
void blur(InputArray src, OutputArray dst,
Size ksize, Point anchor, int borderType)
{
CV_INSTRUMENT_REGION();

boxFilter( src, dst, -1, ksize, anchor, true, borderType );
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
void boxFilter(InputArray _src, OutputArray _dst, int ddepth,
Size ksize, Point anchor,
bool normalize, int borderType)
{
CV_INSTRUMENT_REGION();

CV_Assert(!_src.empty());

CV_OCL_RUN(_dst.isUMat() &&
(borderType == BORDER_REPLICATE || borderType == BORDER_CONSTANT ||
borderType == BORDER_REFLECT || borderType == BORDER_REFLECT_101),
ocl_boxFilter3x3_8UC1(_src, _dst, ddepth, ksize, anchor, borderType, normalize))

CV_OCL_RUN(_dst.isUMat(), ocl_boxFilter(_src, _dst, ddepth, ksize, anchor, borderType, normalize))

Mat src = _src.getMat();
int stype = src.type(), sdepth = CV_MAT_DEPTH(stype), cn = CV_MAT_CN(stype);
if( ddepth < 0 )
ddepth = sdepth;
_dst.create( src.size(), CV_MAKETYPE(ddepth, cn) );
Mat dst = _dst.getMat();
if( borderType != BORDER_CONSTANT && normalize && (borderType & BORDER_ISOLATED) != 0 )
{
if( src.rows == 1 )
ksize.height = 1;
if( src.cols == 1 )
ksize.width = 1;
}

Point ofs;
Size wsz(src.cols, src.rows);
if(!(borderType&BORDER_ISOLATED))
src.locateROI( wsz, ofs );

CALL_HAL(boxFilter, cv_hal_boxFilter, src.ptr(), src.step, dst.ptr(), dst.step, src.cols, src.rows, sdepth, ddepth, cn,
ofs.x, ofs.y, wsz.width - src.cols - ofs.x, wsz.height - src.rows - ofs.y, ksize.width, ksize.height,
anchor.x, anchor.y, normalize, borderType&~BORDER_ISOLATED);

CV_OVX_RUN(true,
openvx_boxfilter(src, dst, ddepth, ksize, anchor, normalize, borderType))

//CV_IPP_RUN_FAST(ipp_boxfilter(src, dst, ksize, anchor, normalize, borderType));

borderType = (borderType&~BORDER_ISOLATED);

Ptr<FilterEngine> f = createBoxFilter( src.type(), dst.type(),
ksize, anchor, normalize, borderType );

f->apply( src, dst, wsz, ofs );
}

medianBlur

中值滤波使用中位数来代替模板中原点对应位置的值。

1
medianBlur ( src, dst, 3 );

源码位于modules/imgproc/src/median_blur.dispatch.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
void medianBlur( InputArray _src0, OutputArray _dst, int ksize )
{
CV_INSTRUMENT_REGION();

CV_Assert(!_src0.empty());

CV_Assert( (ksize % 2 == 1) && (_src0.dims() <= 2 ));

if( ksize <= 1 || _src0.empty() )
{
_src0.copyTo(_dst);
return;
}

// 尝试使用 OpenCL
CV_OCL_RUN(_dst.isUMat(),
ocl_medianFilter(_src0,_dst, ksize))

Mat src0 = _src0.getMat();
_dst.create( src0.size(), src0.type() );
Mat dst = _dst.getMat();

// 尝试使用 cv_hal_medianBlur,目前这是个空方法
CALL_HAL(medianBlur, cv_hal_medianBlur, src0.data, src0.step, dst.data, dst.step, src0.cols, src0.rows, src0.depth(),
src0.channels(), ksize);

// 尝试使用 OpenVX
CV_OVX_RUN(true,
openvx_medianFilter(_src0, _dst, ksize))

//CV_IPP_RUN_FAST(ipp_medianFilter(src0, dst, ksize));

CV_CPU_DISPATCH(medianBlur, (src0, dst, ksize),
CV_CPU_DISPATCH_MODES_ALL);
}

最后的 OpenCV 中内置的中值滤波的实现方法:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
void medianBlur(const Mat& src0, /*const*/ Mat& dst, int ksize)
{
CV_INSTRUMENT_REGION();

bool useSortNet = ksize == 3 || (ksize == 5);

Mat src;
if( useSortNet )
{
if( dst.data != src0.data )
src = src0;
else
src0.copyTo(src);

if( src.depth() == CV_8U )
medianBlur_SortNet<MinMax8u, MinMaxVec8u>( src, dst, ksize );
else if( src.depth() == CV_16U )
medianBlur_SortNet<MinMax16u, MinMaxVec16u>( src, dst, ksize );
else if( src.depth() == CV_16S )
medianBlur_SortNet<MinMax16s, MinMaxVec16s>( src, dst, ksize );
else if( src.depth() == CV_32F )
medianBlur_SortNet<MinMax32f, MinMaxVec32f>( src, dst, ksize );
else
CV_Error(CV_StsUnsupportedFormat, "");

return;
}
else
{
// TODO AVX guard (external call)
cv::copyMakeBorder( src0, src, 0, 0, ksize/2, ksize/2, BORDER_REPLICATE|BORDER_ISOLATED);

int cn = src0.channels();
CV_Assert( src.depth() == CV_8U && (cn == 1 || cn == 3 || cn == 4) );

double img_size_mp = (double)(src0.total())/(1 << 20);
if( ksize <= 3 + (img_size_mp < 1 ? 12 : img_size_mp < 4 ? 6 : 2)*
(CV_SIMD ? 1 : 3))
medianBlur_8u_Om( src, dst, ksize );
else
medianBlur_8u_O1( src, dst, ksize );
}
}

可以看到共有 3 个出口,如果是核的大小是 3 或 5,则使用 SortNet。如果是较小的核使用 Om 的算法,较大的核使用 O1 的算法。

再看下medianBlur_SortNet()方法:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
// Op 用于比较两个值的大小,VecOp 用于比较两个向量的大小
template<class Op, class VecOp>
static void
medianBlur_SortNet( const Mat& _src, Mat& _dst, int m )
{
CV_INSTRUMENT_REGION();

typedef typename Op::value_type T;
typedef typename Op::arg_type WT;
typedef typename VecOp::arg_type VT;
#if CV_SIMD_WIDTH > 16
typedef typename VecOp::warg_type WVT;
#endif

const T* src = _src.ptr<T>();
T* dst = _dst.ptr<T>();
int sstep = (int)(_src.step/sizeof(T));
int dstep = (int)(_dst.step/sizeof(T));
Size size = _dst.size();
int i, j, k, cn = _src.channels();
Op op;
VecOp vop;

if( m == 3 )
{
if( size.width == 1 || size.height == 1 )
{
// 处理特殊情况
return;
}

size.width *= cn;
// 行遍历
for( i = 0; i < size.height; i++, dst += dstep )
{
const T* row0 = src + std::max(i - 1, 0)*sstep;
const T* row1 = src + i*sstep;
const T* row2 = src + std::min(i + 1, size.height-1)*sstep;
int limit = cn;

// 列遍历
for(j = 0;; )
{
// 遍历 cn 个像素点,对每个通道都进行比较得到中值
for( ; j < limit; j++ )
{
int j0 = j >= cn ? j - cn : j;
int j2 = j < size.width - cn ? j + cn : j;
WT p0 = row0[j0], p1 = row0[j], p2 = row0[j2];
WT p3 = row1[j0], p4 = row1[j], p5 = row1[j2];
WT p6 = row2[j0], p7 = row2[j], p8 = row2[j2];

// 当 p1 < p2 时交换两个数字
op(p1, p2); op(p4, p5); op(p7, p8); op(p0, p1);
op(p3, p4); op(p6, p7); op(p1, p2); op(p4, p5);
op(p7, p8); op(p0, p3); op(p5, p8); op(p4, p7);
op(p3, p6); op(p1, p4); op(p2, p5); op(p4, p7);
op(p4, p2); op(p6, p4); op(p4, p2);
dst[j] = (T)p4;
}

if( limit == size.width )
break;

for( ; j <= size.width - VecOp::SIZE - cn; j += VecOp::SIZE )
{
VT p0 = vop.load(row0+j-cn), p1 = vop.load(row0+j), p2 = vop.load(row0+j+cn);
VT p3 = vop.load(row1+j-cn), p4 = vop.load(row1+j), p5 = vop.load(row1+j+cn);
VT p6 = vop.load(row2+j-cn), p7 = vop.load(row2+j), p8 = vop.load(row2+j+cn);

vop(p1, p2); vop(p4, p5); vop(p7, p8); vop(p0, p1);
vop(p3, p4); vop(p6, p7); vop(p1, p2); vop(p4, p5);
vop(p7, p8); vop(p0, p3); vop(p5, p8); vop(p4, p7);
vop(p3, p6); vop(p1, p4); vop(p2, p5); vop(p4, p7);
vop(p4, p2); vop(p6, p4); vop(p4, p2);
vop.store(dst+j, p4);
}

limit = size.width;
}
}
}
else if( m == 5 )
{
// 类似于 m == 3
}
}

形态学

有关形态学操作的代码在modules/imgproc/src/morph.dispatch.cpp中。

其中getStructuringElement()方法用于构造形态学操作的核。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
Mat getStructuringElement(int shape, Size ksize, Point anchor)
{
int i, j;
int r = 0, c = 0;
double inv_r2 = 0;

anchor = normalizeAnchor(anchor, ksize);
if( ksize == Size(1,1) )
shape = MORPH_RECT;
if( shape == MORPH_ELLIPSE )
{
r = ksize.height/2;
c = ksize.width/2;
inv_r2 = r ? 1./((double)r*r) : 0;
}

Mat elem(ksize, CV_8U);
for( i = 0; i < ksize.height; i++ )
{
uchar* ptr = elem.ptr(i);
int j1 = 0, j2 = 0;
/// 形状为矩形或者十字刚好处于横着的那一行时,直接填充一整行
if( shape == MORPH_RECT || (shape == MORPH_CROSS && i == anchor.y) )
j2 = ksize.width;
else if( shape == MORPH_CROSS )
j1 = anchor.x, j2 = j1 + 1;
else
{
int dy = i - r;
if( std::abs(dy) <= r )
{
/// 计算近似的椭圆的宽度
int dx = saturate_cast<int>(c*std::sqrt((r*r - dy*dy)*inv_r2));
j1 = std::max( c - dx, 0 );
j2 = std::min( c + dx + 1, ksize.width );
}
}
/// 填充 j1 ~ j2 区间内的数字
for( j = 0; j < j1; j++ )
ptr[j] = 0;
for( ; j < j2; j++ )
ptr[j] = 1;
for( ; j < ksize.width; j++ )
ptr[j] = 0;
}

return elem;
}

形态学中常见的腐蚀、膨胀操作:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
void erode( InputArray src, OutputArray dst, InputArray kernel,
Point anchor, int iterations,
int borderType, const Scalar& borderValue )
{
CV_INSTRUMENT_REGION();

CV_Assert(!src.empty());

morphOp( MORPH_ERODE, src, dst, kernel, anchor, iterations, borderType, borderValue );
}


void dilate( InputArray src, OutputArray dst, InputArray kernel,
Point anchor, int iterations,
int borderType, const Scalar& borderValue )
{
CV_INSTRUMENT_REGION();

CV_Assert(!src.empty());

morphOp( MORPH_DILATE, src, dst, kernel, anchor, iterations, borderType, borderValue );
}

morphOp() 方法中会提供一个空的 cv_hal_morph 方法供用户自行定义 morph 方法的实现。若没有自己定义的实现,则调用 opencv 内置提供的 ocvMorph 方法。和其他滤波器类似,在该方法中,调用了 createMorphologyFilter 得到一个 FilterEngine,最后调用 apply 方法进行计算。最后实际进行图形学滤波运算的是 MorphFilter 这样一个模板类:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

template<class Op, class VecOp> struct MorphFilter : BaseFilter
{
typedef typename Op::rtype T;
MorphFilter( const Mat& _kernel, Point _anchor ) { ... }

void operator()(const uchar** src, uchar* dst, int dststep, int count, int width, int cn) CV_OVERRIDE
{
...
width *= cn;
/// 遍历每一行
for( ; count > 0; count--, dst += dststep, src++ )
{
...
/// 遍历每一列
for( ; i < width; i++ )
{
T s0 = kp[0][i];
/// 滤波操作
for( k = 1; k < nz; k++ )
s0 = op(s0, kp[k][i]);
D[i] = s0;
}
}
}
};

对于 erode 和 dilate 两种操作,只需要分别传入 MinOp(返回值更小的那个) 和 MaxOp(返回值更大的那个) 即可。以 erode 为例,传入 MinOp 之后,对于核上每一个为 1 的点,覆盖到图像上的对应位置也必须为 1,否则由于 min 操作的特性,只要有一个是 0 最后的结果就会是 0,这个操作的结果就是,将核中心放在结果图像上任意一个为 1 的点,都能够被原图像包裹,即结果图像是源图像的腐蚀。膨胀则使用最大值,分析类似。

形态学中的其它操作基本都转化为 erode 和 dilate 操作,由 morphologyEx() 方法可见:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
switch( op )
{
case MORPH_ERODE: /// 腐蚀操作
erode( src, dst, kernel, anchor, iterations, borderType, borderValue );
break;
case MORPH_DILATE: /// 扩张操作
dilate( src, dst, kernel, anchor, iterations, borderType, borderValue );
break;
case MORPH_OPEN: /// 开操作
erode( src, dst, kernel, anchor, iterations, borderType, borderValue );
dilate( dst, dst, kernel, anchor, iterations, borderType, borderValue );
break;
case MORPH_CLOSE: /// 闭操作
dilate( src, dst, kernel, anchor, iterations, borderType, borderValue );
erode( dst, dst, kernel, anchor, iterations, borderType, borderValue );
break;
case MORPH_GRADIENT: /// 梯度计算操作
erode( src, temp, kernel, anchor, iterations, borderType, borderValue );
dilate( src, dst, kernel, anchor, iterations, borderType, borderValue );
dst -= temp;
break;
case MORPH_TOPHAT: /// 顶帽操作
if( src.data != dst.data )
temp = dst;
erode( src, temp, kernel, anchor, iterations, borderType, borderValue );
dilate( temp, temp, kernel, anchor, iterations, borderType, borderValue );
dst = src - temp;
break;
case MORPH_BLACKHAT: /// 黑帽操作
if( src.data != dst.data )
temp = dst;
dilate( src, temp, kernel, anchor, iterations, borderType, borderValue );
erode( temp, temp, kernel, anchor, iterations, borderType, borderValue );
dst = temp - src;
break;
...
}

图像分割

Canny 边缘检测

OpenCV 中实施 Canny 边缘检测的一般形式:

1
Canny( detected_edges, detected_edges, lowThreshold, lowThreshold*ratio, kernel_size );

源代码位于modules/imgproc/src/canny.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
// image	8 bit 输入图像
// edges 输出的图像边缘,单通道 8 bit 图像,和原图像有同样大小
// threshold1 低阈值
// threshold2 高阈值
// apertureSize Sobel 孔径大小
// L2gradient flag,是否采用更精确的 L2 梯度
void Canny( InputArray _src, OutputArray _dst,
double low_thresh, double high_thresh,
int aperture_size, bool L2gradient )
{

// 省略一些检查和初始化工作

const Size size = _src.size();

_dst.create(size, CV_8U);

Mat src0 = _src.getMat(), dst = _dst.getMat();
Mat src(src0.size(), src0.type(), src0.data, src0.step);

/// 尝试是否可以通过 OpenCL、HAL、OpenVX、IPP 完成 Canny 算法的计算
CV_OCL_RUN(...)
CALL_HAL(...);
CV_OVX_RUN(...)
CV_IPP_RUN_FAST(...)

// 如果使用 L2 梯度,需要修正 thresh
if (L2gradient)
{
low_thresh = std::min(32767.0, low_thresh);
high_thresh = std::min(32767.0, high_thresh);

if (low_thresh > 0) low_thresh *= low_thresh;
if (high_thresh > 0) high_thresh *= high_thresh;
}
int low = cvFloor(low_thresh);
int high = cvFloor(high_thresh);

// If Scharr filter: aperture size is 3, ksize2 is 1
int ksize2 = aperture_size < 0 ? 1 : aperture_size / 2;
// 计算可用的线程数
int numOfThreads = std::max(1, std::min(getNumThreads(), getNumberOfCPUs()));
// Make a fallback for pictures with too few rows.
int grainSize = src.rows / numOfThreads;
int minGrainSize = 2 * (ksize2 + 1);
if (grainSize < minGrainSize)
numOfThreads = std::max(1, src.rows / minGrainSize);

Mat map;
std::deque<uchar*> stack;

// 并行计算 Canny 算法
parallel_for_(Range(0, src.rows), parallelCanny(src, map, stack, low, high, aperture_size, L2gradient), numOfThreads);

CV_TRACE_REGION("global_hysteresis");
// 全局进行 edge track
ptrdiff_t mapstep = map.cols;

while (!stack.empty())
{
uchar* m = stack.back();
stack.pop_back();

if (!m[-mapstep-1]) CANNY_PUSH((m-mapstep-1), stack);
if (!m[-mapstep]) CANNY_PUSH((m-mapstep), stack);
if (!m[-mapstep+1]) CANNY_PUSH((m-mapstep+1), stack);
if (!m[-1]) CANNY_PUSH((m-1), stack);
if (!m[1]) CANNY_PUSH((m+1), stack);
if (!m[mapstep-1]) CANNY_PUSH((m+mapstep-1), stack);
if (!m[mapstep]) CANNY_PUSH((m+mapstep), stack);
if (!m[mapstep+1]) CANNY_PUSH((m+mapstep+1), stack);
}

CV_TRACE_REGION_NEXT("finalPass");
// 转换 map 中的点,显示出边缘
parallel_for_(Range(0, src.rows), finalPass(map, dst), src.total()/(double)(1<<16));
}

主要函数为 parallelCanny(src, map, stack, low, high, aperture_size, L2gradient),源码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
if(needGradient)
{
if (aperture_size == 7)
{
scale = 1 / 16.0;
}
// 如果需要梯度,用 Sobel 算子计算
Sobel(src.rowRange(rowStart, rowEnd), dx, CV_16S, 1, 0, aperture_size, scale, 0, BORDER_REPLICATE);
Sobel(src.rowRange(rowStart, rowEnd), dy, CV_16S, 0, 1, aperture_size, scale, 0, BORDER_REPLICATE);
}
else
{
dx = src.rowRange(rowStart, rowEnd);
dy = src2.rowRange(rowStart, rowEnd);
}

计算强度和角度梯度,应用非最大抑制:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92

// 用值填充 map
// 0 - 像素可能属于边界
// 1 - 像素不属于边界
// 2 - 像素属于边界
for (int i = rowStart; i <= boundaries.end; ++i)
{
// 梯度计算

if(i < rowEnd)
{
// 计算下一行
_dx = dx.ptr<short>(i - rowStart);
_dy = dy.ptr<short>(i - rowStart);

// 使用 L2 梯度
if (L2gradient)
{
int j = 0, width = src.cols * cn;
// 省略 SIMD 相关代码
for ( ; j < width; ++j)
_mag_n[j] = int(_dx[j])*_dx[j] + int(_dy[j])*_dy[j];
}
else // 使用 L1 梯度
{
int j = 0, width = src.cols * cn;
for ( ; j < width; ++j)
_mag_n[j] = std::abs(int(_dx[j])) + std::abs(int(_dy[j]));
}
...
}
else
{
...
}

...

// 非极大值抑制
const int TG22 = 13573;
int j = 0;

for (; j < src.cols; j++)
{
int m = _mag_a[j];

if (m > low)
{
short xs = _dx[j];
short ys = _dy[j];
int x = (int)std::abs(xs);
int y = (int)std::abs(ys) << 15;

int tg22x = x * TG22;

// 水平方向梯度
if (y < tg22x)
{
if (m > _mag_a[j - 1] && m >= _mag_a[j + 1])
{
// 检查是否大于 high 阈值,是入栈且 pmap 设为 0
CANNY_CHECK(m, high, (_pmap+j), stack);
continue;
}
}
else
{
// 垂直方向梯度
int tg67x = tg22x + (x << 16);
if (y > tg67x)
{
if (m > _mag_p[j] && m >= _mag_n[j])
{
CANNY_CHECK(m, high, (_pmap+j), stack);
continue;
}
}
else
{
// 斜方向梯度
int s = (xs ^ ys) < 0 ? -1 : 1;
if(m > _mag_p[j - s] && m > _mag_n[j + s])
{
CANNY_CHECK(m, high, (_pmap+j), stack);
continue;
}
}
}
}
_pmap[j] = 1;
}
}

进行 egde track,即使用滞后阈值进行改进:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
// now track the edges (hysteresis thresholding)
CV_TRACE_REGION_NEXT("hysteresis");
while (!stack.empty())
{
uchar *m = stack.back();
stack.pop_back();
// 如果不是位于边界
if((unsigned)(m - pmapLower) < pmapDiff)
{
if (!m[-mapstep-1]) CANNY_PUSH((m-mapstep-1), stack);
if (!m[-mapstep]) CANNY_PUSH((m-mapstep), stack);
if (!m[-mapstep+1]) CANNY_PUSH((m-mapstep+1), stack);
if (!m[-1]) CANNY_PUSH((m-1), stack);
if (!m[1]) CANNY_PUSH((m+1), stack);
if (!m[mapstep-1]) CANNY_PUSH((m+mapstep-1), stack);
if (!m[mapstep]) CANNY_PUSH((m+mapstep), stack);
if (!m[mapstep+1]) CANNY_PUSH((m+mapstep+1), stack);
}
else
{
// 处理边界
borderPeaksLocal.push_back(m);
ptrdiff_t mapstep2 = m < pmapLower ? mapstep : -mapstep;

if (!m[-1]) CANNY_PUSH((m-1), stack);
if (!m[1]) CANNY_PUSH((m+1), stack);
if (!m[mapstep2-1]) CANNY_PUSH((m+mapstep2-1), stack);
if (!m[mapstep2]) CANNY_PUSH((m+mapstep2), stack);
if (!m[mapstep2+1]) CANNY_PUSH((m+mapstep2+1), stack);
}
}

最后一步,将边缘点(标记为2)映射为灰度值 255(白),其它映射到 0:

1
2
3
4
5
6
7
8
9
10
11
12
// the final pass, form the final image
for (int i = boundaries.start; i < boundaries.end; i++)
{
int j = 0;
uchar *pdst = dst.ptr<uchar>(i);
const uchar *pmap = map.ptr<uchar>(i + 1);
pmap += 1;
for (; j < dst.cols; j++)
{
pdst[j] = (uchar)-(pmap[j] >> 1);
}
}

>>1后只有 2 变为了 1,uchar - 1 变成了 255,其它的 0 和 1,则映射到了 0。这样就得到了边缘图像。

OTSU

彩色图像处理

OpenCV-Python

OpenCV 中,所有的算法都是通过 C++ 实现。但也可以通过其它的语言调用,如 Java,Python。绑定生成器使这成为可能。这些生成器在C ++和Python之间建立了桥梁,使用户能够从Python调用C ++函数。

Python 官方提供了将 C++ 扩展到 Python 的示例。但 OpenCV 在modules/python/src2中的一些脚本,根据 C++ 头文件自动生成包装器函数。

modules/python/CMakeFiles.txt是一个CMake脚本,用于检查要扩展到Python的模块。它将自动检查所有要扩展的模块并获取其头文件。这些头文件包含该特定模块的所有类,函数,常量等的列表。

其次,将这些头文件传递到Python脚本modules/python/src2/gen2.py。这是Python Binding生成器脚本。它调用另一个Python脚本module/python/src2/hdr_parser.py。这是标头解析器脚本。此标头解析器将完整的标头文件拆分为较小的Python列表。因此,这些列表包含有关特定函数,类等的所有详细信息。例如,将对一个函数进行解析以获取一个包含函数名称,返回类型,输入参数,参数类型等的列表。最终列表包含所有函数,枚举的详细信息,头文件中的structs,classs等。

因此头解析器将返回已解析函数的最终大列表。我们的生成器脚本(gen2.py)将为标头解析器解析的所有函数/类/枚举/结构创建包装函数(你可以在编译期间在build/modules/python/文件夹中以pyopencv_genic_*.h文件找到这些标头文件)。但是可能会有一些基本的OpenCV数据类型,例如Mat,Vec4i,Size。它们需要手动扩展。例如,Mat类型应扩展为Numpy数组,Size应扩展为两个整数的元组,等等。类似地,可能会有一些复杂的结构/类/函数等需要手动扩展。所有这些手动包装函数都放在modules/python/src2/cv2.cpp中。

所以现在剩下的就是这些包装文件的编译了,这给了我们cv2模块。因此,当你调用函数时,例如在Python中说res = equalizeHist(img1,img2),你将传递两个numpy数组,并期望另一个numpy数组作为输出。因此,将这些numpy数组转换为cv::Mat,然后在C++中调用equalizeHist()函数。最终结果将res转换回Numpy数组。简而言之,几乎所有操作都是在 C++ 中完成的,这给了我们几乎与C++相同的速度。

参考

OpenCV: OpenCV Tutorials

OpenCV: OpenCV modules

opencv 源码初探 | Little csd’s blog

OpenCV: How OpenCV-Python Bindings Works?

OpenCVTutorials/11_1_OpenCV-Python Bindings.md at master · fendouai/OpenCVTutorials

数字图像处理(第五版) 电子工业出版社