metafilter3.cc

Go to the documentation of this file.

/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2015 - 2023 by Kay F. Jahnke                    */
/*                                                                      */
/*    Permission is hereby granted, free of charge, to any person       */
/*    obtaining a copy of this software and associated documentation    */
/*    files (the "Software"), to deal in the Software without           */
/*    restriction, including without limitation the rights to use,      */
/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
/*    sell copies of the Software, and to permit persons to whom the    */
/*    Software is furnished to do so, subject to the following          */
/*    conditions:                                                       */
/*                                                                      */
/*    The above copyright notice and this permission notice shall be    */
/*    included in all copies or substantial portions of the             */
/*    Software.                                                         */
/*                                                                      */
/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
/*                                                                      */
/************************************************************************/
 
/// \file metafilter3.cc
///
/// \brief implementing a locally adapted filter
///
/// taking the method introduced in metafilter.cc one step further,
/// this file uses a filter which varies with the locus of it's
/// application. Since the loci of the pickup points, relative to the
/// filter's center, are decoupled from the filter weights, we can
/// 'reshape' the filter by manipulating them. This program adapts the
/// filter so that it is applied to a perspective-corrected view of
/// the image - or, to express it differently (and this is what happens
/// technically) - the filter is reshaped with the locus of it's application.
/// So the filter is applied to 'what was seen' rather than
/// 'what the image holds'. This is a subtle difference, which is not
/// easy to spot - you can provoke a clearly visible result by specifying
/// a large hfov, and also by applying the program to images with sharp
/// contrasts and single off-coloured pixels; with ordinary photographs
/// and 'normal' viewing angles you'll notice a blur increasing towards
/// the edges of the image. In a wide-angle shot, small circular
/// objects near the edges (especially the corners) of the image appear
/// like ellipses. Filtering the image with a normal convolution with
/// a symmetric filter (like the binomial we're using here) will result
/// in the ellipses being surrounded by a blurred halo which has the same
/// width everywhere. Using the filter implemented here, we see an image
/// of a blurred circular object: the blurred halo is wider in radial
/// direction. And this is precisely what 'should' be seen: Filtering
/// the image with a static filter is a simplification which produces
/// results that look reasonably convincing given a reasonable field 
/// of view, but to 'properly' model some change in viewing conditions
/// (like haze in the atmosphere) we need to model what was seen.
/// While this distinction may look like hair-splitting when it comes to
/// photography, it may be more convincing when you consider feature
/// detection. If you have a filter detecting circular patterns, the
/// filter's response to the elliptical shapes occurring in a wide-angle
/// shot near the margins will be sub-optimal: the feature detector will
/// not respond maximally, because, after all, it's input is an ellipse
/// and not a circle. But if you use the technique I implement in this
/// program, the detector will adapt to the place it 'looks at' and detect
/// representations of circular shapes in the image with proper full
/// response. In fact, such a detector will react with lessened response if
/// the image, near the border, shows circular structures, because these
/// clearly can't have originated from circular objects in the view.
/// This program works with rectilinear input, but the implications are
/// even greater for, say, full spherical panoramas. In such images, there
/// is intense distortion 'near the poles', and using ordinary filters on
/// such images does not produce truly correct effects. With filters which
/// adapt to the locus of their application, such images can be filtered
/// adequately.
/// A word of warning is in order here: you still have to obey the
/// sampling theorem and make sure that the pickup points ar not further
/// than half the minimal wave length captured in the image (expressed as
/// an angle in spherical coordinates). Otherwise you will get aliasing,
/// which may produce visible artifacts. If in doubt, use a larger kernel
/// and scale it down. 'Kernel' size can be adapted freely with the technique
/// proposed here. To see the filter 'at work' process an image with a few
/// isolated black dots on white background and inspect the output, which
/// shows a different result, depending on the distance from the image's
/// center, rather than a uniform impulse response everywhere.
/// All of this sounds like a lot of CPU cycles to be gotten through,
/// but in fact due to vspline's use of SIMD and multi-threading it only
/// takes a fraction of a second for a full HD image.
/// In this variant of the program, we add a few 'nice-to-have' features
/// which are easy to implement with the given tool set: we add the option
/// to scale the filter - we can do this steplessly because the positions
/// of the pickup points are decoupled from the grid. We also use a larger
/// kernel (9X9) which we can use over a wider range of scales - if we
/// scale it so that it 'covers' just about a single pixel in the center,
/// it will be larger towards the image margins, the increase depending on
/// the FOV. If we use it on scaled versions of the image, we can scale it
/// so that it covers an equally sized patch (e.g. x arch seconds in
/// diameter) on every scale. We could be tempted to use it unscaled on
/// the image with the most detail, but we have to keep in mind that in
/// the corners of the image, this will put the pickup points further
/// apart than one unit step, which will only produce correct results
/// if excessively high frequencies are absent from the spectrum in these
/// parts of the image.
/// This suggests that using metafilters can be useful for work with
/// image pyramids: given two adjacent layers of the pyramid L and M,
/// where M is a scaled-down version of L with a scaling factor of
/// l (in 'conventional' pyramids, l = 1/2), we can choose a metakernel
/// so that it will not be too 'scattered' when projected to the edges
/// of L (so as to avoid aliasing) and scale the metakernel with l
/// when we apply it to M. Then we expect similar filter response from
/// both cases at corresponding loci - we have killed two birds with
/// one stone: the filter response is independent of scale and reflects
/// 'what is looked at' rather than it's projected image. With the free
/// scalability of the metafiler, we can use other l than 0.5 which
/// allows 'steeper' and 'shallower' pyramids (ref. to spline pyramids)
/// with the equivalence of applying two differently-scaled versions of
/// a metafilter to differently-scaled versions of an image, we can
/// also 'extrapolate' filters: rather than applying a metafilter with
/// a large number of coefficients to a detailed image, we can apply
/// a smaller metafilter to a scaled-down version of the image. (note
/// that this puts constraints on the filter, it should be akin to a
/// gaussian).
 
