形态学图像处理

1 工具

1.1 灰度腐蚀和膨胀

当平坦结构元b的原点是(x,y)时，它在(x,y)处对图像f的灰度腐蚀定义为，图像f与b重合区域中的最小值。结构元b在位置(x,y)处对图像f的腐蚀写为：

类似地，当b的反射的原点是(x,y)时，平坦结构元b在(x,y)处对图像f的膨胀，定义为图像f与b的反射重合区域中的最大值，即

示例：灰度X射线图像分别使用半径为2像素的圆形结构元腐蚀、膨胀图像

Mat src = imread("./7.png", 0);
Mat erodeImg, dilateImg;
Mat element = getStructuringElement(MORPH_ELLIPSE, Size(4, 4));
morphologyEx(src, erodeImg, MORPH_ERODE, element);
Mat element1 = getStructuringElement(MORPH_ELLIPSE, Size(4, 4));
morphologyEx(src, dilateImg, MORPH_DILATE, element1);

1.2 灰度开运算和闭运算

灰度图像开运算和闭运算的公式，形式上与二值图像开运算和闭运算的公式相同。结构元b对图像f的开运算是

开运算照例首先用b腐蚀f然后用b膨胀得到结果。类似地，b对f的闭运算是

示例：灰度X射线图像分别使用a）半径为3像素的圆形结构元进行开运算，b）使用半径为5像素的结构元进行闭运算

Mat src = imread("./7.png", 0);
Mat openImg, closeImg;
Mat element = getStructuringElement(MORPH_ELLIPSE, Size(6, 6));
morphologyEx(src, openImg, MORPH_OPEN, element);
Mat element1 = getStructuringElement(MORPH_ELLIPSE, Size(10, 10));

2 灰度图像形态学算法

2.1 形态学平滑

天鹅星座环超新星图像使用半径为5像素的圆形结构元进行开运算+闭运算

Mat src = imread("./8.png", 0);
Mat openImg, closeImg;
Mat element = getStructuringElement(MORPH_ELLIPSE, Size(5, 5));
morphologyEx(src, openImg, MORPH_OPEN, element);
Mat element1 = getStructuringElement(MORPH_ELLIPSE, Size(5, 5));
morphologyEx(openImg, closeImg, MORPH_CLOSE, element1);

2.2 形态学梯度

脑部CT图像膨胀和腐蚀结合图像相减

Mat src = imread("./9.png", 0);
Mat erodeImg, dilateImg;
Mat element = getStructuringElement(MORPH_ELLIPSE, Size(4, 4));
morphologyEx(src, erodeImg, MORPH_ERODE, element);
Mat element1 = getStructuringElement(MORPH_ELLIPSE, Size(4, 4));
morphologyEx(src, dilateImg, MORPH_DILATE, element1);
Mat sub = dilateImg - erodeImg;

2.3 顶帽变换和底帽变换

图像相减结合开运算和闭运算，可以得到所谓的顶帽变换和底帽变换。灰度图像f的顶帽变换定义为f减去其开运算：

类似地，f的底帽变换定义为f的闭运算减去f：

这些变换的主要应用之一是，在开运算或闭运算中用一个结构元从图像中删除目标，而不是拟合将被删除的目标。然后，差运算得到一幅仅保留已删除分量的图像。顶帽变换用于暗背景上的亮目标，而底帽变换则用于亮背景上的暗目标。因此我们通常将这两个变换称为白顶帽变换和黑底帽变换。

使用顶帽变换校正阴影

Mat src = imread("./10.png", 0);
Mat binImg, tophatImg, blackhatImg, dst;
// 1. Otsu最优阈值处理方法
threshold(src, binImg, 127, 255, THRESH_OTSU);
 
// 2. 顶帽变换校正阴影后Otsu
Mat element = getStructuringElement(MORPH_CROSS, Size(41, 41));
morphologyEx(src, tophatImg, MORPH_TOPHAT, element);
threshold(tophatImg, dst, 127, 255, THRESH_OTSU);

2.4 粒度测定

粒度测定是指确定图像中颗粒的大小分布。由于颗粒急剧无法整齐的分开，因此采用逐个识别颗粒的方法来计算颗粒的属两个非常困难。形态学可间接估计颗粒的大小分布，而不需要识别和测量各个颗粒。对于比背景亮且形状规则的颗粒，这种方法是用逐渐增大的结构元对图像进行开运算。基本思想是，某个特殊大小的开运算会对包含类似大小颗粒的输入图像的那些区域产生最大影响。对于开运算得到的每幅图像，我们计算像素之和。这个和值称为表面区域，它随结构元的增大而减小，因为开运算会减小图像的亮特征。这一过程得到一个一维阵列，阵列中的每个元素都是对应大小的结构元开运算后的像素之和。为了强调两个连续开运算之间的变化，我们计算一维阵列中相邻两个元素的差。画出差值的图像，曲线中的峰值就会指明图像中主要大小颗粒的分布。

