| /* 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., 59 Temple Place, Suite 330, Boston, MA |
| 02111-1307 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 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 <vips/intl.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( PEL* restrict pout, \ |
| const PEL* 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 = top * left_p_cent; \ |
| 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 + bot * cent_p_rite; \ |
| 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, |
| 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 + .5 ); |
| const int iy = FAST_PSEUDO_FLOOR( absolute_y + .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 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 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_vsqbs_class_init( VipsInterpolateVsqbsClass *klass ) |
| { |
| GObjectClass *gobject_class = G_OBJECT_CLASS( klass ); |
| VipsObjectClass *object_class = VIPS_OBJECT_CLASS( klass ); |
| VipsInterpolateClass *interpolate_class = VIPS_INTERPOLATE_CLASS( klass ); |
| |
| gobject_class->set_property = vips_object_set_property; |
| gobject_class->get_property = vips_object_get_property; |
| |
| object_class->nickname = "vsqbs"; |
| object_class->description = _( "B-Splines with antialiasing smoothing" ); |
| |
| interpolate_class->interpolate = vips_interpolate_vsqbs_interpolate; |
| interpolate_class->window_size = 3; |
| interpolate_class->window_offset = 1; |
| } |
| |
| static void |
| vips_interpolate_vsqbs_init( VipsInterpolateVsqbs *vsqbs ) |
| { |
| } |
