Allink  v0.1
ElPolyDrawSkinCGAL.cpp
1 /***********************************************************************
2 ElPoly:This progam provide a graphical visualisation of the data
3 opend by VarData using openGL glut. The most important option are
4 the possibility of changing the backfold of the polymers with 'c',
5 see the subsequent file in the list with '>', see the bond with 'b'.
6 Copyright (C) 2008-2010 by Giovanni Marelli <sabeiro@virgilio.it>
7 
8 
9 This program is free software; you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation; either version 2 of the License, or
12 (at your option) any later version.
13 
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18 
19 You should have received a copy of the GNU General Public License
20 along with this program; If not, see <http://www.gnu.org/licenses/>.
21 ***********************************************************************/
22 #include "ElPoly.h"
23 
24 #ifdef __glut_h__
25 extern Draw *Dr;
26 #ifdef USE_CGAL
27 #ifndef __CGAL_h__
28 #include <fstream>
29 #include <list>
30 // #include <CGAL/basic.h>
31 // #include <CGAL/Cartesian.h>
32 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
33 //#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
34 //#include <CGAL/make_skin_surface_mesh_3.h>
35 #include <CGAL/Skin_surface_3.h>
36 #include <CGAL/mesh_skin_surface_3.h>
37 #include <CGAL/subdivide_skin_surface_mesh_3.h>
38 //#include <CGAL/Skin_surface_3.h>
39 #include <CGAL/Polyhedron_3.h>
40 #include "skin_surface_writer.h"
41 #include <CGAL/IO/Polyhedron_iostream.h>
42 #include <CGAL/subdivide_skin_surface_mesh_3.h>
43 #include <CGAL/Skin_surface_refinement_policy_3.h>
44 //#include <CGAL/mesh_skin_surface_3.h>
45 //#include <CGAL/subdivide_skin_surface_mesh_3.h>
46 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
47 typedef CGAL::Skin_surface_traits_3<K> Traits;
48 typedef CGAL::Skin_surface_3<Traits> Skin_surface_3;
49 typedef Skin_surface_3::FT FT;
50 typedef Skin_surface_3::Weighted_point Weighted_point;
51 //typedef CGAL::Weighted_point<Bare_point,K::RT> Weighted_point;
52 typedef Weighted_point::Point Bare_point;
53 // typedef K::Point_3 Bare_point;
54 typedef CGAL::Polyhedron_3<K,
55  CGAL::Skin_surface_polyhedral_items_3<Skin_surface_3> > Polyhedron;
56 //typedef CGAL::Polyhedron_3<K> Polyhedron;
57 typedef Polyhedron::Traits::Vector_3 Vector;
58 typedef Polyhedron::Vertex_iterator Vertex_iterator;
59 typedef Polyhedron::Halfedge_iterator Halfedge_iterator;
60 typedef Polyhedron::Facet_iterator Facet_iterator;
61 typedef Polyhedron::Halfedge_around_facet_circulator HFC;
62 typedef Polyhedron::Vertex_handle Vertex_handle;
63 // CGAL::Skin_surface_refinement_policy_3<SkinSurface, Polyhedron> policy(skin);
64 
65 
66 
67 void ElPoly::DefineSkin(int NSample){
68  std::list<Weighted_point> l;
69  FT shrinkfactor = 0.5;
70  double *Plot = new double[pNType()*CUBE(NSample)];
71  double *Count = new double[CUBE(NSample)];
72  double Thre = 10.;
73  double Radius = pEdge(0)/(double)NSample;
74  for(int p=0;p<pNPart();p++){
75  int t = pType(p);
76  int vx = (int)(pPos(p,0)/pEdge(0)*NSample);
77  int vy = (int)(pPos(p,1)/pEdge(1)*NSample);
78  int vz = (int)(pPos(p,2)/pEdge(2)*NSample);
79  int vTot = (vz*NSample+vy)*NSample+vx;
80  Plot[vTot*pNType()+t] += 1.;
81  }
82  double *Norm = (double *)calloc(pNType(),sizeof(double));
83  for(int t=0;t<pNType();t++){
84  for(int v=0;v<CUBE(NSample);v++){
85  if(Norm[t] < Plot[v*pNType()+t])
86  Norm[t] = Plot[v*pNType()+t];
87  }
88  Norm[t] = Norm[t] <= 0. ? 1. : Norm[t];
89  }
90  for(int vx=0;vx<NSample;vx++){
91  double x = vx*pEdge(0)/(double)NSample;
92  for(int vy=0;vy<NSample;vy++){
93  double y = vy*pEdge(1)/(double)NSample;
94  for(int vz=0;vz<NSample;vz++){
95  double z = vz*pEdge(2)/(double)NSample;
96  int vTot = (vz*NSample+vy)*NSample+vx;
97  if(Plot[vTot*pNType()] > Thre){
98  l.push_front(Weighted_point(Bare_point(x,y,z),Radius));
99  }
100  }
101  }
102  }
103 
104  Polyhedron Polyhe;
105 
106  Skin_surface_3 skin_surface(l.begin(), l.end(), shrinkfactor);
107  CGAL::mesh_skin_surface_3(skin_surface, Polyhe);
108 
109  // CGAL::subdivide_skin_surface_mesh_3(skin_surface, Polyhe);
110 
111  // std::ofstream out("mesh.off");
112  // out << Polyhe;
113 
114  glDeleteLists(Dr->Particles,1);
115  Dr->Particles = glGenLists(1);
116  glNewList(Dr->Particles,GL_COMPILE);
117 
118  // Polyhedron::Facet_iterator fcUp = Polyhe.facets_begin();
119  // for(;fcUp != Polyhe.facets_end(); ++fcUp){
120  // Polyhedron::Supports_facet_halfedge = fcUp.halfedge();
121  // //Halfedge_around_facet_circulator heUp = fcUp.halfedge();
122  // }
123 
124  // for (Vertex_iterator vit = Polyhe.vertices_begin();vit != Polyhe.vertices_end(); vit++){
125  // // Vector n = policy.normal(vit);
126  // // n = n/sqrt(n*n);
127  // cout << vit->point() << std::endl;
128  // Halfedge_iterator heUp = Polyhe.halfedges_begin();
129  // for(;heUp != Polyhe.halfedges_end(); ++heUp){
130  // //Polyhedron::Halfedge_handle Half = *heUp;
131  // Vertex_handle veUp = heUp->vertex();
132  // K::Point_3 pf1 = vit->point();
133  // }
134  // }
135  CGAL::Inverse_index<Vertex_handle> index(Polyhe.vertices_begin(),
136  Polyhe.vertices_end());
137 
138  for(Facet_iterator fi = Polyhe.facets_begin();fi != Polyhe.facets_end(); ++fi) {
139  HFC hc = fi->facet_begin();
140  HFC hc_end = hc;
141  Polyhedron::Vertex_handle vf1 = (*hc).vertex();
142  hc++;
143  Polyhedron::Vertex_handle vf2 = (*hc).vertex();
144  hc++;
145  Polyhedron::Vertex_handle vf3 = (*hc).vertex();
146  hc++;
147  K::Point_3 pf1 = vf1->point();
148  K::Point_3 pf2 = vf2->point();
149  K::Point_3 pf3 = vf3->point();
150  Vettore v1(pf1.x(),pf1.y(),pf1.z());
151  Vettore v2(pf2.x(),pf2.y(),pf2.z());
152  Vettore v3(pf3.x(),pf3.y(),pf3.z());
153  Vettore vN(3);
154  v1.Mult(InvScaleUn);
155  v2.Mult(InvScaleUn);
156  v3.Mult(InvScaleUn);
157  vN = (v1-v2) ^ (v3-v2);
158  //if(vN.Norm() > 2.*AreaMean) continue;
159  double Sfumatura = .3*Mat->Casuale();
160  glColor4f(0.1,.4+Sfumatura,0.2,1.);
161  //glColor4f(HueUp[p].r,HueUp[p].g,HueUp[p].b,HueUp[p].a);
162  DrTria(&v1,&v2,&v3,&vN);
163  glColor4f(1.,.0,0.,1.);
164  DrTriaContour(&v1,&v2,&v3);
165 
166 
167  // glPushMatrix();//Particle
168  // glBegin(GL_LINES);
169  // do {
170  // Polyhedron::Vertex_handle vh = (*hc).vertex();
171  // K::Point_3 pf1 = vh->point();
172  // glVertex3d(pf1.x(),pf1.y(),pf1.z());
173  // } while (++hc != hc_end);
174  // glEnd();
175  // glPopMatrix();//Particle
176  }
177 
178  glEndList();
179 
180  // Tr tr; // 3D-Delaunay triangulation
181  // C2t3 c2t3 (tr); // 2D-complex in 3D-Delaunay triangulation
182 
183  // Surface_3 surface(sphere_function,Sphere_3(CGAL::ORIGIN, 2.));
184  // Surface_mesh_default_criteria_3<Tr> criteria(30., 0.1,0.1);
185  // // meshing surface
186  // make_surface_mesh(c2t3, surface, criteria, CGAL::Non_manifold_tag());
187 
188  // std::cout << "Final number of points: " << tr.number_of_vertices() << "\n";
189 
190  // DT dt;
191  // for(int c=0;c<Gen->NChain;c++){
192  // if(CHAIN_IF_TYPE(Ch[c].Type,CHAIN_UP) )continue;
193  // Point_3<K> ChPos(pPos(c,CLat1),pPos(c,CLat2),pPos(c,CNorm));
194  // dt.insert(ChPos);
195  // }
196  // Face_iterator fcTr = dt.finite_faces_begin();
197  // glDeleteLists(Dr->Particles,1);
198  // Dr->Particles = glGenLists(1);
199  // glNewList(Dr->Particles,GL_COMPILE);
200  // for(;fcTr != dt.faces_end(); ++fcTr){
201  // Vertex_handle vf1 = fcTr->vertex(0),
202  // vf2 = fcTr->vertex(1),vf3 = fcTr->vertex(2);
203  // Point pf1 = vf1->point();
204  // Point pf2 = vf2->point();
205  // Point pf3 = vf3->point();
206  // Vettore v1(pf1.x() - pf2.x(),pf1.y() - pf2.y(),pf1.z()-pf2.z());
207  // Vettore v2(pf3.x() - pf2.x(),pf3.y() - pf2.y(),pf1.z()-pf2.z());
208  // Vettore vN(3);
209  // vN = v1 ^ v2;
210  // DrTira(v1,v2,v3,vN);
211 
212  // }
213  // glEndList();
214  }
215 #else
216 void ElPoly::DefineSkin(){
217  printf("DefSurface without CGAL not implemented\n");
218 }
219 #endif // USE_CGAL
220 #endif // __CGAL_h__
221 #endif //__glut_h__
double InvScaleUn
Define the shrink factor between the box edges.
Definition: ElPoly.h:457
Draw provides the basic configuration of the openGL libraries used in every derived program...
Definition: Draw.h:15
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
void DefineSkin(int NSample)
Find the covering surface for given points.
void DrTriaContour(Vettore *v00, Vettore *v01, Vettore *v11)
Draw the contour of three vectors.
int pType(int p)
Return the type.
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
double Casuale()
Random uniform number.
GLuint Particles
Refers to the list of the total position of the particles which will be generated in another program...
Definition: Draw.h:200
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
int pNType()
of types of the particle
double pPos(int p, int d)
Return back folded position.
void DrTria(Vettore *v00, Vettore *v01, Vettore *v11, Vettore *vN)
single triangle defined by three vectors
int pNPart()
Number of particle.