diff --git a/VolTriZmapGraphics.cpp b/VolTriZmapGraphics.cpp index f77fb23..b535ddc 100644 --- a/VolTriZmapGraphics.cpp +++ b/VolTriZmapGraphics.cpp @@ -2598,224 +2598,6 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) } } - // 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) { - - 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 ; - - // 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)) ; - - double dParInt1, dParInt2 ; - Point3d ptVoxMin( ( i * N_DEXVOXRATIO + 0.5) * m_dStep, - ( j * N_DEXVOXRATIO + 0.5) * m_dStep, ( k * N_DEXVOXRATIO + 0.5) * m_dStep) ; - Point3d ptVoxMax( ( ( i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, - ( ( j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( k + 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( 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[nCompCount - 1] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; - // Il triangolo è pronto - Triangle3d CurrentTriangle ; - 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 = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ; - // 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 ; - - if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) { - - Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; - - double dParInt1, dParInt2 ; - Point3d ptVoxMin( ( i * N_DEXVOXRATIO + 0.5) * m_dStep, - ( j * N_DEXVOXRATIO + 0.5) * m_dStep, ( k * N_DEXVOXRATIO + 0.5) * m_dStep) ; - Point3d ptVoxMax( ( ( i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, - ( ( j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( k + 1) * N_DEXVOXRATIO + 0.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 ; - - 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 - Triangle3d CurrentTriangle ; - 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) ; - } - } - - else { - int nCouple = 0 ; - int nCoupleIndex[4] ; - - // Valuto il numero di coppie di vettori - // quasi coincidenti e per ogni coppia salvo gli - // indici dei vettori - for ( int nNI = 0 ; nNI < 3 ; ++ nNI) { - for ( int nNJ = nNI + 1 ; nNJ < 4 ; ++ nNJ) { - if ( AreSameVectorApprox( CompoVert[nCompCount - 1][nNI].vtNorm, - CompoVert[nCompCount - 1][nNJ].vtNorm)) { - ++ nCouple ; - if ( nCouple == 1) { - nCoupleIndex[0] = nNI ; - nCoupleIndex[1] = nNJ ; - } - else if ( nCouple == 2) { - nCoupleIndex[2] = nNI ; - nCoupleIndex[3] = nNJ ; - } - } - } - } - - // caso due coppie - if ( nCouple == 2) { - // vedo se c'è un triangolo con normale invertita rispetto a quelle - // si entrambi i vertici, se esiste si passerà a std MC - for ( int ni = 0 ; ni < 4 ; ++ ni) { - int nj = ( ni == 3 ? 0 : ni + 1) ; - if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - 0.95 && - triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - 0.95) { - bDangerInversion = true ; - break ; - } - } - } - // caso una coppia - else if ( nCouple == 1) { - // cerco gli indici dei vettori non appartenenti alla coppia - for ( int ni = 0 ; ni < 4 ; ++ ni) { - if ( ni != nCoupleIndex[0] && ni != nCoupleIndex[1]) { - nCoupleIndex[2] = ni ; - break ; - } - } - for ( int ni = 0 ; ni < 4 ; ++ ni) { - if ( ni != nCoupleIndex[0] && ni != nCoupleIndex[1] && ni != nCoupleIndex[2]) { - nCoupleIndex[3] = ni ; - break ; - } - } - // Media dei vettori coppia - Vector3d vtAv01 = 0.5 * ( CompoVert[nCompCount - 1][nCoupleIndex[0]].vtNorm + - CompoVert[nCompCount - 1][nCoupleIndex[1]].vtNorm) ; - // vettore nello spazio genenrato dai due non appartenenti alla coppia - Vector3d vtAv23 = 0.5 * ( CompoVert[nCompCount - 1][nCoupleIndex[2]].vtNorm + - CompoVert[nCompCount - 1][nCoupleIndex[3]].vtNorm) ; - - vtAv01.Normalize() ; - vtAv23.Normalize() ; - double dDAvAV = vtAv01 * vtAv23 ; - // se angolo grande si esegue std MC - if ( abs( vtAv01 * vtAv23) < EPS_SMALL) { - for ( int ni = 0 ; ni < 4 ; ++ ni) { - int nj = ni == 3 ? 0 : ni + 1 ; - if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - 0.95 && - triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - 0.95) { - bDangerInversion = true ; - break ; - } - } - } - else { - - double dD23 = CompoVert[nCompCount - 1][nCoupleIndex[2]].vtNorm * - CompoVert[nCompCount - 1][nCoupleIndex[3]].vtNorm ; - if ( dD23 > 0.7 && dDAvAV < 0.7) { - - for ( int ni = 0 ; ni < 4 ; ++ ni) { - int nj = ni == 3 ? 0 : ni + 1 ; - if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - 0.9 && - triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - 0.9) { - bDangerInversion = true ; - break ; - } - } - } - } - } - } - } - } - } - } - // Controllo sulle normali, se risultano in tutti i punti // approssimativamente uguali passiamo alla routine standard int nContSize = int( triContainer.size()) ;