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

details vigra/boundarytensor.hxx VIGRA

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.3.2 (27 Jan 2005)