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

Ongoing Rearrangement of filepath 

delete old trimesh content
This commit is contained in:
ganovelli 2011-04-01 17:17:15 +00:00
parent f4a5512500
commit 3c7efa7bff
10 changed files with 0 additions and 7555 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

551
vcg/complex/trimesh/create/advancing_front.h
View File

@ -1,551 +0,0 @@
#ifndef MLS_ADVANCE_H
#define MLS_ADVANCE_H
#include <iostream>
#include <list>
#include <wrap/callback.h>
#include <vcg/complex/trimesh/update/topology.h>
#include <vcg/complex/trimesh/update/flag.h>
#include <map>
namespace vcg {
namespace tri {
class FrontEdge {
public:
int v0, v1, v2; //v0, v1 represent the FrontEdge, v2 the other vertex
//in the face this FrontEdge belongs to
int face; //index of the face
bool active; //keep tracks of wether it is in front or in deads
//the loops in the front are mantained as a double linked list
std::list<FrontEdge>::iterator next;
std::list<FrontEdge>::iterator previous;
FrontEdge() {}
FrontEdge(int _v0, int _v1, int _v2, int _face):
v0(_v0), v1(_v1), v2(_v2), face(_face), active(true) {
assert(v0 != v1 && v1 != v2 && v0 != v2);
}
const bool operator==(const FrontEdge& f) const
{
return ((v0 == f.v0) && (v1 == f.v1) && (v2 == f.v2) && (face == f.face));
}
};
template <class MESH> class AdvancingFront {
public:
typedef typename MESH::VertexType VertexType;
typedef typename MESH::FaceType FaceType;
typedef typename MESH::ScalarType ScalarType;
typedef typename MESH::VertexType::CoordType Point3x;
//class FrontEdgeLists
//{
//};
// protected:
std::list<FrontEdge> front;
std::list<FrontEdge> deads;
std::vector<int> nb; //number of fronts a vertex is into,
//this is used for the Visited and Border flags
//but adding topology may not be needed anymore
public:
MESH &mesh; //this structure will be filled by the algorithm
AdvancingFront(MESH &_mesh): mesh(_mesh) {
UpdateFlags<MESH>::FaceBorderFromNone(mesh);
UpdateFlags<MESH>::VertexBorderFromFace(mesh);
nb.clear();
nb.resize(mesh.vert.size(), 0);
CreateLoops();
}
virtual ~AdvancingFront() {}
virtual ScalarType radi() { return 0; }
void BuildMesh(CallBackPos call = NULL, int interval = 512) {
float finalfacesext = mesh.vert.size() * 2.0f;
if(call) call(0, "Advancing front");
while(1) {
for(int i = 0; i < interval; i++) {
if(!front.size() && !SeedFace()) return;
AddFace();
if(call)
{
float rap = float(mesh.face.size()) / finalfacesext;
int perc = (int) (100.0f * rap);
(*call)(perc,"Adding Faces");
}
}
}
}
protected:
//Implement these functions in your subclass
enum ListID {FRONT,DEADS};
typedef std::pair< ListID,std::list<FrontEdge>::iterator > ResultIterator;
virtual bool Seed(int &v0, int &v1, int &v2) = 0;
virtual int Place(FrontEdge &e, ResultIterator &touch) = 0;
bool CheckFrontEdge(int v0, int v1) {
int tot = 0;
//HACK to speed up things until i can use a seach structure
// int i = mesh.face.size() - 4*(front.size());
// if(front.size() < 100) i = mesh.face.size() - 100;
int i = 0;
if(i < 0) i = 0;
for(; i < (int)mesh.face.size(); i++) {
FaceType &f = mesh.face[i];
for(int k = 0; k < 3; k++) {
if(v1== (int)f.V(k) && v0 == (int)f.V((k+1)%3)) ++tot;
else if(v0 == (int)f.V(k) && v1 == (int)f.V((k+1)%3)) { //orientation non constistent
return false;
}
}
if(tot >= 2) { //non manifold
return false;
}
}
return true;
}
//create the FrontEdge loops from seed faces
void CreateLoops() {
VertexType *start = &*mesh.vert.begin();
for(int i = 0; i < (int)mesh.face.size(); i++) {
FaceType &f = mesh.face[i];
if(f.IsD()) continue;
for(int k = 0; k < 3; k++) {
if(f.IsB(k)) {
NewEdge(FrontEdge(f.V0(k) - start, f.V1(k) - start, f.V2(k) - start, i));
nb[f.V0(k)-start]++;
}
}
}
for(std::list<FrontEdge>::iterator s = front.begin(); s != front.end(); s++) {
(*s).previous = front.end();
(*s).next = front.end();
}
//now create loops:
for(std::list<FrontEdge>::iterator s = front.begin(); s != front.end(); s++) {
for(std::list<FrontEdge>::iterator j = front.begin(); j != front.end(); j++) {
if(s == j) continue;
if((*s).v1 != (*j).v0) continue;
if((*j).previous != front.end()) continue;
(*s).next = j;
(*j).previous = s;
break;
}
}
for(std::list<FrontEdge>::iterator s = front.begin(); s != front.end(); s++) {
assert((*s).next != front.end());
assert((*s).previous != front.end());
}
}
bool SeedFace() {
int v[3];
bool success = Seed(v[0], v[1], v[2]);
if(!success) return false;
nb.resize(mesh.vert.size(), 0);
//create the border of the first face
std::list<FrontEdge>::iterator e = front.end();
std::list<FrontEdge>::iterator last = e;
std::list<FrontEdge>::iterator first;
for(int i = 0; i < 3; i++) {
int v0 = v[i];
int v1 = v[((i+1)%3)];
int v2 = v[((i+2)%3)];
mesh.vert[v0].SetB();
nb[v[i]]++;
e = front.insert(front.begin(), FrontEdge(v0, v1, v2, mesh.face.size()));
if(i != 0) {
(*last).next = e;
(*e).previous = last;
} else
first = e;
last = e;
}
//connect last and first
(*last).next = first;
(*first).previous = last;
AddFace(v[0], v[1], v[2]);
return true;
}
public:
bool AddFace() {
if(!front.size()) return false;
std::list<FrontEdge>::iterator ei = front.begin();
FrontEdge &current = *ei;
FrontEdge &previous = *current.previous;
FrontEdge &next = *current.next;
int v0 = current.v0, v1 = current.v1;
assert(nb[v0] < 10 && nb[v1] < 10);
ResultIterator touch;
touch.first = FRONT;
touch.second = front.end();
int v2 = Place(current, touch);
if(v2 == -1) {
KillEdge(ei);
return false;
}
assert(v2 != v0 && v2 != v1);
if ((touch.first == FRONT) && (touch.second != front.end()) ||
(touch.first == DEADS) && (touch.second != deads.end()))
{
//check for orientation and manifoldness
//touch == current.previous?
if(v2 == previous.v0) {
if(!CheckEdge(v2, v1)) {
KillEdge(ei);
return false;
}
/*touching previous FrontEdge (we reuse previous)
next
------->v2 -----> v1------>
\ /
\ /
previous \ / current
\ /
v0 */
Detach(v0);
std::list<FrontEdge>::iterator up = NewEdge(FrontEdge(v2, v1, v0, mesh.face.size()));
MoveFront(up);
(*up).previous = previous.previous;
(*up).next = current.next;
(*previous.previous).next = up;
next.previous = up;
Erase(current.previous);
Erase(ei);
Glue(up);
//touch == (*current.next).next
} else if(v2 == next.v1) {
if(!CheckEdge(v0, v2)) {
KillEdge(ei);
return false;
}
/*touching next FrontEdge (we reuse next)
previous
------->v0 -----> v2------>
\ /
\ /
\ / next
\ /
v1 */
Detach(v1);
std::list<FrontEdge>::iterator up = NewEdge(FrontEdge(v0, v2, v1, mesh.face.size()));
MoveFront(up);
(*up).previous = current.previous;
(*up).next = (*current.next).next;
previous.next = up;
(*next.next).previous = up;
Erase(current.next);
Erase(ei);
Glue(up);
} else {
if(!CheckEdge(v0, v2) || !CheckEdge(v2, v1)) {
KillEdge(ei);
return false;
}
//touching some loop: split (or merge it is local does not matter.
//like this
/*
left right
<--------v2-<------
/|\
/ \
up / \ down
/ \
/ V
----v0 - - - > v1---------
current */
std::list<FrontEdge>::iterator left = touch.second;
std::list<FrontEdge>::iterator right = (*touch.second).previous;
//this would be a really bad join
if(v1 == (*right).v0 || v0 == (*left).v1) {
KillEdge(ei);
return false;
}
nb[v2]++;
std::list<FrontEdge>::iterator down = NewEdge(FrontEdge(v2, v1, v0, mesh.face.size()));
std::list<FrontEdge>::iterator up = NewEdge(FrontEdge(v0, v2, v1, mesh.face.size()));
(*right).next = down;
(*down).previous = right;
(*down).next = current.next;
next.previous = down;
(*left).previous = up;
(*up).next = left;
(*up).previous = current.previous;
previous.next = up;
Erase(ei);
}
}
else if ((touch.first == FRONT) && (touch.second == front.end()) ||
(touch.first == DEADS) && (touch.second == deads.end()))
{
// assert(CheckEdge(v0, v2));
// assert(CheckEdge(v2, v1));
/* adding a new vertex
v2
/|\
/ \
up / \ down
/ \
/ V
----v0 - - - > v1--------- */
assert(!mesh.vert[v2].IsB()); //fatal error! a new point is already a border?
nb[v2]++;
mesh.vert[v2].SetB();
std::list<FrontEdge>::iterator down = NewEdge(FrontEdge(v2, v1, v0, mesh.face.size()));
std::list<FrontEdge>::iterator up = NewEdge(FrontEdge(v0, v2, v1, mesh.face.size()));
(*down).previous = up;
(*up).next = down;
(*down).next = current.next;
next.previous = down;
(*up).previous = current.previous;
previous.next = up;
Erase(ei);
}
AddFace(v0, v2, v1);
return false;
}
protected:
void AddFace(int v0, int v1, int v2) {
assert(v0 < (int)mesh.vert.size() && v1 < (int)mesh.vert.size() && v2 < (int)mesh.vert.size());
FaceType face;
face.V(0) = &mesh.vert[v0];
face.V(1) = &mesh.vert[v1];
face.V(2) = &mesh.vert[v2];
ComputeNormalizedNormal(face);
mesh.face.push_back(face);
mesh.fn++;
}
void AddVertex(VertexType &vertex) {
VertexType *oldstart = NULL;
if(mesh.vert.size()) oldstart = &*mesh.vert.begin();
mesh.vert.push_back(vertex);
mesh.vn++;
VertexType *newstart = &*mesh.vert.begin();
if(oldstart && oldstart != newstart) {
for(int i = 0; i < mesh.face.size(); i++) {
FaceType &face = mesh.face[i];
for(int k = 0; k < 3; k++)
face.V(k) = newstart + (face.V(k) - oldstart);
}
}
nb.push_back(0);
}
bool CheckEdge(int v0, int v1) {
int tot = 0;
//HACK to speed up things until i can use a seach structure
/* int i = mesh.face.size() - 4*(front.size());
if(front.size() < 100) i = mesh.face.size() - 100;
if(i < 0) i = 0;*/
VertexType *vv0 = &(mesh.vert[v0]);
VertexType *vv1 = &(mesh.vert[v1]);
for(int i = 0; i < (int)mesh.face.size(); i++) {
FaceType &f = mesh.face[i];
for(int k = 0; k < 3; k++) {
if(vv0 == f.V0(k) && vv1 == f.V1(k)) //orientation non constistent
return false;
else if(vv1 == f.V0(k) && vv0 == f.V1(k)) ++tot;
}
if(tot >= 2) { //non manifold
return false;
}
}
return true;
}
//front management:
//Add a new FrontEdge to the back of the queue
std::list<FrontEdge>::iterator NewEdge(FrontEdge e) {
return front.insert(front.end(), e);
}
//move an Edge among the dead ones
void KillEdge(std::list<FrontEdge>::iterator e)
{
if (e->active)
{
(*e).active = false;
//std::list<FrontEdge>::iterator res = std::find(front.begin(),front.end(),e);
FrontEdge tmp = *e;
deads.splice(deads.end(), front, e);
std::list<FrontEdge>::iterator newe = std::find(deads.begin(),deads.end(),tmp);
tmp.previous->next = newe;
tmp.next->previous = newe;
}
}
void Erase(std::list<FrontEdge>::iterator e) {
if((*e).active) front.erase(e);
else deads.erase(e);
}
//move an FrontEdge to the back of the queue
void MoveBack(std::list<FrontEdge>::iterator e) {
front.splice(front.end(), front, e);
}
void MoveFront(std::list<FrontEdge>::iterator e) {
front.splice(front.begin(), front, e);
}
//check if e can be sewed with one of oits neighbours
bool Glue(std::list<FrontEdge>::iterator e) {
return Glue((*e).previous, e) || Glue(e, (*e).next);
}
//Glue toghether a and b (where a.next = b
bool Glue(std::list<FrontEdge>::iterator a, std::list<FrontEdge>::iterator b) {
if((*a).v0 != (*b).v1) return false;
std::list<FrontEdge>::iterator previous = (*a).previous;
std::list<FrontEdge>::iterator next = (*b).next;
(*previous).next = next;
(*next).previous = previous;
Detach((*a).v1);
Detach((*a).v0);
Erase(a);
Erase(b);
return true;
}
void Detach(int v) {
assert(nb[v] > 0);
if(--nb[v] == 0) {
mesh.vert[v].ClearB();
}
}
};
template <class MESH> class AdvancingTest: public AdvancingFront<MESH> {
public:
typedef typename MESH::VertexType VertexType;
typedef typename MESH::VertexIterator VertexIterator;
typedef typename MESH::FaceType FaceType;
typedef typename MESH::FaceIterator FaceIterator;
typedef typename MESH::ScalarType ScalarType;
typedef typename MESH::VertexType::CoordType Point3x;
AdvancingTest(MESH &_mesh): AdvancingFront<MESH>(_mesh) {}
bool Seed(int &v0, int &v1, int &v2) {
VertexType v[3];
v[0].P() = Point3x(0, 0, 0);
v[1].P() = Point3x(1, 0, 0);
v[2].P() = Point3x(0, 1, 0);
v[0].ClearFlags();
v[1].ClearFlags();
v[2].ClearFlags();
v0 = this->mesh.vert.size();
AddVertex(v[0]);
v1 = this->mesh.vert.size();
AddVertex(v[1]);
v2 = this->mesh.vert.size();
AddVertex(v[2]);
return true;
}
int Place(FrontEdge &e, typename AdvancingFront<MESH>::ResultIterator &touch)
{
Point3f p[3];
p[0] = this->mesh.vert[e.v0].P();
p[1] = this->mesh.vert[e.v1].P();
p[2] = this->mesh.vert[e.v2].P();
Point3f point = p[0] + p[1] - p[2];
int vn = this->mesh.vert.size();
for(int i = 0; i < this->mesh.vert.size(); i++)
{
if((this->mesh.vert[i].P() - point).Norm() < 0.1)
{
vn = i;
//find the border
assert(this->mesh.vert[i].IsB());
for(std::list<FrontEdge>::iterator k = this->front.begin(); k != this->front.end(); k++)
if((*k).v0 == i)
{
touch.first = AdvancingFront<MESH>::FRONT;
touch.second = k;
}
for(std::list<FrontEdge>::iterator k = this->deads.begin(); k != this->deads.end(); k++)
if((*k).v0 == i)
if((*k).v0 == i)
{
touch.first = AdvancingFront<MESH>::FRONT;
touch.second = k;
}
break;
}
}
if(vn == this->mesh.vert.size()) {
VertexType v;
v.P() = point;
v.ClearFlags();
AddVertex(v);
}
return vn;
}
};
}//namespace tri
}//namespace vcg
#endif

418
vcg/complex/trimesh/create/ball_pivoting.h
View File

