OpenCV: 如何进行Delaunay三角剖分并返回邻接矩阵？

Question

我正在尝试使用OpenCV对一组点进行Delaunay三角剖分，但遇到了问题。

该函数接收一个坐标矩阵，并返回一个邻接矩阵。（如果有连接点i和点j的边，则adj（i，j）= 1，否则为0。）

我无法使其正常工作。以下代码给出奇怪的结果。

请问您能提供帮助吗？

这里提供了Delaunay Triangulation的示例：链接。

非常感谢您的帮助。

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

using namespace std;
using namespace cv;

Mat delaunay(const Mat& points, int imRows, int imCols)
/// Return the Delaunay triangulation, under the form of an adjacency matrix
/// points is a Nx2 mat containing the coordinates (x, y) of the points
{
    Mat adj(points.rows, points.rows, CV_32S, Scalar(0));

    /// Create subdiv and insert the points to it
    Subdiv2D subdiv(Rect(0,0,imCols,imRows));
    for(int p = 0; p < points.rows; p++)
    {
        float xp = points.at<float>(p, 0);
        float yp = points.at<float>(p, 1);
        Point2f fp(xp, yp);
        subdiv.insert(fp);
    }

    /// Get the number of edges
    vector<Vec4f> edgeList;
    subdiv.getEdgeList(edgeList);
    int nE = edgeList.size();

    /// Check adjacency
    for(int e = 1; e <= nE; e++)
    {
        int p = subdiv.edgeOrg(e); // Edge's origin
        int q = subdiv.edgeDst(e); // Edge's destination

        if(p < points.rows && q < points.rows)
            adj.at<int>(p, q) = 1;
//        else
//        {
//            cout<<p<<", "<<q<<endl;
//            assert(p < points.rows && q < points.rows);
//        }
    }
    return adj;
}



int main()
{
    Mat points = Mat(100, 2, CV_32F);
    randu(points, 0, 99);

    int rows = 100, cols = 100;
    Mat im(rows, cols, CV_8UC3, Scalar::all(0));
    Mat adj = delaunay(points, rows, cols);

    for(int i = 0; i < points.rows; i++)
    {
        int xi = points.at<float>(i,0); 
        int yi = points.at<float>(i,1);
        /// Draw the edges
        for(int j = i+1; j < points.rows; j++)
        {
            if(adj.at<int>(i,j) > 0)
            {
                int xj = points.at<float>(j,0); 
                int yj = points.at<float>(j,1);
                line(im, Point(xi,yi), Point(xj,yj), Scalar(255,0,0), 1);
            }
        /// Draw the nodes
        circle(im, Point(xi, yi), 1, Scalar(0,0,255), -1);
        }
    }
    namedWindow("im", CV_WINDOW_NORMAL);
    imshow("im",im);
    waitKey();
    return 0;
}

- f10w

Answer 1

你正在向邻接矩阵中插入Subdiv2d边缘的索引，这些索引不对应于点的索引。

你可以通过将点及其索引存储到std::map中来解决此问题。当从Subdiv2d检索边缘时，您需要检查该边缘是否由您的点而非Subdiv2d添加的边界点组成。通过存储点的索引，您现在能够正确构建邻接矩阵。

请查看以下代码：

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

struct lessPoint2f
{
    bool operator()(const Point2f& lhs, const Point2f& rhs) const
    {
        return (lhs.x == rhs.x) ? (lhs.y < rhs.y) : (lhs.x < rhs.x);
    }
};

Mat delaunay(const Mat1f& points, int imRows, int imCols)
/// Return the Delaunay triangulation, under the form of an adjacency matrix
/// points is a Nx2 mat containing the coordinates (x, y) of the points
{
    map<Point2f, int, lessPoint2f> mappts;

    Mat1b adj(points.rows, points.rows, uchar(0));

    /// Create subdiv and insert the points to it
    Subdiv2D subdiv(Rect(0, 0, imCols, imRows));
    for (int p = 0; p < points.rows; p++)
    {
        float xp = points(p, 0);
        float yp = points(p, 1);
        Point2f fp(xp, yp);

        // Don't add duplicates
        if (mappts.count(fp) == 0)
        {
            // Save point and index
            mappts[fp] = p;

            subdiv.insert(fp);
        }
    }

    /// Get the number of edges
    vector<Vec4f> edgeList;
    subdiv.getEdgeList(edgeList);
    int nE = edgeList.size();

    /// Check adjacency
    for (int i = 0; i < nE; i++)
    {
        Vec4f e = edgeList[i];
        Point2f pt0(e[0], e[1]);
        Point2f pt1(e[2], e[3]);

        if (mappts.count(pt0) == 0 || mappts.count(pt1) == 0) {
            // Not a valid point
            continue;
        }

        int idx0 = mappts[pt0];
        int idx1 = mappts[pt1];

        // Symmetric matrix
        adj(idx0, idx1) = 1;
        adj(idx1, idx0) = 1;
    }
    return adj;
}



int main()
{
    Mat1f points(10, 2);
    randu(points, 0, 99);

    int rows = 100, cols = 100;
    Mat3b im(rows, cols, Vec3b(0,0,0));
    Mat1b adj = delaunay(points, rows, cols);

    for (int i = 0; i < points.rows; i++)
    {
        int xi = points.at<float>(i, 0);
        int yi = points.at<float>(i, 1);

        /// Draw the edges
        for (int j = i + 1; j < points.rows; j++)
        {
            if (adj(i, j))
            {
                int xj = points(j, 0);
                int yj = points(j, 1);
                line(im, Point(xi, yi), Point(xj, yj), Scalar(255, 0, 0), 1);
            }
        }
    }

    for (int i = 0; i < points.rows; i++)
    {
        int xi = points(i, 0);
        int yi = points(i, 1);

        /// Draw the nodes
        circle(im, Point(xi, yi), 1, Scalar(0, 0, 255), -1);
    }

    imshow("im", im);
    waitKey();
    return 0;
}

- Miki