cvMorphology形态学原理解析及源码分析

⑴ 图像形态学处理的概念...1

⑵ 二值图像的逻辑运算...3

⑶ 膨胀和腐蚀...4

(4) 高级形态学变换...8

(5) 细化...10

 

⑴ 图像形态学处理的概念

数字图像处理中的形态学处理是指将数字形态学作为工具从图像中提取对于表达和描绘区域形状有用处的图像分量，比如边界、骨架以及凸壳，还包括用于预处理或后处理的形态学过滤、细化和修剪等。图像形态学处理中我们感兴趣的主要是二值图像。

在二值图像中，所有黑色像素的集合是图像完整的形态学描述，二值图像的各个分量是Z2的元素。假定二值图像A和形态学处理的结构元素B是定义在笛卡儿网格上的集合，网格中值为1的点是集合的元素，当结构元素的原点移到点(x,y)时，记为Sxy，为简单起见，结构元素为3x3，且全都为1，在这种限制下，决定输出结果的是逻辑运算。

IplConvKernel结构元素

typedef struct _IplConvKernel

{

    int  nCols; //结构元素的行宽

    int  nRows; //列高

    int anchorX; //结构原点位置水平坐标

    int anchorY; //结构原点位置垂直坐标

int *values; //当nShiftR为自定义时，value是指向结构元素数据的指针

//如果结构元素的大小定义为8*6，那么values为48长的int数组，值为0或1。

    int nShiftR;// 用于表示结构元素的形状类型

}IplConvKernel;

cvCreateStructuringElementEx创建结构元素 

IplConvKernel* cvCreateStructuringElementEx( int cols, int rows, intanchor_x, int anchor_y,

                                                    int shape, int*values=NULL );

cols 结构元素的列数目 rows 结构元素的行数目

anchor_x 锚点的相对水平偏移量 anchor_y 锚点的相对垂直偏移量 

shape 结构元素的形状，可以是下列值： 

CV_SHAPE_RECT, 长方形元素; 

CV_SHAPE_CROSS, 十字交叉型，交错元素 a cross-shaped element; 

CV_SHAPE_ELLIPSE, 椭圆元素; 

CV_SHAPE_CUSTOM, 用户自定义元素。这种情况下参数 values 在封闭矩形内定义核的形状,即象素的那个邻域必须考虑。

values 指向结构元素的指针，它是一个平面数组，表示对元素矩阵逐行扫描。(非零点表示该点属于结构元)。如果指针为空，则表示平面数组中的所有元素都是非零的，即结构元是一个长方形(该参数仅仅当shape参数是 CV_SHAPE_CUSTOM 时才予以考虑)。 

形态核与卷积核不同，不需要任何数值填充核。当核在图像上移动时，核的元素只需要简单表明应该在哪个范围内计算最大值和最小值，参考点制定核与源图像的位置关系，同时也锁定了计算结果在目标图像中的位置。行和列确定了所构造的矩形的大小(结构元素在矩形内)，anchor_x和anchor_y是核的封闭矩形内的参考点坐标。

cvReleaseStructuringElement删除结构元素 

void cvReleaseStructuringElement( IplConvKernel** element );

element 被删除的结构元素的指针,函数 cvReleaseStructuringElement 释放结构 IplConvKernel 。如果 *element 为 NULL, 则函数不作用。

CV_IMPL IplConvKernel*

cvCreateStructuringElementEx( int cols, introws,

                              intanchorX,int anchorY,

                              intshape,int *values )

{

    cv::Sizeksize =cv::Size(cols,rows);

    cv::Pointanchor =cv::Point(anchorX,anchorY);

       // 检测输入数据，当用户自定义的时候value不能为空，value默认为NULL

    CV_Assert(cols > 0 &&rows > 0 &&anchor.inside(cv::Rect(0,0,cols,rows)) &&

               (shape!= CV_SHAPE_CUSTOM || values != 0));

 

    int i, size = rows *cols;

    int element_size = sizeof(IplConvKernel) +size*sizeof(int);

       // 为什么创建的内存要比实际的大呢？大了size*sizeof(int)+32

    IplConvKernel*element = (IplConvKernel*)cvAlloc(element_size+ 32);

 

    element->nCols =cols;

    element->nRows =rows;

    element->anchorX =anchorX;

    element->anchorY =anchorY;

//   enum     {            CV_SHAPE_RECT      =0,             CV_SHAPE_CROSS     =1,           

//          CV_SHAPE_ELLIPSE   =2,             CV_SHAPE_CUSTOM    =100       };

    element->nShiftR =shape<CV_SHAPE_ELLIPSE ?shape :CV_SHAPE_CUSTOM;

       // element指向结构的首地址

    element->values = (int*)(element + 1);

       // 如果为用户自定义的类型，从values中取值

    if( shape == CV_SHAPE_CUSTOM)

    {

        for( i = 0; i < size; i++ )

            element->values[i] =values[i];

    }

    else

    {

              // 根据不同的结构类型获得不同的数值

        cv::Matelem =cv::getStructuringElement(shape,ksize,anchor);

        for( i = 0; i < size; i++ )

            element->values[i] =elem.data[i];

    }

 

    return element;

}

