diff --git a/include/boost/gil/color_base_algorithm.hpp b/include/boost/gil/color_base_algorithm.hpp
index 6746f83002..1999a73295 100644
--- a/include/boost/gil/color_base_algorithm.hpp
+++ b/include/boost/gil/color_base_algorithm.hpp
@@ -363,6 +363,13 @@ struct element_recursion
         op2(semantic_at_c<N-1>(p1), semantic_at_c<N-1>(p2), semantic_at_c<N-1>(p3));
         return op2;
     }
+    //static_for_each with 5 sources
+    template <typename P1,typename P2,typename P3,typename P4,typename P5,typename Op>
+    static Op static_for_each(const P1& p1, const P2& p2, const P3& p3, const P4& p4, P5& p5, Op op) {
+        Op op2(element_recursion<N-1>::static_for_each(p1,p2,p3,p4,p5,op));
+        op2(semantic_at_c<N-1>(p1), semantic_at_c<N-1>(p2), semantic_at_c<N-1>(p3), semantic_at_c<N-1>(p4), semantic_at_c<N-1>(p5));
+        return op2;
+    }
     //static_transform with one source
     template <typename P1,typename Dst,typename Op>
     static Op static_transform(P1& src, Dst& dst, Op op) {
@@ -401,6 +408,8 @@ struct element_recursion
         semantic_at_c<N-1>(dst)=op2(semantic_at_c<N-1>(src1), semantic_at_c<N-1>(src2));
         return op2;
     }
+
+    
 };
 
 // Termination condition of the compile-time recursion for element operations on a color base
@@ -426,12 +435,16 @@ template<> struct element_recursion<0> {
     //static_for_each with three sources
     template <typename P1,typename P2,typename P3,typename Op>
     static Op static_for_each(const P1&,const P2&,const P3&,Op op){return op;}
+    //static_for_each with 5 sources
+    template <typename P1,typename P2,typename P3,typename P4,typename P5,typename Op>
+    static Op static_for_each(const P1&,const P2&,const P3&,const P4&,const P5&,Op op){return op;}
     //static_transform with one source
     template <typename P1,typename Dst,typename Op>
     static Op static_transform(const P1&,const Dst&,Op op){return op;}
     //static_transform with two sources
     template <typename P1,typename P2,typename Dst,typename Op>
     static Op static_transform(const P1&,const P2&,const Dst&,Op op){return op;}
+
 };
 
 // std::min and std::max don't have the mutable overloads...
@@ -706,6 +719,9 @@ Op static_for_each(const P1& p1,const P2& p2,P3& p3,Op op) { return detail::elem
 template <typename P1,typename P2,typename P3,typename Op>
 BOOST_FORCEINLINE
 Op static_for_each(const P1& p1,const P2& p2,const P3& p3,Op op) { return detail::element_recursion<size<P1>::value>::static_for_each(p1,p2,p3,op); }
+template <typename P1,typename P2,typename P3,typename P4,typename P5,typename Op>
+BOOST_FORCEINLINE
+Op static_for_each(const P1& p1,const P2& p2,const P3& p3,const P4& p4,P5& p5,Op op) { return detail::element_recursion<size<P1>::value>::static_for_each(p1,p2,p3,p4,p5,op); }
 ///\}
 
 } }  // namespace boost::gil
