Virtual milling additivo

EgtGeomKernel.vcxproj

@@ -394,6 +394,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" />
@@ -557,6 +559,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" />
@@ -604,6 +607,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

@@ -453,6 +453,12 @@
     <ClCompile Include="CDeRectPrismoidTria.cpp">
       <Filter>File di origine\GeoCollision</Filter>
     </ClCompile>
+    <ClCompile Include="tpp_assert.cpp">
+      <Filter>File di origine</Filter>
+    </ClCompile>
+    <ClCompile Include="tpp_impl.cpp">
+      <Filter>File di origine</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="stdafx.h">
@@ -1076,6 +1082,21 @@
     <ClInclude Include="CDeRectPrismoidTria.h">
       <Filter>File di intestazione</Filter>
     </ClInclude>
+    <ClInclude Include="dpoint.hpp">
+      <Filter>File di intestazione</Filter>
+    </ClInclude>
+    <ClInclude Include="tpp_assert.hpp">
+      <Filter>File di intestazione</Filter>
+    </ClInclude>
+    <ClInclude Include="tpp_interface.hpp">
+      <Filter>File di intestazione</Filter>
+    </ClInclude>
+    <ClInclude Include="triangle.h">
+      <Filter>File di intestazione</Filter>
+    </ClInclude>
+    <ClInclude Include="triangle_impl.hpp">
+      <Filter>File di intestazione</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="EgtGeomKernel.rc">

GdbExecutor.cpp

@@ -2775,7 +2775,448 @@ GdbExecutor::ExecuteVolZmap( const string& sCmd2, const STRVECTOR& vsParams)
    //else if ( sCmd2 == "EDGES") {
    //   return ExecuteVolZmapEdges( vsParams) ;
    //}
-   return false ;
+   ///////////////////////////if ( sCmd2 == "ELLCYL") {
+   ///////////////////////////   // Parametri: nParentId, nLineId1, nLineId2, ptLine, vtLine, ptCirc, vtAx, dRad, vtSweptVec, nInOut, ptPP1, vtPV1, ptPP2, vtPV2, ptPP3, vtPV3
+   ///////////////////////////   if (vsParams.size() != 16)
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero il riferimento in cui è immerso
+   ///////////////////////////   Frame3d frRef;
+   ///////////////////////////   if (!m_pGDB->GetGroupGlobFrame(GetIdParam(vsParams[0]), frRef))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero punto iniziale retta
+   ///////////////////////////   Point3d ptLine;
+   ///////////////////////////   if (!GetPointParam(vsParams[3], frRef, ptLine))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore retta
+   ///////////////////////////   Vector3d vtLine;
+   ///////////////////////////   if (!GetVectorParam(vsParams[4], frRef, vtLine))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtLine.Normalize();
+   ///////////////////////////   // recupero centro circonferenza
+   ///////////////////////////   Point3d ptCirc;
+   ///////////////////////////   if (!GetPointParam(vsParams[5], frRef, ptCirc))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore asse circonferenza
+   ///////////////////////////   Vector3d vtAx;
+   ///////////////////////////   if (!GetVectorParam(vsParams[6], frRef, vtAx))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtAx.Normalize();
+   ///////////////////////////   // recupero raggio della circonferenza
+   ///////////////////////////   double dRad;
+   ///////////////////////////   if (!FromString(vsParams[7], dRad))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore di traslazione
+   ///////////////////////////   Vector3d vtSweptVec;
+   ///////////////////////////   if (!GetVectorParam(vsParams[8], frRef, vtSweptVec))
+   ///////////////////////////      return false;
+   ///////////////////////////   // intero per recupero flag bInOut
+   ///////////////////////////   int nInOut;
+   ///////////////////////////   if (!FromString(vsParams[9], nInOut))
+   ///////////////////////////      return false;
+   ///////////////////////////   bool bInOut = nInOut == 0 ? false : true;
+   ///////////////////////////   // recupero punto piano 1
+   ///////////////////////////   Point3d ptPP1;
+   ///////////////////////////   if (!GetPointParam(vsParams[10], frRef, ptPP1))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 1
+   ///////////////////////////   Vector3d vtPV1;
+   ///////////////////////////   if (!GetVectorParam(vsParams[11], frRef, vtPV1))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV1.Normalize();
+   ///////////////////////////   // recupero punto piano 2
+   ///////////////////////////   Point3d ptPP2;
+   ///////////////////////////   if (!GetPointParam(vsParams[12], frRef, ptPP2))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 2
+   ///////////////////////////   Vector3d vtPV2;
+   ///////////////////////////   if (!GetVectorParam(vsParams[13], frRef, vtPV2))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV2.Normalize();
+   ///////////////////////////   // recupero punto piano 3
+   ///////////////////////////   Point3d ptPP3;
+   ///////////////////////////   if (!GetPointParam(vsParams[14], frRef, ptPP3))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 3
+   ///////////////////////////   Vector3d vtPV3;
+   ///////////////////////////   if (!GetVectorParam(vsParams[15], frRef, vtPV3))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV3.Normalize();
+   ///////////////////////////   // definisco i piani
+   ///////////////////////////   Plane3d plPlane1;
+   ///////////////////////////   plPlane1.Set(ptPP1, vtPV1);
+   ///////////////////////////   Plane3d plPlane2;
+   ///////////////////////////   plPlane2.Set(ptPP2, vtPV2);
+   ///////////////////////////   Plane3d plPlane3;
+   ///////////////////////////   plPlane3.Set(ptPP3, vtPV3);
+   ///////////////////////////   vector<Plane3d> vPlanesVec;
+   ///////////////////////////   /*vPlanesVec.emplace_back( plPlane1) ;
+   ///////////////////////////   vPlanesVec.emplace_back( plPlane2) ;
+   ///////////////////////////   vPlanesVec.emplace_back( plPlane3) ;*/
+   ///////////////////////////   // Intersezione
+   ///////////////////////////   VolZmap MyVol;
+   ///////////////////////////   Point3d ptInt1;
+   ///////////////////////////   Vector3d vtN1;
+   ///////////////////////////   Point3d ptInt2;
+   ///////////////////////////   Vector3d vtN2;
+   ///////////////////////////   int nSol = MyVol.IntersLineCircSweptSurfCuttedByPlanes(ptLine, vtLine, ptCirc, vtAx, dRad, vtSweptVec, bInOut, vPlanesVec,
+   ///////////////////////////      ptInt1, vtN1, ptInt2, vtN2);
+   ///////////////////////////   if (nSol == 1) {
+   ///////////////////////////      PtrOwner<CurveLine> pLine(CreateBasicCurveLine());
+   ///////////////////////////      pLine->Set(ptInt1, ptInt1 + vtN1);
+   ///////////////////////////      // inserisco  nel DB
+   ///////////////////////////      return AddGeoObj(vsParams[1], vsParams[0], Release(pLine));
+   ///////////////////////////   }
+   ///////////////////////////   else if (nSol == 2) {
+   ///////////////////////////      PtrOwner<CurveLine> pLine1(CreateBasicCurveLine());
+   ///////////////////////////      pLine1->Set(ptInt1, ptInt1 + vtN1);
+   ///////////////////////////      PtrOwner<CurveLine> pLine2(CreateBasicCurveLine());
+   ///////////////////////////      pLine2->Set(ptInt2, ptInt2 + vtN2);
+   ///////////////////////////      // inserisco  nel DB
+   ///////////////////////////      return AddGeoObj(vsParams[1], vsParams[0], Release(pLine1)) &&
+   ///////////////////////////         AddGeoObj(vsParams[2], vsParams[0], Release(pLine2));
+   ///////////////////////////   }
+   ///////////////////////////}
+   ///////////////////////////else if (sCmd2 == "STDCYL") {
+   ///////////////////////////   //  Parametri : nParentId, nLineId1, nLineId2, ptLineP, vtLineDir, ptBaseCen, vtAx, dRad, dH, nInOut, ptPP1, vtPV1, ptPP2, vtPV2, ptPP3, vtPV3 
+   ///////////////////////////   if (vsParams.size() != 16)
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero il riferimento in cui è immerso
+   ///////////////////////////   Frame3d frRef;
+   ///////////////////////////   if (!m_pGDB->GetGroupGlobFrame(GetIdParam(vsParams[0]), frRef))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero punto iniziale retta
+   ///////////////////////////   Point3d ptLineP;
+   ///////////////////////////   if (!GetPointParam(vsParams[3], frRef, ptLineP))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore retta
+   ///////////////////////////   Vector3d vtLineDir;
+   ///////////////////////////   if (!GetVectorParam(vsParams[4], frRef, vtLineDir))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtLineDir.Normalize();
+   ///////////////////////////   // recupero centro circonferenza
+   ///////////////////////////   Point3d ptBaseCen;
+   ///////////////////////////   if (!GetPointParam(vsParams[5], frRef, ptBaseCen))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore asse circonferenza
+   ///////////////////////////   Vector3d vtAx;
+   ///////////////////////////   if (!GetVectorParam(vsParams[6], frRef, vtAx))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtAx.Normalize();
+   ///////////////////////////   // recupero raggio della circonferenza
+   ///////////////////////////   double dRad;
+   ///////////////////////////   if (!FromString(vsParams[7], dRad))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero raggio della circonferenza
+   ///////////////////////////   double dH;
+   ///////////////////////////   if (!FromString(vsParams[8], dH))
+   ///////////////////////////      return false;
+   ///////////////////////////   // intero per recupero flag bInOut
+   ///////////////////////////   int nInOut;
+   ///////////////////////////   if (!FromString(vsParams[9], nInOut))
+   ///////////////////////////      return false;
+   ///////////////////////////   bool bInOut = nInOut == 0 ? false : true;
+   ///////////////////////////   // recupero punto piano 1
+   ///////////////////////////   Point3d ptPP1;
+   ///////////////////////////   if (!GetPointParam(vsParams[10], frRef, ptPP1))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 1
+   ///////////////////////////   Vector3d vtPV1;
+   ///////////////////////////   if (!GetVectorParam(vsParams[11], frRef, vtPV1))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV1.Normalize();
+   ///////////////////////////   // recupero punto piano 2
+   ///////////////////////////   Point3d ptPP2;
+   ///////////////////////////   if (!GetPointParam(vsParams[12], frRef, ptPP2))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 2
+   ///////////////////////////   Vector3d vtPV2;
+   ///////////////////////////   if (!GetVectorParam(vsParams[13], frRef, vtPV2))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV2.Normalize();
+   ///////////////////////////   // recupero punto piano 3
+   ///////////////////////////   Point3d ptPP3;
+   ///////////////////////////   if (!GetPointParam(vsParams[14], frRef, ptPP3))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 3
+   ///////////////////////////   Vector3d vtPV3;
+   ///////////////////////////   if (!GetVectorParam(vsParams[15], frRef, vtPV3))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV3.Normalize();
+   ///////////////////////////   // definisco i piani
+   ///////////////////////////   Plane3d plPlane1;
+   ///////////////////////////   plPlane1.Set(ptPP1, vtPV1);
+   ///////////////////////////   Plane3d plPlane2;
+   ///////////////////////////   plPlane2.Set(ptPP2, vtPV2);
+   ///////////////////////////   Plane3d plPlane3;
+   ///////////////////////////   plPlane3.Set(ptPP3, vtPV3);
+   ///////////////////////////   vector<Plane3d> vPlanesVec;
+   ///////////////////////////   /*vPlanesVec.emplace_back( plPlane1) ;
+   ///////////////////////////   vPlanesVec.emplace_back( plPlane2) ;
+   ///////////////////////////   vPlanesVec.emplace_back( plPlane3) ;*/
+   ///////////////////////////   // Intersezione
+   ///////////////////////////   VolZmap MyVol;
+   ///////////////////////////   Point3d ptInt1;
+   ///////////////////////////   Vector3d vtN1;
+   ///////////////////////////   Point3d ptInt2;
+   ///////////////////////////   Vector3d vtN2;
+   ///////////////////////////   int nSol = MyVol.IntersLineCylinderCuttedByPlanes(ptLineP, vtLineDir, ptBaseCen, vtAx, dRad, dH, bInOut, vPlanesVec, ptInt1, vtN1, ptInt2, vtN2);
+   ///////////////////////////   if (nSol == 1) {
+   ///////////////////////////      PtrOwner<CurveLine> pLine(CreateBasicCurveLine());
+   ///////////////////////////      pLine->Set(ptInt1, ptInt1 + vtN1);
+   ///////////////////////////      // inserisco  nel DB
+   ///////////////////////////      return AddGeoObj(vsParams[1], vsParams[0], Release(pLine));
+   ///////////////////////////   }
+   ///////////////////////////   else if (nSol == 2) {
+   ///////////////////////////      PtrOwner<CurveLine> pLine1(CreateBasicCurveLine());
+   ///////////////////////////      pLine1->Set(ptInt1, ptInt1 + vtN1);
+   ///////////////////////////      PtrOwner<CurveLine> pLine2(CreateBasicCurveLine());
+   ///////////////////////////      pLine2->Set(ptInt2, ptInt2 + vtN2);
+   ///////////////////////////      // inserisco  nel DB
+   ///////////////////////////      return AddGeoObj(vsParams[1], vsParams[0], Release(pLine1)) &&
+   ///////////////////////////         AddGeoObj(vsParams[2], vsParams[0], Release(pLine2));
+   ///////////////////////////   }
+   ///////////////////////////}
+   ///////////////////////////else if ( sCmd2 == "PLANE") {
+   ///////////////////////////   //  Parametri : nParentId, nLineId, ptLineP, vtLineDir, ptParOrig, vtSeg1, vtSeg2, nExtNorm
+   ///////////////////////////   if (vsParams.size() != 8)
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero il riferimento in cui è immerso
+   ///////////////////////////   Frame3d frRef;
+   ///////////////////////////   if (!m_pGDB->GetGroupGlobFrame(GetIdParam(vsParams[0]), frRef))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero punto iniziale retta
+   ///////////////////////////   Point3d ptLineP;
+   ///////////////////////////   if (!GetPointParam(vsParams[2], frRef, ptLineP))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore retta
+   ///////////////////////////   Vector3d vtLineDir;
+   ///////////////////////////   if (!GetVectorParam(vsParams[3], frRef, vtLineDir))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtLineDir.Normalize();
+   ///////////////////////////   // recupero origine parallelogramma
+   ///////////////////////////   Point3d ptParOrig;
+   ///////////////////////////   if (!GetPointParam(vsParams[4], frRef, ptParOrig))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore segmento 1
+   ///////////////////////////   Vector3d vtSeg1;
+   ///////////////////////////   if (!GetVectorParam(vsParams[5], frRef, vtSeg1))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore segmento 2
+   ///////////////////////////   Vector3d vtSeg2;
+   ///////////////////////////   if (!GetVectorParam(vsParams[6], frRef, vtSeg2))
+   ///////////////////////////      return false;
+   ///////////////////////////   // intero per recupero flag nExtNorm
+   ///////////////////////////   int nExtNorm;
+   ///////////////////////////   if (!FromString(vsParams[7], nExtNorm))
+   ///////////////////////////      return false;
+   ///////////////////////////   bool bExtNorm = nExtNorm == 0 ? false : true;
+   ///////////////////////////   // Intersezione
+   ///////////////////////////   VolZmap MyVol;
+   ///////////////////////////   Point3d ptInt;
+   ///////////////////////////   Vector3d vtN;
+   ///////////////////////////   int nSol = MyVol.IntersLineParallelogram( ptLineP, vtLineDir, ptParOrig, vtSeg1, vtSeg2, bExtNorm, ptInt, vtN) ;
+   ///////////////////////////   if (nSol == 1) {
+   ///////////////////////////      PtrOwner<CurveLine> pLine(CreateBasicCurveLine());
+   ///////////////////////////      pLine->Set(ptInt, ptInt + vtN);
+   ///////////////////////////      // inserisco  nel DB
+   ///////////////////////////      return AddGeoObj(vsParams[1], vsParams[0], Release(pLine));
+   ///////////////////////////   }
+   ///////////////////////////}
+   ///////////////////////////else if (sCmd2 == "CONE") {
+   ///////////////////////////   // Parametri : nParentId, nLineId1, nLineId2, ptLineP, vtLineDir, ptVert, vtAx, dRad, dH, nInOut, ptPP1, vtPV1, ptPP2, vtPV2, ptPP3, vtPV3
+   ///////////////////////////   if (vsParams.size() != 16)
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero il riferimento in cui è immerso
+   ///////////////////////////   Frame3d frRef;
+   ///////////////////////////   if (!m_pGDB->GetGroupGlobFrame(GetIdParam(vsParams[0]), frRef))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero punto iniziale retta
+   ///////////////////////////   Point3d ptLineP;
+   ///////////////////////////   if (!GetPointParam(vsParams[3], frRef, ptLineP))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore retta
+   ///////////////////////////   Vector3d vtLineDir;
+   ///////////////////////////   if (!GetVectorParam(vsParams[4], frRef, vtLineDir))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtLineDir.Normalize();
+   ///////////////////////////   // recupero vertice cono
+   ///////////////////////////   Point3d ptVert;
+   ///////////////////////////   if (!GetPointParam(vsParams[5], frRef, ptVert))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore asse cono
+   ///////////////////////////   Vector3d vtAx;
+   ///////////////////////////   if (!GetVectorParam(vsParams[6], frRef, vtAx))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtAx.Normalize();
+   ///////////////////////////   // recupero raggio della base del cono
+   ///////////////////////////   double dRad;
+   ///////////////////////////   if (!FromString(vsParams[7], dRad))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero altezza del cono
+   ///////////////////////////   double dH;
+   ///////////////////////////   if (!FromString(vsParams[8], dH))
+   ///////////////////////////      return false;
+   ///////////////////////////   // intero per recupero flag nExtNorm
+   ///////////////////////////   int nInOut;
+   ///////////////////////////   if (!FromString(vsParams[9], nInOut))
+   ///////////////////////////      return false;
+   ///////////////////////////   bool bInOut = nInOut == 0 ? false : true;
+   ///////////////////////////   // recupero punto piano 1
+   ///////////////////////////   Point3d ptPP1;
+   ///////////////////////////   if (!GetPointParam(vsParams[10], frRef, ptPP1))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 1
+   ///////////////////////////   Vector3d vtPV1;
+   ///////////////////////////   if (!GetVectorParam(vsParams[11], frRef, vtPV1))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV1.Normalize();
+   ///////////////////////////   // recupero punto piano 2
+   ///////////////////////////   Point3d ptPP2;
+   ///////////////////////////   if (!GetPointParam(vsParams[12], frRef, ptPP2))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 2
+   ///////////////////////////   Vector3d vtPV2;
+   ///////////////////////////   if (!GetVectorParam(vsParams[13], frRef, vtPV2))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV2.Normalize();
+   ///////////////////////////   // recupero punto piano 3
+   ///////////////////////////   Point3d ptPP3;
+   ///////////////////////////   if (!GetPointParam(vsParams[14], frRef, ptPP3))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 3
+   ///////////////////////////   Vector3d vtPV3;
+   ///////////////////////////   if (!GetVectorParam(vsParams[15], frRef, vtPV3))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV3.Normalize();
+   ///////////////////////////   // definisco i piani
+   ///////////////////////////   Plane3d plPlane1;
+   ///////////////////////////   plPlane1.Set(ptPP1, vtPV1);
+   ///////////////////////////   Plane3d plPlane2;
+   ///////////////////////////   plPlane2.Set(ptPP2, vtPV2);
+   ///////////////////////////   Plane3d plPlane3;
+   ///////////////////////////   plPlane3.Set(ptPP3, vtPV3);
+   ///////////////////////////   vector<Plane3d> vPlanesVec;
+   ///////////////////////////   /*vPlanesVec.emplace_back( plPlane1) ;
+   ///////////////////////////   vPlanesVec.emplace_back( plPlane2) ;
+   ///////////////////////////   vPlanesVec.emplace_back( plPlane3) ;*/
+   ///////////////////////////   // Intersezione
+   ///////////////////////////   VolZmap MyVol;
+   ///////////////////////////   Point3d ptInt1;
+   ///////////////////////////   Vector3d vtN1;
+   ///////////////////////////   Point3d ptInt2;
+   ///////////////////////////   Vector3d vtN2;
+   ///////////////////////////   int nSol = MyVol.IntersLineConeCuttedByPlanes( ptLineP, vtLineDir, ptVert, vtAx, dRad, dH, bInOut, vPlanesVec, ptInt1, vtN1, ptInt2, vtN2);
+   ///////////////////////////   if (nSol == 1) {
+   ///////////////////////////      PtrOwner<CurveLine> pLine(CreateBasicCurveLine());
+   ///////////////////////////      pLine->Set(ptInt1, ptInt1 + vtN1);
+   ///////////////////////////      // inserisco  nel DB
+   ///////////////////////////      return AddGeoObj(vsParams[1], vsParams[0], Release(pLine));
+   ///////////////////////////   }
+   ///////////////////////////   else if (nSol == 2) {
+   ///////////////////////////      PtrOwner<CurveLine> pLine1(CreateBasicCurveLine());
+   ///////////////////////////      pLine1->Set(ptInt1, ptInt1 + vtN1);
+   ///////////////////////////      PtrOwner<CurveLine> pLine2(CreateBasicCurveLine());
+   ///////////////////////////      pLine2->Set(ptInt2, ptInt2 + vtN2);
+   ///////////////////////////      // inserisco  nel DB
+   ///////////////////////////      return AddGeoObj(vsParams[1], vsParams[0], Release(pLine1)) &&
+   ///////////////////////////         AddGeoObj(vsParams[2], vsParams[0], Release(pLine2));
+   ///////////////////////////   }
+   ///////////////////////////}
+   ///////////////////////////else if (sCmd2 == "SPHERE") {
+   ///////////////////////////   // Parametri: nParentId, nLineId1, nLineId2, ptLineP, vtLineD, ptCen, dRad, nInOut, ptPP1, vtPV1, ptPP2, vtPV2, ptPP3, vtPV3
+   ///////////////////////////   if (vsParams.size() != 14)
+   ///////////////////////////         return false;
+   ///////////////////////////   // recupero il riferimento in cui è immerso
+   ///////////////////////////   Frame3d frRef;
+   ///////////////////////////   if (!m_pGDB->GetGroupGlobFrame(GetIdParam(vsParams[0]), frRef))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero punto iniziale retta
+   ///////////////////////////   Point3d ptLineP;
+   ///////////////////////////   if (!GetPointParam(vsParams[3], frRef, ptLineP))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore retta
+   ///////////////////////////   Vector3d vtLineD;
+   ///////////////////////////   if (!GetVectorParam(vsParams[4], frRef, vtLineD))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtLineD.Normalize();
+   ///////////////////////////   // recupero centro sfera
+   ///////////////////////////   Point3d ptCen;
+   ///////////////////////////   if (!GetPointParam(vsParams[5], frRef, ptCen))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero raggio della sfera
+   ///////////////////////////   double dRad;
+   ///////////////////////////   if (!FromString(vsParams[6], dRad))
+   ///////////////////////////      return false;
+   ///////////////////////////   // intero per recupero flag nExtNorm
+   ///////////////////////////   int nInOut;
+   ///////////////////////////   if (!FromString(vsParams[7], nInOut))
+   ///////////////////////////      return false;
+   ///////////////////////////   bool bInOut = nInOut == 0 ? false : true;
+   ///////////////////////////   // recupero punto piano 1
+   ///////////////////////////   Point3d ptPP1;
+   ///////////////////////////   if (!GetPointParam(vsParams[8], frRef, ptPP1))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 1
+   ///////////////////////////   Vector3d vtPV1;
+   ///////////////////////////   if (!GetVectorParam(vsParams[9], frRef, vtPV1))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV1.Normalize();
+   ///////////////////////////   // recupero punto piano 2
+   ///////////////////////////   Point3d ptPP2;
+   ///////////////////////////   if (!GetPointParam(vsParams[10], frRef, ptPP2))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 2
+   ///////////////////////////   Vector3d vtPV2;
+   ///////////////////////////   if (!GetVectorParam(vsParams[11], frRef, vtPV2))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV2.Normalize();
+   ///////////////////////////   // recupero punto piano 3
+   ///////////////////////////   Point3d ptPP3;
+   ///////////////////////////   if (!GetPointParam(vsParams[12], frRef, ptPP3))
+   ///////////////////////////      return false;
+   ///////////////////////////   // recupero vettore piano 3
+   ///////////////////////////   Vector3d vtPV3;
+   ///////////////////////////   if (!GetVectorParam(vsParams[13], frRef, vtPV3))
+   ///////////////////////////      return false;
+   ///////////////////////////   vtPV3.Normalize();
+   ///////////////////////////   // definisco i piani
+   ///////////////////////////   Plane3d plPlane1;
+   ///////////////////////////   plPlane1.Set(ptPP1, vtPV1);
+   ///////////////////////////   Plane3d plPlane2;
+   ///////////////////////////   plPlane2.Set(ptPP2, vtPV2);
+   ///////////////////////////   Plane3d plPlane3;
+   ///////////////////////////   plPlane3.Set(ptPP3, vtPV3);
+   ///////////////////////////   vector<Plane3d> vPlanesVec;
+   ///////////////////////////   /*vPlanesVec.emplace_back( plPlane1) ;
+   ///////////////////////////   vPlanesVec.emplace_back( plPlane2) ;
+   ///////////////////////////   vPlanesVec.emplace_back( plPlane3) ;*/
+   ///////////////////////////   // Intersezione
+   ///////////////////////////   VolZmap MyVol;
+   ///////////////////////////   Point3d ptInt1;
+   ///////////////////////////   Vector3d vtN1;
+   ///////////////////////////   Point3d ptInt2;
+   ///////////////////////////   Vector3d vtN2;
+   ///////////////////////////   int nSol = MyVol.IntersLineSphereCuttedByPlanes(ptLineP, vtLineD, ptCen, dRad, bInOut, vPlanesVec, ptInt1, vtN1, ptInt2, vtN2);
+   ///////////////////////////   if (nSol == 1) {
+   ///////////////////////////      PtrOwner<CurveLine> pLine(CreateBasicCurveLine());
+   ///////////////////////////      pLine->Set(ptInt1, ptInt1 + vtN1);
+   ///////////////////////////      // inserisco  nel DB
+   ///////////////////////////      return AddGeoObj(vsParams[1], vsParams[0], Release(pLine));
+   ///////////////////////////   }
+   ///////////////////////////   else if (nSol == 2) {
+   ///////////////////////////      PtrOwner<CurveLine> pLine1(CreateBasicCurveLine());
+   ///////////////////////////      pLine1->Set(ptInt1, ptInt1 + vtN1);
+   ///////////////////////////      PtrOwner<CurveLine> pLine2(CreateBasicCurveLine());
+   ///////////////////////////      pLine2->Set(ptInt2, ptInt2 + vtN2);
+   ///////////////////////////      // inserisco  nel DB
+   ///////////////////////////      return AddGeoObj(vsParams[1], vsParams[0], Release(pLine1)) &&
+   ///////////////////////////         AddGeoObj(vsParams[2], vsParams[0], Release(pLine2));
+   ///////////////////////////   }
+   ///////////////////////////}
+   ///////////////////////////else
+   ///////////////////////////   return false ;
+   ///////////////////////////return true ;
+   return true;
 }
 /*
 //----------------------------------------------------------------------------

GeoConst.h

@@ -49,7 +49,7 @@ const double BEZARC_ANG_CEN_MAX = 90 ;
 
 //----------------- Costanti per superfici TriMesh ---------------------------
 // tolleranza lineare standard
-const double STM_STD_LIN_TOL = 0.1 ;
+const double STM_STD_LIN_TOL = EPS_SMALL/*0.1*/ ;
 // angolo limite per definire un edge che è contorno di poligono
 const double STM_STD_BOUNDARY_ANG = 0.1 ;
 // angolo limite per mediare le normali in un vertice

