From 793cd334fde88f89e8899be63b15ab0433518546 Mon Sep 17 00:00:00 2001 From: Scramjet911 Date: Sat, 13 Mar 2021 09:05:33 +0530 Subject: [PATCH 1/3] Bicubic sampler with bugs --- .../boost/gil/extension/numeric/sampler.hpp | 340 ++++++++++++++++++ 1 file changed, 340 insertions(+) diff --git a/include/boost/gil/extension/numeric/sampler.hpp b/include/boost/gil/extension/numeric/sampler.hpp index 0f71a986d7..95fd41d84f 100644 --- a/include/boost/gil/extension/numeric/sampler.hpp +++ b/include/boost/gil/extension/numeric/sampler.hpp @@ -85,6 +85,24 @@ struct add_dst_mul_src { // dst); } }; + + +template +struct bicubic_interpolate +{ + void operator()(const SrcP& p0, const SrcP& p1, const SrcP& p2, const SrcP& p3, const Weight x, DstP& dst) const + { + using add_dst_mul_src_current = add_dst_mul_src::value> > >; + + add_dst_mul_src_current()(p0, ((-0.5 * x) + (x*x) - (0.5 * x*x*x)), dst); + add_dst_mul_src_current()(p1, (1 - (2.5 * x*x) + (1.5 * x*x*x)), dst); + add_dst_mul_src_current()(p2, ((0.5 * x) + (2 * x*x) - (1.5 * x*x*x)), dst); + add_dst_mul_src_current()(p3, ((-0.5 * x*x) + (0.5 * x*x*x)), dst); + } +}; + + + } // namespace detail /// \brief A sampler that sets the destination pixel as the bilinear interpolation of the four closest pixels from the source. @@ -183,6 +201,328 @@ bool sample(bilinear_sampler, SrcView const& src, point const& p, DstP& resul return true; } +/// \brief A sampler that sets the destination pixel as the bicubic interpolation of the four closest pixels from the source. +/// If outside the bounds, it doesn't change the destination +/// \ingroup ImageAlgorithms +struct bicubic_sampler {}; + +template +bool sample(bicubic_sampler, SrcView const& src, point const& p, DstP& result) +{ + using SrcP = typename SrcView::value_type; + using bicubic_interpolate_current = typename detail::bicubic_interpolate::value> > >; + + point_t p0(ifloor(p.x), ifloor(p.y)); // the closest integer coordinate top left from p + point 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::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::value> > >()(loc.x()[1], 1, mpres); + } + else if(p0.x == 0) + { + // Top left after corner pixel + bicubic_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::value> > >()(*loc, 1, mpres); + } + else if(p0.x+2 == src.width()) + { + // Top right pixel before corner pixel + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mpres); + } + else + { + // All other top row pixels + bicubic_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(); + bicubic_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 + bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + mp0 = mp1; + ++loc.y(); + bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + ++loc.y(); + bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); + + bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + } + else if(p0.x+1 == src.width()) + { + // Right corner pixel of 2nd row + bicubic_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 + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1); + mp0 = mp1; + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp3); + + bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + } + else + { + // All other pixels of 2nd row + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + mp0 = mp1; + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); + + bicubic_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::value> > >()(loc.x()[1], 1, mpres); + } + else if(p0.x == 0) + { + // Pixel after left corner of last row + bicubic_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::value> > >()(*loc, 1, mpres); + } + else if(p0.x+2 == src.width()) + { + // Pixel before right corner of last row + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mpres); + } + else + { + // All other pixels of last row + bicubic_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(); + bicubic_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(); + bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); + ++loc.y(); + bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + ++loc.y(); + bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + mp3 = mp2; + + bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + } + else if(p0.x+1 == src.width()) + { + // Right corner of second last row + bicubic_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(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp0); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2); + mp3 = mp2; + + bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + } + else + { + // All other pixels in second last row + --loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + mp3 = mp2; + + bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + } + } + else + { + // All other rows + if(p0.x == -1) + { + // First column + ++loc.x(); + bicubic_interpolate_current()(*(loc.y()-1), *loc, *(loc.y()+1), *(loc.y()+2), frac.y, mpres); + } + else if(p0.x == 0) + { + // 2nd column + --loc.y(); + bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); + ++loc.y(); + bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + ++loc.y(); + bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + ++loc.y(); + bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); + + bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + } + else if(p0.x+1 == src.width()) + { + // Last column + bicubic_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(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp0); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp3); + + bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + } + else + { + // All general cases, center pixels + --loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + ++loc.y(); + bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); + + bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + } + } + + // if(p0.y == -1) + // { + // ++loc.y(); + // detail::add_dst_mul_src::value> > >()(*loc, 1, mpres); + // } + // else if(p0.y == 0 || p0.y+1 == src.height() || p0.y+2 == src.height()) + // { + // detail::add_dst_mul_src::value> > >()(*loc, 1, mpres); + // } + // else if (p0.y+2 < src.height()) + // { + // if(p0.x == -1) + // { + // ++loc.x(); + // detail::add_dst_mul_src::value> > >()(*loc, 1, mpres); + // } + // if(p0.x+1 == src.width()) + // { + // detail::add_dst_mul_src::value> > >()(*loc, 1, mpres); + // } + // // if(p0.x == -1 || p0.x+1 == src.width()) + // // { + // // // First or last column pixels + // // bicubic_interpolate_current()(*(loc.y()-1), *loc, *(loc.y()+1), *(loc.y()+2), frac.y, mpres); + // // } + // else if(p0.x == 0) + // { + // // Most common case, inside image not any corner or last rows or columns + // --loc.y(); + // bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); + // ++loc.y(); + // bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + // ++loc.y(); + // bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + // ++loc.y(); + // bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); + + // bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + // } + // else if(p0.x+2 == src.width()) + // { + // // + // --loc.y(); + // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp0); + // ++loc.y(); + // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1); + // ++loc.y(); + // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2); + // ++loc.y(); + // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp3); + + // bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + // } + // else if(p0.x+2 < src.width()) + // { + // // Most common case, inside image not any corner or last rows or columns + // --loc.y(); + // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); + // ++loc.y(); + // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + // ++loc.y(); + // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + // ++loc.y(); + // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); + + // bicubic_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 From 7979ec6686c841138185a4a066042ae19759c9f8 Mon Sep 17 00:00:00 2001 From: Scramjet911 Date: Sat, 13 Mar 2021 10:51:58 +0530 Subject: [PATCH 2/3] Cleaned bicubic sampler with bugs --- .../boost/gil/extension/numeric/sampler.hpp | 69 ------------------- 1 file changed, 69 deletions(-) diff --git a/include/boost/gil/extension/numeric/sampler.hpp b/include/boost/gil/extension/numeric/sampler.hpp index 95fd41d84f..c228bb98d8 100644 --- a/include/boost/gil/extension/numeric/sampler.hpp +++ b/include/boost/gil/extension/numeric/sampler.hpp @@ -446,75 +446,6 @@ bool sample(bicubic_sampler, SrcView const& src, point const& p, DstP& result bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); } } - - // if(p0.y == -1) - // { - // ++loc.y(); - // detail::add_dst_mul_src::value> > >()(*loc, 1, mpres); - // } - // else if(p0.y == 0 || p0.y+1 == src.height() || p0.y+2 == src.height()) - // { - // detail::add_dst_mul_src::value> > >()(*loc, 1, mpres); - // } - // else if (p0.y+2 < src.height()) - // { - // if(p0.x == -1) - // { - // ++loc.x(); - // detail::add_dst_mul_src::value> > >()(*loc, 1, mpres); - // } - // if(p0.x+1 == src.width()) - // { - // detail::add_dst_mul_src::value> > >()(*loc, 1, mpres); - // } - // // if(p0.x == -1 || p0.x+1 == src.width()) - // // { - // // // First or last column pixels - // // bicubic_interpolate_current()(*(loc.y()-1), *loc, *(loc.y()+1), *(loc.y()+2), frac.y, mpres); - // // } - // else if(p0.x == 0) - // { - // // Most common case, inside image not any corner or last rows or columns - // --loc.y(); - // bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); - // ++loc.y(); - // bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); - // ++loc.y(); - // bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); - // ++loc.y(); - // bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); - - // bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); - // } - // else if(p0.x+2 == src.width()) - // { - // // - // --loc.y(); - // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp0); - // ++loc.y(); - // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1); - // ++loc.y(); - // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2); - // ++loc.y(); - // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp3); - - // bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); - // } - // else if(p0.x+2 < src.width()) - // { - // // Most common case, inside image not any corner or last rows or columns - // --loc.y(); - // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); - // ++loc.y(); - // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); - // ++loc.y(); - // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); - // ++loc.y(); - // bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); - - // bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); - // } - // } SrcP src_result; cast_pixel(mpres, src_result); From b81d3cc75c4b73c6b5349e646d934efff560f109 Mon Sep 17 00:00:00 2001 From: Scramjet911 Date: Tue, 30 Mar 2021 14:18:14 +0530 Subject: [PATCH 3/3] Bicubic Sampler and static_for_each with 5 sources --- include/boost/gil/color_base_algorithm.hpp | 16 ++ .../boost/gil/extension/numeric/sampler.hpp | 170 +++++++++++------- test/extension/numeric/resample.cpp | 52 ++++++ 3 files changed, 174 insertions(+), 64 deletions(-) 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(p1), semantic_at_c(p2), semantic_at_c(p3)); return op2; } + //static_for_each with 5 sources + template + 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::static_for_each(p1,p2,p3,p4,p5,op)); + op2(semantic_at_c(p1), semantic_at_c(p2), semantic_at_c(p3), semantic_at_c(p4), semantic_at_c(p5)); + return op2; + } //static_transform with one source template static Op static_transform(P1& src, Dst& dst, Op op) { @@ -401,6 +408,8 @@ struct element_recursion semantic_at_c(dst)=op2(semantic_at_c(src1), semantic_at_c(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 static Op static_for_each(const P1&,const P2&,const P3&,Op op){return op;} + //static_for_each with 5 sources + template + 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 static Op static_transform(const P1&,const Dst&,Op op){return op;} //static_transform with two sources template 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 BOOST_FORCEINLINE Op static_for_each(const P1& p1,const P2& p2,const P3& p3,Op op) { return detail::element_recursion::value>::static_for_each(p1,p2,p3,op); } +template +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::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 c228bb98d8..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 void cast_pixel(const SrcPixel& src, DstPixel& dst) { static_for_each(src,dst,cast_channel_fn()); } +// template +// void cast_pixel_bicubic(const SrcPixel& src, DstPixel& dst) { +// static_for_each(src, dst, cubic_interpolate_check()); +// } + namespace detail { template @@ -86,23 +93,55 @@ struct add_dst_mul_src { } }; +struct cubic_interpolate_check { + template + void operator()(const SrcChannel& , DstChannel& dst) { + // using dst_value_t = typename channel_traits::value_type; + using src_value_t = typename channel_traits::value_type; + + if(dst > channel_traits::max_value()){ + dst = src_value_t(channel_traits::max_value()); + } + else if(dst < channel_traits::min_value()){ + dst = src_value_t(channel_traits::min_value()); + } + else{ + dst = src_value_t(dst); + } + } +}; + +template +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 + 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 -struct bicubic_interpolate +struct cubic_interpolate { void operator()(const SrcP& p0, const SrcP& p1, const SrcP& p2, const SrcP& p3, const Weight x, DstP& dst) const { - using add_dst_mul_src_current = add_dst_mul_src::value> > >; + 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)); - add_dst_mul_src_current()(p0, ((-0.5 * x) + (x*x) - (0.5 * x*x*x)), dst); - add_dst_mul_src_current()(p1, (1 - (2.5 * x*x) + (1.5 * x*x*x)), dst); - add_dst_mul_src_current()(p2, ((0.5 * x) + (2 * x*x) - (1.5 * x*x*x)), dst); - add_dst_mul_src_current()(p3, ((-0.5 * x*x) + (0.5 * x*x*x)), dst); + static_for_each(p0, p1, p2, p3, dst, cubic_interpolate_channel(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. @@ -201,8 +240,9 @@ bool sample(bilinear_sampler, SrcView const& src, point const& p, DstP& resul return true; } -/// \brief A sampler that sets the destination pixel as the bicubic interpolation of the four closest pixels from the source. -/// If outside the bounds, it doesn't change the destination +/// \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 {}; @@ -210,7 +250,7 @@ template bool sample(bicubic_sampler, SrcView const& src, point const& p, DstP& result) { using SrcP = typename SrcView::value_type; - using bicubic_interpolate_current = typename detail::bicubic_interpolate::value> > >; + using cubic_interpolate_current = typename detail::cubic_interpolate::value> > >; point_t p0(ifloor(p.x), ifloor(p.y)); // the closest integer coordinate top left from p point frac(p.x-p0.x, p.y-p0.y); @@ -223,7 +263,7 @@ bool sample(bicubic_sampler, SrcView const& src, point const& p, DstP& result pixel::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 @@ -236,7 +276,7 @@ bool sample(bicubic_sampler, SrcView const& src, point const& p, DstP& result else if(p0.x == 0) { // Top left after corner pixel - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres); } else if(p0.x+1 == src.width()) { @@ -246,12 +286,12 @@ bool sample(bicubic_sampler, SrcView const& src, point const& p, DstP& result else if(p0.x+2 == src.width()) { // Top right pixel before corner pixel - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mpres); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mpres); } else { // All other top row pixels - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres); } } else if(p0.y == 0) @@ -261,48 +301,48 @@ bool sample(bicubic_sampler, SrcView const& src, point const& p, DstP& result { // Left corner pixel of 2nd row ++loc.x(); - bicubic_interpolate_current()(*loc, *loc, *(loc.y()+1), *(loc.y()+2), frac.y, mpres); + 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 - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); mp0 = mp1; ++loc.y(); - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); ++loc.y(); - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); - bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); } else if(p0.x+1 == src.width()) { // Right corner pixel of 2nd row - bicubic_interpolate_current()(*loc, *loc, *(loc.y()+1), *(loc.y()+2), frac.y, mpres); + 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 - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1); mp0 = mp1; ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp3); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp3); - bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); } else { // All other pixels of 2nd row - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); mp0 = mp1; ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); - bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); } } @@ -317,7 +357,7 @@ bool sample(bicubic_sampler, SrcView const& src, point const& p, DstP& result else if(p0.x == 0) { // Pixel after left corner of last row - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres); } else if(p0.x+1 == src.width()) { @@ -327,12 +367,12 @@ bool sample(bicubic_sampler, SrcView const& src, point const& p, DstP& result else if(p0.x+2 == src.width()) { // Pixel before right corner of last row - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mpres); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mpres); } else { // All other pixels of last row - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mpres); } } else if(p0.y+2 == src.height()) @@ -342,51 +382,51 @@ bool sample(bicubic_sampler, SrcView const& src, point const& p, DstP& result { // Leftmost corner of 2nd last row ++loc.x(); - bicubic_interpolate_current()(*(loc.y()-1), *loc, *(loc.y()+1), *(loc.y()+1), frac.y, mpres); + 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(); - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); ++loc.y(); - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); ++loc.y(); - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); mp3 = mp2; - bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); } else if(p0.x+1 == src.width()) { // Right corner of second last row - bicubic_interpolate_current()(*(loc.y()-1), *loc, *(loc.y()+1), *(loc.y()+1), frac.y, mpres); + 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(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp0); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp0); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2); mp3 = mp2; - bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); } else { // All other pixels in second last row --loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); mp3 = mp2; - bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); } } else @@ -396,58 +436,60 @@ bool sample(bicubic_sampler, SrcView const& src, point const& p, DstP& result { // First column ++loc.x(); - bicubic_interpolate_current()(*(loc.y()-1), *loc, *(loc.y()+1), *(loc.y()+2), frac.y, mpres); + 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(); - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); ++loc.y(); - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); ++loc.y(); - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); ++loc.y(); - bicubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); + cubic_interpolate_current()(*loc, *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); - bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); } else if(p0.x+1 == src.width()) { // Last column - bicubic_interpolate_current()(*(loc.y()-1), *loc, *(loc.y()+1), *(loc.y()+2), frac.y, mpres); + 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(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp0); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp0); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp1); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp2); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp3); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+1), frac.x, mp3); - bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); } else { // All general cases, center pixels --loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp0); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp1); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp2); ++loc.y(); - bicubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); + cubic_interpolate_current()(*(loc.x()-1), *loc, *(loc.x()+1), *(loc.x()+2), frac.x, mp3); - bicubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + cubic_interpolate_current()(mp0, mp1, mp2, mp3, frac.y, mpres); + } } SrcP src_result; + cast_pixel(mpres, src_result); color_convert(src_result, result); 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 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(); }