diff --git a/AABB_tree/include/CGAL/AABB_traits.h b/AABB_tree/include/CGAL/AABB_traits.h
index 470420162a7..613d08c174b 100644
--- a/AABB_tree/include/CGAL/AABB_traits.h
+++ b/AABB_tree/include/CGAL/AABB_traits.h
@@ -26,7 +26,6 @@
 #include <CGAL/internal/AABB_tree/Primitive_helper.h>
 
 #include <boost/optional.hpp>
-#include <boost/bind.hpp>
 
 /// \file AABB_traits.h
 
@@ -274,13 +273,13 @@ public:
         switch(Traits::longest_axis(bbox))
         {
         case AT::CGAL_AXIS_X: // sort along x
-          std::nth_element(first, middle, beyond, boost::bind(Traits::less_x,_1,_2,m_traits));
+          std::nth_element(first, middle, beyond, [this](const Primitive& p1, const Primitive& p2){ return Traits::less_x(p1, p2, this->m_traits); });
           break;
         case AT::CGAL_AXIS_Y: // sort along y
-          std::nth_element(first, middle, beyond, boost::bind(Traits::less_y,_1,_2,m_traits));
+          std::nth_element(first, middle, beyond, [this](const Primitive& p1, const Primitive& p2){ return Traits::less_y(p1, p2, this->m_traits); });
           break;
         case AT::CGAL_AXIS_Z: // sort along z
-          std::nth_element(first, middle, beyond, boost::bind(Traits::less_z,_1,_2,m_traits));
+          std::nth_element(first, middle, beyond, [this](const Primitive& p1, const Primitive& p2){ return Traits::less_z(p1, p2, this->m_traits); });
           break;
         default:
           CGAL_error();
diff --git a/Arrangement_on_surface_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h b/Arrangement_on_surface_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h
index 984ea6170a6..47afb1f2f68 100644
--- a/Arrangement_on_surface_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h
+++ b/Arrangement_on_surface_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h
@@ -295,13 +295,6 @@ public:
         return y;
     }
 
-    //! \brief the same as \c evaluate but arguments are passed by value
-    //! (needed to substitute variables in bivariate polynomial)
-    inline static NT binded_eval(Poly_1 poly, NT x)
-    {
-        return evaluate(poly, x);
-    }
-
     //! \brief evalutates a polynomial at certain x-coordinate
     static NT evaluate(const Poly_1& poly, const NT& x,
         bool *error_bounds_ = nullptr)
