ExtDelaunayTetrahedralizer.h

// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of NVIDIA CORPORATION nor the names of its
//    contributors may be used to endorse or promote products derived
//    from this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ''AS IS'' AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
// OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Copyright (c) 2008-2022 NVIDIA Corporation. All rights reserved.
 
#ifndef EXT_DELAUNAY_TETRAHEDRALIZER_H
#define EXT_DELAUNAY_TETRAHEDRALIZER_H
 
#include "foundation/PxArray.h"
#include "foundation/PxVec3.h"
#include "foundation/PxHashSet.h"
#include "GuTetrahedron.h"
#include "ExtVec3.h"
 
namespace physx
{   
 
namespace Ext
{
    using Edge = PxPair<PxI32, PxI32>;
    using Tetrahedron = Gu::TetrahedronT<PxI32>;
    using Tetrahedron16 = Gu::TetrahedronT<PxI16>;
 
    void buildNeighborhood(const PxArray<Tetrahedron>& tets, PxArray<PxI32>& result);
    void buildNeighborhood(const PxI32* tets, PxU32 numTets, PxArray<PxI32>& result);
 
    PX_FORCE_INLINE PxF64 tetVolume(const Vec3& a, const Vec3& b, const Vec3& c, const Vec3& d)
    {
        return (-1.0 / 6.0) * (a - d).dot((b - d).cross(c - d));
    }
 
    PX_FORCE_INLINE PxF64 tetVolume(const Tetrahedron& tet, const PxArray<Vec3>& points)
    {
        return tetVolume(points[tet[0]], points[tet[1]], points[tet[2]], points[tet[3]]);
    }
 
    //Returns the intersection (set operation) of two sorted lists
    PX_FORCE_INLINE void intersectionOfSortedLists(const PxArray<PxI32>& sorted1, const PxArray<PxI32>& sorted2, PxArray<PxI32>& result)
    {
        PxU32 a = 0;
        PxU32 b = 0;
        result.clear();
        while (a < sorted1.size() && b < sorted2.size())
        {
            if (sorted1[a] == sorted2[b])
            {
                result.pushBack(sorted1[a]);
                ++a;
                ++b;
            }
            else if (sorted1[a] > sorted2[b])
                ++b;
            else
                ++a;
        }
    }
    
    PX_FORCE_INLINE bool intersectionOfSortedListsContainsElements(const PxArray<PxI32>& sorted1, const PxArray<PxI32>& sorted2)
    {
        PxU32 a = 0;
        PxU32 b = 0;
        while (a < sorted1.size() && b < sorted2.size())
        {
            if (sorted1[a] == sorted2[b])           
                return true;            
            else if (sorted1[a] > sorted2[b])
                ++b;
            else
                ++a;
        }
        return false;
    }
 
 
    class BaseTetAnalyzer
    {
    public:
        virtual PxF64 quality(const PxArray<PxI32> tetIndices) const = 0;
 
        virtual PxF64 quality(const PxArray<Tetrahedron> tetrahedraToCheck) const = 0;
 
        virtual bool improved(PxF64 previousQuality, PxF64 newQuality) const = 0;
 
        virtual ~BaseTetAnalyzer() {}
    };
 
    class MinimizeMaxAmipsEnergy : public BaseTetAnalyzer
    {
        const PxArray<Vec3>& points;
        const PxArray<Tetrahedron>& tetrahedra;
 
    public:
        MinimizeMaxAmipsEnergy(const PxArray<Vec3>& points_, const PxArray<Tetrahedron>& tetrahedra_) : points(points_), tetrahedra(tetrahedra_)
        {}
 
        PxF64 quality(const PxArray<PxI32> tetIndices) const;
 
        PxF64 quality(const PxArray<Tetrahedron> tetrahedraToCheck) const;
 
        bool improved(PxF64 previousQuality, PxF64 newQuality) const;
 
        virtual ~MinimizeMaxAmipsEnergy() {}
 
    private:
        PX_NOCOPY(MinimizeMaxAmipsEnergy)
    };
 
    class MaximizeMinTetVolume : public BaseTetAnalyzer
    {
        const PxArray<Vec3>& points;
        const PxArray<Tetrahedron>& tetrahedra;
 
    public:
        MaximizeMinTetVolume(const PxArray<Vec3>& points_, const PxArray<Tetrahedron>& tetrahedra_) : points(points_), tetrahedra(tetrahedra_)
        {}
 
        PxF64 quality(const PxArray<PxI32> tetIndices) const;
 
        PxF64 quality(const PxArray<Tetrahedron> tetrahedraToCheck) const;
 
        bool improved(PxF64 previousQuality, PxF64 newQuality) const;
 
        virtual ~MaximizeMinTetVolume() {}
 
    private:
        PX_NOCOPY(MaximizeMinTetVolume)
    };
 
 
 
    //Helper class to extract surface triangles from a tetmesh
    struct SortedTriangle
    {
    public:
        PxI32 A;
        PxI32 B;
        PxI32 C;
        bool Flipped;
 
