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

vigra/boundarytensor.hxx VIGRA

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.8.0 (20 Sep 2011)