【图像处理】——距离变换算法_图像操作 计算区域的距离变换

前言

通过本文可以了解到

• 什么是图像的距离？
• 什么是距离变换
• 距离变换的计算
• OpenCV中距离变换的实现

什么是图像的距离？

距离（distance）是描述图像两点像素之间的远近关系的度量，常见的度量距离有欧式距离（Euchildean distance）、城市街区距离（City block distance）、棋盘距离（Chessboard distance）。

• 欧式距离

欧式距离的定义源于经典的几何学，与我们数学中所学的简单几何的两点之间的距离一致,为两个像素点坐标值得平方根。欧式距离的优点在于其定义非常地直观，是显而易见的，但缺点在于平方根的计算是非常耗时的。

De = sqrt(((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1)));
1

• 城市街区距离

城市街区距离也称D4距离，此距离描述的是只允许像素坐标系平面中横向和纵向的移动距离，4表示在这种定义下，像素点是4邻接的，即每个点只与它的上、下、左、右相邻的4个点之间的距离为1。其定义为

D4 = abs(x1 - x2) + abs(y1 - y2);
1

• 棋盘距离

如果允许在图像坐标系中像素点的对角线方向的移动，就可以得到棋盘距离，也称D8距离，8表示在这种定义下，像素点是8邻接的，即每个点只与它的上、下、左、右、四个对角线方向相邻的8个点之间的距离为1。

D8 = max(abs(x1 - x2), (y1 - y2));
1

图1表示了三个距离之间的关系

​ 图1 三种距离示意图

什么是距离变换？

距离变换也叫作距离函数或者斜切算法。它是距离概念的一个应用，图像处理的一些算法以距离变换为基础。距离变换描述的是图像中像素点与某个区域块的距离，区域块中的像素点值为0，临近区域块的像素点的值较小，离它越远值越大。

图2是一个简单距离转换的例子，使用的棋盘距离

​ 图2 棋盘距离变换例子

距离变换算法步骤

直观计算

确定了距离公式之后，根据直观的理解应该怎么计算呢？还是拿上面的例子进行说明，B代表着黑色背景的集合，里面有b₁, b₂, b₃……等一些点；W代表着白色物体的集合，里面有w₁，w₂，w₃……等一些点。现在我们就是要计算W集合里面的每一个点，到B集合的距离，那么我们就每一次取一个W里面的w_i，然后根据距离公式，计算w_i到B集合每一个点b_j的距离D_ij，取其中距离D_ij 最小的值，作为最后的度量结果，直到W集合的每一个点都计算一遍。实际上，在进行计算距离这一步之前，还要遍历图像把像素点分成两个集合，这种直观的计算相对来说，较复杂且花费时间较长，因此前人发明一种快速计算距离变换的方法。

快速计算

对于欧式距离变换算法，相关学者研究了速度更快的倒角距离变换算法，来近似欧式距离变换的效果。具体过程如下：

• 使用前向模板如图3中左边3*3模板，对图像从上到下，从左到右进行扫描，模板中心0点对应的像素值如果为0则跳过，如果为1则计算模板中每个元素与其对应的像素值的和，分别为Sum1,Sum2,Sum3,Sum4，Sum5，而中心像素值为这五个和值中的最小值。
• 使用后向模板如图Fig.3中右边的3*3模板，对图像从下到上，从右到左进行扫描。
• 一般我们使用的倒角距离变换模板为3 * 3和5 * 5，分别如下图所示：

​ 图3 倒角距离变换模板

示例

AL AL
AL P          mask1
AL
   BR
P  BR         mask2
BR BR
1
2
3
4
5
6

• 将图像进行二值化，子图像值为0，背景为255；
• 利用Mask 1从左向右，从上到下扫描，p点是当前像素点，q点是Mask 1中AL邻域中的点，D为距离计算，包括棋盘距离、城市距离和欧式距离。Fp为p点的像素值，计算

