diff --git a/EgtGeomKernel.vcxproj.filters b/EgtGeomKernel.vcxproj.filters
index c24cae1..b22fe9d 100644
--- a/EgtGeomKernel.vcxproj.filters
+++ b/EgtGeomKernel.vcxproj.filters
@@ -273,6 +273,9 @@
File di origine\GeoInters
+
+ File di origine\GeoInters
+
File di origine\GeoCreate
diff --git a/IntersLineSurfBez.cpp b/IntersLineSurfBez.cpp
index 3637788..34e45fe 100644
--- a/IntersLineSurfBez.cpp
+++ b/IntersLineSurfBez.cpp
@@ -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 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 ;
}
\ No newline at end of file
diff --git a/SurfTriMesh.cpp b/SurfTriMesh.cpp
index 19157d9..8bac766 100644
--- a/SurfTriMesh.cpp
+++ b/SurfTriMesh.cpp
@@ -505,13 +505,6 @@ SurfTriMesh::GetFirstVertex( Point3d& ptP) const
return GetNextVertex( SVT_NULL, ptP) ;
}
-//----------------------------------------------------------------------------
-int
-SurfTriMesh::GetFirstVertexParam( int nId, double& dU, double& dV) const
-{
- return GetNextVertexParam( SVT_NULL, dU, dV) ;
-}
-
//----------------------------------------------------------------------------
int
SurfTriMesh::GetNextVertex( int nId, Point3d& ptP) const