/* vertex-split subdivision followed by quadratic b-spline smoothing

 *

 * C. Racette 23-28/05/2010 based on code by N. Robidoux and J. Cupitt

 *

 * N. Robidoux 29-30/05/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

 */

/*

 * 2010 (c) Chantal Racette, Nicolas Robidoux, John Cupitt.

 *

 * Nicolas Robidoux thanks Adam Turcotte, Geert Jordaens, Ralf Meyer,

 * Øyvind Kolås, Minglun Gong and Eric Daoust for useful comments and

 * code.

 *

 * Chantal Racette's image resampling research and programming funded

 * in part by a NSERC Discovery Grant awarded to Julien Dompierre

 * (20-61098).

 */

/*

 * Vertex-Split Quadratic B-Splines (VSQBS) is a brand new method

 * which consists of vertex-split subdivision, a subdivision method

 * with the (as yet unknown?) property that data which is (locally)

 * constant on diagonals is subdivided into data which is (locally)

 * constant on diagonals, followed by quadratic B-Spline smoothing.

 * Because both methods are linear, their combination can be

 * implemented as if there is no subdivision.

 *

 * At high enlargement ratios, VSQBS is very effective at "masking"

 * that the original has pixels uniformly distributed on a grid. In

 * particular, VSQBS produces resamples with only very mild

 * staircasing. Like cubic B-Spline smoothing, however, VSQBS is not

 * an interpolatory method. For example, using VSQBS to perform the

 * identity geometric transformation (enlargement by a scaling factor

 * equal to 1) on an image does not return the original: VSQBS

 * effectively smooths out the image with the convolution mask

 *

 *     1/8

 * 1/8 1/2 1/8

 *     1/8

 *

 * which is a fairly moderate blur (although the checkerboard mode is

 * in its nullspace).

 *

 * By blending VSQBS with an interpolatory method (bilinear, say) in a

 * transformation adaptive environment (current GEGL, for example), it

 * is quite easy to restore that resampling for identity geometric

 * transformation is equivalent to the identity, and rotations are not

 * affected by the above, implicit, blur. Contact N. Robidoux for

 * details.

 *

 * An article on VSQBS is forthcoming.

 */

#ifdef HAVE_CONFIG_H

#include <config.h>

#endif /*HAVE_CONFIG_H*/

#include <glib/gi18n-lib.h>

#include <cstdio>

#include <cstdlib>

#include <vips/vips.h>

#include <vips/internal.h>

#include "templates.h"

#define VIPS_TYPE_INTERPOLATE_VSQBS \

  (vips_interpolate_vsqbs_get_type())

#define VIPS_INTERPOLATE_VSQBS(obj) \

  (G_TYPE_CHECK_INSTANCE_CAST((obj), \

    VIPS_TYPE_INTERPOLATE_VSQBS, VipsInterpolateVsqbs))

#define VIPS_INTERPOLATE_VSQBS_CLASS(klass) \

  (G_TYPE_CHECK_CLASS_CAST((klass), \

    VIPS_TYPE_INTERPOLATE_VSQBS, VipsInterpolateVsqbsClass))

#define VIPS_IS_INTERPOLATE_VSQBS(obj) \

  (G_TYPE_CHECK_INSTANCE_TYPE((obj), VIPS_TYPE_INTERPOLATE_VSQBS))

#define VIPS_IS_INTERPOLATE_VSQBS_CLASS(klass) \

  (G_TYPE_CHECK_CLASS_TYPE((klass), VIPS_TYPE_INTERPOLATE_VSQBS))

#define VIPS_INTERPOLATE_VSQBS_GET_CLASS(obj) \

  (G_TYPE_INSTANCE_GET_CLASS((obj), \

    VIPS_TYPE_INTERPOLATE_VSQBS, VipsInterpolateVsqbsClass))

typedef struct _VipsInterpolateVsqbs {

  VipsInterpolate parent_object;

} VipsInterpolateVsqbs;

