diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index 58c2def..19b6ff2 100644 Binary files a/EgtGeomKernel.rc and b/EgtGeomKernel.rc differ diff --git a/PolyLine.cpp b/PolyLine.cpp index f94795c..b3ca97f 100644 --- a/PolyLine.cpp +++ b/PolyLine.cpp @@ -686,7 +686,7 @@ PolyLine::RemoveAlignedPoints( double dToler) dToler = max( dToler, LIN_TOL_MIN) ; // coefficiente deduzione tolleranza const double COEFF_TOL = 0.90 ; - // si analizza la distanza di un punto dal segmento che unisce precente e successivo + // si analizza la distanza di un punto dal segmento che unisce precedente e successivo // punto precedente auto precP = m_lUPoints.begin() ; // punto corrente diff --git a/VolTriZmapCalculus.cpp b/VolTriZmapCalculus.cpp index 258f074..e353f2f 100644 --- a/VolTriZmapCalculus.cpp +++ b/VolTriZmapCalculus.cpp @@ -19,6 +19,8 @@ #include "IntersLineSurfTm.h" #include "/EgtDev/Include/EgtNumUtils.h" #include "/EgtDev/Include/EGkIntersLineTria.h" +#include "/EgtDev/Include/EGkIntersLinePlane.h" +#include "/EgtDev/Include/EGkChainCurves.h" using namespace std ; @@ -1399,4 +1401,335 @@ VolZmap::GetPartVolume( int nPart, double& dVol) const return true ; } +//---------------------------------------------------------------------------- +// Il piano è espresso nel sistema locale, la funzione lo esprime nel sistema Zmap. +// I loop della flat region ottenuta dall'intersezione vengono salvati nello +// standard vector vLoop passato per riferimento. +// Se il processo va a buon fine viene restituito true, false altrimenti. +bool +VolZmap::GetPlaneIntersection( const Plane3d& plPlane, POCRVVECTOR& vpLoop) const +{ + // Verifico validità parametri + if ( ! plPlane.IsValid()) + return false ; + // Porto il piano nel sistema Zmap + Plane3d plPlaneL = plPlane ; + plPlaneL.ToLoc( m_MapFrame[0]) ; + // Vettore del piano + Vector3d vtNL = plPlaneL.GetVersN() ; + // Se il piano non interseca il bounding box del solido, abbiamo finito + if ( ! TestIntersPlaneZmapBBox( plPlaneL)) + return true ; + + // Cerco la mappa con i dexel più perpendicolari al piano e valuto + // se il versore del piano è equiverso o controverso ai dexel. + int nGrid = 0 ; + int nSign = ( vtNL.z > 0 ? 1 : - 1) ; + if ( abs( vtNL.x) > abs( vtNL.y) && abs( vtNL.x) > abs( vtNL.z)) { + nGrid = 1 ; + nSign = ( vtNL.x > 0 ? 1 : - 1) ; + } + else if ( abs( vtNL.y) > abs( vtNL.x) && abs( vtNL.y) > abs( vtNL.z)) { + nGrid = 2 ; + nSign = ( vtNL.y > 0 ? 1 : - 1) ; + } + + // Ciclo sulle celle + vector vLine ; + for ( int ni = - 1 ; ni < int( m_nNx[nGrid]) ; ++ ni) { + for ( int nj = - 1 ; nj < int( m_nNy[nGrid]) ; ++ nj) { + ProcessCell( nGrid, ni, nj, plPlaneL, vLine) ; + } + } + + // Creo i loop + ChainCurves LoopCreator ; + LoopCreator.Init( false, EPS_SMALL, int( vLine.size())) ; + // Carico le curve per concatenarle + for ( int nCv = 0 ; nCv < int( vLine.size()) ; ++ nCv) { + Point3d ptSt = vLine[nCv].GetStart() ; + Point3d ptEn = vLine[nCv].GetEnd() ; + Vector3d vtDir ; vLine[nCv].GetStartDir( vtDir) ; + LoopCreator.AddCurve( nCv + 1, ptSt, vtDir, ptEn, vtDir) ; + } + // Recupero i concatenamenti + INTVECTOR vIds ; + Point3d ptNearStart ; + while ( LoopCreator.GetChainFromNear( ptNearStart, false, vIds)) { + PtrOwner pLoop( CreateBasicCurveComposite()) ; + if ( IsNull( pLoop)) + return false ; + for ( auto i : vIds) { + // Aggiungo la linea alla curva composta. + if ( ! pLoop->AddCurve( vLine[i-1], true, EPS_SMALL)) + return false ; + } + pLoop->SetExtrusion( vtNL) ; + pLoop->ToGlob( m_MapFrame[0]) ; + // Se normali controverse, devo invertire il verso del loop. + if ( nSign < 0) + pLoop->Invert() ; + // eseguo approssimazione + PolyLine PL ; + if ( ! pLoop->ApproxWithLines( m_dStep, ANG_TOL_STD_DEG, ICurve::APL_RIGHT, PL) || + ! pLoop->Clear() || ! pLoop->FromPolyLine( PL)) + return false ; + pLoop->MergeCurves( m_dStep, ANG_TOL_STD_DEG) ; + // Inserisco la curva composita nella raccolta da ritornare + vpLoop.emplace_back( Release( pLoop)) ; + } + + return true ; +} + +//---------------------------------------------------------------------------- +// Data una cella, identificata dal numero della sua griglia e dai suoi indici i,j +// all'interno della griglia, genera gli eventuali segmenti che costituiscono la curva +// corrispondente alla frontiera della regione di piano interna la solido. +// Se la griglia o la cella non è valida, la funzione restituisce false, +// altrimenti valuta la tipologia della cella e crea gli eventuali segmenti. +// Se il processo va a buon fine la funzione restituisce true. +bool +VolZmap::ProcessCell( int nGrid, int nCellI, int nCellJ, const Plane3d& plPlane, vector& vLine) const +{ + // Se la griglia non esiste vi è un errore + if ( nGrid < 0 || nGrid > 2) + return false ; + // Se la cella non esiste vi è un errore + if ( nCellI < - 1 || nCellI >= int( m_nNx[nGrid]) || + nCellJ < - 1 || nCellJ >= int( m_nNy[nGrid])) + return false ; + + // Determino la configurazione della cella + int nIndex = CalcIndexForPlaneCells( plPlane, nGrid, nCellI, nCellJ) ; + + // Tabella segmenti + static int nLineTable[16][5] = { + { -1, -1, -1, -1, -1}, + { 0, 3, -1, -1, -1}, + { 1, 0, -1, -1, -1}, + { 1, 3, -1, -1, -1}, + { 2, 1, -1, -1, -1}, + { 0, 1, 2, 3, -1}, + { 2, 0, -1, -1, -1}, + { 2, 3, -1, -1, -1}, + { 3, 2, -1, -1, -1}, + { 0, 2, -1, -1, -1}, + { 1, 2, 3, 0, -1}, + { 1, 2, -1, -1, -1}, + { 3, 1, -1, -1, -1}, + { 0, 1, -1, -1, -1}, + { 3, 0, -1, -1, -1}, + { -1, -1, -1, -1, -1} + } ; + // Tabella dei punti medi dei segmenti: le righe rappresentano + // il numero di segmento (0, 1, 2, 3) le colonne rappresentano + // i e j. + static double dDeltaIJ[4][2] = { + { 0.5, 0}, + { 1, 0.5}, + { 0.5, 1}, + { 0, 0.5} + } ; + // Se la cella è di frontiera costruisco i segmenti della curva + if ( nIndex != 0 && nIndex != 15) { + // Ciclo su tutti gli spigoli della cella che vengono attraversati dalla curva. + for ( int nEdge = 0 ; nLineTable[nIndex][nEdge] != -1 ; nEdge += 2) { + + int nStEdge = nLineTable[nIndex][nEdge] ; + int nEnEdge = nLineTable[nIndex][nEdge + 1] ; + double nStDeltaI = dDeltaIJ[nStEdge][0] ; + double nStDeltaJ = dDeltaIJ[nStEdge][1] ; + double nEnDeltaI = dDeltaIJ[nEnEdge][0] ; + double nEnDeltaJ = dDeltaIJ[nEnEdge][1] ; + + if ( nGrid == 0) { + // Punti sulla griglia + Point3d ptGrSt( ( nCellI + nStDeltaI + 0.5) * m_dStep, ( nCellJ + nStDeltaJ + 0.5) * m_dStep, 0) ; + Point3d ptGrEn( ( nCellI + nEnDeltaI + 0.5) * m_dStep, ( nCellJ + nEnDeltaJ + 0.5) * m_dStep, 0) ; + // Corrispondenti punti sul piano + Point3d ptSt, ptEn ; + if ( IntersLinePlane( ptGrSt, Z_AX, 10, plPlane, ptSt, false) && + IntersLinePlane( ptGrEn, Z_AX, 10, plPlane, ptEn, false)) { + // Costruisco il tratto di curva + CurveLine cvLine ; + cvLine.Set( ptSt, ptEn) ; + vLine.emplace_back( cvLine) ; + } + } + + else if ( nGrid == 1) { + // Punti sulla griglia + Point3d ptGrSt( 0, ( nCellI + nStDeltaI + 0.5) * m_dStep, ( nCellJ + nStDeltaJ + 0.5) * m_dStep) ; + Point3d ptGrEn( 0, ( nCellI + nEnDeltaI + 0.5) * m_dStep, ( nCellJ + nEnDeltaJ + 0.5) * m_dStep) ; + // Corrispondenti punti sul piano + Point3d ptSt, ptEn ; + if ( IntersLinePlane( ptGrSt, X_AX, 10, plPlane, ptSt, false) && + IntersLinePlane( ptGrEn, X_AX, 10, plPlane, ptEn, false)) { + // Costruisco il tratto di curva + CurveLine cvLine ; + cvLine.Set( ptSt, ptEn) ; + vLine.emplace_back( cvLine) ; + } + } + + else { + // Punti sulla griglia + Point3d ptGrSt( ( nCellJ + nStDeltaJ + 0.5) * m_dStep, 0, ( nCellI + nStDeltaI + 0.5) * m_dStep) ; + Point3d ptGrEn( ( nCellJ + nEnDeltaJ + 0.5) * m_dStep, 0, ( nCellI + nEnDeltaI + 0.5) * m_dStep) ; + // Corrispondenti punti sul piano + Point3d ptSt, ptEn ; + if ( IntersLinePlane( ptGrSt, Y_AX, 10, plPlane, ptSt, false) && + IntersLinePlane( ptGrEn, Y_AX, 10, plPlane, ptEn, false)) { + // Costruisco il tratto di curva + CurveLine cvLine ; + cvLine.Set( ptSt, ptEn) ; + vLine.emplace_back( cvLine) ; + } + } + } + } + + return true ; +} + +//---------------------------------------------------------------------------- +// Dato un dexel, identificato dal numero della sua griglia e dai suoi indici, +// determina se un suo intervallo attraversa il piano dato. +// Il numero di griglia e gli indici del dexel devono essere validi, se +// non lo sono, la funzione restituisce false, altrimenti cerca un +// tratto di dexel che intersechi il piano, se lo trova restituisce true, +// false altrimenti. +bool +VolZmap::InOut( const Plane3d& plPlane, int nGrid, int nI, int nJ) const +{ + // Se la griglia non esiste, vi è un errore. + if ( nGrid < 0 || nGrid > 2) + return false ; + // Se gli indici sono di frontiera per lo + // Zmap o non esistono, non sono interni. + if ( nI <= - 1 || nI >= int( m_nNx[nGrid]) || + nJ <= - 1 || nJ >= int( m_nNy[nGrid])) + return false ; + // Numero del dexel + int nDex = nJ * m_nNx[nGrid] + nI ; + + // Ciclo sui segmenti del dexel + bool bNotFound = true ; + for ( int nK = 0 ; nK < int( m_Values[nGrid][nDex].size()) ; ++ nK) { + if ( nGrid == 0) { + // Punti estremi del segmento + Point3d ptSt( ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMin) ; + Point3d ptEn( ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMax) ; + // Se il segmento interseca il piano abbiamo finito + Point3d ptInt ; + if ( IntersLinePlane( ptSt, ptEn, plPlane, ptInt)) + return true ; + } + else if ( nGrid == 1) { + // Punti estremi del segmento + Point3d ptSt( m_Values[nGrid][nDex][nK].dMin, ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep) ; + Point3d ptEn( m_Values[nGrid][nDex][nK].dMax, ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep) ; + // Se il segmento interseca il piano abbiamo finito + Point3d ptInt ; + if ( IntersLinePlane( ptSt, ptEn, plPlane, ptInt)) + return true ; + } + else { + // Punti estremi del segmento + Point3d ptSt( ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMin, ( nI + 0.5) * m_dStep) ; + Point3d ptEn( ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMax, ( nI + 0.5) * m_dStep) ; + // Se il segmento interseca il piano abbiamo finito + Point3d ptInt ; + if ( IntersLinePlane( ptSt, ptEn, plPlane, ptInt)) + return true ; + } + } + return false ; +} + +//---------------------------------------------------------------------------- +int +VolZmap::CalcIndexForPlaneCells( const Plane3d& plPlane, int nGrid, int nCellI, int nCellJ) const +{ + int nIndex = 0 ; + + if ( InOut( plPlane, nGrid, nCellI, nCellJ)) + nIndex |= ( 1 << 0) ; + if ( InOut( plPlane, nGrid, nCellI + 1, nCellJ)) + nIndex |= ( 1 << 1) ; + if ( InOut( plPlane, nGrid, nCellI + 1, nCellJ + 1)) + nIndex |= ( 1 << 2) ; + if ( InOut( plPlane, nGrid, nCellI, nCellJ + 1)) + nIndex |= ( 1 << 3) ; + + return nIndex ; +} + +//---------------------------------------------------------------------------- +// dZ è una quota relativa alla mappa nGrid, se il dexel nDex esiste e +// dZ è interna a un suo intervallo restituisce true, false altrimenti. +bool +VolZmap::IsZInsideInterval( int nGrid, int nDex, double dZ) const +{ + // Se la griglia non esiste vi è un errore. + if ( nGrid < 0 || nGrid > 2) + return false ; + // Se il dexel non esiste vi è un errore. + if ( nDex < 0 && nDex > int( m_Values[nGrid].size() - 1)) + return false ; + // Valuto se dZ è interna a un intervallo. + for ( int nk = 0 ; nk < int( m_Values[nGrid][nDex].size()) ; ++ nk) { + double dZ1 = m_Values[nGrid][nDex][nk].dMin ; + double dZ2 = m_Values[nGrid][nDex][nk].dMax ; + // Se troviamo dZ in un intervallo abbiamo finito. + if ( dZ > dZ1 && dZ < dZ2) + return true ; + } + // dZ non sta in un intervallo. + return false ; +} + +//---------------------------------------------------------------------------- +// Box e piano devono essere nello stesso riferimento. Il box deve essere axis aligned. +bool +TestIntersPlaneBox( const BBox3d& b3Box, const Plane3d& plPlane) +{ + // Calcolo le distanze con segno dei punti dal piano + Point3d ptE0 = b3Box.GetMin() ; + double dDist0 = DistPointPlane( ptE0, plPlane) ; + Point3d ptE1( b3Box.GetMax().x, b3Box.GetMin().y, b3Box.GetMin().z) ; + double dDist1 = DistPointPlane( ptE1, plPlane) ; + Point3d ptE2( b3Box.GetMax().x, b3Box.GetMax().y, b3Box.GetMin().z) ; + double dDist2 = DistPointPlane( ptE2, plPlane) ; + Point3d ptE3( b3Box.GetMin().x, b3Box.GetMax().y, b3Box.GetMin().z) ; + double dDist3 = DistPointPlane( ptE3, plPlane) ; + Point3d ptE4( b3Box.GetMin().x, b3Box.GetMin().y, b3Box.GetMax().z) ; + double dDist4 = DistPointPlane( ptE4, plPlane) ; + Point3d ptE5( b3Box.GetMax().x, b3Box.GetMin().y, b3Box.GetMax().z) ; + double dDist5 = DistPointPlane( ptE5, plPlane) ; + Point3d ptE6 = b3Box.GetMax() ; + double dDist6 = DistPointPlane( ptE6, plPlane) ; + Point3d ptE7( b3Box.GetMin().x, b3Box.GetMax().y, b3Box.GetMax().z) ; + double dDist7 = DistPointPlane( ptE7, plPlane) ; + // Distanze tutte positive + if ( dDist0 > EPS_SMALL && dDist1 > EPS_SMALL && dDist2 > EPS_SMALL && dDist3 > EPS_SMALL && + dDist4 > EPS_SMALL && dDist5 > EPS_SMALL && dDist6 > EPS_SMALL && dDist7 > EPS_SMALL) + return false ; + // Distanze tutte negative + if ( dDist0 < - EPS_SMALL && dDist1 < - EPS_SMALL && dDist2 < - EPS_SMALL && dDist3 < - EPS_SMALL && + dDist4 < - EPS_SMALL && dDist5 < - EPS_SMALL && dDist6 < - EPS_SMALL && dDist7 < - EPS_SMALL) + return false ; + // Il piano interseca il box + return true ; +} + +//---------------------------------------------------------------------------- +// Il piano deve essere espresso nel sistema di riferimento intrinseco dello Zmap. +bool +VolZmap::TestIntersPlaneZmapBBox( const Plane3d& plPlane) const +{ + BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ; + return TestIntersPlaneBox( b3Zmap, plPlane) ; +} diff --git a/VolZmap.h b/VolZmap.h index efa9fe8..00e73b8 100644 --- a/VolZmap.h +++ b/VolZmap.h @@ -17,6 +17,7 @@ #include "GeoObjRW.h" #include "CurveComposite.h" #include "SurfFlatRegion.h" +#include "CurveLine.h" #include "/EgtDev/Include/EGkVolZmap.h" #include "/EgtDev/Include/EGkPoint3d.h" #include "/EgtDev/Include/EGkVector3d.h" @@ -99,6 +100,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool GetDepth( const Point3d& ptP, const Vector3d& vtD, double& dInLength, double& dOutLength, bool bExact) const override ; bool GetDepthWithDexel( const Point3d& ptP, const Vector3d& vtDir, double& dInLength, double& dOutLength) const ; bool GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtDir, double& dInLength, double& dOutLength, bool bEnh) const ; + bool GetPlaneIntersection( const Plane3d& plPlane, POCRVVECTOR& vpLoop) const override ; bool AvoidBox( const Frame3d& frBox, const Vector3d& vtDiag) const override ; VolZmap* ClonePart( int nPart) const override ; bool RemovePart( int nPart) override ; @@ -146,15 +148,19 @@ class VolZmap : public IVolZmap, public IGeoObjRW const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DLIST& lstTria) const ; bool ProcessCube( int nVoxI, int nVoxJ, int nVoxK, int& nCubeIndex, TRIA3DLIST& lstTria) const ; bool ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DLIST& lstTria, TriHolder& triHold, bool bExt) const ; ///// ora + bool ProcessCell( int nGrid, int nCellI, int nCellJ, const Plane3d& plPlane, std::vector& vLine) const ; bool ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) const ; bool ExtMarchingCubes( std::vector& vVox, TRIA3DLIST& lstTria, bool bEnh) const ; bool FlipEdgesII( TriHolder& TriHold, bool bGraph) const ; bool FlipEdgesBB( TriaMatrix& InterTria) const ; bool IsThereMat( int nI, int nJ, int nK) const ; - int CalcIndex( int nI, int nJ, int nK) const ; + bool InOut( const Plane3d& plPlane, int nGrid, int nI, int nJ) const ; + int CalcIndex( int nI, int nJ, int nK) const ; + int CalcIndexForPlaneCells( const Plane3d& plPlane, int nGrid, int nCellI, int nCellJ) const ; bool IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, Vector3d& vtNormal) const ; bool IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP, double dPrec = EPS_SMALL) const ; bool GetPointVoxel( const Point3d& ptP, int& nVoxI, int& nVoxJ, int& nVoxK) const ; + bool IsZInsideInterval( int nGrid, int nDex, double dZ) const ; // OPERAZIONI SU INTERVALLI bool SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ, double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax) ; @@ -271,7 +277,8 @@ class VolZmap : public IVolZmap, public IGeoObjRW Point3d& ptInt1, Point3d& ptInt2) ; bool IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir, const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaX, - Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2) ; + Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2) ; + bool TestIntersPlaneZmapBBox( const Plane3d& plPlane) const ; // Voxel: esistenza e passaggio da N a ijk per i voxel bool IsValidVoxel( int nN) const ; bool IsValidVoxel( int nI, int nJ, int nK) const ;