• Create Account

Banner advertising on our site currently available from just \$5!

# Zouflain

Member Since 12 Jan 2007
Offline Last Active Jan 09 2015 09:45 AM

### [C++] Can't find error with Delaunay Triangulation

31 December 2014 - 07:49 AM

I'm more rusty than I care to admit, both in the math and programming department. Coming back from a very long hiatus, I'm having serious trouble identifying my mistake in implementing the Bowyer-Watson algorithm.

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);
}
#endif
```
Delaunay.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;

const int INFINITY = 10000;
const Triangle super_triangle(Point(-INFINITY,INFINITY),Point(0,-INFINITY),Point(INFINITY,INFINITY));
triangles.push_back(super_triangle);

for(auto& point : points){

//Find non-delaunay triangle
for(auto& triangle : triangles){
if(triangle.circumcircles(point))
}

//Get good edges from bad triangles
EdgeList good_edges;
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;
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);
}
}
}

}

//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.

### [MinGW] Permanent Msys environment variables

30 December 2014 - 09:23 PM

Eons ago when I did windows programming I used Microsoft's Visual Express, but then I migrated to Linux for a few years and fell in love with GCC and makefiles. It's a much more sensible, portable toolchain than Visual Express and it's hard to give that up.

Well I've since started working on a Windows system for various reasons and I'm fighting with MinGW/Msys over environment variables.

I have a clear directory for my own includes/libs (namely G:/MinGW/usr/local/) and I can get a project to compile if I export the CPLUS_INCLUDE_PATH and LIBRARY_PATH, but this is an irritating thing to have to repeatedly configure every time I open a MSYS shell. Is there any way to make these environment variables permanent?

I've already tried hardcoding GCC's specs file, but despite following the MinGW instructions, GCC seems to ignore the specs file. I don't really know another workaround.

### Seeking Advice for Scripting Objects with Data Oriented Design

17 February 2014 - 03:36 PM

Ordinarily, when performing game logic on instances of a generic Entity class, after spacial partitioning and assembling a list of nearby, relevant entities, I would simply toss the list of entities and the relevant entities into a script (usually lua) and do something like:
```
--Below: A wolf hunts for rabbits to eat

TYPES = { WOLF=0,RABBIT=1}

...

function UpdateEntities(entities_to_update,relevant_entities)

for entity in entities_to_update do

if GetEntityType(entity)==TYPES.WOLF then

DoWolfAi(entity,relevant_entities);

end

end

end

ACTIONS = {HUNT=0,FLEE=1,PRAY=1029};

function DoWolfAi(entity,relevant_entities)

--Below: look for food

local actions = {};

for other in relevant_entities do

if other~=entity then

if GetEntityType(other)==TYPES.RABBIT then

table.insert(a,{action=ACTIONS.HUNT,target=other,priority=GetEntityHunger(entity)*GetDistanceBetween(entity,other)});

end

...

end

end

--Ommited: sort actions by member.priority

--Ommited: pick top action (hunt) and do hunt function

end

