[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]
![]() |
vigra/stdconvolution.hxx | ![]() |
---|
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2002 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* ( Version 1.3.3, Aug 18 2005 ) */ 00008 /* You may use, modify, and distribute this software according */ 00009 /* to the terms stated in the LICENSE file included in */ 00010 /* the VIGRA distribution. */ 00011 /* */ 00012 /* The VIGRA Website is */ 00013 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00014 /* Please direct questions, bug reports, and contributions to */ 00015 /* koethe@informatik.uni-hamburg.de */ 00016 /* */ 00017 /* THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR */ 00018 /* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED */ 00019 /* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */ 00020 /* */ 00021 /************************************************************************/ 00022 00023 00024 #ifndef VIGRA_STDCONVOLUTION_HXX 00025 #define VIGRA_STDCONVOLUTION_HXX 00026 00027 #include <cmath> 00028 #include "vigra/stdimage.hxx" 00029 #include "vigra/bordertreatment.hxx" 00030 #include "vigra/separableconvolution.hxx" 00031 #include "vigra/utilities.hxx" 00032 00033 namespace vigra { 00034 00035 template <class SrcIterator, class SrcAccessor, 00036 class DestIterator, class DestAccessor, 00037 class KernelIterator, class KernelAccessor, 00038 class KSumType> 00039 void internalPixelEvaluationByClip(int x, int y, int w, int h, SrcIterator xs, 00040 SrcAccessor src_acc, DestIterator xd, DestAccessor dest_acc, 00041 KernelIterator ki, Diff2D kul, Diff2D klr, KernelAccessor ak, 00042 KSumType norm) 00043 { 00044 typedef typename 00045 NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType; 00046 typedef 00047 NumericTraits<typename DestAccessor::value_type> DestTraits; 00048 00049 // calculate width and height of the kernel 00050 int kernel_width = klr.x - kul.x + 1; 00051 int kernel_height = klr.y - kul.y + 1; 00052 00053 SumType sum = NumericTraits<SumType>::zero(); 00054 int xx, yy; 00055 int x0, y0, x1, y1; 00056 00057 y0 = (y<klr.y) ? -y : -klr.y; 00058 y1 = (h-y-1<-kul.y) ? h-y-1 : -kul.y; 00059 00060 x0 = (x<klr.x) ? -x : -klr.x; 00061 x1 = (w-x-1<-kul.x) ? w-x-1 : -kul.x; 00062 00063 SrcIterator yys = xs + Diff2D(x0, y0); 00064 KernelIterator yk = ki - Diff2D(x0, y0); 00065 00066 KSumType ksum = NumericTraits<KSumType>::zero(); 00067 kernel_width = x1 - x0 + 1; 00068 kernel_height = y1 - y0 + 1; 00069 00070 //es wird zuerst abgeschnitten und dann gespigelt! 00071 00072 for(yy=0; yy<kernel_height; ++yy, ++yys.y, --yk.y) 00073 { 00074 SrcIterator xxs = yys; 00075 KernelIterator xk = yk; 00076 00077 for(xx=0; xx<kernel_width; ++xx, ++xxs.x, --xk.x) 00078 { 00079 sum += ak(xk) * src_acc(xxs); 00080 ksum += ak(xk); 00081 } 00082 } 00083 00084 // store average in destination pixel 00085 dest_acc.set(DestTraits::fromRealPromote((norm / ksum) * sum), xd); 00086 00087 } 00088 00089 00090 #if 0 00091 00092 template <class SrcIterator, class SrcAccessor, 00093 class DestIterator, class DestAccessor, 00094 class KernelIterator, class KernelAccessor> 00095 void internalPixelEvaluationByWrapReflectRepeat(int x, int y, int src_width, int src_height, SrcIterator xs, 00096 SrcAccessor src_acc, DestIterator xd, DestAccessor dest_acc, 00097 KernelIterator ki, Diff2D kul, Diff2D klr, KernelAccessor ak, 00098 BorderTreatmentMode border) 00099 { 00100 00101 typedef typename 00102 NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType; 00103 typedef 00104 NumericTraits<typename DestAccessor::value_type> DestTraits; 00105 00106 SumType sum = NumericTraits<SumType>::zero(); 00107 00108 SrcIterator src_ul = xs - Diff2D(x, y); 00109 SrcIterator src_lr = src_ul + Diff2D(src_width, src_height); 00110 00111 SrcIterator yys = xs; 00112 KernelIterator yk = ki; 00113 00114 // calculate width and height of the kernel 00115 int kernel_width = klr.x - kul.x + 1; 00116 int kernel_height = klr.y - kul.y + 1; 00117 00118 //Zeigt an wo der Kernel ’ber die Grenzen hinausgeht 00119 bool top_to_much = (y<klr.y) ? true : false; 00120 bool down_to_much = (src_height-y-1<-kul.y)? true : false; 00121 bool left_to_much = (x<klr.x)? true : false; 00122 bool right_to_much = (src_width-x-1<-kul.x)? true : false; 00123 00124 //Die Richtung x und y !!! 00125 //in der bei der Iteration ’ber das aktuelle Bereich im Bild 00126 //iteriert wird. Also wenn von ur->ll dann (-1, +1) und wenn lr->ul 00127 //dann (-1, -1). 00128 Diff2D way_increment; 00129 00130 /* iteriert wird immer aus dem g’ltigen in den ung’ltigen 00131 Bereich! dieser Tupel setzt sich wie folgt zusammen: 00132 1. Wird bei der Iteration in X-Richtung ung’ltiger Bereich 00133 erreicht so wird mit border_increment.first gesprungen und 00134 mit border_increment.third weiter iteriert. 00135 2. Wird bei der Iteration in Y-Richtung ung’ltiger Bereich 00136 erreicht so wird mit border_increment.second gesprungen und 00137 mit border_increment.fourth weiter iteriert. 00138 */ 00139 tuple4<int, int, int, int> border_increment; 00140 if (border == BORDER_TREATMENT_REPEAT){ 00141 border_increment = tuple4<int, int, int, int>(1, 1, 0, 0); 00142 }else if (border == BORDER_TREATMENT_REFLECT){ 00143 border_increment = tuple4<int, int, int, int>(2, 2, -1, -1); 00144 }else{ // BORDER_TREATMENT_WRAP 00145 border_increment = tuple4<int, int, int, int>(src_width, src_height, 1, 1); 00146 } 00147 00148 pair<int, int> valid_step_count; 00149 00150 if(left_to_much && !top_to_much && !down_to_much) 00151 { 00152 yys += klr; 00153 yk += kul; 00154 way_increment = Diff2D(-1, -1); 00155 border_increment.third = -border_increment.third; 00156 border_increment.fourth = -border_increment.fourth; 00157 valid_step_count = std::make_pair((yys - src_ul).x + 1, kernel_height); 00158 } 00159 else if(top_to_much && !left_to_much && !right_to_much) 00160 { 00161 yys += klr; 00162 yk += kul; 00163 way_increment = Diff2D(-1, -1); 00164 border_increment.third = -border_increment.third; 00165 border_increment.fourth = -border_increment.fourth; 00166 valid_step_count = std::make_pair(kernel_width, (yys - src_ul).y + 1); 00167 } 00168 else if(right_to_much && !top_to_much && !down_to_much) 00169 { 00170 yys += kul; 00171 yk += klr; 00172 way_increment = Diff2D(1, 1); 00173 border_increment.first = -border_increment.first; 00174 border_increment.second = -border_increment.second; 00175 valid_step_count = std::make_pair((src_lr - yys).x, kernel_height); 00176 } 00177 else if(down_to_much && !left_to_much && !right_to_much) 00178 { 00179 yys += kul; 00180 yk += klr; 00181 way_increment = Diff2D(1, 1); 00182 border_increment.first = -border_increment.first; 00183 border_increment.second = -border_increment.second; 00184 valid_step_count = std::make_pair(kernel_width, (src_lr - yys).y); 00185 } 00186 else if(down_to_much && left_to_much) 00187 { 00188 yys += kul + Diff2D(kernel_width - 1, 0); 00189 yk += kul + Diff2D(0, kernel_height - 1); 00190 way_increment = Diff2D(-1, 1); 00191 border_increment.second = -border_increment.second; 00192 border_increment.third = -border_increment.third; 00193 valid_step_count = std::make_pair((yys - src_ul).x + 1, (src_lr - yys).y); 00194 } 00195 else if(down_to_much && right_to_much) 00196 { 00197 yys += kul; 00198 yk += klr; 00199 way_increment = Diff2D(1, 1); 00200 border_increment.first = -border_increment.first; 00201 border_increment.second = -border_increment.second; 00202 valid_step_count = std::make_pair((src_lr - yys).x, (src_lr - yys).y); 00203 } 00204 else if(top_to_much && left_to_much) 00205 { 00206 yys += klr; 00207 yk += kul; 00208 way_increment = Diff2D(-1, -1); 00209 border_increment.third = -border_increment.third; 00210 border_increment.fourth = -border_increment.fourth; 00211 valid_step_count = std::make_pair((yys - src_ul).x + 1, (yys - src_ul).y + 1); 00212 } 00213 else 00214 { //top_to_much && right_to_much 00215 yys += kul + Diff2D(0, kernel_height - 1); 00216 yk += kul + Diff2D(kernel_width - 1, 0); 00217 way_increment = Diff2D(1, -1); 00218 border_increment.first = -border_increment.first; 00219 border_increment.fourth = -border_increment.fourth; 00220 valid_step_count = std::make_pair((src_lr - yys).x, (yys - src_ul).y + 1); 00221 } 00222 00223 int yy = 0, xx; 00224 00225 //laeuft den zul„ssigen Bereich in y-Richtung durch 00226 for(; yy < valid_step_count.second; ++yy, yys.y += way_increment.y, yk.y -= way_increment.y ) 00227 { 00228 SrcIterator xxs = yys; 00229 KernelIterator xk = yk; 00230 00231 //laeuft den zul„ssigen Bereich in x-Richtung durch 00232 for(xx = 0; xx < valid_step_count.first; ++xx, xxs.x += way_increment.x, xk.x -= way_increment.x) 00233 { 00234 sum += ak(xk) * src_acc(xxs); 00235 } 00236 00237 //N„chstes ++xxs.x wuerde in unzul„ssigen Bereich 00238 //bringen => Sprung in zulaessigen Bereich 00239 xxs.x += border_increment.first; 00240 00241 for( ; xx < kernel_width; ++xx, xxs.x += border_increment.third, xk.x -= way_increment.x ) 00242 { 00243 sum += ak(xk) * src_acc(xxs); 00244 } 00245 } 00246 00247 //N„chstes ++yys.y wuerde in unzul„ssigen Bereich 00248 //bringen => Sprung in zulaessigen Bereich 00249 yys.y += border_increment.second; 00250 00251 for( ; yy < kernel_height; ++yy, yys.y += border_increment.third, yk.y -= way_increment.y) 00252 { 00253 SrcIterator xxs = yys; 00254 KernelIterator xk = yk; 00255 00256 for(xx=0; xx < valid_step_count.first; ++xx, xxs.x += way_increment.x, xk.x -= way_increment.x) 00257 { 00258 sum += ak(xk) * src_acc(xxs); 00259 } 00260 00261 //Sprung in den zulaessigen Bereich 00262 xxs.x += border_increment.first; 00263 00264 for( ; xx < kernel_width; ++xx, xxs.x += border_increment.third, xk.x -= way_increment.x ) 00265 { 00266 sum += ak(xk) * src_acc(xxs); 00267 } 00268 } 00269 00270 // store average in destination pixel 00271 dest_acc.set(DestTraits::fromRealPromote(sum), xd); 00272 00273 }// end of internalPixelEvaluationByWrapReflectRepeat 00274 #endif /* #if 0 */ 00275 00276 00277 template <class SrcIterator, class SrcAccessor, 00278 class KernelIterator, class KernelAccessor, 00279 class SumType> 00280 void 00281 internalPixelEvaluationByWrapReflectRepeat(SrcIterator xs, SrcAccessor src_acc, 00282 KernelIterator xk, KernelAccessor ak, 00283 int left, int right, int kleft, int kright, 00284 int borderskipx, int borderinc, SumType & sum) 00285 { 00286 SrcIterator xxs = xs + left; 00287 KernelIterator xxk = xk - left; 00288 00289 for(int xx = left; xx <= right; ++xx, ++xxs, --xxk) 00290 { 00291 sum += ak(xxk) * src_acc(xxs); 00292 } 00293 00294 xxs = xs + left - borderskipx; 00295 xxk = xk - left + 1; 00296 for(int xx = left - 1; xx >= -kright; --xx, xxs -= borderinc, ++xxk) 00297 { 00298 sum += ak(xxk) * src_acc(xxs); 00299 } 00300 00301 xxs = xs + right + borderskipx; 00302 xxk = xk - right - 1; 00303 for(int xx = right + 1; xx <= -kleft; ++xx, xxs += borderinc, --xxk) 00304 { 00305 sum += ak(xxk) * src_acc(xxs); 00306 } 00307 } 00308 00309 00310 /** \addtogroup StandardConvolution Two-dimensional convolution functions 00311 00312 Perform 2D non-separable convolution, with and without ROI mask. 00313 00314 These generic convolution functions implement 00315 the standard 2D convolution operation for images that fit 00316 into the required interface. Arbitrary ROI's are supported 00317 by the mask version of the algorithm. 00318 The functions need a suitable 2D kernel to operate. 00319 */ 00320 //@{ 00321 00322 /** \brief Performs a 2 dimensional convolution of the source image using the given 00323 kernel. 00324 00325 The KernelIterator must point to the center of the kernel, and 00326 the kernel's size is given by its upper left (x and y of distance <= 0) and 00327 lower right (distance >= 0) corners. The image must always be larger than the 00328 kernel. At those positions where the kernel does not completely fit 00329 into the image, the specified \ref BorderTreatmentMode is 00330 applied. You can choice between following BorderTreatmentModes: 00331 <ul> 00332 <li>BORDER_TREATMENT_CLIP</li> 00333 <li>BORDER_TREATMENT_AVOID</li> 00334 <li>BORDER_TREATMENT_WRAP</li> 00335 <li>BORDER_TREATMENT_REFLECT</li> 00336 <li>BORDER_TREATMENT_REPEAT</li> 00337 </ul><br> 00338 The images's pixel type (SrcAccessor::value_type) must be a 00339 linear space over the kernel's value_type (KernelAccessor::value_type), 00340 i.e. addition of source values, multiplication with kernel values, 00341 and NumericTraits must be defined. 00342 The kernel's value_type must be an algebraic field, 00343 i.e. the arithmetic operations (+, -, *, /) and NumericTraits must 00344 be defined. 00345 00346 <b> Declarations:</b> 00347 00348 pass arguments explicitly: 00349 \code 00350 namespace vigra { 00351 template <class SrcIterator, class SrcAccessor, 00352 class DestIterator, class DestAccessor, 00353 class KernelIterator, class KernelAccessor> 00354 void convolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc, 00355 DestIterator dest_ul, DestAccessor dest_acc, 00356 KernelIterator ki, KernelAccessor ak, 00357 Diff2D kul, Diff2D klr, BorderTreatmentMode border); 00358 } 00359 \endcode 00360 00361 00362 use argument objects in conjunction with \ref ArgumentObjectFactories: 00363 \code 00364 namespace vigra { 00365 template <class SrcIterator, class SrcAccessor, 00366 class DestIterator, class DestAccessor, 00367 class KernelIterator, class KernelAccessor> 00368 void convolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00369 pair<DestIterator, DestAccessor> dest, 00370 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, 00371 BorderTreatmentMode> kernel); 00372 } 00373 \endcode 00374 00375 <b> Usage:</b> 00376 00377 <b>\#include</b> "<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>"<br> 00378 Namespace: vigra 00379 00380 00381 \code 00382 vigra::FImage src(w,h), dest(w,h); 00383 ... 00384 00385 // define horizontal Sobel filter 00386 vigra::Kernel2D<float> sobel; 00387 00388 sobel.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = // upper left and lower right 00389 0.125, 0.0, -0.125, 00390 0.25, 0.0, -0.25, 00391 0.125, 0.0, -0.125; 00392 00393 vigra::convolveImage(srcImageRange(src), destImage(dest), kernel2d(sobel)); 00394 \endcode 00395 00396 <b> Required Interface:</b> 00397 00398 \code 00399 ImageIterator src_ul, src_lr; 00400 ImageIterator dest_ul; 00401 ImageIterator ik; 00402 00403 SrcAccessor src_accessor; 00404 DestAccessor dest_accessor; 00405 KernelAccessor kernel_accessor; 00406 00407 NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(src_ul); 00408 00409 s = s + s; 00410 s = kernel_accessor(ik) * s; 00411 s -= s; 00412 00413 dest_accessor.set( 00414 NumericTraits<DestAccessor::value_type>::fromRealPromote(s), dest_ul); 00415 00416 NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik); 00417 00418 k += k; 00419 k -= k; 00420 k = k / k; 00421 00422 \endcode 00423 00424 <b> Preconditions:</b> 00425 00426 \code 00427 kul.x <= 0 00428 kul.y <= 0 00429 klr.x >= 0 00430 klr.y >= 0 00431 src_lr.x - src_ul.x >= klr.x + kul.x + 1 00432 src_lr.y - src_ul.y >= klr.y + kul.y + 1 00433 \endcode 00434 00435 If border == BORDER_TREATMENT_CLIP: Sum of kernel elements must be 00436 != 0. 00437 00438 */ 00439 template <class SrcIterator, class SrcAccessor, 00440 class DestIterator, class DestAccessor, 00441 class KernelIterator, class KernelAccessor> 00442 void convolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc, 00443 DestIterator dest_ul, DestAccessor dest_acc, 00444 KernelIterator ki, KernelAccessor ak, 00445 Diff2D kul, Diff2D klr, BorderTreatmentMode border) 00446 { 00447 vigra_precondition((border == BORDER_TREATMENT_CLIP || 00448 border == BORDER_TREATMENT_AVOID || 00449 border == BORDER_TREATMENT_REFLECT || 00450 border == BORDER_TREATMENT_REPEAT || 00451 border == BORDER_TREATMENT_WRAP), 00452 "convolveImage():\n" 00453 " Border treatment must be one of follow treatments:\n" 00454 " - BORDER_TREATMENT_CLIP\n" 00455 " - BORDER_TREATMENT_AVOID\n" 00456 " - BORDER_TREATMENT_REFLECT\n" 00457 " - BORDER_TREATMENT_REPEAT\n" 00458 " - BORDER_TREATMENT_WRAP\n"); 00459 00460 vigra_precondition(kul.x <= 0 && kul.y <= 0, 00461 "convolveImage(): coordinates of " 00462 "kernel's upper left must be <= 0."); 00463 vigra_precondition(klr.x >= 0 && klr.y >= 0, 00464 "convolveImage(): coordinates of " 00465 "kernel's lower right must be >= 0."); 00466 00467 // use traits to determine SumType as to prevent possible overflow 00468 typedef typename 00469 NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType; 00470 typedef typename 00471 NumericTraits<typename KernelAccessor::value_type>::RealPromote KernelSumType; 00472 typedef 00473 NumericTraits<typename DestAccessor::value_type> DestTraits; 00474 00475 // calculate width and height of the image 00476 int w = src_lr.x - src_ul.x; 00477 int h = src_lr.y - src_ul.y; 00478 00479 // calculate width and height of the kernel 00480 int kernel_width = klr.x - kul.x + 1; 00481 int kernel_height = klr.y - kul.y + 1; 00482 00483 vigra_precondition(w >= kernel_width && h >= kernel_height, 00484 "convolveImage(): kernel larger than image."); 00485 00486 int x,y; 00487 00488 KernelSumType norm = NumericTraits<KernelSumType>::zero(); 00489 if(border == BORDER_TREATMENT_CLIP) 00490 { 00491 // caluclate the sum of the kernel elements for renormalization 00492 KernelIterator yk = ki + klr; 00493 00494 //Die Summe der Punkte im Kernel wird ermittelt (= norm) 00495 for(y=0; y<kernel_height; ++y, --yk.y) 00496 { 00497 KernelIterator xk = yk; 00498 for(x=0; x<kernel_width; ++x, --xk.x) 00499 { 00500 norm += ak(xk); 00501 } 00502 } 00503 vigra_precondition(norm != NumericTraits<KernelSumType>::zero(), 00504 "convolveImage(): Cannot use BORDER_TREATMENT_CLIP with a DC-free kernel"); 00505 } 00506 00507 // create iterators for the interior part of the image (where the kernel always fits into the image) 00508 DestIterator yd = dest_ul + Diff2D(klr.x, klr.y); 00509 SrcIterator ys = src_ul + Diff2D(klr.x, klr.y); 00510 SrcIterator send = src_lr + Diff2D(kul.x, kul.y); 00511 00512 // iterate over the interior part 00513 for(; ys.y < send.y; ++ys.y, ++yd.y) 00514 { 00515 // create x iterators 00516 DestIterator xd(yd); 00517 SrcIterator xs(ys); 00518 00519 for(; xs.x < send.x; ++x, ++xs.x, ++xd.x) 00520 { 00521 // init the sum 00522 SumType sum = NumericTraits<SumType>::zero(); 00523 00524 SrcIterator yys = xs - klr; 00525 SrcIterator yyend = xs - kul; 00526 KernelIterator yk = ki + klr; 00527 00528 for(; yys.y <= yyend.y; ++yys.y, --yk.y) 00529 { 00530 typename SrcIterator::row_iterator xxs = yys.rowIterator(); 00531 typename SrcIterator::row_iterator xxe = xxs + kernel_width; 00532 typename KernelIterator::row_iterator xk = yk.rowIterator(); 00533 00534 for(; xxs < xxe; ++xxs, --xk) 00535 { 00536 sum += ak(xk) * src_acc(xxs); 00537 } 00538 } 00539 00540 // store convolution result in destination pixel 00541 dest_acc.set(DestTraits::fromRealPromote(sum), xd); 00542 } 00543 } 00544 00545 if(border == BORDER_TREATMENT_AVOID) 00546 return; // skip processing near the border 00547 00548 int interiorskip = w + kul.x - klr.x - 1; 00549 int borderskipx; 00550 int borderskipy; 00551 int borderinc; 00552 if(border == BORDER_TREATMENT_REPEAT) 00553 { 00554 borderskipx = 0; 00555 borderskipy = 0; 00556 borderinc = 0; 00557 } 00558 else if(border == BORDER_TREATMENT_REFLECT) 00559 { 00560 borderskipx = -1; 00561 borderskipy = -1; 00562 borderinc = -1; 00563 } 00564 else if(border == BORDER_TREATMENT_WRAP) 00565 { 00566 borderskipx = -w+1; 00567 borderskipy = -h+1; 00568 borderinc = 1; 00569 } 00570 00571 // create iterators for the entire image 00572 yd = dest_ul; 00573 ys = src_ul; 00574 00575 // go over the entire image (but skip the already computed points in the loop) 00576 for(y=0; y < h; ++y, ++ys.y, ++yd.y) 00577 { 00578 int top = std::max(-klr.y, src_ul.y - ys.y); 00579 int bottom = std::min(-kul.y, src_lr.y - ys.y - 1); 00580 00581 // create x iterators 00582 DestIterator xd(yd); 00583 SrcIterator xs(ys); 00584 00585 for(x=0; x < w; ++x, ++xs.x, ++xd.x) 00586 { 00587 // check if we are away from the border 00588 if(y >= klr.y && y < h+kul.y && x == klr.x) 00589 { 00590 // yes => skip the already computed points 00591 x += interiorskip; 00592 xs.x += interiorskip; 00593 xd.x += interiorskip; 00594 continue; 00595 } 00596 if (border == BORDER_TREATMENT_CLIP) 00597 { 00598 internalPixelEvaluationByClip(x, y, w, h, xs, src_acc, xd, dest_acc, ki, kul, klr, ak, norm); 00599 } 00600 else 00601 { 00602 int left = std::max(-klr.x, src_ul.x - xs.x); 00603 int right = std::min(-kul.x, src_lr.x - xs.x - 1); 00604 00605 // init the sum 00606 SumType sum = NumericTraits<SumType>::zero(); 00607 00608 // create iterators for the part of the kernel that fits into the image 00609 SrcIterator yys = xs + Size2D(0, top); 00610 KernelIterator yk = ki - Size2D(0, top); 00611 00612 int yy; 00613 for(yy = top; yy <= bottom; ++yy, ++yys.y, --yk.y) 00614 { 00615 internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak, 00616 left, right, kul.x, klr.x, borderskipx, borderinc, sum); 00617 } 00618 yys = xs + Size2D(0, top - borderskipy); 00619 yk = ki - Size2D(0, top - 1); 00620 for(yy = top - 1; yy >= -klr.y; --yy, yys.y -= borderinc, ++yk.y) 00621 { 00622 internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak, 00623 left, right, kul.x, klr.x, borderskipx, borderinc, sum); 00624 } 00625 yys = xs + Size2D(0, bottom + borderskipy); 00626 yk = ki - Size2D(0, bottom + 1); 00627 for(yy = bottom + 1; yy <= -kul.y; ++yy, yys.y += borderinc, --yk.y) 00628 { 00629 internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak, 00630 left, right, kul.x, klr.x, borderskipx, borderinc, sum); 00631 } 00632 00633 // store convolution result in destination pixel 00634 dest_acc.set(DestTraits::fromRealPromote(sum), xd); 00635 00636 // internalPixelEvaluationByWrapReflectRepeat(x, y, w, h, xs, src_acc, xd, dest_acc, ki, kul, klr, ak, border); 00637 } 00638 } 00639 } 00640 } 00641 00642 00643 template <class SrcIterator, class SrcAccessor, 00644 class DestIterator, class DestAccessor, 00645 class KernelIterator, class KernelAccessor> 00646 inline 00647 void convolveImage( 00648 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00649 pair<DestIterator, DestAccessor> dest, 00650 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, 00651 BorderTreatmentMode> kernel) 00652 { 00653 convolveImage(src.first, src.second, src.third, 00654 dest.first, dest.second, 00655 kernel.first, kernel.second, kernel.third, 00656 kernel.fourth, kernel.fifth); 00657 } 00658 00659 00660 /** \brief Performs a 2-dimensional normalized convolution, i.e. convolution with a mask image. 00661 00662 This functions computes 00663 <a href ="http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/PIRODDI1/NormConv/NormConv.html">normalized 00664 convolution</a> as defined in 00665 Knutsson, H. and Westin, C-F.: <i>Normalized and differential convolution: 00666 Methods for Interpolation and Filtering of incomplete and uncertain data</i>. 00667 Proc. of the IEEE Conf. on Computer Vision and Pattern Recognition, 1993, 515-523. 00668 00669 The mask image must be binary and encodes which pixels of the original image 00670 are valid. It is used as follows: 00671 Only pixel under the mask are used in the calculations. Whenever a part of the 00672 kernel lies outside the mask, it is ignored, and the kernel is renormalized to its 00673 original norm (analogous to the CLIP \ref BorderTreatmentMode). Thus, a useful convolution 00674 result is computed whenever <i>at least one valid pixel is within the current window</i> 00675 Thus, destination pixels not under the mask still receive a value if they are <i>near</i> 00676 the mask. Therefore, this algorithm is useful as an interpolator of sparse input data. 00677 If you are only interested in the destination values under the mask, you can perform 00678 a subsequent \ref copyImageIf(). 00679 00680 The KernelIterator must point to the center of the kernel, and 00681 the kernel's size is given by its upper left (x and y of distance <= 0) and 00682 lower right (distance >= 0) corners. The image must always be larger than the 00683 kernel. At those positions where the kernel does not completely fit 00684 into the image, the specified \ref BorderTreatmentMode is 00685 applied. Only BORDER_TREATMENT_CLIP and BORDER_TREATMENT_AVOID are currently 00686 supported. 00687 00688 The images's pixel type (SrcAccessor::value_type) must be a 00689 linear space over the kernel's value_type (KernelAccessor::value_type), 00690 i.e. addition of source values, multiplication with kernel values, 00691 and NumericTraits must be defined. 00692 The kernel's value_type must be an algebraic field, 00693 i.e. the arithmetic operations (+, -, *, /) and NumericTraits must 00694 be defined. 00695 00696 <b> Declarations:</b> 00697 00698 pass arguments explicitly: 00699 \code 00700 namespace vigra { 00701 template <class SrcIterator, class SrcAccessor, 00702 class MaskIterator, class MaskAccessor, 00703 class DestIterator, class DestAccessor, 00704 class KernelIterator, class KernelAccessor> 00705 void 00706 normalizedConvolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc, 00707 MaskIterator mul, MaskAccessor am, 00708 DestIterator dest_ul, DestAccessor dest_acc, 00709 KernelIterator ki, KernelAccessor ak, 00710 Diff2D kul, Diff2D klr, BorderTreatmentMode border); 00711 } 00712 \endcode 00713 00714 00715 use argument objects in conjunction with \ref ArgumentObjectFactories: 00716 \code 00717 namespace vigra { 00718 template <class SrcIterator, class SrcAccessor, 00719 class MaskIterator, class MaskAccessor, 00720 class DestIterator, class DestAccessor, 00721 class KernelIterator, class KernelAccessor> 00722 inline 00723 void normalizedConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00724 pair<MaskIterator, MaskAccessor> mask, 00725 pair<DestIterator, DestAccessor> dest, 00726 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, 00727 BorderTreatmentMode> kernel); 00728 } 00729 \endcode 00730 00731 <b> Usage:</b> 00732 00733 <b>\#include</b> "<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>"<br> 00734 Namespace: vigra 00735 00736 00737 \code 00738 vigra::FImage src(w,h), dest(w,h); 00739 vigra::CImage mask(w,h); 00740 ... 00741 00742 // define 3x3 binomial filter 00743 vigra::Kernel2D<float> binom; 00744 00745 binom.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = // upper left and lower right 00746 0.0625, 0.125, 0.0625, 00747 0.125, 0.25, 0.125, 00748 0.0625, 0.125, 0.0625; 00749 00750 vigra::normalizedConvolveImage(srcImageRange(src), maskImage(mask), destImage(dest), kernel2d(binom)); 00751 \endcode 00752 00753 <b> Required Interface:</b> 00754 00755 \code 00756 ImageIterator src_ul, src_lr; 00757 ImageIterator mul; 00758 ImageIterator dest_ul; 00759 ImageIterator ik; 00760 00761 SrcAccessor src_accessor; 00762 MaskAccessor mask_accessor; 00763 DestAccessor dest_accessor; 00764 KernelAccessor kernel_accessor; 00765 00766 NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(src_ul); 00767 00768 s = s + s; 00769 s = kernel_accessor(ik) * s; 00770 s -= s; 00771 00772 if(mask_accessor(mul)) ...; 00773 00774 dest_accessor.set( 00775 NumericTraits<DestAccessor::value_type>::fromRealPromote(s), dest_ul); 00776 00777 NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik); 00778 00779 k += k; 00780 k -= k; 00781 k = k / k; 00782 00783 \endcode 00784 00785 <b> Preconditions:</b> 00786 00787 \code 00788 kul.x <= 0 00789 kul.y <= 0 00790 klr.x >= 0 00791 klr.y >= 0 00792 src_lr.x - src_ul.x >= klr.x + kul.x + 1 00793 src_lr.y - src_ul.y >= klr.y + kul.y + 1 00794 border == BORDER_TREATMENT_CLIP || border == BORDER_TREATMENT_AVOID 00795 \endcode 00796 00797 Sum of kernel elements must be != 0. 00798 00799 */ 00800 template <class SrcIterator, class SrcAccessor, 00801 class DestIterator, class DestAccessor, 00802 class MaskIterator, class MaskAccessor, 00803 class KernelIterator, class KernelAccessor> 00804 void 00805 normalizedConvolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc, 00806 MaskIterator mul, MaskAccessor am, 00807 DestIterator dest_ul, DestAccessor dest_acc, 00808 KernelIterator ki, KernelAccessor ak, 00809 Diff2D kul, Diff2D klr, BorderTreatmentMode border) 00810 { 00811 vigra_precondition((border == BORDER_TREATMENT_CLIP || 00812 border == BORDER_TREATMENT_AVOID), 00813 "normalizedConvolveImage(): " 00814 "Border treatment must be BORDER_TREATMENT_CLIP or BORDER_TREATMENT_AVOID."); 00815 00816 vigra_precondition(kul.x <= 0 && kul.y <= 0, 00817 "normalizedConvolveImage(): left borders must be <= 0."); 00818 vigra_precondition(klr.x >= 0 && klr.y >= 0, 00819 "normalizedConvolveImage(): right borders must be >= 0."); 00820 00821 // use traits to determine SumType as to prevent possible overflow 00822 typedef typename 00823 NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType; 00824 typedef typename 00825 NumericTraits<typename KernelAccessor::value_type>::RealPromote KSumType; 00826 typedef 00827 NumericTraits<typename DestAccessor::value_type> DestTraits; 00828 00829 // calculate width and height of the image 00830 int w = src_lr.x - src_ul.x; 00831 int h = src_lr.y - src_ul.y; 00832 int kernel_width = klr.x - kul.x + 1; 00833 int kernel_height = klr.y - kul.y + 1; 00834 00835 int x,y; 00836 int ystart = (border == BORDER_TREATMENT_AVOID) ? klr.y : 0; 00837 int yend = (border == BORDER_TREATMENT_AVOID) ? h+kul.y : h; 00838 int xstart = (border == BORDER_TREATMENT_AVOID) ? klr.x : 0; 00839 int xend = (border == BORDER_TREATMENT_AVOID) ? w+kul.x : w; 00840 00841 // create y iterators 00842 DestIterator yd = dest_ul + Diff2D(xstart, ystart); 00843 SrcIterator ys = src_ul + Diff2D(xstart, ystart); 00844 MaskIterator ym = mul + Diff2D(xstart, ystart); 00845 00846 KSumType norm = ak(ki); 00847 int xx, yy; 00848 KernelIterator yk = ki + klr; 00849 for(yy=0; yy<kernel_height; ++yy, --yk.y) 00850 { 00851 KernelIterator xk = yk; 00852 00853 for(xx=0; xx<kernel_width; ++xx, --xk.x) 00854 { 00855 norm += ak(xk); 00856 } 00857 } 00858 norm -= ak(ki); 00859 00860 00861 for(y=ystart; y < yend; ++y, ++ys.y, ++yd.y, ++ym.y) 00862 { 00863 // create x iterators 00864 DestIterator xd(yd); 00865 SrcIterator xs(ys); 00866 MaskIterator xm(ym); 00867 00868 for(x=xstart; x < xend; ++x, ++xs.x, ++xd.x, ++xm.x) 00869 { 00870 // how much of the kernel fits into the image ? 00871 int x0, y0, x1, y1; 00872 00873 y0 = (y<klr.y) ? -y : -klr.y; 00874 y1 = (h-y-1<-kul.y) ? h-y-1 : -kul.y; 00875 x0 = (x<klr.x) ? -x : -klr.x; 00876 x1 = (w-x-1<-kul.x) ? w-x-1 : -kul.x; 00877 00878 bool first = true; 00879 // init the sum 00880 SumType sum; 00881 KSumType ksum; 00882 00883 SrcIterator yys = xs + Diff2D(x0, y0); 00884 MaskIterator yym = xm + Diff2D(x0, y0); 00885 KernelIterator yk = ki - Diff2D(x0, y0); 00886 00887 int xx, kernel_width, kernel_height; 00888 kernel_width = x1 - x0 + 1; 00889 kernel_height = y1 - y0 + 1; 00890 for(yy=0; yy<kernel_height; ++yy, ++yys.y, --yk.y, ++yym.y) 00891 { 00892 typename SrcIterator::row_iterator xxs = yys.rowIterator(); 00893 typename SrcIterator::row_iterator xxend = xxs + kernel_width; 00894 typename MaskIterator::row_iterator xxm = yym.rowIterator(); 00895 typename KernelIterator::row_iterator xk = yk.rowIterator(); 00896 00897 for(xx=0; xxs < xxend; ++xxs.x, --xk.x, ++xxm.x) 00898 { 00899 if(!am(xxm)) continue; 00900 00901 if(first) 00902 { 00903 sum = ak(xk) * src_acc(xxs); 00904 ksum = ak(xk); 00905 first = false; 00906 } 00907 else 00908 { 00909 sum += ak(xk) * src_acc(xxs); 00910 ksum += ak(xk); 00911 } 00912 } 00913 } 00914 // store average in destination pixel 00915 if(!first && 00916 ksum != NumericTraits<KSumType>::zero()) 00917 { 00918 dest_acc.set(DestTraits::fromRealPromote((norm / ksum) * sum), xd); 00919 } 00920 } 00921 } 00922 } 00923 00924 00925 template <class SrcIterator, class SrcAccessor, 00926 class DestIterator, class DestAccessor, 00927 class MaskIterator, class MaskAccessor, 00928 class KernelIterator, class KernelAccessor> 00929 inline 00930 void normalizedConvolveImage( 00931 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00932 pair<MaskIterator, MaskAccessor> mask, 00933 pair<DestIterator, DestAccessor> dest, 00934 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, 00935 BorderTreatmentMode> kernel) 00936 { 00937 normalizedConvolveImage(src.first, src.second, src.third, 00938 mask.first, mask.second, 00939 dest.first, dest.second, 00940 kernel.first, kernel.second, kernel.third, 00941 kernel.fourth, kernel.fifth); 00942 } 00943 00944 /** \brief Deprecated name of 2-dimensional normalized convolution, i.e. convolution with a mask image. 00945 00946 See \ref normalizedConvolveImage() for documentation. 00947 00948 <b> Declarations:</b> 00949 00950 pass arguments explicitly: 00951 \code 00952 namespace vigra { 00953 template <class SrcIterator, class SrcAccessor, 00954 class MaskIterator, class MaskAccessor, 00955 class DestIterator, class DestAccessor, 00956 class KernelIterator, class KernelAccessor> 00957 void 00958 convolveImageWithMask(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc, 00959 MaskIterator mul, MaskAccessor am, 00960 DestIterator dest_ul, DestAccessor dest_acc, 00961 KernelIterator ki, KernelAccessor ak, 00962 Diff2D kul, Diff2D klr, BorderTreatmentMode border); 00963 } 00964 \endcode 00965 00966 00967 use argument objects in conjunction with \ref ArgumentObjectFactories: 00968 \code 00969 namespace vigra { 00970 template <class SrcIterator, class SrcAccessor, 00971 class MaskIterator, class MaskAccessor, 00972 class DestIterator, class DestAccessor, 00973 class KernelIterator, class KernelAccessor> 00974 inline 00975 void convolveImageWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00976 pair<MaskIterator, MaskAccessor> mask, 00977 pair<DestIterator, DestAccessor> dest, 00978 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, 00979 BorderTreatmentMode> kernel); 00980 } 00981 \endcode 00982 */ 00983 template <class SrcIterator, class SrcAccessor, 00984 class DestIterator, class DestAccessor, 00985 class MaskIterator, class MaskAccessor, 00986 class KernelIterator, class KernelAccessor> 00987 inline void 00988 convolveImageWithMask(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc, 00989 MaskIterator mul, MaskAccessor am, 00990 DestIterator dest_ul, DestAccessor dest_acc, 00991 KernelIterator ki, KernelAccessor ak, 00992 Diff2D kul, Diff2D klr, BorderTreatmentMode border) 00993 { 00994 normalizedConvolveImage(src_ul, src_lr, src_acc, 00995 mul, am, 00996 dest_ul, dest_acc, 00997 ki, ak, kul, klr, border); 00998 } 00999 01000 template <class SrcIterator, class SrcAccessor, 01001 class DestIterator, class DestAccessor, 01002 class MaskIterator, class MaskAccessor, 01003 class KernelIterator, class KernelAccessor> 01004 inline 01005 void convolveImageWithMask( 01006 triple<SrcIterator, SrcIterator, SrcAccessor> src, 01007 pair<MaskIterator, MaskAccessor> mask, 01008 pair<DestIterator, DestAccessor> dest, 01009 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, 01010 BorderTreatmentMode> kernel) 01011 { 01012 normalizedConvolveImage(src.first, src.second, src.third, 01013 mask.first, mask.second, 01014 dest.first, dest.second, 01015 kernel.first, kernel.second, kernel.third, 01016 kernel.fourth, kernel.fifth); 01017 } 01018 01019 //@} 01020 01021 /********************************************************/ 01022 /* */ 01023 /* Kernel2D */ 01024 /* */ 01025 /********************************************************/ 01026 01027 /** \brief Generic 2 dimensional convolution kernel. 01028 01029 This kernel may be used for convolution of 2 dimensional signals. 01030 01031 Convolution functions access the kernel via an ImageIterator 01032 which they get by calling \ref center(). This iterator 01033 points to the center of the kernel. The kernel's size is given by its upperLeft() 01034 (upperLeft().x <= 0, upperLeft().y <= 0) 01035 and lowerRight() (lowerRight().x >= 0, lowerRight().y >= 0) methods. 01036 The desired border treatment mode is returned by borderTreatment(). 01037 (Note that the \ref StandardConvolution "2D convolution functions" don't currently 01038 support all modes.) 01039 01040 The different init functions create a kernel with the specified 01041 properties. The requirements for the kernel's value_type depend 01042 on the init function used. At least NumericTraits must be defined. 01043 01044 The kernel defines a factory function kernel2d() to create an argument object 01045 (see \ref KernelArgumentObjectFactories). 01046 01047 <b> Usage:</b> 01048 01049 <b>\#include</b> "<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>"<br> 01050 Namespace: vigra 01051 01052 \code 01053 vigra::FImage src(w,h), dest(w,h); 01054 ... 01055 01056 // define horizontal Sobel filter 01057 vigra::Kernel2D<float> sobel; 01058 01059 sobel.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = // upper left and lower right 01060 0.125, 0.0, -0.125, 01061 0.25, 0.0, -0.25, 01062 0.125, 0.0, -0.125; 01063 01064 vigra::convolveImage(srcImageRange(src), destImage(dest), kernel2d(sobel)); 01065 \endcode 01066 01067 <b> Required Interface:</b> 01068 01069 \code 01070 value_type v = NumericTraits<value_type>::one(); 01071 \endcode 01072 01073 See also the init functions. 01074 */ 01075 template <class ARITHTYPE> 01076 class Kernel2D 01077 { 01078 public: 01079 /** the kernel's value type 01080 */ 01081 typedef ARITHTYPE value_type; 01082 01083 /** 2D random access iterator over the kernel's values 01084 */ 01085 typedef typename BasicImage<value_type>::traverser Iterator; 01086 01087 /** const 2D random access iterator over the kernel's values 01088 */ 01089 typedef typename BasicImage<value_type>::const_traverser ConstIterator; 01090 01091 /** the kernel's accessor 01092 */ 01093 typedef typename BasicImage<value_type>::Accessor Accessor; 01094 01095 /** the kernel's const accessor 01096 */ 01097 typedef typename BasicImage<value_type>::ConstAccessor ConstAccessor; 01098 01099 struct InitProxy 01100 { 01101 typedef typename 01102 BasicImage<value_type>::ScanOrderIterator Iterator; 01103 01104 InitProxy(Iterator i, int count, value_type & norm) 01105 : iter_(i), base_(i), 01106 count_(count), sum_(count), 01107 norm_(norm) 01108 {} 01109 01110 ~InitProxy() 01111 { 01112 vigra_precondition(count_ == 1 || count_ == sum_, 01113 "Kernel2D::initExplicitly(): " 01114 "Too few init values."); 01115 } 01116 01117 InitProxy & operator,(value_type const & v) 01118 { 01119 if(count_ == sum_) norm_ = *iter_; 01120 01121 --count_; 01122 vigra_precondition(count_ > 0, 01123 "Kernel2D::initExplicitly(): " 01124 "Too many init values."); 01125 01126 norm_ += v; 01127 01128 ++iter_; 01129 *iter_ = v; 01130 01131 return *this; 01132 } 01133 01134 Iterator iter_, base_; 01135 int count_, sum_; 01136 value_type & norm_; 01137 }; 01138 01139 static value_type one() { return NumericTraits<value_type>::one(); } 01140 01141 /** Default constructor. 01142 Creates a kernel of size 1x1 which would copy the signal 01143 unchanged. 01144 */ 01145 Kernel2D() 01146 : kernel_(1, 1, one()), 01147 left_(0, 0), 01148 right_(0, 0), 01149 norm_(one()), 01150 border_treatment_(BORDER_TREATMENT_CLIP) 01151 {} 01152 01153 /** Copy constructor. 01154 */ 01155 Kernel2D(Kernel2D const & k) 01156 : kernel_(k.kernel_), 01157 left_(k.left_), 01158 right_(k.right_), 01159 norm_(k.norm_), 01160 border_treatment_(k.border_treatment_) 01161 {} 01162 01163 /** Copy assignment. 01164 */ 01165 Kernel2D & operator=(Kernel2D const & k) 01166 { 01167 if(this != &k) 01168 { 01169 kernel_ = k.kernel_; 01170 left_ = k.left_; 01171 right_ = k.right_; 01172 norm_ = k.norm_; 01173 border_treatment_ = k.border_treatment_; 01174 } 01175 return *this; 01176 } 01177 01178 /** Initialization. 01179 This initializes the kernel with the given constant. The norm becomes 01180 v*width()*height(). 01181 01182 Instead of a single value an initializer list of length width()*height() 01183 can be used like this: 01184 01185 \code 01186 vigra::Kernel2D<float> binom; 01187 01188 binom.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = 01189 0.0625, 0.125, 0.0625, 01190 0.125, 0.25, 0.125, 01191 0.0625, 0.125, 0.0625; 01192 \endcode 01193 01194 In this case, the norm will be set to the sum of the init values. 01195 An initializer list of wrong length will result in a run-time error. 01196 */ 01197 InitProxy operator=(value_type const & v) 01198 { 01199 int size = (right_.x - left_.x + 1) * 01200 (right_.y - left_.y + 1); 01201 kernel_ = v; 01202 norm_ = (double)size*v; 01203 01204 return InitProxy(kernel_.begin(), size, norm_); 01205 } 01206 01207 /** Destructor. 01208 */ 01209 ~Kernel2D() 01210 {} 01211 01212 /** Init the 2D kernel as the cartesian product of two 1D kernels 01213 of type \ref Kernel1D. The norm becomes the product of the two original 01214 norms. 01215 01216 <b> Required Interface:</b> 01217 01218 The kernel's value_type must be a linear algebra. 01219 01220 \code 01221 vigra::Kernel2D<...>::value_type v; 01222 v = v * v; 01223 \endcode 01224 */ 01225 void initSeparable(Kernel1D<value_type> & kx, 01226 Kernel1D<value_type> & ky) 01227 { 01228 left_ = Diff2D(kx.left(), ky.left()); 01229 right_ = Diff2D(kx.right(), ky.right()); 01230 int w = right_.x - left_.x + 1; 01231 int h = right_.y - left_.y + 1; 01232 kernel_.resize(w, h); 01233 01234 norm_ = kx.norm() * ky.norm(); 01235 01236 typedef typename Kernel1D<value_type>::Iterator KIter; 01237 typename Kernel1D<value_type>::Accessor ka; 01238 01239 KIter kiy = ky.center() + left_.y; 01240 Iterator iy = center() + left_; 01241 01242 for(int y=left_.y; y<=right_.y; ++y, ++kiy, ++iy.y) 01243 { 01244 KIter kix = kx.center() + left_.x; 01245 Iterator ix = iy; 01246 for(int x=left_.x; x<=right_.x; ++x, ++kix, ++ix.x) 01247 { 01248 *ix = ka(kix) * ka(kiy); 01249 } 01250 } 01251 } 01252 01253 /** Init the 2D kernel as the cartesian product of two 1D kernels 01254 given explicitly by iterators and sizes. The norm becomes the 01255 sum of the resulting kernel values. 01256 01257 <b> Required Interface:</b> 01258 01259 The kernel's value_type must be a linear algebra. 01260 01261 \code 01262 vigra::Kernel2D<...>::value_type v; 01263 v = v * v; 01264 v += v; 01265 \endcode 01266 01267 <b> Preconditions:</b> 01268 01269 \code 01270 xleft <= 0; 01271 xright >= 0; 01272 yleft <= 0; 01273 yright >= 0; 01274 \endcode 01275 */ 01276 template <class KernelIterator> 01277 void initSeparable(KernelIterator kxcenter, int xleft, int xright, 01278 KernelIterator kycenter, int yleft, int yright) 01279 { 01280 vigra_precondition(xleft <= 0 && yleft <= 0, 01281 "Kernel2D::initSeparable(): left borders must be <= 0."); 01282 vigra_precondition(xright >= 0 && yright >= 0, 01283 "Kernel2D::initSeparable(): right borders must be >= 0."); 01284 01285 left_ = Point2D(xleft, yleft); 01286 right_ = Point2D(xright, yright); 01287 01288 int w = right_.x - left_.x + 1; 01289 int h = right_.y - left_.y + 1; 01290 kernel_.resize(w, h); 01291 01292 KernelIterator kiy = kycenter + left_.y; 01293 Iterator iy = center() + left_; 01294 01295 for(int y=left_.y; y<=right_.y; ++y, ++kiy, ++iy.y) 01296 { 01297 KernelIterator kix = kxcenter + left_.x; 01298 Iterator ix = iy; 01299 for(int x=left_.x; x<=right_.x; ++x, ++kix, ++ix.x) 01300 { 01301 *ix = *kix * *kiy; 01302 } 01303 } 01304 01305 typename BasicImage<value_type>::iterator i = kernel_.begin(); 01306 typename BasicImage<value_type>::iterator iend = kernel_.end(); 01307 norm_ = *i; 01308 ++i; 01309 01310 for(; i!= iend; ++i) 01311 { 01312 norm_ += *i; 01313 } 01314 } 01315 01316 /** Init the 2D kernel as a circular averaging filter. The norm will be 01317 calculated as 01318 <TT>NumericTraits<value_type>::one() / (number of non-zero kernel values)</TT>. 01319 The kernel's value_type must be a linear space. 01320 01321 <b> Required Interface:</b> 01322 01323 \code 01324 value_type v = vigra::NumericTraits<value_type>::one(); 01325 01326 double d; 01327 v = d * v; 01328 \endcode 01329 01330 <b> Precondition:</b> 01331 01332 \code 01333 radius > 0; 01334 \endcode 01335 */ 01336 void initDisk(int radius) 01337 { 01338 vigra_precondition(radius > 0, 01339 "Kernel2D::initDisk(): radius must be > 0."); 01340 01341 left_ = Point2D(-radius, -radius); 01342 right_ = Point2D(radius, radius); 01343 int w = right_.x - left_.x + 1; 01344 int h = right_.y - left_.y + 1; 01345 kernel_.resize(w, h); 01346 norm_ = NumericTraits<value_type>::one(); 01347 01348 kernel_ = NumericTraits<value_type>::zero(); 01349 double count = 0.0; 01350 01351 Iterator k = center(); 01352 double r2 = (double)radius*radius; 01353 01354 int i; 01355 for(i=0; i<= radius; ++i) 01356 { 01357 double r = (double) i - 0.5; 01358 int w = (int)(VIGRA_CSTD::sqrt(r2 - r*r) + 0.5); 01359 for(int j=-w; j<=w; ++j) 01360 { 01361 k(j, i) = NumericTraits<value_type>::one(); 01362 k(j, -i) = NumericTraits<value_type>::one(); 01363 count += (i != 0) ? 2.0 : 1.0; 01364 } 01365 } 01366 01367 count = 1.0 / count; 01368 01369 for(int y=-radius; y<=radius; ++y) 01370 { 01371 for(int x=-radius; x<=radius; ++x) 01372 { 01373 k(x,y) = count * k(x,y); 01374 } 01375 } 01376 } 01377 01378 /** Init the kernel by an explicit initializer list. 01379 The upper left and lower right corners of the kernel must be passed. 01380 A comma-separated initializer list is given after the assignment operator. 01381 This function is used like this: 01382 01383 \code 01384 // define horizontal Sobel filter 01385 vigra::Kernel2D<float> sobel; 01386 01387 sobel.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = 01388 0.125, 0.0, -0.125, 01389 0.25, 0.0, -0.25, 01390 0.125, 0.0, -0.125; 01391 \endcode 01392 01393 The norm is set to the sum of the initialzer values. If the wrong number of 01394 values is given, a run-time error results. It is, however, possible to give 01395 just one initializer. This creates an averaging filter with the given constant: 01396 01397 \code 01398 vigra::Kernel2D<float> average3x3; 01399 01400 average3x3.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = 1.0/9.0; 01401 \endcode 01402 01403 Here, the norm is set to value*width()*height(). 01404 01405 <b> Preconditions:</b> 01406 01407 \code 01408 1. upperleft.x <= 0; 01409 2. upperleft.y <= 0; 01410 3. lowerright.x >= 0; 01411 4. lowerright.y >= 0; 01412 5. the number of values in the initializer list 01413 is 1 or equals the size of the kernel. 01414 \endcode 01415 */ 01416 Kernel2D & initExplicitly(Diff2D upperleft, Diff2D lowerright) 01417 { 01418 vigra_precondition(upperleft.x <= 0 && upperleft.y <= 0, 01419 "Kernel2D::initExplicitly(): left borders must be <= 0."); 01420 vigra_precondition(lowerright.x >= 0 && lowerright.y >= 0, 01421 "Kernel2D::initExplicitly(): right borders must be >= 0."); 01422 01423 left_ = Point2D(upperleft); 01424 right_ = Point2D(lowerright); 01425 01426 int w = right_.x - left_.x + 1; 01427 int h = right_.y - left_.y + 1; 01428 kernel_.resize(w, h); 01429 01430 return *this; 01431 } 01432 01433 /** Coordinates of the upper left corner of the kernel. 01434 */ 01435 Point2D upperLeft() const { return left_; } 01436 01437 /** Coordinates of the lower right corner of the kernel. 01438 */ 01439 Point2D lowerRight() const { return right_; } 01440 01441 /** Width of the kernel. 01442 */ 01443 int width() const { return right_.x - left_.x + 1; } 01444 01445 /** Height of the kernel. 01446 */ 01447 int height() const { return right_.y - left_.y + 1; } 01448 01449 /** ImageIterator that points to the center of the kernel (coordinate (0,0)). 01450 */ 01451 Iterator center() { return kernel_.upperLeft() - left_; } 01452 01453 /** ImageIterator that points to the center of the kernel (coordinate (0,0)). 01454 */ 01455 ConstIterator center() const { return kernel_.upperLeft() - left_; } 01456 01457 /** Access kernel entry at given position. 01458 */ 01459 value_type & operator()(int x, int y) 01460 { return kernel_[Diff2D(x,y) - left_]; } 01461 01462 /** Read kernel entry at given position. 01463 */ 01464 value_type operator()(int x, int y) const 01465 { return kernel_[Diff2D(x,y) - left_]; } 01466 01467 /** Access kernel entry at given position. 01468 */ 01469 value_type & operator[](Diff2D const & d) 01470 { return kernel_[d - left_]; } 01471 01472 /** Read kernel entry at given position. 01473 */ 01474 value_type operator[](Diff2D const & d) const 01475 { return kernel_[d - left_]; } 01476 01477 /** Norm of the kernel (i.e. sum of its elements). 01478 */ 01479 value_type norm() const { return norm_; } 01480 01481 /** The kernels default accessor. 01482 */ 01483 Accessor accessor() { return Accessor(); } 01484 01485 /** The kernels default const accessor. 01486 */ 01487 ConstAccessor accessor() const { return ConstAccessor(); } 01488 01489 /** Normalize the kernel to the given value. (The norm is the sum of all kernel 01490 elements.) The kernel's value_type must be a division algebra or 01491 algebraic field. 01492 01493 <b> Required Interface:</b> 01494 01495 \code 01496 value_type v = vigra::NumericTraits<value_type>::one(); // if norm is not 01497 // given explicitly 01498 01499 v += v; 01500 v = v * v; 01501 v = v / v; 01502 \endcode 01503 */ 01504 void normalize(value_type norm) 01505 { 01506 typename BasicImage<value_type>::iterator i = kernel_.begin(); 01507 typename BasicImage<value_type>::iterator iend = kernel_.end(); 01508 typename NumericTraits<value_type>::RealPromote sum = *i; 01509 ++i; 01510 01511 for(; i!= iend; ++i) 01512 { 01513 sum += *i; 01514 } 01515 01516 sum = norm / sum; 01517 i = kernel_.begin(); 01518 for(; i != iend; ++i) 01519 { 01520 *i = *i * sum; 01521 } 01522 01523 norm_ = norm; 01524 } 01525 01526 /** Normalize the kernel to norm 1. 01527 */ 01528 void normalize() 01529 { 01530 normalize(one()); 01531 } 01532 01533 /** current border treatment mode 01534 */ 01535 BorderTreatmentMode borderTreatment() const 01536 { return border_treatment_; } 01537 01538 /** Set border treatment mode. 01539 Only <TT>BORDER_TREATMENT_CLIP</TT> and <TT>BORDER_TREATMENT_AVOID</TT> are currently 01540 allowed. 01541 */ 01542 void setBorderTreatment( BorderTreatmentMode new_mode) 01543 { 01544 vigra_precondition((new_mode == BORDER_TREATMENT_CLIP || 01545 new_mode == BORDER_TREATMENT_AVOID || 01546 new_mode == BORDER_TREATMENT_REFLECT || 01547 new_mode == BORDER_TREATMENT_REPEAT || 01548 new_mode == BORDER_TREATMENT_WRAP), 01549 "convolveImage():\n" 01550 " Border treatment must be one of follow treatments:\n" 01551 " - BORDER_TREATMENT_CLIP\n" 01552 " - BORDER_TREATMENT_AVOID\n" 01553 " - BORDER_TREATMENT_REFLECT\n" 01554 " - BORDER_TREATMENT_REPEAT\n" 01555 " - BORDER_TREATMENT_WRAP\n"); 01556 01557 border_treatment_ = new_mode; 01558 } 01559 01560 01561 private: 01562 BasicImage<value_type> kernel_; 01563 Point2D left_, right_; 01564 value_type norm_; 01565 BorderTreatmentMode border_treatment_; 01566 }; 01567 01568 /**************************************************************/ 01569 /* */ 01570 /* Argument object factories for Kernel2D */ 01571 /* */ 01572 /* (documentation: see vigra/convolution.hxx) */ 01573 /* */ 01574 /**************************************************************/ 01575 01576 template <class KernelIterator, class KernelAccessor> 01577 inline 01578 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, BorderTreatmentMode> 01579 kernel2d(KernelIterator ik, KernelAccessor ak, Diff2D kul, Diff2D klr, 01580 BorderTreatmentMode border) 01581 01582 { 01583 return 01584 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, BorderTreatmentMode> ( 01585 ik, ak, kul, klr, border); 01586 } 01587 01588 template <class T> 01589 inline 01590 tuple5<typename Kernel2D<T>::ConstIterator, 01591 typename Kernel2D<T>::ConstAccessor, 01592 Diff2D, Diff2D, BorderTreatmentMode> 01593 kernel2d(Kernel2D<T> const & k) 01594 01595 { 01596 return 01597 tuple5<typename Kernel2D<T>::ConstIterator, 01598 typename Kernel2D<T>::ConstAccessor, 01599 Diff2D, Diff2D, BorderTreatmentMode>( 01600 k.center(), 01601 k.accessor(), 01602 k.upperLeft(), k.lowerRight(), 01603 k.borderTreatment()); 01604 } 01605 01606 template <class T> 01607 inline 01608 tuple5<typename Kernel2D<T>::ConstIterator, 01609 typename Kernel2D<T>::ConstAccessor, 01610 Diff2D, Diff2D, BorderTreatmentMode> 01611 kernel2d(Kernel2D<T> const & k, BorderTreatmentMode border) 01612 01613 { 01614 return 01615 tuple5<typename Kernel2D<T>::ConstIterator, 01616 typename Kernel2D<T>::ConstAccessor, 01617 Diff2D, Diff2D, BorderTreatmentMode>( 01618 k.center(), 01619 k.accessor(), 01620 k.upperLeft(), k.lowerRight(), 01621 border); 01622 } 01623 01624 01625 } // namespace vigra 01626 01627 #endif // VIGRA_STDCONVOLUTION_HXX
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|