CGAL users discussion list

[Manuel Caroli <>]

Hi Issha,

Issha Kayo wrote:
> [...]
>
> Assume a point at (0.01 0.5 0.5) might be in a cell, one of whose
> vertices is out of the original cube. In this case I expect
> point(Periodic_point) to return the position of the vertex, say, (-0.05,
> 0.5, 0.5). The actual return, however, seems to be (0.95, 0.5, 0.5).
As you have noted, vertices always store their representative in the
original cube. In the case a cell "wraps around" its representation is
transposed in positive x-, y-, z-direction, respectively. That is, in
your case we do not transpose
        (0.95, 0.5, 0.5) to (-0.05, 0.5, 0.5)     (*)
but rather
        (0.01, 0.5, 0.5) to (1.01, 0.5,0.5)
That choice is somewhat arbitrary, but it was necessary to define a
consistent behavior.

If you really need the output to behave like in (*) you could iterate
through all four periodic_points of a cell and subtract 1 from the
offset (second entry) of each vertex in each direction the cell wraps
around (example below). Beware that this works only if you are computing
in 1-sheeted covering space.

hope this helps

Manuel


============================================================

Example:

// Additionally required typedefs
typedef PDT::Periodic_point         Periodic_point;
typedef PDT::Offset                 Offset;

// Extract the Periodic_points from the data structure
Periodic_point pp[4];
for (int i=0 ; i<4 ; i++) {
 pp[i] = T.periodic_point(c,i);
}

// Figure out in which directions the cell wraps around
int off[3];
for (int d=0 ; d<3 ; d++) {
 bool does_wrap=false;
 for (int i=0 ; i<4 ; i++) {
   does_wrap = does_wrap || (pp[i].second[d]!=0);
 }
 off[d] = (does_wrap ? -1 : 0);
}

// wrap_offset contains the required translation
Offset wrap_offset(off[0],off[1],off[2]);

// Apply the translation in wrap_offset to the Periodic_points
for (int i=0 ; i<4 ; i++) {
 pp[i] = std::make_pair(pp[i].first, pp[i].second + wrap_offset);
 std::cout<< T.point(pp[i])<<std::endl;
}

> As an example, please see the following code. I can easily imagine I'm
> doing something very silly ....
>
> Thank you,
> Issha
>
> -------------
> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
> #include <CGAL/Periodic_3_triangulation_traits_3.h>
> #include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
>
> #include <iostream>
>
> typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
> typedef CGAL::Periodic_3_triangulation_traits_3<K> GT;
> typedef CGAL::Periodic_3_Delaunay_triangulation_3<GT> PDT;
> typedef PDT::Point                  Point;
>
> int main()
> {
>   PDT T;
>
>   // Input point grid (27 points)
>   for (double x=0.1 ; x < .9 ; x += 0.33) {
>     for (double y=0.1 ; y < .9 ; y += 0.33) {
>       for (double z=0.1 ; z < .9 ; z += 0.33) {
>         T.insert(Point(x,y,z)); 
>   } } }
>
>   Point p(0.05, 0.5, 0.5);
>   PDT::Cell_handle c;
>   c=T.locate(p);
>
>   T.convert_to_1_sheeted_covering();
>
>   std::cout << T.point(T.periodic_point(c,0)) << std::endl;
>   std::cout << T.point(T.periodic_point(c,1)) << std::endl;
>   std::cout << T.point(T.periodic_point(c,2)) << std::endl;
>   std::cout << T.point(T.periodic_point(c,3)) << std::endl;
>
>   return 0;
> }
> -----------------------------
>
> The resulting output is
>
> 1.1 0.43 0.43   <--- If 0.1 0.43 0.43, I'm happy.
> 1.1 0.76 0.43   <--- If 0.1 0.76 0.43, I'm happ..
> 0.76 0.76 0.76  <--- If -0.24 0.76 0.76, I'm hap...
> 1.1 0.43 0.76   <--- If 0.1 0.43 0.76, ..