// We use the following strategy:
// codify the image's metrics in an object which provides methods to
// convert an image coordinate into spherical coordinates and ray
// coordinates into image coordinates.
// Then convert each batch of incoming 2D image coordinates to spherical
// coordinates. The spherical coordinates are interpreted as a rotation,
// and this rotation is applied to each of the kernel coefficients' loci
// in turn, yielding 'pickup coordinates' where the image is evaluated.
// The evaluations for all coefficients are weighted and summed up,
// yielding the filter response, which is stored in the target image.
// So the mathematics are quite involved: For each image point, we
// have to calculate the rotational quaternion, apply it to all the
// kernel coefficients, evaluate at these locations, and form the weighted
// sum. The SIMD code used here does that with batches of 16 points at a
// time, which is a great speed-up.
// Note how this strategy is SIMD-friendly, whereas following the concept
// of 'reorient a ray bundle to each source coordinate' is not so easily
// calculated with SIMD arithmetic. The outcome of both schemes is the
// same, but changing the sequence around allows us to go 'in SIMD'
// from a SIMD vector of target coordinates to a SIMD vector of pickup
// coordinates to a series of SIMD vectors of weighted partial results
// corresponding to kernel coefficients, and finally, when this series
// is summed up, to a SIMD vector of result pixels.
 
#include <iostream>
 
#include <vspline/vspline.h>
 
// we want to use some optional headers from vspline, which make our
// lives easier, allowing more concise notation.
 
#define XEL_TYPE vigra::TinyVector
 
#define VECTOR_TYPE vspline::simd_type
#include "vspline/opt/vector_promote.h"
#include "vspline/opt/xel_of_vector.h"
#undef VECTOR_TYPE
 
