Skip to content

OpenCV 源码分析

Posted on:January 12, 2021 at 18:42:18 GMT+8

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 类的定义:

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 的构造函数之一:

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() 函数来创建对象:

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);
}

读取、写入、显示图像

从文件中读取图像:

Mat img = imread(filename);

灰度图像:

Mat img = imread(filename, IMREAD_GRAYSCALE);

将图像写入文件:

imwrite(filename, img);

展示图像:

imshow( "Title", img );

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

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_() 方法实现:

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

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() 源码:

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,函数定义如下:

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

内部调用了:

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++ 中:

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

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

void blur(InputArray src, OutputArray dst,
          Size ksize, Point anchor, int borderType)
{
    CV_INSTRUMENT_REGION();

    boxFilter( src, dst, -1, ksize, anchor, true, borderType );
}
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

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

medianBlur ( src, dst, 3 );

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

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 中内置的中值滤波的实现方法:

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()方法:

// 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()方法用于构造形态学操作的核。

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;
}

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

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 这样一个模板类:


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() 方法可见:

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 边缘检测的一般形式:

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

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

// 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),源码:

        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);
        }

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


        //   用值填充 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,即使用滞后阈值进行改进:

        // 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:

        // 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/文件夹中以pyopencvgenic*.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

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