@@ -913,10 +906,9 @@ void get_precached_poly(int var, const NT& key, int /* level */, Poly_1& poly)
 //     }
 
     if(not_cached||not_found) {
-        poly = Poly_1(::boost::make_transform_iterator(coeffs->begin(),
-                            boost::bind2nd(std::ptr_fun(binded_eval), key1)),
-                      ::boost::make_transform_iterator(coeffs->end(),
-                            boost::bind2nd(std::ptr_fun(binded_eval), key1)));
+        auto fn = [&key1](const Poly_1& poly){ return evaluate(poly, key1); };
+        poly = Poly_1(::boost::make_transform_iterator(coeffs->begin(), fn),
+                      ::boost::make_transform_iterator(coeffs->end(), fn));
         if(not_cached)
             return;
     // all available space consumed: drop the least recently used entry
diff --git a/Arrangement_on_surface_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_traits.h b/Arrangement_on_surface_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_traits.h
index 752389a0f35..56595cd3835 100644
--- a/Arrangement_on_surface_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_traits.h
+++ b/Arrangement_on_surface_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_traits.h
@@ -16,7 +16,6 @@
 
 #include <CGAL/basic.h>
 #include <CGAL/function_objects.h>
-#include <boost/functional.hpp>
 
 /*! \file CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_traits.h
  * \brief
@@ -107,11 +106,13 @@ struct Transform {
 
     template <class X>
     OutputPoly_2 operator()(const CGAL::Polynomial<X>& p, Op op = Op()) const {
-
-        Transform<typename OutputPoly_2::NT, typename InputPoly_2::NT, Op> tr;
+        typedef typename InputPoly_2::NT NT_in;
+        typedef typename OutputPoly_2::NT NT_out;
+        Transform<NT_out, NT_in, Op> tr;
+        auto fn = [&op, &tr](const NT_in& v){ return tr(v, op); };
         return OutputPoly_2(
-            ::boost::make_transform_iterator(p.begin(), boost::bind2nd(tr, op)),
-            ::boost::make_transform_iterator(p.end(), boost::bind2nd(tr, op)));
+            ::boost::make_transform_iterator(p.begin(), fn),
+            ::boost::make_transform_iterator(p.end(), fn));
     }
 
     OutputPoly_2 operator()(
diff --git a/Bounding_volumes/include/CGAL/Min_annulus_d.h b/Bounding_volumes/include/CGAL/Min_annulus_d.h
index 70f45488c0b..0be72c1b427 100644
--- a/Bounding_volumes/include/CGAL/Min_annulus_d.h
+++ b/Bounding_volumes/include/CGAL/Min_annulus_d.h
@@ -32,9 +32,10 @@
 #include <CGAL/QP_solver/QP_full_exact_pricing.h>
 #include <CGAL/boost/iterator/counting_iterator.hpp>
 #include <CGAL/boost/iterator/transform_iterator.hpp>
-#include <boost/functional.hpp>
 #include <CGAL/NT_converter.h>
 
+#include <functional>
+
 // here is how it works. We have d+2 variables:
 // R (big radius), r (small radius), c (center). The problem is
 //
@@ -258,9 +259,6 @@ private:
   typedef  QP_access_by_index
   <typename std::vector<Point>::const_iterator, int> Point_by_index;
 
-  typedef  boost::binder2nd< std::divides<int> >
-  Divide;
-
   typedef  std::vector<int>           Index_vector;
 
   typedef  std::vector<NT>            NT_vector;
@@ -272,7 +270,7 @@ public:
 
   typedef  CGAL::Join_input_iterator_1<
     Basic_variable_index_iterator,
-    CGAL::Unary_compose_1<Point_by_index,Divide> >
+    std::function<Point(int)> >
   Support_point_iterator;
 
 
@@ -331,9 +329,7 @@ public:
                                     "support_points_begin: not enough points");
     return Support_point_iterator(
                                   solver->basic_original_variable_indices_begin(),
-                                  CGAL::compose1_1(
-                                                   Point_by_index( points.begin()),
-                                                   boost::bind2nd( std::divides<int>(), 2)));
+                                  [this](int i){ return Point_by_index(this->points.begin())(i/2); });
   }
 
   Support_point_iterator
@@ -342,9 +338,7 @@ public:
                                     "support_points_begin: not enough points");
     return Support_point_iterator(
                                   solver->basic_original_variable_indices_end(),
-                                  CGAL::compose1_1(
-                                                   Point_by_index( points.begin()),
-                                                   boost::bind2nd( std::divides<int>(), 2)));
+                                  [this](int i){ return Point_by_index(this->points.begin())(i/2); });
   }
 
   int  number_of_inner_support_points() const { return static_cast<int>(inner_indices.size());}
@@ -592,9 +586,7 @@ private:
   bool
   check_dimension( std::size_t  offset = 0)
   { return ( std::find_if( points.begin()+offset, points.end(),
-                           CGAL::compose1_1( boost::bind2nd(
-                                                          std::not_equal_to<int>(), d),
-                                             tco.access_dimension_d_object()))
+                          [this](const Point& p){ return this->d != this->tco.access_dimension_d_object()(p); })
              == points.end()); }
 
   // compute smallest enclosing annulus
diff --git a/Bounding_volumes/include/CGAL/Rectangular_p_center_traits_2.h b/Bounding_volumes/include/CGAL/Rectangular_p_center_traits_2.h
index 4233123e225..bc93c75bf74 100644
--- a/Bounding_volumes/include/CGAL/Rectangular_p_center_traits_2.h
+++ b/Bounding_volumes/include/CGAL/Rectangular_p_center_traits_2.h
@@ -290,8 +290,8 @@ bounding_box_2(ForwardIterator f, ForwardIterator l, const Traits& t)
   return rect(v(rect(*xmin, *ymin), 0), v(rect(*xmax, *ymax), 2));
 } // bounding_box_2(f, l, t)
 template < class ForwardIterator >
-inline typename
-std::iterator_traits< ForwardIterator >::value_type::R::Iso_rectangle_2
+inline
+auto
 bounding_box_2(ForwardIterator f, ForwardIterator l)
 // PRE: f != l.
 {
diff --git a/Bounding_volumes/include/CGAL/min_quadrilateral_2.h b/Bounding_volumes/include/CGAL/min_quadrilateral_2.h
index 867593b2d1d..dfa0eb1e30c 100644
--- a/Bounding_volumes/include/CGAL/min_quadrilateral_2.h
+++ b/Bounding_volumes/include/CGAL/min_quadrilateral_2.h
@@ -20,8 +20,7 @@
 #include <CGAL/basic.h>
 #include <CGAL/Optimisation/assertions.h>
 #include <iterator>
-#include <boost/bind.hpp>
-#include <boost/function.hpp>
+#include <functional>
 
 #ifdef CGAL_OPTIMISATION_EXPENSIVE_PRECONDITION_TAG
 #include <CGAL/Polygon_2_algorithms.h>
@@ -70,13 +69,16 @@ convex_bounding_box_2(
   typedef typename Traits::Point_2                Point_2;
   typedef typename Traits::Less_xy_2              Less_xy_2;
   typedef typename Traits::Less_yx_2              Less_yx_2;
-  typedef boost::function2<bool,Point_2,Point_2>  Greater_xy_2;
-  typedef boost::function2<bool,Point_2,Point_2>  Greater_yx_2;
+  typedef std::function<bool(const Point_2&,
+                             const Point_2&)>     Greater_xy_2;
+  typedef Greater_xy_2                            Greater_yx_2;
 
   Less_xy_2    less_xy_2    = t.less_xy_2_object();
   Less_yx_2    less_yx_2    = t.less_yx_2_object();
-  Greater_xy_2 greater_xy_2 = boost::bind(less_xy_2, _2, _1);
-  Greater_yx_2 greater_yx_2 = boost::bind(less_yx_2, _2, _1);
+  Greater_xy_2 greater_xy_2 = [&less_xy_2](const Point_2& p1, const Point_2& p2)
+                              { return less_xy_2(p2, p1); };
+  Greater_yx_2 greater_yx_2 = [&less_yx_2](const Point_2& p1, const Point_2& p2)
+                              { return less_yx_2(p2, p1); };
 
   if (less_xy_2(*minx, *f) ||
       (less_yx_2(*minx, *f) && !less_xy_2(*f, *minx)))
@@ -268,21 +270,21 @@ namespace Optimisation {
     // ---------------------------------------------------------------
     // Right_of_implicit_line_2
     // ---------------------------------------------------------------
-    typedef boost::function3<bool,Point_2,Point_2,Direction_2>
+    typedef std::function<bool(const Point_2&, const Point_2& , const Direction_2&)>
       Right_of_implicit_line_2;
 
     Right_of_implicit_line_2 right_of_implicit_line_2_object() const {
-      return boost::bind(has_on_negative_side_2_object(),
-                         boost::bind(construct_line_2_object(), _2, _3),
-                         _1);
+      return [this](const Point_2& p1, const Point_2& p2, const Direction_2& d)
+             {  return this->has_on_negative_side_2_object()(this->construct_line_2_object()(p2, d), p1); };
     }
 
-    typedef boost::function2<Direction_2,Point_2,Point_2>
+    typedef std::function<Direction_2(const Point_2&, const Point_2&)>
       Construct_direction_2;
 
     Construct_direction_2 construct_direction_2_object() const {
-      return boost::bind(Base::construct_direction_2_object(),
-                         boost::bind(construct_vector_2_object(), _1, _2));
+      return [this](const Point_2& p1, const Point_2& p2)
+             { return this->Base::construct_direction_2_object()(
+                        this->construct_vector_2_object()(p1, p2)); };
     }
 
     template < class Kernel >
@@ -324,13 +326,11 @@ namespace Optimisation {
     rotate_direction_by_multiple_of_pi_2_object() const
     { return Rotate_direction_by_multiple_of_pi_2(*this); }
 
-    typedef boost::function2<bool,Direction_2,Direction_2>
+    typedef std::function<bool(const Direction_2&, const Direction_2&)>
       Less_angle_with_x_axis_2;
     Less_angle_with_x_axis_2 less_angle_with_x_axis_2_object() const {
-      return boost::bind(std::equal_to<Comparison_result>(),
-                         boost::bind(compare_angle_with_x_axis_2_object(),
-                                     _1, _2),
-                         SMALLER);
+      return [this](const Direction_2& d1, const Direction_2& d2)
+             { return this->compare_angle_with_x_axis_2_object()(d1, d2) == SMALLER; };
     }
 
   };
diff --git a/Bounding_volumes/include/CGAL/pierce_rectangles_2.h b/Bounding_volumes/include/CGAL/pierce_rectangles_2.h
index 9a3bfb91883..ebaf8081e1c 100644
--- a/Bounding_volumes/include/CGAL/pierce_rectangles_2.h
+++ b/Bounding_volumes/include/CGAL/pierce_rectangles_2.h
@@ -22,7 +22,6 @@
 #include <algorithm>
 #include <iterator>
 #include <vector>
-#include <boost/bind.hpp>
 
 #if defined(BOOST_MSVC)
 #  pragma warning(push)
@@ -235,14 +234,16 @@ struct Staircases : public Loc_domain< Traits_ > {
     do {
       brstc.push_back(*i++);
       i = find_if(i, ysort.end(),
-                  boost::bind(this->traits.less_x_2_object(), brstc.back(), _1));
+                  [this](const Point_2& p)
+                  { return this->traits.less_x_2_object()(this->brstc.back(), p);});
     } while (i != ysort.end());
     // top-left
     Riterator j = ysort.rbegin();
     do {
       tlstc.push_back(*j++);
       j = find_if(j, ysort.rend(),
-                  boost::bind(this->traits.less_x_2_object(), _1, tlstc.back()));
+                     [this](const Point_2& p)
+                     { return this->traits.less_x_2_object()(p, this->tlstc.back());});
     } while (j != ysort.rend());
 
     // build left-bottom and right-top staircases
@@ -252,14 +253,16 @@ struct Staircases : public Loc_domain< Traits_ > {
     do {
       lbstc.push_back(*i++);
       i = find_if(i, xsort.end(),
-                  boost::bind(this->traits.less_y_2_object(), _1, lbstc.back()));
+                  [this](const Point_2& p)
+                  { return this->traits.less_y_2_object()(p, this->lbstc.back());});
     } while (i != xsort.end());
     // right-top
     j = xsort.rbegin();
     do {
       rtstc.push_back(*j++);
       j = find_if(j, xsort.rend(),
-                  boost::bind(this->traits.less_y_2_object(), rtstc.back(), _1));
+                  [this](const Point_2& p)
+                  { return this->traits.less_y_2_object()(this->rtstc.back(), p);});
     } while (j != xsort.rend());
   } // Staircases(b, e, t)
 
@@ -300,16 +303,20 @@ struct Staircases : public Loc_domain< Traits_ > {
       min_element_if(
         this->pts.begin(), this->pts.end(),
         this->traits.less_x_2_object(),
-        boost::bind(std::logical_and< bool >(),
-                    boost::bind(this->traits.less_x_2_object(), p, _1),
-                    boost::bind(this->traits.less_y_2_object(), p, _1)));
+        [&p, this](const Point_2& pt)
+        {
+          return this->traits.less_x_2_object()(p, pt) &&
+                 this->traits.less_y_2_object()(p, pt);
+        });
     Citerator j =
       max_element_if(
         this->pts.begin(), this->pts.end(),
         this->traits.less_x_2_object(),
-        boost::bind(std::logical_and< bool >(),
-                    boost::bind(this->traits.less_x_2_object(), _1, q),
-                    boost::bind(this->traits.less_y_2_object(), q, _1)));
+        [&q, this](const Point_2& pt)
+        {
+          return this->traits.less_x_2_object()(pt, q) &&
+                 this->traits.less_y_2_object()(q, pt);
+        });
     return Intervall(i == this->pts.end() ? this->maxx : *i,
                      j == this->pts.end() ? this->minx : *j);
   } // top_intervall()
@@ -326,16 +333,20 @@ struct Staircases : public Loc_domain< Traits_ > {
       min_element_if(
         this->pts.begin(), this->pts.end(),
         this->traits.less_x_2_object(),
-        boost::bind(std::logical_and< bool >(),
-                    boost::bind(this->traits.less_x_2_object(), p, _1),
-                    boost::bind(this->traits.less_y_2_object(), _1, p)));
+        [&p, this](const Point_2& pt)
+        {
+          return this->traits.less_x_2_object()(p, pt) &&
+                 this->traits.less_y_2_object()(pt, p);
+        });
     Citerator j =
       max_element_if(
         this->pts.begin(), this->pts.end(),
         this->traits.less_x_2_object(),
-        boost::bind(std::logical_and< bool >(),
-                    boost::bind(this->traits.less_x_2_object(), _1, q),
-                    boost::bind(this->traits.less_y_2_object(), _1, q)));
+        [&q, this](const Point_2& pt)
+        {
+          return this->traits.less_x_2_object()(pt, q) &&
+                 this->traits.less_y_2_object()(pt, q);
+        });
     return Intervall(i == this->pts.end() ? this->maxx : *i,
                      j == this->pts.end() ? this->minx : *j);
   } // bottom_intervall()
@@ -352,16 +363,20 @@ struct Staircases : public Loc_domain< Traits_ > {
       min_element_if(
         this->pts.begin(), this->pts.end(),
         this->traits.less_y_2_object(),
-        boost::bind(std::logical_and< bool >(),
-                    boost::bind(this->traits.less_x_2_object(), _1, p),
-                    boost::bind(this->traits.less_y_2_object(), p, _1)));
+        [&p, this](const Point_2& pt)
+        {
+          return this->traits.less_x_2_object()(pt, p) &&
+                 this->traits.less_y_2_object()(p, pt);
+        });
     Citerator j =
       max_element_if(
         this->pts.begin(), this->pts.end(),
         this->traits.less_y_2_object(),
-        boost::bind(std::logical_and< bool >(),
-                    boost::bind(this->traits.less_x_2_object(), _1, q),
-                    boost::bind(this->traits.less_y_2_object(), _1, q)));
+        [&q, this](const Point_2& pt)
+        {
+          return this->traits.less_x_2_object()(pt, q) &&
+                 this->traits.less_y_2_object()(pt, q);
+        });
     return Intervall(i == this->pts.end() ? this->maxy : *i,
                      j == this->pts.end() ? this->miny : *j);
   } // left_intervall()
@@ -378,16 +393,20 @@ struct Staircases : public Loc_domain< Traits_ > {
       min_element_if(
         this->pts.begin(), this->pts.end(),
         this->traits.less_y_2_object(),
-        boost::bind(std::logical_and< bool >(),
-                    boost::bind(this->traits.less_x_2_object(), p, _1),
-                    boost::bind(this->traits.less_y_2_object(), p, _1)));
+        [&p, this](const Point_2& pt)
+        {
+          return this->traits.less_x_2_object()(p, pt) &&
+                 this->traits.less_y_2_object()(p, pt);
+        });
     Citerator j =
       max_element_if(
         this->pts.begin(), this->pts.end(),
         this->traits.less_y_2_object(),
-        boost::bind(std::logical_and< bool >(),
-                    boost::bind(this->traits.less_x_2_object(), q, _1),
-                    boost::bind(this->traits.less_y_2_object(), _1, q)));
+        [&q, this](const Point_2& pt)
+        {
+          return this->traits.less_x_2_object()(q, pt) &&
+                 this->traits.less_y_2_object()(pt, q);
+        });
     return Intervall(i == this->pts.end() ? this->maxy : *i,
                      j == this->pts.end() ? this->miny : *j);
   } // right_intervall()
@@ -487,6 +506,7 @@ two_cover_points(
   using std::less;
 
   typedef typename Traits::FT           FT;
+  typedef typename Traits::Point_2           Point_2;
   typename Traits::Infinity_distance_2 dist =
     d.traits.infinity_distance_2_object();
   typename Traits::Signed_infinity_distance_2 sdist =
@@ -504,11 +524,8 @@ two_cover_points(
     if (d.end() ==
         find_if(d.begin(),
                 d.end(),
-                boost::bind(less<FT>(),
-                     d.r,
-                     boost::bind(Min<FT>(),
-                          boost::bind(dist, d[0], _1),
-                          boost::bind(dist, d[2], _1)))))
+                [&dist,&d](const Point_2& p)
+                { return d.r < Min<FT>()(dist(d[0], p), dist(d[2],p)); }))
       {
         *o++ = d[0];
         *o++ = d[2];
@@ -520,11 +537,8 @@ two_cover_points(
     if (d.end() ==
         find_if(d.begin(),
                 d.end(),
-                boost::bind(less<FT>(),
-                     d.r,
-                     boost::bind(Min<FT>(),
-                          boost::bind(dist, d[1], _1),
-                          boost::bind(dist, d[3], _1)))))
+                [&dist,&d](const Point_2& p)
+                { return d.r < Min<FT>()(dist(d[1], p), dist(d[3],p)); }))
       {
         *o++ = d[1];
         *o++ = d[3];
@@ -551,7 +565,6 @@ three_cover_points(
   CGAL_optimisation_precondition(!d.empty());
 
   // typedefs:
-  typedef typename Traits::FT                      FT;
   typedef typename Traits::Point_2                 Point_2;
   typedef typename Loc_domain< Traits >::Iterator  Iterator;
   typename Traits::Infinity_distance_2 dist =
@@ -565,7 +578,8 @@ three_cover_points(
 
     // find first point not covered by the rectangle at d[k]
     Iterator i = find_if(d.begin(), d.end(),
-                         boost::bind(less<FT>(), d.r, boost::bind(dist, corner, _1)));
+                         [&d,&dist, &corner](const Point_2& p)
+                         { return d.r < dist(corner, p); });
 
     // are all points already covered?
     if (i == d.end()) {
@@ -608,12 +622,12 @@ three_cover_points(
 
     CGAL_optimisation_expensive_assertion(
       save_end == find_if(d.end(), save_end,
-                          boost::bind(less<FT>(), d.r, boost::bind(dist, corner, _1))));
+                         [&d, &dist, &corner](const Point_2& p)
+                         { return d.r < dist(corner, p); }));
     CGAL_optimisation_expensive_assertion(
       d.end() == find_if(d.begin(), d.end(),
-                         boost::bind(std::greater_equal<FT>(),
-                              d.r,
-                              boost::bind(dist, corner, _1))));
+                         [&d,&dist, &corner](const Point_2& p)
+                         { return d.r >= dist(corner, p); }));
 
 
     two_cover_points(d, o, ok);
@@ -702,7 +716,8 @@ four_cover_points(Staircases< Traits >& d, OutputIterator o, bool& ok)
 
     // find first point not covered by the rectangle at d[k]
     Iterator i = find_if(d.begin(), d.end(),
-                         boost::bind(less<FT>(), d.r, boost::bind(dist, corner, _1)));
+                         [&d,&dist,&corner](const Point_2& p)
+                         { return d.r < dist(corner, p); });
 
     // are all points already covered?
     if (i == d.end()) {
@@ -745,12 +760,12 @@ four_cover_points(Staircases< Traits >& d, OutputIterator o, bool& ok)
 
     CGAL_optimisation_expensive_assertion(
       save_end == find_if(d.end(), save_end,
-                          boost::bind(less<FT>(), d.r, boost::bind(dist, corner, _1))));
+                          [&d,&dist,&corner](const Point_2& p)
+                          { return d.r < dist(corner, p); }));
     CGAL_optimisation_expensive_assertion(
       d.end() == find_if(d.begin(), d.end(),
-                         boost::bind(std::greater_equal<FT>(),
-                              d.r,
-                              boost::bind(dist, corner, _1))));
+                         [&d,&dist,&corner](const Point_2& p)
+                         { return d.r >= dist(corner, p); }));
 
 
     three_cover_points(d, o, ok);
diff --git a/Bounding_volumes/include/CGAL/rectangular_3_center_2.h b/Bounding_volumes/include/CGAL/rectangular_3_center_2.h
index 531fe486fd4..cd3a1f14bee 100644
--- a/Bounding_volumes/include/CGAL/rectangular_3_center_2.h
+++ b/Bounding_volumes/include/CGAL/rectangular_3_center_2.h
@@ -22,8 +22,6 @@
 #include <CGAL/Rectangular_p_center_traits_2.h>
 #include <algorithm>
 #include <vector>
-#include <boost/bind.hpp>
-#include <boost/function.hpp>
 
 namespace CGAL {
 
@@ -38,8 +36,6 @@ rectangular_2_center_2(
   Traits& t)
 {
   using std::pair;
-  using std::greater;
-  using std::less;
 
   typedef typename Traits::Iso_rectangle_2        Rectangle;
   typedef typename Traits::Point_2                Point;
@@ -53,7 +49,6 @@ rectangular_2_center_2(
     P_below_right;
   typedef typename Traits::Construct_point_2_below_left_implicit_point_2
     P_below_left;
-  typedef boost::function1<FT, Point>              Gamma;
 
   // fetch function objects from traits class
   CVertex       v      = t.construct_vertex_2_object();
@@ -72,15 +67,15 @@ rectangular_2_center_2(
 
   // two cases: top-left & bottom-right or top-right & bottom-left
   Min< FT > minft;
-  Gamma gamma1 =
-    boost::bind(minft, boost::bind(dist, v(bb, 0), _1), boost::bind(dist, v(bb, 2), _1));
-  Gamma gamma2 =
-    boost::bind(minft, boost::bind(dist, v(bb, 1), _1), boost::bind(dist, v(bb, 3), _1));
+  auto gamma1 =
+    [&minft, &bb, &dist, &v](const Point& p){ return minft(dist( v(bb, 0), p), dist(v(bb, 2),p));};
+  auto gamma2 =
+    [&minft, &bb, &dist, &v](const Point& p){ return minft(dist( v(bb, 1), p), dist(v(bb, 3),p));};
 
   pair< ForwardIterator, ForwardIterator > cand =
     min_max_element(f, l,
-                    boost::bind(greater<FT>(), boost::bind(gamma1, _1), boost::bind(gamma1, _2)),
-                    boost::bind(less<FT>(), boost::bind(gamma2, _1), boost::bind(gamma2, _2)));
+                    [&gamma1](const Point& p1, const Point& p2){ return std::greater<FT>()(gamma1(p1), gamma1(p2)); },
+                    [&gamma2](const Point& p1, const Point& p2){ return std::less<FT>()(gamma2(p1), gamma2(p2)); });
 
   // return the result
   if (gamma1(*cand.first) < gamma2(*cand.second)) {
@@ -106,9 +101,6 @@ rectangular_3_center_2_type1(
   typename Traits::FT& rad,
   Traits& t)
 {
-  using std::max;
-  using std::less;
-
   typedef typename Traits::FT                         FT;
   typedef typename Traits::Iso_rectangle_2            Rectangle;
   typedef typename Traits::Point_2                    Point;
@@ -124,7 +116,6 @@ rectangular_3_center_2_type1(
     P_below_right;
   typedef typename Traits::Construct_point_2_below_left_implicit_point_2
     P_below_left;
-  typedef boost::function1<FT, Point>                 Gamma;
 
   // fetch function objects from traits class
   Rect          rect   = t.construct_iso_rectangle_2_object();
@@ -158,14 +149,14 @@ rectangular_3_center_2_type1(
     RandomAccessIterator e = l;
     bool b_empty = true;
     Min< FT > minft;
-    Gamma gamma = boost::bind(minft,
-                       boost::bind(dist, v(r, i), _1),
-                       boost::bind(dist, v(r, 2 + i), _1));
+    auto gamma = [&minft, &dist, &v, &r, i](const Point& p)
+                  { return minft(dist(v(r, i), p), dist(v(r, 2 + i), p)); };
 
     while (e - s > 1) {
       // step (a)
       RandomAccessIterator m = s + (e - s - 1) / 2;
-      std::nth_element(s, m, e, boost::bind(less<FT>(), boost::bind(gamma, _1), boost::bind(gamma, _2)));
+      std::nth_element(s, m, e, [&gamma](const Point& p1, const Point& p2)
+                                {return gamma(p1) < gamma(p2);});
 
       // step (b)
       Rectangle b_prime = bounding_box_2(m + 1, e, t);
@@ -220,8 +211,9 @@ struct Rectangular_3_center_2_type2_operations_base {
   typedef typename R::Infinity_distance_2        Infinity_distance_2;
   typedef typename R::Less_x_2                   Less_x_2;
   typedef typename R::Less_y_2                   Less_y_2;
-  typedef boost::function2<bool,Point_2,Point_2> Greater_x_2;
-  typedef boost::function2<bool,Point_2,Point_2> Greater_y_2;
+  typedef std::function<bool(const Point_2& ,
+                             const Point_2&)>    Greater_x_2;
+  typedef Greater_x_2                            Greater_y_2;
   typedef Min< Point_2, Less_x_2 >               Min_x_2;
   typedef Max< Point_2, Less_x_2 >               Max_x_2;
   typedef Min< Point_2, Less_y_2 >               Min_y_2;
@@ -236,15 +228,15 @@ struct Rectangular_3_center_2_type2_operations_base {
     Construct_point_2_below_right_implicit_point_2;
   typedef typename R::Construct_point_2_below_left_implicit_point_2
     Construct_point_2_below_left_implicit_point_2;
-  typedef boost::function1<FT,Point_2>           Delta;
+  typedef std::function<FT(const Point_2&)>           Delta;
 
   Delta  delta() const { return delta_; }
   Less_x_2  less_x_2_object() const { return r_.less_x_2_object(); }
   Less_y_2  less_y_2_object() const { return r_.less_y_2_object(); }
   Greater_x_2  greater_x_2_object() const
-  { return boost::bind(less_x_2_object(),_2,_1); }
+  { return [this](const Point_2& p1, const Point_2& p2){ return this->less_x_2_object()(p2, p1); }; }
   Greater_y_2  greater_y_2_object() const
-  { return boost::bind(less_y_2_object(),_2,_1); }
+  { return [this](const Point_2& p1, const Point_2& p2){ return this->less_y_2_object()(p2, p1); }; }
   Infinity_distance_2  distance() const
   { return r_.infinity_distance_2_object(); }
   Construct_vertex_2  construct_vertex_2_object() const
@@ -277,7 +269,7 @@ struct Rectangular_3_center_2_type2_operations_base {
 public:
 
   Rectangular_3_center_2_type2_operations_base(R& r, const Point_2& p)
-    : r_(r), delta_(boost::bind(r.infinity_distance_2_object(), p, _1))
+    : r_(r), delta_([&r, &p](const Point_2& q){ return r.infinity_distance_2_object()(p, q); })
   {}
 
 };
@@ -846,10 +838,6 @@ rectangular_3_center_2_type2(
   Operations op)
 {
   BOOST_USING_STD_MAX();
-  using std::less;
-  using std::greater;
-  using std::greater_equal;
-  using std::not_equal_to;
   using std::logical_and;
   using std::max_element;
   using std::find_if;
@@ -894,7 +882,7 @@ rectangular_3_center_2_type2(
   {
     // First try whether the best radius so far can be reached at all
     RandomAccessIterator m =
-      partition(f, l, boost::bind(greater< FT >(), rad, boost::bind(op.delta(), _1)));
+      partition(f, l, [&rad, &op](const Point& p){ return rad > op.delta()(p); });
     IP pos = min_max_element(m, l, op.compare_x(), op.compare_y());
     // extreme points of the two other squares
     Point q_t =
@@ -905,11 +893,11 @@ rectangular_3_center_2_type2(
       op.place_y_square(op.place_y_square(Q_r_empty, Q_r, *pos.second, r),
                         r,
                         rad);
-    boost::function1<bool,FT> le_rad = boost::bind(greater_equal<FT>(), rad, _1);
+    auto le_rad = [&rad](const FT& v){return rad >= v;};
     RandomAccessIterator b1 =
-      partition(m, l, boost::bind(le_rad, boost::bind(op.distance(), q_t, _1)));
+      partition(m, l, [&le_rad, &op, q_t](const Point& p){ return le_rad(op.distance()(q_t, p)); });
     RandomAccessIterator b2 =
-      partition(b1, l, boost::bind(le_rad, boost::bind(op.distance(), q_r, _1)));
+      partition(b1, l, [&le_rad, &op, q_r](const Point& p){ return le_rad(op.distance()(q_r, p)); });
 
     if (b2 != l)
       return o;
@@ -921,8 +909,7 @@ rectangular_3_center_2_type2(
   while (e - s > 6) {
     std::ptrdiff_t cutoff = (e - s) / 2;
     RandomAccessIterator m = s + cutoff - 1;
-    std::nth_element(s, m, e,
-                boost::bind(less<FT>(), boost::bind(op.delta(), _1), boost::bind(op.delta(), _2)));
+    std::nth_element(s, m, e, [&op](const Point& p1, const Point& p2){ return op.delta()(p1) < op.delta()(p2); });
 
     // step (b)
     IP pos = min_max_element(m + 1, e, op.compare_x(), op.compare_y());
@@ -935,13 +922,11 @@ rectangular_3_center_2_type2(
     Point q_r = op.place_y_square(q_r_afap, r, op.delta()(*m));
 
     // check for covering
-    boost::function1<bool,FT>
-      le_delta_m = boost::bind(greater_equal<FT>(), op.delta()(*m), _1);
+    auto le_delta_m = [&op, m](const FT& v){ return op.delta()(*m) >= v; };
     RandomAccessIterator b1 =
-      partition(m + 1, e,
-                boost::bind(le_delta_m, boost::bind(op.distance(), q_t, _1)));
+      partition(m + 1, e, [&le_delta_m, &op, & q_t](const Point& p){ return le_delta_m(op.distance()(q_t, p)); });
     RandomAccessIterator b2 =
-      partition(b1, e, boost::bind(le_delta_m, boost::bind(op.distance(), q_r, _1)));
+      partition(b1, e, [&le_delta_m, &op, & q_r](const Point& p){ return le_delta_m(op.distance()(q_r, p)); });
 
     if (b2 != e)
       s = m;
@@ -960,7 +945,7 @@ rectangular_3_center_2_type2(
     std::ptrdiff_t cutoff = (e - s) / fraction;
     RandomAccessIterator m = s + cutoff - 1;
     std::nth_element(s, m, e,
-                boost::bind(less<FT>(), boost::bind(op.delta(), _1), boost::bind(op.delta(), _2)));
+                [&op](const Point& p1, const Point& p2){ return op.delta()(p1) < op.delta()(p2); });
 
     // step (b)
     IP pos = min_max_element(m + 1, e, op.compare_x(), op.compare_y());
@@ -1007,16 +992,16 @@ rectangular_3_center_2_type2(
     // partition the range [m+1, e) into ranges
     // [m+1, b1), [b1, b2),   [b2, b3) and [b3, e)
     //     R      G cap q_t  G cap q_r      none
-    boost::function1<bool,FT>
-      le_delta_m = boost::bind(greater_equal<FT>(), op.delta()(*m), _1);
+    auto le_delta_m = [&op, m](const FT& v){ return op.delta()(*m) >= v; };
     RandomAccessIterator b2 =
-      partition(m + 1, e, boost::bind(le_delta_m, boost::bind(op.distance(), q_t, _1)));
+      partition(m + 1, e,
+                [&le_delta_m, &op, &q_t](const Point& p) { return le_delta_m(op.distance()(q_t, p)); });
     RandomAccessIterator b1 =
       partition(m + 1, b2,
-                boost::bind(le_delta_m, boost::bind(op.distance(), q_r, _1)));
+                [&le_delta_m, &op, &q_r](const Point& p) { return le_delta_m(op.distance()(q_r, p)); });
     RandomAccessIterator b3 =
-      partition(b2, e, boost::bind(le_delta_m, boost::bind(op.distance(), q_r, _1)));
-
+      partition(b2, e,
+                [&le_delta_m, &op, &q_r](const Point& p) { return le_delta_m(op.distance()(q_r, p)); });
 
     // step (c)
     if (b3 != e ||
@@ -1100,9 +1085,8 @@ rectangular_3_center_2_type2(
       // step 1
       RandomAccessIterator s_m = s_b + (s_e - s_b - 1) / 2;
       std::nth_element(s_b, s_m, s_e,
-                  boost::bind(less<FT>(),
-                       boost::bind(op.delta(), _1),
-                       boost::bind(op.delta(), _2)));
+                       [&op](const Point& p1, const Point& p2)
+                       { return op.delta()(p1) < op.delta()(p2); });
 
       // step 2 (as above)
       Point q_t_m = q_t_afap;
@@ -1138,14 +1122,17 @@ CGAL_3CENTER_REPEAT_CHECK:
     // partition the range [s_b+1, e) into ranges
     // [s_b+1, b1), [b1, b2),   [b2, b3) and [b3, e)
     //     R      G cap q_t  G cap q_r      none
-    boost::function1<bool,FT>
-      le_delta_sb = boost::bind(greater_equal<FT>(), op.delta()(*s_b), _1);
-    b2 = partition(s_b + 1, e, boost::bind(le_delta_sb,
-                                    boost::bind(op.distance(), q_t, _1)));
-    b1 = partition(s_b + 1, b2, boost::bind(le_delta_sb,
-                                     boost::bind(op.distance(), q_r, _1)));
+    auto le_delta_sb = [&op, s_b](const FT& v){ return op.delta()(*s_b) >= v;} ;
+
+    b2 = partition(s_b + 1, e,
+                   [&le_delta_sb, &op, &q_t](const Point& p)
+                   { return le_delta_sb(op.distance()(q_t, p)); });
+    b1 = partition(s_b + 1, b2,
+                   [&le_delta_sb, &op, &q_r](const Point& p)
+                   { return le_delta_sb(op.distance()(q_r, p)); });
     b3 = partition(b2, e,
-                   boost::bind(le_delta_sb, boost::bind(op.distance(), q_r, _1)));
+                   [&le_delta_sb, &op, &q_r](const Point& p)
+                   { return le_delta_sb(op.distance()(q_r, p)); });
 
     if (b3 != e ||
         (!Q_t_empty && op.compute_x_distance(q_t, Q_t) > op.delta()(*s_b)) ||
@@ -1183,7 +1170,7 @@ CGAL_3CENTER_REPEAT_CHECK:
         std::vector< Point > tmppts(f, l);
         RandomAccessIterator ii =
           partition(tmppts.begin(), tmppts.end(),
-                    boost::bind(le_delta_sb, boost::bind(op.delta(), _1)));
+                    [&le_delta_sb, &op](const FT& v){ return le_delta_sb(op.delta()(v)); });
         IP tmppos = min_max_element(ii, tmppts.end(),
                                     op.compare_x(), op.compare_y());
         )
@@ -1228,12 +1215,8 @@ CGAL_3CENTER_REPEAT_CHECK:
     // we have to take the next smaller radius
     RandomAccessIterator next =
       max_element_if(s, s_b,
-                     boost::bind(less<FT>(),
-                          boost::bind(op.delta(), _1),
-                          boost::bind(op.delta(), _2)),
-                     boost::bind(not_equal_to<FT>(),
-                          op.delta()(*s_b),
-                          boost::bind(op.delta(), _1)));
+                     [&op](const Point& p1, const Point& p2){ return op.delta()(p1) < op.delta()(p2); },
+                     [&op, s_b](const Point& p){ return op.delta()(*s_b) != op.delta()(p); });
     rho_max = op.delta()(*s_b);
     q_t_at_rho_max = q_t, q_r_at_rho_max = q_r;
     CGAL_optimisation_assertion(op.delta()(*next) < op.delta()(*s_b));
@@ -1243,14 +1226,13 @@ CGAL_3CENTER_REPEAT_CHECK:
     q_r = op.place_y_square(q_r_afap, r, op.delta()(*next));
 
     // again check for covering
-    boost::function1<bool,FT>
-      le_delta_next = boost::bind(greater_equal<FT>(), op.delta()(*next), _1);
+    auto le_delta_next = [&op, next](const FT& v){ return op.delta()(*next) >= v; };
     b2 = partition(s_b, e,
-                   boost::bind(le_delta_next, boost::bind(op.distance(), q_t, _1)));
+                   [&op, &le_delta_next, &q_t](const Point& p){ return le_delta_next( op.distance()(q_t,p) ); });
     b1 = partition(s_b, b2,
-                   boost::bind(le_delta_next, boost::bind(op.distance(), q_r, _1)));
+                   [&op, &le_delta_next, &q_r](const Point& p){ return le_delta_next( op.distance()(q_r,p) ); });
     b3 = partition(b2, e,
-                   boost::bind(le_delta_next, boost::bind(op.distance(), q_r, _1)));
+                   [&op, &le_delta_next, &q_r](const Point& p){ return le_delta_next( op.distance()(q_r,p) ); });
 
     if (b3 != e ||
         (!Q_t_empty && op.compute_x_distance(q_t, Q_t) > op.delta()(*next)) ||
@@ -1304,7 +1286,7 @@ CGAL_3CENTER_REPEAT_CHECK:
     Point q_t_afap = op.place_x_square(Q_t_empty, Q_t, r);
     Point q_r_afap = op.place_y_square(Q_r_empty, Q_r, r);
     if (s != e) {
-      sort(s, e, boost::bind(less<FT>(), boost::bind(op.delta(), _1), boost::bind(op.delta(), _2)));
+      sort(s, e, [&op](const Point& p1, const Point& p2){ return op.delta()(p1) < op.delta()(p2);});
       rho_max = op.delta()(*--t);
     } else
       rho_max = rho_min;
@@ -1347,16 +1329,15 @@ CGAL_3CENTER_REPEAT_CHECK:
       q_r = op.place_y_square(q_r_afap, r, try_rho);
 
       // check for covering
-      boost::function1<bool,FT>
-        greater_rho_max = boost::bind(less<FT>(), try_rho, _1);
+      auto greater_rho_max = [&try_rho](const FT& v){ return try_rho < v; };
       if ((!Q_t_empty && op.compute_x_distance(q_t, Q_t) > try_rho) ||
           (!Q_r_empty && op.compute_y_distance(q_r, Q_r) > try_rho) ||
           e != find_if(
             t + 1,
             e,
-            boost::bind(logical_and<bool>(),
-                 boost::bind(greater_rho_max, boost::bind(op.distance(), q_t, _1)),
-                 boost::bind(greater_rho_max, boost::bind(op.distance(), q_r, _1)))))
+            [&greater_rho_max, &q_t, &q_r, &op](const Point& p)
+            { return greater_rho_max(op.distance()(q_t,p)) &&
+                     greater_rho_max(op.distance()(q_r,p)); }))
         {
           rho_min = try_rho;
           q_t_q_r_cover_at_rho_min = 0;
@@ -1393,17 +1374,15 @@ CGAL_3CENTER_REPEAT_CHECK:
   CGAL_optimisation_assertion(rho_min >= 0);
   FT rad_2 = q_t_q_r_cover_at_rho_min;
   if (s_at_rho_min != e_at_rho_min) {
-    boost::function1<FT,Point>
-      mydist = boost::bind(Min<FT>(),
-                    boost::bind(op.distance(), q_t_at_rho_min, _1),
-                    boost::bind(op.distance(), q_r_at_rho_min, _1));
+    auto mydist = [&q_t_at_rho_min, &q_r_at_rho_min, &op](const Point& p)
+                  { return Min<FT>()( op.distance()(q_t_at_rho_min, p),
+                                      op.distance()(q_r_at_rho_min, p)); };
     rad_2 =
       max BOOST_PREVENT_MACRO_SUBSTITUTION (
         rad_2,
         mydist(*max_element(s_at_rho_min, e_at_rho_min,
-                            boost::bind(less< FT >(),
-                                 boost::bind(mydist, _1),
-                                 boost::bind(mydist, _2)))));
+                            [&mydist](const Point& p1, const Point& p2)
+                            { return mydist(p1) < mydist(p2); })));
   }
   CGAL_optimisation_assertion(rad_2 == 0 || rad_2 > rho_min);
 
diff --git a/Bounding_volumes/include/CGAL/rectangular_p_center_2.h b/Bounding_volumes/include/CGAL/rectangular_p_center_2.h
index 6034a94f173..1d680e9572f 100644
--- a/Bounding_volumes/include/CGAL/rectangular_p_center_2.h
+++ b/Bounding_volumes/include/CGAL/rectangular_p_center_2.h
@@ -316,7 +316,6 @@ rectangular_p_center_2_matrix_search(
   const Traits& t)
 {
   typedef typename Traits::FT FT;
-  using std::minus;
 
   return rectangular_p_center_2_matrix_search(
     f,
@@ -325,7 +324,9 @@ rectangular_p_center_2_matrix_search(
     r,
     pf,
     t,
-    boost::bind(Max<FT>(), 0, boost::bind(minus<FT>(), _1, _2)));
+    std::function<FT(const FT&, const FT&)>(
+      [](const FT& a, const FT& b) { return Max<FT>()(0, std::minus<FT>()(a,b)); }
+    ));
 
 } // Pcenter_matrix_search( ... )
 
diff --git a/Convex_hull_2/examples/Convex_hull_2/ch_graham_anderson.cpp b/Convex_hull_2/examples/Convex_hull_2/ch_graham_anderson.cpp
index bc98d8c91ed..019ede1ed24 100644
--- a/Convex_hull_2/examples/Convex_hull_2/ch_graham_anderson.cpp
+++ b/Convex_hull_2/examples/Convex_hull_2/ch_graham_anderson.cpp
@@ -2,7 +2,6 @@
 #include <CGAL/ch_graham_andrew.h>
 #include <vector>
 #include <algorithm>
-#include <boost/bind.hpp>
 
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
@@ -23,7 +22,8 @@ ch_graham_anderson( InputIterator  first, InputIterator  beyond,
   std::vector< Point_2 >  V (first, beyond);
   typename std::vector< Point_2 >::iterator it =
                std::min_element(V.begin(), V.end(), Less_xy_2());
-  std::sort( V.begin(), V.end(), boost::bind(Less_rotate_ccw_2(), *it, _1, _2) );
+  const Point_2 p = *it;
+  std::sort( V.begin(), V.end(), [&p](const Point_2& p1, const Point_2& p2){return Less_rotate_ccw_2()(p, p1, p2);} );
   if ( *(V.begin()) != *(V.rbegin()) )
   {
     result = CGAL::ch_graham_andrew_scan( V.begin(), V.end(), result, ch_traits);
diff --git a/Convex_hull_2/include/CGAL/Convex_hull_2/ch_akl_toussaint_impl.h b/Convex_hull_2/include/CGAL/Convex_hull_2/ch_akl_toussaint_impl.h
index c06435efcef..0fdc1a0367b 100644
--- a/Convex_hull_2/include/CGAL/Convex_hull_2/ch_akl_toussaint_impl.h
+++ b/Convex_hull_2/include/CGAL/Convex_hull_2/ch_akl_toussaint_impl.h
@@ -25,7 +25,6 @@
 #include <CGAL/ch_graham_andrew.h>
 #include <CGAL/algorithm.h>
 #include <CGAL/IO/Tee_for_output_iterator.h>
-#include <boost/bind.hpp>
 #include <CGAL/tuple.h>
 #include <CGAL/utility.h>
 #include <iterator>
@@ -296,9 +295,11 @@ ch_akl_toussaint(ForwardIterator first, ForwardIterator last,
   std::sort( std::next(region2.begin() ), region2.end(),
              ch_traits.less_xy_2_object() );
   std::sort( std::next(region3.begin() ), region3.end(),
-             boost::bind(ch_traits.less_xy_2_object(), _2, _1) );
+             [&ch_traits](const Point_2& p1, const Point_2& p2)
+             { return ch_traits.less_xy_2_object()(p2, p1); });
   std::sort( std::next(region4.begin() ), region4.end(),
-             boost::bind(ch_traits.less_xy_2_object(), _2, _1) );
+             [&ch_traits](const Point_2& p1, const Point_2& p2)
+             { return ch_traits.less_xy_2_object()(p2, p1); });
 
   if (! equal_points(*w,*s) )
   {
diff --git a/Convex_hull_2/include/CGAL/Convex_hull_2/ch_bykat_impl.h b/Convex_hull_2/include/CGAL/Convex_hull_2/ch_bykat_impl.h
index 32334e4e7f9..bcd0f1b13c0 100644
--- a/Convex_hull_2/include/CGAL/Convex_hull_2/ch_bykat_impl.h
+++ b/Convex_hull_2/include/CGAL/Convex_hull_2/ch_bykat_impl.h
@@ -26,7 +26,6 @@
 #include <CGAL/algorithm.h>
 #include <list>
 #include <algorithm>
-#include <boost/bind.hpp>
 
 namespace CGAL {
 template <class InputIterator, class OutputIterator, class Traits>
@@ -43,7 +42,6 @@ ch_bykat(InputIterator first, InputIterator last,
   typedef typename Traits::Equal_2                         Equal_2;
 
   Left_turn_2 left_turn    = ch_traits.left_turn_2_object();
-  Less_dist   less_dist    = ch_traits.less_signed_distance_to_line_2_object();
   Equal_2     equal_points = ch_traits.equal_2_object();
 
   if (first == last) return result;
@@ -77,18 +75,23 @@ ch_bykat(InputIterator first, InputIterator last,
   H.push_back( a );
   L.push_back( P.begin() );
   R.push_back( l = std::partition(P.begin(), P.end(),
-                                  boost::bind(left_turn, boost::cref(a), boost::cref(b), _1)));
-  r = std::partition( l, P.end(), boost::bind(left_turn, boost::cref(b), boost::cref(a), _1));
+                                  [&left_turn, &a, &b](const Point_2& p){ return left_turn(a,b,p); }) );
+  r = std::partition( l, P.end(), [&left_turn, &a, &b](const Point_2& p){ return left_turn(b,a,p); });
 
   for (;;)
   {
+    //  This functor must be in the for loop so that the Convex_hull_constructive traits_2 works correctly
+    Less_dist less_dist = ch_traits.less_signed_distance_to_line_2_object();
       if ( l != r)
       {
-        Point_2 c = *std::min_element( l, r, boost::bind(less_dist, boost::cref(a), boost::cref(b), _1, _2));
+        Point_2 c = *std::min_element( l, r, [&less_dist,&a,&b](const Point_2&p1, const Point_2& p2)
+                                             { return less_dist(a, b, p1, p2); });
           H.push_back( b );
           L.push_back( l );
-          R.push_back( l = std::partition(l, r, boost::bind(left_turn, boost::cref(b), boost::cref(c), _1)));
-          r = std::partition(l, r, boost::bind(left_turn, boost::cref(c), boost::cref(a), _1));
+          R.push_back( l = std::partition(l, r, [&left_turn,&c,&b](const Point_2&p)
+                                                { return left_turn(b, c, p); }));
+          r = std::partition(l, r, [&left_turn,&c,&a](const Point_2&p)
+                                   { return left_turn(c, a, p); });
           b = c;
       }
       else
@@ -173,8 +176,10 @@ ch_bykat_with_threshold(InputIterator   first, InputIterator last,
   H.push_back( a );
   L.push_back( Pbegin );
   Left_turn_2 left_turn = ch_traits.left_turn_2_object();
-  R.push_back( l = std::partition( Pbegin, Pend,  boost::bind(left_turn, boost::cref(a), boost::cref(b), _1)));
-  r = std::partition( l, Pend, boost::bind(left_turn, boost::cref(b), boost::cref(a), _1));
+  R.push_back( l = std::partition( Pbegin, Pend,  [&left_turn,&a,&b](const Point_2&p)
+                                                  { return left_turn(a, b, p); }));
+  r = std::partition( l, Pend, [&left_turn,&a,&b](const Point_2&p)
+                               { return left_turn(b, a, p); });
 
   Less_dist less_dist = ch_traits.less_signed_distance_to_line_2_object();
   for (;;)
@@ -183,11 +188,14 @@ ch_bykat_with_threshold(InputIterator   first, InputIterator last,
       {
           if ( r-l > CGAL_ch_THRESHOLD )
           {
-              Point_2 c = *std::min_element( l, r, boost::bind(less_dist, boost::cref(a), boost::cref(b), _1, _2));
+              Point_2 c = *std::min_element( l, r, [&less_dist,&a,&b](const Point_2&p1, const Point_2& p2)
+                                                   { return less_dist(a, b, p1, p2); });
               H.push_back( b );
               L.push_back( l );
-              R.push_back( l = std::partition(l, r, boost::bind(left_turn, boost::cref(b), boost::cref(c), _1)));
-              r = std::partition(l, r, boost::bind(left_turn, boost::cref(c), boost::cref(a), _1));
+              R.push_back( l = std::partition(l, r, [&left_turn,&c,&b](const Point_2&p)
+                                                    { return left_turn(b, c, p); }));
+              r = std::partition(l, r, [&left_turn,&a,&c](const Point_2&p)
+                                       { return left_turn(c, a, p); });
               b = c;
           }
           else
@@ -201,8 +209,8 @@ ch_bykat_with_threshold(InputIterator   first, InputIterator last,
               }
               else
               {
-                std::sort(std::next(l), r,
-                            boost::bind(ch_traits.less_xy_2_object(), _2, _1) );
+                std::sort(std::next(l), r, [&ch_traits](const Point_2&p1, const Point_2& p2)
+                                           { return ch_traits.less_xy_2_object()(p2, p1); });
               }
               ch__ref_graham_andrew_scan(l, std::next(r), res, ch_traits);
               std::swap( a, *l);
diff --git a/Convex_hull_2/include/CGAL/Convex_hull_2/ch_eddy_impl.h b/Convex_hull_2/include/CGAL/Convex_hull_2/ch_eddy_impl.h
index 2a7b29c2ddd..19df2467657 100644
--- a/Convex_hull_2/include/CGAL/Convex_hull_2/ch_eddy_impl.h
+++ b/Convex_hull_2/include/CGAL/Convex_hull_2/ch_eddy_impl.h
@@ -25,7 +25,6 @@
 #include <CGAL/algorithm.h>
 #include <list>
 #include <algorithm>
-#include <boost/bind.hpp>
 
 namespace CGAL {
 
@@ -45,7 +44,8 @@ ch__recursive_eddy(List& L,
 
   CGAL_ch_precondition( \
     std::find_if(a_it, b_it, \
-                 boost::bind(left_turn, *b_it, *a_it, _1)) \
+                 [&left_turn, a_it, b_it](const Point_2& p)
+                 { return left_turn(*b_it, *a_it, p); }) \
     != b_it );
 
 
@@ -53,11 +53,14 @@ ch__recursive_eddy(List& L,
   Less_dist less_dist = ch_traits.less_signed_distance_to_line_2_object();
   ListIterator
       c_it = std::min_element( f_it, b_it,  // max before
-                               boost::bind(less_dist, *a_it, *b_it, _1, _2));
+                               [&less_dist, a_it, b_it](const Point_2& p1, const Point_2& p2)
+                               { return less_dist(*a_it, *b_it, p1, p2); });
   Point_2 c = *c_it;
 
-  c_it = std::partition(f_it, b_it, boost::bind(left_turn, c, *a_it, _1));
-  f_it = std::partition(c_it, b_it, boost::bind(left_turn, *b_it, c, _1));
+  c_it = std::partition(f_it, b_it, [&left_turn, &c, a_it](const Point_2& p)
+                                    {return left_turn(c, *a_it, p);});
+  f_it = std::partition(c_it, b_it, [&left_turn, &c, b_it](const Point_2& p)
+                                    {return left_turn(*b_it, c, p);});
   c_it = L.insert(c_it, c);
   L.erase( f_it, b_it );
 
@@ -104,7 +107,8 @@ ch_eddy(InputIterator first, InputIterator last,
   L.erase(e);
 
   e = std::partition(L.begin(), L.end(),
-                     boost::bind(left_turn, ep, wp, _1) );
+                     [&left_turn, &wp, &ep](const Point_2& p)
+                     {return left_turn(ep, wp, p);} );
   L.push_front(wp);
   e = L.insert(e, ep);
 
@@ -112,7 +116,8 @@ ch_eddy(InputIterator first, InputIterator last,
   {
       ch__recursive_eddy( L, L.begin(), e, ch_traits);
   }
-  w = std::find_if( e, L.end(), boost::bind(left_turn, wp, ep, _1) );
+  w = std::find_if( e, L.end(), [&left_turn, &wp, &ep](const Point_2& p)
+                                { return left_turn(wp, ep, p); });
   if ( w == L.end() )
   {
       L.erase( ++e, L.end() );
diff --git a/Convex_hull_2/include/CGAL/Convex_hull_2/ch_jarvis_impl.h b/Convex_hull_2/include/CGAL/Convex_hull_2/ch_jarvis_impl.h
index c5d3309cbd0..65b2d1fc29c 100644
--- a/Convex_hull_2/include/CGAL/Convex_hull_2/ch_jarvis_impl.h
+++ b/Convex_hull_2/include/CGAL/Convex_hull_2/ch_jarvis_impl.h
@@ -24,7 +24,6 @@
 #include <CGAL/Convex_hull_2/ch_assertions.h>
 #include <CGAL/ch_selected_extreme_points_2.h>
 #include <algorithm>
-#include <boost/bind.hpp>
 
 namespace CGAL {
 
@@ -65,7 +64,8 @@ ch_jarvis_march(ForwardIterator first, ForwardIterator last,
       Point previous_point = start_p; )
 
   ForwardIterator it = std::min_element( first, last,
-                                         boost::bind(rotation_predicate, boost::cref(start_p), _1, _2) );
+                                         [&start_p, &rotation_predicate](const Point& p1, const Point& p2)
+                                         {return rotation_predicate(start_p, p1, p2);} );
   while (! equal_points(*it, stop_p) )
   {
       CGAL_ch_exactness_assertion( \
@@ -80,7 +80,8 @@ ch_jarvis_march(ForwardIterator first, ForwardIterator last,
           constructed_points <= count_points + 1 );
 
       it = std::min_element( first, last,
-                             boost::bind(rotation_predicate, *it, _1, _2) );
+                             [it, &rotation_predicate](const Point& p1, const Point& p2)
+                             {return rotation_predicate(*it, p1, p2);} );
   }
   CGAL_ch_postcondition( \
       is_ccw_strongly_convex_2( res.output_so_far_begin(), \
diff --git a/Convex_hull_3/include/CGAL/convex_hull_3.h b/Convex_hull_3/include/CGAL/convex_hull_3.h
index 995c978abd8..cf35a6d78af 100644
--- a/Convex_hull_3/include/CGAL/convex_hull_3.h
+++ b/Convex_hull_3/include/CGAL/convex_hull_3.h
@@ -43,7 +43,6 @@
 #include <CGAL/boost/graph/named_params_helper.h>
 #include <CGAL/is_iterator.h>
 
-#include <boost/bind.hpp>
 #include <boost/next_prior.hpp>
 #include <boost/type_traits/is_floating_point.hpp>
 #include <boost/type_traits/is_same.hpp>
@@ -576,7 +575,8 @@ farthest_outside_point(Face_handle f, std::list<Point>& outside_set,
    Outside_set_iterator farthest_it =
           std::max_element(outside_set.begin(),
                            outside_set.end(),
-                           boost::bind(less_dist_to_plane, plane, _1, _2));
+                           [&less_dist_to_plane,&plane](const Point& p1, const Point& p2)
+                           { return less_dist_to_plane(plane, p1, p2); });
    return farthest_it;
 }
 
@@ -795,8 +795,10 @@ ch_quickhull_face_graph(std::list<typename Traits::Point_3>& points,
   // plane.
   std::pair<P3_iterator, P3_iterator> min_max;
   min_max = CGAL::min_max_element(points.begin(), points.end(),
-                                  boost::bind(compare_dist, plane, _1, _2),
-                                  boost::bind(compare_dist, plane, _1, _2));
+                                  [&compare_dist, &plane](const Point_3& p1, const Point_3& p2)
+                                  { return compare_dist(plane, p1, p2); },
+                                  [&compare_dist, &plane](const Point_3& p1, const Point_3& p2)
+                                  { return compare_dist(plane, p1, p2); });
   P3_iterator max_it;
   if (coplanar(*point1_it, *point2_it, *point3_it, *min_max.second))
   {
@@ -929,8 +931,10 @@ convex_hull_3(InputIterator first, InputIterator beyond,
      Less_dist less_dist = traits.less_distance_to_point_3_object();
      P3_iterator_pair endpoints =
       min_max_element(points.begin(), points.end(),
-                      boost::bind(less_dist, *points.begin(), _1, _2),
-                      boost::bind(less_dist, *points.begin(), _1, _2));
+                      [&points, &less_dist](const Point_3& p1, const Point_3& p2)
+                      { return less_dist(*points.begin(), p1, p2); },
+                      [&points, &less_dist](const Point_3& p1, const Point_3& p2)
+                      { return less_dist(*points.begin(), p1, p2); });
 
      typename Traits::Construct_segment_3 construct_segment =
             traits.construct_segment_3_object();
@@ -1032,8 +1036,10 @@ void convex_hull_3(InputIterator first, InputIterator beyond,
     Less_dist less_dist = traits.less_distance_to_point_3_object();
     P3_iterator_pair endpoints =
     min_max_element(points.begin(), points.end(),
-                    boost::bind(less_dist, *points.begin(), _1, _2),
-                    boost::bind(less_dist, *points.begin(), _1, _2));
+                    [&points, &less_dist](const Point_3& p1, const Point_3& p2)
+                    { return less_dist(*points.begin(), p1, p2); },
+                    [&points, &less_dist](const Point_3& p1, const Point_3& p2)
+                    { return less_dist(*points.begin(), p1, p2); });
     Convex_hull_3::internal::add_isolated_points(*endpoints.first, polyhedron);
     Convex_hull_3::internal::add_isolated_points(*endpoints.second, polyhedron);
     return;
diff --git a/Generator/include/CGAL/random_convex_set_2.h b/Generator/include/CGAL/random_convex_set_2.h
index a118503b4c3..f8189fcea1c 100644
--- a/Generator/include/CGAL/random_convex_set_2.h
+++ b/Generator/include/CGAL/random_convex_set_2.h
@@ -23,7 +23,6 @@
 #include <algorithm>
 #include <numeric>
 #include <CGAL/Random_convex_set_traits_2.h>
-#include <boost/functional.hpp>
 
 namespace CGAL {
 
@@ -80,7 +79,7 @@ random_convex_set_2( std::size_t n,
     points.begin(),
     points.end(),
     points.begin(),
-    boost::bind2nd( Sum(), scale( centroid, FT( -1))));
+    [&centroid, &sum, &scale](const Point_2& p) { return sum(p, scale(centroid, FT( -1))); });
 
   // sort them according to their direction's angle
   // w.r.t. the positive x-axis:
@@ -102,8 +101,9 @@ random_convex_set_2( std::size_t n,
     points.begin(),
     points.end(),
     points.begin(),
-    boost::bind2nd( Sum(), sum( centroid,
-                         scale( new_centroid, FT( -1)))));
+    [&centroid, &new_centroid, &sum, &scale](const Point_2& p)
+    {return sum(p, sum( centroid, scale(new_centroid, FT( -1)))); }
+  );
 
   // compute maximal coordinate:
   FT maxcoord( max_coordinate(
@@ -118,7 +118,7 @@ random_convex_set_2( std::size_t n,
     points.begin(),
     points.end(),
     o,
-    boost::bind2nd( Scale(), FT( pg.range()) / maxcoord));
+    [&pg, &maxcoord, &scale](const Point_2& p){ return scale(p, FT( pg.range()) / maxcoord); });
 
 } // random_convex_set_2( n, o, pg, t)
 
diff --git a/Inscribed_areas/include/CGAL/Extremal_polygon_traits_2.h b/Inscribed_areas/include/CGAL/Extremal_polygon_traits_2.h
index 64f3ac2dd1d..a85e9a76f47 100644
--- a/Inscribed_areas/include/CGAL/Extremal_polygon_traits_2.h
+++ b/Inscribed_areas/include/CGAL/Extremal_polygon_traits_2.h
@@ -19,7 +19,6 @@
 #include <CGAL/Optimisation/assertions.h>
 #include <CGAL/squared_distance_2.h>
 #include <CGAL/Polygon_2.h>
-#include <boost/bind.hpp>
 #include <boost/function.hpp>
 
 namespace CGAL {
@@ -47,7 +46,7 @@ struct Extremal_polygon_area_traits_2 {
   };
 
   typedef Kgon_triangle_area                         Baseop;
-  typedef boost::function2<FT,Point_2,Point_2>       Operation;
+  typedef std::function<FT(const Point_2& , const Point_2&)>       Operation;
 
   Extremal_polygon_area_traits_2() {}
   Extremal_polygon_area_traits_2(const K& k_) : k(k_) {}
@@ -60,7 +59,7 @@ struct Extremal_polygon_area_traits_2 {
 
   Operation
   operation( const Point_2& p) const
-  { return boost::bind(Baseop(k), _1, _2, p); }
+  { return [&p, this](const Point_2& p1, const Point_2& p2){ return Baseop(this->k)(p1, p2, p);}; }
 
   template < class RandomAccessIC, class OutputIterator >
   OutputIterator
@@ -184,7 +183,7 @@ struct Extremal_polygon_perimeter_traits_2 {
   };
 
   typedef Kgon_triangle_perimeter                    Baseop;
-  typedef boost::function2<FT,Point_2,Point_2>       Operation;
+  typedef std::function<FT(const Point_2&, const Point_2&)>       Operation;
 
   Extremal_polygon_perimeter_traits_2() {}
   Extremal_polygon_perimeter_traits_2(const K& k_) : k(k_) {}
@@ -197,7 +196,7 @@ struct Extremal_polygon_perimeter_traits_2 {
 
   Operation
   operation( const Point_2& p) const
-  { return boost::bind(Baseop(k), _1, _2, p); }
+  { return [this,&p](const Point_2& p1, const Point_2& p2){ return Baseop(this->k)(p1,p2,p); }; }
 
   template < class RandomAccessIC, class OutputIterator >
   OutputIterator
@@ -234,10 +233,9 @@ struct Extremal_polygon_perimeter_traits_2 {
       max_element(
         points_begin + 1,
         points_end,
-        boost::bind(
-          less< FT >(),
-          boost::bind(operation(points_begin[0]), _1, points_begin[0]),
-          boost::bind(operation(points_begin[0]), _2, points_begin[0]))));
+        [points_begin, this](const Point_2& p1, const Point_2& p2)
+        { return less< FT >()( this->operation(points_begin[0])(p1, points_begin[0]),
+                               this->operation(points_begin[0])(p2, points_begin[0])); } ));
 
     // give result:
     max_perimeter = operation(*points_begin)(*maxi, *points_begin);
diff --git a/Inscribed_areas/include/CGAL/extremal_polygon_2.h b/Inscribed_areas/include/CGAL/extremal_polygon_2.h
index 33c6a126f9a..2cd09fdd7c0 100644
--- a/Inscribed_areas/include/CGAL/extremal_polygon_2.h
+++ b/Inscribed_areas/include/CGAL/extremal_polygon_2.h
@@ -25,8 +25,6 @@
 #include <functional>
 #include <algorithm>
 #include <CGAL/Extremal_polygon_traits_2.h>
-#include <boost/bind.hpp>
-#include <boost/bind/make_adaptable.hpp>
 namespace CGAL {
 
 //!!! This will eventually be integrated into function_objects.h
@@ -443,8 +441,8 @@ extremal_polygon_2(
     k,
     CGAL::transform_iterator(
       o,
-      boost::make_adaptable<Point_2, int>(boost::bind(Index_operator< RandomAccessIC, int, Point_2 >(),
-                                                      points_begin, _1))),
+      std::function<Point_2(int)>([points_begin](int i)
+      { return Index_operator< RandomAccessIC, int, Point_2 >()(points_begin, i); })),
     t);
 }
 
diff --git a/Matrix_search/examples/Matrix_search/sorted_matrix_search.cpp b/Matrix_search/examples/Matrix_search/sorted_matrix_search.cpp
index 6fce8f6dd10..a9c261a494c 100644
--- a/Matrix_search/examples/Matrix_search/sorted_matrix_search.cpp
+++ b/Matrix_search/examples/Matrix_search/sorted_matrix_search.cpp
@@ -38,7 +38,7 @@ int main()
   CGAL::sorted_matrix_search(
     &M, &M + 1,
     CGAL::sorted_matrix_search_traits_adaptor(
-      boost::bind2nd(std::greater_equal<Value>(), bound), M));
+      [&bound](const auto& m){ return std::greater_equal<Value>()(m, bound); }, M));
   std::cout << "Upper bound for " << bound << " is "
             << upper_bound << "." << std::endl;
 
diff --git a/Matrix_search/include/CGAL/Cartesian_matrix.h b/Matrix_search/include/CGAL/Cartesian_matrix.h
index 1bad69fc2d4..d3ad276a715 100644
--- a/Matrix_search/include/CGAL/Cartesian_matrix.h
+++ b/Matrix_search/include/CGAL/Cartesian_matrix.h
@@ -18,6 +18,7 @@
 
 #include <CGAL/basic.h>
 #include <CGAL/Optimisation/assertions.h>
+#include <type_traits>
 
 namespace CGAL {
 
@@ -26,7 +27,12 @@ template < class Operation,
            class RandomAccessIC_column >
 class Cartesian_matrix {
 public:
+
+#if CGAL_CXX17 && __has_cpp_attribute(nodiscard)
+  typedef typename std::invoke_result<Operation, typename std::iterator_traits<RandomAccessIC_row>::value_type, typename std::iterator_traits<RandomAccessIC_column>::value_type>::type Value;
+#else
   typedef typename Operation::result_type           Value;
+#endif
 
   Cartesian_matrix(RandomAccessIC_row r_f,
                    RandomAccessIC_row r_l,
diff --git a/Matrix_search/include/CGAL/Transform_iterator.h b/Matrix_search/include/CGAL/Transform_iterator.h
index 4797271123d..c3ab6e8c437 100644
--- a/Matrix_search/include/CGAL/Transform_iterator.h
+++ b/Matrix_search/include/CGAL/Transform_iterator.h
@@ -34,7 +34,6 @@ struct Transform_iterator {
   typedef std::_Unchecked_iterator_tag _Checked_iterator_category;
   typedef std::output_iterator_tag             iterator_category;
   typedef Transform_iterator< OutputIterator, Operation >   self;
-  typedef typename Operation::argument_type        argument_type;
 
   typedef typename std::iterator_traits<OutputIterator>::difference_type difference_type;
   typedef typename std::iterator_traits<OutputIterator>::value_type      value_type;
@@ -54,6 +53,7 @@ struct Transform_iterator {
 
   self& operator++( int) { return *this; }
 
+  template <typename argument_type>
   self& operator=( const argument_type& a) {
     *(o_++) = op_( a);
     return *this;
diff --git a/Matrix_search/include/CGAL/sorted_matrix_search.h b/Matrix_search/include/CGAL/sorted_matrix_search.h
index 8d2cf105ffb..dd21234b4f5 100644
--- a/Matrix_search/include/CGAL/sorted_matrix_search.h
+++ b/Matrix_search/include/CGAL/sorted_matrix_search.h
@@ -18,7 +18,6 @@
 
 #include <CGAL/basic.h>
 #include <CGAL/Optimisation/assertions.h>
-#include <boost/bind.hpp>
 #include <algorithm>
 #include <vector>
 #include <CGAL/Sorted_matrix_search_traits_adaptor.h>
@@ -254,10 +253,10 @@ sorted_matrix_search(InputIterator f, InputIterator l, Traits t)
     std::nth_element(active_cells.begin(),
                 active_cells.begin() + upper_median_rank,
                 active_cells.end(),
-                boost::bind(
-                  t.compare_strictly(),
-                  boost::bind(Cell_min<Cell>(), _1),
-                  boost::bind(Cell_min<Cell>(), _2)));
+                [&t](const Cell& c1, const Cell& c2)
+                {
+                  return t.compare_strictly()(Cell_min<Cell>()(c1), Cell_min<Cell>()(c2));
+                });
 
     Cell_iterator lower_median_cell =
       active_cells.begin() + upper_median_rank;
@@ -267,10 +266,10 @@ sorted_matrix_search(InputIterator f, InputIterator l, Traits t)
     std::nth_element(active_cells.begin(),
                 active_cells.begin() + lower_median_rank,
                 active_cells.end(),
-                boost::bind(
-                  t.compare_strictly(),
-                  boost::bind(Cell_max< Cell >(ccd), _1),
-                  boost::bind(Cell_max< Cell >(ccd), _2)));
+                [&t, &ccd](const Cell& c1, const Cell& c2)
+                {
+                  return t.compare_strictly()(Cell_max< Cell >(ccd)(c1), Cell_max< Cell >(ccd)(c2));
+                });
 
     Cell_iterator upper_median_cell =
       active_cells.begin() + lower_median_rank;
@@ -282,10 +281,10 @@ sorted_matrix_search(InputIterator f, InputIterator l, Traits t)
       lower_median_cell =
         find_if(active_cells.begin(),
                 active_cells.end(),
-                boost::bind(
-                  equal_to< Value >(),
-                  lower_median,
-                  boost::bind(Cell_min< Cell >(), _1)));
+                [&lower_median](const Cell& c)
+                {
+                  return equal_to< Value >()(lower_median, Cell_min< Cell >()(c));
+                });
     CGAL_optimisation_assertion(lower_median_cell != active_cells.end());
     // ------------------------------------------------------
     // test feasibility of medians and remove cells accordingly:
@@ -318,10 +317,10 @@ sorted_matrix_search(InputIterator f, InputIterator l, Traits t)
           remove_if(
             active_cells.begin() + 1,
             active_cells.end(),
-            boost::bind(
-              t.compare_non_strictly(),
-              min_median,
-              boost::bind(Cell_min< Cell >(), _1)));
+            [&t, &min_median](const Cell& c)
+            {
+              return t.compare_non_strictly()(min_median, Cell_min< Cell >()(c));
+            });
 
       } // lower_median and upper_median are feasible
       else { // lower_median is feasible, but upper_median is not
@@ -337,16 +336,11 @@ sorted_matrix_search(InputIterator f, InputIterator l, Traits t)
           remove_if(
             active_cells.begin() + 1,
             active_cells.end(),
-            boost::bind(
-              logical_or< bool >(),
-              boost::bind(
-                t.compare_non_strictly(),
-                lower_median,
-                boost::bind(Cell_min< Cell >(), _1)),
-              boost::bind(
-                  t.compare_non_strictly(),
-                  boost::bind(Cell_max< Cell >( ccd), _1),
-                  upper_median)));
+            [&t, &lower_median, &upper_median, &ccd](const Cell& c)
+            {
+              return t.compare_non_strictly()(lower_median, Cell_min< Cell >()(c)) ||
+                     t.compare_non_strictly()(Cell_max< Cell >(ccd)(c), upper_median);
+            });
 
       } // lower_median is feasible, but upper_median is not
     else
@@ -364,16 +358,11 @@ sorted_matrix_search(InputIterator f, InputIterator l, Traits t)
           remove_if(
             active_cells.begin() + 1,
             active_cells.end(),
-            boost::bind(
-              logical_or< bool >(),
-              boost::bind(
-                t.compare_non_strictly(),
-                upper_median,
-                boost::bind(Cell_min< Cell >(), _1)),
-              boost::bind(
-                t.compare_non_strictly(),
-                boost::bind(Cell_max< Cell >( ccd), _1),
-                lower_median)));
+            [&t, &lower_median, &upper_median, &ccd](const Cell& c)
+            {
+              return t.compare_non_strictly()(upper_median, Cell_min<Cell>()(c)) ||
+                     t.compare_non_strictly()(Cell_max<Cell>(ccd)(c),lower_median);
+            });
 
       } // upper_median is feasible, but lower_median is not
       else { // both upper_median and lower_median are infeasible
@@ -385,11 +374,11 @@ sorted_matrix_search(InputIterator f, InputIterator l, Traits t)
           remove_if(
             active_cells.begin(),
             active_cells.end(),
-            boost::bind(
-              t.compare_non_strictly(),
-              boost::bind(Cell_max< Cell >( ccd), _1),
-              max BOOST_PREVENT_MACRO_SUBSTITUTION ( lower_median, upper_median)));
-
+            [&t, &ccd, &lower_median, &upper_median](const Cell& c)
+            {
+              return t.compare_non_strictly()(Cell_max<Cell>( ccd)(c),
+                                              max BOOST_PREVENT_MACRO_SUBSTITUTION ( lower_median, upper_median));
+            });
       } // both upper_median and lower_median are infeasible
 
     active_cells.erase( new_end, active_cells.end());
diff --git a/Matrix_search/test/Matrix_search/sorted_matrix_search_test.cpp b/Matrix_search/test/Matrix_search/sorted_matrix_search_test.cpp
index af9d062c74f..22b2e16b82d 100644
--- a/Matrix_search/test/Matrix_search/sorted_matrix_search_test.cpp
+++ b/Matrix_search/test/Matrix_search/sorted_matrix_search_test.cpp
@@ -181,7 +181,7 @@ main( int argc, char* argv[])
       matrices.begin(),
       matrices.end(),
       sorted_matrix_search_traits_adaptor(
-        boost::bind( greater_equal< Value >(), _1, bound),
+        [&bound](const auto& m){ return greater_equal< Value >()(m, bound); },
         *(matrices.begin()))));
 
   #ifdef OUTPUT
diff --git a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h b/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h
index f810e9fc044..ab9a4844e5f 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h
@@ -50,7 +50,6 @@
 
 #include <CGAL/boost/iterator/transform_iterator.hpp>
 
-#include <boost/bind.hpp>
 #include <boost/iterator/function_output_iterator.hpp>
 #ifndef CGAL_NO_ASSERTIONS
 #  include <boost/math/special_functions/next.hpp> // for float_prior
@@ -1879,7 +1878,7 @@ next_vertex_along_curve(const Vertex_handle& start,
   adjacent_vertices.erase
     (std::remove_if(adjacent_vertices.begin(),
                     adjacent_vertices.end(),
-                    boost::bind(&Adjacent_vertices::value_type::second, _1) != curve_index),
+                    [curve_index](const auto& p){ return p.second != curve_index; }),
      adjacent_vertices.end());
   CGAL_assertion(adjacent_vertices.size() == 2);
 
diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h
index 8df5d276c14..8b0f3b1fb6c 100644
--- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h
+++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h
@@ -49,7 +49,6 @@
 
 #include <CGAL/boost/iterator/transform_iterator.hpp>
 
-#include <boost/bind.hpp>
 #include <boost/iterator/function_output_iterator.hpp>
 #ifndef CGAL_NO_ASSERTIONS
 #  include <boost/math/special_functions/next.hpp> // for float_prior
@@ -2892,7 +2891,7 @@ next_vertex_along_curve(const Vertex_handle& start,
   adjacent_vertices.erase
     (std::remove_if(adjacent_vertices.begin(),
                     adjacent_vertices.end(),
-                    boost::bind(&Adjacent_vertices::value_type::second, _1) != curve_index),
+                    [curve_index](const auto& p){ return p.second != curve_index;}),
      adjacent_vertices.end());
 
 //  typename Adjacent_vertices::const_iterator iv = adjacent_vertices.begin();
diff --git a/Polytope_distance_d/include/CGAL/Polytope_distance_d.h b/Polytope_distance_d/include/CGAL/Polytope_distance_d.h
index 3f6a4b74f93..b52d1ca7219 100644
--- a/Polytope_distance_d/include/CGAL/Polytope_distance_d.h
+++ b/Polytope_distance_d/include/CGAL/Polytope_distance_d.h
@@ -665,9 +665,8 @@ private:
   check_dimension( InputIterator first, InputIterator last)
   { return ( std::find_if
              ( first, last,
-               CGAL::compose1_1
-               ( boost::bind2nd(std::not_equal_to<int>(), d),
-                 tco.access_dimension_d_object()))
+               [this](const auto& o)
+               { return this->tco.access_dimension_d_object()(o) != this->d; })
              == last); }
 
   // compute (squared) distance
diff --git a/Polytope_distance_d/include/CGAL/all_furthest_neighbors_2.h b/Polytope_distance_d/include/CGAL/all_furthest_neighbors_2.h
index ad05b31358a..d2927f43d6b 100644
--- a/Polytope_distance_d/include/CGAL/all_furthest_neighbors_2.h
+++ b/Polytope_distance_d/include/CGAL/all_furthest_neighbors_2.h
@@ -22,7 +22,6 @@
 #include <CGAL/monotone_matrix_search.h>
 #include <CGAL/Polygon_2_algorithms.h>
 #include <algorithm>
-#include <boost/bind.hpp>
 
 namespace CGAL {
 template < class Operation, class RandomAccessIC >
@@ -116,7 +115,7 @@ all_furthest_neighbors_2( RandomAccessIC points_begin,
   return transform(v.begin(),
                    v.end(),
                    o,
-                   boost::bind(modulus<int>(), _1, number_of_points));
+                   [number_of_points](int i){ return i % number_of_points;} );
 } // all_furthest_neighbors_2( ... )
 
 
diff --git a/QP_solver/include/CGAL/QP_solution.h b/QP_solver/include/CGAL/QP_solution.h
index f1908434cc5..722acf3e7dc 100644
--- a/QP_solver/include/CGAL/QP_solution.h
+++ b/QP_solver/include/CGAL/QP_solution.h
@@ -27,8 +27,6 @@
 #include <CGAL/function_objects.h>
 #include <CGAL/Algebraic_structure_traits.h>
 #include <CGAL/QP_solver/assertions.h>
-#include <boost/bind.hpp>
-#include <boost/function.hpp>
 
 #include <CGAL/boost/iterator/counting_iterator.hpp>
 #include <CGAL/boost/iterator/transform_iterator.hpp>
@@ -75,7 +73,7 @@ public:
   typedef QP_solution_detail::Quotient_normalizer<ET>
   Quotient_normalizer; // normalizer (ET, ET) -> (ET, ET)
 
-  typedef boost::function1< Quotient<ET>, ET >
+  typedef std::function< Quotient<ET>(const ET&) >
   Quotient_maker;
 
   typedef std::vector<int>
@@ -130,10 +128,7 @@ public:
     ET n = solution_numerator();
     ET d = solution_denominator();
     return
-      boost::bind
-      (Quotient_normalizer(), boost::bind
-       (U_Quotient_creator(), _1, _2))
-      (n, d);
+      Quotient_normalizer()( U_Quotient_creator()(n,d) );
       // (solution_numerator(), solution_denominator());
   }
   virtual Quadratic_program_status status() const = 0;
@@ -164,20 +159,22 @@ public:
   original_variable_values_begin( ) const
   { return Variable_value_iterator
       (original_variables_numerator_begin(),
-       boost::bind
-       (boost::bind
-        (Quotient_normalizer(), boost::bind
-         (U_Quotient_creator(), _1, _2)), _1, variables_common_denominator()));
+       [this](const ET& n)
+       {
+         return Quotient_normalizer()(
+           U_Quotient_creator()(n, this->variables_common_denominator()));
+       });
   }
 
   Variable_value_iterator
   original_variable_values_end  ( ) const
   { return Variable_value_iterator
       (original_variables_numerator_end(),
-       boost::bind
-       (boost::bind
-        (Quotient_normalizer(), boost::bind
-         (U_Quotient_creator(), _1, _2)), _1, variables_common_denominator()));
+       [this](const ET& n)
+       {
+         return Quotient_normalizer()(
+           U_Quotient_creator()(n, this->variables_common_denominator()));
+       });
   }
 
   // Basic variables and constraints
@@ -234,10 +231,11 @@ public:
   {
     return Lambda_iterator
      (lambda_numerator_begin(),
-      boost::bind
-      (boost::bind
-       (Quotient_normalizer(), boost::bind
-        (U_Quotient_creator(), _1, _2)), _1, variables_common_denominator()));
+       [this](const ET& n)
+       {
+         return Quotient_normalizer()(
+           U_Quotient_creator()(n, this->variables_common_denominator()));
+       });
   }
 
   Lambda_iterator
@@ -245,10 +243,11 @@ public:
   {
     return Lambda_iterator
      (lambda_numerator_end(),
-      boost::bind
-      (boost::bind
-       (Quotient_normalizer(), boost::bind
-        (U_Quotient_creator(), _1, _2)), _1, variables_common_denominator()));
+       [this](const ET& n)
+       {
+         return Quotient_normalizer()(
+           U_Quotient_creator()(n, this->variables_common_denominator()));
+       });
   }
 
   // destruction
diff --git a/QP_solver/include/CGAL/QP_solver/QP_basis_inverse.h b/QP_solver/include/CGAL/QP_solver/QP_basis_inverse.h
index 91aaa9972f8..b9f3906da7e 100644
--- a/QP_solver/include/CGAL/QP_solver/QP_basis_inverse.h
+++ b/QP_solver/include/CGAL/QP_solver/QP_basis_inverse.h
@@ -410,7 +410,7 @@ class QP_basis_inverse {
         m_it1 = M.begin()+l;
         for ( row = 0; row < s; ++row, ++m_it1) {
             std::transform( m_it1->begin(), m_it1->begin()+s, m_it1->begin(),
-                            boost::bind2nd( std::multiplies<ET>(), d));
+                            [this](const ET& v){return v * this->d;});
         }
 
         // new denominator: |det(A_B)|^2
diff --git a/QP_solver/include/CGAL/QP_solver/QP_basis_inverse_impl.h b/QP_solver/include/CGAL/QP_solver/QP_basis_inverse_impl.h
index 10873b0db55..04e8f3f36d4 100644
--- a/QP_solver/include/CGAL/QP_solver/QP_basis_inverse_impl.h
+++ b/QP_solver/include/CGAL/QP_solver/QP_basis_inverse_impl.h
@@ -178,13 +178,10 @@ z_replace_original_by_original(ForwardIterator y_l_it,
 
     // tmp_l -part
     std::transform(y_l_it, (y_l_it+s), x_l.begin(), tmp_l.begin(),
-        compose2_2(std::plus<ET>(), Identity<ET>(),
-        boost::bind1st(std::multiplies<ET>(), s_delta)));
-
+        [&s_delta](const ET& v1, const ET& v2){ return std::plus<ET>()(v1, s_delta * v2); });
     // tmp_x -part
     std::transform(y_x_it, (y_x_it+b), x_x.begin(), tmp_x.begin(),
-        compose2_2(std::plus<ET>(), Identity<ET>(),
-        boost::bind1st(std::multiplies<ET>(), s_delta)));
+        [&s_delta](const ET& v1, const ET& v2){ return std::plus<ET>()(v1, s_delta * v2); });
     tmp_x[k_i] -= d;
 
     // prepare \hat{k}_{2} -scalar
diff --git a/QP_solver/include/CGAL/QP_solver/QP_solver.h b/QP_solver/include/CGAL/QP_solver/QP_solver.h
index 74c7bb367e1..889fdf25a1b 100644
--- a/QP_solver/include/CGAL/QP_solver/QP_solver.h
+++ b/QP_solver/include/CGAL/QP_solver/QP_solver.h
@@ -38,9 +38,6 @@
 
 #include <CGAL/IO/Verbose_ostream.h>
 
-#include <boost/bind.hpp>
-#include <boost/function.hpp>
-
 #include <CGAL/boost/iterator/transform_iterator.hpp>
 
 #include <vector>
@@ -262,7 +259,7 @@ private:
 
   typedef  QP_matrix_accessor< A_iterator, false, true, false, false>
   A_accessor;
-  typedef  boost::function1<typename A_accessor::result_type, int>
+  typedef  std::function<typename A_accessor::result_type(int)>
   A_row_by_index_accessor;
   typedef  boost::transform_iterator
   < A_row_by_index_accessor, Index_iterator >
@@ -1500,9 +1497,12 @@ transition( Tag_false)
   inv_M_B.transition
     (boost::make_transform_iterator
      (B_O.begin(),
-      boost::bind
-      (D_transition_creator_iterator(), B_O.begin(),
-       boost::bind (D_transition_creator_accessor(), qp_D, _1))));
+      [this](int i)
+      {
+        return D_transition_creator_iterator()(
+          this->B_O.begin(),D_transition_creator_accessor()(this->qp_D, i));
+      })
+    );
 }
 
 template < typename Q, typename ET, typename Tags >  inline                                 // LP case
@@ -1597,12 +1597,10 @@ ratio_test_1__q_x_S( Tag_false)
                     A_by_index_iterator( S_B.begin(),
                                          A_by_index_accessor( *(qp_A + j))),
                     q_x_S.begin(),
-                    boost::bind(std::minus<ET>(),
-                                _1,
-                                boost::bind(std::multiplies<ET>(), d,
-                                  boost::bind(
-                                    NT_converter<RT,ET>(),
-                                    _2))));
+                    [this](const ET& n1, const RT& n2)
+                    {
+                      return std::minus<ET>()(n1,this->d * NT_converter<RT,ET>() (n2));
+                    });
   }
 
   // q_x_S = -+ ( A_S_BxB_O * q_x_O - A_S_Bxj)
@@ -1883,7 +1881,8 @@ basis_matrix_stays_regular()
   if ( has_ineq && (i >= qp_n)) {        // slack variable
     new_row = slack_A[ i-qp_n].first;
     A_row_by_index_accessor  a_accessor =
-      boost::bind (A_accessor( qp_A, 0, qp_n), _1, new_row);
+      [new_row, this](int i){ return A_accessor( this->qp_A, 0, this->qp_n)(i, new_row); };
+
     typedef typename std::iterator_traits<A_row_by_index_iterator>::value_type RT;
     std::transform(A_row_by_index_iterator( B_O.begin(), a_accessor),
                    A_row_by_index_iterator( B_O.end  (), a_accessor),
@@ -1940,17 +1939,14 @@ compute__x_B_S( Tag_false /*has_equalities_only_and_full_rank*/,
                   B_by_index_iterator( S_B.end  (), b_accessor),
                   x_B_S.begin(),
                   x_B_S.begin(),
-                  compose2_2( std::minus<ET>(),
-                              boost::bind1st( std::multiplies<ET>(), d),
-                              Identity<ET>()));
+                  [this](const ET& n1, const ET& n2)
+                  { return std::minus<ET>()(this->d * n1, n2); });
 
   // b_S_B - ( A_S_BxB_O * x_B_O) - r_S_B
   std::transform(x_B_S.begin(), x_B_S.begin()+S_B.size(),
                  r_S_B.begin(), x_B_S.begin(),
-                 compose2_2(std::minus<ET>(),
-                            Identity<ET>(),
-                            boost::bind1st( std::multiplies<ET>(), d)));
-
+                 [this](const ET& n1, const ET& n2)
+                 { return std::minus<ET>()(n1, this->d * n2); });
 
   // x_B_S = +- ( b_S_B - A_S_BxB_O * x_B_O)
   Value_iterator  x_it = x_B_S.begin();
