/* nohalo subdivision followed by lbb (locally bounded bicubic)
 * interpolation resampler
 *
 * Nohalo level 1 with bilinear finishing scheme hacked for VIPS by
 * J. Cupitt based on code by N. Robidoux, 20/1/09
 *
 * N. Robidoux and J. Cupitt, 4-17/3/09
 *
 * N. Robidoux, 1/4-29/5/2009
 *
 * Nohalo level 2 with bilinear finishing scheme by N. Robidoux based
 * on code by N. Robidoux, A. Turcotte and J. Cupitt, 27/1/2010
 *
 * Nohalo level 1 with LBB finishing scheme by N. Robidoux and
 * C. Racette, 11-18/5/2010
 */

/*

    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., 51 Franklin Street, Fifth Floor, Boston, MA
    02110-1301 USA

 */

/*

    These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk

 */

/*
 * 2009-2010 (c) Nicolas Robidoux, Chantal Racette, John Cupitt and
 * Adam Turcotte
 *
 * N. Robidoux thanks Geert Jordaens, Ralf Meyer, Øyvind Kolås,
 * Minglun Gong, Eric Daoust and Sven Neumann for useful comments and
 * code.
 *
 * N. Robidoux's early research on Nohalo funded in part by an NSERC
 * (National Science and Engineering Research Council of Canada)
 * Discovery Grant awarded to him (298424--2004).
 *
 * C. Racette's image resampling research and programming funded in
 * part by an NSERC (National Science and Engineering Research Council
 * of Canada) Alexander Graham Bell Canada Graduate Scholarship, by an
 * NSERC Discovery Grant awarded to Julien Dompierre (grant number
 * 20-61098) and by N. Robidoux's Laurentian University professional
 * allowance.
 *
 * A. Turcotte's image resampling research on reduced halo funded in
 * part by an NSERC Alexander Graham Bell Canada Graduate Scholarhip
 * awarded to him and by a Google Summer of Code 2010 award awarded to
 * GIMP (Gnu Image Manipulation Program).
 *
 * Nohalo with LBB finishing scheme was developed by N. Robidoux and
 * C. Racette at the Department of Mathematics and Computer Science of
 * Laurentian University in the course of C. Racette's Masters thesis
 * in Computational Sciences. Preliminary work on Nohalo and monotone
 * interpolation was performed by C. Racette and N. Robidoux in the
 * course of her honours thesis, by N. Robidoux, A. Turcotte and
 * E. Daoust during Google Summer of Code 2009 (through two awards
 * made to GIMP to improve GEGL), and, earlier, by N. Robidoux,
 * A. Turcotte, J. Cupitt, M. Gong and K. Martinez.
 */

/*
 * Nohalo with LBB as finishing scheme has two versions, which are
 * only different in the way LBB is implemented:
 *
 *   A "soft" version, which shows a little less staircasing and a
 *   little more haloing, and which is a little more expensive to
 *   compute. We recommend this as the default.
 *
 *   A "sharp" version, which shows a little more staircasing and a
 *   little less haloing, and which is a little cheaper (it uses 6
 *   less comparisons and 12 less "? :").
 *
 * The only difference between the two is that the "soft" versions
 * uses local minima and maxima computed over 3x3 square blocks, and
 * the "sharp" version uses local minima and maxima computed over 3x3
 * crosses.
 *
 * The "sharp" version is (a little) faster. We don't know yet for
 * sure, but it appears that the "soft" version gives marginally
 * better results.
 *
 * If you want to use the "sharp" (cheaper) version, uncomment the
 * following three pre-processor code lines:
 */

/*
#ifndef __NOHALO_CHEAP_H__
#define __NOHALO_CHEAP_H__
#endif
 */

/*
 * ================
 * NOHALO RESAMPLER
 * ================
 *
 * "Nohalo" is a resampler with a mission: smoothly straightening
 * oblique lines without undesirable side-effects. In particular,
 * without much blurring and with no added haloing.
 *
 * In this code, one Nohalo subdivision is performed. The
 * interpolation is finished with LBB (Locally Bounded Bicubic).
 *
 * Key properties:
 *
 * =======================
 * Nohalo is interpolatory
 * =======================
 *
 * That is, Nohalo preserves point values: If asked for the value at
 * the center of an input pixel, the sampler returns the corresponding
 * value, unchanged. In addition, because Nohalo is continuous, if
 * asked for a value at a location "very close" to the center of an
 * input pixel, then the sampler returns a value "very close" to
 * it. (Nohalo is not smoothing like, say, B-Spline
 * pseudo-interpolation.)
 *
 * ====================================================================
 * Nohalo subdivision is co-monotone (this is why it's called "no-halo")
 * ====================================================================
 *
 * One consequence of monotonicity is that additional subdivided
 * values are in the range of the four closest input values, which is
 * a form of local boundedness.  (Note: plain vanilla bilinear and
 * nearest neighbour are also co-monotone.) LBB is also locally
 * bounded. Consequently, Nohalo subdivision followed by LBB is
 * locally bounded. When used as a finishing scheme for Nohalo, the
 * standard LBB bounds imply that the final interpolated value is in
 * the range of the nine closest input values. This property is why
 * there is very little added haloing, even when a finishing scheme
 * which is not strictly monotone. Another consequence of local
 * boundedness is that clamping is unnecessary (provided abyss values
 * are within the range of acceptable values, which is "always" the
 * case).
 *
 * Note: If the abyss policy is an extrapolating one---for example,
 * linear or bilinear extrapolation---clamping is still unnecessary
 * UNLESS one attempts to resample outside of the convex hull of the
 * input pixel positions. Consequence: the "corner" image size
 * convention does not require clamping when using linear
 * extrapolation abyss policy when performing image resizing, but the
 * "center" one does, when upscaling, at locations very close to the
 * boundary. If computing values at locations outside of the convex
 * hull of the pixel locations of the input image, nearest neighbour
 * abyss policy is most likely better anyway, because linear
 * extrapolation produces "streaks" if positions far outside the
 * original image boundary are resampled.
 *
 * ========================
 * Nohalo is a local method
 * ========================
 *
 * The interpolated pixel value when using Nohalo subdivision followed
 * by LBB only depends on the 21 (5x5 minus the four corners) closest
 * input values.
 *
 * ===============================
 * Nohalo is second order accurate
 * ===============================
 *
 * (Except possibly near the boundary: it is easy to make this
 * property carry over everywhere but this requires a tuned abyss
 * policy---linear extrapolation, say---or building the boundary
 * conditions inside the sampler.)  Nohalo+LBB is exact on linear
 * intensity profiles, meaning that if the input pixel values (in the
 * stencil) are obtained from a function of the form f(x,y) = a + b*x
 * + c*y (a, b, c constants), then the computed pixel value is exactly
 * the value of f(x,y) at the asked-for sampling location. The
 * boundary condition which is emulated by VIPS through the "extend"
 * extension of the input image---this corresponds to the nearest
 * neighbour abyss policy---does NOT make this resampler exact on
 * linears near the boundary. It does, however, guarantee that no
 * clamping is required even when resampled values are computed at
 * positions outside of the extent of the input image (when
 * extrapolation is required).
 *
 * ===================
 * Nohalo is nonlinear
 * ===================
 *
 * Both Nohalo and LBB are nonlinear, consequently their composition
 * is nonlinear.  In particular, resampling a sum of images may not be
 * the same as summing the resamples. (This occurs even without taking
 * into account over and underflow issues: images can only take values
 * within a banded range, and consequently no sampler is truly
 * linear.)
 *
 * ====================
 * Weaknesses of Nohalo
 * ====================
 *
 * In some cases, the initial subdivision computation is wasted:
 *
 * If a region is bi-chromatic, the nonlinear component of Nohalo
 * subdivision is zero in the interior of the region, and consequently
 * Nohalo subdivision boils down to bilinear. For such images, LBB is
 * probably a better choice.
 *
 * =========================
 * Bibliographical reference
 * =========================
 *
 * For more information about Nohalo (a prototype version with
 * bilinear finish instead of LBB), see
 *
 * CPU, SMP and GPU implementations of Nohalo level 1, a fast
 * co-convex antialiasing image resampler by Nicolas Robidoux, Minglun
 * Gong, John Cupitt, Adam Turcotte, and Kirk Martinez, in C3S2E '09:
 * Proceedings of the 2nd Canadian Conference on Computer Science and
 * Software Engineering, p. 185--195, ACM, New York, NY, USA, 2009.
 * http://doi.acm.org/10.1145/1557626.1557657.
 */

