[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
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ö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) |
html generated using doxygen and Python
|