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