Files
EgtGeomKernel/CurveAux.cpp
T
Daniele Bariletti 563697f840 EgtGeomKernel :
- correzione al box del volume spazzato VM 5assi
- modificata gestione conversione curve a bezier.
2025-09-26 10:47:06 +02:00

2198 lines
80 KiB
C++

//----------------------------------------------------------------------------
// 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/Extern/Eigen/Dense"
using namespace std ;
static bool FindSpan( double dU, int nDeg, const DBLVECTOR& vKnots, int& nSpan) ;
static bool CalcBasisFunc( double dU, 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<int> vProp{ crvC.GetTempProp( 0), crvC.GetTempProp( 1)} ;
if ( vProp[0] != 0 || vProp[1] != 0)
sOut += string( "TmpProp=") + ToString( vProp) + szNewLine ;
// eventuali parametri temporanei
vector<double> 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<ICurve> pBezierForm ;
switch ( pCrv->GetType()) {
case CRV_LINE :{
const ICurveLine* pCrvLine = static_cast<const ICurveLine*>( pCrv) ;
int nDegWanted = nDeg == -1 ? 3 : nDeg ;
pBezierForm.Set( LineToBezierCurve( pCrvLine, nDegWanted, bMakeRatOrNot)) ;
break ;
}
case CRV_ARC : {
const ICurveArc* pCrvArc = static_cast<const ICurveArc*>( pCrv) ;
int nDegWanted = nDeg == -1 ? 3 : nDeg ;
pBezierForm.Set( ArcToBezierCurve( pCrvArc, nDegWanted, bMakeRatOrNot)) ;
break ;
}
case CRV_COMPO : {
const ICurveComposite* pCrvCompo = static_cast<const ICurveComposite*>( pCrv) ;
pBezierForm.Set( CompositeToBezierCurve( pCrvCompo, nDeg, bMakeRatOrNot)) ;
break ;
}
case CRV_BEZIER : {
const ICurveBezier* pCrvBezier = static_cast<const ICurveBezier*>( 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<ICurveBezier> 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<ICurveBezier> 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<ICurveComposite> 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<ICurveBezier> 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<ICurveComposite> pCCBezier( CreateCurveComposite()) ;
for ( int i = 0 ; i < int( pCC->GetCurveCount()) ; ++i) {
PtrOwner<ICurve> 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<ICurveBezier> 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<ICurveComposite> pBezCubics( CreateCurveComposite()) ;
pBezCubics->AddCurve( ApproxBezierWithCubics(pCrvBezier, dTol)) ;
if ( IsNull( pBezCubics))
return nullptr ;
// adatto ogni sottocurva cubica
PtrOwner<ICurveComposite> 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<ICurveComposite> pBezCubics( CreateCurveComposite()) ;
pBezCubics->AddCurve( ApproxBezierWithCubics(pCrvBezier, dTol)) ;
if ( IsNull( pBezCubics) || ! pBezCubics->IsValid())
return nullptr ;
// adatto ogni sottocurva cubica
PtrOwner<ICurveComposite> 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<ICurveBezier> 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<ICurveBezier> 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
CalcBezierApproxError( 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<ICurveLine> 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<CurveLine> pCL1( CreateBasicCurveLine()) ;
pCL1->Set( ptStart, ptAux1) ;
PtrOwner<CurveLine> 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<ICurveComposite> pCC( CreateBasicCurveComposite()) ;
// approssimo l'ultimo tratto
int nCount = 0 ;
while ( dErr > dTol && nCount < 100) {
PtrOwner<ICurveBezier> pCrvPart( pCrvBezier->Clone()) ;
pCrvPart->TrimStartEndAtParam( double( nParts - 1) / nParts, 1) ;
PtrOwner<ICurveBezier> 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)) ;
}
}
CalcBezierApproxError( 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<ICurveBezier> pCrvPart( pCrvBezier->Clone()) ;
pCrvPart->TrimStartEndAtParam( double( i - 1) / nParts, double(i) / nParts) ;
PtrOwner<ICurveBezier> 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<ICurveBezier> 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<ICurveBezier> 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) ;
}
//----------------------------------------------------------------------------
ICurve*
ApproxCurveWithBezier( const ICurve*, double dTol)
{
// campiono punti lungo la curva e poi li interpolo
// oppure faccio la fat curve e poi calcolo una bezier che stia all'interno di quella regione
PtrOwner<ICurveComposite> pCC( CreateBasicCurveComposite()) ;
return Release( pCC) ;
}
//----------------------------------------------------------------------------
ICurve*
InterpolatePointSetWithBezier( const PNTVECTOR& vPnt, double dTol)
{
PtrOwner<ICurve> pCrvInt ;
// pag 364 del Piegl
// scelgo il parametro associato ad ogni punto in modo che rispecchi la distanza tra i punti
int nPoints = vPnt.size() ;
if ( nPoints < 4)
return nullptr ;
int nDeg = 3 ;
DBLVECTOR vLen ;
double dLenTot = 0 ;
for ( int i = 0 ; i < int( nPoints -1) ; ++i) {
double dLen = Dist(vPnt[i], vPnt[i+1]) ;
vLen.push_back( dLen) ;
dLenTot += dLen ;
}
DBLVECTOR vPntParam ;
vPntParam.resize( nPoints) ;
vPntParam[0] = 0 ;
vPntParam.back() = 1 ;
for ( int i = 1 ; i < int( 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 ; j < i + nDeg ; ++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) {
int nSpan = 0 ; FindSpan(vPntParam[i], nDeg, vKnots, nSpan) ;
if ( i == 0)
mA.row(0).col(0) << 1 ;
else if ( i == nPoints - 1)
mA.row(i).col(nPoints - 1) << 1 ;
else {
DBLVECTOR vBasis ; vBasis.resize( nDeg + 1) ;
CalcBasisFunc( vPntParam[i], nDeg, vKnots, vBasis) ;
for( int j = nSpan - nDeg ; j < nSpan ; ++j)
mA.row(i).col(j) << vBasis[j - nSpan + nDeg] ;
}
}
int nDim = 3 ;
Eigen::MatrixXd mb ; mb.resize( nPoints, nDim) ;
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)) ;
return Release( pCrvInt) ;
}
//----------------------------------------------------------------------------
static bool
FindSpan( double dU, int nDeg, const DBLVECTOR& vKnots, int& nSpan)
{
if ( dU < 0)
return false ;
else if ( dU < EPS_ZERO) {
nSpan = nDeg ;
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 ;
return true ;
}
int nLow = nDeg ;
int nHigh = nKnots ;
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)
break ;
}
nSpan = nMid ;
return true ;
}
//----------------------------------------------------------------------------
static bool
CalcBasisFunc( double dU, 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 ;
int nSpan = 0 ;
FindSpan( dU, nDeg, vKnots, nSpan) ;
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 dTemp = vBasis[r] / ( vRight[r+1] + vLeft[j-r]) ;
vBasis[r] = dSaved + vRight[r+1] * dTemp ;
dSaved = vLeft[j-r] * dTemp ;
}
vBasis[j] = dSaved ;
}
return true ;
}
//----------------------------------------------------------------------------
bool
CalcBezierApproxError( const ICurveBezier* pCrvOri, const ICurveBezier* pCrvNew, double& dErr, int nPoints)
{
if ( pCrvOri->GetType() != CRV_BEZIER || pCrvNew->GetType() != CRV_BEZIER)
return false ;
// 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<CurveComposite> 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<CurveComposite> 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<double>( 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<double>( 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<double>( cnData.vU.begin() + 2, cnData.vU.end()) ;
}
// sennò ne tolgo uno all'inizio e uno alla fine
else
cnData.vU = vector<double>( 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<double>( 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<double&>( 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<double&>( 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<CurveBezier> 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<CurveComposite> 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<CurveBezier> 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<CurveBezier> 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<ICurve> 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<ICurve> 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<ICurve> 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<ICurve> 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<Voronoi> 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<Voronoi> 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<Voronoi> 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<Voronoi> 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 ;
}
}