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

details vigra/orientedtensorfilters.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 2002-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.3, Aug 18 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 #ifndef VIGRA_ORIENTEDTENSORFILTERS_HXX
00024 #define VIGRA_ORIENTEDTENSORFILTERS_HXX
00025 
00026 #include <cmath>
00027 #include "vigra/utilities.hxx"
00028 #include "vigra/initimage.hxx"
00029 #include "vigra/stdconvolution.hxx"
00030 
00031 namespace vigra {
00032 
00033 /** \addtogroup TensorImaging Tensor Image Processing
00034 */
00035 //@{
00036 
00037 /********************************************************/
00038 /*                                                      */
00039 /*                     hourGlassFilter                  */
00040 /*                                                      */
00041 /********************************************************/
00042 
00043 /** \brief Anisotropic tensor smoothing with the hourglass filter.
00044 
00045     This function implements anisotropic tensor smoothing by an
00046     hourglass-shaped filters as described in
00047     
00048     U. Köthe: <a href="http://kogs-www.informatik.uni-hamburg.de/~koethe/papers/abstracts/structureTensor.html">
00049     <i>"Edge and Junction Detection with an Improved Structure Tensor"</i></a>, 
00050      in: Proc. of 25th DAGM Symposium, Magdeburg 2003, Lecture Notes in Computer Science 2781, 
00051      pp. 25-32, Heidelberg: Springer, 2003
00052      
00053     It is closely related to the structure tensor (see 
00054     \link CommonConvolutionFilters#structureTensor structureTensor\endlink()), but
00055     replaces the linear tensor smoothing with a smoothing along edges only. 
00056     Smoothing accross edges is largely suppressed. This means that the
00057     image structure is preserved much better because nearby features
00058     such as parallel edges are not blended into each other. 
00059     
00060     The hourglass filter is typically applied to a gradient tensor, i.e. the 
00061     Euclidean product of the gradient with itself, which can be obtained by a
00062     gradient operator followed with \ref vectorToTensor(), see example below. 
00063     The hourglass shape of the filter can be interpreted as indicating the likely 
00064     continuations of a local edge element. The parameter <tt>sigma</tt> determines
00065     the radius of the hourglass (i.e. how far the influence of the edge element 
00066     reaches), and <tt>rho</tt> controls its opening angle (i.e. how narrow the 
00067     edge orientation os followed). Recommended values are <tt>sigma = 1.4</tt>
00068     (or, more generally, two to three times the scale of the gradient operator
00069     used in the first step), and <tt>rho = 0.4</tt> which corresponds to an 
00070     opening angle of 22.5 degrees to either side of the edge.
00071     
00072     <b> Declarations:</b>
00073 
00074     pass arguments explicitly:
00075     \code
00076     namespace vigra {
00077         template <class SrcIterator, class SrcAccessor,
00078                   class DestIterator, class DestAccessor>
00079         void hourGlassFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
00080                              DestIterator dul, DestAccessor dest,
00081                              double sigma, double rho);
00082     }
00083     \endcode
00084 
00085 
00086     use argument objects in conjunction with \ref ArgumentObjectFactories:
00087     \code
00088     namespace vigra {
00089         template <class SrcIterator, class SrcAccessor,
00090                   class DestIterator, class DestAccessor>
00091         inline
00092         void hourGlassFilter(triple<SrcIterator, SrcIterator, SrcAccessor> s,
00093                              pair<DestIterator, DestAccessor> d,
00094                              double sigma, double rho);
00095     }
00096     \endcode
00097 
00098     <b> Usage:</b>
00099 
00100     <b>\#include</b> "<a href="orientedtensorfilters_8hxx-source.html">vigra/orientedtensorfilters.hxx</a>"
00101 
00102     \code
00103     FImage img(w,h);
00104     FVector2Image gradient(w,h);
00105     FVector3Image tensor(w,h), smoothedTensor(w,h);
00106     
00107     gaussianGradient(srcImageRange(img), destImage(gradient), 1.0);
00108     vectorToTensor(srcImageRange(gradient), destImage(tensor));
00109     hourGlassFilter(srcImageRange(tensor), destImage(smoothedTensor), 2.0, 0.4);
00110     \endcode
00111     
00112     \see vectorToTensor()
00113 */
00114 template <class SrcIterator, class SrcAccessor,
00115           class DestIterator, class DestAccessor>
00116 void hourGlassFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
00117                      DestIterator dul, DestAccessor dest,
00118                      double sigma, double rho)
00119 {
00120     vigra_precondition(sigma >= 0.0 && rho >= 0.0,
00121                        "hourGlassFilter(): sigma and rho must be >= 0.0");
00122     vigra_precondition(src.size(sul) == 3,
00123                        "hourGlassFilter(): input image must have 3 bands.");
00124     vigra_precondition(dest.size(dul) == 3,
00125                        "hourGlassFilter(): output image must have 3 bands.");
00126 
00127     // TODO: normalization
00128 
00129     int w = slr.x - sul.x;
00130     int h = slr.y - sul.y;
00131 
00132     double radius = VIGRA_CSTD::floor(3.0*sigma + 0.5);
00133     double sigma2 = -0.5 / sigma / sigma;
00134     double rho2 = -0.5 / rho / rho;
00135     double norm = 1.0 / (2.0 * M_PI * sigma * sigma);
00136 
00137     initImage(dul, dul+Diff2D(w,h), dest, NumericTraits<typename DestAccessor::value_type>::zero());
00138 
00139     for(int y=0; y<h; ++y, ++sul.y, ++dul.y)
00140     {
00141         SrcIterator s = sul;
00142         DestIterator d = dul;
00143         for(int x=0; x<w; ++x, ++s.x, ++d.x)
00144         {
00145             double phi = 0.5 * VIGRA_CSTD::atan2(
00146                                      2.0*src.getComponent(s,1),
00147                                      (double)src.getComponent(s,0) - src.getComponent(s,2));
00148             double u = -VIGRA_CSTD::sin(phi);
00149             double v = VIGRA_CSTD::cos(phi);
00150 
00151             double x0 = x - radius < 0 ? -x : -radius;
00152             double y0 = y - radius < 0 ? -y : -radius;
00153             double x1 = x + radius >= w ? w - x - 1 : radius;
00154             double y1 = y + radius >= h ? h - y - 1 : radius;
00155 
00156             DestIterator dwul = d + Diff2D((int)x0, (int)y0);
00157 
00158             for(double yy=y0; yy <= y1; ++yy, ++dwul.y)
00159             {
00160                 typename DestIterator::row_iterator dw = dwul.rowIterator();
00161                 for(double xx=x0; xx <= x1; ++xx, ++dw)
00162                 {
00163                     double r2 = xx*xx + yy*yy;
00164                     double p  = u*xx - v*yy;
00165                     double q  = v*xx + u*yy;
00166                     double kernel = (p == 0.0) ?
00167                                       (q == 0.0) ?
00168                                        norm :
00169                                        0.0 :
00170                                        norm * VIGRA_CSTD::exp(sigma2*r2 + rho2*q*q/p/p);
00171                     dest.set(dest(dw) + kernel*src(s), dw);
00172                 }
00173             }
00174         }
00175     }
00176 }
00177 
00178 template <class SrcIterator, class SrcAccessor,
00179           class DestIterator, class DestAccessor>
00180 inline
00181 void hourGlassFilter(triple<SrcIterator, SrcIterator, SrcAccessor> s,
00182                      pair<DestIterator, DestAccessor> d,
00183                      double sigma, double rho)
00184 {
00185     hourGlassFilter(s.first, s.second, s.third, d.first, d.second, sigma, rho);
00186 }
00187 
00188 /********************************************************/
00189 /*                                                      */
00190 /*                    ellipticGaussian                  */
00191 /*                                                      */
00192 /********************************************************/
00193 
00194 template <class SrcIterator, class SrcAccessor,
00195           class DestIterator, class DestAccessor>
00196 void ellipticGaussian(SrcIterator sul, SrcIterator slr, SrcAccessor src,
00197                       DestIterator dul, DestAccessor dest,
00198                       double sigmax, double sigmin)
00199 {
00200     vigra_precondition(sigmax >= sigmin && sigmin >= 0.0,
00201                        "ellipticGaussian(): "
00202                        "sigmax >= sigmin and sigmin >= 0.0 required");
00203 
00204     int w = slr.x - sul.x;
00205     int h = slr.y - sul.y;
00206 
00207     double radius = VIGRA_CSTD::floor(3.0*sigmax + 0.5);
00208     double sigmin2 = -0.5 / sigmin / sigmin;
00209     double norm = 1.0 / (2.0 * M_PI * sigmin * sigmax);
00210 
00211     initImage(dul, dul+Diff2D(w,h), dest, NumericTraits<typename DestAccessor::value_type>::zero());
00212 
00213     for(int y=0; y<h; ++y, ++sul.y, ++dul.y)
00214     {
00215         SrcIterator s = sul;
00216         DestIterator d = dul;
00217         for(int x=0; x<w; ++x, ++s.x, ++d.x)
00218         {
00219             typedef typename 
00220                NumericTraits<typename SrcAccessor::component_type>::RealPromote TmpType;
00221             TmpType d1 = src.getComponent(s,0) + src.getComponent(s,2);
00222             TmpType d2 = src.getComponent(s,0) - src.getComponent(s,2);
00223             TmpType d3 = 2.0 * src.getComponent(s,1);
00224             TmpType d4 = VIGRA_CSTD::sqrt(sq(d2) + sq(d3));
00225             TmpType excentricity = 1.0 - (d1 - d4) / (d1 + d4);
00226             double sigmax2 = -0.5 / sq((sigmax - sigmin)*excentricity + sigmin);
00227             
00228             double phi = 0.5 * VIGRA_CSTD::atan2(d3, d2);
00229             double u = -VIGRA_CSTD::sin(phi);
00230             double v = VIGRA_CSTD::cos(phi);
00231 
00232             double x0 = x - radius < 0 ? -x : -radius;
00233             double y0 = y - radius < 0 ? -y : -radius;
00234             double x1 = x + radius >= w ? w - x - 1 : radius;
00235             double y1 = y + radius >= h ? h - y - 1 : radius;
00236 
00237             DestIterator dwul = d + Diff2D((int)x0, (int)y0);
00238 
00239             for(double yy=y0; yy <= y1; ++yy, ++dwul.y)
00240             {
00241                 typename DestIterator::row_iterator dw = dwul.rowIterator();
00242                 for(double xx=x0; xx <= x1; ++xx, ++dw)
00243                 {
00244                     double p  = u*xx - v*yy;
00245                     double q  = v*xx + u*yy;
00246                     double kernel = norm * VIGRA_CSTD::exp(sigmax2*p*p + sigmin2*q*q);
00247                     dest.set(dest(dw) + kernel*src(s), dw);
00248                 }
00249             }
00250         }
00251     }
00252 }
00253 
00254 template <class SrcIterator, class SrcAccessor,
00255           class DestIterator, class DestAccessor>
00256 inline
00257 void ellipticGaussian(triple<SrcIterator, SrcIterator, SrcAccessor> s,
00258                       pair<DestIterator, DestAccessor> d,
00259                       double sigmax, double sigmin)
00260 {
00261     ellipticGaussian(s.first, s.second, s.third, d.first, d.second, sigmax, sigmin);
00262 }
00263 
00264 /********************************************************/
00265 /*                                                      */
00266 /*         kernels for orientedTrigonometricFilter      */
00267 /*                                                      */
00268 /********************************************************/
00269 
00270 class FoerstnerKernelBase
00271 {
00272   public:
00273     typedef double ResultType;
00274     typedef double WeightType;
00275     typedef DVector2Image::value_type VectorType;
00276   
00277     FoerstnerKernelBase(double scale, bool ringShaped = false)
00278     : radius_((int)(3.0*scale+0.5)),
00279       weights_(2*radius_+1, 2*radius_+1),
00280       vectors_(2*radius_+1, 2*radius_+1)
00281     {
00282         double norm = 1.0 / (2.0 * M_PI * scale * scale);
00283         double s2 = -0.5 / scale / scale;
00284         
00285         for(int y = -radius_; y <= radius_; ++y)
00286         {
00287             for(int x = -radius_; x <= radius_; ++x)
00288             {
00289                 double d2 = x*x + y*y;
00290                 double d = VIGRA_CSTD::sqrt(d2);
00291                 vectors_(x+radius_,y+radius_) = d != 0.0 ?
00292                                                   VectorType(x/d, -y/d) :
00293                                                   VectorType(0, 0);
00294                 weights_(x+radius_,y+radius_) = ringShaped ? 
00295                                        norm * d2 * VIGRA_CSTD::exp(d2 * s2) :
00296                                        norm * VIGRA_CSTD::exp(d2 * s2);
00297             }
00298         }
00299     }   
00300     
00301     ResultType operator()(int x, int y, VectorType const &) const
00302     {
00303         // isotropic filtering
00304         return weights_(radius_, radius_);
00305     }
00306 
00307     int radius_;
00308     DImage weights_;
00309     DVector2Image vectors_;
00310 };
00311 
00312 class FoerstnerRingKernelBase
00313 : public FoerstnerKernelBase
00314 {
00315   public:
00316     FoerstnerRingKernelBase(double scale)
00317     : FoerstnerKernelBase(scale, true)
00318     {}
00319 };
00320 
00321 class Cos2RingKernel
00322 : public FoerstnerKernelBase
00323 {
00324   public:
00325     typedef double ResultType;
00326     typedef double WeightType;
00327     typedef DVector2Image::value_type VectorType;
00328   
00329     Cos2RingKernel(double scale)
00330     : FoerstnerKernelBase(scale, true)
00331     {}
00332     
00333     ResultType operator()(int x, int y, VectorType const & v) const
00334     {
00335         if(x == 0 && y == 0)
00336             return weights_(radius_, radius_);
00337         double d = dot(vectors_(x+radius_, y+radius_), v);
00338         return d * d * weights_(x+radius_, y+radius_);
00339     }
00340 };
00341 
00342 class Cos2Kernel
00343 : public FoerstnerKernelBase
00344 {
00345   public:
00346     typedef double ResultType;
00347     typedef double WeightType;
00348     typedef DVector2Image::value_type VectorType;
00349   
00350     Cos2Kernel(double scale)
00351     : FoerstnerKernelBase(scale, false)
00352     {}
00353     
00354     ResultType operator()(int x, int y, VectorType const & v) const
00355     {
00356         if(x == 0 && y == 0)
00357             return weights_(radius_, radius_);
00358         double d = dot(vectors_(x+radius_, y+radius_), v);
00359         return d * d * weights_(x+radius_, y+radius_);
00360     }
00361 };
00362 
00363 class Sin2RingKernel
00364 : public FoerstnerKernelBase
00365 {
00366   public:
00367     typedef double ResultType;
00368     typedef double WeightType;
00369     typedef DVector2Image::value_type VectorType;
00370   
00371     Sin2RingKernel(double scale)
00372     : FoerstnerKernelBase(scale, true)
00373     {}
00374     
00375     ResultType operator()(int x, int y, VectorType const & v) const
00376     {
00377         if(x == 0 && y == 0)
00378             return weights_(radius_, radius_);
00379         double d = dot(vectors_(x+radius_, y+radius_), v);
00380         return (1.0 - d * d) * weights_(x+radius_, y+radius_);
00381     }
00382 };
00383 
00384 class Sin2Kernel
00385 : public FoerstnerKernelBase
00386 {
00387   public:
00388     typedef double ResultType;
00389     typedef double WeightType;
00390     typedef DVector2Image::value_type VectorType;
00391   
00392     Sin2Kernel(double scale)
00393     : FoerstnerKernelBase(scale, false)
00394     {}
00395     
00396     ResultType operator()(int x, int y, VectorType const & v) const
00397     {
00398         if(x == 0 && y == 0)
00399             return weights_(radius_, radius_);
00400         double d = dot(vectors_(x+radius_, y+radius_), v);
00401         return (1.0 - d * d) * weights_(x+radius_, y+radius_);
00402     }
00403 };
00404 
00405 class Sin6RingKernel
00406 : public FoerstnerKernelBase
00407 {
00408   public:
00409     typedef double ResultType;
00410     typedef double WeightType;
00411     typedef DVector2Image::value_type VectorType;
00412   
00413     Sin6RingKernel(double scale)
00414     : FoerstnerKernelBase(scale, true)
00415     {}
00416     
00417     ResultType operator()(int x, int y, VectorType const & v) const
00418     {
00419         if(x == 0 && y == 0)
00420             return weights_(radius_, radius_);
00421         double d = dot(vectors_(x+radius_, y+radius_), v);
00422         return VIGRA_CSTD::pow(1.0 - d * d, 3) * weights_(x+radius_, y+radius_);
00423     }
00424 };
00425 
00426 class Sin6Kernel
00427 : public FoerstnerKernelBase
00428 {
00429   public:
00430     typedef double ResultType;
00431     typedef double WeightType;
00432     typedef DVector2Image::value_type VectorType;
00433   
00434     Sin6Kernel(double scale)
00435     : FoerstnerKernelBase(scale, false)
00436     {}
00437     
00438     ResultType operator()(int x, int y, VectorType const & v) const
00439     {
00440         if(x == 0 && y == 0)
00441             return weights_(radius_, radius_);
00442         double d = dot(vectors_(x+radius_, y+radius_), v);
00443         return VIGRA_CSTD::pow(1.0 - d * d, 3) * weights_(x+radius_, y+radius_);
00444     }
00445 };
00446 
00447 class Cos6RingKernel
00448 : public FoerstnerKernelBase
00449 {
00450   public:
00451     typedef double ResultType;
00452     typedef double WeightType;
00453     typedef DVector2Image::value_type VectorType;
00454   
00455     Cos6RingKernel(double scale)
00456     : FoerstnerKernelBase(scale, true)
00457     {}
00458     
00459     ResultType operator()(int x, int y, VectorType const & v) const
00460     {
00461         if(x == 0 && y == 0)
00462             return weights_(radius_, radius_);
00463         double d = dot(vectors_(x+radius_, y+radius_), v);
00464         return (1.0 - VIGRA_CSTD::pow(1.0 - d * d, 3)) * weights_(x+radius_, y+radius_);
00465     }
00466 };
00467 
00468 class Cos6Kernel
00469 : public FoerstnerKernelBase
00470 {
00471   public:
00472     typedef double ResultType;
00473     typedef double WeightType;
00474     typedef DVector2Image::value_type VectorType;
00475   
00476     Cos6Kernel(double scale)
00477     : FoerstnerKernelBase(scale, false)
00478     {}
00479     
00480     ResultType operator()(int x, int y, VectorType const & v) const
00481     {
00482         if(x == 0 && y == 0)
00483             return weights_(radius_, radius_);
00484         double d = dot(vectors_(x+radius_, y+radius_), v);
00485         return (1.0 - VIGRA_CSTD::pow(1.0 - d * d, 3)) * weights_(x+radius_, y+radius_);
00486     }
00487 };
00488 
00489 /********************************************************/
00490 /*                                                      */
00491 /*              orientedTrigonometricFilter             */
00492 /*                                                      */
00493 /********************************************************/
00494 
00495 template <class SrcIterator, class SrcAccessor,
00496           class DestIterator, class DestAccessor,
00497           class Kernel>
00498 void orientedTrigonometricFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
00499                     DestIterator dul, DestAccessor dest,
00500                     Kernel const & kernel)
00501 {
00502     vigra_precondition(src.size(sul) == 2,
00503                        "orientedTrigonometricFilter(): input image must have 2 bands.");
00504     vigra_precondition(dest.size(dul) == 3,
00505                        "orientedTrigonometricFilter(): output image must have 3 bands.");
00506 
00507     int w = slr.x - sul.x;
00508     int h = slr.y - sul.y;
00509     int radius = kernel.radius_;
00510     
00511     typedef typename SrcAccessor::value_type VectorType;
00512     typedef typename DestAccessor::value_type TensorType;
00513 
00514     initImage(dul, dul+Diff2D(w,h), dest, NumericTraits<TensorType>::zero());
00515 
00516     for(int y=0; y<h; ++y, ++sul.y, ++dul.y)
00517     {
00518         SrcIterator s = sul;
00519         DestIterator d = dul;
00520         for(int x=0; x<w; ++x, ++s.x, ++d.x)
00521         {
00522             int x0 = x - radius < 0 ? -x : -radius;
00523             int y0 = y - radius < 0 ? -y : -radius;
00524             int x1 = x + radius >= w ? w - x - 1 : radius;
00525             int y1 = y + radius >= h ? h - y - 1 : radius;
00526 
00527             VectorType v(src(s));
00528             TensorType t(sq(v[0]), v[0]*v[1], sq(v[1]));
00529             double sqMag = t[0] + t[2];
00530             double mag = VIGRA_CSTD::sqrt(sqMag);
00531             if(mag != 0.0)
00532                 v /= mag;
00533             else
00534                 v *= 0.0;
00535             Diff2D dd;
00536             for(dd.y = y0; dd.y <= y1; ++dd.y)
00537             {
00538                 for(dd.x = x0; dd.x <= x1; ++dd.x)
00539                 {
00540                     dest.set(dest(d, dd) + kernel(dd.x, dd.y, v) * t, d, dd);
00541                 }
00542             }
00543         }
00544     }
00545 }
00546 
00547 template <class SrcIterator, class SrcAccessor,
00548           class DestIterator, class DestAccessor,
00549           class Kernel>
00550 inline
00551 void orientedTrigonometricFilter(triple<SrcIterator, SrcIterator, SrcAccessor> s,
00552                       pair<DestIterator, DestAccessor> d,
00553                       Kernel const & kernel)
00554 {
00555     orientedTrigonometricFilter(s.first, s.second, s.third, d.first, d.second, kernel);
00556 }
00557 
00558 //@}
00559 
00560 } // namespace vigra
00561 
00562 #endif /* VIGRA_ORIENTEDTENSORFILTERS_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.3 (18 Aug 2005)