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

| vigra/orientedtensorfilters.hxx |
--- | --- | ---

    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.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 #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.2 (27 Jan 2005)