EgtGeomKernel :

- correzione errori compliazione con Clang-cl/LLVM.

EgtGeomKernel.vcxproj

@@ -114,6 +114,7 @@
       <MultiProcessorCompilation>true</MultiProcessorCompilation>
       <DebugInformationFormat>ProgramDatabase</DebugInformationFormat>
       <EnablePREfast>false</EnablePREfast>
+      <AdditionalIncludeDirectories>\EgtDev\Extern\GeometricToolsEngine</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Windows</SubSystem>
@@ -148,7 +149,7 @@ copy $(TargetPath) \EgtProg\DllD32</Command>
       <MinimalRebuild>false</MinimalRebuild>
       <MultiProcessorCompilation>true</MultiProcessorCompilation>
       <DebugInformationFormat>ProgramDatabase</DebugInformationFormat>
-      <LanguageStandard>stdcpp17</LanguageStandard>
+      <AdditionalIncludeDirectories>\EgtDev\Extern\GeometricToolsEngine</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Windows</SubSystem>
@@ -195,6 +196,7 @@ copy $(TargetPath) \EgtProg\DllD64</Command>
       <EnableParallelCodeGeneration>true</EnableParallelCodeGeneration>
       <WholeProgramOptimization>false</WholeProgramOptimization>
       <DebugInformationFormat>ProgramDatabase</DebugInformationFormat>
+      <AdditionalIncludeDirectories>\EgtDev\Extern\GeometricToolsEngine</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Windows</SubSystem>
@@ -240,7 +242,7 @@ copy $(TargetPath) \EgtProg\Dll32</Command>
       <EnableFiberSafeOptimizations>true</EnableFiberSafeOptimizations>
       <WholeProgramOptimization>false</WholeProgramOptimization>
       <DebugInformationFormat>None</DebugInformationFormat>
-      <LanguageStandard>stdcpp17</LanguageStandard>
+      <AdditionalIncludeDirectories>\EgtDev\Extern\GeometricToolsEngine</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Windows</SubSystem>
@@ -397,6 +399,8 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
     <ClCompile Include="SurfTriMeshBooleans.cpp" />
     <ClCompile Include="TextureData.cpp" />
     <ClCompile Include="Tool.cpp" />
+    <ClCompile Include="tpp_assert.cpp" />
+    <ClCompile Include="tpp_impl.cpp" />
     <ClCompile Include="UserObjDefault.cpp" />
     <ClCompile Include="UserObjFactory.cpp" />
     <ClCompile Include="OutTsc.cpp" />
@@ -560,6 +564,7 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
     <ClInclude Include="DistPointCrvComposite.h" />
     <ClInclude Include="DistPointLine.h" />
     <ClInclude Include="DllMain.h" />
+    <ClInclude Include="dpoint.hpp" />
     <ClInclude Include="ExtDimension.h" />
     <ClInclude Include="ExtText.h" />
     <ClInclude Include="FontAux.h" />
@@ -607,6 +612,10 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
     <ClInclude Include="TextureData.h" />
     <ClInclude Include="Tool.h" />
     <ClInclude Include="CAvToolTriangle.h" />
+    <ClInclude Include="tpp_assert.hpp" />
+    <ClInclude Include="tpp_interface.hpp" />
+    <ClInclude Include="triangle.h" />
+    <ClInclude Include="triangle_impl.hpp" />
     <ClInclude Include="UserObjDefault.h" />
     <ClInclude Include="OutTsc.h" />
     <ClInclude Include="PointsPCA.h" />

EgtGeomKernel.vcxproj.filters

@@ -46,6 +46,12 @@
     <Filter Include="File di origine\GeoCollision">
       <UniqueIdentifier>{865b76ee-b10d-41fc-861c-b48ce52fa277}</UniqueIdentifier>
     </Filter>
+    <Filter Include="File di intestazione\tpp">
+      <UniqueIdentifier>{9596e38a-8880-4ab5-bd65-e1db39c92ac5}</UniqueIdentifier>
+    </Filter>
+    <Filter Include="File di origine\tpp">
+      <UniqueIdentifier>{636e60fa-4b24-4a7b-b2e9-74e1c1e8fb31}</UniqueIdentifier>
+    </Filter>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="Vector3d.cpp">
@@ -453,6 +459,12 @@
     <ClCompile Include="CDeRectPrismoidTria.cpp">
       <Filter>File di origine\GeoCollision</Filter>
     </ClCompile>
+    <ClCompile Include="tpp_assert.cpp">
+      <Filter>File di origine\tpp</Filter>
+    </ClCompile>
+    <ClCompile Include="tpp_impl.cpp">
+      <Filter>File di origine\tpp</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="stdafx.h">
@@ -1076,6 +1088,21 @@
     <ClInclude Include="CDeRectPrismoidTria.h">
       <Filter>File di intestazione</Filter>
     </ClInclude>
+    <ClInclude Include="dpoint.hpp">
+      <Filter>File di intestazione\tpp</Filter>
+    </ClInclude>
+    <ClInclude Include="tpp_assert.hpp">
+      <Filter>File di intestazione\tpp</Filter>
+    </ClInclude>
+    <ClInclude Include="tpp_interface.hpp">
+      <Filter>File di intestazione\tpp</Filter>
+    </ClInclude>
+    <ClInclude Include="triangle.h">
+      <Filter>File di intestazione\tpp</Filter>
+    </ClInclude>
+    <ClInclude Include="triangle_impl.hpp">
+      <Filter>File di intestazione\tpp</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="EgtGeomKernel.rc">

PolyLine.cpp

@@ -22,6 +22,8 @@
 #include "/EgtDev/Include/EGkPolyLine.h"
 #include "/EgtDev/Include/EGkPlane3d.h"
 #include "/EgtDev/Include/EGnStringUtils.h"
+#include "/EgtDev/Include/EgtNumUtils.h"
+#include "/EgtDev/Include/EGkPolygon3d.h"
 
 using namespace std ;
 
@@ -1296,3 +1298,326 @@ PolyLine::Trim( const Plane3d& plPlane, bool bInVsOut)
 
    return true ;
 }
+
+
+//----------------------------------------------------------------------------
+/*static*/ bool
+ChangePolyLineStart( const Point3d& ptNewStart, PolyLine& Loop, double dTol)
+{
+   // Rinomino la lista di punti della PolyLine.
+   PNTULIST& LoopList = Loop.GetUPointList() ;
+   // Ciclo sui segmenti del loop per cercare il tratto del loop chiuso più vicino al punto.   
+   double dMinSqDist = DBL_MAX ;
+   auto itMinDist = LoopList.end() ;
+   auto itSt = LoopList.begin() ;
+   auto itEn = itSt ;
+   ++ itEn ; 
+   for ( ; itSt != LoopList.end()  &&  itEn != LoopList.end() ; ++ itSt, ++ itEn) {
+      // Estremi del segmento corrente del loop
+      Point3d ptSegSt = itSt->first ;
+      Point3d ptSegEn = itEn->first ;
+      // Distanza del punto dal segmento del loop
+      DistPointLine dDistCalc( ptNewStart, ptSegSt, ptSegEn) ;
+      double dSqDist ;
+      dDistCalc.GetSqDist( dSqDist) ;
+      if ( dSqDist < dMinSqDist) {
+         dMinSqDist = dSqDist ;
+         itMinDist = itSt ;
+      }
+   }
+   // Se il punto non sta sul loop, errore
+   if ( dMinSqDist > dTol * dTol)
+      return false ;
+   // Se il punto non sta su un vertice del segmento, lo aggiungo. Altrimenti non devo fare nulla.
+   auto itNewPointSt = LoopList.begin() ;
+   auto itNext = itMinDist ;
+   ++ itNext ;
+   bool bOnStart = AreSamePointApprox( ptNewStart, itMinDist->first) ;
+   bool bOnEnd = AreSamePointApprox( ptNewStart, itNext->first) ;
+   itNewPointSt = LoopList.emplace( itNext, ptNewStart, 0) ;
+   // Sposto i punti precedenti in coda.
+   bool bStartRemoved = false ;
+   auto it = LoopList.begin() ;
+   while ( it != itNewPointSt) {
+      if ( bStartRemoved) {
+         LoopList.emplace_back( it->first, it->second) ;  
+      }
+      bStartRemoved = true ;
+      it = LoopList.erase( it) ;
+   }
+   // Se il punto inserito non coincide con l'inizio del segmento chiudo il loop.
+   if ( ! bOnStart) {
+      LoopList.emplace_back( ptNewStart, 0) ;
+      // Se coincide con la fine tolgo il punto di fine che diviene inutile.
+      if ( bOnEnd) {
+         //LoopList.erase( itNext) ;
+         auto itNewStart = LoopList.begin() ;
+         ++ itNewStart ;
+         LoopList.erase( itNewStart) ;
+      }
+   }
+
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+// nSegNum 0-based
+/*static*/ bool
+PointPositionOnPolyLine( const Point3d& ptPoint, /*const*/ PolyLine& Loop, int& nSegNum, double& dParOnSeg, double dTol)
+{
+   // Rinomino la lista di punti della PolyLine.
+   /*const*/ PNTULIST& LoopList = Loop.GetUPointList() ;
+   // Ciclo sui segmenti del loop per cercare il tratto del loop chiuso più vicino al punto.
+   nSegNum = - 1 ;
+   double dMinSqDist = DBL_MAX ;
+   int nS = 0 ;
+   auto itMinDistSt = LoopList.end() ;
+   auto itMinDistEn = itMinDistSt ;
+   auto itSt = LoopList.begin() ;
+   auto itEn = itSt ;
+   ++ itEn ; 
+   for ( ; itSt != LoopList.end()  &&  itEn != LoopList.end() ; ++ itSt, ++ itEn, ++ nS) {
+      // Estremi del segmento corrente del loop
+      Point3d ptSegSt = itSt->first ;
+      Point3d ptSegEn = itEn->first ;
+      // Distanza del punto dal segmento del loop
+      DistPointLine dDistCalc( ptPoint, ptSegSt, ptSegEn) ;
+      double dSqDist ;
+      dDistCalc.GetSqDist( dSqDist) ;
+      if ( dSqDist < dMinSqDist) {
+         nSegNum = nS ;
+         dMinSqDist = dSqDist ;
+         itMinDistSt = itSt ;
+         itMinDistEn = itEn ;
+      }
+   }
+   // Se il punto non sta sul loop, lo segnalo.
+   if ( dMinSqDist > dTol * dTol)
+      return false ;
+   // Calcolo il parametro lungo il segmento.
+   Vector3d vtSeg = itMinDistEn->first - itMinDistSt->first ;
+   double dSegLen = vtSeg.Len() ;
+   if ( dSegLen < EPS_SMALL)
+      return false ;
+   vtSeg /= dSegLen ;
+   dParOnSeg = Clamp( ( ptPoint - itMinDistSt->first) * vtSeg, 0., dSegLen) ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+/*static*/ bool
+IsPointInsidePolyLine( const Point3d& ptP, const PolyLine& plPoly)
+{
+   // Se la PolyLine non è chiusa, il punto non può essere interno.
+   if ( ! plPoly.IsClosed())
+      return false ;
+   // Lista dei punti
+   PolyLine& MyPlPoly = const_cast< PolyLine&>( plPoly) ;
+   const PNTULIST& List = MyPlPoly.GetUPointList() ;
+   // Ciclo sui segmenti della PolyLine per cercarne il tratto più vicino al punto.
+   double dMinSqDist = DBL_MAX ;
+   Point3d ptMinDist ;
+   auto itMinDistSt = List.end() ;
+   auto itSt = List.begin() ;
+   auto itEn = itSt ;
+   ++ itEn ;
+   for ( ; itSt != List.end()  &&  itEn != List.end() ; ++ itSt, ++ itEn) {
+      // Estremi del segmento corrente del loop
+      Point3d ptSegSt = itSt->first ;
+      Point3d ptSegEn = itEn->first ;
+      // Distanza del punto dal segmento del loop
+      DistPointLine dDistCalc( ptP, ptSegSt, ptSegEn) ;
+      double dSqDist ;
+      dDistCalc.GetSqDist( dSqDist) ;
+      if ( dSqDist < dMinSqDist) {
+         dMinSqDist = dSqDist ;
+         dDistCalc.GetMinDistPoint( ptMinDist) ;
+         itMinDistSt = itSt ;
+      }
+   }
+   // Termine del segmento di minima distanza.
+   auto itMinDistEn = itMinDistSt ;
+   ++ itMinDistEn ;
+   // Punto di minima distanza nell'estremo iniziale del segento
+   if ( AreSamePointApprox( ptMinDist, itMinDistSt->first)) {
+      auto itPrevSt = List.begin() ; 
+      if ( itMinDistSt == List.begin()) {
+         auto itAuxNext = itPrevSt ;
+         ++ ( ++ itAuxNext) ; 
+         for ( ; itAuxNext != List.end() ; ++ itPrevSt, ++ itAuxNext)
+            ;
+      }
+      else {
+         auto itAuxNext = itPrevSt ;
+         ++ itAuxNext ; 
+         for (; itAuxNext != itMinDistSt; ++itPrevSt, ++itAuxNext)
+            ;
+      }
+      Vector3d vtPrevTan = itMinDistSt->first - itPrevSt->first ;
+      vtPrevTan.Normalize() ;
+      Vector3d vtTan = itMinDistEn->first - itMinDistSt->first ;
+      vtTan.Normalize() ;
+      Polygon3d AuxPolygon ;
+      AuxPolygon.FromPolyLine( plPoly) ;
+      Vector3d vtPolyNorm = AuxPolygon.GetVersN() ;
+      Vector3d vtPrevOut = vtPrevTan ^ vtPolyNorm ;
+      Vector3d vtOut = vtTan ^ vtPolyNorm ;
+      // Caso concavo
+      if ( vtTan * vtPrevOut > 0) {
+         Vector3d vtTest = ptP - ptMinDist ;
+         if ( vtTest * vtPrevOut < 0  ||  vtTest * vtOut < 0)
+            return true ;
+      }
+      // Caso convesso
+      else {
+         Vector3d vtTest = ptP - ptMinDist ;
+         if ( vtTest * vtPrevOut < 0  &&  vtTest * vtOut < 0)
+            return true ;
+      }
+   }
+   // Punto di minima distanza nell'estremo finale del segento
+   else if ( AreSamePointApprox( ptMinDist, itMinDistEn->first)) {
+      auto itNextEn = itMinDistEn ;
+      ++ itNextEn ;
+      if ( itNextEn == List.end()) {
+         itNextEn = List.begin() ;
+         ++ itNextEn ;
+      }
+      Vector3d vtTan = itMinDistEn->first - itMinDistSt->first ;
+      vtTan.Normalize() ;
+      Vector3d vtNextTan = itNextEn->first - itMinDistEn->first ;
+      vtNextTan.Normalize() ;
+      Polygon3d AuxPolygon ;
+      AuxPolygon.FromPolyLine( plPoly) ;
+      Vector3d vtPolyNorm = AuxPolygon.GetVersN() ;
+      Vector3d vtOut =  vtTan ^ vtPolyNorm ;
+      Vector3d vtNextOut =  vtNextTan ^ vtPolyNorm ;
+      // Caso concavo
+      if ( vtNextTan * vtOut > 0) {
+         Vector3d vtTest = ptP - ptMinDist ;
+         if ( vtTest * vtOut < 0  ||  vtTest * vtNextOut < 0)
+            return true ;
+      }
+      // Caso convesso
+      else {
+         Vector3d vtTest = ptP - ptMinDist ;
+         if ( vtTest * vtOut < 0  &&  vtTest * vtNextOut < 0)
+            return true ;
+      }
+   }
+   // Punto di minima distanza interno al segmeno
+   else {
+      Vector3d vtP = ptP - itMinDistSt->first ;
+      Vector3d vtTan = itMinDistEn->first - itMinDistSt->first ;
+      vtTan.Normalize() ;
+      Polygon3d AuxPolygon ;
+      AuxPolygon.FromPolyLine( plPoly) ;
+      Vector3d vtPolyNorm = AuxPolygon.GetVersN() ;
+      Vector3d vtOut = vtTan ^ vtPolyNorm ;
+      vtP -= ( vtP * vtTan) * vtTan ;
+      if ( vtP * vtOut < - EPS_SMALL)
+         return true ;
+   }
+   return false ;
+}
+
+//----------------------------------------------------------------------------
+bool
+DistPointPolyLine( const Point3d& ptP, const PolyLine& plPoly, double& dPointPolyLineDist)
+{
+   if ( plPoly.GetPointNbr() == 0)
+      return false ;
+   dPointPolyLineDist = DBL_MAX ;
+   Point3d ptSt, ptEn ;
+   bool bContinue = plPoly.GetFirstPoint( ptSt)  &&  plPoly.GetNextPoint( ptEn) ;
+   while ( bContinue) {
+      double dPoinLineDist ;
+      DistPointLine PointLineDistCalc( ptP, ptSt, ptEn) ;
+      PointLineDistCalc.GetDist( dPoinLineDist) ;
+      if ( dPoinLineDist < dPointPolyLineDist)
+         dPointPolyLineDist = dPoinLineDist ;
+      ptSt = ptEn ;
+      bContinue = plPoly.GetNextPoint( ptEn) ;
+   }
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+/*static*/ bool
+SplitPolyLineAtPoint( const Point3d& ptPoint, /*const*/ PolyLine& Loop, PolyLine& Loop1, PolyLine& Loop2, double dTol)
+{
+   // Rinomino la lista di punti della PolyLine.
+   /*const*/ PNTULIST& LoopList = Loop.GetUPointList() ;
+   // Ciclo sui segmenti del loop per cercare il tratto del loop chiuso più vicino al punto.
+   double dMinSqDist = DBL_MAX ;
+   auto itMinDistSt = LoopList.end() ;
+   auto itSt = LoopList.begin() ;
+   auto itEn = itSt ;
+   ++ itEn ; 
+   for ( ; itSt != LoopList.end()  &&  itEn != LoopList.end() ; ++ itSt, ++ itEn) {
+      // Estremi del segmento corrente del loop
+      Point3d ptSegSt = itSt->first ;
+      Point3d ptSegEn = itEn->first ;
+      // Distanza del punto dal segmento del loop
+      DistPointLine dDistCalc( ptPoint, ptSegSt, ptSegEn) ;
+      double dSqDist ;
+      dDistCalc.GetSqDist( dSqDist) ;
+      if ( dSqDist < dMinSqDist) {
+         dMinSqDist = dSqDist ;
+         itMinDistSt = itSt ;
+      }
+   }
+   // Se il punto non sta sul loop, lo segnalo.
+   if ( dMinSqDist > dTol * dTol)
+      return false ;
+   // Se il punto di stop sta su un vertice non devo aggiungerlo e il
+   // punto di stop sarà uno degli estremi del segmento su cui giace.
+   auto itStop = itMinDistSt ;
+   auto itNext = itMinDistSt ;
+   ++ itNext ;
+   if ( AreSamePointApprox( ptPoint, itStop->first))
+      ;
+   else if ( AreSamePointApprox( ptPoint, itNext->first))
+      itStop = itNext ;
+   else {
+      itStop = LoopList.emplace( itNext, ptPoint, 0.) ;
+   }
+   // Creo i due loop
+   PNTULIST& LoopList1 = Loop1.GetUPointList() ;
+   PNTULIST& LoopList2 = Loop2.GetUPointList() ;
+   for ( auto it = LoopList.begin() ; it != itStop ; ++ it) {
+      LoopList1.emplace_back( it->first, it->second) ;
+   }
+   LoopList1.emplace_back( itStop->first, itStop->second) ;
+   for ( auto it = itStop ; it != LoopList.end() ; ++ it) {
+      LoopList2.emplace_back( it->first, it->second) ;
+   }
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+/*static*/ bool
+AddPolyLineToPolyLine( PolyLine& Poly, PolyLine& PolyToAdd, double dTol)
+{
+   // Se la PolyLine a cui devo aggiungere l'altra è chiusa, non posso aggiungere nulla.
+   if ( Poly.IsClosed())
+      return false ;
+   // Se la PolyLina che devo aggiungere è vuota, ho finito.
+   PNTULIST& PolyToAddList = PolyToAdd.GetUPointList() ;
+   if ( int( PolyToAddList.size()) == 0)
+      return true ;
+   // Se Poly non è vuota e la sua fine non coincide con l'inizio di PolyToAdd, non è possibile aggiungere nulla.
+   Point3d ptLast ;
+   Poly.GetLastPoint( ptLast) ;
+   auto it = PolyToAddList.begin() ;
+   if ( Poly.GetPointNbr() != 0  &&  ! AreSamePointEpsilon( it->first, ptLast, dTol))
+      return false ;
+   /*if ( Poly.GetPointNbr() == 0)
+   Poly.AddUPoint( 0., it->first) ;
+   ++ it ;*/
+   // Aggiungo i punti.
+   for ( ; it != PolyToAddList.end() ; ++ it) {
+      Poly.AddUPoint( 0., it->first) ;
+   }
+   return true ;
+}

SurfTriMeshBooleans.cpp

@@ -2804,66 +2804,6 @@ SurfTriMesh::IntersFacetFacet( const SurfFlatRegion& RegionA, const PolyLine& Ex
 //   return bOk ;
 //}
 
-//----------------------------------------------------------------------------
-static bool
-ChangePolyLineStart( const Point3d& ptNewStart, PolyLine& Loop)
-{
-   // Rinomino la lista di punti della PolyLine.
-   PNTULIST& LoopList = Loop.GetUPointList() ;
-   // Ciclo sui segmenti del loop per cercare il tratto del loop chiuso più vicino al punto.   
-   double dMinSqDist = DBL_MAX ;
-   auto itMinDist = LoopList.end() ;
-   auto itSt = LoopList.begin() ;
-   auto itEn = itSt ;
-   ++ itEn ; 
-   for ( ; itSt != LoopList.end()  &&  itEn != LoopList.end() ; ++ itSt, ++ itEn) {
-      // Estremi del segmento corrente del loop
-      Point3d ptSegSt = itSt->first ;
-      Point3d ptSegEn = itEn->first ;
-      // Distanza del punto dal segmento del loop
-      DistPointLine dDistCalc( ptNewStart, ptSegSt, ptSegEn) ;
-      double dSqDist ;
-      dDistCalc.GetSqDist( dSqDist) ;
-      if ( dSqDist < dMinSqDist) {
-         dMinSqDist = dSqDist ;
-         itMinDist = itSt ;
-      }
-   }
-   // Se il punto non sta sul loop, errore
-   if ( dMinSqDist > 100 * SQ_EPS_SMALL)
-      return false ;
-   // Se il punto non sta su un vertice del segmento, lo aggiungo. Altrimenti non devo fare nulla.
-   auto itNewPointSt = LoopList.begin() ;
-   auto itNext = itMinDist ;
-   ++ itNext ;
-   bool bOnStart = AreSamePointApprox( ptNewStart, itMinDist->first) ;
-   bool bOnEnd = AreSamePointApprox( ptNewStart, itNext->first) ;
-   itNewPointSt = LoopList.emplace( itNext, ptNewStart, 0) ;
-   // Sposto i punti precedenti in coda.
-   bool bStartRemoved = false ;
-   auto it = LoopList.begin() ;
-   while ( it != itNewPointSt) {
-      if ( bStartRemoved) {
-         LoopList.emplace_back( it->first, it->second) ;  
-      }
-      bStartRemoved = true ;
-      it = LoopList.erase( it) ;
-   }
-   // Se il punto inserito non coincide con l'inizio del segmento chiudo il loop.
-   if ( ! bOnStart) {
-      LoopList.emplace_back( ptNewStart, 0) ;
-      // Se coincide con la fine tolgo il punto di fine che diviene inutile.
-      if ( bOnEnd) {
-         //LoopList.erase( itNext) ;
-         auto itNewStart = LoopList.begin() ;
-         ++ itNewStart ;
-         LoopList.erase( itNewStart) ;
-      }
-   }
-      
-   return true ;
-}
-
 //----------------------------------------------------------------------------
 // nSegNum 0-based
 static bool
@@ -3024,27 +2964,6 @@ IsPointInsidePolyLine( const Point3d& ptP, /*const*/ PolyLine& plPoly)
    return false ;
 }
 
-//----------------------------------------------------------------------------
-bool
-DistPointPolyLine( const Point3d& ptP, const PolyLine& plPoly, double& dPointPolyLineDist)
-{
-   if ( plPoly.GetPointNbr() == 0)
-      return false ;
-   dPointPolyLineDist = DBL_MAX ;
-   Point3d ptSt, ptEn ;
-   bool bContinue = plPoly.GetFirstPoint( ptSt)  &&  plPoly.GetNextPoint( ptEn) ;
-   while ( bContinue) {
-      double dPoinLineDist ;
-      DistPointLine PointLineDistCalc( ptP, ptSt, ptEn) ;
-      PointLineDistCalc.GetDist( dPoinLineDist) ;
-      if ( dPoinLineDist < dPointPolyLineDist)
-         dPointPolyLineDist = dPoinLineDist ;
-      ptSt = ptEn ;
-      bContinue = plPoly.GetNextPoint( ptEn) ;
-   }
-   return true ;
-}
-
 //----------------------------------------------------------------------------
 // Una faccia di una trimesh ha una sola componente connessa.
 // Si assume che i loop siano corretti e rispettino tale proprietà.
@@ -3088,98 +3007,6 @@ DistPointFacet( const Point3d& ptP, /*const*/ POLYLINEVECTOR& vPolyVec, double&
    return true ;
 }
 
-//----------------------------------------------------------------------------
-static bool
-SplitPolyLineAtPoint( const Point3d& ptPoint, /*const*/ PolyLine& Loop, PolyLine& Loop1, PolyLine& Loop2)
-{
-   // Rinomino la lista di punti della PolyLine.
-   /*const*/ PNTULIST& LoopList = Loop.GetUPointList() ;
-   // Ciclo sui segmenti del loop per cercare il tratto del loop chiuso più vicino al punto.
-   double dMinSqDist = DBL_MAX ;
-   auto itMinDistSt = LoopList.end() ;
-   auto itSt = LoopList.begin() ;
-   auto itEn = itSt ;
-   ++ itEn ; 
-   for ( ; itSt != LoopList.end()  &&  itEn != LoopList.end() ; ++ itSt, ++ itEn) {
-      // Estremi del segmento corrente del loop
-      Point3d ptSegSt = itSt->first ;
-      Point3d ptSegEn = itEn->first ;
-      // Distanza del punto dal segmento del loop
-      DistPointLine dDistCalc( ptPoint, ptSegSt, ptSegEn) ;
-      double dSqDist ;
-      dDistCalc.GetSqDist( dSqDist) ;
-      if ( dSqDist < dMinSqDist) {
-         dMinSqDist = dSqDist ;
-         itMinDistSt = itSt ;
-      }
-   }
-   // Se il punto non sta sul loop, lo segnalo.
-   if ( dMinSqDist > 100 * SQ_EPS_SMALL)
-      return false ;
-   // Se il punto di stop sta su un vertice non devo aggiungerlo e il
-   // punto di stop sarà uno degli estremi del segmento su cui giace.
-   auto itStop = itMinDistSt ;
-   auto itNext = itMinDistSt ;
-   ++ itNext ;
-   if ( AreSamePointApprox( ptPoint, itStop->first))
-      ;
-   else if ( AreSamePointApprox( ptPoint, itNext->first))
-      itStop = itNext ;
-   else {
-      itStop = LoopList.emplace( itNext, ptPoint, 0.) ;
-   }
-   // Creo i due loop
-   PNTULIST& LoopList1 = Loop1.GetUPointList() ;
-   PNTULIST& LoopList2 = Loop2.GetUPointList() ;
-   for ( auto it = LoopList.begin() ; it != itStop ; ++ it) {
-      LoopList1.emplace_back( it->first, it->second) ;
-   }
-   LoopList1.emplace_back( itStop->first, itStop->second) ;
-   for ( auto it = itStop ; it != LoopList.end() ; ++ it) {
-      LoopList2.emplace_back( it->first, it->second) ;
-   }
-   return true ;
-}
-
-//----------------------------------------------------------------------------
-static bool
-AddPolyLineToPolyLine( PolyLine& Poly, PolyLine& PolyToAdd)
-{
-   // Se la PolyLine a cui devo aggiungere l'altra è chiusa, non posso aggiungere nulla.
-   if ( Poly.IsClosed())
-      return false ;
-   // Se la PolyLina che devo aggiungere è vuota, ho finito.
-   PNTULIST& PolyToAddList = PolyToAdd.GetUPointList() ;
-   if ( int( PolyToAddList.size()) == 0)
-      return true ;
-   // Se Poly non è vuota e la sua fine non coincide con l'inizio di PolyToAdd, non è possibile aggiungere nulla.
-   Point3d ptLast ;
-   Poly.GetLastPoint( ptLast) ;
-   auto it = PolyToAddList.begin() ;
-   if ( Poly.GetPointNbr() != 0  &&  ! AreSamePointEpsilon( it->first, ptLast, 10 * EPS_SMALL))
-      return false ;
-   /*if ( Poly.GetPointNbr() == 0)
-      Poly.AddUPoint( 0., it->first) ;
-   ++ it ;*/
-   // Aggiungo i punti.
-   for ( ; it != PolyToAddList.end() ; ++ it) {
-      Poly.AddUPoint( 0., it->first) ;
-   }
-   return true ;
-}
-
-//----------------------------------------------------------------------------
-struct PositionOnPolyLine {
-   int nIndexInVec ;
-   int nSegNum ;
-   double dParOnSeg ;
-   PositionOnPolyLine( int nIndex, int nSeg, double dPar) {
-      nIndexInVec = nIndex ;
-      nSegNum = nSeg ;
-      dParOnSeg = dPar ;
-   }
-} ;
-
 //----------------------------------------------------------------------------
 bool
 SurfTriMesh::SplitFacet( const INTERSCHAINMAP& IntersLineMap, PieceMap& NewFacet)

Triangulate.cpp

@@ -16,12 +16,21 @@
 #include "DllMain.h"
 #include "Triangulate.h"
 #include "ProjPlane.h"
+#include "CurveComposite.h"
+#include "CurveLine.h"
+#include "/EgtDev/Include/EGkIntersCurves.h"
+#include "/EgtDev/Include/EGkDistPointCurve.h"
 #include "/EgtDev/Include/EGkPolyLine.h"
 #include "/EgtDev/Include/EGkPlane3d.h"
 #include "/EgtDev/Include/EGkStringUtils3d.h"
+#include "Mathematics/TriangulateCDT.h"
 #include <algorithm>
 
+#include "tpp_interface.hpp"
+
 using namespace std ;
+using namespace tpp ;
+using namespace gte ;
 
 //----------------------------------------------------------------------------
 enum EarStatus{ EAS_NULL = -1, EAS_NO = 0, EAS_OK = 1} ;
@@ -35,18 +44,22 @@ static bool ChangeStartPntVector( int nNewStart, PNTVECTOR& vPi) ;
 //       INTVECTOR (int Vector) : 3*T indices of above points for T triangles   
 //----------------------------------------------------------------------------
 bool
-Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr)
+Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr, TrgType trgType)
 {
-  // verifico che la polilinea sia chiusa e piana e calcolo il piano medio del poligono
+   // verifico che la polilinea sia chiusa e piana e calcolo il piano medio del poligono
    double dArea ;
    Plane3d plPlane ;
    if ( ! PL.IsClosedAndFlat( plPlane, dArea, 50 * EPS_SMALL))
       return false ;
-  // determino il piano ottimale di proiezione e il relativo senso di rotazione
+
+   if ( trgType != TRG_STANDARD)
+      return Make( POLYLINEVECTOR{ PL}, vPt, vTr, trgType) ;
+
+   // determino il piano ottimale di proiezione e il relativo senso di rotazione
    bool bCCW ;
    if ( ! CalcProjPlane( plPlane.GetVersN(), m_nPlane, bCCW))
       return false ;
-  // riempio il vettore con i vertici del poligono da triangolare
+   // riempio il vettore con i vertici del poligono da triangolare
    vPt.clear() ;
    vPt.reserve( PL.GetPointNbr() - 1) ;
    // salto il primo punto (coincide con l'ultimo)
@@ -56,24 +69,24 @@ Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr)
    // inserisco i punti
    while ( PL.GetNextPoint( ptP))
       vPt.push_back( ptP) ;
-  // se non CCW inverto il vettore dei punti
+   // se non CCW inverto il vettore dei punti
    if ( ! bCCW)
       reverse( vPt.begin(), vPt.end()) ;
-  // creo il vettore degli indici del Poligono
+   // creo il vettore degli indici del Poligono
    INTVECTOR vPol ;
    int n = int( vPt.size()) ;
    vPol.reserve( n) ;
    for ( int i = 0 ; i < n ; ++ i)
       vPol.push_back( i) ;
 
-  // eseguo la triangolazione
+   // eseguo la triangolazione
    if ( ! MakeByEC2( vPt, vPol, vTr)  &&
-        ! MakeByEC3( vPt, vPol, vTr)) {
+      ! MakeByEC3( vPt, vPol, vTr)) {
       LOG_ERROR( GetEGkLogger(), "Error in MakeByEC23(1)")
-      return false ;
+         return false ;
    }
 
-  // se era CW, devo invertire il senso dei triangoli
+   // se era CW, devo invertire il senso dei triangoli
    if ( ! bCCW)
       reverse( vTr.begin(), vTr.end()) ;
 
@@ -86,7 +99,7 @@ Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr)
 //       INTVECTOR (int Vector) : 3*T indices of above points for T triangles   
 //----------------------------------------------------------------------------
 bool
-Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
+Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr, TrgType trgType)
 {
   // pulisco i vettori di ritorno
    vPt.clear() ;
@@ -115,6 +128,38 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
            ! AreOppositeVectorApprox( plExtPlane.GetVersN(), plPlane.GetVersN()))
          return false ;
    }
