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

EgtGeomKernel :

Browse Source 
- miglioramenti e ottimizzazioni in CD su Zmap per cilindri e tronchi di cono.
This commit is contained in:
DarioS
2023-05-16 20:26:35 +02:00
parent 8420edecb5
commit b78212c3a1
9 changed files with 474 additions and 205 deletions
Download Patch File Download Diff File Expand all files Collapse all files
VolZmapCalculus.cpp
+118 -196
View File
@@ -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<bool>> 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 ;
}
}