        PX_FORCE_INLINE SortedTriangle(PxI32 a, PxI32 b, PxI32 c)
        {
            A = a; B = b; C = c; Flipped = false;
            if (A > B) { PxSwap(A, B); Flipped = !Flipped; }
            if (B > C) { PxSwap(B, C); Flipped = !Flipped; }
            if (A > B) { PxSwap(A, B); Flipped = !Flipped; }
        }
    };
 
    struct TriangleHash
    {
        PX_FORCE_INLINE std::size_t operator()(const SortedTriangle& k) const
        {
            return k.A ^ k.B ^ k.C;
        }
 
        PX_FORCE_INLINE bool equal(const SortedTriangle& first, const SortedTriangle& second) const
        {
            return first.A == second.A && first.B == second.B && first.C == second.C;
        }
    };
 
 
    struct RecoverEdgeMemoryCache
    {
        PxArray<PxI32> facesStart;
        PxHashSet<PxI32> tetsDoneStart;
        PxArray<PxI32> resultStart;
        PxArray<PxI32> facesEnd;
        PxHashSet<PxI32> tetsDoneEnd;
        PxArray<PxI32> resultEnd;
    };
 
    class StackMemory
    {
    public:
 
        void clear()
        {
            faces.forceSize_Unsafe(0);
            hashSet.clear();
        }
 
        PxArray<PxI32> faces;
        PxHashSet<PxI32> hashSet;
    
    };
 
    //Incremental delaunay tetrahedralizer
    class DelaunayTetrahedralizer
    {
    public:
        //The bounds specified must contain all points that will get inserted by calling insertPoints.
        DelaunayTetrahedralizer(const Vec3& min, const Vec3& max);
 
        DelaunayTetrahedralizer(PxArray<Vec3>& points, PxArray<Tetrahedron>& tets);
 
        void initialize(PxArray<Vec3>& points, PxArray<Tetrahedron>& tets);
 
        //Inserts a bunch of new points into the tetrahedralization and keeps the delaunay condition satisfied. The new result will
        //get stored in the tetrahedra array. Points to insert must already be present in inPoints, the indices of the points to insert
        //can be controlled with start and end index (end index is exclusive, start index is inclusive)
        void insertPoints(const PxArray<Vec3>& inPoints, PxI32 start, PxI32 end, PxArray<Tetrahedron>& tetrahedra);
        
        bool insertPoints(const PxArray<Vec3>& inPoints, PxI32 start, PxI32 end);
 
        void exportTetrahedra(PxArray<Tetrahedron>& tetrahedra);
 
        bool canCollapseEdge(PxI32 edgeVertexToKeep, PxI32 edgeVertexToRemove, PxF64 volumeChangeThreshold = 0.1, BaseTetAnalyzer* tetAnalyzer = NULL);
        bool canCollapseEdge(PxI32 edgeVertexToKeep, PxI32 edgeVertexToRemove, const PxArray<PxI32>& tetsConnectedToA, const PxArray<PxI32>& tetsConnectedToB,
            PxF64& qualityAfterCollapse, PxF64 volumeChangeThreshold = 0.1, BaseTetAnalyzer* tetAnalyzer = NULL);
 
        void collapseEdge(PxI32 edgeVertexToKeep, PxI32 edgeVertexToRemove);
 
        void collapseEdge(PxI32 edgeVertexAToKeep, PxI32 edgeVertexBToRemove, const PxArray<PxI32>& tetsConnectedToA, const PxArray<PxI32>& tetsConnectedToB);
 
        void collectTetsConnectedToVertex(PxI32 vertexIndex, PxArray<PxI32>& tetIds);
 
        void collectTetsConnectedToVertex(PxArray<PxI32>& faces, PxHashSet<PxI32>& tetsDone, PxI32 vertexIndex, PxArray<PxI32>& tetIds);
 
        void collectTetsConnectedToEdge(PxI32 edgeStart, PxI32 edgeEnd, PxArray<PxI32>& tetIds);
 
        PX_FORCE_INLINE const Vec3& point(PxI32 index) const { return centeredNormalizedPoints[index]; }
 
        PX_FORCE_INLINE PxU32 numPoints() const { return centeredNormalizedPoints.size(); }
 
        PX_FORCE_INLINE PxArray<Vec3>& points() { return centeredNormalizedPoints; }
        PX_FORCE_INLINE const PxArray<Vec3>& points() const { return centeredNormalizedPoints; }
 
        PxU32 addPoint(const Vec3& p)
        {
            centeredNormalizedPoints.pushBack(p);
            vertexToTet.pushBack(-1);
            return centeredNormalizedPoints.size() - 1;
        }
 
        PX_FORCE_INLINE const Tetrahedron& tetrahedron(PxI32 index) const { return result[index]; }
 
        PX_FORCE_INLINE PxU32 numTetrahedra() const { return result.size(); }
 
