Good catch, std10093. I'll investigate.
Edit 1: There was a typo when defining v3. The last component should be p3.z-p1.z, not p3.z-p3.z.
Edit 2: For the purposes of implementing this algorithm, we need to distinguish not only between "outside" and "inside" a half-space, but also "on the dividing plane". My code in #12 requires origin (0,0,0) to be inside the polyhedron, but you put one of the vertices there. That is not "inside", that's "on the dividing plane". That's why it fails.
In the mean time, you can use a much more robust half-space logic instead. This only requires the vertices to form a convex polyhedron, and that's it. (The origin does not need to lie within the polyhedron).
The corresponding halfspaces() function isCode:Loop over every triplet of vertex points p1, p2, p3:
Calculate n' = (p2 - p1) × (p3 - p2)
Calculate nsqr = n' · n'
If nsqr < eps*eps, skip to next triplet (as p1, p2, p3 are [almost] collinear)
Calculate n = n' / sqrt(nsqr) ( == n' / || n' || )
Calculate d = n · p1 ( == n · p2 == n · p3 )
If half-space (n, d) is already known, skip to next triplet.
Calculate the number of points inside, at, and outside the half-space.
For each vertex p, p != p1, p != p2, p != p3:
If (p · n - d > eps):
p is outside the half-space
Else, if (p · n - d < -eps):
p is inside the half-space
Else:
p is on the space-dividing plane.
If there was:
at least one point inside the half-space, and
no points outside the half-space, and
at least three points on the space-dividing plane,
then this half-space is one that defines the convex polyhedron;
add it to the list of known half-spaces.
End loop
Edit 3: Here is the final version of halfspaces() function, using the same logic as outlined earlier in this message, but only using N3/6-N2/2+N/3 iterations (instead of N3).Code:/* Given a list of 'points' vertices in 'point',
* defining a convex polyhedron,
* such that origin (0,0,0) is inside the polyhedron,
* return the number of half-spaces needed to
* define the polyhedron.
*
* Return 0 with errno set if an error occurs.
*/
size_t halfspaces(halfspace_t **const listptr,
size_t *const sizeptr,
const point_t *const point,
const size_t points,
const double epsilon)
{
const double esqr = epsilon * epsilon;
halfspace_t *list, curr;
size_t size, used;
size_t i, i1, i2, i3;
/* Sanity check. */
if (!listptr || !sizeptr || !point || points < 3) {
errno = EINVAL;
return (size_t)0;
}
/* If listptr points to NULL, then size is 0. */
list = *listptr;
if (list)
size = *sizeptr;
else
size = 0;
/* No halfspaces defined yet. */
used = 0;
for (i1 = 0; i1 < points; i1++) {
const point_t p1 = point[i1];
for (i2 = 0; i2 < points; i2++) {
if (i2 != i1) {
const point_t p2 = point[i2];
for (i3 = 0; i3 < points; i3++) {
if (i3 != i1 && i3 != i2) {
const point_t p3 = point[i3];
const point_t v2 = { p2.x - p1.x, p2.y - p1.y, p2.z - p1.z };
const point_t v3 = { p3.x - p1.x, p3.y - p1.y, p3.z - p1.z };
const point_t n = { v2.y * v3.z - v2.z * v3.y,
v2.z * v3.x - v2.x * v3.z,
v2.x * v3.y - v2.y * v3.x };
const double nn = n.x * n.x + n.y * n.y + n.z * n.z;
if (nn > esqr) {
const double nlen = sqrt(nn);
size_t n_in = 0; /* Vertices inside the half-space */
size_t n_at = 0; /* Vertices on the half-space plane */
size_t n_out = 0; /* Vertices outside the half-space */
curr.x = n.x / nlen;
curr.y = n.y / nlen;
curr.z = n.z / nlen;
curr.d = p1.x * curr.x + p1.y * curr.y + p1.z * curr.z;
/* Already included in the list? */
for (i = 0; i < used; i++)
if (fabs(curr.x - list[i].x) <= epsilon &&
fabs(curr.y - list[i].y) <= epsilon &&
fabs(curr.z - list[i].z) <= epsilon &&
fabs(curr.d - list[i].d) <= epsilon)
break;
if (i < used)
continue;
/* Count points inside and outside the half-space. */
for (i = 0; i < points; i++) {
const double d = curr.x * point[i].x
+ curr.y * point[i].y
+ curr.z * point[i].z
- curr.d;
if (d > epsilon)
n_out++;
else
if (d < -epsilon)
n_in++;
else
n_at++;
}
/* There must be:
* - at least one point inside,
* - at least three points on the plane
* - no points outside
* for this to be a valid half-space.
*/
if (n_in < 1 || n_at < 3 || n_out > 0)
continue;
/* Make room for a new halfspace. */
if (used >= size) {
size = (used | 15) + 17;
list = realloc(list, size * sizeof *list);
if (!list) {
errno = ENOMEM;
return (size_t)0;
}
*listptr = list;
*sizeptr = size;
}
list[used++] = curr;
}
}
}
}
}
}
errno = 0;
return (size_t)used;
}
The difference between this and the code in post #12 in this thread, is that this code does not rely on any extra rules for origin.Code:/* Given a list of 'points' vertices in 'point',
* defining a convex polyhedron,
* such that origin (0,0,0) is inside the polyhedron,
* return the number of half-spaces needed to
* define the polyhedron.
*
* Return 0 with errno set if an error occurs.
*/
size_t halfspaces(halfspace_t **const listptr,
size_t *const sizeptr,
const point_t *const point,
const size_t points,
const double epsilon)
{
const double esqr = epsilon * epsilon;
halfspace_t *list, curr;
size_t size, used;
size_t i, i1, i2, i3;
/* Sanity check. */
if (!listptr || !sizeptr || !point || points < 3) {
errno = EINVAL;
return (size_t)0;
}
/* If listptr points to NULL, then size is 0. */
list = *listptr;
if (list)
size = *sizeptr;
else
size = 0;
/* No halfspaces defined yet. */
used = 0;
for (i1 = 0; i1 < points - 2; i1++) {
const point_t p1 = point[i1];
for (i2 = i1 + 1; i2 < points - 1; i2++) {
const point_t p2 = point[i2];
for (i3 = i2 + 1; i3 < points; i3++) {
const point_t p3 = point[i3];
const point_t v2 = { p2.x - p1.x, p2.y - p1.y, p2.z - p1.z };
const point_t v3 = { p3.x - p1.x, p3.y - p1.y, p3.z - p1.z };
const point_t n = { v2.y * v3.z - v2.z * v3.y,
v2.z * v3.x - v2.x * v3.z,
v2.x * v3.y - v2.y * v3.x };
const double nn = n.x * n.x + n.y * n.y + n.z * n.z;
if (nn > esqr) {
const double nlen = sqrt(nn);
size_t n_in = 0; /* Vertices inside the half-space */
size_t n_at = 0; /* Vertices on the dividing plane */
size_t n_out = 0; /* Vertices outside the half-space */
curr.x = n.x / nlen;
curr.y = n.y / nlen;
curr.z = n.z / nlen;
curr.d = p1.x * curr.x + p1.y * curr.y + p1.z * curr.z;
/* Already included in the list? Inverse included? */
for (i = 0; i < used; i++)
if (fabs(curr.x - list[i].x) <= epsilon &&
fabs(curr.y - list[i].y) <= epsilon &&
fabs(curr.z - list[i].z) <= epsilon &&
fabs(curr.d - list[i].d) <= epsilon)
break;
else
if (fabs(curr.x + list[i].x) <= epsilon &&
fabs(curr.y + list[i].y) <= epsilon &&
fabs(curr.z + list[i].z) <= epsilon &&
fabs(curr.d + list[i].d) <= epsilon)
break;
if (i < used)
continue;
/* Count points inside and outside the half-space. */
for (i = 0; i < points; i++) {
const double d = curr.x * point[i].x
+ curr.y * point[i].y
+ curr.z * point[i].z
- curr.d;
if (d > epsilon)
n_out++;
else
if (d < -epsilon)
n_in++;
else
n_at++;
}
/* We need at least three points on the half-space plane,
* so that they form a true face. */
if (n_at < 3)
continue;
/* If there are points both inside and outside
* the half-space, it does not define the polyhedron. */
if (n_in > 0 && n_out > 0)
continue;
/* If we had points only outside the half-space,
* we need to reverse the half-space. Then,
* all those points were actually inside. */
if (n_out > 0 && n_in == 0) {
curr.x = -curr.x;
curr.y = -curr.y;
curr.z = -curr.z;
curr.d = -curr.d;
n_in = n_out;
n_out = 0;
}
/* If there are points outside the half-space,
* or no points inside the half-space,
* then this cannot be one that defines the polyhedron. */
if (n_in < 1 || n_out > 0)
continue;
/* Make room for a new halfspace. */
if (used >= size) {
size = (used | 15) + 17;
list = realloc(list, size * sizeof *list);
if (!list) {
errno = ENOMEM;
return (size_t)0;
}
*listptr = list;
*sizeptr = size;
}
list[used++] = curr;
}
}
}
}
errno = 0;
return (size_t)used;
}
If you fix the line defining v3 in the code in post #12, that code will work; just remember that that code requires you to specify vertices that form a polyhedron with origin inside (and not on the face, edge, or vertex) of the polyhedron.