@ -1,418 +0,0 @@
#ifndef BALL_PIVOTING_H
#define BALL_PIVOTING_H
#include "advancing_front.h"
#include <vcg/space/index/grid_static_ptr.h>
#include <vcg/complex/trimesh/closest.h>
/* Ball pivoting algorithm:
1) the vertices used in the new mesh are marked as visited
2) the border vertices of the new mesh are marked as border
3) the vector nb is used to keep track of the number of borders a vertex belongs to
4) usedBit flag is used to select the points in the mesh already processed
*/
namespace vcg {
namespace tri {
template <class MESH> class BallPivoting: public AdvancingFront<MESH> {
public:
typedef typename MESH::VertexType VertexType;
typedef typename MESH::FaceType FaceType;
typedef typename MESH::ScalarType ScalarType;
typedef typename MESH::VertexIterator VertexIterator;
typedef typename MESH::VertexType::CoordType Point3x;
typedef GridStaticPtr<typename MESH::VertexType, typename MESH::ScalarType > StaticGrid;
float radius; //radius of the ball
float min_edge; //min lenght of an edge
float max_edge; //min lenght of an edge
float max_angle; //max angle between 2 faces (cos(angle) actually)
public:
ScalarType radi() { return radius; }
// if radius ==0 an autoguess for the ball pivoting radius is attempted
// otherwise the passed value (in absolute mesh units) is used.
BallPivoting(MESH &_mesh, float _radius = 0,
float minr = 0.2, float angle = M_PI/2):
AdvancingFront<MESH>(_mesh), radius(_radius),
min_edge(minr), max_edge(1.8), max_angle(cos(angle)),
last_seed(-1) {
//compute bbox
baricenter = Point3x(0, 0, 0);
UpdateBounding<MESH>::Box(_mesh);
for(VertexIterator vi=this->mesh.vert.begin();vi!=this->mesh.vert.end();++vi)
if( !(*vi).IsD() ) baricenter += (*vi).P();
baricenter /= this->mesh.vn;
assert(this->mesh.vn > 3);
if(radius == 0) // radius ==0 means that an auto guess should be attempted.
radius = sqrt((this->mesh.bbox.Diag()*this->mesh.bbox.Diag())/this->mesh.vn);
min_edge *= radius;
max_edge *= radius;
//enlarging the bbox for out-of-box queries
Box3<ScalarType> BPbbox=this->mesh.bbox;
BPbbox.Offset(4*radius);
grid.Set(this->mesh.vert.begin(), this->mesh.vert.end(), BPbbox);
//mark visited points
std::vector<VertexType *> targets;
std::vector<Point3x> points;
std::vector<ScalarType> dists;
usedBit = VertexType::NewBitFlag();
for(int i = 0; i < (int)this->mesh.vert.size(); i++)
this->mesh.vert[i].ClearUserBit(usedBit);
UpdateFlags<MESH>::VertexClearV(this->mesh);
for(int i = 0; i < (int)this->mesh.face.size(); i++) {
FaceType &f = this->mesh.face[i];
if(f.IsD()) continue;
for(int k = 0; k < 3; k++) {
f.V(k)->SetV();
int n = tri::GetInSphereVertex(this->mesh, grid, f.V(k)->P(), min_edge, targets, dists, points);
for(int t = 0; t < n; t++) {
targets[t]->SetUserBit(usedBit);
assert(targets[t]->IsUserBit(usedBit));
}
assert(f.V(k)->IsUserBit(usedBit));
}
}
}
~BallPivoting() {
VertexType::DeleteBitFlag(usedBit);
}
bool Seed(int &v0, int &v1, int &v2) {
//get a sphere of neighbours
std::vector<VertexType *> targets;
std::vector<Point3x> points;
std::vector<ScalarType> dists;
while(++last_seed < (int)(this->mesh.vert.size())) {
VertexType &seed = this->mesh.vert[last_seed];
if(seed.IsD() || seed.IsUserBit(usedBit)) continue;
seed.SetUserBit(usedBit);
int n = tri::GetInSphereVertex(this->mesh, grid, seed.P(), 2*radius, targets, dists, points);
if(n < 3) {
continue;
}
bool success = true;
//find the closest visited or boundary
for(int i = 0; i < n; i++) {
VertexType &v = *(targets[i]);
if(v.IsV()) {
success = false;
break;
}
}
if(!success) continue;
VertexType *vv0, *vv1, *vv2;
success = false;
//find a triplet that does not contains any other point
Point3x center;
for(int i = 0; i < n; i++) {
vv0 = targets[i];
if(vv0->IsD()) continue;
Point3x &p0 = vv0->P();
for(int k = i+1; k < n; k++) {
vv1 = targets[k];
if(vv1->IsD()) continue;
Point3x &p1 = vv1->P();
float d2 = (p1 - p0).Norm();
if(d2 < min_edge || d2 > max_edge) continue;
for(int j = k+1; j < n; j++) {
vv2 = targets[j];
if(vv2->IsD()) continue;
Point3x &p2 = vv2->P();
float d1 = (p2 - p0).Norm();
if(d1 < min_edge || d1 > max_edge) continue;
float d0 = (p2 - p1).Norm();
if(d0 < min_edge || d0 > max_edge) continue;
Point3x normal = (p1 - p0)^(p2 - p0);
if(normal.dot(p0 - baricenter) < 0) continue;
/* if(use_normals) {
if(normal * vv0->N() < 0) continue;
if(normal * vv1->N() < 0) continue;
if(normal * vv2->N() < 0) continue;
}*/
if(!FindSphere(p0, p1, p2, center)) {
continue;
}
//check no other point inside
int t;
for(t = 0; t < n; t++) {
if((center - targets[t]->P()).Norm() <= radius)
break;
}
if(t < n) {
continue;
}
//check on the other side there is not a surface
Point3x opposite = center + normal*(((center - p0).dot(normal))*2/normal.SquaredNorm());
for(t = 0; t < n; t++) {
VertexType &v = *(targets[t]);
if((v.IsV()) && (opposite - v.P()).Norm() <= radius)
break;
}
if(t < n) {
continue;
}
success = true;
i = k = j = n;
}
}
}
if(!success) { //see bad luck above
continue;
}
Mark(vv0);
Mark(vv1);
Mark(vv2);
v0 = vv0 - &*this->mesh.vert.begin();
v1 = vv1 - &*this->mesh.vert.begin();
v2 = vv2 - &*this->mesh.vert.begin();
return true;
}
return false;
}
//select a new vertex, mark as Visited and mark as usedBit all neighbours (less than min_edge)
int Place(FrontEdge &edge,typename AdvancingFront<MESH>::ResultIterator &touch) {
Point3x v0 = this->mesh.vert[edge.v0].P();
Point3x v1 = this->mesh.vert[edge.v1].P();
Point3x v2 = this->mesh.vert[edge.v2].P();
/* TODO why using the face normals everything goes wrong? should be
exactly the same................................................
Point3x &normal = mesh.face[edge.face].N(); ?
*/
Point3x normal = ((v1 - v0)^(v2 - v0)).Normalize();
Point3x middle = (v0 + v1)/2;
Point3x center;
if(!FindSphere(v0, v1, v2, center)) {
// assert(0);
return -1;
}
Point3x start_pivot = center - middle;
Point3x axis = (v1 - v0);
ScalarType axis_len = axis.SquaredNorm();
if(axis_len > 4*radius*radius) {
return -1;
}
axis.Normalize();
// r is the radius of the thorus of all possible spheres passing throug v0 and v1
ScalarType r = sqrt(radius*radius - axis_len/4);
std::vector<VertexType *> targets;
std::vector<ScalarType> dists;
std::vector<Point3x> points;
tri::GetInSphereVertex(this->mesh, grid, middle, r + radius, targets, dists, points);
if(targets.size() == 0) {
return -1; //this really would be strange but one never knows.
}
VertexType *candidate = NULL;
ScalarType min_angle = M_PI;
for(int i = 0; i < static_cast<int>(targets.size()); i++) {
VertexType *v = targets[i];
int id = v - &*this->mesh.vert.begin();
if(v->IsD()) continue;
// this should always be true IsB => IsV , IsV => IsU
if(v->IsB()) assert(v->IsV());
if(v->IsV()) assert(v->IsUserBit(usedBit));
if(v->IsUserBit(usedBit) && !(v->IsB())) continue;
if(id == edge.v0 || id == edge.v1 || id == edge.v2) continue;
Point3x p = this->mesh.vert[id].P();
/* Find the sphere through v0, p, v1 (store center on end_pivot */
if(!FindSphere(v0, p, v1, center)) {
continue;
}
/* Angle between old center and new center */
ScalarType alpha = Angle(start_pivot, center - middle, axis);
/* adding a small bias to already chosen vertices.
doesn't solve numerical problems, but helps. */
// if(this->mesh.vert[id].IsB()) alpha -= 0.001;
/* Sometimes alpha might be little less then M_PI while it should be 0,
by numerical errors: happens for example pivoting
on the diagonal of a square. */
/* if(alpha > 2*M_PI - 0.8) {
// Angle between old center and new *point*
//TODO is this really overshooting? shouldbe enough to alpha -= 2*M_PI
Point3x proj = p - axis * (axis * p - axis * middle);
ScalarType beta = angle(start_pivot, proj - middle, axis);
if(alpha > beta) alpha -= 2*M_PI;
} */
if(candidate == NULL || alpha < min_angle) {
candidate = v;
min_angle = alpha;
}
}
if(min_angle >= M_PI - 0.1) {
return -1;
}
if(candidate == NULL) {
return -1;
}
if(!candidate->IsB()) {
assert((candidate->P() - v0).Norm() > min_edge);
assert((candidate->P() - v1).Norm() > min_edge);
}
int id = candidate - &*this->mesh.vert.begin();
assert(id != edge.v0 && id != edge.v1);
Point3x newnormal = ((candidate->P() - v0)^(v1 - v0)).Normalize();
if(normal.dot(newnormal) < max_angle || this->nb[id] >= 2) {
return -1;
}
//test if id is in some border (to return touch
for(std::list<FrontEdge>::iterator k = this->front.begin(); k != this->front.end(); k++)
{
if((*k).v0 == id)
{
touch.first = AdvancingFront<MESH>::FRONT;
touch.second = k;
}
}
for(std::list<FrontEdge>::iterator k = this->deads.begin(); k != this->deads.end(); k++)
{
if((*k).v0 == id)
{
touch.first = AdvancingFront<MESH>::DEADS;
touch.second = k;
}
}
//mark vertices close to candidate
Mark(candidate);
return id;
}
private:
int last_seed; //used for new seeds when front is empty
int usedBit; //use to detect if a vertex has been already processed.
Point3x baricenter;//used for the first seed.
StaticGrid grid; //lookup grid for points
/* returns the sphere touching p0, p1, p2 of radius r such that
the normal of the face points toward the center of the sphere */
bool FindSphere(Point3x &p0, Point3x &p1, Point3x &p2, Point3x &center) {
//we want p0 to be always the smallest one.
Point3x p[3];
if(p0 < p1 && p0 < p2) {
p[0] = p0;
p[1] = p1;
p[2] = p2;
} else if(p1 < p0 && p1 < p2) {
p[0] = p1;
p[1] = p2;
p[2] = p0;
} else {
p[0] = p2;
p[1] = p0;
p[2] = p1;
}
Point3x q1 = p[1] - p[0];
Point3x q2 = p[2] - p[0];
Point3x up = q1^q2;
ScalarType uplen = up.Norm();
//the three points are aligned
if(uplen < 0.001*q1.Norm()*q2.Norm()) {
return false;
}
up /= uplen;
ScalarType a11 = q1.dot(q1);
ScalarType a12 = q1.dot(q2);
ScalarType a22 = q2.dot(q2);
ScalarType m = 4*(a11*a22 - a12*a12);
ScalarType l1 = 2*(a11*a22 - a22*a12)/m;
ScalarType l2 = 2*(a11*a22 - a12*a11)/m;
center = q1*l1 + q2*l2;
ScalarType circle_r = center.Norm();
if(circle_r > radius) {
return false; //need too big a sphere
}
ScalarType height = sqrt(radius*radius - circle_r*circle_r);
center += p[0] + up*height;
return true;
}
/* compute angle from p to q, using axis for orientation */
ScalarType Angle(Point3x p, Point3x q, Point3x &axis) {
p.Normalize();
q.Normalize();
Point3x vec = p^q;
ScalarType angle = acos(p.dot(q));
if(vec.dot(axis) < 0) angle = -angle;
if(angle < 0) angle += 2*M_PI;
return angle;
}
void Mark(VertexType *v) {
std::vector<VertexType *> targets;
std::vector<Point3x> points;
std::vector<ScalarType> dists;
int n = tri::GetInSphereVertex(this->mesh, grid, v->P(), min_edge, targets, dists, points);
for(int t = 0; t < n; t++)
targets[t]->SetUserBit(usedBit);
v->SetV();
}
};
} //namespace tri
} //namespace vcg
#endif

928
vcg/complex/trimesh/create/emc_lookup_table.h
View File

