CGAL users discussion list

["Sebastien Loriot (GeometryFactory)" <>]

Here are two versions compiling and working,
but with different output!

main.cpp: since the constructor with range of alpha shape is used (which
is calling then the constructor with range of the Delaunay triangulation
the order in which the points are inserted in not the same than the one
in the sequence). This is why I introduced a map which
If you want to keep track of the original insertion order, look at this
thread:
https://lists-sop.inria.fr/sympa/arc/cgal-discuss/2010-01/msg00003.html

main-esbtl.cpp:This example uses the ESBTL library 
(http://esbtl.sourceforge.net/) to read your input pdbfile.
It uses a kernel where points are atoms and thus you can access their
atom serial number and print it.

S.


MrVH wrote:
> Thank you for your anwser. 
> It's really unclear for me.
> it doesn't work :( To answer you about point, I just want to now the
> belonging of my point to tetraedra/facets/Edge
> input
> coord
> 1 x y z w
> 2 x y z w
> ....
>
> output
> Results (results of Alpha Shape)
> edges 1 2
> edges 4 24
> ...
> Tetrehedron 1 4 6 2
> ...
> Triangles 78 4 3
> ...
>
>
>
> 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;
>
>
> /*TEST AREA*/
> typedefTriple<Cell_handle, int, int> Edge;  // does not compile
> typedef std::pair<Cell_handle, int> Facet; 
>
>
> int main()
> {
>   
>   std::list<Weighted_point> lwp;
>   
>   ifstream fpdb;
>   fpdb.open("1GXD.pdb");
>   string line;
>   while (getline(fpdb, line)){
>         
>         /* Get coordinates*/
>         string type(line.substr(0,4));
>         if (type != "ATOM"){continue;}
>         
>         string xstr(line.substr(30,8));
>         string ystr(line.substr(38,8));
>         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));
>         
>   }
>   fpdb.close();
>   
>   
>   //build alpha_shape  in GENERAL mode and set alpha=0
>   Alpha_shape_3  as(lwp.begin(), lwp.end(), 0, Alpha_shape_3::GENERAL);
>   //explore the 0-shape - It is dual to the boundary of the union.
>   //CGAL::Triangulation_data_structure_3::Cell_handle cells;
>   std::list<Cell_handle>       cells;
>   std::list<Facet>       facets;
>   std::list<Edge>        edges;
>   //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);
>   
>   
>   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++);  }
>   
>   for (int k=0;k<4;++k) std::cout << index[ &( *(cells->vertex(k)) ) ];
>   std::cout <<std::endl; 
>
>   for (k=0;k<3;++k)
>    std::cout << index[ &(*(facet.first->vertex( (facets.second+k)%4 )))];
>   std::cout <<std::endl; 
>
>   std::cout << index[ &(*(edges.first->vertex( edges.second )))] <<" ";
>   std::cout << index[ &(*(edges.first->vertex( edges.third )))];
>   std::cout <<std::endl; 
>   
>   /*
>   std::list<Edge>::iterator fbeg = edges.begin();
>   std::list<Edge>::iterator fend = edges.end();
>   
>   for (fbeg; fbeg != fend; fbeg++){
>     std::cout << fbeg.first << std::endl;
>
>
>
>     }*/
>
>   
>   /*
> for (Alpha_Shape_Edges_Iterator edge_alpha_it=edges.alpha_begin();
> edge_alpha_it!=edges.alpha_end(); edge_alpha_it++)
>     {
>       Alpha_Edge& edge = *edge_alpha_it;
>       
>       const Point_2& p1 =
> edge.first->vertex(edge.first->ccw(edge.second))->point();
>       const Point_2& p2 = edge.first->vertex(edge.first->cw
> (edge.second))->point(); 
>       std::cout << p1 << p2 << endl;
>       }*/
>
> //  std::ifstream iFileT("output",std::ios::in);
>   //iFileT >> as;
>   
>   
>   std::cout << "TOTO" << std::endl;
>   
>   //std::cout << cells.neighbor(1) << 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>


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;

