Getting the 3D facets of a 4D Triangulation

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Getting the 3D facets of a 4D Triangulation

GSBicalho
I apologize if the question is obvious, but I am a beginner to CGAL.

Given a 3d cube with vertices 1 to 8, I would obtain a shape made of 6 squares, each of which would define a face of the cube.

If I wanted render this cube, then it would improve my performance if I transformed all the square faces into triangles.

I would like to take this logic and increase one dimension on it. That means that, given the 16 vertices of a hypercube, I'd like to find the pyramids that divide its 3d faces.

I've looked at CGAL and it looks like it would be able to do what I want with its dD Triangulations, by first triangulating it into 4D pyramids and finding the 3d faces of those.

For testing, the code below does the triangulation with a 4d pyramid, a 5-cell.

---CODE---

#include <CGAL/Epick_d.h>
#include <CGAL/Triangulation.h>
#include <CGAL/algorithm.h>
#include <CGAL/assertions.h>
#include <CGAL/Kernel_d/Point_d.h>
#include <CGAL/Delaunay_triangulation.h>
#include <iostream>
#include <iterator>
#include <vector>
#include <fstream>

const int D = 4; // we work in Euclidean 4-space
typedef CGAL::Dimension_tag< D > Dim_tag;
typedef CGAL::Epick_d< Dim_tag >  Kernel;
typedef CGAL::Epick_d< Dim_tag >::Point_d Point_4;
typedef CGAL::Triangulation< Kernel > Triangulation;
typedef CGAL::Delaunay_triangulation<Kernel> DenTriangulation;

int main()
{
        std::vector<std::vector<double>> coords = {
                { 0.0, 0.0, 0.0, 0.0 },
                { 0.0, 0.0, 0.0, 1.0 },
                { 0.0, 0.0, 1.0, 0.0 },
                { 0.0, 1.0, 0.0, 0.0 },
                { 1.0, 0.0, 0.0, 0.0 },
        };

        std::vector<Point_4> points;
        for (int i = 0; i < coords.size(); i++) {
                Point_4 p = Point_4(D, coords[i].begin(), coords[i].end());
                points.push_back(p);
        }
       
        Triangulation t(D);                      // create triangulation
        CGAL_assertion(t.empty());
        t.insert(points.begin(), points.end());  // compute triangulation
        CGAL_assertion(t.is_valid());

        std::cout << "Tn " << t.number_of_finite_full_cells() << std::endl;

        int i = 0;
        typedef Triangulation::Full_cell_iterator Full_cell_iterator;
        typedef Triangulation::Facet Facet;

        for (Full_cell_iterator cit = t.full_cells_begin();
                cit != t.full_cells_end(); ++cit){
                if (!t.is_infinite(cit))
                        continue;
                Facet ft(cit, cit->index(t.infinite_vertex()));

                std::cout << "V " << std::endl;
                std::cout << ft.full_cell()->vertex(0)->point() << std::endl;
                std::cout << ft.full_cell()->vertex(1)->point() << std::endl;
                std::cout << ft.full_cell()->vertex(2)->point() << std::endl;
                std::cout << ft.full_cell()->vertex(3)->point() << std::endl;
                std::cout << ft.full_cell()->vertex(4)->point() << std::endl;

                ++i;// |ft| is a facet of the convex hull
        }
        std::cout << "There are " << i << " facets on the convex hull." << std::endl;
       
        i = 0;
        for (auto cit = t.finite_full_cells_begin();
        cit != t.finite_full_cells_end(); ++cit)
        {
        std::cout << "V " << std::endl;
        std::cout << cit->vertex(0)->point() << std::endl;
        std::cout << cit->vertex(1)->point() << std::endl;
        std::cout << cit->vertex(2)->point() << std::endl;
        std::cout << cit->vertex(3)->point() << std::endl;
        std::cout << cit->vertex(4)->point() << std::endl;

        ++i;// |ft| is a facet of the convex hull
        }
        std::cout << "There are " << i << " full_cells on the convex hull." << std::endl;

        return 0;
}

------

It produces the following output:

---OUTPUT---