#ifdef USE_VC
 
#define VECTOR_TYPE vspline::vc_simd_type
#include "vspline/opt/vector_promote.h"
#include "vspline/opt/xel_of_vector.h"
#undef VECTOR_TYPE
 
#endif
 
#ifdef USE_HWY
 
#define VECTOR_TYPE vspline::hwy_simd_type
#include "vspline/opt/vector_promote.h"
#include "vspline/opt/xel_of_vector.h"
#undef VECTOR_TYPE
 
#endif
 
#undef XEL_TYPE
 
#include <vigra/stdimage.hxx>
#include <vigra/imageinfo.hxx>
#include <vigra/impex.hxx>
#include <vigra/quaternion.hxx>
 
// we silently assume we have a colour image
 
typedef vigra::RGBValue < float , 0 , 1 , 2 > pixel_type ; 
 
// coordinate_type has a 2D coordinate
 
typedef vigra::TinyVector < float , 2 > coordinate_type ;
 
// ray_type has a 3D coordinate
 
typedef vigra::TinyVector < float , 3 > ray_type ;
 
// type of b-spline object
 
typedef vspline::bspline < pixel_type , 2 > spline_type ;
 
// target_type is a 2D array of pixels 
 
typedef vigra::MultiArray < 2 , pixel_type > target_type ;
 
// filter coefficients are TinyVectors of three values:
// the horizontal and vertical distance from the filter's
// center, and the corresponding weight.
 
typedef vigra::TinyVector < float , 3 > coefficient_type ;
 
// kernel_type is a 2D MultiArray of coefficients.
 
typedef vigra::MultiArray < 2 , coefficient_type > kernel_type ;
 
// image_base_type has the basic metrics of an image both in pixel
// units and model space units.
 
struct image_base_type
{
  // metrics of the image in pixel units
 
  size_t px_width ;
  size_t px_height ;
  double hfov ;
 
  // corresponding width and height values in unit sphere radii
  // to clarify: the rectilinear source image is 'mounted' in model
  // space at one unit distance from model center, and it's horizontal
  // extent in model space, u_width, is equal to the twice the tangens
  // of half the horizontal field of view. u_height is cacluclated
  // from u_width via the aspect ratio. px_per_u and u_per_px are
  // simply linear scaling factors.
 
  double u_width ;
  double u_height ;
 
  // ratio of pixel units to unit sphere radii (u) and inverse
 
  double px_per_u ;
  double u_per_px ;
 
  image_base_type ( const size_t & _px_width ,
                    const size_t & _px_height ,
                    const double & _hfov )
  : px_width ( _px_width ) ,
    px_height ( _px_height ) ,
    hfov ( _hfov )
  {
    assert ( hfov >= 0.0 ) ;
    u_width = 2.0 * tan ( hfov / 2.0 ) ;
    u_height = ( u_width * px_height ) / px_width ;
    px_per_u = px_width / u_width ;
    u_per_px = u_width / px_width ;
  }
} ;
 
// image_type provides conversion from 2D rectilinear image coordinates
// to spherical coordinates, and from 3D cartesian to 2D rectilinear.
// It uses two vspline::domains to encode the shift and scale operation,
// and the reprojection.
 
template < std::size_t vsize = vspline::vector_traits < float > :: size >
struct image_type
{
  const vspline::domain_type < coordinate_type , vsize > scale_to_model ;
  const vspline::domain_type < coordinate_type , vsize > scale_to_image ;
  const image_base_type image_base ;
 
private:
 
  // scale_to_model transforms image coordinates to model space coordinates.
  // This is a simple scale-and-shift operation which is implemented with
  // a vspline::domain_type. scale_to_image performs the inverse operation.
 
