Polyhedron Vertex ID

[petrasvestartas]

How can I get Polyhedron vertex id? I would like to iterate over faces and get vertex id of those face vertices. In code below I am iterating over faces, but I do not know how to get index of a vertex.

 ```
   for ( auto it = output_mesh.facets_begin (); it != output_mesh.facets_end (); ++it ) {
            auto j = it->facet_begin ();
            CGAL_assertion (CGAL::circulator_size (j) >= 3);

            do {
                //j->vertex()->point(). //How can I get an index ?
                //std::cout << ' ' << std::distance (P.vertices_begin (), j->vertex ());
            } while ( ++j != it->facet_begin () );

        }

[raphaelsulzer]

See here: http://cgal-discuss.949826.n4.nabble.com/Convert-Polyhedron-3-Vertex-handles-to-integers-tp4664108p4664110.html

So for the first solution use CGAL::Polyhedron_items_with_id_3. Loop over your mesh vertices and set vertex->id() = i++, which you can then access with j->vertex()->id() in your while loop.

[petrasvestartas]

I added correct header and Polyhedron3 is indexed. I am able to add vertex indices and in other loop print them to console. However, I still cannot fet indices for mrom mesh faces.

The face circulator I am doing is simply wrong. I am missing something basic. Could you advice me why I do not get id when circulating over faces vertices?

I did this:

           int counter = 0;
           for ( Polyhedron3::Vertex_iterator it = output_mesh.vertices_begin (); it < output_mesh.vertices_end (); ++it ) {
               it->id () = counter;
               counter++;
           }

Then I am iterating over face vertex indices:

            int b = 0;
            for ( auto it = output_mesh.facets_begin (); it != output_mesh.facets_end (); ++it ) {
                auto j = it->facet_begin ();

                do {
                    n_o[b++] = j->vertex()->id ();
                } while ( ++j != it->facet_begin () );

            }

But this outputs strange information. What I am doing wrong?

18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615
18446744073709551615

This is the full code:

    Polyhedron3 output_mesh;
        double average_spacing = CGAL::compute_average_spacing<CGAL::Sequential_tag>     (points, 6, CGAL::parameters::point_map (CGAL::First_of_pair_property_map<PointVectorPair> ()));
        if ( CGAL::poisson_surface_reconstruction_delaunay  (points.begin (), points.end (), CGAL::First_of_pair_property_map<PointVectorPair> (),   CGAL::Second_of_pair_property_map<PointVectorPair> (),  output_mesh, average_spacing) ) {

            ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
            //Run Open3D Method
            ////////////////////////////////////////////////////////////////////////////////////////////////////////////////

            p_c_o = output_mesh.size_of_vertices () * 3;
            p_o = new double[p_c_o];

            n_c_o = output_mesh.size_of_facets () * 3;
            n_o = new double[n_c_o];

            c_c_o = output_mesh.size_of_vertices () * 3;
            c_o = new double[c_c_o];

            std::ofstream myfile;
           myfile.open ("C:\\libs\\Cockroach\\CockroachPInvoke\\Cockroach_CSHARP_DLL\\bin\\x64\\PoissonStep1.txt");

           std::size_t counter = 0;
           for ( Polyhedron3::Vertex_iterator it = output_mesh.vertices_begin (); it < output_mesh.vertices_end (); ++it ) {
               it->id () = counter;
               counter++;
           }

           std::size_t counter2 = 0;
           for ( Polyhedron3::Facet_iterator it = output_mesh.facets_begin (); it < output_mesh.facets_end (); ++it ) {
               it->id() = counter2;
               counter2++;
           }

            int i = 0;
            int a = 0;
            int c = 0;
            for ( Polyhedron3::Vertex_iterator it = output_mesh.vertices_begin (); it < output_mesh.vertices_end (); ++it ){

               // it->id () = i++;
                p_o[a++] = it->point ().x ();
                p_o[a++] = it->point ().y ();
                p_o[a++] = it->point ().z ();
                c_o[a++] = 0;
                c_o[a++] = 0;
                c_o[a++] = 0;
                myfile << it->point ().x () << "\n";
                myfile << it->point ().y () << "\n";
                myfile << it->point ().z () << "\n";
            }
            myfile.close ();

            //i = 0;
            //for ( Polyhedron3::Facet_iterator it = output_mesh.facets_begin (); it != output_mesh.facets_end (); ++it )
            //    it->id () = i++;

            //
            std::ofstream myfile2;

            int b = 0;
            myfile2.open ("C:\\libs\\Cockroach\\CockroachPInvoke\\Cockroach_CSHARP_DLL\\bin\\x64\\PoissonStep2.txt");

            for ( Polyhedron3::Face_iterator it = output_mesh.facets_begin (); it != output_mesh.facets_end (); ++it ) {
                Polyhedron3::Halfedge_around_facet_circulator circ =         it->facet_begin ();
                //cout << "Vertex indices of facet: " << it->id ();

                do {
                    n_o[b++] = circ->vertex ()->id ();
                    myfile2 << circ->vertex ()->id () << "\n";
                } while ( ++circ != it->facet_begin () );

            }

            myfile2.close ();

        } 

[sloriot]

Try

for (Polyhedron3::Facet_handle fh : output_mesh.facet_handles())
{
  Polyhedron3::Halfedge_handle start = fh->halfedge(), h = start;
  do{
    std::cout << h->vertex()->id() << "\n";
    h = h->next();
  }while(h!=start);
}

[sloriot]

The issue in your code was it < output_mesh.vertices_end () you should use !=.

[petrasvestartas]

Thank you @sloriot . It worked after changing < to != in the iterator. I added the source code for the PInvoke I use incase somebody else have a similar issue.

However I could not manage to write output_mesh.facet_handles(). facet_handles is missing. Do I need to reference a specific header file?

SourceCode.txt

[afabri]

Note also that it might make sense to use Surface_mesh instead of Polyhedron and the BGL API. While the former gives you the index for free, as a vertex is identified by an index, and the latter gives you more elegant loops.

for(boost::graph_traits<Poyhedron3>::face_descriptor fd : faces(output_mesh)){
  for(boost::graph_traits<Poyhedron3>::vertex_descriptor vd :  vertices_around_face(halfedge(fd, output_mesh), output_mesh){
    int id  = vd; 
  }
}

[afabri]

The function facet_handles() got introduced in 5.2, so you must either upgrade, or use the _begin()/_end() or switch to the BGL API.

[petrasvestartas]

Thank you, @afabri I did not know that it is possible to use Surface mesh as well. I will try to use it. And it seems I am using CGAL 5.1.1

Is there any possibility to prevent the code from crashing? If I input less than 150 points or if points are lying on a plane like the screenshot attach, it crashes.

I tried to use try/catch statement, but it does not prevent from crash. Is there any way to prevent this? (For me it is fine to get error, but I would like to avoid restarting application)

try { 
        double average_spacing = CGAL::compute_average_spacing<CGAL::Sequential_tag>     (points, 6, CGAL::parameters::point_map (CGAL::First_of_pair_property_map<PointVectorPair> ()));
        if ( CGAL::poisson_surface_reconstruction_delaunay  (points.begin (), points.end (), CGAL::First_of_pair_property_map<PointVectorPair> (),   CGAL::Second_of_pair_property_map<PointVectorPair> (),  output_mesh, average_spacing) ) {
...
   }

    } catch ( const std::exception& e ) {

        std::ofstream myfile2;
        myfile2.open ("C:\\libs\\Cockroach\\CockroachPInvoke\\Cockroach_CSHARP_DLL\\bin\\x64\\PoissonDebug.txt");
        myfile2<<e.what();
        myfile2.close();
    }

[afabri]

We should work on this.
And let me make a remark for full code examples. You better put them on gist.github.com so that they are really full, that is with #include and typedef. As a bonus this gets versioned.

[petrasvestartas]

These are two examples files inlcuding header and cpp:

https://gist.github.com/petrasvestartas/7d84b73170fafa38c51b23ad699af5d4 https://gist.github.com/petrasvestartas/10c521c7b2a82ed65536cfe1ebfec200

I know it is PInvoke but at least you can see code better.

I assume, the issue of crashing when input is not good, cannot be solved at a current stage?

[afabri]

You have to check if all points are coplanar. In this case Poisson does not really make sense.

[petrasvestartas]

It is not only about planarity. For instance if input pointset is downsampled to really small number of points i.e. 100 it crashes as well.

I am wondering if there is any way to track if input pointsets are totally irrational like this to avoid running Poisson:

[afabri]

Can you just send me a small point set with normals for which Poisson crashes so that I do not have to invent it myself.

[petrasvestartas]

These are two pointclouds that crashes. If you need other than .ply format let me know.

Clouds.zip

[petrasvestartas]

With normals CloudsWithNormals.zip

[petrasvestartas]

Dear @afabri maybe there was a chance to check the pointclouds?

[afabri]

Sorry @petrasvestartas for coming back to you so late. Do you have a .cpp with your call that leads to a crash, so that I see the parameters etc.

[petrasvestartas]

Hi,

The Poisson Surface Reconstruction starts at line 786: https://gist.github.com/petrasvestartas/c37b9d8089c2c0532cd4a5ac358505df

It is a PInvoke method. I believe it would be straight forward to understand.

[afabri]

I just created the self contained poisson.cpp and it works for Cloud1.ply

[petrasvestartas]

Thank you @afabri for the example file. I am really grateful that CGAL is such an active and responsive library.

I am trying to compare the example I took before and yours. To make your method work do I need to install the recent CGAL version? It seems on older version I cannot call average spacing with two variables and for poisson reconstruction point_map / normal_map does not exist. Would it work with older version as well?

[afabri]

Ideally you cloned what I've put on gist.github.com, so that we have the same code base and can eventually reproduce. It would also be good if for that you used the latest public release (or even the master branch). So it is not so much an example but a minimal self contained test case.

[petrasvestartas]

Hi @afabri ,

I installed CGAL 5.2, seems that the methods / constructor overloads you are using,, are also working on my pc. I linked header files exactly as you did.

Since I am using PInvoke, I would to ask a few question due inputs and outputs. They are simple, but I strictly would like to follow your example using the same data-structures.

1) In your example you are reading a file, in my case I need to add points and vectors to the point set. Previously I was pushing PointVectorPair , can I somehow add points and vectors to the points object of your example? Because there is no push_back method in CGAL::Point_set_3

    CGAL::Point_set_3<Kernel2::Point_3, Kernel2::Vector_3> points;

    for ( int i = 0; i < p_c; i++ ) {

        //points.push_back (PointVectorPair (
        //    Kernel2::Point_3 (p[3 * i + 0], p[3 * i + 1], p[3 * i + 2]),
        //    Kernel2::Vector_3 (n[3 * i + 0], n[3 * i + 1], n[3 * i + 2])
        //));
    }

2) This question about reading vertices / faces from Polyhedron2. The commented line before worked now there not index. How can I iterate through vertices/faces and get each of their id?

Polyhedron2 output_mesh;

for ( auto it = output_mesh.vertices_begin (); it != output_mesh.vertices_end (); ++it ) {
  // it->id () = counter++;
}

std::size_t counter2 = 0;
for ( auto it = output_mesh.facets_begin (); it != output_mesh.facets_end (); ++it ) {
   //it->id () = counter2++;
 }

int b = 0;
for ( auto it = output_mesh.facets_begin (); it != output_mesh.facets_end (); ++it ) {

    auto circ = it->facet_begin ();

    do {
            // n_o[b++] = circ->vertex ()->id ();
    } while ( ++circ != it->facet_begin () );

}

[petrasvestartas]

Dear @afabri ,

Would it be possible to get a reply for this issue?

[sloriot]

1. There is the insert method.
2. You need to use Polyhedron_items_with_id_3 as template parameter of Polyhedron_3.

[afabri]

The class Point_set_3 has a method to insert() taking a point and a normal as arguments. As the documentation of the method says, you must add the normal property map once before.

[petrasvestartas]

Works perfectly, thank you @afabri and @sloriot (I consider this issue closed).

Can you also explain me why mesh properties are normally not indexed? Is it faster to have it that way? Before I was using OpenMesh as part of EPFL courses I was taking and they were not indexed either. Since my tutors did not answer this question, maybe you could?

Here is a working version from Rhino:

[sloriot]

If by index property you mean a property that can be accessed in constant time (using an index), then obviously it requires an index on the considered item. Point_set_3 properties is considering points as indices to properties associated to points (like coordinates) are in a vector. The same hold for Surface_mesh and OpenMesh. For Polyhedron, if you use the item with id mentioned above, only dynamic properties will be indexed.