typedef struct _VipsInterpolateVsqbsClass {

  VipsInterpolateClass parent_class;

} VipsInterpolateVsqbsClass;

/*

 * THE STENCIL OF INPUT VALUES:

 *

 * Pointer arithmetic is used to implicitly reflect the input stencil

 * about dos_two---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 dos_two.

 *

 * The following code and picture assumes that the stencil reflexion

 * has already been performed. (X is the sampling location.)

 *

 *

 *               (ix,iy-1)    (ix+1,iy-1)

 *               = uno_two    = uno_thr

 *

 *

 *

 *  (ix-1,iy)    (ix,iy)      (ix+1,iy)

 *  = dos_one    = dos_two    = dos_thr

 *                       X

 *

 *

 *  (ix-1,iy+1)  (ix,iy+1)    (ix+1,iy+1)

 *  = tre_one    = tre_two    = tre_thr

 *

 *

 * The above input pixel values are the ones needed in order to

 * IMPLICITLY make available the following values, needed by quadratic

 * B-Splines, which is performed on (shifted) double density data:

 *

 *

 *  uno_one_1 =      uno_two_1 =      uno_thr_1 =

 *  (ix-1/4,iy-1/4)  (ix+1/4,iy-1/4)  (ix+3/4,iy-1/4)

 *

 *

 *

 *                 X            or X

 *  dos_one_1 =      dos_two_1 =      dos_thr_1 =

 *  (ix-1/4,iy+1/4)  (ix+1/4,iy+1/4)  (ix+3/4,iy+1/4)

 *              or X            or X

 *

 *

 *

 *  tre_one_1 =      tre_two_1 =      tre_thr_1 =

 *  (ix-1/4,iy+3/4)  (ix+1/4,iy+3/4)  (ix+3/4,iy+3/4)

 *

 *

 * In the coefficient computations, we fix things so that coordinates

 * are relative to dos_two_1, and so that distances are rescaled so

 * that double density pixel locations are at a distance of 1.

 */