cv::Matcv::getStructuringElement(intshape,Sizeksize, Pointanchor)

{

    int i, j;

    int r = 0, c = 0;

    double inv_r2 = 0;

      

    CV_Assert(shape ==MORPH_RECT||shape ==MORPH_CROSS||shape ==MORPH_ELLIPSE);

       //ifanchor.x=-1,anchor.x=ksize.width/2; if anchor.y=-1,anchor.y=ksize.height/2

       // 并判断是否在rect(0, 0, ksize.width, ksize.height)内

    anchor =normalizeAnchor(anchor,ksize);

       // 当只有一个结构元素的时候cols=1 rows=1，长方形结构元素

    if( ksize == Size(1,1))

        shape= MORPH_RECT;

       // 如果为椭圆形的结构元素

    if( shape == MORPH_ELLIPSE)

    {

              //r  c分别为椭圆的半径

        r = ksize.height/2;

        c = ksize.width/2;

              // 如果r!=0，inv_r2=1/(r*r)

        inv_r2= r ? 1./((double)r*r) : 0;

    }

      

    Mat elem(ksize, CV_8U);

 

    for( i = 0; i < ksize.height; i++ )

    {

        uchar*ptr =elem.data +i*elem.step;

        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

        {

            intdy =i -r;

            if(std::abs(dy) <=r )

            {

                intdx =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的赋，大于j1小于j2的赋，其余的赋0

        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;

}

⑵ 二值图像的逻辑运算

逻辑运算尽管本质上很简单，但对于实现以形态学为基础的图像处理算法是一种有力的补充手段。在图像处理中用到的主要逻辑运算是：与、或和非（求补），它们可以互相组合形成其他逻辑运算。

⑶ 膨胀和腐蚀

膨胀和腐蚀这两种操作是形态学处理的基础，许多形态学算法都是以这两种运算为基础的。

① 膨胀

结构元素B可以看作一个卷积模板，区别在于膨胀是以集合运算为基础的，卷积是以算术运算为基础的，但两者的处理过程是相似的。

⑴ 用结构元素B，扫描图像A的每一个像素

⑵ 用结构元素与其覆盖的二值图像做“与”操作

⑶ 如果都为0，结果图像的该像素为0。否则为1

cvDilate使用任意结构元素膨胀图像 

voidcvDilate( const CvArr* src, CvArr* dst, IplConvKernel* element=NULL, intiterations=1 );

src输入图像. dst 输出图像. element 用于膨胀的结构元素。

若为 NULL,则使用(1+iterations*2)*(1+iterations*2)的长方形的结构元素 .iterations膨胀的次数 

函数 cvDilate 对输入图像使用指定的结构元进行膨胀，该结构决定每个具有最小值象素点的邻域形状： dst=dilate(src,element): dst(x,y)=max((x',y') in element))src(x+x',y+y')

函数支持（in-place）模式。膨胀可以重复进行 (iterations) 次. 对彩色图像，每个彩色通道单独处理。

在试图找到连通分支（即具有相似的颜色或强度的像素点的大块互相分离的区域）时通常使用膨胀操作。

CV_IMPL void

cvDilate( const CvArr* srcarr, CvArr* dstarr, IplConvKernel* element,int iterations)

{

    cv::Matsrc =cv::cvarrToMat(srcarr),dst =cv::cvarrToMat(dstarr),kernel;

       // 输入输出必须是同等尺寸、同类型的

    CV_Assert(src.size()==dst.size()&&src.type()==dst.type());

    cv::Pointanchor;

       // 若没有结构元素输入，kernel=NULL，anchor=(1,1)

       // 否则将结构元素中的值读入kernel和anchor

    convertConvKernel(element,kernel,anchor );

       // 边界差值方法采用边界复制

    cv::dilate(src,dst, kernel,anchor, iterations,cv::BORDER_REPLICATE);

}

static voidconvertConvKernel(constIplConvKernel*src,cv::Mat&dst,cv::Point&anchor)

{

       // 若没有输入结构元素

    if(!src)

    {

        anchor= cv::Point(1,1);

        dst.release();

        return;

    }

       // 获取结构原点的坐标

    anchor =cv::Point(src->anchorX,src->anchorY);

       // 读取结构元素的值

    dst.create(src->nRows,src->nCols,CV_8U);

    int i, size = src->nRows*src->nCols;

    for( i = 0; i < size; i++ )

        dst.data[i] = (uchar)src->values[i];

}

void cv::dilate( InputArraysrc, OutputArraydst, InputArraykernel,

                 Pointanchor,intiterations,

                 intborderType,constScalar&borderValue)

{

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

}

static voidmorphOp(int op, InputArray _src, OutputArray_dst,

                     InputArray_kernel,

                     Pointanchor,intiterations,

                     intborderType,constScalar&borderValue)

{

    Mat src = _src.getMat(),kernel=_kernel.getMat();

       // 如果输入的时候不输入kernel，则kernel.data=NULL，那么ksize=(3,3)

    Size ksize = kernel.data ?kernel.size() :Size(3,3);

       // ifanchor.x=-1,anchor.x=ksize.width/2; if anchor.y=-1,anchor.y=ksize.height/2

       // 并判断是否在rect(0, 0, ksize.width, ksize.height)内

    anchor =normalizeAnchor(anchor,ksize);

       // 这一句是多余的，因为在上面normalizeAnchor已经判断了

    CV_Assert(anchor.inside(Rect(0, 0,ksize.width,ksize.height)) );

 

    _dst.create(src.size(),src.type() );

    Mat dst = _dst.getMat();

       // 如果迭代步数为或者结构元素的尺寸为，不进行处理，直接输出

    if( iterations == 0 || kernel.rows*kernel.cols == 1 )

    {

        src.copyTo(dst);

        return;

    }

       // 如果没有输入结构元素，那么创建(1+iterations*2)*(1+iterations*2)的长方形结构元素

       // 结构元素的中心点为(iterations, iterations)，并将迭代步数设置为

    if( !kernel.data )

    {

        kernel= getStructuringElement(MORPH_RECT, Size(1+iterations*2,1+iterations*2));

        anchor= Point(iterations,iterations);

        iterations= 1;

    }

       // 如果结构步数大于的话并且kernel为长方形的结构元素，重新创建结构元素

    else if( iterations> 1 && countNonZero(kernel) == kernel.rows*kernel.cols )

    {

        anchor= Point(anchor.x*iterations,anchor.y*iterations);

        kernel= getStructuringElement(MORPH_RECT,

                                       Size(ksize.width + (iterations-1)*(ksize.width-1),

                                            ksize.height +(iterations-1)*(ksize.height-1)),

                                       anchor);

        iterations= 1;

    }

 

int nStripes = 1;

// Tegra是NVIDIA公司于2008年推出的基于ARM构架通用处理器品牌（即CPU，NVIDIA称为“Computer on a chip”片上计算机），能够为便携设备提供高性能、低功耗体验。

#if definedHAVE_TEGRA_OPTIMIZATION

    if (src.data != dst.data &&iterations == 1 &&  //NOTE:threads are not used for inplace processing

        (borderType & BORDER_ISOLATED) == 0&& //TODO: check border types

        src.rows >= 64 ) //NOTE: justheuristics

        nStripes = 4;

#endif

 

    parallel_for_(Range(0,nStripes),

                  MorphologyRunner(src,dst,nStripes,iterations,op,kernel,anchor,borderType,borderType,borderValue));

 

    //Ptr<FilterEngine>f = createMorphologyFilter(op, src.type(),

    //                                            kernel, anchor, borderType, borderType, borderValue );

 

    //f->apply(src, dst );

    //for( int i = 1;i < iterations; i++ )

    //    f->apply( dst, dst );

}

// 是否采用并行处理

void cv::parallel_for_(constcv::Range&range,constcv::ParallelLoopBody&body,doublenstripes=-1)

{

       // 大部分代码省略，如果定义了并行框架，可以采用并行处理，一般不定义

    (void)nstripes;

    body(range);

}

class MorphologyRunner: public ParallelLoopBody

{

public:

    MorphologyRunner(Mat_src,Mat _dst, int _nStripes,int_iterations,

                     int_op,Mat_kernel,Point_anchor,

                     int_rowBorderType,int_columnBorderType,constScalar&_borderValue):

                                  borderValue(_borderValue)

    {

        src= _src;

        dst= _dst;

 

        nStripes= _nStripes;

        iterations= _iterations;

 

        op =_op;

        kernel= _kernel;

        anchor= _anchor;

        rowBorderType= _rowBorderType;

        columnBorderType= _columnBorderType;

    }

       // ()操作符，最主要的运算符号

    void operator () ( const Range& range) const

    {

        int row0 = min(cvRound(range.start *src.rows / nStripes),src.rows);

        int row1 = min(cvRound(range.end *src.rows /nStripes),src.rows);

 

        /*if(0)

            printf("Size = (%d, %d),range[%d,%d), row0 = %d, row1 = %d\n",

                   src.rows, src.cols,range.start, range.end, row0, row1);*/

 

        Mat srcStripe = src.rowRange(row0,row1);

        Mat dstStripe = dst.rowRange(row0,row1);

               // 创建形态学滤波器

        Ptr<FilterEngine>f=createMorphologyFilter(op,src.type(),kernel,anchor,

                                                    rowBorderType,columnBorderType,borderValue);

              // 主要的处理步骤在这里面，还未解读

        f->apply(srcStripe,dstStripe );

        for( int i = 1; i < iterations;i++ )

            f->apply(dstStripe,dstStripe );

    }

 

private:

    Mat src;

    Mat dst;

    int nStripes;

    int iterations;

 

    int op;

    Mat kernel;

    Point anchor;

    int rowBorderType;

    int columnBorderType;

    Scalar borderValue;

};

② 腐蚀

对Z中的集合A和B，B对A进行腐蚀的整个过程如下：

⑴ 用结构元素B，扫描图像A的每一个像素

⑵ 用结构元素与其覆盖的二值图像做“与”操作

⑶ 如果都为1，结果图像的该像素为1。否则为0

腐蚀处理的结果是使原来的二值图像减小一圈。

cvErode使用任意结构元素腐蚀图像 

void cvErode( const CvArr* src, CvArr* dst,IplConvKernel* element=NULL, int iterations=1 );

src 输入图像. dst 输出图像.element 用于腐蚀的结构元素。

若为 NULL, 则使用 (1+iterations*2)*(1+iterations*2)的长方形的结构元素iterations 腐蚀的次数 

函数 cvErode 对输入图像使用指定的结构元素进行腐蚀，该结构元素决定每个具有最小值象素点的邻域形状： dst=erode(src,element):  dst(x,y)=min((x',y') inelement))src(x+x',y+y')

函数可能是本地操作支持in-place，不需另外开辟存储空间的意思。腐蚀可以重复进行 (iterations) 次. 对彩色图像，每个彩色通道单独处理。 

腐蚀操作通常消除图像中的斑点噪声，确保图像中较大的区域仍然存在。

cvErode的源代码与cvDialte的源代码相似，在此不再对其进行解读

CV_IMPL void

cvErode( const CvArr* srcarr, CvArr* dstarr, IplConvKernel* element,int iterations)

{

    cv::Matsrc =cv::cvarrToMat(srcarr),dst =cv::cvarrToMat(dstarr),kernel;

    CV_Assert(src.size()==dst.size()&&src.type()==dst.type());

    cv::Pointanchor;

    convertConvKernel(element,kernel,anchor );

    cv::erode(src,dst, kernel,anchor, iterations,cv::BORDER_REPLICATE);

}

void cv::erode( InputArraysrc, OutputArraydst, InputArraykernel,

                Pointanchor,intiterations,

                intborderType,constScalar&borderValue)

{

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

}

(4) 高级形态学变换

开操作是先腐蚀、后膨胀处理。

闭操作是先膨胀、后腐蚀处理。

cvMorphologyEx高级形态学变换 

void cvMorphologyEx( const CvArr* src, CvArr* dst, CvArr* temp,

                    IplConvKernel* element, int operation, int iterations=1 );

src 输入图像. dst 输出图像. temp 临时图像，某些情况下需要，与源图像同样大小。临时图像 temp 在形态梯度以及对“顶帽”和“黑帽”操作时的 in-place 模式下需要。element 结构元素，如果没有输入，则使用3*3的长方形结构元素。 iterations 迭代次数. 

operation 形态操作的类型: CV_MOP_OPEN - 开运算 CV_MOP_CLOSE - 闭运算 CV_MOP_GRADIENT - 形态梯度 CV_MOP_TOPHAT - "顶帽" CV_MOP_BLACKHAT - "黑帽" 

函数 cvMorphologyEx 在膨胀和腐蚀基本操作的基础上，完成一些高级的形态变换： 

开运算 dst=open(src,element)=dilate(erode(src,element),element)

开运算通常可以用来统计二值图像中的区域数。

闭运算 dst=close(src,element)=erode(dilate(src,element),element)

在多数连通域分析方法中用闭运算去除噪声区域

对于连通域分析，通常先采用腐蚀或者闭运算来消除纯粹噪声引起的部分，然后用开运算来连接邻近的区域。闭运算消除低于其邻近点的孤立点，开运算消除高于其邻近点的孤立点。对于iterations=2，就开运算而言其实是腐蚀->腐蚀->膨胀->膨胀这样的过程。

形态梯度dst=morph_grad(src,element)=dilate(src,element)-erode(src,element)

对图像进行这一操作，可以将团块blob的边缘以高亮区域突出出来，保留完整的外围边缘。

"顶帽" dst=tophat(src,element)=src-open(src,element) 

"黑帽" dst=blackhat(src,element)=close(src,element)-src 

当试图孤立的部分相对于其邻近的部分有亮度变化时可以使用，分离比邻近的点亮或暗的一些斑块。开运算带来的结果是放大裂缝或局部低亮度区域，因此顶帽操作可以突出与核大小相关的比源图像周围的区域更明亮的区域。黑帽操作突出比源图像周围的区域黑暗的区域。

CV_IMPL void

cvMorphologyEx( const void* srcarr,void*dstarr, void*,

                IplConvKernel*element,intop,int iterations )

{

    cv::Matsrc =cv::cvarrToMat(srcarr),dst =cv::cvarrToMat(dstarr),kernel;

    CV_Assert(src.size()==dst.size()&&src.type()==dst.type());

    cv::Pointanchor;

    IplConvKernel*temp_element =NULL;

       // 如果没有给定结构元素，则定义*3的长方形元素，元素原点为(1,1)

    if (!element)

    {

        temp_element= cvCreateStructuringElementEx(3, 3, 1, 1, CV_SHAPE_RECT);

    } else {

        temp_element= element;

    }

       // 读取结构元素中的值

    convertConvKernel(temp_element,kernel,anchor );

       // 释放定义的结构元素

    if (!element)

    {

        cvReleaseStructuringElement(&temp_element);

    }

       // 执行形态学操作

    cv::morphologyEx(src,dst,op, kernel,anchor,iterations,cv::BORDER_REPLICATE );

}

void cv::morphologyEx( InputArray_src, OutputArray_dst, int op,

                       InputArraykernel,Pointanchor,intiterations,

                       intborderType,constScalar&borderValue)

{

    Mat src = _src.getMat(),temp;

    _dst.create(src.size(),src.type());

    Mat dst = _dst.getMat();

 

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

        dilate(src,dst,kernel,anchor,iterations,borderType,borderValue );

        erode(dst,dst,kernel,anchor,iterations,borderType,borderValue );

        break;

    case CV_MOP_GRADIENT:

        erode(src,temp,kernel,anchor,iterations,borderType,borderValue );

        dilate(src,dst,kernel,anchor,iterations,borderType,borderValue );

        dst-= temp;

        break;

    case CV_MOP_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 CV_MOP_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;

    default:

        CV_Error(CV_StsBadArg,"unknownmorphological operation" );

    }

}

