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