CGAL users discussion list

[Quentin Merigot <>]

Hi Evan,

I attach a function that does exactly what Olivier suggest in order to compute the intersection of halfspaces in 3D, assuming that the origin is contained inside the intersection. Note that the code makes some geometric constructions and is not robust to degenerate situations. (The right way to do this would probably be to write modified predicates for convex_hull_3 in order to work to work implicitely with the points that are the dual of the planes.)

PS: despite the CGAL namespace, this code is not part of CGAL, all bugs are mine.

Best,

Quentin

2013/4/30 Olivier Devillers <>

Le 4/29/13 5:58 PM, Evan Behar a écrit :

Hi,

I have computed many Plane_3 from points and vectors that I am given, and I want to combine them into a convex polyhedron. I have been trying to do this by constructing Nef_polyhedron_3 from the Plane_3 to create the halfspace and then iteratively intersecting the Nef_polyhedron_3, however this becomes quite slow after a while.

It is not important that I specifically have a Nef_polyhedron_3, a Polyhedron_3 will work just fine. Given that, is there a canonical, efficient way to compute a convex polyhedron from the intersection of half-spaces?

If you know a point inside the intersection then, by duality, an intersection of half spaces corresponds to a 3D convex hull of points.

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

#ifndef CGAL_HALFSPACES_INTERSECTION_HPP
#define CGAL_HALFSPACES_INTERSECTION_HPP

#include <CGAL/Convex_hull_traits_3.h>
#include <CGAL/convex_hull_3.h>
#include <CGAL/Polyhedron_incremental_builder_3.h>

CGAL_BEGIN_NAMESPACE
namespace internal
{
  template <class Polyhedron>
  class Build_dual_polyhedron :
    public CGAL::Modifier_base<typename Polyhedron::HalfedgeDS>
  { 
    typedef typename Polyhedron::HalfedgeDS HDS;
    const Polyhedron &_primal;
      
  public:
    Build_dual_polyhedron (const Polyhedron & primal):
      _primal (primal)
    {}
      
    void operator () (HDS &hds)
    {
      typedef typename Polyhedron::Facet Facet;
      typedef typename Polyhedron::Facet_const_handle Facet_const_handle;
      typedef typename Polyhedron::Facet_const_iterator Facet_const_iterator;
      typedef typename Polyhedron::Vertex_const_iterator Vertex_const_iterator;
      typename CGAL::Polyhedron_incremental_builder_3<HDS> B( hds, true);
      
      B.begin_surface(_primal.size_of_facets(),
		      _primal.size_of_vertices(),
		      _primal.size_of_vertices());
      
      // compute coordinates of extreme vertices in the dual polyhedron
      // from primal faces
      std::map<Facet_const_handle, size_t> extreme_points;
      size_t n = 0;

      for (Facet_const_iterator it = _primal.facets_begin();
	   it != _primal.facets_end(); ++it, ++n)
	{
	  typename Facet::Halfedge_const_handle h = it->halfedge();
	  typename Facet::Plane_3 p ( h->vertex()->point(),
				      h->next()->vertex()->point(),
				      h->next()->next()->vertex()->point());
	  B.add_vertex(CGAL::ORIGIN + p.orthogonal_vector () / (-p.d()));
	  extreme_points[it] = n;
	}
      
      // build faces
      for (Vertex_const_iterator it = _primal.vertices_begin();
	   it != _primal.vertices_end(); ++it)
	{
	  assert (it->is_bivalent() == false);

	  typename Polyhedron::Halfedge_around_vertex_const_circulator
	    h0 = it->vertex_begin(), hf = h0;

	  B.begin_facet();
	  do
	    {
	      B.add_vertex_to_facet(extreme_points[hf->facet()]);
	    } while (++hf != h0);
	  B.end_facet();
	}

      B.end_surface();
    }
  };

  template <class PlaneIterator, class Polyhedron, class K>
  void
  halfspaces_intersection(PlaneIterator pbegin, PlaneIterator pend,
			  Polyhedron &P, const K &k)
  {
    typedef typename CGAL::Convex_hull_traits_3<K> Traits;
    typedef typename Polyhedron::HalfedgeDS HalfedgeDS;
    typedef typename K::Point_3 Point;
    typedef typename CGAL::internal::Build_dual_polyhedron<Polyhedron> Builder;

    // construct dual points to apply the convex hull
    std::vector<Point> dual_points;
    for (PlaneIterator p = pbegin; p != pend; ++p)
      dual_points.push_back(CGAL::ORIGIN + p->orthogonal_vector () / (-p->d()));

    Polyhedron ch;
    CGAL::convex_hull_3(dual_points.begin(), dual_points.end(), ch);    

    Builder build_dual (ch);
    P.delegate(build_dual);
  }
}

CGAL_END_NAMESPACE

#endif