CGAL users discussion list

[Pierre Alliez <>]

hi all,

inspired by http://www.qhull.org/html/qhalf.htm
(see qhalf notes)

I am using the QP solver (used as LP solver) to compute the kernel of a 
polyon in 2D. as I obtain strange results - I would appreciate if 
someone with experience on the QP solve could take a look at the 
following code. the input is a set of 2D segments.


// LP solver to compute polygon kernels
#include <CGAL/QP_functions.h>
#include <CGAL/QP_models.h>

template <class Kernel, class ET>
class CPolygon_kernel
{
private:

  // basic types
  typedef typename Kernel::FT FT;
  typedef typename Kernel::Line_2 Line;
  typedef typename Kernel::Point_2 Point;
  typedef typename Kernel::Vector_2 Vector;
  typedef typename Kernel::Segment_2 Segment;

  // program and solution types
  typedef CGAL::Quadratic_program<double> LP;
  typedef typename CGAL::Quadratic_program_solution<typename ET> Solution;
  typedef typename Solution::Variable_value_iterator 
Variable_value_iterator;
  Solution m_solution;
  Point m_inside_point;

  typedef CGAL::Real_embeddable_traits<typename 
Variable_value_iterator::value_type> RE_traits;
  typename RE_traits::To_double to_double;

public:
  CPolygon_kernel() {}
  ~CPolygon_kernel() {}

public:
  Point& inside_point() { return m_inside_point; }
  const Point& inside_point() const { return m_inside_point; }
  unsigned int number_of_iterations() { return 
m_solution.number_of_iterations(); }

public:

  // assume polygon given per segment and CCW
  template < class InputIterator >
  bool solve(InputIterator begin, InputIterator end)
  {
    // solve linear program
    LP lp(CGAL::LARGER); // with constraints Ax >= b

    // column indices
    const int index_x1 = 0;
    const int index_x2 = 1;
    const int index_x3 = 2;

    // assemble linear program
    int j = 0; // row index
    InputIterator it;
    for(it = begin; it != end; it++,j++) // iterate over segments
    {
      const Segment& segment = *it;
      const Point a = segment.source();
      const Point b = segment.target();
      Vector ab = b - a;
      Vector normal = unit_normal(ab); // CW 90 deg rotation
      const FT aj = normal.x();
      const FT bj = normal.y();

      // constraint (line j): aj * x1 + bj * x2 - x3 >= -cj
      lp.set_a(index_x1, j,   aj);
      lp.set_a(index_x2, j,   bj);
      lp.set_a(index_x3, j, -1.0);

      // right hand side
      Line line(a,b);
      FT c = to_double(line.c());
      FT cj = distance_to_origin(line) * CGAL::sign(c);
      lp.set_b(j, -cj);

      // or equivalently...
      // Line line(a,normalize(ab));
      // FT cj = to_double(line.c());
      // lp.set_b(j, -cj);
    }

    // objective function -> max x3 (negative sign set because
    // the lp solver always minimizes an objective function)
    lp.set_c(index_x3,-1.0);

    // add constraints over x3 (>= 0). somewhat redundant
    // constraint as we later test for x3 > 0
    lp.set_l(index_x3,true,0.0);

    // solve the linear program
    m_solution = CGAL::solve_linear_program(lp, ET());

    if(m_solution.is_infeasible())
      return false;

    if(!m_solution.is_optimal())
      return false;

    // get variables
    Variable_value_iterator X = m_solution.variable_values_begin();

    // solution if x3 > 0
    double x3 = to_double(X[2]);
    if(x3 <= 0.0)
      return false;

    // define inside point as (x1;x2)
    double x1 = to_double(X[0]);
    double x2 = to_double(X[1]);
    m_inside_point = Point(x1,x2);

    return true;
  }

  Vector normalize(Vector v)
  {
    return v / std::sqrt(v*v);
  }

  Vector unit_normal(const Vector& v)
  {
    Vector normal(v.y(),-v.x()); // 90 deg CW rotation
    return normalize(normal);
  }

  FT distance_to_origin(const Line& line)
  {
    return std::sqrt(CGAL::squared_distance(Point(CGAL::ORIGIN),line));
  }

};