CGAL users discussion list

["Santosh Tiwari" <>]

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
>


-- 
Santosh Tiwari
Fluor Daniel EIB 326
Clemson University, Clemson, SC 29634
http://www.clemson.edu/~stiwari/#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;
}