【图像处理】——距离变换算法

前言

通过本文可以了解到

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

什么是图像的距离？

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

欧式距离

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

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

城市街区距离

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

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

棋盘距离

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

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

图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
ALBR
P  BR         mask2
BR BR

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

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

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

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

代码实现距离变换算法

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//ALstd::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);//Fqstd::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  BRstd::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);//Fqstd::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));
}

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);

参数详解：

InputArray src：输入图像，一般为二值图像；
OutputArray dst：输出的图像，距离变换结果，灰度图；
int distanceType：用于距离变换的距离类型（欧氏距离：DIST_L2 = 2；D4D_4D4距离：DIST_L1 = 1；D8D_8D8距离：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调换，使得原点距离为0mat = 1 - mat;//显示原始图像(显示为黑色)imshow("原始图片", mat);//分别利用欧式距离、D4距离和D8距离作距离变换，将结果存入transMatD4、transMatD8和transMatEdistanceTransform(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;
}

运行结果

参考资料

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