CGAL users discussion list

["Sebastien Loriot (GeometryFactory)" <>]

Here is one way to compute the volume you are looking for.
Note that it is not robust as the periodic translation of points
as well as circumcenter computation require exact constructions.

It should be working for both periodic and non-periodic triangulations.


S.

Maarten Moesen wrote:
> Dear all,
>
> I'm looking for a quick way to compute or estimate the volume of a 
> Voronoi cell around a vertex v in a periodic triangulation. Can one just 
> compute the volume of the incident tetrahedra and divide this by four? 
> Or does one really need to construct the Voronoi cell (and partition it 
> somehow e.g. into pyramids around v?
>
> Many thanks,
>
> Maarten
> --
> Maarten Moesen, PHD
> Department of Metallurgy and Materials Engineering (MTM)
> K.U.Leuven
> Kasteelpark Arenberg 44 - Bus 02450
> B-3001 Heverlee, Belgium
> tel. +32 (0)16 32 13 17
> fax. +32 (0)16 32 19 90

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_triangulation_traits_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>



template <class Triangulation>
typename Triangulation::Facet get_first_vertex(
    Triangulation& T,
    const typename Triangulation::Facet& f,
    const typename Triangulation::Edge& e)
{
  int r=0;
  int c=f.first->index( e.first->vertex(e.second) );
  int l=f.first->index( e.first->vertex(e.third) );
  while (r==c || r==l || r==f.second)
    ++r;
  //if orientation(l,c,r) seen from f.second
  //is counterclockwise then return f else return f.mirror_facet()
  //if f.second =         0    2     1     3
  //then we              213  103   302   201
  //observe              321  310   023   012
  //if ccw               132  031   230   120
  //MAYBE SHOULD USE FUNCTION ccw in Triangulation_utils_3
  int a[3]={l,c,r};
  int im=0;
  if (f.second%2==0){
    while (a[im]!=3) ++im;
    if(a[(im+2)%3] < a[(im+1)%3])
      return f;
  }
  else{
    while (a[im]!=0) ++im;
    if (a[(im+2)%3] > a[(im+1)%3])
      return f;
  }
  return T.mirror_facet(f);
}


template <class Triangulation>
bool is_infinite_vertex(const Triangulation& T, typename Triangulation::Vertex_handle vh)
{
  return vh==T.infinite_vertex();
}

template <class PT,class TDS,class Vertex_handle>
bool is_infinite_vertex(const CGAL::Periodic_3_Delaunay_triangulation_3<PT,TDS>&,Vertex_handle)
{
  return false;
}


template <class Triangulation>
typename Triangulation::Point
point(const Triangulation&, typename Triangulation::Vertex_handle vh)
{
  return vh->point();
}

template <class PT,class TDS,class Vertex_handle>
typename CGAL::Periodic_3_Delaunay_triangulation_3<PT,TDS>::Points
point(const CGAL::Periodic_3_Delaunay_triangulation_3<PT,TDS>& T,Vertex_handle vh)
{
  return T.point(T.periodic_point(vh));
}





template <class Kernel,class Triangulation,class Circumcenter>
typename Kernel::FT 
get_area_of_dual_delaunay_edge( const Triangulation& T,
                                const typename Triangulation::Triangulation_data_structure::Edge& edge,
                                const typename Kernel::Point_3& p1,
                                const typename Kernel::Point_3& p2,
                                const Circumcenter& circumcenter)
{
  typename Kernel::FT area(0);
   
  //vertices of the edge
  typename Triangulation::Vertex_handle vh1=edge.first->vertex( edge.second );
  typename Triangulation::Vertex_handle vh2=edge.first->vertex( edge.third  );
  typename Triangulation::Facet_circulator f_current=T.incident_facets(edge);
  typename Triangulation::Facet_circulator f_start=f_current;
  
  std::vector<typename Kernel::Point_3> boundary;
  
  do{
    typename Triangulation::Facet oriented_facet=get_first_vertex(T,*f_current,edge);
    typename Triangulation::Vertex_handle vh3=oriented_facet.first->vertex( oriented_facet.second );
    //test if the dual facet is infinite.
    if (is_infinite_vertex(T,vh3))
      return 0;
    typename Triangulation::Vertex_handle vh4;
    for (int i=1;i<4;++i){
      vh4=oriented_facet.first->vertex( (oriented_facet.second+i)%4 );
      if (vh4!=vh1 && vh4!=vh2) break;
    }

    //test if the dual facet is infinite.
    if (is_infinite_vertex(T,vh4))
      return 0;
    
    assert (vh1!=vh2 && vh1!=vh3 && vh1!=vh4 && vh2!=vh3 && vh2!=vh4 && vh3!=vh4 );
    
    typename Kernel::Point_3 p3=point(T,vh3);
    typename Kernel::Point_3 p4=point(T,vh4);
    
    boundary.push_back( circumcenter( p1,p2,p3,p4 ) );
    --f_current;
  }while(f_start!=f_current);
  
  int n=boundary.size();
  //compute area using a triangulation of convex polygon
  for (int i=1; i<n-1;++i ){
    typename Kernel::Vector_3 vect1=boundary[i]-boundary[0];
    typename Kernel::Vector_3 vect2=boundary[i+1]-boundary[0];
    area+=sqrt(CGAL::cross_product(vect1,vect2).squared_length()) / typename Kernel::FT(2);
  }  
  return area;
}


template <class Triangulation>
double
volume_of_finite_dual_cell(const Triangulation& T,typename Triangulation::Vertex_handle v){
  typedef typename Triangulation::Triangulation_data_structure TDS;
  typedef typename Triangulation::Geom_traits                  Kernel;
  
  typename Kernel::Construct_circumcenter_3 circumcenter=Kernel().construct_circumcenter_3_object();
  
  double volume_t=0;
  typename std::list<typename Triangulation::Edge> edgelist;
  T.incident_edges(v,std::back_inserter(edgelist));
  
  for (typename std::list<typename Triangulation::Edge>::iterator 
          e_current=edgelist.begin();
          e_current!=edgelist.end();
          ++e_current)
  {
    //vertices of the edge
    typename Triangulation::Vertex_handle vh1=e_current->first->vertex( e_current->second );
    typename Triangulation::Vertex_handle vh2=e_current->first->vertex( e_current->third  );    
    //if the vertex is infinite or on the convex hull
    if (is_infinite_vertex(T,vh1) || is_infinite_vertex(T,vh2) )
      return 0;
    typename Kernel::Point_3  p1=point(T,vh1);
    typename Kernel::Point_3  p2=point(T,vh2);
    double area_t=get_area_of_dual_delaunay_edge<Kernel>(T,*e_current,p1,p2,circumcenter);
    double height=CGAL::sqrt( squared_distance( p1,p2 ) ) / double(2);    
    volume_t+=area_t*height/double(3);
  }
  return volume_t;
}


typedef CGAL::Exact_predicates_inexact_constructions_kernel     Kernel;
typedef CGAL::Delaunay_triangulation_3<Kernel>                  Triangulation;
typedef CGAL::Periodic_3_triangulation_traits_3<Kernel>         PK;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<PK>           P3DT3;

int main(){
  Triangulation T;
  volume_of_finite_dual_cell(T,T.finite_vertices_begin());
  
  P3DT3 pdt3;
  volume_of_finite_dual_cell(pdt3,pdt3.finite_vertices_begin());
  
}