        PX_FORCE_INLINE PxArray<Tetrahedron>& tetrahedra() { return result; }
        PX_FORCE_INLINE const PxArray<Tetrahedron>& tetrahedra() const { return result; }
 
        void copyInternalPointsTo(PxArray<Vec3>& points) { points = centeredNormalizedPoints; }
 
        bool optimizeByFlipping(PxArray<PxI32>& affectedFaces, const BaseTetAnalyzer& qualityAnalyzer);
 
        void insertPointIntoEdge(PxI32 newPointIndex, PxI32 edgeA, PxI32 edgeB, PxArray<PxI32>& affectedTets, BaseTetAnalyzer* qualityAnalyzer = NULL);
 
        bool removeEdgeByFlip(PxI32 edgeA, PxI32 edgeB, PxArray<PxI32>& tetIndices, BaseTetAnalyzer* qualityAnalyzer = NULL);
 
        void addLockedEdges(const PxArray<Gu::IndexedTriangleT<PxI32>>& triangles);
 
        void addLockedTriangles(const PxArray<Gu::IndexedTriangleT<PxI32>>& triangles);
 
        void clearLockedEdges() { lockedEdges.clear(); }
 
        void clearLockedTriangles() { lockedTriangles.clear(); }
 
        bool recoverEdgeByFlip(PxI32 eStart, PxI32 eEnd, RecoverEdgeMemoryCache& cache);
 
        void generateTetmeshEnforcingEdges(const PxArray<Vec3>& trianglePoints, const PxArray<Gu::IndexedTriangleT<PxI32>>& triangles, PxArray<PxArray<PxI32>>& allEdges,
            PxArray<PxArray<PxI32>>& pointToOriginalTriangle, 
            PxArray<Vec3>& points, PxArray<Tetrahedron>& finalTets);
 
    private:
        PxArray<Vec3> centeredNormalizedPoints;
        PxArray<PxI32> neighbors;
        PxArray<PxI32> unusedTets;
        PxArray<PxI32> vertexToTet;
        PxArray<Tetrahedron> result;
        PxI32 numAdditionalPointsAtBeginning = 4;
 
        PxHashSet<SortedTriangle, TriangleHash> lockedTriangles;
        PxHashSet<PxU64> lockedEdges;
 
        StackMemory     stackMemory;
        PX_NOCOPY(DelaunayTetrahedralizer)
    };
 
    struct EdgeWithLength
    {
        PxI32 A;
        PxI32 B;
        PxF64 Length;
 
        EdgeWithLength(PxI32 a_, PxI32 b_, PxF64 length_)
        {
            A = a_;
            B = b_;
            Length = length_;
        }       
    };
 
    PX_FORCE_INLINE bool operator <(const EdgeWithLength& lhs, const EdgeWithLength& rhs)
    {
        return lhs.Length < rhs.Length;
    }
 
    struct SplitEdge
    {
        PxI32 A;
        PxI32 B;
        PxF64 Q;
        PxF64 L;
        bool InteriorEdge;
 
        SplitEdge(PxI32 a, PxI32 b, PxF64 q, PxF64 l, bool interiorEdge)
        {
            A = a;
            B = b;
            Q = q;
            L = l;
            InteriorEdge = interiorEdge;
        }       
    };
 
    PX_FORCE_INLINE bool operator >(const SplitEdge& lhs, const SplitEdge& rhs)
    {
        if (lhs.Q == rhs.Q)
            return lhs.L > rhs.L;
        return lhs.Q > rhs.Q;
    }
    
    
    bool optimizeByCollapsing(DelaunayTetrahedralizer& del, const PxArray<EdgeWithLength>& edges,
        PxArray<PxArray<PxI32>>& pointToOriginalTriangle, PxI32 numFixPoints, BaseTetAnalyzer* qualityAnalyzer = NULL);
 
    bool optimizeBySwapping(DelaunayTetrahedralizer& del, const PxArray<EdgeWithLength>& edges,
        const PxArray<PxArray<PxI32>>& pointToOriginalTriangle, BaseTetAnalyzer* qualityAnalyzer);
 
    bool optimizeBySplitting(DelaunayTetrahedralizer& del, const PxArray<EdgeWithLength>& edges, const PxArray<PxArray<PxI32>>& pointToOriginalTriangle,
        PxI32 maxPointsToInsert = -1, bool sortByQuality = false, BaseTetAnalyzer* qualityAnalyzer = NULL, PxF64 qualityThreshold = 10);
 
    //Modified tetmesh quality improvement implementation of the method described in https://cs.nyu.edu/~yixinhu/tetwild.pdf Section 3.2 Mesh Improvement
    void optimize(DelaunayTetrahedralizer& del, PxArray<PxArray<PxI32>>& pointToOriginalTriangle, PxI32 numFixPoints,
        PxArray<Vec3>& optimizedPoints, PxArray<Tetrahedron>& optimizedTets, PxI32 numPasses = 10);
}
}
 
#endif