diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index 666d9cc..d7d74f8 100644 Binary files a/EgtGeomKernel.rc and b/EgtGeomKernel.rc differ diff --git a/VolZmapGraphics.cpp b/VolZmapGraphics.cpp index 21808c7..bf9cd8b 100644 --- a/VolZmapGraphics.cpp +++ b/VolZmapGraphics.cpp @@ -1648,7 +1648,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& 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) { - + // Classificazione dei vertici: interni o esterni al materiale int nIndex = CalcIndex( i, j, k) ; @@ -1857,21 +1857,39 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold else if ( nAllConfig[nIndex] == 6) { // Procedura analoga a quella della configurazione 3 - Vector3d vtCmpAvg0, vtCmpAvg1 ; - bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.7) ; + /*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 ; + bool bFirstTest = ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > 0.7) ;*/ + // Altro metodo + double dDotSum = 0 ; + for ( int nFV = 0 ; nFV < 4 ; ++ nFV) { + for ( int nTV = 0 ; nTV < 3 ; ++ nTV) { + dDotSum += CompoVert[0][nFV].vtNorm * CompoVert[1][nTV].vtNorm ; + } + } + for ( int nFVI = 0 ; nFVI < 3 ; ++ nFVI) { + for ( int nFVJ = nFVI + 1 ; nFVJ < 4 ; ++ nFVJ) { + dDotSum -= CompoVert[0][nFVI].vtNorm * CompoVert[0][nFVJ].vtNorm ; + } + } + for ( int nTVI = 0 ; nTVI < 2 ; ++ nTVI) { + for ( int nTVJ = nTVI + 1 ; nTVJ < 3 ; ++ nTVJ) { + dDotSum -= CompoVert[1][nTVI].vtNorm * CompoVert[1][nTVJ].vtNorm ; + } + } + bool bTestOnSum = dDotSum > - 4 ; - if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > 0.7)) { + if ( bTestOnSum) { 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) { // Indici vertici @@ -1885,10 +1903,37 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold CurrentTriangle.SetVertexNorm( 0, VecField[i0].vtNorm) ; CurrentTriangle.SetVertexNorm( 1, VecField[i1].vtNorm) ; CurrentTriangle.SetVertexNorm( 2, VecField[i2].vtNorm) ; + // Setto il numero di utensile (conta solo positivo, nullo o negativo) + int nTool0 = Clamp( VecField[i0].nToolFlag, -1, 1) ; + int nTool1 = Clamp( VecField[i1].nToolFlag, -1, 1) ; + int nTool2 = Clamp( VecField[i2].nToolFlag, -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].vtNorm ; + if ( CurrentTriangle.GetN() * vtVertNorm > 0.6) + CurrentTriangle.SetVertexNorm( 0, vtVertNorm) ; + else + CurrentTriangle.SetVertexNorm( 0, CurrentTriangle.GetN()) ; + vtVertNorm = VecField[i1].vtNorm ; + if ( CurrentTriangle.GetN() * vtVertNorm > 0.6) + CurrentTriangle.SetVertexNorm( 1, vtVertNorm) ; + else + CurrentTriangle.SetVertexNorm( 1, CurrentTriangle.GetN()) ; + vtVertNorm = VecField[i2].vtNorm ; + if ( CurrentTriangle.GetN() * vtVertNorm > 0.6) + CurrentTriangle.SetVertexNorm( 2, vtVertNorm) ; + else + CurrentTriangle.SetVertexNorm( 2, CurrentTriangle.GetN()) ; + } + // Riporto il triangolo nel sistema in cui è immerso quello dello Zmap CurrentTriangle.ToGlob( m_MapFrame) ; - // Aggiungo alla lista + // Aggiungo alla lista lstTria.emplace_back( CurrentTriangle) ; - } + } continue ; } } @@ -2079,26 +2124,20 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold // 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) { - 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 ; 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[nCompCount - 1] - 1 ; ++ ni) { @@ -2107,7 +2146,6 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold nNegDot ++ ; } } - 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 @@ -2122,14 +2160,11 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold break ; } } - // Se tale normale esiste if ( bInversione2) { - // 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, ptSol)) { - Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; double dParInt1, dParInt2 ; Point3d ptVoxMin( ( i * N_DEXVOXRATIO + 0.5) * m_dStep, @@ -2139,14 +2174,10 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold // 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( ptSol, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { - triContainer.resize( 0) ; - double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 : dParInt2 + ( dParInt1 - dParInt2) / 100 ; - ptSol += dPar * vtNullSpace ; - // Costruisco triangoli di prova for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; @@ -2175,18 +2206,14 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold 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( ptSol, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) { - // Classificazione del voxel adiacente int nAdjIndex = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ; - // Se il voxel è pieno if ( EdgeTable[nAdjIndex] != 0) bDangerInversion = true ; @@ -2197,84 +2224,133 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold } } } - - // Ventaglio con base a quattro vertici + // 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. + // Controllo preliminare sulle normali bool bInversione = false ; - for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { - - int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; - + 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 ; + double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ; + if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) { + bInversione = true ; + break ; + } } - // Se tale triangolo esiste continuo i test + // Si necessita di un ulteriore controllo if ( bInversione) { - - // Se la feature non cade nel suo voxel - if ( ! IsPointInsideVoxelApprox( i, j, k, ptSol)) { - - 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( ptSol, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { - - triContainer.resize( 0) ; - - double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 : - dParInt2 + ( dParInt1 - dParInt2) / 100 ; - - ptSol += dPar * vtNullSpace ; - // 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 - Triangle3dEx CurrentTriangle ; - CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; - CurrentTriangle.SetAttrib( 1, CompoVert[nCompCount - 1][nj].nToolFlag) ; - CurrentTriangle.SetAttrib( 2, CompoVert[nCompCount - 1][ni].nToolFlag) ; - CurrentTriangle.SetVertexNorm( 1, CompoVert[nCompCount - 1][nj].vtNorm) ; - CurrentTriangle.SetVertexNorm( 2, CompoVert[nCompCount - 1][ni].vtNorm) ; - if ( CurrentTriangle.GetVertexNorm( 1) * CurrentTriangle.GetVertexNorm( 2) > 0.5) - CurrentTriangle.SetVertexNorm( 0, 0.5 * ( CurrentTriangle.GetVertexNorm( 1) + - CurrentTriangle.GetVertexNorm( 2))) ; - CurrentTriangle.Validate( true) ; - // Assegno i colori - if ( CurrentTriangle.GetAttrib( 1) < 0 || - CurrentTriangle.GetAttrib( 2) < 0) - CurrentTriangle.SetGrade( - 1) ; - else if ( CurrentTriangle.GetAttrib( 1) > 0 || - CurrentTriangle.GetAttrib( 2) > 0) - CurrentTriangle.SetGrade( 1) ; - else - CurrentTriangle.SetGrade( 0) ; - // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; - } + PolygonPlane Polygon ; + Point3d ptP0 = CompoVert[nCompCount - 1][0].ptInt ; + Point3d ptP1 = CompoVert[nCompCount - 1][1].ptInt ; + Point3d ptP2 = CompoVert[nCompCount - 1][2].ptInt ; + Point3d ptP3 = CompoVert[nCompCount - 1][3].ptInt ; + Polygon.AddPoint( ptP0) ; + Polygon.AddPoint( ptP1) ; + Polygon.AddPoint( ptP2) ; + Polygon.AddPoint( ptP3) ; + Point3d ptCentr ; + Vector3d vtPolyN ; + bool bCentFound = Polygon.GetCentroid( ptCentr) ; + bool bNormFound = Polygon.GetNormal( vtPolyN) ; + if ( bCentFound && bNormFound) { + Point3d ptPoly0 = ptP0 + ( ptP0 - ptCentr) - ( ptP0 - ptCentr) * vtPolyN * vtPolyN ; + Point3d ptPoly1 = ptP1 + ( ptP1 - ptCentr) - ( ptP1 - ptCentr) * vtPolyN * vtPolyN ; + Point3d ptPoly2 = ptP2 + ( ptP2 - ptCentr) - ( ptP2 - ptCentr) * vtPolyN * vtPolyN ; + Point3d ptPoly3 = ptP3 + ( ptP3 - ptCentr) - ( ptP3 - ptCentr) * vtPolyN * vtPolyN ; + Point3d ptSolP = ptSol + ( ptSol - ptCentr) - ( ptSol - ptCentr) * vtPolyN * vtPolyN ; + Vector3d vtSeg0 = ptPoly1 - ptPoly0 ; + Vector3d vtSeg1 = ptPoly2 - ptPoly1 ; + Vector3d vtSeg2 = ptPoly3 - ptPoly2 ; + Vector3d vtSeg3 = ptPoly0 - ptPoly3 ; + Vector3d vtSol0 = ptSolP - ptPoly0 ; + Vector3d vtSol1 = ptSolP - ptPoly1 ; + Vector3d vtSol2 = ptSolP - ptPoly2 ; + Vector3d vtSol3 = ptSolP - ptPoly3 ; + vtSeg0.Normalize() ; + vtSeg1.Normalize() ; + vtSeg2.Normalize() ; + vtSeg3.Normalize() ; + /*vtSol0.Normalize() ; + vtSol1.Normalize() ; + vtSol2.Normalize() ; + vtSol3.Normalize() ;*/ + double dIn[4] ; + dIn[0] = ( vtSeg0 ^ vtSol0) * vtPolyN ; + dIn[1] = ( vtSeg1 ^ vtSol1) * vtPolyN ; + dIn[2] = ( vtSeg2 ^ vtSol2) * vtPolyN ; + dIn[3] = ( vtSeg3 ^ vtSol3) * vtPolyN ; + double dThr = 0.1 * N_DEXVOXRATIO * m_dStep ; + bool bSpecial = false ; + for ( int nE = 0 ; nE < 4 ; ++ nE) { + if ( dIn[nE] < dThr) { + int nF = nE + 1 < 4 ? nE + 1 : 0 ; + double dDE = triContainer[nE].GetN() * CompoVert[nCompCount - 1][nE].vtNorm ; + double dDF = triContainer[nE].GetN() * CompoVert[nCompCount - 1][nF].vtNorm ; + double dLim = - 0.5/*EPS_SMALL*/ ; + if ( dDE < dLim || dDF < dLim) { + bSpecial = true ; + break ; + } + } } - - // 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 ; - } + if ( bSpecial) { + // Se la feature non cade nel suo voxel + if ( ! IsPointInsideVoxelApprox( i, j, k, ptSol)) { + // Vettore di direzione feature + 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( ptSol, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { + triContainer.resize( 0) ; + double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 : + dParInt2 + ( dParInt1 - dParInt2) / 100 ; + ptSol += dPar * vtNullSpace ; + // 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 + Triangle3dEx CurrentTriangle ; + CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; + CurrentTriangle.SetAttrib( 1, CompoVert[nCompCount - 1][nj].nToolFlag) ; + CurrentTriangle.SetAttrib( 2, CompoVert[nCompCount - 1][ni].nToolFlag) ; + CurrentTriangle.SetVertexNorm( 1, CompoVert[nCompCount - 1][nj].vtNorm) ; + CurrentTriangle.SetVertexNorm( 2, CompoVert[nCompCount - 1][ni].vtNorm) ; + if ( CurrentTriangle.GetVertexNorm( 1) * CurrentTriangle.GetVertexNorm( 2) > 0.5) + CurrentTriangle.SetVertexNorm( 0, 0.5 * ( CurrentTriangle.GetVertexNorm( 1) + + CurrentTriangle.GetVertexNorm( 2))) ; + CurrentTriangle.Validate( true) ; + // Assegno i colori + if ( CurrentTriangle.GetAttrib( 1) < 0 || + CurrentTriangle.GetAttrib( 2) < 0) + CurrentTriangle.SetGrade( - 1) ; + else if ( CurrentTriangle.GetAttrib( 1) > 0 || + CurrentTriangle.GetAttrib( 2) > 0) + CurrentTriangle.SetGrade( 1) ; + else + CurrentTriangle.SetGrade( 0) ; + // 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 ; + } + } + } } } } @@ -2408,7 +2484,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold } // ExtMC non confermato, si passa a MC - else { + else { vector vTria ; // Costruzione dei triangoli for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) { @@ -2459,7 +2535,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold } // Standard MC - else { + else { vector vTria ; // Costruzione dei triangoli for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {