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

vigra/multi_convolution.hxx VIGRA

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 #include "tinyvector.hxx"
00051 #include "algorithm.hxx"
00052 
00053 namespace vigra
00054 {
00055 
00056 namespace detail
00057 {
00058 
00059 struct DoubleYielder
00060 {
00061     const double value;
00062     DoubleYielder(double v, unsigned, const char *const) : value(v) {}
00063     DoubleYielder(double v)                              : value(v) {}
00064     void operator++() {}
00065     double operator*() const { return value; }
00066 };
00067 
00068 template <typename X>
00069 struct IteratorDoubleYielder
00070 {
00071     X it;
00072     IteratorDoubleYielder(X i, unsigned, const char *const) : it(i) {}
00073     IteratorDoubleYielder(X i)                              : it(i) {}
00074     void operator++() { ++it; }
00075     double operator*() const { return *it; }
00076 };
00077 
00078 template <typename X>
00079 struct SequenceDoubleYielder
00080 {
00081     typename X::const_iterator it;
00082     SequenceDoubleYielder(const X & seq, unsigned dim,
00083                           const char *const function_name = "SequenceDoubleYielder")
00084         : it(seq.begin())
00085     {
00086         if (seq.size() == dim)
00087             return;
00088         std::string msg = "(): Parameter number be equal to the number of spatial dimensions.";
00089         vigra_precondition(false, function_name + msg);
00090     }
00091     void operator++() { ++it; }
00092     double operator*() const { return *it; }
00093 };
00094 
00095 template <typename X>
00096 struct WrapDoubleIterator
00097 {
00098     typedef
00099     typename IfBool< IsConvertibleTo<X, double>::value,
00100         DoubleYielder,
00101         typename IfBool< IsIterator<X>::value || IsArray<X>::value,
00102             IteratorDoubleYielder<X>,
00103             SequenceDoubleYielder<X>
00104         >::type
00105     >::type type;
00106 };
00107 
00108 template <class Param1, class Param2, class Param3>
00109 struct WrapDoubleIteratorTriple
00110 {
00111     typename WrapDoubleIterator<Param1>::type sigma_eff_it;
00112     typename WrapDoubleIterator<Param2>::type sigma_d_it;
00113     typename WrapDoubleIterator<Param3>::type step_size_it;
00114     WrapDoubleIteratorTriple(Param1 sigma_eff, Param2 sigma_d, Param3 step_size)
00115         : sigma_eff_it(sigma_eff), sigma_d_it(sigma_d), step_size_it(step_size) {}
00116     void operator++()
00117     {
00118         ++sigma_eff_it;
00119         ++sigma_d_it;
00120         ++step_size_it;
00121     }
00122     double sigma_eff() const { return *sigma_eff_it; }
00123     double sigma_d() const { return *sigma_d_it; }
00124     double step_size() const { return *step_size_it; }
00125     static double sqr(double x) { return x * x; }
00126     static void sigma_precondition(double sigma, const char *const function_name)
00127     {
00128         if (sigma < 0.0)
00129         {
00130              std::string msg = "(): Scale must be positive.";
00131              vigra_precondition(false, function_name + msg);
00132         }
00133     }
00134     double sigma_scaled(const char *const function_name = "unknown function ") const
00135     {
00136         sigma_precondition(sigma_eff(), function_name);
00137         sigma_precondition(sigma_d(), function_name);
00138         double sigma_squared = sqr(sigma_eff()) - sqr(sigma_d());
00139         if (sigma_squared > 0.0)
00140         {
00141             return std::sqrt(sigma_squared) / step_size();
00142         }
00143         else
00144         {
00145              std::string msg = "(): Scale would be imaginary or zero.";
00146              vigra_precondition(false, function_name + msg);
00147              return 0;
00148         }
00149     }
00150 };
00151 
00152 template <unsigned dim>
00153 struct multiArrayScaleParam
00154 {
00155     typedef TinyVector<double, dim> p_vector;
00156     typedef typename p_vector::const_iterator return_type;
00157     p_vector vec;
00158 
00159     template <class Param>
00160     multiArrayScaleParam(Param val, const char *const function_name = "multiArrayScaleParam")
00161     {
00162         typename WrapDoubleIterator<Param>::type in(val, dim, function_name);
00163         for (unsigned i = 0; i != dim; ++i, ++in)
00164             vec[i] = *in;
00165     }
00166     return_type operator()() const
00167     {
00168         return vec.begin();
00169     }
00170     static void precondition(unsigned n_par, const char *const function_name = "multiArrayScaleParam")
00171     {
00172         char n[3] = "0.";
00173         n[0] += dim;
00174         std::string msg = "(): dimension parameter must be ";
00175         vigra_precondition(dim == n_par, function_name + msg + n);
00176     }
00177     multiArrayScaleParam(double v0, double v1, const char *const function_name = "multiArrayScaleParam")
00178     {
00179         precondition(2, function_name);
00180         vec = p_vector(v0, v1);
00181     }
00182     multiArrayScaleParam(double v0, double v1, double v2, const char *const function_name = "multiArrayScaleParam")
00183     {
00184         precondition(3, function_name);
00185         vec = p_vector(v0, v1, v2);
00186     }
00187     multiArrayScaleParam(double v0, double v1, double v2,  double v3, const char *const function_name = "multiArrayScaleParam")
00188     {
00189         precondition(4, function_name);
00190         vec = p_vector(v0, v1, v2, v3);
00191     }
00192     multiArrayScaleParam(double v0, double v1, double v2,  double v3, double v4, const char *const function_name = "multiArrayScaleParam")
00193     {
00194         precondition(5, function_name);
00195         vec = p_vector(v0, v1, v2, v3, v4);
00196     }
00197 };
00198 
00199 } // namespace detail
00200 
00201 #define VIGRA_CONVOLUTION_OPTIONS(function_name, default_value, member_name) \
00202     template <class Param> \
00203     ConvolutionOptions & function_name(const Param & val) \
00204     { \
00205         member_name = ParamVec(val, "ConvolutionOptions::" #function_name); \
00206         return *this; \
00207     } \
00208     ConvolutionOptions & function_name() \
00209     { \
00210         member_name = ParamVec(default_value, "ConvolutionOptions::" #function_name); \
00211         return *this; \
00212     } \
00213     ConvolutionOptions & function_name(double v0, double v1) \
00214     { \
00215         member_name = ParamVec(v0, v1, "ConvolutionOptions::" #function_name); \
00216         return *this; \
00217     } \
00218     ConvolutionOptions & function_name(double v0, double v1, double v2) \
00219     { \
00220         member_name = ParamVec(v0, v1, v2, "ConvolutionOptions::" #function_name); \
00221         return *this; \
00222     } \
00223     ConvolutionOptions & function_name(double v0, double v1, double v2, double v3) \
00224     { \
00225         member_name = ParamVec(v0, v1, v2, v3, "ConvolutionOptions::" #function_name); \
00226         return *this; \
00227     } \
00228     ConvolutionOptions & function_name(double v0, double v1, double v2, double v3, double v4) \
00229     { \
00230         member_name = ParamVec(v0, v1, v2, v3, v4, "ConvolutionOptions::" #function_name); \
00231         return *this; \
00232     }
00233 
00234 
00235 /** \brief  Options class template for convolutions.
00236  
00237   <b>\#include</b> <vigra/multi_convolution.hxx>
00238   
00239   This class enables the calculation of scale space convolutions
00240   such as \ref gaussianGradientMultiArray() on data with anisotropic
00241   discretization. For these, the result of the ordinary calculation
00242   has to be multiplied by factors of \f$1/w^{n}\f$ for each dimension,
00243   where \f$w\f$ is the step size of the grid in said dimension and
00244   \f$n\f$ is the differential order of the convolution, e.g., 1 for
00245   gaussianGradientMultiArray(), and 0 for gaussianSmoothMultiArray(),
00246   respectively. Also for each dimension in turn, the convolution's scale
00247   parameter \f$\sigma\f$ has to be replaced by
00248   \f$\sqrt{\sigma_\mathrm{eff}^2 - \sigma_\mathrm{D}^2}\Big/w\f$,
00249   where \f$\sigma_\mathrm{eff}\f$ is the resulting effective filtering
00250   scale. The data is assumed to be already filtered by a 
00251   gaussian smoothing with the scale parameter \f$\sigma_\mathrm{D}\f$
00252   (such as by measuring equipment). All of the above changes are
00253   automatically employed by the convolution functions for <tt>MultiArray</tt>s
00254   if a corresponding options object is provided.
00255 
00256   The <tt>ConvolutionOptions</tt> class must be parameterized by the dimension
00257   <tt>dim</tt>
00258   of the <tt>MultiArray</tt>s on which it is used. The actual per-axis
00259   options are set by (overloaded) member functions explained below,
00260   or else default to neutral values corresponding to the absence of the
00261   particular option.
00262   
00263   All member functions set <tt>dim</tt> values of the respective convolution
00264   option, one for each dimension. They may be set explicitly by multiple
00265   arguments for up to five dimensions, or by a single argument to the same
00266   value for all dimensions. For the general case, a single argument that is
00267   either a C-syle array, an iterator, or a C++ standard library style
00268   sequence (such as <tt>std::vector</tt>, with member functions <tt>begin()</tt>
00269   and <tt>size()</tt>) supplies the option values for any number of dimensions.
00270   
00271   Note that the return value of all member functions is <tt>*this</tt>, which
00272   provides the mechanism for concatenating member function calls as shown below.
00273 
00274   <b>usage with explicit parameters:</b>
00275 
00276   \code
00277   ConvolutionOptions<2> opt = ConvolutionOptions<2>().stepSize(1, 2.3);
00278   \endcode
00279  
00280   <b>usage with arrays:</b>
00281  
00282   \code
00283   const double step_size[3] = { x_scale, y_scale, z_scale };
00284   ConvolutionOptions<3> opt = ConvolutionOptions<3>().stepSize(step_size);
00285   \endcode
00286 
00287   <b>usage with C++ standard library style sequences:</b>
00288  
00289   \code
00290   TinyVector<double, 4> step_size(1, 1, 2.0, 1.5);
00291   TinyVector<double, 4>  r_sigmas(1, 1, 2.3, 3.2);
00292   ConvolutionOptions<4> opt = ConvolutionOptions<4>().stepSize(step_size).resolutionStdDev(r_sigmas);
00293   \endcode
00294 
00295   <b>usage with iterators:</b>
00296 
00297   \code
00298   ArrayVector<double> step_size;
00299   step_size.push_back(0);
00300   step_size.push_back(3);
00301   step_size.push_back(4);
00302   ArrayVector<double>::iterator i = step_size.begin();
00303   ++i;
00304   ConvolutionOptions<2> opt = ConvolutionOptions<2>().stepSize(i);
00305   \endcode
00306 
00307   <b>general usage in a convolution function call:</b>
00308 
00309   \code
00310   MultiArray<3, double> test_image;
00311   MultiArray<3, double> out_image;
00312   gaussianSmoothMultiArray(srcMultiArrayRange(test_image),
00313                            destMultiArray(out_image),
00314                            5.0,
00315                            ConvolutionOptions<3>()
00316                               .stepSize        (1, 1, 3.2)
00317                               .resolutionStdDev(1, 1, 4)
00318                           );
00319   \endcode
00320  
00321 */
00322 template <unsigned dim>
00323 class ConvolutionOptions
00324 {
00325   public:
00326     typedef typename MultiArrayShape<dim>::type Shape;
00327     typedef detail::multiArrayScaleParam<dim> ParamVec;
00328     typedef typename ParamVec::return_type    ParamIt;
00329 
00330     ParamVec sigma_eff;
00331     ParamVec sigma_d;
00332     ParamVec step_size;
00333     ParamVec outer_scale;
00334     double window_ratio;
00335     Shape from_point, to_point;
00336      
00337     ConvolutionOptions()
00338     : sigma_eff(0.0),
00339       sigma_d(0.0),
00340       step_size(1.0),
00341       outer_scale(0.0),
00342       window_ratio(0.0)
00343     {}
00344 
00345     typedef typename detail::WrapDoubleIteratorTriple<ParamIt, ParamIt, ParamIt>
00346         ScaleIterator;
00347     typedef typename detail::WrapDoubleIterator<ParamIt>::type
00348         StepIterator;
00349 
00350     ScaleIterator scaleParams() const
00351     {
00352         return ScaleIterator(sigma_eff(), sigma_d(), step_size());
00353     }
00354     StepIterator stepParams() const
00355     {
00356         return StepIterator(step_size());
00357     }
00358 
00359     ConvolutionOptions outerOptions() const
00360     {
00361         ConvolutionOptions outer = *this;
00362         // backward-compatible values:
00363         return outer.stdDev(outer_scale()).resolutionStdDev(0.0);
00364     }
00365 
00366     // Step size per axis.
00367     // Default: dim values of 1.0
00368     VIGRA_CONVOLUTION_OPTIONS(stepSize, 1.0, step_size)
00369 #ifdef DOXYGEN
00370         /** Step size(s) per axis, i.e., the distance between two
00371             adjacent pixels. Required for <tt>MultiArray</tt>
00372             containing anisotropic data.
00373  
00374             Note that a convolution containing a derivative operator
00375             of order <tt>n</tt> results in a multiplication by 
00376             \f${\rm stepSize}^{-n}\f$ for each axis.
00377             Also, the above standard deviations
00378             are scaled according to the step size of each axis.
00379             Default value for the options object if this member function is not
00380             used: A value of 1.0 for each dimension.
00381         */
00382     ConvolutionOptions<dim> & stepSize(...);
00383 #endif
00384 
00385     // Resolution standard deviation per axis.
00386     // Default: dim values of 0.0
00387     VIGRA_CONVOLUTION_OPTIONS(resolutionStdDev, 0.0, sigma_d)
00388 #ifdef DOXYGEN
00389         /** Resolution standard deviation(s) per axis, i.e., a supposed
00390             pre-existing gaussian filtering by this value.
00391        
00392             The standard deviation actually used by the convolution operators
00393             is \f$\sqrt{{\rm sigma}^{2} - {\rm resolutionStdDev}^{2}}\f$ for each
00394             axis.
00395             Default value for the options object if this member function is not
00396             used: A value of 0.0 for each dimension.
00397         */
00398     ConvolutionOptions<dim> & resolutionStdDev(...);
00399 #endif
00400 
00401     // Standard deviation of scale space operators.
00402     // Default: dim values of 0.0
00403     VIGRA_CONVOLUTION_OPTIONS(stdDev, 0.0, sigma_eff)
00404     VIGRA_CONVOLUTION_OPTIONS(innerScale, 0.0, sigma_eff)
00405 
00406 #ifdef DOXYGEN
00407         /** Standard deviation(s) of scale space operators, or inner scale(s) for \ref structureTensorMultiArray().
00408         
00409             Usually not
00410             needed, since a single value for all axes may be specified as a parameter
00411             <tt>sigma</tt> to the call of
00412             an convolution operator such as \ref gaussianGradientMultiArray(), and
00413             anisotropic data requiring the use of the stepSize() member function.
00414             Default value for the options object if this member function is not
00415             used: A value of 0.0 for each dimension.
00416         */
00417     ConvolutionOptions<dim> & stdDev(...);
00418 
00419         /** Standard deviation(s) of scale space operators, or inner scale(s) for \ref structureTensorMultiArray().
00420         
00421             Usually not
00422             needed, since a single value for all axes may be specified as a parameter
00423             <tt>sigma</tt> to the call of
00424             an convolution operator such as \ref gaussianGradientMultiArray(), and
00425             anisotropic data requiring the use of the stepSize() member function.
00426             Default value for the options object if this member function is not
00427             used: A value of 0.0 for each dimension.
00428         */
00429     ConvolutionOptions<dim> & innerScale(...);
00430 #endif
00431 
00432     // Outer scale, for structure tensor.
00433     // Default: dim values of 0.0
00434     VIGRA_CONVOLUTION_OPTIONS(outerScale, 0.0, outer_scale)
00435 #ifdef DOXYGEN
00436         /** Standard deviation(s) of the second convolution of the
00437             structure tensor. 
00438 
00439             Usually not needed, since a single value for
00440             all axes may be specified as a parameter <tt>outerScale</tt> to
00441             the call of \ref structureTensorMultiArray(), and
00442             anisotropic data requiring the use of the stepSize() member
00443             function.
00444             Default value for the options object if this member function is not
00445             used: A value of 0.0 for each dimension.
00446         */
00447     ConvolutionOptions<dim> & outerScale(...);
00448 #endif
00449 
00450         /** Size of the filter window as a multiple of the scale parameter. 
00451 
00452             This option is only used for Gaussian filters and their derivatives.
00453             By default, the window size of a Gaussian filter is automatically 
00454             determined such that the error resulting from restricting the 
00455             infinitely large Gaussian function to a finite size is minimized. 
00456             In particular, the window radius is determined as
00457             <tt>radius = round(3.0 * sigma + 0.5 * order)</tt>, where 'order' is the 
00458             desired derivative order. In some cases, it is desirable to trade off 
00459             accuracy for speed, and this function can be used to request a smaller
00460             window radius.
00461             
00462             Default: <tt>0.0</tt> (i.e. determine the window size automatically)
00463         */
00464     ConvolutionOptions<dim> & filterWindowSize(double ratio)
00465     {
00466         vigra_precondition(ratio >= 0.0,
00467             "ConvolutionOptions::filterWindowSize(): ratio must not be negative.");
00468         window_ratio = ratio;
00469         return *this;
00470     }
00471 
00472         /** Restrict the filter to a subregion of the input array. 
00473 
00474             This is useful for speeding up computations by ignoring irrelevant 
00475             areas in the array. <b>Note:</b> It is assumed that the output array
00476             of the convolution has the size given in this function.
00477             
00478             Default: <tt>from = Shape(), to = Shape()</tt> (i.e. use entire array)
00479         */
00480     ConvolutionOptions<dim> & subarray(Shape const & from, Shape const & to)
00481     {
00482         from_point = from;
00483         to_point = to;
00484         return *this;
00485     }
00486 };
00487 
00488 namespace detail
00489 {
00490 
00491 /********************************************************/
00492 /*                                                      */
00493 /*        internalSeparableConvolveMultiArray           */
00494 /*                                                      */
00495 /********************************************************/
00496 
00497 template <class SrcIterator, class SrcShape, class SrcAccessor,
00498           class DestIterator, class DestAccessor, class KernelIterator>
00499 void
00500 internalSeparableConvolveMultiArrayTmp(
00501                       SrcIterator si, SrcShape const & shape, SrcAccessor src,
00502                       DestIterator di, DestAccessor dest, KernelIterator kit)
00503 {
00504     enum { N = 1 + SrcIterator::level };
00505 
00506     typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
00507     typedef typename AccessorTraits<TmpType>::default_accessor TmpAcessor;
00508 
00509     // temporary array to hold the current line to enable in-place operation
00510     ArrayVector<TmpType> tmp( shape[0] );
00511 
00512     typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
00513     typedef MultiArrayNavigator<DestIterator, N> DNavigator;
00514     
00515     TmpAcessor acc;
00516 
00517     {
00518         // only operate on first dimension here
00519         SNavigator snav( si, shape, 0 );
00520         DNavigator dnav( di, shape, 0 );
00521         
00522         for( ; snav.hasMore(); snav++, dnav++ )
00523         {
00524              // first copy source to tmp for maximum cache efficiency
00525              copyLine(snav.begin(), snav.end(), src, tmp.begin(), acc);
00526 
00527              convolveLine(srcIterRange(tmp.begin(), tmp.end(), acc),
00528                           destIter( dnav.begin(), dest ),
00529                           kernel1d( *kit ) );
00530         }
00531         ++kit;
00532     }
00533 
00534     // operate on further dimensions
00535     for( int d = 1; d < N; ++d, ++kit )
00536     {
00537         DNavigator dnav( di, shape, d );
00538 
00539         tmp.resize( shape[d] );
00540 
00541         for( ; dnav.hasMore(); dnav++ )
00542         {
00543              // first copy source to tmp since convolveLine() cannot work in-place
00544              copyLine(dnav.begin(), dnav.end(), dest, tmp.begin(), acc);
00545 
00546              convolveLine(srcIterRange(tmp.begin(), tmp.end(), acc),
00547                           destIter( dnav.begin(), dest ),
00548                           kernel1d( *kit ) );
00549         }
00550     }
00551 }
00552 
00553 /********************************************************/
00554 /*                                                      */
00555 /*         internalSeparableConvolveSubarray            */
00556 /*                                                      */
00557 /********************************************************/
00558 
00559 template <class SrcIterator, class SrcShape, class SrcAccessor,
00560           class DestIterator, class DestAccessor, class KernelIterator>
00561 void
00562 internalSeparableConvolveSubarray(
00563                       SrcIterator si, SrcShape const & shape, SrcAccessor src,
00564                       DestIterator di, DestAccessor dest, KernelIterator kit,
00565                       SrcShape const & start, SrcShape const & stop)
00566 {
00567     enum { N = 1 + SrcIterator::level };
00568 
00569     typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
00570     typedef MultiArray<N, TmpType> TmpArray;
00571     typedef typename TmpArray::traverser TmpIterator;
00572     typedef typename AccessorTraits<TmpType>::default_accessor TmpAcessor;
00573     
00574     SrcShape sstart, sstop, axisorder, tmpshape;
00575     TinyVector<double, N> overhead;
00576     for(int k=0; k<N; ++k)
00577     {
00578         sstart[k] = start[k] - kit[k].right();
00579         if(sstart[k] < 0)
00580             sstart[k] = 0;
00581         sstop[k] = stop[k] - kit[k].left();
00582         if(sstop[k] > shape[k])
00583             sstop[k] = shape[k];
00584         overhead[k] = double(sstop[k] - sstart[k]) / (stop[k] - start[k]);
00585     }
00586     
00587     indexSort(overhead.begin(), overhead.end(), axisorder.begin(), std::greater<double>());
00588     
00589     SrcShape dstart, dstop(sstop - sstart);
00590     dstop[axisorder[0]]  = stop[axisorder[0]] - start[axisorder[0]];
00591     
00592     // temporary array to hold the current line to enable in-place operation
00593     MultiArray<N, TmpType> tmp(dstop);
00594 
00595     typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
00596     typedef MultiArrayNavigator<TmpIterator, N> TNavigator;
00597     typedef MultiArrayNavigator<DestIterator, N> DNavigator;
00598     
00599     TmpAcessor acc;
00600 
00601     {
00602         // only operate on first dimension here
00603         SNavigator snav( si, sstart, sstop, axisorder[0]);
00604         TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[0]);
00605         
00606         ArrayVector<TmpType> tmpline(sstop[axisorder[0]] - sstart[axisorder[0]]);
00607         
00608         int lstart = start[axisorder[0]] - sstart[axisorder[0]];
00609         int lstop  = lstart + (stop[axisorder[0]] - start[axisorder[0]]);
00610 
00611         for( ; snav.hasMore(); snav++, tnav++ )
00612         {
00613             // first copy source to tmp for maximum cache efficiency
00614             copyLine(snav.begin(), snav.end(), src, tmpline.begin(), acc);
00615             
00616             convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
00617                          destIter(tnav.begin(), acc),
00618                          kernel1d( kit[axisorder[0]] ), lstart, lstop);
00619         }
00620     }
00621     
00622     // operate on further dimensions
00623     for( int d = 1; d < N; ++d)
00624     {
00625         TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[d]);
00626         
00627         ArrayVector<TmpType> tmpline(dstop[axisorder[d]] - dstart[axisorder[d]]);
00628         
00629         int lstart = start[axisorder[d]] - sstart[axisorder[d]];
00630         int lstop  = lstart + (stop[axisorder[d]] - start[axisorder[d]]);
00631 
00632         for( ; tnav.hasMore(); tnav++ )
00633         {
00634             // first copy source to tmp because convolveLine() cannot work in-place
00635             copyLine(tnav.begin(), tnav.end(), acc, tmpline.begin(), acc );
00636 
00637             convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
00638                          destIter( tnav.begin() + lstart, acc ),
00639                          kernel1d( kit[axisorder[d]] ), lstart, lstop);
00640         }
00641         
00642         dstart[axisorder[d]] = lstart;
00643         dstop[axisorder[d]] = lstop;
00644     }
00645     
00646     copyMultiArray(tmp.traverser_begin()+dstart, stop-start, acc, di, dest);              
00647 }
00648 
00649 
00650 template <class K>
00651 void 
00652 scaleKernel(K & kernel, double a)
00653 {
00654     for(int i = kernel.left(); i <= kernel.right(); ++i)
00655         kernel[i] = detail::RequiresExplicitCast<typename K::value_type>::cast(kernel[i] * a);
00656 }
00657 
00658 
00659 } // namespace detail
00660 
00661 /** \addtogroup MultiArrayConvolutionFilters Convolution filters for multi-dimensional arrays.
00662 
00663     These functions realize a separable convolution on an arbitrary dimensional
00664     array that is specified by iterators (compatible to \ref MultiIteratorPage)
00665     and shape objects. It can therefore be applied to a wide range of data structures
00666     (\ref vigra::MultiArrayView, \ref vigra::MultiArray etc.).
00667 */
00668 //@{
00669 
00670 /********************************************************/
00671 /*                                                      */
00672 /*             separableConvolveMultiArray              */
00673 /*                                                      */
00674 /********************************************************/
00675 
00676 /** \brief Separated convolution on multi-dimensional arrays.
00677 
00678     This function computes a separated convolution on all dimensions
00679     of the given multi-dimensional array. Both source and destination
00680     arrays are represented by iterators, shape objects and accessors.
00681     The destination array is required to already have the correct size.
00682 
00683     There are two variants of this functions: one takes a single kernel
00684     of type \ref vigra::Kernel1D which is then applied to all dimensions,
00685     whereas the other requires an iterator referencing a sequence of
00686     \ref vigra::Kernel1D objects, one for every dimension of the data.
00687     Then the first kernel in this sequence is applied to the innermost
00688     dimension (e.g. the x-dimension of an image), while the last is applied to the
00689     outermost dimension (e.g. the z-dimension in a 3D image).
00690 
00691     This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
00692     A full-sized internal array is only allocated if working on the destination
00693     array directly would cause round-off errors (i.e. if
00694     <tt>typeid(typename NumericTraits<typename DestAccessor::value_type>::RealPromote)
00695     != typeid(typename DestAccessor::value_type)</tt>.
00696     
00697     If <tt>start</tt> and <tt>stop</tt> have non-default values, they must represent
00698     a valid subarray of the input array. The convolution is then restricted to that 
00699     subarray, and it is assumed that the output array only refers to the
00700     subarray (i.e. <tt>diter</tt> points to the element corresponding to 
00701     <tt>start</tt>). 
00702 
00703     <b> Declarations:</b>
00704 
00705     pass arguments explicitly:
00706     \code
00707     namespace vigra {
00708         // apply the same kernel to all dimensions
00709         template <class SrcIterator, class SrcShape, class SrcAccessor,
00710                   class DestIterator, class DestAccessor, class T>
00711         void
00712         separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00713                                     DestIterator diter, DestAccessor dest,
00714                                     Kernel1D<T> const & kernel,
00715                                     SrcShape const & start = SrcShape(),
00716                                     SrcShape const & stop = SrcShape());
00717 
00718         // apply each kernel from the sequence 'kernels' in turn
00719         template <class SrcIterator, class SrcShape, class SrcAccessor,
00720                   class DestIterator, class DestAccessor, class KernelIterator>
00721         void
00722         separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00723                                     DestIterator diter, DestAccessor dest,
00724                                     KernelIterator kernels,
00725                                     SrcShape const & start = SrcShape(),
00726                                     SrcShape const & stop = SrcShape());
00727     }
00728     \endcode
00729 
00730     use argument objects in conjunction with \ref ArgumentObjectFactories :
00731     \code
00732     namespace vigra {
00733         // apply the same kernel to all dimensions
00734         template <class SrcIterator, class SrcShape, class SrcAccessor,
00735                   class DestIterator, class DestAccessor, class T>
00736         void
00737         separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00738                                     pair<DestIterator, DestAccessor> const & dest,
00739                                     Kernel1D<T> const & kernel,
00740                                     SrcShape const & start = SrcShape(),
00741                                     SrcShape const & stop = SrcShape());
00742 
00743         // apply each kernel from the sequence 'kernels' in turn
00744         template <class SrcIterator, class SrcShape, class SrcAccessor,
00745                   class DestIterator, class DestAccessor, class KernelIterator>
00746         void
00747         separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00748                                     pair<DestIterator, DestAccessor> const & dest,
00749                                     KernelIterator kernels,
00750                                     SrcShape const & start = SrcShape(),
00751                                     SrcShape const & stop = SrcShape());
00752     }
00753     \endcode
00754 
00755     <b> Usage:</b>
00756 
00757     <b>\#include</b> <vigra/multi_convolution.hxx>
00758 
00759     \code
00760     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
00761     MultiArray<3, unsigned char> source(shape);
00762     MultiArray<3, float> dest(shape);
00763     ...
00764     Kernel1D<float> gauss;
00765     gauss.initGaussian(sigma);
00766     // create 3 Gauss kernels, one for each dimension
00767     ArrayVector<Kernel1D<float> > kernels(3, gauss);
00768 
00769     // perform Gaussian smoothing on all dimensions
00770     separableConvolveMultiArray(srcMultiArrayRange(source), destMultiArray(dest), 
00771                                 kernels.begin());
00772     \endcode
00773 
00774     <b> Required Interface:</b>
00775 
00776     see \ref separableConvolveMultiArray(), in addition:
00777 
00778     \code
00779     int dimension = 0;
00780     VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
00781     \endcode
00782 
00783     \see vigra::Kernel1D, convolveLine()
00784 */
00785 doxygen_overloaded_function(template <...> void separableConvolveMultiArray)
00786 
00787 template <class SrcIterator, class SrcShape, class SrcAccessor,
00788           class DestIterator, class DestAccessor, class KernelIterator>
00789 void
00790 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
00791                              DestIterator d, DestAccessor dest, 
00792                              KernelIterator kernels,
00793                              SrcShape const & start = SrcShape(),
00794                              SrcShape const & stop = SrcShape())
00795 {
00796     typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
00797 
00798     if(stop != SrcShape())
00799     {
00800         enum { N = 1 + SrcIterator::level };
00801         for(int k=0; k<N; ++k)
00802             vigra_precondition(0 <= start[k] && start[k] < stop[k] && stop[k] <= shape[k],
00803               "separableConvolveMultiArray(): invalid subarray shape.");
00804 
00805         detail::internalSeparableConvolveSubarray(s, shape, src, d, dest, kernels, start, stop);
00806     }
00807     else if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
00808     {
00809         // need a temporary array to avoid rounding errors
00810         MultiArray<SrcShape::static_size, TmpType> tmpArray(shape);
00811         detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
00812              tmpArray.traverser_begin(), typename AccessorTraits<TmpType>::default_accessor(), kernels );
00813         copyMultiArray(srcMultiArrayRange(tmpArray), destIter(d, dest));
00814     }
00815     else
00816     {
00817         // work directly on the destination array
00818         detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
00819     }
00820 }
00821 
00822 template <class SrcIterator, class SrcShape, class SrcAccessor,
00823           class DestIterator, class DestAccessor, class KernelIterator>
00824 inline
00825 void separableConvolveMultiArray(
00826     triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00827     pair<DestIterator, DestAccessor> const & dest, 
00828     KernelIterator kit,
00829     SrcShape const & start = SrcShape(),
00830     SrcShape const & stop = SrcShape())
00831 {
00832     separableConvolveMultiArray( source.first, source.second, source.third,
00833                                  dest.first, dest.second, kit, start, stop );
00834 }
00835 
00836 template <class SrcIterator, class SrcShape, class SrcAccessor,
00837           class DestIterator, class DestAccessor, class T>
00838 inline void
00839 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
00840                              DestIterator d, DestAccessor dest,
00841                              Kernel1D<T> const & kernel,
00842                              SrcShape const & start = SrcShape(),
00843                              SrcShape const & stop = SrcShape())
00844 {
00845     ArrayVector<Kernel1D<T> > kernels(shape.size(), kernel);
00846 
00847     separableConvolveMultiArray( s, shape, src, d, dest, kernels.begin(), start, stop);
00848 }
00849 
00850 template <class SrcIterator, class SrcShape, class SrcAccessor,
00851           class DestIterator, class DestAccessor, class T>
00852 inline void
00853 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00854                             pair<DestIterator, DestAccessor> const & dest,
00855                             Kernel1D<T> const & kernel,
00856                             SrcShape const & start = SrcShape(),
00857                             SrcShape const & stop = SrcShape())
00858 {
00859     ArrayVector<Kernel1D<T> > kernels(source.second.size(), kernel);
00860 
00861     separableConvolveMultiArray( source.first, source.second, source.third,
00862                                  dest.first, dest.second, kernels.begin(), start, stop);
00863 }
00864 
00865 /********************************************************/
00866 /*                                                      */
00867 /*            convolveMultiArrayOneDimension            */
00868 /*                                                      */
00869 /********************************************************/
00870 
00871 /** \brief Convolution along a single dimension of a multi-dimensional arrays.
00872 
00873     This function computes a convolution along one dimension (specified by
00874     the parameter <tt>dim</tt> of the given multi-dimensional array with the given
00875     <tt>kernel</tt>. Both source and destination arrays are represented by
00876     iterators, shape objects and accessors. The destination array is required to
00877     already have the correct size.
00878 
00879     If <tt>start</tt> and <tt>stop</tt> have non-default values, they must represent
00880     a valid subarray of the input array. The convolution is then restricted to that 
00881     subarray, and it is assumed that the output array only refers to the
00882     subarray (i.e. <tt>diter</tt> points to the element corresponding to 
00883     <tt>start</tt>). 
00884 
00885     This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
00886 
00887     <b> Declarations:</b>
00888 
00889     pass arguments explicitly:
00890     \code
00891     namespace vigra {
00892         template <class SrcIterator, class SrcShape, class SrcAccessor,
00893                   class DestIterator, class DestAccessor, class T>
00894         void
00895         convolveMultiArrayOneDimension(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00896                                        DestIterator diter, DestAccessor dest,
00897                                        unsigned int dim, vigra::Kernel1D<T> const & kernel,
00898                                        SrcShape const & start = SrcShape(),
00899                                        SrcShape const & stop = SrcShape());
00900     }
00901     \endcode
00902 
00903     use argument objects in conjunction with \ref ArgumentObjectFactories :
00904     \code
00905     namespace vigra {
00906         template <class SrcIterator, class SrcShape, class SrcAccessor,
00907                   class DestIterator, class DestAccessor, class T>
00908         void
00909         convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00910                                        pair<DestIterator, DestAccessor> const & dest,
00911                                        unsigned int dim, vigra::Kernel1D<T> const & kernel,
00912                                        SrcShape const & start = SrcShape(),
00913                                        SrcShape const & stop = SrcShape());
00914     }
00915     \endcode
00916 
00917     <b> Usage:</b>
00918 
00919     <b>\#include</b> <vigra/multi_convolution.hxx>
00920 
00921     \code
00922     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
00923     MultiArray<3, unsigned char> source(shape);
00924     MultiArray<3, float> dest(shape);
00925     ...
00926     Kernel1D<float> gauss;
00927     gauss.initGaussian(sigma);
00928 
00929     // perform Gaussian smoothing along dimensions 1 (height)
00930     convolveMultiArrayOneDimension(srcMultiArrayRange(source), destMultiArray(dest), 1, gauss);
00931     \endcode
00932 
00933     \see separableConvolveMultiArray()
00934 */
00935 doxygen_overloaded_function(template <...> void convolveMultiArrayOneDimension)
00936 
00937 template <class SrcIterator, class SrcShape, class SrcAccessor,
00938           class DestIterator, class DestAccessor, class T>
00939 void
00940 convolveMultiArrayOneDimension(SrcIterator s, SrcShape const & shape, SrcAccessor src,
00941                                DestIterator d, DestAccessor dest,
00942                                unsigned int dim, vigra::Kernel1D<T> const & kernel,
00943                                SrcShape const & start = SrcShape(),
00944                                SrcShape const & stop = SrcShape())
00945 {
00946     enum { N = 1 + SrcIterator::level };
00947     vigra_precondition( dim < N,
00948                         "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller "
00949                         "than the data dimensionality" );
00950 
00951     typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
00952     typedef typename AccessorTraits<TmpType>::default_const_accessor TmpAccessor;
00953     ArrayVector<TmpType> tmp( shape[dim] );
00954 
00955     typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
00956     typedef MultiArrayNavigator<DestIterator, N> DNavigator;
00957     
00958     SrcShape sstart, sstop(shape), dstart, dstop(shape);
00959     
00960     if(stop != SrcShape())
00961     {
00962         sstart = start;
00963         sstop  = stop;
00964         sstart[dim] = 0;
00965         sstop[dim]  = shape[dim];
00966         dstop = stop - start;
00967     }
00968 
00969     SNavigator snav( s, sstart, sstop, dim );
00970     DNavigator dnav( d, dstart, dstop, dim );
00971 
00972     for( ; snav.hasMore(); snav++, dnav++ )
00973     {
00974          // first copy source to temp for maximum cache efficiency
00975          copyLine( snav.begin(), snav.end(), src,
00976            tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
00977 
00978          convolveLine(srcIterRange( tmp.begin(), tmp.end(), TmpAccessor()),
00979                       destIter( dnav.begin(), dest ),
00980                       kernel1d( kernel), start[dim], stop[dim]);
00981     }
00982 }
00983 
00984 template <class SrcIterator, class SrcShape, class SrcAccessor,
00985           class DestIterator, class DestAccessor, class T>
00986 inline void
00987 convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00988                                pair<DestIterator, DestAccessor> const & dest,
00989                                unsigned int dim, vigra::Kernel1D<T> const & kernel,
00990                                SrcShape const & start = SrcShape(),
00991                                SrcShape const & stop = SrcShape())
00992 {
00993     convolveMultiArrayOneDimension(source.first, source.second, source.third,
00994                                    dest.first, dest.second, dim, kernel, start, stop);
00995 }
00996 
00997 /********************************************************/
00998 /*                                                      */
00999 /*             gaussianSmoothMultiArray                 */
01000 /*                                                      */
01001 /********************************************************/
01002 
01003 /** \brief Isotropic Gaussian smoothing of a multi-dimensional arrays.
01004 
01005     This function computes an isotropic convolution of the given N-dimensional
01006     array with a Gaussian filter at the given standard deviation <tt>sigma</tt>.
01007     Both source and destination arrays are represented by
01008     iterators, shape objects and accessors. The destination array is required to
01009     already have the correct size. This function may work in-place, which means
01010     that <tt>siter == diter</tt> is allowed. It is implemented by a call to
01011     \ref separableConvolveMultiArray() with the appropriate kernel.
01012 
01013     Anisotropic data should be passed with appropriate
01014     \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional
01015     unless the parameter <tt>sigma</tt> is left out.
01016 
01017     <b> Declarations:</b>
01018 
01019     pass arguments explicitly:
01020     \code
01021     namespace vigra {
01022         template <class SrcIterator, class SrcShape, class SrcAccessor,
01023                   class DestIterator, class DestAccessor>
01024         void
01025         gaussianSmoothMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
01026                                  DestIterator diter, DestAccessor dest,
01027                                  double sigma, const ConvolutionOptions<N> & opt);
01028     }
01029     \endcode
01030 
01031     use argument objects in conjunction with \ref ArgumentObjectFactories :
01032     \code
01033     namespace vigra {
01034         template <class SrcIterator, class SrcShape, class SrcAccessor,
01035                   class DestIterator, class DestAccessor>
01036         void
01037         gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01038                                  pair<DestIterator, DestAccessor> const & dest,
01039                                  double sigma, const ConvolutionOptions<N> & opt);
01040     }
01041     \endcode
01042 
01043     <b> Usage:</b>
01044 
01045     <b>\#include</b> <vigra/multi_convolution.hxx>
01046 
01047     \code
01048     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
01049     MultiArray<3, unsigned char> source(shape);
01050     MultiArray<3, float> dest(shape);
01051     ...
01052     // perform isotropic Gaussian smoothing at scale 'sigma'
01053     gaussianSmoothMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
01054     \endcode
01055 
01056     <b> Usage with anisotropic data:</b>
01057 
01058     <b>\#include</b> <vigra/multi_convolution.hxx>
01059 
01060     \code
01061     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
01062     MultiArray<3, unsigned char> source(shape);
01063     MultiArray<3, float> dest(shape);
01064     TinyVector<float, 3> step_size;
01065     TinyVector<float, 3> resolution_sigmas;
01066     ...
01067     // perform anisotropic Gaussian smoothing at scale 'sigma'
01068     gaussianSmoothMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma,
01069                              ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
01070     \endcode
01071 
01072     \see separableConvolveMultiArray()
01073 */
01074 doxygen_overloaded_function(template <...> void gaussianSmoothMultiArray)
01075 
01076 template <class SrcIterator, class SrcShape, class SrcAccessor,
01077           class DestIterator, class DestAccessor>
01078 void
01079 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
01080                    DestIterator d, DestAccessor dest,
01081                    const ConvolutionOptions<SrcShape::static_size> & opt,
01082                    const char *const function_name = "gaussianSmoothMultiArray" )
01083 {
01084     typedef typename DestAccessor::value_type DestType;
01085 
01086     static const int N = SrcShape::static_size;
01087 
01088     typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams();
01089     ArrayVector<Kernel1D<double> > kernels(N);
01090 
01091     for (int dim = 0; dim < N; ++dim, ++params)
01092         kernels[dim].initGaussian(params.sigma_scaled(function_name), 1.0, opt.window_ratio);
01093 
01094     separableConvolveMultiArray(s, shape, src, d, dest, kernels.begin(), opt.from_point, opt.to_point);
01095 }
01096 
01097 template <class SrcIterator, class SrcShape, class SrcAccessor,
01098           class DestIterator, class DestAccessor>
01099 inline void
01100 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
01101                    DestIterator d, DestAccessor dest, double sigma,
01102                    const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01103 {
01104     ConvolutionOptions<SrcShape::static_size> par = opt;
01105     gaussianSmoothMultiArray(s, shape, src, d,  dest, par.stdDev(sigma));
01106 }
01107 
01108 template <class SrcIterator, class SrcShape, class SrcAccessor,
01109           class DestIterator, class DestAccessor>
01110 inline void
01111 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01112                   pair<DestIterator, DestAccessor> const & dest,
01113                   const ConvolutionOptions<SrcShape::static_size> & opt)
01114 {
01115     gaussianSmoothMultiArray( source.first, source.second, source.third,
01116                               dest.first, dest.second, opt );
01117 }
01118 
01119 template <class SrcIterator, class SrcShape, class SrcAccessor,
01120           class DestIterator, class DestAccessor>
01121 inline void
01122 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01123                   pair<DestIterator, DestAccessor> const & dest, double sigma,
01124                   const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01125 {
01126     gaussianSmoothMultiArray( source.first, source.second, source.third,
01127                               dest.first, dest.second, sigma, opt );
01128 }
01129 
01130 
01131 /********************************************************/
01132 /*                                                      */
01133 /*             gaussianGradientMultiArray               */
01134 /*                                                      */
01135 /********************************************************/
01136 
01137 /** \brief Calculate Gaussian gradient of a multi-dimensional arrays.
01138 
01139     This function computes the Gaussian gradient of the given N-dimensional
01140     array with a sequence of first-derivative-of-Gaussian filters at the given
01141     standard deviation <tt>sigma</tt> (differentiation is applied to each dimension
01142     in turn, starting with the innermost dimension). Both source and destination arrays
01143     are represented by iterators, shape objects and accessors. The destination array is
01144     required to have a vector valued pixel type with as many elements as the number of
01145     dimensions. This function is implemented by calls to
01146     \ref separableConvolveMultiArray() with the appropriate kernels.
01147 
01148     Anisotropic data should be passed with appropriate
01149     \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional
01150     unless the parameter <tt>sigma</tt> is left out.
01151 
01152     <b> Declarations:</b>
01153 
01154     pass arguments explicitly:
01155     \code
01156     namespace vigra {
01157         template <class SrcIterator, class SrcShape, class SrcAccessor,
01158                   class DestIterator, class DestAccessor>
01159         void
01160         gaussianGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
01161                                    DestIterator diter, DestAccessor dest,
01162                                    double sigma, const ConvolutionOptions<N> & opt);
01163     }
01164     \endcode
01165 
01166     use argument objects in conjunction with \ref ArgumentObjectFactories :
01167     \code
01168     namespace vigra {
01169         template <class SrcIterator, class SrcShape, class SrcAccessor,
01170                   class DestIterator, class DestAccessor>
01171         void
01172         gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01173                                    pair<DestIterator, DestAccessor> const & dest,
01174                                    double sigma, const ConvolutionOptions<N> & opt);
01175     }
01176     \endcode
01177 
01178     <b> Usage:</b>
01179 
01180     <b>\#include</b> <vigra/multi_convolution.hxx>
01181 
01182     \code
01183     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
01184     MultiArray<3, unsigned char> source(shape);
01185     MultiArray<3, TinyVector<float, 3> > dest(shape);
01186     ...
01187     // compute Gaussian gradient at scale sigma
01188     gaussianGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
01189     \endcode
01190 
01191     <b> Usage with anisotropic data:</b>
01192 
01193     <b>\#include</b> <vigra/multi_convolution.hxx>
01194 
01195     \code
01196     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
01197     MultiArray<3, unsigned char> source(shape);
01198     MultiArray<3, TinyVector<float, 3> > dest(shape);
01199     TinyVector<float, 3> step_size;
01200     TinyVector<float, 3> resolution_sigmas;
01201     ...
01202     // compute Gaussian gradient at scale sigma
01203     gaussianGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma,
01204                                ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
01205     \endcode
01206 
01207     <b> Required Interface:</b>
01208 
01209     see \ref separableConvolveMultiArray(), in addition:
01210 
01211     \code
01212     int dimension = 0;
01213     VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
01214     \endcode
01215 
01216     \see separableConvolveMultiArray()
01217 */
01218 doxygen_overloaded_function(template <...> void gaussianGradientMultiArray)
01219 
01220 template <class SrcIterator, class SrcShape, class SrcAccessor,
01221           class DestIterator, class DestAccessor>
01222 void
01223 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
01224                            DestIterator di, DestAccessor dest,
01225                            ConvolutionOptions<SrcShape::static_size> const & opt,
01226                            const char *const function_name = "gaussianGradientMultiArray")
01227 {
01228     typedef typename DestAccessor::value_type DestType;
01229     typedef typename DestType::value_type     DestValueType;
01230     typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
01231    
01232     static const int N = SrcShape::static_size;
01233     typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
01234 
01235     for(int k=0; k<N; ++k)
01236         if(shape[k] <=0)
01237             return;
01238 
01239     vigra_precondition(N == (int)dest.size(di),
01240         "gaussianGradientMultiArray(): Wrong number of channels in output array.");
01241 
01242     ParamType params = opt.scaleParams();
01243     ParamType params2(params);
01244 
01245     ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
01246     for (int dim = 0; dim < N; ++dim, ++params)
01247     {
01248         double sigma = params.sigma_scaled(function_name);
01249         plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
01250     }
01251 
01252     typedef VectorElementAccessor<DestAccessor> ElementAccessor;
01253 
01254     // compute gradient components
01255     for (int dim = 0; dim < N; ++dim, ++params2)
01256     {
01257         ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
01258         kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 1, 1.0, opt.window_ratio);
01259         detail::scaleKernel(kernels[dim], 1.0 / params2.step_size());
01260         separableConvolveMultiArray(si, shape, src, di, ElementAccessor(dim, dest), kernels.begin(), 
01261                                     opt.from_point, opt.to_point);
01262     }
01263 }
01264 
01265 template <class SrcIterator, class SrcShape, class SrcAccessor,
01266           class DestIterator, class DestAccessor>
01267 void
01268 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
01269                            DestIterator di, DestAccessor dest, double sigma,
01270                            const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01271 {
01272     ConvolutionOptions<SrcShape::static_size> par = opt;
01273     gaussianGradientMultiArray(si, shape, src, di, dest, par.stdDev(sigma));
01274 }
01275 
01276 template <class SrcIterator, class SrcShape, class SrcAccessor,
01277           class DestIterator, class DestAccessor>
01278 inline void
01279 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01280                            pair<DestIterator, DestAccessor> const & dest,
01281                            ConvolutionOptions<SrcShape::static_size> const & opt )
01282 {
01283     gaussianGradientMultiArray( source.first, source.second, source.third,
01284                                 dest.first, dest.second, opt );
01285 }
01286 
01287 template <class SrcIterator, class SrcShape, class SrcAccessor,
01288           class DestIterator, class DestAccessor>
01289 inline void
01290 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01291                            pair<DestIterator, DestAccessor> const & dest,
01292                            double sigma,
01293                            const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01294 {
01295     gaussianGradientMultiArray( source.first, source.second, source.third,
01296                                 dest.first, dest.second, sigma, opt );
01297 }
01298 
01299 /********************************************************/
01300 /*                                                      */
01301 /*             symmetricGradientMultiArray              */
01302 /*                                                      */
01303 /********************************************************/
01304 
01305 /** \brief Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
01306 
01307     This function computes the gradient of the given N-dimensional
01308     array with a sequence of symmetric difference filters a (differentiation is applied
01309     to each dimension in turn, starting with the innermost dimension). Both source and
01310     destination arrays are represented by iterators, shape objects and accessors.
01311     The destination array is required to have a vector valued pixel type with as many
01312     elements as the number of dimensions. This function is implemented by calls to
01313     \ref convolveMultiArrayOneDimension() with the symmetric difference kernel.
01314 
01315     Anisotropic data should be passed with appropriate
01316     \ref ConvolutionOptions, the parameter <tt>opt</tt> is optional
01317     otherwise.
01318     
01319     <b> Declarations:</b>
01320 
01321     pass arguments explicitly:
01322     \code
01323     namespace vigra {
01324         template <class SrcIterator, class SrcShape, class SrcAccessor,
01325                   class DestIterator, class DestAccessor>
01326         void
01327         symmetricGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
01328                                     DestIterator diter, DestAccessor dest,
01329                                     const ConvolutionOptions<N> & opt);
01330     }
01331     \endcode
01332 
01333     use argument objects in conjunction with \ref ArgumentObjectFactories :
01334     \code
01335     namespace vigra {
01336         template <class SrcIterator, class SrcShape, class SrcAccessor,
01337                   class DestIterator, class DestAccessor>
01338         void
01339         symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01340                                     pair<DestIterator, DestAccessor> const & dest,
01341                                     const ConvolutionOptions<N> & opt);
01342     }
01343     \endcode
01344 
01345     <b> Usage:</b>
01346 
01347     <b>\#include</b> <vigra/multi_convolution.hxx>
01348 
01349     \code
01350     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
01351     MultiArray<3, unsigned char> source(shape);
01352     MultiArray<3, TinyVector<float, 3> > dest(shape);
01353     ...
01354     // compute gradient
01355     symmetricGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest));
01356     \endcode
01357 
01358     <b> Usage with anisotropic data:</b>
01359 
01360     <b>\#include</b> <vigra/multi_convolution.hxx>
01361 
01362     \code
01363     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
01364     MultiArray<3, unsigned char> source(shape);
01365     MultiArray<3, TinyVector<float, 3> > dest(shape);
01366     TinyVector<float, 3> step_size;
01367     ...
01368     // compute gradient
01369     symmetricGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest),
01370                                 ConvolutionOptions<3>().stepSize(step_size));
01371     \endcode
01372 
01373     <b> Required Interface:</b>
01374 
01375     see \ref convolveMultiArrayOneDimension(), in addition:
01376 
01377     \code
01378     int dimension = 0;
01379     VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
01380     \endcode
01381 
01382     \see convolveMultiArrayOneDimension()
01383 */
01384 doxygen_overloaded_function(template <...> void symmetricGradientMultiArray)
01385 
01386 template <class SrcIterator, class SrcShape, class SrcAccessor,
01387           class DestIterator, class DestAccessor>
01388 void
01389 symmetricGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
01390                             DestIterator di, DestAccessor dest,
01391                             const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01392 {
01393     typedef typename DestAccessor::value_type DestType;
01394     typedef typename DestType::value_type     DestValueType;
01395     typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
01396 
01397     static const int N = SrcShape::static_size;
01398     typedef typename ConvolutionOptions<N>::StepIterator StepType;
01399 
01400     for(int k=0; k<N; ++k)
01401         if(shape[k] <=0)
01402             return;
01403 
01404     vigra_precondition(N == (int)dest.size(di),
01405         "symmetricGradientMultiArray(): Wrong number of channels in output array.");
01406 
01407     Kernel1D<KernelType> filter;
01408     filter.initSymmetricGradient();
01409 
01410     StepType step_size_it = opt.stepParams();
01411 
01412     typedef VectorElementAccessor<DestAccessor> ElementAccessor;
01413 
01414     // compute gradient components
01415     for (int d = 0; d < N; ++d, ++step_size_it)
01416     {
01417         Kernel1D<KernelType> symmetric(filter);
01418         detail::scaleKernel(symmetric, 1 / *step_size_it);
01419         convolveMultiArrayOneDimension(si, shape, src,
01420                                        di, ElementAccessor(d, dest),
01421                                        d, symmetric, opt.from_point, opt.to_point);
01422     }
01423 }
01424 
01425 template <class SrcIterator, class SrcShape, class SrcAccessor,
01426           class DestIterator, class DestAccessor>
01427 inline void
01428 symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01429                             pair<DestIterator, DestAccessor> const & dest,
01430                             const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01431 {
01432     symmetricGradientMultiArray(source.first, source.second, source.third,
01433                                 dest.first, dest.second, opt);
01434 }
01435 
01436 /********************************************************/
01437 /*                                                      */
01438 /*            laplacianOfGaussianMultiArray             */
01439 /*                                                      */
01440 /********************************************************/
01441 
01442 /** \brief Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
01443 
01444     This function computes the Laplacian the given N-dimensional
01445     array with a sequence of second-derivative-of-Gaussian filters at the given
01446     standard deviation <tt>sigma</tt>. Both source and destination arrays
01447     are represented by iterators, shape objects and accessors. Both source and destination 
01448     arrays must have scalar value_type. This function is implemented by calls to
01449     \ref separableConvolveMultiArray() with the appropriate kernels, followed by summation.
01450 
01451     Anisotropic data should be passed with appropriate
01452     \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional
01453     unless the parameter <tt>sigma</tt> is left out.
01454 
01455     <b> Declarations:</b>
01456 
01457     pass arguments explicitly:
01458     \code
01459     namespace vigra {
01460         template <class SrcIterator, class SrcShape, class SrcAccessor,
01461                   class DestIterator, class DestAccessor>
01462         void
01463         laplacianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
01464                                       DestIterator diter, DestAccessor dest,
01465                                       double sigma, const ConvolutionOptions<N> & opt);
01466     }
01467     \endcode
01468 
01469     use argument objects in conjunction with \ref ArgumentObjectFactories :
01470     \code
01471     namespace vigra {
01472         template <class SrcIterator, class SrcShape, class SrcAccessor,
01473                   class DestIterator, class DestAccessor>
01474         void
01475         laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01476                                       pair<DestIterator, DestAccessor> const & dest,
01477                                       double sigma, const ConvolutionOptions<N> & opt);
01478     }
01479     \endcode
01480 
01481     <b> Usage:</b>
01482 
01483     <b>\#include</b> <vigra/multi_convolution.hxx>
01484 
01485     \code
01486     MultiArray<3, float> source(shape);
01487     MultiArray<3, float> laplacian(shape);
01488     ...
01489     // compute Laplacian at scale sigma
01490     laplacianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(laplacian), sigma);
01491     \endcode
01492 
01493     <b> Usage with anisotropic data:</b>
01494 
01495     <b>\#include</b> <vigra/multi_convolution.hxx>
01496 
01497     \code
01498     MultiArray<3, float> source(shape);
01499     MultiArray<3, float> laplacian(shape);
01500     TinyVector<float, 3> step_size;
01501     TinyVector<float, 3> resolution_sigmas;
01502     ...
01503     // compute Laplacian at scale sigma
01504     laplacianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(laplacian), sigma,
01505                                   ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
01506     \endcode
01507 
01508     <b> Required Interface:</b>
01509 
01510     see \ref separableConvolveMultiArray(), in addition:
01511 
01512     \code
01513     int dimension = 0;
01514     VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
01515     \endcode
01516 
01517     \see separableConvolveMultiArray()
01518 */
01519 doxygen_overloaded_function(template <...> void laplacianOfGaussianMultiArray)
01520 
01521 template <class SrcIterator, class SrcShape, class SrcAccessor,
01522           class DestIterator, class DestAccessor>
01523 void
01524 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
01525                               DestIterator di, DestAccessor dest,
01526                               ConvolutionOptions<SrcShape::static_size> const & opt )
01527 { 
01528     using namespace functor;
01529     
01530     typedef typename DestAccessor::value_type DestType;
01531     typedef typename NumericTraits<DestType>::RealPromote KernelType;
01532     typedef typename AccessorTraits<KernelType>::default_accessor DerivativeAccessor;
01533 
01534     static const int N = SrcShape::static_size;
01535     typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
01536     
01537     ParamType params = opt.scaleParams();
01538     ParamType params2(params);
01539 
01540     ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
01541     for (int dim = 0; dim < N; ++dim, ++params)
01542     {
01543         double sigma = params.sigma_scaled("laplacianOfGaussianMultiArray");
01544         plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
01545     }
01546     
01547     SrcShape dshape(shape);
01548     if(opt.to_point != SrcShape())
01549         dshape = opt.to_point - opt.from_point;
01550     
01551     MultiArray<N, KernelType> derivative(dshape);
01552 
01553     // compute 2nd derivatives and sum them up
01554     for (int dim = 0; dim < N; ++dim, ++params2)
01555     {
01556         ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
01557         kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 2, 1.0, opt.window_ratio);
01558         detail::scaleKernel(kernels[dim], 1.0 / sq(params2.step_size()));
01559 
01560         if (dim == 0)
01561         {
01562             separableConvolveMultiArray( si, shape, src, 
01563                                          di, dest, kernels.begin(), opt.from_point, opt.to_point);
01564         }
01565         else
01566         {
01567             separableConvolveMultiArray( si, shape, src, 
01568                                          derivative.traverser_begin(), DerivativeAccessor(), 
01569                                          kernels.begin(), opt.from_point, opt.to_point);
01570             combineTwoMultiArrays(di, dshape, dest, derivative.traverser_begin(), DerivativeAccessor(), 
01571                                   di, dest, Arg1() + Arg2() );
01572         }
01573     }
01574 }
01575 
01576 template <class SrcIterator, class SrcShape, class SrcAccessor,
01577           class DestIterator, class DestAccessor>
01578 void
01579 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
01580                               DestIterator di, DestAccessor dest, double sigma,
01581                               const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01582 {
01583     ConvolutionOptions<SrcShape::static_size> par = opt;
01584     laplacianOfGaussianMultiArray(si, shape, src, di, dest, par.stdDev(sigma));
01585 }
01586 
01587 template <class SrcIterator, class SrcShape, class SrcAccessor,
01588           class DestIterator, class DestAccessor>
01589 inline void
01590 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01591                               pair<DestIterator, DestAccessor> const & dest,
01592                               ConvolutionOptions<SrcShape::static_size> const & opt )
01593 {
01594     laplacianOfGaussianMultiArray( source.first, source.second, source.third,
01595                                    dest.first, dest.second, opt );
01596 }
01597 
01598 template <class SrcIterator, class SrcShape, class SrcAccessor,
01599           class DestIterator, class DestAccessor>
01600 inline void
01601 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01602                               pair<DestIterator, DestAccessor> const & dest,
01603                               double sigma,
01604                               const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01605 {
01606     laplacianOfGaussianMultiArray( source.first, source.second, source.third,
01607                                    dest.first, dest.second,  sigma, opt );
01608 }
01609 
01610 /********************************************************/
01611 /*                                                      */
01612 /*              hessianOfGaussianMultiArray             */
01613 /*                                                      */
01614 /********************************************************/
01615 
01616 /** \brief Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
01617 
01618     This function computes the Hessian matrix the given scalar N-dimensional
01619     array with a sequence of second-derivative-of-Gaussian filters at the given
01620     standard deviation <tt>sigma</tt>. Both source and destination arrays
01621     are represented by iterators, shape objects and accessors. The destination array must 
01622     have a vector valued element type with N*(N+1)/2 elements (it represents the
01623     upper triangular part of the symmetric Hessian matrix). This function is implemented by calls to
01624     \ref separableConvolveMultiArray() with the appropriate kernels.
01625 
01626     Anisotropic data should be passed with appropriate
01627     \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional
01628     unless the parameter <tt>sigma</tt> is left out.
01629 
01630     <b> Declarations:</b>
01631 
01632     pass arguments explicitly:
01633     \code
01634     namespace vigra {
01635         template <class SrcIterator, class SrcShape, class SrcAccessor,
01636                   class DestIterator, class DestAccessor>
01637         void
01638         hessianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
01639                                     DestIterator diter, DestAccessor dest,
01640                                     double sigma, const ConvolutionOptions<N> & opt);
01641     }
01642     \endcode
01643 
01644     use argument objects in conjunction with \ref ArgumentObjectFactories :
01645     \code
01646     namespace vigra {
01647         template <class SrcIterator, class SrcShape, class SrcAccessor,
01648                   class DestIterator, class DestAccessor>
01649         void
01650         hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01651                                     pair<DestIterator, DestAccessor> const & dest,
01652                                     double sigma, const ConvolutionOptions<N> & opt);
01653     }
01654     \endcode
01655 
01656     <b> Usage:</b>
01657 
01658     <b>\#include</b> <vigra/multi_convolution.hxx>
01659 
01660     \code
01661     MultiArray<3, float> source(shape);
01662     MultiArray<3, TinyVector<float, 6> > dest(shape);
01663     ...
01664     // compute Hessian at scale sigma
01665     hessianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
01666     \endcode
01667 
01668     <b> Usage with anisotropic data:</b>
01669 
01670     <b>\#include</b> <vigra/multi_convolution.hxx>
01671 
01672     \code
01673     MultiArray<3, float> source(shape);
01674     MultiArray<3, TinyVector<float, 6> > dest(shape);
01675     TinyVector<float, 3> step_size;
01676     TinyVector<float, 3> resolution_sigmas;
01677     ...
01678     // compute Hessian at scale sigma
01679     hessianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma,
01680                                 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
01681     \endcode
01682 
01683     <b> Required Interface:</b>
01684 
01685     see \ref separableConvolveMultiArray(), in addition:
01686 
01687     \code
01688     int dimension = 0;
01689     VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
01690     \endcode
01691 
01692     \see separableConvolveMultiArray(), vectorToTensorMultiArray()
01693 */
01694 doxygen_overloaded_function(template <...> void hessianOfGaussianMultiArray)
01695 
01696 template <class SrcIterator, class SrcShape, class SrcAccessor,
01697           class DestIterator, class DestAccessor>
01698 void
01699 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
01700                             DestIterator di, DestAccessor dest,
01701                             ConvolutionOptions<SrcShape::static_size> const & opt )
01702 { 
01703     typedef typename DestAccessor::value_type DestType;
01704     typedef typename DestType::value_type     DestValueType;
01705     typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
01706 
01707     static const int N = SrcShape::static_size;
01708     static const int M = N*(N+1)/2;
01709     typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
01710     
01711     for(int k=0; k<N; ++k)
01712         if(shape[k] <=0)
01713             return;
01714 
01715     vigra_precondition(M == (int)dest.size(di),
01716         "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
01717 
01718     ParamType params_init = opt.scaleParams();
01719 
01720     ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
01721     ParamType params(params_init);
01722     for (int dim = 0; dim < N; ++dim, ++params)
01723     {
01724         double sigma = params.sigma_scaled("hessianOfGaussianMultiArray");
01725         plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
01726     }
01727 
01728     typedef VectorElementAccessor<DestAccessor> ElementAccessor;
01729 
01730     // compute elements of the Hessian matrix
01731     ParamType params_i(params_init);
01732     for (int b=0, i=0; i<N; ++i, ++params_i)
01733     {
01734         ParamType params_j(params_i);
01735         for (int j=i; j<N; ++j, ++b, ++params_j)
01736         {
01737             ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
01738             if(i == j)
01739             {
01740                 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 2, 1.0, opt.window_ratio);
01741             }
01742             else
01743             {
01744                 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 1, 1.0, opt.window_ratio);
01745                 kernels[j].initGaussianDerivative(params_j.sigma_scaled(), 1, 1.0, opt.window_ratio);
01746             }
01747             detail::scaleKernel(kernels[i], 1 / params_i.step_size());
01748             detail::scaleKernel(kernels[j], 1 / params_j.step_size());
01749             separableConvolveMultiArray(si, shape, src, di, ElementAccessor(b, dest),
01750                                         kernels.begin(), opt.from_point, opt.to_point);
01751         }
01752     }
01753 }
01754 
01755 template <class SrcIterator, class SrcShape, class SrcAccessor,
01756           class DestIterator, class DestAccessor>
01757 inline void
01758 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
01759                             DestIterator di, DestAccessor dest, double sigma,
01760                             const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01761 {
01762     ConvolutionOptions<SrcShape::static_size> par = opt;
01763     hessianOfGaussianMultiArray(si, shape, src, di, dest, par.stdDev(sigma));
01764 }
01765 
01766 template <class SrcIterator, class SrcShape, class SrcAccessor,
01767           class DestIterator, class DestAccessor>
01768 inline void
01769 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01770                             pair<DestIterator, DestAccessor> const & dest,
01771                             ConvolutionOptions<SrcShape::static_size> const & opt )
01772 {
01773     hessianOfGaussianMultiArray( source.first, source.second, source.third,
01774                                  dest.first, dest.second, opt );
01775 }
01776 
01777 template <class SrcIterator, class SrcShape, class SrcAccessor,
01778           class DestIterator, class DestAccessor>
01779 inline void
01780 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01781                             pair<DestIterator, DestAccessor> const & dest,
01782                             double sigma,
01783                             const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01784 {
01785     hessianOfGaussianMultiArray( source.first, source.second, source.third,
01786                                  dest.first, dest.second, sigma, opt );
01787 }
01788 
01789 namespace detail {
01790 
01791 template<int N, class VectorType>
01792 struct StructurTensorFunctor
01793 {
01794     typedef VectorType result_type;
01795     typedef typename VectorType::value_type ValueType;
01796     
01797     template <class T>
01798     VectorType operator()(T const & in) const
01799     {
01800         VectorType res;
01801         for(int b=0, i=0; i<N; ++i)
01802         {
01803             for(int j=i; j<N; ++j, ++b)
01804             {
01805                 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
01806             }
01807         }
01808         return res;
01809     }
01810 };
01811 
01812 } // namespace detail
01813 
01814 /********************************************************/
01815 /*                                                      */
01816 /*               structureTensorMultiArray              */
01817 /*                                                      */
01818 /********************************************************/
01819 
01820 /** \brief Calculate th structure tensor of a multi-dimensional arrays.
01821 
01822     This function computes the gradient (outer product) tensor for each element
01823     of the given N-dimensional array with first-derivative-of-Gaussian filters at 
01824     the given <tt>innerScale</tt>, followed by Gaussian smoothing at <tt>outerScale</tt>.
01825     Both source and destination arrays are represented by iterators, shape objects and 
01826     accessors. The destination array must have a vector valued pixel type with 
01827     N*(N+1)/2 elements (it represents the upper triangular part of the symmetric 
01828     structure tensor matrix). If the source array is also vector valued, the 
01829     resulting structure tensor is the sum of the individual tensors for each channel.
01830     This function is implemented by calls to
01831     \ref separableConvolveMultiArray() with the appropriate kernels.
01832 
01833     Anisotropic data should be passed with appropriate
01834     \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional
01835     unless the parameters <tt>innerScale</tt> and <tt>outerScale</tt> are
01836     both left out.
01837 
01838     <b> Declarations:</b>
01839 
01840     pass arguments explicitly:
01841     \code
01842     namespace vigra {
01843         template <class SrcIterator, class SrcShape, class SrcAccessor,
01844                   class DestIterator, class DestAccessor>
01845         void
01846         structureTensorMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
01847                                   DestIterator diter, DestAccessor dest,
01848                                   double innerScale, double outerScale,
01849                                   const ConvolutionOptions<N> & opt);
01850     }
01851     \endcode
01852 
01853     use argument objects in conjunction with \ref ArgumentObjectFactories :
01854     \code
01855     namespace vigra {
01856         template <class SrcIterator, class SrcShape, class SrcAccessor,
01857                   class DestIterator, class DestAccessor>
01858         void
01859         structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01860                                   pair<DestIterator, DestAccessor> const & dest,
01861                                   double innerScale, double outerScale,
01862                                   const ConvolutionOptions<N> & opt);
01863     }
01864     \endcode
01865 
01866     <b> Usage:</b>
01867 
01868     <b>\#include</b> <vigra/multi_convolution.hxx>
01869 
01870     \code
01871     MultiArray<3, RGBValue<float> > source(shape);
01872     MultiArray<3, TinyVector<float, 6> > dest(shape);
01873     ...
01874     // compute structure tensor at scales innerScale and outerScale
01875     structureTensorMultiArray(srcMultiArrayRange(source), destMultiArray(dest), innerScale, outerScale);
01876     \endcode
01877 
01878     <b> Usage with anisotropic data:</b>
01879 
01880     <b>\#include</b> <vigra/multi_convolution.hxx>
01881 
01882     \code
01883     MultiArray<3, RGBValue<float> > source(shape);
01884     MultiArray<3, TinyVector<float, 6> > dest(shape);
01885     TinyVector<float, 3> step_size;
01886     TinyVector<float, 3> resolution_sigmas;
01887     ...
01888     // compute structure tensor at scales innerScale and outerScale
01889     structureTensorMultiArray(srcMultiArrayRange(source), destMultiArray(dest), innerScale, outerScale,
01890                               ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
01891     \endcode
01892 
01893     <b> Required Interface:</b>
01894 
01895     see \ref separableConvolveMultiArray(), in addition:
01896 
01897     \code
01898     int dimension = 0;
01899     VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
01900     \endcode
01901 
01902     \see separableConvolveMultiArray(), vectorToTensorMultiArray()
01903 */
01904 doxygen_overloaded_function(template <...> void structureTensorMultiArray)
01905 
01906 template <class SrcIterator, class SrcShape, class SrcAccessor,
01907           class DestIterator, class DestAccessor>
01908 void
01909 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
01910                           DestIterator di, DestAccessor dest, 
01911                           ConvolutionOptions<SrcShape::static_size> const & opt)
01912 {
01913     static const int N = SrcShape::static_size;
01914     static const int M = N*(N+1)/2;
01915     
01916     typedef typename DestAccessor::value_type DestType;
01917     typedef typename DestType::value_type     DestValueType;
01918     typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
01919     typedef TinyVector<KernelType, N> GradientVector;
01920     typedef typename AccessorTraits<GradientVector>::default_accessor GradientAccessor;
01921     typedef typename AccessorTraits<DestType>::default_accessor GradientTensorAccessor;
01922 
01923     for(int k=0; k<N; ++k)
01924         if(shape[k] <=0)
01925             return;
01926 
01927     vigra_precondition(M == (int)dest.size(di),
01928         "structureTensorMultiArray(): Wrong number of channels in output array.");
01929         
01930     ConvolutionOptions<N> innerOptions = opt;
01931     ConvolutionOptions<N> outerOptions = opt.outerOptions();
01932     typename ConvolutionOptions<N>::ScaleIterator params = outerOptions.scaleParams();
01933     
01934     SrcShape gradientShape(shape);
01935     if(opt.to_point != SrcShape())
01936     {
01937         for(int k=0; k<N; ++k, ++params)
01938         {
01939             Kernel1D<double> gauss;
01940             gauss.initGaussian(params.sigma_scaled("structureTensorMultiArray"), 1.0, opt.window_ratio);
01941             int dilation = gauss.right();
01942             innerOptions.from_point[k] = std::max<MultiArrayIndex>(0, opt.from_point[k] - dilation);
01943             innerOptions.to_point[k] = std::min<MultiArrayIndex>(shape[k], opt.to_point[k] + dilation);
01944         }
01945         outerOptions.from_point -= innerOptions.from_point;
01946         outerOptions.to_point -= innerOptions.from_point;
01947         gradientShape = innerOptions.to_point - innerOptions.from_point;
01948     }
01949 
01950     MultiArray<N, GradientVector> gradient(gradientShape);
01951     MultiArray<N, DestType> gradientTensor(gradientShape);
01952     gaussianGradientMultiArray(si, shape, src, 
01953                                gradient.traverser_begin(), GradientAccessor(), 
01954                                innerOptions,
01955                                "structureTensorMultiArray");
01956 
01957     transformMultiArray(gradient.traverser_begin(), gradientShape, GradientAccessor(), 
01958                         gradientTensor.traverser_begin(), GradientTensorAccessor(), 
01959                         detail::StructurTensorFunctor<N, DestType>());
01960 
01961     gaussianSmoothMultiArray(gradientTensor.traverser_begin(), gradientShape, GradientTensorAccessor(), 
01962                              di, dest, outerOptions,
01963                              "structureTensorMultiArray");
01964 }
01965 
01966 template <class SrcIterator, class SrcShape, class SrcAccessor,
01967           class DestIterator, class DestAccessor>
01968 inline void
01969 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
01970                           DestIterator di, DestAccessor dest,
01971                           double innerScale, double outerScale,
01972                           const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01973 {
01974     ConvolutionOptions<SrcShape::static_size> par = opt;
01975     structureTensorMultiArray(si, shape, src, di, dest,
01976                               par.stdDev(innerScale).outerScale(outerScale));
01977 }
01978 
01979 template <class SrcIterator, class SrcShape, class SrcAccessor,
01980           class DestIterator, class DestAccessor>
01981 inline void
01982 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01983                           pair<DestIterator, DestAccessor> const & dest, 
01984                           ConvolutionOptions<SrcShape::static_size> const & opt )
01985 {
01986     structureTensorMultiArray( source.first, source.second, source.third,
01987                                dest.first, dest.second, opt );
01988 }
01989 
01990 
01991 template <class SrcIterator, class SrcShape, class SrcAccessor,
01992           class DestIterator, class DestAccessor>
01993 inline void
01994 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
01995                           pair<DestIterator, DestAccessor> const & dest,
01996                           double innerScale, double outerScale,
01997                           const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
01998 {
01999     structureTensorMultiArray( source.first, source.second, source.third,
02000                                dest.first, dest.second,
02001                                innerScale, outerScale, opt);
02002 }
02003 
02004 //@}
02005 
02006 } //-- namespace vigra
02007 
02008 
02009 #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.8.0 (20 Sep 2011)