CGAL users discussion list

["J.L.M." <>]

John,

I am interested in how you implemented your Minkowski difference using 
the bounding box method. My code throws an assertion in the 
minkowski_sum_3 function on my attached shapes:

terminate called after throwing an instance of 'CGAL::Assertion_exception'
  what():  CGAL ERROR: assertion violation!
Expr: p != 0
File: 
/import/lightcycle/home/jlm/school/masters_thesis/code/CGAL-3.7/include/CGAL/Union_find.h
Line: 135
Aborted

Could you run my shapes through your Minkowski difference code? For
A ominus B, set A=l_shape.off and B=cube.off.

How are you creating your bounding box? I am getting the minimum and 
maximum X,Y,Z values for both shapes, adding them together, and adding a 
little bit more for padding. These then give me two coordinates at 
opposite corners of a cube. I enter these values into a function that 
generates a Polyhedron_3 as I pulled out of the CGAL manual, and then 
convert that into a Nef_polyhedron_3. Beyond that, my code is pretty 
much the same as what you describe. In the attached C++ code, I've 
removed my invert function, but that doesn't eliminate the assertion.

If you happen to have some time and can review why my code causes 
assertions and yours doesn't I would greatly appreciate it. It should be 
self-explanatory with comments, but if it isn't don't waist too much 
time figuring it out.

Thanks.

