CGAL users discussion list

["Tardugno, Angelo" <>]

Hi CGAL experts,

 

I've been trying to compute the intersection between a Delaunay facet and ist dual (Voronoi Edge). This operation should result in a point, but, surprisingly sometimes (actually very frequently) it returns something different.

Here's part of my code:

 

Finite_cells_iterator ci=T.finite_cells_begin();
  out_stream.open("Vvertex");
  o_stream.open("s3points");
  o1stream.open("triangle");
  o2stream.open("segment");

 

  for (;ci !=T.finite_cells_end(); ++ci)
  {
   voronoi_vertices = (T.dual(ci));
   out_stream << 4 << " " << voronoi_vertices << "\n";

 

   // Finding the dual of each facet of the delaunay triangulation (Voronoi edge)
   // and compute its intersection with the facet in case the triangle
   // is acute

 

 

 

   for (int i=0; i<4; i++)
   {
    CGAL::Object o = T.dual(ci,i);
    K::Triangle_3 t = T.triangle(ci,i);
    CGAL::Object o2;
       K::Segment_3 s,s2;
       K::Ray_3     r;
    K::Plane_3 plane;
    K::Point_3 point;

 

    K::Construct_supporting_plane_3 construct_plane;

 

    

 

    if (is_acute(t))
    {
        if (CGAL::assign(s,o))
        {
      //cout << i << "\n";
      o1stream << t << "\n";
      o2stream << s << "\n";
      plane = construct_plane(t);        // Building the supporting plane for the triangle
      o2 = intersection(s,plane);        // Evaluating intersection between voronoi edge and plane
      if (assign(point,o2))         // Assigning object intersection to point
      {
       xx = (point.x());
       yy = (point.y());
       zz = (point.z());
       out_stream << 3 << " " << xx << " " << yy  << " " << zz << "\n";

 

       n_of_s3_points++;
      } else
      {
      cout << "Something different\n";
      }
     }
        if (CGAL::assign(r,o))
        {
      //cout << i << "\n";
      plane = construct_plane(t);        // Building the supporting plane for the triangle
      o2 = intersection(s,plane);        // Evaluating intersection between voronoi edge and plane
      if (assign(point,o2))         // Assigning object intersection to point
      {
       xx = (point.x());
       yy = (point.y());
       zz = (point.z());
       o_stream <<  3 << " " << xx << " " << yy  << " " << zz << "\n";

 

       n_of_s3_points++;
      } else
      {
            cout << "Something different\n";
      }

 

     }
    }
    else cout << "\nObtuse Angle\n";

 

   }

 

  }

 

  o1stream.close();
  o2stream.close();
  o_stream.close();
  out_stream.close();
 

After running this program, I checked the files triangle and segment but when I check for the first triangle and the first segment it seems like the segment isn't the actual dual of the Delaunay facet, since there is no intersection between the two. Maybe I've made a mistake when I created the dual of the facet or the triangle describing the facet. Anyway, I've tried my best to understand what's happening but I have no clue.

 

Anybody can  help?

 

Thanks

 

Angelo Tardugno