OpenCV4入门教程048：距离变换

索引地址：系列索引

距离变换于1966年被学者首次提出,目前已被广泛应用于图像分析、计算机视觉、模式识别等领域，人们利用它来实现目标细化、骨架提取、形状插值及匹配、粘连物体的分离等。距离变换是针对二值图像的一种变换。在二维空间中，一幅二值图像可以认为仅仅包含目标和背景两种像素，目标的像素值为1，背景的像素值为0；距离变换的结果不是另一幅二值图像，而是一幅灰度级图像，即距离图像，图像中每个像素的灰度值为该像素与距其最近的背景像素间的距离。

距离变换按照距离的类型可以分为欧式距离变换(Eudlidean Distance Transfrom)和非欧式距离变换两种，其中，非欧式距离变换又包括棋盘距离变换(Chessboard Distance Transform)，城市街区距离变换(Cityblock Distance Transform)，倒角距离变换(Chamfer Distance Transform)等；

距离变换的主要过程：

假设一幅二值图像I，包含一个连通区域S，其中有目标集O和背景集B，距离图为D，则距离变换的定义如公式：

$d(p)=Min(disf(p,q)) \ p\in O,q\in B$

具体步骤如下：

1，将图像中的目标像素点分类，分为内部点，外部点和孤立点。

以中心像素的四邻域为例，如果中心像素为目标像素(值为1)且四邻域都为目标像素(值为1)，则该点为内部点。如果该中心像素为目标像素，四邻域为背景像素(值为0)，则该中心点为孤立点。除了内部点和孤立点之外的目标区域点为边界点。

内部点

	1
1	1	1
	1

孤立点

	0
0	1	0
	0

2，计算图像中所有的内部点和非内部点，点集分别为S1，S2。

3，对于S1中的每一个内部点(x,y)，使用距离公式disf()计算骑在S2中的最小距离，这些最小距离构成集合S3。

4，计算S3中的最大最小值Max,Min。

5，对于每一个内部点，转换后的灰度值G计算如下公式所示：

$G(x,y)=\frac{255 * |S3(x,y)-Min|}{Max-Min}$

其中，S3(x,y)表示S1中的点(x,y)在S2中的最短距离

6，对于孤立点保持不变。

在以上距离变换的过程中，距离函数disf()的选取如果是欧式距离，则该距离变换称为欧式距离变换，依次类推。对于距离的求取，目前主要的距离公式如下：

欧式距离：

$disf(p(x_1,y_1),q(x_2,y_2))=\sqrt{(x_1-x_2)^2+(y_1-y_2)^2}$

棋盘距离：

$disf(p(x_1,y_1),q(x_2,y_2))=|x_1-x_2|+|y_1-y_2|$

城市街区距离：

$disf(p(x_1,y_1),q(x_2,y_2))=Max(|x_1-x_2|,|y_1-y_2|)$

对于欧式距离变换，由于其结果准确，而计算相比非欧式距离变换较为复杂，因此，出现了较多的快速欧式距离变换算法，这里笔者介绍一种基于 $3 * 3$ 模板的快速欧式距离变换算法(文献2)，具体过程如下：

1、按照从上到下，从左到右的顺序，使用模板如下，依次循环遍历图像I，此过程称为前向循环。

$\begin{array}{|c|c|c|} \hline & \text{快速欧式距离变换模板} & \\ \hline q_2 & q_3 & q_4 ~\\ \hline q_1 & p & q_5 \\ \hline q_8 & q_7 & q_6 \\ \hline \end{array}$

对于p对应的像素(x,y)，我们计算五个距离：d0,d1,d2,d3,d4：

$d0=p(x,y) \\ d1=p(x-1,y)+disf((x-1,y),(x,y)) \\ d2=p(x-1,y-1)+disf((x-1,y-1),(x,y)) \\ d3=p(x,y-1)+disf((x,y-1),(x,y)) \\ d4=p(x+1,y-1)+disf((x-1,y-1),(x,y))$

则p(x,y)变换后的像素值为：

$p(x,y)=Min(d0,d1,d2,d3,d4);$

使用上述算法得到图像I'。

2、按照从下到上，从右到左的顺序，使用上图所示模板依次循环遍历图像I’，此过程称为后向循环。

对于p对应的像素(x,y)，我们计算五个距离：d0,d5,d6,d7,d8：

