EgtGeomKernel/MedialAxis.cpp

//----------------------------------------------------------------------------
//  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"

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()) ;
   for ( int i = 0 ; i < nCount ; ++ i) {
      if ( ccClass[i].nClass == CRVC_IN) {
        // 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 ;
         }
        // 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) ;
         }
      }
   }
   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 ;
}