Dario Sassi
e7b2c55cbb
- inserite le modifiche di Lorenzo a Zmap - aggiunta a Zmap la possibilità di visualizzare le normali.
2511 lines
94 KiB
C++
2511 lines
94 KiB
C++
//----------------------------------------------------------------------------
|
|
// EgalTech 2015-2016
|
|
//----------------------------------------------------------------------------
|
|
// File : VolZmap.cpp Data : 22.01.15 Versione : 1.6a4
|
|
// Contenuto : Implementazione della classe Volume Zmap (tre griglie)
|
|
//
|
|
//
|
|
//
|
|
// Modifiche : 22.01.15 DS Creazione modulo.
|
|
//
|
|
//
|
|
//----------------------------------------------------------------------------
|
|
|
|
//--------------------------- Include ----------------------------------------
|
|
|
|
#include "stdafx.h"
|
|
#include "CurveLine.h"
|
|
#include "VolZmap.h"
|
|
#include "GeoConst.h"
|
|
#include "IntersLineSurfTm.h"
|
|
#include "MC_Tables.h"
|
|
#include "\EgtDev\Include\EGkIntervals.h"
|
|
#include "\EgtDev\Include\EgtNumUtils.h"
|
|
#include "\EgtDev\Include\EGkStringUtils3d.h"
|
|
#include "\EgtDev\Extern\Eigen\Core"
|
|
#include "\EgtDev\Extern\Eigen\SVD"
|
|
|
|
|
|
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 { NoFeature = 0, Corner = 1, Edge = 2} ;
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
TestOnNormal( const VectorField CompoVert[], int nCompoElem, int& FeatureType)
|
|
{
|
|
int nI, nJ ;
|
|
double dMinCosTheta = 1.001 ;
|
|
const double dCosThetaSharp = 0.9 ;
|
|
|
|
for ( int i = 0 ; i < nCompoElem ; ++ i) {
|
|
|
|
for ( int j = i + 1 ; j < nCompoElem ; ++ j) {
|
|
|
|
double dCurrentCos = CompoVert[i].vtNorm * CompoVert[j].vtNorm ;
|
|
|
|
if ( dCurrentCos < dMinCosTheta) {
|
|
|
|
nI = i ;
|
|
nJ = j ;
|
|
dMinCosTheta = dCurrentCos ;
|
|
}
|
|
}
|
|
}
|
|
|
|
if ( dMinCosTheta >= dCosThetaSharp) {
|
|
|
|
FeatureType = NoFeature ;
|
|
return false ;
|
|
}
|
|
|
|
Vector3d vtK = CompoVert[nI].vtNorm ^ CompoVert[nJ].vtNorm ;
|
|
|
|
const double dCosPhiCorner = 0.5 ;
|
|
|
|
for ( int i = 0 ; i < nCompoElem ; ++ i) {
|
|
|
|
double dAbsCurrentCos = abs( CompoVert[i].vtNorm * vtK) ;
|
|
|
|
if ( dAbsCurrentCos > dCosPhiCorner) {
|
|
|
|
FeatureType = Corner ;
|
|
return true ;
|
|
}
|
|
}
|
|
|
|
FeatureType = Edge ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
TestOnSlice( const VectorField CompoVert[], int nLimArrey,
|
|
int i,
|
|
int j,
|
|
int k, double dStep)
|
|
{
|
|
// Coordinate dei piani
|
|
double dX0 = ( i + 0.5) * dStep ;
|
|
double dX1 = ( i + 1.5) * dStep ;
|
|
double dY0 = ( j + 0.5) * dStep ;
|
|
double dY1 = ( j + 1.5) * dStep ;
|
|
double dZ0 = ( k + 0.5) * dStep ;
|
|
double dZ1 = ( k + 1.5) * dStep ;
|
|
|
|
|
|
// Analisi faccia 1
|
|
int nCount = 0 ;
|
|
double dDelta = abs( CompoVert[0].ptInt.x - dX0) ;
|
|
|
|
while ( dDelta < EPS_SMALL && nCount < nLimArrey) {
|
|
|
|
++ nCount ;
|
|
dDelta = abs( CompoVert[nCount].ptInt.x - dX0) ;
|
|
}
|
|
|
|
if ( nCount == nLimArrey)
|
|
return false ;
|
|
|
|
// Analisi faccia 2
|
|
nCount = 0 ;
|
|
dDelta = abs( CompoVert[0].ptInt.x - dX1) ;
|
|
|
|
while ( dDelta < EPS_SMALL && nCount < nLimArrey) {
|
|
|
|
++ nCount ;
|
|
dDelta = abs( CompoVert[nCount].ptInt.x - dX1) ;
|
|
}
|
|
|
|
if ( nCount == nLimArrey)
|
|
return false ;
|
|
|
|
// Analisi faccia 3
|
|
nCount = 0 ;
|
|
dDelta = abs( CompoVert[0].ptInt.y - dY0) ;
|
|
|
|
while ( dDelta < EPS_SMALL && nCount < nLimArrey) {
|
|
|
|
++ nCount ;
|
|
dDelta = abs( CompoVert[nCount].ptInt.y - dY0) ;
|
|
}
|
|
|
|
if ( nCount == nLimArrey)
|
|
return false ;
|
|
|
|
// Analisi faccia 4
|
|
nCount = 0 ;
|
|
dDelta = abs( CompoVert[0].ptInt.y - dY1) ;
|
|
|
|
while ( dDelta < EPS_SMALL && nCount < nLimArrey) {
|
|
|
|
++ nCount ;
|
|
dDelta = abs( CompoVert[nCount].ptInt.y - dY1) ;
|
|
}
|
|
|
|
if ( nCount == nLimArrey)
|
|
return false ;
|
|
|
|
// Analisi faccia 5
|
|
nCount = 0 ;
|
|
dDelta = abs( CompoVert[0].ptInt.z - dZ0) ;
|
|
|
|
while ( dDelta < EPS_SMALL && nCount < nLimArrey) {
|
|
|
|
++ nCount ;
|
|
dDelta = abs( CompoVert[nCount].ptInt.z - dZ0) ;
|
|
}
|
|
|
|
if ( nCount == nLimArrey)
|
|
return false ;
|
|
|
|
// Analisi faccia 6
|
|
nCount = 0 ;
|
|
dDelta = abs( CompoVert[0].ptInt.z - dZ1) ;
|
|
|
|
while ( dDelta < EPS_SMALL && nCount < nLimArrey) {
|
|
|
|
++ nCount ;
|
|
dDelta = abs( CompoVert[nCount].ptInt.z - dZ1) ;
|
|
}
|
|
|
|
if ( nCount == nLimArrey)
|
|
return false ;
|
|
|
|
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 = 1 ; j < int( m_Values[nDir][nPos].size()) ; j += 2) {
|
|
// aggiungo polilinea a lista
|
|
lstPL.emplace_back() ;
|
|
// inserisco punti estremi
|
|
lstPL.back().AddUPoint( 0, ptP + m_Values[nDir][nPos][j-1].dZVal * m_MapFrame[nDir].VersZ()) ;
|
|
lstPL.back().AddUPoint( 1, ptP + m_Values[nDir][nPos][j].dZVal * 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 = 1 ; j < int( m_Values[nDir][nPos].size()) ; j += 2) {
|
|
// aggiungo polilinea a lista
|
|
lstPL.emplace_back() ;
|
|
// calcolo e inserisco punto inizio spillone
|
|
Point3d ptQ = ptP + m_Values[nDir][nPos][j-1].dZVal * m_MapFrame[nDir].VersZ() ;
|
|
lstPL.back().AddUPoint( 0, ptQ) ;
|
|
// calcolo e inserisco punto su termine sua normale
|
|
Vector3d vtV = m_Values[nDir][nPos][j-1].vtN ;
|
|
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].dZVal * m_MapFrame[nDir].VersZ() ;
|
|
lstPL.back().AddUPoint( 0, ptR) ;
|
|
// calcolo e inserisco punto su termine sua normale
|
|
Vector3d vtW = m_Values[nDir][nPos][j].vtN ;
|
|
vtW.ToGlob( m_MapFrame[0]) ;
|
|
lstPL.back().AddUPoint( 1, ptR + vtW * m_dStep / 4) ;
|
|
}
|
|
return true ;
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetAllTriangles( TRIA3DLIST& lstTria) const
|
|
{
|
|
if ( m_nMapNum == 1) {
|
|
const int MAX_DIM_CHUNK = 128 ;
|
|
for ( int i = 0 ; i < int( m_nNx[0]) ; i += MAX_DIM_CHUNK) {
|
|
int nDimChunkX = min( MAX_DIM_CHUNK, int( m_nNx[0]) - i) ;
|
|
for ( int j = 0 ; j < int( m_nNy[0]) ; j += MAX_DIM_CHUNK) {
|
|
int nDimChunkY = min( MAX_DIM_CHUNK, int( m_nNy[0]) - j) ;
|
|
GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, MAX_DIM_CHUNK, lstTria) ;
|
|
}
|
|
}
|
|
}
|
|
//else {
|
|
//
|
|
// //std::vector <TriHolder> vecTria ;
|
|
// //vecTria.resize( int( m_BlockToUpdate.size())) ;
|
|
|
|
// //for ( int i = 0 ; i < int( m_BlockToUpdate.size()) ; ++ i) {
|
|
//
|
|
// //if ( m_BlockToUpdate[i])
|
|
// // ExtMarchingCubes( i, lstTria, vecTria[i]) ; }
|
|
|
|
// TriHolder triHold ;
|
|
|
|
// ExtMarchingCubes( 0, lstTria, triHold) ;
|
|
// FlipEdges( triHold) ;
|
|
//
|
|
// for ( int i = 0 ; i < int( triHold.size()) ; ++ i)
|
|
// for ( int j = 0 ; j < int( triHold[i].vecTria.size()) ; ++ j)
|
|
|
|
// lstTria.emplace_back( triHold[i].vecTria[j]) ;
|
|
//}
|
|
else
|
|
MarchingCubes( lstTria) ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetBlockTriangles( int nBlock, TRIA3DLIST& lstTria) const
|
|
{
|
|
// Controllo sulla validità del blocco
|
|
if ( nBlock < 0 || nBlock >= int( m_BlockToUpdate.size()))
|
|
return false ;
|
|
|
|
// Caso di singola mappa
|
|
if ( m_nMapNum == 1) {
|
|
|
|
const int MAX_DIM_CHUNK = 128 ;
|
|
|
|
// Calcolo posizione del blocco nella griglia
|
|
int nIBlock = nBlock % int( m_nFracLin[0]) ;
|
|
int nJBlock = nBlock / 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, lstTria) ;
|
|
}
|
|
}
|
|
}
|
|
// Caso con tre mappe
|
|
else {
|
|
|
|
// Calcolo i limiti sugli indici dei voxel del blocco
|
|
// Vettore indici i,j,k del blocco
|
|
int nIJK[3] ;
|
|
GetBlockIJK( nIJK, nBlock) ;
|
|
// Vettore limiti sugli indici dei voxel del blocco
|
|
int LimitIndexes[6] ;
|
|
GetBlockLimitsIJK( nIJK, LimitIndexes) ;
|
|
|
|
TriHolder triHold ;
|
|
|
|
ExtMarchingCubes( LimitIndexes, lstTria, triHold) ;
|
|
FlipEdges( triHold) ;
|
|
|
|
// Valuto se esistono voxel contenenti sharp feature
|
|
if ( ! triHold.size())
|
|
|
|
;
|
|
|
|
else {
|
|
|
|
// Cerco fra i voxel con sharp feature quelli
|
|
// di frontiera e quelli ad essi adiacenti.
|
|
|
|
// Ciclo sui voxel con sharp feature
|
|
for ( size_t t = 0 ; t < triHold.size() ; ++ t) {
|
|
|
|
// Voxel di frontiera
|
|
if ( triHold[t].i == LimitIndexes[0] || triHold[t].i == LimitIndexes[1] - 1 ||
|
|
triHold[t].j == LimitIndexes[2] || triHold[t].j == LimitIndexes[3] - 1 ||
|
|
triHold[t].k == LimitIndexes[4] || triHold[t].k == LimitIndexes[5] - 1) {
|
|
|
|
// Ciclo sui voxel adiacenti
|
|
for ( int nI = triHold[t].i - 1 ; nI <= triHold[t].i + 1 ; ++ nI) {
|
|
for ( int nJ = triHold[t].j - 1 ; nJ <= triHold[t].j + 1 ; ++ nJ) {
|
|
for ( int nK = triHold[t].k - 1 ; nK <= triHold[t].k + 1 ; ++ nK) {
|
|
|
|
// Voxel adiacente appartenente allo
|
|
// Zmap e non appartenente al blocco
|
|
if ( ( nI >= -1 && nI <= int( m_nNx[0]) - 1 &&
|
|
nJ >= -1 && nJ <= int( m_nNy[0]) - 1 &&
|
|
nK >= -1 && nK <= int( m_nNy[1]) - 1) &&
|
|
! ( nI >= LimitIndexes[0] && nI < LimitIndexes[1] &&
|
|
nJ >= LimitIndexes[2] && nJ < LimitIndexes[3] &&
|
|
nK >= LimitIndexes[4] && nK < LimitIndexes[5])) {
|
|
|
|
std::vector <int> AdjVoxel ;
|
|
|
|
AdjVoxel.emplace_back( nI) ;
|
|
AdjVoxel.emplace_back( nJ) ;
|
|
AdjVoxel.emplace_back( nK) ;
|
|
|
|
TriHolder NewTriHold ;
|
|
|
|
ExtMarchingCubes( AdjVoxel, NewTriHold) ;
|
|
|
|
if ( NewTriHold.size())
|
|
FlipEdgesLocalFlipEdges( triHold[t], NewTriHold[0]) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
for ( size_t t1 = 0 ; t1 < triHold.size() ; ++ t1)
|
|
for ( size_t t2 = 0 ; t2 < triHold[t1].vCompoTria.size() ; ++ t2)
|
|
for ( size_t t3 = 0 ; t3 < triHold[t1].vCompoTria[t2].size() ; ++ t3)
|
|
lstTria.emplace_back( triHold[t1].vCompoTria[t2][t3]) ;
|
|
|
|
}
|
|
|
|
m_BlockToUpdate[nBlock] = false ;
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetBlockInfo( std::vector<bool> & bModified) const
|
|
{
|
|
bModified = m_BlockToUpdate ;
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int
|
|
VolZmap::GetBlockCount( void) const
|
|
{
|
|
return m_nNumBlock ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
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()) != 2)
|
|
bIsSimple = false ;
|
|
else if ( i == 0 && j == 0) {
|
|
dBotZ = m_Values[0][nPos][0].dZVal ;
|
|
dTopZ = m_Values[0][nPos][1].dZVal ;
|
|
}
|
|
else if ( abs( m_Values[0][nPos][0].dZVal - dBotZ) > EPS_SMALL ||
|
|
abs( m_Values[0][nPos][1].dZVal - 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][1].dZVal * m_MapFrame[0].VersZ() ;
|
|
Vector3d vtDZb = m_Values[0][nPos][0].dZVal * 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 = 1 ; i < int( m_Values[0][nPos].size()) ; i += 2) {
|
|
Vector3d vtDZt = m_Values[0][nPos][i].dZVal * m_MapFrame[0].VersZ() ;
|
|
Vector3d vtDZb = m_Values[0][nPos][i-1].dZVal * 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 = 1 ; i < int( m_Values[0][nPos].size()) ; i += 2)
|
|
intFace.Add( m_Values[0][nPos][i-1].dZVal, m_Values[0][nPos][i].dZVal) ;
|
|
if ( nPosAdj > 0) {
|
|
for ( int i = 1 ; i < int( m_Values[0][nPosAdj].size()) ; i += 2)
|
|
intFace.Subtract( m_Values[0][nPosAdj][i-1].dZVal, m_Values[0][nPosAdj][i].dZVal) ;
|
|
}
|
|
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( TRIA3DLIST& lstTria) const
|
|
{
|
|
// Ciclo su tutti i voxel dello Zmap
|
|
for ( int i = - 1 ; i < int( m_nNx[0]) ; ++ i) {
|
|
for ( int j = - 1 ; j < int( m_nNy[0]) ; ++ j) {
|
|
for ( int k = - 1 ; k < int( m_nNy[1]) ; ++ 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 }
|
|
} ;
|
|
|
|
Point3d ptIntPoint[12] ;
|
|
|
|
// 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] ;
|
|
|
|
// 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] ;
|
|
|
|
// Il triangolo è pronto
|
|
Triangle3d CurrentTriangle ;
|
|
CurrentTriangle.Set( ptIntPoint[i0], ptIntPoint[i1], ptIntPoint[i2]) ;
|
|
CurrentTriangle.Validate() ;
|
|
|
|
// Aggiungo triangolo
|
|
lstTria.emplace_back( CurrentTriangle) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const
|
|
{
|
|
if ( nBlock < 0 || nBlock >= int( m_BlockToUpdate.size()))
|
|
return false ;
|
|
|
|
Point3d ptMapOrig = m_MapFrame[0].Orig() ;
|
|
|
|
// 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( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& triHold) const
|
|
{
|
|
Point3d ptMapOrig = m_MapFrame[0].Orig() ;
|
|
|
|
// 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) {
|
|
|
|
// 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}
|
|
} ;
|
|
|
|
int nIndex = 0 ;
|
|
|
|
// Classificazione dei vertici: interni o esterni al materiale
|
|
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 }
|
|
} ;
|
|
|
|
// Arrey di strutture punto di intersezione
|
|
// e normale alla superficie in esso.
|
|
VectorField VecField[12] ;
|
|
|
|
// 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.
|
|
NewIntersPos( IndexCorner[n1], IndexCorner[n2], bN1,
|
|
VecField[EdgeIndex].ptInt,
|
|
VecField[EdgeIndex].vtNorm) ;
|
|
|
|
VecField[EdgeIndex].ptInt.ToGlob( m_MapFrame[0]) ;
|
|
VecField[EdgeIndex].vtNorm.ToGlob( m_MapFrame[0]) ;
|
|
}
|
|
|
|
// Determino il numero di componenti connesse
|
|
int nComponents = TriangleTableEn[nIndex][1][0] ;
|
|
|
|
// Serve nel ciclo che salva i punti e vettori di
|
|
// una componente nell'arrey di compentenza: La tabella
|
|
// fornisce numero di componenti, numero di vertici per
|
|
// componenti per OGNUNA delle componenti e in fine
|
|
// elenca i vertici della prima componente, seguiti da quelli
|
|
// della seconda e così via.
|
|
int nTableOffset = nComponents ;
|
|
|
|
// Numero di feature nel voxel: al più vi è una feature per componente connessa
|
|
int nFeatureInVoxel = 0 ;
|
|
|
|
// Ciclo sulle componenti
|
|
for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) {
|
|
|
|
// Numero vertici per componenti
|
|
int nVertComp = TriangleTableEn[nIndex][1][nCompCount] ;
|
|
|
|
// Vettore di Vector3d
|
|
VectorField CompoVert[12] ;
|
|
|
|
// Riempio il vettore
|
|
for ( int nVertCount = 0 ; nVertCount < nVertComp ; ++ nVertCount)
|
|
// Nota che il primo elemento dell'array (0-esimo) non viene inizializzato
|
|
CompoVert[nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nTableOffset + 1]] ;
|
|
|
|
int nFeatureType ;
|
|
|
|
// Valuto le relazioni reciproche fra le normali e
|
|
// se i punti sono su un piano di griglia.
|
|
bool bNormal = TestOnNormal( CompoVert, nVertComp, nFeatureType) ;
|
|
bool bExt = bNormal ;
|
|
|
|
// Extended MC
|
|
if ( bExt) {
|
|
|
|
// Passo al sistema di riferimento del baricentro
|
|
Point3d ptGravityCenter( 0, 0, 0) ;
|
|
|
|
for ( int ni = 0 ; ni < nVertComp ; ++ ni)
|
|
ptGravityCenter += CompoVert[ni].ptInt ;
|
|
|
|
ptGravityCenter /= nVertComp ;
|
|
|
|
Vector3d vtO = ptGravityCenter - ORIG ;
|
|
|
|
Point3d ptTrasf[12] ;
|
|
|
|
for ( int ni = 0 ; ni < nVertComp ; ++ ni)
|
|
ptTrasf[ni] = CompoVert[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, 3) ;
|
|
dKnownVector.resize( nVertComp, 1) ;
|
|
dUnknownVector.resize( 3, 1) ;
|
|
|
|
// Studio del caso 4 punti su un piano
|
|
int nEqual = 0 ;
|
|
int nPosD ;
|
|
Vector3d vtD, vtE ;
|
|
|
|
if ( nVertComp == 4 && nFeatureType == 2) {
|
|
int nPosEq ;
|
|
for ( int ni = 0 ; ni < 2 ; ++ ni) {
|
|
for ( int nj = ni + 1 ; nj < nVertComp ; ++ nj) {
|
|
|
|
if ( AreSameVectorApprox( CompoVert[ni].vtNorm,
|
|
CompoVert[nj].vtNorm)) {
|
|
nEqual ++ ;
|
|
nPosEq = ni ;
|
|
}
|
|
}
|
|
if ( nEqual == 2)
|
|
break ;
|
|
else
|
|
nEqual = 0 ;
|
|
}
|
|
if ( nEqual == 2) {
|
|
|
|
for ( int ni = 0 ; ni < nVertComp ; ++ ni)
|
|
|
|
if ( ! AreSameVectorApprox( CompoVert[ni].vtNorm,
|
|
CompoVert[nPosEq].vtNorm)) {
|
|
nPosD = ni ;
|
|
vtD = CompoVert[ni].vtNorm ;
|
|
vtE = CompoVert[nPosEq].vtNorm ;
|
|
}
|
|
}
|
|
}
|
|
|
|
double dDot = abs( ( CompoVert[1].ptInt - CompoVert[0].ptInt) *
|
|
( ( CompoVert[2].ptInt - CompoVert[1].ptInt) ^
|
|
( CompoVert[3].ptInt - CompoVert[2].ptInt))) ;
|
|
|
|
// Caso quattro punti su un piano
|
|
if ( nVertComp == 4 && nEqual == 2 && dDot < EPS_SMALL) {
|
|
|
|
for ( int ni = 0 ; ni < nVertComp ; ++ ni) {
|
|
if ( ni != nPosD) {
|
|
dMatrixN( ni, 0) = CompoVert[ni].vtNorm.x ;
|
|
dMatrixN( ni, 1) = CompoVert[ni].vtNorm.y ;
|
|
dMatrixN( ni, 2) = CompoVert[ni].vtNorm.z ;
|
|
dKnownVector( ni) = CompoVert[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 ; ++ ni) {
|
|
dMatrixN( ni, 0) = CompoVert[ni].vtNorm.x ;
|
|
dMatrixN( ni, 1) = CompoVert[ni].vtNorm.y ;
|
|
dMatrixN( ni, 2) = CompoVert[ni].vtNorm.z ;
|
|
dKnownVector( ni) = CompoVert[ni].vtNorm * ( ptTrasf[ni] - ORIG) ;
|
|
}
|
|
}
|
|
|
|
DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ;
|
|
|
|
dMatrixU = svd.matrixU() ;
|
|
dMatrixV = svd.matrixV() ;
|
|
dSingularValue = svd.singularValues() ;
|
|
|
|
double s0 = dSingularValue( 0) ;
|
|
double s1 = dSingularValue( 1) ;
|
|
double s2 = dSingularValue( 2) ;
|
|
|
|
if ( nFeatureType == 2 )
|
|
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 ( abs( dSingularValue( ni)) > EPS_SMALL) {
|
|
|
|
for ( int nj = 0 ; nj < nVertComp ; ++ 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 ;
|
|
}
|
|
|
|
// Limito la feature entro una distanza di 3
|
|
// volte la diagonale del voxel dal baricentro.
|
|
Vector3d vtFeature( dUnknownVector( 0),
|
|
dUnknownVector( 1),
|
|
dUnknownVector( 2)) ;
|
|
|
|
double dDistFeature = vtFeature.Len() ;
|
|
|
|
const double dMaxDist = sqrt( 3) * m_dStep ;
|
|
|
|
if ( dDistFeature > dMaxDist) {
|
|
|
|
// 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 ;
|
|
|
|
// Il triangolo è pronto
|
|
CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
|
|
// Se il triangolo non è degenere lo aggiungo alla lista
|
|
if ( ! ( AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) ||
|
|
AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) ||
|
|
AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0))))
|
|
lstTria.emplace_back( CurrentTriangle) ;
|
|
}
|
|
|
|
continue ;
|
|
}
|
|
|
|
// Esprimo la soluzione nel sistema di riferimento dello z-map.
|
|
Point3d ptSol( dUnknownVector( 0) + vtO.x,
|
|
dUnknownVector( 1) + vtO.y,
|
|
dUnknownVector( 2) + vtO.z) ;
|
|
|
|
Triangle3d CurrentTriangle ;
|
|
TRIA3DVECTOR triContainer ;
|
|
|
|
bool bInvNormal = false ;
|
|
|
|
for ( int ni = 0 ; ni < nVertComp ; ++ ni) {
|
|
|
|
// Il triangolo è pronto
|
|
int nj = ( ni + 1 < nVertComp) ? ni + 1 : 0 ;
|
|
CurrentTriangle.Set( ptSol, CompoVert[nj].ptInt, CompoVert[ni].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
|
|
// Test sulla degenerazione
|
|
if ( CurrentTriangle.GetArea() < EPS_SMALL)
|
|
continue ;
|
|
|
|
// Controllo sull'inversione delle normali
|
|
if ( CurrentTriangle.GetN() * CompoVert[nj].vtNorm < - 0.7 ||
|
|
CurrentTriangle.GetN() * CompoVert[ni].vtNorm < - 0.7) {
|
|
//string sInfo ;
|
|
//sInfo = "TriaN=" + ToString( CurrentTriangle.GetN()) +
|
|
// "VertJ=" + ToString( CompoVert[nj].vtNorm) +
|
|
// "VertI=" + ToString( CompoVert[ni].vtNorm) ;
|
|
//LOG_INFO( GetEGkLogger(), sInfo.c_str())
|
|
//bInvNormal = true ;
|
|
}
|
|
|
|
// Aggiungo triangolo al vettore temporaneo
|
|
triContainer.emplace_back( CurrentTriangle) ;
|
|
}
|
|
|
|
// Valuto normali
|
|
int nContSize = 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 || bInvNormal)) {
|
|
|
|
// Aggiorno il numero di feature.
|
|
++ nFeatureInVoxel ;
|
|
|
|
// Se siamo in presenza della prima feature del
|
|
// voxel, ridimensiono il vettore che contiene
|
|
// la struttura dati del voxel.
|
|
if ( nFeatureInVoxel == 1) {
|
|
|
|
triHold.resize( triHold.size() + 1) ;
|
|
|
|
// Questi dati dipendono solo dal voxel,
|
|
// quindi sono validi 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 = 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 < nContSize ; ++ ni)
|
|
|
|
triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( triContainer[ni]) ;
|
|
}
|
|
|
|
else {
|
|
// 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 ;
|
|
|
|
// Il triangolo è pronto
|
|
CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
|
|
// Se il triangolo non è degenere lo aggiungo alla lista
|
|
if ( ! ( AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) ||
|
|
AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) ||
|
|
AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0))))
|
|
lstTria.emplace_back( CurrentTriangle) ;
|
|
}
|
|
}
|
|
}
|
|
// Standard MC
|
|
else {
|
|
|
|
// 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 ;
|
|
|
|
// Il triangolo è pronto
|
|
CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
|
|
// Se il triangolo non è degenere lo aggiungo alla lista
|
|
if ( ! ( AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) ||
|
|
AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) ||
|
|
AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0))))
|
|
lstTria.emplace_back( CurrentTriangle) ;
|
|
}
|
|
}
|
|
|
|
nTableOffset += nVertComp ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::ExtMarchingCubes( const std::vector <int> VoxelsIndexes, TriHolder& triHold) const
|
|
{
|
|
// Controllo sulla validità della dimensione del
|
|
// vector: essa deve essere un multiplo di 3
|
|
if ( ( VoxelsIndexes.size() % 3))
|
|
return false ;
|
|
|
|
Point3d ptMapOrig = m_MapFrame[0].Orig() ;
|
|
|
|
// Ciclo sui voxel
|
|
for ( size_t t = 0 ; t < VoxelsIndexes.size() - 2 ; t += 3) {
|
|
|
|
// Indici del voxel
|
|
int i = VoxelsIndexes[t] ;
|
|
int j = VoxelsIndexes[t + 1] ;
|
|
int k = VoxelsIndexes[t + 2] ;
|
|
|
|
// 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}
|
|
} ;
|
|
|
|
int nIndex = 0 ;
|
|
|
|
// Classificazione dei vertici: interni o esterni al materiale
|
|
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 }, { 0, 3 }, { 4, 5 }, { 5, 6 },
|
|
{ 6, 7 }, { 4, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }
|
|
} ;
|
|
|
|
// Arrey di strutture punto di intersezione
|
|
// e normale alla superficie in esso.
|
|
VectorField VecField[12] ;
|
|
|
|
// 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] ;
|
|
|
|
// Determino con precisione il punto di intersezione sullo spigolo.
|
|
IntersPos( IndexCorner[n1], IndexCorner[n2],
|
|
VecField[EdgeIndex].ptInt,
|
|
VecField[EdgeIndex].vtNorm) ;
|
|
|
|
VecField[EdgeIndex].ptInt.ToGlob( m_MapFrame[0]) ;
|
|
VecField[EdgeIndex].vtNorm.ToGlob( m_MapFrame[0]) ;
|
|
}
|
|
|
|
// Determino il numero di componenti connesse
|
|
int nComponents = TriangleTableEn[nIndex][1][0] ;
|
|
|
|
// Serve nel ciclo che salva i punti e vettori di
|
|
// una componente nell'arrey di compentenza: La tabella
|
|
// fornisce numero di componenti, numero di vertici per
|
|
// componenti per OGNUNA delle componenti e in fine
|
|
// elenca i vertici della prima componente, seguiti da quelli
|
|
// della seconda e così via.
|
|
int nTableOffset = nComponents ;
|
|
|
|
// Numero di feature nel voxel: al più vi
|
|
// è una feature per componente connessa
|
|
int nFeatureInVoxel = 0 ;
|
|
|
|
// Ciclo sulle componenti
|
|
for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) {
|
|
|
|
// Numero vertici per componenti
|
|
int nVertComp = TriangleTableEn[nIndex][1][nCompCount] ;
|
|
|
|
// Vettore di Vector3d
|
|
VectorField CompoVert[12] ;
|
|
|
|
// Riempio il vettore
|
|
for ( int nVertCount = 0 ; nVertCount < nVertComp ; ++ nVertCount)
|
|
// Nota che il primo elemento dell'array (0-esimo) non viene inizializzato
|
|
CompoVert[nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nTableOffset + 1]] ;
|
|
|
|
int nFeatureType ;
|
|
|
|
// Valuto le relazioni reciproche fra le normali
|
|
bool bExt = TestOnNormal( CompoVert, nVertComp, nFeatureType) ;
|
|
|
|
// Extended MC
|
|
if ( bExt) {
|
|
|
|
// Aggiorno il numero di feature.
|
|
++ nFeatureInVoxel ;
|
|
|
|
// Se siamo in presenza della prima feature del
|
|
// voxel, ridimensiono il vettore che contiene
|
|
// la struttura dati del voxel.
|
|
if ( nFeatureInVoxel == 1) {
|
|
|
|
triHold.resize( triHold.size() + 1) ;
|
|
|
|
// Questi dati dipendono solo dal voxel,
|
|
// quindi sono validi 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 ;
|
|
|
|
// Passo al sistema di riferimento del baricentro
|
|
Point3d ptGravityCenter( 0, 0, 0) ;
|
|
|
|
for ( int ni = 0 ; ni < nVertComp ; ++ ni)
|
|
|
|
ptGravityCenter += CompoVert[ni].ptInt ;
|
|
|
|
ptGravityCenter /= nVertComp ;
|
|
|
|
Vector3d vtO = ptGravityCenter - ORIG ;
|
|
|
|
Point3d ptTrasf[12] ;
|
|
|
|
for ( int ni = 0 ; ni < nVertComp ; ++ ni)
|
|
|
|
ptTrasf[ni] = CompoVert[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, 3) ;
|
|
dKnownVector.resize( nVertComp, 1) ;
|
|
dUnknownVector.resize( 3, 1) ;
|
|
|
|
// Studio del caso 4 punti su un piano
|
|
int nEqual = 0 ;
|
|
int nPosD ;
|
|
Vector3d vtD, vtE ;
|
|
|
|
if ( nVertComp == 4 && nFeatureType == 2) {
|
|
|
|
int nPosEq ;
|
|
|
|
for ( int ni = 0 ; ni < 2 ; ++ ni) {
|
|
|
|
for ( int nj = ni + 1 ; nj < nVertComp ; ++ nj) {
|
|
|
|
if ( AreSameVectorApprox( CompoVert[ni].vtNorm,
|
|
CompoVert[nj].vtNorm)) {
|
|
|
|
nEqual ++ ;
|
|
nPosEq = ni ;
|
|
}
|
|
}
|
|
|
|
if ( nEqual == 2)
|
|
break ;
|
|
else
|
|
nEqual = 0 ;
|
|
}
|
|
|
|
if ( nEqual == 2) {
|
|
|
|
for ( int ni = 0 ; ni < nVertComp ; ++ ni)
|
|
|
|
if ( ! AreSameVectorApprox( CompoVert[ni].vtNorm,
|
|
CompoVert[nPosEq].vtNorm)) {
|
|
|
|
nPosD = ni ;
|
|
vtD = CompoVert[ni].vtNorm ;
|
|
vtE = CompoVert[nPosEq].vtNorm ;
|
|
}
|
|
}
|
|
}
|
|
|
|
double dDot = abs( ( CompoVert[1].ptInt - CompoVert[0].ptInt) *
|
|
( ( CompoVert[2].ptInt - CompoVert[1].ptInt) ^
|
|
( CompoVert[3].ptInt - CompoVert[2].ptInt))) ;
|
|
|
|
|
|
|
|
// Caso quattro punti su un piano
|
|
if ( nVertComp == 4 && nEqual == 2 && dDot < EPS_SMALL) {
|
|
|
|
for ( int ni = 0 ; ni < nVertComp ; ++ ni) {
|
|
|
|
if ( ni != nPosD) {
|
|
|
|
dMatrixN( ni, 0) = CompoVert[ni].vtNorm.x ;
|
|
dMatrixN( ni, 1) = CompoVert[ni].vtNorm.y ;
|
|
dMatrixN( ni, 2) = CompoVert[ni].vtNorm.z ;
|
|
|
|
dKnownVector( ni) = CompoVert[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 ; ++ ni) {
|
|
|
|
dMatrixN( ni, 0) = CompoVert[ni].vtNorm.x ;
|
|
dMatrixN( ni, 1) = CompoVert[ni].vtNorm.y ;
|
|
dMatrixN( ni, 2) = CompoVert[ni].vtNorm.z ;
|
|
|
|
dKnownVector( ni) = CompoVert[ni].vtNorm * ( ptTrasf[ni] - ORIG) ;
|
|
}
|
|
}
|
|
|
|
DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ;
|
|
|
|
dMatrixU = svd.matrixU() ;
|
|
dMatrixV = svd.matrixV() ;
|
|
dSingularValue = svd.singularValues() ;
|
|
|
|
double s0 = dSingularValue( 0) ;
|
|
double s1 = dSingularValue( 1) ;
|
|
double s2 = dSingularValue( 2) ;
|
|
|
|
if ( nFeatureType == 2 )
|
|
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 ( abs( dSingularValue( ni)) > EPS_SMALL) {
|
|
|
|
for ( int nj = 0 ; nj < nVertComp ; ++ 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 ;
|
|
}
|
|
|
|
// Limito la feature entro una distanza di 3
|
|
// volte la diagonale del voxel dal baricentro.
|
|
Vector3d vtFeature( dUnknownVector( 0),
|
|
dUnknownVector( 1),
|
|
dUnknownVector( 2)) ;
|
|
|
|
double dDistFeature = vtFeature.Len() ;
|
|
|
|
const double dMaxDist = 2 * sqrt( 3) * m_dStep / 2 ;
|
|
|
|
if ( dDistFeature > dMaxDist) {
|
|
|
|
vtFeature = ( dMaxDist / dDistFeature) * vtFeature ;
|
|
|
|
dUnknownVector( 0) = vtFeature.x ;
|
|
dUnknownVector( 1) = vtFeature.y ;
|
|
dUnknownVector( 2) = vtFeature.z ;
|
|
}
|
|
|
|
// Esprimo la soluzione nel sistema di riferimento dello z-map
|
|
Point3d ptSol( dUnknownVector( 0) + vtO.x,
|
|
dUnknownVector( 1) + vtO.y,
|
|
dUnknownVector( 2) + vtO.z) ;
|
|
|
|
// Aggiungo vertice della componente
|
|
// connessa alla lista dei vertici.
|
|
triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ;
|
|
|
|
int nOldFeatureNum = triHold[nCurrent].vCompoTria.size() ;
|
|
int nNewFeatureNum = nOldFeatureNum + 1 ;
|
|
|
|
// Aggiungo una componente per il vettore
|
|
// dei triangoli della componente connessa.
|
|
triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ;
|
|
|
|
Triangle3d CurrentTriangle ;
|
|
TRIA3DVECTOR triContainer ;
|
|
|
|
for ( int ni = 0 ; ni < nVertComp - 1 ; ++ ni) {
|
|
|
|
// Il triangolo è pronto
|
|
CurrentTriangle.Set( ptSol, CompoVert[ni+1].ptInt, CompoVert[ni].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
|
|
// Test sulla degenerazione
|
|
if ( CurrentTriangle.GetArea() < EPS_SMALL /*AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) ||
|
|
AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) ||
|
|
AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0))*/)
|
|
|
|
continue ;
|
|
|
|
// Aggiungo triangolo
|
|
triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( CurrentTriangle) ;
|
|
}
|
|
|
|
// Ultimo triangolo
|
|
CurrentTriangle.Set( ptSol, CompoVert[0].ptInt, CompoVert[nVertComp - 1].ptInt) ;
|
|
CurrentTriangle.Validate( true) ;
|
|
|
|
// Test sulla degenerazione
|
|
if ( CurrentTriangle.GetArea() < EPS_SMALL /*AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) ||
|
|
AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) ||
|
|
AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0))*/)
|
|
|
|
continue ;
|
|
|
|
// Aggiungo ultimo triangolo
|
|
triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( CurrentTriangle) ;
|
|
}
|
|
|
|
nTableOffset += nVertComp ;
|
|
}
|
|
}
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::FlipEdges( TriHolder& triHold) const
|
|
{
|
|
// Numero di voxel in cui si presentano sharp feature
|
|
int nVoxelNum = int( triHold.size()) ;
|
|
|
|
// Ciclo su tali voxel
|
|
for ( int n = 0 ; n < nVoxelNum ; ++ n) {
|
|
for ( int m = n ; m < nVoxelNum ; ++ m) {
|
|
|
|
// Voxel adiacenti o coincidenti
|
|
if ( m == n ||
|
|
( ( triHold[m].i < int( m_nNx[0]) &&
|
|
triHold[m].j < int( m_nNy[0]) &&
|
|
triHold[m].k < int( m_nNy[1])) &&
|
|
( triHold[m].i == triHold[n].i + 1 ||
|
|
triHold[m].j == triHold[n].j + 1 ||
|
|
triHold[m].k == triHold[n].k + 1))) {
|
|
|
|
// Numero delle componenti connesse nei due voxel
|
|
int nNumCompoN = triHold[n].ptCompoVert.size() ;
|
|
int nNumCompoM = triHold[m].ptCompoVert.size() ;
|
|
|
|
int nCompoN = 0 ;
|
|
|
|
// Ciclo sulle componenti
|
|
for ( ; nCompoN < nNumCompoN ; ++ nCompoN) {
|
|
|
|
int nCompoM = ( m == n ? nCompoN + 1 : 0) ;
|
|
|
|
for ( ; nCompoM < nNumCompoM ; ++ nCompoM) {
|
|
|
|
// Numero di triangoli per le componenti connesse
|
|
int nNumN = int( triHold[n].vCompoTria[nCompoN].size()) ;
|
|
int nNumM = int( triHold[m].vCompoTria[nCompoM].size()) ;
|
|
|
|
for ( int triN = 0 ; triN < nNumN ; ++ triN) {
|
|
|
|
bool bModified = false ;
|
|
|
|
for ( int triM = 0 ; triM < nNumM ; ++ triM) {
|
|
|
|
std::vector <int> SharedIndex ;
|
|
|
|
for ( int vertN = 0 ; vertN < 3 ; ++ vertN) {
|
|
|
|
for ( int vertM = 0 ; vertM < 3 ; ++ vertM) {
|
|
|
|
Point3d ptN = triHold[n].vCompoTria[nCompoN][triN].GetP( vertN) ;
|
|
Point3d ptM = triHold[m].vCompoTria[nCompoM][triM].GetP( vertM) ;
|
|
|
|
if ( SqDist( ptN, ptM) < EPS_SMALL * EPS_SMALL) {
|
|
|
|
Point3d ptVertN = triHold[n].ptCompoVert[nCompoN] ;
|
|
Point3d ptVertM = triHold[m].ptCompoVert[nCompoM] ;
|
|
|
|
if ( ! ( AreSamePointApprox( ptN, ptVertN) ||
|
|
AreSamePointApprox( ptM, ptVertM))) {
|
|
|
|
SharedIndex.emplace_back( vertN) ;
|
|
SharedIndex.emplace_back( vertM) ;
|
|
}
|
|
}
|
|
|
|
|
|
if ( SharedIndex.size() > 2)
|
|
break ;
|
|
}
|
|
|
|
if ( SharedIndex.size() > 2)
|
|
break ;
|
|
}
|
|
// Si deve operare la modifica dei triangoli
|
|
if ( SharedIndex.size() > 2) {
|
|
|
|
// Controllo sulle normali
|
|
Point3d ptMFeature = triHold[m].ptCompoVert[nCompoM] ;
|
|
Point3d ptNFeature = triHold[n].ptCompoVert[nCompoN] ;
|
|
|
|
Vector3d vtFeature = ptMFeature - ptNFeature ;
|
|
|
|
vtFeature.Normalize() ;
|
|
|
|
double dDot = vtFeature * triHold[m].vCompoTria[nCompoM][triM].GetN() ;
|
|
|
|
// Ulteriore controllo sulle normali
|
|
Point3d ptP0 = ptNFeature ;
|
|
Point3d ptP1 = triHold[n].vCompoTria[nCompoN][triN].GetP( SharedIndex[0]) ;
|
|
Point3d ptP2 = triHold[n].vCompoTria[nCompoN][triN].GetP( SharedIndex[2]) ;
|
|
Point3d ptP3 = ptMFeature ;
|
|
|
|
double dPlane = ( ( ptP1 - ptP0) ^ ( ptP2 - ptP0)) * ( ptP3 - ptP0) ;
|
|
|
|
Vector3d vtN = triHold[n].vCompoTria[nCompoN][triN].GetN() ;
|
|
Vector3d vtM = triHold[m].vCompoTria[nCompoM][triM].GetN() ;
|
|
|
|
double dDotNM = vtN * vtM ;
|
|
|
|
int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ;
|
|
|
|
// ---
|
|
if ( nProd != 1 && nProd != - 2 && nProd != 4) {
|
|
|
|
triHold[n].vCompoTria[nCompoN][triN].SetP( SharedIndex[0],
|
|
triHold[m].ptCompoVert[nCompoM]) ;
|
|
triHold[m].vCompoTria[nCompoM][triM].SetP( SharedIndex[3],
|
|
triHold[n].ptCompoVert[nCompoN]) ;
|
|
|
|
triHold[n].vCompoTria[nCompoN][triN].Validate( true) ;
|
|
triHold[m].vCompoTria[nCompoM][triM].Validate( true) ;
|
|
|
|
bModified = true ;
|
|
break ;
|
|
}
|
|
}
|
|
}
|
|
|
|
if ( bModified)
|
|
break ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::FlipEdgesLocalFlipEdges( TriaStruct& triStrCurr, TriaStruct& triStrAdj) const
|
|
{
|
|
// Numero delle componenti connesse nei due voxel
|
|
size_t nCurVoxCompNum = triStrCurr.ptCompoVert.size() ;
|
|
size_t nAdjVoxCompNum = triStrAdj.ptCompoVert.size() ;
|
|
|
|
// Ciclo sulle componenti connesse
|
|
for ( size_t nCVCurComp = 0 ; nCVCurComp < nCurVoxCompNum ; ++ nCVCurComp) {
|
|
for ( size_t nAVCurComp = 0 ; nAVCurComp < nAdjVoxCompNum ; ++ nAVCurComp) {
|
|
|
|
// Numero di triangoli per le componenti connesse
|
|
size_t nCVTriaNum = triStrCurr.vCompoTria[nCVCurComp].size() ;
|
|
size_t nAVTriaNum = triStrAdj.vCompoTria[nAVCurComp].size() ;
|
|
|
|
for ( size_t nCVT = 0 ; nCVT < nCVTriaNum ; ++ nCVT) {
|
|
|
|
bool bModified = false ;
|
|
|
|
for ( size_t nAVT = 0 ; nAVT < nAVTriaNum ; ++ nAVT) {
|
|
|
|
std::vector<int> SharedIndex ;
|
|
|
|
for ( size_t nCurVoxVert = 0 ; nCurVoxVert < 3 ; ++ nCurVoxVert) {
|
|
for ( size_t nAdjVoxVert = 0 ; nAdjVoxVert < 3 ; ++ nAdjVoxVert) {
|
|
|
|
Point3d ptC = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetP( nCurVoxVert) ;
|
|
Point3d ptA = triStrAdj.vCompoTria[nAVCurComp][nAVT].GetP( nAdjVoxVert) ;
|
|
|
|
if ( SqDist( ptC, ptA) < EPS_SMALL * EPS_SMALL) {
|
|
|
|
SharedIndex.emplace_back( int( nCurVoxVert)) ;
|
|
SharedIndex.emplace_back( int( nAdjVoxVert)) ;
|
|
}
|
|
|
|
if ( SharedIndex.size() > 2)
|
|
break ;
|
|
}
|
|
|
|
if ( SharedIndex.size() > 2)
|
|
break ;
|
|
}
|
|
// Si deve operare la modifica dei triangoli
|
|
if ( SharedIndex.size() > 2) {
|
|
|
|
// Controllo sulle normali
|
|
Point3d ptCurFeature = triStrCurr.ptCompoVert[nCVCurComp] ;
|
|
Point3d ptAdjFeature = triStrAdj.ptCompoVert[nAVCurComp] ;
|
|
|
|
Vector3d vtFeature = ptAdjFeature - ptCurFeature ;
|
|
|
|
vtFeature.Normalize() ;
|
|
|
|
double dDot = vtFeature * triStrAdj.vCompoTria[nAVCurComp][nAVT].GetN() ;
|
|
|
|
// Ulteriore controllo sulle normali
|
|
Point3d ptP0 = ptCurFeature ;
|
|
Point3d ptP1 = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetP( SharedIndex[0]) ;
|
|
Point3d ptP2 = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetP( SharedIndex[2]) ;
|
|
Point3d ptP3 = ptAdjFeature ;
|
|
|
|
double dPlane = ( ( ptP1 - ptP0) ^ ( ptP2 - ptP0)) * ( ptP3 - ptP0) ;
|
|
|
|
Vector3d vtN = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetN() ;
|
|
Vector3d vtM = triStrAdj.vCompoTria[nAVCurComp][nAVT].GetN() ;
|
|
|
|
double dDotNM = vtN * vtM ;
|
|
|
|
int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ;
|
|
|
|
if ( nProd != 1 && nProd != - 2 && nProd != 4 &&
|
|
( dDot < - EPS_SMALL || ( dPlane < EPS_SMALL && dDotNM < - 0.95))) {
|
|
|
|
triStrCurr.vCompoTria[nCVCurComp][nCVT].SetP( SharedIndex[0],
|
|
triStrAdj.ptCompoVert[nAVCurComp]) ;
|
|
triStrAdj.vCompoTria[nAVCurComp][nAVT].SetP( SharedIndex[3],
|
|
triStrCurr.ptCompoVert[nCVCurComp]) ;
|
|
|
|
triStrCurr.vCompoTria[nCVCurComp][nCVT].Validate( true) ;
|
|
triStrAdj.vCompoTria[nAVCurComp][nAVT].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].dZVal - dEps &&
|
|
dZ[nGrid] < m_Values[nGrid][nPos][nIndex + 1].dZVal + dEps) {
|
|
|
|
++ nCount ;
|
|
break ;
|
|
}
|
|
nIndex += 2 ;
|
|
}
|
|
}
|
|
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 += 2) {
|
|
|
|
double dx1 = m_Values[1][nDexel][i].dZVal ;
|
|
double dx2 = m_Values[1][nDexel][i+1].dZVal ;
|
|
|
|
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 += 2) {
|
|
|
|
double dy1 = m_Values[2][nDexel][j].dZVal ;
|
|
double dy2 = m_Values[2][nDexel][j+1].dZVal ;
|
|
|
|
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 += 2) {
|
|
|
|
double dz1 = m_Values[0][nDexel][k].dZVal ;
|
|
double dz2 = m_Values[0][nDexel][k+1].dZVal ;
|
|
|
|
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[], Point3d& ptInt, Vector3d& vtNormal) const
|
|
{
|
|
double Eps = EPS_SMALL ;
|
|
|
|
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 ;
|
|
|
|
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 += 2) {
|
|
|
|
double dx1 = m_Values[1][nDexel][i].dZVal ;
|
|
double dx2 = m_Values[1][nDexel][i+1].dZVal ;
|
|
|
|
if ( dx1 < dMinX - Eps && dx2 > dMinX - Eps && dx2 < dMaxX + Eps) {
|
|
ptInt.x = dx2 ;
|
|
vtNormal = m_Values[1][nDexel][i+1].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
else if ( dx1 > dMinX - Eps && dx1 < dMaxX + Eps && dx2 > dMaxX + Eps) {
|
|
ptInt.x = dx1 ;
|
|
vtNormal = m_Values[1][nDexel][i].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}/*
|
|
if ( dx1 < dMinX + Eps && dx2 > dMinX && dx2 < dMaxX - Eps) {
|
|
ptInt.x = dx2 ;
|
|
vtNormal = m_Values[1][nDexel][i+1].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
else if ( dx1 > dMinX + Eps && dx1 < dMaxX + Eps && dx2 > dMaxX + Eps) {
|
|
ptInt.x = dx1 ;
|
|
vtNormal = m_Values[1][nDexel][i].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}*/
|
|
}
|
|
if ( ! bFound) {
|
|
|
|
ptInt.x = ( dMinX + dMaxX) / 2 ;
|
|
/*
|
|
if ( IsThereMat( nVec1[0], nVec1[1], nVec1[2])) {
|
|
|
|
double dY = ( nVec2[1] + 0.5) * m_dStep ;
|
|
double dZ = ( nVec2[2] + 0.5) * m_dStep ;
|
|
|
|
unsigned int YDirMinIndex = 0 ;
|
|
unsigned int ZDirMinIndex = 0 ;
|
|
|
|
double dYMinDist = DBL_MAX ;
|
|
double dZMinDist = DBL_MAX ;
|
|
|
|
unsigned int nYDirGrid = nVec2[0] * m_nNx[2] + nVec2[2] ;
|
|
unsigned int nZDirGrid = nVec2[1] * m_nNx[0] + nVec2[0] ;
|
|
|
|
for ( unsigned int t = 0 ; t < m_Values[0][nZDirGrid].size() ; ++ t) {
|
|
|
|
double dMZ1 = m_Values[0][nZDirGrid][t].dZVal ;
|
|
|
|
if ( abs( dMZ1 - dZ) < dZMinDist) {
|
|
|
|
ZDirMinIndex = t ;
|
|
dZMinDist = abs( dMZ1 - dZ) ;
|
|
}
|
|
}
|
|
|
|
for ( unsigned int t = 0 ; t < m_Values[2][nYDirGrid].size() ; ++ t) {
|
|
|
|
double dMY1 = m_Values[2][nYDirGrid][t].dZVal ;
|
|
|
|
if ( abs( dMY1 - dY) < dYMinDist) {
|
|
|
|
YDirMinIndex = t ;
|
|
dYMinDist = abs( dMY1 - dY) ;
|
|
}
|
|
}
|
|
|
|
unsigned int AbsoluteMinIndex ;
|
|
|
|
if ( dZMinDist < dYMinDist)
|
|
|
|
vtNormal = m_Values[0][nZDirGrid][ZDirMinIndex].vtN ;
|
|
else
|
|
vtNormal = m_Values[0][nYDirGrid][YDirMinIndex].vtN ;
|
|
}
|
|
else {
|
|
|
|
}*/
|
|
}
|
|
}
|
|
|
|
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 ;
|
|
|
|
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 += 2) {
|
|
|
|
double dy1 = m_Values[2][nDexel][j].dZVal ;
|
|
double dy2 = m_Values[2][nDexel][j+1].dZVal ;
|
|
|
|
if ( dy1 < dMinY - Eps && dy2 > dMinY - Eps && dy2 < dMaxY + Eps) {
|
|
ptInt.y = dy2 ;
|
|
vtNormal = m_Values[2][nDexel][j+1].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
else if ( dy1 > dMinY - Eps && dy1 < dMaxY + Eps && dy2 > dMaxY + Eps) {
|
|
ptInt.y = dy1 ;
|
|
vtNormal = m_Values[2][nDexel][j].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}/*
|
|
if ( dy1 < dMinY + Eps && dy2 > dMinY && dy2 < dMaxY - Eps) {
|
|
ptInt.y = dy2 ;
|
|
vtNormal = m_Values[2][nDexel][j+1].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
else if ( dy1 > dMinY + Eps && dy1 < dMaxY + Eps && dy2 > dMaxY + Eps) {
|
|
ptInt.y = dy1 ;
|
|
vtNormal = m_Values[2][nDexel][j].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}*/
|
|
}
|
|
if ( ! bFound) {
|
|
ptInt.y = ( dMinY + dMaxY) / 2 ;
|
|
// Versore Normale ???
|
|
}
|
|
}
|
|
|
|
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 ;
|
|
|
|
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 += 2) {
|
|
|
|
double dz1 = m_Values[0][nDexel][k].dZVal ;
|
|
double dz2 = m_Values[0][nDexel][k+1].dZVal ;
|
|
|
|
if ( dz1 < dMinZ - Eps && dz2 > dMinZ - Eps && dz2 < dMaxZ + Eps) {
|
|
ptInt.z = dz2 ;
|
|
vtNormal = m_Values[0][nDexel][k+1].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
else if ( dz1 > dMinZ - Eps && dz1 < dMaxZ + Eps && dz2 > dMaxZ + Eps) {
|
|
ptInt.z = dz1 ;
|
|
vtNormal = m_Values[0][nDexel][k].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}/*
|
|
if ( dz1 < dMinZ + Eps && dz2 > dMinZ && dz2 < dMaxZ - Eps) {
|
|
ptInt.z = dz2 ;
|
|
vtNormal = m_Values[0][nDexel][k+1].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}
|
|
else if ( dz1 > dMinZ + Eps && dz1 < dMaxZ + Eps && dz2 > dMaxZ + Eps) {
|
|
ptInt.z = dz1 ;
|
|
vtNormal = m_Values[0][nDexel][k].vtN ;
|
|
bFound = true ;
|
|
break ;
|
|
}*/
|
|
}
|
|
if ( ! bFound) {
|
|
ptInt.z = ( dMinZ + dMaxZ) / 2 ;
|
|
// Versore Normale ???
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, Vector3d& vtNormal) const
|
|
{
|
|
double dEps = 2 * EPS_SMALL ;
|
|
|
|
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) {
|
|
|
|
size_t n = nSize - 1 ;
|
|
double dX = m_Values[1][nDexel][n].dZVal ;
|
|
|
|
while ( n > 0 && dX > dMinX - dEps) {
|
|
|
|
if ( dX < dMaxX + dEps) {
|
|
|
|
ptInt.x = dX ;
|
|
vtNormal = m_Values[1][nDexel][n].vtN ;
|
|
break ;
|
|
}
|
|
|
|
if ( n == 1)
|
|
break ;
|
|
|
|
n -= 2 ;
|
|
dX = m_Values[1][nDexel][n].dZVal ;
|
|
}
|
|
}
|
|
else {
|
|
|
|
size_t n = 0 ;
|
|
double dX = m_Values[1][nDexel][0].dZVal ;
|
|
|
|
while ( n <= nSize - 2 && dX < dMaxX + dEps) {
|
|
|
|
if ( dX > dMinX - dEps) {
|
|
|
|
ptInt.x = dX ;
|
|
vtNormal = m_Values[1][nDexel][n].vtN ;
|
|
break ;
|
|
}
|
|
|
|
if ( n == nSize - 2)
|
|
break ;
|
|
|
|
n += 2 ;
|
|
dX = m_Values[1][nDexel][n].dZVal ;
|
|
}
|
|
}
|
|
}
|
|
|
|
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) {
|
|
|
|
size_t n = nSize - 1 ;
|
|
double dY = m_Values[2][nDexel][n].dZVal ;
|
|
|
|
while ( n > 0 && dY > dMinY - dEps) {
|
|
|
|
if ( dY < dMaxY + dEps) {
|
|
|
|
ptInt.y = dY ;
|
|
vtNormal = m_Values[2][nDexel][n].vtN ;
|
|
break ;
|
|
}
|
|
|
|
if ( n == 1)
|
|
break ;
|
|
|
|
n -= 2 ;
|
|
dY = m_Values[2][nDexel][n].dZVal ;
|
|
}
|
|
}
|
|
else {
|
|
|
|
size_t n = 0 ;
|
|
double dY = m_Values[2][nDexel][0].dZVal ;
|
|
|
|
while ( n <= nSize - 2 && dY < dMaxY + dEps) {
|
|
|
|
if ( dY > dMinY - dEps) {
|
|
|
|
ptInt.y = dY ;
|
|
vtNormal = m_Values[2][nDexel][n].vtN ;
|
|
break ;
|
|
}
|
|
|
|
if ( n == nSize - 2)
|
|
break ;
|
|
|
|
n += 2 ;
|
|
dY = m_Values[2][nDexel][n].dZVal ;
|
|
}
|
|
}
|
|
}
|
|
|
|
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) {
|
|
|
|
size_t n = nSize - 1 ;
|
|
double dZ = m_Values[0][nDexel][n].dZVal ;
|
|
|
|
while ( n > 0 && dZ > dMinZ - dEps) {
|
|
|
|
if ( dZ < dMaxZ + dEps) {
|
|
|
|
ptInt.z = dZ ;
|
|
vtNormal = m_Values[0][nDexel][n].vtN ;
|
|
break ;
|
|
}
|
|
|
|
if ( n == 1)
|
|
break ;
|
|
|
|
n -= 2 ;
|
|
dZ = m_Values[0][nDexel][n].dZVal ;
|
|
}
|
|
}
|
|
else {
|
|
|
|
size_t n = 0 ;
|
|
double dZ = m_Values[0][nDexel][0].dZVal ;
|
|
|
|
while ( n <= nSize - 2 && dZ < dMaxZ + dEps) {
|
|
|
|
if ( dZ > dMinZ - dEps) {
|
|
|
|
ptInt.z = dZ ;
|
|
vtNormal = m_Values[0][nDexel][n].vtN ;
|
|
break ;
|
|
}
|
|
|
|
if ( n == nSize - 2)
|
|
break ;
|
|
|
|
n += 2 ;
|
|
dZ = m_Values[0][nDexel][n].dZVal ;
|
|
}
|
|
}
|
|
}
|
|
|
|
return true ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
bool
|
|
VolZmap::GetBlockIJK( int nIJK[], int nBlock) 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::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 ;
|
|
}
|