diff --git a/docs/filter.cpp b/docs/filter.cpp new file mode 100644 index 0000000..3cc19f7 --- /dev/null +++ b/docs/filter.cpp @@ -0,0 +1,846 @@ +/* + Copyright 2008 Larry Gritz and the other authors and contributors. + All Rights Reserved. + Based on BSD-licensed software Copyright 2004 NVIDIA Corp. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are + met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the software's owners nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + (This is the Modified BSD License) +*/ + + +// Filter results comparison +// pow2 downres (4K -> 128 tested) +// - katana 2.7.12 +// - ImageMagick 6.3.6 +// - prman 16.0 txmake +// +// box: oiio, prman, katana, imagemagick match +// lanczos3: oiio, katana, imagemagick match. prman is far sharper (perhaps lanczos2?) +// sinc: oiio, prman match. Katana is slighly softer. imagemagick is much softer +// blackman harris: all differ. In order of decreasing sharpness... imagemagick, oiio, prman +// catrom: oiio, imagemagick, prman match +// gaussian: prman, katana match (gaussian3). imgmagick, oiio are sharper (gaussian2) + + + +#include +#include +#include +#include + +#include +#include +#include + + +OIIO_NAMESPACE_BEGIN + +// Below are the implementations of several 2D filters. They all +// inherit their interface from Filter2D. Each must redefine two +// virtual functions: +// char *name() Return the filter name +// float operator(float,float) Evaluate the filter +// + +class FilterBox1D : public Filter1D { +public: + FilterBox1D (float width) : Filter1D(width) { } + ~FilterBox1D (void) { } + float operator() (float x) const { + return (fabsf(x) <= m_w*0.5f) ? 1.0f : 0.0f; + } + string_view name (void) const { return "box"; } +}; + + + +class FilterBox2D : public Filter2D { +public: + FilterBox2D (float width, float height) : Filter2D(width,height) { } + ~FilterBox2D (void) { } + float operator() (float x, float y) const { + return (fabsf(x) <= m_w*0.5f && fabsf(y) <= m_h*0.5f) ? 1.0f : 0.0f; + } + bool separable (void) const { return true; } + float xfilt (float x) const { return fabsf(x) <= m_w*0.5f ? 1.0f : 0.0f; } + float yfilt (float y) const { return fabsf(y) <= m_h*0.5f ? 1.0f : 0.0f; } + string_view name (void) const { return "box"; } +}; + + + +class FilterTriangle1D : public Filter1D { +public: + FilterTriangle1D (float width) : Filter1D(width), m_rad_inv(2.0f/width) { } + ~FilterTriangle1D (void) { } + float operator() (float x) const { + return tri1d (x * m_rad_inv); + } + string_view name (void) const { return "triangle"; } + + static float tri1d (float x) { + x = fabsf(x); + return (x < 1.0f) ? (1.0f - x) : 0.0f; + } +private: + float m_rad_inv; +}; + + + +class FilterTriangle2D : public Filter2D { +public: + FilterTriangle2D (float width, float height) + : Filter2D(width,height), m_wrad_inv(2.0f/width), + m_hrad_inv(2.0f/height) { } + ~FilterTriangle2D (void) { } + float operator() (float x, float y) const { + return FilterTriangle1D::tri1d (x * m_wrad_inv) + * FilterTriangle1D::tri1d (y * m_hrad_inv); + } + bool separable (void) const { return true; } + float xfilt (float x) const { + return FilterTriangle1D::tri1d (x * m_wrad_inv); + } + float yfilt (float y) const { + return FilterTriangle1D::tri1d (y * m_hrad_inv); + } + string_view name (void) const { return "triangle"; } +private: + float m_wrad_inv, m_hrad_inv; +}; + + + +class FilterGaussian1D : public Filter1D { +public: + FilterGaussian1D (float width) : Filter1D(width), m_rad_inv(2.0f/width) { } + ~FilterGaussian1D (void) { } + float operator() (float x) const { + return gauss1d (x * m_rad_inv); + } + static float gauss1d (float x) { + x = fabsf(x); + return (x < 1.0f) ? fast_exp (-2.0f * (x*x)) : 0.0f; + } + string_view name (void) const { return "gaussian"; } +private: + float m_rad_inv; +}; + + + +class FilterGaussian2D : public Filter2D { +public: + FilterGaussian2D (float width, float height) + : Filter2D(width,height), m_wrad_inv(2.0f/width), + m_hrad_inv(2.0f/height) { } + ~FilterGaussian2D (void) { } + float operator() (float x, float y) const { + return FilterGaussian1D::gauss1d (x * m_wrad_inv) + * FilterGaussian1D::gauss1d (y * m_hrad_inv); + } + bool separable (void) const { return true; } + float xfilt (float x) const { + return FilterGaussian1D::gauss1d (x * m_wrad_inv); + } + float yfilt (float y) const { + return FilterGaussian1D::gauss1d (y * m_hrad_inv); + } + string_view name (void) const { return "gaussian"; } +private: + float m_wrad_inv, m_hrad_inv; +}; + + + +class FilterSharpGaussian1D : public Filter1D { +public: + FilterSharpGaussian1D (float width) : Filter1D(width), m_rad_inv(2.0f/width) { } + ~FilterSharpGaussian1D (void) { } + float operator() (float x) const { + return gauss1d (x * m_rad_inv); + } + static float gauss1d (float x) { + x = fabsf(x); + return (x < 1.0f) ? fast_exp (-4.0f * (x*x)) : 0.0f; + } + string_view name (void) const { return "gaussian"; } +private: + float m_rad_inv; +}; + + + +class FilterSharpGaussian2D : public Filter2D { +public: + FilterSharpGaussian2D (float width, float height) + : Filter2D(width,height), m_wrad_inv(2.0f/width), + m_hrad_inv(2.0f/height) { } + ~FilterSharpGaussian2D (void) { } + float operator() (float x, float y) const { + return FilterSharpGaussian1D::gauss1d (x * m_wrad_inv) + * FilterSharpGaussian1D::gauss1d (y * m_hrad_inv); + } + bool separable (void) const { return true; } + float xfilt (float x) const { + return FilterSharpGaussian1D::gauss1d (x * m_wrad_inv); + } + float yfilt (float y) const { + return FilterSharpGaussian1D::gauss1d (y * m_hrad_inv); + } + string_view name (void) const { return "gaussian"; } +private: + float m_wrad_inv, m_hrad_inv; +}; + + + +class FilterCatmullRom1D : public Filter1D { +public: + FilterCatmullRom1D (float width) + : Filter1D(4.0f), m_scale(4.0f/width) { } + ~FilterCatmullRom1D (void) { } + float operator() (float x) const { return catrom1d(x * m_scale); } + string_view name (void) const { return "catmull-rom"; } + + static float catrom1d (float x) { + x = fabsf(x); + float x2 = x * x; + float x3 = x * x2; + return (x >= 2.0f) ? 0.0f : ((x < 1.0f) ? + (3.0f * x3 - 5.0f * x2 + 2.0f) : + (-x3 + 5.0f * x2 - 8.0f * x + 4.0f) ); + } +private: + float m_scale; +}; + + + +class FilterCatmullRom2D : public Filter2D { +public: + FilterCatmullRom2D (float width, float height) + : Filter2D(width,height), m_wscale(4.0f/width), + m_hscale(4.0f/height) { } + ~FilterCatmullRom2D (void) { } + float operator() (float x, float y) const { + return FilterCatmullRom1D::catrom1d(x * m_wscale) + * FilterCatmullRom1D::catrom1d(y * m_hscale); + } + bool separable (void) const { return true; } + float xfilt (float x) const { return FilterCatmullRom1D::catrom1d(x * m_wscale); } + float yfilt (float y) const { return FilterCatmullRom1D::catrom1d(y * m_hscale); } + string_view name (void) const { return "catmull-rom"; } +private: + float m_wscale, m_hscale; +}; + + + +class FilterBlackmanHarris1D : public Filter1D { +public: + FilterBlackmanHarris1D (float width) + : Filter1D(width), m_rad_inv(2.0f/width) { } + ~FilterBlackmanHarris1D (void) { } + float operator() (float x) const { + return bh1d (x * m_rad_inv); + } + string_view name (void) const { return "blackman-harris"; } + static float bh1d (float x) { + if (x < -1.0f || x > 1.0f) // Early out if outside filter range + return 0.0f; + // Compute BH. Straight from classic BH paper, but the usual + // formula assumes that the filter is centered at 0.5, so scale: + x = (x + 1.0f) * 0.5f; + const float A0 = 0.35875f; + const float A1 = -0.48829f; + const float A2 = 0.14128f; + const float A3 = -0.01168f; + const float m_pi = float (M_PI); +#if 0 + // original -- three cos calls! + return A0 + A1 * cosf(2.f * m_pi * x) + + A2 * cosf(4.f * m_pi * x) + A3 * cosf(6.f * m_pi * x); +#else + // Use trig identintities to reduce to just one cos. + // https://en.wikipedia.org/wiki/List_of_trigonometric_identities + // cos(2x) = 2 cos^2(x) - 1 + // cos(3x) = 4 cos^3(x) − 3 cos(x) + float cos2pix = cosf(2.f * m_pi * x); + float cos4pix = 2.0f * cos2pix * cos2pix - 1.0f; + float cos6pix = cos2pix * (2.0f * cos4pix - 1.0f); + return A0 + A1 * cos2pix + A2 * cos4pix + A3 * cos6pix; +#endif + } +private: + float m_rad_inv; +}; + + + +class FilterBlackmanHarris2D : public Filter2D { +public: + FilterBlackmanHarris2D (float width, float height) + : Filter2D(width,height), + m_wrad_inv(2.0f/width), m_hrad_inv(2.0f/height) { } + ~FilterBlackmanHarris2D (void) { } + float operator() (float x, float y) const { + return FilterBlackmanHarris1D::bh1d (x*m_wrad_inv) + * FilterBlackmanHarris1D::bh1d (y*m_hrad_inv); + } + bool separable (void) const { return true; } + float xfilt (float x) const { return FilterBlackmanHarris1D::bh1d(x*m_wrad_inv); } + float yfilt (float y) const { return FilterBlackmanHarris1D::bh1d(y*m_hrad_inv); } + string_view name (void) const { return "blackman-harris"; } +private: + float m_wrad_inv, m_hrad_inv; +}; + + + +class FilterSinc1D : public Filter1D { +public: + FilterSinc1D (float width) : Filter1D(width), m_rad(width/2.0f) { } + ~FilterSinc1D (void) { } + float operator() (float x) const { return sinc1d (x, m_rad); } + string_view name (void) const { return "sinc"; } + + static float sinc1d (float x, float rad) { + x = fabsf(x); + if (x > rad) + return 0.0f; + const float m_pi = float (M_PI); + return (x < 0.0001f) ? 1.0f : sinf (m_pi*x)/(m_pi*x); + } +private: + float m_rad; +}; + + + +class FilterSinc2D : public Filter2D { +public: + FilterSinc2D (float width, float height) + : Filter2D(width,height), + m_wrad(width/2.0f), m_hrad(height/2.0f) { } + ~FilterSinc2D (void) { } + float operator() (float x, float y) const { + return FilterSinc1D::sinc1d(x,m_wrad) * FilterSinc1D::sinc1d(y,m_hrad); + } + bool separable (void) const { return true; } + float xfilt (float x) const { return FilterSinc1D::sinc1d(x,m_wrad); } + float yfilt (float y) const { return FilterSinc1D::sinc1d(y,m_hrad); } + string_view name (void) const { return "sinc"; } +private: + float m_wrad, m_hrad; +}; + + + +class FilterLanczos3_1D : public Filter1D { +public: + FilterLanczos3_1D (float width) + : Filter1D(width), m_scale(6.0f/width) { } + ~FilterLanczos3_1D (void) { } + float operator() (float x) const { + return lanczos3 (x * m_scale); + } + string_view name (void) const { return "lanczos3"; } + + static float lanczos3 (float x) { + const float a = 3.0f; // Lanczos 3 lobe + const float ainv = 1.0f / a; + const float m_pi = float (M_PI); + x = fabsf(x); + if (x > a) + return 0.0f; + if (x < 0.0001f) + return 1.0f; +#if 0 + // Full precision, for reference: + float pix = m_pi * x; + return a/(x*x*(m_pi*m_pi)) * sinf(pix)*sinf(pix*ainv); +#elif 0 + // Use approximate fast_sinphi -- about 0.1% absolute error, but + // around 2.5x times faster. BUT when graphed, looks very icky near + // f(0). + return a/(x*x*(m_pi*m_pi)) * fast_sinpi(x)*fast_sinpi(x*ainv); +#else + // Compromise: full-precision sin(), but use the trig identity + // sin(3x) = -4 sin^3(x) + 3 sin(x) + // to make it so only one call to sin is sufficient. This is still + // about 1.5x the speed of the reference implementation, but with + // no precision compromises. + float s1 = sinf(x*ainv*m_pi); // sin(x*pi/a) + float s3 = (-4.0f * s1*s1 + 3.0f) * s1; // sin(3*x*pi/a) == sin(x*pi) + return a/(x*x*(m_pi*m_pi)) * s1 * s3; +#endif + } +private: + float m_scale; +}; + + + +class FilterLanczos3_2D : public Filter2D { +public: + FilterLanczos3_2D (float width, float height) + : Filter2D(width,height), m_wscale(6.0f/width), + m_hscale(6.0f/height) { } + ~FilterLanczos3_2D (void) { } + float operator() (float x, float y) const { + return FilterLanczos3_1D::lanczos3 (x * m_wscale) + * FilterLanczos3_1D::lanczos3 (y * m_hscale); + } + bool separable (void) const { return true; } + float xfilt (float x) const { return FilterLanczos3_1D::lanczos3(x * m_wscale); } + float yfilt (float y) const { return FilterLanczos3_1D::lanczos3(y * m_hscale); } + string_view name (void) const { return "lanczos3"; } +protected: + float m_wscale, m_hscale; +}; + + + +class FilterRadialLanczos3_2D : public FilterLanczos3_2D { +public: + FilterRadialLanczos3_2D (float width, float height) + : FilterLanczos3_2D(width,height) { } + float operator() (float x, float y) const { + x *= m_wscale; + y *= m_hscale; + return FilterLanczos3_1D::lanczos3(sqrtf(x*x + y*y)); + } + bool separable (void) const { return false; } + string_view name (void) const { return "radial-lanczos3"; } +}; + + + +class FilterMitchell1D : public Filter1D { +public: + FilterMitchell1D (float width) : Filter1D(width), m_rad_inv(2.0f/width) { } + ~FilterMitchell1D (void) { } + float operator() (float x) const { + return mitchell1d (x * m_rad_inv); + } + string_view name (void) const { return "mitchell"; } + + static float mitchell1d (float x) { + x = fabsf (x); + if (x > 1.0f) + return 0.0f; + // Computation stright out of the classic Mitchell paper. + // In the paper, the range is -2 to 2, so we rescale: + x *= 2.0f; + float x2 = x*x; + const float B = 1.0f/3.0f; + const float C = 1.0f/3.0f; + const float SIXTH = 1.0f/6.0f; + if (x >= 1.0f) + return ((-B - 6.0f*C)*x*x2 + (6.0f*B + 30.0f*C)*x2 + + (-12.0f*B - 48.0f*C)*x + (8.0f*B + 24.0f*C)) * SIXTH; + else + return ((12.0f - 9.0f*B - 6.0f*C)*x*x2 + + (-18.0f + 12.0f*B + 6.0f*C)*x2 + (6.0f - 2.0f*B)) * SIXTH; + } +private: + float m_rad_inv; +}; + + + +class FilterMitchell2D : public Filter2D { +public: + FilterMitchell2D (float width, float height) + : Filter2D(width,height), + m_wrad_inv(2.0f/width), m_hrad_inv(2.0f/height) { } + ~FilterMitchell2D (void) { } + float operator() (float x, float y) const { + return FilterMitchell1D::mitchell1d (x * m_wrad_inv) + * FilterMitchell1D::mitchell1d (y * m_hrad_inv); + } + bool separable (void) const { return true; } + float xfilt (float x) const { + return FilterMitchell1D::mitchell1d (x * m_wrad_inv); + } + float yfilt (float y) const { + return FilterMitchell1D::mitchell1d (y * m_hrad_inv); + } + string_view name (void) const { return "mitchell"; } +private: + float m_wrad_inv, m_hrad_inv; +}; + + + +// B-spline filter from Stark et al, JGT 10(1) +class FilterBSpline1D : public Filter1D { +public: + FilterBSpline1D (float width) + : Filter1D(width), m_wscale(4.0f/width) + { } + ~FilterBSpline1D (void) { } + float operator() (float x) const { + return bspline1d (x*m_wscale); + } + string_view name (void) const { return "b-spline"; } + + static float bspline1d (float x) { + x = fabsf (x); + if (x <= 1.0f) + return b1 (1.0f-x); + else if (x < 2.0f) + return b0 (2.0f-x); + else return 0.0f; + } +private: + float m_wscale; // width scale factor + static float b0 (float t) { return t*t*t / 6.0f; } + static float b1 (float t) { + return 0.5f * t * (t * (1.0f - t) + 1.0f) + 1.0f/6.0f; + } +}; + + + +class FilterBSpline2D : public Filter2D { +public: + FilterBSpline2D (float width, float height) + : Filter2D(width,height), m_wscale(4.0f/width), m_hscale(4.0f/height) + { } + ~FilterBSpline2D (void) { } + float operator() (float x, float y) const { + return FilterBSpline1D::bspline1d (x * m_wscale) + * FilterBSpline1D::bspline1d (y * m_hscale); + } + bool separable (void) const { return true; } + float xfilt (float x) const { + return FilterBSpline1D::bspline1d(x*m_wscale); + } + float yfilt (float y) const { + return FilterBSpline1D::bspline1d(y*m_hscale); + } + string_view name (void) const { return "b-spline"; } +private: + float m_wscale, m_hscale; +}; + + + +class FilterDisk2D : public Filter2D { +public: + FilterDisk2D (float width, float height) : Filter2D(width,height) { } + ~FilterDisk2D (void) { } + float operator() (float x, float y) const { + x /= (m_w*0.5f); + y /= (m_h*0.5f); + return ((x*x+y*y) < 1.0f) ? 1.0f : 0.0f; + } + string_view name (void) const { return "disk"; } +}; + + +class FilterCubic1D : public Filter1D { +public: + FilterCubic1D (float width) + : Filter1D(width), m_a(0.0f), m_rad_inv(2.0f / width) + { } + ~FilterCubic1D (void) { } + float operator() (float x) const { + return cubic(x * m_rad_inv, m_a); + } + + static float cubic (float x, float a) { + x = fabsf (x); + if (x > 1.0f) + return 0.0f; + // Because range is -2 to 2, we rescale + x *= 2.0f; + + if (x >= 1.0f) + return a * (x * (x * (x - 5.0f) + 8.0f) - 4.0f); + // return a * x*x*x - 5.0f * a * x*x + 8.0f * a * x - 4.0f * a; + else + return x*x*((a + 2.0f) * x - (a + 3.0f)) + 1.0f; + // return (a + 2.0f) * x*x*x - (a + 3.0f) * x*x + 1.0f; + } + + virtual string_view name (void) const { return "cubic"; } +protected: + float m_a; + float m_rad_inv; +}; + + + +class FilterCubic2D : public Filter2D { +public: + FilterCubic2D (float width, float height) + : Filter2D(width,height), m_a(0.0f) + , m_wrad_inv(2.0f/width), m_hrad_inv(2.0f/height) { } + ~FilterCubic2D (void) { } + float operator() (float x, float y) const { + return FilterCubic1D::cubic(x * m_wrad_inv, m_a) + * FilterCubic1D::cubic(y * m_hrad_inv, m_a); + } + bool separable (void) const { return true; } + float xfilt (float x) const { + return FilterCubic1D::cubic (x * m_wrad_inv, m_a); + } + float yfilt (float y) const { + return FilterCubic1D::cubic (y * m_hrad_inv, m_a); + } + virtual string_view name (void) const { return "cubic"; } +protected: + float m_a; + float m_wrad_inv, m_hrad_inv; +}; + + + +class FilterKeys1D : public FilterCubic1D { +public: + FilterKeys1D (float width) : FilterCubic1D(width) { m_a = -0.5f; } + ~FilterKeys1D (void) { } + virtual string_view name (void) const { return "keys"; } +}; + + +class FilterKeys2D : public FilterCubic2D { +public: + FilterKeys2D (float width, float height) : FilterCubic2D(width,height) { m_a = -0.5f; } + ~FilterKeys2D (void) { } + virtual string_view name (void) const { return "keys"; } +}; + + + +class FilterSimon1D : public FilterCubic1D { +public: + FilterSimon1D (float width) : FilterCubic1D(width) { m_a = -0.75f; } + ~FilterSimon1D (void) { } + virtual string_view name (void) const { return "simon"; } +}; + + +class FilterSimon2D : public FilterCubic2D { +public: + FilterSimon2D (float width, float height) : FilterCubic2D(width,height) { m_a = -0.75f; } + ~FilterSimon2D (void) { } + virtual string_view name (void) const { return "simon"; } +}; + + + +class FilterRifman1D : public FilterCubic1D { +public: + FilterRifman1D (float width) : FilterCubic1D(width) { m_a = -1.0f; } + ~FilterRifman1D (void) { } + virtual string_view name (void) const { return "rifman"; } +}; + + +class FilterRifman2D : public FilterCubic2D { +public: + FilterRifman2D (float width, float height) : FilterCubic2D(width,height) { m_a = -1.0f; } + ~FilterRifman2D (void) { } + virtual string_view name (void) const { return "rifman"; } +}; + + + +namespace { +FilterDesc filter1d_list[] = { + // name dim width fixedwidth scalable separable + { "box", 1, 1, false, true, true }, + { "triangle", 1, 2, false, true, true }, + { "gaussian", 1, 3, false, true, true }, + { "sharp-gaussian", 1, 2, false, true, true }, + { "catmull-rom", 1, 4, false, true, true }, + { "blackman-harris", 1, 3, false, true, true }, + { "sinc", 1, 4, false, true, true }, + { "lanczos3", 1, 6, false, true, true }, + { "mitchell", 1, 4, false, true, true }, + { "bspline", 1, 4, false, true, true }, + { "cubic", 1, 4, false, true, true }, + { "keys", 1, 4, false, true, true }, + { "simon", 1, 4, false, true, true }, + { "rifman", 1, 4, false, true, true } +}; +} + +int +Filter1D::num_filters () +{ + return sizeof(filter1d_list)/sizeof(filter1d_list[0]); +} + +void +Filter1D::get_filterdesc (int filternum, FilterDesc *filterdesc) +{ + ASSERT (filternum >= 0 && filternum < num_filters()); + *filterdesc = filter1d_list[filternum]; +} + + + +// Filter1D::create is the static method that, given a filter name, +// width, and height, returns an allocated and instantiated filter of +// the correct implementation. If the name is not recognized, return +// NULL. +Filter1D * +Filter1D::create (string_view filtername, float width) +{ + if (filtername == "box") + return new FilterBox1D (width); + if (filtername == "triangle") + return new FilterTriangle1D (width); + if (filtername == "gaussian") + return new FilterGaussian1D (width); + if (filtername == "sharp-gaussian") + return new FilterSharpGaussian1D (width); + if (filtername == "catmull-rom" || filtername == "catrom") + return new FilterCatmullRom1D (width); + if (filtername == "blackman-harris") + return new FilterBlackmanHarris1D (width); + if (filtername == "sinc") + return new FilterSinc1D (width); + if (filtername == "lanczos3" || filtername == "lanczos") + return new FilterLanczos3_1D (width); + if (filtername == "mitchell") + return new FilterMitchell1D (width); + if (filtername == "b-spline" || filtername == "bspline") + return new FilterBSpline1D (width); + if (filtername == "cubic") + return new FilterCubic1D (width); + if (filtername == "keys") + return new FilterKeys1D (width); + if (filtername == "simon") + return new FilterSimon1D (width); + if (filtername == "rifman") + return new FilterRifman1D (width); + return NULL; +} + + + +void +Filter1D::destroy (Filter1D *filt) +{ + delete filt; +} + + + +static FilterDesc filter2d_list[] = { + // name dim width fixedwidth scalable separable + { "box", 2, 1, false, true, true }, + { "triangle", 2, 2, false, true, true }, + { "gaussian", 2, 3, false, true, true }, + { "sharp-gaussian", 2, 2, false, true, true }, + { "catmull-rom", 2, 4, false, true, true }, + { "blackman-harris", 2, 3, false, true, true }, + { "sinc", 2, 4, false, true, true }, + { "lanczos3", 2, 6, false, true, true }, + { "radial-lanczos3", 2, 6, false, true, false }, + { "mitchell", 2, 4, false, true, true }, + { "bspline", 2, 4, false, true, true }, + { "disk", 2, 1, false, true, false }, + { "cubic", 2, 4, false, true, true }, + { "keys", 2, 4, false, true, true }, + { "simon", 2, 4, false, true, true }, + { "rifman", 2, 4, false, true, true } +}; + +int +Filter2D::num_filters () +{ + return sizeof(filter2d_list)/sizeof(filter2d_list[0]); +} + +void +Filter2D::get_filterdesc (int filternum, FilterDesc *filterdesc) +{ + ASSERT (filternum >= 0 && filternum < num_filters()); + *filterdesc = filter2d_list[filternum]; +} + + + +// Filter2D::create is the static method that, given a filter name, +// width, and height, returns an allocated and instantiated filter of +// the correct implementation. If the name is not recognized, return +// NULL. + +Filter2D * +Filter2D::create (string_view filtername, float width, float height) +{ + if (filtername == "box") + return new FilterBox2D (width, height); + if (filtername == "triangle") + return new FilterTriangle2D (width, height); + if (filtername == "gaussian") + return new FilterGaussian2D (width, height); + if (filtername == "sharp-gaussian") + return new FilterSharpGaussian2D (width, height); + if (filtername == "catmull-rom" || filtername == "catrom") + return new FilterCatmullRom2D (width, height); + if (filtername == "blackman-harris") + return new FilterBlackmanHarris2D (width, height); + if (filtername == "sinc") + return new FilterSinc2D (width, height); + if (filtername == "lanczos3" || filtername == "lanczos") + return new FilterLanczos3_2D (width, height); + if (filtername == "radial-lanczos3" || filtername == "radial-lanczos") + return new FilterRadialLanczos3_2D (width, height); + if (filtername == "mitchell") + return new FilterMitchell2D (width, height); + if (filtername == "b-spline" || filtername == "bspline") + return new FilterBSpline2D (width, height); + if (filtername == "disk") + return new FilterDisk2D (width, height); + if (filtername == "cubic") + return new FilterCubic2D (width, height); + if (filtername == "keys") + return new FilterKeys2D (width, height); + if (filtername == "simon") + return new FilterSimon2D (width, height); + if (filtername == "rifman") + return new FilterRifman2D (width, height); + return NULL; +} + + + +void +Filter2D::destroy (Filter2D *filt) +{ + delete filt; +} + + +OIIO_NAMESPACE_END diff --git a/docs/imagebufalgo_xform.cpp b/docs/imagebufalgo_xform.cpp new file mode 100644 index 0000000..77223a2 --- /dev/null +++ b/docs/imagebufalgo_xform.cpp @@ -0,0 +1,1102 @@ +/* + Copyright 2008 Larry Gritz and the other authors and contributors. + All Rights Reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are + met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the software's owners nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + (This is the Modified BSD License) +*/ + +/// \file +/// ImageBufAlgo functions for filtered transformations + + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include "imageio_pvt.h" + +OIIO_NAMESPACE_BEGIN + + +namespace { + +// Poor man's Dual2 makes it easy to compute with differentials. For +// a rich man's implementation and full documentation, see +// OpenShadingLanguage (dual2.h). +class Dual2 { +public: + float val() const { return m_val; } + float dx() const { return m_dx; } + float dy() const { return m_dy; } + Dual2 (float val) : m_val(val), m_dx(0.0f), m_dy(0.0f) {} + Dual2 (float val, float dx, float dy) : m_val(val), m_dx(dx), m_dy(dy) {} + Dual2& operator= (float f) { m_val = f; m_dx = 0.0f; m_dy = 0.0f; return *this; } + friend Dual2 operator+ (const Dual2 &a, const Dual2 &b) { + return Dual2 (a.m_val+b.m_val, a.m_dx+b.m_dx, a.m_dy+b.m_dy); + } + friend Dual2 operator+ (const Dual2 &a, float b) { + return Dual2 (a.m_val+b, a.m_dx, a.m_dy); + } + friend Dual2 operator* (const Dual2 &a, float b) { + return Dual2 (a.m_val*b, a.m_dx*b, a.m_dy*b); + } + friend Dual2 operator* (const Dual2 &a, const Dual2 &b) { + // Use the chain rule + return Dual2 (a.m_val*b.m_val, + a.m_val*b.m_dx + a.m_dx*b.m_val, + a.m_val*b.m_dy + a.m_dy*b.m_val); + } + friend Dual2 operator/ (const Dual2 &a, const Dual2 &b) { + float bvalinv = 1.0f / b.m_val; + float aval_bval = a.m_val * bvalinv; + return Dual2 (aval_bval, + bvalinv * (a.m_dx - aval_bval * b.m_dx), + bvalinv * (a.m_dy - aval_bval * b.m_dy)); + } +private: + float m_val, m_dx, m_dy; +}; + +/// Transform a 2D point (x,y) with derivatives by a 3x3 affine matrix to +/// obtain a transformed point with derivatives. +inline void +robust_multVecMatrix (const Imath::M33f &M, + const Dual2 &x, const Dual2 &y, + Dual2 &outx, Dual2 &outy) +{ + Dual2 a = x * M[0][0] + y * M[1][0] + M[2][0]; + Dual2 b = x * M[0][1] + y * M[1][1] + M[2][1]; + Dual2 w = x * M[0][2] + y * M[1][2] + M[2][2]; + + if (w.val() != 0.0f) { + outx = a / w; + outy = b / w; + } else { + outx = 0.0f; + outy = 0.0f; + } +} + + + +// Transform an ROI by an affine matrix. +ROI +transform (const Imath::M33f &M, ROI roi) +{ + Imath::V2f ul (roi.xbegin+0.5f, roi.ybegin+0.5f); + Imath::V2f ur (roi.xend-0.5f, roi.ybegin+0.5f); + Imath::V2f ll (roi.xbegin+0.5f, roi.yend-0.5f); + Imath::V2f lr (roi.xend-0.5f, roi.yend-0.5f); + M.multVecMatrix (ul, ul); + M.multVecMatrix (ur, ur); + M.multVecMatrix (ll, ll); + M.multVecMatrix (lr, lr); + Imath::Box2f box (ul); + box.extendBy (ll); + box.extendBy (ur); + box.extendBy (lr); + int xmin = int (floorf(box.min.x)); + int ymin = int (floorf(box.min.y)); + int xmax = int (floorf(box.max.x)) + 1; + int ymax = int (floorf(box.max.y)) + 1; + return ROI (xmin, xmax, ymin, ymax, roi.zbegin, roi.zend, roi.chbegin, roi.chend); +} + + + +// Given s,t image space coordinates and their derivatives, compute a +// filtered sample using the derivatives to guide the size of the filter +// footprint. +template +inline void +filtered_sample (const ImageBuf &src, float s, float t, + float dsdx, float dtdx, float dsdy, float dtdy, + const Filter2D *filter, ImageBuf::WrapMode wrap, + float *result) +{ + DASSERT (filter); + // Just use isotropic filtering + float ds = std::max (1.0f, std::max (fabsf(dsdx), fabsf(dsdy))); + float dt = std::max (1.0f, std::max (fabsf(dtdx), fabsf(dtdy))); + float ds_inv = 1.0f / ds; + float dt_inv = 1.0f / dt; + float filterrad_s = 0.5f * ds * filter->width(); + float filterrad_t = 0.5f * dt * filter->width(); + ImageBuf::ConstIterator samp (src, + (int)floorf(s-filterrad_s), (int)ceilf(s+filterrad_s), + (int)floorf(t-filterrad_t), (int)ceilf(t+filterrad_t), + 0, 1, wrap); + int nc = src.nchannels(); + float *sum = ALLOCA (float, nc); + memset (sum, 0, nc*sizeof(float)); + float total_w = 0.0f; + for ( ; ! samp.done(); ++samp) { + float w = (*filter) (ds_inv*(samp.x()+0.5f-s), dt_inv*(samp.y()+0.5f-t)); + for (int c = 0; c < nc; ++c) + sum[c] += w * samp[c]; + total_w += w; + } + if (total_w != 0.0f) + for (int c = 0; c < nc; ++c) + result[c] = sum[c] / total_w; + else + for (int c = 0; c < nc; ++c) + result[c] = 0.0f; +} + +} // end anon namespace + + + + +template +static bool +resize_ (ImageBuf &dst, const ImageBuf &src, + Filter2D *filter, ROI roi, int nthreads) +{ + ImageBufAlgo::parallel_image (roi, nthreads, [&](ROI roi){ + + const ImageSpec &srcspec (src.spec()); + const ImageSpec &dstspec (dst.spec()); + int nchannels = dstspec.nchannels; + + // Local copies of the source image window, converted to float + float srcfx = srcspec.full_x; + float srcfy = srcspec.full_y; + float srcfw = srcspec.full_width; + float srcfh = srcspec.full_height; + + // Ratios of dst/src size. Values larger than 1 indicate that we + // are maximizing (enlarging the image), and thus want to smoothly + // interpolate. Values less than 1 indicate that we are minimizing + // (shrinking the image), and thus want to properly filter out the + // high frequencies. + float xratio = float(dstspec.full_width) / srcfw; // 2 upsize, 0.5 downsize + float yratio = float(dstspec.full_height) / srcfh; + + float dstfx = float (dstspec.full_x); + float dstfy = float (dstspec.full_y); + float dstfw = float (dstspec.full_width); + float dstfh = float (dstspec.full_height); + float dstpixelwidth = 1.0f / dstfw; + float dstpixelheight = 1.0f / dstfh; + float *pel = ALLOCA (float, nchannels); + float filterrad = filter->width() / 2.0f; + + // radi,radj is the filter radius, as an integer, in source pixels. We + // will filter the source over [x-radi, x+radi] X [y-radj,y+radj]. + int radi = (int) ceilf (filterrad/xratio); + int radj = (int) ceilf (filterrad/yratio); + int xtaps = 2*radi + 1; + int ytaps = 2*radj + 1; + bool separable = filter->separable(); + float *yfiltval = ALLOCA (float, ytaps); + float *xfiltval_all = NULL; + if (separable) { + // For separable filters, horizontal tap weights will be the same + // for every column. So we precompute all the tap weights for every + // x position we'll need. We do the same thing in y, but row by row + // inside the loop (since we never revisit a y row). This + // substantially speeds up resize. + xfiltval_all = ALLOCA (float, xtaps * roi.width()); + for (int x = roi.xbegin; x < roi.xend; ++x) { + float *xfiltval = xfiltval_all + (x-roi.xbegin) * xtaps; + float s = (x-dstfx+0.5f)*dstpixelwidth; + float src_xf = srcfx + s * srcfw; + int src_x; + float src_xf_frac = floorfrac (src_xf, &src_x); + for (int c = 0; c < nchannels; ++c) + pel[c] = 0.0f; + float totalweight_x = 0.0f; + for (int i = 0; i < xtaps; ++i) { + float w = filter->xfilt (xratio * (i-radi-(src_xf_frac-0.5f))); + xfiltval[i] = w; + totalweight_x += w; + } + if (totalweight_x != 0.0f) + for (int i = 0; i < xtaps; ++i) // normalize x filter + xfiltval[i] /= totalweight_x; // weights + } + } + +#if 0 + std::cerr << "Resizing " << srcspec.full_width << "x" << srcspec.full_height + << " to " << dstspec.full_width << "x" << dstspec.full_height << "\n"; + std::cerr << "ratios = " << xratio << ", " << yratio << "\n"; + std::cerr << "examining src filter " << filter->name() + << " support radius of " << radi << " x " << radj << " pixels\n"; + std::cout << " " << xtaps << "x" << ytaps << " filter taps\n"; + std::cerr << "dst range " << roi << "\n"; + std::cerr << "separable filter\n"; +#endif + +#define USE_SPECIAL 0 +#if USE_SPECIAL + // Special case: src and dst are local memory, float buffers, and we're + // operating on all channels, <= 4. + bool special = + ( (is_same::value || is_same::value) + && (is_same::value || is_same::value) + // && dst.localpixels() // has to be, because it's writeable + && src.localpixels() + // && R.contains_roi(roi) // has to be, because IBAPrep + && src.contains_roi(roi) + && roi.chbegin == 0 && roi.chend == dst.nchannels() + && roi.chend == src.nchannels() && roi.chend <= 4 + && separable); +#endif + + // We're going to loop over all output pixels we're interested in. + // + // (s,t) = NDC space coordinates of the output sample we are computing. + // This is the "sample point". + // (src_xf, src_xf) = source pixel space float coordinates of the + // sample we're computing. We want to compute the weighted sum + // of all the source image pixels that fall under the filter when + // centered at that location. + // (src_x, src_y) = image space integer coordinates of the floor, + // i.e., the closest pixel in the source image. + // src_xf_frac and src_yf_frac are the position within that pixel + // of our sample. + // + // Separate cases for separable and non-separable filters. + if (separable) { + ImageBuf::Iterator out (dst, roi); + ImageBuf::ConstIterator srcpel (src, ImageBuf::WrapClamp); + for (int y = roi.ybegin; y < roi.yend; ++y) { + float t = (y-dstfy+0.5f)*dstpixelheight; + float src_yf = srcfy + t * srcfh; + int src_y; + float src_yf_frac = floorfrac (src_yf, &src_y); + // If using separable filters, our vertical set of filter tap + // weights will be the same for the whole scanline we're on. Just + // compute and normalize them once. + float totalweight_y = 0.0f; + for (int j = 0; j < ytaps; ++j) { + float w = filter->yfilt (yratio * (j-radj-(src_yf_frac-0.5f))); + yfiltval[j] = w; + totalweight_y += w; + } + if (totalweight_y != 0.0f) + for (int i = 0; i < ytaps; ++i) + yfiltval[i] /= totalweight_y; + + for (int x = roi.xbegin; x < roi.xend; ++x, ++out) { + float s = (x-dstfx+0.5f)*dstpixelwidth; + float src_xf = srcfx + s * srcfw; + int src_x = ifloor (src_xf); + for (int c = 0; c < nchannels; ++c) + pel[c] = 0.0f; + const float *xfiltval = xfiltval_all + (x-roi.xbegin) * xtaps; + float totalweight_x = 0.0f; + for (int i = 0; i < xtaps; ++i) + totalweight_x += xfiltval[i]; + if (totalweight_x != 0.0f) { + srcpel.rerange (src_x-radi, src_x+radi+1, + src_y-radj, src_y+radj+1, + 0, 1, ImageBuf::WrapClamp); + for (int j = -radj; j <= radj; ++j) { + float wy = yfiltval[j+radj]; + if (wy == 0.0f) { + // 0 weight for this y tap -- move to next line + srcpel.pos (srcpel.x(), srcpel.y()+1, srcpel.z()); + continue; + } + for (int i = 0; i < xtaps; ++i, ++srcpel) { + float w = wy * xfiltval[i]; + if (w) + for (int c = 0; c < nchannels; ++c) + pel[c] += w * srcpel[c]; + } + } + } + // Copy the pixel value (already normalized) to the output. + DASSERT (out.x() == x && out.y() == y); + if (totalweight_y == 0.0f) { + // zero it out + for (int c = 0; c < nchannels; ++c) + out[c] = 0.0f; + } else { + for (int c = 0; c < nchannels; ++c) + out[c] = pel[c]; + } + } + } + + } else { + // Non-separable filter + ImageBuf::Iterator out (dst, roi); + ImageBuf::ConstIterator srcpel (src, ImageBuf::WrapClamp); + for (int y = roi.ybegin; y < roi.yend; ++y) { + float t = (y-dstfy+0.5f)*dstpixelheight; + float src_yf = srcfy + t * srcfh; + int src_y; + float src_yf_frac = floorfrac (src_yf, &src_y); + for (int x = roi.xbegin; x < roi.xend; ++x, ++out) { + float s = (x-dstfx+0.5f)*dstpixelwidth; + float src_xf = srcfx + s * srcfw; + int src_x; + float src_xf_frac = floorfrac (src_xf, &src_x); + for (int c = 0; c < nchannels; ++c) + pel[c] = 0.0f; + float totalweight = 0.0f; + srcpel.rerange (src_x-radi, src_x+radi+1, + src_y-radi, src_y+radi+1, + 0, 1, ImageBuf::WrapClamp); + for (int j = -radj; j <= radj; ++j) { + for (int i = -radi; i <= radi; ++i, ++srcpel) { + DASSERT (! srcpel.done()); + float w = (*filter)(xratio * (i-(src_xf_frac-0.5f)), + yratio * (j-(src_yf_frac-0.5f))); + if (w) { + totalweight += w; + for (int c = 0; c < nchannels; ++c) + pel[c] += w * srcpel[c]; + } + } + } + DASSERT (srcpel.done()); + // Rescale pel to normalize the filter and write it to the + // output image. + DASSERT (out.x() == x && out.y() == y); + if (totalweight == 0.0f) { + // zero it out + for (int c = 0; c < nchannels; ++c) + out[c] = 0.0f; + } else { + for (int c = 0; c < nchannels; ++c) + out[c] = pel[c] / totalweight; + } + } + } + } + + }); // end of parallel_image + return true; +} + + + +static std::shared_ptr +get_resize_filter (string_view filtername, float fwidth, + ImageBuf& dst, float wratio, float hratio) +{ + // Set up a shared pointer with custom deleter to make sure any + // filter we allocate here is properly destroyed. + std::shared_ptr filter ((Filter2D*)nullptr, Filter2D::destroy); + if (filtername.empty()) { + // No filter name supplied -- pick a good default + if (wratio > 1.0f || hratio > 1.0f) + filtername = "blackman-harris"; + else + filtername = "lanczos3"; + } + for (int i = 0, e = Filter2D::num_filters(); i < e; ++i) { + FilterDesc fd; + Filter2D::get_filterdesc (i, &fd); + if (fd.name == filtername) { + float w = fwidth > 0.0f ? fwidth : fd.width * std::max (1.0f, wratio); + float h = fwidth > 0.0f ? fwidth : fd.width * std::max (1.0f, hratio); + filter.reset (Filter2D::create (filtername, w, h)); + break; + } + } + if (! filter) { + dst.error ("Filter \"%s\" not recognized", filtername); + } + return filter; +} + + + +bool +ImageBufAlgo::resize (ImageBuf &dst, const ImageBuf &src, + Filter2D *filter, ROI roi, int nthreads) +{ + pvt::LoggedTimer logtime("IBA::resize"); + if (! IBAprep (roi, &dst, &src, + IBAprep_REQUIRE_SAME_NCHANNELS | IBAprep_NO_SUPPORT_VOLUME | + IBAprep_NO_COPY_ROI_FULL)) + return false; + + // Set up a shared pointer with custom deleter to make sure any + // filter we allocate here is properly destroyed. + std::shared_ptr filterptr ((Filter2D*)NULL, Filter2D::destroy); + if (! filter) { + // If no filter was provided, punt and just linearly interpolate. + const ImageSpec &srcspec (src.spec()); + const ImageSpec &dstspec (dst.spec()); + float wratio = float(dstspec.full_width) / float(srcspec.full_width); + float hratio = float(dstspec.full_height) / float(srcspec.full_height); + float w = 2.0f * std::max (1.0f, wratio); + float h = 2.0f * std::max (1.0f, hratio); + filter = Filter2D::create ("triangle", w, h); + filterptr.reset (filter); + } + + bool ok; + OIIO_DISPATCH_COMMON_TYPES2 (ok, "resize", resize_, + dst.spec().format, src.spec().format, + dst, src, filter, roi, nthreads); + return ok; +} + + + +bool +ImageBufAlgo::resize (ImageBuf &dst, const ImageBuf &src, + string_view filtername, float fwidth, + ROI roi, int nthreads) +{ + pvt::LoggedTimer logtime("IBA::resize"); + if (! IBAprep (roi, &dst, &src, + IBAprep_REQUIRE_SAME_NCHANNELS | IBAprep_NO_SUPPORT_VOLUME | + IBAprep_NO_COPY_ROI_FULL)) + return false; + const ImageSpec &srcspec (src.spec()); + const ImageSpec &dstspec (dst.spec()); + + // Resize ratios + float wratio = float(dstspec.full_width) / float(srcspec.full_width); + float hratio = float(dstspec.full_height) / float(srcspec.full_height); + + // Set up a shared pointer with custom deleter to make sure any + // filter we allocate here is properly destroyed. + auto filter = get_resize_filter (filtername, fwidth, dst, wratio, hratio); + if (! filter) + return false; // error issued in get_resize_filter + + logtime.stop(); // it will be picked up again by the next call... + return resize (dst, src, filter.get(), roi, nthreads); +} + + + +ImageBuf +ImageBufAlgo::resize (const ImageBuf &src, + Filter2D *filter, ROI roi, int nthreads) +{ + ImageBuf result; + bool ok = resize (result, src, filter, roi, nthreads); + if (!ok && !result.has_error()) + result.error ("ImageBufAlgo::resize() error"); + return result; +} + + +ImageBuf +ImageBufAlgo::resize (const ImageBuf &src, + string_view filtername, float filterwidth, + ROI roi, int nthreads) +{ + ImageBuf result; + bool ok = resize (result, src, filtername, filterwidth, roi, nthreads); + if (!ok && !result.has_error()) + result.error ("ImageBufAlgo::resize() error"); + return result; +} + + + +bool +ImageBufAlgo::fit (ImageBuf &dst, const ImageBuf &src, + Filter2D *filter, bool exact, + ROI roi, int nthreads) +{ + // No time logging, it will be accounted in the underlying warp/resize + if (! IBAprep (roi, &dst, &src, + IBAprep_REQUIRE_SAME_NCHANNELS | IBAprep_NO_SUPPORT_VOLUME | + IBAprep_NO_COPY_ROI_FULL)) + return false; + + const ImageSpec &srcspec (src.spec()); + + // Compute scaling factors and use action_resize to do the heavy lifting + int fit_full_width = roi.width(); + int fit_full_height = roi.height(); + int fit_full_x = roi.xbegin; + int fit_full_y = roi.ybegin; + float oldaspect = float(srcspec.full_width) / srcspec.full_height; + float newaspect = float(fit_full_width) / fit_full_height; + int resize_full_width = fit_full_width; + int resize_full_height = fit_full_height; + int xoffset = 0, yoffset = 0; + float xoff = 0.0f, yoff = 0.0f; + float scale = 1.0f; + + if (newaspect >= oldaspect) { // same or wider than original + resize_full_width = int(resize_full_height * oldaspect + 0.5f); + xoffset = (fit_full_width - resize_full_width) / 2; + scale = float(fit_full_height) / float(srcspec.full_height); + xoff = float(fit_full_width - scale * srcspec.full_width) / 2.0f; + } else { // narrower than original + resize_full_height = int(resize_full_width / oldaspect + 0.5f); + yoffset = (fit_full_height - resize_full_height) / 2; + scale = float(fit_full_width) / float(srcspec.full_width); + yoff = float(fit_full_height - scale * srcspec.full_height) / 2.0f; + } + + ROI newroi (fit_full_x, fit_full_x+fit_full_width, + fit_full_y, fit_full_y+fit_full_height, + 0, 1, 0, srcspec.nchannels); + // std::cout << " Fitting " << srcspec.roi() + // << " into " << newroi << "\n"; + // std::cout << " Fit scale factor " << scale << "\n"; + + // Set up a shared pointer with custom deleter to make sure any + // filter we allocate here is properly destroyed. + std::shared_ptr filterptr ((Filter2D*)NULL, Filter2D::destroy); + if (! filter) { + // If no filter was provided, punt and just linearly interpolate. + float wratio = float(resize_full_width) / float(srcspec.full_width); + float hratio = float(resize_full_height) / float(srcspec.full_height); + float w = 2.0f * std::max (1.0f, wratio); + float h = 2.0f * std::max (1.0f, hratio); + filter = Filter2D::create ("triangle", w, h); + filterptr.reset (filter); + } + + bool ok = true; + if (exact) { + // Full partial-pixel filtered resize -- exactly preserves aspect + // ratio and exactly centers the padded image, but might make the + // edges of the resized area blurry because it's not a whole number + // of pixels. + Imath::M33f M (scale, 0.0f, 0.0f, + 0.0f, scale, 0.0f, + xoff, yoff, 1.0f); + // std::cout << " Fit performing warp with " << M << "\n"; + ImageSpec newspec = srcspec; + newspec.set_roi (newroi); + newspec.set_roi_full (newroi); + dst.reset (newspec); + ImageBuf::WrapMode wrap = ImageBuf::WrapMode_from_string ("black"); + ok &= ImageBufAlgo::warp (dst, src, M, filter, false, wrap, {}, nthreads); + } else { + // Full pixel resize -- gives the sharpest result, but for odd-sized + // destination resolution, may not be exactly centered and will only + // preserve the aspect ratio to the nearest integer pixel size. + if (resize_full_width != srcspec.full_width || + resize_full_height != srcspec.full_height || + fit_full_x != srcspec.full_x || fit_full_y != srcspec.full_y) { + ROI resizeroi (fit_full_x, fit_full_x+resize_full_width, + fit_full_y, fit_full_y+resize_full_height, + 0, 1, 0, srcspec.nchannels); + ImageSpec newspec = srcspec; + newspec.set_roi (resizeroi); + newspec.set_roi_full (resizeroi); + dst.reset (newspec); + ok &= ImageBufAlgo::resize (dst, src, filter, resizeroi, nthreads); + } else { + ok &= dst.copy (src); // no resize is necessary + } + dst.specmod().full_width = fit_full_width; + dst.specmod().full_height = fit_full_height; + dst.specmod().full_x = fit_full_x; + dst.specmod().full_y = fit_full_y; + dst.specmod().x = xoffset; + dst.specmod().y = yoffset; + } + return ok; +} + + + +bool +ImageBufAlgo::fit (ImageBuf &dst, const ImageBuf &src, + string_view filtername, float fwidth, bool exact, + ROI roi, int nthreads) +{ + pvt::LoggedTimer logtime("IBA::fit"); + if (! IBAprep (roi, &dst, &src, + IBAprep_REQUIRE_SAME_NCHANNELS | IBAprep_NO_SUPPORT_VOLUME | + IBAprep_NO_COPY_ROI_FULL)) + return false; + const ImageSpec &srcspec (src.spec()); + const ImageSpec &dstspec (dst.spec()); + // Resize ratios + float wratio = float(dstspec.full_width) / float(srcspec.full_width); + float hratio = float(dstspec.full_height) / float(srcspec.full_height); + + // Set up a shared pointer with custom deleter to make sure any + // filter we allocate here is properly destroyed. + auto filter = get_resize_filter (filtername, fwidth, dst, wratio, hratio); + if (! filter) + return false; // error issued in get_resize_filter + + logtime.stop(); // it will be picked up again by the next call... + return fit (dst, src, filter.get(), exact, roi, nthreads); +} + + + +ImageBuf +ImageBufAlgo::fit (const ImageBuf &src, Filter2D *filter, bool exact, + ROI roi, int nthreads) +{ + ImageBuf result; + bool ok = fit (result, src, filter, exact, roi, nthreads); + if (!ok && !result.has_error()) + result.error ("ImageBufAlgo::fit() error"); + return result; +} + + +ImageBuf +ImageBufAlgo::fit (const ImageBuf &src, + string_view filtername, float filterwidth, bool exact, + ROI roi, int nthreads) +{ + ImageBuf result; + bool ok = fit (result, src, filtername, filterwidth, exact, roi, nthreads); + if (!ok && !result.has_error()) + result.error ("ImageBufAlgo::fit() error"); + return result; +} + + + +template +static bool +resample_ (ImageBuf &dst, const ImageBuf &src, bool interpolate, + ROI roi, int nthreads) +{ + ImageBufAlgo::parallel_image (roi, nthreads, [&](ROI roi){ + const ImageSpec &srcspec (src.spec()); + const ImageSpec &dstspec (dst.spec()); + int nchannels = src.nchannels(); + bool deep = src.deep(); + ASSERT (deep == dst.deep()); + + // Local copies of the source image window, converted to float + float srcfx = srcspec.full_x; + float srcfy = srcspec.full_y; + float srcfw = srcspec.full_width; + float srcfh = srcspec.full_height; + + float dstfx = dstspec.full_x; + float dstfy = dstspec.full_y; + float dstfw = dstspec.full_width; + float dstfh = dstspec.full_height; + float dstpixelwidth = 1.0f / dstfw; + float dstpixelheight = 1.0f / dstfh; + float *pel = ALLOCA (float, nchannels); + + ImageBuf::Iterator out (dst, roi); + ImageBuf::ConstIterator srcpel (src); + for (int y = roi.ybegin; y < roi.yend; ++y) { + // s,t are NDC space + float t = (y-dstfy+0.5f)*dstpixelheight; + // src_xf, src_xf are image space float coordinates + float src_yf = srcfy + t * srcfh; + // src_x, src_y are image space integer coordinates of the floor + int src_y = ifloor (src_yf); + for (int x = roi.xbegin; x < roi.xend; ++x, ++out) { + float s = (x-dstfx+0.5f)*dstpixelwidth; + float src_xf = srcfx + s * srcfw; + int src_x = ifloor (src_xf); + if (deep) { + srcpel.pos (src_x, src_y, 0); + int nsamps = srcpel.deep_samples(); + ASSERT (nsamps == out.deep_samples()); + if (! nsamps) + continue; + for (int c = 0; c < nchannels; ++c) { + if (dstspec.channelformat(c) == TypeDesc::UINT32) + for (int samp = 0; samp < nsamps; ++samp) + out.set_deep_value (c, samp, srcpel.deep_value_uint(c, samp)); + else + for (int samp = 0; samp < nsamps; ++samp) + out.set_deep_value (c, samp, srcpel.deep_value(c, samp)); + } + } else if (interpolate) { + // Non-deep image, bilinearly interpolate + src.interppixel (src_xf, src_yf, pel); + for (int c = roi.chbegin; c < roi.chend; ++c) + out[c] = pel[c]; + } else { + // Non-deep image, just copy closest pixel + srcpel.pos (src_x, src_y, 0); + for (int c = roi.chbegin; c < roi.chend; ++c) + out[c] = srcpel[c]; + } + } + } + }); + return true; +} + + + +bool +ImageBufAlgo::resample (ImageBuf &dst, const ImageBuf &src, + bool interpolate, ROI roi, int nthreads) +{ + pvt::LoggedTimer logtime("IBA::resample"); + if (! IBAprep (roi, &dst, &src, + IBAprep_REQUIRE_SAME_NCHANNELS | IBAprep_NO_SUPPORT_VOLUME | + IBAprep_NO_COPY_ROI_FULL | IBAprep_SUPPORT_DEEP)) + return false; + + if (dst.deep()) { + // If it's deep, figure out the sample allocations first, because + // it's not thread-safe to do that simultaneously with copying the + // values. + const ImageSpec &srcspec (src.spec()); + const ImageSpec &dstspec (dst.spec()); + float srcfx = srcspec.full_x; + float srcfy = srcspec.full_y; + float srcfw = srcspec.full_width; + float srcfh = srcspec.full_height; + float dstpixelwidth = 1.0f / dstspec.full_width; + float dstpixelheight = 1.0f / dstspec.full_height; + ImageBuf::ConstIterator srcpel (src, roi); + ImageBuf::Iterator dstpel (dst, roi); + for ( ; !dstpel.done(); ++dstpel, ++srcpel) { + float s = (dstpel.x()-dstspec.full_x+0.5f)*dstpixelwidth; + float t = (dstpel.y()-dstspec.full_y+0.5f)*dstpixelheight; + int src_y = ifloor (srcfy + t * srcfh); + int src_x = ifloor (srcfx + s * srcfw); + srcpel.pos (src_x, src_y, 0); + dstpel.set_deep_samples (srcpel.deep_samples ()); + } + } + + bool ok; + OIIO_DISPATCH_COMMON_TYPES2 (ok, "resample", resample_, + dst.spec().format, src.spec().format, + dst, src, interpolate, roi, nthreads); + return ok; +} + + + +ImageBuf +ImageBufAlgo::resample (const ImageBuf &src, bool interpolate, + ROI roi, int nthreads) +{ + ImageBuf result; + bool ok = resample (result, src, interpolate, roi, nthreads); + if (!ok && !result.has_error()) + result.error ("ImageBufAlgo::resample() error"); + return result; +} + + + + +#if 0 +template +static bool +affine_resample_ (ImageBuf &dst, const ImageBuf &src, const Imath::M33f &Minv, + ROI roi, int nthreads) +{ + ImageBufAlgo::parallel_image (roi, nthreads, [&](ROI roi){ + ImageBuf::Iterator d (dst, roi); + ImageBuf::ConstIterator s (src); + for ( ; ! d.done(); ++d) { + Imath::V2f P (d.x()+0.5f, d.y()+0.5f); + Minv.multVecMatrix (P, P); + s.pos (int(floorf(P.x)), int(floorf(P.y)), d.z()); + for (int c = roi.chbegin; c < roi.chend; ++c) + d[c] = s[c]; + } + }); + return true; +} +#endif + + + +template +static bool +warp_ (ImageBuf &dst, const ImageBuf &src, const Imath::M33f &M, + const Filter2D *filter, ImageBuf::WrapMode wrap, + ROI roi, int nthreads) +{ + ImageBufAlgo::parallel_image (roi, nthreads, [&](ROI roi){ + int nc = dst.nchannels(); + float *pel = ALLOCA (float, nc); + memset (pel, 0, nc*sizeof(float)); + Imath::M33f Minv = M.inverse(); + ImageBuf::Iterator out (dst, roi); + for ( ; ! out.done(); ++out) { + Dual2 x (out.x()+0.5f, 1.0f, 0.0f); + Dual2 y (out.y()+0.5f, 0.0f, 1.0f); + robust_multVecMatrix (Minv, x, y, x, y); + filtered_sample (src, x.val(), y.val(), + x.dx(), y.dx(), x.dy(), y.dy(), + filter, wrap, pel); + for (int c = roi.chbegin; c < roi.chend; ++c) + out[c] = pel[c]; + + } + }); + return true; +} + + + +bool +ImageBufAlgo::warp (ImageBuf &dst, const ImageBuf &src, + const Imath::M33f &M, + const Filter2D *filter, + bool recompute_roi, ImageBuf::WrapMode wrap, + ROI roi, int nthreads) +{ + pvt::LoggedTimer logtime("IBA::warp"); + ROI src_roi_full = src.roi_full(); + ROI dst_roi, dst_roi_full; + if (dst.initialized()) { + dst_roi = roi.defined() ? roi : dst.roi(); + dst_roi_full = dst.roi_full(); + } else { + dst_roi = roi.defined() ? roi : (recompute_roi ? transform (M, src.roi()) : src.roi()); + dst_roi_full = src_roi_full; + } + dst_roi.chend = std::min (dst_roi.chend, src.nchannels()); + dst_roi_full.chend = std::min (dst_roi_full.chend, src.nchannels()); + + if (! IBAprep (dst_roi, &dst, &src, IBAprep_NO_SUPPORT_VOLUME)) + return false; + + // Set up a shared pointer with custom deleter to make sure any + // filter we allocate here is properly destroyed. + std::shared_ptr filterptr ((Filter2D*)NULL, Filter2D::destroy); + if (filter == NULL) { + // If no filter was provided, punt and just linearly interpolate. + filterptr.reset (Filter2D::create ("lanczos3", 6.0f, 6.0f)); + filter = filterptr.get(); + } + + bool ok; + OIIO_DISPATCH_COMMON_TYPES2 (ok, "warp", warp_, + dst.spec().format, src.spec().format, + dst, src, M, filter, wrap, dst_roi, nthreads); + return ok; +} + + + +bool +ImageBufAlgo::warp (ImageBuf &dst, const ImageBuf &src, + const Imath::M33f &M, + string_view filtername_, float filterwidth, + bool recompute_roi, ImageBuf::WrapMode wrap, + ROI roi, int nthreads) +{ + // Set up a shared pointer with custom deleter to make sure any + // filter we allocate here is properly destroyed. + std::shared_ptr filter ((Filter2D*)NULL, Filter2D::destroy); + std::string filtername = filtername_.size() ? filtername_ : "lanczos3"; + for (int i = 0, e = Filter2D::num_filters(); i < e; ++i) { + FilterDesc fd; + Filter2D::get_filterdesc (i, &fd); + if (fd.name == filtername) { + float w = filterwidth > 0.0f ? filterwidth : fd.width; + float h = filterwidth > 0.0f ? filterwidth : fd.width; + filter.reset (Filter2D::create (filtername, w, h)); + break; + } + } + if (! filter) { + dst.error ("Filter \"%s\" not recognized", filtername); + return false; + } + + return warp (dst, src, M, filter.get(), recompute_roi, wrap, roi, nthreads); +} + + + +ImageBuf +ImageBufAlgo::warp (const ImageBuf &src, + const Imath::M33f &M, const Filter2D *filter, + bool recompute_roi, ImageBuf::WrapMode wrap, + ROI roi, int nthreads) +{ + ImageBuf result; + bool ok = warp (result, src, M, filter, recompute_roi, wrap, roi, nthreads); + if (!ok && !result.has_error()) + result.error ("ImageBufAlgo::warp() error"); + return result; +} + + + +ImageBuf +ImageBufAlgo::warp (const ImageBuf &src, const Imath::M33f &M, + string_view filtername, float filterwidth, + bool recompute_roi, ImageBuf::WrapMode wrap, + ROI roi, int nthreads) +{ + ImageBuf result; + bool ok = warp (result, src, M, filtername, filterwidth, + recompute_roi, wrap, roi, nthreads); + if (!ok && !result.has_error()) + result.error ("ImageBufAlgo::warp() error"); + return result; +} + + + +bool +ImageBufAlgo::rotate (ImageBuf &dst, const ImageBuf &src, + float angle, float center_x, float center_y, + Filter2D *filter, bool recompute_roi, + ROI roi, int nthreads) +{ + // Calculate the rotation matrix + Imath::M33f M; + M.translate(Imath::V2f(-center_x, -center_y)); + M.rotate(angle); + M *= Imath::M33f().translate(Imath::V2f(center_x, center_y)); + return ImageBufAlgo::warp (dst, src, M, filter, + recompute_roi, ImageBuf::WrapBlack, + roi, nthreads); +} + + + +bool +ImageBufAlgo::rotate (ImageBuf &dst, const ImageBuf &src, + float angle, float center_x, float center_y, + string_view filtername, float filterwidth, + bool recompute_roi, ROI roi, int nthreads) +{ + // Calculate the rotation matrix + Imath::M33f M; + M.translate(Imath::V2f(-center_x, -center_y)); + M.rotate(angle); + M *= Imath::M33f().translate(Imath::V2f(center_x, center_y)); + return ImageBufAlgo::warp (dst, src, M, filtername, filterwidth, + recompute_roi, ImageBuf::WrapBlack, + roi, nthreads); +} + + + +bool +ImageBufAlgo::rotate (ImageBuf &dst, const ImageBuf &src, float angle, + Filter2D *filter, + bool recompute_roi, ROI roi, int nthreads) +{ + ROI src_roi_full = src.roi_full(); + float center_x = 0.5f * (src_roi_full.xbegin + src_roi_full.xend); + float center_y = 0.5f * (src_roi_full.ybegin + src_roi_full.yend); + return ImageBufAlgo::rotate (dst, src, angle, center_x, center_y, + filter, recompute_roi, roi, nthreads); +} + + + +bool +ImageBufAlgo::rotate (ImageBuf &dst, const ImageBuf &src, float angle, + string_view filtername, float filterwidth, + bool recompute_roi, + ROI roi, int nthreads) +{ + ROI src_roi_full = src.roi_full(); + float center_x = 0.5f * (src_roi_full.xbegin + src_roi_full.xend); + float center_y = 0.5f * (src_roi_full.ybegin + src_roi_full.yend); + return ImageBufAlgo::rotate (dst, src, angle, center_x, center_y, + filtername, filterwidth, recompute_roi, + roi, nthreads); +} + + + +ImageBuf +ImageBufAlgo::rotate (const ImageBuf &src, + float angle, float center_x, float center_y, + Filter2D *filter, bool recompute_roi, + ROI roi, int nthreads) +{ + ImageBuf result; + bool ok = rotate (result, src, angle, center_x, center_y, + filter, recompute_roi, roi, nthreads); + if (!ok && !result.has_error()) + result.error ("ImageBufAlgo::rotate() error"); + return result; +} + + + +ImageBuf +ImageBufAlgo::rotate (const ImageBuf &src, + float angle, float center_x, float center_y, + string_view filtername, float filterwidth, + bool recompute_roi, ROI roi, int nthreads) +{ + ImageBuf result; + bool ok = rotate (result, src, angle, center_x, center_y, + filtername, filterwidth, recompute_roi, roi, nthreads); + if (!ok && !result.has_error()) + result.error ("ImageBufAlgo::rotate() error"); + return result; +} + + + +ImageBuf +ImageBufAlgo::rotate (const ImageBuf &src, float angle, + Filter2D *filter, + bool recompute_roi, ROI roi, int nthreads) +{ + ImageBuf result; + bool ok = rotate (result, src, angle, filter, recompute_roi, roi, nthreads); + if (!ok && !result.has_error()) + result.error ("ImageBufAlgo::rotate() error"); + return result; +} + + + +ImageBuf +ImageBufAlgo::rotate (const ImageBuf &src, float angle, + string_view filtername, float filterwidth, + bool recompute_roi, + ROI roi, int nthreads) +{ + ImageBuf result; + bool ok = rotate (result, src, angle, filtername, filterwidth, + recompute_roi, roi, nthreads); + if (!ok && !result.has_error()) + result.error ("ImageBufAlgo::rotate() error"); + return result; +} + + +OIIO_NAMESPACE_END