int main()
{
  
  std::list<Weighted_point> lwp;
  
  std::ifstream fpdb;
  fpdb.open("1GXD.pdb");
  std::string line;
  while (std::getline(fpdb, 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));
	
  }
  fpdb.close();
  
  std::cout << "read " << lwp.size() << " atoms" << std::endl;
  //build alpha_shape  in GENERAL mode and set alpha=0
  Alpha_shape_3  as(lwp.begin(), lwp.end(), 0, Alpha_shape_3::GENERAL);
  //explore the 0-shape - It is dual to the boundary of the union.
  std::list<Cell_handle>       cells;
  std::list<Facet>       facets;
  std::list<Edge>        edges;
  
  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);
  
  
  std::map<Alpha_shape_3::Vertex*,unsigned int> index;
  
  //set index of each vertex
  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++));
  
  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 (std::list<Facet>::iterator facet=facets.begin();facet!=facets.end();++facet){
    std::cout << "facet ";
    for (int k=0;k<3;++k)
      std::cout << index[ &(*(facet->first->vertex( (facet->second+k)%4 )))] << " ";
    std::cout <<std::endl;
  }

  
  for (std::list<Edge>::iterator edge=edges.begin();edge!=edges.end();++edge)
  {
    std::cout << "edge ";
    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>

//ESBTL include
#include <ESBTL/default.h>
#include <ESBTL/weighted_atom_iterator.h>
#include <ESBTL/CGAL/EPIC_kernel_with_atom.h>
#include <ESBTL/atom_classifier.h>

#include <fstream>
#include <iostream>

//Kernel from ESBTL
typedef ESBTL::CGAL::EPIC_kernel_with_atom                   K;
//CGAL stuff
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;

//ESBTL stuff
typedef ESBTL::CGAL::Default_system                                                 System;
typedef ESBTL::Generic_classifier<ESBTL::Radius_of_atom<double,System::Atom> >      T_Atom_classifier;
typedef ESBTL::Accept_none_occupancy_policy<ESBTL::PDB::Line_format<> >             Accept_none_occupancy_policy;
typedef ESBTL::Weighted_atom_iterator<System::Model,
                                      CGAL::Weighted_point<K::Point_3,double>,
                                      ESBTL::Weight_of_atoms<T_Atom_classifier> >   Weighted_atom_iterator;

int main(int argc,char** argv)
{
  if (argc != 2 ){
    std::cerr << "Please provide a filename"  << std::endl;  
    return EXIT_FAILURE;
  }
  
  ESBTL::PDB_line_selector_two_systems sel;
  
  std::vector<System> systems;
  ESBTL::All_atom_system_builder<System> builder(systems,sel.max_nb_systems());
  T_Atom_classifier atom_classifier;
    //read pdb file
  ESBTL::read_a_pdb_file(argv[1],sel,builder,Accept_none_occupancy_policy());
  if ( systems.empty() || systems[0].has_no_model() ){
    std::cerr << "No atoms found" << std::endl;
    return EXIT_FAILURE;
  }  
  
  
  //Consider only the first model of the first system
  const System::Model& model=*systems[0].models_begin();
  

  //build alpha_shape  in GENERAL mode and set alpha=0
  Alpha_shape_3  as(  Weighted_atom_iterator(model.atoms_begin(),&atom_classifier),
                      Weighted_atom_iterator(model.atoms_end(),&atom_classifier),
                      0, Alpha_shape_3::GENERAL );
  //explore the 0-shape - It is dual to the boundary of the union.
  std::list<Cell_handle>        cells;
  std::list<Facet>              facets;
  std::list<Edge>               edges;
  
  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);
  
  
  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)->point().atom_serial_number() << " ";
    std::cout <<std::endl; 
  }

  
  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().atom_serial_number() << " ";
    std::cout <<std::endl;
  }

  
  for (std::list<Edge>::iterator edge=edges.begin();edge!=edges.end();++edge)
  {
    std::cout << "edge ";
    std::cout << edge->first->vertex( edge->second )->point().atom_serial_number() <<" ";
    std::cout << edge->first->vertex( edge->third )->point().atom_serial_number() <<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;
  

}
# CMAKE generated file: DO NOT EDIT!
# Generated by "Unix Makefiles" Generator, CMake Version 2.8

# Default target executed when no arguments are given to make.
default_target: all
.PHONY : default_target

#=============================================================================
# Special targets provided by cmake.

# Disable implicit rules so canoncical targets will work.
.SUFFIXES:

# Remove some rules from gmake that .SUFFIXES does not remove.
SUFFIXES =

.SUFFIXES: .hpux_make_needs_suffix_list

# Suppress display of executed commands.
$(VERBOSE).SILENT:

# A target that is always out of date.
cmake_force:
.PHONY : cmake_force

#=============================================================================
# Set environment variables for the build.

# The shell in which to execute make rules.
SHELL = /bin/sh

# The CMake executable.
CMAKE_COMMAND = /usr/bin/cmake

# The command to remove a file.
RM = /usr/bin/cmake -E remove -f

# The program to use to edit the cache.
CMAKE_EDIT_COMMAND = /usr/bin/ccmake

# The top-level source directory on which CMake was run.
CMAKE_SOURCE_DIR = /tmp/as

# The top-level build directory on which CMake was run.
CMAKE_BINARY_DIR = /tmp/as

