CGAL users discussion list

[Muthukumaran Chandrasekaran <>]

Hi,

Please find my attempt at doing what Marc suggested. Code is attached.

1) I create a full triangulation of points with 1 finite cell.

2) Rewrote a variant of the barycentric_subdivide that now uses the full triangulation

3) Had to reimplement CGAL::barycenter() because it wasn't defined for Point_d

However, I am not sure if what I am doing is 100% correct.

I couldn't test this because I kept getting an error at compile time at line 48 (v = v + it - CGAL::ORIGIN;)

The error:
 "no match for ‘operator+’ (operand types are ‘Point_d {aka CGAL::Wrap::Point_d<CGAL::Epick_d<CGAL::Dynamic_dimension_tag> >}’ and ‘Point_d {aka CGAL::Wrap::Point_d<CGAL::Epick_d<CGAL::Dynamic_dimension_tag> >}’)"

Can someone please tell me why I am getting this error?

Also, can you please check the rest of the code to see if the implementation seems right?

Thanks

Muthu Chandrasekaran

PhD Candidate, THINC Lab

Department of Comp. Science

University of Georgia

On Wed, Feb 15, 2017 at 2:14 PM, Muthukumaran Chandrasekaran <> wrote:

Dear Marc,

Thank you for getting back to me so quickly!

You could build a triangulation of your points, which has just one finite cell, then apply a variant of barycentric_subdivide from the example to it (the example works on the TDS, you would want to work on a full triangulation, with points with coordinates), and iterate on the resulting finite cells.

Is there an example/demo/code of how I could build a full triangulation of points (which has 1 finite cell)?

I am trying to read up on all this as quickly as I can. Apparently there is a LOT that I don't know about geometry :P. Its getting a little overwhelming. Any direction would immensely help.

Sincerely,

Muthu

--

Muthukumaran Chandrasekaran
Department of Computer Science
University of Georgia

Phone: 7062472873

Website: http://muthuchandrasekaran.com

Email:
         

--

Muthukumaran Chandrasekaran
Department of Computer Science
University of Georgia

Phone: 7062472873

Website: http://muthuchandrasekaran.com

Email:
         

/*
 * barycentric_subdivision_test.cpp
 *
 *  Created on: Feb 17, 2017
 *      Author: AIUGA\mkran
 */

#define CGAL_EIGEN3_ENABLED
#if defined(__GNUC__) && defined(__GNUC_MINOR__) && (__GNUC__ <= 4) && (__GNUC_MINOR__ < 4)
#include <iostream>
int main()
{
  std::cerr << "NOTICE: This test requires G++ >= 4.4, and will not be compiled." << std::endl;
}
#else

#include <CGAL/Triangulation_data_structure.h>
#include <CGAL/internal/Combination_enumerator.h>
#include <CGAL/Epick_d.h>
#include <CGAL/point_generators_d.h>
#include <CGAL/Triangulation.h>
#include <CGAL/algorithm.h>
#include <CGAL/assertions.h>
#include <iostream>
#include <iterator>
#include <vector>
#include <CGAL/Cartesian_d.h>
using namespace std;

typedef float											FT;
typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag>		K;
typedef K::Point_d 										Point_d;
typedef CGAL::Triangulation<K>							Triangulation;

Point_d barycenter(int current_dimension, const vector<typename Triangulation::Vertex_handle> vertices){

	typedef pair<Point_d,FT> Point_Pair;
	vector<FT> dims;

	for(int i = 0; i < current_dimension; i++){
		dims.push_back(0.0);
	}
	Point_d v = Point_d(current_dimension,dims.begin(),dims.end());

	FT norm = 0;
	for( int i = 0; i < current_dimension; i++){
		Point_d it = vertices.at(i)->point();
		v = v + it - CGAL::ORIGIN;
		norm += 1;
	}

	vector<FT> dims2;
	for(int i=0;i<current_dimension;i++){
		dims2.push_back(v[i]/norm);
	}
	Point_d cd = Point_d(current_dimension,dims2.begin(),dims2.end());

	return cd;
}

/*
 * The main goal of this function is to find a full cell that
 * contains a given set of vertices |face_vertices|. Then, it
 * builds a corresponding |face|.
 */