+
+  // triangolazione Delaunay 
+   if ( trgType == TRG_DEL_CONFORMING) {
+      if ( ! MakeByDelaunay( vPL, vPt, vTr, true, true)) {
+         LOG_ERROR( GetEGkLogger(), "Error in MakeByDelaunay ( conforming)") ;
+         return false ;
+      } 
+      return true ;
+   }
+   else if ( trgType == TRG_DEL_QUALITY) {
+      if ( ! MakeByDelaunay( vPL, vPt, vTr, false, true)) {
+         LOG_ERROR( GetEGkLogger(), "Error in MakeByDelaunay ( quality)") ;
+         return false ;
+      } 
+      return true ;   
+   }
+   else if ( trgType == TRG_DEL_NOQUALITY) {
+      if ( ! MakeByDelaunay( vPL, vPt, vTr, false, false)) {
+         LOG_ERROR( GetEGkLogger(), "Error in MakeByDelaunay ( no quality)") ;
+         return false ;
+      } 
+      return true ;  
+   }
+   else if ( trgType == TRG_DEL_GTE) {
+      if ( ! MakeByDelaunayGTE( vPL, vPt, vTr)) {
+         LOG_ERROR( GetEGkLogger(), "Error in MakeByDelaunay ( Geometric Tools Engine)") ;
+         return false ;
+      } 
+      return true ;  
+   }
+
+  // ear clipping 
   // se non CCW inverto tutte le polilinee
    if ( ! bCCW) {
       for ( int i = 0 ; i < int( vPL.size()) ; ++i)
@@ -172,7 +217,7 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
    INTVECTOR vPol ;
    int n = int( vPt.size()) ;
    vPol.reserve( n) ;
-   // non devo gestire separatamente CCW perchè ho già invertito i punti
+   // non devo gestire separatamente CCW perch� ho gi� invertito i punti
    for ( int i = 0 ; i < n ; ++ i)
       vPol.push_back( i) ;
 
@@ -190,6 +235,362 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
    return true ;
 }
 
+//----------------------------------------------------------------------------
+// Delaunay Triangulation ( Triangle library)
+//----------------------------------------------------------------------------
+bool 
+Triangulate::MakeByDelaunay( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr, bool bConforming, bool bQuality) 
+{
+  // sistema di riferimento locale sul piano della polyline  
+   Plane3d plPlane ;
+   double dArea;
+   if ( ! vPL[0].IsClosedAndFlat( plPlane, dArea, 50 * EPS_SMALL))
+      return false ;
+   Frame3d frLoc ;
+   frLoc.Set( plPlane.GetPoint(), plPlane.GetVersN()) ;
+
+   POLYLINEVECTOR vMyPL = vPL ;
+   for ( size_t t = 0 ; t < vPL.size() ; ++ t) {
+      vMyPL[t].ToLoc( frLoc) ; 
+      if ( vMyPL[t].GetPointNbr() < 4)   // per poter triangolare devo avere almeno tre punti distinti
+         return false ;
+   }
+
+  // vertici e constraint per la triangolazione 
+   vector<Delaunay::Point> vDelaunayVert ;
+   INTVECTOR vDelaunayConstr ; 
+
+   for ( size_t i = 0 ; i < vMyPL.size() ; i++) {      
+      Point3d ptFirst, pt ;
+      vMyPL[i].GetFirstPoint( ptFirst) ;
+      vDelaunayVert.push_back( Delaunay::Point( ptFirst.x, ptFirst.y)) ;
+      int nFirst = ( int)vDelaunayVert.size() - 1 ;  // indice del primo punto della polyline fra i vertici della triangolazione
+      vDelaunayConstr.push_back( nFirst) ;
+
+      while ( vMyPL[i].GetNextPoint( pt)) {
+         vDelaunayVert.push_back( Delaunay::Point( pt.x, pt.y)) ;
+        // nei constraint gli indici vanno inseriti due volte perchè vengono letti a due a due per definire un segmento
+         vDelaunayConstr.push_back( ( int)vDelaunayVert.size() - 1) ;
+         vDelaunayConstr.push_back( ( int)vDelaunayVert.size() - 1) ;        
+      }
+
+     // impongo che l'ultimo vertice coincida con il primo ( curva chiusa)
+      vDelaunayVert.pop_back() ;
+      vDelaunayConstr.erase(vDelaunayConstr.end() - 2, vDelaunayConstr.end()) ;
+      vDelaunayConstr.push_back( nFirst) ;      
+   }
+
+  // holes : sono definiti da un punto interno al buco  
+   std::vector<Delaunay::Point> vDelaunayHoles ;
+   for ( size_t i = 1 ; i < vMyPL.size() ; i++) {
+      Point3d pt ;      
+      CurveComposite * pCrvHole = CreateBasicCurveComposite() ;
+      pCrvHole->FromPolyLine( vMyPL[i]) ;
+      pCrvHole->GetCentroid( pt) ;
+     // se il centroide fosse esterno, cerco per tentativi un punto qualsiasi all'interno della curva 
+      double dPar = 0.5 ;
+      while ( ! IsPointInsidePolyLine( pt, vMyPL[i])) {
+         double dParS, dParE ;
+         pCrvHole->GetDomain( dParS, dParE) ;
+         if ( dPar > dParE)
+            return false ;         
+         Vector3d vtDir ;
+         pCrvHole->GetPointTang( dPar, ICurve::FROM_MINUS, pt, vtDir) ;
+         vtDir.Rotate( Z_AX, 0, -1) ;
+         pt += 2 * EPS_SMALL * vtDir ; 
+        // tento con un nuovo punto
+         dPar += 10 * EPS_SMALL ;
+      }
+
+      vDelaunayHoles.push_back( Delaunay::Point( pt.x, pt.y)) ;
+   }
+
+  // parti concave 
+   PNTVECTOR vPtConvexHull ;
+   vMyPL[0].GetConvexHullXY( vPtConvexHull) ;
+   CurveComposite* pCrv = CreateBasicCurveComposite() ;
+   pCrv->FromPolyLine( vMyPL[0]) ;
+
+   for ( size_t i = 0 ; i < vPtConvexHull.size() ; i ++) {
+      size_t NextIdx = ( i == vPtConvexHull.size() - 1) ? 0 : i + 1 ;
+      CurveLine * pLine = CreateBasicCurveLine() ;
+      pLine->Set( vPtConvexHull[i], vPtConvexHull[NextIdx]) ;
+
+     // verifico se ho regioni da eliminare 
+      IntersCurveCurve intCC( *pLine, *pCrv) ;
+      CRVCVECTOR ccClass ;
+      intCC.GetCurveClassification( 0, ccClass) ;
+      for ( size_t j = 0 ; j < ccClass.size() ; j ++) {
+         if ( ccClass[j].nClass == CRVC_OUT) {
+           // cerco per tentativi un punto a caso all'interno della regione da eliminare 
+            double dPar = ( ccClass[j].dParS + ccClass[j].dParE) / 2 ;
+            double dDist = 0 ;
+            Point3d pt ;
+            Vector3d vtDir ;
+            while ( dDist < EPS_SMALL) {
+               if ( dPar > ccClass[j].dParE)
+                  return false ;
+               pLine->GetPointTang( dPar, ICurve::FROM_MINUS, pt, vtDir) ; 
+               vtDir.Rotate( Z_AX, 0, 1) ;
+               DistPointCurve distPtCrv( pt, *pCrv) ;
+               dDist = 0.0 ;
+               distPtCrv.GetDist( dDist) ;
+               if ( dDist > EPS_SMALL) {
+                  pt += EPS_SMALL * ( vtDir) ;    
+                  vDelaunayHoles.push_back( Delaunay::Point( pt.x , pt.y)) ;      
+               }
+              // tento con nuovo punto 
+               dPar += 10 * EPS_SMALL ;
+            }
+         }
+      }
+   }
+
+  // triangolazione
+   Delaunay trGenerator( vDelaunayVert) ;
+   trGenerator.setSegmentConstraint( vDelaunayConstr) ;
+   trGenerator.setHolesConstraint( vDelaunayHoles) ;
+   if ( bConforming)
+      trGenerator.TriangulateConf( true) ;
+   else
+      trGenerator.Triangulate( bQuality) ;
+
+  // se non ho generato triangoli, errore
+   if ( trGenerator.ntriangles() == 0)
+      return false ;
+
+  // vertici
+   std::vector< Delaunay::Point> vVertex = trGenerator.MyVertexTraverse() ;
+   for (size_t i = 0; i < vVertex.size(); i++) {
+      Point3d pt( vVertex[i][0],  vVertex[i][1]) ;
+      pt.ToGlob( frLoc) ;
+      vPt.push_back( pt) ;
+   }
+
+  // triangoli
+   vTr = trGenerator.MyTriangleTraverse() ;
+
+   return true ;
+}
+
+
+//----------------------------------------------------------------------------
+// Delaunay Triangulation ( Geometric Tools Engine library)
+//----------------------------------------------------------------------------
+bool 
+Triangulate::MakeByDelaunayGTE( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr) 
+{
+  // sistema di riferimento locale sul piano della polyline
+   Plane3d plPlane ;
+   double dArea;
+   if ( ! vPL[0].IsClosedAndFlat( plPlane, dArea, 50 * EPS_SMALL))
+      return false ;
+   Frame3d frLoc ;
+   frLoc.Set( plPlane.GetPoint(), plPlane.GetVersN()) ;
+
+   POLYLINEVECTOR vMyPL = vPL ;
+   for ( size_t t = 0 ; t < vPL.size() ; ++ t) {
+      vMyPL[t].ToLoc( frLoc) ;
+      if ( vMyPL[t].GetPointNbr() < 4)
+         return false ;
+   }
+
+   CurveComposite * pCrvPL = CreateBasicCurveComposite() ;
+   pCrvPL->FromPolyLine( vMyPL[0]) ;
+
+  // controllo che holes siano contenuti nel loop esterno
+   for ( size_t i = 1 ; i < vMyPL.size() ; i ++) {      
+      CurveComposite * pCrvHole = CreateBasicCurveComposite() ;
+      pCrvHole->FromPolyLine( vMyPL[i]) ;
+
+      IntersCurveCurve intCC( *pCrvPL, *pCrvHole) ;
+      CRVCVECTOR ccClass ;
+      intCC.GetCurveClassification( 1, ccClass) ;
+
+      if ( ccClass.size() > 1 || ccClass[0].nClass != CRVC_IN) 
+         return false ; 
+   }
+
+   // Verifico non ci siano autointersezioni
+   SelfIntersCurve Inters( *pCrvPL) ;
+   if ( Inters.GetIntersCount() > 0 && Inters.GetOverlaps() == 0) 
+      return false ;
+
+   // Gestione autointersezione   
+   // if ( Inters.GetIntersCount() > 0 && Inters.GetOverlaps() == 0) {
+   //    
+   //   // vertici della triangolazione  
+   //    vector<Vector2<double>> vPoints ;   
+   //    Point3d pt ;
+   //    vMyPL[0].GetFirstPoint( pt) ;
+   //    vPoints.push_back( {pt.x, pt.y}) ;
+   //    pt.ToGlob( frLoc) ;
+   //    vPt.push_back( pt) ;
+   //    while ( vMyPL[0].GetNextPoint( pt) ) {
+   //       vPoints.push_back( {pt.x, pt.y}) ;
+   //       pt.ToGlob( frLoc) ;
+   //       vPt.push_back( pt) ;
+   //    }
+   //    vPt.pop_back() ; 
+   //    vPoints.pop_back() ;
+   //
+   //   // i punti di autointersezione P sono salvati come pair: parametro di P, indice posizione di P in vPoints
+   //    vector<pair<double,int>> dParams ; 
+   //    for ( int i = 0 ; i < Inters.GetIntersCount() ; i++) {
+   //       
+   //       IntCrvCrvInfo aInfo ;
+   //       Inters.GetIntCrvCrvInfo( i, aInfo) ;
+   //               
+   //       Point3d pt ; 
+   //       pt = aInfo.IciA[0].ptI ;
+   //       pt.ToGlob( frLoc) ;
+   //       
+   //       auto it = find_if( vPt.begin(), vPt.end(), [pt]( const Point3d& ptCrv){ return AreSamePointApprox( ptCrv, pt) ;} ) ;
+   //      // se il punto di interezione è uno dei vertici 
+   //       if ( it != vPt.end()) {
+   //          dParams.push_back( pair<double, int> (aInfo.IciA[0].dU, it - vPt.begin())) ;
+   //          dParams.push_back( pair<double, int> (aInfo.IciB[0].dU, it - vPt.begin())) ;      
+   //       }
+   //      // se non è uno dei vertici, lo aggiungo in vPt e in vPoints
+   //       else {
+   //          vPt.push_back( pt) ;
+   //          pt.ToLoc( frLoc) ;
+   //          vPoints.push_back( { pt.x, pt.y}) ;
+   //          dParams.push_back( pair<double, int> (aInfo.IciA[0].dU, vPoints.size() - 1)) ;
+   //          dParams.push_back( pair<double, int> (aInfo.IciB[0].dU, vPoints.size() - 1)) ;
+   //       }      
+   //    }
+   //    sort( dParams.begin(), dParams.end(), []( const pair<double, int>& a, const pair<double, int>& b){ return a.first < b.first ; }) ;
+   //    
+   //   // lati del poligono :
+   //    INTVECTOR vPolygonIdx ; 
+   //    for ( int i = 0 ; i < vMyPL[0].GetPointNbr() ; i++)
+   //       vPolygonIdx.push_back( i) ;
+   //     
+   //   // inserisco i nuovi lati determinati dai punti di intersezione
+   //    for ( int i = 0 ; i < dParams.size() ; i++) {
+   //      // se il parametro c'è già, non lo inserisco 
+   //       if ( find_if( vPolygonIdx.begin(), vPolygonIdx.end(),  
+   //                     [dParams, i]( const int& val) { return abs( val - dParams[i].first) < 0 ; }) != vPolygonIdx.end())
+   //           break ;
+   //      // altrimenti lo inserisco in modo ordinato nel vettore 
+   //       auto it = lower_bound( vPolygonIdx.begin(), vPolygonIdx.end(), dParams[i].first, 
+   //                              []( const int& a , const int& b) { return a < b ; }) ;
+   //       vPolygonIdx.insert( it, dParams[i].second) ;          
+   //    }
+   //
+   //    auto pPolygon = std::make_shared<PolygonTree>() ;         
+   //    pPolygon->polygon = vPolygonIdx ;
+   //    
+   //   // holes
+   //    int k = vPoints.size() ;
+   //    for ( size_t i = 1 ; i < vMyPL.size() ; i++) {
+   //       
+   //       auto pHole = std::make_shared<PolygonTree>() ; 
+   //       Point3d pt ;         
+   //       vMyPL[i].GetFirstPoint( pt) ;
+   //       vPoints.push_back( { pt.x, pt.y}) ;
+   //       pHole->polygon.push_back( k) ;
+   //       k ++ ;
+   //       while ( vMyPL[i].GetNextPoint( pt)) {
+   //          vPoints.push_back( { pt.x, pt.y}) ;
+   //          pHole->polygon.push_back( k) ;
+   //          k ++ ;
+   //       }
+   //    
+   //       vPoints.pop_back() ;
+   //       pHole->polygon.pop_back() ;
+   //       k -- ;       
+   //
+   //       pPolygon->child.push_back( pHole) ;
+   //    
+   //    }
+   //   
+   //   // Triangolazione 
+   //    PolygonTreeEx PTOutput ;
+   //    TriangulateCDT<double> triangulator ;
+   //    triangulator( vPoints, pPolygon, PTOutput) ;
+   //    
+   //    for ( auto const& tri : PTOutput.interiorTriangles) {
+   //       vTr.push_back( tri[0]) ;
+   //       vTr.push_back( tri[1]) ;
+   //       vTr.push_back( tri[2]) ;
+   //    } 
+   //    
+   //    return true ;
+   //
+   // }
+
+
+  // i vertici della triangolazione vPt sono i punti delle polyline nel sistema globale  
+   for ( size_t i = 0 ; i < vPL.size() ; i++) {
+      Point3d pt ;
+      vPL[i].GetFirstPoint( pt) ;
+      vPt.push_back( Point3d( pt.x, pt.y, pt.z)) ;
+      while ( vPL[i].GetNextPoint( pt ) ) 
+         vPt.push_back( Point3d( pt.x, pt.y, pt.z)) ; 
+     // ultimo punto della polyline coincide con il primo, quindi non serve 
+      vPt.pop_back() ;  
+   }
+
+   vector<Vector2<double>> vPoints ;                   // vertici per la triangolazione 
+   auto pPolygon = std::make_shared<PolygonTree>() ;   // poligono da triangolare
+
+   Point3d pt ;
+   int k = 0 ;   
+   vMyPL[0].GetFirstPoint( pt) ;
+   vPoints.push_back( { pt.x, pt.y}) ;
+   pPolygon->polygon.push_back( k) ;
+   k ++ ;
+
+   while ( vMyPL[0].GetNextPoint( pt) ) {
+      vPoints.push_back( { pt.x, pt.y}) ;
+      pPolygon->polygon.push_back( k) ;
+      k ++ ;
+   }
+
+  // ultimo punto della polyline coincide con il primo, quindi non serve
+   vPoints.pop_back() ;
+   pPolygon->polygon.pop_back() ;
+   k -- ;
+
+  // holes
+   for ( size_t i = 1 ; i < vMyPL.size() ; i++) {
+
+      auto pHole = std::make_shared<PolygonTree>() ; 
+      vMyPL[i].GetFirstPoint( pt) ;
+      vPoints.push_back( { pt.x, pt.y}) ;
+      pHole->polygon.push_back( k) ;
+      k ++ ;
+      while ( vMyPL[i].GetNextPoint( pt)) {
+         vPoints.push_back( { pt.x, pt.y}) ;
+         pHole->polygon.push_back( k) ;
+         k ++ ;
+      }
+
+      vPoints.pop_back() ;
+      pHole->polygon.pop_back() ;
+      k -- ;        
+
+      pPolygon->child.push_back( pHole) ;
+
+   }
+
+  // Triangolazione 
+   PolygonTreeEx PTOutput ;
+   TriangulateCDT<double> triangulator ;
+   triangulator( vPoints, pPolygon, PTOutput) ;
+
+   for ( auto const & tri : PTOutput.interiorTriangles)  {  
+      vTr.push_back( tri[0]) ;
+      vTr.push_back( tri[1]) ;
+      vTr.push_back( tri[2]) ;
+   } 
+
+   return true ;
+}
+
+
 //----------------------------------------------------------------------------
 bool
 Triangulate::PrepareGrid( const PNTVECTOR& vPt, const INTVECTOR& vPol,
@@ -273,7 +674,7 @@ Triangulate::MakeByEC( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& v
              vTr.push_back( vPol[i]) ;
              vTr.push_back( vPol[vNext[i]]) ;
           }
-         // ‘Delete’ vertex v[i] by redirecting next and previous links
+         // �Delete� vertex v[i] by redirecting next and previous links
          // of neighboring verts past it. Decrement vertex count
           vNext[vPrev[i]] = vNext[i] ;
           vPrev[vNext[i]] = vPrev[i] ;
@@ -393,7 +794,7 @@ Triangulate::MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR&
         // Reset earity of diagonal endpoints
          vEar[vPrev[i]] = EAS_NULL ;
          vEar[vNext[i]] = EAS_NULL ;
-        // ‘Delete’ vertex v[i] by redirecting next and previous links
+        // �Delete� vertex v[i] by redirecting next and previous links
         // of neighboring verts past it. Decrement vertex count
          vNext[vPrev[i]] = vNext[i] ;
          vPrev[vNext[i]] = vPrev[i] ;
@@ -506,7 +907,7 @@ Triangulate::MakeByEC3( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR&
         // Reset earity of diagonal endpoints
          vEar[vPrev[i]] = EAS_NULL ;
          vEar[vNext[i]] = EAS_NULL ;
-        // ‘Delete’ vertex v[i] by redirecting next and previous links
+        // �Delete� vertex v[i] by redirecting next and previous links
         // of neighboring verts past it. Decrement vertex count
          vNext[vPrev[i]] = vNext[i] ;
          vPrev[vNext[i]] = vPrev[i] ;
@@ -570,7 +971,7 @@ Triangulate::TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol,
       }
    }
    else {
-     // The ‘ear’ triangle is clockwise so v[i] is not an ear
+     // The �ear� triangle is clockwise so v[i] is not an ear
       bIsEar = false ;
    }
 
