To describe the issue, I'd have to say that somewhere along the line the wrong edges are being removed from the triangulation (aka "std::list<Triangle> triangles") and I cannot figure out why. In the attached example.png I highlighted two problem areas: in red you see an edge that clearly should not exist, and in mauve you see a point that ought to have edges but doesn't. The greenish lines are the remaining edges and the rectangles are centered about the points fed into the bowyers watson algorithm.

I tried to follow the algorithm as tightly as possible. Anyone willing to point me in the right direction would be greatly appreciated. I've been analyzing and reanalyzing line by line for 3 straight days to figure out what the problem is to no avail.

The code is below. I've tried to keep it as clean and readable as possible, but the algorithm is verbose. Most of it is boiler plate. The last function is the actual algorithm. I've included the boiler plate just in case the error is there.

Delaunay.hpp

#ifndef DELAUNAY_HPP #define DELAUNAY_HPP #include "types.hpp" #include <list> namespace Delaunay{ struct Point; struct Triangle; struct Edge; typedef std::list<Point> PointList; typedef std::list<Triangle> TriangleList; typedef std::list<Edge> EdgeList; struct Mat3i{ S64 matrix[9];//COLUMN MAJOR order inline S64 determinant(void) const{ return matrix[0]*(matrix[4]*matrix[8]-matrix[7]*matrix[5]) -matrix[3]*(matrix[1]*matrix[8]-matrix[7]*matrix[2]) +matrix[6]*(matrix[1]*matrix[5]-matrix[4]*matrix[2]); } inline S64& operator[](int n) { return matrix[n]; } Mat3i(S64 i_x,S64 j_x,S64 k_x,S64 i_y,S64 j_y,S64 k_y,S64 i_z,S64 j_z,S64 k_z);//ROW MAJOR input, for easy reading }; struct Mat4i{ S64 matrix[16];//COLUMN MAJOR order inline S64& operator[](int n) { return matrix[n]; } S64 determinant(void) const; Mat4i(S64 i_x,S64 j_x,S64 k_x,S64 l_x,S64 i_y,S64 j_y,S64 k_y,S64 l_y, S64 i_z,S64 j_z,S64 k_z,S64 l_z,S64 i_w,S64 j_w,S64 k_w,S64 l_w );//ROW MAJOR input }; struct Point{ S32 x,y,room_id; Point& operator=(const Point& point); inline bool operator==(const Point& point) const{ return(x==point.x&&y==point.y); } inline bool operator!=(const Point& point) const { return !(*this == point); } Point(const Point& point); Point(S32 x=0,S32 y=0,S32 room_id=-1); }; struct Edge{ Point points[2]; inline bool hasPoint(const Point& point) const{ return points[0]==point||points[1]==point; } Edge& operator=(const Edge& edge); inline bool operator==(const Edge& edge) const{ return hasPoint(edge.points[0])&&hasPoint(edge.points[1])&&points[0]!=points[1]; } Edge(const Edge& edge); Edge(const Point& p1=Point(),const Point& p2=Point()); }; struct Triangle{ Point points[3];//ALWAYS in counter clockwise order (after initialization) inline bool hasEdge(const Edge& edge) const{ return Edge(points[0],points[1])==edge||Edge(points[1],points[2])==edge||Edge(points[2],points[0])==edge; } inline bool hasPoint(const Point& point) const{ return points[0]==point||points[1]==point||points[2]==point; } inline bool sharesPoint(const Triangle& triangle) const{ return hasPoint(triangle.points[0])||hasPoint(triangle.points[1])||hasPoint(triangle.points[2]); } bool circumcircles(const Point& p) const; void makeCounterClockwise(void);//used in each initialization path, keeps points in CCW order Triangle& operator=(const Triangle& triangle); inline bool operator==(const Triangle& triangle) const{ return sharesPoint(triangle.points[0])&&sharesPoint(triangle.points[1])&&sharesPoint(triangle.points[2]); } inline bool operator!=(const Triangle& triangle) const{ return !(*this==triangle); } Triangle(const Triangle& triangle); Triangle(const Point& p1=Point(),const Point& p2=Point(),const Point& p3=Point()); Triangle(const Point& point,const Edge& edge); }; EdgeList bowyerWatson(const PointList& points); } #endifDelaunay.cpp