Mat src = imread("./11.png", 0);
Mat openImg, dst;
int tempGray = sum(src)[0];
for (size_t i = 1; i < 36; i=i+2)
{
    Mat element = getStructuringElement(MORPH_ELLIPSE, Size(i, i));
    morphologyEx(src, openImg, MORPH_OPEN, element);
    int sumGray = sum(openImg)[0];
    int d = tempGray - sumGray;
    cout << d << endl;
    tempGray = sumGray;
    if (i%10==5)
        imshow("openImg(ksize="+to_string(i)+")", openImg);
}

2.5 纹理分割

　　形态学的纹理分割是以纹理内容为基础，找到两个区域的边界，将图像分割为不同的区域。下图是一幅在亮背景上叠加了暗斑点的噪声图像。图像有两个纹理区域，左侧区域包括一些较小的斑点，而右侧区域包括一些较大的斑点。由于目标斑点比背景暗，可以用一个尺寸大于较小斑点的圆形结构元对图像进行闭运算，删除较小的斑点，就得到只有大斑点的图像。再用尺寸大于较大斑点的圆形结构元对图像进行开运算，可以删除较大斑点之间的亮间距，整个图像形成左侧亮色和右侧暗色两个区域。通过形态学梯度运算，就得到两个区域的边界。最后将形态学梯度获得的边界叠加到原图像上，就实现了左右两种不同纹理的区域分割。

纹理分割

Mat src = imread("./12.png", 0);
Mat closeImg, dst, openImg, gradImg;
Mat element = getStructuringElement(MORPH_ELLIPSE, Size(19, 19));
morphologyEx(src, closeImg, MORPH_CLOSE, element);
Mat element1 = getStructuringElement(MORPH_ELLIPSE, Size(61, 61));
morphologyEx(closeImg, openImg, MORPH_OPEN, element1);
Mat element2 = getStructuringElement(MORPH_RECT, Size(3, 3));
morphologyEx(openImg, gradImg, MORPH_GRADIENT, element2);
bitwise_or(gradImg, src, dst);

3二值图像形态学算法

3.1 边界提取

前景像素集合A的边界ß(A)可按如下方式得到：首先使用合适的结构元B腐蚀A，然后求A和腐蚀结果的差集。也就是说，

边界提取

Mat src = imread("./4.png",0);
Mat thr, dst, ero;
threshold(src, thr, 128, 255, THRESH_BINARY);
Mat kernel = getStructuringElement(MORPH_RECT, Size(3, 3));
morphologyEx(thr, ero, MORPH_ERODE, kernel);

3.2 孔洞填充

Mat I = imread("./2.png",0);
threshold(I, I, 128, 255, THRESH_BINARY);
Mat Ic,thr, dst, tmp, dilImg;
bitwise_not(I, Ic);
Mat mask = Mat::zeros(I.size(), CV_8UC1);
tmp = mask.clone();
mask.at<uchar>(0, 0) = 255;
Mat kernel = getStructuringElement(MORPH_CROSS, Size(3, 3));
int diff = -1;
while (diff!=0) {
    tmp = mask.clone();
    morphologyEx(mask, dilImg, MORPH_DILATE, kernel);
    bitwise_and(dilImg, Ic, mask);
    diff = sum(mask - tmp)[0];
}

3.3 提取连通分量

　　从二值图像中提取连通分量是自动图像分析的核心步骤。冈萨雷斯《数字图像处理（第四版）》提供了一种提取连通分量的形态学算法，构造一个元素为0的阵列X0，其中对应连通分量的像素值为1，采用迭代过程可以得到所有的连通分量：

该算法与约束膨胀孔洞填充的思路相同，使用条件膨胀来限制膨胀的增长，但用I代替Ic以寻找前景点。

　　对于内含多个连通分量的图像A，从仅为连通分量A1内部的某个像素B开始，用3*3的结构元不断进行膨胀。由于其它连通分量与A1之间至少有一条像素宽度的空隙，每次膨胀都不会产生位于其它连通区域内的点。用每次膨胀后的图像与原始图像A取交集，就把膨胀限制在A1内部。随着集合B的不断膨胀，B的区域不断生长，但又被限制在连通分量A1的内部，最终就会充满整个连通分量A1，从而实现对连通分量A1的提取。

Mat dilImg, tmp1, sub, I;
Mat src = imread("./5.png",0);
threshold(src, I, 200, 255, THRESH_BINARY);
Mat tmp = I.clone();
Mat markImg = Mat::zeros(I.size(), CV_8UC3);
Mat mask = Mat::zeros(I.size(), CV_8UC1);
Mat kernel = getStructuringElement(MORPH_RECT, Size(3, 3));
vector<vector<Point>> cons;
while (countNonZero(tmp) != 0) {
    vector<Point> idx;
    findNonZero(tmp, idx);
    mask.at<uchar>(idx[0].y, idx[0].x) = 255;
    tmp1 = mask.clone();
    sub = mask.clone();
    while (countNonZero(sub) != 0)
    {
        morphologyEx(mask, dilImg, MORPH_DILATE, kernel);
        bitwise_and(dilImg, I, mask);
        sub = mask - tmp1;   
        tmp1 = mask.clone();
     }
     tmp -= mask;
}

3.4 凸壳

当且仅当数字集合A的欧式凸壳只包含于属于A的数字点，该数字集合A是凸的。一种简单的可视化方法是，用直（连续的）欧几里得线段连接其边界点，如果只有前景点包含于由这些线段形成的集合，那么集合是凸的。令Bi,i=1,2,3,4表示下图中的4个结构元。凸壳程序由如下形态学公式实现：

其中，X0i=I。分别用4个结构元做击中-击不中变换直至收敛，再求并集就是A的凸壳。为了防止凸壳增长到超出保证凸性所需的最小尺寸，设定限制凸壳不超过集合A的垂直和水平尺寸。

3.5 细化

图像细化（Image Thinning），一般指二值图像的骨架化（Image Skeletonization）的一种操作运算。所谓的细化就是经过一层层的剥离，从原来的图中去掉一些点，但仍要保持原来的形状，直到得到图像的骨架。骨架，可以理解为图象的中轴。

Mat src = imread("./6.png", 0);
Mat edge = src.clone();
Mat element = getStructuringElement(MORPH_CROSS, Size(3, 3));
Mat dst = Mat::zeros(src.size(), CV_8UC1);
Mat openImg, tmp;
int n = -1;
while (countNonZero(edge) != 0) {
    morphologyEx(edge, openImg, MORPH_OPEN, element); // 开运算
    subtract(edge, openImg, tmp); // 获得骨架子集
    bitwise_or(dst, tmp, dst); // 将删除的像素添加到骨架图
    morphologyEx(edge, edge, MORPH_ERODE, element); // 腐蚀，用于下一次迭代
}

opencv实现

#include <opencv2/opencv.hpp>
#include <algorithm>
#include <iostream>
using namespace std;
using namespace cv;
 
void ThinOnce(Mat& pSrc, Mat& pDst, int flag) {
    int rows = pSrc.rows;
    int cols = pSrc.cols;
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            if (pSrc.at<float>(i, j) == 1.0f) {
                /// get 8 neighbors
                /// calculate C(p)
                int P2 = (int)pSrc.at<float>(i - 1, j);
                int P3 = (int)pSrc.at<float>(i - 1, j + 1);
                int P4 = (int)pSrc.at<float>(i, j + 1);
                int P5 = (int)pSrc.at<float>(i + 1, j + 1);
                int P6 = (int)pSrc.at<float>(i + 1, j);
                int P7 = (int)pSrc.at<float>(i + 1, j - 1);
                int P8 = (int)pSrc.at<float>(i, j - 1);
                int P9 = (int)pSrc.at<float>(i - 1, j - 1);
                int sum = P2 + P3 + P4 + P5 + P6 + P7 + P8 + P9;
                int arr[9]{ P2, P3, P4, P5, P6, P7, P8, P9, P2 };
                int count = 0;
                for (int i = 0; i < 8; i++) {
                    if (arr[i] == 0 && arr[i + 1] == 1)
                        count++;
                } 
                if (flag == 0 && sum < 7 && sum > 1 && P2 * P4 * P6 == 0 && P4 * P6 * P8 == 0 && count == 1)
                    pDst.at<float>(i, j) = 0.0f;                
                if (flag == 1 && sum < 7 && sum > 1 && P2 * P4 * P8 == 0 && P2 * P6 * P8 == 0 && count == 1)
                    pDst.at<float>(i, j) = 0.0f;
            }
        }
    }
}
void thin(Mat& src, Mat& dst) {
    bool bDone = false;
    int rows = src.rows;
    int cols = src.cols;
 
    /// pad source and dst
    Mat p_enlarged_src, p_enlarged_dst;
    copyMakeBorder(src, p_enlarged_src, 1, 1, 1, 1, BORDER_CONSTANT, Scalar(0));
    p_enlarged_src.convertTo(p_enlarged_src, CV_32F, 1/255.0);
    p_enlarged_dst = p_enlarged_src.clone();<br>
    Mat p_cmp = Mat::zeros(rows + 2, cols + 2, CV_8UC1);
    int iter = 0;
    while (bDone != true) {
        // sub-iteration
        int i = iter % 2;
        ThinOnce(p_enlarged_src, p_enlarged_dst, i);
        // compare
        compare(p_enlarged_src, p_enlarged_dst, p_cmp, CMP_EQ);
        // check
        int num_non_zero = countNonZero(p_cmp);
        if (num_non_zero == (rows + 2) * (cols + 2)) {
            bDone = true;
        }
        // copy
        p_enlarged_dst.copyTo(p_enlarged_src);
        iter++;
    }
    // copy result
    dst = p_enlarged_dst({ 1, 1, cols, rows }).clone();
}
 
int main(int argc, char* argv[])
{
    Mat src = imread("./6.png", 0);
    Mat dst;
    thin(src, dst);
    imshow("src", src);
    imshow("dst", dst);
    waitKey(0);
    return 0;
}