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