Tn 1
V
4 0 0 1 0
4 0 0 0 0
4 1 0 0 0
4 0 0 1 0
4 0 1 0 0
V
4 0 0 1 0
4 0 0 0 1
4 0 0 0 0
4 0 0 1 0
4 0 1 0 0
V
4 0 0 1 0
4 1 0 0 0
4 0 0 0 1
4 0 0 1 0
4 0 1 0 0
V
4 0 0 1 0
4 0 0 0 1
4 1 0 0 0
4 0 0 1 0
4 0 0 0 0
V
4 0 0 1 0
4 0 0 0 1
4 1 0 0 0
4 0 0 0 0
4 0 1 0 0
There are 5 facets on the convex hull.
V
4 0 0 0 1
4 1 0 0 0
4 0 0 0 0
4 0 0 1 0
4 0 1 0 0
There are 1 full_cells on the convex hull.

------

As I understand, it should print the vertexes of 5 pyramids, which it indeed does, however these pyramids make no sense. They should have 4 vertexes, but have 5. On the first four, one of the vertex is repeated. On the fifth, it is the same vertexes of the original pyramid!

I believe that the error may be on the "ft.full_cell()", which I use to get the vertexes. This seems to get me the full_cell the facet belongs to, which would be wrong. However, I see not other way to access any information regarding the facet.

The number of full_cells appears to be right, but since the previous result, which I made following the tutorial, didn't work out, I am hesitant to trust that it would work correctly for larger numbers.

Am I going at this the wrong way? Am I not accessing the vertexes of the facet correctly? Or have I misundertood what the dDTriangulation does?
Reply | Threaded
Open this post in threaded view
|

Re: Getting the 3D facets of a 4D Triangulation

Marc Glisse
On Fri, 21 Apr 2017, GSBicalho wrote:

