This topic is 2646 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

I have developed convex hull 3D algorithm using gift wrapping at work. What I'm currently trying to do is to make the algorithm always return the correct result.
The plan to do that is:
1) Transform the algorithm to use sign tests
2) Make sign tests return always correct result - no matter what is input

As the sign tests I'm using:
1) For points A, B, C, D (3D) return whether point D is laying on, below, above the oriented plane A, B, C
2) For points A, B, C (2D) return whether point C is laying on, on the left, on the right of oriented line A, B
3) For points A, B, C (2D or 3D) determine whether point C is laying on infinite line A, B

I'm not intending to use more sign tests than these. Maybe even sign tests under 3 are not needed.

The algorithm I have developed is not just returning the points on the convex hull, but complete topology faces. If more than 3 points are laying on the plane (create the face of convex hull) then not just triangles, but whole face topology is returned as 1 face. The algorithm takes N points in floating point precision (float in C++) as input at start.

However I have following problems transforming the algorithm to sign tests:
Having 3D points P1 to Pn, n > 3, laying on the plane (or almost on the plane in case the sign tests can sometimes still fail) I need some way of how to transfer (possibly project 3D points to 2D points) in order to compute the 2D convex hull - just forming 1 face of 3D convex hull.
There are various ways of how to accomplish this I know, but none which I find quite satisfying considering numerical errors:
1) Find view point V (in 3D) which is not laying on the plane consisting of points P1 to Pn - then simulate all sign tests in 2D using the sign tests in 3D including point V - asking for whether point C is on, on the left or on the right of oriented line A, B is done so that it's asked whether point C is on, below, above the oriented plane A, B, V
2) Compute the normal N of plane consisting of points P1 to Pn - then simulate all sign tests in 2D using the sign tests in 3D where the N normal (3D vector) is used as one row in the matrix - variant of sign test 1) with just 3 vectors forming parallel-piped
3) Compute the normal N of plane consisting of points P1 to Pn - use N to denote 2 tangent vectors on plane - U, V and project the points to 2D using U, V
4) Compute the normal N of plane consisting of points P1 to Pn - use highest absolute value coordinate of normalized N to denote the coordinate to omit, omit that coordinate from 3D points to project them to 2D points

I'm currently using the method 3) as this is most straightforward.
However all of these methods are having various problems when you consider the numerical stability of the algorithm.
There are 3 criteria which I use to consider the method:
C1) Is the method correct
C2) Is the method introducing some numerical problems
C3) Is the method still going to work in case the points P1 to Pn are not exactly on the plane (due to numerical round-offs in sign tests - this is usually that points are returned to be on plane instead of being below or above plane)

For the method 1) there is basically problem of how to find good point V - one way is to compute the normal N and add the normal to arbitrary (e.g. first) chosen point from P1 to Pn. This is having the problem of computing normal N - badly numerically stable (I will explain below). Another methods exist - e.g. computing AA bounding box and randomly choosing the point V inside it - using the sign test 1) then test that point is not laying on the plane. This is omitting the computation of normal however the found point can still be quite near the plane just making it hard if the sign tests are not always absolutely correct - violating C3).
Or just randomly choose V and test whether it's not laying on the plane - same problem - the sign tests can fail quite easily - violating C3)

Question 1) Any other suggestions here for method of how to find point V?

For method 2) The method 2) is having the advantage that instead of point V, the normal N is used as sort of view direction and therefore can be chosen so that the algorithm is quite numerically stable against condition of C3). However the problem here is the computation of normal N - as this basically can't be done by sign tests - and therefore violating C2). E.g. imagine that P1 to Pn are having all coordinates in denormalized floating point (basically in like 10^-30 or even less). Then you are going to end up with normal (0, 0, 0) quite easily and it's just not much numerically stable.

Question 2) Can it be done so that it's numerically stable in all cases? How? (Considering floating points underflow and overflow)

For method 3) This is having same problem as 2) + projection is another problem as it requires another dot product

For method 4) This is most promising as one coordinate is omitted. However if we wanna to avoid the computation of normal N. Then the problem is how to choose which coordinate (x, y or z) to omit?

Question 3) Is there a way of how to find which coordinate to omit in some more elegant way? Also considering the case when P1 to Pn may not lay exactly on plane.

The another problem is how to make the algorithm to at least work and return some approximately correct output in case the sign tests are not going to be always right due to numerical round-offs.

The trick I'm currently using with sign tests is casting all coordinates from floats to doubles before evaluating the sign test - this is basically able to prevent some underflows and some overflows although still obviously not completely safe method for evaluating sign tests. However this method is fast to say and uses HW float and double support well.

