diff --git a/DistLineLine.cpp b/DistLineLine.cpp new file mode 100644 index 0000000..8451373 --- /dev/null +++ b/DistLineLine.cpp @@ -0,0 +1,220 @@ +// File : DistLineLine.h Data : Versione : +// Contenuto : Dichiarazione della classe distanza fra elementi lineari. +// +// +// +// Modifiche : +// +// +//---------------------------------------------------------------------------- + +#include "stdafx.h" +#include "DistLineLine.h" +#include "/EgtDev/Include/EgtNumUtils.h" +#include "/EgtDev/Include/EGkGeoCollection.h" +#include "/EgtDev/Include/EGkGeoConst.h" +#include + +//---------------------------------------------------------------------------- +DistLineLine::DistLineLine( const Point3d& ptSt1, const Point3d& ptEn1, + const Point3d& ptSt2, const Point3d& ptEn2, + bool bIsSegment1, bool bIsSegment2) +{ + Vector3d vtD1 = ptEn1 - ptSt1 ; + double dLen1 = vtD1.Len() ; + vtD1 /= dLen1 ; + Vector3d vtD2 = ptEn2 - ptSt2 ; + double dLen2 = vtD1.Len() ; + vtD2 /= dLen2 ; + DistLineLine( ptSt1, vtD1, dLen1, ptSt2, vtD2, dLen2, bIsSegment1, bIsSegment2) ; +} + +//---------------------------------------------------------------------------- +// I vettori devono essere normalizzati +DistLineLine::DistLineLine( const Point3d& ptSt1, const Vector3d& vtD1, double dLen1, + const Point3d& ptSt2, const Vector3d& vtD2, double dLen2, + bool bIsSegment1, bool bIsSegment2) +{ + if ( dLen1 < EPS_SMALL || dLen2 < EPS_SMALL) { + m_dSqDist = - 1 ; + m_dDist = - 1 ; + return ; + } + Calculate( ptSt1, vtD1, dLen1, ptSt2, vtD2, dLen2, bIsSegment1, bIsSegment2) ; +} + +//---------------------------------------------------------------------------- +bool +DistLineLine::GetSqDist( double& dSqDist) +{ + if ( m_dSqDist < 0) + return false ; + dSqDist = m_dSqDist ; + return true ; +} + +//---------------------------------------------------------------------------- +bool +DistLineLine::GetDist( double& dDist) +{ + if ( m_dSqDist < 0) + return false ; + dDist = sqrt( m_dSqDist) ; + return true ; +} + +//---------------------------------------------------------------------------- +bool +DistLineLine::GetMinDistPoints( Point3d& ptMinDist1, Point3d& ptMinDist2) +{ + if ( m_dSqDist < 0) + return false ; + ptMinDist1 = m_ptMinDist1 ; + ptMinDist2 = m_ptMinDist2 ; + return true ; +} + +//---------------------------------------------------------------------------- +bool +DistLineLine::GetParamsAtMinDistPoints( double& dPar1, double& dPar2) +{ + if ( m_dSqDist < 0) + return false ; + dPar1 = m_dPar1 ; + dPar2 = m_dPar2 ; + return true ; +} + +//---------------------------------------------------------------------------- +// Calcola la distanza fra i due elemnti lineari, i punti di minima distanza e +// i loro rispettivi parametri. +// Se la coppia di punti di minima distanza non è unica ne viene scelta una +// in base a comodità di calcolo. +void +DistLineLine::Calculate( const Point3d& ptSt1, const Vector3d& vtD1, double dLen1, + const Point3d& ptSt2, const Vector3d& vtD2, double dLen2, + bool bIsSegment1, bool bIsSegment2) +{ + // Caso di elementi lineari paralleli/antiparalleli + if ( AreSameOrOppositeVectorExact( vtD1, vtD2)) { + // Almeno un elemento è una retta + if ( ! ( bIsSegment1 && bIsSegment2)) { + // Il primo elemento è segmento, quindi deve essere una retta il secondo + if ( bIsSegment1) { + Vector3d vtStSt = ptSt1 - ptSt2 ; + double dLong = vtStSt * vtD2 ; + Vector3d vtDist = vtStSt - dLong * vtD2 ; + m_dSqDist = vtDist.SqLen() ; + m_dDist = sqrt( m_dSqDist) ; + m_dPar1 = 0 ; + m_dPar2 = dLong ; + m_ptMinDist1 = ptSt1 ; + m_ptMinDist2 = ptSt2 + dLong * vtD2 ; + } + // Il primo elemento è una retta + else { + Vector3d vtStSt = ptSt2 - ptSt1 ; + double dLong = vtStSt * vtD1 ; + Vector3d vtDist = vtStSt - dLong * vtD1 ; + m_dSqDist = vtDist.SqLen() ; + m_dDist = sqrt( m_dSqDist) ; + m_dPar1 = dLong ; + m_dPar2 = 0 ; + m_ptMinDist1 = ptSt1 + dLong * vtD1 ; + m_ptMinDist2 = ptSt2 ; + } + } + // Entrambi gli elementi sono segmenti + else { + Point3d ptEn1 = ptSt1 + dLen1 * vtD1 ; + Point3d ptEn2 = ptSt2 + dLen2 * vtD2 ; + Vector3d vtStSt = ptSt2 - ptSt1 ; + Vector3d vtStEn = ptEn2 - ptSt1 ; + double dStU = vtStSt * vtD1 ; + double dEnU = vtStEn * vtD1 ; + // Classifico i punti del segmento segmento in base alla loro + // coordinata rispetto all'ordinamento generato dal primo. + double dMinPar, dMaxPar ; + Point3d ptMinPar, ptMaxPar ; + if ( dStU < dEnU) { + ptMinPar = ptSt2 ; + ptMaxPar = ptEn2 ; + dMinPar = dStU ; + dMaxPar = dEnU ; + } + else { + ptMinPar = ptEn2 ; + ptMaxPar = ptSt2 ; + dMinPar = dEnU ; + dMaxPar = dStU ; + } + // Possibili posizioni reciproche dei segmenti + if ( dMinPar > dLen1) { + m_dSqDist = SqDist( ptEn1, ptMinPar) ; + m_dDist = sqrt( m_dSqDist) ; + m_ptMinDist1 = ptEn1 ; + m_ptMinDist2 = ptMinPar ; + m_dPar1 = dLen1 ; + m_dPar2 = Clamp( ( m_ptMinDist2 - ptSt2) * vtD2, 0., dLen2) ; + } + else if ( dMinPar > 0) { + m_dSqDist = std::max( vtStSt * vtStSt - dStU * dStU, 0.) ; + m_dDist = sqrt( m_dSqDist) ; + m_ptMinDist1 = ptSt1 + dMinPar * vtD1 ; + m_ptMinDist2 = ptMinPar ; + m_dPar1 = dMinPar ; + m_dPar2 = Clamp( ( m_ptMinDist2 - ptSt2) * vtD2, 0., dLen2) ; + } + else if ( dMaxPar > 0) { + m_dSqDist = std::max( vtStSt * vtStSt - dStU * dStU, 0.) ; + m_dDist = sqrt( m_dSqDist) ; + m_ptMinDist1 = ptSt1 ; + m_ptMinDist2 = ptSt2 + ( ptSt1 - ptSt2) * vtD2 * vtD2 ; + m_dPar1 = 0 ; + m_dPar2 = Clamp( ( m_ptMinDist2 - ptSt2) * vtD2, 0., dLen2) ; + } + else { + m_dSqDist = SqDist( ptSt1, ptMaxPar) ; + m_dDist = sqrt( m_dSqDist) ; + m_ptMinDist1 = ptSt1 ; + m_ptMinDist2 = ptMaxPar ; + m_dPar1 = 0 ; + m_dPar2 = Clamp( ( m_ptMinDist2 - ptSt2) * vtD2, 0., dLen2) ; + } + } + return ; + } + // Caso generale + Vector3d vtDist0 = ptSt2 - ptSt1 ; + double dDist01 = vtDist0 * vtD1 ; + double dDist02 = vtDist0 * vtD2 ; + double dDotD1D2 = vtD1 * vtD2 ; + double dT1 = dDist01 + ( ( dDist01 * dDotD1D2 - dDist02) * dDotD1D2) / ( 1 - dDotD1D2 * dDotD1D2) ; + double dT2 = ( dDist01 * dDotD1D2 - dDist02) / ( 1 - dDotD1D2 * dDotD1D2) ; + double dMin1 = - INFINITO ; + double dMax1 = INFINITO ; + double dMin2 = - INFINITO ; + double dMax2 = INFINITO ; + double dSt1On2 = ( ptSt1 - ptSt2) * vtD2 ; + double dEn1On2 = ( ptSt1 + dLen1 * vtD1 - ptSt2) * vtD2 ; + if ( bIsSegment1) { + dMin1 = 0 ; + dMax1 = dLen1 ; + dMin2 = std::min( dSt1On2, dEn1On2) ; + dMax2 = std::max( dSt1On2, dEn1On2) ; + } + if ( bIsSegment2) { + double dSt2On1 = ( ptSt2 - ptSt1) * vtD1 ; + double dEn2On1 = ( ptSt2 + dLen2 * vtD2 - ptSt1) * vtD1 ; + dMin1 = std::max( dMin1, std::min( dSt2On1, dEn2On1)) ; + dMax1 = std::min( dMax1, std::max( dSt2On1, dEn2On1)) ; + dMin2 = std::max( dMin2, 0.) ; + dMax2 = std::min( dMax2, dLen2) ; + } + m_dPar1 = Clamp( dT1, dMin1, dMax1) ; + m_dPar2 = Clamp( dT2, dMin2, dMax2) ; + m_ptMinDist1 = ptSt1 + m_dPar1 * vtD1 ; + m_ptMinDist2 = ptSt2 + m_dPar2 * vtD2 ; + m_dSqDist = SqDist( m_ptMinDist1, m_ptMinDist2) ; + m_dDist = sqrt( m_dSqDist) ; +} \ No newline at end of file diff --git a/DistLineLine.h b/DistLineLine.h new file mode 100644 index 0000000..8d14a9d --- /dev/null +++ b/DistLineLine.h @@ -0,0 +1,57 @@ +// File : DistLineLine.h Data : Versione : +// Contenuto : Dichiarazione della classe distanza fra elementi lineari. +// +// +// +// Modifiche : +// +// +//---------------------------------------------------------------------------- + +#pragma once + +#include "/EgtDev/Include/EGkVector3d.h" +#include "/EgtDev/Include/EGkPoint3d.h" + + +class DistLineLine +{ + public : + DistLineLine( const Point3d& ptSt1, const Point3d& ptEn1, + const Point3d& ptSt2, const Point3d& ptEn2, + bool bIsSegment1 = true, bool bIsSegment2 = true) ; + DistLineLine( const Point3d& ptSt1, const Vector3d& vtD1, double dLen1, + const Point3d& ptSt2, const Vector3d& vtD2, double dLen2, + bool bIsSegment1 = true, bool bIsSegment2 = true) ; + + public : + bool GetSqDist( double& dSqDist) ; + bool GetDist( double& dDist) ; + bool IsEpsilon( double dTol) + { + double dSqDist ; + return ( GetSqDist( dSqDist) && ( dSqDist < SQ_EPS_ZERO || dSqDist < dTol * dTol)) ; + } + bool IsSmall( void) + { + return IsEpsilon( EPS_SMALL) ; + } + bool IsZero( void) + { + return IsEpsilon( EPS_ZERO) ; + } + bool GetMinDistPoints( Point3d& ptMinDist1, Point3d& ptMinDist2) ; + bool GetParamsAtMinDistPoints( double& dPar1, double& dPar2) ; + + private : + void Calculate( const Point3d& ptSt1, const Vector3d& vtD1, double dLen1, + const Point3d& ptSt2, const Vector3d& vtD2, double dLen2, + bool bIsSegment1 = true, bool bIsSegment2 = true) ; + private: + double m_dSqDist ; + double m_dDist ; + double m_dPar1 ; + double m_dPar2 ; + Point3d m_ptMinDist1 ; + Point3d m_ptMinDist2 ; +} ; \ No newline at end of file diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index aacde69..cf92b40 100644 Binary files a/EgtGeomKernel.rc and b/EgtGeomKernel.rc differ diff --git a/EgtGeomKernel.vcxproj b/EgtGeomKernel.vcxproj index 467e55d..c959fba 100644 --- a/EgtGeomKernel.vcxproj +++ b/EgtGeomKernel.vcxproj @@ -286,6 +286,7 @@ copy $(TargetPath) \EgtProg\Dll64 + @@ -526,6 +527,7 @@ copy $(TargetPath) \EgtProg\Dll64 + diff --git a/EgtGeomKernel.vcxproj.filters b/EgtGeomKernel.vcxproj.filters index 24128ef..aa29634 100644 --- a/EgtGeomKernel.vcxproj.filters +++ b/EgtGeomKernel.vcxproj.filters @@ -420,6 +420,9 @@ File di origine\GeoInters + + File di origine\GeoDist + @@ -1010,6 +1013,9 @@ File di intestazione\Include + + File di intestazione + diff --git a/IntersLineTria.cpp b/IntersLineTria.cpp index 3df4db5..c902512 100644 --- a/IntersLineTria.cpp +++ b/IntersLineTria.cpp @@ -17,6 +17,8 @@ #include "CurveLine.h" #include "IntersLineLine.h" #include "IntersLineTria.h" +#include "DistPointLine.h" +#include "DistLineLine.h" #include "/EgtDev/Include/EGkIntersLinePlane.h" #include "/EgtDev/Include/EGkFrame3d.h" #include @@ -53,15 +55,59 @@ IntersLineTria( const Point3d& ptL, const Vector3d& vtL, double dLen, const Tria // se c'è una intersezione else if ( nRes == ILPT_START || nRes == ILPT_END || nRes == ILPT_YES) { int nPTT = PointInTria( ptInt, trTria) ; - switch ( nPTT) { - case PTT_OUT : - return ILTT_NO ; - case PTT_VERT : - return ILTT_VERT ; - case PTT_EDGE : - return ILTT_EDGE ; - default : + if ( nPTT == PTT_IN) return ILTT_IN ; + else if ( nPTT == PTT_OUT) + return ILTT_NO ; + else if ( nPTT == PTT_VERT) { + for ( int nV = 0 ; nV < 3 ; ++ nV) { + int nW = ( nV + 1) % 3 ; + int nX = ( nW + 1) % 3 ; + Point3d ptPrev = trTria.GetP( nV) ; + Point3d ptCent = trTria.GetP( nW) ; + Point3d ptNext = trTria.GetP( nX) ; + Vector3d vtPrevCent = ptCent - ptPrev ; + vtPrevCent.Normalize() ; + Vector3d vtCentNext = ptNext - ptCent ; + vtCentNext.Normalize() ; + Vector3d vtPrevP = ptInt - ptPrev ; + vtPrevP.Normalize() ; + Vector3d vtCentP = ptInt - ptCent ; + vtCentP.Normalize() ; + if ( ( vtPrevCent ^ vtPrevP) * trTria.GetN() < 0 || ( vtCentNext ^ vtCentP) * trTria.GetN() < 0) { + DistPointLine DistCalc( ptCent, ptL, vtL, dLen, bFinite) ; + double dSqDist ; + DistCalc.GetSqDist( dSqDist) ; + if ( dSqDist < EPS_SMALL * EPS_SMALL) { + DistCalc.GetMinDistPoint( ptInt) ; + break ; + } + } + } + return ILTT_VERT ; + } + else { + for ( int nV = 0 ; nV < 3 ; ++ nV) { + int nW = ( nV + 1) % 3 ; + Point3d ptVertSt = trTria.GetP( nV) ; + Point3d ptVertEn = trTria.GetP( nW) ; + Vector3d vtSeg = ptVertEn - ptVertSt ; + double dSegLen = vtSeg.Len() ; + vtSeg /= dSegLen ; + Vector3d vtIntP = ptInt - ptVertSt ; + vtIntP.Normalize() ; + if ( ( vtSeg ^ vtIntP) * trTria.GetN() < 0) { + DistLineLine DistCalc( ptL, vtL, dLen, ptVertSt, vtSeg, dSegLen, bFinite, true) ; + double dSqDist ; + DistCalc.GetSqDist( dSqDist) ; + if ( dSqDist < EPS_SMALL * EPS_SMALL) { + Point3d ptMinDistOnSeg ; + DistCalc.GetMinDistPoints( ptInt, ptMinDistOnSeg) ; + break ; + } + } + } + return ILTT_EDGE ; } } // se la linea giace nel piano del triangolo diff --git a/VolZmapCreation.cpp b/VolZmapCreation.cpp index 1941259..597b30d 100644 --- a/VolZmapCreation.cpp +++ b/VolZmapCreation.cpp @@ -579,6 +579,22 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex // Determino le intersezioni della retta con la TriMesh ILSIVECTOR IntersectionResults ; intPLSTM.GetInters( ptP0, vtLen.v[(g+2)%3], IntersectionResults) ; + + for ( int nI = 0 ; nI < int( IntersectionResults.size()) - 3 ; ++ nI) { + int nJ = nI + 1 ; + int nK = nJ + 1 ; + int nT = nK + 1 ; + int nSgnI = IntersectionResults[nI].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nI].dCosDN > -EPS_SMALL ? 0 : - 1 ; + int nSgnJ = IntersectionResults[nJ].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nJ].dCosDN > -EPS_SMALL ? 0 : - 1 ; + int nSgnK = IntersectionResults[nK].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nK].dCosDN > -EPS_SMALL ? 0 : - 1 ; + int nSgnT = IntersectionResults[nT].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nT].dCosDN > -EPS_SMALL ? 0 : - 1 ; + double dUJ = IntersectionResults[nJ].dU ; + double dUK = IntersectionResults[nK].dU ; + if ( nSgnI != 0 && nSgnI == nSgnJ && nSgnK != 0 && nSgnK == nSgnT && nSgnI == - nSgnT && abs( dUJ - dUK) < EPS_SMALL) { + IntersectionResults.erase( IntersectionResults.begin() + nK) ; + IntersectionResults.erase( IntersectionResults.begin() + nJ) ; + } + } int nInt = int( IntersectionResults.size()) ;