OpenCV:如何进行 Delaunay 三角测量并返回邻接矩阵

OpenCV: How to do Delaunay Triangulation and return an adjacency matrix?

本文关键字:测量 返回 邻接矩阵 三角 Delaunay 何进行 OpenCV      更新时间:2023-10-16

我正在尝试对OpenCV中的一组点进行Delaunay 三角测量,但遇到了一个问题。

该函数采用坐标矩阵并返回邻接矩阵。(如果有 和 边连接点 i 和点 j,则 adj(i,j( = 1,否则为 0。

我没有让它工作。下面的代码给出了奇怪的结果。

你能帮忙吗?

这里给出了德劳奈三角测量的一个例子。

提前谢谢你。

#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;
}

您将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;
}