CGAL users discussion list

["Sebastien Loriot (GeometryFactory)" <>]

Several problems:
*you need to use a vertex with info in the TDS of the triangulation.
*the function build_triangulation_with_indices must be called on a 
triangulation.
*you must then build the alpha_shape from the triangulation.

S.

MrVH wrote:
> I'm trying to understand your mapping code and to modify it. But I've got
> some compilation errors and I don't catch some lines. I really have to use
> this code cause to avoid to recode AlphaShape Weighted from scratch in
> Python :(
> I don't see where the Delaunay function is launched. 
>
> In main I'm running:
> Alpha_shape_3  as;
> build_triangulation_with_indices(lwp.begin(),lwp.end(),as);
>
> The build_triangulation_with_indices function take now an Alpha_Shape_3
> object:
> //function inserting points into a triangulation
> //and setting the info field to the order in the input list.
> template <class Alpha_Shape_3,class Point_iterator>
> void
> build_triangulation_with_indices(
>         Point_iterator begin,
>         Point_iterator end,
>         Alpha_Shape_3& T)
> {
> // The vector is now a Weighted_point vector cause it is regular
> triangulation
>   std::vector<std::pair<const typename Gt::Weighted_point*,int> > points;
>   int index=-1;
>   for (Point_iterator it=begin;it!=end;++it){
>     points.push_back(std::make_pair(&(*it),++index));
>   }
>   std::random_shuffle (points.begin(), points.end()); // shuffle point but
> why ???
>   spatial_sort (points.begin(),points.end(),
> Traits_for_spatial_sort<Alpha_Shape_3>()); // build a structure struct
> Traits_for_spatial_sort:public Alpha_Shape_3::Geom_traits
>
> // from here I do not understand. Where Trigulation is launched ???
>   typename Alpha_Shape_3::Cell_handle hint;
>   for (typename std::vector< std::pair<const typename
> Gt::Weighted_point*,int> >::const_iterator
>         p = points.begin();p != points.end(); ++p)
>   {
>     typename Alpha_Shape_3::Locate_type lt;
>     typename Alpha_Shape_3::Cell_handle c;
>     int li, lj;
>     c = T.locate (*(p->first), lt, li, lj, hint);
>
>     typename Alpha_Shape_3::Vertex_handle v = T.insert (*(p->first), lt, c,
> li, lj);
>     if ( v==typename Alpha_Shape_3::Vertex_handle() )
>       hint=c;
>     else{
>       v->info() = p->second ;
>       hint=v->cell();
>     }
>   }
> }  
>
> the structure
> //typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
>
> typedef CGAL::Triangulation_vertex_base_with_info_3<int, K> Wb; // avoid
> conflict so Wb instead of Vb cause casting is different than the previous
> <int, K>
> typedef CGAL::Triangulation_data_structure_3<Wb>                    Tds2;
> typedef CGAL::Delaunay_triangulation_3<K, Tds2>                     
> Delaunay; // avoir conflict
>
> typedef Delaunay::Point   Point;
>
> //spatial sort traits to use with a pair of point pointers and integer.
> template<class Alpha_Shape_3>
> struct Traits_for_spatial_sort:public Alpha_Shape_3::Geom_traits{
>   typedef typename Alpha_Shape_3::Geom_traits Gt;
>   typedef std::pair<const typename Gt::Weighted_point*,int> Point_3;
>   
>   struct Less_x_3{
>     bool operator()(const Point_3& p,const Point_3& q) const {
>       return typename Gt::Less_x_3()(*(p.first),*(q.first));
>     }
>   };
>
>   struct Less_y_3{
>     bool operator()(const Point_3& p,const Point_3& q) const {
>       return typename Gt::Less_y_3()(*(p.first),*(q.first));
>     }
>   };
>
>   struct Less_z_3{
>     bool operator()(const Point_3& p,const Point_3& q) const {
>       return typename Gt::Less_z_3()(*(p.first),*(q.first));
>     }
>   };  
>   
>   Less_x_3      less_x_3_object () const {return Less_x_3();}
>   Less_y_3      less_y_3_object () const {return Less_y_3();}
>   Less_z_3      less_z_3_object () const {return Less_z_3();}
> };
>
>
> Thank you for your help and your explanations,
>
>
>
>
> ############################################################
> The full code:
> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
> #include <CGAL/Regular_triangulation_euclidean_traits_3.h>
> #include <CGAL/Regular_triangulation_3.h>
> #include <CGAL/Alpha_shape_3.h>
> #include <fstream>
> #include <iostream>
> #include <cstdio>
> #include <cstring>
> #include <string>
> #include <map>
> #include <set>
> #include <vector>
> #include <cmath>
> #include <stdio.h>
> #include <stdlib.h>
>
>
> //#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
> #include <CGAL/Delaunay_triangulation_3.h>
> #include <CGAL/Triangulation_vertex_base_with_info_3.h>
> //#include <iostream>
>
>
> typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
> typedef CGAL::Regular_triangulation_euclidean_traits_3<K> Gt;
> typedef CGAL::Alpha_shape_vertex_base_3<Gt>         Vb;
> typedef CGAL::Alpha_shape_cell_base_3<Gt>           Fb;
> typedef CGAL::Triangulation_data_structure_3<Vb,Fb> Tds;
> typedef CGAL::Regular_triangulation_3<Gt,Tds>       Triangulation_3;
> typedef CGAL::Alpha_shape_3<Triangulation_3>        Alpha_shape_3;
>
> typedef Alpha_shape_3::Cell_handle          Cell_handle;
> typedef Alpha_shape_3::Vertex_handle        Vertex_handle;
> typedef Alpha_shape_3::Facet                Facet;
> typedef Alpha_shape_3::Edge                 Edge;
>
> typedef Gt::Weighted_point                  Weighted_point;
> typedef Gt::Bare_point                      Bare_point;
>
>
>
>
>
>
>
>
>
>
>
>
> //typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
>
> typedef CGAL::Triangulation_vertex_base_with_info_3<int, K> Wb;
> typedef CGAL::Triangulation_data_structure_3<Wb>                    Tds2;
> typedef CGAL::Delaunay_triangulation_3<K, Tds2>                     
> Delaunay;
>
> typedef Delaunay::Point   Point;
>
> //spatial sort traits to use with a pair of point pointers and integer.
> template<class Alpha_Shape_3>
> struct Traits_for_spatial_sort:public Alpha_Shape_3::Geom_traits{
>   typedef typename Alpha_Shape_3::Geom_traits Gt;
>   typedef std::pair<const typename Gt::Weighted_point*,int> Point_3;
>   
>   struct Less_x_3{
>     bool operator()(const Point_3& p,const Point_3& q) const {
>       return typename Gt::Less_x_3()(*(p.first),*(q.first));
>     }
>   };
>
>   struct Less_y_3{
>     bool operator()(const Point_3& p,const Point_3& q) const {
>       return typename Gt::Less_y_3()(*(p.first),*(q.first));
>     }
>   };
>
>   struct Less_z_3{
>     bool operator()(const Point_3& p,const Point_3& q) const {
>       return typename Gt::Less_z_3()(*(p.first),*(q.first));
>     }
>   };  
>   
>   Less_x_3      less_x_3_object () const {return Less_x_3();}
>   Less_y_3      less_y_3_object () const {return Less_y_3();}
>   Less_z_3      less_z_3_object () const {return Less_z_3();}
> };
>
>
> //function inserting points into a triangulation
> //and setting the info field to the order in the input list.
> template <class Alpha_Shape_3,class Point_iterator>
> void
> build_triangulation_with_indices(
>         Point_iterator begin,
>         Point_iterator end,
>         Alpha_Shape_3& T)
> {
>   std::vector<std::pair<const typename Gt::Weighted_point*,int> > points;
>   int index=-1;
>   for (Point_iterator it=begin;it!=end;++it){
>     points.push_back(std::make_pair(&(*it),++index));
>   }
>   std::random_shuffle (points.begin(), points.end());
>   spatial_sort (points.begin(),points.end(),
> Traits_for_spatial_sort<Alpha_Shape_3>());
>
>   typename Alpha_Shape_3::Cell_handle hint;
>   for (typename std::vector< std::pair<const typename
> Gt::Weighted_point*,int> >::const_iterator
>         p = points.begin();p != points.end(); ++p)
>   {
>     typename Alpha_Shape_3::Locate_type lt;
>     typename Alpha_Shape_3::Cell_handle c;
>     int li, lj;
>     c = T.locate (*(p->first), lt, li, lj, hint);
>
>     typename Alpha_Shape_3::Vertex_handle v = T.insert (*(p->first), lt, c,
> li, lj);
>     if ( v==typename Alpha_Shape_3::Vertex_handle() )
>       hint=c;
>     else{
>       v->info() = p->second ;
>       hint=v->cell();
>     }
>   }
> }  
>
>
>
>
>
>
>
>
>
>
>
>
>
> int main(int argc, char *argv[]) 
> {
>  
>   std::list<Weighted_point> lwp;
>      
>   /*Get arguments*/
>   char *filename = argv[1];    // Format of file %8.3f%8.3f%8.3f%8.3f, x, y,
> z, weight 
>   float sizeas = atof(argv[2]);  // size of AlphaSphere
>   /*----------------*/
>   
>     
>   std::cout << "# File" << filename << std::endl;
>   
>   /*Read the file*/
>   std::ifstream fcoords;
>   fcoords.open(filename);
>   std::string line;
>   while (std::getline(fcoords, line)){
>     
>         /* Get coordinates*/
>         std::string type(line.substr(0,1));
>         if (type == "#" || type == "\n"){continue;} // comment line
>         
>         std::string xstr(line.substr(0,8));
>         std::string ystr(line.substr(8,8));
>         std::string zstr(line.substr(16,8));
>         std::string weightstr(line.substr(24,8));
>         //std::cout << xstr << ystr << zstr << weightstr << std::endl;
>
>         float x = atof(xstr.c_str());
>         float y = atof(ystr.c_str());
>         float z = atof(zstr.c_str());
>         float weight = atof(weightstr.c_str());
>         
>         /* fill each point with its associate weight */
>         lwp.push_back(Weighted_point(Bare_point( x, y, z), weight));
>        
>   }
>   fcoords.close();
>   
>   
>   std::cout << "# read " << lwp.size() << " atoms" << std::endl;
>   
>   /* Launch AlphaShape*/
>   //Alpha_shape_3  as(lwp.begin(), lwp.end(), sizeas,
> Alpha_shape_3::GENERAL);
>   Alpha_shape_3  as;
>   build_triangulation_with_indices(lwp.begin(),lwp.end(),as);
>   
>   
>   /* Create empty storage*/
>   std::list<Cell_handle>       cells;
>   std::list<Facet>       facets;
>   std::list<Edge>        edges;
>   
>   /*Get information 
>     Cell -> Tetraedron
>     Facet -> Triangle
>     Edge -> segment
>   */
>  
> as.get_alpha_shape_cells(std::back_inserter(cells),Alpha_shape_3::INTERIOR);
>  
> as.get_alpha_shape_facets(std::back_inserter(facets),Alpha_shape_3::REGULAR);
>  
> as.get_alpha_shape_facets(std::back_inserter(facets),Alpha_shape_3::SINGULAR);
>  
> as.get_alpha_shape_edges(std::back_inserter(edges),Alpha_shape_3::SINGULAR);
>   
>   
>   /* Enumeration */
>   std::map<Alpha_shape_3::Vertex*,unsigned int> index;
>   
>   
>   unsigned int i=0;
>   for (Alpha_shape_3::Finite_vertices_iterator
> vit=as.finite_vertices_begin();vit!=as.finite_vertices_end();++vit)
>     index.insert(std::make_pair(&(*vit),i++));
>   
>   /* Tetrahedron*/
>   for (std::list<Cell_handle>::iterator
> cell=cells.begin();cell!=cells.end();++cell){
>     std::cout << "cell ";
>     //for (int k=0;k<4;++k) std::cout << index[ &( *((*cell)->vertex(k)) ) ]
> << " ";
>     //std::cout <<std::endl;
>     for (int k=0;k<4;++k) std::cout <<  (*cell)->vertex(k)->point() << " ";
>     std::cout <<std::endl; 
>   }
>   
>   /* Triangle */
>   for (std::list<Facet>::iterator
> facet=facets.begin();facet!=facets.end();++facet){
>     std::cout << "facet ";
>     for (int k=0;k<3;++k)
>       std::cout << (facet->first->vertex( (facet->second+k)%4 ))->point() <<
> " ";
>       //std::cout << index[ &(*(facet->first->vertex( (facet->second+k)%4
> )))] << " ";
>     std::cout <<std::endl;
>   }
>
>   /* Edge */
>   for (std::list<Edge>::iterator
> edge=edges.begin();edge!=edges.end();++edge)
>     {
>     std::cout << "edge ";
>     std::cout << (edge->first->vertex( edge->second ))->point() <<" ";
>     std::cout << (edge->first->vertex( edge->third ))->point() << std::endl;
>     //std::cout << index[ &(*(edge->first->vertex( edge->second )))] <<" ";
>     //std::cout << index[ &(*(edge->first->vertex( edge->third )))]
> <<std::endl;
>     }
>   
>   std::cout << "# The 0-shape has : " << std::endl;
>   std::cout << "# " << cells.size() << " interior tetrahedra" << std::endl;
>   std::cout << "# " << facets.size() << " boundary facets" << std::endl;
>   std::cout << "# " << edges.size()  << " singular edges" << std::endl;
>   return 1;
>  
>
> }

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_euclidean_traits_3.h>
#include <CGAL/Regular_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <fstream>
#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <map>
#include <set>
#include <vector>
#include <cmath>
#include <stdio.h>
#include <stdlib.h>