diff --git a/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h b/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h
index bc3788745b2..2d3a66015ef 100644
--- a/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h
+++ b/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h
@@ -14,7 +14,6 @@
 //                 Kaspar Fischer
 
 #include <CGAL/QP_solver/Initialization.h>
-#include <boost/bind.hpp>
 #include <CGAL/NT_converter.h>
 
 namespace CGAL {
@@ -64,9 +63,7 @@ transition( )
     std::transform( C_by_index_iterator( B_O.begin(), c_accessor),
                     C_by_index_iterator( B_O.end  (), c_accessor),
                     minus_c_B.begin(),
-                    boost::bind(
-                      NT_converter<RT,ET>(),
-                      boost::bind(std::negate<RT>(), _1)));
+                    [](const RT& n){return NT_converter<RT,ET>()(-n);} );
 
     // compute initial solution of phase II
     compute_solution(Is_nonnegative());
@@ -1438,7 +1435,7 @@ replace_variable_slack_slack( )
 
     // update basis inverse
     A_row_by_index_accessor  a_accessor =
-      boost::bind( A_accessor( qp_A, 0, qp_n), _1, new_row);
+      [new_row, this](int i){ return A_accessor( this->qp_A, 0, this->qp_n)(i, new_row); };
     typedef typename std::iterator_traits<A_row_by_index_iterator>::value_type RT;
     std::transform(A_row_by_index_iterator( B_O.begin(), a_accessor),
                    A_row_by_index_iterator( B_O.end  (), a_accessor),
@@ -1591,7 +1588,7 @@ replace_variable_original_slack( )
 
     // update basis inverse
     A_row_by_index_accessor  a_accessor =
-      boost::bind (A_accessor( qp_A, 0, qp_n), _1, new_row);
+      [new_row, this](int i){ return A_accessor( this->qp_A, 0, this->qp_n)(i, new_row); };
     typedef typename std::iterator_traits<A_row_by_index_iterator>::value_type RT;
     std::transform(A_row_by_index_iterator( B_O.begin(), a_accessor),
                    A_row_by_index_iterator( B_O.end  (), a_accessor),
@@ -1946,7 +1943,7 @@ leave_variable( )
 
         // update basis inverse
         A_row_by_index_accessor  a_accessor =
-          boost::bind (A_accessor( qp_A, 0, qp_n), _1, new_row);
+          [new_row, this](int i){ return A_accessor( this->qp_A, 0, this->qp_n)(i, new_row); };
         std::copy( A_row_by_index_iterator( B_O.begin(), a_accessor),
                    A_row_by_index_iterator( B_O.end  (), a_accessor),
                    tmp_x.begin());
@@ -2283,7 +2280,7 @@ z_replace_variable_slack_by_original( )
 
     // prepare u
     A_row_by_index_accessor  a_accessor =
-      boost::bind (A_accessor( qp_A, 0, qp_n), _1, new_row);
+      [new_row, this](int i){ return A_accessor( this->qp_A, 0, this->qp_n)(i, new_row); };
     std::copy( A_row_by_index_iterator( B_O.begin(), a_accessor),
                A_row_by_index_iterator( B_O.end  (), a_accessor),
                tmp_x.begin());
@@ -2370,7 +2367,7 @@ z_replace_variable_slack_by_slack( )
     // update basis inverse
     // --------------------
     A_row_by_index_accessor  a_accessor =
-      boost::bind ( A_accessor( qp_A, 0, qp_n), _1, new_row);
+      [new_row, this](int i){ return A_accessor( this->qp_A, 0, this->qp_n)(i, new_row); };
     std::copy( A_row_by_index_iterator( B_O.begin(), a_accessor),
                A_row_by_index_iterator( B_O.end  (), a_accessor),
                tmp_x.begin());
@@ -2937,7 +2934,7 @@ check_basis_inverse( Tag_false)
             }
         }
         v_it = std::find_if( q_x_O.begin(), q_x_O.begin()+cols,
-                             boost::bind2nd( std::not_equal_to<ET>(), et0));
+                             [this](const ET& v){ return v != this->et0; });
         if ( v_it != q_x_O.begin()+cols) {
             if ( ! vout4.verbose()) {
                 std::cerr << std::endl << "basis-inverse check: ";
@@ -2978,7 +2975,7 @@ check_basis_inverse( Tag_false)
         }
 
         v_it = std::find_if( q_lambda.begin(), q_lambda.begin()+rows,
-                             boost::bind2nd( std::not_equal_to<ET>(), et0));
+                             [this](const ET& v){ return v != this->et0; });
         if ( v_it != q_lambda.begin()+rows) {
             if ( ! vout4.verbose()) {
                 std::cerr << std::endl << "basis-inverse check: ";
diff --git a/QP_solver/test/QP_solver/test_bind.cpp b/QP_solver/test/QP_solver/test_bind.cpp
index 3733f359200..ce828b73068 100644
--- a/QP_solver/test/QP_solver/test_bind.cpp
+++ b/QP_solver/test/QP_solver/test_bind.cpp
@@ -2,7 +2,7 @@
 
 #include <CGAL/Quotient.h>
 #include <CGAL/function_objects.h>
-#include <boost/bind.hpp>
+#include <functional>
 
 // functor int x int -> Quotient<int>, (a,b) -> a/b
 // ------------------------------------------------
@@ -20,14 +20,16 @@ struct Quotient_inverter
   }
 };
 
