00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
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
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
00071
00072
00073 vigra_precondition((levels >= 1 && levels <= 29),
00074 "filterHalfWidth: levels outside of range [1,29]");
00075
00076
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
00198
00199 vigra_precondition(src_w > 1 && src_h > 1,
00200 "src image too small in reduce");
00201
00202
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
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
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
00236 int dstx = 0;
00237
00238
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
00278 if (!evenX) {
00279
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
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
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
00334
00335
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
00345 ++sx.x;
00346 ++ax.x;
00347 dx = dy;
00348 dax = day;
00349
00350
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
00389 if (!evenX) {
00390
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
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
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
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
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
00466 ++sx.x;
00467 ++ax.x;
00468
00469
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
00490 if (!evenX) {
00491
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
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
00521 {
00522 if (!evenY) {
00523
00524
00525
00526
00527
00528
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
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
00544
00545
00546
00547
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
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
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
00611
00612 vigra_precondition(src_w > 1 && src_h > 1,
00613 "src image too small in reduce");
00614
00615
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
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
00634 int dstx = 0;
00635
00636
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
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
00665 if (!evenX) {
00666
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
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
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
00706
00707
00708 sx = sy;
00709 isr1 = isr0 + isrp;
00710 isr0 = SKIPSMImagePixelType(sa(sx));
00711
00712 ++sx.x;
00713 dx = dy;
00714
00715
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
00737 if (!evenX) {
00738
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
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
00771 sx = sy;
00772 isr1 = isr0 + isrp;
00773 isr0 = SKIPSMImagePixelType(sa(sx));
00774
00775 ++sx.x;
00776
00777
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
00792 if (!evenX) {
00793
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
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
00817 {
00818 if (!evenY) {
00819
00820
00821
00822
00823
00824
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
00832
00833
00834
00835
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
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
00862
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
00893
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
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
00945
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
00973
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
01001
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
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
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
01125
01126
01127
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
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
01160
01161 ++dyy.y;
01162
01163 srcy = 1;
01164 ++sy.y;
01165
01166
01167 if (src_h > 1) {
01168
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
01178
01179 srcx = 1;
01180 ++sx.x;
01181 dx = dy;
01182 dxx = dyy;
01183
01184
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
01193 for (srcx = 2, ++sx.x; srcx < src_w; ++srcx, ++sx.x) {
01194 SKIPSM_EXPAND(56, 56, 16, 16)
01195 }
01196
01197
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
01206 SKIPSM_EXPAND_COLUMN_END(42, 28, 12, 8)
01207 }
01208
01209 }
01210 else {
01211
01212
01213 srcx = 0;
01214 sr0 = SKIPSMImageZero;
01215 sr1 = SKIPSMImageZero;
01216
01217 dx = dy;
01218 dxx = dyy;
01219
01220 if (src_w > 1) {
01221
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
01230 for (srcx = 2; srcx < src_w; ++srcx) {
01231 SKIPSM_EXPAND_ROW_END(48, 48, 8, 8)
01232 }
01233
01234
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
01243
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
01256
01257 dy.y += 2;
01258 dyy.y += 2;
01259
01260 srcy = 2;
01261 ++sy.y;
01262
01263
01264 for (srcy = 2, sx = sy; srcy < src_h; ++srcy, ++sy.y, dy.y += 2, dyy.y += 2) {
01265
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
01275
01276 srcx = 1;
01277 ++sx.x;
01278 dx = dy;
01279 dxx = dyy;
01280
01281
01282 if (src_w > 1) {
01283 if (wraparound) {
01284 SKIPSM_EXPAND_SHIFT
01285 } else {
01286 SKIPSM_EXPAND(56, 64, 14, 16)
01287 }
01288
01289
01290 for (srcx = 2, ++sx.x; srcx < src_w; ++srcx, ++sx.x) {
01291
01292 SKIPSM_EXPAND_SHIFT
01293 }
01294
01295
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
01304
01305
01306 SKIPSM_EXPAND_COLUMN_END(48, 32, 12, 8)
01307 }
01308 }
01309
01310
01311 {
01312 srcx = 0;
01313 sr0 = SKIPSMImageZero;
01314 sr1 = SKIPSMImageZero;
01315
01316 dx = dy;
01317 dxx = dyy;
01318
01319 if (src_w > 1) {
01320
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
01329 for (srcx = 2; srcx < src_w; ++srcx) {
01330 SKIPSM_EXPAND_ROW_END(56, 56, 8, 8)
01331 }
01332
01333
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
01342
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
01355
01356
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
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
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
01407 int w = src_lowerright.x - src_upperleft.x;
01408 int h = src_lowerright.y - src_upperleft.y;
01409
01410
01411 PyramidImageType *gp0 = new PyramidImageType(w, h);
01412
01413
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
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
01434 w = (w + 1) >> 1;
01435 h = (h + 1) >> 1;
01436
01437
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
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
01496 int w = src_lowerright.x - src_upperleft.x;
01497 int h = src_lowerright.y - src_upperleft.y;
01498
01499
01500 PyramidImageType *gp0 = new PyramidImageType(w, h);
01501
01502
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
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
01523 w = (w + 1) >> 1;
01524 h = (h + 1) >> 1;
01525
01526
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
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
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
01578
01579 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
01580 cout << "Generating Laplacian pyramid:";
01581 cout.flush();
01582 }
01583
01584
01585
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
01603
01604 return gp;
01605 };
01606
01607
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
01635
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
01655 template <typename PyramidImageType>
01656 void exportPyramid(vector<PyramidImageType*> *v, const char *prefix, VigraTrueType) {
01657 typedef typename PyramidImageType::value_type PyramidValueType;
01658
01659
01660
01661
01662
01663
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
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
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
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
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
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 }
01719
01720 #endif