manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Include header cleaning and reordering.

This commit is contained in:
Paolo Cignoni 2013-11-25 10:32:41 +00:00
parent f2bbdb787a
commit cc72b3e3e1
4 changed files with 1399 additions and 1413 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

715
vcg/complex/algorithms/create/extended_marching_cubes.h
View File

@ -27,342 +27,337 @@
#define __VCG_EXTENDED_MARCHING_CUBES #define __VCG_EXTENDED_MARCHING_CUBES
#include <float.h> #include <float.h>
#include <assert.h>
#include <vector>
#include <vcg/math/base.h>
#include <vcg/simplex/face/topology.h> #include <vcg/simplex/face/topology.h>
#include <vcg/complex/algorithms/update/normal.h> #include <vcg/complex/algorithms/update/normal.h>
#include <vcg/complex/algorithms/update/topology.h> #include <vcg/complex/algorithms/update/topology.h>
#include <vcg/complex/allocate.h>
#include <vcg/complex/algorithms/clean.h> #include <vcg/complex/algorithms/clean.h>
#include <vcg/space/point3.h>
#include "emc_lookup_table.h" #include "emc_lookup_table.h"
#include <eigenlib/Eigen/SVD> #include <eigenlib/Eigen/SVD>
namespace vcg namespace vcg
{ {
namespace tri namespace tri
{ {
// Doxygen documentation // Doxygen documentation
/** \addtogroup trimesh */ /** \addtogroup trimesh */
/*@{*/ /*@{*/
/* /*
* Cube description: * Cube description:
* 3 ________ 2 _____2__ * 3 ________ 2 _____2__
* /| /| / | /| * /| /| / | /|
* / | / | 11/ 3 10/ | * / | / | 11/ 3 10/ |
* 7 /_______ / | /__6_|__ / |1 * 7 /_______ / | /__6_|__ / |1
* | | |6 | | | | * | | |6 | | | |
* | 0|__|_____|1 | |__|_0|__| * | 0|__|_____|1 | |__|_0|__|
* | / | / 7 8/ 5 / * | / | / 7 8/ 5 /
* | / | / | / | /9 * | / | / | / | /9
* |/_______|/ |/___4___|/ * |/_______|/ |/___4___|/
* 4 5 * 4 5
*/ */
//! This class implements the Extended Marching Cubes algorithm. //! This class implements the Extended Marching Cubes algorithm.
/*! /*!
* The implementation is enough generic: this class works only on one volume cell for each * The implementation is enough generic: this class works only on one volume cell for each
* call to <CODE>ProcessCell</CODE>. Using the field value at the cell corners, it adds to the * call to <CODE>ProcessCell</CODE>. Using the field value at the cell corners, it adds to the
* mesh the triangles set approximating the surface that cross that cell. * mesh the triangles set approximating the surface that cross that cell.
* @param TRIMESH_TYPE (Template parameter) the mesh type that will be constructed * @param TRIMESH_TYPE (Template parameter) the mesh type that will be constructed
* @param WALKER_TYPE (Template parameter) the class that implements the traversal ordering of the volume. * @param WALKER_TYPE (Template parameter) the class that implements the traversal ordering of the volume.
**/ **/
template<class TRIMESH_TYPE, class WALKER_TYPE> template<class TRIMESH_TYPE, class WALKER_TYPE>
class ExtendedMarchingCubes class ExtendedMarchingCubes
{ {
public: public:
#if defined(__GNUC__) #if defined(__GNUC__)
typedef unsigned int size_t; typedef unsigned int size_t;
#else #else
#ifdef _WIN64 #ifdef _WIN64
typedef unsigned __int64 size_t; typedef unsigned __int64 size_t;
#else #else
typedef _W64 unsigned int size_t; typedef _W64 unsigned int size_t;
#endif #endif
#endif #endif
typedef typename vcg::tri::Allocator< TRIMESH_TYPE > AllocatorType; typedef typename vcg::tri::Allocator< TRIMESH_TYPE > AllocatorType;
typedef typename TRIMESH_TYPE::ScalarType ScalarType; typedef typename TRIMESH_TYPE::ScalarType ScalarType;
typedef typename TRIMESH_TYPE::VertexType VertexType; typedef typename TRIMESH_TYPE::VertexType VertexType;
typedef typename TRIMESH_TYPE::VertexPointer VertexPointer; typedef typename TRIMESH_TYPE::VertexPointer VertexPointer;
typedef typename TRIMESH_TYPE::VertexIterator VertexIterator; typedef typename TRIMESH_TYPE::VertexIterator VertexIterator;
typedef typename TRIMESH_TYPE::FaceType FaceType; typedef typename TRIMESH_TYPE::FaceType FaceType;
typedef typename TRIMESH_TYPE::FacePointer FacePointer; typedef typename TRIMESH_TYPE::FacePointer FacePointer;
typedef typename TRIMESH_TYPE::FaceIterator FaceIterator; typedef typename TRIMESH_TYPE::FaceIterator FaceIterator;
typedef typename TRIMESH_TYPE::CoordType CoordType; typedef typename TRIMESH_TYPE::CoordType CoordType;
typedef typename TRIMESH_TYPE::CoordType* CoordPointer; typedef typename TRIMESH_TYPE::CoordType* CoordPointer;
struct LightEdge struct LightEdge
{ {
LightEdge(size_t _face, size_t _edge):face(_face), edge(_edge) { } LightEdge(size_t _face, size_t _edge):face(_face), edge(_edge) { }
size_t face, edge; size_t face, edge;
}; };
/*! /*!
* Constructor * Constructor
* \param mesh The mesh that will be constructed * \param mesh The mesh that will be constructed
* \param volume The volume describing the field * \param volume The volume describing the field
* \param walker The class implementing the traversal policy * \param walker The class implementing the traversal policy
* \param angle The feature detection threshold misuring the sharpness of a feature(default is 30 degree) * \param angle The feature detection threshold misuring the sharpness of a feature(default is 30 degree)
*/ */
ExtendedMarchingCubes(TRIMESH_TYPE &mesh, WALKER_TYPE &walker, ScalarType angle=30) ExtendedMarchingCubes(TRIMESH_TYPE &mesh, WALKER_TYPE &walker, ScalarType angle=30)
{ {
_mesh = &mesh; _mesh = &mesh;
_walker = &walker; _walker = &walker;
_featureAngle = vcg::math::ToRad(angle); _featureAngle = vcg::math::ToRad(angle);
_initialized = _finalized = false; _initialized = _finalized = false;
}; };
/*! /*!
* Execute the initialiazation. * Execute the initialiazation.
* This method must be executed before the first call to <CODE>ApplyEMC</CODE> * This method must be executed before the first call to <CODE>ApplyEMC</CODE>
*/ */
void Initialize() void Initialize()
{ {
assert(!_initialized && !_finalized); assert(!_initialized && !_finalized);
_featureFlag = VertexType::NewBitFlag(); _featureFlag = VertexType::NewBitFlag();
_initialized = true; _initialized = true;
}; };
/*! /*!
* *
* This method must be executed after the last call to <CODE>ApplyEMC</CODE> * This method must be executed after the last call to <CODE>ApplyEMC</CODE>
*/ */
void Finalize() void Finalize()
{ {
assert(_initialized && !_finalized); assert(_initialized && !_finalized);
FlipEdges(); FlipEdges();
VertexIterator v_iter = _mesh->vert.begin(); VertexIterator v_iter = _mesh->vert.begin();
VertexIterator v_end = _mesh->vert.end(); VertexIterator v_end = _mesh->vert.end();
for ( ; v_iter!=v_end; v_iter++) for ( ; v_iter!=v_end; v_iter++)
v_iter->ClearUserBit( _featureFlag ); v_iter->ClearUserBit( _featureFlag );
VertexType::DeleteBitFlag( _featureFlag ); VertexType::DeleteBitFlag( _featureFlag );
_featureFlag = 0; _featureFlag = 0;
_mesh = NULL; _mesh = NULL;
_walker = NULL; _walker = NULL;
_finalized = true; _finalized = true;
}; };
/*! /*!
* Apply the <I>extended marching cubes</I> algorithm to the volume cell identified by the two points <CODE>min</CODE> and <CODE>max</CODE>. * Apply the <I>extended marching cubes</I> algorithm to the volume cell identified by the two points <CODE>min</CODE> and <CODE>max</CODE>.
* All the three coordinates of the first point must be smaller than the respectives three coordinatas of the second point. * All the three coordinates of the first point must be smaller than the respectives three coordinatas of the second point.
* \param min the first point * \param min the first point
* \param max the second point * \param max the second point
*/ */
void ProcessCell(const vcg::Point3i &min, const vcg::Point3i &max) void ProcessCell(const vcg::Point3i &min, const vcg::Point3i &max)
{ {
assert(_initialized && !_finalized); assert(_initialized && !_finalized);
assert(min[0]<max[0] && min[1]<max[1] && min[2]<max[2]); assert(min[0]<max[0] && min[1]<max[1] && min[2]<max[2]);
_corners[0].X()=min.X(); _corners[0].Y()=min.Y(); _corners[0].Z()=min.Z(); _corners[0].X()=min.X(); _corners[0].Y()=min.Y(); _corners[0].Z()=min.Z();
_corners[1].X()=max.X(); _corners[1].Y()=min.Y(); _corners[1].Z()=min.Z(); _corners[1].X()=max.X(); _corners[1].Y()=min.Y(); _corners[1].Z()=min.Z();
_corners[2].X()=max.X(); _corners[2].Y()=max.Y(); _corners[2].Z()=min.Z(); _corners[2].X()=max.X(); _corners[2].Y()=max.Y(); _corners[2].Z()=min.Z();
_corners[3].X()=min.X(); _corners[3].Y()=max.Y(); _corners[3].Z()=min.Z(); _corners[3].X()=min.X(); _corners[3].Y()=max.Y(); _corners[3].Z()=min.Z();
_corners[4].X()=min.X(); _corners[4].Y()=min.Y(); _corners[4].Z()=max.Z(); _corners[4].X()=min.X(); _corners[4].Y()=min.Y(); _corners[4].Z()=max.Z();
_corners[5].X()=max.X(); _corners[5].Y()=min.Y(); _corners[5].Z()=max.Z(); _corners[5].X()=max.X(); _corners[5].Y()=min.Y(); _corners[5].Z()=max.Z();
_corners[6].X()=max.X(); _corners[6].Y()=max.Y(); _corners[6].Z()=max.Z(); _corners[6].X()=max.X(); _corners[6].Y()=max.Y(); _corners[6].Z()=max.Z();
_corners[7].X()=min.X(); _corners[7].Y()=max.Y(); _corners[7].Z()=max.Z(); _corners[7].X()=min.X(); _corners[7].Y()=max.Y(); _corners[7].Z()=max.Z();
unsigned char cubetype = 0; unsigned char cubetype = 0;
if ((_field[0] = _walker->V(_corners[0].X(), _corners[0].Y(), _corners[0].Z())) >= 0) cubetype+= 1; if ((_field[0] = _walker->V(_corners[0].X(), _corners[0].Y(), _corners[0].Z())) >= 0) cubetype+= 1;
if ((_field[1] = _walker->V(_corners[1].X(), _corners[1].Y(), _corners[1].Z())) >= 0) cubetype+= 2; if ((_field[1] = _walker->V(_corners[1].X(), _corners[1].Y(), _corners[1].Z())) >= 0) cubetype+= 2;
if ((_field[2] = _walker->V(_corners[2].X(), _corners[2].Y(), _corners[2].Z())) >= 0) cubetype+= 4; if ((_field[2] = _walker->V(_corners[2].X(), _corners[2].Y(), _corners[2].Z())) >= 0) cubetype+= 4;
if ((_field[3] = _walker->V(_corners[3].X(), _corners[3].Y(), _corners[3].Z())) >= 0) cubetype+= 8; if ((_field[3] = _walker->V(_corners[3].X(), _corners[3].Y(), _corners[3].Z())) >= 0) cubetype+= 8;
if ((_field[4] = _walker->V(_corners[4].X(), _corners[4].Y(), _corners[4].Z())) >= 0) cubetype+= 16; if ((_field[4] = _walker->V(_corners[4].X(), _corners[4].Y(), _corners[4].Z())) >= 0) cubetype+= 16;
if ((_field[5] = _walker->V(_corners[5].X(), _corners[5].Y(), _corners[5].Z())) >= 0) cubetype+= 32; if ((_field[5] = _walker->V(_corners[5].X(), _corners[5].Y(), _corners[5].Z())) >= 0) cubetype+= 32;
if ((_field[6] = _walker->V(_corners[6].X(), _corners[6].Y(), _corners[6].Z())) >= 0) cubetype+= 64; if ((_field[6] = _walker->V(_corners[6].X(), _corners[6].Y(), _corners[6].Z())) >= 0) cubetype+= 64;
if ((_field[7] = _walker->V(_corners[7].X(), _corners[7].Y(), _corners[7].Z())) >= 0) cubetype+=128; if ((_field[7] = _walker->V(_corners[7].X(), _corners[7].Y(), _corners[7].Z())) >= 0) cubetype+=128;
if (cubetype==0 || cubetype==255) if (cubetype==0 || cubetype==255)
return; return;
size_t vertices_idx[12]; size_t vertices_idx[12];
memset(vertices_idx, -1, 12*sizeof(size_t)); memset(vertices_idx, -1, 12*sizeof(size_t));
int code = EMCLookUpTable::EdgeTable(cubetype); int code = EMCLookUpTable::EdgeTable(cubetype);
VertexPointer vp = NULL; VertexPointer vp = NULL;
if ( 1&code ) { _walker->GetXIntercept(_corners[0], _corners[1], vp); vertices_idx[ 0] = vp - &_mesh->vert[0]; } if ( 1&code ) { _walker->GetXIntercept(_corners[0], _corners[1], vp); vertices_idx[ 0] = vp - &_mesh->vert[0]; }
if ( 2&code ) { _walker->GetYIntercept(_corners[1], _corners[2], vp); vertices_idx[ 1] = vp - &_mesh->vert[0]; } if ( 2&code ) { _walker->GetYIntercept(_corners[1], _corners[2], vp); vertices_idx[ 1] = vp - &_mesh->vert[0]; }
if ( 4&code ) { _walker->GetXIntercept(_corners[3], _corners[2], vp); vertices_idx[ 2] = vp - &_mesh->vert[0]; } if ( 4&code ) { _walker->GetXIntercept(_corners[3], _corners[2], vp); vertices_idx[ 2] = vp - &_mesh->vert[0]; }
if ( 8&code ) { _walker->GetYIntercept(_corners[0], _corners[3], vp); vertices_idx[ 3] = vp - &_mesh->vert[0]; } if ( 8&code ) { _walker->GetYIntercept(_corners[0], _corners[3], vp); vertices_idx[ 3] = vp - &_mesh->vert[0]; }
if ( 16&code ) { _walker->GetXIntercept(_corners[4], _corners[5], vp); vertices_idx[ 4] = vp - &_mesh->vert[0]; } if ( 16&code ) { _walker->GetXIntercept(_corners[4], _corners[5], vp); vertices_idx[ 4] = vp - &_mesh->vert[0]; }
if ( 32&code ) { _walker->GetYIntercept(_corners[5], _corners[6], vp); vertices_idx[ 5] = vp - &_mesh->vert[0]; } if ( 32&code ) { _walker->GetYIntercept(_corners[5], _corners[6], vp); vertices_idx[ 5] = vp - &_mesh->vert[0]; }
if ( 64&code ) { _walker->GetXIntercept(_corners[7], _corners[6], vp); vertices_idx[ 6] = vp - &_mesh->vert[0]; } if ( 64&code ) { _walker->GetXIntercept(_corners[7], _corners[6], vp); vertices_idx[ 6] = vp - &_mesh->vert[0]; }
if ( 128&code ) { _walker->GetYIntercept(_corners[4], _corners[7], vp); vertices_idx[ 7] = vp - &_mesh->vert[0]; } if ( 128&code ) { _walker->GetYIntercept(_corners[4], _corners[7], vp); vertices_idx[ 7] = vp - &_mesh->vert[0]; }
if ( 256&code ) { _walker->GetZIntercept(_corners[0], _corners[4], vp); vertices_idx[ 8] = vp - &_mesh->vert[0]; } if ( 256&code ) { _walker->GetZIntercept(_corners[0], _corners[4], vp); vertices_idx[ 8] = vp - &_mesh->vert[0]; }
if ( 512&code ) { _walker->GetZIntercept(_corners[1], _corners[5], vp); vertices_idx[ 9] = vp - &_mesh->vert[0]; } if ( 512&code ) { _walker->GetZIntercept(_corners[1], _corners[5], vp); vertices_idx[ 9] = vp - &_mesh->vert[0]; }
if (1024&code ) { _walker->GetZIntercept(_corners[2], _corners[6], vp); vertices_idx[10] = vp - &_mesh->vert[0]; } if (1024&code ) { _walker->GetZIntercept(_corners[2], _corners[6], vp); vertices_idx[10] = vp - &_mesh->vert[0]; }
if (2048&code ) { _walker->GetZIntercept(_corners[3], _corners[7], vp); vertices_idx[11] = vp - &_mesh->vert[0]; } if (2048&code ) { _walker->GetZIntercept(_corners[3], _corners[7], vp); vertices_idx[11] = vp - &_mesh->vert[0]; }
int m, n, vertices_num; int m, n, vertices_num;
int components = EMCLookUpTable::TriTable(cubetype, 1)[0]; //unsigned int components = triTable[cubetype][1][0]; int components = EMCLookUpTable::TriTable(cubetype, 1)[0]; //unsigned int components = triTable[cubetype][1][0];
int *indices = &EMCLookUpTable::TriTable(cubetype, 1)[components+1]; //int *indices = &EMCLookUpTable::TriTable(cubetype, 1, components+1); int *indices = &EMCLookUpTable::TriTable(cubetype, 1)[components+1]; //int *indices = &EMCLookUpTable::TriTable(cubetype, 1, components+1);
std::vector< size_t > vertices_list; std::vector< size_t > vertices_list;
for (m=1; m<=components; m++) for (m=1; m<=components; m++)
{ {
// current sheet contains vertices_num vertices // current sheet contains vertices_num vertices
vertices_num = EMCLookUpTable::TriTable(cubetype, 1)[m]; //vertices_num = triTable[cubetype][1][m]; vertices_num = EMCLookUpTable::TriTable(cubetype, 1)[m]; //vertices_num = triTable[cubetype][1][m];
// collect vertices // collect vertices
vertices_list.clear(); vertices_list.clear();
for (n=0; n<vertices_num; ++n) for (n=0; n<vertices_num; ++n)
vertices_list.push_back( vertices_idx[ indices[n] ] ); vertices_list.push_back( vertices_idx[ indices[n] ] );
VertexPointer feature = FindFeature( vertices_list ); VertexPointer feature = FindFeature( vertices_list );
if (feature != NULL) // i.e. is a valid vertex if (feature != NULL) // i.e. is a valid vertex
{ {
// feature -> create triangle fan around feature vertex // feature -> create triangle fan around feature vertex
size_t feature_idx = feature - &_mesh->vert[0]; size_t feature_idx = feature - &_mesh->vert[0];
size_t face_idx = _mesh->face.size(); size_t face_idx = _mesh->face.size();
vertices_list.push_back( vertices_list[0] ); vertices_list.push_back( vertices_list[0] );
AllocatorType::AddFaces(*_mesh, (int) vertices_num); AllocatorType::AddFaces(*_mesh, (int) vertices_num);
for (int j=0; j<vertices_num; ++j, face_idx++) for (int j=0; j<vertices_num; ++j, face_idx++)
{ {
_mesh->face[face_idx].V(0) = &_mesh->vert[ vertices_list[j ] ]; _mesh->face[face_idx].V(0) = &_mesh->vert[ vertices_list[j ] ];
_mesh->face[face_idx].V(1) = &_mesh->vert[ vertices_list[j+1] ]; _mesh->face[face_idx].V(1) = &_mesh->vert[ vertices_list[j+1] ];
_mesh->face[face_idx].V(2) = &_mesh->vert[ feature_idx ]; _mesh->face[face_idx].V(2) = &_mesh->vert[ feature_idx ];
} }
} }
else else
{ {
// no feature -> old marching cubes triangle table // no feature -> old marching cubes triangle table
for (int j=0; EMCLookUpTable::PolyTable(vertices_num, j) != -1; j+=3) //for (int j=0; polyTable[vertices_num][j] != -1; j+=3) for (int j=0; EMCLookUpTable::PolyTable(vertices_num, j) != -1; j+=3) //for (int j=0; polyTable[vertices_num][j] != -1; j+=3)
{ {
size_t face_idx = _mesh->face.size(); size_t face_idx = _mesh->face.size();
AllocatorType::AddFaces(*_mesh, 1); AllocatorType::AddFaces(*_mesh, 1);
//_mesh->face[ face_idx].V(0) = &_mesh->vert[ vertices_idx[ indices[ polyTable[vertices_num][j ] ] ] ]; //_mesh->face[ face_idx].V(0) = &_mesh->vert[ vertices_idx[ indices[ polyTable[vertices_num][j ] ] ] ];
//_mesh->face[ face_idx].V(1) = &_mesh->vert[ vertices_idx[ indices[ polyTable[vertices_num][j+1] ] ] ]; //_mesh->face[ face_idx].V(1) = &_mesh->vert[ vertices_idx[ indices[ polyTable[vertices_num][j+1] ] ] ];
//_mesh->face[ face_idx].V(2) = &_mesh->vert[ vertices_idx[ indices[ polyTable[vertices_num][j+2] ] ] ]; //_mesh->face[ face_idx].V(2) = &_mesh->vert[ vertices_idx[ indices[ polyTable[vertices_num][j+2] ] ] ];
_mesh->face[ face_idx].V(0) = &_mesh->vert[ vertices_idx[ indices[ EMCLookUpTable::PolyTable(vertices_num, j ) ] ] ]; _mesh->face[ face_idx].V(0) = &_mesh->vert[ vertices_idx[ indices[ EMCLookUpTable::PolyTable(vertices_num, j ) ] ] ];
_mesh->face[ face_idx].V(1) = &_mesh->vert[ vertices_idx[ indices[ EMCLookUpTable::PolyTable(vertices_num, j+1) ] ] ]; _mesh->face[ face_idx].V(1) = &_mesh->vert[ vertices_idx[ indices[ EMCLookUpTable::PolyTable(vertices_num, j+1) ] ] ];
_mesh->face[ face_idx].V(2) = &_mesh->vert[ vertices_idx[ indices[ EMCLookUpTable::PolyTable(vertices_num, j+2) ] ] ]; _mesh->face[ face_idx].V(2) = &_mesh->vert[ vertices_idx[ indices[ EMCLookUpTable::PolyTable(vertices_num, j+2) ] ] ];
} }
} }
indices += vertices_num; indices += vertices_num;
} }
}; // end of ApplyEMC }; // end of ApplyEMC
private: private:
/*! /*!
*/ */
WALKER_TYPE *_walker; WALKER_TYPE *_walker;
/*! /*!
*/ */
TRIMESH_TYPE *_mesh; TRIMESH_TYPE *_mesh;
/*! /*!
*/ */
bool _initialized;; bool _initialized;;
/*! /*!
*/ */
bool _finalized; bool _finalized;
/*! /*!
* The feature detection threshold misuring the sharpness of a feature * The feature detection threshold misuring the sharpness of a feature
*/ */
ScalarType _featureAngle; ScalarType _featureAngle;
/*! /*!
* The flag used for marking the feature vertices. * The flag used for marking the feature vertices.
*/ */
int _featureFlag; int _featureFlag;
/*! /*!
* Array of the 8 corners of the volume cell being processed * Array of the 8 corners of the volume cell being processed
*/ */
vcg::Point3i _corners[8]; vcg::Point3i _corners[8];
/*! /*!
* The field value at the cell corners * The field value at the cell corners
*/ */
ScalarType _field[8]; ScalarType _field[8];
/*! /*!
* Tests if the surface patch crossing the current cell contains a sharp feature * Tests if the surface patch crossing the current cell contains a sharp feature
* \param vertices_idx The list of vertex indices intersecting the edges of the current cell * \param vertices_idx The list of vertex indices intersecting the edges of the current cell
* \return The pointer to the new Vertex if a feature is detected; NULL otherwise. * \return The pointer to the new Vertex if a feature is detected; NULL otherwise.
*/ */
VertexPointer FindFeature(const std::vector<size_t> &vertices_idx) VertexPointer FindFeature(const std::vector<size_t> &vertices_idx)
{ {
unsigned int i, j, rank; unsigned int i, j, rank;
size_t vertices_num = (size_t) vertices_idx.size(); size_t vertices_num = (size_t) vertices_idx.size();
CoordType *points = new CoordType[ vertices_num ]; CoordType *points = new CoordType[ vertices_num ];
CoordType *normals = new CoordType[ vertices_num ]; CoordType *normals = new CoordType[ vertices_num ];
Box3<ScalarType> bb; Box3<ScalarType> bb;
for (i=0; i<vertices_num; i++) for (i=0; i<vertices_num; i++)
{ {
points[i] = _mesh->vert[ vertices_idx[i] ].P(); points[i] = _mesh->vert[ vertices_idx[i] ].P();
normals[i].Import(_mesh->vert[ vertices_idx[i] ].N()); normals[i].Import(_mesh->vert[ vertices_idx[i] ].N());
bb.Add(points[i]); bb.Add(points[i]);
} }
// move barycenter of points into (0, 0, 0) // move barycenter of points into (0, 0, 0)
CoordType center((ScalarType) 0.0, (ScalarType) 0.0, (ScalarType) 0.0); CoordType center((ScalarType) 0.0, (ScalarType) 0.0, (ScalarType) 0.0);
for (i=0; i<vertices_num; ++i) for (i=0; i<vertices_num; ++i)
center += points[i]; center += points[i];
center /= (ScalarType) vertices_num; center /= (ScalarType) vertices_num;
for (i=0; i<vertices_num; ++i) for (i=0; i<vertices_num; ++i)
points[i] -= center; points[i] -= center;
// normal angle criterion // normal angle criterion
double c, minC, maxC; double c, minC, maxC;
CoordType axis; CoordType axis;
for (minC=1.0, i=0; i<vertices_num-1; ++i) for (minC=1.0, i=0; i<vertices_num-1; ++i)
{ {
for (j=i+1; j<vertices_num; ++j) for (j=i+1; j<vertices_num; ++j)
{ {
c = normals[i]*normals[j]; c = normals[i]*normals[j];
if (c < minC) if (c < minC)
{ {
minC = c; minC = c;
axis = normals[i] ^ normals[j]; axis = normals[i] ^ normals[j];
} }
} }
} //end for (minC=1.0, i=0; i<vertNumber; ++i) } //end for (minC=1.0, i=0; i<vertNumber; ++i)
if (minC > cos(_featureAngle)) if (minC > cos(_featureAngle))
return NULL; // invalid vertex return NULL; // invalid vertex
// ok, we have a feature: is it edge or corner, i.e. rank 2 or 3 ? // ok, we have a feature: is it edge or corner, i.e. rank 2 or 3 ?
axis.Normalize(); axis.Normalize();
for (minC=1.0, maxC=-1.0, i=0; i<vertices_num; ++i) for (minC=1.0, maxC=-1.0, i=0; i<vertices_num; ++i)
{ {
c = axis * normals[i]; c = axis * normals[i];
if (c < minC) minC = c; if (c < minC) minC = c;
if (c > maxC) maxC = c; if (c > maxC) maxC = c;
} }
c = std::max< double >(fabs(minC), fabs(maxC)); c = std::max< double >(fabs(minC), fabs(maxC));
c = sqrt(1.0-c*c); c = sqrt(1.0-c*c);
rank = (c > cos(_featureAngle) ? 2 : 3); rank = (c > cos(_featureAngle) ? 2 : 3);
// setup linear system (find intersection of tangent planes) // setup linear system (find intersection of tangent planes)
//--vcg::ndim::Matrix<double> A((unsigned int) vertices_num, 3); //--vcg::ndim::Matrix<double> A((unsigned int) vertices_num, 3);
Eigen::MatrixXd A(vertices_num,3); Eigen::MatrixXd A(vertices_num,3);
//--double *b = new double[ vertices_num ]; //--double *b = new double[ vertices_num ];
Eigen::MatrixXd b(vertices_num,1); Eigen::MatrixXd b(vertices_num,1);
for (i=0; i<vertices_num; ++i) for (i=0; i<vertices_num; ++i)
{ {
//--A[i][0] = normals[i][0]; //--A[i][0] = normals[i][0];
//--A[i][1] = normals[i][1]; //--A[i][1] = normals[i][1];
//--A[i][2] = normals[i][2]; //--A[i][2] = normals[i][2];
//--b[i] = (points[i] * normals[i]); //--b[i] = (points[i] * normals[i]);
A(i,0) = normals[i][0]; A(i,0) = normals[i][0];
A(i,1) = normals[i][1]; A(i,1) = normals[i][1];
A(i,2) = normals[i][2]; A(i,2) = normals[i][2];
b(i) = (points[i] * normals[i]); b(i) = (points[i] * normals[i]);
} }
// SVD of matrix A // SVD of matrix A
Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV); Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
Eigen::MatrixXd sol(3,1); Eigen::MatrixXd sol(3,1);
sol=svd.solve(b); sol=svd.solve(b);
// vcg::ndim::Matrix<double> V(3, 3); // vcg::ndim::Matrix<double> V(3, 3);
// double *w = new double[vertices_num]; // double *w = new double[vertices_num];
// vcg::SingularValueDecomposition< typename vcg::ndim::Matrix<double> > (A, w, V, LeaveUnsorted, 100); // vcg::SingularValueDecomposition< typename vcg::ndim::Matrix<double> > (A, w, V, LeaveUnsorted, 100);
// rank == 2 -> suppress smallest singular value // rank == 2 -> suppress smallest singular value
// if (rank == 2) // if (rank == 2)
// { // {
// double smin = DBL_MAX; // the max value, as defined in <float.h> // double smin = DBL_MAX; // the max value, as defined in <float.h>
@ -384,91 +379,91 @@ namespace vcg
// double *x = new double[3]; // double *x = new double[3];
// vcg::SingularValueBacksubstitution< vcg::ndim::Matrix<double> >(A, w, V, x, b); // vcg::SingularValueBacksubstitution< vcg::ndim::Matrix<double> >(A, w, V, x, b);
// transform x to world coords // transform x to world coords
//--CoordType point((ScalarType) x[0], (ScalarType) x[1], (ScalarType) x[2]); //--CoordType point((ScalarType) x[0], (ScalarType) x[1], (ScalarType) x[2]);
CoordType point((ScalarType) sol(0), (ScalarType) sol(1), (ScalarType) sol(2)); CoordType point((ScalarType) sol(0), (ScalarType) sol(1), (ScalarType) sol(2));
point += center; point += center;
// Safety check if the feature point found by svd is // Safety check if the feature point found by svd is
// out of the bbox of the vertices perhaps it is better to put it back in the center... // out of the bbox of the vertices perhaps it is better to put it back in the center...
if(!bb.IsIn(point)) point = center; if(!bb.IsIn(point)) point = center;
// insert the feature-point // insert the feature-point
VertexPointer mean_point = &*AllocatorType::AddVertices( *_mesh, 1); VertexPointer mean_point = &*AllocatorType::AddVertices( *_mesh, 1);
mean_point->SetUserBit(_featureFlag); mean_point->SetUserBit(_featureFlag);
mean_point->P() = point; mean_point->P() = point;
mean_point->N().SetZero(); mean_point->N().SetZero();
// delete []x; // delete []x;
delete []points; delete []points;
delete []normals; delete []normals;
return mean_point; return mean_point;
} // end of FindFeature } // end of FindFeature
/*! /*!
* Postprocessing step performed during the finalization tha flip some of the mesh edges. * Postprocessing step performed during the finalization tha flip some of the mesh edges.
* The flipping criterion is quite simple: each edge is flipped if it will connect two * The flipping criterion is quite simple: each edge is flipped if it will connect two
* feature samples after the flip. * feature samples after the flip.
*/ */
void FlipEdges() void FlipEdges()
{ {
std::vector< LightEdge > edges; std::vector< LightEdge > edges;
for (FaceIterator fi = _mesh->face.begin(); fi!=_mesh->face.end(); fi++) for (FaceIterator fi = _mesh->face.begin(); fi!=_mesh->face.end(); fi++)
{ {
size_t i = tri::Index(*_mesh,*fi); size_t i = tri::Index(*_mesh,*fi);
if (fi->V(1) > fi->V(0)) edges.push_back( LightEdge(i,0) ); if (fi->V(1) > fi->V(0)) edges.push_back( LightEdge(i,0) );
if (fi->V(2) > fi->V(1)) edges.push_back( LightEdge(i,1) ); if (fi->V(2) > fi->V(1)) edges.push_back( LightEdge(i,1) );
if (fi->V(0) > fi->V(2)) edges.push_back( LightEdge(i,2) ); if (fi->V(0) > fi->V(2)) edges.push_back( LightEdge(i,2) );
} }
vcg::tri::UpdateTopology< TRIMESH_TYPE >::FaceFace( *_mesh ); vcg::tri::UpdateTopology< TRIMESH_TYPE >::FaceFace( *_mesh );
// Select all the triangles that has a vertex shared with a non manifold edge. // Select all the triangles that has a vertex shared with a non manifold edge.
int nonManifEdge = tri::Clean< TRIMESH_TYPE >::CountNonManifoldEdgeFF(*_mesh,true); int nonManifEdge = tri::Clean< TRIMESH_TYPE >::CountNonManifoldEdgeFF(*_mesh,true);
if(nonManifEdge >0) if(nonManifEdge >0)
tri::UpdateSelection< TRIMESH_TYPE >::FaceFromVertexLoose(*_mesh); tri::UpdateSelection< TRIMESH_TYPE >::FaceFromVertexLoose(*_mesh);
//qDebug("Got %i non manif edges",nonManifEdge); //qDebug("Got %i non manif edges",nonManifEdge);
typename std::vector< LightEdge >::iterator e_it = edges.begin(); typename std::vector< LightEdge >::iterator e_it = edges.begin();
typename std::vector< LightEdge >::iterator e_end = edges.end(); typename std::vector< LightEdge >::iterator e_end = edges.end();
FacePointer g, f; FacePointer g, f;
int w, z; int w, z;
for( ; e_it!=e_end; e_it++) for( ; e_it!=e_end; e_it++)
{ {
f = &_mesh->face[e_it->face]; f = &_mesh->face[e_it->face];
z = (int) e_it->edge; z = (int) e_it->edge;
// v2------v1 swap the diagonal only if v2 and v3 are feature and v0 and v1 are not. // v2------v1 swap the diagonal only if v2 and v3 are feature and v0 and v1 are not.
// | / | // | / |
// | / | // | / |
// v0------v3 // v0------v3
if (!(f->IsS()) && vcg::face::CheckFlipEdge< FaceType >(*f, z)) if (!(f->IsS()) && vcg::face::CheckFlipEdge< FaceType >(*f, z))
{ {
VertexPointer v0, v1, v2, v3; VertexPointer v0, v1, v2, v3;
v0 = f->V(z); v0 = f->V(z);
v1 = f->V1(z); v1 = f->V1(z);
v2 = f->V2(z); v2 = f->V2(z);
g = f->FFp(z); g = f->FFp(z);
w = f->FFi(z); w = f->FFi(z);
v3 = g->V2(w); v3 = g->V2(w);
bool b0, b1, b2, b3; bool b0, b1, b2, b3;
b0 = !v0->IsUserBit(_featureFlag) ; b0 = !v0->IsUserBit(_featureFlag) ;
b1 = !v1->IsUserBit(_featureFlag) ; b1 = !v1->IsUserBit(_featureFlag) ;
b2 = v2->IsUserBit(_featureFlag) ; b2 = v2->IsUserBit(_featureFlag) ;
b3 = v3->IsUserBit(_featureFlag) ; b3 = v3->IsUserBit(_featureFlag) ;
if( b0 && b1 && b2 && b3) if( b0 && b1 && b2 && b3)
vcg::face::FlipEdge< FaceType >(*f, z); vcg::face::FlipEdge< FaceType >(*f, z);
} // end if (vcg::face::CheckFlipEdge< _Face >(*f, z)) } // end if (vcg::face::CheckFlipEdge< _Face >(*f, z))
} // end for( ; e_it!=e_end; e_it++) } // end for( ; e_it!=e_end; e_it++)
} //end of FlipEdges } //end of FlipEdges
}; // end of class ExtendedMarchingCubes }; // end of class ExtendedMarchingCubes
// /*! @} */ // /*! @} */
// end of Doxygen documentation // end of Doxygen documentation
} // end of namespace tri } // end of namespace tri
}; // end of namespace vcg }; // end of namespace vcg
#endif // __VCG_EXTENDED_MARCHING_CUBES #endif // __VCG_EXTENDED_MARCHING_CUBES

1299
vcg/complex/algorithms/create/marching_cubes.h
View File

File diff suppressed because it is too large Load Diff

395
vcg/complex/algorithms/create/platonic.h
View File

@ -25,7 +25,6 @@
#define __VCGLIB_PLATONIC #define __VCGLIB_PLATONIC
#include<vcg/math/base.h> #include<vcg/math/base.h>
#include<vcg/complex/allocate.h>
#include<vcg/complex/algorithms/refine.h> #include<vcg/complex/algorithms/refine.h>
#include<vcg/complex/algorithms/update/flag.h> #include<vcg/complex/algorithms/update/flag.h>
#include<vcg/complex/algorithms/update/position.h> #include<vcg/complex/algorithms/update/position.h>
@ -42,9 +41,9 @@ namespace tri {
that represent surfaces of platonic solids, that represent surfaces of platonic solids,
and other simple shapes. and other simple shapes.
The 1st parameter is the mesh that will The 1st parameter is the mesh that will
be filled with the solid. be filled with the solid.
*/ */
template <class TetraMeshType> template <class TetraMeshType>
void Tetrahedron(TetraMeshType &in) void Tetrahedron(TetraMeshType &in)
{ {
@ -135,31 +134,31 @@ void Dodecahedron(DodMeshType & in)
int h,i,j,m=0; int h,i,j,m=0;
bool used[N_points]; bool used[N_points];
for (i=0; i<N_points; i++) used[i]=false; for (i=0; i<N_points; i++) used[i]=false;
int reindex[20+12 *10]; int reindex[20+12 *10];
ScalarType xx,yy,zz, sx,sy,sz; ScalarType xx,yy,zz, sx,sy,sz;
int order[5]={0,1,8,6,2}; int order[5]={0,1,8,6,2};
int added[12]; int added[12];
VertexIterator vi=in.vert.begin(); VertexIterator vi=in.vert.begin();
for (i=0; i<12; i++) { for (i=0; i<12; i++) {
sx=sy=sz=0; sx=sy=sz=0;
for (int j=0; j<5; j++) { for (int j=0; j<5; j++) {
h= penta[ i*9 + order[j] ]-1; h= penta[ i*9 + order[j] ]-1;
xx=vv[h*3];yy=vv[h*3+1];zz=vv[h*3+2]; sx+=xx; sy+=yy; sz+=zz; xx=vv[h*3];yy=vv[h*3+1];zz=vv[h*3+2]; sx+=xx; sy+=yy; sz+=zz;
if (!used[h]) { if (!used[h]) {
(*vi).P()=CoordType( xx, yy, zz ); vi++; (*vi).P()=CoordType( xx, yy, zz ); vi++;
used[h]=true; used[h]=true;
reindex[ h ] = m++; reindex[ h ] = m++;
} }
} }
(*vi).P()=CoordType( sx/5.0, sy/5.0, sz/5.0 ); vi++; (*vi).P()=CoordType( sx/5.0, sy/5.0, sz/5.0 ); vi++;
added[ i ] = m++; added[ i ] = m++;
} }
std::vector<VertexPointer> index(in.vn); std::vector<VertexPointer> index(in.vn);
@ -167,19 +166,19 @@ void Dodecahedron(DodMeshType & in)
FaceIterator fi=in.face.begin(); FaceIterator fi=in.face.begin();
for (i=0; i<12; i++) { for (i=0; i<12; i++) {
for (j=0; j<5; j++){ for (j=0; j<5; j++){
(*fi).V(0)=index[added[i] ]; (*fi).V(0)=index[added[i] ];
(*fi).V(1)=index[reindex[penta[i*9 + order[j ] ] -1 ] ]; (*fi).V(1)=index[reindex[penta[i*9 + order[j ] ] -1 ] ];
(*fi).V(2)=index[reindex[penta[i*9 + order[(j+1)%5] ] -1 ] ]; (*fi).V(2)=index[reindex[penta[i*9 + order[(j+1)%5] ] -1 ] ];
if (HasPerFaceFlags(in)) { if (HasPerFaceFlags(in)) {
// tag faux edges // tag faux edges
(*fi).SetF(0); (*fi).SetF(0);
(*fi).SetF(2); (*fi).SetF(2);
} }
fi++; fi++;
} }
} }
} }
template <class OctMeshType> template <class OctMeshType>
@ -233,24 +232,24 @@ void Icosahedron(IcoMeshType &in)
CoordType ( 0,-L, 1), CoordType ( 0,-L, 1),
CoordType ( 0,-L,-1), CoordType ( 0,-L,-1),
CoordType ( L, 1, 0), CoordType ( L, 1, 0),
CoordType ( L,-1, 0), CoordType ( L,-1, 0),
CoordType (-L, 1, 0), CoordType (-L, 1, 0),
CoordType (-L,-1, 0), CoordType (-L,-1, 0),
CoordType ( 1, 0, L), CoordType ( 1, 0, L),
CoordType (-1, 0, L), CoordType (-1, 0, L),
CoordType ( 1, 0,-L), CoordType ( 1, 0,-L),
CoordType (-1, 0,-L) CoordType (-1, 0,-L)
}; };
int ff[20][3]={ int ff[20][3]={
{1,0,4},{0,1,6},{2,3,5},{3,2,7}, {1,0,4},{0,1,6},{2,3,5},{3,2,7},
{4,5,10},{5,4,8},{6,7,9},{7,6,11}, {4,5,10},{5,4,8},{6,7,9},{7,6,11},
{8,9,2},{9,8,0},{10,11,1},{11,10,3}, {8,9,2},{9,8,0},{10,11,1},{11,10,3},
{0,8,4},{0,6,9},{1,4,10},{1,11,6}, {0,8,4},{0,6,9},{1,4,10},{1,11,6},
{2,5,8},{2,9,7},{3,10,5},{3,7,11} {2,5,8},{2,9,7},{3,10,5},{3,7,11}
}; };
in.Clear(); in.Clear();
@ -368,33 +367,33 @@ void Sphere(MeshType &in, const int subdiv = 3 )
typedef typename MeshType::FaceIterator FaceIterator; typedef typename MeshType::FaceIterator FaceIterator;
if(in.vn==0 && in.fn==0) Icosahedron(in); if(in.vn==0 && in.fn==0) Icosahedron(in);
VertexIterator vi; VertexIterator vi;
for(vi = in.vert.begin(); vi!=in.vert.end();++vi) for(vi = in.vert.begin(); vi!=in.vert.end();++vi)
vi->P().Normalize(); vi->P().Normalize();
tri::UpdateFlags<MeshType>::FaceBorderFromNone(in); tri::UpdateFlags<MeshType>::FaceBorderFromNone(in);
tri::UpdateTopology<MeshType>::FaceFace(in); tri::UpdateTopology<MeshType>::FaceFace(in);
size_t lastsize = 0; size_t lastsize = 0;
for(int i = 0 ; i < subdiv; ++i) for(int i = 0 ; i < subdiv; ++i)
{ {
Refine< MeshType, MidPoint<MeshType> >(in, MidPoint<MeshType>(&in), 0); Refine< MeshType, MidPoint<MeshType> >(in, MidPoint<MeshType>(&in), 0);
for(vi = in.vert.begin() + lastsize; vi != in.vert.end(); ++vi) for(vi = in.vert.begin() + lastsize; vi != in.vert.end(); ++vi)
vi->P().Normalize(); vi->P().Normalize();
lastsize = in.vert.size(); lastsize = in.vert.size();
} }
} }
/// r1 = raggio 1, r2 = raggio2, h = altezza (asse y) /// r1 = raggio 1, r2 = raggio2, h = altezza (asse y)
template <class MeshType> template <class MeshType>
void Cone( MeshType& in, void Cone( MeshType& in,
const typename MeshType::ScalarType r1, const typename MeshType::ScalarType r1,
const typename MeshType::ScalarType r2, const typename MeshType::ScalarType r2,
const typename MeshType::ScalarType h, const typename MeshType::ScalarType h,
const int SubDiv = 36 ) const int SubDiv = 36 )
{ {
typedef typename MeshType::ScalarType ScalarType; typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType; typedef typename MeshType::CoordType CoordType;
@ -433,14 +432,14 @@ void Cone( MeshType& in,
b2 += SubDiv; b2 += SubDiv;
} }
if(r2!=0) if(r2!=0)
{ {
for(i=0;i<SubDiv;++i) for(i=0;i<SubDiv;++i)
{ {
double a = math::ToRad(i*360.0/SubDiv); double a = math::ToRad(i*360.0/SubDiv);
ivp[cnt]=&*vi; (*vi).P()= CoordType( r2*cos(a), h/2.0, r2*sin(a)); ++vi;++cnt; ivp[cnt]=&*vi; (*vi).P()= CoordType( r2*cos(a), h/2.0, r2*sin(a)); ++vi;++cnt;
} }
} }
FaceIterator fi=in.face.begin(); FaceIterator fi=in.face.begin();
@ -456,29 +455,29 @@ void Cone( MeshType& in,
(*fi).V(1)=ivp[b2+(i+1)%SubDiv]; (*fi).V(1)=ivp[b2+(i+1)%SubDiv];
} }
if(r1==0) for(i=0;i<SubDiv;++i,++fi) if(r1==0) for(i=0;i<SubDiv;++i,++fi)
{ {
(*fi).V(0)=ivp[0]; (*fi).V(0)=ivp[0];
(*fi).V(1)=ivp[b2+i]; (*fi).V(1)=ivp[b2+i];
(*fi).V(2)=ivp[b2+(i+1)%SubDiv]; (*fi).V(2)=ivp[b2+(i+1)%SubDiv];
} }
if(r2==0) for(i=0;i<SubDiv;++i,++fi){ if(r2==0) for(i=0;i<SubDiv;++i,++fi){
(*fi).V(0)=ivp[1]; (*fi).V(0)=ivp[1];
(*fi).V(2)=ivp[b1+i]; (*fi).V(2)=ivp[b1+i];
(*fi).V(1)=ivp[b1+(i+1)%SubDiv]; (*fi).V(1)=ivp[b1+(i+1)%SubDiv];
} }
if(r1!=0 && r2!=0)for(i=0;i<SubDiv;++i) if(r1!=0 && r2!=0)for(i=0;i<SubDiv;++i)
{ {
(*fi).V(0)=ivp[b1+i]; (*fi).V(0)=ivp[b1+i];
(*fi).V(1)=ivp[b2+i]; (*fi).V(1)=ivp[b2+i];
(*fi).V(2)=ivp[b2+(i+1)%SubDiv]; (*fi).V(2)=ivp[b2+(i+1)%SubDiv];
++fi; ++fi;
(*fi).V(0)=ivp[b1+i]; (*fi).V(0)=ivp[b1+i];
(*fi).V(1)=ivp[b2+(i+1)%SubDiv]; (*fi).V(1)=ivp[b2+(i+1)%SubDiv];
(*fi).V(2)=ivp[b1+(i+1)%SubDiv]; (*fi).V(2)=ivp[b1+(i+1)%SubDiv];
++fi; ++fi;
} }
} }
@ -654,10 +653,10 @@ void Grid(MeshType & in, int w, int h, float wl, float hl, float *data=0)
template <class MeshType> template <class MeshType>
void FaceGrid(MeshType & in, int w, int h) void FaceGrid(MeshType & in, int w, int h)
{ {
assert(in.vn == (int)in.vert.size()); // require a compact vertex vector assert(in.vn == (int)in.vert.size()); // require a compact vertex vector
assert(in.vn >= w*h); // the number of vertices should match the number of expected grid vertices assert(in.vn >= w*h); // the number of vertices should match the number of expected grid vertices
Allocator<MeshType>::AddFaces(in,(w-1)*(h-1)*2); Allocator<MeshType>::AddFaces(in,(w-1)*(h-1)*2);
// i+0,j+0 -- i+0,j+1 // i+0,j+0 -- i+0,j+1
// | \ | // | \ |
@ -694,8 +693,8 @@ void FaceGrid(MeshType & in, int w, int h)
template <class MeshType> template <class MeshType>
void FaceGrid(MeshType & in, const std::vector<int> &grid, int w, int h) void FaceGrid(MeshType & in, const std::vector<int> &grid, int w, int h)
{ {
assert(in.vn == (int)in.vert.size()); // require a compact vertex vector assert(in.vn == (int)in.vert.size()); // require a compact vertex vector
assert(in.vn <= w*h); // the number of vertices should match the number of expected grid vertices assert(in.vn <= w*h); // the number of vertices should match the number of expected grid vertices
// V0 V1 // V0 V1
// i+0,j+0 -- i+0,j+1 // i+0,j+0 -- i+0,j+1
@ -715,27 +714,27 @@ void FaceGrid(MeshType & in, const std::vector<int> &grid, int w, int h)
int V2i= grid[(i+1)*w+j+0]; int V2i= grid[(i+1)*w+j+0];
int V3i= grid[(i+1)*w+j+1]; int V3i= grid[(i+1)*w+j+1];
int ndone=0; int ndone=0;
bool quad = (V0i>=0 && V1i>=0 && V2i>=0 && V3i>=0 ) && tri::HasPerFaceFlags(in); bool quad = (V0i>=0 && V1i>=0 && V2i>=0 && V3i>=0 ) && tri::HasPerFaceFlags(in);
if(V0i>=0 && V2i>=0 && V3i>=0 ) if(V0i>=0 && V2i>=0 && V3i>=0 )
{ {
typename MeshType::FaceIterator f= Allocator<MeshType>::AddFaces(in,1); typename MeshType::FaceIterator f= Allocator<MeshType>::AddFaces(in,1);
f->V(0)=&(in.vert[V3i]); f->V(0)=&(in.vert[V3i]);
f->V(1)=&(in.vert[V2i]); f->V(1)=&(in.vert[V2i]);
f->V(2)=&(in.vert[V0i]); f->V(2)=&(in.vert[V0i]);
if (quad) f->SetF(2); if (quad) f->SetF(2);
ndone++; ndone++;
} }
if(V0i>=0 && V1i>=0 && V3i>=0 ) if(V0i>=0 && V1i>=0 && V3i>=0 )
{ {
typename MeshType::FaceIterator f= Allocator<MeshType>::AddFaces(in,1); typename MeshType::FaceIterator f= Allocator<MeshType>::AddFaces(in,1);
f->V(0)=&(in.vert[V0i]); f->V(0)=&(in.vert[V0i]);
f->V(1)=&(in.vert[V1i]); f->V(1)=&(in.vert[V1i]);
f->V(2)=&(in.vert[V3i]); f->V(2)=&(in.vert[V3i]);
if (quad) f->SetF(2); if (quad) f->SetF(2);
ndone++; ndone++;
} }
if (ndone==0) { // try diag the other way if (ndone==0) { // try diag the other way
if(V2i>=0 && V0i>=0 && V1i>=0 ) if(V2i>=0 && V0i>=0 && V1i>=0 )
@ -850,49 +849,49 @@ void OrientedDisk(MeshType &m, int slices, Point3f center, Point3f norm, float r
template <class MeshType> template <class MeshType>
void Cylinder(int slices, int stacks, MeshType & m){ void Cylinder(int slices, int stacks, MeshType & m){
typename MeshType::VertexIterator vi = vcg::tri::Allocator<MeshType>::AddVertices(m,slices*(stacks+1)); typename MeshType::VertexIterator vi = vcg::tri::Allocator<MeshType>::AddVertices(m,slices*(stacks+1));
for ( int i = 0; i < stacks+1; ++i) for ( int i = 0; i < stacks+1; ++i)
for ( int j = 0; j < slices; ++j) for ( int j = 0; j < slices; ++j)
{ {
float x,y,h; float x,y,h;
x = cos( 2.0 * M_PI / slices * j); x = cos( 2.0 * M_PI / slices * j);
y = sin( 2.0 * M_PI / slices * j); y = sin( 2.0 * M_PI / slices * j);
h = 2 * i / (float)(stacks) - 1; h = 2 * i / (float)(stacks) - 1;
(*vi).P() = typename MeshType::CoordType(x,h,y); (*vi).P() = typename MeshType::CoordType(x,h,y);
++vi; ++vi;
} }
typename MeshType::FaceIterator fi ; typename MeshType::FaceIterator fi ;
for ( int j = 0; j < stacks; ++j) for ( int j = 0; j < stacks; ++j)
for ( int i = 0; i < slices; ++i) for ( int i = 0; i < slices; ++i)
{ {
int a,b,c,d; int a,b,c,d;
a = (j+0)*slices + i; a = (j+0)*slices + i;
b = (j+1)*slices + i; b = (j+1)*slices + i;
c = (j+1)*slices + (i+1)%slices; c = (j+1)*slices + (i+1)%slices;
d = (j+0)*slices + (i+1)%slices; d = (j+0)*slices + (i+1)%slices;
if(((i+j)%2) == 0){ if(((i+j)%2) == 0){
fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1); fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1);
(*fi).V(0) = &m.vert[ a ]; (*fi).V(0) = &m.vert[ a ];
(*fi).V(1) = &m.vert[ b ]; (*fi).V(1) = &m.vert[ b ];
(*fi).V(2) = &m.vert[ c ]; (*fi).V(2) = &m.vert[ c ];
fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1); fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1);
(*fi).V(0) = &m.vert[ c ]; (*fi).V(0) = &m.vert[ c ];
(*fi).V(1) = &m.vert[ d ]; (*fi).V(1) = &m.vert[ d ];
(*fi).V(2) = &m.vert[ a ]; (*fi).V(2) = &m.vert[ a ];
} }
else{ else{
fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1); fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1);
(*fi).V(0) = &m.vert[ b ]; (*fi).V(0) = &m.vert[ b ];
(*fi).V(1) = &m.vert[ c ]; (*fi).V(1) = &m.vert[ c ];
(*fi).V(2) = &m.vert[ d ]; (*fi).V(2) = &m.vert[ d ];
fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1); fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1);
(*fi).V(0) = &m.vert[ d ]; (*fi).V(0) = &m.vert[ d ];
(*fi).V(1) = &m.vert[ a ]; (*fi).V(1) = &m.vert[ a ];
(*fi).V(2) = &m.vert[ b ]; (*fi).V(2) = &m.vert[ b ];
} }
} }
@ -908,47 +907,47 @@ void Cylinder(int slices, int stacks, MeshType & m){
template <class MeshType> template <class MeshType>
void GenerateCameraMesh(MeshType &in){ void GenerateCameraMesh(MeshType &in){
typedef typename MeshType::CoordType MV; typedef typename MeshType::CoordType MV;
MV vv[52]={ MV vv[52]={
MV(-0.000122145 , -0.2 ,0.35), MV(-0.000122145 , -0.2 ,0.35),
MV(0.000122145 , -0.2 ,-0.35),MV(-0.000122145 , 0.2 ,0.35),MV(0.000122145 , 0.2 ,-0.35),MV(0.999878 , -0.2 ,0.350349),MV(1.00012 , -0.2 ,-0.349651),MV(0.999878 , 0.2 ,0.350349),MV(1.00012 , 0.2 ,-0.349651),MV(1.28255 , 0.1 ,0.754205),MV(1.16539 , 0.1 ,1.03705),MV(0.88255 , 0.1 ,1.15421), MV(0.000122145 , -0.2 ,-0.35),MV(-0.000122145 , 0.2 ,0.35),MV(0.000122145 , 0.2 ,-0.35),MV(0.999878 , -0.2 ,0.350349),MV(1.00012 , -0.2 ,-0.349651),MV(0.999878 , 0.2 ,0.350349),MV(1.00012 , 0.2 ,-0.349651),MV(1.28255 , 0.1 ,0.754205),MV(1.16539 , 0.1 ,1.03705),MV(0.88255 , 0.1 ,1.15421),
MV(0.599707 , 0.1 ,1.03705),MV(0.48255 , 0.1 ,0.754205),MV(0.599707 , 0.1 ,0.471362),MV(0.88255 , 0.1 ,0.354205),MV(1.16539 , 0.1 ,0.471362),MV(1.28255 , -0.1 ,0.754205),MV(1.16539 , -0.1 ,1.03705),MV(0.88255 , -0.1 ,1.15421),MV(0.599707 , -0.1 ,1.03705),MV(0.48255 , -0.1 ,0.754205), MV(0.599707 , 0.1 ,1.03705),MV(0.48255 , 0.1 ,0.754205),MV(0.599707 , 0.1 ,0.471362),MV(0.88255 , 0.1 ,0.354205),MV(1.16539 , 0.1 ,0.471362),MV(1.28255 , -0.1 ,0.754205),MV(1.16539 , -0.1 ,1.03705),MV(0.88255 , -0.1 ,1.15421),MV(0.599707 , -0.1 ,1.03705),MV(0.48255 , -0.1 ,0.754205),
MV(0.599707 , -0.1 ,0.471362),MV(1.16539 , -0.1 ,0.471362),MV(0.88255 , -0.1 ,0.354205),MV(3.49164e-005 , 0 ,-0.1),MV(1.74582e-005 , -0.0866025 ,-0.05),MV(-1.74582e-005 , -0.0866025 ,0.05),MV(-3.49164e-005 , 8.74228e-009 ,0.1),MV(-1.74582e-005 , 0.0866025 ,0.05),MV(1.74582e-005 , 0.0866025 ,-0.05),MV(-0.399913 , 1.99408e-022 ,-0.25014), MV(0.599707 , -0.1 ,0.471362),MV(1.16539 , -0.1 ,0.471362),MV(0.88255 , -0.1 ,0.354205),MV(3.49164e-005 , 0 ,-0.1),MV(1.74582e-005 , -0.0866025 ,-0.05),MV(-1.74582e-005 , -0.0866025 ,0.05),MV(-3.49164e-005 , 8.74228e-009 ,0.1),MV(-1.74582e-005 , 0.0866025 ,0.05),MV(1.74582e-005 , 0.0866025 ,-0.05),MV(-0.399913 , 1.99408e-022 ,-0.25014),
MV(-0.399956 , -0.216506 ,-0.12514),MV(-0.400044 , -0.216506 ,0.12486),MV(-0.400087 , 2.18557e-008 ,0.24986),MV(-0.400044 , 0.216506 ,0.12486),MV(-0.399956 , 0.216506 ,-0.12514),MV(0.479764 , 0.1 ,0.754205),MV(0.362606 , 0.1 ,1.03705),MV(0.0797637 , 0.1 ,1.15421),MV(-0.203079 , 0.1 ,1.03705),MV(-0.320236 , 0.1 ,0.754205), MV(-0.399956 , -0.216506 ,-0.12514),MV(-0.400044 , -0.216506 ,0.12486),MV(-0.400087 , 2.18557e-008 ,0.24986),MV(-0.400044 , 0.216506 ,0.12486),MV(-0.399956 , 0.216506 ,-0.12514),MV(0.479764 , 0.1 ,0.754205),MV(0.362606 , 0.1 ,1.03705),MV(0.0797637 , 0.1 ,1.15421),MV(-0.203079 , 0.1 ,1.03705),MV(-0.320236 , 0.1 ,0.754205),
MV(-0.203079 , 0.1 ,0.471362),MV(0.0797637 , 0.1 ,0.354205),MV(0.362606 , 0.1 ,0.471362),MV(0.479764 , -0.1 ,0.754205),MV(0.362606 , -0.1 ,1.03705),MV(0.0797637 , -0.1 ,1.15421),MV(-0.203079 , -0.1 ,1.03705),MV(-0.320236 , -0.1 ,0.754205),MV(0.0797637 , -0.1 ,0.354205),MV(0.362606 , -0.1 ,0.471362), MV(-0.203079 , 0.1 ,0.471362),MV(0.0797637 , 0.1 ,0.354205),MV(0.362606 , 0.1 ,0.471362),MV(0.479764 , -0.1 ,0.754205),MV(0.362606 , -0.1 ,1.03705),MV(0.0797637 , -0.1 ,1.15421),MV(-0.203079 , -0.1 ,1.03705),MV(-0.320236 , -0.1 ,0.754205),MV(0.0797637 , -0.1 ,0.354205),MV(0.362606 , -0.1 ,0.471362),
MV(-0.203079 , -0.1 ,0.471362), }; MV(-0.203079 , -0.1 ,0.471362), };
int ff[88][3]={ int ff[88][3]={
{0,2,3}, {0,2,3},
{3,1,0},{4,5,7},{7,6,4},{0,1,5},{5,4,0},{1,3,7},{7,5,1},{3,2,6},{6,7,3},{2,0,4}, {3,1,0},{4,5,7},{7,6,4},{0,1,5},{5,4,0},{1,3,7},{7,5,1},{3,2,6},{6,7,3},{2,0,4},
{4,6,2},{10,9,8},{10,12,11},{10,13,12},{10,14,13},{10,15,14},{10,8,15},{8,17,16},{8,9,17},{9,18,17}, {4,6,2},{10,9,8},{10,12,11},{10,13,12},{10,14,13},{10,15,14},{10,8,15},{8,17,16},{8,9,17},{9,18,17},
{9,10,18},{10,19,18},{10,11,19},{11,20,19},{11,12,20},{12,21,20},{12,13,21},{13,23,21},{13,14,23},{14,22,23}, {9,10,18},{10,19,18},{10,11,19},{11,20,19},{11,12,20},{12,21,20},{12,13,21},{13,23,21},{13,14,23},{14,22,23},
{14,15,22},{15,16,22},{15,8,16},{23,16,17},{23,17,18},{23,18,19},{23,19,20},{23,20,21},{23,22,16},{25,27,26}, {14,15,22},{15,16,22},{15,8,16},{23,16,17},{23,17,18},{23,18,19},{23,19,20},{23,20,21},{23,22,16},{25,27,26},
{25,28,27},{25,29,28},{25,24,29},{24,31,30},{24,25,31},{25,32,31},{25,26,32},{26,33,32},{26,27,33},{27,34,33}, {25,28,27},{25,29,28},{25,24,29},{24,31,30},{24,25,31},{25,32,31},{25,26,32},{26,33,32},{26,27,33},{27,34,33},
{27,28,34},{28,35,34},{28,29,35},{29,30,35},{29,24,30},{35,30,31},{35,31,32},{35,32,33},{35,33,34},{42,37,36}, {27,28,34},{28,35,34},{28,29,35},{29,30,35},{29,24,30},{35,30,31},{35,31,32},{35,32,33},{35,33,34},{42,37,36},
{42,38,37},{42,39,38},{42,40,39},{42,41,40},{42,36,43},{36,45,44},{36,37,45},{37,46,45},{37,38,46},{38,47,46}, {42,38,37},{42,39,38},{42,40,39},{42,41,40},{42,36,43},{36,45,44},{36,37,45},{37,46,45},{37,38,46},{38,47,46},
{38,39,47},{39,48,47},{39,40,48},{40,51,48},{40,41,51},{41,49,51},{41,42,49},{42,50,49},{42,43,50},{43,44,50}, {38,39,47},{39,48,47},{39,40,48},{40,51,48},{40,41,51},{41,49,51},{41,42,49},{42,50,49},{42,43,50},{43,44,50},
{43,36,44},{51,44,45},{51,45,46},{51,46,47},{51,47,48},{51,49,50},{51,50,44}, {43,36,44},{51,44,45},{51,45,46},{51,46,47},{51,47,48},{51,49,50},{51,50,44},
}; };
in.Clear(); in.Clear();
Allocator<MeshType>::AddVertices(in,52); Allocator<MeshType>::AddVertices(in,52);
Allocator<MeshType>::AddFaces(in,88); Allocator<MeshType>::AddFaces(in,88);
in.vn=52;in.fn=88; in.vn=52;in.fn=88;
int i,j; int i,j;
for(i=0;i<in.vn;i++) for(i=0;i<in.vn;i++)
in.vert[i].P()=vv[i];; in.vert[i].P()=vv[i];;
std::vector<typename MeshType::VertexPointer> index(in.vn); std::vector<typename MeshType::VertexPointer> index(in.vn);
typename MeshType::VertexIterator vi; typename MeshType::VertexIterator vi;
for(j=0,vi=in.vert.begin();j<in.vn;++j,++vi) index[j] = &*vi; for(j=0,vi=in.vert.begin();j<in.vn;++j,++vi) index[j] = &*vi;
for(j=0;j<in.fn;++j) for(j=0;j<in.fn;++j)
{ {
in.face[j].V(0)=index[ff[j][0]]; in.face[j].V(0)=index[ff[j][0]];
in.face[j].V(1)=index[ff[j][1]]; in.face[j].V(1)=index[ff[j][1]];
in.face[j].V(2)=index[ff[j][2]]; in.face[j].V(2)=index[ff[j][2]];
} }
} }
template <class MeshType> template <class MeshType>

399
vcg/complex/algorithms/create/zonohedron.h
View File

@ -24,33 +24,28 @@
#ifndef __VCGLIB_ZONOHEDRON #ifndef __VCGLIB_ZONOHEDRON
#define __VCGLIB_ZONOHEDRON #define __VCGLIB_ZONOHEDRON
#include<vcg/complex/allocate.h>
#include<map>
typedef unsigned int uint;
namespace vcg { namespace vcg {
namespace tri { namespace tri {
/** \addtogroup trimesh */ /** \addtogroup trimesh */
//@{ //@{
/** /**
A class to build a Zonohedron. A class to build a Zonohedron.
Given a set of input vectors, a zonohedron is defined Given a set of input vectors, a zonohedron is defined
as the convex hull of all the points which can be costructed by summing as the convex hull of all the points which can be costructed by summing
together any subset of input vectors. together any subset of input vectors.
The surface closing this solid is composed only of flat parallelograms, The surface closing this solid is composed only of flat parallelograms,
(which have the input vectors as sides). (which have the input vectors as sides).
It is always point-symmetric. It is always point-symmetric.
Mesh created by this class are pure-quad meshes (triangular bit-quad), Mesh created by this class are pure-quad meshes (triangular bit-quad),
(when coplanar vectors are fed, then planar groups of quads can be seen as (when coplanar vectors are fed, then planar groups of quads can be seen as
forming planar faces with more than 4 vertices). forming planar faces with more than 4 vertices).
USAGE: USAGE:
1) Instantiate a Zonohedron. 1) Instantiate a Zonohedron.
2) Add input vectors at will to it, with addVector(s) 2) Add input vectors at will to it, with addVector(s)
3) When you are done, call createMesh. 3) When you are done, call createMesh.
*/ */
@ -59,185 +54,185 @@ template <class Scalar>
class Zonohedron{ class Zonohedron{
public: public:
typedef Point3<Scalar> Vec3; typedef Point3<Scalar> Vec3;
Zonohedron(){} Zonohedron(){}
void addVector(Scalar x, Scalar y, Scalar z); void addVector(Scalar x, Scalar y, Scalar z);
void addVector(Vec3 v); void addVector(Vec3 v);
void addVectors(const std::vector< Vec3 > ); void addVectors(const std::vector< Vec3 > );
const std::vector< Vec3 >& vectors() const { const std::vector< Vec3 >& vectors() const {
return vec; return vec;
} }
template<class MeshType> template<class MeshType>
void createMesh( MeshType& output ); void createMesh( MeshType& output );
private: private:
/* classes for internal use */ /* classes for internal use */
/****************************/ /****************************/
typedef int VecIndex; // a number in [0..n) typedef int VecIndex; // a number in [0..n)
/* the signature of a vertex (a 0 or 1 per input vector) */ /* the signature of a vertex (a 0 or 1 per input vector) */
struct Signature { struct Signature {
std::vector< bool > v; std::vector< bool > v;
Signature(){} Signature(){}
Signature(int n){ v.resize(n,false); } Signature(int n){ v.resize(n,false); }
bool operator == (const Signature & b) const { bool operator == (const Signature & b) const {
return (b.v == v); return (b.v == v);
} }
bool operator < (const Signature & b) const { bool operator < (const Signature & b) const {
return (b.v < v); return (b.v < v);
} }
Signature& set(VecIndex i, bool value){ Signature& set(VecIndex i, bool value){
v[i] = value; v[i] = value;
return *this; return *this;
} }
Signature& set(VecIndex i, bool valueI, VecIndex j, bool valueJ){ Signature& set(VecIndex i, bool valueI, VecIndex j, bool valueJ){
v[i] = valueI; v[i] = valueI;
v[j] = valueJ; v[j] = valueJ;
return *this; return *this;
} }
}; };
struct Face { struct Face {
int vert[4]; // index to vertex array int vert[4]; // index to vertex array
}; };
/* precomputed cross products for all pairs of vectors */ /* precomputed cross products for all pairs of vectors */
std::vector< Vec3 > precomputedCross; std::vector< Vec3 > precomputedCross;
void precompteAllCrosses(){ void precompteAllCrosses(){
precomputedCross.resize(n*n); precomputedCross.resize(n*n);
for (int i=0; i<n; i++) for (int j=0; j<n; j++) { for (int i=0; i<n; i++) for (int j=0; j<n; j++) {
precomputedCross[i*n+j] = vec[i] ^ vec[j] ; precomputedCross[i*n+j] = vec[i] ^ vec[j] ;
} }
} }
Vec3 cross(VecIndex i, VecIndex j){ Vec3 cross(VecIndex i, VecIndex j){
return precomputedCross[i*n+j]; return precomputedCross[i*n+j];
} }
// given a vector, returns a copy pointing a unique verse // given a vector, returns a copy pointing a unique verse
static Vec3 uniqueVerse(Vec3 v){ static Vec3 uniqueVerse(Vec3 v){
if (v.X()>0) return v; if (v.X()>0) return v;
else if (v.X()<0) return -v; else if (v.X()<0) return -v;
else if (v.Y()>0) return v; else if (v.Y()>0) return v;
else if (v.Y()<0) return -v; else if (v.Y()<0) return -v;
else if (v.Z()>0) return v; else if (v.Z()>0) return v;
return -v; return -v;
} }
static Vec3 altVec(int i) { static Vec3 altVec(int i) {
return Vec3(1, i, i*i); return Vec3(1, i, i*i);
} }
static Scalar tripleProduct( const Vec3 &a, const Vec3 &b, const Vec3 & c){ static Scalar tripleProduct( const Vec3 &a, const Vec3 &b, const Vec3 & c){
return ( a ^ b ) * c; return ( a ^ b ) * c;
} }
// returns signof: (i x j) * k // returns signof: (i x j) * k
bool signOf_IxJoK(VecIndex i, VecIndex j, VecIndex k){ bool signOf_IxJoK(VecIndex i, VecIndex j, VecIndex k){
const float EPSILON_SQUARED = 1e-12; const float EPSILON_SQUARED = 1e-12;
bool invert = false; bool invert = false;
// sort i,j,k // sort i,j,k
if (i<j) { std::swap(i,j); invert = !invert; } if (i<j) { std::swap(i,j); invert = !invert; }
if (j<k) { std::swap(j,k); invert = !invert; if (j<k) { std::swap(j,k); invert = !invert;
if (i<j) { std::swap(i,j); invert = !invert; } if (i<j) { std::swap(i,j); invert = !invert; }
} }
//Scalar res = Vec3::dot( Vec3::cross( vec[i] , vec[j] ) , vec[k] ); //Scalar res = Vec3::dot( Vec3::cross( vec[i] , vec[j] ) , vec[k] );
Scalar res = cross( i , j ) * vec[k] ; Scalar res = cross( i , j ) * vec[k] ;
if (res*res<=EPSILON_SQUARED) { if (res*res<=EPSILON_SQUARED) {
// three coplanar vectors! // three coplanar vectors!
// use derivative... // use derivative...
//res = uniqueVerse( cross(i,j) ) * cross(j,k) ; //res = uniqueVerse( cross(i,j) ) * cross(j,k) ;
res = tripleProduct( altVec(i), vec[j], vec[k]) + res = tripleProduct( altVec(i), vec[j], vec[k]) +
tripleProduct( vec[i], altVec(j), vec[k]) + tripleProduct( vec[i], altVec(j), vec[k]) +
tripleProduct( vec[i], vec[j], altVec(k)) ; tripleProduct( vec[i], vec[j], altVec(k)) ;
if (res*res<=EPSILON_SQUARED) { if (res*res<=EPSILON_SQUARED) {
// zero derivative (happens, if three colinear vectors, or...) // zero derivative (happens, if three colinear vectors, or...)
res = tripleProduct( vec[i], altVec(j), altVec(k)) + res = tripleProduct( vec[i], altVec(j), altVec(k)) +
tripleProduct( altVec(i), vec[j], altVec(k)) + tripleProduct( altVec(i), vec[j], altVec(k)) +
tripleProduct( altVec(i), altVec(j), vec[k]) ; tripleProduct( altVec(i), altVec(j), vec[k]) ;
} }
if (res*res<=EPSILON_SQUARED) { if (res*res<=EPSILON_SQUARED) {
// zero second derivative (happens if three zero-vectors, i.e. never? or...) // zero second derivative (happens if three zero-vectors, i.e. never? or...)
res = tripleProduct( altVec(i), altVec(j), altVec(k) ); res = tripleProduct( altVec(i), altVec(j), altVec(k) );
} }
} }
return ( (res>=0) != invert ); // XOR return ( (res>=0) != invert ); // XOR
} }
int n; // number of input vectors int n; // number of input vectors
std::vector<Vec3> vec; // input vectors std::vector<Vec3> vec; // input vectors
int vertCount; int vertCount;
std::vector<Face> _face; std::vector<Face> _face;
typedef std::map< Signature, int > VertexMap; typedef std::map< Signature, int > VertexMap;
VertexMap vertexMap; VertexMap vertexMap;
// given a vertex signature, returns index of vert (newly created or not) // given a vertex signature, returns index of vert (newly created or not)
VecIndex vertexIndex(const Signature &s){ VecIndex vertexIndex(const Signature &s){
typename VertexMap::iterator i; typename VertexMap::iterator i;
//Vec3 pos = s; //toPos(s); //Vec3 pos = s; //toPos(s);
i = vertexMap.find( s ); i = vertexMap.find( s );
if (i!= vertexMap.end() ) return i->second; if (i!= vertexMap.end() ) return i->second;
else { else {
int newVertex = vertCount++; int newVertex = vertCount++;
//vertexMap.insert(s) //vertexMap.insert(s)
vertexMap[s] = newVertex; vertexMap[s] = newVertex;
return newVertex; return newVertex;
} }
} }
// given two index of vectors, returns face // given two index of vectors, returns face
Face& face(VecIndex i, VecIndex j){ Face& face(VecIndex i, VecIndex j){
assert(i!=j); assert(i!=j);
assert( i*n + j < (int) _face.size() ); assert( i*n + j < (int) _face.size() );
return _face[i*n + j]; return _face[i*n + j];
} }
Vec3 toPos(const Signature &s) const{ Vec3 toPos(const Signature &s) const{
Vec3 res(0,0,0); Vec3 res(0,0,0);
for (int i=0; i<n; i++) for (int i=0; i<n; i++)
if (s.v[i]) res += vec[i]; if (s.v[i]) res += vec[i];
return res; return res;
} }
void createInternalMesh() { void createInternalMesh() {
n = vec.size(); n = vec.size();
precompteAllCrosses(); precompteAllCrosses();
// allocate faces // allocate faces
_face.resize( n*n ); _face.resize( n*n );
vertCount = 0; vertCount = 0;
vertexMap.clear(); vertexMap.clear();
for (int i=0; i<n; i++) { for (int i=0; i<n; i++) {
//::showProgress(i,n); //::showProgress(i,n);
for (int j=0; j<n; j++) if(i!=j) { for (int j=0; j<n; j++) if(i!=j) {
Signature s(n); Signature s(n);
for (int k=0; k<n; k++) if ((k!=j) && (k!=i)) for (int k=0; k<n; k++) if ((k!=j) && (k!=i))
{ {
s.set( k , signOf_IxJoK( i,j,k ) ); s.set( k , signOf_IxJoK( i,j,k ) );
} }
face(i,j).vert[0] = vertexIndex( s.set(i,false, j,false) ); face(i,j).vert[0] = vertexIndex( s.set(i,false, j,false) );
face(i,j).vert[1] = vertexIndex( s.set(i,false, j,true ) ); face(i,j).vert[1] = vertexIndex( s.set(i,false, j,true ) );
face(i,j).vert[2] = vertexIndex( s.set(i,true, j,true ) ); face(i,j).vert[2] = vertexIndex( s.set(i,true, j,true ) );
face(i,j).vert[3] = vertexIndex( s.set(i,true, j,false) ); face(i,j).vert[3] = vertexIndex( s.set(i,true, j,false) );
} }
} }
} }
}; };
@ -245,64 +240,64 @@ private:
template<class Scalar> template<class Scalar>
void Zonohedron<Scalar>::addVectors(std::vector< Zonohedron<Scalar>::Vec3 > input){ void Zonohedron<Scalar>::addVectors(std::vector< Zonohedron<Scalar>::Vec3 > input){
for (uint i=0; i<input.size(); i++) { for (size_t i=0; i<input.size(); i++) {
addVector( input[i]); addVector( input[i]);
} }
} }
template<class Scalar> template<class Scalar>
void Zonohedron<Scalar>::addVector(Scalar x, Scalar y, Scalar z) { void Zonohedron<Scalar>::addVector(Scalar x, Scalar y, Scalar z) {
addVector( Vec3(x,y,z) ); addVector( Vec3(x,y,z) );
} }
template<class Scalar> template<class Scalar>
void Zonohedron<Scalar>::addVector(Zonohedron<Scalar>::Vec3 v){ void Zonohedron<Scalar>::addVector(Zonohedron<Scalar>::Vec3 v){
vec.push_back(v); vec.push_back(v);
} }
template<class Scalar> template<class Scalar>
template<class MeshType> template<class MeshType>
void Zonohedron<Scalar>::createMesh(MeshType &m){ void Zonohedron<Scalar>::createMesh(MeshType &m){
typedef MeshType Mesh; typedef MeshType Mesh;
typedef typename Mesh::VertexPointer MeshVertexPointer; typedef typename Mesh::VertexPointer MeshVertexPointer;
typedef typename Mesh::VertexIterator MeshVertexIterator; typedef typename Mesh::VertexIterator MeshVertexIterator;
typedef typename Mesh::FaceIterator MeshFaceIterator; typedef typename Mesh::FaceIterator MeshFaceIterator;
typedef typename Mesh::FaceType MeshFace; typedef typename Mesh::FaceType MeshFace;
createInternalMesh(); createInternalMesh();
m.Clear(); m.Clear();
Allocator<MeshType>::AddVertices(m,vertexMap.size()); Allocator<MeshType>::AddVertices(m,vertexMap.size());
Allocator<MeshType>::AddFaces(m,n*(n-1) * 2); Allocator<MeshType>::AddFaces(m,n*(n-1) * 2);
// assign vertex positions // assign vertex positions
MeshVertexIterator vi=m.vert.begin(); MeshVertexIterator vi=m.vert.begin();
for (typename VertexMap::iterator i=vertexMap.begin(); i!=vertexMap.end(); i++){ for (typename VertexMap::iterator i=vertexMap.begin(); i!=vertexMap.end(); i++){
(vi + i->second )->P() = toPos( i->first ); (vi + i->second )->P() = toPos( i->first );
} }
// assegn FV connectivity // assegn FV connectivity
MeshFaceIterator fi=m.face.begin(); MeshFaceIterator fi=m.face.begin();
for (int i=0; i<n; i++) { for (int i=0; i<n; i++) {
for (int j=0; j<n; j++) if (i!=j) { for (int j=0; j<n; j++) if (i!=j) {
const Face &f( face(i,j) ); const Face &f( face(i,j) );
for (int k=0; k<2; k++) { // two tri faces per quad for (int k=0; k<2; k++) { // two tri faces per quad
for (int w=0; w<3; w++) { for (int w=0; w<3; w++) {
fi->V(w) = &* (vi + f.vert[(w+k*2)%4] ); fi->V(w) = &* (vi + f.vert[(w+k*2)%4] );
} }
if (tri::HasPerFaceNormal(m)) { if (tri::HasPerFaceNormal(m)) {
fi->N() = cross(i,j).normalized(); fi->N() = cross(i,j).normalized();
} }
if (tri::HasPerFaceFlags(m)) { if (tri::HasPerFaceFlags(m)) {
fi->SetF(2); // quad diagonals are faux fi->SetF(2); // quad diagonals are faux
} }
fi++; fi++;
} }
} }
} }
} }