diff --git a/include/boost/gil/extension/numeric/sampler.hpp b/include/boost/gil/extension/numeric/sampler.hpp
index 0f71a986d7..c0560d3f57 100644
--- a/include/boost/gil/extension/numeric/sampler.hpp
+++ b/include/boost/gil/extension/numeric/sampler.hpp
@@ -13,7 +13,7 @@
 
 namespace boost { namespace gil {
 
-// Nearest-neighbor and bilinear image samplers.
+// Nearest-neighbor, bilinear and bicubic image samplers.
 // NOTE: The code is for example use only. It is not optimized for performance
 
 ///////////////////////////////////////////////////////////////////////////
@@ -55,11 +55,18 @@ struct cast_channel_fn {
     }
 };
 
+
+
 template <typename SrcPixel, typename DstPixel>
 void cast_pixel(const SrcPixel& src, DstPixel& dst) {
     static_for_each(src,dst,cast_channel_fn());
 }
 
+// template <typename SrcPixel, typename DstPixel>
+// void cast_pixel_bicubic(const SrcPixel& src, DstPixel& dst) {
+//     static_for_each(src, dst, cubic_interpolate_check());
+// }
+
 namespace detail {
 
 template <typename Weight>
@@ -85,6 +92,56 @@ struct add_dst_mul_src {
 //            dst);
     }
 };
+
+struct cubic_interpolate_check {
+    template <typename SrcChannel, typename DstChannel>
+    void operator()(const SrcChannel& , DstChannel& dst) {
+        // using dst_value_t = typename channel_traits<DstChannel>::value_type;
+        using src_value_t = typename channel_traits<SrcChannel>::value_type;
+
+        if(dst > channel_traits<SrcChannel>::max_value()){
+            dst = src_value_t(channel_traits<SrcChannel>::max_value());
+        }
+        else if(dst < channel_traits<SrcChannel>::min_value()){
+            dst = src_value_t(channel_traits<SrcChannel>::min_value());
+        }
+        else{
+            dst = src_value_t(dst);
+        }
+    }
+};
+
+template <typename Weight>
+struct cubic_interpolate_channel {
+    Weight _w1;
+    Weight _w2;
+    Weight _w3;
+    Weight _w4;
+    cubic_interpolate_channel(Weight w1, Weight w2, Weight w3, Weight w4) : _w1(w1), _w2(w2), _w3(w3), _w4(w4) {}
+
+    template <typename SrcChannel, typename DstChannel>
+    void operator()(const SrcChannel& src1, const SrcChannel& src2, const SrcChannel& src3, const SrcChannel& src4, DstChannel& dst) const {
+        dst = DstChannel(src1*_w1 + src2*_w2 + src3*_w3 + src4*_w4);
+    }
+};
+
+template <typename SrcView, typename SrcP, typename Weight, typename DstP>
+struct cubic_interpolate
+{
+    void operator()(const SrcP& p0, const SrcP& p1, const SrcP& p2, const SrcP& p3, const Weight x, DstP& dst) const 
+    {
+        Weight w0 = ((-0.5 * x) + (x*x) - (0.5 * x*x*x));
+        Weight w1 = (1 - (2.5 * x*x) + (1.5 * x*x*x));
+        Weight w2 = ((0.5 * x) + (2 * x*x) - (1.5 * x*x*x));
+        Weight w3 = ((-0.5 * x*x) + (0.5 * x*x*x));
+
+        static_for_each(p0, p1, p2, p3, dst, cubic_interpolate_channel<Weight>(w0, w1, w2, w3));
+        
+        // Overflow check. Should optimally be after all cubic calculations
+        static_for_each(p0, dst, cubic_interpolate_check());
+    }
+};
+
 } // namespace detail
 
 /// \brief A sampler that sets the destination pixel as the bilinear interpolation of the four closest pixels from the source.
@@ -183,6 +240,262 @@ bool sample(bilinear_sampler, SrcView const& src, point<F> const& p, DstP& resul
 	return true;
 }
 