PolyLine.cpp

@@ -19,9 +19,11 @@
 #include "PolygonPlane.h"
 #include "PointsPCA.h"
 #include "GeoConst.h"
+#include "/EgtDev/Include/EGkPolygon3d.h"
 #include "/EgtDev/Include/EGkPolyLine.h"
 #include "/EgtDev/Include/EGkPlane3d.h"
 #include "/EgtDev/Include/EGnStringUtils.h"
+#include "/EgtDev/Include/EGtNumUtils.h"
 
 using namespace std ;
 
@@ -53,7 +55,7 @@ PolyLine::AddUPoint( double dPar, const Point3d& ptP, bool bEndOrStart)
 {
   // se da aggiungere in coda
    if ( bEndOrStart) {
-     // se il punto è uguale all'ultimo (ignoro parametro), non lo inserisco ma ok
+     // se il punto � uguale all'ultimo (ignoro parametro), non lo inserisco ma ok
       if ( m_lUPoints.size() > 0  &&  AreSamePointApprox( ptP, m_lUPoints.back().first)) {
          ++ m_nRejected ;
          return true ;
@@ -68,7 +70,7 @@ PolyLine::AddUPoint( double dPar, const Point3d& ptP, bool bEndOrStart)
    }
   // altrimenti si aggiunge in testa
    else {
-     // se il punto è uguale al primo (ignoro parametro), non lo inserisco ma ok
+     // se il punto � uguale al primo (ignoro parametro), non lo inserisco ma ok
       if ( m_lUPoints.size() > 0  &&  AreSamePointApprox( ptP, m_lUPoints.front().first)) {
          ++ m_nRejected ;
          return true ;
@@ -92,7 +94,7 @@ PolyLine::Close( void)
   // ci devono essere almeno 2 punti
    if ( m_lUPoints.size() < 2)
       return false ;
-  // verifico non sia già chiuso
+  // verifico non sia gi� chiuso
    if ( AreSamePointApprox( m_lUPoints.front().first, m_lUPoints.back().first))
       return false ;
   // aggiungo un punto uguale al primo in coda
@@ -227,7 +229,7 @@ PolyLine::ToLoc( const Frame3d& frRef)
 bool
 PolyLine::LocToLoc( const Frame3d& frOri, const Frame3d& frDest)
 {
-  // se i due riferimenti coincidono, non devo fare alcunché
+  // se i due riferimenti coincidono, non devo fare alcunch�
    if ( AreSameFrame( frOri, frDest))
       return true ;
   // ciclo sui punti
@@ -242,7 +244,7 @@ PolyLine::LocToLoc( const Frame3d& frOri, const Frame3d& frDest)
 bool
 PolyLine::Join( PolyLine& PL, double dOffsetPar)
 {
-  // se l'altra polilinea non contiene alcunchè, esco con ok
+  // se l'altra polilinea non contiene alcunch�, esco con ok
    if ( PL.m_lUPoints.size() == 0)
       return true ;
   // verifico che l'ultimo punto di questa polilinea coincida con il primo dell'altra
@@ -394,7 +396,7 @@ PolyLine::GetPrevUPoint( double* pdPar, Point3d* pptP, bool bNotFirst) const
 bool
 PolyLine::GetCurrUPoint( double* pdPar, Point3d* pptP) const
 {
-  // verifico validità punto corrente
+  // verifico validit� punto corrente
    if ( m_iter == m_lUPoints.end())
       return false ;
 
@@ -435,7 +437,7 @@ PolyLine::GetFirstULine( double* pdIni, Point3d* pptIni, double* pdFin, Point3d*
 bool
 PolyLine::GetNextULine( double* pdIni, Point3d* pptIni, double* pdFin, Point3d* pptFin) const
 {
-  // parametro e punto iniziali (è il precedente finale)
+  // parametro e punto iniziali (� il precedente finale)
    if ( m_iter == m_lUPoints.end())
       return false ;
    if ( pdIni != nullptr)
@@ -516,19 +518,19 @@ PolyLine::IsFlat( int& nRank, Point3d& ptCen, Vector3d& vtDir, double dToler) co
    PointsPCA ptsPCA ;
    for ( bool bFound = GetFirstLine( ptP1, ptP2) ; bFound ; bFound = GetNextLine( ptP1, ptP2))
       ptsPCA.AddPoint( Media( ptP1, ptP2), Dist( ptP1, ptP2)) ;
-  // recupero il rango, ovvero la dimensionalità dell'insieme di punti
+  // recupero il rango, ovvero la dimensionalit� dell'insieme di punti
    nRank = ptsPCA.GetRank() ;
   // se dimensione nulla, o non ci sono punti o sono tutti praticamente coincidenti
    if ( nRank == 0)
       return ptsPCA.GetCenter( ptCen) ;
   // se dimensione 1, allora i punti sono distribuiti su una linea
    if ( nRank == 1) {
-     // assegno il centro e la direzione della linea (il verso è indifferente)
+     // assegno il centro e la direzione della linea (il verso � indifferente)
       ptsPCA.GetCenter( ptCen) ;
       ptsPCA.GetPrincipalComponent( 0, vtDir) ;
       return true ;
    }
-  // altrimenti dimensione 2 o 3, allora è determinato un piano principale, verifico se tutti i punti vi giacciono
+  // altrimenti dimensione 2 o 3, allora � determinato un piano principale, verifico se tutti i punti vi giacciono
    // Center and normal vector
    ptsPCA.GetCenter( ptCen) ;
    Vector3d vtX, vtY ;
@@ -536,9 +538,9 @@ PolyLine::IsFlat( int& nRank, Point3d& ptCen, Vector3d& vtDir, double dToler) co
    ptsPCA.GetPrincipalComponent( 1, vtY) ;
    vtDir = vtX ^ vtY ;
    if ( ! vtDir.Normalize()) {
-     // riduco la dimensionalità a lineare
+     // riduco la dimensionalit� a lineare
       nRank = 1 ;
-     // assegno il centro e la direzione della linea (il verso è indifferente)
+     // assegno il centro e la direzione della linea (il verso � indifferente)
       ptsPCA.GetCenter( ptCen) ;
       vtDir = vtX ;
       return true ;
@@ -567,12 +569,12 @@ PolyLine::IsFlat( Plane3d& plPlane, double dToler) const
       plPlane.Reset() ;
       return false ;
    }
-  // recupero dati sulla planarità della polilinea
+  // recupero dati sulla planarit� della polilinea
    int nRank ;
    Point3d ptCen ;
    Vector3d vtDir ;
    bool bFlat = IsFlat( nRank, ptCen, vtDir, dToler) ;
-   // imposto il piano a seconda della dimensionalità
+   // imposto il piano a seconda della dimensionalit�
    switch ( nRank) {
    case 0 :      // punto
       plPlane.Set( ptCen, Z_AX) ;
@@ -643,7 +645,7 @@ PolyLine::GetAreaXY( double& dArea) const
   // verifico sia chiusa
    if ( ! IsClosed())
       return false ;
-  // calcolo l'area considerando solo XY (è la Z di Newell)
+  // calcolo l'area considerando solo XY (� la Z di Newell)
    dArea = 0 ;
    Point3d ptIni, ptFin ;
    for ( bool bFound = GetFirstLine( ptIni, ptFin) ; bFound ; bFound = GetNextLine( ptIni, ptFin)) {
@@ -777,10 +779,10 @@ PolyLine::RemoveAlignedPoints( double dToler)
    }
   // se curva chiusa con almeno 4 punti, devo analizzare il terzetto attorno alla chiusura
    if ( IsClosed()  &&  m_lUPoints.size() >= 4) {
-     // precP e currP sono già corretti
+     // precP e currP sono gi� corretti
      // il primo punto ripete l'ultimo (geometricamente coincide con currP)
       auto firstP = m_lUPoints.begin() ;
-     // questo è il vero successivo
+     // questo � il vero successivo
       nextP = next( firstP) ;
      // distanza del punto corrente dal segmento che unisce gli adiacenti
       DistPointLine dPL( currP->first, precP->first, nextP->first) ;
@@ -831,7 +833,7 @@ PolyLine::MyChangeStart( int nPos)
       return false ;
   // cancello ultimo punto ( coincide con primo)
    m_lUPoints.pop_back() ;
-  // sposto la metà iniziale dei punti alla fine
+  // sposto la met� iniziale dei punti alla fine
    for ( int i = 0 ; i < nPos ; ++ i)
       m_lUPoints.splice( m_lUPoints.end(), m_lUPoints, m_lUPoints.begin()) ;
   // aggiungo punto finale come copia dell'iniziale
@@ -936,7 +938,7 @@ PolyLine::MyApproxOnSide( const Vector3d& vtN,  bool bLeftSide, double dToler)
             }
          }
       }
-     // non è stato eliminato alcunché
+     // non � stato eliminato alcunch�
      // ripristino la tolleranza corrente
       dCurrToler = dToler ;
      // avanzo il terzetto di uno step
@@ -977,7 +979,7 @@ PolyLine::MakeConvex( const Vector3d& vtN,  bool bLeftSide)
 bool
 PolyLine::MyMakeConvex( const Vector3d& vtN,  bool bLeftSide)
 {
-  // ciclo i controlli finchè non ci sono rimozioni
+  // ciclo i controlli finch� non ci sono rimozioni
    bool bRemoved = true ;
    while ( bRemoved) {
       bRemoved = false ;
@@ -1005,7 +1007,7 @@ PolyLine::MyMakeConvex( const Vector3d& vtN,  bool bLeftSide)
             bRemoved = true ;
             continue ;
          }
-        // non è stato eliminato alcunché : avanzo il terzetto di uno step
+        // non � stato eliminato alcunch� : avanzo il terzetto di uno step
          precP = currP ;
          currP = nextP ;
          ++ nextP ;
@@ -1032,7 +1034,7 @@ PolyLine::Invert( bool bInvertU)
    m_lUPoints.reverse() ;
   // se richiesto, inverto anche il parametro U
    if ( bInvertU) {
-     // recupero il primo valore di U che è il vecchio finale ed è il riferimento di inversione
+     // recupero il primo valore di U che � il vecchio finale ed � il riferimento di inversione
       double dUfin = m_lUPoints.front().second ;
      // ciclo su tutti gli elementi
       for ( auto& UPoint : m_lUPoints) {
@@ -1247,7 +1249,7 @@ PolyLine::GetMinAreaRectangleXY( Point3d& ptCen, Vector3d& vtAx, double& dLen, d
 bool
 PolyLine::Trim( const Plane3d& plPlane, bool bInVsOut)
 {
-  // se vuota non faccio alcunché
+  // se vuota non faccio alcunch�
    if ( m_lUPoints.size() == 0)
       return false ;
 
@@ -1296,3 +1298,329 @@ PolyLine::Trim( const Plane3d& plPlane, bool bInVsOut)
 
    return true ;
 }
+
+//----------------------------------------------------------------------------
+bool
+PolyLine::ChangePolyLineStart( const Point3d& ptNewStart, double dTol)
+{
+   // Rinomino la lista di punti della PolyLine.
+   PNTULIST& LoopList = 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
+bool
+PolyLine::PointPositionOnPolyLine( const Point3d& ptPoint, int& nSegNum, double& dParOnSeg, double dTol) const
+{
+   // Rinomino la lista di punti della PolyLine.
+   const PNTULIST& LoopList = 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 ;
+}
+
+//----------------------------------------------------------------------------
+bool
+PolyLine::IsPointInsidePolyLine( const Point3d& ptP) const
+{
+   // Se la PolyLine non � chiusa, il punto non pu� essere interno.
+   if ( ! IsClosed())
+      return false ;
+   // Lista dei punti
+   const PNTULIST& List = 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( *this) ;
+      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( *this) ;
+      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( *this) ;
+      Vector3d vtPolyNorm = AuxPolygon.GetVersN() ;
+      Vector3d vtOut = vtTan ^ vtPolyNorm ;
+      vtP -= ( vtP * vtTan) * vtTan ;
+      if ( vtP * vtOut < - EPS_SMALL)
+         return true ;
+   }
+   return false ;
+}
+
+//----------------------------------------------------------------------------
+bool
+PolyLine::DistPointPolyLine( const Point3d& ptP, double& dPointPolyLineDist) const
+{
+   if ( GetPointNbr() == 0)
+      return false ;
+   dPointPolyLineDist = DBL_MAX ;
+   Point3d ptSt, ptEn ;
+   bool bContinue = GetFirstPoint( ptSt)  &&  GetNextPoint( ptEn) ;
+   while ( bContinue) {
+      double dPoinLineDist ;
+      DistPointLine PointLineDistCalc( ptP, ptSt, ptEn) ;
+      PointLineDistCalc.GetDist( dPoinLineDist) ;
+      if ( dPoinLineDist < dPointPolyLineDist)
+         dPointPolyLineDist = dPoinLineDist ;
+      ptSt = ptEn ;
+      bContinue = GetNextPoint( ptEn) ;
+   }
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+PolyLine::SplitPolyLineAtPoint( const Point3d& ptPoint, PolyLine& Loop1, PolyLine& Loop2, double dTol) const
+{
+   // Rinomino la lista di punti della PolyLine.
+   const PNTULIST& LoopList = 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.
+   bool bPointOnEnd = false ;
+   auto itStop = itMinDistSt ;
+   auto itNext = itMinDistSt ;
+   ++ itNext ;
+   if ( AreSamePointApprox( ptPoint, itStop->first))
+      ;
+   else if ( AreSamePointApprox( ptPoint, itNext->first))
+      itStop = itNext ;
+   else {
+      bPointOnEnd = true ;
+      itStop = itNext ;
+   }
+   // 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( ptPoint, 0) ;
+   if ( bPointOnEnd) {
+      LoopList2.emplace_back( ptPoint, 0) ;
+   }
+   for ( auto it = itStop ; it != LoopList.end() ; ++ it) {
+      LoopList2.emplace_back( it->first, it->second) ;
+   }
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+PolyLine::AddPolyLineToPolyLine( PolyLine& PolyToAdd, double dTol)
+{
+   // Se la PolyLine a cui devo aggiungere l'altra � chiusa, non posso aggiungere nulla.
+   if ( 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 ;
+   GetLastPoint( ptLast) ;
+   auto it = PolyToAddList.begin() ;
+   if ( 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) {
+      AddUPoint( 0., it->first) ;
+   }
+   return true ;
+}

SurfTriMesh.cpp

@@ -41,7 +41,7 @@ GEOOBJ_REGISTER( SRF_TRIMESH, NGE_S_TRM, SurfTriMesh) ;
 SurfTriMesh::SurfTriMesh( void)
    : m_nStatus( TO_VERIFY), m_dLinTol( STM_STD_LIN_TOL), m_dBoundaryAng( STM_STD_BOUNDARY_ANG),
      m_dSmoothAng( STM_STD_SMOOTH_ANG), m_bOriented( false), m_bClosed( false), m_bFaceted( false),
-     m_nTimeStamp( 0), m_nMaxTFlag( 0), m_nParts( -1),  m_pHGrd3d( nullptr)
+     m_nTimeStamp( 0), m_nMaxTFlag( 0), m_nParts( -1),  m_pHGrd3d( nullptr), m_pPGrd3d( nullptr)
 {
    m_dCosBndAng = cos( m_dBoundaryAng * DEGTORAD) ;
    m_dCosSmAng = cos( m_dSmoothAng * DEGTORAD) ;
@@ -53,19 +53,21 @@ SurfTriMesh::SurfTriMesh( void)
 SurfTriMesh::~SurfTriMesh( void)
 {
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 }
 
 //----------------------------------------------------------------------------
 bool
 SurfTriMesh::Init( int nNumVert, int nNumTria, int nNumFacet)
 {
-  // imposto ricalcolo della grafica e di hashgrids3d
+  // imposto ricalcolo della grafica e di hashgrids3d e pointgrid3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
   // se superficie vuota
    if ( nNumVert == 0  &&  nNumTria == 0  &&  nNumFacet == 0)
       return true ;
-  // verifico validità parametri
+  // verifico validit� parametri
    if ( nNumVert < 3  ||  nNumTria < 1)
       return false ;
   // prealloco la memoria
@@ -99,6 +101,7 @@ SurfTriMesh::Clear( void)
    m_vTria.clear() ;
    m_vFacet.clear() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
    m_nMaxTFlag = 0 ;
    m_nParts = -1 ;
    return true ;
@@ -112,10 +115,19 @@ SurfTriMesh::AddVertex( const Point3d& ptVert)
    m_nStatus = TO_VERIFY ;
    m_nParts = - 1 ;
    m_OGrMgr.Reset() ;
-   ResetHashGrids3d() ;
+   // Verifico che il punto non sia distante meno di epsilon da un altro gi� esistente.
+   Point3d ptAddingPoint ;
+   VerifyPointGrid3d() ;
+   int nNearestId ;
+   if ( m_pPGrd3d->FindNearest( ptVert, EPS_SMALL, nNearestId))
+      GetVertex( nNearestId, ptAddingPoint) ;
+   else
+      ptAddingPoint = ptVert ; 
   // inserisco il vertice
    try { m_vVert.emplace_back( ptVert) ;}
    catch(...) { return SVT_NULL ;}
+   // Inserisco il vertice nella griglia
+   m_pPGrd3d->InsertPoint( ptVert, int( m_vVert.size() - 1)) ;
   // ne determino l'indice
    return int( m_vVert.size() - 1) ;
 }
@@ -182,7 +194,8 @@ SurfTriMesh::AddTriangle( const int nIdVert[3], int nTFlag)
    m_nStatus = TO_VERIFY ;
    m_nParts = - 1 ;
    m_OGrMgr.Reset() ;
-   ResetHashGrids3d() ;
+   //ResetHashGrids3d() ;
+   //ResetPointGrid3d() ;
   // inserisco il triangolo
    try { m_vTria.emplace_back( nIdVert, nTFlag) ;}
    catch(...) { return SVT_NULL ;}
@@ -200,6 +213,7 @@ SurfTriMesh::SetTriangle( int nInd, const StmTria& tT)
    m_nStatus = TO_VERIFY ;
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
   // recupero la dimensione originale
    int nPrevSize = int( m_vTria.size()) ;
   // determino la dimensione necessaria
@@ -231,14 +245,14 @@ SurfTriMesh::RemoveTriangle( int nId)
   // verifico esistenza del triangolo 
    if ( nId < 0  ||  nId >= GetTriangleSize())
       return false ;
-  // verifico se già cancellato
+  // verifico se gi� cancellato
    if ( m_vTria[nId].nIdVert[0] == SVT_DEL)
       return true ;
   // aggiorno eventuali riferimenti dei vertici
    for ( int i = 0 ; i < 3 ; ++ i) {
      // indice vertice
       int nV = m_vTria[nId].nIdVert[i] ;
-     // se vertice non c'è passo al prossimo
+     // se vertice non c'� passo al prossimo
       if ( nV < 0  ||  nV >= int( m_vTria.size()))
          continue ;
       if ( m_vVert[nV].nIdTria == nId) {
@@ -256,7 +270,7 @@ SurfTriMesh::RemoveTriangle( int nId)
    for ( int i = 0 ; i < 3 ; ++ i) {
      // indice triangolo adiacente
       int nAdjT = m_vTria[nId].nIdAdjac[i] ;
-     // se triangolo adiacente non c'è passo al prossimo
+     // se triangolo adiacente non c'� passo al prossimo
       if ( nAdjT == SVT_NULL  ||  m_vTria[nAdjT].nIdVert[0] == SVT_DEL)
          continue ;
      // ne sistemo la contro-adiacenza
@@ -269,7 +283,7 @@ SurfTriMesh::RemoveTriangle( int nId)
    m_vTria[nId].nIdVert[0] = SVT_DEL ;
   // invalido calcolo facce
    m_bFaceted = false ;
-  // invalido calcolo connettività
+  // invalido calcolo connettivit�
    m_nParts = - 1 ;
    return true ;
 }
@@ -335,7 +349,7 @@ SurfTriMesh::GetCentroid( Point3d& ptCen) const
   // la superficie deve essere validata
    if ( m_nStatus != OK)
       return false ;
-  // se la superficie è chiusa, calcolo il centroide del solido
+  // se la superficie � chiusa, calcolo il centroide del solido
    if ( IsClosed()) {
      // applico le formule di R. Nurnberg Imperial College London ad ogni faccia
       Triangle3d Tria ;
@@ -688,13 +702,13 @@ SurfTriMesh::GetTriangleBoundaryEdges( int nId, TriFlags3d& TFlags) const
    for ( int i = 0 ; i < 3 ; ++ i) {
      // indice triangolo adiacente al lato
       int nT = m_vTria[nId].nIdAdjac[i] ;
-     // se già definite le facce, verifico indice faccia
+     // se gi� definite le facce, verifico indice faccia
       if ( m_bFaceted) {
          TFlags.bFlag[i] = ( nT == SVT_NULL  ||  m_vTria[nId].nIdFacet != m_vTria[nT].nIdFacet) ;
       }
      // altrimenti verifico con le normali
       else
-        // se non c'è triangolo adiacente o se forma un angolo oltre il limite, il lato è un contorno
+        // se non c'� triangolo adiacente o se forma un angolo oltre il limite, il lato � un contorno
          TFlags.bFlag[i] = ( nT == SVT_NULL  ||  m_vTria[nId].vtN * m_vTria[nT].vtN < m_dCosBndAng) ;
    }
 
@@ -778,7 +792,7 @@ SurfTriMesh::GetTriangleSmoothNormal( int nT, int nV, Vector3d& vtN) const
    if ( nPos == -1)
       return false ;
 
-  // medio le normali, finché non incontro degli spigoli
+  // medio le normali, finch� non incontro degli spigoli
    vtN = m_vTria[nT].vtN ;
    // parto dal triangolo e vado in direzione positiva
    int nLim = nPos ;
@@ -881,7 +895,7 @@ SurfTriMesh::GetLoops( POLYLINEVECTOR& vPL) const
             if ( ! MarchAlongLoop( nT, 1, m_nTimeStamp, vPL.back()))
                return false ;
          }
-        // se il lato 0 è di contorno
+        // se il lato 0 � di contorno
          else if ( nAdjT[0] == SVT_NULL) {
            // ho trovato l'inizio di un loop
             vPL.emplace_back() ;
@@ -893,7 +907,7 @@ SurfTriMesh::GetLoops( POLYLINEVECTOR& vPL) const
             if ( ! MarchAlongLoop( nT, 1, m_nTimeStamp, vPL.back()))
                return false ;
          }
-        // se il lato 1 è di contorno
+        // se il lato 1 � di contorno
          else if ( nAdjT[1] == SVT_NULL) {
            // ho trovato l'inizio di un loop
             vPL.emplace_back() ;
@@ -905,7 +919,7 @@ SurfTriMesh::GetLoops( POLYLINEVECTOR& vPL) const
             if ( ! MarchAlongLoop( nT, 2, m_nTimeStamp, vPL.back()))
                return false ;
          }
-        // se il lato 2 è di contorno
+        // se il lato 2 � di contorno
          else if ( nAdjT[2] == SVT_NULL) {
            // ho trovato l'inizio di un loop
             vPL.emplace_back() ;
@@ -917,7 +931,7 @@ SurfTriMesh::GetLoops( POLYLINEVECTOR& vPL) const
             if ( ! MarchAlongLoop( nT, 0, m_nTimeStamp, vPL.back()))
                return false ;
          }
-        // altrimenti non c'è contorno
+        // altrimenti non c'� contorno
          else {
            // marco il triangolo come verificato
             m_vTria[nT].nTemp = m_nTimeStamp ;
@@ -962,10 +976,10 @@ SurfTriMesh::MarchOneTria( int& nT, int& nV, int nTimeStamp,
       return false ;
   // vertice di fine adiacenza e indice del successivo lato
    int nAdjV = Next( nAdjS) ;
-  // verifico se il lato successivo è un bordo
+  // verifico se il lato successivo � un bordo
    int nNextT = m_vTria[nAdjT].nIdAdjac[nAdjV] ;
    if ( nNextT == SVT_NULL) {
-     // se già recuperato
+     // se gi� recuperato
       if ( m_vTria[nAdjT].nTemp == nTimeStamp) {
          bEnd = true ;
          return true ;
@@ -1022,7 +1036,7 @@ SurfTriMesh::GetSilhouette( const Vector3d& vtDir, double dTol, POLYLINEVECTOR&
 
   // Elimino i triangoli con normale non equiversa alla direzione scelta
    for ( int i = 0 ; i < int( pStm->m_vTria.size()) ; ++ i) {
-     // se già cancellato, passo oltre
+     // se gi� cancellato, passo oltre
       if ( pStm->m_vTria[i].nIdVert[0] == SVT_DEL)
          continue ;
      // verifico la normale
@@ -1208,7 +1222,7 @@ SurfTriMesh::Dump( string& sOut, bool bMM, const char* szNewLine) const
   // segnalo eventuale incongruenza di orientamento
    if ( ! m_bOriented)
       sOut += string( "Inconsistent Orientation") + szNewLine ;
-  // segnalo numero di parti se più di una
+  // segnalo numero di parti se pi� di una
    int nParts = GetPartCount() ;
    if ( nParts > 1)
       sOut += string( "Parts =") + ToString( nParts) + szNewLine ;
@@ -1297,11 +1311,12 @@ SurfTriMesh::Save( NgeWriter& ngeOut) const
 bool
 SurfTriMesh::Load( NgeReader& ngeIn)
 {
-  // imposto ricalcolo della grafica, della connessione e di hashgrids3d
+  // imposto ricalcolo della grafica, della connessione e di hashgrids3d e pointgrid3d
    m_OGrMgr.Reset() ;
    m_nMaxTFlag = 0 ;
    m_nParts = -1 ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
   // leggo la prossima linea ( 2 parametri : dLinTol e dSmoothAng)
   // tolleranza lineare di costruzione
    double dLinTol ;
@@ -1438,7 +1453,7 @@ SurfTriMesh::Validate( bool bCorrect)
 
   // Verifico che i triangoli riferiti dalle facce esistano
    for ( int i = 0 ; i < GetFacetSize() && m_nStatus == OK && m_bFaceted ; ++ i) {
-     // verifico validità triangolo riferito 
+     // verifico validit� triangolo riferito 
       if ( m_vFacet[i] <= SVT_NULL  ||
            m_vFacet[i] >= GetTriangleSize()  ||
            m_vTria[ m_vFacet[i]].nIdVert[0] == SVT_DEL)
@@ -1593,7 +1608,7 @@ SurfTriMesh::AdjustAdjacencies( void)
 bool
 SurfTriMesh::AdjustOrientations( void)
 {
-  // se non ci sono almeno 2 triangoli è inutile fare test
+  // se non ci sono almeno 2 triangoli � inutile fare test
    if ( m_vTria.size() < 2) {
       m_bOriented = true ;
       return true ;
@@ -1637,7 +1652,7 @@ SurfTriMesh::AdjustTriaOrientation( TRINTDEQUE& S3iQ)
    S3iQ.pop_front() ;
   // assegno time stamp al triangolo
    m_vTria[nT].nTemp = m_nTimeStamp ;
-  // se c'è triangolo di riferimento, devo verificare se da invertire
+  // se c'� triangolo di riferimento, devo verificare se da invertire
    if ( nRefT != SVT_NULL) {
      // cerco indice half-edge in comune
       int nE = 0 ;
@@ -1655,7 +1670,7 @@ SurfTriMesh::AdjustTriaOrientation( TRINTDEQUE& S3iQ)
    bool bOk = true ;
    for ( int j = 0 ; j < 3 ; ++ j) {
       int nAdjT = m_vTria[nT].nIdAdjac[j] ;
-     // se non c'è adiacenza o va sul triangolo di provenienza
+     // se non c'� adiacenza o va sul triangolo di provenienza
       if ( nAdjT == SVT_NULL  ||  nAdjT == nRefT)
          ;
      // la verifico
@@ -1708,7 +1723,7 @@ SurfTriMesh::TestSealing( void)
 bool
 SurfTriMesh::AdjustTopology( void)
 {
-  // se non è rimasto alcunché di valido, pulisco tutto ed esco
+  // se non � rimasto alcunch� di valido, pulisco tutto ed esco
    if ( GetVertexCount() < 3  ||  GetTriangleCount() < 1) {
       Clear() ;
       m_bOriented = true ;
@@ -1725,7 +1740,7 @@ SurfTriMesh::AdjustTopology( void)
   // verifica adiacenze
    if ( ! AdjustAdjacencies())
       return false ;
-  // verifica continuità orientazione
+  // verifica continuit� orientazione
    if ( ! AdjustOrientations())
       return false ;
   // verifica chiusura
@@ -1763,7 +1778,7 @@ SurfTriMesh::PackVertices( void)
          vVId.push_back( SVT_DEL) ;
       }
    }
-  // se non c'è stata compattazione, esco
+  // se non c'� stata compattazione, esco
    if ( nFirstFree == SVT_NULL)
       return true ;
   // lunghezza vettore indici vertici
@@ -1778,7 +1793,7 @@ SurfTriMesh::PackVertices( void)
      // salto i triangoli cancellati
       if ( vOId[0] == SVT_DEL)
          continue ;
-     // verifico la validità degli indici
+     // verifico la validit� degli indici
       if ( vOId[0] < 0  ||  vOId[0] >= nVIdSize  ||
            vOId[1] < 0  ||  vOId[1] >= nVIdSize  ||
            vOId[2] < 0  ||  vOId[2] >= nVIdSize)
@@ -1818,11 +1833,12 @@ SurfTriMesh::PackTriangles( void)
          vTId.push_back( SVT_DEL) ;
       }
    }
-  // se non c'è stata compattazione, esco
+  // se non c'� stata compattazione, esco
    if ( nFirstFree == SVT_NULL)
       return true ;
-  // Invalido HashGrid
+  // Invalido HashGrid e Pointgrid
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
   // lunghezza vettore indici triangoli
    int nTIdSize = int( vTId.size()) ;
   // aggiorno gli indici ai triangoli dai vertici
@@ -1832,7 +1848,7 @@ SurfTriMesh::PackTriangles( void)
      // salto vertice cancellato
       if ( nOId == SVT_DEL)
          continue ;
-     // verifico la validità dell'indice
+     // verifico la validit� dell'indice
       if ( nOId < 0  ||  nOId >= nTIdSize)
          return false ;
      // aggiorno    
@@ -1862,7 +1878,7 @@ SurfTriMesh::PackTriangles( void)
      // salto le facets non valide 
       if ( m_vFacet[nId] == SVT_DEL)
          continue ;
-     // verifico validità indice a triangolo
+     // verifico validit� indice a triangolo
       if ( m_vFacet[nId] < 0  ||  m_vFacet[nId] >= nTIdSize)
          return false ;
      // aggiorno
@@ -1963,8 +1979,9 @@ SurfTriMesh::CreateByExtrusion( const PolyLine& PL, const Vector3d& vtExtr)
    m_nStatus = ERR ;
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
-  // verifico se la polilinea è chiusa
+  // verifico se la polilinea � chiusa
    bool bClosed = PL.IsClosed() ;
 
   // costruisco la mesh
@@ -1997,7 +2014,7 @@ SurfTriMesh::CreateByExtrusion( const PolyLine& PL, const Vector3d& vtExtr)
    }
   // se curva chiusa, aggiungo gli ultimi due triangoli
    if ( bClosed) {
-     // non devo aggiungere i vertici, perchè coincidono con quelli iniziali
+     // non devo aggiungere i vertici, perch� coincidono con quelli iniziali
      // aggiungo i due triangoli relativi
       nIdV[0] = nV ;
       nIdV[1] = nV - 1 ;
@@ -2016,7 +2033,7 @@ SurfTriMesh::CreateByExtrusion( const PolyLine& PL, const Vector3d& vtExtr)
 bool
 SurfTriMesh::CreateByPointCurve( const Point3d& ptP, const PolyLine& PL)
 {
-  // verifico validità punto/polilinea
+  // verifico validit� punto/polilinea
    bool bClosed = PL.IsClosed() ;
   // se chiusa, la polilinea deve avere almeno 3 punti
    if ( bClosed) {
@@ -2033,6 +2050,7 @@ SurfTriMesh::CreateByPointCurve( const Point3d& ptP, const PolyLine& PL)
    m_nStatus = ERR ;
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // costruisco la mesh
    int nVertNbr =  1 + PL.GetPointNbr() ;
@@ -2060,7 +2078,7 @@ SurfTriMesh::CreateByPointCurve( const Point3d& ptP, const PolyLine& PL)
      // aggiorno indice punto precedente su curva
       nIdV[1] = nIdV[2] ;
    }
-  // se chiusa aggiungo l'ultimo triangolo (non il vertice perchè è il primo della curva)
+  // se chiusa aggiungo l'ultimo triangolo (non il vertice perch� � il primo della curva)
    if ( bClosed) {
       nIdV[2] = 1 ;
      // inserisco il triangolo A2p -> A1p -> A1s
@@ -2080,6 +2098,7 @@ SurfTriMesh::CreateByTwoCurves( const PolyLine& PL1, const PolyLine& PL2, int nR
    m_nStatus = ERR ;
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // se rigata a minima distanza tra le due curve
    if ( nRuledType == RLT_MINDIST) {
@@ -2205,7 +2224,7 @@ SurfTriMesh::CreateByTwoCurves( const PolyLine& PL1, const PolyLine& PL2, int nR
       }
      // ciclo sui punti
       while ( bNext1  ||  bNext2) {
-        // se c'è nuovo V1s e la diagonale più corta è V2p -> V1s oppure non c'è V2s
+        // se c'� nuovo V1s e la diagonale pi� corta � V2p -> V1s oppure non c'� V2s
          if ( ! bNext2  ||  ( bNext1  &&  ( nP1s == vPnt2[nP2p].second  ||  vPnt1[nP1s].second == nP2p))) {
            // inserisco il vertice V1s (se ultimo e curve chiuse, prendo il primo)
             if ( nP1s == nTotP1 - 1  &&  bClosed)
@@ -2227,7 +2246,7 @@ SurfTriMesh::CreateByTwoCurves( const PolyLine& PL1, const PolyLine& PL2, int nR
             nV1p = nV1s ; nP1p = nP1s ; ++ nP1s ;
             bNext1 = ( nP1s < nTotP1) ;
          }
-        // altrimenti è V1p -> V2s
+        // altrimenti � V1p -> V2s
          else {
            // inserisco il vertice V2s (se ultimo e curve chiuse, prendo il primo)
             if ( nP2s == nTotP2 - 1  &&  bClosed)
@@ -2254,7 +2273,7 @@ SurfTriMesh::CreateByTwoCurves( const PolyLine& PL1, const PolyLine& PL2, int nR
    
   // altrimenti rigata con parametrizzazione sincrona sulle due curve
    else {
-     // verifico validità polilinee (devono avere almeno 2 punti e non coincidere se non agli estremi aperti)
+     // verifico validit� polilinee (devono avere almeno 2 punti e non coincidere se non agli estremi aperti)
       if ( ! VerifyPolylinesForTwoCurves( PL1, PL2))
          return false ;
 
@@ -2329,7 +2348,7 @@ SurfTriMesh::CreateByTwoCurves( const PolyLine& PL1, const PolyLine& PL2, int nR
       }
      // ciclo sui punti
       while ( bNext1  ||  bNext2) {
-        // se c'è nuovo dA1s e la diagonale più corta è dA2p -> dA1s oppure non c'è dA2s
+        // se c'� nuovo dA1s e la diagonale pi� corta � dA2p -> dA1s oppure non c'� dA2s
          if ( ( bNext1  &&  ( dA1s - dA2p) <= ( dA2s - dA1p) + EPS_PARAM)  ||  ! bNext2) {
            // inserisco il vertice A1s
             if ( ( nV1s = AddVertex( ptP1s)) == SVT_NULL)
@@ -2353,7 +2372,7 @@ SurfTriMesh::CreateByTwoCurves( const PolyLine& PL1, const PolyLine& PL2, int nR
             if ( bNext1)
                dA1s = ( dU1s - dU1F) / dDeltaU1 ;
          }
-        // altrimenti è dA1p -> dA2s
+        // altrimenti � dA1p -> dA2s
          else {
            // inserisco il vertice A2s
             if ( ( nV2s = AddVertex( ptP2s)) == SVT_NULL)
@@ -2382,7 +2401,7 @@ SurfTriMesh::CreateByTwoCurves( const PolyLine& PL1, const PolyLine& PL2, int nR
       if ( bClosed) {
          dA1s = 1 ;
          dA2s = 1 ;
-        // se la diagonale più corta è dA2p -> dA1s = 0
+        // se la diagonale pi� corta � dA2p -> dA1s = 0
          if ( ( dA1s - dA2p) <= ( dA2s - dA1p) + EPS_PARAM) {
            // inserisco il triangolo A2p -> A1p -> A1s = 0
             nIdV[0] = nV2p ;
@@ -2397,7 +2416,7 @@ SurfTriMesh::CreateByTwoCurves( const PolyLine& PL1, const PolyLine& PL2, int nR
             if ( AddTriangle( nIdV) == SVT_NULL)
                return false ;
          }
-        // altrimenti è dA1p -> dA2s = 1
+        // altrimenti � dA1p -> dA2s = 1
          else {
            // inserisco il triangolo A2p -> A1p -> A2s = 1
             nIdV[0] = nV2p ;
@@ -2469,7 +2488,7 @@ SurfTriMesh::VerifyPolylinesForTwoCurves( const PolyLine& PL1, const PolyLine& P
       return false ;
   // verifiche sui punti successivi (non sugli ultimi)
    while ( bNext1  ||  bNext2) {
-     // se c'è nuovo dA1s e la diagonale più corta è dA2p -> dA1s oppure non c'è dA2s
+     // se c'� nuovo dA1s e la diagonale pi� corta � dA2p -> dA1s oppure non c'� dA2s
       if ( ( bNext1  &&  ( dA1s - dA2p) <= ( dA2s - dA1p) + EPS_PARAM)  ||  ! bNext2) {
         // verifico se coincidono
          if ( AreSamePointApprox( ptP2p, ptP1s))
@@ -2480,7 +2499,7 @@ SurfTriMesh::VerifyPolylinesForTwoCurves( const PolyLine& PL1, const PolyLine& P
          if ( bNext1)
             dA1s = ( dU1s - dU1F) / dDeltaU1 ;
       }
-     // altrimenti è dA1p -> dA2s = 1
+     // altrimenti � dA1p -> dA2s = 1
       else {
         // verifico se coincidono
          if ( AreSamePointApprox( ptP1p, ptP2s))
@@ -2584,7 +2603,7 @@ AdjustPolylineForRevolution( PolyLine& PL, const Point3d& ptAx, const Vector3d&
    plTrim.Translate( DIST_SIC * vtTrN) ;
    PL.Trim( plTrim, false) ;
 
-  // Se polilinea risultante è aperta con estremità molto vicine all'asse le porto su questo
+  // Se polilinea risultante � aperta con estremit� molto vicine all'asse le porto su questo
    if ( ! PL.IsClosed()) {
      // verifico l'inizio
       double  dUStart ;
@@ -2661,8 +2680,9 @@ SurfTriMesh::CreateByScrewing( const PolyLine& PL, const Point3d& ptAx, const Ve
    m_nStatus = ERR ;
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
-  // verifico se la polilinea è chiusa
+  // verifico se la polilinea � chiusa
    bool bClosed = MyPL.IsClosed() ;
 
   // costruisco la mesh
@@ -2683,7 +2703,7 @@ SurfTriMesh::CreateByScrewing( const PolyLine& PL, const Point3d& ptAx, const Ve
    // verifico se il punto giace sull'asse e vi sta fisso
    bool bPrevOnAx = bOnlyRev && DistPointLine( ptP, ptAx, vtAx, 1, false).IsSmall() ;
    int  nVPrevOnAx = nV ;
-   // se non è fisso sull'asse, inserisco le copie ruotate
+   // se non � fisso sull'asse, inserisco le copie ruotate
    if ( ! bPrevOnAx) {
       for ( int i = 1 ; i <= nStep ; ++i) {
          ptP.Rotate( ptAx, vtAx, dCosStepRot, dSinStepRot) ;
@@ -2705,7 +2725,7 @@ SurfTriMesh::CreateByScrewing( const PolyLine& PL, const Point3d& ptAx, const Ve
       bool bOnAx = bOnlyRev && DistPointLine( ptP, ptAx, vtAx, 1, false).IsSmall() ;
       // ciclo sugli step
       for ( int i = 1 ; i <= nStep ; ++i) {
-        // se non è fisso sull'asse, inserisco le copie ruotate
+        // se non � fisso sull'asse, inserisco le copie ruotate
          if ( ! bOnAx) { 
             ptP.Rotate( ptAx, vtAx, dCosStepRot, dSinStepRot) ;
             if ( ! bOnlyRev)
@@ -2714,8 +2734,8 @@ SurfTriMesh::CreateByScrewing( const PolyLine& PL, const Point3d& ptAx, const Ve
                return false ;
             ++ nV ;
          }
-        // per i controlli già fatti non è possibile avere contemp. prec e corr su asse
-        // se il precedente è sull'asse, aggiungo un solo triangolo
+        // per i controlli gi� fatti non � possibile avere contemp. prec e corr su asse
+        // se il precedente � sull'asse, aggiungo un solo triangolo
          if ( bPrevOnAx) {
             nIdV[0] = nVPrevOnAx ;
             nIdV[1] = nV ;
@@ -2723,7 +2743,7 @@ SurfTriMesh::CreateByScrewing( const PolyLine& PL, const Point3d& ptAx, const Ve
             if ( AddTriangle( nIdV) == SVT_NULL)
                return false ;
          }
-        // se il corrente è sull'asse, aggiungo un solo triangolo
+        // se il corrente � sull'asse, aggiungo un solo triangolo
          else if ( bOnAx) {
             nIdV[0] = nV - ( nStep + 2) + i ;
             nIdV[1] = nIdV[0] + 1 ;
@@ -2743,8 +2763,8 @@ SurfTriMesh::CreateByScrewing( const PolyLine& PL, const Point3d& ptAx, const Ve
       }
      // se rivoluzione completa, aggiungo i due triangoli di chiusura
       if ( bFullRev) {
-        // per i controlli già fatti non è possibile avere contemp. prec e corr su asse
-        // se il precedente è sull'asse, aggiungo un solo triangolo
+        // per i controlli gi� fatti non � possibile avere contemp. prec e corr su asse
+        // se il precedente � sull'asse, aggiungo un solo triangolo
          if ( bPrevOnAx) {
             nIdV[0] = nVPrevOnAx ;
             nIdV[1] = nV - nStep ;
@@ -2752,7 +2772,7 @@ SurfTriMesh::CreateByScrewing( const PolyLine& PL, const Point3d& ptAx, const Ve
             if ( AddTriangle( nIdV) == SVT_NULL)
                return false ;
          }
-        // se il corrente è sull'asse, aggiungo un solo triangolo
+        // se il corrente � sull'asse, aggiungo un solo triangolo
          else if ( bOnAx) {
             nIdV[0] = nV - 1 ;
             nIdV[1] = nV - ( nStep + 1) ;
@@ -2775,7 +2795,7 @@ SurfTriMesh::CreateByScrewing( const PolyLine& PL, const Point3d& ptAx, const Ve
   // altrimenti ultimo punto di polilinea chiusa
    if ( bClosed) {
       for ( int i = 1 ; i <= nStep ; ++i) {
-        // non devo aggiungere i vertici, perchè coincidono con quelli iniziali
+        // non devo aggiungere i vertici, perch� coincidono con quelli iniziali
         // aggiungo triangolo in basso a sinistra
          nIdV[0] = nV - nStep + i - 1 ; nIdV[1] = nV - nStep + i ; nIdV[2] = i ;
          if ( AddTriangle( nIdV) == SVT_NULL)
@@ -2811,7 +2831,7 @@ SurfTriMesh::AddBiTriangle( const int nIdVert[4])
   //    |    |
   //    1 -> 2
   int nIdV[3] ;
-  // se la diagonale 0->2 è uguale o più corta della 1->3
+  // se la diagonale 0->2 � uguale o pi� corta della 1->3
    if ( SqDist( m_vVert[nIdVert[0]].ptP, m_vVert[nIdVert[2]].ptP) <=
         SqDist( m_vVert[nIdVert[1]].ptP, m_vVert[nIdVert[3]].ptP) + EPS_SMALL) {
      // triangolo 0->1->2
@@ -2854,6 +2874,7 @@ SurfTriMesh::DoCompacting( double dTol)
    m_nStatus = ERR ;
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // definisco un Grid per i vertici della superficie
    PointGrid3d VertGrid ;
@@ -2872,14 +2893,14 @@ SurfTriMesh::DoCompacting( double dTol)
       }
      // recupero la posizione geometrica del vertice
       Point3d ptP = m_vVert[nId].ptP ;
-     // se non c'è già un vertice con la stessa posizione lo inserisco nel grid
+     // se non c'� gi� un vertice con la stessa posizione lo inserisco nel grid
       int nAliasId ;
       if ( ! VertGrid.Find( ptP, dTol, nAliasId)) {
          VertGrid.InsertPoint( ptP, nId) ;
         // salvo l'Id nel vettore di reindirizzo
          vVId.push_back( nId) ;
       }
-     // c'è un vertice coincidente
+     // c'� un vertice coincidente
       else {
         // salvo l'Id alias nel vettore di reindirizzo
          vVId.push_back( nAliasId) ;
@@ -2899,7 +2920,7 @@ SurfTriMesh::DoCompacting( double dTol)
      // salto i triangoli cancellati
       if ( vOId[0] == SVT_DEL)
          continue ;
-     // verifico la validità degli indici
+     // verifico la validit� degli indici
       if ( vOId[0] < 0  ||  vOId[0] >= nVIdSize  ||
            vOId[1] < 0  ||  vOId[1] >= nVIdSize  ||
            vOId[2] < 0  ||  vOId[2] >= nVIdSize)
@@ -2941,6 +2962,7 @@ SurfTriMesh::DoSewing( const ISurfTriMesh& stmOther, const Frame3d& frOther, dou
    m_nStatus = ERR ;
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // definisco un Grid per i vertici delle due superfici
    PointGrid3d VertGrid ;
@@ -2970,7 +2992,7 @@ SurfTriMesh::DoSewing( const ISurfTriMesh& stmOther, const Frame3d& frOther, dou
       Point3d ptOP = pOther->m_vVert[nOId].ptP ;
      // la porto nel riferimento della prima superficie
       ptOP.ToGlob( frOther) ;
-     // se non c'è già un vertice con la stessa posizione lo inserisco
+     // se non c'� gi� un vertice con la stessa posizione lo inserisco
       int nNewId ;
       if ( ! VertGrid.Find( ptOP, dTol, nNewId)) {
          if ( ( nNewId = AddVertex( ptOP)) == SVT_NULL)
@@ -2992,7 +3014,7 @@ SurfTriMesh::DoSewing( const ISurfTriMesh& stmOther, const Frame3d& frOther, dou
      // salto i triangoli cancellati
       if ( vOId[0] == SVT_DEL)
          continue ;
-     // verifico la validità degli indici
+     // verifico la validit� degli indici
       if ( vOId[0] < 0  ||  vOId[0] >= nVIdSize  ||
            vOId[1] < 0  ||  vOId[1] >= nVIdSize  ||
            vOId[2] < 0  ||  vOId[2] >= nVIdSize)
@@ -3044,7 +3066,7 @@ SurfTriMesh::GetBBox( const Frame3d& frRef, BBox3d& b3Ref, int nFlag) const
   // verifico lo stato
    if ( m_nStatus != OK)
       return false ;
-  // verifico validità del frame
+  // verifico validit� del frame
    if ( frRef.GetType() == Frame3d::ERR)
       return false ;
   // assegno il box nel riferimento
@@ -3068,9 +3090,10 @@ SurfTriMesh::Translate( const Vector3d& vtMove)
    if ( m_nStatus != OK)
       return false ;
 
-  // imposto ricalcolo della grafica e di hashgrids3d
+  // imposto ricalcolo della grafica e di hashgrids3d e pointgrid3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // traslo i vertici
    for ( int i = 0 ; i < GetVertexSize() ; ++ i) {
@@ -3089,13 +3112,14 @@ SurfTriMesh::Rotate( const Point3d& ptAx, const Vector3d& vtAx, double dCosAng,
    if ( m_nStatus != OK)
       return false ;
 
-  // verifico validità dell'asse di rotazione
+  // verifico validit� dell'asse di rotazione
    if ( vtAx.IsSmall())
       return false ;
 
-  // imposto ricalcolo della grafica e di hashgrids3d
+  // imposto ricalcolo della grafica e di hashgrids3d e pointgrid3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // ruoto i vertici
    for ( int i = 0 ; i < GetVertexSize() ; ++ i) {
@@ -3155,7 +3179,7 @@ SurfTriMesh::Scale( const Frame3d& frRef, double dCoeffX, double dCoeffY, double
    bool bRecalc = ( abs( dCoeffX) < EPS_ZERO  ||  abs( dCoeffY) < EPS_ZERO  ||  abs( dCoeffZ) < EPS_ZERO) ;
    for ( int i = 0 ; i < GetTriangleSize() ; ++ i) {
       if ( m_vTria[i].nIdVert[0] != SVT_DEL) {
-        // se c'è mirror, devo invertire la faccia
+        // se c'� mirror, devo invertire la faccia
          if ( bMirror)
             InvertTriangle( i) ;
         // aggiorno la normale
@@ -3177,7 +3201,7 @@ SurfTriMesh::Scale( const Frame3d& frRef, double dCoeffX, double dCoeffY, double
 bool 
 SurfTriMesh::InvertTriangle( int nT)
 {
-  // controllo validità triangolo
+  // controllo validit� triangolo
    if ( m_vTria[nT].nIdVert[0] == SVT_DEL)
       return true ;
   // scambio di due vertici
@@ -3193,10 +3217,10 @@ SurfTriMesh::InvertTriangle( int nT)
 bool 
 SurfTriMesh::CalcTriangleNormal( int nT)
 {
-  // controllo validità triangolo
+  // controllo validit� triangolo
    if ( m_vTria[nT].nIdVert[0] == SVT_DEL)
       return true ;
-  // controllo validità vertici riferiti dal triangolo
+  // controllo validit� vertici riferiti dal triangolo
    if ( m_vTria[nT].nIdVert[0] < 0  ||
         m_vTria[nT].nIdVert[0] >= GetVertexSize()  ||
         m_vVert[m_vTria[nT].nIdVert[0]].nIdTria == SVT_DEL  ||
@@ -3226,13 +3250,14 @@ SurfTriMesh::Mirror( const Point3d& ptOn, const Vector3d& vtNorm)
    if ( m_nStatus != OK)
       return false ;
 
-  // verifico validità del piano di specchiatura
+  // verifico validit� del piano di specchiatura
    if ( vtNorm.IsSmall())
       return false ;
 
-  // imposto ricalcolo della grafica e di hashgrids3d
+  // imposto ricalcolo della grafica e di hashgrids3d e pointgrid3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // specchio i vertici
    for ( int i = 0 ; i < GetVertexSize() ; ++ i)
@@ -3259,13 +3284,14 @@ SurfTriMesh::Shear( const Point3d& ptOn, const Vector3d& vtNorm, const Vector3d&
    if ( m_nStatus != OK)
       return false ;
 
-  // verifico validità dei parametri
+  // verifico validit� dei parametri
    if ( vtNorm.IsSmall()  ||  vtDir.IsSmall())
       return false ;
 
-  // imposto ricalcolo della grafica e di hashgrids3d
+  // imposto ricalcolo della grafica e di hashgrids3d e pointgrid3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // eseguo scorrimento dei vertici
    for ( int i = 0 ; i < GetVertexSize() ; ++ i) {
@@ -3291,13 +3317,14 @@ SurfTriMesh::ToGlob( const Frame3d& frRef)
    if ( m_nStatus != OK)
       return false ;
 
-  // verifico validità del frame
+  // verifico validit� del frame
    if ( frRef.GetType() == Frame3d::ERR)
       return false ;
 
-  // imposto ricalcolo della grafica e di hashgrids3d
+  // imposto ricalcolo della grafica e di hashgrids3d e pointgrid3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // trasformo i vertici
    for ( int i = 0 ; i < GetVertexSize() ; ++ i) {
@@ -3322,13 +3349,14 @@ SurfTriMesh::ToLoc( const Frame3d& frRef)
    if ( m_nStatus != OK)
       return false ;
 
-  // verifico validità del frame
+  // verifico validit� del frame
    if ( frRef.GetType() == Frame3d::ERR)
       return false ;
 
-  // imposto ricalcolo della grafica e di hashgrids3d
+  // imposto ricalcolo della grafica e di hashgrids3d e pointgrid3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // trasformo i vertici
    for ( int i = 0 ; i < GetVertexSize() ; ++ i) {
@@ -3353,17 +3381,18 @@ SurfTriMesh::LocToLoc( const Frame3d& frOri, const Frame3d& frDest)
    if ( m_nStatus != OK)
       return false ;
 
-  // verifico validità dei frame
+  // verifico validit� dei frame
    if ( frOri.GetType() == Frame3d::ERR  ||  frDest.GetType() == Frame3d::ERR)
       return false ;
 
-  // se i due riferimenti coincidono, non devo fare alcunché
+  // se i due riferimenti coincidono, non devo fare alcunch�
    if ( AreSameFrame( frOri, frDest))
       return true ;
 
-  // imposto ricalcolo della grafica e di hashgrids3d
+  // imposto ricalcolo della grafica e di hashgrids3d e pointgrid3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // trasformo i vertici
    for ( int i = 0 ; i < GetVertexSize() ; ++ i) {
@@ -3388,9 +3417,10 @@ SurfTriMesh::Invert( void)
    if ( m_nStatus != OK)
       return false ;
 
-   // imposto ricalcolo della grafica e di hashgrids3d
+   // imposto ricalcolo della grafica e di hashgrids3d e pointgrid3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
   // inverto i triangoli
    for ( int i = 0 ; i < GetTriangleSize() ; ++ i)
@@ -3403,7 +3433,7 @@ SurfTriMesh::Invert( void)
 bool
 SurfTriMesh::Cut( const Plane3d& plPlane, bool bSaveOnEq)
 {
-#define UseTria 1
+#define UseTria 0
 #if UseTria
    // la superficie deve essere validata
    if ( m_nStatus != OK)
@@ -3619,10 +3649,10 @@ SurfTriMesh::Cut( const Plane3d& plPlane, bool bSaveOnEq)
    return true ;
 
 #else
- 
+   
    bool bModif = false ;
 
-  // Setto i triangoli come né fuori né dentro.
+  // Setto i triangoli come n� fuori n� dentro.
    int nTriaNum = GetTriangleSize() ;
    for ( int nT = 0 ; nT < nTriaNum ; ++ nT) {
       m_vTria[nT].nTempPart = 0 ;
@@ -3641,8 +3671,11 @@ SurfTriMesh::Cut( const Plane3d& plPlane, bool bSaveOnEq)
       PtrOwner<SurfFlatRegion> pReg( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLineVector( vLoopVec))) ;
       if ( IsNull( pReg))
          return false ;
+      double dArea;
+      Plane3d plFacetPlane;
+      vLoopVec[0].IsClosedAndFlat( plFacetPlane, dArea) ;
       LineFacetClassVector IntersLinePart ;
-      int nIntType = IntersFacetPlane( *pReg, vLoopVec[0], plPlane, IntersLinePart) ;
+      int nIntType = IntersFacetPlane( vLoopVec, plPlane, IntersLinePart) ; /*IntersFacetPlane(*pReg, vLoopVec[0], plPlane, IntersLinePart);*/
       if ( nIntType == FacetPlaneIntersType::FPI_CUT) {
          for ( int nPart = 0 ; nPart < int( IntersLinePart.size()) ; ++ nPart) {
             // Salvo intersezione per la faccia.
@@ -3653,13 +3686,13 @@ SurfTriMesh::Cut( const Plane3d& plPlane, bool bSaveOnEq)
             else {
                IntersLineMap.emplace( nF, IntersInnChain( 1, IntersInnSeg( IntersLinePart[nPart].ptSt, IntersLinePart[nPart].ptEn))) ;
             }
-         }
+         } 
       }
-      else if ( nIntType == FacetPlaneIntersType::FPI_ON) { 
+      else if ( nIntType == FacetPlaneIntersType::FPI_ON) {
          INTVECTOR vT ;
          GetAllTriaInFacet( nF, vT) ;
          for ( auto& nT : vT)
-            m_vTria[nT].nTempPart = pReg->GetNormVersor() * plPlane.GetVersN() > 0. ? 2 : - 2 ;
+            m_vTria[nT].nTempPart = plFacetPlane.GetVersN() * plPlane.GetVersN() > 0. ? 2 : - 2 ; /*pReg->GetNormVersor() * plPlane.GetVersN() > 0. ? 2 : - 2 ;*/
       }
       else if ( nIntType == FacetPlaneIntersType::FPI_INN) {
          INTVECTOR vT ;
@@ -3677,18 +3710,24 @@ SurfTriMesh::Cut( const Plane3d& plPlane, bool bSaveOnEq)
          ;
    }
 
-  // Divido le facce.
+  // Divido e ritriangolo le facce.
    PieceMap NewFacet ;
    SplitFacet( IntersLineMap, NewFacet) ;
    INTERSEDGEMAP EdgeInterLineMap, EdgeEdgeLineMap ;
    RetriangulateFacetPieces( NewFacet, EdgeInterLineMap, EdgeEdgeLineMap) ;
-
+   if ( ! AdjustVertices()  ||  ! DoCompacting())
+      return false ;
+   // Identifico le parti
+   IdentifyParts() ;
+   // Elimino i triangoli fuori dal semi-spazio d'interesse.
    int nNumTria = GetTriangleSize() ;
    for ( int nT = 0 ; nT < nNumTria ; ++ nT)
       if ( m_vTria[nT].nTempPart == - 1  ||  m_vTria[nT].nTempPart == - 2  ||  ( ! bSaveOnEq  &&  m_vTria[nT].nTempPart == 2))
          RemoveTriangle( nT) ;
 
-   return AdjustVertices()  &&  DoCompacting() ;
+   if ( AdjustVertices()  &&  DoCompacting() )
+      RemoveTJunctions();
+   return AdjustVertices()  &&  DoCompacting() ; 
 #endif
 }
 
@@ -3707,7 +3746,7 @@ SurfTriMesh::ResetHashGrids3d( void) const
 bool
 SurfTriMesh::VerifyHashGrids3d( void) const
 {
-  // se già calcolato, non devo fare altro
+  // se gi� calcolato, non devo fare altro
    if ( m_pHGrd3d != nullptr)
       return true ;
   // alloco
@@ -3736,6 +3775,45 @@ SurfTriMesh::VerifyHashGrids3d( void) const
    return true ;
 }
 
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+//----------------------------------------------------------------------------
+void
+SurfTriMesh::ResetPointGrid3d( void) const
+{
+   if ( m_pPGrd3d != nullptr) {
+      delete m_pPGrd3d ; 
+      m_pPGrd3d = nullptr ;
+   } 
+}
+
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::VerifyPointGrid3d( void) const
+{
+   // Se gi� calcolato, non devo fare altro
+   if ( m_pPGrd3d != nullptr)
+      return true ;
+   // Alloco
+   m_pPGrd3d = new PointGrid3d ;
+   if ( m_pPGrd3d == nullptr)
+      return false ;
+   // Riempio
+   m_pPGrd3d->Init( 100) ;
+   Point3d ptPV ;
+   int nV = GetFirstVertex( ptPV) ;
+   while ( nV != SVT_NULL)
+   {
+      if ( ! m_pPGrd3d->InsertPoint( ptPV, nV)) {
+         ResetPointGrid3d() ;
+         return false ;
+      }
+      nV = GetNextVertex( nV, ptPV) ;   
+   }
+
+   return true ; 
+}
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
 //----------------------------------------------------------------------------
 bool
 SurfTriMesh::GetAllTriaOverlapBox( const BBox3d& b3Box, INTVECTOR& vT) const
@@ -3778,7 +3856,7 @@ SurfTriMesh::VerifyConnection( void) const
   // ciclo sui triangoli
    m_nParts = 0 ;
    for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i) {
-     // salto triangoli cancellati o già assegnati
+     // salto triangoli cancellati o gi� assegnati
       if ( m_vTria[i].nIdVert[0] == SVT_DEL  ||
            m_vTria[i].nPart != SVT_NULL)
          continue ;
@@ -3831,9 +3909,10 @@ SurfTriMesh::RemovePart( int nPart)
   // Aggiorno il numero di componenti
    m_nParts = nPartsOld - 1 ;
 
-  // imposto ricalcolo della grafica e di hashgrids3d
+  // imposto ricalcolo della grafica e di hashgrids3d e poingrid3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
    return true ;
 }

SurfTriMesh.h

@@ -18,6 +18,7 @@
 #include "GeoObjRW.h"
 #include "/EgtDev/Include/EGkSurfTriMesh.h"
 #include "/EgtDev/Include/EGkHashGrids3d.h"
+#include "/EgtDev/Include/EGkPointGrid3d.h"
 #include <deque>
 #include <set>
 
@@ -133,6 +134,10 @@ typedef std::vector<LineFacetClass> LineFacetClassVector ;
 //----------------------------------------------------------------------------
 enum FacetPlaneIntersType { FPI_ERROR = 0, FPI_CUT = 1, FPI_INN = 2, FPI_OUT = 3, FPI_ON = 4 } ;
 
+
+//----------------------------------------------------------------------------
+enum LineClassificationWithExtPolygon { LCP_ERROR = 0, LCP_INN = 1, LCP_OUT = 2, LCP_OVERLAP = 3 } ;
+
 //----------------------------------------------------------------------------
 // Definizione strutture e contenitori
 // Contatti interno-interno
@@ -226,7 +231,7 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       bool CopyFrom( const IGeoObj* pGObjSrc) override ;
       bool Init( int nNumVert, int nNumTria, int nNumFacet = 0) override ;
       void SetLinearTolerance( double dLinTol) override
-              { m_dLinTol = std::max( dLinTol, EPS_SMALL) ; }
+              { m_dLinTol = EPS_SMALL/*std::max( dLinTol, EPS_SMALL)*/ ; } ///////////////////////////////////////////////////////////////
       void SetBoundaryAngle( double dBoundaryAngDeg) override
               { m_dBoundaryAng = std::max( dBoundaryAngDeg, EPS_ANG_SMALL) ;
                 m_dCosBndAng = cos( m_dBoundaryAng * DEGTORAD) ;
@@ -380,6 +385,8 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       bool MarchOneFacetTria( int nF, int& nT, int& nV, int nTimeStamp, PolyLine& PL, bool& bEnd) const ;
       void ResetHashGrids3d( void) const ;
       bool VerifyHashGrids3d( void) const ;
+      void ResetPointGrid3d( void) const ;
+      bool VerifyPointGrid3d( void) const ;
       bool VerifyConnection( void) const ;
       bool DecomposeLoop( CHAINVECTOR& cvOpenChain, INTVECTOR& vnDegVec, PNTMATRIX& cvBoundClosedLoopVec, BOOLVECTOR& vbInOut) ;
       bool RetriangulationForBooleanOperation( CHAINMAP& LoopLines, TRIA3DVECTORMAP& Ambiguos, SurfTriMesh& Surf, bool& bModif) ;
@@ -387,10 +394,12 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       bool IntersectTriMeshTriangle( SurfTriMesh& Other) ;
       int IntersFacetPlane( const SurfFlatRegion& Region, const PolyLine& ExtLoop, const Plane3d& plCutPlane,
                             LineFacetClassVector& IntersLinePart) ;
+      int IntersFacetPlane( const POLYLINEVECTOR& vFacetLoopsVec, const Plane3d& plCutPlane,
+                            LineFacetClassVector& IntersLinePart) ;
       bool IntersFacetFacet( const SurfFlatRegion& RegionA, const PolyLine& ExtLoopA,
                              const SurfFlatRegion& RegionB, const PolyLine& ExtLoopB,
                              LineFacetClassVector& IntersLinePart) ;
-      bool ItersectTriMeshFacets( SurfTriMesh& Other) ;
+      bool IntersectTriMeshFacets( SurfTriMesh& Other) ;
       bool RetriangulateFacetPieces( const PieceMap& NewFacet,
                                      const INTERSEDGEMAP& EdgeInterLineMap,
                                      const INTERSEDGEMAP& EdgeEdgeLineMap) ;
@@ -416,20 +425,29 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       double         m_dCosBndAng ;    // coseno dell'angolo limite per considerare un lato un contorno  
       double         m_dSmoothAng ;    // angolo limite per mediare le normali (in gradi)
       double         m_dCosSmAng ;     // coseno dell'angolo limite per mediare le normali
-      bool           m_bOriented ;     // la superficie è orientata consistentemente in tutte le sue parti
+      bool           m_bOriented ;     // la superficie � orientata consistentemente in tutte le sue parti
       bool           m_bClosed ;       // la superficie racchiude un volume
-      bool           m_bFaceted ;      // flag di validità della sfaccettatura
+      bool           m_bFaceted ;      // flag di validit� della sfaccettatura
       VERTVECTOR     m_vVert ;         // vettore dei vertici
       TRIAVECTOR     m_vTria ;         // vettore dei triangoli
       INTVECTOR      m_vFacet ;        // vettore delle sfaccettature
       mutable int    m_nTimeStamp ;    // orologio locale 
-      int            m_nTempProp[2] ;  // vettore proprietà temporanee
+      int            m_nTempProp[2] ;  // vettore propriet� temporanee
       int            m_nMaxTFlag ;     // massimo valore dei TFlag dei triangoli
       mutable int    m_nParts ;        // numero di parti connesse (-1 se da calcolare)
       mutable HashGrids3d* m_pHGrd3d ; // Hash Grid 3d nel suo riferimento
-      mutable BBox3d m_b3HGrd3d ;      // Box3d collegato a Hash Grid 3d 
+      mutable BBox3d m_b3HGrd3d ;      // Box3d collegato a Hash Grid 3d
+      mutable PointGrid3d* m_pPGrd3d ; // Point Grid 3d nel suo riferimento
 } ;
 
+//----------------------------------------------------------------------------
+//static bool ChangePolyLineStart( const Point3d& ptNewStart, PolyLine& Loop) ;
+//// nSegNum 0-based
+//static bool PointPositionOnPolyLine( const Point3d& ptPoint, /*const*/ PolyLine& Loop, int& nSegNum, double& dParOnSeg) ;
+//static bool IsPointInsidePolyLine( const Point3d& ptP, /*const*/ PolyLine& plPoly) ;
+//static bool SplitPolyLineAtPoint( const Point3d& ptPoint, /*const*/ PolyLine& Loop, PolyLine& Loop1, PolyLine& Loop2) ;
+//static bool AddPolyLineToPolyLine(PolyLine& Poly, PolyLine& PolyToAdd) ;
+
 //-----------------------------------------------------------------------------
 inline SurfTriMesh* CreateBasicSurfTriMesh( void)
                { return ( static_cast<SurfTriMesh*>( CreateGeoObj( SRF_TRIMESH))) ; }

SurfTriMeshFaceting.cpp

@@ -114,7 +114,7 @@ SurfTriMesh::UpdateTriaFaceting( int nRefT, int nFacet, const Plane3d& plPlane,
    }
    if ( nV == SVT_NULL)
       return false ;
-   double dTol = max( min( m_dLinTol, 20 * EPS_SMALL), 2 * EPS_SMALL) ;
+   double dTol = max( min( m_dLinTol, 20 * EPS_SMALL), /*2 **/ EPS_SMALL) ;     ////////////////////////////////////////////////////
    if ( ! PointInPlaneEpsilon( m_vVert[m_vTria[nT].nIdVert[nV]].ptP, plPlane, dTol))
       return true ;
   // il triangolo fa parte della faccia
@@ -413,7 +413,7 @@ SurfTriMesh::GetFacetLoops( int nF, POLYLINEVECTOR& vPL) const
    for ( int i = 0 ; i < int( vPL.size()) ; ++ i) {
       Plane3d plPlane ;
       double dArea ;
-      if ( ! vPL[i].IsClosedAndFlat( plPlane, dArea, 100 * EPS_SMALL))
+      if ( ! vPL[i].IsClosedAndFlat( plPlane, dArea, /*100 **/ EPS_SMALL))
          return false ;
      // se loop esterno
       if ( vtN * plPlane.GetVersN() > 0) {
@@ -712,9 +712,10 @@ SurfTriMesh::RemoveFacet( int nF)
    if ( ! DoCompacting())
       return false ;
 
-  // dichiaro necessità ricalcolo della grafica e di hashgrids3d
+  // dichiaro necessità ricalcolo della grafica e di hashgrids3d e pointgrid3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
+   ResetPointGrid3d() ;
 
    return true ;
 }

Tool.cpp

@@ -29,7 +29,8 @@ using namespace std ;
 Tool::Tool( bool bApproxWithLines) 
          : m_bApproxWithLines( bApproxWithLines), m_dLinTol( LIN_TOL_STD), m_dAngTolDeg( ANG_TOL_APPROX_DEG),
            m_nType( UNDEF), m_nCurrentNum( 0), m_dHeight( 0), m_dTipHeight( 0), m_dRadius( 0), m_dRCorner( 0),
-           m_dTipRadius( 0), m_dRefRadius( 0), m_dCutterHeight( 0), m_dMrtChsWidth( 0), m_dMrtChsThickness( 0)
+           m_dTipRadius( 0), m_dRefRadius( 0), m_dCutterHeight( 0), m_dMrtChsWidth( 0), m_dMrtChsThickness( 0)/*,
+           m_bAllPartCut( true)*/
 {
 }
 
@@ -73,6 +74,10 @@ Tool::SetTolerances( double dLinTol, double dAngTolDeg)
 bool
 Tool::SetStdTool( const string& sToolName, double dH, double dR, double dCornR, double dCutterH, int nToolNum)
 {
+   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+   //////////////////////////////////////// Per test additivi //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+   //return SetAdditiveTool( sToolName, dH, dR, dCornR, nToolNum) ;
+   /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
   // Impostazioni generali
    m_sName = sToolName ;
    m_nCurrentNum = nToolNum ;
@@ -100,9 +105,13 @@ Tool::SetStdTool( const string& sToolName, double dH, double dR, double dCornR,
       Point3d pt3( m_dRadius, - m_dHeight, 0) ;
       Point3d pt4( 0, - m_dHeight, 0) ;
       m_Outline.AddPoint( pt0) ;
-      m_Outline.AddLine( pt1);
+      m_Outline.AddLine( pt1) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
       m_Outline.AddLine( pt3) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       m_Outline.AddLine( pt4) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
    }
   // utensile naso di toro
    else if ( dCornR < dR - EPS_SMALL) {
@@ -120,10 +129,15 @@ Tool::SetStdTool( const string& sToolName, double dH, double dR, double dCornR,
       Point3d pt3( m_dTipRadius, - m_dHeight, 0) ;
       Point3d pt4( 0, - m_dHeight, 0) ;
       m_Outline.AddPoint( pt0) ;
-      m_Outline.AddLine( pt1);
-      m_Outline.AddLine( pt2);
+      m_Outline.AddLine( pt1) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
+      m_Outline.AddLine( pt2) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       m_Outline.AddArcTg( pt3) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
       m_Outline.AddLine( pt4) ;
+      m_Outline.SetCurveTempProp( 3, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
      // se da approosimare
       if ( m_bApproxWithLines)
         return SetGenTool( sToolName, &m_Outline, nToolNum) ;
@@ -143,9 +157,13 @@ Tool::SetStdTool( const string& sToolName, double dH, double dR, double dCornR,
       Point3d pt2( m_dRadius, - m_dHeight + m_dTipHeight, 0) ;
       Point3d pt4( 0, - m_dHeight, 0) ;
       m_Outline.AddPoint( pt0) ;
-      m_Outline.AddLine( pt1);
-      m_Outline.AddLine( pt2);
+      m_Outline.AddLine( pt1) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
+      m_Outline.AddLine( pt2) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       m_Outline.AddArcTg( pt4) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
    }
   // impossibile
    else
@@ -201,9 +219,14 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
      // profilo
       m_Outline.AddPoint( pt0) ;
       m_Outline.AddLine( pt1) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
       m_Outline.AddLine( pt2) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       m_Outline.AddLine( Point3d( m_dTipRadius, - m_dHeight, 0)) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
       m_Outline.AddLine( pt5) ;
+      m_Outline.SetCurveTempProp( 3, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
      // eventuali sistemazioni per altezza tagliente
       if ( ModifyForCutterHeight())
          return SetGenTool( sToolName, &m_Outline, nToolNum) ;
@@ -213,7 +236,7 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
 
   // Altrimenti utensile generico
 
-  // se Tip a punta ( TipRadius è raggio teorico della parte finale dell'utensile senza raccordo)
+  // se Tip a punta ( TipRadius � raggio teorico della parte finale dell'utensile senza raccordo)
    if ( m_dTipRadius < m_dRadius) {
      // se punta a sfera
       if ( m_dTipRadius < EPS_SMALL) {
@@ -242,10 +265,16 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
         // creazione curva composita
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt1) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt2) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
          m_Outline.AddLine( pt3) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
          m_Outline.AddCurve( cvArc) ;
+         m_Outline.SetCurveTempProp( 3, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 4, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
      // altrimenti punta a naso di toro
       else {
@@ -271,16 +300,22 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
         // creazione curva composita
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt1) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt2) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
          m_Outline.AddLine( pt3) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
          m_Outline.AddCurve( Release( pArc)) ;
+         m_Outline.SetCurveTempProp( 3, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 4, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
    }
 
-   // altrimenti Tip a coda di rondine ( TipRadius è raggio misurabile della parte finale dell'utensile)
+   // altrimenti Tip a coda di rondine ( TipRadius � raggio misurabile della parte finale dell'utensile)
    else {
-     // il raggio della punta non può essere inferiore al raggio corner
+     // il raggio della punta non pu� essere inferiore al raggio corner
       if ( m_dTipRadius < m_dRCorner)
          return false ;
      // Assegno il raggio di riferimento
@@ -305,10 +340,16 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
      // creazione curva composita
       m_Outline.AddPoint( pt0) ;
       m_Outline.AddLine( pt1) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
       m_Outline.AddLine( pt2) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       m_Outline.AddLine( pt3) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
       m_Outline.AddCurve( cvArc) ;
+      m_Outline.SetCurveTempProp( 3, 1, 1) ;
       m_Outline.AddLine( pt5) ;
+      m_Outline.SetCurveTempProp( 4, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
    }
 
   // eventuali sistemazioni per altezza tagliente
@@ -321,7 +362,7 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
 bool
 Tool::ModifyForCutterHeight( void)
 {
-  // Se altezza tagliente non definita o superiore alla altezza utensile non devo fare alcunché
+  // Se altezza tagliente non definita o superiore alla altezza utensile non devo fare alcunch�
    if ( m_dCutterHeight < EPS_SMALL  ||  m_dCutterHeight > m_dHeight - EPS_SMALL)
       return false ;
   // quota di taglio
@@ -340,6 +381,8 @@ Tool::ModifyForCutterHeight( void)
       m_Outline.TrimStartAtParam( dU) ;
       m_Outline.AddLine( Point3d( 0, dYtrim, 0), false) ;
       m_Outline.AddLine( ORIG, false) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       return true ;
    } 
    return false ;
@@ -375,22 +418,36 @@ Tool::SetSawTool( const string& sToolName, double dH, double dR,
         // creazione profilo
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt1) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt2) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
          m_Outline.AddLine( pt3) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
          m_Outline.AddLine( pt4) ;
+         m_Outline.SetCurveTempProp( 3, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 4, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
      // altrimenti con raggio corner
       else {
         // creazione profilo
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt1) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt2) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
          m_Outline.AddLine( pt3 - X_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddArcTg( pt3 - Y_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
          m_Outline.AddLine( pt4 + Y_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 3, 1, 1) ;
          m_Outline.AddArcTg( pt4 - X_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 4, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 5, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
    }
   // altrimenti senza gambo
@@ -405,18 +462,28 @@ Tool::SetSawTool( const string& sToolName, double dH, double dR,
         // creazione profilo
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt3) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt4) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
      // altrimenti con raggio corner
       else {
         // creazione profilo
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt3 - X_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddArcTg( pt3 - Y_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt4 + Y_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddArcTg( pt4 - X_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
    }
   // Assegno il raggio di riferimento
@@ -436,17 +503,30 @@ Tool::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline,
    if ( pToolOutline != &m_Outline)
       m_Outline.Clear() ;
    m_ArcLineApprox.Clear() ;
-
+                                                           
   // Copio il profilo e garantisco sia di soli archi e rette (converto eventuali curve di Bezier)
    if ( ! m_Outline.CopyFrom( pToolOutline)  ||
         ! m_Outline.ArcsBezierCurvesToArcsPerpExtr( m_dLinTol, m_dAngTolDeg))
       return false ;
 
-  // Ciclo sulle curve componenti 
+  // Valuto se tutte le curve tagliano
+   m_ArcLineApprox.SetTempProp( 1, 1) ;
    const ICurve* pCurve = m_Outline.GetFirstCurve() ;
+   while ( pCurve != nullptr  &&  m_ArcLineApprox.GetTempProp( 1) == 1) {
+      int nCutTempProp = pCurve->GetTempProp( 1) ;
+      if ( nCutTempProp == 0) {
+         m_ArcLineApprox.SetTempProp( 0, 1) ;
+      }
+      pCurve = m_Outline.GetNextCurve() ;
+   }
+
+  // Ciclo sulle curve componenti 
+   pCurve = m_Outline.GetFirstCurve() ;
    while ( pCurve != nullptr) {
 
-     // Se la curva è un arco ed è richiesto la verifica per l'approssimazione,
+      int nCutTempProp = pCurve->GetTempProp( 1) ;
+
+     // Se la curva � un arco ed � richiesto la verifica per l'approssimazione,
      // verifico se approssimarlo
       if ( m_bApproxWithLines  &&  pCurve->GetType() == CRV_ARC) {
 
@@ -454,11 +534,11 @@ Tool::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline,
          Point3d ptO = GetBasicCurveArc( pCurve)->GetCenter() ;
          double dRadius = GetBasicCurveArc( pCurve)->GetRadius() ;
 
-        // Se il centro è fuori dall'asse devo approssimare
+        // Se il centro � fuori dall'asse devo approssimare
          bool bCurrApprox = ( abs( ptO.x) > EPS_SMALL) ;
 
         // Se una delle altre curve dista dal centro meno del raggio, devo approssimare
-         for ( int nI = 0 ; ! bCurrApprox  &&  nI < m_Outline.GetCurveCount() ; ++ nI) {
+         for ( int nI = 0 ; ! bCurrApprox  &&  nI < m_Outline.GetCurveCount()  &&  m_ArcLineApprox.GetTempProp( 1) == 1 ; ++ nI) {
             const ICurve* pOtherCrv = m_Outline.GetCurve( nI) ; 
             if ( pOtherCrv != pCurve) {
                DistPointCurve CalcDist( ptO, *pOtherCrv) ;
@@ -493,6 +573,7 @@ Tool::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline,
                m_vArcNormals.emplace_back( vtExtN) ;
                int nCvCount = m_ArcLineApprox.GetCurveCount() ;
                m_ArcLineApprox.SetCurveTempProp( nCvCount - 1, int( m_vArcNormals.size()) - 1) ;
+               m_ArcLineApprox.SetCurveTempProp( nCvCount - 1, nCutTempProp, 1) ;
             }
          }
         // altrimenti lo aggiungo semplicemente 
@@ -500,16 +581,20 @@ Tool::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline,
             m_ArcLineApprox.AddCurve( *pCurve, true) ;
             int nCvCount = m_ArcLineApprox.GetCurveCount() ;
             m_ArcLineApprox.SetCurveTempProp( nCvCount - 1, - 1) ;
+            m_ArcLineApprox.SetCurveTempProp( nCvCount - 1, nCutTempProp, 1) ;
          }
       }
-     // altrimenti è segmento e lo aggiungo semplicemente
-      else
+     // altrimenti � segmento e lo aggiungo semplicemente
+      else {
          m_ArcLineApprox.AddCurve( *pCurve, true) ;
+         int nCvCount = m_ArcLineApprox.GetCurveCount() ;
+         m_ArcLineApprox.SetCurveTempProp( nCvCount - 1, nCutTempProp, 1) ;
+      }
 
       pCurve = m_Outline.GetNextCurve() ;
    }
 
-  // Il profilo dell'utensile deve stare nel 1° e 4° quadrante del piano XY
+  // Il profilo dell'utensile deve stare nel 1� e 4� quadrante del piano XY
    BBox3d Bounding ;
    m_Outline.GetLocalBBox( Bounding) ;
    if ( Bounding.GetMin().x < - 10 * EPS_SMALL)
@@ -521,7 +606,7 @@ Tool::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline,
   // Assegno le dimensioni dell'utensile
    m_dHeight = - Bounding.GetMin().y ;
    m_dRadius = Bounding.GetMax().x ;
-  // Assegno il raggio di riferimento se non già assegnato
+  // Assegno il raggio di riferimento se non gi� assegnato
    if ( m_dRefRadius < EPS_SMALL)
       m_dRefRadius = 0.25 * m_dRadius ;
  
@@ -579,3 +664,85 @@ Tool::SetChiselTool( const string& sToolName, double dH, double dW, double dTh,
 
    return true ;
 }
+
+//----------------------------------------------------------------------------
+bool
+Tool::SetAdditiveTool( const std::string& sToolName, double dH, double dR, double dRC, int nToolNum)
+{
+   // Impostazioni generali
+   m_sName = sToolName ;
+   m_nCurrentNum = nToolNum ;
+   m_nType = UNDEF ;
+   m_Outline.Clear() ;
+   m_ArcLineApprox.Clear() ;
+
+   // verifica sulle minime dimensioni globali
+   if ( dH < EPS_SMALL  ||  dR < EPS_SMALL  ||  dRC < - EPS_SMALL)
+      return false ;
+
+   m_nType = ADDITIVE ;
+   m_dHeight = dH ;
+   m_dRadius = dR ;
+   m_dRCorner = dRC ;
+   m_dTipHeight = 0 ;
+   m_dTipRadius = 0 ;
+   m_dRefRadius = 0 ;
+   m_dCutterHeight = dH ;
+
+   double dSquareCornerRadProj = m_dRCorner * m_dRCorner - 0.25 * m_dHeight * m_dHeight ;
+   // Utensile sfiancato
+   if ( dSquareCornerRadProj > 0) {
+      double dCenX = m_dRadius - m_dRCorner ;
+      double dCylRad = dCenX + sqrt( dSquareCornerRadProj) ;
+      // Utensile mal definito
+      if ( dCylRad < 0)
+         return false ;
+      // Profilo
+      m_Outline.AddPoint( Point3d( 0, 0, 0)) ;
+      m_Outline.AddLine( Point3d( dCylRad, 0, 0)) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
+      CurveArc cvArc ;
+      cvArc.SetC2P( Point3d( dCenX, - 0.5 * m_dHeight, 0), Point3d( dCylRad, 0, 0), Point3d( dCylRad, - m_dHeight, 0)) ;
+      m_Outline.AddCurve( cvArc) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
+      m_Outline.AddLine( Point3d( 0, - m_dHeight, 0)) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
+   }
+   // Utensile cilindrico con eventuale raggio corner 
+   else {
+      // Utensile mal definito
+      if ( m_dRadius - m_dRCorner < EPS_SMALL)
+         return false ;
+      // Raggio corner nullo
+      if ( m_dRadius < EPS_SMALL) {
+         // Profilo
+         m_Outline.AddPoint( Point3d( 0, 0, 0)) ;
+         m_Outline.AddLine( Point3d( m_dRadius, 0, 0)) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
+         m_Outline.AddLine( Point3d( m_dRadius, - m_dHeight, 0)) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
+         m_Outline.AddLine( Point3d( 0, - m_dHeight, 0)) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
+      }
+      else {
+         // Profilo
+         m_Outline.AddPoint( Point3d( 0, 0, 0)) ;
+         m_Outline.AddLine( Point3d( m_dRadius - m_dRCorner, 0, 0)) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
+         CurveArc cvArc ;
+         cvArc.SetC2P( Point3d( m_dRadius - m_dRCorner, - m_dRCorner, 0), Point3d( m_dRadius - m_dRCorner, 0, 0), Point3d( m_dRadius, - m_dRCorner, 0)) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
+         m_Outline.AddLine( Point3d( m_dRadius, - m_dHeight + m_dRCorner, 0)) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
+         cvArc.SetC2P( Point3d( m_dRadius - m_dRCorner, - m_dHeight + m_dRCorner, 0), Point3d( m_dRadius, - m_dHeight + m_dRCorner, 0), Point3d( m_dRadius - m_dRCorner, - m_dHeight, 0)) ;
+         m_Outline.AddCurve( cvArc) ;
+         m_Outline.SetCurveTempProp( 3, 1, 1) ;
+         m_Outline.AddLine( Point3d( 0, - m_dHeight, 0)) ;
+         m_Outline.SetCurveTempProp( 4, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
+      }
+   } 
+
+   return SetGenTool( sToolName, &m_Outline, nToolNum) ;  
+}

Tool.h

@@ -16,7 +16,7 @@
 #include "CurveComposite.h"
 
 //----------------------------------------------------------------------------
-class Tool
+class Tool  
 {
    public :
       Tool( bool bApproxWithLines = false) ;
@@ -33,6 +33,7 @@ class Tool
       bool SetGenTool( const std::string& sToolName, const ICurveComposite* pToolOutline, int nToolNum) ;
       bool SetMortiserTool( const std::string& sToolName, double dH, double dW, double dTh, double dRc, int nToolNum) ;
       bool SetChiselTool( const std::string& sToolName, double dH, double dW, double dTh, int nToolNum) ;
+      bool SetAdditiveTool( const std::string& sToolName, double dH, double dR, double dRC, int nToolNum) ;
       bool SetToolNum( int nToolNum)
                { m_nCurrentNum = nToolNum ; return true ; }
       int GetType() const
@@ -61,6 +62,8 @@ class Tool
                { return ( m_ArcLineApprox.GetCurveCount() == 0 ? m_Outline : m_ArcLineApprox) ; }
       const VCT3DVECTOR& GetArcNormalVec( void) const
                { return m_vArcNormals ; }
+      bool GetCuttingFlag() const
+               { /*return m_bAllPartCut ;*/ return ( m_nType == GEN ? m_ArcLineApprox.GetTempProp( 1) == 1 : true ) ; }
 
    public :
       enum ToolType { UNDEF = 0,         // Utensile indefinito
@@ -70,13 +73,15 @@ class Tool
                       BULLNOSEMILL = 4,  // Naso di toro
                       CONEMILL = 5,      // Con parte terminale conica
                       MORTISER = 6,      // Mortasatrice
-                      CHISEL = 7} ;      // Scalpello
+                      CHISEL = 7,        // Scalpello
+                      ADDITIVE = 8} ;    // Additivo
 
    private :
        bool ModifyForCutterHeight( void) ;
 
    private :
       bool           m_bApproxWithLines ;
+      //bool           m_bAllPartCut ;
       double         m_dLinTol ;
       double         m_dAngTolDeg ;
       std::string    m_sName ;

Triangulate.cpp

@@ -12,17 +12,34 @@
 //----------------------------------------------------------------------------
 
 //--------------------------- Include ----------------------------------------
+//#include "stdafx.h"
+//#include "DllMain.h"
+//#include "Triangulate.h"
+//#include "ProjPlane.h"
+//#include "/EgtDev/Include/EGkPolyLine.h"
+//#include "/EgtDev/Include/EGkPlane3d.h"
+//#include "/EgtDev/Include/EGkStringUtils3d.h"
+//#include <algorithm>
 #include "stdafx.h"
 #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 "tpp_interface.hpp"
+
 #include <algorithm>
 
 using namespace std ;
 
+using namespace tpp ;
+
+
 //----------------------------------------------------------------------------
 enum EarStatus{ EAS_NULL = -1, EAS_NO = 0, EAS_OK = 1} ;
 
@@ -35,13 +52,17 @@ 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
    double dArea ;
    Plane3d plPlane ;
    if ( ! PL.IsClosedAndFlat( plPlane, dArea, 50 * EPS_SMALL))
       return false ;
+
+   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))
@@ -82,12 +103,13 @@ Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr)
 
 //----------------------------------------------------------------------------
 // In : POLYLINEVECTOR : vector of polylines, the first outer, the others inner
+//      trgType : triangulation type
 // Out : PNTVECTOR (Point3d Vector) : points of the polyline 
 //       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() ;
    vTr.clear() ;
@@ -95,7 +117,7 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
    if ( &vPL == nullptr  ||  vPL.empty())
       return false ;
   // se una sola polilinea mi riconduco al caso precedente
-   if ( vPL.size() == 1)
+   if ( vPL.size() == 1 && trgType == TRG_STANDARD)
       return Make( vPL[0], vPt, vTr) ;
   // verifico che la polilinea esterna sia chiusa e piana e calcolo il piano medio del poligono
    double dArea ;
@@ -115,6 +137,31 @@ 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 ;  
+   }
+
+  // ear clipping 
   // se non CCW inverto tutte le polilinee
    if ( ! bCCW) {
       for ( int i = 0 ; i < int( vPL.size()) ; ++i)
@@ -172,7 +219,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 +237,142 @@ 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) 
+{
+   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) ; 
+   }
+
+   // 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 = 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( vDelaunayVert.size() - 1) ;
+         vDelaunayConstr.push_back( 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 ( ! vMyPL[i].IsPointInsidePolyLine( pt)) {
+         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 ;
+}
+
+
+
 //----------------------------------------------------------------------------
 bool
 Triangulate::PrepareGrid( const PNTVECTOR& vPt, const INTVECTOR& vPol,
@@ -273,7 +456,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 +576,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 +689,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] ;
@@ -534,7 +717,10 @@ Triangulate::TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol,
    bool bIsEar = true ;
   // An ear must be convex (here counterclockwise)
    if ( TriangleIsCCW( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]])  &&
-        ! Collinear( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]])) {
+        ! Collinear( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]])  /*&&
+        ! AreSamePoint( vPt[vPol[vPrev[i]]], vPt[vPol[i]])  &&
+        ! AreSamePoint( vPt[vPol[i]], vPt[vPol[vNext[i]]])  &&
+        ! AreSamePoint( vPt[vPol[vNext[i]]], vPt[vPol[vPrev[i]]])*/  ) {
      // Loop over all vertices not part of the tentative ear
       BBox3d b3Tria ;
       b3Tria.Add( vPt[vPol[vPrev[i]]]) ;
@@ -570,7 +756,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 +1082,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 +1105,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 +1141,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 +1156,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,19 @@
 #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  
+             } ;
+
 //----------------------------------------------------------------------------
 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 +36,7 @@ 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 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) ;
@@ -50,3 +58,5 @@ class Triangulate
       PointGrid3d m_VertGrid ;
       INTVECTOR m_vVert ;
 } ;
+
+

VolZmap.cpp

@@ -35,7 +35,7 @@ GEOOBJ_REGISTER( VOL_ZMAP, NGE_V_ZMP, VolZmap) ;
 //----------------------------------------------------------------------------
 VolZmap::VolZmap(void)
    : m_nStatus( TO_VERIFY), m_dStep( 10.0), m_nMapNum( 0), m_nShape( GENERIC), m_nVoxNumPerBlock( N_VOXBLOCK),
-     m_nDexVoxRatio( 1), m_nNumBlock( 0), m_nConnectedCompoCount( 0), m_Tool( true)
+     m_nDexVoxRatio( 1), m_nNumBlock( 0), m_nConnectedCompoCount( 0), m_nCurrTool( - 1)/*, m_Tool( true)*/
 {
    for ( int i = 0 ; i < N_MAPS ; ++ i) {
       m_nNx[i] = 0 ;
@@ -47,7 +47,7 @@ VolZmap::VolZmap(void)
    }
    m_nTempProp[0] = 0 ;
    m_nTempProp[1] = 0 ;
-   m_Tool.SetTolerances( LIN_TOL_STD, ANG_TOL_APPROX_DEG) ;
+   //m_Tool.SetTolerances( LIN_TOL_STD, ANG_TOL_APPROX_DEG) ;
 }
 
 //----------------------------------------------------------------------------
@@ -76,7 +76,7 @@ VolZmap::Clear( void)
    m_dStep = EPS_SMALL ;
    m_nTempProp[0] = 0 ;
    m_nTempProp[1] = 0 ;
-   m_Tool.Clear() ;
+   ResetAllTools() ;
   // imposto ricalcolo della grafica
    m_OGrMgr.Reset() ;
 
@@ -1680,8 +1680,10 @@ VolZmap::Cut( const Plane3d& plPlane)
    Plane3d plMyPlane = plPlane ;
    plMyPlane.ToLoc( m_MapFrame) ;
   // Imposto numero fittizio di utensile per avere il colore di sezione opportuno
-   int nToolNumOld = m_Tool.GetCurrentToolNum() ;
-   m_Tool.SetToolNum( 1) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   int nToolNumOld = m_vTool[m_nCurrTool].GetCurrentToolNum() ;
+   m_vTool[m_nCurrTool].SetToolNum( 1) ;
   // Interseco lo Zmap col piano, ciclando sulle griglie
    bool bModified = false ;
    for ( int nMap = 0 ; nMap < int( m_nMapNum) ; ++ nMap) {
@@ -1735,7 +1737,7 @@ VolZmap::Cut( const Plane3d& plPlane)
    }
 
   // Ripristino numero utensile
-   m_Tool.SetToolNum( nToolNumOld) ;
+   m_vTool[m_nCurrTool].SetToolNum( nToolNumOld) ;
 
    if ( bModified == true) {
      // Imposto forma generica
@@ -1970,14 +1972,18 @@ VolZmap::CalcBlockNum( void)
 bool 
 VolZmap::SetToolTolerances( double dLinTol, double dAngTolDeg) 
 {
-   return m_Tool.SetTolerances( dLinTol, dAngTolDeg) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetTolerances( dLinTol, dAngTolDeg) ;
 }
 
 //----------------------------------------------------------------------------
 bool
 VolZmap::SetStdTool( const string& sToolName, double dH, double dR, double dCornR, double dCutterH, int nFlag) 
 {
-   return m_Tool.SetStdTool( sToolName, dH, dR, dCornR, dCutterH, nFlag) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetStdTool( sToolName, dH, dR, dCornR, dCutterH, nFlag) ;
 }
 
 //----------------------------------------------------------------------------
@@ -1985,7 +1991,9 @@ bool
 VolZmap::SetAdvTool( const string& sToolName, 
                      double dH, double dR, double dTipH, double dTipR, double dCornR, double dCutterH, int nFlag)
 {
-   return m_Tool.SetAdvTool( sToolName, dH, dR, dTipH, dTipR, dCornR, dCutterH, nFlag) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetAdvTool( sToolName, dH, dR, dTipH, dTipR, dCornR, dCutterH, nFlag) ;
 }
 
 //----------------------------------------------------------------------------
@@ -1993,33 +2001,90 @@ bool
 VolZmap::SetSawTool( const string& sToolName,
                      double dH, double dR, double dThick, double dStemR, double dCornR, int nFlag)
 {
-   return m_Tool.SetSawTool( sToolName, dH, dR, dThick, dStemR, dCornR, nFlag) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetSawTool( sToolName, dH, dR, dThick, dStemR, dCornR, nFlag) ;
 }
 
 //----------------------------------------------------------------------------
 bool 
 VolZmap::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline, int nFlag)
 {
-   return m_Tool.SetGenTool( sToolName, pToolOutline, nFlag) ; 
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetGenTool( sToolName, pToolOutline, nFlag) ; 
 }
 
 //----------------------------------------------------------------------------
 bool 
 VolZmap::SetMortiserTool( const string& sToolName, double dH, double dW, double dTh, double dRc, int nFlag)
 {
-   return m_Tool.SetMortiserTool( sToolName, dH, dW, dTh, dRc, nFlag) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetMortiserTool( sToolName, dH, dW, dTh, dRc, nFlag) ;
 }
 
 //----------------------------------------------------------------------------
 bool 
 VolZmap::SetChiselTool( const string& sToolName, double dH, double dW, double dTh, int nFlag)
 {
-   return m_Tool.SetChiselTool( sToolName, dH, dW, dTh, nFlag) ; 
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetChiselTool( sToolName, dH, dW, dTh, nFlag) ; 
+}
+
+//----------------------------------------------------------------------------
+bool
+VolZmap::SetCurrTool( int nCurrTool)
+{
+   if ( nCurrTool < 0  ||  nCurrTool >= int( m_vTool.size()))
+      return false ;
+   m_nCurrTool = nCurrTool ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+int
+VolZmap::GetToolCount( void) const
+{
+   return int( m_vTool.size()) ;
+}
+
+//----------------------------------------------------------------------------
+bool
+VolZmap::AddTool( void)
+{
+   m_vTool.emplace_back( true) ;
+   m_vTool.back().SetTolerances( LIN_TOL_STD, ANG_TOL_APPROX_DEG) ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+VolZmap::ResetAllTools( void)
+{
+   bool bOk = true ;
+   for ( int n = 0 ; n < int( m_vTool.size()) ; ++ n) {
+      bOk = bOk  ||  m_vTool[n].Clear() ;
+   }
+   m_vTool.clear() ;
+   return bOk ;
 }
 
 //----------------------------------------------------------------------------
 bool 
 VolZmap::ResetTool( void)
 {
-   return m_Tool.Clear() ; 
+   bool bOk = m_vTool.back().Clear() ;
+   m_vTool.erase( m_vTool.begin() + int( m_vTool.size()) - 1) ; 
+   return bOk ; 
+}
+
+//----------------------------------------------------------------------------
+const ICurveComposite&
+VolZmap::GetToolOutline( bool bApprox) const
+{
+   if ( m_nCurrTool < 0)
+      return cvEmptyOutline ;
+   return ( bApprox ? m_vTool[m_nCurrTool].GetApproxOutline() : m_vTool[m_nCurrTool].GetOutline()) ;
 }

VolZmap.h

@@ -72,6 +72,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool Create( const Point3d& ptO, double dDimX, double dDimY, double dDimZ, double dStep, bool bTriDex) override ;
       bool CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double dStep, bool bTriDex) override ;
       bool CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex) override ;
+      bool CreateEmptyMap( const Point3d& ptO, double dLengthX, double dLengthY, double dLengthZ, double dStep, bool bTriDex) override ;
       int  GetBlockCount( void) const override ;
       int  GetBlockUpdatingCounter( int nBlock) const override ;
       bool GetBlockTriangles( int nBlock, TRIA3DEXVECTOR& vTria) const override ;
@@ -90,9 +91,14 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool SetGenTool( const std::string& sToolName, const ICurveComposite* pToolOutline, int nFlag) override ;
       bool SetMortiserTool( const std::string& sToolName, double dH, double dW, double dTh, double dRc, int nFlag) override ;
       bool SetChiselTool( const std::string& sToolName, double dH, double dW, double dTh, int nFlag) override ;
+      bool SetCurrTool( int nCurrTool) override ;
+      int GetToolCount( void) const override ;
+      int GetCurrTool( void) const override
+             { return m_nCurrTool ; }
+      bool AddTool( void) override ;
+      bool ResetAllTools( void) override ;
       bool ResetTool( void) override ;
-      const ICurveComposite& GetToolOutline( bool bApprox = false) const override
-                  { return ( bApprox ? m_Tool.GetApproxOutline() : m_Tool.GetOutline()) ;}
+      const ICurveComposite& GetToolOutline( bool bApprox = false) const override ;
       bool MillingStep( const Point3d& ptPs, const Vector3d& vtDs, const Point3d& ptPe, const Vector3d& vtDe) override ;
       bool MillingStep( const Point3d& ptPs, const Vector3d& vtDs, const Vector3d& vtAs,
                         const Point3d& ptPe, const Vector3d& vtDe, const Vector3d& vtAe) override ;
@@ -303,6 +309,38 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       inline bool TestParaBBox( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtD, const Vector3d& vtA,
                                 double dLenX, double dLenY, double dLenZ,
                                 int& nStI, int& nStJ, int& nEnI, int& nEnJ) ;
+
+    // Asportazioni superfici elementari vuote
+      bool SurfCircCrown_ZDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, double dMaxRad, double dMinRad) ;
+      bool SurfCircCrown_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, double dMaxRad, double dMinRad) ;
+      bool SurfCircCrown_ZMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, double dMaxRad, double dMinRad) ;
+      bool SurfCircCrown_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, double dMaxRad, double dMinRad) ;
+      /*bool SurfCyl_ZDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir,
+                              double dHei, double dRad, bool bTapB, bool bTapT) ;
+      bool SurfCyl_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir,
+                             double dHei, double dRad, bool bTapB, bool bTapT) ;*/
+      bool SurfCyl_ZMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir,
+                             double dHei, double dRad, bool bOuterCutter) ;
+      bool SurfCyl_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE,
+                            const Vector3d& vtToolDir, double dHei, double dRad, bool bOuterCutter, bool bTapB, bool bTapT) ;
+      bool SurfConus_ZDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, double dHei, double dMaxRad, double dMinRad,
+                                bool bOuterCutter, const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR) ;
+      bool SurfConus_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, 
+                               double dHei, double dMaxRad, double dMinRad, bool bOuterCutter, bool bTapB, bool bTapT,
+                               const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR) ;
+      bool SurfConus_ZMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, double dHei, double dMaxRad, double dMinRad,
+                               bool bOuterCutter, const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR) ;
+      bool SurfConus_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, 
+                              double dHei, double dMaxRad, double dMinRad, bool bOuterCutter, bool bTapB, bool bTapT,
+                              const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR) ;
+      bool SurfSphericalShellPart_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx,
+                                           double dRad, double dInfH, double dSupH, bool bOuterCutter) ;
+      bool SurfSphericalShell_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx,
+                                       double dRad, double dHei, bool bOuterCutter) ;
+    // Additivo
+      bool AddingMotion( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx/*, double dHei, double dRad, double dCornerRad*/) ;
+
+
     // Intersezioni
       bool IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax) const ;
       bool IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax,
@@ -332,6 +370,34 @@ class VolZmap : public IVolZmap, public IGeoObjRW
                                        const Frame3d& frTruncPyramFrame, double dSegMin, double dSegMax, double dHeight,
                                        Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
       bool TestIntersPlaneZmapBBox( const Plane3d& plPlane) const ;
+
+     // Intersezioni per asportazioni avanzate
+      int IntersLineCircCrown( const Point3d& ptLineP, const Vector3d& vtLineDir,
+                               const Point3d& ptCen, const Vector3d& vtAx, double dMaxRad, double dMinRad,
+                               Point3d& ptInt, Vector3d& vtN) const ;
+      int IntersLineParallelogram( const Point3d& ptLineP, const Vector3d& vtLineDir,
+                                   const Point3d& ptParOrig, const Vector3d& vtSeg1, const Vector3d& vtSeg2,
+                                   bool bExtNorm, Point3d& ptInt, Vector3d& vtN) const ;
+      int IntersLineCylinderCuttedByPlanes( const Point3d& ptLineP, const Vector3d& vtLineDir,
+                                            const Point3d& ptBaseCen, const Vector3d& vtAx, double dRad, double dH, bool bInOut,
+                                            const std::vector<Plane3d>& vPlanesVec,
+                                            Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
+      int IntersLineCircSweptSurfCuttedByPlanes( const Point3d& ptLineP, const Vector3d& vtLineDir,
+                                                 const Point3d& ptCen, const Vector3d& vtAx, double dRad, const Vector3d& vtSweptVec, bool bInOut,
+                                                 const std::vector<Plane3d>& vPlanesVec,
+                                                 Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
+      int IntersLineConeCuttedByPlanes( const Point3d& ptLineP, const Vector3d& vtLineDir,
+                                        const Point3d& ptVert, const Vector3d& vtAx, double dRad, double dH, bool bInOut,
+                                        const std::vector<Plane3d>& vPlanesVec,
+                                        Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
+      int IntersLineSphereCuttedByPlanes( const Point3d& ptLineP, const Vector3d& vtLineD,
+                                          const Point3d& ptCen, double dRad, bool bInOut, 
+                                          const std::vector<Plane3d>& vPlanesVec,
+                                          Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
+      int IntersLineCatTongue( const Point3d& ptLineP, const Vector3d& vtLineD,
+                               const Point3d& ptCenSt, const Point3d& ptCenEn, const Vector3d& vtNorm, double dRad,
+                               Point3d& ptInt, Vector3d& vtN) const ;
+     
      // Voxel: esistenza e passaggio da N a ijk per i voxel
       bool IsValidVoxel( int nN) const ;
       bool IsValidVoxel( int nI, int nJ, int nK) const ;
@@ -452,7 +518,10 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       mutable std::atomic<bool>           m_bBreak ;
       std::atomic<bool>                   m_bIsBox ;
 
-      Tool m_Tool ;
+      // Utensili
+      std::vector<Tool> m_vTool ;
+      int m_nCurrTool ;
+      CurveComposite cvEmptyOutline ;
 } ;
 
 

VolZmapCreation.cpp

@@ -28,6 +28,10 @@ using namespace std ;
 bool 
 VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dLengthZ, double dStep, bool bTriDex)
 {
+   /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+   //////////////////////////////////////// Per test additivi //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+   //return CreateEmptyMap( ptO, dLengthX, dLengthY, dLengthZ, dStep, bTriDex) ;
+   /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
   // Controlli sull'ammissibilità delle dimensioni lineari del grezzo e del passo
    if ( dStep < EPS_SMALL  ||  dLengthX < EPS_SMALL  ||  dLengthY < EPS_SMALL  ||  dLengthZ < EPS_SMALL)
       return false ;
@@ -665,3 +669,77 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex
 
    return bCompleted ;
 }
+
+//----------------------------------------------------------------------------
+bool 
+VolZmap::CreateEmptyMap( const Point3d& ptO, double dLengthX, double dLengthY, double dLengthZ, double dStep, bool bTriDex)
+{
+  // Controlli sull'ammissibilità delle dimensioni lineari del grezzo e del passo
+   if ( dStep < EPS_SMALL  ||  dLengthX < EPS_SMALL  ||  dLengthY < EPS_SMALL  ||  dLengthZ < EPS_SMALL)
+      return false ;
+
+  // Il passo di discretizzazione non può essere inferiore a 100 * EPS_SMALL
+   m_dStep = max( dStep, 100 * EPS_SMALL) ;
+
+  // Aggiorno la dimensione della mappa 1 o 3
+   m_nMapNum = ( bTriDex ? 3 : 1) ;
+
+  // Disponendo i sistemi di riferimento in una successione, le coordinate x,y,z
+  // di uno si ottengono da una permutazione ciclica di quelle del precedente sistema.
+  // es: X(n) = Z(n-1), Y(n) = X(n-1), Z(n) = Y(n-1)
+
+  // Definisco il sistema di riferimento intrinseco
+   m_MapFrame.Set( ptO, X_AX, Y_AX, Z_AX) ;
+
+  // Definisco i vettori dei limiti su indici
+   m_nNx[0] = max( int( ( dLengthX + EPS_SMALL) / m_dStep + 0.5), 1) ;
+   m_nNy[0] = max( int( ( dLengthY + EPS_SMALL) / m_dStep + 0.5), 1) ;
+
+  // Numero di componenti connesse 
+   m_nConnectedCompoCount = 1 ;
+
+  // Se tridexel
+   if ( bTriDex) {
+      m_nNx[1] = m_nNy[0] ;
+      m_nNy[1] = max( int( ( dLengthZ + EPS_SMALL) / m_dStep + 0.5), 1) ;
+      m_nNx[2] = m_nNy[1] ;
+      m_nNy[2] = m_nNx[0] ;
+   }
+
+  // altrimenti mono dexel
+   else {
+      m_nNx[1] = 0 ;
+      m_nNy[1] = 0 ;
+      m_nNx[2] = 0 ;
+      m_nNy[2] = 0 ; 
+   }
+
+   // Definisco il numero di blocchi lungo x,y e z
+   if ( ! CalcBlockNum())
+      return false ;
+
+  // Creazione delle mappe
+  // Calcolo del numero di celle per ogni mappa
+   for ( int i = 0 ; i < m_nMapNum ; ++ i)
+      m_nDim[i] = m_nNx[i] * m_nNy[i] ;
+
+  // Creazione delle celle per ogni mappa
+   for ( int i = 0 ; i < m_nMapNum ; ++ i)
+      m_Values[i].resize( m_nDim[i]) ;
+
+  // Definizione delle limitazioni iniziali in Z per ogni mappa
+   m_dMinZ[0] = 0 ;
+   m_dMaxZ[0] = dLengthZ ;
+   m_dMinZ[1] = 0 ;
+   m_dMaxZ[1] = ( bTriDex ? dLengthX : 0) ;
+   m_dMinZ[2] = 0 ;  
+   m_dMaxZ[2] = ( bTriDex ? dLengthY : 0) ;
+
+   // Tipologia
+   m_nShape = BOX ;
+
+  // Aggiornamento dello stato
+   m_nStatus = OK ;
+
+   return true ;
+}

VolZmapGraphics.cpp

@@ -844,7 +844,7 @@ VolZmap::UpdateTripleMapGraphics( void) const
    VecTriHold.resize( m_nNumBlock) ;
 
   // Ciclo sui blocchi per eliminare le slice fra blocchi da aggiornare
-   for ( int t = 0 ; t < m_nNumBlock ; ++ t) {
+   for ( int t = 0 ; t < m_nNumBlock ; ++ t) { 
       for ( auto it = m_SliceXY[t].begin() ; it != m_SliceXY[t].end() ;) {
          int nSlIJK[3] ;
          if ( GetVoxIJKFromN( it->first, nSlIJK[0], nSlIJK[1], nSlIJK[2])) {
@@ -1083,7 +1083,7 @@ VolZmap::ExtMarchingCubes( int nBlock, VoxelContainer& vVox) const
 
            // Flag di regolarità dei campi scalare e vettoriale
             bool bReg = true ;
-
+            
            // Ciclo sui segmenti
             for ( int EdgeIndex = 0 ; EdgeIndex < 12 ; ++ EdgeIndex) {
               // Se il segmento non attraversa la superficie passo al successivo

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 */
+
+