```

It's pretty straightforward, and works just fine for a small number of instances. In fact, it's served me well for the admitedly hackish projects I've thrown together. But, I'm moving on to larger projects and running into issues of scalability. After getting a massive wakeup call to DoD (Data Oriented Design) and things like cache misses (which didn't matter so much when it was 50 or so objects with simple scripts, but really matter when it's 100k complex objects), I've come to realize that my old way of doing things is grossly insufficient.

The way this handles data seems attrocious, with what I know now. First there's all those "GetBlahBlah(entity)" calls, which usually refer back to an map (std::map or boost::container::flat_map) to fetch data from a generic Entity class. Each one is a lot of overhead, because it has to look up the map to find the Entity, then find Entity.BlahBlah, then pass it all the way back to the script. I shudder to think what's happening to the cache during all this.

That wont work for hundreds or worse thousands of objects. It's also hard to envision a model that works in a cache-friendly manner. Each "wolf" entity needs to fetch data from every nearby "rabbit" entity (namely its location and type), and thus it's dependent on other instances. I also wont know which wolf will need which rabbit when the wolf and rabbit are allocated, so even though flat_map has contiguous memory, I'm forced to use random access rather than sequential (read: cache friendly).

My current project has a bit of a challenge mixed in: It has to run decently on a specific machine with an Intel Atom N570 (read: crap). 4 cores, 1.42ghz, ancient technology. That means I have to dance carefully.

So, are there any examples of a DoD (Data Oriented Design) approach to the same issue? Every book, every paper, every article I've read uses something more complex but fundamentally similar to the above. From gamedev to gamasutra, very little seems to use Data Oriented Design. I could really use some reference material, or just informative discussion. Thank you.

Note: this forum (and only this forum) likes to erase whitspaces and newlines. Not sure why, and it seems unique to me. If it's illegible I appologize, I'll try to edit and fix it but the edits often times erase more...

13 February 2014 - 09:23 PM

I'm implementing a quadtree for the spacial partitioning of a very large group of objects (target is 400k and trying to get it as fast as I can (target is no more than 50ms for all 400k). My initial naive approach was only able to handle about 1,000 objects in 27 seconds, but I've since been able to get it down to 10,000 objects in 113ms mostly by paying close attention to unnecessary move, copy, and temporary operations. I'm comfortably certain that there are no more such operations that can be eliminated.

The quadtree is a structure with two std::unordered_set: the first contains the id#'s of all objects inside the tree node aka "owned" objects, and the second contains the first plus any object close enough to the node's border to be included aka "known" objects. It also has a bounding rectangle defined by its top-left and bottom-right points (can be instantiated Rect(x1,y1,x2,y1) to avoid unecessary Point() constructors).

At every update, every known object is inserted into the quadtree. The quadtree is updated with quadtree::insert(const Point&amp;amp;amp; where,const uint64_t&amp;amp;amp; which).
--The method checks first if the quadtree has any children AND if it's depth is &amp;amp;lt; a maximum AND if it has &amp;amp;gt; max objecs within it. If that's true, it creates 4 children and dumps both of it's id containers into the children.
--The method then checks if it has any children - if so, call quadtree::insert for each child. If not, then check if the bounding rect contains "where"
----If the bounding rect contains "where" insert "which" into both containers
----If not, check if the bounding rect collides with a rect formed around where Rect(where.x-width/2,where.y-width/2,where.x+width/2,where.y+width/2). If so, add it to just the "known" container.

Unfortunately, after profiling 10 iterations of 10,000 randomly distributed objects, this leads to about ~408,000 calls to quadtree::insert (~37% execution time), 3,570,000 rect constructors (~6.6%), 7,256,000 point constructors (~10.38%), 3,586,000 stl emplacements (~13.21%).

The only really meaningful way to further optimize this that I can see would be to find a less naive way of splitting the quadtree, so that the top node needn't dump every "known" object into every child, which is what is dramatically multiplying the number of constructors. Any advice in this regard would be seriously appreciated.

Barely Annotated Code:
```
bool LogicTree::insert(const QuadPoint&amp;amp; where,const EntityID&amp;amp; which){

bool within = false;

//check if split needed

if(children.size()==0 &amp;amp;&amp;amp; depth &amp;lt; MAX_DEPTH &amp;amp;&amp;amp; owned.size()&amp;gt;MAX_ENTITIES){

Coord width = bounds.rightBottom.x-bounds.leftTop.x,height = bounds.rightBottom.y-bounds.leftTop.y;

c2(bounds.leftTop.x+width/2,bounds.leftTop.y,bounds.leftTop.x+width,bounds.leftTop.y+height/2),

c3(bounds.leftTop.x,bounds.leftTop.y+height/2,bounds.leftTop.x+width/2,bounds.leftTop.y+height),

c4(bounds.leftTop.x+width/2,bounds.leftTop.y+height/2,bounds.rightBottom.x,bounds.rightBottom.y);

children.emplace_back(new LogicTree(c1,depth+1));

children.emplace_back(new LogicTree(c2,depth+1));

children.emplace_back(new LogicTree(c3,depth+1));

children.emplace_back(new LogicTree(c4,depth+1));

for(auto&amp;amp; child : children){

massInsertChild(child);

}

}

//Pass to children or insert

if(children.size()&amp;gt;0){

for(auto&amp;amp; child : children){

if(child-&amp;gt;insert(where,which)){

within = true;

}

}

}else{

if(bounds.contains(where)){

owned.emplace(which);

known.emplace(which);

within = true;

known.emplace(which);

}

}

return within;

}

```
P.S&amp;gt;If this comes out as a solid block of text, my appologies - for some reason this forum likes to delete my newline characters.

Edit: Attempt #1 to fix the whitespace nuke. Also fixed a mistake (400k not 3mil, lol).