diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index f55c585..75c3f2d 100644 Binary files a/EgtGeomKernel.rc and b/EgtGeomKernel.rc differ diff --git a/VolTriZmapCalculus.cpp b/VolTriZmapCalculus.cpp index 3c25fff..2715de0 100644 --- a/VolTriZmapCalculus.cpp +++ b/VolTriZmapCalculus.cpp @@ -480,7 +480,7 @@ VolZmap::AvoidBox( const Frame3d& frBox, const Vector3d& vtDiag) for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) { if ( ! ( dZmax < m_Values[0][nPos][nIndex].dMin - EPS_SMALL || - dZmin > m_Values[0][nPos][nIndex + 1].dMax + EPS_SMALL)) + dZmin > m_Values[0][nPos][nIndex].dMax + EPS_SMALL)) return false ; } } diff --git a/VolTriZmapGraphics.cpp b/VolTriZmapGraphics.cpp index 8e2acf9..268ba29 100644 --- a/VolTriZmapGraphics.cpp +++ b/VolTriZmapGraphics.cpp @@ -713,7 +713,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) 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) { - + // Riconoscimento dei voxel di frontiera int nVoxIndexes[3] = { i, j, k} ; bool bBoundary = IsAVoxelOnBoundary( nLimits, nVoxIndexes, true) ; @@ -851,169 +851,173 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) nStdTabOff += 3 * ( nVertComp[nCompCount - 1] - 2) ; } - // Test sulla topologia - if ( nAllConfig[nIndex] == 3) { + // 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. + if ( bReg) { + if ( nAllConfig[nIndex] == 3) { - Vector3d vtCmpAvg0, vtCmpAvg1 ; + Vector3d vtCmpAvg0, vtCmpAvg1 ; - // 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). - bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ; - bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ; + // 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). + 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 ; + // 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 ; + if ( bTest0 && bTest1) + dScProd = vtCmpAvg0 * vtCmpAvg1 ; - double dThreshold = 0.7 ; + double dThreshold = 0.7 ; - if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) { + if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) { - int nt = 0 ; + int nt = 0 ; - while ( nIndexVsIndex3[nt][0] != nIndex) - ++ nt ; + while ( nIndexVsIndex3[nt][0] != nIndex) + ++ nt ; - int nRotCase = nIndexVsIndex3[nt][1] ; + int nRotCase = nIndexVsIndex3[nt][1] ; - nComponents = Cases3Plus[nRotCase][1][0] ; + nComponents = Cases3Plus[nRotCase][1][0] ; - // Riaggiorno gli offsets - nExtTabOff = nComponents ; - nStdTabOff = 0 ; + // Riaggiorno gli offsets + nExtTabOff = nComponents ; + nStdTabOff = 0 ; - // Modifico le matrici - for ( int nC = 1 ; nC <= nComponents ; ++ nC) { + // Modifico le matrici + for ( int nC = 1 ; nC <= nComponents ; ++ nC) { - // Numero vertici per componenti - nVertComp[nC - 1] = Cases3Plus[nRotCase][1][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) + // 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]] ; + CompoVert[nC - 1][nFanVert] = VecField[Cases3Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ; - // Matrici dei vertici dei triangoli in assenza di sharp feature - for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) { + // 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]] ; - } + 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]] ; + } - // Aggiorno gli offsets per raggiungere i - // vertici della componente successiva. - nExtTabOff += nVertComp[nC - 1] ; - nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; - } + // Aggiorno gli offsets per raggiungere i + // vertici della componente successiva. + nExtTabOff += nVertComp[nC - 1] ; + nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; + } - } - } - 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 ; - - if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) { - - int nt = 0 ; - - while ( nIndexVsIndex6[nt][0] != nIndex) - ++ nt ; - - int nRotCase = nIndexVsIndex6[nt][1] ; - - - // Costruzione dei triangoli - for ( int TriIndex = 0 ; TriIndex < 15 ; TriIndex += 3) { - - // Costruzione triangolo - int i0 = Cases6Plus[nRotCase][TriIndex + 2] ; - int i1 = Cases6Plus[nRotCase][TriIndex + 1] ; - int i2 = Cases6Plus[nRotCase][TriIndex] ; - - Triangle3d CurrentTriangle ; - - // Il triangolo è pronto - CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ; - bool bV = CurrentTriangle.Validate( true) ; - - // Aggiungo alla lista - lstTria.emplace_back( CurrentTriangle) ; - } - continue ; + } } - } - else if ( nAllConfig[nIndex] == 10) { + else if ( nAllConfig[nIndex] == 6) { + + // Procedura analoga a quella della configurazione 3 + Vector3d vtCmpAvg0, vtCmpAvg1 ; - Vector3d vtCmpAvg0, vtCmpAvg1 ; + bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.7) ; + bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ; + + double dScProd = - 2 ; - // 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 - bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ; - bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ; + if ( bTest0 && bTest1) + dScProd = vtCmpAvg0 * vtCmpAvg1 ; - if ( ! ( bTest0 && bTest1)) { + double dThreshold = 0.7 ; - int nt = 0 ; + if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) { + + int nt = 0 ; + + while ( nIndexVsIndex6[nt][0] != nIndex) + ++ nt ; + + int nRotCase = nIndexVsIndex6[nt][1] ; - while ( nIndexVsIndex10[nt][0] != nIndex) - ++ nt ; - int nRotCase = nIndexVsIndex10[nt][1] ; + // Costruzione dei triangoli + for ( int TriIndex = 0 ; TriIndex < 15 ; TriIndex += 3) { - // Riaggiorno gli offsets - nExtTabOff = 2 ; - nStdTabOff = 0 ; + // Costruzione triangolo + int i0 = Cases6Plus[nRotCase][TriIndex + 2] ; + int i1 = Cases6Plus[nRotCase][TriIndex + 1] ; + int i2 = Cases6Plus[nRotCase][TriIndex] ; - // Modifico le matrici - for ( int nC = 1 ; nC <= 2 ; ++ nC) { - - // Numero vertici per componenti - nVertComp[nC - 1] = Cases10Plus[nRotCase][1][nC] ; + Triangle3d CurrentTriangle ; - // Matrice dei vertici della base del fan - for ( int nFanVert = 0 ; nFanVert < 4 ; ++ nFanVert) + // Il triangolo è pronto + CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ; + bool bV = CurrentTriangle.Validate( true) ; - 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]] ; - CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ; - CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert]] ; + // Aggiungo alla lista + lstTria.emplace_back( CurrentTriangle) ; } + continue ; + } + } + else if ( nAllConfig[nIndex] == 10) { + + Vector3d vtCmpAvg0, vtCmpAvg1 ; + + // 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 + bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ; + bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ; + + if ( ! ( bTest0 && bTest1)) { + + int nt = 0 ; + + while ( nIndexVsIndex10[nt][0] != nIndex) + ++ nt ; - // Aggiorno gli offsets per raggiungere i - // vertici della componente successiva. - nExtTabOff += nVertComp[nC - 1] ; - nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; - } - } + int nRotCase = nIndexVsIndex10[nt][1] ; + + // Riaggiorno gli offsets + nExtTabOff = 2 ; + nStdTabOff = 0 ; + + // Modifico le matrici + 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]] ; + + // 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]] ; + CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ; + CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert]] ; + } + + // Aggiorno gli offsets per raggiungere i + // vertici della componente successiva. + nExtTabOff += nVertComp[nC - 1] ; + nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; + } + } + } } @@ -1288,7 +1292,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) int nPosD ; Vector3d vtD, vtE ; - if ( nVertComp[nCompCount - 1] == 4 && nFeatureType == 2) { + if ( nVertComp[nCompCount - 1] == 4 && nFeatureType == EDGE) { int nPosEq ; for ( int ni = 0 ; ni < 2 ; ++ ni) { for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) { @@ -1410,17 +1414,16 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) dUnknownVector( 2)) ; double dDistFeature = vtFeature.Len() ; - const double dMaxDist = sqrt( 3) * m_dStep ; + const double MAX_DIST = sqrt( 3) * m_dStep ; // Limito la feature a una distanza dal // baricentro pari alla diagonale del voxel. - if ( dDistFeature > dMaxDist) + if ( dDistFeature > MAX_DIST) bOutside = true ; Triangle3d CurrentTriangle ; TRIA3DVECTOR triContainer ; - - + // Costruisco triangoli di prova for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { @@ -1432,106 +1435,242 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // Aggiungo triangolo al vettore temporaneo triContainer.emplace_back( CurrentTriangle) ; } - - // Controllo delle inversioni dei triangoli - bool bInversione = false ; + + // Controllo delle inversioni dei triangoli bool bDangerInversion = false ; + + // Caso ventaglio con tre vertici di base + if ( nVertComp[nCompCount - 1] == 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[nCompCount - 1] ; ++ ni) { - for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; - int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; - - if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - EPS_SMALL || - triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL) + double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; + double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ; + + if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) bInversione = true ; - } - - if ( bInversione && - ( nVertComp[nCompCount - 1] == 3 || - nVertComp[nCompCount - 1] == 4)) { - - int nNegDot = 0 ; - - for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] - 1 ; ++ ni) { - for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) { - if ( CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL) - nNegDot ++ ; - } } - if ( nNegDot == nVertComp[nCompCount - 1] - 1) { + + // Se tale triangolo esiste procedo + if ( bInversione) { - Point3d ptSolZMapFrame = ptSol ; - ptSolZMapFrame.ToLoc( m_MapFrame[0]) ; - - - if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) { - - Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; - vtNullSpace.ToLoc( m_MapFrame[0]) ; - double dParInt1, dParInt2 ; - Point3d ptVoxMin( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, ( k + 0.5) * m_dStep) ; - Point3d ptVoxMax( ( i + 1.5) * m_dStep, ( j + 1.5) * m_dStep, ( k + 1.5) * m_dStep) ; - - 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 ; - - ptNewSol.ToGlob( m_MapFrame[0]) ; - - ptSol = ptNewSol ; - - // Costruisco triangoli di prova - for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { - - int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; - // Il triangolo è pronto - CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; - CurrentTriangle.Validate( true) ; - - // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; - } + // Conto le coppie di normali con angolo compreso maggiore di 90° + int nNegDot = 0 ; + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] - 1 ; ++ ni) { + for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) { + if ( CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL) + nNegDot ++ ; } - else { + } + + if ( nNegDot == nVertComp[nCompCount - 1] - 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[nCompCount - 1] ; ++ ni) { + + int nj = ( ni == 0 ? nVertComp[nCompCount - 1] - 1 : ni - 1) ; + + double dDLast = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; + double dDPrev = triContainer[nj].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; + + if ( dDLast < EPS_SMALL && dDPrev < EPS_SMALL) + bInversione2 = true ; + } + + // Se tale normale esiste + if ( bInversione2) { + + // Soluzione del sistema nel sistema Zmap + Point3d ptSolZMapFrame = ptSol ; + ptSolZMapFrame.ToLoc( m_MapFrame[0]) ; + + // Se la soluzione non cade nel voxel di appartenenza vedo se può + // essere riportata dentro muovendosi lungo la linea di feature. + if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) { + + Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; + vtNullSpace.ToLoc( m_MapFrame[0]) ; + double dParInt1, dParInt2 ; + Point3d ptVoxMin( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, ( k + 0.5) * m_dStep) ; + Point3d ptVoxMax( ( i + 1.5) * m_dStep, ( j + 1.5) * m_dStep, ( k + 1.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 ; + + ptNewSol.ToGlob( m_MapFrame[0]) ; + + ptSol = ptNewSol ; + + // Costruisco triangoli di prova + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { + + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + // Il triangolo è pronto + CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; + 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)) { + int nAdjVoxI, nAdjVoxJ, nAdjVoxK ; + if ( GetPointVoxel( ptSolZMapFrame, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) { - // Classificazione del voxel adiacente - int nAdjIndex = 0 ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK)) - nAdjIndex |= ( 1 << 0) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK)) - nAdjIndex |= ( 1 << 1) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK)) - nAdjIndex |= ( 1 << 2) ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK)) - nAdjIndex |= ( 1 << 3) ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 4) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 5) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 6) ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 7) ; + // Classificazione del voxel adiacente + int nAdjIndex = 0 ; + if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK)) + nAdjIndex |= ( 1 << 0) ; + if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK)) + nAdjIndex |= ( 1 << 1) ; + if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK)) + nAdjIndex |= ( 1 << 2) ; + if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK)) + nAdjIndex |= ( 1 << 3) ; + if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK + 1)) + nAdjIndex |= ( 1 << 4) ; + if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK + 1)) + nAdjIndex |= ( 1 << 5) ; + if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK + 1)) + nAdjIndex |= ( 1 << 6) ; + if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK + 1)) + nAdjIndex |= ( 1 << 7) ; - // Se il voxel è pieno - if ( EdgeTable[nAdjIndex] != 0) - bDangerInversion = true ; - } - } - } - } - else { + // Se il voxel è pieno + if ( EdgeTable[nAdjIndex] != 0) + bDangerInversion = true ; + } + } + } + } + } + } + } + // Ventaglio con base a quattro vertici + else if ( nVertComp[nCompCount - 1] == 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[nCompCount - 1] ; ++ ni) { - if ( nVertComp[nCompCount - 1] == 4) { + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + + double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; + double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ; + + if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) + bInversione = true ; + } + // Se tale triangolo esiste continuo i test + if ( bInversione) { + + // Conto il numero di coppie di normali con prodotto scalare negativo + int nNegDot = 0 ; + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] - 1 ; ++ ni) { + for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) { + double dDot = CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm ; + if ( dDot < - EPS_SMALL) + nNegDot ++ ; + } + } + + // Caso in cui tale numero è 3 + if ( nNegDot == nVertComp[nCompCount - 1] - 1) { + + Point3d ptSolZMapFrame = ptSol ; + ptSolZMapFrame.ToLoc( m_MapFrame[0]) ; + + // Se la feature non cade nel suo voxel + if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) { + + Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; + vtNullSpace.ToLoc( m_MapFrame[0]) ; + double dParInt1, dParInt2 ; + Point3d ptVoxMin( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, ( k + 0.5) * m_dStep) ; + Point3d ptVoxMax( ( i + 1.5) * m_dStep, ( j + 1.5) * m_dStep, ( k + 1.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( 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 ; + + ptNewSol.ToGlob( m_MapFrame[0]) ; + + ptSol = ptNewSol ; + + // Costruisco triangoli di prova + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { + + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + // Il triangolo è pronto + CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; + 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( ptSolZMapFrame, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) { + + // Classificazione del voxel adiacente + int nAdjIndex = 0 ; + if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK)) + nAdjIndex |= ( 1 << 0) ; + if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK)) + nAdjIndex |= ( 1 << 1) ; + if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK)) + nAdjIndex |= ( 1 << 2) ; + if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK)) + nAdjIndex |= ( 1 << 3) ; + if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK + 1)) + nAdjIndex |= ( 1 << 4) ; + if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK + 1)) + nAdjIndex |= ( 1 << 5) ; + if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK + 1)) + nAdjIndex |= ( 1 << 6) ; + if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK + 1)) + nAdjIndex |= ( 1 << 7) ; + + // Se il voxel è pieno + if ( EdgeTable[nAdjIndex] != 0) + bDangerInversion = true ; + } + } + } + } + // Caso in cui il numero di coppie di normali con prodotto + // scalare negativo non è 3 + else { Point3d ptSolZMapFrame = ptSol ; ptSolZMapFrame.ToLoc( m_MapFrame[0]) ; @@ -1571,10 +1710,10 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) } } } - } - } + } + } } - + // Valuto normali: questo è ancora un controllo // sulle normali, se risultano in tutti i punti // approssimativamente uguali passiamo alla @@ -1598,7 +1737,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) if ( ! bPlane) break ; } - + // Se la feature non è fuori dai limiti // e i triangoli formano una superficie // non piana confermo ExtMC @@ -1717,7 +1856,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) } // ExtMC non confermato, si passa a MC else { - + // Costruzione dei triangoli for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) { @@ -1736,7 +1875,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) } // Standard MC else { - + // Costruzione dei triangoli for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) { @@ -2796,5 +2935,3 @@ VolZmap::IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert, else return false ; } - - diff --git a/VolTriZmapVolume.cpp b/VolTriZmapVolume.cpp index 8c2cf4a..ce7c9c4 100644 --- a/VolTriZmapVolume.cpp +++ b/VolTriZmapVolume.cpp @@ -3272,7 +3272,7 @@ VolZmap::CompCyl_ZDrilling( unsigned int nGrid, const Point3d & ptS, const Point // Parametri geometrici dell'utensile double dSqRad = dRad * dRad ; - double dSqfeSqRad = dSqRad - 2 * dRad * EPS_SMALL ; + double dSafeSqRad = dSqRad - 2 * dRad * EPS_SMALL ; // Punte del gambo Point3d ptTStemS = ptS - vtToolDir * dHei ; @@ -3293,7 +3293,7 @@ VolZmap::CompCyl_ZDrilling( unsigned int nGrid, const Point3d & ptS, const Point double dSqLen = vtC.SqLen() ; // Se il punto si trova dentro il cerchio taglio - if ( dSqLen < dSqfeSqRad) + if ( dSqLen < dSafeSqRad) SubtractIntervals( nGrid, i, j, dMinStemZ, dMaxStemZ, Z_AX, - Z_AX) ; } } @@ -3318,8 +3318,9 @@ VolZmap::CompConus_ZDrilling( unsigned int nGrid, const Point3d & ptS, const Poi double dAngC = dHei / ( dMaxRad - dMinRad) ; double dSqMinRad = dMinRad * dMinRad ; double dSqMaxRad = dMaxRad * dMaxRad ; - double dDeltaR = dMaxRad - dMinRad ; - + double dSafeSqMaxRad = dSqMaxRad - 2 * dMaxRad * EPS_SMALL ; // Questa variabile è sperimentale: serve per evitare il taglio di un dexel dalla parte cilindrica del volume spazzato dalla traslazione del cono. + double dDeltaR = dMaxRad - dMinRad ; // Per tornare alla versione precedente basta sostituire dSafeSqMaxRad con dSqMaxRad. Per risolvere il problema in modo forse più sicuro, ma + // computazionalmente più pesante è sottrarre prima il cilindro con dSafeSqMaxRad e dopo il cono con dSqMaxRad. // Studio delle simmetrie if ( vtToolDir.z > 0) { @@ -3354,7 +3355,7 @@ VolZmap::CompConus_ZDrilling( unsigned int nGrid, const Point3d & ptS, const Poi SubtractIntervals( nGrid, i, j, dZMin, dZMax, Z_AX, -Z_AX) ; - else if ( dSqDist < dSqMaxRad) { + else if ( dSqDist < dSafeSqMaxRad) { // dSafeSqMaxRad è sperimentale double dr = sqrt( dSqDist) ;