@ -1,928 +0,0 @@
/*===========================================================================*\
* *
* IsoEx *
* Copyright (C) 2002 by Computer Graphics Group, RWTH Aachen *
* www.rwth-graphics.de *
* *
*---------------------------------------------------------------------------*
* *
* License *
* *
* This library is free software; you can redistribute it and/or modify it *
* under the terms of the GNU Library General Public License as published *
* by the Free Software Foundation, version 2. *
* *
* This library is distributed in the hope that it will be useful, but *
* WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
* Library General Public License for more details. *
* *
* You should have received a copy of the GNU Library General Public *
* License along with this library; if not, write to the Free Software *
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *
* *
\*===========================================================================*/
//== TABLES ==================================================================
#ifndef __VCG_EMC_LOOK_UP_TABLE
#define __VCG_EMC_LOOK_UP_TABLE
namespace vcg
{
namespace tri
{
class EMCLookUpTable
{
public:
static const int EdgeTable(unsigned char cubetype)
{
static const int edgeTable[256]=
{
0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0
};
return edgeTable[cubetype];
}; // end of EdgeTable
//-----------------------------------------------------------------------------
static int* TriTable(unsigned char cubetype, int u)
{
static int triTable[256][2][17] =
{{{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{ 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 1, 9, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 2, 10, 9, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 10, 9, 8, 3, 2 , -1,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 0, 8, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 10 */
{{1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 1, 9, 0, 2, 3,11, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 9, 8, 11, 2, 1,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 3, 11,10, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 8, 11, 10, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 11,10, 9, 0, 3, -1,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 15 */
{{9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 8, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 4, 7, 3, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1}},
{{4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 7, 3, 1, 9, 4, -1,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 20 */
{{1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 1, 2,10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 3, 0, 4, 7, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1}},
{{9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 2,10, 9, 0, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}},
{{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1, -1},
{1, 6, 7, 3, 2,10, 9, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 25 */
{{11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 5, 2, 0, 4, 7,11,-1,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 3, 3, 3, 9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1}},
{{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1, -1},
{2, 4, 4, 2, 1, 9, 11, 11,9,4,7, -1, -1, -1, -1, -1 ,-1}},
{{3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 3, 11,10, 1, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1}},
{{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1, -1},
{1, 6, 1, 0, 4, 7,11,10, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 30 */
{{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1, -1},
{2, 3, 5, 4, 7, 8, 0, 3, 11, 10, 9, -1, -1, -1, -1, -1, -1}},
{{4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 4, 7,11,10, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 0, 1, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 35 */
{{8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 3, 1, 5, 4, 8,-1,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 1, 2,10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 3, 3, 3, 3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1}},
{{5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 4, 0, 2,10, 5,-1,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1, -1},
{2, 4, 4, 2, 10, 5, 3, 4, 8, 3, 5, -1, -1, -1, -1, -1, -1}},
/* 40 */
{{9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 0, 8, 11, 2, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}},
{{0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 0, 1, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}},
{{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1, -1},
{1, 6, 2, 1, 5, 4, 8,11, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1},
{ 2, 4, 3, 3,11,10, 1, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1}},
/* 45 */
{{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1, -1},
{2, 3, 5, 4, 9, 5, 1, 0, 8,11, 10, -1, -1, -1, -1, -1, -1}},
{{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1, -1},
{1, 6, 5, 4, 0, 3,11, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 5, 4, 8, 11, 10,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 7, 8, 9, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 5, 7, 3, 0, 9,-1,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 50 */
{{0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 1, 5, 7, 8, 0,-1,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 3, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 7, 8, 9, 5,10, 1, 2, -1, -1, -1, -1, -1, -1, -1}},
{{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1, -1},
{ 2, 3, 5,10, 1, 2, 0, 9, 5, 7, 3,-1, -1, -1, -1, -1, -1}},
{{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1, -1},
{1, 6, 2,10, 5, 7, 8, 0,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 55 */
{{2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 2,10, 5, 7, 3,-1,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 7, 8, 9, 5, 3,11, 2, -1, -1, -1, -1, -1, -1, -1}},
{{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1, -1},
{1, 6, 2, 0, 9, 5, 7,11, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1, -1},
{2, 3, 5, 2, 3,11, 8, 0, 1, 5, 7, -1, -1, -1, -1, -1, -1}},
{{11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5,11, 2, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 60 */
{{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1, -1},
{2, 4, 4, 3,11, 10, 1, 5, 7, 8, 9, -1, -1, -1, -1, -1, -1}},
{{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1, -1},
{1, 7, 5, 7, 11,10, 1, 0, 9, -1, -1, -1, -1, -1, -1, -1, -1}},
{{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1, -1},
{1, 7, 11,10,5, 7, 8, 0,3, -1, -1, -1, -1, -1, -1, -1, -1}},
{{11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 4, 5, 7, 11,10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 3,10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 65 */
{{0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 1, 9, 8, 3, 5,10, 6, -1, -1, -1, -1, -1, -1, -1}},
{{1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 1, 2, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 1, 2, 6, 5, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1}},
/* 70 */
{{9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 0, 2, 6, 5, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1, -1},
{1, 6, 2, 6, 5, 9, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 2, 3,11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1}},
{{11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1},
{ 2, 4, 3, 0, 8, 11, 2, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}},
{{0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 3, 3, 3, 0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1}},
/* 75 */
{{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1, -1},
{2, 3, 5, 5,10, 6, 2, 1, 9, 8,11, -1, -1, -1, -1, -1, -1}},
{{6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 5, 1, 3, 11,6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1, -1},
{1, 6, 5, 1, 0, 8,11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1, -1},
{2, 4, 4, 3, 11, 6, 0, 5, 9, 0, 6, -1, -1, -1, -1}},
{{6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 6, 5, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 80 */
{{5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 5,10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1}},
{{4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 4, 7, 3, 0, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1}},
{{1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 3, 3, 3, 1, 9, 0, 5,10, 6, 8, 4, 7, -1, -1, -1, -1}},
{{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1, -1},
{ 2, 3, 5,10, 6, 5, 9, 4, 7, 3, 1,-1, -1, -1, -1, -1, -1}},
{{6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 1, 2, 6, 5, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1}},
/* 85 */
{{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1, -1},
{2, 4, 4, 2, 6, 5, 1, 3, 0, 4, 7, -1, -1, -1, -1, -1, -1}},
{{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1, -1},
{2, 3, 5, 8, 4, 7, 5, 9, 0, 2, 6, -1, -1, -1, -1, -1, -1}},
{{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1, -1},
{1, 7, 7, 3, 2, 6, 5, 9, 4,-1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 3, 3, 3, 3,11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1}},
{{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1, -1},
{2, 3, 5, 5,10, 6, 7,11, 2, 0, 4, -1, -1, -1, -1, -1, -1}},
/* 90 */
{{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1},
{4, 3, 3, 3, 3, 0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6}},
{{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1, -1},
{3, 4, 4, 3, 2, 1, 9,11, 4, 7, 11, 9, 5, 10, 6, -1, -1}},
{{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1, -1},
{2, 3, 5, 8, 4, 7, 11, 6, 5, 1, 3, -1, -1, -1, -1, -1, -1}},
{{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1, -1},
{1, 7, 5, 1, 0, 4, 7,11, 6, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1, -1},
{3, 4, 4, 3, 0, 6, 5, 9, 3, 11, 6, 0, 8, 4, 7, -1, -1}},
/* 95 */
{{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1, -1},
{2, 4, 4, 9, 4, 7, 11, 6, 5, 9, 11,-1, -1, -1, -1, -1, -1}},
{{10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 4, 4, 9, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 4, 9,10, 6, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1}},
{{10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 5, 6, 4, 0, 1, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1, -1},
{1, 6, 1,10, 6, 4, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 100 */
{{1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 2, 6, 4, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1, -1},
{2, 3, 5, 3, 0, 8, 9, 1, 2, 6, 4, -1, -1, -1, -1, -1, -1}},
{{0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 2, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 8, 3, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1, -1},
{ 2, 4, 3, 10, 6, 4, 9,11, 2, 3, -1, -1, -1, -1, -1, -1, -1}},
/* 105 */
{{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1, -1},
{2, 4, 4, 2, 11, 8, 0, 10, 6, 4, 9, -1, -1, -1, -1, -1, -1}},
{{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1, -1},
{2, 3, 5, 3,11, 2, 1, 10,6, 4, 0, -1, -1, -1, -1, -1, -1}},
{{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1, -1},
{1, 7, 6, 4, 8,11, 2, 1,10, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1, -1},
{1, 6, 3,11, 6, 4, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1, -1},
{1, 7, 8,11, 6, 4, 9, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 110 */
{{3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 3,11, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 8, 11, 6, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 8, 9,10, 6, 7,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1, -1},
{1, 6, 0, 9, 10, 6, 7, 3, -1,-1,-1, -1, -1, -1, -1, -1, -1}},
{{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1, -1},
{ 2, 4, 4, 8, 0, 1, 7, 10, 6, 7, 1,-1, -1, -1, -1, -1, -1}},
/* 115 */
{{10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 5, 10, 6, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1, -1},
{1, 6, 1, 2, 6, 7, 8, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1, -1},
{1, 7, 2, 6, 7, 3, 0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 7, 8, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 7, 3, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 120 */
{{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1, -1},
{2, 3, 5, 2, 3,11, 6, 7, 8, 9,10, -1, -1, -1, -1, -1, -1}},
{{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1, -1},
{1, 7, 2, 0, 9,10,6, 7, 11, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1, -1},
{3, 4, 4, 3, 8, 0, 1, 7, 10, 6, 7, 1, 11, 2, 3, -1, -1}},
{{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1, -1},
{ 2, 4, 4, 11, 2, 1,7, 1, 10, 6, 7,-1, -1, -1, -1, -1, -1}},
{{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1, -1},
{1, 7, 8, 9, 1, 3, 11, 6, 7,-1, -1, -1, -1, -1, -1, -1, -1}},
/* 125 */
{{0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1}},
{{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1, -1},
{2, 4, 4, 0, 3,11, 6, 7, 8, 0, 6, -1, -1, -1, -1, -1, -1}},
{{7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 130 */
{{0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1}},
{{8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 1, 9, 8, 3,11, 7, 6, -1, -1, -1, -1, -1, -1, -1}},
{{10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{ 2, 3, 3,10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 3, 3, 3, 1, 2,10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1}},
{{2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 2, 10, 9, 0, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}},
/* 135 */
{{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1, -1},
{2, 3, 5, 6, 11, 7, 3, 2,10, 9, 8, -1, -1, -1, -1, -1, -1}},
{{7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 2, 3, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 6, 2, 0, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 2, 3, 7, 6, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1}},
{{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1, -1},
{1, 6, 6, 2, 1, 9, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 140 */
{{10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 5, 1, 3, 7, 6,10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1, -1},
{ 2, 4, 4, 10, 1, 7, 6, 8, 7, 1, 0,-1, -1, -1, -1, -1, -1}},
{{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1, -1},
{1, 6,10, 9, 0, 3, 7, 6,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 7, 6, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 6, 11, 8, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 145 */
{{3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 0, 4, 6,11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 6,11, 8, 4, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1, -1},
{1, 6, 6,11, 3, 1, 9, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 6, 11, 8, 4, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1, -1},
{2, 3, 5, 1, 2, 10,11, 3,0,4, 6, -1, -1, -1, -1, -1, -1}},
/* 150 */
{{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1, -1},
{2, 4, 4, 4, 6, 11, 8, 2,10, 9, 0, -1, -1, -1, -1, -1, -1}},
{{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1, -1},
{1, 7, 10,9, 4, 6, 11, 3, 2, -1, -1, -1, -1, -1, -1, -1, -1}},
{{8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 4, 6, 2, 3, 8,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 4, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1, -1},
{2, 3, 5, 1, 9, 0, 3, 8, 4, 6, 2, -1, -1, -1, -1, -1, -1}},
/* 155 */
{{1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 1, 9, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1, -1},
{1, 6, 1, 3, 8, 4, 6,10, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 5,10, 1,0,4,6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1, -1},
{1, 7, 4, 6, 10, 9, 0,3, 8, -1, -1, -1, -1, -1, -1, -1, -1}},
{{10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 4, 4, 6, 10, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 160 */
{{4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 3, 3, 3, 0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1}},
{{5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 0, 1, 5, 4, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}},
{{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1, -1},
{ 2, 3, 5,11, 7, 6, 4, 8, 3, 1, 5,-1, -1, -1, -1, -1, -1}},
{{9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 3, 3, 3, 9, 5, 4,10, 1, 2, 7, 6, 11, -1, -1, -1, -1}},
/* 165 */
{{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1, -1},
{4, 3, 3, 3, 3, 6,11, 7, 1, 2,10, 0, 8, 3, 4, 9, 5}},
{{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1, -1},
{2, 3, 5, 7, 6, 11, 10, 5, 4, 0, 2,-1, -1, -1, -1, -1, -1}},
{{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1, -1},
{3, 4, 4, 3, 5, 3, 2,10, 4, 8, 3, 5, 6, 11, 7, 6, -1}},
{{7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 4, 3, 2, 3, 7, 6, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1}},
{{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1, -1},
{2, 3, 5, 9, 5, 4, 8, 7, 6, 2, 0, -1, -1, -1, -1, -1, -1}},
/* 170 */
{{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1, -1},
{2, 4, 4, 3, 7, 6, 2, 0, 1, 5, 4, -1, -1, -1, -1, -1, -1}},
{{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1, -1},
{1, 7, 6, 2, 1, 5, 4, 8, 7,-1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1, -1},
{2, 3, 5, 9, 5, 4, 6,10, 1, 3, 7,-1, -1, -1, -1, -1, -1}},
{{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1, -1},
{3, 4, 4, 3, 0, 8, 7, 1, 6, 10, 1, 7, 9, 5, 4, -1, -1}},
{{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1, -1},
{1, 7, 4, 0, 3, 7, 6, 10, 5, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 175 */
{{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1, -1},
{2, 4, 4, 4, 8, 10, 5, 7, 6,10, 8, -1, -1, -1, -1, -1, -1}},
{{6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5,11, 8, 9, 5, 6,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1, -1},
{2, 4, 4, 0, 9, 5, 6, 6,11, 3, 0, -1, -1, -1, -1, -1, -1}},
{{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1, -1},
{1, 6, 0, 1, 5, 6,11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 6,11, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/*180 */
{{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1, -1},
{2, 3, 5, 1, 2, 10, 5, 6,11, 8, 9, -1, -1, -1, -1, -1, -1}},
{{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1, -1},
{3, 4, 4, 3, 11, 3,0, 6, 9, 5, 6, 0, 2, 10, 1, 2, 10}},
{{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1, -1},
{ 1, 7,11, 8, 0, 2,10, 5, 6,-1, -1, -1, -1, -1, -1, -1, -1}},
{{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1, -1},
{2, 4, 4, 6,11, 3, 5, 10, 5, 3, 2, -1, -1, -1, -1, -1, -1}},
{{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1, -1},
{1, 6, 2, 3, 8, 9, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 185 */
{{9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 9, 5, 6, 2, 0,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1, -1},
{1, 7, 1, 5, 6, 2, 3, 8, 0, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 1, 5, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1, -1},
{1, 7, 1, 3, 8, 9, 5, 6,10, -1, -1, -1, -1, -1, -1, -1, -1}},
{{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1, -1},
{ 2, 4, 4, 5, 6, 0, 9, 10, 1, 0, 6, -1, -1, -1, -1, -1, -1}},
/* 190 */
{{0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1}},
{{10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 3,10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 4, 5,10, 11, 7,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1, -1},
{ 2, 4, 3, 5,10,11, 7, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1}},
{{5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1, -1},
{ 2, 4, 3, 5, 10, 11, 7, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1}},
/* 195 */
{{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1, -1},
{ 2, 4, 4, 10, 11, 7, 5, 1, 9, 8, 3, -1, -1, -1, -1, -1, -1}},
{{11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 5, 7, 5, 1, 2,11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1, -1},
{2, 3, 5, 0, 8, 3, 2,11, 7, 5,1, -1, -1, -1, -1, -1, -1}},
{{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1, -1},
{1, 6, 2,11, 7, 5, 9, 0,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1, -1},
{1, 7, 7, 5, 9, 8, 3, 2,11,-1, -1, -1, -1, -1, -1, -1, -1}},
/* 200 */
{{2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 3, 7, 5,10, 2,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1, -1},
{1, 6, 5,10, 2, 0, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1, -1},
{2, 3, 5, 9, 0, 1, 10, 2, 3, 7, 5, -1, -1, -1, -1, -1, -1}},
{{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1, -1},
{1, 7, 9, 8, 7, 5,10, 2, 1,-1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 3, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 205 */
{{0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 0, 8, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 9, 0, 3, 7, 5,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 7, 5, 9, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 10,11, 8, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1, -1},
{1, 6, 0, 4, 5,10,11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 210 */
{{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1, -1},
{2, 3, 5, 0, 1, 9, 4, 5, 10, 11, 8, -1, -1, -1, -1, -1, -1}},
{{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1, -1},
{ 1, 7,10, 11, 3, 1, 9,4, 5,-1, -1, -1, -1, -1, -1, -1}},
{{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1, -1},
{1, 6, 2,11, 8, 4, 5, 1,-1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1, -1},
{1, 7, 0, 4, 5, 1, 2, 11, 3,-1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1, -1},
{1, 7, 0, 2,11, 8, 4, 5, 9, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 215 */
{{9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1}},
{{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1, -1},
{2, 4, 4, 2, 3, 5, 10, 4, 5, 3, 8,-1, -1, -1, -1, -1, -1}},
{{5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 5,10, 2, 0, 4,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1, -1},
{3, 4, 4, 3, 3, 5, 10, 2, 8, 4, 5, 3, 0, 1, 9, -1, -1}},
{{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1, -1},
{1, 6,10, 2, 1, 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 220 */
{{8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 8, 4, 5, 1, 3,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 0, 4, 5, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1, -1},
{2, 4, 4, 0, 3, 5, 9, 8, 4, 5, 3, -1, -1, -1, -1, -1, -1}},
{{9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 9,10, 11, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 225 */
{{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1, -1},
{2, 3, 5, 0, 8, 3, 7, 4, 9, 10, 11, -1, -1, -1, -1, -1, -1}},
{{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1, -1},
{1, 6, 1, 10,11, 7, 4, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1, -1},
{1, 7, 3, 1,10,11, 7, 4, 8, -1, -1, -1, -1, -1, -1, -1, -1}},
{{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1, -1},
{2, 4, 4, 2, 11, 9, 1, 4, 9, 11, 7, -1, -1, -1, -1, -1, -1}},
{{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1, -1},
{3, 4, 4, 3, 1, 2, 11, 9, 7, 4, 9,11, 8, 3, 0, 8, 3}},
/* 230 */
{{11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1, -1},
{ 1, 5, 11, 7, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1, -1},
{ 2, 4, 4, 11, 7, 4, 2, 3, 2, 4, 8,-1, -1, -1, -1, -1, -1}},
{{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1, -1},
{1, 6, 2, 3, 7, 4, 9,10, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1, -1},
{1, 7, 9,10, 2, 0, 8, 7, 4,-1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1, -1},
{1, 7, 3, 7, 4, 0, 1,10, 2, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 235 */
{{1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 3, 1,10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1}},
{{4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 4, 9, 1, 3, 7,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1, -1},
{2, 4, 4, 8, 7, 1, 0, 4, 9, 1, 7, -1, -1, -1, -1, -1, -1}},
{{4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 3, 7, 4, 0,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 240 */
{{9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 8, 9, 10,11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 3, 0, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 0, 1, 10,11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 3, 1,10, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 1, 2, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 245 */
{{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1, -1},
{2, 4, 4, 2,11, 9, 1, 3, 0, 9, 11, -1, -1, -1, -1, -1,-1}},
{{0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 0, 2,11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 2, 3, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 2, 0, 9,10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
/* 250 */
{{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1, -1},
{2, 4, 4, 2, 3, 8, 10, 1, 10, 8, 0, -1, -1, -1, -1, -1, -1}},
{{1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 4, 1, 3, 8, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}},
{{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{ 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}
};
return &triTable[cubetype][u][0];
}; // end of TriTable
//-----------------------------------------------------------------------------
static const int PolyTable(unsigned int cubetype, int u)
{
static const int polyTable[8][16] =
{
{-1},
{-1},
{-1},
{0, 1, 2, -1},
{0, 1, 2, 2, 3, 0, -1},
{0, 1, 2, 0, 2, 4, 4, 2, 3, -1},
{0, 1, 2, 2, 3, 4, 4, 5, 0, 0, 2, 4, -1},
{0, 1, 5, 0, 5, 6, 1, 2, 5, 4, 5, 3, 2, 3, 5, -1}
};
return polyTable[cubetype][u];
}; // end of PolyTable
//=============================================================================
}; //end of class EMCLookUpTable
}; // end of namespace tri
}; // end of namespace vcg
#endif // __VCG_EMC_LOOK_UP_TABLE

462
vcg/complex/trimesh/create/extended_marching_cubes.h
View File

@ -1,462 +0,0 @@
/****************************************************************************
* VCGLib o o *
* Visual and Computer Graphics Library o o *
* _ O _ *
* Copyright(C) 2004 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
/***************************************************************************/
#ifndef __VCG_EXTENDED_MARCHING_CUBES
#define __VCG_EXTENDED_MARCHING_CUBES
#include <float.h>
#include <assert.h>
#include <vector>
#include <vcg/math/base.h>
#include <vcg/math/matrix.h>
#include <vcg/math/lin_algebra.h>
#include <vcg/simplex/face/topology.h>
#include <vcg/complex/trimesh/update/edges.h>
#include <vcg/complex/trimesh/update/normal.h>
#include <vcg/complex/trimesh/update/topology.h>
#include <vcg/complex/trimesh/allocate.h>
#include <vcg/complex/trimesh/clean.h>
#include <vcg/space/point3.h>
#include "emc_lookup_table.h"
namespace vcg
{
namespace tri
{
// Doxygen documentation
/** \addtogroup trimesh */
/*@{*/
/*
* Cube description:
* 3 ________ 2 _____2__
* /| /| / | /|
* / | / | 11/ 3 10/ |
* 7 /_______ / | /__6_|__ / |1
* | | |6 | | | |
* | 0|__|_____|1 | |__|_0|__|
* | / | / 7 8/ 5 /
* | / | / | / | /9
* |/_______|/ |/___4___|/
* 4 5
*/
//! This class implements the Extended Marching Cubes algorithm.
/*!
* 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
* mesh the triangles set approximating the surface that cross that cell.
* @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.
**/
template<class TRIMESH_TYPE, class WALKER_TYPE>
class ExtendedMarchingCubes
{
public:
#if defined(__GNUC__)
typedef unsigned int size_t;
#else
#ifdef _WIN64
typedef unsigned __int64 size_t;
#else
typedef _W64 unsigned int size_t;
#endif
#endif
typedef typename vcg::tri::Allocator< TRIMESH_TYPE > AllocatorType;
typedef typename TRIMESH_TYPE::ScalarType ScalarType;
typedef typename TRIMESH_TYPE::VertexType VertexType;
typedef typename TRIMESH_TYPE::VertexPointer VertexPointer;
typedef typename TRIMESH_TYPE::VertexIterator VertexIterator;
typedef typename TRIMESH_TYPE::FaceType FaceType;
typedef typename TRIMESH_TYPE::FacePointer FacePointer;
typedef typename TRIMESH_TYPE::FaceIterator FaceIterator;
typedef typename TRIMESH_TYPE::CoordType CoordType;
typedef typename TRIMESH_TYPE::CoordType* CoordPointer;
struct LightEdge
{
LightEdge(size_t _face, size_t _edge):face(_face), edge(_edge) { }
size_t face, edge;
};
/*!
* Constructor
* \param mesh The mesh that will be constructed
* \param volume The volume describing the field
* \param walker The class implementing the traversal policy
* \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)
{
_mesh = &mesh;
_walker = &walker;
_featureAngle = vcg::math::ToRad(angle);
_initialized = _finalized = false;
};
/*!
* Execute the initialiazation.
* This method must be executed before the first call to <CODE>ApplyEMC</CODE>
*/
void Initialize()
{
assert(!_initialized && !_finalized);
_featureFlag = VertexType::NewBitFlag();
_initialized = true;
};
/*!
*
* This method must be executed after the last call to <CODE>ApplyEMC</CODE>
*/
void Finalize()
{
assert(_initialized && !_finalized);
FlipEdges();
VertexIterator v_iter = _mesh->vert.begin();
VertexIterator v_end = _mesh->vert.end();
for ( ; v_iter!=v_end; v_iter++)
v_iter->ClearUserBit( _featureFlag );
VertexType::DeleteBitFlag( _featureFlag );
_featureFlag = 0;
_mesh = NULL;
_walker = NULL;
_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>.
* 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 max the second point
*/
void ProcessCell(const vcg::Point3i &min, const vcg::Point3i &max)
{
assert(_initialized && !_finalized);
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[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[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[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[7].X()=min.X(); _corners[7].Y()=max.Y(); _corners[7].Z()=max.Z();
unsigned char cubetype = 0;
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[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[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[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 (cubetype==0 || cubetype==255)
return;
size_t vertices_idx[12];
memset(vertices_idx, -1, 12*sizeof(size_t));
int code = EMCLookUpTable::EdgeTable(cubetype);
VertexPointer vp = NULL;
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 ( 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 ( 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 ( 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 ( 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 (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]; }
int m, n, vertices_num;
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);
std::vector< size_t > vertices_list;
for (m=1; m<=components; m++)
{
// current sheet contains vertices_num vertices
vertices_num = EMCLookUpTable::TriTable(cubetype, 1)[m]; //vertices_num = triTable[cubetype][1][m];
// collect vertices
vertices_list.clear();
for (n=0; n<vertices_num; ++n)
vertices_list.push_back( vertices_idx[ indices[n] ] );
VertexPointer feature = FindFeature( vertices_list );
if (feature != NULL) // i.e. is a valid vertex
{
// feature -> create triangle fan around feature vertex
size_t feature_idx = feature - &_mesh->vert[0];
size_t face_idx = _mesh->face.size();
vertices_list.push_back( vertices_list[0] );
AllocatorType::AddFaces(*_mesh, (int) vertices_num);
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(1) = &_mesh->vert[ vertices_list[j+1] ];
_mesh->face[face_idx].V(2) = &_mesh->vert[ feature_idx ];
}
}
else
{
// 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)
{
size_t face_idx = _mesh->face.size();
AllocatorType::AddFaces(*_mesh, 1);
//_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(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(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) ] ] ];
}
}
indices += vertices_num;
}
}; // end of ApplyEMC
private:
/*!
*/
WALKER_TYPE *_walker;
/*!
*/
TRIMESH_TYPE *_mesh;
/*!
*/
bool _initialized;;
/*!
*/
bool _finalized;
/*!
* The feature detection threshold misuring the sharpness of a feature
*/
ScalarType _featureAngle;
/*!
* The flag used for marking the feature vertices.
*/
int _featureFlag;
/*!
* Array of the 8 corners of the volume cell being processed
*/
vcg::Point3i _corners[8];
/*!
* The field value at the cell corners
*/
ScalarType _field[8];
/*!
* 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
* \return The pointer to the new Vertex if a feature is detected; NULL otherwise.
*/
VertexPointer FindFeature(const std::vector<size_t> &vertices_idx)
{
unsigned int i, j, rank;
size_t vertices_num = (size_t) vertices_idx.size();
CoordType *points = new CoordType[ vertices_num ];
CoordType *normals = new CoordType[ vertices_num ];
Box3<ScalarType> bb;
for (i=0; i<vertices_num; i++)
{
points[i] = _mesh->vert[ vertices_idx[i] ].P();
normals[i].Import(_mesh->vert[ vertices_idx[i] ].N());
bb.Add(points[i]);
}
// move barycenter of points into (0, 0, 0)
CoordType center((ScalarType) 0.0, (ScalarType) 0.0, (ScalarType) 0.0);
for (i=0; i<vertices_num; ++i)
center += points[i];
center /= (ScalarType) vertices_num;
for (i=0; i<vertices_num; ++i)
points[i] -= center;
// normal angle criterion
double c, minC, maxC;
CoordType axis;
for (minC=1.0, i=0; i<vertices_num-1; ++i)
{
for (j=i+1; j<vertices_num; ++j)
{
c = normals[i]*normals[j];
if (c < minC)
{
minC = c;
axis = normals[i] ^ normals[j];
}
}
} //end for (minC=1.0, i=0; i<vertNumber; ++i)
if (minC > cos(_featureAngle))
return NULL; // invalid vertex
// ok, we have a feature: is it edge or corner, i.e. rank 2 or 3 ?
axis.Normalize();
for (minC=1.0, maxC=-1.0, i=0; i<vertices_num; ++i)
{
c = axis * normals[i];
if (c < minC) minC = c;
if (c > maxC) maxC = c;
}
c = std::max< double >(fabs(minC), fabs(maxC));
c = sqrt(1.0-c*c);
rank = (c > cos(_featureAngle) ? 2 : 3);
// setup linear system (find intersection of tangent planes)
vcg::ndim::Matrix<double> A((unsigned int) vertices_num, 3);
double *b = new double[ vertices_num ];
for (i=0; i<vertices_num; ++i)
{
A[i][0] = normals[i][0];
A[i][1] = normals[i][1];
A[i][2] = normals[i][2];
b[i] = (points[i] * normals[i]);
}
// SVD of matrix A
vcg::ndim::Matrix<double> V(3, 3);
double *w = new double[vertices_num];
vcg::SingularValueDecomposition< typename vcg::ndim::Matrix<double> > (A, w, V, LeaveUnsorted, 100);
// rank == 2 -> suppress smallest singular value
if (rank == 2)
{
double smin = DBL_MAX; // the max value, as defined in <float.h>
unsigned int sminid = 0;
unsigned int srank = std::min< unsigned int >(vertices_num, 3u);
for (i=0; i<srank; ++i)
{
if (w[i] < smin)
{
smin = w[i];
sminid = i;
}
}
w[sminid] = 0.0;
}
// SVD backsubstitution -> least squares, least norm solution x
double *x = new double[3];
vcg::SingularValueBacksubstitution< vcg::ndim::Matrix<double> >(A, w, V, x, b);
// transform x to world coords
CoordType point((ScalarType) x[0], (ScalarType) x[1], (ScalarType) x[2]);
point += center;
// 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...
if(!bb.IsIn(point)) point = center;
// insert the feature-point
VertexPointer mean_point = &*AllocatorType::AddVertices( *_mesh, 1);
mean_point->SetUserBit(_featureFlag);
mean_point->P() = point;
mean_point->N().SetZero();
delete []x;
delete []points;
delete []normals;
return mean_point;
} // end of FindFeature
/*!
* 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
* feature samples after the flip.
*/
void FlipEdges()
{
size_t i;
std::vector< LightEdge > edges;
FaceIterator f_iter = _mesh->face.begin();
FaceIterator f_end = _mesh->face.end();
for (i=0; f_iter!=f_end; f_iter++, i++)
{
if (f_iter->V(1) > f_iter->V(0)) edges.push_back( LightEdge(i,0) );
if (f_iter->V(2) > f_iter->V(1)) edges.push_back( LightEdge(i,1) );
if (f_iter->V(0) > f_iter->V(2)) edges.push_back( LightEdge(i,2) );
}
vcg::tri::UpdateTopology< TRIMESH_TYPE >::FaceFace( *_mesh );
// Select all the triangles that has a vertex shared with a non manifold edge.
int nonManifEdge = tri::Clean< TRIMESH_TYPE >::CountNonManifoldEdgeFF(*_mesh,true);
if(nonManifEdge >0)
tri::UpdateSelection< TRIMESH_TYPE >::FaceFromVertexLoose(*_mesh);
//qDebug("Got %i non manif edges",nonManifEdge);
typename std::vector< LightEdge >::iterator e_it = edges.begin();
typename std::vector< LightEdge >::iterator e_end = edges.end();
FacePointer g, f;
int w, z;
for( ; e_it!=e_end; e_it++)
{
f = &_mesh->face[e_it->face];
z = (int) e_it->edge;
// v2------v1 swap the diagonal only if v2 and v3 are feature and v0 and v1 are not.
// | / |
// | / |
// v0------v3
if (!(f->IsS()) && vcg::face::CheckFlipEdge< FaceType >(*f, z))
{
VertexPointer v0, v1, v2, v3;
v0 = f->V(z);
v1 = f->V1(z);
v2 = f->V2(z);
g = f->FFp(z);
w = f->FFi(z);
v3 = g->V2(w);
bool b0, b1, b2, b3;
b0 = !v0->IsUserBit(_featureFlag) ;
b1 = !v1->IsUserBit(_featureFlag) ;
b2 = v2->IsUserBit(_featureFlag) ;
b3 = v3->IsUserBit(_featureFlag) ;
if( b0 && b1 && b2 && b3)
vcg::face::FlipEdge< FaceType >(*f, z);
} // end if (vcg::face::CheckFlipEdge< _Face >(*f, z))
} // end for( ; e_it!=e_end; e_it++)
}; //end of FlipEdges
}; // end of class ExtendedMarchingCubes
// /*! @} */
// end of Doxygen documentation
} // end of namespace tri
}; // end of namespace vcg
#endif // __VCG_EXTENDED_MARCHING_CUBES

730
vcg/complex/trimesh/create/marching_cubes.h
View File

@ -1,730 +0,0 @@
/****************************************************************************
* VCGLib o o *
* Visual and Computer Graphics Library o o *
* _ O _ *
* Copyright(C) 2004 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
/***************************************************************************/
#ifndef __VCG_MARCHING_CUBES
#define __VCG_MARCHING_CUBES
#include <assert.h>
#include <vcg/space/point3.h>
#include <vcg/complex/trimesh/allocate.h>
#include "mc_lookup_table.h"
namespace vcg
{
namespace tri
{
// Doxygen documentation
/** \addtogroup trimesh */
/*@{*/
/*
* Cube description:
* 3 ________ 2 _____2__
* /| /| / | /|
* / | / | 11/ 3 10/ |
* 7 /_______ / | /__6_|__ / |1
* | | |6 | | | | |
* | 0|__|_____|1 | |__|__0__|
* | / | / 7 8/ 5 /
* | / | / | / | /9
* |/_______|/ |/___4___|/
* 4 5
*/
//! This class implements the Marching Cubes algorithm.
/*!
* 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
* mesh the triangles set approximating the surface that cross that cell. The ambiguities
* are resolved using an enhanced topologically controlled lookup table.
* @param TRIMESH_TYPE (Template parameter) the mesh type that will be constructed
* @param WALKER_TYPE (Template parameter) the class that implement the traversal ordering of the volume
**/
template<class TRIMESH_TYPE, class WALKER_TYPE>
class MarchingCubes
{
public:
enum Dimension {X, Y, Z};
#if defined(__GNUC__)
typedef unsigned int size_t;
#else
#ifdef _WIN64
typedef unsigned __int64 size_t;
#else
typedef _W64 unsigned int size_t;
#endif
#endif
typedef typename vcg::tri::Allocator< TRIMESH_TYPE > AllocatorType;
typedef typename TRIMESH_TYPE::ScalarType ScalarType;
typedef typename TRIMESH_TYPE::VertexType VertexType;
typedef typename TRIMESH_TYPE::VertexPointer VertexPointer;
typedef typename TRIMESH_TYPE::VertexIterator VertexIterator;
typedef typename TRIMESH_TYPE::FaceType FaceType;
typedef typename TRIMESH_TYPE::FacePointer FacePointer;
typedef typename TRIMESH_TYPE::FaceIterator FaceIterator;
typedef typename TRIMESH_TYPE::CoordType CoordType;
typedef typename TRIMESH_TYPE::CoordType* CoordPointer;
/*!
* Constructor
* \param mesh the mesh that will be constructed
* \param walker the class implementing the traversal policy
*/
MarchingCubes(TRIMESH_TYPE &mesh, WALKER_TYPE &walker)
{
_mesh = &mesh;
_walker = &walker;
};
/*!
* Execute the initialiazation.
* This method must be executed before the first call to <CODE>ApplyMC</CODE>
*/
void Initialize()
{
_mesh->Clear();
}; // end of Initialize()
/*!
*
* This method must be executed after the last call to <CODE>ApplyMC</CODE>
*/
void Finalize()
{
_mesh = NULL;
_walker = NULL;
}; // end of Finalize()
/*!
* Apply the <I>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.
* \param min the first point
* \param max the second point
*/
void ProcessCell(const vcg::Point3i &min, const vcg::Point3i &max)
{
_case = _subconfig = _config = -1;
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[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[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[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[7].X()=min.X(); _corners[7].Y()=max.Y(); _corners[7].Z()=max.Z();
for (int i=0; i<8; i++)
_field[i] = _walker->V( _corners[i].X(), _corners[i].Y(), _corners[i].Z() );
unsigned char cubetype = 0;
for (int i=0; i<8; i++)
if (_field[i]>0) cubetype += 1<<i;
_case = MCLookUpTable::Cases(cubetype, 0); //_case = cases[cubetype][0];
_config = MCLookUpTable::Cases(cubetype, 1); //_config = cases[cubetype][1];
_subconfig = 0;
VertexPointer v12 = NULL;
switch( _case )
{
case 0 : { break ; }
case 1 : { AddTriangles( MCLookUpTable::Tiling1(_config), 1 ); break; } //case 1 : { AddTriangles( tiling1[_config], 1 ); break; }
case 2 : { AddTriangles( MCLookUpTable::Tiling2(_config), 2 ); break; } //case 2 : { AddTriangles( tiling2[_config], 2 ); break; }
case 3 :
{
//if( TestFace( test3[_config]) ) AddTriangles( tiling3_2[_config], 4 ) ; // 3.2
if( TestFace( MCLookUpTable::Test3(_config)) )
AddTriangles( MCLookUpTable::Tiling3_2(_config), 4 ) ; // 3.2
else
AddTriangles( MCLookUpTable::Tiling3_1(_config), 2 ) ; // 3.1
break ;
}
case 4 :
{
//if( TestInterior( test4[_config]) ) AddTriangles( tiling4_1[_config], 2 ) ; // 4.1.1
if( TestInterior( MCLookUpTable::Test4(_config) ) )
AddTriangles( MCLookUpTable::Tiling4_1(_config), 2 ) ; // 4.1.1
else
AddTriangles( MCLookUpTable::Tiling4_2(_config), 6 ) ; // 4.1.2
break ;
}
case 5 : { AddTriangles( MCLookUpTable::Tiling5(_config), 3 ); break; }
case 6 :
{
//if( TestFace( test6[_config][0]) )
if( TestFace( MCLookUpTable::Test6(_config, 0)) )
AddTriangles( MCLookUpTable::Tiling6_2(_config), 5 ) ; // 6.2
else
{
if( TestInterior( MCLookUpTable::Test6(_config, 1)) )
AddTriangles( MCLookUpTable::Tiling6_1_1(_config), 3 ) ; // 6.1.1
else
AddTriangles( MCLookUpTable::Tiling6_1_2(_config), 7 ) ; // 6.1.2
}
break ;
}
case 7 :
{
//if( TestFace( test7[_config][0] ) ) _subconfig += 1 ;
//if( TestFace( test7[_config][1] ) ) _subconfig += 2 ;
//if( TestFace( test7[_config][2] ) ) _subconfig += 4 ;
if( TestFace( MCLookUpTable::Test7(_config, 0) ) ) _subconfig += 1 ;
if( TestFace( MCLookUpTable::Test7(_config, 1) ) ) _subconfig += 2 ;
if( TestFace( MCLookUpTable::Test7(_config, 2) ) ) _subconfig += 4 ;
switch( _subconfig )
{
case 0 : { AddTriangles( MCLookUpTable::Tiling7_1(_config), 3 ) ; break; }
case 1 : { AddTriangles( MCLookUpTable::Tiling7_2(_config,0), 5 ) ; break; }
case 2 : { AddTriangles( MCLookUpTable::Tiling7_2(_config,1), 5 ) ; break; }
case 3 : { ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling7_3(_config,0), 9, v12 ) ; break ; }
case 4 : { AddTriangles( MCLookUpTable::Tiling7_2(_config, 2), 5 ) ; break ;}
case 5 : { ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling7_3(_config,1), 9, v12 ) ; break ; }
case 6 : { ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling7_3(_config,2), 9, v12 ) ; break ; }
case 7 :
{
if( TestInterior( MCLookUpTable::Test7(_config, 3) ) )
AddTriangles( MCLookUpTable::Tiling7_4_2(_config), 9 ) ;
else
AddTriangles( MCLookUpTable::Tiling7_4_1(_config), 5 ) ;
break ;
}
};
break ;
} // end of case 7
case 8 : { AddTriangles( MCLookUpTable::Tiling8(_config), 2 ) ; break ;}
case 9 : { AddTriangles( MCLookUpTable::Tiling9(_config), 4 ) ; break ;}
case 10 :
{
if( TestFace( MCLookUpTable::Test10(_config, 0)) ) //if( TestFace( test10[_config][0]) )
{
if( TestFace( MCLookUpTable::Test10(_config,1) ) )
AddTriangles( MCLookUpTable::Tiling10_1_1_(_config), 4 ) ; // 10.1.1
else
{
ComputeCVertex(v12);
AddTriangles( MCLookUpTable::Tiling10_2(_config), 8, v12 ) ; // 10.2
}
}
else
{
if( TestFace( MCLookUpTable::Test10(_config, 1) ) )
{
ComputeCVertex(v12) ;
AddTriangles( MCLookUpTable::Tiling10_2_(_config), 8, v12 ) ; // 10.2
}
else
{
if( TestInterior( MCLookUpTable::Test10(_config, 2) ) )
AddTriangles( MCLookUpTable::Tiling10_1_1(_config), 4 ) ; // 10.1.1
else
AddTriangles( MCLookUpTable::Tiling10_1_2(_config), 8 ) ; // 10.1.2
}
}
break ;
} // end of case 10
case 11 : { AddTriangles( MCLookUpTable::Tiling11(_config), 4 ) ; break ; }
case 12 :
{
if( TestFace( MCLookUpTable::Test12(_config, 0) ) ) //if( TestFace( test12[_config][0]) )
{
if( TestFace( MCLookUpTable::Test12(_config, 1) ) )
AddTriangles( MCLookUpTable::Tiling12_1_1_(_config), 4 ) ; // 12.1.1
else
{
ComputeCVertex(v12) ;
AddTriangles( MCLookUpTable::Tiling12_2(_config), 8, v12 ) ; // 12.2
}
}
else
{
if( TestFace( MCLookUpTable::Test12(_config, 1) ) )
{
ComputeCVertex(v12) ;
AddTriangles( MCLookUpTable::Tiling12_2_(_config), 8, v12 ) ; // 12.2
}
else
{
if( TestInterior( MCLookUpTable::Test12(_config, 2) ) )
AddTriangles( MCLookUpTable::Tiling12_1_1(_config), 4 ) ; // 12.1.1
else
AddTriangles( MCLookUpTable::Tiling12_1_2(_config), 8 ) ; // 12.1.2
}
}
break ;
} // end of case 12
case 13 :
{
//if( TestFace( test13[_config][0] ) ) _subconfig += 1 ;
//if( TestFace( test13[_config][1] ) ) _subconfig += 2 ;
//if( TestFace( test13[_config][2] ) ) _subconfig += 4 ;
//if( TestFace( test13[_config][3] ) ) _subconfig += 8 ;
//if( TestFace( test13[_config][4] ) ) _subconfig += 16 ;
//if( TestFace( test13[_config][5] ) ) _subconfig += 32 ;
if( TestFace( MCLookUpTable::Test13(_config, 0) ) ) _subconfig += 1 ;
if( TestFace( MCLookUpTable::Test13(_config, 1) ) ) _subconfig += 2 ;
if( TestFace( MCLookUpTable::Test13(_config, 2) ) ) _subconfig += 4 ;
if( TestFace( MCLookUpTable::Test13(_config, 3) ) ) _subconfig += 8 ;
if( TestFace( MCLookUpTable::Test13(_config, 4) ) ) _subconfig += 16 ;
if( TestFace( MCLookUpTable::Test13(_config, 5) ) ) _subconfig += 32 ;
switch( MCLookUpTable::Subconfig13(_subconfig) ) //switch( subconfig13[_subconfig] )
{
case 0 : { /* 13.1 */ AddTriangles( MCLookUpTable::Tiling13_1(_config) , 4 ) ; break ; }
case 1 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 0), 6 ) ; break ; }
case 2 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 1), 6 ) ; break ; }
case 3 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 2), 6 ) ; break ; }
case 4 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 3), 6 ) ; break ; }
case 5 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 4), 6 ) ; break ; }
case 6 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 5), 6 ) ; break ; }
case 7 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 0), 10, v12 ) ; break ; }
case 8 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 1), 10, v12 ) ; break ; }
case 9 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 2), 10, v12 ) ; break ; }
case 10 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 3), 10, v12 ) ; break ; }
case 11 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 4), 10, v12 ) ; break ; }
case 12 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 5), 10, v12 ) ; break ; }
case 13 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 6), 10, v12 ) ; break ; }
case 14 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 7), 10, v12 ) ; break ; }
case 15 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 8), 10, v12 ) ; break ; }
case 16 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 9), 10, v12 ) ; break ; }
case 17 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config,10), 10, v12 ) ; break ; }
case 18 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config,11), 10, v12 ) ; break ; }
case 19 : { /* 13.4 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_4(_config, 0), 12, v12 ) ; break ; }
case 20 : { /* 13.4 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_4(_config, 1), 12, v12 ) ; break ; }
case 21 : { /* 13.4 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_4(_config, 2), 12, v12 ) ; break ; }
case 22 : { /* 13.4 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_4(_config, 3), 12, v12 ) ; break ; }
case 23 :
{ /* 13.5 */
_subconfig = 0 ;
if( TestInterior( MCLookUpTable::Test13(_config, 6) ) )
AddTriangles( MCLookUpTable::Tiling13_5_1(_config, 0), 6 ) ;
else
AddTriangles( MCLookUpTable::Tiling13_5_2(_config, 0), 10 ) ;
break ;
}
case 24 :
{ /* 13.5 */
_subconfig = 1 ;
if( TestInterior( MCLookUpTable::Test13(_config, 6) ) )
AddTriangles( MCLookUpTable::Tiling13_5_1(_config, 1), 6 ) ;
else
AddTriangles( MCLookUpTable::Tiling13_5_2(_config, 1), 10 ) ;
break ;
}
case 25 :
{/* 13.5 */
_subconfig = 2 ;
if( TestInterior( MCLookUpTable::Test13(_config, 6) ) )
AddTriangles( MCLookUpTable::Tiling13_5_1(_config, 2), 6 ) ;
else
AddTriangles( MCLookUpTable::Tiling13_5_2(_config, 2), 10 ) ;
break ;
}
case 26 :
{/* 13.5 */
_subconfig = 3 ;
if( TestInterior( MCLookUpTable::Test13(_config, 6) ) )
AddTriangles( MCLookUpTable::Tiling13_5_1(_config, 3), 6 ) ;
else
AddTriangles( MCLookUpTable::Tiling13_5_2(_config, 3), 10 ) ;
break ;
}
case 27 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 0), 10, v12 ) ; break ; }
case 28 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 1), 10, v12 ) ; break ; }
case 29 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 2), 10, v12 ) ; break ; }
case 30 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 3), 10, v12 ) ; break ; }
case 31 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 4), 10, v12 ) ; break ; }
case 32 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 5), 10, v12 ) ; break ; }
case 33 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 6), 10, v12 ) ; break ; }
case 34 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 7), 10, v12 ) ; break ; }
case 35 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 8), 10, v12 ) ; break ; }
case 36 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 9), 10, v12 ) ; break ; }
case 37 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config,10), 10, v12 ) ; break ; }
case 38 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config,11), 10, v12 ) ; break ; }
case 39 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,0), 6 ) ; break ; }
case 40 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,1), 6 ) ; break ; }
case 41 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,2), 6 ) ; break ; }
case 42 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,3), 6 ) ; break ; }
case 43 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,4), 6 ) ; break ; }
case 44 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,5), 6 ) ; break ; }
case 45 : { /* 13.1 */ AddTriangles( MCLookUpTable::Tiling13_1_(_config) , 4 ) ; break ; }
default : { /*Impossible case 13*/ assert(false); }
}
break ;
} // end of case 13
case 14 : { AddTriangles( MCLookUpTable::Tiling14(_config), 4 ) ; }
break ;
} //end of switch (_case)
}; // end of ApplyMC
private:
/*!
*/
WALKER_TYPE *_walker;
/*!
*/
TRIMESH_TYPE *_mesh;
/*!
* The field value at the cell corners
*/
ScalarType _field[8];
/*!
* Array of the 8 corners of the volume cell being processed
*/
vcg::Point3i _corners[8];
/*!
* Case of the volume cell being processed
*/
unsigned char _case;
/*!
* Configuration of the volume cell being processed
*/
unsigned char _config;
/*!
* Subconfiguration of the volume cell being processed
*/
unsigned char _subconfig;
/*!
* Tests if the components of the tesselation of the cube should be connected
* by the interior of an ambiguous face
*/
inline bool TestFace(signed char face)
{
ScalarType A,B,C,D ;
switch( face )
{
case -1 : case 1 : A = _field[0] ; B = _field[4] ; C = _field[5] ; D = _field[1] ; break ;
case -2 : case 2 : A = _field[1] ; B = _field[5] ; C = _field[6] ; D = _field[2] ; break ;
case -3 : case 3 : A = _field[2] ; B = _field[6] ; C = _field[7] ; D = _field[3] ; break ;
case -4 : case 4 : A = _field[3] ; B = _field[7] ; C = _field[4] ; D = _field[0] ; break ;
case -5 : case 5 : A = _field[0] ; B = _field[3] ; C = _field[2] ; D = _field[1] ; break ;
case -6 : case 6 : A = _field[4] ; B = _field[7] ; C = _field[6] ; D = _field[5] ; break ;
default : assert(false); // Invalid face code
};
return face * A * ( A*C - B*D ) >= 0 ; // face and A invert signs
}; // end of TestFace
/*!
* Tests if the components of the tesselation of the cube should be connected
* through the interior of the cube
*/
inline bool TestInterior(signed char s)
{
ScalarType t, At=0, Bt=0, Ct=0, Dt=0, a, b ;
char test = 0 ;
char edge = -1 ; // reference edge of the triangulation
switch( _case )
{
case 4 :
case 10 :
{
a = (_field[4]-_field[0])*(_field[6]-_field[2]) - (_field[7]-_field[3])*(_field[5]-_field[1]);
b = _field[2]*(_field[4]-_field[0])+_field[0]*(_field[6]-_field[2])-_field[1]*(_field[7]-_field[3])-_field[3]*(_field[5]-_field[1]);
t = - b / (2*a) ;
if( t<0 || t>1 )
return s>0 ;
At = _field[0] + ( _field[4] - _field[0] ) * t ;
Bt = _field[3] + ( _field[7] - _field[3] ) * t ;
Ct = _field[2] + ( _field[6] - _field[2] ) * t ;
Dt = _field[1] + ( _field[5] - _field[1] ) * t ;
break ;
}
case 6 :
case 7 :
case 12 :
case 13 :
switch( _case )
{
case 6 : edge = MCLookUpTable::Test6 (_config, 2) ; break ;
case 7 : edge = MCLookUpTable::Test7 (_config, 4) ; break ;
case 12 : edge = MCLookUpTable::Test12(_config, 3) ; break ;
case 13 : edge = MCLookUpTable::Tiling13_5_1(_config, _subconfig)[0] ; break ;
}
switch( edge )
{
case 0 :
t = _field[0] / ( _field[0] - _field[1] ) ;
At = 0 ;
Bt = _field[3] + ( _field[2] - _field[3] ) * t ;
Ct = _field[7] + ( _field[6] - _field[7] ) * t ;
Dt = _field[4] + ( _field[5] - _field[4] ) * t ;
break ;
case 1 :
t = _field[1] / ( _field[1] - _field[2] ) ;
At = 0 ;
Bt = _field[0] + ( _field[3] - _field[0] ) * t ;
Ct = _field[4] + ( _field[7] - _field[4] ) * t ;
Dt = _field[5] + ( _field[6] - _field[5] ) * t ;
break ;
case 2 :
t = _field[2] / ( _field[2] - _field[3] ) ;
At = 0 ;
Bt = _field[1] + ( _field[0] - _field[1] ) * t ;
Ct = _field[5] + ( _field[4] - _field[5] ) * t ;
Dt = _field[6] + ( _field[7] - _field[6] ) * t ;
break ;
case 3 :
t = _field[3] / ( _field[3] - _field[0] ) ;
At = 0 ;
Bt = _field[2] + ( _field[1] - _field[2] ) * t ;
Ct = _field[6] + ( _field[5] - _field[6] ) * t ;
Dt = _field[7] + ( _field[4] - _field[7] ) * t ;
break ;
case 4 :
t = _field[4] / ( _field[4] - _field[5] ) ;
At = 0 ;
Bt = _field[7] + ( _field[6] - _field[7] ) * t ;
Ct = _field[3] + ( _field[2] - _field[3] ) * t ;
Dt = _field[0] + ( _field[1] - _field[0] ) * t ;
break ;
case 5 :
t = _field[5] / ( _field[5] - _field[6] ) ;
At = 0 ;
Bt = _field[4] + ( _field[7] - _field[4] ) * t ;
Ct = _field[0] + ( _field[3] - _field[0] ) * t ;
Dt = _field[1] + ( _field[2] - _field[1] ) * t ;
break ;
case 6 :
t = _field[6] / ( _field[6] - _field[7] ) ;
At = 0 ;
Bt = _field[5] + ( _field[4] - _field[5] ) * t ;
Ct = _field[1] + ( _field[0] - _field[1] ) * t ;
Dt = _field[2] + ( _field[3] - _field[2] ) * t ;
break ;
case 7 :
t = _field[7] / ( _field[7] - _field[4] ) ;
At = 0 ;
Bt = _field[6] + ( _field[5] - _field[6] ) * t ;
Ct = _field[2] + ( _field[1] - _field[2] ) * t ;
Dt = _field[3] + ( _field[0] - _field[3] ) * t ;
break ;
case 8 :
t = _field[0] / ( _field[0] - _field[4] ) ;
At = 0 ;
Bt = _field[3] + ( _field[7] - _field[3] ) * t ;
Ct = _field[2] + ( _field[6] - _field[2] ) * t ;
Dt = _field[1] + ( _field[5] - _field[1] ) * t ;
break ;
case 9 :
t = _field[1] / ( _field[1] - _field[5] ) ;
At = 0 ;
Bt = _field[0] + ( _field[4] - _field[0] ) * t ;
Ct = _field[3] + ( _field[7] - _field[3] ) * t ;
Dt = _field[2] + ( _field[6] - _field[2] ) * t ;
break ;
case 10 :
t = _field[2] / ( _field[2] - _field[6] ) ;
At = 0 ;
Bt = _field[1] + ( _field[5] - _field[1] ) * t ;
Ct = _field[0] + ( _field[4] - _field[0] ) * t ;
Dt = _field[3] + ( _field[7] - _field[3] ) * t ;
break ;
case 11 :
t = _field[3] / ( _field[3] - _field[7] ) ;
At = 0 ;
Bt = _field[2] + ( _field[6] - _field[2] ) * t ;
Ct = _field[1] + ( _field[5] - _field[1] ) * t ;
Dt = _field[0] + ( _field[4] - _field[0] ) * t ;
break ;
default: { assert(false); /* Invalid edge */ break ; }
}
break ;
default : assert(false); /* Invalid ambiguous case */ break;
}
if( At >= 0 ) test ++ ;
if( Bt >= 0 ) test += 2 ;
if( Ct >= 0 ) test += 4 ;
if( Dt >= 0 ) test += 8 ;
switch( test )
{
case 0 : return s>0 ;
case 1 : return s>0 ;
case 2 : return s>0 ;
case 3 : return s>0 ;
case 4 : return s>0 ;
case 5 : if( At * Ct < Bt * Dt ) return s>0 ; break ;
case 6 : return s>0 ;
case 7 : return s<0 ;
case 8 : return s>0 ;
case 9 : return s>0 ;
case 10 : if( At * Ct >= Bt * Dt ) return s>0 ; break ;
case 11 : return s<0 ;
case 12 : return s>0 ;
case 13 : return s<0 ;
case 14 : return s<0 ;
case 15 : return s<0 ;
}
return s<0 ;
}; //end of TestInterior
/*!
* Adds a vertex inside the current cube
* \param v The pointer to the new vertex along the edge
*/
inline void ComputeCVertex(VertexPointer &v12)
{
v12 = &*AllocatorType::AddVertices(*_mesh, 1);
v12->P() = CoordType(0.0, 0.0, 0.0);
unsigned int count = 0;
VertexPointer v = NULL;
if (_walker->Exist(_corners[0], _corners[1], v) )
{
count++;
v12->P() += v->P();
}
if (_walker->Exist(_corners[1], _corners[2], v) )
{
count++;
v12->P() += v->P();
}
if (_walker->Exist(_corners[3], _corners[2], v) )
{
count++;
v12->P() += v->P();
}
if (_walker->Exist(_corners[0], _corners[3], v) )
{
count++;
v12->P() += v->P();
}
if (_walker->Exist(_corners[4], _corners[5], v) )
{
count++;
v12->P() += v->P();
}
if (_walker->Exist(_corners[5], _corners[6], v) )
{
count++;
v12->P() += v->P();
}
if (_walker->Exist(_corners[7], _corners[6], v) )
{
count++;
v12->P() += v->P();
}
if (_walker->Exist(_corners[4], _corners[7], v) )
{
count++;
v12->P() += v->P();
}
if (_walker->Exist(_corners[0], _corners[4], v) )
{
count++;
v12->P() += v->P();
}
if (_walker->Exist(_corners[1], _corners[5], v) )
{
count++;
v12->P() += v->P();
}
if (_walker->Exist(_corners[2], _corners[6], v) )
{
count++;
v12->P() += v->P();
}
if (_walker->Exist(_corners[3], _corners[7], v) )
{
count++;
v12->P() += v->P();
}
v12->P() /= (float) count;
} // end of AddCVertex
/*!
* Adds new triangles to the mesh
* \param vertices_list The list of vertex indices
* \param n The number of triangles that will be added to the mesh
* \param v12 The pointer to the vertex inside the current cell
*/
inline void AddTriangles(const char *vertices_list, char n, VertexPointer v12=NULL)
{
VertexPointer vp = NULL;
size_t face_idx = _mesh->face.size();
size_t v12_idx = -1;
size_t vertices_idx[3];
if (v12 != NULL) v12_idx = v12 - &_mesh->vert[0];
AllocatorType::AddFaces(*_mesh, (int) n);
for (int trig=0; trig<3*n; face_idx++ )
{
vp = NULL;
memset(vertices_idx, -1, 3*sizeof(size_t));
for (int vert=0; vert<3; vert++, trig++) //ok
{
switch ( vertices_list[trig] )
{
case 0: { _walker->GetXIntercept(_corners[0], _corners[1], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 1: { _walker->GetYIntercept(_corners[1], _corners[2], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 2: { _walker->GetXIntercept(_corners[3], _corners[2], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 3: { _walker->GetYIntercept(_corners[0], _corners[3], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 4: { _walker->GetXIntercept(_corners[4], _corners[5], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 5: { _walker->GetYIntercept(_corners[5], _corners[6], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 6: { _walker->GetXIntercept(_corners[7], _corners[6], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 7: { _walker->GetYIntercept(_corners[4], _corners[7], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 8: { _walker->GetZIntercept(_corners[0], _corners[4], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 9: { _walker->GetZIntercept(_corners[1], _corners[5], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 10: { _walker->GetZIntercept(_corners[2], _corners[6], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 11: { _walker->GetZIntercept(_corners[3], _corners[7], vp); vertices_idx[vert] = vp - &_mesh->vert[0]; break; }
case 12: { assert(v12 != NULL); vertices_idx[vert] = v12_idx; break; }
default: { assert(false); /* Invalid edge identifier */ }
} // end of switch
// Note that vp can be zero if we are in case 12 and that vertices_idx is surely >0 so the following assert has to be corrected as below.
// assert((vp - &_mesh->vert[0])>=0 && vertices_idx[vert]<_mesh->vert.size());
assert(vertices_idx[vert]<_mesh->vert.size());
} // end for (int vert=0 ...)
_mesh->face[face_idx].V(0) = &_mesh->vert[vertices_idx[0]];
_mesh->face[face_idx].V(1) = &_mesh->vert[vertices_idx[1]];
_mesh->face[face_idx].V(2) = &_mesh->vert[vertices_idx[2]];
} // end for (int trig=0...)
}; // end of AddTriangles
}; // end of class MarchingCubes
/*! @} */
//end of Doxygen documentation
}; // end of namespace tri
}; // end of namespace vcg
#endif //__VCG_MARCHING_CUBES

2560
vcg/complex/trimesh/create/mc_lookup_table.h
View File

File diff suppressed because it is too large Load Diff

346
vcg/complex/trimesh/create/mc_trivial_walker.h
View File

@ -1,346 +0,0 @@
/****************************************************************************
* VCGLib o o *
* Visual and Computer Graphics Library o o *
* _ O _ *
* Copyright(C) 2004-2009 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
#ifndef __VCG_TRIVIAL_WALKER
#define __VCG_TRIVIAL_WALKER
#include <wrap/callback.h>
namespace vcg {
// Very simple volume class.
// just an example of the interface that the trivial walker expects
template <class VOX_TYPE>
class SimpleVolume
{
public:
typedef VOX_TYPE VoxelType;
std::vector<VoxelType> Vol;
Point3i sz; /// Dimensioni griglia come numero di celle per lato
const Point3i &ISize() {return sz;}; /// Dimensioni griglia come numero di celle per lato
void Init(Point3i _sz)
{
sz=_sz;
Vol.resize(sz[0]*sz[1]*sz[2]);
}
float Val(const int &x,const int &y,const int &z) const {
return cV(x,y,z).V();
//else return numeric_limits<float>::quiet_NaN( );
}
float &Val(const int &x,const int &y,const int &z) {
return V(x,y,z).V();
//else return numeric_limits<float>::quiet_NaN( );
}
VOX_TYPE &V(const int &x,const int &y,const int &z) {
return Vol[x+y*sz[0]+z*sz[0]*sz[1]];
}
const VOX_TYPE &cV(const int &x,const int &y,const int &z) const {
return Vol[x+y*sz[0]+z*sz[0]*sz[1]];
}
typedef enum { XAxis=0,YAxis=1,ZAxis=2} VolumeAxis;
template < class VertexPointerType, VolumeAxis AxisVal >
void GetIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr)
{
float f1 = Val(p1.X(), p1.Y(), p1.Z())-thr;
float f2 = Val(p2.X(), p2.Y(), p2.Z())-thr;
float u = (float) f1/(f1-f2);
if(AxisVal==XAxis) v->P().X() = (float) p1.X()*(1-u) + u*p2.X();
else v->P().X() = (float) p1.X();
if(AxisVal==YAxis) v->P().Y() = (float) p1.Y()*(1-u) + u*p2.Y();
else v->P().Y() = (float) p1.Y();
if(AxisVal==ZAxis) v->P().Z() = (float) p1.Z()*(1-u) + u*p2.Z();
else v->P().Z() = (float) p1.Z();
}
template < class VertexPointerType >
void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr)
{ GetIntercept<VertexPointerType,XAxis>(p1,p2,v,thr); }
template < class VertexPointerType >
void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr)
{ GetIntercept<VertexPointerType,YAxis>(p1,p2,v,thr); }
template < class VertexPointerType >
void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr)
{ GetIntercept<VertexPointerType,ZAxis>(p1,p2,v,thr); }
};
template <class VolumeType>
class RawVolumeImporter
{
public:
enum DataType
{
// Funzioni superiori
UNDEF=0,
BYTE=1,
SHORT=2,
FLOAT=3
};
static bool Open(const char *filename, VolumeType &V, Point3i sz, DataType d)
{
return true;
}
};
class SimpleVoxel
{
private:
float _v;
public:
float &V() {return _v;};
float V() const {return _v;};
};
namespace tri {
// La classe Walker implementa la politica di visita del volume; conoscendo l'ordine di visita del volume
// Ë conveniente che il Walker stesso si faccia carico del caching dei dati utilizzati durante l'esecuzione
// degli algoritmi MarchingCubes ed ExtendedMarchingCubes, in particolare il calcolo del volume ai vertici
// delle celle e delle intersezioni della superficie con le celle. In questo esempio il volume da processare
// viene suddiviso in fette; in questo modo se il volume ha dimensione h*l*w (rispettivamente altezza,
// larghezza e profondit‡), lo spazio richiesto per il caching dei vertici gi‡ allocati passa da O(h*l*w)
// a O(h*l).
template <class MeshType, class VolumeType>
class TrivialWalker
{
private:
typedef int VertexIndex;
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexPointer VertexPointer;
public:
// bbox is the portion of the volume to be computed
// resolution determine the sampling step:
// should be a divisor of bbox size (e.g. if bbox size is 256^3 resolution could be 128,64, etc)
void Init(VolumeType &volume)
{
_bbox = Box3i(Point3i(0,0,0),volume.ISize());
_slice_dimension = _bbox.DimX()*_bbox.DimZ();
_x_cs = new VertexIndex[ _slice_dimension ];
_y_cs = new VertexIndex[ _slice_dimension ];
_z_cs = new VertexIndex[ _slice_dimension ];
_x_ns = new VertexIndex[ _slice_dimension ];
_z_ns = new VertexIndex[ _slice_dimension ];
};
~TrivialWalker()
{_thr=0;}
template<class EXTRACTOR_TYPE>
void BuildMesh(MeshType &mesh, VolumeType &volume, EXTRACTOR_TYPE &extractor, const float threshold, vcg::CallBackPos * cb=0)
{
Init(volume);
_volume = &volume;
_mesh = &mesh;
_mesh->Clear();
_thr=threshold;
vcg::Point3i p1, p2;
Begin();
extractor.Initialize();
for (int j=_bbox.min.Y(); j<(_bbox.max.Y()-1)-1; j+=1)
{
if(cb && ((j%10)==0) ) cb(j*_bbox.DimY()/100.0,"Marching volume");
for (int i=_bbox.min.X(); i<(_bbox.max.X()-1)-1; i+=1)
{
for (int k=_bbox.min.Z(); k<(_bbox.max.Z()-1)-1; k+=1)
{
p1.X()=i; p1.Y()=j; p1.Z()=k;
p2.X()=i+1; p2.Y()=j+1; p2.Z()=k+1;
extractor.ProcessCell(p1, p2);
}
}
NextSlice();
}
extractor.Finalize();
_volume = NULL;
_mesh = NULL;
};
float V(int pi, int pj, int pk)
{
return _volume->Val(pi, pj, pk)-_thr;
}
bool Exist(const vcg::Point3i &p0, const vcg::Point3i &p1, VertexPointer &v)
{
int pos = p0.X()+p0.Z()*_bbox.max.X();
int vidx;
if (p0.X()!=p1.X()) // punti allineati lungo l'asse X
vidx = (p0.Y()==_current_slice) ? _x_cs[pos] : _x_ns[pos];
else if (p0.Y()!=p1.Y()) // punti allineati lungo l'asse Y
vidx = _y_cs[pos];
else if (p0.Z()!=p1.Z()) // punti allineati lungo l'asse Z
vidx = (p0.Y()==_current_slice)? _z_cs[pos] : _z_ns[pos];
else
assert(false);
v = (vidx!=-1)? &_mesh->vert[vidx] : NULL;
return v!=NULL;
}
void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
{
int i = p1.X() - _bbox.min.X();
int z = p1.Z() - _bbox.min.Z();
VertexIndex index = i+z*_bbox.max.X();
VertexIndex pos;
if (p1.Y()==_current_slice)
{
if ((pos=_x_cs[index])==-1)
{
_x_cs[index] = (VertexIndex) _mesh->vert.size();
pos = _x_cs[index];
Allocator<MeshType>::AddVertices( *_mesh, 1 );
v = &_mesh->vert[pos];
_volume->GetXIntercept(p1, p2, v, _thr);
return;
}
}
if (p1.Y()==_current_slice+1)
{
if ((pos=_x_ns[index])==-1)
{
_x_ns[index] = (VertexIndex) _mesh->vert.size();
pos = _x_ns[index];
Allocator<MeshType>::AddVertices( *_mesh, 1 );
v = &_mesh->vert[pos];
_volume->GetXIntercept(p1, p2, v,_thr);
return;
}
}
assert(pos >=0 && size_t(pos)< _mesh->vert.size());
v = &_mesh->vert[pos];
}
void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
{
int i = p1.X() - _bbox.min.X();
int z = p1.Z() - _bbox.min.Z();
VertexIndex index = i+z*_bbox.max.X();
VertexIndex pos;
if ((pos=_y_cs[index])==-1)
{
_y_cs[index] = (VertexIndex) _mesh->vert.size();
pos = _y_cs[index];
Allocator<MeshType>::AddVertices( *_mesh, 1);
v = &_mesh->vert[ pos ];
_volume->GetYIntercept(p1, p2, v,_thr);
}
v = &_mesh->vert[pos];
}
void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
{
int i = p1.X() - _bbox.min.X();
int z = p1.Z() - _bbox.min.Z();
VertexIndex index = i+z*_bbox.max.X();
VertexIndex pos;
if (p1.Y()==_current_slice)
{
if ((pos=_z_cs[index])==-1)
{
_z_cs[index] = (VertexIndex) _mesh->vert.size();
pos = _z_cs[index];
Allocator<MeshType>::AddVertices( *_mesh, 1 );
v = &_mesh->vert[pos];
_volume->GetZIntercept(p1, p2, v,_thr);
return;
}
}
if (p1.Y()==_current_slice+1)
{
if ((pos=_z_ns[index])==-1)
{
_z_ns[index] = (VertexIndex) _mesh->vert.size();
pos = _z_ns[index];
Allocator<MeshType>::AddVertices( *_mesh, 1 );
v = &_mesh->vert[pos];
_volume->GetZIntercept(p1, p2, v,_thr);
return;
}
}
v = &_mesh->vert[pos];
}
protected:
Box3i _bbox;
int _slice_dimension;
int _current_slice;
VertexIndex *_x_cs; // indici dell'intersezioni della superficie lungo gli Xedge della fetta corrente
VertexIndex *_y_cs; // indici dell'intersezioni della superficie lungo gli Yedge della fetta corrente
VertexIndex *_z_cs; // indici dell'intersezioni della superficie lungo gli Zedge della fetta corrente
VertexIndex *_x_ns; // indici dell'intersezioni della superficie lungo gli Xedge della prossima fetta
VertexIndex *_z_ns; // indici dell'intersezioni della superficie lungo gli Zedge della prossima fetta
MeshType *_mesh;
VolumeType *_volume;
float _thr;
void NextSlice()
{
memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex));
memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex));
memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex));
std::swap(_x_cs, _x_ns);
std::swap(_z_cs, _z_ns);
_current_slice += 1;
}
void Begin()
{
_current_slice = _bbox.min.Y();
memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex));
memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex));
memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex));
memset(_x_ns, -1, _slice_dimension*sizeof(VertexIndex));
memset(_z_ns, -1, _slice_dimension*sizeof(VertexIndex));
}
};
} // end namespace
} // end namespace
#endif // __VCGTEST_WALKER

829
vcg/complex/trimesh/create/platonic.h
View File

@ -1,829 +0,0 @@
/****************************************************************************
* VCGLib o o *
* Visual and Computer Graphics Library o o *
* _ O _ *
* Copyright(C) 2004 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
#ifndef __VCGLIB_PLATONIC
#define __VCGLIB_PLATONIC
#include<vcg/math/base.h>
#include<vcg/complex/trimesh/allocate.h>
#include<vcg/complex/trimesh/refine.h>
#include<vcg/complex/trimesh/update/flag.h>
namespace vcg {
namespace tri {
/** \addtogroup trimesh */
//@{
/**
A set of functions that builds meshes
that represent surfaces of platonic solids,
and other simple shapes.
The 1st parameter is the mesh that will
be filled with the solid.
*/
template <class TetraMeshType>
void Tetrahedron(TetraMeshType &in)
{
typedef TetraMeshType MeshType;
typedef typename TetraMeshType::CoordType CoordType;
typedef typename TetraMeshType::VertexPointer VertexPointer;
typedef typename TetraMeshType::VertexIterator VertexIterator;
typedef typename TetraMeshType::FaceIterator FaceIterator;
in.Clear();
Allocator<TetraMeshType>::AddVertices(in,4);
Allocator<TetraMeshType>::AddFaces(in,4);
VertexPointer ivp[4];
VertexIterator vi=in.vert.begin();
ivp[0]=&*vi;(*vi).P()=CoordType ( 1.0, 1.0, 1.0); ++vi;
ivp[1]=&*vi;(*vi).P()=CoordType (-1.0, 1.0,-1.0); ++vi;
ivp[2]=&*vi;(*vi).P()=CoordType (-1.0,-1.0, 1.0); ++vi;
ivp[3]=&*vi;(*vi).P()=CoordType ( 1.0,-1.0,-1.0);
FaceIterator fi=in.face.begin();
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[1]; (*fi).V(2)=ivp[2]; ++fi;
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[2]; (*fi).V(2)=ivp[3]; ++fi;
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[3]; (*fi).V(2)=ivp[1]; ++fi;
(*fi).V(0)=ivp[3]; (*fi).V(1)=ivp[2]; (*fi).V(2)=ivp[1];
}
/// builds a Dodecahedron,
/// (each pentagon is composed of 5 triangles)
template <class DodMeshType>
void Dodecahedron(DodMeshType & in)
{
typedef DodMeshType MeshType;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FaceIterator FaceIterator;
typedef typename MeshType::ScalarType ScalarType;
const int N_penta=12;
const int N_points=62;
int penta[N_penta*3*3]=
{20,11, 18, 18, 11, 8, 8, 11, 4,
13,23, 4, 4, 23, 8, 8, 23, 16,
13, 4, 30, 30, 4, 28, 28, 4, 11,
16,34, 8, 8, 34, 18, 18, 34, 36,
11,20, 28, 28, 20, 45, 45, 20, 38,
13,30, 23, 23, 30, 41, 41, 30, 47,
16,23, 34, 34, 23, 50, 50, 23, 41,
20,18, 38, 38, 18, 52, 52, 18, 36,
30,28, 47, 47, 28, 56, 56, 28, 45,
50,60, 34, 34, 60, 36, 36, 60, 52,
45,38, 56, 56, 38, 60, 60, 38, 52,
50,41, 60, 60, 41, 56, 56, 41, 47 };
//A B E D C
const ScalarType p=(1.0 + math::Sqrt(5.0)) / 2.0;
const ScalarType p2=p*p;
const ScalarType p3=p*p*p;
ScalarType vv[N_points*3]=
{
0, 0, 2*p2, p2, 0, p3, p, p2, p3,
0, p, p3, -p, p2, p3, -p2, 0, p3,
-p, -p2, p3, 0, -p, p3, p, -p2, p3,
p3, p, p2, p2, p2, p2, 0, p3, p2,
-p2, p2, p2, -p3, p, p2, -p3, -p, p2,
-p2, -p2, p2, 0, -p3, p2, p2, -p2, p2,
p3, -p, p2, p3, 0, p, p2, p3, p,
-p2, p3, p, -p3, 0, p, -p2, -p3, p,
p2, -p3, p, 2*p2, 0, 0, p3, p2, 0,
p, p3, 0, 0, 2*p2, 0, -p, p3, 0,
-p3, p2, 0, -2*p2, 0, 0, -p3, -p2, 0,
-p, -p3, 0, 0, -2*p2, 0, p, -p3, 0,
p3, -p2, 0, p3, 0, -p, p2, p3, -p,
-p2, p3, -p, -p3, 0, -p, -p2, -p3, -p,
p2, -p3, -p, p3, p, -p2, p2, p2, -p2,
0, p3, -p2, -p2, p2, -p2, -p3, p, -p2,
-p3, -p, -p2, -p2, -p2, -p2, 0, -p3, -p2,
p2, -p2, -p2, p3, -p, -p2, p2, 0, -p3,
p, p2, -p3, 0, p, -p3, -p, p2, -p3,
-p2, 0, -p3, -p, -p2, -p3, 0, -p, -p3,
p, -p2, -p3, 0, 0, -2*p2
};
in.Clear();
//in.face.clear();
Allocator<DodMeshType>::AddVertices(in,20+12);
Allocator<DodMeshType>::AddFaces(in, 5*12); // five pentagons, each made by 5 tri
int h,i,j,m=0;
bool used[N_points];
for (i=0; i<N_points; i++) used[i]=false;
int reindex[20+12 *10];
ScalarType xx,yy,zz, sx,sy,sz;
int order[5]={0,1,8,6,2};
int added[12];
VertexIterator vi=in.vert.begin();
for (i=0; i<12; i++) {
sx=sy=sz=0;
for (int j=0; j<5; j++) {
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;
if (!used[h]) {
(*vi).P()=CoordType( xx, yy, zz ); vi++;
used[h]=true;
reindex[ h ] = m++;
}
}
(*vi).P()=CoordType( sx/5.0, sy/5.0, sz/5.0 ); vi++;
added[ i ] = m++;
}
std::vector<VertexPointer> index(in.vn);
for(j=0,vi=in.vert.begin();j<in.vn;++j,++vi) index[j] = &(*vi);
FaceIterator fi=in.face.begin();
for (i=0; i<12; i++) {
for (j=0; j<5; j++){
(*fi).V(0)=index[added[i] ];
(*fi).V(1)=index[reindex[penta[i*9 + order[j ] ] -1 ] ];
(*fi).V(2)=index[reindex[penta[i*9 + order[(j+1)%5] ] -1 ] ];
if (in.HasPerFaceFlags()) {
// tag faux edges
(*fi).SetF(0);
(*fi).SetF(2);
}
fi++;
}
}
}
template <class OctMeshType>
void Octahedron(OctMeshType &in)
{
typedef OctMeshType MeshType;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FaceIterator FaceIterator;
in.Clear();
Allocator<OctMeshType>::AddVertices(in,6);
Allocator<OctMeshType>::AddFaces(in,8);
VertexPointer ivp[6];
VertexIterator vi=in.vert.begin();
ivp[0]=&*vi;(*vi).P()=CoordType ( 1, 0, 0); ++vi;
ivp[1]=&*vi;(*vi).P()=CoordType ( 0, 1, 0); ++vi;
ivp[2]=&*vi;(*vi).P()=CoordType ( 0, 0, 1); ++vi;
ivp[3]=&*vi;(*vi).P()=CoordType (-1, 0, 0); ++vi;
ivp[4]=&*vi;(*vi).P()=CoordType ( 0,-1, 0); ++vi;
ivp[5]=&*vi;(*vi).P()=CoordType ( 0, 0,-1);
FaceIterator fi=in.face.begin();
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[1]; (*fi).V(2)=ivp[2]; ++fi;
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[2]; (*fi).V(2)=ivp[4]; ++fi;
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[4]; (*fi).V(2)=ivp[5]; ++fi;
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[5]; (*fi).V(2)=ivp[1]; ++fi;
(*fi).V(0)=ivp[3]; (*fi).V(1)=ivp[1]; (*fi).V(2)=ivp[5]; ++fi;
(*fi).V(0)=ivp[3]; (*fi).V(1)=ivp[5]; (*fi).V(2)=ivp[4]; ++fi;
(*fi).V(0)=ivp[3]; (*fi).V(1)=ivp[4]; (*fi).V(2)=ivp[2]; ++fi;
(*fi).V(0)=ivp[3]; (*fi).V(1)=ivp[2]; (*fi).V(2)=ivp[1];
}
template <class IcoMeshType>
void Icosahedron(IcoMeshType &in)
{
typedef IcoMeshType MeshType;
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FaceIterator FaceIterator;
ScalarType L=ScalarType((math::Sqrt(5.0)+1.0)/2.0);
CoordType vv[12]={
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 ( 1, 0, L),
CoordType (-1, 0, L),
CoordType ( 1, 0,-L),
CoordType (-1, 0,-L)
};
int ff[20][3]={
{1,0,4},{0,1,6},{2,3,5},{3,2,7},
{4,5,10},{5,4,8},{6,7,9},{7,6,11},
{8,9,2},{9,8,0},{10,11,1},{11,10,3},
{0,8,4},{0,6,9},{1,4,10},{1,11,6},
{2,5,8},{2,9,7},{3,10,5},{3,7,11}
};
in.Clear();
Allocator<IcoMeshType>::AddVertices(in,12);
Allocator<IcoMeshType>::AddFaces(in,20);
VertexPointer ivp[12];
VertexIterator vi;
int i;
for(i=0,vi=in.vert.begin();vi!=in.vert.end();++i,++vi){
(*vi).P()=vv[i];
ivp[i]=&*vi;
}
FaceIterator fi;
for(i=0,fi=in.face.begin();fi!=in.face.end();++i,++fi){
(*fi).V(0)=ivp[ff[i][0]];
(*fi).V(1)=ivp[ff[i][1]];
(*fi).V(2)=ivp[ff[i][2]];
}
}
template <class MeshType>
void Hexahedron(MeshType &in)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FaceIterator FaceIterator;
in.Clear();
Allocator<MeshType>::AddVertices(in,8);
Allocator<MeshType>::AddFaces(in,12);
VertexPointer ivp[8];
VertexIterator vi=in.vert.begin();
ivp[7]=&*vi;(*vi).P()=CoordType (-1,-1,-1); ++vi;
ivp[6]=&*vi;(*vi).P()=CoordType ( 1,-1,-1); ++vi;
ivp[5]=&*vi;(*vi).P()=CoordType (-1, 1,-1); ++vi;
ivp[4]=&*vi;(*vi).P()=CoordType ( 1, 1,-1); ++vi;
ivp[3]=&*vi;(*vi).P()=CoordType (-1,-1, 1); ++vi;
ivp[2]=&*vi;(*vi).P()=CoordType ( 1,-1, 1); ++vi;
ivp[1]=&*vi;(*vi).P()=CoordType (-1, 1, 1); ++vi;
ivp[0]=&*vi;(*vi).P()=CoordType ( 1, 1, 1);
FaceIterator fi=in.face.begin();
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[1]; (*fi).V(2)=ivp[2]; ++fi;
(*fi).V(0)=ivp[3]; (*fi).V(1)=ivp[2]; (*fi).V(2)=ivp[1]; ++fi;
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[2]; (*fi).V(2)=ivp[4]; ++fi;
(*fi).V(0)=ivp[6]; (*fi).V(1)=ivp[4]; (*fi).V(2)=ivp[2]; ++fi;
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[4]; (*fi).V(2)=ivp[1]; ++fi;
(*fi).V(0)=ivp[5]; (*fi).V(1)=ivp[1]; (*fi).V(2)=ivp[4]; ++fi;
(*fi).V(0)=ivp[7]; (*fi).V(1)=ivp[5]; (*fi).V(2)=ivp[6]; ++fi;
(*fi).V(0)=ivp[4]; (*fi).V(1)=ivp[6]; (*fi).V(2)=ivp[5]; ++fi;
(*fi).V(0)=ivp[7]; (*fi).V(1)=ivp[6]; (*fi).V(2)=ivp[3]; ++fi;
(*fi).V(0)=ivp[2]; (*fi).V(1)=ivp[3]; (*fi).V(2)=ivp[6]; ++fi;
(*fi).V(0)=ivp[7]; (*fi).V(1)=ivp[3]; (*fi).V(2)=ivp[5]; ++fi;
(*fi).V(0)=ivp[1]; (*fi).V(1)=ivp[5]; (*fi).V(2)=ivp[3];
if (in.HasPerFaceFlags()) {
FaceIterator fi=in.face.begin();
for (int k=0; k<12; k++) {
(*fi).SetF(1); fi++;
}
}
}
template <class MeshType>
void Square(MeshType &in)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FaceIterator FaceIterator;
in.Clear();
Allocator<MeshType>::AddVertices(in,4);
Allocator<MeshType>::AddFaces(in,2);
VertexPointer ivp[4];
VertexIterator vi=in.vert.begin();
ivp[0]=&*vi;(*vi).P()=CoordType ( 1, 0, 0); ++vi;
ivp[1]=&*vi;(*vi).P()=CoordType ( 0, 1, 0); ++vi;
ivp[2]=&*vi;(*vi).P()=CoordType (-1, 0, 0); ++vi;
ivp[3]=&*vi;(*vi).P()=CoordType ( 0,-1, 0);
FaceIterator fi=in.face.begin();
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[1]; (*fi).V(2)=ivp[2]; ++fi;
(*fi).V(0)=ivp[2]; (*fi).V(1)=ivp[3]; (*fi).V(2)=ivp[0];
if (in.HasPerFaceFlags()) {
FaceIterator fi=in.face.begin();
for (int k=0; k<2; k++) {
(*fi).SetF(2); fi++;
}
}
}
// this function build a sphere starting from a eventually not empty mesh.
// If the mesh is not empty it is 'spherified' and used as base for the subdivision process.
// otherwise an icosahedron is used.
template <class MeshType>
void Sphere(MeshType &in, const int subdiv = 3 )
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FaceIterator FaceIterator;
if(in.vn==0 && in.fn==0) Icosahedron(in);
VertexIterator vi;
for(vi = in.vert.begin(); vi!=in.vert.end();++vi)
vi->P().Normalize();
tri::UpdateFlags<MeshType>::FaceBorderFromNone(in);
tri::UpdateTopology<MeshType>::FaceFace(in);
size_t lastsize = 0;
for(int i = 0 ; i < subdiv; ++i)
{
Refine< MeshType, MidPoint<MeshType> >(in, MidPoint<MeshType>(&in), 0);
for(vi = in.vert.begin() + lastsize; vi != in.vert.end(); ++vi)
vi->P().Normalize();
lastsize = in.vert.size();
}
}
/// r1 = raggio 1, r2 = raggio2, h = altezza (asse y)
template <class MeshType>
void Cone( MeshType& in,
const typename MeshType::ScalarType r1,
const typename MeshType::ScalarType r2,
const typename MeshType::ScalarType h,
const int SubDiv = 36 )
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FaceIterator FaceIterator;
int i,b1,b2;
in.Clear();
int VN,FN;
if(r1==0 || r2==0) {
VN=SubDiv+2;
FN=SubDiv*2;
} else {
VN=SubDiv*2+2;
FN=SubDiv*4;
}
Allocator<MeshType>::AddVertices(in,VN);
Allocator<MeshType>::AddFaces(in,FN);
VertexPointer *ivp = new VertexPointer[VN];
VertexIterator vi=in.vert.begin();
ivp[0]=&*vi;(*vi).P()=CoordType ( 0,-h/2.0,0 ); ++vi;
ivp[1]=&*vi;(*vi).P()=CoordType ( 0, h/2.0,0 ); ++vi;
b1 = b2 = 2;
int cnt=2;
if(r1!=0)
{
for(i=0;i<SubDiv;++i)
{
double a = math::ToRad(i*360.0/SubDiv);
ivp[cnt]=&*vi; (*vi).P()= CoordType(r1*cos(a), -h/2.0, r1*sin(a)); ++vi;++cnt;
}
b2 += SubDiv;
}
if(r2!=0)
{
for(i=0;i<SubDiv;++i)
{
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;
}
}
FaceIterator fi=in.face.begin();
if(r1!=0) for(i=0;i<SubDiv;++i,++fi) {
(*fi).V(0)=ivp[0];
(*fi).V(1)=ivp[b1+i];
(*fi).V(2)=ivp[b1+(i+1)%SubDiv];
}
if(r2!=0) for(i=0;i<SubDiv;++i,++fi) {
(*fi).V(0)=ivp[1];
(*fi).V(2)=ivp[b2+i];
(*fi).V(1)=ivp[b2+(i+1)%SubDiv];
}
if(r1==0) for(i=0;i<SubDiv;++i,++fi)
{
(*fi).V(0)=ivp[0];
(*fi).V(1)=ivp[b2+i];
(*fi).V(2)=ivp[b2+(i+1)%SubDiv];
}
if(r2==0) for(i=0;i<SubDiv;++i,++fi){
(*fi).V(0)=ivp[1];
(*fi).V(2)=ivp[b1+i];
(*fi).V(1)=ivp[b1+(i+1)%SubDiv];
}
if(r1!=0 && r2!=0)for(i=0;i<SubDiv;++i)
{
(*fi).V(0)=ivp[b1+i];
(*fi).V(1)=ivp[b2+i];
(*fi).V(2)=ivp[b2+(i+1)%SubDiv];
++fi;
(*fi).V(0)=ivp[b1+i];
(*fi).V(1)=ivp[b2+(i+1)%SubDiv];
(*fi).V(2)=ivp[b1+(i+1)%SubDiv];
++fi;
}
}
template <class MeshType >
void Box(MeshType &in, const typename MeshType::BoxType & bb )
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FaceIterator FaceIterator;
in.Clear();
Allocator<MeshType>::AddVertices(in,8);
Allocator<MeshType>::AddFaces(in,12);
VertexPointer ivp[8];
VertexIterator vi=in.vert.begin();
ivp[0]=&*vi;(*vi).P()=CoordType (bb.min[0],bb.min[1],bb.min[2]); ++vi;
ivp[1]=&*vi;(*vi).P()=CoordType (bb.max[0],bb.min[1],bb.min[2]); ++vi;
ivp[2]=&*vi;(*vi).P()=CoordType (bb.min[0],bb.max[1],bb.min[2]); ++vi;
ivp[3]=&*vi;(*vi).P()=CoordType (bb.max[0],bb.max[1],bb.min[2]); ++vi;
ivp[4]=&*vi;(*vi).P()=CoordType (bb.min[0],bb.min[1],bb.max[2]); ++vi;
ivp[5]=&*vi;(*vi).P()=CoordType (bb.max[0],bb.min[1],bb.max[2]); ++vi;
ivp[6]=&*vi;(*vi).P()=CoordType (bb.min[0],bb.max[1],bb.max[2]); ++vi;
ivp[7]=&*vi;(*vi).P()=CoordType (bb.max[0],bb.max[1],bb.max[2]);
FaceIterator fi=in.face.begin();
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[1]; (*fi).V(2)=ivp[2]; ++fi;
(*fi).V(0)=ivp[3]; (*fi).V(1)=ivp[2]; (*fi).V(2)=ivp[1]; ++fi;
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[2]; (*fi).V(2)=ivp[4]; ++fi;
(*fi).V(0)=ivp[6]; (*fi).V(1)=ivp[4]; (*fi).V(2)=ivp[2]; ++fi;
(*fi).V(0)=ivp[0]; (*fi).V(1)=ivp[4]; (*fi).V(2)=ivp[1]; ++fi;
(*fi).V(0)=ivp[5]; (*fi).V(1)=ivp[1]; (*fi).V(2)=ivp[4]; ++fi;
(*fi).V(0)=ivp[7]; (*fi).V(1)=ivp[5]; (*fi).V(2)=ivp[6]; ++fi;
(*fi).V(0)=ivp[4]; (*fi).V(1)=ivp[6]; (*fi).V(2)=ivp[5]; ++fi;
(*fi).V(0)=ivp[7]; (*fi).V(1)=ivp[6]; (*fi).V(2)=ivp[3]; ++fi;
(*fi).V(0)=ivp[2]; (*fi).V(1)=ivp[3]; (*fi).V(2)=ivp[6]; ++fi;
(*fi).V(0)=ivp[7]; (*fi).V(1)=ivp[3]; (*fi).V(2)=ivp[5]; ++fi;
(*fi).V(0)=ivp[1]; (*fi).V(1)=ivp[5]; (*fi).V(2)=ivp[3];
if (in.HasPerFaceFlags()) {
FaceIterator fi=in.face.begin();
for (int k=0; k<12; k++) {
(*fi).SetF(1); fi++;
}
}
}
// this function build a mesh starting from a vector of generic coords (objects having a triple of float at their beginning)
// and a vector of faces (objects having a triple of ints at theri beginning).
template <class MeshType,class V, class F >
void Build( MeshType & in, const V & v, const F & f)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FaceIterator FaceIterator;
Allocator<MeshType>::AddVertices(in,v.size());
Allocator<MeshType>::AddFaces(in,f.size());
typename V::const_iterator vi;
typename MeshType::VertexType tv;
for(int i=0;i<v.size();++i)
{
float *vv=(float *)(&v[i]);
in.vert[i].P() = CoordType( vv[0],vv[1],vv[2]);
}
std::vector<VertexPointer> index(in.vn);
VertexIterator j;
int k;
for(k=0,j=in.vert.begin();j!=in.vert.end();++j,++k)
index[k] = &*j;
typename F::const_iterator fi;
typename MeshType::FaceType ft;
for(int i=0;i<f.size();++i)
{
int * ff=(int *)(&f[i]);
assert( ff[0]>=0 );
assert( ff[1]>=0 );
assert( ff[2]>=0 );
assert( ff[0]<in.vn );
assert( ff[1]<in.vn );
assert( ff[2]<in.vn );
in.face[i].V(0) = &in.vert[ ff[0] ];
in.face[i].V(1) = &in.vert[ ff[0] ];
in.face[i].V(2) = &in.vert[ ff[0] ];
}
}
// Build a regular grid mesh as a typical height field mesh
// x y are the position on the grid scaled by wl and hl (at the end x is in the range 0..wl and y is in 0..hl)
// z is taken from the <data> array
// Once generated the vertex positions it uses the FaceGrid function to generate the faces;
template <class MeshType>
void Grid(MeshType & in, int w, int h, float wl, float hl, float *data)
{
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FaceIterator FaceIterator;
in.Clear();
Allocator<MeshType>::AddVertices(in,w*h);
float wld=wl/float(w);
float hld=hl/float(h);
for(int i=0;i<h;++i)
for(int j=0;j<w;++j)
in.vert[i*w+j].P()=CoordType ( j*wld, i*hld, data[i*w+j]);
FaceGrid(in,w,h);
}
// Build a regular grid mesh of faces as a typical height field mesh
// Vertexes are assumed to be already be allocated.
template <class MeshType>
void FaceGrid(MeshType & in, int w, int h)
{
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
Allocator<MeshType>::AddFaces(in,(w-1)*(h-1)*2);
// i+0,j+0 -- i+0,j+1
// | \ |
// | \ |
// | \ |
// | \ |
// i+1,j+0 -- i+1,j+1
//
for(int i=0;i<h-1;++i)
for(int j=0;j<w-1;++j)
{
in.face[2*(i*(w-1)+j)+0].V(0) = &(in.vert[(i+1)*w+j+1]);
in.face[2*(i*(w-1)+j)+0].V(1) = &(in.vert[(i+0)*w+j+1]);
in.face[2*(i*(w-1)+j)+0].V(2) = &(in.vert[(i+0)*w+j+0]);
in.face[2*(i*(w-1)+j)+1].V(0) = &(in.vert[(i+0)*w+j+0]);
in.face[2*(i*(w-1)+j)+1].V(1) = &(in.vert[(i+1)*w+j+0]);
in.face[2*(i*(w-1)+j)+1].V(2) = &(in.vert[(i+1)*w+j+1]);
}
if (in.HasPerFaceFlags()) {
for (int k=0; k<(h-1)*(w-1)*2; k++) {
in.face[k].SetF(2);
}
}
}
// Build a regular grid mesh of faces as a typical height field mesh
// Vertexes are assumed to be already be allocated, but not oll the grid vertexes are present.
// For this purpos a grid of indexes is also passed. negative indexes means that there is no vertex.
template <class MeshType>
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 <= w*h); // the number of vertices should match the number of expected grid vertices
// V0 V1
// i+0,j+0 -- i+0,j+1
// | \ |
// | \ |
// | \ |
// | \ |
// i+1,j+0 -- i+1,j+1
// V2 V3
for(int i=0;i<h-1;++i)
for(int j=0;j<w-1;++j)
{
int V0i= grid[(i+0)*w+j+0];
int V1i= grid[(i+0)*w+j+1];
int V2i= grid[(i+1)*w+j+0];
int V3i= grid[(i+1)*w+j+1];
int ndone=0;
bool quad = (V0i>=0 && V1i>=0 && V2i>=0 && V3i>=0 ) && in.HasPerFaceFlags();
if(V0i>=0 && V2i>=0 && V3i>=0 )
{
typename MeshType::FaceIterator f= Allocator<MeshType>::AddFaces(in,1);
f->V(0)=&(in.vert[V3i]);
f->V(1)=&(in.vert[V2i]);
f->V(2)=&(in.vert[V0i]);
if (quad) f->SetF(2);
ndone++;
}
if(V0i>=0 && V1i>=0 && V3i>=0 )
{
typename MeshType::FaceIterator f= Allocator<MeshType>::AddFaces(in,1);
f->V(0)=&(in.vert[V0i]);
f->V(1)=&(in.vert[V1i]);
f->V(2)=&(in.vert[V3i]);
if (quad) f->SetF(2);
ndone++;
}
if (ndone==0) { // try diag the other way
if(V2i>=0 && V0i>=0 && V1i>=0 )
{
typename MeshType::FaceIterator f= Allocator<MeshType>::AddFaces(in,1);
f->V(0)=&(in.vert[V2i]);
f->V(1)=&(in.vert[V0i]);
f->V(2)=&(in.vert[V1i]);
ndone++;
}
if(V1i>=0 && V3i>=0 && V2i>=0 )
{
typename MeshType::FaceIterator f= Allocator<MeshType>::AddFaces(in,1);
f->V(0)=&(in.vert[V1i]);
f->V(1)=&(in.vert[V3i]);
f->V(2)=&(in.vert[V2i]);
ndone++;
}
}
}
}
template <class MeshType>
void Cylinder(int slices, int stacks, MeshType & m){
typename MeshType::VertexIterator vi = vcg::tri::Allocator<MeshType>::AddVertices(m,slices*(stacks+1));
for ( int i = 0; i < stacks+1; ++i)
for ( int j = 0; j < slices; ++j)
{
float x,y,h;
x = cos( 2.0 * M_PI / slices * j);
y = sin( 2.0 * M_PI / slices * j);
h = 2 * i / (float)(stacks) - 1;
(*vi).P() = typename MeshType::CoordType(x,h,y);
++vi;
}
typename MeshType::FaceIterator fi ;
for ( int j = 0; j < stacks; ++j)
for ( int i = 0; i < slices; ++i)
{
int a,b,c,d;
a = (j+0)*slices + i;
b = (j+1)*slices + i;
c = (j+1)*slices + (i+1)%slices;
d = (j+0)*slices + (i+1)%slices;
if(((i+j)%2) == 0){
fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1);
(*fi).V(0) = &m.vert[ a ];
(*fi).V(1) = &m.vert[ b ];
(*fi).V(2) = &m.vert[ c ];
fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1);
(*fi).V(0) = &m.vert[ c ];
(*fi).V(1) = &m.vert[ d ];
(*fi).V(2) = &m.vert[ a ];
}
else{
fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1);
(*fi).V(0) = &m.vert[ b ];
(*fi).V(1) = &m.vert[ c ];
(*fi).V(2) = &m.vert[ d ];
fi = vcg::tri::Allocator<MeshType>::AddFaces(m,1);
(*fi).V(0) = &m.vert[ d ];
(*fi).V(1) = &m.vert[ a ];
(*fi).V(2) = &m.vert[ b ];
}
}
if (m.HasPerFaceFlags()) {
for (typename MeshType::FaceIterator fi=m.face.begin(); fi!=m.face.end(); fi++) {
(*fi).SetF(2);
}
}
}
template <class MeshType>
void GenerateCameraMesh(MeshType &in){
typedef typename MeshType::CoordType MV;
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.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 ,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.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), };
int ff[88][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},
{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},
{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},
{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},
{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},
};
in.Clear();
Allocator<MeshType>::AddVertices(in,52);
Allocator<MeshType>::AddFaces(in,88);
in.vn=52;in.fn=88;
int i,j;
for(i=0;i<in.vn;i++)
in.vert[i].P()=vv[i];;
std::vector<typename MeshType::VertexPointer> index(in.vn);
typename MeshType::VertexIterator vi;
for(j=0,vi=in.vert.begin();j<in.vn;++j,++vi) index[j] = &*vi;
for(j=0;j<in.fn;++j)
{
in.face[j].V(0)=index[ff[j][0]];
in.face[j].V(1)=index[ff[j][1]];
in.face[j].V(2)=index[ff[j][2]];
}
}
//@}
} // End Namespace TriMesh
} // End Namespace vcg
#endif

102
vcg/complex/trimesh/create/readme.txt
View File

@ -1,102 +0,0 @@
MARCHING CUBES & EXTENDED MARCHING CUBES
===================================================================================
In breve le classi coinvolte sono 3 e sono:
* MerchingCubes ed ExtendedMarchingCubes
processano una cella alla volta, aggiungendo per ogni chiamata a ProcessCell
l'insieme di triangoli approssimante la superficie che interseca la cella
* Walker
gestisce l'attraversamento del volume, servendo le chiamate effettuate dagli
algoritmi di estrazione della superficie al volume e cachandone i risultato
* Volume
conosce come calcolare il campo scalare all'interno del volume da processare
e come calcolare le intersezioni superficie/segmenti.
DESCRIZIONE
====================================================================================
Le classi che implementano gli algoritmi MarchingCubes ed ExtendedMarchingCubes
sono state implementate così da risultare quanto più generiche possibile: ogni chiamata
al metodo ProcessCell(Point3i p1, Point3i p2) analizza la cella del volume individuata
dai due punti p1 e p2 e l'analisi di questa cella si conclude esattamente al ritorno da questa
chiamata: nel caso infatti la superficie da estrarre attraversi questa cella, all'interno
di questa stessa chiamata la mesh viene aggiornata con un opportuno insieme di triangoli.
L'assunzione alla base di questa astrazione è l'esistenza di altri due entità, il Walker
ed il Volume; sebbene sulla loro implementazione è lasciata la più completa libertà, è utile
precisare in quale relazione essi stiano rispetto agli algoritmi di estrazione di superfici.
Un esempio che riassume quanto qui esposto è incluso nella libreria: vd. vcg/apps/test/extractors.
VOLUME
====================================================================================
Gli algoritmi di estrazione di superfici risalgono alla superficie utilizzando i valori
di un campo scalare definito sul volume da processare campionato sui vertici di una qualche
griglia. Questo campo scalare sarà generalmente diverso a seconda del tipo di applicazione.
Il Volume è appunto quella classe che racchiude il campo scalare e di cui ne conosce le proprietà.
In realtà, all'interno dell'algoritmo di estrazione di superfici, non esiste alcun collegamento
esplicito con il Volume: tutte le sue chiamate sono rivolte al Walker. Questo perché
(per motivi che saranno esposti successivamente a proposito del Walker) il Walker potrebbe già
possedere il valore del campo calcolato in un dato punto, che evita così di farselo ricalcolare
nuovamente dal Volume: Similmente, quando viene chiesto di calcolare il punto di intersezione
della superficie con un segmento, se il Walker dispone già di questa informazione, la restituisce
all'algoritmo di estrazione di superfici, altrimenti il calcolo verrà effettuato dal Volume: il
punto ottenuto verrà restituito al Walker che potrà quindi soddisfare la richiesta iniziale.
Il motivo per cui si è scelto di frapporre un Walker fra gli algoritmi di estrazione di superfici
ed il Volume è esclusivamente di ottimizzazione. L'idea di fondo è che il Walker è quell'oggetto
attrverso cui avviene la visita del volume: per ogni cella del volume, esso effettua la chiamata
ProcessCell. Conoscendo l'ordine di visita, il Walker è anche la classe candidata ad implementare
le politiche di caching, in quanto sa esattamente da che momento può essere utile una certa
informazione e per quanto a lungo può essere conveniente mantenerla prima di liberarsene.
WALKER
====================================================================================
Poiché la politica di visita del volume è realizzata all'interno del walker,
è opportuno che sempre all'interno del walker vengano realizzate le politiche
di caching rivolte ad ottimizzare l'esecuzione degli algoritmi MC ed EMC. Durante
il processing di ogni cella questi algoritmi possono chiamare le seguenti funzioni
del Walker:
MC EMC
------------------------------------------
V(i, j, k) X X
GetXIntercept(p1, p2, v) X X
GetYIntercept(p1, p2, v) X X
GetZIntercept(p1, p2, v) X X
Exist(p1, p2, v) X
const float V(int i, int j, int k) const
La superficie che attraversa ogni cella viene ricavata dall'algoritmo di estrazione
analizzando il valore del campo sugli otto spigoli di ogni voxel del volume;
per ogni voxel, il valore del campo sui suoi otto spigoli vengono richiesti
dall'algoritmo di estrazione al walker: se questo valore è già stato calcolato
e cachato, il walker restituisce direttamente tale valore; altrimenti il valore
del campo in questo spigolo viene calcolato (eventualmente cachato) e restituito
al walker. In questo modo il valoro del campo ad ogni punto viene calcolato una
sola volta anziché 8, questo puo' essere molto utile nel caso si utilizzi dataset
volumetrici mantenuti implicitamente.
void GetXIntercept(Point3i p1, Point3i p2, VertexPointer v)
void GetYIntercept(Point3i p1, Point3i p2, VertexPointer v)
void GetZIntercept(Point3i p1, Point3i p2, VertexPointer v)
Dall'analisi del valore del campo agli spigoli di un dato voxel, l'algoritmo di
estrazione ha rilevato che la superficie interseca lo spigolo avente estremi p1
e p2(a seconda dell'orientazione di questo spigolo, viene chiamato uno dei tre
metodi): al termine di una di queste chiamate, v deve puntare al vertice della
mesh (di coordinate comprese fra p1 e p2) attraverso cui passa la superficie.
Se questo vertice è stato già calcolato ed inserito nella mesh, il walker deve
risalire a tale vertice e memorizzarne in v il suo puntatore. Altrimenti deve provvedere
ad aggiugnere un nuovo vertice ed a calcolare la sua posizione; v deve puntare
in questo caso al vertice appena inserito.
Il motivo per cui questo calcolo non viene implementato direttamente negli algoritmi
di estrazione (possibile per es. attraverso interpolazione lineare del valore
del campo nei punti p1 e p2) è che questo calcolo può essere fatto in maniere
molto più precisa conoscendo come il campo viene calcolato, essendo infatti ciò
dipendente dall'applicazione.
bool Exist(Point3i p1, Point3i p2, VertexPointer v)
Questo metodo viene chiamato solamente all'interno dell'algoritmo MarchingCubes.
A differenza dei tre motodi precedenti, in questo caso si vuole sapere solamente
se esiste già un vertice tra i punti p1 e p2: nel caso tale vertice esista, Exist
deve resituire true ed v deve puntare a tale vertice; se invece tale vertice non
esiste, Exist deve restituire false e v deve prendere il valore NULL.
NB: nel caso in cui il vertice non esiste, alla mesh non deve essere
aggiunto alcun nuovo vertice.

629
vcg/complex/trimesh/create/resampler.h
View File

@ -1,629 +0,0 @@
/****************************************************************************
* VCGLib o o *
* Visual and Computer Graphics Library o o *
* _ O _ *
* Copyright(C) 2004 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
#ifndef __VCG_MESH_RESAMPLER
#define __VCG_MESH_RESAMPLER
#include <vcg/complex/trimesh/update/normal.h>
#include <vcg/complex/trimesh/update/flag.h>
#include <vcg/complex/trimesh/update/bounding.h>
#include <vcg/complex/trimesh/update/edges.h>
#include <vcg/complex/trimesh/create/marching_cubes.h>
#include <vcg/space/index/grid_static_ptr.h>
#include <vcg/complex/trimesh/closest.h>
#include <vcg/space/box3.h>
namespace vcg {
namespace tri {
/** \addtogroup trimesh */
/*@{*/
/*@{*/
/** Class Resampler.
This is class reasmpling a mesh using marching cubes methods
@param OLD_MESH_TYPE (Template Parameter) Specifies the type of mesh to be resampled
@param NEW_MESH_TYPE (Template Parameter) Specifies the type of output mesh.
*/
template <class OLD_MESH_TYPE,class NEW_MESH_TYPE, class FLT, class DISTFUNCTOR = vcg::face::PointDistanceBaseFunctor<typename OLD_MESH_TYPE::ScalarType > >
class Resampler : public BasicGrid<FLT>
{
typedef OLD_MESH_TYPE Old_Mesh;
typedef NEW_MESH_TYPE New_Mesh;
//template <class OLD_MESH_TYPE,class NEW_MESH_TYPE>
class Walker : BasicGrid<float>
{
private:
typedef int VertexIndex;
typedef OLD_MESH_TYPE Old_Mesh;
typedef NEW_MESH_TYPE New_Mesh;
typedef typename New_Mesh::CoordType NewCoordType;
typedef typename New_Mesh::VertexType* VertexPointer;
typedef typename Old_Mesh::FaceContainer FaceCont;
typedef typename vcg::GridStaticPtr<typename Old_Mesh::FaceType> GridType;
protected:
int SliceSize;
int CurrentSlice;
typedef tri::FaceTmark<Old_Mesh> MarkerFace;
MarkerFace markerFunctor;
VertexIndex *_x_cs; // indici dell'intersezioni della superficie lungo gli Xedge della fetta corrente
VertexIndex *_y_cs; // indici dell'intersezioni della superficie lungo gli Yedge della fetta corrente
VertexIndex *_z_cs; // indici dell'intersezioni della superficie lungo gli Zedge della fetta corrente
VertexIndex *_x_ns; // indici dell'intersezioni della superficie lungo gli Xedge della prossima fetta
VertexIndex *_z_ns; // indici dell'intersezioni della superficie lungo gli Zedge della prossima fetta
//float *_v_cs;///values of distance fields for each direction in current slice
//float *_v_ns;///values of distance fields for each direction in next slice
typedef typename std::pair<bool,float> field_value;
field_value* _v_cs;
field_value* _v_ns;
New_Mesh *_newM;
Old_Mesh *_oldM;
GridType _g;
public:
float max_dim; // the limit value of the search (that takes into account of the offset)
float offset; // an offset value that is always added to the returned value. Useful for extrarting isosurface at a different threshold
bool DiscretizeFlag; // if the extracted surface should be discretized or not.
bool MultiSampleFlag;
bool AbsDistFlag; // if true the Distance Field computed is no more a signed one.
Walker(const Box3f &_bbox, Point3i _siz )
{
this->bbox= _bbox;
this->siz=_siz;
ComputeDimAndVoxel();
SliceSize = (this->siz.X()+1)*(this->siz.Z()+1);
CurrentSlice = 0;
offset=0;
DiscretizeFlag=false;
MultiSampleFlag=false;
AbsDistFlag=false;
_x_cs = new VertexIndex[ SliceSize ];
_y_cs = new VertexIndex[ SliceSize ];
_z_cs = new VertexIndex[ SliceSize ];
_x_ns = new VertexIndex[ SliceSize ];
_z_ns = new VertexIndex[ SliceSize ];
_v_cs= new field_value[(this->siz.X()+1)*(this->siz.Z()+1)];
_v_ns= new field_value[(this->siz.X()+1)*(this->siz.Z()+1)];
};
~Walker()
{}
float V(const Point3i &p)
{
return V(p.V(0),p.V(1),p.V(2));
}
std::pair<bool,float> VV(int x,int y,int z)
{
assert ((y==CurrentSlice)||(y==(CurrentSlice+1)));
//test if it is outside the bb of the mesh
//vcg::Point3f test=vcg::Point3f((float)x,(float)y,(float)z);
/*if (!_oldM->bbox.IsIn(test))
return (1.f);*/
int index=GetSliceIndex(x,z);
if (y==CurrentSlice) return _v_cs[index];
else return _v_ns[index];
}
float V(int x,int y,int z)
{
if(DiscretizeFlag) return VV(x,y,z).second+offset<0?-1:1;
return VV(x,y,z).second+offset;
}
///return true if the distance form the mesh is less than maxdim and return distance
field_value DistanceFromMesh(Point3f &pp,Old_Mesh */*mesh*/)
{
float dist;
typename Old_Mesh::FaceType *f=NULL;
const float max_dist = max_dim;
vcg::Point3f testPt;
this->IPfToPf(pp,testPt);
vcg::Point3f closestNormV,closestNormF;
vcg::Point3f closestPt;
vcg::Point3f pip(-1,-1,-1);
// Note that PointDistanceBaseFunctor does not require the edge and plane precomptued.
// while the PointDistanceFunctor requires them.
DISTFUNCTOR PDistFunct;
f = _g.GetClosest(PDistFunct,markerFunctor,testPt,max_dist,dist,closestPt);
if (f==NULL) return field_value(false,0);
if(AbsDistFlag) return field_value(true,dist);
assert(!f->IsD());
bool retIP;
// To compute the interpolated normal we use the more robust function that require to know what is the most orhogonal direction of the face.
if((*f).Flags() & Old_Mesh::FaceType::NORMX) retIP=InterpolationParameters(*f,0,closestPt, pip);
else if((*f).Flags() & Old_Mesh::FaceType::NORMY) retIP=InterpolationParameters(*f,1,closestPt, pip);
else if((*f).Flags() & Old_Mesh::FaceType::NORMZ) retIP=InterpolationParameters(*f,2,closestPt, pip);
else assert(0);
assert(retIP); // this should happen only if the starting mesh has degenerate faces.
const float InterpolationEpsilon = 0.00001f;
int zeroCnt=0;
if(pip[0]<InterpolationEpsilon) ++zeroCnt;
if(pip[1]<InterpolationEpsilon) ++zeroCnt;
if(pip[2]<InterpolationEpsilon) ++zeroCnt;
assert(zeroCnt<3);
Point3f dir=(testPt-closestPt).Normalize();
// Note that the two signs could be discordant.
// Always choose the best one according to where the nearest point falls.
float signBest;
// Compute test if the point see the surface normal from inside or outside
// Surface normal for improved robustness is computed both by face and interpolated from vertices.
if(zeroCnt>0) // we Not are in the middle of the face so the face normal is NOT reliable.
{
closestNormV = (f->V(0)->cN())*pip[0] + (f->V(1)->cN())*pip[1] + (f->V(2)->cN())*pip[2] ;
signBest = dir.dot(closestNormV) ;
}
else
{
closestNormF = f->cN() ;
signBest = dir.dot(closestNormF) ;
}
if(signBest<0) dist=-dist;
return field_value(true,dist);
}
field_value MultiDistanceFromMesh(Point3f &pp, Old_Mesh */*mesh*/)
{
float distSum=0;
int positiveCnt=0; // positive results counter
const int MultiSample=7;
const Point3f delta[7]={Point3f(0,0,0),
Point3f( 0.2, -0.01, -0.02),
Point3f(-0.2, 0.01, 0.02),
Point3f( 0.01, 0.2, 0.01),
Point3f( 0.03, -0.2, -0.03),
Point3f(-0.02, -0.03, 0.2 ),
Point3f(-0.01, 0.01, -0.2 )};
for(int qq=0;qq<MultiSample;++qq)
{
Point3f pp2=pp+delta[qq];
field_value ff= DistanceFromMesh(pp2,_oldM);
if(ff.first==false) return field_value(false,0);
distSum += fabs(ff.second);
if(ff.second>0) positiveCnt ++;
}
if(positiveCnt<=MultiSample/2) distSum = -distSum;
return field_value(true, distSum/MultiSample);
}
/// compute the values if an entire slice (per y) distances>dig of a cell are signed with double of
/// the distance of the bb
void ComputeSliceValues(int slice,field_value *slice_values)
{
for (int i=0; i<=this->siz.X(); i++)
{
for (int k=0; k<=this->siz.Z(); k++)
{
int index=GetSliceIndex(i,k);
Point3f pp(i,slice,k);
if(this->MultiSampleFlag) slice_values[index] = MultiDistanceFromMesh(pp,_oldM);
else slice_values[index] = DistanceFromMesh(pp,_oldM);
}
}
//ComputeConsensus(slice,slice_values);
}
/*
For some reasons it can happens that the sign of the computed distance could not correct.
this function tries to correct these issues by flipping the isolated voxels with discordant sign
*/
void ComputeConsensus(int slice, field_value *slice_values)
{
float max_dist = min(min(this->voxel[0],this->voxel[1]),this->voxel[2]);
int flippedCnt=0;
int flippedTot=0;
int flippedTimes=0;
do
{
flippedCnt=0;
for (int i=0; i<=this->siz.X(); i++)
{
for (int k=0; k<=this->siz.Z(); k++)
{
int goodCnt=0;
int badCnt=0;
int index=GetSliceIndex(i,k);
int index_l,index_r,index_u,index_d;
if(slice_values[index].first)
{
float curVal= slice_values[index].second;
if(i > 0 ) index_l=GetSliceIndex(i-1,k); else index_l = index;
if(i < this->siz.X() ) index_r=GetSliceIndex(i+1,k); else index_r = index;
if(k > 0 ) index_d=GetSliceIndex(i,k-1); else index_d = index;
if(k < this->siz.Z() ) index_u=GetSliceIndex(i,k+1); else index_u = index;
if(slice_values[index_l].first) { goodCnt++; if(fabs(slice_values[index_l].second - curVal) > max_dist) badCnt++; }
if(slice_values[index_r].first) { goodCnt++; if(fabs(slice_values[index_r].second - curVal) > max_dist) badCnt++; }
if(slice_values[index_u].first) { goodCnt++; if(fabs(slice_values[index_u].second - curVal) > max_dist) badCnt++; }
if(slice_values[index_d].first) { goodCnt++; if(fabs(slice_values[index_d].second - curVal) > max_dist) badCnt++; }
if(badCnt >= goodCnt) {
slice_values[index].second *=-1.0f;
//slice_values[index].first = false;
flippedCnt++;
}
}
}
}
flippedTot+=flippedCnt;
flippedTimes++;
} while(flippedCnt>0);
#ifndef NO_QT
if(flippedTot>0)
qDebug("Flipped %i values in %i times",flippedTot,flippedTimes);
#endif
}
template<class EXTRACTOR_TYPE>
void ProcessSlice(EXTRACTOR_TYPE &extractor)
{
for (int i=0; i<this->siz.X(); i++)
{
for (int k=0; k<this->siz.Z(); k++)
{
bool goodCell=true;
Point3i p1(i,CurrentSlice,k);
Point3i p2=p1+Point3i(1,1,1);
for(int ii=0;ii<2;++ii)
for(int jj=0;jj<2;++jj)
for(int kk=0;kk<2;++kk)
goodCell &= VV(p1[0]+ii,p1[1]+jj,p1[2]+kk).first;
if(goodCell) extractor.ProcessCell(p1, p2);
}
}
}
template<class EXTRACTOR_TYPE>
void BuildMesh(Old_Mesh &old_mesh,New_Mesh &new_mesh,EXTRACTOR_TYPE &extractor,vcg::CallBackPos *cb)
{
_newM=&new_mesh;
_oldM=&old_mesh;
// the following two steps are required to be sure that the point-face distance without precomputed data works well.
tri::UpdateNormals<Old_Mesh>::PerFaceNormalized(old_mesh);
tri::UpdateNormals<Old_Mesh>::PerVertexAngleWeighted(old_mesh);
tri::UpdateFlags<Old_Mesh>::FaceProjection(old_mesh);
int _size=(int)old_mesh.fn*100;
_g.Set(_oldM->face.begin(),_oldM->face.end(),_size);
markerFunctor.SetMesh(&old_mesh);
_newM->Clear();
Begin();
extractor.Initialize();
for (int j=0; j<=this->siz.Y(); j++)
{
cb((100*j)/this->siz.Y(),"Marching ");
ProcessSlice<EXTRACTOR_TYPE>(extractor);//find cells where there is the isosurface and examine it
NextSlice();
}
extractor.Finalize();
typename New_Mesh::VertexIterator vi;
for(vi=new_mesh.vert.begin();vi!=new_mesh.vert.end();++vi)
if(!(*vi).IsD())
{
IPfToPf((*vi).cP(),(*vi).P());
}
}
//return the index of a vertex in slide as it was stored
int GetSliceIndex(int x,int z)
{
VertexIndex index = x+z*(this->siz.X()+1);
return (index);
}
//swap slices , the initial value of distance fields ids set as double of bbox of space
void NextSlice()
{
memset(_x_cs, -1, SliceSize*sizeof(VertexIndex));
memset(_y_cs, -1, SliceSize*sizeof(VertexIndex));
memset(_z_cs, -1, SliceSize*sizeof(VertexIndex));
std::swap(_x_cs, _x_ns);
std::swap(_z_cs, _z_ns);
std::swap(_v_cs, _v_ns);
CurrentSlice ++;
ComputeSliceValues(CurrentSlice + 1,_v_ns);
}
//initialize data strucures , the initial value of distance fields ids set as double of bbox of space
void Begin()
{
CurrentSlice = 0;
memset(_x_cs, -1, SliceSize*sizeof(VertexIndex));
memset(_y_cs, -1, SliceSize*sizeof(VertexIndex));
memset(_z_cs, -1, SliceSize*sizeof(VertexIndex));
memset(_x_ns, -1, SliceSize*sizeof(VertexIndex));
memset(_z_ns, -1, SliceSize*sizeof(VertexIndex));
ComputeSliceValues(CurrentSlice,_v_cs);
ComputeSliceValues(CurrentSlice+1,_v_ns);
}
bool Exist(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
{
int i = p1.X();// - _bbox.min.X())/_cell_size.X();
int z = p1.Z();// - _bbox.min.Z())/_cell_size.Z();
VertexIndex index = i+z*this->siz.X();
//VertexIndex index =GetSliceIndex(//
int v_ind = 0;
if (p1.X()!=p2.X()) //intersezione della superficie con un Xedge
{
if (p1.Y()==CurrentSlice)
{
if (_x_cs[index]!=-1)
{
v_ind = _x_cs[index];
v = &_newM->vert[v_ind];
assert(!v->IsD());
return true;
}
}
else
{
if (_x_ns[index]!=-1)
{
v_ind = _x_ns[index];
v = &_newM->vert[v_ind];
assert(!v->IsD());
return true;
}
}
v = NULL;
return false;
}
else if (p1.Y()!=p2.Y()) //intersezione della superficie con un Yedge
{
if (_y_cs[index]!=-1)
{
v_ind =_y_cs[index];
v = &_newM->vert[v_ind];
assert(!v->IsD());
return true;
}
else
{
v = NULL;
return false;
}
}
else if (p1.Z()!=p2.Z())
//intersezione della superficie con un Zedge
{
if (p1.Y()==CurrentSlice)
{
if ( _z_cs[index]!=-1)
{
v_ind = _z_cs[index];
v = &_newM->vert[v_ind];
assert(!v->IsD());
return true;
}
}
else
{
if (_z_ns[index]!=-1)
{
v_ind = _z_ns[index];
v = &_newM->vert[v_ind];
assert(!v->IsD());
return true;
}
}
v = NULL;
return false;
}
assert (0);
return false;
}
///interpolate
NewCoordType Interpolate(const vcg::Point3i &p1, const vcg::Point3i &p2,int dir)
{
float f1 = (float)V(p1);
float f2 = (float)V(p2);
float u = (float) f1/(f1-f2);
NewCoordType ret=vcg::Point3f((float)p1.V(0),(float)p1.V(1),(float)p1.V(2));
ret.V(dir) = (float) p1.V(dir)*(1.f-u) + u*(float)p2.V(dir);
return (ret);
}
///if there is a vertex in z axis of a cell return the vertex or create it
void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
{
assert(p1.X()+1 == p2.X());
assert(p1.Y() == p2.Y());
assert(p1.Z() == p2.Z());
int i = p1.X();// (p1.X() - _bbox.min.X())/_cell_size.X();
int z = p1.Z();//(p1.Z() - _bbox.min.Z())/_cell_size.Z();
VertexIndex index = i+z*this->siz.X();
VertexIndex pos=-1;
if (p1.Y()==CurrentSlice)
{
if ((pos=_x_cs[index])==-1)
{
_x_cs[index] = (VertexIndex) _newM->vert.size();
pos = _x_cs[index];
Allocator<New_Mesh>::AddVertices( *_newM, 1 );
v = &_newM->vert[pos];
v->P()=Interpolate(p1,p2,0);
return;
}
}
if (p1.Y()==CurrentSlice+1)
{
if ((pos=_x_ns[index])==-1)
{
_x_ns[index] = (VertexIndex) _newM->vert.size();
pos = _x_ns[index];
Allocator<New_Mesh>::AddVertices( *_newM, 1 );
v = &_newM->vert[pos];
v->P()=Interpolate(p1,p2,0);
return;
}
}
assert(pos>=0);
v = &_newM->vert[pos];
}
///if there is a vertex in y axis of a cell return the vertex or create it
void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
{
assert(p1.X() == p2.X());
assert(p1.Y()+1 == p2.Y());
assert(p1.Z() == p2.Z());
int i = p1.X(); // (p1.X() - _bbox.min.X())/_cell_size.X();
int z = p1.Z(); // (p1.Z() - _bbox.min.Z())/_cell_size.Z();
VertexIndex index = i+z*this->siz.X();
VertexIndex pos=-1;
if ((pos=_y_cs[index])==-1)
{
_y_cs[index] = (VertexIndex) _newM->vert.size();
pos = _y_cs[index];
Allocator<New_Mesh>::AddVertices( *_newM, 1);
v = &_newM->vert[ pos ];
v->P()=Interpolate(p1,p2,1);
}
assert(pos>=0);
v = &_newM->vert[pos];
}
///if there is a vertex in z axis of a cell return the vertex or create it
void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
{
assert(p1.X() == p2.X());
assert(p1.Y() == p2.Y());
assert(p1.Z()+1 == p2.Z());
int i = p1.X(); //(p1.X() - _bbox.min.X())/_cell_size.X();
int z = p1.Z(); //(p1.Z() - _bbox.min.Z())/_cell_size.Z();
VertexIndex index = i+z*this->siz.X();
VertexIndex pos=-1;
if (p1.Y()==CurrentSlice)
{
if ((pos=_z_cs[index])==-1)
{
_z_cs[index] = (VertexIndex) _newM->vert.size();
pos = _z_cs[index];
Allocator<New_Mesh>::AddVertices( *_newM, 1 );
v = &_newM->vert[pos];
v->P()=Interpolate(p1,p2,2);
return;
}
}
if (p1.Y()==CurrentSlice+1)
{
if ((pos=_z_ns[index])==-1)
{
_z_ns[index] = (VertexIndex) _newM->vert.size();
pos = _z_ns[index];
Allocator<New_Mesh>::AddVertices( *_newM, 1 );
v = &_newM->vert[pos];
v->P()=Interpolate(p1,p2,2);
return;
}
}
assert(pos>=0);
v = &_newM->vert[pos];
}
};//end class walker
public:
typedef Walker /*< Old_Mesh,New_Mesh>*/ MyWalker;
typedef vcg::tri::MarchingCubes<New_Mesh, MyWalker> MyMarchingCubes;
///resample the mesh using marching cube algorithm ,the accuracy is the dimension of one cell the parameter
static void Resample(Old_Mesh &old_mesh,New_Mesh &new_mesh, Box3f volumeBox, vcg::Point3<int> accuracy,float max_dist, float thr=0, bool DiscretizeFlag=false, bool MultiSampleFlag=false, bool AbsDistFlag=false, vcg::CallBackPos *cb=0 )
{
///be sure that the bounding box is updated
vcg::tri::UpdateBounding<Old_Mesh>::Box(old_mesh);
MyWalker walker(volumeBox,accuracy);
walker.max_dim=max_dist+fabs(thr);
walker.offset = - thr;
walker.DiscretizeFlag = DiscretizeFlag;
walker.MultiSampleFlag = MultiSampleFlag;
walker.AbsDistFlag = AbsDistFlag;
MyMarchingCubes mc(new_mesh, walker);
walker.BuildMesh(old_mesh,new_mesh,mc,cb);
}
};//end class resampler
};//end namespace tri
};//end namespace vcg
#endif