【图像处理】灰度图、亮度峰值极值查找

引言

图像分割中，采用腐蚀、膨胀等形态学操作一般可以对图像进行一定效果的分割。腐蚀采用的kernel往往决定了分割的效果，核太小则难以分割开来，核太大则会将某些部分腐蚀掉同时对原本图像形状产生较大的影响。所以，换了一个思路想，从外边分割不如从里面扩张。所以分割问题的重点放在了图像质心的查找。

图像处理中，针对二值化分割图可以使用距离变幻得到一张灰度图，而亮度中心则是质心，查找质心既可以实现简单的分割。Matlab中查找局部峰值的函数有imregionalmax，而在python中模块skimage中有peak_local_max函数同样实现了对灰度图的峰值查找。但是，在C/C++中，笔者还没有找到相关的图像处理库有这一方面的功能，所以只能对着论文写了。

1、距离变换

距离变换在图像处理应用十分广泛，特别在图像分割方面。简单来说，距离变换是针对二值化图像，其实就是计算非零点到最近零点的距离。所以，图像经过距离变换后，分割图像中心则呈现出峰值。距离的定义一般有以下3种类型：

OpenCV中，则有7中类型的距离定义，除此之外用于还可以使用自己定义的距离。（所以，有无限制种吧）

/** Distance types for Distance Transform and M-estimators */
enum
{
    CV_DIST_USER    =-1,  /**< User defined distance */
    CV_DIST_L1      =1,   /**< distance = |x1-x2| + |y1-y2| */
    CV_DIST_L2      =2,   /**< the simple euclidean distance */
    CV_DIST_C       =3,   /**< distance = max(|x1-x2|,|y1-y2|) */
    CV_DIST_L12     =4,   /**< L1-L2 metric: distance = 2(sqrt(1+x*x/2) - 1)) */
    CV_DIST_FAIR    =5,   /**< distance = c^2(|x|/c-log(1+|x|/c)), c = 1.3998 */
    CV_DIST_WELSCH  =6,   /**< distance = c^2/2(1-exp(-(x/c)^2)), c = 2.9846 */
    CV_DIST_HUBER   =7    /**< distance = |x|<c ? x^2/2 : c(|x|-c/2), c=1.345 */
};

以下是针对两个圆形粘连情况下采用CV_DIST_L2（常用欧式距离）作为距离进行距离变换：

2、形态学灰度重建（Morphological Grayscale Reconstruction）

一般对于局部峰值极值点的查找很容易联想到使用窗口判断的方法。但是，窗口的大小决定着算法的好坏。更致命的是，窗口这种方法很难应对出现连续峰值的问题（峰值平坦？不太会描述）。所以这种方法很难解决峰值查找的问题。于是，笔者在各种论坛中查找各种峰值查找的方法，最终还是回到Matlab参考文献。

Matlab中的imregionalmax函数是使用形态学灰度重建的方法实现了局部极大值的查找，具体参考论文：Morphological Grayscale Reconstruction in Image Analysis: Applications and Effcient Algorithms。这里只进行简单描述。

 从上图可知，对于图像I，用其偏置I-1进行grayscale reconstruction灰度重建得到重建后的图像。然后在使用原图与重建后的图像相减即可得到局部极大值部分。上图也可以解释Matlab中imregionalmax等于原图与其imreconstruct相减。接下来，则是如何使用算法实现这个过程呢？论文末尾处已经给出了许多算法实现的伪代码，其中包括并行算法、串行算法等等。本文使用队列对算法进行实现，该算法在论文中已经给出伪代码：

特别说明，笔者取J=I-1，R(J)是遍历每个点取八邻近最大值来初始化队列。NG(p)也是取得p点的八邻近。部分实现代码如下：

Mat findpeeks(Mat I,Mat J)
{   
    queue<Point> que;    
    /*最大值滤波*/
    J = regionalmaximaImg(J);
    /*对满足条件的点入队列，初始化队列*/
    for(int i=0;i<I.rows;i++)
        for(int j=0;j<I.cols;j++)
            if(J.at<uchar>(i,j) && isexistNBzero(J,Point(j,i)))
                que.push(Point(j,i));
    /*循环队列，直至队列为空*/
    while(que.empty())
    {
        Point p = que.front();
        que.pop();
        
        for(int i=0;i<NEGHTBOR;i++)
        {
            /*get the neighbor point*/
            Point q;
            q.x = p.x+x_offset[i];
            q.y = p.y+y_offset[i];
            if(q.x>I.cols-1) q.x = I.cols-1;
            if(q.x<0)       q.x = 0;
            if(q.y>I.rows-1) q.y = I.rows-1;
            if(q.y<0)       q.y = 0;
            
            if(J.at<uchar>(p.y,p.x)<J.at<uchar>(q.y,q.x) &&
                    I.at<uchar>(q.y,q.x)!=I.at<uchar>(q.y,q.x))
            {
                    J.at<uchar>(p.y,p.x)=min(J.at<uchar>(p.y,p.x),
                                                I.at<uchar>(q.y,q.x));
                    que.push(p);
            }
        }
    }
    return (I-J)>0;
}

实现效果如下：

由图像可见峰值点可以找出来，但是不一定为质心，因为连接处也可能存在峰值，这里需要后续的优化。需要注意的是峰值点不是单点，可能出现连续或间断点；笔者认为后续使用区域生长、主动轮廓模型可以完成实现分割。

结语

形态学灰度重建提供了一种峰值查找的解决方法，论文比较长，我也看得不太懂，由于前阶段毕业设计急着需要用就直接按照后面的伪代码实现了，后续有时间再将毕设后续的内容总结上来。

另外：该算法的并行算法实现详见我的简书：OpenCV CUDA 图像峰值查找，我总觉得搬过来不太好。

本文源码下载：01_findpeeks

参考文献：

Morphological Grayscale Reconstruction in Image Analysis: Applications and Effcient Algorithms