Files
EgtGeomKernel/MedialAxis.cpp
T
Riccardo Elitropi bd80f16812 EgtGeomKernel :
- Modifica Medial Axis.
2022-12-15 16:24:58 +01:00

238 lines
7.4 KiB
C++

//----------------------------------------------------------------------------
// EgalTech 2022-2022
//----------------------------------------------------------------------------
// File : MedialAxis.cpp Data : 22.08.22 Versione : 2.4h2
// Contenuto : Funzione calcolo approssimato MedialAxis principale.
//
//
//
// Modifiche : 22.08.22 DS Creazione modulo.
//
//
//
//----------------------------------------------------------------------------
//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "GeoConst.h"
#include "CurveLine.h"
#include "CurveArc.h"
#include "CurveBezier.h"
#include "CurveComposite.h"
#include "/EgtDev/Include/EGkBBox3d.h"
#include "/EgtDev/Include/EGkIntersCurves.h"
#include "/EgtDev/Include/EGkMedialAxis.h"
#include <EgtDev/Include/EGnStringUtils.h>
using namespace std ;
//----------------------------------------------------------------------------
static bool
ProcessCircleCurve( const Point3d& ptCen, double dDiam, const CurveComposite& crvCompo,
PNTVECTOR& vPnt, DBLVECTOR& vDiam)
{
// pulizia risultati
vPnt.clear() ;
vDiam.clear() ;
// creazione circonferenza
CurveArc crvCirc ;
crvCirc.Set( ptCen, Z_AX, dDiam / 2) ;
// classificazioni parti di circonferenza rispetto alla curva
IntersCurveCurve LnCmpInt( crvCirc, crvCompo) ;
CRVCVECTOR ccClass ;
if ( ! LnCmpInt.GetCurveClassification( 0, 10 * EPS_SMALL, ccClass))
return false ;
int nCount = int( ccClass.size()) ;
// ... per evitare percorsi di ritorno
int nInparts = 0 ;
for ( int i = 0 ; i < nCount ; ++ i) {
if ( ccClass[i].nClass == CRVC_IN) {
nInparts ++ ;
// gestione unione intervallo a cavallo del parametro 0==1
if ( i == nCount - 1 && ccClass[0].nClass == CRVC_IN && ccClass[0].dParS < EPS_ZERO && ccClass[i].dParE > 1 - EPS_ZERO) {
double dPar = ( ccClass[i].dParS + 1 + ccClass[0].dParE) / 2 ;
if ( dPar > 1)
dPar -= 1 ;
Point3d ptP ;
crvCirc.GetPointD1D2( dPar, ICurve::FROM_MINUS, ptP) ;
vPnt[0] = ptP ;
double dDelta = ccClass[i].dParE - ccClass[i].dParS + ccClass[0].dParE - ccClass[0].dParS ;
vDiam[0] = PIGRECO * dDiam * dDelta ;
nInparts-- ;
}
// altrimenti caso standard
else {
double dPar = ( ccClass[i].dParS + ccClass[i].dParE) / 2 ;
Point3d ptP ;
crvCirc.GetPointD1D2( dPar, ICurve::FROM_MINUS, ptP) ;
vPnt.emplace_back( ptP) ;
double dDelta = ccClass[i].dParE - ccClass[i].dParS ;
vDiam.emplace_back( PIGRECO * dDiam * dDelta) ;
}
}
}
if ( nInparts == 1)
vPnt.clear() ;
return true ;
}
//----------------------------------------------------------------------------
bool
CurveSimpleMedialAxis( const ICurve* pCrv, PolyLine& PL)
{
// pulizia della polilinea
PL.Clear() ;
// verifico che la curva esista
if ( pCrv == nullptr)
return false ;
// verifico che la curva sia chiusa
if ( ! pCrv->IsClosed())
return false ;
// verifico che la curva sia piana
Plane3d plPlane ;
if ( ! pCrv->IsFlat( plPlane, false, 10 * EPS_SMALL))
return false ;
// eventuali trasformazioni in composita
CurveComposite crvCompo ;
if ( pCrv->GetType() == CRV_ARC) {
const CurveArc* pArc = GetBasicCurveArc( pCrv) ;
crvCompo.AddCurve( *pArc) ;
}
else if ( pCrv->GetType() == CRV_BEZIER) {
const CurveBezier* pBez = GetBasicCurveBezier( pCrv) ;
PolyArc PA ;
if ( ! pBez->ApproxWithArcs( 10 * EPS_SMALL, ANG_TOL_STD_DEG, PA) || ! crvCompo.FromPolyArc( PA))
return false ;
}
else if ( pCrv->GetType() == CRV_COMPO) {
crvCompo.AddCurve( *pCrv) ;
}
else
return false ;
// se è una circonferenza
Point3d ptCen ; Vector3d vtN ; double dRad ; bool bCCW ;
if ( crvCompo.IsACircle( 10 * EPS_SMALL, ptCen, vtN, dRad, bCCW)) {
PL.AddUPoint( 0, ptCen) ;
return true ;
}
// la porto nel riferimento del piano
Frame3d frLoc ;
frLoc.Set( ORIG, plPlane.GetVersN()) ;
crvCompo.ToLoc( frLoc) ;
// lo oriento in senso antiorario
double dArea ;
if ( ! crvCompo.GetAreaXY( dArea))
return false ;
if ( dArea < 0)
crvCompo.Invert() ;
// ne calcolo il box
BBox3d b3Compo ;
if ( ! crvCompo.GetLocalBBox( b3Compo, BBF_STANDARD))
return false ;
// ricerca di un punto di partenza
CurveLine crvLine ;
if ( b3Compo.GetDimX() > b3Compo.GetDimY()) {
Point3d ptP1 = Media( b3Compo.GetMin(), b3Compo.GetMax()) - 0.6 * b3Compo.GetDimY() * Y_AX ;
Point3d ptP2 = ptP1 + 1.2 * b3Compo.GetDimY() * Y_AX ;
crvLine.Set( ptP1, ptP2) ;
}
else {
Point3d ptP1 = Media( b3Compo.GetMin(), b3Compo.GetMax()) - 0.6 * b3Compo.GetDimX() * X_AX ;
Point3d ptP2 = ptP1 + 1.2 * b3Compo.GetDimX() * X_AX ;
crvLine.Set( ptP1, ptP2) ;
}
IntersCurveCurve LnCmpInt( crvLine, crvCompo) ;
if ( LnCmpInt.GetCrossIntersCount() < 2)
return false ;
IntCrvCrvInfo aInfo[2] ;
LnCmpInt.GetIntCrvCrvInfo( 0, aInfo[0]) ;
LnCmpInt.GetIntCrvCrvInfo( 1, aInfo[1]) ;
Point3d ptOrig = Media( aInfo[0].IciA->ptI, aInfo[1].IciA->ptI) ;
double dDiamOrig = 1.2 * DistXY( aInfo[0].IciA->ptI, aInfo[1].IciA->ptI) ;
Vector3d vtTg0 ; Point3d ptP0 ; crvCompo.GetPointTang( aInfo[0].IciB->dU, ICurve::FROM_MINUS, ptP0, vtTg0) ;
Vector3d vtTg1 ; Point3d ptP1 ; crvCompo.GetPointTang( aInfo[1].IciB->dU, ICurve::FROM_MINUS, ptP1, vtTg1) ;
Vector3d vtDirOrig = Media( vtTg0, -vtTg1) ; vtDirOrig.Normalize() ;
// movimento lungo un ramo
Point3d ptStart = ptOrig ;
double dDiam = dDiamOrig ;
Vector3d vtDir = vtDirOrig ;
int i = 0 ;
while ( i < 1000) {
PNTVECTOR vPnt ;
DBLVECTOR vDiam ;
if ( ! ProcessCircleCurve( ptStart, dDiam, crvCompo, vPnt, vDiam))
return false ;
if ( vPnt.empty())
break ;
int nInd = -1 ;
double dCosMax = -0.5 ;
for ( int j = 0 ; j < int( vPnt.size()) ; ++ j) {
Vector3d vtNewDir = vPnt[j] - ptStart ;
if ( vtNewDir.Normalize()) {
double dCosA = vtDir * vtNewDir ;
if ( dCosA > dCosMax) {
nInd = j ;
dCosMax = dCosA ;
}
}
}
if ( nInd < 0)
break ;
vtDir = vPnt[nInd] - ptStart ;
vtDir.Normalize() ;
ptStart = vPnt[nInd] ;
dDiam = 1.2 * vDiam[nInd] ;
PL.AddUPoint( 0, GetToGlob( ptStart, frLoc)) ;
++ i ;
}
// movimento lungo l'altro ramo
ptStart = ptOrig ;
dDiam = dDiamOrig ;
vtDir = -vtDirOrig ;
i = 0 ;
while ( i < 1000) {
PNTVECTOR vPnt ;
DBLVECTOR vDiam ;
if ( ! ProcessCircleCurve( ptStart, dDiam, crvCompo, vPnt, vDiam))
return false ;
if ( vPnt.empty())
break ;
int nInd = -1 ;
double dCosMax = -0.5 ;
for ( int j = 0 ; j < int( vPnt.size()) ; ++ j) {
Vector3d vtNewDir = vPnt[j] - ptStart ;
if ( vtNewDir.Normalize()) {
double dCosA = vtDir * vtNewDir ;
if ( dCosA > dCosMax) {
nInd = j ;
dCosMax = dCosA ;
}
}
}
if ( nInd < 0)
break ;
vtDir = vPnt[nInd] - ptStart ;
vtDir.Normalize() ;
ptStart = vPnt[nInd] ;
dDiam = 1.2 * vDiam[nInd] ;
PL.AddUPoint( 0, GetToGlob( ptStart, frLoc), false) ;
++ i ;
}
return true ;
}