diff --git a/CDeClosedSurfTmClosedSurfTm.cpp b/CDeClosedSurfTmClosedSurfTm.cpp
index 81bdc5a..50c7158 100644
--- a/CDeClosedSurfTmClosedSurfTm.cpp
+++ b/CDeClosedSurfTmClosedSurfTm.cpp
@@ -103,7 +103,7 @@ CDeClosedSurfTmClosedSurfTm( const SurfTriMesh& SurfA, const SurfTriMesh& SurfB,
// Processo il triangolo: se i due triangoli collidono ho finito.
if ( CDeTriaTria( trTriaA, trTriaB))
return true ;
- // Segno il triangolo come processato: nTFlag = 1
+ // Segno il triangolo come processato: nTemp = 1
SurfB.SetTempInt( nTB, 1) ;
}
}
diff --git a/EgtGeomKernel.vcxproj b/EgtGeomKernel.vcxproj
index f1ddf40..07e3f01 100644
--- a/EgtGeomKernel.vcxproj
+++ b/EgtGeomKernel.vcxproj
@@ -331,6 +331,8 @@ copy $(TargetPath) \EgtProg\Dll64
+
+
@@ -599,6 +601,8 @@ copy $(TargetPath) \EgtProg\Dll64
+
+
diff --git a/EgtGeomKernel.vcxproj.filters b/EgtGeomKernel.vcxproj.filters
index f6a5e4c..6e688e2 100644
--- a/EgtGeomKernel.vcxproj.filters
+++ b/EgtGeomKernel.vcxproj.filters
@@ -471,6 +471,12 @@
File di origine\GeoOffset
+
+ File di origine\GeoInters
+
+
+ File di origine\GeoInters
+
@@ -1103,6 +1109,12 @@
File di intestazione
+
+ File di intestazione
+
+
+ File di intestazione
+
diff --git a/IntersLineCone.cpp b/IntersLineCone.cpp
new file mode 100644
index 0000000..97f5d61
--- /dev/null
+++ b/IntersLineCone.cpp
@@ -0,0 +1,138 @@
+//----------------------------------------------------------------------------
+// EgalTech 2023-2023
+//----------------------------------------------------------------------------
+// File : IntersLineCone.cpp Data : 16.05.23 Versione : 2.5e3
+// Contenuto : Implementazione della intersezione linea/tronco di cono.
+//
+//
+//
+// Modifiche : 16.05.23 DS Creazione modulo.
+//
+//
+//----------------------------------------------------------------------------
+
+//--------------------------- Include ----------------------------------------
+#include "stdafx.h"
+#include "IntersLineCone.h"
+#include "IntersLineCyl.h"
+#include "/EgtDev/Include/ENkPolynomialRoots.h"
+
+using namespace std ;
+
+//----------------------------------------------------------------------------
+// Linea e tronco di cono sono nel medesimo riferimento.
+// Il tronco di cono è centrato sull'asse Z e appoggiato con RMin sul piano XY.
+// In caso di intersezione viene restituito true e i parametri in dU1 e dU2.
+//----------------------------------------------------------------------------
+bool
+IntersLineCone( const Point3d& ptL, const Vector3d& vtL,
+ double dRadMin, double dRadMax, double dHeight,
+ double& dU1, double& dU2)
+{
+ // Verifico il versore
+ if ( vtL.IsSmall())
+ return false ;
+
+ // Verifico il tronco di cono
+ if ( ( dRadMin < EPS_SMALL && dRadMax < EPS_SMALL) || dHeight < EPS_SMALL)
+ return false ;
+
+ // Se è un cilindro, rimando a questo
+ if ( abs( dRadMax - dRadMin) < EPS_SMALL)
+ return IntersLineCyl( ptL, vtL, ( dRadMin + dRadMax) / 2, dHeight, dU1, dU2) ;
+
+ // Se raggi invertiti, li scambio
+ if ( dRadMin > dRadMax)
+ swap( dRadMin, dRadMax) ;
+
+ // Tangente dell'angolo di semi-apertura del cono
+ double dTanTheta = ( dRadMax - dRadMin) / dHeight ;
+ double dSqTanTheta = dTanTheta * dTanTheta ;
+ double dSqCosTheta = 1 / ( 1 + dSqTanTheta) ;
+
+ // Determino le eventuali intersezioni con le due basi a quota minima e massima (solo se linea non giace sul cono)
+ int nBasInt = 0 ;
+ if ( abs( vtL.z) > EPS_ZERO) {
+ // le linee tangenti al cono non sono considerate intersecanti
+ bool bSameHAng = ( abs( abs( vtL.x) - abs( vtL.y)) < EPS_SMALL && abs( dSqCosTheta - vtL.z * vtL.z) < 2 * abs( vtL.z) * EPS_SMALL) ;
+ double EpsRad = ( bSameHAng ? - EPS_SMALL : EPS_SMALL) ;
+ Point3d ptInt1 = ptL + ( ( 0 - ptL.z) / vtL.z) * vtL ;
+ if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y < dRadMin * dRadMin + 2 * dRadMin * EpsRad) {
+ dU1 = ( ptInt1 - ptL) * vtL ;
+ nBasInt += 1 ;
+ }
+ Point3d ptInt2 = ptL + ( ( dHeight - ptL.z) / vtL.z) * vtL ;
+ if ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y < dRadMax * dRadMax + 2 * dRadMax * EpsRad) {
+ dU2 = ( ptInt2 - ptL) * vtL ;
+ nBasInt += 2 ;
+ }
+ }
+
+ // Se la linea interseca entrambe le basi, si sono trovate le due intersezioni
+ if ( nBasInt == 3) {
+ if ( dU1 > dU2)
+ swap( dU1, dU2) ;
+ // Trovate intersezioni
+ return true ;
+ }
+
+ // Posizione del vertice del cono
+ double dDeltaH = ( dRadMin < EPS_SMALL ? 0 : dRadMin / dTanTheta) ;
+ // Sposto il punto di passaggio della linea di conseguenza
+ Point3d ptMyL = ptL + Z_AX * dDeltaH ;
+
+ // Determino le intersezioni con la superficie laterale del cono
+ DBLVECTOR vdCoeff{ ptMyL.x * ptMyL.x + ptMyL.y * ptMyL.y - ptMyL.z * ptMyL.z * dSqTanTheta,
+ 2 * ( ptMyL.x * vtL.x + ptMyL.y * vtL.y - ptMyL.z * vtL.z * dSqTanTheta),
+ vtL.x * vtL.x + vtL.y * vtL.y - vtL.z * vtL.z * dSqTanTheta} ;
+ DBLVECTOR vdRoots ;
+ int nRoot = PolynomialRoots( 2, vdCoeff, vdRoots) ;
+
+ // Elimino le soluzioni cha danno intersezioni fuori dai limiti in Z del tronco di cono
+ if ( nRoot == 2) {
+ double dIntZ2 = ptL.z + vdRoots[1] * vtL.z ;
+ if ( dIntZ2 < 0 - EPS_SMALL || dIntZ2 > dHeight + EPS_SMALL)
+ -- nRoot ;
+ }
+ if ( nRoot >= 1) {
+ double dIntZ1 = ptL.z + vdRoots[0] * vtL.z ;
+ if ( dIntZ1 < 0 - EPS_SMALL || dIntZ1 > dHeight + EPS_SMALL) {
+ if ( nRoot == 2)
+ vdRoots[0] = vdRoots[1] ;
+ -- nRoot ;
+ }
+ }
+
+ // Due soluzioni: la retta interseca due volte la superficie laterale
+ if ( nRoot == 2) {
+ dU1 = vdRoots[0] ;
+ dU2 = vdRoots[1] ;
+ if ( dU1 > dU2)
+ swap( dU1, dU2) ;
+ // Trovate intersezioni
+ return true ;
+ }
+
+ // Una soluzione : la retta interseca la superficie laterale e un piano
+ else if ( nRoot == 1) {
+ // Se piano superiore
+ if ( nBasInt == 2) {
+ dU1 = vdRoots[0] ;
+ }
+ // altrimenti piano inferiore
+ else if ( nBasInt == 1) {
+ dU2 = vdRoots[0] ;
+ }
+ // altrimenti niente
+ else
+ return false ;
+ if ( dU1 > dU2)
+ swap( dU1, dU2) ;
+ // Trovate intersezioni
+ return true ;
+ }
+
+ // Nessuna soluzione : nessuna intersezione
+ else
+ return false ;
+}
diff --git a/IntersLineCone.h b/IntersLineCone.h
new file mode 100644
index 0000000..8ffa732
--- /dev/null
+++ b/IntersLineCone.h
@@ -0,0 +1,35 @@
+//----------------------------------------------------------------------------
+// EgalTech 2023-2023
+//----------------------------------------------------------------------------
+// File : IntersLineCone.h Data : 16.05.23 Versione : 2.5e3
+// Contenuto : Dichiarazione funzioni base per intersezione linea/cono tronco.
+//
+//
+//
+// Modifiche : 16.05.23 DS Creazione modulo.
+//
+//
+//----------------------------------------------------------------------------
+
+#pragma once
+
+#include "/EgtDev/Include/EGkPoint3d.h"
+
+//----------------------------------------------------------------------------
+// Linea e tronco di cono sono nel medesimo riferimento.
+// Il tronco di cono è centrato sull'asse Z e appoggiato con RMin sul piano XY.
+// Con intersezione viene restituito true e i parametri in dU1 e dU2.
+//----------------------------------------------------------------------------
+bool
+IntersLineCone( const Point3d& ptL, const Vector3d& vtL,
+ double dRadMin, double dRadMax, double dHeight,
+ double& dU1, double& dU2) ;
+
+//----------------------------------------------------------------------------
+inline bool
+TestIntersLineCone( const Point3d& ptL, const Vector3d& vtL,
+ double dRadMin, double dRadMax, double dHeight)
+{
+ double dU1, dU2 ;
+ return IntersLineCone( ptL, vtL, dRadMin, dRadMax, dHeight, dU1, dU2) ;
+}
diff --git a/IntersLineCyl.cpp b/IntersLineCyl.cpp
new file mode 100644
index 0000000..c78b3b7
--- /dev/null
+++ b/IntersLineCyl.cpp
@@ -0,0 +1,118 @@
+//----------------------------------------------------------------------------
+// EgalTech 2023-2023
+//----------------------------------------------------------------------------
+// File : IntersLineCyl.cpp Data : 16.05.23 Versione : 2.5e3
+// Contenuto : Implementazione della intersezione linea/cilindro.
+//
+//
+//
+// Modifiche : 16.05.23 DS Creazione modulo.
+//
+//
+//----------------------------------------------------------------------------
+
+//--------------------------- Include ----------------------------------------
+#include "stdafx.h"
+#include "IntersLineCyl.h"
+#include "/EgtDev/Include/ENkPolynomialRoots.h"
+
+using namespace std ;
+
+//----------------------------------------------------------------------------
+// Linea e cilindro sono nel medesimo riferimento.
+// Il cilindro è centrato sull'asse Z e appoggiato sul piano XY.
+// In caso di intersezione viene restituito true e i parametri in dU1 e dU2.
+//----------------------------------------------------------------------------
+bool
+IntersLineCyl( const Point3d& ptL, const Vector3d& vtL,
+ double dRad, double dHeight,
+ double& dU1, double& dU2)
+{
+ // Verifico il versore
+ if ( vtL.IsSmall())
+ return false ;
+
+ // Verifico il cilindro
+ if ( dRad < EPS_SMALL || dHeight < EPS_SMALL)
+ return false ;
+
+ // Determino le eventuali intersezioni con le due basi a quota minima e massima (solo se linea non parallela ad esse)
+ int nBasInt = 0 ;
+ if ( abs( vtL.z) > EPS_ZERO) {
+ // le linee tangenti al cilindro non sono considerate intersecanti
+ double EpsRad = ( vtL.IsZeroXY() ? - EPS_SMALL : EPS_SMALL) ;
+ Point3d ptInt1 = ptL + ( ( 0 - ptL.z) / vtL.z) * vtL ;
+ if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y < dRad * dRad + 2 * dRad * EpsRad) {
+ dU1 = ( ptInt1 - ptL) * vtL ;
+ nBasInt += 1 ;
+ }
+ Point3d ptInt2 = ptL + ( ( dHeight - ptL.z) / vtL.z) * vtL ;
+ if ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y < dRad * dRad + 2 * dRad * EpsRad) {
+ dU2 = ( ptInt2 - ptL) * vtL ;
+ nBasInt += 2 ;
+ }
+ }
+
+ // Se la linea interseca entrambe le basi, si sono trovate le due intersezioni
+ if ( nBasInt == 3) {
+ if ( dU1 > dU2)
+ swap( dU1, dU2) ;
+ // Trovate intersezioni
+ return true ;
+ }
+
+ // Determino le intersezioni con la superficie laterale del cilindro
+ DBLVECTOR vdCoeff{ ptL.x * ptL.x + ptL.y * ptL.y - dRad * dRad,
+ 2 * ( ptL.x * vtL.x + ptL.y * vtL.y),
+ vtL.x * vtL.x + vtL.y * vtL.y} ;
+ DBLVECTOR vdRoots ;
+ int nRoot = PolynomialRoots( 2, vdCoeff, vdRoots) ;
+
+ // Elimino le soluzioni cha danno intersezioni fuori dai limiti in Z del cilindro
+ if ( nRoot == 2) {
+ double dIntZ2 = ptL.z + vdRoots[1] * vtL.z ;
+ if ( dIntZ2 < 0 - EPS_SMALL || dIntZ2 > dHeight + EPS_SMALL)
+ -- nRoot ;
+ }
+ if ( nRoot >= 1) {
+ double dIntZ1 = ptL.z + vdRoots[0] * vtL.z ;
+ if ( dIntZ1 < 0 - EPS_SMALL || dIntZ1 > dHeight + EPS_SMALL) {
+ if ( nRoot == 2)
+ vdRoots[0] = vdRoots[1] ;
+ -- nRoot ;
+ }
+ }
+
+ // Due soluzioni: la retta interseca due volte la superficie laterale
+ if ( nRoot == 2) {
+ dU1 = vdRoots[0] ;
+ dU2 = vdRoots[1] ;
+ if ( dU1 > dU2)
+ swap( dU1, dU2) ;
+ // Trovate intersezioni
+ return true ;
+ }
+
+ // Una soluzione : la retta interseca la superficie laterale e un piano
+ else if ( nRoot == 1) {
+ // Se piano superiore
+ if ( nBasInt == 2) {
+ dU1 = vdRoots[0] ;
+ }
+ // altrimenti piano inferiore
+ else if ( nBasInt == 1) {
+ dU2 = vdRoots[0] ;
+ }
+ // altrimenti niente
+ else
+ return false ;
+ if ( dU1 > dU2)
+ swap( dU1, dU2) ;
+ // Trovate intersezioni
+ return true ;
+ }
+
+ // Nessuna soluzione : nessuna intersezione
+ else
+ return false ;
+}
diff --git a/IntersLineCyl.h b/IntersLineCyl.h
new file mode 100644
index 0000000..4448739
--- /dev/null
+++ b/IntersLineCyl.h
@@ -0,0 +1,35 @@
+//----------------------------------------------------------------------------
+// EgalTech 2023-2023
+//----------------------------------------------------------------------------
+// File : IntersLineCyl.h Data : 16.05.23 Versione : 2.5e3
+// Contenuto : Dichiarazione funzioni base per intersezione linea/cilindro.
+//
+//
+//
+// Modifiche : 16.05.23 DS Creazione modulo.
+//
+//
+//----------------------------------------------------------------------------
+
+#pragma once
+
+#include "/EgtDev/Include/EGkPoint3d.h"
+
+//----------------------------------------------------------------------------
+// Linea e cilindro sono nel medesimo riferimento.
+// Il cilindro è centrato sull'asse Z e appoggiato sul piano XY.
+// Con intersezione viene restituito true e i parametri in dU1 e dU2.
+//----------------------------------------------------------------------------
+bool
+IntersLineCyl( const Point3d& ptL, const Vector3d& vtL,
+ double dRad, double dHeight,
+ double& dU1, double& dU2) ;
+
+//----------------------------------------------------------------------------
+inline bool
+TestIntersLineCyl( const Point3d& ptL, const Vector3d& vtL,
+ double dRad, double dHeight)
+{
+ double dU1, dU2 ;
+ return IntersLineCyl( ptL, vtL, dRad, dHeight, dU1, dU2) ;
+}
diff --git a/IntersLineSphere.cpp b/IntersLineSphere.cpp
index befeab1..6d2629b 100644
--- a/IntersLineSphere.cpp
+++ b/IntersLineSphere.cpp
@@ -26,20 +26,25 @@ IntersLineSphere( const Point3d& ptL, const Vector3d& vtL, const Point3d& ptCen,
return ILST_NO ;
// Proiezione del centro della sfera sulla linea
Point3d ptP = ptL + (( ptCen - ptL) * vtL) * vtL ;
- // Distanza di questo punto di proiezione dal centro della sfera
- double dDist = Dist( ptCen, ptP) ;
+ // Quadrato della distanza di questo punto di proiezione dal centro della sfera
+ double dSqDist = SqDist( ptCen, ptP) ;
+ // Differenza tra quadrato del raggio e quadrato della distanza
+ double dSqDelta = dRad * dRad - dSqDist ;
+
+ // Se distanza superiore al raggio, nessuna intersezione
+ if ( dSqDelta < - 2 * dRad * EPS_SMALL)
+ return ILST_NO ;
+
// Se distanza uguale al raggio, intersezione tangente
- if ( abs( dDist - dRad) < EPS_SMALL) {
+ if ( dSqDelta < EPS_SMALL * EPS_SMALL) {
ptI1 = ptP ;
ptI2 = ptP ;
return ILST_TG ;
}
- // Se distanza superiore al raggio, nessuna intersezione
- if ( dDist > dRad)
- return ILST_NO ;
+
// Distanza inferiore al raggio, due intersezioni secanti
- double dDist2 = sqrt( dRad * dRad - dDist * dDist) ;
+ double dDist2 = sqrt( dSqDelta) ;
ptI1 = ptP - dDist2 * vtL ;
ptI2 = ptP + dDist2 * vtL ;
return ILST_SEC ;
-}
\ No newline at end of file
+}
diff --git a/VolZmapCalculus.cpp b/VolZmapCalculus.cpp
index da84395..1443894 100644
--- a/VolZmapCalculus.cpp
+++ b/VolZmapCalculus.cpp
@@ -17,6 +17,8 @@
#include "VolZmap.h"
#include "GeoConst.h"
#include "IntersLineBox.h"
+#include "IntersLineCyl.h"
+#include "IntersLineCone.h"
#include "IntersLineSurfStd.h"
#include "/EgtDev/Include/EGkIntersLineTria.h"
#include "/EgtDev/Include/EGkIntersLinePlane.h"
@@ -453,16 +455,13 @@ VolZmap::AvoidSimpleBox( const Frame3d& frBox, const Vector3d& vtDiag, bool bPre
continue ;
for ( int k = 0 ; k < 5 ; ++ k) {
Point3d ptT = ptO + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ;
- if ( k == 0)
- ;
- else if ( k == 1)
- ptT += - 0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ;
- else if ( k == 2)
- ptT += + 0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ;
- else if ( k == 3)
- ptT += + 0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ;
- else if ( k == 4)
- ptT += - 0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ;
+ switch ( k) {
+ case 0 : break ;
+ case 1 : ptT += -0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ; break ;
+ case 2 : ptT += +0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ; break ;
+ case 3 : ptT += +0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ; break ;
+ case 4 : ptT += -0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ; break ;
+ }
double dZmin, dZmax ;
if ( IntersLineBox( ptT, vtK, ORIG, ORIG + vtDiag, dZmin, dZmax)) {
for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
@@ -829,6 +828,14 @@ VolZmap::AvoidSimpleCylinder( const Frame3d& frCyl, double dR, double dH, bool b
int nStJ = Clamp( int( b3Int.GetMin().y / m_dStep), 0, m_nNy[0] - 1) ;
int nEnJ = Clamp( int( b3Int.GetMax().y / m_dStep), 0, m_nNy[0] - 1) ;
+ // Vettore direzione dei dexel nel riferimento del Box
+ Vector3d vtK = GetToLoc( Z_AX, frCylInt) ;
+
+ // Riferimento intrinseco dei dexel nel riferimento del box
+ Point3d ptO = GetToLoc( ORIG, frCylInt) ;
+ Vector3d vtX = GetToLoc( X_AX, frCylInt) ;
+ Vector3d vtY = GetToLoc( Y_AX, frCylInt) ;
+
// Ciclo di intersezione dei dexel con il cilindro (nel riferimento intrinseco)
for ( int i = nStI ; i <= nEnI ; ++ i) {
for ( int j = nStJ ; j <= nEnJ ; ++ j) {
@@ -839,19 +846,16 @@ VolZmap::AvoidSimpleCylinder( const Frame3d& frCyl, double dR, double dH, bool b
if ( m_Values[0][nPos][nSize-1].dMax < b3Int.GetMin().z || m_Values[0][nPos][0].dMin > b3Int.GetMax().z)
continue ;
for ( int k = 0 ; k < 5 ; ++ k) {
- Point3d ptT = ORIG + ( i + 0.5) * m_dStep * X_AX + ( j + 0.5) * m_dStep * Y_AX ;
+ Point3d ptT = ptO + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ;
switch ( k) {
case 0 : break ;
- case 1 : ptT += -0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
- case 2 : ptT += +0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
- case 3 : ptT += +0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
- case 4 : ptT += -0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
+ case 1 : ptT += -0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ; break ;
+ case 2 : ptT += +0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ; break ;
+ case 3 : ptT += +0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ; break ;
+ case 4 : ptT += -0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ; break ;
}
- Point3d ptI1, ptI2 ;
- Vector3d vtN1, vtN2 ;
- if ( IntersLineCylinder( ptT, Z_AX, frCylInt, dH, dR, true, true, ptI1, vtN1, ptI2, vtN2)) {
- double dZmin = min( ptI1.z, ptI2.z) ;
- double dZmax = max( ptI1.z, ptI2.z) ;
+ double dZmin, dZmax ;
+ if ( IntersLineCyl( ptT, vtK, dR, dH, dZmin, dZmax)) {
for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
if ( dZmax > m_Values[0][nPos][nIndex].dMin - EPS_SMALL &&
dZmin < m_Values[0][nPos][nIndex].dMax + EPS_SMALL)
@@ -867,6 +871,7 @@ VolZmap::AvoidSimpleCylinder( const Frame3d& frCyl, double dR, double dH, bool b
else {
// Ciclo di intersezione dei dexel con il cilindro (nel riferimento intrinseco)
for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) {
+ Point3d ptO = ORIG ;
Vector3d vtX = X_AX ;
Vector3d vtY = Y_AX ;
Vector3d vtK = Z_AX ;
@@ -891,6 +896,11 @@ VolZmap::AvoidSimpleCylinder( const Frame3d& frCyl, double dR, double dH, bool b
vtY = X_AX ;
vtK = Y_AX ;
}
+ // Passo da riferimento intrinseco griglia a riferimento cilindro
+ ptO.ToLoc( frCylInt) ;
+ vtX.ToLoc( frCylInt) ;
+ vtY.ToLoc( frCylInt) ;
+ vtK.ToLoc( frCylInt) ;
// Limiti su indici
int nStI = Clamp( int( ptBoxInf.x / m_dStep), 0, m_nNx[nMap] - 1) ;
int nEnI = Clamp( int( ptBoxSup.x / m_dStep), 0, m_nNx[nMap] - 1) ;
@@ -905,27 +915,13 @@ VolZmap::AvoidSimpleCylinder( const Frame3d& frCyl, double dR, double dH, bool b
continue ;
if ( m_Values[nMap][nPos][nSize-1].dMax < b3Int.GetMin().z || m_Values[nMap][nPos][0].dMin > b3Int.GetMax().z)
continue ;
- Point3d ptT = ORIG + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ;
- Point3d ptI1, ptI2 ;
- Vector3d vtN1, vtN2 ;
- // La linea del dexel interseca il cilindro.
- if ( IntersLineCylinder( ptT, vtK, frCylInt, dH, dR, true, true, ptI1, vtN1, ptI2, vtN2)) {
- double dMinU, dMaxU ;
- if ( nMap == 0) {
- dMinU = min( ptI1.z, ptI2.z) ;
- dMaxU = max( ptI1.z, ptI2.z) ;
- }
- else if ( nMap == 1) {
- dMinU = min( ptI1.x, ptI2.x) ;
- dMaxU = max( ptI1.x, ptI2.x) ;
- }
- else {
- dMinU = min( ptI1.y, ptI2.y) ;
- dMaxU = max( ptI1.y, ptI2.y) ;
- }
- // Ciclo sui segmenti del dexel.
+ Point3d ptT = ptO + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ;
+ double dMinU, dMaxU ;
+ // La retta associata al dexel interseca il cilindro.
+ if ( IntersLineCyl( ptT, vtK, dR, dH, dMinU, dMaxU)) {
+ // Ciclo sui segmenti del dexel
for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
- // Se il segmento è interno all'intervallo d'intersezione, ho finito.
+ // Se il segmento è interno all'intervallo d'intersezione, ho finito.
if ( dMaxU > m_Values[nMap][nPos][nIndex].dMin - EPS_SMALL &&
dMinU < m_Values[nMap][nPos][nIndex].dMax + EPS_SMALL)
return false ;
@@ -1126,177 +1122,103 @@ VolZmap::AvoidSimpleConeFrustum( const Frame3d& frCone, double dMinRad, double d
int nEnI = Clamp( int( b3Int.GetMax().x / m_dStep), 0, m_nNx[0] - 1) ;
int nStJ = Clamp( int( b3Int.GetMin().y / m_dStep), 0, m_nNy[0] - 1) ;
int nEnJ = Clamp( int( b3Int.GetMax().y / m_dStep), 0, m_nNy[0] - 1) ;
- // Limiti su Z
- double dZmin = b3Int.GetMin().z ;
- double dZmax = b3Int.GetMax().z ;
- // Numero massimo di thread
- int nThreadMax = max( 1, int( thread::hardware_concurrency()) - 1) ;
- // se un solo thread
- if ( nThreadMax == 1) {
- m_bBreak = false ;
- bool bCollision = SingleMapDexelConeCollision( nStI, nEnI, nStJ, nEnJ, ptRefPoint, vtRefAx,
- dMinRad, dMaxRad, dHeight, dZmin, dZmax) ;
- return ( ! bCollision) ;
- }
- // altrimenti esecuzione in parallelo dei calcoli
- else {
- //string sOut = "I=" + ToString( nStI) + "," + ToString( nEnI) + " J=" + ToString( nStJ) + "," + ToString( nEnJ) ;
- //LOG_INFO( GetEGkLogger(), sOut.c_str())
- m_bBreak = false ;
- // lancio dei thread
- int nSpanI = ( nEnI - nStI + 1) / nThreadMax + 1 ;
- int nSpanJ = ( nEnJ - nStJ + 1) / nThreadMax + 1 ;
- bool bOnI = ( nSpanI >= nSpanJ) ;
- int nThreadTot = 0 ;
- vector< future> vRes( nThreadMax) ;
- for ( int nT = 0 ; nT < nThreadMax ; ++ nT) {
- int nMyStI = ( bOnI ? nStI + nT * nSpanI : nStI) ;
- int nMyEnI = ( bOnI ? min( nMyStI + nSpanI - 1, nEnI) : nEnI) ;
- int nMyStJ = ( bOnI ? nStJ : nStJ + nT * nSpanJ) ;
- int nMyEnJ = ( bOnI ? nEnJ : min( nMyStJ + nSpanJ - 1, nEnJ)) ;
- if ( nMyStI > nEnI || nMyStJ > nEnJ)
- break ;
- //string sOut = "MyI=" + ToString( nMyStI) + "," + ToString( nMyEnI) + " MyJ=" + ToString( nMyStJ) + "," + ToString( nMyEnJ) ;
- //LOG_INFO( GetEGkLogger(), sOut.c_str())
- vRes[nT] = async( launch::async, &VolZmap::SingleMapDexelConeCollision, this,
- nMyStI, nMyEnI, nMyStJ, nMyEnJ, cref( ptRefPoint), cref( vtRefAx),
- dMinRad, dMaxRad, dHeight, dZmin, dZmax) ;
- ++ nThreadTot ;
- }
- // recupero i risultati dei thread alla loro terminazione
- bool bCollision = false ;
- int nTerminated = 0 ;
- while ( nTerminated < nThreadTot) {
- for ( int nT = 0 ; nT < nThreadTot ; ++ nT) {
- // Async terminato
- if ( vRes[nT].valid() && vRes[nT].wait_for( chrono::nanoseconds{ 1}) == future_status::ready) {
- ++ nTerminated ;
- // Se c'è collisione ...
- if ( vRes[nT].get()) {
- bCollision = true ;
- m_bBreak = true ;
+
+ // Vettore direzione dei dexel nel riferimento del tronco di cono
+ Vector3d vtK = GetToLoc( Z_AX, frConeInt) ;
+
+ // Riferimento intrinseco dei dexel nel riferimento del tronco di cono
+ Point3d ptO = GetToLoc( ORIG, frConeInt) ;
+ Vector3d vtX = GetToLoc( X_AX, frConeInt) ;
+ Vector3d vtY = GetToLoc( Y_AX, frConeInt) ;
+
+ // Ciclo di intersezione dei dexel con il tronco di cono (nel riferimento intrinseco)
+ for ( int i = nStI ; i <= nEnI ; ++ i) {
+ for ( int j = nStJ ; j <= nEnJ ; ++ j) {
+ int nPos = j * m_nNx[0] + i ;
+ int nSize = int( m_Values[0][nPos].size()) ;
+ if ( nSize == 0)
+ continue ;
+ if ( m_Values[0][nPos][nSize-1].dMax < b3Int.GetMin().z || m_Values[0][nPos][0].dMin > b3Int.GetMax().z)
+ continue ;
+ for ( int k = 0 ; k < 5 ; ++ k) {
+ Point3d ptT = ptO + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ;
+ switch ( k) {
+ case 0 : break ;
+ case 1 : ptT += -0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ; break ;
+ case 2 : ptT += +0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ; break ;
+ case 3 : ptT += +0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ; break ;
+ case 4 : ptT += -0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ; break ;
+ }
+ double dZmin, dZmax ;
+ if ( IntersLineCone( ptT, vtK, dMinRad, dMaxRad, dHeight, dZmin, dZmax)) {
+ for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
+ if ( dZmax > m_Values[0][nPos][nIndex].dMin - EPS_SMALL &&
+ dZmin < m_Values[0][nPos][nIndex].dMax + EPS_SMALL)
+ return false ;
}
}
}
}
- return ( ! bCollision) ;
}
}
// Uso tutte le mappe
else {
- // Ciclo sulle mappe.
+ // Ciclo di intersezione dei dexel con il cilindro (nel riferimento intrinseco)
for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) {
- Point3d ptInfIntBox = b3Int.GetMin() ;
- Point3d ptSupIntBox = b3Int.GetMax() ;
- // Dal sistema intrinseco al sistema griglia (per la prima griglia coincidono).
+ Point3d ptO = ORIG ;
+ Vector3d vtX = X_AX ;
+ Vector3d vtY = Y_AX ;
+ Vector3d vtK = Z_AX ;
+ // Estremi del box
+ Point3d ptBoxInf = b3Int.GetMin() ;
+ Point3d ptBoxSup = b3Int.GetMax() ;
if ( nMap == 1) {
- swap( ptInfIntBox.x, ptInfIntBox.z) ;
- swap( ptInfIntBox.x, ptInfIntBox.y) ;
- swap( ptSupIntBox.x, ptSupIntBox.z) ;
- swap( ptSupIntBox.x, ptSupIntBox.y) ;
+ swap( ptBoxInf.x, ptBoxInf.z) ;
+ swap( ptBoxInf.x, ptBoxInf.y) ;
+ swap( ptBoxSup.x, ptBoxSup.z) ;
+ swap( ptBoxSup.x, ptBoxSup.y) ;
+ vtX = Y_AX ;
+ vtY = Z_AX ;
+ vtK = X_AX ;
}
else if ( nMap == 2) {
- swap( ptInfIntBox.y, ptInfIntBox.z) ;
- swap( ptInfIntBox.x, ptInfIntBox.y) ;
- swap( ptSupIntBox.y, ptSupIntBox.z) ;
- swap( ptSupIntBox.x, ptSupIntBox.y) ;
+ swap( ptBoxInf.y, ptBoxInf.z) ;
+ swap( ptBoxInf.x, ptBoxInf.y) ;
+ swap( ptBoxSup.y, ptBoxSup.z) ;
+ swap( ptBoxSup.x, ptBoxSup.y) ;
+ vtX = Z_AX ;
+ vtY = X_AX ;
+ vtK = Y_AX ;
}
+ // Passo da riferimento intrinseco griglia a riferimento cilindro
+ ptO.ToLoc( frConeInt) ;
+ vtX.ToLoc( frConeInt) ;
+ vtY.ToLoc( frConeInt) ;
+ vtK.ToLoc( frConeInt) ;
// Limiti su indici
- int nStI = Clamp( int( ptInfIntBox.x / m_dStep), 0, m_nNx[nMap] - 1) ;
- int nEnI = Clamp( int( ptSupIntBox.x / m_dStep), 0, m_nNx[nMap] - 1) ;
- int nStJ = Clamp( int( ptInfIntBox.y / m_dStep), 0, m_nNy[nMap] - 1) ;
- int nEnJ = Clamp( int( ptSupIntBox.y / m_dStep), 0, m_nNy[nMap] - 1) ;
- // Ciclo sui dexel.
- for ( int nDex = 0 ; nDex < int( m_Values[nMap].size()) ; ++ nDex) {
- int nDexSize = (int)m_Values[nMap][nDex].size() ;
- if ( nDexSize == 0 ||
- m_Values[nMap][nDex][ nDexSize- 1].dMax < ptInfIntBox.z ||
- m_Values[nMap][nDex][0].dMin > ptSupIntBox.z)
- continue ;
- // Indici del dexel.
- int nI = nDex % m_nNx[nMap] ;
- int nJ = nDex / m_nNx[nMap] ;
- // Se fuori dalla regione ammissibile salto l'iterazione
- if ( nI < nStI || nI > nEnI || nJ < nStJ || nJ > nEnJ)
- continue ;
- // Posizione del dexel.
- double dX = ( nI + 0.5) * m_dStep ;
- double dY = ( nJ + 0.5) * m_dStep ;
- Point3d ptLineSt( dX, dY, 0.) ;
- Vector3d vtLineDir( 0., 0., 1.) ;
- // Dal sistema griglia al sistema intrinseco (per la prima griglia coincidono).
- if ( nMap == 1) {
- swap( ptLineSt.x, ptLineSt.y) ;
- swap( ptLineSt.x, ptLineSt.z) ;
- swap( vtLineDir.x, vtLineDir.y) ;
- swap( vtLineDir.x, vtLineDir.z) ;
- }
- else if ( nMap == 2) {
- swap( ptLineSt.x, ptLineSt.y) ;
- swap( ptLineSt.y, ptLineSt.z) ;
- swap( vtLineDir.x, vtLineDir.y) ;
- swap( vtLineDir.y, vtLineDir.z) ;
- }
- // Cono proprio
- if ( dMinRad < EPS_SMALL) {
- double dMinPar = m_Values[nMap][nDex][0].dMin ;
- double dMaxPar = m_Values[nMap][nDex][nDexSize - 1].dMax ;
- Point3d ptSegSt = ptLineSt + dMinPar * vtLineDir ;
- double dU1, dU2 ;
- int nIntType = SegmentCone( ptSegSt, vtLineDir, dMaxPar - dMinPar, ptRefPoint, vtRefAx,
- dMaxRad, dHeight, dU1, dU2) ;
- if ( nIntType == LinCompCCIntersType::CC_ERROR_INT)
- return false ;
- else if ( nIntType != LinCompCCIntersType::CC_NO_INTERS) {
- for ( int nIndex = 0 ; nIndex < nDexSize ; nIndex += 1) {
- if ( m_Values[nMap][nDex][nIndex].dMax >= dU1 && m_Values[nMap][nDex][nIndex].dMin <= dU2)
- return false ;
- }
- }
- nIntType = SegmentDisc( ptSegSt, vtLineDir, dMaxPar - dMinPar, ptRefPoint + dHeight * vtRefAx,
- vtRefAx, dMaxRad, dU1, dU2) ;
- if ( nIntType == LinCompDiscIntersType::D_ERROR_INT)
- return false ;
- else if ( nIntType != LinCompDiscIntersType::D_NO_INTERS) {
- for ( int nIndex = 0 ; nIndex < nDexSize ; nIndex += 1) {
- if ( m_Values[nMap][nDex][nIndex].dMax >= dU1 && m_Values[nMap][nDex][nIndex].dMin <= dU2)
- return false ;
- }
- }
- }
- // Tronco di cono
- else {
- double dMinPar = m_Values[nMap][nDex][0].dMin ;
- double dMaxPar = m_Values[nMap][nDex][nDexSize - 1].dMax ;
- Point3d ptSegSt = ptLineSt + dMinPar * vtLineDir ;
- double dU1, dU2 ;
- int nIntType = SegmentConeFrustum( ptSegSt, vtLineDir, dMaxPar - dMinPar, ptRefPoint, vtRefAx,
- dMinRad, dMaxRad, dHeight, dU1, dU2) ;
- if ( nIntType == LinCompCCIntersType::CC_ERROR_INT)
- return false ;
- else if ( nIntType != LinCompCCIntersType::CC_NO_INTERS) {
- for ( int nIndex = 0 ; nIndex < nDexSize ; nIndex += 1) {
- if ( m_Values[nMap][nDex][nIndex].dMax >= dU1 && m_Values[nMap][nDex][nIndex].dMin <= dU2)
- return false ;
- }
- }
- nIntType = SegmentDisc( ptSegSt, vtLineDir, dMaxPar - dMinPar, ptRefPoint + dHeight * vtRefAx,
- vtRefAx, dMaxRad, dU1, dU2) ;
- if ( nIntType == LinCompDiscIntersType::D_ERROR_INT)
- return false ;
- else if ( nIntType != LinCompDiscIntersType::D_NO_INTERS) {
- for ( int nIndex = 0 ; nIndex < nDexSize ; nIndex += 1) {
- if ( m_Values[nMap][nDex][nIndex].dMax >= dU1 && m_Values[nMap][nDex][nIndex].dMin <= dU2)
- return false ;
- }
- }
- nIntType = SegmentDisc( ptSegSt, vtLineDir, dMaxPar - dMinPar, ptRefPoint, vtRefAx, dMinRad, dU1, dU2) ;
- if ( nIntType == LinCompDiscIntersType::D_ERROR_INT)
- return false ;
- else if ( nIntType != LinCompDiscIntersType::D_NO_INTERS) {
- for ( int nIndex = 0 ; nIndex < nDexSize ; nIndex += 1) {
- if ( m_Values[nMap][nDex][nIndex].dMax >= dU1 && m_Values[nMap][nDex][nIndex].dMin <= dU2)
+ int nStI = Clamp( int( ptBoxInf.x / m_dStep), 0, m_nNx[nMap] - 1) ;
+ int nEnI = Clamp( int( ptBoxSup.x / m_dStep), 0, m_nNx[nMap] - 1) ;
+ int nStJ = Clamp( int( ptBoxInf.y / m_dStep), 0, m_nNy[nMap] - 1) ;
+ int nEnJ = Clamp( int( ptBoxSup.y / m_dStep), 0, m_nNy[nMap] - 1) ;
+ // Ciclo sui dexel.
+ for ( int i = nStI ; i <= nEnI ; ++ i) {
+ for ( int j = nStJ ; j <= nEnJ ; ++ j) {
+ int nPos = j * m_nNx[nMap] + i ;
+ int nSize = int( m_Values[nMap][nPos].size()) ;
+ if ( nSize == 0)
+ continue ;
+ if ( m_Values[nMap][nPos][nSize-1].dMax < b3Int.GetMin().z || m_Values[nMap][nPos][0].dMin > b3Int.GetMax().z)
+ continue ;
+ Point3d ptT = ptO + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ;
+ double dMinU, dMaxU ;
+ // La retta associata al dexel interseca il tronco di cono
+ if ( IntersLineCone( ptT, vtK, dMinRad, dMaxRad, dHeight, dMinU, dMaxU)) {
+ // Ciclo sui segmenti del dexel
+ for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
+ // Se il segmento è interno all'intervallo d'intersezione, ho finito.
+ if ( dMaxU > m_Values[nMap][nPos][nIndex].dMin - EPS_SMALL &&
+ dMinU < m_Values[nMap][nPos][nIndex].dMax + EPS_SMALL)
return false ;
}
}