//#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
//#include <iostream>


typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Regular_triangulation_euclidean_traits_3<K> Gt;
typedef CGAL::Alpha_shape_vertex_base_3<Gt>         Vbase;
typedef CGAL::Triangulation_vertex_base_with_info_3<int, Gt,Vbase> Vb;
typedef CGAL::Alpha_shape_cell_base_3<Gt>           Fb;
typedef CGAL::Triangulation_data_structure_3<Vb,Fb> Tds;
typedef CGAL::Regular_triangulation_3<Gt,Tds>       Triangulation_3;
typedef CGAL::Alpha_shape_3<Triangulation_3>        Alpha_shape_3;

typedef Alpha_shape_3::Cell_handle          Cell_handle;
typedef Alpha_shape_3::Vertex_handle        Vertex_handle;
typedef Alpha_shape_3::Facet                Facet;
typedef Alpha_shape_3::Edge                 Edge;

typedef Gt::Weighted_point                  Weighted_point;
typedef Gt::Bare_point                      Bare_point;



//spatial sort traits to use with a pair of point pointers and integer.
template<class Alpha_Shape_3>
struct Traits_for_spatial_sort:public Alpha_Shape_3::Geom_traits{
  typedef typename Alpha_Shape_3::Geom_traits Gt;
  typedef std::pair<const typename Gt::Weighted_point*,int> Point_3;
  
