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