From a4609dfa0f17665fb714e8785b0bf8570d5d8135 Mon Sep 17 00:00:00 2001 From: Dario Sassi Date: Wed, 2 Oct 2019 06:27:18 +0000 Subject: [PATCH] EgtGeomKernel 2.1.10.1 : - nuova versione delle operazioni booleane sulle superfici trimesh. --- EgtGeomKernel.rc | Bin 11710 -> 11718 bytes GdbExecutor.cpp | 215 +++-- GdbExecutor.h | 13 +- IntersTriaTria.cpp | 27 +- SurfTriMesh.h | 20 + SurfTriMeshBooleans.cpp | 1863 +++++++++++++++++++++++++++++++-------- VolZmap.cpp | 107 +-- VolZmapVolume.cpp | 4 +- 8 files changed, 1751 insertions(+), 498 deletions(-) diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index e0c4a5cda281d7a211312473ea42cf05aa4318ad..64bc7b98f5265da45c71d3dd35c1d4fdf74155c0 100644 GIT binary patch delta 126 zcmdlNeJpyzH#T-d1_K5i2E)z2+2om-a~KRKD{|{@_TuVcfy>Wy)SyhH#SyF1|0_D&G-4vGfl1&(wY2&Q-wK`!FaNwuGeneralizedCut( *pCurve, true) ; } +//---------------------------------------------------------------------------- +bool +GdbExecutor::SurfTrimeshAdd( const STRVECTOR& vsParams) +{ + // Parametri: nSurfId, nOtherSurfId + if ( vsParams.size() != 2) + return false ; + // Recupero Id superficie + int nIdSurf = GetIdParam( vsParams[0]) ; + ISurfTriMesh* pStmM = GetSurfTriMesh( m_pGDB->GetGeoObj( nIdSurf)) ; + SurfTriMesh* pStmF = static_cast ( pStmM) ; + int nIdOtherSurf = GetIdParam( vsParams[1]) ; + ISurfTriMesh* pStmOtherM = GetSurfTriMesh( m_pGDB->GetGeoObj( nIdOtherSurf)) ; + SurfTriMesh* pStmOtherF = static_cast ( pStmOtherM) ; + return pStmF->Add( *pStmOtherF) ; +} + +//---------------------------------------------------------------------------- +bool +GdbExecutor::SurfTrimeshInters( const STRVECTOR& vsParams) +{ + // Parametri: nSurfId, nOtherSurfId + if ( vsParams.size() != 2) + return false ; + // Recupero Id superficie + int nIdSurf = GetIdParam( vsParams[0]) ; + ISurfTriMesh* pStmM = GetSurfTriMesh( m_pGDB->GetGeoObj( nIdSurf)) ; + SurfTriMesh* pStmF = static_cast ( pStmM) ; + int nIdOtherSurf = GetIdParam( vsParams[1]) ; + ISurfTriMesh* pStmOtherM = GetSurfTriMesh( m_pGDB->GetGeoObj( nIdOtherSurf)) ; + SurfTriMesh* pStmOtherF = static_cast ( pStmOtherM) ; + return pStmF->Intersect( *pStmOtherF) ; +} + +//---------------------------------------------------------------------------- +bool +GdbExecutor::SurfTrimeshSubtract( const STRVECTOR& vsParams) +{ + // Parametri: nSurfId, nOtherSurfId + if ( vsParams.size() != 2) + return false ; + // Recupero Id superficie + int nIdSurf = GetIdParam( vsParams[0]) ; + ISurfTriMesh* pStmM = GetSurfTriMesh( m_pGDB->GetGeoObj( nIdSurf)) ; + SurfTriMesh* pStmF = static_cast ( pStmM) ; + int nIdOtherSurf = GetIdParam( vsParams[1]) ; + SurfTriMesh* pStmOther = GetBasicSurfTriMesh( m_pGDB->GetGeoObj( nIdOtherSurf)) ; + if ( pStmOther == nullptr) + return false ; + return pStmF->Subtract( *pStmOther) ; +} + //---------------------------------------------------------------------------- bool GdbExecutor::ExecuteVolZmap( const string& sCmd2, const STRVECTOR& vsParams) @@ -2705,15 +2766,15 @@ GdbExecutor::ExecuteVolZmap( const string& sCmd2, const STRVECTOR& vsParams) //else if ( sCmd2 == "TRIMESHZMAP") { // return VolZmapCreateFromTriMesh( vsParams) ; //} - if ( sCmd2 == "CUT") { - return ExecuteVolZmapCut( vsParams) ; - } - else if ( sCmd2 == "COMP") { - return ExecuteVolZmapCompact( vsParams) ; - } - else if ( sCmd2 == "EDGES") { - return ExecuteVolZmapEdges( vsParams) ; - } + //if ( sCmd2 == "CUT") { + // return ExecuteVolZmapCut( vsParams) ; + //} + //else if ( sCmd2 == "COMP") { + // return ExecuteVolZmapCompact( vsParams) ; + //} + //else if ( sCmd2 == "EDGES") { + // return ExecuteVolZmapEdges( vsParams) ; + //} return false ; } /* @@ -3024,76 +3085,76 @@ GdbExecutor::VolZmapBBoxZmapIntersection( const STRVECTOR& vsParams) }*/ //---------------------------------------------------------------------------- -bool -GdbExecutor::ExecuteVolZmapCut( const STRVECTOR& vsParams) -{ - // Parametri : ZmapId, ParentId, ptP, vtN - if ( vsParams.size() != 4) - return false ; - // Recupero Zmap - int nZmapId = GetIdParam( vsParams[0]) ; - VolZmap* pZmap = GetBasicVolZmap( m_pGDB->GetGeoObj( nZmapId)) ; - if ( pZmap == nullptr) - return false ; - // recupero il riferimento in cui è immerso lo Zmap - Frame3d frRef ; - if ( ! m_pGDB->GetGroupGlobFrame( GetIdParam( vsParams[1]), frRef)) - return false ; - // Recupero punto piano - Point3d ptP ; - if ( ! GetPointParam( vsParams[2], frRef, ptP)) - return false ; - // primo vettore frBBoxFrame - Vector3d vtN ; - if ( ! GetVectorParam( vsParams[3], frRef, vtN)) - return false ; - vtN.Normalize() ; - Plane3d plPlane ; - plPlane.Set( ptP, vtN) ; - pZmap->Cut( plPlane) ; - return true ; -} +//bool +//GdbExecutor::ExecuteVolZmapCut( const STRVECTOR& vsParams) +//{ +// // Parametri : ZmapId, ParentId, ptP, vtN +// if ( vsParams.size() != 4) +// return false ; +// // Recupero Zmap +// int nZmapId = GetIdParam( vsParams[0]) ; +// VolZmap* pZmap = GetBasicVolZmap( m_pGDB->GetGeoObj( nZmapId)) ; +// if ( pZmap == nullptr) +// return false ; +// // recupero il riferimento in cui è immerso lo Zmap +// Frame3d frRef ; +// if ( ! m_pGDB->GetGroupGlobFrame( GetIdParam( vsParams[1]), frRef)) +// return false ; +// // Recupero punto piano +// Point3d ptP ; +// if ( ! GetPointParam( vsParams[2], frRef, ptP)) +// return false ; +// // primo vettore frBBoxFrame +// Vector3d vtN ; +// if ( ! GetVectorParam( vsParams[3], frRef, vtN)) +// return false ; +// vtN.Normalize() ; +// Plane3d plPlane ; +// plPlane.Set( ptP, vtN) ; +// pZmap->Cut( plPlane) ; +// return true ; +//} //---------------------------------------------------------------------------- -bool -GdbExecutor::ExecuteVolZmapCompact( const STRVECTOR& vsParams) -{ - // Parametri : ZmapId - if ( vsParams.size() != 1) - return false ; - // Recupero Zmap - int nZmapId = GetIdParam( vsParams[0]) ; - VolZmap* pZmap = GetBasicVolZmap( m_pGDB->GetGeoObj( nZmapId)) ; - if ( pZmap == nullptr) - return false ; - pZmap->Compact() ; - return true ; -} +//bool +//GdbExecutor::ExecuteVolZmapCompact( const STRVECTOR& vsParams) +//{ +// // Parametri : ZmapId +// if ( vsParams.size() != 1) +// return false ; +// // Recupero Zmap +// int nZmapId = GetIdParam( vsParams[0]) ; +// VolZmap* pZmap = GetBasicVolZmap( m_pGDB->GetGeoObj( nZmapId)) ; +// if ( pZmap == nullptr) +// return false ; +// pZmap->Compact() ; +// return true ; +//} //---------------------------------------------------------------------------- -bool -GdbExecutor::ExecuteVolZmapEdges(const STRVECTOR& vsParams) -{ - // Parametri : ZmapId, ParentId - if ( vsParams.size() != 2) - return false ; - // Recupero Zmap - int nZmapId = GetIdParam( vsParams[0]) ; - VolZmap* pZmap = GetBasicVolZmap( m_pGDB->GetGeoObj( nZmapId)) ; - if ( pZmap == nullptr) - return false; - // recupero il riferimento in cui è immerso lo Zmap - Frame3d frRef ; - if ( ! m_pGDB->GetGroupGlobFrame( GetIdParam( vsParams[1]), frRef)) - return false ; - ICURVEPOVECTOR vpLoop ; - pZmap->GetEdges( vpLoop) ; - bool bOk = true ; - for ( int n = 0 ; n < int( vpLoop.size()) && bOk ; ++ n) { - bOk = bOk && AddGeoObj( "$NN", vsParams[1], Release( vpLoop[n])) ; - } - return bOk ; -} +//bool +//GdbExecutor::ExecuteVolZmapEdges(const STRVECTOR& vsParams) +//{ +// // Parametri : ZmapId, ParentId +// if ( vsParams.size() != 2) +// return false ; +// // Recupero Zmap +// int nZmapId = GetIdParam( vsParams[0]) ; +// VolZmap* pZmap = GetBasicVolZmap( m_pGDB->GetGeoObj( nZmapId)) ; +// if ( pZmap == nullptr) +// return false; +// // recupero il riferimento in cui è immerso lo Zmap +// Frame3d frRef ; +// if ( ! m_pGDB->GetGroupGlobFrame( GetIdParam( vsParams[1]), frRef)) +// return false ; +// ICURVEPOVECTOR vpLoop ; +// pZmap->GetEdges( vpLoop) ; +// bool bOk = true ; +// for ( int n = 0 ; n < int( vpLoop.size()) && bOk ; ++ n) { +// bOk = bOk && AddGeoObj( "$NN", vsParams[1], Release( vpLoop[n])) ; +// } +// return bOk ; +//} //---------------------------------------------------------------------------- bool diff --git a/GdbExecutor.h b/GdbExecutor.h index bd6d0a7..7447b5f 100644 --- a/GdbExecutor.h +++ b/GdbExecutor.h @@ -1,7 +1,7 @@ //---------------------------------------------------------------------------- -// EgalTech 2013-2013 +// EgalTech 2013-2019 //---------------------------------------------------------------------------- -// File : GdbExecutor.h Data : 25.11.13 Versione : 1.3a1 +// File : GdbExecutor.h Data : 01.10.19 Versione : 2.1i3 // Contenuto : Dichiarazione della classe GdbExecutor. // // @@ -118,6 +118,9 @@ class GdbExecutor : public IGdbExecutor bool SurfTrimeshRemovePart( const STRVECTOR& vsParams) ; bool SurfTrimeshClonePart( const STRVECTOR& vsParams) ; bool SurfTrimeshGenCut( const STRVECTOR& vsParams) ; + bool SurfTrimeshAdd( const STRVECTOR& vsParams) ; + bool SurfTrimeshInters( const STRVECTOR& vsParams) ; + bool SurfTrimeshSubtract( const STRVECTOR& vsParams) ; bool ExecuteVolZmap( const std::string& sCmd2, const STRVECTOR& vsParams) ; //bool VolZmapCreate( const STRVECTOR& vsParams) ; //bool VolZmapCreateFromFlatRegion( const STRVECTOR& vsParams) ; @@ -127,9 +130,9 @@ class GdbExecutor : public IGdbExecutor //bool VolZmapSetGenTool( const STRVECTOR& vsParams) ; //bool VolZmapDeepnessMeasure( const STRVECTOR& vsParams) ; //bool VolZmapBBoxZmapIntersection( const STRVECTOR& vsParams) ; - bool ExecuteVolZmapCut( const STRVECTOR& vsParams) ; - bool ExecuteVolZmapCompact( const STRVECTOR& vsParams) ; - bool ExecuteVolZmapEdges( const STRVECTOR& vsParams) ; + //bool ExecuteVolZmapCut( const STRVECTOR& vsParams) ; + //bool ExecuteVolZmapCompact( const STRVECTOR& vsParams) ; + //bool ExecuteVolZmapEdges( const STRVECTOR& vsParams) ; bool ExecuteIntersection( const std::string& sCmd2, const STRVECTOR& vsParams) ; bool LineDiscInters( const STRVECTOR& vsParams) ; bool RayDiscInters( const STRVECTOR& vsParams) ; diff --git a/IntersTriaTria.cpp b/IntersTriaTria.cpp index 3aa70be..8cc3e92 100644 --- a/IntersTriaTria.cpp +++ b/IntersTriaTria.cpp @@ -25,6 +25,14 @@ using namespace std ; +//---------------------------------------------------------------------------- +int nTriaTriaIntersCases[6][6] = { { ITTTS_NO, ITTTS_NO, ITTTS_NO, ITTTS_NO, ITTTS_NO, ITTTS_NO }, + { ITTTS_NO, ITTTS_INT_INT_SEG, ITTTS_INT_EDGE, ITTTS_NO, ITTTS_NO, ITTTS_NO }, + { ITTTS_NO, ITTTS_EDGE_INT, ITTTS_EDGE_EDGE_SEG, ITTTS_NO, ITTTS_NO, ITTTS_NO }, + { ITTTS_NO, ITTTS_NO, ITTTS_NO, ITTTS_VERT_VERT, ITTTS_VERT_EDGE, ITTTS_VERT_INT }, + { ITTTS_NO, ITTTS_NO, ITTTS_NO, ITTTS_EDGE_VERT, ITTTS_EDGE_EDGE_PNT, ITTTS_NO }, + { ITTTS_NO, ITTTS_NO, ITTTS_NO, ITTTS_INT_VERT, ITTTS_NO, ITTTS_NO } } ; + //---------------------------------------------------------------------------- static int IntersCoplanarTriaTria( const Triangle3d& trTria1, const Triangle3d& trTria2, TRIA3DVECTOR& vTria) ; @@ -81,15 +89,6 @@ IntersTriaTria( const Triangle3d& trTria1, const Triangle3d& trTria2, Point3d& p // limito la linea di intersezione con il primo triangolo int nRes1 = IntersCoplanarLineTria( ptL1, vtL1, 100.0, trTria1, ptInt, ptInt2, false) ; - switch ( nRes1) { - case ILTT_SEGM : - case ILTT_SEGM_ON_EDGE : - break ; - case ILTT_VERT : - return ( PointInTria( ptInt, trTria2) ? ITTT_VERT : ITTT_NO) ; - default : - return ITTT_NO ; - } // il segmento calcolato va limitato col secondo triangolo Point3d ptL2 = ptInt ; @@ -97,14 +96,8 @@ IntersTriaTria( const Triangle3d& trTria1, const Triangle3d& trTria2, Point3d& p double dLen = vtL2.Len() ; vtL2 /= dLen ; int nRes2 = IntersCoplanarLineTria( ptL2, vtL2, dLen, trTria2, ptInt, ptInt2, true) ; - switch( nRes2) { - case ILTT_NO : return ITTT_NO ; - case ILTT_SEGM : return ITTT_YES ; - case ILTT_SEGM_ON_EDGE : return ITTT_EDGE ; - case ILTT_VERT : return ITTT_VERT ; - case ILTT_EDGE : return ITTT_PNT ; - default : return ITTT_NO ; - } + + return nTriaTriaIntersCases[nRes1][nRes2] ; } //---------------------------------------------------------------------------- diff --git a/SurfTriMesh.h b/SurfTriMesh.h index 6f1d6ba..91530e0 100644 --- a/SurfTriMesh.h +++ b/SurfTriMesh.h @@ -59,6 +59,7 @@ class StmTria int nEFlag ; mutable int nPart ; mutable int nTemp ; + mutable int nTempPart ; } ; //---------------------------------------------------------------------------- @@ -68,6 +69,19 @@ struct Stm3Int { int nI1 ; int nI2 ; int nI3 ; } ; +//---------------------------------------------------------------------------- +// Struttura segmento di intersezione +struct IntSegment { + Point3d ptSt ; + Point3d ptEn ; + Vector3d vtOuter ; + bool bDegenerate ; +} ; +// Tipo chain +typedef std::vector Chain ; +// Tipo vettore di Chain +typedef std::vector CHAINVECTOR ; + //---------------------------------------------------------------------------- //---------------------------------------------------------------------------- class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW @@ -189,6 +203,9 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW bool GetFacetBBox( int nF, const Frame3d& frRef, BBox3d& b3Ref, int nFlag = BBF_STANDARD) const override ; bool Cut( const Plane3d& plPlane, bool bSaveOnEq) override ; bool GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) override ; + bool Add( const ISurfTriMesh& Other) override ; + bool Intersect( const ISurfTriMesh& Other) override ; + bool Subtract( const ISurfTriMesh& Other) override ; bool GetAllTriaOverlapBox( const BBox3d& b3Box, INTVECTOR& vT) const override ; const BBox3d& GetAllTriaBox( void) const override ; int GetPartCount( void) const override ; @@ -255,6 +272,9 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW void ResetHashGrids3d( void) const ; bool VerifyHashGrids3d( void) const ; bool VerifyConnection( void) const ; + bool DecomposeLoop( CHAINVECTOR& cvOpenChain, INTVECTOR& vnDegVec, std::vector& cvBoundClosedLoopVec, std::vector& vbInOut) ; + bool IntersectTriMeshTriangle( SurfTriMesh& Other) ; + bool IdentifyParts( void) const ; private : ObjGraphicsMgr m_OGrMgr ; // gestore grafica dell'oggetto diff --git a/SurfTriMeshBooleans.cpp b/SurfTriMeshBooleans.cpp index 33ebac6..96298a1 100644 --- a/SurfTriMeshBooleans.cpp +++ b/SurfTriMeshBooleans.cpp @@ -21,27 +21,17 @@ #include "GeoConst.h" #include "/EgtDev/Include/EgkCurve.h" #include "/EgtDev/Include/EgkDistPointCurve.h" +#include "/EgtDev/Include/EgkDistPointTria.h" #include "/EgtDev/Include/EgkIntersLineTria.h" #include "/EgtDev/Include/EgkIntersTriaTria.h" #include "/EgtDev/Include/EGkChainCurves.h" #include "/EgtDev/Include/EGkGeoCollection.h" +#include "/EgtDev/Include/EGkPolygon3d.h" #include +#include using namespace std ; -//---------------------------------------------------------------------------- -// Struttura segmento di intersezione -struct IntSegment { - Point3d ptSt ; - Point3d ptEn ; - Vector3d vtOuter ; - bool bDegenerate ; -} ; -// Tipo chain -typedef vector Chain ; -// Tipo vettore di Chain -typedef vector CHAINVECTOR ; - //---------------------------------------------------------------------------- int IntersRectangleTriangle( const Point3d& ptP, const Vector3d& vtL1, const Vector3d& vtL2, @@ -61,9 +51,9 @@ IntersRectangleTriangle( const Point3d& ptP, const Vector3d& vtL1, const Vector3 Point3d ptIntA1, ptIntA2 ; TRIA3DVECTOR vTriaA ; int nIntTypeA = IntersTriaTria( trTria, trTriaA, ptIntA1, ptIntA2, vTriaA) ; - if ( nIntTypeA == ITTT_PNT || nIntTypeA == ITTT_VERT) + if ( nFromSpecialToNormal[nIntTypeA] == ITTT_PNT || nFromSpecialToNormal[nIntTypeA] == ITTT_VERT) nIntA = 1 ; - else if ( nIntTypeA == ITTT_YES || nIntTypeA == ITTT_EDGE) { + else if ( nFromSpecialToNormal[nIntTypeA] == ITTT_YES || nFromSpecialToNormal[nIntTypeA] == ITTT_EDGE) { nIntA = 2 ; } // Interseco il triangolo con il secondo dei due triangoli del rettangolo @@ -71,9 +61,9 @@ IntersRectangleTriangle( const Point3d& ptP, const Vector3d& vtL1, const Vector3 Point3d ptIntB1, ptIntB2 ; TRIA3DVECTOR vTriaB ; int nIntTypeB = IntersTriaTria( trTria, trTriaB, ptIntB1, ptIntB2, vTriaB) ; - if ( nIntTypeB == ITTT_PNT || nIntTypeB == ITTT_VERT) + if ( nFromSpecialToNormal[nIntTypeB] == ITTT_PNT || nFromSpecialToNormal[nIntTypeB] == ITTT_VERT) nIntB = 1 ; - else if ( nIntTypeB == ITTT_YES || nIntTypeB == ITTT_EDGE) { + else if ( nFromSpecialToNormal[nIntTypeB] == ITTT_YES || nFromSpecialToNormal[nIntTypeB] == ITTT_EDGE) { nIntB = 2 ; } // Unisco le due intersezioni @@ -117,15 +107,171 @@ IntersRectangleTriangle( const Point3d& ptP, const Vector3d& vtL1, const Vector3 } else if ( nIntTot == 1) { if ( nIntA == 1) - ptStSeg = ptIntA1 ; + ptStSeg = ptIntA1 ; else - ptStSeg = ptIntB1 ; + ptStSeg = ptIntB1 ; return 1 ; } else return 0 ; } +//---------------------------------------------------------------------------- +bool +ChangeStart( const Point3d& ptNewStart, PNTVECTOR& Loop) +{ + // Cerco se esiste un tratto del loop chiuso su cui giace il punto + int nSeg = - 1 ; + for ( int nPt = 0 ; nPt < int( Loop.size()) && nSeg == - 1 ; ++ nPt) { + // Estremi del segmento corrente del loop + Point3d ptSegSt = Loop[nPt] ; + Point3d ptSegEn = Loop[( nPt + 1) % int( Loop.size())] ; + // Vedo se il punto giace sul segmento del loop + DistPointLine dDistCalc( ptNewStart, ptSegSt, ptSegEn) ; + double dSqDist ; + dDistCalc.GetSqDist( dSqDist) ; + if ( dSqDist < SQ_EPS_SMALL) { + nSeg = nPt ; + } + } + // Se il punto non sta sul loop, errore + if ( nSeg == - 1) + return false ; + // Verifico che il punto stia su un vertice, in tal caso non devo fare nulla + bool bOnStart = AreSamePointApprox( Loop[nSeg], ptNewStart) ; + bool bOnEnd = AreSamePointApprox( Loop[( nSeg + 1) % int( Loop.size())], ptNewStart) ; + if ( bOnStart || bOnEnd) { + if ( bOnEnd) { + ++ nSeg ; + if ( nSeg % int( Loop.size()) == 0) + return true ; + } + PNTVECTOR vTempVec ; + for ( int nPt = 0 ; nPt < nSeg ; ++ nPt) + vTempVec.emplace_back( Loop[nPt]) ; + int nSize = int( Loop.size()) ; + for ( int nPt = 0 ; nPt < nSize - nSeg ; ++ nPt) { + Loop[nPt] = Loop[nPt + nSeg] ; + } + for ( int nPt = 0 ; nPt < int( vTempVec.size()) ; ++ nPt) { + Loop[nPt + nSize - nSeg] = vTempVec[nPt] ; + } + return true ; + } + // Ridimensiono il loop + Loop.resize( Loop.size() + 1) ; + // Copio i primi punti + PNTVECTOR LoopTemp ; + for ( int nPt = 0 ; nPt <= nSeg ; ++ nPt) + LoopTemp.emplace_back( Loop[nPt]) ; + // Aggiungo il nuovo punto all'inizio + Loop[0] = ptNewStart ; + // Sposto gli ultimi in testa + int nLastPointNum = int( Loop.size()) - 1 - nSeg ; + for ( int nPt = 1 ; nPt <= nLastPointNum ; ++ nPt) { + Loop[nPt] = Loop[nPt + nSeg] ; + } + // Porto i primi in fondo + for ( int nPt = 0 ; nPt < int( LoopTemp.size()) ; ++ nPt) { + Loop[nPt + nLastPointNum] = LoopTemp[nPt] ; + } + return true ; +} + +//---------------------------------------------------------------------------- +bool +SplitAtPoint( const Point3d& ptStop, const PNTVECTOR& Loop, PNTVECTOR& Loop1, PNTVECTOR& Loop2) +{ + // Cerco se esiste un tratto del loop chiuso su cui giace il punto + int nSeg = - 1 ; + for ( int nPt = 0 ; nPt < int( Loop.size()) && nSeg == - 1 ; ++ nPt) { + // Estremi del segmento corrente del loop + Point3d ptSegSt = Loop[nPt] ; + Point3d ptSegEn = Loop[( nPt + 1) % int( Loop.size())] ; + // Vedo se il punto giace sul segmento del loop + DistPointLine dDistCalc( ptStop, ptSegSt, ptSegEn) ; + double dSqDist ; + dDistCalc.GetSqDist( dSqDist) ; + if ( dSqDist < SQ_EPS_SMALL) { + nSeg = nPt ; + } + } + // Se il punto non sta sul loop, errore + if ( nSeg == - 1) + return false ; + // Verifico che il punto stia su un vertice, in tal caso non devo aggiungerlo + bool bFirst = AreSamePointApprox( Loop[nSeg], ptStop) ; + bool bLast = AreSamePointApprox( Loop[( nSeg + 1) % int( Loop.size())], ptStop) ; + // Se il punto è sul vertice finale del segmento, aggiungo il vertice alla lista da inglobare al primo loop + if ( bLast) + ++ nSeg ; + // Inglobo fino a nSeg nel primo loop + for ( int nPt = 0 ; nPt <= nSeg ; ++ nPt) + Loop1.emplace_back( Loop[nPt]) ; + // Se il punto è interno al segmento, lo inglobo in entrambi i loop + if ( ! ( bFirst || bLast)) { + Loop1.emplace_back( ptStop) ; + Loop2.emplace_back( ptStop) ; + } + else { + Loop2.emplace_back( Loop[nSeg]) ; + } + // Inglobo gli ultimi vertici in Loop2 + for ( int nPt = nSeg + 1 ; nPt < int( Loop.size()) ; ++ nPt) + Loop2.emplace_back( Loop[nPt]) ; + Loop2.emplace_back( Loop[0]) ; + + return true ; +} + +//---------------------------------------------------------------------------- +bool +AddChainToChain( const Chain& ChainToAdd, PNTVECTOR& OrigChain) +{ + // Se la catena da aggiungere è vuota, non devo fare alcunchè + if ( ChainToAdd.size() == 0) + return true ; + // Se la catena originale è vuota, non è possibile aggiungere nulla + if ( OrigChain.size() == 0) + return false ; + // Se la catena originale è chiusa non posso aggiungere nulla + int nLastOrig = max( int( OrigChain.size()) - 1, 0) ; + if ( AreSamePointApprox( OrigChain[0], OrigChain[nLastOrig])) + return false ; + int nLastToAdd = max( int( ChainToAdd.size()) - 1, 0) ; + if ( AreSamePointApprox( OrigChain[nLastOrig], ChainToAdd[0].ptSt)) { + for ( int nPt = 1 ; nPt <= nLastToAdd ; ++ nPt) { + if ( nPt == nLastToAdd) { + if ( ! AreSamePointApprox(OrigChain[0], ChainToAdd[nPt].ptSt)) + OrigChain.emplace_back( ChainToAdd[nPt].ptSt) ; + } + else if ( nPt == 1) { + if ( ! AreSamePointApprox( OrigChain[nLastOrig], ChainToAdd[nPt].ptSt)) + OrigChain.emplace_back( ChainToAdd[nPt].ptSt) ; + } + else + OrigChain.emplace_back( ChainToAdd[nPt].ptSt) ; + } + return true ; + } + /*else if ( AreSamePointApprox( OrigChain[0], ChainToAdd[nLastToAdd].ptEn)) { + int nOldDim = OrigChain.size() ; + int nNewDim = nOldDim + ChainToAdd.size() + 1 ; + OrigChain.resize( nNewDim) ; + int nShift = std::max( nNewDim - nOldDim, 0) ; + for ( int nPt = nOldDim - 1 ; nPt >= 0 ; -- nPt) { + OrigChain[nPt + nShift] = OrigChain[nPt] ; + } + for ( int nPt = 0 ; nPt < int( ChainToAdd.size()) ; ++ nPt) { + OrigChain[nPt] = ChainToAdd[nPt].ptSt ; + } + OrigChain[int( ChainToAdd.size())] = ChainToAdd[int(ChainToAdd.size()) - 1].ptEn ; + return true ; + }*/ + else + return false ; +} + //---------------------------------------------------------------------------- bool SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) @@ -173,9 +319,7 @@ SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) // Recupero il triangolo Triangle3d trTria ; if ( ! GetTriangle( nT, trTria)) - continue ; - // Ne calcolo una copia proiettata sul piano della curva - + continue ; // Box del triangolo nel riferimento locale della curva BBox3d b3Tria ; trTria.GetLocalBBox( b3Tria) ; @@ -396,10 +540,7 @@ SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) PNTVECTOR cvFirstLoop ; cvFirstLoop.emplace_back( trTria.GetP( 0)) ; cvFirstLoop.emplace_back( trTria.GetP( 1)) ; - cvFirstLoop.emplace_back( trTria.GetP( 1)) ; cvFirstLoop.emplace_back( trTria.GetP( 2)) ; - cvFirstLoop.emplace_back( trTria.GetP( 2)) ; - cvFirstLoop.emplace_back( trTria.GetP( 0)) ; vector cvBoundClosedLoopVec ; cvBoundClosedLoopVec.emplace_back(cvFirstLoop); @@ -409,291 +550,46 @@ SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) while ( cvOpenChain.size() > 0) { int nLastOpenLoopN = int( cvOpenChain.size()) - 1 ; if ( vnDegVec[nLastOpenLoopN] == 1) { - for ( int nLoop = 0 ; nLoop < int( cvBoundClosedLoopVec.size()) ; ++ nLoop) { // Estremi del loop aperto int nLastOpenLoopPoint = max( int( cvOpenChain[nLastOpenLoopN].size()) - 1, 0) ; Point3d ptOpenLoopStP = cvOpenChain[nLastOpenLoopN][0].ptSt ; Point3d ptOpenLoopEnP = cvOpenChain[nLastOpenLoopN][nLastOpenLoopPoint].ptEn ; - // Cerco se esistono dei tratti del loop chiuso corrente che sono - // toccati dagli estremi del loop aperto corrente - int nCvSt = -1 ; - int nCvEn = -1 ; - for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size()) - 1 && ( nCvSt == -1 || nCvEn == -1) ; nLine += 2) { - // Estremi del segmento corrente del loop chiuso corrente - Point3d ptSegSt = cvBoundClosedLoopVec[nLoop][nLine] ; - Point3d ptSegEn = cvBoundClosedLoopVec[nLoop][nLine + 1] ; - // Vedo se gli estremi del loop aperto stanno su un segmento del chiuso - DistPointLine dStDistCalc( ptOpenLoopStP, ptSegSt, ptSegEn) ; - DistPointLine dEnDistCalc( ptOpenLoopEnP, ptSegSt, ptSegEn) ; - double dSqDistSt ; - dStDistCalc.GetSqDist(dSqDistSt) ; - if ( dSqDistSt < SQ_EPS_SMALL) { - nCvSt = nLine ; - } - double dSqDistEn ; - dEnDistCalc.GetSqDist(dSqDistEn) ; - if ( dSqDistEn < SQ_EPS_SMALL) { - nCvEn = nLine ; - } + PNTVECTOR Loop1, Loop2 ; + bool bChangedStart = ChangeStart( ptOpenLoopStP, cvBoundClosedLoopVec[nLoop]) ; + bool bSplitted = SplitAtPoint( ptOpenLoopEnP, cvBoundClosedLoopVec[nLoop], Loop1, Loop2) ; + if ( ! ( bChangedStart && bSplitted)) + continue ; + Chain cvCounterChain ; + for ( int nPt = int( cvOpenChain[nLastOpenLoopN].size()) - 1 ; nPt >= 0 ; -- nPt) { + IntSegment CurSeg ; + CurSeg.ptSt = cvOpenChain[nLastOpenLoopN][nPt].ptEn ; + CurSeg.ptEn = cvOpenChain[nLastOpenLoopN][nPt].ptSt ; + CurSeg.vtOuter = - cvOpenChain[nLastOpenLoopN][nPt].vtOuter ; + CurSeg.bDegenerate = cvOpenChain[nLastOpenLoopN][nPt].bDegenerate ; + cvCounterChain.emplace_back( CurSeg) ; } - // Se entrambi gli estremi del loop aperto sono su un segmento del loop chiuso devo dividere il loop in due - if ( nCvSt != - 1 && nCvEn != - 1) { - // Entrambi gli estremi del loop aperto sono su uno stesso segmento del loop chiuso - if ( nCvSt == nCvEn) { - bool bFirstInside ; - bool bFirstSt ; - PNTVECTOR cvSplitLoop1, cvSplitLoop2 ; - // Creo primo loop - for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size()) - 1 ; nLine += 2) { - // Segmenti loop chiuso - if ( nLine != nCvSt && nLine != nCvEn) { - cvSplitLoop1.emplace_back( cvBoundClosedLoopVec[nLoop][nLine]) ; - cvSplitLoop1.emplace_back( cvBoundClosedLoopVec[nLoop][nLine+1]) ; - } - // Dal chiuso all'aperto - else if ( nLine == nCvSt) { - // Distanze degli estremi del loop aperto dal punto iniziale del segmento corrente - Point3d ptSegSt, ptSegEn ; - ptSegSt = cvBoundClosedLoopVec[nLoop][nLine] ; - double dDistStSt = SqDist( ptSegSt, ptOpenLoopStP) ; - double dDistStEn = SqDist( ptSegSt, ptOpenLoopEnP) ; - // Devo percorrere il loop aperto nel suo verso - if ( dDistStSt < dDistStEn) { - // Dall'inizio del segmento corrente fino all'inizio del loop aperto - ptSegEn = ptOpenLoopStP ; - cvSplitLoop1.emplace_back( ptSegSt) ; - cvSplitLoop1.emplace_back( ptSegEn) ; - // Loop aperto - for ( int nOpenLine = 0 ; nOpenLine < int( cvOpenChain[nLastOpenLoopN].size()) ; ++ nOpenLine) { - Point3d ptOpenSt = cvOpenChain[nLastOpenLoopN][nOpenLine].ptSt ; - Point3d ptOpenEn = cvOpenChain[nLastOpenLoopN][nOpenLine].ptEn ; - cvSplitLoop1.emplace_back( ptOpenSt) ; - cvSplitLoop1.emplace_back( ptOpenEn) ; - if ( nOpenLine == 0) { - Vector3d vtLast = ptSegSt - ptSegEn ; - vtLast.Normalize() ; - bFirstInside = vtLast * cvOpenChain[nLastOpenLoopN][nOpenLine].vtOuter < 0. ; - } - } - // Dalla fine del loop aperto alla fine del segmento corrente - ptSegSt = ptOpenLoopEnP ; - ptSegEn = cvBoundClosedLoopVec[nLoop][nLine+1] ; - cvSplitLoop1.emplace_back( ptSegSt) ; - cvSplitLoop1.emplace_back( ptSegEn) ; - bFirstSt = true ; - } - // Devo percorrere il loop aperto contro il suo verso - else { - // Dall'inizio del segmento corrente fino alla fine del loop aperto - ptSegEn = ptOpenLoopEnP ; - cvSplitLoop1.emplace_back( ptSegSt) ; - cvSplitLoop1.emplace_back( ptSegEn) ; - // Loop aperto - for ( int nOpenLine = int( cvOpenChain[nLastOpenLoopN].size()) - 1 ; nOpenLine >= 0 ; -- nOpenLine) { - Point3d ptOpenSt = cvOpenChain[nLastOpenLoopN][nOpenLine].ptEn ; - Point3d ptOpenEn = cvOpenChain[nLastOpenLoopN][nOpenLine].ptSt ; - cvSplitLoop1.emplace_back( ptOpenSt) ; - cvSplitLoop1.emplace_back( ptOpenEn) ; - if ( nOpenLine == int( cvOpenChain[nLastOpenLoopN].size()) - 1) { - Vector3d vtLast = ptSegSt - ptSegEn ; - vtLast.Normalize() ; - bFirstInside = vtLast * cvOpenChain[nLastOpenLoopN][nOpenLine].vtOuter < 0. ; - } - } - // Dall'inizio del loop aperto fino alla fine del segmento corrente - ptSegSt = ptOpenLoopStP ; - ptSegEn = cvBoundClosedLoopVec[nLoop][nLine+1] ; - cvSplitLoop1.emplace_back( ptSegSt) ; - cvSplitLoop1.emplace_back( ptSegEn) ; - bFirstSt = false ; - } - } - } - // Creo secondo loop - if ( bFirstSt) { - // Tratto segmento ove cadono gli estremi del loop aperto - cvSplitLoop2.emplace_back( ptOpenLoopStP) ; - cvSplitLoop2.emplace_back( ptOpenLoopEnP) ; - // Loop aperto - for ( int nOpenLine = int( cvOpenChain[nLastOpenLoopN].size()) - 1 ; nOpenLine >= 0 ; -- nOpenLine) { - Point3d ptOpenSt = cvOpenChain[nLastOpenLoopN][nOpenLine].ptEn ; - Point3d ptOpenEn = cvOpenChain[nLastOpenLoopN][nOpenLine].ptSt ; - cvSplitLoop2.emplace_back( ptOpenSt) ; - cvSplitLoop2.emplace_back( ptOpenEn) ; - } - } - else { - // Tratto segmento ove cadono gli estremi del loop aperto - cvSplitLoop2.emplace_back( ptOpenLoopEnP) ; - cvSplitLoop2.emplace_back( ptOpenLoopStP) ; - // Loop aperto - for ( int nOpenLine = 0 ; nOpenLine < int( cvOpenChain[nLastOpenLoopN].size()) ; ++ nOpenLine) { - Point3d ptOpenSt = cvOpenChain[nLastOpenLoopN][nOpenLine].ptSt ; - Point3d ptOpenEn = cvOpenChain[nLastOpenLoopN][nOpenLine].ptEn ; - cvSplitLoop2.emplace_back( ptOpenSt) ; - cvSplitLoop2.emplace_back( ptOpenEn) ; - } - } - // Aggiungo i nuovi loop nel vettore - int nCurSize = int(cvBoundClosedLoopVec.size()); - cvBoundClosedLoopVec.resize(nCurSize + 1); - vbInOut.resize(nCurSize + 1); - for (int nCL = nCurSize - 1; nCL > nLoop; --nCL) { - cvBoundClosedLoopVec[nCL + 1] = cvBoundClosedLoopVec[nCL]; - vbInOut[nCL + 1] = vbInOut[nCL]; - } - cvBoundClosedLoopVec[nLoop] = cvSplitLoop1 ; - cvBoundClosedLoopVec[nLoop + 1] = cvSplitLoop2; - vbInOut[nLoop] = bFirstInside; - vbInOut[nLoop + 1] = !bFirstInside; - ++ nLoop ; - } - // Estremi dell'aperto su due diversi segmenti del chiuso - else { - bool bFirstInside ; - bool bFirstSt ; - PNTVECTOR cvSplitLoop1, cvSplitLoop2 ; - for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size()) - 1 ; nLine += 2) { - // Segmento su cui non incide nessun estremo del loop aperto - if ( nLine != nCvSt && nLine != nCvEn) { - cvSplitLoop1.emplace_back( cvBoundClosedLoopVec[nLoop][nLine]) ; - cvSplitLoop1.emplace_back( cvBoundClosedLoopVec[nLoop][nLine + 1]) ; - } - // Segmento su cui incide l'estremo iniziale del loop aperto - else if ( nLine == nCvSt) { - // Dall'inizio del segmento corrente fino all'inizio del loop aperto - Point3d ptSegSt = cvBoundClosedLoopVec[nLoop][nLine] ; - Point3d ptSegEn = ptOpenLoopStP ; - cvSplitLoop1.emplace_back( ptSegSt) ; - cvSplitLoop1.emplace_back( ptSegEn) ; - // Loop aperto - for ( int nOpenLine = 0 ; nOpenLine < int( cvOpenChain[nLastOpenLoopN].size()) ; ++ nOpenLine) { - // Valuto se cvSplitLoop1 è interno o esterno - Point3d ptOpenSt = cvOpenChain[nLastOpenLoopN][nOpenLine].ptSt ; - Point3d ptOpenEn = cvOpenChain[nLastOpenLoopN][nOpenLine].ptEn ; - cvSplitLoop1.emplace_back( ptOpenSt) ; - cvSplitLoop1.emplace_back( ptOpenEn) ; - if ( nOpenLine == 0) { - Vector3d vtLast = ptSegSt - ptOpenLoopStP ; - if ( vtLast.Normalize()) { - bFirstInside = vtLast * cvOpenChain[nLastOpenLoopN][nOpenLine].vtOuter < 0. ; - } - else { - vtLast = cvBoundClosedLoopVec[nLoop][nLine + 1] - ptOpenLoopStP ; - vtLast.Normalize() ; - bFirstInside = vtLast * cvOpenChain[nLastOpenLoopN][nOpenLine].vtOuter > 0. ; - } - } - } - // Dalla fine del loop aperto fino alla fine del segmento ove arriva il loop aperto - nLine = nCvEn ; - ptSegEn = cvBoundClosedLoopVec[nLoop][nLine+1] ; - cvSplitLoop1.emplace_back( ptOpenLoopEnP) ; - cvSplitLoop1.emplace_back( ptSegEn) ; - bFirstSt = true ; - } - // Segmento su cui incide l'estremo finale del loop aperto - else { - // Dall'inizio del segmento corrente fino alla fine del loop aperto - Point3d ptSegSt = cvBoundClosedLoopVec[nLoop][nLine] ; - cvSplitLoop1.emplace_back( ptSegSt) ; - cvSplitLoop1.emplace_back( ptOpenLoopEnP) ; - // Loop aperto - for ( int nOpenLine = int( cvOpenChain[nLastOpenLoopN].size()) - 1 ; nOpenLine >= 0 ; -- nOpenLine) { - Point3d ptOpenSt = cvOpenChain[nLastOpenLoopN][nOpenLine].ptEn ; - Point3d ptOpenEn = cvOpenChain[nLastOpenLoopN][nOpenLine].ptSt ; - cvSplitLoop1.emplace_back( ptOpenSt) ; - cvSplitLoop1.emplace_back( ptOpenEn) ; - if ( nOpenLine == int( cvOpenChain[nLastOpenLoopN].size()) - 1) { - Vector3d vtLast = ptSegSt - ptOpenLoopEnP ; - if ( vtLast.Normalize()) { - bFirstInside = vtLast * cvOpenChain[nLastOpenLoopN][nOpenLine].vtOuter < 0. ; - } - else { - bFirstInside = vtLast * cvOpenChain[nLastOpenLoopN][nOpenLine].vtOuter > 0. ; - } - } - } - // Dalla fine del loop aperto fino alla fine del segmento ove arriva il loop aperto - nLine = nCvSt ; - Point3d ptSegEn = cvBoundClosedLoopVec[nLoop][nLine + 1] ; - cvSplitLoop1.emplace_back( ptOpenLoopStP) ; - cvSplitLoop1.emplace_back( ptSegEn) ; - bFirstSt = false; - } - } - // Creo secondo loop - if ( bFirstSt) { - // Dall'estremo inizinale del loop aperto alla fine del segmento su cui giace - Point3d ptSegSt, ptSegEn ; - ptSegSt = ptOpenLoopStP ; - ptSegEn = cvBoundClosedLoopVec[nLoop][nCvSt+1] ; - cvSplitLoop2.emplace_back( ptSegSt) ; - cvSplitLoop2.emplace_back( ptSegEn) ; - // Loop chiuso fino a segmento ove cade l'estremo finale dell'aperto - for (int nLine = (nCvSt + 2) % int(cvBoundClosedLoopVec[nLoop].size()); nLine != nCvEn; nLine += 2) { - ptSegSt = cvBoundClosedLoopVec[nLoop][nLine % int(cvBoundClosedLoopVec[nLoop].size())]; - ptSegEn = cvBoundClosedLoopVec[nLoop][(nLine + 1) % int(cvBoundClosedLoopVec[nLoop].size())]; - cvSplitLoop2.emplace_back(ptSegSt); - cvSplitLoop2.emplace_back(ptSegEn); - } - // Dall'inizio del segmento del loop chiuso ove cade l'estremo finale del loop aperto - // a quest'ultimo - ptSegSt = cvBoundClosedLoopVec[nLoop][nCvEn] ; - ptSegEn = ptOpenLoopEnP ; - cvSplitLoop2.emplace_back( ptSegSt) ; - cvSplitLoop2.emplace_back( ptSegEn) ; - // Loop aperto - for ( int nOpenLine = int( cvOpenChain[nLastOpenLoopN].size()) - 1 ; nOpenLine >= 0 ; -- nOpenLine) { - Point3d ptOpenSt = cvOpenChain[nLastOpenLoopN][nOpenLine].ptEn ; - Point3d ptOpenEn = cvOpenChain[nLastOpenLoopN][nOpenLine].ptSt ; - cvSplitLoop2.emplace_back( ptOpenSt) ; - cvSplitLoop2.emplace_back( ptOpenEn) ; - } - } - else { - // Dall'estremo finale del loop aperto alla fine del segmento su cui giace - Point3d ptSegSt = ptOpenLoopEnP ; - Point3d ptSegEn = cvBoundClosedLoopVec[nLoop][nCvEn+1] ; - cvSplitLoop2.emplace_back( ptSegSt) ; - cvSplitLoop2.emplace_back( ptSegEn) ; - // Loop chiuso fino al segmento ove cade l'inizio del loop aperto - for (int nLine = (nCvEn + 2) % int(cvBoundClosedLoopVec[nLoop].size()); nLine != nCvSt; nLine += 2) { - ptSegSt = cvBoundClosedLoopVec[nLoop][nLine % int(cvBoundClosedLoopVec[nLoop].size())]; - ptSegEn = cvBoundClosedLoopVec[nLoop][(nLine + 1) % int(cvBoundClosedLoopVec[nLoop].size())]; - cvSplitLoop2.emplace_back(ptSegSt); - cvSplitLoop2.emplace_back(ptSegEn); - } - // Dall'inizio del segmento del loop chiuso ove cade l'estremo iniziale del loop aperto - // a quest'ultimo - ptSegSt = cvBoundClosedLoopVec[nLoop][nCvSt] ; - ptSegEn = ptOpenLoopStP ; - cvSplitLoop2.emplace_back( ptSegSt) ; - cvSplitLoop2.emplace_back( ptSegEn) ; - // Loop aperto - for ( int nOpenLoop = 0 ; nOpenLoop < int( cvOpenChain[nLastOpenLoopN].size()) ; ++ nOpenLoop) { - Point3d ptOpenSt = cvOpenChain[nLastOpenLoopN][nOpenLoop].ptSt ; - Point3d ptOpenEn = cvOpenChain[nLastOpenLoopN][nOpenLoop].ptEn ; - cvSplitLoop2.emplace_back( ptOpenSt) ; - cvSplitLoop2.emplace_back( ptOpenEn) ; - } - } - // Aggiungo i nuovi loop nel vettore - int nCurSize = int( cvBoundClosedLoopVec.size()) ; - cvBoundClosedLoopVec.resize( nCurSize + 1) ; - vbInOut.resize( nCurSize + 1) ; - for ( int nCL = nCurSize - 1 ; nCL > nLoop ; -- nCL) { - cvBoundClosedLoopVec[nCL + 1] = cvBoundClosedLoopVec[nCL] ; - vbInOut[nCL + 1] = vbInOut[nCL] ; - } - cvBoundClosedLoopVec[nLoop] = cvSplitLoop1 ; - cvBoundClosedLoopVec[nLoop + 1] = cvSplitLoop2 ; - vbInOut[nLoop] = bFirstInside ; - vbInOut[nLoop + 1] = ! bFirstInside ; - ++ nLoop ; - } - } + bool bAdded1 = AddChainToChain( cvCounterChain, Loop1) ; + bool bAdded2 = AddChainToChain( cvOpenChain[nLastOpenLoopN], Loop2) ; + if ( ! ( bAdded1 && bAdded2)) + continue ; + // Aggiungo i nuovi loop nel vettore + int nCurSize = int( cvBoundClosedLoopVec.size()) ; + cvBoundClosedLoopVec.resize( nCurSize + 1) ; + vbInOut.resize( nCurSize + 1) ; + for ( int nCL = nCurSize - 1 ; nCL > nLoop ; -- nCL) { + cvBoundClosedLoopVec[nCL + 1] = cvBoundClosedLoopVec[nCL] ; + vbInOut[nCL + 1] = vbInOut[nCL] ; + } + int nLastPointLoop2 = int( Loop2.size()) - 1 ; + Vector3d vtTest = Loop1[1] - Loop1[0] ; + vtTest.Normalize() ; + bool bSecondInside = vtTest * cvOpenChain[nLastOpenLoopN][0].vtOuter < 0. ; + cvBoundClosedLoopVec[nLoop] = Loop1 ; + cvBoundClosedLoopVec[nLoop + 1] = Loop2 ; + vbInOut[nLoop] = bSecondInside ; + vbInOut[nLoop + 1] = ! bSecondInside ; + ++ nLoop ; } } // Degenere @@ -707,13 +603,13 @@ SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) Point3d ptOpenLoopStP = cvOpenChain[nLastOpenLoopN][0].ptSt ; Point3d ptOpenLoopEnP = cvOpenChain[nLastOpenLoopN][0].ptEn ; // Cerco se esistono dei tratti del loop chiuso corrente che sono - // toccati dagli estremi del loop aperto corrente + // toccati dagli estremi del loop aperto corrente int nCvFirst = -1 ; int nCvSecond = -1 ; - for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size()) - 1 && nCvSecond == -1 ; nLine += 2) { + for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size()) && nCvSecond == -1 ; ++ nLine) { // Estremi del segmento corrente del loop chiuso corrente Point3d ptSegSt = cvBoundClosedLoopVec[nLoop][nLine] ; - Point3d ptSegEn = cvBoundClosedLoopVec[nLoop][nLine + 1] ; + Point3d ptSegEn = cvBoundClosedLoopVec[nLoop][( nLine + 1) % int( cvBoundClosedLoopVec[nLoop].size())] ; // Vettore congiungente i su definiti punti Vector3d vtClosedLoopSeg = ptSegEn - ptSegSt ; vtClosedLoopSeg.Normalize() ; @@ -728,6 +624,7 @@ SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) nCvSecond = nLine ; } } + if ( nCvFirst != nCvSecond && nCvSecond != -1) { // li ordino in senso crescente if ( nCvFirst > nCvSecond) @@ -735,55 +632,68 @@ SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) // punto medio tra primo e secondo int nCount = 0 ; Point3d ptM12 ; - for ( int i = nCvFirst + 1 ; i < nCvSecond ; i += 2) { + for ( int i = nCvFirst + 1 ; i <= nCvSecond ; ++ i) { ptM12 += cvBoundClosedLoopVec[nLoop][i] ; ++ nCount ; } ptM12 /= nCount ; - // verifico se questo punto è dalla parte valida o no - bool bC12 = ( ( ptM12 - ptProva) * vtVecProva < 0) ; - // numero totale di punti - int nPntTot = int( cvBoundClosedLoopVec[nLoop].size()) ; + // Distanza quadrata media dei punti tra primo e secondo dal baricentro + double dVar12 = 0. ; + for ( int i = nCvFirst + 1 ; i <= nCvSecond ; ++ i) { + dVar12 += ( cvBoundClosedLoopVec[nLoop][i] - ptM12) * ( cvBoundClosedLoopVec[nLoop][i] - ptM12) ; + } + dVar12 /= nCount ; + // punto medio fra secondo e primo + nCount = 0 ; + Point3d ptM21 ; + for ( int i = nCvSecond + 1 ; i % int( cvBoundClosedLoopVec[nLoop].size()) ; ++ i) { + ptM21 += cvBoundClosedLoopVec[nLoop][i] ; + ++ nCount ; + } + for ( int i = 0 ; i <= nCvFirst ; ++ i) { + ptM21 += cvBoundClosedLoopVec[nLoop][i] ; + ++ nCount ; + } + ptM21 /= nCount ; + // Distanza quadrata media dei punti tra secondo e primo dal baricentro + double dVar21 = 0. ; + for ( int i = nCvSecond ; i < i % int( cvBoundClosedLoopVec[nLoop].size()) ; ++ i) { + dVar21 += ( cvBoundClosedLoopVec[nLoop][i] - ptM21) * ( cvBoundClosedLoopVec[nLoop][i] - ptM21) ; + ++ nCount ; + } + for ( int i = 0 ; i <= nCvFirst ; ++ i) { + dVar21 += ( cvBoundClosedLoopVec[nLoop][i] - ptM21) * ( cvBoundClosedLoopVec[nLoop][i] - ptM21) ; + ++ nCount ; + } + dVar21 /= nCount ; // elimino i punti dalla parte non valida - if ( bC12) { + if ( dVar12 > dVar21) { // assegno i nuovi valori cvBoundClosedLoopVec[nLoop][nCvFirst] = ptProva ; - cvBoundClosedLoopVec[nLoop][nCvSecond+1] = ptProva ; + cvBoundClosedLoopVec[nLoop][( nCvSecond + 1) % int(cvBoundClosedLoopVec[nLoop].size())] = ptProva ; + // numero totale di punti + int nPntTot = int( cvBoundClosedLoopVec[nLoop].size()); // elimino i punti superflui dopo - for ( int i = nPntTot - 1 ; i > nCvSecond+1 ; -- i) + for ( int i = nPntTot - 1 ; i > nCvSecond + 1 ; -- i) cvBoundClosedLoopVec[nLoop].pop_back() ; // elimino i punti superflui prima for ( int i = 0 ; i < nCvFirst ; ++ i) cvBoundClosedLoopVec[nLoop].erase( cvBoundClosedLoopVec[nLoop].begin()) ; + // verifico se questo punto è dalla parte valida o no + bool bC12 = ( ( ptM12 - ptProva) * vtVecProva < 0) ; + vbInOut[nLoop] = bC12 ; } else { // assegno i nuovi valori - cvBoundClosedLoopVec[nLoop][nCvFirst+1] = ptProva ; + cvBoundClosedLoopVec[nLoop][nCvFirst + 1] = ptProva ; cvBoundClosedLoopVec[nLoop][nCvSecond] = ptProva ; // elimino i punti superflui intermedi - for ( int i = nCvFirst+2 ; i < nCvSecond ; ++ i) - cvBoundClosedLoopVec[nLoop].erase( cvBoundClosedLoopVec[nLoop].begin()+nCvFirst+2) ; + for ( int i = nCvFirst + 2 ; i < nCvSecond ; ++ i) + cvBoundClosedLoopVec[nLoop].erase( cvBoundClosedLoopVec[nLoop].begin() + nCvFirst + 2) ; + // verifico se questo punto è dalla parte valida o no + bool bC21 = ( ( ptM21 - ptProva) * vtVecProva < 0) ; + vbInOut[nLoop] = bC21 ; } - vbInOut[nLoop] = true ; - - // PNTVECTOR cvSplitLoop; - // for (int nLine = 0; nLine < nCvFirst; nLine += 2) { - // cvSplitLoop.emplace_back(cvBoundClosedLoopVec[nLoop][nLine]); - // cvSplitLoop.emplace_back(cvBoundClosedLoopVec[nLoop][nLine + 1]); - // } - // cvSplitLoop.emplace_back(cvBoundClosedLoopVec[nLoop][nCvFirst]); - // cvSplitLoop.emplace_back(ptProva); - // cvSplitLoop.emplace_back(ptProva); - // cvSplitLoop.emplace_back(cvBoundClosedLoopVec[nLoop][nCvSecond + 1]); - // for (int nLine = nCvSecond + 2; nLine < int(cvBoundClosedLoopVec[nLoop].size()) - 1; nLine += 2) { - // cvSplitLoop.emplace_back(cvBoundClosedLoopVec[nLoop][nLine]); - // cvSplitLoop.emplace_back(cvBoundClosedLoopVec[nLoop][nLine + 1]); - // } - // cvBoundClosedLoopVec[nLoop] = cvSplitLoop; - // if ((cvBoundClosedLoopVec[nLoop][nCvFirst] - ptProva) * vtVecProva < 0.) - // vbInOut[nLoop] = true; - // else - // vbInOut[nLoop] = false; } } } @@ -795,11 +705,11 @@ SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) // Trasformo i loop compositi in loop polyline POLYLINEVECTOR vplPolyVec ; vplPolyVec.resize( cvBoundClosedLoopVec.size()) ; - for ( int nLoop = 0 ; nLoop < int( vplPolyVec.size()) ; ++ nLoop) { - for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size()) - 1 ; nLine += 2) { + for ( int nLoop = 0 ; nLoop < int( vplPolyVec.size()) ; ++ nLoop) { + for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size()) ; ++ nLine) { vplPolyVec[nLoop].AddUPoint( 0., cvBoundClosedLoopVec[nLoop][nLine]) ; - vplPolyVec[nLoop].AddUPoint( 0., cvBoundClosedLoopVec[nLoop][nLine + 1]) ; } + vplPolyVec[nLoop].AddUPoint( 0., cvBoundClosedLoopVec[nLoop][0]) ; if ( vbInOut[nLoop]) { // Eseguo triangolazione Triangulate CreateTriangulation ; @@ -813,14 +723,14 @@ SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) AddVertex( vPt[nNewTriaVertId[1]]), AddVertex( vPt[nNewTriaVertId[2]]) } ; AddTriangle( nNewId) ; - bModif = true ; + bModif = true ; } } } - } + } } else if ( nVertInside == 0) - RemoveTriangle( nT) ; + RemoveTriangle( nT) ; } // Se avvenuta modifica, aggiorno tutto @@ -829,3 +739,1268 @@ SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq) return true ; } + +//---------------------------------------------------------------------------- +bool +SurfTriMesh::DecomposeLoop( CHAINVECTOR& cvOpenChain, INTVECTOR& vnDegVec, vector& cvBoundClosedLoopVec, vector& vbInOut) +{ + // Divido il loop di partenza in sotto-loop + while ( cvOpenChain.size() > 0) { + int nLastOpenLoopN = int( cvOpenChain.size()) - 1 ; + if ( vnDegVec[nLastOpenLoopN] == 1) { + for ( int nLoop = 0 ; nLoop < int( cvBoundClosedLoopVec.size()) ; ++ nLoop) { + // Estremi del loop aperto + int nLastOpenLoopPoint = max( int( cvOpenChain[nLastOpenLoopN].size()) - 1, 0) ; + Point3d ptOpenLoopStP = cvOpenChain[nLastOpenLoopN][0].ptSt ; + Point3d ptOpenLoopEnP = cvOpenChain[nLastOpenLoopN][nLastOpenLoopPoint].ptEn ; + PNTVECTOR Loop1, Loop2 ; + bool bChangedStart = ChangeStart( ptOpenLoopStP, cvBoundClosedLoopVec[nLoop]) ; + bool bSplitted = SplitAtPoint( ptOpenLoopEnP, cvBoundClosedLoopVec[nLoop], Loop1, Loop2) ; + if ( ! ( bChangedStart && bSplitted)) + continue ; + Chain cvCounterChain ; + for ( int nPt = int( cvOpenChain[nLastOpenLoopN].size()) - 1 ; nPt >= 0 ; -- nPt) { + IntSegment CurSeg ; + CurSeg.ptSt = cvOpenChain[nLastOpenLoopN][nPt].ptEn ; + CurSeg.ptEn = cvOpenChain[nLastOpenLoopN][nPt].ptSt ; + CurSeg.vtOuter = - cvOpenChain[nLastOpenLoopN][nPt].vtOuter ; + CurSeg.bDegenerate = cvOpenChain[nLastOpenLoopN][nPt].bDegenerate ; + cvCounterChain.emplace_back( CurSeg) ; + } + bool bAdded1 = AddChainToChain( cvCounterChain, Loop1) ; + bool bAdded2 = AddChainToChain( cvOpenChain[nLastOpenLoopN], Loop2) ; + if ( ! ( bAdded1 && bAdded2)) + continue ; + // Aggiungo i nuovi loop nel vettore + int nCurSize = int( cvBoundClosedLoopVec.size()) ; + cvBoundClosedLoopVec.resize( nCurSize + 1) ; + vbInOut.resize( nCurSize + 1) ; + for ( int nCL = nCurSize - 1 ; nCL > nLoop ; -- nCL) { + cvBoundClosedLoopVec[nCL + 1] = cvBoundClosedLoopVec[nCL] ; + vbInOut[nCL + 1] = vbInOut[nCL] ; + } + int nLastPointLoop2 = int( Loop2.size()) - 1 ; + Vector3d vtTest = Loop1[1] - Loop1[0] ; + vtTest.Normalize() ; + bool bSecondInside = vtTest * cvOpenChain[nLastOpenLoopN][0].vtOuter < 0. ; + cvBoundClosedLoopVec[nLoop] = Loop1 ; + cvBoundClosedLoopVec[nLoop + 1] = Loop2 ; + vbInOut[nLoop] = bSecondInside ; + vbInOut[nLoop + 1] = ! bSecondInside ; + ++ nLoop ; + } + } + // Degenere + else { + Point3d ptProva = 0.5 * ( cvOpenChain[nLastOpenLoopN][0].ptSt + cvOpenChain[nLastOpenLoopN][0].ptEn) ; + Vector3d vtVecProva = cvOpenChain[nLastOpenLoopN][0].vtOuter ; + vtVecProva.Normalize( EPS_ZERO) ; + for ( int nLoop = 0 ; nLoop < int( cvBoundClosedLoopVec.size()) ; ++ nLoop) { + // Estremi del loop aperto + int nLastOpenLoopPoint = max( int( cvOpenChain[nLastOpenLoopN].size()) - 1, 0) ; + Point3d ptOpenLoopStP = cvOpenChain[nLastOpenLoopN][0].ptSt ; + Point3d ptOpenLoopEnP = cvOpenChain[nLastOpenLoopN][0].ptEn ; + // Cerco se esistono dei tratti del loop chiuso corrente che sono + // toccati dagli estremi del loop aperto corrente + int nCvFirst = - 1 ; + int nCvSecond = - 1 ; + for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size()) && nCvSecond == - 1 ; ++ nLine) { + // Estremi del segmento corrente del loop chiuso corrente + Point3d ptSegSt = cvBoundClosedLoopVec[nLoop][nLine] ; + Point3d ptSegEn = cvBoundClosedLoopVec[nLoop][( nLine + 1) % int(cvBoundClosedLoopVec[nLoop].size())] ; + // Vettore congiungente i su definiti punti + Vector3d vtClosedLoopSeg = ptSegEn - ptSegSt ; + vtClosedLoopSeg.Normalize() ; + // Vedo se gli estremi del loop aperto stanno su un segmento del chiuso + DistPointLine DistCalc( ptProva, ptSegSt, ptSegEn) ; + double dSqDist ; + DistCalc.GetSqDist( dSqDist) ; + if ( dSqDist < 2 * SQ_EPS_SMALL) { + if ( nCvFirst == - 1) + nCvFirst = nLine ; + else + nCvSecond = nLine ; + } + } + + if ( nCvFirst != nCvSecond && nCvSecond != - 1) { + // li ordino in senso crescente + if ( nCvFirst > nCvSecond) + swap( nCvFirst, nCvSecond) ; + // punto medio tra primo e secondo + int nCount = 0 ; + Point3d ptM12 ; + for ( int i = nCvFirst + 1 ; i <= nCvSecond ; ++ i) { + ptM12 += cvBoundClosedLoopVec[nLoop][i] ; + ++ nCount ; + } + ptM12 /= nCount ; + // Distanza quadrata media dei punti tra primo e secondo dal baricentro + double dVar12 = 0. ; + for ( int i = nCvFirst + 1 ; i <= nCvSecond ; ++ i) { + dVar12 += ( cvBoundClosedLoopVec[nLoop][i] - ptM12) * ( cvBoundClosedLoopVec[nLoop][i] - ptM12) ; + } + dVar12 /= nCount ; + // punto medio fra secondo e primo + nCount = 0 ; + Point3d ptM21 ; + for ( int i = nCvSecond + 1 ; i % int( cvBoundClosedLoopVec[nLoop].size()) ; ++ i) { + ptM21 += cvBoundClosedLoopVec[nLoop][i] ; + ++ nCount ; + } + for ( int i = 0 ; i <= nCvFirst ; ++ i) { + ptM21 += cvBoundClosedLoopVec[nLoop][i] ; + ++ nCount ; + } + ptM21 /= nCount ; + // Distanza quadrata media dei punti tra secondo e primo dal baricentro + double dVar21 = 0. ; + for ( int i = nCvSecond ; i < i % int( cvBoundClosedLoopVec[nLoop].size()) ; ++ i) { + dVar21 += ( cvBoundClosedLoopVec[nLoop][i] - ptM21) * ( cvBoundClosedLoopVec[nLoop][i] - ptM21) ; + ++ nCount ; + } + for ( int i = 0 ; i <= nCvFirst ; ++ i) { + dVar21 += ( cvBoundClosedLoopVec[nLoop][i] - ptM21) * ( cvBoundClosedLoopVec[nLoop][i] - ptM21) ; + ++ nCount ; + } + dVar21 /= nCount ; + // elimino i punti dalla parte non valida + if ( dVar12 > dVar21) { + // assegno i nuovi valori + cvBoundClosedLoopVec[nLoop][nCvFirst] = ptProva ; + cvBoundClosedLoopVec[nLoop][( nCvSecond + 1) % int( cvBoundClosedLoopVec[nLoop].size())] = ptProva ; + // numero totale di punti + int nPntTot = int( cvBoundClosedLoopVec[nLoop].size()) ; + // elimino i punti superflui dopo + for ( int i = nPntTot - 1 ; i > nCvSecond + 1 ; -- i) + cvBoundClosedLoopVec[nLoop].pop_back() ; + // elimino i punti superflui prima + for ( int i = 0 ; i < nCvFirst ; ++ i) + cvBoundClosedLoopVec[nLoop].erase( cvBoundClosedLoopVec[nLoop].begin()) ; + // verifico se questo punto è dalla parte valida o no + bool bC12 = ( ( ptM12 - ptProva) * vtVecProva < 0) ; + vbInOut[nLoop] = bC12 ; + } + else { + // assegno i nuovi valori + cvBoundClosedLoopVec[nLoop][nCvFirst + 1] = ptProva ; + cvBoundClosedLoopVec[nLoop][nCvSecond] = ptProva ; + // elimino i punti superflui intermedi + for ( int i = nCvFirst + 2 ; i < nCvSecond ; ++ i) + cvBoundClosedLoopVec[nLoop].erase( cvBoundClosedLoopVec[nLoop].begin() + nCvFirst + 2) ; + // verifico se questo punto è dalla parte valida o no + bool bC21 = ( ( ptM21 - ptProva) * vtVecProva < 0) ; + vbInOut[nLoop] = bC21 ; + } + } + } + } + vnDegVec.resize( nLastOpenLoopN) ; + cvOpenChain.resize( nLastOpenLoopN) ; + } + return true ; +} + +//---------------------------------------------------------------------------- +bool +SurfTriMesh::IntersectTriMeshTriangle( SurfTriMesh& Other) +{ + bool bModif = false ; + SurfTriMesh& SurfB = Other ; + // Le superfici devono essere valide + if ( m_nStatus != OK || ! SurfB.IsValid()) + return false ; + // Unordered map dei segmenti di intersezione + std::unordered_map LineMapA ; + std::unordered_map LineMapB ; + // Ciclo sui triangoli delle mesh + int nTriaNumA = GetTriangleSize() ; + int nTriaNumB = SurfB.GetTriangleSize() ; + // Setto il triangolo come né fuori né dentro + for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) + m_vTria[nTA].nTempPart = 0 ; + for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) + SurfB.m_vTria[nTB].nTempPart = 0 ; + for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) { + // Se il triangolo A non è valido, continuo + Triangle3d trTriaA ; + if ( ! ( GetTriangle( nTA, trTriaA) && trTriaA.Validate( true))) + continue ; + // Box del triangolo A + BBox3d b3dTriaA ; + trTriaA.GetLocalBBox( b3dTriaA) ; + bool bNewTriaA = true ; + for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) { + // Se il triangolo B non è valido, continuo + Triangle3d trTriaB ; + if ( ! ( SurfB.GetTriangle( nTB, trTriaB) && trTriaB.Validate( true))) + continue ; + // Box del triangolo B + BBox3d b3dTriaB ; + trTriaB.GetLocalBBox( b3dTriaB) ; + // Se i box non si sovrappongono, contiuno + if ( ! b3dTriaA.Overlaps( b3dTriaB)) + continue ; + // Interseco i triangoli + if ( abs( trTriaA.GetN() * trTriaB.GetN()) < 1 - EPS_ZERO) { + Point3d ptSegSt, ptSegEn ; + TRIA3DVECTOR vTria ; + int nIntType = IntersTriaTria( trTriaA, trTriaB, ptSegSt, ptSegEn, vTria) ; + if ( ! ( nFromSpecialToNormal[nIntType] == ITTT_NO || + nFromSpecialToNormal[nIntType] == ITTT_OVERLAPS || + nFromSpecialToNormal[nIntType] == ITTTS_VERT_VERT)) { + // Assegno i dati di intersezione + IntSegment CurInters ; + if ( nFromSpecialToNormal[nIntType] == ITTT_EDGE || nFromSpecialToNormal[nIntType] == ITTT_YES) { + CurInters.ptSt = ptSegSt ; + CurInters.ptEn = ptSegEn ; + CurInters.bDegenerate = false ; + } + else { + CurInters.ptSt = ptSegSt ; + CurInters.ptEn = ptSegSt ; + CurInters.bDegenerate = true ; + } + CurInters.vtOuter = trTriaB.GetN() ; + CurInters.vtOuter -= ( ( CurInters.vtOuter * trTriaA.GetN()) * trTriaA.GetN()) ; + CurInters.vtOuter.Normalize() ; + // Salvo intersezione per superficie A + bool bIntOnEndgeA = false ; + if ( ! ( nIntType == ITTTS_EDGE_EDGE_SEG || nIntType == ITTTS_EDGE_INT)) { + auto itA = LineMapA.find( nTA) ; + if ( itA != LineMapA.end()) { + itA->second.emplace_back( CurInters) ; + } + else { + Chain chTemp ; + chTemp.emplace_back( CurInters) ; + LineMapA.emplace( nTA, chTemp) ; + } + } + else + bIntOnEndgeA = true ; + + CurInters.vtOuter = trTriaA.GetN() ; + CurInters.vtOuter -= ( ( CurInters.vtOuter * trTriaB.GetN()) * trTriaB.GetN()) ; + CurInters.vtOuter.Normalize() ; + + // Salvo intersezione per superficie B + bool bIntOnEndgeB = false ; + if ( ! ( nIntType == ITTTS_EDGE_EDGE_SEG || nIntType == ITTTS_INT_EDGE)) { + auto itB = LineMapB.find( nTB) ; + if ( itB != LineMapB.end()) { + itB->second.emplace_back( CurInters) ; + } + else { + Chain chTemp ; + chTemp.emplace_back( CurInters) ; + LineMapB.emplace( nTB, chTemp) ; + } + } + else + bIntOnEndgeB = true ; + + if ( bIntOnEndgeA && ! bIntOnEndgeB) { + double dMaxDist = 0. ; + int nSegMaxDist = - 1 ; + for ( int nVA = 0 ; nVA < 3 ; ++ nVA) { + double dDist = abs( ( trTriaA.GetP( nVA) - trTriaB.GetP( 0)) * trTriaB.GetN()) ; + if ( dMaxDist < dDist) { + nSegMaxDist = nVA ; + dMaxDist = dDist ; + } + } + if ( nSegMaxDist >= 0) { + m_vTria[nTA].nTempPart = ( ( trTriaA.GetP( nSegMaxDist) - trTriaB.GetP( 0)) * trTriaB.GetN() < - EPS_SMALL ? 1 : - 1) ; + } + } + else if ( ! bIntOnEndgeA && bIntOnEndgeB) { + double dMaxDist = 0. ; + int nSegMaxDist = - 1 ; + for ( int nVB = 0 ; nVB < 3 ; ++ nVB) { + double dDist = abs( ( trTriaB.GetP( nVB) - trTriaA.GetP( 0)) * trTriaA.GetN()) ; + if ( dMaxDist < dDist) { + nSegMaxDist = nVB ; + dMaxDist = dDist ; + } + } + if ( nSegMaxDist >= 0) { + SurfB.m_vTria[nTB].nTempPart = ( ( trTriaB.GetP( nSegMaxDist) - trTriaA.GetP( 0)) * trTriaA.GetN() < - EPS_SMALL ? 1 : - 1) ; + } + } + } + else { + ; + } + } + } + } + + + // Ritriangolarizzo i triangoli della superficie A + for ( auto it = LineMapA.begin() ; it != LineMapA.end() ; ++ it) { + for ( int nS1 = 0 ; nS1 < int( it->second.size()) - 1 ; ++ nS1) { + for ( int nS2 = nS1 + 1 ; nS2 < int( it->second.size()) ; ++ nS2) { + if ( AreSamePointApprox( it->second[nS1].ptSt, it->second[nS2].ptEn) && + AreSamePointApprox( it->second[nS1].ptEn, it->second[nS2].ptSt) && + it->second[nS1].vtOuter * it->second[nS2].vtOuter < - EPS_SMALL) { + it->second.erase( it->second.begin() + nS2) ; + it->second.erase( it->second.begin() + nS1) ; + -- nS1 ; + -- nS2 ; + } + } + } + if ( int( it->second.size()) == 0) + continue ; + + // Recupero il triangolo + Triangle3d trTriaA ; + GetTriangle( it->first, trTriaA) ; + // Lo rimuovo dalla mesh + RemoveTriangle( it->first) ; + bModif = true ; + CHAINVECTOR vChain ; + // Creo i loop + ChainCurves LoopCreator ; + LoopCreator.Init( false, EPS_SMALL, int( it->second.size())) ; + // Carico le curve per concatenarle + for ( int nCv = 0 ; nCv < int( it->second.size()) ; ++ nCv) { + Point3d ptSt = it->second[nCv].ptSt ; + Point3d ptEn = it->second[nCv].ptEn ; + Vector3d vtDir = ptEn - ptSt ; + vtDir.Normalize() ; + LoopCreator.AddCurve( nCv + 1, ptSt, vtDir, ptEn, vtDir) ; + } + // Recupero i concatenamenti + INTVECTOR vIds ; + Point3d ptNearStart ; + while ( LoopCreator.GetChainFromNear( ptNearStart, false, vIds)) { + Chain chTemp ; + for ( auto i : vIds) { + // Aggiungo la linea alla curva composta. + chTemp.emplace_back( it->second[i - 1]) ; + } + vChain.emplace_back( chTemp) ; + } + // Lavoro su loop e catene per regolarizzarle + int nChainCnt = int( vChain.size()) ; + // unisco eventuali catene estreme che sono parte di una stessa catena + if ( nChainCnt > 1) { + if ( AreSamePointApprox( vChain[0].front().ptSt, vChain[nChainCnt - 1].back().ptEn)) { + vChain[0].insert( vChain[0].begin(), vChain[nChainCnt - 1].begin(), vChain[nChainCnt - 1].end()) ; + vChain.pop_back() ; + -- nChainCnt ; + } + else if ( AreSamePointApprox(vChain[0].back().ptEn, vChain[nChainCnt - 1].front().ptSt)) { + vChain[0].insert(vChain[0].end(), vChain[nChainCnt - 1].begin(), vChain[nChainCnt - 1].end()) ; + vChain.pop_back() ; + -- nChainCnt ; + } + } + // semplifico catene formate da punti degeneri + for ( int nCh = 0 ; nCh < nChainCnt ; ++ nCh) { + if ( vChain[nCh].size() == 2 && ( vChain[nCh][0].bDegenerate || vChain[nCh][1].bDegenerate)) { + vChain[nCh][0].ptEn = vChain[nCh][1].ptEn ; + vChain[nCh][0].vtOuter = ( vChain[nCh][0].bDegenerate ? vChain[nCh][1].vtOuter : vChain[nCh][0].vtOuter) ; + vChain[nCh][0].bDegenerate = AreSamePointApprox( vChain[nCh][0].ptSt, vChain[nCh][0].ptEn) ; + vChain[nCh].resize( 1) ; + } + } + // Elimino la seconda copia di catene doppie + for ( int nI = 0 ; nI < nChainCnt ; ++ nI) { + for ( int nJ = nI + 1 ; nJ < nChainCnt ; ++ nJ) { + if ( vChain[nI].size() == vChain[nJ].size()) { + bool bSame = true ; + for ( int nK = 0 ; nK < int( vChain[nI].size()) ; ++ nK) { + if ( ! AreSamePointApprox( vChain[nI][nK].ptSt, vChain[nJ][nK].ptSt) || + ! AreSamePointApprox( vChain[nI][nK].ptEn, vChain[nJ][nK].ptEn)) { + bSame = false ; + break ; + } + } + if ( bSame) { + vChain.erase( vChain.begin() + nJ) ; + -- nChainCnt ; + -- nJ ; + } + } + } + } + + // Fra le catene trovate separo le aperte dalle chiuse + int nDegenerateChainNum = 0 ; + INTVECTOR vnDegVec ; + CHAINVECTOR cvClosedChain ; + CHAINVECTOR cvOpenChain ; + for ( int nL = 0 ; nL < int( vChain.size()) ; ++ nL) { + bool bChainDegenerate = false ; + if ( vChain[nL].size() == 1 && AreSamePointApprox( vChain[nL][0].ptSt, vChain[nL][0].ptEn)) { + bChainDegenerate = true ; + } + + if ( bChainDegenerate) + ++ nDegenerateChainNum ; + int nCurLoopLast = max( int( vChain[nL].size()) - 1, 0) ; + if ( ( ! bChainDegenerate) && AreSamePointApprox( vChain[nL][0].ptSt, vChain[nL][nCurLoopLast].ptEn)) + cvClosedChain.emplace_back( vChain[nL]) ; + else { + cvOpenChain.emplace_back( vChain[nL]) ; + if ( bChainDegenerate) + vnDegVec.emplace_back( 0) ; + else + vnDegVec.emplace_back( 1) ; + } + } + + // Creo il loop chiuso padre di tutti, il perimetro del triangolo. + // Questo viene diviso in sotto-loop chiusi mediante quelli aperti. + // I loop chiusi trovati precedentemente sono interni a uno dei sotto-loop + // chiusi di cui è formato il perimetro. + PNTVECTOR cvFirstLoop ; + cvFirstLoop.emplace_back( trTriaA.GetP( 0)) ; + cvFirstLoop.emplace_back( trTriaA.GetP( 1)) ; + cvFirstLoop.emplace_back( trTriaA.GetP( 2)) ; + + vector cvBoundClosedLoopVec ; + cvBoundClosedLoopVec.emplace_back( cvFirstLoop) ; + vector vbInOut ; + vbInOut.push_back( true) ; + // Divido il loop usando le catene + DecomposeLoop( cvOpenChain, vnDegVec, cvBoundClosedLoopVec, vbInOut) ; + + // Rimuovo il triangolo corrente + RemoveTriangle( it->first) ; + // Trasformo i loop compositi in loop polyline + POLYLINEVECTOR vplPolyVec ; + vplPolyVec.resize( cvBoundClosedLoopVec.size()) ; + for ( int nLoop = 0 ; nLoop < int( vplPolyVec.size()) ; ++ nLoop) { + for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size()) ; ++ nLine) { + vplPolyVec[nLoop].AddUPoint( 0., cvBoundClosedLoopVec[nLoop][nLine]) ; + } + vplPolyVec[nLoop].AddUPoint( 0., cvBoundClosedLoopVec[nLoop][0]) ; + + // Assegno ai loop trovati i rispettivi interni + // Assumo che i loop interni a uno dei loop creati fino ad'ora siano tutti sullo stesso livello. + // Il caso generale si risolve con una struttura ad albero in cui il nodi corrispondente a un + // loop è figlio del nodo corrispondente al loop che lo contiene. + vector vInnerLoop ; + for ( int nCLI = 0 ; nCLI < int( cvClosedChain.size()) ; ++ nCLI) { + Point3d ptLoopStart = cvClosedChain[nCLI][0].ptSt ; + double dMinDist = DBL_MAX ; + Point3d ptMinDist ; + bool bPointOnSt = false ; + bool bPointOnEn = false ; + int nSegNum = 0 ; + int nSegMin ; + Point3d ptS, ptE ; + bool bContinueS = vplPolyVec[nLoop].GetFirstPoint( ptS) ; + bool bContinueE = vplPolyVec[nLoop].GetNextPoint( ptE) ; + while ( bContinueS && bContinueE) { + ++ nSegNum ; + DistPointLine DistCalculator( ptLoopStart, ptS, ptE) ; + double dDist ; + DistCalculator.GetDist( dDist) ; + if ( dDist < dMinDist) { + DistCalculator.GetMinDistPoint( ptMinDist) ; + bPointOnSt = AreSamePointExact( ptMinDist, ptS) ; + bPointOnEn = AreSamePointExact( ptMinDist, ptE) ; + dMinDist = dDist ; + nSegMin = nSegNum ; + } + ptS = ptE ; + bContinueS = bContinueE ; + bContinueE = vplPolyVec[nLoop].GetNextPoint( ptE) ; + } + if ( ! ( bPointOnSt || bPointOnEn)) { + vplPolyVec[nLoop].GetFirstPoint( ptS) ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + for ( int nSeg = 1 ; nSeg < nSegMin ; ++ nSeg) { + ptS = ptE ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + } + Vector3d vtTan = ptE - ptS ; + vtTan.Normalize() ; + Vector3d vtOut = vtTan ^ trTriaA.GetN() ; + Point3d ptMinDist ; + DistPointLine DistCalculator( ptLoopStart, ptS, ptE) ; + DistCalculator.GetMinDistPoint( ptMinDist) ; + double dMinDistDot = ( ptLoopStart - ptMinDist) * vtOut ; + if ( dMinDistDot < 0.) + vInnerLoop.emplace_back( nCLI) ; + } + else if ( bPointOnSt) { + Point3d ptPrevS, ptPrevE ; + if ( nSegMin == 1) { + vplPolyVec[nLoop].GetFirstPoint( ptS) ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + vplPolyVec[nLoop].GetLastPoint( ptPrevE) ; + vplPolyVec[nLoop].GetPrevPoint( ptPrevS) ; + } + else { + -- nSegMin ; + vplPolyVec[nLoop].GetFirstPoint( ptPrevS) ; + vplPolyVec[nLoop].GetNextPoint( ptPrevE) ; + for ( int nSeg = 1 ; nSeg < nSegMin ; ++ nSeg) { + ptPrevS = ptPrevE ; + vplPolyVec[nLoop].GetNextPoint( ptPrevE) ; + } + ptS = ptPrevE ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + } + Vector3d vtTan = ptE - ptS ; + vtTan.Normalize() ; + Vector3d vtTanPrev = ptPrevE - ptPrevS ; + vtTanPrev.Normalize() ; + Vector3d vtBisector = 0.5 * ( vtTan + vtTanPrev) ^ trTriaA.GetN() ; + vtBisector.Normalize() ; + double dMinDistDot = ( ptLoopStart - ptMinDist) * vtBisector ; + if ( dMinDistDot < 0.) + vInnerLoop.emplace_back(nCLI) ; + } + else if ( bPointOnEn) { + Point3d ptLast ; + vplPolyVec[nLoop].GetLastPoint( ptLast) ; + vplPolyVec[nLoop].GetFirstPoint( ptS) ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + for ( int nSeg = 1 ; nSeg < nSegMin ; ++ nSeg) { + ptS = ptE ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + } + Point3d ptNextS, ptNextE ; + if ( AreSamePointExact( ptE, ptLast)) { + vplPolyVec[nLoop].GetFirstPoint( ptNextS) ; + vplPolyVec[nLoop].GetNextPoint( ptNextE) ; + } + else { + vplPolyVec[nLoop].GetNextPoint( ptNextS) ; + vplPolyVec[nLoop].GetNextPoint( ptNextE) ; + } + Vector3d vtTan = ptE - ptS ; + vtTan.Normalize() ; + Vector3d vtTanNext = ptNextE - ptNextS ; + vtTanNext.Normalize() ; + Vector3d vtBisector = 0.5 * ( vtTan + vtTanNext) ^ trTriaA.GetN() ; + vtBisector.Normalize() ; + double dMinDistDot = ( ptLoopStart - ptMinDist) * vtBisector ; + if ( dMinDistDot < 0.) + vInnerLoop.emplace_back( nCLI) ; + } + } + + if ( vInnerLoop.size() == 0) { + // Eseguo triangolazione + Triangulate CreateTriangulation ; + PNTVECTOR vPt ; + INTVECTOR vTr ; + if ( Triangulate().Make( vplPolyVec[nLoop], vPt, vTr)) { + // Inserisco i nuovi triangoli + for ( int n = 0 ; n < int( vTr.size()) - 2 ; n += 3) { + int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] } ; + int nNewId[3] = { AddVertex( vPt[nNewTriaVertId[0]]), + AddVertex( vPt[nNewTriaVertId[1]]), + AddVertex( vPt[nNewTriaVertId[2]]) } ; + int nNewTriaNum = AddTriangle( nNewId) ; + if ( vbInOut[nLoop]) + m_vTria[nNewTriaNum].nTempPart = 1 ; + else + m_vTria[nNewTriaNum].nTempPart = - 1 ; + bModif = true ; + } + } + } + else { + + POLYLINEVECTOR vPolygons ; + vPolygons.emplace_back( vplPolyVec[nLoop]) ; + + for ( int nL = 0 ; nL < int( vInnerLoop.size()) ; ++ nL) { + PolyLine CurLoop ; + for ( int nV = 0 ; nV < int( cvClosedChain[vInnerLoop[nL]].size()) ; ++ nV) { + CurLoop.AddUPoint( 0., cvClosedChain[vInnerLoop[nL]][nV].ptSt) ; + } + CurLoop.AddUPoint( 0., cvClosedChain[vInnerLoop[nL]][0].ptSt) ; + vPolygons.emplace_back( CurLoop) ; + } + + Polygon3d pgPol ; + pgPol.FromPolyLine( vPolygons[1]) ; + + if ( trTriaA.GetN() * pgPol.GetVersN() > 0.) { + for (int nL = 1; nL < int(vPolygons.size()); ++nL) { + vPolygons[nL].Invert(); + } + + PNTVECTOR vPt; + INTVECTOR vTr; + if (Triangulate().Make(vPolygons, vPt, vTr)) { + // Inserisco i nuovi triangoli + for (int n = 0; n < int(vTr.size()) - 2; n += 3) { + int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] }; + int nNewId[3] = { AddVertex(vPt[nNewTriaVertId[0]]), + AddVertex(vPt[nNewTriaVertId[1]]), + AddVertex(vPt[nNewTriaVertId[2]]) }; + int nNewTriaNum = AddTriangle(nNewId); + m_vTria[nNewTriaNum].nTempPart = -1; + bModif = true; + } + } + + for (int nL = 1; nL < int(vPolygons.size()); ++nL) { + vPolygons[nL].Invert(); + if (Triangulate().Make(vPolygons[nL], vPt, vTr)) { + // Inserisco i nuovi triangoli + for (int n = 0; n < int(vTr.size()) - 2; n += 3) { + int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] }; + int nNewId[3] = { AddVertex(vPt[nNewTriaVertId[0]]), + AddVertex(vPt[nNewTriaVertId[1]]), + AddVertex(vPt[nNewTriaVertId[2]]) }; + int nNewTriaNum = AddTriangle(nNewId); + m_vTria[nNewTriaNum].nTempPart = 1; + bModif = true; + } + } + } + } + else { + PNTVECTOR vPt; + INTVECTOR vTr; + if (Triangulate().Make(vPolygons, vPt, vTr)) { + // Inserisco i nuovi triangoli + for (int n = 0; n < int(vTr.size()) - 2; n += 3) { + int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] }; + int nNewId[3] = { AddVertex(vPt[nNewTriaVertId[0]]), + AddVertex(vPt[nNewTriaVertId[1]]), + AddVertex(vPt[nNewTriaVertId[2]]) }; + int nNewTriaNum = AddTriangle(nNewId); + m_vTria[nNewTriaNum].nTempPart = 1; + bModif = true; + } + } + + for (int nL = 1; nL < int(vPolygons.size()); ++nL) { + vPolygons[nL].Invert(); + if (Triangulate().Make(vPolygons[nL], vPt, vTr)) { + // Inserisco i nuovi triangoli + for (int n = 0; n < int(vTr.size()) - 2; n += 3) { + int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] }; + int nNewId[3] = { AddVertex(vPt[nNewTriaVertId[0]]), + AddVertex(vPt[nNewTriaVertId[1]]), + AddVertex(vPt[nNewTriaVertId[2]]) }; + int nNewTriaNum = AddTriangle(nNewId); + m_vTria[nNewTriaNum].nTempPart = -1; + bModif = true; + } + } + } + } + } + vInnerLoop.resize( 0) ; + } + } + + // Ritriangolarizzo i triangoli della superficie B + for ( auto it = LineMapB.begin() ; it != LineMapB.end() ; ++ it) { + for ( int nS1 = 0 ; nS1 < int( it->second.size()) - 1 ; ++ nS1) { + for ( int nS2 = nS1 + 1 ; nS2 < int( it->second.size()) ; ++ nS2) { + if ( AreSamePointApprox( it->second[nS1].ptSt, it->second[nS2].ptEn) && + AreSamePointApprox( it->second[nS1].ptEn, it->second[nS2].ptSt) && + it->second[nS1].vtOuter * it->second[nS2].vtOuter < - EPS_SMALL) { + it->second.erase( it->second.begin() + nS2) ; + it->second.erase( it->second.begin() + nS1) ; + -- nS1 ; + -- nS2 ; + } + } + } + if ( int( it->second.size()) == 0) + continue ; + + // Recupero il triangolo + Triangle3d trTriaB ; + SurfB.GetTriangle( it->first, trTriaB) ; + // Lo rimuovo dalla mesh + SurfB.RemoveTriangle( it->first) ; + bModif = true ; + CHAINVECTOR vChain ; + // Creo i loop + ChainCurves LoopCreator ; + LoopCreator.Init( false, EPS_SMALL, int( it->second.size())) ; + // Carico le curve per concatenarle + for ( int nCv = 0 ; nCv < int( it->second.size()) ; ++ nCv) { + Point3d ptSt = it->second[nCv].ptSt ; + Point3d ptEn = it->second[nCv].ptEn ; + Vector3d vtDir = ptEn - ptSt ; + vtDir.Normalize() ; + LoopCreator.AddCurve( nCv + 1, ptSt, vtDir, ptEn, vtDir) ; + } + // Recupero i concatenamenti + INTVECTOR vIds ; + Point3d ptNearStart ; + while ( LoopCreator.GetChainFromNear( ptNearStart, false, vIds)) { + Chain chTemp ; + for ( auto i : vIds) { + // Aggiungo la linea alla curva composta. + chTemp.emplace_back( it->second[i - 1]) ; + } + vChain.emplace_back( chTemp) ; + } + // Lavoro su loop e catene per regolarizzarle + int nChainCnt = int( vChain.size()) ; + // unisco eventuali catene estreme che sono parte di una stessa catena + if ( nChainCnt > 1) { + if ( AreSamePointApprox( vChain[0].front().ptSt, vChain[nChainCnt - 1].back().ptEn)) { + vChain[0].insert( vChain[0].begin(), vChain[nChainCnt - 1].begin(), vChain[nChainCnt - 1].end()) ; + vChain.pop_back() ; + -- nChainCnt ; + } + else if ( AreSamePointApprox( vChain[0].back().ptEn, vChain[nChainCnt - 1].front().ptSt)) { + vChain[0].insert( vChain[0].end(), vChain[nChainCnt - 1].begin(), vChain[nChainCnt - 1].end()) ; + vChain.pop_back() ; + -- nChainCnt ; + } + } + // semplifico catene formate da punti degeneri + for ( int nCh = 0 ; nCh < nChainCnt ; ++ nCh) { + if ( vChain[nCh].size() == 2 && ( vChain[nCh][0].bDegenerate || vChain[nCh][1].bDegenerate)) { + vChain[nCh][0].ptEn = vChain[nCh][1].ptEn ; + vChain[nCh][0].vtOuter = ( vChain[nCh][0].bDegenerate ? vChain[nCh][1].vtOuter : vChain[nCh][0].vtOuter) ; + vChain[nCh][0].bDegenerate = AreSamePointApprox( vChain[nCh][0].ptSt, vChain[nCh][0].ptEn) ; + vChain[nCh].resize( 1) ; + } + } + // Elimino la seconda copia di catene doppie + for ( int nI = 0 ; nI < nChainCnt ; ++ nI) { + for ( int nJ = nI + 1 ; nJ < nChainCnt ; ++ nJ) { + if ( vChain[nI].size() == vChain[nJ].size()) { + bool bSame = true ; + for ( int nK = 0 ; nK < int( vChain[nI].size()) ; ++ nK) { + if ( ! AreSamePointApprox( vChain[nI][nK].ptSt, vChain[nJ][nK].ptSt) || + ! AreSamePointApprox( vChain[nI][nK].ptEn, vChain[nJ][nK].ptEn)) { + bSame = false ; + break ; + } + } + if ( bSame) { + vChain.erase( vChain.begin() + nJ) ; + -- nChainCnt ; + -- nJ ; + } + } + } + } + + // Fra le catene trovate separo le aperte dalle chiuse + int nDegenerateChainNum = 0 ; + INTVECTOR vnDegVec ; + CHAINVECTOR cvClosedChain ; + CHAINVECTOR cvOpenChain ; + for ( int nL = 0 ; nL < int( vChain.size()) ; ++ nL) { + bool bChainDegenerate = false ; + if ( vChain[nL].size() == 1 && AreSamePointApprox( vChain[nL][0].ptSt, vChain[nL][0].ptEn)) { + bChainDegenerate = true ; + } + + if ( bChainDegenerate) + ++ nDegenerateChainNum ; + int nCurLoopLast = max( int( vChain[nL].size()) - 1, 0) ; + if ( ( ! bChainDegenerate) && AreSamePointApprox( vChain[nL][0].ptSt, vChain[nL][nCurLoopLast].ptEn)) + cvClosedChain.emplace_back( vChain[nL]) ; + else { + cvOpenChain.emplace_back( vChain[nL]) ; + if ( bChainDegenerate) + vnDegVec.emplace_back( 0) ; + else + vnDegVec.emplace_back( 1) ; + } + } + + // Creo il loop chiuso padre di tutti, il perimetro del triangolo. + // Questo viene diviso in sotto-loop chiusi mediante quelli aperti. + // I loop chiusi trovati precedentemente sono interni a uno dei sotto-loop + // chiusi di cui è formato il perimetro. + PNTVECTOR cvFirstLoop ; + cvFirstLoop.emplace_back( trTriaB.GetP( 0)) ; + cvFirstLoop.emplace_back( trTriaB.GetP( 1)) ; + cvFirstLoop.emplace_back( trTriaB.GetP( 2)) ; + + vector cvBoundClosedLoopVec ; + cvBoundClosedLoopVec.emplace_back( cvFirstLoop) ; + vector vbInOut ; + vbInOut.push_back( true) ; + // Divido il loop usando le catene + DecomposeLoop( cvOpenChain, vnDegVec, cvBoundClosedLoopVec, vbInOut) ; + // Rimuovo il triangolo corrente + SurfB.RemoveTriangle( it->first) ; + // Trasformo i loop compositi in loop polyline + POLYLINEVECTOR vplPolyVec ; + vplPolyVec.resize( cvBoundClosedLoopVec.size()) ; + for ( int nLoop = 0 ; nLoop < int( vplPolyVec.size()) ; ++ nLoop) { + for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size()) ; ++ nLine) { + vplPolyVec[nLoop].AddUPoint( 0., cvBoundClosedLoopVec[nLoop][nLine]) ; + } + vplPolyVec[nLoop].AddUPoint( 0., cvBoundClosedLoopVec[nLoop][0]) ; + + // Assegno ai loop trovati i rispettivi interni + // Assumo che i loop interni a uno dei loop creati fino ad'ora siano tutti sullo stesso livello. + // Il caso generale si risolve con una struttura ad albero in cui il nodi corrispondente a un + // loop è figlio del nodo corrispondente al loop che lo contiene. + vector vInnerLoop ; + for ( int nCLI = 0 ; nCLI < int( cvClosedChain.size()) ; ++ nCLI) { + Point3d ptLoopStart = cvClosedChain[nCLI][0].ptSt ; + double dMinDist = DBL_MAX ; + Point3d ptMinDist ; + bool bPointOnSt = false ; + bool bPointOnEn = false ; + int nSegNum = 0 ; + int nSegMin ; + Point3d ptS, ptE ; + bool bContinueS = vplPolyVec[nLoop].GetFirstPoint( ptS) ; + bool bContinueE = vplPolyVec[nLoop].GetNextPoint( ptE) ; + while ( bContinueS && bContinueE) { + ++ nSegNum ; + DistPointLine DistCalculator( ptLoopStart, ptS, ptE) ; + double dDist ; + DistCalculator.GetDist( dDist) ; + if ( dDist < dMinDist) { + DistCalculator.GetMinDistPoint( ptMinDist) ; + bPointOnSt = AreSamePointExact( ptMinDist, ptS) ; + bPointOnEn = AreSamePointExact( ptMinDist, ptE) ; + dMinDist = dDist ; + nSegMin = nSegNum ; + } + ptS = ptE ; + bContinueS = bContinueE ; + bContinueE = vplPolyVec[nLoop].GetNextPoint( ptE) ; + } + if ( ! ( bPointOnSt || bPointOnEn)) { + vplPolyVec[nLoop].GetFirstPoint( ptS) ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + for ( int nSeg = 1 ; nSeg < nSegMin ; ++ nSeg) { + ptS = ptE ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + } + Vector3d vtTan = ptE - ptS ; + vtTan.Normalize() ; + Vector3d vtOut = vtTan ^ trTriaB.GetN() ; + Point3d ptMinDist ; + DistPointLine DistCalculator( ptLoopStart, ptS, ptE) ; + DistCalculator.GetMinDistPoint( ptMinDist) ; + double dMinDistDot = ( ptLoopStart - ptMinDist) * vtOut ; + if ( dMinDistDot < 0.) + vInnerLoop.emplace_back( nCLI) ; + } + else if ( bPointOnSt) { + Point3d ptPrevS, ptPrevE ; + if ( nSegMin == 1) { + vplPolyVec[nLoop].GetFirstPoint( ptS) ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + vplPolyVec[nLoop].GetLastPoint( ptPrevE) ; + vplPolyVec[nLoop].GetPrevPoint( ptPrevS) ; + } + else { + -- nSegMin ; + vplPolyVec[nLoop].GetFirstPoint( ptPrevS) ; + vplPolyVec[nLoop].GetNextPoint( ptPrevE) ; + for ( int nSeg = 1 ; nSeg < nSegMin ; ++ nSeg) { + ptPrevS = ptPrevE ; + vplPolyVec[nLoop].GetNextPoint( ptPrevE) ; + } + ptS = ptPrevE ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + } + Vector3d vtTan = ptE - ptS ; + vtTan.Normalize() ; + Vector3d vtTanPrev = ptPrevE - ptPrevS ; + vtTanPrev.Normalize() ; + Vector3d vtBisector = 0.5 * ( vtTan + vtTanPrev) ^ trTriaB.GetN() ; + vtBisector.Normalize() ; + double dMinDistDot = ( ptLoopStart - ptMinDist) * vtBisector ; + if ( dMinDistDot < 0.) + vInnerLoop.emplace_back( nCLI) ; + } + else if ( bPointOnEn) { + Point3d ptLast ; + vplPolyVec[nLoop].GetLastPoint( ptLast) ; + vplPolyVec[nLoop].GetFirstPoint( ptS) ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + for ( int nSeg = 1 ; nSeg < nSegMin ; ++ nSeg) { + ptS = ptE ; + vplPolyVec[nLoop].GetNextPoint( ptE) ; + } + Point3d ptNextS, ptNextE ; + if ( AreSamePointExact( ptE, ptLast)) { + vplPolyVec[nLoop].GetFirstPoint( ptNextS) ; + vplPolyVec[nLoop].GetNextPoint( ptNextE) ; + } + else { + vplPolyVec[nLoop].GetNextPoint( ptNextS) ; + vplPolyVec[nLoop].GetNextPoint( ptNextE) ; + } + Vector3d vtTan = ptE - ptS ; + vtTan.Normalize() ; + Vector3d vtTanNext = ptNextE - ptNextS ; + vtTanNext.Normalize() ; + Vector3d vtBisector = 0.5 * ( vtTan + vtTanNext) ^ trTriaB.GetN() ; + vtBisector.Normalize() ; + double dMinDistDot = ( ptLoopStart - ptMinDist) * vtBisector ; + if ( dMinDistDot < 0.) + vInnerLoop.emplace_back( nCLI) ; + } + } + + if ( vInnerLoop.size() == 0) { + // Eseguo triangolazione + Triangulate CreateTriangulation ; + PNTVECTOR vPt ; + INTVECTOR vTr ; + if ( Triangulate().Make( vplPolyVec[nLoop], vPt, vTr)) { + // Inserisco i nuovi triangoli + for ( int n = 0 ; n < int( vTr.size()) - 2 ; n += 3) { + int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] } ; + int nNewId[3] = { SurfB.AddVertex(vPt[nNewTriaVertId[0]]), + SurfB.AddVertex(vPt[nNewTriaVertId[1]]), + SurfB.AddVertex(vPt[nNewTriaVertId[2]]) } ; + int nNewTriaNum = SurfB.AddTriangle( nNewId) ; + if ( vbInOut[nLoop]) + SurfB.m_vTria[nNewTriaNum].nTempPart = 1 ; + else + SurfB.m_vTria[nNewTriaNum].nTempPart = - 1 ; + bModif = true ; + } + } + } + else { + POLYLINEVECTOR vPolygons ; + vPolygons.emplace_back( vplPolyVec[nLoop]) ; + if ( vbInOut[nLoop]) { + for ( int nL = 0 ; nL < int( vInnerLoop.size()) ; ++ nL) { + PolyLine CurLoop ; + for ( int nV = 0 ; nV < int( cvClosedChain[vInnerLoop[nL]].size()) ; ++ nV) { + CurLoop.AddUPoint( 0., cvClosedChain[vInnerLoop[nL]][nV].ptSt) ; + } + CurLoop.AddUPoint( 0., cvClosedChain[vInnerLoop[nL]][0].ptSt) ; + vPolygons.emplace_back( CurLoop) ; + } + + PNTVECTOR vPt ; + INTVECTOR vTr ; + if ( Triangulate().Make( vPolygons, vPt, vTr)) { + // Inserisco i nuovi triangoli + for ( int n = 0 ; n < int( vTr.size()) - 2 ; n += 3) { + int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] } ; + int nNewId[3] = { SurfB.AddVertex(vPt[nNewTriaVertId[0]]), + SurfB.AddVertex(vPt[nNewTriaVertId[1]]), + SurfB.AddVertex(vPt[nNewTriaVertId[2]]) } ; + int nNewTriaNum = SurfB.AddTriangle( nNewId) ; + SurfB.m_vTria[nNewTriaNum].nTempPart = 1 ; + bModif = true ; + } + } + + for ( int nL = 1 ; nL < int( vPolygons.size()) ; ++ nL) { + vPolygons[nL].Invert() ; + if ( Triangulate().Make(vPolygons[nL], vPt, vTr)) { + // Inserisco i nuovi triangoli + for ( int n = 0 ; n < int( vTr.size()) - 2 ; n += 3) { + int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] } ; + int nNewId[3] = { SurfB.AddVertex(vPt[nNewTriaVertId[0]]), + SurfB.AddVertex(vPt[nNewTriaVertId[1]]), + SurfB.AddVertex(vPt[nNewTriaVertId[2]]) } ; + int nNewTriaNum = SurfB.AddTriangle( nNewId) ; + SurfB.m_vTria[nNewTriaNum].nTempPart = - 1 ; + bModif = true ; + } + } + } + + } + else { + for ( int nL = 0 ; nL < int( vInnerLoop.size()); ++ nL) { + PolyLine CurLoop ; + for ( int nV = 0 ; nV < int( cvClosedChain[vInnerLoop[nL]].size()) ; ++ nV) { + CurLoop.AddUPoint( 0., cvClosedChain[vInnerLoop[nL]][nV].ptSt) ; + } + CurLoop.AddUPoint( 0., cvClosedChain[vInnerLoop[nL]][0].ptSt) ; + CurLoop.Invert() ; + vPolygons.emplace_back( CurLoop) ; + } + + PNTVECTOR vPt ; + INTVECTOR vTr ; + + if ( Triangulate().Make( vPolygons, vPt, vTr)) { + // Inserisco i nuovi triangoli + for ( int n = 0 ; n < int( vTr.size()) - 2 ; n += 3) { + int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] } ; + int nNewId[3] = { SurfB.AddVertex(vPt[nNewTriaVertId[0]]), + SurfB.AddVertex(vPt[nNewTriaVertId[1]]), + SurfB.AddVertex(vPt[nNewTriaVertId[2]]) } ; + int nNewTriaNum = SurfB.AddTriangle( nNewId) ; + SurfB.m_vTria[nNewTriaNum].nTempPart = 1 ; + bModif = true ; + } + } + + for ( int nL = 1 ; nL < int( vPolygons.size()) ; ++ nL) { + vPolygons[nL].Invert() ; + if ( Triangulate().Make( vPolygons[nL], vPt, vTr)) { + // Inserisco i nuovi triangoli + for ( int n = 0 ; n < int( vTr.size()) - 2 ; n += 3) { + int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] } ; + int nNewId[3] = { SurfB.AddVertex(vPt[nNewTriaVertId[0]]), + SurfB.AddVertex(vPt[nNewTriaVertId[1]]), + SurfB.AddVertex(vPt[nNewTriaVertId[2]]) } ; + int nNewTriaNum = SurfB.AddTriangle( nNewId) ; + SurfB.m_vTria[nNewTriaNum].nTempPart = - 1 ; + bModif = true ; + } + } + } + } + } + vInnerLoop.resize( 0) ; + } + } + // Se i triangoli delle superfici non si intersecano, una delle due è totalmente interna o esterna all'altra. + if ( ! bModif) { + int nVertNum = 0 ; + Point3d ptFirstV ; + int nCurVert = GetFirstVertex( ptFirstV) ; + int nInOutNum = 0 ; + while ( nInOutNum == 0 && nCurVert != SVT_NULL) { + int nTriaNum = - 1 ; + double dMinDist = DBL_MAX ; + for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) { + // Se il triangolo B non è valido, continuo + Triangle3d trTriaB ; + if ( ! ( SurfB.GetTriangle( nTB, trTriaB) && trTriaB.Validate( true))) + continue ; + DistPointTriangle DistCalculator( ptFirstV, trTriaB) ; + double dDist ; + DistCalculator.GetDist( dDist) ; + if ( dDist < dMinDist) { + nTriaNum = nTB ; + dMinDist = dDist ; + } + } + if ( nTriaNum >= 0) { + Triangle3d trTriaB ; + SurfB.GetTriangle( nTriaNum, trTriaB) ; + + if ( ( ptFirstV - trTriaB.GetP(0)) * trTriaB.GetN() < - EPS_SMALL) { + nInOutNum = 1 ; + } + else if ( ( ptFirstV - trTriaB.GetP(0)) * trTriaB.GetN() > EPS_SMALL) { + nInOutNum = - 1 ; + } + } + if ( nInOutNum == 0) { + nCurVert = GetNextVertex( nVertNum, ptFirstV) ; + ++ nVertNum ; + } + } + for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) { + m_vTria[nTA].nTempPart = nInOutNum ; + } + nVertNum = 0 ; + ptFirstV ; + nCurVert = SurfB.GetFirstVertex( ptFirstV) ; + nInOutNum = 0 ; + while ( nInOutNum == 0 && nCurVert != SVT_NULL) { + int nTriaNum = - 1 ; + double dMinDist = DBL_MAX ; + for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) { + // Se il triangolo A non è valido, continuo + Triangle3d trTriaA ; + if ( ! ( GetTriangle( nTA, trTriaA) && trTriaA.Validate( true))) + continue ; + DistPointTriangle DistCalculator( ptFirstV, trTriaA) ; + double dDist ; + DistCalculator.GetDist( dDist) ; + if ( dDist < dMinDist) { + nTriaNum = nTA ; + dMinDist = dDist ; + } + } + if ( nTriaNum >= 0) { + Triangle3d trTriaA ; + GetTriangle( nTriaNum, trTriaA) ; + + if ( ( ptFirstV - trTriaA.GetP( 0)) * trTriaA.GetN() < - EPS_SMALL) { + nInOutNum = 1 ; + } + else if ( ( ptFirstV - trTriaA.GetP(0)) * trTriaA.GetN() > EPS_SMALL) { + nInOutNum = - 1 ; + } + } + if ( nInOutNum == 0) { + nCurVert = SurfB.GetNextVertex( nVertNum, ptFirstV) ; + ++ nVertNum ; + } + } + for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) { + SurfB.m_vTria[nTB].nTempPart = nInOutNum ; + } + } + + bool bContinue = true ; + // Se avvenuta modifica, aggiorno tutto + if ( bModif) + bContinue = ( AdjustVertices() && DoCompacting() && SurfB.AdjustVertices() && SurfB.DoCompacting()) ; + // Triangoli sovrapposti + if ( bContinue) { + int nTriaNumA = GetTriangleSize() ; + int nTriaNumB = SurfB.GetTriangleSize() ; + for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) { + // Se il triangolo A non è valido, continuo + Triangle3d trTriaA ; + if ( ! ( GetTriangle( nTA, trTriaA) && trTriaA.Validate( true))) + continue ; + // Box del triangolo A + BBox3d b3dTriaA ; + trTriaA.GetLocalBBox( b3dTriaA) ; + bool bNewTriaA = true ; + for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) { + // Se il triangolo B non è valido, continuo + Triangle3d trTriaB ; + if ( ! ( SurfB.GetTriangle( nTB, trTriaB) && trTriaB.Validate( true))) + continue ; + // Box del triangolo B + BBox3d b3dTriaB ; + trTriaB.GetLocalBBox( b3dTriaB) ; + // Se i box non si sovrappongono, contiuno + if ( ! b3dTriaA.Overlaps( b3dTriaB)) + continue ; + // Se i triangoli sono sovrapposti + TRIA3DVECTOR vTriaAB ; + Point3d ptTempA, ptTempB ; + int nIntTypeAB = IntersTriaTria( trTriaA, trTriaB, ptTempA, ptTempB, vTriaAB) ; + if ( nIntTypeAB == ITTTS_OVERLAPS) { + bool bInvertB = trTriaA.GetN() * trTriaB.GetN() < 0. ; + m_vTria[nTA].nTempPart = ( m_vTria[nTA].nTempPart == 1 ? ( bInvertB ? -2 : 2) : m_vTria[nTA].nTempPart) ; + SurfB.m_vTria[nTB].nTempPart = ( SurfB.m_vTria[nTB].nTempPart == 1 ? ( bInvertB ? - 2 : 2) : SurfB.m_vTria[nTB].nTempPart) ; + } + } + } + return ( AdjustVertices() && DoCompacting() && SurfB.AdjustVertices() && SurfB.DoCompacting()) ; + } + + return true ; +} + +//---------------------------------------------------------------------------- +bool +SurfTriMesh::IdentifyParts( void) const +{ + for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i) { + // salto triangoli cancellati o già assegnati + if ( m_vTria[i].nIdVert[0] == SVT_DEL || + abs( m_vTria[i].nTempPart) != 1) + continue ; + // set di triangoli da aggiornare + set stTria ; + stTria.insert( i) ; + while ( ! stTria.empty()) { + // tolgo un triangolo dal set + const auto iIt = stTria.begin() ; + int nT = *iIt ; + stTria.erase( iIt) ; + // aggiorno i triangoli adiacenti + for ( int j = 0 ; j < 3 ; ++ j) { + int nAdjT = m_vTria[nT].nIdAdjac[j] ; + if ( nAdjT != SVT_NULL && m_vTria[nAdjT].nTempPart == 0) { + m_vTria[nAdjT].nTempPart = m_vTria[nT].nTempPart ; + stTria.insert( nAdjT) ; + } + } + } + } + + return true ; +} + +//---------------------------------------------------------------------------- +bool +SurfTriMesh::Add( const ISurfTriMesh& Other) +{ + SurfTriMesh SurfB ; + SurfB.CopyFrom( &Other) ; + IntersectTriMeshTriangle( SurfB) ; + IdentifyParts() ; + SurfB.IdentifyParts() ; + int nTriaNumA = GetTriangleSize() ; + for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) { + if ( m_vTria[nTA].nTempPart == 1 || m_vTria[nTA].nTempPart == - 2) + RemoveTriangle( nTA) ; + } + int nTriaNumB = SurfB.GetTriangleSize() ; + for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) { + if ( SurfB.m_vTria[nTB].nTempPart == - 1) { + int nNewVert[3] ; + for ( int nV = 0 ; nV < 3 ; ++ nV) { + nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ; + } + AddTriangle( nNewVert) ; + } + } + + return ( AdjustVertices() && DoCompacting()) ; +} + +//---------------------------------------------------------------------------- +bool +SurfTriMesh::Intersect( const ISurfTriMesh& Other) +{ + SurfTriMesh SurfB ; + SurfB.CopyFrom( &Other) ; + IntersectTriMeshTriangle( SurfB) ; + IdentifyParts() ; + SurfB.IdentifyParts() ; + int nTriaNumA = GetTriangleSize() ; + for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) { + if ( m_vTria[nTA].nTempPart == - 1 || m_vTria[nTA].nTempPart == - 2) + RemoveTriangle( nTA) ; + } + int nTriaNumB = SurfB.GetTriangleSize() ; + for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) { + if ( SurfB.m_vTria[nTB].nTempPart == 1) { + int nNewVert[3] ; + for ( int nV = 0 ; nV < 3 ; ++ nV) { + nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ; + } + AddTriangle( nNewVert) ; + } + } + + return ( AdjustVertices() && DoCompacting()) ; +} + +//---------------------------------------------------------------------------- +bool +SurfTriMesh::Subtract( const ISurfTriMesh& Other) +{ + SurfTriMesh SurfB ; + SurfB.CopyFrom( &Other) ; + IntersectTriMeshTriangle( SurfB) ; + IdentifyParts() ; + SurfB.IdentifyParts() ; + int nTriaNumA = GetTriangleSize() ; + for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) { + if ( m_vTria[nTA].nTempPart == 1 || m_vTria[nTA].nTempPart == 2) + RemoveTriangle( nTA) ; + } + int nTriaNumB = SurfB.GetTriangleSize() ; + for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) { + if ( SurfB.m_vTria[nTB].nTempPart == 1) { + int nNewVert[3] ; + for ( int nV = 0 ; nV < 3 ; ++ nV) { + nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ; + } + swap( nNewVert[1], nNewVert[2]) ; + AddTriangle( nNewVert) ; + } + } + + return ( AdjustVertices() && DoCompacting()) ; +} diff --git a/VolZmap.cpp b/VolZmap.cpp index 88ef179..1979695 100644 --- a/VolZmap.cpp +++ b/VolZmap.cpp @@ -622,40 +622,41 @@ VolZmap::CheckMapConnection( void) for ( int tInt = 0 ; tInt < int( m_Values[tMap][tDex].size()) ; ++ tInt) { m_Values[tMap][tDex][tInt].nCompo = 0 ; - // Controlli sui tratti di dexel non incidenti su nodi del reticolo - if ( tMap == 0) { - // Z degli estremi del segmento - double dZMin = m_Values[tMap][tDex][tInt].dMin ; - double dZMax = m_Values[tMap][tDex][tInt].dMax ; - // Indici k dei voxels in cui cadono le Z - int nKmin = int( floor( ( dZMin - 2 * EPS_SMALL) / m_dStep - 0.5)) ; - int nKmax = int( floor( ( dZMax + 2 * EPS_SMALL) / m_dStep - 0.5)) ; - // Se cadono nello stesso voxel imposto a -1 il valore della componente connessa - if ( nKmax - nKmin == 0) - m_Values[tMap][tDex][tInt].nCompo = -1 ; - } - else if ( tMap == 1) { - // X degli estremi del segmento - double dXMin = m_Values[tMap][tDex][tInt].dMin ; - double dXMax = m_Values[tMap][tDex][tInt].dMax ; - // Indici i dei voxels in cui cadono le X - int nImin = int( floor( ( dXMin - 2 * EPS_SMALL) / m_dStep - 0.5)) ; - int nImax = int( floor( ( dXMax + 2 * EPS_SMALL) / m_dStep - 0.5)) ; - // Se cadono nello stesso voxel imposto a -1 il valore della componente connessa - if ( nImax - nImin == 0) - m_Values[tMap][tDex][tInt].nCompo = -1 ; - } - else { - // Y degli estremi del segmento - double dYMin = m_Values[tMap][tDex][tInt].dMin ; - double dYMax = m_Values[tMap][tDex][tInt].dMax ; - // Indici j dei voxels in cui cadono le X - int nJmin = int( floor( ( dYMin - 2 * EPS_SMALL) / m_dStep - 0.5)) ; - int nJmax = int( floor( ( dYMax + 2 * EPS_SMALL) / m_dStep - 0.5)) ; - // Se cadono nello stesso voxel imposto a -1 il valore della componente connessa - if ( nJmax - nJmin == 0) - m_Values[tMap][tDex][tInt].nCompo = -1 ; - } + + //// Controlli sui tratti di dexel non incidenti su nodi del reticolo + // if ( tMap == 0) { + // // Z degli estremi del segmento + // double dZMin = m_Values[tMap][tDex][tInt].dMin ; + // double dZMax = m_Values[tMap][tDex][tInt].dMax ; + // // Indici k dei voxels in cui cadono le Z + // int nKmin = int( floor( ( dZMin - 2 * EPS_SMALL) / m_dStep - 0.5)) ; + // int nKmax = int( floor( ( dZMax + 2 * EPS_SMALL) / m_dStep - 0.5)) ; + // // Se cadono nello stesso voxel imposto a -1 il valore della componente connessa + // if ( nKmax - nKmin == 0) + // m_Values[tMap][tDex][tInt].nCompo = -1 ; + // } + // else if ( tMap == 1) { + // // X degli estremi del segmento + // double dXMin = m_Values[tMap][tDex][tInt].dMin ; + // double dXMax = m_Values[tMap][tDex][tInt].dMax ; + // // Indici i dei voxels in cui cadono le X + // int nImin = int( floor( ( dXMin - 2 * EPS_SMALL) / m_dStep - 0.5)) ; + // int nImax = int( floor( ( dXMax + 2 * EPS_SMALL) / m_dStep - 0.5)) ; + // // Se cadono nello stesso voxel imposto a -1 il valore della componente connessa + // if ( nImax - nImin == 0) + // m_Values[tMap][tDex][tInt].nCompo = -1 ; + // } + // else { + // // Y degli estremi del segmento + // double dYMin = m_Values[tMap][tDex][tInt].dMin ; + // double dYMax = m_Values[tMap][tDex][tInt].dMax ; + // // Indici j dei voxels in cui cadono le X + // int nJmin = int( floor( ( dYMin - 2 * EPS_SMALL) / m_dStep - 0.5)) ; + // int nJmax = int( floor( ( dYMax + 2 * EPS_SMALL) / m_dStep - 0.5)) ; + // // Se cadono nello stesso voxel imposto a -1 il valore della componente connessa + // if ( nJmax - nJmin == 0) + // m_Values[tMap][tDex][tInt].nCompo = -1 ; + // } } } } @@ -725,8 +726,8 @@ VolZmap::ExpandFromXInterval( IntContaier& IntCont) // Quote estreme del segmento lungo X double dMinX = m_Values[1][tDex][tInt].dMin ; double dMaxX = m_Values[1][tDex][tInt].dMax ; - double dMinDX = max( floor( ( dMinX - EPS_SMALL) / m_dStep - 0.5), 0.) ; - double dMaxDX = max( floor( ( dMaxX + EPS_SMALL) / m_dStep - 0.5), 0.) ; + double dMinDX = max( floor( ( dMinX - 2 * EPS_SMALL) / m_dStep - 0.5), 0.) ; + double dMaxDX = max( floor( ( dMaxX + 2 * EPS_SMALL) / m_dStep - 0.5), 0.) ; // Indici estremi dei dei dexel ortogonali // che possono intersecare il segmento di partenza int tStartI = min( int( dMinDX), ( m_nNx[0] - 1)) ; @@ -746,8 +747,8 @@ VolZmap::ExpandFromXInterval( IntContaier& IntCont) // Se i segmenti si incrociano e il nuovo trovato non // ha già un indice assegnato, assegno l'indice e // aggiungo l'intervallo trovato allo stack. - if ( dZmin - EPS_SMALL < dZ && - dZmax + EPS_SMALL > dZ && + if ( dZmin - 2 * EPS_SMALL < dZ && + dZmax + 2 * EPS_SMALL > dZ && m_Values[0][tGrIndex1 * m_nNx[0] + tI][tIntZ].nCompo == 0) { m_Values[0][tGrIndex1 * m_nNx[0] + tI][tIntZ].nCompo = m_Values[1][tDex][tInt].nCompo ; IntervalIndexes NewInterval ; @@ -766,8 +767,8 @@ VolZmap::ExpandFromXInterval( IntContaier& IntCont) // Se i segmenti si incrociano e il nuovo trovato non // ha già un indice assegnato, assegno l'indice e // aggiungo l'intervallo trovato allo stack. - if ( dYmin - EPS_SMALL < dY && - dYmax + EPS_SMALL > dY && + if ( dYmin - 2 * EPS_SMALL < dY && + dYmax + 2 * EPS_SMALL > dY && m_Values[2][tI * m_nNx[2] + tGrIndex2][tIntY].nCompo == 0) { m_Values[2][tI * m_nNx[2] + tGrIndex2][tIntY].nCompo = m_Values[1][tDex][tInt].nCompo ; IntervalIndexes NewInterval ; @@ -796,8 +797,8 @@ VolZmap::ExpandFromYInterval( IntContaier& IntCont) // Quote estreme del segmento lungo Y double dMinY = m_Values[2][tDex][tInt].dMin ; double dMaxY = m_Values[2][tDex][tInt].dMax ; - double dMinDY = max( floor( ( dMinY - EPS_SMALL) / m_dStep - 0.5), 0.) ; - double dMaxDY = max( floor( ( dMaxY + EPS_SMALL) / m_dStep - 0.5), 0.) ; + double dMinDY = max( floor( ( dMinY - 2 * EPS_SMALL) / m_dStep - 0.5), 0.) ; + double dMaxDY = max( floor( ( dMaxY + 2 * EPS_SMALL) / m_dStep - 0.5), 0.) ; // Indici estremi dei dei dexel ortogonali // che possono intersecare il segmento di partenza int tStartJ = min( int( dMinDY), ( m_nNy[0] - 1)) ; @@ -817,8 +818,8 @@ VolZmap::ExpandFromYInterval( IntContaier& IntCont) // Se i segmenti si incrociano e il nuovo trovato non // ha già un indice assegnato, assegno l'indice e // aggiungo l'intervallo trovato allo stack. - if ( dZmin - EPS_SMALL < dZ && - dZmax + EPS_SMALL > dZ && + if ( dZmin - 2 * EPS_SMALL < dZ && + dZmax + 2 * EPS_SMALL > dZ && m_Values[0][tJ * m_nNx[0] + tGrIndex2][tIntZ].nCompo == 0) { m_Values[0][tJ * m_nNx[0] + tGrIndex2][tIntZ].nCompo = m_Values[2][tDex][tInt].nCompo ; IntervalIndexes NewInterval ; @@ -837,8 +838,8 @@ VolZmap::ExpandFromYInterval( IntContaier& IntCont) // Se i segmenti si incrociano e il nuovo trovato non // ha già un indice assegnato, assegno l'indice e // aggiungo l'intervallo trovato allo stack. - if ( dXmin - EPS_SMALL < dX && - dXmax + EPS_SMALL > dX && + if ( dXmin - 2 * EPS_SMALL < dX && + dXmax + 2 * EPS_SMALL > dX && m_Values[1][tGrIndex1 * m_nNx[1] + tJ][tIntX].nCompo == 0) { m_Values[1][tGrIndex1 * m_nNx[1] + tJ][tIntX].nCompo = m_Values[2][tDex][tInt].nCompo ; IntervalIndexes NewInterval ; @@ -867,8 +868,8 @@ VolZmap::ExpandFromZInterval( IntContaier& IntCont) // Quote estreme del segmento lungo Z double dMinZ = m_Values[0][tDex][tInt].dMin ; double dMaxZ = m_Values[0][tDex][tInt].dMax ; - double dMinDZ = max( floor( ( dMinZ - EPS_SMALL) / m_dStep - 0.5), 0.) ; - double dMaxDZ = max( floor( ( dMaxZ + EPS_SMALL) / m_dStep - 0.5), 0.) ; + double dMinDZ = max( floor( ( dMinZ - 2 * EPS_SMALL) / m_dStep - 0.5), 0.) ; + double dMaxDZ = max( floor( ( dMaxZ + 2 * EPS_SMALL) / m_dStep - 0.5), 0.) ; // Indici estremi dei dexel ortogonali // che possono intersecare il segmento di partenza int tStartK = min( int( dMinDZ), ( m_nNy[1] - 1)) ; @@ -888,9 +889,9 @@ VolZmap::ExpandFromZInterval( IntContaier& IntCont) // Se i segmenti si incrociano e il nuovo trovato non // ha già un indice assegnato, assegno l'indice e // aggiungo l'intervallo trovato allo stack. - if ( dXmin - EPS_SMALL < dX && - dXmax + EPS_SMALL > dX && - m_Values[1][tK * m_nNx[1] + tGrIndex2][tIntX].nCompo == 0) { + if ( dXmin - 2 * EPS_SMALL < dX && + dXmax + 2 * EPS_SMALL > dX && + m_Values[1][tK * m_nNx[1] + tGrIndex2][tIntX].nCompo == 0) { m_Values[1][tK * m_nNx[1] + tGrIndex2][tIntX].nCompo = m_Values[0][tDex][tInt].nCompo ; IntervalIndexes NewInterval ; NewInterval.tMap = 1 ; @@ -908,8 +909,8 @@ VolZmap::ExpandFromZInterval( IntContaier& IntCont) // Se i segmenti si incrociano e il nuovo trovato non // ha già un indice assegnato, assegno l'indice e // aggiungo l'intervallo trovato allo stack. - if ( dYmin - EPS_SMALL < dY && - dYmax + EPS_SMALL > dY && + if ( dYmin - 2 * EPS_SMALL < dY && + dYmax + 2 * EPS_SMALL > dY && m_Values[2][tGrIndex1 * m_nNx[2] + tK][tIntY].nCompo == 0) { m_Values[2][tGrIndex1 * m_nNx[2] + tK][tIntY].nCompo = m_Values[0][tDex][tInt].nCompo ; IntervalIndexes NewInterval ; diff --git a/VolZmapVolume.cpp b/VolZmapVolume.cpp index 9abd4c3..9d0cdb7 100644 --- a/VolZmapVolume.cpp +++ b/VolZmapVolume.cpp @@ -120,8 +120,8 @@ VolZmap::SubtractIntervals( int nGrid, int nI, int nJ, // Elimino residui di intervalli inutili for ( int i = 0 ; i < int( vDexel.size()) ; ++ i) { - int nMinN = int( floor( ( vDexel[i].dMax + 2 * EPS_SMALL - 0.5 * m_dStep) / m_dStep)) ; - int nMaxN = int( floor( ( vDexel[i].dMin - 2 * EPS_SMALL - 0.5 * m_dStep) / m_dStep)) ; + int nMaxN = int( floor( ( vDexel[i].dMax + 2 * EPS_SMALL - 0.5 * m_dStep) / m_dStep)) ; + int nMinN = int( floor( ( vDexel[i].dMin - 2 * EPS_SMALL - 0.5 * m_dStep) / m_dStep)) ; if ( nMinN == nMaxN) { vDexel.erase( vDexel.begin() + i) ; -- i ;