//---------------------------------------------------------------------------- // EgalTech 2013-2013 //---------------------------------------------------------------------------- // File : CurveAux.cpp Data : 22.11.13 Versione : 1.3a1 // Contenuto : Implementazione di alcune funzioni di utilità per le curve. // // // // Modifiche : 22.11.13 DS Creazione modulo. // // //---------------------------------------------------------------------------- //--------------------------- Include ---------------------------------------- #include "stdafx.h" #include "CurveAux.h" #include "GeoConst.h" #include "CurveLine.h" #include "CurveArc.h" #include "CurveBezier.h" #include "CurveComposite.h" #include "Voronoi.h" #include "IntersLineLine.h" #include "/EgtDev/Include/EGkDistPointCurve.h" #include "/EgtDev/Include/EGkStringUtils3d.h" #include "/EgtDev/Include/EGkUiUnits.h" #include "/EgtDev/Include/EgtPointerOwner.h" #include "/EgtDev/Include/EGkCurveByInterp.h" #define EIGEN_NO_IO #include "/EgtDev/Extern/Eigen/Dense" using namespace std ; static bool FindSpan( double dU, int nDeg, const DBLVECTOR& vKnots, int& nSpan) ; static bool CalcBasisFunc( double dU, int nSpan, int nDeg, const DBLVECTOR& vKnots, DBLVECTOR& vBasis) ; //---------------------------------------------------------------------------- bool IsClosed( const ICurve& crvC) { Point3d ptStart, ptEnd ; return ( crvC.GetStartPoint( ptStart) && crvC.GetEndPoint( ptEnd) && AreSamePointApprox( ptStart, ptEnd)) ; } //---------------------------------------------------------------------------- bool IsValidParam( const ICurve& crvC, double dPar, ICurve::Side nSide) { // recupero l'intervallo di validità del parametro double dStart, dEnd ; if ( ! crvC.GetDomain( dStart, dEnd)) return false ; // il parametro deve stare nei limiti if ( dPar < dStart - EPS_PARAM || dPar > dEnd + EPS_PARAM) return false ; // se curva chiusa non sono necessari altri controlli if ( crvC.IsClosed()) return true ; // per curve aperte non posso studiare la curva prima dell'inizio e dopo la fine if ( ( dPar < dStart + EPS_PARAM && nSide == ICurve::FROM_MINUS) || ( dPar > dEnd - EPS_PARAM && nSide == ICurve::FROM_PLUS)) return false ; return true ; } //---------------------------------------------------------------------------- bool IsStartParam( const ICurve& crvC, double dPar) { // recupero l'intervallo di validità del parametro double dStart, dEnd ; if ( ! crvC.GetDomain( dStart, dEnd)) return false ; // se il parametro non è nell'intorno dell'inizio if ( abs( dPar - dStart) > EPS_PARAM) return false ; return true ; } //---------------------------------------------------------------------------- bool IsEndParam( const ICurve& crvC, double dPar) { // recupero l'intervallo di validità del parametro double dStart, dEnd ; if ( ! crvC.GetDomain( dStart, dEnd)) return false ; // se il parametro non è nell'intorno della fine if ( abs( dPar - dEnd) > EPS_PARAM) return false ; return true ; } //---------------------------------------------------------------------------- bool GetNearestExtremityToPoint( const Point3d& ptP, const ICurve& Curve, bool& bStart) { // recupero punto iniziale Point3d ptStart ; if ( ! Curve.GetStartPoint( ptStart)) return false ; // recupero punto finale Point3d ptEnd ; if ( ! Curve.GetEndPoint( ptEnd)) return false ; // se curva aperta if ( ! AreSamePointApprox( ptStart, ptEnd)) { // è il più vicino tra inizio e fine bStart = ( SqDist( ptP, ptStart) <= SqDist( ptP, ptEnd)) ; return true ; } // altrimenti curva chiusa, devo proiettare il punto sulla curva e misurare la distanza su di essa DistPointCurve distPC( ptP, Curve) ; int nFlag ; double dLen, dU, dULen ; if ( Curve.GetLength( dLen) && distPC.GetParamAtMinDistPoint( 0, dU, nFlag) && Curve.GetLengthAtParam( dU, dULen)) { bStart = ( dULen <= ( dLen - dULen)) ; return true ; } else return false ; } //---------------------------------------------------------------------------- bool MoveParamToAvoidTg( double& dU, ICurve::Side nSide, const ICurve& Curve) { // verifico che il parametro sia accettabile if ( ! Curve.IsValidParam( dU, nSide)) return false ; // determino un incremento piccolo ma valido del parametro U double dDeltaU = 1000 * EPS_PARAM ; Point3d ptPos ; Vector3d vtDer1 ; if ( Curve.GetPointD1D2( dU, nSide, ptPos, &vtDer1)) { double dDer1 = vtDer1.LenXY() ; if ( dDer1 > EPS_ZERO) dDeltaU = 2 * EPS_SMALL / dDer1 ; } // cerco di spostarmi per evitare eventuali problemi di tangenza if ( nSide == ICurve::FROM_MINUS) { double dUm = dU - dDeltaU ; if ( ! Curve.IsValidParam( dUm, nSide)) { if ( Curve.IsClosed()) { double dStart, dEnd ; Curve.GetDomain( dStart, dEnd) ; dUm = ( dEnd - dStart) - dDeltaU ; } else dUm = dU ; } dU = dUm ; return true ; } else { // nSide == ICurve::FROM_PLUS double dUm = dU + dDeltaU ; if ( ! Curve.IsValidParam( dUm, nSide)) { if ( Curve.IsClosed()) { double dStart, dEnd ; Curve.GetDomain( dStart, dEnd) ; dUm = dStart + dDeltaU ; } else dUm = dU ; } dU = dUm ; return true ; } } //---------------------------------------------------------------------------- bool GetTang( const ICurve& crvC, double dU, ICurve::Side nS, Vector3d& vtTang) { Point3d ptPos ; return GetPointTang( crvC, dU, nS, ptPos, vtTang) ; } //---------------------------------------------------------------------------- bool GetPointTang( const ICurve& crvC, double dU, ICurve::Side nS, Point3d& ptPos, Vector3d& vtTang) { if ( ! crvC.GetPointD1D2( dU, nS, ptPos, &vtTang)) return false ; if ( vtTang.Normalize( EPS_ZERO)) return true ; // nel caso la derivata prima sia nulla, utilizziamo la seconda Vector3d vtDummy ; if ( ! crvC.GetPointD1D2( dU, nS, ptPos, &vtDummy, &vtTang)) return false ; double dUmod = dU + ( nS == ICurve::FROM_MINUS ? -1 : +1) * 100 * EPS_PARAM ; Point3d ptDummy ; // se anche la derivata seconda è nulla, provo a spostarmi di poco if ( ! vtTang.Normalize( EPS_ZERO)) { if ( ! crvC.GetPointD1D2( dUmod, nS, ptDummy, &vtDummy, &vtTang)) return false ; if ( ! vtTang.Normalize( EPS_ZERO)) return false ; } // per il verso, lo confrontiamo con quello della derivata prima un poco discosta Vector3d vtD1near ; if ( ! crvC.GetPointD1D2( dUmod, nS, ptDummy, &vtD1near)) return false ; if ( ( vtTang * vtD1near) < 0) vtTang.Invert() ; return true ; } //---------------------------------------------------------------------------- bool GetPointDiffGeom( const ICurve& crvC, double dU, ICurve::Side nS, CrvPointDiffGeom& oDiffG) { // calcolo punto e derivate Point3d ptPos ; Vector3d vtDer1, vtDer2 ; if ( ! crvC.GetPointD1D2( dU, nS, ptPos, &vtDer1, &vtDer2)) return false ; // assegno parametro e posizione oDiffG.dU = dU ; oDiffG.ptP = ptPos ; // se esiste la derivata prima non nulla double dSqLenD1 = vtDer1.SqLen() ; if ( vtDer1.Normalize()) { oDiffG.vtT = vtDer1 ; // del vettore deriv2^ tengo la sola componente perpendicolare al vettore tangente oDiffG.vtN = vtDer2 - ( vtDer2 * vtDer1) * vtDer1 ; if ( oDiffG.vtN.Normalize()) { oDiffG.dCurv = ( vtDer1 ^ vtDer2).Len() / dSqLenD1 ; oDiffG.nStatus = CrvPointDiffGeom::NCRV ; } else { oDiffG.dCurv = 0 ; oDiffG.nStatus = CrvPointDiffGeom::TANG ; } } // se esiste almeno derivata seconda non nulla, definisce la tangente else if ( vtDer2.Normalize()) { oDiffG.vtT = vtDer2 ; oDiffG.vtN.Set( 0, 0, 0) ; oDiffG.dCurv = 0 ; // per il verso, lo confrontiamo con quello della derivata prima un poco discosta Point3d ptDummy ; Vector3d vtD1near ; dU += ( nS == ICurve::FROM_MINUS ? -1 : +1) * 100 * EPS_PARAM ; if ( crvC.GetPointD1D2( dU, nS, ptDummy, &vtD1near)) { if ( ( oDiffG.vtT * vtD1near) < 0) oDiffG.vtT.Invert() ; oDiffG.nStatus = CrvPointDiffGeom::TANG ; } else { oDiffG.vtT.Set( 0, 0, 0) ; oDiffG.nStatus = CrvPointDiffGeom::POS ; } } // altrimenti non sono definite la tangente e la normale else { oDiffG.vtT.Set( 0, 0, 0) ; oDiffG.vtN.Set( 0, 0, 0) ; oDiffG.dCurv = 0 ; oDiffG.nStatus = CrvPointDiffGeom::POS ; } return true ; } //---------------------------------------------------------------------------- bool ImproveCurveParamAtPoint( double& dU, const Point3d& ptP, const ICurve* pCrv) { // Da usare quando il parametro è già molto vicino a quello esatto // calcolo il punto della curva in corrispondenza al parametro Point3d ptQ ; Vector3d vtD ; if ( ! pCrv->GetPointD1D2( dU, ICurve::FROM_MINUS, ptQ, &vtD)) return false ; // se sono uguali, è già tutto ok if ( AreSamePointExact( ptP, ptQ)) return true ; // se derivata nulla, non posso migliorare if ( vtD.IsZero()) return false ; // incremento parametro su linea tangente (uso la derivata) double dDeltaU = ( ptP - ptQ) * vtD / vtD.SqLen() ; // se migliorato, tengo nuovo valore Point3d ptR ; if ( pCrv->GetPointD1D2( dU + dDeltaU, ICurve::FROM_MINUS, ptR) && ( ptP - ptR).SqLen() < ( ptP - ptQ).SqLen()) dU += dDeltaU ; return true ; } //---------------------------------------------------------------------------- bool CurveGetAreaXY( const ICurve& crvC, double& dArea) { // verifico sia chiusa if ( ! crvC.IsClosed()) return false ; double dAreaXY = 0 ; int nType = crvC.GetType() ; if ( nType == CRV_ARC) { // se circonferenza const CurveArc* pArc = GetBasicCurveArc( &crvC) ; if ( pArc == nullptr) return false ; double dAng = pArc->GetAngCenter() ; if ( abs( abs( dAng) - ANG_FULL) > EPS_ANG_SMALL) return false ; double dRad = pArc->GetRadius() ; Vector3d vtN = pArc->GetNormVersor() ; dAreaXY = PIGRECO * dRad * dRad * vtN * Z_AX * ( dAng > 0 ? 1 : -1) ; } else if ( nType == CRV_BEZIER) { // se Bezier approssimo la curva con una polilinea PolyLine PL ; crvC.ApproxWithLines( LIN_TOL_STD, ANG_TOL_STD_DEG, ICurve::APL_SPECIAL_INT, PL) ; PL.GetAreaXY( dAreaXY) ; } else if ( nType == CRV_COMPO) { // se composita calcolo il contributo di ogni sottotratto const CurveComposite* pCompo = GetBasicCurveComposite( &crvC) ; for ( int i = 0 ; i < pCompo->GetCurveCount() ; i++) { const ICurve* pCrv = pCompo->GetCurve( i) ; if ( pCrv == nullptr) return false ; if ( pCrv->GetType() == CRV_LINE) { Point3d ptS, ptE ; if ( ! pCrv->GetStartPoint( ptS) || ! pCrv->GetEndPoint( ptE)) return false ; dAreaXY += ( ptS.x - ptE.x) * ( ptS.y + ptE.y) ; } else if ( pCrv->GetType() == CRV_ARC) { const CurveArc* pArc = GetBasicCurveArc( pCrv) ; if ( pArc == nullptr) return false ; Point3d ptS, ptE ; if ( ! pArc->GetStartPoint( ptS) || ! pArc->GetEndPoint( ptE)) return false ; // recupero i dati dell'arco Point3d ptC = pArc->GetCenter() ; Vector3d vtN = pArc->GetNormVersor() ; double dRad = pArc->GetRadius() ; double dAng = pArc->GetAngCenter() ; dAreaXY += ( ptS.x - ptC.x) * ( ptS.y + ptC.y) + ( ptC.x - ptE.x) * ( ptE.y + ptC.y) + dRad * dRad * dAng * DEGTORAD * vtN * Z_AX ; } else if ( pCrv->GetType() == CRV_BEZIER) { // approssimo la sottocurva con polilinea PolyLine PL ; pCrv->ApproxWithLines( LIN_TOL_STD, ANG_TOL_STD_DEG, ICurve::APL_SPECIAL_INT, PL) ; Point3d ptS, ptE ; for ( bool bFound = PL.GetFirstLine( ptS, ptE) ; bFound ; bFound = PL.GetNextLine( ptS, ptE)) { dAreaXY += ( ptS.x - ptE.x) * ( ptS.y + ptE.y) ; } } } dAreaXY *= 0.5 ; } else return false ; // restituisco il valore if ( &dArea != nullptr) dArea = dAreaXY ; return true ; } //---------------------------------------------------------------------------- bool CurveGetArea( const ICurve& crvC, Plane3d& plPlane, double& dArea) { // verifico sia chiusa if ( ! crvC.IsClosed()) return false ; // approssimo la curva con una polilinea PolyLine PL ; crvC.ApproxWithLines( LIN_TOL_STD, ANG_TOL_STD_DEG, ICurve::APL_SPECIAL_INT, PL) ; // calcolo l'area Plane3d plMyPlane ; double dMyArea = 0 ; if ( ! PL.IsClosedAndFlat( plMyPlane, dMyArea, 50 * EPS_SMALL)) return false ; // restituisco i valori if ( &plPlane != nullptr) plPlane = plMyPlane ; if ( &dArea != nullptr) dArea = dMyArea ; return true ; } //---------------------------------------------------------------------------- bool CurveDump( const ICurve& crvC, string& sOut, bool bMM, const char* szNewLine) { // verifico validità curva if ( ! crvC.IsValid()) return false ; // punti iniziale e finale Point3d ptPS ; crvC.GetStartPoint( ptPS) ; sOut += "PS(" + ToString( GetInUiUnits( ptPS, bMM), 3) + ")" + szNewLine ; Point3d ptPE ; crvC.GetEndPoint( ptPE) ; sOut += "PE(" + ToString( GetInUiUnits( ptPE, bMM), 3) + ")" + szNewLine ; // direzioni iniziale e finale Vector3d vtDirS ; crvC.GetStartDir( vtDirS) ; sOut += "VS(" + ToString( vtDirS, 3) + ")" + szNewLine ; Vector3d vtDirE ; crvC.GetEndDir( vtDirE) ; sOut += "VE(" + ToString( vtDirE, 3) + ")" + szNewLine ; // eventuali direzione di estrusione e spessore Vector3d vtExtr ; crvC.GetExtrusion( vtExtr) ; if ( ! vtExtr.IsSmall()) { double dThick ; crvC.GetThickness( dThick) ; sOut += "Extr(" + ToString( vtExtr, 3) + ") Th=" + ToString( GetInUiUnits( dThick, bMM), 3) + szNewLine ; } // lunghezza double dLen = 0 ; crvC.GetLength( dLen) ; sOut += "Len=" + ToString( GetInUiUnits( dLen, bMM), 3) + szNewLine ; // altri dati per curva chiusa double dAreaXY ; if ( CurveGetAreaXY( crvC, dAreaXY)) { bool bCCW = ( dAreaXY > 0) ; double dAreaUi = GetAreaInUiUnits( abs( dAreaXY), bMM) ; int nDec = ( dAreaUi > 100 ? 1 : ( dAreaUi > 0.1 ? 3 : 6)) ; sOut += string( "Closed") + ( bCCW ? " CCW" : " CW") + " AreaXY=" + ToString( dAreaUi, nDec) + szNewLine ; } // eventuali proprietà temporanee vector vProp{ crvC.GetTempProp( 0), crvC.GetTempProp( 1)} ; if ( vProp[0] != 0 || vProp[1] != 0) sOut += string( "TmpProp=") + ToString( vProp) + szNewLine ; // eventuali parametri temporanei vector vParam{ crvC.GetTempParam( 0), crvC.GetTempParam( 1)} ; if ( abs( vParam[0]) > EPS_ZERO || abs( vParam[1]) > EPS_ZERO) sOut += string( "TmpParam=") + ToString( vParam) + szNewLine ; return true ; } //---------------------------------------------------------------------------- bool CopyExtrusion( const ICurve* pSouCrv, ICurve* pDestCrv) { // verifico puntatori if ( pSouCrv == nullptr || pDestCrv == nullptr) return false ; // eseguo copia Vector3d vtExtr ; pSouCrv->GetExtrusion( vtExtr) ; pDestCrv->SetExtrusion( vtExtr) ; return true ; } //---------------------------------------------------------------------------- bool CopyThickness( const ICurve* pSouCrv, ICurve* pDestCrv) { // verifico puntatori if ( pSouCrv == nullptr || pDestCrv == nullptr) return false ; // eseguo copia double dThick = 0 ; pSouCrv->GetThickness( dThick) ; pDestCrv->SetThickness( dThick) ; return true ; } //---------------------------------------------------------------------------- ICurve* CurveToBezierCurve( const ICurve* pCrv, int nDeg, bool bMakeRatOrNot) { PtrOwner pBezierForm ; switch ( pCrv->GetType()) { case CRV_LINE :{ const ICurveLine* pCrvLine = static_cast( pCrv) ; int nDegWanted = nDeg == -1 ? 3 : nDeg ; pBezierForm.Set( LineToBezierCurve( pCrvLine, nDegWanted, bMakeRatOrNot)) ; break ; } case CRV_ARC : { const ICurveArc* pCrvArc = static_cast( pCrv) ; int nDegWanted = nDeg == -1 ? 3 : nDeg ; pBezierForm.Set( ArcToBezierCurve( pCrvArc, nDegWanted, bMakeRatOrNot)) ; break ; } case CRV_COMPO : { const ICurveComposite* pCrvCompo = static_cast( pCrv) ; pBezierForm.Set( CompositeToBezierCurve( pCrvCompo, nDeg, bMakeRatOrNot)) ; break ; } case CRV_BEZIER : { const ICurveBezier* pCrvBezier = static_cast( pCrv) ; pBezierForm.Set( EditBezierCurve( pCrvBezier, nDeg, bMakeRatOrNot)) ; break ; } default : return nullptr ; break ; } return Release( pBezierForm) ; } //---------------------------------------------------------------------------- ICurveBezier* LineToBezierCurve( const ICurveLine* pCrvLine, int nDeg, bool bMakeRatOrNot) { // verifico che esista la linea if ( pCrvLine == nullptr) return nullptr ; PtrOwner pCrvBezier( CreateCurveBezier()) ; // rendo tutte le curve di grado 2 e razionali così posso convertire anche archi e avere tutte curve dello stesso grado e razionali pCrvBezier->Init( nDeg, true) ; pCrvBezier->FromLine( *pCrvLine) ; if ( bMakeRatOrNot) pCrvBezier->MakeRational() ; return Release( pCrvBezier) ; } //---------------------------------------------------------------------------- ICurve* ArcToBezierCurve( const ICurveArc* pArc, int nDeg, bool bMakeRatOrNot) { // una spirale non può essere forzata al grado 2 // verifico che esista l'arco if ( pArc == nullptr) return nullptr ; // se angolo al centro sotto il limite, basta una curva if ( abs( pArc->GetAngCenter()) < BEZARC_ANG_CEN_MAX + EPS_ANG_SMALL) { // creo la curva di Bezier equivalente all'arco PtrOwner pCrvBez( CreateBasicCurveBezier()) ; if ( IsNull( pCrvBez) || ! pCrvBez->FromArc( *pArc)) return nullptr ; if ( ! bMakeRatOrNot) { Point3d ptCen = pArc->GetCenter() ; Vector3d vtN = pArc->GetNormVersor() ; pCrvBez.Set( ApproxArcCurveBezierWithSingleCubic( pCrvBez, ptCen, vtN)) ; } if ( IsNull( pCrvBez)) return nullptr ; // aumento il grado della curva come richiesto while ( pCrvBez->GetDegree() < nDeg) pCrvBez.Set( BezierIncreaseDegree( pCrvBez)) ; // restituisco la curva return Release( pCrvBez) ; } // altrimenti curva composita di Bezier else { // creo la curva composita PtrOwner pCrvCompo( CreateBasicCurveComposite()) ; if ( IsNull( pCrvCompo)) return nullptr ; // inserisco nella CC le curve di Bezier equivalenti alle parti dell'arco int nParts = (int) ceil( abs( pArc->GetAngCenter()) / ( BEZARC_ANG_CEN_MAX + EPS_ANG_SMALL)) ; nParts = max( nParts, 2) ; for ( int i = 0 ; i < nParts ; ++ i) { // copio l'arco originale CurveArc cArc = *GetBasicCurveArc( pArc) ; // lo limito alla parte di interesse cArc.TrimStartEndAtParam( double( i) / nParts, double( i + 1) / nParts) ; // creo la curva di Bezier equivalente PtrOwner pCrvBez( CreateBasicCurveBezier()) ; if ( IsNull( pCrvBez) || ! pCrvBez->FromArc( cArc)) return nullptr ; // aumento il grado della curva come richiesto while ( pCrvBez->GetDegree() < nDeg) pCrvBez.Set( BezierIncreaseDegree( pCrvBez)) ; // aggiungo la curva di Bezier a quella composita if ( ! pCrvCompo->AddCurve( Release( pCrvBez))) return nullptr ; } // copio estrusione e spessore CopyExtrusion( pArc, pCrvCompo) ; CopyThickness( pArc, pCrvCompo) ; // restituisco la curva return Release( pCrvCompo) ; } } //---------------------------------------------------------------------------- ICurve* CompositeToBezierCurve( const ICurveComposite* pCC, int nDeg, bool bMakeRatOrNot) { // verifico che esista la curva if ( pCC == nullptr) return nullptr ; // converto tutte le curve in bezier razionali di grado 2 PtrOwner pCCBezier( CreateCurveComposite()) ; for ( int i = 0 ; i < int( pCC->GetCurveCount()) ; ++i) { PtrOwner pCrvNew ; if ( pCC->GetCurve(i)->GetType() == CRV_ARC) { const CurveArc* crArc = GetBasicCurveArc( pCC->GetCurve(i)) ; ICurve* pCrvBezier = ArcToBezierCurve( crArc) ; if ( pCrvBezier == nullptr) return nullptr ; // se la curva è di grado superiore al secondo allora devo ricondurla al secondo grado pCrvNew.Set( pCrvBezier) ; } else if ( pCC->GetCurve(i)->GetType() == CRV_LINE) { const CurveLine* crLine = GetBasicCurveLine( pCC->GetCurve(i)) ; ICurve* pCrvBezier = LineToBezierCurve( crLine, nDeg, bMakeRatOrNot) ; if ( pCrvBezier == nullptr) return nullptr ; pCrvNew.Set( pCrvBezier) ; } else if ( pCC->GetCurve(i)->GetType() == CRV_BEZIER ) { const CurveBezier* crvBezier = GetBasicCurveBezier( pCC->GetCurve(i)) ; ICurve* pCrvBezier = EditBezierCurve( crvBezier, nDeg, bMakeRatOrNot) ; if ( pCrvBezier == nullptr) return nullptr ; pCrvNew.Set( pCrvBezier) ; } pCCBezier->AddCurve( Release( pCrvNew)) ; } return Release( pCCBezier) ; } //---------------------------------------------------------------------------- ICurve* EditBezierCurve( const ICurveBezier* pCrvBezier, int nDeg, bool bMakeRatOrNot, double dTol) { // se nDeg == -1 allora viene mantenuto il grado della curva originale // verifico sia una bezier if ( pCrvBezier == nullptr) return nullptr ; if ( nDeg == 2 || nDeg == 1) return nullptr ; PtrOwner pCrvNew( pCrvBezier->Clone()) ; int nDegCurr = pCrvNew->GetDegree() ; bool bRat = pCrvNew->IsRational() ; int nDegWanted = nDeg == -1 ? nDegCurr : nDeg ; if ( ! bMakeRatOrNot) { if ( ! pCrvNew->MakeNonRational( dTol)) { // 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)), nDegWanted, bMakeRatOrNot, dTol)) ) return nullptr ; } return Release( pCCEdited) ; } } // se la curva è già nella forma giusta la restituisco if ( nDegCurr == nDegWanted && ( ( ! bRat && ! bMakeRatOrNot) || (bRat && bMakeRatOrNot))) return Release( pCrvNew) ; // sennò mi porto al grado giusto if ( nDegCurr < nDegWanted) { while ( nDegCurr < nDegWanted) { pCrvNew.Set( BezierIncreaseDegree( pCrvNew)) ; ++ nDegCurr ; } } else if ( nDegCurr > nDegWanted) { while ( nDegCurr > nDegWanted) { 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 ; } } if ( bMakeRatOrNot) pCrvNew->MakeRational() ; return Release( pCrvNew) ; } //---------------------------------------------------------------------------- ICurveBezier* BezierIncreaseDegree( const ICurveBezier* pCrvBezier) { if ( pCrvBezier == nullptr) return nullptr ; // creo la versione con grado aumentato PtrOwner pNewBezier( CreateCurveBezier()) ; if ( IsNull( pNewBezier)) return nullptr ; int nDeg = pCrvBezier->GetDegree() + 1 ; bool bRat = pCrvBezier->IsRational() ; pNewBezier->Init( nDeg, bRat) ; // prev e curr sono riferiti alla curva di partenza // salvo il primo punto Point3d ptCtrlPrev = pCrvBezier->GetControlPoint( 0) ; double dWprev = 1 ; if ( bRat) { dWprev = pCrvBezier->GetControlWeight( 0) ; pNewBezier->SetControlPoint( 0, ptCtrlPrev, dWprev) ; } else pNewBezier->SetControlPoint( 0, ptCtrlPrev) ; // ciclo sui punti di controllo intermedi per calcolare quelli nuovi Point3d ptCtrlCurr ; double dWcurr ; for ( int i = 1 ; i < nDeg ; ++i) { ptCtrlCurr = pCrvBezier->GetControlPoint( i) ; double dAlpha = double( i) / nDeg ; if ( bRat) { dWcurr = pCrvBezier->GetControlWeight( i) ; double dWnew = dAlpha * dWprev + ( 1 - dAlpha) * dWcurr ; Point3d ptNew = dAlpha * ptCtrlPrev * dWprev + ( 1 - dAlpha) * ptCtrlCurr * dWcurr; ptNew /= dWnew ; pNewBezier->SetControlPoint( i, ptNew, dWnew) ; dWprev = dWcurr ; } else { Point3d ptNew = dAlpha * ptCtrlPrev + ( 1 - dAlpha) * ptCtrlCurr ; pNewBezier->SetControlPoint( i, ptNew) ; } ptCtrlPrev = ptCtrlCurr ; } // salvo l'ultimo punto ptCtrlCurr = pCrvBezier->GetControlPoint( nDeg - 1) ; if ( bRat) { dWcurr = pCrvBezier->GetControlWeight( nDeg - 1) ; pNewBezier->SetControlPoint( nDeg, ptCtrlCurr, dWcurr) ; } else pNewBezier->SetControlPoint( nDeg, ptCtrlCurr) ; return Release( pNewBezier) ; } //---------------------------------------------------------------------------- ICurveBezier* BezierDecreaseDegree(const ICurveBezier* pCrvBezier, double dTol) { if ( pCrvBezier == nullptr) return nullptr ; PtrOwner pNewBezier( CreateCurveBezier()) ; int nDeg = pCrvBezier->GetDegree() - 1 ; if ( nDeg == 0) return nullptr ; bool bRat = pCrvBezier->IsRational() ; pNewBezier->Init( nDeg, bRat) ; // prev e curr sono riferiti alla nuova curva //salvo il primo punto Point3d ptCtrlPrev = pCrvBezier->GetControlPoint( 0) ; double dWprev = 1 ; if ( bRat ) { dWprev = pCrvBezier->GetControlWeight( 0) ; pNewBezier->SetControlPoint( 0, ptCtrlPrev, dWprev) ; } else pNewBezier->SetControlPoint( 0, ptCtrlPrev) ; // ciclo sui punti intermedi per trovare i nuovi punti di controllo // equazioni da NURBS book Point3d ptCtrlCurr ; double dWcurr ; int r = int (nDeg / 2) ; // prima metà for ( double i = 1 ; i < r ; ++i) { double dAlpha = i / ( nDeg + 1) ; // old è riferito alla curva originale Point3d ptOld = pCrvBezier->GetControlPoint( int( i)) ; if ( bRat) { double dWold = pCrvBezier->GetControlWeight( int( i)) ; dWcurr = ( dWold + (- dAlpha * dWprev)) / ( 1 - dAlpha) ; ptCtrlCurr = ( ptOld * dWold + (- dAlpha * ptCtrlPrev * dWprev)) / ( 1 - dAlpha) ; ptCtrlCurr /= dWcurr ; pNewBezier->SetControlPoint( int( i), ptCtrlCurr, dWcurr) ; dWprev = dWcurr ; } else { ptCtrlCurr = ( ptOld + (- dAlpha * ptCtrlPrev)) / ( 1 - dAlpha) ; pNewBezier->SetControlPoint( int( i), ptCtrlCurr) ; } ptCtrlPrev = ptCtrlCurr ; } // risolvo il punto di mezzo per il caso nDeg pari if ( ( nDeg + 1) % 2 == 0) { // pari double dAlpha = r / ( nDeg + 1.) ; Point3d ptOld = pCrvBezier->GetControlPoint( r) ; if ( bRat) { double dWold = pCrvBezier->GetControlWeight( r) ; dWcurr = ( dWold + (- dAlpha * dWprev)) / ( 1 - dAlpha) ; ptCtrlCurr = ( ptOld * dWold + (- dAlpha * ptCtrlPrev * dWprev)) / ( 1 - dAlpha) ; ptCtrlCurr /= dWcurr ; pNewBezier->SetControlPoint( r, ptCtrlCurr, dWcurr) ; dWprev = dWcurr ; } else { ptCtrlCurr = ( ptOld + (- dAlpha * ptCtrlPrev)) / ( 1 - dAlpha) ; pNewBezier->SetControlPoint( r, ptCtrlCurr) ; } } // salvo l'ultimo punto ptCtrlCurr = pCrvBezier->GetControlPoint( nDeg + 1) ; if ( bRat ) { dWcurr = pCrvBezier->GetControlWeight( nDeg + 1) ; pNewBezier->SetControlPoint( nDeg, ptCtrlCurr, dWcurr) ; } else pNewBezier->SetControlPoint( nDeg, ptCtrlCurr) ; ptCtrlPrev = ptCtrlCurr ; if ( bRat) dWprev = dWcurr ; // seconda metà for ( double i = nDeg - 1 ; i > r ; --i) { double dAlpha = ( i + 1) / ( nDeg + 1) ; Point3d ptOld = pCrvBezier->GetControlPoint( int( i) + 1) ; if ( bRat) { double dWold = pCrvBezier->GetControlWeight( int( i) + 1) ; dWcurr = ( dWold - ( ( 1 - dAlpha) * dWprev)) / dAlpha ; ptCtrlCurr = ( ptOld * dWold + ( - ( 1 - dAlpha) * ptCtrlPrev * dWprev)) / dAlpha ; ptCtrlCurr /= dWcurr ; pNewBezier->SetControlPoint( int( i), ptCtrlCurr, dWcurr) ; dWprev = dWcurr ; } else { ptCtrlCurr = ( ptOld + ( - ( 1 - dAlpha) * ptCtrlPrev)) / dAlpha ; pNewBezier->SetControlPoint( int( i), ptCtrlCurr) ; } ptCtrlPrev = ptCtrlCurr ; } // risolvo il punto di mezzo per il caso nDeg dispari // calcolo anche l'errore di approssimazione double dErr = 0 ; if ( (nDeg + 1) % 2 != 0) { // dispari // puntdo di sinistra double dAlpha = r / ( nDeg + 1.) ; ptCtrlPrev = pNewBezier->GetControlPoint( r - 1) ; Point3d ptOld = pCrvBezier->GetControlPoint( r) ; Point3d ptL ; double dWL = 1 ; if ( bRat) { dWprev = pNewBezier->GetControlWeight( r - 1) ; double dWold =pCrvBezier->GetControlWeight( r) ; dWL = ( dWold - ( dAlpha * dWprev)) / ( 1 - dAlpha) ; ptL = ( ptOld * dWold + ( - dAlpha * ptCtrlPrev * dWprev)) / ( 1 - dAlpha) ; } else ptL = ( ptOld + ( - dAlpha * ptCtrlPrev)) / ( 1 - dAlpha) ; // punto di destra dAlpha = ( r + 1.) / ( nDeg + 1.) ; ptCtrlPrev = pNewBezier->GetControlPoint( r + 1) ; ptOld = pCrvBezier->GetControlPoint( r + 1) ; double dWR = 1 ; Point3d ptR ; if ( bRat) { dWprev = pNewBezier->GetControlWeight( r + 1) ; double dWold =pCrvBezier->GetControlWeight( r + 1) ; dWR = ( dWold - ( ( 1 - dAlpha) * dWprev)) / dAlpha ; ptR = ( ptOld * dWold + ( - ( 1 - dAlpha) * ptCtrlPrev * dWprev)) / dAlpha ; } else ptR = ( ptOld + ( - ( 1 - dAlpha) * ptCtrlPrev)) / dAlpha ; // calcolo il punto di mezzo e lo aggiungo Point3d ptNew = ( ptL + ptR) / 2 ; if ( bRat) { double dWnew = ( dWL + dWR) / 2 ; ptNew /= dWnew ; pNewBezier->SetControlPoint( r, ptNew, dWnew) ; } else pNewBezier->SetControlPoint( r, ptNew) ; dErr = Dist( ptL, ptR) ; } else { ptCtrlCurr = pNewBezier->GetControlPoint( r + 1) ; ptCtrlPrev = pNewBezier->GetControlPoint( r) ; Point3d ptOld = pCrvBezier->GetControlPoint( r + 1) ; dErr = Dist( ptOld, 0.5 * ( ptCtrlPrev + ptCtrlCurr)) ; } // se l'approssimazione dà un errore troppo alto allora annullo tutto // ricalcolo l'errore a mano, per avere un valore più attendibile CalcApproxError( pCrvBezier, pNewBezier, dErr) ; if ( dErr > dTol) return nullptr ; return Release( pNewBezier) ; } //---------------------------------------------------------------------------- ICurve* ApproxBezierWithCubics(const ICurve* pCrv, double dTol) { // verifico sia una bezier const CurveBezier* pCrvBezier = GetBasicCurveBezier( pCrv) ; if ( pCrvBezier == nullptr) return nullptr ; // cerco di stimare quanti cambi di concavità ho // tiro una linea tra il primo punto di controllo e l'ultimo e poi scorro gli altri punti di controllo // controllando quante volte salto da un lato all'altro della linea PtrOwner pCL( CreateCurveLine()) ; Point3d ptStart, ptEnd ; pCrvBezier->GetStartPoint( ptStart) ; pCrvBezier->GetEndPoint( ptEnd) ; pCL->Set( ptStart, ptEnd) ; Point3d ptCtrl = pCrvBezier->GetControlPoint( 1) ; Vector3d vtN = (ptCtrl - ptStart) ^ ( ptCtrl - ptEnd) ; // controllo se la curva di bezier è un arco perfetto // in questo caso userò i numeri approssimati int nDeg = pCrvBezier->GetDegree() ; bool bRat = pCrvBezier->IsRational() ; bool bIsArc = false ; Point3d ptCen ; if ( nDeg == 2 && bRat) { Vector3d vtStartDir, vtEndDir ; pCrvBezier->GetPointD1D2( 0, ICurve::FROM_MINUS, ptStart, &vtStartDir) ; pCrvBezier->GetPointD1D2( 1, ICurve::FROM_PLUS, ptEnd, &vtEndDir) ; vtStartDir.Rotate( vtN, 90) ; vtEndDir.Rotate( vtN, 90) ; Point3d ptAux1 = ptStart + vtStartDir ; Point3d ptAux2 = ptEnd + vtEndDir ; PtrOwner pCL1( CreateBasicCurveLine()) ; pCL1->Set( ptStart, ptAux1) ; PtrOwner pCL2( CreateBasicCurveLine()) ; pCL2->Set( ptEnd, ptAux2) ; IntersLineLine ill( *pCL1, *pCL2, false) ; if ( ill.GetNumInters() != 0) { IntCrvCrvInfo iccInfo ; ill.GetIntCrvCrvInfo( iccInfo) ; ptCen = iccInfo.IciA[0].ptI ; double dDist1 = Dist( ptStart, ptCen) ; double dDist2 = Dist( ptEnd, ptCen) ; if ( abs(dDist1 - dDist2) < EPS_SMALL ) { CurveArc cArc ; cArc.SetC2PN( ptCen, ptStart, ptEnd, vtN) ; // controllo se il raggio della circonferenza risultante coincide con quello calcolato prima // impongo inoltre che l'angolo al centro sia minore di 90 gradi if ( abs( cArc.GetRadius() - dDist1) < EPS_SMALL && cArc.GetAngCenter() < 90 + EPS_SMALL) bIsArc = true ; } } } int nSidePrec = -2 ; int nCross = 0 ; // se non ho un arco cerco di capire se ho cambi di concavità if ( ! bIsArc) { for ( int i = 1 ; i < pCrvBezier->GetDegree() - 1; ++i) { ptCtrl = pCrvBezier->GetControlPoint( i) ; DistPointCurve dpc( ptCtrl, *pCL) ; int nSide = -2 ; dpc.GetSideAtMinDistPoint( 0, vtN, nSide) ; if ( i > 0 && nSide * nSidePrec < 1) ++nCross ; if ( nSide != 0) nSidePrec = nSide ; } } // determino in quanti tratti spezzare la curva originale guardando il grado e i cambi di concavità int nParts = nCross > nDeg / 2 ? nCross : nDeg / 2 + 1 ; bool bOneIsEnough = false ; if ((( nDeg == 2 || nDeg == 3) && ! bRat) || bIsArc) { nParts = 1 ; bOneIsEnough = true ; } double dErr = INFINITO ; PtrOwner pCC( CreateBasicCurveComposite()) ; // approssimo l'ultimo tratto int nCount = 0 ; while ( dErr > dTol && nCount < 100) { PtrOwner pCrvPart( pCrvBezier->Clone()) ; pCrvPart->TrimStartEndAtParam( double( nParts - 1) / nParts, 1) ; PtrOwner pCrvCubic ; if ( ! bIsArc) pCrvCubic.Set( ApproxCurveBezierWithSingleCubic( pCrvPart)) ; else { pCrvCubic.Set( ApproxArcCurveBezierWithSingleCubic( pCrvPart, ptCen, vtN)) ; if ( IsNull( pCrvCubic)) { // se fallisce allora riprovo usando più di una bezier bIsArc = false ; bOneIsEnough = false ; nParts = 2 ; pCrvPart->TrimStartEndAtParam( double( nParts - 1) / nParts, 1) ; pCrvCubic.Set( ApproxCurveBezierWithSingleCubic( pCrvPart)) ; } } CalcApproxError( pCrvPart, pCrvCubic, dErr) ; if ( dErr > dTol && ! bOneIsEnough) nParts += 2 ; else { pCC->AddCurve( Release( pCrvCubic)) ; if ( bOneIsEnough) break ; } ++ nCount ; } if ( nCount < 100) { // approssimo annche tutti gli altri tratti con una cubica for ( int i = nParts - 1 ; i > 0 ; --i) { PtrOwner pCrvPart( pCrvBezier->Clone()) ; pCrvPart->TrimStartEndAtParam( double( i - 1) / nParts, double(i) / nParts) ; PtrOwner pCrvCubic( ApproxCurveBezierWithSingleCubic( pCrvPart)) ; if ( ! pCC->AddCurve( Release( pCrvCubic), false)) return nullptr ; } return Release( pCC) ; } else return nullptr ; } //---------------------------------------------------------------------------- ICurveBezier* ApproxCurveBezierWithSingleCubic( const ICurve* pCrv) { // verifico sia una bezier const CurveBezier* pCrvBez = GetBasicCurveBezier( pCrv) ; if ( pCrvBez == nullptr) return nullptr ; Point3d ptStart, ptEnd ; pCrvBez->GetStartPoint( ptStart) ; pCrvBez->GetEndPoint( ptEnd) ; Vector3d vtStartDir, vtEndDir ; pCrvBez->GetStartDir( vtStartDir) ; pCrvBez->GetEndDir( vtEndDir) ; Vector3d vtStartD1 ; Vector3d vtEndD1 ; Point3d ptCalc ; pCrvBez->GetPointD1D2(0, ICurve::FROM_PLUS, ptCalc, &vtStartD1) ; pCrvBez->GetPointD1D2(1, ICurve::FROM_MINUS, ptCalc, &vtEndD1) ; PtrOwner pCrvCubic( CreateCurveBezier()) ; pCrvCubic->Init( 3, false) ; // approssimo per tentativi questo tratto di curva con una cubica // come vincolo tengo i punti e le direzioni di start e di end pCrvCubic->SetControlPoint( 0, ptStart) ; pCrvCubic->SetControlPoint( 3, ptEnd) ; // calcolo i punti intermedi per approssimare la curva Point3d ptInter1 = ptStart + 1./3 * vtStartD1 ; Point3d ptInter2 = ptEnd - 1./3 * vtEndD1 ; pCrvCubic->SetControlPoint( 1, ptInter1) ; pCrvCubic->SetControlPoint( 2, ptInter2) ; return Release( pCrvCubic) ; } //---------------------------------------------------------------------------- ICurveBezier* ApproxArcCurveBezierWithSingleCubic( const ICurve* pCrv, const Point3d& ptCen, const Vector3d& vtN) { // verifico sia una bezier const CurveBezier* pCrvBez = GetBasicCurveBezier( pCrv) ; if ( pCrvBez == nullptr) return nullptr ; // mi metto nel frame della curva ( lavoro nel piano XY) Point3d ptStart = pCrvBez->GetControlPoint( 0) ; Point3d ptEnd = pCrvBez->GetControlPoint( 2) ; Frame3d frCrv; frCrv.Set( ptStart, vtN) ; ptStart = ORIG ; ptEnd.ToLoc( frCrv) ; Point3d ptCenLoc = ptCen ; ptCenLoc.ToLoc( frCrv) ; // converto una curva di bezier che definisce un arco perfetto di circonferenza // dato il centro( della circonferenza), inizio e fine // N.B. : per archi con angolo al centro < 90 PtrOwner pCrvCubic( CreateCurveBezier()) ; int nDeg = pCrvBez->GetDegree() ; bool bRat = pCrvBez->IsRational() ; if ( nDeg != 2 || ! bRat) return nullptr ; if ( AreSamePointEpsilon( ptStart, ptEnd, EPS_SMALL * 50)) return nullptr ; nDeg = 3 ; bRat = false ; pCrvCubic->Init( nDeg, bRat) ; pCrvCubic->SetControlPoint( 0, ptStart) ; pCrvCubic->SetControlPoint( 3, ptEnd) ; // from the article of Aleksas Riškus "Approximation of a cubic bezier curve by circular arcs and vice versa" // with corrections from "Hans Muller's Flex Blog", page "More About Approximating Circular Arcs With a Cubic Bezier Path" double ax = ptStart.x - ptCenLoc.x ; double ay = ptStart.y - ptCenLoc.y ; double bx = ptEnd.x - ptCenLoc.x ; double by = ptEnd.y - ptCenLoc.y ; double q1 = ax * ax + ay * ay ; double q2 = q1 + ax * bx + ay * by ; double k2 = 4. / 3 * ( sqrt( 2 * q1 * q2) - q2) / ( ax * by - ay * bx) ; Point3d ptCtrl1, ptCtrl2 ; ptCtrl1.x = ptStart.x - k2 * ay ; ptCtrl1.y = ptStart.y + k2 * ax ; ptCtrl2.x = ptEnd.x + k2 * by ; ptCtrl2.y = ptEnd.y - k2 * bx ; pCrvCubic->SetControlPoint( 1, ptCtrl1) ; pCrvCubic->SetControlPoint( 2, ptCtrl2) ; // ritorno nel frame originale pCrvCubic->ToGlob( frCrv) ; return Release( pCrvCubic) ; } //---------------------------------------------------------------------------- static bool FindSpan( double dU, int nDeg, const DBLVECTOR& vKnots, int& nSpan) { if ( dU < 0) return false ; else if ( dU < EPS_ZERO) { nSpan = nDeg - 1 ; return true ; } // trovo a quale span appartiene il parametro dU int nKnots = int( vKnots.size()) ; if ( abs( dU - vKnots[nKnots-1]) < EPS_SMALL) { nSpan = nKnots - 1 - nDeg ; return true ; } int nLow = nDeg - 1 ; int nHigh = nKnots - 1 ; int nMid = ( nLow + nHigh) / 2 ; while ( dU < vKnots[nMid] || dU >= vKnots[nMid + 1] ) { if ( dU < vKnots[nMid]) nHigh = nMid ; else nLow = nMid ; nMid = ( nLow + nHigh) / 2 ; if ( nMid == nDeg - 1) break ; } nSpan = nMid ; return true ; } //---------------------------------------------------------------------------- static bool CalcBasisFunc( double dU, int nSpan, int nDeg, const DBLVECTOR& vKnots, DBLVECTOR& vBasis) { // mi aspetto che il vettore vBasis sia di lunghezza nDeg + 1 if ( vBasis.size() != nDeg + 1) return false ; vBasis[0] = 1 ; DBLVECTOR vLeft ; vLeft.resize( nDeg + 1) ; DBLVECTOR vRight ; vRight.resize( nDeg + 1) ; for ( int j = 1 ; j <= nDeg ; ++j) { vLeft[j] = dU - vKnots[nSpan + 1 -j] ; vRight[j] = vKnots[nSpan + j] - dU ; double dSaved = 0 ; for ( int r = 0 ; r < j ; ++r) { double dSum = vRight[r+1] + vLeft[j-r] ; double dTemp = 0 ; if ( dSum > EPS_ZERO) dTemp = vBasis[r] / dSum ; vBasis[r] = dSaved + vRight[r+1] * dTemp ; dSaved = vLeft[j-r] * dTemp ; } vBasis[j] = dSaved ; } return true ; } //---------------------------------------------------------------------------- ICurve* InterpolatePointSetWithBezierNoIntermedLines( const PNTVECTOR& vPnt, int nStart, int nEnd, int nDeg, const DBLVECTOR& vLen, double dLenTot) { PtrOwner pCrvInt ; int nPoints = nEnd - nStart + 1 ; if ( nPoints < 2) return nullptr ; else if ( nPoints == 2) { // se ho solo due punti restituisco una linea CurveLine cl ; cl.Set( vPnt[nStart], vPnt[nEnd]) ; pCrvInt.Set( LineToBezierCurve( &cl, nDeg, false)) ; if ( ! IsNull( pCrvInt) && pCrvInt->IsValid()) return Release( pCrvInt) ; else return nullptr ; } else if (nPoints == 3) { // se ho solo tre punti uso un altro algoritmo CurveByInterp cbi ; for ( int i = nStart ; i <= nEnd ; ++i) cbi.AddPoint( vPnt[i]) ; pCrvInt.Set( cbi.GetCurve( CurveByInterp::AKIMA_CORNER, CurveByInterp::CUBIC_BEZIERS)) ; if ( ! IsNull( pCrvInt) && pCrvInt->IsValid()) return Release( pCrvInt) ; else return nullptr ; } DBLVECTOR vPntParam ; vPntParam.resize( nPoints) ; vPntParam[0] = 0 ; vPntParam.back() = 1 ; for ( int i = 1 ; i < nPoints - 1 ; ++i) vPntParam[i] = vPntParam[i-1] + vLen[i-1] / dLenTot ; DBLVECTOR vKnots ; vKnots.resize( nPoints + nDeg - 1) ; for ( int i = 0 ; i < nDeg ; ++i) { vKnots[i] = 0 ; vKnots.end()[-i-1] = 1 ; } for ( int i = nDeg ; i < nPoints - 1 ; ++i) { double dKnot = 0 ; for ( int j = i + 1 ; j < i + nDeg + 1 ; ++j) dKnot += vPntParam[j - nDeg] ; dKnot /= nDeg ; vKnots[i] = dKnot ; } Eigen::MatrixXd mA( nPoints, nPoints) ; mA.fill( 0) ; for ( int i = 0 ; i < nPoints ; ++i) { if ( i == 0) mA.row(0).col(0) << 1 ; else if ( i == nPoints - 1) mA.row(i).col(nPoints - 1) << 1 ; else { int nSpan = 0 ; FindSpan( vPntParam[i], nDeg, vKnots, nSpan) ; DBLVECTOR vBasis ; vBasis.resize( nDeg + 1) ; CalcBasisFunc( vPntParam[i], nSpan, nDeg, vKnots, vBasis) ; for ( int j = nSpan - nDeg + 1 ; j <= nSpan + 1 ; ++j) mA.row(i).col(j) << vBasis[j - nSpan + nDeg - 1] ; } } int nDim = 3 ; Eigen::MatrixXd mb ; mb.resize( nPoints, nDim) ; for ( int i = 0 ; i < nPoints ; ++i) { mb.row(i).col(0) << vPnt[nStart + i].x ; mb.row(i).col(1) << vPnt[nStart + i].y ; mb.row(i).col(2) << vPnt[nStart + i].z ; } Eigen::MatrixXd mX = mA.fullPivLu().solve(mb) ; PNTVECTOR vPntCtrl ; for ( int i = 0 ; i < nPoints ; ++i) vPntCtrl.emplace_back( mX(i,0), mX(i,1), mX(i,2)) ; CNurbsData cNurbs ; cNurbs.nDeg = nDeg ; cNurbs.vCP = vPntCtrl ; cNurbs.vU = vKnots ; cNurbs.bRat = false ; pCrvInt.Set( NurbsToBezierCurve( cNurbs)) ; if ( ! IsNull(pCrvInt) && pCrvInt->IsValid()) return Release( pCrvInt) ; else return nullptr ; } //---------------------------------------------------------------------------- ICurve* InterpolatePointSetWithBezier( const PNTVECTOR& vPnt, double dLinTol, double dMaxLen) { PolyLine pl ; pl.FromPointVector( vPnt) ; PtrOwner pCrvOri( CreateCurveComposite()) ; pCrvOri->FromPolyLine( pl) ; PtrOwner pCrvInt( CreateCurveComposite()) ; // pag 364 del Piegl // scelgo il parametro associato ad ogni punto in modo che rispecchi la distanza tra i punti double dErr = INFINITO ; int nItCount = 0 ; while ( dErr > dLinTol && nItCount < 10) { pCrvInt->Clear() ; int nPoints = int( vPnt.size()) ; int nDeg = 3 ; if ( nPoints < 2) return nullptr ; else if ( nPoints == 2) { // se ho solo due punti restituisco una linea CurveLine cl ; cl.Set( vPnt[0], vPnt[1]) ; pCrvInt->AddCurve( LineToBezierCurve( &cl, nDeg, false)) ; if ( ! IsNull( pCrvInt) && pCrvInt->IsValid()) return Release( pCrvInt) ; else return nullptr ; } else if ( nPoints == 3) { // se ho solo tre punti uso un altro algoritmo CurveByInterp cbi ; for ( int i = 0 ; i < int( vPnt.size()) ; ++i) cbi.AddPoint( vPnt[i]) ; pCrvInt->AddCurve( cbi.GetCurve( CurveByInterp::AKIMA_CORNER, CurveByInterp::CUBIC_BEZIERS)) ; if ( ! IsNull( pCrvInt) && pCrvInt->IsValid()) return Release( pCrvInt) ; else return nullptr ; } // scorro il vettore di punti passato in input e interpolo separatamente gruppi di punti se incontro // tratti più lunghi della distanza dMaxLen int nStart = 0, nEnd = 0 ; bool bOk = true ; bool bFoundLine = false ; while ( nStart < nPoints - 1 && bOk) { bFoundLine = false ; nEnd = 0 ; DBLVECTOR vLen ; double dLenTot = 0 ; for ( int i = nStart ; i < nPoints - 1 ; ++i) { double dLen = Dist( vPnt[i], vPnt[i+1]) ; if ( dLen < dMaxLen) { vLen.push_back( dLen) ; dLenTot += dLen ; } else { nEnd = i ; bFoundLine = true ; break ; } } if ( vLen.size() != 0) { if ( nEnd == 0) nEnd = nPoints - 1 ; pCrvInt->AddCurve( InterpolatePointSetWithBezierNoIntermedLines( vPnt, nStart, nEnd, nDeg, vLen, dLenTot)) ; } if ( bFoundLine) { CurveLine cl ; cl.Set( vPnt[nEnd], vPnt[nEnd+1]) ; pCrvInt->AddCurve( LineToBezierCurve( &cl, nDeg, false)) ; ++ nEnd ; } bOk = pCrvInt->IsValid() ; nStart = nEnd ; } CalcApproxError( pCrvOri, pCrvInt, dErr) ; if ( dErr > dLinTol && dMaxLen > 200 * EPS_SMALL) dMaxLen /= 2 ; ++ nItCount ; } if ( ! IsNull( pCrvInt) && pCrvInt->IsValid() && dErr < dLinTol) return Release( pCrvInt) ; else return nullptr ; } //---------------------------------------------------------------------------- ICurve* ApproxCurveWithBezier( const ICurve* pCrv , double dTol, int nType) { PolyLine plApprox ; double dAngTolFine = 2 ; pCrv->ApproxWithLines( dTol, dAngTolFine, ICurve::APL_STD, plApprox) ; PNTVECTOR vPnt ; Point3d pt ; plApprox.GetFirstPoint( pt) ; do { vPnt.push_back( pt) ; } while ( plApprox.GetNextPoint( pt)) ; // campiono punti lungo la curva e poi li interpolo PtrOwner pCC( InterpolatePointSetWithBezier( vPnt, dTol, 100)) ; if ( ! IsNull( pCC) && pCC->IsValid()) return Release( pCC) ; else return nullptr ; } //---------------------------------------------------------------------------- ICurve* ApproxPointSetWithBezier( const ICurve* pCrv, double dTol) { // campiono punti lungo la curva e poi li interpolo PtrOwner pCC( CreateBasicCurveComposite()) ; return Release( pCC) ; } //---------------------------------------------------------------------------- bool CalcApproxError( const ICurve* pCrvOri, const ICurve* pCrvNew, double& dErr, int nPoints) { // controllo l'errore effettivo campionando più finemente double dLenOri = 0 ; pCrvOri->GetLength( dLenOri) ; double dLenNew = 0 ; pCrvNew->GetLength( dLenNew) ; dErr = 0 ; for ( int i = 1 ; i < nPoints ; ++i) { Point3d ptOri, ptNew ; double dParOri, dParNew ; pCrvOri->GetParamAtLength( double (i) / nPoints * dLenOri, dParOri) ; pCrvOri->GetPointD1D2(dParOri, ICurve::FROM_MINUS, ptOri) ; pCrvNew->GetParamAtLength( double (i) / nPoints * dLenNew, dParNew) ; pCrvNew->GetPointD1D2(dParNew, ICurve::FROM_MINUS, ptNew) ; ; double dErrLoc = Dist( ptOri, ptNew) ; if ( dErrLoc > dErr) dErr = dErrLoc ; } return true ; } //---------------------------------------------------------------------------- ICurve* CurveToNoArcsCurve( const ICurve* pCrv) { // verifico validità curva if ( pCrv == nullptr) return nullptr ; // se arco, devo trasformarlo in curva di Bezier (semplice o composta) if ( pCrv->GetType() == CRV_ARC) { return ArcToBezierCurve( GetCurveArc( pCrv)) ; } // se curva composita, devo trasformarla in composita senza archi else if ( pCrv->GetType() == CRV_COMPO) { PtrOwner pCopyCC( GetBasicCurveComposite( pCrv->Clone())) ; if ( IsNull( pCopyCC) || ! pCopyCC->ArcsToBezierCurves()) return nullptr ; return Release( pCopyCC) ; } // altrimenti devo solo copiarla else { return pCrv->Clone() ; } } //---------------------------------------------------------------------------- ICurve* CurveToArcsPerpExtrCurve( const ICurve* pCrv, double dLinTol, double dAngTolDeg) { // verifico validità curva if ( pCrv == nullptr) return nullptr ; // se arco in piano non perpendicolare ad estrusione, curva composita o curva di Bezier trasformo if ( ( pCrv->GetType() == CRV_ARC && ! GetBasicCurveArc(pCrv)->IsInPlanePerpExtr()) || pCrv->GetType() == CRV_COMPO || pCrv->GetType() == CRV_BEZIER) { // creo un poliarco PolyArc PA ; if ( ! pCrv->ApproxWithArcs( dLinTol, dAngTolDeg, PA)) return nullptr ; // creo una nuova curva composita a partire dal poliarco PtrOwner pCopyCC( CreateBasicCurveComposite()) ; if ( IsNull( pCopyCC)) return nullptr ; if ( ! pCopyCC->FromPolyArc( PA)) return nullptr ; // copio estrusione e spessore CopyExtrusion( pCrv, pCopyCC) ; CopyThickness( pCrv, pCopyCC) ; return Release( pCopyCC) ; } // altrimenti devo solo copiarla else { return pCrv->Clone() ; } } //---------------------------------------------------------------------------- bool NurbsCurveCanonicalize( CNurbsData& cnData) { // se con nodi extra if ( cnData.bExtraKnotes) { int nKnotesNbr = int( cnData.vU.size()) ; if ( nKnotesNbr < 4) return false ; cnData.bExtraKnotes = false ; for ( int i = 0 ; i < nKnotesNbr - 2 ; ++ i) cnData.vU[i] = cnData.vU[i+1] ; cnData.vU.resize( nKnotesNbr - 2) ; } // se periodica if ( cnData.bPeriodic || ! cnData.bClamped) { bool bAlreadyChecked = false ; // se la curva è peridica verifco che effettivamente ci sia un numero di punti ripetituti uguale al grado della curva // wrap della curva su se stessa if ( cnData.bPeriodic && (int(cnData.vU.size()) > int(cnData.vCP.size()) + cnData.nDeg - 1)) { bool bRepeated = true ; for ( int i = 0 ; i < cnData.nDeg ; ++i) { if ( ! AreSamePointApprox( cnData.vCP[i], cnData.vCP.end()[-cnData.nDeg + i]) ) { bRepeated = false ; break ; } } bool bFirstAddedAtEnd = false ; if ( ! bRepeated || (bRepeated && AreSamePointApprox( cnData.vCP[0],cnData.vCP[cnData.nDeg]))){ // salvo il vettore dei nodi in caso mi accorga di avere tra le mani una curva unclamped DBLVECTOR vU = cnData.vU ; // se effettivamente ho dei nodi in più da togliere allora li tolgo ed eventualmente aggiungo punti di controllo if ( int(cnData.vU.size()) > int(cnData.vCP.size()) + cnData.nDeg - 1 ) { // se il primo e l'ultimo punto non coincidono allora aggiungo il primo punto in fondo al vettore dei punti di controllo if ( ! AreSamePointApprox( cnData.vCP[0], cnData.vCP.back())) { bFirstAddedAtEnd = true ; cnData.vCP.push_back( cnData.vCP[0]) ; if ( cnData.bRat) cnData.vW.push_back( cnData.vW[0]) ; } // devo poi anche togliere i nodi di troppo // presuppongo che la convenzione sia che i nodi di troppo siano alla fine del vettore dei nodi cnData.vU = DBLVECTOR( cnData.vU.begin(), cnData.vU.end() - cnData.nDeg) ; // controllo eventualmente anche i nodi extra // se ne ho due in più ne tolgo uno in cima e uno in fondo if ( cnData.vU.size() == int( cnData.vCP.size()) + cnData.nDeg + 1 ) { // significa che ci sono due nodi extra, uno all'inizio e uno alla fine, da togliere cnData.vU = vector( cnData.vU.begin() + 1, cnData.vU.end() - 1) ; } // se ne ho solo uno in più lo tolgo in cima else if ( cnData.vU.size() == int( cnData.vCP.size()) + cnData.nDeg) { cnData.vU = vector( cnData.vU.begin() + 1, cnData.vU.end()) ; } } // controllo se il vettore dei nodi ha la giusta molteplicità all'inizio e alla fine, sennò ha comunque bisogno di essere resa non periodica double dU0 = cnData.vU[0] ; double dULast = cnData.vU.back() ; bool bSame = true ; for ( int i = 1 ; i < cnData.nDeg ; ++i ) { bSame = bSame && abs(cnData.vU[i] - dU0) < EPS_SMALL ; bSame = bSame && abs(cnData.vU.end()[-( i+ 1)] - dULast) < EPS_SMALL ; } if ( bSame) { cnData.bPeriodic = false ; return true ; } else { // aggiungo i punti ripetuti ( controllando se il primo l'ho già aggiunto o c'è già) bFirstAddedAtEnd = bFirstAddedAtEnd || AreSamePointApprox( cnData.vCP[0], cnData.vCP.back()) ; for ( int i = bFirstAddedAtEnd ? 1 : 0 ; i < cnData.nDeg ; ++i ) { cnData.vCP.push_back( cnData.vCP[i]) ; if ( cnData.bRat) cnData.vW.push_back( cnData.vW[i]) ; } // recupero il vettore dei nodi cnData.vU = vU ; // verifico se ho nodi extra if ( cnData.vU.size() == int( cnData.vCP.size()) + cnData.nDeg + 1 ) { // significa che ci sono due nodi extra: // se la curva ha grado maggiore di 1 e i primi due nodi sono uguali allora tolgo quelli if ( cnData.nDeg > 1 && abs(cnData.vU[1] - cnData.vU[0]) < EPS_SMALL) { cnData.vU = vector( cnData.vU.begin() + 2, cnData.vU.end()) ; } // sennò ne tolgo uno all'inizio e uno alla fine else cnData.vU = vector( cnData.vU.begin() + 1, cnData.vU.end() - 1) ; } // se ne ho solo uno in più lo tolgo in cima else if ( cnData.vU.size() == int( cnData.vCP.size()) + cnData.nDeg) cnData.vU = vector( cnData.vU.begin() + 1, cnData.vU.end()) ; } bAlreadyChecked = true ; } } if ( ! bAlreadyChecked) { // se non ho già controllato guardo se ho già la giusta molteplicità all'inizio e alla fine del vettore dei nodi double dU0 = cnData.vU[0] ; double dULast = cnData.vU.back() ; bool bSame = true ; for ( int i = 1 ; i < cnData.nDeg ; ++i ) { bSame = bSame && abs(cnData.vU[i] - dU0) < EPS_SMALL ; bSame = bSame && abs(cnData.vU.end()[-( i+ 1)] - dULast) < EPS_SMALL ; } if ( bSame) { cnData.bPeriodic = false ; return true ; } } // qui aggiungo un controllo se la curva è collassata in un punto ( ho un polo), lascio stare bool bCollapsed = true ; Point3d ptFirst = cnData.vCP.front() ; for ( int i = 1 ; i < int( cnData.vCP.size()) ; ++i) { if ( ! AreSamePointApprox( ptFirst, cnData.vCP[i])) { bCollapsed = false ; break ; } } if ( bCollapsed) return false ; // va trasformata in non-periodica (clamped) // bisogna aumentare la molteplicità dei nodi u_p-1 e u_(m-p+1) fino ad arrivare al grado della nurbs // e poi scartare nodi e punti fuori dalla regione clamped ( al di fuori della regione u_p-1 -> u_(m-p+1)) // l'agoritmo per l'inserimento dei nodi l' A5.1 del libro delle Nurbs ( Piegl e Tiller), con qualche modifica // agli indici perché uso u_p-1 e u_(m-p+1), anziché u_p e u_m-p // comincio ad aumentare la molteplictià del nodo u_m-p+1 int nCP = int( cnData.vCP.size()) ; int nU = nCP + cnData.nDeg - 1 ; int nDeg = cnData.nDeg ; PNTVECTOR vBC ; vBC.resize( nDeg + 1) ; DBLVECTOR vBW ; vBW.resize( nDeg + 1) ; // trovo il nodo di cui aumentare la molteplicità e ne calcolo la molteplicità int b = nU - nDeg - 1 + 1 ; int i = b ; //int c = b ; while ( b > 0 && abs( cnData.vU[b] - cnData.vU[b - 1]) < EPS_ZERO) -- b ; int mult = min( i - b + 1, nDeg) ; // mi aspetto che sia 1, ma comunque sarà < nDeg //while ( i > 0 && abs( cnData.vU[i] - cnData.vU[i - 1]) < EPS_ZERO) // -- i ; //int mult = min( b - i + 1, nDeg) ; // mi aspetto che sia 1, ma comunque sarà < nDeg //// devo controllare anche i nodi successivi! //while ( c < nU - 1 && abs( cnData.vU[c + 1] - cnData.vU[c]) < EPS_ZERO) // ++ c ; //int mult = min( c - i + 1, nDeg) ; // mi aspetto che sia 1, ma comunque sarà < nDeg // recupero i punti da modificare if ( ! cnData.bRat) { for ( int i = 0 ; i <= nDeg - mult ; ++ i) vBC[i] = cnData.vCP[b - nDeg + 1 + i] ; } else { for ( int i = 0 ; i <= nDeg - mult ; ++ i) { vBC[i] = cnData.vCP[b - nDeg + 1 + i] * cnData.vW[b - nDeg + 1 + i] ; vBW[i] = cnData.vW[b - nDeg + 1 + i] ; } } // salvo i punti inalterati int r = nDeg - mult ; // numero di volte che dovrò inserire il nodo cnData.vCP.resize( nCP + r) ; for ( int p = nCP - 1 ; p > b - mult ; --p) { cnData.vCP[r + p] = cnData.vCP[p] ; } if ( cnData.bRat ) { cnData.vW.resize( nCP + r) ; for ( int p = nCP - 1 ; p > b - mult ; --p) { cnData.vW[r + p] = cnData.vW[p] ; } } // procedo all'inserimento int L = 0 ; double alpha ; if ( mult < nDeg) { // inserisco il nodo r volte for ( int j = 1 ; j <= r ; ++ j) { L = b - nDeg + j ; for ( int i = 0; i <= r - j ; ++i) { alpha = (cnData.vU[b] - cnData.vU[L + i])/ ( cnData.vU[i + b + 1] - cnData.vU[L + i]) ; vBC[i] = alpha * vBC[i +1 ] + ( 1 - alpha) * vBC[i] ; if ( cnData.bRat) { vBW[i] = alpha * vBW[i + 1] + ( 1 - alpha) * vBW[i] ; } } cnData.vCP[L + 1] = vBC[0] ; cnData.vCP[b + nDeg - j - mult] = vBC[r - j] ; if ( cnData.bRat ) { cnData.vW[L + 1] = vBW[0] ; cnData.vW[b + nDeg - j - mult] = vBW[r-j] ; } } } // allungo il vettore dei nodi e sposto gli ultimi nodi cnData.vU.resize( nU + r) ; for ( int p = nU - 1 ; p > b ; --p) cnData.vU[p + r] = cnData.vU[p] ; // aggiungo i nodi nuovi for ( int p = 0 ; p < r ; ++p) cnData.vU[b + 1 + p] = cnData.vU[b] ; nU = nU + r ; nCP = nCP + r ; // aumento la molteplicità del punto u_p-1 b = nDeg - 1 ; i = b ; //c = b ; while ( b > 0 && abs( cnData.vU[b] - cnData.vU[b - 1]) < EPS_ZERO) -- b ; mult = min( i - b + 1, nDeg) ; // mi aspetto che sia 1, ma comunque sarà < cnData.nDeg //while ( i > 0 && abs( cnData.vU[i] - cnData.vU[i - 1]) < EPS_ZERO) // -- i ; //mult = min( b - i + 1, nDeg) ; // mi aspetto che sia 1, ma comunque sarà < cnData.nDeg //// devo controllare anche i nodi successivi! //while ( c < nU -1 && abs(cnData.vU[c + 1] - cnData.vU[c]) < EPS_ZERO ) // ++ c ; //mult = min( c - i + 1, nDeg) ; // mi aspetto che sia 1, ma comunque sarà < cnData.nDeg // recupero i punti da modificare if ( ! cnData.bRat) { for ( int i = 0 ; i <= nDeg - mult ; ++ i) vBC[i] = cnData.vCP[i] ; } else { for ( int i = 0 ; i <= nDeg - mult ; ++ i) { vBC[i] = cnData.vCP[i] * cnData.vW[i] ; vBW[i] = cnData.vW[i] ; } } r = nDeg - mult ; // salvo i punti inalterati cnData.vCP.resize( nCP + r) ; for ( int p = nCP - 1 ; p > b - mult ; --p) { cnData.vCP[r + p] = cnData.vCP[p] ; } if ( cnData.bRat ) { cnData.vW.resize( nCP + r) ; for ( int p = nCP - 1 ; p > b - mult ; --p) { cnData.vW[r + p] = cnData.vW[p] ; } } // procedo all'inserimento L = 0 ; if ( mult < nDeg) { // inserisco il nodo r volte for ( int j = 1 ; j <= r ; ++ j) { L = b - nDeg + j ; for ( int i = 0; i <= r - j ; ++i) { alpha = (cnData.vU[b] - cnData.vU[L + i])/ ( cnData.vU[i + b + 1] - cnData.vU[L + i]) ; vBC[i] = alpha * vBC[i + 1] + ( 1 - alpha) * vBC[i] ; if ( cnData.bRat) { vBW[i] = alpha * vBW[i + 1] + ( 1 - alpha) * vBW[i] ; } } cnData.vCP[L + 1] = vBC[0] ; cnData.vCP[b + nDeg - j - mult] = vBC[r - j] ; if ( cnData.bRat ) { cnData.vW[L + 1] = vBW[0] ; cnData.vW[b + nDeg - j - mult] = vBW[r - j] ; } } } // allungo il vettore dei nodi e sposto gli ultimi nodi cnData.vU.resize(nU + r) ; for ( int p = nU - 1 ; p > b ; --p) cnData.vU[p+r] = cnData.vU[p] ; // aggiungo i nodi nuovi for ( int p = 0 ; p < r ; ++p) cnData.vU[b + 1 + p] = cnData.vU[b] ; nU = nU + r ; nCP = nCP + r ; // rendo la curva chiusa e non periodica eliminando i primi e gli ultimi nDeg punti e nodi cnData.bPeriodic = false ; nCP = nCP - 2 * ( nDeg - 1); nU = nU - 2 * ( nDeg - 1); PNTVECTOR vCP_clamped ; vCP_clamped.resize( nCP) ; DBLVECTOR vW_clamped ; vW_clamped.resize( nCP) ; DBLVECTOR vU_clamped ; vU_clamped.resize( nU) ; for ( int i = 0 ; i < nCP ; ++i) { if ( ! cnData.bRat) vCP_clamped[i] = cnData.vCP[i + nDeg - 1] ; else { vW_clamped[i] = cnData.vW[i + nDeg - 1] ; vCP_clamped[i] = cnData.vCP[i + nDeg - 1] / cnData.vW[i + nDeg - 1] ; } } cnData.vCP = vCP_clamped ; cnData.vW = vW_clamped ; for ( int i = 0 ; i < nU ; ++i) { vU_clamped[i] = cnData.vU[i + nDeg - 1] ; } cnData.vU = vU_clamped ; } return true ; } //---------------------------------------------------------------------------- ICurve* NurbsToBezierCurve( const CNurbsData& cnData) { // la curva Nurbs deve essere in forma canonica if ( cnData.bPeriodic || cnData.bExtraKnotes) return nullptr ; // numero dei nodi int nU = int( cnData.vCP.size()) + cnData.nDeg - 1 ; // controllo relazione nodi - punti di controllo if ( nU != int( cnData.vU.size())) return nullptr ; // numero degli intervalli int nInt = nU - 2 * cnData.nDeg + 1 ; // sistemo le condizioni agli estremi sui nodi (i primi nDeg nodi e gli ultimi nDeg nodi devono essere uguali tra loro) for ( int i = 0 ; i < cnData.nDeg - 1 ; ++ i) { if ( abs( cnData.vU[i] - cnData.vU[cnData.nDeg - 1]) >= EPS_ZERO) const_cast( cnData.vU[i]) = cnData.vU[cnData.nDeg - 1] ; } for ( int i = 0 ; i < cnData.nDeg - 1 ; ++ i) { if ( abs( cnData.vU[nU - 1 - i] - cnData.vU[nU - cnData.nDeg]) >= EPS_ZERO) const_cast( cnData.vU[nU - 1 - i]) = cnData.vU[nU - cnData.nDeg] ; } // se 1 solo intervallo, la Nurbs è già una curva di Bezier if ( nInt == 1) { // creo la curva di Bezier PtrOwner pCrvBez( CreateBasicCurveBezier()) ; if ( IsNull( pCrvBez)) return nullptr ; // la inizializzo if ( ! pCrvBez->Init( cnData.nDeg, cnData.bRat)) return nullptr ; for ( int i = 0 ; i <= cnData.nDeg ; ++ i) { if ( ! cnData.bRat) { if ( ! pCrvBez->SetControlPoint( i, cnData.vCP[i])) return nullptr ; } else { if ( ! pCrvBez->SetControlPoint( i, cnData.vCP[i], cnData.vW[i])) return nullptr ; } } // se non è una curva ma un punto, la invalido if ( pCrvBez->IsAPoint()) pCrvBez->Init( cnData.nDeg, cnData.bRat) ; // restituisco la curva return Release( pCrvBez) ; } // altrimenti è equivalente ad una curva composita, la creo PtrOwner pCrvCompo( CreateBasicCurveComposite()) ; if ( IsNull( pCrvCompo)) return nullptr ; // vettore dei punti di controllo della curva di Bezier PNTVECTOR vBC ; vBC.resize( cnData.nDeg + 1) ; DBLVECTOR vBW ; vBW.resize( cnData.nDeg + 1) ; if ( ! cnData.bRat) { for ( int i = 0 ; i <= cnData.nDeg ; ++ i) vBC[i] = cnData.vCP[i] ; } else { for ( int i = 0 ; i <= cnData.nDeg ; ++ i) { vBC[i] = cnData.vCP[i] * cnData.vW[i] ; vBW[i] = cnData.vW[i] ; } } // primi coefficienti della successiva PNTVECTOR vNextBC ; vNextBC.resize( cnData.nDeg - 1) ; DBLVECTOR vNextBW ; vNextBW.resize( cnData.nDeg - 1) ; // ... DBLVECTOR vAlfa ; vAlfa.resize( cnData.nDeg - 1) ; int a = cnData.nDeg - 1 ; int b = cnData.nDeg ; bool bPrevRejected = false ; //algoritmo A5.6 di Piegl e Tiller // ciclo while ( b < nU - 1) { int i = b ; while ( b < nU - 1 && abs( cnData.vU[b+1] - cnData.vU[b]) < EPS_ZERO) ++ b ; int mult = min( b - i + 1, cnData.nDeg) ; if ( mult < cnData.nDeg) { // numeratore di alfa double numer = cnData.vU[b] - cnData.vU[a] ; // calcola e salva gli alfa for ( int j = cnData.nDeg ; j > mult ; -- j) vAlfa[j-mult-1] = numer / ( cnData.vU[a+j] - cnData.vU[a]) ; // inserisco il nodo r volte int r = cnData.nDeg - mult ; for ( int j = 1 ; j <= r ; ++ j) { int save = r - j ; int s = mult + j ; for ( int k = cnData.nDeg ; k >= s ; -- k) vBC[k] = vAlfa[k-s] * vBC[k] + ( 1 - vAlfa[k-s]) * vBC[k-1] ; if ( cnData.bRat) { for ( int k = cnData.nDeg ; k >= s ; -- k) vBW[k] = vAlfa[k-s] * vBW[k] + ( 1 - vAlfa[k-s]) * vBW[k-1] ; } if ( b < nU - 1) { vNextBC[save] = vBC[cnData.nDeg] ; if ( cnData.bRat) vNextBW[save] = vBW[cnData.nDeg] ; } } } // costruisco la curva di Bezier e la inserisco nella curva composita PtrOwner pCrvBez( CreateBasicCurveBezier()) ; if ( IsNull( pCrvBez)) return nullptr ; // se precedente saltata if ( bPrevRejected) { // prendo l'ultimo punto della curva composita per garantire la continuità Point3d ptEnd ; if ( pCrvCompo->GetEndPoint( ptEnd)) vBC[0] = ptEnd ; } // la inizializzo if ( ! pCrvBez->Init( cnData.nDeg, cnData.bRat)) return nullptr ; if ( ! cnData.bRat) { for ( int i = 0 ; i <= cnData.nDeg ; ++ i) { if ( ! pCrvBez->SetControlPoint( i, vBC[i])) return nullptr ; } } else { for ( int i = 0 ; i <= cnData.nDeg ; ++ i) { if ( ! pCrvBez->SetControlPoint( i, vBC[i] / vBW[i], vBW[i])) return nullptr ; } } // se è una vera curva, la aggiungo alla curva composita if ( ! pCrvBez->IsAPoint()) { if ( ! pCrvCompo->AddCurve( Release( pCrvBez))) return nullptr ; bPrevRejected = false ; } // altrimenti è un punto, la cancello else { pCrvBez.Reset() ; bPrevRejected = true ; } // inizializzazioni per la prossima curva di Bezier if ( b < nU - 1) { if ( ! cnData.bRat) { for ( int i = 0 ; i < cnData.nDeg - 1 ; ++ i) vBC[i] = vNextBC[i] ; for ( int i = cnData.nDeg - mult ; i <= cnData.nDeg ; ++ i) vBC[i] = cnData.vCP[b-cnData.nDeg+i+1] ; } else { for ( int i = 0 ; i < cnData.nDeg - 1 ; ++ i) { vBC[i] = vNextBC[i] ; vBW[i] = vNextBW[i] ; } for ( int i = cnData.nDeg - mult ; i <= cnData.nDeg ; ++ i) { vBC[i] = cnData.vCP[b-cnData.nDeg+i+1] * cnData.vW[b-cnData.nDeg+i+1] ; vBW[i] = cnData.vW[b-cnData.nDeg+i+1] ; } } a = b ; ++ b ; } } // se la curva ha grado 1, manca da aggiungere l'ultimo tratto if ( cnData.nDeg == 1 ) { // costruisco la curva di Bezier e la inserisco nella curva composita PtrOwner pCrvBez( CreateBasicCurveBezier()) ; if ( ! pCrvBez->Init( cnData.nDeg, cnData.bRat)) return nullptr ; if ( ! cnData.bRat) { for ( int i = 0 ; i <= cnData.nDeg ; ++ i) { if ( ! pCrvBez->SetControlPoint( i, vBC[i])) return nullptr ; } } else { for ( int i = 0 ; i <= cnData.nDeg ; ++ i) { if ( ! pCrvBez->SetControlPoint( i, vBC[i] / vBW[i], vBW[i])) return nullptr ; } } if ( ! pCrvBez->IsAPoint()) { if ( ! pCrvCompo->AddCurve( Release( pCrvBez))) return nullptr ; } } // restituisco la curva composita return Release( pCrvCompo) ; } //---------------------------------------------------------------------------- ICurve* FlattenCurve( const ICurve& crCrv, double dToler, double dAngToler, int nFlag) { // Determino il piano medio della curva e verifico se scostamento accettabile Plane3d plMid ; if ( ! crCrv.IsFlat( plMid, ( nFlag == FLTCRV_USE_EXTR), dToler)) return nullptr ; // Recupero estrusione rispetto alla normale al piano medio Vector3d vtExtr ; crCrv.GetExtrusion( vtExtr) ; if ( nFlag == FLTCRV_USE_EXTR) plMid.Set( plMid.GetPoint(), vtExtr) ; // Determino il suo centro geometrico Point3d ptCen ; if ( ! crCrv.GetCentroid( ptCen)) return nullptr ; // Verifico se curva già piatta PolyLine PL ; if ( ! crCrv.ApproxWithLines( LIN_TOL_FINE, ANG_TOL_STD_DEG, ICurve::APL_SPECIAL_INT, PL)) return nullptr ; bool bFlat = true ; Plane3d plFlat ; plFlat.Set( ptCen, plMid.GetVersN()) ; Point3d ptP ; bool bPoint = PL.GetFirstPoint( ptP) ; while ( bPoint && bFlat) { if ( DistPointPlane( ptP, plFlat) > 0.1 * EPS_SMALL) bFlat = false ; bPoint = PL.GetNextPoint( ptP) ; } // Se curva già piatta, la copio ed esco if ( bFlat) { PtrOwner pCrv( crCrv.Clone()) ; if ( IsNull( pCrv)) return nullptr ; // assegno estrusione if ( nFlag != FLTCRV_SET_EXTR) pCrv->SetExtrusion( vtExtr) ; else { if ( vtExtr * plMid.GetVersN() >= 0) pCrv->SetExtrusion( plMid.GetVersN()) ; else pCrv->SetExtrusion( - plMid.GetVersN()) ; } // restituisco la copia della curva return Release( pCrv) ; } // altrimenti la appiattisco else { // mi assicuro che la curva non contenga archi PtrOwner pCrv( CurveToNoArcsCurve( &crCrv)) ; if ( IsNull( pCrv)) return nullptr ; // calcolo un riferimento con piano XY coincidente con il piano di proiezione Frame3d frRef ; if ( ! frRef.Set( ptCen, plMid.GetVersN())) return nullptr ; // eseguo scalatura con fattori X e Y unitari e fattore Z nullo if ( ! pCrv->Scale( frRef, 1, 1, 0)) return nullptr ; // assegno estrusione if ( nFlag != FLTCRV_SET_EXTR) pCrv->SetExtrusion( vtExtr) ; else { if ( vtExtr * plMid.GetVersN() >= 0) pCrv->SetExtrusion( plMid.GetVersN()) ; else pCrv->SetExtrusion( - plMid.GetVersN()) ; } // restituisco la nuova curva return Release( pCrv) ; } } //---------------------------------------------------------------------------- ICurve* ProjectCurveOnPlane( const ICurve& crCrv, const Plane3d& plPlane) { // determino se curva piana e suo eventuale piano Plane3d plCrv ; if ( crCrv.IsFlat( plCrv, false, EPS_SMALL / 2)) { // se il piano della curva è parallelo a quello di proiezione if ( AreSameOrOppositeVectorExact( plCrv.GetVersN(), plPlane.GetVersN())) { // copio la curva PtrOwner pCrv( crCrv.Clone()) ; if ( IsNull( pCrv)) return nullptr ; // se non coincidenti, basta eseguire una traslazione Point3d ptOC = ORIG + plCrv.GetDist() * plCrv.GetVersN() ; Point3d ptOP = ORIG + plPlane.GetDist() * plPlane.GetVersN() ; if ( ! AreSamePointApprox( ptOC, ptOP)) pCrv->Translate( ptOP - ptOC) ; // restituisco la nuova curva return Release( pCrv) ; } } // mi assicuro che la curva non contenga archi PtrOwner pCrv( CurveToNoArcsCurve( &crCrv)) ; if ( IsNull( pCrv)) return nullptr ; // calcolo un riferimento con piano XY coincidente con il piano di proiezione Frame3d frRef ; if ( ! frRef.Set( ORIG + plPlane.GetDist() * plPlane.GetVersN(), plPlane.GetVersN())) return nullptr ; // eseguo scalatura con fattori X e Y unitari e fattore Z nullo if ( ! pCrv->Scale( frRef, 1, 1, 0)) return nullptr ; // restituisco la nuova curva return Release( pCrv) ; } //---------------------------------------------------------------------------- bool AdjustCurveSlope( ICurveComposite* pCrv, double dNini, double dNfin) { // verifico curva if ( pCrv == nullptr) return false ; // determino versore estrusione Vector3d vtN ; if ( ! pCrv->GetExtrusion( vtN) || vtN.IsSmall()) vtN = Z_AX ; // assegno la corretta pendenza int i = 0 ; double dCurrLen = 0 ; double dLen ; pCrv->GetLength( dLen) ; const ICurve* pSCrv = pCrv->GetFirstCurve() ; while ( pSCrv != nullptr) { double dCrvLen ; pSCrv->GetLength( dCrvLen) ; double dCoeff = dCurrLen / dLen ; double dCurrN = dNini * ( 1.0 - dCoeff) + dNfin * dCoeff ; Point3d ptJoin ; pSCrv->GetStartPoint( ptJoin) ; pCrv->ModifyJoint( i, ptJoin + vtN * ( dCurrN - ( ptJoin - ORIG) * vtN)) ; // passo al successivo dCurrLen += dCrvLen ; pSCrv = pCrv->GetNextCurve() ; ++ i ; } Point3d ptFin ; pCrv->GetEndPoint( ptFin) ; ptFin += vtN * ( dNfin - ( ptFin - ORIG) * vtN) ; pCrv->ModifyEnd( ptFin) ; return true ; } //---------------------------------------------------------------------------- Voronoi* GetCurveVoronoi( const ICurve& crvC) { switch ( crvC.GetType()) { case CRV_LINE : return GetBasicCurveLine( &crvC)->GetVoronoiObject() ; case CRV_ARC : return GetBasicCurveArc( &crvC)->GetVoronoiObject() ; case CRV_BEZIER : return GetBasicCurveBezier( &crvC)->GetVoronoiObject() ; case CRV_COMPO : return GetBasicCurveComposite( &crvC)->GetVoronoiObject() ; default : return nullptr ; } return nullptr ; } //---------------------------------------------------------------------------- bool CalcCurveVoronoiDiagram( const ICurve& crvC, ICURVEPOVECTOR& vCrvs, int nBound) { Voronoi* pVoronoiObj = GetCurveVoronoi( crvC) ; if ( pVoronoiObj == nullptr) return false ; return pVoronoiObj->CalcVoronoiDiagram( vCrvs, nBound) ; } //---------------------------------------------------------------------------- bool CalcCurvesVoronoiDiagram( const CICURVEPVECTOR& vCrvC, ICURVEPOVECTOR& vCrvs, int nBound) { // creo oggetto Voronoi con le curve passate PtrOwner pVoronoiObj( new( std::nothrow) Voronoi()) ; if ( pVoronoiObj == nullptr) return false ; for ( int i = 0 ; i < int( vCrvC.size()) ; i ++) { if ( ! pVoronoiObj->AddCurve( vCrvC[i])) return false ; } return pVoronoiObj->CalcVoronoiDiagram( vCrvs, nBound) ; } //---------------------------------------------------------------------------- bool CalcCurveMedialAxis( const ICurve& crvC, ICURVEPOVECTOR& vCrvs, int nSide) { Voronoi* pVoronoiObj = GetCurveVoronoi( crvC) ; if ( pVoronoiObj == nullptr) return false ; return pVoronoiObj->CalcMedialAxis( vCrvs, nSide) ; } //---------------------------------------------------------------------------- bool CalcCurvesMedialAxis( const CICURVEPVECTOR& vCrvC, ICURVEPOVECTOR& vCrvs, int nSide) { // creo oggetto Voronoi con le curve passate PtrOwner pVoronoiObj( new( std::nothrow) Voronoi()) ; if ( pVoronoiObj == nullptr) return false ; for ( int i = 0 ; i < int( vCrvC.size()) ; i ++) { if ( ! pVoronoiObj->AddCurve( vCrvC[i])) return false ; } return pVoronoiObj->CalcMedialAxis( vCrvs, nSide) ; } //---------------------------------------------------------------------------- bool CalcCurveFatCurve( const ICurve& crvC, ICURVEPOVECTOR& vCrvs, double dRadius, bool bSquareEnds, bool bSquareMids, bool bMergeOnlySameProps) { Voronoi* pVoronoiObj = GetCurveVoronoi( crvC) ; if ( pVoronoiObj == nullptr) return false ; return pVoronoiObj->CalcFatCurve( vCrvs, dRadius, bSquareEnds, bSquareMids, bMergeOnlySameProps) ; } //---------------------------------------------------------------------------- bool CalcCurveLimitOffset( const ICurve& crvC, double& dOffs) { // se curva aperta errore if ( ! crvC.IsClosed()) return false ; Voronoi* pVoronoiObj = GetCurveVoronoi( crvC) ; if ( pVoronoiObj == nullptr) return false ; // verifico l'orientamento della curva rispetto al piano di vroni per capire se l'offset critico è quello a destra ( offset > 0) // oppure quello a sinistra ( offset < 0) Plane3d plPlane, plPlaneVroni ; double dArea ; if ( ! CurveGetArea( crvC, plPlane, dArea) || ! pVoronoiObj->GetVroniPlane( plPlaneVroni)) return false ; bool bLeft = AreSameVectorApprox( plPlane.GetVersN(), plPlaneVroni.GetVersN()) ; return pVoronoiObj->CalcLimitOffset( 0, bLeft, dOffs) ; } //---------------------------------------------------------------------------- bool CalcCurveSingleCurvesOffset( const ICurve& crvC, ICURVEPOVECTOR& vCrvs, double dOffs) { Voronoi* pVoronoiObj = GetCurveVoronoi( crvC) ; if ( pVoronoiObj == nullptr) return false ; return pVoronoiObj->CalcSingleCurvesOffset( vCrvs, dOffs) ; } //---------------------------------------------------------------------------- bool CalcOffsetCurves( const ICURVEPVECTOR& vpCrvs, ICURVEPOVECTOR& vCrvs, double dOffs, int nType) { // creo oggetto Voronoi con le curve passate PtrOwner pVoronoiObj( new( std::nothrow) Voronoi()) ; if ( pVoronoiObj == nullptr) return false ; for ( int i = 0 ; i < int( vpCrvs.size()) ; i ++) { if ( ! pVoronoiObj->AddCurve( vpCrvs[i])) return false ; } return pVoronoiObj->CalcOffset( vCrvs, dOffs, nType) ; } //---------------------------------------------------------------------------- bool CalcFatOffsetCurves( const ICURVEPVECTOR& vpCrvs, ICURVEPOVECTOR& vCrvs, double dOffs, bool bSquareEnds, bool bSquareMids, bool bMergeOnlySameProps) { // controllo validità delle curve for ( auto& pCrv : vpCrvs) { if ( pCrv == nullptr) return false ; } // se offset nullo restituisco direttamente le curve if ( abs( dOffs) < EPS_SMALL) { for ( auto& pCrv : vpCrvs) { if ( ! vCrvs.emplace_back( pCrv->Clone())) return false ; } return true ; } // creo oggetto Voronoi con le curve passate PtrOwner pVoronoiObj( new( std::nothrow) Voronoi()) ; if ( pVoronoiObj == nullptr) return false ; for ( int i = 0 ; i < int( vpCrvs.size()) ; i ++) { if ( ! pVoronoiObj->AddCurve( vpCrvs[i])) return false ; } return pVoronoiObj->CalcFatCurve( vCrvs, dOffs, bSquareEnds, bSquareMids, bMergeOnlySameProps) ; } //---------------------------------------------------------------------------- void ResetCurveVoronoi( const ICurve& crvC) { switch ( crvC.GetType()) { case CRV_LINE : GetBasicCurveLine( &crvC)->ResetVoronoiObject() ; break ; case CRV_ARC : GetBasicCurveArc( &crvC)->ResetVoronoiObject() ; break ; case CRV_BEZIER : GetBasicCurveBezier( &crvC)->ResetVoronoiObject() ; break ; case CRV_COMPO : GetBasicCurveComposite( &crvC)->ResetVoronoiObject() ; break ; default : break ; } }