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

数字图像处理（第五版） 电子工业出版社