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

EgtGeomKernel :

Browse Source 
- migliorata GetDepth di VolZmap, con flag per scelta algoritmo.
This commit is contained in:
Dario Sassi
2017-10-23 06:47:22 +00:00
parent 9770d57793
commit 2ecd82c61f
13 changed files with 1306 additions and 440 deletions
Download Patch File Download Diff File Expand all files Collapse all files
VolTriZmapCalculus.cpp
+149 -113
View File
@@ -79,6 +79,19 @@ VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV,
return ( dU2 >= dU1) ;
}
//----------------------------------------------------------------------------
// La retta è rappresentata da un punto e dal versore espressi nel riferimento dello Zmap.
// Se non vi è intersezione, viene restituito falso
bool
VolZmap::IntersLineZMapLattice( const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) const
{
// Punti estremi della diagonale del box del reticolo
Point3d ptMin( - 0.5 * m_dStep, - 0.5 * m_dStep, - 0.5 * m_dStep) ;
Point3d ptMax( ( m_nNx[0] + 0.5) * m_dStep, ( m_nNy[0] + 0.5) * m_dStep, ( m_nNy[1] + 0.5) * m_dStep) ;
return IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2) ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IntersLineZMapBBox( const Point3d& ptP, const Vector3d& vtV, unsigned int nGrid, double& dU1, double& dU2) const
@@ -130,7 +143,7 @@ VolZmap::IntersLineDexel( const Point3d& ptP, const Vector3d& vtV, unsigned int
//----------------------------------------------------------------------------
bool
VolZmap::IntersRayDexel( const Point3d& ptP, const Vector3d& vtV, unsigned int nGrid, unsigned int nI, unsigned int nJ,
VolZmap::IntersRayDexel( const Point3d& ptP, const Vector3d& vtV, unsigned int nGrid, unsigned int nI, unsigned int nJ,
double& dU1, double& dU2) const
{
// Determino l'indice del dexel e il numero di suoi intervalli
@@ -172,25 +185,34 @@ VolZmap::IntersRayDexel( const Point3d& ptP, const Vector3d& vtV, unsigned int n
}
//----------------------------------------------------------------------------
// Calcola la profondità del materiale lungo un raggio, identificato da un punto e un versoreore.
// Punto e versore devono essere espressi nel sistema di riferimento locale, in cui è immerso quello dello Zmap.
// InLength = distanza di ingresso (se -1 il punto è interno, se -2 il punto è esterno e il raggio non interseca lo Zmap)
// OutLength = distanza di uscita
bool
VolZmap::GetDepth( const Point3d& ptPGlob, const Vector3d& vtDir, double& dInLength, double& dOutLength) const
VolZmap::GetDepth( const Point3d& ptPLoc, const Vector3d& vtDLoc, double& dInLength, double& dOutLength, bool bExact) const
{
#if 1
GetDepthWithVoxel( ptPGlob, vtDir, dInLength, dOutLength) ;
return true ;
#else
Point3d ptP = ptPLoc ;
Vector3d vtD = vtDLoc ;
ptP.ToLoc( m_MapFrame[0]) ;
vtD.ToLoc( m_MapFrame[0]) ;
if ( bExact)
return GetDepthWithVoxel( ptP, vtD, dInLength, dOutLength, true) ;
else
return GetDepthWithDexel( ptP, vtD, dInLength, dOutLength) ;
}
//----------------------------------------------------------------------------
// Calcola la profondità del materiale usando i parallelepipedi corrispondenti ai segmenti di dexel di una griglia.
// Punto e versore devono essere espressi nel sistema di riferimento Zmap. Il vettore deve essere normalizzato.
// InLength = distanza di ingresso (se -1 il punto è interno, se -2 il punto è esterno e il raggio non interseca lo Zmap)
// OutLength = distanza di uscita
bool
VolZmap::GetDepthWithDexel( const Point3d& ptP, const Vector3d& vtV, double& dInLength, double& dOutLength) const
{
int nGrid = 0 ;
// Porto il raggio nel riferimento intrinseco
Point3d ptP = ptPGlob ;
ptP.ToLoc( m_MapFrame[nGrid]) ;
Vector3d vtV = vtDir ;
vtV.ToLoc( m_MapFrame[nGrid]) ;
vtV.Normalize() ;
// Intersezione fra semiretta e BBox dello Zmap
double dU1, dU2 ;
bool bTest = IntersLineZMapBBox( ptP, vtV, nGrid, dU1, dU2) && ( dU1 > 0 || dU2 > 0) ;
@@ -286,23 +308,17 @@ VolZmap::GetDepth( const Point3d& ptPGlob, const Vector3d& vtDir, double& dInLen
dInLength = - 1 ;
return true ;
#endif
}
//----------------------------------------------------------------------------
// Punti e vettori devono essere espressi nel sistema locale (quello in cui è immerso lo Zmap)
// Calcola la profondità del materiale usando i triangoli della rappresentazione grafica.
// Punto e versore devono essere espressi nel sistema di riferimento Zmap.
// Se si vuole che vengano utilizzati i triangoli con sharp-featue bEnh deve valere true, false altrimenti.
// InLength = distanza di ingresso (se -1 il punto è interno, se -2 il punto è esterno e il raggio non interseca lo Zmap)
// OutLength = distanza di uscita
// OutLength = distanza di uscita.
bool
VolZmap::GetDepthWithVoxel( const Point3d& ptPLoc, const Vector3d& vtDir, double& dInLength, double& dOutLength) const
VolZmap::GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtD, double& dInLength, double& dOutLength, bool bEnh) const
{
// Porto punto e vettore della retta nel sistema Zmap
Point3d ptP = ptPLoc ;
ptP.ToLoc( m_MapFrame[0]) ;
Vector3d vtD = vtDir ;
vtD.ToLoc( m_MapFrame[0]) ;
vtD.Normalize() ;
// Parametri di intersezione retta BBox dello Zmap
double dU1, dU2 ;
@@ -322,106 +338,120 @@ VolZmap::GetDepthWithVoxel( const Point3d& ptPLoc, const Vector3d& vtDir, double
if ( ! GetPointVoxel( ptP + dU1 * vtD, nVoxI, nVoxJ, nVoxK))
return false ;
}
// Indici dell'ultimo voxel
int nVxEndI, nVxEndJ, nVxEndK ;
// Determino il voxel finale
if ( ! GetPointVoxel( ptP + dU2 * vtD, nVxEndI, nVxEndJ, nVxEndK))
return false ;
// Vettore di voxel
vector<Voxel> vVox ;
// Struttura studio dell'intersezione
// Incrementi degli indici per sguire la retta
int nDeltaI = ( vtD.x > 0 ? 1 : ( vtD.x < 0 ? - 1 : 0)) ;
int nDeltaJ = ( vtD.y > 0 ? 1 : ( vtD.y < 0 ? - 1 : 0)) ;
int nDeltaK = ( vtD.z > 0 ? 1 : ( vtD.z < 0 ? - 1 : 0)) ;
// Scelgo i piani del voxel con cui intersecare la retta,
// per determinare il voxel successivo
int nPlaneI = ( vtD.x >= 0 ? 1 : 0) ;
int nPlaneJ = ( vtD.y >= 0 ? 1 : 0) ;
int nPlaneK = ( vtD.z >= 0 ? 1 : 0) ;
// Ciclo sui voxel
while ( IsValidVoxel( nVoxI, nVoxJ, nVoxK)) {
// Estremi della diagonale del voxel corrente
Point3d ptMin( ( nVoxI + 0.5) * m_dStep,
( nVoxJ + 0.5) * m_dStep,
( nVoxK + 0.5) * m_dStep) ;
Point3d ptMax( ( nVoxI + 1.5) * m_dStep,
( nVoxJ + 1.5) * m_dStep,
( nVoxK + 1.5) * m_dStep) ;
// Studio il voxel corrente
if ( IntersLineBox( ptP, vtD, ptMin, ptMax)) {
Voxel NewVox ;
NewVox.nI = nVoxI ;
NewVox.nJ = nVoxJ ;
NewVox.nK = nVoxK ;
vVox.emplace_back( NewVox) ;
}
// Interseco la retta con i piani frontiera del voxel
double dMaxTX = ( abs( vtD.x) > EPS_ZERO ? abs( ( ( nVoxI + nPlaneI + 0.5) * m_dStep - ptP.x) / vtD.x) : INFINITO) ;
double dMaxTY = ( abs( vtD.y) > EPS_ZERO ? abs( ( ( nVoxJ + nPlaneJ + 0.5) * m_dStep - ptP.y) / vtD.y) : INFINITO) ;
double dMaxTZ = ( abs( vtD.z) > EPS_ZERO ? abs( ( ( nVoxK + nPlaneK + 0.5) * m_dStep - ptP.z) / vtD.z) : INFINITO) ;
if ( dMaxTX < dMaxTY) {
if ( dMaxTX < dMaxTZ)
nVoxI += nDeltaI ;
else
nVoxK += nDeltaK ;
}
else {
if ( dMaxTY < dMaxTZ)
nVoxJ += nDeltaJ ;
else
nVoxK += nDeltaK ;
}
}
TRIA3DLIST lstTria ;
ExtMarchingCubes( vVox, lstTria, bEnh) ;
// Struttura studio dell'intersezione
struct LineTriaInt {
double dPar ;
double dDot ;
} ;
std::vector<std::vector<LineTriaInt>> IntMatrix ;
vector<vector<LineTriaInt>> IntMatrix ;
int nStI = min( nVoxI, nVxEndI) ;
int nStJ = min( nVoxJ, nVxEndJ) ;
int nStK = min( nVoxK, nVxEndK) ;
int nEnI = max( nVoxI, nVxEndI) ;
int nEnJ = max( nVoxJ, nVxEndJ) ;
int nEnK = max( nVoxK, nVxEndK) ;
for ( int nI = nStI ; nI <= nEnI ; ++ nI) {
for ( int nJ = nStJ ; nJ <= nEnJ ; ++ nJ) {
for ( int nK = nStK ; nK <= nEnK ; ++ nK) {
Point3d ptMin( ( nI + 0.5) * m_dStep,
( nJ + 0.5) * m_dStep,
( nK + 0.5) * m_dStep) ;
Point3d ptMax( ( nI + 1.5) * m_dStep,
( nJ + 1.5) * m_dStep,
( nK + 1.5) * m_dStep) ;
if ( IsValidVoxel( nI, nJ, nK) &&
IntersLineBox( ptP, vtD, ptMin, ptMax)) {
// Analisi del voxel
int nCbType ;
TRIA3DLIST lstTria ;
ProcessCube( nI, nJ, nK, nCbType, lstTria) ;
// Se il voxel contiene triangoli
if ( nCbType == VOX_ON_BOUNDARY) {
// Ciclo sui triangoli del voxel
for ( auto it = lstTria.begin() ; it != lstTria.end() ; ++it ) {
Triangle3d CurrTria = *it ;
Point3d ptLineTria1, ptLineTria2 ;
// Studio dell'intersezione della retta con il triangolo corrente
int nIntType = IntersLineTria( ptP, vtD, 1.5 * dU2, CurrTria, ptLineTria1, ptLineTria2) ;
// Se non ci sono intersezioni passo al prossimo triangolo
if ( nIntType == ILTT_NO)
continue ;
// se altrimenti c'è una sola intersezione
else if ( nIntType == ILTT_VERT ||
nIntType == ILTT_EDGE ||
nIntType == ILTT_IN) {
LineTriaInt NewInt ;
NewInt.dPar = ( ptLineTria1 - ptP) * vtD ;
NewInt.dDot = vtD * CurrTria.GetN() ;
std::vector<LineTriaInt> vSing ;
vSing.emplace_back( NewInt) ;
// Ciclo sui triangoli del voxel
for ( auto it = lstTria.begin() ; it != lstTria.end() ; ++it ) {
// Triangolo corrente e suoi punti di intersezione con la retta
Triangle3d CurrTria = *it ;
Point3d ptLineTria1, ptLineTria2 ;
// Studio dell'intersezione della retta con il triangolo corrente
int nIntType = IntersLineTria( ptP, vtD, 1.5 * dU2, CurrTria, ptLineTria1, ptLineTria2) ;
// Se non ci sono intersezioni passo al prossimo triangolo
if ( nIntType == ILTT_NO)
continue ;
// se altrimenti c'è una sola intersezione
else if ( nIntType == ILTT_VERT ||
nIntType == ILTT_EDGE ||
nIntType == ILTT_IN) {
LineTriaInt NewInt ;
NewInt.dPar = ( ptLineTria1 - ptP) * vtD ;
NewInt.dDot = vtD * CurrTria.GetN() ;
vector<LineTriaInt> vSing ;
vSing.emplace_back( NewInt) ;
IntMatrix.emplace_back( vSing) ;
}
// altrimenti ci sono due intersezioni
else {
LineTriaInt NewInt1, NewInt2 ;
std::vector<LineTriaInt> vCouple ;
double dP1 = ( ptLineTria1 - ptP) * vtD ;
double dP2 = ( ptLineTria2 - ptP) * vtD ;
NewInt1.dPar = ( dP1 < dP2 ? dP1 : dP2) ;
NewInt2.dPar = ( dP1 < dP2 ? dP2 : dP1) ;
NewInt1.dDot = vtD * CurrTria.GetN() ;
NewInt2.dDot = NewInt1.dDot ;
IntMatrix.emplace_back( vSing) ;
}
// altrimenti ci sono due intersezioni
else {
LineTriaInt NewInt1, NewInt2 ;
vector<LineTriaInt> vCouple ;
double dP1 = ( ptLineTria1 - ptP) * vtD ;
double dP2 = ( ptLineTria2 - ptP) * vtD ;
NewInt1.dPar = ( dP1 < dP2 ? dP1 : dP2) ;
NewInt2.dPar = ( dP1 < dP2 ? dP2 : dP1) ;
NewInt1.dDot = vtD * CurrTria.GetN() ;
NewInt2.dDot = NewInt1.dDot ;
vCouple.emplace_back( NewInt1) ;
vCouple.emplace_back( NewInt2) ;
IntMatrix.emplace_back( vCouple) ;
}
}
}
}
}
vCouple.emplace_back( NewInt1) ;
vCouple.emplace_back( NewInt2) ;
IntMatrix.emplace_back( vCouple) ;
}
}
// Ordino le intersezioni in base al parametro distanza con segno da ptP
for ( int nN = 0 ; nN < int( IntMatrix.size()) - 1 ; ++ nN) {
for ( int nM = nN + 1 ; nM < int( IntMatrix.size()) ; ++ nM) {
double dFP = ( IntMatrix[nN].size() == 2 ? 0.5 * ( IntMatrix[nN][0].dPar + IntMatrix[nN][1].dPar) :
IntMatrix[nN][0].dPar) ;
double dLP = ( IntMatrix[nN+1].size() == 2 ? 0.5 * ( IntMatrix[nN+1][0].dPar + IntMatrix[nN+1][1].dPar) :
IntMatrix[nN+1][0].dPar) ;
double dFP = ( IntMatrix[nN].size() == 2 ? 0.5 * ( IntMatrix[nN][0].dPar + IntMatrix[nN][1].dPar) :
IntMatrix[nN][0].dPar) ;
double dLP = ( IntMatrix[nM].size() == 2 ? 0.5 * ( IntMatrix[nM][0].dPar + IntMatrix[nM][1].dPar) :
IntMatrix[nM][0].dPar) ;
if ( dFP > dLP)
swap( IntMatrix[nN], IntMatrix[nN+1]) ;
if ( dFP > dLP)
swap( IntMatrix[nN], IntMatrix[nM]) ;
}
}
std::vector<LineTriaInt> vInt ;
vector<LineTriaInt> vInt ;
for ( int nN = 0 ; nN < int( IntMatrix.size()) ; ++ nN) {
vInt.emplace_back( IntMatrix[nN][0]) ;
@@ -437,7 +467,7 @@ VolZmap::GetDepthWithVoxel( const Point3d& ptPLoc, const Vector3d& vtDir, double
int nFirstPosN ;
int nN = 0 ;
for ( ; nN < int( vInt.size()) ; ++ nN) {
if ( vInt[nN].dPar >= 0) {
if ( /*vInt[nN].dPar >= 0*/vInt[nN].dPar > - EPS_SMALL) { // SCEGLIERE FRA LE DUE OPZIONI
nFirstPosN = nN ;
break ;
}
@@ -454,8 +484,14 @@ VolZmap::GetDepthWithVoxel( const Point3d& ptPLoc, const Vector3d& vtDir, double
dInLength = -1 ;
}
else if ( nFirstPosN == 0) {
if ( vInt[nFirstPosN].dDot > - EPS_ZERO)
dInLength = -1 ;
if ( vInt[nFirstPosN].dDot > EPS_ZERO) //
dInLength = -1 ;
else if ( GetPointVoxel( ptP, nVoxI, nVoxJ, nVoxK)) {
int nCubeType = CalcIndex( nVoxI, nVoxJ, nVoxK) ;
if ( nCubeType == 255)
dInLength = -1 ;
}
}
for ( int nN = nFirstPosN ; nN < int( vInt.size()) ; ++ nN) {
@@ -1001,8 +1037,8 @@ VolZmap::IntersLineEllipticalCylinder( const Vector3d& vtLineDir, const Point3d&
ptP.ToLoc( CircFrame) ;
vtV.ToLoc( CircFrame) ;
std::vector <double> vdCoef(3) ;
std::vector <double> vdRoots ;
vector <double> vdCoef(3) ;
vector <double> vdRoots ;
vdCoef[0] = dSqCoef * ptP.x * ptP.x + ptP.y * ptP.y + ptP.z * ptP.z - 2 * dObCoef * ptP.x * ptP.y - dSqRad ;
vdCoef[1] = 2 * ( dSqCoef * vtV.x * ptP.x + vtV.y * ptP.y + vtV.z * ptP.z - dObCoef * ( vtV.x * ptP.y + vtV.y * ptP.x)) ;