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