/* Uncomment to enable bounds checking for VIPS_REGION_ADDR().
 */
#define DEBUG

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /*HAVE_CONFIG_H*/
#include <glib/gi18n-lib.h>

#include <stdio.h>
#include <stdlib.h>

#include <vips/vips.h>
#include <vips/internal.h>

#include "templates.h"

#define VIPS_TYPE_INTERPOLATE_NOHALO \
	(vips_interpolate_nohalo_get_type())
#define VIPS_INTERPOLATE_NOHALO( obj ) \
	(G_TYPE_CHECK_INSTANCE_CAST( (obj), \
	VIPS_TYPE_INTERPOLATE_NOHALO, VipsInterpolateNohalo ))
#define VIPS_INTERPOLATE_NOHALO_CLASS( klass ) \
	(G_TYPE_CHECK_CLASS_CAST( (klass), \
	VIPS_TYPE_INTERPOLATE_NOHALO, VipsInterpolateNohaloClass))
#define VIPS_IS_INTERPOLATE_NOHALO( obj ) \
	(G_TYPE_CHECK_INSTANCE_TYPE( (obj), VIPS_TYPE_INTERPOLATE_NOHALO ))
#define VIPS_IS_INTERPOLATE_NOHALO_CLASS( klass ) \
	(G_TYPE_CHECK_CLASS_TYPE( (klass), VIPS_TYPE_INTERPOLATE_NOHALO ))
#define VIPS_INTERPOLATE_NOHALO_GET_CLASS( obj ) \
	(G_TYPE_INSTANCE_GET_CLASS( (obj), \
	VIPS_TYPE_INTERPOLATE_NOHALO, VipsInterpolateNohaloClass ))

typedef struct _VipsInterpolateNohalo {
	VipsInterpolate parent_object;

} VipsInterpolateNohalo;

typedef struct _VipsInterpolateNohaloClass {
	VipsInterpolateClass parent_class;

} VipsInterpolateNohaloClass;

/*
 * NOHALO_MINMOD is an implementation of the minmod function which
 * only needs two "conditional moves."
 * NOHALO_MINMOD(a,b,a_times_a,a_times_b) "returns"
 * minmod(a,b). The macro parameter ("input") a_times_a is assumed to
 * contain the square of a; a_times_b, the product of a and 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), or chips with good branch prediction, the following
 * version of the minmod function may work well:
 *
 * ( (a_times_b)>=0. ? ( (a_times_b)<(a_times_a) ? (b) : (a) ) : 0. )
 *
 * In this 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 the above version is not
 * as effective for images with regions with constant pixel values (or
 * regions with pixel values which vary linearly or bilinearly) since
 * we apply minmod to pairs of differences.
 *
 * The following version is more 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. This
 * may be preferable if "branch flag look ahead" does not work so
 * well.
 *
 * ( (a_times_b)>=0. ? ( (a_times_a)<=(a_times_b) ? (a) : (b) ) : 0. )
 *
 * This last version appears to be slightly better than the former in
 * speed tests performed on a recent multicore Intel chip, especially
 * when enlarging a sharp image by a large factor, hence the choice.
 */
#define NOHALO_MINMOD(a,b,a_times_a,a_times_b) \
  ( ( (a_times_b)>=0. ) ? ( (a_times_a)<=(a_times_b) ? (a) : (b) ) : 0. )

/*
 * Absolute value and sign macros:
 */
#define NOHALO_ABS(x)  ( ((x)>=0.) ? (x) : -(x) )
#define NOHALO_SIGN(x) ( ((x)>=0.) ? 1.  : -1.  )

/*
 * MIN and MAX macros set up so that I can put the likely winner in
 * the first argument (forward branch likely blah blah blah):
 */
#define NOHALO_MIN(x,y) ( ((x)<=(y)) ? (x) : (y) )
#define NOHALO_MAX(x,y) ( ((x)>=(y)) ? (x) : (y) )


