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

vigra/multi_convolution.hxx |

    00001 //-- -*- c++ -*-
    00002 /************************************************************************/
    00003 /*                                                                      */
    00004 /*               Copyright 2003 by Christian-Dennis Rahn                */
    00005 /*                        and Ullrich Koethe                            */
    00006 /*                                                                      */
    00007 /*    This file is part of the VIGRA computer vision library.           */
    00008 /*    The VIGRA Website is                                              */
    00009 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
    00010 /*    Please direct questions, bug reports, and contributions to        */
    00011 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00012 /*        vigra@informatik.uni-hamburg.de                               */
    00013 /*                                                                      */
    00014 /*    Permission is hereby granted, free of charge, to any person       */
    00015 /*    obtaining a copy of this software and associated documentation    */
    00016 /*    files (the "Software"), to deal in the Software without           */
    00017 /*    restriction, including without limitation the rights to use,      */
    00018 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
    00019 /*    sell copies of the Software, and to permit persons to whom the    */
    00020 /*    Software is furnished to do so, subject to the following          */
    00021 /*    conditions:                                                       */
    00022 /*                                                                      */
    00023 /*    The above copyright notice and this permission notice shall be    */
    00024 /*    included in all copies or substantial portions of the             */
    00025 /*    Software.                                                         */
    00026 /*                                                                      */
    00027 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
    00028 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
    00029 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
    00030 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
    00031 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
    00032 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
    00033 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
    00034 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
    00035 /*                                                                      */
    00036 /************************************************************************/
    00037
    00038 #ifndef VIGRA_MULTI_CONVOLUTION_H
    00039 #define VIGRA_MULTI_CONVOLUTION_H
    00040
    00041 #include "separableconvolution.hxx"
    00042 #include "array_vector.hxx"
    00043 #include "multi_array.hxx"
    00044 #include "accessor.hxx"
    00045 #include "numerictraits.hxx"
    00046 #include "navigator.hxx"
    00047 #include "metaprogramming.hxx"
    00048 #include "multi_pointoperators.hxx"
    00049 #include "functorexpression.hxx"
    00050 #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.9.0 (Tue Nov 6 2012)