vigra - vigra/boundarytensor.hxx Source File

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

vigra/boundarytensor.hxx |

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

© 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)