Allink  v0.1
VarDataCGAL.cpp
1 /* colored_face.C */
2 /* ---------------- */
3 #include <VarData.h>
4 //#ifndef __CGAL_h__
5 #ifdef USE_CGAL
6 #include <CGAL/basic.h>
7 #include <CGAL/Cartesian.h>
8 #include <CGAL/IO/Color.h>
9 #include <CGAL/Triangulation_euclidean_traits_2.h>
10 #include <CGAL/Triangulation_default_data_structure_2.h>
11 #include <CGAL/Triangulation_2.h>
12 #include <CGAL/Delaunay_triangulation_2.h>
13 #include <CGAL/Voronoi_diagram_2.h>
14 #include <CGAL/Delaunay_triangulation_adaptation_traits_2.h>
15 #include <CGAL/Delaunay_triangulation_adaptation_policies_2.h>
16 #include <vector>
17 #include <iostream>
18 #include <fstream>
19 #include <cassert>
20 
21 using namespace CGAL;
22 using namespace std;
23 
24 /* A facet with a color member variable. */
25 template < class Gt >
26 class My_face_base : public Triangulation_face_base_2<Gt>
27 {
28 public:
29  Color color;
30  My_face_base() :
31  Triangulation_face_base_2<Gt>() {}
32  My_face_base(void* v0, void* v1, void* v2) :
33  Triangulation_face_base_2<Gt>(v0,v1,v2) {}
34  My_face_base(void* v0, void* v1, void* v2, void* n0, void* n1, void* n2) :
35  Triangulation_face_base_2<Gt>(v0,v1,v2,n0,n1,n2) {}
36 };
37 
38 typedef Cartesian<double> Rp;
39 typedef Triangulation_euclidean_traits_2<Rp> Gt;
40 typedef Triangulation_vertex_base_2<Gt> Vb;
41 typedef My_face_base<Gt> Fb;
42 typedef Triangulation_default_data_structure_2<Gt,Vb,Fb > Tds;
43 typedef Triangulation_2<Gt,Tds> Triangulation;
44 //typedef Delaunay_graph DG;
45 typedef Delaunay_triangulation_2<Gt,Tds> DT;
46 typedef Delaunay_triangulation_adaptation_traits_2<DT> AT;
47 typedef Delaunay_triangulation_caching_degeneracy_removal_policy_2 <DT> AP;
48 typedef Voronoi_diagram_2<DT,AT,AP> VD;
49 typedef Point_2<Rp> Point;
50 typedef Triangulation::Face_handle Face_handle;
51 typedef Triangulation::Face_iterator Face_iterator;
52 typedef Triangulation::Vertex Vertex;
53 typedef Triangulation::Vertex_handle Vertex_handle;
54 typedef Triangulation::Vertex_iterator Vertex_iterator;
55 // typedef VD::Face_handle Face_handle;
56 // typedef VD::Face_iterator Face_iterator;
57 // typedef VD::Vertex Vertex;
58 // typedef VD::Vertex_handle Vertex_handle;
59 // typedef VD::Vertex_iterator Vertex_iterator;
60 typedef Triangulation::Finite_vertices_iterator Finite_vertices_iterator;
61 
62 void VarData::AreaDistr(double *Distr,double *RadDistr,int Values) {
63  double BoxRadInv = pInvEdge(3);
64  // VD vd;
65  double AreaMean = 3.*.5*Gen->NChain/(Gen->Edge[CLat1]*Gen->Edge[CLat2]);
66  DT dtUp;
67  DT dtDown;
68  //ofstream File2Write("AreaSnapshot.dat");
69  //Triangulation dt;
70  //Point p;
71  for(int c=0;c<Gen->NChain;c++){
72  if(CHAIN_IF_TYPE(Ch[c].Type,CHAIN_UP) ){
73  Point ChPos(Ch[c].Pos[CLat1],Ch[c].Pos[CLat2]);
74  dtUp.insert(ChPos);
75  }
76  else{
77  Point ChPos(Ch[c].Pos[CLat1],Ch[c].Pos[CLat2]);
78  dtDown.insert(ChPos);
79  }
80  }
81  // Face_iterator fc = t.faces_begin();
82  Face_iterator fcTrUp = dtUp.finite_faces_begin();
83  Face_iterator fcTrDown = dtDown.finite_faces_begin();
84  //VD::Face_iterator fcVor = vd.faces_begin();
85  //VD::Face_iterator fcVor = vd.faces_begin();
86  double *Norma = new double[Values];
87  memset(RadDistr,0.,Values*sizeof(double));
88  vector <double> Areas;
89  for(;fcTrUp != dtUp.faces_end(); ++fcTrUp){
90  Vertex_handle vf1 = fcTrUp->vertex(0),
91  vf2 = fcTrUp->vertex(1),vf3 = fcTrUp->vertex(2);
92  Point pf1 = vf1->point();
93  Point pf2 = vf2->point();
94  Point pf3 = vf3->point();
95  Vettore v1(pf1.x() - pf2.x(),pf1.y() - pf2.y(),0.);
96  Vettore v2(pf3.x() - pf2.x(),pf3.y() - pf2.y(),0.);
97  double Area = v1.VetV(&v2);
98  double Distx = Nano->Pos[CLat1] - (pf1.x() + pf2.x() + pf3.x())/3.;
99  double Disty = Nano->Pos[CLat2] - (pf1.y() + pf2.y() + pf3.y())/3.;
100  double Rad = sqrt( SQR(Distx) + SQR(Disty) );
101  if(Rad < Nano->Rad){continue;}
102  int v = (int)(Rad*BoxRadInv*Values);
103  //assert(v >= 0 && v < Values);
104  if(v < 0 || v >= Values){continue;}
105  RadDistr[v] += Area;
106  Areas.push_back(Area);
107  Norma[v] += 1.;
108  //File2Write << vf1->point() << endl << vf2->point()<< endl << vf3->point() << endl << endl << endl ;
109  }
110  for(;fcTrDown != dtDown.faces_end(); ++fcTrDown){
111  Vertex_handle vf1 = fcTrDown->vertex(0),
112  vf2 = fcTrDown->vertex(1),vf3 = fcTrDown->vertex(2);
113  Point pf1 = vf1->point();
114  Point pf2 = vf2->point();
115  Point pf3 = vf3->point();
116  Vettore v1(pf1.x() - pf2.x(),pf1.y() - pf2.y(),0.);
117  Vettore v2(pf3.x() - pf2.x(),pf3.y() - pf2.y(),0.);
118  double Area = v1.VetV(&v2);
119  double Distx = Nano->Pos[CLat1] - (pf1.x() + pf2.x() + pf3.x())/3.;
120  double Disty = Nano->Pos[CLat2] - (pf1.y() + pf2.y() + pf3.y())/3.;
121  double Rad = sqrt( SQR(Distx) + SQR(Disty) );
122  //if(Rad < Nano->Rad){continue;}
123  int v = (int)(Rad*BoxRadInv*Values);
124  if(v < 0 || v >= Values){continue;}
125  RadDistr[v] += Area;
126  Areas.push_back(Area);
127  Norma[v] += 1.;
128  //File2Write << vf1->point() << endl << vf2->point()<< endl << vf3->point() << endl << endl << endl ;
129  }
130  // for(int v=0;v<Values;v++){
131  // RadDistr[v] /= Norma[v] > 0 ? Norma[v] : 1.;
132  // }
133  vector <double>::iterator iArea;
134  for(iArea=Areas.begin();iArea!=Areas.end();iArea++){
135  int v = (int)(*iArea/AreaMean*Values);
136  if(v < 0 || v >= Values){continue;}
137  Distr[v] += 1.;
138  }
139  delete [] Norma;
140  //File2Write.close();
141  return;
142  //Finite_vertices_iterator vc = dt.finite_vertices_begin();
143  // while(vc != dtUp.finite_vertices_end()){
144  // Face_handle fv = vc->face();
145  // Vertex_handle vf1 = fv->vertex(0),
146  // vf2 = fv->vertex(1),vf3 = fv->vertex(2);
147  // //File2Write << vf1->point() << endl << vf2->point() << endl << vf3->point() << endl << endl << endl;
148  // ++vc;
149  // }
150  // while(fcVor != vd.faces_end()){
151  // // VD::Face fc1 = *fcVor;
152  // //VD::Vertex_handle vf1 = fcVor->vertex(0),//fcVor->vertex(0),
153  // //vf2 = fcVor->vertex(1),vf3 = fcVor->vertex(2);
154  // //File2Write << vf1->point() << endl << vf2->point() << endl << vf3->point() << endl << endl << endl;
155  // ++fcVor;
156  // }
157 
158  // Finite_vertices_iterator it1 = dt.finite_vertices_begin(),it2(it1), it3(it1);
159  // ++it2;
160  // ++it3; ++it3;
161  // // Vertex_handle vf0 = fc.vertex(0);
162  // while( it3 != dt.finite_vertices_end()) {
163  // //File2Write << it1->point() << endl << it2->point() << endl << it3->point() << endl << endl << endl;
164  // ++it1; ++it2; ++it3;
165  // }
166 
167  // return 0;
168 }
169 #else
170 void VarData::AreaDistr(double *Distr,double *RadDistr,int Values){
171  printf("AreaDistr without CGAL not implemented\n");
172 }
173 
174 #endif //USE_CGAL
void AreaDistr(double *Distr, double *RadDistr, int NSample)
Calculate the (temporal/radial) area distribution.
Definition: VarDataCGAL.cpp:62
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
STL namespace.
double VetV(const Vettore *u, const Vettore *v)
Vectorial product between two Vettore returns the area.