[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | 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.6.0, Aug 13 2008 )                                    */
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 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00012 /*        vigra@informatik.uni-hamburg.de                               */
00013 /*                                                                      */
00014 /*    Permission is hereby granted, free of charge, to any person       */
00015 /*    obtaining a copy of this software and associated documentation    */
00016 /*    files (the "Software"), to deal in the Software without           */
00017 /*    restriction, including without limitation the rights to use,      */
00018 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00019 /*    sell copies of the Software, and to permit persons to whom the    */
00020 /*    Software is furnished to do so, subject to the following          */
00021 /*    conditions:                                                       */
00022 /*                                                                      */
00023 /*    The above copyright notice and this permission notice shall be    */
00024 /*    included in all copies or substantial portions of the             */
00025 /*    Software.                                                         */
00026 /*                                                                      */
00027 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00028 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00029 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00030 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00031 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00032 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00033 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00034 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00035 /*                                                                      */
00036 /************************************************************************/
00037 
00038 
00039 #ifndef VIGRA_BOUNDARYTENSOR_HXX
00040 #define VIGRA_BOUNDARYTENSOR_HXX
00041 
00042 #include <cmath>
00043 #include <functional>
00044 #include "utilities.hxx"
00045 #include "array_vector.hxx"
00046 #include "basicimage.hxx"
00047 #include "combineimages.hxx"
00048 #include "numerictraits.hxx"
00049 #include "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 > 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 <= 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 doxygen_overloaded_function(template <...> void rieszTransformOfLOG)
00419 
00420 template <class SrcIterator, class SrcAccessor,
00421           class DestIterator, class DestAccessor>
00422 void rieszTransformOfLOG(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
00423                          DestIterator dupperleft, DestAccessor dest,
00424                          double scale, unsigned int xorder, unsigned int yorder)
00425 {
00426     unsigned int order = xorder + yorder;
00427 
00428     vigra_precondition(order <= 2,
00429             "rieszTransformOfLOG(): can only compute Riesz transforms up to order 2.");
00430     vigra_precondition(scale > 0.0,
00431             "rieszTransformOfLOG(): scale must be positive.");
00432 
00433     int w = slowerright.x - supperleft.x;
00434     int h = slowerright.y - supperleft.y;
00435 
00436     typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
00437     typedef BasicImage<TmpType> TmpImage;
00438 
00439     switch(order)
00440     {
00441         case 0:
00442         {
00443             detail::KernelArray k2;
00444             detail::initGaussianPolarFilters2(scale, k2);
00445 
00446             TmpImage t1(w, h), t2(w, h);
00447 
00448             convolveImage(srcIterRange(supperleft, slowerright, src),
00449                           destImage(t1), k2[2], k2[0]);
00450             convolveImage(srcIterRange(supperleft, slowerright, src),
00451                           destImage(t2), k2[0], k2[2]);
00452             combineTwoImages(srcImageRange(t1), srcImage(t2),
00453                              destIter(dupperleft, dest), std::plus<TmpType>());
00454             break;
00455         }
00456         case 1:
00457         {
00458             detail::KernelArray k1;
00459             detail::initGaussianPolarFilters1(scale, k1);
00460 
00461             TmpImage t1(w, h), t2(w, h);
00462 
00463             if(xorder == 1)
00464             {
00465                 convolveImage(srcIterRange(supperleft, slowerright, src),
00466                             destImage(t1), k1[3], k1[0]);
00467                 convolveImage(srcIterRange(supperleft, slowerright, src),
00468                             destImage(t2), k1[1], k1[2]);
00469             }
00470             else
00471             {
00472                 convolveImage(srcIterRange(supperleft, slowerright, src),
00473                             destImage(t1), k1[0], k1[3]);
00474                 convolveImage(srcIterRange(supperleft, slowerright, src),
00475                             destImage(t2), k1[2], k1[1]);
00476             }
00477             combineTwoImages(srcImageRange(t1), srcImage(t2),
00478                              destIter(dupperleft, dest), std::plus<TmpType>());
00479             break;
00480         }
00481         case 2:
00482         {
00483             detail::KernelArray k2;
00484             detail::initGaussianPolarFilters2(scale, k2);
00485 
00486             convolveImage(srcIterRange(supperleft, slowerright, src),
00487                           destIter(dupperleft, dest), k2[xorder], k2[yorder]);
00488             break;
00489         }
00490         /* for test purposes only: compute 3rd order polar filters */
00491         case 3:
00492         {
00493             detail::KernelArray k3;
00494             detail::initGaussianPolarFilters3(scale, k3);
00495 
00496             TmpImage t1(w, h), t2(w, h);
00497 
00498             if(xorder == 3)
00499             {
00500                 convolveImage(srcIterRange(supperleft, slowerright, src),
00501                             destImage(t1), k3[3], k3[0]);
00502                 convolveImage(srcIterRange(supperleft, slowerright, src),
00503                             destImage(t2), k3[1], k3[2]);
00504             }
00505             else
00506             {
00507                 convolveImage(srcIterRange(supperleft, slowerright, src),
00508                             destImage(t1), k3[0], k3[3]);
00509                 convolveImage(srcIterRange(supperleft, slowerright, src),
00510                             destImage(t2), k3[2], k3[1]);
00511             }
00512             combineTwoImages(srcImageRange(t1), srcImage(t2),
00513                              destIter(dupperleft, dest), std::minus<TmpType>());
00514             break;
00515         }
00516     }
00517 }
00518 
00519 template <class SrcIterator, class SrcAccessor,
00520           class DestIterator, class DestAccessor>
00521 inline
00522 void rieszTransformOfLOG(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00523                          pair<DestIterator, DestAccessor> dest,
00524                          double scale, unsigned int xorder, unsigned int yorder)
00525 {
00526     rieszTransformOfLOG(src.first, src.second, src.third, dest.first, dest.second,
00527                         scale, xorder, yorder);
00528 }
00529 //@}
00530 
00531 /** \addtogroup TensorImaging Tensor Image Processing
00532 */
00533 //@{
00534 
00535 /********************************************************/
00536 /*                                                      */
00537 /*                     boundaryTensor                   */
00538 /*                                                      */
00539 /********************************************************/
00540 
00541 /** \brief Calculate the boundary tensor for a scalar valued image.
00542 
00543     These functions calculate a spatial domain approximation of
00544     the boundary tensor as described in
00545 
00546     U. K&ouml;the: <a href="http://kogs-www.informatik.uni-hamburg.de/~koethe/papers/abstracts/polarfilters.html">
00547     <i>"Integrated Edge and Junction Detection with the Boundary Tensor"</i></a>,
00548      in: ICCV 03, Proc. of 9th Intl. Conf. on Computer Vision, Nice 2003, vol. 1,
00549      pp. 424-431, Los Alamitos: IEEE Computer Society, 2003
00550 
00551     with the Laplacian of Gaussian as the underlying bandpass filter (see
00552     \ref rieszTransformOfLOG()). The output image must have 3 bands which will hold the
00553     tensor components in the order t11, t12 (== t21), t22. The function
00554     \ref boundaryTensor1() with the same interface implements a variant of the
00555     boundary tensor where the 0th-order Riesz transform has been dropped, so that the
00556     tensor is no longer sensitive to blobs.
00557 
00558     <b> Declarations:</b>
00559 
00560     pass arguments explicitly:
00561     \code
00562     namespace vigra {
00563         template <class SrcIterator, class SrcAccessor,
00564                   class DestIterator, class DestAccessor>
00565         void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
00566                             DestIterator dupperleft, DestAccessor dest,
00567                             double scale);
00568     }
00569     \endcode
00570 
00571     use argument objects in conjunction with \ref ArgumentObjectFactories :
00572     \code
00573     namespace vigra {
00574         template <class SrcIterator, class SrcAccessor,
00575                   class DestIterator, class DestAccessor>
00576         void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00577                             pair<DestIterator, DestAccessor> dest,
00578                             double scale);
00579     }
00580     \endcode
00581 
00582     <b> Usage:</b>
00583 
00584     <b>\#include</b> <<a href="boundarytensor_8hxx-source.html">vigra/boundarytensor.hxx</a>>
00585 
00586     \code
00587     FImage img(w,h);
00588     FVector3Image bt(w,h);
00589     ...
00590     boundaryTensor(srcImageRange(img), destImage(bt), 2.0);
00591     \endcode
00592 
00593 */
00594 doxygen_overloaded_function(template <...> void boundaryTensor)
00595 
00596 template <class SrcIterator, class SrcAccessor,
00597           class DestIterator, class DestAccessor>
00598 void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
00599                     DestIterator dupperleft, DestAccessor dest,
00600                     double scale)
00601 {
00602     vigra_precondition(dest.size(dupperleft) == 3,
00603                        "boundaryTensor(): image for even output must have 3 bands.");
00604     vigra_precondition(scale > 0.0,
00605                        "boundaryTensor(): scale must be positive.");
00606 
00607     detail::evenPolarFilters(supperleft, slowerright, src,
00608                              dupperleft, dest, scale, false);
00609     detail::oddPolarFilters(supperleft, slowerright, src,
00610                              dupperleft, dest, scale, true);
00611 }
00612 
00613 template <class SrcIterator, class SrcAccessor,
00614           class DestIterator, class DestAccessor>
00615 inline
00616 void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00617                     pair<DestIterator, DestAccessor> dest,
00618                     double scale)
00619 {
00620     boundaryTensor(src.first, src.second, src.third,
00621                    dest.first, dest.second, scale);
00622 }
00623 
00624 /** \brief Boundary tensor variant.
00625 
00626     This function implements a variant of the boundary tensor where the 
00627     0th-order Riesz transform has been dropped, so that the tensor is no 
00628     longer sensitive to blobs. See \ref boundaryTensor() for more detailed 
00629     documentation.
00630 
00631     <b> Declarations:</b>
00632 
00633     <b>\#include</b> <<a href="boundarytensor_8hxx-source.html">vigra/boundarytensor.hxx</a>>
00634 
00635     pass arguments explicitly:
00636     \code
00637     namespace vigra {
00638         template <class SrcIterator, class SrcAccessor,
00639                   class DestIterator, class DestAccessor>
00640         void boundaryTensor1(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
00641                              DestIterator dupperleft, DestAccessor dest,
00642                              double scale);
00643     }
00644     \endcode
00645 
00646     use argument objects in conjunction with \ref ArgumentObjectFactories :
00647     \code
00648     namespace vigra {
00649         template <class SrcIterator, class SrcAccessor,
00650                   class DestIterator, class DestAccessor>
00651         void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00652                              pair<DestIterator, DestAccessor> dest,
00653                              double scale);
00654     }
00655     \endcode
00656 */
00657 doxygen_overloaded_function(template <...> void boundaryTensor1)
00658 
00659 template <class SrcIterator, class SrcAccessor,
00660           class DestIterator, class DestAccessor>
00661 void boundaryTensor1(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
00662                     DestIterator dupperleft, DestAccessor dest,
00663                     double scale)
00664 {
00665     vigra_precondition(dest.size(dupperleft) == 3,
00666                        "boundaryTensor1(): image for even output must have 3 bands.");
00667     vigra_precondition(scale > 0.0,
00668                        "boundaryTensor1(): scale must be positive.");
00669 
00670     detail::evenPolarFilters(supperleft, slowerright, src,
00671                              dupperleft, dest, scale, true);
00672     detail::oddPolarFilters(supperleft, slowerright, src,
00673                              dupperleft, dest, scale, true);
00674 }
00675 
00676 template <class SrcIterator, class SrcAccessor,
00677           class DestIterator, class DestAccessor>
00678 inline
00679 void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00680                      pair<DestIterator, DestAccessor> dest,
00681                      double scale)
00682 {
00683     boundaryTensor1(src.first, src.second, src.third,
00684                     dest.first, dest.second, scale);
00685 }
00686 
00687 /********************************************************/
00688 /*                                                      */
00689 /*                    boundaryTensor3                   */
00690 /*                                                      */
00691 /********************************************************/
00692 
00693 /*  Add 3rd order Riesz transform to boundary tensor
00694     ??? Does not work -- bug or too coarse approximation for 3rd order ???
00695 */
00696 template <class SrcIterator, class SrcAccessor,
00697           class DestIteratorEven, class DestAccessorEven,
00698           class DestIteratorOdd, class DestAccessorOdd>
00699 void boundaryTensor3(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor sa,
00700                      DestIteratorEven dupperleft_even, DestAccessorEven even,
00701                      DestIteratorOdd dupperleft_odd, DestAccessorOdd odd,
00702                      double scale)
00703 {
00704     vigra_precondition(even.size(dupperleft_even) == 3,
00705                        "boundaryTensor3(): image for even output must have 3 bands.");
00706     vigra_precondition(odd.size(dupperleft_odd) == 3,
00707                        "boundaryTensor3(): image for odd output must have 3 bands.");
00708 
00709     detail::evenPolarFilters(supperleft, slowerright, sa,
00710                              dupperleft_even, even, scale, false);
00711 
00712     int w = slowerright.x - supperleft.x;
00713     int h = slowerright.y - supperleft.y;
00714 
00715     typedef typename
00716        NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
00717     typedef BasicImage<TinyVector<TmpType, 4> > TmpImage;
00718     TmpImage t1(w, h), t2(w, h);
00719 
00720     detail::KernelArray k1, k3;
00721     detail::initGaussianPolarFilters1(scale, k1);
00722     detail::initGaussianPolarFilters3(scale, k3);
00723 
00724     // calculate filter responses for odd filters
00725     VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t1.accessor());
00726     convolveImage(srcIterRange(supperleft, slowerright, sa),
00727                   destImage(t1, tmpBand), k1[3], k1[0]);
00728     tmpBand.setIndex(1);
00729     convolveImage(srcIterRange(supperleft, slowerright, sa),
00730                   destImage(t1, tmpBand), k1[1], k1[2]);
00731     tmpBand.setIndex(2);
00732     convolveImage(srcIterRange(supperleft, slowerright, sa),
00733                   destImage(t1, tmpBand), k3[3], k3[0]);
00734     tmpBand.setIndex(3);
00735     convolveImage(srcIterRange(supperleft, slowerright, sa),
00736                   destImage(t1, tmpBand), k3[1], k3[2]);
00737 
00738     tmpBand.setIndex(0);
00739     convolveImage(srcIterRange(supperleft, slowerright, sa),
00740                   destImage(t2, tmpBand), k1[0], k1[3]);
00741     tmpBand.setIndex(1);
00742     convolveImage(srcIterRange(supperleft, slowerright, sa),
00743                   destImage(t2, tmpBand), k1[2], k1[1]);
00744     tmpBand.setIndex(2);
00745     convolveImage(srcIterRange(supperleft, slowerright, sa),
00746                   destImage(t2, tmpBand), k3[0], k3[3]);
00747     tmpBand.setIndex(3);
00748     convolveImage(srcIterRange(supperleft, slowerright, sa),
00749                   destImage(t2, tmpBand), k3[2], k3[1]);
00750 
00751     // create odd tensor from filter responses
00752     typedef typename TmpImage::traverser TmpTraverser;
00753     TmpTraverser tul1(t1.upperLeft());
00754     TmpTraverser tlr1(t1.lowerRight());
00755     TmpTraverser tul2(t2.upperLeft());
00756     for(; tul1.y != tlr1.y; ++tul1.y, ++tul2.y, ++dupperleft_odd.y)
00757     {
00758         typename TmpTraverser::row_iterator tr1 = tul1.rowIterator();
00759         typename TmpTraverser::row_iterator trend1 = tr1 + w;
00760         typename TmpTraverser::row_iterator tr2 = tul2.rowIterator();
00761         typename DestIteratorOdd::row_iterator o = dupperleft_odd.rowIterator();
00762         for(; tr1 != trend1; ++tr1, ++tr2, ++o)
00763         {
00764             TmpType d11 =  (*tr1)[0] + (*tr1)[2];
00765             TmpType d12 = -(*tr1)[1] - (*tr1)[3];
00766             TmpType d31 =  (*tr2)[0] - (*tr2)[2];
00767             TmpType d32 =  (*tr2)[1] - (*tr2)[3];
00768             TmpType d111 = 0.75 * d11 + 0.25 * d31;
00769             TmpType d112 = 0.25 * (d12 + d32);
00770             TmpType d122 = 0.25 * (d11 - d31);
00771             TmpType d222 = 0.75 * d12 - 0.25 * d32;
00772             TmpType d2 = sq(d112);
00773             TmpType d3 = sq(d122);
00774 
00775             odd.setComponent(0.25 * (sq(d111) + 2.0*d2 + d3), o, 0);
00776             odd.setComponent(0.25 * (d111*d112 + 2.0*d112*d122 + d122*d222), o, 1);
00777             odd.setComponent(0.25 * (d2 + 2.0*d3 + sq(d222)), o, 2);
00778         }
00779     }
00780 }
00781 
00782 template <class SrcIterator, class SrcAccessor,
00783           class DestIteratorEven, class DestAccessorEven,
00784           class DestIteratorOdd, class DestAccessorOdd>
00785 inline
00786 void boundaryTensor3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00787                      pair<DestIteratorEven, DestAccessorEven> even,
00788                      pair<DestIteratorOdd, DestAccessorOdd> odd,
00789                      double scale)
00790 {
00791     boundaryTensor3(src.first, src.second, src.third,
00792                     even.first, even.second, odd.first, odd.second, scale);
00793 }
00794 
00795 //@}
00796 
00797 } // namespace vigra
00798 
00799 #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.6.0 (13 Aug 2008)