diff --git a/EgtGeomKernel.vcxproj b/EgtGeomKernel.vcxproj index 6fe399c..e776a02 100644 --- a/EgtGeomKernel.vcxproj +++ b/EgtGeomKernel.vcxproj @@ -309,6 +309,7 @@ copy $(TargetPath) \EgtProg\Dll64 + @@ -450,6 +451,7 @@ copy $(TargetPath) \EgtProg\Dll64 + diff --git a/EgtGeomKernel.vcxproj.filters b/EgtGeomKernel.vcxproj.filters index 4802df7..5b7bed1 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 new file mode 100644 index 0000000..34e45fe --- /dev/null +++ b/IntersLineSurfBez.cpp @@ -0,0 +1,235 @@ +//---------------------------------------------------------------------------- +// EgalTech 2024 +//---------------------------------------------------------------------------- +// File : IntersLineSurfBez.cpp Data : 06.02.24 Versione : 2.6b1 +// Contenuto : Implementazione della intersezione linea/superficie bezier. +// +// +// +// Modifiche : 06.02.24 DB Creazione modulo. +// +// +//---------------------------------------------------------------------------- + +//--------------------------- Include ---------------------------------------- +#include "stdafx.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 ; + +//---------------------------------------------------------------------------- +// raffino i punti di intersezione e recupero le rispettive coordinate nello spazio parametrico +bool +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 ; + pSurfTm->GetVertexParam( nVert[0], dU0, dV0) ; + pSurfTm->GetVertexParam( nVert[1], dU1, dV1) ; + pSurfTm->GetVertexParam( nVert[2], dU2, dV2) ; + Point3d pt0, pt1, pt2 ; + pSurfTm->GetVertex( nVert[0], pt0) ; + pSurfTm->GetVertex( nVert[1], pt1) ; + pSurfTm->GetVertex( nVert[2], pt2) ; + // 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 = 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/IntersLineSurfBez.h b/IntersLineSurfBez.h new file mode 100644 index 0000000..6f70f09 --- /dev/null +++ b/IntersLineSurfBez.h @@ -0,0 +1 @@ +#pragma once