Dario Sassi
4644 lines
223 KiB
C++
4644 lines
223 KiB
C++
//----------------------------------------------------------------------------
|
|
// EgalTech 2015-2019
|
|
//----------------------------------------------------------------------------
|
|
// File : VolZmap.cpp Data : 07.03.19 Versione : 2.1c2
|
|
// 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 "MC_Tables.h"
|
|
#include "PolygonPlane.h"
|
|
#include "/EgtDev/Include/EGkIntervals.h"
|
|
#include "/EgtDev/Include/EGkStringUtils3d.h"
|
|
#include "/EgtDev/Include/EGkChainCurves.h"
|
|
#include "/EgtDev/Include/EgtNumUtils.h"
|
|
#include "/EgtDev/Extern/Eigen/Core"
|
|
#include "/EgtDev/Extern/Eigen/SVD"
|
|
|
|
using namespace std ;
|
|
|
|
// ------------------------- Tipi per SVD con Eigen ------------------------------------------------------------------------------
|
|
static const int MAX_FAN_BASE_VERTS = 7 ;
|
|
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor, MAX_FAN_BASE_VERTS, 3> SvdMatrix ;
|
|
typedef Eigen::Matrix<double, Eigen::Dynamic, 1, Eigen::ColMajor, MAX_FAN_BASE_VERTS, 1> SvdVector ;
|
|
typedef Eigen::JacobiSVD<SvdMatrix, Eigen::QRPreconditioners::ColPivHouseholderQRPreconditioner> SvdDecomposer ;
|
|
|
|
// ------------------------- 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} ;
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
Config2VertOrder( int nInd)
|
|
{
|
|
if ( nInd == 111 || nInd == 119 || nInd == 159 || nInd == 187 || nInd == 207 || nInd == 221 ||
|
|
nInd == 238 || nInd == 243 || nInd == 246 || nInd == 249 || nInd == 205)
|
|
return true ;
|
|
return false ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
Config3Duality(int nInd)
|
|
{
|
|
if (nInd == 5 || nInd == 10 || nInd == 18 || nInd == 24 || nInd == 33 || nInd == 36 ||
|
|
nInd == 66 || nInd == 72 || nInd == 80 || nInd == 129 || nInd == 132 || nInd == 160)
|
|
return true;
|
|
return false;
|
|
}
|
|
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
Config6Duality( int nInd)
|
|
{
|
|
if ( nInd == 21 || nInd == 22 || nInd == 28 || nInd == 41 || nInd == 42 || nInd == 44 ||
|
|
nInd == 52 || nInd == 56 || nInd == 67 || nInd == 69 || nInd == 73 || nInd == 81 ||
|
|
nInd == 84 || nInd == 97 || nInd == 104 || nInd == 131 || nInd == 134 || nInd == 138 ||
|
|
nInd == 146 || nInd == 148 || nInd == 162 || nInd == 168 || nInd == 193 || nInd == 194)
|
|
return true ;
|
|
return false ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
Config7Duality( int nInd)
|
|
{
|
|
if ( nInd == 26 || nInd == 37 || nInd == 74 || nInd == 82 ||
|
|
nInd == 88 || nInd == 133 || nInd == 161 || nInd == 164)
|
|
return true ;
|
|
return false ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int
|
|
TestOnNormal( const AppliedVector CompoField[], 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 = CompoField[i].vtVec * CompoField[j].vtVec ;
|
|
if ( dCurrCos < dMinCosTheta) {
|
|
dMinCosTheta = dCurrCos ;
|
|
nI = i ;
|
|
nJ = j ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Se la massima deviazione non supera il limite non è feature
|
|
const double SHARP_COS_SUP = 0.866 ;
|
|
const double SHARP_COS_INF = - 0.985 ;
|
|
if ( dMinCosTheta >= SHARP_COS_SUP || dMinCosTheta <= SHARP_COS_INF)
|
|
return NO_FEATURE ;
|
|
|
|
// Verifico se Edge o Corner
|
|
// direzione perpendicolare alle normali con massima differenza (potenziale edge)
|
|
Vector3d vtK = CompoField[nI].vtVec ^ CompoField[nJ].vtVec ;
|
|
vtK.Normalize() ;
|
|
// cerco normale con massima vicinanza al potenziale edge
|
|
double dMaxAbsCos = 0 ;
|
|
for ( int i = 0 ; i < nCompoElem ; ++ i) {
|
|
double dAbsCurrentCos = abs( CompoField[i].vtVec * vtK) ;
|
|
if ( dAbsCurrentCos > dMaxAbsCos)
|
|
dMaxAbsCos = dAbsCurrentCos ;
|
|
}
|
|
// se esiste normale diretta quasi come potenziale edge, allora corner
|
|
const double CORNER_COS = 0.39 ;
|
|
if ( dMaxAbsCos > CORNER_COS)
|
|
return CORNER ;
|
|
else
|
|
return EDGE ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
CanonicPlaneTest( const AppliedVector CompoField[], int nDir, double& dPos, int& nTool)
|
|
{
|
|
// 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 = CompoField[0].ptPApp.v[nI] ;
|
|
double dPos1 = CompoField[1].ptPApp.v[nI] ;
|
|
double dPos2 = CompoField[2].ptPApp.v[nI] ;
|
|
double dPos3 = CompoField[3].ptPApp.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 ;
|
|
}
|
|
int nDifferent = 0 ;
|
|
for ( int i = 0 ; i < 4 ; ++ i) {
|
|
if ( CompoField[i].vtVec * vtN < 0.999)
|
|
return false ;
|
|
for ( int j = i + 1 ; j < 4 ; ++ j) {
|
|
if ( CompoField[i].nPropIndex != CompoField[j].nPropIndex)
|
|
return false ;
|
|
}
|
|
}
|
|
nTool = CompoField[0].nPropIndex ;
|
|
// Superati tutti i test
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
DotTest( const AppliedVector CompoField[], 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 = CompoField[i].vtVec * CompoField[j].vtVec ;
|
|
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 += CompoField[i].vtVec ;
|
|
vtAvg /= nCompoElem ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
CreateBigTriangleXY( double dMinX, double dMaxX, double dMinY, double dMaxY, double dZ,
|
|
bool bNormZPlus, int nGrade, Triangle3dEx& trTria1, Triangle3dEx& trTria2)
|
|
{
|
|
// Punti che definiscono il perimetro del rettangolo
|
|
Point3d ptP0( dMinX, dMinY, dZ) ;
|
|
Point3d ptP1( dMinX, dMaxY, dZ) ;
|
|
Point3d ptP2( dMaxX, dMaxY, dZ) ;
|
|
Point3d ptP3( dMaxX, dMinY, dZ) ;
|
|
|
|
// Versore normale diretto come Z+
|
|
if ( bNormZPlus) {
|
|
trTria1.Set( ptP0, ptP2, ptP1) ;
|
|
trTria2.Set( ptP0, ptP3, ptP2) ;
|
|
}
|
|
// Versore normale diretto come Z-
|
|
else {
|
|
trTria1.Set( ptP0, ptP1, ptP2) ;
|
|
trTria2.Set( ptP0, ptP2, ptP3) ;
|
|
}
|
|
|
|
trTria1.SetGrade( nGrade) ;
|
|
trTria2.SetGrade( nGrade) ;
|
|
|
|
return trTria1.Validate( true) && trTria2.Validate( true) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
CreateBigTriangleXZ( double dMinX, double dMaxX, double dMinZ, double dMaxZ, double dY,
|
|
bool bNormYPlus, int nGrade, Triangle3dEx& trTria1, Triangle3dEx& trTria2)
|
|
{
|
|
// Punti che definiscono il perimetro del rettangolo
|
|
Point3d ptP0( dMinX, dY, dMinZ) ;
|
|
Point3d ptP1( dMinX, dY, dMaxZ) ;
|
|
Point3d ptP2( dMaxX, dY, dMaxZ) ;
|
|
Point3d ptP3( dMaxX, dY, dMinZ) ;
|
|
|
|
// Versore normale diretto come Y+
|
|
if ( bNormYPlus) {
|
|
trTria1.Set( ptP0, ptP1, ptP2) ;
|
|
trTria2.Set( ptP0, ptP2, ptP3) ;
|
|
}
|
|
// Versore normale diretto come Y-
|
|
else {
|
|
trTria1.Set( ptP0, ptP2, ptP1) ;
|
|
trTria2.Set( ptP0, ptP3, ptP2) ;
|
|
}
|
|
|
|
trTria1.SetGrade( nGrade) ;
|
|
trTria2.SetGrade( nGrade) ;
|
|
|
|
return trTria1.Validate( true) && trTria2.Validate( true) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
CreateBigTriangleYZ( double dMinY, double dMaxY, double dMinZ, double dMaxZ, double dX,
|
|
bool bNormXPlus, int nGrade, Triangle3dEx& trTria1, Triangle3dEx& trTria2)
|
|
{
|
|
// Punti che definiscono il perimetro del rettangolo
|
|
Point3d ptP0( dX, dMinY, dMinZ) ;
|
|
Point3d ptP1( dX, dMinY, dMaxZ) ;
|
|
Point3d ptP2( dX, dMaxY, dMaxZ) ;
|
|
Point3d ptP3( dX, dMaxY, dMinZ) ;
|
|
|
|
// Versore normale diretto come X+
|
|
if ( bNormXPlus) {
|
|
trTria1.Set( ptP0, ptP2, ptP1) ;
|
|
trTria2.Set( ptP0, ptP3, ptP2) ;
|
|
}
|
|
// Versore normale diretto come X-
|
|
else {
|
|
trTria1.Set( ptP0, ptP1, ptP2) ;
|
|
trTria2.Set( ptP0, ptP2, ptP3) ;
|
|
}
|
|
|
|
trTria1.SetGrade( nGrade) ;
|
|
trTria2.SetGrade( nGrade) ;
|
|
|
|
return trTria1.Validate( true) && trTria2.Validate( true) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::FindAdjComp( const vector<VoxelContainer>& vVecVox, int nCurBlock, int nCurVox, int nCurComp,
|
|
INTVECTOR& vAdjBlockVoxComp, INTVECTOR& vAdjBordBlockVoxComp) const
|
|
{
|
|
// Controllo sulla validità del blocco
|
|
if ( nCurBlock < 0 || nCurBlock >= int( m_nNumBlock))
|
|
return false ;
|
|
// Calcolo gli indici ijk minimi e massimi dei voxel del blocco
|
|
// Vettore indici i,j,k del blocco
|
|
int nIJK[3] ;
|
|
GetBlockIJKFromN( nCurBlock, nIJK) ;
|
|
// Vettore limiti sugli indici dei voxel del blocco
|
|
int nLimits[6] ;
|
|
GetBlockLimitsIJK( nIJK, nLimits) ;
|
|
// Voxel corrente
|
|
// Determino gli indici ijk del voxel corrente
|
|
int nCurIJK[3] ;
|
|
GetVoxIJKFromN( nCurVox, nCurIJK[0], nCurIJK[1], nCurIJK[2]) ;
|
|
bool bInnerVox = ! IsAVoxelOnBoundary( nLimits, nCurIJK, true) ;
|
|
Voxel CurVox = bInnerVox ? ( * vVecVox[nCurBlock].find( nCurVox)).second :
|
|
( * m_InterBlockVox[nCurBlock].find( nCurVox)).second ;
|
|
// Ciclo su tutti i voxel adiacenti
|
|
for ( int nI = - 1 ; nI <= 1 ; ++ nI) {
|
|
for ( int nJ = - 1 ; nJ <= 1 ; ++ nJ) {
|
|
for ( int nK = - 1 ; nK <= 1 ; ++ nK) {
|
|
// Salto il voxel stesso
|
|
if ( nI == 0 && nJ == 0 && nK == 0)
|
|
continue ;
|
|
// Indici del voxel adiacente
|
|
int nAdjVoxIJK[3] = { nCurIJK[0] + nI, nCurIJK[1] + nJ, nCurIJK[2] + nK} ;
|
|
// Risalgo all'indice del voxel adiacente (chiave) a partire dagli indici ijk
|
|
int nAdjVoxInd ;
|
|
if ( GetVoxNFromIJK( nAdjVoxIJK[0], nAdjVoxIJK[1], nAdjVoxIJK[2], nAdjVoxInd)) {
|
|
// Risalgo al blocco del voxel adiacente
|
|
int nAdjBlockIJK[3] ;
|
|
// Se tale blocco esiste
|
|
if ( GetVoxelBlockIJK( nAdjVoxIJK, nAdjBlockIJK)) {
|
|
// Risalgo all'indice del blocco adiacente dai suoi indici ijk
|
|
int nAdjBlockInd ;
|
|
GetBlockNFromIJK( nAdjBlockIJK, nAdjBlockInd) ;
|
|
// Se tale blocco è stato aggiornato, il voxel adiacente contiene componenti feature ed è interno
|
|
if ( m_BlockToUpdate[nAdjBlockInd] &&
|
|
vVecVox[nAdjBlockInd].find( nAdjVoxInd) != vVecVox[nAdjBlockInd].end()) {
|
|
// Accedo al voxel adiacente
|
|
Voxel AdjVox = ( * vVecVox[nAdjBlockInd].find( nAdjVoxInd)).second ;
|
|
// Valuto la connessione fra le componenti dei due voxels
|
|
bool bConnected = false ;
|
|
for ( int nAdjComp = 0 ; nAdjComp < AdjVox.nNumComp ; ++ nAdjComp) {
|
|
for ( int nCurV = 0 ; nCurV < CurVox.Compo[nCurComp].nVertNum ; ++ nCurV) {
|
|
int nNextCurV = ( nCurV + 1) % CurVox.Compo[nCurComp].nVertNum ;
|
|
for ( int nAdjV = 0 ; nAdjV < AdjVox.Compo[nAdjComp].nVertNum ; ++ nAdjV) {
|
|
int nVextAdjV = ( nAdjV + 1) % AdjVox.Compo[nAdjComp].nVertNum ;
|
|
Point3d ptV = CurVox.Compo[nCurComp].CompVecField[nCurV].ptPApp ;
|
|
Point3d ptVN = CurVox.Compo[nCurComp].CompVecField[nNextCurV].ptPApp ;
|
|
Point3d ptAdV = AdjVox.Compo[nAdjComp].CompVecField[nAdjV].ptPApp ;
|
|
Point3d ptAdVN = AdjVox.Compo[nAdjComp].CompVecField[nVextAdjV].ptPApp ;
|
|
if ( AreSamePointExact( ptV, ptAdVN) &&/*||*/ AreSamePointExact( ptVN, ptAdV)) {
|
|
bConnected = true ;
|
|
break ;
|
|
}
|
|
}
|
|
if ( bConnected)
|
|
break ;
|
|
}
|
|
if ( bConnected) {
|
|
vAdjBlockVoxComp.emplace_back( nAdjBlockInd) ;
|
|
vAdjBlockVoxComp.emplace_back( nAdjVoxInd) ;
|
|
vAdjBlockVoxComp.emplace_back( nAdjComp) ;
|
|
}
|
|
}
|
|
}
|
|
// A prescindere dal fatto che il blocco sia stato aggiornato,
|
|
// se il voxel è di frontiera lo analizzo.
|
|
if ( m_InterBlockVox[nAdjBlockInd].find( nAdjVoxInd) != m_InterBlockVox[nAdjBlockInd].end()) {
|
|
// Accedo al voxel adiacente
|
|
Voxel AdjVox = ( * m_InterBlockVox[nAdjBlockInd].find( nAdjVoxInd)).second ;
|
|
// Valuto la connessione fra le componenti dei due voxels
|
|
bool bConnected = false ;
|
|
for ( int nAdjComp = 0 ; nAdjComp < AdjVox.nNumComp ; ++ nAdjComp) {
|
|
for ( int nCurV = 0 ; nCurV < CurVox.Compo[nCurComp].nVertNum ; ++ nCurV) {
|
|
int nNextCurV = ( nCurV + 1) % CurVox.Compo[nCurComp].nVertNum ;
|
|
for ( int nAdjV = 0 ; nAdjV < AdjVox.Compo[nAdjComp].nVertNum ; ++ nAdjV) {
|
|
int nVextAdjV = ( nAdjV + 1) % AdjVox.Compo[nAdjComp].nVertNum ;
|
|
Point3d ptV = CurVox.Compo[nCurComp].CompVecField[nCurV].ptPApp ;
|
|
Point3d ptVN = CurVox.Compo[nCurComp].CompVecField[nNextCurV].ptPApp ;
|
|
Point3d ptAdV = AdjVox.Compo[nAdjComp].CompVecField[nAdjV].ptPApp ;
|
|
Point3d ptAdVN = AdjVox.Compo[nAdjComp].CompVecField[nVextAdjV].ptPApp ;
|
|
if ( AreSamePointExact( ptV, ptAdVN) &&/*||*/ AreSamePointExact( ptVN, ptAdV)) {
|
|
bConnected = true ;
|
|
break ;
|
|
}
|
|
}
|
|
if ( bConnected)
|
|
break ;
|
|
}
|
|
if ( bConnected) {
|
|
vAdjBordBlockVoxComp.emplace_back( nAdjBlockInd) ;
|
|
vAdjBordBlockVoxComp.emplace_back( nAdjVoxInd) ;
|
|
vAdjBordBlockVoxComp.emplace_back( nAdjComp) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
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::GetBlockTriangles( int nBlock, TRIA3DEXVECTOR& vTria) const
|
|
{
|
|
// Dexel (singola mappa)
|
|
if ( m_nMapNum == 1) {
|
|
// se il blocco non esiste, errore
|
|
if ( nBlock < 0 || nBlock > m_nNumBlock - 1)
|
|
return false ;
|
|
// lancio eventuale aggiornamento pendente della grafica
|
|
UpdateSingleMapGraphics() ;
|
|
// copio i triangoli
|
|
vTria = m_SingleMapTria[nBlock] ;
|
|
return true ;
|
|
}
|
|
|
|
// Tridexel (tre mappe)
|
|
else {
|
|
// se il blocco non esiste, errore
|
|
if ( nBlock < 0 || nBlock > ( m_nNumBlock + 1) - 1)
|
|
return false ;
|
|
// lancio eventuale aggiornamento pendente della grafica
|
|
UpdateTripleMapGraphics() ;
|
|
// copio i triangoli del blocco
|
|
if ( nBlock != m_nNumBlock) {
|
|
vTria.clear() ;
|
|
vTria.reserve( 10000) ;
|
|
// triangoli smooth
|
|
for ( int nVx = 0 ; nVx < int( m_BlockSmoothTria[nBlock].size()) ; ++ nVx) {
|
|
for ( int nTr = 0 ; nTr < int( m_BlockSmoothTria[nBlock][nVx].vTria.size()) ; ++ nTr)
|
|
vTria.emplace_back( m_BlockSmoothTria[nBlock][nVx].vTria[nTr]) ;
|
|
}
|
|
// triangoli grandi piatti
|
|
for ( int tBl = 0 ; tBl < int( m_BlockBigTria[nBlock].size()) ; ++ tBl) {
|
|
vTria.emplace_back( m_BlockBigTria[nBlock][tBl]) ;
|
|
}
|
|
// triangoli di feature nel blocco (ciclo sui voxel del blocco)
|
|
for ( int t1 = 0 ; t1 < int( m_BlockSharpTria[nBlock].size()) ; ++ t1) {
|
|
// ciclo sulle componenti connesse del voxel
|
|
for ( int t2 = 0 ; t2 < int( m_BlockSharpTria[nBlock][t1].vCompoTria.size()) ; ++ t2) {
|
|
// ciclo sui triangoli delle componenti connesse
|
|
for ( int t3 = 0 ; t3 < int( m_BlockSharpTria[nBlock][t1].vCompoTria[t2].size()) ; ++ t3) {
|
|
// aggiungo triangolo alla lista
|
|
vTria.emplace_back( m_BlockSharpTria[nBlock][t1].vCompoTria[t2][t3]) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
// blocco speciale aggiunto per i triangoli feature a cavallo dei blocchi
|
|
else {
|
|
vTria.clear() ;
|
|
vTria.reserve( 1000) ;
|
|
for ( int t = 0 ; t < int( m_InterBlockSharpTria.size()) ; ++ t) {
|
|
for ( int t1 = 0 ; t1 < int( m_InterBlockSharpTria[t].size()) ; ++ t1) {
|
|
for ( int t2 = 0 ; t2 < int( m_InterBlockSharpTria[t][t1].vCompoTria.size()) ; ++ t2) {
|
|
for ( int t3 = 0 ; t3 < int( m_InterBlockSharpTria[t][t1].vCompoTria[t2].size()) ; ++ t3) {
|
|
if ( m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3].GetArea() > SQ_EPS_SMALL) {
|
|
// aggiungo triangolo alla lista
|
|
vTria.emplace_back( m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3]) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return true ;
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int
|
|
VolZmap::GetBlockCount( void) const
|
|
{
|
|
return ( m_nNumBlock + ( m_nMapNum == 1 ? 0 : 1)) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int
|
|
VolZmap::GetBlockUpdatingCounter( int nBlock) const
|
|
{
|
|
// se il blocco non esiste, errore
|
|
if ( nBlock < 0 || nBlock >= int( m_BlockUpdatingCounter.size()))
|
|
return -1 ;
|
|
// Lancio eventuale aggiornamento pendente della grafica
|
|
if ( m_nMapNum == 1)
|
|
UpdateSingleMapGraphics() ;
|
|
else
|
|
UpdateTripleMapGraphics() ;
|
|
// restituisco il suo indice di aggiornamento
|
|
return m_BlockUpdatingCounter[nBlock] ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::UpdateSingleMapGraphics( void) const
|
|
{
|
|
const int MAX_DIM_CHUNK = 128 ;
|
|
|
|
// Ciclo sui blocchi
|
|
for ( int t = 0 ; t < m_nNumBlock ; ++ t) {
|
|
|
|
// Se il blocco deve essere aggiornato, eseguo
|
|
if ( m_BlockToUpdate[t]) {
|
|
|
|
m_SingleMapTria[t].clear() ;
|
|
|
|
// 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, t) ;
|
|
}
|
|
}
|
|
|
|
m_BlockToUpdate[t] = false ;
|
|
++ m_BlockUpdatingCounter[t] ;
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, int nDimChk, int nBlock) 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, nBlock) ;
|
|
}
|
|
// 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, nBlock) ;
|
|
}
|
|
}
|
|
}
|
|
// 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, nBlock) ;
|
|
}
|
|
}
|
|
}
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::CalcChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, int nBlock) 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
|
|
m_SingleMapTria[nBlock].emplace_back() ;
|
|
m_SingleMapTria[nBlock].back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_MapFrame.VersZ()) ;
|
|
m_SingleMapTria[nBlock].emplace_back() ;
|
|
m_SingleMapTria[nBlock].back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_MapFrame.VersZ()) ;
|
|
// faccia inferiore P1b->P4b->P3b->P2b : sempre visibile
|
|
m_SingleMapTria[nBlock].emplace_back() ;
|
|
m_SingleMapTria[nBlock].back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_MapFrame.VersZ()) ;
|
|
m_SingleMapTria[nBlock].emplace_back() ;
|
|
m_SingleMapTria[nBlock].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(), nBlock) ;
|
|
}
|
|
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(), nBlock) ;
|
|
}
|
|
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(), nBlock) ;
|
|
}
|
|
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(), nBlock) ;
|
|
}
|
|
//
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::CalcDexelPrisms( int nPos1, int nPos2, int nBlock) 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
|
|
m_SingleMapTria[nBlock].emplace_back() ;
|
|
m_SingleMapTria[nBlock].back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_MapFrame.VersZ()) ;
|
|
m_SingleMapTria[nBlock].emplace_back() ;
|
|
m_SingleMapTria[nBlock].back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_MapFrame.VersZ()) ;
|
|
// faccia inferiore P1b->P4b->P3b->P2b : sempre visibile
|
|
m_SingleMapTria[nBlock].emplace_back() ;
|
|
m_SingleMapTria[nBlock].back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_MapFrame.VersZ()) ;
|
|
m_SingleMapTria[nBlock].emplace_back() ;
|
|
m_SingleMapTria[nBlock].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(), nBlock) ;
|
|
int nPosNord = ( nPos2 < int( m_nNy[0] - 1) ? nPos + m_nNx[0] : - 1) ;
|
|
AddDexelSideFace( nPos, nPosNord, ptP3, ptP4, m_MapFrame.VersZ(), m_MapFrame.VersY(), nBlock) ;
|
|
int nPosWest = ( nPos1 > 0 ? nPos - 1 : - 1) ;
|
|
AddDexelSideFace( nPos, nPosWest, ptP4, ptP1, m_MapFrame.VersZ(), - m_MapFrame.VersX(), nBlock) ;
|
|
int nPosSud = ( nPos2 > 0 ? nPos - m_nNx[0] : - 1) ;
|
|
AddDexelSideFace( nPos, nPosSud, ptP1, ptP2, m_MapFrame.VersZ(), - m_MapFrame.VersY(), nBlock) ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Point3d& ptQ,
|
|
const Vector3d& vtZ, const Vector3d& vtNorm, int nBlock) 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 ;
|
|
m_SingleMapTria[nBlock].emplace_back() ;
|
|
m_SingleMapTria[nBlock].back().Set( ptP + vtDZb, ptQ + vtDZb, ptQ + vtDZt, vtNorm) ;
|
|
m_SingleMapTria[nBlock].emplace_back() ;
|
|
m_SingleMapTria[nBlock].back().Set( ptQ + vtDZt, ptP + vtDZt, ptP + vtDZb, vtNorm) ;
|
|
bFound = intFace.GetNext( dMin, dMax) ;
|
|
}
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::UpdateTripleMapGraphics( void) const
|
|
{
|
|
// se non ci sono blocchi da aggiornare, esco
|
|
if ( find( begin( m_BlockToUpdate), end( m_BlockToUpdate), true) == end( m_BlockToUpdate))
|
|
return true ;
|
|
|
|
// contenitori temporanei
|
|
vector<VoxelContainer> vVoxContainerVec ;
|
|
vVoxContainerVec.resize( m_nNumBlock) ;
|
|
SharpTriaMatrix VecTriHold ;
|
|
VecTriHold.resize( m_nNumBlock) ;
|
|
|
|
// Ciclo sui blocchi per eliminare le slice fra blocchi da aggiornare
|
|
for ( int t = 0 ; t < m_nNumBlock ; ++ t) {
|
|
for ( auto it = m_SliceXY[t].begin() ; it != m_SliceXY[t].end() ;) {
|
|
int nSlIJK[3] ;
|
|
if ( GetVoxIJKFromN( it->first, nSlIJK[0], nSlIJK[1], nSlIJK[2])) {
|
|
int nBlockIJK[3] ;
|
|
if ( GetVoxelBlockIJK( nSlIJK, nBlockIJK)) {
|
|
int nLimits[6] ;
|
|
int nDeltaIndex[3] ;
|
|
if ( GetBlockLimitsIJK( nBlockIJK, nLimits) &&
|
|
IsAVoxelOnBoundary( nLimits, nSlIJK, nDeltaIndex)) {
|
|
for ( int nInd = 0 ; nInd < 3 ; ++ nInd)
|
|
nSlIJK[nInd] += nDeltaIndex[nInd] ;
|
|
int nAdBlockIJK[3] ;
|
|
int nAdBlockNum ;
|
|
if ( GetVoxelBlockIJK( nSlIJK, nAdBlockIJK) &&
|
|
GetBlockNFromIJK( nAdBlockIJK, nAdBlockNum) &&
|
|
m_BlockToUpdate[nAdBlockNum]) {
|
|
it = m_SliceXY[t].erase( it) ;
|
|
continue ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
++ it ;
|
|
}
|
|
for ( auto it = m_SliceXZ[t].begin() ; it != m_SliceXZ[t].end() ;) {
|
|
int nSlIJK[3] ;
|
|
if ( GetVoxIJKFromN( it->first, nSlIJK[0], nSlIJK[1], nSlIJK[2])) {
|
|
int nBlockIJK[3] ;
|
|
if ( GetVoxelBlockIJK( nSlIJK, nBlockIJK)) {
|
|
int nLimits[6] ;
|
|
int nDeltaIndex[3] ;
|
|
if ( GetBlockLimitsIJK( nBlockIJK, nLimits) &&
|
|
IsAVoxelOnBoundary( nLimits, nSlIJK, nDeltaIndex)) {
|
|
for ( int nInd = 0 ; nInd < 3 ; ++ nInd)
|
|
nSlIJK[nInd] += nDeltaIndex[nInd] ;
|
|
int nAdBlockIJK[3] ;
|
|
int nAdBlockNum ;
|
|
if ( GetVoxelBlockIJK( nSlIJK, nAdBlockIJK) &&
|
|
GetBlockNFromIJK( nAdBlockIJK, nAdBlockNum) &&
|
|
m_BlockToUpdate[nAdBlockNum]) {
|
|
it = m_SliceXZ[t].erase( it) ;
|
|
continue ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
++ it ;
|
|
}
|
|
for ( auto it = m_SliceYZ[t].begin() ; it != m_SliceYZ[t].end() ;) {
|
|
int nSlIJK[3] ;
|
|
if ( GetVoxIJKFromN( it->first, nSlIJK[0], nSlIJK[1], nSlIJK[2])) {
|
|
int nBlockIJK[3] ;
|
|
if ( GetVoxelBlockIJK( nSlIJK, nBlockIJK)) {
|
|
int nLimits[6] ;
|
|
int nDeltaIndex[3] ;
|
|
if ( GetBlockLimitsIJK( nBlockIJK, nLimits) &&
|
|
IsAVoxelOnBoundary( nLimits, nSlIJK, nDeltaIndex)) {
|
|
for ( int nInd = 0 ; nInd < 3 ; ++ nInd)
|
|
nSlIJK[nInd] += nDeltaIndex[nInd] ;
|
|
int nAdBlockIJK[3] ;
|
|
int nAdBlockNum ;
|
|
if ( GetVoxelBlockIJK( nSlIJK, nAdBlockIJK) &&
|
|
GetBlockNFromIJK( nAdBlockIJK, nAdBlockNum) &&
|
|
m_BlockToUpdate[nAdBlockNum]) {
|
|
it = m_SliceYZ[t].erase( it) ;
|
|
continue ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
++ it ;
|
|
}
|
|
}
|
|
|
|
// Calcolo i triangoli sui blocchi
|
|
for ( int t = 0 ; t < m_nNumBlock ; ++ t) {
|
|
// Se il blocco deve essere processato
|
|
if ( m_BlockToUpdate[t]) {
|
|
// processo ...
|
|
ExtMarchingCubes( int( t), vVoxContainerVec[t]) ;
|
|
}
|
|
}
|
|
|
|
// Regolarizzo la catena
|
|
RegulateFeaturesChain( vVoxContainerVec) ;
|
|
|
|
// Costruisco i triangoli di feature
|
|
bool bCalcInterBlock = false ;
|
|
for ( int t = 0 ; t < m_nNumBlock ; ++ t) {
|
|
// Se il blocco è stato processato
|
|
if ( m_BlockToUpdate[t]) {
|
|
CreateSharpFeatureTriangle( t, vVoxContainerVec[t]) ;
|
|
// Flipping fra voxel interni
|
|
FlipEdgesII( t) ;
|
|
bCalcInterBlock = true ;
|
|
m_BlockToUpdate[t] = false ;
|
|
++ m_BlockUpdatingCounter[t] ;
|
|
// Sistemo le normali ai vertici (ciclo sui voxel)
|
|
for ( int t1 = 0 ; t1 < int( m_BlockSharpTria[t].size()) ; ++ t1) {
|
|
// ciclo sulle componenti connesse del voxel
|
|
for ( int t2 = 0 ; t2 < int( m_BlockSharpTria[t][t1].vCompoTria.size()) ; ++ t2) {
|
|
// ciclo sui triangoli delle componenti connesse
|
|
for ( int t3 = 0 ; t3 < int( m_BlockSharpTria[t][t1].vCompoTria[t2].size()) ; ++ t3) {
|
|
// Controllo normali
|
|
Vector3d vtN = m_BlockSharpTria[t][t1].vCompoTria[t2][t3].GetN() ;
|
|
bool bNormN = vtN.IsNormalized() ;
|
|
for ( int nV = 0 ; nV < 3 ; ++ nV) {
|
|
Vector3d vtNV = m_BlockSharpTria[t][t1].vCompoTria[t2][t3].GetVertexNorm( nV) ;
|
|
bool bNormV = vtNV.IsNormalized() ;
|
|
if ( bNormN && bNormV && vtN * vtNV < 0.7)
|
|
m_BlockSharpTria[t][t1].vCompoTria[t2][t3].SetVertexNorm( nV, vtN) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Calcolo i triangoli di frontiera tra feature di blocchi diversi (lavoro su copia e conservo gli originali)
|
|
SharpTriaMatrix InterBlockTria ;
|
|
if ( bCalcInterBlock) {
|
|
// Eseguo
|
|
m_InterBlockSharpTria = m_InterBlockOriginalSharpTria ;
|
|
FlipEdgesBB() ;
|
|
++ m_BlockUpdatingCounter[m_nNumBlock] ;
|
|
// Sistemo le normali ai vertici
|
|
for ( int t = 0 ; t < int( m_InterBlockSharpTria.size()) ; ++ t) {
|
|
for ( int t1 = 0 ; t1 < int( m_InterBlockSharpTria[t].size()) ; ++ t1) {
|
|
for ( int t2 = 0 ; t2 < int( m_InterBlockSharpTria[t][t1].vCompoTria.size()) ; ++ t2) {
|
|
for ( int t3 = 0 ; t3 < int( m_InterBlockSharpTria[t][t1].vCompoTria[t2].size()) ; ++ t3) {
|
|
if ( m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3].GetArea() > SQ_EPS_SMALL) {
|
|
// Controllo normali
|
|
Vector3d vtN = m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3].GetN() ;
|
|
bool bNormN = vtN.IsNormalized() ;
|
|
for ( int nV = 0 ; nV < 3 ; ++ nV) {
|
|
Vector3d vtNV = m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3].GetVertexNorm( nV) ;
|
|
bool bNormV = vtNV.IsNormalized() ;
|
|
if ( bNormN && bNormV && vtN * vtNV < 0.7)
|
|
m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3].SetVertexNorm( nV, vtN) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::ExtMarchingCubes( int nBlock, VoxelContainer& vVox) const
|
|
{
|
|
// Controllo sulla validità del blocco
|
|
if ( nBlock < 0 || nBlock >= int( m_nNumBlock))
|
|
return false ;
|
|
|
|
// Calcolo i limiti sugli indici dei voxel del blocco
|
|
// Vettore indici i,j,k del blocco
|
|
int nIJK[3] ;
|
|
GetBlockIJKFromN( nBlock, nIJK) ;
|
|
// Vettore limiti sugli indici dei voxel del blocco
|
|
int nLimits[6] ;
|
|
GetBlockLimitsIJK( nIJK, nLimits) ;
|
|
|
|
// Pulisco il contenitore dei voxel di frontiera
|
|
m_InterBlockVox[nBlock].clear() ;
|
|
m_InterBlockOriginalSharpTria[nBlock].clear() ;
|
|
m_BlockSharpTria[nBlock].clear() ;
|
|
m_BlockSmoothTria[nBlock].clear() ;
|
|
m_BlockBigTria[nBlock].clear() ;
|
|
|
|
// Unordered Map per la riduzione del numero di triangoli
|
|
int nDim = m_nVoxNumPerBlock * m_nVoxNumPerBlock ;
|
|
FlatVoxelContainer VoxContXZInf( nDim) ;
|
|
FlatVoxelContainer VoxContXZSup( nDim) ;
|
|
FlatVoxelContainer VoxContXYInf( nDim) ;
|
|
FlatVoxelContainer VoxContXYSup( nDim) ;
|
|
FlatVoxelContainer VoxContYZInf( nDim) ;
|
|
FlatVoxelContainer VoxContYZSup( nDim) ;
|
|
|
|
// Unordered map per la coerenza topologica nel blocco
|
|
InterVoxMatter SliceXY( 200) ;
|
|
InterVoxMatter SliceXZ( 200) ;
|
|
InterVoxMatter SliceYZ( 200) ;
|
|
|
|
// Ciclo su tutti i voxel del blocco
|
|
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) {
|
|
|
|
if ( m_nShape == BOX && ! IsVoxelOnBoxEdge( i, j, k))
|
|
continue ;
|
|
// 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.
|
|
AppliedVector 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]))
|
|
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.
|
|
AppliedVector CompoVert[4][7] ;
|
|
AppliedVector CompoTriVert[4][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[4] ;
|
|
|
|
int nExtTabOff = nComponents ;
|
|
int nStdTabOff = 0 ;
|
|
|
|
// Carico le matrici CompoVert e CompoTriVert
|
|
for ( int nComp = 0 ; nComp < nComponents ; ++ nComp) {
|
|
// Numero vertici per componenti
|
|
nVertComp[nComp] = TriangleTableEn[nIndex][1][nComp+1] ;
|
|
|
|
// Riempio il nCompCount-esimo vettore di vertici della base del fan
|
|
for ( int nVert = 0 ; nVert < nVertComp[nComp] ; ++ nVert)
|
|
CompoVert[nComp][nVert] = VecField[TriangleTableEn[nIndex][1][nVert + 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 nVert = 0 ; nVert < 3 * ( nVertComp[nComp] - 2) ; nVert += 3) {
|
|
CompoTriVert[nComp][nVert] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVert + 2]] ;
|
|
CompoTriVert[nComp][nVert+1] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVert + 1]] ;
|
|
CompoTriVert[nComp][nVert+2] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVert]] ;
|
|
}
|
|
|
|
// Aggiorno gli offsets per raggiungere i vertici della componente successiva.
|
|
nExtTabOff += nVertComp[nComp] ;
|
|
nStdTabOff += 3 * ( nVertComp[nComp] - 2) ;
|
|
}
|
|
|
|
// Controllo se il voxel ha una sola faccia che giace in un piano canonico e quindi ha gestione speciale
|
|
if ( m_nShape != BOX) {
|
|
// Faccia XY normale Z+
|
|
if ( nIndex == 15) {
|
|
int nTool ; double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], Z_PLUS, dPos, nTool)) {
|
|
int nN ; GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContXYSup.emplace( nN, HeigthAndColor( nTool, dPos)) ;
|
|
continue ;
|
|
}
|
|
}
|
|
// Faccia YZ normale X+
|
|
else if ( nIndex == 153) {
|
|
int nTool ; double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], X_PLUS, dPos, nTool)) {
|
|
int nN ; GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContYZSup.emplace( nN, HeigthAndColor( nTool, dPos)) ;
|
|
continue ;
|
|
}
|
|
}
|
|
// Faccia ZX normale Y+
|
|
else if ( nIndex == 51) {
|
|
int nTool ; double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], Y_PLUS, dPos, nTool)) {
|
|
int nN ; GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContXZSup.emplace( nN, HeigthAndColor( nTool, dPos)) ;
|
|
continue ;
|
|
}
|
|
}
|
|
// Faccia YX normale Z-
|
|
else if ( nIndex == 240) {
|
|
int nTool ; double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], Z_MINUS, dPos, nTool)) {
|
|
int nN ; GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContXYInf.emplace( nN, HeigthAndColor( nTool, dPos)) ;
|
|
continue ;
|
|
}
|
|
}
|
|
// Faccia ZY normale X-
|
|
else if ( nIndex == 102) {
|
|
int nTool ; double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], X_MINUS, dPos, nTool)) {
|
|
int nN ; GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContYZInf.emplace( nN, HeigthAndColor( nTool, dPos)) ;
|
|
continue ;
|
|
}
|
|
}
|
|
// Faccia XZ normale Y-
|
|
else if ( nIndex == 204) {
|
|
int nTool ; double dPos ;
|
|
if ( CanonicPlaneTest( CompoVert[0], Y_MINUS, dPos, nTool)) {
|
|
int nN ; GetVoxNFromIJK( i, j, k, nN) ;
|
|
VoxContXZInf.emplace( nN, HeigthAndColor( nTool, dPos)) ;
|
|
continue ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Configurazione 3
|
|
if ( nAllConfig[nIndex] == 3) {
|
|
// Test sulla topologia
|
|
bool bDefTopology = false ;
|
|
bool bMatOnSlice = false ;
|
|
int nCount = 0 ;
|
|
while ( nIndexConfig3[nCount] != nIndex)
|
|
++ nCount ;
|
|
// Vedo se la topologia è definita: se sì uso l'informazione
|
|
// passata dall'altro voxel, altrimenti la calcolo
|
|
int nIJKSl[3] = { ( nAdjVox3[nCount] != 1 ? i : i + 1),
|
|
( nAdjVox3[nCount] != 2 ? j : j + 1),
|
|
( nAdjVox3[nCount] != 3 ? k : k + 1) } ;
|
|
int nSliceN ;
|
|
int nSlBlockN ;
|
|
if ( GetVoxNFromIJK(nIJKSl[0], nIJKSl[1], nIJKSl[2], nSliceN)) {
|
|
int nSlBlockIJK[3] ;
|
|
GetVoxelBlockIJK(nIJKSl, nSlBlockIJK) ;
|
|
if ( abs( nAdjVox3[nCount]) == 1) {
|
|
auto it = SliceYZ.find( nSliceN) ;
|
|
if ( it != SliceYZ.end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
|
|
auto it = m_SliceYZ[nSlBlockN].find( nSliceN) ;
|
|
if ( it != m_SliceYZ[nSlBlockN].end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
}
|
|
}
|
|
else if ( abs( nAdjVox3[nCount]) == 2) {
|
|
auto it = SliceXZ.find( nSliceN) ;
|
|
if ( it != SliceXZ.end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
|
|
auto it = m_SliceXZ[nSlBlockN].find( nSliceN) ;
|
|
if ( it != m_SliceXZ[nSlBlockN].end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
}
|
|
}
|
|
else if ( abs( nAdjVox3[nCount]) == 3) {
|
|
auto it = SliceXY.find( nSliceN) ;
|
|
if ( it != SliceXY.end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
|
|
auto it = m_SliceXY[nSlBlockN].find( nSliceN) ;
|
|
if ( it != m_SliceXY[nSlBlockN].end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
// La topologia è indefinita: la calcolo
|
|
if ( ! bDefTopology && bReg) {
|
|
|
|
Point3d ptFirstBar = ( CompoVert[0][0].ptPApp + CompoVert[0][1].ptPApp + CompoVert[0][2].ptPApp) / 3 ;
|
|
Point3d ptSecondBar = ( CompoVert[1][0].ptPApp + CompoVert[1][1].ptPApp + CompoVert[1][2].ptPApp) / 3 ;
|
|
Vector3d vtDiagBar = ptSecondBar - ptFirstBar ;
|
|
vtDiagBar.Normalize() ;
|
|
int nSum = 0 ;
|
|
for ( int nVec = 0 ; nVec < 3 ; ++ nVec) {
|
|
if ( CompoVert[0][nVec].vtVec * vtDiagBar > 0.7)
|
|
++ nSum ;
|
|
else if ( CompoVert[0][nVec].vtVec * vtDiagBar < -0.7)
|
|
-- nSum;
|
|
if ( CompoVert[1][nVec].vtVec * vtDiagBar < -0.7)
|
|
++ nSum ;
|
|
else if ( CompoVert[1][nVec].vtVec * vtDiagBar > 0.7)
|
|
-- nSum ;
|
|
}
|
|
|
|
if ( nSum < - 2)
|
|
bMatOnSlice = true ;
|
|
if ( nSum <= 2) {
|
|
Vector3d vtVoxCentre ;
|
|
if ( nAdjVox3[nCount] == - 1)
|
|
vtVoxCentre = X_AX ;
|
|
else if ( nAdjVox3[nCount] == 1)
|
|
vtVoxCentre = - X_AX ;
|
|
else if ( nAdjVox3[nCount] == - 2)
|
|
vtVoxCentre = Y_AX ;
|
|
else if ( nAdjVox3[nCount] == 2)
|
|
vtVoxCentre = - Y_AX ;
|
|
else if ( nAdjVox3[nCount] == - 3)
|
|
vtVoxCentre = Z_AX ;
|
|
else if ( nAdjVox3[nCount] == 3)
|
|
vtVoxCentre = - Z_AX ;
|
|
|
|
int nPlusNum = 0 ;
|
|
int nMinusNum = 0 ;
|
|
for ( int nVec = 0 ; nVec < 3 ; ++ nVec) {
|
|
if ( CompoVert[0][nVec].vtVec * vtVoxCentre > EPS_SMALL)
|
|
++ nPlusNum ;
|
|
else if ( CompoVert[0][nVec].vtVec * vtVoxCentre < - EPS_SMALL)
|
|
-- nMinusNum ;
|
|
if ( CompoVert[1][nVec].vtVec * vtVoxCentre > EPS_SMALL)
|
|
++ nPlusNum ;
|
|
else if ( CompoVert[1][nVec].vtVec * vtVoxCentre < - EPS_SMALL)
|
|
-- nMinusNum ;
|
|
}
|
|
|
|
if ( nPlusNum >= nMinusNum)
|
|
bMatOnSlice = true ;
|
|
}
|
|
else
|
|
;
|
|
}
|
|
// Conservo l'informazione per i voxel successivi
|
|
if ( GetVoxNFromIJK( nIJKSl[0], nIJKSl[1], nIJKSl[2], nSliceN)) {
|
|
if ( abs(nAdjVox3[nCount]) == 1) {
|
|
if ( nSlBlockN == nBlock)
|
|
SliceYZ.emplace( nSliceN, bMatOnSlice) ;
|
|
else
|
|
m_SliceYZ[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
|
|
}
|
|
else if ( abs(nAdjVox3[nCount]) == 2) {
|
|
if ( nSlBlockN == nBlock)
|
|
SliceXZ.emplace( nSliceN, bMatOnSlice) ;
|
|
else
|
|
m_SliceXZ[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
|
|
}
|
|
else if ( abs(nAdjVox3[nCount]) == 3) {
|
|
if ( nSlBlockN == nBlock)
|
|
SliceXY.emplace(nSliceN, bMatOnSlice) ;
|
|
else
|
|
m_SliceXY[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
|
|
}
|
|
}
|
|
|
|
bool bNewTopology = Config3Duality( nIndex) ? bMatOnSlice : ! bMatOnSlice ;
|
|
|
|
// Si passa alla seconda topologia
|
|
if ( bNewTopology) {
|
|
// Ricerca del caso corrispondente della nuova topologia
|
|
int nt = 0 ;
|
|
while ( nIndexVsIndex3[nt][0] != nIndex)
|
|
++ nt;
|
|
int nRotCase = nIndexVsIndex3[nt][1] ;
|
|
// Aggiorno numero di componenti
|
|
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) ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Configurazione 6
|
|
else if ( nAllConfig[nIndex] == 6) {
|
|
// Test sulla topologia
|
|
bool bDefTopology = false ;
|
|
bool bMatOnSlice = false ;
|
|
int nCount = 0 ;
|
|
while ( nIndexConfig6[nCount] != nIndex)
|
|
++ nCount ;
|
|
// Vedo se la topologia è definita: altrimenti la ricalcolo
|
|
int nIJKSl[3] = { ( nAdjVox6[nCount] != 1 ? i : i + 1),
|
|
( nAdjVox6[nCount] != 2 ? j : j + 1),
|
|
( nAdjVox6[nCount] != 3 ? k : k + 1)} ;
|
|
int nSliceN ;
|
|
int nSlBlockN ;
|
|
if ( GetVoxNFromIJK( nIJKSl[0], nIJKSl[1], nIJKSl[2], nSliceN)) {
|
|
int nSlBlockIJK[3] ;
|
|
GetVoxelBlockIJK( nIJKSl, nSlBlockIJK) ;
|
|
if ( abs( nAdjVox6[nCount]) == 1) {
|
|
auto it = SliceYZ.find( nSliceN) ;
|
|
if ( it != SliceYZ.end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
|
|
auto it = m_SliceYZ[nSlBlockN].find( nSliceN) ;
|
|
if ( it != m_SliceYZ[nSlBlockN].end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
}
|
|
}
|
|
else if ( abs( nAdjVox6[nCount]) == 2) {
|
|
auto it = SliceXZ.find( nSliceN) ;
|
|
if ( it != SliceXZ.end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
|
|
auto it = m_SliceXZ[nSlBlockN].find( nSliceN) ;
|
|
if ( it != m_SliceXZ[nSlBlockN].end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
}
|
|
}
|
|
else if ( abs( nAdjVox6[nCount]) == 3) {
|
|
auto it = SliceXY.find( nSliceN) ;
|
|
if ( it != SliceXY.end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
|
|
auto it = m_SliceXY[nSlBlockN].find( nSliceN) ;
|
|
if ( it != m_SliceXY[nSlBlockN].end()) {
|
|
bMatOnSlice = it->second ;
|
|
bDefTopology = true ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
// Topologia indefinita: la calcolo
|
|
if ( ! bDefTopology && bReg) {
|
|
// Test sulla topologia
|
|
Point3d ptFirstBar = ( CompoVert[0][0].ptPApp + CompoVert[0][1].ptPApp +
|
|
CompoVert[0][2].ptPApp + CompoVert[0][3].ptPApp) / 4 ;
|
|
Point3d ptSecondBar = ( CompoVert[1][0].ptPApp + CompoVert[1][1].ptPApp +
|
|
CompoVert[1][2].ptPApp) / 3 ;
|
|
Vector3d vtDiagBar = ptSecondBar - ptFirstBar ;
|
|
vtDiagBar.Normalize() ;
|
|
int nSum = 0 ;
|
|
for ( int nVec = 0 ; nVec < 4 ; ++ nVec) {
|
|
if ( CompoVert[0][nVec].vtVec * vtDiagBar > 0.7)
|
|
++ nSum ;
|
|
else if ( CompoVert[0][nVec].vtVec * vtDiagBar < - 0.7)
|
|
-- nSum ;
|
|
}
|
|
for ( int nVec = 0 ; nVec < 3 ; ++ nVec) {
|
|
if ( CompoVert[1][nVec].vtVec * vtDiagBar < - 0.7)
|
|
++ nSum ;
|
|
else if ( CompoVert[1][nVec].vtVec * vtDiagBar > 0.7)
|
|
-- nSum ;
|
|
}
|
|
|
|
if ( nSum < - 3)
|
|
bMatOnSlice = true ;
|
|
else if ( nSum <= 3) {
|
|
Vector3d vtVoxCentre ;
|
|
if ( nAdjVox6[nCount] == -1)
|
|
vtVoxCentre = X_AX ;
|
|
else if ( nAdjVox6[nCount] == 1)
|
|
vtVoxCentre = - X_AX ;
|
|
else if ( nAdjVox6[nCount] == -2)
|
|
vtVoxCentre = Y_AX ;
|
|
else if ( nAdjVox6[nCount] == 2)
|
|
vtVoxCentre = - Y_AX ;
|
|
else if ( nAdjVox6[nCount] == -3)
|
|
vtVoxCentre = Z_AX ;
|
|
else if ( nAdjVox6[nCount] == 3)
|
|
vtVoxCentre = - Z_AX ;
|
|
|
|
int nPlusNum = 0 ;
|
|
int nMinusNum = 0 ;
|
|
for ( int nVec = 0 ; nVec < 4 ; ++ nVec) {
|
|
if ( CompoVert[0][nVec].vtVec * vtVoxCentre > EPS_SMALL)
|
|
++ nPlusNum ;
|
|
else if ( CompoVert[0][nVec].vtVec * vtVoxCentre < -EPS_SMALL)
|
|
-- nMinusNum ;
|
|
}
|
|
for ( int nVec = 0 ; nVec < 3 ; ++ nVec) {
|
|
if ( CompoVert[1][nVec].vtVec * vtVoxCentre > EPS_SMALL)
|
|
++ nPlusNum ;
|
|
else if ( CompoVert[1][nVec].vtVec * vtVoxCentre < -EPS_SMALL)
|
|
-- nMinusNum ;
|
|
}
|
|
|
|
if ( nPlusNum >= nMinusNum)
|
|
bMatOnSlice = true ;
|
|
|
|
}
|
|
else
|
|
;
|
|
}
|
|
// Conservo l'informazione
|
|
if ( GetVoxNFromIJK( nIJKSl[0], nIJKSl[1], nIJKSl[2], nSliceN)) {
|
|
if ( abs(nAdjVox6[nCount]) == 1) {
|
|
if ( nSlBlockN == nBlock)
|
|
SliceYZ.emplace( nSliceN, bMatOnSlice) ;
|
|
else
|
|
m_SliceYZ[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
|
|
}
|
|
else if ( abs(nAdjVox6[nCount]) == 2) {
|
|
if ( nSlBlockN == nBlock)
|
|
SliceXZ.emplace( nSliceN, bMatOnSlice) ;
|
|
else
|
|
m_SliceXZ[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
|
|
}
|
|
else if (abs(nAdjVox6[nCount]) == 3) {
|
|
if ( nSlBlockN == nBlock)
|
|
SliceXY.emplace( nSliceN, bMatOnSlice) ;
|
|
else
|
|
m_SliceXY[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
|
|
}
|
|
}
|
|
bool bNewTopology = ( Config6Duality( nIndex) ? bMatOnSlice : ! bMatOnSlice) ;
|
|
// Si deve passare alla seconda topologia
|
|
if ( bNewTopology) {
|
|
// Ricerca del caso corrispondente della nuova topologia
|
|
int nt = 0 ;
|
|
while ( nIndexVsIndex6[nt][0] != nIndex)
|
|
++ nt ;
|
|
int nRotCase = nIndexVsIndex6[nt][1] ;
|
|
// Aggiorno numero di componenti
|
|
nComponents = Cases6Plus[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] = Cases6Plus[nRotCase][1][nC] ;
|
|
// Matrice dei vertici della base del fan
|
|
for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert)
|
|
CompoVert[nC - 1][nFanVert] = VecField[Cases6Plus[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[Cases6Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
|
|
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases6Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
|
|
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases6Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
|
|
}
|
|
// Aggiorno gli offsets per raggiungere i vertici della componente successiva.
|
|
nExtTabOff += nVertComp[nC - 1] ;
|
|
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Configurazione 7
|
|
else if ( nAllConfig[nIndex] == 7) {
|
|
// !!! Versione provvisoria, deve essere riveduta e semplificata !!!
|
|
// Test sulla topologia
|
|
bool bMatOnSliceYZ = false ;
|
|
bool bMatOnSliceXZ = false ;
|
|
bool bMatOnSliceXY = false ;
|
|
bool bDefSliceYZ = false ;
|
|
bool bDefSliceXZ = false ;
|
|
bool bDefSliceXY = false ;
|
|
int nCount = 0 ;
|
|
while ( nIndexConfig7[nCount] != nIndex)
|
|
++ nCount ;
|
|
|
|
// Vedo se la topologia è definita: se sì uso l'informazione già posseduta,
|
|
// altrimenti devo calcolare la topologia
|
|
|
|
// Numerazione delle facce del voxel: 0: YZ- 1: XZ- 2: YZ+ 3: XZ+ 4: XY- 5: XY+
|
|
// - 1 faccia non in gioco, 0 faccia in gioco ma con topologia non definita,
|
|
// 1 faccia in gioco con topologia definita
|
|
int nFace[6] = { -1, -1, -1, -1, -1, -1} ;
|
|
// Facce in gioco
|
|
int nCurFaceYZ = 0 ;
|
|
int nCurFaceXZ = 1 ;
|
|
int nCurFaceXY = 4 ;
|
|
|
|
int nIJKSlYZ[3] = { i, j, k} ;
|
|
if ( nAdjVox7[nCount][0] == 1) {
|
|
++ nIJKSlYZ[0] ;
|
|
nCurFaceYZ = 2 ;
|
|
}
|
|
int nIJKSlXZ[3] = { i, j, k} ;
|
|
if ( nAdjVox7[nCount][1] == 2) {
|
|
++ nIJKSlXZ[1] ;
|
|
nCurFaceXZ = 3 ;
|
|
}
|
|
int nIJKSlXY[3] = { i, j, k} ;
|
|
if ( nAdjVox7[nCount][2] == 3) {
|
|
++ nIJKSlXY[2] ;
|
|
nCurFaceXY = 5 ;
|
|
}
|
|
// Assegno l'indice zero alle facce in gioco. L'indice zero è per le facce in gioco
|
|
// su cui non si è certi se ci sia o meno materiale.
|
|
nFace[nCurFaceYZ] = 0 ;
|
|
nFace[nCurFaceXZ] = 0 ;
|
|
nFace[nCurFaceXY] = 0 ;
|
|
|
|
// Per ogni faccia in gioco vediamo se c'è materiale: se sì assegnamo l'indice 1
|
|
int nSlYZN ;
|
|
int nSlYZBlockN ;
|
|
if ( GetVoxNFromIJK( nIJKSlYZ[0], nIJKSlYZ[1], nIJKSlYZ[2], nSlYZN)) {
|
|
// Slice interna al blocco
|
|
auto it = SliceYZ.find( nSlYZN) ;
|
|
// Topologia definita
|
|
if ( it != SliceYZ.end()) {
|
|
bMatOnSliceYZ = it->second ;
|
|
bDefSliceYZ = true ;
|
|
nFace[nCurFaceYZ] = 1 ;
|
|
}
|
|
// Slice sulla frontiera
|
|
int nSlBlockIJK[3] ;
|
|
GetVoxelBlockIJK( nIJKSlYZ, nSlBlockIJK) ;
|
|
if ( GetBlockNFromIJK( nSlBlockIJK, nSlYZBlockN)) {
|
|
auto it = m_SliceYZ[nSlYZBlockN].find( nSlYZN) ;
|
|
// Topologia definita
|
|
if ( it != m_SliceYZ[nSlYZBlockN].end()) {
|
|
bMatOnSliceYZ = it->second ;
|
|
bDefSliceYZ = true ;
|
|
nFace[nCurFaceYZ] = 1 ;
|
|
}
|
|
}
|
|
}
|
|
int nSlXZN ;
|
|
int nSlXZBlockN ;
|
|
if ( GetVoxNFromIJK( nIJKSlXZ[0], nIJKSlXZ[1], nIJKSlXZ[2], nSlXZN)) {
|
|
// Slice interna al blocco
|
|
auto it = SliceXZ.find( nSlXZN) ;
|
|
// Topologia definita
|
|
if ( it != SliceXZ.end()) {
|
|
bMatOnSliceXZ = it->second ;
|
|
bDefSliceXZ = true ;
|
|
nFace[nCurFaceXZ] = 1 ;
|
|
}
|
|
// Slice sulla frontiera
|
|
int nSlBlockIJK[3] ;
|
|
GetVoxelBlockIJK( nIJKSlXZ, nSlBlockIJK) ;
|
|
if ( GetBlockNFromIJK( nSlBlockIJK, nSlXZBlockN)) {
|
|
auto it = m_SliceYZ[nSlXZBlockN].find( nSlXZN) ;
|
|
// Topologia definita
|
|
if ( it != m_SliceYZ[nSlXZBlockN].end()) {
|
|
bMatOnSliceXZ = it->second ;
|
|
bDefSliceXZ = true ;
|
|
nFace[nCurFaceXZ] = 1 ;
|
|
}
|
|
}
|
|
}
|
|
int nSlXYN ;
|
|
int nSlXYBlockN ;
|
|
if ( GetVoxNFromIJK( nIJKSlXY[0], nIJKSlXY[1], nIJKSlXY[2], nSlXYN)) {
|
|
// Slice interna al blocco
|
|
auto it = SliceXY.find( nSlXYN) ;
|
|
// Topologia definita
|
|
if ( it != SliceXY.end()) {
|
|
bMatOnSliceXY = it->second ;
|
|
bDefSliceXY = true ;
|
|
nFace[nCurFaceXY] = 1 ;
|
|
}
|
|
// Slice sulla frontiera
|
|
int nSlBlockIJK[3] ;
|
|
GetVoxelBlockIJK( nIJKSlXY, nSlBlockIJK) ;
|
|
if ( GetBlockNFromIJK( nSlBlockIJK, nSlXYBlockN)) {
|
|
auto it = m_SliceYZ[nSlXYBlockN].find( nSlXYN) ;
|
|
// Topologia definita
|
|
if ( it != m_SliceYZ[nSlXYBlockN].end()) {
|
|
bMatOnSliceXY = it->second ;
|
|
bDefSliceXY = true ;
|
|
nFace[nCurFaceXY] = 1 ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Numerazione delle facce del voxel 0: YZ- 1: XZ- 2: YZ+ 3: XZ+ 4: XY- 5: XY+
|
|
// Gli spigoli sono ordinati in senso antiorario dal punto di vista di un osservatore esterno al voxel
|
|
static int nSliceEdges[6][4] = { { 3, 8, 7, 11}, { 0, 9, 4, 8}, { 1, 10, 5, 9}, { 2, 11, 6, 10},
|
|
{ 0, 3, 2, 1}, { 4, 5, 6, 7} } ;
|
|
// Le facce in gioco hano indici nCurFaceXY, nCurFaceXZ e nCurFaceYZ. Se per almeno una faccia non si sa se c'è materiale,
|
|
// ovvero almeno una faccia ha indice 0, e il campo vettoriale normale è regolare verifico la presenza di materiale.
|
|
if ( ( nFace[nCurFaceXY] == 0 || nFace[nCurFaceXZ] == 0 || nFace[nCurFaceYZ] == 0) && bReg) {
|
|
// Ciclo sulle facce
|
|
for ( int nFaceN = 0 ; nFaceN < 6 ; ++ nFaceN) {
|
|
// Faccia da analizzare
|
|
if ( nFace[nFaceN] == 0) {
|
|
// Dalla tabella determino le due componenti connesse che si appoggiano alla faccia
|
|
int nFaceCompo1 = 0 ;
|
|
int nFaceCompo2 = 0 ;
|
|
// Ciclo sulle componenti connesse del voxel
|
|
for ( int nCurComp = 1 ; nCurComp <= 3 ; ++ nCurComp) {
|
|
int nMatchEdge = 0 ;
|
|
// Ciclo sugli edge della componente
|
|
for ( int nEdge = 0 ; nEdge < 3 ; ++ nEdge) {
|
|
// Ciclo sugli edge della faccia
|
|
for ( int nFaceEdge = 0 ; nFaceEdge < 4 ; ++ nFaceEdge) {
|
|
// Edge della componente coincide con quello della faccia
|
|
if ( nSliceEdges[nFaceN][nFaceEdge] == TriangleTableEn[nIndex][1][3 + nCurComp - 1 + nEdge]) {
|
|
++ nMatchEdge ;
|
|
break ;
|
|
}
|
|
}
|
|
// La componente corrente ha due edge sulla faccia,
|
|
// interrompiamo la ricerca di edge corrispondenti
|
|
if ( nMatchEdge == 2)
|
|
break ;
|
|
}
|
|
// Abbiamo trovato una nuova componente con due edge sulla faccia
|
|
if ( nMatchEdge == 2) {
|
|
if ( nFaceCompo1 == 0)
|
|
nFaceCompo1 = nCurComp ;
|
|
else
|
|
nFaceCompo2 = nCurComp ;
|
|
}
|
|
// Se le componenti con due edge sulla faccia sono due, interrompiamo la ricerca
|
|
if ( nFaceCompo2 != 0)
|
|
break ;
|
|
}
|
|
// Valuto la topologia VecField[EdgeIndex] edgeIndex = 0, 1, ..., 11
|
|
double dDotSum = 0 ;
|
|
for ( int nV = 0 ; nV < 3 ; ++ nV) {
|
|
Vector3d vtV1 = VecField[TriangleTableEn[nIndex][1][3 + nFaceCompo1 - 1 + nV]].vtVec ;
|
|
Vector3d vtV2 = VecField[TriangleTableEn[nIndex][1][3 + nFaceCompo2 - 1 + nV]].vtVec ;
|
|
dDotSum += ( vtV1 * vtV2) ;
|
|
}
|
|
if ( nFace[nFaceN] == 0 || nFace[nFaceN] == 2)
|
|
bMatOnSliceXZ = dDotSum > 0. ;
|
|
else if ( nFace[nFaceN] == 1 || nFace[nFaceN] == 3)
|
|
bMatOnSliceYZ = dDotSum > 0. ;
|
|
else
|
|
bMatOnSliceXY = dDotSum > 0. ;
|
|
}
|
|
}
|
|
}
|
|
// Uso le informazioni per scegliere la topologia
|
|
int nFaceWithMatNum = 0 ;
|
|
if ( bMatOnSliceXZ)
|
|
++ nFaceWithMatNum ;
|
|
if ( bMatOnSliceYZ)
|
|
++ nFaceWithMatNum ;
|
|
if ( bMatOnSliceXY)
|
|
++ nFaceWithMatNum ;
|
|
if ( nFaceWithMatNum == 1) {
|
|
int nFaceCase = ( bMatOnSliceYZ ? 0 : (bMatOnSliceXZ ? 1 : 2)) ;
|
|
// Aggiorno numero di componenti
|
|
nComponents = Cases7Plus[nCount][nFaceCase][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] = Cases7Plus[nCount][nFaceCase][1][nC] ;
|
|
// Matrice dei vertici della base del fan
|
|
for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert)
|
|
CompoVert[nC - 1][nFanVert] = VecField[Cases7Plus[nCount][nFaceCase][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[Cases7Plus[nCount][nFaceCase][0][nStdTabOff + nTriVert + 2]] ;
|
|
CompoTriVert[nC - 1][nTriVert + 1] = VecField[Cases7Plus[nCount][nFaceCase][0][nStdTabOff + nTriVert + 1]] ;
|
|
CompoTriVert[nC - 1][nTriVert + 2] = VecField[Cases7Plus[nCount][nFaceCase][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 ( nFaceWithMatNum == 2) {
|
|
;
|
|
}
|
|
// Conservo l'informazione
|
|
if ( GetVoxNFromIJK( nIJKSlYZ[0], nIJKSlYZ[1], nIJKSlYZ[2], nSlYZN)) {
|
|
if ( nSlYZBlockN == nBlock)
|
|
SliceYZ.emplace( nSlYZN, bMatOnSliceYZ) ;
|
|
else
|
|
m_SliceYZ[nSlYZBlockN].emplace( nSlYZN, bMatOnSliceYZ) ;
|
|
}
|
|
if ( GetVoxNFromIJK( nIJKSlXZ[0], nIJKSlXZ[1], nIJKSlXZ[2], nSlXZN)) {
|
|
if ( nSlXZBlockN == nBlock)
|
|
SliceXZ.emplace( nSlXZN, bMatOnSliceXZ) ;
|
|
else
|
|
m_SliceXZ[nSlXZBlockN].emplace( nSlXZN, bMatOnSliceXZ) ;
|
|
}
|
|
if ( GetVoxNFromIJK( nIJKSlXY[0], nIJKSlXY[1], nIJKSlXY[2], nSlXYN)) {
|
|
if ( nSlXYBlockN == nBlock)
|
|
SliceXY.emplace( nSlXYN, bMatOnSliceXY) ;
|
|
else
|
|
m_SliceXY[nSlXYBlockN].emplace( nSlXYN, bMatOnSliceXY) ;
|
|
}
|
|
}
|
|
|
|
// Configurazione 10
|
|
else if ( nAllConfig[nIndex] == 10) {
|
|
// Test sulla topologia
|
|
bool bDefStTopology = false ;
|
|
bool bMatOnSlice = false ;
|
|
int nCount = 0 ;
|
|
while ( nIndexConfig10[nCount] != nIndex)
|
|
++ nCount ;
|
|
// Vedo se la topologia è definita: se sì uso l'informazione già posseduta,
|
|
// altrimenti devo calcolare la topologia
|
|
int nIJKSlSt[3] = { i, j, k} ;
|
|
int nIJKSlEn[3] = { ( nAdjVox10[nCount] != 1 ? i : i + 1),
|
|
( nAdjVox10[nCount] != 2 ? j : j + 1),
|
|
( nAdjVox10[nCount] != 3 ? k : k + 1)} ;
|
|
int nSliceStN, nSliceEnN ;
|
|
int nSlBlockEnN ;
|
|
if ( GetVoxNFromIJK( nIJKSlSt[0], nIJKSlSt[1], nIJKSlSt[2], nSliceStN) &&
|
|
GetVoxNFromIJK( nIJKSlEn[0], nIJKSlEn[1], nIJKSlEn[2], nSliceEnN)) {
|
|
if ( abs( nAdjVox10[nCount]) == 1) {
|
|
auto itSt = SliceYZ.find( nSliceStN) ;
|
|
if ( itSt != SliceYZ.end()) {
|
|
bMatOnSlice = itSt->second ;
|
|
bDefStTopology = true ;
|
|
}
|
|
}
|
|
else if ( abs( nAdjVox10[nCount]) == 2) {
|
|
auto itSt = SliceXZ.find( nSliceStN) ;
|
|
if ( itSt != SliceXZ.end()) {
|
|
bMatOnSlice = itSt->second ;
|
|
bDefStTopology = true ;
|
|
}
|
|
}
|
|
else if ( abs( nAdjVox10[nCount]) == 3) {
|
|
auto itSt = SliceXY.find( nSliceStN) ;
|
|
if ( itSt != SliceXY.end()) {
|
|
bMatOnSlice = itSt->second ;
|
|
bDefStTopology = true ;
|
|
}
|
|
}
|
|
}
|
|
// La topologia non è definita, la calcolo
|
|
if ( ! bDefStTopology && bReg) {
|
|
// Verifico concordanza tra i versori di una stessa componente
|
|
// (ogni coppia di vettori di una medesima componente deve avere prodotto scalare non inferiore a 0.0)
|
|
Vector3d vtCmpAvg0, vtCmpAvg1 ;
|
|
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.0) ;
|
|
bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1, 0.0) ;
|
|
bMatOnSlice = ( ! bTest0 || ! bTest1) ;
|
|
}
|
|
// Conservo l'informazioe e la trasmetto al voxel successivo
|
|
if ( GetVoxNFromIJK( nIJKSlSt[0], nIJKSlSt[1], nIJKSlSt[2], nSliceStN) &&
|
|
GetVoxNFromIJK( nIJKSlEn[0], nIJKSlEn[1], nIJKSlEn[2], nSliceEnN)) {
|
|
if ( abs( nAdjVox6[nCount]) == 1) {
|
|
if ( GetBlockNFromIJK( nIJKSlEn, nSlBlockEnN)) {
|
|
auto it = m_SliceYZ[nSlBlockEnN].find( nSliceEnN) ;
|
|
if ( it != m_SliceYZ[nSlBlockEnN].end()) {
|
|
if ( it->second != bMatOnSlice)
|
|
m_BlockToUpdate[nSlBlockEnN] = true ;
|
|
it->second = bMatOnSlice ;
|
|
}
|
|
else
|
|
m_SliceYZ[nSlBlockEnN].emplace( nSliceEnN, bMatOnSlice) ;
|
|
}
|
|
}
|
|
else if ( abs( nAdjVox6[nCount]) == 2) {
|
|
if ( GetBlockNFromIJK( nIJKSlEn, nSlBlockEnN)) {
|
|
auto it = m_SliceXZ[nSlBlockEnN].find( nSliceEnN) ;
|
|
if ( it != m_SliceXZ[nSlBlockEnN].end()) {
|
|
if ( it->second != bMatOnSlice)
|
|
m_BlockToUpdate[nSlBlockEnN] = true ;
|
|
it->second = bMatOnSlice ;
|
|
}
|
|
else
|
|
m_SliceXZ[nSlBlockEnN].emplace( nSliceEnN, bMatOnSlice) ;
|
|
}
|
|
}
|
|
else if ( abs( nAdjVox6[nCount]) == 3) {
|
|
if ( GetBlockNFromIJK( nIJKSlEn, nSlBlockEnN)) {
|
|
auto it = m_SliceXY[nSlBlockEnN].find( nSliceEnN) ;
|
|
if ( it != m_SliceXY[nSlBlockEnN].end()) {
|
|
if ( it->second != bMatOnSlice)
|
|
m_BlockToUpdate[nSlBlockEnN] = true ;
|
|
it->second = bMatOnSlice ;
|
|
}
|
|
else
|
|
m_SliceXY[nSlBlockEnN].emplace( nSliceEnN, bMatOnSlice) ;
|
|
}
|
|
}
|
|
}
|
|
bool bNewTopology = bMatOnSlice ;
|
|
// Si passa alla seconda topologia
|
|
if ( bNewTopology) {
|
|
// Ricerca del caso corrispondente della nuova topologia
|
|
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) ;
|
|
}
|
|
}
|
|
}
|
|
|
|
Voxel VoxConf ;
|
|
VoxConf.nNumComp = 0 ;
|
|
|
|
SmoothTriaStruct VoxSmoothTria ;
|
|
int nVoxSmootSizePrev = int( VoxSmoothTria.vTria.size()) ;
|
|
int nVoxSmootSize = nVoxSmootSizePrev ;
|
|
|
|
// Numero di feature nel voxel: al più vi è una feature per componente connessa.
|
|
int nFeatureInVoxel = 0 ;
|
|
|
|
// Ciclo sulle componenti
|
|
for ( int nComp = 0 ; nComp < nComponents ; ++ nComp) {
|
|
|
|
int nFeatureType = NO_FEATURE ;
|
|
// Se i componenti sono regolari valuto le normali per stabilire se eseguire ExtMC o MC
|
|
if ( bReg)
|
|
nFeatureType = TestOnNormal( CompoVert[nComp], nVertComp[nComp]) ;
|
|
|
|
// Extended MC
|
|
if ( nFeatureType != NO_FEATURE) {
|
|
|
|
// Passo al sistema di riferimento del baricentro
|
|
Point3d ptGravityCenter( 0, 0, 0) ;
|
|
for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni)
|
|
ptGravityCenter += CompoVert[nComp][ni].ptPApp ;
|
|
ptGravityCenter = ptGravityCenter / nVertComp[nComp] ;
|
|
|
|
Vector3d vtTrasf[12] ;
|
|
for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni)
|
|
vtTrasf[ni] = CompoVert[nComp][ni].ptPApp - ptGravityCenter ;
|
|
|
|
// Preparo le matrici per il sistema
|
|
SvdMatrix dMatrixN( nVertComp[nComp], 3) ;
|
|
SvdVector dKnownVector( nVertComp[nComp]) ;
|
|
|
|
// medio le normali adiacenti molto vicine (delta angolare inferiore a 22.5 deg)
|
|
Vector3d vtNorm[12] ;
|
|
for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni)
|
|
vtNorm[ni] = CompoVert[nComp][ni].vtVec ;
|
|
for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) {
|
|
int nj = ( ni + 1) % nVertComp[nComp] ;
|
|
if ( vtNorm[ni] * vtNorm[nj] > 0.92) {
|
|
Vector3d vtNI = ( 0.6 * vtNorm[ni] + 0.4 * vtNorm[nj]) ;
|
|
Vector3d vtNJ = ( 0.4 * vtNorm[ni] + 0.6 * vtNorm[nj]) ;
|
|
vtNorm[ni] = vtNI ;
|
|
vtNorm[ni].Normalize() ;
|
|
vtNorm[nj] = vtNJ ;
|
|
vtNorm[nj].Normalize() ;
|
|
++ ni ;
|
|
}
|
|
}
|
|
|
|
// Definisco la matrice del sistema
|
|
for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) {
|
|
dMatrixN( ni, 0) = vtNorm[ni].x ;
|
|
dMatrixN( ni, 1) = vtNorm[ni].y ;
|
|
dMatrixN( ni, 2) = vtNorm[ni].z ;
|
|
dKnownVector( ni) = vtNorm[ni] * vtTrasf[ni] ;
|
|
}
|
|
|
|
// calcolo SVD
|
|
SvdDecomposer svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ;
|
|
auto dMatrixV = svd.matrixV() ;
|
|
auto 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
|
|
auto 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 ;
|
|
|
|
// Gestisco i ribaltamenti dei triangoli
|
|
bool bExtConfirmed = true ;
|
|
Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
|
|
if ( nFeatureType == EDGE && vtNullSpace.Normalize() && ! IsPointInsideVoxelApprox( i, j, k, ptSol)) {
|
|
int nNewAdjVoxIJK[3] ;
|
|
if ( GetPointVoxel( ptSol, nNewAdjVoxIJK[0], nNewAdjVoxIJK[1], nNewAdjVoxIJK[2])) {
|
|
int nNewAdjVoxIndex = CalcIndex( nNewAdjVoxIJK[0], nNewAdjVoxIJK[1], nNewAdjVoxIJK[2]) ;
|
|
if ( nNewAdjVoxIndex != 0) {
|
|
bool bOverTurning = false ;
|
|
for ( int m = 0 ; m < nVertComp[nComp] ; ++ m) {
|
|
int l = ( m + 1) % nVertComp[nComp] ;
|
|
Vector3d vtNlm = ( CompoVert[nComp][l].ptPApp - ptSol) ^ ( CompoVert[nComp][m].ptPApp - ptSol) ;
|
|
vtNlm.Normalize() ;
|
|
double dDotl = vtNlm * CompoVert[nComp][l].vtVec ;
|
|
double dDotm = vtNlm * CompoVert[nComp][m].vtVec ;
|
|
if ( dDotl < - 0.85 || dDotm < - 0.85) {
|
|
bOverTurning = true ;
|
|
break ;
|
|
}
|
|
}
|
|
if ( bOverTurning) {
|
|
Point3d ptMinDiag( ( i * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( j * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( k * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
Point3d ptMaxDiag( ( ( i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( k + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
double dNewU1, dNewU2 ;
|
|
if ( IntersLineBox( ptSol, vtNullSpace, ptMinDiag, ptMaxDiag, dNewU1, dNewU2)) {
|
|
if ( dNewU1 > 0.) {
|
|
ptSol += ( dNewU1 + EPS_SMALL * ( dNewU2 - dNewU1)) * vtNullSpace ;
|
|
}
|
|
else if ( dNewU2 < 0.) {
|
|
ptSol += ( dNewU2 + EPS_SMALL * ( dNewU1 - dNewU2)) * vtNullSpace ;
|
|
}
|
|
}
|
|
else {
|
|
|
|
if ( nVertComp[nComp] == 5) {
|
|
double dDotAvarage = 0.;
|
|
for (int m = 0; m < nVertComp[nComp] - 1; ++m) {
|
|
for (int l = m + 1 ; l < nVertComp[nComp]; ++l) {
|
|
dDotAvarage += CompoVert[nComp][m].vtVec * CompoVert[nComp][l].vtVec;
|
|
}
|
|
}
|
|
dDotAvarage /= (( nVertComp[nComp] * ( nVertComp[nComp] - 1)) / 2) ;
|
|
|
|
int nNumPar = 0 ;
|
|
for ( int m = 0 ; m < nVertComp[nComp] - 1 ; ++ m) {
|
|
for ( int l = m + 1 ; l < nVertComp[nComp] ; ++ l) {
|
|
if ( AreSameVectorExact( CompoVert[nComp][m].vtVec, CompoVert[nComp][l].vtVec))
|
|
++ nNumPar ;
|
|
}
|
|
}
|
|
bExtConfirmed = nNumPar == 3 ;
|
|
}
|
|
else if (nVertComp[nComp] < 5) {
|
|
bExtConfirmed = false ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Gestione casi speciali: tutti i vettori quasi equiversi tranne uno
|
|
if ( bExtConfirmed && ( nVertComp[nComp] == 4 || nVertComp[nComp] == 5)) {
|
|
int nParNum = 0 ;
|
|
for ( int m = 0 ; m < nVertComp[nComp] - 1 ; ++ m) {
|
|
for ( int l = m + 1 ; l < nVertComp[nComp] ; ++ l) {
|
|
if ( CompoVert[nComp][m].vtVec * CompoVert[nComp][l].vtVec > 0.9) {
|
|
++ nParNum ;
|
|
}
|
|
}
|
|
}
|
|
if ( nParNum == ( nVertComp[nComp] - 1) * ( nVertComp[nComp] - 2) / 2) {
|
|
int nNotParIndex = - 1 ;
|
|
for ( int m = 0 ; m < nVertComp[nComp] ; ++ m) {
|
|
int nNumDiff = 0 ;
|
|
for ( int l = 0 ; l < nVertComp[nComp] ; ++ l) {
|
|
if ( CompoVert[nComp][m].vtVec * CompoVert[nComp][l].vtVec <= 0.9) {
|
|
++ nNumDiff ;
|
|
}
|
|
}
|
|
if ( nNumDiff == nVertComp[nComp] - 1) {
|
|
nNotParIndex = m ;
|
|
break ;
|
|
}
|
|
}
|
|
int nPrevL = nNotParIndex > 0 ? nNotParIndex - 1 : nVertComp[nComp] - 1 ;
|
|
int nNextL = nNotParIndex < nVertComp[nComp] - 1 ? nNotParIndex + 1 : 0 ;
|
|
Vector3d vtPrevN = ( CompoVert[nComp][nNotParIndex].ptPApp - ptSol) ^ ( CompoVert[nComp][nPrevL].ptPApp - ptSol) ;
|
|
Vector3d vtNextN = ( CompoVert[nComp][nNextL].ptPApp - ptSol) ^ ( CompoVert[nComp][nNotParIndex].ptPApp - ptSol) ;
|
|
vtPrevN.Normalize() ;
|
|
vtNextN.Normalize() ;
|
|
if ( nNotParIndex >= 0 &&
|
|
CompoVert[nComp][nNotParIndex].vtVec * vtPrevN < - EPS_SMALL &&
|
|
CompoVert[nComp][nNotParIndex].vtVec * vtNextN < - EPS_SMALL) {
|
|
bExtConfirmed = false ;
|
|
}
|
|
}
|
|
}
|
|
|
|
// ExtMC confermato
|
|
if ( bExtConfirmed) {
|
|
int tOldCompo = VoxConf.nNumComp ;
|
|
++ VoxConf.nNumComp ;
|
|
ConComp& CurFan = VoxConf.Compo[tOldCompo] ;
|
|
CurFan.nVertNum = nVertComp[nComp] ;
|
|
for ( int nV = 0 ; nV < nVertComp[nComp] ; ++ nV) {
|
|
CurFan.CompVecField[nV] = CompoVert[nComp][nV] ;
|
|
}
|
|
for ( int nV = 0 ; nV < CurFan.nVertNum && CurFan.nVertNum >= 2 ; ++ nV) {
|
|
int nW = ( nV + 1) % CurFan.nVertNum ;
|
|
if ( AreSamePointApprox( CurFan.CompVecField[nV].ptPApp, CurFan.CompVecField[nW].ptPApp) &&
|
|
CurFan.CompVecField[nV].vtVec * CurFan.CompVecField[nW].vtVec > 0.9) {
|
|
Vector3d vtAvNorm = 0.5 * ( CurFan.CompVecField[nV].vtVec + CurFan.CompVecField[nW].vtVec) ;
|
|
vtAvNorm.Normalize() ;
|
|
CurFan.CompVecField[nV].ptPApp = 0.5 * ( CurFan.CompVecField[nV].ptPApp + CurFan.CompVecField[nW].ptPApp) ;
|
|
for ( int nL = nW ; nL < CurFan.nVertNum - 1 ; ++ nL) {
|
|
CurFan.CompVecField[nV].ptPApp = CurFan.CompVecField[nL + 1].ptPApp ;
|
|
CurFan.CompVecField[nL].vtVec = CurFan.CompVecField[nL + 1].vtVec ;
|
|
CurFan.CompVecField[nL].nPropIndex = CurFan.CompVecField[nL + 1].nPropIndex;
|
|
}
|
|
-- CurFan.nVertNum ;
|
|
}
|
|
}
|
|
CurFan.ptVert = ptSol ;
|
|
CurFan.vtNullSpace = vtNullSpace ;
|
|
CurFan.bInside = IsPointInsideVoxelApprox( i, j, k, ptSol, EPS_SMALL) ;
|
|
CurFan.bCorner = ( nFeatureType == CORNER) ;
|
|
}
|
|
// ExtMC non confermato
|
|
else
|
|
CreateSmoothTriangle( nIndex, nVertComp[nComp], CompoTriVert[nComp], true, VoxSmoothTria) ;
|
|
}
|
|
// Standard MC
|
|
else if ( m_nShape != BOX) {
|
|
CreateSmoothTriangle( nIndex, nVertComp[nComp], CompoTriVert[nComp], false, VoxSmoothTria) ;
|
|
}
|
|
}
|
|
// Se nel voxel abbiamo trovato feature
|
|
// aggiorniamo i contenitori.
|
|
if ( VoxConf.nNumComp > 0) {
|
|
VoxConf.i = i ;
|
|
VoxConf.j = j ;
|
|
VoxConf.k = k ;
|
|
int nIJK[3] = { i, j, k} ;
|
|
int nKey ;
|
|
GetVoxNFromIJK( i, j, k, nKey) ;
|
|
if ( IsAVoxelOnBoundary( nLimits, nIJK, true))
|
|
m_InterBlockVox[nBlock].emplace( nKey, VoxConf) ;
|
|
else
|
|
vVox.emplace( nKey, VoxConf) ;
|
|
}
|
|
// Se nel voxel abbiamo trovato componenti smooth
|
|
// aggiorno i contenitori
|
|
nVoxSmootSize = int( VoxSmoothTria.vTria.size()) ;
|
|
if ( nVoxSmootSize > nVoxSmootSizePrev) {
|
|
VoxSmoothTria.i = i ;
|
|
VoxSmoothTria.j = j ;
|
|
VoxSmoothTria.k = k ;
|
|
m_BlockSmoothTria[nBlock].emplace_back( VoxSmoothTria) ;
|
|
nVoxSmootSizePrev = nVoxSmootSize ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Se il solido è un parallelepipedo creo direttamente i triangoli grandi
|
|
if ( m_nShape == BOX) {
|
|
int nBlockIJK[3] ;
|
|
GetBlockIJKFromN( nBlock, nBlockIJK) ;
|
|
// Determino il primo nodo pieno della mappa
|
|
int nFirstVoxI, nFirstVoxJ, nFirstVoxK ;
|
|
GetFirstVoxIJK( nFirstVoxI, nFirstVoxJ, nFirstVoxK) ;
|
|
// Determino il primo nodo pieno della mappa
|
|
int nLastVoxI, nLastVoxJ, nLastVoxK ;
|
|
GetLastVoxIJK( nLastVoxI, nLastVoxJ, nLastVoxK) ;
|
|
// Costruisco i triangoli paralleli al piano YZ
|
|
if ( nBlockIJK[0] == 0 || nBlockIJK[0] + 1 == m_nFracLin[0]) {
|
|
// Indici dei dexel corrispondenti al confine inferiore fra voxel piatti e feature
|
|
int nDexMinJ = ( nBlockIJK[1] == 0 ? N_DEXVOXRATIO * ( nFirstVoxJ + 1) :
|
|
N_DEXVOXRATIO * nBlockIJK[1] * m_nVoxNumPerBlock) ;
|
|
int nDexMinK = ( nBlockIJK[2] == 0 ? N_DEXVOXRATIO * ( nFirstVoxK + 1) :
|
|
N_DEXVOXRATIO * nBlockIJK[2] * m_nVoxNumPerBlock) ;
|
|
// Indici dei dexel corrispondenti al confine superiore fra voxel piatti e feature
|
|
int nDexMaxJ = ( nBlockIJK[1] + 1 < int( m_nFracLin[1]) ? N_DEXVOXRATIO * ( nBlockIJK[1] + 1) * m_nVoxNumPerBlock :
|
|
N_DEXVOXRATIO * nLastVoxJ) ;
|
|
int nDexMaxK = ( nBlockIJK[2] + 1 < int( m_nFracLin[2]) ? N_DEXVOXRATIO * ( nBlockIJK[2] + 1) * m_nVoxNumPerBlock :
|
|
N_DEXVOXRATIO * nLastVoxK) ;
|
|
// Determino coordinate minime e massime dei punti dei triangoli
|
|
double dYMin = ( nDexMinJ + 0.5) * m_dStep ;
|
|
double dZMin = ( nDexMinK + 0.5) * m_dStep ;
|
|
double dYMax = ( nDexMaxJ + 0.5) * m_dStep ;
|
|
double dZMax = ( nDexMaxK + 0.5) * m_dStep ;
|
|
// Determino colore dei triangoli
|
|
int nFaceInfGrade = m_Values[1][nDexMinK * m_nNx[1] + nDexMinJ][0].nToolMin ;
|
|
int nFaceSupGrade = m_Values[1][nDexMinK * m_nNx[1] + nDexMinJ][0].nToolMax ;
|
|
// Piano di coordinata x inferiore: versore normale rivolto come X-
|
|
if ( nBlockIJK[0] == 0) {
|
|
Triangle3dEx trTria1, trTria2 ;
|
|
CreateBigTriangleYZ( dYMin, dYMax, dZMin, dZMax, m_dMinZ[1], false, nFaceInfGrade, trTria1, trTria2) ;
|
|
trTria1.ToGlob( m_MapFrame) ;
|
|
trTria2.ToGlob( m_MapFrame) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria1) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria2) ;
|
|
}
|
|
// Piano di coordinata x superiore: versore normale rivolto come X+
|
|
if ( nBlockIJK[0] + 1 == m_nFracLin[0]) {
|
|
Triangle3dEx trTria1, trTria2 ;
|
|
CreateBigTriangleYZ( dYMin, dYMax, dZMin, dZMax, m_dMaxZ[1], true, nFaceSupGrade, trTria1, trTria2) ;
|
|
trTria1.ToGlob( m_MapFrame) ;
|
|
trTria2.ToGlob( m_MapFrame) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria1) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria2) ;
|
|
}
|
|
}
|
|
// Costruisco i triangoli paralleli al piano XZ
|
|
if ( nBlockIJK[1] == 0 || nBlockIJK[1] + 1 == m_nFracLin[1]) {
|
|
// Indici dei dexel corrispondenti al confine inferiore fra voxel piatti e feature
|
|
int nDexMinI = ( nBlockIJK[0] == 0 ? N_DEXVOXRATIO * ( nFirstVoxI + 1) :
|
|
N_DEXVOXRATIO * nBlockIJK[0] * m_nVoxNumPerBlock) ;
|
|
int nDexMinK = ( nBlockIJK[2] == 0 ? N_DEXVOXRATIO * ( nFirstVoxK + 1) :
|
|
N_DEXVOXRATIO * nBlockIJK[2] * m_nVoxNumPerBlock) ;
|
|
// Indici dei dexel corrispondenti al confine superiore fra voxel piatti e feature
|
|
int nDexMaxI = ( nBlockIJK[0] + 1 < int( m_nFracLin[0]) ? N_DEXVOXRATIO * ( nBlockIJK[0] + 1) * m_nVoxNumPerBlock :
|
|
N_DEXVOXRATIO * nLastVoxI) ;
|
|
int nDexMaxK = ( nBlockIJK[2] + 1 < int( m_nFracLin[2]) ? N_DEXVOXRATIO * ( nBlockIJK[2] + 1) * m_nVoxNumPerBlock :
|
|
N_DEXVOXRATIO * nLastVoxK) ;
|
|
// Determino coordinate minime e massime dei punti dei triangoli
|
|
double dXMin = ( nDexMinI + 0.5) * m_dStep ;
|
|
double dZMin = ( nDexMinK + 0.5) * m_dStep ;
|
|
double dXMax = ( nDexMaxI + 0.5) * m_dStep ;
|
|
double dZMax = ( nDexMaxK + 0.5) * m_dStep ;
|
|
// Determino colore dei triangoli
|
|
int nFaceInfGrade = m_Values[2][nDexMinI * m_nNx[2] + nDexMinK][0].nToolMin ;
|
|
int nFaceSupGrade = m_Values[2][nDexMinI * m_nNx[2] + nDexMinK][0].nToolMax ;
|
|
// Piano di coordinata y inferiore: versore normale rivolto come Y-
|
|
if ( nBlockIJK[1] == 0) {
|
|
Triangle3dEx trTria1, trTria2 ;
|
|
CreateBigTriangleXZ( dXMin, dXMax, dZMin, dZMax, m_dMinZ[2], false, nFaceInfGrade, trTria1, trTria2) ;
|
|
trTria1.ToGlob( m_MapFrame) ;
|
|
trTria2.ToGlob( m_MapFrame) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria1) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria2) ;
|
|
}
|
|
// Piano di coordinata y superiore: versore normale rivolto come Y+
|
|
if ( nBlockIJK[1] + 1 == m_nFracLin[1]) {
|
|
Triangle3dEx trTria1, trTria2 ;
|
|
CreateBigTriangleXZ( dXMin, dXMax, dZMin, dZMax, m_dMaxZ[2], true, nFaceSupGrade, trTria1, trTria2) ;
|
|
trTria1.ToGlob( m_MapFrame) ;
|
|
trTria2.ToGlob( m_MapFrame) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria1) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria2) ;
|
|
}
|
|
}
|
|
// Costruisco i triangoli paralleli al piano XY
|
|
if ( nBlockIJK[2] == 0 || nBlockIJK[2] + 1 == m_nFracLin[2]) {
|
|
// Indici dei dexel corrispondenti al confine inferiore fra voxel piatti e feature
|
|
int nDexMinI = ( nBlockIJK[0] == 0 ? N_DEXVOXRATIO * ( nFirstVoxI + 1) :
|
|
N_DEXVOXRATIO * nBlockIJK[0] * m_nVoxNumPerBlock) ;
|
|
int nDexMinJ = ( nBlockIJK[1] == 0 ? N_DEXVOXRATIO * ( nFirstVoxJ + 1) :
|
|
N_DEXVOXRATIO * nBlockIJK[1] * m_nVoxNumPerBlock) ;
|
|
// Indici dei dexel corrispondenti al confine superiore fra voxel piatti e feature
|
|
int nDexMaxI = ( nBlockIJK[0] + 1 < int( m_nFracLin[0]) ? N_DEXVOXRATIO * ( nBlockIJK[0] + 1) * m_nVoxNumPerBlock :
|
|
N_DEXVOXRATIO * nLastVoxI) ;
|
|
int nDexMaxJ = ( nBlockIJK[1] + 1 < int( m_nFracLin[1]) ? N_DEXVOXRATIO * ( nBlockIJK[1] + 1) * m_nVoxNumPerBlock :
|
|
N_DEXVOXRATIO * nLastVoxJ) ;
|
|
// Determino coordinate minime e massime dei punti dei triangoli
|
|
double dXMin = ( nDexMinI + 0.5) * m_dStep ;
|
|
double dYMin = ( nDexMinJ + 0.5) * m_dStep ;
|
|
double dXMax = ( nDexMaxI + 0.5) * m_dStep ;
|
|
double dYMax = ( nDexMaxJ + 0.5) * m_dStep ;
|
|
// Determino colore dei triangoli
|
|
int nFaceInfGrade = m_Values[0][nDexMinJ * m_nNx[0] + nDexMinI][0].nToolMin ;
|
|
int nFaceSupGrade = m_Values[0][nDexMinJ * m_nNx[0] + nDexMinI][0].nToolMax ;
|
|
// Piano di coordinata Z inferiore: versore normale rivolto come Z-
|
|
if ( nBlockIJK[2] == 0) {
|
|
Triangle3dEx trTria1, trTria2 ;
|
|
CreateBigTriangleXY( dXMin, dXMax, dYMin, dYMax, m_dMinZ[0], false, nFaceInfGrade, trTria1, trTria2) ;
|
|
trTria1.ToGlob( m_MapFrame) ;
|
|
trTria2.ToGlob( m_MapFrame) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria1) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria2) ;
|
|
}
|
|
// Piano di coordinata Z superiore: versore normale rivolto come Z+
|
|
if ( nBlockIJK[2] + 1 == m_nFracLin[2]) {
|
|
Triangle3dEx trTria1, trTria2 ;
|
|
CreateBigTriangleXY( dXMin, dXMax, dYMin, dYMax, m_dMaxZ[0], true, nFaceSupGrade, trTria1, trTria2) ;
|
|
trTria1.ToGlob( m_MapFrame) ;
|
|
trTria2.ToGlob( m_MapFrame) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria1) ;
|
|
m_BlockBigTria[nBlock].emplace_back( trTria2) ;
|
|
}
|
|
}
|
|
}
|
|
// Processo i Voxel con possibile superficie piana
|
|
else {
|
|
ProcessVoxContXY( VoxContXYInf, nBlock, false) ;
|
|
ProcessVoxContXY( VoxContXYSup, nBlock, true) ;
|
|
ProcessVoxContYZ( VoxContYZInf, nBlock, false) ;
|
|
ProcessVoxContYZ( VoxContYZSup, nBlock, true) ;
|
|
ProcessVoxContXZ( VoxContXZInf, nBlock, false) ;
|
|
ProcessVoxContXZ( VoxContXZSup, nBlock, true) ;
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::RegulateFeaturesChain( vector<VoxelContainer>& vVecVox) const
|
|
{
|
|
// Ciclo sui blocchi
|
|
for ( int nBlock = 0 ; nBlock < int( m_nNumBlock) ; ++ nBlock) {
|
|
// Se il blocco è da aggiornare
|
|
if ( m_BlockToUpdate[nBlock]) {
|
|
// Ciclo sui voxel interni
|
|
for ( auto itVox = vVecVox[nBlock].begin() ; itVox != vVecVox[nBlock].end() ; ++ itVox) {
|
|
int nVox ;
|
|
Voxel& CurVox = itVox->second ;
|
|
GetVoxNFromIJK( CurVox.i, CurVox.j, CurVox.k, nVox) ;
|
|
// Ciclo sulle componenti
|
|
for ( int nComp = 0 ; nComp < CurVox.nNumComp ; ++ nComp) {
|
|
// Vertice fuori dal suo voxel
|
|
if ( ! CurVox.Compo[nComp].bInside) {
|
|
if ( IsPointInsideVoxelApprox( CurVox.i, CurVox.j, CurVox.k, CurVox.Compo[nComp].ptVert, 1.5 * EPS_SMALL))
|
|
continue ;
|
|
// Caso corner
|
|
if ( CurVox.Compo[nComp].bCorner) {
|
|
// Cerco i primi vicini
|
|
INTVECTOR vNearFirst, vBordNearFirst ;
|
|
FindAdjComp( vVecVox, nBlock, nVox, nComp, vNearFirst, vBordNearFirst) ;
|
|
// Ciclo sui primi vicini
|
|
int nInnSizeF = int( vNearFirst.size()) ;
|
|
int nBorSizeF = int( vBordNearFirst.size()) ;
|
|
for ( int nF = 0 ; nF < nInnSizeF + nBorSizeF ; nF += 3) {
|
|
// Indice e vettore corrente dei primi vicini
|
|
int nFCur = nF < nInnSizeF ? nF : nF - nInnSizeF ;
|
|
INTVECTOR vVecNFCur = nF < nInnSizeF ? vNearFirst : vBordNearFirst ;
|
|
// Cerco i secondi vicini
|
|
INTVECTOR vNearSecond, vBordNearSecond ;
|
|
FindAdjComp( vVecVox, vVecNFCur[nFCur], vVecNFCur[nFCur+1], vVecNFCur[nFCur+2],
|
|
vNearSecond, vBordNearSecond) ;
|
|
// Ciclo sui secondi vicini
|
|
int nInnSizeS = int( vNearSecond.size()) ;
|
|
int nBorSizeS = int( vBordNearSecond.size()) ;
|
|
for ( int nS = 0 ; nS < nInnSizeS + nBorSizeS ; nS += 3) {
|
|
// Indice e vettore corrente dei secondi vicini
|
|
int nNSCur = nS < nInnSizeS ? nS : nS - nInnSizeS ;
|
|
INTVECTOR vVecNSCur = nS < nInnSizeS ? vNearSecond : vBordNearSecond ;
|
|
// Escludo dai secondi i primi vicini
|
|
bool bFirst = false ;
|
|
if ( vVecNSCur[nNSCur] == nBlock &&
|
|
vVecNSCur[nNSCur+1] == nVox &&
|
|
vVecNSCur[nNSCur+2] == nComp)
|
|
bFirst = true ;
|
|
for ( int nOldF = 0 ; nOldF < nInnSizeF + nBorSizeF ; nOldF += 3) {
|
|
if ( nOldF < nInnSizeF) {
|
|
if ( vVecNSCur[nNSCur] == vNearFirst[nOldF] &&
|
|
vVecNSCur[nNSCur+1] == vNearFirst[nOldF+1] &&
|
|
vVecNSCur[nNSCur+2] == vNearFirst[nOldF+2])
|
|
bFirst = true ;
|
|
}
|
|
else {
|
|
if ( vVecNSCur[nNSCur] == vBordNearFirst[nOldF-nInnSizeF] &&
|
|
vVecNSCur[nNSCur+1] == vBordNearFirst[nOldF-nInnSizeF+1] &&
|
|
vVecNSCur[nNSCur+2] == vBordNearFirst[nOldF-nInnSizeF+2])
|
|
bFirst = true ;
|
|
}
|
|
}
|
|
// Se trovo un secondo fra i primi salto l'iterazione
|
|
if ( bFirst)
|
|
continue ;
|
|
// Se necessario regolarizzo la catena
|
|
Voxel& VoxNearFirst = nF < nInnSizeF ? vVecVox[vVecNFCur[nFCur]].find( vVecNFCur[nFCur+1])->second :
|
|
m_InterBlockVox[vVecNFCur[nFCur]].find( vVecNFCur[nFCur+1])->second ;
|
|
Voxel& VoxNearSecond = nS < nInnSizeS ? vVecVox[vVecNSCur[nNSCur]].find( vVecNSCur[nNSCur+1])->second :
|
|
m_InterBlockVox[vVecNSCur[nNSCur]].find( vVecNSCur[nNSCur+1])->second ;
|
|
Point3d ptCurV = CurVox.Compo[nComp].ptVert ;
|
|
Point3d ptFirst = VoxNearFirst.Compo[vVecNFCur[nFCur+2]].ptVert ;
|
|
Point3d ptSecond = VoxNearSecond.Compo[vVecNSCur[nNSCur+2]].ptVert ;
|
|
Vector3d vtF = ptFirst - ptCurV ;
|
|
Vector3d vtS = ptSecond - ptCurV ;
|
|
vtF.Normalize() ;
|
|
vtS.Normalize() ;
|
|
// Calcolo i baricentri dei vertici dei fan
|
|
Point3d ptCurBar ;
|
|
for ( int n = 0 ; n < CurVox.Compo[nComp].nVertNum ; ++ n)
|
|
ptCurBar += CurVox.Compo[nComp].CompVecField[n].ptPApp ;
|
|
ptCurBar /= CurVox.Compo[nComp].nVertNum ;
|
|
Point3d ptFirstBar ;
|
|
for ( int n = 0 ; n < VoxNearFirst.Compo[vVecNFCur[nFCur+2]].nVertNum ; ++ n)
|
|
ptFirstBar += VoxNearFirst.Compo[vVecNFCur[nFCur+2]].CompVecField[n].ptPApp ;
|
|
ptFirstBar /= VoxNearFirst.Compo[vVecNFCur[nFCur+2]].nVertNum ;
|
|
Point3d ptSecondBar ;
|
|
for ( int n = 0 ; n < VoxNearSecond.Compo[vVecNSCur[nNSCur+2]].nVertNum ; ++ n)
|
|
ptSecondBar += VoxNearSecond.Compo[vVecNSCur[nNSCur+2]].CompVecField[n].ptPApp ;
|
|
ptSecondBar /= VoxNearSecond.Compo[vVecNSCur[nNSCur+2]].nVertNum ;
|
|
Vector3d vtBarCF = ptFirstBar - ptCurBar ;
|
|
Vector3d vtBarCS = ptSecondBar - ptCurBar ;
|
|
Vector3d vtBarSF = ptSecondBar - ptFirstBar ;
|
|
vtBarCF.Normalize() ;
|
|
vtBarCS.Normalize() ;
|
|
vtBarSF.Normalize() ;
|
|
if ( vtBarCF * vtBarCS > 0.5 && vtBarCF * vtBarSF > 0.5 && vtBarCS * vtBarSF > 0.5 &&
|
|
vtF * vtS < - 0.9 && ! VoxNearFirst.Compo[vVecNFCur[nFCur+2]].bCorner) {
|
|
VoxNearFirst.Compo[vVecNFCur[nFCur+2]].ptVert = 0.5 * ( ptCurV + ptSecond) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
// Caso feature
|
|
else {
|
|
INTVECTOR vNearInn, vNearBord ;
|
|
FindAdjComp( vVecVox, nBlock, nVox, nComp, vNearInn, vNearBord) ;
|
|
int nSizeInn = int( vNearInn.size()) ;
|
|
int nSizeBord = int( vNearBord.size() );
|
|
if ( nSizeInn + nSizeBord == 6) {
|
|
const Voxel* pVoxSt = nullptr ;
|
|
const Voxel* pVoxEn = nullptr ;
|
|
if ( nSizeInn == 6) {
|
|
pVoxSt = &( vVecVox[vNearInn[0]].find( vNearInn[1])->second) ;
|
|
pVoxEn = &( vVecVox[vNearInn[3]].find( vNearInn[4])->second) ;
|
|
}
|
|
else if ( nSizeBord == 6) {
|
|
pVoxSt = &( m_InterBlockVox[vNearBord[0]].find( vNearBord[1])->second) ;
|
|
pVoxEn = &( m_InterBlockVox[vNearBord[3]].find( vNearBord[4])->second) ;
|
|
}
|
|
else {
|
|
pVoxSt = &( vVecVox[vNearInn[0]].find( vNearInn[1])->second) ;
|
|
pVoxEn = &( m_InterBlockVox[vNearBord[0]].find( vNearBord[1])->second) ;
|
|
}
|
|
Point3d ptPCur = CurVox.Compo[nComp].ptVert ;
|
|
Point3d ptPSt ;
|
|
Point3d ptPEn ;
|
|
if ( nSizeInn == 6) {
|
|
ptPSt = pVoxSt->Compo[vNearInn[2]].ptVert ;
|
|
ptPEn = pVoxEn->Compo[vNearInn[5]].ptVert ;
|
|
}
|
|
else if ( nSizeBord == 6) {
|
|
ptPSt = pVoxSt->Compo[vNearBord[2]].ptVert ;
|
|
ptPEn = pVoxEn->Compo[vNearBord[5]].ptVert ;
|
|
}
|
|
else {
|
|
ptPSt = pVoxSt->Compo[vNearInn[2]].ptVert ;
|
|
ptPEn = pVoxEn->Compo[vNearBord[2]].ptVert ;
|
|
}
|
|
Vector3d vtStCurr = ptPCur - ptPSt ;
|
|
Vector3d vtStEn = ptPEn - ptPSt ;
|
|
Vector3d vtCurrEn = ptPEn - ptPCur ;
|
|
vtStCurr.Normalize() ;
|
|
vtStEn.Normalize() ;
|
|
vtCurrEn.Normalize() ;
|
|
Point3d ptMid = 0.5 * ( ptPSt + ptPEn) ;
|
|
double dMidU = ( ptMid - ptPSt) * vtStEn ;
|
|
double dCurU = ( ptPCur - ptPSt) * vtStEn ;
|
|
Point3d ptNew ;
|
|
Point3d ptPLine ;
|
|
Vector3d vtDLine ;
|
|
if ( dMidU < dCurU) {
|
|
ptPLine = ptPEn ;
|
|
vtDLine = - vtCurrEn ;
|
|
ptNew = ptPEn + ( ptMid - ptPEn) * vtCurrEn * vtCurrEn ;
|
|
}
|
|
else {
|
|
ptPLine = ptPSt ;
|
|
vtDLine = vtStCurr ;
|
|
ptNew = ptPSt + ( ptMid - ptPSt) * vtStCurr * vtStCurr ;
|
|
}
|
|
Point3d ptCubeInf( ( CurVox.i * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( CurVox.j * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( CurVox.k * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
Point3d ptCubeSup( ( ( CurVox.i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( CurVox.j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( CurVox.k + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
double dU1, dU2 ;
|
|
if ( 1 - abs( vtStCurr * vtCurrEn ) < EPS_ZERO &&
|
|
IntersLineBox( ptPLine, vtDLine, ptCubeInf, ptCubeSup, dU1, dU2)) {
|
|
double dU = abs( dU1) < abs( dU2) ? dU1 + ( dU2 - dU1) / 2 : dU2 + ( dU1 - dU2) / 2 ;
|
|
ptNew = ptPLine + dU * vtDLine ;
|
|
}
|
|
bool bNewInside = IsPointInsideVoxelApprox( CurVox.i, CurVox.j, CurVox.k, ptNew, 0) ;
|
|
if ( abs( vtStCurr * vtStEn) > 0.95 && abs( vtStCurr * vtCurrEn) > 0.95 &&
|
|
abs( vtStEn * vtCurrEn) > 0.95) {
|
|
CurVox.Compo[nComp].ptVert = ptNew ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
// Ciclo sui voxel di frontiera
|
|
for ( auto itVox = m_InterBlockVox[nBlock].begin() ; itVox != m_InterBlockVox[nBlock].end() ; ++ itVox) {
|
|
int nVox ;
|
|
Voxel& CurVox = itVox->second ;
|
|
GetVoxNFromIJK( CurVox.i, CurVox.j, CurVox.k, nVox) ;
|
|
// Ciclo sulle componenti
|
|
for ( int nComp = 0 ; nComp < CurVox.nNumComp ; ++ nComp) {
|
|
// Vertice fuori dal suo voxel
|
|
if ( ! CurVox.Compo[nComp].bInside) {
|
|
if ( IsPointInsideVoxelApprox( CurVox.i, CurVox.j, CurVox.k, CurVox.Compo[nComp].ptVert, 1.5 * EPS_SMALL))
|
|
continue ;
|
|
// Caso feature
|
|
if ( ! CurVox.Compo[nComp].bCorner) {
|
|
INTVECTOR vNearInn, vNearBord ;
|
|
FindAdjComp( vVecVox, nBlock, nVox, nComp, vNearInn, vNearBord) ;
|
|
int nSizeInn = int( vNearInn.size()) ;
|
|
int nSizeBord = int( vNearBord.size() );
|
|
if ( nSizeInn + nSizeBord == 6) {
|
|
const Voxel* pVoxSt = nullptr ;
|
|
const Voxel* pVoxEn = nullptr ;
|
|
if ( nSizeInn == 6) {
|
|
pVoxSt = &( vVecVox[vNearInn[0]].find( vNearInn[1])->second) ;
|
|
pVoxEn = &( vVecVox[vNearInn[3]].find( vNearInn[4])->second) ;
|
|
}
|
|
else if ( nSizeBord == 6) {
|
|
pVoxSt = &( m_InterBlockVox[vNearBord[0]].find( vNearBord[1])->second) ;
|
|
pVoxEn = &( m_InterBlockVox[vNearBord[3]].find( vNearBord[4])->second) ;
|
|
}
|
|
else {
|
|
pVoxSt = &( vVecVox[vNearInn[0]].find( vNearInn[1])->second) ;
|
|
pVoxEn = &( m_InterBlockVox[vNearBord[0]].find( vNearBord[1])->second) ;
|
|
}
|
|
Point3d ptPCur = CurVox.Compo[nComp].ptVert ;
|
|
Point3d ptPSt ;
|
|
Point3d ptPEn ;
|
|
if ( nSizeInn == 6) {
|
|
ptPSt = pVoxSt->Compo[vNearInn[2]].ptVert ;
|
|
ptPEn = pVoxEn->Compo[vNearInn[5]].ptVert ;
|
|
}
|
|
else if ( nSizeBord == 6) {
|
|
ptPSt = pVoxSt->Compo[vNearBord[2]].ptVert ;
|
|
ptPEn = pVoxEn->Compo[vNearBord[5]].ptVert ;
|
|
}
|
|
else {
|
|
ptPSt = pVoxSt->Compo[vNearInn[2]].ptVert ;
|
|
ptPEn = pVoxEn->Compo[vNearBord[2]].ptVert ;
|
|
}
|
|
Vector3d vtStCurr = ptPCur - ptPSt ;
|
|
Vector3d vtStEn = ptPEn - ptPSt ;
|
|
Vector3d vtCurrEn = ptPEn - ptPCur ;
|
|
vtStCurr.Normalize() ;
|
|
vtStEn.Normalize() ;
|
|
vtCurrEn.Normalize() ;
|
|
Point3d ptMid = 0.5 * ( ptPSt + ptPEn) ;
|
|
double dMidU = ( ptMid - ptPSt) * vtStEn ;
|
|
double dCurU = ( ptPCur - ptPSt) * vtStEn ;
|
|
Point3d ptNew ;
|
|
Point3d ptPLine ;
|
|
Vector3d vtDLine ;
|
|
if ( dMidU < dCurU) {
|
|
ptPLine = ptPEn ;
|
|
vtDLine = - vtCurrEn ;
|
|
ptNew = ptPEn + ( ptMid - ptPEn) * vtCurrEn * vtCurrEn ;
|
|
}
|
|
else {
|
|
ptPLine = ptPSt ;
|
|
vtDLine = vtStCurr ;
|
|
ptNew = ptPSt + ( ptMid - ptPSt) * vtStCurr * vtStCurr ;
|
|
}
|
|
Point3d ptCubeInf( ( CurVox.i * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( CurVox.j * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( CurVox.k * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
Point3d ptCubeSup( ( ( CurVox.i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( CurVox.j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( CurVox.k + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
double dU1, dU2 ;
|
|
if ( 1 - abs( vtStCurr * vtCurrEn ) < EPS_ZERO &&
|
|
IntersLineBox( ptPLine, vtDLine, ptCubeInf, ptCubeSup, dU1, dU2)) {
|
|
double dU = abs( dU1) < abs( dU2) ? dU1 + ( dU2 - dU1) / 2 : dU2 + ( dU1 - dU2) / 2 ;
|
|
ptNew = ptPLine + dU * vtDLine ;
|
|
}
|
|
bool bNewInside = IsPointInsideVoxelApprox( CurVox.i, CurVox.j, CurVox.k, ptNew, 0) ;
|
|
if ( abs( vtStCurr * vtStEn) > 0.95 && abs( vtStCurr * vtCurrEn) > 0.95 &&
|
|
abs( vtStEn * vtCurrEn) > 0.95 /*&& bNewInside*/) {
|
|
CurVox.Compo[nComp].ptVert = ptNew ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// A partire dalla struttura dati di tipo voxel con sharp features crea i triangoli corrispondenti.
|
|
// Accetta come parametri il numero del blocco, l'unordered map di voxel e l'apposito contenitore
|
|
// ove salva i triangoli.
|
|
bool
|
|
VolZmap::CreateSharpFeatureTriangle( int nBlock, const VoxelContainer& vVoxel) const
|
|
{
|
|
// Calcolo i limiti sugli indici dei voxel del blocco
|
|
// Vettore indici i,j,k del blocco
|
|
int nBlockIJK[3] ;
|
|
GetBlockIJKFromN( nBlock, nBlockIJK) ;
|
|
// Vettore limiti sugli indici dei voxel del blocco
|
|
int nLimits[6] ;
|
|
GetBlockLimitsIJK( nBlockIJK, nLimits) ;
|
|
|
|
// Ciclo sui voxel interni
|
|
for ( auto it = vVoxel.begin() ; it != vVoxel.end() ; ++ it) {
|
|
int tOldSize = int( m_BlockSharpTria[nBlock].size()) ;
|
|
m_BlockSharpTria[nBlock].resize( tOldSize + 1) ;
|
|
m_BlockSharpTria[nBlock][tOldSize].i = it->second.i ;
|
|
m_BlockSharpTria[nBlock][tOldSize].j = it->second.j ;
|
|
m_BlockSharpTria[nBlock][tOldSize].k = it->second.k ;
|
|
// Ciclo sulle componenti connesse del voxel
|
|
for ( int nComp = 0 ; nComp < it->second.nNumComp ; ++ nComp) {
|
|
m_BlockSharpTria[nBlock][tOldSize].ptCompoVert.emplace_back( it->second.Compo[nComp].ptVert) ;
|
|
m_BlockSharpTria[nBlock][tOldSize].ptCompoVert.back().ToGlob( m_MapFrame) ;
|
|
int tOldCompNum = int( m_BlockSharpTria[nBlock][tOldSize].vCompoTria.size()) ;
|
|
m_BlockSharpTria[nBlock][tOldSize].vCompoTria.resize( tOldCompNum + 1) ;
|
|
m_BlockSharpTria[nBlock][tOldSize].vbFlipped.resize( tOldCompNum + 1) ;
|
|
// ciclo sui vertici della componente connessa
|
|
int nNumVert = it->second.Compo[nComp].nVertNum ;
|
|
for ( int nVert = 0 ; nVert < nNumVert ; ++ nVert) {
|
|
int nNextVert = ( nVert + 1 < nNumVert ? nVert + 1 : 0) ;
|
|
// Definisco il triangolo
|
|
Triangle3dEx CurrTri ;
|
|
CurrTri.Set( it->second.Compo[nComp].ptVert,
|
|
it->second.Compo[nComp].CompVecField[nNextVert].ptPApp,
|
|
it->second.Compo[nComp].CompVecField[nVert].ptPApp) ;
|
|
// Setto il numero di utensile ai vertici di base del fan
|
|
CurrTri.SetAttrib( 1, it->second.Compo[nComp].CompVecField[nNextVert].nPropIndex) ;
|
|
CurrTri.SetAttrib( 2, it->second.Compo[nComp].CompVecField[nVert].nPropIndex) ;
|
|
// Setto il numero di utensile al triangolo nel complesso
|
|
if ( CurrTri.GetAttrib( 1) < 0 ||
|
|
CurrTri.GetAttrib( 2) < 0)
|
|
CurrTri.SetGrade( - 1) ;
|
|
else if ( CurrTri.GetAttrib( 1) > 0 ||
|
|
CurrTri.GetAttrib( 2) > 0)
|
|
CurrTri.SetGrade( 1) ;
|
|
else
|
|
CurrTri.SetGrade( 0) ;
|
|
// Setto le normali a ogni vertice
|
|
CurrTri.SetVertexNorm( 1, it->second.Compo[nComp].CompVecField[nNextVert].vtVec) ;
|
|
CurrTri.SetVertexNorm( 2, it->second.Compo[nComp].CompVecField[nVert].vtVec) ;
|
|
// Setto i flag a false
|
|
CurrTri.SetEdgeFlag( 0, false) ;
|
|
CurrTri.SetEdgeFlag( 1, false) ;
|
|
CurrTri.SetEdgeFlag( 2, false) ;
|
|
// Valido il triangolo
|
|
if ( ! CurrTri.Validate( true))
|
|
CurrTri.SetVertexNorm( 0, 0.5 * ( CurrTri.GetVertexNorm( 1) + CurrTri.GetVertexNorm( 2))) ;
|
|
m_BlockSharpTria[nBlock][tOldSize].vCompoTria[tOldCompNum].emplace_back( CurrTri) ;
|
|
m_BlockSharpTria[nBlock][tOldSize].vCompoTria[tOldCompNum].back().ToGlob( m_MapFrame) ;
|
|
int nTri = int( m_BlockSharpTria[nBlock][tOldSize].vbFlipped[tOldCompNum].size()) ;
|
|
m_BlockSharpTria[nBlock][tOldSize].vbFlipped[tOldCompNum].resize( nTri + 1) ;
|
|
m_BlockSharpTria[nBlock][tOldSize].vbFlipped[tOldCompNum][nTri] = false ;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Pulisco il contenitore dei voxel di frontiera
|
|
m_InterBlockOriginalSharpTria[nBlock].clear() ;
|
|
// Ciclo sui voxel di frontiera
|
|
for ( auto itVox = m_InterBlockVox[nBlock].begin() ; itVox != m_InterBlockVox[nBlock].end() ; ++ itVox) {
|
|
// Indici del voxel
|
|
int nVoxIJK[3] = { itVox->second.i,
|
|
itVox->second.j,
|
|
itVox->second.k} ;
|
|
// Ridimensiono il contenitore dei triangoli interni
|
|
int nOldSizeInn = int( m_BlockSharpTria[nBlock].size()) ;
|
|
m_BlockSharpTria[nBlock].resize( nOldSizeInn + 1) ;
|
|
m_BlockSharpTria[nBlock][nOldSizeInn].i = nVoxIJK[0] ;
|
|
m_BlockSharpTria[nBlock][nOldSizeInn].j = nVoxIJK[1] ;
|
|
m_BlockSharpTria[nBlock][nOldSizeInn].k = nVoxIJK[2] ;
|
|
// Ridimensiono il contenitore dei triangoli di frontiera
|
|
int nOldSizeBor = int( m_InterBlockOriginalSharpTria[nBlock].size()) ;
|
|
m_InterBlockOriginalSharpTria[nBlock].resize( nOldSizeBor + 1) ;
|
|
m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].i = nVoxIJK[0] ;
|
|
m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].j = nVoxIJK[1] ;
|
|
m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].k = nVoxIJK[2] ;
|
|
// Ciclo sulle componenti connesse del voxel
|
|
for ( int nComp = 0 ; nComp < itVox->second.nNumComp ; ++ nComp) {
|
|
bool bNewCompInn = true ;
|
|
bool bNewCompBor = true ;
|
|
// ciclo sui vertici della componente connessa
|
|
int nNumVert = itVox->second.Compo[nComp].nVertNum ;
|
|
for ( int nVert = 0 ; nVert < nNumVert ; ++ nVert) {
|
|
int nNextVert = ( nVert + 1 < nNumVert ? nVert + 1 : 0) ;
|
|
// Definisco il triangolo
|
|
Triangle3dEx CurrTri ;
|
|
CurrTri.Set( itVox->second.Compo[nComp].ptVert,
|
|
itVox->second.Compo[nComp].CompVecField[nNextVert].ptPApp,
|
|
itVox->second.Compo[nComp].CompVecField[nVert].ptPApp) ;
|
|
// Setto il numero di utensile ai vertici di base del fan
|
|
CurrTri.SetAttrib( 1, itVox->second.Compo[nComp].CompVecField[nNextVert].nPropIndex) ;
|
|
CurrTri.SetAttrib( 2, itVox->second.Compo[nComp].CompVecField[nVert].nPropIndex) ;
|
|
// Setto il numero di utensile al triangolo nel complesso
|
|
if ( CurrTri.GetAttrib( 1) < 0 ||
|
|
CurrTri.GetAttrib( 2) < 0)
|
|
CurrTri.SetGrade( - 1) ;
|
|
else if ( CurrTri.GetAttrib( 1) > 0 ||
|
|
CurrTri.GetAttrib( 2) > 0)
|
|
CurrTri.SetGrade( 1) ;
|
|
else
|
|
CurrTri.SetGrade( 0) ;
|
|
// Setto le normali a ogni vertice
|
|
CurrTri.SetVertexNorm( 1, itVox->second.Compo[nComp].CompVecField[nNextVert].vtVec) ;
|
|
CurrTri.SetVertexNorm( 2, itVox->second.Compo[nComp].CompVecField[nVert].vtVec) ;
|
|
// Setto i flag a false
|
|
CurrTri.SetEdgeFlag( 0, false) ;
|
|
CurrTri.SetEdgeFlag( 1, false) ;
|
|
CurrTri.SetEdgeFlag( 2, false) ;
|
|
// Valido il triangolo
|
|
if ( ! CurrTri.Validate( true))
|
|
CurrTri.SetVertexNorm( 0, 0.5 * ( CurrTri.GetVertexNorm( 1) + CurrTri.GetVertexNorm( 2)));
|
|
// Smisto i triangoli fra di frontiera e interni
|
|
bool bTriOnBorder = IsTriangleOnBorder( CurrTri, nLimits, nVoxIJK) ;
|
|
// Triangolo di frontiera
|
|
if ( bTriOnBorder) {
|
|
if ( bNewCompBor) {
|
|
Point3d ptVert = itVox->second.Compo[nComp].ptVert ;
|
|
ptVert.ToGlob( m_MapFrame) ;
|
|
m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].ptCompoVert.emplace_back( ptVert) ;
|
|
int tOldComp = int( m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vCompoTria.size()) ;
|
|
m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vCompoTria.resize( tOldComp + 1) ;
|
|
m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vbFlipped.resize( tOldComp + 1) ;
|
|
bNewCompBor = false ;
|
|
}
|
|
int tCurrSz = int( m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vCompoTria.size()) ;
|
|
m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vCompoTria[tCurrSz - 1].emplace_back( CurrTri) ;
|
|
m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vCompoTria[tCurrSz - 1].back().ToGlob( m_MapFrame) ;
|
|
int tTri =int( m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vbFlipped[tCurrSz - 1].size()) ;
|
|
m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vbFlipped[tCurrSz - 1].resize( tTri + 1) ;
|
|
m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vbFlipped[tCurrSz - 1][tTri] = false ;
|
|
}
|
|
// Triangolo interno
|
|
else {
|
|
if ( bNewCompInn) {
|
|
Point3d ptVert = itVox->second.Compo[nComp].ptVert ;
|
|
ptVert.ToGlob( m_MapFrame) ;
|
|
m_BlockSharpTria[nBlock][nOldSizeInn].ptCompoVert.emplace_back( ptVert) ;
|
|
int tOldComp = int( m_BlockSharpTria[nBlock][nOldSizeInn].vCompoTria.size()) ;
|
|
m_BlockSharpTria[nBlock][nOldSizeInn].vCompoTria.resize( tOldComp + 1) ;
|
|
m_BlockSharpTria[nBlock][nOldSizeInn].vbFlipped.resize( tOldComp + 1) ;
|
|
bNewCompInn = false ;
|
|
}
|
|
int tCurrSz = int( m_BlockSharpTria[nBlock][nOldSizeInn].vCompoTria.size()) ;
|
|
m_BlockSharpTria[nBlock][nOldSizeInn].vCompoTria[tCurrSz - 1].emplace_back( CurrTri) ;
|
|
m_BlockSharpTria[nBlock][nOldSizeInn].vCompoTria[tCurrSz - 1].back().ToGlob( m_MapFrame) ;
|
|
int tTri = int( m_BlockSharpTria[nBlock][nOldSizeInn].vbFlipped[tCurrSz - 1].size()) ;
|
|
m_BlockSharpTria[nBlock][nOldSizeInn].vbFlipped[tCurrSz - 1].resize( tTri + 1) ;
|
|
m_BlockSharpTria[nBlock][nOldSizeInn].vbFlipped[tCurrSz - 1][tTri] = false ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::CreateSmoothTriangle(int nIndex, int nVertNum, AppliedVector TriVert[], bool bWasSharp, SmoothTriaStruct& VoxSmoothTria) const
|
|
{
|
|
TRIA3DEXVECTOR vTria ;
|
|
// Costruzione dei triangoli
|
|
for ( int TriIndex = 0 ; TriIndex < (nVertNum - 2) * 3 ; TriIndex += 3) {
|
|
// Il triangolo è pronto
|
|
Triangle3dEx CurrentTriangle ;
|
|
CurrentTriangle.Set( TriVert[TriIndex].ptPApp,
|
|
TriVert[TriIndex + 1].ptPApp,
|
|
TriVert[TriIndex + 2].ptPApp) ;
|
|
// Setto il numero di utensile (conta solo positivo, nullo o negativo)
|
|
int nTool0 = Clamp( TriVert[TriIndex].nPropIndex, - 1, 1) ;
|
|
int nTool1 = Clamp( TriVert[TriIndex + 1].nPropIndex, - 1, 1) ;
|
|
int nTool2 = Clamp( TriVert[TriIndex + 2].nPropIndex, - 1, 1) ;
|
|
int nGrade = 0 ;
|
|
if ( nTool0 == -1 || nTool1 == -1 || nTool2 == -1)
|
|
nGrade = -1 ;
|
|
else if ( nTool0 == 1 || nTool1 == 1 || nTool2 == 1)
|
|
nGrade = 1 ;
|
|
CurrentTriangle.SetGrade( nGrade) ;
|
|
CurrentTriangle.SetAttrib( 0, nTool0) ;
|
|
CurrentTriangle.SetAttrib( 1, nTool1) ;
|
|
CurrentTriangle.SetAttrib( 2, nTool2) ;
|
|
// Valido il triangolo e setto le normali del campo vettoriale ai corrispondenti vertici
|
|
if ( CurrentTriangle.Validate( true)) {
|
|
double dCosAngThreshold = ( bWasSharp ? 0.7 : 0.5) ;
|
|
for ( int nV = 0 ; nV < 3 ; ++ nV) {
|
|
const Vector3d& vtVertNorm = TriVert[TriIndex + nV].vtVec ;
|
|
if ( CurrentTriangle.GetN() * vtVertNorm > dCosAngThreshold)
|
|
CurrentTriangle.SetVertexNorm( nV, vtVertNorm) ;
|
|
else
|
|
CurrentTriangle.SetVertexNorm( nV, CurrentTriangle.GetN()) ;
|
|
}
|
|
}
|
|
// Riporto le coordinate nel sistema in cui è immerso lo Zmap
|
|
CurrentTriangle.ToGlob( m_MapFrame) ;
|
|
// Aggiungo alla lista
|
|
vTria.emplace_back( CurrentTriangle) ;
|
|
}
|
|
// Controllo i colori di configurazioni 2 e 8
|
|
if ( nAllConfig[nIndex] == 2 || nAllConfig[nIndex] == 8) {
|
|
int nFst1 = - 1 ;
|
|
int nFst2 = - 1 ;
|
|
int nSec1 = - 1 ;
|
|
int nSec2 = - 1 ;
|
|
for ( int nF = 0 ; nF < 3 ; ++ nF) {
|
|
for ( int nS = 0 ; nS < 3 ; ++ nS) {
|
|
if ( AreSamePointExact( vTria[0].GetP( nF), vTria[1].GetP( nS))) {
|
|
if ( nFst1 < 0)
|
|
nFst1 = nF ;
|
|
else
|
|
nFst2 = nF ;
|
|
if ( nSec1 < 0)
|
|
nSec1 = nS ;
|
|
else
|
|
nSec2 = nS ;
|
|
}
|
|
}
|
|
}
|
|
int nFst3 = -1 ;
|
|
int nSec3 = -1 ;
|
|
for ( int n = 0 ; n < 3 ; ++n) {
|
|
if ( n != nFst1 && n != nFst2)
|
|
nFst3 = n ;
|
|
if ( n != nSec1 && n != nSec2)
|
|
nSec3 = n ;
|
|
}
|
|
if ( vTria[0].GetAttrib( nFst1) != vTria[0].GetAttrib( nFst2) &&
|
|
vTria[0].GetAttrib( nFst3) == vTria[1].GetAttrib( nSec3)) {
|
|
vTria[0].SetP( nFst1, vTria[1].GetP( nSec3)) ;
|
|
vTria[0].SetAttrib( nFst1, vTria[1].GetAttrib( nSec3)) ;
|
|
vTria[0].SetVertexNorm( nFst1, vTria[1].GetVertexNorm( nSec3)) ;
|
|
int nFirstGrade = 0 ;
|
|
if ( vTria[0].GetAttrib(0) == -1 || vTria[0].GetAttrib(1) == -1 || vTria[0].GetAttrib(2) == -1)
|
|
nFirstGrade = -1 ;
|
|
else if ( vTria[0].GetAttrib(0) == 1 || vTria[0].GetAttrib(1) == 1 || vTria[0].GetAttrib(2) == 1)
|
|
nFirstGrade = 1 ;
|
|
vTria[0].SetGrade( nFirstGrade) ;
|
|
vTria[0].Validate( true) ;
|
|
vTria[1].SetP( nSec2, vTria[0].GetP( nFst3)) ;
|
|
vTria[1].SetAttrib( nSec2, vTria[0].GetAttrib( nFst3)) ;
|
|
vTria[1].SetVertexNorm( nSec2, vTria[0].GetVertexNorm( nFst3)) ;
|
|
int nSecondGrade = 0 ;
|
|
if ( vTria[1].GetAttrib(0) == -1 || vTria[1].GetAttrib(1) == -1 || vTria[1].GetAttrib(2) == -1)
|
|
nSecondGrade = -1 ;
|
|
else if ( vTria[1].GetAttrib(0) == 1 || vTria[1].GetAttrib(1) == 1 || vTria[1].GetAttrib(2) == 1)
|
|
nSecondGrade = 1 ;
|
|
vTria[1].SetGrade( nSecondGrade) ;
|
|
vTria[1].Validate( true) ;
|
|
}
|
|
else if ( vTria[0].GetAttrib(nFst1) == vTria[0].GetAttrib(nFst2) &&
|
|
vTria[0].GetAttrib(nFst3) != vTria[1].GetAttrib(nSec3)) {
|
|
int nFirstGrade = 0 ;
|
|
if ( vTria[0].GetAttrib(0) == -1 || vTria[0].GetAttrib(1) == -1 || vTria[0].GetAttrib(2) == -1)
|
|
nFirstGrade = -1 ;
|
|
else if ( vTria[0].GetAttrib(0) == 1 || vTria[0].GetAttrib(1) == 1 || vTria[0].GetAttrib(2) == 1)
|
|
nFirstGrade = 1 ;
|
|
vTria[0].SetGrade( nFirstGrade) ;
|
|
int nSecondGrade = 0 ;
|
|
if ( vTria[1].GetAttrib(0) == -1 || vTria[1].GetAttrib(1) == -1 || vTria[1].GetAttrib(2) == -1)
|
|
nSecondGrade = -1 ;
|
|
else if ( vTria[1].GetAttrib(0) == 1 || vTria[1].GetAttrib(1) == 1 || vTria[1].GetAttrib(2) == 1)
|
|
nSecondGrade = 1 ;
|
|
vTria[1].SetGrade( nSecondGrade) ;
|
|
}
|
|
}
|
|
for ( int nT = 0 ; nT < int( vTria.size()) ; ++nT)
|
|
VoxSmoothTria.vTria.emplace_back( vTria[nT]) ;
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Esegue il flipping dei triangoli contenuti nel TriHolder,
|
|
// bGraph indica se chiamata per grafica o per calcolo di profondità.
|
|
bool
|
|
VolZmap::FlipEdgesII( int nBlock) const
|
|
{
|
|
// Numero di voxel in cui si presentano sharp feature
|
|
int nVoxelNum = int( m_BlockSharpTria[nBlock].size()) ;
|
|
|
|
// Ciclo sui voxel con sharp feature
|
|
for ( int n1 = 0 ; n1 < nVoxelNum ; ++ n1) {
|
|
SharpTriaStruct& SharpTria1 = m_BlockSharpTria[nBlock][n1] ;
|
|
for ( int n2 = n1 ; n2 < nVoxelNum ; ++ n2) {
|
|
SharpTriaStruct& SharpTria2 = m_BlockSharpTria[nBlock][n2] ;
|
|
// Se non adiacenti o coincidenti vado oltre
|
|
if ( abs( SharpTria2.i - SharpTria1.i) > 1 ||
|
|
abs( SharpTria2.j - SharpTria1.j) > 1 ||
|
|
abs( SharpTria2.k - SharpTria1.k) > 1 )
|
|
continue ;
|
|
// Ciclo sulle componenti connesse del primo voxel
|
|
int nNumCompo1 = int( SharpTria1.ptCompoVert.size()) ;
|
|
for ( int nCompo1 = 0 ; nCompo1 < nNumCompo1 ; ++ nCompo1) {
|
|
TRIA3DEXVECTOR& vTria1 = SharpTria1.vCompoTria[nCompo1] ;
|
|
// Numero di triangoli della componente connessa
|
|
int nTriNum1 = int( vTria1.size()) ;
|
|
// Ciclo sulle componenti connesse del secondo voxel
|
|
int nNumCompo2 = int( SharpTria2.ptCompoVert.size()) ;
|
|
int nCompo2 = ( n1 == n2 ? nCompo1 + 1 : 0) ;
|
|
for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) {
|
|
TRIA3DEXVECTOR& vTria2 = SharpTria2.vCompoTria[nCompo2] ;
|
|
// Numero di triangoli della componente connessa
|
|
int nTriNum2 = int( vTria2.size()) ;
|
|
for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) {
|
|
bool bModified = false ;
|
|
for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) {
|
|
// Punti che devono essere in comune fra i due triangoli
|
|
const Point3d& ptP11 = vTria1[nTri1].GetP( 1) ;
|
|
const Point3d& ptP12 = vTria1[nTri1].GetP( 2) ;
|
|
const Point3d& ptP21 = vTria2[nTri2].GetP( 1) ;
|
|
const Point3d& ptP22 = vTria2[nTri2].GetP( 2) ;
|
|
// I triangoli sono da flippare
|
|
if ( AreSamePointEpsilon( ptP11, ptP22, EPS_ZERO) &&
|
|
AreSamePointEpsilon( ptP12, ptP21, EPS_ZERO) &&
|
|
! ( SharpTria1.vbFlipped[nCompo1][nTri1] ||
|
|
SharpTria2.vbFlipped[nCompo2][nTri2])) {
|
|
// Assegno l'array dei punti di contorno
|
|
Point3d vPnt[4] ;
|
|
vPnt[0] = vTria1[nTri1].GetP( 1) ;
|
|
vPnt[1] = SharpTria1.ptCompoVert[nCompo1] ;
|
|
vPnt[2] = vTria1[nTri1].GetP( 2) ;
|
|
vPnt[3] = SharpTria2.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 = vTria1[nTri1] ;
|
|
Triangle3d trT2 = vTria2[nTri2] ;
|
|
int nVert1, nVert2 ;
|
|
// Determino gli indici dei punti sharp-feature
|
|
for ( int nP = 0 ; nP < 3 ; ++ nP) {
|
|
if ( nP != 1 && nP != 2)
|
|
nVert1 = nP ;
|
|
if ( nP != 2 && nP != 1)
|
|
nVert2 = nP ;
|
|
}
|
|
trT1.SetP( 1, trT2.GetP( nVert2)) ;
|
|
trT2.SetP( 1, trT1.GetP( nVert1)) ;
|
|
trT1.Validate( true) ;
|
|
trT2.Validate( true) ;
|
|
bInv = ( trT1.GetN() * trT2.GetN() < 0) ;
|
|
}
|
|
// Se non vi è inversione eseguo il flipping
|
|
if ( ! bInv) {
|
|
// Vertice condiviso fra nTri1 e quello del suo fan
|
|
int nCol1 = vTria1[nTri1].GetAttrib( 2) ;
|
|
int nCol2 = vTria2[nTri2].GetAttrib( 2) ;
|
|
// Modifico i punti e gli indici
|
|
vTria1[nTri1].SetP( 1, SharpTria2.ptCompoVert[nCompo2]) ;
|
|
vTria2[nTri2].SetP( 1, SharpTria1.ptCompoVert[nCompo1]) ;
|
|
vTria1[nTri1].SetGrade( nCol1) ;
|
|
vTria2[nTri2].SetGrade( nCol2) ;
|
|
// Setto le normali
|
|
Vector3d vtN1 = vTria1[nTri1].GetVertexNorm( 2) ;
|
|
Vector3d vtN2 = vTria2[nTri2].GetVertexNorm( 2) ;
|
|
vTria1[nTri1].SetVertexNorm( 0, V_NULL) ;
|
|
vTria1[nTri1].SetVertexNorm( 1, V_NULL) ;
|
|
vTria2[nTri2].SetVertexNorm( 0, V_NULL) ;
|
|
vTria2[nTri2].SetVertexNorm( 1, V_NULL) ;
|
|
// Valido i triangoli
|
|
vTria1[nTri1].Validate( true) ;
|
|
vTria2[nTri2].Validate( true) ;
|
|
// Setto i bFlags
|
|
if ( ! AreSameVectorApprox( vTria1[nTri1].GetN(), vTria2[nTri2].GetN())) {
|
|
vTria1[nTri1].SetEdgeFlag( 0, true) ;
|
|
}
|
|
// Setto i triangoli come flippati
|
|
SharpTria1.vbFlipped[nCompo1][nTri1] = true ;
|
|
SharpTria2.vbFlipped[nCompo2][nTri2] = true ;
|
|
// Avvenuto flipping
|
|
bModified = true ;
|
|
break ;
|
|
}
|
|
else {
|
|
// In questo caso i due triangoli sono necessariamente su un piano,
|
|
// quindi hanno normali concordi. A entrambi assegno il colore del vertice con
|
|
// normale più concorde a quella del triangolo
|
|
double dDotVec[4] ;
|
|
dDotVec[0] = vTria1[nTri1].GetN() * vTria1[nTri1].GetVertexNorm( 1) ;
|
|
dDotVec[1] = vTria2[nTri2].GetN() * vTria2[nTri2].GetVertexNorm( 2) ;
|
|
dDotVec[2] = vTria1[nTri1].GetN() * vTria1[nTri1].GetVertexNorm( 2) ;
|
|
dDotVec[3] = vTria2[nTri2].GetN() * vTria2[nTri2].GetVertexNorm( 1) ;
|
|
// Cerco il massimo dei prodotti scalari
|
|
int nMaxPos = 0 ;
|
|
double dMaxDot = - 1 ;
|
|
for ( int nPos = 0 ; nPos < 4 && dMaxDot < 1 ; ++ nPos) {
|
|
if ( dDotVec[nPos] > dMaxDot) {
|
|
dMaxDot = dDotVec[nPos] ;
|
|
nMaxPos = nPos ;
|
|
}
|
|
}
|
|
// Trovo il colore associato al vertice di massimo prodotto scalare
|
|
int nCol ;
|
|
switch ( nMaxPos) {
|
|
case 0 :
|
|
nCol = vTria1[nTri1].GetAttrib( 1) ;
|
|
break ;
|
|
case 1 :
|
|
nCol = vTria2[nTri2].GetAttrib( 2) ;
|
|
break ;
|
|
case 2 :
|
|
nCol = vTria1[nTri1].GetAttrib( 2) ;
|
|
break ;
|
|
case 3 :
|
|
nCol = vTria2[nTri2].GetAttrib( 1) ;
|
|
break ;
|
|
}
|
|
// Assegno il colore ai triangoli
|
|
vTria1[nTri1].SetGrade( nCol) ;
|
|
vTria2[nTri2].SetGrade( nCol) ;
|
|
// Setto i flag dei vertici
|
|
double dDist02 = ( vTria1[nTri1].GetP( 0) - vTria1[nTri1].GetP( 2)).Len() +
|
|
( vTria2[nTri2].GetP( 1) - vTria2[nTri2].GetP( 0)).Len() ;
|
|
double dDist20 = ( vTria1[nTri1].GetP( 1) - vTria1[nTri1].GetP( 0)).Len() +
|
|
( vTria2[nTri2].GetP( 0) - vTria2[nTri2].GetP( 2)).Len() ;
|
|
if ( dDist02 > dDist20) {
|
|
vTria1[nTri1].SetEdgeFlag( 0, true) ;
|
|
vTria2[nTri2].SetEdgeFlag( 2, true) ;
|
|
}
|
|
else {
|
|
vTria1[nTri1].SetEdgeFlag( 2, true) ;
|
|
vTria2[nTri2].SetEdgeFlag( 0, true) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if ( bModified)
|
|
break ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Ciclo sui voxel con sharp feature
|
|
for ( int nVox = 0 ; nVox < nVoxelNum ; ++nVox) {
|
|
SharpTriaStruct& SharpTria = m_BlockSharpTria[nBlock][nVox] ;
|
|
// Ciclo sulle componenti connesse del voxel
|
|
int nCompNum = int( SharpTria.ptCompoVert.size()) ;
|
|
for ( int nComp = 0 ; nComp < nCompNum ; ++ nComp) {
|
|
TRIA3DEXVECTOR& vTria = SharpTria.vCompoTria[nComp] ;
|
|
// Ciclo sui triangoli della componente connessa
|
|
int nTriaNum = int( vTria.size()) ;
|
|
for ( int nTria = 0 ; nTria < nTriaNum ; ++ nTria) {
|
|
// Lunghezza del segmento opposto al vertice
|
|
double dDist = Dist( vTria[nTria].GetP( 1), vTria[nTria].GetP( 2)) ;
|
|
// Se è quasi nullo
|
|
if ( dDist < 2 * EPS_SMALL && ! SharpTria.vbFlipped[nComp][nTria]) {
|
|
// Cerco triangolo adiacente con normale più vicina
|
|
int nPrevTria = ( nTria - 1 >= 0 ? nTria - 1 : nTriaNum - 1) ;
|
|
int nNextTria = ( nTria + 1 < nTriaNum ? nTria + 1 : 0) ;
|
|
double dDotPrev = vTria[nTria].GetN() * vTria[nPrevTria].GetN() ;
|
|
double dDotNext = vTria[nTria].GetN() * vTria[nNextTria].GetN() ;
|
|
if ( dDotPrev > dDotNext) {
|
|
double dAdjDist = Dist( vTria[nPrevTria].GetP( 1), vTria[nPrevTria].GetP( 2)) ;
|
|
// Se normale sufficientemente vicina e lato opposto al vertice abbastanza lungo assegno normale al vertice del triangolo corrente
|
|
if ( dDotPrev > 0.5 && dAdjDist > 2 * EPS_SMALL)
|
|
vTria[nTria].SetVertexNorm( 0, vTria[nPrevTria].GetVertexNorm( 0)) ;
|
|
}
|
|
else {
|
|
double dAdjDist = Dist( vTria[nNextTria].GetP( 1), vTria[nNextTria].GetP( 2)) ;
|
|
// Se normale sufficientemente vicina e lato opposto al vertice abbastanza lungo assegno normale al vertice del triangolo corrente
|
|
if ( dDotNext > 0.5 && dAdjDist > 2 * EPS_SMALL)
|
|
vTria[nTria].SetVertexNorm( 0, vTria[nNextTria].GetVertexNorm( 0)) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Ciclo sui voxel con sharp feature
|
|
for ( int nVox = 0 ; nVox < nVoxelNum ; ++ nVox) {
|
|
SharpTriaStruct& SharpTria = m_BlockSharpTria[nBlock][nVox] ;
|
|
// Ciclo sulle componenti connesse del voxel
|
|
int nCompNum = int( SharpTria.ptCompoVert.size()) ;
|
|
for ( int nComp = 0 ; nComp < nCompNum ; ++ nComp) {
|
|
TRIA3DEXVECTOR& vTria = SharpTria.vCompoTria[nComp] ;
|
|
// Ciclo sui triangoli della componente connessa
|
|
int nTriaNum = int( vTria.size()) ;
|
|
for ( int nTria = 0 ; nTria < nTriaNum ; ++ nTria) {
|
|
// Se il triangolo ha partecipato a un flipping
|
|
if ( vTria[nTria].GetVertexNorm( 0).IsSmall() &&
|
|
vTria[nTria].GetVertexNorm( 1).IsSmall()) {
|
|
const Point3d& ptP0 = vTria[nTria].GetP( 0) ;
|
|
const Point3d& ptP1 = vTria[nTria].GetP( 1) ;
|
|
const Point3d& ptP2 = vTria[nTria].GetP( 2) ;
|
|
int nPrevTria = ( nTria - 1 >= 0 ? nTria - 1 : nTriaNum - 1) ;
|
|
const Point3d& ptPrevP0 = vTria[nPrevTria].GetP( 0) ;
|
|
const Point3d& ptPrevP1 = vTria[nPrevTria].GetP( 1) ;
|
|
if ( AreSamePointApprox( ptP0, ptPrevP0) && AreSamePointApprox( ptP2, ptPrevP1) &&
|
|
vTria[nPrevTria].GetVertexNorm( 0).IsNormalized())
|
|
vTria[nTria].SetVertexNorm( 0, vTria[nPrevTria].GetVertexNorm( 0)) ;
|
|
else {
|
|
vTria[nTria].Validate( true) ;
|
|
vTria[nTria].SetVertexNorm( 1, V_NULL) ;
|
|
}
|
|
// Ciclo su voxel adiacenti
|
|
bool bFound = false ;
|
|
for ( int nAdjVox = 0 ; nAdjVox < nVoxelNum ; ++ nAdjVox) {
|
|
SharpTriaStruct& AdjSharpTria = m_BlockSharpTria[nBlock][nAdjVox] ;
|
|
// Se coincidenti o non adiacenti vado oltre
|
|
if ( nAdjVox == nVox ||
|
|
abs( AdjSharpTria.i - SharpTria.i) > 1 ||
|
|
abs( AdjSharpTria.j - SharpTria.j) > 1 ||
|
|
abs( AdjSharpTria.k - SharpTria.k) > 1)
|
|
continue ;
|
|
// Ciclo sulle componenti connesse del voxel adiacente
|
|
int nAdjCompNum = int( AdjSharpTria.ptCompoVert.size()) ;
|
|
for ( int nAdjComp = 0 ; nAdjComp < nAdjCompNum ; ++ nAdjComp) {
|
|
TRIA3DEXVECTOR& vAdjTria = AdjSharpTria.vCompoTria[nAdjComp] ;
|
|
// Ciclo sui triangoli della componente connessa
|
|
int nAdjTriaNum = int( vAdjTria.size()) ;
|
|
for ( int nAdjTria = 0 ; nAdjTria < nAdjTriaNum ; ++ nAdjTria) {
|
|
const Point3d& ptAdjP0 = vAdjTria[nAdjTria].GetP( 0) ;
|
|
const Point3d& ptAdjP2 = vAdjTria[nAdjTria].GetP( 2) ;
|
|
if ( AreSamePointApprox( ptP1, ptAdjP0) && AreSamePointApprox( ptP2, ptAdjP2)) {
|
|
if ( vAdjTria[nAdjTria].GetVertexNorm( 0).IsNormalized())
|
|
vTria[nTria].SetVertexNorm( 1, vAdjTria[nAdjTria].GetVertexNorm( 0)) ;
|
|
else
|
|
vTria[nTria].Validate( true) ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
}
|
|
if ( bFound)
|
|
break ;
|
|
}
|
|
if ( bFound)
|
|
break ;
|
|
}
|
|
vTria[nTria].Validate( true) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::FlipEdgesBB() const
|
|
{
|
|
// Numero di blocchi
|
|
int nBlocksNum = int( m_InterBlockSharpTria.size()) ;
|
|
// ciclo sui blocchi
|
|
for ( int tFB = 0 ; tFB < nBlocksNum ; ++ tFB) {
|
|
int nFBijk[3] ;
|
|
GetBlockIJKFromN( int( tFB), nFBijk) ;
|
|
for ( int 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
|
|
int nVoxelNumFB = int( m_InterBlockSharpTria[tFB].size()) ;
|
|
int nVoxelNumLB = int( m_InterBlockSharpTria[tLB].size()) ;
|
|
// Ciclo sui voxel dei due blocchi
|
|
for ( int tVFB = 0 ; tVFB < nVoxelNumFB ; ++ tVFB) {
|
|
for ( int tVLB = 0 ; tVLB < nVoxelNumLB ; ++ tVLB) {
|
|
// Se i voxel non sono adiacenti salto l'iterazione
|
|
if ( ! ( abs( m_InterBlockSharpTria[tFB][tVFB].i - m_InterBlockSharpTria[tLB][tVLB].i) <= 1 &&
|
|
abs( m_InterBlockSharpTria[tFB][tVFB].j - m_InterBlockSharpTria[tLB][tVLB].j) <= 1 &&
|
|
abs( m_InterBlockSharpTria[tFB][tVFB].k - m_InterBlockSharpTria[tLB][tVLB].k) <= 1))
|
|
continue ;
|
|
// Numero di componenti connesse dei voxel
|
|
int nCompoVFBNum = int( m_InterBlockSharpTria[tFB][tVFB].ptCompoVert.size()) ;
|
|
int nCompoVLBNum = int( m_InterBlockSharpTria[tLB][tVLB].ptCompoVert.size()) ;
|
|
// Ciclo sulle componenti connesse
|
|
for ( int tCmpF = 0 ; tCmpF < nCompoVFBNum ; ++ tCmpF) {
|
|
for ( int tCmpL = 0 ; tCmpL < nCompoVLBNum ; ++ tCmpL) {
|
|
// Numero di triangoli delle componenti connesse
|
|
int nTriFBNum = int( m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF].size()) ;
|
|
int nTriLBNum = int( m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL].size()) ;
|
|
// Ciclo sui triangoli
|
|
for ( int tTriFB = 0 ; tTriFB < nTriFBNum ; ++ tTriFB) {
|
|
bool bModified = false ;
|
|
for ( int tTriLB = 0 ; tTriLB < nTriLBNum ; ++ tTriLB) {
|
|
// Punti che devono essere in comune fra i due triangoli
|
|
Point3d ptPF1 = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 1) ;
|
|
Point3d ptPF2 = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 2) ;
|
|
Point3d ptPL1 = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 1) ;
|
|
Point3d ptPL2 = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 2) ;
|
|
// Si deve operare la modifica dei triangoli
|
|
if ( AreSamePointEpsilon( ptPF1, ptPL2, EPS_ZERO) &&
|
|
AreSamePointEpsilon( ptPF2, ptPL1, EPS_ZERO) &&
|
|
! ( m_InterBlockSharpTria[tFB][tVFB].vbFlipped[tCmpF][tTriFB] ||
|
|
m_InterBlockSharpTria[tLB][tVLB].vbFlipped[tCmpL][tTriLB])) {
|
|
// Assegno l'array dei punti di contorno
|
|
Point3d vPnt[4] ;
|
|
vPnt[0] = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 1) ;
|
|
vPnt[1] = m_InterBlockSharpTria[tFB][tVFB].ptCompoVert[tCmpF] ;
|
|
vPnt[2] = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 2) ;
|
|
vPnt[3] = m_InterBlockSharpTria[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) {
|
|
Triangle3dEx trTF = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB] ;
|
|
Triangle3dEx trTL = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB] ;
|
|
int nVertF, nVertL ;
|
|
// Determino gli indici dei punti sharp-feature
|
|
for ( int nP = 0 ; nP < 3 ; ++ nP) {
|
|
if ( nP != 1 && nP != 2)
|
|
nVertF = nP ;
|
|
if ( nP != 2 && nP != 1)
|
|
nVertL = nP ;
|
|
}
|
|
trTF.SetP( 1, trTL.GetP( nVertL)) ;
|
|
trTF.Validate( true) ;
|
|
trTL.SetP( 1, trTF.GetP( nVertF)) ;
|
|
trTL.Validate( true) ;
|
|
bInv = ( trTF.GetN() * trTL.GetN() < 0) ;
|
|
}
|
|
// Se non vi è inversione eseguo il flipping
|
|
if ( ! bInv) {
|
|
// Vertice condiviso fra nTri1 e quello del suo fan
|
|
int nColF = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetAttrib( 2) ;
|
|
int nColL = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetAttrib( 2) ;
|
|
// modifico punti e colori
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetP( 1, m_InterBlockSharpTria[tLB][tVLB].ptCompoVert[tCmpL]) ;
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetP( 1, m_InterBlockSharpTria[tFB][tVFB].ptCompoVert[tCmpF]) ;
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetGrade( nColF) ;
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetGrade( nColL) ;
|
|
// Valido i triangoli
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].Validate( true) ;
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].Validate( true) ;
|
|
// Setto le normali a ogni vertice
|
|
Vector3d vtNormF = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetVertexNorm( 2) ;
|
|
Vector3d vtNormL = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetVertexNorm( 2) ;
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetVertexNorm( 1, vtNormF) ;
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetVertexNorm( 0, vtNormF) ;
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetVertexNorm( 1, vtNormL) ;
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetVertexNorm( 0, vtNormL) ;
|
|
// Valido i triangoli
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].Validate( true) ;
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].Validate( true) ;
|
|
// Setto i bFlags
|
|
if ( ! AreSameVectorApprox( m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetN(),
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetN())) {
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetEdgeFlag( 0, true) ;
|
|
}
|
|
// Setto i triangoli come flippati
|
|
m_InterBlockSharpTria[tFB][tVFB].vbFlipped[tCmpF][tTriFB] = true ;
|
|
m_InterBlockSharpTria[tLB][tVLB].vbFlipped[tCmpL][tTriLB] = true ;
|
|
bModified = true ;
|
|
break ;
|
|
}
|
|
else {
|
|
// In questo caso i due triangoli sono necessariamente su un piano,
|
|
// quindi hanno normali concordi. A entrambi assegno il colore del vertice con
|
|
// normale più concorde a quella del triangolo
|
|
double dDotVec[4] ;
|
|
dDotVec[0] = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetN() *
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetVertexNorm( 1) ;
|
|
dDotVec[1] = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetN() *
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetVertexNorm( 2) ;
|
|
dDotVec[2] = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetN() *
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetVertexNorm( 2) ;
|
|
dDotVec[3] = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetN() *
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetVertexNorm( 1) ;
|
|
// Cerco il massimo dei prodotti scalari
|
|
int nMaxPos = 0 ;
|
|
double dMaxDot = - 1 ;
|
|
for ( int nPos = 0 ; nPos < 4 && dMaxDot < 1 ; ++ nPos) {
|
|
if ( dDotVec[nPos] > dMaxDot) {
|
|
dMaxDot = dDotVec[nPos] ;
|
|
nMaxPos = nPos ;
|
|
}
|
|
}
|
|
// Trovo il colore associato al vertice di massimo prodotto scalare
|
|
int nCol ;
|
|
switch ( nMaxPos) {
|
|
case 0 :
|
|
nCol = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetAttrib( 1) ;
|
|
break ;
|
|
case 1 :
|
|
nCol = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetAttrib( 2) ;
|
|
break ;
|
|
case 2 :
|
|
nCol = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetAttrib( 2) ;
|
|
break ;
|
|
case 3 :
|
|
nCol = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetAttrib( 1) ;
|
|
break ;
|
|
}
|
|
// Assegno il colore ai triangoli
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetGrade( nCol) ;
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetGrade( nCol) ;
|
|
// Setto i flag dei vertici
|
|
double dDist02 = ( m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 0) -
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 2)).Len() +
|
|
( m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 1) -
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 0)).Len() ;
|
|
double dDist20 = ( m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 1) -
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 0)).Len() +
|
|
( m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 0) -
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 2)).Len() ;
|
|
if ( dDist02 > dDist20) {
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetEdgeFlag( 0, true) ;
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetEdgeFlag( 2, true) ;
|
|
}
|
|
else {
|
|
m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetEdgeFlag( 2, true) ;
|
|
m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetEdgeFlag( 0, true) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if ( bModified)
|
|
break ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IsThereMat( int nI, int nJ, int nK) const
|
|
{
|
|
// Trasformo gli 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
|
|
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
|
|
int nIndex = 0 ;
|
|
int nPos = nGrJ * m_nNx[nGrid] + nGrI ;
|
|
int nDexSize = int( 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, AppliedVector& vfField) const
|
|
{
|
|
const double dEps = 2 * EPS_SMALL ;
|
|
const double dEpsClamp = EPS_SMALL ; // EPS_SMALL ; 0
|
|
|
|
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 ;
|
|
|
|
vfField.ptPApp.y = ( nVec1[1] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
vfField.ptPApp.z = ( nVec1[2] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
|
|
int 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) {
|
|
vfField.ptPApp.x = Clamp( dX, dMinX + dEpsClamp, dMaxX - dEpsClamp) ;
|
|
vfField.vtVec = m_Values[1][nDexel][n].vtMaxN ;
|
|
vfField.nPropIndex = m_Values[1][nDexel][n].nToolMax ;
|
|
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) {
|
|
vfField.ptPApp.x = Clamp( dX, dMinX + dEpsClamp, dMaxX - dEpsClamp) ;
|
|
vfField.vtVec = m_Values[1][nDexel][n].vtMinN ;
|
|
vfField.nPropIndex = m_Values[1][nDexel][n].nToolMin ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
n += 1 ;
|
|
if ( n < nSize)
|
|
dX = m_Values[1][nDexel][n].dMin ;
|
|
}
|
|
}
|
|
|
|
if ( ! bFound)
|
|
vfField.ptPApp.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 ;
|
|
|
|
vfField.ptPApp.x = ( nVec1[0] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
vfField.ptPApp.z = ( nVec1[2] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
|
|
int 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) {
|
|
vfField.ptPApp.y = Clamp( dY, dMinY + dEpsClamp, dMaxY - dEpsClamp) ;
|
|
vfField.vtVec = m_Values[2][nDexel][n].vtMaxN ;
|
|
vfField.nPropIndex = m_Values[2][nDexel][n].nToolMax;
|
|
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) {
|
|
vfField.ptPApp.y = Clamp( dY, dMinY + dEpsClamp, dMaxY - dEpsClamp) ;
|
|
vfField.vtVec = m_Values[2][nDexel][n].vtMinN ;
|
|
vfField.nPropIndex = m_Values[2][nDexel][n].nToolMin ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
n += 1 ;
|
|
if ( n < nSize)
|
|
dY = m_Values[2][nDexel][n].dMin ;
|
|
}
|
|
}
|
|
|
|
if ( ! bFound)
|
|
vfField.ptPApp.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 ;
|
|
|
|
vfField.ptPApp.x = ( nVec1[0] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
vfField.ptPApp.y = ( nVec1[1] * N_DEXVOXRATIO + 0.5) * m_dStep ;
|
|
|
|
int 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) {
|
|
vfField.ptPApp.z = Clamp( dZ, dMinZ + dEpsClamp, dMaxZ - dEpsClamp) ;
|
|
vfField.vtVec = m_Values[0][nDexel][n].vtMaxN ;
|
|
vfField.nPropIndex = m_Values[0][nDexel][n].nToolMax ;
|
|
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) {
|
|
vfField.ptPApp.z = Clamp( dZ, dMinZ + dEpsClamp, dMaxZ - dEpsClamp) ;
|
|
vfField.vtVec = m_Values[0][nDexel][n].vtMinN ;
|
|
vfField.nPropIndex = m_Values[0][nDexel][n].nToolMin ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
n += 1 ;
|
|
if ( n < nSize)
|
|
dZ = m_Values[0][nDexel][n].dMin ;
|
|
}
|
|
}
|
|
|
|
if ( ! bFound)
|
|
vfField.ptPApp.z = 0.5 * ( dMinZ + dMaxZ) ;
|
|
}
|
|
|
|
return bFound ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IsValidVoxel( int nN) const
|
|
{
|
|
// Calcolo il numero di voxel lungo X,Y e Z
|
|
int nVoxNumX = m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2) ;
|
|
int nVoxNumY = m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] % N_DEXVOXRATIO == 0 ? 1 : 2) ;
|
|
int nVoxNumZ = m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2) ;
|
|
// Verifico la validità del voxel
|
|
return ( nN >= 0 && nN < ( 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 - 2 &&
|
|
nJ >= -1 && nJ <= nVoxNumY - 2 &&
|
|
nK >= -1 && nK <= nVoxNumZ - 2) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
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 - 2 &&
|
|
nJ >= -1 && nJ <= nVoxNumY - 2 &&
|
|
nK >= -1 && nK <= nVoxNumZ - 2) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
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 < -1 || nI > nVoxNumX - 2 ||
|
|
nJ < -1 || nJ > nVoxNumY - 2 ||
|
|
nK < -1 || nK > nVoxNumZ - 2)
|
|
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::GetVoxelBox( int nI, int nJ, int nK, BBox3d& b3VoxBox) 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 < -1 || nI > nVoxNumX - 2 ||
|
|
nJ < -1 || nJ > nVoxNumY - 2 ||
|
|
nK < -1 || nK > nVoxNumZ - 2)
|
|
return false ;
|
|
|
|
// Punti della diagonale del voxel
|
|
Point3d ptCubeInf( ( nI * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nJ * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nK * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
Point3d ptCubeSup( ( ( nI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( nJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( ( nK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
b3VoxBox.Set( ptCubeInf, ptCubeSup) ;
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
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] >= m_nFracLin[0] ||
|
|
nIJK[1] < 0 || nIJK[1] >= m_nFracLin[1] ||
|
|
nIJK[2] < 0 || nIJK[2] >= 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] * m_nVoxNumPerBlock) ;
|
|
nLimits[1] = ( nIJK[0] + 1 == m_nFracLin[0] ?
|
|
nVoxNumX - 1 : ( nIJK[0] + 1) * m_nVoxNumPerBlock) ;
|
|
|
|
// Calcolo limiti per l'indice j
|
|
nLimits[2] = ( nIJK[1] == 0 ? - 1 : nIJK[1] * m_nVoxNumPerBlock) ;
|
|
nLimits[3] = ( nIJK[1] + 1 == m_nFracLin[1] ?
|
|
nVoxNumY - 1 : ( nIJK[1] + 1) * m_nVoxNumPerBlock) ;
|
|
|
|
// Calcolo limiti per l'indice k
|
|
nLimits[4] = ( nIJK[2] == 0 ? - 1 : nIJK[2] * m_nVoxNumPerBlock) ;
|
|
nLimits[5] = ( nIJK[2] + 1 == m_nFracLin[2] ?
|
|
nVoxNumZ - 1 : ( nIJK[2] + 1) * m_nVoxNumPerBlock) ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetBlockBox( const int nIJK[], BBox3d& b3VoxBox) const
|
|
{
|
|
// Calcolo limiti sugli indici del blocco; se il blocco non è valido, errore
|
|
int nLimits[6] ;
|
|
if ( ! GetBlockLimitsIJK( nIJK, nLimits))
|
|
return false ;
|
|
// Costruisco il bounding-bx del blocco
|
|
Point3d ptMinDiag( ( nLimits[0] * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nLimits[2] * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nLimits[4] * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
Point3d ptMaxDiag( ( nLimits[1] * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nLimits[3] * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nLimits[5] * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
b3VoxBox.Set( ptMinDiag, ptMaxDiag) ;
|
|
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 < -1 || nI > nVoxNumX - 2 ||
|
|
nJ < -1 || nJ > nVoxNumY - 2 ||
|
|
nK < -1 || nK > nVoxNumZ - 2)
|
|
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 - 2 &&
|
|
nVoxJ >= - 1 && nVoxJ <= nVoxNumY - 2 &&
|
|
nVoxK >= - 1 && nVoxK <= nVoxNumZ - 2) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
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] < -1 || nVoxIJK[0] > nVoxNumX - 2 ||
|
|
nVoxIJK[1] < -1 || nVoxIJK[1] > nVoxNumY - 2 ||
|
|
nVoxIJK[2] < -1 || nVoxIJK[2] > nVoxNumZ - 2)
|
|
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::GetVoxelBlockN( int nVox, int& nBlock) const
|
|
{
|
|
int nVoxIJK[3] ;
|
|
int nBlockIJK[3] ;
|
|
return ( GetVoxIJKFromN( nVox, nVoxIJK[0], nVoxIJK[1], nVoxIJK[2]) &&
|
|
GetVoxelBlockIJK( nVoxIJK, nBlockIJK) &&
|
|
GetBlockNFromIJK( nBlockIJK, nBlock)) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
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] < -1 || nIJK[0] > nVoxNumX - 2 ||
|
|
nIJK[1] < -1 || nIJK[1] > nVoxNumY - 2 ||
|
|
nIJK[2] < -1 || 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] < m_nFracLin[0] &&
|
|
nAdjBlockIJK[1] > -1 && nAdjBlockIJK[1] < m_nFracLin[1] &&
|
|
nAdjBlockIJK[2] > -1 && nAdjBlockIJK[2] < m_nFracLin[2] &&
|
|
( nCurrentBlockIJK[0] != nAdjBlockIJK[0] ||
|
|
nCurrentBlockIJK[1] != nAdjBlockIJK[1] ||
|
|
nCurrentBlockIJK[2] != nAdjBlockIJK[2])) {
|
|
return true ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return false ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], int nDeltaIndex[]) 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] < -1 || nIJK[0] > nVoxNumX - 2 ||
|
|
nIJK[1] < -1 || nIJK[1] > nVoxNumY - 2 ||
|
|
nIJK[2] < -1 || nIJK[2] > nVoxNumZ - 2)
|
|
return false ;
|
|
nDeltaIndex[0] = 0 ;
|
|
nDeltaIndex[1] = 0 ;
|
|
nDeltaIndex[2] = 0 ;
|
|
if ( nIJK[0] == nLimits[0])
|
|
-- nDeltaIndex[0] ;
|
|
else if ( nIJK[0] == nLimits[1] - 1)
|
|
++ nDeltaIndex[0] ;
|
|
if ( nIJK[0] == nLimits[0])
|
|
-- nDeltaIndex[1];
|
|
else if ( nIJK[0] == nLimits[1] - 1)
|
|
++ nDeltaIndex[1];
|
|
if ( nIJK[0] == nLimits[0])
|
|
-- nDeltaIndex[2];
|
|
else if ( nIJK[0] == nLimits[1] - 1)
|
|
++ nDeltaIndex[2];
|
|
return ( nDeltaIndex[0] != 0 || nDeltaIndex[1] != 0 || nDeltaIndex[2] != 0) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Per riferimento restituisce i minimi degli indici ijk riferiti ai Voxel.
|
|
// Si assume che il volume sia un parallelepipedo.
|
|
bool
|
|
VolZmap::GetFirstVoxIJK( int& i, int& j, int& k) const
|
|
{
|
|
// Se non è un parallelepipedo, errore
|
|
if ( m_nShape != BOX)
|
|
return false ;
|
|
// Calcolo indici
|
|
int ni, nj ;
|
|
for ( ni = 0 ; ni < 2 ; ++ ni) {
|
|
bool bNotEmpty = false ;
|
|
for ( nj = 0 ; nj < 2 ; ++ nj) {
|
|
int nDex = nj * int( m_nNx[0]) + ni ;
|
|
if ( m_Values[0][nDex].size() > 0) {
|
|
bNotEmpty = true ;
|
|
break ;
|
|
}
|
|
}
|
|
if ( bNotEmpty)
|
|
break ;
|
|
}
|
|
int mi, mj ;
|
|
for ( mi = 0 ; mi < 2 ; ++ mi) {
|
|
bool bNotEmpty = false ;
|
|
for ( mj = 0 ; mj < 2 ; ++ mj) {
|
|
int nDex = mj * int( m_nNx[1]) + mi ;
|
|
if ( m_Values[1][nDex].size() > 0) {
|
|
bNotEmpty = true ;
|
|
break ;
|
|
}
|
|
}
|
|
if ( bNotEmpty)
|
|
break ;
|
|
}
|
|
i = ni / N_DEXVOXRATIO - ( ni % N_DEXVOXRATIO) - 1 ;
|
|
j = nj / N_DEXVOXRATIO - ( nj % N_DEXVOXRATIO) - 1 ;
|
|
k = mj / N_DEXVOXRATIO - ( mj % N_DEXVOXRATIO) - 1 ;
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Per riferimento restituisce i massimi degli indici ijk riferiti ai Voxel.
|
|
// Si assume che il volume sia un parallelepipedo.
|
|
bool
|
|
VolZmap::GetLastVoxIJK( int& i, int& j, int& k) const
|
|
{
|
|
// Se non è un parallelepipedo, errore
|
|
if ( m_nShape != BOX)
|
|
return false ;
|
|
// Calcolo indici
|
|
int ni, nj ;
|
|
for ( ni = int( m_nNx[0]) - 1 ; ni > int( m_nNx[0]) - 3 ; -- ni) {
|
|
bool bNotEmpty = false ;
|
|
for ( nj = int( m_nNy[0]) - 1 ; nj > int( m_nNy[0]) - 3 ; -- nj) {
|
|
int nDex = nj * int( m_nNx[0]) + ni ;
|
|
if ( m_Values[0][nDex].size() > 0) {
|
|
bNotEmpty = true ;
|
|
break ;
|
|
}
|
|
}
|
|
if ( bNotEmpty)
|
|
break ;
|
|
}
|
|
int mi, mj ;
|
|
for ( mi = int( m_nNx[1]) - 1 ; mi > int( m_nNx[1]) - 3 ; -- mi) {
|
|
bool bNotEmpty = false ;
|
|
for ( mj = int( m_nNy[1]) - 1 ; mj > int( m_nNy[1]) - 3 ; -- mj) {
|
|
int nDex = mj * int( m_nNx[1]) + mi ;
|
|
if ( m_Values[1][nDex].size() > 0) {
|
|
bNotEmpty = true ;
|
|
break ;
|
|
}
|
|
}
|
|
if ( bNotEmpty)
|
|
break ;
|
|
}
|
|
i = ni / N_DEXVOXRATIO ;
|
|
j = nj / N_DEXVOXRATIO ;
|
|
k = mj / N_DEXVOXRATIO ;
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Il volume deve essere un parallelepipedo.
|
|
// Verifica se il voxel contiene uno dei dodici spigoli del parallelepipedo.
|
|
bool
|
|
VolZmap::IsVoxelOnBoxEdge( int i, int j, int k) const
|
|
{
|
|
// Se non è un parallelepipedo, errore
|
|
if ( m_nShape != BOX)
|
|
return false ;
|
|
// Determino il primo nodo pieno della mappa
|
|
int nFirstVoxI, nFirstVoxJ, nFirstVoxK ;
|
|
GetFirstVoxIJK( nFirstVoxI, nFirstVoxJ, nFirstVoxK) ;
|
|
// Determino il primo nodo pieno della mappa
|
|
int nLastVoxI, nLastVoxJ, nLastVoxK ;
|
|
GetLastVoxIJK( nLastVoxI, nLastVoxJ, nLastVoxK) ;
|
|
|
|
// Determino se il voxel è su un edge del box
|
|
int nIndexOnLimitNum = 0 ;
|
|
if ( i == nFirstVoxI || i == nLastVoxI)
|
|
++ nIndexOnLimitNum ;
|
|
if ( j == nFirstVoxJ || j == nLastVoxJ)
|
|
++ nIndexOnLimitNum ;
|
|
if ( k == nFirstVoxK || k == nLastVoxK)
|
|
++ nIndexOnLimitNum ;
|
|
return ( nIndexOnLimitNum >= 2) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::IsTriangleOnBorder( const Triangle3dEx& trTria, const int nBlockLimits[], const int nVoxIJK[]) const
|
|
{
|
|
// Punti del triangolo sulla griglia
|
|
Point3d ptFirstGrPt = trTria.GetP( 1) ;
|
|
Point3d ptSecondGrPt = trTria.GetP( 2) ;
|
|
|
|
// 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::CheckForFanNodeInterferance( AppliedVector BaseVecField[], Point3d& ptFanVert,
|
|
int nBasePointNum, int nVoxConfig, int nVoxI, int nVoxJ, int nVoxK) const
|
|
{
|
|
bool bInterferance = false ;
|
|
int nFilterIndex = 1 ;
|
|
for ( int n = 0 ; n < 7 ; ++ n) {
|
|
bool bNodeInside = true ;
|
|
if ( ( nVoxConfig & nFilterIndex) != 0) {
|
|
int nNodeI = nVoxI ;
|
|
int nNodeJ = nVoxJ ;
|
|
int nNodeK = nVoxK ;
|
|
if ( ( n + 1) % 4 >= 2)
|
|
++ nNodeI ;
|
|
if ( n % 4 >= 2)
|
|
++ nNodeJ ;
|
|
if ( n / 4 == 1)
|
|
++ nNodeK ;
|
|
Point3d ptNodeP( ( nNodeI * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nNodeJ * N_DEXVOXRATIO + 0.5) * m_dStep,
|
|
( nNodeK * N_DEXVOXRATIO + 0.5) * m_dStep) ;
|
|
for ( int nV = 0 ; nV < nBasePointNum ; ++ nV) {
|
|
int nW = ( nV + 1) % nBasePointNum ;
|
|
Plane3d plPlane ;
|
|
Vector3d vtN = ( BaseVecField[nW].ptPApp - ptFanVert) ^ ( BaseVecField[nV].ptPApp - ptFanVert) ;
|
|
if ( plPlane.Set( ptFanVert, vtN) && DistPointPlane( ptNodeP, plPlane) > -EPS_SMALL) {
|
|
bNodeInside = false ;
|
|
break ;
|
|
}
|
|
}
|
|
if ( bNodeInside) {
|
|
for ( int nV = 2 ; nV < nBasePointNum ; ++ nV) {
|
|
Plane3d plPlane ;
|
|
Vector3d vtN = ( BaseVecField[nV - 1].ptPApp - BaseVecField[0].ptPApp) ^ ( BaseVecField[nV].ptPApp - BaseVecField[0].ptPApp) ;
|
|
if ( plPlane.Set( BaseVecField[0].ptPApp, vtN) && DistPointPlane( ptNodeP, plPlane) > -EPS_SMALL) {
|
|
bNodeInside = false ;
|
|
break ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else
|
|
bNodeInside = false ;
|
|
|
|
bInterferance = bInterferance || bNodeInside ;
|
|
|
|
nFilterIndex *= 2 ;
|
|
}
|
|
|
|
return bInterferance ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::ProcessVoxContXY( FlatVoxelContainer& VoxContXY, int nBlock, bool bPlus) const
|
|
{
|
|
for ( auto it = VoxContXY.cbegin() ;
|
|
it != VoxContXY.end() ;
|
|
it = VoxContXY.begin()) {
|
|
|
|
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 ;
|
|
int nToolNum = it->second.nTool ;
|
|
double dCordZ = it->second.dHeigth ;
|
|
|
|
// 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, nToolNum))
|
|
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, nToolNum))
|
|
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, nToolNum))
|
|
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, nToolNum))
|
|
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) ;
|
|
|
|
Triangle3dEx 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) ;
|
|
}
|
|
Tria0.SetGrade( nToolNum) ;
|
|
Tria1.SetGrade( nToolNum) ;
|
|
bool bV0 = Tria0.Validate( true) ;
|
|
bool bV1 = Tria1.Validate( true) ;
|
|
// Aggiungo alla lista
|
|
m_BlockBigTria[nBlock].emplace_back( Tria0) ;
|
|
m_BlockBigTria[nBlock].emplace_back( Tria1) ;
|
|
|
|
// Elimino il voxel da cui sono partito a ingrandire.
|
|
VoxContXY.erase( nN) ;
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::ProcessVoxContYZ( FlatVoxelContainer& VoxContYZ, int nBlock, bool bPlus) const
|
|
{
|
|
for ( auto it = VoxContYZ.begin() ;
|
|
it != VoxContYZ.end() ;
|
|
it = VoxContYZ.begin()) {
|
|
|
|
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 ;
|
|
int nToolNum = it->second.nTool ;
|
|
double dCordX = it->second.dHeigth ;
|
|
|
|
// 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, nToolNum))
|
|
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, nToolNum))
|
|
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, nToolNum))
|
|
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, nToolNum))
|
|
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) ;
|
|
|
|
Triangle3dEx 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) ;
|
|
}
|
|
Tria0.SetGrade( nToolNum) ;
|
|
Tria1.SetGrade( nToolNum) ;
|
|
bool bV0 = Tria0.Validate( true) ;
|
|
bool bV1 = Tria1.Validate( true) ;
|
|
// Aggiungo alla lista
|
|
m_BlockBigTria[nBlock].emplace_back( Tria0) ;
|
|
m_BlockBigTria[nBlock].emplace_back( Tria1) ;
|
|
|
|
// Elimino il voxel da cui sono partito a ingrandire.
|
|
VoxContYZ.erase( nN) ;
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::ProcessVoxContXZ( FlatVoxelContainer& VoxContXZ, int nBlock, bool bPlus) const
|
|
{
|
|
for ( auto it = VoxContXZ.begin() ;
|
|
it != VoxContXZ.end() ;
|
|
it = VoxContXZ.begin()) {
|
|
|
|
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 ;
|
|
int nToolNum = it->second.nTool ;
|
|
double dCordY = it->second.dHeigth ;
|
|
|
|
// 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, nToolNum))
|
|
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, nToolNum))
|
|
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, nToolNum))
|
|
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, nToolNum))
|
|
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) ;
|
|
|
|
Triangle3dEx 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) ;
|
|
}
|
|
Tria0.SetGrade( nToolNum) ;
|
|
Tria1.SetGrade( nToolNum) ;
|
|
bool bV0 = Tria0.Validate( true) ;
|
|
bool bV1 = Tria1.Validate( true) ;
|
|
// Aggiungo alla lista
|
|
m_BlockBigTria[nBlock].emplace_back( Tria0) ;
|
|
m_BlockBigTria[nBlock].emplace_back( Tria1) ;
|
|
|
|
// Elimino il voxel da cui sono partito a ingrandire.
|
|
VoxContXZ.erase( nN) ;
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::Find( const FlatVoxelContainer& VoxCont, int nI, int nJ, int nK, double dPos, int nTool) 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.dHeigth) < EPS_SMALL &&
|
|
nTool == iter->second.nTool) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::Remove( FlatVoxelContainer& VoxCont, int nI, int nJ, int nK) const
|
|
{
|
|
int nN ;
|
|
if ( ! GetVoxNFromIJK( nI, nJ, nK, nN))
|
|
return false ;
|
|
return ( VoxCont.erase( nN) > 0) ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetEdges( ICURVEPOVECTOR& vpCurve) const
|
|
{
|
|
// Garantisco grafica aggiornata
|
|
UpdateTripleMapGraphics() ;
|
|
|
|
// Vettore dei segmenti di feature (coppie di punti)
|
|
BIPNTVECTOR vBpt ;
|
|
|
|
// Ciclo sui triangoli feature all'interno dei blocchi
|
|
for ( int nBlock = 0 ; nBlock < m_nNumBlock ; ++ nBlock) {
|
|
// Numero di voxel in cui si presentano sharp feature
|
|
int nVoxelNum = int( m_BlockSharpTria[nBlock].size()) ;
|
|
// Ciclo sui voxel con sharp feature
|
|
for ( int n = 0 ; n < nVoxelNum ; ++ n) {
|
|
const SharpTriaStruct& SharpTria = m_BlockSharpTria[nBlock][n] ;
|
|
// Ciclo sulle componenti connesse del voxel
|
|
int nNumCompo = int( SharpTria.ptCompoVert.size()) ;
|
|
for ( int nCompo = 0 ; nCompo < nNumCompo ; ++ nCompo) {
|
|
const TRIA3DEXVECTOR& vTria = SharpTria.vCompoTria[nCompo] ;
|
|
// Numero di triangoli della componente connessa
|
|
int nTriNum = int( vTria.size()) ;
|
|
for ( int nTri = 0 ; nTri < nTriNum ; ++ nTri) {
|
|
// Il triangolo ha un lato sulla sharp-feature
|
|
if ( vTria[nTri].GetEdgeFlag( 0)) {
|
|
const Point3d& ptPS = vTria[nTri].GetP( 0) ;
|
|
const Point3d& ptPE = vTria[nTri].GetP( 1) ;
|
|
// Segmento sharp-feature
|
|
vBpt.emplace_back( ptPS, ptPE) ;
|
|
}
|
|
else if ( vTria[nTri].GetEdgeFlag( 1)) {
|
|
const Point3d& ptPS = vTria[nTri].GetP( 1) ;
|
|
const Point3d& ptPE = vTria[nTri].GetP( 2) ;
|
|
// Segmento sharp-feature
|
|
vBpt.emplace_back( ptPS, ptPE) ;
|
|
}
|
|
else if ( vTria[nTri].GetEdgeFlag( 2)) {
|
|
const Point3d& ptPS = vTria[nTri].GetP( 2) ;
|
|
const Point3d& ptPE = vTria[nTri].GetP( 0) ;
|
|
// Segmento sharp-feature
|
|
vBpt.emplace_back( ptPS, ptPE) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Ciclo sui triangoli feature al confine tra blocchi
|
|
for ( int tB = 0 ; tB < m_nNumBlock ; ++ tB) {
|
|
// Numero di voxel nel blocco corrente
|
|
int nVoxelNum = int( m_InterBlockSharpTria[tB].size()) ;
|
|
// Ciclo sui voxel del blocco
|
|
for ( int tV = 0 ; tV < nVoxelNum ; ++ tV) {
|
|
// Numero di componenti connesse del voxel
|
|
int nCompoNum = int( m_InterBlockSharpTria[tB][tV].ptCompoVert.size()) ;
|
|
// Ciclo sulle componenti connesse
|
|
for ( int tC = 0 ; tC < nCompoNum ; ++ tC) {
|
|
// Numero di triangoli della componente connessa
|
|
int nTriNum = int( m_InterBlockSharpTria[tB][tV].vCompoTria[tC].size()) ;
|
|
// Ciclo sui triangoli
|
|
for ( int tT = 0 ; tT < nTriNum ; ++ tT) {
|
|
// Il triangolo ha un lato sulla sharp-feature
|
|
if ( m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetEdgeFlag( 0)) {
|
|
const Point3d& ptPS = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 0) ;
|
|
const Point3d& ptPE = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 1) ;
|
|
// Segmento sharp-feature
|
|
vBpt.emplace_back( ptPS, ptPE) ;
|
|
}
|
|
else if ( m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetEdgeFlag( 1)) {
|
|
const Point3d& ptPS = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 1) ;
|
|
const Point3d& ptPE = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 2) ;
|
|
// Segmento sharp-feature
|
|
vBpt.emplace_back( ptPS, ptPE) ;
|
|
}
|
|
else if ( m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetEdgeFlag( 2)) {
|
|
const Point3d& ptPS = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 2) ;
|
|
const Point3d& ptPE = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 0) ;
|
|
// Segmento sharp-feature
|
|
vBpt.emplace_back( ptPS, ptPE) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Concateno le curve rispettando le biforcazioni
|
|
const double dToler = 20 * EPS_SMALL ;
|
|
ChainCurves chainC ;
|
|
chainC.Init( true, dToler, int( vBpt.size())) ;
|
|
for ( int i = 0 ; i < int( vBpt.size()) ; ++ i) {
|
|
Vector3d vtDir = vBpt[i].second - vBpt[i].first ;
|
|
vtDir.Normalize() ;
|
|
if ( ! chainC.AddCurve( i + 1, vBpt[i].first, vtDir, vBpt[i].second, vtDir))
|
|
return false ;
|
|
}
|
|
// recupero i percorsi concatenati
|
|
Point3d ptNear = ( vBpt.empty() ? ORIG : vBpt[0].first) ;
|
|
INTVECTOR vId ;
|
|
while ( chainC.GetChainFromNear( ptNear, true, vId)) {
|
|
// creo una curva composita
|
|
PtrOwner<ICurveComposite> pCrvCompo( CreateCurveComposite()) ;
|
|
if ( IsNull( pCrvCompo))
|
|
return false ;
|
|
// recupero gli estremi dei segmenti, creo le linee e le inserisco nella composita
|
|
for ( int i = 0 ; i < int( vId.size()) ; ++ i) {
|
|
// creo un segmento di retta
|
|
int nInd = abs( vId[i]) - 1 ;
|
|
bool bInvert = ( vId[i] < 0) ;
|
|
PtrOwner<ICurveLine> pLine( CreateCurveLine()) ;
|
|
if ( IsNull( pLine) || ! pLine->Set( vBpt[nInd].first, vBpt[nInd].second))
|
|
return false ;
|
|
if ( bInvert)
|
|
pLine->Invert() ;
|
|
// lo accodo alla composita
|
|
if ( ! pCrvCompo->AddCurve( Release( pLine), true, 1.1 * dToler) &&
|
|
! AreSamePointApprox( ptNear, ( bInvert ? vBpt[nInd].first : vBpt[nInd].second)))
|
|
return false ;
|
|
// aggiorno il prossimo near
|
|
ptNear = ( bInvert ? vBpt[nInd].first : vBpt[nInd].second) ;
|
|
}
|
|
// se lunghezza curva inferiore a 5 volte la tolleranza, la salto
|
|
double dCrvLen ;
|
|
if ( ! pCrvCompo->GetLength( dCrvLen) || dCrvLen < 5 * dToler)
|
|
continue ;
|
|
// unisco segmenti allineati
|
|
pCrvCompo->MergeCurves( dToler, ANG_TOL_STD_DEG) ;
|
|
// la salvo
|
|
vpCurve.emplace_back( Release( pCrvCompo)) ;
|
|
}
|
|
|
|
return true ;
|
|
}
|