$d0=p(x,y) \\ d5=p(x+1,y)+disf((x+1,y),(x,y)) \\ d6=p(x+1,y+1)+disf((x+1,y+1),(x,y)) \\ d7=p(x,y+1)+disf((x,y+1),(x,y)) \\ d8=p(x-1,y+1)+disf((x-1,y+1),(x,y))$

则p(x,y)后向变换后的像素值为：

$p(x,y)=Min(d0,d5,d6,d7,d8);$

使用上述算法得到的图像即为距离变换得到的灰度图像。

以上过程即文献2所示快速欧式距离变换算法。如果我们需要非欧氏距离变换的快速算法，只需要修改文献2算法中的欧式距离公式disf()为非欧式距离公式，如棋盘距离，城市街区距离等，过程依次类推。

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

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

OpenCV中函数原型为：

void cv::distanceTransform  (   InputArray      src,
        OutputArray     dst,
        int     distanceType,
        int     maskSize,
        int     dstType = CV_32F 
    );

参数说明：

src：8位单通道（二值）源图像
dst：保存了每一个点与最近的零点的距离信息，图像上越亮的点，代表了离零点的距离越远。是个8位或者32位浮点单通道与源图像同大小的图像
distanceType：距离类型，详见DistanceTypes
maskSize：距离变换掩码的大小，请参阅距离变换掩码。在DIST_L1或DIST_C距离类型的情况下，参数被强制为3，因为3×3掩码给出的结果与5或任何更大的孔径相同。
dstType：输出图像的类型。可以CV_8U也可以CV_32F。类型CV_8U只能用于函数的第一个变体和距离类型DIST_L1。

可以根据距离变换的这个性质，经过简单的运算，用于细化字符的轮廓和查找物体质心（中心）。

测试代码：

// 功能：代码 3-19 distanceTransform 距离变换
// 作者：朱伟 zhu1988wei@163.com
// 来源：《OpenCV图像处理编程实例》
// 博客：http://blog.csdn.net/zhuwei1988
// 更新：2016-8-1
// 说明：版权所有，引用或摘录请联系作者，并按照上面格式注明出处，谢谢。//

#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>

int main() {
    cv::Mat srcImage = cv::imread("lena.jpg");
    if (!srcImage.data)
        return -1;
    // 转换为灰度图像
    cv::Mat srcGray;
    cvtColor(srcImage, srcGray, cv::COLOR_BGR2GRAY);
    // 转换为二值图像
    cv::Mat srcBinary;
    threshold(srcGray, srcBinary, 160, 255, cv::THRESH_BINARY);
    // 距离变换
    cv::Mat dstImage;
    cv::distanceTransform(srcBinary, dstImage, cv::DIST_L2, cv::DIST_MASK_PRECISE);
    // 归一化矩阵
    cv::normalize(dstImage, dstImage, 0, 1., cv::NORM_MINMAX);
    cv::imshow("srcBinary", srcBinary);
    cv::imshow("dstImage", dstImage);
    cv::waitKey(0);

    return 0;
}

测试结果：

distance

动态调整参数测试代码：

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <stdio.h>
using namespace cv;

int maskSize0   = DIST_MASK_5;
int voronoiType = -1;
int edgeThresh  = 100;
int distType0   = DIST_L1;

// The output and temporary images
Mat gray;

