diff --git a/CurveAux.cpp b/CurveAux.cpp index ec1416d..5bbb6eb 100644 --- a/CurveAux.cpp +++ b/CurveAux.cpp @@ -555,17 +555,28 @@ CompositeToBezierCurve( const ICurveComposite* pCC, int nDeg, bool bMakeRatOrNot //---------------------------------------------------------------------------- ICurve* -EditBezierCurve( const ICurveBezier* pCrvBezier, int nDeg, bool bMakeRatOrNot) +EditBezierCurve( const ICurveBezier* pCrvBezier, int nDeg, bool bMakeRatOrNot, double dTol) { - // resta da calcolare un errore sull'approssimazione oppure usare la tecnica di spezzare la curva originale in sottocurve e approssimarle con bezier cubiche - // per ridurre molto l'errore - + if( nDeg < 3) + return nullptr ; // dovrei restituire una bezier di grado 2, razionale per poter essere uniforme con le altre curve trasmorate in bezier PtrOwner pCrvNew( pCrvBezier->Clone()) ; - double dErr = 0 ; - if( ! bMakeRatOrNot) - if( ! pCrvNew->MakeNonRational( 10 * EPS_SMALL, dErr)) - return nullptr ; + if( ! bMakeRatOrNot) { + if( ! pCrvNew->MakeNonRational( 10 * EPS_SMALL)) { + // se ho fallito la conversione diretta in curva non razionale allora la spezzo in bezier cubiche + PtrOwner pBezCubics( CreateCurveComposite()) ; + pBezCubics->AddCurve( ApproxBezierWithCubics(pCrvBezier, dTol)) ; + if( IsNull( pBezCubics)) + return nullptr ; + // adatto ogni sottocurva cubica + PtrOwner pCCEdited( CreateCurveComposite()) ; + for ( int i = 0 ; i < pBezCubics->GetCurveCount() ; ++i) { + if( ! pCCEdited->AddCurve( EditBezierCurve( GetCurveBezier( pBezCubics->GetCurve( i)), nDeg, bMakeRatOrNot, dTol)) ) + return nullptr ; + } + return Release( pCCEdited) ; + } + } int nDegCurr = pCrvNew->GetDegree() ; bool bRat = pCrvNew->IsRational() ; // se la curva è già nella forma giusta la restituisco @@ -581,9 +592,22 @@ EditBezierCurve( const ICurveBezier* pCrvBezier, int nDeg, bool bMakeRatOrNot) } else if ( nDegCurr > nDegWanted) { while ( nDegCurr > nDegWanted) { - pCrvNew.Set( BezierDecreaseDegree( pCrvNew, 1)) ; - if ( IsNull( pCrvNew) || ! pCrvNew->IsValid()) - return nullptr ; + ICurveBezier* pCrvDec = BezierDecreaseDegree( pCrvNew, dTol) ; + if( pCrvDec == nullptr || ! pCrvDec->IsValid()) { + // se ho fallito la riduzione di grado entro la tolleranza richiesta allora la spezzo in bezier cubiche prima di adattare + PtrOwner pBezCubics( CreateCurveComposite()) ; + pBezCubics->AddCurve( ApproxBezierWithCubics(pCrvBezier, dTol)) ; + if ( IsNull( pBezCubics) || ! pBezCubics->IsValid()) + return nullptr ; + // adatto ogni sottocurva cubica + PtrOwner pCCEdited( CreateCurveComposite()) ; + for ( int i = 0 ; i < pBezCubics->GetCurveCount() ; ++i) { + if( ! pCCEdited->AddCurve( EditBezierCurve( GetCurveBezier( pBezCubics->GetCurve( i)), nDeg, bMakeRatOrNot, dTol)) ) + return nullptr ; + } + return Release( pCCEdited) ; + } + pCrvNew.Set( pCrvDec) ; -- nDegCurr ; } } @@ -784,11 +808,12 @@ BezierDecreaseDegree(const ICurveBezier* pCrvBezier, double dTol) Point3d ptOld = pCrvBezier->GetControlPoint( r + 1) ; dErr = Dist( ptOld, 0.5 * ( ptCtrlPrev + ptCtrlCurr)) ; } - //// se l'approssimazione dà un errore troppo alto allora annullo tutto // da controllare - //if ( dErr > dTol) - // return nullptr ; - + // se l'approssimazione dà un errore troppo alto allora annullo tutto + // ricalcolo l'errore a mano, per avere un valore più attendibile CalcBezierApproxError( pCrvBezier, pNewBezier, dErr) ; + if ( dErr > dTol) + return nullptr ; + return Release( pNewBezier) ; } diff --git a/CurveBezier.cpp b/CurveBezier.cpp index dda98fa..7ce4149 100644 --- a/CurveBezier.cpp +++ b/CurveBezier.cpp @@ -2309,52 +2309,68 @@ CurveBezier::MakeRationalStandardForm( void) //---------------------------------------------------------------------------- bool -CurveBezier::MakeNonRational( double dTol, double& dErr) +CurveBezier::MakeNonRational( double dTol) { if( ! m_bRat) return true ; - // provo ad approssimare la curva di bezier con una controparte non razionale + // controllo se i pesi sono tutti == 1 allora è una finta razionale e mi basta fare una copia dei punti di controllo + bool bIsActualRat = true ; + for ( int i = 0 ; i < m_nDeg && bIsActualRat ; ++i) + bIsActualRat = bIsActualRat && abs(m_vWeCtrl[i] - 1) > EPS_SMALL ; + bool bOk = true ; - PtrOwner pNewBez( CreateBasicCurveBezier()) ; - int nDeg = m_nDeg + 2 ; - pNewBez->Init( nDeg, false) ; - PNTVECTOR vPntCtrl ; - PNTVECTOR vPntSampling ; - for ( int p = 0 ; p < nDeg + 1; ++p) { - Point3d pt ; GetPointD1D2( double(p) / nDeg, pt) ; - pNewBez->SetControlPoint( p, pt) ; - vPntCtrl.push_back( pt) ; - } - vPntSampling = vPntCtrl ; - int c = 0 ; - dErr = INFINITO ; - while ( dErr > dTol / 100 && c < 100) { - double dErrMax = 0 ; - // calcolo le differenze tra i punti di sampling sulla nuova curva e quelli sulla curva originale - for ( int p = 0 ; p < nDeg + 1; ++p) { - Point3d pt ; pNewBez->GetPointD1D2( double(p) / nDeg, pt) ; - Vector3d vDiff = vPntSampling[p] - pt ; - double dErrLoc = vDiff.Len() ; - if( dErrLoc > dErrMax) - dErrMax = dErrLoc ; - // aggiorno il vettore dei punti di controllo della nuova curva - vPntCtrl[p] += vDiff ; + if ( ! bIsActualRat ) { + PtrOwner pNewBez( CreateBasicCurveBezier()) ; + for ( int p = 0 ; p < m_nDeg ; ++p) { + Point3d pt = GetControlPoint( p) ; + pNewBez->SetControlPoint( p, pt) ; } - dErr = dErrMax ; - // aggiorno i punti di controllo della nuova curva - for ( int i = 0 ; i < nDeg + 1 ; ++i) - pNewBez->SetControlPoint( i, vPntCtrl[i]) ; - ++c ; } - bOk = dErr < 5 * dTol ; - if( bOk) { + else { + // provo ad approssimare la curva di bezier con una controparte non razionale + PtrOwner pNewBez( CreateBasicCurveBezier()) ; + int nDeg = m_nDeg + 2 ; + pNewBez->Init( nDeg, false) ; + PNTVECTOR vPntCtrl ; + PNTVECTOR vPntSampling ; + for ( int p = 0 ; p < nDeg + 1; ++p) { + Point3d pt ; GetPointD1D2( double(p) / nDeg, pt) ; + pNewBez->SetControlPoint( p, pt) ; + vPntCtrl.push_back( pt) ; + } + vPntSampling = vPntCtrl ; + int c = 0 ; + double dErr = INFINITO ; + while ( dErr > dTol && c < 100) { + double dErrMax = 0 ; + // calcolo le differenze tra i punti di sampling sulla nuova curva e quelli sulla curva originale + for ( int p = 0 ; p < nDeg + 1; ++p) { + Point3d pt ; pNewBez->GetPointD1D2( double(p) / nDeg, pt) ; + Vector3d vDiff = vPntSampling[p] - pt ; + double dErrLoc = vDiff.Len() ; + if( dErrLoc > dErrMax) + dErrMax = dErrLoc ; + // aggiorno il vettore dei punti di controllo della nuova curva + vPntCtrl[p] += vDiff ; + } + dErr = dErrMax ; + // aggiorno i punti di controllo della nuova curva + for ( int i = 0 ; i < nDeg + 1 ; ++i) + pNewBez->SetControlPoint( i, vPntCtrl[i]) ; + ++c ; + } + + // calcolo l'errore di approssimazione sulla curva CalcBezierApproxError( this, pNewBez, dErr) ; - // aggiorno la curva di bezier originale con quella approssimata - Init( nDeg, false) ; - for( int i = 0 ; i < nDeg + 1 ; ++i) { - SetControlPoint( i, pNewBez->GetControlPoint( i)) ; - SetControlWeight( i, pNewBez->GetControlWeight( i)) ; + bOk = dErr < dTol ; + if( bOk) { + // aggiorno la curva di bezier originale con quella approssimata + Init( nDeg, false) ; + for( int i = 0 ; i < nDeg + 1 ; ++i) { + SetControlPoint( i, pNewBez->GetControlPoint( i)) ; + SetControlWeight( i, pNewBez->GetControlWeight( i)) ; + } } } diff --git a/CurveBezier.h b/CurveBezier.h index 91d8598..2bee20c 100644 --- a/CurveBezier.h +++ b/CurveBezier.h @@ -151,7 +151,7 @@ class CurveBezier : public ICurveBezier, public IGeoObjRW int GetSingularParam( double& dPar) const override ; bool MakeRational( void) override ; bool MakeRationalStandardForm( void) override ; - bool MakeNonRational( double dTol, double& dErr) override ; + bool MakeNonRational( double dTol) override ; public : // IGeoObjRW int GetNgeId( void) const override ; diff --git a/SurfBezier.cpp b/SurfBezier.cpp index 164c94b..ed01926 100644 --- a/SurfBezier.cpp +++ b/SurfBezier.cpp @@ -4754,11 +4754,12 @@ ParametrizeByLen( const ICurveComposite* pCurve, DBLVECTOR& vParam) int nSpanU = pCurve->GetCurveCount() ; DBLVECTOR vLen ; double dLenTot = 0 ; + vParam.push_back( 0) ; for( int i = 0 ; i < nSpanU ; ++i) { const ICurve* pSubCrv = pCurve->GetCurve( i) ; double dLen ; pSubCrv->GetLength( dLen) ; dLenTot += dLen ; - vLen.push_back( dLen) ; + vLen.push_back( dLenTot) ; } // determino il parametro di ogni curva rispetto alla lunghezza totale for ( int i = 0 ; i < nSpanU ; ++i) @@ -4771,53 +4772,49 @@ ParametrizeByLen( const ICurveComposite* pCurve, DBLVECTOR& vParam) static bool BuildCommonParam( const DBLMATRIX& mParam, DBLVECTOR& vCommonParam) { - // aggiungo tutti gli end delle sottocurve - Intervals iInt ; - iInt.Set( 0, 1) ; - for ( int i = 0 ; i < int( mParam.size()) ; ++i) { - for ( int j = 0 ; j < int(mParam[i].size()) - 1 ; ++j) { - iInt.Add( mParam[i][j], mParam[i][j+1]) ; + vCommonParam = mParam[0] ; + for ( int i = 1 ; i < int( mParam.size()) ; ++i) { + for ( int j = 0 ; j < int( mParam[i].size()) ; ++j) + vCommonParam.push_back( mParam[i][j]) ; + } + //riordino ed elimino i doppioni entro una certa tolleranza + sort( vCommonParam.begin(), vCommonParam.end()) ; + for ( int i = 0 ; i < int(vCommonParam.size()) - 1 ; ++i) { + for( int j = i + 1 ; j < int(vCommonParam.size()) ; ++j) { + if ( vCommonParam[j] - vCommonParam[i] < EPS_SMALL){ + vCommonParam.erase(vCommonParam.begin() + j) ; + --j ; + } + else + break ; } } - //// NON FUNZIONANO COSì GLI INTERVALSS!!!!!!/////////////////////////// - // - //// controllo che non ce ne siano di troppo vicini - //// in tal caso elimino il secondo - //for ( ) { - // iInt. - //} - return true ; } -////---------------------------------------------------------------------------- -//bool -//SurfBezier::FindMatchByParam( const DBLVECTOR& vParam0, const DBLVECTOR& vParam1, INTVECTOR& vMatch0, INTVECTOR& vMatch1) const -//{ -// // per la modalità parametrata sulla lunghezza -// int nPoints0 = vParam0.size() ; -// int nPoints1 = vParam1.size() ; -// DBLVECTOR vSplit ; -// for ( int i = 0 ; i < nPoints0 - 1 ; ++i) { -// // valuto il limite dell'area di influenza dei punti della curva con meno punti -// // per i punti intermedi il limite è la metà tra i punti -// // per il primo e l'ultimo punto aumento il peso di quest'ultimi nel calcolo della media -// if ( i == 0) -// vSplit.push_back( (i * (1./3.) + ( i + 1) * ( 2./3.)) / ( nPoints0 - 1)) ; -// else if ( i == nPoints0 - 2) -// vSplit.push_back( (i * (2./3.) + ( i + 1) * ( 1./3.)) / ( nPoints0 - 1)) ; -// else -// vSplit.push_back( ( 2. * i + 1) / (2 * ( nPoints0 - 1))) ; -// } -// int nCount = 0 ; -// for ( double j = 0 ; j < nPoints1 ; ++j) { -// if ( nCount < int(vSplit.size()) && j / ( nPoints1 - 1) > vSplit[nCount] ) -// ++nCount ; -// vMatch0.push_back( nCount) ; -// } -// return true ; -//} +//---------------------------------------------------------------------------- +static bool +SplitByCommonParam( ICURVEPOVECTOR& vCrvBezUnif, DBLVECTOR& vCommonParam, DBLMATRIX& mParam) +{ + for ( int i = 0 ; i < int( vCrvBezUnif.size()) ; ++i) { + ICurveComposite* pCC = GetCurveComposite(vCrvBezUnif[i]) ; + int c = mParam[i].size() - 1 ; + DBLVECTOR vParam = mParam[i] ; + for( int j = vCommonParam.size() - 1 ; j >= 0 ; --j) { + // capisco su quale sottocurva devo fare lo split e riconverto il parametro rispetto al numero di sottocurve + while ( vCommonParam[j] < vParam[c]) + --c ; + if( vCommonParam[j] - vParam[c] > EPS_SMALL) { + double dSplit = (vCommonParam[j] - vParam[c]) / (vParam[c+1] - vParam[c]) + c ; + pCC->AddJoint( dSplit) ; + vParam.insert( vParam.begin() + c + 1, vCommonParam[j]) ; + ++c ; + } + } + } + return true ; +} static bool ChangeStartForClosed( PolyLine& plU0, PolyLine& plU1, ICurveComposite* pCrvU0, ICurveComposite* pCrvU1) { @@ -4869,22 +4866,53 @@ ChangeStartForClosed( PolyLine& plU0, PolyLine& plU1, ICurveComposite* pCrvU0, I bool SurfBezier::CreateBySetOfCurves( const ICURVEPOVECTOR& vCrvBez) { - //uniformo le curve e determino il grado e il numero di span condiviso - //... - //... + // uniformo le curve e determino il grado e il numero di span condiviso + bool bRat = false ; int nDegU = 3 ; - int nSpanU = 2 ; + ICURVEPOVECTOR vCrvBezUnif ; + DBLMATRIX mParam ; + // scorro le curve e uniformo tutte in curve di grado 3 non razionale + for ( int i = 0 ; i < int( vCrvBez.size()) ; ++i) { + const ICurveComposite* pCCFromSet = GetCurveComposite( vCrvBez[i]) ; + if( pCCFromSet == nullptr) + return false ; + PtrOwner pCC( CreateCurveComposite()) ; + for( int j = 0 ; j < int( pCCFromSet->GetCurveCount()) ; ++j) { + const ICurveBezier* pCrvBez = GetCurveBezier( pCCFromSet->GetCurve( j)) ; + if( pCrvBez == nullptr) + return false ; + if( ! pCC->AddCurve( EditBezierCurve( pCrvBez, nDegU, bRat))) + return false ; + } + mParam.emplace_back() ; + ParametrizeByLen( pCC, mParam.back()) ; + vCrvBezUnif.emplace_back( Release( pCC)) ; + } + // ora vado a splittare le curve in modo da avere una divisione condivisa + DBLVECTOR vdCommonPar ; + BuildCommonParam( mParam, vdCommonPar) ; + SplitByCommonParam( vCrvBezUnif, vdCommonPar, mParam) ; + int nSpanU = GetCurveComposite( vCrvBezUnif[0])->GetCurveCount() ; + + //// debug + //int nDegU = 3 ; + //bool bRat = false ; + //int nSpanU = 2 ; + //ICURVEPOVECTOR vCrvBezUnif ; + //for ( int i= 0 ; i < int( vCrvBez.size()) ; ++i) + // vCrvBezUnif.emplace_back( vCrvBez[i]->Clone()) ; + //// debug // calcolo i punti di controllo in V int nDegV = 3 ; - int nSpanV = int( vCrvBez.size()) - 1 ; + int nSpanV = int( vCrvBezUnif.size()) - 1 ; PNTMATRIX vPntCrvs ; for ( int j = 0 ; j < nSpanU ; ++j ) { PNTVECTOR vPntCtrl0 ; PNTVECTOR vPntCtrl1 ; - for ( int i = 0 ; i < int( vCrvBez.size()) ; ++i) { - const ICurveComposite* pCC = GetCurveComposite(vCrvBez[i]) ; + for ( int i = 0 ; i < int( vCrvBezUnif.size()) ; ++i) { + const ICurveComposite* pCC = GetCurveComposite(vCrvBezUnif[i]) ; const ICurveBezier* pCrvBez = GetCurveBezier(pCC->GetCurve( j)) ; if ( j == 0) vPntCtrl0.push_back( pCrvBez->GetControlPoint(0)) ; @@ -4906,9 +4934,9 @@ SurfBezier::CreateBySetOfCurves( const ICURVEPOVECTOR& vCrvBez) Vector3d vtDirXGeneral = ptMean - vPntCrvs[s][0] ; // prendo le curve a gruppi di 3 per costruire la parabola per trovare la pendenza "intuitiva" della superficie in V for ( int g = 0 ; g < nSpanV - 1 ; ++g) { - const ICurveComposite* pCC0 = GetCurveComposite( vCrvBez[g]) ; - const ICurveComposite* pCC1 = GetCurveComposite( vCrvBez[g + 1]) ; - const ICurveComposite* pCC2 = GetCurveComposite( vCrvBez[g + 2]) ; + const ICurveComposite* pCC0 = GetCurveComposite( vCrvBezUnif[g]) ; + const ICurveComposite* pCC1 = GetCurveComposite( vCrvBezUnif[g + 1]) ; + const ICurveComposite* pCC2 = GetCurveComposite( vCrvBezUnif[g + 2]) ; // su ognuna di queste 3 curve prendo un punto per ogni punto di controllo ( semplicemente dividendo uniformemente il parametrico) for ( int n = s == 0 ? 0 : 1 ; n < nDegU + 1 ; ++n) { const ICurveBezier* pCrv0 = GetCurveBezier( pCC0->GetCurve( s)) ; @@ -4967,7 +4995,7 @@ SurfBezier::CreateBySetOfCurves( const ICURVEPOVECTOR& vCrvBez) mA.col(0) << pow(pt0.x,2), pow(pt1.x,2) , pow(pt2.x,2) ; mA.col(1) << pt0.x, pt1.x , pt2.x ; mA.col(2) << 1, 1, 1 ; - if( abs( mA.determinant()) < EPS_SMALL) + if( abs( mA.determinant()) < EPS_SMALL) return false ; Eigen::Vector3d b ( pt0.y, pt1.y, pt2.y) ; Eigen::Vector3d coeff = mA.fullPivLu().solve(b) ; diff --git a/Tree.cpp b/Tree.cpp index 4a8bb23..2388616 100644 --- a/Tree.cpp +++ b/Tree.cpp @@ -801,7 +801,7 @@ Tree::BuildTree( double dLinTol, double dSideMin, double dSideMax) Plane3d plPlane ; double dArea = 0 ; bFlat = plCell.IsClosedAndFlat( plPlane, dArea, 10 * EPS_SMALL) ; } - if( ! bFlat) { + if( ! bFlat && dDist > 5 * EPS_SMALL) { bTwist = true ; // devo decidere in quale direzione splittare // dovrei capire in quale delle due direzioni è più torta la superficie