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

vigra/orientedtensorfilters.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 2002-2004 by Ullrich Koethe                  */
00004 /*                                                                      */
00005 /*    This file is part of the VIGRA computer vision library.           */
00006 /*    The VIGRA Website is                                              */
00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00008 /*    Please direct questions, bug reports, and contributions to        */
00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00010 /*        vigra@informatik.uni-hamburg.de                               */
00011 /*                                                                      */
00012 /*    Permission is hereby granted, free of charge, to any person       */
00013 /*    obtaining a copy of this software and associated documentation    */
00014 /*    files (the "Software"), to deal in the Software without           */
00015 /*    restriction, including without limitation the rights to use,      */
00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00017 /*    sell copies of the Software, and to permit persons to whom the    */
00018 /*    Software is furnished to do so, subject to the following          */
00019 /*    conditions:                                                       */
00020 /*                                                                      */
00021 /*    The above copyright notice and this permission notice shall be    */
00022 /*    included in all copies or substantial portions of the             */
00023 /*    Software.                                                         */
00024 /*                                                                      */
00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */                
00033 /*                                                                      */
00034 /************************************************************************/
00035 
00036 #ifndef VIGRA_ORIENTEDTENSORFILTERS_HXX
00037 #define VIGRA_ORIENTEDTENSORFILTERS_HXX
00038 
00039 #include <cmath>
00040 #include "utilities.hxx"
00041 #include "initimage.hxx"
00042 #include "stdconvolution.hxx"
00043 
00044 namespace vigra {
00045 
00046 /** \addtogroup TensorImaging Tensor Image Processing
00047 */
00048 //@{
00049 
00050 /********************************************************/
00051 /*                                                      */
00052 /*                     hourGlassFilter                  */
00053 /*                                                      */
00054 /********************************************************/
00055 
00056 /** \brief Anisotropic tensor smoothing with the hourglass filter.
00057 
00058     This function implements anisotropic tensor smoothing by an
00059     hourglass-shaped filters as described in
00060     
00061     U. K&ouml;the: <a href="http://hci.iwr.uni-heidelberg.de/people/ukoethe/papers/index.php#cite_structureTensor">
00062     <i>"Edge and Junction Detection with an Improved Structure Tensor"</i></a>, 
00063      in: Proc. of 25th DAGM Symposium, Magdeburg 2003, Lecture Notes in Computer Science 2781, 
00064      pp. 25-32, Heidelberg: Springer, 2003
00065      
00066     It is closely related to the structure tensor (see \ref structureTensor()), but
00067     replaces the linear tensor smoothing with a smoothing along edges only. 
00068     Smoothing across edges is largely suppressed. This means that the
00069     image structure is preserved much better because nearby features
00070     such as parallel edges are not blended into each other. 
00071     
00072     The hourglass filter is typically applied to a gradient tensor, i.e. the 
00073     Euclidean product of the gradient with itself, which can be obtained by a
00074     gradient operator followed with \ref vectorToTensor(), see example below. 
00075     The hourglass shape of the filter can be interpreted as indicating the likely 
00076     continuations of a local edge element. The parameter <tt>sigma</tt> determines
00077     the radius of the hourglass (i.e. how far the influence of the edge element 
00078     reaches), and <tt>rho</tt> controls its opening angle (i.e. how narrow the 
00079     edge orientation os followed). Recommended values are <tt>sigma = 1.4</tt>
00080     (or, more generally, two to three times the scale of the gradient operator
00081     used in the first step), and <tt>rho = 0.4</tt> which corresponds to an 
00082     opening angle of 22.5 degrees to either side of the edge.
00083     
00084     <b> Declarations:</b>
00085 
00086     pass arguments explicitly:
00087     \code
00088     namespace vigra {
00089         template <class SrcIterator, class SrcAccessor,
00090                   class DestIterator, class DestAccessor>
00091         void hourGlassFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
00092                              DestIterator dul, DestAccessor dest,
00093                              double sigma, double rho);
00094     }
00095     \endcode
00096 
00097 
00098     use argument objects in conjunction with \ref ArgumentObjectFactories :
00099     \code
00100     namespace vigra {
00101         template <class SrcIterator, class SrcAccessor,
00102                   class DestIterator, class DestAccessor>
00103         inline
00104         void hourGlassFilter(triple<SrcIterator, SrcIterator, SrcAccessor> s,
00105                              pair<DestIterator, DestAccessor> d,
00106                              double sigma, double rho);
00107     }
00108     \endcode
00109 
00110     <b> Usage:</b>
00111 
00112     <b>\#include</b> <vigra/orientedtensorfilters.hxx>
00113 
00114     \code
00115     FImage img(w,h);
00116     FVector2Image gradient(w,h);
00117     FVector3Image tensor(w,h), smoothedTensor(w,h);
00118     
00119     gaussianGradient(srcImageRange(img), destImage(gradient), 1.0);
00120     vectorToTensor(srcImageRange(gradient), destImage(tensor));
00121     hourGlassFilter(srcImageRange(tensor), destImage(smoothedTensor), 2.0, 0.4);
00122     \endcode
00123     
00124     \see vectorToTensor()
00125 */
00126 doxygen_overloaded_function(template <...> void hourGlassFilter)
00127 
00128 template <class SrcIterator, class SrcAccessor,
00129           class DestIterator, class DestAccessor>
00130 void hourGlassFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
00131                      DestIterator dul, DestAccessor dest,
00132                      double sigma, double rho)
00133 {
00134     vigra_precondition(sigma >= 0.0 && rho >= 0.0,
00135                        "hourGlassFilter(): sigma and rho must be >= 0.0");
00136     vigra_precondition(src.size(sul) == 3,
00137                        "hourGlassFilter(): input image must have 3 bands.");
00138     vigra_precondition(dest.size(dul) == 3,
00139                        "hourGlassFilter(): output image must have 3 bands.");
00140 
00141     // TODO: normalization
00142 
00143     int w = slr.x - sul.x;
00144     int h = slr.y - sul.y;
00145 
00146     double radius = VIGRA_CSTD::floor(3.0*sigma + 0.5);
00147     double sigma2 = -0.5 / sigma / sigma;
00148     double rho2 = -0.5 / rho / rho;
00149     double norm = 1.0 / (2.0 * M_PI * sigma * sigma);
00150 
00151     initImage(dul, dul+Diff2D(w,h), dest, NumericTraits<typename DestAccessor::value_type>::zero());
00152 
00153     for(int y=0; y<h; ++y, ++sul.y, ++dul.y)
00154     {
00155         SrcIterator s = sul;
00156         DestIterator d = dul;
00157         for(int x=0; x<w; ++x, ++s.x, ++d.x)
00158         {
00159             double phi = 0.5 * VIGRA_CSTD::atan2(
00160                                      2.0*src.getComponent(s,1),
00161                                      (double)src.getComponent(s,0) - src.getComponent(s,2));
00162             double u = VIGRA_CSTD::sin(phi);
00163             double v = VIGRA_CSTD::cos(phi);
00164 
00165             double x0 = x - radius < 0 ? -x : -radius;
00166             double y0 = y - radius < 0 ? -y : -radius;
00167             double x1 = x + radius >= w ? w - x - 1 : radius;
00168             double y1 = y + radius >= h ? h - y - 1 : radius;
00169 
00170             DestIterator dwul = d + Diff2D((int)x0, (int)y0);
00171 
00172             for(double yy=y0; yy <= y1; ++yy, ++dwul.y)
00173             {
00174                 typename DestIterator::row_iterator dw = dwul.rowIterator();
00175                 for(double xx=x0; xx <= x1; ++xx, ++dw)
00176                 {
00177                     double r2 = xx*xx + yy*yy;
00178                     double p  = u*xx - v*yy;
00179                     double q  = v*xx + u*yy;
00180                     double kernel = (p == 0.0) ?
00181                                       (q == 0.0) ?
00182                                        norm :
00183                                        0.0 :
00184                                        norm * VIGRA_CSTD::exp(sigma2*r2 + rho2*q*q/p/p);
00185                     dest.set(dest(dw) + kernel*src(s), dw);
00186                 }
00187             }
00188         }
00189     }
00190 }
00191 
00192 template <class SrcIterator, class SrcAccessor,
00193           class DestIterator, class DestAccessor>
00194 inline
00195 void hourGlassFilter(triple<SrcIterator, SrcIterator, SrcAccessor> s,
00196                      pair<DestIterator, DestAccessor> d,
00197                      double sigma, double rho)
00198 {
00199     hourGlassFilter(s.first, s.second, s.third, d.first, d.second, sigma, rho);
00200 }
00201 
00202 /********************************************************/
00203 /*                                                      */
00204 /*                    ellipticGaussian                  */
00205 /*                                                      */
00206 /********************************************************/
00207 
00208 template <class SrcIterator, class SrcAccessor,
00209           class DestIterator, class DestAccessor>
00210 void ellipticGaussian(SrcIterator sul, SrcIterator slr, SrcAccessor src,
00211                       DestIterator dul, DestAccessor dest,
00212                       double sigmax, double sigmin)
00213 {
00214     vigra_precondition(sigmax >= sigmin && sigmin >= 0.0,
00215                        "ellipticGaussian(): "
00216                        "sigmax >= sigmin and sigmin >= 0.0 required");
00217 
00218     int w = slr.x - sul.x;
00219     int h = slr.y - sul.y;
00220 
00221     double radius = VIGRA_CSTD::floor(3.0*sigmax + 0.5);
00222     double sigmin2 = -0.5 / sigmin / sigmin;
00223     double norm = 1.0 / (2.0 * M_PI * sigmin * sigmax);
00224 
00225     initImage(dul, dul+Diff2D(w,h), dest, NumericTraits<typename DestAccessor::value_type>::zero());
00226 
00227     for(int y=0; y<h; ++y, ++sul.y, ++dul.y)
00228     {
00229         SrcIterator s = sul;
00230         DestIterator d = dul;
00231         for(int x=0; x<w; ++x, ++s.x, ++d.x)
00232         {
00233             typedef typename 
00234                NumericTraits<typename SrcAccessor::component_type>::RealPromote TmpType;
00235             TmpType d1 = src.getComponent(s,0) + src.getComponent(s,2);
00236             TmpType d2 = src.getComponent(s,0) - src.getComponent(s,2);
00237             TmpType d3 = 2.0 * src.getComponent(s,1);
00238             TmpType d4 = VIGRA_CSTD::sqrt(sq(d2) + sq(d3));
00239             TmpType excentricity = 1.0 - (d1 - d4) / (d1 + d4);
00240             double sigmax2 = -0.5 / sq((sigmax - sigmin)*excentricity + sigmin);
00241             
00242             double phi = 0.5 * VIGRA_CSTD::atan2(d3, d2);
00243             double u = VIGRA_CSTD::sin(phi);
00244             double v = VIGRA_CSTD::cos(phi);
00245 
00246             double x0 = x - radius < 0 ? -x : -radius;
00247             double y0 = y - radius < 0 ? -y : -radius;
00248             double x1 = x + radius >= w ? w - x - 1 : radius;
00249             double y1 = y + radius >= h ? h - y - 1 : radius;
00250 
00251             DestIterator dwul = d + Diff2D((int)x0, (int)y0);
00252 
00253             for(double yy=y0; yy <= y1; ++yy, ++dwul.y)
00254             {
00255                 typename DestIterator::row_iterator dw = dwul.rowIterator();
00256                 for(double xx=x0; xx <= x1; ++xx, ++dw)
00257                 {
00258                     double p  = u*xx - v*yy;
00259                     double q  = v*xx + u*yy;
00260                     double kernel = norm * VIGRA_CSTD::exp(sigmax2*p*p + sigmin2*q*q);
00261                     dest.set(dest(dw) + kernel*src(s), dw);
00262                 }
00263             }
00264         }
00265     }
00266 }
00267 
00268 template <class SrcIterator, class SrcAccessor,
00269           class DestIterator, class DestAccessor>
00270 inline
00271 void ellipticGaussian(triple<SrcIterator, SrcIterator, SrcAccessor> s,
00272                       pair<DestIterator, DestAccessor> d,
00273                       double sigmax, double sigmin)
00274 {
00275     ellipticGaussian(s.first, s.second, s.third, d.first, d.second, sigmax, sigmin);
00276 }
00277 
00278 /********************************************************/
00279 /*                                                      */
00280 /*         kernels for orientedTrigonometricFilter      */
00281 /*                                                      */
00282 /********************************************************/
00283 
00284 class FoerstnerKernelBase
00285 {
00286   public:
00287     typedef double ResultType;
00288     typedef double WeightType;
00289     typedef DVector2Image::value_type VectorType;
00290     typedef VectorType::value_type    ValueType;
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(ValueType(0), ValueType(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.8.0 (20 Sep 2011)