void find_face_from_vertices(const Triangulation & T,
        const std::vector<typename Triangulation::Vertex_handle> & face_vertices,
        typename Triangulation::Face & face){

    typedef typename Triangulation::Vertex_handle           Vertex_handle;
    typedef typename Triangulation::Full_cell_handle        Full_cell_handle;
    typedef typename Triangulation::Full_cell::Vertex_handle_iterator Vertex_h_iterator;

    // get the dimension of the face we want to build
    std::size_t fdim(face_vertices.size() - 1);
    if( fdim <= 0) exit(-1);

    // find all full cells incident to the first vertex of |face|
    typedef std::vector<Full_cell_handle> Cells;
    Cells cells;
    std::back_insert_iterator<Cells> out(cells);
    T.incident_full_cells(face_vertices[0], out);

    // Iterate over the cells to find one which contains the face_vertices
    for( typename Cells::iterator cit = cells.begin(); cit != cells.end(); ++cit){
        // find if the cell *cit contains the Face |face|
        std::size_t i = 0;
        for( ; i <= fdim; ++i ) {
            Vertex_handle face_v(face_vertices[i]);
            bool found(false);
            Vertex_h_iterator vit = (*cit)->vertices_begin();
            for( ; vit != (*cit)->vertices_end(); ++vit ) {
                if( *vit == face_v ) {
                    found = true;
                    break;
                }
            }
            if( ! found )
                break;
        }
        if( i > fdim ) {// the full cell |*cit| contains |face|
            face.set_full_cell(*cit);
            for(size_t i = 0; i <= fdim; ++i){
              face.set_index(static_cast<int>(i),(*cit)->index(face_vertices[i]));
            }
            return;
        }
    }
    cerr << "Could not build a face from vertices"<<std::endl;
    CGAL_assertion(false);
}

/*
 * This function builds the barycentric subdivision of a single
 * finite full cell |fc| from a full triangulation of points with coordinates.
 */
void barycentric_subdivision(Triangulation & T, Triangulation::Full_cell_handle fc){
	typedef typename Triangulation::Vertex_handle Vertex_handle;
    typedef typename Triangulation::Face Face;
    const int dim = T.current_dimension();

    // First, read handles to the cell's vertices
    vector<Vertex_handle> vertices;
    vector<Vertex_handle> face_vertices;
    for( int i = 0; i <= dim; ++i ) vertices.push_back(fc->vertex(i));

    //First find the barycenter of the finite full cell
    Point_d bcenter = barycenter(dim,vertices);
    cout << bcenter << endl;

    // Then, subdivide the cell |fc| once by inserting one vertex
    // (and associate the vertex to the barycenter).
    T.insert_in_full_cell(bcenter,fc);

    // From now on, we can't use the variable |fc|...
    // Then, subdivide faces of |fc| in order of decreasing dimension
    for( int d = dim-1; d > 0; --d )
    {
        face_vertices.resize(d+1);
        // The following class
        // enumerates all (d+1)-tuple of the set {0, 1, ..., dim}
        CGAL::internal::Combination_enumerator combi(d+1, 0, dim);
        while( ! combi.end() )
        {
            for( int i = 0; i <= d; ++i )
                face_vertices[i] = vertices[combi[i]];
            // we need to find a face with face_vertices
            Face face(dim);
            find_face_from_vertices(T, face_vertices, face);

            //Given the vertices of face in dimension d, find its barycenter
            Point_d bcenter_face = barycenter(d,face_vertices);
            cout << bcenter_face << endl;

            //Associate this point to a vertex and insert this vertex into the face
            T.insert_in_face(bcenter_face,face);
            ++combi;
        }
    }
}

int main(){
    vector<Point_d > points_d;
    const int D = 5;

    //Initializing the points in the convex hull of the original simplex
	for(int kk = 0 ; kk<D; kk++){
		std::vector<float> dims;
		for(int kkk=0;kkk<D;kkk++){
			dims.push_back( kk==kkk ? 1.0 :0.0);
		}
		points_d.push_back(Point_d(dims.begin(),dims.end()));
	}

	//Create a full Triangulation and insert the points defined previously
	vector<Triangulation::Point> points = points_d;
    Triangulation t(D);                      // create triangulation
    CGAL_assertion(t.empty());
    t.insert(points.begin(), points.end());  // compute triangulation
    CGAL_assertion( t.is_valid() );

    //Sanity Check
    cout<<"[Before Subdivision] Sanity Check: "<<endl;
    std::cout<<t.number_of_vertices()<<std::endl;
    std::cout<<t.number_of_full_cells()<<std::endl;
    std::cout<<t.number_of_finite_full_cells()<<std::endl;

    //Perform barycentric subdivision using the full triangulation of points
    cout<<"Doing subdivision..."<<endl;

    Triangulation::Full_cell_handle fc;
    fc = t.finite_full_cells_begin()->vertex(0)->full_cell();

    barycentric_subdivision(t, fc);

    cout<<"Done subdivision..."<<endl;

    cout<<"[After Subdivision] Sanity Check: "<<endl;
    std::cout<<t.number_of_vertices()<<std::endl;
    std::cout<<t.number_of_full_cells()<<std::endl;
    std::cout<<t.number_of_finite_full_cells()<<std::endl;

    cout<<"Output the coordinates of "<<t.number_of_finite_full_cells()<<" finite full cells resulting from the subdivision:"<<endl;
    for(Triangulation::Finite_full_cell_iterator it=t.finite_full_cells_begin(); it!=t.finite_full_cells_end(); it++){
		for(int k=0;k<D;k++){
			cout<<it->vertex(k)->point()<<endl;
		}
	}

    return 0;
}
#endif