  struct Less_x_3{
    bool operator()(const Point_3& p,const Point_3& q) const {
      return typename Gt::Less_x_3()(*(p.first),*(q.first));
    }
  };

  struct Less_y_3{
    bool operator()(const Point_3& p,const Point_3& q) const {
      return typename Gt::Less_y_3()(*(p.first),*(q.first));
    }
  };

  struct Less_z_3{
    bool operator()(const Point_3& p,const Point_3& q) const {
      return typename Gt::Less_z_3()(*(p.first),*(q.first));
    }
  };  
  
  Less_x_3      less_x_3_object () const {return Less_x_3();}
  Less_y_3 	less_y_3_object () const {return Less_y_3();}
  Less_z_3 	less_z_3_object () const {return Less_z_3();}
};


//function inserting points into a triangulation
//and setting the info field to the order in the input list.
template <class Triangulation,class Point_iterator>
void
build_triangulation_with_indices(
	Point_iterator begin,
	Point_iterator end,
	Triangulation& T)
{
  std::vector<std::pair<const typename Gt::Weighted_point*,int> > points;
  int index=-1;
  for (Point_iterator it=begin;it!=end;++it){
    points.push_back(std::make_pair(&(*it),++index));
  }
  std::random_shuffle (points.begin(), points.end());
  spatial_sort (points.begin(),points.end(),Traits_for_spatial_sort<Triangulation>());

  typename Triangulation::Cell_handle hint;
  for (typename std::vector< std::pair<const typename Gt::Weighted_point*,int> >::const_iterator
        p = points.begin();p != points.end(); ++p)
  {
    typename Triangulation::Locate_type lt;
    typename Triangulation::Cell_handle c;
    int li, lj;
    c = T.locate (*(p->first), lt, li, lj, hint);

    typename Triangulation::Vertex_handle v = T.insert (*(p->first), lt, c,li, lj);
    if ( v==typename Triangulation::Vertex_handle() )
      hint=c;
    else{
      v->info() = p->second ;
      hint=v->cell();
    }
  }
}  

