VTK-m  2.0
WaveletDWT.h
Go to the documentation of this file.
1 //============================================================================
2 // Copyright (c) Kitware, Inc.
3 // All rights reserved.
4 // See LICENSE.txt for details.
5 //
6 // This software is distributed WITHOUT ANY WARRANTY; without even
7 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 // PURPOSE. See the above copyright notice for more information.
9 //============================================================================
10 
11 #ifndef vtk_m_worklet_wavelets_waveletdwt_h
12 #define vtk_m_worklet_wavelets_waveletdwt_h
13 
14 #include <vtkm/worklet/wavelets/WaveletBase.h>
15 
16 #include <vtkm/worklet/wavelets/WaveletTransforms.h>
17 
18 #include <vtkm/cont/ArrayHandleConcatenate.h>
19 #include <vtkm/cont/ArrayHandleCounting.h>
20 #include <vtkm/cont/ArrayHandlePermutation.h>
21 
22 #include <vtkm/Math.h>
23 #include <vtkm/cont/Timer.h>
24 
25 namespace vtkm
26 {
27 namespace worklet
28 {
29 namespace wavelets
30 {
31 
32 class WaveletDWT : public WaveletBase
33 {
34 public:
35  // Constructor
36  WaveletDWT(WaveletName name)
37  : WaveletBase(name)
38  {
39  }
40 
41  // Function: extend a cube in X direction
42  template <typename SigInArrayType, typename ExtensionArrayType>
43  vtkm::Id Extend3DLeftRight(const SigInArrayType& sigIn, // input
44  vtkm::Id sigDimX,
45  vtkm::Id sigDimY,
46  vtkm::Id sigDimZ,
47  vtkm::Id sigStartX,
48  vtkm::Id sigStartY,
49  vtkm::Id sigStartZ,
50  vtkm::Id sigPretendDimX,
51  vtkm::Id sigPretendDimY,
52  vtkm::Id sigPretendDimZ,
53  ExtensionArrayType& ext1, // output
54  ExtensionArrayType& ext2, // output
55  vtkm::Id addLen,
56  vtkm::worklet::wavelets::DWTMode ext1Method,
57  vtkm::worklet::wavelets::DWTMode ext2Method,
58  bool pretendSigPaddedZero,
59  bool padZeroAtExt2)
60  {
61  // pretendSigPaddedZero and padZeroAtExt2 cannot happen at the same time
62  VTKM_ASSERT(!pretendSigPaddedZero || !padZeroAtExt2);
63 
64  if (addLen == 0) // Haar kernel
65  {
66  ext1.Allocate(0); // No extension on the left side
67  if (pretendSigPaddedZero || padZeroAtExt2) // plane of size 1*dimY*dimZ
68  {
69  ext2.Allocate(sigPretendDimY * sigPretendDimZ);
70  WaveletBase::DeviceAssignZero3DPlaneX(ext2, 1, sigPretendDimY, sigPretendDimZ, 0);
71  }
72  else
73  {
74  ext2.Allocate(0); // No extension on the right side
75  }
76  return 0;
77  }
78 
79  using ValueType = typename SigInArrayType::ValueType;
80  using ExtendArrayType = vtkm::cont::ArrayHandle<ValueType>;
81  using ExtensionWorklet = vtkm::worklet::wavelets::ExtensionWorklet3D;
82  using DispatcherType = typename vtkm::worklet::DispatcherMapField<ExtensionWorklet>;
83  vtkm::Id extDimX, extDimY, extDimZ;
84  vtkm::worklet::wavelets::ExtensionDirection dir;
85 
86  { // First work on left extension
87  dir = LEFT;
88  extDimX = addLen;
89  extDimY = sigPretendDimY;
90  extDimZ = sigPretendDimZ;
91 
92  ext1.Allocate(extDimX * extDimY * extDimZ);
93  ExtensionWorklet worklet(extDimX,
94  extDimY,
95  extDimZ,
96  sigDimX,
97  sigDimY,
98  sigDimZ,
99  sigStartX,
100  sigStartY,
101  sigStartZ,
102  sigPretendDimX,
103  sigPretendDimY,
104  sigPretendDimZ,
105  ext1Method,
106  dir,
107  false); // not treating input signal as having zeros
108  DispatcherType dispatcher(worklet);
109  dispatcher.Invoke(ext1, sigIn);
110  }
111 
112  // Then work on right extension
113  dir = RIGHT;
114  extDimY = sigPretendDimY;
115  extDimZ = sigPretendDimZ;
116  if (!pretendSigPaddedZero && !padZeroAtExt2)
117  {
118  extDimX = addLen;
119  ext2.Allocate(extDimX * extDimY * extDimZ);
120  ExtensionWorklet worklet(extDimX,
121  extDimY,
122  extDimZ,
123  sigDimX,
124  sigDimY,
125  sigDimZ,
126  sigStartX,
127  sigStartY,
128  sigStartZ,
129  sigPretendDimX,
130  sigPretendDimY,
131  sigPretendDimZ,
132  ext2Method,
133  dir,
134  false);
135  DispatcherType dispatcher(worklet);
136  dispatcher.Invoke(ext2, sigIn);
137  }
138  else if (!pretendSigPaddedZero && padZeroAtExt2)
139  { // This case is not exactly padding a zero at the end of Ext2.
140  // Rather, it is to increase extension length by one and fill it
141  // to be whatever mirrored.
142  extDimX = addLen + 1;
143  ext2.Allocate(extDimX * extDimY * extDimZ);
144  ExtensionWorklet worklet(extDimX,
145  extDimY,
146  extDimZ,
147  sigDimX,
148  sigDimY,
149  sigDimZ,
150  sigStartX,
151  sigStartY,
152  sigStartZ,
153  sigPretendDimX,
154  sigPretendDimY,
155  sigPretendDimZ,
156  ext2Method,
157  dir,
158  false);
159  DispatcherType dispatcher(worklet);
160  dispatcher.Invoke(ext2, sigIn);
161  }
162  else // pretendSigPaddedZero
163  {
164  ExtendArrayType ext2Temp;
165  extDimX = addLen;
166  ext2Temp.Allocate(extDimX * extDimY * extDimZ);
167  ExtensionWorklet worklet(extDimX,
168  extDimY,
169  extDimZ,
170  sigDimX,
171  sigDimY,
172  sigDimZ,
173  sigStartX,
174  sigStartY,
175  sigStartZ,
176  sigPretendDimX,
177  sigPretendDimY,
178  sigPretendDimZ,
179  ext2Method,
180  dir,
181  true); // pretend sig is padded a zero
182  DispatcherType dispatcher(worklet);
183  dispatcher.Invoke(ext2Temp, sigIn);
184 
185  // Give ext2 one layer thicker to hold the pretend zeros from signal.
186  ext2.Allocate((extDimX + 1) * extDimY * extDimZ);
187  WaveletBase::DeviceCubeCopyTo(
188  ext2Temp, extDimX, extDimY, extDimZ, ext2, extDimX + 1, extDimY, extDimZ, 1, 0, 0);
189  WaveletBase::DeviceAssignZero3DPlaneX(ext2, extDimX + 1, extDimY, extDimZ, 0);
190  }
191  return 0;
192  }
193 
194  // Function: extend a cube in Y direction
195  template <typename SigInArrayType, typename ExtensionArrayType>
196  vtkm::Id Extend3DTopDown(const SigInArrayType& sigIn, // input
197  vtkm::Id sigDimX,
198  vtkm::Id sigDimY,
199  vtkm::Id sigDimZ,
200  vtkm::Id sigStartX,
201  vtkm::Id sigStartY,
202  vtkm::Id sigStartZ,
203  vtkm::Id sigPretendDimX,
204  vtkm::Id sigPretendDimY,
205  vtkm::Id sigPretendDimZ,
206  ExtensionArrayType& ext1, // output
207  ExtensionArrayType& ext2, // output
208  vtkm::Id addLen,
209  vtkm::worklet::wavelets::DWTMode ext1Method,
210  vtkm::worklet::wavelets::DWTMode ext2Method,
211  bool pretendSigPaddedZero,
212  bool padZeroAtExt2)
213  {
214  // pretendSigPaddedZero and padZeroAtExt2 cannot happen at the same time
215  VTKM_ASSERT(!pretendSigPaddedZero || !padZeroAtExt2);
216 
217  if (addLen == 0) // Haar kernel
218  {
219  ext1.Allocate(0); // No extension on the top side
220  if (pretendSigPaddedZero || padZeroAtExt2) // plane of size dimX*dimZ
221  {
222  ext2.Allocate(sigPretendDimX * 1 * sigPretendDimZ);
223  WaveletBase::DeviceAssignZero3DPlaneY(ext2, sigPretendDimX, 1, sigPretendDimZ, 0);
224  }
225  else
226  {
227  ext2.Allocate(0); // No extension on the right side
228  }
229  return 0;
230  }
231 
232  using ValueType = typename SigInArrayType::ValueType;
233  using ExtendArrayType = vtkm::cont::ArrayHandle<ValueType>;
234  using ExtensionWorklet = vtkm::worklet::wavelets::ExtensionWorklet3D;
235  using DispatcherType = typename vtkm::worklet::DispatcherMapField<ExtensionWorklet>;
236  vtkm::Id extDimX, extDimY, extDimZ;
237  vtkm::worklet::wavelets::ExtensionDirection dir;
238 
239  { // First work on top extension
240  dir = TOP;
241  extDimX = sigPretendDimX;
242  extDimY = addLen;
243  extDimZ = sigPretendDimZ;
244 
245  ext1.Allocate(extDimX * extDimY * extDimZ);
246  ExtensionWorklet worklet(extDimX,
247  extDimY,
248  extDimZ,
249  sigDimX,
250  sigDimY,
251  sigDimZ,
252  sigStartX,
253  sigStartY,
254  sigStartZ,
255  sigPretendDimX,
256  sigPretendDimY,
257  sigPretendDimZ,
258  ext1Method,
259  dir,
260  false); // not treating input signal as having zeros
261  DispatcherType dispatcher(worklet);
262  dispatcher.Invoke(ext1, sigIn);
263  }
264 
265  // Then work on bottom extension
266  dir = BOTTOM;
267  extDimX = sigPretendDimX;
268  extDimZ = sigPretendDimZ;
269  if (!pretendSigPaddedZero && !padZeroAtExt2)
270  {
271  extDimY = addLen;
272  ext2.Allocate(extDimX * extDimY * extDimZ);
273  ExtensionWorklet worklet(extDimX,
274  extDimY,
275  extDimZ,
276  sigDimX,
277  sigDimY,
278  sigDimZ,
279  sigStartX,
280  sigStartY,
281  sigStartZ,
282  sigPretendDimX,
283  sigPretendDimY,
284  sigPretendDimZ,
285  ext2Method,
286  dir,
287  false);
288  DispatcherType dispatcher(worklet);
289  dispatcher.Invoke(ext2, sigIn);
290  }
291  else if (!pretendSigPaddedZero && padZeroAtExt2)
292  { // This case is not exactly padding a zero at the end of Ext2.
293  // Rather, it is to increase extension length by one and fill it
294  // to be whatever mirrored.
295  extDimY = addLen + 1;
296  ext2.Allocate(extDimX * extDimY * extDimZ);
297  ExtensionWorklet worklet(extDimX,
298  extDimY,
299  extDimZ,
300  sigDimX,
301  sigDimY,
302  sigDimZ,
303  sigStartX,
304  sigStartY,
305  sigStartZ,
306  sigPretendDimX,
307  sigPretendDimY,
308  sigPretendDimZ,
309  ext2Method,
310  dir,
311  false);
312  DispatcherType dispatcher(worklet);
313  dispatcher.Invoke(ext2, sigIn);
314  }
315  else // pretendSigPaddedZero
316  {
317  ExtendArrayType ext2Temp;
318  extDimY = addLen;
319  ext2Temp.Allocate(extDimX * extDimY * extDimZ);
320  ExtensionWorklet worklet(extDimX,
321  extDimY,
322  extDimZ,
323  sigDimX,
324  sigDimY,
325  sigDimZ,
326  sigStartX,
327  sigStartY,
328  sigStartZ,
329  sigPretendDimX,
330  sigPretendDimY,
331  sigPretendDimZ,
332  ext2Method,
333  dir,
334  true); // pretend sig is padded a zero
335  DispatcherType dispatcher(worklet);
336  dispatcher.Invoke(ext2Temp, sigIn);
337 
338  // Give ext2 one layer thicker to hold the pretend zeros from signal.
339  ext2.Allocate(extDimX * (extDimY + 1) * extDimZ);
340  WaveletBase::DeviceCubeCopyTo(
341  ext2Temp, extDimX, extDimY, extDimZ, ext2, extDimX, extDimY + 1, extDimZ, 0, 1, 0);
342  WaveletBase::DeviceAssignZero3DPlaneY(ext2, extDimX, extDimY + 1, extDimZ, 0);
343  }
344  return 0;
345  }
346 
347  // Function: extend a cube in Z direction
348  template <typename SigInArrayType, typename ExtensionArrayType>
349  vtkm::Id Extend3DFrontBack(const SigInArrayType& sigIn, // input
350  vtkm::Id sigDimX,
351  vtkm::Id sigDimY,
352  vtkm::Id sigDimZ,
353  vtkm::Id sigStartX,
354  vtkm::Id sigStartY,
355  vtkm::Id sigStartZ,
356  vtkm::Id sigPretendDimX,
357  vtkm::Id sigPretendDimY,
358  vtkm::Id sigPretendDimZ,
359  ExtensionArrayType& ext1, // output
360  ExtensionArrayType& ext2, // output
361  vtkm::Id addLen,
362  vtkm::worklet::wavelets::DWTMode ext1Method,
363  vtkm::worklet::wavelets::DWTMode ext2Method,
364  bool pretendSigPaddedZero,
365  bool padZeroAtExt2)
366  {
367  // pretendSigPaddedZero and padZeroAtExt2 cannot happen at the same time
368  VTKM_ASSERT(!pretendSigPaddedZero || !padZeroAtExt2);
369 
370  if (addLen == 0) // Haar kernel
371  {
372  ext1.Allocate(0); // No extension on the front side
373  if (pretendSigPaddedZero || padZeroAtExt2) // plane of size dimX * dimY
374  {
375  ext2.Allocate(sigPretendDimX * sigPretendDimY * 1);
376  WaveletBase::DeviceAssignZero3DPlaneZ(ext2, sigPretendDimX, sigPretendDimY, 1, 0);
377  }
378  else
379  {
380  ext2.Allocate(0); // No extension on the right side
381  }
382  return 0;
383  }
384 
385  using ValueType = typename SigInArrayType::ValueType;
386  using ExtendArrayType = vtkm::cont::ArrayHandle<ValueType>;
387  using ExtensionWorklet = vtkm::worklet::wavelets::ExtensionWorklet3D;
388  using DispatcherType = typename vtkm::worklet::DispatcherMapField<ExtensionWorklet>;
389  vtkm::Id extDimX, extDimY, extDimZ;
390  vtkm::worklet::wavelets::ExtensionDirection dir;
391 
392  { // First work on front extension
393  dir = FRONT;
394  extDimX = sigPretendDimX;
395  extDimY = sigPretendDimY;
396  extDimZ = addLen;
397 
398  ext1.Allocate(extDimX * extDimY * extDimZ);
399  ExtensionWorklet worklet(extDimX,
400  extDimY,
401  extDimZ,
402  sigDimX,
403  sigDimY,
404  sigDimZ,
405  sigStartX,
406  sigStartY,
407  sigStartZ,
408  sigPretendDimX,
409  sigPretendDimY,
410  sigPretendDimZ,
411  ext1Method,
412  dir,
413  false); // not treating input signal as having zeros
414  DispatcherType dispatcher(worklet);
415  dispatcher.Invoke(ext1, sigIn);
416  }
417 
418  // Then work on back extension
419  dir = BACK;
420  extDimX = sigPretendDimX;
421  extDimY = sigPretendDimY;
422  if (!pretendSigPaddedZero && !padZeroAtExt2)
423  {
424  extDimZ = addLen;
425  ext2.Allocate(extDimX * extDimY * extDimZ);
426  ExtensionWorklet worklet(extDimX,
427  extDimY,
428  extDimZ,
429  sigDimX,
430  sigDimY,
431  sigDimZ,
432  sigStartX,
433  sigStartY,
434  sigStartZ,
435  sigPretendDimX,
436  sigPretendDimY,
437  sigPretendDimZ,
438  ext2Method,
439  dir,
440  false);
441  DispatcherType dispatcher(worklet);
442  dispatcher.Invoke(ext2, sigIn);
443  }
444  else if (!pretendSigPaddedZero && padZeroAtExt2)
445  { // This case is not exactly padding a zero at the end of Ext2.
446  // Rather, it is to increase extension length by one and fill it
447  // to be whatever mirrored.
448  extDimZ = addLen + 1;
449  ext2.Allocate(extDimX * extDimY * extDimZ);
450  ExtensionWorklet worklet(extDimX,
451  extDimY,
452  extDimZ,
453  sigDimX,
454  sigDimY,
455  sigDimZ,
456  sigStartX,
457  sigStartY,
458  sigStartZ,
459  sigPretendDimX,
460  sigPretendDimY,
461  sigPretendDimZ,
462  ext2Method,
463  dir,
464  false);
465  DispatcherType dispatcher(worklet);
466  dispatcher.Invoke(ext2, sigIn);
467  }
468  else // pretendSigPaddedZero
469  {
470  ExtendArrayType ext2Temp;
471  extDimZ = addLen;
472  ext2Temp.Allocate(extDimX * extDimY * extDimZ);
473  ExtensionWorklet worklet(extDimX,
474  extDimY,
475  extDimZ,
476  sigDimX,
477  sigDimY,
478  sigDimZ,
479  sigStartX,
480  sigStartY,
481  sigStartZ,
482  sigPretendDimX,
483  sigPretendDimY,
484  sigPretendDimZ,
485  ext2Method,
486  dir,
487  true); // pretend sig is padded a zero
488  DispatcherType dispatcher(worklet);
489  dispatcher.Invoke(ext2Temp, sigIn);
490 
491  // Give ext2 one layer thicker to hold the pretend zeros from signal.
492  ext2.Allocate(extDimX * extDimY * (extDimZ + 1));
493  WaveletBase::DeviceCubeCopyTo(
494  ext2Temp, extDimX, extDimY, extDimZ, ext2, extDimX, extDimY, extDimZ + 1, 0, 0, 1);
495  WaveletBase::DeviceAssignZero3DPlaneZ(ext2, extDimX, extDimY, extDimZ + 1, 0);
496  }
497  return 0;
498  }
499 
500  // L[3] L[15]
501  // -----------------------
502  // / / /|
503  // L[5] / / / |
504  // / LLH / HLH / |
505  // / / / | L[16]
506  // ----------------------- |
507  // / / /| |
508  // L[2] / / / | /|
509  // / / / | / |
510  // /___L[0]___/___L[12]__/ | / | L[22]
511  // | | | |/ |
512  // L[1] | | | /HHH /
513  // | LLL | HLL | /| /
514  // | | | / | / L[23]
515  // |---------------------|/ | /
516  // | | | |/
517  // | | | /
518  // L[7] | LHL | HHL | /
519  // | | | / L[20]
520  // |__________|__________|/
521  // L[6] L[18]
522  //
523  // Performs one level of 3D discrete wavelet transform on a small cube of input array
524  // The output has the same size as the small cube
525  template <typename ArrayInType, typename ArrayOutType>
526  vtkm::Float64 DWT3D(ArrayInType& sigIn,
527  vtkm::Id sigDimX,
528  vtkm::Id sigDimY,
529  vtkm::Id sigDimZ,
530  vtkm::Id sigStartX,
531  vtkm::Id sigStartY,
532  vtkm::Id sigStartZ,
533  vtkm::Id sigPretendDimX,
534  vtkm::Id sigPretendDimY,
535  vtkm::Id sigPretendDimZ,
536  ArrayOutType& coeffOut,
537  bool discardSigIn) // discard sigIn on devices
538  {
539  std::vector<vtkm::Id> L(27, 0);
540 
541  // LLL
542  L[0] = WaveletBase::GetApproxLength(sigPretendDimX);
543  L[1] = WaveletBase::GetApproxLength(sigPretendDimY);
544  L[2] = WaveletBase::GetApproxLength(sigPretendDimZ);
545  // LLH
546  L[3] = L[0];
547  L[4] = L[1];
548  L[5] = WaveletBase::GetDetailLength(sigPretendDimZ);
549  // LHL
550  L[6] = L[0];
551  L[7] = WaveletBase::GetDetailLength(sigPretendDimY);
552  L[8] = L[2];
553  // LHH
554  L[9] = L[0];
555  L[10] = L[7];
556  L[11] = L[5];
557  // HLL
558  L[12] = WaveletBase::GetDetailLength(sigPretendDimX);
559  L[13] = L[1];
560  L[14] = L[2];
561  // HLH
562  L[15] = L[12];
563  L[16] = L[1];
564  L[17] = L[5];
565  // HHL
566  L[18] = L[12];
567  L[19] = L[7];
568  L[20] = L[2];
569  // HHH
570  L[21] = L[12];
571  L[22] = L[7];
572  L[23] = L[5];
573 
574  L[24] = sigPretendDimX;
575  L[25] = sigPretendDimY;
576  L[26] = sigPretendDimZ;
577 
578  vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
579  bool oddLow = true;
580  if (filterLen % 2 != 0)
581  {
582  oddLow = false;
583  }
584  vtkm::Id addLen = filterLen / 2;
585 
586  using ValueType = typename ArrayInType::ValueType;
587  using ArrayType = vtkm::cont::ArrayHandle<ValueType>;
588  using LeftRightXFormType = vtkm::worklet::wavelets::ForwardTransform3DLeftRight;
589  using TopDownXFormType = vtkm::worklet::wavelets::ForwardTransform3DTopDown;
590  using FrontBackXFormType = vtkm::worklet::wavelets::ForwardTransform3DFrontBack;
591  using LeftRightDispatcherType = vtkm::worklet::DispatcherMapField<LeftRightXFormType>;
592  using TopDownDispatcherType = vtkm::worklet::DispatcherMapField<TopDownXFormType>;
593  using FrontBackDispatcherType = vtkm::worklet::DispatcherMapField<FrontBackXFormType>;
594 
595  vtkm::cont::Timer timer;
596  vtkm::Float64 computationTime = 0.0;
597 
598  // First transform in X direction
599  ArrayType afterX;
600  afterX.Allocate(sigPretendDimX * sigPretendDimY * sigPretendDimZ);
601  {
602  ArrayType leftExt, rightExt;
603  this->Extend3DLeftRight(sigIn,
604  sigDimX,
605  sigDimY,
606  sigDimZ,
607  sigStartX,
608  sigStartY,
609  sigStartZ,
610  sigPretendDimX,
611  sigPretendDimY,
612  sigPretendDimZ,
613  leftExt,
614  rightExt,
615  addLen,
616  WaveletBase::wmode,
617  WaveletBase::wmode,
618  false,
619  false);
620  LeftRightXFormType worklet(filterLen,
621  L[0],
622  oddLow,
623  addLen,
624  sigPretendDimY,
625  sigPretendDimZ,
626  sigDimX,
627  sigDimY,
628  sigDimZ,
629  sigStartX,
630  sigStartY,
631  sigStartZ,
632  sigPretendDimX,
633  sigPretendDimY,
634  sigPretendDimZ,
635  addLen,
636  sigPretendDimY,
637  sigPretendDimZ);
638  LeftRightDispatcherType dispatcher(worklet);
639  timer.Start();
640  dispatcher.Invoke(leftExt,
641  sigIn,
642  rightExt,
643  WaveletBase::filter.GetLowDecomposeFilter(),
644  WaveletBase::filter.GetHighDecomposeFilter(),
645  afterX);
646  computationTime += timer.GetElapsedTime();
647  }
648 
649  if (discardSigIn)
650  {
651  sigIn.ReleaseResourcesExecution();
652  }
653 
654  // Then do transform in Y direction
655  ArrayType afterY;
656  afterY.Allocate(sigPretendDimX * sigPretendDimY * sigPretendDimZ);
657  {
658  ArrayType topExt, bottomExt;
659  this->Extend3DTopDown(afterX,
660  sigPretendDimX,
661  sigPretendDimY,
662  sigPretendDimZ,
663  0,
664  0,
665  0,
666  sigPretendDimX,
667  sigPretendDimY,
668  sigPretendDimZ,
669  topExt,
670  bottomExt,
671  addLen,
672  WaveletBase::wmode,
673  WaveletBase::wmode,
674  false,
675  false);
676  TopDownXFormType worklet(filterLen,
677  L[1],
678  oddLow,
679  sigPretendDimX,
680  addLen,
681  sigPretendDimZ,
682  sigPretendDimX,
683  sigPretendDimY,
684  sigPretendDimZ,
685  0,
686  0,
687  0,
688  sigPretendDimX,
689  sigPretendDimY,
690  sigPretendDimZ,
691  sigPretendDimX,
692  addLen,
693  sigPretendDimZ);
694  TopDownDispatcherType dispatcher(worklet);
695  timer.Start();
696  dispatcher.Invoke(topExt,
697  afterX,
698  bottomExt,
699  WaveletBase::filter.GetLowDecomposeFilter(),
700  WaveletBase::filter.GetHighDecomposeFilter(),
701  afterY);
702  computationTime += timer.GetElapsedTime();
703  }
704 
705  // Then do transform in Z direction
706  afterX.ReleaseResources(); // release afterX
707  {
708  ArrayType frontExt, backExt;
709  coeffOut.Allocate(sigPretendDimX * sigPretendDimY * sigPretendDimZ);
710  this->Extend3DFrontBack(afterY,
711  sigPretendDimX,
712  sigPretendDimY,
713  sigPretendDimZ,
714  0,
715  0,
716  0,
717  sigPretendDimX,
718  sigPretendDimY,
719  sigPretendDimZ,
720  frontExt,
721  backExt,
722  addLen,
723  WaveletBase::wmode,
724  WaveletBase::wmode,
725  false,
726  false);
727  FrontBackXFormType worklet(filterLen,
728  L[1],
729  oddLow,
730  sigPretendDimX,
731  sigPretendDimY,
732  addLen,
733  sigPretendDimX,
734  sigPretendDimY,
735  sigPretendDimZ,
736  0,
737  0,
738  0,
739  sigPretendDimX,
740  sigPretendDimY,
741  sigPretendDimZ,
742  sigPretendDimX,
743  sigPretendDimY,
744  addLen);
745  FrontBackDispatcherType dispatcher(worklet);
746  timer.Start();
747  dispatcher.Invoke(frontExt,
748  afterY,
749  backExt,
750  WaveletBase::filter.GetLowDecomposeFilter(),
751  WaveletBase::filter.GetHighDecomposeFilter(),
752  coeffOut);
753  computationTime += timer.GetElapsedTime();
754  }
755 
756  return computationTime;
757  }
758 
759  // Performs one level of IDWT on a small cube of a big cube
760  // The output array has the same dimensions as the small cube.
761  template <typename ArrayInType, typename ArrayOutType>
762  vtkm::Float64 IDWT3D(ArrayInType& coeffIn,
763  vtkm::Id inDimX,
764  vtkm::Id inDimY,
765  vtkm::Id inDimZ,
766  vtkm::Id inStartX,
767  vtkm::Id inStartY,
768  vtkm::Id inStartZ,
769  const std::vector<vtkm::Id>& L,
770  ArrayOutType& sigOut,
771  bool discardCoeffIn) // can we discard coeffIn?
772  {
773  //VTKM_ASSERT( L.size() == 27 );
774  //VTKM_ASSERT( inDimX * inDimY * inDimZ == coeffIn.GetNumberOfValues() );
775  vtkm::Id inPretendDimX = L[0] + L[12];
776  vtkm::Id inPretendDimY = L[1] + L[7];
777  vtkm::Id inPretendDimZ = L[2] + L[5];
778 
779  vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
780  using BasicArrayType = vtkm::cont::ArrayHandle<typename ArrayInType::ValueType>;
781  using LeftRightXFormType = vtkm::worklet::wavelets::InverseTransform3DLeftRight;
782  using TopDownXFormType = vtkm::worklet::wavelets::InverseTransform3DTopDown;
783  using FrontBackXFormType = vtkm::worklet::wavelets::InverseTransform3DFrontBack;
784  using LeftRightDispatcherType = vtkm::worklet::DispatcherMapField<LeftRightXFormType>;
785  using TopDownDispatcherType = vtkm::worklet::DispatcherMapField<TopDownXFormType>;
786  using FrontBackDispatcherType = vtkm::worklet::DispatcherMapField<FrontBackXFormType>;
787 
788  vtkm::cont::Timer timer;
789  vtkm::Float64 computationTime = 0.0;
790 
791  // First, inverse transform in Z direction
792  BasicArrayType afterZ;
793  afterZ.Allocate(inPretendDimX * inPretendDimY * inPretendDimZ);
794  {
795  BasicArrayType ext1, ext2, ext3, ext4;
796  vtkm::Id extDimX = inPretendDimX;
797  vtkm::Id extDimY = inPretendDimY;
798  vtkm::Id ext1DimZ = 0, ext2DimZ = 0, ext3DimZ = 0, ext4DimZ = 0;
799  this->IDWTHelper3DFrontBack(coeffIn,
800  inDimX,
801  inDimY,
802  inDimZ,
803  inStartX,
804  inStartY,
805  inStartZ,
806  inPretendDimX,
807  inPretendDimY,
808  inPretendDimZ,
809  L[2],
810  L[5],
811  ext1,
812  ext2,
813  ext3,
814  ext4,
815  ext1DimZ,
816  ext2DimZ,
817  ext3DimZ,
818  ext4DimZ,
819  filterLen,
820  wmode);
821  FrontBackXFormType worklet(filterLen,
822  extDimX,
823  extDimY,
824  ext1DimZ, // ext1
825  extDimX,
826  extDimY,
827  ext2DimZ, // ext2
828  extDimX,
829  extDimY,
830  ext3DimZ, // ext3
831  extDimX,
832  extDimY,
833  ext4DimZ, // ext4
834  inPretendDimX,
835  inPretendDimY,
836  L[2], // cA
837  inPretendDimX,
838  inPretendDimY,
839  L[5], // cD
840  inDimX,
841  inDimY,
842  inDimZ, // coeffIn
843  inStartX,
844  inStartY,
845  inStartZ); // coeffIn
846  FrontBackDispatcherType dispatcher(worklet);
847  timer.Start();
848  dispatcher.Invoke(ext1,
849  ext2,
850  ext3,
851  ext4,
852  coeffIn,
853  WaveletBase::filter.GetLowReconstructFilter(),
854  WaveletBase::filter.GetHighReconstructFilter(),
855  afterZ);
856  computationTime += timer.GetElapsedTime();
857  }
858 
859  if (discardCoeffIn)
860  {
861  coeffIn.ReleaseResourcesExecution();
862  }
863 
864  // Second, inverse transform in Y direction
865  BasicArrayType afterY;
866  afterY.Allocate(inPretendDimX * inPretendDimY * inPretendDimZ);
867  {
868  BasicArrayType ext1, ext2, ext3, ext4;
869  vtkm::Id extDimX = inPretendDimX;
870  vtkm::Id extDimZ = inPretendDimZ;
871  vtkm::Id ext1DimY = 0, ext2DimY = 0, ext3DimY = 0, ext4DimY = 0;
872  this->IDWTHelper3DTopDown(afterZ,
873  inPretendDimX,
874  inPretendDimY,
875  inPretendDimZ,
876  0,
877  0,
878  0,
879  inPretendDimX,
880  inPretendDimY,
881  inPretendDimZ,
882  L[1],
883  L[7],
884  ext1,
885  ext2,
886  ext3,
887  ext4,
888  ext1DimY,
889  ext2DimY,
890  ext3DimY,
891  ext4DimY,
892  filterLen,
893  wmode);
894  TopDownXFormType worklet(filterLen,
895  extDimX,
896  ext1DimY,
897  extDimZ, // ext1
898  extDimX,
899  ext2DimY,
900  extDimZ, // ext2
901  extDimX,
902  ext3DimY,
903  extDimZ, // ext3
904  extDimX,
905  ext4DimY,
906  extDimZ, // ext4
907  inPretendDimX,
908  L[1],
909  inPretendDimZ, // cA
910  inPretendDimX,
911  L[7],
912  inPretendDimZ, // cD
913  inPretendDimX,
914  inPretendDimY,
915  inPretendDimZ, // actual signal
916  0,
917  0,
918  0);
919  TopDownDispatcherType dispatcher(worklet);
920  timer.Start();
921  dispatcher.Invoke(ext1,
922  ext2,
923  ext3,
924  ext4,
925  afterZ,
926  WaveletBase::filter.GetLowReconstructFilter(),
927  WaveletBase::filter.GetHighReconstructFilter(),
928  afterY);
929  computationTime += timer.GetElapsedTime();
930  }
931 
932  // Lastly, inverse transform in X direction
933  afterZ.ReleaseResources();
934  {
935  BasicArrayType ext1, ext2, ext3, ext4;
936  vtkm::Id extDimY = inPretendDimY;
937  vtkm::Id extDimZ = inPretendDimZ;
938  vtkm::Id ext1DimX = 0, ext2DimX = 0, ext3DimX = 0, ext4DimX = 0;
939  this->IDWTHelper3DLeftRight(afterY,
940  inPretendDimX,
941  inPretendDimY,
942  inPretendDimZ,
943  0,
944  0,
945  0,
946  inPretendDimX,
947  inPretendDimY,
948  inPretendDimZ,
949  L[0],
950  L[12],
951  ext1,
952  ext2,
953  ext3,
954  ext4,
955  ext1DimX,
956  ext2DimX,
957  ext3DimX,
958  ext4DimX,
959  filterLen,
960  wmode);
961  sigOut.Allocate(inPretendDimX * inPretendDimY * inPretendDimZ);
962  LeftRightXFormType worklet(filterLen,
963  ext1DimX,
964  extDimY,
965  extDimZ, // ext1
966  ext2DimX,
967  extDimY,
968  extDimZ, // ext2
969  ext3DimX,
970  extDimY,
971  extDimZ, // ext3
972  ext4DimX,
973  extDimY,
974  extDimZ, // ext4
975  L[0],
976  inPretendDimY,
977  inPretendDimZ, // cA
978  L[12],
979  inPretendDimY,
980  inPretendDimZ, // cD
981  inPretendDimX,
982  inPretendDimY,
983  inPretendDimZ, // actual signal
984  0,
985  0,
986  0);
987  LeftRightDispatcherType dispatcher(worklet);
988  timer.Start();
989  dispatcher.Invoke(ext1,
990  ext2,
991  ext3,
992  ext4,
993  afterY,
994  WaveletBase::filter.GetLowReconstructFilter(),
995  WaveletBase::filter.GetHighReconstructFilter(),
996  sigOut);
997  computationTime += timer.GetElapsedTime();
998  }
999 
1000  return computationTime;
1001  }
1002 
1003  //=============================================================================
1004 
1005  template <typename SigInArrayType, typename ExtensionArrayType>
1006  vtkm::Id Extend2D(const SigInArrayType& sigIn, // Input
1007  vtkm::Id sigDimX,
1008  vtkm::Id sigDimY,
1009  vtkm::Id sigStartX,
1010  vtkm::Id sigStartY,
1011  vtkm::Id sigPretendDimX,
1012  vtkm::Id sigPretendDimY,
1013  ExtensionArrayType& ext1, // left/top extension
1014  ExtensionArrayType& ext2, // right/bottom extension
1015  vtkm::Id addLen,
1016  vtkm::worklet::wavelets::DWTMode ext1Method,
1017  vtkm::worklet::wavelets::DWTMode ext2Method,
1018  bool pretendSigPaddedZero,
1019  bool padZeroAtExt2,
1020  bool modeLR) // true = left-right, false = top-down
1021  {
1022  // pretendSigPaddedZero and padZeroAtExt2 cannot happen at the same time
1023  VTKM_ASSERT(!pretendSigPaddedZero || !padZeroAtExt2);
1024 
1025  if (addLen == 0) // Haar kernel
1026  {
1027  ext1.Allocate(0); // no need to extend on left/top
1028  if (pretendSigPaddedZero || padZeroAtExt2)
1029  {
1030  if (modeLR) // right extension
1031  {
1032  ext2.Allocate(sigPretendDimY);
1033  WaveletBase::DeviceAssignZero2DColumn(ext2, 1, sigPretendDimY, 0);
1034  }
1035  else // bottom extension
1036  {
1037  ext2.Allocate(sigPretendDimX);
1038  WaveletBase::DeviceAssignZero2DRow(ext2, sigPretendDimX, 1, 0);
1039  }
1040  }
1041  else
1042  {
1043  ext2.Allocate(0);
1044  }
1045  return 0;
1046  }
1047 
1048  using ValueType = typename SigInArrayType::ValueType;
1049  using ExtendArrayType = vtkm::cont::ArrayHandle<ValueType>;
1050  using ExtensionWorklet = vtkm::worklet::wavelets::ExtensionWorklet2D;
1051  using DispatcherType = typename vtkm::worklet::DispatcherMapField<ExtensionWorklet>;
1052  vtkm::Id extDimX, extDimY;
1053  vtkm::worklet::wavelets::ExtensionDirection dir;
1054 
1055  { // Work on left/top extension
1056  if (modeLR)
1057  {
1058  dir = LEFT;
1059  extDimX = addLen;
1060  extDimY = sigPretendDimY;
1061  }
1062  else
1063  {
1064  dir = TOP;
1065  extDimX = sigPretendDimX;
1066  extDimY = addLen;
1067  }
1068  ext1.Allocate(extDimX * extDimY);
1069  ExtensionWorklet worklet(extDimX,
1070  extDimY,
1071  sigDimX,
1072  sigDimY,
1073  sigStartX,
1074  sigStartY,
1075  sigPretendDimX,
1076  sigPretendDimY, // use this area
1077  ext1Method,
1078  dir,
1079  false); // not treating sigIn as having zeros
1080  DispatcherType dispatcher(worklet);
1081  dispatcher.Invoke(ext1, sigIn);
1082  }
1083 
1084  // Work on right/bottom extension
1085  if (!pretendSigPaddedZero && !padZeroAtExt2)
1086  {
1087  if (modeLR)
1088  {
1089  dir = RIGHT;
1090  extDimX = addLen;
1091  extDimY = sigPretendDimY;
1092  }
1093  else
1094  {
1095  dir = BOTTOM;
1096  extDimX = sigPretendDimX;
1097  extDimY = addLen;
1098  }
1099  ext2.Allocate(extDimX * extDimY);
1100  ExtensionWorklet worklet(extDimX,
1101  extDimY,
1102  sigDimX,
1103  sigDimY,
1104  sigStartX,
1105  sigStartY,
1106  sigPretendDimX,
1107  sigPretendDimY, // use this area
1108  ext2Method,
1109  dir,
1110  false);
1111  DispatcherType dispatcher(worklet);
1112  dispatcher.Invoke(ext2, sigIn);
1113  }
1114  else if (!pretendSigPaddedZero && padZeroAtExt2)
1115  {
1116  if (modeLR)
1117  {
1118  dir = RIGHT;
1119  extDimX = addLen + 1;
1120  extDimY = sigPretendDimY;
1121  }
1122  else
1123  {
1124  dir = BOTTOM;
1125  extDimX = sigPretendDimX;
1126  extDimY = addLen + 1;
1127  }
1128  ext2.Allocate(extDimX * extDimY);
1129  ExtensionWorklet worklet(extDimX,
1130  extDimY,
1131  sigDimX,
1132  sigDimY,
1133  sigStartX,
1134  sigStartY,
1135  sigPretendDimX,
1136  sigPretendDimY,
1137  ext2Method,
1138  dir,
1139  false);
1140  DispatcherType dispatcher(worklet);
1141  dispatcher.Invoke(ext2, sigIn);
1142  /* Pad a zero at the end of cDTemp, when cDTemp is forced to have the same
1143  length as cATemp. For example, with odd length signal, cA is 1 element
1144  longer than cD.
1145  */
1146  /* Update 10/24/2016: the extra element of cD shouldn't be zero, just be
1147  * whatever it extends to be.
1148  * if( modeLR )
1149  * WaveletBase::DeviceAssignZero2DColumn( ext2, extDimX, extDimY,
1150  * extDimX-1 );
1151  * else
1152  * WaveletBase::DeviceAssignZero2DRow( ext2, extDimX, extDimY,
1153  * extDimY-1 );
1154  */
1155  }
1156  else // pretendSigPaddedZero
1157  {
1158  ExtendArrayType ext2Temp;
1159  if (modeLR)
1160  {
1161  dir = RIGHT;
1162  extDimX = addLen;
1163  extDimY = sigPretendDimY;
1164  }
1165  else
1166  {
1167  dir = BOTTOM;
1168  extDimX = sigPretendDimX;
1169  extDimY = addLen;
1170  }
1171  ext2Temp.Allocate(extDimX * extDimY);
1172  ExtensionWorklet worklet(extDimX,
1173  extDimY,
1174  sigDimX,
1175  sigDimY,
1176  sigStartX,
1177  sigStartY,
1178  sigPretendDimX,
1179  sigPretendDimY,
1180  ext2Method,
1181  dir,
1182  true); // pretend sig is padded a zero
1183  DispatcherType dispatcher(worklet);
1184  dispatcher.Invoke(ext2Temp, sigIn);
1185 
1186  if (modeLR)
1187  {
1188  ext2.Allocate((extDimX + 1) * extDimY);
1189  WaveletBase::DeviceRectangleCopyTo(
1190  ext2Temp, extDimX, extDimY, ext2, extDimX + 1, extDimY, 1, 0);
1191  WaveletBase::DeviceAssignZero2DColumn(ext2, extDimX + 1, extDimY, 0);
1192  }
1193  else
1194  {
1195  ext2.Allocate(extDimX * (extDimY + 1));
1196  WaveletBase::DeviceRectangleCopyTo(
1197  ext2Temp, extDimX, extDimY, ext2, extDimX, extDimY + 1, 0, 1);
1198  WaveletBase::DeviceAssignZero2DRow(ext2, extDimX, extDimY + 1, 0);
1199  }
1200  }
1201  return 0;
1202  }
1203 
1204  // Extend 1D signal
1205  template <typename SigInArrayType, typename SigExtendedArrayType>
1206  vtkm::Id Extend1D(const SigInArrayType& sigIn, // Input
1207  SigExtendedArrayType& sigOut, // Output
1208  vtkm::Id addLen,
1209  vtkm::worklet::wavelets::DWTMode leftExtMethod,
1210  vtkm::worklet::wavelets::DWTMode rightExtMethod,
1211  bool attachZeroRightLeft,
1212  bool attachZeroRightRight)
1213  {
1214  // "right extension" can be attached a zero on either end, but not both ends.
1215  VTKM_ASSERT(!attachZeroRightRight || !attachZeroRightLeft);
1216 
1217  using ValueType = typename SigInArrayType::ValueType;
1218  using ExtensionArrayType = vtkm::cont::ArrayHandle<ValueType>;
1219  using ArrayConcat = vtkm::cont::ArrayHandleConcatenate<ExtensionArrayType, SigInArrayType>;
1220 
1221  ExtensionArrayType leftExtend, rightExtend;
1222 
1223  if (addLen == 0) // Haar kernel
1224  {
1225  if (attachZeroRightLeft || attachZeroRightRight)
1226  {
1227  leftExtend.Allocate(0);
1228  rightExtend.Allocate(1);
1229  WaveletBase::DeviceAssignZero(rightExtend, 0);
1230  }
1231  else
1232  {
1233  leftExtend.Allocate(0);
1234  rightExtend.Allocate(0);
1235  }
1236  ArrayConcat leftOn(leftExtend, sigIn);
1237  sigOut = vtkm::cont::make_ArrayHandleConcatenate(leftOn, rightExtend);
1238  return 0;
1239  }
1240 
1241  leftExtend.Allocate(addLen);
1242  vtkm::Id sigInLen = sigIn.GetNumberOfValues();
1243 
1244  using LeftSYMH = vtkm::worklet::wavelets::LeftSYMHExtentionWorklet;
1245  using LeftSYMW = vtkm::worklet::wavelets::LeftSYMWExtentionWorklet;
1246  using RightSYMH = vtkm::worklet::wavelets::RightSYMHExtentionWorklet;
1247  using RightSYMW = vtkm::worklet::wavelets::RightSYMWExtentionWorklet;
1248  using LeftASYMH = vtkm::worklet::wavelets::LeftASYMHExtentionWorklet;
1249  using LeftASYMW = vtkm::worklet::wavelets::LeftASYMWExtentionWorklet;
1250  using RightASYMH = vtkm::worklet::wavelets::RightASYMHExtentionWorklet;
1251  using RightASYMW = vtkm::worklet::wavelets::RightASYMWExtentionWorklet;
1252 
1253  switch (leftExtMethod)
1254  {
1255  case SYMH:
1256  {
1257  LeftSYMH worklet(addLen);
1258  vtkm::worklet::DispatcherMapField<LeftSYMH> dispatcher(worklet);
1259  dispatcher.Invoke(leftExtend, sigIn);
1260  break;
1261  }
1262  case SYMW:
1263  {
1264  LeftSYMW worklet(addLen);
1265  vtkm::worklet::DispatcherMapField<LeftSYMW> dispatcher(worklet);
1266  dispatcher.Invoke(leftExtend, sigIn);
1267  break;
1268  }
1269  case ASYMH:
1270  {
1271  LeftASYMH worklet(addLen);
1272  vtkm::worklet::DispatcherMapField<LeftASYMH> dispatcher(worklet);
1273  dispatcher.Invoke(leftExtend, sigIn);
1274  break;
1275  }
1276  case ASYMW:
1277  {
1278  LeftASYMW worklet(addLen);
1279  vtkm::worklet::DispatcherMapField<LeftASYMW> dispatcher(worklet);
1280  dispatcher.Invoke(leftExtend, sigIn);
1281  break;
1282  }
1283  default:
1284  {
1285  vtkm::cont::ErrorInternal("Left extension mode not supported!");
1286  return 1;
1287  }
1288  }
1289 
1290  if (!attachZeroRightLeft) // no attach zero, or only attach on RightRight
1291  {
1292  // Allocate memory
1293  if (attachZeroRightRight)
1294  {
1295  rightExtend.Allocate(addLen + 1);
1296  }
1297  else
1298  {
1299  rightExtend.Allocate(addLen);
1300  }
1301 
1302  switch (rightExtMethod)
1303  {
1304  case SYMH:
1305  {
1306  RightSYMH worklet(sigInLen);
1307  vtkm::worklet::DispatcherMapField<RightSYMH> dispatcher(worklet);
1308  dispatcher.Invoke(rightExtend, sigIn);
1309  break;
1310  }
1311  case SYMW:
1312  {
1313  RightSYMW worklet(sigInLen);
1314  vtkm::worklet::DispatcherMapField<RightSYMW> dispatcher(worklet);
1315  dispatcher.Invoke(rightExtend, sigIn);
1316  break;
1317  }
1318  case ASYMH:
1319  {
1320  RightASYMH worklet(sigInLen);
1321  vtkm::worklet::DispatcherMapField<RightASYMH> dispatcher(worklet);
1322  dispatcher.Invoke(rightExtend, sigIn);
1323  break;
1324  }
1325  case ASYMW:
1326  {
1327  RightASYMW worklet(sigInLen);
1328  vtkm::worklet::DispatcherMapField<RightASYMW> dispatcher(worklet);
1329  dispatcher.Invoke(rightExtend, sigIn);
1330  break;
1331  }
1332  default:
1333  {
1334  vtkm::cont::ErrorInternal("Right extension mode not supported!");
1335  return 1;
1336  }
1337  }
1338  if (attachZeroRightRight)
1339  {
1340  WaveletBase::DeviceAssignZero(rightExtend, addLen);
1341  }
1342  }
1343  else // attachZeroRightLeft mode
1344  {
1345  using ConcatArray = vtkm::cont::ArrayHandleConcatenate<SigInArrayType, ExtensionArrayType>;
1346  // attach a zero at the end of sigIn
1347  ExtensionArrayType singleValArray;
1348  singleValArray.Allocate(1);
1349  WaveletBase::DeviceAssignZero(singleValArray, 0);
1350  ConcatArray sigInPlusOne(sigIn, singleValArray);
1351 
1352  // allocate memory for extension
1353  rightExtend.Allocate(addLen);
1354 
1355  switch (rightExtMethod)
1356  {
1357  case SYMH:
1358  {
1359  RightSYMH worklet(sigInLen + 1);
1360  vtkm::worklet::DispatcherMapField<RightSYMH> dispatcher(worklet);
1361  dispatcher.Invoke(rightExtend, sigInPlusOne);
1362  break;
1363  }
1364  case SYMW:
1365  {
1366  RightSYMW worklet(sigInLen + 1);
1367  vtkm::worklet::DispatcherMapField<RightSYMW> dispatcher(worklet);
1368  dispatcher.Invoke(rightExtend, sigInPlusOne);
1369  break;
1370  }
1371  case ASYMH:
1372  {
1373  RightASYMH worklet(sigInLen + 1);
1374  vtkm::worklet::DispatcherMapField<RightASYMH> dispatcher(worklet);
1375  dispatcher.Invoke(rightExtend, sigInPlusOne);
1376  break;
1377  }
1378  case ASYMW:
1379  {
1380  RightASYMW worklet(sigInLen + 1);
1381  vtkm::worklet::DispatcherMapField<RightASYMW> dispatcher(worklet);
1382  dispatcher.Invoke(rightExtend, sigInPlusOne);
1383  break;
1384  }
1385  default:
1386  {
1387  vtkm::cont::ErrorInternal("Right extension mode not supported!");
1388  return 1;
1389  }
1390  }
1391 
1392  // make a copy of rightExtend with a zero attached to the left
1393  ExtensionArrayType rightExtendPlusOne;
1394  rightExtendPlusOne.Allocate(addLen + 1);
1395  WaveletBase::DeviceCopyStartX(rightExtend, rightExtendPlusOne, 1);
1396  WaveletBase::DeviceAssignZero(rightExtendPlusOne, 0);
1397  rightExtend = rightExtendPlusOne;
1398  }
1399 
1400  ArrayConcat leftOn(leftExtend, sigIn);
1401  sigOut = vtkm::cont::make_ArrayHandleConcatenate(leftOn, rightExtend);
1402 
1403  return 0;
1404  }
1405 
1406  // Performs one level of 1D discrete wavelet transform
1407  // It takes care of boundary conditions, etc.
1408  template <typename SignalArrayType, typename CoeffArrayType>
1409  vtkm::Float64 DWT1D(const SignalArrayType& sigIn, // Input
1410  CoeffArrayType& coeffOut, // Output: cA followed by cD
1411  std::vector<vtkm::Id>& L) // Output: how many cA and cD.
1412  {
1413  vtkm::Id sigInLen = sigIn.GetNumberOfValues();
1414  if (GetWaveletMaxLevel(sigInLen) < 1)
1415  {
1416  vtkm::cont::ErrorInternal("Signal is too short to perform DWT!");
1417  return -1;
1418  }
1419 
1420  //VTKM_ASSERT( L.size() == 3 );
1421  L[0] = WaveletBase::GetApproxLength(sigInLen);
1422  L[1] = WaveletBase::GetDetailLength(sigInLen);
1423  L[2] = sigInLen;
1424 
1425  //VTKM_ASSERT( L[0] + L[1] == L[2] );
1426 
1427  vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
1428 
1429  bool doSymConv = false;
1430  if (WaveletBase::filter.isSymmetric())
1431  {
1432  if ((WaveletBase::wmode == SYMW && (filterLen % 2 != 0)) ||
1433  (WaveletBase::wmode == SYMH && (filterLen % 2 == 0)))
1434  {
1435  doSymConv = true;
1436  }
1437  }
1438 
1439  vtkm::Id addLen; // for extension
1440  bool oddLow = true;
1441  bool oddHigh = true;
1442  if (filterLen % 2 != 0)
1443  {
1444  oddLow = false;
1445  }
1446  if (doSymConv)
1447  {
1448  addLen = filterLen / 2;
1449  }
1450  else
1451  {
1452  addLen = filterLen - 1;
1453  }
1454 
1455  vtkm::Id sigExtendedLen = sigInLen + 2 * addLen;
1456 
1457  using SigInValueType = typename SignalArrayType::ValueType;
1458  using SigInBasic = vtkm::cont::ArrayHandle<SigInValueType>;
1459 
1460  using ConcatType1 = vtkm::cont::ArrayHandleConcatenate<SigInBasic, SignalArrayType>;
1461  using ConcatType2 = vtkm::cont::ArrayHandleConcatenate<ConcatType1, SigInBasic>;
1462 
1463  ConcatType2 sigInExtended;
1464 
1465  this->Extend1D(
1466  sigIn, sigInExtended, addLen, WaveletBase::wmode, WaveletBase::wmode, false, false);
1467  //VTKM_ASSERT( sigInExtended.GetNumberOfValues() == sigExtendedLen );
1468 
1469  // initialize a worklet for forward transform
1470  vtkm::worklet::wavelets::ForwardTransform forwardTransform(
1471  filterLen, L[0], L[1], oddLow, oddHigh);
1472 
1473  coeffOut.Allocate(sigExtendedLen);
1474  vtkm::worklet::DispatcherMapField<vtkm::worklet::wavelets::ForwardTransform> dispatcher(
1475  forwardTransform);
1476  // put a timer
1477  vtkm::cont::Timer timer;
1478  timer.Start();
1479  dispatcher.Invoke(sigInExtended,
1480  WaveletBase::filter.GetLowDecomposeFilter(),
1481  WaveletBase::filter.GetHighDecomposeFilter(),
1482  coeffOut);
1483  vtkm::Float64 elapsedTime = timer.GetElapsedTime();
1484 
1485  //VTKM_ASSERT( L[0] + L[1] <= coeffOut.GetNumberOfValues() );
1486  coeffOut.Allocate(L[0] + L[1], vtkm::CopyFlag::On);
1487 
1488  return elapsedTime;
1489  }
1490 
1491  // Performs one level of inverse wavelet transform
1492  // It takes care of boundary conditions, etc.
1493  template <typename CoeffArrayType, typename SignalArrayType>
1494  vtkm::Float64 IDWT1D(const CoeffArrayType& coeffIn, // Input, cA followed by cD
1495  std::vector<vtkm::Id>& L, // Input, how many cA and cD
1496  SignalArrayType& sigOut) // Output
1497  {
1498  vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
1499  bool doSymConv = false;
1500  vtkm::worklet::wavelets::DWTMode cALeftMode = WaveletBase::wmode;
1501  vtkm::worklet::wavelets::DWTMode cARightMode = WaveletBase::wmode;
1502  vtkm::worklet::wavelets::DWTMode cDLeftMode = WaveletBase::wmode;
1503  vtkm::worklet::wavelets::DWTMode cDRightMode = WaveletBase::wmode;
1504 
1505  if (WaveletBase::filter.isSymmetric()) // this is always true with the 1st 4 filters.
1506  {
1507  if ((WaveletBase::wmode == SYMW && (filterLen % 2 != 0)) ||
1508  (WaveletBase::wmode == SYMH && (filterLen % 2 == 0)))
1509  {
1510  doSymConv = true; // doSymConv is always true with the 1st 4 filters.
1511 
1512  if (WaveletBase::wmode == SYMH)
1513  {
1514  cDLeftMode = ASYMH;
1515  if (L[2] % 2 != 0)
1516  {
1517  cARightMode = SYMW;
1518  cDRightMode = ASYMW;
1519  }
1520  else
1521  {
1522  cDRightMode = ASYMH;
1523  }
1524  }
1525  else
1526  {
1527  cDLeftMode = SYMH;
1528  if (L[2] % 2 != 0)
1529  {
1530  cARightMode = SYMW;
1531  cDRightMode = SYMH;
1532  }
1533  else
1534  {
1535  cARightMode = SYMH;
1536  }
1537  }
1538  }
1539  }
1540 
1541  vtkm::Id cATempLen, cDTempLen; //, reconTempLen;
1542  vtkm::Id addLen = 0;
1543  vtkm::Id cDPadLen = 0;
1544  if (doSymConv) // extend cA and cD
1545  {
1546  addLen = filterLen / 4; // addLen == 0 for Haar kernel
1547  if ((L[0] > L[1]) && (WaveletBase::wmode == SYMH))
1548  {
1549  cDPadLen = L[0];
1550  }
1551  cATempLen = L[0] + 2 * addLen;
1552  cDTempLen = cATempLen; // same length
1553  }
1554  else // not extend cA and cD
1555  { // (biorthogonal kernels won't come into this case)
1556  cATempLen = L[0];
1557  cDTempLen = L[1];
1558  }
1559 
1560  using IdArrayType = vtkm::cont::ArrayHandleCounting<vtkm::Id>;
1561  using PermutArrayType = vtkm::cont::ArrayHandlePermutation<IdArrayType, CoeffArrayType>;
1562 
1563  // Separate cA and cD
1564  IdArrayType approxIndices(0, 1, L[0]);
1565  IdArrayType detailIndices(L[0], 1, L[1]);
1566  PermutArrayType cA(approxIndices, coeffIn);
1567  PermutArrayType cD(detailIndices, coeffIn);
1568 
1569  using CoeffValueType = typename CoeffArrayType::ValueType;
1570  using ExtensionArrayType = vtkm::cont::ArrayHandle<CoeffValueType>;
1571  using Concat1 = vtkm::cont::ArrayHandleConcatenate<ExtensionArrayType, PermutArrayType>;
1572  using Concat2 = vtkm::cont::ArrayHandleConcatenate<Concat1, ExtensionArrayType>;
1573 
1574  Concat2 cATemp, cDTemp;
1575 
1576  if (doSymConv) // Actually extend cA and cD
1577  {
1578  // first extend cA to be cATemp
1579  this->Extend1D(cA, cATemp, addLen, cALeftMode, cARightMode, false, false);
1580 
1581  // Then extend cD based on extension needs
1582  if (cDPadLen > 0)
1583  {
1584  // Add back the missing final cD, 0.0, before doing extension
1585  this->Extend1D(cD, cDTemp, addLen, cDLeftMode, cDRightMode, true, false);
1586  }
1587  else
1588  {
1589  vtkm::Id cDTempLenWouldBe = L[1] + 2 * addLen;
1590  if (cDTempLenWouldBe == cDTempLen)
1591  {
1592  this->Extend1D(cD, cDTemp, addLen, cDLeftMode, cDRightMode, false, false);
1593  }
1594  else if (cDTempLenWouldBe == cDTempLen - 1)
1595  {
1596  this->Extend1D(cD, cDTemp, addLen, cDLeftMode, cDRightMode, false, true);
1597  }
1598  else
1599  {
1600  vtkm::cont::ErrorInternal("cDTemp Length not match!");
1601  return 1;
1602  }
1603  }
1604  }
1605  else
1606  {
1607  // make cATemp
1608  ExtensionArrayType dummyArray;
1609  dummyArray.Allocate(0);
1610  Concat1 cALeftOn(dummyArray, cA);
1611  cATemp =
1612  vtkm::cont::make_ArrayHandleConcatenate<Concat1, ExtensionArrayType>(cALeftOn, dummyArray);
1613 
1614  // make cDTemp
1615  Concat1 cDLeftOn(dummyArray, cD);
1616  cDTemp =
1617  vtkm::cont::make_ArrayHandleConcatenate<Concat1, ExtensionArrayType>(cDLeftOn, dummyArray);
1618  }
1619 
1620  vtkm::cont::ArrayHandleConcatenate<Concat2, Concat2> coeffInExtended(cATemp, cDTemp);
1621 
1622  // Allocate memory for sigOut
1623  sigOut.Allocate(cATempLen + cDTempLen);
1624 
1625  vtkm::Float64 elapsedTime = 0;
1626  if (filterLen % 2 != 0)
1627  {
1628  vtkm::worklet::wavelets::InverseTransformOdd inverseXformOdd(filterLen, L[0], cATempLen);
1629  vtkm::worklet::DispatcherMapField<vtkm::worklet::wavelets::InverseTransformOdd> dispatcher(
1630  inverseXformOdd);
1631  // use a timer
1632  vtkm::cont::Timer timer;
1633  timer.Start();
1634  dispatcher.Invoke(coeffInExtended,
1635  WaveletBase::filter.GetLowReconstructFilter(),
1636  WaveletBase::filter.GetHighReconstructFilter(),
1637  sigOut);
1638  elapsedTime = timer.GetElapsedTime();
1639  }
1640  else
1641  {
1642  vtkm::worklet::wavelets::InverseTransformEven inverseXformEven(
1643  filterLen, L[0], cATempLen, !doSymConv);
1644  vtkm::worklet::DispatcherMapField<vtkm::worklet::wavelets::InverseTransformEven> dispatcher(
1645  inverseXformEven);
1646  // use a timer
1647  vtkm::cont::Timer timer;
1648  timer.Start();
1649  dispatcher.Invoke(coeffInExtended,
1650  WaveletBase::filter.GetLowReconstructFilter(),
1651  WaveletBase::filter.GetHighReconstructFilter(),
1652  sigOut);
1653  elapsedTime = timer.GetElapsedTime();
1654  }
1655 
1656  sigOut.Allocate(L[2], vtkm::CopyFlag::On);
1657 
1658  return elapsedTime;
1659  }
1660 
1661  // Performs one level of 2D discrete wavelet transform
1662  // It takes care of boundary conditions, etc.
1663  // N.B.
1664  // L[0] == L[2]
1665  // L[1] == L[5]
1666  // L[3] == L[7]
1667  // L[4] == L[6]
1668  //
1669  // ____L[0]_______L[4]____
1670  // | | |
1671  // L[1] | cA | cDv | L[5]
1672  // | (LL) | (HL) |
1673  // | | |
1674  // |---------------------|
1675  // | | |
1676  // | cDh | cDd | L[7]
1677  // L[3] | (LH) | (HH) |
1678  // | | |
1679  // |__________|__________|
1680  // L[2] L[6]
1681  //
1682  // Performs one level of 2D discrete wavelet transform on a small rectangle of input array
1683  // The output has the same size as the small rectangle
1684  template <typename ArrayInType, typename ArrayOutType>
1685  vtkm::Float64 DWT2D(const ArrayInType& sigIn,
1686  vtkm::Id sigDimX,
1687  vtkm::Id sigDimY,
1688  vtkm::Id sigStartX,
1689  vtkm::Id sigStartY,
1690  vtkm::Id sigPretendDimX,
1691  vtkm::Id sigPretendDimY,
1692  ArrayOutType& coeffOut,
1693  std::vector<vtkm::Id>& L)
1694  {
1695  L[0] = WaveletBase::GetApproxLength(sigPretendDimX);
1696  L[2] = L[0];
1697  L[1] = WaveletBase::GetApproxLength(sigPretendDimY);
1698  L[5] = L[1];
1699  L[3] = WaveletBase::GetDetailLength(sigPretendDimY);
1700  L[7] = L[3];
1701  L[4] = WaveletBase::GetDetailLength(sigPretendDimX);
1702  L[6] = L[4];
1703  L[8] = sigPretendDimX;
1704  L[9] = sigPretendDimY;
1705 
1706  vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
1707  bool oddLow = true;
1708  if (filterLen % 2 != 0)
1709  {
1710  oddLow = false;
1711  }
1712  vtkm::Id addLen = filterLen / 2;
1713 
1714  using ValueType = typename ArrayInType::ValueType;
1715  using ArrayType = vtkm::cont::ArrayHandle<ValueType>;
1716  using ForwardXForm = vtkm::worklet::wavelets::ForwardTransform2D;
1717  using DispatcherType = typename vtkm::worklet::DispatcherMapField<ForwardXForm>;
1718 
1719  vtkm::cont::Timer timer;
1720  vtkm::Float64 computationTime = 0.0;
1721 
1722  ArrayType afterX;
1723  afterX.Allocate(sigPretendDimX * sigPretendDimY);
1724 
1725  // First transform on rows
1726  {
1727  ArrayType leftExt, rightExt;
1728  this->Extend2D(sigIn,
1729  sigDimX,
1730  sigDimY,
1731  sigStartX,
1732  sigStartY,
1733  sigPretendDimX,
1734  sigPretendDimY,
1735  leftExt,
1736  rightExt,
1737  addLen,
1738  WaveletBase::wmode,
1739  WaveletBase::wmode,
1740  false,
1741  false,
1742  true); // Extend in left-right direction
1743  ForwardXForm worklet(filterLen,
1744  L[0],
1745  oddLow,
1746  true, // left-right
1747  addLen,
1748  sigPretendDimY,
1749  sigDimX,
1750  sigDimY,
1751  sigStartX,
1752  sigStartY,
1753  sigPretendDimX,
1754  sigPretendDimY,
1755  addLen,
1756  sigPretendDimY);
1757  DispatcherType dispatcher(worklet);
1758  timer.Start();
1759  dispatcher.Invoke(leftExt,
1760  sigIn,
1761  rightExt,
1762  WaveletBase::filter.GetLowDecomposeFilter(),
1763  WaveletBase::filter.GetHighDecomposeFilter(),
1764  afterX);
1765  computationTime += timer.GetElapsedTime();
1766  }
1767 
1768  // Then do transform in Y direction
1769  {
1770  ArrayType topExt, bottomExt;
1771  coeffOut.Allocate(sigPretendDimX * sigPretendDimY);
1772  this->Extend2D(afterX,
1773  sigPretendDimX,
1774  sigPretendDimY,
1775  0,
1776  0,
1777  sigPretendDimX,
1778  sigPretendDimY,
1779  topExt,
1780  bottomExt,
1781  addLen,
1782  WaveletBase::wmode,
1783  WaveletBase::wmode,
1784  false,
1785  false,
1786  false); // Extend in top-down direction
1787  ForwardXForm worklet(filterLen,
1788  L[1],
1789  oddLow,
1790  false, // top-down
1791  sigPretendDimX,
1792  addLen,
1793  sigPretendDimX,
1794  sigPretendDimY,
1795  0,
1796  0,
1797  sigPretendDimX,
1798  sigPretendDimY,
1799  sigPretendDimX,
1800  addLen);
1801  DispatcherType dispatcher(worklet);
1802  timer.Start();
1803  dispatcher.Invoke(topExt,
1804  afterX,
1805  bottomExt,
1806  WaveletBase::filter.GetLowDecomposeFilter(),
1807  WaveletBase::filter.GetHighDecomposeFilter(),
1808  coeffOut);
1809  computationTime += timer.GetElapsedTime();
1810  }
1811 
1812  return computationTime;
1813  }
1814 
1815  // Performs one level of IDWT.
1816  // The output array has the same dimensions as the small rectangle.
1817  template <typename ArrayInType, typename ArrayOutType>
1818  vtkm::Float64 IDWT2D(const ArrayInType& coeffIn,
1819  vtkm::Id inDimX,
1820  vtkm::Id inDimY,
1821  vtkm::Id inStartX,
1822  vtkm::Id inStartY,
1823  const std::vector<vtkm::Id>& L,
1824  ArrayOutType& sigOut)
1825  {
1826  vtkm::Id inPretendDimX = L[0] + L[4];
1827  vtkm::Id inPretendDimY = L[1] + L[3];
1828 
1829  vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
1830  using BasicArrayType = vtkm::cont::ArrayHandle<typename ArrayInType::ValueType>;
1831  using IDWT2DWorklet = vtkm::worklet::wavelets::InverseTransform2D;
1832  using Dispatcher = vtkm::worklet::DispatcherMapField<IDWT2DWorklet>;
1833  vtkm::cont::Timer timer;
1834  vtkm::Float64 computationTime = 0.0;
1835 
1836  // First inverse transform on columns
1837  BasicArrayType afterY;
1838  {
1839  BasicArrayType ext1, ext2, ext3, ext4;
1840  vtkm::Id extDimX = inPretendDimX;
1841  vtkm::Id ext1DimY = 0, ext2DimY = 0, ext3DimY = 0, ext4DimY = 0;
1842  this->IDWTHelper2DTopDown(coeffIn,
1843  inDimX,
1844  inDimY,
1845  inStartX,
1846  inStartY,
1847  inPretendDimX,
1848  inPretendDimY,
1849  L[1],
1850  L[3],
1851  ext1,
1852  ext2,
1853  ext3,
1854  ext4,
1855  ext1DimY,
1856  ext2DimY,
1857  ext3DimY,
1858  ext4DimY,
1859  filterLen,
1860  wmode);
1861 
1862  afterY.Allocate(inPretendDimX * inPretendDimY);
1863  IDWT2DWorklet worklet(filterLen,
1864  extDimX,
1865  ext1DimY, // ext1
1866  inPretendDimX,
1867  L[1], // cA
1868  extDimX,
1869  ext2DimY, // ext2
1870  extDimX,
1871  ext3DimY, // ext3
1872  inPretendDimX,
1873  L[3], // cD
1874  extDimX,
1875  ext4DimY, // ext4
1876  inDimX,
1877  inDimY, // coeffIn
1878  inStartX,
1879  inStartY, // coeffIn
1880  false); // top-down
1881  Dispatcher dispatcher(worklet);
1882  timer.Start();
1883  dispatcher.Invoke(ext1,
1884  ext2,
1885  ext3,
1886  ext4,
1887  coeffIn,
1888  WaveletBase::filter.GetLowReconstructFilter(),
1889  WaveletBase::filter.GetHighReconstructFilter(),
1890  afterY);
1891  computationTime += timer.GetElapsedTime();
1892  }
1893 
1894  // Then inverse transform on rows
1895  {
1896  BasicArrayType ext1, ext2, ext3, ext4;
1897  vtkm::Id extDimY = inPretendDimY;
1898  vtkm::Id ext1DimX = 0, ext2DimX = 0, ext3DimX = 0, ext4DimX = 0;
1899  this->IDWTHelper2DLeftRight(afterY,
1900  inPretendDimX,
1901  inPretendDimY,
1902  0,
1903  0,
1904  inPretendDimX,
1905  inPretendDimY,
1906  L[0],
1907  L[4],
1908  ext1,
1909  ext2,
1910  ext3,
1911  ext4,
1912  ext1DimX,
1913  ext2DimX,
1914  ext3DimX,
1915  ext4DimX,
1916  filterLen,
1917  wmode);
1918  sigOut.Allocate(inPretendDimX * inPretendDimY);
1919  IDWT2DWorklet worklet(filterLen,
1920  ext1DimX,
1921  extDimY, // ext1
1922  L[0],
1923  inPretendDimY, // cA
1924  ext2DimX,
1925  extDimY, // ext2
1926  ext3DimX,
1927  extDimY, // ext3
1928  L[4],
1929  inPretendDimY, // cA
1930  ext4DimX,
1931  extDimY, // ext4
1932  inPretendDimX,
1933  inPretendDimY,
1934  0,
1935  0,
1936  true); // left-right
1937  Dispatcher dispatcher(worklet);
1938  timer.Start();
1939  dispatcher.Invoke(ext1,
1940  ext2,
1941  ext3,
1942  ext4,
1943  afterY,
1944  WaveletBase::filter.GetLowReconstructFilter(),
1945  WaveletBase::filter.GetHighReconstructFilter(),
1946  sigOut);
1947  computationTime += timer.GetElapsedTime();
1948  }
1949 
1950  return computationTime;
1951  }
1952 
1953  // decides the correct extension modes for cA and cD separately,
1954  // and fill the extensions (2D matrices)
1955  template <typename ArrayInType, typename ArrayOutType>
1956  void IDWTHelper2DLeftRight(const ArrayInType& coeffIn,
1957  vtkm::Id inDimX,
1958  vtkm::Id inDimY,
1959  vtkm::Id inStartX,
1960  vtkm::Id inStartY,
1961  vtkm::Id inPretendDimX,
1962  vtkm::Id inPretendDimY,
1963  vtkm::Id cADimX,
1964  vtkm::Id cDDimX,
1965  ArrayOutType& ext1,
1966  ArrayOutType& ext2, // output
1967  ArrayOutType& ext3,
1968  ArrayOutType& ext4, // output
1969  vtkm::Id& ext1DimX,
1970  vtkm::Id& ext2DimX, // output
1971  vtkm::Id& ext3DimX,
1972  vtkm::Id& ext4DimX, // output
1973  vtkm::Id filterLen,
1974  DWTMode mode)
1975  {
1976  VTKM_ASSERT(inPretendDimX == (cADimX + cDDimX));
1977 
1978  // determine extension modes
1979  DWTMode cALeft, cARight, cDLeft, cDRight;
1980  cALeft = cARight = cDLeft = cDRight = mode;
1981  if (mode == SYMH)
1982  {
1983  cDLeft = ASYMH;
1984  if (inPretendDimX % 2 != 0)
1985  {
1986  cARight = SYMW;
1987  cDRight = ASYMW;
1988  }
1989  else
1990  {
1991  cDRight = ASYMH;
1992  }
1993  }
1994  else // mode == SYMW
1995  {
1996  cDLeft = SYMH;
1997  if (inPretendDimX % 2 != 0)
1998  {
1999  cARight = SYMW;
2000  cDRight = SYMH;
2001  }
2002  else
2003  {
2004  cARight = SYMH;
2005  }
2006  }
2007  // determine length after extension
2008  vtkm::Id cAExtendedDimX, cDExtendedDimX;
2009  vtkm::Id cDPadLen = 0;
2010  vtkm::Id addLen = filterLen / 4; // addLen == 0 for Haar kernel
2011  if ((cADimX > cDDimX) && (mode == SYMH))
2012  {
2013  cDPadLen = cADimX;
2014  }
2015  cAExtendedDimX = cADimX + 2 * addLen;
2016  cDExtendedDimX = cAExtendedDimX;
2017 
2018  // extend cA
2019  vtkm::Id cADimY = inPretendDimY;
2020  this->Extend2D(coeffIn,
2021  inDimX,
2022  inDimY,
2023  inStartX,
2024  inStartY,
2025  cADimX,
2026  cADimY,
2027  ext1,
2028  ext2,
2029  addLen,
2030  cALeft,
2031  cARight,
2032  false,
2033  false,
2034  true);
2035 
2036  ext1DimX = ext2DimX = addLen;
2037 
2038  // extend cD
2039  vtkm::Id cDDimY = inPretendDimY;
2040  if (cDPadLen > 0)
2041  {
2042  this->Extend2D(coeffIn,
2043  inDimX,
2044  inDimY,
2045  inStartX + cADimX,
2046  inStartY,
2047  cDDimX,
2048  cDDimY,
2049  ext3,
2050  ext4,
2051  addLen,
2052  cDLeft,
2053  cDRight,
2054  true,
2055  false,
2056  true);
2057  ext3DimX = addLen;
2058  ext4DimX = addLen + 1;
2059  }
2060  else
2061  {
2062  vtkm::Id cDExtendedWouldBe = cDDimX + 2 * addLen;
2063  if (cDExtendedWouldBe == cDExtendedDimX)
2064  {
2065  this->Extend2D(coeffIn,
2066  inDimX,
2067  inDimY,
2068  inStartX + cADimX,
2069  inStartY,
2070  cDDimX,
2071  cDDimY,
2072  ext3,
2073  ext4,
2074  addLen,
2075  cDLeft,
2076  cDRight,
2077  false,
2078  false,
2079  true);
2080  ext3DimX = ext4DimX = addLen;
2081  }
2082  else if (cDExtendedWouldBe == cDExtendedDimX - 1)
2083  {
2084  this->Extend2D(coeffIn,
2085  inDimX,
2086  inDimY,
2087  inStartX + cADimX,
2088  inStartY,
2089  cDDimX,
2090  cDDimY,
2091  ext3,
2092  ext4,
2093  addLen,
2094  cDLeft,
2095  cDRight,
2096  false,
2097  true,
2098  true);
2099  ext3DimX = addLen;
2100  ext4DimX = addLen + 1;
2101  }
2102  else
2103  {
2104  vtkm::cont::ErrorInternal("cDTemp Length not match!");
2105  }
2106  }
2107  }
2108 
2109  // decides the correct extension modes for cA and cD separately,
2110  // and fill the extensions (2D matrices)
2111  template <typename ArrayInType, typename ArrayOutType>
2112  void IDWTHelper2DTopDown(const ArrayInType& coeffIn,
2113  vtkm::Id inDimX,
2114  vtkm::Id inDimY,
2115  vtkm::Id inStartX,
2116  vtkm::Id inStartY,
2117  vtkm::Id inPretendDimX,
2118  vtkm::Id inPretendDimY,
2119  vtkm::Id cADimY,
2120  vtkm::Id cDDimY,
2121  ArrayOutType& ext1,
2122  ArrayOutType& ext2, // output
2123  ArrayOutType& ext3,
2124  ArrayOutType& ext4, // output
2125  vtkm::Id& ext1DimY,
2126  vtkm::Id& ext2DimY, // output
2127  vtkm::Id& ext3DimY,
2128  vtkm::Id& ext4DimY, // output
2129  vtkm::Id filterLen,
2130  DWTMode mode)
2131  {
2132  VTKM_ASSERT(inPretendDimY == (cADimY + cDDimY));
2133 
2134  // determine extension modes
2135  DWTMode cATopMode, cADownMode, cDTopMode, cDDownMode;
2136  cATopMode = cADownMode = cDTopMode = cDDownMode = mode;
2137  if (mode == SYMH)
2138  {
2139  cDTopMode = ASYMH;
2140  if (inPretendDimY % 2 != 0)
2141  {
2142  cADownMode = SYMW;
2143  cDDownMode = ASYMW;
2144  }
2145  else
2146  {
2147  cDDownMode = ASYMH;
2148  }
2149  }
2150  else // mode == SYMW
2151  {
2152  cDTopMode = SYMH;
2153  if (inPretendDimY % 2 != 0)
2154  {
2155  cADownMode = SYMW;
2156  cDDownMode = SYMH;
2157  }
2158  else
2159  {
2160  cADownMode = SYMH;
2161  }
2162  }
2163  // determine length after extension
2164  vtkm::Id cAExtendedDimY, cDExtendedDimY;
2165  vtkm::Id cDPadLen = 0;
2166  vtkm::Id addLen = filterLen / 4; // addLen == 0 for Haar kernel
2167  if ((cADimY > cDDimY) && (mode == SYMH))
2168  cDPadLen = cADimY;
2169  cAExtendedDimY = cADimY + 2 * addLen;
2170  cDExtendedDimY = cAExtendedDimY;
2171 
2172  // extend cA
2173  vtkm::Id cADimX = inPretendDimX;
2174  this->Extend2D(coeffIn,
2175  inDimX,
2176  inDimY,
2177  inStartX,
2178  inStartY,
2179  cADimX,
2180  cADimY,
2181  ext1,
2182  ext2,
2183  addLen,
2184  cATopMode,
2185  cADownMode,
2186  false,
2187  false,
2188  false);
2189  ext1DimY = ext2DimY = addLen;
2190 
2191  // extend cD
2192  vtkm::Id cDDimX = inPretendDimX;
2193  if (cDPadLen > 0)
2194  {
2195  this->Extend2D(coeffIn,
2196  inDimX,
2197  inDimY,
2198  inStartX,
2199  inStartY + cADimY,
2200  cDDimX,
2201  cDDimY,
2202  ext3,
2203  ext4,
2204  addLen,
2205  cDTopMode,
2206  cDDownMode,
2207  true,
2208  false,
2209  false);
2210  ext3DimY = addLen;
2211  ext4DimY = addLen + 1;
2212  }
2213  else
2214  {
2215  vtkm::Id cDExtendedWouldBe = cDDimY + 2 * addLen;
2216  if (cDExtendedWouldBe == cDExtendedDimY)
2217  {
2218  this->Extend2D(coeffIn,
2219  inDimX,
2220  inDimY,
2221  inStartX,
2222  inStartY + cADimY,
2223  cDDimX,
2224  cDDimY,
2225  ext3,
2226  ext4,
2227  addLen,
2228  cDTopMode,
2229  cDDownMode,
2230  false,
2231  false,
2232  false);
2233  ext3DimY = ext4DimY = addLen;
2234  }
2235  else if (cDExtendedWouldBe == cDExtendedDimY - 1)
2236  {
2237  this->Extend2D(coeffIn,
2238  inDimX,
2239  inDimY,
2240  inStartX,
2241  inStartY + cADimY,
2242  cDDimX,
2243  cDDimY,
2244  ext3,
2245  ext4,
2246  addLen,
2247  cDTopMode,
2248  cDDownMode,
2249  false,
2250  true,
2251  false);
2252  ext3DimY = addLen;
2253  ext4DimY = addLen + 1;
2254  }
2255  else
2256  {
2257  vtkm::cont::ErrorInternal("cDTemp Length not match!");
2258  }
2259  }
2260  }
2261 
2262  // decides the correct extension modes for cA and cD separately,
2263  // and fill the extensions (3D cubes)
2264  template <typename ArrayInType, typename ArrayOutType>
2265  void IDWTHelper3DLeftRight(const ArrayInType& coeffIn,
2266  vtkm::Id inDimX,
2267  vtkm::Id inDimY,
2268  vtkm::Id inDimZ,
2269  vtkm::Id inStartX,
2270  vtkm::Id inStartY,
2271  vtkm::Id inStartZ,
2272  vtkm::Id inPretendDimX,
2273  vtkm::Id inPretendDimY,
2274  vtkm::Id inPretendDimZ,
2275  vtkm::Id cADimX,
2276  vtkm::Id cDDimX,
2277  ArrayOutType& ext1,
2278  ArrayOutType& ext2, // output
2279  ArrayOutType& ext3,
2280  ArrayOutType& ext4, // output
2281  vtkm::Id& ext1DimX,
2282  vtkm::Id& ext2DimX, // output
2283  vtkm::Id& ext3DimX,
2284  vtkm::Id& ext4DimX, // output
2285  vtkm::Id filterLen,
2286  DWTMode mode)
2287  {
2288  VTKM_ASSERT(inPretendDimX == (cADimX + cDDimX));
2289 
2290  // determine extension modes
2291  DWTMode cALeftMode, cARightMode, cDLeftMode, cDRightMode;
2292  cALeftMode = cARightMode = cDLeftMode = cDRightMode = mode;
2293  if (mode == SYMH)
2294  {
2295  cDLeftMode = ASYMH;
2296  if (inPretendDimX % 2 != 0)
2297  {
2298  cARightMode = SYMW;
2299  cDRightMode = ASYMW;
2300  }
2301  else
2302  {
2303  cDRightMode = ASYMH;
2304  }
2305  }
2306  else // mode == SYMW
2307  {
2308  cDLeftMode = SYMH;
2309  if (inPretendDimX % 2 != 0)
2310  {
2311  cARightMode = SYMW;
2312  cDRightMode = SYMH;
2313  }
2314  else
2315  {
2316  cARightMode = SYMH;
2317  }
2318  }
2319 
2320  // determine length after extension
2321  vtkm::Id cAExtendedDimX, cDExtendedDimX;
2322  vtkm::Id cDPadLen = 0;
2323  vtkm::Id addLen = filterLen / 4; // addLen == 0 for Haar kernel
2324  if ((cADimX > cDDimX) && (mode == SYMH))
2325  {
2326  cDPadLen = cADimX;
2327  }
2328  cAExtendedDimX = cADimX + 2 * addLen;
2329  cDExtendedDimX = cAExtendedDimX;
2330 
2331  // extend cA
2332  vtkm::Id cADimY = inPretendDimY;
2333  vtkm::Id cADimZ = inPretendDimZ;
2334  this->Extend3DLeftRight(coeffIn,
2335  inDimX,
2336  inDimY,
2337  inDimZ,
2338  inStartX,
2339  inStartY,
2340  inStartZ,
2341  cADimX,
2342  cADimY,
2343  cADimZ,
2344  ext1,
2345  ext2,
2346  addLen,
2347  cALeftMode,
2348  cARightMode,
2349  false,
2350  false);
2351  ext1DimX = ext2DimX = addLen;
2352 
2353  // extend cD
2354  vtkm::Id cDDimY = inPretendDimY;
2355  vtkm::Id cDDimZ = inPretendDimZ;
2356  bool pretendSigPaddedZero, padZeroAtExt2;
2357  if (cDPadLen > 0)
2358  {
2359  ext3DimX = addLen;
2360  ext4DimX = addLen + 1;
2361  pretendSigPaddedZero = true;
2362  padZeroAtExt2 = false;
2363  }
2364  else
2365  {
2366  vtkm::Id cDExtendedWouldBe = cDDimX + 2 * addLen;
2367  if (cDExtendedWouldBe == cDExtendedDimX)
2368  {
2369  ext3DimX = ext4DimX = addLen;
2370  pretendSigPaddedZero = false;
2371  padZeroAtExt2 = false;
2372  }
2373  else if (cDExtendedWouldBe == cDExtendedDimX - 1)
2374  {
2375  ext3DimX = addLen;
2376  ext4DimX = addLen + 1;
2377  pretendSigPaddedZero = false;
2378  padZeroAtExt2 = true;
2379  }
2380  else
2381  {
2382  pretendSigPaddedZero = padZeroAtExt2 = false; // so the compiler doesn't complain
2383  vtkm::cont::ErrorInternal("cDTemp Length not match!");
2384  }
2385  }
2386  this->Extend3DLeftRight(coeffIn,
2387  inDimX,
2388  inDimY,
2389  inDimZ,
2390  inStartX + cADimX,
2391  inStartY,
2392  inStartZ,
2393  cDDimX,
2394  cDDimY,
2395  cDDimZ,
2396  ext3,
2397  ext4,
2398  addLen,
2399  cDLeftMode,
2400  cDRightMode,
2401  pretendSigPaddedZero,
2402  padZeroAtExt2);
2403  }
2404 
2405  template <typename ArrayInType, typename ArrayOutType>
2406  void IDWTHelper3DTopDown(const ArrayInType& coeffIn,
2407  vtkm::Id inDimX,
2408  vtkm::Id inDimY,
2409  vtkm::Id inDimZ,
2410  vtkm::Id inStartX,
2411  vtkm::Id inStartY,
2412  vtkm::Id inStartZ,
2413  vtkm::Id inPretendDimX,
2414  vtkm::Id inPretendDimY,
2415  vtkm::Id inPretendDimZ,
2416  vtkm::Id cADimY,
2417  vtkm::Id cDDimY,
2418  ArrayOutType& ext1,
2419  ArrayOutType& ext2, // output
2420  ArrayOutType& ext3,
2421  ArrayOutType& ext4, // output
2422  vtkm::Id& ext1DimY,
2423  vtkm::Id& ext2DimY, // output
2424  vtkm::Id& ext3DimY,
2425  vtkm::Id& ext4DimY, // output
2426  vtkm::Id filterLen,
2427  DWTMode mode)
2428  {
2429  VTKM_ASSERT(inPretendDimY == (cADimY + cDDimY));
2430 
2431  // determine extension modes
2432  DWTMode cATopMode, cADownMode, cDTopMode, cDDownMode;
2433  cATopMode = cADownMode = cDTopMode = cDDownMode = mode;
2434  if (mode == SYMH)
2435  {
2436  cDTopMode = ASYMH;
2437  if (inPretendDimY % 2 != 0)
2438  {
2439  cADownMode = SYMW;
2440  cDDownMode = ASYMW;
2441  }
2442  else
2443  {
2444  cDDownMode = ASYMH;
2445  }
2446  }
2447  else // mode == SYMW
2448  {
2449  cDTopMode = SYMH;
2450  if (inPretendDimY % 2 != 0)
2451  {
2452  cADownMode = SYMW;
2453  cDDownMode = SYMH;
2454  }
2455  else
2456  {
2457  cADownMode = SYMH;
2458  }
2459  }
2460 
2461  // determine length after extension
2462  vtkm::Id cAExtendedDimY, cDExtendedDimY;
2463  vtkm::Id cDPadLen = 0;
2464  vtkm::Id addLen = filterLen / 4; // addLen == 0 for Haar kernel
2465  if ((cADimY > cDDimY) && (mode == SYMH))
2466  {
2467  cDPadLen = cADimY;
2468  }
2469  cAExtendedDimY = cADimY + 2 * addLen;
2470  cDExtendedDimY = cAExtendedDimY;
2471 
2472  // extend cA
2473  vtkm::Id cADimX = inPretendDimX;
2474  vtkm::Id cADimZ = inPretendDimZ;
2475  this->Extend3DTopDown(coeffIn,
2476  inDimX,
2477  inDimY,
2478  inDimZ,
2479  inStartX,
2480  inStartY,
2481  inStartZ,
2482  cADimX,
2483  cADimY,
2484  cADimZ,
2485  ext1,
2486  ext2,
2487  addLen,
2488  cATopMode,
2489  cADownMode,
2490  false,
2491  false);
2492  ext1DimY = ext2DimY = addLen;
2493 
2494  // extend cD
2495  vtkm::Id cDDimX = inPretendDimX;
2496  vtkm::Id cDDimZ = inPretendDimZ;
2497  bool pretendSigPaddedZero, padZeroAtExt2;
2498  if (cDPadLen > 0)
2499  {
2500  ext3DimY = addLen;
2501  ext4DimY = addLen + 1;
2502  pretendSigPaddedZero = true;
2503  padZeroAtExt2 = false;
2504  }
2505  else
2506  {
2507  vtkm::Id cDExtendedWouldBe = cDDimY + 2 * addLen;
2508  if (cDExtendedWouldBe == cDExtendedDimY)
2509  {
2510  ext3DimY = ext4DimY = addLen;
2511  pretendSigPaddedZero = false;
2512  padZeroAtExt2 = false;
2513  }
2514  else if (cDExtendedWouldBe == cDExtendedDimY - 1)
2515  {
2516  ext3DimY = addLen;
2517  ext4DimY = addLen + 1;
2518  pretendSigPaddedZero = false;
2519  padZeroAtExt2 = true;
2520  }
2521  else
2522  {
2523  pretendSigPaddedZero = padZeroAtExt2 = false; // so the compiler doesn't complain
2524  vtkm::cont::ErrorInternal("cDTemp Length not match!");
2525  }
2526  }
2527  this->Extend3DTopDown(coeffIn,
2528  inDimX,
2529  inDimY,
2530  inDimZ,
2531  inStartX,
2532  inStartY + cADimY,
2533  inStartZ,
2534  cDDimX,
2535  cDDimY,
2536  cDDimZ,
2537  ext3,
2538  ext4,
2539  addLen,
2540  cDTopMode,
2541  cDDownMode,
2542  pretendSigPaddedZero,
2543  padZeroAtExt2);
2544  }
2545 
2546  template <typename ArrayInType, typename ArrayOutType>
2547  void IDWTHelper3DFrontBack(const ArrayInType& coeffIn,
2548  vtkm::Id inDimX,
2549  vtkm::Id inDimY,
2550  vtkm::Id inDimZ,
2551  vtkm::Id inStartX,
2552  vtkm::Id inStartY,
2553  vtkm::Id inStartZ,
2554  vtkm::Id inPretendDimX,
2555  vtkm::Id inPretendDimY,
2556  vtkm::Id inPretendDimZ,
2557  vtkm::Id cADimZ,
2558  vtkm::Id cDDimZ,
2559  ArrayOutType& ext1,
2560  ArrayOutType& ext2, // output
2561  ArrayOutType& ext3,
2562  ArrayOutType& ext4, // output
2563  vtkm::Id& ext1DimZ,
2564  vtkm::Id& ext2DimZ, // output
2565  vtkm::Id& ext3DimZ,
2566  vtkm::Id& ext4DimZ, // output
2567  vtkm::Id filterLen,
2568  DWTMode mode)
2569  {
2570  VTKM_ASSERT(inPretendDimZ == (cADimZ + cDDimZ));
2571 
2572  // determine extension modes
2573  DWTMode cAFrontMode, cABackMode, cDFrontMode, cDBackMode;
2574  cAFrontMode = cABackMode = cDFrontMode = cDBackMode = mode;
2575  if (mode == SYMH)
2576  {
2577  cDFrontMode = ASYMH;
2578  if (inPretendDimZ % 2 != 0)
2579  {
2580  cABackMode = SYMW;
2581  cDBackMode = ASYMW;
2582  }
2583  else
2584  {
2585  cDBackMode = ASYMH;
2586  }
2587  }
2588  else // mode == SYMW
2589  {
2590  cDFrontMode = SYMH;
2591  if (inPretendDimZ % 2 != 0)
2592  {
2593  cABackMode = SYMW;
2594  cDBackMode = SYMH;
2595  }
2596  else
2597  {
2598  cABackMode = SYMH;
2599  }
2600  }
2601 
2602  // determine length after extension
2603  vtkm::Id cAExtendedDimZ, cDExtendedDimZ;
2604  vtkm::Id cDPadLen = 0;
2605  vtkm::Id addLen = filterLen / 4; // addLen == 0 for Haar kernel
2606  if ((cADimZ > cDDimZ) && (mode == SYMH))
2607  {
2608  cDPadLen = cADimZ;
2609  }
2610  cAExtendedDimZ = cADimZ + 2 * addLen;
2611  cDExtendedDimZ = cAExtendedDimZ;
2612 
2613  // extend cA
2614  vtkm::Id cADimX = inPretendDimX;
2615  vtkm::Id cADimY = inPretendDimY;
2616  this->Extend3DFrontBack(coeffIn,
2617  inDimX,
2618  inDimY,
2619  inDimZ,
2620  inStartX,
2621  inStartY,
2622  inStartZ,
2623  cADimX,
2624  cADimY,
2625  cADimZ,
2626  ext1,
2627  ext2,
2628  addLen,
2629  cAFrontMode,
2630  cABackMode,
2631  false,
2632  false);
2633  ext1DimZ = ext2DimZ = addLen;
2634 
2635  // extend cD
2636  vtkm::Id cDDimX = inPretendDimX;
2637  vtkm::Id cDDimY = inPretendDimY;
2638  bool pretendSigPaddedZero, padZeroAtExt2;
2639  if (cDPadLen > 0)
2640  {
2641  ext3DimZ = addLen;
2642  ext4DimZ = addLen + 1;
2643  pretendSigPaddedZero = true;
2644  padZeroAtExt2 = false;
2645  }
2646  else
2647  {
2648  vtkm::Id cDExtendedWouldBe = cDDimZ + 2 * addLen;
2649  if (cDExtendedWouldBe == cDExtendedDimZ)
2650  {
2651  ext3DimZ = ext4DimZ = addLen;
2652  pretendSigPaddedZero = false;
2653  padZeroAtExt2 = false;
2654  }
2655  else if (cDExtendedWouldBe == cDExtendedDimZ - 1)
2656  {
2657  ext3DimZ = addLen;
2658  ext4DimZ = addLen + 1;
2659  pretendSigPaddedZero = false;
2660  padZeroAtExt2 = true;
2661  }
2662  else
2663  {
2664  pretendSigPaddedZero = padZeroAtExt2 = false; // so the compiler doesn't complain
2665  vtkm::cont::ErrorInternal("cDTemp Length not match!");
2666  }
2667  }
2668  this->Extend3DFrontBack(coeffIn,
2669  inDimX,
2670  inDimY,
2671  inDimZ,
2672  inStartX,
2673  inStartY,
2674  inStartZ + cADimZ,
2675  cDDimX,
2676  cDDimY,
2677  cDDimZ,
2678  ext3,
2679  ext4,
2680  addLen,
2681  cDFrontMode,
2682  cDBackMode,
2683  pretendSigPaddedZero,
2684  padZeroAtExt2);
2685  }
2686 };
2687 
2688 } // namespace wavelets
2689 } // namespace worklet
2690 } // namespace vtkm
2691 
2692 #endif
ArrayHandleConcatenate.h
vtkm::worklet::wavelets::RightSYMWExtentionWorklet
Definition: WaveletTransforms.h:3366
vtkm::worklet::wavelets::WaveletBase::GetDetailLength
vtkm::Id GetDetailLength(vtkm::Id sigInLen)
Definition: WaveletBase.h:62
vtkm::cont::ArrayHandle
Manages an array-worth of data.
Definition: ArrayHandle.h:283
vtkm::worklet::wavelets::InverseTransform3DTopDown
Definition: WaveletTransforms.h:1537
vtkm::worklet::wavelets::ASYMH
@ ASYMH
Definition: WaveletTransforms.h:30
vtkm::worklet::wavelets::WaveletBase::DeviceCubeCopyTo
void DeviceCubeCopyTo(const SmallArrayType &smallCube, vtkm::Id smallX, vtkm::Id smallY, vtkm::Id smallZ, BigArrayType &bigCube, vtkm::Id bigX, vtkm::Id bigY, vtkm::Id bigZ, vtkm::Id startX, vtkm::Id startY, vtkm::Id startZ)
Definition: WaveletBase.h:295
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
vtkm::cont::ArrayHandleConcatenate
Definition: ArrayHandleConcatenate.h:241
vtkm::worklet::wavelets::WaveletDWT::Extend3DFrontBack
vtkm::Id Extend3DFrontBack(const SigInArrayType &sigIn, vtkm::Id sigDimX, vtkm::Id sigDimY, vtkm::Id sigDimZ, vtkm::Id sigStartX, vtkm::Id sigStartY, vtkm::Id sigStartZ, vtkm::Id sigPretendDimX, vtkm::Id sigPretendDimY, vtkm::Id sigPretendDimZ, ExtensionArrayType &ext1, ExtensionArrayType &ext2, vtkm::Id addLen, vtkm::worklet::wavelets::DWTMode ext1Method, vtkm::worklet::wavelets::DWTMode ext2Method, bool pretendSigPaddedZero, bool padZeroAtExt2)
Definition: WaveletDWT.h:349
VTKM_ASSERT
#define VTKM_ASSERT(condition)
Definition: Assert.h:43
vtkm::cont::ArrayHandle::Allocate
VTKM_CONT void Allocate(vtkm::Id numberOfValues, vtkm::CopyFlag preserve, vtkm::cont::Token &token) const
Allocates an array large enough to hold the given number of values.
Definition: ArrayHandle.h:465
vtkm::cont::Timer::Start
VTKM_CONT void Start()
vtkm::worklet::wavelets::InverseTransformEven
Definition: WaveletTransforms.h:3035
vtkm::worklet::wavelets::InverseTransformOdd
Definition: WaveletTransforms.h:2953
WaveletTransforms.h
vtkm::worklet::wavelets::ExtensionWorklet2D
Definition: WaveletTransforms.h:2239
vtkm::worklet::wavelets::WaveletDWT::WaveletDWT
WaveletDWT(WaveletName name)
Definition: WaveletDWT.h:36
vtkm::worklet::wavelets::LeftASYMWExtentionWorklet
Definition: WaveletTransforms.h:3310
vtkm::worklet::wavelets::WaveletDWT::IDWTHelper2DLeftRight
void IDWTHelper2DLeftRight(const ArrayInType &coeffIn, vtkm::Id inDimX, vtkm::Id inDimY, vtkm::Id inStartX, vtkm::Id inStartY, vtkm::Id inPretendDimX, vtkm::Id inPretendDimY, vtkm::Id cADimX, vtkm::Id cDDimX, ArrayOutType &ext1, ArrayOutType &ext2, ArrayOutType &ext3, ArrayOutType &ext4, vtkm::Id &ext1DimX, vtkm::Id &ext2DimX, vtkm::Id &ext3DimX, vtkm::Id &ext4DimX, vtkm::Id filterLen, DWTMode mode)
Definition: WaveletDWT.h:1956
vtkm::worklet::wavelets::ForwardTransform2D
Definition: WaveletTransforms.h:2383
vtkm::worklet::wavelets::WaveletBase
Definition: WaveletBase.h:30
vtkm::worklet::wavelets::WaveletDWT::Extend3DLeftRight
vtkm::Id Extend3DLeftRight(const SigInArrayType &sigIn, vtkm::Id sigDimX, vtkm::Id sigDimY, vtkm::Id sigDimZ, vtkm::Id sigStartX, vtkm::Id sigStartY, vtkm::Id sigStartZ, vtkm::Id sigPretendDimX, vtkm::Id sigPretendDimY, vtkm::Id sigPretendDimZ, ExtensionArrayType &ext1, ExtensionArrayType &ext2, vtkm::Id addLen, vtkm::worklet::wavelets::DWTMode ext1Method, vtkm::worklet::wavelets::DWTMode ext2Method, bool pretendSigPaddedZero, bool padZeroAtExt2)
Definition: WaveletDWT.h:43
vtkm::worklet::wavelets::InverseTransform3DFrontBack
Definition: WaveletTransforms.h:1768
vtkm::worklet::wavelets::WaveletFilter::GetFilterLength
vtkm::Id GetFilterLength()
Definition: WaveletFilter.h:108
vtkm::worklet::wavelets::WaveletDWT::IDWT2D
vtkm::Float64 IDWT2D(const ArrayInType &coeffIn, vtkm::Id inDimX, vtkm::Id inDimY, vtkm::Id inStartX, vtkm::Id inStartY, const std::vector< vtkm::Id > &L, ArrayOutType &sigOut)
Definition: WaveletDWT.h:1818
vtkm::worklet::wavelets::WaveletDWT::DWT1D
vtkm::Float64 DWT1D(const SignalArrayType &sigIn, CoeffArrayType &coeffOut, std::vector< vtkm::Id > &L)
Definition: WaveletDWT.h:1409
vtkm::Id
vtkm::Int32 Id
Represents an ID (index into arrays).
Definition: Types.h:191
vtkm::worklet::wavelets::WaveletBase::DeviceAssignZero3DPlaneY
void DeviceAssignZero3DPlaneY(ArrayType &array, vtkm::Id dimX, vtkm::Id dimY, vtkm::Id dimZ, vtkm::Id zeroY)
Definition: WaveletBase.h:159
vtkm::worklet::wavelets::WaveletDWT::IDWTHelper3DTopDown
void IDWTHelper3DTopDown(const ArrayInType &coeffIn, vtkm::Id inDimX, vtkm::Id inDimY, vtkm::Id inDimZ, vtkm::Id inStartX, vtkm::Id inStartY, vtkm::Id inStartZ, vtkm::Id inPretendDimX, vtkm::Id inPretendDimY, vtkm::Id inPretendDimZ, vtkm::Id cADimY, vtkm::Id cDDimY, ArrayOutType &ext1, ArrayOutType &ext2, ArrayOutType &ext3, ArrayOutType &ext4, vtkm::Id &ext1DimY, vtkm::Id &ext2DimY, vtkm::Id &ext3DimY, vtkm::Id &ext4DimY, vtkm::Id filterLen, DWTMode mode)
Definition: WaveletDWT.h:2406
vtkm::worklet::wavelets::WaveletDWT::Extend3DTopDown
vtkm::Id Extend3DTopDown(const SigInArrayType &sigIn, vtkm::Id sigDimX, vtkm::Id sigDimY, vtkm::Id sigDimZ, vtkm::Id sigStartX, vtkm::Id sigStartY, vtkm::Id sigStartZ, vtkm::Id sigPretendDimX, vtkm::Id sigPretendDimY, vtkm::Id sigPretendDimZ, ExtensionArrayType &ext1, ExtensionArrayType &ext2, vtkm::Id addLen, vtkm::worklet::wavelets::DWTMode ext1Method, vtkm::worklet::wavelets::DWTMode ext2Method, bool pretendSigPaddedZero, bool padZeroAtExt2)
Definition: WaveletDWT.h:196
vtkm::worklet::wavelets::WaveletBase::DeviceRectangleCopyTo
void DeviceRectangleCopyTo(const SmallArrayType &smallRect, vtkm::Id smallX, vtkm::Id smallY, BigArrayType &bigRect, vtkm::Id bigX, vtkm::Id bigY, vtkm::Id startX, vtkm::Id startY)
Definition: WaveletBase.h:278
vtkm::worklet::wavelets::WaveletBase::GetApproxLength
vtkm::Id GetApproxLength(vtkm::Id sigInLen)
Definition: WaveletBase.h:49
vtkm::worklet::wavelets::LEFT
@ LEFT
Definition: WaveletTransforms.h:36
vtkm::worklet::wavelets::InverseTransform3DLeftRight
Definition: WaveletTransforms.h:1306
vtkm::cont::ErrorInternal
This class is thrown when VTKm detects an internal state that should never be reached.
Definition: ErrorInternal.h:26
vtkm::worklet::wavelets::BACK
@ BACK
Definition: WaveletTransforms.h:41
vtkm::worklet::DispatcherMapField
Dispatcher for worklets that inherit from WorkletMapField.
Definition: DispatcherMapField.h:25
Math.h
ArrayHandlePermutation.h
Timer.h
vtkm::worklet::wavelets::ASYMW
@ ASYMW
Definition: WaveletTransforms.h:31
vtkm::worklet::wavelets::WaveletBase::DeviceAssignZero2DRow
void DeviceAssignZero2DRow(ArrayType &array, vtkm::Id dimX, vtkm::Id dimY, vtkm::Id rowIdx)
Definition: WaveletBase.h:119
vtkm::worklet::wavelets::ForwardTransform3DLeftRight
Definition: WaveletTransforms.h:834
vtkm::worklet::wavelets::WaveletBase::DeviceAssignZero
void DeviceAssignZero(ArrayType &array, vtkm::Id index)
Definition: WaveletBase.h:109
vtkm::worklet::wavelets::ForwardTransform3DTopDown
Definition: WaveletTransforms.h:987
vtkm::cont::ArrayHandlePermutation
Implicitly permutes the values in an array.
Definition: ArrayHandlePermutation.h:227
vtkm::cont::ArrayHandleCounting< vtkm::Id >
vtkm::worklet::wavelets::SYMH
@ SYMH
Definition: WaveletTransforms.h:28
vtkm::worklet::wavelets::SYMW
@ SYMW
Definition: WaveletTransforms.h:29
vtkm::worklet::contourtree_augmented::IdArrayType
vtkm::cont::ArrayHandle< vtkm::Id > IdArrayType
Definition: filter/scalar_topology/worklet/contourtree_augmented/Types.h:90
vtkm::worklet::wavelets::ForwardTransform3DFrontBack
Definition: WaveletTransforms.h:1140
vtkm::worklet::wavelets::WaveletDWT::DWT2D
vtkm::Float64 DWT2D(const ArrayInType &sigIn, vtkm::Id sigDimX, vtkm::Id sigDimY, vtkm::Id sigStartX, vtkm::Id sigStartY, vtkm::Id sigPretendDimX, vtkm::Id sigPretendDimY, ArrayOutType &coeffOut, std::vector< vtkm::Id > &L)
Definition: WaveletDWT.h:1685
vtkm::worklet::wavelets::FRONT
@ FRONT
Definition: WaveletTransforms.h:40
vtkm::worklet::wavelets::WaveletBase::wmode
DWTMode wmode
Definition: WaveletBase.h:330
vtkm::worklet::wavelets::DWTMode
DWTMode
Definition: WaveletTransforms.h:26
vtkm::worklet::wavelets::LeftSYMHExtentionWorklet
Definition: WaveletTransforms.h:3226
vtkm::worklet::wavelets::WaveletDWT::Extend2D
vtkm::Id Extend2D(const SigInArrayType &sigIn, vtkm::Id sigDimX, vtkm::Id sigDimY, vtkm::Id sigStartX, vtkm::Id sigStartY, vtkm::Id sigPretendDimX, vtkm::Id sigPretendDimY, ExtensionArrayType &ext1, ExtensionArrayType &ext2, vtkm::Id addLen, vtkm::worklet::wavelets::DWTMode ext1Method, vtkm::worklet::wavelets::DWTMode ext2Method, bool pretendSigPaddedZero, bool padZeroAtExt2, bool modeLR)
Definition: WaveletDWT.h:1006
vtkm::cont::Timer
A class that can be used to time operations in VTK-m that might be occuring in parallel.
Definition: Timer.h:43
vtkm::cont::make_ArrayHandleConcatenate
VTKM_CONT ArrayHandleConcatenate< ArrayHandleType1, ArrayHandleType2 > make_ArrayHandleConcatenate(const ArrayHandleType1 &array1, const ArrayHandleType2 &array2)
Definition: ArrayHandleConcatenate.h:266
vtkm::CopyFlag::On
@ On
vtkm::worklet::wavelets::RightASYMHExtentionWorklet
Definition: WaveletTransforms.h:3394
vtkm::worklet::wavelets::WaveletName
WaveletName
Definition: WaveletFilter.h:28
vtkm::worklet::wavelets::LeftSYMWExtentionWorklet
Definition: WaveletTransforms.h:3254
vtkm::worklet::wavelets::WaveletDWT::IDWTHelper3DFrontBack
void IDWTHelper3DFrontBack(const ArrayInType &coeffIn, vtkm::Id inDimX, vtkm::Id inDimY, vtkm::Id inDimZ, vtkm::Id inStartX, vtkm::Id inStartY, vtkm::Id inStartZ, vtkm::Id inPretendDimX, vtkm::Id inPretendDimY, vtkm::Id inPretendDimZ, vtkm::Id cADimZ, vtkm::Id cDDimZ, ArrayOutType &ext1, ArrayOutType &ext2, ArrayOutType &ext3, ArrayOutType &ext4, vtkm::Id &ext1DimZ, vtkm::Id &ext2DimZ, vtkm::Id &ext3DimZ, vtkm::Id &ext4DimZ, vtkm::Id filterLen, DWTMode mode)
Definition: WaveletDWT.h:2547
vtkm::worklet::wavelets::RightSYMHExtentionWorklet
Definition: WaveletTransforms.h:3338
vtkm::worklet::wavelets::WaveletBase::DeviceAssignZero2DColumn
void DeviceAssignZero2DColumn(ArrayType &array, vtkm::Id dimX, vtkm::Id dimY, vtkm::Id colIdx)
Definition: WaveletBase.h:132
vtkm::worklet::wavelets::RIGHT
@ RIGHT
Definition: WaveletTransforms.h:37
vtkm::worklet::wavelets::ForwardTransform
Definition: WaveletTransforms.h:2874
vtkm::worklet::wavelets::BOTTOM
@ BOTTOM
Definition: WaveletTransforms.h:39
vtkm::worklet::wavelets::WaveletBase::DeviceAssignZero3DPlaneX
void DeviceAssignZero3DPlaneX(ArrayType &array, vtkm::Id dimX, vtkm::Id dimY, vtkm::Id dimZ, vtkm::Id zeroX)
Definition: WaveletBase.h:145
vtkm::worklet::wavelets::ExtensionWorklet3D
Definition: WaveletTransforms.h:46
vtkm::worklet::wavelets::WaveletDWT::IDWT1D
vtkm::Float64 IDWT1D(const CoeffArrayType &coeffIn, std::vector< vtkm::Id > &L, SignalArrayType &sigOut)
Definition: WaveletDWT.h:1494
vtkm::worklet::wavelets::RightASYMWExtentionWorklet
Definition: WaveletTransforms.h:3422
vtkm::worklet::wavelets::WaveletDWT::DWT3D
vtkm::Float64 DWT3D(ArrayInType &sigIn, vtkm::Id sigDimX, vtkm::Id sigDimY, vtkm::Id sigDimZ, vtkm::Id sigStartX, vtkm::Id sigStartY, vtkm::Id sigStartZ, vtkm::Id sigPretendDimX, vtkm::Id sigPretendDimY, vtkm::Id sigPretendDimZ, ArrayOutType &coeffOut, bool discardSigIn)
Definition: WaveletDWT.h:526
ArrayHandleCounting.h
vtkm::worklet::wavelets::WaveletDWT::IDWTHelper3DLeftRight
void IDWTHelper3DLeftRight(const ArrayInType &coeffIn, vtkm::Id inDimX, vtkm::Id inDimY, vtkm::Id inDimZ, vtkm::Id inStartX, vtkm::Id inStartY, vtkm::Id inStartZ, vtkm::Id inPretendDimX, vtkm::Id inPretendDimY, vtkm::Id inPretendDimZ, vtkm::Id cADimX, vtkm::Id cDDimX, ArrayOutType &ext1, ArrayOutType &ext2, ArrayOutType &ext3, ArrayOutType &ext4, vtkm::Id &ext1DimX, vtkm::Id &ext2DimX, vtkm::Id &ext3DimX, vtkm::Id &ext4DimX, vtkm::Id filterLen, DWTMode mode)
Definition: WaveletDWT.h:2265
vtkm::Float64
double Float64
Definition: Types.h:155
vtkm::worklet::wavelets::ExtensionDirection
ExtensionDirection
Definition: WaveletTransforms.h:34
vtkm::worklet::wavelets::WaveletDWT::IDWTHelper2DTopDown
void IDWTHelper2DTopDown(const ArrayInType &coeffIn, vtkm::Id inDimX, vtkm::Id inDimY, vtkm::Id inStartX, vtkm::Id inStartY, vtkm::Id inPretendDimX, vtkm::Id inPretendDimY, vtkm::Id cADimY, vtkm::Id cDDimY, ArrayOutType &ext1, ArrayOutType &ext2, ArrayOutType &ext3, ArrayOutType &ext4, vtkm::Id &ext1DimY, vtkm::Id &ext2DimY, vtkm::Id &ext3DimY, vtkm::Id &ext4DimY, vtkm::Id filterLen, DWTMode mode)
Definition: WaveletDWT.h:2112
WaveletBase.h
vtkm::worklet::wavelets::WaveletDWT::IDWT3D
vtkm::Float64 IDWT3D(ArrayInType &coeffIn, vtkm::Id inDimX, vtkm::Id inDimY, vtkm::Id inDimZ, vtkm::Id inStartX, vtkm::Id inStartY, vtkm::Id inStartZ, const std::vector< vtkm::Id > &L, ArrayOutType &sigOut, bool discardCoeffIn)
Definition: WaveletDWT.h:762
vtkm::worklet::wavelets::WaveletDWT
Definition: WaveletDWT.h:32
vtkm::worklet::wavelets::LeftASYMHExtentionWorklet
Definition: WaveletTransforms.h:3282
vtkm::worklet::wavelets::TOP
@ TOP
Definition: WaveletTransforms.h:38
vtkm::worklet::wavelets::WaveletBase::filter
WaveletFilter filter
Definition: WaveletBase.h:331
vtkm::worklet::wavelets::InverseTransform2D
Definition: WaveletTransforms.h:2569
vtkm::cont::Timer::GetElapsedTime
VTKM_CONT vtkm::Float64 GetElapsedTime() const
Get the elapsed time measured by the given device adapter.
vtkm::worklet::wavelets::WaveletBase::GetWaveletMaxLevel
vtkm::Id GetWaveletMaxLevel(vtkm::Id sigInLen)
Definition: WaveletBase.h:89
vtkm::worklet::wavelets::WaveletBase::DeviceCopyStartX
void DeviceCopyStartX(const ArrayType1 &srcArray, ArrayType2 &dstArray, vtkm::Id startIdx)
Definition: WaveletBase.h:99
vtkm::worklet::wavelets::WaveletBase::DeviceAssignZero3DPlaneZ
void DeviceAssignZero3DPlaneZ(ArrayType &array, vtkm::Id dimX, vtkm::Id dimY, vtkm::Id dimZ, vtkm::Id zeroZ)
Definition: WaveletBase.h:173
vtkm::worklet::wavelets::WaveletDWT::Extend1D
vtkm::Id Extend1D(const SigInArrayType &sigIn, SigExtendedArrayType &sigOut, vtkm::Id addLen, vtkm::worklet::wavelets::DWTMode leftExtMethod, vtkm::worklet::wavelets::DWTMode rightExtMethod, bool attachZeroRightLeft, bool attachZeroRightRight)
Definition: WaveletDWT.h:1206
Generated by   doxygen 1.8.17