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