[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/multi_convolution.hxx
00001 //-- -*- c++ -*-
00002 /************************************************************************/
00003 /*                                                                      */
00004 /*               Copyright 2003 by Christian-Dennis Rahn                */
00005 /*                        and Ullrich Koethe                            */
00006 /*                                                                      */
00007 /*    This file is part of the VIGRA computer vision library.           */
00008 /*    The VIGRA Website is                                              */
00009 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00010 /*    Please direct questions, bug reports, and contributions to        */
00011 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00012 /*        vigra@informatik.uni-hamburg.de                               */
00013 /*                                                                      */
00014 /*    Permission is hereby granted, free of charge, to any person       */
00015 /*    obtaining a copy of this software and associated documentation    */
00016 /*    files (the "Software"), to deal in the Software without           */
00017 /*    restriction, including without limitation the rights to use,      */
00018 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00019 /*    sell copies of the Software, and to permit persons to whom the    */
00020 /*    Software is furnished to do so, subject to the following          */
00021 /*    conditions:                                                       */
00022 /*                                                                      */
00023 /*    The above copyright notice and this permission notice shall be    */
00024 /*    included in all copies or substantial portions of the             */
00025 /*    Software.                                                         */
00026 /*                                                                      */
00027 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00028 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00029 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00030 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00031 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00032 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00033 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00034 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */                
00035 /*                                                                      */
00036 /************************************************************************/
00037 
00038 #ifndef VIGRA_MULTI_CONVOLUTION_H
00039 #define VIGRA_MULTI_CONVOLUTION_H
00040 
00041 #include "separableconvolution.hxx"
00042 #include "array_vector.hxx"
00043 #include "multi_array.hxx"
00044 #include "accessor.hxx"
00045 #include "numerictraits.hxx"
00046 #include "navigator.hxx"
00047 #include "metaprogramming.hxx"
00048 #include "multi_pointoperators.hxx"
00049 #include "functorexpression.hxx"
00050 
00051 namespace vigra
00052 {
00053 
00054 
00055 namespace detail
00056 {
00057 
00058 /********************************************************/
00059 /*                                                      */
00060 /*        internalSeparableConvolveMultiArray           */
00061 /*                                                      */
00062 /********************************************************/
00063 
00064 template <class SrcIterator, class SrcShape, class SrcAccessor,
00065           class DestIterator, class DestAccessor, class KernelIterator>
00066 void
00067 internalSeparableConvolveMultiArrayTmp(
00068                       SrcIterator si, SrcShape const & shape, SrcAccessor src,
00069                       DestIterator di, DestAccessor dest, KernelIterator kit)
00070 {
00071     enum { N = 1 + SrcIterator::level };
00072 
00073     typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
00074 
00075     // temporay array to hold the current line to enable in-place operation
00076     ArrayVector<TmpType> tmp( shape[0] );
00077 
00078     typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
00079     typedef MultiArrayNavigator<DestIterator, N> DNavigator;
00080 
00081     { // only operate on first dimension here
00082         SNavigator snav( si, shape, 0 );
00083         DNavigator dnav( di, shape, 0 );
00084 
00085         for( ; snav.hasMore(); snav++, dnav++ )
00086         {
00087              // first copy source to temp for maximum cache efficiency
00088              copyLine( snav.begin(), snav.end(), src,
00089                        tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
00090 
00091              convolveLine( srcIterRange(tmp.begin(), tmp.end(),
00092                                         typename AccessorTraits<TmpType>::default_const_accessor()),
00093                            destIter( dnav.begin(), dest ),
00094                            kernel1d( *kit ) );
00095         }
00096         ++kit;
00097     }
00098 
00099     // operate on further dimensions
00100     for( int d = 1; d < N; ++d, ++kit )
00101     {
00102         DNavigator dnav( di, shape, d );
00103 
00104         tmp.resize( shape[d] );
00105 
00106         for( ; dnav.hasMore(); dnav++ )
00107         {
00108              // first copy source to temp for maximum cache efficiency
00109              copyLine( dnav.begin(), dnav.end(), dest,
00110                        tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
00111 
00112              convolveLine( srcIterRange(tmp.begin(), tmp.end(),
00113                                         typename AccessorTraits<TmpType>::default_const_accessor()),
00114                            destIter( dnav.begin(), dest ),
00115                            kernel1d( *kit ) );
00116         }
00117     }
00118 }
00119 
00120 
00121 } // namespace detail
00122 
00123 /** \addtogroup MultiArrayConvolutionFilters Convolution filters for multi-dimensional arrays.
00124 
00125     These functions realize a separable convolution on an arbitrary dimensional
00126     array that is specified by iterators (compatible to \ref MultiIteratorPage)
00127     and shape objects. It can therefore be applied to a wide range of data structures
00128     (\ref vigra::MultiArrayView, \ref vigra::MultiArray etc.).
00129 */
00130 //@{
00131 
00132 /********************************************************/
00133 /*                                                      */
00134 /*             separableConvolveMultiArray              */
00135 /*                                                      */
00136 /********************************************************/
00137 
00138 /** \brief Separated convolution on multi-dimensional arrays.
00139 
00140     This function computes a separated convolution on all dimensions
00141     of the given multi-dimensional array. Both source and destination
00142     arrays are represented by iterators, shape objects and accessors.
00143     The destination array is required to already have the correct size.
00144 
00145     There are two variants of this functions: one takes a single kernel
00146     of type \ref vigra::Kernel1D which is then applied to all dimensions,
00147     whereas the other requires an iterator referencing a sequence of
00148     \ref vigra::Kernel1D objects, one for every dimension of the data.
00149     Then the first kernel in this sequence is applied to the innermost
00150     dimension (e.g. the x-dimension of an image), while the last is applied to the
00151     outermost dimension (e.g. the z-dimension in a 3D image).
00152 
00153     This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
00154     A full-sized internal array is only allocated if working on the destination
00155     array directly would cause round-off errors (i.e. if
00156     <tt>typeid(typename NumericTraits<typename DestAccessor::value_type>::RealPromote)
00157     != typeid(typename DestAccessor::value_type)</tt>.
00158 
00159     <b> Declarations:</b>
00160 
00161     pass arguments explicitly:
00162     \code
00163     namespace vigra {
00164         // apply the same kernel to all dimensions
00165         template <class SrcIterator, class SrcShape, class SrcAccessor,
00166                   class DestIterator, class DestAccessor, class T>
00167         void
00168         separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00169                                     DestIterator diter, DestAccessor dest,
00170                                     Kernel1D<T> const & kernel);
00171 
00172         // apply each kernel from the sequence 'kernels' in turn
00173         template <class SrcIterator, class SrcShape, class SrcAccessor,
00174                   class DestIterator, class DestAccessor, class KernelIterator>
00175         void
00176         separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00177                                     DestIterator diter, DestAccessor dest,
00178                                     KernelIterator kernels);
00179     }
00180     \endcode
00181 
00182     use argument objects in conjunction with \ref ArgumentObjectFactories :
00183     \code
00184     namespace vigra {
00185         // apply the same kernel to all dimensions
00186         template <class SrcIterator, class SrcShape, class SrcAccessor,
00187                   class DestIterator, class DestAccessor, class T>
00188         void
00189         separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00190                                     pair<DestIterator, DestAccessor> const & dest,
00191                                     Kernel1D<T> const & kernel);
00192 
00193         // apply each kernel from the sequence 'kernels' in turn
00194         template <class SrcIterator, class SrcShape, class SrcAccessor,
00195                   class DestIterator, class DestAccessor, class KernelIterator>
00196         void
00197         separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00198                                     pair<DestIterator, DestAccessor> const & dest,
00199                                     KernelIterator kernels);
00200     }
00201     \endcode
00202 
00203     <b> Usage:</b>
00204 
00205     <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
00206 
00207     \code
00208     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
00209     MultiArray<3, unsigned char> source(shape);
00210     MultiArray<3, float> dest(shape);
00211     ...
00212     Kernel1D<float> gauss;
00213     gauss.initGaussian(sigma);
00214     // create 3 Gauss kernels, one for each dimension
00215     ArrayVector<Kernel1D<float> > kernels(3, gauss);
00216 
00217     // perform Gaussian smoothing on all dimensions
00218     separableConvolveMultiArray(srcMultiArrayRange(source), destMultiArray(dest), 
00219                                 kernels.begin());
00220     \endcode
00221 
00222     \see vigra::Kernel1D, convolveLine()
00223 */
00224 doxygen_overloaded_function(template <...> void separableConvolveMultiArray)
00225 
00226 template <class SrcIterator, class SrcShape, class SrcAccessor,
00227           class DestIterator, class DestAccessor, class KernelIterator>
00228 void
00229 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
00230                              DestIterator d, DestAccessor dest, KernelIterator kernels )
00231 {
00232     typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
00233 
00234     if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
00235     {
00236         // need a temporary array to avoid rounding errors
00237         MultiArray<SrcShape::static_size, TmpType> tmpArray(shape);
00238         detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
00239              tmpArray.traverser_begin(), typename AccessorTraits<TmpType>::default_accessor(), kernels );
00240         copyMultiArray(srcMultiArrayRange(tmpArray), destIter(d, dest));
00241     }
00242     else
00243     {
00244         // work directly on the destination array
00245         detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
00246     }
00247 }
00248 
00249 template <class SrcIterator, class SrcShape, class SrcAccessor,
00250           class DestIterator, class DestAccessor, class KernelIterator>
00251 inline
00252 void separableConvolveMultiArray(
00253     triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00254     pair<DestIterator, DestAccessor> const & dest, KernelIterator kit )
00255 {
00256     separableConvolveMultiArray( source.first, source.second, source.third,
00257                                  dest.first, dest.second, kit );
00258 }
00259 
00260 template <class SrcIterator, class SrcShape, class SrcAccessor,
00261           class DestIterator, class DestAccessor, class T>
00262 inline void
00263 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
00264                              DestIterator d, DestAccessor dest,
00265                              Kernel1D<T> const & kernel )
00266 {
00267     ArrayVector<Kernel1D<T> > kernels(shape.size(), kernel);
00268 
00269     separableConvolveMultiArray( s, shape, src, d, dest, kernels.begin() );
00270 }
00271 
00272 template <class SrcIterator, class SrcShape, class SrcAccessor,
00273           class DestIterator, class DestAccessor, class T>
00274 inline void
00275 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00276                             pair<DestIterator, DestAccessor> const & dest,
00277                             Kernel1D<T> const & kernel )
00278 {
00279     ArrayVector<Kernel1D<T> > kernels(source.second.size(), kernel);
00280 
00281     separableConvolveMultiArray( source.first, source.second, source.third,
00282                                  dest.first, dest.second, kernels.begin() );
00283 }
00284 
00285 /********************************************************/
00286 /*                                                      */
00287 /*            convolveMultiArrayOneDimension            */
00288 /*                                                      */
00289 /********************************************************/
00290 
00291 /** \brief Convolution along a single dimension of a multi-dimensional arrays.
00292 
00293     This function computes a convolution along one dimension (specified by
00294     the parameter <tt>dim</tt> of the given multi-dimensional array with the given
00295     <tt>kernel</tt>. Both source and destination arrays are represented by
00296     iterators, shape objects and accessors. The destination array is required to
00297     already have the correct size.
00298 
00299     This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
00300 
00301     <b> Declarations:</b>
00302 
00303     pass arguments explicitly:
00304     \code
00305     namespace vigra {
00306         template <class SrcIterator, class SrcShape, class SrcAccessor,
00307                   class DestIterator, class DestAccessor, class T>
00308         void
00309         convolveMultiArrayOneDimension(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00310                                        DestIterator diter, DestAccessor dest,
00311                                        unsigned int dim, vigra::Kernel1D<T> const & kernel);
00312     }
00313     \endcode
00314 
00315     use argument objects in conjunction with \ref ArgumentObjectFactories :
00316     \code
00317     namespace vigra {
00318         template <class SrcIterator, class SrcShape, class SrcAccessor,
00319                   class DestIterator, class DestAccessor, class T>
00320         void
00321         convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00322                                        pair<DestIterator, DestAccessor> const & dest,
00323                                        unsigned int dim, vigra::Kernel1D<T> const & kernel);
00324     }
00325     \endcode
00326 
00327     <b> Usage:</b>
00328 
00329     <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
00330 
00331     \code
00332     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
00333     MultiArray<3, unsigned char> source(shape);
00334     MultiArray<3, float> dest(shape);
00335     ...
00336     Kernel1D<float> gauss;
00337     gauss.initGaussian(sigma);
00338 
00339     // perform Gaussian smoothing along dimensions 1 (height)
00340     convolveMultiArrayOneDimension(srcMultiArrayRange(source), destMultiArray(dest), 1, gauss);
00341     \endcode
00342 
00343     \see separableConvolveMultiArray()
00344 */
00345 doxygen_overloaded_function(template <...> void convolveMultiArrayOneDimension)
00346 
00347 template <class SrcIterator, class SrcShape, class SrcAccessor,
00348           class DestIterator, class DestAccessor, class T>
00349 void
00350 convolveMultiArrayOneDimension(SrcIterator s, SrcShape const & shape, SrcAccessor src,
00351                                DestIterator d, DestAccessor dest,
00352                                unsigned int dim, vigra::Kernel1D<T> const & kernel )
00353 {
00354     enum { N = 1 + SrcIterator::level };
00355     vigra_precondition( dim < N,
00356                         "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller "
00357                         "than the data dimensionality" );
00358 
00359     typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
00360     ArrayVector<TmpType> tmp( shape[dim] );
00361 
00362     typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
00363     typedef MultiArrayNavigator<DestIterator, N> DNavigator;
00364 
00365     SNavigator snav( s, shape, dim );
00366     DNavigator dnav( d, shape, dim );
00367 
00368     for( ; snav.hasMore(); snav++, dnav++ )
00369     {
00370          // first copy source to temp for maximum cache efficiency
00371          copyLine( snav.begin(), snav.end(), src,
00372            tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
00373 
00374          convolveLine( srcIterRange( tmp.begin(), tmp.end(), typename AccessorTraits<TmpType>::default_const_accessor()),
00375                        destIter( dnav.begin(), dest ),
00376                        kernel1d( kernel ) );
00377     }
00378 }
00379 
00380 template <class SrcIterator, class SrcShape, class SrcAccessor,
00381           class DestIterator, class DestAccessor, class T>
00382 inline void
00383 convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00384                                pair<DestIterator, DestAccessor> const & dest,
00385                                unsigned int dim, vigra::Kernel1D<T> const & kernel )
00386 {
00387     convolveMultiArrayOneDimension( source.first, source.second, source.third,
00388                                    dest.first, dest.second, dim, kernel );
00389 }
00390 
00391 /********************************************************/
00392 /*                                                      */
00393 /*             gaussianSmoothMultiArray                 */
00394 /*                                                      */
00395 /********************************************************/
00396 
00397 /** \brief Isotropic Gaussian smoothing of a multi-dimensional arrays.
00398 
00399     This function computes an isotropic convolution of the given multi-dimensional
00400     array with a Gaussian filter at the given standard deviation <tt>sigma</tt>.
00401     Both source and destination arrays are represented by
00402     iterators, shape objects and accessors. The destination array is required to
00403     already have the correct size. This function may work in-place, which means
00404     that <tt>siter == diter</tt> is allowed. It is implemented by a call to
00405     \ref separableConvolveMultiArray() with the appropriate kernel.
00406     If the data are anisotropic (different pixel size along different dimensions)
00407     you should call \ref separableConvolveMultiArray() directly with the appropriate
00408     anisotropic Gaussians.
00409 
00410     <b> Declarations:</b>
00411 
00412     pass arguments explicitly:
00413     \code
00414     namespace vigra {
00415         template <class SrcIterator, class SrcShape, class SrcAccessor,
00416                   class DestIterator, class DestAccessor>
00417         void
00418         gaussianSmoothMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00419                                  DestIterator diter, DestAccessor dest,
00420                                  double sigma);
00421     }
00422     \endcode
00423 
00424     use argument objects in conjunction with \ref ArgumentObjectFactories :
00425     \code
00426     namespace vigra {
00427         template <class SrcIterator, class SrcShape, class SrcAccessor,
00428                   class DestIterator, class DestAccessor>
00429         void
00430         gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00431                                  pair<DestIterator, DestAccessor> const & dest,
00432                                  double sigma);
00433     }
00434     \endcode
00435 
00436     <b> Usage:</b>
00437 
00438     <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
00439 
00440     \code
00441     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
00442     MultiArray<3, unsigned char> source(shape);
00443     MultiArray<3, float> dest(shape);
00444     ...
00445     // perform isotropic Gaussian smoothing at scale 'sigma'
00446     gaussianSmoothMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
00447     \endcode
00448 
00449     \see separableConvolveMultiArray()
00450 */
00451 doxygen_overloaded_function(template <...> void gaussianSmoothMultiArray)
00452 
00453 template <class SrcIterator, class SrcShape, class SrcAccessor,
00454           class DestIterator, class DestAccessor>
00455 void
00456 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
00457                    DestIterator d, DestAccessor dest, double sigma )
00458 {
00459     Kernel1D<double> gauss;
00460     gauss.initGaussian( sigma );
00461 
00462     separableConvolveMultiArray( s, shape, src, d, dest, gauss);
00463 }
00464 
00465 template <class SrcIterator, class SrcShape, class SrcAccessor,
00466           class DestIterator, class DestAccessor>
00467 inline void
00468 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00469                   pair<DestIterator, DestAccessor> const & dest,
00470                   double sigma )
00471 {
00472     gaussianSmoothMultiArray( source.first, source.second, source.third,
00473                               dest.first, dest.second, sigma );
00474 }
00475 
00476 /********************************************************/
00477 /*                                                      */
00478 /*             gaussianGradientMultiArray               */
00479 /*                                                      */
00480 /********************************************************/
00481 
00482 /** \brief Calculate Gaussian gradient of a multi-dimensional arrays.
00483 
00484     This function computes the Gaussian gradient of the given multi-dimensional
00485     array with a sequence of first-derivative-of-Gaussian filters at the given
00486     standard deviation <tt>sigma</tt> (differentiation is applied to each dimension
00487     in turn, starting with the innermost dimension). Both source and destination arrays
00488     are represented by iterators, shape objects and accessors. The destination array is
00489     required to have a vector valued pixel type with as many elements as the number of
00490     dimensions. This function is implemented by calls to
00491     \ref separableConvolveMultiArray() with the appropriate kernels.
00492     If the data are anisotropic (different pixel size along different dimensions)
00493     you should call \ref separableConvolveMultiArray() directly with the appropriate
00494     anisotropic Gaussian derivatives.
00495 
00496     <b> Declarations:</b>
00497 
00498     pass arguments explicitly:
00499     \code
00500     namespace vigra {
00501         template <class SrcIterator, class SrcShape, class SrcAccessor,
00502                   class DestIterator, class DestAccessor>
00503         void
00504         gaussianGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00505                                    DestIterator diter, DestAccessor dest,
00506                                    double sigma);
00507     }
00508     \endcode
00509 
00510     use argument objects in conjunction with \ref ArgumentObjectFactories :
00511     \code
00512     namespace vigra {
00513         template <class SrcIterator, class SrcShape, class SrcAccessor,
00514                   class DestIterator, class DestAccessor>
00515         void
00516         gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00517                                    pair<DestIterator, DestAccessor> const & dest,
00518                                    double sigma);
00519     }
00520     \endcode
00521 
00522     <b> Usage:</b>
00523 
00524     <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
00525 
00526     \code
00527     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
00528     MultiArray<3, unsigned char> source(shape);
00529     MultiArray<3, TinyVector<float, 3> > dest(shape);
00530     ...
00531     // compute Gaussian gradient at scale sigma
00532     gaussianGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
00533     \endcode
00534 
00535     <b> Required Interface:</b>
00536 
00537     see \ref separableConvolveMultiArray(), in addition:
00538 
00539     \code
00540     int dimension = 0;
00541     VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
00542     \endcode
00543 
00544     \see separableConvolveMultiArray()
00545 */
00546 doxygen_overloaded_function(template <...> void gaussianGradientMultiArray)
00547 
00548 template <class SrcIterator, class SrcShape, class SrcAccessor,
00549           class DestIterator, class DestAccessor>
00550 void
00551 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
00552                            DestIterator di, DestAccessor dest, double sigma )
00553 {
00554     typedef typename DestAccessor::value_type DestType;
00555     typedef typename DestType::value_type     DestValueType;
00556     typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
00557    
00558     static const int N = SrcShape::static_size;
00559 
00560     for(int k=0; k<N; ++k)
00561         if(shape[k] <=0)
00562             return;
00563 
00564     vigra_precondition(N == dest.size(di),
00565         "gaussianGradientMultiArray(): Wrong number of channels in output array.");
00566 
00567     vigra_precondition(sigma > 0.0, "gaussianGradientMultiArray(): Scale must be positive.");
00568 
00569     Kernel1D<KernelType> gauss, derivative;
00570     gauss.initGaussian(sigma);
00571 
00572     typedef VectorElementAccessor<DestAccessor> ElementAccessor;
00573 
00574     // compute gradient components
00575     for(int d = 0; d < N; ++d )
00576     {
00577         ArrayVector<Kernel1D<KernelType> > kernels(N, gauss);
00578         kernels[d].initGaussianDerivative(sigma, 1);
00579         separableConvolveMultiArray( si, shape, src, di, ElementAccessor(d, dest), kernels.begin());
00580     }
00581 }
00582 
00583 template <class SrcIterator, class SrcShape, class SrcAccessor,
00584           class DestIterator, class DestAccessor>
00585 inline void
00586 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00587                            pair<DestIterator, DestAccessor> const & dest, double sigma )
00588 {
00589     gaussianGradientMultiArray( source.first, source.second, source.third,
00590                                 dest.first, dest.second, sigma );
00591 }
00592 
00593 /********************************************************/
00594 /*                                                      */
00595 /*             symmetricGradientMultiArray              */
00596 /*                                                      */
00597 /********************************************************/
00598 
00599 /** \brief Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
00600 
00601     This function computes the gradient of the given multi-dimensional
00602     array with a sequence of symmetric difference filters a (differentiation is applied
00603     to each dimension in turn, starting with the innermost dimension). Both source and
00604     destination arrays are represented by iterators, shape objects and accessors.
00605     The destination array is required to have a vector valued pixel type with as many
00606     elements as the number of dimensions. This function is implemented by calls to
00607     \ref convolveMultiArrayOneDimension() with the symmetric difference kernel.
00608 
00609     <b> Declarations:</b>
00610 
00611     pass arguments explicitly:
00612     \code
00613     namespace vigra {
00614         template <class SrcIterator, class SrcShape, class SrcAccessor,
00615                   class DestIterator, class DestAccessor>
00616         void
00617         symmetricGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00618                                     DestIterator diter, DestAccessor dest);
00619     }
00620     \endcode
00621 
00622     use argument objects in conjunction with \ref ArgumentObjectFactories :
00623     \code
00624     namespace vigra {
00625         template <class SrcIterator, class SrcShape, class SrcAccessor,
00626                   class DestIterator, class DestAccessor>
00627         void
00628         symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00629                                     pair<DestIterator, DestAccessor> const & dest);
00630     }
00631     \endcode
00632 
00633     <b> Usage:</b>
00634 
00635     <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
00636 
00637     \code
00638     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
00639     MultiArray<3, unsigned char> source(shape);
00640     MultiArray<3, TinyVector<float, 3> > dest(shape);
00641     ...
00642     // compute gradient
00643     symmetricGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest));
00644     \endcode
00645 
00646     <b> Required Interface:</b>
00647 
00648     see \ref convolveMultiArrayOneDimension(), in addition:
00649 
00650     \code
00651     int dimension = 0;
00652     VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
00653     \endcode
00654 
00655     \see convolveMultiArrayOneDimension()
00656 */
00657 doxygen_overloaded_function(template <...> void symmetricGradientMultiArray)
00658 
00659 template <class SrcIterator, class SrcShape, class SrcAccessor,
00660           class DestIterator, class DestAccessor>
00661 void
00662 symmetricGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
00663                             DestIterator di, DestAccessor dest)
00664 {
00665     typedef typename DestAccessor::value_type DestType;
00666     typedef typename DestType::value_type     DestValueType;
00667     typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
00668 
00669     static const int N = SrcShape::static_size;
00670 
00671     for(int k=0; k<N; ++k)
00672         if(shape[k] <=0)
00673             return;
00674 
00675     vigra_precondition(N == dest.size(di),
00676         "symmetricGradientMultiArray(): Wrong number of channels in output array.");
00677 
00678     Kernel1D<KernelType> filter;
00679     filter.initSymmetricGradient();
00680 
00681     typedef VectorElementAccessor<DestAccessor> ElementAccessor;
00682 
00683     // compute gradient components
00684     for(int d = 0; d < N; ++d )
00685     {
00686         convolveMultiArrayOneDimension(si, shape, src,
00687                                        di, ElementAccessor(d, dest),
00688                                        d, filter);
00689     }
00690 }
00691 
00692 template <class SrcIterator, class SrcShape, class SrcAccessor,
00693           class DestIterator, class DestAccessor>
00694 inline void
00695 symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00696                             pair<DestIterator, DestAccessor> const & dest )
00697 {
00698     symmetricGradientMultiArray(source.first, source.second, source.third,
00699                                 dest.first, dest.second);
00700 }
00701 
00702 
00703 /********************************************************/
00704 /*                                                      */
00705 /*            laplacianOfGaussianMultiArray             */
00706 /*                                                      */
00707 /********************************************************/
00708 
00709 /** \brief Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
00710 
00711     This function computes the Laplacian the given N-dimensional
00712     array with a sequence of second-derivative-of-Gaussian filters at the given
00713     standard deviation <tt>sigma</tt>. Both source and destination arrays
00714     are represented by iterators, shape objects and accessors. Both source and destination 
00715     arrays must have scalar value_type. This function is implemented by calls to
00716     \ref separableConvolveMultiArray() with the appropriate kernels, followed by summation.
00717 
00718     <b> Declarations:</b>
00719 
00720     pass arguments explicitly:
00721     \code
00722     namespace vigra {
00723         template <class SrcIterator, class SrcShape, class SrcAccessor,
00724                   class DestIterator, class DestAccessor>
00725         void
00726         laplacianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00727                                       DestIterator diter, DestAccessor dest,
00728                                       double sigma);
00729     }
00730     \endcode
00731 
00732     use argument objects in conjunction with \ref ArgumentObjectFactories :
00733     \code
00734     namespace vigra {
00735         template <class SrcIterator, class SrcShape, class SrcAccessor,
00736                   class DestIterator, class DestAccessor>
00737         void
00738         laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00739                                       pair<DestIterator, DestAccessor> const & dest,
00740                                       double sigma);
00741     }
00742     \endcode
00743 
00744     <b> Usage:</b>
00745 
00746     <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
00747 
00748     \code
00749     MultiArray<3, float> source(shape);
00750     MultiArray<3, float> laplacian(shape);
00751     ...
00752     // compute Laplacian at scale sigma
00753     laplacianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(laplacian), sigma);
00754     \endcode
00755 
00756     <b> Required Interface:</b>
00757 
00758     see \ref separableConvolveMultiArray(), in addition:
00759 
00760     \code
00761     int dimension = 0;
00762     VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
00763     \endcode
00764 
00765     \see separableConvolveMultiArray()
00766 */
00767 doxygen_overloaded_function(template <...> void laplacianOfGaussianMultiArray)
00768 
00769 template <class SrcIterator, class SrcShape, class SrcAccessor,
00770           class DestIterator, class DestAccessor>
00771 void
00772 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
00773                               DestIterator di, DestAccessor dest, double sigma )
00774 { 
00775     using namespace functor;
00776     
00777     typedef typename DestAccessor::value_type DestType;
00778     typedef typename NumericTraits<DestType>::RealPromote KernelType;
00779     typedef typename AccessorTraits<KernelType>::default_accessor DerivativeAccessor;
00780 
00781     static const int N = SrcShape::static_size;
00782     
00783     vigra_precondition(sigma > 0.0, "laplacianOfGaussianMultiArray(): Scale must be positive.");
00784 
00785     Kernel1D<KernelType> gauss;
00786     gauss.initGaussian(sigma);
00787     
00788     MultiArray<N, KernelType> derivative(shape);
00789 
00790     // compute 2nd derivatives and sum them up
00791     for(int d = 0; d < N; ++d )
00792     {
00793         ArrayVector<Kernel1D<KernelType> > kernels(N, gauss);
00794         kernels[d].initGaussianDerivative(sigma, 2);
00795         if(d == 0)
00796         {
00797             separableConvolveMultiArray( si, shape, src, 
00798                                          di, dest, kernels.begin());
00799         }
00800         else
00801         {
00802             separableConvolveMultiArray( si, shape, src, 
00803                                          derivative.traverser_begin(), DerivativeAccessor(), 
00804                                          kernels.begin());
00805             combineTwoMultiArrays(di, shape, dest, derivative.traverser_begin(), DerivativeAccessor(), 
00806                                   di, dest, Arg1() + Arg2() );
00807         }
00808     }
00809 }
00810 
00811 template <class SrcIterator, class SrcShape, class SrcAccessor,
00812           class DestIterator, class DestAccessor>
00813 inline void
00814 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00815                             pair<DestIterator, DestAccessor> const & dest, double sigma )
00816 {
00817     laplacianOfGaussianMultiArray( source.first, source.second, source.third,
00818                                    dest.first, dest.second, sigma );
00819 }
00820 
00821 /********************************************************/
00822 /*                                                      */
00823 /*              hessianOfGaussianMultiArray             */
00824 /*                                                      */
00825 /********************************************************/
00826 
00827 /** \brief Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
00828 
00829     This function computes the Hessian matrix the given scalar N-dimensional
00830     array with a sequence of second-derivative-of-Gaussian filters at the given
00831     standard deviation <tt>sigma</tt>. Both source and destination arrays
00832     are represented by iterators, shape objects and accessors. The destination array must 
00833     have a vector valued element type with N*(N+1)/2 elements (it represents the
00834     upper triangular part of the symmetric Hessian matrix). This function is implemented by calls to
00835     \ref separableConvolveMultiArray() with the appropriate kernels.
00836 
00837     <b> Declarations:</b>
00838 
00839     pass arguments explicitly:
00840     \code
00841     namespace vigra {
00842         template <class SrcIterator, class SrcShape, class SrcAccessor,
00843                   class DestIterator, class DestAccessor>
00844         void
00845         hessianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00846                                     DestIterator diter, DestAccessor dest,
00847                                     double sigma);
00848     }
00849     \endcode
00850 
00851     use argument objects in conjunction with \ref ArgumentObjectFactories :
00852     \code
00853     namespace vigra {
00854         template <class SrcIterator, class SrcShape, class SrcAccessor,
00855                   class DestIterator, class DestAccessor>
00856         void
00857         hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00858                                     pair<DestIterator, DestAccessor> const & dest,
00859                                     double sigma);
00860     }
00861     \endcode
00862 
00863     <b> Usage:</b>
00864 
00865     <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
00866 
00867     \code
00868     MultiArray<3, float> source(shape);
00869     MultiArray<3, TinyVector<float, 6> > dest(shape);
00870     ...
00871     // compute Hessian at scale sigma
00872     hessianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
00873     \endcode
00874 
00875     <b> Required Interface:</b>
00876 
00877     see \ref separableConvolveMultiArray(), in addition:
00878 
00879     \code
00880     int dimension = 0;
00881     VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
00882     \endcode
00883 
00884     \see separableConvolveMultiArray(), vectorToTensorMultiArray()
00885 */
00886 doxygen_overloaded_function(template <...> void hessianOfGaussianMultiArray)
00887 
00888 template <class SrcIterator, class SrcShape, class SrcAccessor,
00889           class DestIterator, class DestAccessor>
00890 void
00891 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
00892                             DestIterator di, DestAccessor dest, double sigma )
00893 { 
00894     typedef typename DestAccessor::value_type DestType;
00895     typedef typename DestType::value_type     DestValueType;
00896     typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
00897 
00898     static const int N = SrcShape::static_size;
00899     static const int M = N*(N+1)/2;
00900     
00901     for(int k=0; k<N; ++k)
00902         if(shape[k] <=0)
00903             return;
00904 
00905     vigra_precondition(M == dest.size(di),
00906         "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
00907 
00908     vigra_precondition(sigma > 0.0, "hessianOfGaussianMultiArray(): Scale must be positive.");
00909 
00910     Kernel1D<KernelType> gauss;
00911     gauss.initGaussian(sigma);
00912 
00913     typedef VectorElementAccessor<DestAccessor> ElementAccessor;
00914 
00915     // compute elements of the Hessian matrix
00916     for(int b=0, i=0; i<N; ++i)
00917     {
00918         for(int j=i; j<N; ++j, ++b)
00919         {
00920             ArrayVector<Kernel1D<KernelType> > kernels(N, gauss);
00921             if(i == j)
00922             {
00923                 kernels[i].initGaussianDerivative(sigma, 2);
00924             }
00925             else
00926             {
00927                 kernels[i].initGaussianDerivative(sigma, 1);
00928                 kernels[j].initGaussianDerivative(sigma, 1);
00929             }
00930             separableConvolveMultiArray(si, shape, src, di, ElementAccessor(b, dest),
00931                                         kernels.begin());
00932         }
00933     }
00934 }
00935 
00936 template <class SrcIterator, class SrcShape, class SrcAccessor,
00937           class DestIterator, class DestAccessor>
00938 inline void
00939 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00940                             pair<DestIterator, DestAccessor> const & dest, double sigma )
00941 {
00942     hessianOfGaussianMultiArray( source.first, source.second, source.third,
00943                                  dest.first, dest.second, sigma );
00944 }
00945 
00946 namespace detail {
00947 
00948 template<int N, class VectorType>
00949 struct StructurTensorFunctor
00950 {
00951     typedef VectorType result_type;
00952     typedef typename VectorType::value_type ValueType;
00953     
00954     template <class T>
00955     VectorType operator()(T const & in) const
00956     {
00957         VectorType res;
00958         for(int b=0, i=0; i<N; ++i)
00959         {
00960             for(int j=i; j<N; ++j, ++b)
00961             {
00962                 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
00963             }
00964         }
00965         return res;
00966     }
00967 };
00968 
00969 } // namespace detail
00970 
00971 /********************************************************/
00972 /*                                                      */
00973 /*               structureTensorMultiArray              */
00974 /*                                                      */
00975 /********************************************************/
00976 
00977 /** \brief Calculate th structure tensor of a multi-dimensional arrays.
00978 
00979     This function computes the gradient (outer product) tensor for each element
00980     of the given N-dimensional array with first-derivative-of-Gaussian filters at 
00981     the given <tt>innerScale</tt>, followed by Gaussian smoothing at <tt>outerScale</tt>.
00982     Both source and destination arrays are represented by iterators, shape objects and 
00983     accessors. The destination array must have a vector valued pixel type with 
00984     N*(N+1)/2 elements (it represents the upper triangular part of the symmetric 
00985     structure tensor matrix). If the source array is also vector valued, the 
00986     resulting structure tensor is the sum of the individual tensors for each channel.
00987     This function is implemented by calls to
00988     \ref separableConvolveMultiArray() with the appropriate kernels.
00989 
00990     <b> Declarations:</b>
00991 
00992     pass arguments explicitly:
00993     \code
00994     namespace vigra {
00995         template <class SrcIterator, class SrcShape, class SrcAccessor,
00996                   class DestIterator, class DestAccessor>
00997         void
00998         structureTensorMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00999                                   DestIterator diter, DestAccessor dest,
01000                                   double innerScale, double outerScale);
01001     }
01002     \endcode
01003 
01004     use argument objects in conjunction with \ref ArgumentObjectFactories :
01005     \code
01006     namespace vigra {
01007         template <class SrcIterator, class SrcShape, class SrcAccessor,
01008                   class DestIterator, class DestAccessor>
01009         void
01010         structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01011                                   pair<DestIterator, DestAccessor> const & dest,
01012                                   double innerScale, double outerScale);
01013     }
01014     \endcode
01015 
01016     <b> Usage:</b>
01017 
01018     <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
01019 
01020     \code
01021     MultiArray<3, RGBValue<float> > source(shape);
01022     MultiArray<3, TinyVector<float, 6> > dest(shape);
01023     ...
01024     // compute structure tensor at scales innerScale and outerScale
01025     structureTensorMultiArray(srcMultiArrayRange(source), destMultiArray(dest), innerScale, outerScale);
01026     \endcode
01027 
01028     <b> Required Interface:</b>
01029 
01030     see \ref separableConvolveMultiArray(), in addition:
01031 
01032     \code
01033     int dimension = 0;
01034     VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
01035     \endcode
01036 
01037     \see separableConvolveMultiArray(), vectorToTensorMultiArray()
01038 */
01039 doxygen_overloaded_function(template <...> void structureTensorMultiArray)
01040 
01041 template <class SrcIterator, class SrcShape, class SrcAccessor,
01042           class DestIterator, class DestAccessor>
01043 void
01044 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
01045                           DestIterator di, DestAccessor dest, 
01046                           double innerScale, double outerScale)
01047 { 
01048     static const int N = SrcShape::static_size;
01049     static const int M = N*(N+1)/2;
01050     
01051     typedef typename DestAccessor::value_type DestType;
01052     typedef typename DestType::value_type     DestValueType;
01053     typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
01054     typedef TinyVector<KernelType, N> GradientVector;
01055     typedef typename AccessorTraits<GradientVector>::default_accessor GradientAccessor;
01056 
01057     for(int k=0; k<N; ++k)
01058         if(shape[k] <=0)
01059             return;
01060 
01061     vigra_precondition(M == dest.size(di),
01062         "structureTensorMultiArray(): Wrong number of channels in output array.");
01063 
01064     vigra_precondition(innerScale > 0.0 && outerScale >= 0.0,
01065          "structureTensorMultiArray(): Scale must be positive.");
01066 
01067     MultiArray<N, GradientVector> gradient(shape);
01068     gaussianGradientMultiArray(si, shape, src, 
01069                                gradient.traverser_begin(), GradientAccessor(), 
01070                                innerScale);
01071 
01072     transformMultiArray(gradient.traverser_begin(), shape, GradientAccessor(), 
01073                         di, dest, 
01074                         detail::StructurTensorFunctor<N, DestType>());
01075 
01076     gaussianSmoothMultiArray(di, shape, dest, di, dest, outerScale);
01077 }
01078 
01079 template <class SrcIterator, class SrcShape, class SrcAccessor,
01080           class DestIterator, class DestAccessor>
01081 inline void
01082 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01083                           pair<DestIterator, DestAccessor> const & dest, 
01084                           double innerScale, double outerScale)
01085 {
01086     structureTensorMultiArray( source.first, source.second, source.third,
01087                                dest.first, dest.second, innerScale, outerScale );
01088 }
01089 
01090 //@}
01091 
01092 } //-- namespace vigra
01093 
01094 
01095 #endif        //-- VIGRA_MULTI_CONVOLUTION_H

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.7.0 (Thu Aug 25 2011)