From f2006a6f1932096764e56b88a0c1d85cc4d58b36 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, which can be enabled via IMeshTools_Parameters::MeshAlgo option or CSF_MeshAlgo environment variable. Do not fill cirles filter upon explicit initialization. Call base postProcessMesh functionality after initialization of circles in BRepMesh_CustomDelaunayBaseMeshAlgo. Added Vsprintf() wrapper for vsprintf() preserving C locale. --- dox/overview/overview.md | 4 + .../modeling_algos/modeling_algos.md | 34 + src/BRepMesh/BRepMesh_BaseMeshAlgo.hxx | 4 +- .../BRepMesh_ConstrainedBaseMeshAlgo.hxx | 4 +- src/BRepMesh/BRepMesh_Context.cxx | 44 +- src/BRepMesh/BRepMesh_Context.hxx | 2 +- src/BRepMesh/BRepMesh_CustomBaseMeshAlgo.hxx | 29 +- .../BRepMesh_CustomDelaunayBaseMeshAlgo.hxx | 11 +- .../BRepMesh_DelabellaBaseMeshAlgo.cxx | 193 +++ .../BRepMesh_DelabellaBaseMeshAlgo.hxx | 46 + .../BRepMesh_DelabellaMeshAlgoFactory.cxx | 145 +++ .../BRepMesh_DelabellaMeshAlgoFactory.hxx | 44 + src/BRepMesh/BRepMesh_Delaun.cxx | 82 +- src/BRepMesh/BRepMesh_Delaun.hxx | 29 +- .../BRepMesh_DelaunayBaseMeshAlgo.hxx | 2 +- ...Mesh_DelaunayDeflectionControlMeshAlgo.hxx | 2 +- ...BRepMesh_DelaunayNodeInsertionMeshAlgo.hxx | 2 +- src/BRepMesh/BRepMesh_IncrementalMesh.cxx | 2 +- src/BRepMesh/BRepMesh_MeshAlgoFactory.hxx | 2 +- src/BRepMesh/FILES | 6 + src/BRepMesh/delabella.cpp | 1043 +++++++++++++++++ src/BRepMesh/delabella.pxx | 90 ++ src/IMeshTools/FILES | 1 + src/IMeshTools/IMeshTools_MeshAlgoType.hxx | 25 + src/IMeshTools/IMeshTools_Parameters.hxx | 5 + src/MeshTest/MeshTest.cxx | 74 +- src/Standard/Standard_CString.cxx | 6 + src/Standard/Standard_CString.hxx | 9 + tests/bugs/mesh/bug28089_1 | 24 + tests/bugs/mesh/bug28089_2 | 27 + 30 files changed, 1923 insertions(+), 68 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.cpp create mode 100644 src/BRepMesh/delabella.pxx create mode 100644 src/IMeshTools/IMeshTools_MeshAlgoType.hxx create mode 100644 tests/bugs/mesh/bug28089_1 create mode 100644 tests/bugs/mesh/bug28089_2 diff --git a/dox/overview/overview.md b/dox/overview/overview.md index 6eb6c1b2b9..24d666b656 100644 --- a/dox/overview/overview.md +++ b/dox/overview/overview.md @@ -156,6 +156,10 @@ RapidJSON is optionally used by OCCT for reading glTF files (https://rapidjson.o **DejaVu** fonts are a font family based on the Vera Fonts under a permissive license (MIT-like, https://dejavu-fonts.github.io/License.html). DejaVu Sans (basic Latin sub-set) is used by OCCT as fallback font when no system font is available. +**Delabella** is an open-source, cross-platform implementation of the Newton Apple Wrapper algorithm producing 2D Delaunay triangulation. +Delabella is used by BRepMesh as one of alternative 2D triangulation algorithms. +Delabella is licensed under the MIT license (https://github.com/msokalski/delabella). + Adobe Systems, Inc. provides **Adobe Reader**, which can be used to view files in Portable Document Format (PDF). @section OCCT_OVW_SECTION_3 Documentation diff --git a/dox/user_guides/modeling_algos/modeling_algos.md b/dox/user_guides/modeling_algos/modeling_algos.md index df0cdcd104..9977a562e8 100644 --- a/dox/user_guides/modeling_algos/modeling_algos.md +++ b/dox/user_guides/modeling_algos/modeling_algos.md @@ -3251,6 +3251,40 @@ In general, face meshing algorithms have the following structure: * Classes *BRepMesh_DelaunayNodeInsertionMeshAlgo* and *BRepMesh_SweepLineNodeInsertionMeshAlgo* implement algorithm-specific functionality related to addition of internal nodes supplementing functionality provided by *BRepMesh_NodeInsertionMeshAlgo*; * *BRepMesh_DelaunayDeflectionControlMeshAlgo* extends functionality of *BRepMesh_DelaunayNodeInsertionMeshAlgo* by additional procedure controlling deflection of generated triangles. +BRepMesh provides user a way to switch default triangulation algorithm to a custom one, either implemented by user or available worldwide. +There are three base classes that can be currently used to integrate 3rd-party algorithms: + * *BRepMesh_ConstrainedBaseMeshAlgo* base class for tools providing generation of triangulations with constraints requiring no common processing by BRepMesh; + * *BRepMesh_CustomBaseMeshAlgo* provides the entry point for generic algorithms without support of constraints and supposed for generation of base mesh only. + Constraint edges are processed using standard functionality provided by the component itself upon background mesh produced by 3rd-party solver; + * *BRepMesh_CustomDelaunayBaseMeshAlgo* contains initialization part for tools used by BRepMesh for checks or optimizations using results of 3rd-party algorithm. + +Meshing algorithms could be provided by implemeting IMeshTools_MeshAlgoFactory with related interfaces and passing it to BRepMesh_Context::SetFaceDiscret(). +OCCT comes with two base 2D meshing algorithms: BRepMesh_MeshAlgoFactory (used by default) and BRepMesh_DelabellaMeshAlgoFactory. + +The following example demonstrates how it could be done from Draw environment: +~~~~~ +psphere s 10 + +### Default Algo ### +incmesh s 0.0001 -algo default + +### Delabella Algo ### +incmesh s 0.0001 -algo delabella +~~~~~ + +The code snippet below shows passing a custom mesh factory to BRepMesh_IncrementalMesh: +~~~~~ +IMeshTools_Parameters aMeshParams; +Handle(IMeshTools_Context) aContext = new BRepMesh_Context(); +aContext->SetFaceDiscret (new BRepMesh_FaceDiscret (new BRepMesh_DelabellaMeshAlgoFactory())); + +BRepMesh_IncrementalMesh aMesher; +aMesher.SetShape (aShape); +aMesher.ChangeParameters() = aMeshParams; + +aMesher.Perform (aContext); +~~~~~ + #### Range splitter Range splitter tools provide functionality to generate internal surface nodes defined within the range computed using discrete model data. The base functionality is provided by *BRepMesh_DefaultRangeSplitter* which can be used without modifications in case of planar surface. The default splitter does not generate any internal node. diff --git a/src/BRepMesh/BRepMesh_BaseMeshAlgo.hxx b/src/BRepMesh/BRepMesh_BaseMeshAlgo.hxx index 1ab210fe8e..ddf4419933 100644 --- a/src/BRepMesh/BRepMesh_BaseMeshAlgo.hxx +++ b/src/BRepMesh/BRepMesh_BaseMeshAlgo.hxx @@ -25,7 +25,7 @@ class BRepMesh_DataStructureOfDelaun; class BRepMesh_Delaun; -//! Class provides base fuctionality for algorithms building face triangulation. +//! Class provides base functionality for algorithms building face triangulation. //! Performs initialization of BRepMesh_DataStructureOfDelaun and nodes map structures. class BRepMesh_BaseMeshAlgo : public IMeshTools_MeshAlgo { @@ -108,7 +108,7 @@ protected: private: - //! If the given edge has another pcurve for current face coinsiding with specified one, + //! If the given edge has another pcurve for current face coinciding with specified one, //! returns TopAbs_INTERNAL flag. Elsewhere returns orientation of specified pcurve. TopAbs_Orientation fixSeamEdgeOrientation( const IMeshData::IEdgeHandle& theDEdge, diff --git a/src/BRepMesh/BRepMesh_ConstrainedBaseMeshAlgo.hxx b/src/BRepMesh/BRepMesh_ConstrainedBaseMeshAlgo.hxx index 38a63aab73..59e9836027 100644 --- a/src/BRepMesh/BRepMesh_ConstrainedBaseMeshAlgo.hxx +++ b/src/BRepMesh/BRepMesh_ConstrainedBaseMeshAlgo.hxx @@ -23,7 +23,7 @@ class BRepMesh_DataStructureOfDelaun; class BRepMesh_Delaun; -//! Class provides base fuctionality to build face triangulation using Dealunay approach. +//! Class provides base functionality to build face triangulation using Dealunay approach. //! Performs generation of mesh using raw data from model. class BRepMesh_ConstrainedBaseMeshAlgo : public BRepMesh_BaseMeshAlgo { @@ -49,7 +49,7 @@ protected: return std::pair (-1, -1); } - //! Perfroms processing of generated mesh. + //! Performs processing of generated mesh. //! By default does nothing. //! Expected to be called from method generateMesh() in successor classes. virtual void postProcessMesh (BRepMesh_Delaun& /*theMesher*/, diff --git a/src/BRepMesh/BRepMesh_Context.cxx b/src/BRepMesh/BRepMesh_Context.cxx index e48d038467..fe9b6520d7 100644 --- a/src/BRepMesh/BRepMesh_Context.cxx +++ b/src/BRepMesh/BRepMesh_Context.cxx @@ -20,19 +20,59 @@ #include #include #include + #include +#include +#include +#include //======================================================================= // Function: Constructor // Purpose : //======================================================================= -BRepMesh_Context::BRepMesh_Context () +BRepMesh_Context::BRepMesh_Context (IMeshTools_MeshAlgoType theMeshType) { + if (theMeshType == IMeshTools_MeshAlgoType_DEFAULT) + { + TCollection_AsciiString aValue = OSD_Environment ("CSF_MeshAlgo").Value(); + aValue.LowerCase(); + if (aValue == "watson" + || aValue == "0") + { + theMeshType = IMeshTools_MeshAlgoType_Watson; + } + else if (aValue == "delabella" + || aValue == "1") + { + theMeshType = IMeshTools_MeshAlgoType_Delabella; + } + else + { + if (!aValue.IsEmpty()) + { + Message::SendWarning (TCollection_AsciiString("BRepMesh_Context, ignore unknown algorithm '") + aValue + "' specified in CSF_MeshAlgo variable"); + } + theMeshType = IMeshTools_MeshAlgoType_Watson; + } + } + + Handle (IMeshTools_MeshAlgoFactory) aAlgoFactory; + switch (theMeshType) + { + case IMeshTools_MeshAlgoType_DEFAULT: + case IMeshTools_MeshAlgoType_Watson: + aAlgoFactory = new BRepMesh_MeshAlgoFactory(); + break; + case IMeshTools_MeshAlgoType_Delabella: + aAlgoFactory = new BRepMesh_DelabellaMeshAlgoFactory(); + 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_Context.hxx b/src/BRepMesh/BRepMesh_Context.hxx index a802a6b44b..e2b004d1d6 100644 --- a/src/BRepMesh/BRepMesh_Context.hxx +++ b/src/BRepMesh/BRepMesh_Context.hxx @@ -25,7 +25,7 @@ class BRepMesh_Context : public IMeshTools_Context public: //! Constructor. - Standard_EXPORT BRepMesh_Context (); + Standard_EXPORT BRepMesh_Context (IMeshTools_MeshAlgoType theMeshType = IMeshTools_MeshAlgoType_DEFAULT); //! Destructor. Standard_EXPORT virtual ~BRepMesh_Context (); diff --git a/src/BRepMesh/BRepMesh_CustomBaseMeshAlgo.hxx b/src/BRepMesh/BRepMesh_CustomBaseMeshAlgo.hxx index e175091de9..74d6b0c81f 100644 --- a/src/BRepMesh/BRepMesh_CustomBaseMeshAlgo.hxx +++ b/src/BRepMesh/BRepMesh_CustomBaseMeshAlgo.hxx @@ -25,7 +25,7 @@ class BRepMesh_DataStructureOfDelaun; -//! Class provides base fuctionality to build face triangulation using custom triangulation algorithm. +//! Class provides base functionality to build face triangulation using custom triangulation algorithm. //! Performs generation of mesh using raw data from model. class BRepMesh_CustomBaseMeshAlgo : public BRepMesh_ConstrainedBaseMeshAlgo { @@ -46,19 +46,42 @@ public: protected: //! Generates mesh for the contour stored in data structure. - Standard_EXPORT virtual void generateMesh () Standard_OVERRIDE + Standard_EXPORT virtual void generateMesh (const Message_ProgressRange& theRange) 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 (); - postProcessMesh (aMesher); + postProcessMesh (aMesher, theRange); } protected: diff --git a/src/BRepMesh/BRepMesh_CustomDelaunayBaseMeshAlgo.hxx b/src/BRepMesh/BRepMesh_CustomDelaunayBaseMeshAlgo.hxx index 6bde7b90a8..2143db6359 100644 --- a/src/BRepMesh/BRepMesh_CustomDelaunayBaseMeshAlgo.hxx +++ b/src/BRepMesh/BRepMesh_CustomDelaunayBaseMeshAlgo.hxx @@ -19,7 +19,7 @@ class BRepMesh_DataStructureOfDelaun; class BRepMesh_Delaun; -//! Class provides base fuctionality to build face triangulation using custom +//! Class provides base functionality to build face triangulation using custom //! triangulation algorithm with possibility to modify final mesh. //! Performs generation of mesh using raw data from model. template @@ -39,14 +39,15 @@ public: protected: - //! Perfroms processing of generated mesh. - virtual void postProcessMesh(BRepMesh_Delaun& theMesher) + //! Performs processing of generated mesh. + virtual void postProcessMesh (BRepMesh_Delaun& theMesher, + const Message_ProgressRange& theRange) { - BaseAlgo::postProcessMesh (theMesher); - const Handle(BRepMesh_DataStructureOfDelaun)& aStructure = this->getStructure(); std::pair aCellsCount = this->getCellsCount (aStructure->NbNodes()); theMesher.InitCirclesTool (aCellsCount.first, aCellsCount.second); + + BaseAlgo::postProcessMesh (theMesher, theRange); } }; diff --git a/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.cxx b/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.cxx new file mode 100644 index 0000000000..40bbb361e9 --- /dev/null +++ b/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.cxx @@ -0,0 +1,193 @@ +// 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 "delabella.pxx" + +IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_DelabellaBaseMeshAlgo, BRepMesh_CustomBaseMeshAlgo) + +namespace +{ + //! Redirect algorithm messages to OCCT messenger. + static int logDelabella2Occ (void* theStream, const char* theFormat, ...) + { + (void )theStream; + char aBuffer[1024]; // should be more than enough for Delabella messages + + va_list anArgList; + va_start(anArgList, theFormat); + Vsprintf(aBuffer, theFormat, anArgList); + va_end(anArgList); + + Message_Gravity aGravity = Message_Warning; + switch ((int )theFormat[1]) + { + case int('E'): aGravity = Message_Fail; break; // [ERR] + case int('W'): aGravity = Message_Trace; break; // [WRN] + case int('N'): aGravity = Message_Trace; break; // [NFO] + } + Message::Send (aBuffer, aGravity); + return 0; + } +} + +//======================================================================= +// 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(); + if (aTriangulator == NULL) // should never happen + { + throw Standard_ProgramError ("BRepMesh_DelabellaBaseMeshAlgo::buildBaseTriangulation: unable creating a triangulation algorithm"); + } + + aTriangulator->SetErrLog (logDelabella2Occ, NULL); + try + { + 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 (); + aTriangulator = NULL; + } + catch (Standard_Failure const& theException) + { + if (aTriangulator != NULL) + { + aTriangulator->Destroy (); + aTriangulator = NULL; + } + + throw Standard_Failure (theException); + } + catch (...) + { + if (aTriangulator != NULL) + { + aTriangulator->Destroy (); + aTriangulator = NULL; + } + + 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..9d6f2e156a --- /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 functionality 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_RTTIEXT(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..8b02c4a46c --- /dev/null +++ b/src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.cxx @@ -0,0 +1,145 @@ +// 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; + }; +} + +IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_DelabellaMeshAlgoFactory, IMeshTools_MeshAlgoFactory) + +//======================================================================= +// 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..486e0d51ae --- /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 complexity 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_RTTIEXT(BRepMesh_DelabellaMeshAlgoFactory, IMeshTools_MeshAlgoFactory) +}; + +#endif diff --git a/src/BRepMesh/BRepMesh_Delaun.cxx b/src/BRepMesh/BRepMesh_Delaun.cxx index 5a8d30ffe5..6ca5cee049 100644 --- a/src/BRepMesh/BRepMesh_Delaun.cxx +++ b/src/BRepMesh/BRepMesh_Delaun.cxx @@ -90,9 +90,10 @@ 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), + myInitCircles (Standard_False) { - memset (mySupVert, 0, sizeof (mySupVert)); if (isFillCircles) { InitCirclesTool (theCellsCountU, theCellsCountV); @@ -105,9 +106,10 @@ 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), + myInitCircles (Standard_False) { - memset (mySupVert, 0, sizeof (mySupVert)); if ( theVertices.Length() > 2 ) { myMeshData = new BRepMesh_DataStructureOfDelaun( @@ -126,9 +128,10 @@ 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), + myInitCircles (Standard_False) { - memset (mySupVert, 0, sizeof (mySupVert)); if ( theVertices.Length() > 2 ) { Init( theVertices ); @@ -144,9 +147,10 @@ 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), + myInitCircles (Standard_False) { - memset (mySupVert, 0, sizeof (mySupVert)); perform(theVertexIndices); } @@ -160,9 +164,10 @@ 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), + myInitCircles (Standard_False) { - memset (mySupVert, 0, sizeof (mySupVert)); perform (theVertexIndices, theCellsCountU, theCellsCountV); } @@ -240,6 +245,8 @@ void BRepMesh_Delaun::initCirclesTool (const Bnd_Box2d& theBox, myCircles.SetMinMaxSize( gp_XY( aMinX, aMinY ), gp_XY( aMaxX, aMaxY ) ); myCircles.SetCellSize ( aDeltaX / Max (theCellsCountU, aScaler), aDeltaY / Max (theCellsCountV, aScaler)); + + myInitCircles = Standard_True; } //======================================================================= @@ -290,14 +297,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]; @@ -324,7 +331,7 @@ void BRepMesh_Delaun::superMesh(const Bnd_Box2d& theBox) void BRepMesh_Delaun::deleteTriangle(const Standard_Integer theIndex, IMeshData::MapOfIntegerInteger& theLoopEdges ) { - if (!myCircles.IsEmpty()) + if (myInitCircles) { myCircles.Delete (theIndex); } @@ -373,12 +380,25 @@ void BRepMesh_Delaun::compute(IMeshData::VectorOfInteger& theVertexIndexes) createTrianglesOnNewVertices (theVertexIndexes, Message_ProgressRange()); } + 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) + BRepMesh_SelectorOfDataStructureOfDelaun aSelector (myMeshData); + for (Standard_Integer aSupVertId = 0; aSupVertId < mySupVert.Size(); ++aSupVertId) aSelector.NeighboursOfNode( mySupVert[aSupVertId] ); - - aLoopEdges.Clear(); + IMeshData::IteratorOfMapOfInteger aFreeTriangles( aSelector.Elements() ); for ( ; aFreeTriangles.More(); aFreeTriangles.Next() ) deleteTriangle( aFreeTriangles.Key(), aLoopEdges ); @@ -393,7 +413,7 @@ void BRepMesh_Delaun::compute(IMeshData::VectorOfInteger& theVertexIndexes) } // The tops of the super triangle are destroyed - for (Standard_Integer aSupVertId = 0; aSupVertId < 3; ++aSupVertId) + for (Standard_Integer aSupVertId = 0; aSupVertId < mySupVert.Size (); ++aSupVertId) myMeshData->RemoveNode( mySupVert[aSupVertId] ); } @@ -786,9 +806,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; } @@ -1246,11 +1264,15 @@ void BRepMesh_Delaun::addTriangle( const Standard_Integer (&theEdgesId)[3], myMeshData->AddElement(BRepMesh_Triangle(theEdgesId, theEdgesOri, BRepMesh_Free)); - Standard_Boolean isAdded = myCircles.Bind( - aNewTriangleId, - GetVertex( theNodesId[0] ).Coord(), - GetVertex( theNodesId[1] ).Coord(), - GetVertex( theNodesId[2] ).Coord() ); + Standard_Boolean isAdded = Standard_True; + if (myInitCircles) + { + isAdded = myCircles.Bind( + aNewTriangleId, + GetVertex( theNodesId[0] ).Coord(), + GetVertex( theNodesId[1] ).Coord(), + GetVertex( theNodesId[2] ).Coord() ); + } if ( !isAdded ) myMeshData->RemoveElement( aNewTriangleId ); diff --git a/src/BRepMesh/BRepMesh_Delaun.hxx b/src/BRepMesh/BRepMesh_Delaun.hxx index d60b76dc18..f0812d48d1 100755 --- a/src/BRepMesh/BRepMesh_Delaun.hxx +++ b/src/BRepMesh/BRepMesh_Delaun.hxx @@ -149,6 +149,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 @@ -353,13 +364,27 @@ private: //! Performs insertion of internal edges into mesh. void insertInternalEdges(); + //! Checks whether the given vertex id relates to super contour. + Standard_Boolean isSupVertex (const Standard_Integer theVertexIdx) const + { + 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; + Standard_Boolean myInitCircles; BRepMesh_Triangle mySupTrian; - }; #endif diff --git a/src/BRepMesh/BRepMesh_DelaunayBaseMeshAlgo.hxx b/src/BRepMesh/BRepMesh_DelaunayBaseMeshAlgo.hxx index 60793d7f4e..627d489cf0 100644 --- a/src/BRepMesh/BRepMesh_DelaunayBaseMeshAlgo.hxx +++ b/src/BRepMesh/BRepMesh_DelaunayBaseMeshAlgo.hxx @@ -23,7 +23,7 @@ class BRepMesh_DataStructureOfDelaun; class BRepMesh_Delaun; -//! Class provides base fuctionality to build face triangulation using Dealunay approach. +//! Class provides base functionality to build face triangulation using Dealunay approach. //! Performs generation of mesh using raw data from model. class BRepMesh_DelaunayBaseMeshAlgo : public BRepMesh_ConstrainedBaseMeshAlgo { diff --git a/src/BRepMesh/BRepMesh_DelaunayDeflectionControlMeshAlgo.hxx b/src/BRepMesh/BRepMesh_DelaunayDeflectionControlMeshAlgo.hxx index 128901f558..ea355e797f 100644 --- a/src/BRepMesh/BRepMesh_DelaunayDeflectionControlMeshAlgo.hxx +++ b/src/BRepMesh/BRepMesh_DelaunayDeflectionControlMeshAlgo.hxx @@ -46,7 +46,7 @@ public: protected: - //! Perfroms processing of generated mesh. Generates surface nodes and inserts them into structure. + //! Performs processing of generated mesh. Generates surface nodes and inserts them into structure. virtual void postProcessMesh (BRepMesh_Delaun& theMesher, const Message_ProgressRange& theRange) Standard_OVERRIDE { diff --git a/src/BRepMesh/BRepMesh_DelaunayNodeInsertionMeshAlgo.hxx b/src/BRepMesh/BRepMesh_DelaunayNodeInsertionMeshAlgo.hxx index d9f46e9c46..a2e6d5037e 100644 --- a/src/BRepMesh/BRepMesh_DelaunayNodeInsertionMeshAlgo.hxx +++ b/src/BRepMesh/BRepMesh_DelaunayNodeInsertionMeshAlgo.hxx @@ -84,7 +84,7 @@ protected: &this->getRangeSplitter()); } - //! Perfroms processing of generated mesh. Generates surface nodes and inserts them into structure. + //! Performs processing of generated mesh. Generates surface nodes and inserts them into structure. virtual void postProcessMesh (BRepMesh_Delaun& theMesher, const Message_ProgressRange& theRange) Standard_OVERRIDE { diff --git a/src/BRepMesh/BRepMesh_IncrementalMesh.cxx b/src/BRepMesh/BRepMesh_IncrementalMesh.cxx index cfae699e78..eb00e303eb 100644 --- a/src/BRepMesh/BRepMesh_IncrementalMesh.cxx +++ b/src/BRepMesh/BRepMesh_IncrementalMesh.cxx @@ -88,7 +88,7 @@ BRepMesh_IncrementalMesh::~BRepMesh_IncrementalMesh() //======================================================================= void BRepMesh_IncrementalMesh::Perform(const Message_ProgressRange& theRange) { - Handle(BRepMesh_Context) aContext = new BRepMesh_Context; + Handle(BRepMesh_Context) aContext = new BRepMesh_Context (myParameters.MeshAlgo); Perform (aContext, theRange); } diff --git a/src/BRepMesh/BRepMesh_MeshAlgoFactory.hxx b/src/BRepMesh/BRepMesh_MeshAlgoFactory.hxx index 72a27c5587..ed6b9b1db5 100644 --- a/src/BRepMesh/BRepMesh_MeshAlgoFactory.hxx +++ b/src/BRepMesh/BRepMesh_MeshAlgoFactory.hxx @@ -22,7 +22,7 @@ #include //! Default implementation of IMeshTools_MeshAlgoFactory providing algorithms -//! of different compexity depending on type of target surface. +//! of different complexity depending on type of target surface. class BRepMesh_MeshAlgoFactory : public IMeshTools_MeshAlgoFactory { public: diff --git a/src/BRepMesh/FILES b/src/BRepMesh/FILES index a37a4504e4..44333fcb2f 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.pxx +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..e88ccf297f --- /dev/null +++ b/src/BRepMesh/delabella.cpp @@ -0,0 +1,1043 @@ +/* + +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. + +*/ + +#include +#include +#include + +#if (defined(__APPLE__)) +#include +#else +#include +#endif + +#include +#include "delabella.pxx" // 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.pxx b/src/BRepMesh/delabella.pxx new file mode 100644 index 0000000000..a7ff00b56a --- /dev/null +++ b/src/BRepMesh/delabella.pxx @@ -0,0 +1,90 @@ +/* + +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. + +*/ + +#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 diff --git a/src/IMeshTools/FILES b/src/IMeshTools/FILES index e582bc1d1b..5c4537202c 100644 --- a/src/IMeshTools/FILES +++ b/src/IMeshTools/FILES @@ -2,6 +2,7 @@ IMeshTools_Context.hxx IMeshTools_CurveTessellator.hxx IMeshTools_MeshAlgo.hxx IMeshTools_MeshAlgoFactory.hxx +IMeshTools_MeshAlgoType.hxx IMeshTools_MeshBuilder.hxx IMeshTools_MeshBuilder.cxx IMeshTools_ModelAlgo.hxx diff --git a/src/IMeshTools/IMeshTools_MeshAlgoType.hxx b/src/IMeshTools/IMeshTools_MeshAlgoType.hxx new file mode 100644 index 0000000000..ca71c6c873 --- /dev/null +++ b/src/IMeshTools/IMeshTools_MeshAlgoType.hxx @@ -0,0 +1,25 @@ +// Copyright (c) 2020 OPEN CASCADE SAS +// +// 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 _IMeshTools_MeshAlgoType_HeaderFile +#define _IMeshTools_MeshAlgoType_HeaderFile + +//! Enumerates built-in meshing algorithms factories implementing IMeshTools_MeshAlgoFactory interface. +enum IMeshTools_MeshAlgoType +{ + IMeshTools_MeshAlgoType_DEFAULT = -1, //!< use global default (IMeshTools_MeshAlgoType_Watson or CSF_MeshAlgo) + IMeshTools_MeshAlgoType_Watson = 0, //!< generate 2D Delaunay triangulation based on Watson algorithm (BRepMesh_MeshAlgoFactory) + IMeshTools_MeshAlgoType_Delabella, //!< generate 2D Delaunay triangulation based on Delabella algorithm (BRepMesh_DelabellaMeshAlgoFactory) +}; + +#endif diff --git a/src/IMeshTools/IMeshTools_Parameters.hxx b/src/IMeshTools/IMeshTools_Parameters.hxx index 17c3591113..4bbff6a902 100644 --- a/src/IMeshTools/IMeshTools_Parameters.hxx +++ b/src/IMeshTools/IMeshTools_Parameters.hxx @@ -16,6 +16,7 @@ #ifndef _IMeshTools_Parameters_HeaderFile #define _IMeshTools_Parameters_HeaderFile +#include #include //! Structure storing meshing parameters @@ -24,6 +25,7 @@ struct IMeshTools_Parameters { //! Default constructor IMeshTools_Parameters () : + MeshAlgo (IMeshTools_MeshAlgoType_DEFAULT), Angle(0.5), Deflection(0.001), AngleInterior(-1.0), @@ -47,6 +49,9 @@ struct IMeshTools_Parameters { return 0.1; } + //! 2D Delaunay triangulation algorithm factory to use + IMeshTools_MeshAlgoType MeshAlgo; + //! Angular deflection used to tessellate the boundary edges Standard_Real Angle; diff --git a/src/MeshTest/MeshTest.cxx b/src/MeshTest/MeshTest.cxx index 5d337dee44..d02cfcd917 100644 --- a/src/MeshTest/MeshTest.cxx +++ b/src/MeshTest/MeshTest.cxx @@ -43,6 +43,11 @@ #include #include +#include +#include +#include +#include + //epa Memory leaks test //OAN: for triepoints #ifdef _WIN32 @@ -94,7 +99,8 @@ options:\n\ -adjust_min enables local adjustment of min size depending on edge size (switched off by default)\n\ -force_face_def disables usage of shape tolerances for computing face deflection (switched off by default)\n\ -decrease enforces the meshing of the shape even if current mesh satisfies the new criteria\ - (switched off by default).\n"; + (switched off by default).\n\ + -algo {watson|delabella} changes core triangulation algorithm to one with specified id (watson is used by default)\n"; return 0; } @@ -109,6 +115,7 @@ options:\n\ aMeshParams.Deflection = aMeshParams.DeflectionInterior = Max(Draw::Atof(argv[2]), Precision::Confusion()); + Handle (IMeshTools_Context) aContext = new BRepMesh_Context; if (nbarg > 3) { Standard_Integer i = 3; @@ -135,23 +142,54 @@ options:\n\ aMeshParams.AllowQualityDecrease = Standard_True; else if (i < nbarg) { - Standard_Real aVal = Draw::Atof(argv[i++]); - if (aOpt == "-a") + if (aOpt == "-algo") { - aMeshParams.Angle = aVal * M_PI / 180.; - } - else if (aOpt == "-ai") - { - aMeshParams.AngleInterior = aVal * M_PI / 180.; - } - else if (aOpt == "-min") - aMeshParams.MinSize = aVal; - else if (aOpt == "-di") - { - aMeshParams.DeflectionInterior = aVal; + TCollection_AsciiString anAlgoStr (argv[i++]); + anAlgoStr.LowerCase(); + if (anAlgoStr == "watson" + || anAlgoStr == "0") + { + aMeshParams.MeshAlgo = IMeshTools_MeshAlgoType_Watson; + aContext->SetFaceDiscret (new BRepMesh_FaceDiscret (new BRepMesh_MeshAlgoFactory)); + } + else if (anAlgoStr == "delabella" + || anAlgoStr == "1") + { + aMeshParams.MeshAlgo = IMeshTools_MeshAlgoType_Delabella; + aContext->SetFaceDiscret (new BRepMesh_FaceDiscret (new BRepMesh_DelabellaMeshAlgoFactory)); + } + else if (anAlgoStr == "-1" + || anAlgoStr == "default") + { + // already handled by BRepMesh_Context constructor + //aMeshParams.MeshAlgo = IMeshTools_MeshAlgoType_DEFAULT; + } + else + { + di << "Syntax error at " << anAlgoStr; + return 1; + } } else - --i; + { + Standard_Real aVal = Draw::Atof(argv[i++]); + if (aOpt == "-a") + { + aMeshParams.Angle = aVal * M_PI / 180.; + } + else if (aOpt == "-ai") + { + aMeshParams.AngleInterior = aVal * M_PI / 180.; + } + else if (aOpt == "-min") + aMeshParams.MinSize = aVal; + else if (aOpt == "-di") + { + aMeshParams.DeflectionInterior = aVal; + } + else + --i; + } } } } @@ -160,7 +198,11 @@ options:\n\ << (aMeshParams.InParallel ? "ON" : "OFF") << "\n"; Handle(Draw_ProgressIndicator) aProgress = new Draw_ProgressIndicator(di, 1); - BRepMesh_IncrementalMesh aMesher (aShape, aMeshParams, aProgress->Start()); + BRepMesh_IncrementalMesh aMesher; + aMesher.SetShape (aShape); + aMesher.ChangeParameters() = aMeshParams; + + aMesher.Perform (aContext, aProgress->Start()); di << "Meshing statuses: "; const Standard_Integer aStatus = aMesher.GetStatusFlags(); diff --git a/src/Standard/Standard_CString.cxx b/src/Standard/Standard_CString.cxx index 9e70d38166..b9a40d69d2 100755 --- a/src/Standard/Standard_CString.cxx +++ b/src/Standard/Standard_CString.cxx @@ -133,3 +133,9 @@ int Sprintf (char* theBuffer, const char* theFormat, ...) va_end(argp); return result; } + +int Vsprintf (char* theBuffer, const char* theFormat, va_list theArgList) +{ + SAVE_TL(); + return vsprintf_l(theBuffer, Standard_CLocaleSentry::GetCLocale(), theFormat, theArgList); +} diff --git a/src/Standard/Standard_CString.hxx b/src/Standard/Standard_CString.hxx index 319a63d22e..451fc2f534 100644 --- a/src/Standard/Standard_CString.hxx +++ b/src/Standard/Standard_CString.hxx @@ -85,6 +85,15 @@ Standard_EXPORT int Fprintf (FILE* theFile, const char* theFormat, ...); //! Equivalent of standard C function sprintf() that always uses C locale Standard_EXPORT int Sprintf (char* theBuffer, const char* theFormat, ...); +//! Equivalent of standard C function vsprintf() that always uses C locale. +//! Note that this function does not check buffer bounds and should be used with precaution measures +//! (only with format fitting into the buffer of known size). +//! @param theBuffer [in] [out] string buffer to fill +//! @param theFormat [in] format to apply +//! @param theArgList [in] argument list for specified format +//! @return the total number of characters written, or a negative number on error +Standard_EXPORT int Vsprintf (char* theBuffer, const char* theFormat, va_list theArgList); + #ifdef __cplusplus } #endif /* __cplusplus */ diff --git a/tests/bugs/mesh/bug28089_1 b/tests/bugs/mesh/bug28089_1 new file mode 100644 index 0000000000..7e731c05b1 --- /dev/null +++ b/tests/bugs/mesh/bug28089_1 @@ -0,0 +1,24 @@ +puts "========" +puts "0028089: Mesh - New algorithm for triangulation of 2d polygons" +puts "========" +puts "" + +pload MODELING VISUALIZATION +restore [locate_data_file terrain.brep] t1 +tclean t1 +#tscale t1 0 0 0 0.01 +tcopy t1 t2 +ttranslate t1 0 0 0 +ttranslate t2 8000 0 0 +incmesh t1 1000 -algo 0 -a 50 +trinfo t1 +incmesh t2 1000 -algo 1 -a 50 +trinfo t2 +vclear +vinit View1 +vtop +vdefaults -autotriang 0 +vdisplay -dispMode 1 t1 t2 +vfit +vaspects t1 -drawEdges 1 +vaspects t2 -drawEdges 1 diff --git a/tests/bugs/mesh/bug28089_2 b/tests/bugs/mesh/bug28089_2 new file mode 100644 index 0000000000..02e3d1324f --- /dev/null +++ b/tests/bugs/mesh/bug28089_2 @@ -0,0 +1,27 @@ +puts "========" +puts "0028089: Mesh - New algorithm for triangulation of 2d polygons" +puts "========" +puts "" + +pload MODELING VISUALIZATION +psphere s1 1 +psphere s2 1 +box b1 0 2 0 1 2 2 +box b2 2 2 0 1 2 2 +ttranslate s2 2 0 0 +incmesh s1 0.5 -algo 0 -a 50 +incmesh s2 0.5 -algo 1 -a 50 +incmesh b1 0.5 -algo 0 +incmesh b2 0.5 -algo 1 + +vclear +vinit View1 +vdefaults -autotriang 0 +vdisplay -dispMode 1 s1 s2 b1 b2 +vfit +vaspects s1 -drawEdges 1 +vaspects s2 -drawEdges 1 +vaspects b1 -drawEdges 1 +vaspects b2 -drawEdges 1 + +vdump ${imagedir}/${casename}.png