static void inline
nohalo_subdivision (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           dos_fiv,
                    const double           tre_one,
                    const double           tre_two,
                    const double           tre_thr,
                    const double           tre_fou,
                    const double           tre_fiv,
                    const double           qua_one,
                    const double           qua_two,
                    const double           qua_thr,
                    const double           qua_fou,
                    const double           qua_fiv,
                    const double           cin_two,
                    const double           cin_thr,
                    const double           cin_fou,
                          double* restrict uno_one_1,
                          double* restrict uno_two_1,
                          double* restrict uno_thr_1,
                          double* restrict uno_fou_1,
                          double* restrict dos_one_1,
                          double* restrict dos_two_1,
                          double* restrict dos_thr_1,
                          double* restrict dos_fou_1,
                          double* restrict tre_one_1,
                          double* restrict tre_two_1,
                          double* restrict tre_thr_1,
                          double* restrict tre_fou_1,
                          double* restrict qua_one_1,
                          double* restrict qua_two_1,
                          double* restrict qua_thr_1,
                          double* restrict qua_fou_1)
{
  /*
   * nohalo_subdivision calculates the missing twelve double density
   * pixel values, and also returns the "already known" four, so that
   * the sixteen values which make up the stencil of LBB are
   * available.
   */
  /*
   * THE STENCIL OF INPUT VALUES:
   *
   * Pointer arithmetic is used to implicitly reflect the input
   * stencil about tre_thr---assumed closer to the sampling location
   * than other pixels (ties are OK)---in such a way that after
   * reflection the sampling point is to the bottom right of tre_thr.
   *
   * The following code and picture assumes that the stencil reflexion
   * has already been performed.
   *
   *               (ix-1,iy-2)  (ix,iy-2)    (ix+1,iy-2)
   *               =uno_two     = uno_thr    = uno_fou
   *
   *
   *
   *  (ix-2,iy-1)  (ix-1,iy-1)  (ix,iy-1)    (ix+1,iy-1)  (ix+2,iy-1)
   *  = dos_one    = dos_two    = dos_thr    = dos_fou    = dos_fiv
   *
   *
   *
   *  (ix-2,iy)    (ix-1,iy)    (ix,iy)      (ix+1,iy)    (ix+2,iy)
   *  = tre_one    = tre_two    = tre_thr    = tre_fou    = tre_fiv
   *                                    X
   *
   *
   *  (ix-2,iy+1)  (ix-1,iy+1)  (ix,iy+1)    (ix+1,iy+1)  (ix+2,iy+1)
   *  = qua_one    = qua_two    = qua_thr    = qua_fou    = qua_fiv
   *
   *
   *
   *               (ix-1,iy+2)  (ix,iy+2)    (ix+1,iy+2)
   *               = cin_two    = cin_thr    = cin_fou
   *
   *
   * The above input pixel values are the ones needed in order to make
   * available the following values, needed by LBB:
   *
   *  uno_one_1 =      uno_two_1 =  uno_thr_1 =      uno_fou_1 =
   *  (ix-1/2,iy-1/2)  (ix,iy-1/2)  (ix+1/2,iy-1/2)  (ix+1,iy-1/2)
   *
   *
   *
   *
   *  dos_one_1 =      dos_two_1 =  dos_thr_1 =      dos_fou_1 =
   *  (ix-1/2,iy)      (ix,iy)      (ix+1/2,iy)      (ix+1,iy)
   *
   *                             X
   *
   *
   *  tre_one_1 =      tre_two_1 =  tre_thr_1 =      tre_fou_1 =
   *  (ix-1/2,iy+1/2)  (ix,iy+1/2)  (ix+1/2,iy+1/2)  (ix+1,iy+1/2)
   *
   *
   *
   *
   *  qua_one_1 =      qua_two_1 =  qua_thr_1 =      qua_fou_1 =
   *  (ix-1/2,iy+1)    (ix,iy+1)    (ix+1/2,iy+1)    (ix+1,iy+1)
   *
   */

  /*
   * Computation of the nonlinear slopes: If two consecutive pixel
   * value differences have the same sign, the smallest one (in
   * absolute value) is taken to be the corresponding slope; if the
   * two consecutive pixel value differences don't have the same sign,
   * the corresponding slope is set to 0.
   *
   * In other words: Apply minmod to consecutive differences.
   */
  /*
   * Two vertical simple differences:
   */
  const double d_unodos_two = dos_two - uno_two;
  const double d_dostre_two = tre_two - dos_two;
  const double d_trequa_two = qua_two - tre_two;
  const double d_quacin_two = cin_two - qua_two;
  /*
   * Thr(ee) vertical differences:
   */
  const double d_unodos_thr = dos_thr - uno_thr;
  const double d_dostre_thr = tre_thr - dos_thr;
  const double d_trequa_thr = qua_thr - tre_thr;
  const double d_quacin_thr = cin_thr - qua_thr;
  /*
   * Fou(r) vertical differences:
   */
  const double d_unodos_fou = dos_fou - uno_fou;
  const double d_dostre_fou = tre_fou - dos_fou;
  const double d_trequa_fou = qua_fou - tre_fou;
  const double d_quacin_fou = cin_fou - qua_fou;
  /*
   * Dos horizontal differences:
   */
  const double d_dos_onetwo = dos_two - dos_one;
  const double d_dos_twothr = dos_thr - dos_two;
  const double d_dos_thrfou = dos_fou - dos_thr;
  const double d_dos_foufiv = dos_fiv - dos_fou;
  /*
   * Tre(s) horizontal differences:
   */
  const double d_tre_onetwo = tre_two - tre_one;
  const double d_tre_twothr = tre_thr - tre_two;
  const double d_tre_thrfou = tre_fou - tre_thr;
  const double d_tre_foufiv = tre_fiv - tre_fou;
  /*
   * Qua(ttro) horizontal differences:
   */
  const double d_qua_onetwo = qua_two - qua_one;
  const double d_qua_twothr = qua_thr - qua_two;
  const double d_qua_thrfou = qua_fou - qua_thr;
  const double d_qua_foufiv = qua_fiv - qua_fou;

  /*
   * Recyclable vertical products and squares:
   */
  const double d_unodos_times_dostre_two = d_unodos_two * d_dostre_two;
  const double d_dostre_two_sq           = d_dostre_two * d_dostre_two;
  const double d_dostre_times_trequa_two = d_dostre_two * d_trequa_two;
  const double d_trequa_times_quacin_two = d_quacin_two * d_trequa_two;
  const double d_quacin_two_sq           = d_quacin_two * d_quacin_two;

  const double d_unodos_times_dostre_thr = d_unodos_thr * d_dostre_thr;
  const double d_dostre_thr_sq           = d_dostre_thr * d_dostre_thr;
  const double d_dostre_times_trequa_thr = d_trequa_thr * d_dostre_thr;
  const double d_trequa_times_quacin_thr = d_trequa_thr * d_quacin_thr;
  const double d_quacin_thr_sq           = d_quacin_thr * d_quacin_thr;

  const double d_unodos_times_dostre_fou = d_unodos_fou * d_dostre_fou;
  const double d_dostre_fou_sq           = d_dostre_fou * d_dostre_fou;
  const double d_dostre_times_trequa_fou = d_trequa_fou * d_dostre_fou;
  const double d_trequa_times_quacin_fou = d_trequa_fou * d_quacin_fou;
  const double d_quacin_fou_sq           = d_quacin_fou * d_quacin_fou;
  /*
   * Recyclable horizontal products and squares:
   */
  const double d_dos_onetwo_times_twothr = d_dos_onetwo * d_dos_twothr;
  const double d_dos_twothr_sq           = d_dos_twothr * d_dos_twothr;
  const double d_dos_twothr_times_thrfou = d_dos_twothr * d_dos_thrfou;
  const double d_dos_thrfou_times_foufiv = d_dos_thrfou * d_dos_foufiv;
  const double d_dos_foufiv_sq           = d_dos_foufiv * d_dos_foufiv;

  const double d_tre_onetwo_times_twothr = d_tre_onetwo * d_tre_twothr;
  const double d_tre_twothr_sq           = d_tre_twothr * d_tre_twothr;
  const double d_tre_twothr_times_thrfou = d_tre_thrfou * d_tre_twothr;
  const double d_tre_thrfou_times_foufiv = d_tre_thrfou * d_tre_foufiv;
  const double d_tre_foufiv_sq           = d_tre_foufiv * d_tre_foufiv;

  const double d_qua_onetwo_times_twothr = d_qua_onetwo * d_qua_twothr;
  const double d_qua_twothr_sq           = d_qua_twothr * d_qua_twothr;
  const double d_qua_twothr_times_thrfou = d_qua_thrfou * d_qua_twothr;
  const double d_qua_thrfou_times_foufiv = d_qua_thrfou * d_qua_foufiv;
  const double d_qua_foufiv_sq           = d_qua_foufiv * d_qua_foufiv;

  /*
   * Minmod slopes and first level pixel values:
   */
  const double dos_thr_y = NOHALO_MINMOD( d_dostre_thr, d_unodos_thr,
                                          d_dostre_thr_sq,
                                          d_unodos_times_dostre_thr );
  const double tre_thr_y = NOHALO_MINMOD( d_dostre_thr, d_trequa_thr,
                                          d_dostre_thr_sq,
                                          d_dostre_times_trequa_thr );

  const double newval_uno_two =
    .5 * ( dos_thr + tre_thr )
    +
    .25 * ( dos_thr_y - tre_thr_y );

  const double qua_thr_y = NOHALO_MINMOD( d_quacin_thr, d_trequa_thr,
                                          d_quacin_thr_sq,
                                          d_trequa_times_quacin_thr );

  const double newval_tre_two =
    .5 * ( tre_thr + qua_thr )
    +
    .25 * ( tre_thr_y - qua_thr_y );

  const double tre_fou_y = NOHALO_MINMOD( d_dostre_fou, d_trequa_fou,
                                          d_dostre_fou_sq,
                                          d_dostre_times_trequa_fou );
  const double qua_fou_y = NOHALO_MINMOD( d_quacin_fou, d_trequa_fou,
                                          d_quacin_fou_sq,
                                          d_trequa_times_quacin_fou );

  const double newval_tre_fou =
    .5 * ( tre_fou + qua_fou )
    +
    .25 * ( tre_fou_y - qua_fou_y );

  const double dos_fou_y = NOHALO_MINMOD( d_dostre_fou, d_unodos_fou,
                                          d_dostre_fou_sq,
                                          d_unodos_times_dostre_fou );

  const double newval_uno_fou =
     .5 * ( dos_fou + tre_fou )
     +
     .25 * (dos_fou_y - tre_fou_y );

  const double tre_two_x = NOHALO_MINMOD( d_tre_twothr, d_tre_onetwo,
                                          d_tre_twothr_sq,
                                          d_tre_onetwo_times_twothr );
  const double tre_thr_x = NOHALO_MINMOD( d_tre_twothr, d_tre_thrfou,
                                          d_tre_twothr_sq,
                                          d_tre_twothr_times_thrfou );

  const double newval_dos_one =
    .5 * ( tre_two + tre_thr )
    +
    .25 * ( tre_two_x - tre_thr_x );

  const double tre_fou_x = NOHALO_MINMOD( d_tre_foufiv, d_tre_thrfou,
                                          d_tre_foufiv_sq,
                                          d_tre_thrfou_times_foufiv );

  const double tre_thr_x_minus_tre_fou_x =
    tre_thr_x - tre_fou_x;

  const double newval_dos_thr =
    .5 * ( tre_thr + tre_fou )
    +
    .25 * tre_thr_x_minus_tre_fou_x;

  const double qua_thr_x = NOHALO_MINMOD( d_qua_twothr, d_qua_thrfou,
                                          d_qua_twothr_sq,
                                          d_qua_twothr_times_thrfou );
  const double qua_fou_x = NOHALO_MINMOD( d_qua_foufiv, d_qua_thrfou,
                                          d_qua_foufiv_sq,
                                          d_qua_thrfou_times_foufiv );

  const double qua_thr_x_minus_qua_fou_x =
    qua_thr_x - qua_fou_x;

  const double newval_qua_thr =
    .5 * ( qua_thr + qua_fou )
    +
    .25 * qua_thr_x_minus_qua_fou_x;

  const double qua_two_x = NOHALO_MINMOD( d_qua_twothr, d_qua_onetwo,
                                          d_qua_twothr_sq,
                                          d_qua_onetwo_times_twothr );

  const double newval_qua_one =
    .5 * ( qua_two + qua_thr )
    +
    .25 * ( qua_two_x - qua_thr_x );

  const double newval_tre_thr =
    .125 * ( tre_thr_x_minus_tre_fou_x + qua_thr_x_minus_qua_fou_x )
    +
    .5 * ( newval_tre_two + newval_tre_fou );

  const double dos_thr_x = NOHALO_MINMOD( d_dos_twothr, d_dos_thrfou,
                                          d_dos_twothr_sq,
                                          d_dos_twothr_times_thrfou );
  const double dos_fou_x = NOHALO_MINMOD( d_dos_foufiv, d_dos_thrfou,
                                          d_dos_foufiv_sq,
                                          d_dos_thrfou_times_foufiv );

  const double newval_uno_thr =
    .25 * ( dos_fou - tre_thr )
    +
    .125 * ( dos_fou_y - tre_fou_y + dos_thr_x - dos_fou_x )
    +
    .5 * ( newval_uno_two + newval_dos_thr );

  const double tre_two_y = NOHALO_MINMOD( d_dostre_two, d_trequa_two,
                                          d_dostre_two_sq,
                                          d_dostre_times_trequa_two );
  const double qua_two_y = NOHALO_MINMOD( d_quacin_two, d_trequa_two,
                                          d_quacin_two_sq,
                                          d_trequa_times_quacin_two );

  const double newval_tre_one =
    .25 * ( qua_two - tre_thr )
    +
    .125 * ( qua_two_x - qua_thr_x + tre_two_y - qua_two_y )
    +
    .5 * ( newval_dos_one + newval_tre_two );

  const double dos_two_x = NOHALO_MINMOD( d_dos_twothr, d_dos_onetwo,
                                          d_dos_twothr_sq,
                                          d_dos_onetwo_times_twothr );

  const double dos_two_y = NOHALO_MINMOD( d_dostre_two, d_unodos_two,
                                          d_dostre_two_sq,
                                          d_unodos_times_dostre_two );

  const double newval_uno_one =
    .25 * ( dos_two + dos_thr + tre_two + tre_thr )
    +
    .125 * ( dos_two_x - dos_thr_x + tre_two_x - tre_thr_x
             +
             dos_two_y + dos_thr_y - tre_two_y - tre_thr_y );

  /*
   * Return the sixteen LBB stencil values:
   */
  *uno_one_1 = newval_uno_one;
  *uno_two_1 = newval_uno_two;
  *uno_thr_1 = newval_uno_thr;
  *uno_fou_1 = newval_uno_fou;
  *dos_one_1 = newval_dos_one;
  *dos_two_1 =        tre_thr;
  *dos_thr_1 = newval_dos_thr;
  *dos_fou_1 =        tre_fou;
  *tre_one_1 = newval_tre_one;
  *tre_two_1 = newval_tre_two;
  *tre_thr_1 = newval_tre_thr;
  *tre_fou_1 = newval_tre_fou;
  *qua_one_1 = newval_qua_one;
  *qua_two_1 =        qua_thr;
  *qua_thr_1 = newval_qua_thr;
  *qua_fou_1 =        qua_fou;
}

