Dario Sassi
67cd4ef923
- correzione a creazione arco con Set2PNB quando deltaN < epsilon - dump di CurveComposite limitato a 1000 entità - correzione in Zmap a calcolo spillone con movimento ortogonale di sfera.
1716 lines
63 KiB
C++
1716 lines
63 KiB
C++
//----------------------------------------------------------------------------
|
|
// 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 "IntersLineSurfTm.h"
|
|
#include "/EgtDev/Include/EgtNumUtils.h"
|
|
#include "/EgtDev/Include/EGkIntersLineTria.h"
|
|
#include "/EgtDev/Include/EGkIntersLinePlane.h"
|
|
#include "/EgtDev/Include/EGkIntersLineSphere.h"
|
|
#include "/EgtDev/Include/EGkChainCurves.h"
|
|
|
|
using namespace std ;
|
|
|
|
//----------------------------------------------------------------------------
|
|
// 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
|
|
{
|
|
// Il box è allineato agli assi
|
|
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, 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
|
|
unsigned int nDexelPos = nJ * m_nNx[nGrid] + nI ;
|
|
unsigned int nDexelSize = unsigned 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 ( unsigned 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, 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
|
|
unsigned int nDexelPos = nJ * m_nNx[nGrid] + nI ;
|
|
unsigned int nDexelSize = unsigned 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 ( unsigned 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
|
|
{
|
|
Point3d ptP = ptPLoc ;
|
|
Vector3d vtD = vtDLoc ;
|
|
ptP.ToLoc( m_MapFrame) ;
|
|
vtD.ToLoc( m_MapFrame) ;
|
|
|
|
if ( bExact && m_nMapNum == 3)
|
|
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
|
|
{
|
|
// Intersezione fra semiretta e BBox dello Zmap
|
|
double dU1, dU2 ;
|
|
bool bTest = IntersLineZMapBBox( ptP, vtV, 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[3] ;
|
|
double dOutLen[3] ;
|
|
|
|
// Ciclo sulle griglie
|
|
for ( int nGrid = 0 ; nGrid < int( m_nMapNum) ; ++ nGrid) {
|
|
|
|
Point3d ptP0 = ptP ;
|
|
Vector3d vtV0 = vtV ;
|
|
|
|
Point3d ptI, ptF ;
|
|
// Una sola intersezione valida ( punto interno, intersezione valida 2)
|
|
if ( dU1 < 0 && dU2 > 0) {
|
|
ptI = ptP ;
|
|
ptF = ptP + dU2 * vtV ;
|
|
}
|
|
// due soluzioni valide ( punto esterno)
|
|
else {
|
|
ptI = ptP + dU1 * vtV ;
|
|
ptF = ptP + dU2 * vtV ;
|
|
}
|
|
// 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.
|
|
// 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.
|
|
bool
|
|
VolZmap::GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtD, double& dInLength, double& dOutLength, bool bEnh) const
|
|
{
|
|
// Intersezione fra semiretta e BBox dello Zmap
|
|
double dU1, dU2 ;
|
|
bool bLineBBoxInters = IntersLineZMapBBox( ptP, vtD, dU1, dU2) && ( dU1 > 0 || dU2 > 0) ;
|
|
|
|
// Semiretta esterna al box dello Zmap quindi esterna anche allo Zmap
|
|
if ( ! bLineBBoxInters) {
|
|
dInLength = - 2 ;
|
|
dOutLength = - 2 ;
|
|
return true ;
|
|
}
|
|
|
|
// Determino il voxel di partenza
|
|
int nVoxI, nVoxJ, nVoxK ;
|
|
if ( ! GetPointVoxel( ptP, nVoxI, nVoxJ, nVoxK)) {
|
|
if ( ! GetPointVoxel( ptP + dU1 * vtD, nVoxI, nVoxJ, nVoxK))
|
|
return false ;
|
|
}
|
|
|
|
// 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
|
|
vector<VoxelIndexes> vVox ;
|
|
while ( IsValidVoxel( nVoxI, nVoxJ, nVoxK)) {
|
|
// Estremi della diagonale del voxel corrente
|
|
Point3d ptMin( ( nVoxI * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nVoxJ * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nVoxK * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
Point3d ptMax( ( ( nVoxI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( nVoxJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( nVoxK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
// Studio il voxel corrente
|
|
if ( IntersLineBox( ptP, vtD, ptMin, ptMax)) {
|
|
VoxelIndexes 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) * N_DEXVOXRATIO + 0.5) * m_dStep - ptP.x) / vtD.x) : INFINITO) ;
|
|
double dMaxTY = ( abs( vtD.y) > EPS_ZERO ?
|
|
abs( ( ( ( nVoxJ + nPlaneJ) * N_DEXVOXRATIO + 0.5) * m_dStep - ptP.y) / vtD.y) : INFINITO) ;
|
|
double dMaxTZ = ( abs( vtD.z) > EPS_ZERO ?
|
|
abs( ( ( ( nVoxK + nPlaneK) * N_DEXVOXRATIO + 0.5) * m_dStep - ptP.z) / vtD.z) : INFINITO) ;
|
|
// Determino gli incrementi
|
|
if ( dMaxTX < dMaxTY) {
|
|
if ( dMaxTX < dMaxTZ)
|
|
nVoxI += nDeltaI ;
|
|
else
|
|
nVoxK += nDeltaK ;
|
|
}
|
|
else {
|
|
if ( dMaxTY < dMaxTZ)
|
|
nVoxJ += nDeltaJ ;
|
|
else
|
|
nVoxK += nDeltaK ;
|
|
}
|
|
}
|
|
|
|
// Triangoli di frontiera dei voxel
|
|
TRIA3DEXLIST lstTria ;
|
|
ExtMarchingCubes( vVox, lstTria, bEnh) ;
|
|
|
|
// Dati dell'intersezione
|
|
struct LineTriaInt {
|
|
int nNum ;
|
|
double dPar1 ;
|
|
double dPar2 ;
|
|
double dDot ;
|
|
LineTriaInt( void) : nNum( 0) {}
|
|
LineTriaInt( double dP, double dD) : nNum( 1), dPar1( dP), dDot( dD) {}
|
|
LineTriaInt( double dP1, double dP2, double dD)
|
|
: nNum( 2), dPar1( dP1), dPar2( dP2), dDot( dD) {}
|
|
} ;
|
|
vector<LineTriaInt> vInt ;
|
|
|
|
// Ciclo sui triangoli dei voxel
|
|
for ( auto it = lstTria.begin() ; it != lstTria.end() ; ++it ) {
|
|
// Triangolo corrente e suoi punti di intersezione con la retta
|
|
const 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) {
|
|
vInt.emplace_back( ( ptLineTria1 - ptP) * vtD, vtD * CurrTria.GetN()) ;
|
|
}
|
|
// altrimenti ci sono due intersezioni
|
|
else {
|
|
double dP1 = ( ptLineTria1 - ptP) * vtD ;
|
|
double dP2 = ( ptLineTria2 - ptP) * vtD ;
|
|
double dD = vtD * CurrTria.GetN() ;
|
|
vInt.emplace_back( ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1), dD) ;
|
|
}
|
|
}
|
|
|
|
// Ordino le intersezioni in base al parametro distanza con segno da ptP
|
|
sort( vInt.begin(), vInt.end(),
|
|
[]( const LineTriaInt& a, const LineTriaInt& b)
|
|
{ double dFP = ( a.nNum == 2 ? 0.5 * ( a.dPar1 + a.dPar2) : a.dPar1) ;
|
|
double dLP = ( b.nNum == 2 ? 0.5 * ( b.dPar1 + b.dPar2) : b.dPar1) ;
|
|
return ( dLP > dFP) ; }) ;
|
|
|
|
// Inizializzo le distanze di ingresso e uscita:
|
|
// dInLength diminuisce, dOutLength aumenta.
|
|
dInLength = INFINITO ;
|
|
dOutLength = - INFINITO ;
|
|
|
|
int nFirstPosN ;
|
|
int nN = 0 ;
|
|
for ( ; nN < int( vInt.size()) ; ++ nN) {
|
|
if ( vInt[nN].dPar1 > - EPS_SMALL) {
|
|
nFirstPosN = nN ;
|
|
break ;
|
|
}
|
|
}
|
|
|
|
if ( nN == int( vInt.size())) {
|
|
dInLength = - 2 ;
|
|
dOutLength = - 2 ;
|
|
return true ;
|
|
}
|
|
|
|
if ( nFirstPosN > 0) {
|
|
if ( vInt[nFirstPosN - 1].dDot < EPS_ZERO)
|
|
dInLength = -1 ;
|
|
}
|
|
else if ( nFirstPosN == 0) {
|
|
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) {
|
|
if ( vInt[nN].dDot > - EPS_ZERO) {
|
|
if ( vInt[nN].nNum == 2 && dOutLength < vInt[nN].dPar2)
|
|
dOutLength = vInt[nN].dPar2 ;
|
|
else if ( dOutLength < vInt[nN].dPar1)
|
|
dOutLength = vInt[nN].dPar1 ;
|
|
}
|
|
if ( vInt[nN].dDot < EPS_ZERO &&
|
|
dInLength > vInt[nN].dPar1)
|
|
dInLength = vInt[nN].dPar1 ;
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::AvoidBox( const Frame3d& frBox, const Vector3d& vtDiag) const
|
|
{
|
|
// BBox
|
|
BBox3d b3Box( ORIG, ORIG + vtDiag) ;
|
|
|
|
// lo porto nel riferimento intrinseco dello Zmap
|
|
b3Box.LocToLoc( frBox, 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, frBox) ;
|
|
|
|
// Riferimento intrinseco dei dexel nel riferimento del box
|
|
Point3d ptO = ORIG ; ptO.LocToLoc( m_MapFrame, frBox) ;
|
|
Vector3d vtX = X_AX ; vtX.LocToLoc( m_MapFrame, frBox) ;
|
|
Vector3d vtY = Y_AX ; vtY.LocToLoc( m_MapFrame, frBox) ;
|
|
|
|
// 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 + vtDiag, 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) const
|
|
{
|
|
// Porto la sfera nel riferimento intrinseco dello Zmap
|
|
Point3d ptC = ptCenter ;
|
|
ptC.ToLoc( m_MapFrame) ;
|
|
|
|
// 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 = ptI1.z ;
|
|
double dZmax = 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 dL, double dR) const
|
|
{
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// NB: L'origine del sistema di riferimento è nel centro della circonferenza di base
|
|
// e l'asse di simmetria coincide 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 dR, bool bTapB, bool bTapT,
|
|
Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2)
|
|
{
|
|
Point3d ptP = ptLineSt ;
|
|
Vector3d vtV = vtLineDir ;
|
|
|
|
// Trasformazione delle coordinate:
|
|
// l'asse del cilindro corrisponde con
|
|
// l'asse z del sistema di riferimento
|
|
ptP.ToLoc( CylFrame) ;
|
|
vtV.ToLoc( CylFrame) ;
|
|
|
|
DBLVECTOR vdCoef(3) ;
|
|
DBLVECTOR vdRoots ;
|
|
// Non vogliamo che i dexel a filo vengano tagliati
|
|
double dSqRad = dR * dR ;
|
|
double dSqRadSafe = dSqRad - 2 * dR * EPS_SMALL ;
|
|
|
|
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 ;
|
|
|
|
// Computo radici
|
|
int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ;
|
|
|
|
// Nessuna soluzione
|
|
if ( nRoot == 0 || nRoot == 1) {
|
|
if ( abs( vtV.z) > EPS_ZERO) {
|
|
// Intersezioni con i piani che limitano il cilindro in altezza
|
|
ptInt1 = ptP - ( ptP.z / vtV.z) * vtV ;
|
|
ptInt2 = ptP + ( ( dH - ptP.z) / vtV.z) * vtV ;
|
|
// Normali nei punti di interseione
|
|
vtN1 = Z_AX ;
|
|
vtN2 = - Z_AX ;
|
|
// Se le soluzioni sono all'interno delle circonferenze
|
|
if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y <= dSqRadSafe &&
|
|
ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y <= dSqRadSafe) {
|
|
// Trasformiamo le coordinate nel sistema Zmap e abbiamo finito
|
|
ptInt1.ToGlob( CylFrame) ;
|
|
ptInt2.ToGlob( CylFrame) ;
|
|
vtN1.ToGlob( CylFrame) ;
|
|
vtN2.ToGlob( CylFrame) ;
|
|
return true ;
|
|
}
|
|
// Nessuna intersezione
|
|
else
|
|
return false ;
|
|
}
|
|
// Nessuna intersezione
|
|
else
|
|
return false ;
|
|
}
|
|
|
|
// L'equazione ammette o due soluzioni (eventualmente
|
|
// coincidenti) oppure nessuna o infinite se la la retta
|
|
// appartiene alla superficie
|
|
|
|
if ( nRoot == 2) {
|
|
// Ordino i parametri di intersezione
|
|
double dUmin = vdRoots[0] ;
|
|
double dUmax = vdRoots[1] ;
|
|
if ( dUmin > dUmax)
|
|
swap( dUmin, dUmax) ;
|
|
// Calcolo i punti d'intersezione (ordinati secondo Z crescente)
|
|
ptInt1 = ptP + vdRoots[0] * vtV ;
|
|
ptInt2 = ptP + vdRoots[1] * vtV ;
|
|
if ( ptInt1.z > ptInt2.z)
|
|
swap( ptInt1, ptInt2) ;
|
|
// Quote limitazione, dipendenti dalla tappatura estremità
|
|
double dZbot = ( bTapB ? -EPS_SMALL : EPS_SMALL) ;
|
|
double dZtop = ( bTapT ? dH + EPS_SMALL : dH - EPS_SMALL) ;
|
|
// Se intersezioni entrambe fuori dal cilindro limitato, non vanno considerate
|
|
if ( ptInt2.z < dZbot || ptInt1.z > dZtop)
|
|
return false ;
|
|
// Calcolo le normali
|
|
vtN1.Set( ( ORIG - ptInt1).x, ( ORIG - ptInt1).y, 0) ;
|
|
vtN1.Normalize() ;
|
|
vtN2.Set( ( ORIG - ptInt2).x, ( ORIG - ptInt2).y, 0) ;
|
|
vtN2.Normalize() ;
|
|
// Limitazioni per intersezione con piano basso
|
|
if ( ptInt1.z < dZbot - EPS_ZERO) {
|
|
ptInt1 = ptP + Clamp( ( dZbot - ptP.z / vtV.z), dUmin, dUmax) * vtV ;
|
|
vtN1.Set( 0, 0, 1) ;
|
|
}
|
|
// Limitazioni per intersezione con piano alto
|
|
if ( ptInt2.z > dZtop + EPS_ZERO) {
|
|
ptInt2 = ptP + Clamp( ( ( dZtop - ptP.z) / vtV.z), dUmin, dUmax) * vtV ;
|
|
vtN2.Set( 0, 0, -1) ;
|
|
}
|
|
// Riporto le coordinate nel sistema di riferimento griglia
|
|
ptInt1.ToGlob( CylFrame) ;
|
|
ptInt2.ToGlob( CylFrame) ;
|
|
vtN1.ToGlob( CylFrame) ;
|
|
vtN2.ToGlob( CylFrame) ;
|
|
}
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IntersLineConus( const Point3d& ptLineSt, const Vector3d& vtLineDir,
|
|
const Frame3d& ConusFrame, double dTan, double dh, double dH, bool bTapLow, bool bTapUp,
|
|
Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2)
|
|
{
|
|
// NB: L'origine del sistema di riferimento deve essere
|
|
// nel vertice del cono e l'asse di simmetria deve coincidere con l'asse x.
|
|
// La funzione restituisce true in caso di intersezione, false altrimenti.
|
|
|
|
Point3d ptP = ptLineSt ;
|
|
Vector3d vtV = vtLineDir ;
|
|
|
|
// Trasformazione delle coordinate
|
|
ptP.ToLoc( ConusFrame) ;
|
|
vtV.ToLoc( ConusFrame) ;
|
|
|
|
DBLVECTOR vdCoef(3) ;
|
|
DBLVECTOR vdRoots ;
|
|
|
|
double dSqTan = dTan * dTan ;
|
|
double dMinRad = dTan * dh ;
|
|
double dMaxRad = dTan * dH ;
|
|
double dDeltaR = dMaxRad - dMinRad ;
|
|
double dHei = dH - dh ;
|
|
|
|
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 ;
|
|
|
|
// Computo radici
|
|
int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ;
|
|
|
|
// Nessuna soluzione
|
|
if ( nRoot == 0)
|
|
return false ;
|
|
|
|
double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ;
|
|
double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ;
|
|
|
|
// Una soluzione: la retta iterseca superficie
|
|
// laterale e un piano
|
|
if ( nRoot == 1) {
|
|
|
|
ptInt1 = ptP + vdRoots[0] * vtV ;
|
|
|
|
Vector3d vtU = ( ptInt1 - ORIG) - ( ptInt1 - ORIG).z * Z_AX ;
|
|
|
|
vtU.Normalize() ;
|
|
|
|
vtN1 = dDeltaR * Z_AX - dHei * vtU ;
|
|
|
|
vtN1.Normalize() ;
|
|
|
|
if ( ptInt1.z < dH + dEpsUp) {
|
|
|
|
if ( ptInt1.z > dh + dEpsLow) {
|
|
|
|
ptInt2 = ptP + ( ( dH - ptP.z) / vtV.z) * vtV ;
|
|
|
|
vtN2 = - Z_AX ;
|
|
|
|
}
|
|
else if ( ptInt1.z > - EPS_SMALL) {
|
|
|
|
ptInt1 = ptP + ( ( dh - ptP.z) / vtV.z) * vtV ;
|
|
ptInt2 = ptP + ( ( dH - ptP.z) / vtV.z) * vtV ;
|
|
|
|
vtN1 = Z_AX ;
|
|
vtN2 = - Z_AX ;
|
|
|
|
if ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y > dMaxRad * dMaxRad)
|
|
|
|
return false ;
|
|
}
|
|
else
|
|
return false ;
|
|
|
|
ptInt1.ToGlob( ConusFrame) ;
|
|
ptInt2.ToGlob( ConusFrame) ;
|
|
|
|
vtN1.ToGlob( ConusFrame) ;
|
|
vtN2.ToGlob( ConusFrame) ;
|
|
|
|
return true ;
|
|
}
|
|
else
|
|
return false ;
|
|
|
|
}
|
|
// Due soluzioni: la retta interseca due volte la
|
|
// superficie laterale
|
|
else if ( nRoot == 2) {
|
|
|
|
ptInt1 = ptP + vdRoots[0] * vtV ;
|
|
ptInt2 = ptP + vdRoots[1] * vtV ;
|
|
|
|
if ( ptInt1.z > ptInt2.z)
|
|
swap( ptInt1, ptInt2) ;
|
|
|
|
Vector3d vtU1 = ( ptInt1 - ORIG) - ( ptInt1 - ORIG).z * Z_AX ;
|
|
Vector3d vtU2 = ( ptInt2 - ORIG) - ( ptInt2 - ORIG).z * Z_AX ;
|
|
|
|
vtU1.Normalize() ;
|
|
vtU2.Normalize() ;
|
|
|
|
vtN1 = dDeltaR * Z_AX - dHei * vtU1 ;
|
|
vtN2 = dDeltaR * Z_AX - dHei * vtU2 ;
|
|
|
|
vtN1.Normalize() ;
|
|
vtN2.Normalize() ;
|
|
|
|
if ( abs( vtV.z) < EPS_ZERO) {
|
|
|
|
if ( ptInt1.z > dh + dEpsLow && ptInt1.z < dH + dEpsUp) {
|
|
|
|
ptInt1.ToGlob( ConusFrame) ;
|
|
ptInt2.ToGlob( ConusFrame) ;
|
|
|
|
vtN1.ToGlob( ConusFrame) ;
|
|
vtN2.ToGlob( ConusFrame) ;
|
|
|
|
vtN1.Normalize() ;
|
|
vtN2.Normalize() ;
|
|
|
|
return true ;
|
|
}
|
|
else
|
|
return false ;
|
|
}
|
|
|
|
|
|
if ( ptInt1.z < dH + dEpsUp) {
|
|
|
|
if ( ptInt1.z > dh + dEpsLow) {
|
|
|
|
if ( ptInt2.z > dH + dEpsUp) {
|
|
|
|
ptInt2 = ptP + ( ( dH - ptP.z) / vtV.z) * vtV ;
|
|
|
|
vtN2 = - Z_AX ;
|
|
}
|
|
}
|
|
else if ( ptInt1.z > - EPS_SMALL) {
|
|
|
|
if ( ptInt2.z > dH + dEpsUp) {
|
|
|
|
ptInt1 = ptP + ( ( dh - ptP.z) / vtV.z) * vtV ;
|
|
ptInt2 = ptP + ( ( dH - ptP.z) / vtV.z) * vtV ;
|
|
|
|
vtN1 = Z_AX ;
|
|
vtN2 = - Z_AX ;
|
|
}
|
|
else if ( ptInt2.z > dh + dEpsLow) {
|
|
|
|
ptInt1 = ptP + ( ( dh - ptP.z) / vtV.z) * vtV ;
|
|
vtN1 = Z_AX ;
|
|
}
|
|
else
|
|
return false ;
|
|
}
|
|
else {
|
|
|
|
if ( ptInt2.z < 0)
|
|
|
|
return false ;
|
|
|
|
else if ( ptInt2.z < dh + dEpsLow) {
|
|
|
|
ptInt1 = ptP + ( ( dh - ptP.z) / vtV.z) * vtV ;
|
|
ptInt2 = ptP + ( ( dH - ptP.z) / vtV.z) * vtV ;
|
|
|
|
vtN1 = Z_AX ;
|
|
vtN2 = - Z_AX ;
|
|
}
|
|
else if ( ptInt2.z < dH + dEpsUp) {
|
|
|
|
ptInt1 = ptP + ( ( dH - ptP.z) / vtV.z) * vtV ;
|
|
vtN1 = - Z_AX ;
|
|
}
|
|
else
|
|
return false ;
|
|
}
|
|
|
|
ptInt1.ToGlob( ConusFrame) ;
|
|
ptInt2.ToGlob( ConusFrame) ;
|
|
|
|
vtN1.ToGlob( ConusFrame) ;
|
|
vtN2.ToGlob( ConusFrame) ;
|
|
|
|
vtN1.Normalize() ;
|
|
vtN2.Normalize() ;
|
|
|
|
return true ;
|
|
}
|
|
else
|
|
return false ;
|
|
}
|
|
return false ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// NB: 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 Vector3d& vtLineDir, const Point3d& ptLineSt,
|
|
const Frame3d& CircFrame, double dRad, double dLongMvLen, double dOrtMvLen,
|
|
bool bTapLow, bool bTapUp,
|
|
Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2)
|
|
{
|
|
// Quadrato del raggio
|
|
double dSqRad = dRad * dRad ;
|
|
// Punto e vettore individuanti la retta
|
|
Point3d ptP = ptLineSt ;
|
|
Vector3d vtV = vtLineDir ;
|
|
// Asse cilindro ellittico
|
|
Vector3d vtAx( dOrtMvLen, 0, dLongMvLen) ;
|
|
vtAx.Normalize() ;
|
|
// Se il cilindro ellittico degenera in una superficie,
|
|
// non bisogna tagliare
|
|
if ( abs( vtAx.z) < EPS_SMALL)
|
|
return false ;
|
|
// Trasformazione delle coordinate
|
|
ptP.ToLoc( CircFrame) ;
|
|
vtV.ToLoc( CircFrame) ;
|
|
// 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 ;
|
|
}
|
|
|
|
vector <double> vdCoef(3) ;
|
|
vector <double> vdRoots ;
|
|
// 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
|
|
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
|
|
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 ci sono intersezioni con i tappi o l'equazione
|
|
// degenera in una di primo grado, le eventuali
|
|
// soluzioni sono già state trovate.
|
|
if ( nRoot == 0 || nRoot == 1)
|
|
return false ;
|
|
// Due soluzioni trovate
|
|
else if ( nRoot == 2) {
|
|
// Flag per i tappi
|
|
double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ;
|
|
double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ;
|
|
// 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
|
|
Vector3d vtTest1 = ( ptInt1 - ORIG) - ( ptInt1 - ORIG) * vtAx * vtAx ;
|
|
Vector3d vtTest2 = ( ptInt2 - ORIG) - ( ptInt2 - ORIG) * vtAx * vtAx ;
|
|
double dX0_1, dX0_2 ;
|
|
if ( vtTest1.x > 0) {
|
|
dX0_1 = ( dSqRad - ptInt1.y * ptInt1.y > 0 ? sqrt( dSqRad - ptInt1.y * ptInt1.y) : 0) ;
|
|
}
|
|
else {
|
|
dX0_1 = ( dSqRad - ptInt1.y * ptInt1.y > 0 ? - sqrt( dSqRad - ptInt1.y * ptInt1.y) : 0) ;
|
|
}
|
|
Vector3d vtCirc1( - dX0_1, - ptInt1.y, 0) ;
|
|
Vector3d vtTan1( vtCirc1.y, - vtCirc1.x, 0) ;
|
|
Vector3d vtCross1 = vtTan1 ^ vtMv ;
|
|
// Vettore 1
|
|
vtN1 = ( vtCross1 * vtCirc1 > - EPS_ZERO ? vtCross1 : - vtCross1) ;
|
|
if ( vtTest2.x > 0) {
|
|
dX0_2 = ( dSqRad - ptInt2.y * ptInt2.y > 0 ? sqrt( dSqRad - ptInt2.y * ptInt2.y) : 0) ;
|
|
}
|
|
else {
|
|
dX0_2 = ( dSqRad - ptInt2.y * ptInt2.y > 0 ? - sqrt( dSqRad - ptInt2.y * ptInt2.y) : 0) ;
|
|
}
|
|
Vector3d vtCirc2( - dX0_2, - ptInt2.y, 0) ;
|
|
Vector3d vtTan2( vtCirc2.y, - vtCirc2.x, 0) ;
|
|
Vector3d vtCross2 = vtTan2 ^ vtMv ;
|
|
// Vettore 2
|
|
vtN2 = ( vtCross2 * vtCirc2 > - EPS_ZERO ? vtCross2 : - vtCross2) ;
|
|
// Normalizzo i vettori
|
|
vtN1.Normalize() ;
|
|
vtN2.Normalize() ;
|
|
|
|
// Studio le soluzioni: se ua è 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 ;
|
|
|
|
// Riporto le coordinate nel sistema
|
|
// di riferimento griglia.
|
|
ptInt1.ToGlob( CircFrame) ;
|
|
ptInt2.ToGlob( CircFrame) ;
|
|
vtN1.ToGlob( CircFrame) ;
|
|
vtN2.ToGlob( CircFrame) ;
|
|
return true ;
|
|
}
|
|
return false ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir,
|
|
const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaZ,
|
|
Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2)
|
|
{
|
|
Point3d ptP = ptLineSt ;
|
|
Vector3d vtV = vtLineDir ;
|
|
|
|
// Trasformazione delle coordinate
|
|
ptP.ToLoc( PolyFrame) ;
|
|
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) ;
|
|
|
|
// Servono per descrivere i piani obliqui
|
|
Vector3d vtFacet5 = ptFacet135 - ptP ;
|
|
Vector3d vtFacet6 = ptFacet246 - ptP ;
|
|
|
|
Vector3d vtOb( - dDeltaZ, 0, dLenX) ;
|
|
|
|
vtOb.Normalize() ;
|
|
|
|
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 ;
|
|
|
|
// 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) > dLenY / 2 - EPS_SMALL)
|
|
return false ;
|
|
// Controllo sulle facce 3 e 4
|
|
if ( abs( vtV.x) < EPS_ZERO &&
|
|
( ptP.x < EPS_SMALL ||
|
|
ptP.x > dLenX - EPS_SMALL))
|
|
return false ;
|
|
// Controllo sulle facce 5 e 6
|
|
/*Vector3d vtW( 0, dLenX, dDeltaZ) ;
|
|
vtW.Normalize() ;
|
|
Vector3d vtU = vtV - vtV.y * Y_AX - vtV * vtW * vtW ;
|
|
if ( vtU.Len() < EPS_ZERO &&
|
|
( ptP.z * dLenX < dDeltaZ * ptP.x + dLenX * EPS_SMALL ||
|
|
ptP.z * dLenX > dDeltaZ * ptP.x + dLenX * ( dLenY - EPS_SMALL)))
|
|
return false ;*/
|
|
double dDotObV = abs( vtV * vtOb) ;
|
|
Vector3d vtP1 = ptFacet135 - ptP ;
|
|
Vector3d vtP2 = ptFacet246 - ptP ;
|
|
double dP1 = abs ( vtP1 * vtOb) ;
|
|
double dP2 = abs ( vtP2 * vtOb) ;
|
|
if ( dDotObV < EPS_ZERO &&
|
|
( dP1 < EPS_SMALL || dP2 < EPS_SMALL))
|
|
return false ;
|
|
|
|
// 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) {
|
|
ptInt1.ToGlob( PolyFrame) ;
|
|
ptInt2.ToGlob( PolyFrame) ;
|
|
vtN1.ToGlob( PolyFrame) ;
|
|
vtN2.ToGlob( PolyFrame) ;
|
|
return true ;
|
|
}
|
|
else
|
|
return false ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
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 ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// 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 plPlaneL = plPlane ;
|
|
plPlaneL.ToLoc( m_MapFrame) ;
|
|
// Vettore del piano
|
|
Vector3d vtNL = plPlaneL.GetVersN() ;
|
|
|
|
// Se il piano non interseca il bounding box del solido, abbiamo finito
|
|
if ( ! TestIntersPlaneZmapBBox( plPlaneL))
|
|
return true ;
|
|
|
|
// Cerco la mappa con i dexel più perpendicolari al piano e valuto
|
|
// se il versore del piano è equiverso o controverso ai dexel.
|
|
int nGrid = 0 ;
|
|
int nSign = ( vtNL.z > 0 ? 1 : - 1) ;
|
|
if ( abs( vtNL.x) > abs( vtNL.y) && abs( vtNL.x) > abs( vtNL.z)) {
|
|
nGrid = 1 ;
|
|
nSign = ( vtNL.x > 0 ? 1 : - 1) ;
|
|
}
|
|
else if ( abs( vtNL.y) > abs( vtNL.x) && abs( vtNL.y) > abs( vtNL.z)) {
|
|
nGrid = 2 ;
|
|
nSign = ( vtNL.y > 0 ? 1 : - 1) ;
|
|
}
|
|
|
|
// Ciclo sulle celle
|
|
vector<CurveLine> vLine ;
|
|
for ( int ni = - 1 ; ni < int( m_nNx[nGrid]) ; ++ ni) {
|
|
for ( int nj = - 1 ; nj < int( m_nNy[nGrid]) ; ++ nj) {
|
|
ProcessCell( nGrid, ni, nj, plPlaneL, vLine) ;
|
|
}
|
|
}
|
|
|
|
// 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, EPS_SMALL))
|
|
return false ;
|
|
}
|
|
pLoop->SetExtrusion( vtNL) ;
|
|
pLoop->ToGlob( m_MapFrame) ;
|
|
// Se normali controverse, devo invertire il verso del loop.
|
|
if ( nSign < 0)
|
|
pLoop->Invert() ;
|
|
// eseguo approssimazione
|
|
PolyLine PL ;
|
|
if ( ! pLoop->ApproxWithLines( m_dStep, ANG_TOL_STD_DEG, ICurve::APL_RIGHT, PL) ||
|
|
! pLoop->Clear() || ! pLoop->FromPolyLine( PL))
|
|
return false ;
|
|
pLoop->MergeCurves( m_dStep, ANG_TOL_STD_DEG) ;
|
|
// Inserisco la curva composita nella raccolta da ritornare
|
|
vpLoop.emplace_back( Release( pLoop)) ;
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Data una cella, identificata dal numero della sua griglia e dai suoi indici i,j
|
|
// all'interno della griglia, genera gli eventuali segmenti che costituiscono la curva
|
|
// corrispondente alla frontiera della regione di piano interna la solido.
|
|
// Se la griglia o la cella non è valida, la funzione restituisce false,
|
|
// altrimenti valuta la tipologia della cella e crea gli eventuali segmenti.
|
|
// Se il processo va a buon fine la funzione restituisce true.
|
|
bool
|
|
VolZmap::ProcessCell( int nGrid, int nCellI, int nCellJ, const Plane3d& plPlane, vector<CurveLine>& vLine) const
|
|
{
|
|
// Se la griglia non esiste vi è un errore
|
|
if ( nGrid < 0 || nGrid > 2)
|
|
return false ;
|
|
// Se la cella non esiste vi è un errore
|
|
if ( nCellI < - 1 || nCellI >= int( m_nNx[nGrid]) ||
|
|
nCellJ < - 1 || nCellJ >= int( m_nNy[nGrid]))
|
|
return false ;
|
|
|
|
// Determino la configurazione della cella
|
|
int nIndex = CalcIndexForPlaneCells( plPlane, nGrid, nCellI, nCellJ) ;
|
|
|
|
// Tabella segmenti
|
|
static int nLineTable[16][5] = {
|
|
{ -1, -1, -1, -1, -1},
|
|
{ 0, 3, -1, -1, -1},
|
|
{ 1, 0, -1, -1, -1},
|
|
{ 1, 3, -1, -1, -1},
|
|
{ 2, 1, -1, -1, -1},
|
|
{ 0, 1, 2, 3, -1},
|
|
{ 2, 0, -1, -1, -1},
|
|
{ 2, 3, -1, -1, -1},
|
|
{ 3, 2, -1, -1, -1},
|
|
{ 0, 2, -1, -1, -1},
|
|
{ 1, 2, 3, 0, -1},
|
|
{ 1, 2, -1, -1, -1},
|
|
{ 3, 1, -1, -1, -1},
|
|
{ 0, 1, -1, -1, -1},
|
|
{ 3, 0, -1, -1, -1},
|
|
{ -1, -1, -1, -1, -1}
|
|
} ;
|
|
// Tabella dei punti medi dei segmenti: le righe rappresentano
|
|
// il numero di segmento (0, 1, 2, 3) le colonne rappresentano
|
|
// i e j.
|
|
static double dDeltaIJ[4][2] = {
|
|
{ 0.5, 0},
|
|
{ 1, 0.5},
|
|
{ 0.5, 1},
|
|
{ 0, 0.5}
|
|
} ;
|
|
// Se la cella è di frontiera costruisco i segmenti della curva
|
|
if ( nIndex != 0 && nIndex != 15) {
|
|
// Ciclo su tutti gli spigoli della cella che vengono attraversati dalla curva.
|
|
for ( int nEdge = 0 ; nLineTable[nIndex][nEdge] != -1 ; nEdge += 2) {
|
|
|
|
int nStEdge = nLineTable[nIndex][nEdge] ;
|
|
int nEnEdge = nLineTable[nIndex][nEdge + 1] ;
|
|
double nStDeltaI = dDeltaIJ[nStEdge][0] ;
|
|
double nStDeltaJ = dDeltaIJ[nStEdge][1] ;
|
|
double nEnDeltaI = dDeltaIJ[nEnEdge][0] ;
|
|
double nEnDeltaJ = dDeltaIJ[nEnEdge][1] ;
|
|
|
|
if ( nGrid == 0) {
|
|
// Punti sulla griglia
|
|
Point3d ptGrSt( ( nCellI + nStDeltaI + 0.5) * m_dStep, ( nCellJ + nStDeltaJ + 0.5) * m_dStep, 0) ;
|
|
Point3d ptGrEn( ( nCellI + nEnDeltaI + 0.5) * m_dStep, ( nCellJ + nEnDeltaJ + 0.5) * m_dStep, 0) ;
|
|
// Corrispondenti punti sul piano
|
|
Point3d ptSt, ptEn ;
|
|
if ( IntersLinePlane( ptGrSt, Z_AX, 10, plPlane, ptSt, false) &&
|
|
IntersLinePlane( ptGrEn, Z_AX, 10, plPlane, ptEn, false)) {
|
|
// Costruisco il tratto di curva
|
|
CurveLine cvLine ;
|
|
cvLine.Set( ptSt, ptEn) ;
|
|
vLine.emplace_back( cvLine) ;
|
|
}
|
|
}
|
|
|
|
else if ( nGrid == 1) {
|
|
// Punti sulla griglia
|
|
Point3d ptGrSt( 0, ( nCellI + nStDeltaI + 0.5) * m_dStep, ( nCellJ + nStDeltaJ + 0.5) * m_dStep) ;
|
|
Point3d ptGrEn( 0, ( nCellI + nEnDeltaI + 0.5) * m_dStep, ( nCellJ + nEnDeltaJ + 0.5) * m_dStep) ;
|
|
// Corrispondenti punti sul piano
|
|
Point3d ptSt, ptEn ;
|
|
if ( IntersLinePlane( ptGrSt, X_AX, 10, plPlane, ptSt, false) &&
|
|
IntersLinePlane( ptGrEn, X_AX, 10, plPlane, ptEn, false)) {
|
|
// Costruisco il tratto di curva
|
|
CurveLine cvLine ;
|
|
cvLine.Set( ptSt, ptEn) ;
|
|
vLine.emplace_back( cvLine) ;
|
|
}
|
|
}
|
|
|
|
else {
|
|
// Punti sulla griglia
|
|
Point3d ptGrSt( ( nCellJ + nStDeltaJ + 0.5) * m_dStep, 0, ( nCellI + nStDeltaI + 0.5) * m_dStep) ;
|
|
Point3d ptGrEn( ( nCellJ + nEnDeltaJ + 0.5) * m_dStep, 0, ( nCellI + nEnDeltaI + 0.5) * m_dStep) ;
|
|
// Corrispondenti punti sul piano
|
|
Point3d ptSt, ptEn ;
|
|
if ( IntersLinePlane( ptGrSt, Y_AX, 10, plPlane, ptSt, false) &&
|
|
IntersLinePlane( ptGrEn, Y_AX, 10, plPlane, ptEn, false)) {
|
|
// Costruisco il tratto di curva
|
|
CurveLine cvLine ;
|
|
cvLine.Set( ptSt, ptEn) ;
|
|
vLine.emplace_back( cvLine) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Dato un dexel, identificato dal numero della sua griglia e dai suoi indici,
|
|
// determina se un suo intervallo attraversa il piano dato.
|
|
// Il numero di griglia e gli indici del dexel devono essere validi, se
|
|
// non lo sono, la funzione restituisce false, altrimenti cerca un
|
|
// tratto di dexel che intersechi il piano, se lo trova restituisce true,
|
|
// false altrimenti.
|
|
bool
|
|
VolZmap::InOut( const Plane3d& plPlane, int nGrid, int nI, int nJ) const
|
|
{
|
|
// Se la griglia non esiste, vi è un errore.
|
|
if ( nGrid < 0 || nGrid > 2)
|
|
return false ;
|
|
// Se gli indici sono di frontiera per lo
|
|
// Zmap o non esistono, non sono interni.
|
|
if ( nI <= - 1 || nI >= int( m_nNx[nGrid]) ||
|
|
nJ <= - 1 || nJ >= int( m_nNy[nGrid]))
|
|
return false ;
|
|
// Numero del dexel
|
|
int nDex = nJ * m_nNx[nGrid] + nI ;
|
|
|
|
// Ciclo sui segmenti del dexel
|
|
bool bNotFound = true ;
|
|
for ( int nK = 0 ; nK < int( m_Values[nGrid][nDex].size()) ; ++ nK) {
|
|
if ( nGrid == 0) {
|
|
// Punti estremi del segmento
|
|
Point3d ptSt( ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMin) ;
|
|
Point3d ptEn( ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMax) ;
|
|
// Se il segmento interseca il piano abbiamo finito
|
|
Point3d ptInt ;
|
|
if ( IntersLinePlane( ptSt, ptEn, plPlane, ptInt))
|
|
return true ;
|
|
}
|
|
else if ( nGrid == 1) {
|
|
// Punti estremi del segmento
|
|
Point3d ptSt( m_Values[nGrid][nDex][nK].dMin, ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep) ;
|
|
Point3d ptEn( m_Values[nGrid][nDex][nK].dMax, ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep) ;
|
|
// Se il segmento interseca il piano abbiamo finito
|
|
Point3d ptInt ;
|
|
if ( IntersLinePlane( ptSt, ptEn, plPlane, ptInt))
|
|
return true ;
|
|
}
|
|
else {
|
|
// Punti estremi del segmento
|
|
Point3d ptSt( ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMin, ( nI + 0.5) * m_dStep) ;
|
|
Point3d ptEn( ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMax, ( nI + 0.5) * m_dStep) ;
|
|
// Se il segmento interseca il piano abbiamo finito
|
|
Point3d ptInt ;
|
|
if ( IntersLinePlane( ptSt, ptEn, plPlane, ptInt))
|
|
return true ;
|
|
}
|
|
}
|
|
return false ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int
|
|
VolZmap::CalcIndexForPlaneCells( const Plane3d& plPlane, int nGrid, int nCellI, int nCellJ) const
|
|
{
|
|
int nIndex = 0 ;
|
|
|
|
if ( InOut( plPlane, nGrid, nCellI, nCellJ))
|
|
nIndex |= ( 1 << 0) ;
|
|
if ( InOut( plPlane, nGrid, nCellI + 1, nCellJ))
|
|
nIndex |= ( 1 << 1) ;
|
|
if ( InOut( plPlane, nGrid, nCellI + 1, nCellJ + 1))
|
|
nIndex |= ( 1 << 2) ;
|
|
if ( InOut( plPlane, nGrid, nCellI, nCellJ + 1))
|
|
nIndex |= ( 1 << 3) ;
|
|
|
|
return nIndex ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// dZ è una quota relativa alla mappa nGrid, se il dexel nDex esiste e
|
|
// dZ è interna a un suo intervallo restituisce true, false altrimenti.
|
|
bool
|
|
VolZmap::IsZInsideInterval( int nGrid, int nDex, double dZ) const
|
|
{
|
|
// Se la griglia non esiste vi è un errore.
|
|
if ( nGrid < 0 || nGrid > 2)
|
|
return false ;
|
|
// Se il dexel non esiste vi è un errore.
|
|
if ( nDex < 0 && nDex > int( m_Values[nGrid].size() - 1))
|
|
return false ;
|
|
// Valuto se dZ è interna a un intervallo.
|
|
for ( int nk = 0 ; nk < int( m_Values[nGrid][nDex].size()) ; ++ nk) {
|
|
double dZ1 = m_Values[nGrid][nDex][nk].dMin ;
|
|
double dZ2 = m_Values[nGrid][nDex][nk].dMax ;
|
|
// Se troviamo dZ in un intervallo abbiamo finito.
|
|
if ( dZ > dZ1 && dZ < dZ2)
|
|
return true ;
|
|
}
|
|
// dZ non sta in un intervallo.
|
|
return false ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Box e piano devono essere nello stesso riferimento. Il box deve essere axis aligned.
|
|
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 ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// 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) ;
|
|
}
|