Path: blob/master/thirdparty/jolt_physics/Jolt/Physics/Collision/Shape/ConvexHullShape.cpp
9913 views
// Jolt Physics Library (https://github.com/jrouwe/JoltPhysics)1// SPDX-FileCopyrightText: 2021 Jorrit Rouwe2// SPDX-License-Identifier: MIT34#include <Jolt/Jolt.h>56#include <Jolt/Physics/Collision/Shape/ConvexHullShape.h>7#include <Jolt/Physics/Collision/Shape/ScaleHelpers.h>8#include <Jolt/Physics/Collision/Shape/PolyhedronSubmergedVolumeCalculator.h>9#include <Jolt/Physics/Collision/RayCast.h>10#include <Jolt/Physics/Collision/CastResult.h>11#include <Jolt/Physics/Collision/CollidePointResult.h>12#include <Jolt/Physics/Collision/TransformedShape.h>13#include <Jolt/Physics/Collision/CollideSoftBodyVertexIterator.h>14#include <Jolt/Geometry/ConvexHullBuilder.h>15#include <Jolt/Geometry/ClosestPoint.h>16#include <Jolt/ObjectStream/TypeDeclarations.h>17#include <Jolt/Core/StringTools.h>18#include <Jolt/Core/StreamIn.h>19#include <Jolt/Core/StreamOut.h>20#include <Jolt/Core/UnorderedMap.h>2122JPH_NAMESPACE_BEGIN2324JPH_IMPLEMENT_SERIALIZABLE_VIRTUAL(ConvexHullShapeSettings)25{26JPH_ADD_BASE_CLASS(ConvexHullShapeSettings, ConvexShapeSettings)2728JPH_ADD_ATTRIBUTE(ConvexHullShapeSettings, mPoints)29JPH_ADD_ATTRIBUTE(ConvexHullShapeSettings, mMaxConvexRadius)30JPH_ADD_ATTRIBUTE(ConvexHullShapeSettings, mMaxErrorConvexRadius)31JPH_ADD_ATTRIBUTE(ConvexHullShapeSettings, mHullTolerance)32}3334ShapeSettings::ShapeResult ConvexHullShapeSettings::Create() const35{36if (mCachedResult.IsEmpty())37Ref<Shape> shape = new ConvexHullShape(*this, mCachedResult);38return mCachedResult;39}4041ConvexHullShape::ConvexHullShape(const ConvexHullShapeSettings &inSettings, ShapeResult &outResult) :42ConvexShape(EShapeSubType::ConvexHull, inSettings, outResult),43mConvexRadius(inSettings.mMaxConvexRadius)44{45using BuilderFace = ConvexHullBuilder::Face;46using Edge = ConvexHullBuilder::Edge;47using Faces = Array<BuilderFace *>;4849// Check convex radius50if (mConvexRadius < 0.0f)51{52outResult.SetError("Invalid convex radius");53return;54}5556// Build convex hull57const char *error = nullptr;58ConvexHullBuilder builder(inSettings.mPoints);59ConvexHullBuilder::EResult result = builder.Initialize(cMaxPointsInHull, inSettings.mHullTolerance, error);60if (result != ConvexHullBuilder::EResult::Success && result != ConvexHullBuilder::EResult::MaxVerticesReached)61{62outResult.SetError(error);63return;64}65const Faces &builder_faces = builder.GetFaces();6667// Check the consistency of the resulting hull if we fully built it68if (result == ConvexHullBuilder::EResult::Success)69{70ConvexHullBuilder::Face *max_error_face;71float max_error_distance, coplanar_distance;72int max_error_idx;73builder.DetermineMaxError(max_error_face, max_error_distance, max_error_idx, coplanar_distance);74if (max_error_distance > 4.0f * max(coplanar_distance, inSettings.mHullTolerance)) // Coplanar distance could be bigger than the allowed tolerance if the points are far apart75{76outResult.SetError(StringFormat("Hull building failed, point %d had an error of %g (relative to tolerance: %g)", max_error_idx, (double)max_error_distance, double(max_error_distance / inSettings.mHullTolerance)));77return;78}79}8081// Calculate center of mass and volume82builder.GetCenterOfMassAndVolume(mCenterOfMass, mVolume);8384// Calculate covariance matrix85// See:86// - Why the inertia tensor is the inertia tensor - Jonathan Blow (http://number-none.com/blow/inertia/deriving_i.html)87// - How to find the inertia tensor (or other mass properties) of a 3D solid body represented by a triangle mesh (Draft) - Jonathan Blow, Atman J Binstock (http://number-none.com/blow/inertia/bb_inertia.doc)88Mat44 covariance_canonical(Vec4(1.0f / 60.0f, 1.0f / 120.0f, 1.0f / 120.0f, 0), Vec4(1.0f / 120.0f, 1.0f / 60.0f, 1.0f / 120.0f, 0), Vec4(1.0f / 120.0f, 1.0f / 120.0f, 1.0f / 60.0f, 0), Vec4(0, 0, 0, 1));89Mat44 covariance_matrix = Mat44::sZero();90for (BuilderFace *f : builder_faces)91{92// Fourth point of the tetrahedron is at the center of mass, we subtract it from the other points so we get a tetrahedron with one vertex at zero93// The first point on the face will be used to form a triangle fan94Edge *e = f->mFirstEdge;95Vec3 v1 = inSettings.mPoints[e->mStartIdx] - mCenterOfMass;9697// Get the 2nd point98e = e->mNextEdge;99Vec3 v2 = inSettings.mPoints[e->mStartIdx] - mCenterOfMass;100101// Loop over the triangle fan102for (e = e->mNextEdge; e != f->mFirstEdge; e = e->mNextEdge)103{104Vec3 v3 = inSettings.mPoints[e->mStartIdx] - mCenterOfMass;105106// Affine transform that transforms a unit tetrahedron (with vertices (0, 0, 0), (1, 0, 0), (0, 1, 0) and (0, 0, 1) to this tetrahedron107Mat44 a(Vec4(v1, 0), Vec4(v2, 0), Vec4(v3, 0), Vec4(0, 0, 0, 1));108109// Calculate covariance matrix for this tetrahedron110float det_a = a.GetDeterminant3x3();111Mat44 c = det_a * (a * covariance_canonical * a.Transposed());112113// Add it114covariance_matrix += c;115116// Prepare for next triangle117v2 = v3;118}119}120121// Calculate inertia matrix assuming density is 1, note that element (3, 3) is garbage122mInertia = Mat44::sIdentity() * (covariance_matrix(0, 0) + covariance_matrix(1, 1) + covariance_matrix(2, 2)) - covariance_matrix;123124// Convert polygons from the builder to our internal representation125using VtxMap = UnorderedMap<int, uint8>;126VtxMap vertex_map;127vertex_map.reserve(VtxMap::size_type(inSettings.mPoints.size()));128for (BuilderFace *builder_face : builder_faces)129{130// Determine where the vertices go131JPH_ASSERT(mVertexIdx.size() <= 0xFFFF);132uint16 first_vertex = (uint16)mVertexIdx.size();133uint16 num_vertices = 0;134135// Loop over vertices in face136Edge *edge = builder_face->mFirstEdge;137do138{139// Remap to new index, not all points in the original input set are required to form the hull140uint8 new_idx;141int original_idx = edge->mStartIdx;142VtxMap::iterator m = vertex_map.find(original_idx);143if (m != vertex_map.end())144{145// Found, reuse146new_idx = m->second;147}148else149{150// This is a new point151// Make relative to center of mass152Vec3 p = inSettings.mPoints[original_idx] - mCenterOfMass;153154// Update local bounds155mLocalBounds.Encapsulate(p);156157// Add to point list158JPH_ASSERT(mPoints.size() <= 0xff);159new_idx = (uint8)mPoints.size();160mPoints.push_back({ p });161vertex_map[original_idx] = new_idx;162}163164// Append to vertex list165JPH_ASSERT(mVertexIdx.size() < 0xffff);166mVertexIdx.push_back(new_idx);167num_vertices++;168169edge = edge->mNextEdge;170} while (edge != builder_face->mFirstEdge);171172// Add face173mFaces.push_back({ first_vertex, num_vertices });174175// Add plane176Plane plane = Plane::sFromPointAndNormal(builder_face->mCentroid - mCenterOfMass, builder_face->mNormal.Normalized());177mPlanes.push_back(plane);178}179180// Test if GetSupportFunction can support this many points181if (mPoints.size() > cMaxPointsInHull)182{183outResult.SetError(StringFormat("Internal error: Too many points in hull (%u), max allowed %d", (uint)mPoints.size(), cMaxPointsInHull));184return;185}186187for (int p = 0; p < (int)mPoints.size(); ++p)188{189// For each point, find faces that use the point190Array<int> faces;191for (int f = 0; f < (int)mFaces.size(); ++f)192{193const Face &face = mFaces[f];194for (int v = 0; v < face.mNumVertices; ++v)195if (mVertexIdx[face.mFirstVertex + v] == p)196{197faces.push_back(f);198break;199}200}201202if (faces.size() < 2)203{204outResult.SetError("A point must be connected to 2 or more faces!");205return;206}207208// Find the 3 normals that form the largest tetrahedron209// The largest tetrahedron we can get is ((1, 0, 0) x (0, 1, 0)) . (0, 0, 1) = 1, if the volume is only 5% of that,210// the three vectors are too coplanar and we fall back to using only 2 plane normals211float biggest_volume = 0.05f;212int best3[3] = { -1, -1, -1 };213214// When using 2 normals, we get the two with the biggest angle between them with a minimal difference of 1 degree215// otherwise we fall back to just using 1 plane normal216float smallest_dot = Cos(DegreesToRadians(1.0f));217int best2[2] = { -1, -1 };218219for (int face1 = 0; face1 < (int)faces.size(); ++face1)220{221Vec3 normal1 = mPlanes[faces[face1]].GetNormal();222for (int face2 = face1 + 1; face2 < (int)faces.size(); ++face2)223{224Vec3 normal2 = mPlanes[faces[face2]].GetNormal();225Vec3 cross = normal1.Cross(normal2);226227// Determine the 2 face normals that are most apart228float dot = normal1.Dot(normal2);229if (dot < smallest_dot)230{231smallest_dot = dot;232best2[0] = faces[face1];233best2[1] = faces[face2];234}235236// Determine the 3 face normals that form the largest tetrahedron237for (int face3 = face2 + 1; face3 < (int)faces.size(); ++face3)238{239Vec3 normal3 = mPlanes[faces[face3]].GetNormal();240float volume = abs(cross.Dot(normal3));241if (volume > biggest_volume)242{243biggest_volume = volume;244best3[0] = faces[face1];245best3[1] = faces[face2];246best3[2] = faces[face3];247}248}249}250}251252// If we didn't find 3 planes, use 2, if we didn't find 2 use 1253if (best3[0] != -1)254faces = { best3[0], best3[1], best3[2] };255else if (best2[0] != -1)256faces = { best2[0], best2[1] };257else258faces = { faces[0] };259260// Copy the faces to the points buffer261Point &point = mPoints[p];262point.mNumFaces = (int)faces.size();263for (int i = 0; i < (int)faces.size(); ++i)264point.mFaces[i] = faces[i];265}266267// If the convex radius is already zero, there's no point in further reducing it268if (mConvexRadius > 0.0f)269{270// Find out how thin the hull is by walking over all planes and checking the thickness of the hull in that direction271float min_size = FLT_MAX;272for (const Plane &plane : mPlanes)273{274// Take the point that is furthest away from the plane as thickness of this hull275float max_dist = 0.0f;276for (const Point &point : mPoints)277{278float dist = -plane.SignedDistance(point.mPosition); // Point is always behind plane, so we need to negate279if (dist > max_dist)280max_dist = dist;281}282min_size = min(min_size, max_dist);283}284285// We need to fit in 2x the convex radius in min_size, so reduce the convex radius if it's bigger than that286mConvexRadius = min(mConvexRadius, 0.5f * min_size);287}288289// Now walk over all points and see if we have to further reduce the convex radius because of sharp edges290if (mConvexRadius > 0.0f)291{292for (const Point &point : mPoints)293if (point.mNumFaces != 1) // If we have a single face, shifting back is easy and we don't need to reduce the convex radius294{295// Get first two planes296Plane p1 = mPlanes[point.mFaces[0]];297Plane p2 = mPlanes[point.mFaces[1]];298Plane p3;299Vec3 offset_mask;300301if (point.mNumFaces == 3)302{303// Get third plane304p3 = mPlanes[point.mFaces[2]];305306// All 3 planes will be offset by the convex radius307offset_mask = Vec3::sReplicate(1);308}309else310{311// Third plane has normal perpendicular to the other two planes and goes through the vertex position312JPH_ASSERT(point.mNumFaces == 2);313p3 = Plane::sFromPointAndNormal(point.mPosition, p1.GetNormal().Cross(p2.GetNormal()));314315// Only the first and 2nd plane will be offset, the 3rd plane is only there to guide the intersection point316offset_mask = Vec3(1, 1, 0);317}318319// Plane equation: point . normal + constant = 0320// Offsetting the plane backwards with convex radius r: point . normal + constant + r = 0321// To find the intersection 'point' of 3 planes we solve:322// |n1x n1y n1z| |x| | r + c1 |323// |n2x n2y n2z| |y| = - | r + c2 | <=> n point = -r (1, 1, 1) - (c1, c2, c3)324// |n3x n3y n3z| |z| | r + c3 |325// Where point = (x, y, z), n1x is the x component of the first plane, c1 = plane constant of plane 1, etc.326// The relation between how much the intersection point shifts as a function of r is: -r * n^-1 (1, 1, 1) = r * offset327// Where offset = -n^-1 (1, 1, 1) or -n^-1 (1, 1, 0) in case only the first 2 planes are offset328// The error that is introduced by a convex radius r is: error = r * |offset| - r329// So the max convex radius given error is: r = error / (|offset| - 1)330Mat44 n = Mat44(Vec4(p1.GetNormal(), 0), Vec4(p2.GetNormal(), 0), Vec4(p3.GetNormal(), 0), Vec4(0, 0, 0, 1)).Transposed();331float det_n = n.GetDeterminant3x3();332if (det_n == 0.0f)333{334// If the determinant is zero, the matrix is not invertible so no solution exists to move the point backwards and we have to choose a convex radius of zero335mConvexRadius = 0.0f;336break;337}338Mat44 adj_n = n.Adjointed3x3();339float offset = ((adj_n * offset_mask) / det_n).Length();340JPH_ASSERT(offset > 1.0f);341float max_convex_radius = inSettings.mMaxErrorConvexRadius / (offset - 1.0f);342mConvexRadius = min(mConvexRadius, max_convex_radius);343}344}345346// Calculate the inner radius by getting the minimum distance from the origin to the planes of the hull347mInnerRadius = FLT_MAX;348for (const Plane &p : mPlanes)349mInnerRadius = min(mInnerRadius, -p.GetConstant());350mInnerRadius = max(0.0f, mInnerRadius); // Clamp against zero, this should do nothing as the shape is centered around the center of mass but for flat convex hulls there may be numerical round off issues351352outResult.Set(this);353}354355MassProperties ConvexHullShape::GetMassProperties() const356{357MassProperties p;358359float density = GetDensity();360361// Calculate mass362p.mMass = density * mVolume;363364// Calculate inertia matrix365p.mInertia = density * mInertia;366p.mInertia(3, 3) = 1.0f;367368return p;369}370371Vec3 ConvexHullShape::GetSurfaceNormal(const SubShapeID &inSubShapeID, Vec3Arg inLocalSurfacePosition) const372{373JPH_ASSERT(inSubShapeID.IsEmpty(), "Invalid subshape ID");374375const Plane &first_plane = mPlanes[0];376Vec3 best_normal = first_plane.GetNormal();377float best_dist = abs(first_plane.SignedDistance(inLocalSurfacePosition));378379// Find the face that has the shortest distance to the surface point380for (Array<Face>::size_type i = 1; i < mFaces.size(); ++i)381{382const Plane &plane = mPlanes[i];383Vec3 plane_normal = plane.GetNormal();384float dist = abs(plane.SignedDistance(inLocalSurfacePosition));385if (dist < best_dist)386{387best_dist = dist;388best_normal = plane_normal;389}390}391392return best_normal;393}394395class ConvexHullShape::HullNoConvex final : public Support396{397public:398explicit HullNoConvex(float inConvexRadius) :399mConvexRadius(inConvexRadius)400{401static_assert(sizeof(HullNoConvex) <= sizeof(SupportBuffer), "Buffer size too small");402JPH_ASSERT(IsAligned(this, alignof(HullNoConvex)));403}404405virtual Vec3 GetSupport(Vec3Arg inDirection) const override406{407// Find the point with the highest projection on inDirection408float best_dot = -FLT_MAX;409Vec3 best_point = Vec3::sZero();410411for (Vec3 point : mPoints)412{413// Check if its support is bigger than the current max414float dot = point.Dot(inDirection);415if (dot > best_dot)416{417best_dot = dot;418best_point = point;419}420}421422return best_point;423}424425virtual float GetConvexRadius() const override426{427return mConvexRadius;428}429430using PointsArray = StaticArray<Vec3, cMaxPointsInHull>;431432inline PointsArray & GetPoints()433{434return mPoints;435}436437const PointsArray & GetPoints() const438{439return mPoints;440}441442private:443float mConvexRadius;444PointsArray mPoints;445};446447class ConvexHullShape::HullWithConvex final : public Support448{449public:450explicit HullWithConvex(const ConvexHullShape *inShape) :451mShape(inShape)452{453static_assert(sizeof(HullWithConvex) <= sizeof(SupportBuffer), "Buffer size too small");454JPH_ASSERT(IsAligned(this, alignof(HullWithConvex)));455}456457virtual Vec3 GetSupport(Vec3Arg inDirection) const override458{459// Find the point with the highest projection on inDirection460float best_dot = -FLT_MAX;461Vec3 best_point = Vec3::sZero();462463for (const Point &point : mShape->mPoints)464{465// Check if its support is bigger than the current max466float dot = point.mPosition.Dot(inDirection);467if (dot > best_dot)468{469best_dot = dot;470best_point = point.mPosition;471}472}473474return best_point;475}476477virtual float GetConvexRadius() const override478{479return 0.0f;480}481482private:483const ConvexHullShape * mShape;484};485486class ConvexHullShape::HullWithConvexScaled final : public Support487{488public:489HullWithConvexScaled(const ConvexHullShape *inShape, Vec3Arg inScale) :490mShape(inShape),491mScale(inScale)492{493static_assert(sizeof(HullWithConvexScaled) <= sizeof(SupportBuffer), "Buffer size too small");494JPH_ASSERT(IsAligned(this, alignof(HullWithConvexScaled)));495}496497virtual Vec3 GetSupport(Vec3Arg inDirection) const override498{499// Find the point with the highest projection on inDirection500float best_dot = -FLT_MAX;501Vec3 best_point = Vec3::sZero();502503for (const Point &point : mShape->mPoints)504{505// Calculate scaled position506Vec3 pos = mScale * point.mPosition;507508// Check if its support is bigger than the current max509float dot = pos.Dot(inDirection);510if (dot > best_dot)511{512best_dot = dot;513best_point = pos;514}515}516517return best_point;518}519520virtual float GetConvexRadius() const override521{522return 0.0f;523}524525private:526const ConvexHullShape * mShape;527Vec3 mScale;528};529530const ConvexShape::Support *ConvexHullShape::GetSupportFunction(ESupportMode inMode, SupportBuffer &inBuffer, Vec3Arg inScale) const531{532// If there's no convex radius, we don't need to shrink the hull533if (mConvexRadius == 0.0f)534{535if (ScaleHelpers::IsNotScaled(inScale))536return new (&inBuffer) HullWithConvex(this);537else538return new (&inBuffer) HullWithConvexScaled(this, inScale);539}540541switch (inMode)542{543case ESupportMode::IncludeConvexRadius:544case ESupportMode::Default:545if (ScaleHelpers::IsNotScaled(inScale))546return new (&inBuffer) HullWithConvex(this);547else548return new (&inBuffer) HullWithConvexScaled(this, inScale);549550case ESupportMode::ExcludeConvexRadius:551if (ScaleHelpers::IsNotScaled(inScale))552{553// Create support function554HullNoConvex *hull = new (&inBuffer) HullNoConvex(mConvexRadius);555HullNoConvex::PointsArray &transformed_points = hull->GetPoints();556JPH_ASSERT(mPoints.size() <= cMaxPointsInHull, "Not enough space, this should have been caught during shape creation!");557558for (const Point &point : mPoints)559{560Vec3 new_point;561562if (point.mNumFaces == 1)563{564// Simply shift back by the convex radius using our 1 plane565new_point = point.mPosition - mPlanes[point.mFaces[0]].GetNormal() * mConvexRadius;566}567else568{569// Get first two planes and offset inwards by convex radius570Plane p1 = mPlanes[point.mFaces[0]].Offset(-mConvexRadius);571Plane p2 = mPlanes[point.mFaces[1]].Offset(-mConvexRadius);572Plane p3;573574if (point.mNumFaces == 3)575{576// Get third plane and offset inwards by convex radius577p3 = mPlanes[point.mFaces[2]].Offset(-mConvexRadius);578}579else580{581// Third plane has normal perpendicular to the other two planes and goes through the vertex position582JPH_ASSERT(point.mNumFaces == 2);583p3 = Plane::sFromPointAndNormal(point.mPosition, p1.GetNormal().Cross(p2.GetNormal()));584}585586// Find intersection point between the three planes587if (!Plane::sIntersectPlanes(p1, p2, p3, new_point))588{589// Fallback: Just push point back using the first plane590new_point = point.mPosition - p1.GetNormal() * mConvexRadius;591}592}593594// Add point595transformed_points.push_back(new_point);596}597598return hull;599}600else601{602// Calculate scaled convex radius603float convex_radius = ScaleHelpers::ScaleConvexRadius(mConvexRadius, inScale);604605// Create new support function606HullNoConvex *hull = new (&inBuffer) HullNoConvex(convex_radius);607HullNoConvex::PointsArray &transformed_points = hull->GetPoints();608JPH_ASSERT(mPoints.size() <= cMaxPointsInHull, "Not enough space, this should have been caught during shape creation!");609610// Precalculate inverse scale611Vec3 inv_scale = inScale.Reciprocal();612613for (const Point &point : mPoints)614{615// Calculate scaled position616Vec3 pos = inScale * point.mPosition;617618// Transform normals for plane 1 with scale619Vec3 n1 = (inv_scale * mPlanes[point.mFaces[0]].GetNormal()).Normalized();620621Vec3 new_point;622623if (point.mNumFaces == 1)624{625// Simply shift back by the convex radius using our 1 plane626new_point = pos - n1 * convex_radius;627}628else629{630// Transform normals for plane 2 with scale631Vec3 n2 = (inv_scale * mPlanes[point.mFaces[1]].GetNormal()).Normalized();632633// Get first two planes and offset inwards by convex radius634Plane p1 = Plane::sFromPointAndNormal(pos, n1).Offset(-convex_radius);635Plane p2 = Plane::sFromPointAndNormal(pos, n2).Offset(-convex_radius);636Plane p3;637638if (point.mNumFaces == 3)639{640// Transform last normal with scale641Vec3 n3 = (inv_scale * mPlanes[point.mFaces[2]].GetNormal()).Normalized();642643// Get third plane and offset inwards by convex radius644p3 = Plane::sFromPointAndNormal(pos, n3).Offset(-convex_radius);645}646else647{648// Third plane has normal perpendicular to the other two planes and goes through the vertex position649JPH_ASSERT(point.mNumFaces == 2);650p3 = Plane::sFromPointAndNormal(pos, n1.Cross(n2));651}652653// Find intersection point between the three planes654if (!Plane::sIntersectPlanes(p1, p2, p3, new_point))655{656// Fallback: Just push point back using the first plane657new_point = pos - n1 * convex_radius;658}659}660661// Add point662transformed_points.push_back(new_point);663}664665return hull;666}667}668669JPH_ASSERT(false);670return nullptr;671}672673void ConvexHullShape::GetSupportingFace(const SubShapeID &inSubShapeID, Vec3Arg inDirection, Vec3Arg inScale, Mat44Arg inCenterOfMassTransform, SupportingFace &outVertices) const674{675JPH_ASSERT(inSubShapeID.IsEmpty(), "Invalid subshape ID");676677Vec3 inv_scale = inScale.Reciprocal();678679// Need to transform the plane normals using inScale680// Transforming a direction with matrix M is done through multiplying by (M^-1)^T681// In this case M is a diagonal matrix with the scale vector, so we need to multiply our normal by 1 / scale and renormalize afterwards682Vec3 plane0_normal = inv_scale * mPlanes[0].GetNormal();683float best_dot = plane0_normal.Dot(inDirection) / plane0_normal.Length();684int best_face_idx = 0;685686for (Array<Plane>::size_type i = 1; i < mPlanes.size(); ++i)687{688Vec3 plane_normal = inv_scale * mPlanes[i].GetNormal();689float dot = plane_normal.Dot(inDirection) / plane_normal.Length();690if (dot < best_dot)691{692best_dot = dot;693best_face_idx = (int)i;694}695}696697// Get vertices698const Face &best_face = mFaces[best_face_idx];699const uint8 *first_vtx = mVertexIdx.data() + best_face.mFirstVertex;700const uint8 *end_vtx = first_vtx + best_face.mNumVertices;701702// If we have more than 1/2 the capacity of outVertices worth of vertices, we start skipping vertices (note we can't fill the buffer completely since extra edges will be generated by clipping).703// TODO: This really needs a better algorithm to determine which vertices are important!704int max_vertices_to_return = outVertices.capacity() / 2;705int delta_vtx = (int(best_face.mNumVertices) + max_vertices_to_return) / max_vertices_to_return;706707// Calculate transform with scale708Mat44 transform = inCenterOfMassTransform.PreScaled(inScale);709710if (ScaleHelpers::IsInsideOut(inScale))711{712// Flip winding of supporting face713for (const uint8 *v = end_vtx - 1; v >= first_vtx; v -= delta_vtx)714outVertices.push_back(transform * mPoints[*v].mPosition);715}716else717{718// Normal winding of supporting face719for (const uint8 *v = first_vtx; v < end_vtx; v += delta_vtx)720outVertices.push_back(transform * mPoints[*v].mPosition);721}722}723724void ConvexHullShape::GetSubmergedVolume(Mat44Arg inCenterOfMassTransform, Vec3Arg inScale, const Plane &inSurface, float &outTotalVolume, float &outSubmergedVolume, Vec3 &outCenterOfBuoyancy JPH_IF_DEBUG_RENDERER(, RVec3Arg inBaseOffset)) const725{726// Trivially calculate total volume727Vec3 abs_scale = inScale.Abs();728outTotalVolume = mVolume * abs_scale.GetX() * abs_scale.GetY() * abs_scale.GetZ();729730// Check if shape has been scaled inside out731bool is_inside_out = ScaleHelpers::IsInsideOut(inScale);732733// Convert the points to world space and determine the distance to the surface734int num_points = int(mPoints.size());735PolyhedronSubmergedVolumeCalculator::Point *buffer = (PolyhedronSubmergedVolumeCalculator::Point *)JPH_STACK_ALLOC(num_points * sizeof(PolyhedronSubmergedVolumeCalculator::Point));736PolyhedronSubmergedVolumeCalculator submerged_vol_calc(inCenterOfMassTransform * Mat44::sScale(inScale), &mPoints[0].mPosition, sizeof(Point), num_points, inSurface, buffer JPH_IF_DEBUG_RENDERER(, inBaseOffset));737738if (submerged_vol_calc.AreAllAbove())739{740// We're above the water741outSubmergedVolume = 0.0f;742outCenterOfBuoyancy = Vec3::sZero();743}744else if (submerged_vol_calc.AreAllBelow())745{746// We're fully submerged747outSubmergedVolume = outTotalVolume;748outCenterOfBuoyancy = inCenterOfMassTransform.GetTranslation();749}750else751{752// Calculate submerged volume753int reference_point_idx = submerged_vol_calc.GetReferencePointIdx();754for (const Face &f : mFaces)755{756const uint8 *first_vtx = mVertexIdx.data() + f.mFirstVertex;757const uint8 *end_vtx = first_vtx + f.mNumVertices;758759// If any of the vertices of this face are the reference point, the volume will be zero so we can skip this face760bool degenerate = false;761for (const uint8 *v = first_vtx; v < end_vtx; ++v)762if (*v == reference_point_idx)763{764degenerate = true;765break;766}767if (degenerate)768continue;769770// Triangulate the face771int i1 = *first_vtx;772if (is_inside_out)773{774// Reverse winding775for (const uint8 *v = first_vtx + 2; v < end_vtx; ++v)776{777int i2 = *(v - 1);778int i3 = *v;779submerged_vol_calc.AddFace(i1, i3, i2);780}781}782else783{784// Normal winding785for (const uint8 *v = first_vtx + 2; v < end_vtx; ++v)786{787int i2 = *(v - 1);788int i3 = *v;789submerged_vol_calc.AddFace(i1, i2, i3);790}791}792}793794// Get the results795submerged_vol_calc.GetResult(outSubmergedVolume, outCenterOfBuoyancy);796}797798#ifdef JPH_DEBUG_RENDERER799// Draw center of buoyancy800if (sDrawSubmergedVolumes)801DebugRenderer::sInstance->DrawWireSphere(inBaseOffset + outCenterOfBuoyancy, 0.05f, Color::sRed, 1);802#endif // JPH_DEBUG_RENDERER803}804805#ifdef JPH_DEBUG_RENDERER806void ConvexHullShape::Draw(DebugRenderer *inRenderer, RMat44Arg inCenterOfMassTransform, Vec3Arg inScale, ColorArg inColor, bool inUseMaterialColors, bool inDrawWireframe) const807{808if (mGeometry == nullptr)809{810Array<DebugRenderer::Triangle> triangles;811for (const Face &f : mFaces)812{813const uint8 *first_vtx = mVertexIdx.data() + f.mFirstVertex;814const uint8 *end_vtx = first_vtx + f.mNumVertices;815816// Draw first triangle of polygon817Vec3 v0 = mPoints[first_vtx[0]].mPosition;818Vec3 v1 = mPoints[first_vtx[1]].mPosition;819Vec3 v2 = mPoints[first_vtx[2]].mPosition;820Vec3 uv_direction = (v1 - v0).Normalized();821triangles.push_back({ v0, v1, v2, Color::sWhite, v0, uv_direction });822823// Draw any other triangles in this polygon824for (const uint8 *v = first_vtx + 3; v < end_vtx; ++v)825triangles.push_back({ v0, mPoints[*(v - 1)].mPosition, mPoints[*v].mPosition, Color::sWhite, v0, uv_direction });826}827mGeometry = new DebugRenderer::Geometry(inRenderer->CreateTriangleBatch(triangles), GetLocalBounds());828}829830// Test if the shape is scaled inside out831DebugRenderer::ECullMode cull_mode = ScaleHelpers::IsInsideOut(inScale)? DebugRenderer::ECullMode::CullFrontFace : DebugRenderer::ECullMode::CullBackFace;832833// Determine the draw mode834DebugRenderer::EDrawMode draw_mode = inDrawWireframe? DebugRenderer::EDrawMode::Wireframe : DebugRenderer::EDrawMode::Solid;835836// Draw the geometry837Color color = inUseMaterialColors? GetMaterial()->GetDebugColor() : inColor;838RMat44 transform = inCenterOfMassTransform.PreScaled(inScale);839inRenderer->DrawGeometry(transform, color, mGeometry, cull_mode, DebugRenderer::ECastShadow::On, draw_mode);840841// Draw the outline if requested842if (sDrawFaceOutlines)843for (const Face &f : mFaces)844{845const uint8 *first_vtx = mVertexIdx.data() + f.mFirstVertex;846const uint8 *end_vtx = first_vtx + f.mNumVertices;847848// Draw edges of face849inRenderer->DrawLine(transform * mPoints[*(end_vtx - 1)].mPosition, transform * mPoints[*first_vtx].mPosition, Color::sGrey);850for (const uint8 *v = first_vtx + 1; v < end_vtx; ++v)851inRenderer->DrawLine(transform * mPoints[*(v - 1)].mPosition, transform * mPoints[*v].mPosition, Color::sGrey);852}853}854855void ConvexHullShape::DrawShrunkShape(DebugRenderer *inRenderer, RMat44Arg inCenterOfMassTransform, Vec3Arg inScale) const856{857// Get the shrunk points858SupportBuffer buffer;859const HullNoConvex *support = mConvexRadius > 0.0f? static_cast<const HullNoConvex *>(GetSupportFunction(ESupportMode::ExcludeConvexRadius, buffer, inScale)) : nullptr;860861RMat44 transform = inCenterOfMassTransform * Mat44::sScale(inScale);862863for (int p = 0; p < (int)mPoints.size(); ++p)864{865const Point &point = mPoints[p];866RVec3 position = transform * point.mPosition;867RVec3 shrunk_point = support != nullptr? transform * support->GetPoints()[p] : position;868869// Draw difference between shrunk position and position870inRenderer->DrawLine(position, shrunk_point, Color::sGreen);871872// Draw face normals that are contributing873for (int i = 0; i < point.mNumFaces; ++i)874inRenderer->DrawLine(position, position + 0.1f * mPlanes[point.mFaces[i]].GetNormal(), Color::sYellow);875876// Draw point index877inRenderer->DrawText3D(position, ConvertToString(p), Color::sWhite, 0.1f);878}879}880#endif // JPH_DEBUG_RENDERER881882bool ConvexHullShape::CastRayHelper(const RayCast &inRay, float &outMinFraction, float &outMaxFraction) const883{884if (mFaces.size() == 2)885{886// If we have only 2 faces, we're a flat convex hull and we need to test edges instead of planes887888// Check if plane is parallel to ray889const Plane &p = mPlanes.front();890Vec3 plane_normal = p.GetNormal();891float direction_projection = inRay.mDirection.Dot(plane_normal);892if (abs(direction_projection) >= 1.0e-12f)893{894// Calculate intersection point895float distance_to_plane = inRay.mOrigin.Dot(plane_normal) + p.GetConstant();896float fraction = -distance_to_plane / direction_projection;897if (fraction < 0.0f || fraction > 1.0f)898{899// Does not hit plane, no hit900outMinFraction = 0.0f;901outMaxFraction = 1.0f + FLT_EPSILON;902return false;903}904Vec3 intersection_point = inRay.mOrigin + fraction * inRay.mDirection;905906// Test all edges to see if point is inside polygon907const Face &f = mFaces.front();908const uint8 *first_vtx = mVertexIdx.data() + f.mFirstVertex;909const uint8 *end_vtx = first_vtx + f.mNumVertices;910Vec3 p1 = mPoints[*end_vtx].mPosition;911for (const uint8 *v = first_vtx; v < end_vtx; ++v)912{913Vec3 p2 = mPoints[*v].mPosition;914if ((p2 - p1).Cross(intersection_point - p1).Dot(plane_normal) < 0.0f)915{916// Outside polygon, no hit917outMinFraction = 0.0f;918outMaxFraction = 1.0f + FLT_EPSILON;919return false;920}921p1 = p2;922}923924// Inside polygon, a hit925outMinFraction = fraction;926outMaxFraction = fraction;927return true;928}929else930{931// Parallel ray doesn't hit932outMinFraction = 0.0f;933outMaxFraction = 1.0f + FLT_EPSILON;934return false;935}936}937else938{939// Clip ray against all planes940int fractions_set = 0;941bool all_inside = true;942float min_fraction = 0.0f, max_fraction = 1.0f + FLT_EPSILON;943for (const Plane &p : mPlanes)944{945// Check if the ray origin is behind this plane946Vec3 plane_normal = p.GetNormal();947float distance_to_plane = inRay.mOrigin.Dot(plane_normal) + p.GetConstant();948bool is_outside = distance_to_plane > 0.0f;949all_inside &= !is_outside;950951// Check if plane is parallel to ray952float direction_projection = inRay.mDirection.Dot(plane_normal);953if (abs(direction_projection) >= 1.0e-12f)954{955// Get intersection fraction between ray and plane956float fraction = -distance_to_plane / direction_projection;957958// Update interval of ray that is inside the hull959if (direction_projection < 0.0f)960{961min_fraction = max(fraction, min_fraction);962fractions_set |= 1;963}964else965{966max_fraction = min(fraction, max_fraction);967fractions_set |= 2;968}969}970else if (is_outside)971return false; // Outside the plane and parallel, no hit!972}973974// Test if both min and max have been set975if (fractions_set == 3)976{977// Output fractions978outMinFraction = min_fraction;979outMaxFraction = max_fraction;980981// Test if the infinite ray intersects with the hull (the length will be checked later)982return min_fraction <= max_fraction && max_fraction >= 0.0f;983}984else985{986// Degenerate case, either the ray is parallel to all planes or the ray has zero length987outMinFraction = 0.0f;988outMaxFraction = 1.0f + FLT_EPSILON;989990// Return if the origin is inside the hull991return all_inside;992}993}994}995996bool ConvexHullShape::CastRay(const RayCast &inRay, const SubShapeIDCreator &inSubShapeIDCreator, RayCastResult &ioHit) const997{998// Determine if ray hits the shape999float min_fraction, max_fraction;1000if (CastRayHelper(inRay, min_fraction, max_fraction)1001&& min_fraction < ioHit.mFraction) // Check if this is a closer hit1002{1003// Better hit than the current hit1004ioHit.mFraction = min_fraction;1005ioHit.mSubShapeID2 = inSubShapeIDCreator.GetID();1006return true;1007}1008return false;1009}10101011void ConvexHullShape::CastRay(const RayCast &inRay, const RayCastSettings &inRayCastSettings, const SubShapeIDCreator &inSubShapeIDCreator, CastRayCollector &ioCollector, const ShapeFilter &inShapeFilter) const1012{1013// Test shape filter1014if (!inShapeFilter.ShouldCollide(this, inSubShapeIDCreator.GetID()))1015return;10161017// Determine if ray hits the shape1018float min_fraction, max_fraction;1019if (CastRayHelper(inRay, min_fraction, max_fraction)1020&& min_fraction < ioCollector.GetEarlyOutFraction()) // Check if this is closer than the early out fraction1021{1022// Better hit than the current hit1023RayCastResult hit;1024hit.mBodyID = TransformedShape::sGetBodyID(ioCollector.GetContext());1025hit.mSubShapeID2 = inSubShapeIDCreator.GetID();10261027// Check front side hit1028if (inRayCastSettings.mTreatConvexAsSolid || min_fraction > 0.0f)1029{1030hit.mFraction = min_fraction;1031ioCollector.AddHit(hit);1032}10331034// Check back side hit1035if (inRayCastSettings.mBackFaceModeConvex == EBackFaceMode::CollideWithBackFaces1036&& max_fraction < ioCollector.GetEarlyOutFraction())1037{1038hit.mFraction = max_fraction;1039ioCollector.AddHit(hit);1040}1041}1042}10431044void ConvexHullShape::CollidePoint(Vec3Arg inPoint, const SubShapeIDCreator &inSubShapeIDCreator, CollidePointCollector &ioCollector, const ShapeFilter &inShapeFilter) const1045{1046// Test shape filter1047if (!inShapeFilter.ShouldCollide(this, inSubShapeIDCreator.GetID()))1048return;10491050// Check if point is behind all planes1051for (const Plane &p : mPlanes)1052if (p.SignedDistance(inPoint) > 0.0f)1053return;10541055// Point is inside1056ioCollector.AddHit({ TransformedShape::sGetBodyID(ioCollector.GetContext()), inSubShapeIDCreator.GetID() });1057}10581059void ConvexHullShape::CollideSoftBodyVertices(Mat44Arg inCenterOfMassTransform, Vec3Arg inScale, const CollideSoftBodyVertexIterator &inVertices, uint inNumVertices, int inCollidingShapeIndex) const1060{1061Mat44 inverse_transform = inCenterOfMassTransform.InversedRotationTranslation();10621063Vec3 inv_scale = inScale.Reciprocal();1064bool is_not_scaled = ScaleHelpers::IsNotScaled(inScale);1065float scale_flip = ScaleHelpers::IsInsideOut(inScale)? -1.0f : 1.0f;10661067for (CollideSoftBodyVertexIterator v = inVertices, sbv_end = inVertices + inNumVertices; v != sbv_end; ++v)1068if (v.GetInvMass() > 0.0f)1069{1070Vec3 local_pos = inverse_transform * v.GetPosition();10711072// Find most facing plane1073float max_distance = -FLT_MAX;1074Vec3 max_plane_normal = Vec3::sZero();1075uint max_plane_idx = 0;1076if (is_not_scaled)1077{1078// Without scale, it is trivial to calculate the distance to the hull1079for (const Plane &p : mPlanes)1080{1081float distance = p.SignedDistance(local_pos);1082if (distance > max_distance)1083{1084max_distance = distance;1085max_plane_normal = p.GetNormal();1086max_plane_idx = uint(&p - mPlanes.data());1087}1088}1089}1090else1091{1092// When there's scale we need to calculate the planes first1093for (uint i = 0; i < (uint)mPlanes.size(); ++i)1094{1095// Calculate plane normal and point by scaling the original plane1096Vec3 plane_normal = (inv_scale * mPlanes[i].GetNormal()).Normalized();1097Vec3 plane_point = inScale * mPoints[mVertexIdx[mFaces[i].mFirstVertex]].mPosition;10981099float distance = plane_normal.Dot(local_pos - plane_point);1100if (distance > max_distance)1101{1102max_distance = distance;1103max_plane_normal = plane_normal;1104max_plane_idx = i;1105}1106}1107}1108bool is_outside = max_distance > 0.0f;11091110// Project point onto that plane1111Vec3 closest_point = local_pos - max_distance * max_plane_normal;11121113// Check edges if we're outside the hull (when inside we know the closest face is also the closest point to the surface)1114if (is_outside)1115{1116// Loop over edges1117float closest_point_dist_sq = FLT_MAX;1118const Face &face = mFaces[max_plane_idx];1119for (const uint8 *v_start = &mVertexIdx[face.mFirstVertex], *v1 = v_start, *v_end = v_start + face.mNumVertices; v1 < v_end; ++v1)1120{1121// Find second point1122const uint8 *v2 = v1 + 1;1123if (v2 == v_end)1124v2 = v_start;11251126// Get edge points1127Vec3 p1 = inScale * mPoints[*v1].mPosition;1128Vec3 p2 = inScale * mPoints[*v2].mPosition;11291130// Check if the position is outside the edge (if not, the face will be closer)1131Vec3 edge_normal = (p2 - p1).Cross(max_plane_normal);1132if (scale_flip * edge_normal.Dot(local_pos - p1) > 0.0f)1133{1134// Get closest point on edge1135uint32 set;1136Vec3 closest = ClosestPoint::GetClosestPointOnLine(p1 - local_pos, p2 - local_pos, set);1137float distance_sq = closest.LengthSq();1138if (distance_sq < closest_point_dist_sq)1139closest_point = local_pos + closest;1140}1141}1142}11431144// Check if this is the largest penetration1145Vec3 normal = local_pos - closest_point;1146float normal_length = normal.Length();1147float penetration = normal_length;1148if (is_outside)1149penetration = -penetration;1150else1151normal = -normal;1152if (v.UpdatePenetration(penetration))1153{1154// Calculate contact plane1155normal = normal_length > 0.0f? normal / normal_length : max_plane_normal;1156Plane plane = Plane::sFromPointAndNormal(closest_point, normal);11571158// Store collision1159v.SetCollision(plane.GetTransformed(inCenterOfMassTransform), inCollidingShapeIndex);1160}1161}1162}11631164class ConvexHullShape::CHSGetTrianglesContext1165{1166public:1167CHSGetTrianglesContext(Mat44Arg inTransform, bool inIsInsideOut) : mTransform(inTransform), mIsInsideOut(inIsInsideOut) { }11681169Mat44 mTransform;1170bool mIsInsideOut;1171size_t mCurrentFace = 0;1172};11731174void ConvexHullShape::GetTrianglesStart(GetTrianglesContext &ioContext, const AABox &inBox, Vec3Arg inPositionCOM, QuatArg inRotation, Vec3Arg inScale) const1175{1176static_assert(sizeof(CHSGetTrianglesContext) <= sizeof(GetTrianglesContext), "GetTrianglesContext too small");1177JPH_ASSERT(IsAligned(&ioContext, alignof(CHSGetTrianglesContext)));11781179new (&ioContext) CHSGetTrianglesContext(Mat44::sRotationTranslation(inRotation, inPositionCOM) * Mat44::sScale(inScale), ScaleHelpers::IsInsideOut(inScale));1180}11811182int ConvexHullShape::GetTrianglesNext(GetTrianglesContext &ioContext, int inMaxTrianglesRequested, Float3 *outTriangleVertices, const PhysicsMaterial **outMaterials) const1183{1184static_assert(cGetTrianglesMinTrianglesRequested >= 12, "cGetTrianglesMinTrianglesRequested is too small");1185JPH_ASSERT(inMaxTrianglesRequested >= cGetTrianglesMinTrianglesRequested);11861187CHSGetTrianglesContext &context = (CHSGetTrianglesContext &)ioContext;11881189int total_num_triangles = 0;1190for (; context.mCurrentFace < mFaces.size(); ++context.mCurrentFace)1191{1192const Face &f = mFaces[context.mCurrentFace];11931194const uint8 *first_vtx = mVertexIdx.data() + f.mFirstVertex;1195const uint8 *end_vtx = first_vtx + f.mNumVertices;11961197// Check if there is still room in the output buffer for this face1198int num_triangles = f.mNumVertices - 2;1199inMaxTrianglesRequested -= num_triangles;1200if (inMaxTrianglesRequested < 0)1201break;1202total_num_triangles += num_triangles;12031204// Get first triangle of polygon1205Vec3 v0 = context.mTransform * mPoints[first_vtx[0]].mPosition;1206Vec3 v1 = context.mTransform * mPoints[first_vtx[1]].mPosition;1207Vec3 v2 = context.mTransform * mPoints[first_vtx[2]].mPosition;1208v0.StoreFloat3(outTriangleVertices++);1209if (context.mIsInsideOut)1210{1211// Store first triangle in this polygon flipped1212v2.StoreFloat3(outTriangleVertices++);1213v1.StoreFloat3(outTriangleVertices++);12141215// Store other triangles in this polygon flipped1216for (const uint8 *v = first_vtx + 3; v < end_vtx; ++v)1217{1218v0.StoreFloat3(outTriangleVertices++);1219(context.mTransform * mPoints[*v].mPosition).StoreFloat3(outTriangleVertices++);1220(context.mTransform * mPoints[*(v - 1)].mPosition).StoreFloat3(outTriangleVertices++);1221}1222}1223else1224{1225// Store first triangle in this polygon1226v1.StoreFloat3(outTriangleVertices++);1227v2.StoreFloat3(outTriangleVertices++);12281229// Store other triangles in this polygon1230for (const uint8 *v = first_vtx + 3; v < end_vtx; ++v)1231{1232v0.StoreFloat3(outTriangleVertices++);1233(context.mTransform * mPoints[*(v - 1)].mPosition).StoreFloat3(outTriangleVertices++);1234(context.mTransform * mPoints[*v].mPosition).StoreFloat3(outTriangleVertices++);1235}1236}1237}12381239// Store materials1240if (outMaterials != nullptr)1241{1242const PhysicsMaterial *material = GetMaterial();1243for (const PhysicsMaterial **m = outMaterials, **m_end = outMaterials + total_num_triangles; m < m_end; ++m)1244*m = material;1245}12461247return total_num_triangles;1248}12491250void ConvexHullShape::SaveBinaryState(StreamOut &inStream) const1251{1252ConvexShape::SaveBinaryState(inStream);12531254inStream.Write(mCenterOfMass);1255inStream.Write(mInertia);1256inStream.Write(mLocalBounds.mMin);1257inStream.Write(mLocalBounds.mMax);1258inStream.Write(mPoints);1259inStream.Write(mFaces);1260inStream.Write(mPlanes);1261inStream.Write(mVertexIdx);1262inStream.Write(mConvexRadius);1263inStream.Write(mVolume);1264inStream.Write(mInnerRadius);1265}12661267void ConvexHullShape::RestoreBinaryState(StreamIn &inStream)1268{1269ConvexShape::RestoreBinaryState(inStream);12701271inStream.Read(mCenterOfMass);1272inStream.Read(mInertia);1273inStream.Read(mLocalBounds.mMin);1274inStream.Read(mLocalBounds.mMax);1275inStream.Read(mPoints);1276inStream.Read(mFaces);1277inStream.Read(mPlanes);1278inStream.Read(mVertexIdx);1279inStream.Read(mConvexRadius);1280inStream.Read(mVolume);1281inStream.Read(mInnerRadius);1282}12831284Shape::Stats ConvexHullShape::GetStats() const1285{1286// Count number of triangles1287uint triangle_count = 0;1288for (const Face &f : mFaces)1289triangle_count += f.mNumVertices - 2;12901291return Stats(1292sizeof(*this)1293+ mPoints.size() * sizeof(Point)1294+ mFaces.size() * sizeof(Face)1295+ mPlanes.size() * sizeof(Plane)1296+ mVertexIdx.size() * sizeof(uint8),1297triangle_count);1298}12991300void ConvexHullShape::sRegister()1301{1302ShapeFunctions &f = ShapeFunctions::sGet(EShapeSubType::ConvexHull);1303f.mConstruct = []() -> Shape * { return new ConvexHullShape; };1304f.mColor = Color::sGreen;1305}13061307JPH_NAMESPACE_END130813091310