F(p) = min( F(p),  F(q)+D(p,q) ), 其中，q属于AL.
1

• 再利用Mask 2从右向左，从下向上扫描，计算

F(p) = min( F(p),  F(q)+D(p,q) ), 其中，q属于BR.
1

代码实现距离变换算法

void imageprocess::DistanceTransform(const cv::Mat& src_image, cv::Mat& dst_image)
{
    //step1: check the input parameters: 检查输入参数
    assert(!src_image.empty());
    assert(src_image.channels() == 1);

    //step2: initialize dst_image : 初始化目标图像
    cv::threshold(src_image, dst_image, 100, 255, cv::THRESH_BINARY);


    //step3: pass throuth from top to bottom, left to right: 从上到下，从做到右遍历
    for (size_t i = 1; i < dst_image.rows - 1; i++)
    {
        for (size_t j = 1; j < dst_image.cols; j++)
        {
            //AL  AL    
            //AL  P
            //AL
            std::array<cv::Point, 4> AL;
            AL.at(0) =  cv::Point( i - 1, j - 1 );
            AL.at(1) =  cv::Point( i - 1, j );
            AL.at(2) =  cv::Point( i, j - 1 );
            AL.at(3) =  cv::Point(i + 1, j - 1 );

            int Fp = dst_image.at<uchar>(i, j);

            //Fq
            std::array<int, 4> Fq = { 0 };
            Fq.at(0) = dst_image.at<uchar>(i - 1, j - 1);
            Fq.at(1) = dst_image.at<uchar>(i - 1, j);
            Fq.at(2) = dst_image.at<uchar>(i, j - 1);
            Fq.at(3) = dst_image.at<uchar>(i + 1, j - 1);

            std::array<int, 4> Dpq = { 0 };
            std::array<int, 4> DpqAddFq = { 0 };
        
            for (size_t k = 0; k < 4; k++)
            {
                //D(p, q)
                Dpq.at(k) = D4(i, AL.at(k).x, j, AL.at(k).y);
                //D(p,q) + F(q)
                DpqAddFq.at(k) = Dpq.at(k) + Fq.at(k);
            }
            //F(p) = min[F(p), D(p,q) + F(q)]
            std::sort(DpqAddFq.begin(), DpqAddFq.end());
            
            auto min = DpqAddFq.at(0);
            Fp = std::min(Fp, min);

            dst_image.at<uchar>(i, j) = Fp;

        }
    }

    
    //step4: pass throuth from bottom to top, right to left： 从下到上，从右到左遍历

    for (int i = dst_image.rows - 2; i > 0; i--)
    {
        for (int j = dst_image.cols -  2; j >= 0 ; j--)
        {
            //        BR
            //  P   BR
            //  BR  BR
            std::array<cv::Point, 4> BR;
            BR.at(0) = cv::Point( i - 1, j + 1 );
            BR.at(1) = cv::Point( i , j + 1 );
            BR.at(2) = cv::Point( i + 1, j + 1 );
            BR.at(3) = cv::Point( i + 1, j );

            int Fp = dst_image.at<uchar>(i, j);

            //Fq
            std::array<int, 4> Fq = { 0 };
            Fq.at(0) = dst_image.at<uchar>(i - 1, j + 1);
            Fq.at(1) = dst_image.at<uchar>(i, j + 1);
            Fq.at(2) = dst_image.at<uchar>(i + 1, j + 1);
            Fq.at(3) = dst_image.at<uchar>(i + 1, j);

            std::array<int, 4> Dpq = { 0 };
            std::array<int, 4> DpqAddFq = { 0 };

            for (size_t k = 0; k < 4; k++)
            {
                //D(p, q)
                Dpq.at(k) = D4(i, BR.at(k).x, j, BR.at(k).y);
                //D(p,q) + F(q)
                DpqAddFq.at(k) = Dpq.at(k) + Fq.at(k);
            }

            //F(p) = min[F(p), D(p,q) + F(q)]
            std::sort(DpqAddFq.begin(), DpqAddFq.end());

            auto min = DpqAddFq.at(0);
            Fp = std::min(Fp, min);

            dst_image.at<uchar>(i, j) = static_cast<uchar>(Fp);

        }
    }
    
}

