CGAL users discussion list

[Josip Dzolonga <>]

Thank you for your reply.

On Thu, Nov 12, 2009 at 9:46 AM, Manuel Caroli <> wrote:

By the way, I don't understand why you need the set seen but anyway this should not cause the problem.

Because  if I remove the same vertex from multiple facets it belongs to I would get a segfault. This way, I am ensuring that each vertex is removed only once. I have a set of Facets whose vertices I do not need anymore.
Since this is part of a bigger project, it's not easy to extract a compilable example that fails. Pasting the whole would certainly be confusing. So, I will outline the main idea. We are performing region growing on the Facets on the Delaunay
triangulation structure. The growing procedure is iterated until all points are included, and after each step we remove all points that were reached from the seed (i.e. facet vertices). These are hold into grown->facets (has type std::list<Facet>).
The toRemove list holds all vertices that belong to some Facet in the grown set, however every one is included once by the reasons explained above. I will include the code here for completeness.

std::set<Vertex_handle> toRemove;

std::set<Point> seen;

for(std::list<Facet>::iterator iter = grown->facets->begin();

        iter != grown->facets->end(); ++iter) {

       Cell_handle c = (*iter).first;

   int index = (*iter).second;

   Point p1 = c->vertex((index+1)%4)->point();

   Point p2 = c->vertex((index+2)%4)->point();

            Point p3 = c->vertex((index+3)%4)->point();

   // Do not remove points already removed from a neighboring cell

   if(seen.insert(p1).second) {

   toRemove.insert(c->vertex((index+1)%4));

   }

   if(seen.insert(p2).second) {

   toRemove.insert(c->vertex((index+2)%4));

   }

   if(seen.insert(p3).second) {

   toRemove.insert(c->vertex((index+3)%4));

   }

}

T.remove(toRemove.begin(), toRemove.end());

At the moment I can't provide you with an example that crashes since this is part of a bigger applications and it would be a lot of work to get it

What exactly does it mean that the call to circumcenter "fails"? The only issue I see right now would be that the triangle is very flat. Then you might get a division by 0 or an overflow because the circumcenter is very far away...

That is what I'm getting, the triangles are "very flat". The triangle in question is

x: 0.74084 -4.76403 -0.5324
y: 0 0 0
z: 8.44014e-317 8.25717e-317 8.25554e-317

I guess there are some rounding errors or faults in the code since z is equal to y. Anything else that can cause this? The points are extracted from a Facet (namely T.triangle(facet)[0/1/2]).

Thank you,
Josip
 

Could you check how the triangle looks like, for which circumcenter fails ?

best

Manuel

The call to the circumcenter functions fails, and I suppose that this is because I'm using an inexact_constructions kernel. Is there any trick around this, i.e. some condition/heuristic that will skip these triangles? Changing the kernel itself is of course a possibility, however it will introduce a lot of changes in the whole project (because it is assumed in most of the code that some CGAL functions return something castable to double).

Thank you very much,
Josip

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss