CGAL users discussion list

[Jonas Schild <>]

Hi again,

with all your help I could gladly solve the issue. I first tried 
Santosh's code which worked. Though this was quite inconvenient, by 
writing and reading streams, it brought me to looking how the alpha 
shape is put into the ostream in alpha_shape_3.h. There I could find the 
same solution Monique and Mariette kindly proposed. So thank you all!
For anyone interested here is the code:

 int number_of_vertices = 0;
 std::vector<CgalPoint3> vertices_list;
 std::map<CgalVertexHandle,int> vertices_map;

 CgalVerticesIterator vit = mAlphaShape->alpha_shape_vertices_begin();
 for( ; vit != mAlphaShape->alpha_shape_vertices_end(); ++vit)
 {
   vertices_list.push_back((*vit)->point());
   vertices_map[*vit] = number_of_vertices++;
 }

 CgalCellHandle c;
 CgalFacetsIterator fit = mAlphaShape->alpha_shape_facets_begin();
 std::vector<osg::Vec3> triangle(3);
 CgalPoint3 v1, v2, v3;
 int i, i0, i1, i2, j;

 for( ; fit != mAlphaShape->alpha_shape_facets_end(); ++fit)
 {
   c = fit->first;
   i = fit->second;
   // the following ensures that regular facets are output
   // in ccw order
   if ((mAlphaShape->classify(*fit) == CgalAlphaShape3::REGULAR) &&
       (mAlphaShape->classify(c) == CgalAlphaShape3::INTERIOR))
   {
     c = c->neighbor(i);
     i = c->index(fit->first);
   }
   i0 = CgalAlphaShape3::Triangulation_utils_3::vertex_triple_index(i,0);
   i1 = CgalAlphaShape3::Triangulation_utils_3::vertex_triple_index(i,1);
   i2 = CgalAlphaShape3::Triangulation_utils_3::vertex_triple_index(i,2);

   v1 = vertices_list[vertices_map[c->vertex(i0)]];
   v2 = vertices_list[vertices_map[c->vertex(i1)]];
   v3 = vertices_list[vertices_map[c->vertex(i2)]];

   triangle[0] = osg::Vec3(v1[0], v1[1], v1[2]);
   triangle[1] = osg::Vec3(v2[0], v2[1], v2[2]);
   triangle[2] = osg::Vec3(v3[0], v3[1], v3[2]);

   osg::Vec3 a = triangle[1] - triangle[0];
   osg::Vec3 b = triangle[2] - triangle[1];
   osg::Vec3 n = a ^ b;
   n.normalize();

 }

Best regards
Jonas

