距离变换及其数学推导

距离变换

假设二值图像I，其中有前景点集F，背景点集B。则由通过距离变换可以得到距离图D。
其中

$$D(p)=min(d(p,q)),p\in F,q\in B$$
由这个简单的想法我们可以得到距离变换的实现代码

double ComputeEuDistance(Point2i x1, Point2i x2)
{
    return norm(x1 - x2);
}
Mat BruteForceDistanceTransform(Mat binary_img)
{
    vector<Point2i> foreth;
    vector<Point2i> background;

    for (size_t i = 0; i < binary_img.rows; i++)
    {
        for (size_t j = 0; j < binary_img.cols; j++)
        {
            if (binary_img.at<uchar>(j, i))
                foreth.emplace_back(j, i);
            else
                background.emplace_back(j, i);
        }
    }
    Mat result = Mat::zeros(binary_img.rows, binary_img.cols, CV_8U);
    for (auto& pf : foreth)
    {
        double mind = INFINITY;
        for (auto& bf : background)
        {
            auto d = ComputeEuDistance(pf, bf);
            mind = min(mind, d);
        }
        result.at<uchar>(pf) = mind;
    }
    double minimal, maximum;
    minMaxLoc(result, &minimal, &maximum);
    result = (result - minimal) * (255.0 / (maximum - minimal));
    return result;
}

L1-norm和L∞-norm下的加速

显然，上述算法时间复杂度过高，消耗的时间往往是不能接受的。

实际实现的时候，我们可以利用图像在空间上的连续性和在存储上的离散性进行加速。

在1-范数和∞-范数下（也就是所谓的街区距离和棋盘距离），我们可以通过相邻像素点的空间关系得到线性时间复杂度的算法。

1-范数下，我们可以通过对图像进行一次从左上到右下的遍历和一次相反方向的遍历来计算距离变换。第一次遍历的时候对当前前景点p，和它左边与上边的点$q_0,q_1$，我们有D(p)=min(D(q0),D(q1))+1。

第二次遍历的时候对当前前景点p，和它右边与下边的点$q_2,q_3$，我们有D(p) = min(D(p),min(D(q2),D(q3))+1)。

由此得到如下代码：

Mat L1DistanceTransform(Mat binary_img)
{
    Mat result = binary_img.clone()/255;
    for (size_t i = 1; i < result.rows; i++)
    {
        for (size_t j = 1; j < result.cols; j++)
        {
            if (result.at<uchar>(j, i))
            {
                result.at<uchar>(j, i) = min(result.at<uchar>(j - 1, i), result.at<uchar>(j, i - 1)) + 1;
            }
        }
    }
    for (int i = result.rows - 2; i >=0; i--)
    {
        for (int j = result.cols - 2; j >=0; j--)
        {
            if (result.at<uchar>(j, i))
            {
                result.at<uchar>(j, i) = min((int)result.at<uchar>(j, i),
                    (min(result.at<uchar>(j + 1, i), result.at<uchar>(j, i + 1))+1) );
            }
        }
    }
    double minimal, maximum;
    minMaxLoc(result, &minimal, &maximum);
    result = (result - minimal) * (255.0 / (maximum - minimal));
    return result;
}

对于$L\infty$-norm可以用类似的方法得到距离变换。在此就不展开。

L2-norm下的加速

使用棋盘距离和街区距离可以得到线性复杂度的算法，然而问题是这两种距离都没有旋转不变性，也不符合人类直觉。在许多场合欧氏距离进行变换可以得到更好的效果。

好在使用类似的思路，我们也可以对欧氏距离下的距离变换进行加速。假设前景一点p和它周围的一点q，假设已经算出了q的距离变换且离q最近的点为r，d(p,r)>d(q,r)，我们有$D(p)^2 = d(p,r)^2=d(q,r)^2+\Delta$

于是有

$$D(p)^2 = \min_i d(q_i,r)^2+\Delta$$

其中

