[ 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.4.0, Dec 21 2005 ) */ 00008 /* The VIGRA Website is */ 00009 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00010 /* Please direct questions, bug reports, and contributions to */ 00011 /* koethe@informatik.uni-hamburg.de or */ 00012 /* vigra@kogs1.informatik.uni-hamburg.de */ 00013 /* */ 00014 /* Permission is hereby granted, free of charge, to any person */ 00015 /* obtaining a copy of this software and associated documentation */ 00016 /* files (the "Software"), to deal in the Software without */ 00017 /* restriction, including without limitation the rights to use, */ 00018 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00019 /* sell copies of the Software, and to permit persons to whom the */ 00020 /* Software is furnished to do so, subject to the following */ 00021 /* conditions: */ 00022 /* */ 00023 /* The above copyright notice and this permission notice shall be */ 00024 /* included in all copies or substantial portions of the */ 00025 /* Software. */ 00026 /* */ 00027 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00028 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00029 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00030 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00031 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00032 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00033 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00034 /* OTHER DEALINGS IN THE SOFTWARE. */ 00035 /* */ 00036 /************************************************************************/ 00037 00038 00039 #ifndef VIGRA_STDCONVOLUTION_HXX 00040 #define VIGRA_STDCONVOLUTION_HXX 00041 00042 #include <cmath> 00043 #include "vigra/stdimage.hxx" 00044 #include "vigra/bordertreatment.hxx" 00045 #include "vigra/separableconvolution.hxx" 00046 #include "vigra/utilities.hxx" 00047 #include "vigra/sized_int.hxx" 00048 00049 namespace vigra { 00050 00051 template <class SrcIterator, class SrcAccessor, 00052 class DestIterator, class DestAccessor, 00053 class KernelIterator, class KernelAccessor, 00054 class KSumType> 00055 void internalPixelEvaluationByClip(int x, int y, int w, int h, SrcIterator xs, 00056 SrcAccessor src_acc, DestIterator xd, DestAccessor dest_acc, 00057 KernelIterator ki, Diff2D kul, Diff2D klr, KernelAccessor ak, 00058 KSumType norm) 00059 { 00060 typedef typename 00061 NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType; 00062 typedef 00063 NumericTraits<typename DestAccessor::value_type> DestTraits; 00064 00065 // calculate width and height of the kernel 00066 int kernel_width = klr.x - kul.x + 1; 00067 int kernel_height = klr.y - kul.y + 1; 00068 00069 SumType sum = NumericTraits<SumType>::zero(); 00070 int xx, yy; 00071 int x0, y0, x1, y1; 00072 00073 y0 = (y<klr.y) ? -y : -klr.y; 00074 y1 = (h-y-1<-kul.y) ? h-y-1 : -kul.y; 00075 00076 x0 = (x<klr.x) ? -x : -klr.x; 00077 x1 = (w-x-1<-kul.x) ? w-x-1 : -kul.x; 00078 00079 SrcIterator yys = xs + Diff2D(x0, y0); 00080 KernelIterator yk = ki - Diff2D(x0, y0); 00081 00082 KSumType ksum = NumericTraits<KSumType>::zero(); 00083 kernel_width = x1 - x0 + 1; 00084 kernel_height = y1 - y0 + 1; 00085 00086 //es wird zuerst abgeschnitten und dann gespigelt! 00087 00088 for(yy=0; yy<kernel_height; ++yy, ++yys.y, --yk.y) 00089 { 00090 SrcIterator xxs = yys; 00091 KernelIterator xk = yk; 00092 00093 for(xx=0; xx<kernel_width; ++xx, ++xxs.x, --xk.x) 00094 { 00095 sum += ak(xk) * src_acc(xxs); 00096 ksum += ak(xk); 00097 } 00098 } 00099 00100 // store average in destination pixel 00101 dest_acc.set(DestTraits::fromRealPromote((norm / ksum) * sum), xd); 00102 00103 } 00104 00105 00106 #if 0 00107 00108 template <class SrcIterator, class SrcAccessor, 00109 class DestIterator, class DestAccessor, 00110 class KernelIterator, class KernelAccessor> 00111 void internalPixelEvaluationByWrapReflectRepeat(int x, int y, int src_width, int src_height, SrcIterator xs, 00112 SrcAccessor src_acc, DestIterator xd, DestAccessor dest_acc, 00113 KernelIterator ki, Diff2D kul, Diff2D klr, KernelAccessor ak, 00114 BorderTreatmentMode border) 00115 { 00116 00117 typedef typename 00118 NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType; 00119 typedef 00120 NumericTraits<typename DestAccessor::value_type> DestTraits; 00121 00122 SumType sum = NumericTraits<SumType>::zero(); 00123 00124 SrcIterator src_ul = xs - Diff2D(x, y); 00125 SrcIterator src_lr = src_ul + Diff2D(src_width, src_height); 00126 00127 SrcIterator yys = xs; 00128 KernelIterator yk = ki; 00129 00130 // calculate width and height of the kernel 00131 int kernel_width = klr.x - kul.x + 1; 00132 int kernel_height = klr.y - kul.y + 1; 00133 00134 // where the kernel is beyond the borders: 00135 bool top_to_much = (y<klr.y) ? true : false; 00136 bool down_to_much = (src_height-y-1<-kul.y)? true : false; 00137 bool left_to_much = (x<klr.x)? true : false; 00138 bool right_to_much = (src_width-x-1<-kul.x)? true : false; 00139 00140 // direction of iteration, 00141 // e.g. (-1, +1) for ll->ur or (-1, -1) for lr->ul 00142 Diff2D way_increment; 00143 00144 /* Iteration is always done from valid to invalid range. 00145 The following tuple is composed as such: 00146 - If an invalid range is reached while iterating in X, 00147 a jump of border_increment.first is performed and 00148 border_increment.third is used for further iterating. 00149 - If an invalid range is reached while iterating in Y, 00150 a jump of border_increment.second is performed and 00151 border_increment.fourth is used for further iterating. 00152 */ 00153 tuple4<int, int, int, int> border_increment; 00154 if (border == BORDER_TREATMENT_REPEAT){ 00155 border_increment = tuple4<int, int, int, int>(1, 1, 0, 0); 00156 }else if (border == BORDER_TREATMENT_REFLECT){ 00157 border_increment = tuple4<int, int, int, int>(2, 2, -1, -1); 00158 }else{ // BORDER_TREATMENT_WRAP 00159 border_increment = tuple4<int, int, int, int>(src_width, src_height, 1, 1); 00160 } 00161 00162 pair<int, int> valid_step_count; 00163 00164 if(left_to_much && !top_to_much && !down_to_much) 00165 { 00166 yys += klr; 00167 yk += kul; 00168 way_increment = Diff2D(-1, -1); 00169 border_increment.third = -border_increment.third; 00170 border_increment.fourth = -border_increment.fourth; 00171 valid_step_count = std::make_pair((yys - src_ul).x + 1, kernel_height); 00172 } 00173 else if(top_to_much && !left_to_much && !right_to_much) 00174 { 00175 yys += klr; 00176 yk += kul; 00177 way_increment = Diff2D(-1, -1); 00178 border_increment.third = -border_increment.third; 00179 border_increment.fourth = -border_increment.fourth; 00180 valid_step_count = std::make_pair(kernel_width, (yys - src_ul).y + 1); 00181 } 00182 else if(right_to_much && !top_to_much && !down_to_much) 00183 { 00184 yys += kul; 00185 yk += klr; 00186 way_increment = Diff2D(1, 1); 00187 border_increment.first = -border_increment.first; 00188 border_increment.second = -border_increment.second; 00189 valid_step_count = std::make_pair((src_lr - yys).x, kernel_height); 00190 } 00191 else if(down_to_much && !left_to_much && !right_to_much) 00192 { 00193 yys += kul; 00194 yk += klr; 00195 way_increment = Diff2D(1, 1); 00196 border_increment.first = -border_increment.first; 00197 border_increment.second = -border_increment.second; 00198 valid_step_count = std::make_pair(kernel_width, (src_lr - yys).y); 00199 } 00200 else if(down_to_much && left_to_much) 00201 { 00202 yys += kul + Diff2D(kernel_width - 1, 0); 00203 yk += kul + Diff2D(0, kernel_height - 1); 00204 way_increment = Diff2D(-1, 1); 00205 border_increment.second = -border_increment.second; 00206 border_increment.third = -border_increment.third; 00207 valid_step_count = std::make_pair((yys - src_ul).x + 1, (src_lr - yys).y); 00208 } 00209 else if(down_to_much && right_to_much) 00210 { 00211 yys += kul; 00212 yk += klr; 00213 way_increment = Diff2D(1, 1); 00214 border_increment.first = -border_increment.first; 00215 border_increment.second = -border_increment.second; 00216 valid_step_count = std::make_pair((src_lr - yys).x, (src_lr - yys).y); 00217 } 00218 else if(top_to_much && left_to_much) 00219 { 00220 yys += klr; 00221 yk += kul; 00222 way_increment = Diff2D(-1, -1); 00223 border_increment.third = -border_increment.third; 00224 border_increment.fourth = -border_increment.fourth; 00225 valid_step_count = std::make_pair((yys - src_ul).x + 1, (yys - src_ul).y + 1); 00226 } 00227 else 00228 { //top_to_much && right_to_much 00229 yys += kul + Diff2D(0, kernel_height - 1); 00230 yk += kul + Diff2D(kernel_width - 1, 0); 00231 way_increment = Diff2D(1, -1); 00232 border_increment.first = -border_increment.first; 00233 border_increment.fourth = -border_increment.fourth; 00234 valid_step_count = std::make_pair((src_lr - yys).x, (yys - src_ul).y + 1); 00235 } 00236 00237 int yy = 0, xx; 00238 00239 //laeuft den zulässigen Bereich in y-Richtung durch 00240 for(; yy < valid_step_count.second; ++yy, yys.y += way_increment.y, yk.y -= way_increment.y ) 00241 { 00242 SrcIterator xxs = yys; 00243 KernelIterator xk = yk; 00244 00245 //laeuft den zulässigen Bereich in x-Richtung durch 00246 for(xx = 0; xx < valid_step_count.first; ++xx, xxs.x += way_increment.x, xk.x -= way_increment.x) 00247 { 00248 sum += ak(xk) * src_acc(xxs); 00249 } 00250 00251 //Nächstes ++xxs.x wuerde in unzulässigen Bereich 00252 //bringen => Sprung in zulaessigen Bereich 00253 xxs.x += border_increment.first; 00254 00255 for( ; xx < kernel_width; ++xx, xxs.x += border_increment.third, xk.x -= way_increment.x ) 00256 { 00257 sum += ak(xk) * src_acc(xxs); 00258 } 00259 } 00260 00261 //Nächstes ++yys.y wuerde in unzulässigen Bereich 00262 //bringen => Sprung in zulaessigen Bereich 00263 yys.y += border_increment.second; 00264 00265 for( ; yy < kernel_height; ++yy, yys.y += border_increment.third, yk.y -= way_increment.y) 00266 { 00267 SrcIterator xxs = yys; 00268 KernelIterator xk = yk; 00269 00270 for(xx=0; xx < valid_step_count.first; ++xx, xxs.x += way_increment.x, xk.x -= way_increment.x) 00271 { 00272 sum += ak(xk) * src_acc(xxs); 00273 } 00274 00275 //Sprung in den zulaessigen Bereich 00276 xxs.x += border_increment.first; 00277 00278 for( ; xx < kernel_width; ++xx, xxs.x += border_increment.third, xk.x -= way_increment.x ) 00279 { 00280 sum += ak(xk) * src_acc(xxs); 00281 } 00282 } 00283 00284 // store average in destination pixel 00285 dest_acc.set(DestTraits::fromRealPromote(sum), xd); 00286 00287 }// end of internalPixelEvaluationByWrapReflectRepeat 00288 #endif /* #if 0 */ 00289 00290 00291 template <class SrcIterator, class SrcAccessor, 00292 class KernelIterator, class KernelAccessor, 00293 class SumType> 00294 void 00295 internalPixelEvaluationByWrapReflectRepeat(SrcIterator xs, SrcAccessor src_acc, 00296 KernelIterator xk, KernelAccessor ak, 00297 int left, int right, int kleft, int kright, 00298 int borderskipx, int borderinc, SumType & sum) 00299 { 00300 SrcIterator xxs = xs + left; 00301 KernelIterator xxk = xk - left; 00302 00303 for(int xx = left; xx <= right; ++xx, ++xxs, --xxk) 00304 { 00305 sum += ak(xxk) * src_acc(xxs); 00306 } 00307 00308 xxs = xs + left - borderskipx; 00309 xxk = xk - left + 1; 00310 for(int xx = left - 1; xx >= -kright; --xx, xxs -= borderinc, ++xxk) 00311 { 00312 sum += ak(xxk) * src_acc(xxs); 00313 } 00314 00315 xxs = xs + right + borderskipx; 00316 xxk = xk - right - 1; 00317 for(int xx = right + 1; xx <= -kleft; ++xx, xxs += borderinc, --xxk) 00318 { 00319 sum += ak(xxk) * src_acc(xxs); 00320 } 00321 } 00322 00323 00324 /** \addtogroup StandardConvolution Two-dimensional convolution functions 00325 00326 Perform 2D non-separable convolution, with and without ROI mask. 00327 00328 These generic convolution functions implement 00329 the standard 2D convolution operation for images that fit 00330 into the required interface. Arbitrary ROI's are supported 00331 by the mask version of the algorithm. 00332 The functions need a suitable 2D kernel to operate. 00333 */ 00334 //@{ 00335 00336 /** \brief Performs a 2 dimensional convolution of the source image using the given 00337 kernel. 00338 00339 The KernelIterator must point to the center of the kernel, and 00340 the kernel's size is given by its upper left (x and y of distance <= 0) and 00341 lower right (distance >= 0) corners. The image must always be larger than the 00342 kernel. At those positions where the kernel does not completely fit 00343 into the image, the specified \ref BorderTreatmentMode is 00344 applied. You can choice between following BorderTreatmentModes: 00345 <ul> 00346 <li>BORDER_TREATMENT_CLIP</li> 00347 <li>BORDER_TREATMENT_AVOID</li> 00348 <li>BORDER_TREATMENT_WRAP</li> 00349 <li>BORDER_TREATMENT_REFLECT</li> 00350 <li>BORDER_TREATMENT_REPEAT</li> 00351 </ul><br> 00352 The images's pixel type (SrcAccessor::value_type) must be a 00353 linear space over the kernel's value_type (KernelAccessor::value_type), 00354 i.e. addition of source values, multiplication with kernel values, 00355 and NumericTraits must be defined. 00356 The kernel's value_type must be an algebraic field, 00357 i.e. the arithmetic operations (+, -, *, /) and NumericTraits must 00358 be defined. 00359 00360 <b> Declarations:</b> 00361 00362 pass arguments explicitly: 00363 \code 00364 namespace vigra { 00365 template <class SrcIterator, class SrcAccessor, 00366 class DestIterator, class DestAccessor, 00367 class KernelIterator, class KernelAccessor> 00368 void convolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc, 00369 DestIterator dest_ul, DestAccessor dest_acc, 00370 KernelIterator ki, KernelAccessor ak, 00371 Diff2D kul, Diff2D klr, BorderTreatmentMode border); 00372 } 00373 \endcode 00374 00375 00376 use argument objects in conjunction with \ref ArgumentObjectFactories: 00377 \code 00378 namespace vigra { 00379 template <class SrcIterator, class SrcAccessor, 00380 class DestIterator, class DestAccessor, 00381 class KernelIterator, class KernelAccessor> 00382 void convolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00383 pair<DestIterator, DestAccessor> dest, 00384 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, 00385 BorderTreatmentMode> kernel); 00386 } 00387 \endcode 00388 00389 <b> Usage:</b> 00390 00391 <b>\#include</b> "<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>"<br> 00392 Namespace: vigra 00393 00394 00395 \code 00396 vigra::FImage src(w,h), dest(w,h); 00397 ... 00398 00399 // define horizontal Sobel filter 00400 vigra::Kernel2D<float> sobel; 00401 00402 sobel.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = // upper left and lower right 00403 0.125, 0.0, -0.125, 00404 0.25, 0.0, -0.25, 00405 0.125, 0.0, -0.125; 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 int x,y; 00502 00503 KernelSumType norm = NumericTraits<KernelSumType>::zero(); 00504 if(border == BORDER_TREATMENT_CLIP) 00505 { 00506 // caluclate the sum of the kernel elements for renormalization 00507 KernelIterator yk = ki + klr; 00508 00509 //Die Summe der Punkte im Kernel wird ermittelt (= norm) 00510 for(y=0; y<kernel_height; ++y, --yk.y) 00511 { 00512 KernelIterator xk = yk; 00513 for(x=0; x<kernel_width; ++x, --xk.x) 00514 { 00515 norm += ak(xk); 00516 } 00517 } 00518 vigra_precondition(norm != NumericTraits<KernelSumType>::zero(), 00519 "convolveImage(): Cannot use BORDER_TREATMENT_CLIP with a DC-free kernel"); 00520 } 00521 00522 // create iterators for the interior part of the image (where the kernel always fits into the image) 00523 DestIterator yd = dest_ul + Diff2D(klr.x, klr.y); 00524 SrcIterator ys = src_ul + Diff2D(klr.x, klr.y); 00525 SrcIterator send = src_lr + Diff2D(kul.x, kul.y); 00526 00527 // iterate over the interior part 00528 for(; ys.y < send.y; ++ys.y, ++yd.y) 00529 { 00530 // create x iterators 00531 DestIterator xd(yd); 00532 SrcIterator xs(ys); 00533 00534 for(; xs.x < send.x; ++x, ++xs.x, ++xd.x) 00535 { 00536 // init the sum 00537 SumType sum = NumericTraits<SumType>::zero(); 00538 00539 SrcIterator yys = xs - klr; 00540 SrcIterator yyend = xs - kul; 00541 KernelIterator yk = ki + klr; 00542 00543 for(; yys.y <= yyend.y; ++yys.y, --yk.y) 00544 { 00545 typename SrcIterator::row_iterator xxs = yys.rowIterator(); 00546 typename SrcIterator::row_iterator xxe = xxs + kernel_width; 00547 typename KernelIterator::row_iterator xk = yk.rowIterator(); 00548 00549 for(; xxs < xxe; ++xxs, --xk) 00550 { 00551 sum += ak(xk) * src_acc(xxs); 00552 } 00553 } 00554 00555 // store convolution result in destination pixel 00556 dest_acc.set(DestTraits::fromRealPromote(sum), xd); 00557 } 00558 } 00559 00560 if(border == BORDER_TREATMENT_AVOID) 00561 return; // skip processing near the border 00562 00563 int interiorskip = w + kul.x - klr.x - 1; 00564 int borderskipx; 00565 int borderskipy; 00566 int borderinc; 00567 if(border == BORDER_TREATMENT_REPEAT) 00568 { 00569 borderskipx = 0; 00570 borderskipy = 0; 00571 borderinc = 0; 00572 } 00573 else if(border == BORDER_TREATMENT_REFLECT) 00574 { 00575 borderskipx = -1; 00576 borderskipy = -1; 00577 borderinc = -1; 00578 } 00579 else if(border == BORDER_TREATMENT_WRAP) 00580 { 00581 borderskipx = -w+1; 00582 borderskipy = -h+1; 00583 borderinc = 1; 00584 } 00585 00586 // create iterators for the entire image 00587 yd = dest_ul; 00588 ys = src_ul; 00589 00590 // go over the entire image (but skip the already computed points in the loop) 00591 for(y=0; y < h; ++y, ++ys.y, ++yd.y) 00592 { 00593 int top = std::max(static_cast<IntBiggest>(-klr.y), 00594 static_cast<IntBiggest>(src_ul.y - ys.y)); 00595 int bottom = std::min(static_cast<IntBiggest>(-kul.y), 00596 static_cast<IntBiggest>(src_lr.y - ys.y - 1)); 00597 00598 // create x iterators 00599 DestIterator xd(yd); 00600 SrcIterator xs(ys); 00601 00602 for(x=0; x < w; ++x, ++xs.x, ++xd.x) 00603 { 00604 // check if we are away from the border 00605 if(y >= klr.y && y < h+kul.y && x == klr.x) 00606 { 00607 // yes => skip the already computed points 00608 x += interiorskip; 00609 xs.x += interiorskip; 00610 xd.x += interiorskip; 00611 continue; 00612 } 00613 if (border == BORDER_TREATMENT_CLIP) 00614 { 00615 internalPixelEvaluationByClip(x, y, w, h, xs, src_acc, xd, dest_acc, ki, kul, klr, ak, norm); 00616 } 00617 else 00618 { 00619 int left = std::max(-klr.x, src_ul.x - xs.x); 00620 int right = std::min(-kul.x, src_lr.x - xs.x - 1); 00621 00622 // init the sum 00623 SumType sum = NumericTraits<SumType>::zero(); 00624 00625 // create iterators for the part of the kernel that fits into the image 00626 SrcIterator yys = xs + Size2D(0, top); 00627 KernelIterator yk = ki - Size2D(0, top); 00628 00629 int yy; 00630 for(yy = top; yy <= bottom; ++yy, ++yys.y, --yk.y) 00631 { 00632 internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak, 00633 left, right, kul.x, klr.x, borderskipx, borderinc, sum); 00634 } 00635 yys = xs + Size2D(0, top - borderskipy); 00636 yk = ki - Size2D(0, top - 1); 00637 for(yy = top - 1; yy >= -klr.y; --yy, yys.y -= borderinc, ++yk.y) 00638 { 00639 internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak, 00640 left, right, kul.x, klr.x, borderskipx, borderinc, sum); 00641 } 00642 yys = xs + Size2D(0, bottom + borderskipy); 00643 yk = ki - Size2D(0, bottom + 1); 00644 for(yy = bottom + 1; yy <= -kul.y; ++yy, yys.y += borderinc, --yk.y) 00645 { 00646 internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak, 00647 left, right, kul.x, klr.x, borderskipx, borderinc, sum); 00648 } 00649 00650 // store convolution result in destination pixel 00651 dest_acc.set(DestTraits::fromRealPromote(sum), xd); 00652 00653 // internalPixelEvaluationByWrapReflectRepeat(x, y, w, h, xs, src_acc, xd, dest_acc, ki, kul, klr, ak, border); 00654 } 00655 } 00656 } 00657 } 00658 00659 00660 template <class SrcIterator, class SrcAccessor, 00661 class DestIterator, class DestAccessor, 00662 class KernelIterator, class KernelAccessor> 00663 inline 00664 void convolveImage( 00665 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00666 pair<DestIterator, DestAccessor> dest, 00667 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, 00668 BorderTreatmentMode> kernel) 00669 { 00670 convolveImage(src.first, src.second, src.third, 00671 dest.first, dest.second, 00672 kernel.first, kernel.second, kernel.third, 00673 kernel.fourth, kernel.fifth); 00674 } 00675 00676 00677 /** \brief Performs a 2-dimensional normalized convolution, i.e. convolution with a mask image. 00678 00679 This functions computes 00680 <a href ="http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/PIRODDI1/NormConv/NormConv.html">normalized 00681 convolution</a> as defined in 00682 Knutsson, H. and Westin, C-F.: <i>Normalized and differential convolution: 00683 Methods for Interpolation and Filtering of incomplete and uncertain data</i>. 00684 Proc. of the IEEE Conf. on Computer Vision and Pattern Recognition, 1993, 515-523. 00685 00686 The mask image must be binary and encodes which pixels of the original image 00687 are valid. It is used as follows: 00688 Only pixel under the mask are used in the calculations. Whenever a part of the 00689 kernel lies outside the mask, it is ignored, and the kernel is renormalized to its 00690 original norm (analogous to the CLIP \ref BorderTreatmentMode). Thus, a useful convolution 00691 result is computed whenever <i>at least one valid pixel is within the current window</i> 00692 Thus, destination pixels not under the mask still receive a value if they are <i>near</i> 00693 the mask. Therefore, this algorithm is useful as an interpolator of sparse input data. 00694 If you are only interested in the destination values under the mask, you can perform 00695 a subsequent \ref copyImageIf(). 00696 00697 The KernelIterator must point to the center of the kernel, and 00698 the kernel's size is given by its upper left (x and y of distance <= 0) and 00699 lower right (distance >= 0) corners. The image must always be larger than the 00700 kernel. At those positions where the kernel does not completely fit 00701 into the image, the specified \ref BorderTreatmentMode is 00702 applied. Only BORDER_TREATMENT_CLIP and BORDER_TREATMENT_AVOID are currently 00703 supported. 00704 00705 The images's pixel type (SrcAccessor::value_type) must be a 00706 linear space over the kernel's value_type (KernelAccessor::value_type), 00707 i.e. addition of source values, multiplication with kernel values, 00708 and NumericTraits must be defined. 00709 The kernel's value_type must be an algebraic field, 00710 i.e. the arithmetic operations (+, -, *, /) and NumericTraits must 00711 be defined. 00712 00713 <b> Declarations:</b> 00714 00715 pass arguments explicitly: 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 void 00723 normalizedConvolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc, 00724 MaskIterator mul, MaskAccessor am, 00725 DestIterator dest_ul, DestAccessor dest_acc, 00726 KernelIterator ki, KernelAccessor ak, 00727 Diff2D kul, Diff2D klr, BorderTreatmentMode border); 00728 } 00729 \endcode 00730 00731 00732 use argument objects in conjunction with \ref ArgumentObjectFactories: 00733 \code 00734 namespace vigra { 00735 template <class SrcIterator, class SrcAccessor, 00736 class MaskIterator, class MaskAccessor, 00737 class DestIterator, class DestAccessor, 00738 class KernelIterator, class KernelAccessor> 00739 inline 00740 void normalizedConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00741 pair<MaskIterator, MaskAccessor> mask, 00742 pair<DestIterator, DestAccessor> dest, 00743 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, 00744 BorderTreatmentMode> kernel); 00745 } 00746 \endcode 00747 00748 <b> Usage:</b> 00749 00750 <b>\#include</b> "<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>"<br> 00751 Namespace: vigra 00752 00753 00754 \code 00755 vigra::FImage src(w,h), dest(w,h); 00756 vigra::CImage mask(w,h); 00757 ... 00758 00759 // define 3x3 binomial filter 00760 vigra::Kernel2D<float> binom; 00761 00762 binom.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = // upper left and lower right 00763 0.0625, 0.125, 0.0625, 00764 0.125, 0.25, 0.125, 00765 0.0625, 0.125, 0.0625; 00766 00767 vigra::normalizedConvolveImage(srcImageRange(src), maskImage(mask), destImage(dest), kernel2d(binom)); 00768 \endcode 00769 00770 <b> Required Interface:</b> 00771 00772 \code 00773 ImageIterator src_ul, src_lr; 00774 ImageIterator mul; 00775 ImageIterator dest_ul; 00776 ImageIterator ik; 00777 00778 SrcAccessor src_accessor; 00779 MaskAccessor mask_accessor; 00780 DestAccessor dest_accessor; 00781 KernelAccessor kernel_accessor; 00782 00783 NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(src_ul); 00784 00785 s = s + s; 00786 s = kernel_accessor(ik) * s; 00787 s -= s; 00788 00789 if(mask_accessor(mul)) ...; 00790 00791 dest_accessor.set( 00792 NumericTraits<DestAccessor::value_type>::fromRealPromote(s), dest_ul); 00793 00794 NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik); 00795 00796 k += k; 00797 k -= k; 00798 k = k / k; 00799 00800 \endcode 00801 00802 <b> Preconditions:</b> 00803 00804 \code 00805 kul.x <= 0 00806 kul.y <= 0 00807 klr.x >= 0 00808 klr.y >= 0 00809 src_lr.x - src_ul.x >= klr.x + kul.x + 1 00810 src_lr.y - src_ul.y >= klr.y + kul.y + 1 00811 border == BORDER_TREATMENT_CLIP || border == BORDER_TREATMENT_AVOID 00812 \endcode 00813 00814 Sum of kernel elements must be != 0. 00815 00816 */ 00817 template <class SrcIterator, class SrcAccessor, 00818 class DestIterator, class DestAccessor, 00819 class MaskIterator, class MaskAccessor, 00820 class KernelIterator, class KernelAccessor> 00821 void 00822 normalizedConvolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc, 00823 MaskIterator mul, MaskAccessor am, 00824 DestIterator dest_ul, DestAccessor dest_acc, 00825 KernelIterator ki, KernelAccessor ak, 00826 Diff2D kul, Diff2D klr, BorderTreatmentMode border) 00827 { 00828 vigra_precondition((border == BORDER_TREATMENT_CLIP || 00829 border == BORDER_TREATMENT_AVOID), 00830 "normalizedConvolveImage(): " 00831 "Border treatment must be BORDER_TREATMENT_CLIP or BORDER_TREATMENT_AVOID."); 00832 00833 vigra_precondition(kul.x <= 0 && kul.y <= 0, 00834 "normalizedConvolveImage(): left borders must be <= 0."); 00835 vigra_precondition(klr.x >= 0 && klr.y >= 0, 00836 "normalizedConvolveImage(): right borders must be >= 0."); 00837 00838 // use traits to determine SumType as to prevent possible overflow 00839 typedef typename 00840 NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType; 00841 typedef typename 00842 NumericTraits<typename KernelAccessor::value_type>::RealPromote KSumType; 00843 typedef 00844 NumericTraits<typename DestAccessor::value_type> DestTraits; 00845 00846 // calculate width and height of the image 00847 int w = src_lr.x - src_ul.x; 00848 int h = src_lr.y - src_ul.y; 00849 int kernel_width = klr.x - kul.x + 1; 00850 int kernel_height = klr.y - kul.y + 1; 00851 00852 int x,y; 00853 int ystart = (border == BORDER_TREATMENT_AVOID) ? klr.y : 0; 00854 int yend = (border == BORDER_TREATMENT_AVOID) ? h+kul.y : h; 00855 int xstart = (border == BORDER_TREATMENT_AVOID) ? klr.x : 0; 00856 int xend = (border == BORDER_TREATMENT_AVOID) ? w+kul.x : w; 00857 00858 // create y iterators 00859 DestIterator yd = dest_ul + Diff2D(xstart, ystart); 00860 SrcIterator ys = src_ul + Diff2D(xstart, ystart); 00861 MaskIterator ym = mul + Diff2D(xstart, ystart); 00862 00863 KSumType norm = ak(ki); 00864 int xx, yy; 00865 KernelIterator yk = ki + klr; 00866 for(yy=0; yy<kernel_height; ++yy, --yk.y) 00867 { 00868 KernelIterator xk = yk; 00869 00870 for(xx=0; xx<kernel_width; ++xx, --xk.x) 00871 { 00872 norm += ak(xk); 00873 } 00874 } 00875 norm -= ak(ki); 00876 00877 00878 for(y=ystart; y < yend; ++y, ++ys.y, ++yd.y, ++ym.y) 00879 { 00880 // create x iterators 00881 DestIterator xd(yd); 00882 SrcIterator xs(ys); 00883 MaskIterator xm(ym); 00884 00885 for(x=xstart; x < xend; ++x, ++xs.x, ++xd.x, ++xm.x) 00886 { 00887 // how much of the kernel fits into the image ? 00888 int x0, y0, x1, y1; 00889 00890 y0 = (y<klr.y) ? -y : -klr.y; 00891 y1 = (h-y-1<-kul.y) ? h-y-1 : -kul.y; 00892 x0 = (x<klr.x) ? -x : -klr.x; 00893 x1 = (w-x-1<-kul.x) ? w-x-1 : -kul.x; 00894 00895 bool first = true; 00896 // init the sum 00897 SumType sum; 00898 KSumType ksum; 00899 00900 SrcIterator yys = xs + Diff2D(x0, y0); 00901 MaskIterator yym = xm + Diff2D(x0, y0); 00902 KernelIterator yk = ki - Diff2D(x0, y0); 00903 00904 int xx, kernel_width, kernel_height; 00905 kernel_width = x1 - x0 + 1; 00906 kernel_height = y1 - y0 + 1; 00907 for(yy=0; yy<kernel_height; ++yy, ++yys.y, --yk.y, ++yym.y) 00908 { 00909 typename SrcIterator::row_iterator xxs = yys.rowIterator(); 00910 typename SrcIterator::row_iterator xxend = xxs + kernel_width; 00911 typename MaskIterator::row_iterator xxm = yym.rowIterator(); 00912 typename KernelIterator::row_iterator xk = yk.rowIterator(); 00913 00914 for(xx=0; xxs < xxend; ++xxs.x, --xk.x, ++xxm.x) 00915 { 00916 if(!am(xxm)) continue; 00917 00918 if(first) 00919 { 00920 sum = ak(xk) * src_acc(xxs); 00921 ksum = ak(xk); 00922 first = false; 00923 } 00924 else 00925 { 00926 sum += ak(xk) * src_acc(xxs); 00927 ksum += ak(xk); 00928 } 00929 } 00930 } 00931 // store average in destination pixel 00932 if(!first && 00933 ksum != NumericTraits<KSumType>::zero()) 00934 { 00935 dest_acc.set(DestTraits::fromRealPromote((norm / ksum) * sum), xd); 00936 } 00937 } 00938 } 00939 } 00940 00941 00942 template <class SrcIterator, class SrcAccessor, 00943 class DestIterator, class DestAccessor, 00944 class MaskIterator, class MaskAccessor, 00945 class KernelIterator, class KernelAccessor> 00946 inline 00947 void normalizedConvolveImage( 00948 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00949 pair<MaskIterator, MaskAccessor> mask, 00950 pair<DestIterator, DestAccessor> dest, 00951 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, 00952 BorderTreatmentMode> kernel) 00953 { 00954 normalizedConvolveImage(src.first, src.second, src.third, 00955 mask.first, mask.second, 00956 dest.first, dest.second, 00957 kernel.first, kernel.second, kernel.third, 00958 kernel.fourth, kernel.fifth); 00959 } 00960 00961 /** \brief Deprecated name of 2-dimensional normalized convolution, i.e. convolution with a mask image. 00962 00963 See \ref normalizedConvolveImage() for documentation. 00964 00965 <b> Declarations:</b> 00966 00967 pass arguments explicitly: 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 void 00975 convolveImageWithMask(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc, 00976 MaskIterator mul, MaskAccessor am, 00977 DestIterator dest_ul, DestAccessor dest_acc, 00978 KernelIterator ki, KernelAccessor ak, 00979 Diff2D kul, Diff2D klr, BorderTreatmentMode border); 00980 } 00981 \endcode 00982 00983 00984 use argument objects in conjunction with \ref ArgumentObjectFactories: 00985 \code 00986 namespace vigra { 00987 template <class SrcIterator, class SrcAccessor, 00988 class MaskIterator, class MaskAccessor, 00989 class DestIterator, class DestAccessor, 00990 class KernelIterator, class KernelAccessor> 00991 inline 00992 void convolveImageWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00993 pair<MaskIterator, MaskAccessor> mask, 00994 pair<DestIterator, DestAccessor> dest, 00995 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, 00996 BorderTreatmentMode> kernel); 00997 } 00998 \endcode 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 (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|