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