Egalware/EgtGeomKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

EgtGeomKernel :

Browse Source 
- implementata l'intersezione tra linea e sup bezier.
This commit is contained in:
Daniele Bariletti
2024-02-12 12:52:56 +01:00
parent 2c159d7ce6
commit 0f8012ce61
3 changed files with 199 additions and 104 deletions
Download Patch File Download Diff File Expand all files Collapse all files
IntersLineSurfBez.cpp
+196 -97
View File
@@ -13,108 +13,21 @@
//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "IntersLineBox.h"
#include "/EgtDev/Include/EGkIntersLineTria.h"
#include "/EgtDev/Include/EGkIntersLineSurfTm.h"
#include "/EgtDev/Include/EGkIntersLineSurfBez.h"
#include "/EgtDev/Include/EGkSurfBezier.h"
#include "/EgtDev/Extern/Eigen/Dense"
#include "DistPointLine.h"
#include "CurveLine.h"
using namespace std ;
//----------------------------------------------------------------------------
// Intersezione di una linea con una superficie TriMesh
//----------------------------------------------------------------------------
bool
IntersLineSurfBz( const Point3d& ptL, const Vector3d& vtL, double dLen, const ISurfBezier& SBz,
ILSBIVECTOR& vInfo, bool bFinite)
{
// verifico linea
Vector3d vtDir = vtL ;
if ( ! vtDir.Normalize( EPS_ZERO))
return false ;
// verifico superficie
if ( &SBz == nullptr)
return false ;
// verifico parametro di ritorno
if ( &vInfo == nullptr)
return false ;
vInfo.clear() ;
// trovo le intersezioni con la trimesh ausiliaria
const ISurfTriMesh* pSurfTm = SBz.GetAuxSurf() ;
ILSIVECTOR vInfoTm ;
if ( ! IntersLineSurfTm( ptL, vtL, dLen, *pSurfTm, vInfoTm, bFinite))
return false ;
// ricavo le intersezioni con la superficie di Bezier
for ( IntLinStmInfo InfoTm : vInfoTm ) {
// devo raffinare i parametri lungo la curva, l'angolo e i punti di intersezione
double dUU, dUU2, dCos ;
Point3d ptI, ptI2 ;
// devo trovare le intersezioni
Point3d ptSP, ptSP2 ;
RefinePoints( pSurfTm, InfoTm.nILTT, InfoTm.nT, ptI, ptI2,ptSP, ptSP2) ;
IntLinSbzInfo InfoBz( InfoTm.nILTT, dUU, dUU2, InfoTm.nT, dCos, InfoTm.ptI, InfoTm.ptI2, ptSP, ptSP2) ;
vInfo.emplace_back( InfoBz) ;
}
//////////////////////////////////////////////////////// interlineSurfTm
//// limito la linea al box dei triangoli della superficie
// BBox3d b3Stm = Stm.GetAllTriaBox() ;
// if ( b3Stm.IsEmpty())
// return false ;
//// lo ingrandisco per non avere problemi con faccia piana su piani canonici
// b3Stm.Expand( 10 * EPS_SMALL) ;
// double dU1, dU2 ;
// if ( ! IntersLineBox( ptL, vtL, b3Stm.GetMin() , b3Stm.GetMax(), dU1, dU2))
// return true ;
// if ( bFinite) {
// dU1 = max( dU1, 0.) ;
// dU2 = min( dU2, dLen) ;
// if ( dU2 - dU1 < EPS_SMALL)
// return true ;
// }
// Point3d ptStart = ptL + dU1 * vtL ;
// double dLenEff = dU2 - dU1 ;
//// cerco i triangoli intersecati dalla linea
// const double BOX_STEP = 10 ;
// int nStep = int( ceil( dLenEff / BOX_STEP)) ;
// Vector3d vtStep = dLenEff / nStep * vtL ;
// INTVECTOR vPrevT ;
// for ( int i = 0 ; i < nStep ; ++ i) {
// BBox3d b3Box( ptStart + i * vtStep, ptStart + ( i + 1) * vtStep) ;
// INTVECTOR vT ;
// if ( Stm.GetAllTriaOverlapBox( b3Box, vT)) {
// for ( auto nT : vT) {
// if ( find( vPrevT.begin(), vPrevT.end(), nT) == vPrevT.end()) {
// vPrevT.emplace_back( nT) ;
// Triangle3d Tria ;
// Stm.GetTriangle( nT, Tria) ;
// // aggiorno info con intersezione
// UpdateInfoIntersLineSurfTm( ptL, vtDir, dLen, nT, Tria, vInfo, bFinite) ;
// }
// }
// }
// }
//// ordino il vettore delle eventuali intersezioni secondo il senso crescente del parametro di linea
// OrderInfoIntersLineSurfTm( vInfo) ;
// return true ;
//////////////////////////////////////////////////////// interlineSurfTm
return true ;
}
//----------------------------------------------------------------------------
// raffino i punti di intersezione e recupero le rispettive coordinate nello spazio parametrico
bool
RefinePoints( const ISurfTriMesh* pSurfTm ,int nIL, int nT, Point3d& ptI, Point3d& ptI2, Point3d& ptSP, Point3d& ptSP2) {
FindParametricCoord( const ISurfTriMesh* pSurfTm ,int nIL, int nT, const Point3d& ptI, Point3d& ptSP) {
// recupero i dati dei vertici del triangolo che fa intersezione
int nVert[3] ;
pSurfTm->GetTriangle( nT, nVert) ;
double dU0, dV0, dU1, dV1,dU2, dV2 ;
@@ -125,12 +38,198 @@ RefinePoints( const ISurfTriMesh* pSurfTm ,int nIL, int nT, Point3d& ptI, Point3
pSurfTm->GetVertex( nVert[0], pt0) ;
pSurfTm->GetVertex( nVert[1], pt1) ;
pSurfTm->GetVertex( nVert[2], pt2) ;
Eigen::Matrix3d m3x3 ;
//m3x3 << pt0.x, pt1.x, pt2.x, pt0.y, pt1.y, pt2.y, pt0.z, pt1.z, pt2.z ;
m3x3.col(0) << pt0.x, pt0.y , pt0.z ;
m3x3.col(1) << pt1.x, pt1.y , pt1.z ;
m3x3.col(2) << pt2.x, pt2.y , pt2.z ;
// calcolo approssimativamente le coordinate nello spazio parametrico del punto di intersezione
// quindi prima calcolo la composizione lineare tra i vertici del triangolo per ottenere il punto di intersezione
Eigen::Matrix3d mA ;
mA.col(0) << pt0.x, pt0.y , pt0.z ;
mA.col(1) << pt1.x, pt1.y , pt1.z ;
mA.col(2) << pt2.x, pt2.y , pt2.z ;
Eigen::Vector3d b ( ptI.x, ptI.y, ptI.z) ;
Eigen::Vector3d x = m3x3.fullPivLu().solve(b) ;
Eigen::Vector3d x = mA.fullPivLu().solve(b) ;
// applico questa composizione alle loro coordinate parametriche
Eigen::Matrix3d mB ;
mB.col(0) << dU0, dV0, 0 ;
mB.col(1) << dU1, dV1, 0 ;
mB.col(2) << dU2, dV2, 0 ;
Eigen::Vector3d ptParam = mB * x ;
ptSP.x = ptParam.x() ;
ptSP.y = ptParam.y() ;
return true ;
}
//----------------------
bool
RefineIntersNewton( const Point3d& ptL, const Vector3d& vtL, double dLen, bool bFinite,
const ISurfBezier* pSurfBz, Point3d& ptSP, Point3d& ptIBz) {
// la funzione raffina la posisione del punto ptSP, minimizzando la distanza dalla retta e restituisce il punto di intersezione ptIBz
pSurfBz->GetPointD1D2( ptSP.x / SBZ_TREG_COEFF, ptSP.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBz) ;
// usando un algoritmo di newton cerco di avvicinarmi il più possibile alla retta
DistPointLine dpl( ptIBz, ptL, vtL, dLen, bFinite) ;
double dDistNew = 0, dDistPre = 0 ;
dpl.GetDist(dDistNew) ;
int nCount = 0 ;
double dh = EPS_SMALL ;
pSurfBz->GetPointD1D2( ptSP.x, ptSP.y, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBz) ;
// metodo di newton in più dimensioni
// vario sia il parametro U che il parametro V e verifico se la distanza dalla retta diminuisce per scostamenti positivi o negativi.
while ( dDistNew > EPS_SMALL && nCount < 100) {
dDistPre = dDistNew ;
Point3d ptIBzNew1 ;
pSurfBz->GetPointD1D2( ( ptSP.x + dh) / SBZ_TREG_COEFF, ptSP.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew1) ;
DistPointLine dplNewU( ptIBzNew1, ptL, vtL, dLen, bFinite) ;
dplNewU.GetDist( dDistNew) ;
double dfdU = ( dDistNew - dDistPre) / dh ;
Point3d ptIBzNew2 ;
pSurfBz->GetPointD1D2( ptSP.x / SBZ_TREG_COEFF, ( ptSP.y + dh) / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew2) ;
DistPointLine dplNewV( ptIBzNew2, ptL, vtL, dLen, bFinite) ;
dplNewV.GetDist( dDistNew) ;
double dfdV = ( dDistNew - dDistPre) / dh ;
//// opzione 0
////scelgo h1 e h2 separatamente e in modo da annullare f(x)
//// opzione 1
//// valore fisso
//double dr = EPS_SMALL ;
//if ( dDistPre > 1)
// dr = 1 ;
//else if ( dDistPre > 0.1)
// dr = 0.1 ;
//else if ( dDistPre > 0.01)
// dr = 0.01 ;
//// opzione 2
//// valore direttamente vincolato
//double dr = dDistPre ;
//// opzione 3
//// valuto la deformazione locale in base allo spostamento del punto sulla bezier // non serve
//double dh1 = Dist( ptIBz, ptIBzNew1) ;
//double dh2 = Dist( ptIBz, ptIBzNew2) ;
// potrei valutare il nuovo spostamento in base all'ultima variazione di dDist
// potrei anche vedere se sto uscendo dal triangolo ( definito nello spazio parametrico)
// mi avvicino cercando di annullare la distanza in un colpo solo
double dr = - dDistPre / ( dfdU + dfdV) ;
pSurfBz->GetPointD1D2(( ptSP.x + dr * dfdU) / SBZ_TREG_COEFF, ( ptSP.y + dr * dfdV) / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBz) ;
DistPointLine dplNew( ptIBz, ptL, vtL, dLen, bFinite) ;
dplNew.GetDist( dDistNew) ;
++nCount ;
}
return nCount != 99 ;
}
////----------------------------------------------------------------------------
//bool
//RefineIntersBisec( const Point3d& ptL, const Vector3d& vtL, double dLen, bool bFinite,
// const ISurfBezier* pSurfBz, Point3d& ptSP, Point3d& ptIBz) {
//
//}
//----------------------------------------------------------------------------
void
UpdateInfoIntersLineSurfBz( const Point3d& ptL, const Vector3d& vtDir, int nILT, int nT, const Point3d& ptSP, const Point3d& ptIBz, double dCos,
const Point3d& ptSP2, const Point3d& ptIBz2, double dCos2, ILSBIVECTOR& vInfo)
{
if ( nILT == ILTT_IN || nILT == ILTT_EDGE || nILT == ILTT_VERT) {
double dU = ( ptIBz - ptL) * vtDir ;
vInfo.emplace_back( nILT, dU, nT, dCos, ptIBz, ptSP) ;
}
else if ( nILT == ILTT_SEGM || nILT == ILTT_SEGM_ON_EDGE) {
double dU = ( ptIBz - ptL) * vtDir ;
double dU2 = ( ptIBz2 - ptL) * vtDir ;
vInfo.emplace_back( nILT, dU, dU2, nT, dCos2, ptIBz, ptIBz2, ptSP, ptSP2) ;
}
}
//----------------------------------------------------------------------------
void
OrderInfoIntersLineSurfBz( ILSBIVECTOR& vInfo)
{
// se non trovati, esco
if ( vInfo.size() == 0)
return ;
// ordino il vettore delle intersezioni secondo il senso crescente del parametro di linea
sort( vInfo.begin(), vInfo.end(),
[]( const IntLinSbzInfo& a, const IntLinSbzInfo& b)
{ double dUa = ( ( a.nILTT == ILTT_SEGM || a.nILTT == ILTT_SEGM_ON_EDGE) ? ( a.dU + a.dU2) / 2 : a.dU) ;
double dUb = ( ( b.nILTT == ILTT_SEGM || b.nILTT == ILTT_SEGM_ON_EDGE) ? ( b.dU + b.dU2) / 2 : b.dU) ;
return ( dUa < dUb) ; }) ;
}
//----------------------------------------------------------------------------
// Intersezione di una linea con una superficie TriMesh
//----------------------------------------------------------------------------
bool
IntersLineSurfBz( const Point3d& ptL, const Vector3d& vtL, double dLen, const ISurfBezier* pSurfBz,
ILSBIVECTOR& vInfo, bool bFinite)
{
PtrOwner<ICurveLine> pCL( CreateCurveLine()) ;
pCL->SetPVL(ptL, vtL, dLen) ;
// verifico linea
Vector3d vtDir = vtL ;
if ( ! vtDir.Normalize( EPS_ZERO))
return false ;
// verifico superficie
if ( pSurfBz == nullptr)
return false ;
// verifico parametro di ritorno
if ( &vInfo == nullptr)
return false ;
vInfo.clear() ;
// trovo le intersezioni con la trimesh ausiliaria
const ISurfTriMesh* pSurfTm = pSurfBz->GetAuxSurf() ;
ILSIVECTOR vInfoTm ;
if ( ! IntersLineSurfTm( ptL, vtL, dLen, *pSurfTm, vInfoTm, bFinite))
return false ;
// ricavo le intersezioni con la superficie di Bezier
for ( IntLinStmInfo InfoTm : vInfoTm ) {
// devo raffinare i parametri lungo la curva, l'angolo e i punti di intersezione
Point3d ptI, ptI2 ;
// devo trovare le intersezioni
Point3d ptSP, ptSP2 ; // coordinate parametriche delle soluzioni
FindParametricCoord( pSurfTm, InfoTm.nILTT, InfoTm.nT, InfoTm.ptI, ptSP) ;
Point3d ptIBz, ptIBz2 ;
if ( ! RefineIntersNewton( ptL, vtL, dLen, bFinite, pSurfBz, ptSP, ptIBz)) {
/////// posso provare anche a rilanciare newton con un punto di partenza diverso oppure con una direzione di avvicinamento diversa///////////////////////////////////
// per restare nel triangolo mi sposto verso un vertice
int nVert[3] ;
pSurfTm->GetTriangle( InfoTm.nT, nVert) ;
double dU0, dV0 ;
pSurfTm->GetVertexParam( nVert[0], dU0, dV0) ;
ptSP = ptSP + Point3d(dU0, dV0, 0) ;
if ( ! RefineIntersNewton( ptL,vtL, dLen, bFinite, pSurfBz, ptSP, ptIBz))
return false ;
}
Vector3d vtN ;
pSurfBz->GetPointNrmD1D2(ptSP.x / SBZ_TREG_COEFF, ptSP.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBz, vtN) ;
double dCos = vtN * vtL ;
double dCos2 = 0 ;
// eventualmente ripeto tutto per ptI2 ( se ho un'intersezione con sovrapposizione)
if ( InfoTm.nILTT == ILTT_SEGM || InfoTm.nILTT == ILTT_SEGM_ON_EDGE ) {
FindParametricCoord( pSurfTm, InfoTm.nILTT, InfoTm.nT, InfoTm.ptI2, ptSP2) ;
if ( !RefineIntersNewton(ptL, vtL, dLen, bFinite, pSurfBz, ptSP2, ptIBz2) ) {
int nVert[3] ;
pSurfTm->GetTriangle( InfoTm.nT, nVert) ;
double dU0, dV0 ;
pSurfTm->GetVertexParam( nVert[0], dU0, dV0) ;
ptSP = ptSP + Point3d(dU0, dV0, 0) ;
if ( ! RefineIntersNewton( ptL,vtL, dLen, bFinite, pSurfBz, ptSP, ptIBz))
return false ;
}
pSurfBz->GetPointNrmD1D2( ptSP2.x / SBZ_TREG_COEFF, ptSP2.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBz2, vtN) ;
dCos2 = vtN * vtL ;
}
//IntLinSbzInfo InfoBz( InfoTm.nILTT, dUU, dUU2, InfoTm.nT, dCos, ptIBz, ptIBz2, ptSP, ptSP2) ;
//vInfo.emplace_back( InfoBz) ;
UpdateInfoIntersLineSurfBz( ptL, vtL, InfoTm.nILTT, InfoTm.nT, ptSP, ptIBz, dCos, ptSP2, ptIBz2, dCos2, vInfo) ;
}
OrderInfoIntersLineSurfBz( vInfo) ;
//////////////////////////////////////////////////////// interlineSurfTm
//UpdateInfoIntersLineSurfTm( ptL, vtDir, dLen, nT, Tria, vInfo, bFinite) ;
// OrderInfoIntersLineSurfTm( vInfo) ;
//////////////////////////////////////////////////////// interlineSurfTm
return true ;
}