## Friday, December 10, 2010

### Computational Geometry Sucks!

Computational Geometry is hard. Most of the examples out there are crap and the good stuff is without exception hard to understand.

You can usually whip up a 100 liner to solve a problems, if your input is cloud of random points. This is what pretty much all the Java code out there does, but the input data for your game is not a random point cloud!

I hate computation geometry, but yet it is so important! It drives me crazy.

I was recently looking for a robust 3D convex hull algorithm which would work for a couple of hundred points in realtime. That is, something sub 2.0 ms on a 2GHz Core Duo.

I wrote a couple of solutions myself, including one which was O(n^4) and surprisingly did not scale beyond 25 points. Finally I was hinted towards this patch for Bullet Physics, and it turned out to be a real gem!

It is fast, scales well (it's O(n log h)), and is robust too. My litmus test is a grid of 4x4x4 points. It passes that and even handles coplanar input. Pretty much all the Java code out there fails on that input.

So I took that patch, regexp'd it, copy pasted some Bullet code, deleted stuff I did not need, shuffled things around a bit, and made it output triangles and I'm finally a happy man!

1. Been five years since I studied computational geometry, but if I remember, there were some very simple, elegant, fast randomized algorithms which sacrificed absolute complexity for expected complexity (i.e., runtime is O(n^4) but *expected* runtime depending on random processes is O(n log n)). Don't remember how robust they were w/r/t non-general position.

2. You're probably referring to the incremental algorithm. It is simple to implement, but does not handle coplanar faces well (i.e. That tesselated cube i mentioned).

3. Let me recommend CGAL as the go-to library for computational geometry. It's pretty heavyweight to get started, but well worth the effort. There might be licensing issues for you, check to make sure before you get carried away.

For doing geometric tests ("are these four points coplanar?") you might be interested in Shewchuk's robust geometric predicates -- Google "shewchuk predicates.c" If you're not interested, either you've thought extremely carefully and you have formal proofs that it's OK to be inexact, or you're writing bugs. You've already hit this fact with the tesselated cube.

Shewchuk's approach is pretty fast: compared to using doubles, it's a factor of two for the incircle test, and about 20% overall for a 2d Delaunay triangulation algorithm. You could also use an exact arithmetic package like GMP but then you lose an order of magnitude overall.

This makes is so that all your tests are exact; no more issues there! But of course when you do a construction ("add a point on the plane defined by these three points") it's still a tough life we live, with these computers that can't handle real numbers. When constructions are involved, you're stuck either praying or doing higher math.

For books, de Berg et al ("Computational Geometry") is great for 2d geometry. I don't have a great text to recommend for 3d.

Cheers!

4. Finally added this convex hull computation to the Bullet library:

The implementation has been tested in production, and is part of MAXON Cinema4D package.

5. Erwin, great! I'd like to make a simple standalone library from the code, but I have been swamped with my other projects.

6. Unlike CGAL, the Bullet library keeps dependencies low, this snippet depends only on a handful of files. Distibuting such gems as snippet/standalone is a good idea.

7. Needs a 100% public domain C version! Sadly, can't get there from here (if I tried to paraphrase the Bullet one, then it wouldn't be legimitate to public domainify it; and I don't have the time/energy/interest to figure it out from scratch).

8. Shewchuk predicates are not completely stable, at least I wasn't able to make it stable on new HW:

double pointA[3] = {0.12539999f, 0.0016452915f, -0.019413333f};
double pointB[3] = {0.12539999f, 0.0017933375f, -0.019214222f};
double pointC[3] = {0.12539999f, 0.0017933375f, -0.017919999f};
double pointD[3] = {0.11700000f, 0.0017933375f, -0.018710587f};
double pointE[3] = {0.12539999f, 0.0017933375f, -0.019413333f};

std::cout << "On plane: A, C, B, E: " << onPlane(pointA, pointC, pointB, pointE) << std::endl;
std::cout << "On plane: B, C, D, E: " << onPlane(pointB, pointC, pointD, pointE) << std::endl;
std::cout << "On plane: B, C, D, A: " << onPlane(pointB, pointC, pointD, pointA) << std::endl;

This three tests ultimately fail telling you that first 2 tests are true and the last one is false, which is logically impossible and therefore predicates are not working correctly on modern HW at least.

9. @Pavel - It took me ages to figure this out! Your points A, B, C and E are all exactly on a straight line, so your first onPlane test is not meaningful. The 'modern hardware' and precision issues raised in the Shewchuck paper are really not a problem on x86 or x64 processors, as long as the FPU precision is not changed from the 64-bit (53-bit significand) default. I'll have a carefully checked .NET port of the Shewchuck predicates posted soon - look for it on github: https://github.com/govert/RobustGeometry.NET.