@@ -896,14 +1297,14 @@ Triangulate::GetOuterPntToJoin( const PNTVECTOR& vPt, const Point3d& ptP, int& n
          break ;
       }
    }
-  // non ho trovato alcunché, errore
+  // non ho trovato alcunch�, errore
    if ( nI == - 1)
       return false ;
   // se ho trovato un punto esatto del contorno, non devo fare altri controlli
    if ( AreSamePointApprox( ptInt, vPt[nI]))
       return true ;
   // devo ora verificare che il segmento che unisce i punti non intersechi altri lati del contorno esterno
-  // altrimenti tengo il punto con raggio più vicino a X_AX o Y_AX o Z_AX secondo m_nPlane
+  // altrimenti tengo il punto con raggio pi� vicino a X_AX o Y_AX o Z_AX secondo m_nPlane
    int nJ = nI ;
    Point3d ptPa = ptP ;
    Point3d ptPb = vPt[nI] ;
@@ -919,7 +1320,7 @@ Triangulate::GetOuterPntToJoin( const PNTVECTOR& vPt, const Point3d& ptP, int& n
    double dMinTan = INFINITO ;
    double dMinSqDist = SQ_INFINITO ;
    for ( int i = 0 ; i < nNumPt ; ++ i) {
-     // salto il punto già trovato
+     // salto il punto gi� trovato
       if ( i == nJ)
          continue ;
      // verifico se sta nel triangolo
@@ -955,7 +1356,7 @@ Triangulate::GetOuterPntToJoin( const PNTVECTOR& vPt, const Point3d& ptP, int& n
 bool
 Triangulate::PointInSector( const Point3d& ptTest, const Point3d& ptPrev, const Point3d& ptCorn, const Point3d& ptNext)
 {
-  // la parte valida del settore è a sinistra dei segmenti ptPrev --> ptCorn --> ptNext
+  // la parte valida del settore � a sinistra dei segmenti ptPrev --> ptCorn --> ptNext
   // se corner convesso
    if ( TriangleIsCCW( ptPrev, ptCorn, ptNext, 0))
       return ( TriangleIsCCW( ptPrev, ptCorn, ptTest)  &&
@@ -970,10 +1371,10 @@ Triangulate::PointInSector( const Point3d& ptTest, const Point3d& ptPrev, const
 bool
 ChangeStartPntVector( int nNewStart, PNTVECTOR& vPi)
 {
-  // se il nuovo inizio coincide col vecchio, non devo fare alcunché
+  // se il nuovo inizio coincide col vecchio, non devo fare alcunch�
    if ( nNewStart == 0)
       return true ;
-  // se il nuovo indice è oltre la dimensione del vettore, errore
+  // se il nuovo indice � oltre la dimensione del vettore, errore
    if ( nNewStart >= int( vPi.size()))
       return false ;
   // ciclo di aggiustamento

Triangulate.h

@@ -16,12 +16,20 @@
 #include "/EgtDev/Include/EGkPolyLine.h"
 #include "/EgtDev/Include/EGkPointGrid3d.h"
 
+// 
+enum TrgType { TRG_STANDARD,          // ear clipping
+               TRG_DEL_CONFORMING,    // conforming constrained Delaunay ( with quality constraint)
+               TRG_DEL_QUALITY,       // constrained Delaunay with quality constraints ( no angle smaller than 20 degrees)
+               TRG_DEL_NOQUALITY,     // constrained Delaunay without quality constraints 
+               TRG_DEL_GTE            // constrained Delaunay with Geomtric Tools Engine library
+             } ;
+
 //----------------------------------------------------------------------------
 class Triangulate
 {
    public :
-      bool Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr) ;
-      bool Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr) ;
+      bool Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr, TrgType trgType = TRG_STANDARD) ;
+      bool Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr, TrgType trgType = TRG_STANDARD) ;
 
    private :
       bool PrepareGrid( const PNTVECTOR& vPt, const INTVECTOR& vPol,
@@ -29,6 +37,8 @@ class Triangulate
       bool MakeByEC( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ;
       bool MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ;
       bool MakeByEC3( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ;
+      bool MakeByDelaunay( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr, bool bConforming, bool bQuality) ;
+      bool MakeByDelaunayGTE( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr) ;
       bool TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol,
                          const INTVECTOR& vPrev, INTVECTOR& vNext, int i) ;
       double CalcTriangleAspectRatio( const Point3d& ptPa, const Point3d& ptPb, const Point3d& ptPc) ;

dpoint.hpp

@@ -0,0 +1,635 @@
+/*! \file dpoint.hpp
+    \brief d-dimensional point class
+    
+    A d-dimensional point class which is written carefully using templates. It allows for basic
+    operations on points in any dimension. Orientation tests for 2 and 3 dimensional points are 
+    supported using <a href="http://www.cs.berkeley.edu/~jrs/">Jonathan's</a> code. This class 
+    forms the building block of other classes like dplane, dsphere etc.
+
+    \author <a href="www.compgeom.com/~piyush">Piyush Kumar</a>
+    \bug    No known bugs.
+ */
+
+#ifndef REVIVER_POINT_HPP
+#define REVIVER_POINT_HPP
+
+// changed mrkkrj --
+//#include "assert.hpp"
+#include "tpp_assert.hpp"
+// END changed --
+#include <iostream>
+#include <valarray>
+#include <stdio.h>
+#include <limits>
+
+
+//! The reviver namespace signifies the part of the code borrowed from reviver (dpoint.hpp). 
+namespace reviver {
+
+
+// Forward Declaration of the main Point Class
+// Eucledian d-dimensional point. The distance is L_2
+
+template<typename NumType, unsigned D>
+class dpoint;
+
+
+///////////////////////////////////////////////////////
+// Internal number type traits for dpoint
+///////////////////////////////////////////////////////
+
+template<typename T> 
+class InternalNumberType;
+
+template<>
+class InternalNumberType<float>{
+public: 
+  typedef double __INT;
+};
+
+
+template<>
+class InternalNumberType<int>{
+public: 
+  typedef long long __INT;
+};
+
+
+template<>
+class InternalNumberType<double>{
+public: 
+  typedef double __INT;
+};
+
+template<>
+class InternalNumberType<long>{
+public: 
+  typedef long long __INT;
+};
+
+
+///////////////////////////////////////////////////////
+// Origin of d-dimensional point
+///////////////////////////////////////////////////////
+template< typename NumType, unsigned D, unsigned I > struct origin
+{
+   static inline void eval( dpoint<NumType,D>& p )
+   {
+      p[I] = 0.0;
+      origin< NumType, D, I-1 >::eval( p );
+   }
+};
+
+
+// Partial Template Specialization
+template <typename NumType, unsigned D> struct origin<NumType, D, 0>
+{
+   static inline void eval( dpoint<NumType,D>& p )
+   {
+      p[0] = 0.0;
+   }
+};
+
+
+//! A structure to compute squared distances between points
+/*!
+    Uses unrolling of loops using templates.
+*/
+///////////////////////////////////////////////////////
+// Squared Distance of d-dimensional point
+///////////////////////////////////////////////////////
+template< typename NumType, unsigned D, unsigned I > struct Distance
+{
+   static inline double eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
+   {
+     double sum = ((double)p[I] -  (double)q[I] ) *( (double)p[I] - (double)q[I] );
+      return sum + Distance< NumType, D, I-1 >::eval( p,q );
+   }
+};
+
+
+//!  Partial Template Specialization for distance calculations
+template <typename NumType, unsigned D> struct Distance<NumType, D, 0>
+{
+   static inline double eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
+   {
+     return ((double) p[0] - (double)q[0] )*( (double)p[0] - (double)q[0] );
+   }
+};
+
+
+//! A structure to compute dot product between two points or associated vectors
+/*!
+    Uses unrolling of loops using templates.
+*/
+///////////////////////////////////////////////////////
+// Dot Product of two d-dimensional points
+///////////////////////////////////////////////////////
+template< typename __INT, typename NumType, unsigned D, unsigned I > struct DotProd
+{
+   static inline __INT eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
+   {
+      __INT sum = ( ((__INT)p[I]) * ((__INT)q[I]) );
+      return sum + DotProd< __INT, NumType, D, I-1 >::eval( p,q );
+   }
+};
+
+
+//! Partial Template Specialization for dot product calculations
+template < typename __INT, typename NumType, unsigned D> struct DotProd<__INT,NumType, D, 0>
+{
+   static inline __INT eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
+   {
+      return ( ((__INT)p[0]) * ((__INT)q[0]) );
+   }
+};
+
+
+///////////////////////////////////////////////////////
+// Equality of two d-dimensional points
+///////////////////////////////////////////////////////
+template< typename NumType, unsigned D, unsigned I > struct IsEqual
+{
+   static inline bool eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
+   {
+      if( p[I]  != q[I] ) return false;
+      else return IsEqual< NumType, D, I-1 >::eval( p,q );
+   }
+};
+
+
+// Partial Template Specialization
+template <typename NumType, unsigned D> struct IsEqual<NumType, D, 0>
+{
+   static inline NumType eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
+   {
+       return (p[0] == q[0])?1:0;
+   }
+};
+
+
+
+//!   Equate two d-dimensional points. 
+/*!
+    Uses unrolling of loops using templates. 
+    A class used to implement operator= for points. This class also helps in automatic type 
+    conversions of points (with explicit calls for conversion).
+*/
+template< 
+  typename NumType1, 
+  typename NumType2, 
+  unsigned D, 
+  unsigned I 
+  > struct Equate
+{
+   static inline void eval( dpoint<NumType1,D>& p,const dpoint<NumType2,D>& q )
+   {
+      p[I]  = q[I];
+      Equate< NumType1, NumType2, D, I-1 >::eval( p,q );
+   }
+};
+
+
+//! Partial Template Specialization for Equate
+template <
+  typename NumType1, 
+  typename NumType2,
+  unsigned D
+  > struct Equate<NumType1,NumType2, D, 0>
+{
+   static inline void  eval( dpoint<NumType1,D>& p,const dpoint<NumType2,D>& q )
+   {
+       p[0] = q[0];
+   }
+};
+
+
+//! A structure to add two points
+/*!
+    Uses unrolling of loops using templates.
+*/
+///////////////////////////////////////////////////////
+// Add two d-dimensional points
+///////////////////////////////////////////////////////
+template< typename NumType, unsigned D, unsigned I > struct Add
+{
+   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
+   {
+      result[I] = p[I]  + q[I];
+      Add< NumType, D, I-1 >::eval( result,p,q );
+   }
+};
+
+
+//! Partial Template Specialization for Add structure
+template <typename NumType, unsigned D> struct Add<NumType, D, 0>
+{
+   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
+   {
+       result[0] = p[0] + q[0];
+   }
+};
+
+
+///////////////////////////////////////////////////////
+// Subtract two d-dimensional points
+///////////////////////////////////////////////////////
+// Could actually be done using scalar multiplication and addition
+
+
+// What about unsigned types?
+template< typename NumType > 
+inline NumType Subtract_nums(const NumType& x, const NumType& y) {
+  if(!std::numeric_limits<NumType>::is_signed) {
+      std::cerr << "Exception: Can't subtract unsigned types."; exit(1);
+  }
+  return x - y;
+}
+
+
+//!   Subtract two d-dimensional vectors
+/*!
+      Caution: Do not use on unsigned types.
+*/
+template< typename NumType, unsigned D, unsigned I > struct Subtract
+{
+   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
+   {
+     
+          result[I] = Subtract_nums(p[I] , q[I]);
+      Subtract< NumType, D, I-1 >::eval( result,p,q );
+   }
+};
+
+
+//! Partial Template Specialization for subtraction of points (associated vectors)
+template <typename NumType, unsigned D> struct Subtract<NumType, D, 0>
+{
+   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
+   {
+       result[0] = Subtract_nums(p[0] , q[0]);
+   }
+};
+
+
+
+
+
+//!   Mutiply scalar with d-dimensional point
+/*!
+      Scalar mulipltication of d-dimensional point with a number using template unrolling.
+*/
+template< typename NumType, unsigned D, unsigned I > struct Multiply
+{
+   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, NumType k)
+   {
+      result[I] = p[I] * k;
+      Multiply< NumType, D, I-1 >::eval( result,p,k );
+   }
+};
+
+
+//! Partial Template Specialization for scalar multiplication
+template <typename NumType, unsigned D> struct Multiply<NumType, D, 0>
+{
+   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, NumType k )
+   {
+       result[0] = p[0] * k;
+   }
+};
+
+
+
+//!  Main d dimensional Point Class
+/*!
+    -  NumType = Floating Point Type
+    -  D       = Dimension of Point
+*/
+template<typename NumType = double, unsigned D = 3>
+class dpoint {
+
+        // Makes Swap operation fast
+        NumType  x[D];
+
+public:
+        typedef NumType NT;
+        typedef typename InternalNumberType<NumType>::__INT __INT;
+
+    // To be defined in a cpp file
+    //  const MgcVector2 MgcVector2::ZERO(0,0);
+    //  static const dpoint<NumType,D> Zero;
+
+    inline void move2origin(){ origin<NumType, D, D-1>::eval(*this); };
+
+    dpoint(){ 
+        Assert( (D >= 1), "Dimension < 1 not allowed" ); 
+        // move2origin(); 
+    };
+
+    //! 1 D Point
+    dpoint(NumType x0){ x[0] = x0; };
+    //! 2 D Point
+    dpoint(NumType x0,NumType x1){ x[0] = x0;  x[1] = x1; };
+    //! 3 D Point
+    dpoint(NumType x0,NumType x1,NumType x2){  x[0] = x0;  x[1] = x1; x[2] = x2; };
+    //! Array Initialization
+    dpoint(NumType ax[]){ for(int i =0; i < D; ++i) x[i] = ax[i]; };
+    //! Initialization from another point : Copy Constructor
+        dpoint(const dpoint<NumType,D>& p){  Equate<NumType,NumType,D,D-1>::eval((*this),p);	};
+
+         
+    //! Automatic type conversions of points.
+    //! Only allowed if the conversion is specified explicitly by the programmer.
+    template<class OtherNumType>
+        explicit dpoint(const dpoint<OtherNumType,D>& p){ Equate<NumType,OtherNumType,D,D-1>::eval((*this),p); };
+
+    // Destructor
+    ~dpoint(){};
+
+    inline int      dim() const { return D; };
+    inline double   sqr_dist(const dpoint<NumType,D> q) const ;
+    inline double   distance(const dpoint<NumType,D> q) const ;
+    inline __INT    dotprod (const dpoint<NumType,D> q) const ;
+    inline __INT    sqr_length(void)  const;
+    inline void     normalize (void);
+    inline NumType& operator[](int i);
+    inline NumType  operator[](int i) const;
+
+    inline dpoint&  operator= (const dpoint<NumType,D>& q);
+
+    template<typename NT, unsigned __DIM>
+    friend dpoint<NT,__DIM>   operator- (const dpoint<NT,__DIM>& p, const dpoint<NT,__DIM>& q);
+
+    template<typename NT, unsigned __DIM>
+    friend dpoint<NT,__DIM>   operator+ (const dpoint<NT,__DIM>& p, const dpoint<NT,__DIM>& q);
+
+    template<typename NT, unsigned __DIM>
+    friend bool   operator== (const dpoint<NT,__DIM>& p, const dpoint<NT,__DIM>& q);
+
+    template<typename NT, unsigned __DIM>
+    friend bool   operator!= (const dpoint<NT,__DIM>& p, const dpoint<NT,__DIM>& q);
+
+
+//	inline dpoint&  operator= (const valarray<NumType>& v);
+//	inline operator valarray<NumType>() const;
+
+    template<typename __NT,unsigned __DIM>
+    friend void iswap(dpoint<__NT,__DIM>& p,dpoint<__NT,__DIM>& q);
+};
+
+template<typename NumType, unsigned D>
+void dpoint<NumType,D>::normalize (void){
+    double len = sqrt(sqr_length());
+    if (len > 0.00001)
+    for(int i = 0; i < D; ++i){
+        x[i] /= len;
+    }
+}
+
+/*
+template<typename NumType, unsigned D>
+dpoint<NumType,D>::operator valarray<NumType>() const{
+    valarray<NumType> result((*this).x , D);
+    return result;
+}
+
+//Warning : Valarray should be of size D
+//TODO: Unwind this for loop into a template system
+template<typename NumType, unsigned D>
+dpoint<NumType,D>&
+dpoint<NumType,D>::operator= (const valarray<NumType>& v){
+    dpoint<NumType,D> result;
+    for(int i = 0; i < D; i++) (*this).x[i] = v[i];
+    return (*this);
+}
+*/
+
+template<typename NT, unsigned __DIM>
+dpoint<NT,__DIM>
+operator+ (const dpoint<NT,__DIM>& p, const dpoint<NT,__DIM>& q){
+    dpoint<NT,__DIM> result;
+    Add<NT,__DIM,__DIM-1>::eval(result,p,q);	
+    return result;
+}
+
+template<typename NT, unsigned __DIM>
+dpoint<NT,__DIM>
+operator- (const dpoint<NT,__DIM>& p, const dpoint<NT,__DIM>& q){
+    dpoint<NT,__DIM> result;
+    // cout << "Subtracting..." << p << " from " << q << " = ";
+    Subtract<NT,__DIM,__DIM-1>::eval(result,p,q);	
+    // cout << result << endl;	
+    return result;
+}
+
+template<typename NT, unsigned __DIM>
+bool
+operator== (const dpoint<NT,__DIM>& p, const dpoint<NT,__DIM>& q){
+    return IsEqual<NT,__DIM,__DIM-1>::eval(p,q);	
+}
+
+template<typename NT, unsigned __DIM>
+bool
+operator!= (const dpoint<NT,__DIM>& p, const dpoint<NT,__DIM>& q){
+    return !(IsEqual<NT,__DIM,__DIM-1>::eval(p,q));	
+}
+
+template<typename NT, unsigned __DIM>
+dpoint<NT,__DIM>
+operator* (const dpoint<NT,__DIM>& p, const NT k){
+    dpoint<NT,__DIM> result;
+    Multiply<NT,__DIM,__DIM-1>::eval(result,p,k);	
+    return result;
+}
+
+template<typename NT, unsigned __DIM>
+dpoint<NT,__DIM>
+operator/ (const dpoint<NT,__DIM>& p, const NT k){
+    Assert( (k != 0), "Hell division by zero man...\n");
+    dpoint<NT,__DIM> result;
+    Multiply<NT,__DIM,__DIM-1>::eval(result,p,((double)1.0)/k);	
+    return result;
+}
+
+template < typename NumType, unsigned D >
+dpoint<NumType,D>&
+dpoint<NumType,D>::operator=(const dpoint<NumType,D> &q)
+{
+  Assert((this != &q), "Error p = p");
+  Equate<NumType,NumType,D,D-1>::eval(*this,q);	
+  return *this;
+}
+
+template < typename NumType, unsigned D >
+NumType
+dpoint<NumType,D>::operator[](int i) const
+{ return x[i]; }
+
+template < typename NumType, unsigned D >
+NumType&
+dpoint<NumType,D>::operator[](int i)
+{
+  return x[i]; 
+}
+
+
+template<typename NumType, unsigned D>
+double
+dpoint<NumType,D>::sqr_dist (const dpoint<NumType,D> q) const {
+    return Distance<NumType,D,D-1>::eval(*this,q);	
+}
+
+template<typename NumType, unsigned D>
+double 
+dpoint<NumType,D>::distance (const dpoint<NumType,D> q) const {
+    return sqrt(static_cast<double>(Distance<NumType,D,D-1>::eval(*this,q)));	
+}
+
+
+template<typename NumType, unsigned D>
+typename dpoint<NumType,D>::__INT
+dpoint<NumType,D>::dotprod (const dpoint<NumType,D> q) const {
+    return DotProd<__INT,NumType,D,D-1>::eval(*this,q);	
+}
+
+template<typename NumType, unsigned D>
+typename dpoint<NumType,D>::__INT
+dpoint<NumType,D>::sqr_length (void) const {
+#ifdef _DEBUG	
+    if( DotProd<__INT,NumType,D,D-1>::eval(*this,*this) < 0) {
+          std::cerr << "Point that caused error: ";
+      std::cerr << *this << std::endl;
+      std::cerr << DotProd<__INT,NumType,D,D-1>::eval(*this,*this) << std::endl;
+      std::cerr << "Fatal: Hell!\n"; exit(1);
+    }
+#endif
+    return DotProd<__INT,NumType,D,D-1>::eval(*this,*this);	
+    
+}
+
+template < class NumType, unsigned D >
+std::ostream&
+operator<<(std::ostream& os,const dpoint<NumType,D> &p)
+{
+     os << "Point (d=";
+     os << D << ", (";
+     for (unsigned int i=0; i<D-1; ++i)
+        os << p[i] << ", ";
+    return os << p[D-1] << "))";
+    
+};
+
+template < class NumType, unsigned D >
+std::istream&
+operator>>(std::istream& is,dpoint<NumType,D> &p)
+{
+     for (int i=0; i<D; ++i)
+         if(!(is >> p[i])){
+             if(!is.eof()){
+                std::cerr << "Error Reading Point:" 
+                      << is << std::endl;
+                exit(1);
+             }
+         }
+         
+    return is;
+    
+};
+
+/*
+template<typename __NT,unsigned __DIM>
+static inline void iswap(dpoint<__NT,__DIM>& p,dpoint<__NT,__DIM>& q){
+    __NT *y;
+    y = p.x;
+    p.x = q.x;
+    q.x = y;
+}
+*/
+
+
+
+template < typename NumType, unsigned D >
+dpoint<NumType, D> CrossProd(const dpoint<NumType, D>& vector1, 
+                 const dpoint<NumType, D>& vector2) {
+   Assert(D == 3, "Cross product only defined for 3d vectors");
+   dpoint<NumType, D> vector;
+   vector[0] = (vector1[1] * vector2[2]) - (vector2[1] * vector1[2]);
+   vector[1] = (vector2[0] * vector1[2]) - (vector1[0] * vector2[2]);
+   vector[2] = (vector1[0] * vector2[1]) - (vector2[0] * vector1[1]); 
+   return vector;
+}
+
+
+
+
+template < typename __NT, unsigned __DIM >
+int
+orientation(const dpoint<__NT,__DIM> p[__DIM+1])
+{
+    int _sign = + 1;
+    // To be implemented
+    std::cerr << "Not yet implemented\n";
+    exit(1);
+    return _sign;
+    
+}
+
+
+template < typename __NT >
+inline __NT
+orientation(
+        const dpoint<__NT,2>& p,
+        const dpoint<__NT,2>& q,
+        const dpoint<__NT,2>& r
+        )
+{
+   // 2D speaciliazation for orientation
+    std::cout << "FATAL";
+  exit(1);
+  return ((p[0]-r[0])*(q[1]-r[1]))-((q[0]-r[0])*(p[1]-r[1]));
+}
+
+
+extern "C" double orient2d(double *p, double *q, double *r);
+
+template < >
+inline double
+orientation<double>(
+        const dpoint<double,2>& p,
+        const dpoint<double,2>& q,
+        const dpoint<double,2>& r
+        )
+{
+   // 2D speaciliazation for orientation
+  double pp[2] = { p[0], p[1] };
+  double qq[2] = { q[0], q[1] };
+  double rr[2] = { r[0], r[1] };
+  return orient2d(pp,qq,rr);
+}
+
+
+template < >
+inline float
+orientation<float>(
+        const dpoint<float,2>& p,
+        const dpoint<float,2>& q,
+        const dpoint<float,2>& r
+        )
+{
+   // 2D speaciliazation for orientation
+  double pp[2] = { p[0], p[1] };
+  double qq[2] = { q[0], q[1] };
+  double rr[2] = { r[0], r[1] };
+  return (float)orient2d(pp,qq,rr);
+}
+
+
+
+};    // Namespace Ends here
+
+
+
+
+#endif
+
+

tpp_assert.cpp

@@ -0,0 +1,43 @@
+/*! \file assert.cpp
+    \brief Implements a better 'Assert'
+ */
+
+#include "stdafx.h"
+
+
+#include <iostream>
+#include <stdlib.h>
+
+#if _WINDOWS
+#include <cassert>
+#endif
+
+namespace tpp {
+
+/*! \def MyAssertFunction
+    \brief Function used by 'Assert' function in _DEBUG mode.
+   
+    Details.
+*/
+bool MyAssertFunction( bool b, const char* desc, int line, const char* file){
+    // changed mrkkrj --
+#if _WINDOWS
+    (void)desc;
+    (void)line;
+    (void)file;
+    assert(b); // use integration with Visual Studio!
+    (void)b;
+   return true;
+#else
+    // Original:
+    if (b) return true;
+    std::cerr << "\n\nAssertion Failure\n";
+    std::cerr << "Description : " << desc << std::endl;
+    std::cerr << "Filename    : " << file << std::endl;
+    std::cerr << "Line No     : " << line << std::endl;
+    exit(1);
+#endif
+}
+
+
+} // end of namespace

tpp_assert.hpp

@@ -0,0 +1,29 @@
+/*! \file assert.hpp
+    \brief Implements a better 'Assert'. 
+           Used in the reviver::dpoint inplementation.
+ */
+
+namespace tpp {
+
+#ifndef REVIVER_ASSERT_HPP
+#define REVIVER_ASSERT_HPP
+
+/*! \def MyAssertFunction
+    \brief Function used by 'Assert' function in _DEBUG mode.
+*/
+extern bool MyAssertFunction( bool b, const char* desc, int line, const char* file);
+
+
+#if defined( _DEBUG )
+// changed mrkkrj --
+//#define Assert( exp, description ) MyAssertFunction( (int)(exp), description, __LINE__, __FILE__ )
+#define Assert( exp, description ) tpp::MyAssertFunction( (int)(exp), description, __LINE__, __FILE__ )
+// END changed --
+#else
+#define Assert( exp, description )
+#endif
+
+
+#endif
+
+}

tpp_interface.hpp

@@ -0,0 +1,719 @@
+
+/*! \file tpp_interface.hpp
+    \brief The main Delaunay C++ class of the Triangle++ wrapper.
+
+    Use this class to produce Delaunay triangulations.
+
+    The following description pertains to the original version, the current version
+    was ported to VisualStudio. Thus it doesn't need Python scripts, and is supposed
+    to be used *as it is* in your program!
+ */
+
+
+ /*! \mainpage Triangle++
+ \section intro Introduction
+ <table border="0">
+ <tr><td>
+ If you do not know, what a Delaunay triangulation is, you can read more about it
+ <a href="http://www.compgeom.com/~piyush/teach/5930/slides/lecture8.ppt">here</a> and
+ <a href="http://en.wikipedia.org/wiki/Delaunay_triangulation">here</a>.
+ This C++ library module is just a wrapper class on the
+ <a href="http://www.cs.berkeley.edu/~jrs/">Triangle</a>
+ package of <a href="http://www.cs.berkeley.edu/~jrs/">Jonathan</a>.
+
+ Many times I have had to use triangle in C++ code bases of mine and have been forced to use C.
+ At last I thought I would put a wrapper on his cool C code and it seems that this is what I got.
+
+ The design is not perfect and the code was written in a day, but it does compile and run on the
+ machines I tried (cygwin/redhat). The C++ wrapper will certainly slow access down if you want to
+ mess with the triangulation but the basic delaunay triangulation should be as fast as triangle.
+
+ Look at the tpp_interface.hpp file for getting started on what this wrapper can do for you. Also
+ have a look at main.cpp which shows an example of using this class. The class is thread-safe.
+ <p>
+ <b>Requirements</b> : Python, make and C++ compilers.
+ Supported C/C++ Compilers: g++ / icpc (Intel C++).
+ Also needed is doxygen for generating documentation.
+ </p>
+ <p>
+ <b>Compilation</b> : Just type 'make'</p>
+ <p>
+ <b>Testing</b> : Goto the bin directory, and type './test ../data/input.dat' (after compilation of course).
+ </p>
+ </td>
+ <td><img src="http://upload.wikimedia.org/wikipedia/en/9/92/Delaunay_triangulation.png" alt="Delaunay Triangulation Example"></td>
+ </tr>
+ </table>
+
+ \section Downloads
+ You can download the latest version of the source code from <a href="triangle++.tar.gz">here</a>.
+
+ \section authors Authors
+     <ul>
+         <li><a href="http://compgeom.com/~piyush">Piyush Kumar</a></li>
+         <li><a href="http://www.ib-krajewski.de">Marek Krajewski</a></li>
+         <li>Hopefully more to come... (please feel free to extend this wrapper)</li>
+     </ul>
+
+ \section changelog Change Log
+
+ 17/04/20: mrkkrj – added support Voronoi tesselation <br>
+ 22/01/20: mrkkrj – added support for custom constraints (angle and area) <br>
+ 17/09/18: mrkkrj – ported to 64-bit (preliminary, not thorougly tested!) <br>
+ 11/07/11: mrkkrj - bugfix in Triangle's divandconqdelauney() <br>
+ 10/15/11: mrkkrj - added support for the "quality triangulation" option, added some debug support<br>
+ 08/24/11: mrkkrj - Ported to VisualStudio, added comp. operators, reformatted and added some comments<br>
+ 10/21/06: Replaced vertexsort with C++ sort.<br>
+ 10/25/06: Wrapped in tpp namespace for usage with other libraries with similar names.
+           Added some more documentation/small changes. Used doxygen 1.5.0 and dot. Tested compilation with
+           icc 9.0/9.1, gcc-4.1/3.4.6. <br>
+ 11/03/06: Fixed the compilation system. <br>
+
+ \todo
+ <ol>
+     <li> Intel Compiler Warnings with -Wall </li>
+     <ul>
+          <li> remove the compiler warnings for icpc 9.0/9.1</li>
+     </ul>
+
+     <li> Implement vertexmedian() in C++. </li>
+     <li> Implement the flip operator as a member function of Delaunay. </li>
+ </ol>
+ */
+
+ //-----------------------------------------------------------
+
+#ifndef TRPP_INTERFACE
+#define TRPP_INTERFACE
+
+// changed mrkkrj --
+//#include <dpoint.hpp>
+#include "dpoint.hpp"
+// END changed --
+#include <vector>
+#include <string>
+
+//! The main namespace in which the Triangle++ project lives
+namespace tpp {
+
+   // (mrkkrj)
+   enum DebugOutputLevel {
+      None, 
+      Info,    // most useful; it gives information on algorithmic progress and much more detailed statistics
+      Vertex,  // gives vertex-by-vertex details, and prints so much that Triangle runs much more slowly
+      Debug    // gives information only a debugger could love
+   };
+
+   //!  The main Delaunay Class that wraps around Triangle.
+   /*!
+     This is a C++ wrapper of the Triangle package by JRS.
+
+     This class currently uses the dpoint class written by me (the point class is a d-dimensional point
+     class reviver::dpoint (but for this application it only uses the d=2 case).
+     Additionally, the inner helper C++ class Triwrap groups the original Triangle's C functions.
+
+     \author Piyush Kumar, mrkkrj
+
+     \note   (mrkkrj) For for backgroud info on the Triangle's implementation see "Triangle:
+             Engineering a 2D Quality Mesh Generator and Delaunay Triangulator" by JP Shewchuk:
+             www.cs.cmu.edu/~quake-papers/triangle.ps
+   */
+   class Delaunay {
+   public:
+      //! Point Typedef
+      /*! Warning: If you want to use your own point class, you might have to
+          work hard...
+           - mrkkrj: true!!! - spare your time, use an adapter class.
+      */
+      typedef reviver::dpoint <double, 2> Point;
+
+      //! The main constructor.
+      /*!
+        Takes a vector of 2 dimensional points where each of the coordinates is
+        expressed as double.
+      */
+      Delaunay(const std::vector<Point>& points);
+
+      //! The main destructor.
+      /*!
+        Does memory cleanup mostly.
+      */
+      ~Delaunay();
+
+      //! Delaunay Triangulate the input points (Quality)
+      /*!
+        This function calls Triangle.h to Delaunay-triangulate points given as input to the 
+        constructor of this class. Here a Quality triangualtion will be created.
+
+        \param quality enforce minimal angle (default: 20°) and, minimal area (only if explicitely set)
+      */
+      void Triangulate(bool quality = false, DebugOutputLevel = None);
+
+      //! Delaunay Triangulate the input points (Conforming)
+      /*!
+        This function calls Triangle.h to Delaunay-triangulate points given as input to the
+        constructor of this class. Here a Conforming triangualtion will be created.
+      */
+      void TriangulateConf(bool quality = false, DebugOutputLevel = None);
+
+      //! Voronoi-tesselate the input points (added  mrkkrj)
+      /*!
+        This function calls triangle to create a Voronoi diagram with points given as input
+        to the constructor of this class.
+
+        Note that a Voronoi diagram can be only created if the underlying triangulation is convex 
+        and doesn't have holes!
+
+        \param useConformingDelaunay use conforming Delaunay triangulation as base for the Voronoi diagram
+      */
+      void Tesselate(bool useConformingDelaunay = false, DebugOutputLevel traceLvl = None);
+
+      //! Set a quality constraint for the triangulation
+      /*!
+        \param angle min. resulting angle, if angle <= 0, the constraint will be removed.
+      */
+      void setMinAngle(float angle) {
+         m_minAngle = angle;
+      }
+
+      //! Set a quality constraint for the triangulation
+      /*!
+        \param area max. triangle area, if area <= 0, the constraint will be removed.
+      */
+      void setMaxArea(float area) {
+         m_maxArea = area;
+      }
+
+      //! Set the segments to constrain the triangulation
+      /*!
+        Takes a vector of 2 dimensional points where each consecutive pair of points describes a single segment.
+      Both endpoints of every segment are vertices of the input vector, and a segment may intersect other segments 
+      and vertices only at its endpoints.
+
+      \return true if the input is valid, false otherwise 
+      */
+      bool setSegmentConstraint(const std::vector<Point>& segments);
+
+     //! Set the segments to constrain the triangulation
+     /*!
+       Same as above, but usign indexes of the input points!
+
+     \return true if the input is valid, false otherwise
+     */
+     bool setSegmentConstraint(const std::vector<int>& segmentPointIndexes);
+
+     //! Set the holes to constrain the triangulation
+     /*!
+       Takes a vector of 2 dimensional points where each consecutive pair of points describes a single edge of a hole.
+
+     \return true if the input is valid, false otherwise
+     */
+     bool setHolesConstraint(const std::vector<Point>& holes);
+
+      //! Are the quality constrainst sane?
+      /*!
+        \possible true if is highly probable for triangualtion to succeed
+        \return true if triangualtion is guaranteed to succeed
+      */
+      bool checkConstraints(bool& possible) const;
+
+      //! Are the quality constrainst sane, take two
+      /*!
+        \relaxed report highly probable as correct too, as error otherwise
+        \return true if triangualtion is guaranteed to succeed, or at least higly probable to
+      */
+      bool checkConstraintsOpt(bool relaxed) const;
+
+      //! Get minAngle intervals
+      /*!
+        \guaranteed up to this value triangualtion is guaranteed to succeed
+        \possible up to this value it is highly probable for triangualtion to succeed
+      */
+      static void getMinAngleBoundaries(float& guaranteed, float& possible);
+
+      //! Set a user test function for the triangulation
+      /*!
+        OPEN TODO::: (use the -u switch!!!!)
+      */
+      void setUserConstraint(bool (*f)()) {};
+
+      //! Output a geomview .off file containing the delaunay triangulation
+      /*!
+        \param fname output file name.
+      */
+      void writeoff(std::string& fname);
+
+      //! Number of edges in the triangulation
+      /*!
+        \return Number of Edges
+        Remember to call Triangulate before using this function.
+      */
+      int nedges() const;
+
+      //! Number of triangles in the triangulation
+      /*!
+        \return Number of Triangles
+        Remember to call Triangulate before using this function.
+      */
+      int ntriangles() const;
+
+      //! Number of vertices in the triangulation
+      /*!
+        \return Number of Vertices
+        Remember to call Triangulate before using this function.
+      */
+      int nvertices() const;
+
+      //! Number of vertices on the convex hull.
+      /*!
+        \return Number of vertices on the convex hull.
+        Remember to call Triangulate before using this function.
+      */
+      int hull_size() const;
+
+      //! Number of Voronoi points in the tesselation
+      /*!
+        \return Number of Points
+        Remember to call Tesselate before using this function.
+      */
+      int nvpoints() const;
+
+      //! Number of Voronoi edges in the tesselation
+      /*!
+        \return Number of Edges
+        Remember to call Tesselate before using this function.
+      */
+      int nvedges() const;
+
+
+      ///////////////////////////////
+      //
+      // Vertex Iterator
+      //
+      ///////////////////////////////
+
+      //!  The vertex iterator for the Delaunay class
+      class vIterator {
+      private:
+         vIterator(Delaunay* tiangulator);   //! To set container
+
+         Delaunay* MyDelaunay;   //! Which container do I point
+         void* vloop;            //! Triangles Internal data.
+
+      public:
+         vIterator operator++();
+         vIterator() :vloop(nullptr) {};
+         Point& operator*() const;
+         ~vIterator();
+
+         friend class Delaunay;
+         friend bool operator==(vIterator const&, vIterator const&);
+         friend bool operator!=(vIterator const&, vIterator const&);
+      };
+
+      //! Vertex iterator begin function
+      vIterator vbegin() { return vIterator(this); };
+      //! Vertex iterator end function
+      vIterator vend();
+
+      //! Given an iterator, find its index in the input vector of points.
+      int vertexId(vIterator const& vit) const;
+
+      //! Given an index, return the actual double Point
+      const Point& point_at_vertex_id(int i) { return m_PList[i]; };
+
+      //! Return the Point additionally created in quality mesh generation ("q" option)
+      Point added_point_at_vertex_id(int i);
+
+      friend class vIterator;
+
+
+      ///////////////////////////////
+      //
+      // Face Iterator
+      //
+      ///////////////////////////////
+
+      //!  The face iterator for the Delaunay class
+      class fIterator {
+      private:
+
+         struct tdata {
+            double*** tri;
+            int orient;
+         };
+
+         typedef struct tdata  poface;
+
+         fIterator(Delaunay* tiangulator);  //! To set container
+
+         Delaunay* MyDelaunay;   //! Which container do I point
+         //void *floop;          //! Triangles Internal data.
+         poface floop;
+
+      public:
+         void operator++();
+         fIterator() { floop.tri = nullptr; };
+         ~fIterator();
+
+         friend class Delaunay;
+         friend bool operator==(fIterator const&, fIterator const&);
+         friend bool operator!=(fIterator const&, fIterator const&);
+         friend bool operator<(fIterator const&, fIterator const&); // added mrkkrj
+      };
+
+      //! Face iterator begin function
+      fIterator fbegin() { return fIterator(this); };
+      //! Face iterator end function
+      fIterator fend();
+      int faceId(fIterator const&);
+
+
+      //! Access the origin (Org) vertex of a face.
+      /*!
+        \param fit  Face interator.
+        \param point  if specified: the cordinates of the vertex
+        \return Index of the vertex in m_pList,
+                or -1 if quality option was used and a new vertex was created!
+
+        A triangle abc has origin (org) a,destination (dest) b, and apex (apex)
+        c.  These vertices occur in counterclockwise order about the triangle.
+        Remember to call Triangulate before using this function. Do not use it on a null iterator.
+      */
+      int Org(fIterator const& fit, Point* point = 0);
+
+
+      //! Access the destination (Dest) vertex of a face.
+      /*!
+        \param fit  Face interator.
+        \param point  if specified: the cordinates of the vertex
+        \return Index of the vertex in m_pList,
+                or -1 if quality option was used and a new vertex was created!
+
+        A triangle abc has origin (org) a,destination (dest) b, and apex (apex)
+        c.  These vertices occur in counterclockwise order about the triangle.
+        Remember to call Triangulate before using this function. Do not use it on a null iterator.
+      */
+      int Dest(fIterator const& fit, Point* point = 0);
+
+
+      //! Access the apex (Apex) vertex of a face.
+      /*!
+        \param fit  Face interator.
+        \param point  if specified: the cordinates of the vertex
+        \return Index of the vertex in m_pList,
+                or -1 if quality option was used and a new vertex was created!
+
+        A triangle abc has origin (org) a,destination (dest) b, and apex (apex)
+        c.  These vertices occur in counterclockwise order about the triangle.
+        Remember to call Triangulate before using this function. Do not use it on a null iterator.
+      */
+      int Apex(fIterator const& fit, Point* point = 0);
+
+
+      //! Access the triangle adjoining edge i
+      /*!
+        \param fit  Face Iterator
+        \param i  edge number
+        \return The vertex on the opposite face, or -1 (see Org() above)
+
+        A triangle abc has origin (org) a,destination (dest) b, and apex (apex)
+            c.  These vertices occur in counterclockwise order about the triangle.
+        <ul>
+        <li>sym(abc, 0) -> ba*</li>
+        <li>sym(abc, 1) -> cb*</li>
+        <li>sym(abc, 2) -> ac*</li>
+        </ul>
+        * is the farthest vertex on the adjoining triangle whose index
+        is returned. A -1 is returned if the edge is part of the convex hull.
+        Remember to call Triangulate before using this function.
+        Do not use it on a null iterator.
+      */
+      int Sym(fIterator const& fit, char i);
+
+
+      //! Access the triangle opposite to current edge of the face
+      /*!
+        \param fit  Face iterator
+        \return The iterator of the opposite face
+
+        A triangle abc has origin (org) a,destination (dest) b, and apex (apex)
+            c.  These vertices occur in counterclockwise order about the triangle.
+        The iterator
+        to the triangle is returned. The iterator is empty if the edge
+        is on the convex hull.
+            Remember to call Triangulate before using this function.
+        Do not use it on a null iterator.
+      */
+      fIterator Sym(fIterator const& fit);
+
+
+      //! Is the iterator empty?
+      /*!
+        \param fit  Face interator.
+        \return true if the iterator is empty
+      */
+      inline bool empty(fIterator const& fit)
+      {
+         return fit.floop.tri == nullptr;
+      };
+
+
+      //! Is the iterator pointing to the dummy triangle?
+      /*!
+        \param fit  Face interator.
+        \return true if the iterator is of the dummy triangle.
+      */
+      bool isdummy(fIterator const& fit);
+
+
+      //! Find the next edge (counterclockwise) of a triangle.   
+      /*!
+        \param fit  face iterator
+        \return The face iterator corresponding to the next counterclockwise edge of a triangle
+
+        Lnext(abc) -> bca.
+        Remember to call Triangulate before using this function.
+        Do not use it on a null iterator.
+      */
+      fIterator Lnext(fIterator const& fit);
+
+
+      //! Find the previous edge (clockwise) of a triangle.   
+      /*!
+        \param fit  face iterator
+        \return The face iterator corresponding to the previous clockwise edge of a triangle
+
+        Lprev(abc) -> cab.
+        Remember to call Triangulate before using this function.
+        Do not use it on a null iterator.
+      */
+      fIterator Lprev(fIterator const& fit);
+
+
+      //! Find the next edge (counterclockwise) of a triangle with the same origin
+      /*!
+        \param fit  face iterator
+        \return The face iterator corresponding to the next edge counterclockwise with the same origin.
+
+        Onext(abc) -> ac*.
+        Remember to call Triangulate before using this function.
+        Do not use it on a null iterator.
+      */
+      fIterator Onext(fIterator const& fit);
+
+
+      //! Find the next edge clockwise with the same origin.
+      /*!
+        \param fit  face iterator
+        \return The face iterator corresponding to the next edge clockwise with the same origin.
+
+        Onext(abc) -> a*b.
+        Remember to call Triangulate before using this function.
+        Do not use it on a null iterator.
+      */
+      fIterator Oprev(fIterator const& fit);
+
+
+      // TODO List: (for face iterators)
+      /*  dnext:  Find the next edge counterclockwise with the same destination.   */
+      /*  dnext(abc) -> *ba                                                        */
+      /*                                                                           */
+      /*  dprev:  Find the next edge clockwise with the same destination.          */
+      /*  dprev(abc) -> cb*                                                        */
+      /*                                                                           */
+      /*  rnext:  Find the next edge (counterclockwise) of the adjacent triangle.  */
+      /*  rnext(abc) -> *a*                                                        */
+      /*                                                                           */
+      /*  rprev:  Find the previous edge (clockwise) of the adjacent triangle.     */
+      /*  rprev(abc) -> b**                                                        */
+
+
+      //! Calculate incident triangles around a vertex.
+      /*!
+      \param vertexid The vertex for which you want incident triangles.
+            \param ivv Returns triangles around a vertex in counterclockwise order.
+
+      Note that behaviour is undefined if vertexid is greater than
+      number of vertices - 1. Remember to call Triangulate before using this function.
+      All triangles returned have Org(triangle) = vertexid.
+      All triangles returned are in counterclockwise order.
+      */
+      void trianglesAroundVertex(int vertexid, std::vector<int>& ivv);
+
+
+      //! Calculate the area of a face. 
+      /*!
+        \param fit  Face interator.
+        \return area of the face associated with the iterator.
+
+      */
+      double area(fIterator const& fit);
+
+
+      //! Point locate a vertex v
+      /*!
+        \param vertexid  vertex id
+        \return a face iterator whose origin is v.
+      */
+      fIterator locate(int vertexid); // OPEN:: doesn't seem to be working!
+
+
+      ///////////////////////////////
+      //
+      // Voronoi Points Iterator 
+      //  (added mrkkrj)
+      //
+      ///////////////////////////////
+
+      //!  The Voronoi points iterator for the Delaunay class
+      class vvIterator {
+      public:
+         vvIterator();
+         vvIterator operator++();
+         Point& operator*() const;
+         void advance(int steps);
+
+      private:
+         vvIterator(Delaunay* tiangulator);   //! To set container
+
+         Delaunay* m_delaunay;    //! Which container do I point to
+         void* vvloop;            //! Triangle's Internal data.
+         int vvindex;
+         int vvcount;
+
+         friend class Delaunay;
+         friend bool operator==(vvIterator const&, vvIterator const&);
+         friend bool operator!=(vvIterator const&, vvIterator const&);
+      };
+
+      //! Voronoi Points iterator begin function
+      vvIterator vvbegin() { return vvIterator(this); };
+      //! Voronoi Points iterator end function
+      vvIterator vvend();
+
+
+      ///////////////////////////////
+      //
+      // Voronoi Edges Iterator 
+      //  (added mrkkrj)
+      //
+      ///////////////////////////////
+
+      //!  The Voronoi edges iterator for the Delaunay class
+      class veIterator {
+      public:
+         veIterator();
+         veIterator operator++();
+         int startPointId() const;
+         int endPointId(Point& normvec) const;
+
+      private:
+         veIterator(Delaunay* tiangulator);   //! To set container
+
+         Delaunay* m_delaunay;   //! Which container do I point to
+         void* veloop;           //! Triangle's Internal data.
+         int veindex;
+         int vecount;
+
+         friend class Delaunay;
+         friend bool operator==(veIterator const&, veIterator const&);
+         friend bool operator!=(veIterator const&, veIterator const&);
+      };
+
+      //! Voronoi Points iterator begin function
+      veIterator vebegin() { return veIterator(this); };
+      //! Voronoi Points iterator end function
+      veIterator veend();
+
+
+      //! Access the origin (Org) vertex of an edge. (added mrkkrj)
+      /*!
+        \param eit  Voronoi Edge iterator.
+        \return The start point of the Voronoi edge,
+
+        Remember to call Tesselate before using this function. Do not use it on a null iterator.
+      */
+      const Point& Org(veIterator const& eit);
+
+
+      //! Access the destination (Dest) vertex of an edge. (added mrkkrj)
+      /*!
+        \param eit  Voronoi Edge iterator.
+        \param finiteEdge true for finite edges, false for inifinte rays.
+        \return The end point of the Voronoi edge, for infinite rays the normal vector of the ray
+
+        Remember to call Tesselate before using this function. Do not use it on a null iterator.
+      */
+      Point Dest(veIterator const& eit, bool& finiteEdge);
+
+
+      //--------------------------------------
+      // added mrkkrj - helper for Points 
+      //    OPEN:: compiler cannot instantiate less<> with operator<() for Point, why?!
+      //--------------------------------------
+      struct OrderPoints
+      {
+         bool operator() (const Point& lhs, const Point& rhs) const {
+            // first sort on x, then on y coordinates
+            if (lhs[0] < rhs[0]) {
+               return true;
+            }
+            if (lhs[0] == rhs[0] && lhs[1] < rhs[1]) {
+               return true;
+            }
+            return false;
+         }
+      };
+
+
+     //////////////////////////////////////////////////////////////////////////////////////////////
+      std::vector< Delaunay::Point> MyVertexTraverse( ) ;
+      std::vector< int> Delaunay::MyTriangleTraverse( ) ;
+     //////////////////////////////////////////////////////////////////////////////////////////////
+
+
+   private:
+      void Triangulate(std::string& triswitches);
+
+      // added mrkkrj - helper functions for face iterator access methods 
+      //    HACK:: double* as not to export internal impl.
+      void SetPoint(Point& point, double* vertexptr);
+      int GetVertexIndex(fIterator const& fit, double* vertexptr);
+      int GetFirstIndexNumber() const;
+
+      // added mrkkrj 
+      std::string formatFloatConstraint(float f) const;
+      void setDebugLevelOption(std::string& options, DebugOutputLevel traceLvl);
+      void freeTriangleDataStructs();
+
+      friend class fIterator;
+
+   private:
+      std::vector<Point> m_PList;   /*! Stores the input point list. */
+      void* m_in;             /*! Used for intput to triangle  */
+      void* m_triangleWrap;   /*! Triangle impl. is wrapped in this pointer. */
+      void* m_pmesh;          /*! pointer to triangle mesh */
+      void* m_pbehavior;
+      bool m_Triangulated;
+
+      // added mrkkrj:
+      void* m_vorout;          /*! pointer to Voronoi output */
+
+      // added mrkkrj: quality constraints
+      float m_minAngle;
+      float m_maxArea;
+
+     // added mrkkrj: segment constraints
+     std::vector<int> m_SList;
+     // added mrkkrj: holes 
+     std::vector<Point> m_HList;
+
+   }; // Class Delaunay
+
+} // namespace tpp ends.
+
+
+#endif
+
+

triangle.h

@@ -0,0 +1,300 @@
+
+/*! \file triangle.h
+    \brief Original Triangle package's include file.
+
+    Exports triangulateio structure for use in tpp_impl.hpp. You should not 
+    use struct triangulateio in your application if you are using Triangle++ 
+    wrapper! 
+ */
+
+/*****************************************************************************/
+/*                                                                           */
+/*  (triangle.h)                                                             */
+/*                                                                           */
+/*  Include file for programs that call Triangle.                            */
+/*                                                                           */
+/*  Accompanies Triangle Version 1.6                                         */
+/*  July 28, 2005                                                            */
+/*                                                                           */
+/*  Copyright 1996, 2005                                                     */
+/*  Jonathan Richard Shewchuk                                                */
+/*  2360 Woolsey #H                                                          */
+/*  Berkeley, California  94705-1927                                         */
+/*  jrs@cs.berkeley.edu                                                      */
+/*                                                                           */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/*                                                                           */
+/*  How to call Triangle from another program                                */
+/*                                                                           */
+/*                                                                           */
+/*  If you haven't read Triangle's instructions (run "triangle -h" to read   */
+/*  them), you won't understand what follows.                                */
+/*                                                                           */
+/*  Triangle must be compiled into an object file (triangle.o) with the      */
+/*  TRILIBRARY symbol defined (generally by using the -DTRILIBRARY compiler  */
+/*  switch).  The makefile included with Triangle will do this for you if    */
+/*  you run "make trilibrary".  The resulting object file can be called via  */
+/*  the procedure triangulate().                                             */
+/*                                                                           */
+/*  If the size of the object file is important to you, you may wish to      */
+/*  generate a reduced version of triangle.o.  The REDUCED symbol gets rid   */
+/*  of all features that are primarily of research interest.  Specifically,  */
+/*  the -DREDUCED switch eliminates Triangle's -i, -F, -s, and -C switches.  */
+/*  The CDT_ONLY symbol gets rid of all meshing algorithms above and beyond  */
+/*  constrained Delaunay triangulation.  Specifically, the -DCDT_ONLY switch */
+/*  eliminates Triangle's -r, -q, -a, -u, -D, -Y, -S, and -s switches.       */
+/*                                                                           */
+/*  IMPORTANT:  These definitions (TRILIBRARY, REDUCED, CDT_ONLY) must be    */
+/*  made in the makefile or in triangle.c itself.  Putting these definitions */
+/*  in this file (triangle.h) will not create the desired effect.            */
+/*                                                                           */
+/*                                                                           */
+/*  The calling convention for triangulate() follows.                        */
+/*                                                                           */
+/*      void triangulate(triswitches, in, out, vorout)                       */
+/*      char *triswitches;                                                   */
+/*      struct triangulateio *in;                                            */
+/*      struct triangulateio *out;                                           */
+/*      struct triangulateio *vorout;                                        */
+/*                                                                           */
+/*  `triswitches' is a string containing the command line switches you wish  */
+/*  to invoke.  No initial dash is required.  Some suggestions:              */
+/*                                                                           */
+/*  - You'll probably find it convenient to use the `z' switch so that       */
+/*    points (and other items) are numbered from zero.  This simplifies      */
+/*    indexing, because the first item of any type always starts at index    */
+/*    [0] of the corresponding array, whether that item's number is zero or  */
+/*    one.                                                                   */
+/*  - You'll probably want to use the `Q' (quiet) switch in your final code, */
+/*    but you can take advantage of Triangle's printed output (including the */
+/*    `V' switch) while debugging.                                           */
+/*  - If you are not using the `q', `a', `u', `D', `j', or `s' switches,     */
+/*    then the output points will be identical to the input points, except   */
+/*    possibly for the boundary markers.  If you don't need the boundary     */
+/*    markers, you should use the `N' (no nodes output) switch to save       */
+/*    memory.  (If you do need boundary markers, but need to save memory, a  */
+/*    good nasty trick is to set out->pointlist equal to in->pointlist       */
+/*    before calling triangulate(), so that Triangle overwrites the input    */
+/*    points with identical copies.)                                         */
+/*  - The `I' (no iteration numbers) and `g' (.off file output) switches     */
+/*    have no effect when Triangle is compiled with TRILIBRARY defined.      */
+/*                                                                           */
+/*  `in', `out', and `vorout' are descriptions of the input, the output,     */
+/*  and the Voronoi output.  If the `v' (Voronoi output) switch is not used, */
+/*  `vorout' may be NULL.  `in' and `out' may never be NULL.                 */
+/*                                                                           */
+/*  Certain fields of the input and output structures must be initialized,   */
+/*  as described below.                                                      */
+/*                                                                           */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/*                                                                           */
+/*  The `triangulateio' structure.                                           */
+/*                                                                           */
+/*  Used to pass data into and out of the triangulate() procedure.           */
+/*                                                                           */
+/*                                                                           */
+/*  Arrays are used to store points, triangles, markers, and so forth.  In   */
+/*  all cases, the first item in any array is stored starting at index [0].  */
+/*  However, that item is item number `1' unless the `z' switch is used, in  */
+/*  which case it is item number `0'.  Hence, you may find it easier to      */
+/*  index points (and triangles in the neighbor list) if you use the `z'     */
+/*  switch.  Unless, of course, you're calling Triangle from a Fortran       */
+/*  program.                                                                 */
+/*                                                                           */
+/*  Description of fields (except the `numberof' fields, which are obvious): */
+/*                                                                           */
+/*  `pointlist':  An array of point coordinates.  The first point's x        */
+/*    coordinate is at index [0] and its y coordinate at index [1], followed */
+/*    by the coordinates of the remaining points.  Each point occupies two   */
+/*    REALs.                                                                 */
+/*  `pointattributelist':  An array of point attributes.  Each point's       */
+/*    attributes occupy `numberofpointattributes' REALs.                     */
+/*  `pointmarkerlist':  An array of point markers; one int per point.        */
+/*                                                                           */
+/*  `trianglelist':  An array of triangle corners.  The first triangle's     */
+/*    first corner is at index [0], followed by its other two corners in     */
+/*    counterclockwise order, followed by any other nodes if the triangle    */
+/*    represents a nonlinear element.  Each triangle occupies                */
+/*    `numberofcorners' ints.                                                */
+/*  `triangleattributelist':  An array of triangle attributes.  Each         */
+/*    triangle's attributes occupy `numberoftriangleattributes' REALs.       */
+/*  `trianglearealist':  An array of triangle area constraints; one REAL per */
+/*    triangle.  Input only.                                                 */
+/*  `neighborlist':  An array of triangle neighbors; three ints per          */
+/*    triangle.  Output only.                                                */
+/*                                                                           */
+/*  `segmentlist':  An array of segment endpoints.  The first segment's      */
+/*    endpoints are at indices [0] and [1], followed by the remaining        */
+/*    segments.  Two ints per segment.                                       */
+/*  `segmentmarkerlist':  An array of segment markers; one int per segment.  */
+/*                                                                           */
+/*  `holelist':  An array of holes.  The first hole's x and y coordinates    */
+/*    are at indices [0] and [1], followed by the remaining holes.  Two      */
+/*    REALs per hole.  Input only, although the pointer is copied to the     */
+/*    output structure for your convenience.                                 */
+/*                                                                           */
+/*  `regionlist':  An array of regional attributes and area constraints.     */
+/*    The first constraint's x and y coordinates are at indices [0] and [1], */
+/*    followed by the regional attribute at index [2], followed by the       */
+/*    maximum area at index [3], followed by the remaining area constraints. */
+/*    Four REALs per area constraint.  Note that each regional attribute is  */
+/*    used only if you select the `A' switch, and each area constraint is    */
+/*    used only if you select the `a' switch (with no number following), but */
+/*    omitting one of these switches does not change the memory layout.      */
+/*    Input only, although the pointer is copied to the output structure for */
+/*    your convenience.                                                      */
+/*                                                                           */
+/*  `edgelist':  An array of edge endpoints.  The first edge's endpoints are */
+/*    at indices [0] and [1], followed by the remaining edges.  Two ints per */
+/*    edge.  Output only.                                                    */
+/*  `edgemarkerlist':  An array of edge markers; one int per edge.  Output   */
+/*    only.                                                                  */
+/*  `normlist':  An array of normal vectors, used for infinite rays in       */
+/*    Voronoi diagrams.  The first normal vector's x and y magnitudes are    */
+/*    at indices [0] and [1], followed by the remaining vectors.  For each   */
+/*    finite edge in a Voronoi diagram, the normal vector written is the     */
+/*    zero vector.  Two REALs per edge.  Output only.                        */
+/*                                                                           */
+/*                                                                           */
+/*  Any input fields that Triangle will examine must be initialized.         */
+/*  Furthermore, for each output array that Triangle will write to, you      */
+/*  must either provide space by setting the appropriate pointer to point    */
+/*  to the space you want the data written to, or you must initialize the    */
+/*  pointer to NULL, which tells Triangle to allocate space for the results. */
+/*  The latter option is preferable, because Triangle always knows exactly   */
+/*  how much space to allocate.  The former option is provided mainly for    */
+/*  people who need to call Triangle from Fortran code, though it also makes */
+/*  possible some nasty space-saving tricks, like writing the output to the  */
+/*  same arrays as the input.                                                */
+/*                                                                           */
+/*  Triangle will not free() any input or output arrays, including those it  */
+/*  allocates itself; that's up to you.  You should free arrays allocated by */
+/*  Triangle by calling the trifree() procedure defined below.  (By default, */
+/*  trifree() just calls the standard free() library procedure, but          */
+/*  applications that call triangulate() may replace trimalloc() and         */
+/*  trifree() in triangle.c to use specialized memory allocators.)           */
+/*                                                                           */
+/*  Here's a guide to help you decide which fields you must initialize       */
+/*  before you call triangulate().                                           */
+/*                                                                           */
+/*  `in':                                                                    */
+/*                                                                           */
+/*    - `pointlist' must always point to a list of points; `numberofpoints'  */
+/*      and `numberofpointattributes' must be properly set.                  */
+/*      `pointmarkerlist' must either be set to NULL (in which case all      */
+/*      markers default to zero), or must point to a list of markers.  If    */
+/*      `numberofpointattributes' is not zero, `pointattributelist' must     */
+/*      point to a list of point attributes.                                 */
+/*    - If the `r' switch is used, `trianglelist' must point to a list of    */
+/*      triangles, and `numberoftriangles', `numberofcorners', and           */
+/*      `numberoftriangleattributes' must be properly set.  If               */
+/*      `numberoftriangleattributes' is not zero, `triangleattributelist'    */
+/*      must point to a list of triangle attributes.  If the `a' switch is   */
+/*      used (with no number following), `trianglearealist' must point to a  */
+/*      list of triangle area constraints.  `neighborlist' may be ignored.   */
+/*    - If the `p' switch is used, `segmentlist' must point to a list of     */
+/*      segments, `numberofsegments' must be properly set, and               */
+/*      `segmentmarkerlist' must either be set to NULL (in which case all    */
+/*      markers default to zero), or must point to a list of markers.        */
+/*    - If the `p' switch is used without the `r' switch, then               */
+/*      `numberofholes' and `numberofregions' must be properly set.  If      */
+/*      `numberofholes' is not zero, `holelist' must point to a list of      */
+/*      holes.  If `numberofregions' is not zero, `regionlist' must point to */
+/*      a list of region constraints.                                        */
+/*    - If the `p' switch is used, `holelist', `numberofholes',              */
+/*      `regionlist', and `numberofregions' is copied to `out'.  (You can    */
+/*      nonetheless get away with not initializing them if the `r' switch is */
+/*      used.)                                                               */
+/*    - `edgelist', `edgemarkerlist', `normlist', and `numberofedges' may be */
+/*      ignored.                                                             */
+/*                                                                           */
+/*  `out':                                                                   */
+/*                                                                           */
+/*    - `pointlist' must be initialized (NULL or pointing to memory) unless  */
+/*      the `N' switch is used.  `pointmarkerlist' must be initialized       */
+/*      unless the `N' or `B' switch is used.  If `N' is not used and        */
+/*      `in->numberofpointattributes' is not zero, `pointattributelist' must */
+/*      be initialized.                                                      */
+/*    - `trianglelist' must be initialized unless the `E' switch is used.    */
+/*      `neighborlist' must be initialized if the `n' switch is used.  If    */
+/*      the `E' switch is not used and (`in->numberofelementattributes' is   */
+/*      not zero or the `A' switch is used), `elementattributelist' must be  */
+/*      initialized.  `trianglearealist' may be ignored.                     */
+/*    - `segmentlist' must be initialized if the `p' or `c' switch is used,  */
+/*      and the `P' switch is not used.  `segmentmarkerlist' must also be    */
+/*      initialized under these circumstances unless the `B' switch is used. */
+/*    - `edgelist' must be initialized if the `e' switch is used.            */
+/*      `edgemarkerlist' must be initialized if the `e' switch is used and   */
+/*      the `B' switch is not.                                               */
+/*    - `holelist', `regionlist', `normlist', and all scalars may be ignored.*/
+/*                                                                           */
+/*  `vorout' (only needed if `v' switch is used):                            */
+/*                                                                           */
+/*    - `pointlist' must be initialized.  If `in->numberofpointattributes'   */
+/*      is not zero, `pointattributelist' must be initialized.               */
+/*      `pointmarkerlist' may be ignored.                                    */
+/*    - `edgelist' and `normlist' must both be initialized.                  */
+/*      `edgemarkerlist' may be ignored.                                     */
+/*    - Everything else may be ignored.                                      */
+/*                                                                           */
+/*  After a call to triangulate(), the valid fields of `out' and `vorout'    */
+/*  will depend, in an obvious way, on the choice of switches used.  Note    */
+/*  that when the `p' switch is used, the pointers `holelist' and            */
+/*  `regionlist' are copied from `in' to `out', but no new space is          */
+/*  allocated; be careful that you don't free() the same array twice.  On    */
+/*  the other hand, Triangle will never copy the `pointlist' pointer (or any */
+/*  others); new space is allocated for `out->pointlist', or if the `N'      */
+/*  switch is used, `out->pointlist' remains uninitialized.                  */
+/*                                                                           */
+/*  All of the meaningful `numberof' fields will be properly set; for        */
+/*  instance, `numberofedges' will represent the number of edges in the      */
+/*  triangulation whether or not the edges were written.  If segments are    */
+/*  not used, `numberofsegments' will indicate the number of boundary edges. */
+/*                                                                           */
+/*****************************************************************************/
+
+struct triangulateio {
+  REAL *pointlist;                                               /* In / out */
+  REAL *pointattributelist;                                      /* In / out */
+  int *pointmarkerlist;                                          /* In / out */
+  int numberofpoints;                                            /* In / out */
+  int numberofpointattributes;                                   /* In / out */
+
+  int *trianglelist;                                             /* In / out */
+  REAL *triangleattributelist;                                   /* In / out */
+  REAL *trianglearealist;                                         /* In only */
+  int *neighborlist;                                             /* Out only */
+  int numberoftriangles;                                         /* In / out */
+  int numberofcorners;                                           /* In / out */
+  int numberoftriangleattributes;                                /* In / out */
+
+  int *segmentlist;                                              /* In / out */
+  int *segmentmarkerlist;                                        /* In / out */
+  int numberofsegments;                                          /* In / out */
+
+  REAL *holelist;                        /* In / pointer to array copied out */
+  int numberofholes;                                      /* In / copied out */
+
+  REAL *regionlist;                      /* In / pointer to array copied out */
+  int numberofregions;                                    /* In / copied out */
+
+  int *edgelist;                                                 /* Out only */
+  int *edgemarkerlist;            /* Not used with Voronoi diagram; out only */
+  REAL *normlist;                /* Used only with Voronoi diagram; out only */
+  int numberofedges;                                             /* Out only */
+};
+
+#ifdef ANSI_DECLARATORS
+void triangulate(char *, struct triangulateio *, struct triangulateio *,
+                 struct triangulateio *);
+void trifree(VOID *memptr);
+#else /* not ANSI_DECLARATORS */
+void triangulate();
+void trifree();
+#endif /* not ANSI_DECLARATORS */
+
+