+/// \brief A sampler that sets the destination pixel as the bicubic interpolation of the 16 closest pixels from the source.
+/// If outside the bounds, it doesn't change the destination. Boundary pixels are the cubic interpolation of 4 adjacent
+/// pixels.
+/// \ingroup ImageAlgorithms
+struct bicubic_sampler {};
+
+template <typename DstP, typename SrcView, typename F>
+bool sample(bicubic_sampler, SrcView const& src, point<F> const& p, DstP& result)
+{
+    using SrcP = typename SrcView::value_type;
+    using cubic_interpolate_current = typename detail::cubic_interpolate<SrcView, SrcP, F, pixel<F, devicen_layout_t<num_channels<SrcView>::value> > >;
+
+    point_t p0(ifloor(p.x), ifloor(p.y)); // the closest integer coordinate top left from p
+    point<F> frac(p.x-p0.x, p.y-p0.y);
+
+    if (p0.x < -1 || p0.y < -1 || p0.x>=src.width() || p0.y>=src.height())
+    {
+        return false;
+    }
+
+	pixel<F,devicen_layout_t<num_channels<SrcView>::value> > mp0(0), mp1(0), mp2(0), mp3(0), mpres(0); // suboptimal
+    
+	typename SrcView::xy_locator loc=src.xy_at(p0.x,p0.y);
+    
+    if(p0.y == -1)
+    {
+        // Topmost row pixels
+        ++loc.y();
+        if(p0.x == -1)
+        {
+            // Top left corner
+            detail::add_dst_mul_src<SrcP,F,pixel<F,devicen_layout_t<num_channels<SrcView>::value> > >()(loc.x()[1], 1, mpres); 
+        }
+        else if(p0.x == 0)
+        {
+            // Top left after corner pixel
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres);
+        }
+        else if(p0.x+1 == src.width())
+        {
+            // Top right pixel
+            detail::add_dst_mul_src<SrcP,F,pixel<F,devicen_layout_t<num_channels<SrcView>::value> > >()(*loc, 1, mpres); 
+        }
+        else if(p0.x+2 == src.width())
+        {
+            // Top right pixel before corner pixel
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mpres);
+        }
+        else
+        {
+            // All other top row pixels
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres);
+        }
+    }
+    else if(p0.y == 0)
+    {
+        // Row after topmost row, 2nd row
+        if(p0.x == -1)
+        {
+            // Left corner pixel of 2nd row
+            ++loc.x();
+            cubic_interpolate_current()(*loc, *loc, *(loc.y()+1), *(loc.y()+2), frac.y, mpres);
+        }
+        else if(p0.x == 0)
+        {
+            // Pixel after left corner pixel of 2nd row
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1);
+            mp0 = mp1;
+            ++loc.y();
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2);
+            ++loc.y();
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3);
+            
+            cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres);
+        }
+        else if(p0.x+1 == src.width())
+        {
+            // Right corner pixel of 2nd row
+            cubic_interpolate_current()(*loc, *loc, *(loc.y()+1), *(loc.y()+2), frac.y, mpres);
+        }
+        else if(p0.x+2 == src.width())
+        {
+            // Pixel before right corner pixel of 2nd row
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1);
+            mp0 = mp1;
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp3);
+            
+            cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres);
+        }
+        else
+        {
+            // All other pixels of 2nd row
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1);
+            mp0 = mp1;
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3);
+            
+            cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres);
+        }
+        
+    }
+    else if(p0.y+1 == src.height())
+    {
+        // Bottom most row
+        if(p0.x == -1)
+        {
+            // Left corner pixel of last row
+            detail::add_dst_mul_src<SrcP,F,pixel<F,devicen_layout_t<num_channels<SrcView>::value> > >()(loc.x()[1], 1, mpres); 
+        }
+        else if(p0.x == 0)
+        {
+            // Pixel after left corner of last row
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres);
+        }
+        else if(p0.x+1 == src.width())
+        {
+            // Right corner of last row
+            detail::add_dst_mul_src<SrcP,F,pixel<F,devicen_layout_t<num_channels<SrcView>::value> > >()(*loc, 1, mpres); 
+        }
+        else if(p0.x+2 == src.width())
+        {
+            // Pixel before right corner of last row
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mpres);
+        }
+        else
+        {
+            // All other pixels of last row
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres);
+        }
+    }
+    else if(p0.y+2 == src.height())
+    {
+        // Row before bottom most row
+        if(p0.x == -1)
+        {
+            // Leftmost corner of 2nd last row
+            ++loc.x();
+            cubic_interpolate_current()(*(loc.y()-1), *loc, *(loc.y()+1), *(loc.y()+1), frac.y, mpres);
+        }
+        else if(p0.x == 0)
+        {
+            // Pixel after left corner of second last row
+            --loc.y();
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0);
+            ++loc.y();
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1);
+            ++loc.y();
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2);
+            mp3 = mp2;
+            
+            cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres);
+        }
+        else if(p0.x+1 == src.width())
+        {
+            // Right corner of second last row
+            cubic_interpolate_current()(*(loc.y()-1), *loc, *(loc.y()+1), *(loc.y()+1), frac.y, mpres);
+        }
+        else if(p0.x+2 == src.width())
+        {
+            // Pixel before right corner of second last row
+            --loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp0);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2);
+            mp3 = mp2;
+            
+            cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres);
+        }
+        else
+        {
+            // All other pixels in second last row
+            --loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2);
+            mp3 = mp2;
+
+            cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres);
+        }
+    }
+    else
+    {
+        // All other rows
+        if(p0.x == -1)
+        {
+            // First column
+            ++loc.x();
+            cubic_interpolate_current()(*(loc.y()-1), *loc, *(loc.y()+1), *(loc.y()+2), frac.y, mpres);
+        }
+        else if(p0.x == 0)
+        {
+            // 2nd column
+            --loc.y();
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0);
+            ++loc.y();
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1);
+            ++loc.y();
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2);
+            ++loc.y();
+            cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3);
+            
+            cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres);
+        }
+        else if(p0.x+1 == src.width())
+        {
+            // Last column
+            cubic_interpolate_current()(*(loc.y()-1), *loc, *(loc.y()+1), *(loc.y()+2), frac.y, mpres);
+        }
+        else if(p0.x+2 == src.width())
+        {
+            // 2nd last column
+            --loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp0);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp3);
+            
+            cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres);
+        }
+        else
+        {
+            // All general cases, center pixels
+            --loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2);
+            ++loc.y();
+            cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3);
+            
+            cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres);
+
+        }
+    }
+
+    SrcP src_result;
+    
+	cast_pixel(mpres, src_result);
+
+	color_convert(src_result, result);
+	return true;
+}
+
 }}  // namespace boost::gil
 
 #endif
