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