[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]
![]() |
vigra/separableconvolution.hxx | ![]() |
---|
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2002 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_SEPARABLECONVOLUTION_HXX 00040 #define VIGRA_SEPARABLECONVOLUTION_HXX 00041 00042 #include <cmath> 00043 #include "vigra/utilities.hxx" 00044 #include "vigra/numerictraits.hxx" 00045 #include "vigra/imageiteratoradapter.hxx" 00046 #include "vigra/bordertreatment.hxx" 00047 #include "vigra/gaussians.hxx" 00048 #include "vigra/array_vector.hxx" 00049 00050 namespace vigra { 00051 00052 /********************************************************/ 00053 /* */ 00054 /* internalConvolveLineWrap */ 00055 /* */ 00056 /********************************************************/ 00057 00058 template <class SrcIterator, class SrcAccessor, 00059 class DestIterator, class DestAccessor, 00060 class KernelIterator, class KernelAccessor> 00061 void internalConvolveLineWrap(SrcIterator is, SrcIterator iend, SrcAccessor sa, 00062 DestIterator id, DestAccessor da, 00063 KernelIterator kernel, KernelAccessor ka, 00064 int kleft, int kright) 00065 { 00066 // int w = iend - is; 00067 int w = std::distance( is, iend ); 00068 00069 typedef typename NumericTraits<typename 00070 SrcAccessor::value_type>::RealPromote SumType; 00071 00072 SrcIterator ibegin = is; 00073 00074 for(int x=0; x<w; ++x, ++is, ++id) 00075 { 00076 KernelIterator ik = kernel + kright; 00077 SumType sum = NumericTraits<SumType>::zero(); 00078 00079 if(x < kright) 00080 { 00081 int x0 = x - kright; 00082 SrcIterator iss = iend + x0; 00083 00084 for(; x0; ++x0, --ik, ++iss) 00085 { 00086 sum += ka(ik) * sa(iss); 00087 } 00088 00089 iss = ibegin; 00090 SrcIterator isend = is + (1 - kleft); 00091 for(; iss != isend ; --ik, ++iss) 00092 { 00093 sum += ka(ik) * sa(iss); 00094 } 00095 } 00096 else if(w-x <= -kleft) 00097 { 00098 SrcIterator iss = is + (-kright); 00099 SrcIterator isend = iend; 00100 for(; iss != isend ; --ik, ++iss) 00101 { 00102 sum += ka(ik) * sa(iss); 00103 } 00104 00105 int x0 = -kleft - w + x + 1; 00106 iss = ibegin; 00107 00108 for(; x0; --x0, --ik, ++iss) 00109 { 00110 sum += ka(ik) * sa(iss); 00111 } 00112 } 00113 else 00114 { 00115 SrcIterator iss = is - kright; 00116 SrcIterator isend = is + (1 - kleft); 00117 for(; iss != isend ; --ik, ++iss) 00118 { 00119 sum += ka(ik) * sa(iss); 00120 } 00121 } 00122 00123 da.set(NumericTraits<typename 00124 DestAccessor::value_type>::fromRealPromote(sum), id); 00125 } 00126 } 00127 00128 /********************************************************/ 00129 /* */ 00130 /* internalConvolveLineClip */ 00131 /* */ 00132 /********************************************************/ 00133 00134 template <class SrcIterator, class SrcAccessor, 00135 class DestIterator, class DestAccessor, 00136 class KernelIterator, class KernelAccessor, 00137 class Norm> 00138 void internalConvolveLineClip(SrcIterator is, SrcIterator iend, SrcAccessor sa, 00139 DestIterator id, DestAccessor da, 00140 KernelIterator kernel, KernelAccessor ka, 00141 int kleft, int kright, Norm norm) 00142 { 00143 // int w = iend - is; 00144 int w = std::distance( is, iend ); 00145 00146 typedef typename NumericTraits<typename 00147 SrcAccessor::value_type>::RealPromote SumType; 00148 00149 SrcIterator ibegin = is; 00150 00151 for(int x=0; x<w; ++x, ++is, ++id) 00152 { 00153 KernelIterator ik = kernel + kright; 00154 SumType sum = NumericTraits<SumType>::zero(); 00155 00156 if(x < kright) 00157 { 00158 int x0 = x - kright; 00159 Norm clipped = NumericTraits<Norm>::zero(); 00160 00161 for(; x0; ++x0, --ik) 00162 { 00163 clipped += ka(ik); 00164 } 00165 00166 SrcIterator iss = ibegin; 00167 SrcIterator isend = is + (1 - kleft); 00168 for(; iss != isend ; --ik, ++iss) 00169 { 00170 sum += ka(ik) * sa(iss); 00171 } 00172 00173 sum = norm / (norm - clipped) * sum; 00174 } 00175 else if(w-x <= -kleft) 00176 { 00177 SrcIterator iss = is + (-kright); 00178 SrcIterator isend = iend; 00179 for(; iss != isend ; --ik, ++iss) 00180 { 00181 sum += ka(ik) * sa(iss); 00182 } 00183 00184 Norm clipped = NumericTraits<Norm>::zero(); 00185 00186 int x0 = -kleft - w + x + 1; 00187 00188 for(; x0; --x0, --ik) 00189 { 00190 clipped += ka(ik); 00191 } 00192 00193 sum = norm / (norm - clipped) * sum; 00194 } 00195 else 00196 { 00197 SrcIterator iss = is + (-kright); 00198 SrcIterator isend = is + (1 - kleft); 00199 for(; iss != isend ; --ik, ++iss) 00200 { 00201 sum += ka(ik) * sa(iss); 00202 } 00203 } 00204 00205 da.set(NumericTraits<typename 00206 DestAccessor::value_type>::fromRealPromote(sum), id); 00207 } 00208 } 00209 00210 /********************************************************/ 00211 /* */ 00212 /* internalConvolveLineReflect */ 00213 /* */ 00214 /********************************************************/ 00215 00216 template <class SrcIterator, class SrcAccessor, 00217 class DestIterator, class DestAccessor, 00218 class KernelIterator, class KernelAccessor> 00219 void internalConvolveLineReflect(SrcIterator is, SrcIterator iend, SrcAccessor sa, 00220 DestIterator id, DestAccessor da, 00221 KernelIterator kernel, KernelAccessor ka, 00222 int kleft, int kright) 00223 { 00224 // int w = iend - is; 00225 int w = std::distance( is, iend ); 00226 00227 typedef typename NumericTraits<typename 00228 SrcAccessor::value_type>::RealPromote SumType; 00229 00230 SrcIterator ibegin = is; 00231 00232 for(int x=0; x<w; ++x, ++is, ++id) 00233 { 00234 KernelIterator ik = kernel + kright; 00235 SumType sum = NumericTraits<SumType>::zero(); 00236 00237 if(x < kright) 00238 { 00239 int x0 = x - kright; 00240 SrcIterator iss = ibegin - x0; 00241 00242 for(; x0; ++x0, --ik, --iss) 00243 { 00244 sum += ka(ik) * sa(iss); 00245 } 00246 00247 SrcIterator isend = is + (1 - kleft); 00248 for(; iss != isend ; --ik, ++iss) 00249 { 00250 sum += ka(ik) * sa(iss); 00251 } 00252 } 00253 else if(w-x <= -kleft) 00254 { 00255 SrcIterator iss = is + (-kright); 00256 SrcIterator isend = iend; 00257 for(; iss != isend ; --ik, ++iss) 00258 { 00259 sum += ka(ik) * sa(iss); 00260 } 00261 00262 int x0 = -kleft - w + x + 1; 00263 iss = iend - 2; 00264 00265 for(; x0; --x0, --ik, --iss) 00266 { 00267 sum += ka(ik) * sa(iss); 00268 } 00269 } 00270 else 00271 { 00272 SrcIterator iss = is + (-kright); 00273 SrcIterator isend = is + (1 - kleft); 00274 for(; iss != isend ; --ik, ++iss) 00275 { 00276 sum += ka(ik) * sa(iss); 00277 } 00278 } 00279 00280 da.set(NumericTraits<typename 00281 DestAccessor::value_type>::fromRealPromote(sum), id); 00282 } 00283 } 00284 00285 /********************************************************/ 00286 /* */ 00287 /* internalConvolveLineRepeat */ 00288 /* */ 00289 /********************************************************/ 00290 00291 template <class SrcIterator, class SrcAccessor, 00292 class DestIterator, class DestAccessor, 00293 class KernelIterator, class KernelAccessor> 00294 void internalConvolveLineRepeat(SrcIterator is, SrcIterator iend, SrcAccessor sa, 00295 DestIterator id, DestAccessor da, 00296 KernelIterator kernel, KernelAccessor ka, 00297 int kleft, int kright) 00298 { 00299 // int w = iend - is; 00300 int w = std::distance( is, iend ); 00301 00302 typedef typename NumericTraits<typename 00303 SrcAccessor::value_type>::RealPromote SumType; 00304 00305 SrcIterator ibegin = is; 00306 00307 for(int x=0; x<w; ++x, ++is, ++id) 00308 { 00309 KernelIterator ik = kernel + kright; 00310 SumType sum = NumericTraits<SumType>::zero(); 00311 00312 if(x < kright) 00313 { 00314 int x0 = x - kright; 00315 SrcIterator iss = ibegin; 00316 00317 for(; x0; ++x0, --ik) 00318 { 00319 sum += ka(ik) * sa(iss); 00320 } 00321 00322 SrcIterator isend = is + (1 - kleft); 00323 for(; iss != isend ; --ik, ++iss) 00324 { 00325 sum += ka(ik) * sa(iss); 00326 } 00327 } 00328 else if(w-x <= -kleft) 00329 { 00330 SrcIterator iss = is + (-kright); 00331 SrcIterator isend = iend; 00332 for(; iss != isend ; --ik, ++iss) 00333 { 00334 sum += ka(ik) * sa(iss); 00335 } 00336 00337 int x0 = -kleft - w + x + 1; 00338 iss = iend - 1; 00339 00340 for(; x0; --x0, --ik) 00341 { 00342 sum += ka(ik) * sa(iss); 00343 } 00344 } 00345 else 00346 { 00347 SrcIterator iss = is + (-kright); 00348 SrcIterator isend = is + (1 - kleft); 00349 for(; iss != isend ; --ik, ++iss) 00350 { 00351 sum += ka(ik) * sa(iss); 00352 } 00353 } 00354 00355 da.set(NumericTraits<typename 00356 DestAccessor::value_type>::fromRealPromote(sum), id); 00357 } 00358 } 00359 00360 /********************************************************/ 00361 /* */ 00362 /* internalConvolveLineAvoid */ 00363 /* */ 00364 /********************************************************/ 00365 00366 template <class SrcIterator, class SrcAccessor, 00367 class DestIterator, class DestAccessor, 00368 class KernelIterator, class KernelAccessor> 00369 void internalConvolveLineAvoid(SrcIterator is, SrcIterator iend, SrcAccessor sa, 00370 DestIterator id, DestAccessor da, 00371 KernelIterator kernel, KernelAccessor ka, 00372 int kleft, int kright) 00373 { 00374 // int w = iend - is; 00375 int w = std::distance( is, iend ); 00376 00377 typedef typename NumericTraits<typename 00378 SrcAccessor::value_type>::RealPromote SumType; 00379 00380 is += kright; 00381 id += kright; 00382 00383 for(int x=kright; x<w+kleft; ++x, ++is, ++id) 00384 { 00385 KernelIterator ik = kernel + kright; 00386 SumType sum = NumericTraits<SumType>::zero(); 00387 00388 SrcIterator iss = is + (-kright); 00389 SrcIterator isend = is + (1 - kleft); 00390 for(; iss != isend ; --ik, ++iss) 00391 { 00392 sum += ka(ik) * sa(iss); 00393 } 00394 00395 da.set(NumericTraits<typename 00396 DestAccessor::value_type>::fromRealPromote(sum), id); 00397 } 00398 } 00399 00400 /********************************************************/ 00401 /* */ 00402 /* Separable convolution functions */ 00403 /* */ 00404 /********************************************************/ 00405 00406 /** \addtogroup SeparableConvolution One-dimensional and separable convolution functions 00407 00408 Perform 1D convolution and separable filtering in 2 dimensions. 00409 00410 These generic convolution functions implement 00411 the standard convolution operation for a wide range of images and 00412 signals that fit into the required interface. They need a suitable 00413 kernel to operate. 00414 */ 00415 //@{ 00416 00417 /** \brief Performs a 1 dimensional convolution of the source signal using the given 00418 kernel. 00419 00420 The KernelIterator must point to the center iterator, and 00421 the kernel's size is given by its left (kleft <= 0) and right 00422 (kright >= 0) borders. The signal must always be larger than the kernel. 00423 At those positions where the kernel does not completely fit 00424 into the signal's range, the specified \ref BorderTreatmentMode is 00425 applied. 00426 00427 The signal's value_type (SrcAccessor::value_type) must be a 00428 linear space over the kernel's value_type (KernelAccessor::value_type), 00429 i.e. addition of source values, multiplication with kernel values, 00430 and NumericTraits must be defined. 00431 The kernel's value_type must be an algebraic field, 00432 i.e. the arithmetic operations (+, -, *, /) and NumericTraits must 00433 be defined. 00434 00435 <b> Declarations:</b> 00436 00437 pass arguments explicitly: 00438 \code 00439 namespace vigra { 00440 template <class SrcIterator, class SrcAccessor, 00441 class DestIterator, class DestAccessor, 00442 class KernelIterator, class KernelAccessor> 00443 void convolveLine(SrcIterator is, SrcIterator isend, SrcAccessor sa, 00444 DestIterator id, DestAccessor da, 00445 KernelIterator ik, KernelAccessor ka, 00446 int kleft, int kright, BorderTreatmentMode border) 00447 } 00448 \endcode 00449 00450 00451 use argument objects in conjunction with \ref ArgumentObjectFactories: 00452 \code 00453 namespace vigra { 00454 template <class SrcIterator, class SrcAccessor, 00455 class DestIterator, class DestAccessor, 00456 class KernelIterator, class KernelAccessor> 00457 void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00458 pair<DestIterator, DestAccessor> dest, 00459 tuple5<KernelIterator, KernelAccessor, 00460 int, int, BorderTreatmentMode> kernel) 00461 } 00462 \endcode 00463 00464 <b> Usage:</b> 00465 00466 <b>\#include</b> "<a href="separableconvolution_8hxx-source.html">vigra/separableconvolution.hxx</a>" 00467 00468 00469 \code 00470 std::vector<float> src, dest; 00471 ... 00472 00473 // define binomial filter of size 5 00474 static float kernel[] = 00475 { 1.0/16.0, 4.0/16.0, 6.0/16.0, 4.0/16.0, 1.0/16.0}; 00476 00477 typedef vigra::StandardAccessor<float> FAccessor; 00478 typedef vigra::StandardAccessor<float> KernelAccessor; 00479 00480 00481 vigra::convolveLine(src.begin(), src.end(), FAccessor(), dest.begin(), FAccessor(), 00482 kernel+2, KernelAccessor(), -2, 2, BORDER_TREATMENT_REFLECT); 00483 // ^^^^^^^^ this is the center of the kernel 00484 00485 \endcode 00486 00487 <b> Required Interface:</b> 00488 00489 \code 00490 RandomAccessIterator is, isend; 00491 RandomAccessIterator id; 00492 RandomAccessIterator ik; 00493 00494 SrcAccessor src_accessor; 00495 DestAccessor dest_accessor; 00496 KernelAccessor kernel_accessor; 00497 00498 NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(is); 00499 00500 s = s + s; 00501 s = kernel_accessor(ik) * s; 00502 00503 dest_accessor.set( 00504 NumericTraits<DestAccessor::value_type>::fromRealPromote(s), id); 00505 00506 \endcode 00507 00508 If border == BORDER_TREATMENT_CLIP: 00509 00510 \code 00511 NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik); 00512 00513 k = k + k; 00514 k = k - k; 00515 k = k * k; 00516 k = k / k; 00517 00518 \endcode 00519 00520 <b> Preconditions:</b> 00521 00522 \code 00523 kleft <= 0 00524 kright >= 0 00525 iend - is >= kright + kleft + 1 00526 \endcode 00527 00528 If border == BORDER_TREATMENT_CLIP: Sum of kernel elements must be 00529 != 0. 00530 00531 */ 00532 template <class SrcIterator, class SrcAccessor, 00533 class DestIterator, class DestAccessor, 00534 class KernelIterator, class KernelAccessor> 00535 void convolveLine(SrcIterator is, SrcIterator iend, SrcAccessor sa, 00536 DestIterator id, DestAccessor da, 00537 KernelIterator ik, KernelAccessor ka, 00538 int kleft, int kright, BorderTreatmentMode border) 00539 { 00540 typedef typename KernelAccessor::value_type KernelValue; 00541 00542 vigra_precondition(kleft <= 0, 00543 "convolveLine(): kleft must be <= 0.\n"); 00544 vigra_precondition(kright >= 0, 00545 "convolveLine(): kright must be >= 0.\n"); 00546 00547 // int w = iend - is; 00548 int w = std::distance( is, iend ); 00549 00550 vigra_precondition(w >= kright - kleft + 1, 00551 "convolveLine(): kernel longer than line\n"); 00552 00553 switch(border) 00554 { 00555 case BORDER_TREATMENT_WRAP: 00556 { 00557 internalConvolveLineWrap(is, iend, sa, id, da, ik, ka, kleft, kright); 00558 break; 00559 } 00560 case BORDER_TREATMENT_AVOID: 00561 { 00562 internalConvolveLineAvoid(is, iend, sa, id, da, ik, ka, kleft, kright); 00563 break; 00564 } 00565 case BORDER_TREATMENT_REFLECT: 00566 { 00567 internalConvolveLineReflect(is, iend, sa, id, da, ik, ka, kleft, kright); 00568 break; 00569 } 00570 case BORDER_TREATMENT_REPEAT: 00571 { 00572 internalConvolveLineRepeat(is, iend, sa, id, da, ik, ka, kleft, kright); 00573 break; 00574 } 00575 case BORDER_TREATMENT_CLIP: 00576 { 00577 // find norm of kernel 00578 typedef typename KernelAccessor::value_type KT; 00579 KT norm = NumericTraits<KT>::zero(); 00580 KernelIterator iik = ik + kleft; 00581 for(int i=kleft; i<=kright; ++i, ++iik) norm += ka(iik); 00582 00583 vigra_precondition(norm != NumericTraits<KT>::zero(), 00584 "convolveLine(): Norm of kernel must be != 0" 00585 " in mode BORDER_TREATMENT_CLIP.\n"); 00586 00587 internalConvolveLineClip(is, iend, sa, id, da, ik, ka, kleft, kright, norm); 00588 break; 00589 } 00590 default: 00591 { 00592 vigra_precondition(0, 00593 "convolveLine(): Unknown border treatment mode.\n"); 00594 } 00595 } 00596 } 00597 00598 template <class SrcIterator, class SrcAccessor, 00599 class DestIterator, class DestAccessor, 00600 class KernelIterator, class KernelAccessor> 00601 inline 00602 void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00603 pair<DestIterator, DestAccessor> dest, 00604 tuple5<KernelIterator, KernelAccessor, 00605 int, int, BorderTreatmentMode> kernel) 00606 { 00607 convolveLine(src.first, src.second, src.third, 00608 dest.first, dest.second, 00609 kernel.first, kernel.second, 00610 kernel.third, kernel.fourth, kernel.fifth); 00611 } 00612 00613 /********************************************************/ 00614 /* */ 00615 /* separableConvolveX */ 00616 /* */ 00617 /********************************************************/ 00618 00619 /** \brief Performs a 1 dimensional convolution in x direction. 00620 00621 It calls \link SeparableConvolution#convolveLine convolveLine\endlink() for every row of the 00622 image. See \link SeparableConvolution#convolveLine convolveLine\endlink() for more information about required interfaces 00623 and vigra_preconditions. 00624 00625 <b> Declarations:</b> 00626 00627 pass arguments explicitly: 00628 \code 00629 namespace vigra { 00630 template <class SrcImageIterator, class SrcAccessor, 00631 class DestImageIterator, class DestAccessor, 00632 class KernelIterator, class KernelAccessor> 00633 void separableConvolveX(SrcImageIterator supperleft, 00634 SrcImageIterator slowerright, SrcAccessor sa, 00635 DestImageIterator dupperleft, DestAccessor da, 00636 KernelIterator ik, KernelAccessor ka, 00637 int kleft, int kright, BorderTreatmentMode border) 00638 } 00639 \endcode 00640 00641 00642 use argument objects in conjunction with \ref ArgumentObjectFactories: 00643 \code 00644 namespace vigra { 00645 template <class SrcImageIterator, class SrcAccessor, 00646 class DestImageIterator, class DestAccessor, 00647 class KernelIterator, class KernelAccessor> 00648 void separableConvolveX(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 00649 pair<DestImageIterator, DestAccessor> dest, 00650 tuple5<KernelIterator, KernelAccessor, 00651 int, int, BorderTreatmentMode> kernel) 00652 } 00653 \endcode 00654 00655 <b> Usage:</b> 00656 00657 <b>\#include</b> "<a href="separableconvolution_8hxx-source.html">vigra/separableconvolution.hxx</a>" 00658 00659 00660 \code 00661 vigra::FImage src(w,h), dest(w,h); 00662 ... 00663 00664 // define Gaussian kernel with std. deviation 3.0 00665 vigra::Kernel1D<double> kernel; 00666 kernel.initGaussian(3.0); 00667 00668 vigra::separableConvolveX(srcImageRange(src), destImage(dest), kernel1d(kernel)); 00669 00670 \endcode 00671 00672 */ 00673 template <class SrcIterator, class SrcAccessor, 00674 class DestIterator, class DestAccessor, 00675 class KernelIterator, class KernelAccessor> 00676 void separableConvolveX(SrcIterator supperleft, 00677 SrcIterator slowerright, SrcAccessor sa, 00678 DestIterator dupperleft, DestAccessor da, 00679 KernelIterator ik, KernelAccessor ka, 00680 int kleft, int kright, BorderTreatmentMode border) 00681 { 00682 typedef typename KernelAccessor::value_type KernelValue; 00683 00684 vigra_precondition(kleft <= 0, 00685 "separableConvolveX(): kleft must be <= 0.\n"); 00686 vigra_precondition(kright >= 0, 00687 "separableConvolveX(): kright must be >= 0.\n"); 00688 00689 int w = slowerright.x - supperleft.x; 00690 int h = slowerright.y - supperleft.y; 00691 00692 vigra_precondition(w >= kright - kleft + 1, 00693 "separableConvolveX(): kernel longer than line\n"); 00694 00695 int y; 00696 00697 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y) 00698 { 00699 typename SrcIterator::row_iterator rs = supperleft.rowIterator(); 00700 typename DestIterator::row_iterator rd = dupperleft.rowIterator(); 00701 00702 convolveLine(rs, rs+w, sa, rd, da, 00703 ik, ka, kleft, kright, border); 00704 } 00705 } 00706 00707 template <class SrcIterator, class SrcAccessor, 00708 class DestIterator, class DestAccessor, 00709 class KernelIterator, class KernelAccessor> 00710 inline void 00711 separableConvolveX(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00712 pair<DestIterator, DestAccessor> dest, 00713 tuple5<KernelIterator, KernelAccessor, 00714 int, int, BorderTreatmentMode> kernel) 00715 { 00716 separableConvolveX(src.first, src.second, src.third, 00717 dest.first, dest.second, 00718 kernel.first, kernel.second, 00719 kernel.third, kernel.fourth, kernel.fifth); 00720 } 00721 00722 00723 00724 /********************************************************/ 00725 /* */ 00726 /* separableConvolveY */ 00727 /* */ 00728 /********************************************************/ 00729 00730 /** \brief Performs a 1 dimensional convolution in y direction. 00731 00732 It calls \link SeparableConvolution#convolveLine convolveLine\endlink() for every column of the 00733 image. See \link SeparableConvolution#convolveLine convolveLine\endlink() for more information about required interfaces 00734 and vigra_preconditions. 00735 00736 <b> Declarations:</b> 00737 00738 pass arguments explicitly: 00739 \code 00740 namespace vigra { 00741 template <class SrcImageIterator, class SrcAccessor, 00742 class DestImageIterator, class DestAccessor, 00743 class KernelIterator, class KernelAccessor> 00744 void separableConvolveY(SrcImageIterator supperleft, 00745 SrcImageIterator slowerright, SrcAccessor sa, 00746 DestImageIterator dupperleft, DestAccessor da, 00747 KernelIterator ik, KernelAccessor ka, 00748 int kleft, int kright, BorderTreatmentMode border) 00749 } 00750 \endcode 00751 00752 00753 use argument objects in conjunction with \ref ArgumentObjectFactories: 00754 \code 00755 namespace vigra { 00756 template <class SrcImageIterator, class SrcAccessor, 00757 class DestImageIterator, class DestAccessor, 00758 class KernelIterator, class KernelAccessor> 00759 void separableConvolveY(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 00760 pair<DestImageIterator, DestAccessor> dest, 00761 tuple5<KernelIterator, KernelAccessor, 00762 int, int, BorderTreatmentMode> kernel) 00763 } 00764 \endcode 00765 00766 <b> Usage:</b> 00767 00768 <b>\#include</b> "<a href="separableconvolution_8hxx-source.html">vigra/separableconvolution.hxx</a>" 00769 00770 00771 \code 00772 vigra::FImage src(w,h), dest(w,h); 00773 ... 00774 00775 // define Gaussian kernel with std. deviation 3.0 00776 vigra::Kernel1D kernel; 00777 kernel.initGaussian(3.0); 00778 00779 vigra::separableConvolveY(srcImageRange(src), destImage(dest), kernel1d(kernel)); 00780 00781 \endcode 00782 00783 */ 00784 template <class SrcIterator, class SrcAccessor, 00785 class DestIterator, class DestAccessor, 00786 class KernelIterator, class KernelAccessor> 00787 void separableConvolveY(SrcIterator supperleft, 00788 SrcIterator slowerright, SrcAccessor sa, 00789 DestIterator dupperleft, DestAccessor da, 00790 KernelIterator ik, KernelAccessor ka, 00791 int kleft, int kright, BorderTreatmentMode border) 00792 { 00793 typedef typename KernelAccessor::value_type KernelValue; 00794 00795 vigra_precondition(kleft <= 0, 00796 "separableConvolveY(): kleft must be <= 0.\n"); 00797 vigra_precondition(kright >= 0, 00798 "separableConvolveY(): kright must be >= 0.\n"); 00799 00800 int w = slowerright.x - supperleft.x; 00801 int h = slowerright.y - supperleft.y; 00802 00803 vigra_precondition(h >= kright - kleft + 1, 00804 "separableConvolveY(): kernel longer than line\n"); 00805 00806 int x; 00807 00808 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x) 00809 { 00810 typename SrcIterator::column_iterator cs = supperleft.columnIterator(); 00811 typename DestIterator::column_iterator cd = dupperleft.columnIterator(); 00812 00813 convolveLine(cs, cs+h, sa, cd, da, 00814 ik, ka, kleft, kright, border); 00815 } 00816 } 00817 00818 template <class SrcIterator, class SrcAccessor, 00819 class DestIterator, class DestAccessor, 00820 class KernelIterator, class KernelAccessor> 00821 inline void 00822 separableConvolveY(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00823 pair<DestIterator, DestAccessor> dest, 00824 tuple5<KernelIterator, KernelAccessor, 00825 int, int, BorderTreatmentMode> kernel) 00826 { 00827 separableConvolveY(src.first, src.second, src.third, 00828 dest.first, dest.second, 00829 kernel.first, kernel.second, 00830 kernel.third, kernel.fourth, kernel.fifth); 00831 } 00832 00833 //@} 00834 00835 /********************************************************/ 00836 /* */ 00837 /* Kernel1D */ 00838 /* */ 00839 /********************************************************/ 00840 00841 /** \brief Generic 1 dimensional convolution kernel. 00842 00843 This kernel may be used for convolution of 1 dimensional signals or for 00844 separable convolution of multidimensional signals. 00845 00846 Convlution functions access the kernel via a 1 dimensional random access 00847 iterator which they get by calling \ref center(). This iterator 00848 points to the center of the kernel. The kernel's size is given by its left() (<=0) 00849 and right() (>= 0) methods. The desired border treatment mode is 00850 returned by borderTreatment(). 00851 00852 The different init functions create a kernel with the specified 00853 properties. The kernel's value_type must be a linear space, i.e. it 00854 must define multiplication with doubles and NumericTraits. 00855 00856 00857 The kernel defines a factory function kernel1d() to create an argument object 00858 (see \ref KernelArgumentObjectFactories). 00859 00860 <b> Usage:</b> 00861 00862 <b>\#include</b> "<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>" 00863 00864 \code 00865 vigra::FImage src(w,h), dest(w,h); 00866 ... 00867 00868 // define Gaussian kernel with std. deviation 3.0 00869 vigra::Kernel1D kernel; 00870 kernel.initGaussian(3.0); 00871 00872 vigra::separableConvolveX(srcImageRange(src), destImage(dest), kernel1d(kernel)); 00873 \endcode 00874 00875 <b> Required Interface:</b> 00876 00877 \code 00878 value_type v = vigra::NumericTraits<value_type>::one(); // if norm is not 00879 // given explicitly 00880 double d; 00881 00882 v = d * v; 00883 \endcode 00884 */ 00885 00886 template <class ARITHTYPE> 00887 class Kernel1D 00888 { 00889 public: 00890 typedef ArrayVector<ARITHTYPE> InternalVector; 00891 00892 /** the kernel's value type 00893 */ 00894 typedef typename InternalVector::value_type value_type; 00895 00896 /** the kernel's reference type 00897 */ 00898 typedef typename InternalVector::reference reference; 00899 00900 /** the kernel's const reference type 00901 */ 00902 typedef typename InternalVector::const_reference const_reference; 00903 00904 /** deprecated -- use Kernel1D::iterator 00905 */ 00906 typedef typename InternalVector::iterator Iterator; 00907 00908 /** 1D random access iterator over the kernel's values 00909 */ 00910 typedef typename InternalVector::iterator iterator; 00911 00912 /** const 1D random access iterator over the kernel's values 00913 */ 00914 typedef typename InternalVector::const_iterator const_iterator; 00915 00916 /** the kernel's accessor 00917 */ 00918 typedef StandardAccessor<ARITHTYPE> Accessor; 00919 00920 /** the kernel's const accessor 00921 */ 00922 typedef StandardConstAccessor<ARITHTYPE> ConstAccessor; 00923 00924 struct InitProxy 00925 { 00926 InitProxy(Iterator i, int count, value_type & norm) 00927 : iter_(i), base_(i), 00928 count_(count), sum_(count), 00929 norm_(norm) 00930 {} 00931 00932 ~InitProxy() 00933 { 00934 vigra_precondition(count_ == 1 || count_ == sum_, 00935 "Kernel1D::initExplicitly(): " 00936 "Too few init values."); 00937 } 00938 00939 InitProxy & operator,(value_type const & v) 00940 { 00941 if(sum_ == count_) norm_ = *iter_; 00942 00943 norm_ += v; 00944 00945 --count_; 00946 vigra_precondition(count_ > 0, 00947 "Kernel1D::initExplicitly(): " 00948 "Too many init values."); 00949 00950 ++iter_; 00951 *iter_ = v; 00952 00953 return *this; 00954 } 00955 00956 Iterator iter_, base_; 00957 int count_, sum_; 00958 value_type & norm_; 00959 }; 00960 00961 static value_type one() { return NumericTraits<value_type>::one(); } 00962 00963 /** Default constructor. 00964 Creates a kernel of size 1 which would copy the signal 00965 unchanged. 00966 */ 00967 Kernel1D() 00968 : kernel_(), 00969 left_(0), 00970 right_(0), 00971 border_treatment_(BORDER_TREATMENT_CLIP), 00972 norm_(one()) 00973 { 00974 kernel_.push_back(norm_); 00975 } 00976 00977 /** Copy constructor. 00978 */ 00979 Kernel1D(Kernel1D const & k) 00980 : kernel_(k.kernel_), 00981 left_(k.left_), 00982 right_(k.right_), 00983 border_treatment_(k.border_treatment_), 00984 norm_(k.norm_) 00985 {} 00986 00987 /** Copy assignment. 00988 */ 00989 Kernel1D & operator=(Kernel1D const & k) 00990 { 00991 if(this != &k) 00992 { 00993 left_ = k.left_; 00994 right_ = k.right_; 00995 border_treatment_ = k.border_treatment_; 00996 norm_ = k.norm_; 00997 kernel_ = k.kernel_; 00998 } 00999 return *this; 01000 } 01001 01002 /** Initialization. 01003 This initializes the kernel with the given constant. The norm becomes 01004 v*size(). 01005 01006 Instead of a single value an initializer list of length size() 01007 can be used like this: 01008 01009 \code 01010 vigra::Kernel1D<float> roberts_gradient_x; 01011 01012 roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0; 01013 \endcode 01014 01015 In this case, the norm will be set to the sum of the init values. 01016 An initializer list of wrong length will result in a run-time error. 01017 */ 01018 InitProxy operator=(value_type const & v) 01019 { 01020 int size = right_ - left_ + 1; 01021 for(unsigned int i=0; i<kernel_.size(); ++i) kernel_[i] = v; 01022 norm_ = (double)size*v; 01023 01024 return InitProxy(kernel_.begin(), size, norm_); 01025 } 01026 01027 /** Destructor. 01028 */ 01029 ~Kernel1D() 01030 {} 01031 01032 /** 01033 Init as a sampled Gaussian function. The radius of the kernel is 01034 always 3*std_dev. '<tt>norm</tt>' denotes the sum of all bins of the kernel 01035 (i.e. the kernel is corrected for the normalization error introduced 01036 by windowing the Gaussian to a finite interval). However, 01037 if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic 01038 expression for the Gaussian, and <b>no</b> correction for the windowing 01039 error is performed. 01040 01041 Precondition: 01042 \code 01043 std_dev >= 0.0 01044 \endcode 01045 01046 Postconditions: 01047 \code 01048 1. left() == -(int)(3.0*std_dev + 0.5) 01049 2. right() == (int)(3.0*std_dev + 0.5) 01050 3. borderTreatment() == BORDER_TREATMENT_CLIP 01051 4. norm() == norm 01052 \endcode 01053 */ 01054 void initGaussian(double std_dev, value_type norm); 01055 01056 /** Init as a Gaussian function with norm 1. 01057 */ 01058 void initGaussian(double std_dev) 01059 { 01060 initGaussian(std_dev, one()); 01061 } 01062 01063 01064 /** 01065 Init as Lindeberg's discrete analog of the Gaussian function. The radius of the kernel is 01066 always 3*std_dev. 'norm' denotes the sum of all bins of the kernel. 01067 01068 Precondition: 01069 \code 01070 std_dev >= 0.0 01071 \endcode 01072 01073 Postconditions: 01074 \code 01075 1. left() == -(int)(3.0*std_dev + 0.5) 01076 2. right() == (int)(3.0*std_dev + 0.5) 01077 3. borderTreatment() == BORDER_TREATMENT_REFLECT 01078 4. norm() == norm 01079 \endcode 01080 */ 01081 void initDiscreteGaussian(double std_dev, value_type norm); 01082 01083 /** Init as a LOineberg's discrete analog of the Gaussian function 01084 with norm 1. 01085 */ 01086 void initDiscreteGaussian(double std_dev) 01087 { 01088 initDiscreteGaussian(std_dev, one()); 01089 } 01090 01091 /** 01092 Init as a Gaussian derivative of order '<tt>order</tt>'. 01093 The radius of the kernel is always <tt>3*std_dev + 0.5*order</tt>. 01094 '<tt>norm</tt>' denotes the norm of the kernel so that the 01095 following condition is fulfilled: 01096 01097 \f[ \sum_{i=left()}^{right()} 01098 \frac{(-i)^{order}kernel[i]}{order!} = norm 01099 \f] 01100 01101 Thus, the kernel will be corrected for the error introduced 01102 by windowing the Gaussian to a finite interval. However, 01103 if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic 01104 expression for the Gaussian derivative, and <b>no</b> correction for the 01105 windowing error is performed. 01106 01107 Preconditions: 01108 \code 01109 1. std_dev >= 0.0 01110 2. order >= 1 01111 \endcode 01112 01113 Postconditions: 01114 \code 01115 1. left() == -(int)(3.0*std_dev + 0.5*order + 0.5) 01116 2. right() == (int)(3.0*std_dev + 0.5*order + 0.5) 01117 3. borderTreatment() == BORDER_TREATMENT_REPEAT 01118 4. norm() == norm 01119 \endcode 01120 */ 01121 void initGaussianDerivative(double std_dev, int order, value_type norm); 01122 01123 /** Init as a Gaussian derivative with norm 1. 01124 */ 01125 void initGaussianDerivative(double std_dev, int order) 01126 { 01127 initGaussianDerivative(std_dev, order, one()); 01128 } 01129 01130 /** 01131 Init as a Binomial filter. 'norm' denotes the sum of all bins 01132 of the kernel. 01133 01134 Precondition: 01135 \code 01136 radius >= 0 01137 \endcode 01138 01139 Postconditions: 01140 \code 01141 1. left() == -radius 01142 2. right() == radius 01143 3. borderTreatment() == BORDER_TREATMENT_REFLECT 01144 4. norm() == norm 01145 \endcode 01146 */ 01147 void initBinomial(int radius, value_type norm); 01148 01149 /** Init as a Binomial filter with norm 1. 01150 */ 01151 void initBinomial(int radius) 01152 { 01153 initBinomial(radius, one()); 01154 } 01155 01156 /** 01157 Init as an Averaging filter. 'norm' denotes the sum of all bins 01158 of the kernel. The window size is (2*radius+1) * (2*radius+1) 01159 01160 Precondition: 01161 \code 01162 radius >= 0 01163 \endcode 01164 01165 Postconditions: 01166 \code 01167 1. left() == -radius 01168 2. right() == radius 01169 3. borderTreatment() == BORDER_TREATMENT_CLIP 01170 4. norm() == norm 01171 \endcode 01172 */ 01173 void initAveraging(int radius, value_type norm); 01174 01175 /** Init as a Averaging filter with norm 1. 01176 */ 01177 void initAveraging(int radius) 01178 { 01179 initAveraging(radius, one()); 01180 } 01181 01182 /** 01183 Init as a symmetric gradient filter of the form 01184 <TT>[ 0.5 * norm, 0.0 * norm, -0.5 * norm]</TT> 01185 01186 Postconditions: 01187 \code 01188 1. left() == -1 01189 2. right() == 1 01190 3. borderTreatment() == BORDER_TREATMENT_REPEAT 01191 4. norm() == norm 01192 \endcode 01193 */ 01194 void 01195 initSymmetricGradient(value_type norm ); 01196 01197 /** Init as a symmetric gradient filter with norm 1. 01198 */ 01199 void initSymmetricGradient() 01200 { 01201 initSymmetricGradient(one()); 01202 } 01203 01204 /** Init the kernel by an explicit initializer list. 01205 The left and right boundaries of the kernel must be passed. 01206 A comma-separated initializer list is given after the assignment 01207 operator. This function is used like this: 01208 01209 \code 01210 // define horizontal Roberts filter 01211 vigra::Kernel1D<float> roberts_gradient_x; 01212 01213 roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0; 01214 \endcode 01215 01216 The norm is set to the sum of the initialzer values. If the wrong number of 01217 values is given, a run-time error results. It is, however, possible to give 01218 just one initializer. This creates an averaging filter with the given constant: 01219 01220 \code 01221 vigra::Kernel1D<float> average5x1; 01222 01223 average5x1.initExplicitly(-2, 2) = 1.0/5.0; 01224 \endcode 01225 01226 Here, the norm is set to value*size(). 01227 01228 <b> Preconditions:</b> 01229 01230 \code 01231 01232 1. left <= 0 01233 2. right >= 0 01234 3. the number of values in the initializer list 01235 is 1 or equals the size of the kernel. 01236 \endcode 01237 */ 01238 Kernel1D & initExplicitly(int left, int right) 01239 { 01240 vigra_precondition(left <= 0, 01241 "Kernel1D::initExplicitly(): left border must be <= 0."); 01242 vigra_precondition(right >= 0, 01243 "Kernel1D::initExplicitly(): right border must be <= 0."); 01244 01245 right_ = right; 01246 left_ = left; 01247 01248 kernel_.resize(right - left + 1); 01249 01250 return *this; 01251 } 01252 01253 /** Get iterator to center of kernel 01254 01255 Postconditions: 01256 \code 01257 01258 center()[left()] ... center()[right()] are valid kernel positions 01259 \endcode 01260 */ 01261 iterator center() 01262 { 01263 return kernel_.begin() - left(); 01264 } 01265 01266 const_iterator center() const 01267 { 01268 return kernel_.begin() - left(); 01269 } 01270 01271 /** Access kernel value at specified location. 01272 01273 Preconditions: 01274 \code 01275 01276 left() <= location <= right() 01277 \endcode 01278 */ 01279 reference operator[](int location) 01280 { 01281 return kernel_[location - left()]; 01282 } 01283 01284 const_reference operator[](int location) const 01285 { 01286 return kernel_[location - left()]; 01287 } 01288 01289 /** left border of kernel (inclusive), always <= 0 01290 */ 01291 int left() const { return left_; } 01292 01293 /** right border of kernel (inclusive), always >= 0 01294 */ 01295 int right() const { return right_; } 01296 01297 /** size of kernel (right() - left() + 1) 01298 */ 01299 int size() const { return right_ - left_ + 1; } 01300 01301 /** current border treatment mode 01302 */ 01303 BorderTreatmentMode borderTreatment() const 01304 { return border_treatment_; } 01305 01306 /** Set border treatment mode. 01307 */ 01308 void setBorderTreatment( BorderTreatmentMode new_mode) 01309 { border_treatment_ = new_mode; } 01310 01311 /** norm of kernel 01312 */ 01313 value_type norm() const { return norm_; } 01314 01315 /** set a new norm and normalize kernel, use the normalization formula 01316 for the given </tt>derivativeOrder</tt>. 01317 */ 01318 void 01319 normalize(value_type norm, unsigned int derivativeOrder = 0, double offset = 0.0); 01320 01321 /** normalize kernel to norm 1. 01322 */ 01323 void 01324 normalize() 01325 { 01326 normalize(one()); 01327 } 01328 01329 /** get a const accessor 01330 */ 01331 ConstAccessor accessor() const { return ConstAccessor(); } 01332 01333 /** get an accessor 01334 */ 01335 Accessor accessor() { return Accessor(); } 01336 01337 private: 01338 InternalVector kernel_; 01339 int left_, right_; 01340 BorderTreatmentMode border_treatment_; 01341 value_type norm_; 01342 }; 01343 01344 template <class ARITHTYPE> 01345 void Kernel1D<ARITHTYPE>::normalize(value_type norm, 01346 unsigned int derivativeOrder, 01347 double offset) 01348 { 01349 typedef typename NumericTraits<value_type>::RealPromote TmpType; 01350 01351 // find kernel sum 01352 Iterator k = kernel_.begin(); 01353 TmpType sum = NumericTraits<TmpType>::zero(); 01354 01355 if(derivativeOrder == 0) 01356 { 01357 for(; k < kernel_.end(); ++k) 01358 { 01359 sum += *k; 01360 } 01361 } 01362 else 01363 { 01364 unsigned int faculty = 1; 01365 for(unsigned int i = 2; i <= derivativeOrder; ++i) 01366 faculty *= i; 01367 for(double x = left() + offset; k < kernel_.end(); ++x, ++k) 01368 { 01369 sum += *k * VIGRA_CSTD::pow(-x, int(derivativeOrder)) / faculty; 01370 } 01371 } 01372 01373 vigra_precondition(sum != NumericTraits<value_type>::zero(), 01374 "Kernel1D<ARITHTYPE>::normalize(): " 01375 "Cannot normalize a kernel with sum = 0"); 01376 // normalize 01377 sum = norm / sum; 01378 k = kernel_.begin(); 01379 for(; k != kernel_.end(); ++k) 01380 { 01381 *k = *k * sum; 01382 } 01383 01384 norm_ = norm; 01385 } 01386 01387 /***********************************************************************/ 01388 01389 template <class ARITHTYPE> 01390 void Kernel1D<ARITHTYPE>::initGaussian(double std_dev, 01391 value_type norm) 01392 { 01393 vigra_precondition(std_dev >= 0.0, 01394 "Kernel1D::initGaussian(): Standard deviation must be >= 0."); 01395 01396 if(std_dev > 0.0) 01397 { 01398 Gaussian<ARITHTYPE> gauss(std_dev); 01399 01400 // first calculate required kernel sizes 01401 int radius = (int)(3.0 * std_dev + 0.5); 01402 if(radius == 0) 01403 radius = 1; 01404 01405 // allocate the kernel 01406 kernel_.erase(kernel_.begin(), kernel_.end()); 01407 kernel_.reserve(radius*2+1); 01408 01409 for(ARITHTYPE x = -radius; x <= radius; ++x) 01410 { 01411 kernel_.push_back(gauss(x)); 01412 } 01413 left_ = -radius; 01414 right_ = radius; 01415 } 01416 else 01417 { 01418 kernel_.erase(kernel_.begin(), kernel_.end()); 01419 kernel_.push_back(1.0); 01420 left_ = 0; 01421 right_ = 0; 01422 } 01423 01424 if(norm != 0.0) 01425 normalize(norm); 01426 else 01427 norm_ = 1.0; 01428 01429 // best border treatment for Gaussians is BORDER_TREATMENT_CLIP 01430 border_treatment_ = BORDER_TREATMENT_CLIP; 01431 } 01432 01433 /***********************************************************************/ 01434 01435 template <class ARITHTYPE> 01436 void Kernel1D<ARITHTYPE>::initDiscreteGaussian(double std_dev, 01437 value_type norm) 01438 { 01439 vigra_precondition(std_dev >= 0.0, 01440 "Kernel1D::initDiscreteGaussian(): Standard deviation must be >= 0."); 01441 01442 if(std_dev > 0.0) 01443 { 01444 // first calculate required kernel sizes 01445 int radius = (int)(3.0*std_dev + 0.5); 01446 if(radius == 0) 01447 radius = 1; 01448 01449 double f = 2.0 / std_dev / std_dev; 01450 01451 // allocate the working array 01452 int maxIndex = (int)(2.0 * (radius + 5.0 * VIGRA_CSTD::sqrt((double)radius)) + 0.5); 01453 InternalVector warray(maxIndex+1); 01454 warray[maxIndex] = 0.0; 01455 warray[maxIndex-1] = 1.0; 01456 01457 for(int i = maxIndex-2; i >= radius; --i) 01458 { 01459 warray[i] = warray[i+2] + f * (i+1) * warray[i+1]; 01460 if(warray[i] > 1.0e40) 01461 { 01462 warray[i+1] /= warray[i]; 01463 warray[i] = 1.0; 01464 } 01465 } 01466 01467 // the following rescaling ensures that the numbers stay in a sensible range 01468 // during the rest of the iteration, so no other rescaling is needed 01469 double er = VIGRA_CSTD::exp(-radius*radius / (2.0*std_dev*std_dev)); 01470 warray[radius+1] = er * warray[radius+1] / warray[radius]; 01471 warray[radius] = er; 01472 01473 for(int i = radius-1; i >= 0; --i) 01474 { 01475 warray[i] = warray[i+2] + f * (i+1) * warray[i+1]; 01476 er += warray[i]; 01477 } 01478 01479 double scale = norm / (2*er - warray[0]); 01480 01481 initExplicitly(-radius, radius); 01482 iterator c = center(); 01483 01484 for(int i=0; i<=radius; ++i) 01485 { 01486 c[i] = c[-i] = warray[i] * scale; 01487 } 01488 } 01489 else 01490 { 01491 kernel_.erase(kernel_.begin(), kernel_.end()); 01492 kernel_.push_back(norm); 01493 left_ = 0; 01494 right_ = 0; 01495 } 01496 01497 norm_ = norm; 01498 01499 // best border treatment for Gaussians is BORDER_TREATMENT_REFLECT 01500 border_treatment_ = BORDER_TREATMENT_REFLECT; 01501 } 01502 01503 /***********************************************************************/ 01504 01505 template <class ARITHTYPE> 01506 void 01507 Kernel1D<ARITHTYPE>::initGaussianDerivative(double std_dev, 01508 int order, 01509 value_type norm) 01510 { 01511 vigra_precondition(order >= 0, 01512 "Kernel1D::initGaussianDerivative(): Order must be >= 0."); 01513 01514 if(order == 0) 01515 { 01516 initGaussian(std_dev, norm); 01517 return; 01518 } 01519 01520 vigra_precondition(std_dev > 0.0, 01521 "Kernel1D::initGaussianDerivative(): " 01522 "Standard deviation must be > 0."); 01523 01524 Gaussian<ARITHTYPE> gauss(std_dev, order); 01525 01526 // first calculate required kernel sizes 01527 int radius = (int)(3.0 * std_dev + 0.5 * order + 0.5); 01528 if(radius == 0) 01529 radius = 1; 01530 01531 // allocate the kernels 01532 kernel_.clear(); 01533 kernel_.reserve(radius*2+1); 01534 01535 // fill the kernel and calculate the DC component 01536 // introduced by truncation of the Geussian 01537 ARITHTYPE dc = 0.0; 01538 for(ARITHTYPE x = -radius; x <= radius; ++x) 01539 { 01540 kernel_.push_back(gauss(x)); 01541 dc += kernel_[kernel_.size()-1]; 01542 } 01543 dc /= (2.0*radius + 1.0); 01544 01545 // remove DC, but only if kernel correction is permitted by a non-zero 01546 // value for norm 01547 if(norm != 0.0) 01548 { 01549 for(unsigned int i=0; i < kernel_.size(); ++i) 01550 { 01551 kernel_[i] -= dc; 01552 } 01553 } 01554 01555 left_ = -radius; 01556 right_ = radius; 01557 01558 if(norm != 0.0) 01559 normalize(norm, order); 01560 else 01561 norm_ = 1.0; 01562 01563 // best border treatment for Gaussian derivatives is 01564 // BORDER_TREATMENT_REPEAT 01565 border_treatment_ = BORDER_TREATMENT_REPEAT; 01566 } 01567 01568 /***********************************************************************/ 01569 01570 template <class ARITHTYPE> 01571 void 01572 Kernel1D<ARITHTYPE>::initBinomial(int radius, 01573 value_type norm) 01574 { 01575 vigra_precondition(radius > 0, 01576 "Kernel1D::initBinomial(): Radius must be > 0."); 01577 01578 // allocate the kernel 01579 InternalVector kernel(radius*2+1); 01580 01581 int i,j; 01582 for(i=0; i<radius*2+1; ++i) kernel[i] = 0; 01583 01584 // fill kernel 01585 typename InternalVector::iterator x = kernel.begin() + radius; 01586 x[radius] = 1.0; 01587 01588 for(j=radius-1; j>=-radius; --j) 01589 { 01590 for(i=j; i<radius; ++i) 01591 { 01592 x[i] = (x[i] + x[i+1]) / 2.0; 01593 } 01594 x[radius] /= 2.0; 01595 } 01596 01597 // normalize 01598 kernel_.erase(kernel_.begin(), kernel_.end()); 01599 kernel_.reserve(radius*2+1); 01600 01601 for(i=0; i<=radius*2+1; ++i) 01602 { 01603 kernel_.push_back(kernel[i] * norm); 01604 } 01605 01606 left_ = -radius; 01607 right_ = radius; 01608 norm_ = norm; 01609 01610 // best border treatment for Binomial is BORDER_TREATMENT_REFLECT 01611 border_treatment_ = BORDER_TREATMENT_REFLECT; 01612 } 01613 01614 /***********************************************************************/ 01615 01616 template <class ARITHTYPE> 01617 void Kernel1D<ARITHTYPE>::initAveraging(int radius, 01618 value_type norm) 01619 { 01620 vigra_precondition(radius > 0, 01621 "Kernel1D::initAveraging(): Radius must be > 0."); 01622 01623 // calculate scaling 01624 double scale = 1.0 / (radius * 2 + 1); 01625 01626 // normalize 01627 kernel_.erase(kernel_.begin(), kernel_.end()); 01628 kernel_.reserve(radius*2+1); 01629 01630 for(int i=0; i<=radius*2+1; ++i) 01631 { 01632 kernel_.push_back(scale * norm); 01633 } 01634 01635 left_ = -radius; 01636 right_ = radius; 01637 norm_ = norm; 01638 01639 // best border treatment for Averaging is BORDER_TREATMENT_CLIP 01640 border_treatment_ = BORDER_TREATMENT_CLIP; 01641 } 01642 01643 /***********************************************************************/ 01644 01645 template <class ARITHTYPE> 01646 void 01647 Kernel1D<ARITHTYPE>::initSymmetricGradient(value_type norm) 01648 { 01649 kernel_.erase(kernel_.begin(), kernel_.end()); 01650 kernel_.reserve(3); 01651 01652 kernel_.push_back(0.5 * norm); 01653 kernel_.push_back(0.0 * norm); 01654 kernel_.push_back(-0.5 * norm); 01655 01656 left_ = -1; 01657 right_ = 1; 01658 norm_ = norm; 01659 01660 // best border treatment for SymmetricGradient is 01661 // BORDER_TREATMENT_REPEAT 01662 border_treatment_ = BORDER_TREATMENT_REPEAT; 01663 } 01664 01665 /**************************************************************/ 01666 /* */ 01667 /* Argument object factories for Kernel1D */ 01668 /* */ 01669 /* (documentation: see vigra/convolution.hxx) */ 01670 /* */ 01671 /**************************************************************/ 01672 01673 template <class KernelIterator, class KernelAccessor> 01674 inline 01675 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode> 01676 kernel1d(KernelIterator ik, KernelAccessor ka, 01677 int kleft, int kright, BorderTreatmentMode border) 01678 { 01679 return 01680 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>( 01681 ik, ka, kleft, kright, border); 01682 } 01683 01684 template <class T> 01685 inline 01686 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor, 01687 int, int, BorderTreatmentMode> 01688 kernel1d(Kernel1D<T> const & k) 01689 01690 { 01691 return 01692 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor, 01693 int, int, BorderTreatmentMode>( 01694 k.center(), 01695 k.accessor(), 01696 k.left(), k.right(), 01697 k.borderTreatment()); 01698 } 01699 01700 template <class T> 01701 inline 01702 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor, 01703 int, int, BorderTreatmentMode> 01704 kernel1d(Kernel1D<T> const & k, BorderTreatmentMode border) 01705 01706 { 01707 return 01708 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor, 01709 int, int, BorderTreatmentMode>( 01710 k.center(), 01711 k.accessor(), 01712 k.left(), k.right(), 01713 border); 01714 } 01715 01716 01717 } // namespace vigra 01718 01719 #endif // VIGRA_SEPARABLECONVOLUTION_HXX
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|