From d3d7f94c3aa63a65761f0bde885e70b11db0eb16 Mon Sep 17 00:00:00 2001 From: Daniele Bariletti Date: Tue, 13 Jun 2023 11:55:33 +0200 Subject: [PATCH] =?UTF-8?q?EgtGeomKernel=20:=20-=20risolti=20i=20bug=20not?= =?UTF-8?q?i=20-=20gestiti=20problemi=20numerici=20al=20bordo=20delle=20ce?= =?UTF-8?q?lle=20Manca=20-=20test=20in=20casi=20pi=C3=B9=20complicati?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- SurfBezier.cpp | 10 ++-- Tree.cpp | 156 +++++++++++++++++++++++++++++++++++++------------ Tree.h | 21 ++++++- 3 files changed, 144 insertions(+), 43 deletions(-) diff --git a/SurfBezier.cpp b/SurfBezier.cpp index 36bacaa..c072092 100644 --- a/SurfBezier.cpp +++ b/SurfBezier.cpp @@ -1500,8 +1500,9 @@ SurfBezier::GetAuxSurf( void) const // costruttore della superficie Tree Tree( this, true) ; - Tree.BuildTree_test() ; - //Tree.BuildTree( 0.2, 5) ; + //Tree.BuildTree_test() ; + //Tree.BuildTree( 0.05, 1) ; + Tree.BuildTree( 0.5, 10) ; //Tree.BuildTree( 1, 25) ; vector vvPL ; Tree.GetPolygons( vvPL) ; @@ -1546,8 +1547,9 @@ bool SurfBezier::GetLeaves( std::vector>& vLeaves) const { Tree Tree( this, true) ; - Tree.BuildTree_test() ; - //Tree.BuildTree( 0.2, 5) ; + //Tree.BuildTree_test() ; + //Tree.BuildTree( 0.05, 1) ; + Tree.BuildTree( 0.5, 10) ; //Tree.BuildTree( 1, 25) ; std::vector vCells ; Tree.GetLeaves( vCells) ; diff --git a/Tree.cpp b/Tree.cpp index d5f3f04..23c5e5d 100644 --- a/Tree.cpp +++ b/Tree.cpp @@ -69,7 +69,7 @@ Cell::IsLeaf ( void) const //---------------------------------------------------------------------------- Tree::Tree( void) - : m_pSrfBz( nullptr), m_bTrimmed( false), m_bBilinear( false), m_bMulti( false), m_bClosed( false) + : m_pSrfBz( nullptr), m_bTrimmed( false), m_bBilinear( false), m_bMulti( false), m_bClosed( false), m_bSplitPatches( true) { Point3d ptBl( 0, 0), ptTr ( 1 * SBZ_TREG_COEFF, 1 * SBZ_TREG_COEFF) ; Cell cRoot( ptBl, ptTr) ; @@ -78,9 +78,9 @@ Tree::Tree( void) //---------------------------------------------------------------------------- Tree::Tree( const SurfBezier* pSrfBz, bool bSplitPatches) - : m_bBilinear( false), m_bMulti( false), m_bClosed( false) + : m_bBilinear( false), m_bMulti( false), m_bClosed( false), m_bSplitPatches( bSplitPatches) { - SetSurf( pSrfBz, bSplitPatches) ; + SetSurf( pSrfBz, m_bSplitPatches) ; } //---------------------------------------------------------------------------- @@ -93,6 +93,7 @@ void Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches) { m_pSrfBz = pSrfBz ; + m_bSplitPatches = bSplitPatches ; // le coordinate delle celle sono nello spazio parametrico int nDegU, nDegV, nSpanU, nSpanV ; bool bIsRat, bTrimmed ; @@ -100,6 +101,8 @@ Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches) m_bTrimmed = bTrimmed ; m_nDegU = nDegU ; m_nDegV = nDegV ; + m_nSpanU = nSpanU ; + m_nSpanV = nSpanV ; if ( nDegU == 1 && nDegV == 1) m_bBilinear = true ; if ( nSpanU * nSpanV != 1) @@ -176,7 +179,7 @@ Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches) int nId = -1 ; for ( int i = 1 ; i < nSpanU ; ++i) { m_mTree[nId].SetSplitDirVert( true) ; - Split( nId, i) ; + Split( nId, i * SBZ_TREG_COEFF) ; ++ nId ; ++ nId ; } @@ -185,11 +188,13 @@ Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches) for ( int nId : vLeaves) { for ( int j = nSpanV - 1 ; j > 0 ; --j ) { m_mTree[nId].SetSplitDirVert( false) ; - Split( nId, j) ; + Split( nId, j * SBZ_TREG_COEFF) ; nId = m_mTree[nId].m_nChild2 ; } } - // split preliminari per dividere le patch in modo da triangolarle indipendentemente//////////////////////////////////////////////////////// + vLeaves.clear() ; + GetHeightLeaves( -1, vLeaves) ; + m_vnParents = vLeaves ; } // se non ho già splittato le patches, controllo se la superficie è chiusa. In tal caso la splitto sul parametro su cui è chiusa else { @@ -603,6 +608,8 @@ Tree::BuildTree( double dLinTol, double dSideMin, double dSideMax) m_mTree[nCToSplit].SetProcessed() ; // risalgo i parent finché non trovo il primo Child2 da processare nCToSplit = m_mTree[nCToSplit].m_nParent ; + if ( nCToSplit == -2) + return true ; if ( m_mTree[m_mTree[nCToSplit].m_nChild1].IsProcessed() && m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed()) m_mTree[nCToSplit].SetProcessed() ; while ( m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed()) { @@ -736,6 +743,8 @@ Tree::BuildTree( double dLinTol, double dSideMin, double dSideMax) m_mTree[nCToSplit].SetProcessed() ; // risalgo i parent finché non trovo il primo Child2 da processare nCToSplit = m_mTree[nCToSplit].m_nParent ; + if ( nCToSplit == -2) + return true ; if ( m_mTree[m_mTree[nCToSplit].m_nChild1].IsProcessed() && m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed()) m_mTree[nCToSplit].SetProcessed() ; while ( m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed()) { @@ -1215,7 +1224,8 @@ Tree::GetPolygons( std::vector& vPolygons) { } // trimmata else { - TraceLoopLabelCell() ; + if ( ! TraceLoopLabelCell()) + return false ; POLYLINEVECTOR vPolygonsBasic ; GetPolygonsBasic( vPolygonsBasic) ; // scorro sui poligoni delle celle non trimmate @@ -1416,12 +1426,19 @@ Tree::FindCell( const Point3d& ptToAssign, CurveLine& clTrim) const return nCells ; } - // se lo spazio parametrico non è stato suddiviso ( ho un piano), restituisco la cella root - if ( m_mTree.at(nId).IsLeaf() ) { - nCells.push_back( nId) ; - return nCells ; - } + //// se lo spazio parametrico non è stato suddiviso ( ho un piano), restituisco la cella root + //if ( m_mTree.at(nId).IsLeaf() ) { + // nCells.push_back( nId) ; + // return nCells ; + //} + // se ho diviso preliminarmente le patches e in uno dei due parametri ho un numero dispari di patches devo individuare a mano la cella parent + // in cui individuare la foglia giusta + if ( m_bSplitPatches && ( (m_nSpanU % 2 != 0 && m_nSpanU > 1) || (m_nSpanV % 2 != 0 && m_nSpanV > 1))) { + INTVECTOR nParents = FindCell( ptToAssign, clTrim, m_vnParents) ; + nId = nParents.back() ; + } + // individuo la foglia in cui ho lo start del loop while ( ! m_mTree.at(nId).IsLeaf()) { if ( m_mTree.at(nId).IsSplitVert()) { double dMid = ( m_mTree.at(nId).GetBottomLeft().x + m_mTree.at(nId).GetTopRight().x) / 2 ; @@ -1442,11 +1459,14 @@ Tree::FindCell( const Point3d& ptToAssign, CurveLine& clTrim) const } } } - if ( nId != -1) - nCells.push_back( nId) ; + //if ( nId != -1) + // nCells.push_back( nId) ; + //if ( nId == -1) // nId = -2 ; + nCells.push_back( nId) ; + //// devo verificare se sono sul bordo di una cella: in questo caso dovrei restituire anche il vicino con cui condivido quel lato!! ///////////////////// Point3d ptBr( m_mTree.at(nId).GetTopRight().x , m_mTree.at(nId).GetBottomLeft().y) ; @@ -1482,7 +1502,7 @@ Tree::FindCell( const Point3d& ptToAssign, CurveLine& clTrim) const if ( AreSamePointApprox(ptToAssign, ptTl) || AreSamePointApprox( ptToAssign, m_mTree.at(nId).GetBottomLeft()) || AreSamePointApprox(ptToAssign, ptBr) || - AreSamePointApprox( ptToAssign, ptTl)) + AreSamePointApprox( ptToAssign, m_mTree.at(nId).GetTopRight())) { Point3d ptToAssignPlus ; double dParam ; @@ -1602,6 +1622,7 @@ Tree::TraceLoopLabelCell( void) INTVECTOR :: iterator iter = find( m_vnLeaves.begin(), m_vnLeaves.end(), nId) ; int nIdPolygon = std::distance( m_vnLeaves.begin(), iter) ; bool bVertex = false ; + bool bEraseNextPoint = false ; while ( plLoop.GetNextPoint( ptCurr)) { // sto uscendo dalla cella, quindi cerco l'intersezione Point3d ptTStart, ptTEnd ; @@ -1610,10 +1631,28 @@ Tree::TraceLoopLabelCell( void) CurveLine clTrim ; clTrim.Set( ptTStart, ptTEnd) ; while( ! IsPointInsidePolyLine( ptCurr, vplPolygons[nIdPolygon], dLinTol)) { /// qui devo mettere una tolleranza negativa per poter tener conto anche dei punti che sono SULLA curva + if ( bEraseNextPoint ) { + vptInters.pop_back() ; + bEraseNextPoint = false ; + } bLoopInside = false ; - // trovo l'intersezione e passo alla cella successiva. nId viene aggiornato dalla funzione FindInters + // trovo l'intersezione e passo alla cella successiva. nId viene aggiornato dalla funzione FindInters bVertex = false ; - FindInters( nId, clTrim, vptInters, bVertex) ; + // se non trovo l'intersezione vuol dire che non sono nella cella giusta! + // al precedente FindInters avrei dovuto passare di cella + if ( ! FindInters( nId, clTrim, vptInters, bVertex)) { + // scaterò il punto molto vicino al lato e tengo solo l'intersezione del trim col lato + m_mTree[nId].m_vInters.back().vpt.pop_back() ; + plLoop.GetPrevPoint( ptTEnd) ; + plLoop.GetPrevPoint( ptTStart) ; + plLoop.GetNextPoint( ptCurr) ; + clTrim.Set( ptTStart, ptTEnd) ; + //clTrim.ExtendEndByLen( EPS_SMALL * 2) ; + vptInters.clear() ; + if ( FindInters( nId, clTrim, vptInters, bVertex)) + return false ; + bEraseNextPoint = true ; + } // ricalcolo la posizione di nId nel vettore delle foglie iter = find( m_vnLeaves.begin(), m_vnLeaves.end(), nId) ; nIdPolygon = std::distance( m_vnLeaves.begin(), iter) ; @@ -1644,7 +1683,8 @@ Tree::TraceLoopLabelCell( void) else { // verifico se sono effettivamente nella cella di partenza o nella penultima cella // devo verificare se il loop finisce in un vertice! in questo caso sono nella penultima cella! e devo fare un FindInters di nuovo per sistemare - Point3d ptLast = m_mTree[nId].m_vInters.back().vpt.back() ; + //Point3d ptLast = m_mTree[nId].m_vInters.back().vpt.back() ; + Point3d ptLast = ptCurr ; if ( nId != nFirstCell) { Point3d ptTl( m_mTree[nId].GetBottomLeft().x, m_mTree[nId].GetTopRight().y) ; Point3d ptBr( m_mTree[nId].GetTopRight().x, m_mTree[nId].GetBottomLeft().y) ; @@ -1661,21 +1701,43 @@ Tree::TraceLoopLabelCell( void) vptInters.clear() ; nLastId = nId ; FindInters( nId, clTrim, vptInters, bVertex) ; - m_mTree[nLastId].m_vInters.back().vpt.push_back( vptInters[0]) ; + if ( ! AreSamePointExact( vptInters[0], m_mTree[nLastId].m_vInters.back().vpt.back())) + m_mTree[nLastId].m_vInters.back().vpt.push_back( vptInters[0]) ; } } - int nOut = m_mTree[nId].m_vInters[nPass].nOut ; - m_mTree[nId].m_vInters.back().vpt.insert( m_mTree[nId].m_vInters.back().vpt.end(), - m_mTree[nId].m_vInters[nPass].vpt.begin(), - m_mTree[nId].m_vInters[nPass].vpt.end()) ; - m_mTree[nId].m_vInters[nPass] = m_mTree[nId].m_vInters.back() ; - m_mTree[nId].m_vInters.pop_back() ; - // sistemo il lato d'uscita - m_mTree[nId].m_vInters[nPass].nOut = nOut ; + // sono tornato nella cella iniziale, quindi giunto i due vettori intersezione + if ( nId == nFirstCell) { + int nOut = m_mTree[nId].m_vInters[nPass].nOut ; + m_mTree[nId].m_vInters.back().vpt.insert( m_mTree[nId].m_vInters.back().vpt.end(), + m_mTree[nId].m_vInters[nPass].vpt.begin(), + m_mTree[nId].m_vInters[nPass].vpt.end()) ; + m_mTree[nId].m_vInters[nPass] = m_mTree[nId].m_vInters.back() ; + m_mTree[nId].m_vInters.pop_back() ; + // sistemo il lato d'uscita + m_mTree[nId].m_vInters[nPass].nOut = nOut ; + } + // se non sono tornato alla cella di partenza + else { + m_mTree[nId].m_vInters.pop_back() ; + if ( (int)m_mTree[nId].m_vInters.size() == 0) + m_mTree[nId].m_nFlag = 0 ; + // aggiusto nIn della cella di partenza ( in questo caso dovrei aver avuto un ingresso direttamente in un vertice) + Point3d ptFirst = m_mTree[nFirstCell].m_vInters[nPass].vpt[0] ; + Point3d ptTl( m_mTree[nFirstCell].GetBottomLeft().x, m_mTree[nFirstCell].GetTopRight().y) ; + Point3d ptBr( m_mTree[nFirstCell].GetTopRight().x, m_mTree[nFirstCell].GetBottomLeft().y) ; + if ( AreSamePointApprox( m_mTree[nFirstCell].GetTopRight(), ptFirst)) + m_mTree[nFirstCell].m_vInters[nPass].nIn = 7 ; + else if (AreSamePointApprox( ptTl, ptFirst)) + m_mTree[nFirstCell].m_vInters[nPass].nIn = 4 ; + else if (AreSamePointApprox( m_mTree[nFirstCell].GetBottomLeft(), ptFirst)) + m_mTree[nFirstCell].m_vInters[nPass].nIn = 5 ; + else if (AreSamePointApprox( ptBr, ptFirst)) + m_mTree[nFirstCell].m_vInters[nPass].nIn = 6 ; + } } } - // riordino i vettori di intersezione per ogni cella e setto il flag RightEdgeIn + // riordino i vettori di intersezione per ogni cella e setto il flag RightEdgeIn for ( int nId : m_vnLeaves) { std::sort( m_mTree[nId].m_vInters.begin( ), m_mTree[nId].m_vInters.end()) ; SetRightEdgeIn( nId) ; @@ -1913,16 +1975,21 @@ Tree::FindInters( int& nId, CurveLine& clTrim, PNTVECTOR& vptInters, bool& bVert // intersezione e controlli IntersLineLine illExit( clTrim, clEdge, true) ; IntCrvCrvInfo aInfo, aInfo2 ; - //if ( ! illExit.GetIntCrvCrvInfo( aInfo) ) - // return false ; - illExit.GetIntCrvCrvInfo( aInfo) ; + bool bIntersOn1Found = true ; + if ( ! illExit.GetIntCrvCrvInfo( aInfo) && ! clEdge2.IsValid()) { + bIntersOn1Found = false ; + return false ; + } if ( clEdge2.IsValid()){ IntersLineLine illExit2( clTrim, clEdge2, true) ; //if ( ! illExit2.GetIntCrvCrvInfo( aInfo2) ) // return false ; // verifico su quale dei due lati ho l'intersezione if ( ! illExit2.GetIntCrvCrvInfo( aInfo2)) - ptInters = aInfo.IciA[0].ptI ; + if ( bIntersOn1Found) + ptInters = aInfo.IciA[0].ptI ; + else + return false ; else { // solo intersezione sul lato 2 if ( illExit.GetNumInters() == 0 ) { @@ -1958,7 +2025,10 @@ Tree::FindInters( int& nId, CurveLine& clTrim, PNTVECTOR& vptInters, bool& bVert nEdge = 7 ; } m_mTree[nId].m_vInters.back().nOut = nEdge ; - vptInters.push_back( ptInters) ; + if ( (int)vptInters.size() == 0) + vptInters.push_back( ptInters) ; + else if ( ! AreSamePointExact(ptInters , vptInters.back())) + vptInters.push_back( ptInters) ; // salvo il vettore intersezione per la cella e capisco in quale altra cella passare if ( (int)vptInters.size() == 1) m_mTree[nId].m_vInters.back().vpt.push_back( vptInters[0]) ; @@ -2303,6 +2373,9 @@ Tree::CreateCellPolygons( int nLeafId, std::vector& vPolygons, I int nFirstLoopInPoly = -1 ; INTVECTOR vAddedLoops ; for ( int w = 0 ; w < (int)vToCheckNow.size() ; ++ w) { + if ( m_mTree[nId].m_vInters[w].vpt.size() < 2) { + continue ; + } // indice del loop in m_vInters int j = vToCheckNow[w] ; Inters inA = m_mTree[nId].m_vInters[j] ; @@ -2366,6 +2439,8 @@ Tree::CreateCellPolygons( int nLeafId, std::vector& vPolygons, I // quindi salto al prossimo loop if ( plTrimmedPoly.GetPointNbr() == 1 ) { plTrimmedPoly.Clear() ; + if ( j == nFirstLoopInPoly) + nEdgeIn = -1 ; continue ; } if ( nEdge > 3 && nEdge != 7 ) { @@ -2530,7 +2605,7 @@ Tree::CreateIslandAndHoles( int nLeafId, std::vector& vPolygons, // numero dei loop interni passati Inters inA ; // se ho almeno un loop CW che non è contenuto in un altro poligono o in un loop interno CCW devo aggiungere il bordo - bool bEdgeToAdd = true ; + //bool bEdgeToAdd = true ; bool bAllContained = true ; bool bContained = false ; int nInters = (int) m_mTree[nId].m_vInters.size() ; @@ -2559,7 +2634,8 @@ Tree::CreateIslandAndHoles( int nLeafId, std::vector& vPolygons, break ; } - if ( m_mTree[nId].m_nFlag == 2 && ( bEdgeToAdd || ! bAllContained)) { + //if ( m_mTree[nId].m_nFlag == 2 && ( bEdgeToAdd || ! bAllContained)) { + if ( m_mTree[nId].m_nFlag == 2 && ! bAllContained) { // i loop esterni sono CW, quindi prima dei loop di trim aggiungo il bordo cella Point3d ptVert = m_mTree[nId].GetTopRight() ; plInLoop.AddUPoint( 0, ptVert) ; @@ -2837,7 +2913,11 @@ Tree::AddVertex( int nId, std::vector& vEdgeVertex, PolyLine& plTrimm Point3d ptLast ; plTrimmedPoly.GetLastPoint( ptLast) ; // verifico di essere allineato con un lato, sennò aggiungo e basta - Vector3d vDir = ptToAdd - ptLast ; + Vector3d vDir; + if ( ! AreSamePointExact(ptToAdd, ptLast)) + vDir = ptToAdd - ptLast ; + else + return true ; // se non riesco a normalizzare perché sono troppo vicino ad un vertice allora aggiungo direttamente il vertice if ( ! vDir.Normalize()) { plTrimmedPoly.EraseLastUPoint() ; @@ -3027,7 +3107,7 @@ Tree::SetRightEdgeIn( int nId) m_mTree[nId].m_nRightEdgeIn = 0 ; return true ; } - // se ho inters sul lato destro allora in parte è dentro + // se ho inters sul lato destro ( nEdge == 3) allora in parte è dentro for ( int k = 0 ; k < nPass ; ++ k) { if ( m_mTree[nId].m_vInters[k].nIn == 3 || m_mTree[nId].m_vInters[k].nOut == 3) { m_mTree[nId].m_nRightEdgeIn = 2 ; @@ -3307,7 +3387,7 @@ Tree::CheckIfBetween( Inters& inA, Inters& inB) const } else if ( AreSameEdge(inB.nIn, inA.nIn)) { //devo controllare il loop b sia prima dell'inizio di A - if ( CheckIfBefore(inA.nIn, inB.vpt[0], inA.nIn, inA.vpt[0])) + if ( CheckIfBefore(inA.nIn, inB.vpt[0], inA.vpt[0])) bFound = true ; } else diff --git a/Tree.h b/Tree.h index c1e71ba..39e5333 100644 --- a/Tree.h +++ b/Tree.h @@ -18,6 +18,7 @@ #include "SurfBezier.h" #include "GeoConst.h" #include "CurveLine.h" +#include "/EgtDev/Include/EGkPolyLine.h" struct Inters { int nIn ; @@ -30,7 +31,21 @@ struct Inters { // se ho più intersezioni che entrano in un lato le riordino considerando che percorro i lati in senso antiorario a partire da ptTR bool operator < ( Inters& b) { bool bEqIn = ( nIn == b.nIn) ; - return nIn < b.nIn || + double dAreaA = 0 , dAreaB = 0 ; + if ( bEqIn && nIn == -1) { + PolyLine pl ; + for ( int k = 0 ; k < (int)vpt.size(); ++ k) + pl.AddUPoint( k, vpt[k]) ; + pl.Close() ; + pl.GetAreaXY( dAreaA) ; + pl.Clear() ; + for ( int k = 0 ; k < (int)b.vpt.size(); ++ k) + pl.AddUPoint( k, b.vpt[k]) ; + pl.Close() ; + pl.GetAreaXY( dAreaB) ; + } + return nIn < b.nIn || + ( bEqIn && nIn == -1 && dAreaA > dAreaB) || ( bEqIn && nIn == 0 && vpt[0].x > b.vpt[0].x) || ( bEqIn && nIn == 1 && vpt[0].y > b.vpt[0].y) || ( bEqIn && nIn == 2 && vpt[0].x < b.vpt[0].x) || @@ -153,11 +168,15 @@ private : bool m_bBilinear ; // superficie bilineare bool m_bMulti ; // superficie multi-patch bool m_bClosed ; // superficie chiusa + bool m_bSplitPatches ; // flag che indica se le patches sono state divise prima della creazione dell'albero int m_nDegU ; // grado della superficie nel parametro U int m_nDegV ; // grado della superficie nel parametro V + int m_nSpanU ; + int m_nSpanV ; std::vector m_vPolygons ; // vettore dei poligoni del tree //std::map m_mPolygons ; // mappa dei poligoni std::map m_mTree ; // mappa che contiene tutti i nodi e le foglie dell'albero. -2 è puntatore Null e -1 è root std::map m_mVert ; // mappa che contiene tutti i vertici 3d delle celle del tree. L'Id è lo stesso che la cella ha in m_mTree INTVECTOR m_vnLeaves ; // vettore delle foglie + INTVECTOR m_vnParents ; // vettore delle celle ottenute dalla divisione preliminare in singole patch } ; \ No newline at end of file