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