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

| vigra/boundarytensor.hxx |
--- | --- | ---

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*               Copyright 1998-2004 by Ullrich Koethe                  */
    00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
    00005 /*                                                                      */
    00006 /*    This file is part of the VIGRA computer vision library.           */
    00007 /*    ( Version 1.3.2, Jan 27 2005 )                                    */
    00008 /*    You may use, modify, and distribute this software according       */
    00009 /*    to the terms stated in the LICENSE file included in               */
    00010 /*    the VIGRA distribution.                                           */
    00011 /*                                                                      */
    00012 /*    The VIGRA Website is                                              */
    00013 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
    00014 /*    Please direct questions, bug reports, and contributions to        */
    00015 /*        koethe@informatik.uni-hamburg.de                              */
    00016 /*                                                                      */
    00017 /*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
    00018 /*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
    00019 /*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
    00020 /*                                                                      */
    00021 /************************************************************************/
    00022
    00023
    00024 #ifndef VIGRA_BOUNDARYTENSOR_HXX
    00025 #define VIGRA_BOUNDARYTENSOR_HXX
    00026
    00027 #include <cmath>
    00028 #include <functional>
    00029 #include "vigra/utilities.hxx"
    00030 #include "vigra/array_vector.hxx"
    00031 #include "vigra/basicimage.hxx"
    00032 #include "vigra/combineimages.hxx"
    00033 #include "vigra/numerictraits.hxx"
    00034 #include "vigra/convolution.hxx"
    00035
    00036 namespace vigra {
    00037
    00038 namespace detail {
    00039
    00040 /***********************************************************************/
    00041
    00042 typedef ArrayVector<Kernel1D<double> > KernelArray;
    00043
    00044 void
    00045 initGaussianPolarFilters1(double std_dev, KernelArray & k)
    00046 {
    00047     typedef KernelArray::value_type Kernel;
    00048     typedef Kernel::iterator iterator;
    00049
    00050     vigra_precondition(std_dev >= 0.0,
    00051               "initGaussianPolarFilter1(): "
    00052               "Standard deviation must be >= 0.");
    00053
    00054     k.resize(4);
    00055
    00056     int radius = (int)(4.0*std_dev + 0.5);
    00057     std_dev *= 1.08179074376;
    00058     double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev;  // norm
    00059     double a = 0.558868151788 / VIGRA_CSTD::pow(std_dev, 5);
    00060     double b = -2.04251639729 / VIGRA_CSTD::pow(std_dev, 3);
    00061     double sigma22 = -0.5 / std_dev / std_dev;
    00062
    00063
    00064     for(unsigned int i=0; i<k.size(); ++i)
    00065     {
    00066         k[i].initExplicitly(-radius, radius);
    00067         k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
    00068     }
    00069
    00070     int ix;
    00071     iterator c = k[0].center();
    00072     for(ix=-radius; ix<=radius; ++ix)
    00073     {
    00074         double x = (double)ix;
    00075         c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
    00076     }
    00077
    00078     c = k[1].center();
    00079     for(ix=-radius; ix<=radius; ++ix)
    00080     {
    00081         double x = (double)ix;
    00082         c[ix] = f * x * VIGRA_CSTD::exp(sigma22 * x * x);
    00083     }
    00084
    00085     c = k[2].center();
    00086     double b2 = b / 3.0;
    00087     for(ix=-radius; ix<=radius; ++ix)
    00088     {
    00089         double x = (double)ix;
    00090         c[ix] = f * (b2 + a * x * x) * VIGRA_CSTD::exp(sigma22 * x * x);
    00091     }
    00092
    00093     c = k[3].center();
    00094     for(ix=-radius; ix<=radius; ++ix)
    00095     {
    00096         double x = (double)ix;
    00097         c[ix] = f * x * (b + a * x * x) * VIGRA_CSTD::exp(sigma22 * x * x);
    00098     }
    00099 }
    00100
    00101 void
    00102 initGaussianPolarFilters2(double std_dev, KernelArray & k)
    00103 {
    00104     typedef Kernel1D<double>::iterator iterator;
    00105
    00106     vigra_precondition(std_dev >= 0.0,
    00107               "initGaussianPolarFilter2(): "
    00108               "Standard deviation must be >= 0.");
    00109
    00110     k.resize(3);
    00111
    00112     int radius = (int)(4.0*std_dev + 0.5);
    00113     double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev;  // norm
    00114     double sigma2 = std_dev*std_dev;
    00115     double sigma22 = -0.5 / sigma2;
    00116
    00117     for(unsigned int i=0; i<k.size(); ++i)
    00118     {
    00119         k[i].initExplicitly(-radius, radius);
    00120         k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
    00121     }
    00122
    00123     int ix;
    00124     iterator c = k[0].center();
    00125     for(ix=-radius; ix<=radius; ++ix)
    00126     {
    00127         double x = (double)ix;
    00128         c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
    00129     }
    00130
    00131     c = k[1].center();
    00132     double f1 = f / sigma2;
    00133     for(ix=-radius; ix<=radius; ++ix)
    00134     {
    00135         double x = (double)ix;
    00136         c[ix] = f1 * x * VIGRA_CSTD::exp(sigma22 * x * x);
    00137     }
    00138
    00139     c = k[2].center();
    00140     double f2 = f / (sigma2 * sigma2);
    00141     for(ix=-radius; ix<=radius; ++ix)
    00142     {
    00143         double x = (double)ix;
    00144         c[ix] = f2 * (x * x - sigma2) * VIGRA_CSTD::exp(sigma22 * x * x);
    00145     }
    00146 }
    00147
    00148 void
    00149 initGaussianPolarFilters3(double std_dev, KernelArray & k)
    00150 {
    00151     typedef Kernel1D<double>::iterator iterator;
    00152
    00153     vigra_precondition(std_dev >= 0.0,
    00154               "initGaussianPolarFilter3(): "
    00155               "Standard deviation must be >= 0.");
    00156
    00157     k.resize(4);
    00158
    00159     int radius = (int)(4.0*std_dev + 0.5);
    00160     std_dev *= 1.15470053838;
    00161     double sigma22 = -0.5 / std_dev / std_dev;
    00162     double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev;  // norm
    00163     double a = 0.883887052922 / VIGRA_CSTD::pow(std_dev, 5);
    00164
    00165     for(unsigned int i=0; i<k.size(); ++i)
    00166     {
    00167         k[i].initExplicitly(-radius, radius);
    00168         k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
    00169     }
    00170
    00171     double b = -1.3786348292 / VIGRA_CSTD::pow(std_dev, 3);
    00172
    00173     int ix;
    00174     iterator c = k[0].center();
    00175     for(ix=-radius; ix<=radius; ++ix)
    00176     {
    00177         double x = (double)ix;
    00178         c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
    00179     }
    00180
    00181     c = k[1].center();
    00182     for(ix=-radius; ix<=radius; ++ix)
    00183     {
    00184         double x = (double)ix;
    00185         c[ix] = f * x * VIGRA_CSTD::exp(sigma22 * x * x);
    00186     }
    00187
    00188     c = k[2].center();
    00189     double a2 = 3.0 * a;
    00190     for(ix=-radius; ix<=radius; ++ix)
    00191     {
    00192         double x = (double)ix;
    00193         c[ix] = f * a2 * x * x * VIGRA_CSTD::exp(sigma22 * x * x);
    00194     }
    00195
    00196     c = k[3].center();
    00197     for(ix=-radius; ix<=radius; ++ix)
    00198     {
    00199         double x = (double)ix;
    00200         c[ix] = f * a * x * x * x * VIGRA_CSTD::exp(sigma22 * x * x);
    00201     }
    00202 }
    00203
    00204 template <class SrcIterator, class SrcAccessor,
    00205           class DestIterator, class DestAccessor>
    00206 void
    00207 evenPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
    00208                  DestIterator dupperleft, DestAccessor dest,
    00209                  double scale, bool noLaplacian)
    00210 {
    00211     vigra_precondition(dest.size(dupperleft) == 3,
    00212                        "evenPolarFilters(): image for even output must have 3 bands.");
    00213
    00214     int w = slowerright.x - supperleft.x;
    00215     int h = slowerright.y - supperleft.y;
    00216
    00217     typedef typename
    00218        NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    00219     typedef BasicImage<TinyVector<TmpType, 3> > TmpImage;
    00220     typedef typename TmpImage::traverser TmpTraverser;
    00221     TmpImage t(w, h);
    00222
    00223     KernelArray k2;
    00224     initGaussianPolarFilters2(scale, k2);
    00225
    00226     // calculate filter responses for even filters
    00227     VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor());
    00228     convolveImage(srcIterRange(supperleft, slowerright, src),
    00229                   destImage(t, tmpBand), k2[2], k2[0]);
    00230     tmpBand.setIndex(1);
    00231     convolveImage(srcIterRange(supperleft, slowerright, src),
    00232                   destImage(t, tmpBand), k2[1], k2[1]);
    00233     tmpBand.setIndex(2);
    00234     convolveImage(srcIterRange(supperleft, slowerright, src),
    00235                   destImage(t, tmpBand), k2[0], k2[2]);
    00236
    00237     // create even tensor from filter responses
    00238     TmpTraverser tul(t.upperLeft());
    00239     TmpTraverser tlr(t.lowerRight());
    00240     for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y)
    00241     {
    00242         typename TmpTraverser::row_iterator tr = tul.rowIterator();
    00243         typename TmpTraverser::row_iterator trend = tr + w;
    00244         typename DestIterator::row_iterator d = dupperleft.rowIterator();
    00245         if(noLaplacian)
    00246         {
    00247             for(; tr != trend; ++tr, ++d)
    00248             {
    00249                 TmpType v = 0.5*sq((*tr)[0]-(*tr)[2]) + 2.0*sq((*tr)[1]);
    00250                 dest.setComponent(v, d, 0);
    00251                 dest.setComponent(0, d, 1);
    00252                 dest.setComponent(v, d, 2);
    00253             }
    00254         }
    00255         else
    00256         {
    00257             for(; tr != trend; ++tr, ++d)
    00258             {
    00259                 dest.setComponent(sq((*tr)[0]) + sq((*tr)[1]), d, 0);
    00260                 dest.setComponent(-(*tr)[1] * ((*tr)[0] + (*tr)[2]), d, 1);
    00261                 dest.setComponent(sq((*tr)[1]) + sq((*tr)[2]), d, 2);
    00262             }
    00263         }
    00264     }
    00265 }
    00266
    00267 template <class SrcIterator, class SrcAccessor,
    00268           class DestIterator, class DestAccessor>
    00269 void
    00270 oddPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
    00271                 DestIterator dupperleft, DestAccessor dest,
    00272                 double scale, bool addResult)
    00273 {
    00274     vigra_precondition(dest.size(dupperleft) == 3,
    00275                        "oddPolarFilters(): image for odd output must have 3 bands.");
    00276
    00277     int w = slowerright.x - supperleft.x;
    00278     int h = slowerright.y - supperleft.y;
    00279
    00280     typedef typename
    00281        NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    00282     typedef BasicImage<TinyVector<TmpType, 4> > TmpImage;
    00283     typedef typename TmpImage::traverser TmpTraverser;
    00284     TmpImage t(w, h);
    00285
    00286     detail::KernelArray k1;
    00287     detail::initGaussianPolarFilters1(scale, k1);
    00288
    00289     // calculate filter responses for odd filters
    00290     VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor());
    00291     convolveImage(srcIterRange(supperleft, slowerright, src),
    00292                   destImage(t, tmpBand), k1[3], k1[0]);
    00293     tmpBand.setIndex(1);
    00294     convolveImage(srcIterRange(supperleft, slowerright, src),
    00295                   destImage(t, tmpBand), k1[2], k1[1]);
    00296     tmpBand.setIndex(2);
    00297     convolveImage(srcIterRange(supperleft, slowerright, src),
    00298                   destImage(t, tmpBand), k1[1], k1[2]);
    00299     tmpBand.setIndex(3);
    00300     convolveImage(srcIterRange(supperleft, slowerright, src),
    00301                   destImage(t, tmpBand), k1[0], k1[3]);
    00302
    00303     // create odd tensor from filter responses
    00304     TmpTraverser tul(t.upperLeft());
    00305     TmpTraverser tlr(t.lowerRight());
    00306     for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y)
    00307     {
    00308         typename TmpTraverser::row_iterator tr = tul.rowIterator();
    00309         typename TmpTraverser::row_iterator trend = tr + w;
    00310         typename DestIterator::row_iterator d = dupperleft.rowIterator();
    00311         if(addResult)
    00312         {
    00313             for(; tr != trend; ++tr, ++d)
    00314             {
    00315                 TmpType d0 = (*tr)[0] + (*tr)[2];
    00316                 TmpType d1 = -(*tr)[1] - (*tr)[3];
    00317
    00318                 dest.setComponent(dest.getComponent(d, 0) + sq(d0), d, 0);
    00319                 dest.setComponent(dest.getComponent(d, 1) + d0 * d1, d, 1);
    00320                 dest.setComponent(dest.getComponent(d, 2) + sq(d1), d, 2);
    00321             }
    00322         }
    00323         else
    00324         {
    00325             for(; tr != trend; ++tr, ++d)
    00326             {
    00327                 TmpType d0 = (*tr)[0] + (*tr)[2];
    00328                 TmpType d1 = -(*tr)[1] - (*tr)[3];
    00329
    00330                 dest.setComponent(sq(d0), d, 0);
    00331                 dest.setComponent(d0 * d1, d, 1);
    00332                 dest.setComponent(sq(d1), d, 2);
    00333             }
    00334         }
    00335     }
    00336 }
    00337
    00338 } // namespace detail
    00339
    00340 /** \addtogroup CommonConvolutionFilters Common Filters
    00341 */
    00342 //@{
    00343
    00344 /********************************************************/
    00345 /*                                                      */
    00346 /*                   rieszTransformOfLOG                */
    00347 /*                                                      */
    00348 /********************************************************/
    00349
    00350 /** \brief Calculate Riesz transforms of the Laplacian of Gaussian.
    00351
    00352     The Riesz transforms of the Laplacian of Gaussian have the following transfer
    00353     functions (defined in a polar coordinate representation of the frequency domain):
    00354
    00355     \f[
    00356         F_{\sigma}(r, \phi)=(i \cos \phi)^n (i \sin \phi)^m r^2 e^{-r^2 \sigma^2 / 2}
    00357     \f]
    00358
    00359     where <i>n</i> = <tt>xorder</tt> and <i>m</i> = <tt>yorder</tt> determine th e
    00360     order of the transform, and <tt>sigma &gt; 0</tt> is the scale of the Laplacian
    00361     of Gaussian. This function computes a good spatial domain approximation of
    00362     these transforms for <tt>xorder + yorder &lt;= 2</tt>. The filter responses may be used
    00363     to calculate the monogenic signal or the boundary tensor.
    00364
    00365     <b> Declarations:</b>
    00366
    00367     pass arguments explicitly:
    00368     \code
    00369     namespace vigra {
    00370         template <class SrcIterator, class SrcAccessor,
    00371                 class DestIterator, class DestAccessor>
    00372         void rieszTransformOfLOG(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
    00373                                  DestIterator dupperleft, DestAccessor dest,
    00374                                  double scale, unsigned int xorder, unsigned int yorder);
    00375     }
    00376     \endcode
    00377
    00378
    00379     use argument objects in conjunction with \ref ArgumentObjectFactories:
    00380     \code
    00381     namespace vigra {
    00382         template <class SrcIterator, class SrcAccessor,
    00383                 class DestIterator, class DestAccessor>
    00384         void rieszTransformOfLOG(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00385                                  pair<DestIterator, DestAccessor> dest,
    00386                                  double scale, unsigned int xorder, unsigned int yorder);
    00387     }
    00388     \endcode
    00389
    00390     <b> Usage:</b>
    00391
    00392     <b>\#include</b> "<a href="boundarytensor_8hxx-source.html">vigra/boundarytensor.hxx</a>"
    00393
    00394     \code
    00395     FImage impulse(17,17), res(17, 17);
    00396     impulse(8,8) = 1.0;
    00397
    00398     // calculate the impulse response of the first order Riesz transform in x-direction
    00399     rieszTransformOfLOG(srcImageRange(impulse), destImage(res), 2.0, 1, 0);
    00400     \endcode
    00401
    00402 */
    00403 template <class SrcIterator, class SrcAccessor,
    00404           class DestIterator, class DestAccessor>
    00405 void rieszTransformOfLOG(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
    00406                          DestIterator dupperleft, DestAccessor dest,
    00407                          double scale, unsigned int xorder, unsigned int yorder)
    00408 {
    00409     unsigned int order = xorder + yorder;
    00410
    00411     vigra_precondition(order <= 2,
    00412             "rieszTransformOfLOG(): can only compute Riesz transforms up to order 2.");
    00413     vigra_precondition(scale > 0.0,
    00414             "rieszTransformOfLOG(): scale must be positive.");
    00415
    00416     int w = slowerright.x - supperleft.x;
    00417     int h = slowerright.y - supperleft.y;
    00418
    00419     typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    00420     typedef BasicImage<TmpType> TmpImage;
    00421
    00422     switch(order)
    00423     {
    00424         case 0:
    00425         {
    00426             detail::KernelArray k2;
    00427             detail::initGaussianPolarFilters2(scale, k2);
    00428
    00429             TmpImage t1(w, h), t2(w, h);
    00430
    00431             convolveImage(srcIterRange(supperleft, slowerright, src),
    00432                           destImage(t1), k2[2], k2[0]);
    00433             convolveImage(srcIterRange(supperleft, slowerright, src),
    00434                           destImage(t2), k2[0], k2[2]);
    00435             combineTwoImages(srcImageRange(t1), srcImage(t2),
    00436                              destIter(dupperleft, dest), std::plus<TmpType>());
    00437             break;
    00438         }
    00439         case 1:
    00440         {
    00441             detail::KernelArray k1;
    00442             detail::initGaussianPolarFilters1(scale, k1);
    00443
    00444             TmpImage t1(w, h), t2(w, h);
    00445
    00446             if(xorder == 1)
    00447             {
    00448                 convolveImage(srcIterRange(supperleft, slowerright, src),
    00449                             destImage(t1), k1[3], k1[0]);
    00450                 convolveImage(srcIterRange(supperleft, slowerright, src),
    00451                             destImage(t2), k1[1], k1[2]);
    00452             }
    00453             else
    00454             {
    00455                 convolveImage(srcIterRange(supperleft, slowerright, src),
    00456                             destImage(t1), k1[0], k1[3]);
    00457                 convolveImage(srcIterRange(supperleft, slowerright, src),
    00458                             destImage(t2), k1[2], k1[1]);
    00459             }
    00460             combineTwoImages(srcImageRange(t1), srcImage(t2),
    00461                              destIter(dupperleft, dest), std::plus<TmpType>());
    00462             break;
    00463         }
    00464         case 2:
    00465         {
    00466             detail::KernelArray k2;
    00467             detail::initGaussianPolarFilters2(scale, k2);
    00468
    00469             convolveImage(srcIterRange(supperleft, slowerright, src),
    00470                           destIter(dupperleft, dest), k2[xorder], k2[yorder]);
    00471             break;
    00472         }
    00473         /* for test purposes only: compute 3rd order polar filters */
    00474         case 3:
    00475         {
    00476             detail::KernelArray k3;
    00477             detail::initGaussianPolarFilters3(scale, k3);
    00478
    00479             TmpImage t1(w, h), t2(w, h);
    00480
    00481             if(xorder == 3)
    00482             {
    00483                 convolveImage(srcIterRange(supperleft, slowerright, src),
    00484                             destImage(t1), k3[3], k3[0]);
    00485                 convolveImage(srcIterRange(supperleft, slowerright, src),
    00486                             destImage(t2), k3[1], k3[2]);
    00487             }
    00488             else
    00489             {
    00490                 convolveImage(srcIterRange(supperleft, slowerright, src),
    00491                             destImage(t1), k3[0], k3[3]);
    00492                 convolveImage(srcIterRange(supperleft, slowerright, src),
    00493                             destImage(t2), k3[2], k3[1]);
    00494             }
    00495             combineTwoImages(srcImageRange(t1), srcImage(t2),
    00496                              destIter(dupperleft, dest), std::minus<TmpType>());
    00497             break;
    00498         }
    00499     }
    00500 }
    00501
    00502 template <class SrcIterator, class SrcAccessor,
    00503           class DestIterator, class DestAccessor>
    00504 inline
    00505 void rieszTransformOfLOG(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00506                          pair<DestIterator, DestAccessor> dest,
    00507                          double scale, unsigned int xorder, unsigned int yorder)
    00508 {
    00509     rieszTransformOfLOG(src.first, src.second, src.third, dest.first, dest.second,
    00510                         scale, xorder, yorder);
    00511 }
    00512 //@}
    00513
    00514 /** \addtogroup TensorImaging Tensor Image Processing
    00515 */
    00516 //@{
    00517
    00518 /********************************************************/
    00519 /*                                                      */
    00520 /*                     boundaryTensor                   */
    00521 /*                                                      */
    00522 /********************************************************/
    00523
    00524 /** \brief Calculate the boundary tensor for a scalar valued image.
    00525
    00526     These functions calculates a spatial domain approximation of
    00527     the boundary tensor as described in
    00528
    00529     U. Köthe: <a href="http://kogs-www.informatik.uni-hamburg.de/~koethe/papers/abstracts/polarfilters.html">
    00530     <i>"Integrated Edge and Junction Detection with the Boundary Tensor"</i></a>,
    00531      in: ICCV 03, Proc. of 9th Intl. Conf. on Computer Vision, Nice 2003, vol. 1,
    00532      pp. 424-431, Los Alamitos: IEEE Computer Society, 2003
    00533
    00534     with the Laplacian of Gaussian as the underlying bandpass filter (see
    00535     \ref rieszTransformOfLOG()). The output image must have 3 bands which will hold the
    00536     tensor components in the order t11, t12 (== t21), t22. The function
    00537     \ref boundaryTensor1() with the same interface implements a variant of the
    00538     boundary tensor where the 0th-order Riesz transform has been dropped, so that the
    00539     tensor is no longer sensitive to blobs.
    00540
    00541     <b> Declarations:</b>
    00542
    00543     pass arguments explicitly:
    00544     \code
    00545     namespace vigra {
    00546         template <class SrcIterator, class SrcAccessor,
    00547                 class DestIterator, class DestAccessor>
    00548         void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
    00549                             DestIterator dupperleft, DestAccessor dest,
    00550                             double scale);
    00551     }
    00552     \endcode
    00553
    00554     use argument objects in conjunction with \ref ArgumentObjectFactories:
    00555     \code
    00556     namespace vigra {
    00557         template <class SrcIterator, class SrcAccessor,
    00558                 class DestIterator, class DestAccessor>
    00559         void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00560                             pair<DestIterator, DestAccessor> dest,
    00561                             double scale);
    00562     }
    00563     \endcode
    00564
    00565     <b> Usage:</b>
    00566
    00567     <b>\#include</b> "<a href="boundarytensor_8hxx-source.html">vigra/boundarytensor.hxx</a>"
    00568
    00569     \code
    00570     FImage img(w,h);
    00571     FVector3Image bt(w,h);
    00572     ...
    00573     boundaryTensor(srcImageRange(img), destImage(bt), 2.0);
    00574     \endcode
    00575
    00576 */
    00577 template <class SrcIterator, class SrcAccessor,
    00578           class DestIterator, class DestAccessor>
    00579 void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
    00580                     DestIterator dupperleft, DestAccessor dest,
    00581                     double scale)
    00582 {
    00583     vigra_precondition(dest.size(dupperleft) == 3,
    00584                        "boundaryTensor(): image for even output must have 3 bands.");
    00585     vigra_precondition(scale > 0.0,
    00586                        "boundaryTensor(): scale must be positive.");
    00587
    00588     detail::evenPolarFilters(supperleft, slowerright, src,
    00589                              dupperleft, dest, scale, false);
    00590     detail::oddPolarFilters(supperleft, slowerright, src,
    00591                              dupperleft, dest, scale, true);
    00592 }
    00593
    00594 template <class SrcIterator, class SrcAccessor,
    00595           class DestIterator, class DestAccessor>
    00596 inline
    00597 void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00598                     pair<DestIterator, DestAccessor> dest,
    00599                     double scale)
    00600 {
    00601     boundaryTensor(src.first, src.second, src.third,
    00602                    dest.first, dest.second, scale);
    00603 }
    00604
    00605 template <class SrcIterator, class SrcAccessor,
    00606           class DestIterator, class DestAccessor>
    00607 void boundaryTensor1(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
    00608                     DestIterator dupperleft, DestAccessor dest,
    00609                     double scale)
    00610 {
    00611     vigra_precondition(dest.size(dupperleft) == 3,
    00612                        "boundaryTensor1(): image for even output must have 3 bands.");
    00613     vigra_precondition(scale > 0.0,
    00614                        "boundaryTensor1(): scale must be positive.");
    00615
    00616     detail::evenPolarFilters(supperleft, slowerright, src,
    00617                              dupperleft, dest, scale, true);
    00618     detail::oddPolarFilters(supperleft, slowerright, src,
    00619                              dupperleft, dest, scale, true);
    00620 }
    00621
    00622 template <class SrcIterator, class SrcAccessor,
    00623           class DestIterator, class DestAccessor>
    00624 inline
    00625 void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00626                      pair<DestIterator, DestAccessor> dest,
    00627                      double scale)
    00628 {
    00629     boundaryTensor1(src.first, src.second, src.third,
    00630                     dest.first, dest.second, scale);
    00631 }
    00632
    00633 /********************************************************/
    00634 /*                                                      */
    00635 /*                    boundaryTensor3                   */
    00636 /*                                                      */
    00637 /********************************************************/
    00638
    00639 /*  Add 3rd order Riesz transform to boundary tensor
    00640     ??? Does not work -- bug or too coarse approximation for 3rd order ???
    00641 */
    00642 template <class SrcIterator, class SrcAccessor,
    00643           class DestIteratorEven, class DestAccessorEven,
    00644           class DestIteratorOdd, class DestAccessorOdd>
    00645 void boundaryTensor3(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor sa,
    00646                      DestIteratorEven dupperleft_even, DestAccessorEven even,
    00647                      DestIteratorOdd dupperleft_odd, DestAccessorOdd odd,
    00648                      double scale)
    00649 {
    00650     vigra_precondition(even.size(dupperleft_even) == 3,
    00651                        "boundaryTensor3(): image for even output must have 3 bands.");
    00652     vigra_precondition(odd.size(dupperleft_odd) == 3,
    00653                        "boundaryTensor3(): image for odd output must have 3 bands.");
    00654
    00655     detail::evenPolarFilters(supperleft, slowerright, sa,
    00656                              dupperleft_even, even, scale, false);
    00657
    00658     int w = slowerright.x - supperleft.x;
    00659     int h = slowerright.y - supperleft.y;
    00660
    00661     typedef typename
    00662        NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    00663     typedef BasicImage<TinyVector<TmpType, 4> > TmpImage;
    00664     TmpImage t1(w, h), t2(w, h);
    00665
    00666     detail::KernelArray k1, k3;
    00667     detail::initGaussianPolarFilters1(scale, k1);
    00668     detail::initGaussianPolarFilters3(scale, k3);
    00669
    00670     // calculate filter responses for odd filters
    00671     VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t1.accessor());
    00672     convolveImage(srcIterRange(supperleft, slowerright, sa),
    00673                   destImage(t1, tmpBand), k1[3], k1[0]);
    00674     tmpBand.setIndex(1);
    00675     convolveImage(srcIterRange(supperleft, slowerright, sa),
    00676                   destImage(t1, tmpBand), k1[1], k1[2]);
    00677     tmpBand.setIndex(2);
    00678     convolveImage(srcIterRange(supperleft, slowerright, sa),
    00679                   destImage(t1, tmpBand), k3[3], k3[0]);
    00680     tmpBand.setIndex(3);
    00681     convolveImage(srcIterRange(supperleft, slowerright, sa),
    00682                   destImage(t1, tmpBand), k3[1], k3[2]);
    00683
    00684     tmpBand.setIndex(0);
    00685     convolveImage(srcIterRange(supperleft, slowerright, sa),
    00686                   destImage(t2, tmpBand), k1[0], k1[3]);
    00687     tmpBand.setIndex(1);
    00688     convolveImage(srcIterRange(supperleft, slowerright, sa),
    00689                   destImage(t2, tmpBand), k1[2], k1[1]);
    00690     tmpBand.setIndex(2);
    00691     convolveImage(srcIterRange(supperleft, slowerright, sa),
    00692                   destImage(t2, tmpBand), k3[0], k3[3]);
    00693     tmpBand.setIndex(3);
    00694     convolveImage(srcIterRange(supperleft, slowerright, sa),
    00695                   destImage(t2, tmpBand), k3[2], k3[1]);
    00696
    00697     // create odd tensor from filter responses
    00698     typedef typename TmpImage::traverser TmpTraverser;
    00699     TmpTraverser tul1(t1.upperLeft());
    00700     TmpTraverser tlr1(t1.lowerRight());
    00701     TmpTraverser tul2(t2.upperLeft());
    00702     for(; tul1.y != tlr1.y; ++tul1.y, ++tul2.y, ++dupperleft_odd.y)
    00703     {
    00704         typename TmpTraverser::row_iterator tr1 = tul1.rowIterator();
    00705         typename TmpTraverser::row_iterator trend1 = tr1 + w;
    00706         typename TmpTraverser::row_iterator tr2 = tul2.rowIterator();
    00707         typename DestIteratorOdd::row_iterator o = dupperleft_odd.rowIterator();
    00708         for(; tr1 != trend1; ++tr1, ++tr2, ++o)
    00709         {
    00710             TmpType d11 =  (*tr1)[0] + (*tr1)[2];
    00711             TmpType d12 = -(*tr1)[1] - (*tr1)[3];
    00712             TmpType d31 =  (*tr2)[0] - (*tr2)[2];
    00713             TmpType d32 =  (*tr2)[1] - (*tr2)[3];
    00714             TmpType d111 = 0.75 * d11 + 0.25 * d31;
    00715             TmpType d112 = 0.25 * (d12 + d32);
    00716             TmpType d122 = 0.25 * (d11 - d31);
    00717             TmpType d222 = 0.75 * d12 - 0.25 * d32;
    00718             TmpType d2 = sq(d112);
    00719             TmpType d3 = sq(d122);
    00720
    00721             odd.setComponent(0.25 * (sq(d111) + 2.0*d2 + d3), o, 0);
    00722             odd.setComponent(0.25 * (d111*d112 + 2.0*d112*d122 + d122*d222), o, 1);
    00723             odd.setComponent(0.25 * (d2 + 2.0*d3 + sq(d222)), o, 2);
    00724         }
    00725     }
    00726 }
    00727
    00728 template <class SrcIterator, class SrcAccessor,
    00729           class DestIteratorEven, class DestAccessorEven,
    00730           class DestIteratorOdd, class DestAccessorOdd>
    00731 inline
    00732 void boundaryTensor3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00733                      pair<DestIteratorEven, DestAccessorEven> even,
    00734                      pair<DestIteratorOdd, DestAccessorOdd> odd,
    00735                      double scale)
    00736 {
    00737     boundaryTensor3(src.first, src.second, src.third,
    00738                     even.first, even.second, odd.first, odd.second, scale);
    00739 }
    00740
    00741 //@}
    00742
    00743 } // namespace vigra
    00744
    00745 #endif // VIGRA_BOUNDARYTENSOR_HXX

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)