/*

 * Call vertex-split + quadratic B-splines 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 VSQBS_CONVERSION(conversion) \

  template <typename T> \

  static void inline vsqbs_##conversion(void *restrict pout, \

    const VipsPel *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 uno_two_shift = shift_back_1_row; \

    const int uno_thr_shift = shift_forw_1_pix + shift_back_1_row; \

\

    const int dos_one_shift = shift_back_1_pix; \

    const int dos_two_shift = 0; \

    const int dos_thr_shift = shift_forw_1_pix; \

\

    const int tre_one_shift = shift_back_1_pix + shift_forw_1_row; \

    const int tre_two_shift = shift_forw_1_row; \

    const int tre_thr_shift = shift_forw_1_pix + shift_forw_1_row; \

\

    const double twice_abs_x_0 = (2 * sign_of_x_0) * x_0; \

    const double twice_abs_y_0 = (2 * sign_of_y_0) * y_0; \

    const double x = twice_abs_x_0 + -0.5; \

    const double y = twice_abs_y_0 + -0.5; \

    const double cent = 0.75 - x * x; \

    const double mid = 0.75 - y * y; \

    const double left = -0.5 * (x + cent) + 0.5; \

    const double top = -0.5 * (y + mid) + 0.5; \

    const double left_p_cent = left + cent; \

    const double top_p_mid = top + mid; \

    const double cent_p_rite = 1.0 - left; \

    const double mid_p_bot = 1.0 - top; \

    const double rite = 1.0 - left_p_cent; \

    const double bot = 1.0 - top_p_mid; \

\

    const double four_c_uno_two = left_p_cent * top; \

    const double four_c_dos_one = left * top_p_mid; \

    const double four_c_dos_two = left_p_cent + top_p_mid; \

    const double four_c_dos_thr = cent_p_rite * top_p_mid + rite; \

    const double four_c_tre_two = mid_p_bot * left_p_cent + bot; \

    const double four_c_tre_thr = mid_p_bot * rite + cent_p_rite * bot; \

    const double four_c_uno_thr = top - four_c_uno_two; \

    const double four_c_tre_one = left - four_c_dos_one; \

\

    int band = bands; \

\

    do { \

      const double double_result = \

        (((four_c_uno_two * in[uno_two_shift] + \

            four_c_dos_one * in[dos_one_shift]) + \

           (four_c_dos_two * in[dos_two_shift] + \

             four_c_dos_thr * in[dos_thr_shift])) + \

          ((four_c_tre_two * in[tre_two_shift] + \

             four_c_tre_thr * in[tre_thr_shift]) + \

            (four_c_uno_thr * in[uno_thr_shift] + \

              four_c_tre_one * in[tre_one_shift]))) * \

        0.25; \

\

      const T result = to_##conversion<T>(double_result); \

      in++; \

      *out++ = result; \

\

    } while (--band); \

  }

Unexecuted instantiation: vsqbs.cpp:void vsqbs_nosign<unsigned char>(void*, unsigned char const*, int, int, double, double)

Unexecuted instantiation: vsqbs.cpp:void vsqbs_withsign<signed char>(void*, unsigned char const*, int, int, double, double)

Unexecuted instantiation: vsqbs.cpp:void vsqbs_nosign<unsigned short>(void*, unsigned char const*, int, int, double, double)

Unexecuted instantiation: vsqbs.cpp:void vsqbs_withsign<short>(void*, unsigned char const*, int, int, double, double)

Unexecuted instantiation: vsqbs.cpp:void vsqbs_nosign<unsigned int>(void*, unsigned char const*, int, int, double, double)

Unexecuted instantiation: vsqbs.cpp:void vsqbs_withsign<int>(void*, unsigned char const*, int, int, double, double)

Unexecuted instantiation: vsqbs.cpp:void vsqbs_fptypes<float>(void*, unsigned char const*, int, int, double, double)

Unexecuted instantiation: vsqbs.cpp:void vsqbs_fptypes<double>(void*, unsigned char const*, int, int, double, double)

VSQBS_CONVERSION(fptypes)

VSQBS_CONVERSION(withsign)

VSQBS_CONVERSION(nosign)

#define CALL(T, conversion) \

  vsqbs_##conversion<T>(out, \

    p, \

    bands, \

    lskip, \

    relative_x, \

    relative_y);

/*

 * We need C linkage:

 */

extern "C" {

G_DEFINE_TYPE(VipsInterpolateVsqbs, vips_interpolate_vsqbs,

  VIPS_TYPE_INTERPOLATE);

}

static void

vips_interpolate_vsqbs_interpolate(VipsInterpolate *restrict interpolate,

  void *restrict out,

  VipsRegion *restrict in,

  double absolute_x,

  double absolute_y)

{

  /* absolute_x and absolute_y are always >= 1.0 (see double-check assert

   * below), so we don't need floor().

   *

   * It's 1 not 0 since we ask for a window_offset of 1 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 - 1 >= in->valid.left);

  g_assert(iy - 1 >= in->valid.top);

  g_assert(ix + 1 <= VIPS_RECT_RIGHT(&in->valid));

  g_assert(iy + 1 <= VIPS_RECT_BOTTOM(&in->valid));

  /* Confirm that absolute_x and absolute_y are >= 1, see above.

   */

  g_assert(absolute_x >= 1.0);

  g_assert(absolute_y >= 1.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 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_vsqbs_class_init(VipsInterpolateVsqbsClass *klass)

{

  VipsObjectClass *object_class = VIPS_OBJECT_CLASS(klass);

  VipsInterpolateClass *interpolate_class = VIPS_INTERPOLATE_CLASS(klass);

  object_class->nickname = "vsqbs";

  object_class->description = _("B-Splines with antialiasing smoothing");

  interpolate_class->interpolate = vips_interpolate_vsqbs_interpolate;

  interpolate_class->window_size = 4;

  interpolate_class->window_offset = 1;

}

static void

vips_interpolate_vsqbs_init(VipsInterpolateVsqbs *vsqbs)

{

}