32 #include <vigra/convolution.hxx>
33 #include <vigra/error.hxx>
34 #include <vigra/inspectimage.hxx>
35 #include <vigra/numerictraits.hxx>
36 #include <vigra/rgbvalue.hxx>
37 #include <vigra/sized_int.hxx>
38 #include <vigra/transformimage.hxx>
48 using vigra::linearRangeMapping;
49 using vigra::NumericTraits;
53 using vigra::UInt16Image;
54 using vigra::UInt16RGBImage;
58 #define IMUL6(A) (A * SKIPSMImagePixelType(6))
59 #define IMUL5(A) (A * SKIPSMImagePixelType(5))
60 #define IMUL11(A) (A * SKIPSMImagePixelType(11))
61 #define AMUL6(A) (A * SKIPSMAlphaPixelType(6))
68 template <
typename ImagePixelComponentType>
73 vigra_precondition((levels >= 1 && levels <= 29),
74 "filterHalfWidth: levels outside of range [1,29]");
77 int halfWidth = (levels == 1) ? 0 : ((1 << (levels+1)) - 4);
173 template <
typename SKIPSMImagePixelType,
typename SKIPSMAlphaPixelType,
174 typename SrcImageIterator,
typename SrcAccessor,
175 typename AlphaIterator,
typename AlphaAccessor,
176 typename DestImageIterator,
typename DestAccessor,
177 typename DestAlphaIterator,
typename DestAlphaAccessor>
179 SrcImageIterator src_upperleft,
180 SrcImageIterator src_lowerright,
182 AlphaIterator alpha_upperleft,
184 DestImageIterator dest_upperleft,
185 DestImageIterator dest_lowerright,
187 DestAlphaIterator dest_alpha_upperleft,
188 DestAlphaIterator dest_alpha_lowerright,
189 DestAlphaAccessor daa) {
191 typedef typename DestAccessor::value_type DestPixelType;
192 typedef typename DestAlphaAccessor::value_type DestAlphaPixelType;
194 int src_w = src_lowerright.x - src_upperleft.x;
195 int src_h = src_lowerright.y - src_upperleft.y;
196 int dst_w = dest_lowerright.x - dest_upperleft.x;
199 vigra_precondition(src_w > 1 && src_h > 1,
200 "src image too small in reduce");
203 SKIPSMImagePixelType isr0, isr1, isrp;
204 SKIPSMImagePixelType *isc0 =
new SKIPSMImagePixelType[dst_w + 1];
205 SKIPSMImagePixelType *isc1 =
new SKIPSMImagePixelType[dst_w + 1];
206 SKIPSMImagePixelType *iscp =
new SKIPSMImagePixelType[dst_w + 1];
209 SKIPSMAlphaPixelType asr0, asr1, asrp;
210 SKIPSMAlphaPixelType *asc0 =
new SKIPSMAlphaPixelType[dst_w + 1];
211 SKIPSMAlphaPixelType *asc1 =
new SKIPSMAlphaPixelType[dst_w + 1];
212 SKIPSMAlphaPixelType *ascp =
new SKIPSMAlphaPixelType[dst_w + 1];
215 const SKIPSMImagePixelType SKIPSMImageZero(NumericTraits<SKIPSMImagePixelType>::zero());
216 const SKIPSMAlphaPixelType SKIPSMAlphaZero(NumericTraits<SKIPSMAlphaPixelType>::zero());
217 const SKIPSMAlphaPixelType SKIPSMAlphaOne(NumericTraits<SKIPSMAlphaPixelType>::one());
218 const DestPixelType DestImageZero(NumericTraits<DestPixelType>::zero());
219 const DestAlphaPixelType DestAlphaZero(NumericTraits<DestAlphaPixelType>::zero());
222 DestImageIterator dy = dest_upperleft;
223 DestImageIterator dx = dy;
224 SrcImageIterator sy = src_upperleft;
225 SrcImageIterator sx = sy;
226 AlphaIterator ay = alpha_upperleft;
227 AlphaIterator ax = ay;
228 DestAlphaIterator day = dest_alpha_upperleft;
229 DestAlphaIterator dax = day;
241 asr0 = aa(ay, Diff2D(src_w-2, 0)) ? SKIPSMAlphaOne : SKIPSMAlphaZero;
242 asr1 = SKIPSMAlphaZero;
243 asrp = aa(ay, Diff2D(src_w-1, 0)) ? (SKIPSMAlphaOne * 4) : SKIPSMAlphaZero;
244 isr0 = aa(ay, Diff2D(src_w-2, 0)) ? SKIPSMImagePixelType(sa(sy, Diff2D(src_w-2, 0))) : SKIPSMImageZero;
245 isr1 = SKIPSMImageZero;
246 isrp = aa(ay, Diff2D(src_w-1, 0)) ? (SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1, 0))) * 4) : SKIPSMImageZero;
248 asr0 = SKIPSMAlphaZero;
249 asr1 = SKIPSMAlphaZero;
250 asrp = SKIPSMAlphaZero;
251 isr0 = SKIPSMImageZero;
252 isr1 = SKIPSMImageZero;
253 isrp = SKIPSMImageZero;
256 for (sx = sy, ax = ay, evenX =
true, srcx = 0, dstx = 0; srcx < src_w; ++srcx, ++sx.x, ++ax.x) {
257 SKIPSMAlphaPixelType mcurrent(aa(ax) ? SKIPSMAlphaOne : SKIPSMAlphaZero);
258 SKIPSMImagePixelType icurrent(aa(ax) ? SKIPSMImagePixelType(sa(sx)) : SKIPSMImageZero);
260 asc1[dstx] = SKIPSMAlphaZero;
261 asc0[dstx] = asr1 +
AMUL6(asr0) + asrp + mcurrent;
264 isc1[dstx] = SKIPSMImageZero;
265 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp + icurrent;
282 asc1[dstx] = SKIPSMAlphaZero;
283 asc0[dstx] = asr1 +
AMUL6(asr0) + (aa(ay) ? (SKIPSMAlphaOne * 4) : SKIPSMAlphaZero) + (aa(ay, Diff2D(1,0)) ? SKIPSMAlphaOne : SKIPSMAlphaZero);
284 isc1[dstx] = SKIPSMImageZero;
285 isc0[dstx] = isr1 +
IMUL6(isr0) + (aa(ay) ? (SKIPSMImagePixelType(sa(sy)) * 4) : SKIPSMImageZero)
286 + (aa(ay, Diff2D(1,0)) ? SKIPSMImagePixelType(sa(sy, Diff2D(1,0))) : SKIPSMImageZero);
288 asc1[dstx] = SKIPSMAlphaZero;
289 asc0[dstx] = asr1 +
AMUL6(asr0);
290 isc1[dstx] = SKIPSMImageZero;
291 isc0[dstx] = isr1 +
IMUL6(isr0);
297 asc1[dstx] = SKIPSMAlphaZero;
298 asc0[dstx] = asr1 +
AMUL6(asr0) + asrp + (aa(ay) ? SKIPSMAlphaOne : SKIPSMAlphaZero);
299 isc1[dstx] = SKIPSMImageZero;
300 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp + (aa(ay) ? SKIPSMImagePixelType(sa(sy)) : SKIPSMImageZero);
302 asc1[dstx] = SKIPSMAlphaZero;
303 asc0[dstx] = asr1 +
AMUL6(asr0) + asrp;
304 isc1[dstx] = SKIPSMImageZero;
305 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp;
314 for (evenY =
false, srcy = 1; srcy < src_h; ++srcy, ++sy.y, ++ay.y) {
317 asr0 = aa(ay, Diff2D(src_w-2, 0)) ? SKIPSMAlphaOne : SKIPSMAlphaZero;
318 asr1 = SKIPSMAlphaZero;
319 asrp = aa(ay, Diff2D(src_w-1, 0)) ? (SKIPSMAlphaOne * 4) : SKIPSMAlphaZero;
320 isr0 = aa(ay, Diff2D(src_w-2, 0)) ? SKIPSMImagePixelType(sa(sy, Diff2D(src_w-2,0))) : SKIPSMImageZero;
321 isr1 = SKIPSMImageZero;
322 isrp = aa(ay, Diff2D(src_w-1, 0)) ? (SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1,0))) * 4) : SKIPSMImageZero;
324 asr0 = SKIPSMAlphaZero;
325 asr1 = SKIPSMAlphaZero;
326 asrp = SKIPSMAlphaZero;
327 isr0 = SKIPSMImageZero;
328 isr1 = SKIPSMImageZero;
329 isrp = SKIPSMImageZero;
342 asr0 = aa(ax) ? SKIPSMAlphaOne : SKIPSMAlphaZero;
343 isr0 = aa(ax) ? SKIPSMImagePixelType(sa(sx)) : SKIPSMImageZero;
351 for (evenX =
false, srcx = 1, dstx = 0; srcx < src_w; ++srcx, ++sx.x, ++ax.x) {
352 SKIPSMAlphaPixelType mcurrent(aa(ax) ? SKIPSMAlphaOne : SKIPSMAlphaZero);
353 SKIPSMImagePixelType icurrent(aa(ax) ? SKIPSMImagePixelType(sa(sx)) : SKIPSMImageZero);
355 SKIPSMAlphaPixelType ap = asc1[dstx] +
AMUL6(asc0[dstx]) + ascp[dstx];
356 asc1[dstx] = asc0[dstx] + ascp[dstx];
357 asc0[dstx] = asr1 +
AMUL6(asr0) + asrp + mcurrent;
362 SKIPSMImagePixelType ip = isc1[dstx] +
IMUL6(isc0[dstx]) + iscp[dstx];
363 isc1[dstx] = isc0[dstx] + iscp[dstx];
364 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp + icurrent;
370 da.set(DestPixelType(ip), dx);
371 daa.set(DestAlphaMax, dax);
373 da.set(DestImageZero, dx);
374 daa.set(DestAlphaZero, dax);
393 SKIPSMAlphaPixelType ap = asc1[dstx] +
AMUL6(asc0[dstx]) + ascp[dstx];
394 asc1[dstx] = asc0[dstx] + ascp[dstx];
396 asc0[dstx] = asr1 +
AMUL6(asr0) + (aa(ay) ? (SKIPSMAlphaOne * 4) : SKIPSMAlphaZero) + (aa(ay, Diff2D(1,0)) ? SKIPSMAlphaOne : SKIPSMAlphaZero);
398 asc0[dstx] = asr1 +
AMUL6(asr0);
402 SKIPSMImagePixelType ip = isc1[dstx] +
IMUL6(isc0[dstx]) + iscp[dstx];
403 isc1[dstx] = isc0[dstx] + iscp[dstx];
405 isc0[dstx] = isr1 +
IMUL6(isr0) + (aa(ay) ? (SKIPSMImagePixelType(sa(sy)) * 4) : SKIPSMImageZero)
406 + (aa(ay, Diff2D(1,0)) ? SKIPSMImagePixelType(sa(sy, Diff2D(1,0))) : SKIPSMImageZero);
408 isc0[dstx] = isr1 +
IMUL6(isr0);
414 da.set(DestPixelType(ip), dx);
415 daa.set(DestAlphaMax, dax);
417 da.set(DestImageZero, dx);
418 daa.set(DestAlphaZero, dax);
423 SKIPSMAlphaPixelType ap = asc1[dstx] +
AMUL6(asc0[dstx]) + ascp[dstx];
424 asc1[dstx] = asc0[dstx] + ascp[dstx];
426 asc0[dstx] = asr1 +
AMUL6(asr0) + asrp + (aa(ay) ? SKIPSMAlphaOne : SKIPSMAlphaZero);
428 asc0[dstx] = asr1 +
AMUL6(asr0) + asrp;
432 SKIPSMImagePixelType ip = isc1[dstx] +
IMUL6(isc0[dstx]) + iscp[dstx];
433 isc1[dstx] = isc0[dstx] + iscp[dstx];
435 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp + (aa(ay) ? SKIPSMImagePixelType(sa(sy)) : SKIPSMImageZero);
437 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp;
443 da.set(DestPixelType(ip), dx);
444 daa.set(DestAlphaMax, dax);
446 da.set(DestImageZero, dx);
447 daa.set(DestAlphaZero, dax);
463 asr0 = aa(ax) ? SKIPSMAlphaOne : SKIPSMAlphaZero;
464 isr0 = aa(ax) ? SKIPSMImagePixelType(sa(sx)) : SKIPSMImageZero;
470 for (evenX =
false, srcx = 1, dstx = 0; srcx < src_w; ++srcx, ++sx.x, ++ax.x) {
471 SKIPSMAlphaPixelType mcurrent(aa(ax) ? SKIPSMAlphaOne : SKIPSMAlphaZero);
472 SKIPSMImagePixelType icurrent(aa(ax) ? SKIPSMImagePixelType(sa(sx)) : SKIPSMImageZero);
474 ascp[dstx] = (asr1 +
AMUL6(asr0) + asrp + mcurrent) * 4;
477 iscp[dstx] = (isr1 +
IMUL6(isr0) + isrp + icurrent) * 4;
494 ascp[dstx] = (asr1 +
AMUL6(asr0) + (aa(ay) ? (SKIPSMAlphaOne * 4) : SKIPSMAlphaZero)
495 + (aa(ay, Diff2D(1,0)) ? SKIPSMAlphaOne : SKIPSMAlphaZero)
497 iscp[dstx] = (isr1 +
IMUL6(isr0) + (aa(ay) ? (SKIPSMImagePixelType(sa(sy)) * 4) : SKIPSMImageZero)
498 + (aa(ay, Diff2D(1,0)) ? SKIPSMImagePixelType(sa(sy, Diff2D(1,0))) : SKIPSMImageZero)
501 ascp[dstx] = (asr1 +
AMUL6(asr0)) * 4;
502 iscp[dstx] = (isr1 +
IMUL6(isr0)) * 4;
508 ascp[dstx] = (asr1 +
AMUL6(asr0) + asrp + (aa(ay) ? SKIPSMAlphaOne : SKIPSMAlphaZero)) * 4;
509 iscp[dstx] = (isr1 +
IMUL6(isr0) + isrp + (aa(ay) ? SKIPSMImagePixelType(sa(sy)) : SKIPSMImageZero)) * 4;
511 ascp[dstx] = (asr1 +
AMUL6(asr0) + asrp) * 4;
512 iscp[dstx] = (isr1 +
IMUL6(isr0) + isrp) * 4;
529 for (dstx = 1, dx = dy, dax = day; dstx < (dst_w + 1); ++dstx, ++dx.x, ++dax.x) {
530 SKIPSMAlphaPixelType ap = asc1[dstx] +
AMUL6(asc0[dstx]);
533 SKIPSMImagePixelType ip = (isc1[dstx] +
IMUL6(isc0[dstx])) / ap;
534 da.set(DestPixelType(ip), dx);
535 daa.set(DestAlphaMax, dax);
537 da.set(DestImageZero, dx);
538 daa.set(DestAlphaZero, dax);
548 for (dstx = 1, dx = dy, dax = day; dstx < (dst_w + 1); ++dstx, ++dx.x, ++dax.x) {
549 SKIPSMAlphaPixelType ap = asc1[dstx] +
AMUL6(asc0[dstx]) + ascp[dstx];
552 SKIPSMImagePixelType ip = (isc1[dstx] +
IMUL6(isc0[dstx]) + iscp[dstx]) / ap;
553 da.set(DestPixelType(ip), dx);
554 daa.set(DestAlphaMax, dax);
556 da.set(DestImageZero, dx);
557 daa.set(DestAlphaZero, dax);
574 template <
typename SKIPSMImagePixelType,
typename SKIPSMAlphaPixelType,
575 typename SrcImageIterator,
typename SrcAccessor,
576 typename AlphaIterator,
typename AlphaAccessor,
577 typename DestImageIterator,
typename DestAccessor,
578 typename DestAlphaIterator,
typename DestAlphaAccessor>
580 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
581 pair<AlphaIterator, AlphaAccessor> mask,
582 triple<DestImageIterator, DestImageIterator, DestAccessor> dest,
583 triple<DestAlphaIterator, DestAlphaIterator, DestAlphaAccessor>
destMask) {
584 reduce<SKIPSMImagePixelType, SKIPSMAlphaPixelType>(wraparound,
585 src.first, src.second, src.third,
586 mask.first, mask.second,
587 dest.first, dest.second, dest.third,
588 destMask.first, destMask.second, destMask.third);
594 template <
typename SKIPSMImagePixelType,
595 typename SrcImageIterator,
typename SrcAccessor,
596 typename DestImageIterator,
typename DestAccessor>
598 SrcImageIterator src_upperleft,
599 SrcImageIterator src_lowerright,
601 DestImageIterator dest_upperleft,
602 DestImageIterator dest_lowerright,
605 typedef typename DestAccessor::value_type DestPixelType;
607 int src_w = src_lowerright.x - src_upperleft.x;
608 int src_h = src_lowerright.y - src_upperleft.y;
609 int dst_w = dest_lowerright.x - dest_upperleft.x;
612 vigra_precondition(src_w > 1 && src_h > 1,
613 "src image too small in reduce");
616 SKIPSMImagePixelType isr0, isr1, isrp;
617 SKIPSMImagePixelType *isc0 =
new SKIPSMImagePixelType[dst_w + 1];
618 SKIPSMImagePixelType *isc1 =
new SKIPSMImagePixelType[dst_w + 1];
619 SKIPSMImagePixelType *iscp =
new SKIPSMImagePixelType[dst_w + 1];
622 const SKIPSMImagePixelType SKIPSMImageZero(NumericTraits<SKIPSMImagePixelType>::zero());
624 DestImageIterator dy = dest_upperleft;
625 DestImageIterator dx = dy;
626 SrcImageIterator sy = src_upperleft;
627 SrcImageIterator sx = sy;
639 isr0 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-2, 0)));
640 isr1 = SKIPSMImageZero;
641 isrp = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1, 0))) * 4;
643 isr0 = SKIPSMImagePixelType(sa(sy));
644 isr1 = SKIPSMImageZero;
645 isrp = SKIPSMImagePixelType(sa(sy)) * 4;
649 for (sx = sy, evenX =
true, srcx = 0, dstx = 0; srcx < src_w; ++srcx, ++sx.x) {
650 SKIPSMImagePixelType icurrent(SKIPSMImagePixelType(sa(sx)));
652 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp + icurrent;
653 isc1[dstx] =
IMUL5(isc0[dstx]);
669 isc0[dstx] = isr1 +
IMUL6(isr0) + (SKIPSMImagePixelType(sa(sy)) * 4)
670 + SKIPSMImagePixelType(sa(sy, Diff2D(1,0)));
671 isc1[dstx] =
IMUL5(isc0[dstx]);
673 isc0[dstx] = isr1 +
IMUL11(isr0);
674 isc1[dstx] =
IMUL5(isc0[dstx]);
680 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp + SKIPSMImagePixelType(sa(sy));
681 isc1[dstx] =
IMUL5(isc0[dstx]);
683 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp + (isrp / 4);
684 isc1[dstx] =
IMUL5(isc0[dstx]);
692 for (evenY =
false, srcy = 1; srcy < src_h; ++srcy, ++sy.y) {
695 isr0 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-2,0)));
696 isr1 = SKIPSMImageZero;
697 isrp = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1,0))) * 4;
699 isr0 = SKIPSMImagePixelType(sa(sy));
700 isr1 = SKIPSMImageZero;
701 isrp = SKIPSMImagePixelType(sa(sy)) * 4;
710 isr0 = SKIPSMImagePixelType(sa(sx));
716 for (evenX =
false, srcx = 1, dstx = 0; srcx < src_w; ++srcx, ++sx.x) {
717 SKIPSMImagePixelType icurrent(SKIPSMImagePixelType(sa(sx)));
719 SKIPSMImagePixelType ip = isc1[dstx] +
IMUL6(isc0[dstx]) + iscp[dstx];
720 isc1[dstx] = isc0[dstx] + iscp[dstx];
721 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp + icurrent;
726 da.set(DestPixelType(ip), dx);
741 SKIPSMImagePixelType ip = isc1[dstx] +
IMUL6(isc0[dstx]) + iscp[dstx];
742 isc1[dstx] = isc0[dstx] + iscp[dstx];
744 isc0[dstx] = isr1 +
IMUL6(isr0) + (SKIPSMImagePixelType(sa(sy)) * 4)
745 + SKIPSMImagePixelType(sa(sy, Diff2D(1,0)));
747 isc0[dstx] = isr1 +
IMUL11(isr0);
751 da.set(DestPixelType(ip), dx);
755 SKIPSMImagePixelType ip = isc1[dstx] +
IMUL6(isc0[dstx]) + iscp[dstx];
756 isc1[dstx] = isc0[dstx] + iscp[dstx];
758 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp + SKIPSMImagePixelType(sa(sy));
760 isc0[dstx] = isr1 +
IMUL6(isr0) + isrp + (isrp / 4);
764 da.set(DestPixelType(ip), dx);
773 isr0 = SKIPSMImagePixelType(sa(sx));
778 for (evenX =
false, srcx = 1, dstx = 0; srcx < src_w; ++srcx, ++sx.x) {
779 SKIPSMImagePixelType icurrent(SKIPSMImagePixelType(sa(sx)));
781 iscp[dstx] = (isr1 +
IMUL6(isr0) + isrp + icurrent) * 4;
796 iscp[dstx] = (isr1 +
IMUL6(isr0) + (SKIPSMImagePixelType(sa(sy)) * 4)
797 + SKIPSMImagePixelType(sa(sy, Diff2D(1,0)))
800 iscp[dstx] = (isr1 +
IMUL11(isr0)) * 4;
806 iscp[dstx] = (isr1 +
IMUL6(isr0) + isrp + SKIPSMImagePixelType(sa(sy))) * 4;
808 iscp[dstx] = (isr1 +
IMUL6(isr0) + isrp + (isrp / 4)) * 4;
825 for (dstx = 1, dx = dy; dstx < (dst_w + 1); ++dstx, ++dx.x) {
826 SKIPSMImagePixelType ip = (isc1[dstx] +
IMUL11(isc0[dstx])) / 256;
827 da.set(DestPixelType(ip), dx);
836 for (dstx = 1, dx = dy; dstx < (dst_w + 1); ++dstx, ++dx.x) {
837 SKIPSMImagePixelType ip = (isc1[dstx] +
IMUL6(isc0[dstx]) + iscp[dstx] + (iscp[dstx] / 4)) / 256;
838 da.set(DestPixelType(ip), dx);
850 template <
typename SKIPSMImagePixelType,
851 typename SrcImageIterator,
typename SrcAccessor,
852 typename DestImageIterator,
typename DestAccessor>
854 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
855 triple<DestImageIterator, DestImageIterator, DestAccessor> dest) {
856 reduce<SKIPSMImagePixelType>(wraparound,
857 src.first, src.second, src.third,
858 dest.first, dest.second, dest.third);
863 #define SKIPSM_EXPAND(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11) \
864 current = SKIPSMImagePixelType(sa(sx)); \
865 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \
866 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \
867 out01 = sc0a[srcx]; \
868 out11 = sc0b[srcx]; \
869 sc1a[srcx] = sc0a[srcx]; \
870 sc1b[srcx] = sc0b[srcx]; \
871 sc0a[srcx] = sr1 + IMUL6(sr0) + current; \
872 sc0b[srcx] = (sr0 + current) * 4; \
875 out00 += sc0a[srcx]; \
876 out10 += sc0b[srcx]; \
877 out01 += sc0a[srcx]; \
878 out11 += sc0b[srcx]; \
879 out00 /= SKIPSMImagePixelType(SCALE_OUT00); \
880 out10 /= SKIPSMImagePixelType(SCALE_OUT10); \
881 out01 /= SKIPSMImagePixelType(SCALE_OUT01); \
882 out11 /= SKIPSMImagePixelType(SCALE_OUT11); \
883 da.set(cf(SKIPSMImagePixelType(da(dx)), out00), dx); \
885 da.set(cf(SKIPSMImagePixelType(da(dx)), out10), dx); \
887 da.set(cf(SKIPSMImagePixelType(da(dxx)), out01), dxx); \
889 da.set(cf(SKIPSMImagePixelType(da(dxx)), out11), dxx); \
894 #define SKIPSM_EXPAND_SHIFT \
895 current = SKIPSMImagePixelType(sa(sx)); \
896 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \
897 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \
898 out01 = sc0a[srcx]; \
899 out11 = sc0b[srcx]; \
900 sc1a[srcx] = sc0a[srcx]; \
901 sc1b[srcx] = sc0b[srcx]; \
902 sc0a[srcx] = sr1 + IMUL6(sr0) + current; \
903 sc0b[srcx] = (sr0 + current) * 4; \
906 out00 += sc0a[srcx]; \
907 out10 += sc0b[srcx]; \
908 out01 += sc0a[srcx]; \
909 out11 += sc0b[srcx]; \
914 da.set(cf(SKIPSMImagePixelType(da(dx)), out00), dx); \
916 da.set(cf(SKIPSMImagePixelType(da(dx)), out10), dx); \
918 da.set(cf(SKIPSMImagePixelType(da(dxx)), out01), dxx); \
920 da.set(cf(SKIPSMImagePixelType(da(dxx)), out11), dxx); \
924 #define SKIPSM_EXPAND_ROW_END(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11) \
925 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \
926 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \
927 out00 /= SKIPSMImagePixelType(SCALE_OUT00); \
928 out10 /= SKIPSMImagePixelType(SCALE_OUT10); \
929 da.set(cf(da(dx), out00), dx); \
931 da.set(cf(da(dx), out10), dx); \
933 if ((dst_h & 1) == 0) { \
934 out01 = sc0a[srcx]; \
935 out11 = sc0b[srcx]; \
936 out01 /= SKIPSMImagePixelType(SCALE_OUT01); \
937 out11 /= SKIPSMImagePixelType(SCALE_OUT11); \
938 da.set(cf(da(dxx), out01), dxx); \
940 da.set(cf(da(dxx), out11), dxx); \
946 #define SKIPSM_EXPAND_COLUMN_END(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11) \
947 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \
948 out01 = sc0a[srcx]; \
949 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \
950 out11 = sc0b[srcx]; \
951 sc1a[srcx] = sc0a[srcx]; \
952 sc1b[srcx] = sc0b[srcx]; \
953 sc0a[srcx] = sr1 + IMUL6(sr0); \
954 sc0b[srcx] = sr0 * 4; \
955 out00 += sc0a[srcx]; \
956 out01 += sc0a[srcx]; \
957 out00 /= SKIPSMImagePixelType(SCALE_OUT00); \
958 out01 /= SKIPSMImagePixelType(SCALE_OUT01); \
959 da.set(cf(da(dx), out00), dx); \
960 da.set(cf(da(dxx), out01), dxx); \
961 if ((dst_w & 1) == 0) { \
964 out10 += sc0b[srcx]; \
965 out11 += sc0b[srcx]; \
966 out10 /= SKIPSMImagePixelType(SCALE_OUT10); \
967 out11 /= SKIPSMImagePixelType(SCALE_OUT11); \
968 da.set(cf(da(dx), out10), dx); \
969 da.set(cf(da(dxx), out11), dxx); \
974 #define SKIPSM_EXPAND_COLUMN_END_WRAPAROUND(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11) \
975 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \
976 out01 = sc0a[srcx]; \
977 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \
978 out11 = sc0b[srcx]; \
979 sc1a[srcx] = sc0a[srcx]; \
980 sc1b[srcx] = sc0b[srcx]; \
981 sc0a[srcx] = sr1 + IMUL6(sr0) + SKIPSMImagePixelType(sa(sy)); \
982 sc0b[srcx] = (sr0 + SKIPSMImagePixelType(sa(sy))) * 4; \
983 out00 += sc0a[srcx]; \
984 out01 += sc0a[srcx]; \
985 out00 /= SKIPSMImagePixelType(SCALE_OUT00); \
986 out01 /= SKIPSMImagePixelType(SCALE_OUT01); \
987 da.set(cf(da(dx), out00), dx); \
988 da.set(cf(da(dxx), out01), dxx); \
989 if ((dst_w & 1) == 0) { \
992 out10 += sc0b[srcx]; \
993 out11 += sc0b[srcx]; \
994 out10 /= SKIPSMImagePixelType(SCALE_OUT10); \
995 out11 /= SKIPSMImagePixelType(SCALE_OUT11); \
996 da.set(cf(da(dx), out10), dx); \
997 da.set(cf(da(dxx), out11), dxx); \
1002 #define SKIPSM_EXPAND_ROW_COLUMN_END(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11) \
1003 out00 = sc1a[srcx] + IMUL6(sc0a[srcx]); \
1004 out00 /= SKIPSMImagePixelType(SCALE_OUT00); \
1005 da.set(cf(da(dx), out00), dx); \
1006 if ((dst_w & 1) == 0) { \
1007 out10 = sc1b[srcx] + IMUL6(sc0b[srcx]); \
1008 out10 /= SKIPSMImagePixelType(SCALE_OUT10); \
1010 da.set(cf(da(dx), out10), dx); \
1012 if ((dst_h & 1) == 0) { \
1013 out01 = sc0a[srcx]; \
1014 out01 /= SKIPSMImagePixelType(SCALE_OUT01); \
1015 da.set(cf(da(dxx), out01), dxx); \
1016 if ((dst_w & 1) == 0) { \
1017 out11 = sc0b[srcx]; \
1018 out11 /= SKIPSMImagePixelType(SCALE_OUT11); \
1020 da.set(cf(da(dxx), out11), dxx); \
1085 template <
typename SKIPSMImagePixelType,
1086 typename SrcImageIterator,
typename SrcAccessor,
1087 typename DestImageIterator,
typename DestAccessor,
1088 typename CombineFunctor>
1090 SrcImageIterator src_upperleft,
1091 SrcImageIterator src_lowerright,
1093 DestImageIterator dest_upperleft,
1094 DestImageIterator dest_lowerright,
1096 CombineFunctor cf) {
1098 int src_w = src_lowerright.x - src_upperleft.x;
1099 int src_h = src_lowerright.y - src_upperleft.y;
1100 int dst_w = dest_lowerright.x - dest_upperleft.x;
1101 int dst_h = dest_lowerright.y - dest_upperleft.y;
1104 SKIPSMImagePixelType current;
1105 SKIPSMImagePixelType out00, out10, out01, out11;
1106 SKIPSMImagePixelType sr0, sr1;
1107 SKIPSMImagePixelType *sc0a =
new SKIPSMImagePixelType[src_w + 1];
1108 SKIPSMImagePixelType *sc0b =
new SKIPSMImagePixelType[src_w + 1];
1109 SKIPSMImagePixelType *sc1a =
new SKIPSMImagePixelType[src_w + 1];
1110 SKIPSMImagePixelType *sc1b =
new SKIPSMImagePixelType[src_w + 1];
1113 const SKIPSMImagePixelType SKIPSMImageZero(NumericTraits<SKIPSMImagePixelType>::zero());
1115 DestImageIterator dy = dest_upperleft;
1116 DestImageIterator dyy = dest_upperleft;
1117 DestImageIterator dx = dy;
1118 DestImageIterator dxx = dyy;
1119 SrcImageIterator sy = src_upperleft;
1120 SrcImageIterator sx = sy;
1130 sr0 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1,0)));
1131 sr1 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-2,0)));
1133 sr0 = SKIPSMImageZero;
1134 sr1 = SKIPSMImageZero;
1137 for (sx = sy, srcx = 0; srcx < src_w; ++srcx, ++sx.x) {
1138 current = SKIPSMImagePixelType(sa(sx));
1139 sc0a[srcx] = sr1 +
IMUL6(sr0) + current;
1140 sc0b[srcx] = (sr0 + current) * 4;
1141 sc1a[srcx] = SKIPSMImageZero;
1142 sc1b[srcx] = SKIPSMImageZero;
1149 sc0a[srcx] = sr1 +
IMUL6(sr0) + SKIPSMImagePixelType(sa(sy));
1150 sc0b[srcx] = (sr0 + SKIPSMImagePixelType(sa(sy))) * 4;
1152 sc0a[srcx] = sr1 +
IMUL6(sr0);
1153 sc0b[srcx] = sr0 * 4;
1155 sc1a[srcx] = SKIPSMImageZero;
1156 sc1b[srcx] = SKIPSMImageZero;
1171 sr0 = SKIPSMImagePixelType(sa(sx));
1173 sr1 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1,0)));
1175 sr1 = SKIPSMImageZero;
1193 for (srcx = 2, ++sx.x; srcx < src_w; ++srcx, ++sx.x) {
1214 sr0 = SKIPSMImageZero;
1215 sr1 = SKIPSMImageZero;
1230 for (srcx = 2; srcx < src_w; ++srcx) {
1264 for (srcy = 2, sx = sy; srcy < src_h; ++srcy, ++sy.y, dy.y += 2, dyy.y += 2) {
1268 sr0 = SKIPSMImagePixelType(sa(sx));
1270 sr1 = SKIPSMImagePixelType(sa(sy, Diff2D(src_w-1,0)));
1272 sr1 = SKIPSMImageZero;
1290 for (srcx = 2, ++sx.x; srcx < src_w; ++srcx, ++sx.x) {
1313 sr0 = SKIPSMImageZero;
1314 sr1 = SKIPSMImageZero;
1329 for (srcx = 2; srcx < src_w; ++srcx) {
1357 template<
typename T1,
typename T2,
typename T3>
1360 return NumericTraits<T3>::fromPromote(a + b);
1365 template <
typename SKIPSMImagePixelType,
1366 typename SrcImageIterator,
typename SrcAccessor,
1367 typename DestImageIterator,
typename DestAccessor>
1368 inline void expand(
bool add,
bool wraparound,
1369 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1370 triple<DestImageIterator, DestImageIterator, DestAccessor> dest) {
1372 typedef typename DestAccessor::value_type DestPixelType;
1375 expand<SKIPSMImagePixelType>(add, wraparound,
1376 src.first, src.second, src.third,
1377 dest.first, dest.second, dest.third,
1381 expand<SKIPSMImagePixelType>(add, wraparound,
1382 src.first, src.second, src.third,
1383 dest.first, dest.second, dest.third,
1384 std::minus<SKIPSMImagePixelType>());
1393 template <
typename SrcImageType,
typename AlphaImageType,
typename PyramidImageType,
1394 int PyramidIntegerBits,
int PyramidFractionBits,
1395 typename SKIPSMImagePixelType,
typename SKIPSMAlphaPixelType>
1396 vector<PyramidImageType*> *gaussianPyramid(
unsigned int numLevels,
1398 typename SrcImageType::const_traverser src_upperleft,
1399 typename SrcImageType::const_traverser src_lowerright,
1400 typename SrcImageType::ConstAccessor sa,
1401 typename AlphaImageType::const_traverser alpha_upperleft,
1402 typename AlphaImageType::ConstAccessor aa) {
1404 vector<PyramidImageType*> *gp =
new vector<PyramidImageType*>();
1407 int w = src_lowerright.x - src_upperleft.x;
1408 int h = src_lowerright.y - src_upperleft.y;
1411 PyramidImageType *gp0 =
new PyramidImageType(w, h);
1414 copyToPyramidImage<SrcImageType, PyramidImageType, PyramidIntegerBits, PyramidFractionBits>(
1415 src_upperleft, src_lowerright, sa, gp0->upperLeft(), gp0->accessor());
1419 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1420 cout <<
"Generating Gaussian pyramid: g0";
1424 PyramidImageType *lastGP = gp0;
1425 AlphaImageType *lastA = NULL;
1426 for (
unsigned int l = 1; l < numLevels; l++) {
1428 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1438 PyramidImageType *gpn =
new PyramidImageType(w, h);
1439 AlphaImageType *nextA =
new AlphaImageType(w, h);
1441 if (lastA == NULL) {
1442 reduce<SKIPSMImagePixelType, SKIPSMAlphaPixelType>(wraparound,
1446 reduce<SKIPSMImagePixelType, SKIPSMAlphaPixelType>(wraparound,
1459 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1468 template <
typename SrcImageType,
typename AlphaImageType,
typename PyramidImageType,
1469 int PyramidIntegerBits,
int PyramidFractionBits,
1470 typename SKIPSMImagePixelType,
typename SKIPSMAlphaPixelType>
1471 inline vector<PyramidImageType*> *gaussianPyramid(
unsigned int numLevels,
1473 triple<typename SrcImageType::const_traverser, typename SrcImageType::const_traverser, typename SrcImageType::ConstAccessor> src,
1474 pair<typename AlphaImageType::const_traverser, typename AlphaImageType::ConstAccessor> alpha) {
1475 return gaussianPyramid<SrcImageType, AlphaImageType, PyramidImageType,
1476 PyramidIntegerBits, PyramidFractionBits,
1477 SKIPSMImagePixelType, SKIPSMAlphaPixelType>(
1478 numLevels, wraparound,
1479 src.first, src.second, src.third,
1480 alpha.first, alpha.second);
1484 template <
typename SrcImageType,
typename PyramidImageType,
1485 int PyramidIntegerBits,
int PyramidFractionBits,
1486 typename SKIPSMImagePixelType>
1487 vector<PyramidImageType*> *gaussianPyramid(
unsigned int numLevels,
1489 typename SrcImageType::const_traverser src_upperleft,
1490 typename SrcImageType::const_traverser src_lowerright,
1491 typename SrcImageType::ConstAccessor sa) {
1493 vector<PyramidImageType*> *gp =
new vector<PyramidImageType*>();
1496 int w = src_lowerright.x - src_upperleft.x;
1497 int h = src_lowerright.y - src_upperleft.y;
1500 PyramidImageType *gp0 =
new PyramidImageType(w, h);
1503 copyToPyramidImage<SrcImageType, PyramidImageType, PyramidIntegerBits, PyramidFractionBits>(
1504 src_upperleft, src_lowerright, sa,
1505 gp0->upperLeft(), gp0->accessor());
1509 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1510 cout <<
"Generating Gaussian pyramid: g0";
1514 PyramidImageType *lastGP = gp0;
1515 for (
unsigned int l = 1; l < numLevels; l++) {
1517 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1527 PyramidImageType *gpn =
new PyramidImageType(w, h);
1535 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1543 template <
typename SrcImageType,
typename PyramidImageType,
1544 int PyramidIntegerBits,
int PyramidFractionBits,
1545 typename SKIPSMImagePixelType>
1546 inline vector<PyramidImageType*> *gaussianPyramid(
unsigned int numLevels,
1548 triple<typename SrcImageType::const_traverser, typename SrcImageType::const_traverser, typename SrcImageType::ConstAccessor> src) {
1549 return gaussianPyramid<SrcImageType, PyramidImageType,
1550 PyramidIntegerBits, PyramidFractionBits, SKIPSMImagePixelType>(
1553 src.first, src.second, src.third);
1557 template <
typename SrcImageType,
typename AlphaImageType,
typename PyramidImageType,
1558 int PyramidIntegerBits,
int PyramidFractionBits,
1559 typename SKIPSMImagePixelType,
typename SKIPSMAlphaPixelType>
1560 vector<PyramidImageType*> *laplacianPyramid(
const char* exportName,
unsigned int numLevels,
1562 typename SrcImageType::const_traverser src_upperleft,
1563 typename SrcImageType::const_traverser src_lowerright,
1564 typename SrcImageType::ConstAccessor sa,
1565 typename AlphaImageType::const_traverser alpha_upperleft,
1566 typename AlphaImageType::ConstAccessor aa) {
1569 vector <PyramidImageType*> *gp =
1570 gaussianPyramid<SrcImageType, AlphaImageType, PyramidImageType,
1571 PyramidIntegerBits, PyramidFractionBits,
1572 SKIPSMImagePixelType, SKIPSMAlphaPixelType>(
1573 numLevels, wraparound,
1574 src_upperleft, src_lowerright, sa,
1575 alpha_upperleft, aa);
1579 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1580 cout <<
"Generating Laplacian pyramid:";
1586 for (
unsigned int l = 0; l < (numLevels-1); l++) {
1588 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1593 expand<SKIPSMImagePixelType>(
false, wraparound,
1598 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1599 cout <<
" l" << (numLevels-1) << endl;
1608 template <
typename SrcImageType,
typename AlphaImageType,
typename PyramidImageType,
1609 int PyramidIntegerBits,
int PyramidFractionBits,
1610 typename SKIPSMImagePixelType,
typename SKIPSMAlphaPixelType>
1611 inline vector<PyramidImageType*> *laplacianPyramid(
const char* exportName,
unsigned int numLevels,
1613 triple<typename SrcImageType::const_traverser, typename SrcImageType::const_traverser, typename SrcImageType::ConstAccessor> src,
1614 pair<typename AlphaImageType::const_traverser, typename AlphaImageType::ConstAccessor> alpha) {
1615 return laplacianPyramid<SrcImageType, AlphaImageType, PyramidImageType,
1616 PyramidIntegerBits, PyramidFractionBits,
1617 SKIPSMImagePixelType, SKIPSMAlphaPixelType>(
1619 numLevels, wraparound,
1620 src.first, src.second, src.third,
1621 alpha.first, alpha.second);
1625 template <
typename SKIPSMImagePixelType,
typename Pyram
idImageType>
1626 void collapsePyramid(
bool wraparound, vector<PyramidImageType*> *p) {
1628 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1629 cout <<
"Collapsing Laplacian pyramid: "
1630 <<
"l" << p->size()-1;
1636 for (
int l = (p->size()-2); l >= 0; l--) {
1638 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1643 expand<SKIPSMImagePixelType>(
true, wraparound,
1648 if (Verbose > VERBOSE_PYRAMID_MESSAGES) {
1655 template <
typename Pyram
idImageType>
1656 void exportPyramid(vector<PyramidImageType*> *v,
const char *prefix, VigraTrueType) {
1657 typedef typename PyramidImageType::value_type PyramidValueType;
1665 for (
unsigned int i = 0; i < v->size(); i++) {
1666 char filenameBuf[512];
1667 snprintf(filenameBuf, 512,
"%s%04u.tif", prefix, i);
1670 UInt16Image usPyramid((*v)[i]->width(), (*v)[i]->height());
1677 ImageExportInfo
info(filenameBuf);
1683 template <
typename Pyram
idImageType>
1684 void exportPyramid(vector<PyramidImageType*> *v,
const char *prefix, VigraFalseType) {
1685 typedef typename PyramidImageType::value_type PyramidVectorType;
1686 typedef typename PyramidVectorType::value_type PyramidValueType;
1688 for (
unsigned int i = 0; i < (v->size() - 1); i++) {
1690 initImage(
destImageRange(*((*v)[i])), NumericTraits<PyramidValueType>::zero());
1692 collapsePyramid(
false, v);
1694 for (
unsigned int i = 0; i < v->size(); i++) {
1695 char filenameBuf[512];
1696 snprintf(filenameBuf, 512,
"%s%04u.tif", prefix, i);
1699 UInt16RGBImage usPyramid((*v)[i]->width(), (*v)[i]->height());
1706 ImageExportInfo
info(filenameBuf);
1712 template <
typename Pyram
idImageType>
1713 void exportPyramid(vector<PyramidImageType*> *v,
const char *prefix) {
1714 typedef typename NumericTraits<typename PyramidImageType::value_type>::isScalar pyramid_is_scalar;
1715 exportPyramid(v, prefix, pyramid_is_scalar());
#define SKIPSM_EXPAND_SHIFT
void transformImage(vigra::triple< SrcImageIterator, SrcImageIterator, SrcAccessor > src, vigra::triple< DestImageIterator, DestImageIterator, DestAccessor > dest, std::pair< AlphaImageIterator, AlphaAccessor > alpha, vigra::Diff2D destUL, TRANSFORM &transform, PixelTransform &pixelTransform, bool warparound, Interpolator interpol, AppBase::ProgressDisplay *progress, bool singleThreaded=false)
Transform an image into the panorama.
#define SKIPSM_EXPAND_COLUMN_END(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11)
vigra::pair< typename ROIImage< Image, Alpha >::mask_traverser, typename ROIImage< Image, Alpha >::MaskAccessor > destMask(ROIImage< Image, Alpha > &img)
void reduce(bool wraparound, SrcImageIterator src_upperleft, SrcImageIterator src_lowerright, SrcAccessor sa, AlphaIterator alpha_upperleft, AlphaAccessor aa, DestImageIterator dest_upperleft, DestImageIterator dest_lowerright, DestAccessor da, DestAlphaIterator dest_alpha_upperleft, DestAlphaIterator dest_alpha_lowerright, DestAlphaAccessor daa)
The Burt & Adelson Reduce operation.
#define SKIPSM_EXPAND_ROW_COLUMN_END(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11)
#define SKIPSM_EXPAND_COLUMN_END_WRAPAROUND(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11)
#define SKIPSM_EXPAND_ROW_END(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11)
#define SKIPSM_EXPAND(SCALE_OUT00, SCALE_OUT10, SCALE_OUT01, SCALE_OUT11)
vigra::pair< typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::ImageAccessor > destImage(ROIImage< Image, Alpha > &img)
vigra::triple< typename ROIImage< Image, Mask >::image_const_traverser, typename ROIImage< Image, Mask >::image_const_traverser, typename ROIImage< Image, Mask >::ImageConstAccessor > srcImageRange(const ROIImage< Image, Mask > &img)
helper function for ROIImages
vigra::triple< typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::ImageAccessor > destImageRange(ROIImage< Image, Alpha > &img)
static void info(const char *fmt,...)
void expand(bool add, bool wraparound, SrcImageIterator src_upperleft, SrcImageIterator src_lowerright, SrcAccessor sa, DestImageIterator dest_upperleft, DestImageIterator dest_lowerright, DestAccessor da, CombineFunctor cf)
The Burt & Adelson Expand operation.
unsigned int filterHalfWidth(const unsigned int levels)
Calculate the half-width of a n-level filter.