with all your help I could gladly solve the issue. I first tried

Santosh's code which worked. Though this was quite inconvenient, by

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!

> 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:

[hidden email]
>>> [mailto:

[hidden email]]

>>> Sent: Sunday, May 18, 2008 10:54 AM

>>> To:

[hidden email]
>>> 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;

> }

You are currently subscribed to cgal-discuss.