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