From a58c3fe5af57176bfe783d79fae23a3f5c9b8bf3 Mon Sep 17 00:00:00 2001 From: Dario Sassi Date: Thu, 26 Apr 2018 07:36:37 +0000 Subject: [PATCH] EgtGeomKernel : - correzioni a visualizzazione Zmap. --- VolZmap.h | 4 +- VolZmapGraphics.cpp | 730 +++++++++++++++++++------------------------- 2 files changed, 313 insertions(+), 421 deletions(-) diff --git a/VolZmap.h b/VolZmap.h index edb7b59..349a157 100644 --- a/VolZmap.h +++ b/VolZmap.h @@ -169,13 +169,13 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool CalcDexelPrisms( int nPos1, int nPos2, TRIA3DEXLIST& lstTria) const ; bool AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Point3d& ptQ, const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DEXLIST& lstTria) const ; - bool ProcessCube( int nVoxI, int nVoxJ, int nVoxK, int& nCubeIndex, TRIA3DEXLIST& lstTria) const ; - bool ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, TriHolder& triHold, bool bExt) const ; + bool ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, VoxelContainer& VoxCont, bool bEnh) const ; bool ProcessCell( int nGrid, int nCellI, int nCellJ, const Plane3d& plPlane, std::vector& vLine) const ; bool ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, VoxelContainer& vVox) const ; bool ExtMarchingCubes( std::vector& vVox, TRIA3DEXLIST& lstTria, bool bEnh) const ; bool RegulateFeaturesChain( std::vector& vVecVox) const ; bool CreateSharpFeatureTriangle( int nBlock, const VoxelContainer& vVox, TriHolder& triHold) const ; + bool CreateSharpFeatureTriangle( const VoxelContainer& vVox, TriHolder& triHold) const ; bool FlipEdgesII( TriHolder& TriHold) const ; bool FlipEdgesBB( TriaMatrix& InterTria) const ; bool IsThereMat( int nI, int nJ, int nK) const ; diff --git a/VolZmapGraphics.cpp b/VolZmapGraphics.cpp index 3729038..fbaf5f7 100644 --- a/VolZmapGraphics.cpp +++ b/VolZmapGraphics.cpp @@ -513,7 +513,7 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DEXLIST if ( bNormN && bNormV && vtN * vtNV < 0.7) VecTriHold[t][t1].vCompoTria[t2][t3].SetVertexNorm( nV, vtN) ; } - // aggiungo triangolo alla lista + // aggiungo triangolo alla lista vLstTria[nModifiedBlocks[t]].emplace_back( VecTriHold[t][t1].vCompoTria[t2][t3]) ; } } @@ -758,7 +758,7 @@ VolZmap::AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Poin // il riconoscimento delle sharp-feature bEnh deve valere true. // I triangoli formanti sharp-feature vengono messi nel TriHolder, gli altri nella lista. bool -VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, TriHolder& triHold, bool bEnh) const +VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, VoxelContainer& VoxCont, bool bEnh) const { // Calcolo il numero di voxel lungo X,Y e Z int nVoxNumX = int( m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ; @@ -817,93 +817,85 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr bReg = false ; } - // Determino il numero di componenti connesse nel voxel in caso di configurazione standard + // Determino il numero di componenti connesse nel voxel in caso di configurazione standard int nComponents = TriangleTableEn[nIndex][1][0] ; - // Matrici di campi vettoriali: - // CompoVert[i] ha i vertici della base del triangle fan della (i+1)-esima componente connessa; - // CompoTriVert[i] ha i vertici di tutti i triangoli, nel caso di assenza di sharp feature, - // della (i+1)-esima componente connessa. + // Matrici di campi vettoriali: + // CompoVert[i] ha i vertici della base del triangle fan della (i+1)-esima componente connessa; + // CompoTriVert[i] ha i vertici di tutti i triangoli, nel caso di assenza di sharp feature, + // della (i+1)-esima componente connessa. AppliedVector CompoVert[4][7] ; AppliedVector CompoTriVert[4][17] ; - // Arrey numero di vertici della base del fan per componente - // connessa: nVertComp[i] contiene il numero di vertici - // della base del fan della (i+1)-esima componente connessa. + // Arrey numero di vertici della base del fan per componente + // connessa: nVertComp[i] contiene il numero di vertici + // della base del fan della (i+1)-esima componente connessa. int nVertComp[4] ; int nExtTabOff = nComponents ; int nStdTabOff = 0 ; - // Carico le matrici CompoVert e CompoTriVert + // Carico le matrici CompoVert e CompoTriVert for ( int nComp = 0 ; nComp < nComponents ; ++ nComp) { - // Numero vertici per componenti + // Numero vertici per componenti nVertComp[nComp] = TriangleTableEn[nIndex][1][nComp+1] ; - // Riempio il nCompCount-esimo vettore di vertici della base del fan + // Riempio il nCompCount-esimo vettore di vertici della base del fan for ( int nVertCount = 0 ; nVertCount < nVertComp[nComp] ; ++ nVertCount) CompoVert[nComp][nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1]] ; - // Riempio il nCompCount-esimo vettore di vertici dei triangoli in assenza di - // sharp feature: in una mesh di triangoli con n vertici vi sono n - 2 triangoli. + // Riempio il nCompCount-esimo vettore di vertici dei triangoli in assenza di + // sharp feature: in una mesh di triangoli con n vertici vi sono n - 2 triangoli. for ( int nVert = 0 ; nVert < 3 * ( nVertComp[nComp] - 2) ; nVert += 3) { CompoTriVert[nComp][nVert] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVert+2]] ; CompoTriVert[nComp][nVert+1] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVert+1]] ; CompoTriVert[nComp][nVert+2] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVert]] ; } - // Aggiorno gli offsets per raggiungere i - // vertici della componente successiva. + // Aggiorno gli offsets per raggiungere i + // vertici della componente successiva. nExtTabOff += nVertComp[nComp] ; nStdTabOff += 3 * ( nVertComp[nComp] - 2) ; } - // Test sulla topologia: dal momento che il nostro test - // si fonda sugli angoli compresi fra le normali, esso ha - // senso solo se il campo è regolare. - + // Test sulla topologia: richiede campo regolare perchè si fonda su angoli tra normali if ( bReg) { + // Configurazione 3 if ( nAllConfig[nIndex] == 3) { - // Verifico se i versori delle componenti sono tutti - // più o meno concordi (per esserlo non devono esserci - // due vettori di una medesima componente con prodotto - // scalare inferiore a 0.7). - Vector3d vtCmpAvg0, vtCmpAvg1 ; - bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ; - bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ; - - // Se i versori di entrambe le componenti sono concordi - // ha senso parlare di vettori medi, altrimenti non ha - // senso. Se non ha senso parlare di vettori medi non - // ha senso parlare di prodotti scalari fra loro, - // quindi pongo il loro prodotto a un valore assurdo -2 - // (il prodotto scalare fra versori ha modulo non superiore - // a uno). - double dScProd = - 2 ; - if ( bTest0 && bTest1) - dScProd = vtCmpAvg0 * vtCmpAvg1 ; - - double dThreshold = 0.7 ; + // Test sulla topologia + double dDotSum = 0 ; + for ( int nFV = 0 ; nFV < 3 ; ++ nFV) { + for ( int nTV = 0 ; nTV < 2 ; ++ nTV) { + dDotSum += CompoVert[0][nFV].vtVec * CompoVert[1][nTV].vtVec ; + } + } + for ( int nFVI = 0 ; nFVI < 2 ; ++ nFVI) { + for ( int nFVJ = nFVI + 1 ; nFVJ < 3 ; ++ nFVJ) { + dDotSum -= CompoVert[0][nFVI].vtVec * CompoVert[0][nFVJ].vtVec ; + } + } + for ( int nTVI = 0 ; nTVI < 2 ; ++ nTVI) { + for ( int nTVJ = nTVI + 1 ; nTVJ < 3 ; ++ nTVJ) { + dDotSum -= CompoVert[1][nTVI].vtVec * CompoVert[1][nTVJ].vtVec ; + } + } + bool bTestOnSum = dDotSum > - EPS_SMALL ; // Si passa alla seconda topologia - if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) { - // Ricerca del caso corrispondente della nuova topologia + if ( bTestOnSum) { + // Ricerca del caso corrispondente della nuova topologia int nt = 0 ; - while ( nIndexVsIndex3[nt][0] != nIndex) ++ nt ; - int nRotCase = nIndexVsIndex3[nt][1] ; - // Riaggiorno numero di componenti + // Aggiorno numero di componenti nComponents = Cases3Plus[nRotCase][1][0] ; // Riaggiorno gli offsets nExtTabOff = nComponents ; nStdTabOff = 0 ; // Modifico le matrici for ( int nC = 1 ; nC <= nComponents ; ++ nC) { - // Numero vertici per componenti nVertComp[nC - 1] = Cases3Plus[nRotCase][1][nC] ; - // Matrice dei vertici della base del fan for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert) CompoVert[nC - 1][nFanVert] = VecField[Cases3Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ; @@ -920,72 +912,105 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr } } + // Configurazione 6 else if ( nAllConfig[nIndex] == 6) { - - // Procedura analoga a quella della configurazione 3 - Vector3d vtCmpAvg0, vtCmpAvg1 ; - bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.7) ; - bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ; - double dScProd = - 2 ; - if ( bTest0 && bTest1) - dScProd = vtCmpAvg0 * vtCmpAvg1 ; - - double dThreshold = 0.7 ; - // Si passa alla seconda topologia - if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) { + // Test sulla topologia + double dDotSum = 0 ; + for ( int nFV = 0 ; nFV < 4 ; ++ nFV) { + for ( int nTV = 0 ; nTV < 3 ; ++ nTV) { + dDotSum += CompoVert[0][nFV].vtVec * CompoVert[1][nTV].vtVec ; + } + } + for ( int nFVI = 0 ; nFVI < 3 ; ++ nFVI) { + for ( int nFVJ = nFVI + 1 ; nFVJ < 4 ; ++ nFVJ) { + dDotSum -= CompoVert[0][nFVI].vtVec * CompoVert[0][nFVJ].vtVec ; + } + } + for ( int nTVI = 0 ; nTVI < 2 ; ++ nTVI) { + for ( int nTVJ = nTVI + 1 ; nTVJ < 3 ; ++ nTVJ) { + dDotSum -= CompoVert[1][nTVI].vtVec * CompoVert[1][nTVJ].vtVec ; + } + } + bool bTestOnSum = dDotSum > - 4 ; + // Si deve passare alla seconda topologia + if ( bTestOnSum) { // Ricerca del caso corrispondente della nuova topologia int nt = 0 ; while ( nIndexVsIndex6[nt][0] != nIndex) ++ nt ; int nRotCase = nIndexVsIndex6[nt][1] ; - // Riaggiorno numero di componenti - nComponents = 1 ; - // Riaggiorno gli offsets - nExtTabOff = nComponents ; - nStdTabOff = 0 ; - // Modifico le matrici - for ( int nC = 1 ; nC <= nComponents ; ++ nC) { - // Numero vertici per componenti - nVertComp[nC - 1] = 7 ; - // Matrici dei vertici dei triangoli in assenza di sharp feature - for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) { - CompoTriVert[nC - 1][nTriVert] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ; - CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ; - CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert]] ; + + // Costruzione dei triangoli + for ( int TriIndex = 0 ; TriIndex < 15 ; TriIndex += 3) { + // Indici vertici + int i0 = Cases6Plus[nRotCase][TriIndex + 2] ; + int i1 = Cases6Plus[nRotCase][TriIndex + 1] ; + int i2 = Cases6Plus[nRotCase][TriIndex] ; + // Costruzione triangolo + Triangle3dEx CurrentTriangle ; + CurrentTriangle.Set( VecField[i0].ptPApp, VecField[i1].ptPApp, VecField[i2].ptPApp) ; + CurrentTriangle.Validate( true) ; + CurrentTriangle.SetVertexNorm( 0, VecField[i0].vtVec) ; + CurrentTriangle.SetVertexNorm( 1, VecField[i1].vtVec) ; + CurrentTriangle.SetVertexNorm( 2, VecField[i2].vtVec) ; + // Setto il numero di utensile (conta solo positivo, nullo o negativo) + int nTool0 = Clamp( VecField[i0].nPropIndex, -1, 1) ; + int nTool1 = Clamp( VecField[i1].nPropIndex, -1, 1) ; + int nTool2 = Clamp( VecField[i2].nPropIndex, -1, 1) ; + if ( nTool0 == nTool1 || nTool0 == nTool2) + CurrentTriangle.SetGrade( nTool0) ; + else if ( nTool1 == nTool2) + CurrentTriangle.SetGrade( nTool1) ; + // Valido il triangolo e setto le normali del campo vettoriale ai corrispondenti vertici + if ( CurrentTriangle.Validate( true)) { + Vector3d vtVertNorm = VecField[i0].vtVec ; + if ( CurrentTriangle.GetN() * vtVertNorm > 0.6) + CurrentTriangle.SetVertexNorm( 0, vtVertNorm) ; + else + CurrentTriangle.SetVertexNorm( 0, CurrentTriangle.GetN()) ; + vtVertNorm = VecField[i1].vtVec ; + if ( CurrentTriangle.GetN() * vtVertNorm > 0.6) + CurrentTriangle.SetVertexNorm( 1, vtVertNorm) ; + else + CurrentTriangle.SetVertexNorm( 1, CurrentTriangle.GetN()) ; + vtVertNorm = VecField[i2].vtVec ; + if ( CurrentTriangle.GetN() * vtVertNorm > 0.6) + CurrentTriangle.SetVertexNorm( 2, vtVertNorm) ; + else + CurrentTriangle.SetVertexNorm( 2, CurrentTriangle.GetN()) ; } - // Aggiorno gli offsets per raggiungere i vertici della componente successiva. - nExtTabOff += nVertComp[nC - 1] ; - nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; - } + // Riporto il triangolo nel sistema in cui è immerso quello dello Zmap + CurrentTriangle.ToGlob( m_MapFrame) ; + // Aggiungo alla lista + lstTria.emplace_back( CurrentTriangle) ; + } return true ; } } - + // Configurazione 10 else if ( nAllConfig[nIndex] == 10) { - // Verifico se i versori delle componenti sono tutti - // più o meno concordi (per esserlo non devono esserci - // due vettori di una medesima componente con prodotto - // scalare inferiore a 0). decidere se 0.0 o 0.7 + // Verifico concordanza tra i versori di una stessa componente + // (ogni coppia di vettori di una medesima componente deve avere prodotto scalare non inferiore a 0.0) Vector3d vtCmpAvg0, vtCmpAvg1 ; - bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ; - bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ; + bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.0) ; + bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1, 0.0) ; // Si passa alla seconda topologia - if ( ! ( bTest0 && bTest1)) { + if ( ! bTest0 || ! bTest1) { // Ricerca del caso corrispondente della nuova topologia int nt = 0 ; while ( nIndexVsIndex10[nt][0] != nIndex) ++ nt ; - int nRotCase = nIndexVsIndex10[nt][1] ; // Riaggiorno gli offsets nExtTabOff = 2 ; nStdTabOff = 0 ; // Modifico le matrici + int nRotCase = nIndexVsIndex10[nt][1] ; for ( int nC = 1 ; nC <= 2 ; ++ nC) { // Numero vertici per componenti nVertComp[nC - 1] = Cases10Plus[nRotCase][1][nC] ; // Matrice dei vertici della base del fan for ( int nFanVert = 0 ; nFanVert < 4 ; ++ nFanVert) - CompoVert[nC - 1][nFanVert] = VecField[Cases10Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ; + CompoVert[nC - 1][nFanVert] = VecField[Cases10Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ; // Matrici dei vertici dei triangoli in assenza di sharp feature for ( int nTriVert = 0 ; nTriVert < 6 ; nTriVert += 3) { CompoTriVert[nC - 1][nTriVert] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ; @@ -1000,35 +1025,34 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr } } + Voxel CurrVox ; + CurrVox.nNumComp = 0 ; + // Numero di feature nel voxel: al più vi è una feature per componente connessa. - int nInnerFeatureInVoxel = 0 ; - int nBorderFeatureInVoxel = 0 ; + int nFeatureInVoxel = 0 ; // Ciclo sulle componenti for ( int nComp = 0 ; nComp < nComponents ; ++ nComp) { - + int nFeatureType = NO_FEATURE ; // Se i componenti sono regolari valuto le normali per stabilire se eseguire ExtMC o MC if ( bReg) nFeatureType = TestOnNormal( CompoVert[nComp], nVertComp[nComp]) ; - // Flag ExtMC - bool bExtMC = ( nFeatureType != NO_FEATURE) && bEnh ; - // Extended MC - if ( bExtMC) { + if ( nFeatureType != NO_FEATURE) { // Passo al sistema di riferimento del baricentro Point3d ptGravityCenter( 0, 0, 0) ; for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) - ptGravityCenter = ptGravityCenter + CompoVert[nComp][ni].ptPApp ; + ptGravityCenter += CompoVert[nComp][ni].ptPApp ; ptGravityCenter = ptGravityCenter / nVertComp[nComp] ; Vector3d vtTrasf[12] ; for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) vtTrasf[ni] = CompoVert[nComp][ni].ptPApp - ptGravityCenter ; - // Preparo le matrici per il sistema + // Preparo le matrici per il sistema typedef Eigen::Matrix dSystemMatrix ; typedef Eigen::Matrix dSystemVector ; typedef Eigen::JacobiSVD DecomposerSVD ; @@ -1036,12 +1060,29 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr dSystemMatrix dMatrixN( nVertComp[nComp], 3) ; dSystemVector dKnownVector( nVertComp[nComp], 1) ; - // Caso generale + // medio le normali adiacenti molto vicine (delta angolare inferiore a 22.5 deg) + Vector3d vtNorm[12] ; + for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) + vtNorm[ni] = CompoVert[nComp][ni].vtVec ; for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - dMatrixN( ni, 0) = CompoVert[nComp][ni].vtVec.x ; - dMatrixN( ni, 1) = CompoVert[nComp][ni].vtVec.y ; - dMatrixN( ni, 2) = CompoVert[nComp][ni].vtVec.z ; - dKnownVector( ni) = CompoVert[nComp][ni].vtVec * vtTrasf[ni] ; + int nj = ( ni + 1) % nVertComp[nComp] ; + if ( vtNorm[ni] * vtNorm[nj] > 0.92) { + Vector3d vtNI = ( 0.6 * vtNorm[ni] + 0.4 * vtNorm[nj]) ; + Vector3d vtNJ = ( 0.4 * vtNorm[ni] + 0.6 * vtNorm[nj]) ; + vtNorm[ni] = vtNI ; + vtNorm[ni].Normalize() ; + vtNorm[nj] = vtNJ ; + vtNorm[nj].Normalize() ; + ++ ni ; + } + } + + // Definisco la matrice del sistema + for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { + dMatrixN( ni, 0) = vtNorm[ni].x ; + dMatrixN( ni, 1) = vtNorm[ni].y ; + dMatrixN( ni, 2) = vtNorm[ni].z ; + dKnownVector( ni) = vtNorm[ni] * vtTrasf[ni] ; } // calcolo SVD @@ -1058,318 +1099,97 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // risolvo il sistema con SVD, quindi ai minimi quadrati dSystemVector dUnknownVector( 3, 1) ; dUnknownVector = svd.solve( dKnownVector) ; - // Vettore Baricentro-Feature Vector3d vtFeature( dUnknownVector( 0), dUnknownVector( 1), dUnknownVector( 2)) ; - // Esprimo la soluzione nel sistema di riferimento z-map Point3d ptSol = ptGravityCenter + vtFeature ; - // Limito la feature entro i secondi voxel vicini, - // nel caso essa esca dal reticolo la limito entro una - // distanza dal baricentro pari alla diagonale del voxel. - bool bOutside = false ; - - int nFtVxI, nFtVxJ, nFtVxK ; - if ( GetPointVoxel( ptSol, nFtVxI, nFtVxJ, nFtVxK)) { - - if ( abs( nFtVxI - nVoxI) > 2 || - abs( nFtVxJ - nVoxJ) > 2 || - abs( nFtVxK - nVoxK) > 2) - bOutside = true ; - } - else { - double dDistFeature = vtFeature.Len() ; - const double MAX_DIST = sqrt( 3) * N_DEXVOXRATIO * m_dStep ; - bOutside = ( dDistFeature > MAX_DIST) ; - } - - // Costruisco triangoli del fan - TRIA3DVECTOR triContainer ; - for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - // Il triangolo è pronto - Triangle3d CurrentTriangle ; - CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptPApp, CompoVert[nComp][ni].ptPApp) ; - CurrentTriangle.Validate( true) ; - // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; - } - - // Controllo delle inversioni dei triangoli - bool bDangerInversion = false ; - - // Caso ventaglio con tre vertici di base - if ( nVertComp[nComp] == 3) { - - // Controllo se esiste almeno un triangolo con normale avente prodotto scalare - // negativo con la normale di almeno uno dei vertici di base del ventaglio. - bool bInversione = false ; - for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - - int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - - double dDI = triContainer[ni].GetN() * CompoVert[nComp][ni].vtVec ; - double dDJ = triContainer[ni].GetN() * CompoVert[nComp][nj].vtVec ; - - if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) { - bInversione = true ; + Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; + if ( nFeatureType == EDGE && vtNullSpace.Normalize() && + ! IsPointInsideVoxelApprox( nVoxI, nVoxJ, nVoxK, ptSol)) { + for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { + Vector3d vtBaseVert = CompoVert[nComp][ni].ptPApp - ptSol ; + Vector3d vtDist = vtBaseVert - vtBaseVert * vtNullSpace * vtNullSpace ; + double dMaxDist = 0.005 * N_DEXVOXRATIO * m_dStep ; + double bau = vtDist.Len() ; + if ( vtDist.SqLen() < dMaxDist * dMaxDist) { + ptSol = CompoVert[nComp][ni].ptPApp ; break ; - } - } - - // Se tale triangolo esiste procedo - if ( bInversione) { - - // Conto le coppie di normali con angolo compreso maggiore di 90° - int nNegDot = 0 ; - for ( int ni = 0 ; ni < nVertComp[nComp] - 1 ; ++ ni) { - for ( int nj = ni + 1 ; nj < nVertComp[nComp] ; ++ nj) { - if ( CompoVert[nComp][ni].vtVec * CompoVert[nComp][nj].vtVec < - EPS_SMALL) - nNegDot ++ ; - } - } - - if ( nNegDot == nVertComp[nComp] - 1) { - // Cerco se esiste un punto in cui la normale ha prodotto scalre negativo - // con le normali di entrambi i triangoli che lo contengono - bool bInversione2 = false ; - for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - int nj = ( ni == 0 ? nVertComp[nComp] - 1 : ni - 1) ; - double dDLast = triContainer[ni].GetN() * CompoVert[nComp][ni].vtVec ; - double dDPrev = triContainer[nj].GetN() * CompoVert[nComp][ni].vtVec ; - if ( ( dDLast < EPS_SMALL && dDPrev < EPS_SMALL) || - ( dDLast < - 0.75 || dDPrev < - 0.75)) { - bInversione2 = true ; - break ; - } - } - - // Se tale normale esiste - if ( bInversione2) { - - // Soluzione del sistema nel sistema Zmap - Point3d ptSolZMapFrame = ptSol ; - - // Se la soluzione non cade nel voxel di appartenenza vedo se può - // essere riportata dentro muovendosi lungo la linea di feature. - if ( ! IsPointInsideVoxelApprox( nVoxI, nVoxJ, nVoxK, ptSolZMapFrame)) { - - Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; - double dParInt1, dParInt2 ; - Point3d ptVoxMin( ( nVoxI * N_DEXVOXRATIO + 0.5) * m_dStep, - ( nVoxJ * N_DEXVOXRATIO + 0.5) * m_dStep, ( nVoxK * N_DEXVOXRATIO + 0.5) * m_dStep) ; - Point3d ptVoxMax( ( ( nVoxI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, - ( ( nVoxJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( nVoxK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ; - // Caso in cui può essere riportata dentro: se il voxel in cui cade la feature è pieno - // la riporto muovendola lungo la sua linea - if ( IntersLineBox( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { - - triContainer.resize( 0) ; - - double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 : - dParInt2 + ( dParInt1 - dParInt2) / 100 ; - - Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ; - - ptSol = ptNewSol ; - - // Costruisco triangoli di prova - for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - // Il triangolo è pronto - Triangle3d CurrentTriangle ; - CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptPApp, CompoVert[nComp][ni].ptPApp) ; - CurrentTriangle.Validate( true) ; - // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; - } - } - - // Caso in cui non può essere riportata dentro: - // si passa alla routine standard se il voxel - // in cui cade è pieno. - else { - - int nAdjVoxI, nAdjVoxJ, nAdjVoxK ; - if ( GetPointVoxel( ptSolZMapFrame, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) { - - // Classificazione del voxel adiacente - int nAdjIndex = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ; - // Se il voxel è pieno - if ( EdgeTable[nAdjIndex] != 0) - bDangerInversion = true ; - } - } - } - } - } - } - } - - // Ventaglio con base a quattro vertici - else if ( nVertComp[nComp] == 4) { - // Controllo se esiste almeno un triangolo con normale avente prodotto scalare - // negativo con la normale di almeno uno dei vertici di base del ventaglio. - bool bInversione = false ; - for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - - int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - - double dDI = triContainer[ni].GetN() * CompoVert[nComp][ni].vtVec ; - double dDJ = triContainer[ni].GetN() * CompoVert[nComp][nj].vtVec ; - - if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) - bInversione = true ; - } - - // Se tale triangolo esiste continuo i test - if ( bInversione) { - // Se la feature non cade nel suo voxel - if ( ! IsPointInsideVoxelApprox( nVoxI, nVoxJ, nVoxK, ptSol)) { - - Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; - double dParInt1, dParInt2 ; - Point3d ptVoxMin( ( nVoxI * N_DEXVOXRATIO + 0.5) * m_dStep, - ( nVoxJ * N_DEXVOXRATIO + 0.5) * m_dStep, ( nVoxK * N_DEXVOXRATIO + 0.5) * m_dStep) ; - Point3d ptVoxMax( ( ( nVoxI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, - ( ( nVoxJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( nVoxK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ; - // Se è possibile riportarla dentro e il voxel in cui cade è pieno, la riporto nel suo voxel - // lungo la sua linea - if ( IntersLineBox( ptSol, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { - - triContainer.resize( 0) ; - - double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 : - dParInt2 + ( dParInt1 - dParInt2) / 100 ; - - Point3d ptNewSol = ptSol + dPar * vtNullSpace ; - ptSol = ptNewSol ; - - // Costruisco triangoli di prova - for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - // Il triangolo è pronto - Triangle3d CurrentTriangle ; - CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptPApp, CompoVert[nComp][ni].ptPApp) ; - CurrentTriangle.Validate( true) ; - // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; - } - } - - // Se non è possibile riportarla dentro e il voxel in - // cui cade è pieno passo alla routine standard - else { - int nAdjVoxI, nAdjVoxJ, nAdjVoxK ; - if ( GetPointVoxel( ptSol, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) { - // Classificazione del voxel adiacente - int nAdjIndex = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ; - // Se il voxel è pieno - if ( EdgeTable[nAdjIndex] != 0) - bDangerInversion = true ; - } - } - } - } - } - - // Verifico se il fan di triangoli è piano - bool bPlane = true ; - int nContSize = int( triContainer.size()) ; - if ( nContSize > 0) { - Vector3d vtI = triContainer[0].GetN() ; - for ( int nj = 1 ; nj < nContSize ; ++ nj) { - Vector3d vtJ = triContainer[nj].GetN() ; - if ( ! AreSameVectorApprox( vtI, vtJ)) { - bPlane = false ; - break ; - } - } - } - - // Se feature nei limiti e triangoli non in un piano, confermo ExtMC - if ( ! bOutside && ! bPlane && ! bDangerInversion) { - - TRIA3DVECTOR vInnerTriaTemp ; - - for ( int ni = 0 ; ni < nContSize ; ++ ni) { - - // Se l'area è non nulla aggiungo il triangolo a uno dei due vettori. - if ( triContainer[ni].GetArea() > SQ_EPS_SMALL) { - int nVoxIJK[3] = { nVoxI, nVoxJ, nVoxK} ; - Point3d ptTemp = ptSol ; - Triangle3d trTemp = triContainer[ni] ; - vInnerTriaTemp.emplace_back( triContainer[ni]) ; - } - } - - if ( vInnerTriaTemp.size() > 0) { - - // Aggiorno il numero di feature. - ++ nInnerFeatureInVoxel ; - - // Se siamo in presenza della prima feature del - // voxel, ridimensiono il vettore che contiene - // la struttura dati del voxel. - if ( nInnerFeatureInVoxel == 1) { - triHold.resize( triHold.size() + 1) ; - // Questi dati dipendono solo dal voxel - int nCurrent = int( triHold.size()) - 1 ; - triHold[nCurrent].i = nVoxI ; - triHold[nCurrent].j = nVoxJ ; - triHold[nCurrent].k = nVoxK ; - } - - // Indice che identifica la posizione del voxel nel vector. - int nCurrent = int( triHold.size()) - 1 ; - - // Aggiungo vertice della componente connessa alla lista dei vertici. - triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; - - int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ; - int nNewFeatureNum = nOldFeatureNum + 1 ; - - // Aggiungo una componente per il vettore dei triangoli della componente connessa. - triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ; - for ( int ni = 0 ; ni < int( vInnerTriaTemp.size()) ; ++ ni) { - triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vInnerTriaTemp[ni]) ; - } - } - } - - // ExtMC non confermato, si passa a MC - else { - // Costruzione dei triangoli - for ( int TriIndex = 0; TriIndex < ( nVertComp[nComp] - 2) * 3 ; TriIndex += 3) { - // Il triangolo è pronto - Triangle3d CurrentTriangle ; - CurrentTriangle.Set( CompoTriVert[nComp][TriIndex].ptPApp, - CompoTriVert[nComp][TriIndex+1].ptPApp, - CompoTriVert[nComp][TriIndex+2].ptPApp) ; - bool bV = CurrentTriangle.Validate( true) ; - // Aggiungo alla lista - lstTria.emplace_back( CurrentTriangle) ; - } + } + } } + + size_t tOldCompo = CurrVox.nNumComp ; + ++ CurrVox.nNumComp ; + CurrVox.Compo[tOldCompo].nVertNum = nVertComp[nComp] ; + for ( int nV = 0 ; nV < nVertComp[nComp] ; ++ nV) { + CurrVox.Compo[tOldCompo].CompVecField[nV] = CompoVert[nComp][nV] ; + } + CurrVox.Compo[tOldCompo].ptVert = ptSol ; + CurrVox.Compo[tOldCompo].bInside = IsPointInsideVoxelApprox( nVoxI, nVoxJ, nVoxK, ptSol, EPS_SMALL) ; + CurrVox.Compo[tOldCompo].bCorner = ( nFeatureType == CORNER) ; } - - // Standard MC - else { - // Costruzione dei triangoli + // Standard MC + else { + vector vTria ; + // Costruzione dei triangoli for ( int TriIndex = 0; TriIndex < ( nVertComp[nComp] - 2) * 3 ; TriIndex += 3) { - // Il triangolo è pronto - Triangle3d CurrentTriangle ; + // Il triangolo è pronto + Triangle3dEx CurrentTriangle ; CurrentTriangle.Set( CompoTriVert[nComp][TriIndex].ptPApp, CompoTriVert[nComp][TriIndex+1].ptPApp, CompoTriVert[nComp][TriIndex+2].ptPApp) ; - bool bV = CurrentTriangle.Validate( true) ; - // Aggiungo alla lista - lstTria.emplace_back( CurrentTriangle) ; + // Setto il numero di utensile (conta solo positivo, nullo o negativo) + int nTool0 = Clamp( CompoTriVert[nComp][TriIndex].nPropIndex, -1, 1) ; + int nTool1 = Clamp( CompoTriVert[nComp][TriIndex+1].nPropIndex, -1, 1) ; + int nTool2 = Clamp( CompoTriVert[nComp][TriIndex+2].nPropIndex, -1, 1) ; + if ( nTool0 == nTool1 || nTool0 == nTool2) + CurrentTriangle.SetGrade( nTool0) ; + else if ( nTool1 == nTool2) + CurrentTriangle.SetGrade( nTool1) ; + // Valido il triangolo e setto le normali del campo vettoriale ai corrispondenti vertici + if ( CurrentTriangle.Validate( true)) { + for ( int nV = 0 ; nV < 3 ; ++ nV) { + const Vector3d& vtVertNorm = CompoTriVert[nComp][TriIndex+nV].vtVec ; + if ( CurrentTriangle.GetN() * vtVertNorm > 0.5) + CurrentTriangle.SetVertexNorm( nV, vtVertNorm) ; + else + CurrentTriangle.SetVertexNorm( nV, CurrentTriangle.GetN()) ; + } + } + // Riporto le coordinate nel sistema in cui è immerso lo Zmap + // CurrentTriangle.ToGlob( m_MapFrame) ; + // Aggiungo alla lista + vTria.emplace_back( CurrentTriangle) ; } - } + // Controllo i colori di configurazioni 2 e 8 + if ( ( nAllConfig[nIndex] == 2 || nAllConfig[nIndex] == 8) && + vTria[0].GetN() * vTria[1].GetN() > 0.8) { + if ( vTria[0].GetGrade() < 0 || vTria[1].GetGrade() < 0) { + vTria[0].SetGrade( -1) ; + vTria[1].SetGrade( -1) ; + } + else if ( vTria[0].GetGrade() > 0 || vTria[1].GetGrade() > 0) { + vTria[0].SetGrade( 1) ; + vTria[1].SetGrade( 1) ; + } + } + for ( int nT = 0 ; nT < int( vTria.size()) ; ++ nT) + lstTria.emplace_back( vTria[nT]) ; + } + } + // Se nel voxel abbiamo trovato feature + // aggiorniamo i contenitori. + if ( CurrVox.nNumComp > 0) { + CurrVox.i = nVoxI ; + CurrVox.j = nVoxJ ; + CurrVox.k = nVoxK ; + int nIJK[3] = { nVoxI, nVoxJ, nVoxK} ; + int nKey ; + GetVoxNFromIJK( nVoxI, nVoxJ, nVoxK, nKey) ; + VoxCont.emplace( nKey, CurrVox) ; } - } + } return true ; } @@ -1405,7 +1225,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, VoxelContainer& vV for ( int i = nLimits[0] ; i < nLimits[1] ; ++ i) { for ( int j = nLimits[2] ; j < nLimits[3] ; ++ j) { for ( int k = nLimits[4] ; k < nLimits[5] ; ++ k) { - + // Classificazione dei vertici: interni o esterni al materiale int nIndex = CalcIndex( i, j, k) ; @@ -2239,15 +2059,15 @@ VolZmap::CreateSharpFeatureTriangle( int nBlock, const VoxelContainer& vVoxel, T } } + // Pulisco il contenitore dei voxel di frontiera m_InterBlockTria[nBlock].clear() ; - // Ciclo sui voxel di frontiera for ( auto itVox = m_InterBlockVox[nBlock].begin() ; itVox != m_InterBlockVox[nBlock].end() ; ++ itVox) { // Indici del voxel int nVoxIJK[3] = { itVox->second.i, - itVox->second.j, - itVox->second.k} ; + itVox->second.j, + itVox->second.k} ; // Ridimensiono il contenitore dei triangoli interni int nOldSizeInn = int( triHold.size()) ; triHold.resize( nOldSizeInn + 1) ; @@ -2271,17 +2091,17 @@ VolZmap::CreateSharpFeatureTriangle( int nBlock, const VoxelContainer& vVoxel, T // Definisco il triangolo Triangle3dEx CurrTri ; CurrTri.Set( itVox->second.Compo[nComp].ptVert, - itVox->second.Compo[nComp].CompVecField[nNextVert].ptPApp, - itVox->second.Compo[nComp].CompVecField[nVert].ptPApp) ; + itVox->second.Compo[nComp].CompVecField[nNextVert].ptPApp, + itVox->second.Compo[nComp].CompVecField[nVert].ptPApp) ; // Setto il numero di utensile ai vertici di base del fan CurrTri.SetAttrib( 1, itVox->second.Compo[nComp].CompVecField[nNextVert].nPropIndex) ; CurrTri.SetAttrib( 2, itVox->second.Compo[nComp].CompVecField[nVert].nPropIndex) ; // Setto il numero di utensile al triangolo nel complesso if ( CurrTri.GetAttrib( 1) < 0 || - CurrTri.GetAttrib( 2) < 0) + CurrTri.GetAttrib( 2) < 0) CurrTri.SetGrade( - 1) ; else if ( CurrTri.GetAttrib( 1) > 0 || - CurrTri.GetAttrib( 2) > 0) + CurrTri.GetAttrib( 2) > 0) CurrTri.SetGrade( 1) ; else CurrTri.SetGrade( 0) ; @@ -2290,12 +2110,12 @@ VolZmap::CreateSharpFeatureTriangle( int nBlock, const VoxelContainer& vVoxel, T CurrTri.SetVertexNorm( 2, itVox->second.Compo[nComp].CompVecField[nVert].vtVec) ; if ( CurrTri.GetVertexNorm( 1) * CurrTri.GetVertexNorm( 2) > 0.5) CurrTri.SetVertexNorm( 0, 0.5 * ( CurrTri.GetVertexNorm( 1) + - CurrTri.GetVertexNorm( 2))) ; + CurrTri.GetVertexNorm( 2))) ; // Valido il triangolo CurrTri.Validate( true) ; // Smisto i triangoli fra di frontiera e interni bool bTriOnBorder = IsATriangleOnBorder( CurrTri, itVox->second.Compo[nComp].ptVert, - nLimits, nVoxIJK) ; + nLimits, nVoxIJK) ; // Triangolo di frontiera if ( bTriOnBorder) { if ( bNewCompBor) { @@ -2335,10 +2155,68 @@ VolZmap::CreateSharpFeatureTriangle( int nBlock, const VoxelContainer& vVoxel, T } } } - + return true ; } +//---------------------------------------------------------------------------- +bool +VolZmap::CreateSharpFeatureTriangle( const VoxelContainer& vVoxel, TriHolder& triHold) const +{ + // Ciclo sui voxel interni + for ( auto it = vVoxel.begin() ; it != vVoxel.end() ; ++ it) { + size_t tOldSize = triHold.size() ; + triHold.resize( tOldSize + 1) ; + triHold[tOldSize].i = it->second.i ; + triHold[tOldSize].j = it->second.j ; + triHold[tOldSize].k = it->second.k ; + // Ciclo sulle componenti connesse del voxel + for ( int nComp = 0 ; nComp < it->second.nNumComp ; ++ nComp) { + triHold[tOldSize].ptCompoVert.emplace_back( it->second.Compo[nComp].ptVert) ; + // triHold[tOldSize].ptCompoVert.back().ToGlob( m_MapFrame) ; + size_t tOldCompNum = triHold[tOldSize].vCompoTria.size() ; + triHold[tOldSize].vCompoTria.resize( tOldCompNum + 1) ; + triHold[tOldSize].vbFlipped.resize( tOldCompNum + 1) ; + // ciclo sui vertici della componente connessa + int nNumVert = it->second.Compo[nComp].nVertNum ; + for ( int nVert = 0 ; nVert < nNumVert ; ++ nVert) { + int nNextVert = ( nVert + 1 < nNumVert ? nVert + 1 : 0) ; + // Definisco il triangolo + Triangle3dEx CurrTri ; + CurrTri.Set( it->second.Compo[nComp].ptVert, + it->second.Compo[nComp].CompVecField[nNextVert].ptPApp, + it->second.Compo[nComp].CompVecField[nVert].ptPApp) ; + // Setto il numero di utensile ai vertici di base del fan + CurrTri.SetAttrib( 1, it->second.Compo[nComp].CompVecField[nNextVert].nPropIndex) ; + CurrTri.SetAttrib( 2, it->second.Compo[nComp].CompVecField[nVert].nPropIndex) ; + // Setto il numero di utensile al triangolo nel complesso + if ( CurrTri.GetAttrib( 1) < 0 || + CurrTri.GetAttrib( 2) < 0) + CurrTri.SetGrade( - 1) ; + else if ( CurrTri.GetAttrib( 1) > 0 || + CurrTri.GetAttrib( 2) > 0) + CurrTri.SetGrade( 1) ; + else + CurrTri.SetGrade( 0) ; + // Setto le normali a ogni vertice + CurrTri.SetVertexNorm( 1, it->second.Compo[nComp].CompVecField[nNextVert].vtVec) ; + CurrTri.SetVertexNorm( 2, it->second.Compo[nComp].CompVecField[nVert].vtVec) ; + if ( CurrTri.GetVertexNorm( 1) * CurrTri.GetVertexNorm( 2) > 0.5) + CurrTri.SetVertexNorm( 0, 0.5 * ( CurrTri.GetVertexNorm( 1) + + CurrTri.GetVertexNorm( 2))) ; + // Valido il triangolo + CurrTri.Validate( true) ; + triHold[tOldSize].vCompoTria[tOldCompNum].emplace_back( CurrTri) ; + //triHold[tOldSize].vCompoTria[tOldCompNum].back().ToGlob( m_MapFrame) ; + size_t tTri = triHold[tOldSize].vbFlipped[tOldCompNum].size() ; + triHold[tOldSize].vbFlipped[tOldCompNum].resize( tTri + 1) ; + triHold[tOldSize].vbFlipped[tOldCompNum][tTri] = false ; + } + } + } + + return true ; +} //---------------------------------------------------------------------------- // Calcola i triangoli dei voxel contenuti nel vettore, se si vuole il riconoscimento @@ -2347,8 +2225,8 @@ VolZmap::CreateSharpFeatureTriangle( int nBlock, const VoxelContainer& vVoxel, T bool VolZmap::ExtMarchingCubes( vector& vVox, TRIA3DEXLIST& lstTria, bool bEnh) const { - // Contenitore dei triangoli costituenti sharp-feaure - TriHolder triHold ; + // + VoxelContainer vVoxCont ; int nOriginalSize = int( vVox.size()) ; for ( int nV = 0 ; nV < int( vVox.size()) ; ++ nV) { // I voxel intersecati dalla retta hanno indice nel vettore @@ -2357,9 +2235,9 @@ VolZmap::ExtMarchingCubes( vector& vVox, TRIA3DEXLIST& lstTria, bo if ( nV < nOriginalSize) { // Se a questa iterazione troviamo una sharp-feature // la dimensione del TriHolder aumenta - int nCurrentSize = int( triHold.size()) ; - ProcessCube( vVox[nV].nI, vVox[nV].nJ, vVox[nV].nK, lstTria, triHold, bEnh) ; - int nNewSize = int( triHold.size()) ; + int nCurrentSize = int( vVoxCont.size()) ; + ProcessCube( vVox[nV].nI, vVox[nV].nJ, vVox[nV].nK, lstTria, vVoxCont, bEnh) ; + int nNewSize = int( vVoxCont.size()) ; // Se abbiamo trovato una sharp-feature if ( nNewSize > nCurrentSize) { @@ -2395,8 +2273,12 @@ VolZmap::ExtMarchingCubes( vector& vVox, TRIA3DEXLIST& lstTria, bo // Se il voxel non è attraversato dalla retta, ma fa parte di // quelli adiacenti lo processiamo normalmente else - ProcessCube( vVox[nV].nI, vVox[nV].nJ, vVox[nV].nK, lstTria, triHold, bEnh) ; + ProcessCube( vVox[nV].nI, vVox[nV].nJ, vVox[nV].nK, lstTria, vVoxCont, bEnh) ; } + + // Contenitore dei triangoli costituenti sharp-feaure + TriHolder triHold ; + CreateSharpFeatureTriangle( vVoxCont, triHold) ; // Eseguo il flipping FlipEdgesII( triHold) ; // Aggiungo alla lista di triangoli, tutti quelli che formano una sharp-feature. @@ -2447,7 +2329,9 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const Point3d ptP22 = TriHold[n2].vCompoTria[nCompo2][nTri2].GetP( 2) ; // I triangoli sono da flippare if ( AreSamePointEpsilon( ptP11, ptP22, EPS_ZERO) && - AreSamePointEpsilon( ptP12, ptP21, EPS_ZERO)) { + AreSamePointEpsilon( ptP12, ptP21, EPS_ZERO) && + ! ( TriHold[n1].vbFlipped[nCompo1][nTri1] || + TriHold[n2].vbFlipped[nCompo2][nTri2])) { // Assegno l'array dei punti di contorno Point3d vPnt[4] ; vPnt[0] = TriHold[n1].vCompoTria[nCompo1][nTri1].GetP( 1) ; @@ -2498,6 +2382,9 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const TriHold[n1].vCompoTria[nCompo1][nTri1].SetVertexNorm( 1, vtN1) ; TriHold[n2].vCompoTria[nCompo2][nTri2].SetVertexNorm( 0, vtN2) ; TriHold[n2].vCompoTria[nCompo2][nTri2].SetVertexNorm( 1, vtN2) ; + // Setto i triangoli come flippati + TriHold[n1].vbFlipped[nCompo1][nTri1] = true ; + TriHold[n2].vbFlipped[nCompo2][nTri2] = true ; // Valido i triangoli TriHold[n1].vCompoTria[nCompo1][nTri1].Validate( true) ; TriHold[n2].vCompoTria[nCompo2][nTri2].Validate( true) ; @@ -2609,7 +2496,9 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const Point3d ptPL2 = InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 2) ; // Si deve operare la modifica dei triangoli if ( AreSamePointEpsilon( ptPF1, ptPL2, EPS_ZERO) && - AreSamePointEpsilon( ptPF2, ptPL1, EPS_ZERO)) { + AreSamePointEpsilon( ptPF2, ptPL1, EPS_ZERO) && + ! ( InterTria[tFB][tVFB].vbFlipped[tCmpF][tTriFB] || + InterTria[tLB][tVLB].vbFlipped[tCmpL][tTriLB])) { // Assegno l'array dei punti di contorno Point3d vPnt[4] ; vPnt[0] = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 1) ; @@ -2663,6 +2552,9 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetVertexNorm( 0, vtNormF) ; InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetVertexNorm( 1, vtNormL) ; InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetVertexNorm( 0, vtNormL) ; + // Setto i triangoli come flippati + InterTria[tFB][tVFB].vbFlipped[tCmpF][tTriFB] = true ; + InterTria[tLB][tVLB].vbFlipped[tCmpL][tTriLB] = true ; bModified = true ; break ; }