EgtGeomKernel 1.9l4 :

- fabs sostituito da abs
- in Zmap razionalizzazione operazioni taglio spilloni
- in SurfTriMesh UpdateFaceting senza più chiamate recursive.
This commit is contained in:
Dario Sassi
2018-12-27 11:19:40 +00:00
parent a149384fbb
commit 64c954ad4b
56 changed files with 936 additions and 1408 deletions
+325 -401
View File
@@ -693,312 +693,259 @@ VolZmap::AvoidCylinder( const Frame3d& frCyl, double dL, double dR) const
}
//----------------------------------------------------------------------------
// NB: L'origine del sistema di riferimento è nel centro della circonferenza di base
// e l'asse di simmetria coincide con l'asse z.
// 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 dR, bool bTapB, bool bTapT,
Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2)
const Frame3d& CylFrame, double dH, double dRad, bool bTapLow, bool bTapUp,
Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2)
{
Point3d ptP = ptLineSt ;
Vector3d vtV = vtLineDir ;
// Porto la linea nel riferimento del cilindro
Point3d ptP = ptLineSt ; ptP.ToLoc( CylFrame) ;
Vector3d vtV = vtLineDir ; vtV.ToLoc( CylFrame) ;
// Trasformazione delle coordinate:
// l'asse del cilindro corrisponde con
// l'asse z del sistema di riferimento
ptP.ToLoc( CylFrame) ;
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) ;
DBLVECTOR vdRoots ;
// Non vogliamo che i dexel a filo vengano tagliati
double dSqRad = dR * dR ;
double dSqRadSafe = dSqRad - 2 * dR * EPS_SMALL ;
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 ;
// Computo radici
DBLVECTOR vdRoots ;
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 ;
// 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 ;
}
// Nessuna intersezione
else
return false ;
}
// L'equazione ammette o due soluzioni (eventualmente
// coincidenti) oppure nessuna o infinite se la la retta
// appartiene alla superficie
// Due soluzioni: la retta interseca due volte la superficie laterale
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)
// Punti di intersezione con la superficie del cilindro
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) ;
if ( ptInt1.z > ptInt2.z)
swap( ptInt1, ptInt2) ;
// Determino le normali
vtN1.Set( -ptInt1.x, -ptInt1.y, 0) ;
vtN1.Normalize() ;
vtN2.Set( ( ORIG - ptInt2).x, ( ORIG - ptInt2).y, 0) ;
vtN2.Set( -ptInt2.x, -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
// Porto i punti e i versori nel riferimento globale
ptInt1.ToGlob( CylFrame) ;
ptInt2.ToGlob( CylFrame) ;
vtN1.ToGlob( CylFrame) ;
ptInt2.ToGlob( CylFrame) ;
vtN2.ToGlob( CylFrame) ;
// Trovate intersezioni
return true ;
}
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 dh, double dH, bool bTapLow, bool bTapUp,
Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2)
const Frame3d& ConusFrame, double dTan, double dMinH, double dMaxH, bool bTapLow, bool bTapUp,
Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, 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.
// Porto la linea nel riferimento del cono
Point3d ptP = ptLineSt ; ptP.ToLoc( ConusFrame) ;
Vector3d vtV = vtLineDir ; vtV.ToLoc( ConusFrame) ;
Point3d ptP = ptLineSt ;
Vector3d vtV = vtLineDir ;
// Raggi delle due basi
double dMinRad = dTan * dMinH ;
double dMaxRad = dTan * dMaxH ;
// Trasformazione delle coordinate
ptP.ToLoc( ConusFrame) ;
vtV.ToLoc( ConusFrame) ;
// Epsilon per piani di tappo
double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ;
double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ;
DBLVECTOR vdCoef(3) ;
DBLVECTOR vdRoots ;
// 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 ;
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
DBLVECTOR vdRoots ;
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 ;
// 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
else if ( nRoot == 2) {
// 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) ;
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 ;
// 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() ;
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) ;
// 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() ;
vtN2.Normalize() ;
return true ;
}
// 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 ;
}
return false ;
// Nessuna soluzione : nessuna intersezione
else
return false ;
}
//----------------------------------------------------------------------------
// NB: L'origine del sistema di riferimento deve essere
// 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.
@@ -1008,26 +955,25 @@ VolZmap::IntersLineConus( const Point3d& ptLineSt, const Vector3d& vtLineDir,
// 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,
VolZmap::IntersLineEllipticalCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir,
const Frame3d& CircFrame, double dRad, double dLongMvLen, double dOrtMvLen,
bool bTapLow, bool bTapUp,
Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2)
Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2)
{
// 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 ;
// 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
@@ -1049,14 +995,12 @@ VolZmap::IntersLineEllipticalCylinder( const Vector3d& vtLineDir, const Point3d&
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
// 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
@@ -1064,120 +1008,103 @@ VolZmap::IntersLineEllipticalCylinder( const Vector3d& vtLineDir, const Point3d&
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 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)
// Se numero di soluzioni diverso da due le eventuali intersezioni sono già state trovate
if ( nRoot != 2)
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() ;
// 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() ;
// 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 ;
// 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
return false ;
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 ;
// 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, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2)
Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2)
{
// Controllo sulle dimensioni lineari affinché sia valido il poliedro
if ( dLenX <= 0 || dLenY <= 0 || dLenZ <= 0)
if ( dLenX <= EPS_SMALL || dLenY <= EPS_SMALL || dLenZ <= EPS_SMALL)
return false ;
Point3d ptP = ptLineSt ;
Vector3d vtV = vtLineDir ;
// Trasformazione delle coordinate
ptP.ToLoc( PolyFrame) ;
vtV.ToLoc( PolyFrame) ;
// 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
@@ -1185,13 +1112,27 @@ VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLine
Point3d ptFacet135( 0, dLenY / 2, 0) ;
Point3d ptFacet246( dLenX, - dLenY / 2, dLenZ + dDeltaZ) ;
// Servono per descrivere i piani obliqui
// 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 ;
@@ -1199,37 +1140,20 @@ VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLine
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) > 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 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) {
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 ( 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 ;
@@ -1243,9 +1167,9 @@ VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLine
}
// Intersezione con la terza faccia
if ( nIntNum < 2 &&
ptI3.z >= 0 && ptI3.z <= dLenZ &&
ptI3.y >= - ptFacet135.y && ptI3.y <= ptFacet135.y) {
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 ;
@@ -1259,9 +1183,9 @@ VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLine
}
// 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 < 2 &&
ptI4.z >= dDeltaZ && ptI4.z <= dLenZ + dDeltaZ &&
ptI4.y >= - ptFacet135.y && ptI4.y <= ptFacet135.y) {
if ( nIntNum == 0) {
ptInt1 = ptI4 ;
vtN1 = - X_AX ;
@@ -1275,9 +1199,9 @@ VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLine
}
// Intersezione con la quinta faccia
if ( nIntNum < 2 &&
ptI5.x >= 0 && ptI5.x <= dLenX &&
ptI5.y >= - ptFacet135.y && ptI5.y <= ptFacet135.y) {
if ( nIntNum < 2 &&
ptI5.x >= 0 && ptI5.x <= dLenX &&
ptI5.y >= - ptFacet135.y && ptI5.y <= ptFacet135.y) {
if ( nIntNum == 0) {
ptInt1 = ptI5 ;
vtN1 = vtOb ;
@@ -1291,9 +1215,9 @@ VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLine
}
// Intersezione con la sesta faccia
if ( nIntNum < 2 &&
ptI6.x >= 0 && ptI6.x <= dLenX &&
ptI6.y >= - ptFacet135.y && ptI6.y <= ptFacet135.y) {
if ( nIntNum < 2 &&
ptI6.x >= 0 && ptI6.x <= dLenX &&
ptI6.y >= - ptFacet135.y && ptI6.y <= ptFacet135.y) {
if ( nIntNum == 0) {
ptInt1 = ptI6;
vtN1 = - vtOb ;
@@ -1306,15 +1230,15 @@ VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLine
}
}
if ( nIntNum == 2) {
ptInt1.ToGlob( PolyFrame) ;
ptInt2.ToGlob( PolyFrame) ;
vtN1.ToGlob( PolyFrame) ;
vtN2.ToGlob( PolyFrame) ;
return true ;
}
else
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 ;
}
//----------------------------------------------------------------------------