Santosh Tiwari wrote:
> Attached with this email is the complete code.
>
> I have modified the code and removed unnecessary stuff.
>
> I have not attempted to compile it after the modification. But, you can
> use the code as a blueprint to get the correct normals.
>
> Hope that works.
>
>   
> > Santosh Tiwari wrote:
> >     
> > > This is how I got the correct normals for all the boundary facets.
> > >
> > >
> > >       
> > Hi Santosh,
> > I tried using your code, but failed with correctly setting up the Facet
> > and Vertex data types, so I was not able to create the lists correctly.
> > I then only used your normals to my vertices but they don't seem to be
> > correct. My vertices seem to be in a different order. Could you help me
> > out with your typedefs?
> >
> > Thanks for kindly providing your code!
> > Jonas
> >     
> > > The relevant code.
> > >
> > >
> > > Alpha_shape_3 *as;
> > > list<Facet> facets;
> > > list<Vertex> vertices;
> > > ostringstream facetStream;
> > >
> > > as->get_alpha_shape_facets(back_inserter(facets),
> > > Alpha_shape_3::REGULAR);
> > > as->get_alpha_shape_vertices(back_inserter(vertices),
> > > Alpha_shape_3::REGULAR);
> > >
> > > cout<<"\n Number of boundary facets = "<<facets.size();
> > > cout<<"\n Number of vertices = "<<vertices.size();
> > >
> > > facetStream<<vertices.size()<<endl;
> > > facetStream<<facets.size()<<endl;
> > > facetStream<<*as;
> > >
> > >
> > > //You can now read the facetStream. Here is some code to read it.
> > >
> > > string facetString(facetstream.str());
> > > istringstream in;
> > > in.str(facetString);
> > > in>>numVertices;
> > > in>>numFacets;
> > >
> > > vertices = new double*[numVertices];
> > > for (i=0; i<numVertices; i++) {
> > >         vertices[i] = new double[3];
> > >         in>>x>>y>>z;
> > >         vertices[i][0] = x;
> > >         vertices[i][1] = y;
> > >         vertices[i][2] = z;
> > > }
> > >
> > > facets = new unsigned long int*[numFacets];
> > > for (i=0; i<numFacets; i++) {
> > >         facets[i] = new unsigned long int[3];
> > >         in>>v1>>v2>>v3;
> > >         facets[i][0] = v1;
> > >         facets[i][1] = v2;
> > >         facets[i][2] = v3;
> > > }
> > >
> > > for(i=0; i<numFacets; i++) {
> > >         x1 = vertices[facets[i][0]][0];
> > >         y1 = vertices[facets[i][0]][1];
> > >         z1 = vertices[facets[i][0]][2];
> > >
> > >         x2 = vertices[facets[i][1]][0];
> > >         y2 = vertices[facets[i][1]][1];
> > >         z2 = vertices[facets[i][1]][2];
> > >
> > >         x3 = vertices[facets[i][2]][0];
> > >         y3 = vertices[facets[i][2]][1];
> > >         z3 = vertices[facets[i][2]][2];
> > > }
> > >
> > > //You now have all the facets. You can write them to a STL file.
> > > //You can compute the normals for every facet now.
> > > //Here is some code to compute the normals.
> > >
> > > a1 = x3 - x2;
> > > a2 = y3 - y2;
> > > a3 = z3 - z2;
> > > b1 = x2 - x1;
> > > b2 = y2 - y1;
> > > b3 = z2 - z1;
> > > n1 = a2*b3 - a3*b2;
> > > n2 = a1*b3 - a3*b1;
> > > n3 = a1*b2 - a2*b1;
> > > res = sqrt(n1*n1 + n2*n2 + n3*n3);
> > > n1 /= res;
> > > n2 /= res;
> > > n3 /= res;
> > >
> > > I personally do not prefer accessing the internal data structure of a
> > > library even if it permits that. It is very difficult to debug unless
> > > you
> > > know the data structure really well.
> > >
> > > __
> > > Santosh Tiwari
> > > EIB-326 Clemson University
> > > http://www.clemson.edu/~stiwari/
> > >
> > >
> > >
> > > -----Original Message-----
> > > From: 
> > >
> > > [mailto:]
> > > Sent: Sunday, May 18, 2008 10:54 AM
> > > To: 
> > >
> > > Subject: Re: [cgal-discuss] Normal generation out of 3d alpha shapes
> > >
> > > Andreas Fabri wrote:
> > >
> > >       
> > > > Jonas Schild wrote:
> > > >
> > > >
> > > >         
> > > > > Hi all,
> > > > >
> > > > > I have the following problem:
> > > > >
> > > > > I'm constructing a 3d alpha shape from a set of points using CGAL.
> > > > > After creation of the alpha shape, I'm iterating over all facets to
> > > > > extract the triangles. But:
> > > > > The orientation of the normal (either by using cross product or using
> > > > > CGAL's ortho_vector(point, point, point) of which results are the
> > > > > same) is partially wrong. I've read that the orientation of the
> > > > > triangle is based on orientation of the cell (i.e. pointing into the
> > > > > cell). However, I'm not able to retrieve the normals pointing outside
> > > > > the alpha shape. Part of the problem is, that for certain alpha
> > > > > values, the constructed shape is no solid object. Is it possible to
> > > > > iterate only over alphas where the shape is solid, not necessarily
> > > > > convex?
> > > > >
> > > > > Furthermore, I don't understand how ccw(int) on the triangulation
> > > > > structure is supposed to work regarding finding the correct
> > > > > orientation thus the normals are pointing outside the volume.
> > > > >
> > > > >
> > > > >           
> > > > Hi Jonas,
> > > >
> > > > Here comes a partial answer.
> > > >
> > > > With ccw(int) you refer probably to the base class of all 3D
> > > > traingulations:
> > > >
> > > >
> > > >         
> > > http://www.cgal.org/Manual/3.3/doc_html/cgal_manual/TriangulationDS_3_ref/Cl
> > > ass_Triangulation_utils_3.html#Cross_link_anchor_994
> > >
> > >       
> > > > As the manual says ccw(int) only makes sense when it is 2D.
> > > >
> > > > What you should use is the undocumented funcyion:
> > > >
> > > > static int Triangulation_utils_3<..>::vertex_triple_index(const int i,
> > > > const int j)
> > > >   {
> > > >     // indexes of the  jth vertex  of the facet of a cell
> > > >     // opposite to vertx i
> > > >       CGAL_triangulation_precondition( ( i >= 0 && i < 4 ) &&
> > > >                              ( j >= 0 && j < 3 ) );
> > > >     return tab_vertex_triple_index[i][j];
> > > >   }
> > > >
> > > > The fact that it is undocumented is a bug.
> > > >
> > > >         
> > > really? I would say it is more like a matter of taste...
> > >
> > >         Monique
> > >
> > >       
> > --
> > You are currently subscribed to cgal-discuss.
> > To unsubscribe or access the archives, go to
> > https://lists-sop.inria.fr/wws/info/cgal-discuss
> >
> >     
>
>
>   
> ------------------------------------------------------------------------
>
> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
> #include <CGAL/Delaunay_triangulation_3.h>
> #include <CGAL/Triangulation_hierarchy_3.h>
> #include <CGAL/Alpha_shape_3.h>
> #include <CGAL/IO/io.h>
>
> #include <iostream>
> #include <sstream>
> #include <string>
> #include <fstream>
> #include <list>
>
> using namespace std;
>
> struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
>
> typedef CGAL::Alpha_shape_vertex_base_3<K>                                    
>     Vb;
> typedef CGAL::Triangulation_hierarchy_vertex_base_3<Vb>           Vbh;
> typedef CGAL::Alpha_shape_cell_base_3<K>                                  Fb;
> typedef CGAL::Triangulation_data_structure_3<Vbh,Fb>              Tds;
> typedef CGAL::Delaunay_triangulation_3<K,Tds>                             
> Delaunay;
> typedef CGAL::Triangulation_hierarchy_3<Delaunay>                 
> Delaunay_hierarchy;
> typedef CGAL::Alpha_shape_3<Delaunay_hierarchy>                           
> Alpha_shape_3;
> typedef K::Point_3                                                            
>                           Point;
> typedef Alpha_shape_3::Alpha_iterator                                         
>   Alpha_iterator;
> typedef Alpha_shape_3::NT                                                     
>                   NT;
> typedef Alpha_shape_3::Facet                                                  
>           Facet;
> typedef Alpha_shape_3::Vertex_handle                                          
>                   Vertex;
>
> void writeData(string& facetString);
>
> int main(int argc, const char* const argv[]) {
>
>         if(argc!=2) {
>                 cerr<<"\n Usage: "<<argv[0]<<" point_cloud.dat"<<endl;
>                 exit(1);
>         }
>
>         int mode;
>         mode = static_cast<int>(atoi(argv[1]));
>
>         // Read the points and compute delaunay triangulation
>         Delaunay_hierarchy dt;
>         ifstream is(argv[1]);
>         if(!is) {
>                 cerr<<"\n Could not open file "<<argv[1]<<" for reading, exiting 
> \n";
>                 exit(1);
>         }
>         register unsigned long int n;
>         is>>n;
>         Point p;
>         cout<<"\n "<<n<<" points in the file";
>         for(register unsigned long int i=0; i<n; i++) {
>                 is>>p;
>                 dt.insert(p);
>                 if(i%10000==0) {
>                         cout<<"\n "<<i<<" points read and added to delaunay 
> triangulation";
>                 }
>         }
>         cout<<"\n "<<n<<" points read and added to delaunay triangulation";
>         cout<<"\n Delaunay computed";
>
>         // compute alpha shape
>         Alpha_shape_3 *as = new Alpha_shape_3(dt);
>         cout<<"\n Alpha shape computed";
>
>         // find optimal alpha values
>         Alpha_shape_3::NT alpha_solid = as->find_alpha_solid();
>         Alpha_iterator opt = as->find_optimal_alpha(1);
>         cout<<"\n Smallest alpha value to get a solid through data points is " 
> <<alpha_solid;
>         cout<<"\n Optimal alpha value to get one connected component is 
> "<<*opt;
>         as->set_alpha(*opt);
>         assert(as->number_of_solid_components() == 1);
>         
>         //write the facets
>         list<Facet> facets;
>         list<Vertex> vertices;
>         as->get_alpha_shape_facets(back_inserter(facets), 
> Alpha_shape_3::REGULAR);
>         as->get_alpha_shape_vertices(back_inserter(vertices), 
> Alpha_shape_3::REGULAR);
>
>         cout<<"\n Number of boundary facets = "<<facets.size();
>         cout<<"\n Number of vertices = "<<vertices.size();
>         ostringstream facetStream;
>         facetStream<<vertices.size()<<endl;
>         facetStream<<facets.size()<<endl;
>         facetStream<<*as;
>         writeData(facetStream.str());
>
>         cout<<endl;
>         return 0;
> }
>
> void writeData(string& facetString) {
>         istringstream in;
>         double** vertices;
>         unsigned long int** facets;
>         register unsigned long int i;
>         double x, y, z;
>         unsigned long int v1, v2, v3;
>         double n1, n2, n3, x1, y1, z1, x2, y2, z2, x3, y3, z3;
>         double a1, a2, a3, b1, b2, b3, res;
>         register unsigned long int numVertices;
>         register unsigned long int numFacets;
>
>         in.str(facetString);
>         in>>numVertices;
>         in>>numFacets;
>
>         vertices = new double*[numVertices];
>         for (i=0; i<numVertices; i++) {
>                 vertices[i] = new double[3];
>                 in>>x>>y>>z;
>                 vertices[i][0] = x;
>                 vertices[i][1] = y;
>                 vertices[i][2] = z;
>                 //cout<<vertices[i][0]<<", "<<vertices[i][1]<<", 
> "<<vertices[i][2]<<endl;
>         }
>
>         facets = new unsigned long int*[numFacets];
>         for (i=0; i<numFacets; i++) {
>                 facets[i] = new unsigned long int[3];
>                 in>>v1>>v2>>v3;
>                 facets[i][0] = v1;
>                 facets[i][1] = v2;
>                 facets[i][2] = v3;
>         }
>
>         for(i=0; i<numFacets; i++) {
>                 
>                 x1 = vertices[facets[i][0]][0];
>                 y1 = vertices[facets[i][0]][1];
>                 z1 = vertices[facets[i][0]][2];
>
>                 x2 = vertices[facets[i][1]][0];
>                 y2 = vertices[facets[i][1]][1];
>                 z2 = vertices[facets[i][1]][2];
>
>                 x3 = vertices[facets[i][2]][0];
>                 y3 = vertices[facets[i][2]][1];
>                 z3 = vertices[facets[i][2]][2];
>
>                 a1 = x2 - x1;
>                 a2 = y2 - y1;
>                 a3 = z2 - z1;
>                 b1 = x3 - x2;
>                 b2 = y3 - y2;
>                 b3 = z3 - z2;
>                 n1 = a2*b3 - a3*b2;
>                 n2 = a3*b1 - a1*b3;
>                 n3 = a1*b2 - a2*b1;
>                 
>                 res = sqrt(n1*n1 + n2*n2 + n3*n3);
>                 n1 = n1/res;
>                 n2 = n2/res;
>                 n3 = n3/res;
>                 
>                 //Do whatever you want with facets and normals.
>         }
>
>         return;
> }