> I apologize if the question is obvious, but I am a beginner to CGAL.
>
> Given a 3d cube with vertices 1 to 8, I would obtain a shape made of 6
> squares, each of which would define a face of the cube.
>
> If I wanted render this cube, then it would improve my performance if I
> transformed all the square faces into triangles.
>
> I would like to take this logic and increase one dimension on it. That means
> that, given the 16 vertices of a hypercube, I'd like to find the pyramids
> that divide its 3d faces.
>
> I've looked at CGAL and it looks like it would be able to do what I want
> with its dD Triangulations, by first triangulating it into 4D pyramids and
> finding the 3d faces of those.
>
> For testing, the code below does the triangulation with a 4d pyramid, a
> 5-cell.
>
> ---CODE---
>
> #include <CGAL/Epick_d.h>
> #include <CGAL/Triangulation.h>
> #include <CGAL/algorithm.h>
> #include <CGAL/assertions.h>
> #include <CGAL/Kernel_d/Point_d.h>
> #include <CGAL/Delaunay_triangulation.h>
> #include <iostream>
> #include <iterator>
> #include <vector>
> #include <fstream>
>
> const int D = 4; // we work in Euclidean 4-space
> typedef CGAL::Dimension_tag< D > Dim_tag;
> typedef CGAL::Epick_d< Dim_tag >  Kernel;
> typedef CGAL::Epick_d< Dim_tag >::Point_d Point_4;
> typedef CGAL::Triangulation< Kernel > Triangulation;
> typedef CGAL::Delaunay_triangulation<Kernel> DenTriangulation;
>
> int main()
> {
> std::vector<std::vector&lt;double>> coords = {
> { 0.0, 0.0, 0.0, 0.0 },
> { 0.0, 0.0, 0.0, 1.0 },
> { 0.0, 0.0, 1.0, 0.0 },
> { 0.0, 1.0, 0.0, 0.0 },
> { 1.0, 0.0, 0.0, 0.0 },
> };
>
> std::vector<Point_4> points;
> for (int i = 0; i < coords.size(); i++) {
> Point_4 p = Point_4(D, coords[i].begin(), coords[i].end());
> points.push_back(p);
> }
>
> Triangulation t(D);                      // create triangulation
> CGAL_assertion(t.empty());
> t.insert(points.begin(), points.end());  // compute triangulation
> CGAL_assertion(t.is_valid());
>
> std::cout << "Tn " << t.number_of_finite_full_cells() << std::endl;
>
> int i = 0;
> typedef Triangulation::Full_cell_iterator Full_cell_iterator;
> typedef Triangulation::Facet Facet;
>
> for (Full_cell_iterator cit = t.full_cells_begin();
> cit != t.full_cells_end(); ++cit){
> if (!t.is_infinite(cit))
> continue;
> Facet ft(cit, cit->index(t.infinite_vertex()));
>
> std::cout << "V " << std::endl;
> std::cout << ft.full_cell()->vertex(0)->point() << std::endl;
> std::cout << ft.full_cell()->vertex(1)->point() << std::endl;
> std::cout << ft.full_cell()->vertex(2)->point() << std::endl;
> std::cout << ft.full_cell()->vertex(3)->point() << std::endl;
> std::cout << ft.full_cell()->vertex(4)->point() << std::endl;
>
> ++i;// |ft| is a facet of the convex hull
> }
> std::cout << "There are " << i << " facets on the convex hull." <<
> std::endl;
>
> i = 0;
> for (auto cit = t.finite_full_cells_begin();
> cit != t.finite_full_cells_end(); ++cit)
> {
> std::cout << "V " << std::endl;
> std::cout << cit->vertex(0)->point() << std::endl;
> std::cout << cit->vertex(1)->point() << std::endl;
> std::cout << cit->vertex(2)->point() << std::endl;
> std::cout << cit->vertex(3)->point() << std::endl;
> std::cout << cit->vertex(4)->point() << std::endl;
>
> ++i;// |ft| is a facet of the convex hull
> }
> std::cout << "There are " << i << " full_cells on the convex hull." <<
> std::endl;
>
> return 0;
> }
>
> ------
>
> It produces the following output:
>
> ---OUTPUT---
>
> Tn 1
> V
> 4 0 0 1 0
> 4 0 0 0 0
> 4 1 0 0 0
> 4 0 0 1 0
> 4 0 1 0 0
> V
> 4 0 0 1 0
> 4 0 0 0 1
> 4 0 0 0 0
> 4 0 0 1 0
> 4 0 1 0 0
> V
> 4 0 0 1 0
> 4 1 0 0 0
> 4 0 0 0 1
> 4 0 0 1 0
> 4 0 1 0 0
> V
> 4 0 0 1 0
> 4 0 0 0 1
> 4 1 0 0 0
> 4 0 0 1 0
> 4 0 0 0 0
> V
> 4 0 0 1 0
> 4 0 0 0 1
> 4 1 0 0 0
> 4 0 0 0 0
> 4 0 1 0 0
> There are 5 facets on the convex hull.
> V
> 4 0 0 0 1
> 4 1 0 0 0
> 4 0 0 0 0
> 4 0 0 1 0
> 4 0 1 0 0
> There are 1 full_cells on the convex hull.
>
> ------
>
> As I understand, it should print the vertexes of 5 pyramids, which it indeed
> does, however these pyramids make no sense. They should have 4 vertexes, but
> have 5. On the first four, one of the vertex is repeated. On the fifth, it
> is the same vertexes of the original pyramid!
>
> I believe that the error may be on the "ft.full_cell()", which I use to get
> the vertexes. This seems to get me the full_cell the facet belongs to, which
> would be wrong. However, I see not other way to access any information
> regarding the facet.
>
> The number of full_cells appears to be right, but since the previous result,
> which I made following the tutorial, didn't work out, I am hesitant to trust
> that it would work correctly for larger numbers.
>
> Am I going at this the wrong way? Am I not accessing the vertexes of the
> facet correctly? Or have I misundertood what the dDTriangulation does?

When you ask for a full cell containing your facet, you will get back the
infinite full cell that was used to define the facet. One of its vertices
is the infinite vertex, and asking for the corresponding point returns
random crap. If I understand correctly, for each infinite cell, you want
to print its finite vertices, so something like

for(int i=0;i<5;++i){
   Triangulation::Vertex_const_handle vh=cit->vertex(i);
   if(!t.is_infinite(vh))
     std::cout << vh->point() << std::endl;
}

without needing to create a facet. If you do have a facet, you want to use
index_of_covertex to determine which vertex of the full cell you should
ignore. (it seems that the infinite vertex has index 0, but I don't
remember if that can be relied on)


--
Marc Glisse

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss