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