+using namespace std::placeholders;
+
 int main()
 {
   // create composed functor (a,b) -> b/a...
   // ---------------------------------------
   int three = 3;
   int two = 2;
-  std::cout << boost::bind
-    (Quotient_inverter(), boost::bind
+  std::cout << std::bind
+    (Quotient_inverter(), std::bind
      (Quotient_creator(), _1, _2))
   // ...and apply it to (3, 2)
   // -------------------------
diff --git a/STL_Extension/examples/STL_Extension/min_element_if_example.cpp b/STL_Extension/examples/STL_Extension/min_element_if_example.cpp
index a2a98015533..32e43c1abfb 100644
--- a/STL_Extension/examples/STL_Extension/min_element_if_example.cpp
+++ b/STL_Extension/examples/STL_Extension/min_element_if_example.cpp
@@ -14,8 +14,7 @@ int main()
   std::cout << "min_odd = "
             << *CGAL::min_element_if(v.begin(),
                                      v.end(),
-                                     CGAL::compose1_1(boost::bind2nd(std::greater< int >(), 0),
-                                                      boost::bind2nd(std::modulus< int >(), 2)))
+                                     [](int i){ return (i%2) > 0; })
             << std::endl;
   return 0;
 }
diff --git a/STL_Extension/include/CGAL/Spatial_lock_grid_3.h b/STL_Extension/include/CGAL/Spatial_lock_grid_3.h
index 1189a120a33..f5d2b6df798 100644
--- a/STL_Extension/include/CGAL/Spatial_lock_grid_3.h
+++ b/STL_Extension/include/CGAL/Spatial_lock_grid_3.h
@@ -16,8 +16,6 @@
 
 #include <CGAL/Bbox_3.h>
 
-#include <boost/bind.hpp>
-
 #include <atomic>
 #if TBB_IMPLEMENT_CPP0X
 # include <tbb/compat/thread>
@@ -334,7 +332,7 @@ protected:
   Spatial_lock_grid_base_3(const Bbox_3 &bbox,
                                           int num_grid_cells_per_axis)
     : m_num_grid_cells_per_axis(num_grid_cells_per_axis),
-      m_tls_grids(boost::bind(init_TLS_grid, num_grid_cells_per_axis))
+      m_tls_grids([num_grid_cells_per_axis](){ return init_TLS_grid(num_grid_cells_per_axis); })
   {
     set_bbox(bbox);
   }
diff --git a/STL_Extension/include/CGAL/function_objects.h b/STL_Extension/include/CGAL/function_objects.h
index 09fdf6343a6..82709ab1265 100644
--- a/STL_Extension/include/CGAL/function_objects.h
+++ b/STL_Extension/include/CGAL/function_objects.h
@@ -349,19 +349,16 @@ class Creator_uniform_d {
 
 template < class Op1, class Op2 >
 class Unary_compose_1
-  : public CGAL::cpp98::unary_function< typename Op2::argument_type,
-                                        typename Op1::result_type >
 {
 protected:
   Op1 op1;
   Op2 op2;
 public:
-  typedef typename Op2::argument_type   argument_type;
-  typedef typename Op1::result_type     result_type;
 
   Unary_compose_1(const Op1& x, const Op2& y) : op1(x), op2(y) {}
 
-  result_type
+  template <typename argument_type>
+  auto
   operator()(const argument_type& x) const
   { return op1(op2(x)); }
 };
@@ -373,21 +370,18 @@ compose1_1(const Op1& op1, const Op2& op2)
 
 template < class Op1, class Op2, class Op3 >
 class Binary_compose_1
-  : public CGAL::cpp98::unary_function< typename Op2::argument_type,
-                                        typename Op1::result_type >
 {
 protected:
   Op1 op1;
   Op2 op2;
   Op3 op3;
 public:
-  typedef typename Op2::argument_type  argument_type;
-  typedef typename Op1::result_type    result_type;
 
   Binary_compose_1(const Op1& x, const Op2& y, const Op3& z)
   : op1(x), op2(y), op3(z) {}
 
-  result_type
+  template <typename argument_type>
+  auto
   operator()(const argument_type& x) const
   { return op1(op2(x), op3(x)); }
 };
@@ -399,21 +393,16 @@ compose2_1(const Op1& op1, const Op2& op2, const Op3& op3)
 
 template < class Op1, class Op2 >
 class Unary_compose_2
-  : public CGAL::cpp98::binary_function< typename Op2::first_argument_type,
-                                         typename Op2::second_argument_type,
-                                         typename Op1::result_type >
 {
 protected:
   Op1 op1;
   Op2 op2;
 public:
-  typedef typename Op2::first_argument_type   first_argument_type;
-  typedef typename Op2::second_argument_type  second_argument_type;
-  typedef typename Op1::result_type           result_type;
 
   Unary_compose_2(const Op1& x, const Op2& y) : op1(x), op2(y) {}
 
-  result_type
+  template <typename first_argument_type, typename second_argument_type>
+  auto
   operator()(const first_argument_type& x,
              const second_argument_type& y) const
   { return op1(op2(x, y)); }
@@ -426,18 +415,13 @@ compose1_2(const Op1& op1, const Op2& op2)
 
 template < class Op1, class Op2, class Op3 >
 class Binary_compose_2
-  : public CGAL::cpp98::binary_function< typename Op2::argument_type,
-                                         typename Op3::argument_type,
-                                         typename Op1::result_type >
+
 {
 protected:
   Op1 op1;
   Op2 op2;
   Op3 op3;
 public:
-  typedef typename Op2::argument_type  first_argument_type;
-  typedef typename Op3::argument_type  second_argument_type;
-  typedef typename Op1::result_type    result_type;
 
   Binary_compose_2(const Op1& x, const Op2& y, const Op3& z)
   : op1(x), op2(y), op3(z) {}
@@ -455,7 +439,8 @@ public:
     return *this;
   }
 
-  result_type
+  template <typename first_argument_type, typename second_argument_type>
+  auto
   operator()(const first_argument_type& x,
              const second_argument_type& y) const
   { return op1(op2(x), op3(y)); }
diff --git a/STL_Extension/test/STL_Extension/test_composition.cpp b/STL_Extension/test/STL_Extension/test_composition.cpp
index c7deeaa56f0..a0631437393 100644
--- a/STL_Extension/test/STL_Extension/test_composition.cpp
+++ b/STL_Extension/test/STL_Extension/test_composition.cpp
@@ -12,8 +12,6 @@ using CGAL::compose1_2;
 using CGAL::compose2_1;
 using CGAL::compose2_2;
 using CGAL::compare_to_less;
-using boost::binder1st;
-using boost::bind1st;
 using std::accumulate;
 using std::plus;
 using std::multiplies;
@@ -30,9 +28,8 @@ struct Myf {
 int main()
 {
   plus< int >        pl;
-  multiplies< int >  mu;
-  binder1st< plus< int > >       op1 = bind1st(pl, 1);
-  binder1st< multiplies< int > > op2 = bind1st(mu, 2);
+  std::function<int(int)> op1 = [](int i){ return i+1; };
+  std::function<int(int)> op2 = [](int i){ return i * 2; };
 
   // compose1_2:
   int a[] = {3,5,7,2,4};
@@ -50,7 +47,7 @@ int main()
 
   // compose2_1:
   transform(b, b + 5, a, compose2_1(pl, op2, op2));
-  transform(b, b + 5, b, bind1st(mu, 4));
+  transform(b, b + 5, b, [](int i) { return 4 * i;} );
   assert(equal(a, a + 5, b));
 
   // compare_to_less
diff --git a/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_builder_2_impl.h b/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_builder_2_impl.h
index e2936def38a..5adacb71a4f 100644
--- a/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_builder_2_impl.h
+++ b/Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_builder_2_impl.h
@@ -17,7 +17,6 @@
 #include <CGAL/Real_timer.h>
 #include <CGAL/Unique_hash_map.h>
 
-#include <boost/bind.hpp>
 #include <boost/utility.hpp>
 #include <boost/version.hpp>
 #if BOOST_VERSION == 106000
@@ -2018,12 +2017,12 @@ bool Straight_skeleton_builder_2<Gt,Ss,V>::FinishUp()
 
   std::for_each( mSplitNodes.begin()
                 ,mSplitNodes.end  ()
-                ,boost::bind(&Straight_skeleton_builder_2<Gt,Ss,V>::MergeSplitNodes,this,_1)
+                ,[this](Vertex_handle_pair p){ this->MergeSplitNodes(p); }
                ) ;
 
   std::for_each( mDanglingBisectors.begin()
                 ,mDanglingBisectors.end  ()
-                ,boost::bind(&Straight_skeleton_builder_2<Gt,Ss,V>::EraseBisector,this,_1)
+                ,[this](Halfedge_handle db){ this->EraseBisector(db); }
                ) ;
 
   // MergeCoincidentNodes() locks all extremities of halfedges that have a vertex involved in a multinode.
diff --git a/Surface_mesh/examples/Surface_mesh/sm_do_intersect.cpp b/Surface_mesh/examples/Surface_mesh/sm_do_intersect.cpp
index 1f9f8cc4d95..33ac1ef7dad 100644
--- a/Surface_mesh/examples/Surface_mesh/sm_do_intersect.cpp
+++ b/Surface_mesh/examples/Surface_mesh/sm_do_intersect.cpp
@@ -1,3 +1,7 @@
+#include <algorithm>
+#include <vector>
+#include <fstream>
+
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <CGAL/Surface_mesh.h>
 
@@ -89,16 +93,12 @@ unsigned int intersect(const Mesh& P, const Mesh& Q) {
   Q_box_ptr.reserve(Q.number_of_faces());
 
   // build boxes and pointers to boxes
-  boost::transform(P.faces(),
-                 std::back_inserter(P_boxes),
-                 boost::bind(boost::value_factory<Box>(), _1, boost::cref(P)));
-
-
+  for(auto f : P.faces())
+    P_boxes.push_back( Box(f, P) );
   std::transform(P_boxes.begin(), P_boxes.end(), std::back_inserter(P_box_ptr),
                  &address_of_box);
-  boost::transform(Q.faces(),
-                 std::back_inserter(Q_boxes),
-                 boost::bind(boost::value_factory<Box>(), _1, boost::cref(Q)));
+  for(auto f : Q.faces())
+    Q_boxes.push_back( Box(f, Q) );
   std::transform(Q_boxes.begin(), Q_boxes.end(), std::back_inserter(Q_box_ptr),
                  &address_of_box);
 
@@ -136,5 +136,3 @@ int main(int argc, char* argv[])
 
   return 0;
 }
-
-
diff --git a/Surface_mesh/test/Surface_mesh/surface_mesh_test.cpp b/Surface_mesh/test/Surface_mesh/surface_mesh_test.cpp
index ee0e71417de..949d2538bff 100644
--- a/Surface_mesh/test/Surface_mesh/surface_mesh_test.cpp
+++ b/Surface_mesh/test/Surface_mesh/surface_mesh_test.cpp
@@ -4,9 +4,6 @@
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/Simple_cartesian.h>
 
-#include <boost/bind.hpp>
-#include <boost/range/algorithm.hpp>
-
 #include <CGAL/use.h>
 
 void constructors_test()
@@ -118,15 +115,13 @@ void memory_reuse_test()
   // remove all faces
   std::size_t old_face_size = f.m.number_of_faces();
   std::size_t old_removed_face_size = f.m.number_of_removed_faces();
-  boost::range::for_each(f.m.faces(), boost::bind(&Sm::remove_face, boost::ref(f.m), _1));
+  for(auto face : f.m.faces()) f.m.remove_face(face);
   assert(f.m.number_of_faces()== 0);
   assert(f.m.number_of_removed_faces()== old_face_size + old_removed_face_size);
   // remove all edges
   std::size_t old_edge_size = f.m.number_of_edges();
   std::size_t old_removed_edge_size = f.m.number_of_removed_edges();
-  boost::range::for_each(f.m.edges(),
-                         boost::bind(static_cast<void (Sm::*)(Sm::Edge_index)>(&Sm::remove_edge),
-                                     boost::ref(f.m), _1));
+  for(auto e : f.m.edges()) f.m.remove_edge(e);
   assert(f.m.number_of_faces() == 0);
   assert(f.m.number_of_removed_edges()== old_edge_size + old_removed_edge_size);
 
@@ -151,7 +146,7 @@ void memory_reuse_test()
   std::size_t old_size = f.m.number_of_vertices();
   std::size_t old_removed_size = f.m.number_of_removed_vertices();
 
-  boost::range::for_each(f.m.vertices(), boost::bind(&Sm::remove_vertex, boost::ref(f.m), _1));
+  for(auto v : f.m.vertices()) f.m.remove_vertex(v);
   assert(f.m.number_of_vertices() == 0);
   assert(f.m.number_of_removed_vertices()== old_size + old_removed_size);
 
diff --git a/Surface_mesher/include/CGAL/Surface_mesher/Intersection_data_structure_3.h b/Surface_mesher/include/CGAL/Surface_mesher/Intersection_data_structure_3.h
index 52f158a2918..453d4db1212 100644
--- a/Surface_mesher/include/CGAL/Surface_mesher/Intersection_data_structure_3.h
+++ b/Surface_mesher/include/CGAL/Surface_mesher/Intersection_data_structure_3.h
@@ -20,7 +20,6 @@
 #include <CGAL/Bbox_3.h>
 
 #include <boost/format.hpp>
-#include <boost/bind.hpp>
 #include <boost/ref.hpp>
 #include <CGAL/boost/iterator/transform_iterator.hpp>
 
@@ -144,7 +143,6 @@ public:
 
   void create_data_structure()
   {
-    using boost::bind;
     using boost::make_transform_iterator;
 
     max_width = CGAL_NTS max BOOST_PREVENT_MACRO_SUBSTITUTION
diff --git a/Triangulation_3/include/CGAL/Triangulation_3.h b/Triangulation_3/include/CGAL/Triangulation_3.h
index 4e4c5606545..df69bf1f1b6 100644
--- a/Triangulation_3/include/CGAL/Triangulation_3.h
+++ b/Triangulation_3/include/CGAL/Triangulation_3.h
@@ -53,7 +53,6 @@
 #include <CGAL/Bbox_3.h>
 #include <CGAL/Spatial_lock_grid_3.h>
 
-#include <boost/bind.hpp>
 #include <boost/random/linear_congruential.hpp>
 #include <boost/random/uniform_smallint.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -7294,22 +7293,10 @@ operator==(const Triangulation_3<GT, Tds1, Lds>& t1,
     std::vector<Point> V2 (t2.points_begin(), t2.points_end());
 
     std::sort(V1.begin(), V1.end(),
-              boost::bind<Comparison_result>(
-                cmp1,
-                boost::bind<
-                typename boost::result_of<const Construct_point_3(const Point&)>::type>(cp, _1),
-                boost::bind<
-                typename boost::result_of<const Construct_point_3(const Point&)>::type>(cp, _2))
-              == SMALLER);
+              [&cmp1, &cp](const Point& p1, const Point& p2){ return cmp1(cp(p1), cp(p2))==SMALLER; });
 
     std::sort(V2.begin(), V2.end(),
-              boost::bind<Comparison_result>(
-                cmp2,
-                boost::bind<
-                typename boost::result_of<const Construct_point_3(const Point&)>::type>(cp, _1),
-                boost::bind<
-                typename boost::result_of<const Construct_point_3(const Point&)>::type>(cp, _2))
-              == SMALLER);
+              [&cmp2, &cp](const Point& p1, const Point& p2){ return cmp2(cp(p1), cp(p2))==SMALLER; });
 
     return V1 == V2;
   }