int imageprocess::D4(const int& x1, const int& x2, const int& y1, const int& y2)
{
    return abs(x1 - x2) + abs(y1 - y2);
}

int imageprocess::D8(const int& x1, const int& x2, const int& y1, const int& y2)
{
    return cv::max(abs(x1 - x2), (y1 - y2));
}
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
109
110
111
112

OpenCV distanceTransform函数

• 函数签名

/** @overload
@param src 8-bit, single-channel (binary) source image.
@param dst Output image with calculated distances. It is a 8-bit or 32-bit floating-point,
single-channel image of the same size as src .
@param distanceType Type of distance, see #DistanceTypes
@param maskSize Size of the distance transform mask, see #DistanceTransformMasks. In case of the
#DIST_L1 or #DIST_C distance type, the parameter is forced to 3 because a \f$3\times 3\f$ mask gives
the same result as \f$5\times 5\f$ or any larger aperture.
@param dstType Type of output image. It can be CV_8U or CV_32F. Type CV_8U can be used only for
the first variant of the function and distanceType == #DIST_L1.
*/
CV_EXPORTS_W void distanceTransform( InputArray src, OutputArray dst,
                                     int distanceType, int maskSize, int dstType=CV_32F);
1
2
3
4
5
6
7
8
9
10
11
12
13

参数详解：

• InputArray src：输入图像，一般为二值图像；
• OutputArray dst：输出的图像，距离变换结果，灰度图；
• int distanceType：用于距离变换的距离类型（欧氏距离：DIST_L2 = 2；$D_4$距离：DIST_L1 = 1；$D_8$距离：DIST_C = 3等）；
• int mask_size：距离变换掩模的大小，一般为3或5；
• int dstType：输出图像的数据类型，可以为CV_8U或CV_32F。

测试代码

#include<opencv2/opencv.hpp>
using namespace cv;
using namespace std;

int main()
{
    //初始化输入图像和变换结果图像
    Mat mat(480, 480, CV_8UC1, Scalar(0)), transMatE, transMatD4, transMatD8;

    //给输入图像指定三个像素点作为距离变换原点(区域块)
    mat.at<uchar>(100, 200) = 1;
    mat.at<uchar>(200, 100) = 1;
    mat.at<uchar>(300, 300) = 1;

    //将将输入图像中1和0调换，使得原点距离为0
    mat = 1 - mat;

    //显示原始图像(显示为黑色)
    imshow("原始图片", mat);

    //分别利用欧式距离、D4距离和D8距离作距离变换，将结果存入transMatD4、transMatD8和transMatE
    distanceTransform(mat, transMatE, DIST_L2, 0);
    distanceTransform(mat, transMatD4, DIST_L1, 0, CV_8U);
    distanceTransform(mat, transMatD8, DIST_C, 0);

    //欧式距离与D8距离作变换后，值为32位浮点数，以下代码将其值转为uchar类型
    transMatE.convertTo(transMatE, CV_8U);
    transMatD8.convertTo(transMatD8, CV_8U);

    //显示距离变换结果
    imshow("欧式距离变换后的图片", transMatE);
    imshow("D4距离变换后的图片", transMatD4);
    imshow("D8距离变换后的图片", transMatD8);


    waitKey();

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

运行结果

参考资料

• https://mangoroom.cn/opencv/distance-transfer.html
• https://blog.csdn.net/qq_35535744/article/details/95247601
• https://zhuanlan.zhihu.com/p/38917770
• https://medium.com/on-coding/euclidean-distance-transform-d37e06958216
• https://blog.csdn.net/Trent1985/article/details/18081761