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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)