Egalware/EgtGeomKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
76d94d6194cea232671ab7ca16934339b0127610
EgtGeomKernel/VolTriZmapGraphics.cpp
T
Dario Sassi 76d94d6194 EgtGeomKernel 1.8g1 : 
- modifiche a Zmap
- aggiunto clamp a ratio di font
- migliorata ricerca nomi con *.
2017-08-01 07:59:56 +00:00

2800 lines
116 KiB
C++
Raw Blame History

//----------------------------------------------------------------------------
// EgalTech 2015-2016
//----------------------------------------------------------------------------
// File : VolZmap.cpp Data : 22.01.15 Versione : 1.6a4
// Contenuto : Implementazione della classe Volume Zmap (tre griglie)
//
//
//
// Modifiche : 22.01.15 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------
//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "CurveLine.h"
#include "VolZmap.h"
#include "GeoConst.h"
#include "IntersLineSurfTm.h"
#include "MC_Tables.h"
#include "/EgtDev/Include/EGkIntervals.h"
#include "/EgtDev/Include/EgtNumUtils.h"
#include "/EgtDev/Include/EGkStringUtils3d.h"
#include "/EgtDev/Extern/Eigen\Core"
#include "/EgtDev/Extern/Eigen\SVD"
using namespace std ;
// ------------------------- STRUTTURA VERTICE TRIANGOLO - NORMALE ALLA SUPERFICIE ------------------------------------------------
struct VectorField {
Point3d ptInt ;
Vector3d vtNorm ;
} ;
// ------------------------- 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} ;
//----------------------------------------------------------------------------
int
TestOnNormal( const VectorField CompoVert[], int nCompoElem)
{
// Cerco la massima deviazione tra le normali nei punti della parte connessa
int nI, nJ ;
double dMinCosTheta = 2 ;
for ( int i = 0 ; i < nCompoElem ; ++ i) {
for ( int j = i + 1 ; j < nCompoElem ; ++ j) {
double dCurrCos = CompoVert[i].vtNorm * CompoVert[j].vtNorm ;
if ( dCurrCos < dMinCosTheta) {
dMinCosTheta = dCurrCos ;
nI = i ;
nJ = j ;
}
}
}
// Se la massima deviazione non supera il limite non è feature
const double SHARP_COS = 0.9 ; // 0.8 ;
if ( dMinCosTheta >= SHARP_COS)
return NO_FEATURE ;
// Verifico se Edge o Corner
// direzione perpendicolare alle normali con massima differenza (potenziale edge)
Vector3d vtK = CompoVert[nI].vtNorm ^ CompoVert[nJ].vtNorm ;
// cerco normale con massima vicinanza al potenziale edge
double dMaxAbsCos = 0 ;
for ( int i = 0 ; i < nCompoElem ; ++ i) {
double dAbsCurrentCos = abs( CompoVert[i].vtNorm * vtK) ;
if ( dAbsCurrentCos > dMaxAbsCos)
dMaxAbsCos = dAbsCurrentCos ;
}
// se esiste normale diretta quasi come potenziale edge, allora corner
const double CORNER_COS = 0.7 ; // 0.7 ; // 0.5 ;
if ( dMaxAbsCos > CORNER_COS)
return CORNER ;
else
return EDGE ;
}
//----------------------------------------------------------------------------
bool
DotTest( const VectorField CompoVert[], int nCompoElem, Vector3d& vtAvg, double dThreshold = 0)
{
// Cerco la massima deviazione tra le normali nei punti della parte connessa
double dMinCosTheta = 2 ;
for ( int i = 0 ; i < nCompoElem ; ++ i) {
for ( int j = i + 1 ; j < nCompoElem ; ++ j) {
double dCurrCos = CompoVert[i].vtNorm * CompoVert[j].vtNorm ;
if ( dCurrCos < dMinCosTheta) {
dMinCosTheta = dCurrCos ;
}
}
}
// se normali sparpagliate oltre limite
if ( dMinCosTheta < dThreshold)
return false ;
// determino media delle normali
vtAvg = V_NULL ;
for ( int i = 0 ; i < nCompoElem ; ++ i)
vtAvg += CompoVert[i].vtNorm ;
vtAvg /= nCompoElem ;
return true ;
}
// ------------------------- 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 ;
// 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 = m_MapFrame[nDir].Orig() + dX * m_MapFrame[nDir].VersX() + dY * m_MapFrame[nDir].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 * m_MapFrame[nDir].VersZ()) ;
lstPL.back().AddUPoint( 1, ptP + m_Values[nDir][nPos][j].dMax * m_MapFrame[nDir].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 ;
// 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 = m_MapFrame[nDir].Orig() + dX * m_MapFrame[nDir].VersX() + dY * m_MapFrame[nDir].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 * m_MapFrame[nDir].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[0]) ;
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 * m_MapFrame[nDir].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[0]) ;
lstPL.back().AddUPoint( 1, ptR + vtW * m_dStep / 4) ;
}
return true ;
}
}
//----------------------------------------------------------------------------
bool
VolZmap::GetAllTriangles( TRIA3DLIST& lstTria) const
{
INTVECTOR nModifiedBlocks ;
TRIA3DLISTVECTOR vLstTria ;
if ( ! GetTriangles( true, nModifiedBlocks, vLstTria))
return false ;
lstTria.clear() ;
for ( size_t i = 0 ; i < vLstTria.size() ; ++ i) {
lstTria.splice( lstTria.end(), vLstTria[i]) ;
}
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVECTOR& vLstTria) const
{
// Se nessun blocco modificato, è richiesta esterna e li considero tutti modificati
bool bSomeModif = false ;
for ( size_t i = 0 ; i < m_nNumBlock ; ++ i) {
if ( m_BlockToUpdate[i]) {
bSomeModif = true ;
break ;
}
}
if ( ! bSomeModif)
bAllBlocks = true ;
// Caso di singola mappa
if ( m_nMapNum == 1) {
const int MAX_DIM_CHUNK = 128 ;
nModifiedBlocks.resize( m_nNumBlock) ;
vLstTria.reserve( m_nNumBlock) ;
// Ciclo sui blocchi
for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) {
// Se il blocco deve essere aggiornato, eseguo
if ( bAllBlocks || m_BlockToUpdate[t]) {
// preparo lista
vLstTria.emplace_back() ;
nModifiedBlocks[t] = int( vLstTria.size()) - 1 ;
// Calcolo posizione del blocco nella griglia
int nIBlock = int( t) % int( m_nFracLin[0]) ;
int nJBlock = int( t) / int( m_nFracLin[0]) ;
// Calcolo limiti per l'indice i
int nStartI = nIBlock * int( m_nDexNumPBlock) ;
int nEndI = ( nIBlock + 1 == int( m_nFracLin[0]) ?
int( m_nNx[0]) : ( nIBlock + 1) * int( m_nDexNumPBlock)) ;
// Calcolo limiti per l'indice j
int nStartJ = nJBlock * int( m_nDexNumPBlock) ;
int nEndJ = ( nJBlock + 1 == int( m_nFracLin[1]) ?
int( m_nNy[0]) : ( nJBlock + 1) * int( m_nDexNumPBlock)) ;
// Ciclo su i e j
for ( int i = nStartI ; i < nEndI ; i += MAX_DIM_CHUNK) {
int nDimChunkX = min( MAX_DIM_CHUNK, nEndI - i) ;
for ( int j = nStartJ ; j < nEndJ ; j += MAX_DIM_CHUNK) {
int nDimChunkY = min( MAX_DIM_CHUNK, nEndJ - j) ;
GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, MAX_DIM_CHUNK, vLstTria.back()) ;
}
}
m_BlockToUpdate[t] = false ;
}
else
nModifiedBlocks[t] = -1 ;
}
}
// Caso con tre mappe
else {
nModifiedBlocks.resize( m_nNumBlock + 1) ;
vLstTria.reserve( m_nNumBlock + 1) ;
TriaMatrix VecTriHold ;
VecTriHold.resize( m_nNumBlock) ;
bool bCalcInterBlock = false ;
// Calcolo i triangoli sui blocchi
for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) {
// Se il blocco deve essere processato
if ( bAllBlocks || m_BlockToUpdate[t]) {
// processo ...
vLstTria.emplace_back() ;
nModifiedBlocks[t] = int( vLstTria.size()) - 1 ;
m_InterBlockTria[t].clear() ;
#if 1
ExtMarchingCubes( int( t), vLstTria.back(), VecTriHold[t]) ;
// Flipping fra voxel interni
FlipEdgesII( VecTriHold[t]) ;
bCalcInterBlock = true ;
#else
MarchingCubes( int( t), vLstTria.back()) ;
#endif
m_BlockToUpdate[t] = false ;
}
else
nModifiedBlocks[t] = -1 ;
}
// Calcolo i triangoli di frontiera tra feature di blocchi diversi
// copio i triangoli di frontiera in una matrice gemella
// di m_InterBlockTria per avere sempre a disposizione
// i triangoli non flippati.
TriaMatrix InterBlockTria ;
if ( bCalcInterBlock) {
InterBlockTria = m_InterBlockTria ;
FlipEdgesBB( InterBlockTria) ;
}
// Inserisco in lista i triangoli di feature derivanti dai blocchi
for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) {
if ( nModifiedBlocks[t] >= 0) {
// ciclo sui voxel del blocco
for ( size_t t1 = 0 ; t1 < VecTriHold[t].size() ; ++ t1) {
// ciclo sulle componenti connesse del voxel
for ( size_t t2 = 0 ; t2 < VecTriHold[t][t1].vCompoTria.size() ; ++ t2) {
// ciclo sui triangoli delle componenti connesse
for ( size_t t3 = 0 ; t3 < VecTriHold[t][t1].vCompoTria[t2].size() ; ++ t3) {
// aggiungo triangolo alla lista
vLstTria[nModifiedBlocks[t]].emplace_back( VecTriHold[t][t1].vCompoTria[t2][t3]) ;
}
}
}
}
}
// Inserisco in lista i triangoli di frontiera tra feature di blocchi diversi
if ( bCalcInterBlock) {
vLstTria.resize( vLstTria.size() + 1) ;
size_t nPos = size_t( vLstTria.size() - 1) ;
for ( size_t t = 0 ; t < m_InterBlockTria.size() ; ++ t) {
for ( size_t t1 = 0 ; t1 < m_InterBlockTria[t].size() ; ++ t1) {
for ( size_t t2 = 0 ; t2 < m_InterBlockTria[t][t1].vCompoTria.size() ; ++ t2) {
for ( size_t t3 = 0 ; t3 < m_InterBlockTria[t][t1].vCompoTria[t2].size() ; ++ t3) {
if ( m_InterBlockTria[t][t1].vCompoTria[t2][t3].GetArea() > EPS_SMALL * EPS_SMALL) {
vLstTria[nPos].emplace_back( InterBlockTria[t][t1].vCompoTria[t2][t3]) ;
}
}
}
}
}
nModifiedBlocks.back() = int( nPos) ;
}
else
nModifiedBlocks.back() = - 1 ;
}
return true ;
}
//----------------------------------------------------------------------------
int
VolZmap::GetBlockCount( void) const
{
return m_nNumBlock + ( m_nMapNum == 1 ? 0 : 1) ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, int nDimChk, TRIA3DLIST& lstTria) const
{
// determino se è un semplice parallelepipedo
bool bIsSimple = true ;
double dBotZ ;
double dTopZ ;
for ( int i = 0 ; i < nDim1 && bIsSimple ; ++ i) {
for ( int j = 0 ; j < nDim2 && bIsSimple ; ++ j) {
int nPos = ( nPos1 + i) + ( nPos2 + j) * m_nNx[0] ;
if ( nPos > int( m_nDim[0]) ||
int( m_Values[0][nPos].size()) != 1)
bIsSimple = false ;
else if ( i == 0 && j == 0) {
dBotZ = m_Values[0][nPos][0].dMin ;
dTopZ = m_Values[0][nPos][0].dMax ;
}
else if ( abs( m_Values[0][nPos][0].dMin - dBotZ) > EPS_SMALL ||
abs( m_Values[0][nPos][0].dMax - dTopZ) > EPS_SMALL)
bIsSimple = false ;
}
}
// se semplice parallelepipedo
if ( bIsSimple) {
CalcChunkPrisms( nPos1, nPos2, nDim1, nDim2, lstTria) ;
}
// se chunk di dimensioni accettabili
else if ( nDimChk >= 4) {
int nNewDimChk = nDimChk / 2 ;
for ( int i = nPos1 ; i < int( nPos1 + nDim1) ; i += nNewDimChk) {
int nDimChunkX = min( nNewDimChk, int( nPos1 + nDim1) - i) ;
for ( int j = nPos2 ; j < int( nPos2 + nDim2) ; j += nNewDimChk) {
int nDimChunkY = min( nNewDimChk, int( nPos2 + nDim2) - j) ;
GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, nNewDimChk, lstTria) ;
}
}
}
// altrimenti
else {
// elaboro ogni singolo dexel
for ( int i = 0 ; i < nDim1 ; ++ i) {
for ( int j = 0 ; j < nDim2 ; ++ j) {
CalcDexelPrisms( nPos1 + i, nPos2 + j, lstTria) ;
}
}
}
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::CalcChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, TRIA3DLIST& lstTria) const
{
// verifiche sugli indici
if ( nPos1 < 0 || nPos1 + nDim1 > int( m_nNx[0]) || nPos2 < 0 || nPos2 + nDim2 > int( m_nNy[0]))
return false ;
int nPos = nPos1 + nPos2 * m_nNx[0] ;
if ( nPos < 0 || nPos >= int( m_nDim[0]))
return false ;
// calcolo coordinate punti
double dX = m_dStep * nPos1 ;
double dY = m_dStep * nPos2 ;
Point3d ptP1 = m_MapFrame[0].Orig() + dX * m_MapFrame[0].VersX() + dY * m_MapFrame[0].VersY() ;
Point3d ptP2 = ptP1 + nDim1 * m_dStep * m_MapFrame[0].VersX() ;
Point3d ptP3 = ptP2 + nDim2 * m_dStep * m_MapFrame[0].VersY() ;
Point3d ptP4 = ptP1 + nDim2 * m_dStep * m_MapFrame[0].VersY() ;
// creo le facce sopra e sotto
Vector3d vtDZt = m_Values[0][nPos][0].dMax * m_MapFrame[0].VersZ() ;
Vector3d vtDZb = m_Values[0][nPos][0].dMin * m_MapFrame[0].VersZ() ;
// faccia superiore P1t->P2t->P3t->P4t : sempre visibile
lstTria.emplace_back() ;
lstTria.back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_MapFrame[0].VersZ()) ;
lstTria.emplace_back() ;
lstTria.back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_MapFrame[0].VersZ()) ;
// faccia inferiore P1b->P4b->P3b->P2b : sempre visibile
lstTria.emplace_back() ;
lstTria.back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_MapFrame[0].VersZ()) ;
lstTria.emplace_back() ;
lstTria.back().Set( ptP3 + vtDZb, ptP2 + vtDZb, ptP1 + vtDZb, - m_MapFrame[0].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[0].VersY() ;
Point3d ptP3D = ptP2D + m_dStep * m_MapFrame[0].VersY() ;
AddDexelSideFace( nPosD, nPosEst, ptP2D, ptP3D, m_MapFrame[0].VersZ(), m_MapFrame[0].VersX(), lstTria) ;
}
for ( int i = 0 ; i < nDim1 ; ++ i) {
int nPosD = nPos + ( nDim2 - 1) * m_nNx[0] + i ;
int nPosNord = ( nPos2 + nDim2 - 1 < int( m_nNy[0] - 1) ? nPosD + m_nNx[0] : - 1) ;
Point3d ptP4D = ptP4 + i * m_dStep * m_MapFrame[0].VersX() ;
Point3d ptP3D = ptP4D + m_dStep * m_MapFrame[0].VersX() ;
AddDexelSideFace( nPosD, nPosNord, ptP3D, ptP4D, m_MapFrame[0].VersZ(), m_MapFrame[0].VersY(), lstTria) ;
}
for ( int j = 0 ; j < nDim2 ; ++ j) {
int nPosD = nPos + j * m_nNx[0] ;
int nPosWest = ( nPos1 > 0 ? nPosD - 1 : - 1) ;
Point3d ptP1D = ptP1 + j * m_dStep * m_MapFrame[0].VersY() ;
Point3d ptP4D = ptP1D + m_dStep * m_MapFrame[0].VersY() ;
AddDexelSideFace( nPosD, nPosWest, ptP4D, ptP1D, m_MapFrame[0].VersZ(), - m_MapFrame[0].VersX(), lstTria) ;
}
for ( int i = 0 ; i < nDim1 ; ++ i) {
int nPosD = nPos + i ;
int nPosSud = ( nPos2 > 0 ? nPosD - m_nNx[0] : - 1) ;
Point3d ptP1D = ptP1 + i * m_dStep * m_MapFrame[0].VersX() ;
Point3d ptP2D = ptP1D + m_dStep * m_MapFrame[0].VersX() ;
AddDexelSideFace( nPosD, nPosSud, ptP1D, ptP2D, m_MapFrame[0].VersZ(), - m_MapFrame[0].VersY(), lstTria) ;
}
//
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::CalcDexelPrisms( int nPos1, int nPos2, TRIA3DLIST& lstTria) const
{
// verifiche sugli indici
if ( nPos1 < 0 || nPos1 >= int( m_nNx[0]) || nPos2 < 0 || nPos2 >= int( m_nNy[0]))
return false ;
int nPos = nPos1 + nPos2 * m_nNx[0] ;
if ( nPos < 0 || nPos >= int( m_nDim[0]))
return false ;
// calcolo coordinate punto
double dX = m_dStep * nPos1 ;
double dY = m_dStep * nPos2 ;
Point3d ptP1 = m_MapFrame[0].Orig() + dX * m_MapFrame[0].VersX() + dY * m_MapFrame[0].VersY() ;
Point3d ptP2 = ptP1 + m_dStep * m_MapFrame[0].VersX() ;
Point3d ptP3 = ptP2 + m_dStep * m_MapFrame[0].VersY() ;
Point3d ptP4 = ptP1 + m_dStep * m_MapFrame[0].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[0].VersZ() ;
Vector3d vtDZb = m_Values[0][nPos][i].dMin * m_MapFrame[0].VersZ() ;
// faccia superiore P1t->P2t->P3t->P4t : sempre visibile
lstTria.emplace_back() ;
lstTria.back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_MapFrame[0].VersZ()) ;
lstTria.emplace_back() ;
lstTria.back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_MapFrame[0].VersZ()) ;
// faccia inferiore P1b->P4b->P3b->P2b : sempre visibile
lstTria.emplace_back() ;
lstTria.back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_MapFrame[0].VersZ()) ;
lstTria.emplace_back() ;
lstTria.back().Set( ptP3 + vtDZb, ptP2 + vtDZb, ptP1 + vtDZb, - m_MapFrame[0].VersZ()) ;
}
// creo le facce laterali
int nPosEst = ( nPos1 < int( m_nNx[0] - 1) ? nPos + 1 : - 1) ;
AddDexelSideFace( nPos, nPosEst, ptP2, ptP3, m_MapFrame[0].VersZ(), m_MapFrame[0].VersX(), lstTria) ;
int nPosNord = ( nPos2 < int( m_nNy[0] - 1) ? nPos + m_nNx[0] : - 1) ;
AddDexelSideFace( nPos, nPosNord, ptP3, ptP4, m_MapFrame[0].VersZ(), m_MapFrame[0].VersY(), lstTria) ;
int nPosWest = ( nPos1 > 0 ? nPos - 1 : - 1) ;
AddDexelSideFace( nPos, nPosWest, ptP4, ptP1, m_MapFrame[0].VersZ(), - m_MapFrame[0].VersX(), lstTria) ;
int nPosSud = ( nPos2 > 0 ? nPos - m_nNx[0] : - 1) ;
AddDexelSideFace( nPos, nPosSud, ptP1, ptP2, m_MapFrame[0].VersZ(), - m_MapFrame[0].VersY(), lstTria) ;
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Point3d& ptQ,
const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DLIST& lstTria) const
{
Intervals intFace ;
for ( int i = 0 ; i < int( m_Values[0][nPos].size()) ; i += 1)
intFace.Add( m_Values[0][nPos][i].dMin, m_Values[0][nPos][i].dMax) ;
if ( nPosAdj > 0) {
for ( int i = 0 ; i < int( m_Values[0][nPosAdj].size()) ; i += 1)
intFace.Subtract( m_Values[0][nPosAdj][i].dMin, m_Values[0][nPosAdj][i].dMax) ;
}
double dMin, dMax ;
bool bFound = intFace.GetFirst( dMin, dMax) ;
while ( bFound) {
Vector3d vtDZt = dMax * vtZ ;
Vector3d vtDZb = dMin * vtZ ;
lstTria.emplace_back() ;
lstTria.back().Set( ptP + vtDZb, ptQ + vtDZb, ptQ + vtDZt, vtNorm) ;
lstTria.emplace_back() ;
lstTria.back().Set( ptQ + vtDZt, ptP + vtDZt, ptP + vtDZb, vtNorm) ;
bFound = intFace.GetNext( dMin, dMax) ;
}
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const
{
if ( nBlock < 0 || nBlock >= int( m_nNumBlock))
return false ;
// Calcolo posizione del blocco nel reticolo
int nIBlock = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) % int( m_nFracLin[0]) ;
int nJBlock = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) / int( m_nFracLin[0]) ;
int nKBlock = ( nBlock / int( m_nFracLin[0] * m_nFracLin[1])) ;
// Calcolo limiti per l'indice i
int nStartI = nIBlock * int( m_nDexNumPBlock) - 1 ;
int nEndI = ( nIBlock + 1 == int( m_nFracLin[0]) ?
int( m_nNx[0]) : ( nIBlock + 1) * int( m_nDexNumPBlock)) ;
// Calcolo limiti per l'indice j
int nStartJ = nJBlock * int( m_nDexNumPBlock) - 1 ;
int nEndJ = ( nJBlock + 1 == int( m_nFracLin[1]) ?
int( m_nNy[0]) : ( nJBlock + 1) * int( m_nDexNumPBlock)) ;
// Calcolo limiti per l'indice k
int nStartK = nKBlock * int( m_nDexNumPBlock) - 1 ;
int nEndK = ( nKBlock + 1 == int( m_nFracLin[2]) ?
int( m_nNy[1]) : ( nKBlock + 1) * int( m_nDexNumPBlock)) ;
// Ciclo su tutti i voxel dello Zmap
for ( int i = nStartI ; i < nEndI ; ++ i) {
for ( int j = nStartJ ; j < nEndJ ; ++ j) {
for ( int k = nStartK ; k < nEndK ; ++ k) {
// 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}
} ;
// Classificazione dei vertici: interni o esterni al materiale
int nIndex = 0 ;
if ( IsThereMat( i, j, k))
nIndex |= ( 1 << 0) ;
if ( IsThereMat( i + 1, j, k))
nIndex |= ( 1 << 1) ;
if ( IsThereMat( i + 1, j + 1, k))
nIndex |= ( 1 << 2) ;
if ( IsThereMat( i, j + 1, k))
nIndex |= ( 1 << 3) ;
if ( IsThereMat( i, j, k + 1))
nIndex |= ( 1 << 4) ;
if ( IsThereMat( i + 1, j, k + 1))
nIndex |= ( 1 << 5) ;
if ( IsThereMat( i + 1, j + 1, k + 1))
nIndex |= ( 1 << 6) ;
if ( IsThereMat( i, j + 1, k + 1))
nIndex |= ( 1 << 7) ;
// Se vi è qualche intersezione fra segmenti e superficie
// continuo altrimenti passo al prossimo voxel
if ( EdgeTable[nIndex] == 0)
continue ;
static int intersections[12][2] = {
{ 0, 1 }, { 1, 2 }, { 2, 3 }, { 3, 0 }, { 4, 5 }, { 5, 6 },
{ 6, 7 }, { 7, 4 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }
} ;
// Ciclo sui segmenti
Point3d ptIntPoint[12] ;
for ( int EdgeIndex = 0 ; EdgeIndex < 12 ; ++ EdgeIndex) {
// Se il segmento non attraversa la superficie passo al successivo
if ( ! ( EdgeTable[nIndex] & ( 1 << EdgeIndex)))
continue ;
int n1 = intersections[EdgeIndex][0] ;
int n2 = intersections[EdgeIndex][1] ;
// Determino con precisione il punto di intersezione sullo spigolo
IntersPos( IndexCorner[n1], IndexCorner[n2], ptIntPoint[EdgeIndex]) ;
ptIntPoint[EdgeIndex].ToGlob( m_MapFrame[0]) ;
}
// Costruzione dei triangoli
for ( int TriIndex = 0 ; TriangleTableEn[nIndex][0][TriIndex] != - 1 ; TriIndex += 3) {
// Costruzione triangolo
int i0 = TriangleTableEn[nIndex][0][TriIndex + 2] ;
int i1 = TriangleTableEn[nIndex][0][TriIndex + 1] ;
int i2 = TriangleTableEn[nIndex][0][TriIndex] ;
Triangle3d CurrentTriangle ;
Vector3d vtN = ( ptIntPoint[i1] - ptIntPoint[i0]) ^ ( ptIntPoint[i2] - ptIntPoint[i1]) ;
vtN.Normalize() ;
vtN.ToGlob( m_MapFrame[0]) ;
// Il triangolo è pronto
CurrentTriangle.Set( ptIntPoint[i0], ptIntPoint[i1], ptIntPoint[i2], vtN) ;
// Aggiungo triangolo
lstTria.emplace_back( CurrentTriangle) ;
}
}
}
}
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) const
{
if ( nBlock < 0 || nBlock >= int( m_nNumBlock))
return false ;
// Calcolo i limiti sugli indici dei voxel del blocco
int nIJK[3] ;
// Vettore indici i,j,k del blocco
GetBlockIJKFromN( nBlock, nIJK) ;
// Vettore limiti sugli indici dei voxel del blocco
int nLimits[6] ;
GetBlockLimitsIJK( nIJK, nLimits) ;
// Ciclo su tutti i voxel dello Zmap
for ( int i = nLimits[0] ; i < nLimits[1] ; ++ i) {
for ( int j = nLimits[2] ; j < nLimits[3] ; ++ j) {
for ( int k = nLimits[4] ; k < nLimits[5] ; ++ k) {
// Riconoscimento dei voxel di frontiera
int nVoxIndexes[3] = { i, j, k} ;
bool bBoundary = IsAVoxelOnBoundary( nLimits, nVoxIndexes, true) ;
// 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}
} ;
// Classificazione dei vertici: interni o esterni al materiale
int nIndex = 0 ;
if ( IsThereMat( i, j, k))
nIndex |= ( 1 << 0) ;
if ( IsThereMat( i + 1, j, k))
nIndex |= ( 1 << 1) ;
if ( IsThereMat( i + 1, j + 1, k))
nIndex |= ( 1 << 2) ;
if ( IsThereMat( i, j + 1, k))
nIndex |= ( 1 << 3) ;
if ( IsThereMat( i, j, k + 1))
nIndex |= ( 1 << 4) ;
if ( IsThereMat( i + 1, j, k + 1))
nIndex |= ( 1 << 5) ;
if ( IsThereMat( i + 1, j + 1, k + 1))
nIndex |= ( 1 << 6) ;
if ( IsThereMat( i, j + 1, k + 1))
nIndex |= ( 1 << 7) ;
// Se vi è qualche intersezione fra segmenti e superficie
// continuo altrimenti passo al prossimo voxel.
if ( EdgeTable[nIndex] == 0)
continue ;
static int intersections[12][2] = {
{ 0, 1 }, { 1, 2 }, { 3, 2 }, { 0, 3 }, { 4, 5 }, { 5, 6 },
{ 7, 6 }, { 4, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }
} ;
// Array di strutture punto di intersezione e normale alla superficie in esso.
VectorField VecField[12] ;
// Flag sulla regolatrità dei campi scalare e vettoriale:
// se i campi sono regolari esso resta vero, altrimenti
// assume il valore falso.
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)))
continue ;
int n1 = intersections[EdgeIndex][0] ;
int n2 = intersections[EdgeIndex][1] ;
bool bN1 = ( ( nIndex & ( 1 << n1)) != 0) ;
// Determino con precisione il punto di intersezione sullo spigolo,
// se i campi scalare e vettoriale non sono regolari bReg diviene falso.
if ( ! IntersPos( IndexCorner[n1], IndexCorner[n2], bN1,
VecField[EdgeIndex].ptInt,
VecField[EdgeIndex].vtNorm))
bReg = false ;
// Riporto punti e normali nel sistema locale in cui
// è immerso lo Zmap col suo sistema di riferimento.
VecField[EdgeIndex].ptInt.ToGlob( m_MapFrame[0]) ;
VecField[EdgeIndex].vtNorm.ToGlob( m_MapFrame[0]) ;
}
// 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
// nel caso di assenza di sharp feature, della (i+1)-esima
// componente connessa.
VectorField CompoVert[6][12] ;
VectorField CompoTriVert[6][17] ;
// Arrey numero di vertici della base del fan per componente
// connessa: nVertComp[i] contiene il numero di vertici
// della base del fan della (i+1)-esima componente connessa.
int nVertComp[6] ;
// Matrice di indici dei punti: serve per
// la gestione del caso
int nIndArrey[6][4] ;
int nExtTabOff = nComponents ;
int nStdTabOff = 0 ;
// Carico le matrici CompoVert e CompoTriVert
for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) {
// Numero vertici per componenti
nVertComp[nCompCount - 1] = TriangleTableEn[nIndex][1][nCompCount] ;
// Riempio il nCompCount-esimo vettore di vertici della base del fan
for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount)
CompoVert[nCompCount - 1][nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1]] ;
// Serve per la gestione del caso ...
if ( nVertComp[nCompCount - 1] == 4) {
for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount)
nIndArrey[nCompCount - 1][nVertCount] = TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1] ;
}
// Riempio il nCompCount-esimo vettore di vertici dei triangoli in assenza di
// sharp feature: in una mesh di triangoli con n vertici vi sono n - 2 triangoli.
for ( int nVertCount = 0 ; nVertCount < 3 * ( nVertComp[nCompCount - 1] - 2) ; nVertCount += 3) {
CompoTriVert[nCompCount - 1][nVertCount] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+2]] ;
CompoTriVert[nCompCount - 1][nVertCount+1] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+1]] ;
CompoTriVert[nCompCount - 1][nVertCount+2] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount]] ;
}
// Aggiorno gli offsets per raggiungere i
// vertici della componente successiva.
nExtTabOff += nVertComp[nCompCount - 1] ;
nStdTabOff += 3 * ( nVertComp[nCompCount - 1] - 2) ;
}
// Test sulla topologia
if ( nAllConfig[nIndex] == 3) {
Vector3d vtCmpAvg0, vtCmpAvg1 ;
// Verifico se i versori delle componenti sono tutti
// più o meno concordi (per esserlo non devono esserci
// due vettori di una medesima componente con prodotto
// scalare inferiore a 0.7).
bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ;
bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
// Se i versori di entrambe le componenti sono concordi
// ha senso parlare di vettori medi, altrimenti non ha
// senso. Se non ha senso parlare di vettori medi non
// ha senso parlare di prodotti scalari fra loro,
// quindi pongo il loro prodotto a un valore assurdo -2
// (il prodotto scalare fra versori ha modulo non superiore
// a uno).
double dScProd = - 2 ;
if ( bTest0 && bTest1)
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
double dThreshold = 0.7 ;
if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
int nt = 0 ;
while ( nIndexVsIndex3[nt][0] != nIndex)
++ nt ;
int nRotCase = nIndexVsIndex3[nt][1] ;
nComponents = Cases3Plus[nRotCase][1][0] ;
// Riaggiorno gli offsets
nExtTabOff = nComponents ;
nStdTabOff = 0 ;
// Modifico le matrici
for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
// Numero vertici per componenti
nVertComp[nC - 1] = Cases3Plus[nRotCase][1][nC] ;
// Matrice dei vertici della base del fan
for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert)
CompoVert[nC - 1][nFanVert] = VecField[Cases3Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ;
// Matrici dei vertici dei triangoli in assenza di sharp feature
for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) {
CompoTriVert[nC - 1][nTriVert] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
}
// Aggiorno gli offsets per raggiungere i
// vertici della componente successiva.
nExtTabOff += nVertComp[nC - 1] ;
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
}
}
}
else if ( nAllConfig[nIndex] == 6) {
// Procedura analoga a quella della configurazione 3
Vector3d vtCmpAvg0, vtCmpAvg1 ;
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.7) ;
bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
double dScProd = - 2 ;
if ( bTest0 && bTest1)
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
double dThreshold = 0.7 ;
if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
int nt = 0 ;
while ( nIndexVsIndex6[nt][0] != nIndex)
++ nt ;
int nRotCase = nIndexVsIndex6[nt][1] ;
// Costruzione dei triangoli
for ( int TriIndex = 0 ; TriIndex < 15 ; TriIndex += 3) {
// Costruzione triangolo
int i0 = Cases6Plus[nRotCase][TriIndex + 2] ;
int i1 = Cases6Plus[nRotCase][TriIndex + 1] ;
int i2 = Cases6Plus[nRotCase][TriIndex] ;
Triangle3d CurrentTriangle ;
// Il triangolo è pronto
CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ;
bool bV = CurrentTriangle.Validate( true) ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
}
continue ;
}
}
else if ( nAllConfig[nIndex] == 10) {
Vector3d vtCmpAvg0, vtCmpAvg1 ;
// Verifico se i versori delle componenti sono tutti
// più o meno concordi (per esserlo non devono esserci
// due vettori di una medesima componente con prodotto
// scalare inferiore a 0). decidere se 0.0 o 0.7
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ;
bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ;
if ( ! ( bTest0 && bTest1)) {
int nt = 0 ;
while ( nIndexVsIndex10[nt][0] != nIndex)
++ nt ;
int nRotCase = nIndexVsIndex10[nt][1] ;
// Riaggiorno gli offsets
nExtTabOff = 2 ;
nStdTabOff = 0 ;
// Modifico le matrici
for ( int nC = 1 ; nC <= 2 ; ++ nC) {
// Numero vertici per componenti
nVertComp[nC - 1] = Cases10Plus[nRotCase][1][nC] ;
// Matrice dei vertici della base del fan
for ( int nFanVert = 0 ; nFanVert < 4 ; ++ nFanVert)
CompoVert[nC - 1][nFanVert] = VecField[Cases10Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ;
// Matrici dei vertici dei triangoli in assenza di sharp feature
for ( int nTriVert = 0 ; nTriVert < 6 ; nTriVert += 3) {
CompoTriVert[nC - 1][nTriVert] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
}
// Aggiorno gli offsets per raggiungere i
// vertici della componente successiva.
nExtTabOff += nVertComp[nC - 1] ;
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
}
}
}
// Numero di feature nel voxel: al più vi è una feature per componente connessa.
int nInnerFeatureInVoxel = 0 ;
int nBorderFeatureInVoxel = 0 ;
// Ciclo sulle componenti
for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) {
int nFeatureType = NO_FEATURE ;
// Se i componenti sono regolari valuto le normali per stabilire se eseguire ExtMC o MC
if ( bReg)
nFeatureType = TestOnNormal( CompoVert[nCompCount - 1], nVertComp[nCompCount - 1]) ;
// Controllo per il caso piano su una griglia
// con versori normali a due a due paralleli.
bool bGridControl = true ;
if ( nFeatureType != NO_FEATURE) {
if ( nVertComp[nCompCount - 1] == 4) {
// Ordino i 4 indici in senso crescente
for ( int nSrtInd1 = 0 ; nSrtInd1 < nVertComp[nCompCount - 1] - 1 ; ++ nSrtInd1) {
for ( int nSrtInd2 = nSrtInd1 + 1 ; nSrtInd2 < nVertComp[nCompCount - 1] ; ++ nSrtInd2) {
if ( nIndArrey[nCompCount - 1][nSrtInd1] > nIndArrey[nCompCount - 1][nSrtInd2])
swap( nIndArrey[nCompCount - 1][nSrtInd1], nIndArrey[nCompCount - 1][nSrtInd2]) ;
}
}
if ( ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 2 &&
nIndArrey[nCompCount - 1][2] == 9 && nIndArrey[nCompCount - 1][3] == 10) ||
( nIndArrey[nCompCount - 1][0] == 4 && nIndArrey[nCompCount - 1][1] == 6 &&
nIndArrey[nCompCount - 1][2] == 9 && nIndArrey[nCompCount - 1][3] == 10) ||
( nIndArrey[nCompCount - 1][0] == 4 && nIndArrey[nCompCount - 1][1] == 6 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 2 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 9 ) ||
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 10 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 5 && nIndArrey[nCompCount - 1][1] == 7 &&
nIndArrey[nCompCount - 1][2] == 10 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 5 && nIndArrey[nCompCount - 1][1] == 7 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 9 )) {
VectorField LocVecF[12], LocCompV[12] ;
for ( int LocInd = 0 ; LocInd < 12 ; ++ LocInd) {
LocVecF[LocInd] = VecField[LocInd] ;
LocCompV[LocInd] = CompoVert[nCompCount - 1][LocInd] ;
LocVecF[LocInd].ptInt.ToLoc( m_MapFrame[0]) ;
LocVecF[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ;
LocCompV[LocInd].ptInt.ToLoc( m_MapFrame[0]) ;
LocCompV[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ;
}
if ( ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm *LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) < EPS_SMALL &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) < EPS_SMALL) ||
( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm) < EPS_SMALL &&
abs( LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) < EPS_SMALL)) {
Point3d ptBarycenter = ( LocCompV[0].ptInt + LocCompV[1].ptInt +
LocCompV[2].ptInt + LocCompV[3].ptInt) / 4 ;
if ( abs( LocCompV[0].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) {
double dXBar = ptBarycenter.x / m_dStep - 0.5 ;
int nBarLimSup = int( m_nNx[0]) ;
int nBarInd = 0 ;
while ( nBarInd < nBarLimSup) {
double dXInd = double( nBarInd) ;
if ( abs( dXInd - dXBar) < EPS_SMALL) {
bGridControl = false ;
break ;
}
++ nBarInd ;
}
}
else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) {
double dYBar = ptBarycenter.y / m_dStep - 0.5 ;
int nBarLimSup = int( m_nNy[0]) ;
int nBarInd = 0 ;
while ( nBarInd < nBarLimSup) {
double dYInd = double( nBarInd) ;
if ( abs( dYInd - dYBar) < EPS_SMALL) {
bGridControl = false ;
break ;
}
++ nBarInd ;
}
}
else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[3].ptInt.z - ptBarycenter.z) < EPS_SMALL) {
double dZBar = ptBarycenter.z / m_dStep - 0.5 ;
int nBarLimSup = int( m_nNy[1]) ;
int nBarInd = 0 ;
while ( nBarInd < nBarLimSup) {
double dZInd = double( nBarInd) ;
if ( abs( dZInd - dZBar) < EPS_SMALL) {
bGridControl = false ;
break ;
}
++ nBarInd ;
}
}
}
}
else if ( ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 1 &&
nIndArrey[nCompCount - 1][2] == 4 && nIndArrey[nCompCount - 1][3] == 5) ||
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 2 &&
nIndArrey[nCompCount - 1][2] == 5 && nIndArrey[nCompCount - 1][3] == 6) ||
( nIndArrey[nCompCount - 1][0] == 2 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 6 && nIndArrey[nCompCount - 1][3] == 7) ||
( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 4 && nIndArrey[nCompCount - 1][3] == 7)) {
VectorField LocVecF[12], LocCompV[12] ;
for ( int LocInd = 0 ; LocInd < 12 ; ++ LocInd) {
LocVecF[LocInd] = VecField[LocInd] ;
LocCompV[LocInd] = CompoVert[nCompCount - 1][LocInd] ;
LocVecF[LocInd].ptInt.ToLoc( m_MapFrame[0]) ;
LocVecF[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ;
LocCompV[LocInd].ptInt.ToLoc( m_MapFrame[0]) ;
LocCompV[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ;
}
if ( ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) < EPS_SMALL &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) < EPS_SMALL) ||
( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm) < EPS_SMALL &&
abs( LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) < EPS_SMALL)) {
Point3d ptBarycenter = ( LocCompV[0].ptInt + LocCompV[1].ptInt +
LocCompV[2].ptInt + LocCompV[3].ptInt) / 4 ;
if ( abs( LocCompV[0].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) {
double dXBar = ptBarycenter.x / m_dStep - 0.5 ;
int nBarLimSup = int( m_nNx[0]) ;
int nBarInd = 0 ;
while ( nBarInd < nBarLimSup) {
double dXInd = double( nBarInd) ;
if ( abs( dXInd - dXBar) < EPS_SMALL) {
bGridControl = false ;
break ;
}
++ nBarInd ;
}
}
else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) {
double dYBar = ptBarycenter.y / m_dStep - 0.5 ;
int nBarLimSup = int( m_nNy[0]) ;
int nBarInd = 0 ;
while ( nBarInd < nBarLimSup) {
double dYInd = double( nBarInd) ;
if ( abs( dYInd - dYBar) < EPS_SMALL) {
bGridControl = false ;
break ;
}
++ nBarInd ;
}
}
else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[3].ptInt.z - ptBarycenter.z) < EPS_SMALL) {
double dZBar = ptBarycenter.z / m_dStep - 0.5 ;
int nBarLimSup = int( m_nNy[1]) ;
int nBarInd = 0 ;
while ( nBarInd < nBarLimSup) {
double dZInd = double( nBarInd) ;
if ( abs( dZInd - dZBar) < EPS_SMALL) {
bGridControl = false ;
break ;
}
++ nBarInd ;
}
}
}
}
}
}
// Flag ExtMC
bool bExtMC = ( nFeatureType != NO_FEATURE) && bGridControl ;
// Extended MC
if ( bExtMC) {
// Passo al sistema di riferimento del baricentro
Point3d ptGravityCenter( 0, 0, 0) ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni)
ptGravityCenter = ptGravityCenter + CompoVert[nCompCount - 1][ni].ptInt ;
ptGravityCenter = ptGravityCenter / nVertComp[nCompCount - 1] ;
Vector3d vtO = ptGravityCenter - ORIG ;
Point3d ptTrasf[12] ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni)
ptTrasf[ni] = CompoVert[nCompCount - 1][ni].ptInt - vtO ;
// Preparo le matrici per il sistema
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> dSystemMatrix ;
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> dSystemVector ;
typedef Eigen::JacobiSVD<dSystemMatrix> DecomposerSVD ;
dSystemMatrix dMatrixN, dMatrixU, dMatrixV ;
dSystemVector dKnownVector, dUnknownVector, dSingularValue ;
dMatrixN.resize( nVertComp[nCompCount - 1], 3) ;
dKnownVector.resize( nVertComp[nCompCount - 1], 1) ;
dUnknownVector.resize( 3, 1) ;
// Studio del caso 4 punti su un piano
int nEqual = 0 ;
int nPosD ;
Vector3d vtD, vtE ;
if ( nVertComp[nCompCount - 1] == 4 && nFeatureType == 2) {
int nPosEq ;
for ( int ni = 0 ; ni < 2 ; ++ ni) {
for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) {
if ( AreSameVectorApprox( CompoVert[nCompCount - 1][ni].vtNorm,
CompoVert[nCompCount - 1][nj].vtNorm)) {
nEqual ++ ;
nPosEq = ni ;
}
}
if ( nEqual == 2)
break ;
else
nEqual = 0 ;
}
if ( nEqual == 2) {
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni)
if ( ! AreSameVectorApprox( CompoVert[nCompCount - 1][ni].vtNorm,
CompoVert[nCompCount - 1][nPosEq].vtNorm)) {
nPosD = ni ;
vtD = CompoVert[nCompCount - 1][ni].vtNorm ;
vtE = CompoVert[nCompCount - 1][nPosEq].vtNorm ;
}
}
}
double dDot = abs( ( CompoVert[nCompCount - 1][1].ptInt - CompoVert[nCompCount - 1][0].ptInt) *
( ( CompoVert[nCompCount - 1][2].ptInt - CompoVert[nCompCount - 1][1].ptInt) ^
( CompoVert[nCompCount - 1][3].ptInt - CompoVert[nCompCount - 1][2].ptInt))) ;
// Caso superficie piana
if ( nVertComp[nCompCount - 1] == 4 && nEqual == 2 && dDot < EPS_SMALL) {
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
if ( ni != nPosD) {
dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ;
dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ;
dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ;
dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * ( ptTrasf[ni] - ORIG) ;
}
else {
dMatrixN( ni, 0) = vtE.x ;
dMatrixN( ni, 1) = vtE.y ;
dMatrixN( ni, 2) = vtE.z ;
dKnownVector( ni) = vtE * ( ptTrasf[ni] - ORIG) ;
}
}
}
// Caso generale
else {
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ;
dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ;
dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ;
dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * ( ptTrasf[ni] - ORIG) ;
}
}
DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ;
dMatrixU = svd.matrixU() ;
dMatrixV = svd.matrixV() ;
dSingularValue = svd.singularValues() ;
double dS0 = dSingularValue( 0) ;
double dS1 = dSingularValue( 1) ;
double dS2 = dSingularValue( 2) ;
// Se la feature è un edge annullo il valore singolare minore.
if ( nFeatureType == EDGE) {
dSingularValue( 2) = 0 ;
}
// Back substitution: risolvo il sistema USV*x = b
// Calcolo U^T b
double vTemp[3] ;
for ( int ni = 0 ; ni < 3 ; ++ ni) {
double s = 0 ;
if ( dSingularValue( ni) > 0) {
for ( int nj = 0 ; nj < nVertComp[nCompCount - 1] ; ++ nj)
s += dMatrixU( nj, ni) * dKnownVector( nj) ;
s /= dSingularValue( ni) ;
}
vTemp[ni] = s ;
}
// Moltiplico per V
for ( int ni = 0 ; ni < 3 ; ++ ni) {
double s = 0 ;
for ( int nj = 0 ; nj < 3 ; ++ nj)
s += dMatrixV( ni, nj) * vTemp[nj] ;
dUnknownVector( ni) = s ;
}
// Esprimo la soluzione nel sistema di riferimento
// in cui è immerso quello dello dello z-map.
Point3d ptSol( dUnknownVector( 0) + vtO.x,
dUnknownVector( 1) + vtO.y,
dUnknownVector( 2) + vtO.z) ;
// Flag sulla distanza del vertice dal
// baricentro del sistema di punti
bool bOutside = false ;
// Vettore Baricentro-Feature.
Vector3d vtFeature( dUnknownVector( 0),
dUnknownVector( 1),
dUnknownVector( 2)) ;
double dDistFeature = vtFeature.Len() ;
const double dMaxDist = sqrt( 3) * m_dStep ;
// Limito la feature a una distanza dal
// baricentro pari alla diagonale del voxel.
if ( dDistFeature > dMaxDist)
bOutside = true ;
Triangle3d CurrentTriangle ;
TRIA3DVECTOR triContainer ;
// Costruisco triangoli di prova
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
// Il triangolo è pronto
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
CurrentTriangle.Validate( true) ;
// Aggiungo triangolo al vettore temporaneo
triContainer.emplace_back( CurrentTriangle) ;
}
// Controllo delle inversioni dei triangoli
bool bInversione = false ;
bool bDangerInversion = false ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - EPS_SMALL ||
triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL)
bInversione = true ;
}
if ( bInversione &&
( nVertComp[nCompCount - 1] == 3 ||
nVertComp[nCompCount - 1] == 4)) {
int nNegDot = 0 ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] - 1 ; ++ ni) {
for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) {
if ( CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL)
nNegDot ++ ;
}
}
if ( nNegDot == nVertComp[nCompCount - 1] - 1) {
Point3d ptSolZMapFrame = ptSol ;
ptSolZMapFrame.ToLoc( m_MapFrame[0]) ;
if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) {
Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
vtNullSpace.ToLoc( m_MapFrame[0]) ;
double dParInt1, dParInt2 ;
Point3d ptVoxMin( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, ( k + 0.5) * m_dStep) ;
Point3d ptVoxMax( ( i + 1.5) * m_dStep, ( j + 1.5) * m_dStep, ( k + 1.5) * m_dStep) ;
if ( IntersLineBox( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) {
triContainer.resize( 0) ;
double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 :
dParInt2 + ( dParInt1 - dParInt2) / 100 ;
Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ;
ptNewSol.ToGlob( m_MapFrame[0]) ;
ptSol = ptNewSol ;
// Costruisco triangoli di prova
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
// Il triangolo è pronto
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
CurrentTriangle.Validate( true) ;
// Aggiungo triangolo al vettore temporaneo
triContainer.emplace_back( CurrentTriangle) ;
}
}
else {
int nAdjVoxI, nAdjVoxJ, nAdjVoxK ;
if ( GetPointVoxel( ptSolZMapFrame, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) {
// Classificazione del voxel adiacente
int nAdjIndex = 0 ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK))
nAdjIndex |= ( 1 << 0) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK))
nAdjIndex |= ( 1 << 1) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK))
nAdjIndex |= ( 1 << 2) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK))
nAdjIndex |= ( 1 << 3) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 4) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 5) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 6) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 7) ;
// Se il voxel è pieno
if ( EdgeTable[nAdjIndex] != 0)
bDangerInversion = true ;
}
}
}
}
else {
if ( nVertComp[nCompCount - 1] == 4) {
Point3d ptSolZMapFrame = ptSol ;
ptSolZMapFrame.ToLoc( m_MapFrame[0]) ;
if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) {
Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
vtNullSpace.ToLoc( m_MapFrame[0]) ;
double dParInt1, dParInt2 ;
Point3d ptVoxMin( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, ( k + 0.5) * m_dStep) ;
Point3d ptVoxMax( ( i + 1.5) * m_dStep, ( j + 1.5) * m_dStep, ( k + 1.5) * m_dStep) ;
if ( IntersLineBox( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) {
triContainer.resize( 0) ;
double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100:
dParInt2 + ( dParInt1 - dParInt2) / 100 ;
Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ;
ptNewSol.ToGlob( m_MapFrame[0]) ;
ptSol = ptNewSol ;
// Costruisco triangoli di prova
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
// Il triangolo è pronto
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
CurrentTriangle.Validate( true) ;
// Aggiungo triangolo al vettore temporaneo
triContainer.emplace_back( CurrentTriangle) ;
}
}
}
}
}
}
// Valuto normali: questo è ancora un controllo
// sulle normali, se risultano in tutti i punti
// approssimativamente uguali passiamo alla
// routine standard
int nContSize = int( triContainer.size()) ;
bool bPlane = true ;
for ( int ni = 0 ; ni < nContSize - 1 ; ++ ni) {
for ( int nj = ni + 1 ; nj < nContSize ; ++ nj) {
Vector3d vtI = triContainer[ni].GetN() ;
Vector3d vtJ = triContainer[nj].GetN() ;
if ( ! AreSameVectorApprox( vtI, vtJ)) {
bPlane = false ;
break ;
}
}
if ( ! bPlane)
break ;
}
// Se la feature non è fuori dai limiti
// e i triangoli formano una superficie
// non piana confermo ExtMC
if ( ( ! ( bOutside || bPlane)) && ( ! bDangerInversion)) {
TRIA3DVECTOR vInnerTriaTemp, vBorderTriaTemp ;
for ( int ni = 0 ; ni < nContSize ; ++ ni) {
// Se l'area è non nulla aggiungo il triangolo
// a uno dei due vettori.
if ( triContainer[ni].GetArea() > EPS_SMALL * EPS_SMALL) {
int nVoxIJK[3] = { i, j, k} ;
Point3d ptTemp = ptSol ;
ptTemp.ToLoc( m_MapFrame[0]) ;
Triangle3d trTemp = triContainer[ni] ;
trTemp.ToLoc( m_MapFrame[0]) ;
if ( ! IsATriangleOnBorder( trTemp, ptTemp, nLimits, nVoxIJK, bBoundary))
vInnerTriaTemp.emplace_back( triContainer[ni]) ;
else
vBorderTriaTemp.emplace_back( triContainer[ni]) ;
}
}
// Metto i triangoli negli appositi contenitori:
// Triangoli interni
if ( vInnerTriaTemp.size() > 0) {
// Aggiorno il numero di feature.
++ nInnerFeatureInVoxel ;
// Se siamo in presenza della prima feature del
// voxel, ridimensiono il vettore che contiene
// la struttura dati del voxel.
if ( nInnerFeatureInVoxel == 1) {
triHold.resize( triHold.size() + 1) ;
// Questi dati dipendono solo dal voxel,
// quindi sono validi per tutte le
// componenti che vi appartengono.
int nCurrent = int( triHold.size()) - 1 ;
triHold[nCurrent].i = i ;
triHold[nCurrent].j = j ;
triHold[nCurrent].k = k ;
}
// Indice che identifica la posizione del voxel
// nel vector.
int nCurrent = int( triHold.size()) - 1 ;
// Aggiungo vertice della componente
// connessa alla lista dei vertici.
triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ;
int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ;
int nNewFeatureNum = nOldFeatureNum + 1 ;
// Aggiungo una componente per il vettore
// dei triangoli della componente connessa.
triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ;
for ( int ni = 0 ; ni < int( vInnerTriaTemp.size()) ; ++ ni)
triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vInnerTriaTemp[ni]) ;
}
// Triangoli di frontiera
if ( vBorderTriaTemp.size() > 0) {
// Aggiorno il numero di feature.
++ nBorderFeatureInVoxel ;
// Se siamo in presenza della prima feature del
// voxel, ridimensiono il vettore che contiene
// la struttura dati del voxel.
if ( nBorderFeatureInVoxel == 1) {
m_InterBlockTria[nBlock].resize( m_InterBlockTria[nBlock].size() + 1) ;
// Questi dati dipendono solo dal voxel,
// quindi sono validi per tutte le
// componenti che vi appartengono.
int nCurrent = int( m_InterBlockTria[nBlock].size()) - 1 ;
m_InterBlockTria[nBlock][nCurrent].i = i ;
m_InterBlockTria[nBlock][nCurrent].j = j ;
m_InterBlockTria[nBlock][nCurrent].k = k ;
}
// Indice che identifica la posizione del voxel
// nel vector.
int nCurrent = int( m_InterBlockTria[nBlock].size()) - 1 ;
// Aggiungo vertice della componente
// connessa alla lista dei vertici.
m_InterBlockTria[nBlock][nCurrent].ptCompoVert.emplace_back( ptSol) ;
int nOldFeatureNum = int( m_InterBlockTria[nBlock][nCurrent].vCompoTria.size()) ;
int nNewFeatureNum = nOldFeatureNum + 1 ;
// Aggiungo una componente per il vettore
// dei triangoli della componente connessa.
m_InterBlockTria[nBlock][nCurrent].vCompoTria.resize( nNewFeatureNum) ;
for ( int ni = 0 ; ni < int( vBorderTriaTemp.size()) ; ++ ni)
m_InterBlockTria[nBlock][nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vBorderTriaTemp[ni]) ;
}
}
// ExtMC non confermato, si passa a MC
else {
// Costruzione dei triangoli
for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {
Triangle3d CurrentTriangle ;
// Il triangolo è pronto
CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt,
CompoTriVert[nCompCount - 1][TriIndex+1].ptInt,
CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ;
bool bV = CurrentTriangle.Validate( true) ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
}
}
}
// Standard MC
else {
// Costruzione dei triangoli
for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {
Triangle3d CurrentTriangle ;
// Il triangolo è pronto
CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt,
CompoTriVert[nCompCount - 1][TriIndex+1].ptInt,
CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ;
bool bV = CurrentTriangle.Validate( true) ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
}
}
}
}
}
}
return true ;
}
//
//----------------------------------------------------------------------------
bool
VolZmap::FlipEdgesII( TriHolder& TriHold) const
{
// Numero di voxel in cui si presentano sharp feature
int nVoxelNum = int( TriHold.size()) ;
// Ciclo su tali voxel
for ( int n1 = 0 ; n1 < nVoxelNum ; ++ n1) {
for ( int n2 = n1 ; n2 < nVoxelNum ; ++ n2) {
// Se i voxel sono adiacenti proseguo
if ( ( TriHold[n2].i >= TriHold[n1].i && TriHold[n2].i <= TriHold[n1].i + 1) ||
( TriHold[n2].j >= TriHold[n1].j && TriHold[n2].j <= TriHold[n1].j + 1) ||
( TriHold[n2].k >= TriHold[n1].k && TriHold[n2].k <= TriHold[n1].k + 1)) {
// Numero delle componenti connesse nei due voxel
int nNumCompo1 = int( TriHold[n1].ptCompoVert.size()) ;
int nNumCompo2 = int( TriHold[n2].ptCompoVert.size()) ;
int nCompo1 = 0 ;
// Ciclo sulle componenti
for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) {
int nCompo2 = ( n1 == n2 ? nCompo1 + 1 : 0) ;
for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) {
// Numero di triangoli per le componenti connesse
int nTriNum1 = int( TriHold[n1].vCompoTria[nCompo1].size()) ;
int nTriNum2 = int( TriHold[n2].vCompoTria[nCompo2].size()) ;
for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) {
bool bModified = false ;
for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) {
INTVECTOR SharedIndex ;
for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) {
for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) {
Point3d ptP1 = TriHold[n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ;
Point3d ptP2 = TriHold[n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ;
if ( AreSamePointEpsilon( ptP1, ptP2, EPS_ZERO)) {
Point3d ptVert1 = TriHold[n1].ptCompoVert[nCompo1] ;
Point3d ptVert2 = TriHold[n2].ptCompoVert[nCompo2] ;
if ( ! ( AreSamePointEpsilon( ptP1, ptVert1, EPS_ZERO) ||
AreSamePointEpsilon( ptP2, ptVert2, EPS_ZERO))) {
SharedIndex.emplace_back( nVert1) ;
SharedIndex.emplace_back( nVert2) ;
}
}
if ( SharedIndex.size() > 2)
break ;
}
if ( SharedIndex.size() > 2)
break ;
}
// Si deve operare la modifica dei triangoli
if ( SharedIndex.size() > 2) {
// verifico che i due lati adiacenti siano controversi
int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ;
if ( nProd != 1 && nProd != - 2 && nProd != 4 && nProd != 0) {
TriHold[n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0],
TriHold[n2].ptCompoVert[nCompo2]) ;
TriHold[n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3],
TriHold[n1].ptCompoVert[nCompo1]) ;
TriHold[n1].vCompoTria[nCompo1][nTri1].Validate( true) ;
TriHold[n2].vCompoTria[nCompo2][nTri2].Validate( true) ;
bModified = true ;
break ;
}
}
}
if ( bModified)
break ;
}
}
}
}
}
}
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const
{
// Numero di blocchi
size_t nBlocksNum = InterTria.size() ;
// ciclo sui blocchi
for ( size_t tFB = 0 ; tFB < nBlocksNum - 1 ; ++ tFB) {
int nFBijk[3] ;
GetBlockIJKFromN( int( tFB), nFBijk) ;
for ( size_t tLB = tFB + 1 ; tLB < nBlocksNum ; ++ tLB) {
int nLBijk[3] ;
GetBlockIJKFromN( int( tLB), nLBijk) ;
// Se i blocchi non sono adiacenti salto l'iterazione
if ( ! ( abs( nFBijk[0] - nLBijk[0]) <= 1 &&
abs( nFBijk[1] - nLBijk[1]) <= 1 &&
abs( nFBijk[2] - nLBijk[2]) <= 1))
continue ;
// Numero di voxel nei blocchi correnti
size_t nVoxelNumFB = InterTria[tFB].size() ;
size_t nVoxelNumLB = InterTria[tLB].size() ;
// Ciclo sui voxel dei due blocchi
for ( size_t tVFB = 0 ; tVFB < nVoxelNumFB ; ++ tVFB) {
for ( size_t tVLB = 0 ; tVLB < nVoxelNumLB ; ++ tVLB) {
// Se i voxel non sono adiacenti salto l'iterazione
if ( ! ( abs( InterTria[tFB][tVFB].i - InterTria[tLB][tVLB].i) <= 1 &&
abs( InterTria[tFB][tVFB].j - InterTria[tLB][tVLB].j) <= 1 &&
abs( InterTria[tFB][tVFB].k - InterTria[tLB][tVLB].k) <= 1))
continue ;
// Numero di componenti connesse dei voxel
size_t nCompoVFBNum = InterTria[tFB][tVFB].ptCompoVert.size() ;
size_t nCompoVLBNum = InterTria[tLB][tVLB].ptCompoVert.size() ;
// Ciclo sulle componenti connesse
for ( size_t tCmpF = 0 ; tCmpF < nCompoVFBNum ; ++ tCmpF) {
for ( size_t tCmpL = 0 ; tCmpL < nCompoVLBNum ; ++ tCmpL) {
// Numero di triangoli delle componenti connesse
size_t nTriFBNum = InterTria[tFB][tVFB].vCompoTria[tCmpF].size() ;
size_t nTriLBNum = InterTria[tLB][tVLB].vCompoTria[tCmpL].size() ;
for ( size_t tTriFB = 0 ; tTriFB < nTriFBNum ; ++ tTriFB) {
bool bModified = false ;
for ( size_t tTriLB = 0 ; tTriLB < nTriLBNum ; ++ tTriLB) {
INTVECTOR SharedIndex ;
for ( int nVertF = 0 ; nVertF < 3 ; ++ nVertF) {
for ( int nVertL = 0 ; nVertL < 3 ; ++ nVertL) {
Point3d ptPF = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( nVertF) ;
Point3d ptPL = InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( nVertL) ;
if ( AreSamePointEpsilon( ptPF, ptPL, EPS_ZERO)) {
Point3d ptVertF = InterTria[tFB][tVFB].ptCompoVert[tCmpF] ;
Point3d ptVertL = InterTria[tLB][tVLB].ptCompoVert[tCmpL] ;
if ( ! ( AreSamePointEpsilon( ptPF, ptVertF, EPS_ZERO) ||
AreSamePointEpsilon( ptPL, ptVertL, EPS_ZERO))) {
SharedIndex.emplace_back( nVertF) ;
SharedIndex.emplace_back( nVertL) ;
}
}
if ( SharedIndex.size() > 2)
break ;
}
if ( SharedIndex.size() > 2)
break ;
}
// Si deve operare la modifica dei triangoli
if ( SharedIndex.size() > 2) {
// verifico che i due lai adiacenti siano controversi
int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ;
if ( nProd != 1 && nProd != - 2 && nProd != 4 && nProd != 0) {
InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetP( SharedIndex[0],
InterTria[tLB][tVLB].ptCompoVert[tCmpL]) ;
InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetP( SharedIndex[3],
InterTria[tFB][tVFB].ptCompoVert[tCmpF]) ;
InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].Validate( true) ;
InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].Validate( true) ;
bModified = true ;
break ;
}
}
}
if ( bModified)
break ;
}
}
}
}
}
}
}
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IsThereMat( int nI, int nJ, int nK) const
{
if ( nI == - 1 || nI == m_nNx[0] ||
nJ == - 1 || nJ == m_nNy[0] ||
nK == - 1 || nK == m_nNy[1])
return false ;
double dEps = 2 * EPS_SMALL ;
double dZ[3] ;
dZ[0] = ( nK + 0.5) * m_dStep ;
dZ[1] = ( nI + 0.5) * m_dStep ;
dZ[2] = ( nJ + 0.5) * m_dStep ;
int nCount = 0 ;
for ( int nGrid = 0 ; nGrid < int ( m_nMapNum) ; ++ nGrid) {
unsigned int nGrI, nGrJ ;
if ( nGrid == 0) {
nGrI = nI ;
nGrJ = nJ ;
}
else if ( nGrid == 1) {
nGrI = nJ ;
nGrJ = nK ;
}
else {
nGrI = nK ;
nGrJ = nI ;
}
unsigned int nPos = nGrJ * m_nNx[nGrid] + nGrI ;
size_t nDexSize = m_Values[nGrid][nPos].size() ;
size_t nIndex = 0 ;
while ( nIndex < nDexSize) {
if ( dZ[nGrid] > m_Values[nGrid][nPos][nIndex].dMin - dEps &&
dZ[nGrid] < m_Values[nGrid][nPos][nIndex].dMax + dEps) {
++ nCount ;
break ;
}
nIndex += 1 ;
}
}
return ( nCount == 3) ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d& ptInt) const
{
if ( nVec1[0] != nVec2[0]) {
ptInt.y = ( nVec1[1] + 0.5) * m_dStep ;
ptInt.z = ( nVec1[2] + 0.5) * m_dStep ;
int nMinI = min( nVec1[0], nVec2[0]) ;
int nMaxI = max( nVec1[0], nVec2[0]) ;
double dMinX = ( nMinI + 0.5) * m_dStep ;
double dMaxX = ( nMaxI + 0.5) * m_dStep ;
unsigned int nDexel = nVec1[2] * m_nNx[1] + nVec1[1] ;
size_t nSize = m_Values[1][nDexel].size() ;
bool bFound = false ;
for ( size_t i = 0 ; i < nSize ; i += 1) {
double dx1 = m_Values[1][nDexel][i].dMin ;
double dx2 = m_Values[1][nDexel][i].dMax ;
if ( dx1 < dMinX - EPS_SMALL && dx2 > dMinX - EPS_SMALL && dx2 < dMaxX + EPS_SMALL) {
ptInt.x = dx2 ;
bFound = true ;
break ;
}
else if ( dx1 > dMinX - EPS_SMALL && dx1 < dMaxX + EPS_SMALL && dx2 > dMaxX + EPS_SMALL) {
ptInt.x = dx1 ;
bFound = true ;
break ;
}
}
if ( ! bFound)
ptInt.x = ( dMinX + dMaxX) / 2 ;
}
else if ( nVec1[1] != nVec2[1]) {
ptInt.x = ( nVec1[0] + 0.5) * m_dStep ;
ptInt.z = ( nVec1[2] + 0.5) * m_dStep ;
int nMinJ = min( nVec1[1], nVec2[1]) ;
int nMaxJ = max( nVec1[1], nVec2[1]) ;
double dMinY = ( nMinJ + 0.5) * m_dStep ;
double dMaxY = ( nMaxJ + 0.5) * m_dStep ;
unsigned int nDexel = nVec1[0] * m_nNx[2] + nVec1[2] ;
size_t nSize = m_Values[2][nDexel].size() ;
bool bFound = false ;
for ( size_t j = 0 ; j < nSize ; j += 1) {
double dy1 = m_Values[2][nDexel][j].dMin ;
double dy2 = m_Values[2][nDexel][j].dMax ;
if ( dy1 < dMinY - EPS_SMALL && dy2 > dMinY - EPS_SMALL && dy2 < dMaxY + EPS_SMALL) {
ptInt.y = dy2 ;
bFound = true ;
break ;
}
else if ( dy1 > dMinY - EPS_SMALL && dy1 < dMaxY + EPS_SMALL && dy2 > dMaxY + EPS_SMALL) {
ptInt.y = dy1 ;
bFound = true ;
break ;
}
}
if ( ! bFound)
ptInt.y = ( dMinY + dMaxY) / 2 ;
}
else if ( nVec1[2] != nVec2[2]) {
ptInt.x = ( nVec1[0] + 0.5) * m_dStep ;
ptInt.y = ( nVec1[1] + 0.5) * m_dStep ;
int nMinK = min( nVec1[2], nVec2[2]) ;
int nMaxK = max( nVec1[2], nVec2[2]) ;
double dMinZ = ( nMinK + 0.5) * m_dStep ;
double dMaxZ = ( nMaxK + 0.5) * m_dStep ;
unsigned int nDexel = nVec1[1] * m_nNx[0] + nVec1[0] ;
size_t nSize = m_Values[0][nDexel].size() ;
bool bFound = false ;
for ( size_t k = 0 ; k < nSize ; k += 1) {
double dz1 = m_Values[0][nDexel][k].dMin ;
double dz2 = m_Values[0][nDexel][k].dMax ;
if ( dz1 < dMinZ - EPS_SMALL && dz2 > dMinZ - EPS_SMALL && dz2 < dMaxZ + EPS_SMALL) {
ptInt.z = dz2 ;
bFound = true ;
break ;
}
else if ( dz1 > dMinZ - EPS_SMALL && dz1 < dMaxZ + EPS_SMALL && dz2 > dMaxZ + EPS_SMALL) {
ptInt.z = dz1 ;
bFound = true ;
break ;
}
}
if ( ! bFound)
ptInt.z = ( dMinZ + dMaxZ) / 2 ;
}
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, Vector3d& vtNormal) const
{
double dEps = 2 * EPS_SMALL ;
bool bFound = false ;
if ( nVec1[0] != nVec2[0]) {
int nMinI = min( nVec1[0], nVec2[0]) ;
int nMaxI = max( nVec1[0], nVec2[0]) ;
double dMinX = ( nMinI + 0.5) * m_dStep ;
double dMaxX = ( nMaxI + 0.5) * m_dStep ;
ptInt.y = ( nVec1[1] + 0.5) * m_dStep ;
ptInt.z = ( nVec1[2] + 0.5) * m_dStep ;
size_t nDexel = nVec1[2] * m_nNx[1] + nVec1[1] ;
size_t nSize = m_Values[1][nDexel].size() ;
if ( bFirstCorner) {
int n = int( nSize) - 1 ;
double dX = m_Values[1][nDexel][n].dMax ;
while ( n >= 0 && dX > dMinX - dEps) {
if ( dX < dMaxX + dEps) {
ptInt.x = dX ;
vtNormal = m_Values[1][nDexel][n].vtMaxN ;
bFound = true ;
break ;
}
n -= 1 ;
if ( n >= 0)
dX = m_Values[1][nDexel][n].dMax ;
}
}
else {
size_t n = 0 ;
double dX = m_Values[1][nDexel][n].dMin ;
while ( n < nSize && dX < dMaxX + dEps) {
if ( dX > dMinX - dEps) {
ptInt.x = dX ;
vtNormal = m_Values[1][nDexel][n].vtMinN ;
bFound = true ;
break ;
}
n += 1 ;
if ( n < nSize)
dX = m_Values[1][nDexel][n].dMin ;
}
}
if ( ! bFound)
ptInt.x = 0.5 * ( dMinX + dMaxX) ;
}
else if ( nVec1[1] != nVec2[1]) {
int nMinJ = min( nVec1[1], nVec2[1]) ;
int nMaxJ = max( nVec1[1], nVec2[1]) ;
double dMinY = ( nMinJ + 0.5) * m_dStep ;
double dMaxY = ( nMaxJ + 0.5) * m_dStep ;
ptInt.x = ( nVec1[0] + 0.5) * m_dStep ;
ptInt.z = ( nVec1[2] + 0.5) * m_dStep ;
size_t nDexel = nVec1[0] * m_nNx[2] + nVec1[2] ;
size_t nSize = m_Values[2][nDexel].size() ;
if ( bFirstCorner) {
int n = int( nSize) - 1 ;
double dY = m_Values[2][nDexel][n].dMax ;
while ( n >= 0 && dY > dMinY - dEps) {
if ( dY < dMaxY + dEps) {
ptInt.y = dY ;
vtNormal = m_Values[2][nDexel][n].vtMaxN ;
bFound = true ;
break ;
}
n -= 1 ;
if ( n >= 0)
dY = m_Values[2][nDexel][n].dMax ;
}
}
else {
size_t n = 0 ;
double dY = m_Values[2][nDexel][n].dMin ;
while ( n < nSize && dY < dMaxY + dEps) {
if ( dY > dMinY - dEps) {
ptInt.y = dY ;
vtNormal = m_Values[2][nDexel][n].vtMinN ;
bFound = true ;
break ;
}
n += 1 ;
if ( n < nSize)
dY = m_Values[2][nDexel][n].dMin ;
}
}
if ( ! bFound)
ptInt.y = 0.5 * ( dMinY + dMaxY) ;
}
else if ( nVec1[2] != nVec2[2]) {
int nMinK = min( nVec1[2], nVec2[2]) ;
int nMaxK = max( nVec1[2], nVec2[2]) ;
double dMinZ = ( nMinK + 0.5) * m_dStep ;
double dMaxZ = ( nMaxK + 0.5) * m_dStep ;
ptInt.x = ( nVec1[0] + 0.5) * m_dStep ;
ptInt.y = ( nVec1[1] + 0.5) * m_dStep ;
size_t nDexel = nVec1[1] * m_nNx[0] + nVec1[0] ;
size_t nSize = m_Values[0][nDexel].size() ;
if ( bFirstCorner) {
int n = int( nSize) - 1 ;
double dZ = m_Values[0][nDexel][n].dMax ;
while ( n >= 0 && dZ > dMinZ - dEps) {
if ( dZ < dMaxZ + dEps) {
ptInt.z = dZ ;
vtNormal = m_Values[0][nDexel][n].vtMaxN ;
bFound = true ;
break ;
}
n -= 1 ;
if ( n >= 0)
dZ = m_Values[0][nDexel][n].dMax ;
}
}
else {
size_t n = 0 ;
double dZ = m_Values[0][nDexel][n].dMin ;
while ( n < nSize && dZ < dMaxZ + dEps) {
if ( dZ > dMinZ - dEps) {
ptInt.z = dZ ;
vtNormal = m_Values[0][nDexel][n].vtMinN ;
bFound = true ;
break ;
}
n += 1 ;
if ( n < nSize)
dZ = m_Values[0][nDexel][n].dMin ;
}
}
if ( ! bFound)
ptInt.z = 0.5 * ( dMinZ + dMaxZ) ;
}
return bFound ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetBlockIJKFromN( int nBlock, int nIJK[]) const
{
// Controllo sulla validità del blocco
if ( nBlock < 0 || nBlock >= int( m_nNumBlock))
return false ;
// Calcolo posizione del blocco nel reticolo
nIJK[0] = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) % int( m_nFracLin[0]) ;
nIJK[1] = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) / int( m_nFracLin[0]) ;
nIJK[2] = ( nBlock / int( m_nFracLin[0] * m_nFracLin[1])) ;
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetBlockNFromIJK( int nIJK[], int& nBlock) const
{
// Controllo sulla validità degli indici i, j, k del blocco
if ( nIJK[0] < 0 || nIJK[0] >= int( m_nFracLin[0]) ||
nIJK[1] < 0 || nIJK[1] >= int( m_nFracLin[1]) ||
nIJK[2] < 0 || nIJK[2] >= int( m_nFracLin[2]))
return false ;
// Determino il blocco
nBlock = m_nFracLin[0] * m_nFracLin[1] * nIJK[2] + m_nFracLin[0] * nIJK[1] + nIJK[0] ;
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const
{
// Controllo sulla validità degli indici i, j, k del blocco
if ( nIJK[0] < 0 || nIJK[0] >= int( m_nFracLin[0]) ||
nIJK[1] < 0 || nIJK[1] >= int( m_nFracLin[1]) ||
nIJK[2] < 0 || nIJK[2] >= int( m_nFracLin[2]))
return false ;
// Calcolo limiti per l'indice i
nLimits[0] = ( nIJK[0] == 0 ? - 1 : nIJK[0] * int( m_nDexNumPBlock)) ;
nLimits[1] = ( nIJK[0] + 1 == int( m_nFracLin[0]) ?
int( m_nNx[0]) : ( nIJK[0] + 1) * int( m_nDexNumPBlock)) ;
// Calcolo limiti per l'indice j
nLimits[2] =( nIJK[1] == 0 ? - 1 : nIJK[1] * int( m_nDexNumPBlock)) ;
nLimits[3] = ( nIJK[1] + 1 == int( m_nFracLin[1]) ?
int( m_nNy[0]) : ( nIJK[1] + 1) * int( m_nDexNumPBlock)) ;
// Calcolo limiti per l'indice k
nLimits[4] = ( nIJK[2] == 0 ? - 1 : nIJK[2] * int( m_nDexNumPBlock)) ;
nLimits[5] = ( nIJK[2] + 1 == int( m_nFracLin[2]) ?
int( m_nNy[1]) : ( nIJK[2] + 1) * int( m_nDexNumPBlock)) ;
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IsPointInsideVoxel( int nI, int nJ, int nK, const Point3d& ptP) const
{
// Controllo sull'ammissibilità del voxel
if ( nI < - 1 || nI >= int( m_nNx[0]) ||
nJ < - 1 || nJ >= int( m_nNy[0]) ||
nK < - 1 || nK >= int( m_nNy[1]))
return false ;
int nPointI = int( floor( ( ptP.x - 0.5 * m_dStep) / m_dStep)) ;
int nPointJ = int( floor( ( ptP.y - 0.5 * m_dStep) / m_dStep)) ;
int nPointK = int( floor( ( ptP.z - 0.5 * m_dStep) / m_dStep)) ;
return ( nPointI == nI &&
nPointJ == nJ &&
nPointK == nK) ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP, double dPrec) const
{
// Controllo sull'ammissibilità del voxel
if ( nI < - 1 || nI >= int( m_nNx[0]) ||
nJ < - 1 || nJ >= int( m_nNy[0]) ||
nK < - 1 || nK >= int( m_nNy[1]))
return false ;
if ( dPrec < EPS_ZERO)
dPrec = 0 ;
Point3d ptPZmap = ptP ;
ptPZmap.ToLoc( m_MapFrame[0]) ;
bool bI = ptPZmap.x > ( nI + 0.5) * m_dStep - dPrec &&
ptPZmap.x < ( nI + 1.5) * m_dStep + dPrec ;
bool bJ = ptPZmap.y > ( nJ + 0.5) * m_dStep - dPrec &&
ptPZmap.y < ( nJ + 1.5) * m_dStep + dPrec ;
bool bK = ptPZmap.z > ( nK + 0.5) * m_dStep - dPrec &&
ptPZmap.z < ( nK + 1.5) * m_dStep + dPrec ;
return ( bI && bJ && bK) ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetPointVoxel( const Point3d& ptP, int& nVoxI, int& nVoxJ, int& nVoxK) const
{
nVoxI = int( floor( ( ptP.x - 0.5 * m_dStep) / m_dStep)) ;
nVoxJ = int( floor( ( ptP.y - 0.5 * m_dStep) / m_dStep)) ;
nVoxK = int( floor( ( ptP.z - 0.5 * m_dStep) / m_dStep)) ;
return ( nVoxI >= -1 && nVoxI < int( m_nNx[0])) &&
( nVoxJ >= -1 && nVoxJ < int( m_nNy[0])) &&
( nVoxK >= -1 && nVoxK < int( m_nNy[1])) ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetVoxelBlockIJK( const int nVoxIJK[], int nBlockIJK[]) const
{
// Controllo sull'ammissibilità del voxel
if ( nVoxIJK[0] <= -2 || nVoxIJK[0] >= int( m_nNx[0]) ||
nVoxIJK[1] <= -2 || nVoxIJK[1] >= int( m_nNy[0]) ||
nVoxIJK[2] <= -2 || nVoxIJK[2] >= int( m_nNy[1]))
return false ;
// Divisioni intere
int nIntRatio0 = nVoxIJK[0] / m_nDexNumPBlock ;
int nIntRatio1 = nVoxIJK[1] / m_nDexNumPBlock ;
int nIntRatio2 = nVoxIJK[2] / m_nDexNumPBlock ;
// Resti delle divisioni intere
/*int nMod0 = nVoxIJK[0] % m_nDexNumPBlock ;
int nMod1 = nVoxIJK[1] % m_nDexNumPBlock ;
int nMod2 = nVoxIJK[2] % m_nDexNumPBlock ;*/
// Calcolo indici del blocco
nBlockIJK[0] = ( nVoxIJK[0] == -1 ? 0 : max( 0, nIntRatio0 - ( nIntRatio0 == m_nFracLin[0] ? 1 : 0))) ;
nBlockIJK[1] = ( nVoxIJK[1] == -1 ? 0 : max( 0, nIntRatio1 - ( nIntRatio1 == m_nFracLin[1] ? 1 : 0))) ;
nBlockIJK[2] = ( nVoxIJK[2] == -1 ? 0 : max( 0, nIntRatio2 - ( nIntRatio2 == m_nFracLin[2] ? 1 : 0))) ;
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetAdjBlockToBlock( int nBlockN, int nDeltaI, int nDeltaJ, int nDeltaK, int& nAdjBlockN) const
{
// Test sulla validità degli incrementi su i,j,k
if ( nDeltaI < - 1 || nDeltaI > 1 ||
nDeltaJ < - 1 || nDeltaJ > 1 ||
nDeltaK < - 1 || nDeltaK > 1)
return false ;
// Determino blocco adiacente
nAdjBlockN = nBlockN ;
nAdjBlockN += nDeltaI ;
nAdjBlockN += nDeltaJ * m_nFracLin[0] ;
nAdjBlockN += nDeltaK * m_nFracLin[0] * m_nFracLin[1] ;
// Se il blocco adiacente esiste restituisco vero, altrimenti falso.
return ( nAdjBlockN > -1 && nAdjBlockN < int( m_nNumBlock)) ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], bool bType) const
{
// Test sulla validità dei limiti
if ( nLimits[0] < -1 || nLimits[0] > int( m_nNx[0]) ||
nLimits[1] < -1 || nLimits[1] > int( m_nNx[0]) ||
nLimits[2] < -1 || nLimits[2] > int( m_nNy[0]) ||
nLimits[3] < -1 || nLimits[3] > int( m_nNy[0]) ||
nLimits[4] < -1 || nLimits[4] > int( m_nNy[1]) ||
nLimits[5] < -1 || nLimits[5] > int( m_nNy[1]))
return false ;
// Controllo sull'ammissibilità del voxel
if ( nIJK[0] <= -2 || nIJK[0] >= int( m_nNx[0]) ||
nIJK[1] <= -2 || nIJK[1] >= int( m_nNx[1]) ||
nIJK[2] <= -2 || nIJK[2] >= int( m_nNy[1]))
return false ;
// Se bType è vero cerchiamo i voxel che
// confinano con quelli di atri blocchi,
if ( bType) {
int nCurrentBlockIJK[3] ;
GetVoxelBlockIJK( nIJK, nCurrentBlockIJK) ;
// Condizione necessaria è che il voxel sia sulla frontiera
if ( IsAVoxelOnBoundary( nLimits, nIJK, false)) {
// Ciclo sulle posizioni possibili dei voxel adiacenti
for ( int nDeltaI = -1 ; nDeltaI <= 1 ; ++ nDeltaI) {
for ( int nDeltaJ = -1 ; nDeltaJ <= 1 ; ++ nDeltaJ) {
for ( int nDeltaK = -1 ; nDeltaK <= 1 ; ++ nDeltaK) {
if ( nDeltaI == 0 && nDeltaJ == 0 && nDeltaK == 0)
continue ;
// Indici dei voxel adiacenti
int nAdjIJK[3] = { nIJK[0] + nDeltaI, nIJK[1] + nDeltaJ, nIJK[2] + nDeltaK} ;
int nAdjBlockIJK[3] ;
// Determino (se esiste) il blocco in cui cade il voxel adiacente.
if ( GetVoxelBlockIJK( nAdjIJK, nAdjBlockIJK)) {
if ( nAdjBlockIJK[0] > -1 && nAdjBlockIJK[0] < int( m_nFracLin[0]) &&
nAdjBlockIJK[1] > -1 && nAdjBlockIJK[1] < int( m_nFracLin[1]) &&
nAdjBlockIJK[2] > -1 && nAdjBlockIJK[2] < int( m_nFracLin[2]) &&
( nCurrentBlockIJK[0] != nAdjBlockIJK[0] ||
nCurrentBlockIJK[1] != nAdjBlockIJK[1] ||
nCurrentBlockIJK[2] != nAdjBlockIJK[2])) {
return true ;
}
}
}
}
}
}
}
// altrimenti cerchiamo i voxel che
// sono sulla frontiera del blocco.
else {
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)
return true ;
}
return false ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert,
const int nBlockLimits[], const int nVoxIJK[]) const
{
// Se il triangolo ha area nulla non proseguiamo
if ( trTria.GetArea())
return false ;
// Se il voxel sta sulla frontiera contino,
if ( IsAVoxelOnBoundary( nBlockLimits, nVoxIJK, true)) {
double dSqEps = EPS_SMALL * EPS_SMALL ;
Point3d ptFirstGrPt, ptSecondGrPt ;
int nCount = 0 ;
if ( SqDist( trTria.GetP( 0), ptVert) > dSqEps) {
ptFirstGrPt = trTria.GetP( 0) ;
nCount ++ ;
}
if ( SqDist( trTria.GetP( 1), ptVert) > dSqEps) {
if ( nCount == 0) {
ptFirstGrPt = trTria.GetP( 1) ;
nCount ++ ;
}
else if ( nCount == 1) {
ptSecondGrPt = trTria.GetP( 1) ;
nCount ++ ;
}
}
if ( SqDist( trTria.GetP( 2), ptVert) > dSqEps) {
if ( nCount == 1) {
ptSecondGrPt = trTria.GetP( 2) ;
nCount ++ ;
}
}
/* Altro modo:
std::vector<Point3d> vPt ;
for ( int nC = 0 ; nC < 3 ; ++ nC)
if ( SqDist( trTria.GetP( 0), ptVert) > dSqEps)
vPt.emplace_back( trTria.GetP( 0)) ;
if ( vPt.size() == 2)
......*/
if ( nVoxIJK[0] == nBlockLimits[0]) {
double dXGrid = ( nBlockLimits[0] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.x - dXGrid) < EPS_SMALL &&
abs( ptSecondGrPt.x - dXGrid) < EPS_SMALL)
return true ;
}
if ( nVoxIJK[1] == nBlockLimits[2]) {
double dYGrid = ( nBlockLimits[2] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.y - dYGrid) < EPS_SMALL &&
abs( ptSecondGrPt.y - dYGrid) < EPS_SMALL)
return true ;
}
if ( nVoxIJK[2] == nBlockLimits[4]) {
double dZGrid = ( nBlockLimits[4] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.z - dZGrid) < EPS_SMALL &&
abs( ptSecondGrPt.z - dZGrid) < EPS_SMALL)
return true ;
}
if ( nVoxIJK[0] + 1 == nBlockLimits[1]) {
double dXGrid = ( nBlockLimits[1] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.x - dXGrid) < EPS_SMALL &&
abs( ptSecondGrPt.x - dXGrid) < EPS_SMALL)
return true ;
}
if ( nVoxIJK[1] + 1 == nBlockLimits[3]) {
double dYGrid = ( nBlockLimits[3] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.y - dYGrid) < EPS_SMALL &&
abs( ptSecondGrPt.y - dYGrid) < EPS_SMALL)
return true ;
}
if ( nVoxIJK[2] + 1 == nBlockLimits[5]) {
double dZGrid = ( nBlockLimits[5] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.z - dZGrid) < EPS_SMALL &&
abs( ptSecondGrPt.z - dZGrid) < EPS_SMALL)
return true ;
} return true ;
/*Altro modo
for ( int nC = 0 ; nC < 3 ; ++ nC) {
if ( nVoxIJK[nC] == nBlockLimits[2*nC]) {
double dGrid = ( nBlockLimits[2*nC] + 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] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL &&
abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL)
return true ;
}
}*/
}
// altrimenti non può starci il triangolo
else
return false ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert,
const int nBlockLimits[], const int nVoxIJK[], bool bBorderBox) const
{
// Se il voxel è di frontiera continuo,
if ( bBorderBox) {
// Se il triangolo ha area nulla non proseguiamo
if ( trTria.GetArea() < EPS_SMALL)
return false ;
// Determino i punti del triangolo sulla griglia
double dSqEps = EPS_SMALL * EPS_SMALL ;
int nNotVert[2] ;
int nCount = 0 ;
for ( int nV = 0 ; nV < 3 ; ++ nV) {
if ( SqDist( trTria.GetP( nV), ptVert) > dSqEps) {
if ( nCount == 0) {
nNotVert[nCount] = nV ;
}
else if ( nCount == 1) {
nNotVert[nCount] = nV ;
}
nCount ++ ;
}
}
// Se nCount != 2 deve esserci un errore
if ( nCount != 2)
return false ;
// Punti del triangolo sulla griglia
Point3d ptFirstGrPt = trTria.GetP( nNotVert[0]) ;
Point3d ptSecondGrPt = trTria.GetP( nNotVert[1]) ;
// Verifico se tali punti sono sula griglia
for ( int nC = 0 ; nC < 3 ; ++ nC) {
if ( nVoxIJK[nC] == nBlockLimits[2*nC]) {
double dGrid = ( nBlockLimits[2*nC] + 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] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL &&
abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL)
return true ;
}
}
return false ;
}
// altrimenti non ha senso continuare
else
return false ;
}
Reference in New Issue View Git Blame Copy Permalink