Main Page | Modules | Namespace List | Class Hierarchy | Class List | File List | Namespace Members | Class Members | File Members | Related Pages
hugin_base/vigra_ext/pyramid2.h
Go to the documentation of this file.00001 /* 00002 * Copyright (C) 2004-2007 Andrew Mihal 00003 * 00004 * This file is part of Enblend. 00005 * 00006 * Stripped down and adjust for use in hugin by Pablo d'Angelo 00007 * 00008 * Enblend is free software; you can redistribute it and/or modify 00009 * it under the terms of the GNU General Public License as published by 00010 * the Free Software Foundation; either version 2 of the License, or 00011 * (at your option) any later version. 00012 * 00013 * Enblend is distributed in the hope that it will be useful, 00014 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00015 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00016 * GNU General Public License for more details. 00017 * 00018 * You should have received a copy of the GNU General Public License 00019 * along with Enblend; if not, write to the Free Software 00020 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00021 */ 00022 #ifndef __PYRAMID_H__ 00023 #define __PYRAMID_H__ 00024 00025 #ifdef HAVE_CONFIG_H 00026 #include <config.h> 00027 #endif 00028 00029 #include <functional> 00030 #include <vector> 00031 00032 #include <vigra/convolution.hxx> 00033 #include <vigra/error.hxx> 00034 #include <vigra/inspectimage.hxx> 00035 #include <vigra/numerictraits.hxx> 00036 #include <vigra/rgbvalue.hxx> 00037 #include <vigra/sized_int.hxx> 00038 #include <vigra/transformimage.hxx> 00039 00040 //#include "fixmath.h" 00041 00042 namespace enblend { 00043 00044 using std::cout; 00045 using std::vector; 00046 using std::pair; 00047 00048 using vigra::linearRangeMapping; 00049 using vigra::NumericTraits; 00050 using vigra::transformImage; 00051 using vigra::triple; 00052 using vigra::UInt16; 00053 using vigra::UInt16Image; 00054 using vigra::UInt16RGBImage; 00055 using vigra::Diff2D; 00056 00057 00058 #define IMUL6(A) (A * SKIPSMImagePixelType(6)) 00059 #define IMUL5(A) (A * SKIPSMImagePixelType(5)) 00060 #define IMUL11(A) (A * SKIPSMImagePixelType(11)) 00061 #define AMUL6(A) (A * SKIPSMAlphaPixelType(6)) 00062 00068 template <typename ImagePixelComponentType> 00069 unsigned int filterHalfWidth(const unsigned int levels) { 00070 // For levels >= 30, the full width will just barely fit in int32. 00071 // When this range is added to a bounding box it will certainly 00072 // overflow the Diff2D. 00073 vigra_precondition((levels >= 1 && levels <= 29), 00074 "filterHalfWidth: levels outside of range [1,29]"); 00075 00076 // This is the arithmetic half width. 00077 int halfWidth = (levels == 1) ? 0 : ((1 << (levels+1)) - 4); 00078 00079 return halfWidth; 00080 } 00081 00173 template <typename SKIPSMImagePixelType, typename SKIPSMAlphaPixelType, 00174 typename SrcImageIterator, typename SrcAccessor, 00175 typename AlphaIterator, typename AlphaAccessor, 00176 typename DestImageIterator, typename DestAccessor, 00177 typename DestAlphaIterator, typename DestAlphaAccessor> 00178 inline void reduce(bool wraparound, 00179 SrcImageIterator src_upperleft, 00180 SrcImageIterator src_lowerright, 00181 SrcAccessor sa, 00182 AlphaIterator alpha_upperleft, 00183 AlphaAccessor aa, 00184 DestImageIterator dest_upperleft, 00185 DestImageIterator dest_lowerright, 00186 DestAccessor da, 00187 DestAlphaIterator dest_alpha_upperleft, 00188 DestAlphaIterator dest_alpha_lowerright, 00189 DestAlphaAccessor daa) { 00190 00191 typedef typename DestAccessor::value_type DestPixelType; 00192 typedef typename DestAlphaAccessor::value_type DestAlphaPixelType; 00193 00194 int src_w = src_lowerright.x - src_upperleft.x; 00195 int src_h = src_lowerright.y - src_upperleft.y; 00196 int dst_w = dest_lowerright.x - dest_upperleft.x; 00197 //int dst_h = dest_lowerright.y - dest_upperleft.y; 00198 00199 vigra_precondition(src_w > 1 && src_h > 1, 00200 "src image too small in reduce"); 00201 00202 // State variables for source image pixel values 00203 SKIPSMImagePixelType isr0, isr1, isrp; 00204 SKIPSMImagePixelType *isc0 = new SKIPSMImagePixelType[dst_w + 1]; 00205 SKIPSMImagePixelType *isc1 = new SKIPSMImagePixelType[dst_w + 1]; 00206 SKIPSMImagePixelType *iscp = new SKIPSMImagePixelType[dst_w + 1]; 00207 00208 // State variables for source mask pixel values 00209 SKIPSMAlphaPixelType asr0, asr1, asrp; 00210 SKIPSMAlphaPixelType *asc0 = new SKIPSMAlphaPixelType[dst_w + 1]; 00211 SKIPSMAlphaPixelType *asc1 = new SKIPSMAlphaPixelType[dst_w + 1]; 00212 SKIPSMAlphaPixelType *ascp = new SKIPSMAlphaPixelType[dst_w + 1]; 00213 00214 // Convenient constants 00215 const SKIPSMImagePixelType SKIPSMImageZero(NumericTraits<SKIPSMImagePixelType>::zero()); 00216 const SKIPSMAlphaPixelType SKIPSMAlphaZero(NumericTraits<SKIPSMAlphaPixelType>::zero()); 00217 const SKIPSMAlphaPixelType SKIPSMAlphaOne(NumericTraits<SKIPSMAlphaPixelType>::one()); 00218 const DestPixelType DestImageZero(NumericTraits<DestPixelType>::zero()); 00219 const DestAlphaPixelType DestAlphaZero(NumericTraits<DestAlphaPixelType>::zero()); 00220 const DestAlphaPixelType DestAlphaMax(NumericTraits<DestAlphaPixelType>::max()); 00221 00222 DestImageIterator dy = dest_upperleft; 00223 DestImageIterator dx = dy; 00224 SrcImageIterator sy = src_upperleft; 00225 SrcImageIterator sx = sy; 00226 AlphaIterator ay = alpha_upperleft; 00227 AlphaIterator ax = ay; 00228 DestAlphaIterator day = dest_alpha_upperleft; 00229 DestAlphaIterator dax = day; 00230 00231 bool evenY = true; 00232 bool evenX = true; 00233 int srcy = 0; 00234 int srcx = 0; 00235 //int dsty = 0; 00236 int dstx = 0; 00237 00238 // First row 00239 { 00240 if (wraparound) { 00241 asr0 = aa(ay, Diff2D(src_w-2, 0)) ? SKIPSMAlphaOne : SKIPSMAlphaZero; 00242 asr1 = SKIPSMAlphaZero; 00243 asrp = aa(ay, Diff2D(src_w-1, 0)) ? (SKIPSMAlphaOne * 4) : SKIPSMAlphaZero; 00244 isr0 = aa(ay, Diff2D(src_w-2, 0)) ? SKIPSMImagePixelType(sa(sy, Diff2D(src_w-2, 0))) : SKIPSMImageZero; 00245 isr1 = SKIPSMImageZero; 00246 isrp = aa(ay, Diff2D(src_w-1, 0)) ? (SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1, 0))) * 4) : SKIPSMImageZero; 00247 } else { 00248 asr0 = SKIPSMAlphaZero; 00249 asr1 = SKIPSMAlphaZero; 00250 asrp = SKIPSMAlphaZero; 00251 isr0 = SKIPSMImageZero; 00252 isr1 = SKIPSMImageZero; 00253 isrp = SKIPSMImageZero; 00254 } 00255 00256 for (sx = sy, ax = ay, evenX = true, srcx = 0, dstx = 0; srcx < src_w; ++srcx, ++sx.x, ++ax.x) { 00257 SKIPSMAlphaPixelType mcurrent(aa(ax) ? SKIPSMAlphaOne : SKIPSMAlphaZero); 00258 SKIPSMImagePixelType icurrent(aa(ax) ? SKIPSMImagePixelType(sa(sx)) : SKIPSMImageZero); 00259 if (evenX) { 00260 asc1[dstx] = SKIPSMAlphaZero; 00261 asc0[dstx] = asr1 + AMUL6(asr0) + asrp + mcurrent; 00262 asr1 = asr0 + asrp; 00263 asr0 = mcurrent; 00264 isc1[dstx] = SKIPSMImageZero; 00265 isc0[dstx] = isr1 + IMUL6(isr0) + isrp + icurrent; 00266 isr1 = isr0 + isrp; 00267 isr0 = icurrent; 00268 } 00269 else { 00270 asrp = mcurrent * 4; 00271 isrp = icurrent * 4; 00272 ++dstx; 00273 } 00274 evenX = !evenX; 00275 } 00276 00277 // Last entries in first row 00278 if (!evenX) { 00279 // previous srcx was even 00280 ++dstx; 00281 if (wraparound) { 00282 asc1[dstx] = SKIPSMAlphaZero; 00283 asc0[dstx] = asr1 + AMUL6(asr0) + (aa(ay) ? (SKIPSMAlphaOne * 4) : SKIPSMAlphaZero) + (aa(ay, Diff2D(1,0)) ? SKIPSMAlphaOne : SKIPSMAlphaZero); 00284 isc1[dstx] = SKIPSMImageZero; 00285 isc0[dstx] = isr1 + IMUL6(isr0) + (aa(ay) ? (SKIPSMImagePixelType(sa(sy)) * 4) : SKIPSMImageZero) 00286 + (aa(ay, Diff2D(1,0)) ? SKIPSMImagePixelType(sa(sy, Diff2D(1,0))) : SKIPSMImageZero); 00287 } else { 00288 asc1[dstx] = SKIPSMAlphaZero; 00289 asc0[dstx] = asr1 + AMUL6(asr0); 00290 isc1[dstx] = SKIPSMImageZero; 00291 isc0[dstx] = isr1 + IMUL6(isr0); 00292 } 00293 } 00294 else { 00295 // previous srcx was odd 00296 if (wraparound) { 00297 asc1[dstx] = SKIPSMAlphaZero; 00298 asc0[dstx] = asr1 + AMUL6(asr0) + asrp + (aa(ay) ? SKIPSMAlphaOne : SKIPSMAlphaZero); 00299 isc1[dstx] = SKIPSMImageZero; 00300 isc0[dstx] = isr1 + IMUL6(isr0) + isrp + (aa(ay) ? SKIPSMImagePixelType(sa(sy)) : SKIPSMImageZero); 00301 } else { 00302 asc1[dstx] = SKIPSMAlphaZero; 00303 asc0[dstx] = asr1 + AMUL6(asr0) + asrp; 00304 isc1[dstx] = SKIPSMImageZero; 00305 isc0[dstx] = isr1 + IMUL6(isr0) + isrp; 00306 } 00307 } 00308 } 00309 ++sy.y; 00310 ++ay.y; 00311 00312 // Main Rows 00313 { 00314 for (evenY = false, srcy = 1; srcy < src_h; ++srcy, ++sy.y, ++ay.y) { 00315 00316 if (wraparound) { 00317 asr0 = aa(ay, Diff2D(src_w-2, 0)) ? SKIPSMAlphaOne : SKIPSMAlphaZero; 00318 asr1 = SKIPSMAlphaZero; 00319 asrp = aa(ay, Diff2D(src_w-1, 0)) ? (SKIPSMAlphaOne * 4) : SKIPSMAlphaZero; 00320 isr0 = aa(ay, Diff2D(src_w-2, 0)) ? SKIPSMImagePixelType(sa(sy, Diff2D(src_w-2,0))) : SKIPSMImageZero; 00321 isr1 = SKIPSMImageZero; 00322 isrp = aa(ay, Diff2D(src_w-1, 0)) ? (SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1,0))) * 4) : SKIPSMImageZero; 00323 } else { 00324 asr0 = SKIPSMAlphaZero; 00325 asr1 = SKIPSMAlphaZero; 00326 asrp = SKIPSMAlphaZero; 00327 isr0 = SKIPSMImageZero; 00328 isr1 = SKIPSMImageZero; 00329 isrp = SKIPSMImageZero; 00330 } 00331 00332 if (evenY) { 00333 // Even-numbered row 00334 00335 // First entry in row 00336 sx = sy; 00337 ax = ay; 00338 if (wraparound) { 00339 asr1 = asr0 + asrp; 00340 isr1 = isr0 + isrp; 00341 } 00342 asr0 = aa(ax) ? SKIPSMAlphaOne : SKIPSMAlphaZero; 00343 isr0 = aa(ax) ? SKIPSMImagePixelType(sa(sx)) : SKIPSMImageZero; 00344 // isc*[0] are never used 00345 ++sx.x; 00346 ++ax.x; 00347 dx = dy; 00348 dax = day; 00349 00350 // Main entries in row 00351 for (evenX = false, srcx = 1, dstx = 0; srcx < src_w; ++srcx, ++sx.x, ++ax.x) { 00352 SKIPSMAlphaPixelType mcurrent(aa(ax) ? SKIPSMAlphaOne : SKIPSMAlphaZero); 00353 SKIPSMImagePixelType icurrent(aa(ax) ? SKIPSMImagePixelType(sa(sx)) : SKIPSMImageZero); 00354 if (evenX) { 00355 SKIPSMAlphaPixelType ap = asc1[dstx] + AMUL6(asc0[dstx]) + ascp[dstx]; 00356 asc1[dstx] = asc0[dstx] + ascp[dstx]; 00357 asc0[dstx] = asr1 + AMUL6(asr0) + asrp + mcurrent; 00358 asr1 = asr0 + asrp; 00359 asr0 = mcurrent; 00360 ap += asc0[dstx]; 00361 00362 SKIPSMImagePixelType ip = isc1[dstx] + IMUL6(isc0[dstx]) + iscp[dstx]; 00363 isc1[dstx] = isc0[dstx] + iscp[dstx]; 00364 isc0[dstx] = isr1 + IMUL6(isr0) + isrp + icurrent; 00365 isr1 = isr0 + isrp; 00366 isr0 = icurrent; 00367 if (ap) { 00368 ip += isc0[dstx]; 00369 ip /= ap; 00370 da.set(DestPixelType(ip), dx); 00371 daa.set(DestAlphaMax, dax); 00372 } else { 00373 da.set(DestImageZero, dx); 00374 daa.set(DestAlphaZero, dax); 00375 } 00376 00377 ++dx.x; 00378 ++dax.x; 00379 } 00380 else { 00381 asrp = mcurrent * 4; 00382 isrp = icurrent * 4; 00383 ++dstx; 00384 } 00385 evenX = !evenX; 00386 } 00387 00388 // Last entries in row 00389 if (!evenX) { 00390 // previous srcx was even 00391 ++dstx; 00392 00393 SKIPSMAlphaPixelType ap = asc1[dstx] + AMUL6(asc0[dstx]) + ascp[dstx]; 00394 asc1[dstx] = asc0[dstx] + ascp[dstx]; 00395 if (wraparound) { 00396 asc0[dstx] = asr1 + AMUL6(asr0) + (aa(ay) ? (SKIPSMAlphaOne * 4) : SKIPSMAlphaZero) + (aa(ay, Diff2D(1,0)) ? SKIPSMAlphaOne : SKIPSMAlphaZero); 00397 } else { 00398 asc0[dstx] = asr1 + AMUL6(asr0); 00399 } 00400 ap += asc0[dstx]; 00401 00402 SKIPSMImagePixelType ip = isc1[dstx] + IMUL6(isc0[dstx]) + iscp[dstx]; 00403 isc1[dstx] = isc0[dstx] + iscp[dstx]; 00404 if (wraparound) { 00405 isc0[dstx] = isr1 + IMUL6(isr0) + (aa(ay) ? (SKIPSMImagePixelType(sa(sy)) * 4) : SKIPSMImageZero) 00406 + (aa(ay, Diff2D(1,0)) ? SKIPSMImagePixelType(sa(sy, Diff2D(1,0))) : SKIPSMImageZero); 00407 } else { 00408 isc0[dstx] = isr1 + IMUL6(isr0); 00409 } 00410 if (ap) { 00411 ip += isc0[dstx]; 00412 //ip /= SKIPSMImagePixelType(ap); 00413 ip /= ap; 00414 da.set(DestPixelType(ip), dx); 00415 daa.set(DestAlphaMax, dax); 00416 } else { 00417 da.set(DestImageZero, dx); 00418 daa.set(DestAlphaZero, dax); 00419 } 00420 } 00421 else { 00422 // Previous srcx was odd 00423 SKIPSMAlphaPixelType ap = asc1[dstx] + AMUL6(asc0[dstx]) + ascp[dstx]; 00424 asc1[dstx] = asc0[dstx] + ascp[dstx]; 00425 if (wraparound) { 00426 asc0[dstx] = asr1 + AMUL6(asr0) + asrp + (aa(ay) ? SKIPSMAlphaOne : SKIPSMAlphaZero); 00427 } else { 00428 asc0[dstx] = asr1 + AMUL6(asr0) + asrp; 00429 } 00430 ap += asc0[dstx]; 00431 00432 SKIPSMImagePixelType ip = isc1[dstx] + IMUL6(isc0[dstx]) + iscp[dstx]; 00433 isc1[dstx] = isc0[dstx] + iscp[dstx]; 00434 if (wraparound) { 00435 isc0[dstx] = isr1 + IMUL6(isr0) + isrp + (aa(ay) ? SKIPSMImagePixelType(sa(sy)) : SKIPSMImageZero); 00436 } else { 00437 isc0[dstx] = isr1 + IMUL6(isr0) + isrp; 00438 } 00439 if (ap) { 00440 ip += isc0[dstx]; 00441 //ip /= SKIPSMImagePixelType(ap); 00442 ip /= ap; 00443 da.set(DestPixelType(ip), dx); 00444 daa.set(DestAlphaMax, dax); 00445 } else { 00446 da.set(DestImageZero, dx); 00447 daa.set(DestAlphaZero, dax); 00448 } 00449 } 00450 00451 ++dy.y; 00452 ++day.y; 00453 00454 } 00455 else { 00456 // First entry in odd-numbered row 00457 sx = sy; 00458 ax = ay; 00459 if (wraparound) { 00460 asr1 = asr0 + asrp; 00461 isr1 = isr0 + isrp; 00462 } 00463 asr0 = aa(ax) ? SKIPSMAlphaOne : SKIPSMAlphaZero; 00464 isr0 = aa(ax) ? SKIPSMImagePixelType(sa(sx)) : SKIPSMImageZero; 00465 // isc*[0] are never used 00466 ++sx.x; 00467 ++ax.x; 00468 00469 // Main entries in odd-numbered row 00470 for (evenX = false, srcx = 1, dstx = 0; srcx < src_w; ++srcx, ++sx.x, ++ax.x) { 00471 SKIPSMAlphaPixelType mcurrent(aa(ax) ? SKIPSMAlphaOne : SKIPSMAlphaZero); 00472 SKIPSMImagePixelType icurrent(aa(ax) ? SKIPSMImagePixelType(sa(sx)) : SKIPSMImageZero); 00473 if (evenX) { 00474 ascp[dstx] = (asr1 + AMUL6(asr0) + asrp + mcurrent) * 4; 00475 asr1 = asr0 + asrp; 00476 asr0 = mcurrent; 00477 iscp[dstx] = (isr1 + IMUL6(isr0) + isrp + icurrent) * 4; 00478 isr1 = isr0 + isrp; 00479 isr0 = icurrent; 00480 } 00481 else { 00482 asrp = mcurrent * 4; 00483 isrp = icurrent * 4; 00484 ++dstx; 00485 } 00486 evenX = !evenX; 00487 } 00488 00489 // Last entries in row 00490 if (!evenX) { 00491 // previous srcx was even 00492 ++dstx; 00493 if (wraparound) { 00494 ascp[dstx] = (asr1 + AMUL6(asr0) + (aa(ay) ? (SKIPSMAlphaOne * 4) : SKIPSMAlphaZero) 00495 + (aa(ay, Diff2D(1,0)) ? SKIPSMAlphaOne : SKIPSMAlphaZero) 00496 ) * 4; 00497 iscp[dstx] = (isr1 + IMUL6(isr0) + (aa(ay) ? (SKIPSMImagePixelType(sa(sy)) * 4) : SKIPSMImageZero) 00498 + (aa(ay, Diff2D(1,0)) ? SKIPSMImagePixelType(sa(sy, Diff2D(1,0))) : SKIPSMImageZero) 00499 ) * 4; 00500 } else { 00501 ascp[dstx] = (asr1 + AMUL6(asr0)) * 4; 00502 iscp[dstx] = (isr1 + IMUL6(isr0)) * 4; 00503 } 00504 } 00505 else { 00506 // previous srcx was odd 00507 if (wraparound) { 00508 ascp[dstx] = (asr1 + AMUL6(asr0) + asrp + (aa(ay) ? SKIPSMAlphaOne : SKIPSMAlphaZero)) * 4; 00509 iscp[dstx] = (isr1 + IMUL6(isr0) + isrp + (aa(ay) ? SKIPSMImagePixelType(sa(sy)) : SKIPSMImageZero)) * 4; 00510 } else { 00511 ascp[dstx] = (asr1 + AMUL6(asr0) + asrp) * 4; 00512 iscp[dstx] = (isr1 + IMUL6(isr0) + isrp) * 4; 00513 } 00514 } 00515 } 00516 evenY = !evenY; 00517 } 00518 } 00519 00520 // Last Rows 00521 { 00522 if (!evenY) { 00523 // Last srcy was even 00524 // odd row will set all iscp[] to zero 00525 // even row will do: 00526 //isc0[dstx] = 0; 00527 //isc1[dstx] = isc0[dstx] + 4*iscp[dstx] 00528 //out = isc1[dstx] + 6*isc0[dstx] + 4*iscp[dstx] + newisc0[dstx] 00529 for (dstx = 1, dx = dy, dax = day; dstx < (dst_w + 1); ++dstx, ++dx.x, ++dax.x) { 00530 SKIPSMAlphaPixelType ap = asc1[dstx] + AMUL6(asc0[dstx]); 00531 if (ap) { 00532 //SKIPSMImagePixelType ip = (isc1[dstx] + IMUL6(isc0[dstx])) / SKIPSMImagePixelType(ap); 00533 SKIPSMImagePixelType ip = (isc1[dstx] + IMUL6(isc0[dstx])) / ap; 00534 da.set(DestPixelType(ip), dx); 00535 daa.set(DestAlphaMax, dax); 00536 } else { 00537 da.set(DestImageZero, dx); 00538 daa.set(DestAlphaZero, dax); 00539 } 00540 } 00541 } 00542 else { 00543 // Last srcy was odd 00544 // even row will do: 00545 // isc0[dstx] = 0; 00546 // isc1[dstx] = isc0[dstx] + 4*iscp[dstx] 00547 // out = isc1[dstx] + 6*isc0[dstx] + 4*iscp[dstx] + newisc0[dstx] 00548 for (dstx = 1, dx = dy, dax = day; dstx < (dst_w + 1); ++dstx, ++dx.x, ++dax.x) { 00549 SKIPSMAlphaPixelType ap = asc1[dstx] + AMUL6(asc0[dstx]) + ascp[dstx]; 00550 if (ap) { 00551 //SKIPSMImagePixelType ip = (isc1[dstx] + IMUL6(isc0[dstx]) + iscp[dstx]) / SKIPSMImagePixelType(ap); 00552 SKIPSMImagePixelType ip = (isc1[dstx] + IMUL6(isc0[dstx]) + iscp[dstx]) / ap; 00553 da.set(DestPixelType(ip), dx); 00554 daa.set(DestAlphaMax, dax); 00555 } else { 00556 da.set(DestImageZero, dx); 00557 daa.set(DestAlphaZero, dax); 00558 } 00559 } 00560 } 00561 } 00562 00563 delete[] isc0; 00564 delete[] isc1; 00565 delete[] iscp; 00566 00567 delete[] asc0; 00568 delete[] asc1; 00569 delete[] ascp; 00570 00571 }; 00572 00573 // Version using argument object factories. 00574 template <typename SKIPSMImagePixelType, typename SKIPSMAlphaPixelType, 00575 typename SrcImageIterator, typename SrcAccessor, 00576 typename AlphaIterator, typename AlphaAccessor, 00577 typename DestImageIterator, typename DestAccessor, 00578 typename DestAlphaIterator, typename DestAlphaAccessor> 00579 inline void reduce(bool wraparound, 00580 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 00581 pair<AlphaIterator, AlphaAccessor> mask, 00582 triple<DestImageIterator, DestImageIterator, DestAccessor> dest, 00583 triple<DestAlphaIterator, DestAlphaIterator, DestAlphaAccessor> destMask) { 00584 reduce<SKIPSMImagePixelType, SKIPSMAlphaPixelType>(wraparound, 00585 src.first, src.second, src.third, 00586 mask.first, mask.second, 00587 dest.first, dest.second, dest.third, 00588 destMask.first, destMask.second, destMask.third); 00589 }; 00590 00594 template <typename SKIPSMImagePixelType, 00595 typename SrcImageIterator, typename SrcAccessor, 00596 typename DestImageIterator, typename DestAccessor> 00597 inline void reduce(bool wraparound, 00598 SrcImageIterator src_upperleft, 00599 SrcImageIterator src_lowerright, 00600 SrcAccessor sa, 00601 DestImageIterator dest_upperleft, 00602 DestImageIterator dest_lowerright, 00603 DestAccessor da) { 00604 00605 typedef typename DestAccessor::value_type DestPixelType; 00606 00607 int src_w = src_lowerright.x - src_upperleft.x; 00608 int src_h = src_lowerright.y - src_upperleft.y; 00609 int dst_w = dest_lowerright.x - dest_upperleft.x; 00610 //int dst_h = dest_lowerright.y - dest_upperleft.y; 00611 00612 vigra_precondition(src_w > 1 && src_h > 1, 00613 "src image too small in reduce"); 00614 00615 // State variables for source image pixel values 00616 SKIPSMImagePixelType isr0, isr1, isrp; 00617 SKIPSMImagePixelType *isc0 = new SKIPSMImagePixelType[dst_w + 1]; 00618 SKIPSMImagePixelType *isc1 = new SKIPSMImagePixelType[dst_w + 1]; 00619 SKIPSMImagePixelType *iscp = new SKIPSMImagePixelType[dst_w + 1]; 00620 00621 // Convenient constants 00622 const SKIPSMImagePixelType SKIPSMImageZero(NumericTraits<SKIPSMImagePixelType>::zero()); 00623 00624 DestImageIterator dy = dest_upperleft; 00625 DestImageIterator dx = dy; 00626 SrcImageIterator sy = src_upperleft; 00627 SrcImageIterator sx = sy; 00628 00629 bool evenY = true; 00630 bool evenX = true; 00631 int srcy = 0; 00632 int srcx = 0; 00633 //int dsty = 0; 00634 int dstx = 0; 00635 00636 // First row 00637 { 00638 if (wraparound) { 00639 isr0 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-2, 0))); 00640 isr1 = SKIPSMImageZero; 00641 isrp = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1, 0))) * 4; 00642 } else { 00643 isr0 = SKIPSMImagePixelType(sa(sy)); 00644 isr1 = SKIPSMImageZero; 00645 isrp = SKIPSMImagePixelType(sa(sy)) * 4; 00646 } 00647 00648 // Main pixels in first row 00649 for (sx = sy, evenX = true, srcx = 0, dstx = 0; srcx < src_w; ++srcx, ++sx.x) { 00650 SKIPSMImagePixelType icurrent(SKIPSMImagePixelType(sa(sx))); 00651 if (evenX) { 00652 isc0[dstx] = isr1 + IMUL6(isr0) + isrp + icurrent; 00653 isc1[dstx] = IMUL5(isc0[dstx]); 00654 isr1 = isr0 + isrp; 00655 isr0 = icurrent; 00656 } 00657 else { 00658 isrp = icurrent * 4; 00659 ++dstx; 00660 } 00661 evenX = !evenX; 00662 } 00663 00664 // Last entries in first row 00665 if (!evenX) { 00666 // previous srcx was even 00667 ++dstx; 00668 if (wraparound) { 00669 isc0[dstx] = isr1 + IMUL6(isr0) + (SKIPSMImagePixelType(sa(sy)) * 4) 00670 + SKIPSMImagePixelType(sa(sy, Diff2D(1,0))); 00671 isc1[dstx] = IMUL5(isc0[dstx]); 00672 } else { 00673 isc0[dstx] = isr1 + IMUL11(isr0); 00674 isc1[dstx] = IMUL5(isc0[dstx]); 00675 } 00676 } 00677 else { 00678 // previous srcx was odd 00679 if (wraparound) { 00680 isc0[dstx] = isr1 + IMUL6(isr0) + isrp + SKIPSMImagePixelType(sa(sy)); 00681 isc1[dstx] = IMUL5(isc0[dstx]); 00682 } else { 00683 isc0[dstx] = isr1 + IMUL6(isr0) + isrp + (isrp / 4); 00684 isc1[dstx] = IMUL5(isc0[dstx]); 00685 } 00686 } 00687 } 00688 ++sy.y; 00689 00690 // Main Rows 00691 { 00692 for (evenY = false, srcy = 1; srcy < src_h; ++srcy, ++sy.y) { 00693 00694 if (wraparound) { 00695 isr0 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-2,0))); 00696 isr1 = SKIPSMImageZero; 00697 isrp = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1,0))) * 4; 00698 } else { 00699 isr0 = SKIPSMImagePixelType(sa(sy)); 00700 isr1 = SKIPSMImageZero; 00701 isrp = SKIPSMImagePixelType(sa(sy)) * 4; 00702 } 00703 00704 if (evenY) { 00705 // Even-numbered row 00706 00707 // First entry in row 00708 sx = sy; 00709 isr1 = isr0 + isrp; 00710 isr0 = SKIPSMImagePixelType(sa(sx)); 00711 // isc*[0] are never used 00712 ++sx.x; 00713 dx = dy; 00714 00715 // Main entries in row 00716 for (evenX = false, srcx = 1, dstx = 0; srcx < src_w; ++srcx, ++sx.x) { 00717 SKIPSMImagePixelType icurrent(SKIPSMImagePixelType(sa(sx))); 00718 if (evenX) { 00719 SKIPSMImagePixelType ip = isc1[dstx] + IMUL6(isc0[dstx]) + iscp[dstx]; 00720 isc1[dstx] = isc0[dstx] + iscp[dstx]; 00721 isc0[dstx] = isr1 + IMUL6(isr0) + isrp + icurrent; 00722 isr1 = isr0 + isrp; 00723 isr0 = icurrent; 00724 ip += isc0[dstx]; 00725 ip /= 256; 00726 da.set(DestPixelType(ip), dx); 00727 ++dx.x; 00728 } 00729 else { 00730 isrp = icurrent * 4; 00731 ++dstx; 00732 } 00733 evenX = !evenX; 00734 } 00735 00736 // Last entries in row 00737 if (!evenX) { 00738 // previous srcx was even 00739 ++dstx; 00740 00741 SKIPSMImagePixelType ip = isc1[dstx] + IMUL6(isc0[dstx]) + iscp[dstx]; 00742 isc1[dstx] = isc0[dstx] + iscp[dstx]; 00743 if (wraparound) { 00744 isc0[dstx] = isr1 + IMUL6(isr0) + (SKIPSMImagePixelType(sa(sy)) * 4) 00745 + SKIPSMImagePixelType(sa(sy, Diff2D(1,0))); 00746 } else { 00747 isc0[dstx] = isr1 + IMUL11(isr0); 00748 } 00749 ip += isc0[dstx]; 00750 ip /= 256; 00751 da.set(DestPixelType(ip), dx); 00752 } 00753 else { 00754 // Previous srcx was odd 00755 SKIPSMImagePixelType ip = isc1[dstx] + IMUL6(isc0[dstx]) + iscp[dstx]; 00756 isc1[dstx] = isc0[dstx] + iscp[dstx]; 00757 if (wraparound) { 00758 isc0[dstx] = isr1 + IMUL6(isr0) + isrp + SKIPSMImagePixelType(sa(sy)); 00759 } else { 00760 isc0[dstx] = isr1 + IMUL6(isr0) + isrp + (isrp / 4); 00761 } 00762 ip += isc0[dstx]; 00763 ip /= 256; 00764 da.set(DestPixelType(ip), dx); 00765 } 00766 00767 ++dy.y; 00768 } 00769 else { 00770 // First entry in odd-numbered row 00771 sx = sy; 00772 isr1 = isr0 + isrp; 00773 isr0 = SKIPSMImagePixelType(sa(sx)); 00774 // isc*[0] are never used 00775 ++sx.x; 00776 00777 // Main entries in odd-numbered row 00778 for (evenX = false, srcx = 1, dstx = 0; srcx < src_w; ++srcx, ++sx.x) { 00779 SKIPSMImagePixelType icurrent(SKIPSMImagePixelType(sa(sx))); 00780 if (evenX) { 00781 iscp[dstx] = (isr1 + IMUL6(isr0) + isrp + icurrent) * 4; 00782 isr1 = isr0 + isrp; 00783 isr0 = icurrent; 00784 } 00785 else { 00786 isrp = icurrent * 4; 00787 ++dstx; 00788 } 00789 evenX = !evenX; 00790 } 00791 // Last entries in row 00792 if (!evenX) { 00793 // previous srcx was even 00794 ++dstx; 00795 if (wraparound) { 00796 iscp[dstx] = (isr1 + IMUL6(isr0) + (SKIPSMImagePixelType(sa(sy)) * 4) 00797 + SKIPSMImagePixelType(sa(sy, Diff2D(1,0))) 00798 ) * 4; 00799 } else { 00800 iscp[dstx] = (isr1 + IMUL11(isr0)) * 4; 00801 } 00802 } 00803 else { 00804 // previous srcx was odd 00805 if (wraparound) { 00806 iscp[dstx] = (isr1 + IMUL6(isr0) + isrp + SKIPSMImagePixelType(sa(sy))) * 4; 00807 } else { 00808 iscp[dstx] = (isr1 + IMUL6(isr0) + isrp + (isrp / 4)) * 4; 00809 } 00810 } 00811 } 00812 evenY = !evenY; 00813 } 00814 } 00815 00816 // Last Rows 00817 { 00818 if (!evenY) { 00819 // Last srcy was even 00820 // odd row will set all iscp[] to zero 00821 // even row will do: 00822 //isc0[dstx] = 0; 00823 //isc1[dstx] = isc0[dstx] + 4*iscp[dstx] 00824 //out = isc1[dstx] + 6*isc0[dstx] + 4*iscp[dstx] + newisc0[dstx] 00825 for (dstx = 1, dx = dy; dstx < (dst_w + 1); ++dstx, ++dx.x) { 00826 SKIPSMImagePixelType ip = (isc1[dstx] + IMUL11(isc0[dstx])) / 256; 00827 da.set(DestPixelType(ip), dx); 00828 } 00829 } 00830 else { 00831 // Last srcy was odd 00832 // even row will do: 00833 // isc0[dstx] = 0; 00834 // isc1[dstx] = isc0[dstx] + 4*iscp[dstx] 00835 // out = isc1[dstx] + 6*isc0[dstx] + 4*iscp[dstx] + newisc0[dstx] 00836 for (dstx = 1, dx = dy; dstx < (dst_w + 1); ++dstx, ++dx.x) { 00837 SKIPSMImagePixelType ip = (isc1[dstx] + IMUL6(isc0[dstx]) + iscp[dstx] + (iscp[dstx] / 4)) / 256; 00838 da.set(DestPixelType(ip), dx); 00839 } 00840 } 00841 } 00842 00843 delete[] isc0; 00844 delete[] isc1; 00845 delete[] iscp; 00846 00847 }; 00848 00849 // Version using argument object factories. 00850 template <typename SKIPSMImagePixelType, 00851 typename SrcImageIterator, typename SrcAccessor, 00852 typename DestImageIterator, typename DestAccessor> 00853 inline void reduce(bool wraparound, 00854 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 00855 triple<DestImageIterator, DestImageIterator, DestAccessor> dest) { 00856 reduce<SKIPSMImagePixelType>(wraparound, 00857 src.first, src.second, src.third, 00858 dest.first, dest.second, dest.third); 00859 }; 00860 00861 // SKIPSM update routine used when visiting a pixel in the top two rows 00862 // and the left two rows. 00863 #define SKIPSM_EXPAND(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11) \ 00864 current = SKIPSMImagePixelType(sa(sx)); \ 00865 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \ 00866 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \ 00867 out01 = sc0a[srcx]; \ 00868 out11 = sc0b[srcx]; \ 00869 sc1a[srcx] = sc0a[srcx]; \ 00870 sc1b[srcx] = sc0b[srcx]; \ 00871 sc0a[srcx] = sr1 + IMUL6(sr0) + current; \ 00872 sc0b[srcx] = (sr0 + current) * 4; \ 00873 sr1 = sr0; \ 00874 sr0 = current; \ 00875 out00 += sc0a[srcx]; \ 00876 out10 += sc0b[srcx]; \ 00877 out01 += sc0a[srcx]; \ 00878 out11 += sc0b[srcx]; \ 00879 out00 /= SKIPSMImagePixelType(SCALE_OUT00); \ 00880 out10 /= SKIPSMImagePixelType(SCALE_OUT10); \ 00881 out01 /= SKIPSMImagePixelType(SCALE_OUT01); \ 00882 out11 /= SKIPSMImagePixelType(SCALE_OUT11); \ 00883 da.set(cf(SKIPSMImagePixelType(da(dx)), out00), dx); \ 00884 ++dx.x; \ 00885 da.set(cf(SKIPSMImagePixelType(da(dx)), out10), dx); \ 00886 ++dx.x; \ 00887 da.set(cf(SKIPSMImagePixelType(da(dxx)), out01), dxx); \ 00888 ++dxx.x; \ 00889 da.set(cf(SKIPSMImagePixelType(da(dxx)), out11), dxx); \ 00890 ++dxx.x; 00891 00892 // SKIPSM update routine used when visiting a pixel in the main image body. 00893 // Same as above, but with hard-coded scaling factors. 00894 #define SKIPSM_EXPAND_SHIFT \ 00895 current = SKIPSMImagePixelType(sa(sx)); \ 00896 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \ 00897 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \ 00898 out01 = sc0a[srcx]; \ 00899 out11 = sc0b[srcx]; \ 00900 sc1a[srcx] = sc0a[srcx]; \ 00901 sc1b[srcx] = sc0b[srcx]; \ 00902 sc0a[srcx] = sr1 + IMUL6(sr0) + current; \ 00903 sc0b[srcx] = (sr0 + current) * 4; \ 00904 sr1 = sr0; \ 00905 sr0 = current; \ 00906 out00 += sc0a[srcx]; \ 00907 out10 += sc0b[srcx]; \ 00908 out01 += sc0a[srcx]; \ 00909 out11 += sc0b[srcx]; \ 00910 out00 /= 64; \ 00911 out10 /= 64; \ 00912 out01 /= 16; \ 00913 out11 /= 16; \ 00914 da.set(cf(SKIPSMImagePixelType(da(dx)), out00), dx); \ 00915 ++dx.x; \ 00916 da.set(cf(SKIPSMImagePixelType(da(dx)), out10), dx); \ 00917 ++dx.x; \ 00918 da.set(cf(SKIPSMImagePixelType(da(dxx)), out01), dxx); \ 00919 ++dxx.x; \ 00920 da.set(cf(SKIPSMImagePixelType(da(dxx)), out11), dxx); \ 00921 ++dxx.x; 00922 00923 // SKIPSM update routine used for the extra row under the main image body. 00924 #define SKIPSM_EXPAND_ROW_END(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11) \ 00925 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \ 00926 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \ 00927 out00 /= SKIPSMImagePixelType(SCALE_OUT00); \ 00928 out10 /= SKIPSMImagePixelType(SCALE_OUT10); \ 00929 da.set(cf(da(dx), out00), dx); \ 00930 ++dx.x; \ 00931 da.set(cf(da(dx), out10), dx); \ 00932 ++dx.x; \ 00933 if ((dst_h & 1) == 0) { \ 00934 out01 = sc0a[srcx]; \ 00935 out11 = sc0b[srcx]; \ 00936 out01 /= SKIPSMImagePixelType(SCALE_OUT01); \ 00937 out11 /= SKIPSMImagePixelType(SCALE_OUT11); \ 00938 da.set(cf(da(dxx), out01), dxx); \ 00939 ++dxx.x; \ 00940 da.set(cf(da(dxx), out11), dxx); \ 00941 ++dxx.x; \ 00942 } 00943 00944 // SKIPSM update routine used for the extra column to the right 00945 // of the main image body. 00946 #define SKIPSM_EXPAND_COLUMN_END(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11) \ 00947 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \ 00948 out01 = sc0a[srcx]; \ 00949 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \ 00950 out11 = sc0b[srcx]; \ 00951 sc1a[srcx] = sc0a[srcx]; \ 00952 sc1b[srcx] = sc0b[srcx]; \ 00953 sc0a[srcx] = sr1 + IMUL6(sr0); \ 00954 sc0b[srcx] = sr0 * 4; \ 00955 out00 += sc0a[srcx]; \ 00956 out01 += sc0a[srcx]; \ 00957 out00 /= SKIPSMImagePixelType(SCALE_OUT00); \ 00958 out01 /= SKIPSMImagePixelType(SCALE_OUT01); \ 00959 da.set(cf(da(dx), out00), dx); \ 00960 da.set(cf(da(dxx), out01), dxx); \ 00961 if ((dst_w & 1) == 0) { \ 00962 ++dx.x; \ 00963 ++dxx.x; \ 00964 out10 += sc0b[srcx]; \ 00965 out11 += sc0b[srcx]; \ 00966 out10 /= SKIPSMImagePixelType(SCALE_OUT10); \ 00967 out11 /= SKIPSMImagePixelType(SCALE_OUT11); \ 00968 da.set(cf(da(dx), out10), dx); \ 00969 da.set(cf(da(dxx), out11), dxx); \ 00970 } 00971 00972 // SKIPSM update routine used for the extra column to the right 00973 // of the main image body, with wraparound boundary conditions. 00974 #define SKIPSM_EXPAND_COLUMN_END_WRAPAROUND(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11) \ 00975 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \ 00976 out01 = sc0a[srcx]; \ 00977 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \ 00978 out11 = sc0b[srcx]; \ 00979 sc1a[srcx] = sc0a[srcx]; \ 00980 sc1b[srcx] = sc0b[srcx]; \ 00981 sc0a[srcx] = sr1 + IMUL6(sr0) + SKIPSMImagePixelType(sa(sy)); \ 00982 sc0b[srcx] = (sr0 + SKIPSMImagePixelType(sa(sy))) * 4; \ 00983 out00 += sc0a[srcx]; \ 00984 out01 += sc0a[srcx]; \ 00985 out00 /= SKIPSMImagePixelType(SCALE_OUT00); \ 00986 out01 /= SKIPSMImagePixelType(SCALE_OUT01); \ 00987 da.set(cf(da(dx), out00), dx); \ 00988 da.set(cf(da(dxx), out01), dxx); \ 00989 if ((dst_w & 1) == 0) { \ 00990 ++dx.x; \ 00991 ++dxx.x; \ 00992 out10 += sc0b[srcx]; \ 00993 out11 += sc0b[srcx]; \ 00994 out10 /= SKIPSMImagePixelType(SCALE_OUT10); \ 00995 out11 /= SKIPSMImagePixelType(SCALE_OUT11); \ 00996 da.set(cf(da(dx), out10), dx); \ 00997 da.set(cf(da(dxx), out11), dxx); \ 00998 } 00999 01000 // SKIPSM update routine for the extra column to the right 01001 // of the extra row under the main image body. 01002 #define SKIPSM_EXPAND_ROW_COLUMN_END(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11) \ 01003 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \ 01004 out00 /= SKIPSMImagePixelType(SCALE_OUT00); \ 01005 da.set(cf(da(dx), out00), dx); \ 01006 if ((dst_w & 1) == 0) { \ 01007 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \ 01008 out10 /= SKIPSMImagePixelType(SCALE_OUT10); \ 01009 ++dx.x; \ 01010 da.set(cf(da(dx), out10), dx); \ 01011 } \ 01012 if ((dst_h & 1) == 0) { \ 01013 out01 = sc0a[srcx]; \ 01014 out01 /= SKIPSMImagePixelType(SCALE_OUT01); \ 01015 da.set(cf(da(dxx), out01), dxx); \ 01016 if ((dst_w & 1) == 0) { \ 01017 out11 = sc0b[srcx]; \ 01018 out11 /= SKIPSMImagePixelType(SCALE_OUT11); \ 01019 ++dxx.x; \ 01020 da.set(cf(da(dxx), out11), dxx); \ 01021 } \ 01022 } 01023 01085 template <typename SKIPSMImagePixelType, 01086 typename SrcImageIterator, typename SrcAccessor, 01087 typename DestImageIterator, typename DestAccessor, 01088 typename CombineFunctor> 01089 void expand(bool add, bool wraparound, 01090 SrcImageIterator src_upperleft, 01091 SrcImageIterator src_lowerright, 01092 SrcAccessor sa, 01093 DestImageIterator dest_upperleft, 01094 DestImageIterator dest_lowerright, 01095 DestAccessor da, 01096 CombineFunctor cf) { 01097 01098 int src_w = src_lowerright.x - src_upperleft.x; 01099 int src_h = src_lowerright.y - src_upperleft.y; 01100 int dst_w = dest_lowerright.x - dest_upperleft.x; 01101 int dst_h = dest_lowerright.y - dest_upperleft.y; 01102 01103 // SKIPSM state variables 01104 SKIPSMImagePixelType current; 01105 SKIPSMImagePixelType out00, out10, out01, out11; 01106 SKIPSMImagePixelType sr0, sr1; 01107 SKIPSMImagePixelType *sc0a = new SKIPSMImagePixelType[src_w + 1]; 01108 SKIPSMImagePixelType *sc0b = new SKIPSMImagePixelType[src_w + 1]; 01109 SKIPSMImagePixelType *sc1a = new SKIPSMImagePixelType[src_w + 1]; 01110 SKIPSMImagePixelType *sc1b = new SKIPSMImagePixelType[src_w + 1]; 01111 01112 // Convenient constants 01113 const SKIPSMImagePixelType SKIPSMImageZero(NumericTraits<SKIPSMImagePixelType>::zero()); 01114 01115 DestImageIterator dy = dest_upperleft; 01116 DestImageIterator dyy = dest_upperleft; 01117 DestImageIterator dx = dy; 01118 DestImageIterator dxx = dyy; 01119 SrcImageIterator sy = src_upperleft; 01120 SrcImageIterator sx = sy; 01121 01122 int srcy = 0; 01123 int srcx = 0; 01124 //int dsty = 0; 01125 //int dstx = 0; 01126 01127 // First row 01128 { 01129 if (wraparound) { 01130 sr0 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1,0))); 01131 sr1 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-2,0))); 01132 } else { 01133 sr0 = SKIPSMImageZero; 01134 sr1 = SKIPSMImageZero; 01135 } 01136 01137 for (sx = sy, srcx = 0; srcx < src_w; ++srcx, ++sx.x) { 01138 current = SKIPSMImagePixelType(sa(sx)); 01139 sc0a[srcx] = sr1 + IMUL6(sr0) + current; 01140 sc0b[srcx] = (sr0 + current) * 4; 01141 sc1a[srcx] = SKIPSMImageZero; 01142 sc1b[srcx] = SKIPSMImageZero; 01143 sr1 = sr0; 01144 sr0 = current; 01145 } 01146 01147 // extra column at end of first row 01148 if (wraparound) { 01149 sc0a[srcx] = sr1 + IMUL6(sr0) + SKIPSMImagePixelType(sa(sy)); 01150 sc0b[srcx] = (sr0 + SKIPSMImagePixelType(sa(sy))) * 4; 01151 } else { 01152 sc0a[srcx] = sr1 + IMUL6(sr0); 01153 sc0b[srcx] = sr0 * 4; 01154 } 01155 sc1a[srcx] = SKIPSMImageZero; 01156 sc1b[srcx] = SKIPSMImageZero; 01157 } 01158 01159 // dy = row 0 01160 // dyy = row 1 01161 ++dyy.y; 01162 // sy = row 1 01163 srcy = 1; 01164 ++sy.y; 01165 01166 // Second row 01167 if (src_h > 1) { 01168 // First column 01169 srcx = 0; 01170 sx = sy; 01171 sr0 = SKIPSMImagePixelType(sa(sx)); 01172 if (wraparound) { 01173 sr1 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1,0))); 01174 } else { 01175 sr1 = SKIPSMImageZero; 01176 } 01177 // sc*[0] are irrelevant 01178 01179 srcx = 1; 01180 ++sx.x; 01181 dx = dy; 01182 dxx = dyy; 01183 01184 // Second column 01185 if (src_w > 1) { 01186 if (wraparound) { 01187 SKIPSM_EXPAND(56, 56, 16, 16) 01188 } else { 01189 SKIPSM_EXPAND(49, 56, 14, 16) 01190 } 01191 01192 // Main columns 01193 for (srcx = 2, ++sx.x; srcx < src_w; ++srcx, ++sx.x) { 01194 SKIPSM_EXPAND(56, 56, 16, 16) 01195 } 01196 01197 // extra column at end of second row 01198 if (wraparound) { 01199 SKIPSM_EXPAND_COLUMN_END_WRAPAROUND(56, 56, 16, 16) 01200 } else { 01201 SKIPSM_EXPAND_COLUMN_END(49, 28, 14, 8) 01202 } 01203 } 01204 else { 01205 // Math works out exactly the same for wraparound and no wraparound when src_w ==1 01206 SKIPSM_EXPAND_COLUMN_END(42, 28, 12, 8) 01207 } 01208 01209 } 01210 else { 01211 // No Second Row 01212 // First Column 01213 srcx = 0; 01214 sr0 = SKIPSMImageZero; 01215 sr1 = SKIPSMImageZero; 01216 01217 dx = dy; 01218 dxx = dyy; 01219 01220 if (src_w > 1) { 01221 // Second Column 01222 srcx = 1; 01223 if (wraparound) { 01224 SKIPSM_EXPAND_ROW_END(48, 48, 8, 8) 01225 } else { 01226 SKIPSM_EXPAND_ROW_END(42, 48, 7, 8) 01227 } 01228 01229 // Main columns 01230 for (srcx = 2; srcx < src_w; ++srcx) { 01231 SKIPSM_EXPAND_ROW_END(48, 48, 8, 8) 01232 } 01233 01234 // extra column at end of row 01235 if (wraparound) { 01236 SKIPSM_EXPAND_ROW_COLUMN_END(48, 48, 8, 8) 01237 } else { 01238 SKIPSM_EXPAND_ROW_COLUMN_END(42, 24, 7, 4) 01239 } 01240 } 01241 else { 01242 // No Second Column 01243 // dst_w, dst_h must be at least 2 01244 SKIPSM_EXPAND_ROW_COLUMN_END(36, 24, 6, 4) 01245 } 01246 01247 delete[] sc0a; 01248 delete[] sc0b; 01249 delete[] sc1a; 01250 delete[] sc1b; 01251 01252 return; 01253 } 01254 01255 // dy = row 2 01256 // dyy = row 3 01257 dy.y += 2; 01258 dyy.y += 2; 01259 // sy = row 2 01260 srcy = 2; 01261 ++sy.y; 01262 01263 // Main Rows 01264 for (srcy = 2, sx = sy; srcy < src_h; ++srcy, ++sy.y, dy.y += 2, dyy.y += 2) { 01265 // First column 01266 srcx = 0; 01267 sx = sy; 01268 sr0 = SKIPSMImagePixelType(sa(sx)); 01269 if (wraparound) { 01270 sr1 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1,0))); 01271 } else { 01272 sr1 = SKIPSMImageZero; 01273 } 01274 // sc*[0] are irrelvant 01275 01276 srcx = 1; 01277 ++sx.x; 01278 dx = dy; 01279 dxx = dyy; 01280 01281 // Second column 01282 if (src_w > 1) { 01283 if (wraparound) { 01284 SKIPSM_EXPAND_SHIFT 01285 } else { 01286 SKIPSM_EXPAND(56, 64, 14, 16) 01287 } 01288 01289 // Main columns 01290 for (srcx = 2, ++sx.x; srcx < src_w; ++srcx, ++sx.x) { 01291 //SKIPSM_EXPAND(64, 64, 16, 16) 01292 SKIPSM_EXPAND_SHIFT 01293 } 01294 01295 // extra column at end of row 01296 if (wraparound) { 01297 SKIPSM_EXPAND_COLUMN_END_WRAPAROUND(64, 64, 16, 16) 01298 } else { 01299 SKIPSM_EXPAND_COLUMN_END(56, 32, 14, 8) 01300 } 01301 } 01302 else { 01303 // No second column 01304 // dst_w must be at least 2 01305 // Math works out exactly the same for wraparound and no wraparound when src_w ==1 01306 SKIPSM_EXPAND_COLUMN_END(48, 32, 12, 8) 01307 } 01308 } 01309 01310 // Extra row at end 01311 { 01312 srcx = 0; 01313 sr0 = SKIPSMImageZero; 01314 sr1 = SKIPSMImageZero; 01315 01316 dx = dy; 01317 dxx = dyy; 01318 01319 if (src_w > 1) { 01320 // Second Column 01321 srcx = 1; 01322 if (wraparound) { 01323 SKIPSM_EXPAND_ROW_END(56, 56, 8, 8) 01324 } else { 01325 SKIPSM_EXPAND_ROW_END(49, 56, 7, 8) 01326 } 01327 01328 // Main columns 01329 for (srcx = 2; srcx < src_w; ++srcx) { 01330 SKIPSM_EXPAND_ROW_END(56, 56, 8, 8) 01331 } 01332 01333 // extra column at end of row 01334 if (wraparound) { 01335 SKIPSM_EXPAND_ROW_COLUMN_END(56, 56, 8, 8) 01336 } else { 01337 SKIPSM_EXPAND_ROW_COLUMN_END(49, 28, 7, 4) 01338 } 01339 } 01340 else { 01341 // No Second Column 01342 // dst_w, dst_h must be at least 2 01343 SKIPSM_EXPAND_ROW_COLUMN_END(42, 28, 6, 4) 01344 } 01345 } 01346 01347 delete[] sc0a; 01348 delete[] sc0b; 01349 delete[] sc1a; 01350 delete[] sc1b; 01351 01352 }; 01353 01354 // Functor that adds two values and de-promotes the result. 01355 // Used when collapsing a laplacian pyramid. 01356 // Explict fromPromote necessary to avoid overflow/underflow problems. 01357 template<typename T1, typename T2, typename T3> 01358 struct FromPromotePlusFunctorWrapper : public std::binary_function<T1, T2, T3> { 01359 inline T3 operator()(const T1 &a, const T2 &b) const { 01360 return NumericTraits<T3>::fromPromote(a + b); 01361 } 01362 }; 01363 01364 // Version using argument object factories. 01365 template <typename SKIPSMImagePixelType, 01366 typename SrcImageIterator, typename SrcAccessor, 01367 typename DestImageIterator, typename DestAccessor> 01368 inline void expand(bool add, bool wraparound, 01369 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01370 triple<DestImageIterator, DestImageIterator, DestAccessor> dest) { 01371 01372 typedef typename DestAccessor::value_type DestPixelType; 01373 01374 if (add) { 01375 expand<SKIPSMImagePixelType>(add, wraparound, 01376 src.first, src.second, src.third, 01377 dest.first, dest.second, dest.third, 01378 FromPromotePlusFunctorWrapper<DestPixelType, SKIPSMImagePixelType, DestPixelType>()); 01379 } 01380 else { 01381 expand<SKIPSMImagePixelType>(add, wraparound, 01382 src.first, src.second, src.third, 01383 dest.first, dest.second, dest.third, 01384 std::minus<SKIPSMImagePixelType>()); 01385 } 01386 01387 }; 01388 01389 #if 0 01390 // doesn't compile 01391 01393 template <typename SrcImageType, typename AlphaImageType, typename PyramidImageType, 01394 int PyramidIntegerBits, int PyramidFractionBits, 01395 typename SKIPSMImagePixelType, typename SKIPSMAlphaPixelType> 01396 vector<PyramidImageType*> *gaussianPyramid(unsigned int numLevels, 01397 bool wraparound, 01398 typename SrcImageType::const_traverser src_upperleft, 01399 typename SrcImageType::const_traverser src_lowerright, 01400 typename SrcImageType::ConstAccessor sa, 01401 typename AlphaImageType::const_traverser alpha_upperleft, 01402 typename AlphaImageType::ConstAccessor aa) { 01403 01404 vector<PyramidImageType*> *gp = new vector<PyramidImageType*>(); 01405 01406 // Size of pyramid level 0 01407 int w = src_lowerright.x - src_upperleft.x; 01408 int h = src_lowerright.y - src_upperleft.y; 01409 01410 // Pyramid level 0 01411 PyramidImageType *gp0 = new PyramidImageType(w, h); 01412 01413 // Copy src image into gp0, using fixed-point conversions. 01414 copyToPyramidImage<SrcImageType, PyramidImageType, PyramidIntegerBits, PyramidFractionBits>( 01415 src_upperleft, src_lowerright, sa, gp0->upperLeft(), gp0->accessor()); 01416 01417 gp->push_back(gp0); 01418 01419 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01420 cout << "Generating Gaussian pyramid: g0"; 01421 } 01422 01423 // Make remaining levels. 01424 PyramidImageType *lastGP = gp0; 01425 AlphaImageType *lastA = NULL; 01426 for (unsigned int l = 1; l < numLevels; l++) { 01427 01428 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01429 cout << " g" << l; 01430 cout.flush(); 01431 } 01432 01433 // Size of next level 01434 w = (w + 1) >> 1; 01435 h = (h + 1) >> 1; 01436 01437 // Next pyramid level 01438 PyramidImageType *gpn = new PyramidImageType(w, h); 01439 AlphaImageType *nextA = new AlphaImageType(w, h); 01440 01441 if (lastA == NULL) { 01442 reduce<SKIPSMImagePixelType, SKIPSMAlphaPixelType>(wraparound, 01443 srcImageRange(*lastGP), maskIter(alpha_upperleft, aa), 01444 destImageRange(*gpn), destImageRange(*nextA)); 01445 } else { 01446 reduce<SKIPSMImagePixelType, SKIPSMAlphaPixelType>(wraparound, 01447 srcImageRange(*lastGP), maskImage(*lastA), 01448 destImageRange(*gpn), destImageRange(*nextA)); 01449 } 01450 01451 gp->push_back(gpn); 01452 lastGP = gpn; 01453 delete lastA; 01454 lastA = nextA; 01455 } 01456 01457 delete lastA; 01458 01459 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01460 cout << endl; 01461 } 01462 01463 return gp; 01464 01465 }; 01466 01467 // Version using argument object factories. 01468 template <typename SrcImageType, typename AlphaImageType, typename PyramidImageType, 01469 int PyramidIntegerBits, int PyramidFractionBits, 01470 typename SKIPSMImagePixelType, typename SKIPSMAlphaPixelType> 01471 inline vector<PyramidImageType*> *gaussianPyramid(unsigned int numLevels, 01472 bool wraparound, 01473 triple<typename SrcImageType::const_traverser, typename SrcImageType::const_traverser, typename SrcImageType::ConstAccessor> src, 01474 pair<typename AlphaImageType::const_traverser, typename AlphaImageType::ConstAccessor> alpha) { 01475 return gaussianPyramid<SrcImageType, AlphaImageType, PyramidImageType, 01476 PyramidIntegerBits, PyramidFractionBits, 01477 SKIPSMImagePixelType, SKIPSMAlphaPixelType>( 01478 numLevels, wraparound, 01479 src.first, src.second, src.third, 01480 alpha.first, alpha.second); 01481 }; 01482 01484 template <typename SrcImageType, typename PyramidImageType, 01485 int PyramidIntegerBits, int PyramidFractionBits, 01486 typename SKIPSMImagePixelType> 01487 vector<PyramidImageType*> *gaussianPyramid(unsigned int numLevels, 01488 bool wraparound, 01489 typename SrcImageType::const_traverser src_upperleft, 01490 typename SrcImageType::const_traverser src_lowerright, 01491 typename SrcImageType::ConstAccessor sa) { 01492 01493 vector<PyramidImageType*> *gp = new vector<PyramidImageType*>(); 01494 01495 // Size of pyramid level 0 01496 int w = src_lowerright.x - src_upperleft.x; 01497 int h = src_lowerright.y - src_upperleft.y; 01498 01499 // Pyramid level 0 01500 PyramidImageType *gp0 = new PyramidImageType(w, h); 01501 01502 // Copy src image into gp0, using fixed-point conversions. 01503 copyToPyramidImage<SrcImageType, PyramidImageType, PyramidIntegerBits, PyramidFractionBits>( 01504 src_upperleft, src_lowerright, sa, 01505 gp0->upperLeft(), gp0->accessor()); 01506 01507 gp->push_back(gp0); 01508 01509 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01510 cout << "Generating Gaussian pyramid: g0"; 01511 } 01512 01513 // Make remaining levels. 01514 PyramidImageType *lastGP = gp0; 01515 for (unsigned int l = 1; l < numLevels; l++) { 01516 01517 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01518 cout << " g" << l; 01519 cout.flush(); 01520 } 01521 01522 // Size of next level 01523 w = (w + 1) >> 1; 01524 h = (h + 1) >> 1; 01525 01526 // Next pyramid level 01527 PyramidImageType *gpn = new PyramidImageType(w, h); 01528 01529 reduce<SKIPSMImagePixelType>(wraparound, srcImageRange(*lastGP), destImageRange(*gpn)); 01530 01531 gp->push_back(gpn); 01532 lastGP = gpn; 01533 } 01534 01535 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01536 cout << endl; 01537 } 01538 01539 return gp; 01540 }; 01541 01542 // Version using argument object factories. 01543 template <typename SrcImageType, typename PyramidImageType, 01544 int PyramidIntegerBits, int PyramidFractionBits, 01545 typename SKIPSMImagePixelType> 01546 inline vector<PyramidImageType*> *gaussianPyramid(unsigned int numLevels, 01547 bool wraparound, 01548 triple<typename SrcImageType::const_traverser, typename SrcImageType::const_traverser, typename SrcImageType::ConstAccessor> src) { 01549 return gaussianPyramid<SrcImageType, PyramidImageType, 01550 PyramidIntegerBits, PyramidFractionBits, SKIPSMImagePixelType>( 01551 numLevels, 01552 wraparound, 01553 src.first, src.second, src.third); 01554 }; 01555 01557 template <typename SrcImageType, typename AlphaImageType, typename PyramidImageType, 01558 int PyramidIntegerBits, int PyramidFractionBits, 01559 typename SKIPSMImagePixelType, typename SKIPSMAlphaPixelType> 01560 vector<PyramidImageType*> *laplacianPyramid(const char* exportName, unsigned int numLevels, 01561 bool wraparound, 01562 typename SrcImageType::const_traverser src_upperleft, 01563 typename SrcImageType::const_traverser src_lowerright, 01564 typename SrcImageType::ConstAccessor sa, 01565 typename AlphaImageType::const_traverser alpha_upperleft, 01566 typename AlphaImageType::ConstAccessor aa) { 01567 01568 // First create a Gaussian pyramid. 01569 vector <PyramidImageType*> *gp = 01570 gaussianPyramid<SrcImageType, AlphaImageType, PyramidImageType, 01571 PyramidIntegerBits, PyramidFractionBits, 01572 SKIPSMImagePixelType, SKIPSMAlphaPixelType>( 01573 numLevels, wraparound, 01574 src_upperleft, src_lowerright, sa, 01575 alpha_upperleft, aa); 01576 01577 //exportPyramid(gp, exportName); 01578 01579 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01580 cout << "Generating Laplacian pyramid:"; 01581 cout.flush(); 01582 } 01583 01584 // For each level, subtract the expansion of the next level. 01585 // Stop if there is no next level. 01586 for (unsigned int l = 0; l < (numLevels-1); l++) { 01587 01588 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01589 cout << " l" << l; 01590 cout.flush(); 01591 } 01592 01593 expand<SKIPSMImagePixelType>(false, wraparound, 01594 srcImageRange(*((*gp)[l+1])), 01595 destImageRange(*((*gp)[l]))); 01596 } 01597 01598 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01599 cout << " l" << (numLevels-1) << endl; 01600 } 01601 01602 //exportPyramid(gp, exportName); 01603 01604 return gp; 01605 }; 01606 01607 // Version using argument object factories. 01608 template <typename SrcImageType, typename AlphaImageType, typename PyramidImageType, 01609 int PyramidIntegerBits, int PyramidFractionBits, 01610 typename SKIPSMImagePixelType, typename SKIPSMAlphaPixelType> 01611 inline vector<PyramidImageType*> *laplacianPyramid(const char* exportName, unsigned int numLevels, 01612 bool wraparound, 01613 triple<typename SrcImageType::const_traverser, typename SrcImageType::const_traverser, typename SrcImageType::ConstAccessor> src, 01614 pair<typename AlphaImageType::const_traverser, typename AlphaImageType::ConstAccessor> alpha) { 01615 return laplacianPyramid<SrcImageType, AlphaImageType, PyramidImageType, 01616 PyramidIntegerBits, PyramidFractionBits, 01617 SKIPSMImagePixelType, SKIPSMAlphaPixelType>( 01618 exportName, 01619 numLevels, wraparound, 01620 src.first, src.second, src.third, 01621 alpha.first, alpha.second); 01622 }; 01623 01625 template <typename SKIPSMImagePixelType, typename PyramidImageType> 01626 void collapsePyramid(bool wraparound, vector<PyramidImageType*> *p) { 01627 01628 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01629 cout << "Collapsing Laplacian pyramid: " 01630 << "l" << p->size()-1; 01631 cout.flush(); 01632 } 01633 01634 // For each level, add the expansion of the next level. 01635 // Work backwards from the smallest level to the largest. 01636 for (int l = (p->size()-2); l >= 0; l--) { 01637 01638 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01639 cout << " l" << l; 01640 cout.flush(); 01641 } 01642 01643 expand<SKIPSMImagePixelType>(true, wraparound, 01644 srcImageRange(*((*p)[l+1])), 01645 destImageRange(*((*p)[l]))); 01646 } 01647 01648 if (Verbose > VERBOSE_PYRAMID_MESSAGES) { 01649 cout << endl; 01650 } 01651 01652 }; 01653 01654 // Export a scalar pyramid as a set of UINT16 tiff files. 01655 template <typename PyramidImageType> 01656 void exportPyramid(vector<PyramidImageType*> *v, const char *prefix, VigraTrueType) { 01657 typedef typename PyramidImageType::value_type PyramidValueType; 01658 01659 //for (unsigned int i = 0; i < (v->size() - 1); i++) { 01660 // // Clear all levels except last. 01661 // initImage(destImageRange(*((*v)[i])), NumericTraits<PyramidValueType>::zero()); 01662 //} 01663 //collapsePyramid(false, v); 01664 01665 for (unsigned int i = 0; i < v->size(); i++) { 01666 char filenameBuf[512]; 01667 snprintf(filenameBuf, 512, "%s%04u.tif", prefix, i); 01668 01669 // Rescale the pyramid values to fit in UINT16. 01670 UInt16Image usPyramid((*v)[i]->width(), (*v)[i]->height()); 01671 transformImage(srcImageRange(*((*v)[i])), destImage(usPyramid), 01672 linearRangeMapping(NumericTraits<PyramidValueType>::min(), 01673 NumericTraits<PyramidValueType>::max(), 01674 NumericTraits<UInt16>::min(), 01675 NumericTraits<UInt16>::max())); 01676 01677 ImageExportInfo info(filenameBuf); 01678 exportImage(srcImageRange(usPyramid), info); 01679 } 01680 }; 01681 01682 // Export a vector pyramid as a set of UINT16 tiff files. 01683 template <typename PyramidImageType> 01684 void exportPyramid(vector<PyramidImageType*> *v, const char *prefix, VigraFalseType) { 01685 typedef typename PyramidImageType::value_type PyramidVectorType; 01686 typedef typename PyramidVectorType::value_type PyramidValueType; 01687 01688 for (unsigned int i = 0; i < (v->size() - 1); i++) { 01689 // Clear all levels except last. 01690 initImage(destImageRange(*((*v)[i])), NumericTraits<PyramidValueType>::zero()); 01691 } 01692 collapsePyramid(false, v); 01693 01694 for (unsigned int i = 0; i < v->size(); i++) { 01695 char filenameBuf[512]; 01696 snprintf(filenameBuf, 512, "%s%04u.tif", prefix, i); 01697 01698 // Rescale the pyramid values to fit in UINT16. 01699 UInt16RGBImage usPyramid((*v)[i]->width(), (*v)[i]->height()); 01700 transformImage(srcImageRange(*((*v)[i])), destImage(usPyramid), 01701 linearRangeMapping(PyramidVectorType(NumericTraits<PyramidValueType>::min()), 01702 PyramidVectorType(NumericTraits<PyramidValueType>::max()), 01703 typename UInt16RGBImage::value_type(NumericTraits<UInt16>::min()), 01704 typename UInt16RGBImage::value_type(NumericTraits<UInt16>::max()))); 01705 01706 ImageExportInfo info(filenameBuf); 01707 exportImage(srcImageRange(usPyramid), info); 01708 } 01709 }; 01710 01711 // Export a pyramid as a set of UINT16 tiff files. 01712 template <typename PyramidImageType> 01713 void exportPyramid(vector<PyramidImageType*> *v, const char *prefix) { 01714 typedef typename NumericTraits<typename PyramidImageType::value_type>::isScalar pyramid_is_scalar; 01715 exportPyramid(v, prefix, pyramid_is_scalar()); 01716 }; 01717 #endif 01718 } // namespace enblend 01719 01720 #endif /* __PYRAMID_H__ */
