gsm2.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 gsm2.cc
///
/// \brief  calculating the gradient squared magnitude, derivatives
///
/// alternative implementation of gsm.cc, performing the calculation of the
/// gradient squared magnitude with a functor and transform(), which is faster
/// since the whole operation is multithreaded and potentially vectorized.
///
/// compile with:
/// clang++ -std=c++11 -march=native -o gsm -O3 -pthread -DUSE_VC gsm.cc -lvigraimpex -lVc
/// or clang++ -std=c++11 -march=native -o gsm -O3 -pthread gsm.cc -lvigraimpex
/// (with Vc; use -DUSE_HWY and -lhwy for highway, or -std=c++17 and -DUSE_STDSIMD
/// for the std::simd backend)
///
/// invoke passing an image file. the result will be written to 'gsm2.tif'
 
#include <iostream>
 
#include <vspline/vspline.h>
 
#include <vigra/stdimage.hxx>
#include <vigra/imageinfo.hxx>
#include <vigra/impex.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 ;
 
// 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 ;
 
// b-spline evaluator producing float pixels
typedef vspline::evaluator < coordinate_type , // incoming coordinate's type
                             pixel_type        // singular result data type
                           > ev_type ;
 
/// we build a vspline::unary_functor which calculates the sum of gradient squared
/// magnitudes.
/// Note how the 'compound evaluator' we construct follows a pattern of
/// - derive from vspline::unary_functor
/// - keep const instances of 'inner' types
/// - initialize these in the constructor, yielding a 'pure' functor
/// - if the vector code is identical to the unvectorized code, implement
///   eval() with a template                      
 
// we start out by coding the evaluation functor.
// this class does the actual computations:
 
struct ev_gsm
: public vspline::unary_functor < coordinate_type , pixel_type >
{
  // we create two evaluators for the b-spline, one for the horizontal and
  // one for the vertical gradient. The derivatives for a b-spline are requested
  // by passing a TinyVector with as many elements as the spline's dimension
  // with the desired derivative degree for each dimension. Here we want the
  // first derivatives in x and in y direction:
 
  const vigra::TinyVector < float , 2 > dx1_spec { 1 , 0 } ;
  const vigra::TinyVector < float , 2 > dy1_spec { 0 , 1 } ;
 
  // we keep two 'inner' evaluators, one for each direction
  
  const ev_type xev , yev ;
  
  // which are initialized in the constructor, using the bspline and the
  // derivative specifiers
  
  ev_gsm ( const spline_type & bspl )
  : xev ( bspl , dx1_spec ) ,
    yev ( bspl , dy1_spec )
  { } ;
  
  // since the code is the same for vectorized and unvectorized
  // operation, we can write a template:
  
  template < class IN , class OUT >
  void eval ( const IN & c ,
                    OUT & result ) const
  {
    OUT dx , dy ;
    
    xev.eval ( c , dx ) ; // get the gradient in x direction
    yev.eval ( c , dy ) ; // get the gradient in y direction
    
    // really, we'd like to write:
    // result = dx * dx + dy * dy ;
    // but fail due to problems with type inference: for vectorized
    // types, both dx and dy are vigra::TinyVectors of two simdized
    // values, and multiplying two such objects is not defined.
    // It's possible to activate code performing such operations by
    // defining the relevant traits in namespace vigra, but for this
    // example we'll stick with using multiply-assigns, which are defined.
    
    dx *= dx ;            // square the gradients
    dy *= dy ;
    dx += dy ;            // add them up
    
    result = dx ;         // assign to result
  } 
  
} ;
 
int main ( int argc , char * argv[] )
{
  if ( argc < 2 )
  {
    std::cerr << "pass a colour image file as argument" << std::endl ;
    exit( -1 ) ;
  }
  // get the image file name
  
  vigra::ImageImportInfo imageInfo ( argv[1] ) ;
 
  // we want a b-spline with natural boundary conditions
  
  vigra::TinyVector < vspline::bc_code , 2 > bcv ( vspline::NATURAL ) ;
  
  // create cubic 2D b-spline object containing the image data
  
  spline_type bspl ( imageInfo.shape() , // the shape of the data for the spline
                     3 ,                 // degree 3 == cubic spline
                     bcv                 // specifies natural BCs along both axes
                   ) ;
  
  // load the image data into the b-spline's core. This is a common idiom:
  // the spline's 'core' is a MultiArrayView to that part of the spline's
  // data container which precisely corresponds with the input data.
  // This saves loading the image to some memory first and then transferring
  // the data into the spline. Since the core is a vigra::MultiarrayView,
  // we can pass it to importImage as the desired target for loading the
  // image from disk.
                   
  vigra::importImage ( imageInfo , bspl.core ) ;
  
  // prefilter the b-spline
 
  bspl.prefilter() ;
  
  // now we construct the gsm evaluator
  
  ev_gsm ev ( bspl ) ;
  
  // this is where the result should go:
  
  target_type target ( imageInfo.shape() ) ;
 
  // now we obtain the result by performing a vspline::transform. This function
  // successively passes discrete coordinates into the target to the evaluator
  // it's invoked with, storing the result of the evaluator's evaluation
  // at the self-same coordinates. This is done multithreaded and vectorized
  // automatically, so it's very convenient, if the evaluator is at hand.
  // So here we have invested moderately more coding effort in the evaluator
  // and are rewarded with being able to use the evaluator with vspline's
  // high-level code for a fast implementation of our gsm problem.
  // The difference is quite noticeable. On my system, processing a full HD
  // image, I get:
  
  //   $ time ./gsm image.jpg
  //
  // real    0m0.838s
  // user    0m0.824s
  // sys     0m0.012s
  //
  //   $ time ./gsm2 image.jpg
  // 
  // real    0m0.339s
  // user    0m0.460s
  // sys     0m0.016s
 
  vspline::transform ( ev , target ) ;
 
  // store the result with vigra impex
 
  vigra::ImageExportInfo eximageInfo ( "gsm2.tif" );
  
  std::cout << "storing the target image as 'gsm2.tif'" << std::endl ;
  
  vigra::exportImage ( target ,
                       eximageInfo
                       .setPixelType("UINT8") ) ;
  
  exit ( 0 ) ;
}