diff --git a/include/vips/interpolate.h b/include/vips/interpolate.h index c927d872..1e5670c0 100644 --- a/include/vips/interpolate.h +++ b/include/vips/interpolate.h @@ -147,9 +147,6 @@ VipsInterpolate *vips_interpolate_nearest_static( void ); typedef struct _VipsInterpolateBilinear { VipsInterpolate parent_object; - /* Set this to not use tables ...slightly more accurate. - */ - gboolean slow; } VipsInterpolateBilinear; typedef struct _VipsInterpolateBilinearClass { @@ -164,13 +161,51 @@ typedef struct _VipsInterpolateBilinearClass { } VipsInterpolateBilinearClass; GType vips_interpolate_bilinear_get_type( void ); -void vips_interpolate_bilinear_set_slow( VipsInterpolateBilinear *, gboolean ); VipsInterpolate *vips_interpolate_bilinear_new( void ); /* Convenience: return a static fast bilinear, so no need to free it. */ VipsInterpolate *vips_interpolate_bilinear_static( void ); +/* Slow bilinear class starts. + */ + +#define VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW \ + (vips_interpolate_bilinear_slow_get_type()) +#define VIPS_INTERPOLATE_BILINEAR_SLOW( obj ) \ + (G_TYPE_CHECK_INSTANCE_CAST( (obj), \ + VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW, VipsInterpolateBilinearSlow )) +#define VIPS_INTERPOLATE_BILINEAR_SLOW_CLASS( klass ) \ + (G_TYPE_CHECK_CLASS_CAST( (klass), \ + VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW, \ + VipsInterpolateBilinearSlowClass)) +#define VIPS_IS_INTERPOLATE_BILINEAR_SLOW( obj ) \ + (G_TYPE_CHECK_INSTANCE_TYPE( (obj), \ + VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW )) +#define VIPS_IS_INTERPOLATE_BILINEAR_SLOW_CLASS( klass ) \ + (G_TYPE_CHECK_CLASS_TYPE( (klass), \ + VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW )) +#define VIPS_INTERPOLATE_BILINEAR_SLOW_GET_CLASS( obj ) \ + (G_TYPE_INSTANCE_GET_CLASS( (obj), \ + VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW, VipsInterpolateBilinearSlowClass )) + +typedef struct _VipsInterpolateBilinearSlow { + VipsInterpolate parent_object; + +} VipsInterpolateBilinearSlow; + +typedef struct _VipsInterpolateBilinearSlowClass { + VipsInterpolateClass parent_class; + +} VipsInterpolateBilinearSlowClass; + +GType vips_interpolate_bilinear_slow_get_type( void ); +VipsInterpolate *vips_interpolate_bilinear_slow_new( void ); + +/* Convenience: return a static fast bilinear_slow, so no need to free it. + */ +VipsInterpolate *vips_interpolate_bilinear_slow_static( void ); + /* Yafr class starts. */ @@ -192,6 +227,36 @@ VipsInterpolate *vips_interpolate_bilinear_static( void ); typedef struct _VipsInterpolateYafr { VipsInterpolate parent_object; + /* "sharpening" is a continuous method parameter which is + * proportional to the amount of "diagonal straightening" which the + * nonlinear correction part of the method may add to the underlying + * linear scheme. You may also think of it as a sharpening + * parameter: higher values correspond to more sharpening, and + * negative values lead to strange looking effects. + * + * The default value is sharpening = 29/32 when the scheme being + * "straightened" is Catmull-Rom---as is the case here. This value + * fixes key pixel values near the diagonal boundary between two + * monochrome regions (the diagonal boundary pixel values being set + * to the halfway colour). + * + * If resampling seems to add unwanted texture artifacts, push + * sharpening toward 0. It is not generally not recommended to set + * sharpening to a value larger than 4. + * + * Sharpening is halved because the .5 which has to do with the + * relative coordinates of the evaluation points (which has to do + * with .5*rite_width etc) is folded into the constant to save + * flops. Consequently, the largest recommended value of + * sharpening_over_two is 2=4/2. + * + * In order to simplify interfacing with users, the parameter which + * should be set by the user is normalized so that user_sharpening = + * 1 when sharpening is equal to the recommended value. Consistently + * with the above discussion, values of user_sharpening between 0 + * and about 3.625 give good results. + */ + double sharpening; } VipsInterpolateYafr; typedef struct _VipsInterpolateYafrClass { @@ -199,8 +264,10 @@ typedef struct _VipsInterpolateYafrClass { } VipsInterpolateYafrClass; +GType vips_interpolate_yafr_get_type( void ); VipsInterpolate *vips_interpolate_yafr_new( void ); -void vips_interpolate_yafr_set_thing( VipsInterpolateYafr *, double thing ); +void vips_interpolate_yafr_set_sharpening( VipsInterpolateYafr *, + double sharpening ); /* Convenience: return a static default yafr, so no need to free it. */ diff --git a/libsrc/mosaicing/gegl-sampler-yafr-smooth.c b/libsrc/mosaicing/gegl-sampler-yafr-smooth.c new file mode 100644 index 00000000..5f4b0b8e --- /dev/null +++ b/libsrc/mosaicing/gegl-sampler-yafr-smooth.c @@ -0,0 +1,686 @@ +/* This file is part of GEGL + * + * GEGL 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 3 of the + * License, or (at your option) any later version. + * + * GEGL 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 GEGL; if not, see + * . + * + * 2008 (c) Nicolas Robidoux (developer of Yet Another Fast + * Resampler). + * + * Acknowledgement: N. Robidoux's research on YAFR funded in part by + * an NSERC (National Science and Engineering Research Council of + * Canada) Discovery Grant. + */ + +#include +#include "gegl-types.h" +#include "gegl-buffer-private.h" +#include "gegl-sampler-yafr.h" + +#include + +#ifndef restrict +#ifdef __restrict +#define restrict __restrict +#else +#ifdef __restrict__ +#define restrict __restrict__ +#else +#define restrict +#endif +#endif +#endif + +#ifndef unlikely +#ifdef __builtin_expect +#define unlikely(x) __builtin_expect((x),0) +#else +#define unlikely(x) (x) +#endif +#endif + +enum +{ + PROP_0, + PROP_LAST +}; + +static void gegl_sampler_yafr_get ( GeglSampler *self, + const gdouble x, + const gdouble y, + void *output); + +static void set_property ( GObject *gobject, + guint property_id, + const GValue *value, + GParamSpec *pspec); + +static void get_property (GObject *gobject, + guint property_id, + GValue *value, + GParamSpec *pspec); + +G_DEFINE_TYPE (GeglSamplerYafr, gegl_sampler_yafr, GEGL_TYPE_SAMPLER) + +/* + * YAFR = Yet Another Fast Resampler + * + * Yet Another Fast Resampler is a nonlinear resampler which consists + * of a linear scheme (in this version, Catmull-Rom) plus a nonlinear + * sharpening correction the purpose of which is the straightening of + * diagonal interfaces between flat colour areas. + * + * Key properties: + * + * YAFR (smooth) is interpolatory: + * + * If asked for the value at the center of an input pixel, it will + * return the corresponding value, unchanged. + * + * YAFR (smooth) preserves local averages: + * + * The average of the reconstructed intensity surface over any region + * is the same as the average of the piecewise constant surface with + * values over pixel areas equal to the input pixel values (the + * "nearest neighbour" surface), except for a small amount of blur at + * the boundary of the region. More precicely: YAFR (smooth) is a box + * filtered exact area method. + * + * Main weaknesses of YAFR (smooth): + * + * Weakness 1: YAFR (smooth) improves on Catmull-Rom only for images + * with at least a little bit of smoothness. + * + * Weakness 2: Catmull-Rom introduces a lot of haloing. YAFR (smooth) + * is based on Catmull-Rom, and consequently it too introduces a lot + * of haloing. + * + * More details regarding Weakness 1: + * + * If a portion of the image is such that every pixel has immediate + * neighbours in the horizontal and vertical directions which have + * exactly the same pixel value, then YAFR (smooth) boils down to + * Catmull-Rom, and the computation of the correction is a waste. + * Extreme case: If all the pixels are either pure black or pure white + * in some region, as in some text images (more generally, if the + * region is "bichromatic"), then the YAFR (smooth) correction is 0 in + * the interior of the bichromatic region. + */ + +static void +gegl_sampler_yafr_class_init (GeglSamplerYafrClass *klass) +{ + GeglSamplerClass *sampler_class = GEGL_SAMPLER_CLASS (klass); + GObjectClass *object_class = G_OBJECT_CLASS (klass); + + object_class->set_property = set_property; + object_class->get_property = get_property; + + sampler_class->get = gegl_sampler_yafr_get; + } + +static void +gegl_sampler_yafr_init (GeglSamplerYafr *self) +{ + /* + * The computation stencil is 4x4, and sticks out one column to the + * left and one row above the requested integer position: + */ + GEGL_SAMPLER (self)->context_rect = (GeglRectangle){-1,-1,4,4}; + + GEGL_SAMPLER (self)->interpolate_format = babl_format ("RaGaBaA float"); +} + +static inline gfloat +catrom_yafr (const gfloat cardinal_one, + const gfloat cardinal_two, + const gfloat cardinal_thr, + const gfloat cardinal_fou, + const gfloat cardinal_uno, + const gfloat cardinal_dos, + const gfloat cardinal_tre, + const gfloat cardinal_qua, + const gfloat left_width_times_up__height_times_rite_width, + const gfloat left_width_times_dow_height_times_rite_width, + const gfloat left_width_times_up__height_times_dow_height, + const gfloat rite_width_times_up__height_times_dow_height, + const gfloat* restrict this_channels_uno_one_bptr) +{ + /* + * "sharpening" is a continuous method parameter which is + * proportional to the amount of "diagonal straightening" which the + * nonlinear correction part of the method may add to the underlying + * linear scheme. You may also think of it as a sharpening + * parameter: higher values correspond to more sharpening, and + * negative values lead to strange looking effects. + * + * The default value is sharpening = 29/32 when the scheme being + * "straightened" is Catmull-Rom---as is the case here. This value + * fixes key pixel values near the diagonal boundary between two + * monochrome regions (the diagonal boundary pixel values being set + * to the halfway colour). + * + * If resampling seems to add unwanted texture artifacts, push + * sharpening toward 0. It is not generally not recommended to set + * sharpening to a value larger than 4. + * + * Sharpening is halved because the .5 which has to do with the + * relative coordinates of the evaluation points (which has to do + * with .5*rite_width etc) is folded into the constant to save + * flops. Consequently, the largest recommended value of + * sharpening_over_two is 2=4/2. + * + * In order to simplify interfacing with users, the parameter which + * should be set by the user is normalized so that user_sharpening = + * 1 when sharpening is equal to the recommended value. Consistently + * with the above discussion, values of user_sharpening between 0 + * and about 3.625 give good results. + */ + const gfloat user_sharpening = 1.f; + const gfloat sharpening_over_two = user_sharpening * 0.453125f; + + /* + * The input pixel values are described by the following stencil. + * Spanish abbreviations are used to label positions from top to + * bottom, English ones to label positions from left to right,: + * + * (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 + * + * (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 + */ + + /* + * Load the useful pixel values for the channel under + * consideration. The this_channels_uno_one_bptr pointer is assumed + * to point to uno_one when catrom_yafr is entered. + */ + const gint channels = 4; + const gint pixels_per_buffer_row = 64; + const gfloat uno_one = + this_channels_uno_one_bptr[ 0 ]; + const gfloat uno_two = + this_channels_uno_one_bptr[ channels ]; + const gfloat uno_thr = + this_channels_uno_one_bptr[ 2 * channels ]; + const gfloat uno_fou = + this_channels_uno_one_bptr[ 3 * channels ]; + + const gfloat dos_one = + this_channels_uno_one_bptr[ pixels_per_buffer_row * channels ]; + const gfloat dos_two = + this_channels_uno_one_bptr[ ( 1 + pixels_per_buffer_row ) * channels ]; + const gfloat dos_thr = + this_channels_uno_one_bptr[ ( 2 + pixels_per_buffer_row ) * channels ]; + const gfloat dos_fou = + this_channels_uno_one_bptr[ ( 3 + pixels_per_buffer_row ) * channels ]; + + const gfloat tre_one = + this_channels_uno_one_bptr[ 2 * pixels_per_buffer_row * channels ]; + const gfloat tre_two = + this_channels_uno_one_bptr[ ( 1 + 2 * pixels_per_buffer_row ) * channels ]; + const gfloat tre_thr = + this_channels_uno_one_bptr[ ( 2 + 2 * pixels_per_buffer_row ) * channels ]; + const gfloat tre_fou = + this_channels_uno_one_bptr[ ( 3 + 2 * pixels_per_buffer_row ) * channels ]; + + const gfloat qua_one = + this_channels_uno_one_bptr[ 3 * pixels_per_buffer_row * channels ]; + const gfloat qua_two = + this_channels_uno_one_bptr[ ( 1 + 3 * pixels_per_buffer_row ) * channels ]; + const gfloat qua_thr = + this_channels_uno_one_bptr[ ( 2 + 3 * pixels_per_buffer_row ) * channels ]; + const gfloat qua_fou = + this_channels_uno_one_bptr[ ( 3 + 3 * pixels_per_buffer_row ) * channels ]; + + /* + * Computation of the YAFR correction: + * + * Basically, if two consecutive pixel value differences have the + * same sign, the smallest one (in absolute value) is taken to be + * the corresponding slope. If they don't have the same sign, the + * corresponding slope is set to 0. + * + * Four such pairs (vertical and horizontal) of slopes need to be + * computed, one pair for each of the pixels which potentially + * overlap the unit area centered at the interpolation point. + */ + /* + * Beginning of the computation of the "up" horizontal slopes: + */ + const gfloat prem__up = dos_two - dos_one; + const gfloat deux__up = dos_thr - dos_two; + const gfloat troi__up = dos_fou - dos_thr; + /* + * "down" horizontal slopes: + */ + const gfloat prem_dow = tre_two - tre_one; + const gfloat deux_dow = tre_thr - tre_two; + const gfloat troi_dow = tre_fou - tre_thr; + /* + * "left" vertical slopes: + */ + const gfloat prem_left = dos_two - uno_two; + const gfloat deux_left = tre_two - dos_two; + const gfloat troi_left = qua_two - tre_two; + /* + * "right" vertical slopes: + */ + const gfloat prem_rite = dos_thr - uno_thr; + const gfloat deux_rite = tre_thr - dos_thr; + const gfloat troi_rite = qua_thr - tre_thr; + + /* + * Back to "up": + */ + const gfloat prem__up_squared = prem__up * prem__up; + const gfloat deux__up_squared = deux__up * deux__up; + const gfloat troi__up_squared = troi__up * troi__up; + /* + * Back to "down": + */ + const gfloat prem_dow_squared = prem_dow * prem_dow; + const gfloat deux_dow_squared = deux_dow * deux_dow; + const gfloat troi_dow_squared = troi_dow * troi_dow; + /* + * Back to "left": + */ + const gfloat prem_left_squared = prem_left * prem_left; + const gfloat deux_left_squared = deux_left * deux_left; + const gfloat troi_left_squared = troi_left * troi_left; + /* + * Back to "right": + */ + const gfloat prem_rite_squared = prem_rite * prem_rite; + const gfloat deux_rite_squared = deux_rite * deux_rite; + const gfloat troi_rite_squared = troi_rite * troi_rite; + + /* + * "up": + */ + const gfloat prem__up_times_deux__up = prem__up * deux__up; + const gfloat deux__up_times_troi__up = deux__up * troi__up; + /* + * "down": + */ + const gfloat prem_dow_times_deux_dow = prem_dow * deux_dow; + const gfloat deux_dow_times_troi_dow = deux_dow * troi_dow; + /* + * "left": + */ + const gfloat prem_left_times_deux_left = prem_left * deux_left; + const gfloat deux_left_times_troi_left = deux_left * troi_left; + /* + * "right": + */ + const gfloat prem_rite_times_deux_rite = prem_rite * deux_rite; + const gfloat deux_rite_times_troi_rite = deux_rite * troi_rite; + + /* + * Branching parts of the computation of the YAFR correction (could + * be unbranched using arithmetic branching and C99 math intrinsics, + * although the compiler may be smart enough to remove the branching + * on its own): + */ + /* + * "up": + */ + const gfloat prem__up_vs_deux__up = + prem__up_squared < deux__up_squared ? prem__up : deux__up; + const gfloat deux__up_vs_troi__up = + deux__up_squared < troi__up_squared ? deux__up : troi__up; + /* + * "down": + */ + const gfloat prem_dow_vs_deux_dow = + prem_dow_squared < deux_dow_squared ? prem_dow : deux_dow; + const gfloat deux_dow_vs_troi_dow = + deux_dow_squared < troi_dow_squared ? deux_dow : troi_dow; + /* + * "left": + */ + const gfloat prem_left_vs_deux_left = + prem_left_squared < deux_left_squared ? prem_left : deux_left; + const gfloat deux_left_vs_troi_left = + deux_left_squared < troi_left_squared ? deux_left : troi_left; + /* + * "right": + */ + const gfloat prem_rite_vs_deux_rite = + prem_rite_squared < deux_rite_squared ? prem_rite : deux_rite; + const gfloat deux_rite_vs_troi_rite = + deux_rite_squared < troi_rite_squared ? deux_rite : troi_rite; + /* + * The YAFR correction computation will resume after the computation + * of the Catmull-Rom baseline. + */ + + /* + * Catmull-Rom baseline contribution: + */ + const gfloat catmull_rom = + cardinal_uno * + ( + cardinal_one * uno_one + + + cardinal_two * uno_two + + + cardinal_thr * uno_thr + + + cardinal_fou * uno_fou + ) + + + cardinal_dos * + ( + cardinal_one * dos_one + + + cardinal_two * dos_two + + + cardinal_thr * dos_thr + + + cardinal_fou * dos_fou + ) + + + cardinal_tre * + ( + cardinal_one * tre_one + + + cardinal_two * tre_two + + + cardinal_thr * tre_thr + + + cardinal_fou * tre_fou + ) + + + cardinal_qua * + ( + cardinal_one * qua_one + + + cardinal_two * qua_two + + + cardinal_thr * qua_thr + + + cardinal_fou * qua_fou + ); + + /* + * Computation of the YAFR slopes. + */ + /* + * "up": + */ + const gfloat mx_left__up = + prem__up_times_deux__up < 0.f ? 0.f : prem__up_vs_deux__up; + const gfloat mx_rite__up = + deux__up_times_troi__up < 0.f ? 0.f : deux__up_vs_troi__up; + /* + * "down": + */ + const gfloat mx_left_dow = + prem_dow_times_deux_dow < 0.f ? 0.f : prem_dow_vs_deux_dow; + const gfloat mx_rite_dow = + deux_dow_times_troi_dow < 0.f ? 0.f : deux_dow_vs_troi_dow; + /* + * "left": + */ + const gfloat my_left__up = + prem_left_times_deux_left < 0.f ? 0.f : prem_left_vs_deux_left; + const gfloat my_left_dow = + deux_left_times_troi_left < 0.f ? 0.f : deux_left_vs_troi_left; + /* + * "down": + */ + const gfloat my_rite__up = + prem_rite_times_deux_rite < 0.f ? 0.f : prem_rite_vs_deux_rite; + const gfloat my_rite_dow = + deux_rite_times_troi_rite < 0.f ? 0.f : deux_rite_vs_troi_rite; + + /* + * Assemble the unweighted YAFR correction: + */ + const gfloat unweighted_yafr_correction = + left_width_times_up__height_times_rite_width + * + ( mx_left__up - mx_rite__up ) + + + left_width_times_dow_height_times_rite_width + * + ( mx_left_dow - mx_rite_dow ) + + + left_width_times_up__height_times_dow_height + * + ( my_left__up - my_left_dow ) + + + rite_width_times_up__height_times_dow_height + * + ( my_rite__up - my_rite_dow ); + + /* + * Add the Catmull-Rom baseline and the weighted YAFR correction: + */ + const gfloat newval = + sharpening_over_two * unweighted_yafr_correction + catmull_rom; + + return newval; +} + +static void +gegl_sampler_yafr_get ( GeglSampler *self, + const gdouble x, + const gdouble y, + void *output) +{ + /* + * Note: The computation is structured to foster software + * pipelining. + */ + + /* + * x is understood to increase from left to right, y, from top to + * bottom. Consequently, ix and iy are the indices of the pixel + * located at or to the left, and at or above. the sampling point. + * + * floor is used to make sure that the transition through 0 is + * smooth. If it is known that negative x and y will never be used, + * cast (which truncates) could be used instead. + */ + const gint ix = floorf (x); + const gint iy = floorf (y); + + /* + * Pointer to enlarged input stencil values: + */ + const gfloat* restrict sampler_bptr = gegl_sampler_get_ptr (self, ix, iy); + + /* + * Each (channel's) output pixel value is obtained by combining four + * "pieces," each piece corresponding to the set of points which are + * closest to the four pixels closest to the (x,y) position, pixel + * positions which have coordinates and labels as follows: + * + * (ix,iy) (ix+1,iy) + * =left__up =rite__up + * + * <- (x,y) is somewhere in the convex hull + * + * (ix,iy+1) (ix+1,iy+1) + * =left_dow =rite_dow + */ + /* + * rite_width is the width of the overlaps of the unit averaging box + * (which is centered at the position where an interpolated value is + * desired), with the closest unit pixel areas to the right. + * + * left_width is the width of the overlaps of the unit averaging box + * (which is centered at the position where an interpolated value is + * desired), with the closest unit pixel areas to the left. + */ + const gfloat rite_width = x - ix; + const gfloat dow_height = y - iy; + const gfloat left_width = 1.f - rite_width; + const gfloat up__height = 1.f - dow_height; + /* + * .5*rite_width is the x-coordinate of the center of the overlap of + * the averaging box with the left pixel areas, relative to the + * position of the centers of the left pixels. + * + * -.5*left_width is the x-coordinate ... right pixel areas, + * relative to ... the right pixels. + * + * .5*dow_height is the y-coordinate of the center of the overlap + * of the averaging box with the up pixel areas, relative to the + * position of the centers of the up pixels. + * + * -.5*up__height is the y-coordinate ... down pixel areas, relative + * to ... the down pixels. + */ + const gfloat left_width_times_rite_width = left_width * rite_width; + const gfloat up__height_times_dow_height = up__height * dow_height; + + const gfloat cardinal_two = + left_width_times_rite_width * ( -1.5f * rite_width + 1.f ) + + left_width; + const gfloat cardinal_dos = + up__height_times_dow_height * ( -1.5f * dow_height + 1.f ) + + up__height; + + const gfloat minus_half_left_width_times_rite_width = + -.5f * left_width_times_rite_width; + const gfloat minus_half_up__height_times_dow_height = + -.5f * up__height_times_dow_height; + + const gfloat left_width_times_up__height_times_rite_width = + left_width_times_rite_width * up__height; + const gfloat left_width_times_dow_height_times_rite_width = + left_width_times_rite_width * dow_height; + const gfloat left_width_times_up__height_times_dow_height = + up__height_times_dow_height * left_width; + const gfloat rite_width_times_up__height_times_dow_height = + up__height_times_dow_height * rite_width; + + const gfloat cardinal_one = + minus_half_left_width_times_rite_width * left_width; + const gfloat cardinal_uno = + minus_half_up__height_times_dow_height * up__height; + + const gfloat cardinal_fou = + minus_half_left_width_times_rite_width * rite_width; + const gfloat cardinal_qua = + minus_half_up__height_times_dow_height * dow_height; + + const gfloat cardinal_thr = + 1.f - ( minus_half_left_width_times_rite_width + cardinal_two ); + const gfloat cardinal_tre = + 1.f - ( minus_half_up__height_times_dow_height + cardinal_dos ); + + /* + * The newval array will contain the four (one per channel) + * computed resampled values: + */ + gfloat newval[4]; + + /* + * Set the tile pointer to the first relevant value. Since the + * pointer initially points to dos_two, we need to rewind it one + * tile row, then go back one additional pixel. + */ + const gint channels = 4; + const gint pixels_per_buffer_row = 64; + sampler_bptr -= ( pixels_per_buffer_row + 1 ) * channels; + + newval[0] = catrom_yafr (cardinal_one, + cardinal_two, + cardinal_thr, + cardinal_fou, + cardinal_uno, + cardinal_dos, + cardinal_tre, + cardinal_qua, + left_width_times_up__height_times_rite_width, + left_width_times_dow_height_times_rite_width, + left_width_times_up__height_times_dow_height, + rite_width_times_up__height_times_dow_height, + sampler_bptr++); + newval[1] = catrom_yafr (cardinal_one, + cardinal_two, + cardinal_thr, + cardinal_fou, + cardinal_uno, + cardinal_dos, + cardinal_tre, + cardinal_qua, + left_width_times_up__height_times_rite_width, + left_width_times_dow_height_times_rite_width, + left_width_times_up__height_times_dow_height, + rite_width_times_up__height_times_dow_height, + sampler_bptr++); + newval[2] = catrom_yafr (cardinal_one, + cardinal_two, + cardinal_thr, + cardinal_fou, + cardinal_uno, + cardinal_dos, + cardinal_tre, + cardinal_qua, + left_width_times_up__height_times_rite_width, + left_width_times_dow_height_times_rite_width, + left_width_times_up__height_times_dow_height, + rite_width_times_up__height_times_dow_height, + sampler_bptr++); + newval[3] = catrom_yafr (cardinal_one, + cardinal_two, + cardinal_thr, + cardinal_fou, + cardinal_uno, + cardinal_dos, + cardinal_tre, + cardinal_qua, + left_width_times_up__height_times_rite_width, + left_width_times_dow_height_times_rite_width, + left_width_times_up__height_times_dow_height, + rite_width_times_up__height_times_dow_height, + sampler_bptr); + + /* + * Ship out newval: + */ + babl_process (babl_fish (self->interpolate_format, self->format), + newval, + output, + 1); +} + +static void +set_property ( GObject *gobject, + guint property_id, + const GValue *value, + GParamSpec *pspec) +{ + G_OBJECT_WARN_INVALID_PROPERTY_ID (gobject, property_id, pspec); +} + +static void +get_property (GObject *gobject, + guint property_id, + GValue *value, + GParamSpec *pspec) +{ + G_OBJECT_WARN_INVALID_PROPERTY_ID (gobject, property_id, pspec); +} diff --git a/libsrc/mosaicing/interpolate.c b/libsrc/mosaicing/interpolate.c index ff267598..c718b631 100644 --- a/libsrc/mosaicing/interpolate.c +++ b/libsrc/mosaicing/interpolate.c @@ -39,6 +39,7 @@ #include #include +#include #include #include @@ -54,6 +55,8 @@ static VipsObjectClass *vips_interpolate_parent_class = NULL; static VipsObjectClass *vips_interpolate_nearest_parent_class = NULL; static VipsObjectClass *vips_interpolate_bilinear_parent_class = NULL; +static VipsObjectClass *vips_interpolate_bilinear_slow_parent_class = NULL; +static VipsObjectClass *vips_interpolate_yafr_parent_class = NULL; #ifdef DEBUG static void @@ -154,18 +157,6 @@ vips_interpolate_get_window_size( VipsInterpolate *interpolate ) /* VipsInterpolateNearest class */ -#ifdef DEBUG -static void -vips_interpolate_nearest_finalize( GObject *gobject ) -{ - printf( "vips_interpolate_nearest_finalize: " ); - vips_object_print( VIPS_OBJECT( gobject ) ); - - G_OBJECT_CLASS( vips_interpolate_nearest_parent_class )-> - finalize( gobject ); -} -#endif /*DEBUG*/ - static void vips_interpolate_nearest_interpolate( VipsInterpolate *interpolate, REGION *out, REGION *in, @@ -197,18 +188,12 @@ vips_interpolate_nearest_interpolate( VipsInterpolate *interpolate, static void vips_interpolate_nearest_class_init( VipsInterpolateNearestClass *class ) { -#ifdef DEBUG - GObjectClass *gobject_class = G_OBJECT_CLASS( class ); -#endif /*DEBUG*/ VipsInterpolateClass *interpolate_class = VIPS_INTERPOLATE_CLASS( class ); vips_interpolate_nearest_parent_class = g_type_class_peek_parent( class ); -#ifdef DEBUG - gobject_class->finalize = vips_interpolate_nearest_finalize; -#endif /*DEBUG*/ interpolate_class->interpolate = vips_interpolate_nearest_interpolate; interpolate_class->window_size = 1; } @@ -220,6 +205,7 @@ vips_interpolate_nearest_init( VipsInterpolateNearest *nearest ) printf( "vips_interpolate_nearest_init: " ); vips_object_print( VIPS_OBJECT( nearest ) ); #endif /*DEBUG*/ + } GType @@ -275,18 +261,6 @@ vips_interpolate_nearest_static( void ) * p3 p4 */ -#ifdef DEBUG -static void -vips_interpolate_bilinear_finalize( GObject *gobject ) -{ - printf( "vips_interpolate_bilinear_finalize: " ); - vips_object_print( VIPS_OBJECT( gobject ) ); - - G_OBJECT_CLASS( vips_interpolate_bilinear_parent_class )-> - finalize( gobject ); -} -#endif /*DEBUG*/ - /* Interpolate a section ... int8/16 types. */ #define BILINEAR_INT( TYPE ) { \ @@ -337,21 +311,6 @@ vips_interpolate_bilinear_finalize( GObject *gobject ) c3 * tp3[z] + c4 * tp4[z]; \ } -/* Interpolate a pel ... don't use the pre-calcuated matricies. - */ -#define BILINEAR_SLOW( TYPE ) { \ - TYPE *tq = (TYPE *) q; \ - \ - const TYPE *tp1 = (TYPE *) p1; \ - const TYPE *tp2 = (TYPE *) p2; \ - const TYPE *tp3 = (TYPE *) p3; \ - const TYPE *tp4 = (TYPE *) p4; \ - \ - for( z = 0; z < b; z++ ) \ - tq[z] = c1 * tp1[z] + c2 * tp2[z] + \ - c3 * tp3[z] + c4 * tp4[z]; \ -} - /* Expand for band types. with a fixed-point interpolator and a float * interpolator. */ @@ -385,109 +344,42 @@ vips_interpolate_bilinear_interpolate( VipsInterpolate *interpolate, const int ps = IM_IMAGE_SIZEOF_PEL( in->im ); const int ls = IM_REGION_LSKIP( in ); const int b = in->im->Bands; - int z; + + /* Subtract 0.5 to centre the bilinear. + */ + const double cx = in_x - 0.5; + const double cy = in_y - 0.5; + + /* Now go to scaled int. + */ + const double sx = cx * VIPS_TRANSFORM_SCALE; + const double sy = cy * VIPS_TRANSFORM_SCALE; + const int sxi = FLOOR( sx ); + const int syi = FLOOR( sy ); + + /* Get index into interpolation table and unscaled integer + * position. + */ + const int xi = sxi & (VIPS_TRANSFORM_SCALE - 1); + const int yi = syi & (VIPS_TRANSFORM_SCALE - 1); + const int in_x_int = sxi >> VIPS_TRANSFORM_SHIFT; + const int in_y_int = syi >> VIPS_TRANSFORM_SHIFT; + + const PEL *p1 = (PEL *) IM_REGION_ADDR( in, in_x_int, in_y_int ); + const PEL *p2 = p1 + ps; + const PEL *p3 = p1 + ls; + const PEL *p4 = p3 + ps; PEL *q = (PEL *) IM_REGION_ADDR( out, out_x, out_y ); + int z; - if( bilinear->slow ) { - /* Subtract 0.5 to centre the bilinear. - */ - const double cx = in_x - 0.5; - const double cy = in_y - 0.5; - - /* Top left corner we interpolate from. - */ - const int xi = FLOOR( cx ); - const int yi = FLOOR( cy ); - - /* Fractional part. - */ - const double X = cx - xi; - const double Y = cy - yi; - - /* Residual. - */ - const double Xd = 1.0 - X; - const double Yd = 1.0 - Y; - - /* Weights. - */ - const double c1 = Xd * Yd; - const double c2 = X * Yd; - const double c3 = Xd * Y; - const double c4 = X * Y; - - const PEL *p1 = (PEL *) IM_REGION_ADDR( in, xi, yi ); - const PEL *p2 = p1 + ps; - const PEL *p3 = p1 + ls; - const PEL *p4 = p3 + ps; - - SWITCH_INTERPOLATE( in->im->BandFmt, - BILINEAR_SLOW, BILINEAR_SLOW ); - } - else { - /* Subtract 0.5 to centre the bilinear. - */ - const double cx = in_x - 0.5; - const double cy = in_y - 0.5; - - /* Now go to scaled int. - */ - const double sx = cx * VIPS_TRANSFORM_SCALE; - const double sy = cy * VIPS_TRANSFORM_SCALE; - const int sxi = FLOOR( sx ); - const int syi = FLOOR( sy ); - - /* Get index into interpolation table and unscaled integer - * position. - */ - const int xi = sxi & (VIPS_TRANSFORM_SCALE - 1); - const int yi = syi & (VIPS_TRANSFORM_SCALE - 1); - const int in_x_int = sxi >> VIPS_TRANSFORM_SHIFT; - const int in_y_int = syi >> VIPS_TRANSFORM_SHIFT; - - const PEL *p1 = (PEL *) - IM_REGION_ADDR( in, in_x_int, in_y_int ); - const PEL *p2 = p1 + ps; - const PEL *p3 = p1 + ls; - const PEL *p4 = p3 + ps; - -{ - unsigned char *tq = (unsigned char *) q; - - const int m1 = class->matrix_int[xi][0]; - const int m2 = class->matrix_int[xi][1]; - const int m3 = class->matrix_int[yi][0]; - const int m4 = class->matrix_int[yi][1]; - - const int c1 = (m3 * m1) >> VIPS_INTERPOLATE_SHIFT; - const int c2 = (m3 * m2) >> VIPS_INTERPOLATE_SHIFT; - const int c3 = (m4 * m1) >> VIPS_INTERPOLATE_SHIFT; - const int c4 = (m4 * m2) >> VIPS_INTERPOLATE_SHIFT; - - const unsigned char *tp1 = (unsigned char *) p1; - const unsigned char *tp2 = (unsigned char *) p2; - const unsigned char *tp3 = (unsigned char *) p3; - const unsigned char *tp4 = (unsigned char *) p4; - - for( z = 0; z < b; z++ ) - tq[z] = (c1 * tp1[z] + c2 * tp2[z] + - c3 * tp3[z] + c4 * tp4[z]) >> VIPS_INTERPOLATE_SHIFT; -} - - /* - SWITCH_INTERPOLATE( in->im->BandFmt, - BILINEAR_INT, BILINEAR_FLOAT ); - */ - } + SWITCH_INTERPOLATE( in->im->BandFmt, + BILINEAR_INT, BILINEAR_FLOAT ); } static void vips_interpolate_bilinear_class_init( VipsInterpolateBilinearClass *class ) { -#ifdef DEBUG - GObjectClass *gobject_class = G_OBJECT_CLASS( class ); -#endif /*DEBUG*/ VipsInterpolateClass *interpolate_class = (VipsInterpolateClass *) class; int x; @@ -495,9 +387,6 @@ vips_interpolate_bilinear_class_init( VipsInterpolateBilinearClass *class ) vips_interpolate_bilinear_parent_class = g_type_class_peek_parent( class ); -#ifdef DEBUG - gobject_class->finalize = vips_interpolate_bilinear_finalize; -#endif /*DEBUG*/ interpolate_class->interpolate = vips_interpolate_bilinear_interpolate; interpolate_class->window_size = 2; @@ -525,7 +414,6 @@ vips_interpolate_bilinear_init( VipsInterpolateBilinear *bilinear ) vips_object_print( VIPS_OBJECT( bilinear ) ); #endif /*DEBUG*/ - bilinear->slow = FALSE; } GType @@ -553,13 +441,6 @@ vips_interpolate_bilinear_get_type( void ) return( type ); } -void -vips_interpolate_bilinear_set_slow( VipsInterpolateBilinear *bilinear, - gboolean slow ) -{ - bilinear->slow = slow; -} - VipsInterpolate * vips_interpolate_bilinear_new( void ) { @@ -580,3 +461,777 @@ vips_interpolate_bilinear_static( void ) return( interpolate ); } +/* VipsInterpolateBilinearSlow class + */ + +/* Slow mode is really just for testing ... it doesn't use the pre-calculated + * interpolation factors or the fixed-point arithmetic. + */ + +/* in this class, name vars in the 2x2 grid as eg. + * p1 p2 + * p3 p4 + */ + +/* Interpolate a pel ... don't use the pre-calcuated matricies. + */ +#define BILINEAR_SLOW( TYPE ) { \ + TYPE *tq = (TYPE *) q; \ + \ + const TYPE *tp1 = (TYPE *) p1; \ + const TYPE *tp2 = (TYPE *) p2; \ + const TYPE *tp3 = (TYPE *) p3; \ + const TYPE *tp4 = (TYPE *) p4; \ + \ + for( z = 0; z < b; z++ ) \ + tq[z] = c1 * tp1[z] + c2 * tp2[z] + \ + c3 * tp3[z] + c4 * tp4[z]; \ +} + +static void +vips_interpolate_bilinear_slow_interpolate( VipsInterpolate *interpolate, + REGION *out, REGION *in, + int out_x, int out_y, double in_x, double in_y ) +{ + VipsInterpolateBilinearSlow *bilinear_slow = + VIPS_INTERPOLATE_BILINEAR_SLOW( interpolate ); + VipsInterpolateBilinearSlowClass *class = + VIPS_INTERPOLATE_BILINEAR_SLOW_GET_CLASS( interpolate ); + + /* Pel size and line size. + */ + const int ps = IM_IMAGE_SIZEOF_PEL( in->im ); + const int ls = IM_REGION_LSKIP( in ); + const int b = in->im->Bands; + + /* Subtract 0.5 to centre the bilinear. + */ + const double cx = in_x - 0.5; + const double cy = in_y - 0.5; + + /* Top left corner we interpolate from. + */ + const int xi = FLOOR( cx ); + const int yi = FLOOR( cy ); + + /* Fractional part. + */ + const double X = cx - xi; + const double Y = cy - yi; + + /* Residual. + */ + const double Xd = 1.0 - X; + const double Yd = 1.0 - Y; + + /* Weights. + */ + const double c1 = Xd * Yd; + const double c2 = X * Yd; + const double c3 = Xd * Y; + const double c4 = X * Y; + + const PEL *p1 = (PEL *) IM_REGION_ADDR( in, xi, yi ); + const PEL *p2 = p1 + ps; + const PEL *p3 = p1 + ls; + const PEL *p4 = p3 + ps; + + PEL *q = (PEL *) IM_REGION_ADDR( out, out_x, out_y ); + int z; + + SWITCH_INTERPOLATE( in->im->BandFmt, + BILINEAR_SLOW, BILINEAR_SLOW ); +} + +static void +vips_interpolate_bilinear_slow_class_init( VipsInterpolateBilinearSlowClass *class ) +{ + VipsInterpolateClass *interpolate_class = + (VipsInterpolateClass *) class; + + vips_interpolate_bilinear_slow_parent_class = + g_type_class_peek_parent( class ); + + interpolate_class->interpolate = + vips_interpolate_bilinear_slow_interpolate; + interpolate_class->window_size = 2; +} + +static void +vips_interpolate_bilinear_slow_init( + VipsInterpolateBilinearSlow *bilinear_slow ) +{ +#ifdef DEBUG + printf( "vips_interpolate_bilinear_slow_init: " ); + vips_object_print( VIPS_OBJECT( bilinear_slow ) ); +#endif /*DEBUG*/ + +} + +GType +vips_interpolate_bilinear_slow_get_type( void ) +{ + static GType type = 0; + + if( !type ) { + static const GTypeInfo info = { + sizeof( VipsInterpolateBilinearSlowClass ), + NULL, /* base_init */ + NULL, /* base_finalize */ + (GClassInitFunc) + vips_interpolate_bilinear_slow_class_init, + NULL, /* class_finalize */ + NULL, /* class_data */ + sizeof( VipsInterpolateBilinearSlow ), + 32, /* n_preallocs */ + (GInstanceInitFunc) vips_interpolate_bilinear_slow_init, + }; + + type = g_type_register_static( VIPS_TYPE_INTERPOLATE, + "VipsInterpolateBilinearSlow", &info, 0 ); + } + + return( type ); +} + +VipsInterpolate * +vips_interpolate_bilinear_slow_new( void ) +{ + return( VIPS_INTERPOLATE( g_object_new( + VIPS_TYPE_INTERPOLATE_BILINEAR_SLOW, NULL ) ) ); +} + +/* Convenience: return a static bilinear_slow you don't need to free. + */ +VipsInterpolate * +vips_interpolate_bilinear_slow_static( void ) +{ + static VipsInterpolate *interpolate = NULL; + + if( !interpolate ) + interpolate = vips_interpolate_bilinear_slow_new(); + + return( interpolate ); +} + +/* VipsInterpolateYafr class + */ + +/* Copy-paste of gegl-sampler-yafr-smooth.c starts + */ + +#ifndef restrict +#ifdef __restrict +#define restrict __restrict +#else +#ifdef __restrict__ +#define restrict __restrict__ +#else +#define restrict +#endif +#endif +#endif + +#ifndef unlikely +#ifdef __builtin_expect +#define unlikely(x) __builtin_expect((x),0) +#else +#define unlikely(x) (x) +#endif +#endif + +/* + * YAFR = Yet Another Fast Resampler + * + * Yet Another Fast Resampler is a nonlinear resampler which consists + * of a linear scheme (in this version, Catmull-Rom) plus a nonlinear + * sharpening correction the purpose of which is the straightening of + * diagonal interfaces between flat colour areas. + * + * Key properties: + * + * YAFR (smooth) is interpolatory: + * + * If asked for the value at the center of an input pixel, it will + * return the corresponding value, unchanged. + * + * YAFR (smooth) preserves local averages: + * + * The average of the reconstructed intensity surface over any region + * is the same as the average of the piecewise constant surface with + * values over pixel areas equal to the input pixel values (the + * "nearest neighbour" surface), except for a small amount of blur at + * the boundary of the region. More precicely: YAFR (smooth) is a box + * filtered exact area method. + * + * Main weaknesses of YAFR (smooth): + * + * Weakness 1: YAFR (smooth) improves on Catmull-Rom only for images + * with at least a little bit of smoothness. + * + * Weakness 2: Catmull-Rom introduces a lot of haloing. YAFR (smooth) + * is based on Catmull-Rom, and consequently it too introduces a lot + * of haloing. + * + * More details regarding Weakness 1: + * + * If a portion of the image is such that every pixel has immediate + * neighbours in the horizontal and vertical directions which have + * exactly the same pixel value, then YAFR (smooth) boils down to + * Catmull-Rom, and the computation of the correction is a waste. + * Extreme case: If all the pixels are either pure black or pure white + * in some region, as in some text images (more generally, if the + * region is "bichromatic"), then the YAFR (smooth) correction is 0 in + * the interior of the bichromatic region. + */ + +/* Pointers to write to / read from, how much to add to move right a pixel, + * how much to add to move down a line. + */ + +static inline void +catrom_yafr (float* restrict out, const float* restrict in, + const int channels, + const int pixels_per_buffer_row, + const float sharpening, + + const float cardinal_one, + const float cardinal_two, + const float cardinal_thr, + const float cardinal_fou, + const float cardinal_uno, + const float cardinal_dos, + const float cardinal_tre, + const float cardinal_qua, + const float left_width_times_up__height_times_rite_width, + const float left_width_times_dow_height_times_rite_width, + const float left_width_times_up__height_times_dow_height, + const float rite_width_times_up__height_times_dow_height) +{ + + /* "sharpening" is a continuous method parameter which is + * proportional to the amount of "diagonal straightening" which the + * nonlinear correction part of the method may add to the underlying + * linear scheme. You may also think of it as a sharpening + * parameter: higher values correspond to more sharpening, and + * negative values lead to strange looking effects. + * + * The default value is sharpening = 29/32 when the scheme being + * "straightened" is Catmull-Rom---as is the case here. This value + * fixes key pixel values near the diagonal boundary between two + * monochrome regions (the diagonal boundary pixel values being set + * to the halfway colour). + * + * If resampling seems to add unwanted texture artifacts, push + * sharpening toward 0. It is not generally not recommended to set + * sharpening to a value larger than 4. + * + * Sharpening is halved because the .5 which has to do with the + * relative coordinates of the evaluation points (which has to do + * with .5*rite_width etc) is folded into the constant to save + * flops. Consequently, the largest recommended value of + * sharpening_over_two is 2=4/2. + * + * In order to simplify interfacing with users, the parameter which + * should be set by the user is normalized so that user_sharpening = + * 1 when sharpening is equal to the recommended value. Consistently + * with the above discussion, values of user_sharpening between 0 + * and about 3.625 give good results. + */ + + const float sharpening_over_two = sharpening * 0.453125f; + + /* + * The input pixel values are described by the following stencil. + * Spanish abbreviations are used to label positions from top to + * bottom, English ones to label positions from left to right,: + * + * (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 + * + * (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 + */ + + /* + * Load the useful pixel values for the channel under + * consideration. The in pointer is assumed + * to point to uno_one when catrom_yafr is entered. + */ + const float uno_one = in[ 0 ]; + const float uno_two = in[ channels ]; + const float uno_thr = in[ 2 * channels ]; + const float uno_fou = in[ 3 * channels ]; + + const float dos_one = in[ pixels_per_buffer_row * channels ]; + const float dos_two = in[ ( 1 + pixels_per_buffer_row ) * channels ]; + const float dos_thr = in[ ( 2 + pixels_per_buffer_row ) * channels ]; + const float dos_fou = in[ ( 3 + pixels_per_buffer_row ) * channels ]; + + const float tre_one = in[ 2 * pixels_per_buffer_row * channels ]; + const float tre_two = in[ ( 1 + 2 * pixels_per_buffer_row ) * channels ]; + const float tre_thr = in[ ( 2 + 2 * pixels_per_buffer_row ) * channels ]; + const float tre_fou = in[ ( 3 + 2 * pixels_per_buffer_row ) * channels ]; + + const float qua_one = in[ 3 * pixels_per_buffer_row * channels ]; + const float qua_two = in[ ( 1 + 3 * pixels_per_buffer_row ) * channels ]; + const float qua_thr = in[ ( 2 + 3 * pixels_per_buffer_row ) * channels ]; + const float qua_fou = in[ ( 3 + 3 * pixels_per_buffer_row ) * channels ]; + + /* + * Computation of the YAFR correction: + * + * Basically, if two consecutive pixel value differences have the + * same sign, the smallest one (in absolute value) is taken to be + * the corresponding slope. If they don't have the same sign, the + * corresponding slope is set to 0. + * + * Four such pairs (vertical and horizontal) of slopes need to be + * computed, one pair for each of the pixels which potentially + * overlap the unit area centered at the interpolation point. + */ + /* + * Beginning of the computation of the "up" horizontal slopes: + */ + const float prem__up = dos_two - dos_one; + const float deux__up = dos_thr - dos_two; + const float troi__up = dos_fou - dos_thr; + /* + * "down" horizontal slopes: + */ + const float prem_dow = tre_two - tre_one; + const float deux_dow = tre_thr - tre_two; + const float troi_dow = tre_fou - tre_thr; + /* + * "left" vertical slopes: + */ + const float prem_left = dos_two - uno_two; + const float deux_left = tre_two - dos_two; + const float troi_left = qua_two - tre_two; + /* + * "right" vertical slopes: + */ + const float prem_rite = dos_thr - uno_thr; + const float deux_rite = tre_thr - dos_thr; + const float troi_rite = qua_thr - tre_thr; + + /* + * Back to "up": + */ + const float prem__up_squared = prem__up * prem__up; + const float deux__up_squared = deux__up * deux__up; + const float troi__up_squared = troi__up * troi__up; + /* + * Back to "down": + */ + const float prem_dow_squared = prem_dow * prem_dow; + const float deux_dow_squared = deux_dow * deux_dow; + const float troi_dow_squared = troi_dow * troi_dow; + /* + * Back to "left": + */ + const float prem_left_squared = prem_left * prem_left; + const float deux_left_squared = deux_left * deux_left; + const float troi_left_squared = troi_left * troi_left; + /* + * Back to "right": + */ + const float prem_rite_squared = prem_rite * prem_rite; + const float deux_rite_squared = deux_rite * deux_rite; + const float troi_rite_squared = troi_rite * troi_rite; + + /* + * "up": + */ + const float prem__up_times_deux__up = prem__up * deux__up; + const float deux__up_times_troi__up = deux__up * troi__up; + /* + * "down": + */ + const float prem_dow_times_deux_dow = prem_dow * deux_dow; + const float deux_dow_times_troi_dow = deux_dow * troi_dow; + /* + * "left": + */ + const float prem_left_times_deux_left = prem_left * deux_left; + const float deux_left_times_troi_left = deux_left * troi_left; + /* + * "right": + */ + const float prem_rite_times_deux_rite = prem_rite * deux_rite; + const float deux_rite_times_troi_rite = deux_rite * troi_rite; + + /* + * Branching parts of the computation of the YAFR correction (could + * be unbranched using arithmetic branching and C99 math intrinsics, + * although the compiler may be smart enough to remove the branching + * on its own): + */ + /* + * "up": + */ + const float prem__up_vs_deux__up = + prem__up_squared < deux__up_squared ? prem__up : deux__up; + const float deux__up_vs_troi__up = + deux__up_squared < troi__up_squared ? deux__up : troi__up; + /* + * "down": + */ + const float prem_dow_vs_deux_dow = + prem_dow_squared < deux_dow_squared ? prem_dow : deux_dow; + const float deux_dow_vs_troi_dow = + deux_dow_squared < troi_dow_squared ? deux_dow : troi_dow; + /* + * "left": + */ + const float prem_left_vs_deux_left = + prem_left_squared < deux_left_squared ? prem_left : deux_left; + const float deux_left_vs_troi_left = + deux_left_squared < troi_left_squared ? deux_left : troi_left; + /* + * "right": + */ + const float prem_rite_vs_deux_rite = + prem_rite_squared < deux_rite_squared ? prem_rite : deux_rite; + const float deux_rite_vs_troi_rite = + deux_rite_squared < troi_rite_squared ? deux_rite : troi_rite; + /* + * The YAFR correction computation will resume after the computation + * of the Catmull-Rom baseline. + */ + + /* + * Catmull-Rom baseline contribution: + */ + const float catmull_rom = + cardinal_uno * + ( + cardinal_one * uno_one + + + cardinal_two * uno_two + + + cardinal_thr * uno_thr + + + cardinal_fou * uno_fou + ) + + + cardinal_dos * + ( + cardinal_one * dos_one + + + cardinal_two * dos_two + + + cardinal_thr * dos_thr + + + cardinal_fou * dos_fou + ) + + + cardinal_tre * + ( + cardinal_one * tre_one + + + cardinal_two * tre_two + + + cardinal_thr * tre_thr + + + cardinal_fou * tre_fou + ) + + + cardinal_qua * + ( + cardinal_one * qua_one + + + cardinal_two * qua_two + + + cardinal_thr * qua_thr + + + cardinal_fou * qua_fou + ); + + /* + * Computation of the YAFR slopes. + */ + /* + * "up": + */ + const float mx_left__up = + prem__up_times_deux__up < 0.f ? 0.f : prem__up_vs_deux__up; + const float mx_rite__up = + deux__up_times_troi__up < 0.f ? 0.f : deux__up_vs_troi__up; + /* + * "down": + */ + const float mx_left_dow = + prem_dow_times_deux_dow < 0.f ? 0.f : prem_dow_vs_deux_dow; + const float mx_rite_dow = + deux_dow_times_troi_dow < 0.f ? 0.f : deux_dow_vs_troi_dow; + /* + * "left": + */ + const float my_left__up = + prem_left_times_deux_left < 0.f ? 0.f : prem_left_vs_deux_left; + const float my_left_dow = + deux_left_times_troi_left < 0.f ? 0.f : deux_left_vs_troi_left; + /* + * "right": + */ + const float my_rite__up = + prem_rite_times_deux_rite < 0.f ? 0.f : prem_rite_vs_deux_rite; + const float my_rite_dow = + deux_rite_times_troi_rite < 0.f ? 0.f : deux_rite_vs_troi_rite; + + /* + * Assemble the unweighted YAFR correction: + */ + const float unweighted_yafr_correction = + left_width_times_up__height_times_rite_width + * + ( mx_left__up - mx_rite__up ) + + + left_width_times_dow_height_times_rite_width + * + ( mx_left_dow - mx_rite_dow ) + + + left_width_times_up__height_times_dow_height + * + ( my_left__up - my_left_dow ) + + + rite_width_times_up__height_times_dow_height + * + ( my_rite__up - my_rite_dow ); + + /* + * Add the Catmull-Rom baseline and the weighted YAFR correction: + */ + const float newval = + sharpening_over_two * unweighted_yafr_correction + catmull_rom; + + *out = newval; +} + +static void +vips_interpolate_yafr_interpolate( VipsInterpolate *interpolate, + REGION *out, REGION *in, + int out_x, int out_y, double x, double y ) +{ + VipsInterpolateYafr *yafr = VIPS_INTERPOLATE_YAFR( interpolate ); + VipsInterpolateYafrClass *class = + VIPS_INTERPOLATE_YAFR_GET_CLASS( interpolate ); + + /* + * Note: The computation is structured to foster software + * pipelining. + */ + + /* + * x is understood to increase from left to right, y, from top to + * bottom. Consequently, ix and iy are the indices of the pixel + * located at or to the left, and at or above. the sampling point. + * + * floor is used to make sure that the transition through 0 is + * smooth. If it is known that negative x and y will never be used, + * cast (which truncates) could be used instead. + */ + const gint ix = FLOOR (x); + const gint iy = FLOOR (y); + + /* + * Each (channel's) output pixel value is obtained by combining four + * "pieces," each piece corresponding to the set of points which are + * closest to the four pixels closest to the (x,y) position, pixel + * positions which have coordinates and labels as follows: + * + * (ix,iy) (ix+1,iy) + * =left__up =rite__up + * + * <- (x,y) is somewhere in the convex hull + * + * (ix,iy+1) (ix+1,iy+1) + * =left_dow =rite_dow + */ + /* + * rite_width is the width of the overlaps of the unit averaging box + * (which is centered at the position where an interpolated value is + * desired), with the closest unit pixel areas to the right. + * + * left_width is the width of the overlaps of the unit averaging box + * (which is centered at the position where an interpolated value is + * desired), with the closest unit pixel areas to the left. + */ + const float rite_width = x - ix; + const float dow_height = y - iy; + const float left_width = 1.f - rite_width; + const float up__height = 1.f - dow_height; + /* + * .5*rite_width is the x-coordinate of the center of the overlap of + * the averaging box with the left pixel areas, relative to the + * position of the centers of the left pixels. + * + * -.5*left_width is the x-coordinate ... right pixel areas, + * relative to ... the right pixels. + * + * .5*dow_height is the y-coordinate of the center of the overlap + * of the averaging box with the up pixel areas, relative to the + * position of the centers of the up pixels. + * + * -.5*up__height is the y-coordinate ... down pixel areas, relative + * to ... the down pixels. + */ + const float left_width_times_rite_width = left_width * rite_width; + const float up__height_times_dow_height = up__height * dow_height; + + const float cardinal_two = + left_width_times_rite_width * ( -1.5f * rite_width + 1.f ) + + left_width; + const float cardinal_dos = + up__height_times_dow_height * ( -1.5f * dow_height + 1.f ) + + up__height; + + const float minus_half_left_width_times_rite_width = + -.5f * left_width_times_rite_width; + const float minus_half_up__height_times_dow_height = + -.5f * up__height_times_dow_height; + + const float left_width_times_up__height_times_rite_width = + left_width_times_rite_width * up__height; + const float left_width_times_dow_height_times_rite_width = + left_width_times_rite_width * dow_height; + const float left_width_times_up__height_times_dow_height = + up__height_times_dow_height * left_width; + const float rite_width_times_up__height_times_dow_height = + up__height_times_dow_height * rite_width; + + const float cardinal_one = + minus_half_left_width_times_rite_width * left_width; + const float cardinal_uno = + minus_half_up__height_times_dow_height * up__height; + + const float cardinal_fou = + minus_half_left_width_times_rite_width * rite_width; + const float cardinal_qua = + minus_half_up__height_times_dow_height * dow_height; + + const float cardinal_thr = + 1.f - ( minus_half_left_width_times_rite_width + cardinal_two ); + const float cardinal_tre = + 1.f - ( minus_half_up__height_times_dow_height + cardinal_dos ); + + /* + * Set the tile pointer to the first relevant value. Since the + * pointer initially points to dos_two, we need to rewind it one + * tile row, then go back one additional pixel. + */ + const PEL *p = (PEL *) IM_REGION_ADDR( in, ix - 1, iy - 1 ); + + /* Pel size and line size. + */ + const int channels = in->im->Bands; + const int pixels_per_buffer_row = + IM_REGION_LSKIP( in ) / sizeof( float ); + + /* Where we write the result. + */ + PEL *q = (PEL *) IM_REGION_ADDR( out, out_x, out_y ); + int z; + + for( z = 0; z < channels; z++ ) + catrom_yafr ((float *) q + z, (float *) p + z, + channels, pixels_per_buffer_row, + yafr->sharpening, + cardinal_one, + cardinal_two, + cardinal_thr, + cardinal_fou, + cardinal_uno, + cardinal_dos, + cardinal_tre, + cardinal_qua, + left_width_times_up__height_times_rite_width, + left_width_times_dow_height_times_rite_width, + left_width_times_up__height_times_dow_height, + rite_width_times_up__height_times_dow_height); +} + +static void +vips_interpolate_yafr_class_init( VipsInterpolateYafrClass *class ) +{ + VipsInterpolateClass *interpolate_class = + VIPS_INTERPOLATE_CLASS( class ); + + vips_interpolate_yafr_parent_class = + g_type_class_peek_parent( class ); + + interpolate_class->interpolate = vips_interpolate_yafr_interpolate; + interpolate_class->window_size = 4; +} + +static void +vips_interpolate_yafr_init( VipsInterpolateYafr *yafr ) +{ +#ifdef DEBUG + printf( "vips_interpolate_yafr_init: " ); + vips_object_print( VIPS_OBJECT( yafr ) ); +#endif /*DEBUG*/ + + yafr->sharpening = 1.0; +} + +GType +vips_interpolate_yafr_get_type( void ) +{ + static GType type = 0; + + if( !type ) { + static const GTypeInfo info = { + sizeof( VipsInterpolateYafrClass ), + NULL, /* base_init */ + NULL, /* base_finalize */ + (GClassInitFunc) vips_interpolate_yafr_class_init, + NULL, /* class_finalize */ + NULL, /* class_data */ + sizeof( VipsInterpolateYafr ), + 32, /* n_preallocs */ + (GInstanceInitFunc) vips_interpolate_yafr_init, + }; + + type = g_type_register_static( VIPS_TYPE_INTERPOLATE, + "VipsInterpolateYafr", &info, 0 ); + } + + return( type ); +} + +VipsInterpolate * +vips_interpolate_yafr_new( void ) +{ + return( VIPS_INTERPOLATE( g_object_new( + VIPS_TYPE_INTERPOLATE_YAFR, NULL ) ) ); +} + +void +vips_interpolate_yafr_set_sharpening( VipsInterpolateYafr *yafr, + double sharpening ) +{ + yafr->sharpening = sharpening; +} + +/* Convenience: return a static yafr you don't need to free. + */ +VipsInterpolate * +vips_interpolate_yafr_static( void ) +{ + static VipsInterpolate *interpolate = NULL; + + if( !interpolate ) + interpolate = vips_interpolate_yafr_new(); + + return( interpolate ); +} + + diff --git a/libsrc/mosaicing/mosaicing_dispatch.c b/libsrc/mosaicing/mosaicing_dispatch.c index 1c64f769..cdd7b944 100644 --- a/libsrc/mosaicing/mosaicing_dispatch.c +++ b/libsrc/mosaicing/mosaicing_dispatch.c @@ -565,13 +565,21 @@ affinei_vec( im_object *argv ) break; case 3: - interpolate = vips_interpolate_bilinear_new(); - vips_interpolate_bilinear_set_slow( - VIPS_INTERPOLATE_BILINEAR( interpolate ), TRUE ); + interpolate = vips_interpolate_bilinear_slow_new(); + break; + + case 4: + interpolate = vips_interpolate_yafr_new(); + break; + + case 5: + interpolate = vips_interpolate_yafr_new(); + vips_interpolate_yafr_set_sharpening( + VIPS_INTERPOLATE_YAFR( interpolate ), 2.0 ); break; default: - im_error( "affinei_vec", _( "b ad interpolation" ) ); + im_error( "affinei_vec", _( "bad interpolation" ) ); return( -1 ); }