How to do Delaunay Triangulation and return an adjacency matrix?

edit

I am trying to do Delaunay Triangulation for a set of points in OpenCV, but encountered a problem.

The function takes a matrix of coordinates and return an adjacency matrix. (If there is and edge connecting the point i and the point j, then adj(i,j) = 1, otherwise 0.)

I didn't get it working. The code below give strange results.

Could you please help?

An example of Delaunay Triangulation is given here.

Thank you in advance.

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