// threshold trackbar callback
static void onTrackbar(int, void *) {
    static const Scalar colors[] = {Scalar(0, 0, 0),     Scalar(255, 0, 0), Scalar(255, 128, 0),
                                    Scalar(255, 255, 0), Scalar(0, 255, 0), Scalar(0, 128, 255),
                                    Scalar(0, 255, 255), Scalar(0, 0, 255), Scalar(255, 0, 255)};

    int maskSize = voronoiType >= 0 ? DIST_MASK_5 : maskSize0;
    int distType = voronoiType >= 0 ? DIST_L2 : distType0;

    Mat edge = gray >= edgeThresh, dist, labels, dist8u;

    if (voronoiType < 0)
        distanceTransform(edge, dist, distType, maskSize);
    else
        distanceTransform(edge, dist, labels, distType, maskSize, voronoiType);

    if (voronoiType < 0) {
        // begin "painting" the distance transform result
        dist *= 5000;
        pow(dist, 0.5, dist);

        Mat dist32s, dist8u1, dist8u2;

        dist.convertTo(dist32s, CV_32S, 1, 0.5);
        dist32s &= Scalar::all(255);

        dist32s.convertTo(dist8u1, CV_8U, 1, 0);
        dist32s *= -1;

        dist32s += Scalar::all(255);
        dist32s.convertTo(dist8u2, CV_8U);

        Mat planes[] = {dist8u1, dist8u2, dist8u2};
        merge(planes, 3, dist8u);
    } else {
        dist8u.create(labels.size(), CV_8UC3);
        for (int i = 0; i < labels.rows; i++) {
            const int *  ll = (const int *)labels.ptr(i);
            const float *dd = (const float *)dist.ptr(i);
            uchar *      d  = (uchar *)dist8u.ptr(i);
            for (int j = 0; j < labels.cols; j++) {
                int   idx      = ll[ j ] == 0 || dd[ j ] == 0 ? 0 : (ll[ j ] - 1) % 8 + 1;
                float scale    = 1.f / (1 + dd[ j ] * dd[ j ] * 0.0004f);
                int   b        = round(colors[ idx ][ 0 ] * scale);
                int   g        = round(colors[ idx ][ 1 ] * scale);
                int   r        = round(colors[ idx ][ 2 ] * scale);
                d[ j * 3 ]     = (uchar)b;
                d[ j * 3 + 1 ] = (uchar)g;
                d[ j * 3 + 2 ] = (uchar)r;
            }
        }
    }

    imshow("Distance Map", dist8u);
}

static void help() {
    printf("\n此程序用于演示边缘图像之间的距离变换。\n"
           "\n按键说明:\n"
           "\t【ESC】 -退出程序\n"
           "\t【c】-使用C/Inf度量\n"
           "\t【1】-使用L1度量\n"
           "\t【2】-使用L2度量\n"
           "\t【3】- 使用3 x 3的掩膜\n"
           "\t【5】- 使用5 x 5的掩膜\n"
           "\t【0】- 采用精确的距离变换\n"
           "\t【v】- 切换到Voronoi图（Voronoi diagram）模式\n"
           "\t【p】 - 切换到基于像素的Voronoi图模式\n"
           "\t【SPACE】- 在各种模式间切换\n\n");
}

const char *keys = {"{1| |stuff.jpg|input image file}"};

int main(int argc, const char **argv) {
    help();
    //  CommandLineParser parser(argc, argv, keys);
    //  string filename = parser.get<string>("1");
    gray = imread("lena.jpg", 0);
    //  if(gray.empty())
    //  {
    //      printf("Cannot read image file: %s\n", filename.c_str());
    //      help();
    //      return -1;
    //  }
    voronoiType = 0;
    namedWindow("Distance Map", 1);
    createTrackbar("Brightness Threshold", "Distance Map", &edgeThresh, 255, onTrackbar, 0);

    for (;;) {
        // Call to update the view
        onTrackbar(0, 0);

        int c = waitKey(0) & 255;

        if (c == 27)
            break;

        if (c == 'c' || c == 'C' || c == '1' || c == '2' || c == '3' || c == '5' || c == '0')
            voronoiType = -1;

        if (c == 'c' || c == 'C')
            distType0 = DIST_C;
        else if (c == '1')
            distType0 = DIST_L1;
        else if (c == '2')
            distType0 = DIST_L2;
        else if (c == '3')
            maskSize0 = DIST_MASK_3;
        else if (c == '5')
            maskSize0 = DIST_MASK_5;
        else if (c == '0')
            maskSize0 = DIST_MASK_PRECISE;
        else if (c == 'v')
            voronoiType = 0;
        else if (c == 'p')
            voronoiType = 1;
        else if (c == ' ') {
            if (voronoiType == 0)
                voronoiType = 1;
            else if (voronoiType == 1) {
                voronoiType = -1;
                maskSize0   = DIST_MASK_3;
                distType0   = DIST_C;
            } else if (distType0 == DIST_C)
                distType0 = DIST_L1;
            else if (distType0 == DIST_L1)
                distType0 = DIST_L2;
            else if (maskSize0 == DIST_MASK_3)
                maskSize0 = DIST_MASK_5;
            else if (maskSize0 == DIST_MASK_5)
                maskSize0 = DIST_MASK_PRECISE;
            else if (maskSize0 == DIST_MASK_PRECISE)
                voronoiType = 0;
        }
    }

    return 0;
}