Egalware/EgtGeomKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
5b707f1799bcdb28702a5ccb8450b5c8fc84235d
EgtGeomKernel/VolZmapCalculus.cpp
T
Dario Sassi 5b707f1799 EgtGeomKernel : 
- correzioni a Zmap per GetDepth.
2019-05-14 15:46:42 +00:00

1670 lines
67 KiB
C++
Raw Blame History

//----------------------------------------------------------------------------
// EgalTech 2015-2016
//----------------------------------------------------------------------------
// File : VolZmap.cpp Data : 22.01.15 Versione : 1.6a4
// Contenuto : Implementazione della classe Volume Zmap (tre griglie)
//
//
//
// Modifiche : 22.01.15 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------
//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "CurveLine.h"
#include "VolZmap.h"
#include "GeoConst.h"
#include "/EgtDev/Include/EGkIntersLineTria.h"
#include "/EgtDev/Include/EGkIntersLinePlane.h"
#include "/EgtDev/Include/EGkIntersLineSphere.h"
#include "/EgtDev/Include/EGkIntersPlaneTria.h"
#include "/EgtDev/Include/EGkChainCurves.h"
#include "/EgtDev/Include/ENkPolynomialRoots.h"
#include "/EgtDev/Include/EgtNumUtils.h"
using namespace std ;
//----------------------------------------------------------------------------
// Box e piano devono essere nello stesso riferimento. Il box deve essere axis aligned.
static bool
TestIntersPlaneBox( const BBox3d& b3Box, const Plane3d& plPlane)
{
// Calcolo le distanze con segno dei punti dal piano
Point3d ptE0 = b3Box.GetMin() ;
double dDist0 = DistPointPlane( ptE0, plPlane) ;
Point3d ptE1( b3Box.GetMax().x, b3Box.GetMin().y, b3Box.GetMin().z) ;
double dDist1 = DistPointPlane( ptE1, plPlane) ;
Point3d ptE2( b3Box.GetMax().x, b3Box.GetMax().y, b3Box.GetMin().z) ;
double dDist2 = DistPointPlane( ptE2, plPlane) ;
Point3d ptE3( b3Box.GetMin().x, b3Box.GetMax().y, b3Box.GetMin().z) ;
double dDist3 = DistPointPlane( ptE3, plPlane) ;
Point3d ptE4( b3Box.GetMin().x, b3Box.GetMin().y, b3Box.GetMax().z) ;
double dDist4 = DistPointPlane( ptE4, plPlane) ;
Point3d ptE5( b3Box.GetMax().x, b3Box.GetMin().y, b3Box.GetMax().z) ;
double dDist5 = DistPointPlane( ptE5, plPlane) ;
Point3d ptE6 = b3Box.GetMax() ;
double dDist6 = DistPointPlane( ptE6, plPlane) ;
Point3d ptE7( b3Box.GetMin().x, b3Box.GetMax().y, b3Box.GetMax().z) ;
double dDist7 = DistPointPlane( ptE7, plPlane) ;
// Distanze tutte positive
if ( dDist0 > EPS_SMALL && dDist1 > EPS_SMALL && dDist2 > EPS_SMALL && dDist3 > EPS_SMALL &&
dDist4 > EPS_SMALL && dDist5 > EPS_SMALL && dDist6 > EPS_SMALL && dDist7 > EPS_SMALL)
return false ;
// Distanze tutte negative
if ( dDist0 < - EPS_SMALL && dDist1 < - EPS_SMALL && dDist2 < - EPS_SMALL && dDist3 < - EPS_SMALL &&
dDist4 < - EPS_SMALL && dDist5 < - EPS_SMALL && dDist6 < - EPS_SMALL && dDist7 < - EPS_SMALL)
return false ;
// Il piano interseca il box
return true ;
}
//----------------------------------------------------------------------------
// Punti e vettore devono esse espressi nel medesimo sistema di riferimento.
// Il box è allineato con gli assi di tale sistema di riferimento.
// Viene restituito true in caso di intersezione, false altrimenti.
bool
VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax) const
{
double dU1, dU2 ;
return IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2) ;
}
//----------------------------------------------------------------------------
// Punti e vettore devono esse espressi nel medesimo sistema di riferimento.
// Il box è allineato con gli assi di tale sistema di riferimento.
// In dU1, dU2 vengono restituiti i parametri a cui si trovano le intersezioni.
// Viene restituito true in caso di intersezione, false altrimenti.
bool
VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV,
const Point3d& ptMin, const Point3d& ptMax, double& dU1, double& dU2) const
{
// Il box è allineato agli assi
dU1 = - INFINITO ;
dU2 = INFINITO ;
// confronto con piani YZ (perpendicolari ad asse X)
if ( vtV.x > EPS_ZERO) {
dU1 = max( dU1, ( ptMin.x - ptP.x) / vtV.x) ;
dU2 = min( dU2, ( ptMax.x - ptP.x) / vtV.x) ;
}
else if ( vtV.x < - EPS_ZERO) {
dU1 = max( dU1, ( ptMax.x - ptP.x) / vtV.x) ;
dU2 = min( dU2, ( ptMin.x - ptP.x) / vtV.x) ;
}
else if ( ptP.x < ptMin.x - EPS_SMALL || ptP.x > ptMax.x + EPS_SMALL)
return false ;
// confronto con piani ZX (perpendicolari ad asse Y)
if ( vtV.y > EPS_ZERO) {
dU1 = max( dU1, ( ptMin.y - ptP.y) / vtV.y) ;
dU2 = min( dU2, ( ptMax.y - ptP.y) / vtV.y) ;
}
else if ( vtV.y < - EPS_ZERO) {
dU1 = max( dU1, ( ptMax.y - ptP.y) / vtV.y) ;
dU2 = min( dU2, ( ptMin.y - ptP.y) / vtV.y) ;
}
else if ( ptP.y < ptMin.y - EPS_SMALL || ptP.y > ptMax.y + EPS_SMALL)
return false ;
// confronto con piani XY (perpendicolari ad asse Z)
if ( vtV.z > EPS_ZERO) {
dU1 = max( dU1, ( ptMin.z - ptP.z) / vtV.z) ;
dU2 = min( dU2, ( ptMax.z - ptP.z) / vtV.z) ;
}
else if ( vtV.z < - EPS_ZERO) {
dU1 = max( dU1, ( ptMax.z - ptP.z) / vtV.z) ;
dU2 = min( dU2, ( ptMin.z - ptP.z) / vtV.z) ;
}
else if ( ptP.z < ptMin.z - EPS_SMALL || ptP.z > ptMax.z + EPS_SMALL)
return false ;
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) ;
}
//----------------------------------------------------------------------------
// Trova le intersezioni di una retta, definita tramite un suo punto e il vettore della direzione,
// e il bounding box dello Zmap. Punto e vettore devono essere espressi nel sistema di riferimento
// intrinseco dello Zmap. Se vi è intersezione viene restitutito true, false altrimenti.
// Se vi è intersezione in dU1 e dU2 vengono restituiti i parametri a cui avvengono le intersezioni.
bool
VolZmap::IntersLineZMapBBox( const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) const
{
// Punti estremi del box dello Zmap
Point3d ptMin = ORIG ;
Point3d ptMax = ptMin + Vector3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0]) ;
return IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2) ;
}
//----------------------------------------------------------------------------
// Trova le intersezioni fra una retta e un dexel parallelepipedale (gli intervalli sono
// parallelepipdi). Punto e vettore devono essere espressi nel sistema grilgia che coincide
// con quello intrinseco dello Zmap solo nel caso della prima griglia. Per le altre griglie
// è necessario permutare ciclicamente le coordinate.
bool
VolZmap::IntersLineDexel( const Point3d& ptP, const Vector3d& vtV, int nGrid, int nI, int nJ,
double& dU1, double& dU2) const
{
// Determino l'indice del dexel e il numero di suoi intervalli
int nDexelPos = nJ * m_nNx[nGrid] + nI ;
int nDexelSize = int( m_Values[nGrid][nDexelPos].size()) ;
// Se non c'è materiale non devo fare alcunché
if ( nDexelSize == 0)
return false ;
// Determino estremi nel piano XY intrinseco del dexel
double dXmin = nI * m_dStep ;
double dYmin = nJ * m_dStep ;
double dXmax = ( nI + 1) * m_dStep ;
double dYmax = ( nJ + 1) * m_dStep ;
// ciclo sugli intervalli
dU1 = INFINITO ;
dU2 = - INFINITO ;
bool bInters = false ;
for ( int nIndex = 0 ; nIndex < nDexelSize ; nIndex += 1) {
// estremi del box del singolo intervallo
Point3d ptE1( dXmin, dYmin, m_Values[nGrid][nDexelPos][nIndex].dMin) ;
Point3d ptE2( dXmax, dYmax, m_Values[nGrid][nDexelPos][nIndex].dMax) ;
double dP1, dP2 ;
if ( IntersLineBox( ptP, vtV, ptE1, ptE2, dP1, dP2)) {
bInters = true ;
dU1 = min( dU1, dP1) ;
dU2 = max( dU2, dP2) ;
}
}
return bInters ;
}
//----------------------------------------------------------------------------
// Trova le intersezioni fra una semi-retta e un dexel parallelepipedale (gli intervalli sono
// parallelepipdi). Punto e vettore devono essere espressi nel sistema grilgia che coincide
// con quello intrinseco dello Zmap solo nel caso della prima griglia. Per le altre griglie
// è necessario permutare ciclicamente le coordinate.
bool
VolZmap::IntersRayDexel( const Point3d& ptP, const Vector3d& vtV, int nGrid, int nI, int nJ,
double& dU1, double& dU2) const
{
// Determino l'indice del dexel e il numero di suoi intervalli
int nDexelPos = nJ * m_nNx[nGrid] + nI ;
int nDexelSize = int( m_Values[nGrid][nDexelPos].size()) ;
// Se non c'è materiale non devo fare alcunché
if ( nDexelSize == 0)
return false ;
// Determino estremi nel piano XY intrinseco del dexel
double dXmin = nI * m_dStep ;
double dYmin = nJ * m_dStep ;
double dXmax = ( nI + 1) * m_dStep ;
double dYmax = ( nJ + 1) * m_dStep ;
// ciclo sugli intervalli
dU1 = INFINITO ;
dU2 = - INFINITO ;
bool bInters = false ;
for ( int nIndex = 0 ; nIndex < nDexelSize ; nIndex += 1) {
// estremi del box del singolo intervallo
Point3d ptE1( dXmin, dYmin, m_Values[nGrid][nDexelPos][nIndex].dMin) ;
Point3d ptE2( dXmax, dYmax, m_Values[nGrid][nDexelPos][nIndex].dMax) ;
double dP1, dP2 ;
if ( IntersLineBox( ptP, vtV, ptE1, ptE2, dP1, dP2)) {
if ( dP2 >= 0) {
bInters = true ;
if ( dP1 >= 0)
dU1 = min( dU1, dP1) ;
else
dU1 = - 1 ;
dU2 = max( dU2, dP2) ;
}
}
}
return bInters ;
}
//----------------------------------------------------------------------------
// 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& ptPLoc, const Vector3d& vtDLoc, double& dInLength, double& dOutLength, bool bExact) const
{
if ( bExact && m_nMapNum == 3)
return GetDepthWithVoxel( ptPLoc, vtDLoc, dInLength, dOutLength) ;
else {
return GetDepthWithDexel( ptPLoc, vtDLoc, 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& ptPLoc, const Vector3d& vtDLoc, double& dInLength, double& dOutLength) const
{
// Punto e direzione nel riferimento intrinseco
Point3d ptP = ptPLoc ;
Vector3d vtD = vtDLoc ;
ptP.ToLoc( m_MapFrame) ;
vtD.ToLoc( m_MapFrame) ;
vtD.Normalize() ;
// Intersezione fra semiretta e BBox dello Zmap
double dU1, dU2 ;
bool bTest = IntersLineZMapBBox( ptP, vtD, dU1, dU2) && ( dU1 > 0 || dU2 > 0) ;
// Semiretta esterna al box dello Zmap quindi esterna anche allo Zmap
if ( ! bTest) {
dInLength = -2 ;
dOutLength = -2 ;
return true ;
}
double dInLen[N_MAPS] ;
double dOutLen[N_MAPS] ;
// Ciclo sulle griglie
for ( int nGrid = 0 ; nGrid < m_nMapNum ; ++ nGrid) {
Point3d ptP0 = ptP ;
Vector3d vtV0 = vtD ;
Point3d ptI, ptF ;
// Una sola intersezione valida ( punto interno, intersezione valida 2)
if ( dU1 < 0 && dU2 > 0) {
ptI = ptP ;
ptF = ptP + dU2 * vtD ;
}
// due soluzioni valide ( punto esterno)
else {
ptI = ptP + dU1 * vtD ;
ptF = ptP + dU2 * vtD ;
}
// Passo dal sistema intrinseco alla griglia
if ( nGrid == 1) {
swap( ptP0.y, ptP0.z) ;
swap( ptP0.x, ptP0.z) ;
swap( vtV0.y, vtV0.z) ;
swap( vtV0.x, vtV0.z) ;
swap( ptI.y, ptI.z) ;
swap( ptI.x, ptI.z) ;
swap( ptF.y, ptF.z) ;
swap( ptF.x, ptF.z) ;
}
else if ( nGrid == 2) {
swap( ptP0.x, ptP0.z) ;
swap( ptP0.y, ptP0.z) ;
swap( vtV0.x, vtV0.z) ;
swap( vtV0.y, vtV0.z) ;
swap( ptI.x, ptI.z) ;
swap( ptI.y, ptI.z) ;
swap( ptF.x, ptF.z) ;
swap( ptF.y, ptF.z) ;
}
// Inizializzo distanze
dInLen[nGrid] = INFINITO ;
dOutLen[nGrid] = - INFINITO ;
// Determino asse di spostamento maggiore
bool bOnX = ( abs( ptF.y - ptI.y) <= abs( ptF.x - ptI.x)) ;
// Movimento crescente su asse di movimento principale
if ( ( bOnX && ptF.x < ptI.x) || ( ! bOnX && ptF.y < ptI.y) )
swap( ptI, ptF) ;
// Pendenza
double dDeltaX = ( ptF.x - ptI.x) ;
double dDeltaY = ( ptF.y - ptI.y) ;
double dM = 0 ;
if ( bOnX && abs( dDeltaX) > EPS_SMALL)
dM = dDeltaY / dDeltaX ;
else if ( ! bOnX && abs( dDeltaY) > EPS_SMALL)
dM = dDeltaX / dDeltaY ;
// Determinazione degli indici i j dei punti ptI e ptF
int nIi = Clamp( int( floor( ptI.x / m_dStep)), 0, m_nNx[nGrid] - 1) ;
int nIj = Clamp( int( floor( ptI.y / m_dStep)), 0, m_nNy[nGrid] - 1) ;
int nFi = Clamp( int( floor( ptF.x / m_dStep)), 0, m_nNx[nGrid] - 1) ;
int nFj = Clamp( int( floor( ptF.y / m_dStep)), 0, m_nNy[nGrid] - 1) ;
// mi muovo
int i = nIi ; int j = nIj ;
while ( ( bOnX && i <= nFi) || ( ! bOnX && j <= nFj)) {
// Eseguo controllo
double dU1, dU2 ;
if ( IntersRayDexel( ptP0, vtV0, nGrid, i, j, dU1, dU2)) {
dInLen[nGrid] = min( dInLen[nGrid], dU1) ;
dOutLen[nGrid] = max( dOutLen[nGrid], dU2) ;
}
// Calcolo spostamento (a destra o sopra)
if ( bOnX) {
double dMoveX = ( ( i + 1) * m_dStep - ptI.x) ;
double dMoveY = dMoveX * dM ;
double dY = ptI.y + dMoveY ;
int nNewJ = Clamp( int( floor( dY / m_dStep)), 0, m_nNy[nGrid] - 1) ;
if ( nNewJ != j)
j = nNewJ ;
else
++ i ;
}
else {
double dMoveY = ( ( j + 1) * m_dStep - ptI.y) ;
double dMoveX = dMoveY * dM ;
double dX = ptI.x + dMoveX ;
int nNewI = Clamp( int( floor( dX / m_dStep)), 0, m_nNx[nGrid] - 1) ;
if ( nNewI != i)
i = nNewI ;
else
++ j ;
}
}
// Se non abbiamo incontrato materiale
if ( dInLen[nGrid] > dOutLen[nGrid] - EPS_SMALL) {
dInLen[nGrid] = -2 ;
dOutLen[nGrid] = -2 ;
}
// Se parto dall'interno
else if ( dInLen[nGrid] < - EPS_SMALL)
dInLen[nGrid] = -1 ;
}
if ( m_nMapNum == 1) {
dInLength = dInLen[0] ;
dOutLength = dOutLen[0] ;
return true ;
}
if ( abs( dInLen[0] + 2) < EPS_SMALL &&
abs( dInLen[1] + 2) < EPS_SMALL &&
abs( dInLen[2] + 2) < EPS_SMALL) {
dInLength = -2 ;
dOutLength = -2 ;
return true ;
}
else {
dInLength = INFINITO ;
dOutLength = - INFINITO ;
for ( int nGr = 0 ; nGr < 3 ; ++ nGr) {
if ( dInLen[nGr] > -2 && dInLen[nGr] < dInLength)
dInLength = dInLen[nGr] ;
if ( dOutLen[nGr] > -2 && dOutLen[nGr] > dOutLength)
dOutLength = dOutLen[nGr] ;
}
}
return true ;
}
//----------------------------------------------------------------------------
// Calcola la profondità del materiale usando i triangoli della rappresentazione grafica.
// Punto e versore devono essere espressi nel sistema di riferimento 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::GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtD, double& dInLength, double& dOutLength) const
{
ILZIVECTOR vIntersInfo ;
if ( ! GetLineIntersection( ptP, vtD, vIntersInfo))
return false ;
if ( vIntersInfo.empty()) {
dInLength = -2. ;
dOutLength = -2. ;
return true ;
}
// Inizializzo le distanze di ingresso e uscita:
// dInLength diminuisce, dOutLength aumenta.
dInLength = INFINITO ;
dOutLength = -INFINITO ;
int nFirstPosN ;
for ( int nN = 0; nN < int( vIntersInfo.size()) ; ++ nN) {
if ( vIntersInfo[nN].dU > - EPS_SMALL) {
nFirstPosN = nN ;
break ;
}
}
if ( nFirstPosN == int( vIntersInfo.size())) {
dInLength = -2 ;
dOutLength = -2 ;
return true ;
}
if ( nFirstPosN > 0) {
if ( vIntersInfo[nFirstPosN - 1].trTria.GetN() * vtD < EPS_ZERO)
dInLength = -1 ;
}
else if ( nFirstPosN == 0) {
if ( vIntersInfo[nFirstPosN].trTria.GetN() * vtD > EPS_ZERO)
dInLength = -1 ;
else {
Point3d ptPi = ptP ;
ptPi.ToLoc( m_MapFrame) ;
int nVoxIJK[3] ;
GetBlockIJKFromN( vIntersInfo[0].nVox, nVoxIJK) ;
if ( GetPointVoxel( ptPi, nVoxIJK[0], nVoxIJK[1], nVoxIJK[2])) {
int nCubeType = CalcIndex( nVoxIJK[0], nVoxIJK[1], nVoxIJK[2]) ;
if ( nCubeType == 255)
dInLength = -1 ;
}
}
}
for ( int nN = nFirstPosN ; nN < int( vIntersInfo.size()) ; ++ nN) {
if ( vIntersInfo[nN].trTria.GetN() * vtD > - EPS_ZERO) {
if ( ( vIntersInfo[nN].nILTT == ILTT_SEGM || vIntersInfo[nN].nILTT == ILTT_SEGM_ON_EDGE) &&
dOutLength < vIntersInfo[nN].dU2)
dOutLength = vIntersInfo[nN].dU2 ;
else if ( dOutLength < vIntersInfo[nN].dU)
dOutLength = vIntersInfo[nN].dU ;
}
if ( vIntersInfo[nN].trTria.GetN() * vtD < EPS_ZERO &&
dInLength > vIntersInfo[nN].dU)
dInLength = vIntersInfo[nN].dU ;
}
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::AvoidBox( const Frame3d& frBox, const Vector3d& vtDiag, double dSafeDist) const
{
// BBox
BBox3d b3Box( ORIG, ORIG + vtDiag) ;
// Sistemazioni per distanza di sicurezza
if ( dSafeDist > EPS_SMALL)
b3Box.Expand( dSafeDist) ;
// Riferimento nell'angolo con coordinate minime
Point3d ptMin = b3Box.GetMin() ; ptMin.ToGlob( frBox) ;
Frame3d frB = frBox ; frB.ChangeOrig( ptMin) ;
Vector3d vtDg = b3Box.GetMax() - b3Box.GetMin() ;
b3Box.Set( ORIG, ORIG + vtDg) ;
// Lo porto nel riferimento intrinseco dello Zmap
b3Box.LocToLoc( frB, m_MapFrame) ;
// BBox dello Zmap nel suo riferimento intrinseco
BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;
// Se non interferiscono, posso uscire
BBox3d b3Int ;
if ( ! b3Zmap.FindIntersection( b3Box, b3Int))
return true ;
// Limiti su indici
int nStI = Clamp( int( b3Int.GetMin().x / m_dStep), 0, m_nNx[0] - 1) ;
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) ;
// Vettore direzione dei dexel nel riferimento del Box
Vector3d vtK = Z_AX ; vtK.LocToLoc( m_MapFrame, frB) ;
// Riferimento intrinseco dei dexel nel riferimento del box
Point3d ptO = ORIG ; ptO.LocToLoc( m_MapFrame, frB) ;
Vector3d vtX = X_AX ; vtX.LocToLoc( m_MapFrame, frB) ;
Vector3d vtY = Y_AX ; vtY.LocToLoc( m_MapFrame, frB) ;
// Ciclo di intersezione dei dexel con il BBox
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 ;
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 ;
double dZmin, dZmax ;
if ( IntersLineBox( ptT, vtK, ORIG, ORIG + vtDg, 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 true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::AvoidSphere( const Point3d& ptCenter, double dRad, double dSafeDist) const
{
// Porto la sfera nel riferimento intrinseco dello Zmap
Point3d ptC = ptCenter ;
ptC.ToLoc( m_MapFrame) ;
// Aumento il raggio della distanza di sicurezza
if ( dSafeDist > EPS_SMALL)
dRad += dSafeDist ;
// BBox della sfera
BBox3d b3Box( ptC) ;
b3Box.Expand( dRad) ;
// BBox dello Zmap nel suo riferimento intrinseco
BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;
// Se non interferiscono, posso uscire
BBox3d b3Int ;
if ( ! b3Zmap.FindIntersection( b3Box, b3Int))
return true ;
// Limiti su indici
int nStI = Clamp( int( b3Int.GetMin().x / m_dStep), 0, m_nNx[0] - 1) ;
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) ;
// Ciclo di intersezione dei dexel con la sfera (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 ;
for ( int k = 0 ; k < 5 ; ++ k) {
Point3d ptT = ORIG + ( i + 0.5) * m_dStep * X_AX + ( j + 0.5) * m_dStep * Y_AX ;
if ( k == 0)
;
else if ( k == 1)
ptT += - 0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ;
else if ( k == 2)
ptT += + 0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ;
else if ( k == 3)
ptT += + 0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ;
else if ( k == 4)
ptT += - 0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ;
Point3d ptI1, ptI2 ;
if ( ::IntersLineSphere( ptT, Z_AX, ptC, dRad, ptI1, ptI2) != ILST_NO) {
double dZmin = min( ptI1.z, ptI2.z) ;
double dZmax = max( ptI1.z, ptI2.z) ;
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 true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::AvoidCylinder( const Frame3d& frCyl, double dH, double dR, double dSafeDist) const
{
// Porto il cilindro nel riferimento intrinseco dello Zmap
Frame3d frC = frCyl ;
frC.ToLoc( m_MapFrame) ;
// Se altezza negativa, sposto riferimento da faccia sopra a quella sotto
if ( dH < 0) {
frC.Translate( dH * frC.VersZ()) ;
dH = - dH ;
}
// Sistemazioni per distanza di sicurezza
if ( dSafeDist > EPS_SMALL) {
frC.Translate( - dSafeDist * frC.VersZ()) ;
dH += 2 * dSafeDist ;
dR += dSafeDist ;
}
// BBox del cilindro
Vector3d vtDirL = frC.VersZ() ;
BBox3d b3Box( frC.Orig()) ;
b3Box.Add( frC.Orig() + frC.VersZ() * dH) ;
if ( vtDirL.IsXplus() || vtDirL.IsXminus())
b3Box.Expand( 0, dR, dR) ;
else if ( vtDirL.IsYplus() || vtDirL.IsYminus())
b3Box.Expand( dR, 0, dR) ;
else if ( vtDirL.IsZplus() || vtDirL.IsZminus())
b3Box.Expand( dR, dR, 0) ;
else {
double dExpandX = dR * sqrt( 1 - vtDirL.x * vtDirL.x) ;
double dExpandY = dR * sqrt( 1 - vtDirL.y * vtDirL.y) ;
double dExpandZ = dR * sqrt( 1 - vtDirL.z * vtDirL.z) ;
b3Box.Expand( dExpandX, dExpandY, dExpandZ) ;
}
// BBox dello Zmap nel suo riferimento intrinseco
BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;
// Se non interferiscono, posso uscire
BBox3d b3Int ;
if ( ! b3Zmap.FindIntersection( b3Box, b3Int))
return true ;
// Limiti su indici
int nStI = Clamp( int( b3Int.GetMin().x / m_dStep), 0, m_nNx[0] - 1) ;
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) ;
// 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) {
int nPos = j * m_nNx[0] + i ;
int nSize = int( m_Values[0][nPos].size()) ;
if ( nSize == 0)
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 ;
if ( k == 0)
;
else if ( k == 1)
ptT += - 0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ;
else if ( k == 2)
ptT += + 0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ;
else if ( k == 3)
ptT += + 0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ;
else if ( k == 4)
ptT += - 0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ;
Point3d ptI1, ptI2 ;
Vector3d vtN1, vtN2 ;
if ( IntersLineCylinder( ptT, Z_AX, frC, dH, dR, true, true, ptI1, vtN1, ptI2, vtN2)) {
double dZmin = min( ptI1.z, ptI2.z) ;
double dZmax = max( ptI1.z, ptI2.z) ;
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 true ;
}
//----------------------------------------------------------------------------
// Riferimento con origine nel centro della base e asse di simmetria coincidente con l'asse Z.
// La funzione restituisce true in caso di intersezione, false altrimenti.
//----------------------------------------------------------------------------
bool
VolZmap::IntersLineCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir,
const Frame3d& CylFrame, double dH, double dRad, bool bTapLow, bool bTapUp,
Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const
{
// Porto la linea nel riferimento del cilindro
Point3d ptP = ptLineSt ; ptP.ToLoc( CylFrame) ;
Vector3d vtV = vtLineDir ; vtV.ToLoc( CylFrame) ;
// 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( vtV.z) > EPS_ZERO) {
// le linee tangenti al cilindro non sono considerate intersecanti
double EpsRad = ( vtV.IsZeroXY() ? - EPS_SMALL : EPS_SMALL) ;
ptInt1 = ptP + ( ( 0 - ptP.z) / vtV.z) * vtV ;
if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y < dRad * dRad + 2 * dRad * EpsRad) {
nBasInt += 1 ;
vtN1 = Z_AX ;
}
ptInt2 = ptP + ( ( dH - ptP.z) / vtV.z) * vtV ;
if ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y < dRad * dRad + 2 * dRad * EpsRad) {
nBasInt += 2 ;
vtN2 = - Z_AX ;
}
}
// Se la linea interseca entrambe le basi, si sono trovate le due intersezioni
if ( nBasInt == 3) {
// Porto i punti e i versori nel riferimento globale
ptInt1.ToGlob( CylFrame) ;
vtN1.ToGlob( CylFrame) ;
ptInt2.ToGlob( CylFrame) ;
vtN2.ToGlob( CylFrame) ;
// Trovate intersezioni
return true ;
}
// Determino le intersezioni con la superficie laterale del cilindro
DBLVECTOR vdCoef(3) ;
double dSqRad = dRad * dRad ;
vdCoef[0] = ptP.x * ptP.x + ptP.y * ptP.y - dSqRad ;
vdCoef[1] = 2 * ( ptP.x * vtV.x + ptP.y * vtV.y) ;
vdCoef[2] = vtV.x * vtV.x + vtV.y * vtV.y ;
DBLVECTOR vdRoots ;
int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ;
// Epsilon per piani di tappo
double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ;
double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ;
// Elimino le soluzioni cha danno intersezioni fuori dai limiti in Z del cilindro
if ( nRoot == 2) {
double dIntZ2 = ptP.z + vdRoots[1] * vtV.z ;
if ( dIntZ2 < 0 + dEpsLow || dIntZ2 > dH + dEpsUp)
nRoot = 1 ;
}
if ( nRoot >= 1) {
double dIntZ1 = ptP.z + vdRoots[0] * vtV.z ;
if ( dIntZ1 < 0 + dEpsLow || dIntZ1 > dH + dEpsUp) {
if ( nRoot == 2)
vdRoots[0] = vdRoots[1] ;
-- nRoot ;
}
}
// Due soluzioni: la retta interseca due volte la superficie laterale
if ( nRoot == 2) {
// Punti di intersezione con la superficie del cilindro
ptInt1 = ptP + vdRoots[0] * vtV ;
ptInt2 = ptP + vdRoots[1] * vtV ;
// Determino le normali
vtN1.Set( -ptInt1.x, -ptInt1.y, 0) ;
vtN1.Normalize() ;
vtN2.Set( -ptInt2.x, -ptInt2.y, 0) ;
vtN2.Normalize() ;
// Porto i punti e i versori nel riferimento globale
ptInt1.ToGlob( CylFrame) ;
vtN1.ToGlob( CylFrame) ;
ptInt2.ToGlob( CylFrame) ;
vtN2.ToGlob( CylFrame) ;
// 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) {
// Punto di intersezione
ptInt1 = ptP + vdRoots[0] * vtV ;
// Normale alla superficie del cilindro verso l'interno
vtN1.Set( -ptInt1.x, -ptInt1.y, 0) ;
vtN1.Normalize() ;
}
// altrimenti piano inferiore
else if ( nBasInt == 1) {
// Punto di intersezione
ptInt2 = ptP + vdRoots[0] * vtV ;
// Normale alla superficie del cilindro verso l'interno
vtN2.Set( -ptInt2.x, -ptInt2.y, 0) ;
vtN2.Normalize() ;
}
// altrimenti niente
else
return false ;
// Porto i punti e i versori nel riferimento globale
ptInt1.ToGlob( CylFrame) ;
vtN1.ToGlob( CylFrame) ;
ptInt2.ToGlob( CylFrame) ;
vtN2.ToGlob( CylFrame) ;
// Trovate intersezioni
return true ;
}
// Nessuna soluzione : nessuna intersezione
else
return false ;
}
//----------------------------------------------------------------------------
// Riferimento con origine nel vertice del cono e asse di simmetria coincidente con l'asse Z.
// La funzione restituisce true in caso di intersezione, false altrimenti.
//----------------------------------------------------------------------------
bool
VolZmap::IntersLineConus( const Point3d& ptLineSt, const Vector3d& vtLineDir,
const Frame3d& ConusFrame, double dTan, double dMinH, double dMaxH, bool bTapLow, bool bTapUp,
Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const
{
// Porto la linea nel riferimento del cono
Point3d ptP = ptLineSt ; ptP.ToLoc( ConusFrame) ;
Vector3d vtV = vtLineDir ; vtV.ToLoc( ConusFrame) ;
// Raggi delle due basi
double dMinRad = dTan * dMinH ;
double dMaxRad = dTan * dMaxH ;
// Epsilon per piani di tappo
double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ;
double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ;
// 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( vtV.z) > EPS_ZERO) {
ptInt1 = ptP + ( ( dMinH - ptP.z) / vtV.z) * vtV ;
if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y < dMinRad * dMinRad + 2 * dMinRad * dTan * dEpsLow) {
nBasInt += 1 ;
vtN1 = Z_AX ;
}
ptInt2 = ptP + ( ( dMaxH - ptP.z) / vtV.z) * vtV ;
if ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y < dMaxRad * dMaxRad + 2 * dMaxRad * dTan * dEpsUp) {
nBasInt += 2 ;
vtN2 = - Z_AX ;
}
}
// Se la linea interseca entrambe le basi, si sono trovate le due intersezioni
if ( nBasInt == 3) {
// Porto i punti e i versori nel riferimento globale
ptInt1.ToGlob( ConusFrame) ;
vtN1.ToGlob( ConusFrame) ;
ptInt2.ToGlob( ConusFrame) ;
vtN2.ToGlob( ConusFrame) ;
// Trovate intersezioni
return true ;
}
// Determino le intersezioni con la superficie laterale del cono
DBLVECTOR vdCoef( 3) ;
double dSqTan = dTan * dTan ;
vdCoef[0] = dSqTan * ptP.z * ptP.z - ptP.x * ptP.x - ptP.y * ptP.y ;
vdCoef[1] = 2 * ( dSqTan * ptP.z * vtV.z - ptP.x * vtV.x - ptP.y * vtV.y) ;
vdCoef[2] = dSqTan * vtV.z * vtV.z - vtV.x * vtV.x - vtV.y * vtV.y ;
DBLVECTOR vdRoots ;
int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ;
// Elimino le soluzioni cha danno intersezioni fuori dai limiti in Z del tronco
if ( nRoot == 2) {
double dIntZ2 = ptP.z + vdRoots[1] * vtV.z ;
if ( dIntZ2 < dMinH + dEpsLow || dIntZ2 > dMaxH + dEpsUp)
nRoot = 1 ;
}
if ( nRoot >= 1) {
double dIntZ1 = ptP.z + vdRoots[0] * vtV.z ;
if ( dIntZ1 < dMinH + dEpsLow || dIntZ1 > dMaxH + dEpsUp) {
if ( nRoot == 2)
vdRoots[0] = vdRoots[1] ;
-- nRoot ;
}
}
// Due soluzioni: la retta interseca due volte la superficie laterale
if ( nRoot == 2) {
// Punti di intersezione con la superficie del cono
ptInt1 = ptP + vdRoots[0] * vtV ;
ptInt2 = ptP + vdRoots[1] * vtV ;
if ( ptInt1.z > ptInt2.z)
swap( ptInt1, ptInt2) ;
// Determino le normali
vtN1.Set( -ptInt1.x, -ptInt1.y, ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y) / ptInt1.z) ;
vtN1.Normalize() ;
vtN2.Set( -ptInt2.x, -ptInt2.y, ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y) / ptInt2.z) ;
vtN2.Normalize() ;
// Porto i punti e i versori nel riferimento globale
ptInt1.ToGlob( ConusFrame) ;
vtN1.ToGlob( ConusFrame) ;
ptInt2.ToGlob( ConusFrame) ;
vtN2.ToGlob( ConusFrame) ;
// 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) {
// Punto di intersezione
ptInt1 = ptP + vdRoots[0] * vtV ;
// Normale alla superficie del cono verso l'interno
vtN1.Set( -ptInt1.x, -ptInt1.y, ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y) / ptInt1.z) ;
vtN1.Normalize() ;
}
// altrimenti piano inferiore
else if( nBasInt == 1) {
// Punto di intersezione
ptInt2 = ptP + vdRoots[0] * vtV ;
// Normale alla superficie del cono verso l'interno
vtN2.Set( -ptInt2.x, -ptInt2.y, ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y) / ptInt2.z) ;
vtN2.Normalize() ;
}
// altrimenti niente
else
return false ;
// Porto i punti e i versori nel riferimento globale
ptInt1.ToGlob( ConusFrame) ;
vtN1.ToGlob( ConusFrame) ;
ptInt2.ToGlob( ConusFrame) ;
vtN2.ToGlob( ConusFrame) ;
// Trovate intersezioni
return true ;
}
// Nessuna soluzione : nessuna intersezione
else
return false ;
}
//----------------------------------------------------------------------------
// L'origine del sistema di riferimento deve essere
// nel centro della circonferenza di base, la cui traslazione obliqua
// genera il cilindro ellittico, e l'asse z deve essere l'asse
// di simmetria di tale circonferenza.
// La funzione restituisce true in caso di intersezione, false altrimenti.
// NB: dSqRad è il quadrato del raggio della circonferenza la cui
// traslazione obliqua genera il cilindro ellittico, dLongMvLen e
// dOrtMvLen sono rispettivamente le lunghezze delle proiezioni del
// movimento su z e x del sistema di riferimento CircFrame.
bool
VolZmap::IntersLineEllipticalCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir,
const Frame3d& CircFrame, double dRad, double dLongMvLen, double dOrtMvLen,
bool bTapLow, bool bTapUp,
Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const
{
// Se il cilindrico ellittico degenera in un piano, non bisogna tagliare
if ( abs( dLongMvLen) < EPS_SMALL)
return false ;
// Porto la linea nel riferimento del cilindro
Point3d ptP = ptLineSt ; ptP.ToLoc( CircFrame) ;
Vector3d vtV = vtLineDir ; vtV.ToLoc( CircFrame) ;
// Quadrato del raggio
double dSqRad = dRad * dRad ;
// Asse cilindro ellittico
Vector3d vtAx( dOrtMvLen, 0, dLongMvLen) ;
vtAx.Normalize() ;
// Retta parallela all'asse del cilindro
if ( AreSameOrOppositeVectorApprox( vtV, vtAx)) {
// Interseco la retta con i piani delle circonferenze
Point3d ptOLsCirc( dOrtMvLen, 0, dLongMvLen) ;
ptInt1 = ptP - ( ptP.z / vtV.z) * vtV ;
ptInt2 = ptP - ( ( ptP.z - dLongMvLen) / vtV.z) * vtV ;
double dSafeSqRad = dSqRad - 2 * dRad * EPS_SMALL ;
if ( ( ptInt1 - ORIG).SqLenXY() < dSafeSqRad &&
( ptInt2 - ptOLsCirc).SqLenXY() < dSafeSqRad) {
vtN1 = Z_AX ;
vtN2 = - Z_AX ;
ptInt1.ToGlob( CircFrame) ;
ptInt2.ToGlob( CircFrame) ;
vtN1.ToGlob( CircFrame) ;
vtN2.ToGlob( CircFrame) ;
return true ;
}
else
return false ;
}
// Coefficiente angolare della retta di movimento nel
// piano ZX del sistema di riferimento del movimento e suo quadrato
double dObCoef = dOrtMvLen / dLongMvLen ;
double dSqCoef = dObCoef * dObCoef ;
// Setto i coeficienti dell'equazione
DBLVECTOR vdCoef(3) ;
vdCoef[0] = dSqCoef * ptP.z * ptP.z + ptP.x * ptP.x + ptP.y * ptP.y
- 2 * dObCoef * ptP.z * ptP.x - dSqRad ;
vdCoef[1] = 2 * ( dSqCoef * vtV.z * ptP.z + vtV.x * ptP.x + vtV.y * ptP.y
- dObCoef * ( vtV.z * ptP.x + vtV.x * ptP.z)) ;
vdCoef[2] = dSqCoef * vtV.z * vtV.z + vtV.x * vtV.x + vtV.y * vtV.y
- 2 * dObCoef * vtV.z * vtV.x ;
// Numero di soluzioni
DBLVECTOR vdRoots ;
int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ;
// L'equazione ammette o due soluzioni (eventualmente
// coincidenti) oppure nessuna o infinite se la la retta
// appartiene alla superficie
// Se numero di soluzioni diverso da due le eventuali intersezioni sono già state trovate
if ( nRoot != 2)
return false ;
// Due soluzioni trovate
// Vettore di movimento
Vector3d vtMv( dOrtMvLen, 0, dLongMvLen) ;
// Punti di intersezione
ptInt1 = ptP + vdRoots[0] * vtV ;
ptInt2 = ptP + vdRoots[1] * vtV ;
// Simmetria del problema
if ( ptInt1.z > ptInt2.z)
swap( ptInt1, ptInt2) ;
// Determino le normali alla superficie nei punti d'intersezione
// Vettore 1
Vector3d vtTest1 = ( ptInt1 - ORIG) - ( ptInt1 - ORIG) * vtAx * vtAx ;
double dX0_1 = ( vtTest1.x > 0 ? 1 : -1) * sqrt( max( dSqRad - ptInt1.y * ptInt1.y, 0.)) ;
Vector3d vtCirc1( - dX0_1, - ptInt1.y, 0) ;
Vector3d vtTan1( vtCirc1.y, - vtCirc1.x, 0) ;
Vector3d vtCross1 = vtTan1 ^ vtMv ;
vtN1 = ( vtCross1 * vtCirc1 > - EPS_ZERO ? vtCross1 : - vtCross1) ;
vtN1.Normalize() ;
// Vettore 2
Vector3d vtTest2 = ( ptInt2 - ORIG) - ( ptInt2 - ORIG) * vtAx * vtAx ;
double dX0_2 = ( vtTest2.x > 0 ? 1 : -1) * sqrt( max( dSqRad - ptInt2.y * ptInt2.y, 0.)) ;
Vector3d vtCirc2( - dX0_2, - ptInt2.y, 0) ;
Vector3d vtTan2( vtCirc2.y, - vtCirc2.x, 0) ;
Vector3d vtCross2 = vtTan2 ^ vtMv ;
vtN2 = ( vtCross2 * vtCirc2 > - EPS_ZERO ? vtCross2 : - vtCross2) ;
vtN2.Normalize() ;
// Flag per i tappi
double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ;
double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ;
// Studio le soluzioni: se una è fuori dalla regione
// ammissibile, vuol dire che la retta esce da un tappo.
if ( ptInt1.z < dLongMvLen + dEpsUp) {
if ( ptInt1.z > + dEpsLow) {
// ptInt2 è sul tappo
if ( ptInt2.z > dLongMvLen + dEpsUp) {
ptInt2 = ptP + ( ( dLongMvLen - ptP.z) / vtV.z) * vtV ;
vtN2 = - Z_AX ;
}
}
else {
// Entrambe le soluzioni sono su un tappo
if ( ptInt2.z > dLongMvLen + dEpsUp) {
ptInt1 = ptP - ( ptP.z / vtV.z) * vtV ;
ptInt2 = ptP + ( ( dLongMvLen - ptP.z) / vtV.z) * vtV ;
vtN1.Set( 0, 0, 1) ;
vtN2.Set( 0, 0, -1) ;
if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y > dSqRad &&
ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y > dSqRad)
return false ;
}
// La prima soluzione è sul tappo
else if ( ptInt2.z > dEpsLow) {
ptInt1 = ptP - ( ptP.z / vtV.z) * vtV ;
vtN1.Set( 0, 0, 1) ;
}
else
return false ;
}
}
else
return false ;
// Porto i punti e i versori nel riferimento globale
ptInt1.ToGlob( CircFrame) ;
vtN1.ToGlob( CircFrame) ;
ptInt2.ToGlob( CircFrame) ;
vtN2.ToGlob( CircFrame) ;
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir,
const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaZ,
Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const
{
// Controllo sulle dimensioni lineari affinché sia valido il poliedro
if ( dLenX <= EPS_SMALL || dLenY <= EPS_SMALL || dLenZ <= EPS_SMALL)
return false ;
// Porto la linea nel riferimento del poliedro
Point3d ptP = ptLineSt ; ptP.ToLoc( PolyFrame) ;
Vector3d vtV = vtLineDir ; vtV.ToLoc( PolyFrame) ;
// Facce 1 e 2 parallele a XY
// Facce 3 e 4 parallele a XZ
// Facce 5 e 6 oblique
Point3d ptFacet135( 0, dLenY / 2, 0) ;
Point3d ptFacet246( dLenX, - dLenY / 2, dLenZ + dDeltaZ) ;
// Vettori per descrizione piani obliqui
Vector3d vtFacet5 = ptFacet135 - ptP ;
Vector3d vtFacet6 = ptFacet246 - ptP ;
Vector3d vtOb( - dDeltaZ, 0, dLenX) ;
vtOb.Normalize() ;
// Controlli affinché non vengano tagliati dexel a filo
// con il passaggio dell'utensile:
// Controllo sulle facce 1 e 2
if ( abs( vtV.y) < EPS_ZERO && abs( ptP.y) > ptFacet135.y)
return false ;
// Controllo sulle facce 3 e 4
if ( abs( vtV.x) < EPS_ZERO && ( ptP.x < ptFacet135.x || ptP.x > ptFacet246.x))
return false ;
// Controllo sulle facce 5 e 6
double dP1 = abs ( ( ptFacet135 - ptP) * vtOb) ;
double dP2 = abs ( ( ptFacet246 - ptP) * vtOb) ;
if ( abs( vtV * vtOb) < EPS_ZERO && ( dP1 < EPS_SMALL || dP2 < EPS_SMALL))
return false ;
// Punti notevoli
Point3d ptI1 = ptP + ( ( ptFacet135.y - ptP.y) / vtV.y) * vtV ;
Point3d ptI2 = ptP + ( ( ptFacet246.y - ptP.y) / vtV.y) * vtV ;
Point3d ptI3 = ptP + ( ( ptFacet135.x - ptP.x) / vtV.x) * vtV ;
Point3d ptI4 = ptP + ( ( ptFacet246.x - ptP.x) / vtV.x) * vtV ;
Point3d ptI5 = ptP + ( ( vtFacet5 * vtOb) / ( vtV * vtOb)) * vtV ;
Point3d ptI6 = ptP + ( ( vtFacet6 * vtOb) / ( vtV * vtOb)) * vtV ;
// Ricerca intersezioni con le facce
int nIntNum = 0 ;
// Intersezione con la prima faccia
if ( ptI1.x >= 0 && ptI1.x <= dLenX &&
ptI1.z * dLenX >= dDeltaZ * ptI1.x && ( ptI1.z - dLenZ) * dLenX <= dDeltaZ * ptI1.x) {
ptInt1 = ptI1 ;
vtN1 = - Y_AX ;
++ nIntNum ;
}
// Intersezione con la seconda faccia
if ( ptI2.x >= 0 && ptI2.x <= dLenX &&
ptI2.z * dLenX >= dDeltaZ * ptI2.x && ( ptI2.z - dLenZ) * dLenX <= dDeltaZ * ptI2.x) {
if ( nIntNum == 0) {
ptInt1 = ptI2 ;
vtN1 = Y_AX ;
++ nIntNum ;
}
else if ( ( ptInt1 - ptI2).SqLen() > SQ_EPS_SMALL) {
ptInt2 = ptI2 ;
vtN2 = Y_AX ;
++ nIntNum ;
}
}
// Intersezione con la terza faccia
if ( nIntNum < 2 &&
ptI3.z >= 0 && ptI3.z <= dLenZ &&
ptI3.y >= - ptFacet135.y && ptI3.y <= ptFacet135.y) {
if ( nIntNum == 0) {
ptInt1 = ptI3 ;
vtN1 = X_AX ;
++ nIntNum ;
}
else if ( ( ptInt1 - ptI3).SqLen() > SQ_EPS_SMALL) {
ptInt2 = ptI3 ;
vtN2 = X_AX ;
++ nIntNum ;
}
}
// Intersezione con la quarta faccia
if ( nIntNum < 2 &&
ptI4.z >= dDeltaZ && ptI4.z <= dLenZ + dDeltaZ &&
ptI4.y >= - ptFacet135.y && ptI4.y <= ptFacet135.y) {
if ( nIntNum == 0) {
ptInt1 = ptI4 ;
vtN1 = - X_AX ;
++ nIntNum ;
}
else if ( ( ptInt1 - ptI4).SqLen() > SQ_EPS_SMALL) {
ptInt2 = ptI4 ;
vtN2 = - X_AX ;
++ nIntNum ;
}
}
// Intersezione con la quinta faccia
if ( nIntNum < 2 &&
ptI5.x >= 0 && ptI5.x <= dLenX &&
ptI5.y >= - ptFacet135.y && ptI5.y <= ptFacet135.y) {
if ( nIntNum == 0) {
ptInt1 = ptI5 ;
vtN1 = vtOb ;
++ nIntNum ;
}
else if ( ( ptInt1 - ptI5).SqLen() > SQ_EPS_SMALL) {
ptInt2 = ptI5 ;
vtN2 = vtOb ;
++ nIntNum ;
}
}
// Intersezione con la sesta faccia
if ( nIntNum < 2 &&
ptI6.x >= 0 && ptI6.x <= dLenX &&
ptI6.y >= - ptFacet135.y && ptI6.y <= ptFacet135.y) {
if ( nIntNum == 0) {
ptInt1 = ptI6;
vtN1 = - vtOb ;
++ nIntNum ;
}
else if ( ( ptInt1 - ptI6).SqLen() > SQ_EPS_SMALL) {
ptInt2 = ptI6;
vtN2 = - vtOb ;
++ nIntNum ;
}
}
if ( nIntNum != 2)
return false ;
// Porto i punti e i versori nel riferimento globale
ptInt1.ToGlob( PolyFrame) ;
ptInt2.ToGlob( PolyFrame) ;
vtN1.ToGlob( PolyFrame) ;
vtN2.ToGlob( PolyFrame) ;
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetVolume( double& dVol) const
{
// verifico lo stato
if ( m_nStatus != OK)
return false ;
// Eseguo il calcolo della lunghezza totale degli spilloni
double dLen = 0 ;
for ( int nMap = 0 ; nMap < int( m_nMapNum) ; ++ nMap) {
for ( int nDex = 0 ; nDex < int( m_Values[nMap].size()) ; ++ nDex) {
for ( int nInt = 0 ; nInt < int( m_Values[nMap][nDex].size()) ; ++ nInt) {
dLen += m_Values[nMap][nDex][nInt].dMax - m_Values[nMap][nDex][nInt].dMin ;
}
}
}
// Il volume si trova moltiplicando per l'area di ogni spillone e dividendo per il numero di mappe
dVol = dLen * ( m_dStep * m_dStep) / m_nMapNum ;
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetPartVolume( int nPart, double& dVol) const
{
// verifico lo stato
if ( m_nStatus != OK)
return false ;
// Se una sola mappa o il numero di componenti è indefinito, errore
if ( m_nMapNum == 1 || m_nConnectedCompoCount == -1)
return false ;
// Se la componente richiesta non esiste, errore
if ( nPart < 0 || nPart > m_nConnectedCompoCount - 1)
return false ;
// Eseguo il calcolo della lunghezza totale degli spilloni
double dLen = 0 ;
for ( int nMap = 0 ; nMap < int( m_nMapNum) ; ++ nMap) {
for ( int nDex = 0 ; nDex < int( m_Values[nMap].size()) ; ++ nDex) {
for ( int nInt = 0 ; nInt < int( m_Values[nMap][nDex].size()) ; ++ nInt) {
// se il pezzo di spillone appartiene alla parte, ne aggiungo la lunghezza
if ( m_Values[nMap][nDex][nInt].nCompo == nPart + 1) {
dLen += m_Values[nMap][nDex][nInt].dMax - m_Values[nMap][nDex][nInt].dMin ;
}
}
}
}
// Il volume si trova moltiplicando per l'area di ogni spillone e dividendo per il numero di mappe
dVol = dLen * ( m_dStep * m_dStep) / 3 ;
return true ;
}
//----------------------------------------------------------------------------
// La retta deve essere espressa nel sistema Zmap.
// Per riferimento viene restituito un vettore di intersezioni della retta con i triangoli della superficie del solido.
// Si restituisce false in caso di errore, true altrimenti.
bool
VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTOR& vIntersInfo) const
{
// Direzione normalizzata
Vector3d vtDir = vtD ;
vtDir.Normalize() ;
// Punto e vettore espressi nel riferimento intrinseco dello Zmap
Point3d ptLocP = ptP ;
Vector3d vtDirL = vtDir ;
ptLocP.ToLoc( m_MapFrame) ;
vtDirL.ToLoc( m_MapFrame) ;
// Intersezione fra semiretta e BBox dello Zmap
double dU1, dU2 ;
bool bLineBBoxInters = IntersLineZMapBBox( ptLocP, vtDirL, dU1, dU2) && ( dU1 > 0 || dU2 > 0) ;
// Semiretta esterna al box dello Zmap quindi esterna anche allo Zmap
if ( ! bLineBBoxInters)
return true ;
// Lancio eventuale aggiornamento pendente della grafica
if ( m_nMapNum == 1)
UpdateSingleMapGraphics() ;
else
UpdateTripleMapGraphics() ;
// Ciclo sui blocchi
for ( int nB = 0 ; nB < m_nNumBlock ; ++ nB) {
// Determino indici IJK del blocco; se non trovo il blocco, errore
int nBlockIJK[3] ;
if ( ! GetBlockIJKFromN( nB, nBlockIJK))
return false ;
// Costruisco il bounding-box del blocco; se non è possibile, errore
BBox3d b3BlockBox ;
if ( ! GetBlockBox( nBlockIJK, b3BlockBox))
return false ;
// Se c'è intersezione valuto tutti i voxel interni
if ( IntersLineBox( ptLocP, vtDirL, b3BlockBox.GetMin(), b3BlockBox.GetMax())) {
// Ciclo sui voxel del blocco.
// Triangoli smooth
for ( int nV = 0 ; nV < int( m_BlockSmoothTria[nB].size()) ; ++ nV) {
// Box del voxel
BBox3d b3Vox ;
int nCurVoxIJK[3] = { m_BlockSmoothTria[nB][nV].i,
m_BlockSmoothTria[nB][nV].j,
m_BlockSmoothTria[nB][nV].k } ;
GetVoxelBox( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], b3Vox) ;
// Se non c'è intersezione col voxel, passo al successivo.
if ( ! IntersLineBox( ptLocP, vtDirL, b3Vox.GetMin(), b3Vox.GetMax()))
continue ;
for ( int nT = 0 ; nT < int( m_BlockSmoothTria[nB][nV].vTria.size()) ; ++ nT) {
Triangle3d trTria = m_BlockSmoothTria[nB][nV].vTria[nT] ;
if ( ! trTria.Validate( true))
continue ;
Point3d ptLineTria1, ptLineTria2 ;
// Studio dell'intersezione della retta con il triangolo corrente
int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, 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) {
int nNumVox ;
GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir,
nNumVox, nB, ptLineTria1, trTria) ;
}
// altrimenti ci sono due intersezioni
else {
int nNumVox ;
GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
double dP1 = ( ptLineTria1 - ptP) * vtDir ;
double dP2 = ( ptLineTria2 - ptP) * vtDir ;
vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1),
nNumVox, nB, ptLineTria1, ptLineTria2, trTria) ;
}
}
}
// Triangoli sharp interni al blocco
for ( int nV = 0 ; nV < int( m_BlockSharpTria[nB].size()) ; ++ nV) {
int nCurVoxIJK[3] = { m_BlockSharpTria[nB][nV].i,
m_BlockSharpTria[nB][nV].j,
m_BlockSharpTria[nB][nV].k } ;
// Ciclo sulle componenti connesse
for ( int nC = 0 ; nC < int( m_BlockSharpTria[nB][nV].vCompoTria.size()) ; ++ nC) {
for ( int nT = 0 ; nT < int( m_BlockSharpTria[nB][nV].vCompoTria[nC].size()) ; ++ nT) {
Triangle3d trTria = m_BlockSharpTria[nB][nV].vCompoTria[nC][nT] ;
if ( ! trTria.Validate( true))
continue ;
Point3d ptLineTria1, ptLineTria2 ;
// Studio dell'intersezione della retta con il triangolo corrente
int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, 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) {
int nNumVox ;
GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir,
nNumVox, nB, ptLineTria1, trTria) ;
}
// altrimenti ci sono due intersezioni
else {
int nNumVox ;
GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
double dP1 = ( ptLineTria1 - ptP) * vtDir ;
double dP2 = ( ptLineTria2 - ptP) * vtDir ;
vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1),
nNumVox, nB, ptLineTria1, ptLineTria2, trTria) ;
}
}
}
}
// Triangoli grandi del blocco
for ( int nT = 0 ; nT < int( m_BlockBigTria[nB].size()) ; ++ nT) {
Triangle3d trTria = m_BlockBigTria[nB][nT] ;
Point3d ptLineTria1, ptLineTria2 ;
// Studio dell'intersezione della retta con il triangolo corrente
int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, 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) {
vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir,
-1, nB, ptLineTria1, trTria) ;
}
// altrimenti ci sono due intersezioni
else {
double dP1 = ( ptLineTria1 - ptP) * vtDir ;
double dP2 = ( ptLineTria2 - ptP) * vtDir ;
vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1),
-1, nB, ptLineTria1, ptLineTria2, trTria) ;
}
}
}
// In ogni caso valuto i triangoli sharp fra blocchi
for ( int nV = 0 ; nV < int( m_InterBlockSharpTria[nB].size()) ; ++ nV) {
int nCurVoxIJK[3] = { m_InterBlockSharpTria[nB][nV].i,
m_InterBlockSharpTria[nB][nV].j,
m_InterBlockSharpTria[nB][nV].k } ;
// Ciclo sulle componenti connesse
for ( int nC = 0 ; nC < int( m_InterBlockSharpTria[nB][nV].vCompoTria.size()) ; ++ nC) {
for ( int nT = 0 ; nT < int( m_InterBlockSharpTria[nB][nV].vCompoTria[nC].size()) ; ++ nT) {
Triangle3d trTria = m_InterBlockSharpTria[nB][nV].vCompoTria[nC][nT] ;
if ( ! trTria.Validate( true))
continue ;
Point3d ptLineTria1, ptLineTria2 ;
// Studio dell'intersezione della retta con il triangolo corrente
int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, 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) {
int nNumVox ;
GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir,
nNumVox, nB, ptLineTria1, trTria) ;
}
// altrimenti ci sono due intersezioni
else {
int nNumVox ;
GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
double dP1 = ( ptLineTria1 - ptP) * vtDir ;
double dP2 = ( ptLineTria2 - ptP) * vtDir ;
vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1),
nNumVox, nB, ptLineTria1, ptLineTria2, trTria) ;
}
}
}
}
}
// Ordino le intersezioni in base al parametro distanza con segno da ptP
sort( vIntersInfo.begin(), vIntersInfo.end(),
[]( const IntLineZmapInfo& a, const IntLineZmapInfo& b)
{ double dFP = (( a.nILTT == ILTT_SEGM || a.nILTT == ILTT_SEGM_ON_EDGE) ? 0.5 * ( a.dU + a.dU2) : a.dU) ;
double dLP = (( b.nILTT == ILTT_SEGM || b.nILTT == ILTT_SEGM_ON_EDGE) ? 0.5 * ( b.dU + b.dU2) : b.dU) ;
return ( dLP > dFP) ; }) ;
return true ;
}
//----------------------------------------------------------------------------
// Il piano è espresso nel sistema locale, la funzione lo esprime nel sistema Zmap.
// I loop della flat region ottenuta dall'intersezione vengono salvati nello
// standard vector vLoop passato per riferimento.
// Se il processo va a buon fine viene restituito true, false altrimenti.
bool
VolZmap::GetPlaneIntersection( const Plane3d& plPlane, ICURVEPOVECTOR& vpLoop) const
{
// Verifico validità parametri
if ( ! plPlane.IsValid())
return false ;
// Porto il piano nel sistema Zmap
Plane3d plPlaneLoc = plPlane ;
plPlaneLoc.ToLoc( m_MapFrame) ;
// Se non c'è intersezione fra piano e bounding-box del solido, ho finito.
if ( ! TestIntersPlaneZmapBBox( plPlaneLoc))
return true ;
// Lancio eventuale aggiornamento pendente della grafica
if ( m_nMapNum == 1)
UpdateSingleMapGraphics() ;
else
UpdateTripleMapGraphics() ;
// Vettore di segmenti
vector<CurveLine> vLine ;
// Ciclo sui blocchi
for ( int nB = 0 ; nB < m_nNumBlock ; ++ nB) {
// Determino indici IJK del blocco; Se non trovo il blocco, errore
int nBlockIJK[3] ;
if ( ! GetBlockIJKFromN( nB, nBlockIJK))
return false ;
// Costruisco il bounding-bx del blocco; se non è possiblie, errore
BBox3d b3BlockBox ;
if ( ! GetBlockBox( nBlockIJK, b3BlockBox))
return false ;
// Se c'è intersezione valuto tutti i voxel interni
if ( TestIntersPlaneBox( b3BlockBox, plPlaneLoc)) {
// Ciclo sui voxel del blocco.
// Triangoli smooth
for ( int nV = 0 ; nV < int( m_BlockSmoothTria[nB].size()) ; ++ nV) {
// Box del voxel
BBox3d b3Vox ;
GetVoxelBox( m_BlockSmoothTria[nB][nV].i, m_BlockSmoothTria[nB][nV].j, m_BlockSmoothTria[nB][nV].k, b3Vox) ;
// Se non c'è intersezione col voxel, passo al successivo.
if ( ! TestIntersPlaneBox(b3Vox, plPlaneLoc))
continue ;
for ( int nT = 0 ; nT < int( m_BlockSmoothTria[nB][nV].vTria.size()) ; ++ nT) {
Triangle3d trTria = m_BlockSmoothTria[nB][nV].vTria[nT] ;
Point3d ptSt, ptEn ;
int nIntersType = IntersPlaneTria( plPlane, trTria, ptSt, ptEn) ;
if ( nIntersType == IPTT_EDGE || nIntersType == IPTT_YES) {
// Costruisco il tratto di curva
CurveLine cvLine ;
if ( cvLine.Set( ptSt, ptEn))
vLine.emplace_back( cvLine) ;
}
}
}
// Triangoli sharp interni al blocco
for ( int nV = 0 ; nV < int( m_BlockSharpTria[nB].size()) ; ++ nV) {
// Box del voxel
BBox3d b3Vox ;
GetVoxelBox( m_BlockSharpTria[nB][nV].i, m_BlockSharpTria[nB][nV].j, m_BlockSharpTria[nB][nV].k, b3Vox) ;
// Se non c'è intersezione col voxel, passo al successivo.
// Ciclo sulle componenti connesse
for ( int nC = 0 ; nC < int( m_BlockSharpTria[nB][nV].vCompoTria.size()) ; ++ nC) {
for ( int nT = 0 ; nT < int( m_BlockSharpTria[nB][nV].vCompoTria[nC].size()) ; ++ nT) {
Triangle3d trTria = m_BlockSharpTria[nB][nV].vCompoTria[nC][nT] ;
Point3d ptSt, ptEn ;
int nIntersType = IntersPlaneTria(plPlane, trTria, ptSt, ptEn) ;
if (nIntersType == IPTT_EDGE || nIntersType == IPTT_YES) {
// Costruisco il tratto di curva
CurveLine cvLine ;
if ( cvLine.Set(ptSt, ptEn))
vLine.emplace_back( cvLine) ;
}
}
}
}
// Triangoli grandi del blocco
for ( int nT = 0 ; nT < int( m_BlockBigTria[nB].size()) ; ++ nT) {
Triangle3d trTria = m_BlockBigTria[nB][nT] ;
Point3d ptSt, ptEn ;
int nIntersType = IntersPlaneTria(plPlane, trTria, ptSt, ptEn) ;
if ( nIntersType == IPTT_EDGE || nIntersType == IPTT_YES) {
// Costruisco il tratto di curva
CurveLine cvLine ;
if ( cvLine.Set(ptSt, ptEn))
vLine.emplace_back( cvLine) ;
}
}
}
// In ogni caso valuto i triangoli sharp fra blocchi
for ( int nV = 0 ; nV < int( m_InterBlockSharpTria[nB].size()) ; ++ nV) {
// Ciclo sulle componenti connesse
for ( int nC = 0 ; nC < int( m_InterBlockSharpTria[nB][nV].vCompoTria.size()) ; ++ nC) {
for ( int nT = 0 ; nT < int( m_InterBlockSharpTria[nB][nV].vCompoTria[nC].size()) ; ++ nT) {
Triangle3d trTria = m_InterBlockSharpTria[nB][nV].vCompoTria[nC][nT] ;
Point3d ptSt, ptEn ;
int nIntersType = IntersPlaneTria( plPlane, trTria, ptSt, ptEn) ;
if ( nIntersType == IPTT_EDGE || nIntersType == IPTT_YES) {
// Costruisco il tratto di curva
CurveLine cvLine ;
if ( cvLine.Set(ptSt, ptEn))
vLine.emplace_back( cvLine) ;
}
}
}
}
}
// Creo i loop
ChainCurves LoopCreator ;
LoopCreator.Init( false, EPS_SMALL, int( vLine.size())) ;
// Carico le curve per concatenarle
for ( int nCv = 0 ; nCv < int( vLine.size()) ; ++ nCv) {
Point3d ptSt = vLine[nCv].GetStart() ;
Point3d ptEn = vLine[nCv].GetEnd() ;
Vector3d vtDir; vLine[nCv].GetStartDir(vtDir) ;
LoopCreator.AddCurve( nCv + 1, ptSt, vtDir, ptEn, vtDir) ;
}
// Recupero i concatenamenti
INTVECTOR vIds ;
Point3d ptNearStart ;
while ( LoopCreator.GetChainFromNear( ptNearStart, false, vIds)) {
PtrOwner<CurveComposite> pLoop( CreateBasicCurveComposite()) ;
if ( IsNull( pLoop))
return false ;
for ( auto i : vIds) {
// Aggiungo la linea alla curva composta.
if ( ! pLoop->AddCurve( vLine[i - 1], true, 10 * EPS_SMALL))
return false ;
}
pLoop->SetExtrusion( plPlane.GetVersN()) ;
pLoop->MergeCurves( 10 * EPS_SMALL, ANG_TOL_STD_DEG) ;
// Inserisco la curva composita nella raccolta da ritornare
vpLoop.emplace_back( Release( pLoop)) ;
}
return true ;
}
//----------------------------------------------------------------------------
// Il piano deve essere espresso nel sistema di riferimento intrinseco dello Zmap.
bool
VolZmap::TestIntersPlaneZmapBBox( const Plane3d& plPlane) const
{
BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;
return TestIntersPlaneBox( b3Zmap, plPlane) ;
}
Reference in New Issue View Git Blame Copy Permalink