[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
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ö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 accross 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> <<a href="orientedtensorfilters_8hxx-source.html">vigra/orientedtensorfilters.hxx</a>> 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) |
html generated using doxygen and Python
|