  image_type ( const image_base_type & _ib )
  : image_base ( _ib ) ,
    scale_to_model ( { -.5f ,
                       -.5f } ,
                     { float ( _ib.px_width - .5 ) ,
                       float ( _ib.px_height - 0.5 ) } ,
                     { float ( - _ib.u_width / 2.0 ) ,
                       float ( - _ib.u_height / 2.0 ) } ,
                     { float ( _ib.u_width / 2.0 ) ,
                       float ( _ib.u_height / 2.0 ) } ) ,
    scale_to_image ( { float ( - _ib.u_width / 2.0 ) ,
                       float ( - _ib.u_height / 2.0 ) } ,
                     { float ( _ib.u_width / 2.0 ) ,
                       float ( _ib.u_height / 2.0 ) } ,
                     { -.5f ,
                       -.5f } ,
                     { float ( _ib.px_width - .5 ) ,
                       float ( _ib.px_height - 0.5 ) } )
  { }
 
public:
 
  // public c'tor taking image width, height and hfov
 
  image_type ( const size_t & _px_width ,
               const size_t & _px_height ,
               const double & _hfov )
  : image_type ( image_base_type ( _px_width ,
                                   _px_height ,
                                   _hfov ) )
  { }
 
  // to_spherical calculates azimuth and elevation from a 2D image coordinate
 
  template < typename dtype >
  vigra::TinyVector < dtype , 2 > to_spherical
    ( const vigra::TinyVector < dtype , 2 > & ci ) const
  {
    // cm will hold the coordinate in model space
 
    vigra::TinyVector < dtype , 2 > cm ;
 
    // transform the incoming image coordinate to model space
 
    scale_to_model.eval ( ci , cm ) ;
 
    // find azimuth and elevation
 
    dtype azimuth = atan ( cm[0] ) ;
 
    // distance from origin to point projected to horozontal plane
 
    dtype fl = sqrt ( 1.0 + cm[0] * cm[0] ) ;
 
    dtype elevation = - atan ( cm[1] / fl ) ;
 
    return { azimuth , elevation } ;
  }
 
  // to_image intersects a ray with the image plane, yielding 2D rectilinear.
 
  template < typename dtype >
  vigra::TinyVector < dtype , 2 > to_image
    ( const vigra::TinyVector < dtype , 3 > & cm3 ) const
  {
    vigra::TinyVector < dtype , 2 > cm { cm3[1] / cm3[0] , cm3[2] / cm3[0] } ;
 
    scale_to_image.eval ( cm , cm ) ;
 
    return cm ;
  }
 
} ;
 
// adapted quaternion code taken from lux, omitting the roll component rotation
// which is not used in this program. Note how the inclusion of the headers
// from vspline's opt section introduced traits code for vigra which allows
// us to create vigra::Quaternions with simdized dtype. With this addition
// we can formulate our pixel pipeline in a type-agnostic way so that it
// functions with both scalar and SIMD data.
 
template < typename dtype >
vigra::Quaternion < dtype > get_rotation_q ( dtype yaw , dtype pitch )
{
  // first calculate the component quaternions.
  // for a right-handed system and clockwise rotations, we have this formula:
  // q = cos ( phi / 2 ) + ( ux * i + uy * j + uz * k ) * sin ( phi / 2 )
  // where phi is the angle of rotation  and (ux, uy, uz) is the axis around
  // which we want to rotate some 3D object.
 
  // if we consider a target in the line of sight, we take pitch as moving it upwards,
  // rotating clockwise around the positive y axis (0 , 1, 0)
  vigra::Quaternion < dtype > q_pitch ( cos ( pitch / 2.0 ) , 0.0 , sin ( pitch / 2.0 ) , 0.0 ) ;
 
  // if we consider a target in the line of sight, we take yaw as moving it to the right,
  // which is counterclockwise around the z axis, or clockwise around the negative
  // z axis (0, 0, -1)
  vigra::Quaternion < dtype > q_yaw ( cos ( -yaw / 2.0 ) , 0.0 , 0.0 , - sin ( -yaw / 2.0 ) ) ;
 
  // now produce the concatenated operation by multiplication of the components
  vigra::Quaternion < dtype > total = q_yaw ;
  total *= q_pitch ;
 
  return total ;
}
 
