diff --git a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h index 28123c632be..8037b2f38fc 100644 --- a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h +++ b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h @@ -59,103 +59,103 @@ namespace internal { template class Interior_polyhedron_3 { - // 3D - typedef typename Kernel::FT FT; - typedef typename Kernel::Plane_3 Plane; - typedef typename Kernel::Point_3 Point; + // 3D + typedef typename Kernel::FT FT; + typedef typename Kernel::Plane_3 Plane; + typedef typename Kernel::Point_3 Point; - // program and solution types - typedef CGAL::Quadratic_program LP; - typedef typename CGAL::Quadratic_program_solution Solution; - typedef typename Solution::Variable_value_iterator Variable_value_iterator; - typedef CGAL::Real_embeddable_traits RE_traits; - typename RE_traits::To_double to_double; - Solution m_solution; - Point m_inside_point; - Point m_optimal_point; + // program and solution types + typedef CGAL::Quadratic_program LP; + typedef typename CGAL::Quadratic_program_solution Solution; + typedef typename Solution::Variable_value_iterator Variable_value_iterator; + typedef CGAL::Real_embeddable_traits RE_traits; + typename RE_traits::To_double to_double; + Solution m_solution; + Point m_inside_point; + Point m_optimal_point; - public: - Point& inside_point() { return m_inside_point; } - const Point& inside_point() const { return m_inside_point; } - Point& optimal_point() { return m_optimal_point; } - const Point& optimal_point() const { return m_optimal_point; } +public: + Point& inside_point() { return m_inside_point; } + const Point& inside_point() const { return m_inside_point; } + Point& optimal_point() { return m_optimal_point; } + const Point& optimal_point() const { return m_optimal_point; } - // Determines if a value is infinite or not - template - inline bool isinf(T value) { - return value == std::numeric_limits::infinity(); - } + // Determines if a value is infinite or not + template + inline bool isinf(T value) { + return value == std::numeric_limits::infinity(); + } - // Find a point inside the polyhedron defined by a list of planes - // InputIterator::value_type = Plane - template < class InputIterator > - bool find(InputIterator begin, InputIterator end) { - // solve linear program - LP lp(CGAL::LARGER,false); // with constraints Ax >= b + // Find a point inside the polyhedron defined by a list of planes + // InputIterator::value_type = Plane + template < class InputIterator > + bool find(InputIterator begin, InputIterator end) { + // solve linear program + LP lp(CGAL::LARGER,false); // with constraints Ax >= b - // column indices - const int index_x1 = 0; - const int index_x2 = 1; - const int index_x3 = 2; - const int index_x4 = 3; + // column indices + const int index_x1 = 0; + const int index_x2 = 1; + const int index_x3 = 2; + const int index_x4 = 3; - // assemble linear program - int j = 0; // row index - // iterate over segments - InputIterator it; - for(it = begin; it != end; ++it, j++) { - const Plane& plane = *it; - const double aj = CGAL::to_double(plane.a()); - const double bj = CGAL::to_double(plane.b()); - const double cj = CGAL::to_double(plane.c()); - const double dj = CGAL::to_double(plane.d()); + // assemble linear program + int j = 0; // row index + // iterate over segments + InputIterator it; + for(it = begin; it != end; ++it, j++) { + const Plane& plane = *it; + const double aj = CGAL::to_double(plane.a()); + const double bj = CGAL::to_double(plane.b()); + const double cj = CGAL::to_double(plane.c()); + const double dj = CGAL::to_double(plane.d()); - CGAL_assertion(!isinf(aj)); - CGAL_assertion(!isinf(bj)); - CGAL_assertion(!isinf(cj)); - CGAL_assertion(!isinf(dj)); + CGAL_assertion(!isinf(aj)); + CGAL_assertion(!isinf(bj)); + CGAL_assertion(!isinf(cj)); + CGAL_assertion(!isinf(dj)); - // plane defined the halfspace: aj * x1 + bj * x2 + cj * x3 + dj <= 0 - // <=> - (aj * x1 + bj * x2 + cj * x3 + dj) >= 0 - // j^th constraint: -(aj * x1 + bj * x2 + cj * x3 + x4) >= dj - lp.set_a(index_x1, j, -aj); - lp.set_a(index_x2, j, -bj); - lp.set_a(index_x3, j, -cj); - lp.set_a(index_x4, j, -1.0); + // plane defined the halfspace: aj * x1 + bj * x2 + cj * x3 + dj <= 0 + // <=> - (aj * x1 + bj * x2 + cj * x3 + dj) >= 0 + // j^th constraint: -(aj * x1 + bj * x2 + cj * x3 + x4) >= dj + lp.set_a(index_x1, j, -aj); + lp.set_a(index_x2, j, -bj); + lp.set_a(index_x3, j, -cj); + lp.set_a(index_x4, j, -1.0); - // right hand side - lp.set_b(j, dj); - } + // right hand side + lp.set_b(j, dj); + } - // objective function -> max x4 (negative sign set because - // the lp solver always minimizes an objective function) - lp.set_c(index_x4,-1.0); + // objective function -> max x4 (negative sign set because + // the lp solver always minimizes an objective function) + lp.set_c(index_x4,-1.0); - // solve the linear program - m_solution = CGAL::solve_linear_program(lp, ET()); + // solve the linear program + m_solution = CGAL::solve_linear_program(lp, ET()); - if(m_solution.is_infeasible()) - return false; + if(m_solution.is_infeasible()) + return false; - if(!m_solution.is_optimal()) - return false; + if(!m_solution.is_optimal()) + return false; - // get variables - Variable_value_iterator X = m_solution.variable_values_begin(); + // get variables + Variable_value_iterator X = m_solution.variable_values_begin(); - // solution if x4 > 0 - double x4 = to_double(X[index_x4]); - if(x4 <= 0.0) - return false; + // solution if x4 > 0 + double x4 = to_double(X[index_x4]); + if(x4 <= 0.0) + return false; - // define inside point as (x1;x2;x3) - double x1 = to_double(X[index_x1]); - double x2 = to_double(X[index_x2]); - double x3 = to_double(X[index_x3]); - m_inside_point = Point(x1,x2,x3); + // define inside point as (x1;x2;x3) + double x1 = to_double(X[index_x1]); + double x2 = to_double(X[index_x2]); + double x3 = to_double(X[index_x3]); + m_inside_point = Point(x1,x2,x3); - return true; - } + return true; + } }; } // end of internal namespace @@ -170,7 +170,6 @@ If the intersection is empty, `boost::none` is returned. is \f$ a\, x +b\, y +c\, z + d = 0 \f$ then the corresponding halfspace is defined by \f$ a\, x +b\, y +c\, z + d \le 0 \f$ . \tparam PlaneIterator must be an input iterator with the value type being a `Plane_3` object from \cgal Kernel - */ template boost::optional::value_type>::Kernel::Point_3>