From 007d53b7e0a9df667b5f80c258aa480ea0609f01 Mon Sep 17 00:00:00 2001 From: ponchio Date: Wed, 30 May 2007 15:09:58 +0000 Subject: [PATCH] *** empty log message *** --- vcg/complex/trimesh/create/advancing_front.h | 499 ++++++++ vcg/complex/trimesh/create/ball_pivoting.h | 1172 +++++------------- 2 files changed, 813 insertions(+), 858 deletions(-) create mode 100644 vcg/complex/trimesh/create/advancing_front.h diff --git a/vcg/complex/trimesh/create/advancing_front.h b/vcg/complex/trimesh/create/advancing_front.h new file mode 100644 index 00000000..9efa59f7 --- /dev/null +++ b/vcg/complex/trimesh/create/advancing_front.h @@ -0,0 +1,499 @@ +#ifndef MLS_ADVANCE_H +#define MLS_ADVANCE_H + +#include +#include +#include +#include +#include + +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::iterator next; + std::list::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); + } +}; + +template class AdvancingFront { + public: + + typedef typename MESH::VertexType VertexType; + typedef typename MESH::FaceType FaceType; + typedef typename MESH::ScalarType ScalarType; + typedef typename MESH::VertexType::CoordType Point3x; + + + +// protected: + std::list front; + std::list deads; + std::vector 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) { +// UpdateTopology::VertexFace(mesh); + //UpdateTopology::TestVertexFace(mesh); //odd i would like to return a false not an assert... + UpdateFlags::VertexClear(mesh); +// UpdateFlags::FaceBorderFromVF(mesh); + UpdateFlags::FaceBorderFromNone(mesh); + UpdateFlags::VertexBorderFromFace(mesh); + + nb.clear(); + nb.resize(mesh.vert.size(), 0); + + CreateLoops(); + } + virtual ~AdvancingFront() {} + + void BuildMesh(CallBackPos call = NULL, int interval = 512) { + while(1) { + if(call) call(0, "Advancing front"); + for(int i = 0; i < interval; i++) + if(!AddFace()) return; + } + } + +protected: + //Implement these functions in your subclass + virtual bool Seed(int &v0, int &v1, int &v2) = 0; + virtual int Place(FrontEdge &e, std::list::iterator &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; + // 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::iterator s = front.begin(); s != front.end(); s++) { + (*s).previous = front.end(); + (*s).next = front.end(); + } + //now create loops: + for(std::list::iterator s = front.begin(); s != front.end(); s++) { + for(std::list::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; + + } + } + } + + 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::iterator e = front.end(); + std::list::iterator last = e; + std::list::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 SeedFace(); + } + + std::list::iterator ei = front.begin(); + FrontEdge ¤t = *ei; + FrontEdge &previous = *current.previous; + FrontEdge &next = *current.next; + + int v0 = current.v0, v1 = current.v1; + assert(nb[v0] < 10 && nb[v1] < 10); + + std::list::iterator touch = front.end(); + + int v2 = Place(current, touch); + if(v2 == -1) { + KillEdge(ei); + return true; + } + + assert(v2 != v0 && v2 != v1); + + if(touch != front.end()) { + + //check for orientation and manifoldness + if(!CheckEdge(v0, v2) || !CheckEdge(v2, v1)) { + KillEdge(ei); + return 0; + } + //touch == current.previous? + if(v2 == previous.v0) { + + /*touching previous FrontEdge (we reuse previous) + next + ------->v2 -----> v1------> + \ / + \ / + previous \ / current + \ / + v0 */ + + Detach(v0); + + std::list::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) { + + + /*touching next FrontEdge (we reuse next) + previous + ------->v0 -----> v2------> + \ / + \ / + \ / next + \ / + v1 */ + + Detach(v1); + std::list::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 { + + //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::iterator left = touch; + std::list::iterator right = (*touch).previous; + + //this would be a really bad join + if(v1 == (*right).v0 || v0 == (*left).v1) { + KillEdge(ei); + return 0; + } + + nb[v2]++; + + std::list::iterator down = NewEdge(FrontEdge(v2, v1, v0, mesh.face.size())); + std::list::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 { + + /* 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::iterator down = NewEdge(FrontEdge(v2, v1, v0, mesh.face.size())); + std::list::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 1; + } + +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; + // 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; + } + //front management: + + //Add a new FrontEdge to the back of the queue + std::list::iterator NewEdge(FrontEdge e) { + return front.insert(front.end(), e); + } + + //move an Edge among the dead ones + void KillEdge(std::list::iterator e) { + (*e).active = false; + deads.splice(deads.end(), front, e); + } + + void Erase(std::list::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::iterator e) { + front.splice(front.end(), front, e); + } + + void MoveFront(std::list::iterator e) { + front.splice(front.begin(), front, e); + } + + //check if e can be sewed with one of oits neighbours + bool Glue(std::list::iterator e) { + return Glue((*e).previous, e) || Glue(e, (*e).next); + } + + //Glue toghether a and b (where a.next = b + bool Glue(std::list::iterator a, std::list::iterator b) { + if((*a).v0 != (*b).v1) return false; + + std::list::iterator previous = (*a).previous; + std::list::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 AdvancingTest: public AdvancingFront { + 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) {} + + 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, std::list::iterator &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(list::iterator k = this->front.begin(); k != this->front.end(); k++) + if((*k).v0 == i) touch = k; + for(list::iterator k = this->deads.begin(); k != this->deads.end(); k++) + if((*k).v0 == i) touch = k; + break; + } + } + if(vn == this->mesh.vert.size()) { + VertexType v; + v.P() = point; + v.ClearFlags(); + AddVertex(v); + } + return vn; + } +}; + +}//namespace tri +}//namespace vcg + +#endif diff --git a/vcg/complex/trimesh/create/ball_pivoting.h b/vcg/complex/trimesh/create/ball_pivoting.h index 4fd7971e..fdbf6d28 100644 --- a/vcg/complex/trimesh/create/ball_pivoting.h +++ b/vcg/complex/trimesh/create/ball_pivoting.h @@ -1,905 +1,361 @@ -#ifndef VCG_PIVOT_H -#define VCG_PIVOT_H +#ifndef BALL_PIVOTING_H +#define BALL_PIVOTING_H -#include -#include -#include - -#include "vcg/space/index/grid_static_ptr.h" -#include "vcg/complex/trimesh/closest.h" -#include -using namespace std; +#include "advancing_front.h" +#include +#include +/* 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 { - - + namespace tri { +template class BallPivoting: public AdvancingFront { + public: + typedef typename MESH::VertexType VertexType; + typedef typename MESH::FaceType FaceType; + typedef typename MESH::ScalarType ScalarType; + typedef typename MESH::VertexType::CoordType Point3x; + typedef GridStaticPtr 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: + BallPivoting(MESH &_mesh, float _radius = 0, + float minr = 0.2, float angle = M_PI/2): + + AdvancingFront(_mesh), radius(_radius), + min_edge(minr), max_edge(1.8), max_angle(cos(angle)), + last_seed(-1) { + + //compute bbox + for(int i = 0; i < (int)this->mesh.vert.size(); i++) + this->mesh.bbox.Add(this->mesh.vert[i].P()); -template -class Pivot { - public: - typedef GridStaticPtr StaticGrid; - typedef typename MESH::VertexType CVertex; - typedef typename MESH::FaceType CFace; - typedef typename MESH::ScalarType ScalarType; - typedef typename CVertex::CoordType Point3x; - - ScalarType radius; //default 1 (not meaningful - ScalarType mindist; //minimum distance between points in the mesh (% of radius) - ScalarType crease; // -0.5 - bool normals; //default false - Box3 box; + assert(this->mesh.vn > 3); + if(radius == 0) + radius = sqrt((this->mesh.bbox.Diag()*this->mesh.bbox.Diag())/this->mesh.vn); + else + radius *= this->mesh.bbox.Diag()/100; - MESH &mesh; - StaticGrid grid; - - template class Edge { - public: - int v0, v1, v2; //v0, v1 represent the edge, v2 the other vertex in the face - //this edge belongs to - int face; //corresponding face - Coord center; //center of the sphere touching the face - int count; //test delay touch edges. + min_edge *= radius; + max_edge *= radius; - bool active; //keep tracks of wether it is in front or in deads - float angle; - int candidate; - Coord newcenter; - - //the loops in the front are mantained as a double linked list - typename std::list::iterator next; - typename std::list::iterator previous; + //enlarging the bbox for out-of-box queries + this->mesh.bbox.Offset(4*radius); + grid.Set(this->mesh.vert.begin(), this->mesh.vert.end(), this->mesh.bbox); - Edge() {} - Edge(int _v0, int _v1, int _v2, int _face, Point3f &_center): - v0(_v0), v1(_v1), v2(_v2), - face(_face), center(_center), count(-1), active(true) { - assert(v0 != v1 && v1 != v2 && v0 != v2); - } - }; - typedef Edge Edgex; + //mark visited points + std::vector targets; + std::vector points; + std::vector dists; - /* front edges of the mesh: - to expand the front we get the first edge - if an edge cannot create a new triangle it is marked dead and moved - to the end of the list - the new edges are inserted to the back (before dead_begin) */ - - std::list front; - std::list deads; - std::vector nb; //number of fronts a vertex is into, - //this is used for the Visited and Border flags - //but adding topology may not be needed anymode - int last_seed; + usedBit = VertexType::NewBitFlag(); -// Parameters: -// a standard trimesh. Vertex is required to have normals. -// If some faces are already in the process start from the border edges -// Bounding Box must have been already computed. -// If faces are in border flags must be already computed -// is the radius of the pivoting ball. If set to 0 an automatic guess is used ( sqrt(BBox_Diag^2 / vn) ) -// is the minimal distance between two points -// in the -// -Pivot(MESH &_mesh, ScalarType _radius, ScalarType _mindist = 0.1, ScalarType _crease = -0.5): - mesh(_mesh), radius(_radius), mindist(_mindist), crease(_crease), normals(true), last_seed(0) { - box = mesh.bbox; - //estimate radius if not provided - if(radius <= 0.0f) - radius = sqrt((box.Diag()*box.Diag())/mesh.vn); - - UpdateFlags::VertexClearV(mesh); - - /* we need to enlarge the grid to allow queries from little outside of the box - Someone is a bit lazy... */ - box.Offset(4*radius); - grid.Set(mesh.vert.begin(), mesh.vert.end(), box); - nb.clear(); - nb.resize(mesh.vert.size(), 0); - - if(mesh.face.size()) { - //init border from mesh - Point3x center; - CVertex *start = &*mesh.vert.begin(); - for(int i = 0; i < mesh.face.size(); i++) { - CFace &face = mesh.face[i]; - for(int k = 0; k < 3; k++) { - if(!face.V(k)->IsB()) face.V(k)->SetV(); - cluster(face.V(k) - start); - if(face.IsB(k)) { - //compute center: - findSphere(face.P(k), face.P((k+1)%3), face.P((k+2)%3), center); - newEdge(Edgex(face.V0(k) - start, face.V1(k) - start, face.V2(k) - start, - i, center)); - } - } + 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 = trimesh::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)); } - for(typename std::list::iterator s = front.begin(); s != front.end(); s++) { - (*s).previous = front.end(); - (*s).next = front.end(); - printf("%d %d\n", (*s).v0, (*s).v1); - } - //now create loops: - for(typename std::list::iterator s = front.begin(); s != front.end(); s++) { - for(typename std::list::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; - - } - } -/* for(typename std::list::iterator s = front.begin(); s != front.end(); s++) { - assert((*s).next != front.end()); - assert((*s).previous != front.end()); - } */ -/* for(int i = 0; i < mesh.face.size(); i++) { - CFace &face = mesh.face[i]; - for(int k = 0; k < 3; k++) - face.V(k) = (CVertex *)(face.V(k) - start); - }*/ - - } else { - for(int i = 0; i < mesh.vert.size(); i++) { - mesh.vert[i].ClearFlags(); - } - - } - srand(time(NULL)); - } - - /* return false if you want to stop.\n */ -void buildMesh(CallBackPos *call = NULL, int interval = 512) { - bool done = false; - float estimated_faces = mesh.vn*2; - while(!done) { - //estimating progress - float vdeleted = mesh.vert.size() - mesh.vn; - float vused = mesh.face.size()/2.0f; - float progress = 100*(vdeleted + vused)/mesh.vert.size(); - if(progress > 99) progress = 99; - if(call) call((int)progress, "Pivoting"); - for(int i = 0; i < interval; i++) { - if(addFace() == -1) { - done = true; - break; - } - } - } - } - - /* select a vertex at random, a small group of nearby vertices and looks - for a sphere that touches 3 and contains none. - Use the center of the box to get a sphere inside (or outside) the model - You may be unlucky... */ - -bool seed(bool outside = true, int start = -1) { - - //pick a random point (well...) - if(start == -1) start = 0;//rand()%mesh.vert.size(); - - //get a sphere of neighbours - std::vector targets; - std::vector dists; - int n = getInSphere(mesh.vert[start].P(), 2*radius, targets, dists); - if(n < 3) { - mesh.vert[start].SetD(); - //bad luck. we should call seed again (assuming random pick) up to - //some maximum tries. im lazy. - return false; + assert(f.V(k)->IsUserBit(usedBit)); } + } + } + + ~BallPivoting() { + VertexType::DeleteBitFlag(usedBit); + } + + bool Seed(int &v0, int &v1, int &v2) { + bool use_normals = false; + //get a sphere of neighbours + std::vector targets; + std::vector points; + std::vector 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 = trimesh::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++) { - if(dists[i] < radius) { - int id = targets[i]; - CVertex &v = mesh.vert[id]; - if(v.IsB() || v.IsV()) { - mesh.vert[start].SetD(); - return false; - } + for(int i = 0; i < n; i++) { + VertexType &v = *(targets[i]); + if(v.IsV()) { + success = false; + break; } } - - int v0, v1, v2; - bool found = false; + 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++) { - v0 = targets[i]; - CVertex &vv0 = mesh.vert[v0]; - if(vv0.IsD() || vv0.IsB() || vv0.IsV()) continue; - Point3x &p0 = mesh.vert[v0].P(); - Point3x out = (p0 - box.Center()); - if(!outside) out = -out; + vv0 = targets[i]; + if(vv0->IsD()) continue; + Point3x &p0 = vv0->P(); for(int k = i+1; k < n; k++) { - v1 = targets[k]; - CVertex &vv1 = mesh.vert[v1]; - if(vv1.IsD() || vv1.IsB() || vv1.IsV()) continue; - Point3x &p1 = mesh.vert[v1].P(); - if((p1 - p0).Norm() < mindist*radius) continue; + 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++) { - v2 = targets[j]; - CVertex &vv2 = mesh.vert[v2]; - if(vv2.IsD() || vv2.IsB() || vv2.IsV()) continue; - Point3x &p2 = mesh.vert[v2].P(); - if((p2 - p0).Norm() < mindist*radius) continue; - if((p2 - p1).Norm() < mindist*radius) continue; + 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(!normals) { - //check normal pointing inside - if(normal * out < 0) continue; - } else { - if(normal * vv0.N() < 0) continue; - if(normal * vv1.N() < 0) continue; - if(normal * vv2.N() < 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; - bool failed = false; + if(!FindSphere(p0, p1, p2, center)) { + continue; + } + //check no other point inside - for(int t = 0; t < n; t++) { - Point3x &p = mesh.vert[targets[t]].P(); - if((center - p).Norm() <= radius) { - failed = true; - break; - } + int t; + for(t = 0; t < n; t++) { + if((center - targets[t]->P()).Norm() <= radius) + break; } - //check on the other side there are not a surface - Point3x recenter; - if(!findSphere(p0, p2, p1, recenter)) continue; - for(int t = 0; t < n; t++) { - CVertex &v = mesh.vert[targets[t]]; - Point3x &p = v.P(); - if((center - p).Norm() <= radius && (v.IsV() || v.IsB())) { - failed = true; - break; - } + if(t < n) { + continue; } - if(failed) continue; - found = true; + //check on the other side there is not a surface + Point3x opposite = center + normal*(((center - p0)*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(!found) { //see bad luck above - return false; + if(!success) { //see bad luck above + continue; } + Mark(vv0); + Mark(vv1); + Mark(vv2); - assert(!front.size()); - //TODO: should i check the edgex too? - addFace(v0, v1, v2); - - //create the border of the first face - typename std::list::iterator e = front.end(); - typename std::list::iterator last; - for(int i = 0; i < 3; i++) { - int v0 = mesh.face.back().V0(i) - &*mesh.vert.begin(); - int v1 = mesh.face.back().V1(i) - &*mesh.vert.begin(); - int v2 = mesh.face.back().V2(i) - &*mesh.vert.begin(); - nb[v0] = 1; - assert(!mesh.vert[v0].IsB()); - mesh.vert[v0].SetB(); - Edgex edge(v0, v1, v2, 0, center); - edge.previous = e; - e = front.insert(front.begin(), edge); - if(i == 0) last = e; - (*edge.previous).next = e; - - cluster(v0); - } - - //connect last and first - (*e).next = last; - (*last).previous = e; - return true; + 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, std::list::iterator &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)) return -1; + + Point3x start_pivot = center - middle; + Point3x axis = (v1 - v0); + + ScalarType axis_len = axis.SquaredNorm(); + if(axis_len > 4*radius*radius) return false; + 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 targets; + std::vector dists; + std::vector points; + + int n = trimesh::GetInSphereVertex(this->mesh, grid, middle, r + radius, targets, dists, points); + + if(targets.size() == 0) { + return false; //this really would be strange but one never knows. + } + + VertexType *candidate = NULL; + ScalarType min_angle = M_PI; + + for(int i = 0; i < 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)); - - /* expand the front adding 1 face. Return false on failure (id when - all edges are dead returns: - 1: added a face - 0: added nothing - -1: finished */ - -int addFace() { - - //We try to seed again -/* if(!mesh.face.size()) { - for(int i = 0; i < 100; i++) { - ++last_seed; - CVertex &v = mesh.vert[last_seed-1]; - if(v.IsD() || v.IsV() || v.IsB()) continue; - - printf("seeding new: %i\n", last_seed-1); - if(seed(true, last_seed-1)) return 1; - } - return -1; - }*/ - if(!front.size()) { - //maybe there are unconnected parts of the mesh: - //find a non D, V, B point and try to seed if failed D it. - while(last_seed < mesh.vert.size()) { - ++last_seed; - CVertex &v = mesh.vert[last_seed-1]; - if(v.IsD() || v.IsV() || v.IsB()) continue; - - printf("seeding new: %i\n", last_seed-1); - if(seed(true, last_seed-1)) return 1; - - v.SetD(); - --mesh.vn; - return 0; - } - printf("done\n"); - return -1; - } - if(last_seed > 1) printf("frontsixe: %d\n", front.size()); - typename std::list::iterator ei = front.begin(); - Edgex &e = *ei; - Edgex &previous = *e.previous; - Edgex &next = *e.next; - - int v0 = e.v0, v1 = e.v1; - - assert(nb[v0] < 10 && nb[v1] < 10); - int v2; - Point3x center; - bool success = pivot(e); - v2 = e.candidate; - center = e.newcenter; - - //if no pivoting or we are trying to connect to the inside of the mesh. - if(!success || mesh.vert[v2].IsV()) { - killEdge(ei); - return 0; - } - - //does v2 belongs to a front? (and which?) - typename std::list::iterator touch = touches(ei); - - assert(v2 != v0 && v2 != v1); - - int fn = mesh.face.size(); - if(touch != front.end()) { - - //check for orientation and manifoldness - if(!checkEdge(v0, v2) || !checkEdge(v2, v1)) { - killEdge(ei); - return 0; - } + if(v->IsUserBit(usedBit) && !(v->IsB())) continue; + if(id == edge.v0 || id == edge.v1 || id == edge.v2) continue; - if(v2 == previous.v0) { - - /*touching previous edge (we reuse previous) - next - ------->v2 -----> v1------> - \ / - \ / - previous \ / e - \ / - v0 */ - - detach(v0); - - typename std::list::iterator up = newEdge(Edgex(v2, v1, v0, fn, center)); - (*up).previous = previous.previous; - (*up).next = e.next; - (*previous.previous).next = up; - next.previous = up; - erase(e.previous); - erase(ei); - trovamiunnome(up); - - - } else if(v2 == next.v1) { - - - /*touching next edge (we reuse next) - previous - ------->v0 -----> v2------> - \ / - \ / - \ / next - \ / - v1 */ - - detach(v1); - - typename std::list::iterator up = newEdge(Edgex(v0, v2, v1, fn, center)); - (*up).previous = e.previous; - (*up).next = (*e.next).next; - previous.next = up; - (*next.next).previous = up; - erase(e.next); - erase(ei); - trovamiunnome(up); - - - } else { - - //touching some loop: split (or merge it is local does not matter. - //like this - /* - left right - <--------v2-<------ - /|\ - / \ - up / \ down - / \ - / V - ----v0 - - - > v1--------- - e */ - typename std::list::iterator left = touch; - typename std::list::iterator right = (*touch).previous; - typename std::list::iterator up = ei; - - //this would be a really bad join - if(v1 == (*right).v0 || v0 == (*left).v1) { - killEdge(ei); - return 0; - } - - nb[v2]++; - - typename std::list::iterator down = newEdge(Edgex(v2, v1, v0, fn, center)); - - (*right).next = down; - (*down).previous = right; - - (*down).next = e.next; - next.previous = down; - - (*left).previous = up; - (*up).next = left; - - (*up).v2 = v1; - (*up).v1 = v2; - (*up).face = fn; - (*up).center = center; - moveBack(ei); - } - - - - } else { - - /* adding a new vertex - - v2 - /|\ - / \ - up / \ down - / \ - / V - ----v0 - - - > v1--------- */ - assert(!mesh.vert[v2].IsB()); //fatal error! a new point is already a border? - - //clustering points aroundf v2 - cluster(v2); - - nb[v2]++; - mesh.vert[v2].SetB(); - typename std::list::iterator down = newEdge(Edgex(v2, v1, v0, fn, center)); - (*down).previous = ei; - (*down).next = e.next; - next.previous = down; - - e.v2 = v1; - e.v1 = v2; - e.face = fn; - e.center = center; - e.next = down; - moveBack(ei); - } - addFace(v0, v2, v1); - return 1; - } - - - - /* return new vertex and the center of the new sphere pivoting from edge - if the vertex belongs to another edge, touch points to it. */ - bool pivot(Edgex &edge) { - Point3x v0 = mesh.vert[edge.v0].P(); - Point3x v1 = mesh.vert[edge.v1].P(); - Point3x v2 = mesh.vert[edge.v2].P(); - /* TODO why using the face normals everything goes wrong? should be - exactly the same................................................ - Check if the e.face is correct. - Point3x &normal = mesh.face[edge.face].N(); - */ - - Point3x normal = ((v1 - v0)^(v2 - v0)).Normalize(); - Point3x middle = (v0 + v1)/2; - Point3x start_pivot = edge.center - middle; - Point3x axis = (v1 - v0); - - ScalarType axis_len = axis.SquaredNorm(); - if(axis_len > 4*radius*radius) return false; - 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 targets; - std::vector dists; - getInSphere(middle, r + radius, targets, dists); - - if(targets.size() == 0) { - return false; //this really would be strange but one never knows. - } - - edge.candidate = -1; - ScalarType minangle = 0; - Point3x center; //to be computed for each sample - for(int i = 0; i < targets.size(); i++) { - int id = targets[i]; - - if(id == edge.v0 || id == edge.v1 || id == edge.v2) continue; - - if(mesh.vert[id].IsD()) - continue; - - Point3x p = mesh.vert[id].P(); - + 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); + /* 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(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; - } - //scale alpha by distance: - if(edge.candidate == -1 || alpha < edge.angle) { - edge.candidate = id; - edge.angle = alpha; - edge.newcenter = center; - } - } + /* 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(edge.candidate == -1) { - return false; - } - Point3x n = ((mesh.vert[edge.candidate].P() - v0)^(v1 - v0)).Normalize(); - //found no point suitable. - if(normal * mesh.vert[edge.candidate].N() < 0 || - n * normal < crease || - nb[edge.candidate] >= 2) { - return false; - } - - assert(edge.candidate != edge.v0 && edge.candidate != edge.v1); - return true; - } + if(candidate == NULL) { + return -1; + } + assert((candidate->P() - v0).Norm() > min_edge); + assert((candidate->P() - v1).Norm() > min_edge); - private: - //front management: - //Add a new edge to the back of the queue - typename std::list::iterator newEdge(Edgex e) { - e.active = true; - return front.insert(front.end(), e); - } - //move an Edge among the dead ones - void killEdge(typename std::list::iterator e) { - (*e).active = false; - deads.splice(deads.end(), front, e); - } - - void erase(typename std::list::iterator e) { - if((*e).active) front.erase(e); - else deads.erase(e); - } - //move an Edge to the back of the queue - void moveBack(typename std::list::iterator e) { - front.splice(front.end(), front, e); - } + int id = candidate - &*this->mesh.vert.begin(); + assert(id != edge.v0 && id != edge.v1); + + Point3x newnormal = ((candidate->P() - v0)^(v1 - v0)).Normalize(); + if(normal * newnormal < max_angle || this->nb[id] >= 2) { + return -1; + } + //test if id is in some border (to return touch + for(list::iterator k = this->front.begin(); k != this->front.end(); k++) + if((*k).v0 == id) touch = k; + for(list::iterator k = this->deads.begin(); k != this->deads.end(); k++) + if((*k).v0 == id) touch = 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. + + 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 */ - void moveFront(typename std::list::iterator e) { - front.splice(front.begin(), front, e); - } - - 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; -// i = 0; - if(i < 0) i = 0; - CVertex *start = &*mesh.vert.begin(); - for(; i < mesh.face.size(); i++) { - CFace &f = mesh.face[i]; - for(int k = 0; k < 3; k++) { - if(v1== f.V(k)-start && v0 == f.V((k+1)%3)-start) ++tot; - else if(v0 == f.V(k)-start && v1 == f.V((k+1)%3)-start) { //orientation non constistent - return false; - } - } - if(tot >= 2) { //non manifold - return false; - } - } - return true; - } - - - void cluster(int v) { - /* clean up too close points */ - std::vector targets; - std::vector dists; - getInSphere(mesh.vert[v].P(), mindist*radius, targets, dists); - - for(int i = 0; i < targets.size(); i++) { - int id = targets[i]; - if(id == v) continue; - - CVertex &v = mesh.vert[id]; - if(v.IsD() || v.IsV() || v.IsB()) continue; - v.SetD(); - --mesh.vn; - } - - } - - bool trovamiunnome(typename std::list::iterator e) { - return glue((*e).previous, e) || glue(e, (*e).next); - } - - //glue toghether a and b (where a.next = b - bool glue(typename std::list::iterator a, typename std::list::iterator b) { - if((*a).v0 != (*b).v1) return false; - - typename std::list::iterator previous = (*a).previous; - typename std::list::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].SetV(); - mesh.vert[v].ClearB(); - } - } - - - /* 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*q); - if(vec*axis < 0) angle = -angle; - if(angle < 0) angle += 2*M_PI; - return angle; - } - /* add a new face. compute normals. */ -/* void addFace(int a, int b, int c) { - CFace face; - face.V(0) = (CVertex *)a; - face.V(1) = (CVertex *)b; - face.V(2) = (CVertex *)c; - Point3x &p0 = mesh.vert[a].P(); - Point3x &p1 = mesh.vert[b].P(); - Point3x &p2 = mesh.vert[c].P(); - face.N() = ((p1 - p0)^(p2 - p0)).Normalize(); - - mesh.face.push_back(face); - mesh.fn++; - } */ - - void addFace(int v0, int v1, int v2) { - assert(v0 < mesh.vert.size() && v1 < mesh.vert.size() && v2 < mesh.vert.size()); - CFace 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++; - //update topology? - // FaceIterator f = Allocator::AddFaces(mesh, 1); - } - - void addVertex(CVertex &vertex) { - CVertex *oldstart = &*mesh.vert.begin(); - // VertexIterator v = Allocator::AddVertices(mesh, 1); - mesh.vert.push_back(vertex); - mesh.vn++; - CVertex *newstart = &*mesh.vert.begin(); - if(oldstart != newstart) { - for(int i = 0; i < mesh.face.size(); i++) { - CFace &face = mesh.face[i]; - for(int k = 0; k < 3; k++) - face.V(k) = newstart + (face.V(k) - oldstart); - } - } - nb.push_back(0); - assert(nb.size() == mesh.vert.size()); - } - - - - /* intersects segment [v0, v1] with the sphere of radius radius. */ - bool intersect(int v0, int v1, Point3x ¢er) { - Point3x m = mesh.vert[v1].P() - mesh.vert[v0].P(); - ScalarType t = m*(center - mesh.vert[v0].P()); - if(t < 0) return false; - if(t > m*m) return false; - return true; - } - - - ScalarType distance(int v0, int v1, Point3x ¢er) { - Point3x m = mesh.vert[v1].P() - mesh.vert[v0].P(); - ScalarType t = m*(center - mesh.vert[v0].P())/(m*m); - Point3x p = mesh.vert[v0].P() + m*t; - return (p - center).Norm(); - } - - - /* return all point in a given ball, notice as i want the index - of the vertices not the pointers... this may change in future */ - unsigned int getInSphere(Point3x &p, ScalarType distance, - std::vector &results, - std::vector &dists) { - std::vector ptr; - std::vector points; - int n = trimesh::GetInSphereVertex(mesh, grid, p, distance, ptr, dists, points); - for(int i = 0; i < ptr.size(); i++) - results.push_back(ptr[i] - &(mesh.vert[0])); - return n; - } - - - - /* 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 ¢er) { - Point3x q1 = p1 - p0; - Point3x q2 = p2 - p0; - - 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*q1; - ScalarType a12 = q1*q2; - ScalarType a22 = q2*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 += p0 + up*height; - - return true; - } - - typename std::list::iterator touches(typename std::list::iterator e) { - //TODO what happens when it touches more than one front? - //might still work. - int v = (*e).candidate; - typename std::list::iterator touch = front.end(); - if(mesh.vert[v].IsB()) { - //test nearby Edges: it is faster - typename std::list::iterator p = e; - p = (*e).previous; - if(v == (*p).v0) return p; - e = (*e).next; - if(v == (*e).v0) return e; - - p = (*p).previous; - if(v == (*p).v0) return p; - e = (*e).next; - if(v == (*e).v0) return e; - - //test all. sigh. - - for(typename std::list::iterator k = front.begin(); k != front.end(); k++) { - if(v == (*k).v0) { - touch = k; - break; - } - } - for(typename std::list::iterator k = deads.begin(); k != deads.end(); k++) { - if(v == (*k).v0) { - touch = k; - break; - } - } - assert(touch != front.end()); - } - - return touch; - } - - - public: - - }; - - -}//namespace -}//namespace -/* CODE FOR PIVOTING IN A TOUCH SITUATION not used now. - - //if touch we want to check the ball could really pivot around that point - if(touch != front.end() && touch != (*Edge.next).next && touch != Hinge.previous) { - Point3x &hinge = mesh.vert[min].P(); - Point3x target = (*touch).center - hinge; - float d = (target * start_pivot)/(target.Norm()*start_pivot.Norm()); - if(d < -0.8) { - return false; - } - - - - if(d < 0.5) { //they are far enough so test . - Point3x naxis = (target ^ start_pivot).Normalize(); - Point3x d1 = naxis^start_pivot; - Point3x d2 = target^naxis; - - for(int i = 0; i < targets.size(); i++) { - int id = targets[i]; - if(id == Hinge.v0 || id == Hinge.v1 || id == Hinge.v2 || id == min) continue; - if(mesh.vert[id].IsD()) { - continue; - } - - Point3x intruder = mesh.vert[targets[i]].P() - hinge; - float h = intruder*naxis; - if(fabs(h) > radius) continue; - intruder -= naxis*h; - assert(fabs(intruder *naxis) < 0.01); - float off = radius - intruder.Norm(); //(distance from the center ring of the thorus - if(h*h + off*off > radius*radius) continue; //outside of thorus - if(d1*intruder < 0 || d2*intruder < 0) continue; //ouside of sector - cout << "could not pivot while touching;\n"; - return false; - } - - } - }*/ - + bool FindSphere(Point3x &p0, Point3x &p1, Point3x &p2, Point3x ¢er) { + Point3x q1 = p1 - p0; + Point3x q2 = p2 - p0; + + 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*q1; + ScalarType a12 = q1*q2; + ScalarType a22 = q2*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 += p0 + 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*q); + if(vec*axis < 0) angle = -angle; + if(angle < 0) angle += 2*M_PI; + return angle; + } + + void Mark(VertexType *v) { + std::vector targets; + std::vector points; + std::vector dists; + int n = trimesh::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