songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
cgal/Packages/Mesh_2/Mesh.h

1136 lines
28 KiB
C++
Raw Blame History

#ifndef CGAL_MESH_H
#define CGAL_MESH_H
#include <list>
#include <map>
#include <cmath>
#include <iostream.h>
#include <fstream.h>
#include <qwidget.h>
#include <CGAL/basic.h>
#include <CGAL/Constrained_triangulation_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/predicate_classes_2.h>
#include <CGAL/Triangulation_circulators_2.h>
//CGAL_BEGIN_NAMESPACE
using namespace CGAL;
#define INVALID_EDGE Edge(NULL, -1);
#define UPDATE {viewer->repaint(); qApp->processEvents(); float ff; for(float f=0.0; f<1.2e6; f++) ff += sin(f);}
//#define UPDATE
template <class Gt, class Tds> class Triangulation_finite_vertex_circulator :
public Triangulation_vertex_circulator_2<Gt, Tds>
{
};
template <class Gt, class Tds, class Mtraits = void>
class Mesh:
public Constrained_Delaunay_triangulation_2 <Gt, Tds>
{
public:
QWidget *viewer;
typedef Gt Triangulation_geom_traits;
typedef Tds Triangulation_data_structure;
typedef typename Gt::FT FT;
typedef FT Length;
typedef FT Square_length;
typedef CGAL::Triangulation_2<Gt, Tds> Tr;
typedef CGAL::Constrained_triangulation_2<Gt, Tds> Ct;
typedef CGAL::Constrained_Delaunay_triangulation_2<Gt, Tds> CDt;
// typedef typename CDt::Constraint Constraint;
typedef typename CDt::Vertex Vertex;
typedef typename CDt::Edge Edge;
typedef typename CDt::Face_handle Face_handle;
typedef typename CDt::Vertex_handle Vertex_handle;
typedef typename CDt::Finite_faces_iterator Finite_faces_iterator;
typedef std::pair<Vertex_handle,Vertex_handle> Constrained_edge;
typedef typename Gt::Point_2 Point;
typedef std::pair<Point,Point> Constraint;
typedef std::list<Constraint> List_constraints;
typedef Mtraits Mesh_traits;
enum Cluster_status {REDUCED,NON_REDUCED};
struct Cluster {
Cluster_status status;
Length length;
FT alpha_min;
std::map<Vertex_handle, Length> vertices;
};
void write(ostream &f);
void read(istream &f);
void reset() {
cluster_map.clear();
clear();
}
private:
int current_time;
std::multimap<FT, Face_handle> Bad_faces;
std::list<Constrained_edge> c_edge_queue;
std::map<Vertex_handle, int> insertion_time;
std::map<Vertex_handle, FT> insertion_radius_map;
std::multimap<Vertex_handle, Cluster> cluster_map;
public:
//INSERTION-REMOVAL
void refine_mesh(QWidget *v);
void bounds(FT &xmin, FT &ymin,
FT &xmax, FT &ymax,
FT &xcenter, FT &ycenter);
void bounding_box();
Mesh(const Geom_traits& gt=Geom_traits()):CDt(gt){};
Mesh(List_constraints& lc, const Geom_traits& gt=Geom_traits()):CDt(gt)
{
typename List_constraints::iterator lcit=lc.begin();
for( ; lcit != lc.end(); ++lcit)
{
insert( (*lcit).first, (*lcit).second);
}
CGAL_triangulation_postcondition(is_valid());
}
template <class InputIterator>
Mesh(InputIterator first, InputIterator last, const Geom_traits&
gt=Gt()):CDt(gt)
{
while(first != last){
insert((*first).first, (*first).second);
++first;
}
CGAL_triangulation_postcondition(is_valid());
}
private:
void refine_face(Face_handle f);
void refine_edge(Vertex_handle va, Vertex_handle vb);
void create_clusters();
void create_clusters_of_vertex(Vertex_handle v);
void show_clusters();
Vertex_handle insert_middle(Vertex_handle va, Vertex_handle vb);
void cut_cluster(Vertex_handle va, Vertex_handle vb);
void cut_reduced_cluster(Vertex_handle va, Vertex_handle vb);
void cut_cluster_edge(Vertex_handle va, Vertex_handle vb);
Vertex_handle insert_in_c_edge(Vertex_handle va, Vertex_handle vb, Point p);
void fill_edge_queue();
void fill_facette_map();
void process_edge_queue();
void process_facette_map();
void update_c_edge_queue(Vertex_handle va,
Vertex_handle vb,
Vertex_handle vm);
void update_facette_map(Vertex_handle v);
void update_cluster(Vertex_handle va, Vertex_handle vb, Vertex_handle vm);
inline Edge edge_between(Vertex_handle va, Vertex_handle vb);
//CHECK
bool is_encroached(Vertex_handle va, Vertex_handle vb, Point p);
bool is_encroached(Vertex_handle va, Vertex_handle vb);
bool min_insertion_radius(Vertex_handle v, Cluster &c);
bool is_bad(Face_handle f);
bool is_small_angle(Vertex_handle vleft,
Vertex_handle vmiddle,
Vertex_handle vright);
bool is_cluster_reduced(const Cluster&); //look cluster status
bool is_cluster(Vertex_handle va, Vertex_handle vb);
FT shortest_edge_of_cluster(Vertex_handle v, Cluster &cluster);
void check_cluster_status( Cluster&);
// HELPING functions
Square_length shortest_edge(Face_handle f);
FT angle(Vertex_handle vleft, Vertex_handle vmiddle, Vertex_handle vright);
FT circumradius_to_shortest_edge_ratio(Face_handle f);
FT insertion_radius(Vertex_handle v);
Vertex_handle nearest_incident_vertex(Vertex_handle v);
bool find_cluster(Vertex_handle va, Vertex_handle vb, Cluster &c);
inline Vertex_circulator incr(Vertex_circulator &c) {
c++;
if(Vertex_handle(c) == infinite_vertex()) {
c++;
}
return c;
}
inline Vertex_circulator incr_constraint(Vertex_handle v, Vertex_circulator &c)
{
Vertex_circulator cbegin = c;
incr(c);
while(!edge_between(v, c).first->is_constrained(edge_between(v, c).second)
&& c != cbegin ) {
incr(c);
}
return c;
}
inline Vertex_circulator succ(Vertex_circulator c) { return incr(c); }
inline Vertex_circulator succ_constraint(Vertex_handle v,
Vertex_circulator c)
{
return incr_constraint(v, c);
}
inline bool get_conflicting_edges(const Point &p,
list<Constrained_edge> &cel)
{ //cel = constrained edges list
Face_handle start;
int li;
Locate_type lt;
Face_handle fh = locate(p,lt,li, start);
switch(lt) {
case OUTSIDE_AFFINE_HULL:
case VERTEX:
return false;
case FACE:
case EDGE:
case OUTSIDE_CONVEX_HULL:
propagate_conflicting_edges(p,fh,0,cel);
propagate_conflicting_edges(p,fh,1,cel);
propagate_conflicting_edges(p,fh,2,cel);
return true;
}
CGAL_triangulation_assertion(false);
return false;
}
//private:
inline void propagate_conflicting_edges (const Point &p,
Face_handle fh,
int i,
list<Constrained_edge> &cel)
{
if(is_infinite(make_pair(fh,i))) return;
Face_handle fn = fh->neighbor(i);
Edge e(fh,i);
if(!is_encroached(e.first->vertex(cw(e.second)),
e.first->vertex(ccw(e.second)), p)) {
return;
}
if ( fh->is_constrained(i)) {// || ! test_conflict(p,fn)) {
cel.push_back(Constrained_edge(e.first->vertex(cw(e.second)),
e.first->vertex(ccw(e.second))));//Edge(fn, fn->index(fh));
return;
}
int j = fn->index(fh);
propagate_conflicting_edges(p,fn,ccw(j),cel);
propagate_conflicting_edges(p,fn,cw(j),cel);
return;
}
bool is_infinite(Face_handle fh) {
Vertex_handle va, vb, vc;
va=fh->vertex(0);
vb=fh->vertex(1);
vc=fh->vertex(2);
if(va==infinite_vertex() || vb==infinite_vertex() ||
vc==infinite_vertex()) return true;
//check whether the face still exists...
if(!is_face(va, vb, vc)) return true;
return false;
}
bool is_infinite(Edge e) {
Vertex_handle va, vb;
va=e.first->vertex(cw(e.second));
vb=e.first->vertex(ccw(e.second));
if(va==infinite_vertex() || vb==infinite_vertex()) return true;
//check whether the face still exists...
if(!is_edge(va, vb)) return true;
return false;
}
};
template <class Gt, class Tds, class Mtraits>
bool Mesh<Gt, Tds, Mtraits>::
find_cluster(Vertex_handle va, Vertex_handle vb, Cluster &c)
{
// check if vb is in any cluster of va
pair<multimap<Vertex_handle, Cluster>::iterator,
multimap<Vertex_handle, Cluster>::iterator> range =
cluster_map.equal_range(va);
for(multimap<Vertex_handle, Cluster>::iterator it = range.first;
it != range.second; it++)
{
Cluster &cl = it->second;
if(cl.vertices.find(vb)!=cl.vertices.end()) {
c = it->second;
return true;
}
}
return false;
}
//the function that writes a file
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
write(ostream &f)
{
int nedges = 0;
Edge_iterator eit = edges_begin();
while(eit != edges_end()) {
if((*eit).first->is_constrained((*eit).second))
nedges++;
eit++;
}
f<<nedges<<endl;
eit=edges_begin();
while(eit!=edges_end()) {
if((*eit).first->is_constrained((*eit).second)) {
f<<(*eit).first->vertex(cw((*eit).second))->point().x()<<" ";
f<<(*eit).first->vertex(cw((*eit).second))->point().y()<<" ";
f<<(*eit).first->vertex(ccw((*eit).second))->point().x()<<" ";
f<<(*eit).first->vertex(ccw((*eit).second))->point().y()<<endl;
}
eit++;
}
}
//the function that reads a file
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
read(istream &f)
{
int nedges = 0;
clear();
f>>nedges;
for(int n=0; n<nedges; n++) {
FT x1, y1, x2, y2;
f>>x1>>y1>>x2>>y2;
insert(Point(x1, y1), Point(x2, y2));
}
}
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
fill_edge_queue()
{
Edge_iterator ei = edges_begin();
while(ei != edges_end()) {
Vertex_handle va = (*ei).first->vertex(cw((*ei).second));
Vertex_handle vb = (*ei).first->vertex(ccw((*ei).second));
if((*ei).first->is_constrained((*ei).second) &&
is_encroached(va, vb)&&
va!=infinite_vertex()&&
vb!=infinite_vertex()) {
c_edge_queue.push_back(make_pair(va, vb));
}
ei++;
}
}
//it is necessarry for process_facette_map
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
fill_facette_map()
{
Face_iterator fit = faces_begin();
while(fit != faces_end())
{
if(!is_infinite(fit)) {
if( is_bad(fit))
{
Bad_faces.insert(make_pair(
circumradius_to_shortest_edge_ratio(fit), fit));
}
}
fit++;
}
}
//is used by process_edge_queue
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
process_edge_queue()
{
while(! c_edge_queue.empty() )
{
Constrained_edge ce;
ce=c_edge_queue.front();
c_edge_queue.pop_front();
// if ce is not an edge, choose another. should be more robust
while(!is_edge(ce.first, ce.second) && !c_edge_queue.empty()) {
ce=c_edge_queue.front();
c_edge_queue.pop_front();
}
Face_handle fh;
int i;
if ( is_edge( ce.first, ce.second, fh,i))
{
if(fh->is_constrained(i) &&
is_encroached(ce.first, ce.second))
{
//Vertex_handle va = ce.first;
//Vertex_handle vb = ce.second;
refine_edge(ce.first, ce.second);
//c_edge_queue.pop_front();
}
}
}
}
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
process_facette_map()
{
if(Bad_faces.empty())
{
}
while( !Bad_faces.empty())
{
Face_handle Bf = (*(Bad_faces.begin())).second;
while(is_infinite(Bf)) {
Bad_faces.erase(Bad_faces.begin());
if(Bad_faces.empty()) return;
Bf = (*(Bad_faces.begin())).second;
}
Bad_faces.erase(Bad_faces.begin());
Vertex_handle va, vb, vc;
va = Bf->vertex(0);
vb = Bf->vertex(1);
vc = Bf->vertex(2);
if( is_face(va,vb,vc )&&is_bad(Bf))
{
refine_face(Bf);
}
}
}
//this function split all the segments that are encroached
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
refine_edge(Vertex_handle va, Vertex_handle vb)
{
UPDATE;
Face_handle f;
int i;
is_edge(va, vb, f, i);
if( !is_cluster(va,vb) &&
!is_cluster(vb,va) &&
f->is_constrained(i)) {
Vertex_handle vm = insert_middle(va,vb);
update_c_edge_queue(va, vb, vm);
update_facette_map(vm);
return;
}
else if (is_cluster(va,vb) && is_cluster(vb,va)) {
Vertex_handle vm = insert_middle(va,vb);
update_c_edge_queue(va, vb, vm);
update_facette_map(vm);
if(vm.is_null()) {
exit(-1);
}
update_cluster(va,vb,vm);
update_cluster(vb,va,vm);
return;
}
Vertex_handle vaa = is_cluster(va,vb) ? va : vb;
Vertex_handle vbb = is_cluster(va,vb) ? vb : va;
cut_cluster_edge(vaa,vbb);
}
//split all the bad faces
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
refine_face(Face_handle f)
{
Point pc;
Vertex_handle v;
//check if the face still exists
pc = circumcenter(f);
list<Constrained_edge> conflicts;
list<Constrained_edge>::iterator out_conflicts=c_edge_queue.begin();
if(get_conflicting_edges(pc, c_edge_queue))
{
if(!c_edge_queue.empty()){
while(! c_edge_queue.empty() )
{
Constrained_edge ce;
ce=c_edge_queue.front();
c_edge_queue.pop_front();
while(!is_edge(ce.first, ce.second) && !c_edge_queue.empty()) {
ce=c_edge_queue.front();
c_edge_queue.pop_front();
}
Face_handle fh;
int i;
if ( is_edge( ce.first, ce.second, fh,i))
{
if( fh->is_constrained(i) &&
is_encroached(ce.first, ce.second))
{
Vertex_handle va = ce.first;
Vertex_handle vb = ce.second;
Square_length rg = shortest_edge(f);
Cluster cluster;
Vertex_handle vaa = is_cluster(va,vb) ? va : vb;
Vertex_handle vbb = is_cluster(va,vb) ? vb : va;
find_cluster(vaa, vbb, cluster);
Length rmin = shortest_edge_of_cluster(vaa, cluster);
if(rmin >rg) {
refine_edge(ce.first, ce.second);
} else {
exit(-1);
}
//c_edge_queue.pop_front();
}
}
}
}
else
{
Face_handle start;
int li;
Locate_type lt;
Face_handle fh = locate(pc,lt,li, start);
if(lt!=OUTSIDE_CONVEX_HULL) {
v = insert(pc);
update_facette_map(v);
}
}
}
UPDATE;
current_time++;
insertion_time[v] = current_time;
//flip_around(v);
}
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
create_clusters()
{
Vertex_iterator vit = vertices_begin();
while(vit != vertices_end()) {
create_clusters_of_vertex(vit);
vit++;
}
}
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
create_clusters_of_vertex(Vertex_handle v)
{
// prerequisite: at least too vertices must exist in the triangulation
Vertex_circulator vcirc = incident_vertices(v);
incr_constraint(v, vcirc);
Vertex_circulator vbegin = vcirc;
incr_constraint(v, vcirc);
do {
if(!is_small_angle(vcirc, v, succ_constraint(v, vcirc))) {
break;
}
incr_constraint(v, vcirc);
} while(vcirc!=vbegin);
// now vcirc points to:
// * either the point JUST BEFORE the first large angle
// * or vbegin if such an angle has not been encountered
// begin with the second case: a single cluster exists
if(vcirc==vbegin) {
FT min_angle = 1000.0; // a BIG value for an angle
Cluster cl;
do{
if(is_small_angle(vcirc, v, succ_constraint(v, vcirc))) {
cl.vertices.insert(make_pair(vcirc, squared_distance(v->point(),
vcirc->point())));
min_angle = min(min_angle, angle(vcirc, v, succ_constraint(v, vcirc)));
}
incr_constraint(v, vcirc);
}while(vcirc!=vbegin);
cl.alpha_min = min_angle;
cluster_map.insert(make_pair(v, cl));
} else {
// the first case
vbegin = vcirc;
do{
Cluster cl;
cl.alpha_min = 1000;
do {
cl.vertices.insert(make_pair(vcirc, squared_distance(v->point(),
vcirc->point())));
cl.alpha_min = min(cl.alpha_min, angle(vcirc, v,
succ_constraint(v, vcirc)));
incr_constraint(v, vcirc);
} while(is_small_angle(vcirc, v, succ_constraint(v, vcirc)));
cluster_map.insert(make_pair(v, cl));
} while(vcirc != vbegin);
}
}
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
show_clusters()
{
multimap<Vertex_handle, Cluster>::iterator cmit = cluster_map.begin();
while(cmit != cluster_map.end()) {
Vertex_handle v = (*cmit).first;
map<Vertex_handle, Length>::iterator cit = (*cmit).second.vertices.begin();
while(cit != (*cmit).second.vertices.end()) {
Vertex_handle v1 = (*cit).first;
cit++;
}
cmit++;
}
}
//refine the cluster
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
cut_cluster(Vertex_handle va, Vertex_handle vb)
{
Cluster c;
find_cluster(va, vb, c);
check_cluster_status(c);
if( !is_cluster_reduced(c))
{
cut_cluster_edge(va, vb);
}
else
{
cut_reduced_cluster(va, vb);
}
}
//refine the reduced_cluster
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
cut_reduced_cluster(Vertex_handle va, Vertex_handle vb)
{
Cluster c;
find_cluster(va, vb, c);
Vertex_handle vm;
if(is_cluster_reduced(c))
{
// vm = insert_in_c_edge(va, vb, va->point() + (vb->point() - va->point())/2.0);
}
update_c_edge_queue(va, vb, vm);
update_facette_map(vm);
}
//refine the cluster edges
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
cut_cluster_edge(Vertex_handle va, Vertex_handle vb)
{
//seulement les constraints
Cluster c;
find_cluster(va, vb, c);
Square_length l2 = 2.0; //shortest_edge_of_cluster(va, c)/1.001;
Square_length L2 = squared_distance(va->point(), vb->point());
if(L2 > l2)
{
FT rapport = pow(2.0, rint(0.5*log(L2/l2)/log(2.0)))*::sqrt(l2/(L2*4));
Point pc=va->point()+(vb->point()-va->point())*rapport;
Vertex_handle vc = insert_in_c_edge(va,vb,pc);
flip_around(vc);
update_c_edge_queue(va, vb, vc);
update_facette_map(vc);
update_cluster(va, vb, vc);
}
}
template <class Gt, class Tds, class Mtraits>
Mesh<Gt, Tds, Mtraits>::Vertex_handle Mesh<Gt, Tds, Mtraits>::
insert_middle(Vertex_handle va, Vertex_handle vb)
{
Point midpoint;
midpoint = va->point() + (vb->point() - va->point())/2.0 ;
return insert_in_c_edge(va,vb, midpoint);
}
//insert in constraint edge the middle
template <class Gt, class Tds, class Mtraits>
Mesh<Gt, Tds, Mtraits>::Vertex_handle Mesh<Gt, Tds, Mtraits>::
insert_in_c_edge(Vertex_handle va, Vertex_handle vb, Point p)
{
Face_handle f;
int i;
is_edge(va, vb, f, i);
Vertex_handle v = Ct::special_insert_in_edge(p, f, i);
flip_around(v) ;
update_c_edge_queue(va, vb, v);
return v;
}
//update the encroached segments list
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
update_c_edge_queue(Vertex_handle va, Vertex_handle vb, Vertex_handle vm)
{
//c_edge_queue.pop_back();
Face_circulator fc = incident_faces(vm);
Face_circulator fcbegin = fc;
if(fc.is_empty()) return;
do {
int i;
for(i=0; i<3; i++) {
Edge e=Edge(fc, i);
Vertex_handle v1=e.first->vertex(cw(i));
Vertex_handle v2=e.first->vertex(ccw(i));
if(is_encroached(v1, v2)&&
e.first->is_constrained(i)&&
v1!=infinite_vertex()&&v2!=infinite_vertex())
{
c_edge_queue.push_back(Constrained_edge(v1, v2));
}
}
fc++;
} while(fc != fcbegin);
if(is_encroached(va, vm)) {
c_edge_queue.push_back(Constrained_edge(va, vm));
}
if(is_encroached(vb, vm)) {
c_edge_queue.push_back(Constrained_edge(vb, vm));
}
}
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
update_facette_map(Vertex_handle v)
{
Face_circulator fc = v->incident_faces();
Face_circulator fcbegin = fc;
do {
if(!is_infinite(fc)) {
if(is_bad(fc)) {
Bad_faces.insert(make_pair(
circumradius_to_shortest_edge_ratio(fc), fc));
}
}
fc++;
} while(fc!=fcbegin);
}
//ok
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
update_cluster(Vertex_handle va, Vertex_handle vb, Vertex_handle vm)
{
multimap<Vertex_handle, Cluster>::iterator it_va_cluster =
cluster_map.find(va);
int n_va_cluster = cluster_map.count(va);
for(int i=0; i<n_va_cluster; i++)
{
map<Vertex_handle, Length> &vertices = (*it_va_cluster).second.vertices;
map<Vertex_handle, Length>::iterator it_vertices = vertices.find(vb);
if(it_vertices != vertices.end())
{
vertices.erase(vb);
if(va.is_null()) {}
if(vm.is_null()) {}
vertices[vm] = (squared_distance(va->point(), vm->point()));
break;
}
it_va_cluster++;
}
multimap<Vertex_handle, Cluster>::iterator it_vb_cluster =
cluster_map.find(vb);
int n_vb_cluster = cluster_map.count(vb);
for(int i=0; i<n_vb_cluster; i++)
{
map<Vertex_handle, Length> &vertices = (*it_vb_cluster).second.vertices;
map<Vertex_handle, Length>::iterator it_vertices = vertices.find(va);
if(it_vertices != vertices.end())
{
vertices.erase(va);
vertices[vm] = (squared_distance(vb->point(), vm->point()));
break;
}
it_vb_cluster++;
}
}
template <class Gt, class Tds, class Mtraits>
inline Mesh<Gt, Tds, Mtraits>::Edge
Mesh<Gt, Tds, Mtraits>::edge_between(Vertex_handle va, Vertex_handle vb) {
Edge_circulator ec = va->incident_edges();
Edge_circulator ecbegin = ec;
do {
Edge e = (*ec);
Face_handle f = e.first;
int iedge = e.second;
Vertex_handle v1 = f->vertex(cw(iedge));
Vertex_handle v2 = f->vertex(ccw(iedge));
if((v1 == va && v2 == vb) || (v2 == va && v1 == vb)) {
return e;
}
ec++;
} while(ec != ecbegin);
cerr<<"invalid edge"<<endl;
return INVALID_EDGE;
}
//CHECK
template <class Gt, class Tds, class Mtraits>
bool Mesh<Gt, Tds, Mtraits>::
is_encroached(Vertex_handle va, Vertex_handle vb, Point p)
{
Point pm=midpoint(va->point(), vb->point());
if(2.0*::sqrt(squared_distance(pm, p)) < ::sqrt(squared_distance(va->point(), vb->point())))
{
return true;
}
return false;
}
template <class Gt, class Tds, class Mtraits>
bool Mesh<Gt, Tds, Mtraits>::
is_encroached(Vertex_handle va, Vertex_handle vb)
{
Vertex_iterator vi=vertices_begin();
while(vi!=vertices_end())
{
if(is_encroached(va, vb, vi->point())&&
va!=Vertex_handle(vi)&&
vb!=Vertex_handle(vi))
{
return true;
}
vi++;
}
return false;
}
template <class Gt, class Tds, class Mtraits>
bool Mesh<Gt, Tds, Mtraits>::
min_insertion_radius(Vertex_handle v, Cluster &c)
{
}
//the measure of faces quality
template <class Gt, class Tds, class Mtraits>
bool Mesh<Gt, Tds, Mtraits>::
is_bad(Face_handle f)
{
FT quality = circumradius_to_shortest_edge_ratio(f);
if((quality >1) || (quality < 0.5))
//set_a:(>1.0 || <0.4)
//set_b:(>1.0 || <0.5)
//set_c:(>1.1 || <0.5)
//set_d:(>1.1 || <0.4)
//set_e:(>0.9 || <0.5)
//set_f:(>0.9 || <0.4)//~NOK
//set_g:(>0.8 || <0.4)//NOK
//set_h:(>0.8 || <0.5)
{
return TRUE;
}
else
{
return FALSE;
}
// }
}
//
template <class Gt, class Tds, class Mtraits>
bool Mesh<Gt, Tds, Mtraits>::
is_small_angle(Vertex_handle vleft,
Vertex_handle vmiddle,
Vertex_handle vright)
{
FT cos_alpha;
cos_alpha = angle(vleft, vmiddle, vright);
if(cos_alpha > 1/2)
{
return TRUE; //the same cluster
}
else
{
return FALSE; //another cluster
}
}
//
template <class Gt, class Tds, class Mtraits>
bool Mesh<Gt, Tds, Mtraits>::
is_cluster_reduced(const Cluster& c)
{
return c.status == REDUCED;
}
template <class Gt, class Tds, class Mtraits>
bool Mesh<Gt, Tds, Mtraits>::
is_cluster(Vertex_handle va, Vertex_handle vb)
{
// check if vb is in any cluster of va
pair<multimap<Vertex_handle, Cluster>::iterator,
multimap<Vertex_handle, Cluster>::iterator> range =
cluster_map.equal_range(va);
for(multimap<Vertex_handle, Cluster>::iterator it = range.first;
it != range.second; it++)
{
Cluster &cl = it->second;
if(cl.vertices.find(vb)!=cl.vertices.end()) {
return true;
}
}
return false;
}
template <class Gt, class Tds, class Mtraits>
typename Gt::FT Mesh<Gt, Tds, Mtraits>::
shortest_edge_of_cluster(Vertex_handle v, Cluster &cluster)
{
map<Vertex_handle, Length>::iterator vit = cluster.vertices.begin();
FT min_edge = squared_distance(((*vit).first)->point(), (*v).point());
vit++;
while(vit != cluster.vertices.end())
{
FT temp = squared_distance(((*vit).first)->point(), (*v).point());
if(temp < min_edge)
{
min_edge = temp;
}
vit++;
}
return min_edge;
}
// look if all edges have the same lenght
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
check_cluster_status( Cluster& cluster)
{
Length initl;
map<Vertex_handle, Length>::iterator vi = cluster.vertices.begin();
initl = (*(cluster.vertices.find((*vi).first))).second;
vi++;
for(;vi != cluster.vertices.end(); vi++)
{
if(( (*(cluster.vertices.find((*vi).first))).second) != initl)
{
cluster.status = NON_REDUCED;
return;
}
vi++;
}
cluster.status = REDUCED;
}
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
bounds(FT &xmin, FT &ymin,
FT &xmax, FT &ymax,
FT &xcenter, FT &ycenter)
{
Vertex_iterator vi=vertices_begin();
xmin=xmax=xcenter=vi->point().x();
ymin=ymax=ycenter=vi->point().y();
vi++;
while(vi != vertices_end())
{
xcenter+=vi->point().x();
ycenter+=vi->point().y();
if(vi->point().x() < xmin) xmin=vi->point().x();
if(vi->point().x() > xmax) xmax=vi->point().x();
if(vi->point().y() < ymin) ymin=vi->point().y();
if(vi->point().y() > ymax) ymax=vi->point().y();
vi++;
}
xcenter /= number_of_vertices();
ycenter /= number_of_vertices();
}
//bounding box
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
bounding_box()
{
FT span;
FT xmin, xmax, ymin, ymax, xcenter, ycenter;
bounds(xmin, ymin, xmax, ymax, xcenter, ycenter);
span = max((xmax-xmin), (ymax-ymin));
Point bb1(xcenter - 1.5*span, ycenter - 1.5*span);
Point bb2(xcenter + 1.5*span, ycenter - 1.5*span);
Point bb3(xcenter + 1.5*span, ycenter + 1.5*span);
Point bb4(xcenter - 1.5*span, ycenter + 1.5*span);
insert(bb1);
insert(bb2);
insert(bb3);
insert(bb4);
insert(bb1, bb2);
insert(bb2, bb3);
insert(bb3, bb4);
insert(bb4, bb1);
}
//the angle that are between 2 edges from the triangulation
template <class Gt, class Tds, class Mtraits>
typename Gt::FT Mesh<Gt, Tds, Mtraits>::
angle(Vertex_handle vleft, Vertex_handle vmiddle, Vertex_handle vright)
{
Point
pa = vleft->point(),
pb = vmiddle->point(),
pc = vright->point();
FT a, b, c;
a = ::sqrt(squared_distance(pb, pc));
b = ::sqrt(squared_distance(pc, pa));
c = ::sqrt(squared_distance(pa, pb));
FT cos_alpha = ((a*a+c*c-b*b)/(2*a*c));
return cos_alpha;
}
//the shortest edge that are in a triangle
template <class Gt, class Tds, class Mtraits>
typename Gt::FT Mesh<Gt, Tds, Mtraits>::
shortest_edge(Face_handle f)
{
Point
pa = (f->vertex(0))->point(),
pb = (f->vertex(1))->point(),
pc = (f->vertex(2))->point();
FT a, b, c;
a =(squared_distance(pb, pc));
b =(squared_distance(pc, pa));
c =(squared_distance(pa, pb));
return (min(a, min(b, c))); // regarder la documentation du CGAL sur
// MIN (Developer Manual)
}
//the triangle quality is represented by the
//circumradius_to_shortest_edge_ratio value
template <class Gt, class Tds, class Mtraits>
typename Gt::FT Mesh<Gt, Tds, Mtraits>::
circumradius_to_shortest_edge_ratio(Face_handle f)
{
Point p;
p = circumcenter(f);
Point A = (f->vertex(0))->point();
FT radius =::sqrt(squared_distance(p, A));
FT sh_edge =::sqrt(shortest_edge(f));
return (radius/sh_edge);
}
//insertion radius: the definition
template <class Gt, class Tds, class Mtraits>
typename Gt::FT Mesh<Gt, Tds, Mtraits>::
insertion_radius(Vertex_handle v)
{
Vertex_handle v1 = nearest_incident_vertex(v);
return (squared_distance(v->point(), v1->point()));
}
//this function must compute the vertex that are so close to our vertex.
template <class Gt, class Tds, class Mtraits>
Mesh<Gt, Tds, Mtraits>::Vertex_handle Mesh<Gt, Tds, Mtraits>::
nearest_incident_vertex(Vertex_handle v)
{
Vertex_handle vbegin, vcurrent, vnearest;
Vertex_circulator circ = v->incident_vertices();
vbegin = circ;
vnearest = vbegin;
circ++;
FT dist = (squared_distance(v->point(), vbegin->point()));
while((vcurrent = circ) != vbegin) {
FT d1 = (squared_distance(v->point(), vcurrent->point()));
if(d1 < dist) {
vnearest = vcurrent;
}
circ++;
}
return vnearest;
}
//the mesh refine function
template <class Gt, class Tds, class Mtraits>
void Mesh<Gt, Tds, Mtraits>::
refine_mesh( QWidget *w)
{
viewer = w;
// bounding_box();
create_clusters();
show_clusters();
fill_edge_queue();
fill_facette_map();
process_edge_queue();
process_facette_map();
}
//CGAL_END_NAMESPACE
#endif
Reference in New Issue View Git Blame Copy Permalink