416 lines
14 KiB
C++
416 lines
14 KiB
C++
/* 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 <stdio.h>
|
|
#include <stdlib.h>
|
|
|
|
#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); \
|
|
}
|
|
|
|
|
|
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 )
|
|
{
|
|
}
|