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