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