/// apply the rotation codified in a quaternion to a 3D point
 
template < typename dtype , typename qtype >
vigra::TinyVector < dtype , 3 >
rotate_q ( const vigra::TinyVector < dtype , 3 > & v ,
           vigra::Quaternion < qtype > q )
{
  vigra::TinyVector < dtype , 3 > result ;
  vigra::Quaternion < dtype > vq ( 0.0 , v[0] , v[1] , v[2] ) ;
  vigra::Quaternion < dtype > qq ( q[0] , q[1] , q[2] , q[3] ) ;
  vigra::Quaternion < dtype > qc = conj ( qq ) ;
  qq *= vq ;
  qq *= qc ;
  result[0] = qq[1] ;
  result[1] = qq[2] ;
  result[2] = qq[3] ;
  return result ;
}
 
// the filter is built by functor composition. EV is the
// type of the interpolator producing the pickup values.
// The kernel's coefficients are expected to be scaled to
// their final magnitude, and usually the weight components
// will have been normalized.
// The difference to a normal kernel is that the coefficients'
// loci are not necessarily on a grid, and if they are on a grid,
// they needn't be unit-spaced. If you think of the kernel as a
// set of weighted points on a plane tangent to the unit sphere
// in forward direction, each coefficient's locus defines a ray.
// The filter's output is generated by moving the bundle of rays
// so that their center coincides with the incoming coordinate,
// and the interpolator is evaluated where each of the rays
// hit the image surface. The evaluations are weighted and summed
// up to yield the result.
// With this design, we gain a lot of flexibility: we can scale the
// kernel up and down and use unconventional coefficient patterns.
// And the filter's response is of 'what was looked at' rather than
// of 'what is in the image'. It does in effect have inherent
// perspective correction. So it's a good base for feature detectors,
// which should not be fed 'warped' data for optimal response.
// Note that the funny-sounding comments which mix singular and plural
// try and reflect the fact that the code can process single coordinates
// as well as 'simdized' coordinates.
 
template < typename I , typename O , size_t S ,
           template < typename , typename , size_t > class EV >
