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