diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index 85102c0..864a5d8 100644 Binary files a/EgtGeomKernel.rc and b/EgtGeomKernel.rc differ diff --git a/PolyLine.cpp b/PolyLine.cpp index 027cda5..4498ec9 100644 --- a/PolyLine.cpp +++ b/PolyLine.cpp @@ -734,6 +734,49 @@ PolyLine::RemoveAlignedPoints( double dToler) //---------------------------------------------------------------------------- bool PolyLine::ApproxOnSide( const Vector3d& vtN, bool bLeftSide, double dToler) +{ + // eseguo una prima approssimazione + if ( ! MyApproxOnSide( vtN, bLeftSide, dToler)) + return false ; + + // se polilinea aperta, finito + if ( ! IsClosed()) + return true ; + + // *** polilinea chiusa : devo sistemare a cavallo dell'inizio/fine *** + + // sposto l'inizio dalla parte opposta + MyChangeStartOnMiddle() ; + // rifaccio l'approssimazione con 3/4 della tolleranza ammessa + if ( ! MyApproxOnSide( vtN, bLeftSide, 0.75 * dToler)) + return false ; + // risposto l'inizio dalla parte opposta, quindi lo porto vicino all'originale + MyChangeStartOnMiddle() ; + + return true ; +} + +//---------------------------------------------------------------------------- +bool +PolyLine::MyChangeStartOnMiddle( void) +{ + // solo per polilinee chiuse + if ( ! IsClosed()) + return false ; + // cancello ultimo punto ( coincide con primo) + m_lUPoints.pop_back() ; + // sposto la metà iniziale dei punti alla fine + for ( size_t i = 0 ; i < m_lUPoints.size() / 2 ; ++ i) + m_lUPoints.splice( m_lUPoints.end(), m_lUPoints, m_lUPoints.begin()) ; + // aggiungo punto finale come copia dell'iniziale + m_lUPoints.push_back( m_lUPoints.front()) ; + + return true ; +} + +//---------------------------------------------------------------------------- +bool +PolyLine::MyApproxOnSide( const Vector3d& vtN, bool bLeftSide, double dToler) { // se non ci sono almeno 3 punti, esco subito if ( m_lUPoints.size() < 3) @@ -841,6 +884,31 @@ PolyLine::ApproxOnSide( const Vector3d& vtN, bool bLeftSide, double dToler) //---------------------------------------------------------------------------- bool PolyLine::MakeConvex( const Vector3d& vtN, bool bLeftSide) +{ + // eseguo una prima modifica + if ( ! MyMakeConvex( vtN, bLeftSide)) + return false ; + + // se polilinea aperta, finito + if ( ! IsClosed()) + return true ; + + // *** polilinea chiusa : devo sistemare a cavallo dell'inizio/fine *** + + // sposto l'inizio dalla parte opposta + MyChangeStartOnMiddle() ; + // rifaccio la modifica + if ( ! MyMakeConvex( vtN, bLeftSide)) + return false ; + // risposto l'inizio dalla parte opposta, quindi lo porto vicino all'originale + MyChangeStartOnMiddle() ; + + return true ; +} + +//---------------------------------------------------------------------------- +bool +PolyLine::MyMakeConvex( const Vector3d& vtN, bool bLeftSide) { // ciclo i controlli finchè non ci sono rimozioni bool bRemoved = true ; diff --git a/VolTriZmapCalculus.cpp b/VolTriZmapCalculus.cpp index bac53c0..cabd1c6 100644 --- a/VolTriZmapCalculus.cpp +++ b/VolTriZmapCalculus.cpp @@ -71,170 +71,6 @@ VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV, return ( dU2 >= dU1) ; } -/* -//---------------------------------------------------------------------------- -bool -VolZmap::IntersLineVoxel( const Point3d& ptP, const Vector3d& vtV, int nIndI, int nIndJ, int nIndK, - int& nFace1, int& nFace2, double& dU1, double& dU2) const -{ - // Controllo sull'ammissibilità del voxel - if ( nIndI < - 1 || nIndI >= int( m_nNx[0]) || - nIndJ < - 1 || nIndJ >= int( m_nNy[0]) || - nIndK < - 1 || nIndK >= int( m_nNy[1])) - return false ; - - Point3d ptInt ; - int nIntNum = 0 ; - double dU ; - - // Intersezione con le facce 1 e 3 - if ( abs( vtV.y) > EPS_ZERO) { - - // Intersezione con la prima faccia - dU1 = ( ( nIndJ + 0.5) * m_dStep - ptP.y) / vtV.y ; - ptInt = ptP + dU1 * vtV ; - - if ( ptInt.x > ( nIndI + 0.5) * m_dStep - EPS_SMALL && - ptInt.x < ( nIndI + 1.5) * m_dStep + EPS_SMALL && - ptInt.z > ( nIndK + 0.5) * m_dStep - EPS_SMALL && - ptInt.z < ( nIndK + 1.5) * m_dStep + EPS_SMALL) { - - nFace1 = 1 ; - ++ nIntNum ; - } - - // Intersezione con la terza faccia - dU = ( ( nIndJ + 1.5) * m_dStep - ptP.y) / vtV.y ; - ptInt = ptP + dU * vtV ; - - if ( ptInt.x > ( nIndI + 0.5) * m_dStep - EPS_SMALL && - ptInt.x < ( nIndI + 1.5) * m_dStep + EPS_SMALL && - ptInt.z > ( nIndK + 0.5) * m_dStep - EPS_SMALL && - ptInt.z < ( nIndK + 1.5) * m_dStep + EPS_SMALL) { - - if ( nIntNum == 0) { - - dU1 = dU ; - nFace1 = 3 ; - ++ nIntNum ; - } - else if ( nIntNum == 1) { - - dU2 = dU ; - nFace2 = 3 ; - ++ nIntNum ; - } - } - } - - // Intersezione con le facce 2 e 4 - if ( abs( vtV.x) > EPS_ZERO) { - - // Intersezione con la seconda faccia - dU = ( ( nIndI + 0.5) * m_dStep - ptP.x) / vtV.x ; - ptInt = ptP + dU * vtV ; - - if ( ptInt.y > ( nIndJ + 0.5) * m_dStep - EPS_SMALL && - ptInt.y < ( nIndJ + 1.5) * m_dStep + EPS_SMALL && - ptInt.z > ( nIndK + 0.5) * m_dStep - EPS_SMALL && - ptInt.z < ( nIndK + 1.5) * m_dStep + EPS_SMALL) { - - if ( nIntNum == 0) { - - dU1 = dU ; - nFace1 = 2 ; - ++ nIntNum ; - } - else if ( nIntNum == 1) { - - dU2 = dU ; - nFace2 = 2 ; - ++ nIntNum ; - } - } - - // Intersezione con la quarta faccia - dU = ( ( nIndI + 1.5) * m_dStep - ptP.x) / vtV.x ; - ptInt = ptP + dU * vtV ; - - if ( ptInt.y > ( nIndJ + 0.5) * m_dStep - EPS_SMALL && - ptInt.y < ( nIndJ + 1.5) * m_dStep + EPS_SMALL && - ptInt.z > ( nIndK + 0.5) * m_dStep - EPS_SMALL && - ptInt.z < ( nIndK + 1.5) * m_dStep + EPS_SMALL) { - - if ( nIntNum == 0) { - - dU1 = dU ; - nFace1 = 4 ; - ++ nIntNum ; - } - else if ( nIntNum == 1) { - - dU2 = dU ; - nFace2 = 4 ; - ++ nIntNum ; - } - } - } - - // Intersezione con le facce 5 e 6 - if ( abs( vtV.z) > EPS_ZERO) { - - // Intersezione con la quinta faccia - dU = ( ( nIndK + 0.5) * m_dStep - ptP.z) / vtV.z ; - ptInt = ptP + dU * vtV ; - - if ( ptInt.x > ( nIndI + 0.5) * m_dStep - EPS_SMALL && - ptInt.x < ( nIndI + 1.5) * m_dStep + EPS_SMALL && - ptInt.y > ( nIndJ + 0.5) * m_dStep - EPS_SMALL && - ptInt.y < ( nIndJ + 1.5) * m_dStep + EPS_SMALL) { - - if ( nIntNum == 0) { - - dU1 = dU ; - nFace1 = 5 ; - ++ nIntNum ; - } - else if ( nIntNum == 1) { - - dU2 = dU ; - nFace2 = 5 ; - ++ nIntNum ; - } - } - - // Intersezione con la sesta faccia - dU = ( ( nIndK + 1.5) * m_dStep - ptP.z) / vtV.z ; - ptInt = ptP + dU * vtV ; - - if ( ptInt.x > ( nIndI + 0.5) * m_dStep - EPS_SMALL && - ptInt.x < ( nIndI + 1.5) * m_dStep + EPS_SMALL && - ptInt.y > ( nIndJ + 0.5) * m_dStep - EPS_SMALL && - ptInt.y < ( nIndJ + 1.5) * m_dStep + EPS_SMALL) { - - - if ( nIntNum == 0) { - - dU1 = dU ; - nFace1 = 6 ; - ++ nIntNum ; - } - else if ( nIntNum == 1) { - - dU2 = dU ; - nFace2 = 6 ; - ++ nIntNum ; - } - } - } - - if ( dU1 > dU2) { - - swap( dU1, dU2) ; - swap( nFace1, nFace2) ; - } - return ( nIntNum == 2) ; -} */ //---------------------------------------------------------------------------- bool @@ -888,188 +724,6 @@ VolZmap::IntersZLineCylinder( const Point3d& ptLine, } } -/* -//---------------------------------------------------------------------------- -bool -VolZmap::IntersLineConus( const Point3d& ptLineSt, const Vector3d& vtLineDir, - const Frame3d& ConusFrame, double dTan, double dl, double dL, - Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2) -{ - // NB: L'origine del sistema di riferimento deve essere - // nel vertice del cono e l'asse di simmetria deve coincidere - // con l'asse x. - // La funzione restituisce true in caso di intersezione, - // false altrimenti. - - Point3d ptP = ptLineSt ; - Vector3d vtV = vtLineDir ; - - // Trasformazione delle coordinate - ptP.ToLoc( ConusFrame) ; - vtV.ToLoc( ConusFrame) ; - - DBLVECTOR vdCoef(3) ; - DBLVECTOR vdRoots ; - - double dSqTan = dTan * dTan ; - double dMinRad = dTan * dl ; - double dMaxRad = dTan * dL ; - double dDeltaR = dMaxRad - dMinRad ; - double dHei = dL - dl ; - - vdCoef[0] = dSqTan * ptP.x * ptP.x - ptP.y * ptP.y - ptP.z * ptP.z ; - vdCoef[1] = 2 * ( dSqTan * ptP.x * vtV.x - ptP.y * vtV.y - ptP.z * vtV.z) ; - vdCoef[2] = dSqTan * vtV.x * vtV.x - vtV.y * vtV.y - vtV.z * vtV.z ; - - // Computo radici - int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ; - - // Nessuna soluzione - if ( nRoot == 0) - return false ; - - // Una soluzione: la retta iterseca superficie - // laterale e un piano - if ( nRoot == 1) { - - ptInt1 = ptP + vdRoots[0] * vtV ; - - Vector3d vtU = ( ptInt1 - ORIG) - ( ptInt1 - ORIG).x * X_AX ; - - vtU.Normalize() ; - - vtN1 = dDeltaR * X_AX - dHei * vtU ; - - vtN1.Normalize() ; - - if ( ptInt1.x < dL + EPS_SMALL) { - - if ( ptInt1.x > dl) { - - ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; - - vtN2 = - X_AX ; - - } - else if ( ptInt1.x > - EPS_SMALL) { - - ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ; - ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; - - vtN1 = X_AX ; - vtN2 = - X_AX ; - - if ( ptInt2.y * ptInt2.y + ptInt2.z * ptInt2.z > dMaxRad * dMaxRad) - - return false ; - } - else - return false ; - - ptInt1.ToGlob( ConusFrame) ; - ptInt2.ToGlob( ConusFrame) ; - - vtN1.ToGlob( ConusFrame) ; - vtN2.ToGlob( ConusFrame) ; - - return true ; - } - else - return false ; - - } - // Due soluzioni: la retta interseca due volte la - // superficie laterale - else if ( nRoot == 2) { - - ptInt1 = ptP + vdRoots[0] * vtV ; - ptInt2 = ptP + vdRoots[1] * vtV ; - - if ( ptInt1.x > ptInt2.x) - swap( ptInt1, ptInt2) ; - - Vector3d vtU1 = ( ptInt1 - ORIG) - ( ptInt1 - ORIG).x * X_AX ; - Vector3d vtU2 = ( ptInt2 - ORIG) - ( ptInt2 - ORIG).x * X_AX ; - - vtU1.Normalize() ; - vtU2.Normalize() ; - - vtN1 = dDeltaR * X_AX - dHei * vtU1 ; - vtN2 = dDeltaR * X_AX - dHei * vtU2 ; - - vtN1.Normalize() ; - vtN2.Normalize() ; - - - if ( ptInt1.x < dL + EPS_SMALL) { - - if ( ptInt1.x > dl - EPS_SMALL) { - - if ( ptInt2.x > dL) { - - ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; - - vtN2 = - X_AX ; - } - } - else if ( ptInt1.x > - EPS_SMALL) { - - if ( ptInt2.x > dL) { - - ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ; - ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; - - vtN1 = X_AX ; - vtN2 = - X_AX ; - } - else if ( ptInt2.x > dl - EPS_SMALL) { - - ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ; - vtN1 = X_AX ; - } - else - return false ; - } - else { - - if ( ptInt2.x < 0) - - return false ; - - else if ( ptInt2.x < dl) { - - ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ; - ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; - - vtN1 = X_AX ; - vtN2 = - X_AX ; - } - else if ( ptInt2.x < dL) { - - ptInt1 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; - vtN1 = - X_AX ; - } - else - return false ; - } - - ptInt1.ToGlob( ConusFrame) ; - ptInt2.ToGlob( ConusFrame) ; - - vtN1.ToGlob( ConusFrame) ; - vtN2.ToGlob( ConusFrame) ; - - vtN1.Normalize() ; - vtN2.Normalize() ; - - return true ; - } - else - return false ; - } - return false ; -}*/ - //---------------------------------------------------------------------------- bool VolZmap::IntersLineConus( const Point3d& ptLineSt, const Vector3d& vtLineDir, @@ -1616,150 +1270,3 @@ VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLine else return false ; } -/* -//---------------------------------------------------------------------------- -bool -VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir, - const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaX, - Point3d& ptInt1, Point3d& ptInt2) -{ - - double SqIndet = EPS_SMALL * EPS_SMALL ; - - Point3d ptP = ptLineSt ; - Vector3d vtV = vtLineDir ; - - // Trasformazione delle coordinate - ptP.ToLoc( PolyFrame) ; - vtV.ToLoc( PolyFrame) ; - - // Facce 1 e 2 parallele a XY - // Facce 3 e 4 parallele a XZ - // Facce 5 e 6 oblique - Point3d ptFacet135( 0, 0, dLenZ /2) ; - Point3d ptFacet246( dLenX + dDeltaX, dLenY, - dLenZ / 2) ; - - // Servono per descrivere i piani obliqui - Vector3d vtFacet5 = ptFacet135 - ptP ; - Vector3d vtFacet6 = ptFacet246 - ptP ; - - Vector3d vtOb( dLenY, - dDeltaX, 0) ; - - Point3d ptI1 = ptP + ( ( ptFacet135.z - ptP.z) / vtV.z) * vtV ; - Point3d ptI2 = ptP + ( ( ptFacet246.z - ptP.z) / vtV.z) * vtV ; - Point3d ptI3 = ptP + ( ( ptFacet135.y - ptP.y) / vtV.y) * vtV ; - Point3d ptI4 = ptP + ( ( ptFacet246.y - ptP.y) / vtV.y) * vtV ; - Point3d ptI5 = ptP + ( ( vtFacet5 * vtOb) / ( vtV * vtOb)) * vtV ; - Point3d ptI6 = ptP + ( ( vtFacet6 * vtOb) / ( vtV * vtOb)) * vtV ; - - int nIntNum = 0 ; - - // Intersezione con la prima faccia - if ( ptI1.y > - EPS_SMALL && ptI1.y < dLenY + EPS_SMALL && - ptI1.x * dLenY > dDeltaX * ptI1.y - EPS_SMALL && - ( ptI1.x - dLenX) * dLenY < dDeltaX * ptI1.y + EPS_SMALL) { - - ptInt1 = ptI1 ; - ++ nIntNum ; - } - - // Intersezione con la seconda faccia - if ( ptI2.y > - EPS_SMALL && ptI2.y dDeltaX * ptI2.y - EPS_SMALL && - ( ptI2.x - dLenX) * dLenY < dDeltaX * ptI2.y + EPS_SMALL) { - - if ( nIntNum == 0) { - - ptInt1 = ptI2 ; - ++ nIntNum ; - } - else if ( ( ptInt1 - ptI2).SqLen() > SqIndet) { - - ptInt2 = ptI2 ; - ++ nIntNum ; - } - } - - // Intersezione con la terza faccia - if ( nIntNum < 2 && - ptI3.x > - EPS_SMALL && ptI3.x < dLenX + EPS_SMALL && - ptI3.z > - ptFacet135.z - EPS_SMALL && - ptI3.z < ptFacet135.z + EPS_SMALL) { - - if ( nIntNum == 0) { - - ptInt1 = ptI3 ; - ++ nIntNum ; - } - else if ( ( ptInt1 - ptI3).SqLen() > SqIndet) { - - ptInt2 = ptI3 ; - ++ nIntNum ; - } - } - - // Intersezione con la quarta faccia - if ( nIntNum < 2 && - ptI4.x > dDeltaX - EPS_SMALL && ptI4.x < dLenX + dDeltaX + EPS_SMALL && - ptI4.z > - ptFacet135.z - EPS_SMALL && - ptI4.z < ptFacet135.z + EPS_SMALL) { - - if ( nIntNum == 0) { - - ptInt1 = ptI4 ; - ++ nIntNum ; - } - else if ( ( ptInt1 - ptI4).SqLen() > SqIndet) { - - ptInt2 = ptI4 ; - ++ nIntNum ; - } - } - - // Intersezione con la quinta faccia - if ( nIntNum < 2 && - ptI5.y > - EPS_SMALL && ptI5.y < dLenY + EPS_SMALL && - ptI5.z > - ptFacet135.z - EPS_SMALL && - ptI5.z < ptFacet135.z + EPS_SMALL) { - - if ( nIntNum == 0) { - - ptInt1 = ptI5 ; - ++ nIntNum ; - } - else if ( ( ptInt1 - ptI5).SqLen() > SqIndet) { - - ptInt2 = ptI5 ; - ++ nIntNum ; - } - } - - // Intersezione con la sesta faccia - if ( nIntNum < 2 && - ptI6.y > - EPS_SMALL && ptI6.y < dLenY + EPS_SMALL && - ptI6.z > - ptFacet135.z - EPS_SMALL && - ptI6.z < ptFacet135.z + EPS_SMALL) { - - if ( nIntNum == 0) { - - ptInt1 = ptI6; - ++ nIntNum ; - } - else if ( ( ptInt1 - ptI6).SqLen() > SqIndet) { - - ptInt2 = ptI6; - ++ nIntNum ; - } - } - - - if ( nIntNum == 2) { - - ptInt1.ToGlob( PolyFrame) ; - ptInt2.ToGlob( PolyFrame) ; - - return true ; - } - else - return false ; -} */ \ No newline at end of file diff --git a/VolTriZmapCreation.cpp b/VolTriZmapCreation.cpp index 0bf1aef..0a63b04 100644 --- a/VolTriZmapCreation.cpp +++ b/VolTriZmapCreation.cpp @@ -190,8 +190,6 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double // A partire dalle dimensioni di xy del grezzo determino il numero di colonne e righe // della griglia Zmap e da questi la dimensione del vettore di dexel - //m_nNx[0] = static_cast ( ceil( dLengthX / m_dStep)) ; - //m_nNy[0] = static_cast ( ceil( dLengthY / m_dStep)) ; m_nNx[0] = static_cast ( floor( ( dLengthX + 0.5 * m_dStep + EPS_SMALL) / m_dStep)) ; m_nNy[0] = static_cast ( floor( ( dLengthY + 0.5 * m_dStep + EPS_SMALL) / m_dStep)) ; @@ -216,12 +214,6 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double m_MapFrame[1].Set( ptMapOrig, Y_AX, Z_AX, X_AX) ; m_MapFrame[2].Set( ptMapOrig, Z_AX, X_AX, Y_AX) ; - //m_nNx[1] = static_cast ( ceil( dLengthY / m_dStep)) ; - //m_nNy[1] = static_cast ( ceil( dDimZ / m_dStep)) ; - - //m_nNx[2] = static_cast ( ceil( dDimZ / m_dStep)) ; - //m_nNy[2] = static_cast ( ceil( dLengthX / m_dStep)) ; - m_nNx[1] = static_cast ( floor( ( dLengthY + 0.5 * m_dStep + EPS_SMALL) / m_dStep)) ; m_nNy[1] = static_cast ( floor( ( dDimZ + 0.5 * m_dStep + EPS_SMALL) / m_dStep)) ; @@ -264,9 +256,9 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double // Definisco la retta da intersecare con la regione double dX = ( i + 0.5) * m_dStep ; - Point3d ptP0 = ptMapOrig + Vector3d( dX, 0, 0) ; // CAMBIATO CON EPS + Point3d ptP0 = ptMapOrig + Vector3d( dX, 0, 0) ; CurveLine GridLine ; - GridLine.SetPVL( ptP0, Y_AX, dLengthY) ; // CAMBIATO CON EPS + GridLine.SetPVL( ptP0, Y_AX, dLengthY) ; // Determino le intersezioni della retta con la regione CRVCVECTOR IntersectionResults ; @@ -288,8 +280,8 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double if ( nType == CRVC_IN || nType == CRVC_ON_P || nType == CRVC_ON_M) { // Indici corrispondenti alle coordinate dei punti - int nStartJ = Clamp( int( floor( dt1 * dLengthY / m_dStep - EPS_SMALL + 0.5)), 0, m_nNy[0] - 1) ; // CAMBIATO CON EPS - int nEndJ = Clamp( int( floor( dt2 * dLengthY / m_dStep + EPS_SMALL - 0.5)), 0, m_nNy[0] - 1) ; // CAMBIATO CON EPS + int nStartJ = Clamp( int( floor( dt1 * dLengthY / m_dStep - EPS_SMALL + 0.5)), 0, m_nNy[0] - 1) ; + int nEndJ = Clamp( int( floor( dt2 * dLengthY / m_dStep + EPS_SMALL - 0.5)), 0, m_nNy[0] - 1) ; // Ridimensiono e riempio i dexel for ( int j = nStartJ ; j <= nEndJ ; ++ j) { @@ -319,8 +311,8 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double m_Values[2][nPos2][nCurrentSize].dZVal = dt1 * dLengthY ; m_Values[2][nPos2][nCurrentSize + 1].dZVal = dt2 * dLengthY ; - Point3d ptP1 = ptP0 + dt1 * dLengthY * Y_AX ;//+ ( a + 0.5) * Z_AX ; // CAMBIATO CON EPS - Point3d ptP2 = ptP0 + dt2 * dLengthY * Y_AX ;//+ ( a + 0.5) * Z_AX ; // CAMBIATO CON EPS + Point3d ptP1 = ptP0 + dt1 * dLengthY * Y_AX ;//+ ( a + 0.5) * Z_AX ; + Point3d ptP2 = ptP0 + dt2 * dLengthY * Y_AX ;//+ ( a + 0.5) * Z_AX ; int nChunkNum = Surf.GetChunkCount() ; @@ -373,9 +365,9 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double // Definisco la retta da intersecare con la regione double dX = ( i + 0.5) * m_dStep ; - Point3d ptP0 = ptMapOrig + Vector3d( 0, dX, 0) ; ////// // CAMBIATO CON EPS + Point3d ptP0 = ptMapOrig + Vector3d( 0, dX, 0) ; CurveLine GridLine ; - GridLine.SetPVL( ptP0, X_AX, dLengthX ) ; //// // CAMBIATO CON EPS + GridLine.SetPVL( ptP0, X_AX, dLengthX ) ; // Determino le intersezioni della retta con la regione CRVCVECTOR IntersectionResults ; @@ -404,11 +396,11 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double m_Values[1][nPos1].resize( nCurrentSize + 2) ; - m_Values[1][nPos1][nCurrentSize].dZVal = dt1 * dLengthX ; // CAMBIATO CON EPS - m_Values[1][nPos1][nCurrentSize + 1].dZVal = dt2 * dLengthX ; // CAMBIATO CON EPS + m_Values[1][nPos1][nCurrentSize].dZVal = dt1 * dLengthX ; + m_Values[1][nPos1][nCurrentSize + 1].dZVal = dt2 * dLengthX ; - Point3d ptP1 = ptP0 + dt1 * dLengthX * X_AX ; /// cAMBIO - Point3d ptP2 = ptP0 + dt2 * dLengthX * X_AX ; // CAMBIATO CON EPS + Point3d ptP1 = ptP0 + dt1 * dLengthX * X_AX ; + Point3d ptP2 = ptP0 + dt2 * dLengthX * X_AX ; int nChunkNum = Surf.GetChunkCount() ; diff --git a/VolTriZmapGraphics.cpp b/VolTriZmapGraphics.cpp index 5bb27e9..8181bb2 100644 --- a/VolTriZmapGraphics.cpp +++ b/VolTriZmapGraphics.cpp @@ -52,58 +52,7 @@ static int NeighbourTable[8][4] = { // ------------------------- FUNZIONE TEST SULLE NORMALI -------------------------------------------------------------------------- -enum FatureType { NoFeature = 0, Corner = 1, Edge = 2} ; //, Edge2 = 3} ; - -//---------------------------------------------------------------------------- -//bool -//TestOnNormal( const VectorField CompoVert[], int nCompoElem, int& FeatureType, Vector3d& vtFeature) -//{ -// int nI, nJ ; -// double dMinCosTheta = 1.001 ; -// const double dCosThetaSharp = sqrt( 3) / 2 ; 0.9 ; -// -// for ( int i = 0 ; i < nCompoElem ; ++ i) { -// -// for ( int j = i + 1 ; j < nCompoElem ; ++ j) { -// -// double dCurrentCos = CompoVert[i].vtNorm * CompoVert[j].vtNorm ; -// -// if ( dCurrentCos < dMinCosTheta) { -// -// nI = i ; -// nJ = j ; -// dMinCosTheta = dCurrentCos ; -// } -// } -// } -// -// if ( dMinCosTheta >= dCosThetaSharp) { -// -// FeatureType = NoFeature ; -// return false ; -// } -// -// Vector3d vtK = CompoVert[nI].vtNorm ^ CompoVert[nJ].vtNorm ; -// -// const double dCosPhiCorner = 0.5;//0.7;//0.5 ; -// -// for ( int i = 0 ; i < nCompoElem ; ++ i) { -// -// double dAbsCurrentCos = abs( CompoVert[i].vtNorm * vtK) ; -// -// if ( dAbsCurrentCos > dCosPhiCorner) { -// -// FeatureType = Corner ; -// vtFeature = vtK ; -// return true ; -// } -// } -// -// FeatureType = Edge ; -// vtFeature = vtK ; -// -// return true ; -//} +enum FatureType { NoFeature = 0, Corner = 1, Edge = 2} ; //---------------------------------------------------------------------------- bool @@ -360,62 +309,6 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVE // Flipping fra voxel interni FlipEdgesII( VecTriHold[t]) ; - //for ( int nDeltaI = -1 ; nDeltaI <= 1 ; ++ nDeltaI) { - // for ( int nDeltaJ = -1 ; nDeltaJ <= 1 ; ++ nDeltaJ) { - // for ( int nDeltaK = -1 ; nDeltaK <= 1 ; ++ nDeltaK) { - // - // // Flipping fra voxel dello stesso blocco - // if ( nDeltaI == 0 && nDeltaJ == 0 && nDeltaK == 0) { - // - // // Flipping fra voxel interni - // FlipEdgesII( VecTriHold[t]) ; - // // Flipping fra voxel interni e di frontiera - // //FlipEdgesIB( VecTriHold[t], t) ; - // } - // // Flipping fra voxel adiacenti - // else { - - // int nAdjN ; - // - // if ( GetAdjBlockToBlock( int( t), nDeltaI, nDeltaJ, nDeltaK, nAdjN)) { - // - // // Qui bisogna fare il flipping fra i blocchi: t indica il blocco - // // corrente appena aggiornato. Guardiamo il suo vettore di voxel bislacchi - // // e cicliamo su tali voxel. Per ogni voxel bislacco guardiamo i voxel adiacenti. - // // Se un voxel adiacente è nel blocco adiacente corrente ricalcoliamo i triangoli e flippiamo. - // size_t nBoundaryVoxelNum = m_InterBlockTria[t].size() ; - - // for ( size_t a = 0 ; a < nBoundaryVoxelNum ; ++ a) { - // - // // vedere se il voxel appartiene al blocco adiacente - // int nBVoxI = m_InterBlockTria[t][a].i ; - // int nBVoxJ = m_InterBlockTria[t][a].j ; - // int nBVoxK = m_InterBlockTria[t][a].k ; - - // // Indici blocco adiacente - // int nAdjBlockIJK[3] ; - // GetBlockIJKFromN( nAdjN, nAdjBlockIJK) ; - // // Limiti indici di voxel per il blocco adiacente - // int nLimitAdjBlock[6] ; - - // GetBlockLimitsIJK( nAdjBlockIJK, nLimitAdjBlock) ; - // - // // Se il blocco adiacente esiste eseguo il flipping - // if ( nBVoxI >= nLimitAdjBlock[0] && nBVoxI < nLimitAdjBlock[1] && - // nBVoxJ >= nLimitAdjBlock[2] && nBVoxJ < nLimitAdjBlock[3] && - // nBVoxK >= nLimitAdjBlock[4] && nBVoxK < nLimitAdjBlock[5]) { - // - // m_InterBlockTria[nAdjN].clear() ; - - - - // } - // } - // } - // } - // } - // } - //} m_BlockToUpdate[t] = false ; } else @@ -894,9 +787,6 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& for ( int j = nLimits[2] ; j < nLimits[3] ; ++ j) { for ( int k = nLimits[4] ; k < nLimits[5] ; ++ k) { - if ( !(i == 31 && j == -1 && k == 49)) - ;// continue ; - // Riconoscimento dei voxel di frontiera bool bBoundary = false ; int nVoxIndexes[3] = { i, j, k} ; @@ -1011,164 +901,6 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& // Nota che il primo elemento dell'array (0-esimo) non viene inizializzato CompoVert[nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nTableOffset + 1]] ; - // Controllo per il caso piano su una griglia - // con versori normali a due a due paralleli. - bool bGridControl = true ; - - if ( nVertComp == 4) { - - int nIndArrey[4] ; - // Riempio un vettore con gli indici dei punti - for ( int nCpInd = 0 ; nCpInd < nVertComp ; ++ nCpInd) - nIndArrey[nCpInd] = TriangleTableEn[nIndex][1][nCpInd + nTableOffset + 1] ; - // Ordino i 4 indici in senso crescente - for ( int nSrtInd1 = 0 ; nSrtInd1 < nVertComp - 1 ; ++ nSrtInd1) { - for ( int nSrtInd2 = nSrtInd1 + 1 ; nSrtInd2 < nVertComp ; ++ nSrtInd2) { - if ( nIndArrey[nSrtInd1] > nIndArrey[nSrtInd2]) - swap( nIndArrey[nSrtInd1], nIndArrey[nSrtInd2]) ; - } - } - - if ( ( nIndArrey[0] == 0 && nIndArrey[1] == 2 && nIndArrey[2] == 9 && nIndArrey[3] == 10) || - ( nIndArrey[0] == 4 && nIndArrey[1] == 6 && nIndArrey[2] == 9 && nIndArrey[3] == 10) || - ( nIndArrey[0] == 4 && nIndArrey[1] == 6 && nIndArrey[2] == 8 && nIndArrey[3] == 11) || - ( nIndArrey[0] == 0 && nIndArrey[1] == 2 && nIndArrey[2] == 8 && nIndArrey[3] == 11) || - ( nIndArrey[0] == 1 && nIndArrey[1] == 3 && nIndArrey[2] == 8 && nIndArrey[3] == 9 ) || - ( nIndArrey[0] == 1 && nIndArrey[1] == 3 && nIndArrey[2] == 10 && nIndArrey[3] == 11) || - ( nIndArrey[0] == 5 && nIndArrey[1] == 7 && nIndArrey[2] == 10 && nIndArrey[3] == 11) || - ( nIndArrey[0] == 5 && nIndArrey[1] == 7 && nIndArrey[2] == 8 && nIndArrey[3] == 9 )) { - - if ( AreSameVectorApprox( VecField[nIndArrey[0]].vtNorm, VecField[nIndArrey[1]].vtNorm) && - AreSameVectorApprox( VecField[nIndArrey[2]].vtNorm, VecField[nIndArrey[3]].vtNorm) && - abs( VecField[nIndArrey[0]].vtNorm * VecField[nIndArrey[2]].vtNorm) < EPS_SMALL) { - - Point3d ptBarycenter = ( CompoVert[0].ptInt + CompoVert[1].ptInt + - CompoVert[2].ptInt + CompoVert[3].ptInt) / 4 ; - - if ( abs( CompoVert[0].ptInt.x - ptBarycenter.x) < EPS_SMALL && - abs( CompoVert[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && - abs( CompoVert[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && - abs( CompoVert[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) { - - double dXBar = ptBarycenter.x / m_dStep - 0.5 ; - int nBarLimSup = int( m_nNx[0]) ; - int nBarInd = 0 ; - - while ( nBarInd < nBarLimSup) { - - double dXInd = double( nBarInd) ; - - if ( abs( dXInd - dXBar) < EPS_SMALL) { - - bGridControl = false ; - break ; - } - ++ nBarInd ; - } - } - else if ( abs( CompoVert[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && - abs( CompoVert[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && - abs( CompoVert[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && - abs( CompoVert[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) { - - double dYBar = ptBarycenter.y / m_dStep - 0.5 ; - int nBarLimSup = int( m_nNy[0]) ; - int nBarInd = 0 ; - - while ( nBarInd < nBarLimSup) { - - double dYInd = double( nBarInd) ; - - if ( abs( dYInd - dYBar) < EPS_SMALL) { - - bGridControl = false ; - break ; - } - ++ nBarInd ; - } - } - else if ( abs( CompoVert[0].ptInt.z - ptBarycenter.z) < EPS_SMALL && - abs( CompoVert[1].ptInt.z - ptBarycenter.z) < EPS_SMALL && - abs( CompoVert[2].ptInt.z - ptBarycenter.z) < EPS_SMALL && - abs( CompoVert[3].ptInt.z - ptBarycenter.z) < EPS_SMALL) { - - double dZBar = ptBarycenter.z / m_dStep - 0.5 ; - int nBarLimSup = int( m_nNy[1]) ; - int nBarInd = 0 ; - - while ( nBarInd < nBarLimSup) { - - double dZInd = double( nBarInd) ; - - if ( abs( dZInd - dZBar) < EPS_SMALL) { - - bGridControl = false ; - break ; - } - ++ nBarInd ; - } - } - } - } - else if ( ( nIndArrey[0] == 0 && nIndArrey[1] == 1 && nIndArrey[2] == 4 && nIndArrey[3] == 5) || - ( nIndArrey[0] == 1 && nIndArrey[1] == 2 && nIndArrey[2] == 5 && nIndArrey[3] == 6) || - ( nIndArrey[0] == 2 && nIndArrey[1] == 3 && nIndArrey[2] == 6 && nIndArrey[3] == 7) || - ( nIndArrey[0] == 0 && nIndArrey[1] == 3 && nIndArrey[2] == 4 && nIndArrey[3] == 7)) { - - if ( AreSameVectorApprox( VecField[nIndArrey[0]].vtNorm, VecField[nIndArrey[2]].vtNorm) && - AreSameVectorApprox( VecField[nIndArrey[1]].vtNorm, VecField[nIndArrey[3]].vtNorm) && - abs( VecField[nIndArrey[0]].vtNorm * VecField[nIndArrey[1]].vtNorm) < EPS_SMALL) { - - Point3d ptBarycenter = ( CompoVert[0].ptInt + CompoVert[1].ptInt + - CompoVert[2].ptInt + CompoVert[3].ptInt) / 4 ; - - if ( abs( CompoVert[0].ptInt.x - ptBarycenter.x) < EPS_SMALL && - abs( CompoVert[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && - abs( CompoVert[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && - abs( CompoVert[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) { - - double dXBar = ptBarycenter.x / m_dStep - 0.5 ; - int nBarLimSup = int( m_nNx[0]) ; - int nBarInd = 0 ; - - while ( nBarInd < nBarLimSup) { - - double dXInd = double( nBarInd) ; - - if ( abs( dXInd - dXBar) < EPS_SMALL) { - - bGridControl = false ; - break ; - } - ++ nBarInd ; - } - } - else if ( abs( CompoVert[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && - abs( CompoVert[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && - abs( CompoVert[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && - abs( CompoVert[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) { - - double dYBar = ptBarycenter.y / m_dStep - 0.5 ; - int nBarLimSup = int( m_nNy[0]) ; - int nBarInd = 0 ; - - while ( nBarInd < nBarLimSup) { - - double dYInd = double( nBarInd) ; - - if ( abs( dYInd - dYBar) < EPS_SMALL) { - - bGridControl = false ; - break ; - } - ++ nBarInd ; - } - } - } - } - } - - int nFeatureType ; // Valuto le relazioni reciproche fra le normali e @@ -1186,6 +918,219 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& if ( bReg) bNormal = TestOnNormal( CompoVert, nVertComp, nFeatureType, vtFeatureAxis, nMin1, nMin2) ; + + // Controllo per il caso piano su una griglia + // con versori normali a due a due paralleli. + bool bGridControl = true ; + + if ( bNormal) { + + if ( nVertComp == 4) { + + int nIndArrey[4] ; + // Riempio un vettore con gli indici dei punti + for ( int nCpInd = 0 ; nCpInd < nVertComp ; ++ nCpInd) + nIndArrey[nCpInd] = TriangleTableEn[nIndex][1][nCpInd + nTableOffset + 1] ; + // Ordino i 4 indici in senso crescente + for ( int nSrtInd1 = 0 ; nSrtInd1 < nVertComp - 1 ; ++ nSrtInd1) { + for ( int nSrtInd2 = nSrtInd1 + 1 ; nSrtInd2 < nVertComp ; ++ nSrtInd2) { + if ( nIndArrey[nSrtInd1] > nIndArrey[nSrtInd2]) + swap( nIndArrey[nSrtInd1], nIndArrey[nSrtInd2]) ; + } + } + + if ( ( nIndArrey[0] == 0 && nIndArrey[1] == 2 && nIndArrey[2] == 9 && nIndArrey[3] == 10) || + ( nIndArrey[0] == 4 && nIndArrey[1] == 6 && nIndArrey[2] == 9 && nIndArrey[3] == 10) || + ( nIndArrey[0] == 4 && nIndArrey[1] == 6 && nIndArrey[2] == 8 && nIndArrey[3] == 11) || + ( nIndArrey[0] == 0 && nIndArrey[1] == 2 && nIndArrey[2] == 8 && nIndArrey[3] == 11) || + ( nIndArrey[0] == 1 && nIndArrey[1] == 3 && nIndArrey[2] == 8 && nIndArrey[3] == 9 ) || + ( nIndArrey[0] == 1 && nIndArrey[1] == 3 && nIndArrey[2] == 10 && nIndArrey[3] == 11) || + ( nIndArrey[0] == 5 && nIndArrey[1] == 7 && nIndArrey[2] == 10 && nIndArrey[3] == 11) || + ( nIndArrey[0] == 5 && nIndArrey[1] == 7 && nIndArrey[2] == 8 && nIndArrey[3] == 9 )) { + + VectorField LocVecF[12], LocCompV[12] ; + + for ( int LocInd = 0 ; LocInd < 12 ; ++ LocInd) { + + LocVecF[LocInd] = VecField[LocInd] ; + LocCompV[LocInd] = CompoVert[LocInd] ; + + LocVecF[LocInd].ptInt.ToLoc( m_MapFrame[0]) ; + LocVecF[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ; + LocCompV[LocInd].ptInt.ToLoc( m_MapFrame[0]) ; + LocCompV[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ; + } + + if ( ( AreSameVectorApprox( LocVecF[nIndArrey[0]].vtNorm, LocVecF[nIndArrey[1]].vtNorm) && + abs( LocVecF[nIndArrey[0]].vtNorm *LocVecF[nIndArrey[2]].vtNorm) < EPS_SMALL && + abs( LocVecF[nIndArrey[0]].vtNorm * LocVecF[nIndArrey[3]].vtNorm) < EPS_SMALL) || + ( AreSameVectorApprox( LocVecF[nIndArrey[2]].vtNorm, LocVecF[nIndArrey[3]].vtNorm) && + abs( LocVecF[nIndArrey[2]].vtNorm * LocVecF[nIndArrey[0]].vtNorm) < EPS_SMALL && + abs( LocVecF[nIndArrey[2]].vtNorm * LocVecF[nIndArrey[1]].vtNorm) < EPS_SMALL)) { + + Point3d ptBarycenter = ( LocCompV[0].ptInt + LocCompV[1].ptInt + + LocCompV[2].ptInt + LocCompV[3].ptInt) / 4 ; + + if ( abs( LocCompV[0].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) { + + double dXBar = ptBarycenter.x / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNx[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dXInd = double( nBarInd) ; + + if ( abs( dXInd - dXBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) { + + double dYBar = ptBarycenter.y / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNy[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dYInd = double( nBarInd) ; + + if ( abs( dYInd - dYBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[3].ptInt.z - ptBarycenter.z) < EPS_SMALL) { + + double dZBar = ptBarycenter.z / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNy[1]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dZInd = double( nBarInd) ; + + if ( abs( dZInd - dZBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + } + } + else if ( ( nIndArrey[0] == 0 && nIndArrey[1] == 1 && nIndArrey[2] == 4 && nIndArrey[3] == 5) || + ( nIndArrey[0] == 1 && nIndArrey[1] == 2 && nIndArrey[2] == 5 && nIndArrey[3] == 6) || + ( nIndArrey[0] == 2 && nIndArrey[1] == 3 && nIndArrey[2] == 6 && nIndArrey[3] == 7) || + ( nIndArrey[0] == 0 && nIndArrey[1] == 3 && nIndArrey[2] == 4 && nIndArrey[3] == 7)) { + + VectorField LocVecF[12], LocCompV[12] ; + + for ( int LocInd = 0 ; LocInd < 12 ; ++ LocInd) { + + LocVecF[LocInd] = VecField[LocInd] ; + LocCompV[LocInd] = CompoVert[LocInd] ; + + LocVecF[LocInd].ptInt.ToLoc( m_MapFrame[0]) ; + LocVecF[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ; + LocCompV[LocInd].ptInt.ToLoc( m_MapFrame[0]) ; + LocCompV[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ; + } + + if ( ( AreSameVectorApprox( LocVecF[nIndArrey[0]].vtNorm, LocVecF[nIndArrey[2]].vtNorm) && + abs( LocVecF[nIndArrey[0]].vtNorm * LocVecF[nIndArrey[1]].vtNorm) < EPS_SMALL && + abs( LocVecF[nIndArrey[0]].vtNorm * LocVecF[nIndArrey[3]].vtNorm) < EPS_SMALL) || + ( AreSameVectorApprox( LocVecF[nIndArrey[1]].vtNorm, LocVecF[nIndArrey[3]].vtNorm) && + abs( LocVecF[nIndArrey[1]].vtNorm * LocVecF[nIndArrey[0]].vtNorm) < EPS_SMALL && + abs( LocVecF[nIndArrey[1]].vtNorm * LocVecF[nIndArrey[2]].vtNorm) < EPS_SMALL)) { + + Point3d ptBarycenter = ( LocCompV[0].ptInt + LocCompV[1].ptInt + + LocCompV[2].ptInt + LocCompV[3].ptInt) / 4 ; + + if ( abs( LocCompV[0].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) { + + double dXBar = ptBarycenter.x / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNx[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dXInd = double( nBarInd) ; + + if ( abs( dXInd - dXBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) { + + double dYBar = ptBarycenter.y / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNy[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dYInd = double( nBarInd) ; + + if ( abs( dYInd - dYBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[3].ptInt.z - ptBarycenter.z) < EPS_SMALL) { + + double dZBar = ptBarycenter.z / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNy[1]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dZInd = double( nBarInd) ; + + if ( abs( dZInd - dZBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + } + } + } + } // Flag ExtMC bool bExtMC = bNormal && bGridControl ; @@ -1564,42 +1509,6 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& // non piana confermo ExtMC if ( ! ( bOutside || bPlane)) { - // 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, - // // quindi sono validi per tutte le - // // componenti che vi appartengono. - // int nCurrent = int( triHold.size()) - 1 ; - - // triHold[nCurrent].i = i ; - // triHold[nCurrent].j = j ; - // triHold[nCurrent].k = k ; - // } - - //// 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) ; - TRIA3DVECTOR vInnerTriaTemp, vBorderTriaTemp ; for ( int ni = 0 ; ni < nContSize ; ++ ni) { @@ -1765,442 +1674,6 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& return true ; } -//---------------------------------------------------------------------------- -//bool -//VolZmap::FlipEdges( std::vector& VecTriHold) const -//{ -// double dSqEps = EPS_SMALL * EPS_SMALL ; -// -// // Ciclo sui blocchi -// for ( size_t t1 = 0 ; t1 < m_nNumBlock ; ++ t1) { -// -// // Determino i,j,k del primo blocco -// int nIJK1[3] ; -// GetBlockIJK( int( t1), nIJK1) ; -// -// for ( size_t t2 = t1 ; t2 < m_nNumBlock ; ++ t2) { -// -// // Determino i,j,k del secondo blocco -// int nIJK2[3] ; -// GetBlockIJK( int( t2), nIJK2) ; -// -// // Se i blocchi sono adiacenti o coincidenti proseguo -// if ( ( nIJK2[0] >= nIJK1[0] - 1 && nIJK2[0] <= nIJK1[0] + 1) || -// ( nIJK2[1] >= nIJK1[1] - 1 && nIJK2[1] <= nIJK1[1] + 1) || -// ( nIJK2[2] >= nIJK1[2] - 1 && nIJK2[2] <= nIJK1[2] + 1)) { -// -// // Numero di voxel in cui si presentano sharp feature -// int nVoxelNum1 = int( VecTriHold[t1].size()) ; -// int nVoxelNum2 = int( VecTriHold[t2].size()) ; -// -// // Determino estremi del ciclo sui voxel esterno -// int nSt1 = 0 ; -// int nEn1 = nVoxelNum1 - ( t1 == t2 ? 1 : 0) ; -// -// // Ciclo su tali voxel -// for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { -// -// // Determino estremi del ciclo sui voxel interno -// int nSt2 = ( t1 == t2 ? nSt1 : 0) ; -// int nEn2 = nVoxelNum2 ; -// -// for ( int n2 = n1 ; n2 < nVoxelNum2 ; ++ n2) { -// -// // Se i voxel sono adiacenti proseguo -// if ( VecTriHold[t2][n2].i == VecTriHold[t1][n1].i - 1 || -// VecTriHold[t2][n2].i == VecTriHold[t1][n1].i + 1 || -// VecTriHold[t2][n2].j == VecTriHold[t1][n1].j - 1 || -// VecTriHold[t2][n2].j == VecTriHold[t1][n1].j + 1 || -// VecTriHold[t2][n2].k == VecTriHold[t1][n1].k - 1 || -// VecTriHold[t2][n2].k == VecTriHold[t1][n1].k + 1) { -// -// // Numero delle componenti connesse nei due voxel -// int nNumCompo1 = int( VecTriHold[t1][n1].ptCompoVert.size()) ; -// int nNumCompo2 = int( VecTriHold[t2][n2].ptCompoVert.size()) ; -// -// int nCompo1 = 0 ; -// -// // Ciclo sulle componenti -// for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { -// -// int nCompo2 = ( t1 == t2 && n1 == n2 ? nCompo1 + 1 : 0) ; -// -// for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { -// -// // Numero di triangoli per le componenti connesse -// int nTriNum1 = int( VecTriHold[t1][n1].vCompoTria[nCompo1].size()) ; -// int nTriNum2 = int( VecTriHold[t2][n2].vCompoTria[nCompo2].size()) ; -// -// for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) { -// -// bool bModified = false ; -// -// for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) { -// -// INTVECTOR SharedIndex ; -// -// for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) { -// -// for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) { -// -// Point3d ptP1 = VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; -// Point3d ptP2 = VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; -// -// if ( SqDist( ptP1, ptP2) < dSqEps) { -// -// Point3d ptVert1 = VecTriHold[t1][n1].ptCompoVert[nCompo1] ; -// Point3d ptVert2 = VecTriHold[t2][n2].ptCompoVert[nCompo2] ; -// -// if ( ! ( AreSamePointApprox( ptP1, ptVert1) || -// AreSamePointApprox( ptP2, ptVert2))) { -// -// SharedIndex.emplace_back( nVert1) ; -// SharedIndex.emplace_back( nVert2) ; -// } -// } -// -// if ( SharedIndex.size() > 2) -// break ; -// } -// -// if ( SharedIndex.size() > 2) -// break ; -// } -// -// // Si deve operare la modifica dei triangoli -// if ( SharedIndex.size() > 2) { -// -// int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; -// -// // --- -// if ( nProd != 1 && nProd != - 2 && nProd != 4) { -// -// VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], -// VecTriHold[t2][n2].ptCompoVert[nCompo2]) ; -// VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], -// VecTriHold[t1][n1].ptCompoVert[nCompo1]) ; -// -// VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].Validate( true) ; -// VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].Validate( true) ; -// -// bModified = true ; -// break ; -// } -// } -// } -// -// if ( bModified) -// break ; -// } -// } -// } -// } -// } -// } -// } -// } -// } -// -// return true ; -//} - -//---------------------------------------------------------------------------- -bool -VolZmap::FlipEdges( TriaMatrix& VecTriHold) const -{ - double dSqEps = EPS_SMALL * EPS_SMALL ; - - // Ciclo sui blocchi - for ( size_t t1 = 0 ; t1 < m_nNumBlock ; ++ t1) { - - // Determino i,j,k del primo blocco - int nIJK1[3] ; - GetBlockIJKFromN( int( t1), nIJK1) ; - - for ( size_t t2 = t1 ; t2 < m_nNumBlock ; ++ t2) { - - // Determino i,j,k del secondo blocco - int nIJK2[3] ; - GetBlockIJKFromN( int( t2), nIJK2) ; - - // Se i blocchi sono adiacenti o coincidenti proseguo - if ( ( nIJK2[0] >= nIJK1[0] - 1 && nIJK2[0] <= nIJK1[0] + 1) || - ( nIJK2[1] >= nIJK1[1] - 1 && nIJK2[1] <= nIJK1[1] + 1) || - ( nIJK2[2] >= nIJK1[2] - 1 && nIJK2[2] <= nIJK1[2] + 1)) { - - // Numero di voxel in cui si presentano sharp feature - int nVoxelNum1 = int( VecTriHold[t1].size()) ; - int nVoxelNum2 = int( VecTriHold[t2].size()) ; - - // Determino estremi del ciclo sui voxel esterno - int nSt1 = 0 ; - int nEn1 = nVoxelNum1 - ( t1 == t2 ? 1 : 0) ; - - // Ciclo su tali voxel - for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { - - // Determino estremi del ciclo sui voxel interno - int nSt2 = ( t1 == t2 ? n1 + 1 : 0) ; - int nEn2 = nVoxelNum2 ; - - for ( int n2 = nSt2 ; n2 < nVoxelNum2 ; ++ n2) { - - // Se i voxel sono adiacenti proseguo - if ( VecTriHold[t2][n2].i == VecTriHold[t1][n1].i - 1 || - VecTriHold[t2][n2].i == VecTriHold[t1][n1].i + 1 || - VecTriHold[t2][n2].j == VecTriHold[t1][n1].j - 1 || - VecTriHold[t2][n2].j == VecTriHold[t1][n1].j + 1 || - VecTriHold[t2][n2].k == VecTriHold[t1][n1].k - 1 || - VecTriHold[t2][n2].k == VecTriHold[t1][n1].k + 1) { - - // Numero delle componenti connesse nei due voxel - int nNumCompo1 = int( VecTriHold[t1][n1].ptCompoVert.size()) ; - int nNumCompo2 = int( VecTriHold[t2][n2].ptCompoVert.size()) ; - - int nCompo1 = 0 ; - - // Ciclo sulle componenti - for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { - - int nCompo2 = ( t1 == t2 && n1 == n2 ? nCompo1 + 1 : 0) ; - - for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { - - // Numero di triangoli per le componenti connesse - int nTriNum1 = int( VecTriHold[t1][n1].vCompoTria[nCompo1].size()) ; - int nTriNum2 = int( VecTriHold[t2][n2].vCompoTria[nCompo2].size()) ; - - for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) { - - bool bModified = false ; - - for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) { - - INTVECTOR SharedIndex ; - - for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) { - - for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) { - - Point3d ptP1 = VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; - Point3d ptP2 = VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; - - if ( SqDist( ptP1, ptP2) < dSqEps) { - - Point3d ptVert1 = VecTriHold[t1][n1].ptCompoVert[nCompo1] ; - Point3d ptVert2 = VecTriHold[t2][n2].ptCompoVert[nCompo2] ; - - if ( ! ( AreSamePointApprox( ptP1, ptVert1) || - AreSamePointApprox( ptP2, ptVert2))) { - - SharedIndex.emplace_back( nVert1) ; - SharedIndex.emplace_back( nVert2) ; - } - } - - if ( SharedIndex.size() > 2) - break ; - } - - if ( SharedIndex.size() > 2) - break ; - } - - // Si deve operare la modifica dei triangoli - if ( SharedIndex.size() > 2) { - - int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; - - // --- - if ( nProd != 1 && nProd != - 2 && nProd != 4) { - - VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], - VecTriHold[t2][n2].ptCompoVert[nCompo2]) ; - VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], - VecTriHold[t1][n1].ptCompoVert[nCompo1]) ; - - VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].Validate( true) ; - VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].Validate( true) ; - if ( t1 != t2) { - int bau = 1 ; - } - bModified = true ; - break ; - } - } - } - - if ( bModified) - break ; - } - } - } - } - } - } - } - } - } - - return true ; -} - -//---------------------------------------------------------------------------- -//bool -//VolZmap::FlipEdges( std::vector& VecTriHold) const -//{ -// double dSqEps = EPS_SMALL * EPS_SMALL ; -// size_t nEnd = VecTriHold.size() ; //m_nNumBlock ; -// // Ciclo sui blocchi -// for ( size_t t = 0 ; t < nEnd ; ++ t) { -// -// // Numero di voxel in cui si presentano sharp feature -// int nVoxelNum = int( VecTriHold[t].size()) ; -// -// // Determino estremi del ciclo sui voxel esterno -// int nSt1 = 0 ; -// int nEn1 = nVoxelNum - 1 ; -// -// // Ciclo su tali voxel -// for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { -// -// // Determino estremi del ciclo sui voxel interno -// int nSt2 = nSt1 + 1 ; -// int nEn2 = nVoxelNum ; -// -// for ( int n2 = nSt2 ; n2 < nEn2 ; ++ n2) { -// -// // Se i voxel sono adiacenti proseguo -// if ( VecTriHold[t][n2].i == VecTriHold[t][n1].i + 1 || -// VecTriHold[t][n2].j == VecTriHold[t][n1].j + 1 || -// VecTriHold[t][n2].k == VecTriHold[t][n1].k + 1) { -// -// // Numero delle componenti connesse nei due voxel -// int nNumCompo1 = int( VecTriHold[t][n1].ptCompoVert.size()) ; -// int nNumCompo2 = int( VecTriHold[t][n2].ptCompoVert.size()) ; -// -// int nCompo1 = 0 ; -// -// // Ciclo sulle componenti -// for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { -// -// int nCompo2 = 0 ; -// -// for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { -// -// // Numero di triangoli per le componenti connesse -// int nTriNum1 = int( VecTriHold[t][n1].vCompoTria[nCompo1].size()) ; -// int nTriNum2 = int( VecTriHold[t][n2].vCompoTria[nCompo2].size()) ; -// -// for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) { -// -// bool bModified = false ; -// -// for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) { -// -// INTVECTOR SharedIndex ; -// -// for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) { -// -// for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) { -// -// Point3d ptP1 = VecTriHold[t][n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; -// Point3d ptP2 = VecTriHold[t][n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; -// -// if ( SqDist( ptP1, ptP2) < dSqEps) { -// -// Point3d ptVert1 = VecTriHold[t][n1].ptCompoVert[nCompo1] ; -// Point3d ptVert2 = VecTriHold[t][n2].ptCompoVert[nCompo2] ; -// -// if ( ! ( AreSamePointApprox( ptP1, ptVert1) || -// AreSamePointApprox( ptP2, ptVert2))) { -// -// SharedIndex.emplace_back( nVert1) ; -// SharedIndex.emplace_back( nVert2) ; -// } -// } -// -// if ( SharedIndex.size() > 2) -// break ; -// } -// -// if ( SharedIndex.size() > 2) -// break ; -// } -// -// // Si deve operare la modifica dei triangoli -// if ( SharedIndex.size() > 2) { -// -// int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; -// -// // --- -// if ( nProd != 1 && nProd != - 2 && nProd != 4) { -// -// Triangle3d TriVox1 = VecTriHold[t][n1].vCompoTria[nCompo1][nTri1] ; -// Triangle3d TriVox2 = VecTriHold[t][n2].vCompoTria[nCompo2][nTri2] ; -// -// TriVox1.SetP( SharedIndex[0], VecTriHold[t][n2].ptCompoVert[nCompo2]) ; -// TriVox2.SetP( SharedIndex[3], VecTriHold[t][n1].ptCompoVert[nCompo1]) ; -// -// TriVox1.Validate( true) ; -// TriVox2.Validate( true) ; -// // questo dentro il commento è un'assurdità -// /*bool bContinue = true ; -// -// for ( int nControl = 0 ; nControl < nTriNum1 ; ++ nControl) { -// -// Vector3d vtControl = VecTriHold[t][n1].vCompoTria[nCompo1][nControl].GetN() ; -// Vector3d vtVox = TriVox1.GetN() ; -// -// if ( vtControl * vtVox < - 0.9) { -// bContinue = false ; -// break ; -// } -// } -// -// if ( bContinue) { -// -// for ( int nControl = 0 ; nControl < nTriNum2 ; ++ nControl) { -// -// Vector3d vtControl = VecTriHold[t][n2].vCompoTria[nCompo2][nControl].GetN() ; -// Vector3d vtVox = TriVox2.GetN() ; -// -// if ( vtControl * vtVox < - 0.9) { -// bContinue = false ; -// break ; -// } -// } -// } -// -// if ( bContinue) {*/ -// -// VecTriHold[t][n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], -// VecTriHold[t][n2].ptCompoVert[nCompo2]) ; -// VecTriHold[t][n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], -// VecTriHold[t][n1].ptCompoVert[nCompo1]) ; -// -// VecTriHold[t][n1].vCompoTria[nCompo1][nTri1].Validate( true) ; -// VecTriHold[t][n2].vCompoTria[nCompo2][nTri2].Validate( true) ; -// -// bModified = true ; -// break ; -// /*}*/ -// } -// } -// } -// -// if ( bModified) -// break ; -// } -// } -// } -// } -// } -// } -// } -// return true ; -//} - //---------------------------------------------------------------------------- bool VolZmap::FlipEdgesII( TriHolder& TriHold) const @@ -2286,7 +1759,7 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; - // --- + if ( nProd != 1 && nProd != - 2 && nProd != 4) { Triangle3d TriVox1 = TriHold[n1].vCompoTria[nCompo1][nTri1] ; @@ -2297,47 +1770,17 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const TriVox1.Validate( true) ; TriVox2.Validate( true) ; - // questo dentro il commento è un'assurdità - /*bool bContinue = true ; - for ( int nControl = 0 ; nControl < nTriNum1 ; ++ nControl) { - - Vector3d vtControl = VecTriHold[t][n1].vCompoTria[nCompo1][nControl].GetN() ; - Vector3d vtVox = TriVox1.GetN() ; - - if ( vtControl * vtVox < - 0.9) { - bContinue = false ; - break ; - } - } - - if ( bContinue) { - - for ( int nControl = 0 ; nControl < nTriNum2 ; ++ nControl) { - - Vector3d vtControl = VecTriHold[t][n2].vCompoTria[nCompo2][nControl].GetN() ; - Vector3d vtVox = TriVox2.GetN() ; - - if ( vtControl * vtVox < - 0.9) { - bContinue = false ; - break ; - } - } - } - - if ( bContinue) {*/ - - TriHold[n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], + TriHold[n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], TriHold[n2].ptCompoVert[nCompo2]) ; - TriHold[n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], + TriHold[n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], TriHold[n1].ptCompoVert[nCompo1]) ; - TriHold[n1].vCompoTria[nCompo1][nTri1].Validate( true) ; - TriHold[n2].vCompoTria[nCompo2][nTri2].Validate( true) ; + TriHold[n1].vCompoTria[nCompo1][nTri1].Validate( true) ; + TriHold[n2].vCompoTria[nCompo2][nTri2].Validate( true) ; - bModified = true ; - break ; - /*}*/ + bModified = true ; + break ; } } } @@ -2486,243 +1929,6 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const { return true ; } -//---------------------------------------------------------------------------- -bool -VolZmap::FlipEdgesIB( TriHolder& InnerVox, size_t nBlock) const -{ - double dSqEps = EPS_SMALL * EPS_SMALL ; - - size_t nNumI = InnerVox.size() ; - size_t nNumB = m_InterBlockTria[nBlock].size() ; - - for ( size_t tB = 0 ; tB < nNumB ; ++ tB) { - for ( size_t tI = 0 ; tI < nNumI ; ++ tI) { - - if ( abs( InnerVox[tI].i - m_InterBlockTria[nBlock][tB].i) <= 1 && - abs( InnerVox[tI].j - m_InterBlockTria[nBlock][tB].j) <= 1 && - abs( InnerVox[tI].k - m_InterBlockTria[nBlock][tB].k) <= 1) { - - // Numero delle componenti connessVecTre nei due voxel - int nNumCompoB = int( m_InterBlockTria[nBlock][tB].ptCompoVert.size()) ; - int nNumCompoI = int( InnerVox[tI].ptCompoVert.size()) ; - - int nCompoB = 0 ; - - // Ciclo sulle componenti - for ( ; nCompoB < nNumCompoB ; ++ nCompoB) { - - int nCompoI = 0 ; - - for ( ; nCompoI < nNumCompoI ; ++ nCompoI) { - - // Numero di triangoli per le componenti connesse - int nTriNumB = int( m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB].size()) ; - int nTriNumI = int( InnerVox[tI].vCompoTria[nCompoI].size()) ; - - for ( int nTriB = 0 ; nTriB < nTriNumB ; ++ nTriB) { - - bool bModified = false ; - - for ( int nTriI = 0 ; nTriI < nTriNumI ; ++ nTriI) { - - INTVECTOR SharedIndex ; - - for ( int nVertB = 0 ; nVertB < 3 ; ++ nVertB) { - - for ( int nVertI = 0 ; nVertI < 3 ; ++ nVertI) { - - Point3d ptPB = m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB].GetP( nVertB) ; - Point3d ptPI = InnerVox[tI].vCompoTria[nCompoI][nTriI].GetP( nVertI) ; - - if ( SqDist( ptPB, ptPI) < dSqEps) { - - Point3d ptVertB = m_InterBlockTria[nBlock][tB].ptCompoVert[nCompoB] ; - Point3d ptVertI = InnerVox[tI].ptCompoVert[nCompoI] ; - - if ( ! ( AreSamePointApprox( ptPB, ptVertB) || - AreSamePointApprox( ptPI, ptVertI))) { - - SharedIndex.emplace_back( nVertB) ; - SharedIndex.emplace_back( nVertI) ; - } - } - - if ( SharedIndex.size() > 2) - break ; - } - - if ( SharedIndex.size() > 2) - break ; - } - - // Si deve operare la modifica dei triangoli - if ( SharedIndex.size() > 2) { - - int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; - - // --- - if ( nProd != 1 && nProd != - 2 && nProd != 4) { - - Triangle3d TriVoxB = m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB] ; - Triangle3d TriVoxI = InnerVox[tI].vCompoTria[nCompoI][nTriI] ; - - TriVoxB.SetP( SharedIndex[0], InnerVox[tI].ptCompoVert[nCompoI]) ; - TriVoxI.SetP( SharedIndex[3], m_InterBlockTria[nBlock][tB].ptCompoVert[nCompoB]) ; - - TriVoxB.Validate( true) ; - TriVoxI.Validate( true) ; - - - m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB].SetP( SharedIndex[0], - InnerVox[tI].ptCompoVert[nCompoI]) ; - InnerVox[tI].vCompoTria[nCompoI][nTriI].SetP( SharedIndex[3], - m_InterBlockTria[nBlock][tB].ptCompoVert[nCompoB]) ; - - m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB].Validate( true) ; - InnerVox[tI].vCompoTria[nCompoI][nTriI].Validate( true) ; - - bModified = true ; - break ; - } - } - } - - if ( bModified) - break ; - } - } - } - - } - } - } - return true ; - -} - -//---------------------------------------------------------------------------- -bool -VolZmap::FlipEdgesBB( size_t nCurBlock, size_t nAdjBlock) const -{ - double dSqEps = EPS_SMALL * EPS_SMALL ; - - return true ; - -} - -//---------------------------------------------------------------------------- -bool -VolZmap::FlipEdgesIB( TriHolder& InnerVox, TriHolder& BoundaryVox) const -{ - double dSqEps = EPS_SMALL * EPS_SMALL ; - - size_t nNumI = InnerVox.size() ; - size_t nNumB = BoundaryVox.size() ; - - for ( size_t tB = 0 ; tB < nNumB ; ++ tB) { - for ( size_t tI = 0 ; tI < nNumI ; ++ tI) { - - if ( abs( InnerVox[tI].i - BoundaryVox[tB].i) <= 1 && - abs( InnerVox[tI].j - BoundaryVox[tB].j) <= 1 && - abs( InnerVox[tI].k - BoundaryVox[tB].k) <= 1) { - - // Numero delle componenti connessVecTre nei due voxel - int nNumCompoB = int( BoundaryVox[tB].ptCompoVert.size()) ; - int nNumCompoI = int( InnerVox[tI].ptCompoVert.size()) ; - - int nCompoB = 0 ; - - // Ciclo sulle componenti - for ( ; nCompoB < nNumCompoB ; ++ nCompoB) { - - int nCompoI = 0 ; - - for ( ; nCompoI < nNumCompoI ; ++ nCompoI) { - - // Numero di triangoli per le componenti connesse - int nTriNumB = int( BoundaryVox[tB].vCompoTria[nCompoB].size()) ; - int nTriNumI = int( InnerVox[tI].vCompoTria[nCompoI].size()) ; - - for ( int nTriB = 0 ; nTriB < nTriNumB ; ++ nTriB) { - - bool bModified = false ; - - for ( int nTriI = 0 ; nTriI < nTriNumI ; ++ nTriI) { - - INTVECTOR SharedIndex ; - - for ( int nVertB = 0 ; nVertB < 3 ; ++ nVertB) { - - for ( int nVertI = 0 ; nVertI < 3 ; ++ nVertI) { - - Point3d ptPB = BoundaryVox[tB].vCompoTria[nCompoB][nTriB].GetP( nVertB) ; - Point3d ptPI = InnerVox[tI].vCompoTria[nCompoI][nTriI].GetP( nVertI) ; - - if ( SqDist( ptPB, ptPI) < dSqEps) { - - Point3d ptVertB = BoundaryVox[tB].ptCompoVert[nCompoB] ; - Point3d ptVertI = InnerVox[tI].ptCompoVert[nCompoI] ; - - if ( ! ( AreSamePointApprox( ptPB, ptVertB) || - AreSamePointApprox( ptPI, ptVertI))) { - - SharedIndex.emplace_back( nVertB) ; - SharedIndex.emplace_back( nVertI) ; - } - } - - if ( SharedIndex.size() > 2) - break ; - } - - if ( SharedIndex.size() > 2) - break ; - } - - // Si deve operare la modifica dei triangoli - if ( SharedIndex.size() > 2) { - - int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; - - // --- - if ( nProd != 1 && nProd != - 2 && nProd != 4) { - - Triangle3d TriVoxB = BoundaryVox[tB].vCompoTria[nCompoB][nTriB] ; - Triangle3d TriVoxI = InnerVox[tI].vCompoTria[nCompoI][nTriI] ; - - TriVoxB.SetP( SharedIndex[0], InnerVox[tI].ptCompoVert[nCompoI]) ; - TriVoxI.SetP( SharedIndex[3], BoundaryVox[tB].ptCompoVert[nCompoB]) ; - - TriVoxB.Validate( true) ; - TriVoxI.Validate( true) ; - - - BoundaryVox[tB].vCompoTria[nCompoB][nTriB].SetP( SharedIndex[0], - InnerVox[tI].ptCompoVert[nCompoI]) ; - InnerVox[tI].vCompoTria[nCompoI][nTriI].SetP( SharedIndex[3], - BoundaryVox[tB].ptCompoVert[nCompoB]) ; - - BoundaryVox[tB].vCompoTria[nCompoB][nTriB].Validate( true) ; - InnerVox[tI].vCompoTria[nCompoI][nTriI].Validate( true) ; - - bModified = true ; - break ; - } - } - } - - if ( bModified) - break ; - } - } - } - - } - } - } - return true ; -} - //---------------------------------------------------------------------------- bool VolZmap::IsThereMat( int nI, int nJ, int nK) const diff --git a/VolTriZmapVolume.cpp b/VolTriZmapVolume.cpp index e223569..7e29999 100644 --- a/VolTriZmapVolume.cpp +++ b/VolTriZmapVolume.cpp @@ -767,7 +767,7 @@ VolZmap::CylBall_ZDrilling( unsigned int nGrid, const Point3d & ptS, const Point double dSqLen = vtC.SqLen() ; // Se il punto si trova dentro il cerchio taglio - if ( dSqLen < dSqRad) + if ( dSqLen < dSqRad - 2 * m_dRadius * EPS_SMALL) // utensile cilindrico if ( m_nToolType == CylindricalMill) SubtractIntervals( nGrid, i, j, dMinStemZ, dMaxStemZ, Z_AX, - Z_AX) ; @@ -2812,261 +2812,6 @@ VolZmap::Conus_Milling( unsigned int nGrid, const Point3d & ptS, const Point3d & } // ---------- Utensile generico ---------------------------------------------- - /* -//---------------------------------------------------------------------------- -bool -VolZmap::GenTool_Drilling( unsigned int nGrid, const Point3d & ptS, const Point3d & ptE, const Vector3d & vtToolDir) -{ - // Posizioni iniziale e finale dell'utensile - Point3d ptI = ptS ; - Point3d ptF = ptE ; - - // vettore movimento - Vector3d vtMove = ptE - ptS ; - - // Settaggio profilo - CurveComposite* pToolProfile ; - if ( m_ToolArcLineApprox.GetCurveCount() == 0) - // Se l'utensile non è stato approssimato uso l'originale - pToolProfile = & m_ToolOutline ; - else - // altrimenti usi l'approssimazione - pToolProfile = &m_ToolArcLineApprox ; - - // Ciclo sulle curve - const ICurve* pCurve = pToolProfile->GetFirstCurve() ; - while ( pCurve != nullptr) { - - double dHeight ; - - int nCurveType = pCurve -> GetType() ; - - // Caso di semento - if ( nCurveType == CRV_LINE) { - - Point3d ptStart, ptEnd ; - - pCurve -> GetStartPoint( ptStart) ; - pCurve -> GetEndPoint( ptEnd) ; - - if ( abs( ptStart.y - ptEnd.y) > EPS_SMALL) { - - dHeight = abs( ptStart.y - ptEnd.y) ; - - // Il componente è un cilindro - if ( abs( ptStart.x - ptEnd.x) < EPS_SMALL) { - - double dRadius = ptStart.x ; - - // CompCyl_Drilling( nGrid, ptI, ptF, vtToolDir, dHeight, dRadius) ; - } - // Il componente è un cono con vettore equiverso a quello dell'utensile - else if ( ptStart.x > ptEnd.x) { - - double dMaxRad = ptStart.x ; - double dMinRad = ptEnd.x ; - - CompConus_Drilling( nGrid, ptI, ptF, vtToolDir, dHeight, dMaxRad, dMinRad) ; - } - // Il componente è un cono con vettore opposto a quello dell'utensile - else if ( ptStart.x < ptEnd.x) { - - double dMaxRad = ptEnd.x ; - double dMinRad = ptStart.x ; - - Point3d ptIn = ptI - vtToolDir * dHeight ; - Point3d ptFn = ptIn + vtMove ; - - CompConus_Drilling( nGrid, ptIn, ptFn, - vtToolDir, dHeight, dMaxRad, dMinRad) ; - } - } - else - dHeight = 0 ; - } - - // Caso arco - else if ( nCurveType == CRV_ARC) { - - // Centro e Punti iniziale e finale del cerchio - Point3d ptStart, ptEnd, ptO ; - - pCurve -> GetStartPoint( ptStart) ; - pCurve -> GetEndPoint( ptEnd) ; - pCurve -> GetCenterPoint( ptO) ; - - // Determino il raggio - Vector3d vtStRad = ptStart - ptO ; - Vector3d vtEnRad = ptEnd - ptO ; - - double dRadius = 0.5 * ( vtStRad.LenXY() + vtEnRad.LenXY()) ; - - // Determino le posizioni iniziale e finale del centrodella sfera - Point3d ptOSt = ptI - vtToolDir * ( ptStart.y - ptO.y) ; - Point3d ptOEn = ptOSt + vtMove ; - - // Eseguo l'asportazione del materiale - CompBall_Milling( nGrid, ptOSt, ptOEn, dRadius) ; - - - // aggiorno l'altezza - dHeight = abs( ptStart.y - ptEnd.y) ; - } - - // Determino le posizioni iniziale e finale del componente successivo - ptI = ptI - vtToolDir * dHeight ; - ptF = ptI + vtMove ; - - // Aggiorno il puntatore - pCurve = pToolProfile->GetNextCurve() ; - } - - return true ; -}*/ - /* -//---------------------------------------------------------------------------- -bool -VolZmap::GenTool_Drilling( unsigned int nGrid, const Point3d & ptS, const Point3d & ptE, const Vector3d & vtToolDir) -{ - // Descrizione geometrica del moto - Point3d ptI = ptS ; - Point3d ptF = ptE ; - Vector3d vtMove = ptE - ptS ; - - // Settaggio profilo - CurveComposite* pToolProfile ; - - if ( m_ToolArcLineApprox.GetCurveCount() == 0) - // Se l'utensile non è stato approssimato uso l'originale - pToolProfile = & m_ToolOutline ; - else - // altrimenti uso l'approssimazione - pToolProfile = &m_ToolArcLineApprox ; - - // Ciclo sulle curve - const ICurve* pCurve = pToolProfile->GetFirstCurve() ; - const ICurve* pFirstConst = pCurve ; - - while ( pCurve != nullptr) { - - double dHeight ; - - int nCurveType = pCurve -> GetType() ; - - // Caso di segmento - if ( nCurveType == CRV_LINE) { - - Point3d ptStart, ptEnd ; - - pCurve -> GetStartPoint( ptStart) ; - pCurve -> GetEndPoint( ptEnd) ; - - if ( abs( ptStart.y - ptEnd.y) > EPS_SMALL) { - - dHeight = abs( ptStart.y - ptEnd.y) ; - - bool bTapB, bTapT ; - - CurveComposite* pToolPrev = pToolProfile ; - CurveComposite* ptToolNext = pToolProfile ; - const ICurve* pPrev = pToolPrev -> GetPrevCurve() ; - const ICurve* pNext = ptToolNext -> GetNextCurve() ; - - // Dettagli curva precedente - int nPrevType = pPrev -> GetType() ; - Point3d ptPrevLast ; - pPrev -> GetEndPoint( ptPrevLast) ; - - // Setto le variabili per il tappo superiore - if ( pPrev == pFirstConst || - nPrevType == CRV_ARC || - ptPrevLast.x < ptStart.x) - - bTapB = false ; - else - bTapB = true ; - - // Dettagli curva successiva - int nNextType = pNext -> GetType() ; - Point3d ptNextFirst ; - pNext -> GetEndPoint( ptNextFirst) ; - - if ( pNext == nullptr || - nNextType == CRV_ARC || - ptNextFirst.x < ptEnd.x) - - bTapT = false ; - else - bTapT = true ; - - // Il componente è un cilindro - if ( abs( ptStart.x - ptEnd.x) < EPS_SMALL) { - - double dRadius = ptStart.x ; - - CompCyl_Drilling( nGrid, ptI, ptF, vtToolDir, dHeight, dRadius, bTapB, bTapT) ; - } - // Il componente è un cono con vettore equiverso a quello dell'utensile - else if ( ptStart.x > ptEnd.x) { - - double dMaxRad = ptStart.x ; - double dMinRad = ptEnd.x ; - - CompConus_Drilling( nGrid, ptI, ptF, vtToolDir, dHeight, dMaxRad, dMinRad, bTapB, bTapT) ; - } - // Il componente è un cono con vettore opposto a quello dell'utensile - else if ( ptStart.x < ptEnd.x) { - - double dMaxRad = ptEnd.x ; - double dMinRad = ptStart.x ; - - Point3d ptIn = ptI - vtToolDir * dHeight ; - Point3d ptFn = ptIn + vtMove ; - - CompConus_Drilling( nGrid, ptIn, ptFn, - vtToolDir, dHeight, dMaxRad, dMinRad, bTapT, bTapB) ; - } - } - else - dHeight = 0 ; - } - - // Caso arco - else if ( nCurveType == CRV_ARC) { - - // Centro e Punti iniziale e finale del cerchio - Point3d ptStart, ptEnd, ptO ; - - pCurve -> GetStartPoint( ptStart) ; - pCurve -> GetEndPoint( ptEnd) ; - pCurve -> GetCenterPoint( ptO) ; - - // Determino il raggio - Vector3d vtStRad = ptStart - ptO ; - Vector3d vtEnRad = ptEnd - ptO ; - - double dRadius = 0.5 * ( vtStRad.LenXY() + vtEnRad.LenXY()) ; - - // Determino le posizioni iniziale e finale del centro della sfera - Point3d ptOSt = ptI - vtToolDir * ( ptStart.y - ptO.y) ; - Point3d ptOEn = ptOSt + vtMove ; - - // Eseguo l'asportazione del materiale - CompBall_Milling( nGrid, ptOSt, ptOEn, dRadius) ; - - - // aggiorno l'altezza - dHeight = abs( ptStart.y - ptEnd.y) ; - } - - // Determino le posizioni iniziale e finale del componente successivo - ptI = ptI - vtToolDir * dHeight ; - ptF = ptI + vtMove ; - - // Aggiorno il puntatore - pCurve = pToolProfile->GetNextCurve() ; - } - - return true ; -}*/ //---------------------------------------------------------------------------- bool @@ -3229,116 +2974,6 @@ VolZmap::GenTool_Drilling( unsigned int nGrid, const Point3d & ptS, const Point3 return true ; } -/* -//---------------------------------------------------------------------------- -bool -VolZmap::GenTool_Milling( unsigned int nGrid, const Point3d & ptS, const Point3d & ptE, const Vector3d & vtToolDir) -{ - // Posizioni iniziale e finale dell'utensile - Point3d ptI = ptS ; - Point3d ptF = ptE ; - - // vettore movimento - Vector3d vtMove = ptE - ptS ; - - // Settaggio profilo - CurveComposite* pToolProfile ; - if ( m_ToolArcLineApprox.GetCurveCount() == 0) - // Se l'utensile non è stato approssimato uso l'originale - pToolProfile = &m_ToolOutline ; - else - // altrimenti usi l'approssimazione - pToolProfile = &m_ToolArcLineApprox ; - - // Ciclo sulle curve - const ICurve* pCurve = pToolProfile->GetFirstCurve() ; - while ( pCurve != nullptr) { - - double dHeight ; - - int nCurveType = pCurve -> GetType() ; - - // Caso di semento - if ( nCurveType == CRV_LINE) { - - Point3d ptStart, ptEnd ; - - pCurve -> GetStartPoint( ptStart) ; - pCurve -> GetEndPoint( ptEnd) ; - - if ( abs( ptStart.y - ptEnd.y) > EPS_SMALL) { - - dHeight = abs( ptStart.y - ptEnd.y) ; - - // Il componente è un cilindro - if ( abs( ptStart.x - ptEnd.x) < EPS_SMALL) { - - double dRadius = ptStart.x ; - - CompCyl_Milling( nGrid, ptI, ptF, vtToolDir, dHeight, dRadius) ; - } - // Il componente è un cono con vettore equiverso a quello dell'utensile - else if ( ptStart.x > ptEnd.x) { - - double dMaxRad = ptStart.x ; - double dMinRad = ptEnd.x ; - - CompConus_Milling( nGrid, ptI, ptF, vtToolDir, dHeight, dMaxRad, dMinRad) ; - } - // Il componente è un cono con vettore opposto a quello dell'utensile - else if ( ptStart.x < ptEnd.x) { - - double dMaxRad = ptEnd.x ; - double dMinRad = ptStart.x ; - - Point3d ptIn = ptI - vtToolDir * dHeight ; - Point3d ptFn = ptIn + vtMove ; - - CompConus_Milling( nGrid, ptIn, ptFn, - vtToolDir, dHeight, dMaxRad, dMinRad) ; - } - } - else - dHeight = 0 ; - } - // Caso arco - else if ( nCurveType == CRV_ARC) { - - // Centro e Punti iniziale e finale del cerchio - Point3d ptStart, ptEnd, ptO ; - - pCurve -> GetStartPoint( ptStart) ; - pCurve -> GetEndPoint( ptEnd) ; - pCurve -> GetCenterPoint( ptO) ; - - // Determino il raggio - Vector3d vtStRad = ptStart - ptO ; - Vector3d vtEnRad = ptEnd - ptO ; - - double dRadius = 0.5 * ( vtStRad.LenXY() + vtEnRad.LenXY()) ; - - // Determino le posizioni iniziale e finale del centrodella sfera - Point3d ptOSt = ptI - vtToolDir * ( ptStart.y - ptO.y) ; - Point3d ptOEn = ptOSt + vtMove ; - - // Eseguo l'asportazione del materiale - CompBall_Milling( nGrid, ptOSt, ptOEn, dRadius) ; - - - // aggiorno l'altezza - dHeight = abs( ptStart.y - ptEnd.y) ; - } - - // Determino le posizioni iniziale e finale del componente successivo - ptI = ptI - vtToolDir * dHeight ; - ptF = ptI + vtMove ; - - // Aggiorno il puntatore - pCurve = pToolProfile->GetNextCurve() ; - } - - return true ; -}*/ - //---------------------------------------------------------------------------- bool VolZmap::GenTool_Milling( unsigned int nGrid, const Point3d & ptS, const Point3d & ptE, const Vector3d & vtToolDir) @@ -3784,163 +3419,6 @@ VolZmap::CompCyl_ZMilling( unsigned int nGrid, const Point3d & ptS, const Point3 return true ; } -/* -//---------------------------------------------------------------------------- -bool -VolZmap::CompConus_ZMilling( unsigned int nGrid, const Point3d & ptS, const Point3d & ptE, const Vector3d & vtToolDir, double dHei, double dMaxRad, double dMinRad) -{ - double dMin, dMax, dPLim, dMLim ; - unsigned int nStartI, nStartJ, nEndI, nEndJ ; - - bool Control = BBoxComponent( nGrid, ptS, ptE, vtToolDir, vtToolDir, nStartI, nStartJ, nEndI, nEndJ, dMaxRad, dMinRad, dHei) ; - - if ( ! Control) - return true ; - - Point3d ptI = ( vtToolDir * ( ptE - ptS) > 0 ? ptS : ptE) ; - Point3d ptF = ( vtToolDir * ( ptE - ptS) > 0 ? ptE : ptS) ; - - Point3d ptIxy( ptI.x, ptI.y, 0) ; - Point3d ptFxy( ptF.x, ptF.y, 0) ; - - Vector3d vtMove = ptF - ptI ; double dLen = vtMove.Len() ; - Vector3d vtMLong = ( vtMove * vtToolDir) * vtToolDir ; double dLLong = vtMLong.Len() ; - Vector3d vtMOrt = vtMove - vtMLong ; double dLOrt = vtMOrt.Len() ; - - Vector3d vtV1 = vtToolDir ; - Vector3d vtV2 = vtMOrt ; vtV2.Normalize() ; - Vector3d vtV3 = vtV1 ^ vtV2 ; - - double dZI = ptI.z ; - double dZTI = ptI.z - vtV1.z * dHei ; - double dDeltaZ = ptF.z - ptI.z ; - double dDeltaR = dMaxRad - dMinRad ; - - double dTan = dDeltaR / dHei ; - double dRatio = ( vtMove * vtV1) / ( vtMove * vtV2) ; - - double dCos = dTan * dRatio ; - double dSin = ( 1 - dCos * dCos > 0 ? sqrt( 1 - dCos * dCos) : 0) ; - - double dDen = sqrt( 1 + dTan * dTan) ; - - Point3d ptV = ptI - vtV1 * ( dHei * dMaxRad / dDeltaR) ; - - Vector3d vtNs = - ( dTan / dDen) * vtV1 + ( dCos / dDen) * vtV2 + ( dSin / dDen) * vtV3 ; - Vector3d vtNd = - ( dTan / dDen) * vtV1 + ( dCos / dDen) * vtV2 - ( dSin / dDen) * vtV3 ; - Vector3d vtR0 = ptV - ORIG ; - - double dDots = vtR0 * vtNs ; - double dDotd = vtR0 * vtNd ; - - - for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) { - - for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) { - - double dX = ( i + 0.5) * m_dStep ; double dY = ( j + 0.5) * m_dStep ; - - Point3d ptC( dX, dY, 0) ; - - Vector3d vtCI = ptC - ptIxy ; double dSqDI = vtCI.SqLenXY() ; - Vector3d vtCF = ptC - ptFxy ; double dSqDF = vtCF.SqLenXY() ; - - double dIDO = vtCI * vtV3 ; - double dIDL = vtCI * vtV2 ; - double dIVarCos = dIDL / sqrt( dSqDI) ; - - double dFDL = vtCF * vtV2 ; - double dFVarCos = dFDL / sqrt( dSqDF) ; - - if ( dSqDI < dMaxRad * dMaxRad || dSqDF < dMaxRad * dMaxRad || - (abs( dIDO) < dMaxRad && dIDL > 0 && dIDL < dLOrt)) { - - // Caso dTan > 1 / dRatio - if ( dRatio > 1 / dTan) { - - // Limiti nella direzione positiva di vtV1 - if ( dSqDF < dMaxRad * dMaxRad) - - dPLim = dZI + dDeltaZ ; - - else - - dPLim = dZI + ( dIDL + sqrt( dMaxRad * dMaxRad - dIDO * dIDO)) * dDeltaZ / dLOrt ; - - // Limiti nella direzione negativa di vtV1 - if ( dSqDI < dMinRad * dMinRad) - - dMLim = dZTI ; - - else if ( dSqDI < dMaxRad * dMaxRad) - - dMLim = dZTI + ( sqrt( dSqDI) - dMinRad) * ( dZI - dZTI) / dDeltaR ; - - else - - dMLim = dZI + ( dIDL - sqrt( dMaxRad * dMaxRad - dIDO * dIDO)) * dDeltaZ / dLOrt ; - - } - else { - - // Limiti nella direzione positiva di vtV1 - if ( dSqDF < dMaxRad * dMaxRad) - - dPLim = dZI + dDeltaZ ; - - else - - dPLim = dZI + ( dIDL + sqrt( dMaxRad * dMaxRad - dIDO * dIDO)) * dDeltaZ / dLOrt ; - - // Limiti nella direzione negativa di vtV1 - if ( dSqDI < dMinRad * dMinRad) - - dMLim = dZTI ; - - else if ( dSqDI >= dMinRad * dMinRad && dSqDI < dMaxRad * dMaxRad && dIVarCos < dCos) - - dMLim = dZTI + ( sqrt( dSqDI) - dMinRad) * ( dZI - dZTI) / dDeltaR ; - - else if ( dSqDI >= dMinRad * dMinRad && dIVarCos >= dCos && dFVarCos < dCos && abs( dIDO) < dMaxRad * dSin) { // da qui - - if ( dIDO > - dMaxRad * dSin && dIDO <= - dMinRad * dSin) - - dMLim = ( dDotd - dX * vtNd.x - dY * vtNd.y) / vtNd.z ; - - else if ( dIDO > - dMinRad * dSin && dIDO < dMinRad * dSin) - - dMLim = dZTI + ( dIDL - sqrt( dMinRad * dMinRad - dIDO * dIDO)) * dDeltaZ / dLOrt ; - - else if ( dIDO >= dMinRad * dSin && dIDO < dMaxRad * dSin) - - dMLim = ( dDots - dX * vtNs.x - dY * vtNs.y) / vtNs.z ; // a qui - } - else if ( dFVarCos >= dCos) { - - if ( dSqDF < dMinRad * dMinRad) - - dMLim = dZTI + ( dIDL - sqrt( dMinRad * dMinRad - dIDO * dIDO)) * dDeltaZ / dLOrt ; - - else - - dMLim = dZTI + dDeltaZ + ( sqrt( dSqDF) - dMinRad) * ( dZI - dZTI) / dDeltaR ; - } - else - - dMLim = dZI + ( dIDL - sqrt( dMaxRad * dMaxRad - dIDO * dIDO)) * dDeltaZ / dLOrt ; - } - - dMin = min( dPLim, dMLim) ; - dMax = max( dPLim, dMLim) ; - - SubtractIntervals( nGrid, i, j, dMin, dMax, V_NULL, V_NULL) ; - } - } - } - return true ; -} */ - - //---------------------------------------------------------------------------- bool VolZmap::CompConus_ZMilling( unsigned int nGrid, const Point3d & ptS, const Point3d & ptE, const Vector3d & vtToolDir, double dHei, double dMaxRad, double dMinRad) @@ -4231,227 +3709,6 @@ VolZmap::CompConus_ZMilling( unsigned int nGrid, const Point3d & ptS, const Poin } return true ; } -/* -//---------------------------------------------------------------------------- -bool -VolZmap::CompConus_ZMilling( unsigned int nGrid, const Point3d & ptS, const Point3d & ptE, const Vector3d & vtToolDir, double dHei, double dMaxRad, double dMinRad) -{ - unsigned int nStartI, nStartJ, nEndI, nEndJ ; - - bool Control = BBoxComponent( nGrid, ptS, ptE, vtToolDir, vtToolDir, nStartI, nStartJ, nEndI, nEndJ, dMaxRad, dMinRad, dHei) ; - - if ( ! Control) - return true ; - - Point3d ptI = ( vtToolDir * ( ptE - ptS) > 0 ? ptS : ptE) ; - Point3d ptF = ( vtToolDir * ( ptE - ptS) > 0 ? ptE : ptS) ; - - Point3d ptIT = ptI - vtToolDir * dHei ; - Point3d ptFT = ptF - vtToolDir * dHei ; - - Point3d ptIxy( ptI.x, ptI.y, 0) ; - Point3d ptFxy( ptF.x, ptF.y, 0) ; - - Vector3d vtMove = ptF - ptI ; double dLen = vtMove.Len() ; - Vector3d vtMLong = ( vtMove * vtToolDir) * vtToolDir ; double dLLong = vtMLong.Len() ; - Vector3d vtMOrt = vtMove - vtMLong ; double dLOrt = vtMOrt.Len() ; - - Vector3d vtV1 = vtToolDir ; - Vector3d vtV2 = vtMOrt ; vtV2.Normalize() ; - Vector3d vtV3 = vtV1 ^ vtV2 ; - - double dZI = ptI.z ; - double dZTI = ptI.z - vtV1.z * dHei ; - double dDeltaZ = ptF.z - ptI.z ; - double dDeltaR = dMaxRad - dMinRad ; - - double dTan = dDeltaR / dHei ; - double dRatio = ( vtMove * vtV1) / ( vtMove * vtV2) ; - - bool bCase = ( dRatio * dTan > 1 ? false : true) ; - - double dCos = dTan * dRatio ; - double dSin = ( 1 - dCos * dCos > 0 ? sqrt( 1 - dCos * dCos) : 0) ; - - double dDen = sqrt( 1 + dTan * dTan) ; - - Point3d ptV = ptI - vtV1 * ( dHei * dMaxRad / dDeltaR) ; - Point3d ptVF = ptV + vtMove ; - - Vector3d vtNs = - ( dTan / dDen) * vtV1 + ( dCos / dDen) * vtV2 + ( dSin / dDen) * vtV3 ; - Vector3d vtNd = - ( dTan / dDen) * vtV1 + ( dCos / dDen) * vtV2 - ( dSin / dDen) * vtV3 ; - Vector3d vtR0 = ptV - ORIG ; - - vtNs.Normalize() ; - vtNd.Normalize() ; - - double dDots = vtR0 * vtNs ; - double dDotd = vtR0 * vtNd ; - - double dMin, dMax, dPLim, dMLim ; - Vector3d vtMin, vtMax, vtP, vtM ; - - Vector3d vtUmv = vtMove ; - vtUmv.Normalize() ; - - Point3d ptInt ; - - for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) { - - for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) { - - double dX = ( i + 0.5) * m_dStep ; double dY = ( j + 0.5) * m_dStep ; - - Point3d ptC( dX, dY, 0) ; - - Vector3d vtCI = ptC - ptIxy ; double dSqDI = vtCI.SqLenXY() ; - Vector3d vtCF = ptC - ptFxy ; double dSqDF = vtCF.SqLenXY() ; - - double dIDO = vtCI * vtV3 ; - double dIDL = vtCI * vtV2 ; - double dIVarCos = dIDL / sqrt( dSqDI) ; - - double dFDL = vtCF * vtV2 ; - double dFVarCos = dFDL / sqrt( dSqDF) ; - - if ( dSqDI < dMaxRad * dMaxRad || dSqDF < dMaxRad * dMaxRad || - (abs( dIDO) < dMaxRad && dIDL > 0 && dIDL < dLOrt)) { - - // Caso dTan * dRatio < 1 - if ( bCase) { - - // Limiti nella direzione positiva di vtV1 - if ( dSqDF < dMaxRad * dMaxRad) { - - dPLim = dZI + dDeltaZ ; - - vtP = - vtToolDir ; - } - else { - - dPLim = dZI + ( dIDL + sqrt( dMaxRad * dMaxRad - dIDO * dIDO)) * dDeltaZ / dLOrt ; - - ptInt.Set( dX, dY, dPLim) ; - - vtP = - ( ptInt - ptI) + ( ptInt - ptI) * vtUmv * vtUmv ; - - vtP.Normalize() ; - } - - // Limiti nella direzione negativa di vtV1 - if ( dSqDI < dMinRad * dMinRad) { - - dMLim = dZTI ; - - vtM = vtToolDir ; - } - else if ( dSqDI < dMaxRad * dMaxRad) { - - dMLim = dZTI + ( sqrt( dSqDI) - dMinRad) * ( dZI - dZTI) / dDeltaR ; - - ptInt.Set( dX, dY, dMLim) ; - - Vector3d vtU = ( ptInt - ptV) - ( ptInt - ptV) * vtToolDir * vtToolDir ; - - vtU.Normalize() ; - - vtM = dDeltaR * vtToolDir - dHei * vtU ; - - vtM.Normalize() ; - } - else { - - dMLim = dZI + ( dIDL - sqrt( dMaxRad * dMaxRad - dIDO * dIDO)) * dDeltaZ / dLOrt ; - - ptInt.Set( dX, dY, dMLim) ; - - vtM = - ( ptInt - ptI) + ( ptInt - ptI) * vtUmv * vtUmv ; - - vtM.Normalize() ; - } - } - else { - - // Limiti nella direzione positiva di vtV1 - if ( dSqDF < dMaxRad * dMaxRad) { - - dPLim = dZI + dDeltaZ ; - - vtP = - vtToolDir ; - } - else { - - dPLim = dZI + ( dIDL + sqrt( dMaxRad * dMaxRad - dIDO * dIDO)) * dDeltaZ / dLOrt ; - - ptInt.Set( dX, dY, dPLim) ; - - vtP = - ( ptInt - ptI) + ( ptInt - ptI) * vtUmv * vtUmv ; - - vtP.Normalize() ; - } - - // Limiti nella direzione negativa di vtV1 - if ( dSqDI < dMinRad * dMinRad) { - - dMLim = dZTI ; - - vtM = vtToolDir ; - } - else if ( dSqDI < dMaxRad * dMaxRad) { - - dMLim = dZTI + ( sqrt( dSqDI) - dMinRad) * ( dZI - dZTI) / dDeltaR ; - - ptInt.Set( dX, dY, dMLim) ; - - Vector3d vtU = ( ptInt - ptV) - ( ptInt - ptV) * vtToolDir * vtToolDir ; - - vtU.Normalize() ; - - vtM = dDeltaR * vtToolDir - dHei * vtU ; - - vtM.Normalize() ; - } - else { - - dMLim = dZI + ( dIDL - sqrt( dMaxRad * dMaxRad - dIDO * dIDO)) * dDeltaZ / dLOrt ; - - ptInt.Set( dX, dY, dMLim) ; - - vtM = - ( ptInt - ptI) + ( ptInt - ptI) * vtUmv * vtUmv ; - - vtM.Normalize() ; - } - } - - //dMin = min( dPLim, dMLim) ; - //dMax = max( dPLim, dMLim) ; - - if ( dMLim < dPLim) { - - dMin = dMLim ; - dMax = dPLim ; - - vtMin = vtM ; - vtMax = vtP ; - } - else { - - dMin = dPLim ; - dMax = dMLim ; - - vtMin = vtP ; - vtMax = vtM ; - } - - SubtractIntervals( nGrid, i, j, dMin, dMax, vtMin, vtMax) ; - } - } - } - return true ; -}*/ - - - // Asse di simmetria con orientazione generica: FORATURA diff --git a/VolZmap.h b/VolZmap.h index 25d8158..87e5f9f 100644 --- a/VolZmap.h +++ b/VolZmap.h @@ -37,21 +37,6 @@ typedef std::vector TriHolder ; // il primo indice individua il blocco, il secondo il voxel typedef std::vector TriaMatrix ; -struct IndexStruct { - - int i, j, k ; -} ; - -typedef std::vector> IndexMatrix ; -//struct BoundaryTriangle { -// -// int i, j, k ; -// int nCC ; -// -// Point3d ptVertex ; -// -// bool bInterBlockFlipping ; -//} ; //---------------------------------------------------------------------------- class VolZmap : public IVolZmap, public IGeoObjRW @@ -138,13 +123,9 @@ class VolZmap : public IVolZmap, public IGeoObjRW const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DLIST& lstTria) const ; bool MarchingCubes( TRIA3DLIST& lstTria) const ; bool MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const ; - bool ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& lstTria, TriHolder& triHold) const ; - bool FlipEdges( TriaMatrix& VecTriHold) const ; + bool ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& lstTria, TriHolder& triHold) const ; bool FlipEdgesII( TriHolder& TriHold) const ; - bool FlipEdgesBB( TriaMatrix& InterTria) const ; - bool FlipEdgesIB( TriHolder& InnerVox, size_t nBlock) const ; - bool FlipEdgesBB( size_t nCurBlock, size_t nAdjBlock) const ; - bool FlipEdgesIB( TriHolder& InnerVox, TriHolder& BoundaryVox) const ; + bool FlipEdgesBB( TriaMatrix& InterTria) const ; bool IsThereMat( int nI, int nJ, int nK) const ; bool IsThereMat( const int nMatr[][3], int nNum, double & dHx, double & dHy, double & dHz) const ; bool IntersPos( int nVec1[], int nVec2[], Point3d & ptInt) const ;