[ 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.4.0, Dec 21 2005 ) */ 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 /* koethe@informatik.uni-hamburg.de or */ 00012 /* vigra@kogs1.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 "vigra/utilities.hxx" 00045 #include "vigra/array_vector.hxx" 00046 #include "vigra/basicimage.hxx" 00047 #include "vigra/combineimages.hxx" 00048 #include "vigra/numerictraits.hxx" 00049 #include "vigra/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 template <class SrcIterator, class SrcAccessor, 00419 class DestIterator, class DestAccessor> 00420 void rieszTransformOfLOG(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00421 DestIterator dupperleft, DestAccessor dest, 00422 double scale, unsigned int xorder, unsigned int yorder) 00423 { 00424 unsigned int order = xorder + yorder; 00425 00426 vigra_precondition(order <= 2, 00427 "rieszTransformOfLOG(): can only compute Riesz transforms up to order 2."); 00428 vigra_precondition(scale > 0.0, 00429 "rieszTransformOfLOG(): scale must be positive."); 00430 00431 int w = slowerright.x - supperleft.x; 00432 int h = slowerright.y - supperleft.y; 00433 00434 typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType; 00435 typedef BasicImage<TmpType> TmpImage; 00436 00437 switch(order) 00438 { 00439 case 0: 00440 { 00441 detail::KernelArray k2; 00442 detail::initGaussianPolarFilters2(scale, k2); 00443 00444 TmpImage t1(w, h), t2(w, h); 00445 00446 convolveImage(srcIterRange(supperleft, slowerright, src), 00447 destImage(t1), k2[2], k2[0]); 00448 convolveImage(srcIterRange(supperleft, slowerright, src), 00449 destImage(t2), k2[0], k2[2]); 00450 combineTwoImages(srcImageRange(t1), srcImage(t2), 00451 destIter(dupperleft, dest), std::plus<TmpType>()); 00452 break; 00453 } 00454 case 1: 00455 { 00456 detail::KernelArray k1; 00457 detail::initGaussianPolarFilters1(scale, k1); 00458 00459 TmpImage t1(w, h), t2(w, h); 00460 00461 if(xorder == 1) 00462 { 00463 convolveImage(srcIterRange(supperleft, slowerright, src), 00464 destImage(t1), k1[3], k1[0]); 00465 convolveImage(srcIterRange(supperleft, slowerright, src), 00466 destImage(t2), k1[1], k1[2]); 00467 } 00468 else 00469 { 00470 convolveImage(srcIterRange(supperleft, slowerright, src), 00471 destImage(t1), k1[0], k1[3]); 00472 convolveImage(srcIterRange(supperleft, slowerright, src), 00473 destImage(t2), k1[2], k1[1]); 00474 } 00475 combineTwoImages(srcImageRange(t1), srcImage(t2), 00476 destIter(dupperleft, dest), std::plus<TmpType>()); 00477 break; 00478 } 00479 case 2: 00480 { 00481 detail::KernelArray k2; 00482 detail::initGaussianPolarFilters2(scale, k2); 00483 00484 convolveImage(srcIterRange(supperleft, slowerright, src), 00485 destIter(dupperleft, dest), k2[xorder], k2[yorder]); 00486 break; 00487 } 00488 /* for test purposes only: compute 3rd order polar filters */ 00489 case 3: 00490 { 00491 detail::KernelArray k3; 00492 detail::initGaussianPolarFilters3(scale, k3); 00493 00494 TmpImage t1(w, h), t2(w, h); 00495 00496 if(xorder == 3) 00497 { 00498 convolveImage(srcIterRange(supperleft, slowerright, src), 00499 destImage(t1), k3[3], k3[0]); 00500 convolveImage(srcIterRange(supperleft, slowerright, src), 00501 destImage(t2), k3[1], k3[2]); 00502 } 00503 else 00504 { 00505 convolveImage(srcIterRange(supperleft, slowerright, src), 00506 destImage(t1), k3[0], k3[3]); 00507 convolveImage(srcIterRange(supperleft, slowerright, src), 00508 destImage(t2), k3[2], k3[1]); 00509 } 00510 combineTwoImages(srcImageRange(t1), srcImage(t2), 00511 destIter(dupperleft, dest), std::minus<TmpType>()); 00512 break; 00513 } 00514 } 00515 } 00516 00517 template <class SrcIterator, class SrcAccessor, 00518 class DestIterator, class DestAccessor> 00519 inline 00520 void rieszTransformOfLOG(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00521 pair<DestIterator, DestAccessor> dest, 00522 double scale, unsigned int xorder, unsigned int yorder) 00523 { 00524 rieszTransformOfLOG(src.first, src.second, src.third, dest.first, dest.second, 00525 scale, xorder, yorder); 00526 } 00527 //@} 00528 00529 /** \addtogroup TensorImaging Tensor Image Processing 00530 */ 00531 //@{ 00532 00533 /********************************************************/ 00534 /* */ 00535 /* boundaryTensor */ 00536 /* */ 00537 /********************************************************/ 00538 00539 /** \brief Calculate the boundary tensor for a scalar valued image. 00540 00541 These functions calculates a spatial domain approximation of 00542 the boundary tensor as described in 00543 00544 U. K÷the: <a href="http://kogs-www.informatik.uni-hamburg.de/~koethe/papers/abstracts/polarfilters.html"> 00545 <i>"Integrated Edge and Junction Detection with the Boundary Tensor"</i></a>, 00546 in: ICCV 03, Proc. of 9th Intl. Conf. on Computer Vision, Nice 2003, vol. 1, 00547 pp. 424-431, Los Alamitos: IEEE Computer Society, 2003 00548 00549 with the Laplacian of Gaussian as the underlying bandpass filter (see 00550 \ref rieszTransformOfLOG()). The output image must have 3 bands which will hold the 00551 tensor components in the order t11, t12 (== t21), t22. The function 00552 \ref boundaryTensor1() with the same interface implements a variant of the 00553 boundary tensor where the 0th-order Riesz transform has been dropped, so that the 00554 tensor is no longer sensitive to blobs. 00555 00556 <b> Declarations:</b> 00557 00558 pass arguments explicitly: 00559 \code 00560 namespace vigra { 00561 template <class SrcIterator, class SrcAccessor, 00562 class DestIterator, class DestAccessor> 00563 void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00564 DestIterator dupperleft, DestAccessor dest, 00565 double scale); 00566 } 00567 \endcode 00568 00569 use argument objects in conjunction with \ref ArgumentObjectFactories: 00570 \code 00571 namespace vigra { 00572 template <class SrcIterator, class SrcAccessor, 00573 class DestIterator, class DestAccessor> 00574 void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00575 pair<DestIterator, DestAccessor> dest, 00576 double scale); 00577 } 00578 \endcode 00579 00580 <b> Usage:</b> 00581 00582 <b>\#include</b> "<a href="boundarytensor_8hxx-source.html">vigra/boundarytensor.hxx</a>" 00583 00584 \code 00585 FImage img(w,h); 00586 FVector3Image bt(w,h); 00587 ... 00588 boundaryTensor(srcImageRange(img), destImage(bt), 2.0); 00589 \endcode 00590 00591 */ 00592 template <class SrcIterator, class SrcAccessor, 00593 class DestIterator, class DestAccessor> 00594 void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00595 DestIterator dupperleft, DestAccessor dest, 00596 double scale) 00597 { 00598 vigra_precondition(dest.size(dupperleft) == 3, 00599 "boundaryTensor(): image for even output must have 3 bands."); 00600 vigra_precondition(scale > 0.0, 00601 "boundaryTensor(): scale must be positive."); 00602 00603 detail::evenPolarFilters(supperleft, slowerright, src, 00604 dupperleft, dest, scale, false); 00605 detail::oddPolarFilters(supperleft, slowerright, src, 00606 dupperleft, dest, scale, true); 00607 } 00608 00609 template <class SrcIterator, class SrcAccessor, 00610 class DestIterator, class DestAccessor> 00611 inline 00612 void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00613 pair<DestIterator, DestAccessor> dest, 00614 double scale) 00615 { 00616 boundaryTensor(src.first, src.second, src.third, 00617 dest.first, dest.second, scale); 00618 } 00619 00620 template <class SrcIterator, class SrcAccessor, 00621 class DestIterator, class DestAccessor> 00622 void boundaryTensor1(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src, 00623 DestIterator dupperleft, DestAccessor dest, 00624 double scale) 00625 { 00626 vigra_precondition(dest.size(dupperleft) == 3, 00627 "boundaryTensor1(): image for even output must have 3 bands."); 00628 vigra_precondition(scale > 0.0, 00629 "boundaryTensor1(): scale must be positive."); 00630 00631 detail::evenPolarFilters(supperleft, slowerright, src, 00632 dupperleft, dest, scale, true); 00633 detail::oddPolarFilters(supperleft, slowerright, src, 00634 dupperleft, dest, scale, true); 00635 } 00636 00637 template <class SrcIterator, class SrcAccessor, 00638 class DestIterator, class DestAccessor> 00639 inline 00640 void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00641 pair<DestIterator, DestAccessor> dest, 00642 double scale) 00643 { 00644 boundaryTensor1(src.first, src.second, src.third, 00645 dest.first, dest.second, scale); 00646 } 00647 00648 /********************************************************/ 00649 /* */ 00650 /* boundaryTensor3 */ 00651 /* */ 00652 /********************************************************/ 00653 00654 /* Add 3rd order Riesz transform to boundary tensor 00655 ??? Does not work -- bug or too coarse approximation for 3rd order ??? 00656 */ 00657 template <class SrcIterator, class SrcAccessor, 00658 class DestIteratorEven, class DestAccessorEven, 00659 class DestIteratorOdd, class DestAccessorOdd> 00660 void boundaryTensor3(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor sa, 00661 DestIteratorEven dupperleft_even, DestAccessorEven even, 00662 DestIteratorOdd dupperleft_odd, DestAccessorOdd odd, 00663 double scale) 00664 { 00665 vigra_precondition(even.size(dupperleft_even) == 3, 00666 "boundaryTensor3(): image for even output must have 3 bands."); 00667 vigra_precondition(odd.size(dupperleft_odd) == 3, 00668 "boundaryTensor3(): image for odd output must have 3 bands."); 00669 00670 detail::evenPolarFilters(supperleft, slowerright, sa, 00671 dupperleft_even, even, scale, false); 00672 00673 int w = slowerright.x - supperleft.x; 00674 int h = slowerright.y - supperleft.y; 00675 00676 typedef typename 00677 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType; 00678 typedef BasicImage<TinyVector<TmpType, 4> > TmpImage; 00679 TmpImage t1(w, h), t2(w, h); 00680 00681 detail::KernelArray k1, k3; 00682 detail::initGaussianPolarFilters1(scale, k1); 00683 detail::initGaussianPolarFilters3(scale, k3); 00684 00685 // calculate filter responses for odd filters 00686 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t1.accessor()); 00687 convolveImage(srcIterRange(supperleft, slowerright, sa), 00688 destImage(t1, tmpBand), k1[3], k1[0]); 00689 tmpBand.setIndex(1); 00690 convolveImage(srcIterRange(supperleft, slowerright, sa), 00691 destImage(t1, tmpBand), k1[1], k1[2]); 00692 tmpBand.setIndex(2); 00693 convolveImage(srcIterRange(supperleft, slowerright, sa), 00694 destImage(t1, tmpBand), k3[3], k3[0]); 00695 tmpBand.setIndex(3); 00696 convolveImage(srcIterRange(supperleft, slowerright, sa), 00697 destImage(t1, tmpBand), k3[1], k3[2]); 00698 00699 tmpBand.setIndex(0); 00700 convolveImage(srcIterRange(supperleft, slowerright, sa), 00701 destImage(t2, tmpBand), k1[0], k1[3]); 00702 tmpBand.setIndex(1); 00703 convolveImage(srcIterRange(supperleft, slowerright, sa), 00704 destImage(t2, tmpBand), k1[2], k1[1]); 00705 tmpBand.setIndex(2); 00706 convolveImage(srcIterRange(supperleft, slowerright, sa), 00707 destImage(t2, tmpBand), k3[0], k3[3]); 00708 tmpBand.setIndex(3); 00709 convolveImage(srcIterRange(supperleft, slowerright, sa), 00710 destImage(t2, tmpBand), k3[2], k3[1]); 00711 00712 // create odd tensor from filter responses 00713 typedef typename TmpImage::traverser TmpTraverser; 00714 TmpTraverser tul1(t1.upperLeft()); 00715 TmpTraverser tlr1(t1.lowerRight()); 00716 TmpTraverser tul2(t2.upperLeft()); 00717 for(; tul1.y != tlr1.y; ++tul1.y, ++tul2.y, ++dupperleft_odd.y) 00718 { 00719 typename TmpTraverser::row_iterator tr1 = tul1.rowIterator(); 00720 typename TmpTraverser::row_iterator trend1 = tr1 + w; 00721 typename TmpTraverser::row_iterator tr2 = tul2.rowIterator(); 00722 typename DestIteratorOdd::row_iterator o = dupperleft_odd.rowIterator(); 00723 for(; tr1 != trend1; ++tr1, ++tr2, ++o) 00724 { 00725 TmpType d11 = (*tr1)[0] + (*tr1)[2]; 00726 TmpType d12 = -(*tr1)[1] - (*tr1)[3]; 00727 TmpType d31 = (*tr2)[0] - (*tr2)[2]; 00728 TmpType d32 = (*tr2)[1] - (*tr2)[3]; 00729 TmpType d111 = 0.75 * d11 + 0.25 * d31; 00730 TmpType d112 = 0.25 * (d12 + d32); 00731 TmpType d122 = 0.25 * (d11 - d31); 00732 TmpType d222 = 0.75 * d12 - 0.25 * d32; 00733 TmpType d2 = sq(d112); 00734 TmpType d3 = sq(d122); 00735 00736 odd.setComponent(0.25 * (sq(d111) + 2.0*d2 + d3), o, 0); 00737 odd.setComponent(0.25 * (d111*d112 + 2.0*d112*d122 + d122*d222), o, 1); 00738 odd.setComponent(0.25 * (d2 + 2.0*d3 + sq(d222)), o, 2); 00739 } 00740 } 00741 } 00742 00743 template <class SrcIterator, class SrcAccessor, 00744 class DestIteratorEven, class DestAccessorEven, 00745 class DestIteratorOdd, class DestAccessorOdd> 00746 inline 00747 void boundaryTensor3(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00748 pair<DestIteratorEven, DestAccessorEven> even, 00749 pair<DestIteratorOdd, DestAccessorOdd> odd, 00750 double scale) 00751 { 00752 boundaryTensor3(src.first, src.second, src.third, 00753 even.first, even.second, odd.first, odd.second, scale); 00754 } 00755 00756 //@} 00757 00758 } // namespace vigra 00759 00760 #endif // VIGRA_BOUNDARYTENSOR_HXX
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|