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