BRL-CAD/GeometricTools
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Update ContPointInPolyhedron3.h

Browse Source 
Defect reported by Claude Sonnet5/Fable5. The std::vector member
mProjVertices occurs in a loop over all faces. It stores the projections
onto a plane which are later processed. The member is resized each time the current faceIndices has a larger number of elements that mProjVertices. There is no clear call, so if the number of current faceIndices is smaller than that of mProjVertices, the ensuing code accesses elements of mProjVertices that were set for the previous triangles. I removed the mProjVertices member and now have local storage for the projection vertices that is created for each incoming face. While I was in there, I refactored the code in an attempt to simply
the class interface. There is only one face data structure used for triangles, convex polygons, and simple polygons. Consequently, the class has a single constructor to which you pass a parameter specifying the type of face to be processed. I factored out shared code for ContainsC1C2 and ContainsS1 into SharedContains.
This commit is contained in:
David Eberly
2026-08-13 21:46:47 -07:00
parent 3e43b79905
commit d29e7758ae
1 changed files with 151 additions and 232 deletions
Download Patch File Download Diff File Expand all files Collapse all files
GTE/Mathematics/ContPointInPolyhedron3.h
+151 -232
View File
@@ -3,26 +3,26 @@
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt
// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// File Version: 8.0.2025.05.10
// File Version: 8.0.2026.08.13
#pragma once
// This class contains various implementations for point-in-polyhedron
// queries. The planes stored with the faces are used in all cases to
// queries. The planes stored with the faces are used in all cases to
// reject ray-face intersection tests, a quick culling operation.
//
// The algorithm is to cast a ray from the input point P and test for
// intersection against each face of the polyhedron. If the ray only
// intersection against each face of the polyhedron. If the ray only
// intersects faces at interior points (not vertices, not edge points),
// then the point is inside when the number of intersections is odd and
// the point is outside when the number of intersections is even. If the
// the point is outside when the number of intersections is even. If the
// ray intersects an edge or a vertex, then the counting must be handled
// differently. The details are tedious. As an alternative, the approach
// here is to allow you to specify 2*N+1 rays, where N >= 0. You should
// choose these rays randomly. Each ray reports "inside" or "outside".
// Whichever result occurs N+1 or more times is the "winner". The input
// rayQuantity is 2*N+1. The input array Direction must have rayQuantity
// elements. If you are feeling lucky, choose rayQuantity to be 1.
// differently. The details are tedious. As an alternative, the approach
// here is to allow you to specify 2*N+1 rays, where N >= 0. You should
// choose these rays randomly. Each ray reports "inside" or "outside".
// Whichever result occurs N+1 or more times is the "winner". The input
// rayQuantity is 2*N+1. The input array Direction must have rayQuantity
// elements. If you are feeling lucky, choose rayQuantity to be 1.
#include <Mathematics/Logger.h>
#include <Mathematics/ContPointInPolygon2.h>
@@ -39,54 +39,51 @@ namespace gte
class PointInPolyhedron3
{
public:
// For simple polyhedra with triangle faces.
class TriangleFace
enum FaceType
{
TRIANGLE,
CONVEX,
SIMPLE
};
class Face
{
public:
// The members 'indices' and 'plane' are used for triangle faces,
// for convex polygon faces, and for simple polygon faces. The
// member 'triangles' is used only for simple faces that are not
// convex.
Face()
:
indices{},
plane{},
triangles{}
{
}
// When you view the face from outside, the vertices are
// counterclockwise ordered. The indices array stores the indices
// counterclockwise ordered. The indices array stores the indices
// into the vertex array.
std::array<int32_t, 3> indices;
// The normal vector is unit length and points to the outside of
// the polyhedron.
Plane3<Real> plane;
// Each simple face may be triangulated. The indices are relative
// to the vertex array. Each triple of indices represents a
// triangle in the triangulation.
std::vector<int32_t> triangles;
};
// The Contains query will use ray-triangle intersection queries.
PointInPolyhedron3(int32_t numPoints, Vector3<Real> const* points,
int32_t numFaces, TriangleFace const* faces, int32_t numRays,
Vector3<Real> const* directions)
:
mNumPoints(numPoints),
mPoints(points),
mNumFaces(numFaces),
mTFaces(faces),
mCFaces(nullptr),
mSFaces(nullptr),
mMethod(0),
mNumRays(numRays),
mDirections(directions)
{
}
// For simple polyhedra with convex polygon faces.
class ConvexFace
{
public:
// When you view the face from outside, the vertices are
// counterclockwise ordered. The indices array stores the indices
// into the vertex array.
std::vector<int32_t> indices;
// The normal vector is unit length and points to the outside of
// the polyhedron.
Plane3<Real> plane;
};
// The Contains() query will use ray-convexpolygon intersection
// queries. A ray-convexpolygon intersection query can be implemented
// in many ways. In this context, uiMethod is one of three value:
// The Contains query uses the 'method' parameter to determine which
// point-in-polygon test is used.
//
// Triangle faces.
// 0 : The parameter is unused. Ray-triangle tests are performed
// without projection on the planes of the triangles.
//
// Convex faces.
// 0 : Use a triangle fan and perform a ray-triangle intersection
// query for each triangle.
// 1 : Find the point of intersection of ray and plane of polygon.
@@ -95,89 +92,59 @@ namespace gte
// 2 : Find the point of intersection of ray and plane of polygon.
// Test whether that point is inside the convex polygon using an
// O(log N) test.
PointInPolyhedron3(int32_t numPoints, Vector3<Real> const* points,
int32_t numFaces, ConvexFace const* faces, int32_t numRays,
Vector3<Real> const* directions, uint32_t method)
:
mNumPoints(numPoints),
mPoints(points),
mNumFaces(numFaces),
mTFaces(nullptr),
mCFaces(faces),
mSFaces(nullptr),
mMethod(method),
mNumRays(numRays),
mDirections(directions)
{
}
// For simple polyhedra with simple polygon faces that are generally
// not all convex.
class SimpleFace
{
public:
// When you view the face from outside, the vertices are
// counterclockwise ordered. The Indices array stores the indices
// into the vertex array.
std::vector<int32_t> indices;
// The normal vector is unit length and points to the outside of
// the polyhedron.
Plane3<Real> plane;
// Each simple face may be triangulated. The indices are relative
// to the vertex array. Each triple of indices represents a
// triangle in the triangulation.
std::vector<int32_t> triangles;
};
// The Contains query will use ray-simplepolygon intersection queries.
// A ray-simplepolygon intersection query can be implemented in a
// couple of ways. In this context, uiMethod is one of two value:
//
// Simple faces that are not convex.
// 0 : Iterate over the triangles of each face and perform a
// ray-triangle intersection query for each triangle. This
// requires that the SimpleFace::Triangles array be initialized
// for each face.
// ray-triangle intersection query for each triangle. This
// requires that the Face::Triangles array be initialized for
// each face.
// 1 : Find the point of intersection of ray and plane of polygon.
// Test whether that point is inside the polygon using an O(N)
// test. The SimpleFace::Triangles array is not used for this
// method, so it does not have to be initialized for each face.
PointInPolyhedron3(int32_t numPoints, Vector3<Real> const* points,
int32_t numFaces, SimpleFace const* faces, int32_t numRays,
Vector3<Real> const* directions, uint32_t method)
// test. The Face::Triangles array is not used for this method,
// so it does not have to be initialized for each face.
PointInPolyhedron3(
FaceType type,
int32_t numPoints,
Vector3<Real> const* points,
int32_t numFaces,
Face const* faces,
int32_t numRays,
Vector3<Real> const* directions,
uint32_t method)
:
mType(type),
mNumPoints(numPoints),
mPoints(points),
mNumFaces(numFaces),
mTFaces(nullptr),
mCFaces(nullptr),
mSFaces(faces),
mFaces(faces),
mMethod(method),
mNumRays(numRays),
mDirections(directions)
{
}
// This function will select the actual algorithm based on which
// constructor you used for this class.
// This function will select the actual algorithm based on the type
// of face selected in the constructor.
bool Contains(Vector3<Real> const& p) const
{
if (mTFaces)
if (mType == TRIANGLE)
{
return ContainsT0(p);
}
if (mCFaces)
if (mType == CONVEX)
{
if (mMethod == 0)
{
return ContainsC0(p);
}
return ContainsC1C2(p, mMethod);
else // mMethod is 1 or 2
{
return ContainsC1C2(p, mMethod);
}
}
if (mSFaces)
if (mType == SIMPLE)
{
if (mMethod == 0)
{
@@ -186,7 +153,7 @@ namespace gte
if (mMethod == 1)
{
return ContainsS1(p);
return ContainsS1(p, mMethod);
}
}
@@ -194,9 +161,9 @@ namespace gte
}
private:
// For all types of faces. The ray origin is the test point. The ray
// direction is one of those passed to the constructors. The plane
// origin is a point on the plane of the face. The plane normal is a
// For all types of faces. The ray origin is the test point. The ray
// direction is one of those passed to the constructors. The plane
// origin is a point on the plane of the face. The plane normal is a
// unit-length normal to the face and that points outside the
// polyhedron.
static bool FastNoIntersect(Ray3<Real> const& ray, Plane3<Real> const& plane)
@@ -227,6 +194,49 @@ namespace gte
return false;
}
using ContainsFunction = void (*)(
uint32_t,
PointInPolygon2<Real> const&,
Vector2<Real> const&,
bool&);
static void ContainsPointConvex(
uint32_t method,
PointInPolygon2<Real> const& PIP,
Vector2<Real> const& projIntersect,
bool& odd)
{
if (method == 1)
{
if (PIP.ContainsConvexOrderN(projIntersect))
{
// The ray intersects the triangle.
odd = !odd;
}
}
else
{
if (PIP.ContainsConvexOrderLogN(projIntersect))
{
// The ray intersects the triangle.
odd = !odd;
}
}
}
static void ContainsPointSimple(
uint32_t,
PointInPolygon2<Real> const& PIP,
Vector2<Real> const& projIntersect,
bool& odd)
{
if (PIP.Contains(projIntersect))
{
// The ray intersects the triangle.
odd = !odd;
}
}
// For triangle faces.
bool ContainsT0(Vector3<Real> const& p) const
{
@@ -244,7 +254,7 @@ namespace gte
// Zero intersections to start with.
bool odd = false;
TriangleFace const* face = mTFaces;
Face const* face = mFaces;
for (int32_t i = 0; i < mNumFaces; ++i, ++face)
{
// Attempt to quickly cull the triangle.
@@ -293,7 +303,7 @@ namespace gte
// Zero intersections to start with.
bool odd = false;
ConvexFace const* face = mCFaces;
Face const* face = mFaces;
for (int32_t i = 0; i < mNumFaces; ++i, ++face)
{
// Attempt to quickly cull the triangle.
@@ -327,96 +337,11 @@ namespace gte
return insideCount > mNumRays / 2;
}
bool ContainsC1C2(Vector3<Real> const& p, uint32_t method) const
bool ContainsC1C2(
Vector3<Real> const& p,
uint32_t method) const
{
int32_t insideCount = 0;
FIQuery<Real, Ray3<Real>, Plane3<Real>> rpQuery;
Ray3<Real> ray;
ray.origin = p;
for (int32_t j = 0; j < mNumRays; ++j)
{
ray.direction = mDirections[j];
// Zero intersections to start with.
bool odd = false;
ConvexFace const* face = mCFaces;
for (int32_t i = 0; i < mNumFaces; ++i, ++face)
{
// Attempt to quickly cull the triangle.
if (FastNoIntersect(ray, face->plane))
{
continue;
}
// Compute the ray-plane intersection.
auto result = rpQuery(ray, face->plane);
// If you trigger this assertion, numerical round-off
// errors have led to a discrepancy between
// FastNoIntersect and the Find() result.
LogAssert(result.intersect, "Unexpected condition.");
// Get a coordinate system for the plane. Use vertex 0
// as the origin.
Vector3<Real> const& V0 = mPoints[face->indices[0]];
Vector3<Real> basis[3];
basis[0] = face->plane.normal;
ComputeOrthogonalComplement(1, basis);
// Project the intersection onto the plane.
Vector3<Real> diff = result.point - V0;
Vector2<Real> projIntersect{ Dot(basis[1], diff), Dot(basis[2], diff) };
// Project the face vertices onto the plane of the face.
if (face->indices.size() > mProjVertices.size())
{
mProjVertices.resize(face->indices.size());
}
// Project the remaining vertices. Vertex 0 is always the
// origin.
size_t numIndices = face->indices.size();
mProjVertices[0] = Vector2<Real>::Zero();
for (size_t k = 1; k < numIndices; ++k)
{
diff = mPoints[face->indices[k]] - V0;
mProjVertices[k][0] = Dot(basis[1], diff);
mProjVertices[k][1] = Dot(basis[2], diff);
}
// Test whether the intersection point is in the convex
// polygon.
PointInPolygon2<Real> PIP(static_cast<int32_t>(mProjVertices.size()),
&mProjVertices[0]);
if (method == 1)
{
if (PIP.ContainsConvexOrderN(projIntersect))
{
// The ray intersects the triangle.
odd = !odd;
}
}
else
{
if (PIP.ContainsConvexOrderLogN(projIntersect))
{
// The ray intersects the triangle.
odd = !odd;
}
}
}
if (odd)
{
insideCount++;
}
}
return insideCount > mNumRays / 2;
return SharedContains(p, method, ContainsPointConvex);
}
// For simple faces.
@@ -436,7 +361,7 @@ namespace gte
// Zero intersections to start with.
bool odd = false;
SimpleFace const* face = mSFaces;
Face const* face = mFaces;
for (int32_t i = 0; i < mNumFaces; ++i, ++face)
{
// Attempt to quickly cull the triangle.
@@ -477,7 +402,18 @@ namespace gte
return insideCount > mNumRays / 2;
}
bool ContainsS1(Vector3<Real> const& p) const
bool ContainsS1(
Vector3<Real> const& p,
uint32_t method) const
{
return SharedContains(p, method, ContainsPointSimple);
}
// Shared code for ContainsC1C2 and ContainsS1.
bool SharedContains(
Vector3<Real> const& p,
uint32_t method,
ContainsFunction Contains) const
{
int32_t insideCount = 0;
@@ -492,7 +428,7 @@ namespace gte
// Zero intersections to start with.
bool odd = false;
SimpleFace const* face = mSFaces;
Face const* face = mFaces;
for (int32_t i = 0; i < mNumFaces; ++i, ++face)
{
// Attempt to quickly cull the triangle.
@@ -509,44 +445,35 @@ namespace gte
// FastNoIntersect and the Find() result.
LogAssert(result.intersect, "Unexpected condition.");
// Get a coordinate system for the plane. Use vertex 0
// Get a coordinate system for the plane. Use vertex 0
// as the origin.
Vector3<Real> const& V0 = mPoints[face->indices[0]];
Vector3<Real> basis[3];
std::array<Vector3<Real>, 3> basis{};
basis[0] = face->plane.normal;
ComputeOrthogonalComplement(1, basis);
ComputeOrthogonalComplement(1, basis.data());
// Project the intersection onto the plane.
Vector3<Real> diff = result.point - V0;
Vector2<Real> projIntersect{ Dot(basis[1], diff), Dot(basis[2], diff) };
// Project the face vertices onto the plane of the face.
if (face->indices.size() > mProjVertices.size())
{
mProjVertices.resize(face->indices.size());
}
// Project the remaining vertices. Vertex 0 is always the
// origin.
// Vertex 0 is always the origin.
size_t numIndices = face->indices.size();
mProjVertices[0] = Vector2<Real>::Zero();
std::vector<Vector2<Real>> projVertices(numIndices);
projVertices[0] = Vector2<Real>::Zero();
for (size_t k = 1; k < numIndices; ++k)
{
diff = mPoints[face->indices[k]] - V0;
mProjVertices[k][0] = Dot(basis[1], diff);
mProjVertices[k][1] = Dot(basis[2], diff);
projVertices[k][0] = Dot(basis[1], diff);
projVertices[k][1] = Dot(basis[2], diff);
}
// Test whether the intersection point is in the convex
// polygon.
PointInPolygon2<Real> PIP(static_cast<int32_t>(mProjVertices.size()),
&mProjVertices[0]);
PointInPolygon2<Real> PIP(static_cast<int32_t>(projVertices.size()),
projVertices.data());
if (PIP.Contains(projIntersect))
{
// The ray intersects the triangle.
odd = !odd;
}
Contains(method, PIP, projIntersect, odd);
}
if (odd)
@@ -558,22 +485,14 @@ namespace gte
return insideCount > mNumRays / 2;
}
FaceType mType;
int32_t mNumPoints;
Vector3<Real> const* mPoints;
int32_t mNumFaces;
TriangleFace const* mTFaces;
ConvexFace const* mCFaces;
SimpleFace const* mSFaces;
Face const* mFaces;
uint32_t mMethod;
int32_t mNumRays;
Vector3<Real> const* mDirections;
// Temporary storage for those methods that reduce the problem to 2D
// point-in-polygon queries. The array stores the projections of
// face vertices onto the plane of the face. It is resized as needed.
mutable std::vector<Vector2<Real>> mProjVertices;
};
}