From 057d4b15e3c539cc52e6381f2b61e0d23985e4a0 Mon Sep 17 00:00:00 2001 From: oan Date: Wed, 10 Jul 2019 13:04:25 +0300 Subject: [PATCH] 0028089: Mesh - New algorithm for triangulation of 2d polygons Added custom meshing core algorithm to generate base mesh using Delabella library: https://github.com/msokalski/delabella # Added CSF_MeshAlgo environment variable to select meshing core without necessity of recompilation --- src/BRepMesh/BRepMesh_Context.cxx | 26 +- src/BRepMesh/BRepMesh_CustomBaseMeshAlgo.hxx | 23 + .../BRepMesh_DelabellaBaseMeshAlgo.cxx | 155 +++ .../BRepMesh_DelabellaBaseMeshAlgo.hxx | 46 + .../BRepMesh_DelabellaMeshAlgoFactory.cxx | 143 +++ .../BRepMesh_DelabellaMeshAlgoFactory.hxx | 44 + src/BRepMesh/BRepMesh_Delaun.cxx | 74 +- src/BRepMesh/BRepMesh_Delaun.hxx | 27 +- src/BRepMesh/DELABELLA_LICENSE | 22 + src/BRepMesh/FILES | 6 + src/BRepMesh/delabella.cpp | 1021 +++++++++++++++++ src/BRepMesh/delabella.h | 68 ++ 12 files changed, 1624 insertions(+), 31 deletions(-) create mode 100644 src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.cxx create mode 100644 src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.hxx create mode 100644 src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.cxx create mode 100644 src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.hxx create mode 100644 src/BRepMesh/DELABELLA_LICENSE create mode 100644 src/BRepMesh/delabella.cpp create mode 100644 src/BRepMesh/delabella.h diff --git a/src/BRepMesh/BRepMesh_Context.cxx b/src/BRepMesh/BRepMesh_Context.cxx index e48d038467..188b924175 100644 --- a/src/BRepMesh/BRepMesh_Context.cxx +++ b/src/BRepMesh/BRepMesh_Context.cxx @@ -20,7 +20,9 @@ #include #include #include + #include +#include //======================================================================= // Function: Constructor @@ -28,11 +30,33 @@ //======================================================================= BRepMesh_Context::BRepMesh_Context () { + enum MeshAlgo + { + MeshAlgo_Default = 0x0, + MeshAlgo_Delabella = 0x1 + }; + + char* anAlgoVar; + anAlgoVar = getenv ("CSF_MeshAlgo"); + const Standard_Integer anAlgoId = (anAlgoVar ? atoi (anAlgoVar) : MeshAlgo_Default); + + Handle (IMeshTools_MeshAlgoFactory) aAlgoFactory; + switch (anAlgoId) + { + case MeshAlgo_Delabella: + aAlgoFactory = new BRepMesh_DelabellaMeshAlgoFactory; + break; + + default: + aAlgoFactory = new BRepMesh_MeshAlgoFactory; + break; + } + SetModelBuilder (new BRepMesh_ModelBuilder); SetEdgeDiscret (new BRepMesh_EdgeDiscret); SetModelHealer (new BRepMesh_ModelHealer); SetPreProcessor (new BRepMesh_ModelPreProcessor); - SetFaceDiscret (new BRepMesh_FaceDiscret(new BRepMesh_MeshAlgoFactory)); + SetFaceDiscret (new BRepMesh_FaceDiscret (aAlgoFactory)); SetPostProcessor(new BRepMesh_ModelPostProcessor); } diff --git a/src/BRepMesh/BRepMesh_CustomBaseMeshAlgo.hxx b/src/BRepMesh/BRepMesh_CustomBaseMeshAlgo.hxx index e175091de9..3526eefbd1 100644 --- a/src/BRepMesh/BRepMesh_CustomBaseMeshAlgo.hxx +++ b/src/BRepMesh/BRepMesh_CustomBaseMeshAlgo.hxx @@ -49,12 +49,35 @@ protected: Standard_EXPORT virtual void generateMesh () Standard_OVERRIDE { const Handle (BRepMesh_DataStructureOfDelaun)& aStructure = this->getStructure (); + const Standard_Integer aNodesNb = aStructure->NbNodes (); + buildBaseTriangulation (); std::pair aCellsCount = this->getCellsCount (aStructure->NbNodes ()); BRepMesh_Delaun aMesher (aStructure, aCellsCount.first, aCellsCount.second, Standard_False); + + const Standard_Integer aNewNodesNb = aStructure->NbNodes (); + const Standard_Boolean isRemoveAux = aNewNodesNb > aNodesNb; + if (isRemoveAux) + { + IMeshData::VectorOfInteger aAuxVertices (aNewNodesNb - aNodesNb); + for (Standard_Integer aExtNodesIt = aNodesNb + 1; aExtNodesIt <= aNewNodesNb; ++aExtNodesIt) + { + aAuxVertices.Append (aExtNodesIt); + } + + // Set aux vertices if there are some to clean up mesh correctly. + aMesher.SetAuxVertices (aAuxVertices); + } + aMesher.ProcessConstraints (); + // Destruction of triangles containing aux vertices added (possibly) during base mesh computation. + if (isRemoveAux) + { + aMesher.RemoveAuxElements (); + } + BRepMesh_MeshTool aCleaner (aStructure); aCleaner.EraseFreeLinks (); diff --git a/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.cxx b/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.cxx new file mode 100644 index 0000000000..1a03fb8bf3 --- /dev/null +++ b/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.cxx @@ -0,0 +1,155 @@ +// Created on: 2019-07-05 +// Copyright (c) 2019 OPEN CASCADE SAS +// Created by: Oleg AGASHIN +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include +#include +#include + +#include "delabella.h" + +//======================================================================= +// Function: Constructor +// Purpose : +//======================================================================= +BRepMesh_DelabellaBaseMeshAlgo::BRepMesh_DelabellaBaseMeshAlgo () +{ +} + +//======================================================================= +// Function: Destructor +// Purpose : +//======================================================================= +BRepMesh_DelabellaBaseMeshAlgo::~BRepMesh_DelabellaBaseMeshAlgo () +{ +} + +//======================================================================= +//function : buildBaseTriangulation +//purpose : +//======================================================================= +void BRepMesh_DelabellaBaseMeshAlgo::buildBaseTriangulation() +{ + const Handle(BRepMesh_DataStructureOfDelaun)& aStructure = this->getStructure(); + + Bnd_B2d aBox; + const Standard_Integer aNodesNb = aStructure->NbNodes (); + std::vector aPoints (2 * (aNodesNb + 4)); + for (Standard_Integer aNodeIt = 0; aNodeIt < aNodesNb; ++aNodeIt) + { + const BRepMesh_Vertex& aVertex = aStructure->GetNode (aNodeIt + 1); + + const size_t aBaseIdx = 2 * static_cast (aNodeIt); + aPoints[aBaseIdx + 0] = aVertex.Coord ().X (); + aPoints[aBaseIdx + 1] = aVertex.Coord ().Y (); + + aBox.Add (gp_Pnt2d(aVertex.Coord ())); + } + + aBox.Enlarge (0.1 * (aBox.CornerMax () - aBox.CornerMin ()).Modulus ()); + const gp_XY aMin = aBox.CornerMin (); + const gp_XY aMax = aBox.CornerMax (); + + aPoints[2 * aNodesNb + 0] = aMin.X (); + aPoints[2 * aNodesNb + 1] = aMin.Y (); + aStructure->AddNode (BRepMesh_Vertex ( + aPoints[2 * aNodesNb + 0], + aPoints[2 * aNodesNb + 1], BRepMesh_Free)); + + aPoints[2 * aNodesNb + 2] = aMax.X (); + aPoints[2 * aNodesNb + 3] = aMin.Y (); + aStructure->AddNode (BRepMesh_Vertex ( + aPoints[2 * aNodesNb + 2], + aPoints[2 * aNodesNb + 3], BRepMesh_Free)); + + aPoints[2 * aNodesNb + 4] = aMax.X (); + aPoints[2 * aNodesNb + 5] = aMax.Y (); + aStructure->AddNode (BRepMesh_Vertex ( + aPoints[2 * aNodesNb + 4], + aPoints[2 * aNodesNb + 5], BRepMesh_Free)); + + aPoints[2 * aNodesNb + 6] = aMin.X (); + aPoints[2 * aNodesNb + 7] = aMax.Y (); + aStructure->AddNode (BRepMesh_Vertex ( + aPoints[2 * aNodesNb + 6], + aPoints[2 * aNodesNb + 7], BRepMesh_Free)); + + const Standard_Real aDiffX = (aMax.X () - aMin.X ()); + const Standard_Real aDiffY = (aMax.Y () - aMin.Y ()); + for (size_t i = 0; i < aPoints.size(); i += 2) + { + aPoints[i + 0] = (aPoints[i + 0] - aMin.X ()) / aDiffX - 0.5; + aPoints[i + 1] = (aPoints[i + 1] - aMin.Y ()) / aDiffY - 0.5; + } + + IDelaBella* aTriangulator = IDelaBella::Create(); + try + { + if (aTriangulator != NULL) + { + const int aVerticesNb = aTriangulator->Triangulate ( + static_cast(aPoints.size () / 2), + &aPoints[0], &aPoints[1], 2 * sizeof (Standard_Real)); + + if (aVerticesNb > 0) + { + const DelaBella_Triangle* aTrianglePtr = aTriangulator->GetFirstDelaunayTriangle(); + while (aTrianglePtr != NULL) + { + Standard_Integer aNodes[3] = { + aTrianglePtr->v[0]->i + 1, + aTrianglePtr->v[2]->i + 1, + aTrianglePtr->v[1]->i + 1 + }; + + Standard_Integer aEdges [3]; + Standard_Boolean aOrientations[3]; + for (Standard_Integer k = 0; k < 3; ++k) + { + const BRepMesh_Edge aLink (aNodes[k], aNodes[(k + 1) % 3], BRepMesh_Free); + + const Standard_Integer aLinkInfo = aStructure->AddLink (aLink); + aEdges [k] = Abs (aLinkInfo); + aOrientations[k] = aLinkInfo > 0; + } + + const BRepMesh_Triangle aTriangle (aEdges, aOrientations, BRepMesh_Free); + aStructure->AddElement (aTriangle); + + aTrianglePtr = aTrianglePtr->next; + } + } + + aTriangulator->Destroy (); + } + } + catch (Standard_Failure const& theException) + { + if (aTriangulator != NULL) + { + aTriangulator->Destroy (); + } + + throw Standard_Failure (theException); + } + catch (...) + { + if (aTriangulator != NULL) + { + aTriangulator->Destroy (); + } + + throw Standard_Failure ("BRepMesh_DelabellaBaseMeshAlgo::buildBaseTriangulation: exception in triangulation algorithm"); + } +} diff --git a/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.hxx b/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.hxx new file mode 100644 index 0000000000..43f744fc8e --- /dev/null +++ b/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.hxx @@ -0,0 +1,46 @@ +// Created on: 2019-07-05 +// Copyright (c) 2019 OPEN CASCADE SAS +// Created by: Oleg AGASHIN +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _BRepMesh_DelabellaBaseMeshAlgo_HeaderFile +#define _BRepMesh_DelabellaBaseMeshAlgo_HeaderFile + +#include +#include +#include + +class BRepMesh_DataStructureOfDelaun; +class BRepMesh_Delaun; + +//! Class provides base fuctionality to build face triangulation using Delabella project. +//! Performs generation of mesh using raw data from model. +class BRepMesh_DelabellaBaseMeshAlgo : public BRepMesh_CustomBaseMeshAlgo +{ +public: + + //! Constructor. + Standard_EXPORT BRepMesh_DelabellaBaseMeshAlgo (); + + //! Destructor. + Standard_EXPORT virtual ~BRepMesh_DelabellaBaseMeshAlgo (); + + DEFINE_STANDARD_RTTI_INLINE(BRepMesh_DelabellaBaseMeshAlgo, BRepMesh_CustomBaseMeshAlgo) + +protected: + + //! Builds base triangulation using Delabella project. + Standard_EXPORT virtual void buildBaseTriangulation() Standard_OVERRIDE; +}; + +#endif diff --git a/src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.cxx b/src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.cxx new file mode 100644 index 0000000000..2b218e0f8b --- /dev/null +++ b/src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.cxx @@ -0,0 +1,143 @@ +// Created on: 2019-07-05 +// Copyright (c) 2019 OPEN CASCADE SAS +// Created by: Oleg AGASHIN +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace +{ + struct DefaultBaseMeshAlgo + { + typedef BRepMesh_DelaunayBaseMeshAlgo Type; + }; + + template + struct DefaultNodeInsertionMeshAlgo + { + typedef BRepMesh_DelaunayNodeInsertionMeshAlgo Type; + }; + + struct BaseMeshAlgo + { + typedef BRepMesh_DelabellaBaseMeshAlgo Type; + }; + + template + struct NodeInsertionMeshAlgo + { + typedef BRepMesh_DelaunayNodeInsertionMeshAlgo > Type; + }; + + template + struct DeflectionControlMeshAlgo + { + typedef BRepMesh_DelaunayDeflectionControlMeshAlgo > Type; + }; +} + +//======================================================================= +// Function: Constructor +// Purpose : +//======================================================================= +BRepMesh_DelabellaMeshAlgoFactory::BRepMesh_DelabellaMeshAlgoFactory () +{ +} + +//======================================================================= +// Function: Destructor +// Purpose : +//======================================================================= +BRepMesh_DelabellaMeshAlgoFactory::~BRepMesh_DelabellaMeshAlgoFactory () +{ +} + +//======================================================================= +// Function: GetAlgo +// Purpose : +//======================================================================= +Handle(IMeshTools_MeshAlgo) BRepMesh_DelabellaMeshAlgoFactory::GetAlgo( + const GeomAbs_SurfaceType theSurfaceType, + const IMeshTools_Parameters& theParameters) const +{ + switch (theSurfaceType) + { + case GeomAbs_Plane: + return theParameters.InternalVerticesMode ? + new NodeInsertionMeshAlgo::Type : + new BaseMeshAlgo::Type; + break; + + case GeomAbs_Sphere: + { + NodeInsertionMeshAlgo::Type* aMeshAlgo = + new NodeInsertionMeshAlgo::Type; + aMeshAlgo->SetPreProcessSurfaceNodes (Standard_True); + return aMeshAlgo; + } + break; + + case GeomAbs_Cylinder: + return theParameters.InternalVerticesMode ? + new DefaultNodeInsertionMeshAlgo::Type : + new DefaultBaseMeshAlgo::Type; + break; + + case GeomAbs_Cone: + { + NodeInsertionMeshAlgo::Type* aMeshAlgo = + new NodeInsertionMeshAlgo::Type; + aMeshAlgo->SetPreProcessSurfaceNodes (Standard_True); + return aMeshAlgo; + } + break; + + case GeomAbs_Torus: + { + NodeInsertionMeshAlgo::Type* aMeshAlgo = + new NodeInsertionMeshAlgo::Type; + aMeshAlgo->SetPreProcessSurfaceNodes (Standard_True); + return aMeshAlgo; + } + break; + + case GeomAbs_SurfaceOfRevolution: + { + DeflectionControlMeshAlgo::Type* aMeshAlgo = + new DeflectionControlMeshAlgo::Type; + aMeshAlgo->SetPreProcessSurfaceNodes (Standard_True); + return aMeshAlgo; + } + break; + + default: + { + DeflectionControlMeshAlgo::Type* aMeshAlgo = + new DeflectionControlMeshAlgo::Type; + aMeshAlgo->SetPreProcessSurfaceNodes (Standard_True); + return aMeshAlgo; + } + } +} diff --git a/src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.hxx b/src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.hxx new file mode 100644 index 0000000000..7606b272b7 --- /dev/null +++ b/src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.hxx @@ -0,0 +1,44 @@ +// Created on: 2019-07-05 +// Copyright (c) 2019 OPEN CASCADE SAS +// Created by: Oleg AGASHIN +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _BRepMesh_DelabellaMeshAlgoFactory_HeaderFile +#define _BRepMesh_DelabellaMeshAlgoFactory_HeaderFile + +#include +#include +#include +#include + +//! Implementation of IMeshTools_MeshAlgoFactory providing Delabella-based +//! algorithms of different compexity depending on type of target surface. +class BRepMesh_DelabellaMeshAlgoFactory : public IMeshTools_MeshAlgoFactory +{ +public: + + //! Constructor. + Standard_EXPORT BRepMesh_DelabellaMeshAlgoFactory (); + + //! Destructor. + Standard_EXPORT virtual ~BRepMesh_DelabellaMeshAlgoFactory (); + + //! Creates instance of meshing algorithm for the given type of surface. + Standard_EXPORT virtual Handle(IMeshTools_MeshAlgo) GetAlgo( + const GeomAbs_SurfaceType theSurfaceType, + const IMeshTools_Parameters& theParameters) const Standard_OVERRIDE; + + DEFINE_STANDARD_RTTI_INLINE(BRepMesh_DelabellaMeshAlgoFactory, IMeshTools_MeshAlgoFactory) +}; + +#endif \ No newline at end of file diff --git a/src/BRepMesh/BRepMesh_Delaun.cxx b/src/BRepMesh/BRepMesh_Delaun.cxx index c87bc6ce62..3d02f9d125 100644 --- a/src/BRepMesh/BRepMesh_Delaun.cxx +++ b/src/BRepMesh/BRepMesh_Delaun.cxx @@ -89,7 +89,8 @@ BRepMesh_Delaun::BRepMesh_Delaun ( const Standard_Boolean isFillCircles) : myMeshData ( theOldMesh ), myCircles (new NCollection_IncAllocator( - IMeshData::MEMORY_BLOCK_SIZE_HUGE)) + IMeshData::MEMORY_BLOCK_SIZE_HUGE)), + mySupVert (3) { if (isFillCircles) { @@ -103,7 +104,8 @@ BRepMesh_Delaun::BRepMesh_Delaun ( //======================================================================= BRepMesh_Delaun::BRepMesh_Delaun(IMeshData::Array1OfVertexOfDelaun& theVertices) : myCircles (theVertices.Length(), new NCollection_IncAllocator( - IMeshData::MEMORY_BLOCK_SIZE_HUGE)) + IMeshData::MEMORY_BLOCK_SIZE_HUGE)), + mySupVert (3) { if ( theVertices.Length() > 2 ) { @@ -123,7 +125,8 @@ BRepMesh_Delaun::BRepMesh_Delaun( IMeshData::Array1OfVertexOfDelaun& theVertices) : myMeshData( theOldMesh ), myCircles ( theVertices.Length(), new NCollection_IncAllocator( - IMeshData::MEMORY_BLOCK_SIZE_HUGE)) + IMeshData::MEMORY_BLOCK_SIZE_HUGE)), + mySupVert (3) { if ( theVertices.Length() > 2 ) { @@ -140,7 +143,8 @@ BRepMesh_Delaun::BRepMesh_Delaun( IMeshData::VectorOfInteger& theVertexIndices) : myMeshData( theOldMesh ), myCircles ( theVertexIndices.Length(), new NCollection_IncAllocator( - IMeshData::MEMORY_BLOCK_SIZE_HUGE)) + IMeshData::MEMORY_BLOCK_SIZE_HUGE)), + mySupVert (3) { perform(theVertexIndices); } @@ -155,7 +159,8 @@ BRepMesh_Delaun::BRepMesh_Delaun (const Handle (BRepMesh_DataStructureOfDelaun)& const Standard_Integer theCellsCountV) : myMeshData (theOldMesh), myCircles (theVertexIndices.Length (), new NCollection_IncAllocator( - IMeshData::MEMORY_BLOCK_SIZE_HUGE)) + IMeshData::MEMORY_BLOCK_SIZE_HUGE)), + mySupVert (3) { perform (theVertexIndices, theCellsCountU, theCellsCountV); } @@ -284,14 +289,14 @@ void BRepMesh_Delaun::superMesh(const Bnd_Box2d& theBox) Standard_Real aDeltaMax = Max( aDeltaX, aDeltaY ); Standard_Real aDelta = aDeltaX + aDeltaY; - mySupVert[0] = myMeshData->AddNode( - BRepMesh_Vertex( ( aMinX + aMaxX ) / 2, aMaxY + aDeltaMax, BRepMesh_Free ) ); + mySupVert.Append (myMeshData->AddNode( + BRepMesh_Vertex( ( aMinX + aMaxX ) / 2, aMaxY + aDeltaMax, BRepMesh_Free ) ) ); - mySupVert[1] = myMeshData->AddNode( - BRepMesh_Vertex( aMinX - aDelta, aMinY - aDeltaMin, BRepMesh_Free ) ); + mySupVert.Append (myMeshData->AddNode( + BRepMesh_Vertex( aMinX - aDelta, aMinY - aDeltaMin, BRepMesh_Free ) ) ); - mySupVert[2] = myMeshData->AddNode( - BRepMesh_Vertex( aMaxX + aDelta, aMinY - aDeltaMin, BRepMesh_Free ) ); + mySupVert.Append (myMeshData->AddNode( + BRepMesh_Vertex( aMaxX + aDelta, aMinY - aDeltaMin, BRepMesh_Free ) ) ); Standard_Integer e[3]; Standard_Boolean o[3]; @@ -300,7 +305,7 @@ void BRepMesh_Delaun::superMesh(const Bnd_Box2d& theBox) Standard_Integer aFirstNode = aNodeId; Standard_Integer aLastNode = (aNodeId + 1) % 3; Standard_Integer aLinkIndex = myMeshData->AddLink( BRepMesh_Edge( - mySupVert[aFirstNode], mySupVert[aLastNode], BRepMesh_Free ) ); + mySupVert (aFirstNode), mySupVert (aLastNode), BRepMesh_Free ) ); e[aNodeId] = Abs(aLinkIndex); o[aNodeId] = (aLinkIndex > 0); @@ -367,28 +372,41 @@ void BRepMesh_Delaun::compute(IMeshData::VectorOfInteger& theVertexIndexes) createTrianglesOnNewVertices( theVertexIndexes ); } + RemoveAuxElements (); +} + +//======================================================================= +//function : RemoveAuxElements +//purpose : +//======================================================================= +void BRepMesh_Delaun::RemoveAuxElements () +{ + Handle (NCollection_IncAllocator) aAllocator = new NCollection_IncAllocator ( + IMeshData::MEMORY_BLOCK_SIZE_HUGE); + + IMeshData::MapOfIntegerInteger aLoopEdges (10, aAllocator); + // Destruction of triangles containing a top of the super triangle - BRepMesh_SelectorOfDataStructureOfDelaun aSelector( myMeshData ); - for (Standard_Integer aSupVertId = 0; aSupVertId < 3; ++aSupVertId) - aSelector.NeighboursOfNode( mySupVert[aSupVertId] ); - - aLoopEdges.Clear(); - IMeshData::IteratorOfMapOfInteger aFreeTriangles( aSelector.Elements() ); - for ( ; aFreeTriangles.More(); aFreeTriangles.Next() ) - deleteTriangle( aFreeTriangles.Key(), aLoopEdges ); + BRepMesh_SelectorOfDataStructureOfDelaun aSelector (myMeshData); + for (Standard_Integer aSupVertId = 0; aSupVertId < mySupVert.Size(); ++aSupVertId) + aSelector.NeighboursOfNode (mySupVert (aSupVertId)); + + IMeshData::IteratorOfMapOfInteger aFreeTriangles (aSelector.Elements ()); + for (; aFreeTriangles.More (); aFreeTriangles.Next ()) + deleteTriangle (aFreeTriangles.Key (), aLoopEdges); // All edges that remain free are removed from aLoopEdges; // only the boundary edges of the triangulation remain there - IMeshData::MapOfIntegerInteger::Iterator aFreeEdges( aLoopEdges ); - for ( ; aFreeEdges.More(); aFreeEdges.Next() ) + IMeshData::MapOfIntegerInteger::Iterator aFreeEdges (aLoopEdges); + for (; aFreeEdges.More (); aFreeEdges.Next ()) { - if ( myMeshData->ElementsConnectedTo( aFreeEdges.Key() ).IsEmpty() ) - myMeshData->RemoveLink( aFreeEdges.Key() ); + if (myMeshData->ElementsConnectedTo (aFreeEdges.Key ()).IsEmpty ()) + myMeshData->RemoveLink (aFreeEdges.Key ()); } // The tops of the super triangle are destroyed - for (Standard_Integer aSupVertId = 0; aSupVertId < 3; ++aSupVertId) - myMeshData->RemoveNode( mySupVert[aSupVertId] ); + for (Standard_Integer aSupVertId = 0; aSupVertId < mySupVert.Size (); ++aSupVertId) + myMeshData->RemoveNode (mySupVert (aSupVertId)); } //======================================================================= @@ -771,9 +789,7 @@ void BRepMesh_Delaun::cleanupMesh() myMeshData->ElementNodes (aCurTriangle, v); for (int aNodeIdx = 0; aNodeIdx < 3 && isCanNotBeRemoved; ++aNodeIdx) { - if (v[aNodeIdx] == mySupVert[0] || - v[aNodeIdx] == mySupVert[1] || - v[aNodeIdx] == mySupVert[2]) + if (isSupVertex (v[aNodeIdx])) { isCanNotBeRemoved = Standard_False; } diff --git a/src/BRepMesh/BRepMesh_Delaun.hxx b/src/BRepMesh/BRepMesh_Delaun.hxx index f03a333a11..4d52820428 100755 --- a/src/BRepMesh/BRepMesh_Delaun.hxx +++ b/src/BRepMesh/BRepMesh_Delaun.hxx @@ -147,6 +147,17 @@ public: const Standard_Real theSqTolerance, Standard_Integer& theEdgeOn) const; + //! Explicitly sets ids of auxiliary vertices used to build mesh and used by 3rd-party algorithms. + inline void SetAuxVertices (const IMeshData::VectorOfInteger& theSupVert) + { + mySupVert = theSupVert; + } + + //! Destruction of auxiliary triangles containing the given vertices. + //! Removes auxiliary vertices also. + //! @param theAuxVertices auxiliary vertices to be cleaned up. + Standard_EXPORT void RemoveAuxElements (); + private: enum ReplaceFlag @@ -354,11 +365,25 @@ private: //! Performs insertion of internal edges into mesh. void insertInternalEdges(); + //! Checks whether the given vertex id relates to super contour. + inline Standard_Boolean isSupVertex (const Standard_Integer theVertexIdx) + { + for (IMeshData::VectorOfInteger::Iterator aIt (mySupVert); aIt.More (); aIt.Next ()) + { + if (theVertexIdx == aIt.Value ()) + { + return Standard_True; + } + } + + return Standard_False; + } + private: Handle(BRepMesh_DataStructureOfDelaun) myMeshData; BRepMesh_CircleTool myCircles; - Standard_Integer mySupVert[3]; + IMeshData::VectorOfInteger mySupVert; BRepMesh_Triangle mySupTrian; }; diff --git a/src/BRepMesh/DELABELLA_LICENSE b/src/BRepMesh/DELABELLA_LICENSE new file mode 100644 index 0000000000..5472c6d8a6 --- /dev/null +++ b/src/BRepMesh/DELABELLA_LICENSE @@ -0,0 +1,22 @@ +MIT License + +DELABELLA - Delaunay triangulation library +Copyright (C) 2018 GUMIX - Marcin Sokalski + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/src/BRepMesh/FILES b/src/BRepMesh/FILES index a37a4504e4..d8ad29da70 100755 --- a/src/BRepMesh/FILES +++ b/src/BRepMesh/FILES @@ -86,3 +86,9 @@ BRepMesh_VertexTool.cxx BRepMesh_VertexTool.hxx BRepMesh_CustomBaseMeshAlgo.hxx BRepMesh_CustomDelaunayBaseMeshAlgo.hxx +delabella.h +delabella.cpp +BRepMesh_DelabellaBaseMeshAlgo.hxx +BRepMesh_DelabellaBaseMeshAlgo.cxx +BRepMesh_DelabellaMeshAlgoFactory.hxx +BRepMesh_DelabellaMeshAlgoFactory.cxx diff --git a/src/BRepMesh/delabella.cpp b/src/BRepMesh/delabella.cpp new file mode 100644 index 0000000000..96c7c7517a --- /dev/null +++ b/src/BRepMesh/delabella.cpp @@ -0,0 +1,1021 @@ +/* +DELABELLA - Delaunay triangulation library +Copyright (C) 2018 GUMIX - Marcin Sokalski +*/ + +#include +#include +#include + +#if (defined(__APPLE__)) +#include +#else +#include +#endif + +#include +#include "delabella.h" // we just need LOG() macro + +// assuming BITS is max(X_BITS,Y_BITS) + +typedef double Signed14; // BITS xy coords +typedef double Signed15; // BITS + 1 vect::xy +typedef long double Unsigned28; // 2xBITS z coord +typedef long double Signed29; // 2xBITS + 1 vect::z +typedef long double Signed31; // 2xBITS + 3 norm::z +typedef long double Signed45; // 3xBITS + 3 norm::xy +typedef long double Signed62; // 4xBITS + 6 dot(vect,norm) + +/* +// above typedefs can be used to configure delabella arithmetic types +// in example, EXACT SOLVER (with xy coords limited to 14bits integers in range: +/-8192) +// could be achieved with following configuration: + +typedef int16_t Signed14; // BITS xy coords +typedef int16_t Signed15; // BITS + 1 vect::xy +typedef uint32_t Unsigned28; // 2xBITS z coord +typedef int32_t Signed29; // 2xBITS + 1 vect::z +typedef int32_t Signed31; // 2xBITS + 3 norm::z +typedef int64_t Signed45; // 3xBITS + 3 norm::xy +typedef int64_t Signed62; // 4xBITS + 6 dot(vect,norm) + +// alternatively, one could use some BigInt implementation +// in order to expand xy range +*/ + + +static Unsigned28 s14sqr(const Signed14& s) +{ + return (Unsigned28)((Signed29)s*s); +} + +struct Norm +{ + Signed45 x; + Signed45 y; + Signed31 z; +}; + +struct Vect +{ + Signed15 x, y; + Signed29 z; + + Norm cross (const Vect& v) const // cross prod + { + Norm n; + n.x = (Signed45)y*v.z - (Signed45)z*v.y; + n.y = (Signed45)z*v.x - (Signed45)x*v.z; + n.z = (Signed29)x*v.y - (Signed29)y*v.x; + return n; + } +}; + +struct CDelaBella : IDelaBella +{ + struct Face; + + struct Vert : DelaBella_Vertex + { + Unsigned28 z; + Face* sew; + + Vect operator - (const Vert& v) const // diff + { + Vect d; + d.x = (Signed15)x - (Signed15)v.x; + d.y = (Signed15)y - (Signed15)v.y; + d.z = (Signed29)z - (Signed29)v.z; + return d; + } + + static bool overlap(const Vert* v1, const Vert* v2) + { + return v1->x == v2->x && v1->y == v2->y; + } + + bool operator < (const Vert& v) const + { + return u28cmp(this, &v) < 0; + } + + static int u28cmp(const void* a, const void* b) + { + Vert* va = (Vert*)a; + Vert* vb = (Vert*)b; + if (va->z < vb->z) + return -1; + if (va->z > vb->z) + return 1; + if (va->y < vb->y) + return -1; + if (va->y > vb->y) + return 1; + if (va->x < vb->x) + return -1; + if (va->x > vb->x) + return 1; + if (va->i < vb->i) + return -1; + if (va->i > vb->i) + return 1; + return 0; + } + }; + + struct Face : DelaBella_Triangle + { + Norm n; + + static Face* Alloc(Face** from) + { + Face* f = *from; + *from = (Face*)f->next; + f->next = 0; + return f; + } + + void Free(Face** to) + { + next = *to; + *to = this; + } + + Face* Next(const Vert* p) const + { + // return next face in same direction as face vertices are (cw/ccw) + + if (v[0] == p) + return (Face*)f[1]; + if (v[1] == p) + return (Face*)f[2]; + if (v[2] == p) + return (Face*)f[0]; + return 0; + } + + Signed62 dot(const Vert& p) const // dot + { + Vect d = p - *(Vert*)v[0]; + return (Signed62)n.x * d.x + (Signed62)n.y * d.y + (Signed62)n.z * d.z; + } + + Norm cross() const // cross of diffs + { + return (*(Vert*)v[1] - *(Vert*)v[0]).cross(*(Vert*)v[2] - *(Vert*)v[0]); + } + }; + + Vert* vert_alloc; + Face* face_alloc; + int max_verts; + int max_faces; + + Face* first_dela_face; + Face* first_hull_face; + Vert* first_hull_vert; + + int inp_verts; + int out_verts; + + int(*errlog_proc)(void* file, const char* fmt, ...); + void* errlog_file; + + virtual ~CDelaBella () + { + } + + int Triangulate() + { + int points = inp_verts; + std::sort(vert_alloc, vert_alloc + points); + + // rmove dups + { + int w = 0, r = 1; // skip initial no-dups block + while (r < points && !Vert::overlap(vert_alloc + r, vert_alloc + w)) + { + w++; + r++; + } + + w++; + + while (r < points) + { + r++; + + // skip dups + while (r < points && Vert::overlap(vert_alloc + r, vert_alloc + r - 1)) + r++; + + // copy next no-dups block + while (r < points && !Vert::overlap(vert_alloc + r, vert_alloc + r - 1)) + vert_alloc[w++] = vert_alloc[r++]; + } + + if (points - w) + { + if (errlog_proc) + errlog_proc(errlog_file, "[WRN] detected %d dups in xy array!\n", points - w); + points = w; + } + } + + if (points < 3) + { + if (points == 2) + { + if (errlog_proc) + errlog_proc(errlog_file, "[WRN] all input points are colinear, returning single segment!\n"); + first_hull_vert = vert_alloc + 0; + vert_alloc[0].next = (DelaBella_Vertex*)vert_alloc + 1; + vert_alloc[1].next = 0; + } + else + { + if (errlog_proc) + errlog_proc(errlog_file, "[WRN] all input points are identical, returning signle point!\n"); + first_hull_vert = vert_alloc + 0; + vert_alloc[0].next = 0; + } + + return -points; + } + + int hull_faces = 2 * points - 4; + + if (max_faces < hull_faces) + { + if (max_faces) + free(face_alloc); + max_faces = 0; + face_alloc = (Face*)malloc(sizeof(Face) * hull_faces); + if (face_alloc) + max_faces = hull_faces; + else + { + if (errlog_proc) + errlog_proc(errlog_file, "[ERR] Not enough memory, shop for some more RAM. See you!\n"); + return 0; + } + } + + for (int i = 1; i < hull_faces; i++) + face_alloc[i - 1].next = face_alloc + i; + face_alloc[hull_faces - 1].next = 0; + + Face* cache = face_alloc; + Face* hull = 0; + + Face f; // tmp + f.v[0] = vert_alloc + 0; + f.v[1] = vert_alloc + 1; + f.v[2] = vert_alloc + 2; + f.n = f.cross(); + + bool colinear = f.n.z == 0; + int i = 3; + + ///////////////////////////////////////////////////////////////////////// + // UNTIL INPUT IS COPLANAR, GROW IT IN FORM OF A 2D CONTOUR + /* + . | | after adding . | ________* L + . \ Last points to / Head next point . \ ______/ / + . *____ | -----> .H *____ | + . |\_ \_____ | . |\_ \_____ | + . \ \_ \__* - Tail points to Last . \ \_ \__* T + . \ \_ / . \ \_ / + . \__ \_ __/ . \__ \_ __/ + . \__* - Head points to Tail . \__/ + */ + + Vert* head = (Vert*)f.v[0]; + Vert* tail = (Vert*)f.v[1]; + Vert* last = (Vert*)f.v[2]; + + head->next = tail; + tail->next = last; + last->next = head; + + while (i < points && f.dot(vert_alloc[i]) == 0) + { + Vert* v = vert_alloc + i; + + // it is enough to test just 1 non-zero coord + // but we want also to test stability (assert) + // so we calc all signs... + + // why not testing sign of dot prod of 2 normals? + // that way we'd fall into precission problems + + Norm LvH = (*v - *last).cross(*head - *last); + bool lvh = + (f.n.x > 0 && LvH.x > 0) || (f.n.x < 0 && LvH.x < 0) || + (f.n.y > 0 && LvH.y > 0) || (f.n.y < 0 && LvH.y < 0) || + (f.n.z > 0 && LvH.z > 0) || (f.n.z < 0 && LvH.z < 0); + + Norm TvL = (*v - *tail).cross(*last - *tail); + bool tvl = + (f.n.x > 0 && TvL.x > 0) || (f.n.x < 0 && TvL.x < 0) || + (f.n.y > 0 && TvL.y > 0) || (f.n.y < 0 && TvL.y < 0) || + (f.n.z > 0 && TvL.z > 0) || (f.n.z < 0 && TvL.z < 0); + + if (lvh && !tvl) // insert new f on top of e(2,0) = (last,head) + { + // f.v[0] = head; + f.v[1] = last; + f.v[2] = v; + + last->next = v; + v->next = head; + tail = last; + } + else + if (tvl && !lvh) // insert new f on top of e(1,2) = (tail,last) + { + f.v[0] = last; + //f.v[1] = tail; + f.v[2] = v; + + tail->next = v; + v->next = last; + head = last; + } + else + { + // wtf? dilithium crystals are fucked. + assert(0); + } + + last = v; + i++; + } + + if (i == points) + { + if (colinear) + { + if (errlog_proc) + errlog_proc(errlog_file, "[WRN] all input points are colinear, returning segment list!\n"); + first_hull_vert = head; + last->next = 0; // break contour, make it a list + return -points; + } + else + { + if (points > 3) + { + if (errlog_proc) + errlog_proc(errlog_file, "[NFO] all input points are cocircular.\n"); + } + else + { + if (errlog_proc) + errlog_proc(errlog_file, "[NFO] trivial case of 3 points, thank you.\n"); + + first_dela_face = Face::Alloc(&cache); + first_dela_face->next = 0; + first_hull_face = Face::Alloc(&cache); + first_hull_face->next = 0; + + first_dela_face->f[0] = first_dela_face->f[1] = first_dela_face->f[2] = first_hull_face; + first_hull_face->f[0] = first_hull_face->f[1] = first_hull_face->f[2] = first_dela_face; + + head->sew = tail->sew = last->sew = first_hull_face; + + if (f.n.z < 0) + { + first_dela_face->v[0] = head; + first_dela_face->v[1] = tail; + first_dela_face->v[2] = last; + first_hull_face->v[0] = last; + first_hull_face->v[1] = tail; + first_hull_face->v[2] = head; + + // reverse silhouette + head->next = last; + last->next = tail; + tail->next = head; + + first_hull_vert = last; + } + else + { + first_dela_face->v[0] = last; + first_dela_face->v[1] = tail; + first_dela_face->v[2] = head; + first_hull_face->v[0] = head; + first_hull_face->v[1] = tail; + first_hull_face->v[2] = last; + + first_hull_vert = head; + } + + first_dela_face->n = first_dela_face->cross(); + first_hull_face->n = first_hull_face->cross(); + + return 3; + } + + // retract last point it will be added as a cone's top later + last = head; + head = (Vert*)head->next; + i--; + } + } + + ///////////////////////////////////////////////////////////////////////// + // CREATE CONE HULL WITH TOP AT cloud[i] AND BASE MADE OF CONTOUR LIST + // in 2 ways :( - depending on at which side of the contour a top vertex appears + + Vert* q = vert_alloc + i; + + if (f.dot(*q) > 0) + { + Vert* p = last; + Vert* n = (Vert*)p->next; + + Face* first_side = Face::Alloc(&cache); + first_side->v[0] = p; + first_side->v[1] = n; + first_side->v[2] = q; + first_side->n = first_side->cross(); + hull = first_side; + + p = n; + n = (Vert*)n->next; + + Face* prev_side = first_side; + Face* prev_base = 0; + Face* first_base = 0; + + do + { + Face* base = Face::Alloc(&cache); + base->v[0] = n; + base->v[1] = p; + base->v[2] = last; + base->n = base->cross(); + + Face* side = Face::Alloc(&cache); + side->v[0] = p; + side->v[1] = n; + side->v[2] = q; + side->n = side->cross(); + + side->f[2] = base; + base->f[2] = side; + + side->f[1] = prev_side; + prev_side->f[0] = side; + + base->f[0] = prev_base; + if (prev_base) + prev_base->f[1] = base; + else + first_base = base; + + prev_base = base; + prev_side = side; + + p = n; + n = (Vert*)n->next; + } while (n != last); + + Face* last_side = Face::Alloc(&cache); + last_side->v[0] = p; + last_side->v[1] = n; + last_side->v[2] = q; + last_side->n = last_side->cross(); + + last_side->f[1] = prev_side; + prev_side->f[0] = last_side; + + last_side->f[0] = first_side; + first_side->f[1] = last_side; + + first_base->f[0] = first_side; + first_side->f[2] = first_base; + + last_side->f[2] = prev_base; + prev_base->f[1] = last_side; + + i++; + } + else + { + Vert* p = last; + Vert* n = (Vert*)p->next; + + Face* first_side = Face::Alloc(&cache); + first_side->v[0] = n; + first_side->v[1] = p; + first_side->v[2] = q; + first_side->n = first_side->cross(); + hull = first_side; + + p = n; + n = (Vert*)n->next; + + Face* prev_side = first_side; + Face* prev_base = 0; + Face* first_base = 0; + + do + { + Face* base = Face::Alloc(&cache); + base->v[0] = p; + base->v[1] = n; + base->v[2] = last; + base->n = base->cross(); + + Face* side = Face::Alloc(&cache); + side->v[0] = n; + side->v[1] = p; + side->v[2] = q; + side->n = side->cross(); + + side->f[2] = base; + base->f[2] = side; + + side->f[0] = prev_side; + prev_side->f[1] = side; + + base->f[1] = prev_base; + if (prev_base) + prev_base->f[0] = base; + else + first_base = base; + + prev_base = base; + prev_side = side; + + p = n; + n = (Vert*)n->next; + } while (n != last); + + Face* last_side = Face::Alloc(&cache); + last_side->v[0] = n; + last_side->v[1] = p; + last_side->v[2] = q; + last_side->n = last_side->cross(); + + last_side->f[0] = prev_side; + prev_side->f[1] = last_side; + + last_side->f[1] = first_side; + first_side->f[0] = last_side; + + first_base->f[1] = first_side; + first_side->f[2] = first_base; + + last_side->f[2] = prev_base; + prev_base->f[0] = last_side; + + i++; + } + + ///////////////////////////////////////////////////////////////////////// + // ACTUAL ALGORITHM + + for (; i < points; i++) + { + //ValidateHull(alloc, 2 * i - 4); + + Vert* _q = vert_alloc + i; + Vert* _p = vert_alloc + i - 1; + Face* _f = hull; + + // 1. FIND FIRST VISIBLE FACE + // simply iterate around last vertex using last added triange adjecency info + while (_f->dot(*_q) <= 0) + { + _f = _f->Next(_p); + if (_f == hull) + { + // if no visible face can be located at last vertex, + // let's run through all faces (approximately last to first), + // yes this is emergency fallback and should not ever happen. + _f = face_alloc + 2 * i - 4 - 1; + while (_f->dot(*_q) <= 0) + { + assert(_f != face_alloc); // no face is visible? you must be kidding! + _f--; + } + } + } + + // 2. DELETE VISIBLE FACES & ADD NEW ONES + // (we also build silhouette (vertex loop) between visible & invisible faces) + + int del = 0; + int add = 0; + + // push first visible face onto stack (of visible faces) + Face* stack = _f; + _f->next = _f; // old trick to use list pointers as 'on-stack' markers + while (stack) + { + // pop, take care of last item ptr (it's not null!) + _f = stack; + stack = (Face*)_f->next; + if (stack == _f) + stack = 0; + _f->next = 0; + + // copy parts of old face that we still need after removal + Vert* fv[3] = { (Vert*)_f->v[0],(Vert*)_f->v[1],(Vert*)_f->v[2] }; + Face* ff[3] = { (Face*)_f->f[0],(Face*)_f->f[1],(Face*)_f->f[2] }; + + // delete visible face + _f->Free(&cache); + del++; + + // check all 3 neighbors + for (int e = 0; e < 3; e++) + { + Face* n = ff[e]; + if (n && !n->next) // ensure neighbor is not processed yet & isn't on stack + { + if (n->dot(*_q) <= 0) // if neighbor is not visible we have slihouette edge + { + // build face + add++; + + // ab: given face adjacency [index][], + // it provides [][2] vertex indices on shared edge (CCW order) + const static int ab[3][2] = { { 1,2 },{ 2,0 },{ 0,1 } }; + + Vert* a = fv[ab[e][0]]; + Vert* b = fv[ab[e][1]]; + + Face* s = Face::Alloc(&cache); + s->v[0] = a; + s->v[1] = b; + s->v[2] = _q; + + s->n = s->cross(); + s->f[2] = n; + + // change neighbour's adjacency from old visible face to cone side + if (n->f[0] == _f) + n->f[0] = s; + else + if (n->f[1] == _f) + n->f[1] = s; + else + if (n->f[2] == _f) + n->f[2] = s; + else + assert(0); + + // build silhouette needed for sewing sides in the second pass + a->sew = s; + a->next = b; + } + else + { + // disjoin visible faces + // so they won't be processed more than once + + if (n->f[0] == _f) + n->f[0] = 0; + else + if (n->f[1] == _f) + n->f[1] = 0; + else + if (n->f[2] == _f) + n->f[2] = 0; + else + assert(0); + + // push neighbor face, it's visible and requires processing + n->next = stack ? stack : n; + stack = n; + } + } + } + } + + // if addv[0]; + + Vert* pr = entry; + do + { + // sew pr<->nx + Vert* nx = (Vert*)pr->next; + pr->sew->f[0] = nx->sew; + nx->sew->f[1] = pr->sew; + pr = nx; + } while (pr != entry); + } + + assert(2 * i - 4 == hull_faces); + //ValidateHull(alloc, hull_faces); + + // needed? + for (i = 0; i < points; i++) + { + vert_alloc[i].next = 0; + vert_alloc[i].sew = 0; + } + + i = 0; + Face** prev_dela = &first_dela_face; + Face** prev_hull = &first_hull_face; + for (int j = 0; j < hull_faces; j++) + { + Face* _f = face_alloc + j; + if (_f->n.z < 0) + { + *prev_dela = _f; + prev_dela = (Face**)&_f->next; + i++; + } + else + { + *prev_hull = _f; + prev_hull = (Face**)&_f->next; + if (((Face*)_f->f[0])->n.z < 0) + { + _f->v[1]->next = _f->v[2]; + ((Vert*)_f->v[1])->sew = _f; + } + if (((Face*)_f->f[1])->n.z < 0) + { + _f->v[2]->next = _f->v[0]; + ((Vert*)_f->v[2])->sew = _f; + } + if (((Face*)_f->f[2])->n.z < 0) + { + _f->v[0]->next = _f->v[1]; + ((Vert*)_f->v[0])->sew = _f; + } + } + } + + *prev_dela = 0; + *prev_hull = 0; + + first_hull_vert = (Vert*)first_hull_face->v[0]; + + // todo link slihouette verts into contour + // and sew its edges with hull faces + + return 3*i; + } + + bool ReallocVerts(int points) + { + inp_verts = points; + out_verts = 0; + + first_dela_face = 0; + first_hull_face = 0; + first_hull_vert = 0; + + if (max_verts < points) + { + if (max_verts) + { + free(vert_alloc); + vert_alloc = 0; + max_verts = 0; + } + + vert_alloc = (Vert*)malloc(sizeof(Vert)*points); + if (vert_alloc) + max_verts = points; + else + { + if (errlog_proc) + errlog_proc(errlog_file, "[ERR] Not enough memory, shop for some more RAM. See you!\n"); + return false; + } + } + + return true; + } + + virtual int Triangulate(int points, const float* x, const float* y = 0, int advance_bytes = 0) + { + if (!x) + return 0; + + if (!y) + y = x + 1; + + if (advance_bytes < static_cast(sizeof(float) * 2)) + advance_bytes = static_cast(sizeof(float) * 2); + + if (!ReallocVerts(points)) + return 0; + + for (int i = 0; i < points; i++) + { + Vert* v = vert_alloc + i; + v->i = i; + v->x = (Signed14)*(const float*)((const char*)x + i*advance_bytes); + v->y = (Signed14)*(const float*)((const char*)y + i*advance_bytes); + v->z = s14sqr(v->x) + s14sqr(v->y); + } + + out_verts = Triangulate(); + return out_verts; + } + + virtual int Triangulate(int points, const double* x, const double* y, int advance_bytes) + { + if (!x) + return 0; + + if (!y) + y = x + 1; + + if (advance_bytes < static_cast(sizeof(double) * 2)) + advance_bytes = static_cast(sizeof(double) * 2); + + if (!ReallocVerts(points)) + return 0; + + for (int i = 0; i < points; i++) + { + Vert* v = vert_alloc + i; + v->i = i; + v->x = (Signed14)*(const double*)((const char*)x + i*advance_bytes); + v->y = (Signed14)*(const double*)((const char*)y + i*advance_bytes); + v->z = s14sqr (v->x) + s14sqr (v->y); + } + + out_verts = Triangulate(); + return out_verts; + } + + virtual void Destroy() + { + if (face_alloc) + free(face_alloc); + if (vert_alloc) + free(vert_alloc); + delete this; + } + + // num of points passed to last call to Triangulate() + virtual int GetNumInputPoints() const + { + return inp_verts; + } + + // num of verts returned from last call to Triangulate() + virtual int GetNumOutputVerts() const + { + return out_verts; + } + + virtual const DelaBella_Triangle* GetFirstDelaunayTriangle() const + { + return first_dela_face; + } + + virtual const DelaBella_Triangle* GetFirstHullTriangle() const + { + return first_hull_face; + } + + virtual const DelaBella_Vertex* GetFirstHullVertex() const + { + return first_hull_vert; + } + + virtual void SetErrLog(int(*proc)(void* stream, const char* fmt, ...), void* stream) + { + errlog_proc = proc; + errlog_file = stream; + } +}; + +IDelaBella* IDelaBella::Create() +{ + CDelaBella* db = new CDelaBella; + if (!db) + return 0; + + db->vert_alloc = 0; + db->face_alloc = 0; + db->max_verts = 0; + db->max_faces = 0; + + db->first_dela_face = 0; + db->first_hull_face = 0; + db->first_hull_vert = 0; + + db->inp_verts = 0; + db->out_verts = 0; + + db->errlog_proc = 0; + db->errlog_file = 0; + + return db; +} + +void* DelaBella_Create() +{ + return IDelaBella::Create(); +} + +void DelaBella_Destroy(void* db) +{ + ((IDelaBella*)db)->Destroy(); +} + +void DelaBella_SetErrLog(void* db, int(*proc)(void* stream, const char* fmt, ...), void* stream) +{ + ((IDelaBella*)db)->SetErrLog(proc, stream); +} + +int DelaBella_TriangulateFloat(void* db, int points, float* x, float* y, int advance_bytes) +{ + return ((IDelaBella*)db)->Triangulate(points, x, y, advance_bytes); +} + +int DelaBella_TriangulateDouble(void* db, int points, double* x, double* y, int advance_bytes) +{ + return ((IDelaBella*)db)->Triangulate(points, x, y, advance_bytes); +} + +int DelaBella_GetNumInputPoints(void* db) +{ + return ((IDelaBella*)db)->GetNumInputPoints(); +} + +int DelaBella_GetNumOutputVerts(void* db) +{ + return ((IDelaBella*)db)->GetNumOutputVerts(); +} + +const DelaBella_Triangle* GetFirstDelaunayTriangle(void* db) +{ + return ((IDelaBella*)db)->GetFirstDelaunayTriangle(); +} + +const DelaBella_Triangle* GetFirstHullTriangle(void* db) +{ + return ((IDelaBella*)db)->GetFirstHullTriangle(); +} + +const DelaBella_Vertex* GetFirstHullVertex(void* db) +{ + return ((IDelaBella*)db)->GetFirstHullVertex(); +} + +// depreciated! +int DelaBella(int points, const double* xy, int* abc, int(*errlog)(const char* fmt, ...)) +{ + if (errlog) + errlog("[WRN] Depreciated interface! errlog disabled.\n"); + + if (!xy || points <= 0) + return 0; + + IDelaBella* db = IDelaBella::Create(); + int verts = db->Triangulate(points, xy, 0, 0); + + if (!abc) + return verts; + + if (verts > 0) + { + int tris = verts / 3; + const DelaBella_Triangle* dela = db->GetFirstDelaunayTriangle(); + for (int i = 0; i < tris; i++) + { + for (int j=0; j<3; j++) + abc[3 * i + j] = dela->v[j]->i; + dela = dela->next; + } + } + else + { + int pnts = -verts; + const DelaBella_Vertex* line = db->GetFirstHullVertex(); + for (int i = 0; i < pnts; i++) + { + abc[i] = line->i; + line = line->next; + } + } + + return verts; +} diff --git a/src/BRepMesh/delabella.h b/src/BRepMesh/delabella.h new file mode 100644 index 0000000000..74162f74da --- /dev/null +++ b/src/BRepMesh/delabella.h @@ -0,0 +1,68 @@ +/* +DELABELLA - Delaunay triangulation library +Copyright (C) 2018 GUMIX - Marcin Sokalski +*/ + +#ifndef DELABELLA_H +#define DELABELLA_H + +// returns: positive value: number of triangle indices, negative: number of line segment indices (degenerated input) +// triangle indices in abc array are always returned in clockwise order +// DEPRECIATED. move to new API either extern "C" or IDelaBella (C++) +int DelaBella(int points, const double* xy/*[points][2]*/, int* abc/*[2*points-5][3]*/, int (*errlog)(const char* fmt,...) = printf); + +struct DelaBella_Vertex +{ + int i; // index of original point + double x, y; // coordinates (input copy) + DelaBella_Vertex* next; // next silhouette vertex +}; + +struct DelaBella_Triangle +{ + DelaBella_Vertex* v[3]; // 3 vertices spanning this triangle + DelaBella_Triangle* f[3]; // 3 adjacent faces, f[i] is at the edge opposite to vertex v[i] + DelaBella_Triangle* next; // next triangle (of delaunay set or hull set) +}; + +#ifdef __cplusplus +struct IDelaBella +{ + static IDelaBella* Create(); + + virtual void Destroy() = 0; + + virtual void SetErrLog(int(*proc)(void* stream, const char* fmt, ...), void* stream) = 0; + + // return 0: no output + // negative: all points are colinear, output hull vertices form colinear segment list, no triangles on output + // positive: output hull vertices form counter-clockwise ordered segment contour, delaunay and hull triangles are available + // if 'y' pointer is null, y coords are treated to be located immediately after every x + // if advance_bytes is less than 2*sizeof coordinate type, it is treated as 2*sizeof coordinate type + virtual int Triangulate(int points, const float* x, const float* y = 0, int advance_bytes = 0) = 0; + virtual int Triangulate(int points, const double* x, const double* y = 0, int advance_bytes = 0) = 0; + + // num of points passed to last call to Triangulate() + virtual int GetNumInputPoints() const = 0; + + // num of verts returned from last call to Triangulate() + virtual int GetNumOutputVerts() const = 0; + + virtual const DelaBella_Triangle* GetFirstDelaunayTriangle() const = 0; // valid only if Triangulate() > 0 + virtual const DelaBella_Triangle* GetFirstHullTriangle() const = 0; // valid only if Triangulate() > 0 + virtual const DelaBella_Vertex* GetFirstHullVertex() const = 0; // if Triangulate() < 0 it is list, otherwise closed contour! +}; +#else +void* DelaBella_Create(); +void DelaBella_Destroy(void* db); +void DelaBella_SetErrLog(void* db, int(*proc)(void* stream, const char* fmt, ...), void* stream); +int DelaBella_TriangulateFloat(void* db, int points, float* x, float* y = 0, int advance_bytes = 0); +int DelaBella_TriangulateDouble(void* db, int points, double* x, double* y = 0, int advance_bytes = 0); +int DelaBella_GetNumInputPoints(void* db); +int DelaBella_GetNumOutputVerts(void* db); +const DelaBella_Triangle* GetFirstDelaunayTriangle(void* db); +const DelaBella_Triangle* GetFirstHullTriangle(void* db); +const DelaBella_Vertex* GetFirstHullVertex(void* db); +#endif + +#endif -- 2.39.5