$$\Delta =\begin{cases} 2*(q-r)_x+1&\text{q在p的左右邻域}\\ 2*(q-r)_y+1&\text{q在p的上下邻域}\\ 2*((q-r)_x+(q-r)_y+1)&\text{q在p的对角邻域} \end{cases}$$
同样的，通过两次遍历我们可以得到图像的距离变换，代码如下：

inline int hqp(const int Rx, const int Ry, const int i)
{
    static int G[][2] = { { 1,0 },{ 1,1 },{ 0,1 },{ 1,1 },{ 1,0 },{ 1,1 },{ 0,1 },{ 1,1 } };
    return G[i][0] * (2 * Rx + 1) + G[i][1] * (2 * Ry + 1);
}

//EDT
Mat EuclidienDistanceTransform(const Mat& binary_img)
{
    /*
    q2 q3 q4
    q1 p  q5
    q8 q7 q6
    */
    pair<Mat, Mat> R{ Mat::zeros(binary_img.rows, binary_img.cols, CV_32S) ,
        Mat::zeros(binary_img.rows, binary_img.cols, CV_32S) };
    Mat result = Mat::zeros(binary_img.rows, binary_img.cols,  CV_32S);
    static int G[][2] = { { 1,0 },{ 1,1 },{ 0,1 },{ 1,1 },{ 1,0 },{ 1,1 },{ 0,1 },{ 1,1 } };
    static int q[8][2] = { { -1,0 },{ -1, -1 },{ 0,-1 },{ 1,-1 },{ 1,0 },{ 1,1 },{ 0,1 },{ -1,1 } };
    int qindex[8];
    const int h = result.rows;
    const int w = result.cols;
    for (size_t i = 0; i < 8; i++)
    {
        qindex[i] = q[i][0] + q[i][1] * w;
    }
    int* result_ptr = result.ptr<int>();
    int* Rx_ptr = R.first.ptr<int>();
    int* Ry_ptr = R.second.ptr<int>();
    const unsigned char* binary_img_ptr = binary_img.ptr<uchar>();

    for (size_t i = 1; i < h - 1; i++)
    {
        int index = w*i;
        for (size_t j = 1; j < w - 1; j++)
        {
            index++;
            if (binary_img_ptr[index])
            {
                int& fp = result_ptr[index];
                fp = INT16_MAX;
                int& Rpx = Rx_ptr[index];
                int& Rpy = Ry_ptr[index];
                for (size_t k = 0; k < 4; k++)
                {
                    int absqindex = index + qindex[k];
                    auto Rx = Rx_ptr[absqindex];
                    auto Ry = Ry_ptr[absqindex];
                    auto hp = result_ptr[absqindex] + hqp(Rx, Ry, k);
                    if (hp< fp)
                    {
                        fp = hp;
                        Rpx = Rx + G[k][0];
                        Rpy = Ry + G[k][1];
                    }
                }
            }
        }
    }
    for (int i = h - 2; i >= 1; i--)
    {
        int index = w*i + w - 1;
        for (int j = w - 2; j >= 1; j--)
        {
            index--;
            if (binary_img_ptr[index])
            {
                int& fp = result_ptr[index];
                int& Rpx = Rx_ptr[index];
                int& Rpy = Ry_ptr[index];
                for (int k = 4; k < 8; k++)
                {
                    int absqindex = index + qindex[k];
                    auto Rx = Rx_ptr[absqindex];
                    auto Ry = Ry_ptr[absqindex];
                    auto hp = result_ptr[absqindex] + hqp(Rx, Ry, k);
                    if (hp < fp)
                    {
                        fp = hp;
                        Rpx = Rx + G[k][0];
                        Rpy = Ry + G[k][1];
                    }
                }
            }
        }
    }

    Mat result1;
    result.convertTo(result1, CV_32F);
    sqrt(result1, result1);
    double minimal, maximum;
    minMaxLoc(result1, &minimal, &maximum);
    //cout << minimal << " " << maximum << endl;
    result1 = (result1 - minimal) * (255.0 / (maximum - minimal));
    Mat rr;
    result1.convertTo(rr, CV_8U);
    return rr;
}