Dario Sassi
3812 lines
168 KiB
C++
3812 lines
168 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 "MC_Tables.h"
|
|
#include "/EgtDev/Include/EGkIntervals.h"
|
|
#include "/EgtDev/Include/EgtNumUtils.h"
|
|
#include "/EgtDev/Include/EGkStringUtils3d.h"
|
|
#include "/EgtDev/Extern/Eigen\Core"
|
|
#include "/EgtDev/Extern/Eigen\SVD"
|
|
#include "PolygonPlane.h"
|
|
|
|
|
|
using namespace std ;
|
|
|
|
|
|
|
|
// ------------------------- UNORDERED MAP PER LA GESTIONE DEI TRIANGOLI GRANDI ---------------------------------------------------
|
|
|
|
// ------------------------- TABELLA BLOCCHI ADIACENTI ----------------------------------------------------------------------------
|
|
static int NeighbourTable[8][4] = {
|
|
{0, -1, -1, -1},
|
|
{1, 1, -1, -1},
|
|
{1, 1, 2, -1},
|
|
{2, 1, 2, -1},
|
|
{1, 3, -1, -1},
|
|
{2, 1, 3, -1},
|
|
{2, 2, 3, -1},
|
|
{3, 1, 2, 3}
|
|
} ;
|
|
|
|
// ------------------------- FUNZIONE TEST SULLE NORMALI --------------------------------------------------------------------------
|
|
enum FatureType { NO_FEATURE = 0, CORNER = 1, EDGE = 2} ;
|
|
enum CanonicDir { X_PLUS = 1, X_MINUS = -1, Y_PLUS = 2, Y_MINUS = -2, Z_PLUS = 3, Z_MINUS = -3} ;
|
|
|
|
//----------------------------------------------------------------------------
|
|
int
|
|
TestOnNormal( const VectorField CompoVert[], int nCompoElem)
|
|
{
|
|
// Cerco la massima deviazione tra le normali nei punti della parte connessa
|
|
int nI, nJ ;
|
|
double dMinCosTheta = 2 ;
|
|
for ( int i = 0 ; i < nCompoElem ; ++ i) {
|
|
for ( int j = i + 1 ; j < nCompoElem ; ++ j) {
|
|
double dCurrCos = CompoVert[i].vtNorm * CompoVert[j].vtNorm ;
|
|
if ( dCurrCos < dMinCosTheta) {
|
|
dMinCosTheta = dCurrCos ;
|
|
nI = i ;
|
|
nJ = j ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Se la massima deviazione non supera il limite non è feature
|
|
const double SHARP_COS = 0.9 ;
|
|
if ( dMinCosTheta >= SHARP_COS)
|
|
return NO_FEATURE ;
|
|
|
|
// Verifico se Edge o Corner
|
|
// direzione perpendicolare alle normali con massima differenza (potenziale edge)
|
|
Vector3d vtK = CompoVert[nI].vtNorm ^ CompoVert[nJ].vtNorm ;
|
|
vtK.Normalize() ;
|
|
// cerco normale con massima vicinanza al potenziale edge
|
|
double dMaxAbsCos = 0 ;
|
|
for ( int i = 0 ; i < nCompoElem ; ++ i) {
|
|
double dAbsCurrentCos = abs( CompoVert[i].vtNorm * vtK) ;
|
|
if ( dAbsCurrentCos > dMaxAbsCos)
|
|
dMaxAbsCos = dAbsCurrentCos ;
|
|
}
|
|
// se esiste normale diretta quasi come potenziale edge, allora corner
|
|
const double CORNER_COS = 0.7 ;
|
|
if ( dMaxAbsCos > CORNER_COS)
|
|
return CORNER ;
|
|
else
|
|
return EDGE ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
CanonicPlaneTest( const VectorField CompoVert[], int nDir, double& dPos)
|
|
{
|
|
// Verifico posizione dei punti
|
|
int nI ;
|
|
switch ( nDir) {
|
|
case X_PLUS : case X_MINUS : nI = 0 ; break ;
|
|
case Y_PLUS : case Y_MINUS : nI = 1 ; break ;
|
|
case Z_PLUS : case Z_MINUS : nI = 2 ; break ;
|
|
}
|
|
double dPos0 = CompoVert[0].ptInt.v[nI] ;
|
|
double dPos1 = CompoVert[1].ptInt.v[nI] ;
|
|
double dPos2 = CompoVert[2].ptInt.v[nI] ;
|
|
double dPos3 = CompoVert[3].ptInt.v[nI] ;
|
|
dPos = ( dPos0 + dPos1 + dPos2 + dPos3) / 4 ;
|
|
if ( abs( dPos0 - dPos) > EPS_SMALL || abs( dPos1 - dPos) > EPS_SMALL ||
|
|
abs( dPos2 - dPos) > EPS_SMALL || abs( dPos3 - dPos) > EPS_SMALL)
|
|
return false ;
|
|
// Verifico direzione normali
|
|
Vector3d vtN ;
|
|
switch ( nDir) {
|
|
case X_PLUS : vtN = X_AX ; break ;
|
|
case X_MINUS : vtN = -X_AX ; break ;
|
|
case Y_PLUS : vtN = Y_AX ; break ;
|
|
case Y_MINUS : vtN = - Y_AX ; break ;
|
|
case Z_PLUS : vtN = Z_AX ; break ;
|
|
case Z_MINUS : vtN = - Z_AX ; break ;
|
|
}
|
|
for ( int i = 0 ; i < 4 ; ++ i) {
|
|
if ( CompoVert[i].vtNorm * vtN < 0.999)
|
|
return false ;
|
|
}
|
|
// Superati tutti i test
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
DotTest( const VectorField CompoVert[], int nCompoElem, Vector3d& vtAvg, double dThreshold = 0)
|
|
{
|
|
// Cerco la massima deviazione tra le normali nei punti della parte connessa
|
|
double dMinCosTheta = 2 ;
|
|
for ( int i = 0 ; i < nCompoElem ; ++ i) {
|
|
for ( int j = i + 1 ; j < nCompoElem ; ++ j) {
|
|
double dCurrCos = CompoVert[i].vtNorm * CompoVert[j].vtNorm ;
|
|
if ( dCurrCos < dMinCosTheta) {
|
|
dMinCosTheta = dCurrCos ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// se normali sparpagliate oltre limite
|
|
if ( dMinCosTheta < dThreshold)
|
|
return false ;
|
|
|
|
// determino media delle normali
|
|
vtAvg = V_NULL ;
|
|
for ( int i = 0 ; i < nCompoElem ; ++ i)
|
|
vtAvg += CompoVert[i].vtNorm ;
|
|
vtAvg /= nCompoElem ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GridControl( VectorField VecField[], VectorField CompoVert[][12],
|
|
int nIndArrey[][4], int nVertComp[], int nCompCount, bool& bGridControl) const
|
|
{
|
|
// Ordino i 4 indici in senso crescente
|
|
for ( int nSrtInd1 = 0 ; nSrtInd1 < nVertComp[nCompCount - 1] - 1 ; ++ nSrtInd1) {
|
|
for ( int nSrtInd2 = nSrtInd1 + 1 ; nSrtInd2 < nVertComp[nCompCount - 1] ; ++ nSrtInd2) {
|
|
if ( nIndArrey[nCompCount - 1][nSrtInd1] > nIndArrey[nCompCount - 1][nSrtInd2])
|
|
swap( nIndArrey[nCompCount - 1][nSrtInd1], nIndArrey[nCompCount - 1][nSrtInd2]) ;
|
|
}
|
|
}
|
|
|
|
if ( ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 2 &&
|
|
nIndArrey[nCompCount - 1][2] == 9 && nIndArrey[nCompCount - 1][3] == 10) ||
|
|
( nIndArrey[nCompCount - 1][0] == 4 && nIndArrey[nCompCount - 1][1] == 6 &&
|
|
nIndArrey[nCompCount - 1][2] == 9 && nIndArrey[nCompCount - 1][3] == 10) ||
|
|
( nIndArrey[nCompCount - 1][0] == 4 && nIndArrey[nCompCount - 1][1] == 6 &&
|
|
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 11) ||
|
|
( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 2 &&
|
|
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 11) ||
|
|
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 3 &&
|
|
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 9 ) ||
|
|
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 3 &&
|
|
nIndArrey[nCompCount - 1][2] == 10 && nIndArrey[nCompCount - 1][3] == 11) ||
|
|
( nIndArrey[nCompCount - 1][0] == 5 && nIndArrey[nCompCount - 1][1] == 7 &&
|
|
nIndArrey[nCompCount - 1][2] == 10 && nIndArrey[nCompCount - 1][3] == 11) ||
|
|
( nIndArrey[nCompCount - 1][0] == 5 && nIndArrey[nCompCount - 1][1] == 7 &&
|
|
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 9 )) {
|
|
|
|
|
|
if ( AreSameVectorApprox( VecField[nIndArrey[nCompCount - 1][0]].vtNorm, VecField[nIndArrey[nCompCount - 1][1]].vtNorm) &&
|
|
abs( VecField[nIndArrey[nCompCount - 1][0]].vtNorm * VecField[nIndArrey[nCompCount - 1][2]].vtNorm) < EPS_SMALL &&
|
|
abs( VecField[nIndArrey[nCompCount - 1][0]].vtNorm * VecField[nIndArrey[nCompCount - 1][3]].vtNorm) < EPS_SMALL) {
|
|
|
|
Point3d ptBarycenter = ( CompoVert[nCompCount - 1][0].ptInt + CompoVert[nCompCount - 1][1].ptInt +
|
|
CompoVert[nCompCount - 1][2].ptInt + CompoVert[nCompCount - 1][3].ptInt) / 4 ;
|
|
|
|
Vector3d vtVecField = VecField[nIndArrey[nCompCount - 1][0]].vtNorm ;
|
|
|
|
if ( vtVecField.x * vtVecField.x > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.x - ptBarycenter.x) < EPS_SMALL) {
|
|
|
|
double dXBar = ptBarycenter.x / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNx[0]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dXInd = double( nBarInd) ;
|
|
if ( abs( dXInd - dXBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
else if ( vtVecField.y * vtVecField.y > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.y - ptBarycenter.y) < EPS_SMALL) {
|
|
|
|
double dYBar = ptBarycenter.y / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNy[0]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dYInd = double( nBarInd) ;
|
|
if ( abs( dYInd - dYBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
else if ( vtVecField.z * vtVecField.z > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.z - ptBarycenter.z) < EPS_SMALL) {
|
|
|
|
double dZBar = ptBarycenter.z / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNy[1]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dZInd = double( nBarInd) ;
|
|
if ( abs( dZInd - dZBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else if ( AreSameVectorApprox( VecField[nIndArrey[nCompCount - 1][2]].vtNorm, VecField[nIndArrey[nCompCount - 1][3]].vtNorm) &&
|
|
abs( VecField[nIndArrey[nCompCount - 1][2]].vtNorm * VecField[nIndArrey[nCompCount - 1][0]].vtNorm) < EPS_SMALL &&
|
|
abs( VecField[nIndArrey[nCompCount - 1][2]].vtNorm * VecField[nIndArrey[nCompCount - 1][1]].vtNorm) < EPS_SMALL) {
|
|
|
|
Point3d ptBarycenter = ( CompoVert[nCompCount - 1][0].ptInt + CompoVert[nCompCount - 1][1].ptInt +
|
|
CompoVert[nCompCount - 1][2].ptInt + CompoVert[nCompCount - 1][3].ptInt) / 4 ;
|
|
|
|
Vector3d vtVecField = VecField[nIndArrey[nCompCount - 1][2]].vtNorm ;
|
|
|
|
if ( vtVecField.x * vtVecField.x > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.x - ptBarycenter.x) < EPS_SMALL) {
|
|
|
|
double dXBar = ptBarycenter.x / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNx[0]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dXInd = double( nBarInd) ;
|
|
if ( abs( dXInd - dXBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
else if ( vtVecField.y * vtVecField.y > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.y - ptBarycenter.y) < EPS_SMALL) {
|
|
|
|
double dYBar = ptBarycenter.y / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNy[0]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dYInd = double( nBarInd) ;
|
|
if ( abs( dYInd - dYBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
else if ( vtVecField.z * vtVecField.z > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.z - ptBarycenter.z) < EPS_SMALL) {
|
|
|
|
double dZBar = ptBarycenter.z / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNy[1]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dZInd = double( nBarInd) ;
|
|
if ( abs( dZInd - dZBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
else if ( ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 1 &&
|
|
nIndArrey[nCompCount - 1][2] == 4 && nIndArrey[nCompCount - 1][3] == 5) ||
|
|
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 2 &&
|
|
nIndArrey[nCompCount - 1][2] == 5 && nIndArrey[nCompCount - 1][3] == 6) ||
|
|
( nIndArrey[nCompCount - 1][0] == 2 && nIndArrey[nCompCount - 1][1] == 3 &&
|
|
nIndArrey[nCompCount - 1][2] == 6 && nIndArrey[nCompCount - 1][3] == 7) ||
|
|
( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 3 &&
|
|
nIndArrey[nCompCount - 1][2] == 4 && nIndArrey[nCompCount - 1][3] == 7)) {
|
|
|
|
|
|
if ( AreSameVectorApprox( VecField[nIndArrey[nCompCount - 1][0]].vtNorm, VecField[nIndArrey[nCompCount - 1][2]].vtNorm) &&
|
|
abs( VecField[nIndArrey[nCompCount - 1][0]].vtNorm * VecField[nIndArrey[nCompCount - 1][1]].vtNorm) < EPS_SMALL &&
|
|
abs( VecField[nIndArrey[nCompCount - 1][0]].vtNorm * VecField[nIndArrey[nCompCount - 1][3]].vtNorm) < EPS_SMALL) {
|
|
|
|
Point3d ptBarycenter = ( CompoVert[nCompCount - 1][0].ptInt + CompoVert[nCompCount - 1][1].ptInt +
|
|
CompoVert[nCompCount - 1][2].ptInt + CompoVert[nCompCount - 1][3].ptInt) / 4 ;
|
|
|
|
Vector3d vtVecField = VecField[nIndArrey[nCompCount - 1][0]].vtNorm ;
|
|
|
|
if ( vtVecField.x * vtVecField.x > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.x - ptBarycenter.x) < EPS_SMALL) {
|
|
|
|
double dXBar = ptBarycenter.x / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNx[0]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dXInd = double( nBarInd) ;
|
|
if ( abs( dXInd - dXBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
else if ( vtVecField.y * vtVecField.y > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.y - ptBarycenter.y) < EPS_SMALL) {
|
|
|
|
double dYBar = ptBarycenter.y / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNy[0]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dYInd = double( nBarInd) ;
|
|
if ( abs( dYInd - dYBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
else if ( vtVecField.z * vtVecField.z > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.z - ptBarycenter.z) < EPS_SMALL) {
|
|
|
|
double dZBar = ptBarycenter.z / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNy[1]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dZInd = double( nBarInd) ;
|
|
if ( abs( dZInd - dZBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else if ( AreSameVectorApprox( VecField[nIndArrey[nCompCount - 1][1]].vtNorm, VecField[nIndArrey[nCompCount - 1][3]].vtNorm) &&
|
|
abs( VecField[nIndArrey[nCompCount - 1][1]].vtNorm * VecField[nIndArrey[nCompCount - 1][0]].vtNorm) < EPS_SMALL &&
|
|
abs( VecField[nIndArrey[nCompCount - 1][1]].vtNorm * VecField[nIndArrey[nCompCount - 1][2]].vtNorm) < EPS_SMALL) {
|
|
|
|
Point3d ptBarycenter = ( CompoVert[nCompCount - 1][0].ptInt + CompoVert[nCompCount - 1][1].ptInt +
|
|
CompoVert[nCompCount - 1][2].ptInt + CompoVert[nCompCount - 1][3].ptInt) / 4 ;
|
|
|
|
Vector3d vtVecField = VecField[nIndArrey[nCompCount - 1][0]].vtNorm ;
|
|
|
|
if ( vtVecField.x * vtVecField.x > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.x - ptBarycenter.x) < EPS_SMALL) {
|
|
|
|
double dXBar = ptBarycenter.x / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNx[0]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dXInd = double( nBarInd) ;
|
|
if ( abs( dXInd - dXBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
else if ( vtVecField.y * vtVecField.y > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.y - ptBarycenter.y) < EPS_SMALL) {
|
|
|
|
double dYBar = ptBarycenter.y / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNy[0]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dYInd = double( nBarInd) ;
|
|
if ( abs( dYInd - dYBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
else if ( vtVecField.z * vtVecField.z > 1 - EPS_SMALL * EPS_SMALL) {
|
|
|
|
if ( abs( CompoVert[nCompCount - 1][0].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][1].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][2].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
|
|
abs( CompoVert[nCompCount - 1][3].ptInt.z - ptBarycenter.z) < EPS_SMALL) {
|
|
|
|
double dZBar = ptBarycenter.z / m_dStep - 0.5 ;
|
|
int nBarLimSup = int( m_nNy[1]) ;
|
|
int nBarInd = 0 ;
|
|
|
|
while ( nBarInd < nBarLimSup) {
|
|
double dZInd = double( nBarInd) ;
|
|
if ( abs( dZInd - dZBar) < EPS_SMALL) {
|
|
bGridControl = false ;
|
|
break ;
|
|
}
|
|
++ nBarInd ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return true ;
|
|
}
|
|
|
|
// ------------------------- VISUALIZZAZIONE --------------------------------------------------------------------------------------
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetDexelLines( int nDir, int nPos1, int nPos2, POLYLINELIST& lstPL) const
|
|
{
|
|
// Se richiesti spilloni ( 0 <= nDir < 3)
|
|
if ( nDir < 3) {
|
|
|
|
// Controllo l'ammissibilità della griglia
|
|
if ( nDir < 0 || nDir > 2)
|
|
return false ;
|
|
|
|
// Verifiche sugli indici
|
|
if ( nPos1 < 0 || nPos1 >= int( m_nNx[nDir]) || nPos2 < 0 || nPos2 >= int( m_nNy[nDir]))
|
|
return false ;
|
|
|
|
int nPos = nPos1 + nPos2 * m_nNx[nDir] ;
|
|
|
|
if ( nPos < 0 || nPos >= int( m_Values[nDir].size()))
|
|
return false ;
|
|
|
|
// Definisco un sistema di riferimento opportuno
|
|
Frame3d frMapFrame ;
|
|
if ( nDir == 0)
|
|
frMapFrame = m_MapFrame ;
|
|
else if ( nDir == 1)
|
|
frMapFrame.Set( m_MapFrame.Orig(), m_MapFrame.VersY(), m_MapFrame.VersZ(), m_MapFrame.VersX()) ;
|
|
else if ( nDir == 2)
|
|
frMapFrame.Set( m_MapFrame.Orig(), m_MapFrame.VersZ(), m_MapFrame.VersX(), m_MapFrame.VersY()) ;
|
|
|
|
// Calcolo coordinate punto
|
|
double dX = m_dStep * ( 0.5 + nPos1) ;
|
|
double dY = m_dStep * ( 0.5 + nPos2) ;
|
|
|
|
// Determino il punto di partenza sulla griglia
|
|
Point3d ptP = frMapFrame.Orig() + dX * frMapFrame.VersX() + dY * frMapFrame.VersY() ;
|
|
|
|
// Creo le polilinee
|
|
for ( int j = 0 ; j < int( m_Values[nDir][nPos].size()) ; j += 1) {
|
|
// aggiungo polilinea a lista
|
|
lstPL.emplace_back() ;
|
|
// inserisco punti estremi
|
|
lstPL.back().AddUPoint( 0, ptP + m_Values[nDir][nPos][j].dMin * frMapFrame.VersZ()) ;
|
|
lstPL.back().AddUPoint( 1, ptP + m_Values[nDir][nPos][j].dMax * frMapFrame.VersZ()) ;
|
|
}
|
|
return true ;
|
|
}
|
|
// altrimenti richieste normali ( 3 <= nDir < 6)
|
|
else {
|
|
|
|
// riporto a indice griglia
|
|
nDir -= 3 ;
|
|
|
|
// Controllo l'ammissibilità della griglia
|
|
if ( nDir < 0 || nDir > 2)
|
|
return false ;
|
|
|
|
// Verifiche sugli indici
|
|
if ( nPos1 < 0 || nPos1 >= int( m_nNx[nDir]) || nPos2 < 0 || nPos2 >= int( m_nNy[nDir]))
|
|
return false ;
|
|
|
|
int nPos = nPos1 + nPos2 * m_nNx[nDir] ;
|
|
|
|
if ( nPos < 0 || nPos >= int( m_Values[nDir].size()))
|
|
return false ;
|
|
|
|
// Definisco un sistema di riferimento opportuno
|
|
Frame3d frMapFrame ;
|
|
if ( nDir == 0)
|
|
frMapFrame = m_MapFrame ;
|
|
else if ( nDir == 1)
|
|
frMapFrame.Set( m_MapFrame.Orig(), m_MapFrame.VersY(), m_MapFrame.VersZ(), m_MapFrame.VersX()) ;
|
|
else if ( nDir == 2)
|
|
frMapFrame.Set( m_MapFrame.Orig(), m_MapFrame.VersZ(), m_MapFrame.VersX(), m_MapFrame.VersY()) ;
|
|
|
|
// Calcolo coordinate punto
|
|
double dX = m_dStep * ( 0.5 + nPos1) ;
|
|
double dY = m_dStep * ( 0.5 + nPos2) ;
|
|
|
|
// Determino il punto di partenza sulla griglia
|
|
Point3d ptP = frMapFrame.Orig() + dX * frMapFrame.VersX() + dY * frMapFrame.VersY() ;
|
|
|
|
// Creo le polilinee
|
|
for ( int j = 0 ; j < int( m_Values[nDir][nPos].size()) ; j += 1) {
|
|
// aggiungo polilinea a lista
|
|
lstPL.emplace_back() ;
|
|
// calcolo e inserisco punto inizio spillone
|
|
Point3d ptQ = ptP + m_Values[nDir][nPos][j].dMin * frMapFrame.VersZ() ;
|
|
lstPL.back().AddUPoint( 0, ptQ) ;
|
|
// calcolo e inserisco punto su termine sua normale
|
|
Vector3d vtV = m_Values[nDir][nPos][j].vtMinN ;
|
|
vtV.ToGlob( m_MapFrame) ;
|
|
lstPL.back().AddUPoint( 1, ptQ + vtV * m_dStep / 4) ;
|
|
// aggiungo polilinea a lista
|
|
lstPL.emplace_back() ;
|
|
// calcolo e inserisco punto fine spillone
|
|
Point3d ptR = ptP + m_Values[nDir][nPos][j].dMax * frMapFrame.VersZ() ;
|
|
lstPL.back().AddUPoint( 0, ptR) ;
|
|
// calcolo e inserisco punto su termine sua normale
|
|
Vector3d vtW = m_Values[nDir][nPos][j].vtMaxN ;
|
|
vtW.ToGlob( m_MapFrame) ;
|
|
lstPL.back().AddUPoint( 1, ptR + vtW * m_dStep / 4) ;
|
|
}
|
|
return true ;
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetAllTriangles( TRIA3DLIST& lstTria) const
|
|
{
|
|
INTVECTOR nModifiedBlocks ;
|
|
TRIA3DLISTVECTOR vLstTria ;
|
|
if ( ! GetTriangles( true, nModifiedBlocks, vLstTria))
|
|
return false ;
|
|
lstTria.clear() ;
|
|
for ( size_t i = 0 ; i < vLstTria.size() ; ++ i) {
|
|
lstTria.splice( lstTria.end(), vLstTria[i]) ;
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVECTOR& vLstTria) const
|
|
{
|
|
// Se nessun blocco modificato, è richiesta esterna e li considero tutti modificati
|
|
bool bSomeModif = false ;
|
|
for ( size_t i = 0 ; i < m_nNumBlock ; ++ i) {
|
|
if ( m_BlockToUpdate[i]) {
|
|
bSomeModif = true ;
|
|
break ;
|
|
}
|
|
}
|
|
if ( ! bSomeModif)
|
|
bAllBlocks = true ;
|
|
|
|
// Caso di singola mappa
|
|
if ( m_nMapNum == 1) {
|
|
|
|
const int MAX_DIM_CHUNK = 128 ;
|
|
|
|
nModifiedBlocks.resize( m_nNumBlock) ;
|
|
vLstTria.reserve( m_nNumBlock) ;
|
|
|
|
// Ciclo sui blocchi
|
|
for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) {
|
|
|
|
// Se il blocco deve essere aggiornato, eseguo
|
|
if ( bAllBlocks || m_BlockToUpdate[t]) {
|
|
|
|
// preparo lista
|
|
vLstTria.emplace_back() ;
|
|
nModifiedBlocks[t] = int( vLstTria.size()) - 1 ;
|
|
|
|
// Calcolo posizione del blocco nella griglia
|
|
int nIBlock = int( t) % int( m_nFracLin[0]) ;
|
|
int nJBlock = int( t) / int( m_nFracLin[0]) ;
|
|
|
|
// Calcolo limiti per l'indice i
|
|
int nStartI = nIBlock * int( m_nVoxNumPerBlock) * N_DEXVOXRATIO ;
|
|
int nEndI = ( nIBlock + 1 == int( m_nFracLin[0]) ?
|
|
int( m_nNx[0]) : ( nIBlock + 1) * int( m_nVoxNumPerBlock)) * N_DEXVOXRATIO ;
|
|
|
|
// Calcolo limiti per l'indice j
|
|
int nStartJ = nJBlock * int( m_nVoxNumPerBlock) * N_DEXVOXRATIO ;
|
|
int nEndJ = ( nJBlock + 1 == int( m_nFracLin[1]) ?
|
|
int( m_nNy[0]) : ( nJBlock + 1) * int( m_nVoxNumPerBlock)) * N_DEXVOXRATIO ;
|
|
|
|
// Ciclo su i e j
|
|
for ( int i = nStartI ; i < nEndI ; i += MAX_DIM_CHUNK) {
|
|
int nDimChunkX = min( MAX_DIM_CHUNK, nEndI - i) ;
|
|
for ( int j = nStartJ ; j < nEndJ ; j += MAX_DIM_CHUNK) {
|
|
int nDimChunkY = min( MAX_DIM_CHUNK, nEndJ - j) ;
|
|
GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, MAX_DIM_CHUNK, vLstTria.back()) ;
|
|
}
|
|
}
|
|
|
|
m_BlockToUpdate[t] = false ;
|
|
}
|
|
else
|
|
nModifiedBlocks[t] = -1 ;
|
|
}
|
|
}
|
|
|
|
// Caso con tre mappe
|
|
else {
|
|
|
|
nModifiedBlocks.resize( m_nNumBlock + 1) ;
|
|
vLstTria.reserve( m_nNumBlock + 1) ;
|
|
|
|
TriaMatrix VecTriHold ;
|
|
VecTriHold.resize( m_nNumBlock) ;
|
|
|
|
bool bCalcInterBlock = false ;
|
|
|
|
// Calcolo i triangoli sui blocchi
|
|
for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) {
|
|
|
|
// Se il blocco deve essere processato
|
|
if ( bAllBlocks || m_BlockToUpdate[t]) {
|
|
|
|
// processo ...
|
|
vLstTria.emplace_back() ;
|
|
nModifiedBlocks[t] = int( vLstTria.size()) - 1 ;
|
|
m_InterBlockTria[t].clear() ;
|
|
#if 1
|
|
ExtMarchingCubes( int( t), vLstTria.back(), VecTriHold[t]) ;
|
|
// Flipping fra voxel interni
|
|
FlipEdgesII( VecTriHold[t], true) ;
|
|
bCalcInterBlock = true ;
|
|
#else
|
|
MarchingCubes( int( t), vLstTria.back()) ;
|
|
#endif
|
|
m_BlockToUpdate[t] = false ;
|
|
}
|
|
else
|
|
nModifiedBlocks[t] = -1 ;
|
|
}
|
|
|
|
// Calcolo i triangoli di frontiera tra feature di blocchi diversi
|
|
// copio i triangoli di frontiera in una matrice gemella
|
|
// di m_InterBlockTria per avere sempre a disposizione
|
|
// i triangoli non flippati.
|
|
TriaMatrix InterBlockTria ;
|
|
if ( bCalcInterBlock) {
|
|
InterBlockTria = m_InterBlockTria ;
|
|
FlipEdgesBB( InterBlockTria) ;
|
|
}
|
|
|
|
// Inserisco in lista i triangoli di feature derivanti dai blocchi
|
|
for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) {
|
|
if ( nModifiedBlocks[t] >= 0) {
|
|
// ciclo sui voxel del blocco
|
|
for ( size_t t1 = 0 ; t1 < VecTriHold[t].size() ; ++ t1) {
|
|
// ciclo sulle componenti connesse del voxel
|
|
for ( size_t t2 = 0 ; t2 < VecTriHold[t][t1].vCompoTria.size() ; ++ t2) {
|
|
// ciclo sui triangoli delle componenti connesse
|
|
for ( size_t t3 = 0 ; t3 < VecTriHold[t][t1].vCompoTria[t2].size() ; ++ t3) {
|
|
// aggiungo triangolo alla lista
|
|
vLstTria[nModifiedBlocks[t]].emplace_back( VecTriHold[t][t1].vCompoTria[t2][t3]) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Inserisco in lista i triangoli di frontiera tra feature di blocchi diversi
|
|
if ( bCalcInterBlock) {
|
|
vLstTria.resize( vLstTria.size() + 1) ;
|
|
size_t nPos = size_t( vLstTria.size() - 1) ;
|
|
|
|
for ( size_t t = 0 ; t < m_InterBlockTria.size() ; ++ t) {
|
|
for ( size_t t1 = 0 ; t1 < m_InterBlockTria[t].size() ; ++ t1) {
|
|
for ( size_t t2 = 0 ; t2 < m_InterBlockTria[t][t1].vCompoTria.size() ; ++ t2) {
|
|
for ( size_t t3 = 0 ; t3 < m_InterBlockTria[t][t1].vCompoTria[t2].size() ; ++ t3) {
|
|
if ( m_InterBlockTria[t][t1].vCompoTria[t2][t3].GetArea() > SQ_EPS_SMALL) {
|
|
vLstTria[nPos].emplace_back( InterBlockTria[t][t1].vCompoTria[t2][t3]) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
nModifiedBlocks.back() = int( nPos) ;
|
|
}
|
|
|
|
else
|
|
nModifiedBlocks.back() = - 1 ;
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int
|
|
VolZmap::GetBlockCount( void) const
|
|
{
|
|
return m_nNumBlock + ( m_nMapNum == 1 ? 0 : 1) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, int nDimChk, TRIA3DLIST& lstTria) const
|
|
{
|
|
// determino se è un semplice parallelepipedo
|
|
bool bIsSimple = true ;
|
|
double dBotZ ;
|
|
double dTopZ ;
|
|
for ( int i = 0 ; i < nDim1 && bIsSimple ; ++ i) {
|
|
for ( int j = 0 ; j < nDim2 && bIsSimple ; ++ j) {
|
|
int nPos = ( nPos1 + i) + ( nPos2 + j) * m_nNx[0] ;
|
|
if ( nPos > int( m_nDim[0]) ||
|
|
int( m_Values[0][nPos].size()) != 1)
|
|
bIsSimple = false ;
|
|
else if ( i == 0 && j == 0) {
|
|
dBotZ = m_Values[0][nPos][0].dMin ;
|
|
dTopZ = m_Values[0][nPos][0].dMax ;
|
|
}
|
|
else if ( abs( m_Values[0][nPos][0].dMin - dBotZ) > EPS_SMALL ||
|
|
abs( m_Values[0][nPos][0].dMax - dTopZ) > EPS_SMALL)
|
|
bIsSimple = false ;
|
|
}
|
|
}
|
|
|
|
// se semplice parallelepipedo
|
|
if ( bIsSimple) {
|
|
CalcChunkPrisms( nPos1, nPos2, nDim1, nDim2, lstTria) ;
|
|
}
|
|
// se chunk di dimensioni accettabili
|
|
else if ( nDimChk >= 4) {
|
|
int nNewDimChk = nDimChk / 2 ;
|
|
for ( int i = nPos1 ; i < int( nPos1 + nDim1) ; i += nNewDimChk) {
|
|
int nDimChunkX = min( nNewDimChk, int( nPos1 + nDim1) - i) ;
|
|
for ( int j = nPos2 ; j < int( nPos2 + nDim2) ; j += nNewDimChk) {
|
|
int nDimChunkY = min( nNewDimChk, int( nPos2 + nDim2) - j) ;
|
|
GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, nNewDimChk, lstTria) ;
|
|
}
|
|
}
|
|
}
|
|
// altrimenti
|
|
else {
|
|
// elaboro ogni singolo dexel
|
|
for ( int i = 0 ; i < nDim1 ; ++ i) {
|
|
for ( int j = 0 ; j < nDim2 ; ++ j) {
|
|
CalcDexelPrisms( nPos1 + i, nPos2 + j, lstTria) ;
|
|
}
|
|
}
|
|
}
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::CalcChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, TRIA3DLIST& lstTria) const
|
|
{
|
|
// verifiche sugli indici
|
|
if ( nPos1 < 0 || nPos1 + nDim1 > int( m_nNx[0]) || nPos2 < 0 || nPos2 + nDim2 > int( m_nNy[0]))
|
|
return false ;
|
|
int nPos = nPos1 + nPos2 * m_nNx[0] ;
|
|
if ( nPos < 0 || nPos >= int( m_nDim[0]))
|
|
return false ;
|
|
|
|
// calcolo coordinate punti
|
|
double dX = m_dStep * nPos1 ;
|
|
double dY = m_dStep * nPos2 ;
|
|
Point3d ptP1 = m_MapFrame.Orig() + dX * m_MapFrame.VersX() + dY * m_MapFrame.VersY() ;
|
|
Point3d ptP2 = ptP1 + nDim1 * m_dStep * m_MapFrame.VersX() ;
|
|
Point3d ptP3 = ptP2 + nDim2 * m_dStep * m_MapFrame.VersY() ;
|
|
Point3d ptP4 = ptP1 + nDim2 * m_dStep * m_MapFrame.VersY() ;
|
|
|
|
// creo le facce sopra e sotto
|
|
Vector3d vtDZt = m_Values[0][nPos][0].dMax * m_MapFrame.VersZ() ;
|
|
Vector3d vtDZb = m_Values[0][nPos][0].dMin * m_MapFrame.VersZ() ;
|
|
// faccia superiore P1t->P2t->P3t->P4t : sempre visibile
|
|
lstTria.emplace_back() ;
|
|
lstTria.back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_MapFrame.VersZ()) ;
|
|
lstTria.emplace_back() ;
|
|
lstTria.back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_MapFrame.VersZ()) ;
|
|
// faccia inferiore P1b->P4b->P3b->P2b : sempre visibile
|
|
lstTria.emplace_back() ;
|
|
lstTria.back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_MapFrame.VersZ()) ;
|
|
lstTria.emplace_back() ;
|
|
lstTria.back().Set( ptP3 + vtDZb, ptP2 + vtDZb, ptP1 + vtDZb, - m_MapFrame.VersZ()) ;
|
|
|
|
// creo le facce laterali
|
|
for ( int j = 0 ; j < nDim2 ; ++ j) {
|
|
int nPosD = nPos + nDim1 - 1 + j * m_nNx[0] ;
|
|
int nPosEst = ( nPos1 + nDim1 - 1 < int( m_nNx[0] - 1) ? nPosD + 1 : - 1) ;
|
|
Point3d ptP2D = ptP2 + j * m_dStep * m_MapFrame.VersY() ;
|
|
Point3d ptP3D = ptP2D + m_dStep * m_MapFrame.VersY() ;
|
|
AddDexelSideFace( nPosD, nPosEst, ptP2D, ptP3D, m_MapFrame.VersZ(), m_MapFrame.VersX(), lstTria) ;
|
|
}
|
|
for ( int i = 0 ; i < nDim1 ; ++ i) {
|
|
int nPosD = nPos + ( nDim2 - 1) * m_nNx[0] + i ;
|
|
int nPosNord = ( nPos2 + nDim2 - 1 < int( m_nNy[0] - 1) ? nPosD + m_nNx[0] : - 1) ;
|
|
Point3d ptP4D = ptP4 + i * m_dStep * m_MapFrame.VersX() ;
|
|
Point3d ptP3D = ptP4D + m_dStep * m_MapFrame.VersX() ;
|
|
AddDexelSideFace( nPosD, nPosNord, ptP3D, ptP4D, m_MapFrame.VersZ(), m_MapFrame.VersY(), lstTria) ;
|
|
}
|
|
for ( int j = 0 ; j < nDim2 ; ++ j) {
|
|
int nPosD = nPos + j * m_nNx[0] ;
|
|
int nPosWest = ( nPos1 > 0 ? nPosD - 1 : - 1) ;
|
|
Point3d ptP1D = ptP1 + j * m_dStep * m_MapFrame.VersY() ;
|
|
Point3d ptP4D = ptP1D + m_dStep * m_MapFrame.VersY() ;
|
|
AddDexelSideFace( nPosD, nPosWest, ptP4D, ptP1D, m_MapFrame.VersZ(), - m_MapFrame.VersX(), lstTria) ;
|
|
}
|
|
for ( int i = 0 ; i < nDim1 ; ++ i) {
|
|
int nPosD = nPos + i ;
|
|
int nPosSud = ( nPos2 > 0 ? nPosD - m_nNx[0] : - 1) ;
|
|
Point3d ptP1D = ptP1 + i * m_dStep * m_MapFrame.VersX() ;
|
|
Point3d ptP2D = ptP1D + m_dStep * m_MapFrame.VersX() ;
|
|
AddDexelSideFace( nPosD, nPosSud, ptP1D, ptP2D, m_MapFrame.VersZ(), - m_MapFrame.VersY(), lstTria) ;
|
|
}
|
|
//
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::CalcDexelPrisms( int nPos1, int nPos2, TRIA3DLIST& lstTria) const
|
|
{
|
|
// verifiche sugli indici
|
|
if ( nPos1 < 0 || nPos1 >= int( m_nNx[0]) || nPos2 < 0 || nPos2 >= int( m_nNy[0]))
|
|
return false ;
|
|
int nPos = nPos1 + nPos2 * m_nNx[0] ;
|
|
if ( nPos < 0 || nPos >= int( m_nDim[0]))
|
|
return false ;
|
|
|
|
// calcolo coordinate punto
|
|
double dX = m_dStep * nPos1 ;
|
|
double dY = m_dStep * nPos2 ;
|
|
Point3d ptP1 = m_MapFrame.Orig() + dX * m_MapFrame.VersX() + dY * m_MapFrame.VersY() ;
|
|
Point3d ptP2 = ptP1 + m_dStep * m_MapFrame.VersX() ;
|
|
Point3d ptP3 = ptP2 + m_dStep * m_MapFrame.VersY() ;
|
|
Point3d ptP4 = ptP1 + m_dStep * m_MapFrame.VersY() ;
|
|
|
|
// creo le facce sopra e sotto di ogni intervallo (sempre visibili)
|
|
for ( int i = 0 ; i < int( m_Values[0][nPos].size()) ; i += 1) {
|
|
Vector3d vtDZt = m_Values[0][nPos][i].dMax * m_MapFrame.VersZ() ;
|
|
Vector3d vtDZb = m_Values[0][nPos][i].dMin * m_MapFrame.VersZ() ;
|
|
// faccia superiore P1t->P2t->P3t->P4t : sempre visibile
|
|
lstTria.emplace_back() ;
|
|
lstTria.back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_MapFrame.VersZ()) ;
|
|
lstTria.emplace_back() ;
|
|
lstTria.back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_MapFrame.VersZ()) ;
|
|
// faccia inferiore P1b->P4b->P3b->P2b : sempre visibile
|
|
lstTria.emplace_back() ;
|
|
lstTria.back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_MapFrame.VersZ()) ;
|
|
lstTria.emplace_back() ;
|
|
lstTria.back().Set( ptP3 + vtDZb, ptP2 + vtDZb, ptP1 + vtDZb, - m_MapFrame.VersZ()) ;
|
|
}
|
|
|
|
// creo le facce laterali
|
|
int nPosEst = ( nPos1 < int( m_nNx[0] - 1) ? nPos + 1 : - 1) ;
|
|
AddDexelSideFace( nPos, nPosEst, ptP2, ptP3, m_MapFrame.VersZ(), m_MapFrame.VersX(), lstTria) ;
|
|
int nPosNord = ( nPos2 < int( m_nNy[0] - 1) ? nPos + m_nNx[0] : - 1) ;
|
|
AddDexelSideFace( nPos, nPosNord, ptP3, ptP4, m_MapFrame.VersZ(), m_MapFrame.VersY(), lstTria) ;
|
|
int nPosWest = ( nPos1 > 0 ? nPos - 1 : - 1) ;
|
|
AddDexelSideFace( nPos, nPosWest, ptP4, ptP1, m_MapFrame.VersZ(), - m_MapFrame.VersX(), lstTria) ;
|
|
int nPosSud = ( nPos2 > 0 ? nPos - m_nNx[0] : - 1) ;
|
|
AddDexelSideFace( nPos, nPosSud, ptP1, ptP2, m_MapFrame.VersZ(), - m_MapFrame.VersY(), lstTria) ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Point3d& ptQ,
|
|
const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DLIST& lstTria) const
|
|
{
|
|
Intervals intFace ;
|
|
for ( int i = 0 ; i < int( m_Values[0][nPos].size()) ; i += 1)
|
|
intFace.Add( m_Values[0][nPos][i].dMin, m_Values[0][nPos][i].dMax) ;
|
|
if ( nPosAdj > 0) {
|
|
for ( int i = 0 ; i < int( m_Values[0][nPosAdj].size()) ; i += 1)
|
|
intFace.Subtract( m_Values[0][nPosAdj][i].dMin, m_Values[0][nPosAdj][i].dMax) ;
|
|
}
|
|
double dMin, dMax ;
|
|
bool bFound = intFace.GetFirst( dMin, dMax) ;
|
|
while ( bFound) {
|
|
Vector3d vtDZt = dMax * vtZ ;
|
|
Vector3d vtDZb = dMin * vtZ ;
|
|
lstTria.emplace_back() ;
|
|
lstTria.back().Set( ptP + vtDZb, ptQ + vtDZb, ptQ + vtDZt, vtNorm) ;
|
|
lstTria.emplace_back() ;
|
|
lstTria.back().Set( ptQ + vtDZt, ptP + vtDZt, ptP + vtDZb, vtNorm) ;
|
|
bFound = intFace.GetNext( dMin, dMax) ;
|
|
}
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Calcola i triangoli del voxel con indici nVoxI, nVoxJ, nVoxK; se si vuole
|
|
// il riconoscimento delle sharp-feature bEnh deve valere true.
|
|
// I triangoli formanti sharp-feature vengono messi nel TriHolder, gli altri nella lista.
|
|
bool
|
|
VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DLIST& lstTria, TriHolder& triHold, bool bEnh) const
|
|
{
|
|
// Calcolo il numero di voxel lungo X,Y e Z
|
|
int nVoxNumX = int( m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] / N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
|
|
// Se il voxel non esiste, vi è un errore.
|
|
if ( nVoxI + 1 < 0 || nVoxI + 2 > nVoxNumX ||
|
|
nVoxJ + 1 < 0 || nVoxJ + 2 > nVoxNumY ||
|
|
nVoxK + 1 < 0 || nVoxK + 2 > nVoxNumZ)
|
|
return false ;
|
|
|
|
// Classificazione dei vertici: interni o esterni al materiale
|
|
int nIndex = CalcIndex( nVoxI, nVoxJ , nVoxK) ;
|
|
|
|
// Se vi è qualche intersezione fra segmenti e superficie
|
|
// continuo altrimenti passo al prossimo voxel.
|
|
if ( EdgeTable[nIndex] != 0) {
|
|
|
|
// Indici i,j,k dei vertici
|
|
int IndexCorner[8][3] = {
|
|
{ nVoxI, nVoxJ, nVoxK},
|
|
{ nVoxI + 1, nVoxJ, nVoxK},
|
|
{ nVoxI + 1, nVoxJ + 1, nVoxK},
|
|
{ nVoxI, nVoxJ + 1, nVoxK},
|
|
{ nVoxI, nVoxJ, nVoxK + 1},
|
|
{ nVoxI + 1, nVoxJ, nVoxK + 1},
|
|
{ nVoxI + 1, nVoxJ + 1, nVoxK + 1},
|
|
{ nVoxI, nVoxJ + 1, nVoxK + 1}
|
|
} ;
|
|
|
|
static int intersections[12][2] = {
|
|
{ 0, 1 }, { 1, 2 }, { 3, 2 }, { 0, 3 }, { 4, 5 }, { 5, 6 },
|
|
{ 7, 6 }, { 4, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }
|
|
} ;
|
|
|
|
// Array di strutture punto di intersezione e normale alla superficie in esso.
|
|
VectorField VecField[12] ;
|
|
|
|
// Flag di regolarità dei campi scalare e vettoriale
|
|
bool bReg = true ;
|
|
|
|
// Ciclo sui segmenti
|
|
for ( int EdgeIndex = 0 ; EdgeIndex < 12 ; ++ EdgeIndex) {
|
|
// Se il segmento non attraversa la superficie passo al successivo
|
|
if ( ( EdgeTable[nIndex] & ( 1 << EdgeIndex)) == 0)
|
|
continue ;
|
|
// Indici per linee di griglia sui vertici
|
|
int n1 = intersections[EdgeIndex][0] ;
|
|
int n2 = intersections[EdgeIndex][1] ;
|
|
// Flag posizione corner
|
|
bool bN1 = ( ( nIndex & ( 1 << n1)) != 0) ;
|
|
// Determino con precisione il punto di intersezione sullo spigolo,
|
|
// se i campi scalare e vettoriale non sono regolari bReg diviene falso.
|
|
if ( ! IntersPos( IndexCorner[n1], IndexCorner[n2], bN1,
|
|
VecField[EdgeIndex].ptInt, VecField[EdgeIndex].vtNorm))
|
|
bReg = false ;
|
|
}
|
|
|
|
// Determino il numero di componenti connesse nel voxel in caso di configurazione standard
|
|
int nComponents = TriangleTableEn[nIndex][1][0] ;
|
|
|
|
// Matrici di campi vettoriali:
|
|
// CompoVert[i] ha i vertici della base del triangle fan della (i+1)-esima componente connessa;
|
|
// CompoTriVert[i] ha i vertici di tutti i triangoli, nel caso di assenza di sharp feature,
|
|
// della (i+1)-esima componente connessa.
|
|
VectorField CompoVert[6][12] ;
|
|
VectorField CompoTriVert[6][17] ;
|
|
|
|
// Arrey numero di vertici della base del fan per componente
|
|
// connessa: nVertComp[i] contiene il numero di vertici
|
|
// della base del fan della (i+1)-esima componente connessa.
|
|
int nVertComp[6] ;
|
|
|
|
// Matrice di indici dei punti: serve per la gestione del caso
|
|
int nIndArrey[6][4] ;
|
|
|
|
int nExtTabOff = nComponents ;
|
|
int nStdTabOff = 0 ;
|
|
|
|
// Carico le matrici CompoVert e CompoTriVert
|
|
for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) {
|
|
|
|
// Numero vertici per componenti
|
|
nVertComp[nCompCount - 1] = TriangleTableEn[nIndex][1][nCompCount] ;
|
|
|
|
// Riempio il nCompCount-esimo vettore di vertici della base del fan
|
|
for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount)
|
|
CompoVert[nCompCount - 1][nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1]] ;
|
|
|
|
// Serve per la gestione del caso ...
|
|
if ( nVertComp[nCompCount - 1] == 4) {
|
|
for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount)
|
|
nIndArrey[nCompCount - 1][nVertCount] = TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1] ;
|
|
}
|
|
|
|
// Riempio il nCompCount-esimo vettore di vertici dei triangoli in assenza di
|
|
// sharp feature: in una mesh di triangoli con n vertici vi sono n - 2 triangoli.
|
|
for ( int nVertCount = 0 ; nVertCount < 3 * ( nVertComp[nCompCount - 1] - 2) ; nVertCount += 3) {
|
|
CompoTriVert[nCompCount - 1][nVertCount] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+2]] ;
|
|
CompoTriVert[nCompCount - 1][nVertCount+1] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+1]] ;
|
|
CompoTriVert[nCompCount - 1][nVertCount+2] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount]] ;
|
|
}
|
|
|
|
// Aggiorno gli offsets per raggiungere i
|
|
// vertici della componente successiva.
|
|
nExtTabOff += nVertComp[nCompCount - 1] ;
|
|
nStdTabOff += 3 * ( nVertComp[nCompCount - 1] - 2) ;
|
|
}
|
|
|
|
// Test sulla topologia: dal momento che il nostro test
|
|
// si fonda sugli angoli compresi fra le normali, esso ha
|
|
// senso solo se il campo è regolare.
|
|
// Il flag bSecondConfig6 serve perché nel caso di configurazione 6
|
|
// con una sola componente connessa non si deve eseguire l'Extended MC.
|
|
bool bSecondConfig6 = false ;
|
|
|
|
if ( bReg) {
|
|
if ( nAllConfig[nIndex] == 3) {
|
|
|
|
Vector3d vtCmpAvg0, vtCmpAvg1 ;
|
|
|
|
// Verifico se i versori delle componenti sono tutti
|
|
// più o meno concordi (per esserlo non devono esserci
|
|
// due vettori di una medesima componente con prodotto
|
|
// scalare inferiore a 0.7).
|
|
bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ;
|
|
bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
|
|
|
|
// Se i versori di entrambe le componenti sono concordi
|
|
// ha senso parlare di vettori medi, altrimenti non ha
|
|
// senso. Se non ha senso parlare di vettori medi non
|
|
// ha senso parlare di prodotti scalari fra loro,
|
|
// quindi pongo il loro prodotto a un valore assurdo -2
|
|
// (il prodotto scalare fra versori ha modulo non superiore
|
|
// a uno).
|
|
double dScProd = - 2 ;
|
|
|
|
if ( bTest0 && bTest1)
|
|
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
|
|
|
|
double dThreshold = 0.7 ;
|
|
|
|
if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
|
|
|
|
int nt = 0 ;
|
|
|
|
while ( nIndexVsIndex3[nt][0] != nIndex)
|
|
++ nt ;
|
|
|
|
int nRotCase = nIndexVsIndex3[nt][1] ;
|
|
|
|
nComponents = Cases3Plus[nRotCase][1][0] ;
|
|
|
|
// Riaggiorno gli offsets
|
|
nExtTabOff = nComponents ;
|
|
nStdTabOff = 0 ;
|
|
|
|
// Modifico le matrici
|
|
for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
|
|
|
|
// Numero vertici per componenti
|
|
nVertComp[nC - 1] = Cases3Plus[nRotCase][1][nC] ;
|
|
|
|
// Matrice dei vertici della base del fan
|
|
for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert)
|
|
CompoVert[nC - 1][nFanVert] = VecField[Cases3Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ;
|
|
|
|
// Matrici dei vertici dei triangoli in assenza di sharp feature
|
|
for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) {
|
|
|
|
CompoTriVert[nC - 1][nTriVert] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
|
|
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
|
|
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
|
|
}
|
|
|
|
// Aggiorno gli offsets per raggiungere i
|
|
// vertici della componente successiva.
|
|
nExtTabOff += nVertComp[nC - 1] ;
|
|
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
|
|
}
|
|
}
|
|
}
|
|
|
|
else if ( nAllConfig[nIndex] == 6) {
|
|
|
|
// Procedura analoga a quella della configurazione 3
|
|
Vector3d vtCmpAvg0, vtCmpAvg1 ;
|
|
|
|
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.7) ;
|
|
bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
|
|
|
|
double dScProd = - 2 ;
|
|
|
|
if ( bTest0 && bTest1)
|
|
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
|
|
|
|
double dThreshold = 0.7 ;
|
|
|
|
if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
|
|
|
|
int nt = 0 ;
|
|
|
|
while ( nIndexVsIndex6[nt][0] != nIndex)
|
|
++ nt ;
|
|
|
|
int nRotCase = nIndexVsIndex6[nt][1] ;
|
|
|
|
nComponents = 1 ;
|
|
|
|
// Riaggiorno gli offsets
|
|
nExtTabOff = nComponents ;
|
|
nStdTabOff = 0 ;
|
|
|
|
// Modifico le matrici
|
|
for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
|
|
|
|
// Numero vertici per componenti
|
|
nVertComp[nC - 1] = 7 ;
|
|
|
|
// Matrici dei vertici dei triangoli in assenza di sharp feature
|
|
for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) {
|
|
|
|
CompoTriVert[nC - 1][nTriVert] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
|
|
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
|
|
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
|
|
}
|
|
|
|
// Aggiorno gli offsets per raggiungere i
|
|
// vertici della componente successiva.
|
|
nExtTabOff += nVertComp[nC - 1] ;
|
|
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
|
|
}
|
|
bSecondConfig6 = true ;
|
|
}
|
|
}
|
|
|
|
else if ( nAllConfig[nIndex] == 10) {
|
|
|
|
Vector3d vtCmpAvg0, vtCmpAvg1 ;
|
|
|
|
// Verifico se i versori delle componenti sono tutti
|
|
// più o meno concordi (per esserlo non devono esserci
|
|
// due vettori di una medesima componente con prodotto
|
|
// scalare inferiore a 0). decidere se 0.0 o 0.7
|
|
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ;
|
|
bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ;
|
|
|
|
if ( ! ( bTest0 && bTest1)) {
|
|
|
|
int nt = 0 ;
|
|
while ( nIndexVsIndex10[nt][0] != nIndex)
|
|
++ nt ;
|
|
|
|
int nRotCase = nIndexVsIndex10[nt][1] ;
|
|
|
|
// Riaggiorno gli offsets
|
|
nExtTabOff = 2 ;
|
|
nStdTabOff = 0 ;
|
|
|
|
// Modifico le matrici
|
|
for ( int nC = 1 ; nC <= 2 ; ++ nC) {
|
|
// Numero vertici per componenti
|
|
nVertComp[nC - 1] = Cases10Plus[nRotCase][1][nC] ;
|
|
|
|
// Matrice dei vertici della base del fan
|
|
for ( int nFanVert = 0 ; nFanVert < 4 ; ++ nFanVert)
|
|
CompoVert[nC - 1][nFanVert] = VecField[Cases10Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ;
|
|
|
|
// Matrici dei vertici dei triangoli in assenza di sharp feature
|
|
for ( int nTriVert = 0 ; nTriVert < 6 ; nTriVert += 3) {
|
|
CompoTriVert[nC - 1][nTriVert] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
|
|
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
|
|
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
|
|
}
|
|
|
|
// Aggiorno gli offsets per raggiungere i vertici della componente successiva.
|
|
nExtTabOff += nVertComp[nC - 1] ;
|
|
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Numero di feature nel voxel: al più vi è una feature per componente connessa.
|
|
int nInnerFeatureInVoxel = 0 ;
|
|
int nBorderFeatureInVoxel = 0 ;
|
|
|
|
// Ciclo sulle componenti
|
|
for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) {
|
|
|
|
int nFeatureType = NO_FEATURE ;
|
|
// Se i componenti sono regolari valuto le normali per stabilire se eseguire ExtMC o MC
|
|
if ( bReg)
|
|
nFeatureType = TestOnNormal( CompoVert[nCompCount - 1], nVertComp[nCompCount - 1]) ;
|
|
|
|
// Controllo per il caso piano su una griglia
|
|
// con versori normali a due a due paralleli.
|
|
bool bGridControl = true ;
|
|
if ( nFeatureType != NO_FEATURE) {
|
|
if ( nVertComp[nCompCount - 1] == 4)
|
|
GridControl( VecField, CompoVert, nIndArrey, nVertComp, nCompCount, bGridControl) ;
|
|
}
|
|
|
|
// Flag ExtMC
|
|
bool bExtMC = ( nFeatureType != NO_FEATURE && bGridControl) && ( ! bSecondConfig6) && bEnh ;
|
|
|
|
// Extended MC
|
|
if ( bExtMC) {
|
|
|
|
// Passo al sistema di riferimento del baricentro
|
|
Point3d ptGravityCenter( 0, 0, 0) ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni)
|
|
ptGravityCenter = ptGravityCenter + CompoVert[nCompCount - 1][ni].ptInt ;
|
|
ptGravityCenter = ptGravityCenter / nVertComp[nCompCount - 1] ;
|
|
|
|
Vector3d vtTrasf[12] ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni)
|
|
vtTrasf[ni] = CompoVert[nCompCount - 1][ni].ptInt - ptGravityCenter ;
|
|
|
|
// Preparo le matrici per il sistema
|
|
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> dSystemMatrix ;
|
|
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> dSystemVector ;
|
|
typedef Eigen::JacobiSVD<dSystemMatrix> DecomposerSVD ;
|
|
|
|
dSystemMatrix dMatrixN( nVertComp[nCompCount - 1], 3) ;
|
|
dSystemVector dKnownVector( nVertComp[nCompCount - 1], 1) ;
|
|
|
|
// Caso generale
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ;
|
|
dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ;
|
|
dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ;
|
|
dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * vtTrasf[ni] ;
|
|
}
|
|
|
|
// calcolo SVD
|
|
DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ;
|
|
dSystemMatrix dMatrixV = svd.matrixV() ;
|
|
dSystemVector dSingularValue = svd.singularValues() ;
|
|
|
|
// Se la feature è un edge annullo il valore singolare minore.
|
|
if ( nFeatureType == EDGE) {
|
|
double dThres = 0.5 * ( dSingularValue( 1) + dSingularValue( 2)) / dSingularValue( 0) ;
|
|
svd.setThreshold( dThres) ;
|
|
}
|
|
|
|
// risolvo il sistema con SVD, quindi ai minimi quadrati
|
|
dSystemVector dUnknownVector( 3, 1) ;
|
|
dUnknownVector = svd.solve( dKnownVector) ;
|
|
|
|
// Vettore Baricentro-Feature
|
|
Vector3d vtFeature( dUnknownVector( 0), dUnknownVector( 1), dUnknownVector( 2)) ;
|
|
|
|
// Esprimo la soluzione nel sistema di riferimento z-map
|
|
Point3d ptSol = ptGravityCenter + vtFeature ;
|
|
|
|
// Limito la feature entro i secondi voxel vicini,
|
|
// nel caso essa esca dal reticolo la limito entro una
|
|
// distanza dal baricentro pari alla diagonale del voxel.
|
|
bool bOutside = false ;
|
|
|
|
int nFtVxI, nFtVxJ, nFtVxK ;
|
|
if ( GetPointVoxel( ptSol, nFtVxI, nFtVxJ, nFtVxK)) {
|
|
|
|
if ( abs( nFtVxI - nVoxI) > 2 ||
|
|
abs( nFtVxJ - nVoxJ) > 2 ||
|
|
abs( nFtVxK - nVoxK) > 2)
|
|
bOutside = true ;
|
|
}
|
|
else {
|
|
double dDistFeature = vtFeature.Len() ;
|
|
const double MAX_DIST = sqrt( 3) * N_DEXVOXRATIO * m_dStep ;
|
|
bOutside = ( dDistFeature > MAX_DIST) ;
|
|
}
|
|
|
|
TRIA3DVECTOR triContainer ;
|
|
|
|
// Costruisco triangoli di prova
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
|
|
// Il triangolo è pronto
|
|
Triangle3d CurrentTriangle ;
|
|
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
// Aggiungo triangolo al vettore temporaneo
|
|
triContainer.emplace_back( CurrentTriangle) ;
|
|
}
|
|
|
|
// Controllo delle inversioni dei triangoli
|
|
bool bDangerInversion = false ;
|
|
|
|
// Caso ventaglio con tre vertici di base
|
|
if ( nVertComp[nCompCount - 1] == 3) {
|
|
|
|
// Controllo se esiste almeno un triangolo con normale avente prodotto scalare
|
|
// negativo con la normale di almeno uno dei vertici di base del ventaglio.
|
|
bool bInversione = false ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
|
|
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
|
|
|
|
double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
|
|
double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ;
|
|
|
|
if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) {
|
|
bInversione = true ;
|
|
break ;
|
|
}
|
|
}
|
|
|
|
// Se tale triangolo esiste procedo
|
|
if ( bInversione) {
|
|
|
|
// Conto le coppie di normali con angolo compreso maggiore di 90°
|
|
int nNegDot = 0 ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] - 1 ; ++ ni) {
|
|
for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) {
|
|
if ( CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL)
|
|
nNegDot ++ ;
|
|
}
|
|
}
|
|
|
|
if ( nNegDot == nVertComp[nCompCount - 1] - 1) {
|
|
// Cerco se esiste un punto in cui la normale ha prodotto scalre negativo
|
|
// con le normali di entrambi i triangoli che lo contengono
|
|
bool bInversione2 = false ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
int nj = ( ni == 0 ? nVertComp[nCompCount - 1] - 1 : ni - 1) ;
|
|
double dDLast = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
|
|
double dDPrev = triContainer[nj].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
|
|
if ( ( dDLast < EPS_SMALL && dDPrev < EPS_SMALL) ||
|
|
( dDLast < - 0.75 || dDPrev < - 0.75)) {
|
|
bInversione2 = true ;
|
|
break ;
|
|
}
|
|
}
|
|
|
|
// Se tale normale esiste
|
|
if ( bInversione2) {
|
|
|
|
// Soluzione del sistema nel sistema Zmap
|
|
Point3d ptSolZMapFrame = ptSol ;
|
|
|
|
// Se la soluzione non cade nel voxel di appartenenza vedo se può
|
|
// essere riportata dentro muovendosi lungo la linea di feature.
|
|
if ( ! IsPointInsideVoxelApprox( nVoxI, nVoxJ, nVoxK, ptSolZMapFrame)) {
|
|
|
|
Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
|
|
double dParInt1, dParInt2 ;
|
|
Point3d ptVoxMin( ( nVoxI * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nVoxJ * N_DEXVOXRATIO + 0.5) * m_dStep, ( nVoxK * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
Point3d ptVoxMax( ( ( nVoxI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( nVoxJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( nVoxK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
// Caso in cui può essere riportata dentro: se il voxel in cui cade la feature è pieno
|
|
// la riporto muovendola lungo la sua linea
|
|
if ( IntersLineBox( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) {
|
|
|
|
triContainer.resize( 0) ;
|
|
|
|
double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 :
|
|
dParInt2 + ( dParInt1 - dParInt2) / 100 ;
|
|
|
|
Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ;
|
|
|
|
ptSol = ptNewSol ;
|
|
|
|
// Costruisco triangoli di prova
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
|
|
// Il triangolo è pronto
|
|
Triangle3d CurrentTriangle ;
|
|
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
// Aggiungo triangolo al vettore temporaneo
|
|
triContainer.emplace_back( CurrentTriangle) ;
|
|
}
|
|
}
|
|
|
|
// Caso in cui non può essere riportata dentro:
|
|
// si passa alla routine standard se il voxel
|
|
// in cui cade è pieno.
|
|
else {
|
|
|
|
int nAdjVoxI, nAdjVoxJ, nAdjVoxK ;
|
|
if ( GetPointVoxel( ptSolZMapFrame, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) {
|
|
|
|
// Classificazione del voxel adiacente
|
|
int nAdjIndex = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ;
|
|
|
|
// Se il voxel è pieno
|
|
if ( EdgeTable[nAdjIndex] != 0)
|
|
bDangerInversion = true ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Ventaglio con base a quattro vertici
|
|
else if ( nVertComp[nCompCount - 1] == 4) {
|
|
// Controllo se esiste almeno un triangolo con normale avente prodotto scalare
|
|
// negativo con la normale di almeno uno dei vertici di base del ventaglio.
|
|
bool bInversione = false ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
|
|
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
|
|
|
|
double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
|
|
double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ;
|
|
|
|
if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL)
|
|
bInversione = true ;
|
|
}
|
|
// Se tale triangolo esiste continuo i test
|
|
if ( bInversione) {
|
|
|
|
// Se la feature non cade nel suo voxel
|
|
if ( ! IsPointInsideVoxelApprox( nVoxI, nVoxJ, nVoxK, ptSol)) {
|
|
|
|
Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
|
|
|
|
double dParInt1, dParInt2 ;
|
|
Point3d ptVoxMin( ( nVoxI * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nVoxJ * N_DEXVOXRATIO + 0.5) * m_dStep, ( nVoxK * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
Point3d ptVoxMax( ( ( nVoxI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( nVoxJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( nVoxK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
// Se è possibile riportarla dentro e il voxel in cui cade è pieno, la riporto nel suo voxel
|
|
// lungo la sua linea
|
|
if ( IntersLineBox( ptSol, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) {
|
|
|
|
triContainer.resize( 0) ;
|
|
|
|
double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 :
|
|
dParInt2 + ( dParInt1 - dParInt2) / 100 ;
|
|
|
|
Point3d ptNewSol = ptSol + dPar * vtNullSpace ;
|
|
|
|
ptSol = ptNewSol ;
|
|
|
|
// Costruisco triangoli di prova
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
|
|
// Il triangolo è pronto
|
|
Triangle3d CurrentTriangle ;
|
|
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
// Aggiungo triangolo al vettore temporaneo
|
|
triContainer.emplace_back( CurrentTriangle) ;
|
|
}
|
|
}
|
|
|
|
// Se non è possibile riportarla dentro e il voxel in
|
|
// cui cade è pieno passo alla routine standard
|
|
else {
|
|
int nAdjVoxI, nAdjVoxJ, nAdjVoxK ;
|
|
if ( GetPointVoxel( ptSol, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) {
|
|
// Classificazione del voxel adiacente
|
|
int nAdjIndex = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ;
|
|
// Se il voxel è pieno
|
|
if ( EdgeTable[nAdjIndex] != 0)
|
|
bDangerInversion = true ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Valuto normali: questo è ancora un controllo
|
|
// sulle normali, se risultano in tutti i punti
|
|
// approssimativamente uguali passiamo alla
|
|
// routine standard
|
|
int nContSize = int( triContainer.size()) ;
|
|
|
|
bool bPlane = true ;
|
|
|
|
for ( int ni = 0 ; ni < nContSize - 1 ; ++ ni) {
|
|
for ( int nj = ni + 1 ; nj < nContSize ; ++ nj) {
|
|
|
|
Vector3d vtI = triContainer[ni].GetN() ;
|
|
Vector3d vtJ = triContainer[nj].GetN() ;
|
|
|
|
if ( ! AreSameVectorApprox( vtI, vtJ)) {
|
|
bPlane = false ;
|
|
break ;
|
|
}
|
|
}
|
|
|
|
if ( ! bPlane)
|
|
break ;
|
|
}
|
|
|
|
// Se feature nei limiti e triangoli non in un piano, confermo ExtMC
|
|
if ( ! bOutside && ! bPlane && ! bDangerInversion) {
|
|
|
|
TRIA3DVECTOR vInnerTriaTemp ;
|
|
|
|
for ( int ni = 0 ; ni < nContSize ; ++ ni) {
|
|
|
|
// Se l'area è non nulla aggiungo il triangolo a uno dei due vettori.
|
|
if ( triContainer[ni].GetArea() > SQ_EPS_SMALL) {
|
|
int nVoxIJK[3] = { nVoxI, nVoxJ, nVoxK} ;
|
|
Point3d ptTemp = ptSol ;
|
|
Triangle3d trTemp = triContainer[ni] ;
|
|
vInnerTriaTemp.emplace_back( triContainer[ni]) ;
|
|
}
|
|
}
|
|
|
|
if ( vInnerTriaTemp.size() > 0) {
|
|
|
|
// Aggiorno il numero di feature.
|
|
++ nInnerFeatureInVoxel ;
|
|
|
|
// Se siamo in presenza della prima feature del
|
|
// voxel, ridimensiono il vettore che contiene
|
|
// la struttura dati del voxel.
|
|
if ( nInnerFeatureInVoxel == 1) {
|
|
triHold.resize( triHold.size() + 1) ;
|
|
// Questi dati dipendono solo dal voxel
|
|
int nCurrent = int( triHold.size()) - 1 ;
|
|
triHold[nCurrent].i = nVoxI ;
|
|
triHold[nCurrent].j = nVoxJ ;
|
|
triHold[nCurrent].k = nVoxK ;
|
|
}
|
|
|
|
// Indice che identifica la posizione del voxel nel vector.
|
|
int nCurrent = int( triHold.size()) - 1 ;
|
|
|
|
// Aggiungo vertice della componente connessa alla lista dei vertici.
|
|
triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ;
|
|
|
|
int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ;
|
|
int nNewFeatureNum = nOldFeatureNum + 1 ;
|
|
|
|
// Aggiungo una componente per il vettore dei triangoli della componente connessa.
|
|
triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ;
|
|
for ( int ni = 0 ; ni < int( vInnerTriaTemp.size()) ; ++ ni) {
|
|
triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vInnerTriaTemp[ni]) ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// ExtMC non confermato, si passa a MC
|
|
else {
|
|
// Costruzione dei triangoli
|
|
for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {
|
|
// Il triangolo è pronto
|
|
Triangle3d CurrentTriangle ;
|
|
CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt,
|
|
CompoTriVert[nCompCount - 1][TriIndex+1].ptInt,
|
|
CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ;
|
|
bool bV = CurrentTriangle.Validate( true) ;
|
|
// Aggiungo alla lista
|
|
lstTria.emplace_back( CurrentTriangle) ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Standard MC
|
|
else {
|
|
// Costruzione dei triangoli
|
|
for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {
|
|
// Il triangolo è pronto
|
|
Triangle3d CurrentTriangle ;
|
|
CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt,
|
|
CompoTriVert[nCompCount - 1][TriIndex+1].ptInt,
|
|
CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ;
|
|
bool bV = CurrentTriangle.Validate( true) ;
|
|
// Aggiungo alla lista
|
|
lstTria.emplace_back( CurrentTriangle) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) const
|
|
{
|
|
// Controllo sulla validità del blocco
|
|
if ( nBlock < 0 || nBlock >= int( m_nNumBlock))
|
|
return false ;
|
|
|
|
// Calcolo i limiti sugli indici dei voxel del blocco
|
|
int nIJK[3] ;
|
|
// Vettore indici i,j,k del blocco
|
|
GetBlockIJKFromN( nBlock, nIJK) ;
|
|
// Vettore limiti sugli indici dei voxel del blocco
|
|
int nLimits[6] ;
|
|
GetBlockLimitsIJK( nIJK, nLimits) ;
|
|
|
|
// Unordered Map per la riduzione del numero di triangoli
|
|
int nDim = m_nVoxNumPerBlock * m_nVoxNumPerBlock ;
|
|
VoxelContainer VoxContXZInf( nDim) ;
|
|
VoxelContainer VoxContXZSup( nDim) ;
|
|
VoxelContainer VoxContXYInf( nDim) ;
|
|
VoxelContainer VoxContXYSup( nDim) ;
|
|
VoxelContainer VoxContYZInf( nDim) ;
|
|
VoxelContainer VoxContYZSup( nDim) ;
|
|
|
|
// Ciclo su tutti i voxel dello Zmap
|
|
for ( int i = nLimits[0] ; i < nLimits[1] ; ++ i) {
|
|
for ( int j = nLimits[2] ; j < nLimits[3] ; ++ j) {
|
|
for ( int k = nLimits[4] ; k < nLimits[5] ; ++ k) {
|
|
|
|
// Classificazione dei vertici: interni o esterni al materiale
|
|
int nIndex = CalcIndex( i, j, k) ;
|
|
|
|
// Se vi è qualche intersezione fra segmenti e superficie
|
|
// continuo altrimenti passo al prossimo voxel.
|
|
if ( EdgeTable[nIndex] == 0)
|
|
continue ;
|
|
|
|
// Indici i,j,k dei vertici
|
|
int IndexCorner[8][3] = {
|
|
{ i, j, k},
|
|
{ i + 1, j, k},
|
|
{ i + 1, j + 1, k},
|
|
{ i, j + 1, k},
|
|
{ i, j, k + 1},
|
|
{ i + 1, j, k + 1},
|
|
{ i + 1, j + 1, k + 1},
|
|
{ i, j + 1, k + 1}
|
|
} ;
|
|
|
|
static int intersections[12][2] = {
|
|
{ 0, 1 }, { 1, 2 }, { 3, 2 }, { 0, 3 }, { 4, 5 }, { 5, 6 },
|
|
{ 7, 6 }, { 4, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }
|
|
} ;
|
|
|
|
// Array di strutture punto di intersezione e normale alla superficie in esso.
|
|
VectorField VecField[12] ;
|
|
|
|
// Flag di regolarità dei campi scalare e vettoriale
|
|
bool bReg = true ;
|
|
|
|
// Ciclo sui segmenti
|
|
for ( int EdgeIndex = 0 ; EdgeIndex < 12 ; ++ EdgeIndex) {
|
|
// Se il segmento non attraversa la superficie passo al successivo
|
|
if ( ( EdgeTable[nIndex] & ( 1 << EdgeIndex)) == 0)
|
|
continue ;
|
|
// Indici per linee di griglia sui vertici
|
|
int n1 = intersections[EdgeIndex][0] ;
|
|
int n2 = intersections[EdgeIndex][1] ;
|
|
// Flag posizione corner
|
|
bool bN1 = ( ( nIndex & ( 1 << n1)) != 0) ;
|
|
// Determino con precisione il punto di intersezione sullo spigolo,
|
|
// se i campi scalare e vettoriale non sono regolari bReg diviene falso.
|
|
if ( ! IntersPos( IndexCorner[n1], IndexCorner[n2], bN1,
|
|
VecField[EdgeIndex].ptInt, VecField[EdgeIndex].vtNorm))
|
|
bReg = false ;
|
|
}
|
|
|
|
// Determino il numero di componenti connesse nel voxel in caso di configurazione standard
|
|
int nComponents = TriangleTableEn[nIndex][1][0] ;
|
|
|
|
// Matrici di campi vettoriali:
|
|
// CompoVert[i] ha i vertici della base del triangle fan della (i+1)-esima componente connessa;
|
|
// CompoTriVert[i] ha i vertici di tutti i triangoli, nel caso di assenza di sharp feature,
|
|
// della (i+1)-esima componente connessa.
|
|
VectorField CompoVert[6][12] ;
|
|
VectorField CompoTriVert[6][17] ;
|
|
|
|
// Arrey numero di vertici della base del fan per componente
|
|
// connessa: nVertComp[i] contiene il numero di vertici
|
|
// della base del fan della (i+1)-esima componente connessa.
|
|
int nVertComp[6] ;
|
|
|
|
// Matrice di indici dei punti: serve per la gestione del caso
|
|
int nIndArrey[6][4] ;
|
|
|
|
int nExtTabOff = nComponents ;
|
|
int nStdTabOff = 0 ;
|
|
|
|
// Carico le matrici CompoVert e CompoTriVert
|
|
for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) {
|
|
|
|
// Numero vertici per componenti
|
|
nVertComp[nCompCount - 1] = TriangleTableEn[nIndex][1][nCompCount] ;
|
|
|
|
// Riempio il nCompCount-esimo vettore di vertici della base del fan
|
|
for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount)
|
|
CompoVert[nCompCount - 1][nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1]] ;
|
|
|
|
// Serve per la gestione del caso ...
|
|
if ( nVertComp[nCompCount - 1] == 4) {
|
|
for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount)
|
|
nIndArrey[nCompCount - 1][nVertCount] = TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1] ;
|
|
}
|
|
|
|
// Riempio il nCompCount-esimo vettore di vertici dei triangoli in assenza di
|
|
// sharp feature: in una mesh di triangoli con n vertici vi sono n - 2 triangoli.
|
|
for ( int nVertCount = 0 ; nVertCount < 3 * ( nVertComp[nCompCount - 1] - 2) ; nVertCount += 3) {
|
|
CompoTriVert[nCompCount - 1][nVertCount] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+2]] ;
|
|
CompoTriVert[nCompCount - 1][nVertCount+1] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+1]] ;
|
|
CompoTriVert[nCompCount - 1][nVertCount+2] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount]] ;
|
|
}
|
|
|
|
// Aggiorno gli offsets per raggiungere i
|
|
// vertici della componente successiva.
|
|
nExtTabOff += nVertComp[nCompCount - 1] ;
|
|
nStdTabOff += 3 * ( nVertComp[nCompCount - 1] - 2) ;
|
|
}
|
|
|
|
// Controllo se il voxel ha una sola faccia che giace in un piano canonico e quindi ha gestione speciale
|
|
// Faccia XY normale Z+
|
|
if ( nIndex == 15) {
|
|
double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], Z_PLUS, dPos)) {
|
|
int nN ;
|
|
GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContXYSup.emplace( nN, dPos) ;
|
|
continue ;
|
|
}
|
|
}
|
|
// Faccia YZ normale X+
|
|
else if ( nIndex == 153) {
|
|
double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], X_PLUS, dPos)) {
|
|
int nN ;
|
|
GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContYZSup.emplace( nN, dPos) ;
|
|
continue ;
|
|
}
|
|
}
|
|
// Faccia ZX normale Y+
|
|
else if ( nIndex == 51) {
|
|
double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], Y_PLUS, dPos)) {
|
|
int nN ;
|
|
GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContXZSup.emplace( nN, dPos) ;
|
|
continue ;
|
|
}
|
|
}
|
|
// Faccia YX normale Z-
|
|
else if ( nIndex == 240) {
|
|
double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], Z_MINUS, dPos)) {
|
|
int nN ;
|
|
GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContXYInf.emplace( nN, dPos) ;
|
|
continue ;
|
|
}
|
|
}
|
|
// Faccia ZY normale X-
|
|
else if ( nIndex == 102) {
|
|
double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], X_MINUS, dPos)) {
|
|
int nN ;
|
|
GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContYZInf.emplace( nN, dPos) ;
|
|
continue ;
|
|
}
|
|
}
|
|
// Faccia XZ normale Y-
|
|
else if ( nIndex == 204) {
|
|
double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], Y_MINUS, dPos)) {
|
|
int nN ;
|
|
GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContXZInf.emplace( nN, dPos) ;
|
|
continue ;
|
|
}
|
|
}
|
|
|
|
// Test sulla topologia: dal momento che il nostro test
|
|
// si fonda sugli angoli compresi fra le normali, esso ha
|
|
// senso solo se il campo è regolare.
|
|
if ( bReg) {
|
|
if ( nAllConfig[nIndex] == 3) {
|
|
|
|
Vector3d vtCmpAvg0, vtCmpAvg1 ;
|
|
|
|
// Verifico se i versori delle componenti sono tutti
|
|
// più o meno concordi (per esserlo non devono esserci
|
|
// due vettori di una medesima componente con prodotto
|
|
// scalare inferiore a 0.7).
|
|
bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ;
|
|
bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
|
|
|
|
// Se i versori di entrambe le componenti sono concordi
|
|
// ha senso parlare di vettori medi, altrimenti non ha
|
|
// senso. Se non ha senso parlare di vettori medi non
|
|
// ha senso parlare di prodotti scalari fra loro,
|
|
// quindi pongo il loro prodotto a un valore assurdo -2
|
|
// (il prodotto scalare fra versori ha modulo non superiore
|
|
// a uno).
|
|
double dScProd = - 2 ;
|
|
|
|
if ( bTest0 && bTest1)
|
|
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
|
|
|
|
double dThreshold = 0.7 ;
|
|
|
|
if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
|
|
|
|
int nt = 0 ;
|
|
|
|
while ( nIndexVsIndex3[nt][0] != nIndex)
|
|
++ nt ;
|
|
|
|
int nRotCase = nIndexVsIndex3[nt][1] ;
|
|
|
|
nComponents = Cases3Plus[nRotCase][1][0] ;
|
|
|
|
// Riaggiorno gli offsets
|
|
nExtTabOff = nComponents ;
|
|
nStdTabOff = 0 ;
|
|
|
|
// Modifico le matrici
|
|
for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
|
|
|
|
// Numero vertici per componenti
|
|
nVertComp[nC - 1] = Cases3Plus[nRotCase][1][nC] ;
|
|
|
|
// Matrice dei vertici della base del fan
|
|
for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert)
|
|
|
|
CompoVert[nC - 1][nFanVert] = VecField[Cases3Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ;
|
|
|
|
// Matrici dei vertici dei triangoli in assenza di sharp feature
|
|
for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) {
|
|
|
|
CompoTriVert[nC - 1][nTriVert] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
|
|
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
|
|
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
|
|
}
|
|
|
|
// Aggiorno gli offsets per raggiungere i
|
|
// vertici della componente successiva.
|
|
nExtTabOff += nVertComp[nC - 1] ;
|
|
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
|
|
}
|
|
|
|
|
|
}
|
|
}
|
|
|
|
else if ( nAllConfig[nIndex] == 6) {
|
|
|
|
// Procedura analoga a quella della configurazione 3
|
|
Vector3d vtCmpAvg0, vtCmpAvg1 ;
|
|
|
|
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.7) ;
|
|
bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
|
|
|
|
double dScProd = - 2 ;
|
|
|
|
if ( bTest0 && bTest1)
|
|
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
|
|
|
|
double dThreshold = 0.7 ;
|
|
|
|
if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
|
|
|
|
int nt = 0 ;
|
|
|
|
while ( nIndexVsIndex6[nt][0] != nIndex)
|
|
++ nt ;
|
|
|
|
int nRotCase = nIndexVsIndex6[nt][1] ;
|
|
|
|
// Costruzione dei triangoli
|
|
for ( int TriIndex = 0 ; TriIndex < 15 ; TriIndex += 3) {
|
|
|
|
// Costruzione triangolo
|
|
int i0 = Cases6Plus[nRotCase][TriIndex + 2] ;
|
|
int i1 = Cases6Plus[nRotCase][TriIndex + 1] ;
|
|
int i2 = Cases6Plus[nRotCase][TriIndex] ;
|
|
|
|
Triangle3d CurrentTriangle ;
|
|
|
|
// Il triangolo è pronto
|
|
CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ;
|
|
bool bV = CurrentTriangle.Validate( true) ;
|
|
CurrentTriangle.ToGlob( m_MapFrame) ;
|
|
|
|
// Aggiungo alla lista
|
|
lstTria.emplace_back( CurrentTriangle) ;
|
|
}
|
|
continue ;
|
|
}
|
|
}
|
|
|
|
else if ( nAllConfig[nIndex] == 10) {
|
|
|
|
// Verifico se i versori delle componenti sono tutti
|
|
// più o meno concordi (soglia su prodotto scalare 0, forse 0.7?)
|
|
Vector3d vtCmpAvg0, vtCmpAvg1 ;
|
|
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ;
|
|
bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ;
|
|
|
|
if ( ! bTest0 || ! bTest1) {
|
|
int nt = 0 ;
|
|
while ( nIndexVsIndex10[nt][0] != nIndex)
|
|
++ nt ;
|
|
|
|
// Riaggiorno gli offsets
|
|
nExtTabOff = 2 ;
|
|
nStdTabOff = 0 ;
|
|
|
|
// Modifico le matrici
|
|
int nRotCase = nIndexVsIndex10[nt][1] ;
|
|
for ( int nC = 1 ; nC <= 2 ; ++ nC) {
|
|
// Numero vertici per componenti
|
|
nVertComp[nC - 1] = Cases10Plus[nRotCase][1][nC] ;
|
|
|
|
// Matrice dei vertici della base del fan
|
|
for ( int nFanVert = 0 ; nFanVert < 4 ; ++ nFanVert)
|
|
CompoVert[nC - 1][nFanVert] = VecField[Cases10Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ;
|
|
|
|
// Matrici dei vertici dei triangoli in assenza di sharp feature
|
|
for ( int nTriVert = 0 ; nTriVert < 6 ; nTriVert += 3) {
|
|
CompoTriVert[nC - 1][nTriVert] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
|
|
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
|
|
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
|
|
}
|
|
|
|
// Aggiorno gli offsets per raggiungere i vertici della componente successiva.
|
|
nExtTabOff += nVertComp[nC - 1] ;
|
|
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Numero di feature nel voxel: al più vi è una feature per componente connessa.
|
|
int nInnerFeatureInVoxel = 0 ;
|
|
int nBorderFeatureInVoxel = 0 ;
|
|
|
|
// Ciclo sulle componenti
|
|
for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) {
|
|
|
|
int nFeatureType = NO_FEATURE ;
|
|
// Se i componenti sono regolari valuto le normali per stabilire se eseguire ExtMC o MC
|
|
if ( bReg)
|
|
nFeatureType = TestOnNormal( CompoVert[nCompCount - 1], nVertComp[nCompCount - 1]) ;
|
|
|
|
// Controllo per il caso piano su una griglia
|
|
// con versori normali a due a due paralleli.
|
|
bool bGridControl = true ;
|
|
if ( nFeatureType != NO_FEATURE) {
|
|
if ( nVertComp[nCompCount - 1] == 4)
|
|
GridControl( VecField, CompoVert, nIndArrey, nVertComp, nCompCount, bGridControl) ;
|
|
}
|
|
|
|
// Flag ExtMC
|
|
bool bExtMC = ( nFeatureType != NO_FEATURE && bGridControl) ;
|
|
|
|
// Extended MC
|
|
if ( bExtMC) {
|
|
|
|
// Passo al sistema di riferimento del baricentro
|
|
Point3d ptGravityCenter( 0, 0, 0) ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni)
|
|
ptGravityCenter = ptGravityCenter + CompoVert[nCompCount - 1][ni].ptInt ;
|
|
ptGravityCenter = ptGravityCenter / nVertComp[nCompCount - 1] ;
|
|
|
|
Vector3d vtTrasf[12] ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni)
|
|
vtTrasf[ni] = CompoVert[nCompCount - 1][ni].ptInt - ptGravityCenter ;
|
|
|
|
// Preparo le matrici per il sistema
|
|
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> dSystemMatrix ;
|
|
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> dSystemVector ;
|
|
typedef Eigen::JacobiSVD<dSystemMatrix> DecomposerSVD ;
|
|
|
|
dSystemMatrix dMatrixN( nVertComp[nCompCount - 1], 3) ;
|
|
dSystemVector dKnownVector( nVertComp[nCompCount - 1], 1) ;
|
|
|
|
// Caso generale
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ;
|
|
dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ;
|
|
dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ;
|
|
dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * vtTrasf[ni] ;
|
|
}
|
|
|
|
// calcolo SVD
|
|
DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ;
|
|
dSystemMatrix dMatrixV = svd.matrixV() ;
|
|
dSystemVector dSingularValue = svd.singularValues() ;
|
|
|
|
// Se la feature è un edge annullo il valore singolare minore.
|
|
if ( nFeatureType == EDGE) {
|
|
double dThres = 0.5 * ( dSingularValue( 1) + dSingularValue( 2)) / dSingularValue( 0) ;
|
|
svd.setThreshold( dThres) ;
|
|
}
|
|
|
|
// risolvo il sistema con SVD, quindi ai minimi quadrati
|
|
dSystemVector dUnknownVector( 3, 1) ;
|
|
dUnknownVector = svd.solve( dKnownVector) ;
|
|
|
|
// Vettore Baricentro-Feature
|
|
Vector3d vtFeature( dUnknownVector( 0), dUnknownVector( 1), dUnknownVector( 2)) ;
|
|
|
|
// Esprimo la soluzione nel sistema di riferimento z-map
|
|
Point3d ptSol = ptGravityCenter + vtFeature ;
|
|
|
|
// Limito la feature entro i secondi voxel vicini,
|
|
// nel caso essa esca dal reticolo la limito entro una
|
|
// distanza dal baricentro pari alla diagonale del voxel.
|
|
bool bOutside = false ;
|
|
|
|
int nFtVxI, nFtVxJ, nFtVxK ;
|
|
if ( GetPointVoxel( ptSol, nFtVxI, nFtVxJ, nFtVxK)) {
|
|
|
|
if ( abs( nFtVxI - i) > 2 ||
|
|
abs( nFtVxJ - j) > 2 ||
|
|
abs( nFtVxK - k) > 2)
|
|
bOutside = true ;
|
|
}
|
|
else {
|
|
double dDistFeature = vtFeature.Len() ;
|
|
const double MAX_DIST = sqrt( 3) * N_DEXVOXRATIO * m_dStep ;
|
|
bOutside = ( dDistFeature > MAX_DIST) ;
|
|
}
|
|
|
|
TRIA3DVECTOR triContainer ;
|
|
|
|
// Costruisco triangoli di prova
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
|
|
// Il triangolo è pronto
|
|
Triangle3d CurrentTriangle ;
|
|
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
// Aggiungo triangolo al vettore temporaneo
|
|
triContainer.emplace_back( CurrentTriangle) ;
|
|
}
|
|
|
|
// Controllo delle inversioni dei triangoli
|
|
bool bDangerInversion = false ;
|
|
|
|
// Caso ventaglio con tre vertici di base
|
|
if ( nVertComp[nCompCount - 1] == 3) {
|
|
|
|
// Controllo se esiste almeno un triangolo con normale avente prodotto scalare
|
|
// negativo con la normale di almeno uno dei vertici di base del ventaglio.
|
|
bool bInversione = false ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
|
|
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
|
|
|
|
double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
|
|
double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ;
|
|
|
|
if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) {
|
|
bInversione = true ;
|
|
break ;
|
|
}
|
|
}
|
|
|
|
// Se tale triangolo esiste procedo
|
|
if ( bInversione) {
|
|
|
|
// Conto le coppie di normali con angolo compreso maggiore di 90°
|
|
int nNegDot = 0 ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] - 1 ; ++ ni) {
|
|
for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) {
|
|
if ( CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL)
|
|
nNegDot ++ ;
|
|
}
|
|
}
|
|
|
|
if ( nNegDot == nVertComp[nCompCount - 1] - 1) {
|
|
// Cerco se esiste un punto in cui la normale ha prodotto scalre negativo
|
|
// con le normali di entrambi i triangoli che lo contengono
|
|
bool bInversione2 = false ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
int nj = ( ni == 0 ? nVertComp[nCompCount - 1] - 1 : ni - 1) ;
|
|
double dDLast = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
|
|
double dDPrev = triContainer[nj].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
|
|
if ( ( dDLast < EPS_SMALL && dDPrev < EPS_SMALL) ||
|
|
( dDLast < - 0.75 || dDPrev < - 0.75)) {
|
|
bInversione2 = true ;
|
|
break ;
|
|
}
|
|
}
|
|
|
|
// Se tale normale esiste
|
|
if ( bInversione2) {
|
|
|
|
// Soluzione del sistema nel sistema Zmap
|
|
Point3d ptSolZMapFrame = ptSol ;
|
|
|
|
// Se la soluzione non cade nel voxel di appartenenza vedo se può
|
|
// essere riportata dentro muovendosi lungo la linea di feature.
|
|
if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) {
|
|
|
|
Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
|
|
double dParInt1, dParInt2 ;
|
|
Point3d ptVoxMin( ( i * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( j * N_DEXVOXRATIO + 0.5) * m_dStep, ( k * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
Point3d ptVoxMax( ( ( i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( k + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
// Caso in cui può essere riportata dentro: se il voxel in cui cade la feature è pieno
|
|
// la riporto muovendola lungo la sua linea
|
|
if ( IntersLineBox( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) {
|
|
|
|
triContainer.resize( 0) ;
|
|
|
|
double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 :
|
|
dParInt2 + ( dParInt1 - dParInt2) / 100 ;
|
|
|
|
Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ;
|
|
|
|
ptSol = ptNewSol ;
|
|
|
|
// Costruisco triangoli di prova
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
|
|
// Il triangolo è pronto
|
|
Triangle3d CurrentTriangle ;
|
|
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
// Aggiungo triangolo al vettore temporaneo
|
|
triContainer.emplace_back( CurrentTriangle) ;
|
|
}
|
|
}
|
|
|
|
// Caso in cui non può essere riportata dentro:
|
|
// si passa alla routine standard se il voxel
|
|
// in cui cade è pieno.
|
|
else {
|
|
|
|
int nAdjVoxI, nAdjVoxJ, nAdjVoxK ;
|
|
if ( GetPointVoxel( ptSolZMapFrame, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) {
|
|
|
|
// Classificazione del voxel adiacente
|
|
int nAdjIndex = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ;
|
|
|
|
// Se il voxel è pieno
|
|
if ( EdgeTable[nAdjIndex] != 0)
|
|
bDangerInversion = true ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Ventaglio con base a quattro vertici
|
|
else if ( nVertComp[nCompCount - 1] == 4) {
|
|
|
|
// Controllo se esiste almeno un triangolo con normale avente prodotto scalare
|
|
// negativo con la normale di almeno uno dei vertici di base del ventaglio.
|
|
bool bInversione = false ;
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
|
|
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
|
|
|
|
double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
|
|
double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ;
|
|
|
|
if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL)
|
|
bInversione = true ;
|
|
}
|
|
// Se tale triangolo esiste continuo i test
|
|
if ( bInversione) {
|
|
|
|
// Se la feature non cade nel suo voxel
|
|
if ( ! IsPointInsideVoxelApprox( i, j, k, ptSol)) {
|
|
|
|
Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
|
|
|
|
double dParInt1, dParInt2 ;
|
|
Point3d ptVoxMin( ( i * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( j * N_DEXVOXRATIO + 0.5) * m_dStep, ( k * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
Point3d ptVoxMax( ( ( i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( k + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
// Se è possibile riportarla dentro e il voxel in cui cade è pieno, la riporto nel suo voxel
|
|
// lungo la sua linea
|
|
if ( IntersLineBox( ptSol, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) {
|
|
|
|
triContainer.resize( 0) ;
|
|
|
|
double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 :
|
|
dParInt2 + ( dParInt1 - dParInt2) / 100 ;
|
|
|
|
Point3d ptNewSol = ptSol + dPar * vtNullSpace ;
|
|
|
|
ptSol = ptNewSol ;
|
|
|
|
// Costruisco triangoli di prova
|
|
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
|
|
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
|
|
// Il triangolo è pronto
|
|
Triangle3d CurrentTriangle ;
|
|
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
// Aggiungo triangolo al vettore temporaneo
|
|
triContainer.emplace_back( CurrentTriangle) ;
|
|
}
|
|
}
|
|
|
|
// Se non è possibile riportarla dentro e il voxel in
|
|
// cui cade è pieno passo alla routine standard
|
|
else {
|
|
int nAdjVoxI, nAdjVoxJ, nAdjVoxK ;
|
|
if ( GetPointVoxel( ptSol, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) {
|
|
// Classificazione del voxel adiacente
|
|
int nAdjIndex = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ;
|
|
// Se il voxel è pieno
|
|
if ( EdgeTable[nAdjIndex] != 0)
|
|
bDangerInversion = true ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Controllo sulle normali, se risultano in tutti i punti
|
|
// approssimativamente uguali passiamo alla routine standard
|
|
int nContSize = int( triContainer.size()) ;
|
|
bool bPlane = true ;
|
|
for ( int ni = 0 ; ni < nContSize - 1 ; ++ ni) {
|
|
for ( int nj = ni + 1 ; nj < nContSize ; ++ nj) {
|
|
Vector3d vtI = triContainer[ni].GetN() ;
|
|
Vector3d vtJ = triContainer[nj].GetN() ;
|
|
if ( ! AreSameVectorApprox( vtI, vtJ)) {
|
|
bPlane = false ;
|
|
break ;
|
|
}
|
|
}
|
|
if ( ! bPlane)
|
|
break ;
|
|
}
|
|
|
|
// Se feature nei limiti e triangoli non in un piano, confermo ExtMC
|
|
if ( ! bOutside && ! bPlane && ! bDangerInversion) {
|
|
|
|
// Riconoscimento se voxel di frontiera
|
|
int nVoxIndexes[3] = { i, j, k} ;
|
|
bool bBoundary = IsAVoxelOnBoundary( nLimits, nVoxIndexes, true) ;
|
|
|
|
TRIA3DVECTOR vInnerTriaTemp, vBorderTriaTemp ;
|
|
|
|
for ( int ni = 0 ; ni < nContSize ; ++ ni) {
|
|
|
|
// Se l'area è non nulla aggiungo il triangolo a uno dei due vettori.
|
|
if ( triContainer[ni].GetArea() > SQ_EPS_SMALL) {
|
|
int nVoxIJK[3] = { i, j, k} ;
|
|
Point3d ptTemp = ptSol ;
|
|
Triangle3d trTemp = triContainer[ni] ;
|
|
bool bTriangleOnBorder = IsATriangleOnBorder( trTemp, ptTemp, nLimits, nVoxIJK) ;
|
|
if ( ! bBoundary || ! bTriangleOnBorder)
|
|
vInnerTriaTemp.emplace_back( triContainer[ni]) ;
|
|
else
|
|
vBorderTriaTemp.emplace_back( triContainer[ni]) ;
|
|
}
|
|
}
|
|
|
|
// Porto la soluzione nel sistema in cui è immerso lo Zmap
|
|
ptSol.ToGlob( m_MapFrame) ;
|
|
|
|
// Metto i triangoli negli appositi contenitori:
|
|
|
|
// Triangoli interni
|
|
if ( vInnerTriaTemp.size() > 0) {
|
|
|
|
// Aggiorno il numero di feature.
|
|
++ nInnerFeatureInVoxel ;
|
|
|
|
// Se siamo in presenza della prima feature del
|
|
// voxel, ridimensiono il vettore che contiene
|
|
// la struttura dati del voxel.
|
|
if ( nInnerFeatureInVoxel == 1) {
|
|
triHold.resize( triHold.size() + 1) ;
|
|
// Questi dati dipendono solo dal voxel,
|
|
int nCurrent = int( triHold.size()) - 1 ;
|
|
triHold[nCurrent].i = i ;
|
|
triHold[nCurrent].j = j ;
|
|
triHold[nCurrent].k = k ;
|
|
}
|
|
|
|
// Indice che identifica la posizione del voxel nel vector.
|
|
int nCurrent = int( triHold.size()) - 1 ;
|
|
|
|
// Aggiungo vertice della componente
|
|
// connessa alla lista dei vertici.
|
|
triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ;
|
|
|
|
int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ;
|
|
int nNewFeatureNum = nOldFeatureNum + 1 ;
|
|
|
|
// Aggiungo una componente per il vettore
|
|
// dei triangoli della componente connessa.
|
|
triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ;
|
|
for ( int ni = 0 ; ni < int( vInnerTriaTemp.size()) ; ++ ni) {
|
|
// Riporto le coordinate nel sistema in cui è immerso lo Zmap
|
|
vInnerTriaTemp[ni].ToGlob( m_MapFrame) ;
|
|
triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vInnerTriaTemp[ni]) ;
|
|
}
|
|
}
|
|
|
|
// Triangoli di frontiera
|
|
if ( vBorderTriaTemp.size() > 0) {
|
|
|
|
// Aggiorno il numero di feature.
|
|
++ nBorderFeatureInVoxel ;
|
|
|
|
// Se siamo in presenza della prima feature del
|
|
// voxel, ridimensiono il vettore che contiene
|
|
// la struttura dati del voxel.
|
|
if ( nBorderFeatureInVoxel == 1) {
|
|
m_InterBlockTria[nBlock].resize( m_InterBlockTria[nBlock].size() + 1) ;
|
|
// Questi dati dipendono solo dal voxel
|
|
int nCurrent = int( m_InterBlockTria[nBlock].size()) - 1 ;
|
|
m_InterBlockTria[nBlock][nCurrent].i = i ;
|
|
m_InterBlockTria[nBlock][nCurrent].j = j ;
|
|
m_InterBlockTria[nBlock][nCurrent].k = k ;
|
|
}
|
|
|
|
// Indice che identifica la posizione del voxel nel vector
|
|
int nCurrent = int( m_InterBlockTria[nBlock].size()) - 1 ;
|
|
|
|
// Aggiungo vertice della componente connessa alla lista dei vertici
|
|
m_InterBlockTria[nBlock][nCurrent].ptCompoVert.emplace_back( ptSol) ;
|
|
|
|
int nOldFeatureNum = int( m_InterBlockTria[nBlock][nCurrent].vCompoTria.size()) ;
|
|
int nNewFeatureNum = nOldFeatureNum + 1 ;
|
|
|
|
// Aggiungo una componente per il vettore dei triangoli della componente connessa
|
|
m_InterBlockTria[nBlock][nCurrent].vCompoTria.resize( nNewFeatureNum) ;
|
|
for ( int ni = 0 ; ni < int( vBorderTriaTemp.size()) ; ++ ni) {
|
|
// Riporto le coordinate nel sistema in cui è immerso lo Zmap
|
|
vBorderTriaTemp[ni].ToGlob( m_MapFrame) ;
|
|
m_InterBlockTria[nBlock][nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vBorderTriaTemp[ni]) ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// ExtMC non confermato, si passa a MC
|
|
else {
|
|
// Costruzione dei triangoli
|
|
for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {
|
|
// Il triangolo è pronto
|
|
Triangle3d CurrentTriangle ;
|
|
CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt,
|
|
CompoTriVert[nCompCount - 1][TriIndex+1].ptInt,
|
|
CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ;
|
|
bool bV = CurrentTriangle.Validate( true) ;
|
|
// Riporto le coordinate nel sistema in cui è immerso lo Zmap
|
|
CurrentTriangle.ToGlob( m_MapFrame) ;
|
|
// Aggiungo alla lista
|
|
lstTria.emplace_back( CurrentTriangle) ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Standard MC
|
|
else {
|
|
// Costruzione dei triangoli
|
|
for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {
|
|
// Il triangolo è pronto
|
|
Triangle3d CurrentTriangle ;
|
|
CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt,
|
|
CompoTriVert[nCompCount - 1][TriIndex+1].ptInt,
|
|
CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ;
|
|
bool bV = CurrentTriangle.Validate( true) ;
|
|
// Riporto le coordinate nel sistema in cui è immerso lo Zmap
|
|
CurrentTriangle.ToGlob( m_MapFrame) ;
|
|
// Aggiungo alla lista
|
|
lstTria.emplace_back( CurrentTriangle) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Processo i Voxel con possibile superficie piana
|
|
ProcessVoxContXY( VoxContXYInf, false, lstTria) ;
|
|
ProcessVoxContXY( VoxContXYSup, true, lstTria) ;
|
|
ProcessVoxContYZ( VoxContYZInf, false, lstTria) ;
|
|
ProcessVoxContYZ( VoxContYZSup, true, lstTria) ;
|
|
ProcessVoxContXZ( VoxContXZInf, false, lstTria) ;
|
|
ProcessVoxContXZ( VoxContXZSup, true, lstTria) ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Calcola i triangoli dei voxel contenuti nel vettore, se si vuole il riconoscimento
|
|
// delle sharp-feature passare true in bEhn. In questo caso vengono analizzai anche i
|
|
// voxel adiacenti a quelli contenenti sharp-feature, e viene eseguito il flipping.
|
|
bool
|
|
VolZmap::ExtMarchingCubes( vector<Voxel>& vVox, TRIA3DLIST& lstTria, bool bEnh) const
|
|
{
|
|
// Contenitore dei triangoli costituenti sharp-feaure
|
|
TriHolder triHold ;
|
|
int nOriginalSize = int( vVox.size()) ;
|
|
for ( int nV = 0 ; nV < int( vVox.size()) ; ++ nV) {
|
|
// I voxel intersecati dalla retta hanno indice nel vettore
|
|
// nV < OriginalSize; quelli adiacenti avranno indice
|
|
// nV >= OriginalSize
|
|
if ( nV < nOriginalSize) {
|
|
// Se a questa iterazione troviamo una sharp-feature
|
|
// la dimensione del TriHolder aumenta
|
|
int nCurrentSize = int( triHold.size()) ;
|
|
ProcessCube( vVox[nV].nI, vVox[nV].nJ, vVox[nV].nK, lstTria, triHold, bEnh) ;
|
|
int nNewSize = int( triHold.size()) ;
|
|
|
|
// Se abbiamo trovato una sharp-feature
|
|
if ( nNewSize > nCurrentSize) {
|
|
// Ciclo fra tutti i possibili voxel adiacenti
|
|
for ( int i = - 1 ; i <= 1 ; ++ i) {
|
|
for ( int j = - 1 ; j <= 1 ; ++ j) {
|
|
for ( int k = - 1 ; k <= 1 ; ++ k) {
|
|
if ( i == 0 && j == 0 && k == 0)
|
|
continue ;
|
|
if ( IsValidVoxel( vVox[nV].nI + i, vVox[nV].nJ + j, vVox[nV].nK + k)) {
|
|
// Se il voxel adiacente è già nella lista di quelli attraversati dalla retta
|
|
// lo ignoro, altrimenti lo metto in coda nel TriHolder.
|
|
int nOldV = 0 ;
|
|
for ( ; nOldV < nOriginalSize ; ++ nOldV) {
|
|
|
|
if ( vVox[nOldV].nI == vVox[nV].nI + i ||
|
|
vVox[nOldV].nJ == vVox[nV].nJ + j ||
|
|
vVox[nOldV].nK == vVox[nV].nK + k)
|
|
continue ;
|
|
}
|
|
if ( nOldV == nOriginalSize) {
|
|
Voxel NewVox ;
|
|
NewVox.nI = vVox[nV].nI + i ;
|
|
NewVox.nJ = vVox[nV].nJ + j ;
|
|
NewVox.nK = vVox[nV].nK + k ;
|
|
vVox.emplace_back( NewVox) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
// Se il voxel non è attraversato dalla retta, ma fa parte di
|
|
// quelli adiacenti lo processiamo normalmente
|
|
else
|
|
ProcessCube( vVox[nV].nI, vVox[nV].nJ, vVox[nV].nK, lstTria, triHold, bEnh) ;
|
|
}
|
|
// Eseguo il flipping
|
|
FlipEdgesII( triHold, false) ;
|
|
// Aggiungo alla lista di triangoli, tutti quelli che formano una sharp-feature.
|
|
for ( int nVox = 0 ; nVox < int( triHold.size()) ; ++ nVox) {
|
|
for ( int nCompo = 0 ; nCompo < int( triHold[nVox].vCompoTria.size()) ; ++ nCompo) {
|
|
for ( int nTri = 0 ; nTri < int( triHold[nVox].vCompoTria[nCompo].size()) ; ++ nTri) {
|
|
lstTria.emplace_back( triHold[nVox].vCompoTria[nCompo][nTri]) ;
|
|
}
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Esegue il flipping dei triangoli contenuti nel TriHolder.
|
|
// Se si usa la funzione per la visualizzazione dei triangoli,
|
|
// bGraph deve valere true, se si intende calcolare la profondità
|
|
// del materiale bGraph deve valere false.
|
|
bool
|
|
VolZmap::FlipEdgesII( TriHolder& TriHold, bool bGraph) const
|
|
{
|
|
// Numero di voxel in cui si presentano sharp feature
|
|
int nVoxelNum = int( TriHold.size()) ;
|
|
|
|
// Ciclo su tali voxel
|
|
for ( int n1 = 0 ; n1 < nVoxelNum ; ++ n1) {
|
|
|
|
for ( int n2 = n1 ; n2 < nVoxelNum ; ++ n2) {
|
|
|
|
// Flag voxel adiacenti o meno
|
|
bool bAdj ;
|
|
// Se la funzione è chiamata per la grafica i voxel sono ordinati nel TriHolder,
|
|
// non serve confrontare i triangoli con quelli dei voxel precedenti.
|
|
if ( bGraph) {
|
|
bAdj = ( TriHold[n2].i >= TriHold[n1].i && TriHold[n2].i <= TriHold[n1].i + 1) ||
|
|
( TriHold[n2].j >= TriHold[n1].j && TriHold[n2].j <= TriHold[n1].j + 1) ||
|
|
( TriHold[n2].k >= TriHold[n1].k && TriHold[n2].k <= TriHold[n1].k + 1) ;
|
|
}
|
|
// Se la funzione è chiamata per il calcolo della profondità del materiale i voxel non
|
|
// sono ordinati nel TriHolder, quindi può capitare di dover confrontare i triangoli con
|
|
// quelli dei voxel precedenti.
|
|
else {
|
|
bAdj = abs( TriHold[n2].i - TriHold[n1].i) == 1 ||
|
|
abs( TriHold[n2].j - TriHold[n1].j) == 1 ||
|
|
abs( TriHold[n2].k - TriHold[n1].k) == 1 ;
|
|
}
|
|
|
|
// Se i voxel sono adiacenti proseguo
|
|
if ( ( TriHold[n2].i >= TriHold[n1].i && TriHold[n2].i <= TriHold[n1].i + 1) ||
|
|
( TriHold[n2].j >= TriHold[n1].j && TriHold[n2].j <= TriHold[n1].j + 1) ||
|
|
( TriHold[n2].k >= TriHold[n1].k && TriHold[n2].k <= TriHold[n1].k + 1)) {
|
|
|
|
// Numero delle componenti connesse nei due voxel
|
|
int nNumCompo1 = int( TriHold[n1].ptCompoVert.size()) ;
|
|
int nNumCompo2 = int( TriHold[n2].ptCompoVert.size()) ;
|
|
|
|
int nCompo1 = 0 ;
|
|
|
|
// Ciclo sulle componenti
|
|
for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) {
|
|
|
|
int nCompo2 = ( n1 == n2 ? nCompo1 + 1 : 0) ;
|
|
|
|
for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) {
|
|
|
|
// Numero di triangoli per le componenti connesse
|
|
int nTriNum1 = int( TriHold[n1].vCompoTria[nCompo1].size()) ;
|
|
int nTriNum2 = int( TriHold[n2].vCompoTria[nCompo2].size()) ;
|
|
|
|
for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) {
|
|
|
|
bool bModified = false ;
|
|
|
|
for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) {
|
|
|
|
INTVECTOR SharedIndex ;
|
|
|
|
for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) {
|
|
|
|
for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) {
|
|
|
|
Point3d ptP1 = TriHold[n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ;
|
|
Point3d ptP2 = TriHold[n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ;
|
|
|
|
if ( AreSamePointEpsilon( ptP1, ptP2, EPS_ZERO)) {
|
|
|
|
Point3d ptVert1 = TriHold[n1].ptCompoVert[nCompo1] ;
|
|
Point3d ptVert2 = TriHold[n2].ptCompoVert[nCompo2] ;
|
|
|
|
if ( ! ( AreSamePointEpsilon( ptP1, ptVert1, EPS_ZERO) ||
|
|
AreSamePointEpsilon( ptP2, ptVert2, EPS_ZERO))) {
|
|
|
|
SharedIndex.emplace_back( nVert1) ;
|
|
SharedIndex.emplace_back( nVert2) ;
|
|
}
|
|
}
|
|
|
|
if ( SharedIndex.size() > 2)
|
|
break ;
|
|
}
|
|
|
|
if ( SharedIndex.size() > 2)
|
|
break ;
|
|
}
|
|
|
|
// Si deve operare la modifica dei triangoli
|
|
if ( SharedIndex.size() > 2) {
|
|
|
|
// Verifico che i due lati adiacenti siano controversi
|
|
int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ;
|
|
if ( nProd != 1 && nProd != - 2 && nProd != 4 && nProd != 0) {
|
|
// Assegno l'array dei punti di contorno
|
|
Point3d vPnt[4] ;
|
|
vPnt[0] = TriHold[n1].vCompoTria[nCompo1][nTri1].GetP( SharedIndex[0]) ;
|
|
vPnt[1] = TriHold[n1].ptCompoVert[nCompo1] ;
|
|
vPnt[2] = TriHold[n1].vCompoTria[nCompo1][nTri1].GetP( SharedIndex[2]) ;
|
|
vPnt[3] = TriHold[n2].ptCompoVert[nCompo2] ;
|
|
// Valuto se i triangoli giacciono su un piano
|
|
PolygonPlane Polygon ;
|
|
for ( int i = 0 ; i < 4 ; ++ i)
|
|
Polygon.AddPoint( vPnt[i]) ;
|
|
Plane3d plPlane ;
|
|
bool bOnPlane = Polygon.GetPlane( plPlane) ;
|
|
for ( int i = 0 ; i < 4 && bOnPlane ; ++ i)
|
|
bOnPlane = PointInPlaneApprox( vPnt[i], plPlane) ;
|
|
// Se sono su un piano controllo se avviene inversione
|
|
bool bInv = false ;
|
|
if ( bOnPlane) {
|
|
Triangle3d trT1 = TriHold[n1].vCompoTria[nCompo1][nTri1] ;
|
|
Triangle3d trT2 = TriHold[n2].vCompoTria[nCompo2][nTri2] ;
|
|
int nVert1, nVert2 ;
|
|
// Determino gli indici dei punti sharp-feature
|
|
for ( int nP = 0 ; nP < 3 ; ++ nP) {
|
|
if ( nP != SharedIndex[0] && nP != SharedIndex[2])
|
|
nVert1 = nP ;
|
|
if ( nP != SharedIndex[1] && nP != SharedIndex[3])
|
|
nVert2 = nP ;
|
|
}
|
|
trT1.SetP( SharedIndex[0], trT2.GetP( nVert2)) ;
|
|
trT2.SetP( SharedIndex[3], trT1.GetP( nVert1)) ;
|
|
trT1.Validate( true) ;
|
|
trT2.Validate( true) ;
|
|
bInv = ( trT1.GetN() * trT2.GetN() < 0) ;
|
|
}
|
|
// Se non vi è inversione eseguo il flipping
|
|
if ( ! bInv) {
|
|
TriHold[n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0],
|
|
TriHold[n2].ptCompoVert[nCompo2]) ;
|
|
TriHold[n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3],
|
|
TriHold[n1].ptCompoVert[nCompo1]) ;
|
|
|
|
TriHold[n1].vCompoTria[nCompo1][nTri1].Validate( true) ;
|
|
TriHold[n2].vCompoTria[nCompo2][nTri2].Validate( true) ;
|
|
|
|
bModified = true ;
|
|
break ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if ( bModified)
|
|
break ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const
|
|
{
|
|
// Numero di blocchi
|
|
size_t nBlocksNum = InterTria.size() ;
|
|
|
|
// ciclo sui blocchi
|
|
for ( size_t tFB = 0 ; tFB < nBlocksNum ; ++ tFB) {
|
|
|
|
int nFBijk[3] ;
|
|
GetBlockIJKFromN( int( tFB), nFBijk) ;
|
|
|
|
for ( size_t tLB = tFB ; tLB < nBlocksNum ; ++ tLB) {
|
|
|
|
int nLBijk[3] ;
|
|
GetBlockIJKFromN( int( tLB), nLBijk) ;
|
|
|
|
// Se i blocchi non sono adiacenti salto l'iterazione
|
|
if ( ! ( abs( nFBijk[0] - nLBijk[0]) <= 1 &&
|
|
abs( nFBijk[1] - nLBijk[1]) <= 1 &&
|
|
abs( nFBijk[2] - nLBijk[2]) <= 1))
|
|
continue ;
|
|
|
|
// Numero di voxel nei blocchi correnti
|
|
size_t nVoxelNumFB = InterTria[tFB].size() ;
|
|
size_t nVoxelNumLB = InterTria[tLB].size() ;
|
|
|
|
// Ciclo sui voxel dei due blocchi
|
|
for ( size_t tVFB = 0 ; tVFB < nVoxelNumFB ; ++ tVFB) {
|
|
for ( size_t tVLB = 0 ; tVLB < nVoxelNumLB ; ++ tVLB) {
|
|
|
|
// Se i voxel non sono adiacenti salto l'iterazione
|
|
if ( ! ( abs( InterTria[tFB][tVFB].i - InterTria[tLB][tVLB].i) <= 1 &&
|
|
abs( InterTria[tFB][tVFB].j - InterTria[tLB][tVLB].j) <= 1 &&
|
|
abs( InterTria[tFB][tVFB].k - InterTria[tLB][tVLB].k) <= 1))
|
|
|
|
continue ;
|
|
|
|
// Numero di componenti connesse dei voxel
|
|
size_t nCompoVFBNum = InterTria[tFB][tVFB].ptCompoVert.size() ;
|
|
size_t nCompoVLBNum = InterTria[tLB][tVLB].ptCompoVert.size() ;
|
|
|
|
// Ciclo sulle componenti connesse
|
|
for ( size_t tCmpF = 0 ; tCmpF < nCompoVFBNum ; ++ tCmpF) {
|
|
for ( size_t tCmpL = 0 ; tCmpL < nCompoVLBNum ; ++ tCmpL) {
|
|
|
|
// Numero di triangoli delle componenti connesse
|
|
size_t nTriFBNum = InterTria[tFB][tVFB].vCompoTria[tCmpF].size() ;
|
|
size_t nTriLBNum = InterTria[tLB][tVLB].vCompoTria[tCmpL].size() ;
|
|
|
|
for ( size_t tTriFB = 0 ; tTriFB < nTriFBNum ; ++ tTriFB) {
|
|
|
|
bool bModified = false ;
|
|
|
|
for ( size_t tTriLB = 0 ; tTriLB < nTriLBNum ; ++ tTriLB) {
|
|
|
|
INTVECTOR SharedIndex ;
|
|
|
|
for ( int nVertF = 0 ; nVertF < 3 ; ++ nVertF) {
|
|
for ( int nVertL = 0 ; nVertL < 3 ; ++ nVertL) {
|
|
|
|
Point3d ptPF = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( nVertF) ;
|
|
Point3d ptPL = InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( nVertL) ;
|
|
|
|
if ( AreSamePointEpsilon( ptPF, ptPL, EPS_ZERO)) {
|
|
|
|
Point3d ptVertF = InterTria[tFB][tVFB].ptCompoVert[tCmpF] ;
|
|
Point3d ptVertL = InterTria[tLB][tVLB].ptCompoVert[tCmpL] ;
|
|
|
|
if ( ! ( AreSamePointEpsilon( ptPF, ptVertF, EPS_ZERO) ||
|
|
AreSamePointEpsilon( ptPL, ptVertL, EPS_ZERO))) {
|
|
|
|
SharedIndex.emplace_back( nVertF) ;
|
|
SharedIndex.emplace_back( nVertL) ;
|
|
}
|
|
}
|
|
|
|
if ( SharedIndex.size() > 2)
|
|
break ;
|
|
}
|
|
|
|
if ( SharedIndex.size() > 2)
|
|
break ;
|
|
}
|
|
|
|
// Si deve operare la modifica dei triangoli
|
|
if ( SharedIndex.size() > 2) {
|
|
// Verifico che i due lati adiacenti siano controversi
|
|
int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ;
|
|
if ( nProd != 1 && nProd != - 2 && nProd != 4 && nProd != 0) {
|
|
// Assegno l'array dei punti di contorno
|
|
Point3d vPnt[4] ;
|
|
vPnt[0] = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( SharedIndex[0]) ;
|
|
vPnt[1] = InterTria[tFB][tVFB].ptCompoVert[tCmpF] ;
|
|
vPnt[2] = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( SharedIndex[2]) ;
|
|
vPnt[3] = InterTria[tLB][tVLB].ptCompoVert[tCmpL] ;
|
|
// Valuto se i triangoli giacciono su un piano
|
|
PolygonPlane Polygon ;
|
|
for ( int i = 0 ; i < 4 ; ++ i)
|
|
Polygon.AddPoint( vPnt[i]) ;
|
|
Plane3d plPlane ;
|
|
bool bOnPlane = Polygon.GetPlane( plPlane) ;
|
|
for ( int i = 0 ; i < 4 && bOnPlane ; ++ i)
|
|
bOnPlane = PointInPlaneApprox( vPnt[i], plPlane) ;
|
|
// Se sono su un piano controllo se avviene inversione
|
|
bool bInv = false ;
|
|
if ( bOnPlane) {
|
|
Triangle3d trTF = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB] ;
|
|
Triangle3d trTL = InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB] ;
|
|
int nVertF, nVertL ;
|
|
// Determino gli indici dei punti sharp-feature
|
|
for ( int nP = 0 ; nP < 3 ; ++ nP) {
|
|
if ( nP != SharedIndex[0] && nP != SharedIndex[2])
|
|
nVertF = nP ;
|
|
if ( nP != SharedIndex[1] && nP != SharedIndex[3])
|
|
nVertL = nP ;
|
|
}
|
|
trTF.SetP( SharedIndex[0], trTL.GetP( nVertL)) ;
|
|
trTF.Validate( true) ;
|
|
trTL.SetP( SharedIndex[3], trTF.GetP( nVertF)) ;
|
|
trTL.Validate( true) ;
|
|
bInv = ( trTF.GetN() * trTL.GetN() < 0) ;
|
|
}
|
|
// Se non vi è inversione eseguo il flipping
|
|
if ( ! bInv) {
|
|
InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetP( SharedIndex[0],
|
|
InterTria[tLB][tVLB].ptCompoVert[tCmpL]) ;
|
|
InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetP( SharedIndex[3],
|
|
InterTria[tFB][tVFB].ptCompoVert[tCmpF]) ;
|
|
|
|
InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].Validate( true) ;
|
|
InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].Validate( true) ;
|
|
|
|
bModified = true ;
|
|
break ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if ( bModified)
|
|
break ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IsThereMat( int nI, int nJ, int nK) const
|
|
{
|
|
// Trasformo gl indici della griglia voxel in quelli della griglia dexel
|
|
nI *= N_DEXVOXRATIO ;
|
|
nJ *= N_DEXVOXRATIO ;
|
|
nK *= N_DEXVOXRATIO ;
|
|
|
|
// Se l'indice è alla frontiera del reticolo non vi è materiale
|
|
if ( nI <= - 1 || nI >= int( m_nNx[0]) ||
|
|
nJ <= - 1 || nJ >= int( m_nNy[0]) ||
|
|
nK <= - 1 || nK >= int( m_nNy[1]))
|
|
return false ;
|
|
|
|
// ciclo sulle griglie
|
|
int nCount = 0 ;
|
|
for ( int nGrid = 0 ; nGrid < int ( m_nMapNum) ; ++ nGrid) {
|
|
// assegnazione dati vertice dipendenti dalla griglia
|
|
unsigned int nGrI, nGrJ ;
|
|
double dZ ;
|
|
switch ( nGrid) {
|
|
case 0 :
|
|
nGrI = nI ;
|
|
nGrJ = nJ ;
|
|
dZ = ( nK + 0.5) * m_dStep ;
|
|
break ;
|
|
case 1 :
|
|
nGrI = nJ ;
|
|
nGrJ = nK ;
|
|
dZ = ( nI + 0.5) * m_dStep ;
|
|
break ;
|
|
case 2 :
|
|
nGrI = nK ;
|
|
nGrJ = nI ;
|
|
dZ = ( nJ + 0.5) * m_dStep ;
|
|
break ;
|
|
}
|
|
// verifica spillone su vertice
|
|
size_t nIndex = 0 ;
|
|
unsigned int nPos = nGrJ * m_nNx[nGrid] + nGrI ;
|
|
size_t nDexSize = m_Values[nGrid][nPos].size() ;
|
|
while ( nIndex < nDexSize) {
|
|
if ( dZ > m_Values[nGrid][nPos][nIndex].dMin - 2 * EPS_SMALL &&
|
|
dZ < m_Values[nGrid][nPos][nIndex].dMax + 2 * EPS_SMALL) {
|
|
++ nCount ;
|
|
break ;
|
|
}
|
|
nIndex += 1 ;
|
|
}
|
|
}
|
|
return ( nCount == 3) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int
|
|
VolZmap::CalcIndex( int nI, int nJ, int nK) const
|
|
{
|
|
int nIndex = 0 ;
|
|
if ( IsThereMat( nI, nJ, nK))
|
|
nIndex |= ( 1 << 0) ;
|
|
if ( IsThereMat( nI + 1, nJ, nK))
|
|
nIndex |= ( 1 << 1) ;
|
|
if ( IsThereMat( nI + 1, nJ + 1, nK))
|
|
nIndex |= ( 1 << 2) ;
|
|
if ( IsThereMat( nI, nJ + 1, nK))
|
|
nIndex |= ( 1 << 3) ;
|
|
if ( IsThereMat( nI, nJ, nK + 1))
|
|
nIndex |= ( 1 << 4) ;
|
|
if ( IsThereMat( nI + 1, nJ, nK + 1))
|
|
nIndex |= ( 1 << 5) ;
|
|
if ( IsThereMat( nI + 1, nJ + 1, nK + 1))
|
|
nIndex |= ( 1 << 6) ;
|
|
if ( IsThereMat( nI, nJ + 1, nK + 1))
|
|
nIndex |= ( 1 << 7) ;
|
|
return nIndex ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, Vector3d& vtNormal) const
|
|
{
|
|
const double dEps = 2 * EPS_SMALL ;
|
|
|
|
bool bFound = false ;
|
|
|
|
if ( nVec1[0] != nVec2[0]) {
|
|
|
|
int nMinI = min( nVec1[0], nVec2[0]) * N_DEXVOXRATIO ;
|
|
int nMaxI = max( nVec1[0], nVec2[0]) * N_DEXVOXRATIO ;
|
|
|
|
double dMinX = ( nMinI + 0.5) * m_dStep ;
|
|
double dMaxX = ( nMaxI + 0.5) * m_dStep ;
|
|
|
|
ptInt.y = ( nVec1[1] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
ptInt.z = ( nVec1[2] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
|
|
size_t nDexel = ( nVec1[2] * m_nNx[1] + nVec1[1]) * N_DEXVOXRATIO ;
|
|
int nSize = int( m_Values[1][nDexel].size()) ;
|
|
|
|
if ( bFirstCorner) {
|
|
int n = nSize - 1 ;
|
|
double dX = m_Values[1][nDexel][n].dMax ;
|
|
while ( n >= 0 && dX > dMinX - dEps) {
|
|
if ( dX < dMaxX + dEps) {
|
|
ptInt.x = dX ;
|
|
vtNormal = m_Values[1][nDexel][n].vtMaxN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
n -= 1 ;
|
|
if ( n >= 0)
|
|
dX = m_Values[1][nDexel][n].dMax ;
|
|
}
|
|
}
|
|
|
|
else {
|
|
int n = 0 ;
|
|
double dX = m_Values[1][nDexel][n].dMin ;
|
|
while ( n < nSize && dX < dMaxX + dEps) {
|
|
if ( dX > dMinX - dEps) {
|
|
ptInt.x = dX ;
|
|
vtNormal = m_Values[1][nDexel][n].vtMinN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
n += 1 ;
|
|
if ( n < nSize)
|
|
dX = m_Values[1][nDexel][n].dMin ;
|
|
}
|
|
}
|
|
|
|
if ( ! bFound)
|
|
ptInt.x = 0.5 * ( dMinX + dMaxX) ;
|
|
}
|
|
|
|
else if ( nVec1[1] != nVec2[1]) {
|
|
|
|
int nMinJ = min( nVec1[1], nVec2[1]) * N_DEXVOXRATIO ;
|
|
int nMaxJ = max( nVec1[1], nVec2[1]) * N_DEXVOXRATIO ;
|
|
|
|
double dMinY = ( nMinJ + 0.5) * m_dStep ;
|
|
double dMaxY = ( nMaxJ + 0.5) * m_dStep ;
|
|
|
|
ptInt.x = ( nVec1[0] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
ptInt.z = ( nVec1[2] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
|
|
size_t nDexel = ( nVec1[0] * m_nNx[2] + nVec1[2]) * N_DEXVOXRATIO ;
|
|
int nSize = int( m_Values[2][nDexel].size()) ;
|
|
|
|
if ( bFirstCorner) {
|
|
int n = nSize - 1 ;
|
|
double dY = m_Values[2][nDexel][n].dMax ;
|
|
while ( n >= 0 && dY > dMinY - dEps) {
|
|
if ( dY < dMaxY + dEps) {
|
|
ptInt.y = dY ;
|
|
vtNormal = m_Values[2][nDexel][n].vtMaxN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
n -= 1 ;
|
|
if ( n >= 0)
|
|
dY = m_Values[2][nDexel][n].dMax ;
|
|
}
|
|
}
|
|
|
|
else {
|
|
int n = 0 ;
|
|
double dY = m_Values[2][nDexel][n].dMin ;
|
|
while ( n < nSize && dY < dMaxY + dEps) {
|
|
if ( dY > dMinY - dEps) {
|
|
ptInt.y = dY ;
|
|
vtNormal = m_Values[2][nDexel][n].vtMinN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
n += 1 ;
|
|
if ( n < nSize)
|
|
dY = m_Values[2][nDexel][n].dMin ;
|
|
}
|
|
}
|
|
|
|
if ( ! bFound)
|
|
ptInt.y = 0.5 * ( dMinY + dMaxY) ;
|
|
}
|
|
|
|
else if ( nVec1[2] != nVec2[2]) {
|
|
|
|
int nMinK = min( nVec1[2], nVec2[2]) * N_DEXVOXRATIO ;
|
|
int nMaxK = max( nVec1[2], nVec2[2]) * N_DEXVOXRATIO ;
|
|
|
|
double dMinZ = ( nMinK + 0.5) * m_dStep ;
|
|
double dMaxZ = ( nMaxK + 0.5) * m_dStep ;
|
|
|
|
ptInt.x = ( nVec1[0] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
ptInt.y = ( nVec1[1] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
|
|
size_t nDexel = ( nVec1[1] * m_nNx[0] + nVec1[0]) * N_DEXVOXRATIO ;
|
|
int nSize = int( m_Values[0][nDexel].size()) ;
|
|
|
|
if ( bFirstCorner) {
|
|
int n = nSize - 1 ;
|
|
double dZ = m_Values[0][nDexel][n].dMax ;
|
|
while ( n >= 0 && dZ > dMinZ - dEps) {
|
|
if ( dZ < dMaxZ + dEps) {
|
|
ptInt.z = dZ ;
|
|
vtNormal = m_Values[0][nDexel][n].vtMaxN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
n -= 1 ;
|
|
if ( n >= 0)
|
|
dZ = m_Values[0][nDexel][n].dMax ;
|
|
}
|
|
}
|
|
|
|
else {
|
|
int n = 0 ;
|
|
double dZ = m_Values[0][nDexel][n].dMin ;
|
|
while ( n < nSize && dZ < dMaxZ + dEps) {
|
|
if ( dZ > dMinZ - dEps) {
|
|
ptInt.z = dZ ;
|
|
vtNormal = m_Values[0][nDexel][n].vtMinN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
n += 1 ;
|
|
if ( n < nSize)
|
|
dZ = m_Values[0][nDexel][n].dMin ;
|
|
}
|
|
}
|
|
|
|
if ( ! bFound)
|
|
ptInt.z = 0.5 * ( dMinZ + dMaxZ) ;
|
|
}
|
|
|
|
return bFound ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IsValidVoxel( int nN) const
|
|
{
|
|
// Calcolo il numero di voxel lungo X,Y e Z
|
|
unsigned int nVoxNumX = m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2) ;
|
|
unsigned int nVoxNumY = m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] / N_DEXVOXRATIO == 0 ? 1 : 2) ;
|
|
unsigned int nVoxNumZ = m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2) ;
|
|
// Verifico la validità del voxel
|
|
return ( nN >= 0 && nN < int( nVoxNumX * nVoxNumY * nVoxNumZ)) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IsValidVoxel( int nI, int nJ, int nK) const
|
|
{
|
|
// Calcolo il numero di voxel lungo X,Y e Z
|
|
int nVoxNumX = int( m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] / N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
// Verifico la validità del voxel
|
|
return ( nI >= - 1 && nI < nVoxNumX - 1 &&
|
|
nJ >= - 1 && nJ < nVoxNumY - 1 &&
|
|
nK >= - 1 && nK < nVoxNumZ - 1 ) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetVoxIJKFromN( int nN, int& nI, int& nJ, int& nK) const
|
|
{
|
|
// Calcolo il numero di voxel lungo X,Y e Z
|
|
int nVoxNumX = int( m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] / N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
|
|
// Controllo sulla validità del voxel
|
|
if ( nN < 0 || nN >= nVoxNumX * nVoxNumY * nVoxNumZ)
|
|
return false ;
|
|
|
|
// Calcolo gli indici
|
|
nI = ( nN % ( nVoxNumX * nVoxNumY)) % nVoxNumX - 1 ;
|
|
nJ = ( nN % ( nVoxNumX * nVoxNumY)) / nVoxNumX - 1 ;
|
|
nK = ( nN / ( nVoxNumX * nVoxNumY)) - 1 ;
|
|
|
|
// Controllo la sensatezza del risultato ottenuto
|
|
return ( nI >= - 1 && nI < nVoxNumX - 1 &&
|
|
nJ >= - 1 && nJ < nVoxNumY - 1 &&
|
|
nK >= - 1 && nK < nVoxNumZ - 1 ) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetVoxNFromIJK( int nI, int nJ, int nK, int& nN) const
|
|
{
|
|
// Calcolo il numero di voxel lungo X,Y e Z
|
|
int nVoxNumX = int( m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] / N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
|
|
// Controllo la validità del voxel
|
|
if ( nI <= - 2 || nI >= nVoxNumX - 1 ||
|
|
nJ <= - 2 || nJ >= nVoxNumY - 1 ||
|
|
nK <= - 2 || nK >= nVoxNumZ - 1 )
|
|
return false ;
|
|
|
|
// Calcolo il numero di Voxel
|
|
nN = nVoxNumX * nVoxNumY * ( nK + 1) +
|
|
nVoxNumX * ( nJ + 1) + nI + 1 ;
|
|
|
|
// controllo sulla sensatezza del numero ottenuto
|
|
return ( nN >= 0 && nN < nVoxNumX * nVoxNumY * nVoxNumZ) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetBlockIJKFromN( int nBlock, int nIJK[]) const
|
|
{
|
|
// Controllo sulla validità del blocco
|
|
if ( nBlock < 0 || nBlock >= int( m_nNumBlock))
|
|
return false ;
|
|
|
|
// Calcolo posizione del blocco nel reticolo
|
|
nIJK[0] = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) % int( m_nFracLin[0]) ;
|
|
nIJK[1] = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) / int( m_nFracLin[0]) ;
|
|
nIJK[2] = ( nBlock / int( m_nFracLin[0] * m_nFracLin[1])) ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetBlockNFromIJK( int nIJK[], int& nBlock) const
|
|
{
|
|
// Controllo sulla validità degli indici i, j, k del blocco
|
|
if ( nIJK[0] < 0 || nIJK[0] >= int( m_nFracLin[0]) ||
|
|
nIJK[1] < 0 || nIJK[1] >= int( m_nFracLin[1]) ||
|
|
nIJK[2] < 0 || nIJK[2] >= int( m_nFracLin[2]))
|
|
return false ;
|
|
|
|
// Determino il blocco
|
|
nBlock = m_nFracLin[0] * m_nFracLin[1] * nIJK[2] + m_nFracLin[0] * nIJK[1] + nIJK[0] ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const
|
|
{
|
|
// Controllo sulla validità degli indici i, j, k del blocco
|
|
if ( nIJK[0] < 0 || nIJK[0] >= int( m_nFracLin[0]) ||
|
|
nIJK[1] < 0 || nIJK[1] >= int( m_nFracLin[1]) ||
|
|
nIJK[2] < 0 || nIJK[2] >= int( m_nFracLin[2]))
|
|
return false ;
|
|
|
|
// Calcolo il numero di voxel lungo X,Y e Z
|
|
int nVoxNumX = int( m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] / N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
|
|
// Calcolo limiti per l'indice i
|
|
nLimits[0] = ( nIJK[0] == 0 ? - 1 : nIJK[0] * int( m_nVoxNumPerBlock)) ;
|
|
nLimits[1] = ( nIJK[0] + 1 == int( m_nFracLin[0]) ?
|
|
nVoxNumX - 1 : ( nIJK[0] + 1) * int( m_nVoxNumPerBlock)) ;
|
|
|
|
// Calcolo limiti per l'indice j
|
|
nLimits[2] = ( nIJK[1] == 0 ? - 1 : nIJK[1] * int( m_nVoxNumPerBlock)) ;
|
|
nLimits[3] = ( nIJK[1] + 1 == int( m_nFracLin[1]) ?
|
|
nVoxNumY - 1 : ( nIJK[1] + 1) * int( m_nVoxNumPerBlock)) ;
|
|
|
|
// Calcolo limiti per l'indice k
|
|
nLimits[4] = ( nIJK[2] == 0 ? - 1 : nIJK[2] * int( m_nVoxNumPerBlock)) ;
|
|
nLimits[5] = ( nIJK[2] + 1 == int( m_nFracLin[2]) ?
|
|
nVoxNumZ - 1 : ( nIJK[2] + 1) * int( m_nVoxNumPerBlock)) ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP, double dPrec) const
|
|
{
|
|
// Calcolo il numero di voxel lungo X,Y e Z
|
|
int nVoxNumX = int( m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] / N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
// Controllo sulla validità del voxel
|
|
if ( nI <= - 2 || nI >= nVoxNumX - 1 ||
|
|
nJ <= - 2 || nJ >= nVoxNumY - 1 ||
|
|
nK <= - 2 || nK >= nVoxNumZ - 1 )
|
|
return false ;
|
|
// Se la tolleranza è minore di EPS_SMALL, la considero nulla
|
|
if ( dPrec < EPS_ZERO)
|
|
dPrec = 0 ;
|
|
// Controllo che ogni coordinata stia dentro le relative limitazioni del voxel
|
|
bool bI = ptP.x > ( nI * N_DEXVOXRATIO + 0.5) * m_dStep - dPrec &&
|
|
ptP.x < ( ( nI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep + dPrec ;
|
|
bool bJ = ptP.y > ( nJ * N_DEXVOXRATIO + 0.5) * m_dStep - dPrec &&
|
|
ptP.y < ( ( nJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep + dPrec ;
|
|
bool bK = ptP.z > ( nK * N_DEXVOXRATIO + 0.5) * m_dStep - dPrec &&
|
|
ptP.z < ( ( nK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep + dPrec ;
|
|
|
|
return ( bI && bJ && bK) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetPointVoxel( const Point3d& ptP, int& nVoxI, int& nVoxJ, int& nVoxK) const
|
|
{
|
|
// Calcolo il numero di voxel lungo X,Y e Z
|
|
int nVoxNumX = int( m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] / N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
// Calcolo gli indici del voxel
|
|
nVoxI = int( floor( ( ptP.x - 0.5 * m_dStep) / ( m_dStep * N_DEXVOXRATIO))) ;
|
|
nVoxJ = int( floor( ( ptP.y - 0.5 * m_dStep) / ( m_dStep * N_DEXVOXRATIO))) ;
|
|
nVoxK = int( floor( ( ptP.z - 0.5 * m_dStep) / ( m_dStep * N_DEXVOXRATIO))) ;
|
|
// Controllo la validità del voxel
|
|
return ( nVoxI >= - 1 && nVoxI < nVoxNumX - 1) &&
|
|
( nVoxJ >= - 1 && nVoxJ < nVoxNumY - 1) &&
|
|
( nVoxK >= - 1 && nVoxK < nVoxNumZ - 1) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetVoxelBlockIJK( const int nVoxIJK[], int nBlockIJK[]) const
|
|
{
|
|
// Calcolo il numero di voxel lungo X,Y e Z
|
|
int nVoxNumX = int( m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] / N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
// Controllo sull'ammissibilità del voxel
|
|
if ( nVoxIJK[0] <= - 2 || nVoxIJK[0] >= nVoxNumX - 1 ||
|
|
nVoxIJK[1] <= - 2 || nVoxIJK[1] >= nVoxNumY - 1 ||
|
|
nVoxIJK[2] <= - 2 || nVoxIJK[2] >= nVoxNumZ - 1 )
|
|
return false ;
|
|
// Divisioni intere
|
|
int nIntRatio0 = nVoxIJK[0] / m_nVoxNumPerBlock ;
|
|
int nIntRatio1 = nVoxIJK[1] / m_nVoxNumPerBlock ;
|
|
int nIntRatio2 = nVoxIJK[2] / m_nVoxNumPerBlock ;
|
|
// Calcolo indici del blocco
|
|
nBlockIJK[0] = ( nVoxIJK[0] == -1 ? 0 : max( 0, nIntRatio0 - ( nIntRatio0 == m_nFracLin[0] ? 1 : 0))) ;
|
|
nBlockIJK[1] = ( nVoxIJK[1] == -1 ? 0 : max( 0, nIntRatio1 - ( nIntRatio1 == m_nFracLin[1] ? 1 : 0))) ;
|
|
nBlockIJK[2] = ( nVoxIJK[2] == -1 ? 0 : max( 0, nIntRatio2 - ( nIntRatio2 == m_nFracLin[2] ? 1 : 0))) ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetAdjBlockToBlock( int nBlockN, int nDeltaI, int nDeltaJ, int nDeltaK, int& nAdjBlockN) const
|
|
{
|
|
// Test sulla validità degli incrementi su i,j,k
|
|
if ( nDeltaI < - 1 || nDeltaI > 1 ||
|
|
nDeltaJ < - 1 || nDeltaJ > 1 ||
|
|
nDeltaK < - 1 || nDeltaK > 1)
|
|
return false ;
|
|
|
|
// Determino blocco adiacente
|
|
nAdjBlockN = nBlockN ;
|
|
nAdjBlockN += nDeltaI ;
|
|
nAdjBlockN += nDeltaJ * m_nFracLin[0] ;
|
|
nAdjBlockN += nDeltaK * m_nFracLin[0] * m_nFracLin[1] ;
|
|
|
|
// Se il blocco adiacente esiste restituisco vero, altrimenti falso.
|
|
return ( nAdjBlockN > -1 && nAdjBlockN < int( m_nNumBlock)) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], bool bType) const
|
|
{
|
|
// Calcolo il numero di voxel lungo X,Y e Z
|
|
int nVoxNumX = int( m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] / N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
|
|
|
|
// Test sulla validità dei limiti
|
|
if ( nLimits[0] < - 1 || nLimits[0] > nVoxNumX - 1 ||
|
|
nLimits[1] < - 1 || nLimits[1] > nVoxNumX - 1 ||
|
|
nLimits[2] < - 1 || nLimits[2] > nVoxNumY - 1 ||
|
|
nLimits[3] < - 1 || nLimits[3] > nVoxNumY - 1 ||
|
|
nLimits[4] < - 1 || nLimits[4] > nVoxNumZ - 1 ||
|
|
nLimits[5] < - 1 || nLimits[5] > nVoxNumZ - 1 )
|
|
return false ;
|
|
|
|
// Controllo sull'ammissibilità del voxel
|
|
if ( nIJK[0] <= -2 || nIJK[0] > nVoxNumX - 2 ||
|
|
nIJK[1] <= -2 || nIJK[1] > nVoxNumY - 2 ||
|
|
nIJK[2] <= -2 || nIJK[2] > nVoxNumZ - 2)
|
|
return false ;
|
|
|
|
// Se cerchiamo i voxel che sono sulla frontiera del blocco
|
|
if ( ! bType) {
|
|
return ( nIJK[0] == nLimits[0] || nIJK[0] == nLimits[1] - 1 ||
|
|
nIJK[1] == nLimits[2] || nIJK[1] == nLimits[3] - 1 ||
|
|
nIJK[2] == nLimits[4] || nIJK[2] == nLimits[5] - 1) ;
|
|
}
|
|
|
|
// Altrimenti cerchiamo i voxel che confinano con quelli di altri blocchi
|
|
// Condizione necessaria è che il voxel sia sulla frontiera
|
|
if ( nIJK[0] == nLimits[0] || nIJK[0] == nLimits[1] - 1 ||
|
|
nIJK[1] == nLimits[2] || nIJK[1] == nLimits[3] - 1 ||
|
|
nIJK[2] == nLimits[4] || nIJK[2] == nLimits[5] - 1) {
|
|
// Indici del blocco
|
|
int nCurrentBlockIJK[3] ;
|
|
GetVoxelBlockIJK( nIJK, nCurrentBlockIJK) ;
|
|
// Ciclo sulle possibili posizioni dei voxel adiacenti
|
|
for ( int nDeltaI = -1 ; nDeltaI <= 1 ; ++ nDeltaI) {
|
|
for ( int nDeltaJ = -1 ; nDeltaJ <= 1 ; ++ nDeltaJ) {
|
|
for ( int nDeltaK = -1 ; nDeltaK <= 1 ; ++ nDeltaK) {
|
|
// Evito controllo con se stesso
|
|
if ( nDeltaI == 0 && nDeltaJ == 0 && nDeltaK == 0)
|
|
continue ;
|
|
// Indici del voxel adiacente
|
|
int nAdjIJK[3] = { nIJK[0] + nDeltaI, nIJK[1] + nDeltaJ, nIJK[2] + nDeltaK} ;
|
|
// Determino (se esiste) il blocco in cui cade il voxel adiacente.
|
|
int nAdjBlockIJK[3] ;
|
|
if ( GetVoxelBlockIJK( nAdjIJK, nAdjBlockIJK)) {
|
|
if ( nAdjBlockIJK[0] > -1 && nAdjBlockIJK[0] < int( m_nFracLin[0]) &&
|
|
nAdjBlockIJK[1] > -1 && nAdjBlockIJK[1] < int( m_nFracLin[1]) &&
|
|
nAdjBlockIJK[2] > -1 && nAdjBlockIJK[2] < int( m_nFracLin[2]) &&
|
|
( nCurrentBlockIJK[0] != nAdjBlockIJK[0] ||
|
|
nCurrentBlockIJK[1] != nAdjBlockIJK[1] ||
|
|
nCurrentBlockIJK[2] != nAdjBlockIJK[2])) {
|
|
return true ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return false ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert,
|
|
const int nBlockLimits[], const int nVoxIJK[]) const
|
|
{
|
|
// Se il triangolo ha area nulla, esco
|
|
if ( trTria.GetArea() < EPS_SMALL)
|
|
return false ;
|
|
|
|
// Determino i punti del triangolo sulla griglia
|
|
int nNotVert[2] ;
|
|
int nCount = 0 ;
|
|
for ( int nV = 0 ; nV < 3 ; ++ nV) {
|
|
if ( SqDist( trTria.GetP( nV), ptVert) > SQ_EPS_SMALL) {
|
|
if ( nCount == 0)
|
|
nNotVert[nCount] = nV ;
|
|
else if ( nCount == 1)
|
|
nNotVert[nCount] = nV ;
|
|
nCount ++ ;
|
|
}
|
|
}
|
|
|
|
// Se non ho trovato due punti, errore
|
|
if ( nCount != 2)
|
|
return false ;
|
|
|
|
// Punti del triangolo sulla griglia
|
|
Point3d ptFirstGrPt = trTria.GetP( nNotVert[0]) ;
|
|
Point3d ptSecondGrPt = trTria.GetP( nNotVert[1]) ;
|
|
|
|
// Verifico se tali punti sono sulla griglia
|
|
for ( int nC = 0 ; nC < 3 ; ++ nC) {
|
|
if ( nVoxIJK[nC] == nBlockLimits[2*nC]) {
|
|
double dGrid = ( nBlockLimits[2*nC] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL &&
|
|
abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL)
|
|
return true ;
|
|
}
|
|
if ( nVoxIJK[nC] + 1 == nBlockLimits[2*nC+1]) {
|
|
double dGrid = ( nBlockLimits[2*nC+1] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL &&
|
|
abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL)
|
|
return true ;
|
|
}
|
|
}
|
|
|
|
return false ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::ProcessVoxContXY( VoxelContainer& VoxContXY, bool bPlus, TRIA3DLIST& lstTria) const
|
|
{
|
|
VoxelContainer::iterator it = VoxContXY.begin() ;
|
|
while ( it != VoxContXY.end()) {
|
|
|
|
int nN = ( *it).first ;
|
|
int nI, nJ, nK ;
|
|
GetVoxIJKFromN( nN, nI, nJ, nK) ;
|
|
|
|
// Costruzione del primo rettangolo: un singolo voxel
|
|
int nMinI = nI ;
|
|
int nMinJ = nJ ;
|
|
int nMaxI = nI ;
|
|
int nMaxJ = nJ ;
|
|
double dCordZ = ( *it).second ;
|
|
|
|
// Flag sul ritrovamento di un rettangolo più grande.
|
|
bool bOkI = true ;
|
|
bool bOkJ = true ;
|
|
while ( bOkI || bOkJ) {
|
|
// Se precedente espansione ok e lato I più lungo o non ok J
|
|
if ( bOkI && ( nMaxI - nMinI >= nMaxJ - nMinJ || ! bOkJ)) {
|
|
// Analizzo linea superiore
|
|
bool bSupJ = true ;
|
|
for ( int i = nMinI ; bSupJ && i <= nMaxI ; ++ i) {
|
|
if ( ! Find( VoxContXY, i, nMaxJ + 1, nK, dCordZ))
|
|
bSupJ = false ;
|
|
}
|
|
// Se linea superiore accettata, elimino i voxel
|
|
if ( bSupJ) {
|
|
for ( int i = nMinI ; i <= nMaxI ; ++ i)
|
|
Remove( VoxContXY, i, nMaxJ + 1, nK) ;
|
|
++ nMaxJ ;
|
|
}
|
|
// Analizzo linea inferiore
|
|
bool bInfJ = true ;
|
|
for ( int i = nMinI ; bInfJ && i <= nMaxI ; ++ i) {
|
|
if ( ! Find( VoxContXY, i, nMinJ - 1, nK, dCordZ))
|
|
bInfJ = false ;
|
|
}
|
|
// Se linea inferiore accettata, elimino i voxel
|
|
if ( bInfJ) {
|
|
for ( int i = nMinI ; i <= nMaxI ; ++ i)
|
|
Remove( VoxContXY, i, nMinJ - 1, nK) ;
|
|
-- nMinJ ;
|
|
}
|
|
bOkI = bSupJ || bInfJ ;
|
|
}
|
|
// altrimenti se precedente espansione J ok
|
|
else if ( bOkJ) {
|
|
// Analizzo linea destra
|
|
bool bSupI = true ;
|
|
for ( int j = nMinJ ; bSupI && j <= nMaxJ ; ++ j) {
|
|
if ( ! Find( VoxContXY, nMaxI + 1, j, nK, dCordZ))
|
|
bSupI = false ;
|
|
}
|
|
// Se linea destra accettata, elimino i voxel
|
|
if ( bSupI) {
|
|
for ( int j = nMinJ ; j <= nMaxJ ; ++ j)
|
|
Remove( VoxContXY, nMaxI + 1, j, nK) ;
|
|
++ nMaxI ;
|
|
}
|
|
// Analizzo linea sinistra
|
|
bool bInfI = true ;
|
|
for ( int j = nMinJ ; bInfI && j <= nMaxJ ; ++ j) {
|
|
if ( ! Find( VoxContXY, nMinI - 1, j, nK, dCordZ))
|
|
bInfI = false ;
|
|
}
|
|
// Se linea sinistra accettata, elimino i voxel
|
|
if ( bInfI) {
|
|
for ( int j = nMinJ ; j <= nMaxJ ; ++ j)
|
|
Remove( VoxContXY, nMinI - 1, j, nK) ;
|
|
-- nMinI ;
|
|
}
|
|
bOkJ = bSupI || bInfI ;
|
|
}
|
|
}
|
|
|
|
Point3d ptT0( ( nMinI * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nMinJ * N_DEXVOXRATIO + 0.5) * m_dStep, dCordZ) ;
|
|
Point3d ptT1( ( ( nMaxI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( nMinJ * N_DEXVOXRATIO + 0.5) * m_dStep), dCordZ) ;
|
|
Point3d ptT2( ( ( nMaxI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( ( nMaxJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep), dCordZ) ;
|
|
Point3d ptT3( ( nMinI * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( ( nMaxJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep), dCordZ) ;
|
|
|
|
ptT0.ToGlob( m_MapFrame) ;
|
|
ptT1.ToGlob( m_MapFrame) ;
|
|
ptT2.ToGlob( m_MapFrame) ;
|
|
ptT3.ToGlob( m_MapFrame) ;
|
|
|
|
Triangle3d Tria0, Tria1 ;
|
|
if ( bPlus) {
|
|
Tria0.Set( ptT0, ptT1, ptT3) ;
|
|
Tria1.Set( ptT1, ptT2, ptT3) ;
|
|
}
|
|
else {
|
|
Tria0.Set( ptT0, ptT3, ptT1) ;
|
|
Tria1.Set( ptT1, ptT3, ptT2) ;
|
|
}
|
|
bool bV0 = Tria0.Validate( true) ;
|
|
bool bV1 = Tria1.Validate( true) ;
|
|
// Aggiungo alla lista
|
|
lstTria.emplace_back( Tria0) ;
|
|
lstTria.emplace_back( Tria1) ;
|
|
|
|
// Elimino il voxel da cui sono partito a ingrandire.
|
|
VoxContXY.erase( nN) ;
|
|
|
|
// Passo al primo voxel rimasto
|
|
it = VoxContXY.begin() ;
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::ProcessVoxContYZ( VoxelContainer& VoxContYZ, bool bPlus, TRIA3DLIST& lstTria) const
|
|
{
|
|
VoxelContainer::iterator it = VoxContYZ.begin() ;
|
|
while ( it != VoxContYZ.end()) {
|
|
|
|
int nN = ( *it).first ;
|
|
int nI, nJ, nK ;
|
|
GetVoxIJKFromN( nN, nI, nJ, nK) ;
|
|
|
|
// Costruzione del primo rettangolo: un singolo voxel
|
|
int nMinJ = nJ ;
|
|
int nMinK = nK ;
|
|
int nMaxJ = nJ ;
|
|
int nMaxK = nK ;
|
|
double dCordX = ( *it).second ;
|
|
|
|
// Flag sul ritrovamento di un rettangolo più grande.
|
|
bool bOkJ = true ;
|
|
bool bOkK = true ;
|
|
while ( bOkJ || bOkK) {
|
|
// Se precedente espansione ok e lato J più lungo o non ok K
|
|
if ( bOkJ && ( nMaxJ - nMinJ >= nMaxK - nMinK || ! bOkK)) {
|
|
// Analizzo linea superiore
|
|
bool bSupK = true ;
|
|
for ( int j = nMinJ ; bSupK && j <= nMaxJ ; ++ j) {
|
|
if ( ! Find( VoxContYZ, nI, j, nMaxK + 1, dCordX))
|
|
bSupK = false ;
|
|
}
|
|
// Se linea superiore accettata, elimino i voxel
|
|
if ( bSupK) {
|
|
for ( int j = nMinJ ; j <= nMaxJ ; ++ j)
|
|
Remove( VoxContYZ, nI, j, nMaxK + 1) ;
|
|
++ nMaxK ;
|
|
}
|
|
// Analizzo linea inferiore
|
|
bool bInfK = true ;
|
|
for ( int j = nMinJ ; bInfK && j <= nMaxJ ; ++ j) {
|
|
if ( ! Find( VoxContYZ, nI, j, nMinK - 1, dCordX))
|
|
bInfK = false ;
|
|
}
|
|
// Se linea inferiore accettata, elimino i voxel
|
|
if ( bInfK) {
|
|
for ( int j = nMinJ ; j <= nMaxJ ; ++ j)
|
|
Remove( VoxContYZ, nI, j, nMinK - 1) ;
|
|
-- nMinK ;
|
|
}
|
|
bOkJ = bSupK || bInfK ;
|
|
}
|
|
// altrimenti se precedente espansione K ok
|
|
else if ( bOkK) {
|
|
// Analizzo linea destra
|
|
bool bSupJ = true ;
|
|
for ( int k = nMinK ; bSupJ && k <= nMaxK ; ++ k) {
|
|
if ( ! Find( VoxContYZ, nI, nMaxJ + 1, k, dCordX))
|
|
bSupJ = false ;
|
|
}
|
|
// Se linea destra accettata, elimino i voxel
|
|
if ( bSupJ) {
|
|
for ( int k = nMinK ; k <= nMaxK ; ++ k)
|
|
Remove( VoxContYZ, nI, nMaxJ + 1, k) ;
|
|
++ nMaxJ ;
|
|
}
|
|
// Analizzo linea sinistra
|
|
bool bInfJ = true ;
|
|
for ( int k = nMinK ; bInfJ && k <= nMaxK ; ++ k) {
|
|
if ( ! Find( VoxContYZ, nI, nMinJ - 1, k, dCordX))
|
|
bInfJ = false ;
|
|
}
|
|
// Se linea sinistra accettata, elimino i voxel
|
|
if ( bInfJ) {
|
|
for ( int k = nMinK ; k <= nMaxK ; ++ k)
|
|
Remove( VoxContYZ, nI, nMinJ - 1, k) ;
|
|
-- nMinJ ;
|
|
}
|
|
bOkK = bSupJ || bInfJ ;
|
|
}
|
|
}
|
|
|
|
Point3d ptT0( dCordX, ( nMinJ * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( nMinK * N_DEXVOXRATIO + 0.5) * m_dStep)) ;
|
|
Point3d ptT1( dCordX, ( ( nMaxJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( nMinK * N_DEXVOXRATIO + 0.5) * m_dStep)) ;
|
|
Point3d ptT2( dCordX, ( ( nMaxJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( ( nMaxK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep)) ;
|
|
Point3d ptT3( dCordX, ( nMinJ * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( ( nMaxK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep)) ;
|
|
|
|
ptT0.ToGlob( m_MapFrame) ;
|
|
ptT1.ToGlob( m_MapFrame) ;
|
|
ptT2.ToGlob( m_MapFrame) ;
|
|
ptT3.ToGlob( m_MapFrame) ;
|
|
|
|
Triangle3d Tria0, Tria1 ;
|
|
if ( bPlus) {
|
|
Tria0.Set( ptT0, ptT1, ptT3) ;
|
|
Tria1.Set( ptT1, ptT2, ptT3) ;
|
|
}
|
|
else {
|
|
Tria0.Set( ptT0, ptT3, ptT1) ;
|
|
Tria1.Set( ptT1, ptT3, ptT2) ;
|
|
}
|
|
bool bV0 = Tria0.Validate( true) ;
|
|
bool bV1 = Tria1.Validate( true) ;
|
|
// Aggiungo alla lista
|
|
lstTria.emplace_back( Tria0) ;
|
|
lstTria.emplace_back( Tria1) ;
|
|
|
|
// Elimino il voxel da cui sono partito a ingrandire.
|
|
VoxContYZ.erase( nN) ;
|
|
|
|
// Passo al primo voxel rimasto
|
|
it = VoxContYZ.begin() ;
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::ProcessVoxContXZ( VoxelContainer& VoxContXZ, bool bPlus, TRIA3DLIST& lstTria) const
|
|
{
|
|
VoxelContainer::iterator it = VoxContXZ.begin() ;
|
|
while ( it != VoxContXZ.end()) {
|
|
|
|
int nN = ( *it).first ;
|
|
int nI, nJ, nK ;
|
|
GetVoxIJKFromN( nN, nI, nJ, nK) ;
|
|
|
|
// Costruzione del primo rettangolo: un singolo voxel
|
|
int nMinI = nI ;
|
|
int nMinK = nK ;
|
|
int nMaxI = nI ;
|
|
int nMaxK = nK ;
|
|
double dCordY = ( *it).second ;
|
|
|
|
// Flag sul ritrovamento di un rettangolo più grande.
|
|
bool bOkI = true ;
|
|
bool bOkK = true ;
|
|
while ( bOkI || bOkK) {
|
|
// Se lato I più lungo e precedente espansione ok
|
|
if ( bOkI && ( nMaxI - nMinI >= nMaxK - nMinK || ! bOkK)) {
|
|
// Analizzo linea superiore
|
|
bool bSupK = true ;
|
|
for ( int i = nMinI ; bSupK && i <= nMaxI ; ++ i) {
|
|
if ( ! Find( VoxContXZ, i, nJ, nMaxK + 1, dCordY))
|
|
bSupK = false ;
|
|
}
|
|
// Se linea superiore accettata, elimino i voxel
|
|
if ( bSupK) {
|
|
for ( int i = nMinI ; i <= nMaxI ; ++ i)
|
|
Remove( VoxContXZ, i, nJ, nMaxK + 1) ;
|
|
++ nMaxK ;
|
|
}
|
|
// Analizzo linea inferiore
|
|
bool bInfK = true ;
|
|
for ( int i = nMinI ; bInfK && i <= nMaxI ; ++ i) {
|
|
if ( ! Find( VoxContXZ, i, nJ, nMinK - 1, dCordY))
|
|
bInfK = false ;
|
|
}
|
|
// Se linea inferiore accettata, elimino i voxel
|
|
if ( bInfK) {
|
|
for ( int i = nMinI ; i <= nMaxI ; ++ i)
|
|
Remove( VoxContXZ, i, nJ, nMinK - 1) ;
|
|
-- nMinK ;
|
|
}
|
|
bOkI = bSupK || bInfK ;
|
|
}
|
|
// altrimenti se precedente espansione K ok
|
|
else if ( bOkK) {
|
|
// Analizzo linea destra
|
|
bool bSupI = true ;
|
|
for ( int k = nMinK ; bSupI && k <= nMaxK ; ++ k) {
|
|
if ( ! Find( VoxContXZ, nMaxI + 1, nJ, k, dCordY))
|
|
bSupI = false ;
|
|
}
|
|
// Se linea destra accettata, elimino i voxel
|
|
if ( bSupI) {
|
|
for ( int k = nMinK ; k <= nMaxK ; ++ k)
|
|
Remove( VoxContXZ, nMaxI + 1, nJ, k) ;
|
|
++ nMaxI ;
|
|
}
|
|
// Analizzo linea sinistra
|
|
bool bInfI = true ;
|
|
for ( int k = nMinK ; bInfI && k <= nMaxK ; ++ k) {
|
|
if ( ! Find( VoxContXZ, nMinI - 1, nJ, k, dCordY))
|
|
bInfI = false ;
|
|
}
|
|
// Se linea sinistra accettata, elimino i voxel
|
|
if ( bInfI) {
|
|
for ( int k = nMinK ; k <= nMaxK ; ++ k)
|
|
Remove( VoxContXZ, nMinI - 1, nJ, k) ;
|
|
-- nMinI ;
|
|
}
|
|
bOkK = bSupI || bInfI ;
|
|
}
|
|
}
|
|
|
|
Point3d ptT0( ( nMinI * N_DEXVOXRATIO + 0.5) * m_dStep, dCordY,
|
|
( ( nMinK * N_DEXVOXRATIO + 0.5) * m_dStep)) ;
|
|
Point3d ptT1( ( ( nMaxI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, dCordY,
|
|
( ( nMinK * N_DEXVOXRATIO + 0.5) * m_dStep)) ;
|
|
Point3d ptT2( ( ( nMaxI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, dCordY,
|
|
( ( ( nMaxK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep)) ;
|
|
Point3d ptT3( ( nMinI * N_DEXVOXRATIO + 0.5) * m_dStep, dCordY,
|
|
( ( ( nMaxK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep)) ;
|
|
|
|
ptT0.ToGlob( m_MapFrame) ;
|
|
ptT1.ToGlob( m_MapFrame) ;
|
|
ptT2.ToGlob( m_MapFrame) ;
|
|
ptT3.ToGlob( m_MapFrame) ;
|
|
|
|
Triangle3d Tria0, Tria1 ;
|
|
if ( bPlus) {
|
|
Tria0.Set( ptT0, ptT3, ptT1) ;
|
|
Tria1.Set( ptT1, ptT3, ptT2) ;
|
|
}
|
|
else {
|
|
Tria0.Set( ptT0, ptT1, ptT3) ;
|
|
Tria1.Set( ptT1, ptT2, ptT3) ;
|
|
}
|
|
bool bV0 = Tria0.Validate( true) ;
|
|
bool bV1 = Tria1.Validate( true) ;
|
|
// Aggiungo alla lista
|
|
lstTria.emplace_back( Tria0) ;
|
|
lstTria.emplace_back( Tria1) ;
|
|
|
|
// Elimino il voxel da cui sono partito a ingrandire.
|
|
VoxContXZ.erase( nN) ;
|
|
|
|
// Passo al primo voxel rimasto
|
|
it = VoxContXZ.begin() ;
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::Find( const VoxelContainer& VoxCont, int nI, int nJ, int nK, double dPos) const
|
|
{
|
|
// indice globale del voxel
|
|
int nN ;
|
|
if ( ! GetVoxNFromIJK( nI, nJ, nK, nN))
|
|
return false ;
|
|
// cerco il voxel nel contenitore
|
|
auto iter = VoxCont.find( nN) ;
|
|
return ( iter != VoxCont.end() && abs( dPos - ( *iter).second) < EPS_SMALL) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::Remove( VoxelContainer& VoxCont, int nI, int nJ, int nK) const
|
|
{
|
|
int nN ;
|
|
if ( ! GetVoxNFromIJK( nI, nJ, nK, nN))
|
|
return false ;
|
|
return ( VoxCont.erase( nN) > 0) ;
|
|
}
|