general_filter.h

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                    */
/*                                                                      */
/*    The git repository for this software is at                        */
/*                                                                      */
/*    https://bitbucket.org/kfj/vspline                                 */
/*                                                                      */
/*    Please direct questions, bug reports, and contributions to        */
/*                                                                      */
/*    kfjahnke+vspline@gmail.com                                        */
/*                                                                      */
/*    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 general_filter.h
 
    \brief provides vspline's digital filtering capabilities without the
           b-spline-specific aspects.
 
    The two filters offered by this header are applying 'separated'
    filters: If the signal is nD, a 1D kernel is applied along one or
    all axes of the data. This requires the filter to be 'separable',
    meaning that if, for example, your data is 2D, the 2D kernel
    to be used on the data has to be a cross product of two 1D
    vectors which are applied along the two axes of the data. You
    can't use the code given here for stuff like radial basis functions
    where the kernel is not separable.
 
    vspline uses the code in filter.h for prefiltering and reconstruction.
    In both cases, the filters are set up for the specific use with b-splines.
    But the filtering code is more general, and it can be used for other
    filtering tasks. This header provides convenient access to the two
    filter types implemented in vspline:
 
    - forward-backward n-pole recursive filters
    - convolution with an arbitrary kernel
 
    This code is currently not used in other parts of vspline - vspline uses
    the code in filter.h directly.
    
    The filtering routines will appply the same filter along all axes if
    -1 is passed as 'axis'. Otherwise, the filter will only be applied to
    the specific axis. This way filtering operations with different filters
    along different axes can be achieved by calling the filtering routines
    several times, once per axis, with per-axis arguments.
*/
 
#ifndef VSPLINE_GENERAL_FILTER_H
 
#include "vspline.h"
 
// TODO: now that we have the 'generalized' routines, we might make them part of
// vspline and code 'prefilter' and 'restore' as using them, instead of directly
// calling the code in filter.h.
 
namespace vspline
{
 
/// forward_backward_recursive_filter applies one or more pairs of simple
/// recursive filters to the input. Each pair consists of:
/// 1. the forward filter: x[n] = g * x[n] + p * x[n-1]
/// 2. the backward filter: x[n] = -p * x[n] + p * x[n+1]
/// where 'p' is the 'pole' of the filter (see poles.h) and g is the
/// 'gain', which is computed as g = ( 1 - p ) * ( 1 - 1 / p );
/// multiplication the input with the gain in the first stage
/// makes sure that subsequent convolution with the corresponding
/// reconstruction kernel will restore the original signal.
/// The pair of filters used here, in sum, has no phase shift
/// (the phase shift of the forward filter is precisely matched by
/// an inverse phase shift of the backward filter due to the reversal
/// of the processing direction).
/// input and output must be equally-sized arrays of compatible types
/// (meaning, in_value_type and out_value_type must have the same number
/// of channels). bcv holds a set of boundary conditions, one per
/// dimension. pv is a std::vector holding the filter poles poles,
/// 'tolerance' is the acceptable error, defaulting to
/// a very conservative value. Try and resist from passing zero here,
/// which results in great efforts to produce a 'mathematically correct'
/// result - the result from using the conservative default will differ
/// very little from the one obtained with 'zero tolerance', but is usually
/// faster to compute.
/// Next, a 'boost' value can be passed which will be applied once
/// as a multiplicative factor. If any value apart from -1 is passed for
/// 'axis', filtering will be limited to the indicated axis. Per default,
/// the filters will be applied to all axes in turn.
/// Finally, njobs defines how many jobs will be used to multithread
/// the operation.
/// While b-spline prefiltering uses 'poles' from poles.h, this generalized
/// routine will take any poles you pass. But note that the magnitude of
/// each pole should be below 1. If you use negative poles, the effect will
/// be a high-pass filter, as it is used for b-spline prefiltering. Using
/// positive poles will blur the signal very effectively. Also note that
/// use of small negative poles will amplify the signal strongly. This
/// effect can be seen when prefiltering high-degree b-splines: the resulting
/// signal will have much higher energy. Make sure your data types can deal
/// with the increased dynamic.
 
template < unsigned int dimension ,
           typename in_value_type ,
           typename out_value_type ,
           typename math_ele_type =
                    typename vspline::common_math_ele_type
                             < in_value_type , out_value_type > ,
           size_t vsize =
                  vspline::vector_traits < math_ele_type > :: size
         >
void forward_backward_recursive_filter (
                 const
                 vigra::MultiArrayView
                   < dimension ,
                     in_value_type > & input ,
                 vigra::MultiArrayView
                   < dimension ,
                     out_value_type > & output ,
                 vigra::TinyVector
                   < bc_code ,
                     static_cast < int > ( dimension ) > bcv ,
                 std::vector < vspline::xlf_type > pv ,
                 xlf_type tolerance = -1 ,
                 xlf_type boost = xlf_type ( 1 ) ,
                 int axis = -1 , // -1: apply along all axes
                 int njobs = default_njobs )
{
  if ( output.shape() != input.shape() )
    throw shape_mismatch
     ( "forward_backward_recursive_filter: input and output shape must match" ) ;
 
  if ( tolerance < 0 )
    tolerance = std::numeric_limits < math_ele_type > :: epsilon() ;
 
  int npoles = pv.size() ;
  
  if ( npoles < 1 )
  {
    // if npoles < 1, there is no filter to apply, but we may need
    // to apply 'boost' and/or copy input to output. We use 'amplify'
    // for the purpose, which multithreads the operation (if it is at
    // all necessary). I found this is (slightly) faster than doing the
    // job in a single thread - the process is mainly memory-bound, so
    // the gain is moderate.
 
    amplify < dimension , in_value_type , out_value_type , math_ele_type >
      ( input , output , math_ele_type ( boost ) ) ;
      
    return ;
  }
  
 
  vspline::xlf_type poles [ npoles ] ;
  for ( int i = 0 ; i < npoles ; i++ )
    poles[i] = pv[i] ;
  
  // KFJ 2018-05-08 with the automatic use of vectorization the
  // distinction whether math_ele_type is 'vectorizable' or not
  // is no longer needed: simdized_type will be a Vc::SimdArray
  // if possible, a vspline::simd_type otherwise.
  
  typedef typename vspline::bspl_prefilter
                            < vspline::simdized_type ,
                              math_ele_type ,
                              vsize
                            > filter_type ;
 
  // now call the 'wielding' code in filter.h
 
  if ( axis == -1 )
  {
    // user has passed -1 for 'axis', apply the same filter along all axes
 
    std::vector < vspline::iir_filter_specs > vspecs ;
 
    // package the arguments to the filter; one set of arguments
    // per axis of the data
 
    for ( int axis = 0 ; axis < dimension ; axis++ )
    {
      vspecs.push_back
        ( vspline::iir_filter_specs
          ( bcv [ axis ] , npoles , poles , tolerance , 1 ) ) ;
    }
 
    // 'boost' is only applied to dimension 0, since it is meant to
    // affect the whole data set just once, not once per axis.
 
    vspecs [ 0 ] . boost = boost ;
    
    vspline::filter
    < in_value_type , out_value_type , dimension , filter_type > 
    ( input , output , vspecs , njobs ) ;
  }
  else
  {
    // user has passed a specific axis, apply filter only to this axis
 
    assert ( axis >=0 && axis < dimension ) ;
 
    vspline::filter
    < in_value_type , out_value_type , dimension , filter_type > 
    ( input , output , axis ,
      vspline::iir_filter_specs (
        bcv [ axis ] , npoles , poles , tolerance , boost ) ,
      njobs ) ;
  }
}
 
/// convolution_filter implements convolution of the input with a fixed-size
/// convolution kernel. Note that vspline does *not* follow the DSP convention
/// of using the kernel's coefficients in reverse order. The standard is to
/// calculate sum ( ck * u(n-k) ), vspline uses sum ( ck * u(n+k-h) ) where
/// 'h' is the 'headroom' of the kernel - the number of coefficients which
/// are applied 'to 'past' values, so that for for a kernel with three
/// coefficients and headroom 1, the sum at position n would be
/// c0 * u(n-1) + c1 * u(n) + c2 * u(n+1), and for a kernel with four
/// coefficients and headroom 2, you'd get
/// c0 * u(n-2) + c1 * u(n-1) + c2 * u(n) + c3 * u(n+1)
/// If you use a 'normal' kernel in vspline, you must reverse it (unless it's
/// symmetric, of course), and you must state the 'headroom': if your kernel
/// is odd-sized, this will normally be half the kernel size (integer division!),
/// producing the phase-correct result, and with an even kernel you can't get
/// the phase right, and you have to make up your mind which way the phase
/// should shift. The 'headroom' notation leaves you free to pick your choice.
/// The input, output, bcv, axis and njobs parameters are as in the routine
/// above. Here what's needed here is 'kv', a std::vector holding the filter's
/// coefficients, and the 'headroom', usually kv.size() / 2.
///
/// example: let's say you have data in 'image' and want to convolve in-place
/// with [ 1 , 2 , 1 ], using mirror boundary conditions. Then call:
///
///   vspline::convolution_filter
///    ( image ,
///      image ,
///      { vspline::MIRROR , vspline::MIRROR } ,
///      { 1 , 2 , 1 } ,
///      1 ) ;
 
 
template < unsigned int dimension ,
           typename in_value_type ,
           typename out_value_type ,
           typename math_ele_type =
                    typename vspline::common_math_ele_type
                             < in_value_type , out_value_type > ,
           size_t vsize =
                  vspline::vector_traits < math_ele_type > :: size
         >
void convolution_filter (
                 const
                 vigra::MultiArrayView
                   < dimension ,
                     in_value_type > & input ,
                 vigra::MultiArrayView
                   < dimension ,
                     out_value_type > & output ,
                 vigra::TinyVector
                   < bc_code ,
                     static_cast < int > ( dimension ) > bcv ,
                 std::vector < vspline::xlf_type > kv ,
                 int headroom ,
                 int axis = -1 , // -1: apply along all axes
                 int njobs = default_njobs )
{
  if ( output.shape() != input.shape() )
    throw shape_mismatch
     ( "convolution_filter: input and output shape must match" ) ;
 
  size_t ksize = kv.size() ;
  
  if ( ksize < 1 )
  {
    // we can handle the no-kernel case very efficiently,
    // since we needn't apply a filter at all.
 
    if ( (void*) ( input.data() ) != (void*) ( output.data() ) )
    {
      // operation is not in-place, copy data to output
      output = input ;
    }
    return ;
  }
  
  vspline::xlf_type kernel [ ksize ] ;
  for ( int i = 0 ; i < ksize ; i++ )
    kernel[i] = kv[i] ;
  
  typedef typename vspline::convolve
                            < vspline::simdized_type ,
                              math_ele_type ,
                              vsize
                            > filter_type ;
 
  if ( axis == -1 )
  {
    // user has passed -1 for 'axis', apply the same filter along all axes
 
    std::vector < vspline::fir_filter_specs > vspecs ;
  
    for ( int axis = 0 ; axis < dimension ; axis++ )
    {
      vspecs.push_back 
        ( vspline::fir_filter_specs
          ( bcv [ axis ] , ksize , headroom , kernel ) ) ;
    }
 
    vspline::filter
    < in_value_type , out_value_type , dimension , filter_type > 
    ( input , output , vspecs , njobs ) ;
  }
  else
  {
    // user has passed a specific axis, apply filter only to this axis
 
    assert ( axis >=0 && axis < dimension ) ;
 
    vspline::filter
    < in_value_type , out_value_type , dimension , filter_type > 
    ( input , output , axis ,
      vspline::fir_filter_specs ( bcv [ axis ] , ksize , headroom , kernel ) ,
      njobs ) ;
  }
}
 
} ; // namespace vspline
 
#define VSPLINE_GENERAL_FILTER_H
#endif