Added lbb resampler; Removed mbicubic resampler

libvips/resample/Makefile.am

@@ -2,7 +2,7 @@
 if ENABLE_CXX
 C_SOURCES = \
 	bicubic.cpp \
-	mbicubic.cpp \
+	lbb.cpp \
 	yafrsmooth.cpp \
 	nohalo1.cpp \
 	snohalo1.cpp \
@@ -13,7 +13,7 @@ else
 C_SOURCES = 
 C_DIST = \
 	bicubic.cpp \
-	mbicubic.cpp \
+	lbb.cpp \
 	yafrsmooth.cpp \
 	nohalo1.cpp \
 	snohalo1.cpp \

libvips/resample/interpolate.c

@@ -497,7 +497,7 @@ void
 vips__interpolate_init( void )
 {
 	extern GType vips_interpolate_bicubic_get_type( void );
-	extern GType vips_interpolate_mbicubic_get_type( void );
+	extern GType vips_interpolate_lbb_get_type( void );
 	extern GType vips_interpolate_yafrsmooth_get_type( void );
 	extern GType vips_interpolate_nohalo1_get_type( void );
 	extern GType vips_interpolate_snohalo1_get_type( void );
@@ -508,7 +508,7 @@ vips__interpolate_init( void )
 
 #ifdef ENABLE_CXX
 	vips_interpolate_bicubic_get_type();
-	vips_interpolate_mbicubic_get_type();
+	vips_interpolate_lbb_get_type();
 	vips_interpolate_yafrsmooth_get_type();
 	vips_interpolate_nohalo1_get_type();
 	vips_interpolate_snohalo1_get_type();

libvips/resample/lbb.cpp

@@ -0,0 +1,739 @@
+/* locally bounded bicubic resampler
+ */
+
+/*
+
+    This file is part of VIPS.
+
+    VIPS is free software; you can redistribute it and/or modify it
+    under the terms of the GNU Lesser General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    Lesser General Public License for more details.
+
+    You should have received a copy of the GNU Lesser General Public
+    License along with this program; if not, write to the Free
+    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+    02111-1307 USA
+
+ */
+
+/*
+
+    These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
+
+ */
+
+/*
+ * 2009-2010 (c) Nicolas Robidoux, John Cupitt, Chantal Racette.
+ *
+ * Nicolas Robidoux thanks Ralf Meyer, Minglun Gong, Adam Turcotte,
+ * Eric Daoust, Øyvind Kolås, Geert Jordaens, and Sven Neumann for
+ * useful comments and code.
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif /*HAVE_CONFIG_H*/
+#include <vips/intl.h>
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <vips/vips.h>
+#include <vips/internal.h>
+
+#include "templates.h"
+
+#define VIPS_TYPE_INTERPOLATE_LBB \
+	(vips_interpolate_lbb_get_type())
+#define VIPS_INTERPOLATE_LBB( obj ) \
+	(G_TYPE_CHECK_INSTANCE_CAST( (obj), \
+	VIPS_TYPE_INTERPOLATE_LBB, VipsInterpolateLbb ))
+#define VIPS_INTERPOLATE_LBB_CLASS( klass ) \
+	(G_TYPE_CHECK_CLASS_CAST( (klass), \
+	VIPS_TYPE_INTERPOLATE_LBB, VipsInterpolateLbbClass))
+#define VIPS_IS_INTERPOLATE_LBB( obj ) \
+	(G_TYPE_CHECK_INSTANCE_TYPE( (obj), VIPS_TYPE_INTERPOLATE_LBB ))
+#define VIPS_IS_INTERPOLATE_LBB_CLASS( klass ) \
+	(G_TYPE_CHECK_CLASS_TYPE( (klass), VIPS_TYPE_INTERPOLATE_LBB ))
+#define VIPS_INTERPOLATE_LBB_GET_CLASS( obj ) \
+	(G_TYPE_INSTANCE_GET_CLASS( (obj), \
+	VIPS_TYPE_INTERPOLATE_LBB, VipsInterpolateLbbClass ))
+
+typedef struct _VipsInterpolateLbb {
+	VipsInterpolate parent_object;
+
+} VipsInterpolateLbb;
+
+typedef struct _VipsInterpolateLbbClass {
+	VipsInterpolateClass parent_class;
+
+} VipsInterpolateLbbClass;
+
+static inline double
+lbbicubic( const double c00,
+           const double c10,
+           const double c01,
+           const double c11,
+           const double c00dx,
+           const double c10dx,
+           const double c01dx,
+           const double c11dx,
+           const double c00dy,
+           const double c10dy,
+           const double c01dy,
+           const double c11dy,
+           const double c00dxdy,
+           const double c10dxdy,
+           const double c01dxdy,
+           const double c11dxdy,
+           const double uno_one,
+           const double uno_two,
+           const double uno_thr,
+           const double uno_fou,
+           const double dos_one,
+           const double dos_two,
+           const double dos_thr,
+           const double dos_fou,
+           const double tre_one,
+           const double tre_two,
+           const double tre_thr,
+           const double tre_fou,
+           const double qua_one,
+           const double qua_two,
+           const double qua_thr,
+           const double qua_fou )
+{
+  /*
+   * STENCIL (FOOTPRINT) OF INPUT VALUES:
+   *
+   * The stencil of Symmetrized Monotone Catmull-Rom is the same as
+   * the standard Catmull-Rom's:
+   *
+   *  (ix-1,iy-1)  (ix,iy-1)    (ix+1,iy-1)  (ix+2,iy-1)
+   *  = uno_one    = uno_two    = uno_thr    = uno_fou
+   *
+   *  (ix-1,iy)    (ix,iy)      (ix+1,iy)    (ix+2,iy)
+   *  = dos_one    = dos_two    = dos_thr    = dos_fou
+   *                        X
+   *  (ix-1,iy+1)  (ix,iy+1)    (ix+1,iy+1)  (ix+2,iy+1)
+   *  = tre_one    = tre_two    = tre_thr    = tre_fou
+   *
+   *  (ix-1,iy+2)  (ix,iy+2)    (ix+1,iy+2)  (ix+2,iy+2)
+   *  = qua_one    = qua_two    = qua_thr    = qua_fou
+   *
+   * where ix is the (pseudo-)floor of the requested left-to-right
+   * location ("X"), and iy is the floor of the requested up-to-down
+   * location.
+   */
+
+  /*
+   * Computation of the four min and four max over 3x3 input data
+   * sub-blocks of the 4x4 input stencil (involves 28 flag
+   * computations):
+   */
+  const double m1    = (dos_two <= dos_thr) ? dos_two : dos_thr;
+  const double M1    = (dos_two <= dos_thr) ? dos_thr : dos_two;
+  const double m2    = (tre_two <= tre_thr) ? tre_two : tre_thr;
+  const double M2    = (tre_two <= tre_thr) ? tre_thr : tre_two;
+  const double m3    = (uno_two <= uno_thr) ? uno_two : uno_thr;
+  const double M3    = (uno_two <= uno_thr) ? uno_thr : uno_two;
+  const double m4    = (qua_two <= qua_thr) ? qua_two : qua_thr;
+  const double M4    = (qua_two <= qua_thr) ? qua_thr : qua_two;
+  const double m5    = (m1      <= m2     ) ? m1      : m2     ;
+  const double M5    = (M1      >= M2     ) ? M1      : M2     ;
+  const double m6    = (dos_one <= tre_one) ? dos_one : tre_one;
+  const double M6    = (dos_one <= tre_one) ? tre_one : dos_one;
+  const double m7    = (dos_fou <= tre_fou) ? dos_fou : tre_fou;
+  const double M7    = (dos_fou <= tre_fou) ? tre_fou : dos_fou;
+  const double m8    = (m5      <= m3     ) ? m5      : m3     ;
+  const double M8    = (M5      >= M3     ) ? M5      : M3     ;
+  const double m9    = (m5      <= m4     ) ? m5      : m4     ;
+  const double M9    = (M5      >= M4     ) ? M5      : M4     ;
+  const double m10   = (m6      <= uno_one) ? m6      : uno_one;
+  const double M10   = (M6      >= uno_one) ? M6      : uno_one;
+  const double m11   = (m7      <= uno_fou) ? m7      : uno_fou;
+  const double M11   = (M7      >= uno_fou) ? M7      : uno_fou;
+  const double m12   = (m6      <= qua_one) ? m6      : qua_one;
+  const double M12   = (M6      >= qua_one) ? M6      : qua_one;
+  const double m13   = (m7      <= qua_fou) ? m7      : qua_fou;
+  const double M13   = (M7      >= qua_fou) ? M7      : qua_fou;
+  const double min00 = (m8      <= m10    ) ? m8      : m10    ;
+  const double max00 = (M8      >= M10    ) ? M8      : M10    ;
+  const double min10 = (m8      <= m11    ) ? m8      : m11    ;
+  const double max10 = (M8      >= M11    ) ? M8      : M11    ;
+  const double min01 = (m9      <= m12    ) ? m9      : m12    ;
+  const double max01 = (M9      >= M12    ) ? M9      : M12    ;
+  const double min11 = (m9      <= m13    ) ? m9      : m13    ;
+  const double max11 = (M9      >= M13    ) ? M9      : M13    ;
+
+  /*
+   * Distances to the local min and max:
+   */
+  const double u00 = dos_two - min00;
+  const double v00 = max00 - dos_two;
+  const double u10 = dos_thr - min10;
+  const double v10 = max10 - dos_thr;
+  const double u01 = tre_two - min01;
+  const double v01 = max01 - tre_two;
+  const double u11 = tre_thr - min11;
+  const double v11 = max11 - tre_thr;
+
+  /*
+   * Initial values of the derivatives computed with centered
+   * differences. Factors of 1/2 are left out because they are folded
+   * in later:
+   */
+  const double dble_dzdx00i = dos_thr - dos_one;
+  const double dble_dzdx10i = dos_fou - dos_two;
+  const double dble_dzdx01i = tre_thr - tre_one;
+  const double dble_dzdx11i = tre_fou - tre_two;
+
+  const double dble_dzdy00i = tre_two - uno_two;
+  const double dble_dzdy10i = tre_thr - uno_thr;
+  const double dble_dzdy01i = qua_two - dos_two;
+  const double dble_dzdy11i = qua_thr - dos_thr;
+
+  /*
+   * Signs of the derivatives:
+   */
+  const double sign_dzdx00 = (dble_dzdx00i >= 0.) ? 1. : -1.;
+  const double sign_dzdx10 = (dble_dzdx10i >= 0.) ? 1. : -1.;
+  const double sign_dzdx01 = (dble_dzdx01i >= 0.) ? 1. : -1.;
+  const double sign_dzdx11 = (dble_dzdx11i >= 0.) ? 1. : -1.;
+
+  const double sign_dzdy00 = (dble_dzdy00i >= 0.) ? 1. : -1.;
+  const double sign_dzdy10 = (dble_dzdy10i >= 0.) ? 1. : -1.;
+  const double sign_dzdy01 = (dble_dzdy01i >= 0.) ? 1. : -1.;
+  const double sign_dzdy11 = (dble_dzdy11i >= 0.) ? 1. : -1.;
+
+  /*
+   * Slope limiters. The key multiplier is 3 but we fold a factor of
+   * 2, hence 6:
+   */
+  const double dble_slopelimit_00 = 6.0 * ( (u00 <= v00) ? u00 : v00 );
+  const double dble_slopelimit_10 = 6.0 * ( (u10 <= v10) ? u10 : v10 );
+  const double dble_slopelimit_01 = 6.0 * ( (u01 <= v01) ? u01 : v01 );
+  const double dble_slopelimit_11 = 6.0 * ( (u11 <= v11) ? u11 : v11 );
+
+  /*
+   * Initial values of the cross-derivatives. Factors of 1/4 are left
+   * out because folded in later:
+   */
+  const double quad_d2zdxdy00i = ( uno_one - uno_thr ) + dble_dzdx01i;
+  const double quad_d2zdxdy10i = ( uno_two - uno_fou ) + dble_dzdx11i;
+  const double quad_d2zdxdy01i = ( qua_thr - qua_one ) - dble_dzdx00i;
+  const double quad_d2zdxdy11i = ( qua_fou - qua_two ) - dble_dzdx10i;
+
+  /*
+   * Part of the result which does not need derivatives:
+   */
+  const double newval1 = c00 * dos_two
+                         +
+                         c10 * dos_thr
+                         +
+                         c01 * tre_two
+                         +
+                         c11 * tre_thr;
+
+  /*
+   * Clamped first derivatives:
+   */
+  const double dble_dzdx00 =
+    ( sign_dzdx00 * dble_dzdx00i <= dble_slopelimit_00 )
+    ? dble_dzdx00i :  sign_dzdx00 * dble_slopelimit_00;
+  const double dble_dzdx10 =
+    ( sign_dzdx10 * dble_dzdx10i <= dble_slopelimit_10 )
+    ? dble_dzdx10i :  sign_dzdx10 * dble_slopelimit_10;
+  const double dble_dzdx01 =
+    ( sign_dzdx01 * dble_dzdx01i <= dble_slopelimit_01 )
+    ? dble_dzdx01i :  sign_dzdx01 * dble_slopelimit_01;
+  const double dble_dzdx11 =
+    ( sign_dzdx11 * dble_dzdx11i <= dble_slopelimit_11 )
+    ? dble_dzdx11i :  sign_dzdx11 * dble_slopelimit_11;
+  const double dble_dzdy00 =
+    ( sign_dzdy00 * dble_dzdy00i <= dble_slopelimit_00 )
+    ? dble_dzdy00i :  sign_dzdy00 * dble_slopelimit_00;
+  const double dble_dzdy10 =
+    ( sign_dzdy10 * dble_dzdy10i <= dble_slopelimit_10 )
+    ? dble_dzdy10i :  sign_dzdy10 * dble_slopelimit_10;
+  const double dble_dzdy01 =
+    ( sign_dzdy01 * dble_dzdy01i <= dble_slopelimit_01 )
+    ? dble_dzdy01i :  sign_dzdy01 * dble_slopelimit_01;
+  const double dble_dzdy11 =
+    ( sign_dzdy11 * dble_dzdy11i <= dble_slopelimit_11 )
+    ? dble_dzdy11i :  sign_dzdy11 * dble_slopelimit_11;
+
+  /*
+   * Sums and differences of first derivatives:
+   */
+  const double twelve_sum00 = 6.0 * ( dble_dzdx00 + dble_dzdy00 );
+  const double twelve_dif00 = 6.0 * ( dble_dzdx00 - dble_dzdy00 );
+  const double twelve_sum10 = 6.0 * ( dble_dzdx10 + dble_dzdy10 );
+  const double twelve_dif10 = 6.0 * ( dble_dzdx10 - dble_dzdy10 );
+  const double twelve_sum01 = 6.0 * ( dble_dzdx01 + dble_dzdy01 );
+  const double twelve_dif01 = 6.0 * ( dble_dzdx01 - dble_dzdy01 );
+  const double twelve_sum11 = 6.0 * ( dble_dzdx11 + dble_dzdy11 );
+  const double twelve_dif11 = 6.0 * ( dble_dzdx11 - dble_dzdy11 );
+
+  /*
+   * Part of the result which only needs first derivatives.
+   */
+  const double newval2 = c00dx * dble_dzdx00
+                         +
+                         c10dx * dble_dzdx10
+                         +
+                         c01dx * dble_dzdx01
+                         +
+                         c11dx * dble_dzdx11
+                         +
+                         c00dy * dble_dzdy00
+                         +
+                         c10dy * dble_dzdy10
+                         +
+                         c01dy * dble_dzdy01
+                         +
+                         c11dy * dble_dzdy11;
+
+  /*
+   * Absolute values of the sums:
+   */
+  const double twelve_abs_sum00 =
+    (twelve_sum00 >= 0.0) ? twelve_sum00 : -twelve_sum00;
+  const double twelve_abs_sum10 =
+    (twelve_sum10 >= 0.0) ? twelve_sum10 : -twelve_sum10;
+  const double twelve_abs_sum01 =
+    (twelve_sum01 >= 0.0) ? twelve_sum01 : -twelve_sum01;
+  const double twelve_abs_sum11 =
+    (twelve_sum11 >= 0.0) ? twelve_sum11 : -twelve_sum11;
+
+  /*
+   * Scaled 'u' differences:
+   */
+  const double u00_times_36 = 36. * u00;
+  const double u10_times_36 = 36. * u10;
+  const double u01_times_36 = 36. * u01;
+  const double u11_times_36 = 36. * u11;
+
+  /*
+   * First cross-derivative limiter:
+   */
+  const double first_limit00  = twelve_abs_sum00 - u00_times_36;
+  const double first_limit10  = twelve_abs_sum10 - u10_times_36;
+  const double first_limit01  = twelve_abs_sum01 - u01_times_36;
+  const double first_limit11  = twelve_abs_sum11 - u11_times_36;
+
+  const double quad_d2zdxdy00ii =
+    (quad_d2zdxdy00i >= first_limit00)
+    ? quad_d2zdxdy00i : first_limit00;
+  const double quad_d2zdxdy10ii =
+    (quad_d2zdxdy10i >= first_limit10)
+    ? quad_d2zdxdy10i : first_limit10;
+  const double quad_d2zdxdy01ii =
+    (quad_d2zdxdy01i >= first_limit01)
+    ? quad_d2zdxdy01i : first_limit01;
+  const double quad_d2zdxdy11ii =
+    (quad_d2zdxdy11i >= first_limit11)
+    ? quad_d2zdxdy11i : first_limit11;
+
+  /*
+   * Absolute values of the differences:
+   */
+  const double twelve_abs_dif00 =
+    (twelve_dif00 >= 0.0) ? twelve_dif00 : -twelve_dif00;
+  const double twelve_abs_dif10 =
+    (twelve_dif10 >= 0.0) ? twelve_dif10 : -twelve_dif10;
+  const double twelve_abs_dif01 =
+    (twelve_dif01 >= 0.0) ? twelve_dif01 : -twelve_dif01;
+  const double twelve_abs_dif11 =
+    (twelve_dif11 >= 0.0) ? twelve_dif11 : -twelve_dif11;
+
+  /*
+   * Scaled 'v' differences:
+   */
+  const double v00_times_36 = 36. * v00;
+  const double v10_times_36 = 36. * v10;
+  const double v01_times_36 = 36. * v01;
+  const double v11_times_36 = 36. * v11;
+
+  /*
+   * Second cross-derivative limiter:
+   */
+  const double second_limit00 = v00_times_36 - twelve_abs_sum00;
+  const double second_limit10 = v10_times_36 - twelve_abs_sum10;
+  const double second_limit01 = v01_times_36 - twelve_abs_sum01;
+  const double second_limit11 = v11_times_36 - twelve_abs_sum11;
+
+  const double quad_d2zdxdy00iii =
+    (quad_d2zdxdy00ii <= second_limit00)
+    ? quad_d2zdxdy00ii : second_limit00;
+  const double quad_d2zdxdy10iii =
+    (quad_d2zdxdy10ii <= second_limit10)
+    ? quad_d2zdxdy10ii : second_limit10;
+  const double quad_d2zdxdy01iii =
+    (quad_d2zdxdy01ii <= second_limit01)
+    ? quad_d2zdxdy01ii : second_limit01;
+  const double quad_d2zdxdy11iii =
+    (quad_d2zdxdy11ii <= second_limit11)
+    ? quad_d2zdxdy11ii : second_limit11;
+
+  /*
+   * Third cross-derivative limiter:
+   */
+  const double third_limit00  = u00_times_36 - twelve_abs_dif00;
+  const double third_limit10  = u10_times_36 - twelve_abs_dif10;
+  const double third_limit01  = u01_times_36 - twelve_abs_dif01;
+  const double third_limit11  = u11_times_36 - twelve_abs_dif11;
+
+  const double quad_d2zdxdy00iiii =
+    (quad_d2zdxdy00iii <= third_limit00)
+    ? quad_d2zdxdy00iii : third_limit00;
+  const double quad_d2zdxdy10iiii =
+    (quad_d2zdxdy10iii <= third_limit10)
+    ? quad_d2zdxdy10iii : third_limit10;
+  const double quad_d2zdxdy01iiii =
+    (quad_d2zdxdy01iii <= third_limit01)
+    ? quad_d2zdxdy01iii : third_limit01;
+  const double quad_d2zdxdy11iiii =
+    (quad_d2zdxdy11iii <= third_limit11)
+    ? quad_d2zdxdy11iii : third_limit11;
+
+  /*
+   * Fourth cross-derivative limiter:
+   */
+  const double fourth_limit00 = twelve_abs_dif00 - v00_times_36;
+  const double fourth_limit10 = twelve_abs_dif10 - v10_times_36;
+  const double fourth_limit01 = twelve_abs_dif01 - v01_times_36;
+  const double fourth_limit11 = twelve_abs_dif11 - v11_times_36;
+
+  const double quad_d2zdxdy00 =
+    (quad_d2zdxdy00iiii >= fourth_limit00)
+    ? quad_d2zdxdy00iiii : fourth_limit00;
+  const double quad_d2zdxdy10 =
+    (quad_d2zdxdy10iiii >= fourth_limit10)
+    ? quad_d2zdxdy10iiii : fourth_limit10;
+  const double quad_d2zdxdy01 =
+    (quad_d2zdxdy01iiii >= fourth_limit01)
+    ? quad_d2zdxdy01iiii : fourth_limit01;
+  const double quad_d2zdxdy11 =
+    (quad_d2zdxdy11iiii >= fourth_limit11)
+    ? quad_d2zdxdy11iiii : fourth_limit11;
+
+  /*
+   * Four times the part of the result which only uses cross
+   * derivatives:
+   */
+  const double newval3 = c00dxdy * quad_d2zdxdy00
+                         +
+                         c10dxdy * quad_d2zdxdy10
+                         +
+                         c01dxdy * quad_d2zdxdy01
+                         +
+                         c11dxdy * quad_d2zdxdy11;
+
+  const double newval = newval1 + .5 * newval2 + .25 * newval3;
+
+  return newval;
+}
+
+/*
+ * Call lbb with a type conversion operator as a parameter.
+ *
+ * It would be nice to do this with templates somehow---for one thing
+ * this would allow code comments!---but we can't figure a clean way
+ * to do it.
+ */
+#define LBB_CONVERSION( conversion )                          \
+  template <typename T> static void inline               \
+  lbb_ ## conversion(       PEL*   restrict pout,        \
+                      const PEL*   restrict pin,         \
+                      const int             bands,       \
+                      const int             lskip,       \
+                      const double          relative_x,  \
+                      const double          relative_y ) \
+  { \
+    T* restrict out = (T *) pout; \
+    \
+    const T* restrict in = (T *) pin; \
+    \
+    const int one_shift =  -bands; \
+    const int thr_shift =   bands; \
+    const int fou_shift = 2*bands; \
+    \
+    const int uno_two_shift =  -lskip; \
+    const int dos_two_shift =       0; \
+    const int tre_two_shift =   lskip; \
+    const int qua_two_shift = 2*lskip; \
+    \
+    const int uno_one_shift = uno_two_shift + one_shift; \
+    const int dos_one_shift = dos_two_shift + one_shift; \
+    const int tre_one_shift = tre_two_shift + one_shift; \
+    const int qua_one_shift = qua_two_shift + one_shift; \
+    \
+    const int uno_thr_shift = uno_two_shift + thr_shift; \
+    const int dos_thr_shift = dos_two_shift + thr_shift; \
+    const int tre_thr_shift = tre_two_shift + thr_shift; \
+    const int qua_thr_shift = qua_two_shift + thr_shift; \
+    \
+    const int uno_fou_shift = uno_two_shift + fou_shift; \
+    const int dos_fou_shift = dos_two_shift + fou_shift; \
+    const int tre_fou_shift = tre_two_shift + fou_shift; \
+    const int qua_fou_shift = qua_two_shift + fou_shift; \
+    \
+    const double xp1over2   = relative_x; \
+    const double xm1over2   = xp1over2 - 1.0; \
+    const double onemx      = 1.5 - xp1over2; \
+    const double onepx      = 0.5 + xp1over2; \
+    const double xp1over2sq = xp1over2 * xp1over2; \
+    \
+    const double yp1over2   = relative_y; \
+    const double ym1over2   = yp1over2 - 1.0; \
+    const double onemy      = 1.5 - yp1over2; \
+    const double onepy      = 0.5 + yp1over2; \
+    const double yp1over2sq = yp1over2 * yp1over2; \
+    \
+    const double xm1over2sq = xm1over2 * xm1over2; \
+    const double ym1over2sq = ym1over2 * ym1over2; \
+    \
+    const double twice1mx = onemx + onemx; \
+    const double twice1px = onepx + onepx; \
+    const double twice1my = onemy + onemy; \
+    const double twice1py = onepy + onepy; \
+    \
+    const double xm1over2sq_times_ym1over2sq = xm1over2sq * ym1over2sq; \
+    const double xp1over2sq_times_ym1over2sq = xp1over2sq * ym1over2sq; \
+    const double xp1over2sq_times_yp1over2sq = xp1over2sq * yp1over2sq; \
+    const double xm1over2sq_times_yp1over2sq = xm1over2sq * yp1over2sq; \
+    \
+    const double xm1over2_times_ym1over2  = xm1over2 * ym1over2; \
+    const double xp1over2_times_ym1over2  = xp1over2 * ym1over2; \
+    const double twice_1mx_times_ym1over2 = twice1mx * ym1over2; \
+    const double twice_1px_times_ym1over2 = twice1px * ym1over2; \
+    \
+    const double xm1over2_times_yp1over2  = xm1over2 * yp1over2; \
+    const double xp1over2_times_yp1over2  = xp1over2 * yp1over2; \
+    const double twice_1mx_times_yp1over2 = twice1mx * yp1over2; \
+    const double twice_1px_times_yp1over2 = twice1px * yp1over2; \
+    \
+    const double twice_xm1over2_times_1my = xm1over2 * twice1my; \
+    const double twice_xp1over2_times_1my = xp1over2 * twice1my; \
+    const double four_times_1mx_times_1my = twice1mx * twice1my; \
+    const double four_times_1px_times_1my = twice1px * twice1my; \
+    \
+    const double twice_xm1over2_times_1py = xm1over2 * twice1py; \
+    const double twice_xp1over2_times_1py = xp1over2 * twice1py; \
+    const double four_times_1mx_times_1py = twice1mx * twice1py; \
+    const double four_times_1px_times_1py = twice1px * twice1py; \
+    \
+    const double c00     = \
+      four_times_1px_times_1py * xm1over2sq_times_ym1over2sq; \
+    const double c00dx   = \
+      twice_xp1over2_times_1py * xm1over2sq_times_ym1over2sq; \
+    const double c00dy   = \
+      twice_1px_times_yp1over2 * xm1over2sq_times_ym1over2sq; \
+    const double c00dxdy = \
+       xp1over2_times_yp1over2 * xm1over2sq_times_ym1over2sq; \
+    \
+    const double c10     = \
+      four_times_1mx_times_1py * xp1over2sq_times_ym1over2sq; \
+    const double c10dx   = \
+      twice_xm1over2_times_1py * xp1over2sq_times_ym1over2sq; \
+    const double c10dy   = \
+      twice_1mx_times_yp1over2 * xp1over2sq_times_ym1over2sq; \
+    const double c10dxdy = \
+       xm1over2_times_yp1over2 * xp1over2sq_times_ym1over2sq; \
+    \
+    const double c01     = \
+      four_times_1px_times_1my * xm1over2sq_times_yp1over2sq; \
+    const double c01dx   = \
+      twice_xp1over2_times_1my * xm1over2sq_times_yp1over2sq; \
+    const double c01dy   = \
+      twice_1px_times_ym1over2 * xm1over2sq_times_yp1over2sq; \
+    const double c01dxdy = \
+       xp1over2_times_ym1over2 * xm1over2sq_times_yp1over2sq; \
+    \
+    const double c11     = \
+      four_times_1mx_times_1my * xp1over2sq_times_yp1over2sq; \
+    const double c11dx   = \
+      twice_xm1over2_times_1my * xp1over2sq_times_yp1over2sq; \
+    const double c11dy   = \
+      twice_1mx_times_ym1over2 * xp1over2sq_times_yp1over2sq; \
+    const double c11dxdy = \
+       xm1over2_times_ym1over2 * xp1over2sq_times_yp1over2sq; \
+    \
+    int band = bands; \
+    \
+    do \
+      { \
+        const double double_result =        \
+          lbbicubic( c00,                   \
+                     c10,                   \
+                     c01,                   \
+                     c11,                   \
+                     c00dx,                 \
+                     c10dx,                 \
+                     c01dx,                 \
+                     c11dx,                 \
+                     c00dy,                 \
+                     c10dy,                 \
+                     c01dy,                 \
+                     c11dy,                 \
+                     c00dxdy,               \
+                     c10dxdy,               \
+                     c01dxdy,               \
+                     c11dxdy,               \
+                     in[ uno_one_shift ],   \
+                     in[ uno_two_shift ],   \
+                     in[ uno_thr_shift ],   \
+                     in[ uno_fou_shift ],   \
+                     in[ dos_one_shift ],   \
+                     in[ dos_two_shift ],   \
+                     in[ dos_thr_shift ],   \
+                     in[ dos_fou_shift ],   \
+                     in[ tre_one_shift ],   \
+                     in[ tre_two_shift ],   \
+                     in[ tre_thr_shift ],   \
+                     in[ tre_fou_shift ],   \
+                     in[ qua_one_shift ],   \
+                     in[ qua_two_shift ],   \
+                     in[ qua_thr_shift ],   \
+                     in[ qua_fou_shift ] ); \
+        \
+        const T result = to_ ## conversion<T>( double_result ); \
+        in++; \
+        *out++ = result; \
+      } while (--band); \
+  }
+
+LBB_CONVERSION( fptypes )
+LBB_CONVERSION( withsign )
+LBB_CONVERSION( nosign )
+
+#define CALL( T, conversion )          \
+  lbb_ ## conversion<T>( out,          \
+                         p,            \
+                         bands,        \
+                         lskip,        \
+                         relative_x,   \
+                         relative_y );
+
+/*
+ * We need C linkage:
+ */
+extern "C" {
+G_DEFINE_TYPE( VipsInterpolateLbb, vips_interpolate_lbb,
+	VIPS_TYPE_INTERPOLATE );
+}
+
+static void
+vips_interpolate_lbb_interpolate( VipsInterpolate* restrict interpolate,
+                                  PEL*             restrict out,
+                                  REGION*          restrict in,
+                                  double                    absolute_x,
+                                  double                    absolute_y )
+{
+  /*
+   * Floor's surrogate FAST_PSEUDO_FLOOR is used to make sure that the
+   * transition through 0 is smooth. If it is known that absolute_x
+   * and absolute_y will never be less than 0, plain cast---that is,
+   * const int ix = absolute_x---should be used instead.  Actually,
+   * any function which agrees with floor for non-integer values, and
+   * picks one of the two possibilities for integer values, can be
+   * used. FAST_PSEUDO_FLOOR fits the bill.
+   *
+   * Then, x is the x-coordinate of the sampling point relative to the
+   * position of the top left corner of the convex hull of the 2x2
+   * block of closest pixels. Similarly for y. Range of values: [0,1).
+   */
+  const int ix = FAST_PSEUDO_FLOOR( absolute_x );
+  const int iy = FAST_PSEUDO_FLOOR( absolute_y );
+
+  /*
+   * Move the pointer to (the first band of) the top/left pixel of the
+   * 2x2 group of pixel centers which contains the sampling location
+   * in its convex hull:
+   */
+  const PEL* restrict p = (PEL *) IM_REGION_ADDR( in, ix, iy );
+
+  const double relative_x = absolute_x - ix;
+  const double relative_y = absolute_y - iy;
+
+  /*
+   * VIPS versions of Nicolas's pixel addressing values.
+   */
+  const int actual_bands = in->im->Bands;
+  const int lskip = IM_REGION_LSKIP( in ) / IM_IMAGE_SIZEOF_ELEMENT( in->im );
+  /*
+   * Double the bands for complex images to account for the real and
+   * imaginary parts being computed independently:
+   */
+  const int bands =
+    vips_bandfmt_iscomplex( in->im->BandFmt ) ? 2 * actual_bands : actual_bands;
+
+  switch( in->im->BandFmt ) {
+  case IM_BANDFMT_UCHAR:
+    CALL( unsigned char, nosign );
+    break;
+
+  case IM_BANDFMT_CHAR:
+    CALL( signed char, withsign );
+    break;
+
+  case IM_BANDFMT_USHORT:
+    CALL( unsigned short, nosign );
+    break;
+
+  case IM_BANDFMT_SHORT:
+    CALL( signed short, withsign );
+    break;
+
+  case IM_BANDFMT_UINT:
+    CALL( unsigned int, nosign );
+    break;
+
+  case IM_BANDFMT_INT:
+    CALL( signed int, withsign );
+    break;
+
+  /*
+   * Complex images are handled by doubling of bands.
+   */
+  case IM_BANDFMT_FLOAT:
+  case IM_BANDFMT_COMPLEX:
+    CALL( float, fptypes );
+    break;
+
+  case IM_BANDFMT_DOUBLE:
+  case IM_BANDFMT_DPCOMPLEX:
+    CALL( double, fptypes );
+    break;
+
+  default:
+    g_assert( 0 );
+    break;
+  }
+}
+
+static void
+vips_interpolate_lbb_class_init( VipsInterpolateLbbClass *klass )
+{
+  VipsObjectClass *object_class = VIPS_OBJECT_CLASS( klass );
+  VipsInterpolateClass *interpolate_class =
+    VIPS_INTERPOLATE_CLASS( klass );
+
+  object_class->nickname = "lbb";
+  object_class->description = _( "Reduced halo bicubic" );
+
+  interpolate_class->interpolate = vips_interpolate_lbb_interpolate;
+  interpolate_class->window_size = 4;
+}
+
+static void
+vips_interpolate_lbb_init( VipsInterpolateLbb *lbb )
+{
+}

libvips/resample/mbicubic.cpp

@@ -1,645 +0,0 @@
-/* symmetrized monotone cubic splines
- */
-
-/*
-
-    This file is part of VIPS.
-
-    VIPS is free software; you can redistribute it and/or modify it
-    under the terms of the GNU Lesser General Public License as
-    published by the Free Software Foundation; either version 2 of the
-    License, or (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    Lesser General Public License for more details.
-
-    You should have received a copy of the GNU Lesser General Public
-    License along with this program; if not, write to the Free
-    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
-    02111-1307 USA
-
- */
-
-/*
-
-    These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
-
- */
-
-/*
- * 2009-2010 (c) Nicolas Robidoux, John Cupitt, Eric Daoust and Adam
- * Turcotte
- *
- * E. Daoust and A. Turcotte's symmetrized monotone cubic splines
- * programming funded in part by two Google Summer of Code 2010 awards
- * made to GIMP (Gnu Image Manipulation Program) and its library GEGL.
- *
- * Nicolas Robidoux thanks Chantal Racette, Ralf Meyer, Minglun Gong,
- * Øyvind Kolås, Geert Jordaens and Sven Neumann for useful comments
- * and code.
- */
-
-/*
- * mbicubic is the VIPS name of the symmetrized implementation in 2D
- * of monotone cubic spline interpolation method a.k.a. MP-Quadratic
- * (Monotonicity Preserving with derivative estimated by fitting a
- * parabola (quadratic polynomial)) method, which essentially is
- * Catmull-Rom with derivatives clamped with Fristsh and Carlson's
- * "rule of 3" so as to ensure monotonicity.
- *
- * 1D MP-quadratic (for curve, not surface, interpolation) is
- * described in
- *
- * Accurate Monotone Cubic Interpolation, by Hung T. Huynh, published
- * in the SIAM Journal on Numerical Analysis, Volume 30, Issue 1
- * (February 1993), pages 57-100, 1993. ISSN:0036-1429.
- *
- * and in NASA technical memorandum 103789, which can be downloaded
- * from http://
- * ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19910011517_1991011517.pdf
- *
- * In order to ensure reflexion symmetry about diagonal lines, 1D
- * MP-quadratic is performed two different ways---horizontally then
- * vertically, and vertically then horizontally---and
- * averaged. (Symmetry about 45 degree lines is not automatically
- * respected because MP-quadratic is a nonlinear method: interpolating
- * horizontally then vertically does not necessarily give the same as
- * interpolating vertically then horizontally.)
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif /*HAVE_CONFIG_H*/
-#include <vips/intl.h>
-
-#include <stdio.h>
-#include <stdlib.h>
-
-#include <vips/vips.h>
-#include <vips/internal.h>
-
-#include "templates.h"
-
-#define VIPS_TYPE_INTERPOLATE_MBICUBIC \
-	(vips_interpolate_mbicubic_get_type())
-#define VIPS_INTERPOLATE_MBICUBIC( obj ) \
-	(G_TYPE_CHECK_INSTANCE_CAST( (obj), \
-	VIPS_TYPE_INTERPOLATE_MBICUBIC, VipsInterpolateMbicubic ))
-#define VIPS_INTERPOLATE_MBICUBIC_CLASS( klass ) \
-	(G_TYPE_CHECK_CLASS_CAST( (klass), \
-	VIPS_TYPE_INTERPOLATE_MBICUBIC, VipsInterpolateMbicubicClass))
-#define VIPS_IS_INTERPOLATE_MBICUBIC( obj ) \
-	(G_TYPE_CHECK_INSTANCE_TYPE( (obj), VIPS_TYPE_INTERPOLATE_MBICUBIC ))
-#define VIPS_IS_INTERPOLATE_MBICUBIC_CLASS( klass ) \
-	(G_TYPE_CHECK_CLASS_TYPE( (klass), VIPS_TYPE_INTERPOLATE_MBICUBIC ))
-#define VIPS_INTERPOLATE_MBICUBIC_GET_CLASS( obj ) \
-	(G_TYPE_INSTANCE_GET_CLASS( (obj), \
-	VIPS_TYPE_INTERPOLATE_MBICUBIC, VipsInterpolateMbicubicClass ))
-
-typedef struct _VipsInterpolateMbicubic {
-	VipsInterpolate parent_object;
-
-} VipsInterpolateMbicubic;
-
-typedef struct _VipsInterpolateMbicubicClass {
-	VipsInterpolateClass parent_class;
-
-} VipsInterpolateMbicubicClass;
-
-/*
- * MINMOD is an implementation of the minmod function which only needs
- * two conditional moves.
- *
- * MINMOD(a,b,a_times_a,a_times_b) "returns" minmod(a,b). The
- * parameter ("input") a_times_a is assumed to contain the square of
- * a; the parameter a_times_b, the product of a and b.
- *
- * The version most suitable for images with flat (constant) colour
- * areas, since a, which is a pixel difference, will often be 0, in
- * which case both forward branches are likely:
- *
- * ( (a_times_b)>=0 ? 1. : 0. ) * ( (a_times_a)<=(a_times_b) ? (a) : (b) )
- *
- * For uncompressed natural images in high bit depth (images for which
- * the slopes a and b are unlikely to be equal to zero or be equal to
- * each other), we recommend using
- *
- * ( (a_times_b)>=0. ? 1. : 0. ) * ( (a_times_b)<(a_times_a) ? (b) : (a) )
- *
- * instead. With this second version, the forward branch of the second
- * conditional move is taken when |b|>|a| and when a*b<0. However, the
- * "else" branch is taken when a=0 (or when a=b), which is why this
- * second version is not recommended for images with large regions
- * with constant pixel values (or even, actually, regions with nearby
- * pixel values which vary bilinearly, which may arise from dirt-cheap
- * demosaicing or computer graphics operations).
- *
- * Both of the above use a multiplication instead of a nested
- * "if-then-else" because gcc does not always rewrite the latter using
- * conditional moves.
- *
- * Implementation note: Both of the above are better than FAST_MINMOD
- * (currently found in templates.h and used by the "classic"
- * implementations of the other Nohalo methods). Unfortunately, MINMOD
- * uses different parameters and consequently is not a direct
- * substitute. To be fixed in the future.
- *
- * Note that the two variants differ in whether (a) or (b) is the
- * "out" of the forward branch of the second factor. If there is a
- * difference in likelihood, put the likely one in (a) in the first
- * variant, and the likely one in (b) in the second.
- */
-#define MINMOD(a,b,a_times_a,a_times_b) \
-  ( (a_times_b)>=0. ? 1. : 0. ) * ( (a_times_b)<(a_times_a) ? (b) : (a) )
-
-static inline double
-mcubic( const double one,
-        const double two,
-        const double thr,
-        const double fou,
-        const double coef_thr_point,
-        const double half_coef_two_slope,
-        const double half_coef_thr_slope )
-{
-  /*
-   * Computation of the slopes and slope limiters:
-   *
-   * Differences:
-   */
-  const double prem = two - one;
-  const double deux = thr - two;
-  const double troi = fou - thr;
-
-  const double part = two + deux * coef_thr_point;
-
-  /*
-   * Products useful for the minmod computations:
-   */
-  const double deux_prem = deux * prem;
-  const double deux_deux = deux * deux;
-  const double deux_troi = deux * troi;
-
-  /*
-   * Twice the horizontal limiter slopes (twice_lx) interwoven with
-   * twice the Catmull-Rom slopes (twice_sx).  Because we have twice
-   * the Catmull-Rom slope, we need to use 6 times the minmod slope
-   * instead of the usual 3 (specified by the cited article).
-   */
-  const double twice_lx_two =
-    6. * MINMOD( deux, prem, deux_deux, deux_prem );
-  const double twice_sx_two = deux + prem;
-  const double twice_lx_thr =
-    6. * MINMOD( deux, troi, deux_deux, deux_troi );
-  const double twice_sx_thr = deux + troi;
-
-  const double lx_lx_two = twice_lx_two * twice_lx_two;
-  const double lx_sx_two = twice_lx_two * twice_sx_two;
-  const double lx_lx_thr = twice_lx_thr * twice_lx_thr;
-  const double lx_sx_thr = twice_lx_thr * twice_sx_thr;
-
-  /*
-   * Result of the first interpolations along horizontal lines. Note
-   * that the Catmull-Rom slope almost always satisfies the
-   * monotonicity constraint, hence twice_sx is "likely" to be the one
-   * selected by minmod.
-   */
-  const double newval =
-    part +
-    + half_coef_two_slope
-    * MINMOD( twice_lx_two, twice_sx_two, lx_lx_two, lx_sx_two )
-    + half_coef_thr_slope
-    * MINMOD( twice_lx_thr, twice_sx_thr, lx_lx_thr, lx_sx_thr );
-
-  return newval;
-}
-
-static inline double
-symmetrized_monotone_cubic_splines( const double coef_rite_point,
-                                    const double coef_bot_point,
-                                    const double half_coef_left_slope,
-                                    const double half_coef_rite_slope,
-                                    const double half_coef_top_slope,
-                                    const double half_coef_bot_slope,
-                                    const double uno_one,
-                                    const double uno_two,
-                                    const double uno_thr,
-                                    const double uno_fou,
-                                    const double dos_one,
-                                    const double dos_two,
-                                    const double dos_thr,
-                                    const double dos_fou,
-                                    const double tre_one,
-                                    const double tre_two,
-                                    const double tre_thr,
-                                    const double tre_fou,
-                                    const double qua_one,
-                                    const double qua_two,
-                                    const double qua_thr,
-                                    const double qua_fou )
-{
-  /*
-   * STENCIL (FOOTPRINT) OF INPUT VALUES:
-   *
-   * The stencil of Symmetrized Monotone Catmull-Rom is the same as
-   * the standard Catmull-Rom's:
-   *
-   *  (ix-1,iy-1)  (ix,iy-1)    (ix+1,iy-1)  (ix+2,iy-1)
-   *  = uno_one    = uno_two    = uno_thr    = uno_fou
-   *
-   *  (ix-1,iy)    (ix,iy)      (ix+1,iy)    (ix+2,iy)
-   *  = dos_one    = dos_two    = dos_thr    = dos_fou
-   *                        X
-   *  (ix-1,iy+1)  (ix,iy+1)    (ix+1,iy+1)  (ix+2,iy+1)
-   *  = tre_one    = tre_two    = tre_thr    = tre_fou
-   *
-   *  (ix-1,iy+2)  (ix,iy+2)    (ix+1,iy+2)  (ix+2,iy+2)
-   *  = qua_one    = qua_two    = qua_thr    = qua_fou
-   *
-   * where ix is the (pseudo-)floor of the requested left-to-right
-   * location ("X"), and iy is the floor of the requested up-to-down
-   * location.
-   */
-  /*
-   * Outline of the computation:
-   *
-   * First, four horizontal cubic Hermite interpolations are performed
-   * to get values on the vertical line which passes through X, and
-   * then these four values are used to perform cubic Hermite
-   * interpolation in the vertical direction to get one approximation
-   * of the pixel value at X,
-   *
-   * Then, four vertical cubic Hermite interpolations are performed to
-   * get values on the horizontal line which passes through X, and
-   * then these four values are used to perform cubic Hermite
-   * interpolation in the horizontal direction to get another
-   * approximation of the pixel value at X,
-   *
-   * These two interpolated pixel values are then averaged.
-   */
-
-  /*
-   * Computation of the slopes and slope limiters:
-   *
-   * Uno horizontal differences:
-   */
-  const double uno = mcubic( uno_one,
-                             uno_two,
-                             uno_thr,
-                             uno_fou,
-                             coef_rite_point,
-                             half_coef_left_slope,
-                             half_coef_rite_slope );
-  /*
-   * Do the same with the other three horizontal lines.
-   *
-   * Dos horizontal line:
-   */
-  const double dos = mcubic( dos_one,
-                             dos_two,
-                             dos_thr,
-                             dos_fou,
-                             coef_rite_point,
-                             half_coef_left_slope,
-                             half_coef_rite_slope );
-  /*
-   * Tre(s) horizontal line:
-   */
-  const double tre = mcubic( tre_one,
-                             tre_two,
-                             tre_thr,
-                             tre_fou,
-                             coef_rite_point,
-                             half_coef_left_slope,
-                             half_coef_rite_slope );
-  /*
-   * Qua(ttro) horizontal line:
-   */
-  const double qua = mcubic( qua_one,
-                             qua_two,
-                             qua_thr,
-                             qua_fou,
-                             coef_rite_point,
-                             half_coef_left_slope,
-                             half_coef_rite_slope );
-
-  /*
-   * Perform the interpolation along the one vertical line (filled
-   * with results obtained by interpolating along horizontal lines):
-   */
-  const double partial_y = mcubic( uno,
-                                   dos,
-                                   tre,
-                                   qua,
-                                   coef_bot_point,
-                                   half_coef_top_slope,
-                                   half_coef_bot_slope );
-
-  /*
-   * Redo with four vertical lines (and the corresponding horizontal
-   * one).
-   *
-   * One:
-   */
-  const double one = mcubic( uno_one,
-                             dos_one,
-                             tre_one,
-                             qua_one,
-                             coef_bot_point,
-                             half_coef_top_slope,
-                             half_coef_bot_slope );
-  /*
-   * Two:
-   */
-  const double two = mcubic( uno_two,
-                             dos_two,
-                             tre_two,
-                             qua_two,
-                             coef_bot_point,
-                             half_coef_top_slope,
-                             half_coef_bot_slope );
-  /*
-   * Thr(ee):
-   */
-  const double thr = mcubic( uno_thr,
-                             dos_thr,
-                             tre_thr,
-                             qua_thr,
-                             coef_bot_point,
-                             half_coef_top_slope,
-                             half_coef_bot_slope );
-  /*
-   * Fou(r):
-   */
-  const double fou = mcubic( uno_fou,
-                             dos_fou,
-                             tre_fou,
-                             qua_fou,
-                             coef_bot_point,
-                             half_coef_top_slope,
-                             half_coef_bot_slope );
-
-  /*
-   * Final horizontal line of vertical results:
-   */
-  const double prem_x = two - one;
-  const double deux_x = thr - two;
-  const double troi_x = fou - thr;
-
-  const double partial_newval = partial_y + two + coef_rite_point * deux_x;
-
-  const double deux_prem_x = deux_x * prem_x;
-  const double deux_deux_x = deux_x * deux_x;
-  const double deux_troi_x = deux_x * troi_x;
-
-  const double twice_l_two =
-    6. * MINMOD( deux_x, prem_x, deux_deux_x, deux_prem_x );
-  const double twice_s_two = deux_x + prem_x;
-  const double twice_l_thr =
-    6. * MINMOD( deux_x, troi_x, deux_deux_x, deux_troi_x );
-  const double twice_s_thr = deux_x + troi_x;
-
-  const double l_l_two = twice_l_two * twice_l_two;
-  const double l_s_two = twice_l_two * twice_s_two;
-  const double l_l_thr = twice_l_thr * twice_l_thr;
-  const double l_s_thr = twice_l_thr * twice_s_thr;
-
-  const double newval =
-    (
-      partial_newval
-      +
-      half_coef_left_slope
-      *
-      MINMOD( twice_l_two, twice_s_two, l_l_two, l_s_two )
-      +
-      half_coef_rite_slope
-      *
-      MINMOD( twice_l_thr, twice_s_thr, l_l_thr, l_s_thr )
-    ) * .5;
-
-  return newval;
-}
-
-/*
- * Call Snohalo with an conversion operator as a parameter.
- *
- * It would be nice to do this with templates somehow---for one thing
- * this would allow code comments!---but we can't figure a clean way
- * to do it.
- */
-#define MBICUBIC_CONVERSION( conversion )               \
-  template <typename T> static void inline              \
-  mbicubic_ ## conversion(       PEL*   restrict pout,  \
-                           const PEL*   restrict pin,   \
-                           const int             bands, \
-                           const int             lskip, \
-                           const double          x,     \
-                           const double          y )    \
-  { \
-    T* restrict out = (T *) pout; \
-    \
-    const T* restrict in = (T *) pin; \
-    \
-    const int one_shift =  -bands; \
-    const int thr_shift =   bands; \
-    const int fou_shift = 2*bands; \
-    \
-    const int uno_two_shift =  -lskip; \
-    const int dos_two_shift =       0; \
-    const int tre_two_shift =   lskip; \
-    const int qua_two_shift = 2*lskip; \
-    \
-    const int uno_one_shift = uno_two_shift + one_shift; \
-    const int dos_one_shift = dos_two_shift + one_shift; \
-    const int tre_one_shift = tre_two_shift + one_shift; \
-    const int qua_one_shift = qua_two_shift + one_shift; \
-    \
-    const int uno_thr_shift = uno_two_shift + thr_shift; \
-    const int dos_thr_shift = dos_two_shift + thr_shift; \
-    const int tre_thr_shift = tre_two_shift + thr_shift; \
-    const int qua_thr_shift = qua_two_shift + thr_shift; \
-    \
-    const int uno_fou_shift = uno_two_shift + fou_shift; \
-    const int dos_fou_shift = dos_two_shift + fou_shift; \
-    const int tre_fou_shift = tre_two_shift + fou_shift; \
-    const int qua_fou_shift = qua_two_shift + fou_shift; \
-    \
-    const double x_squared = x * x; \
-    const double y_squared = y * y; \
-    const double twice_x   = x + x; \
-    const double twice_y   = y + y; \
-    const double half_x_squared_minus_x = .5 * ( x_squared - x ); \
-    const double half_y_squared_minus_y = .5 * ( y_squared - y ); \
-    \
-    const double coef_rite_point = x_squared * ( 3. - twice_x ); \
-    const double coef_bot_point  = y_squared * ( 3. - twice_y ); \
-    \
-    const double half_coef_rite_slope = x * half_x_squared_minus_x; \
-    const double half_coef_bot_slope  = y * half_y_squared_minus_y; \
-    const double half_coef_left_slope = \
-      half_coef_rite_slope - half_x_squared_minus_x; \
-    const double half_coef_top_slope  = \
-      half_coef_bot_slope  - half_y_squared_minus_y; \
-    \
-    int band = bands; \
-    \
-    do \
-      { \
-        const double double_result = \
-          symmetrized_monotone_cubic_splines( coef_rite_point,       \
-                                              coef_bot_point,        \
-                                              half_coef_left_slope,  \
-                                              half_coef_rite_slope,  \
-                                              half_coef_top_slope,   \
-                                              half_coef_bot_slope,   \
-                                              in[ uno_one_shift ],   \
-                                              in[ uno_two_shift ],   \
-                                              in[ uno_thr_shift ],   \
-                                              in[ uno_fou_shift ],   \
-                                              in[ dos_one_shift ],   \
-                                              in[ dos_two_shift ],   \
-                                              in[ dos_thr_shift ],   \
-                                              in[ dos_fou_shift ],   \
-                                              in[ tre_one_shift ],   \
-                                              in[ tre_two_shift ],   \
-                                              in[ tre_thr_shift ],   \
-                                              in[ tre_fou_shift ],   \
-                                              in[ qua_one_shift ],   \
-                                              in[ qua_two_shift ],   \
-                                              in[ qua_thr_shift ],   \
-                                              in[ qua_fou_shift ] ); \
-        \
-        const T result = to_ ## conversion<T>( double_result ); \
-        in++; \
-        *out++ = result; \
-      } while (--band); \
-  }
-
-MBICUBIC_CONVERSION( fptypes )
-MBICUBIC_CONVERSION( withsign )
-MBICUBIC_CONVERSION( nosign )
-
-#define CALL( T, conversion )               \
-  mbicubic_ ## conversion<T>( out,          \
-                              p,            \
-                              bands,        \
-                              lskip,        \
-                              relative_x,   \
-                              relative_y );
-
-/*
- * We need C linkage:
- */
-extern "C" {
-G_DEFINE_TYPE( VipsInterpolateMbicubic, vips_interpolate_mbicubic,
-	VIPS_TYPE_INTERPOLATE );
-}
-
-static void
-vips_interpolate_mbicubic_interpolate( VipsInterpolate* restrict interpolate,
-                                       PEL*             restrict out,
-                                       REGION*          restrict in,
-                                       double                    absolute_x,
-                                       double                    absolute_y )
-{
-  /*
-   * Floor's surrogate FAST_PSEUDO_FLOOR is used to make sure that the
-   * transition through 0 is smooth. If it is known that absolute_x
-   * and absolute_y will never be less than 0, plain cast---that is,
-   * const int ix = absolute_x---should be used instead.  Actually,
-   * any function which agrees with floor for non-integer values, and
-   * picks one of the two possibilities for integer values, can be
-   * used. FAST_PSEUDO_FLOOR fits the bill.
-   *
-   * Then, x is the x-coordinate of the sampling point relative to the
-   * position of the center of the convex hull of the 2x2 block of
-   * closest pixels. Similarly for y. Range of values: [-.5,.5).
-   */
-  const int ix = FAST_PSEUDO_FLOOR( absolute_x );
-  const int iy = FAST_PSEUDO_FLOOR( absolute_y );
-
-  /*
-   * Move the pointer to (the first band of) the top/left pixel of the
-   * 2x2 group of pixel centers which contains the sampling location
-   * in its convex hull:
-   */
-  const PEL* restrict p = (PEL *) IM_REGION_ADDR( in, ix, iy );
-
-  const double relative_x = absolute_x - ix;
-  const double relative_y = absolute_y - iy;
-
-  /*
-   * VIPS versions of Nicolas's pixel addressing values.
-   */
-  const int actual_bands = in->im->Bands;
-  const int lskip = IM_REGION_LSKIP( in ) / IM_IMAGE_SIZEOF_ELEMENT( in->im );
-  /*
-   * Double the bands for complex images to account for the real and
-   * imaginary parts being computed independently:
-   */
-  const int bands =
-    vips_bandfmt_iscomplex( in->im->BandFmt ) ? 2 * actual_bands : actual_bands;
-
-  switch( in->im->BandFmt ) {
-  case IM_BANDFMT_UCHAR:
-    CALL( unsigned char, nosign );
-    break;
-
-  case IM_BANDFMT_CHAR:
-    CALL( signed char, withsign );
-    break;
-
-  case IM_BANDFMT_USHORT:
-    CALL( unsigned short, nosign );
-    break;
-
-  case IM_BANDFMT_SHORT:
-    CALL( signed short, withsign );
-    break;
-
-  case IM_BANDFMT_UINT:
-    CALL( unsigned int, nosign );
-    break;
-
-  case IM_BANDFMT_INT:
-    CALL( signed int, withsign );
-    break;
-
-  /*
-   * Complex images are handled by doubling of bands.
-   */
-  case IM_BANDFMT_FLOAT:
-  case IM_BANDFMT_COMPLEX:
-    CALL( float, fptypes );
-    break;
-
-  case IM_BANDFMT_DOUBLE:
-  case IM_BANDFMT_DPCOMPLEX:
-    CALL( double, fptypes );
-    break;
-
-  default:
-    g_assert( 0 );
-    break;
-  }
-}
-
-static void
-vips_interpolate_mbicubic_class_init( VipsInterpolateMbicubicClass *klass )
-{
-  VipsObjectClass *object_class = VIPS_OBJECT_CLASS( klass );
-  VipsInterpolateClass *interpolate_class =
-    VIPS_INTERPOLATE_CLASS( klass );
-
-  object_class->nickname = "mbicubic";
-  object_class->description = _( "Halo-free mbicubic" );
-
-  interpolate_class->interpolate = vips_interpolate_mbicubic_interpolate;
-  interpolate_class->window_size = 4;
-}
-
-static void
-vips_interpolate_mbicubic_init( VipsInterpolateMbicubic *mbicubic )
-{
-}