(5) 细化

图像细化一般作为一种图像预处理技术出现，目的是提取源图像的骨架，即是将原图像中线条宽度大于1个像素的线条细化成只有一个像素宽，形成“骨架”，形成骨架后能比较容易的分析图像，如提取图像的特征。

细化基本思想是“层层剥夺”，即从线条边缘开始一层一层向里剥夺，直到线条剩下一个像素的为止。图像细化大大地压缩了原始图像地数据量，并保持其形状的基本拓扑结构不变，从而为文字识别中的特征抽取等应用奠定了基础。细化算法应满足以下条件：

① 将条形区域变成一条薄线；

② 薄线应位与原条形区域的中心；

③ 薄线应保持原图像的拓扑特性。

细化分成串行细化和并行细化，串行细化即是一边检测满足细化条件的点，一边删除细化点；并行细化即是检测细化点的时候不进行点的删除只进行标记，而在检测完整幅图像后一次性去除要细化的点。

常用的图像细化算法有hilditch算法，pavlidis算法和rosenfeld算法等。注：进行细化算法前要先对图像进行二值化,即图像中只包含“黑”和“白”两种颜色。

cvThin

void cvThin( IplImage* src,IplImage* dst, int iterations=1)

功能：将IPL_DEPTH_8U型二值图像进行细化

参数：src原始IPL_DEPTH_8U型二值图像。dst目标存储空间，必须事先分配好，且和原图像大小类型一致。iterations，迭代次数

在opencv之前的版本中有，后来去除了