From bf3760808e1a3580fedcc87b69c5ace0748df951 Mon Sep 17 00:00:00 2001 From: Daniele Bariletti Date: Tue, 5 Mar 2024 17:17:12 +0100 Subject: [PATCH] EgtGeomKernel : - correzione di errori nella gestione del taglio sul polo e sul bordo di chiusura di una sup. di Bezier. --- SurfBezier.cpp | 544 ++++++++++++++++++++++++++++++++----------------- SurfBezier.h | 1 + 2 files changed, 358 insertions(+), 187 deletions(-) diff --git a/SurfBezier.cpp b/SurfBezier.cpp index d8eca30..39f0e90 100644 --- a/SurfBezier.cpp +++ b/SurfBezier.cpp @@ -1612,6 +1612,42 @@ SurfBezier::ResetTrimRegion( void) m_OGrMgr.Reset() ; } +//---------------------------------------------------------------------------- +bool +SurfBezier::IncreaseUV( double& dUV, double dxy, bool bUOrV) const +{ + dUV += dxy ; + if ( bUOrV) { + if ( dUV < 0) { + if ( m_bClosedU) + dUV += m_nSpanU * SBZ_TREG_COEFF ; + else + dUV = 0 ; + } + if ( dUV > m_nSpanU * SBZ_TREG_COEFF ) { + if ( m_bClosedU) + dUV -= m_nSpanU * SBZ_TREG_COEFF ; + else + dUV = m_nSpanU * SBZ_TREG_COEFF ; + } + } + else { + if ( dUV ) { + if ( m_bClosedV) + dUV += m_nSpanV * SBZ_TREG_COEFF ; + else + dUV = 0 ; + } + if ( dUV > m_nSpanV * SBZ_TREG_COEFF ) { + if ( m_bClosedU) + dUV -= m_nSpanV * SBZ_TREG_COEFF ; + else + dUV = m_nSpanV * SBZ_TREG_COEFF ; + } + } + return true ; +} + //---------------------------------------------------------------------------- bool SurfBezier::UnprojectCurveFromStm( const ICurveComposite* pCC, ICRVCOMPOPVECTOR& vpCC) const @@ -1628,43 +1664,79 @@ SurfBezier::UnprojectCurveFromStm( const ICurveComposite* pCC, ICRVCOMPOPVECTOR& Point3d pt3D, pt2D ; pCrv0->GetStartPoint( pt3D) ; Point3d pt3DEnd ; pCrv0->GetEndPoint( pt3DEnd) ; bool bThroughEdge = false ; + BOOLVECTOR vbThroughEdge ; if ( ! UnprojectPoint(pt3D, pt2D, pt3DEnd, &bThroughEdge)) return false ; + vbThroughEdge.push_back( bThroughEdge) ; //pl.AddUPoint( 0, pt2D) ; // aggiungo tutti i successivi BIPNTVECTOR vBPnt ; - BOOLVECTOR vbThroughEdge ; + bThroughEdge = false ; + int nRejected = 0 ; for ( int i = 0 ; i < int( pCC->GetCurveCount()) ; ++i) { const ICurve* pCrv = pCC->GetCurve( i) ; Point3d pt3DPrev = pt3D ; Point3d pt2DPrev = pt2D ; - if ( bThroughEdge ) { - // devo cambiare le coordinate di pt3DPrev per periodicità + bool bPrevIsPole = false ; + if ( bThroughEdge) { + // devo cambiare le coordinate di pt2DPrev per periodicità // capisco su quale lato è e lo porto sul lato opposto - if ( m_bClosedU ) { + Point3d pt = pt2DPrev ; + if ( m_bClosedU) { if ( pt2DPrev.x < 1) pt2DPrev.x = m_nSpanU * SBZ_TREG_COEFF ; else if ( (m_nSpanU * SBZ_TREG_COEFF - pt2DPrev.x) < 1) pt2DPrev.x = 0 ; } - if ( m_bClosedV ) { + if ( m_bClosedV) { if ( pt2DPrev.y < 1) pt2DPrev.y = m_nSpanV * SBZ_TREG_COEFF ; else if ( (m_nSpanV * SBZ_TREG_COEFF - pt2DPrev.y) < 1) pt2DPrev.y = 0 ; } + bPrevIsPole = AreSamePointApprox( pt, pt2DPrev) ; } pCrv->GetEndPoint( pt3D) ; if ( ! UnprojectPoint( pt3D, pt2D, pt3DPrev, &bThroughEdge)) return false ; - vbThroughEdge.push_back( bThroughEdge) ; + if ( bPrevIsPole) { + // se il punto precedente era di polo allora devo correggere le sue coordinate 2D + if ( ( m_vbPole[0] || m_vbPole[2] )) { + if ( pt2DPrev.y < 1 || m_nSpanV * SBZ_TREG_COEFF - pt2DPrev.y < 1) + pt2DPrev.x = pt2D.x ; + } + if ( ( m_vbPole[1] || m_vbPole[3] )) { + if ( pt2DPrev.x < 1 || m_nSpanU * SBZ_TREG_COEFF - pt2DPrev.x < 1) + pt2DPrev.y = pt2D.y ; + } + } //pl.AddUPoint( i, pt2D) ; Vector3d vtDir = pt2D - pt2DPrev ; vtDir.Normalize() ; + // se mi accorgo che sto per tracciare un taglio lungo un bordo posso semplicmente evitarlo + if ( (1 - abs(vtDir.x) < EPS_SMALL && (pt2D.y < EPS_SMALL || m_nSpanV * SBZ_TREG_COEFF - pt2D.y < 1)) || // parallelo agli edge 0 e 2 e su uno di questi + (1 - abs(vtDir.y) < EPS_SMALL && (pt2D.x < EPS_SMALL || m_nSpanU * SBZ_TREG_COEFF - pt2D.x < 1))) { // parallello agli edge 1 e 3 e su uno di questi + ++ nRejected ; + continue ; + } + if ( bThroughEdge && vbThroughEdge.back()) { + double dParamH, dParamL ; + dParamH = m_nSpanV * SBZ_TREG_COEFF ; + dParamL = m_nSpanU * SBZ_TREG_COEFF ; + // sia questo punto che il precedente sono su un edge, ma il segmento che li unisce non è parallelo ad un edge + // potrei star tracciando un taglio sul bordo di chiusura + // controllo se sto tracciando una linea che unisce due lati di chiusura, allora in realtà sdtarei tracciando un taglio sull'edge e quindi posso non tracciarlo + if ( (abs( vtDir.x) > abs( vtDir.y) && Dist( pt2D, pt2DPrev) > dParamL * 0.5) || + (abs( vtDir.y) > abs( vtDir.x) && Dist( pt2D, pt2DPrev) > dParamH * 0.5)) { + ++ nRejected ; + continue ; + } + } + vbThroughEdge.push_back( bThroughEdge) ; vBPnt.emplace_back( BIPOINT( pt2DPrev, pt2D)) ; - if ( ! chainC.AddCurve( i + 1, pt2DPrev, vtDir, pt2D, vtDir)) + if ( ! chainC.AddCurve( i + 1 - nRejected, pt2DPrev, vtDir, pt2D, vtDir)) return false ; } // ricostruisco le catene in 2D @@ -1686,10 +1758,10 @@ SurfBezier::UnprojectCurveFromStm( const ICurveComposite* pCC, ICRVCOMPOPVECTOR& bAdded = ( Dist( ptStart, ptEnd) > dToler / 2 && pLine->Set( ptStart, ptEnd)) ; bAdded = bAdded && pCC2D->AddCurve( pLine, true, dToler) ; - if ( vbThroughEdge[i]) { - vpCC.emplace_back( Release(pCC2D)) ; - pCC2D.Set( CreateCurveComposite()) ; - } + //if ( vbThroughEdge[i]) { + // vpCC.emplace_back( Release(pCC2D)) ; + // pCC2D.Set( CreateCurveComposite()) ; + //} ptNear = ( bAdded ? ptEnd : ptStart) ; } if ( pCC2D->IsValid()) @@ -2042,11 +2114,13 @@ SurfBezier::Cut( const Plane3d& plPlane, bool bSaveOnEq) if ( icc.GetIntersCount() != 0) { for ( int i = 0 ; i < int( icc.GetIntersCount()); ++i) { IntCrvCrvInfo iccInfo ; icc.GetIntCrvCrvInfo( i, iccInfo) ; - if ( Dist( iccInfo.IciA->ptI, ptStart) < dDistStart) { + if ( ! bStartFound && Dist( iccInfo.IciA->ptI, ptStart) < dDistStart) { + dDistStart = Dist( iccInfo.IciA->ptI, ptStart) ; vICCStart[0] = iccInfo ; tStart = TRINT( c, l, t) ; } - if ( Dist( iccInfo.IciA->ptI, ptEnd) < dDistEnd) { + if ( ! bEndFound && Dist( iccInfo.IciA->ptI, ptEnd) < dDistEnd) { + dDistEnd = Dist( iccInfo.IciA->ptI, ptEnd) ; vICCEnd[0] = iccInfo ; tEnd = TRINT( c, l, t) ; } @@ -2054,23 +2128,47 @@ SurfBezier::Cut( const Plane3d& plPlane, bool bSaveOnEq) } } } + // ricostruisco gli elementi per la mappa mInters - // se ho intersezione con due loop diverse due entry diverse (immediato) - // se ho intersezione con un loop solo allora accorpo i due vettori delle intersezioni - if ( tStart == tEnd) { - vICCStart.emplace_back( vICCEnd[0]) ; - mInters.insert( pair(tStart, vICCStart)) ; + if ( ! bStartFound && ! bEndFound) { + if ( tStart == tEnd) { + // se ho intersezione con un loop solo allora accorpo i due vettori delle intersezioni + vICCStart.emplace_back( vICCEnd[0]) ; + mInters.insert( pair(tStart, vICCStart)) ; + } + else{ + // se ho intersezione con due loop diverse due entry diverse le inserisco + mInters.insert( pair(tStart, vICCStart)) ; + mInters.insert( pair(tEnd, vICCEnd)) ; + } + // ridefinisco il taglio aperto con la sua versione estesa che arriva a toccare i loop dello spazio parametrico + vpCCOpen[t].Set( GetCurveComposite(pCrv->CopyParamRange( vICCStart[0].IciB->dU, vICCEnd[0].IciB->dU))) ; } else { - mInters.insert( pair(tStart, vICCStart)) ; - mInters.insert( pair(tEnd, vICCEnd)) ; + // devo verificare se avevo già trovato una delle due intersezioni e se era sullo stesso loop o no + if ( ! bStartFound) { + if ( mInters.count( tStart) ==1) + mInters[tStart].emplace_back( vICCStart[0]) ; + else + mInters.insert( pair( tStart, vICCStart)) ; + pCrv->TrimStartAtParam( vICCStart[0].IciB->dU) ; + vpCCOpen[t].Set( GetCurveComposite( Release( pCrv))) ; + } + if ( ! bEndFound) { + if ( mInters.count( tEnd) == 1) + mInters[tEnd].emplace_back( vICCEnd[0]) ; + else + mInters.insert( pair( tEnd, vICCEnd)) ; + pCrv->TrimEndAtParam( vICCEnd[0].IciB->dU) ; + vpCCOpen[t].Set( GetCurveComposite( Release( pCrv))) ; + } } } } // vettore di flag che mi indica quali tagli aperti sono stati aggiunti al nuovo bordo BOOLVECTOR vbAdded( vpCCOpen.size()) ; - fill(vbAdded.begin(), vbAdded.end(), false) ; + std::fill( vbAdded.begin(), vbAdded.end(), false) ; PtrOwner pCCNewEdge( CreateCurveComposite()) ; PtrOwner pCL( CreateCurveLine()) ; //bool bAddedAll = false ; @@ -2272,7 +2370,7 @@ SurfBezier::Cut( const Plane3d& plPlane, bool bSaveOnEq) //vtDirS.Rotate( Z_AX, bPos? 90 : -90) ; vtDir = vtDirS + vtDirE ; PtrOwner pCL( CreateCurveLine()) ; - pCL->SetPVL( ptStart, vtDirS, 1e6) ; + pCL->SetPVL( ptStart, vtDir, 1e6) ; IntersCurveCurve icc( *pCL, *pCrv) ; IntCrvCrvInfo iccInfo ; // verifico di guardare verso l'interno ( il numero di intersezioni deve essere pari visto che partivo da un punto sulla curva) @@ -2444,78 +2542,129 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in // se la dimensione maggiore è grande come la dimensione dello spazio parametrico allora potrebbe essere che le coordinate parametriche di un vertice // siano da correggere per periodicità if ( m_bClosedU && abs(vtDir.x) > abs( vtDir.y) && dDist > dParamL * 0.5 ) { - // trovo se dei punti del triangolo sono sul bordo dello spazio parametrico - BOOLVECTOR vbOn(3) ; + // trovo se dei vertici del triangolo sono sul bordo dello spazio parametrico + INTVECTOR vOn(3) ; + fill( vOn.begin(), vOn.end(), -1) ; + int nVertOnPole = -1 ; for ( int p = 0 ; p < 3; ++p ) { - if ( ! m_vCCEdge[1]->IsValid()) { - Point3d pt ; - m_vPLEdge[1].GetFirstPoint( pt) ; - vbOn[p] = AreSamePointApprox( pt, vPT[p]) ; - } - else - vbOn[p] = m_vCCEdge[1]->IsPointOn( vPT[p]) ; - if ( ! m_vCCEdge[3]->IsValid()) { - Point3d pt ; - m_vPLEdge[3].GetFirstPoint( pt) ; - vbOn[p] = vbOn[p] || AreSamePointApprox( pt, vPT[p]) ; - } - else - vbOn[p] = vbOn[p] || m_vCCEdge[3]->IsPointOn( vPT[p]) ; - } - // controllo che almeno un punto sia su un edge - if ( vbOn[0] || vbOn[1] || vbOn[2]) { - double dRightX ; - // tengo per buone le coordinate dei punti che NON sono sul bordo dello spazio parametrico - for ( int p = 0 ; p < 3; ++p) { - if ( !vbOn[p] ) { - dRightX = vPtPa[p].x ; - break ; + for ( int ed = 0 ; ed < 4; ++ed) { + if ( ! m_vCCEdge[ed]->IsValid()) { + Point3d pt ; + m_vPLEdge[ed].GetFirstPoint( pt) ; + if ( AreSamePointApprox( pt, vPT[p])) { + vOn[p] = ed ; + // se un vertice sta su un polo me lo segno + nVertOnPole = p ; + } + } + else { + if (m_vCCEdge[ed]->IsPointOn( vPT[p]) && vOn[p] == -1) + vOn[p] = ed ; } } - for ( int p = 0 ; p < 3; ++p) { - if ( abs(vPtPa[p].x - dRightX) > EPS_SMALL ) { - if ( vPtPa[p].x < EPS_SMALL) - vPtPa[p].x = dParamL ; - else - vPtPa[p].x = 0 ; + } + //// controllo che i vertici che stanno su un edge siano davvero sul lato di chiusura e non sull'eventuale polo + + //// se ho poli + //if ( m_vbPole[0] || m_vbPole[2]) { + // // se un vertice sta su un polo me lo segno + // auto it0 = find(vOn.begin(), vOn.end(), 0) ; + // auto it2 = find(vOn.begin(), vOn.end(), 2) ; + // if ( it0 != vOn.end()) { + // bVertOnPole = true ; + // nVertOnPole = it0 - vOn.begin() ; + // } + // if ( it2 != vOn.end() ) { + // bVertOnPole = true ; + // nVertOnPole = it2 - vOn.begin() ; + // } + //} + + // controllo che almeno un vertice sia su un edge e se è l'unico vertice, che non sia su un polo + if ( vOn[0] > 0 || vOn[1] > 0 || vOn[2] > 0) { + // se ho più un vertice sul lato oppure se ne ho solo uno ma non è sul polo allora procedo alla correzione delle coordinate + if ( vOn[0] * vOn[1] * vOn[2] < 0 || + (vOn[0] > 0 && vOn[0] != nVertOnPole) || + (vOn[1] > 0 && vOn[1] != nVertOnPole) || + (vOn[2] > 0 && vOn[2] != nVertOnPole)) { + double dRightX ; + // tengo per buone le coordinate dei vertici che NON sono sul bordo dello spazio parametrico + for ( int p = 0 ; p < 3; ++p) { + if ( vOn[p] == -1) { + dRightX = vPtPa[p].x ; + break ; + } + } + for ( int p = 0 ; p < 3; ++p) { + if ( abs(vPtPa[p].x - dRightX) > EPS_SMALL ) { + if ( vPtPa[p].x < EPS_SMALL) + vPtPa[p].x = dParamL ; + else + vPtPa[p].x = 0 ; + } } } } } else if ( m_bClosedV && abs(vtDir.y) > abs(vtDir.x) && dDist > dParamH * 0.5) { - BOOLVECTOR vbOn(3) ; + INTVECTOR vOn(3) ; + fill( vOn.begin(), vOn.end(), -1) ; + int nVertOnPole = -1 ; for ( int p = 0 ; p < 3; ++p ) { - if ( ! m_vCCEdge[0]->IsValid()) { - Point3d pt ; - m_vPLEdge[0].GetFirstPoint( pt) ; - vbOn[p] = AreSamePointApprox( pt, vPT[p]) ; - } - else - vbOn[p] = m_vCCEdge[0]->IsPointOn( vPT[p]) ; - if ( ! m_vCCEdge[2]->IsValid()) { - Point3d pt ; - m_vPLEdge[2].GetFirstPoint( pt) ; - vbOn[p] = vbOn[p] || AreSamePointApprox( pt, vPT[p]) ; - } - else - vbOn[p] = vbOn[p] || m_vCCEdge[2]->IsPointOn( vPT[p]) ; - } - // controllo che almeno un punto sia su un edge - if ( vbOn[0] || vbOn[1] || vbOn[2]) { - double dRightY ; - // tengo per buone le coordinate dei punti che NON sono sul bordo dello spazio parametrico - for ( int p = 0 ; p < 3; ++p) { - if ( ! vbOn[p]) { - dRightY = vPtPa[p].y ; - break ; + for (int ed = 0 ; ed < 4 ; ++ed) { + if ( ! m_vCCEdge[ed]->IsValid()) { + Point3d pt ; + m_vPLEdge[ed].GetFirstPoint( pt) ; + if ( AreSamePointApprox( pt, vPT[p])) { + vOn[p] = ed ; + // se un vertice sta su un polo me lo segno + nVertOnPole = p ; + } + } + else { + if( m_vCCEdge[ed]->IsPointOn( vPT[p]) && vOn[p] == -1) + vOn[p] = ed ; } } - for ( int p = 0 ; p < 3; ++p) { - if ( abs(vPtPa[p].y - dRightY) > EPS_SMALL) { - if ( vPtPa[p].y < EPS_SMALL) - vPtPa[p].y = dParamH ; - else - vPtPa[p].y = 0 ; + } + //// controllo che i vertici che stanno su un edge siano davvero sul lato di chiusura e non sull'eventuale polo + //// se ho poli + //if ( m_vbPole[1] || m_vbPole[3]) { + // // se un vertice sta su un polo me lo segno + // auto it1 = find(vOn.begin(), vOn.end(), 1) ; + // auto it3 = find(vOn.begin(), vOn.end(), 3) ; + // if ( it1 != vOn.end()) { + // bVertOnPole = true ; + // nVertOnPole = it1 - vOn.begin() ; + // } + // if ( it3 != vOn.end() ) { + // bVertOnPole = true ; + // nVertOnPole = it3 - vOn.begin() ; + // } + //} + + // controllo che almeno un vertice sia su un edge + if ( vOn[0] > 0 || vOn[1] > 0 || vOn[2] > 0) { + // se ho più un vertice sul lato oppure se ne ho solo uno ma non è sul polo allora procedo alla correzione delle coordinate + if ( vOn[0] * vOn[1] * vOn[2] < 0 || + (vOn[0] > 0 && vOn[0] != nVertOnPole) || + (vOn[1] > 0 && vOn[1] != nVertOnPole) || + (vOn[2] > 0 && vOn[2] != nVertOnPole)) { + double dRightY ; + // tengo per buone le coordinate dei vertici che NON sono sul bordo dello spazio parametrico + for ( int p = 0 ; p < 3; ++p) { + if ( vOn[p] == -1) { + dRightY = vPtPa[p].y ; + break ; + } + } + for ( int p = 0 ; p < 3; ++p) { + if ( abs(vPtPa[p].y - dRightY) > EPS_SMALL) { + if ( vPtPa[p].y < EPS_SMALL) + vPtPa[p].y = dParamH ; + else + vPtPa[p].y = 0 ; + } } } } @@ -2536,7 +2685,7 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in return false ; // trovo quale vertice è sull'edge di polo BOOLVECTOR vbOn(3) ; - fill( vbOn.begin(), vbOn.end(), false) ; + std::fill( vbOn.begin(), vbOn.end(), false) ; for ( int p = 0 ; p < 3; ++p ) { for ( int c = 0 ; c < 4; ++c) { if ( ! m_vCCEdge[c]->IsValid()) { @@ -2599,27 +2748,27 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in ptSP.y = ptParam.y() ; // controllo se le coordinate sono nello spazio parametrico ed eventualmente se posso sistemarle tenendo conto della periodicità ( se la superficie è chiusa) if ( ptSP.x < 0) { - if ( m_bClosedU) - ptSP.x += m_nSpanU * SBZ_TREG_COEFF ; - else + //if ( m_bClosedU) + // ptSP.x += m_nSpanU * SBZ_TREG_COEFF ; + //else ptSP.x = 0 ; } if ( ptSP.x > m_nSpanU * SBZ_TREG_COEFF ) { - if ( m_bClosedU) - ptSP.x -= m_nSpanU * SBZ_TREG_COEFF ; - else + //if ( m_bClosedU) + // ptSP.x -= m_nSpanU * SBZ_TREG_COEFF ; + //else ptSP.x = m_nSpanU * SBZ_TREG_COEFF ; } if ( ptSP.y < 0 ) { - if ( m_bClosedV) - ptSP.y += m_nSpanV * SBZ_TREG_COEFF ; - else + //if ( m_bClosedV) + // ptSP.y += m_nSpanV * SBZ_TREG_COEFF ; + //else ptSP.y = 0 ; } if ( ptSP.y > m_nSpanV * SBZ_TREG_COEFF ) { - if ( m_bClosedU) - ptSP.y -= m_nSpanV * SBZ_TREG_COEFF ; - else + //if ( m_bClosedU) + // ptSP.y -= m_nSpanV * SBZ_TREG_COEFF ; + //else ptSP.y = m_nSpanV * SBZ_TREG_COEFF ; } return true ; @@ -2629,101 +2778,122 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in bool SurfBezier::UnprojectPoint( const Point3d& pt3D, Point3d& ptParam, const Point3d& ptIPrev, bool* bThroughEdge) const { - // dato il punto pt3D sulla superficie di Bezier si cercano le coordinate parametriche ( ptParam) , iterativamente con Newton - // trovato un primo candidato ptParam, ne calcolo l'immagine sulla superficie ( ptBez) e ne calcolo la distanza con il punto pt3D - // ripeto cercando di avvicinarmi il più possibile - // per trovare il primo punto trovo il triangolo della trimesh ausiliaria più vicino e il punto più vicino - DistPointSurfTm dptSurfTm( pt3D, *GetAuxSurf()) ; - Point3d ptI ; dptSurfTm.GetMinDistPoint( ptI) ; - if ( ! UnprojectPointFromStm( -1, ptI, ptParam, 5, ptIPrev, bThroughEdge)) - return false ; - Point3d ptBez ; - GetPointD1D2( ptParam.x / SBZ_TREG_COEFF, ptParam.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptBez) ; - // usando un algoritmo di newton cerco di avvicinarmi il più possibile al punto - double dDistNew = Dist( pt3D, ptBez) ; - double dDistPre ; - double dDist0, dDist1; + // dato il punto pt3D sulla superficie di Bezier si cercano le coordinate parametriche ( ptParam) , iterativamente con Newton + // trovato un primo candidato ptParam, ne calcolo l'immagine sulla superficie ( ptBez) e ne calcolo la distanza con il punto pt3D + // ripeto cercando di avvicinarmi il più possibile + // per trovare il primo punto trovo il triangolo della trimesh ausiliaria più vicino e il punto più vicino + DistPointSurfTm dptSurfTm( pt3D, *GetAuxSurf()) ; + Point3d ptI ; dptSurfTm.GetMinDistPoint( ptI) ; + if ( ! UnprojectPointFromStm( -1, ptI, ptParam, 5, ptIPrev, bThroughEdge)) + return false ; + Point3d ptBez ; + GetPointD1D2( ptParam.x / SBZ_TREG_COEFF, ptParam.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptBez) ; + // usando un algoritmo di newton cerco di avvicinarmi il più possibile al punto + double dDistNew = Dist( pt3D, ptBez) ; + double dDistPre ; + double dDist0, dDist1; - int nCount = 0 ; - double dh = EPS_SMALL ; - // metodo di newton in più dimensioni - // vario sia il parametro U che il parametro V e verifico se la distanza dalla retta diminuisce per scostamenti positivi o negativi. - bool bRetry = false ; - double dApproach = 0.01 ; - bool bDesperate = false ; - double dAng = 1 ; - double dr = 5 ; - double dfdU, dfdV ; - while ( dDistNew > 2 * EPS_SMALL && nCount < 100) { - if ( ! bRetry) { - dDistPre = dDistNew ; - // derivata in U - Point3d ptIBzNew1 ; - GetPointD1D2( ( ptParam.x + dh) / SBZ_TREG_COEFF, ptParam.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew1) ; - dDist0 = Dist( pt3D, ptIBzNew1) ; - dfdU = ( dDist0 - dDistPre) / dh ; - // derivata in V - Point3d ptIBzNew2 ; - GetPointD1D2( ptParam.x / SBZ_TREG_COEFF, ( ptParam.y + dh) / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew2) ; - dDist1 = Dist( pt3D, ptIBzNew2) ; - dfdV = ( dDist1 - dDistPre) / dh ; - } - // calcolo le nuove coordinate - Vector3d vtDir ; - double dASum = abs(dfdU) + abs(dfdV) ; - double dSSum = sqrt( pow(dfdU,2) + pow( dfdV,2)) ; - vtDir.Set( - dfdU, - dfdV, 0) ; - if ( ! vtDir.Normalize() ) - vtDir.Set( - dfdU / dSSum, - dfdV / dSSum, 0) ; - dr = dDistPre / dASum ; - // in modalità Retry riduco lo spostamento - vtDir *= dr * ( bRetry ? 0.1 : 0.5) ; - if ( bDesperate) { - // in depserate mode riduco lo spostamento e comincio a cambiare progressivamente la direzione oscillando tra destra e sinistra della direzione "naturale" + int nCount = 0 ; + double dh = EPS_SMALL ; + // metodo di newton in più dimensioni + // vario sia il parametro U che il parametro V e verifico se la distanza dalla retta diminuisce per scostamenti positivi o negativi. + bool bRetry = false ; + double dApproach = 0.01 ; + bool bDesperate = false ; + double dAng = 1 ; + double dr = 5 ; + double dfdU, dfdV ; + //// se il punto era su un edge allora l'ho già trovato giusto + //if ( bThroughEdge != nullptr && *bThroughEdge) + // return true ; + while ( dDistNew > 2 * EPS_SMALL && nCount < 100) { + if ( ! bRetry) { + dDistPre = dDistNew ; + // derivata in U + Point3d ptIBzNew1 ; + GetPointD1D2( ( ptParam.x + dh) / SBZ_TREG_COEFF, ptParam.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew1) ; + dDist0 = Dist( pt3D, ptIBzNew1) ; + dfdU = ( dDist0 - dDistPre) / dh ; + // derivata in V + Point3d ptIBzNew2 ; + GetPointD1D2( ptParam.x / SBZ_TREG_COEFF, ( ptParam.y + dh) / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew2) ; + dDist1 = Dist( pt3D, ptIBzNew2) ; + dfdV = ( dDist1 - dDistPre) / dh ; + } + // calcolo le nuove coordinate + Vector3d vtDir ; + double dASum = abs(dfdU) + abs(dfdV) ; + double dSSum = sqrt( pow(dfdU,2) + pow( dfdV,2)) ; + vtDir.Set( - dfdU, - dfdV, 0) ; + if ( ! vtDir.Normalize() ) + vtDir.Set( - dfdU / dSSum, - dfdV / dSSum, 0) ; + dr = dDistPre / dASum ; + // in modalità Retry riduco lo spostamento + vtDir *= dr * ( bRetry ? 0.1 : 0.5) ; + if ( bDesperate) { + // in depserate mode riduco lo spostamento e comincio a cambiare progressivamente la direzione oscillando tra destra e sinistra della direzione "naturale" + vtDir *= 0.5 ; + dAng *= -1.25 ; + // riduco ulteriormente lo spostamento dopo qualche giro in desperate mode + if ( abs( dAng) > 5) vtDir *= 0.5 ; - dAng *= -1.25 ; - // riduco ulteriormente lo spostamento dopo qualche giro in desperate mode - if ( abs( dAng) > 5) - vtDir *= 0.5 ; - vtDir.Rotate( Z_AX, dAng) ; - } - ptParam += vtDir ; - // calcolo la nuova distanza tra il punto di partenza e quello che sto trovando con le coordinate attuali - GetPointD1D2( ptParam.x / SBZ_TREG_COEFF, ptParam.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptBez) ; - dDistNew = Dist( pt3D, ptBez) ; - dApproach = dDistPre - dDistNew ; - // se ho peggiorato la situazione rispetto allo step precedente torno indietro e vado in Retry mode - if ( dApproach < EPS_ZERO) { - if ( bRetry) - bDesperate = true ; // entro in desperate mode - ptParam -= vtDir ; - bRetry = true ; - // se ho già fatto 9 giri in desperate mode allora esco - if ( abs( dAng) > 8) - break ; - } - else { - bRetry = false ; - bDesperate = false ; - dAng = 1 ; - } - ++nCount ; + vtDir.Rotate( Z_AX, dAng) ; } - // se il punto era su un edge allora verifico che sia ancora su un edge, sennò ce lo riporto - if ( bThroughEdge != nullptr && *bThroughEdge) { - if ( m_bClosedU ) { - if ( ptParam.x < abs( m_nSpanU * SBZ_TREG_COEFF - ptParam.x)) - ptParam.x = 0 ; - else - ptParam.x = m_nSpanU * SBZ_TREG_COEFF ; - } - if ( m_bClosedV ) { - if ( ptParam.y < abs( m_nSpanV * SBZ_TREG_COEFF - ptParam.y)) - ptParam.y = 0 ; - else - ptParam.y = m_nSpanV * SBZ_TREG_COEFF ; - } + ptParam += vtDir ; + // calcolo la nuova distanza tra il punto di partenza e quello che sto trovando con le coordinate attuali + GetPointD1D2( ptParam.x / SBZ_TREG_COEFF, ptParam.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptBez) ; + dDistNew = Dist( pt3D, ptBez) ; + dApproach = dDistPre - dDistNew ; + // se ho peggiorato la situazione rispetto allo step precedente torno indietro e vado in Retry mode + if ( dApproach < EPS_ZERO) { + if ( bRetry) + bDesperate = true ; // entro in desperate mode + ptParam -= vtDir ; + bRetry = true ; + // se ho già fatto 9 giri in desperate mode allora esco + if ( abs( dAng) > 8) + break ; } + else { + bRetry = false ; + bDesperate = false ; + dAng = 1 ; + } + ++nCount ; + } + //// se il punto era su un edge allora verifico che sia ancora su un edge, sennò ce lo riporto + //if ( bThroughEdge != nullptr && *bThroughEdge) { + // if ( m_bClosedU) { + // if ( abs( m_nSpanU * SBZ_TREG_COEFF - ptParam.x) < 1) + // ptParam.x = 0 ; + // else if ( ptParam.x < 1) + // ptParam.x = m_nSpanU * SBZ_TREG_COEFF ; + // } + // if ( m_bClosedV) { + // if ( abs( m_nSpanV * SBZ_TREG_COEFF - ptParam.y) < 1) + // ptParam.y = 0 ; + // else if ( ptParam.y < 1) + // ptParam.y = m_nSpanV * SBZ_TREG_COEFF ; + // } + //} + + // se il punto era su un edge allora verifico che sia ancora su un edge, sennò ce lo riporto + // guardo a quale dei due lati sono più vicino + // devo distinguere il caso di un triangolo a metà dello spazio, con un vertice su un lato di polo, ma senza vertici su lati di chiusura + if ( bThroughEdge != nullptr && *bThroughEdge) { + if ( m_bClosedU) { + if ( ptParam.x < 1) + ptParam.x = 0 ; + else if ( abs( m_nSpanU * SBZ_TREG_COEFF - ptParam.x) < 1) + ptParam.x = m_nSpanU * SBZ_TREG_COEFF ; + } + if ( m_bClosedV) { + if ( ptParam.y < 1) + ptParam.y = 0 ; + else if ( abs( m_nSpanV * SBZ_TREG_COEFF - ptParam.y) < 1) + ptParam.y = m_nSpanV * SBZ_TREG_COEFF ; + } + } return nCount != 100 || (dDistNew < dDistPre ? dDistNew : dDistPre) < 10 * EPS_SMALL ; } diff --git a/SurfBezier.h b/SurfBezier.h index 195f7fd..f424972 100644 --- a/SurfBezier.h +++ b/SurfBezier.h @@ -119,6 +119,7 @@ class SurfBezier : public ISurfBezier, public IGeoObjRW bool UnprojectCurveFromStm( const ICurveComposite* pCC, ICRVCOMPOPVECTOR& vpCC) const override ; bool Cut( const Plane3d& plPlane, bool bSaveOnEq = false) override ; bool CalcPoles( void) override ; + bool IncreaseUV( double& dU, double dx, bool bUOrV) const override ; public : // IGeoObjRW int GetNgeId( void) const override ;