struct meta_filter
: public vspline::unary_functor < I , O , S >
{
  typedef vspline::unary_functor < I , O , S > base_type ;
  typedef EV < I , O , S > inner_type ;
 
  const inner_type & inner ;
  const kernel_type & kernel ;
  const image_type < S > img ;
 
  meta_filter ( const inner_type & _inner ,
                const kernel_type & _kernel ,
                const coordinate_type & _extent ,
                const float & _hfov )
  : inner ( _inner ) ,
    kernel ( _kernel ) ,
    img ( _extent[0] , _extent[1] , _hfov )
  { } ;
  
  // since the code is the same for vectorized and unvectorized
  // operation, we can write an 'eval' template which can be used
  // both for single-value evaluation and vectorized evaluation.
  // incoming coordinate 'c' is a 2D image coordinate or it's
  // SIMD equivalent, a vigra::TinyVector of two SIMD vectors.
  
  template < class IN , class OUT >
  void eval ( const IN & c ,
                    OUT & result ) const
  {
    // elementary type of IN, may be a scalar or a vector
 
    typedef typename IN::value_type ele_type ;
 
    // clear 'result'
 
    result = 0.0f ;
 
    // convert the incoming coordinate(s) to spherical
 
    auto cc = img.to_spherical ( c ) ;
 
    // - calculate a quaternion containing two rotations: the azimuth
    //   as pitch and the elevation as yaw
 
    auto q = get_rotation_q ( cc[0] , cc[1] ) ;
 
    // iterate over the coefficients
 
    for ( auto const & cf : kernel )
    {
      // cf now holds three values: cf[0] and cf[1] hold the distance
      // of the kernel coefficient from the kernel's center - these two
      // values would be implicit in a 'normal' kernel and result from
      // the coefficient's place in the kernel, but here we allow
      // arbitrary values for maximal flexibility. cf[2] holds the
      // kernel coefficient's value (the 'weight'), which is the only
      // value a normal kernel would hold for each coefficient.
      // Note that the values in cf are not SIMD vectors but simple
      // scalars. We use references for clarity:
 
      const auto & dx ( cf[0] ) ;
      const auto & dy ( cf[1] ) ;
      const auto & weight ( cf[2] ) ;
 
      // cc has spherical coordinates of the current target pixel(s).
      // we have to perform these steps:
 
      // - apply quaternion q to (1, dx, dy) , which yields the pickup
      //   coordinate(s). We need a 3D point, so initially we position it
      //   one unit 'in front':                VVVV
 
      vigra::TinyVector < ele_type , 3 > p3d { 1.0f , dx , dy } ;
 
      // now we apply the quaternion, receiving (a) 3D 'ray' coordinate(s)
 
      p3d = rotate_q ( p3d , q ) ;
 
      // now convert p3d back to image coordinates to yield the
      // pickup coordinate(s) in 2D image coordinates.
 
      auto ci = img.to_image ( p3d ) ;
 
      // there, evaluate the interpolator to yield the pickup value(s)
 
      OUT pickup ; // to hold partial result(s) from pickup point(s)
 
      inner.eval ( ci , pickup ) ;
 
      // apply the current weight to the pickup value(s) and sum up,
      // yielding the filter output(s)
 
      result += pickup * weight ;
    }
  } 
  
} ;
 
// for starters, we provide building a kernel from an ordinary 2D kernel
// with equidistant coefficients, as it would be used for a convolution.
// The loci of the coefficients are calculated by referring to 'img';
// they are given in model space units. Here, their distance is picked
// so that they coincide with the distance of two pixels at the image's
// center. Note how this kernel is for demonstration purposes and
// deliberately simple. It does not really exploit the flexibility which
// the free positioning of the coefficients relative to the center would
// allow.
 
void build_kernel ( const vigra::MultiArrayView < 2 , float > & kernel_2d ,
                    const image_base_type & img ,
                    kernel_type & meta_kernel )
{
  // we get the 2D kernel's extent
 
  auto kernel_shape = kernel_2d.shape() ;
 
  coordinate_type kernel_center ( kernel_shape[0] - 1 , kernel_shape[1] - 1 ) ;
  kernel_center /= 2 ;
 
  assert ( kernel_shape == meta_kernel.shape() ) ;
 
  for ( int x = 0 ; x < kernel_shape[0] ; x++ )
  {
    for ( int y = 0 ; y < kernel_shape[1] ; y++ )
    {
      float dx = img.u_per_px * ( x - kernel_center [ 0 ] ) ;
      float dy = img.u_per_px * ( y - kernel_center [ 1 ] ) ;
      float weight = kernel_2d [ { x , y } ] ;
      
      meta_kernel [ { x , y } ] = { dx , dy , weight } ; 
    }
  }
  // meta_kernel [ { 0 , 0 } ] [ 2 ] = 1.0 ;
}
 
// to scale the kernel, we simply multiply the 'locus' parts of
// the coefficients with a factor. Note how we can exploit vigra's
// convenient capability of multiplying an entire array of T
// (namely k) with one T (cf in this case):
 
void scale_kernel ( kernel_type & k , double sx , double sy )
{
  coefficient_type cf { float ( sx ) , float ( sy ) , 1.0f } ;
  k *= cf ;
}
 