#=============================================================================
# Targets provided globally by CMake.

# Special rule for the target edit_cache
edit_cache:
        @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan 
"Running CMake cache editor..."
        /usr/bin/ccmake -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
.PHONY : edit_cache

# Special rule for the target edit_cache
edit_cache/fast: edit_cache
.PHONY : edit_cache/fast

# Special rule for the target rebuild_cache
rebuild_cache:
        @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan 
"Running CMake to regenerate build system..."
        /usr/bin/cmake -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
.PHONY : rebuild_cache

# Special rule for the target rebuild_cache
rebuild_cache/fast: rebuild_cache
.PHONY : rebuild_cache/fast

# The main all target
all: cmake_check_build_system
        $(CMAKE_COMMAND) -E cmake_progress_start /tmp/as/CMakeFiles 
/tmp/as/CMakeFiles/progress.marks
        $(MAKE) -f CMakeFiles/Makefile2 all
        $(CMAKE_COMMAND) -E cmake_progress_start /tmp/as/CMakeFiles 0
.PHONY : all

# The main clean target
clean:
        $(MAKE) -f CMakeFiles/Makefile2 clean
.PHONY : clean

# The main clean target
clean/fast: clean
.PHONY : clean/fast

# Prepare targets for installation.
preinstall: all
        $(MAKE) -f CMakeFiles/Makefile2 preinstall
.PHONY : preinstall

# Prepare targets for installation.
preinstall/fast:
        $(MAKE) -f CMakeFiles/Makefile2 preinstall
.PHONY : preinstall/fast

# clear depends
depend:
        $(CMAKE_COMMAND) -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) 
--check-build-system CMakeFiles/Makefile.cmake 1
.PHONY : depend

#=============================================================================
# Target rules for targets named main

# Build rule for target.
main: cmake_check_build_system
        $(MAKE) -f CMakeFiles/Makefile2 main
.PHONY : main

# fast build rule for target.
main/fast:
        $(MAKE) -f CMakeFiles/main.dir/build.make CMakeFiles/main.dir/build
.PHONY : main/fast

#=============================================================================
# Target rules for targets named main-esbtl

# Build rule for target.
main-esbtl: cmake_check_build_system
        $(MAKE) -f CMakeFiles/Makefile2 main-esbtl
.PHONY : main-esbtl

# fast build rule for target.
main-esbtl/fast:
        $(MAKE) -f CMakeFiles/main-esbtl.dir/build.make 
CMakeFiles/main-esbtl.dir/build
.PHONY : main-esbtl/fast

# target to build an object file
main-esbtl.o:
        $(MAKE) -f CMakeFiles/main-esbtl.dir/build.make 
CMakeFiles/main-esbtl.dir/main-esbtl.o
.PHONY : main-esbtl.o

# target to preprocess a source file
main-esbtl.i:
        $(MAKE) -f CMakeFiles/main-esbtl.dir/build.make 
CMakeFiles/main-esbtl.dir/main-esbtl.i
.PHONY : main-esbtl.i

# target to generate assembly for a file
main-esbtl.s:
        $(MAKE) -f CMakeFiles/main-esbtl.dir/build.make 
CMakeFiles/main-esbtl.dir/main-esbtl.s
.PHONY : main-esbtl.s

# target to build an object file
main.o:
        $(MAKE) -f CMakeFiles/main.dir/build.make CMakeFiles/main.dir/main.o
.PHONY : main.o

# target to preprocess a source file
main.i:
        $(MAKE) -f CMakeFiles/main.dir/build.make CMakeFiles/main.dir/main.i
.PHONY : main.i

# target to generate assembly for a file
main.s:
        $(MAKE) -f CMakeFiles/main.dir/build.make CMakeFiles/main.dir/main.s
.PHONY : main.s

# Help Target
help:
        @echo "The following are some of the valid targets for this Makefile:"
        @echo "... all (the default if no target is provided)"
        @echo "... clean"
        @echo "... depend"
        @echo "... edit_cache"
        @echo "... main"
        @echo "... main-esbtl"
        @echo "... rebuild_cache"
        @echo "... main-esbtl.o"
        @echo "... main-esbtl.i"
        @echo "... main-esbtl.s"
        @echo "... main.o"
        @echo "... main.i"
        @echo "... main.s"
.PHONY : help



#=============================================================================
# Special targets to cleanup operation of make.

# Special rule to run CMake to check the build system integrity.
# No rule that depends on this can have commands that come from listfiles
# because they might be regenerated.
cmake_check_build_system:
        $(CMAKE_COMMAND) -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) 
--check-build-system CMakeFiles/Makefile.cmake 0
.PHONY : cmake_check_build_system