The method I'm using to tackle the sign test errors is edge stack tracking. Basically every oriented edge (A,B) can get to the stack exactly only once. When the algorithm is working correctly, each oriented edge being the part of result comes in pair of (A, B) and (B, A) - in such case the solution is at least correctly closed - although it can still omit some point from result which should be in result. When the algorithm errors in sign tests then the oriented edge can come more than once and it results in some slightly bad result - through the tracking of edges the algorithm will however end everytime - without the edge tracking the algorithm can end in infinite loop.

Question 4) Any other possible ways of how to make algorithm more stable and more correct? E.g. am I able to guarantee to return always correctly closed hull faces as result even when the sign tests error (for N > 4 number of points)?

Thanks for any suggestions as I'm aware that this is difficult topic.

Also please do not write me to use doubles instead of floats as this is not an argument to make the algorithm really stable - any finite precision floating point representation basically has the round-off, overflow, underflow problem. While infinite floating point representations are slow for computation in practice, so this is also not an option - it's maybe just option for evaluating the sign tests, but nothing else.

##### Share on other sites
We implemented a robust convex hull in our open source Bullet physics library too. There is an option to create a convex hull of a point cloud, used for the separating axis test and contact clipping.

In a nutshell, we first compute the 3D convex hull, using integer arithmetic to avoid numerical issues. This can still generate coplanar faces. The algorithms we use (separating axis and contact clipping) can deal with those coplanar faces, so we could use this and have a robust solution. This implementation is used in the commercial 3D modeler Cinema 3D.

For optimization reasons, I added a second pass that searches coplanar faces and merges them using a 2d convex hull algorithm (Graham scan). To project the 3D faces in the 2D plane, I compute the shortest arc from the xy-plane normal to the face normal, and use that to rotate the 3D vertices into the xy plane, and manually clear the Z coordinate. This second pass is not perfectly robust, so it could create small gaps in the convex hull.

Do you have some unit tests with hard/degenerate cases that you can share?
Thanks,
Erwin

##### Share on other sites

We implemented a robust convex hull in our open source Bullet physics library too. There is an option to create a convex hull of a point cloud, used for the separating axis test and contact clipping.

In a nutshell, we first compute the 3D convex hull, using integer arithmetic to avoid numerical issues. This can still generate coplanar faces. The algorithms we use (separating axis and contact clipping) can deal with those coplanar faces, so we could use this and have a robust solution. This implementation is used in the commercial 3D modeler Cinema 3D.

For optimization reasons, I added a second pass that searches coplanar faces and merges them using a 2d convex hull algorithm (Graham scan). To project the 3D faces in the 2D plane, I compute the shortest arc from the xy-plane normal to the face normal, and use that to rotate the 3D vertices into the xy plane, and manually clear the Z coordinate. This second pass is not perfectly robust, so it could create small gaps in the convex hull.

Do you have some unit tests with hard/degenerate cases that you can share?
Thanks,
Erwin

Hmm, interesting idea to scale the data to some limited domain and then convert it to ints.

I'm using the convex hull currently just for generating the ISO surface faces. That means little points like up to 40 points is sent to convex hull 3D algorithm (gift wrapping).
Basically I have meshes generated by cut-curtisian mesher approach - means a lot of faces parallel to basic axes planes.
See pictures:
[attachment=5422:IsoSurfaceByConvexHull.png]
[attachment=5423:IsoSurfaceByConvexHull2.png]

First I find the ISO value nodes on the edges it goes through on the element - one hopefully closed 3D polyhedron. Then I add all points of element being bigger than currently used ISO value.
And from these points I contruct the convex hull and its faces. Then I go through all the faces and add only faces to solution having all nodes entirely consisting of ISO value nodes.
This process is because this way I separate the plus space (having values bigger than ISO value) and minus space (having values smaller than ISO value).

Particularly I don't know whether the similar approach you are using is useful to me as the ISO value nodes can be sometimes very near the actual nodes of element.
These is some risk that after initial scaling, the nodes will end up in same integer specified location.

Probably the best thing about this integer rescaling solution is however that the result is always correctly enclosed 2-manifold convex hull, which I can't tell about my current algorithm.

>> Do you have some unit tests with hard/degenerate cases that you can share?

I'm tracking the success of my algorithm by checking the correctness of edge usage - triple or more usage means some problem. I can possibly save the floating point coordinates of nodes in case the algorithm reports failed integrity test (the algorithm always ends with some more or less correct results) - particularly in most cases this will be some floating points cases with some points coplanar or almost coplanar. Will this be of some use to you?

1. 1
2. 2
Rutin
16
3. 3
4. 4
5. 5

• 26
• 11
• 9
• 9
• 11
• ### Forum Statistics

• Total Topics
633702
• Total Posts
3013452
×