[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-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 00039 #ifndef VIGRA_BOUNDARYTENSOR_HXX 00040 #define VIGRA_BOUNDARYTENSOR_HXX 00041 00042 #include <cmath> 00043 #include <functional> 00044 #include "utilities.hxx" 00045 #include "array_vector.hxx" 00046 #include "basicimage.hxx" 00047 #include "combineimages.hxx" 00048 #include "numerictraits.hxx" 00049 #include "convolution.hxx" 00050 00051 namespace vigra { 00052 00053 namespace detail { 00054 00055 /***********************************************************************/ 00056 00057 typedef ArrayVector<Kernel1D<double> > KernelArray; 00058 00059 void 00060 initGaussianPolarFilters1(double std_dev, KernelArray & k) 00061 { 00062 typedef KernelArray::value_type Kernel; 00063 typedef Kernel::iterator iterator; 00064 00065 vigra_precondition(std_dev >= 0.0, 00066 "initGaussianPolarFilter1(): " 00067 "Standard deviation must be >= 0."); 00068 00069 k.resize(4); 00070 00071 int radius = (int)(4.0*std_dev + 0.5); 00072 std_dev *= 1.08179074376; 00073 double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev; // norm 00074 double a = 0.558868151788 / VIGRA_CSTD::pow(std_dev, 5); 00075 double b = -2.04251639729 / VIGRA_CSTD::pow(std_dev, 3); 00076 double sigma22 = -0.5 / std_dev / std_dev; 00077 00078 00079 for(unsigned int i=0; i<k.size(); ++i) 00080 { 00081 k[i].initExplicitly(-radius, radius); 00082 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT); 00083 } 00084 00085 int ix; 00086 iterator c = k[0].center(); 00087 for(ix=-radius; ix<=radius; ++ix) 00088 { 00089 double x = (double)ix; 00090 c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x); 00091 } 00092 00093 c = k[1].center(); 00094 for(ix=-radius; ix<=radius; ++ix) 00095 { 00096 double x = (double)ix; 00097 c[ix] = f * x * VIGRA_CSTD::exp(sigma22 * x * x); 00098 } 00099 00100 c = k[2].center(); 00101 double b2 = b / 3.0; 00102 for(ix=-radius; ix<=radius; ++ix) 00103 { 00104 double x = (double)ix; 00105 c[ix] = f * (b2 + a * x * x) * VIGRA_CSTD::exp(sigma22 * x * x); 00106 } 00107 00108 c = k[3].center(); 00109 for(ix=-radius; ix<=radius; ++ix) 00110 { 00111 double x = (double)ix; 00112 c[ix] = f * x * (b + a * x * x) * VIGRA_CSTD::exp(sigma22 * x * x); 00113 } 00114 } 00115 00116 void 00117 initGaussianPolarFilters2(double std_dev, KernelArray & k) 00118 { 00119 typedef Kernel1D<double>::iterator iterator; 00120 00121 vigra_precondition(std_dev >= 0.0, 00122 "initGaussianPolarFilter2(): " 00123 "Standard deviation must be >= 0."); 00124 00125 k.resize(3); 00126 00127 int radius = (int)(4.0*std_dev + 0.5); 00128 double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev; // norm 00129 double sigma2 = std_dev*std_dev; 00130 double sigma22 = -0.5 / sigma2; 00131 00132 for(unsigned int i=0; i<k.size(); ++i) 00133 { 00134 k[i].initExplicitly(-radius, radius); 00135 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT); 00136 } 00137 00138 int ix; 00139 iterator c = k[0].center(); 00140 for(ix=-radius; ix<=radius; ++ix) 00141 { 00142 double x = (double)ix; 00143 c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x); 00144 } 00145 00146 c = k[1].center(); 00147 double f1 = f / sigma2; 00148 for(ix=-radius; ix<=radius; ++ix) 00149 { 00150 double x = (double)ix; 00151 c[ix] = f1 * x * VIGRA_CSTD::exp(sigma22 * x * x); 00152 } 00153 00154 c = k[2].center(); 00155 double f2 = f / (sigma2 * sigma2); 00156 for(ix=-radius; ix<=radius; ++ix) 00157 { 00158 double x = (double)ix; 00159 c[ix] = f2 * (x * x - sigma2) * VIGRA_CSTD::exp(sigma22 * x * x); 00160 } 00161 } 00162 00163 void 00164 initGaussianPolarFilters3(double std_dev, KernelArray & k) 00165 { 00166 typedef Kernel1D<double>::iterator iterator; 00167 00168 vigra_precondition(std_dev >= 0.0, 00169 "initGaussianPolarFilter3(): " 00170 "Standard deviation must be >= 0."); 00171 00172 k.resize(4); 00173 00174 int radius = (int)(4.0*std_dev + 0.5); 00175 std_dev *= 1.15470053838; 00176 double sigma22 = -0.5 / std_dev / std_dev; 00177 double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev; // norm 00178 double a = 0.883887052922 / VIGRA_CSTD::pow(std_dev, 5); 00179 00180 for(unsigned int i=0; i<k.size(); ++i) 00181 { 00182 k[i].initExplicitly(-radius, radius); 00183 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT); 00184 } 00185 00186 //double b = -1.3786348292 / VIGRA_CSTD::pow(std_dev, 3); 00187 00188 int ix; 00189 iterator c = k[0].center(); 00190 for(ix=-radius; ix<=radius; ++ix) 00191 { 00192 double x = (double)ix; 00193 c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x); 00194 } 00195 00196 c = k[1].center(); 00197 for(ix=-radius; ix<=radius; ++ix) 00198 { 00199 double x = (double)ix; 00200 c[ix] = f * x * VIGRA_CSTD::exp(sigma22 * x * x); 00201 } 00202 00203 c = k[2].center(); 00204 double a2 = 3.0 * a; 00205 for(ix=-radius; ix<=radius; ++ix) 00206 { 00207 double x = (double)ix; 00208 c[ix] = f * a2 * x * x * VIGRA_CSTD::exp(sigma22 * x * x); 00209 } 00210 00211 c = k[3].center(); 00212 for(ix=-radius; ix<=radius; ++ix) 00213 { 00214 double x = (double)ix; 00215 c[ix] = f * a * x * x * x * VIGRA_CSTD::exp(sigma22 * x * x); 00216 } 00217 } 00218 00219 template <class SrcIterator, class SrcAccessor, 00220 class DestIterator, class DestAccessor> 00221 void 00222 evenPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00223 DestIterator dupperleft, DestAccessor dest, 00224 double scale, bool noLaplacian) 00225 { 00226 vigra_precondition(dest.size(dupperleft) == 3, 00227 "evenPolarFilters(): image for even output must have 3 bands."); 00228 00229 int w = slowerright.x - supperleft.x; 00230 int h = slowerright.y - supperleft.y; 00231 00232 typedef typename 00233 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType; 00234 typedef BasicImage<TinyVector<TmpType, 3> > TmpImage; 00235 typedef typename TmpImage::traverser TmpTraverser; 00236 TmpImage t(w, h); 00237 00238 KernelArray k2; 00239 initGaussianPolarFilters2(scale, k2); 00240 00241 // calculate filter responses for even filters 00242 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor()); 00243 convolveImage(srcIterRange(supperleft, slowerright, src), 00244 destImage(t, tmpBand), k2[2], k2[0]); 00245 tmpBand.setIndex(1); 00246 convolveImage(srcIterRange(supperleft, slowerright, src), 00247 destImage(t, tmpBand), k2[1], k2[1]); 00248 tmpBand.setIndex(2); 00249 convolveImage(srcIterRange(supperleft, slowerright, src), 00250 destImage(t, tmpBand), k2[0], k2[2]); 00251 00252 // create even tensor from filter responses 00253 TmpTraverser tul(t.upperLeft()); 00254 TmpTraverser tlr(t.lowerRight()); 00255 for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y) 00256 { 00257 typename TmpTraverser::row_iterator tr = tul.rowIterator(); 00258 typename TmpTraverser::row_iterator trend = tr + w; 00259 typename DestIterator::row_iterator d = dupperleft.rowIterator(); 00260 if(noLaplacian) 00261 { 00262 for(; tr != trend; ++tr, ++d) 00263 { 00264 TmpType v = 0.5*sq((*tr)[0]-(*tr)[2]) + 2.0*sq((*tr)[1]); 00265 dest.setComponent(v, d, 0); 00266 dest.setComponent(0, d, 1); 00267 dest.setComponent(v, d, 2); 00268 } 00269 } 00270 else 00271 { 00272 for(; tr != trend; ++tr, ++d) 00273 { 00274 dest.setComponent(sq((*tr)[0]) + sq((*tr)[1]), d, 0); 00275 dest.setComponent(-(*tr)[1] * ((*tr)[0] + (*tr)[2]), d, 1); 00276 dest.setComponent(sq((*tr)[1]) + sq((*tr)[2]), d, 2); 00277 } 00278 } 00279 } 00280 } 00281 00282 template <class SrcIterator, class SrcAccessor, 00283 class DestIterator, class DestAccessor> 00284 void 00285 oddPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00286 DestIterator dupperleft, DestAccessor dest, 00287 double scale, bool addResult) 00288 { 00289 vigra_precondition(dest.size(dupperleft) == 3, 00290 "oddPolarFilters(): image for odd output must have 3 bands."); 00291 00292 int w = slowerright.x - supperleft.x; 00293 int h = slowerright.y - supperleft.y; 00294 00295 typedef typename 00296 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType; 00297 typedef BasicImage<TinyVector<TmpType, 4> > TmpImage; 00298 typedef typename TmpImage::traverser TmpTraverser; 00299 TmpImage t(w, h); 00300 00301 detail::KernelArray k1; 00302 detail::initGaussianPolarFilters1(scale, k1); 00303 00304 // calculate filter responses for odd filters 00305 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor()); 00306 convolveImage(srcIterRange(supperleft, slowerright, src), 00307 destImage(t, tmpBand), k1[3], k1[0]); 00308 tmpBand.setIndex(1); 00309 convolveImage(srcIterRange(supperleft, slowerright, src), 00310 destImage(t, tmpBand), k1[2], k1[1]); 00311 tmpBand.setIndex(2); 00312 convolveImage(srcIterRange(supperleft, slowerright, src), 00313 destImage(t, tmpBand), k1[1], k1[2]); 00314 tmpBand.setIndex(3); 00315 convolveImage(srcIterRange(supperleft, slowerright, src), 00316 destImage(t, tmpBand), k1[0], k1[3]); 00317 00318 // create odd tensor from filter responses 00319 TmpTraverser tul(t.upperLeft()); 00320 TmpTraverser tlr(t.lowerRight()); 00321 for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y) 00322 { 00323 typename TmpTraverser::row_iterator tr = tul.rowIterator(); 00324 typename TmpTraverser::row_iterator trend = tr + w; 00325 typename DestIterator::row_iterator d = dupperleft.rowIterator(); 00326 if(addResult) 00327 { 00328 for(; tr != trend; ++tr, ++d) 00329 { 00330 TmpType d0 = (*tr)[0] + (*tr)[2]; 00331 TmpType d1 = -(*tr)[1] - (*tr)[3]; 00332 00333 dest.setComponent(dest.getComponent(d, 0) + sq(d0), d, 0); 00334 dest.setComponent(dest.getComponent(d, 1) + d0 * d1, d, 1); 00335 dest.setComponent(dest.getComponent(d, 2) + sq(d1), d, 2); 00336 } 00337 } 00338 else 00339 { 00340 for(; tr != trend; ++tr, ++d) 00341 { 00342 TmpType d0 = (*tr)[0] + (*tr)[2]; 00343 TmpType d1 = -(*tr)[1] - (*tr)[3]; 00344 00345 dest.setComponent(sq(d0), d, 0); 00346 dest.setComponent(d0 * d1, d, 1); 00347 dest.setComponent(sq(d1), d, 2); 00348 } 00349 } 00350 } 00351 } 00352 00353 } // namespace detail 00354 00355 /** \addtogroup CommonConvolutionFilters Common Filters 00356 */ 00357 //@{ 00358 00359 /********************************************************/ 00360 /* */ 00361 /* rieszTransformOfLOG */ 00362 /* */ 00363 /********************************************************/ 00364 00365 /** \brief Calculate Riesz transforms of the Laplacian of Gaussian. 00366 00367 The Riesz transforms of the Laplacian of Gaussian have the following transfer 00368 functions (defined in a polar coordinate representation of the frequency domain): 00369 00370 \f[ 00371 F_{\sigma}(r, \phi)=(i \cos \phi)^n (i \sin \phi)^m r^2 e^{-r^2 \sigma^2 / 2} 00372 \f] 00373 00374 where <i>n</i> = <tt>xorder</tt> and <i>m</i> = <tt>yorder</tt> determine th e 00375 order of the transform, and <tt>sigma > 0</tt> is the scale of the Laplacian 00376 of Gaussian. This function computes a good spatial domain approximation of 00377 these transforms for <tt>xorder + yorder <= 2</tt>. The filter responses may be used 00378 to calculate the monogenic signal or the boundary tensor. 00379 00380 <b> Declarations:</b> 00381 00382 pass arguments explicitly: 00383 \code 00384 namespace vigra { 00385 template <class SrcIterator, class SrcAccessor, 00386 class DestIterator, class DestAccessor> 00387 void rieszTransformOfLOG(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00388 DestIterator dupperleft, DestAccessor dest, 00389 double scale, unsigned int xorder, unsigned int yorder); 00390 } 00391 \endcode 00392 00393 00394 use argument objects in conjunction with \ref ArgumentObjectFactories : 00395 \code 00396 namespace vigra { 00397 template <class SrcIterator, class SrcAccessor, 00398 class DestIterator, class DestAccessor> 00399 void rieszTransformOfLOG(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00400 pair<DestIterator, DestAccessor> dest, 00401 double scale, unsigned int xorder, unsigned int yorder); 00402 } 00403 \endcode 00404 00405 <b> Usage:</b> 00406 00407 <b>\#include</b> <<a href="boundarytensor_8hxx-source.html">vigra/boundarytensor.hxx</a>> 00408 00409 \code 00410 FImage impulse(17,17), res(17, 17); 00411 impulse(8,8) = 1.0; 00412 00413 // calculate the impulse response of the first order Riesz transform in x-direction 00414 rieszTransformOfLOG(srcImageRange(impulse), destImage(res), 2.0, 1, 0); 00415 \endcode 00416 00417 */ 00418 doxygen_overloaded_function(template <...> void rieszTransformOfLOG) 00419 00420 template <class SrcIterator, class SrcAccessor, 00421 class DestIterator, class DestAccessor> 00422 void rieszTransformOfLOG(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00423 DestIterator dupperleft, DestAccessor dest, 00424 double scale, unsigned int xorder, unsigned int yorder) 00425 { 00426 unsigned int order = xorder + yorder; 00427 00428 vigra_precondition(order <= 2, 00429 "rieszTransformOfLOG(): can only compute Riesz transforms up to order 2."); 00430 vigra_precondition(scale > 0.0, 00431 "rieszTransformOfLOG(): scale must be positive."); 00432 00433 int w = slowerright.x - supperleft.x; 00434 int h = slowerright.y - supperleft.y; 00435 00436 typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType; 00437 typedef BasicImage<TmpType> TmpImage; 00438 00439 switch(order) 00440 { 00441 case 0: 00442 { 00443 detail::KernelArray k2; 00444 detail::initGaussianPolarFilters2(scale, k2); 00445 00446 TmpImage t1(w, h), t2(w, h); 00447 00448 convolveImage(srcIterRange(supperleft, slowerright, src), 00449 destImage(t1), k2[2], k2[0]); 00450 convolveImage(srcIterRange(supperleft, slowerright, src), 00451 destImage(t2), k2[0], k2[2]); 00452 combineTwoImages(srcImageRange(t1), srcImage(t2), 00453 destIter(dupperleft, dest), std::plus<TmpType>()); 00454 break; 00455 } 00456 case 1: 00457 { 00458 detail::KernelArray k1; 00459 detail::initGaussianPolarFilters1(scale, k1); 00460 00461 TmpImage t1(w, h), t2(w, h); 00462 00463 if(xorder == 1) 00464 { 00465 convolveImage(srcIterRange(supperleft, slowerright, src), 00466 destImage(t1), k1[3], k1[0]); 00467 convolveImage(srcIterRange(supperleft, slowerright, src), 00468 destImage(t2), k1[1], k1[2]); 00469 } 00470 else 00471 { 00472 convolveImage(srcIterRange(supperleft, slowerright, src), 00473 destImage(t1), k1[0], k1[3]); 00474 convolveImage(srcIterRange(supperleft, slowerright, src), 00475 destImage(t2), k1[2], k1[1]); 00476 } 00477 combineTwoImages(srcImageRange(t1), srcImage(t2), 00478 destIter(dupperleft, dest), std::plus<TmpType>()); 00479 break; 00480 } 00481 case 2: 00482 { 00483 detail::KernelArray k2; 00484 detail::initGaussianPolarFilters2(scale, k2); 00485 00486 convolveImage(srcIterRange(supperleft, slowerright, src), 00487 destIter(dupperleft, dest), k2[xorder], k2[yorder]); 00488 break; 00489 } 00490 /* for test purposes only: compute 3rd order polar filters */ 00491 case 3: 00492 { 00493 detail::KernelArray k3; 00494 detail::initGaussianPolarFilters3(scale, k3); 00495 00496 TmpImage t1(w, h), t2(w, h); 00497 00498 if(xorder == 3) 00499 { 00500 convolveImage(srcIterRange(supperleft, slowerright, src), 00501 destImage(t1), k3[3], k3[0]); 00502 convolveImage(srcIterRange(supperleft, slowerright, src), 00503 destImage(t2), k3[1], k3[2]); 00504 } 00505 else 00506 { 00507 convolveImage(srcIterRange(supperleft, slowerright, src), 00508 destImage(t1), k3[0], k3[3]); 00509 convolveImage(srcIterRange(supperleft, slowerright, src), 00510 destImage(t2), k3[2], k3[1]); 00511 } 00512 combineTwoImages(srcImageRange(t1), srcImage(t2), 00513 destIter(dupperleft, dest), std::minus<TmpType>()); 00514 break; 00515 } 00516 } 00517 } 00518 00519 template <class SrcIterator, class SrcAccessor, 00520 class DestIterator, class DestAccessor> 00521 inline 00522 void rieszTransformOfLOG(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00523 pair<DestIterator, DestAccessor> dest, 00524 double scale, unsigned int xorder, unsigned int yorder) 00525 { 00526 rieszTransformOfLOG(src.first, src.second, src.third, dest.first, dest.second, 00527 scale, xorder, yorder); 00528 } 00529 //@} 00530 00531 /** \addtogroup TensorImaging Tensor Image Processing 00532 */ 00533 //@{ 00534 00535 /********************************************************/ 00536 /* */ 00537 /* boundaryTensor */ 00538 /* */ 00539 /********************************************************/ 00540 00541 /** \brief Calculate the boundary tensor for a scalar valued image. 00542 00543 These functions calculate a spatial domain approximation of 00544 the boundary tensor as described in 00545 00546 U. Köthe: <a href="http://kogs-www.informatik.uni-hamburg.de/~koethe/papers/abstracts/polarfilters.html"> 00547 <i>"Integrated Edge and Junction Detection with the Boundary Tensor"</i></a>, 00548 in: ICCV 03, Proc. of 9th Intl. Conf. on Computer Vision, Nice 2003, vol. 1, 00549 pp. 424-431, Los Alamitos: IEEE Computer Society, 2003 00550 00551 with the Laplacian of Gaussian as the underlying bandpass filter (see 00552 \ref rieszTransformOfLOG()). The output image must have 3 bands which will hold the 00553 tensor components in the order t11, t12 (== t21), t22. The function 00554 \ref boundaryTensor1() with the same interface implements a variant of the 00555 boundary tensor where the 0th-order Riesz transform has been dropped, so that the 00556 tensor is no longer sensitive to blobs. 00557 00558 <b> Declarations:</b> 00559 00560 pass arguments explicitly: 00561 \code 00562 namespace vigra { 00563 template <class SrcIterator, class SrcAccessor, 00564 class DestIterator, class DestAccessor> 00565 void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00566 DestIterator dupperleft, DestAccessor dest, 00567 double scale); 00568 } 00569 \endcode 00570 00571 use argument objects in conjunction with \ref ArgumentObjectFactories : 00572 \code 00573 namespace vigra { 00574 template <class SrcIterator, class SrcAccessor, 00575 class DestIterator, class DestAccessor> 00576 void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00577 pair<DestIterator, DestAccessor> dest, 00578 double scale); 00579 } 00580 \endcode 00581 00582 <b> Usage:</b> 00583 00584 <b>\#include</b> <<a href="boundarytensor_8hxx-source.html">vigra/boundarytensor.hxx</a>> 00585 00586 \code 00587 FImage img(w,h); 00588 FVector3Image bt(w,h); 00589 ... 00590 boundaryTensor(srcImageRange(img), destImage(bt), 2.0); 00591 \endcode 00592 00593 */ 00594 doxygen_overloaded_function(template <...> void boundaryTensor) 00595 00596 template <class SrcIterator, class SrcAccessor, 00597 class DestIterator, class DestAccessor> 00598 void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00599 DestIterator dupperleft, DestAccessor dest, 00600 double scale) 00601 { 00602 vigra_precondition(dest.size(dupperleft) == 3, 00603 "boundaryTensor(): image for even output must have 3 bands."); 00604 vigra_precondition(scale > 0.0, 00605 "boundaryTensor(): scale must be positive."); 00606 00607 detail::evenPolarFilters(supperleft, slowerright, src, 00608 dupperleft, dest, scale, false); 00609 detail::oddPolarFilters(supperleft, slowerright, src, 00610 dupperleft, dest, scale, true); 00611 } 00612 00613 template <class SrcIterator, class SrcAccessor, 00614 class DestIterator, class DestAccessor> 00615 inline 00616 void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00617 pair<DestIterator, DestAccessor> dest, 00618 double scale) 00619 { 00620 boundaryTensor(src.first, src.second, src.third, 00621 dest.first, dest.second, scale); 00622 } 00623 00624 /** \brief Boundary tensor variant. 00625 00626 This function implements a variant of the boundary tensor where the 00627 0th-order Riesz transform has been dropped, so that the tensor is no 00628 longer sensitive to blobs. See \ref boundaryTensor() for more detailed 00629 documentation. 00630 00631 <b> Declarations:</b> 00632 00633 <b>\#include</b> <<a href="boundarytensor_8hxx-source.html">vigra/boundarytensor.hxx</a>> 00634 00635 pass arguments explicitly: 00636 \code 00637 namespace vigra { 00638 template <class SrcIterator, class SrcAccessor, 00639 class DestIterator, class DestAccessor> 00640 void boundaryTensor1(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00641 DestIterator dupperleft, DestAccessor dest, 00642 double scale); 00643 } 00644 \endcode 00645 00646 use argument objects in conjunction with \ref ArgumentObjectFactories : 00647 \code 00648 namespace vigra { 00649 template <class SrcIterator, class SrcAccessor, 00650 class DestIterator, class DestAccessor> 00651 void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00652 pair<DestIterator, DestAccessor> dest, 00653 double scale); 00654 } 00655 \endcode 00656 */ 00657 doxygen_overloaded_function(template <...> void boundaryTensor1) 00658 00659 template <class SrcIterator, class SrcAccessor, 00660 class DestIterator, class DestAccessor> 00661 void boundaryTensor1(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00662 DestIterator dupperleft, DestAccessor dest, 00663 double scale) 00664 { 00665 vigra_precondition(dest.size(dupperleft) == 3, 00666 "boundaryTensor1(): image for even output must have 3 bands."); 00667 vigra_precondition(scale > 0.0, 00668 "boundaryTensor1(): scale must be positive."); 00669 00670 detail::evenPolarFilters(supperleft, slowerright, src, 00671 dupperleft, dest, scale, true); 00672 detail::oddPolarFilters(supperleft, slowerright, src, 00673 dupperleft, dest, scale, true); 00674 } 00675 00676 template <class SrcIterator, class SrcAccessor, 00677 class DestIterator, class DestAccessor> 00678 inline 00679 void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00680 pair<DestIterator, DestAccessor> dest, 00681 double scale) 00682 { 00683 boundaryTensor1(src.first, src.second, src.third, 00684 dest.first, dest.second, scale); 00685 } 00686 00687 /********************************************************/ 00688 /* */ 00689 /* boundaryTensor3 */ 00690 /* */ 00691 /********************************************************/ 00692 00693 /* Add 3rd order Riesz transform to boundary tensor 00694 ??? Does not work -- bug or too coarse approximation for 3rd order ??? 00695 */ 00696 template <class SrcIterator, class SrcAccessor, 00697 class DestIteratorEven, class DestAccessorEven, 00698 class DestIteratorOdd, class DestAccessorOdd> 00699 void boundaryTensor3(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor sa, 00700 DestIteratorEven dupperleft_even, DestAccessorEven even, 00701 DestIteratorOdd dupperleft_odd, DestAccessorOdd odd, 00702 double scale) 00703 { 00704 vigra_precondition(even.size(dupperleft_even) == 3, 00705 "boundaryTensor3(): image for even output must have 3 bands."); 00706 vigra_precondition(odd.size(dupperleft_odd) == 3, 00707 "boundaryTensor3(): image for odd output must have 3 bands."); 00708 00709 detail::evenPolarFilters(supperleft, slowerright, sa, 00710 dupperleft_even, even, scale, false); 00711 00712 int w = slowerright.x - supperleft.x; 00713 int h = slowerright.y - supperleft.y; 00714 00715 typedef typename 00716 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType; 00717 typedef BasicImage<TinyVector<TmpType, 4> > TmpImage; 00718 TmpImage t1(w, h), t2(w, h); 00719 00720 detail::KernelArray k1, k3; 00721 detail::initGaussianPolarFilters1(scale, k1); 00722 detail::initGaussianPolarFilters3(scale, k3); 00723 00724 // calculate filter responses for odd filters 00725 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t1.accessor()); 00726 convolveImage(srcIterRange(supperleft, slowerright, sa), 00727 destImage(t1, tmpBand), k1[3], k1[0]); 00728 tmpBand.setIndex(1); 00729 convolveImage(srcIterRange(supperleft, slowerright, sa), 00730 destImage(t1, tmpBand), k1[1], k1[2]); 00731 tmpBand.setIndex(2); 00732 convolveImage(srcIterRange(supperleft, slowerright, sa), 00733 destImage(t1, tmpBand), k3[3], k3[0]); 00734 tmpBand.setIndex(3); 00735 convolveImage(srcIterRange(supperleft, slowerright, sa), 00736 destImage(t1, tmpBand), k3[1], k3[2]); 00737 00738 tmpBand.setIndex(0); 00739 convolveImage(srcIterRange(supperleft, slowerright, sa), 00740 destImage(t2, tmpBand), k1[0], k1[3]); 00741 tmpBand.setIndex(1); 00742 convolveImage(srcIterRange(supperleft, slowerright, sa), 00743 destImage(t2, tmpBand), k1[2], k1[1]); 00744 tmpBand.setIndex(2); 00745 convolveImage(srcIterRange(supperleft, slowerright, sa), 00746 destImage(t2, tmpBand), k3[0], k3[3]); 00747 tmpBand.setIndex(3); 00748 convolveImage(srcIterRange(supperleft, slowerright, sa), 00749 destImage(t2, tmpBand), k3[2], k3[1]); 00750 00751 // create odd tensor from filter responses 00752 typedef typename TmpImage::traverser TmpTraverser; 00753 TmpTraverser tul1(t1.upperLeft()); 00754 TmpTraverser tlr1(t1.lowerRight()); 00755 TmpTraverser tul2(t2.upperLeft()); 00756 for(; tul1.y != tlr1.y; ++tul1.y, ++tul2.y, ++dupperleft_odd.y) 00757 { 00758 typename TmpTraverser::row_iterator tr1 = tul1.rowIterator(); 00759 typename TmpTraverser::row_iterator trend1 = tr1 + w; 00760 typename TmpTraverser::row_iterator tr2 = tul2.rowIterator(); 00761 typename DestIteratorOdd::row_iterator o = dupperleft_odd.rowIterator(); 00762 for(; tr1 != trend1; ++tr1, ++tr2, ++o) 00763 { 00764 TmpType d11 = (*tr1)[0] + (*tr1)[2]; 00765 TmpType d12 = -(*tr1)[1] - (*tr1)[3]; 00766 TmpType d31 = (*tr2)[0] - (*tr2)[2]; 00767 TmpType d32 = (*tr2)[1] - (*tr2)[3]; 00768 TmpType d111 = 0.75 * d11 + 0.25 * d31; 00769 TmpType d112 = 0.25 * (d12 + d32); 00770 TmpType d122 = 0.25 * (d11 - d31); 00771 TmpType d222 = 0.75 * d12 - 0.25 * d32; 00772 TmpType d2 = sq(d112); 00773 TmpType d3 = sq(d122); 00774 00775 odd.setComponent(0.25 * (sq(d111) + 2.0*d2 + d3), o, 0); 00776 odd.setComponent(0.25 * (d111*d112 + 2.0*d112*d122 + d122*d222), o, 1); 00777 odd.setComponent(0.25 * (d2 + 2.0*d3 + sq(d222)), o, 2); 00778 } 00779 } 00780 } 00781 00782 template <class SrcIterator, class SrcAccessor, 00783 class DestIteratorEven, class DestAccessorEven, 00784 class DestIteratorOdd, class DestAccessorOdd> 00785 inline 00786 void boundaryTensor3(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00787 pair<DestIteratorEven, DestAccessorEven> even, 00788 pair<DestIteratorOdd, DestAccessorOdd> odd, 00789 double scale) 00790 { 00791 boundaryTensor3(src.first, src.second, src.third, 00792 even.first, even.second, odd.first, odd.second, scale); 00793 } 00794 00795 //@} 00796 00797 } // namespace vigra 00798 00799 #endif // VIGRA_BOUNDARYTENSOR_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|