Files
j3ml/src/J3ML/Geometry/Polygon.cpp

619 lines
22 KiB
C++

#include <J3ML/Geometry/Polygon.hpp>
#include <J3ML/Geometry/AABB.hpp>
#include <J3ML/Geometry/Triangle.hpp>
#include "J3ML/Geometry/Plane.hpp"
#include "J3ML/Geometry/Line.hpp"
#include <J3ML/Algorithm/GJK.hpp>
namespace J3ML::Geometry {
Vector3 Polygon::AnyPointFast() const { return !vertices.empty() ? vertices[0] : Vector3::NaN; }
AABB Polygon::MinimalEnclosingAABB() const {
AABB aabb;
aabb.SetNegativeInfinity();
for(int i = 0; i < NumVertices(); ++i)
aabb.Enclose(Vertex(i));
return aabb;
}
Vector3 Polygon::ExtremePoint(const Vector3 &direction, float &projectionDistance) const
{
Vector3 mostExtreme = Vector3::NaN;
projectionDistance = -INFINITY;
for(int i = 0; i < NumVertices(); ++i)
{
Vector3 pt = Vertex(i);
float d = Vector3::Dot(direction, pt);
if (d > projectionDistance)
{
projectionDistance = d;
mostExtreme = pt;
}
}
return mostExtreme;
}
Vector3 Polygon::ExtremePoint(const Vector3 &direction) const
{
float projectionDistance;
return ExtremePoint(direction, projectionDistance);
}
void Polygon::ProjectToAxis(const Vector3 &direction, float &outMin, float &outMax) const
{
///\todo Optimize!
Vector3 minPt = ExtremePoint(-direction);
Vector3 maxPt = ExtremePoint(direction);
outMin = Vector3::Dot(minPt, direction);
outMax = Vector3::Dot(maxPt, direction);
}
Vector3 Polygon::Vertex(int vertexIndex) const {
assert(vertexIndex >= 0);
assert(vertexIndex < (int) vertices.size());
return vertices[vertexIndex];
}
int Polygon::NumVertices() const {
return (int)vertices.size();
}
bool Polygon::IsPlanar(float epsilonSq) const
{
if (vertices.empty())
return false;
if (vertices.size() <= 3)
return true;
Vector3 normal = (Vector3(vertices[1])-Vector3(vertices[0])).Cross(Vector3(vertices[2])-Vector3(vertices[0]));
float lenSq = normal.LengthSquared();
for(size_t i = 3; i < vertices.size(); ++i)
{
float d = normal.Dot(Vector3(vertices[i])-Vector3(vertices[0]));
if (d*d > epsilonSq * lenSq)
return false;
}
return true;
}
Vector3 Polygon::BasisU() const
{
if (vertices.size() < 2)
return Vector3::Right;
Vector3 u = (Vector3)vertices[1] - (Vector3)vertices[0];
u = u.Normalized(); // Always succeeds, even if u was zero (generates (1,0,0)).
return u;
}
Vector3 Polygon::BasisV() const
{
if (vertices.size() < 2)
return Vector3::Down;
return Vector3::Cross(PlaneCCW().Normal, BasisU()).Normalized();
}
Plane Polygon::PlaneCCW() const
{
if (vertices.size() > 3)
{
Plane plane;
for(size_t i = 0; i < vertices.size()-2; ++i)
for(size_t j = i+1; j < vertices.size()-1; ++j)
{
Vector3 pij = Vector3(vertices[j])-Vector3(vertices[i]);
for(size_t k = j+1; k < vertices.size(); ++k)
{
plane.Normal = pij.Cross(Vector3(vertices[k])-Vector3(vertices[i]));
float lenSq = plane.Normal.LengthSquared();
if (lenSq > 1e-8f)
{
plane.Normal /= std::sqrt(lenSq);
plane.distance = plane.Normal.Dot(Vector3(vertices[i]));
return plane;
}
}
}
// warn("Polygon contains %d points, but they are all collinear! Cannot form a plane for the Polygon using three points!", (int)p.size());
// Polygon contains multiple points, but they are all collinear.
// Pick an arbitrary plane along the line as the polygon plane (as if the polygon had only two points)
Vector3 dir = (Vector3(vertices[1])-Vector3(vertices[0])).Normalized();
return Plane(Line(vertices[0], dir), dir.Perpendicular());
}
if (vertices.size() == 3)
return Plane(vertices[0], vertices[1], vertices[2]);
if (vertices.size() == 2)
{
Vector3 dir = (Vector3(vertices[1])-Vector3(vertices[0])).Normalized();
return Plane(Line(vertices[0], dir), dir.Perpendicular());
}
if (vertices.size() == 1)
return Plane(vertices[0], Vector3(0,1,0));
return Plane();
}
Vector2 Polygon::MapTo2D(int i) const
{
assert(i >= 0);
assert(i < (int)vertices.size());
return MapTo2D(vertices[i]);
}
Vector2 Polygon::MapTo2D(const Vector3 &point) const
{
assert(!vertices.empty());
Vector3 basisU = BasisU();
Vector3 basisV = BasisV();
Vector3 pt = point - vertices[0];
return Vector2(Vector3::Dot(pt, basisU), Vector3::Dot(pt, basisV));
}
Vector3 Polygon::MapFrom2D(const Vector2 &point) const
{
assert(!vertices.empty());
return (Vector3)vertices[0] + point.x * BasisU() + point.y * BasisV();
}
// A(u) = a1 + u * (a2-a1).
// B(v) = b1 + v * (b2-b1).
// Returns (u,v).
bool IntersectLineLine2D(const Vector2 &a1, const Vector2 &a2, const Vector2 &b1, const Vector2 &b2, Vector2 &out)
{
float u = (b2.x - b1.x)*(a1.y - b1.y) - (b2.y - b1.y)*(a1.x - b1.x);
float v = (a2.x - a1.x)*(a1.y - b1.y) - (a2.y - a1.y)*(a1.x - b1.x);
float det = (b2.y - b1.y)*(a2.x - a1.x) - (b2.x - b1.x)*(a2.y - a1.y);
if (std::abs(det) < 1e-4f)
return false;
det = 1.f / det;
out.x = u * det;
out.y = v * det;
return true;
}
bool IntersectLineSegmentLineSegment2D(const Vector2 &a1, const Vector2 &a2, const Vector2 &b1, const Vector2 &b2, Vector2 &out)
{
bool ret = IntersectLineLine2D(a1, a2, b1, b2, out);
return ret && out.x >= 0.f && out.x <= 1.f && out.y >= 0.f && out.y <= 1.f;
}
/// Returns true if poly[i+1] is an ear.
/// Precondition: i+2 == j (mod poly.size()).
bool IsAnEar(const std::vector<Vector2> &poly, int i, int j)
{
Vector2 dummy;
int x = (int)poly.size()-1;
for(int y = 0; y < i; ++y)
{
if (IntersectLineSegmentLineSegment2D(poly[i], poly[j], poly[x], poly[y], dummy))
return false;
x = y;
}
x = j+1;
for(int y = x+1; y < (int)poly.size(); ++y)
{
if (IntersectLineSegmentLineSegment2D(poly[i], poly[j], poly[x], poly[y], dummy))
return false;
x = y;
}
return true;
}
bool Polygon::Contains(const Polygon &worldSpacePolygon, float polygonThickness) const
{
for(int i = 0; i < worldSpacePolygon.NumVertices(); ++i)
if (!Contains(worldSpacePolygon.Vertex(i), polygonThickness))
return false;
return true;
}
bool Polygon::Contains(const Vector3 &worldSpacePoint, float polygonThicknessSq) const
{
// Implementation based on the description from http://erich.realtimerendering.com/ptinpoly/
if (vertices.size() < 3)
return false;
Vector3 basisU = BasisU();
Vector3 basisV = BasisV();
assert(basisU.IsNormalized());
assert(basisV.IsNormalized());
assert(basisU.IsPerpendicular(basisV));
assert(basisU.IsPerpendicular(PlaneCCW().Normal));
assert(basisV.IsPerpendicular(PlaneCCW().Normal));
Vector3 normal = basisU.Cross(basisV);
// float lenSq = normal.LengthSq(); ///\todo Could we treat basisU and basisV unnormalized here?
float dot = normal.Dot(Vector3(vertices[0]) - worldSpacePoint);
if (dot*dot > polygonThicknessSq)
return false; // The point is not even within the plane of the polygon - can't be contained.
int numIntersections = 0;
const float epsilon = 1e-4f;
// General strategy: transform all points on the polygon onto 2D face plane of the polygon, where the target query point is
// centered to lie in the origin.
// If the test ray (0,0) -> (+inf, 0) intersects exactly an odd number of polygon edge segments, then the query point must have been
// inside the polygon. The test ray is chosen like that to avoid all extra per-edge computations.
// This method works for both simple and non-simple (self-intersecting) polygons.
Vector3 vt = Vector3(vertices.back()) - worldSpacePoint;
Vector2 p0 = Vector2(Vector3::Dot(vt, basisU), Vector3::Dot(vt, basisV));
if (std::abs(p0.y) < epsilon)
p0.y = -epsilon; // Robustness check - if the ray (0,0) -> (+inf, 0) would pass through a vertex, move the vertex slightly.
for(int i = 0; i < (int)vertices.size(); ++i)
{
vt = Vector3(vertices[i]) - worldSpacePoint;
Vector2 p1 = Vector2(Vector3::Dot(vt, basisU), Vector3::Dot(vt, basisV));
if (std::abs(p1.y) < epsilon)
p1.y = -epsilon; // Robustness check - if the ray (0,0) -> (+inf, 0) would pass through a vertex, move the vertex slightly.
if (p0.y * p1.y < 0.f) // If the line segment p0 -> p1 straddles the line x=0, it could intersect the ray (0,0) -> (+inf, 0)
{
if (std::min(p0.x, p1.x) > 0.f) // If both x-coordinates are positive, then there certainly is an intersection with the ray.
++numIntersections;
else if (std::max(p0.x, p1.x) > 0.f) // If one of them is positive, there could be an intersection. (otherwise both are negative and they can't intersect ray)
{
// P = p0 + t*(p1-p0) == (x,0)
// p0.x + t*(p1.x-p0.x) == x
// p0.y + t*(p1.y-p0.y) == 0
// t == -p0.y / (p1.y - p0.y)
// Test whether the lines (0,0) -> (+inf,0) and p0 -> p1 intersect at a positive X-coordinate?
Vector2 d = p1 - p0;
if (d.y != 0.f)
{
float t = -p0.y / d.y; // The line segment parameter, t \in [0,1] forms the line segment p0->p1.
float x = p0.x + t * d.x; // The x-coordinate of intersection with the ray.
if (t >= 0.f && t <= 1.f && x > 0.f)
++numIntersections;
}
}
}
p0 = p1;
}
return numIntersections % 2 == 1;
}
bool Polygon::Contains(const LineSegment &worldSpaceLineSegment, float polygonThickness) const
{
if (vertices.size() < 3)
return false;
Plane plane = PlaneCCW();
if (plane.Distance(worldSpaceLineSegment.A) > polygonThickness ||
plane.Distance(worldSpaceLineSegment.B) > polygonThickness)
return false;
// For robustness, project onto the polygon plane.
LineSegment l = plane.Project(worldSpaceLineSegment);
if (!Contains(l.A) || !Contains(l.B))
return false;
for(int i = 0; i < (int)vertices.size(); ++i)
if (plane.Project(Edge(i)).Intersects(l))
return false;
return true;
}
bool Polygon::Contains(const Triangle &worldSpaceTriangle, float polygonThickness) const
{
return Contains(worldSpaceTriangle.Edge(0), polygonThickness) &&
Contains(worldSpaceTriangle.Edge(1), polygonThickness) &&
Contains(worldSpaceTriangle.Edge(2), polygonThickness);
}
/** The implementation of this function is based on the paper
"Kong, Everett, Toussant. The Graham Scan Triangulates Simple Polygons."
See also p. 772-775 of Geometric Tools for Computer Graphics.
The running time of this function is O(n^2). */
std::vector<Triangle> Polygon::Triangulate() const
{
assert(IsPlanar());
// assume1(IsPlanar(), this->SerializeToString()); // TODO: enable
std::vector<Triangle> t;
// Handle degenerate cases.
if (NumVertices() < 3)
return t;
if (NumVertices() == 3)
{
t.push_back(Triangle(Vertex(0), Vertex(1), Vertex(2)));
return t;
}
std::vector<Vector2> p2d;
std::vector<int> polyIndices;
for(int v = 0; v < NumVertices(); ++v)
{
p2d.push_back(MapTo2D(v));
polyIndices.push_back(v);
}
// Clip ears of the polygon until it has been reduced to a triangle.
int i = 0;
int j = 1;
int k = 2;
size_t numTries = 0; // Avoid creating an infinite loop.
while(p2d.size() > 3 && numTries < p2d.size())
{
if (Vector2::OrientedCCW(p2d[i], p2d[j], p2d[k]) && IsAnEar(p2d, i, k))
{
// The vertex j is an ear. Clip it off.
t.push_back(Triangle(vertices[polyIndices[i]], vertices[polyIndices[j]], vertices[polyIndices[k]]));
p2d.erase(p2d.begin() + j);
polyIndices.erase(polyIndices.begin() + j);
// The previous index might now have become an ear. Move back one index to see if so.
if (i > 0)
{
i = (i + (int)p2d.size() - 1) % p2d.size();
j = (j + (int)p2d.size() - 1) % p2d.size();
k = (k + (int)p2d.size() - 1) % p2d.size();
}
numTries = 0;
}
else
{
// The vertex at j is not an ear. Move to test next vertex.
i = j;
j = k;
k = (k+1) % p2d.size();
++numTries;
}
}
assert(p2d.size() == 3);
if (p2d.size() > 3) // If this occurs, then the polygon is NOT counter-clockwise oriented.
return t;
/*
{
// For conveniency, create a copy that has the winding order fixed, and triangulate that instead.
// (Causes a large performance hit!)
Polygon p2 = *this;
for(size_t i = 0; i < p2.p.size()/2; ++i)
std::swap(p2.p[i], p2.p[p2.p.size()-1-i]);
return p2.Triangulate();
}
*/
// Add the last poly.
t.push_back(Triangle(vertices[polyIndices[0]], vertices[polyIndices[1]], vertices[polyIndices[2]]));
return t;
}
bool Polygon::Intersects(const Capsule &capsule) const {
///@todo Optimize.
std::vector<Triangle> tris = Triangulate();
for(size_t i = 0; i < tris.size(); ++i)
if (tris[i].Intersects(capsule))
return true;
return false;
}
bool Polygon::Intersects(const Line &line) const {
float d;
if (!PlaneCCW().Intersects(line, &d))
return false;
return Contains(line.GetPoint(d));
}
bool Polygon::Intersects(const Ray &ray) const
{
float d;
if (!PlaneCCW().Intersects(ray, &d))
return false;
return Contains(ray.GetPoint(d));
}
bool Polygon::Intersects2D(const LineSegment &localSpaceLineSegment) const
{
if (vertices.size() < 3)
return false;
const Vector3 basisU = BasisU();
const Vector3 basisV = BasisV();
const Vector3 origin = vertices[0];
LineSegment edge;
edge.A = Vector3(Vector3::Dot(vertices.back(), basisU), Vector3::Dot(vertices.back(), basisV), 0); // map to 2D
for (int i = 0; i < (int)vertices.size(); ++i)
{
edge.B = Vector3(Vector3::Dot(vertices[i], basisU), Vector3::Dot(vertices[i], basisV), 0); // map to 2D
if (edge.Intersects(localSpaceLineSegment))
return true;
edge.A = edge.B;
}
// The line segment did not intersect with any of the polygon edges, so either the whole line segment is inside
// the polygon, or it is fully outside the polygon. Test one point of the line segment to determine which.
Vector2 xy = {
localSpaceLineSegment.A.x,
localSpaceLineSegment.A.y
};
return Contains(MapFrom2D(xy));
}
bool Polygon::Intersects(const LineSegment &lineSegment) const
{
Plane plane = PlaneCCW();
// Compute line-plane intersection (unroll Plane::IntersectLinePlane())
float denom = Vector3::Dot(plane.Normal, lineSegment.B - lineSegment.A);
if (std::abs(denom) < 1e-4f) // The plane of the polygon and the line are planar? Do the test in 2D.
return Intersects2D(LineSegment(Vector3(MapTo2D(lineSegment.A), 0), Vector3(MapTo2D(lineSegment.B), 0)));
// The line segment properly intersects the plane of the polygon, so there is exactly one
// point of intersection between the plane of the polygon and the line segment. Test that intersection point against
// the line segment end points.
float t = (plane.distance - Vector3::Dot(plane.Normal, lineSegment.A)) / denom;
if (t < 0.f || t > 1.f)
return false;
return Contains(lineSegment.GetPoint(t));
}
bool Polygon::Intersects(const Plane &plane) const
{
// Project the points of this polygon onto the 1D axis of the plane normal.
// If there are points on both sides of the plane, then the polygon intersects the plane.
float minD = INFINITY;
float maxD = -INFINITY;
for(size_t i = 0; i < vertices.size(); ++i)
{
float d = plane.SignedDistance(vertices[i]);
minD = std::min(minD, d);
maxD = std::max(maxD, d);
}
// Allow a very small epsilon tolerance.
return minD <= 1e-4f && maxD >= -1e-4f;
}
bool Polygon::ConvexIntersects(const AABB &aabb) const
{
return Algorithms::GJKIntersect(*this, aabb);
}
bool Polygon::ConvexIntersects(const OBB &obb) const
{
return Algorithms::GJKIntersect(*this, obb);
}
bool Polygon::ConvexIntersects(const Frustum &frustum) const
{
return Algorithms::GJKIntersect(*this, frustum);
}
template<typename Convex /* = AABB, OBB, Frustum. */>
bool Convex_Intersects_Polygon(const Convex &c, const Polygon &p)
{
LineSegment l;
l.A = p.vertices.back();
for(size_t i = 0; i < p.vertices.size(); ++i)
{
l.B = p.vertices[i];
if (c.Intersects(l))
return true;
l.A = l.B;
}
// Check all the edges of the convex shape against the polygon.
for(int i = 0; i < c.NumEdges(); ++i)
{
l = c.Edge(i);
if (p.Intersects(l))
return true;
}
return false;
}
bool Polygon::Intersects(const AABB &aabb) const
{
// Because GJK test is so fast, use that as an early-out. (computes intersection between the convex hull of this poly, and the aabb)
bool convexIntersects = ConvexIntersects(aabb);
if (!convexIntersects)
return false;
return Convex_Intersects_Polygon(aabb, *this);
}
bool Polygon::Intersects(const OBB &obb) const
{
// Because GJK test is so fast, use that as an early-out. (computes intersection between the convex hull of this poly, and the obb)
bool convexIntersects = ConvexIntersects(obb);
if (!convexIntersects)
return false;
return Convex_Intersects_Polygon(obb, *this);
}
bool Polygon::Intersects(const Frustum &frustum) const
{
// Because GJK test is so fast, use that as an early-out. (computes intersection between the convex hull of this poly, and the frustum)
bool convexIntersects = ConvexIntersects(frustum);
if (!convexIntersects)
return false;
return Convex_Intersects_Polygon(frustum, *this);
}
LineSegment Polygon::Edge(int i) const
{
if (vertices.empty())
return LineSegment(Vector3::NaN, Vector3::NaN);
if (vertices.size() == 1)
return LineSegment(vertices[0], vertices[0]);
return LineSegment(vertices[i], vertices[(i+1)%vertices.size()]);
}
Vector3 Polygon::ClosestPoint(const Vector3 &point) const
{
assert(IsPlanar());
std::vector<Triangle> tris = Triangulate();
Vector3 closestPt = Vector3::NaN;
float closestDist = FLT_MAX;
for(size_t i = 0; i < tris.size(); ++i)
{
Vector3 pt = ((Triangle)tris[i]).ClosestPoint(point);
float d = pt.DistanceSq(point);
if (d < closestDist)
{
closestPt = pt;
closestDist = d;
}
}
return closestPt;
}
Polyhedron Polygon::ToPolyhedron() const {
Polyhedron poly;
poly.v = vertices;
poly.f.push_back(Polyhedron::Face());
poly.f.push_back(Polyhedron::Face());
for(int i = 0; i < NumVertices(); ++i)
{
poly.f[0].v.push_back(i);
poly.f[1].v.push_back(NumVertices()-1-i);
}
return poly;
}
Vector3 Polygon::ClosestPoint(const LineSegment &lineSegment) const
{
return ClosestPoint(lineSegment, 0);
}
Vector3 Polygon::ClosestPoint(const LineSegment &lineSegment, Vector3 *lineSegmentPt) const
{
std::vector<Triangle> tris = Triangulate();
Vector3 closestPt = Vector3::NaN;
Vector3 closestLineSegmentPt = Vector3::NaN;
float closestDist = FLT_MAX;
for(size_t i = 0; i < tris.size(); ++i)
{
Vector3 lineSegPt;
Vector3 pt = ((Triangle)tris[i]).ClosestPoint(lineSegment, &lineSegPt);
float d = pt.DistanceSq(lineSegPt);
if (d < closestDist)
{
closestPt = pt;
closestLineSegmentPt = lineSegPt;
closestDist = d;
}
}
if (lineSegmentPt)
*lineSegmentPt = closestLineSegmentPt;
return closestPt;
}
}