int main ( int argc , char * argv[] )
{
  if ( argc < 5 )
  {
    std::cerr << "pass a colour image file as first argument" << std::endl ;
    std::cerr << "and the hfov as second argument, followed by" << std::endl ;
    std::cerr << "the horizontal and vertical kernel scaling factor" << std::endl ;
    exit( -1 ) ;
  }
 
  float hfov = std::atof ( argv[2] ) * M_PI / 180.0 ; // 90 degree in radians
  float sx = std::atof ( argv[3] ) ; // horizontal kernel scaling factor
  float sy = std::atof ( argv[4] ) ; // vertical kernel scaling factor
 
  // get the image file name
 
  vigra::ImageImportInfo imageInfo ( argv[1] ) ;
 
  image_base_type img ( imageInfo.width() , imageInfo.height() , hfov ) ;
 
  // create cubic 2D b-spline object containing the image data
  // using a cubic b-spline here is quite expensive - after all, for every
  // pickup point this spline has to be evaluated. If the kernel is large,
  // this is overkill, and even a degree-1 spline (bilinear interpolation)
  // yields 'good enough' results.
 
  spline_type bspl ( imageInfo.shape() ) ;
 
  // load the image data into the b-spline's core.
 
  vigra::importImage ( imageInfo , bspl.core ) ;
  
  // prefilter the b-spline
 
  bspl.prefilter() ;
 
  // create a 'safe' evaluator from the b-spline, so we needn't worry
  // about out-of-range access; the default safe evaluator uses mirror
  // boundary conditions and maps out-of-range coordinates into the range
  // by mirroring on the bounds
 
  auto ev = vspline::make_safe_evaluator ( bspl ) ;
 
  // now we set up the meta filter. Here we use a sampling of a quadratic
  // b-spline kernel with a sub-unit-spaced 9*9 grid.
 
  vspline::basis_functor < double > bf ( 2 ) ;
 
  vigra::MultiArray < 2 , float > kernel_2d ( vigra::Shape2 ( 9 , 9 ) ) ;
  double sum = 0.0 ;
 
  for ( int x = -4 ; x < 5 ; x++ )
  {
    for ( int y = -4 ; y < 5 ; y++ )
    {
      double kx = bf ( 0.25 * x ) ;
      double ky = bf ( 0.25 * y ) ;
      double k = kx * ky ;
      std::cout << " " << k ;
      kernel_2d [ { x + 4 , y + 4 } ] = k ;
      sum += k ;
    }
    std::cout << std::endl ;
  }
 
  // normalize to 1.0
 
  kernel_2d /= sum ;
 
  kernel_type meta_kernel ( vigra::Shape2 ( 9 , 9 ) ) ;
 
  build_kernel ( kernel_2d , img , meta_kernel ) ;
  scale_kernel ( meta_kernel , sx , sy ) ;
 
  // now we create the actual meta filter function from the
  // interpolator and the filter data. Note that a large hfov is
  // needed to actually see that the result of the filter looks
  // different in the center and towards the edges and especially
  // the corners.
 
  auto mev = meta_filter < coordinate_type ,
                           pixel_type ,
                           vspline::vector_traits < float > :: size ,
                           vspline::grok_type >
               ( ev , meta_kernel , imageInfo.shape() , hfov ) ;
 
  // this is where the result should go:
 
  target_type target ( imageInfo.shape() ) ;
 
  // finally, we trigger the pixel pipeline by passing the functor 'mev'
  // and the to-be-filled target array to vspline::transform, which
  // takes care of 'wielding' or 'rolling out' the functor, filling
  // the target array with result data.
 
  vspline::transform ( mev , target ) ;
 
  // store the result with vigra impex
 
  vigra::ImageExportInfo eximageInfo ( "metafilter3.tif" );
  
  std::cout << "storing the target image as 'metafilter3.tif'" << std::endl ;
  
  vigra::exportImage ( target ,
                       eximageInfo
                       .setPixelType ( "UINT16" )
                       .setForcedRangeMapping ( 0 , 255 , 0 , 65535 ) ) ;
  
  exit ( 0 ) ;
}