diff --git a/test/extension/numeric/resample.cpp b/test/extension/numeric/resample.cpp
index fa86162535..dd17a52853 100644
--- a/test/extension/numeric/resample.cpp
+++ b/test/extension/numeric/resample.cpp
@@ -95,9 +95,61 @@ void test_bilinear_sampler_test()
     BOOST_TEST_EQ(gil::rgb8_pixel_t(0, 128, 0), dv(3, 3));
 }
 
+void test_bicubic_sampler_test()
+{
+    // R G B R
+    // G W W G
+    // B W W B
+    // R G B R
+    gil::rgb8_image_t img(4, 4);
+    gil::rgb8_view_t v = view(img);
+    v(0, 0) = v(0, 3) = v(3, 0) = v(3, 3) = gil::rgb8_pixel_t(255, 0, 0);
+    v(0, 1) = v(1, 0) = v(1, 3) = v(3, 1) = gil::rgb8_pixel_t(0, 255, 0);
+    v(0, 2) = v(2, 0) = v(2, 3) = v(3, 2) =  gil::rgb8_pixel_t(0, 0, 255);
+    v(1, 1) = v(1, 2) = v(2, 1) = v(2, 2) = gil::rgb8_pixel_t(255, 255, 255);
+
+    gil::rgb8_image_t dims(5, 5);
+    gil::rgb8c_view_t dv = gil::const_view(dims);
+
+    test_map_fn<double, gil::rgb8_image_t::coord_t> mf;
+
+    gil::resample_pixels(gil::const_view(img), gil::view(dims), mf, gil::bicubic_sampler());
+
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(255, 0, 0), dv(0, 0));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 143, 0), dv(0, 1));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(0, 143, 143), dv(0, 2));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 0, 143), dv(0, 3));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(255, 0, 0), dv(0, 4));
+
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 143, 0), dv(1, 0));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 207, 55), dv(1, 1));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 197, 214), dv(1, 2));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 127, 135), dv(1, 3));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 143, 0), dv(1, 4));
+
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(0, 143, 143), dv(2, 0));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 199, 199), dv(2, 1));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(255, 255, 255), dv(2, 2));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 199, 199), dv(2, 3));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(0, 143, 143), dv(2, 4));
+
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 0, 143), dv(3, 0));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 135, 127), dv(3, 1));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 214, 197), dv(3, 2));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 55, 207), dv(3, 3));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 0, 143), dv(3, 4));
+
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(255, 0, 0), dv(4, 0));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 143, 0), dv(4, 1));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(0, 143, 143), dv(4, 2));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(127, 0, 143), dv(4, 3));
+    BOOST_TEST_EQ(gil::rgb8_pixel_t(255, 0, 0), dv(4, 4));
+}
+
 int main()
 {
     test_bilinear_sampler_test();
+    test_bicubic_sampler_test();
 
     return ::boost::report_errors();
 }