5 //---------------------------------------------------------------------------
7 /// namespace for boundary conditions
8 // TODO: could actually be in a separate file... included from here.
13 std::size_t operator()(int i, std::size_t n)
16 if (i>=0) { assert(std::size_t(i) <= n); }
17 if (i == -1) return 3;
18 else if (i == int(n)) return n-4;
19 //else if (i == int(n)) return ( n%2 ? n-2 : n-4 );
26 std::size_t operator()(int i, std::size_t n)
29 if (i>=0) { assert(std::size_t(i) <= n); }
30 if (i == -1) return (n-1);
31 else if (i == int(n)) return 0;
39 //---------------------------------------------------------------------------
41 //----------------------------------------------//
42 // Various anonymous helper functions for DWT //
43 //----------------------------------------------//
47 /// Computes the number of elements reached when hopping in an array of size length with a step of size s.
48 struct SteppedSize : public std::binary_function<std::size_t, std::size_t, std::size_t>
50 std::size_t operator()(std::size_t length, std::size_t step)
51 { return 1 + (length-1) / step; }
54 /// Initial step size -- it should be the largest power of two x so that 1+(4-1)*x is smaller than n.
55 struct StepInit : public std::unary_function<std::size_t, std::size_t>
57 std::size_t operator()(std::size_t n)
59 // This loop may actually be more efficient than a closed form formula for reasonable n.
61 for (res = 1; 1+(4-1)*res <= n; res *= 2);
66 inline std::size_t index_offset(
67 std::size_t x, std::size_t y, std::size_t z,
68 std::size_t dimx, std::size_t dimy)
70 return x + dimx * ( y + dimy * z);
73 inline std::size_t index_offset(
74 std::size_t x, std::size_t y,
82 inline std::size_t wrap(int i, std::size_t n)
85 if (i>=0) { assert(std::size_t(i) <= n); }
86 if (i == -1) return 3;
87 else if (i == int(n)) return n-4;
88 //else if (i == int(n)) return ( n%2 ? n-2 : n-4 );
94 //---------------------------------------------------------------------------
96 template < typename T, typename BC >
97 void DWTCubic<T, BC>::Direct::
98 operator()(T * s, std::size_t n, std::size_t step)
100 for (std::size_t i = 0; i < n; i+=2) s[i*step] -= ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] ) * T(1.0/4.0);
101 for (std::size_t i = 1; i < n; i+=2) s[i*step] -= ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] );
102 for (std::size_t i = 0; i < n; i+=2) s[i*step] += ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] ) * T(3.0/16.0);
104 // This is the "constant L1 norm" normalization
105 //for (std::size_t i = 0; i < n; i+=2) s[i*step] *= T(4);
107 // This is the usual normalization
108 for (std::size_t i = 0; i < n; i+=2) s[i*step] *= T(2);
109 for (std::size_t i = 1; i < n; i+=2) s[i*step] *= T(1.0/2.0);
111 //---------------------------------------------------------------------------
113 template < typename T, typename BC >
114 void DWTCubic<T,BC>::Inverse::
115 operator()(T * s, std::size_t n, std::size_t step)
117 // This is the "constant L1 norm" normalization
118 //for (std::size_t i = 0; i < n; i+=2) s[i*step] *= T(1.0/4.0);
120 std::size_t step2 = 2 * step;
121 std::size_t ns = n * step;
122 // This is the usual normalization
123 for (std::size_t i = 0; i < ns; i += step2) s[i] *= T(1.0/2.0);
124 for (std::size_t i = step; i < ns; i += step2) s[i] *= T(2);
126 for (std::size_t i = 0, j = 0; i < n; i+=2, j += step2) s[j] -= ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] ) * T(3.0/16.0);
127 for (std::size_t i = 1, j = step; i < n; i+=2, j += step2) s[j] += ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] );
128 for (std::size_t i = 0, j = 0; i < n; i+=2, j += step2) s[j] += ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] ) * T(1.0/4.0);
131 //---------------------------------------------------------------------------
133 template < typename T, typename BC >
134 void DWTCubicConjugate< T, BC >::Direct::
135 operator()(T * s, std::size_t n, std::size_t step)
137 for (std::size_t i = 1; i < n; i+=2) s[i*step] += ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] ) * T(1.0/4.0);
138 for (std::size_t i = 0; i < n; i+=2) s[i*step] += ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] );
139 for (std::size_t i = 1; i < n; i+=2) s[i*step] -= ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] ) * T(3.0/16.0);
141 //for (std::size_t i = 0; i < n; i+=2) s[i*step] *= T(1.0/(4.0));
143 //for (std::size_t i = 0; i < n; i+=2) s[i*step] *= T(std::sqrt(2.0)/4);
144 //for (std::size_t i = 1; i < n; i+=2) s[i*step] *= T(std::sqrt(2.0));
146 for (std::size_t i = 0; i < n; i+=2) s[i*step] *= T(1.0/2.0);
147 for (std::size_t i = 1; i < n; i+=2) s[i*step] *= T(2);
150 //---------------------------------------------------------------------------
152 template < typename T, typename BC >
153 void DWTCubicConjugate< T, BC >::Inverse::
154 operator()(T * s, std::size_t n, std::size_t step)
156 for (std::size_t i = 0; i < n; i+=2) s[i*step] *= T(2);
157 for (std::size_t i = 1; i < n; i+=2) s[i*step] *= T(1.0/2.0);
158 for (std::size_t i = 1; i < n; i+=2) s[i*step] += ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] ) * T(3.0/16.0);
159 for (std::size_t i = 0; i < n; i+=2) s[i*step] -= ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] );
160 for (std::size_t i = 1; i < n; i+=2) s[i*step] -= ( s[BC()(i+1,n)*step] + s[BC()(i-1,n)*step] ) * T(1.0/4.0);
163 //---------------------------------------------------------------------------
165 template < typename DWT >
166 void DWTND<DWT, 2>::operator()
168 typename DWT::prec_type * im
169 , numeric_array<std::size_t,2> dim
171 , numeric_array<std::size_t,2> step
172 , numeric_array<std::size_t,2> real_dim
179 for (std::size_t j = 0; j < dim[1]; ++j)
181 m_dwt(im+j*step[1]*real_dim[0], dim[0], step[0]);
186 for (std::size_t i = 0; i < dim[0]; ++i)
188 m_dwt(im+i*step[0], dim[1], real_dim[0]*step[1]);
193 //---------------------------------------------------------------------------
195 template < typename DWT >
196 void DWTND<DWT,3>::operator()(prec_type * im, numeric_array<std::size_t,3> dim, int dir, numeric_array<std::size_t,3> step, numeric_array<std::size_t,3> real_dim)
202 for (std::size_t k = 0; k < dim[2]; ++k)
203 for (std::size_t j = 0; j < dim[1]; ++j)
205 m_dwt(im+index_offset(0, j*step[1], k*step[2], real_dim[0], real_dim[1]), dim[0], step[0]);
210 for (std::size_t k = 0; k < dim[2]; ++k)
211 for (std::size_t i = 0; i < dim[0]; ++i)
213 m_dwt(im+index_offset(i*step[0], 0, k*step[2], real_dim[0], real_dim[1]), dim[1], real_dim[0]*step[1]);
218 for (std::size_t j = 0; j < dim[1]; ++j)
219 for (std::size_t i = 0; i < dim[0]; ++i)
221 m_dwt(im+index_offset(i*step[0], j*step[1], 0, real_dim[0], real_dim[1]), dim[2], real_dim[0]*real_dim[1]*step[2]);
226 //---------------------------------------------------------------------------
228 template < typename DWT >
229 void MultiDWT::Direct<DWT>::operator()
231 prec_type * s, ///< Input array. This also contains the output, the computation being in-place.
232 std::size_t n ///< Length of the input.
239 n = SteppedSize()(n, 2);
244 //---------------------------------------------------------------------------
246 template < typename IDWT >
247 void MultiDWT::Inverse<IDWT>::operator()
249 prec_type * s, ///< The input array. Also contains the output, since the computation is in-place.
250 std::size_t n ///< The input length.
253 for (std::size_t sm = StepInit()(n); sm > 0; sm /= 2)
255 m_idwt(s, SteppedSize()(n, sm), sm);
259 //---------------------------------------------------------------------------
261 template < typename DWT, std::size_t D >
262 void MultiDWTND::Direct<DWT,D>::operator()
265 , const numeric_array<std::size_t,D> dim
268 numeric_array<std::size_t,D> sm;
269 std::fill(sm.begin(), sm.end(), 1);
270 numeric_array<std::size_t,D> sm_bis = sm;
271 numeric_array<std::size_t,D> ssize = dim;
272 numeric_array<std::size_t,D> ssize_bis = dim;
274 numeric_array<bool,D> f;
275 for (std::size_t i = 0; i < D; ++i) f[i] = (dim[i] >= 4);
277 while (or_all(f.begin(), f.end()))
279 // loop in all dimensions, starting from the lowest.
280 for (std::size_t i = 0; i < D; ++i)
284 //std::cout << "dir " << i << " step " << sm[0] << " " << sm[1] << " " << sm[2] << " dim " << ssize[0] << " " << ssize[1] << " " << ssize[2] << std::endl;
285 m_dwt(im, ssize, i, sm, dim);
286 if (ssize_bis[i] >= (4-1)*2+1)
288 //ssize_bis[i] = 1 + (ssize_bis[i]-1) / 2;
289 ssize_bis[i] = SteppedSize()(ssize_bis[i], 2);
300 //---------------------------------------------------------------------------
302 template < typename IDWT, std::size_t D >
304 MultiDWTND::Inverse< IDWT, D >::
305 operator()(prec_type * im, const numeric_array<std::size_t,D> dim)
307 numeric_array<std::size_t,D> sm;
308 numeric_array<std::size_t,D> ssize;
309 std::transform(dim.begin(), dim.end(), sm.begin(), StepInit());
310 // NB: sm_bis has to be initialized here, because it might not be initialized in the
312 numeric_array<std::size_t,D> sm_bis = sm;
316 // loop in all dimensions, starting with the highest
317 for (int i = D-1; i >= 0; --i)
319 if (sm[i] >= 1 && greater_than_others(sm, i))
321 std::transform(dim.begin(), dim.end(), sm.begin(), ssize.begin(), SteppedSize());
322 //for (std::size_t j = 0; j < D; ++j) ssize[j] = SteppedSize()(dim[j], sm[j]);
323 //std::cout << "dir " << i << " step " << sm[0] << " " << sm[1] << " dim " << ssize[0] << " " << ssize[1] << " rdim " << dim[0] << " " << dim[1] << std::endl;
324 //std::cout << "dir " << i << " step " << sm[0] << " " << sm[1] << " " << sm[2] << " dim " << ssize[0] << " " << ssize[1] << " " << ssize[2] << std::endl;
325 m_idwt(im, ssize, i, sm, dim);
326 sm_bis[i] = sm[i] / 2;
330 } while (!all_less(sm, 1));
333 //---------------------------------------------------------------------------
335 template < typename T >
336 void MultiDWTPower< T >::
337 operator()(T * s, std::size_t length, std::size_t hop)
339 for (std::size_t step = 2; SteppedSize()(length, step) >= 4; step *= 2)
341 for (std::size_t i = 0; i < length; i += step)
348 //---------------------------------------------------------------------------
350 template < typename T >
351 void MultiDWTNDPower<T,3>::
352 operator()(T * im, numeric_array<std::size_t,3> dim)
354 MultiDWTPower<T> power;
355 for (std::size_t k = 0; k < dim[2]; ++k)
356 for (std::size_t j = 0; j < dim[1]; ++j)
358 power(im+(k*dim[1]+j)*dim[0], dim[0], 1);
360 for (std::size_t k = 0; k < dim[2]; ++k)
361 for (std::size_t i = 0; i < dim[0]; ++i)
363 power(im+i+ k*dim[0]*dim[1], dim[1], dim[0]);
365 for (std::size_t j = 0; j < dim[1]; ++j)
366 for (std::size_t i = 0; i < dim[0]; ++i)
368 power(im+i+j*dim[0], dim[2], dim[0]*dim[1]);
372 //---------------------------------------------------------------------------
376 /// Reorder sequence by putting the even index first, then the odd.
377 template < typename T >
378 void shuffle_step(T * s, std::size_t length, T * buffer, std::size_t jump)
380 for (std::size_t i = 0; i < length; ++i)
382 buffer[i] = s[i*jump];
384 for (T * p = buffer; p < buffer+length; p+=2, s+=jump)
388 for (T * p = buffer+1; p < buffer+length; p+=2, s+=jump)
394 /// Reorder sequence by dispatching first half on even indexes, second half on odd.
395 template < typename T >
396 void unshuffle_step(T * s, std::size_t length, T * buffer, std::size_t jump)
398 for (std::size_t i = 0, j = 0; i < length; ++i, j += jump)
402 std::size_t jump2 = 2*jump;
403 std::size_t jl = jump * length;
404 for (T * p = s; p < s + jl; p += jump2)
408 for (T * p = s+jump; p < s + jl; p += jump2)
415 //---------------------------------------------------------------------------
417 template < typename T >
418 void MultiDWTShuffle<T>::shuffle(T * s, std::size_t length, std::size_t jump)
420 if (m_buffer.size() < length) m_buffer.resize(length);
421 for (std::size_t step = 1; SteppedSize()(length, step) >= 4; step *= 2)
423 //std::cout << "length " << SteppedSize()(step, length) << " step " << step << std::endl;
424 shuffle_step(s, SteppedSize()(length, step), &m_buffer.front(), jump);
428 //---------------------------------------------------------------------------
430 template < typename T >
431 void MultiDWTShuffle<T>::unshuffle(T * s, std::size_t length, std::size_t jump)
433 if (m_buffer.size() < length) m_buffer.resize(length);
435 for (step = 1; SteppedSize()(length, step) >= 4; step *= 2);
439 unshuffle_step(s, SteppedSize()(length, step), &m_buffer.front(), jump);
443 //---------------------------------------------------------------------------
445 template < typename T >
446 void MultiDWTNDShuffle<T,2>::shuffle(T * im, numeric_array<std::size_t,2> dim)
448 MultiDWTShuffle<T> shuffle;
449 for (std::size_t j = 0; j < dim[1]; ++j)
451 shuffle.shuffle(im+j*dim[0], dim[0], 1);
453 for (std::size_t i = 0; i < dim[0]; ++i)
455 shuffle.shuffle(im+i, dim[1], dim[0]);
459 //---------------------------------------------------------------------------
461 template < typename T >
462 void MultiDWTNDShuffle<T,2>::unshuffle(T * im, numeric_array<std::size_t,2> dim)
464 MultiDWTShuffle<T> shuffle;
465 for (std::size_t j = 0; j < dim[1]; ++j)
467 shuffle.unshuffle(im+j*dim[0], dim[0], 1);
469 for (std::size_t i = 0; i < dim[0]; ++i)
471 shuffle.unshuffle(im+i, dim[1], dim[0]);
475 //---------------------------------------------------------------------------
477 template < typename T >
478 void MultiDWTNDShuffle<T,3>::shuffle(T * im, numeric_array<std::size_t,3> dim)
480 MultiDWTShuffle<T> shuffle;
481 for (std::size_t k = 0; k < dim[2]; ++k)
482 for (std::size_t j = 0; j < dim[1]; ++j)
484 shuffle.shuffle(im+(k*dim[1]+j)*dim[0], dim[0], 1);
486 for (std::size_t k = 0; k < dim[2]; ++k)
487 for (std::size_t i = 0; i < dim[0]; ++i)
489 shuffle.shuffle(im+i+ k*dim[0]*dim[1], dim[1], dim[0]);
491 for (std::size_t j = 0; j < dim[1]; ++j)
492 for (std::size_t i = 0; i < dim[0]; ++i)
494 shuffle.shuffle(im+i+j*dim[0], dim[2], dim[0]*dim[1]);
498 //---------------------------------------------------------------------------
500 template < typename T >
501 void MultiDWTNDShuffle<T,3>::unshuffle(T * im, numeric_array<std::size_t,3> dim)
503 MultiDWTShuffle<T> shuffle;
504 for (std::size_t k = 0; k < dim[2]; ++k)
505 for (std::size_t j = 0; j < dim[1]; ++j)
507 shuffle.unshuffle(im+(k*dim[1]+j)*dim[0], dim[0], 1);
509 for (std::size_t k = 0; k < dim[2]; ++k)
510 for (std::size_t i = 0; i < dim[0]; ++i)
512 shuffle.unshuffle(im+i+ k*dim[0]*dim[1], dim[1], dim[0]);
514 for (std::size_t j = 0; j < dim[1]; ++j)
515 for (std::size_t i = 0; i < dim[0]; ++i)
517 shuffle.unshuffle(im+i+j*dim[0], dim[2], dim[0]*dim[1]);
521 //---------------------------------------------------------------------------
523 /// In-place (4-2) DWT via lifting.
524 template < typename T, typename BC >
525 inline void dwt_cubic(T * s, std::size_t n, std::size_t step, BC)
527 typename DWTCubic<T, BC>::Direct()(s, n, step);
529 /// Same as above with cyclic bc by default
530 template < typename T >
531 inline void dwt_cubic(T * s, std::size_t n, std::size_t step)
533 dwt_cubic(s, n, step, bc::Cyclic());
536 /// In-place inverse (4-2) DWT via lifting.
537 template < typename T, typename BC >
538 inline void idwt_cubic(T * s, std::size_t n, std::size_t step, BC)
540 typename DWTCubic<T, BC>::Inverse()(s, n, step);
542 /// Same as above with cyclic bc by default
543 template < typename T >
544 inline void idwt_cubic(T * s, std::size_t n, std::size_t step)
546 idwt_cubic(s, n, step, bc::Cyclic());
549 template < typename T, std::size_t D >
550 inline void dwtND_cubic(T * im, numeric_array<std::size_t, D> dim)
552 numeric_array<std::size_t, D> step;
553 for (std::size_t i = 0; i < D; ++i)
558 for (std::size_t i = 0; i < D; ++i)
560 DWTND<typename DWTCubic<T, bc::Cyclic>::Direct,D>()(im, dim, i, step, dim);
564 template < typename T, std::size_t D >
565 inline void dwtND_cubicConjugate(T * im, numeric_array<std::size_t, D> dim)
567 numeric_array<std::size_t, D> step;
568 for (std::size_t i = 0; i < D; ++i)
573 for (std::size_t i = 0; i < D; ++i)
575 DWTND<typename DWTCubicConjugate<T, bc::Cyclic>::Direct,D>()(im, dim, i, step, dim);
581 /// In-place multidimensional (4-2) DWT along a given direction.
582 template < typename T, std::size_t D >
583 inline void dwtND_cubic(T * im, numeric_array<std::size_t,D> dim, int dir, numeric_array<std::size_t,D> step, numeric_array<std::size_t,D> real_dim)
585 DWTND<DWTCubic<T>,D>()(im, dim, dir, step, real_dim);
588 /// In-place multidimensional inverse (4-2) DWT along a given direction.
589 template < typename T, std::size_t D >
590 inline void idwtND_cubic(T * im, numeric_array<std::size_t,D> dim, int dir, numeric_array<std::size_t,D> step, numeric_array<std::size_t,D> real_dim)
592 DWTND<IDWTCubic<T>,D>()(im, dim, dir, step, real_dim);
596 template < typename T, typename BC >
597 void multi_dwt_cubic(T * s, std::size_t dim, BC)
599 MultiDWT::Direct<typename DWTCubic<T, BC>::Direct>()(s, dim);
601 template < typename T >
602 void multi_dwt_cubic(T * s, std::size_t dim)
604 multi_dwt_cubic(s, dim, bc::Cyclic());
607 template < typename T, typename BC >
608 void multi_idwt_cubic(T * s, std::size_t dim, BC)
610 MultiDWT::Inverse<typename DWTCubic<T, BC>::Inverse>()(s, dim);
612 template < typename T >
613 void multi_idwt_cubic(T * s, std::size_t dim)
615 multi_idwt_cubic(s, dim, bc::Cyclic());
618 template < typename T, std::size_t D, typename BC >
619 void multi_dwtND_cubic(T * im, numeric_array<std::size_t,D> dim, BC)
621 MultiDWTND::Direct<DWTND<typename DWTCubic<T, BC>::Direct,D>,D>()(im, dim);
623 template < typename T, std::size_t D >
624 void multi_dwtND_cubic(T * im, numeric_array<std::size_t,D> dim)
626 multi_dwtND_cubic(im, dim, bc::Cyclic());
629 template < typename T, std::size_t D, typename BC >
630 void multi_idwtND_cubic(T * im, numeric_array<std::size_t,D> dim, BC)
632 MultiDWTND::Inverse<DWTND<typename DWTCubic<T, BC>::Inverse,D>,D>()(im, dim);
634 template < typename T, std::size_t D >
635 void multi_idwtND_cubic(T * im, numeric_array<std::size_t,D> dim)
637 multi_idwtND_cubic(im, dim, bc::Cyclic());
640 template < typename T >
641 void multi_dwt_shuffle(T * s, std::size_t n)
643 MultiDWTShuffle<T>().shuffle(s,n);
646 template < typename T >
647 void multi_dwt_unshuffle(T * s, std::size_t n)
649 MultiDWTShuffle<T>().unshuffle(s, n);
652 template < typename T, std::size_t D >
653 void multi_dwtND_shuffle(T * im, numeric_array<std::size_t, D> dim)
655 MultiDWTNDShuffle<T,D>().shuffle(im, dim);
658 template < typename T, std::size_t D >
659 void multi_dwtND_unshuffle(T * im, numeric_array<std::size_t, D> dim)
661 MultiDWTNDShuffle<T,D>().unshuffle(im, dim);
664 template < typename T, std::size_t D >
665 void multi_dwtND_power(T * im, numeric_array<std::size_t, D> dim)
667 MultiDWTNDPower<T,D>()(im, dim);