int main(int argc, char *argv[]) 
{
 
  std::list<Weighted_point> lwp;
     
  /*Get arguments*/
  char *filename = argv[1];    // Format of file %8.3f%8.3f%8.3f%8.3f, x, y, z, weight 
  /*----------------*/
  
    
  std::cout << "# File" << filename << std::endl;
  
  /*Read the file*/
  std::ifstream fcoords;
  fcoords.open(filename);
  std::string line;
  while (std::getline(fcoords, line)){
    
        /* Get coordinates*/
        std::string type(line.substr(0,1));
        if (type == "#" || type == "\n"){continue;} // comment line
	
        std::string xstr(line.substr(0,8));
        std::string ystr(line.substr(8,8));
        std::string zstr(line.substr(16,8));
	std::string weightstr(line.substr(24,8));
	//std::cout << xstr << ystr << zstr << weightstr << std::endl;

        float x = atof(xstr.c_str());
        float y = atof(ystr.c_str());
        float z = atof(zstr.c_str());
	float weight = atof(weightstr.c_str());
	
	/* fill each point with its associate weight */
        lwp.push_back(Weighted_point(Bare_point( x, y, z), weight));
       
  }
//  while (std::getline(fcoords, line)){
  	
	/* Get coordinates*/
//	std::string type(line.substr(0,4));
//	if (type != "ATOM"){continue;}
	
//	std::string xstr(line.substr(30,8));
//	std::string ystr(line.substr(38,8));
//	std::string zstr(line.substr(46,8));
	
//	float x = atof(xstr.c_str());
//	float y = atof(ystr.c_str());
//	float z = atof(zstr.c_str());
	
//	lwp.push_back(Weighted_point(Bare_point( x, y, z), 1.80));
	
//  }  
  fcoords.close();
  
  
  std::cout << "# read " << lwp.size() << " atoms" << std::endl;
  
  /* Launch AlphaShape*/
  Triangulation_3 T;
  build_triangulation_with_indices(lwp.begin(),lwp.end(),T);
  std::cout << T.number_of_vertices()<< std::endl;
  Alpha_shape_3 as(T,0.,Alpha_shape_3::GENERAL);
  
  
  /* Create empty storage*/
  std::list<Cell_handle>       cells;
  std::list<Facet>       facets;
  std::list<Edge>        edges;
  
  /*Get information 
    Cell -> Tetraedron
    Facet -> Triangle
    Edge -> segment
  */
  as.get_alpha_shape_cells(std::back_inserter(cells),Alpha_shape_3::INTERIOR);
  as.get_alpha_shape_facets(std::back_inserter(facets),Alpha_shape_3::REGULAR);
  as.get_alpha_shape_facets(std::back_inserter(facets),Alpha_shape_3::SINGULAR);
  as.get_alpha_shape_edges(std::back_inserter(edges),Alpha_shape_3::SINGULAR);
  
  
  /* Enumeration */
  std::map<Alpha_shape_3::Vertex*,unsigned int> index;
  
  
  unsigned int i=0;
  for (Alpha_shape_3::Finite_vertices_iterator vit=as.finite_vertices_begin();vit!=as.finite_vertices_end();++vit)
    index.insert(std::make_pair(&(*vit),i++));
  
  /* Tetrahedron*/
  for (std::list<Cell_handle>::iterator cell=cells.begin();cell!=cells.end();++cell){
    std::cout << "cell ";
    for (int k=0;k<4;++k) std::cout <<  (*cell)->vertex(k)->info() << " ";
    std::cout <<std::endl; 
  }
  
  /* Triangle */
  for (std::list<Facet>::iterator facet=facets.begin();facet!=facets.end();++facet){
    std::cout << "facet ";
    for (int k=0;k<3;++k)
      std::cout << (facet->first->vertex( (facet->second+k)%4 ))->info() << " ";
    std::cout <<std::endl;
  }

  /* Edge */
  for (std::list<Edge>::iterator edge=edges.begin();edge!=edges.end();++edge)
  {
    std::cout << "edge ";
    std::cout << (edge->first->vertex( edge->second ))->info() <<" ";
    std::cout << (edge->first->vertex( edge->third ))->info() << std::endl;
  }
  
  std::cout << "# The 0-shape has : " << std::endl;
  std::cout << "# " << cells.size() << " interior tetrahedra" << std::endl;
  std::cout << "# " << facets.size() << " boundary facets" << std::endl;
  std::cout << "# " << edges.size()  << " singular edges" << std::endl;
  return 0;
}