#include "delaunay.hpp" #include <cstdio> using namespace Delaunay; Delaunay::Mat3i::Mat3i(S64 i_x,S64 j_x,S64 k_x,S64 i_y,S64 j_y,S64 k_y,S64 i_z,S64 j_z,S64 k_z){ //Again this is ROW MAJOR input, COLUMN MAJOR storage matrix[0]=i_x;matrix[3]=j_x;matrix[6]=k_x; matrix[1]=i_y;matrix[4]=j_y;matrix[7]=k_y; matrix[2]=i_z;matrix[5]=j_z;matrix[8]=k_z; } S64 Delaunay::Mat4i::determinant(void) const{ Mat3i a( matrix[5],matrix[9],matrix[13], matrix[6],matrix[10],matrix[14], matrix[7],matrix[11],matrix[15]), b( matrix[1],matrix[9],matrix[13], matrix[2],matrix[10],matrix[14], matrix[3],matrix[11],matrix[15]), c( matrix[1],matrix[5],matrix[13], matrix[2],matrix[6],matrix[14], matrix[3],matrix[7],matrix[15]), d( matrix[1],matrix[5],matrix[9], matrix[2],matrix[6],matrix[10], matrix[3],matrix[7],matrix[11]); return matrix[0]*a.determinant()-matrix[4]*b.determinant()+matrix[8]*c.determinant()-matrix[12]*d.determinant(); } Delaunay::Mat4i::Mat4i(S64 i_x,S64 j_x,S64 k_x,S64 l_x,S64 i_y,S64 j_y,S64 k_y,S64 l_y, S64 i_z,S64 j_z,S64 k_z,S64 l_z,S64 i_w,S64 j_w,S64 k_w,S64 l_w){ //ROW MAJOR input, like mat3 matrix[0]=i_x; matrix[4]=j_x; matrix[8]=k_x; matrix[12]=l_x; matrix[1]=i_y; matrix[5]=j_y; matrix[9]=k_y; matrix[13]=l_y; matrix[2]=i_z; matrix[6]=j_z; matrix[10]=k_z; matrix[14]=l_z; matrix[3]=i_w; matrix[7]=j_w; matrix[11]=k_w; matrix[15]=l_w; } Point& Delaunay::Point::operator=(const Point& point){ room_id = point.room_id; x = point.x; y = point.y; return *this; } Delaunay::Point::Point(const Point& point){ x = point.x; y = point.y; room_id = point.room_id; } Delaunay::Point::Point(S32 x,S32 y,S32 room_id){ this->x = x; this->y = y; this->room_id = room_id; } Edge& Delaunay::Edge::operator=(const Edge& edge){ points[0] = edge.points[0]; points[1] = edge.points[1]; return *this; } Delaunay::Edge::Edge(const Edge& edge){ points[0] = edge.points[0]; points[1] = edge.points[1]; } Delaunay::Edge::Edge(const Point& p1,const Point& p2){ points[0] = p1; points[1] = p2; } bool Delaunay::Triangle::circumcircles(const Point& p) const{ Mat3i geometric_construction( points[0].x-p.x,points[0].y-p.y,(points[0].x*points[0].x-p.x*p.x)+(points[0].y*points[0].y+p.y*p.y), points[1].x-p.x,points[1].y-p.y,(points[1].x*points[1].x-p.x*p.x)+(points[1].y*points[1].y+p.y*p.y), points[2].x-p.x,points[2].y-p.y,(points[2].x*points[2].x-p.x*p.x)+(points[2].y*points[2].y+p.y*p.y) ); return geometric_construction.determinant()>0; }; void Delaunay::Triangle::makeCounterClockwise(void){ //Poor man's vector math S32 vx = points[0].x,vy = points[0].y, ux = points[1].x,uy = points[1].y, wx = points[2].x,wy = points[2].y; S32 vux = vx-ux,vuy = vy-uy, wvx = wx-vx,wvy = wy-vy; //Swap if necessary (dot product > 0) if(vux*wvy-vuy*wvx>0){ Point temp_points[3];//size 3 to make things clear, technically [0] never gets used temp_points[1] = points[2]; temp_points[2] = points[1]; points[1] = temp_points[1]; points[2] = temp_points[2]; } } Triangle& Delaunay::Triangle::operator=(const Triangle& triangle){ points[0] = triangle.points[0]; points[1] = triangle.points[1]; points[2] = triangle.points[2]; return *this; } Delaunay::Triangle::Triangle(const Triangle& triangle){ for(int i=0;i<3;i++){ points[i] = triangle.points[i]; } makeCounterClockwise(); } Delaunay::Triangle::Triangle(const Point& p1,const Point& p2,const Point& p3){ points[0] = p1; points[1] = p2; points[2] = p3; makeCounterClockwise(); } Delaunay::Triangle::Triangle(const Point& point,const Edge& edge){ points[0] = point; points[1] = edge.points[0]; points[2] = edge.points[1]; makeCounterClockwise(); } EdgeList Delaunay::bowyerWatson(const PointList& points){ TriangleList triangles; //add super triangle const int INFINITY = 10000; const Triangle super_triangle(Point(-INFINITY,INFINITY),Point(0,-INFINITY),Point(INFINITY,INFINITY)); triangles.push_back(super_triangle); //add the points incrementally for(auto& point : points){ //Find non-delaunay triangle TriangleList bad_triangles; for(auto& triangle : triangles){ if(triangle.circumcircles(point)) bad_triangles.push_back(triangle); } //Get good edges from bad triangles EdgeList good_edges; for(auto& triangle : bad_triangles){ Edge edges[] = { Edge(triangle.points[0],triangle.points[1]), Edge(triangle.points[1],triangle.points[2]), Edge(triangle.points[2],triangle.points[0]) }; for(auto& edge : edges){ bool keep_edge = true; for(auto& test_triangle : bad_triangles){ if(triangle!=test_triangle && test_triangle.hasEdge(edge)){ keep_edge = false; break;//no need to keep iterating } } if(keep_edge){ good_edges.push_back(edge); } } } //Toss bad triangles for(auto& bad_triangle : bad_triangles){ triangles.remove(bad_triangle); } //Build triangles from good edges for(auto& edge : good_edges){ triangles.push_back(Triangle(point,edge)); } } //Export edges that do not have a point on the super_triangle EdgeList final_edges; for(auto& triangle : triangles){ Edge edges[3] = { Edge(triangle.points[0],triangle.points[1]), Edge(triangle.points[1],triangle.points[2]), Edge(triangle.points[2],triangle.points[0]) }; for(auto& edge : edges){ if(!super_triangle.hasPoint(edge.points[0])&&!super_triangle.hasPoint(edge.points[1])){ final_edges.push_back(edge); } } } return final_edges; }Note: I'll happily supply more info for anyone who's interested in helping. I don't expect anyone to spend too much time wading through somebody else's code, but I included it because I'm not sure what additional information will be needed to find the problem. Thanks again for any pointers.