On 12/28/2010 09:31 PM, johnzjq wrote:
> Hi,
>
> Many thanks for your detailed explaination. And your suggestions are really
> helpful.
>
> I have implmented the code I attached in the post "Any idea about Minkowski
> difference?" and it works fine on many complex CAD models (which should be
> 2-manifold) , but a bit slow.
>
> Regards,
>
> John
OFF
8 6 0
1.000000 -1.000000 -1.000000
1.000000 -1.000000 1.000000
-1.000000 -1.000000 1.000000
-1.000000 -1.000000 -1.000000
1.000000 1.000000 -1.000000
1.000000 1.000000 1.000000
-1.000000 1.000000 1.000000
-1.000000 1.000000 -1.000000
4 0 1 2 3
4 4 7 6 5
4 0 4 5 1
4 1 5 6 2
4 2 6 7 3
4 4 0 3 7
OFF
16 14 0
1.000000 -1.000000 -1.000000
1.000000 -1.000000 1.000000
-1.000000 -1.000000 1.000000
-1.000000 -1.000000 -1.000000
3.000000 3.000000 -1.000000
3.000000 3.000000 1.000000
-1.000000 3.000000 1.000000
-1.000000 3.000000 -1.000000
1.000000 1.000000 1.000000
3.000000 1.000000 1.000000
1.000000 1.000000 -1.000000
3.000000 1.000000 -1.000000
1.000000 3.000000 1.000000
-1.000000 1.000000 1.000000
-1.000000 1.000000 -1.000000
1.000000 3.000000 -1.000000
4 0 1 2 3
4 4 5 9 11
4 13 2 1 8
4 13 14 3 2
4 6 7 14 13
4 6 13 8 12
4 12 8 9 5
4 15 12 5 4
4 7 6 12 15
4 0 10 8 1
4 10 11 9 8
4 0 3 14 10
4 7 15 10 14
4 15 4 11 10
/*Copyright (C) 2010 Jason M'Sadoques


This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

Polyhedron_3 bounding_box(Point_3, Point_3 );

Nef_polyhedron_3 one_nn_boundingbox(const Nef_polyhedron_3& container,
                                    const Nef_polyhedron_3& object)
{

  //Determine a bounding box for both objects. We do this by iterating over each vertex.
  Kernel::FT min_x, min_y, min_z, max_x, max_y, max_z;
  Nef_polyhedron_3::Vertex_const_iterator vci(container.vertices_begin());
  min_x=max_x=vci->point().x();
  min_y=max_y=vci->point().y();
  min_z=max_z=vci->point().z();
  ++vci;
  for(; vci!=container.vertices_end(); ++vci)
  {
    if(vci->point().x()<min_x) min_x=vci->point().x();
    if(vci->point().y()<min_y) min_y=vci->point().y();
    if(vci->point().z()<min_z) min_z=vci->point().z();
    if(vci->point().x()>max_x) max_x=vci->point().x();
    if(vci->point().y()>max_y) max_y=vci->point().y();
    if(vci->point().z()>max_z) max_z=vci->point().z();
  }
  Point_3 container_min(min_x, min_y, min_z);
  Point_3 container_max(max_x, max_y, max_z);

  Nef_polyhedron_3 inverted_object((object));
  vci=inverted_object.vertices_begin();
  min_x=max_x=vci->point().x();
  min_y=max_y=vci->point().y();
  min_z=max_z=vci->point().z();
  ++vci;
  for(; vci!=inverted_object.vertices_end(); ++vci)
  {
    if(vci->point().x()<min_x) min_x=vci->point().x();
    if(vci->point().y()<min_y) min_y=vci->point().y();
    if(vci->point().z()<min_z) min_z=vci->point().z();
    if(vci->point().x()>max_x) max_x=vci->point().x();
    if(vci->point().y()>max_y) max_y=vci->point().y();
    if(vci->point().z()>max_z) max_z=vci->point().z();
  }
  //Add an additional 1% to the min and max values found to give a small bit of extra space for
  //the bounding box. We're assuming that the origin is an element of the object, so the min values
  //will be negative.
  Vector_3 object_min(min_x+0.01*min_x, min_y+0.01*min_y, min_z+0.01*min_z);
  Vector_3 object_max(max_x+0.01*max_x, max_y+0.01*max_y, max_z+0.01*max_z);

  //Create the bounding box corner points.
  Point_3 bounding_box_min(container_min+object_min);
  Point_3 bounding_box_max(container_max+object_max);

  Polyhedron_3 bbox_poly(bounding_box(bounding_box_min, bounding_box_max));
  Nef_polyhedron_3 bbox(bbox_poly);
  Nef_polyhedron_3 bbox_minus_container(bbox-container);

  Nef_polyhedron_3 result=CGAL::minkowski_sum_3(bbox_minus_container, inverted_object);
  //TODO: Try bbox - result, also try container - result.
  return(bbox * !result);
}

//Taken from the CGAL online manual. Creates a cube with opposite corners min_corner and max_corner.
Polyhedron_3 bounding_box(Point_3 min_corner, Point_3 max_corner)
{
  typedef Polyhedron_3::Halfedge_handle Halfedge_handle;
  Polyhedron_3 result;
  Halfedge_handle h=result.make_tetrahedron(Point_3( max_corner.x(), min_corner.y(), min_corner.z()),
                                            Point_3( min_corner.x(), min_corner.y(), max_corner.z()),
                                            Point_3( min_corner.x(), min_corner.y(), min_corner.z()),
                                            Point_3( min_corner.x(), max_corner.y(), min_corner.z()));
  Halfedge_handle g=h->next()->opposite()->next();
  result.split_edge(h->next());
  result.split_edge(g->next());
  result.split_edge(g);
  h->next()->vertex()->point()=Point_3( max_corner.x(), min_corner.y(), max_corner.z());
  g->next()->vertex()->point()=Point_3( min_corner.x(), max_corner.y(), max_corner.z());
  g->opposite()->vertex()->point()=Point_3( max_corner.x(), max_corner.y(), min_corner.z());
  Halfedge_handle f=result.split_facet( g->next(), g->next()->next()->next());
  Halfedge_handle e=result.split_edge(f);
  e->vertex()->point()=Point_3( max_corner.x(), max_corner.y(), max_corner.z());
  result.split_facet(e, f->next()->next());
  CGAL_postcondition(result.is_valid());
  return result;
}