/*
 * LBB (Locally Bounded Bicubic) is a high quality nonlinear variant
 * of Catmull-Rom. Images resampled with LBB have much smaller halos
 * than images resampled with windowed sincs or other interpolatory
 * cubic spline filters. Specifically, LBB halos are narrower and the
 * over/undershoot amplitude is smaller. This is accomplished without
 * a significant reduction in the smoothness of the result (compared
 * to Catmull-Rom).
 *
 * Another important property is that the resampled values are
 * contained within the range of nearby input values. Consequently, no
 * final clamping is needed to stay "in range" (e.g., 0-255 for
 * standard 8-bit images).
 *
 * LBB was developed by N. Robidoux and C. Racette of the Department
 * of Mathematics and Computer Science of Laurentian University in the
 * course of C.'s Masters Thesis in Computational Sciences.
 */

/*
 * LBB is a novel method with the following properties:
 *
 * --LBB is a Hermite bicubic method: The bicubic surface is defined,
 *   one convex hull of four nearby input points at a time, using four
 *   point values, four x-derivatives, four y-derivatives, and four
 *   cross-derivatives.
 *
 * --The stencil for values in a square patch is the usual 4x4.
 *
 * --LBB is interpolatory.
 *
 * --It is C^1 with continuous cross derivatives.
 *
 * --When the limiters are inactive, LBB gives the same results as
 *   Catmull-Rom.
 *
 * --When used on binary images, LBB gives results similar to bicubic
 *   Hermite with all first derivatives---but not necessarily the
 *   cross derivatives--at the input pixel locations set to zero.
 *
 * --The LBB reconstruction is locally bounded: Over each square
 *   patch, the surface is contained between the minimum and the
 *   maximum values among the 16 nearest input pixel values (those in
 *   the stencil).
 *
 * --Consequently, the LBB reconstruction is globally bounded between
 *   the very smallest input pixel value and the very largest input
 *   pixel value. (It is not necessary to clamp results.)
 *
 * The LBB method is based on the method of Ken Brodlie, Petros
 * Mashwama and Sohail Butt for constraining Hermite interpolants
 * between globally defined planes:
 *
 *   Visualization of surface data to preserve positivity and other
 *   simple constraints. Computer & Graphics, Vol. 19, Number 4, pages
 *   585-594, 1995. DOI: 10.1016/0097-8493(95)00036-C.
 *
 * Instead of forcing the reconstructed surface to lie between two
 * GLOBALLY defined planes, LBB constrains one patch at a time to lie
 * between LOCALLY defined planes. This is accomplished by
 * constraining the derivatives (x, y and cross) at each input pixel
 * location so that if the constraint was applied everywhere the
 * surface would fit between the min and max of the values at the 9
 * closest pixel locations. Because this is done with each of the four
 * pixel locations which define the bicubic patch, this forces the
 * reconstructed surface to lie between the min and max of the values
 * at the 16 closest values pixel locations. (Each corner defines its
 * own 3x3 subgroup of the 4x4 stencil. Consequently, the surface is
 * necessarily above the minimum of the four minima, which happens to
 * be the minimum over the 4x4. Similarly with the maxima.)
 *
 * The above paragraph described the "soft" version of LBB. The
 * "sharp" version is similar.
 */

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 LBB is the same as for any standard Hermite
   * bicubic (e.g., Catmull-Rom):
   *
   *  (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.
   */

#if defined (__NOHALO_CHEAP_H__)
  /*
   * Computation of the four min and four max over 3x3 input data
   * sub-crosses of the 4x4 input stencil.
   *
   * We exploit the fact that the data comes from the (co-monotone)
   * method Nohalo so that it is known ahead of time that
   *
   *  dos_thr is between dos_two and dos_fou
   *
   *  tre_two is between dos_two and qua_two
   *
   *  tre_fou is between dos_fou and qua_fou
   *
   *  qua_thr is between qua_two and qua_fou
   *
   *  tre_thr is in the convex hull of dos_two, dos_fou, qua_two and qua_fou
   *
   *  to minimize the number of flags and conditional moves.
   *
   * (The "between" are not strict: "a between b and c" means
   *
   * "min(b,c) <= a <= max(b,c)".)
   *
   * We have, however, succeeded in eliminating one flag computation
   * (one comparison) and one use of an intermediate result. See the
   * two commented out lines below.
   *
   * Overall, only 20 comparisons and 28 "? :" are needed (to compute
   * 4 mins and 4 maxes). If you can figure how to do this more
   * efficiently, let us know.
   */
  const double m1    = (uno_two <= tre_two) ? uno_two : tre_two  ;
  const double M1    = (uno_two <= tre_two) ? tre_two : uno_two  ;
  const double m2    = (dos_thr <= qua_thr) ? dos_thr : qua_thr  ;
  const double M2    = (dos_thr <= qua_thr) ? qua_thr : dos_thr  ;
  const double m3    = (dos_two <= dos_fou) ? dos_two : dos_fou  ;
  const double M3    = (dos_two <= dos_fou) ? dos_fou : dos_two  ;
  const double m4    = (uno_thr <= tre_thr) ? uno_thr : tre_thr  ;
  const double M4    = (uno_thr <= tre_thr) ? tre_thr : uno_thr  ;
  const double m5    = (dos_two <= qua_two) ? dos_two : qua_two  ;
  const double M5    = (dos_two <= qua_two) ? qua_two : dos_two  ;
  const double m6    = (tre_one <= tre_thr) ? tre_one : tre_thr  ;
  const double M6    = (tre_one <= tre_thr) ? tre_thr : tre_one  ;
  const double m7    = (dos_one <= dos_thr) ? dos_one : dos_thr  ;
  const double M7    = (dos_one <= dos_thr) ? dos_thr : dos_one  ;
  const double m8    = (tre_two <= tre_fou) ? tre_two : tre_fou  ;
  const double M8    = (tre_two <= tre_fou) ? tre_fou : tre_two  ;
  const double m9    = NOHALO_MIN(            m1,       dos_two );
  const double M9    = NOHALO_MAX(            M1,       dos_two );
  const double m10   = NOHALO_MIN(            m2,       tre_thr );
  const double M10   = NOHALO_MAX(            M2,       tre_thr );
  const double min10 = NOHALO_MIN(            m3,       m4      );
  const double max10 = NOHALO_MAX(            M3,       M4      );
  const double min01 = NOHALO_MIN(            m5,       m6      );
  const double max01 = NOHALO_MAX(            M5,       M6      );
  const double min00 = NOHALO_MIN(            m9,       m7      );
  const double max00 = NOHALO_MAX(            M9,       M7      );
  const double min11 = NOHALO_MIN(           m10,       m8      );
  const double max11 = NOHALO_MAX(           M10,       M8      );
#else
  /*
   * Computation of the four min and four max over 3x3 input data
   * sub-blocks of the 4x4 input stencil.
   *
   * Surprisingly, we have not succeeded in reducing the number of "?
   * :" needed by using the fact that the data comes from the
   * (co-monotone) method Nohalo so that it is known ahead of time
   * that
   *
   *  dos_thr is between dos_two and dos_fou
   *
   *  tre_two is between dos_two and qua_two
   *
   *  tre_fou is between dos_fou and qua_fou
   *
   *  qua_thr is between qua_two and qua_fou
   *
   *  tre_thr is in the convex hull of dos_two, dos_fou, qua_two and qua_fou
   *
   *  to minimize the number of flags and conditional moves.
   *
   * (The "between" are not strict: "a between b and c" means
   *
   * "min(b,c) <= a <= max(b,c)".)
   *
   * We have, however, succeeded in eliminating one flag computation
   * (one comparison) and one use of an intermediate result. See the
   * two commented out lines below.
   *
   * Overall, only 27 comparisons are needed (to compute 4 mins and 4
   * maxes!). Without the simplification, 28 comparisons would be
   * used. Either way, the number of "? :" used is 34. If you can
   * figure how to do this more efficiently, let us know.
   */
  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 m4    = (qua_two <= qua_thr) ? qua_two : qua_thr  ;
  const double M4    = (qua_two <= qua_thr) ? qua_thr : qua_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 m5    = NOHALO_MIN(            m1,       m2      );
  const double M5    = NOHALO_MAX(            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 m13   = (dos_fou <= qua_fou) ? dos_fou : qua_fou  ;
  const double M13   = (dos_fou <= qua_fou) ? qua_fou : dos_fou  ;
  /*
   * Because the data comes from Nohalo subdivision, the following two
   * lines can be replaced by the above, simpler, two lines without
   * changing the results.
   *
   * const double m13   = NOHALO_MIN(               m7,       qua_fou );
   * const double M13   = NOHALO_MAX(               M7,       qua_fou );
   *
   * This also allows reodering the comparisons to put space between
   * the computation of a result and its use.
   */
  const double m9    = NOHALO_MIN(            m5,       m4      );
  const double M9    = NOHALO_MAX(            M5,       M4      );
  const double m11   = NOHALO_MIN(            m6,       qua_one );
  const double M11   = NOHALO_MAX(            M6,       qua_one );
  const double m10   = NOHALO_MIN(            m6,       uno_one );
  const double M10   = NOHALO_MAX(            M6,       uno_one );
  const double m8    = NOHALO_MIN(            m5,       m3      );
  const double M8    = NOHALO_MAX(            M5,       M3      );
  const double m12   = NOHALO_MIN(            m7,       uno_fou );
  const double M12   = NOHALO_MAX(            M7,       uno_fou );
  const double min11 = NOHALO_MIN(            m9,       m13     );
  const double max11 = NOHALO_MAX(            M9,       M13     );
  const double min01 = NOHALO_MIN(            m9,       m11     );
  const double max01 = NOHALO_MAX(            M9,       M11     );
  const double min00 = NOHALO_MIN(            m8,       m10     );
  const double max00 = NOHALO_MAX(            M8,       M10     );
  const double min10 = NOHALO_MIN(            m8,       m12     );
  const double max10 = NOHALO_MAX(            M8,       M12     );
#endif

  /*
   * The remainder of the "per channel" computation involves the
   * computation of:
   *
   * --8 conditional moves,
   *
   * --8 signs (in which the sign of zero is unimportant),
   *
   * --12 minima of two values,
   *
   * --8 maxima of two values,
   *
   * --8 absolute values,
   *
   * for a grand total of 29 minima, 25 maxima, 8 conditional moves, 8
   * signs, and 8 absolute values. If everything is done with
   * conditional moves, "only" 28+8+8+12+8+8=72 flags are involved
   * (because initial min and max can be computed with one flag).
   *
   * The "per channel" part of the computation also involves 107
   * arithmetic operations (54 *, 21 +, 42 -).
   */

  /*
   * Distances to the local min and max:
   */
  const double u11 = tre_thr - min11;
  const double v11 = max11 - tre_thr;
  const double u01 = tre_two - min01;
  const double v01 = max01 - tre_two;
  const double u00 = dos_two - min00;
  const double v00 = max00 - dos_two;
  const double u10 = dos_thr - min10;
  const double v10 = max10 - dos_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_dzdy11i = qua_thr - dos_thr;
  const double dble_dzdx10i = dos_fou - dos_two;
  const double dble_dzdy01i = qua_two - dos_two;
  const double dble_dzdx01i = tre_thr - tre_one;
  const double dble_dzdy10i = tre_thr - uno_thr;
  const double dble_dzdx11i = tre_fou - tre_two;
  const double dble_dzdy00i = tre_two - uno_two;

  /*
   * Signs of the derivatives. The upcoming clamping does not change
   * them (except if the clamping sends a negative derivative to 0, in
   * which case the sign does not matter anyway).
   */
  const double sign_dzdx00 = NOHALO_SIGN( dble_dzdx00i );
  const double sign_dzdx10 = NOHALO_SIGN( dble_dzdx10i );
  const double sign_dzdx01 = NOHALO_SIGN( dble_dzdx01i );
  const double sign_dzdx11 = NOHALO_SIGN( dble_dzdx11i );

  const double sign_dzdy00 = NOHALO_SIGN( dble_dzdy00i );
  const double sign_dzdy10 = NOHALO_SIGN( dble_dzdy10i );
  const double sign_dzdy01 = NOHALO_SIGN( dble_dzdy01i );
  const double sign_dzdy11 = NOHALO_SIGN( dble_dzdy11i );

  /*
   * 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;

  /*
   * Slope limiters. The key multiplier is 3 but we fold a factor of
   * 2, hence 6:
   */
  const double dble_slopelimit_00 = 6.0 * NOHALO_MIN( u00, v00 );
  const double dble_slopelimit_10 = 6.0 * NOHALO_MIN( u10, v10 );
  const double dble_slopelimit_01 = 6.0 * NOHALO_MIN( u01, v01 );
  const double dble_slopelimit_11 = 6.0 * NOHALO_MIN( u11, v11 );

  /*
   * Clamped first derivatives:
   */
  const double dble_dzdx00 =
    ( sign_dzdx00 * dble_dzdx00i <= dble_slopelimit_00 )
    ? dble_dzdx00i :  sign_dzdx00 * dble_slopelimit_00;
  const double dble_dzdy00 =
    ( sign_dzdy00 * dble_dzdy00i <= dble_slopelimit_00 )
    ? dble_dzdy00i :  sign_dzdy00 * dble_slopelimit_00;
  const double dble_dzdx10 =
    ( sign_dzdx10 * dble_dzdx10i <= dble_slopelimit_10 )
    ? dble_dzdx10i :  sign_dzdx10 * dble_slopelimit_10;
  const double dble_dzdy10 =
    ( sign_dzdy10 * dble_dzdy10i <= dble_slopelimit_10 )
    ? dble_dzdy10i :  sign_dzdy10 * dble_slopelimit_10;
  const double dble_dzdx01 =
    ( sign_dzdx01 * dble_dzdx01i <= dble_slopelimit_01 )
    ? dble_dzdx01i :  sign_dzdx01 * dble_slopelimit_01;
  const double dble_dzdy01 =
    ( sign_dzdy01 * dble_dzdy01i <= dble_slopelimit_01 )
    ? dble_dzdy01i :  sign_dzdy01 * dble_slopelimit_01;
  const double dble_dzdx11 =
    ( sign_dzdx11 * dble_dzdx11i <= dble_slopelimit_11 )
    ? dble_dzdx11i :  sign_dzdx11 * dble_slopelimit_11;
  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 );

  /*
   * Absolute values of the sums:
   */
  const double twelve_abs_sum00 = NOHALO_ABS( twelve_sum00 );
  const double twelve_abs_sum10 = NOHALO_ABS( twelve_sum10 );
  const double twelve_abs_sum01 = NOHALO_ABS( twelve_sum01 );
  const double twelve_abs_sum11 = NOHALO_ABS( twelve_sum11 );

  /*
   * Scaled distances to the min:
   */
  const double u00_times_36 = 36.0 * u00;
  const double u10_times_36 = 36.0 * u10;
  const double u01_times_36 = 36.0 * u01;
  const double u11_times_36 = 36.0 * 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 = NOHALO_MAX( quad_d2zdxdy00i, first_limit00 );
  const double quad_d2zdxdy10ii = NOHALO_MAX( quad_d2zdxdy10i, first_limit10 );
  const double quad_d2zdxdy01ii = NOHALO_MAX( quad_d2zdxdy01i, first_limit01 );
  const double quad_d2zdxdy11ii = NOHALO_MAX( quad_d2zdxdy11i, first_limit11 );

  /*
   * Scaled distances to the max:
   */
  const double v00_times_36 = 36.0 * v00;
  const double v10_times_36 = 36.0 * v10;
  const double v01_times_36 = 36.0 * v01;
  const double v11_times_36 = 36.0 * 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 =
    NOHALO_MIN( quad_d2zdxdy00ii, second_limit00 );
  const double quad_d2zdxdy10iii =
    NOHALO_MIN( quad_d2zdxdy10ii, second_limit10 );
  const double quad_d2zdxdy01iii =
    NOHALO_MIN( quad_d2zdxdy01ii, second_limit01 );
  const double quad_d2zdxdy11iii =
    NOHALO_MIN( quad_d2zdxdy11ii, second_limit11 );

  /*
   * Absolute values of the differences:
   */
  const double twelve_abs_dif00 = NOHALO_ABS( twelve_dif00 );
  const double twelve_abs_dif10 = NOHALO_ABS( twelve_dif10 );
  const double twelve_abs_dif01 = NOHALO_ABS( twelve_dif01 );
  const double twelve_abs_dif11 = NOHALO_ABS( twelve_dif11 );

  /*
   * Third cross-derivative limiter:
   */
  const double third_limit00 = twelve_abs_dif00 - v00_times_36;
  const double third_limit10 = twelve_abs_dif10 - v10_times_36;
  const double third_limit01 = twelve_abs_dif01 - v01_times_36;
  const double third_limit11 = twelve_abs_dif11 - v11_times_36;

  const double quad_d2zdxdy00iiii =
    NOHALO_MAX( quad_d2zdxdy00iii, third_limit00);
  const double quad_d2zdxdy10iiii =
    NOHALO_MAX( quad_d2zdxdy10iii, third_limit10);
  const double quad_d2zdxdy01iiii =
    NOHALO_MAX( quad_d2zdxdy01iii, third_limit01);
  const double quad_d2zdxdy11iiii =
    NOHALO_MAX( quad_d2zdxdy11iii, third_limit11);

  /*
   * Fourth cross-derivative limiter:
   */
  const double fourth_limit00 = u00_times_36 - twelve_abs_dif00;
  const double fourth_limit10 = u10_times_36 - twelve_abs_dif10;
  const double fourth_limit01 = u01_times_36 - twelve_abs_dif01;
  const double fourth_limit11 = u11_times_36 - twelve_abs_dif11;

  const double quad_d2zdxdy00 = NOHALO_MIN( quad_d2zdxdy00iiii, fourth_limit00);
  const double quad_d2zdxdy10 = NOHALO_MIN( quad_d2zdxdy10iiii, fourth_limit10);
  const double quad_d2zdxdy01 = NOHALO_MIN( quad_d2zdxdy01iiii, fourth_limit01);
  const double quad_d2zdxdy11 = NOHALO_MIN( quad_d2zdxdy11iiii, fourth_limit11);

  /*
   * Part of the result which does not need derivatives:
   */
  const double newval1 = c00 * dos_two
                         +
                         c10 * dos_thr
                         +
                         c01 * tre_two
                         +
                         c11 * tre_thr;

  /*
   * Twice the 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;

  /*
   * 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 Nohalo+LBB with a careful type conversion 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 NOHALO_CONVERSION( conversion )               \
  template <typename T> static void inline            \
  nohalo_ ## conversion(       void*  restrict pout,  \
                         const void*  restrict pin,   \
                         const int             bands, \
                         const int             lskip, \
                         const double          x_0,   \
                         const double          y_0 )  \
  { \
    T* restrict out = (T *) pout; \
    \
    const T* restrict in = (T *) pin; \
    \
    \
    const int sign_of_x_0 = 2 * ( x_0 >= 0. ) - 1; \
    const int sign_of_y_0 = 2 * ( y_0 >= 0. ) - 1; \
    \
    \
    const int shift_forw_1_pix = sign_of_x_0 * bands; \
    const int shift_forw_1_row = sign_of_y_0 * lskip; \
    \
    const int shift_back_1_pix = -shift_forw_1_pix; \
    const int shift_back_1_row = -shift_forw_1_row; \
    \
    const int shift_back_2_pix = 2 * shift_back_1_pix; \
    const int shift_back_2_row = 2 * shift_back_1_row; \
    const int shift_forw_2_pix = 2 * shift_forw_1_pix; \
    const int shift_forw_2_row = 2 * shift_forw_1_row; \
    \
    \
    const int uno_two_shift = shift_back_1_pix + shift_back_2_row; \
    const int uno_thr_shift =                    shift_back_2_row; \
    const int uno_fou_shift = shift_forw_1_pix + shift_back_2_row; \
    \
    const int dos_one_shift = shift_back_2_pix + shift_back_1_row; \
    const int dos_two_shift = shift_back_1_pix + shift_back_1_row; \
    const int dos_thr_shift =                    shift_back_1_row; \
    const int dos_fou_shift = shift_forw_1_pix + shift_back_1_row; \
    const int dos_fiv_shift = shift_forw_2_pix + shift_back_1_row; \
    \
    const int tre_one_shift = shift_back_2_pix; \
    const int tre_two_shift = shift_back_1_pix; \
    const int tre_thr_shift = 0;                \
    const int tre_fou_shift = shift_forw_1_pix; \
    const int tre_fiv_shift = shift_forw_2_pix; \
    \
    const int qua_one_shift = shift_back_2_pix + shift_forw_1_row; \
    const int qua_two_shift = shift_back_1_pix + shift_forw_1_row; \
    const int qua_thr_shift =                    shift_forw_1_row; \
    const int qua_fou_shift = shift_forw_1_pix + shift_forw_1_row; \
    const int qua_fiv_shift = shift_forw_2_pix + shift_forw_1_row; \
    \
    const int cin_two_shift = shift_back_1_pix + shift_forw_2_row; \
    const int cin_thr_shift =                    shift_forw_2_row; \
    const int cin_fou_shift = shift_forw_1_pix + shift_forw_2_row; \
    \
    \
    const double xp1over2   = ( 2 * sign_of_x_0 ) * x_0; \
    const double xm1over2   = xp1over2 - 1.0; \
    const double onepx      = 0.5 + xp1over2; \
    const double onemx      = 1.5 - xp1over2; \
    const double xp1over2sq = xp1over2 * xp1over2; \
    \
    const double yp1over2   = ( 2 * sign_of_y_0 ) * y_0; \
    const double ym1over2   = yp1over2 - 1.0; \
    const double onepy      = 0.5 + yp1over2; \
    const double onemy      = 1.5 - yp1over2; \
    const double yp1over2sq = yp1over2 * yp1over2; \
    \
    const double xm1over2sq = xm1over2 * xm1over2; \
    const double ym1over2sq = ym1over2 * ym1over2; \
    \
    const double twice1px = onepx + onepx; \
    const double twice1py = onepy + onepy; \
    const double twice1mx = onemx + onemx; \
    const double twice1my = onemy + onemy; \
    \
    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 four_times_1px_times_1py = twice1px * twice1py; \
    const double four_times_1mx_times_1py = twice1mx * twice1py; \
    const double twice_xp1over2_times_1py = xp1over2 * twice1py; \
    const double twice_xm1over2_times_1py = xm1over2 * twice1py; \
    \
    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_1px_times_ym1over2 = twice1px * ym1over2; \
    const double twice_1mx_times_ym1over2 = twice1mx * ym1over2; \
    const double xp1over2_times_ym1over2  = xp1over2 * ym1over2; \
    const double xm1over2_times_ym1over2  = xm1over2 * 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 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 \
      { \
        double uno_one, uno_two, uno_thr, uno_fou;  \
        double dos_one, dos_two, dos_thr, dos_fou;  \
        double tre_one, tre_two, tre_thr, tre_fou;  \
        double qua_one, qua_two, qua_thr, qua_fou;  \
        \
        nohalo_subdivision( 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[ dos_fiv_shift ], \
                            in[ tre_one_shift ], \
                            in[ tre_two_shift ], \
                            in[ tre_thr_shift ], \
                            in[ tre_fou_shift ], \
                            in[ tre_fiv_shift ], \
                            in[ qua_one_shift ], \
                            in[ qua_two_shift ], \
                            in[ qua_thr_shift ], \
                            in[ qua_fou_shift ], \
                            in[ qua_fiv_shift ], \
                            in[ cin_two_shift ], \
                            in[ cin_thr_shift ], \
                            in[ cin_fou_shift ], \
                            &uno_one,            \
                            &uno_two,            \
                            &uno_thr,            \
                            &uno_fou,            \
                            &dos_one,            \
                            &dos_two,            \
                            &dos_thr,            \
                            &dos_fou,            \
                            &tre_one,            \
                            &tre_two,            \
                            &tre_thr,            \
                            &tre_fou,            \
                            &qua_one,            \
                            &qua_two,            \
                            &qua_thr,            \
                            &qua_fou );          \
        \
        const double double_result =        \
          lbbicubic( c00,                   \
                     c10,                   \
                     c01,                   \
                     c11,                   \
                     c00dx,                 \
                     c10dx,                 \
                     c01dx,                 \
                     c11dx,                 \
                     c00dy,                 \
                     c10dy,                 \
                     c01dy,                 \
                     c11dy,                 \
                     c00dxdy,               \
                     c10dxdy,               \
                     c01dxdy,               \
                     c11dxdy,               \
                     uno_one,               \
                     uno_two,               \
                     uno_thr,               \
                     uno_fou,               \
                     dos_one,               \
                     dos_two,               \
                     dos_thr,               \
                     dos_fou,               \
                     tre_one,               \
                     tre_two,               \
                     tre_thr,               \
                     tre_fou,               \
                     qua_one,               \
                     qua_two,               \
                     qua_thr,               \
                     qua_fou );             \
        \
        {                                                         \
          const T result = to_ ## conversion<T>( double_result ); \
          in++;                                                   \
          *out++ = result;                                        \
        }                                                         \
        \
      } while (--band); \
  }


NOHALO_CONVERSION( fptypes )
NOHALO_CONVERSION( withsign )
NOHALO_CONVERSION( nosign )


#define CALL( T, conversion )             \
  nohalo_ ## conversion<T>( out,          \
                            p,            \
                            bands,        \
                            lskip,        \
                            relative_x,   \
                            relative_y );


/*
 * We need C linkage:
 */
extern "C" {
G_DEFINE_TYPE( VipsInterpolateNohalo, vips_interpolate_nohalo,
	VIPS_TYPE_INTERPOLATE );
}


static void
vips_interpolate_nohalo_interpolate( VipsInterpolate* restrict interpolate,
                                     void*            restrict out,
                                     VipsRegion*      restrict in,
                                     double                    absolute_x,
                                     double                    absolute_y )
{
  /* absolute_x and absolute_y are always >= 2.0 (see double-check assert
   * below), so we don't need floor(). 
   *
   * It's 2 not 0 since we ask for a window_offset of 2 at the bottom.
   */
  const int ix = (int) (absolute_x + 0.5);
  const int iy = (int) (absolute_y + 0.5);

  /*
   * 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 VipsPel* restrict p = VIPS_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 lskip = VIPS_REGION_LSKIP( in ) / 
	  VIPS_IMAGE_SIZEOF_ELEMENT( in->im );

  /*
   * Double the bands for complex images to account for the real and
   * imaginary parts being computed independently:
   */
  const int actual_bands = in->im->Bands;
  const int bands =
    vips_band_format_iscomplex( in->im->BandFmt ) ? 
      2 * actual_bands : actual_bands;

  g_assert( ix - 2 >= in->valid.left );
  g_assert( iy - 2 >= in->valid.top );
  g_assert( ix + 2 <= VIPS_RECT_RIGHT( &in->valid ) );
  g_assert( iy + 2 <= VIPS_RECT_BOTTOM( &in->valid ) );

  /* Confirm that absolute_x and absolute_y are >= 2, see above. 
   */
  g_assert( absolute_x >= 2.0 );
  g_assert( absolute_y >= 2.0 );

  switch( in->im->BandFmt ) {
  case VIPS_FORMAT_UCHAR:
    CALL( unsigned char, nosign );
    break;

  case VIPS_FORMAT_CHAR:
    CALL( signed char, withsign );
    break;

  case VIPS_FORMAT_USHORT:
    CALL( unsigned short, nosign );
    break;

  case VIPS_FORMAT_SHORT:
    CALL( signed short, withsign );
    break;

  case VIPS_FORMAT_UINT:
    CALL( unsigned int, nosign );
    break;

  case VIPS_FORMAT_INT:
    CALL( signed int, withsign );
    break;

  /*
   * Complex images are handled by doubling of bands.
   */
  case VIPS_FORMAT_FLOAT:
  case VIPS_FORMAT_COMPLEX:
    CALL( float, fptypes );
    break;

  case VIPS_FORMAT_DOUBLE:
  case VIPS_FORMAT_DPCOMPLEX:
    CALL( double, fptypes );
    break;

  default:
    g_assert( 0 );
    break;
  }
}

static void
vips_interpolate_nohalo_class_init( VipsInterpolateNohaloClass *klass )
{
  VipsObjectClass *object_class = VIPS_OBJECT_CLASS( klass );
  VipsInterpolateClass *interpolate_class = VIPS_INTERPOLATE_CLASS( klass );

  object_class->nickname    = "nohalo";
  object_class->description =
    _( "edge sharpening resampler with halo reduction" );

  interpolate_class->interpolate   = vips_interpolate_nohalo_interpolate;
  interpolate_class->window_size   = 6;
  interpolate_class->window_offset = 2;
}

static void
vips_interpolate_nohalo_init( VipsInterpolateNohalo *nohalo )
{
}