Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-20 02:55:47 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity
19866b1347
vtk-m/vtkm/worklet/wavelets/WaveletDWT.h

2144 lines
88 KiB
C
Raw Normal View History

move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
IDWT2D seems to be working
2016-08-11 21:24:48 +00:00
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
#ifndef vtk_m_worklet_wavelets_waveletdwt_h
#define vtk_m_worklet_wavelets_waveletdwt_h
#include <vtkm/worklet/wavelets/WaveletBase.h>
#include <vtkm/worklet/wavelets/WaveletTransforms.h>
#include <vtkm/cont/ArrayHandleConcatenate.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/ArrayHandlePermutation.h>
#include <vtkm/Math.h>
performance tunes on IDWT worklet
2016-08-12 23:38:04 +00:00
#include <vtkm/cont/Timer.h>
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
namespace vtkm {
namespace worklet {
namespace wavelets {
class WaveletDWT : public WaveletBase
{
public:
// Constructor
fix a few warnings
2016-08-11 16:39:20 +00:00
WaveletDWT( WaveletName name ) : WaveletBase( name ) {}
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
extension3D function finished; style fixes
2017-02-17 06:23:36 +00:00
assign zero plane, need to implement the function in Base
2017-02-16 06:54:38 +00:00
// Function: extend a cube in X direction
sync
2017-02-16 00:14:33 +00:00
template< typename SigInArrayType, typename ExtensionArrayType, typename DeviceTag >
vtkm::Id Extend3DLeftRight(
assign zero plane, need to implement the function in Base
2017-02-16 06:54:38 +00:00
const SigInArrayType &sigIn, // input
sync
2017-02-16 00:14:33 +00:00
vtkm::Id sigDimX, vtkm::Id sigDimY, vtkm::Id sigDimZ,
extension3D function finished; style fixes
2017-02-17 06:23:36 +00:00
vtkm::Id sigStartX, vtkm::Id sigStartY, vtkm::Id sigStartZ,
sync
2017-02-16 00:14:33 +00:00
vtkm::Id sigPretendDimX, vtkm::Id sigPretendDimY, vtkm::Id sigPretendDimZ,
assign zero plane, need to implement the function in Base
2017-02-16 06:54:38 +00:00
ExtensionArrayType &ext1, // output
ExtensionArrayType &ext2, // output
sync
2017-02-16 00:14:33 +00:00
vtkm::Id addLen,
vtkm::worklet::wavelets::DWTMode ext1Method,
vtkm::worklet::wavelets::DWTMode ext2Method,
bool pretendSigPaddedZero,
bool padZeroAtExt2,
DeviceTag )
{
assign zero plane, need to implement the function in Base
2017-02-16 06:54:38 +00:00
// pretendSigPaddedZero and padZeroAtExt2 cannot happen at the same time
VTKM_ASSERT( !pretendSigPaddedZero || !padZeroAtExt2 );
sync
2017-02-16 00:14:33 +00:00
assign zero plane, need to implement the function in Base
2017-02-16 06:54:38 +00:00
if( addLen == 0 ) // Haar kernel
{
extension3D function finished; style fixes
2017-02-17 06:23:36 +00:00
ext1.PrepareForOutput( 0, DeviceTag() ); // No extension on the left side
if( pretendSigPaddedZero || padZeroAtExt2 ) // plane of size 1*dimY*dimZ
assign zero plane, need to implement the function in Base
2017-02-16 06:54:38 +00:00
{
extension3D function finished; style fixes
2017-02-17 06:23:36 +00:00
ext2.PrepareForOutput( sigPretendDimY * sigPretendDimZ, DeviceTag() );
WaveletBase::DeviceAssignZero3DPlaneX( ext2,
1, sigPretendDimY, sigPretendDimZ,
0,
DeviceTag() );
assign zero plane, need to implement the function in Base
2017-02-16 06:54:38 +00:00
}
else
extension3D function finished; style fixes
2017-02-17 06:23:36 +00:00
ext2.PrepareForOutput( 0, DeviceTag() ); // No extension on the right side
assign zero plane, need to implement the function in Base
2017-02-16 06:54:38 +00:00
return 0;
}
extension3D function finished; style fixes
2017-02-17 06:23:36 +00:00
typedef typename SigInArrayType::ValueType ValueType;
typedef vtkm::cont::ArrayHandle< ValueType > ExtendArrayType;
typedef vtkm::worklet::wavelets::ExtensionWorklet3D ExtensionWorklet;
typedef typename vtkm::worklet::DispatcherMapField< ExtensionWorklet, DeviceTag >
DispatcherType;
vtkm::Id extDimX, extDimY, extDimZ;
vtkm::worklet::wavelets::ExtensionDirection dir;
{ // First work on left extension
dir = LEFT;
extDimX = addLen;
extDimY = sigPretendDimY;
extDimZ = sigPretendDimZ;
ext1.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ, ext1Method,
dir,
false ); // not treating input signal as having zeros
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext1, sigIn );
}
// Then work on right extension
dir = RIGHT;
extDimY = sigPretendDimY;
extDimZ = sigPretendDimZ;
if( !pretendSigPaddedZero && !padZeroAtExt2 )
{
extDimX = addLen;
ext2.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
ext2Method,
dir,
false );
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2, sigIn );
}
else if( !pretendSigPaddedZero && padZeroAtExt2 )
{ // This case is not exactly padding a zero at the end of Ext2.
// Rather, it is to increase extension length by one and fill it
// to be whatever mirrorred.
extDimX = addLen+1;
ext2.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
ext2Method,
dir,
false );
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2, sigIn );
}
else // pretendSigPaddedZero
{
ExtendArrayType ext2Temp;
extDimX = addLen;
ext2Temp.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
ext2Method,
dir,
true ); // pretend sig is padded a zero
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2Temp, sigIn );
// Give ext2 one layer thicker to hold the pretend zeros from signal.
ext2.PrepareForOutput( (extDimX+1) * extDimY * extDimZ, DeviceTag() );
WaveletBase::DeviceCubeCopyTo( ext2Temp,
extDimX, extDimY, extDimZ,
ext2,
extDimX+1, extDimY, extDimZ,
1, 0, 0,
DeviceTag() );
WaveletBase::DeviceAssignZero3DPlaneX( ext2,
extDimX+1, extDimY, extDimZ,
0,
DeviceTag() );
}
return 0;
assign zero plane, need to implement the function in Base
2017-02-16 06:54:38 +00:00
}
extension3D function finished; style fixes
2017-02-17 06:23:36 +00:00
// Function: extend a cube in Y direction
template< typename SigInArrayType, typename ExtensionArrayType, typename DeviceTag >
vtkm::Id Extend3DTopDown(
const SigInArrayType &sigIn, // input
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, // output
ExtensionArrayType &ext2, // output
vtkm::Id addLen,
vtkm::worklet::wavelets::DWTMode ext1Method,
vtkm::worklet::wavelets::DWTMode ext2Method,
bool pretendSigPaddedZero,
bool padZeroAtExt2,
DeviceTag )
{
// pretendSigPaddedZero and padZeroAtExt2 cannot happen at the same time
VTKM_ASSERT( !pretendSigPaddedZero || !padZeroAtExt2 );
if( addLen == 0 ) // Haar kernel
{
ext1.PrepareForOutput( 0, DeviceTag() ); // No extension on the top side
if( pretendSigPaddedZero || padZeroAtExt2 ) // plane of size dimX*dimZ
{
ext2.PrepareForOutput( sigPretendDimX * 1 * sigPretendDimZ, DeviceTag() );
WaveletBase::DeviceAssignZero3DPlaneY( ext2,
sigPretendDimX, 1, sigPretendDimZ,
0,
DeviceTag() );
}
else
ext2.PrepareForOutput( 0, DeviceTag() ); // No extension on the right side
return 0;
}
typedef typename SigInArrayType::ValueType ValueType;
typedef vtkm::cont::ArrayHandle< ValueType > ExtendArrayType;
typedef vtkm::worklet::wavelets::ExtensionWorklet3D ExtensionWorklet;
typedef typename vtkm::worklet::DispatcherMapField< ExtensionWorklet, DeviceTag >
DispatcherType;
vtkm::Id extDimX, extDimY, extDimZ;
vtkm::worklet::wavelets::ExtensionDirection dir;
{ // First work on top extension
dir = TOP;
extDimX = sigPretendDimX;
extDimY = addLen;
extDimZ = sigPretendDimZ;
ext1.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ, ext1Method,
dir,
false ); // not treating input signal as having zeros
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext1, sigIn );
}
// Then work on bottom extension
dir = BOTTOM;
extDimX = sigPretendDimX;
extDimZ = sigPretendDimZ;
if( !pretendSigPaddedZero && !padZeroAtExt2 )
{
extDimY = addLen;
ext2.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
ext2Method,
dir,
false );
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2, sigIn );
}
else if( !pretendSigPaddedZero && padZeroAtExt2 )
{ // This case is not exactly padding a zero at the end of Ext2.
// Rather, it is to increase extension length by one and fill it
// to be whatever mirrorred.
extDimY = addLen + 1;
ext2.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
ext2Method,
dir,
false );
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2, sigIn );
}
else // pretendSigPaddedZero
{
ExtendArrayType ext2Temp;
extDimY = addLen;
ext2Temp.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
ext2Method,
dir,
true ); // pretend sig is padded a zero
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2Temp, sigIn );
// Give ext2 one layer thicker to hold the pretend zeros from signal.
ext2.PrepareForOutput( extDimX * (extDimY+1) * extDimZ, DeviceTag() );
WaveletBase::DeviceCubeCopyTo( ext2Temp,
extDimX, extDimY, extDimZ,
ext2,
extDimX, extDimY + 1, extDimZ,
0, 1, 0,
DeviceTag() );
WaveletBase::DeviceAssignZero3DPlaneY( ext2,
extDimX, extDimY + 1, extDimZ,
0,
DeviceTag() );
}
return 0;
}
// Function: extend a cube in Z direction
template< typename SigInArrayType, typename ExtensionArrayType, typename DeviceTag >
vtkm::Id Extend3DFrontBack(
const SigInArrayType &sigIn, // input
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, // output
ExtensionArrayType &ext2, // output
vtkm::Id addLen,
vtkm::worklet::wavelets::DWTMode ext1Method,
vtkm::worklet::wavelets::DWTMode ext2Method,
bool pretendSigPaddedZero,
bool padZeroAtExt2,
DeviceTag )
{
// pretendSigPaddedZero and padZeroAtExt2 cannot happen at the same time
VTKM_ASSERT( !pretendSigPaddedZero || !padZeroAtExt2 );
if( addLen == 0 ) // Haar kernel
{
ext1.PrepareForOutput( 0, DeviceTag() ); // No extension on the front side
if( pretendSigPaddedZero || padZeroAtExt2 ) // plane of size dimX * dimY
{
ext2.PrepareForOutput( sigPretendDimX * sigPretendDimY * 1, DeviceTag() );
WaveletBase::DeviceAssignZero3DPlaneZ( ext2,
sigPretendDimX, sigPretendDimY, 1,
0,
DeviceTag() );
}
else
ext2.PrepareForOutput( 0, DeviceTag() ); // No extension on the right side
return 0;
}
typedef typename SigInArrayType::ValueType ValueType;
typedef vtkm::cont::ArrayHandle< ValueType > ExtendArrayType;
typedef vtkm::worklet::wavelets::ExtensionWorklet3D ExtensionWorklet;
typedef typename vtkm::worklet::DispatcherMapField< ExtensionWorklet, DeviceTag >
DispatcherType;
vtkm::Id extDimX, extDimY, extDimZ;
vtkm::worklet::wavelets::ExtensionDirection dir;
{ // First work on front extension
dir = FRONT;
extDimX = sigPretendDimX;
extDimY = sigPretendDimY;
extDimZ = addLen;
ext1.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ, ext1Method,
dir,
false ); // not treating input signal as having zeros
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext1, sigIn );
}
// Then work on back extension
dir = BACK;
extDimX = sigPretendDimX;
extDimY = sigPretendDimY;
if( !pretendSigPaddedZero && !padZeroAtExt2 )
{
extDimZ = addLen;
ext2.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
ext2Method,
dir,
false );
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2, sigIn );
}
else if( !pretendSigPaddedZero && padZeroAtExt2 )
{ // This case is not exactly padding a zero at the end of Ext2.
// Rather, it is to increase extension length by one and fill it
// to be whatever mirrorred.
extDimZ = addLen + 1;
ext2.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
ext2Method,
dir,
false );
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2, sigIn );
}
else // pretendSigPaddedZero
{
ExtendArrayType ext2Temp;
extDimZ = addLen;
ext2Temp.PrepareForOutput( extDimX * extDimY * extDimZ, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, extDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
ext2Method,
dir,
true ); // pretend sig is padded a zero
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2Temp, sigIn );
// Give ext2 one layer thicker to hold the pretend zeros from signal.
ext2.PrepareForOutput( extDimX * extDimY * (extDimZ+1), DeviceTag() );
WaveletBase::DeviceCubeCopyTo( ext2Temp,
extDimX, extDimY, extDimZ,
ext2,
extDimX, extDimY, extDimZ + 1,
0, 0, 1,
DeviceTag() );
WaveletBase::DeviceAssignZero3DPlaneZ( ext2,
extDimX, extDimY, extDimZ + 1,
0,
DeviceTag() );
}
return 0;
}
sync
2017-02-16 00:14:33 +00:00
finsih DWT3D
2017-02-17 22:27:30 +00:00
// L[3] L[15]
// -----------------------
// / / /|
// L[5] / / / |
// / LLH / HLH / |
// / / / | L[16]
// ----------------------- |
// / / /| |
// L[2] / / / | /|
// / / / | / |
// /___L[0]___/___L[12]__/ | / | L[22]
// | | | |/ |
// L[1] | | | /HHH /
// | LLL | HLL | /| /
// | | | / | / L[23]
// |---------------------|/ | /
// | | | |/
// | | | /
// L[7] | LHL | HHL | /
// | | | / L[20]
// |__________|__________|/
// L[6] L[18]
//
// Performs one level of 3D discrete wavelet transform on a small cube of input array
// The output has the same size as the small cube
template< typename ArrayInType, typename ArrayOutType, typename DeviceTag >
vtkm::Float64 DWT3D ( const 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,
std::vector<vtkm::Id> &L,
DeviceTag )
{
VTKM_ASSERT( L.size() == 27 );
// LLL
L[0] = WaveletBase::GetApproxLength( sigPretendDimX );
L[1] = WaveletBase::GetApproxLength( sigPretendDimY );
L[2] = WaveletBase::GetApproxLength( sigPretendDimZ );
// LLH
L[3] = L[0];
L[4] = L[1];
L[5] = WaveletBase::GetDetailLength( sigPretendDimZ );
// LHL
L[6] = L[0];
L[7] = WaveletBase::GetDetailLength( sigPretendDimY );
L[8] = L[2];
// LHH
L[9] = L[0];
L[10] = L[7];
L[11] = L[5];
// HLL
L[12] = WaveletBase::GetDetailLength( sigPretendDimX );
L[13] = L[1];
L[14] = L[2];
// HLH
L[15] = L[12];
L[16] = L[1];
L[17] = L[5];
// HHL
L[18] = L[12];
L[19] = L[7];
L[20] = L[2];
// HHH
L[21] = L[12];
L[22] = L[7];
L[23] = L[5];
L[24] = sigPretendDimX;
L[25] = sigPretendDimY;
L[26] = sigPretendDimZ;
vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
bool oddLow = true;
if( filterLen % 2 != 0 )
oddLow = false;
vtkm::Id addLen = filterLen / 2;
typedef typename ArrayInType::ValueType ValueType;
typedef vtkm::cont::ArrayHandle<ValueType> ArrayType;
debug for 3D transforms
2017-02-21 23:56:57 +00:00
typedef vtkm::worklet::wavelets::ForwardTransform3DLeftRight< DeviceTag > LeftRightXFormType;
typedef vtkm::worklet::wavelets::ForwardTransform3DTopDown< DeviceTag > TopDownXFormType;
typedef vtkm::worklet::wavelets::ForwardTransform3DFrontBack< DeviceTag > FrontBackXFormType;
typedef vtkm::worklet::DispatcherMapField< LeftRightXFormType, DeviceTag > LeftRightDispatcherType;
typedef vtkm::worklet::DispatcherMapField< TopDownXFormType, DeviceTag > TopDownDispatcherType;
typedef vtkm::worklet::DispatcherMapField< FrontBackXFormType, DeviceTag > FrontBackDispatcherType;
finsih DWT3D
2017-02-17 22:27:30 +00:00
vtkm::cont::Timer<DeviceTag> timer;
vtkm::Float64 computationTime = 0.0;
// First transform in X direction
ArrayType afterX;
afterX.PrepareForOutput( sigPretendDimX * sigPretendDimY * sigPretendDimZ, DeviceTag() );
{
ArrayType leftExt, rightExt;
this->Extend3DLeftRight( sigIn,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
leftExt,
rightExt,
addLen,
WaveletBase::wmode,
WaveletBase::wmode,
false,
false,
DeviceTag() );
LeftRightXFormType worklet( WaveletBase::filter.GetLowDecomposeFilter(),
WaveletBase::filter.GetHighDecomposeFilter(),
filterLen,
L[0],
oddLow,
addLen, sigPretendDimY, sigPretendDimZ,
sigDimX, sigDimY, sigDimZ,
sigStartX, sigStartY, sigStartZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
addLen, sigPretendDimY, sigPretendDimZ );
LeftRightDispatcherType dispatcher(worklet);
timer.Reset();
dispatcher.Invoke( leftExt, sigIn, rightExt, afterX );
computationTime += timer.GetElapsedTime();
}
// Then do transform in Y direction
ArrayType afterY;
afterY.PrepareForOutput( sigPretendDimX * sigPretendDimY * sigPretendDimZ, DeviceTag() );
{
style optimization
2017-02-19 04:14:20 +00:00
ArrayType topExt, bottomExt;
finsih DWT3D
2017-02-17 22:27:30 +00:00
this->Extend3DTopDown( afterX,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
0, 0, 0,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
topExt,
bottomExt,
addLen,
WaveletBase::wmode,
WaveletBase::wmode,
false,
false,
DeviceTag() );
TopDownXFormType worklet( WaveletBase::filter.GetLowDecomposeFilter(),
WaveletBase::filter.GetHighDecomposeFilter(),
filterLen,
L[1],
oddLow,
sigPretendDimX, addLen, sigPretendDimZ,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
0, 0, 0,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
sigPretendDimX, addLen, sigPretendDimZ );
TopDownDispatcherType dispatcher( worklet );
timer.Reset();
dispatcher.Invoke( topExt, afterX, bottomExt, coeffOut );
computationTime += timer.GetElapsedTime();
}
// Then do transform in Z direction
afterX.ReleaseResources(); // release afterX
{
ArrayType frontExt, backExt;
coeffOut.PrepareForOutput( sigPretendDimX * sigPretendDimY * sigPretendDimZ, DeviceTag() );
this->Extend3DFrontBack( afterY,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
0, 0, 0,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
frontExt,
backExt,
addLen,
WaveletBase::wmode,
WaveletBase::wmode,
false,
false,
DeviceTag() );
FrontBackXFormType worklet( WaveletBase::filter.GetLowDecomposeFilter(),
WaveletBase::filter.GetHighDecomposeFilter(),
filterLen,
L[1],
oddLow,
sigPretendDimX, sigPretendDimY, addLen,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
0, 0, 0,
sigPretendDimX, sigPretendDimY, sigPretendDimZ,
sigPretendDimX, sigPretendDimY, addLen );
FrontBackDispatcherType dispatcher( worklet );
timer.Reset();
dispatcher.Invoke( frontExt, afterY, backExt, coeffOut );
computationTime += timer.GetElapsedTime();
}
return computationTime;
}
debug for 3D transforms
2017-02-21 23:56:57 +00:00
// Performs one level of IDWT on a small cube of a big cube
// The output array has the same dimensions as the small cube.
template< typename ArrayInType, typename ArrayOutType, typename DeviceTag >
vtkm::Float64 IDWT3D( const 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,
DeviceTag )
{
VTKM_ASSERT( L.size() == 27 );
VTKM_ASSERT( inDimX * inDimY * inDimZ == coeffIn.GetNumberOfValues() );
vtkm::Id inPretendDimX = L[0] + L[12];
vtkm::Id inPretendDimY = L[1] + L[7];
vtkm::Id inPretendDimZ = L[2] + L[5];
vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
typedef vtkm::cont::ArrayHandle<typename ArrayInType::ValueType> BasicArrayType;
typedef vtkm::worklet::wavelets::InverseTransform3DLeftRight<DeviceTag> LeftRightXFormType;
typedef vtkm::worklet::wavelets::InverseTransform3DTopDown<DeviceTag> TopDownXFormType;
typedef vtkm::worklet::wavelets::InverseTransform3DFrontBack<DeviceTag> FrontBackXFormType;
typedef vtkm::worklet::DispatcherMapField<LeftRightXFormType, DeviceTag> LeftRightDispatcherType;
typedef vtkm::worklet::DispatcherMapField<TopDownXFormType, DeviceTag> TopDownDispatcherType;
typedef vtkm::worklet::DispatcherMapField<FrontBackXFormType, DeviceTag> FrontBackDispatcherType;
//vtkm::cont::Timer<DeviceTag> timer;
//vtkm::Float64 computationTime = 0.0;
// First inverse transform in Z direction
BasicArrayType afterZ;
{
BasicArrayType ext1, ext2, ext3, ext4;
vtkm::Id extDimX = inPretendDimX;
vtkm::Id extDimY = inPretendDimY;
vtkm::Id ext1DimZ, ext2DimZ, ext3DimZ, ext4DimZ;
this->IDWTHelper3DFrontBack( coeffIn,
inDimX, inDimY, inDimZ,
inStartX, inStartY, inStartZ,
inPretendDimX, inPretendDimY, inPretendDimZ,
L[2], L[5],
ext1, ext2,
ext3, ext4,
ext1DimZ, ext2DimZ,
ext3DimZ, ext4DimZ,
filterLen,
wmode,
DeviceTag() );
afterZ.PrepareForOutput( inPretendDimX * inPretendDimY * inPretendDimZ, DeviceTag() );
FrontBackXFormType worklet(
WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
filterLen,
extDimX, extDimY, ext1DimZ, // ext1
extDimX, extDimY, ext2DimZ, // ext2
extDimX, extDimY, ext3DimZ, // ext3
extDimX, extDimY, ext4DimZ, // ext4
inPretendDimX, inPretendDimY, L[2], // cA
inPretendDimX, inPretendDimY, L[5], // cD
inDimX, inDimY, inDimZ, // coeffIn
inStartX, inStartY, inStartZ ); // coeffIn
FrontBackDispatcherType dispatcher( worklet );
//timer.Reset();
dispatcher.Invoke( ext1, ext2, ext3, ext4, coeffIn, afterZ );
//computationTime += timer.GetElapsedTime();
}
// Second inverse transform in Y direction
BasicArrayType afterY;
{
BasicArrayType ext1, ext2, ext3, ext4;
vtkm::Id extDimX = inPretendDimX;
vtkm::Id extDimZ = inPretendDimZ;
vtkm::Id ext1DimY, ext2DimY, ext3DimY, ext4DimY;
this->IDWTHelper3DTopDown( afterZ,
inPretendDimX, inPretendDimY, inPretendDimZ,
0, 0, 0,
inPretendDimX, inPretendDimY, inPretendDimZ,
L[1], L[7],
ext1, ext2,
ext3, ext4,
ext1DimY, ext2DimY,
ext3DimY, ext4DimY,
filterLen,
wmode,
DeviceTag() );
afterY.PrepareForOutput( inPretendDimX * inPretendDimY * inPretendDimZ, DeviceTag() );
TopDownXFormType worklet(
WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
filterLen,
extDimX, ext1DimY, extDimZ, // ext1
extDimX, ext2DimY, extDimZ, // ext2
extDimX, ext3DimY, extDimZ, // ext3
extDimX, ext4DimY, extDimZ, // ext4
inPretendDimX, L[1], inPretendDimZ, // cA
inPretendDimX, L[7], inPretendDimZ, // cD
inPretendDimX, inPretendDimY, inPretendDimZ, // actual signal
0, 0, 0 );
TopDownDispatcherType dispatcher( worklet );
//timer.Reset();
dispatcher.Invoke( ext1, ext2, ext3, ext4, afterZ, afterY );
//computationTime += timer.GetElapsedTime();
}
// Lastly inverse transform in X direction
afterZ.ReleaseResources();
{
BasicArrayType ext1, ext2, ext3, ext4;
vtkm::Id extDimY = inPretendDimY;
vtkm::Id extDimZ = inPretendDimZ;
vtkm::Id ext1DimX, ext2DimX, ext3DimX, ext4DimX;
this->IDWTHelper3DLeftRight( afterY,
inPretendDimX, inPretendDimY, inPretendDimZ,
0, 0, 0,
inPretendDimX, inPretendDimY, inPretendDimZ,
L[0], L[12],
ext1, ext2,
ext3, ext4,
ext1DimX, ext2DimX,
ext3DimX, ext4DimX,
filterLen,
wmode,
DeviceTag() );
sigOut.PrepareForOutput( inPretendDimX * inPretendDimY * inPretendDimZ, DeviceTag() );
LeftRightXFormType worklet(
WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
filterLen,
ext1DimX, extDimY, extDimZ, // ext1
ext2DimX, extDimY, extDimZ, // ext2
ext3DimX, extDimY, extDimZ, // ext3
ext4DimX, extDimY, extDimZ, // ext4
L[0], inPretendDimY, inPretendDimZ, // cA
L[12], inPretendDimY, inPretendDimZ, // cD
inPretendDimX, inPretendDimY, inPretendDimZ, // actual signal
0, 0, 0 );
LeftRightDispatcherType dispatcher( worklet );
dispatcher.Invoke( ext1, ext2, ext3, ext4, afterY, sigOut );
}
//return computationTime;
return 0.0;
}
//=============================================================================
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
template< typename SigInArrayType, typename ExtensionArrayType, typename DeviceTag >
style fix and code cleaning
2017-02-14 06:54:41 +00:00
vtkm::Id Extend2D (const SigInArrayType &sigIn, // Input
vtkm::Id sigDimX, vtkm::Id sigDimY,
vtkm::Id sigStartX, vtkm::Id sigStartY,
vtkm::Id sigPretendDimX, vtkm::Id sigPretendDimY,
ExtensionArrayType &ext1, // left/top extension
ExtensionArrayType &ext2, // right/bottom extension
vtkm::Id addLen,
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
vtkm::worklet::wavelets::DWTMode ext1Method,
vtkm::worklet::wavelets::DWTMode ext2Method,
bool pretendSigPaddedZero,
bool padZeroAtExt2,
bool modeLR, // true = left-right
// false = top-down
style fix and code cleaning
2017-02-14 06:54:41 +00:00
DeviceTag )
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
{
// pretendSigPaddedZero and padZeroAtExt2 cannot happen at the same time
VTKM_ASSERT( !pretendSigPaddedZero || !padZeroAtExt2 );
if( addLen == 0 ) // Haar kernel
{
assign zero plane, need to implement the function in Base
2017-02-16 06:54:38 +00:00
ext1.PrepareForOutput( 0, DeviceTag() ); // no need to extend on left/top
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
if( pretendSigPaddedZero || padZeroAtExt2 )
{
if( modeLR ) // right extension
{
ext2.PrepareForOutput( sigPretendDimY, DeviceTag() );
WaveletBase::DeviceAssignZero2DColumn( ext2, 1, sigPretendDimY, 0, DeviceTag() );
}
else // bottom extension
{
ext2.PrepareForOutput( sigPretendDimX, DeviceTag() );
WaveletBase::DeviceAssignZero2DRow( ext2, sigPretendDimX, 1, 0, DeviceTag() );
}
}
else
ext2.PrepareForOutput( 0, DeviceTag() );
return 0;
}
style fix and code cleaning
2017-02-14 06:54:41 +00:00
typedef typename SigInArrayType::ValueType ValueType;
typedef vtkm::cont::ArrayHandle< ValueType > ExtendArrayType;
typedef vtkm::worklet::wavelets::ExtensionWorklet2D ExtensionWorklet;
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
typedef typename vtkm::worklet::DispatcherMapField< ExtensionWorklet, DeviceTag >
style fix and code cleaning
2017-02-14 06:54:41 +00:00
DispatcherType;
vtkm::Id extDimX, extDimY;
vtkm::worklet::wavelets::ExtensionDirection dir;
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
style fix and code cleaning
2017-02-14 06:54:41 +00:00
{ // Work on left/top extension
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
if( modeLR )
{
dir = LEFT;
extDimX = addLen;
extDimY = sigPretendDimY;
}
else
{
dir = TOP;
extDimX = sigPretendDimX;
extDimY = addLen;
}
ext1.PrepareForOutput( extDimX * extDimY, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, sigDimX, sigDimY,
sigStartX, sigStartY, sigPretendDimX, sigPretendDimY, // use this area
ext1Method, dir, false ); // not treating sigIn as having zeros
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext1, sigIn );
}
// Work on right/bottom extension
if( !pretendSigPaddedZero && !padZeroAtExt2 )
{
if( modeLR )
{
dir = RIGHT;
extDimX = addLen;
extDimY = sigPretendDimY;
}
else
{
dir = BOTTOM;
extDimX = sigPretendDimX;
extDimY = addLen;
}
ext2.PrepareForOutput( extDimX * extDimY, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, sigDimX, sigDimY,
sigStartX, sigStartY, sigPretendDimX, sigPretendDimY, // use this area
ext2Method, dir, false );
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2, sigIn );
}
else if( !pretendSigPaddedZero && padZeroAtExt2 )
{
if( modeLR )
{
dir = RIGHT;
extDimX = addLen+1;
extDimY = sigPretendDimY;
}
else
{
dir = BOTTOM;
extDimX = sigPretendDimX;
extDimY = addLen+1;
}
ext2.PrepareForOutput( extDimX * extDimY, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, sigDimX, sigDimY,
sigStartX, sigStartY, sigPretendDimX, sigPretendDimY,
ext2Method, dir, false );
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2, sigIn );
update 2D idwt to not assign zero to the last element of cDTemp
2016-10-25 02:11:05 +00:00
/* Pad a zero at the end of cDTemp, when cDTemp is forced to have the same
length as cATemp. For example, with odd length signal, cA is 1 element
longer than cD.
*/
/* Update 10/24/2016: the extra element of cD shouldn't be zero, just be
* whatever it extends to be.
* if( modeLR )
* WaveletBase::DeviceAssignZero2DColumn( ext2, extDimX, extDimY,
* extDimX-1, DeviceTag() );
* else
* WaveletBase::DeviceAssignZero2DRow( ext2, extDimX, extDimY,
* extDimY-1, DeviceTag() );
*/
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
}
else // pretendSigPaddedZero
{
ExtendArrayType ext2Temp;
if( modeLR )
{
dir = RIGHT;
extDimX = addLen;
extDimY = sigPretendDimY;
}
else
{
dir = BOTTOM;
extDimX = sigPretendDimX;
extDimY = addLen;
}
ext2Temp.PrepareForOutput( extDimX * extDimY, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, sigDimX, sigDimY,
sigStartX, sigStartY, sigPretendDimX, sigPretendDimY,
ext2Method, dir, true ); // pretend sig is padded a zero
DispatcherType dispatcher( worklet );
dispatcher.Invoke( ext2Temp, sigIn );
if( modeLR )
{
ext2.PrepareForOutput( (extDimX+1) * extDimY, DeviceTag() );
WaveletBase::DeviceRectangleCopyTo( ext2Temp, extDimX, extDimY,
ext2, extDimX+1, extDimY,
1, 0, DeviceTag() );
WaveletBase::DeviceAssignZero2DColumn( ext2, extDimX+1, extDimY,
0, DeviceTag() );
}
else
{
ext2.PrepareForOutput( extDimX * (extDimY+1), DeviceTag() );
WaveletBase::DeviceRectangleCopyTo( ext2Temp, extDimX, extDimY,
ext2, extDimX, extDimY+1,
0, 1, DeviceTag() );
WaveletBase::DeviceAssignZero2DRow( ext2, extDimX, extDimY+1,
0, DeviceTag() );
}
}
return 0;
}
move to linux to debug
2016-09-11 22:14:01 +00:00
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
move to linux to debug
2016-09-11 22:14:01 +00:00
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
// Extend 1D signal
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
template< typename SigInArrayType, typename SigExtendedArrayType, typename DeviceTag >
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
vtkm::Id Extend1D( const SigInArrayType &sigIn, // Input
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
SigExtendedArrayType &sigOut, // Output
vtkm::Id addLen,
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
vtkm::worklet::wavelets::DWTMode leftExtMethod,
vtkm::worklet::wavelets::DWTMode rightExtMethod,
half done even length filters implementation
2016-08-15 19:46:35 +00:00
bool attachZeroRightLeft,
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
bool attachZeroRightRight,
DeviceTag )
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
{
half done even length filters implementation
2016-08-15 19:46:35 +00:00
// "right extension" can be attached a zero on either end, but not both ends.
VTKM_ASSERT( !attachZeroRightRight || !attachZeroRightLeft );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
typedef typename SigInArrayType::ValueType ValueType;
typedef vtkm::cont::ArrayHandle< ValueType > ExtensionArrayType;
fix a cuda memory bug
2016-09-07 23:43:10 +00:00
typedef vtkm::cont::ArrayHandleConcatenate< ExtensionArrayType, SigInArrayType>
ArrayConcat;
ExtensionArrayType leftExtend, rightExtend;
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
fix a cuda memory bug
2016-09-07 23:43:10 +00:00
if( addLen == 0 ) // Haar kernel
{
if( attachZeroRightLeft || attachZeroRightRight )
{
leftExtend.PrepareForOutput( 0, DeviceTag() );
rightExtend.PrepareForOutput(1, DeviceTag() );
WaveletBase::DeviceAssignZero( rightExtend, 0, DeviceTag() );
}
else
{
leftExtend.PrepareForOutput( 0, DeviceTag() );
rightExtend.PrepareForOutput(0, DeviceTag() );
}
ArrayConcat leftOn( leftExtend, sigIn );
sigOut = vtkm::cont::make_ArrayHandleConcatenate( leftOn, rightExtend );
return 0;
}
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
fix a cuda memory bug
2016-09-07 23:43:10 +00:00
leftExtend.PrepareForOutput( addLen, DeviceTag() );
use a worklet to perform array extension
2016-08-03 20:34:02 +00:00
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
use a worklet to perform array extension
2016-08-03 20:34:02 +00:00
typedef vtkm::worklet::wavelets::LeftSYMHExtentionWorklet LeftSYMH;
typedef vtkm::worklet::wavelets::LeftSYMWExtentionWorklet LeftSYMW;
typedef vtkm::worklet::wavelets::RightSYMHExtentionWorklet RightSYMH;
typedef vtkm::worklet::wavelets::RightSYMWExtentionWorklet RightSYMW;
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
typedef vtkm::worklet::wavelets::LeftASYMHExtentionWorklet LeftASYMH;
typedef vtkm::worklet::wavelets::LeftASYMWExtentionWorklet LeftASYMW;
typedef vtkm::worklet::wavelets::RightASYMHExtentionWorklet RightASYMH;
typedef vtkm::worklet::wavelets::RightASYMWExtentionWorklet RightASYMW;
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
switch( leftExtMethod )
{
working on even length filters; need ASYM* support in Extend1D()
2016-08-15 23:54:15 +00:00
case SYMH:
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
{
use a worklet to perform array extension
2016-08-03 20:34:02 +00:00
LeftSYMH worklet( addLen );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< LeftSYMH, DeviceTag > dispatcher( worklet );
use a worklet to perform array extension
2016-08-03 20:34:02 +00:00
dispatcher.Invoke( leftExtend, sigIn );
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
break;
}
working on even length filters; need ASYM* support in Extend1D()
2016-08-15 23:54:15 +00:00
case SYMW:
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
{
use a worklet to perform array extension
2016-08-03 20:34:02 +00:00
LeftSYMW worklet( addLen );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< LeftSYMW, DeviceTag > dispatcher( worklet );
use a worklet to perform array extension
2016-08-03 20:34:02 +00:00
dispatcher.Invoke( leftExtend, sigIn );
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
break;
}
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
case ASYMH:
{
LeftASYMH worklet( addLen );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< LeftASYMH, DeviceTag > dispatcher( worklet );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
dispatcher.Invoke( leftExtend, sigIn );
break;
}
case ASYMW:
{
LeftASYMW worklet( addLen );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< LeftASYMW, DeviceTag > dispatcher( worklet );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
dispatcher.Invoke( leftExtend, sigIn );
break;
}
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
default:
{
Add throw errors
2016-08-05 23:08:25 +00:00
vtkm::cont::ErrorControlInternal("Left extension mode not supported!");
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
return 1;
}
}
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
if( !attachZeroRightLeft ) // no attach zero, or only attach on RightRight
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
{
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
// Allocate memory
if( attachZeroRightRight )
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
rightExtend.PrepareForOutput( addLen+1, DeviceTag() );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
else
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
rightExtend.PrepareForOutput( addLen, DeviceTag() );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
switch( rightExtMethod )
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
{
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
case SYMH:
{
RightSYMH worklet( sigInLen );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< RightSYMH, DeviceTag > dispatcher( worklet );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
dispatcher.Invoke( rightExtend, sigIn );
break;
}
case SYMW:
{
RightSYMW worklet( sigInLen );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< RightSYMW, DeviceTag > dispatcher( worklet );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
dispatcher.Invoke( rightExtend, sigIn );
break;
}
case ASYMH:
{
RightASYMH worklet( sigInLen );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< RightASYMH, DeviceTag > dispatcher( worklet );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
dispatcher.Invoke( rightExtend, sigIn );
break;
}
case ASYMW:
{
RightASYMW worklet( sigInLen );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< RightASYMW, DeviceTag > dispatcher( worklet );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
dispatcher.Invoke( rightExtend, sigIn );
break;
}
default:
{
vtkm::cont::ErrorControlInternal("Right extension mode not supported!");
return 1;
}
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
}
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
if( attachZeroRightRight )
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
WaveletBase::DeviceAssignZero( rightExtend, addLen, DeviceTag() );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
}
else // attachZeroRightLeft mode
{
typedef vtkm::cont::ArrayHandleConcatenate<SigInArrayType, ExtensionArrayType>
ConcatArray;
// attach a zero at the end of sigIn
ExtensionArrayType singleValArray;
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
singleValArray.PrepareForOutput(1, DeviceTag() );
WaveletBase::DeviceAssignZero( singleValArray, 0, DeviceTag() );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
ConcatArray sigInPlusOne( sigIn, singleValArray );
// allocate memory for extension
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
rightExtend.PrepareForOutput( addLen, DeviceTag() );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
switch( rightExtMethod )
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
{
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
case SYMH:
{
RightSYMH worklet( sigInLen + 1 );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< RightSYMH, DeviceTag > dispatcher( worklet );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
dispatcher.Invoke( rightExtend, sigInPlusOne );
break;
}
case SYMW:
{
RightSYMW worklet( sigInLen + 1 );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< RightSYMW, DeviceTag > dispatcher( worklet );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
dispatcher.Invoke( rightExtend, sigInPlusOne );
break;
}
case ASYMH:
{
RightASYMH worklet( sigInLen + 1 );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< RightASYMH, DeviceTag > dispatcher( worklet );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
dispatcher.Invoke( rightExtend, sigInPlusOne );
break;
}
case ASYMW:
{
RightASYMW worklet( sigInLen + 1 );
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
vtkm::worklet::DispatcherMapField< RightASYMW, DeviceTag > dispatcher( worklet );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
dispatcher.Invoke( rightExtend, sigInPlusOne );
break;
}
default:
{
vtkm::cont::ErrorControlInternal("Right extension mode not supported!");
return 1;
}
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
}
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
// make a copy of rightExtend with a zero attached to the left
half done even length filters implementation
2016-08-15 19:46:35 +00:00
ExtensionArrayType rightExtendPlusOne;
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
rightExtendPlusOne.PrepareForOutput( addLen + 1, DeviceTag() );
WaveletBase::DeviceCopyStartX( rightExtend, rightExtendPlusOne, 1, DeviceTag() );
WaveletBase::DeviceAssignZero( rightExtendPlusOne, 0, DeviceTag() );
half done even length filters implementation
2016-08-15 19:46:35 +00:00
rightExtend = rightExtendPlusOne ;
use a worklet to perform array extension
2016-08-03 20:34:02 +00:00
}
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
ArrayConcat leftOn( leftExtend, sigIn );
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
sigOut = vtkm::cont::make_ArrayHandleConcatenate( leftOn, rightExtend );
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
return 0;
}
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
// Performs one level of 1D discrete wavelet transform
// It takes care of boundary conditions, etc.
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
template< typename SignalArrayType, typename CoeffArrayType, typename DeviceTag>
take off a typedef; some formality changes
2017-02-15 19:54:53 +00:00
vtkm::Float64 DWT1D( const SignalArrayType &sigIn, // Input
CoeffArrayType &coeffOut, // Output: cA followed by cD
std::vector<vtkm::Id> &L, // Output: how many cA and cD.
DeviceTag )
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
{
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
if( GetWaveletMaxLevel( sigInLen ) < 1 )
{
Add throw errors
2016-08-05 23:08:25 +00:00
vtkm::cont::ErrorControlInternal( "Signal is too short to perform DWT!" );
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
return -1;
}
use std::vector<vtkm::Id> instead of raw arrays for bookkeeping
2016-08-08 23:59:15 +00:00
VTKM_ASSERT( L.size() == 3 );
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
L[0] = WaveletBase::GetApproxLength( sigInLen );
L[1] = WaveletBase::GetDetailLength( sigInLen );
L[2] = sigInLen;
VTKM_ASSERT( L[0] + L[1] == L[2] );
Change on class signatures of Filter and Base
2016-08-10 23:24:22 +00:00
vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
bool doSymConv = false;
Change on class signatures of Filter and Base
2016-08-10 23:24:22 +00:00
if( WaveletBase::filter.isSymmetric() )
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
{
if( ( WaveletBase::wmode == SYMW && ( filterLen % 2 != 0 ) ) ||
( WaveletBase::wmode == SYMH && ( filterLen % 2 == 0 ) ) )
doSymConv = true;
}
vtkm::Id sigConvolvedLen = L[0] + L[1]; // approx + detail coeffs
vtkm::Id addLen; // for extension
bool oddLow = true;
bool oddHigh = true;
if( filterLen % 2 != 0 )
oddLow = false;
if( doSymConv )
{
addLen = filterLen / 2;
if( sigInLen % 2 != 0 )
sigConvolvedLen += 1;
}
else
addLen = filterLen - 1;
vtkm::Id sigExtendedLen = sigInLen + 2 * addLen;
DWT worklet uses WholeArrayOut to avoid a memory copy
2016-08-02 18:09:34 +00:00
typedef typename SignalArrayType::ValueType SigInValueType;
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
typedef vtkm::cont::ArrayHandle<SigInValueType> SigInBasic;
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
typedef vtkm::cont::ArrayHandleConcatenate< SigInBasic, SignalArrayType >
ConcatType1;
typedef vtkm::cont::ArrayHandleConcatenate< ConcatType1, SigInBasic >
ConcatType2;
ConcatType2 sigInExtended;
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
this->Extend1D( sigIn, sigInExtended, addLen, WaveletBase::wmode,
WaveletBase::wmode, false, false, DeviceTag() );
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
VTKM_ASSERT( sigInExtended.GetNumberOfValues() == sigExtendedLen );
half done even length filters implementation
2016-08-15 19:46:35 +00:00
// initialize a worklet for forward transform
put filter coeffs into worklet
2016-09-08 21:49:23 +00:00
vtkm::worklet::wavelets::ForwardTransform<DeviceTag> forwardTransform
( WaveletBase::filter.GetLowDecomposeFilter(),
WaveletBase::filter.GetHighDecomposeFilter(),
filterLen, L[0], L[1], oddLow, oddHigh );
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
coeffOut.PrepareForOutput( sigExtendedLen, DeviceTag() );
put filter coeffs into worklet
2016-09-08 21:49:23 +00:00
vtkm::worklet::DispatcherMapField<vtkm::worklet::wavelets::ForwardTransform<DeviceTag>, DeviceTag>
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
dispatcher(forwardTransform);
add timers to report worklet computation time
2016-08-23 23:52:40 +00:00
// put a timer
fix a cuda memory bug
2016-09-07 23:43:10 +00:00
vtkm::cont::Timer<DeviceTag> timer;
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
dispatcher.Invoke( sigInExtended,
DWT worklet uses WholeArrayOut to avoid a memory copy
2016-08-02 18:09:34 +00:00
coeffOut );
add timers to report worklet computation time
2016-08-23 23:52:40 +00:00
vtkm::Float64 elapsedTime = timer.GetElapsedTime();
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
DWT worklet uses WholeArrayOut to avoid a memory copy
2016-08-02 18:09:34 +00:00
VTKM_ASSERT( L[0] + L[1] <= coeffOut.GetNumberOfValues() );
coeffOut.Shrink( L[0] + L[1] );
add timers to report worklet computation time
2016-08-23 23:52:40 +00:00
return elapsedTime;
}
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
// Performs one level of inverse wavelet transform
// It takes care of boundary conditions, etc.
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
template< typename CoeffArrayType, typename SignalArrayType, typename DeviceTag >
take off a typedef; some formality changes
2017-02-15 19:54:53 +00:00
vtkm::Float64 IDWT1D( const CoeffArrayType &coeffIn, // Input, cA followed by cD
std::vector<vtkm::Id> &L, // Input, how many cA and cD
SignalArrayType &sigOut, // Output
DeviceTag )
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
{
use std::vector<vtkm::Id> instead of raw arrays for bookkeeping
2016-08-08 23:59:15 +00:00
VTKM_ASSERT( L.size() == 3 );
half done even length filters implementation
2016-08-15 19:46:35 +00:00
VTKM_ASSERT( L[2] == coeffIn.GetNumberOfValues() );
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
Change on class signatures of Filter and Base
2016-08-10 23:24:22 +00:00
vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
bool doSymConv = false;
vtkm::worklet::wavelets::DWTMode cALeftMode = WaveletBase::wmode;
vtkm::worklet::wavelets::DWTMode cARightMode = WaveletBase::wmode;
vtkm::worklet::wavelets::DWTMode cDLeftMode = WaveletBase::wmode;
vtkm::worklet::wavelets::DWTMode cDRightMode = WaveletBase::wmode;
tried a new idea on transform worklets, without success. Need to discuss with Chris.
2016-09-01 19:20:22 +00:00
if( WaveletBase::filter.isSymmetric() ) // this is always true with the 1st 4 filters.
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
{
if(( WaveletBase::wmode == SYMW && (filterLen % 2 != 0) ) ||
( WaveletBase::wmode == SYMH && (filterLen % 2 == 0) ) )
{
tried a new idea on transform worklets, without success. Need to discuss with Chris.
2016-09-01 19:20:22 +00:00
doSymConv = true; // doSymConv is always true with the 1st 4 filters.
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
if( WaveletBase::wmode == SYMH )
{
cDLeftMode = ASYMH;
if( L[2] % 2 != 0 )
{
cARightMode = SYMW;
cDRightMode = ASYMW;
}
else
cDRightMode = ASYMH;
}
else
{
cDLeftMode = SYMH;
if( L[2] % 2 != 0 )
{
cARightMode = SYMW;
cDRightMode = SYMH;
}
else
cARightMode = SYMH;
}
}
}
working on even length filters; need ASYM* support in Extend1D()
2016-08-15 23:54:15 +00:00
vtkm::Id cATempLen, cDTempLen; //, reconTempLen;
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
vtkm::Id addLen = 0;
vtkm::Id cDPadLen = 0;
if( doSymConv ) // extend cA and cD
{
working on even length filters; need ASYM* support in Extend1D()
2016-08-15 23:54:15 +00:00
addLen = filterLen / 4; // addLen == 0 for Haar kernel
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
if( (L[0] > L[1]) && (WaveletBase::wmode == SYMH) )
half done even length filters implementation
2016-08-15 19:46:35 +00:00
cDPadLen = L[0];
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
cATempLen = L[0] + 2 * addLen;
cDTempLen = cATempLen; // same length
}
else // not extend cA and cD
half done even length filters implementation
2016-08-15 19:46:35 +00:00
{ // (biorthogonal kernels won't come into this case)
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
cATempLen = L[0];
cDTempLen = L[1];
}
IDWT worklet uses WholeArrayOut now
2016-08-02 22:14:03 +00:00
typedef vtkm::cont::ArrayHandleCounting< vtkm::Id > IdArrayType;
typedef vtkm::cont::ArrayHandlePermutation< IdArrayType, CoeffArrayType >
PermutArrayType;
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
// Separate cA and cD
IdArrayType approxIndices( 0, 1, L[0] );
IdArrayType detailIndices( L[0], 1, L[1] );
PermutArrayType cA( approxIndices, coeffIn );
PermutArrayType cD( detailIndices, coeffIn );
typedef typename CoeffArrayType::ValueType CoeffValueType;
typedef vtkm::cont::ArrayHandle< CoeffValueType > ExtensionArrayType;
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
typedef vtkm::cont::ArrayHandleConcatenate< ExtensionArrayType, PermutArrayType >
Concat1;
typedef vtkm::cont::ArrayHandleConcatenate< Concat1, ExtensionArrayType >
Concat2;
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
Concat2 cATemp, cDTemp;
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
if( doSymConv ) // Actually extend cA and cD
{
half done even length filters implementation
2016-08-15 19:46:35 +00:00
// first extend cA to be cATemp
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
this->Extend1D( cA, cATemp, addLen, cALeftMode, cARightMode, false, false, DeviceTag() );
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
// Then extend cD based on extension needs
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
if( cDPadLen > 0 )
{
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
// Add back the missing final cD, 0.0, before doing extension
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
this->Extend1D( cD, cDTemp, addLen, cDLeftMode, cDRightMode, true, false, DeviceTag() );
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
}
else
{
half done even length filters implementation
2016-08-15 19:46:35 +00:00
vtkm::Id cDTempLenWouldBe = L[1] + 2 * addLen;
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
if( cDTempLenWouldBe == cDTempLen )
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
{
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
this->Extend1D( cD, cDTemp, addLen, cDLeftMode, cDRightMode, false, false, DeviceTag());
working on even length filters; need ASYM* support in Extend1D()
2016-08-15 23:54:15 +00:00
}
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
else if( cDTempLenWouldBe == cDTempLen - 1 )
two new wavelet filters, HAAR and CDF8/4 supported now
2016-08-16 21:20:19 +00:00
{
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
this->Extend1D( cD, cDTemp, addLen, cDLeftMode, cDRightMode, false, true , DeviceTag());
working on even length filters; need ASYM* support in Extend1D()
2016-08-15 23:54:15 +00:00
}
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
else
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
{
ready for review, but suspect boost is messing up with smart pointers
2016-08-06 02:45:30 +00:00
vtkm::cont::ErrorControlInternal("cDTemp Length not match!");
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
return 1;
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
}
}
extend1D returns ArrayHandlePermut
2016-08-03 03:41:46 +00:00
}
fix a cuda memory bug
2016-09-07 23:43:10 +00:00
else
half done even length filters implementation
2016-08-15 19:46:35 +00:00
{
// make cATemp
ExtensionArrayType dummyArray;
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
dummyArray.PrepareForOutput(0, DeviceTag() );
half done even length filters implementation
2016-08-15 19:46:35 +00:00
Concat1 cALeftOn( dummyArray, cA );
cATemp = vtkm::cont::make_ArrayHandleConcatenate< Concat1, ExtensionArrayType >
( cALeftOn, dummyArray );
// make cDTemp
Concat1 cDLeftOn( dummyArray, cD );
cDTemp = vtkm::cont::make_ArrayHandleConcatenate< Concat1, ExtensionArrayType >
( cDLeftOn, dummyArray );
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
}
working on even length filters; need ASYM* support in Extend1D()
2016-08-15 23:54:15 +00:00
// make sure signal extension went as expected
VTKM_ASSERT( cATemp.GetNumberOfValues() == cATempLen );
VTKM_ASSERT( cDTemp.GetNumberOfValues() == cDTempLen );
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
working on even length filters; need ASYM* support in Extend1D()
2016-08-15 23:54:15 +00:00
vtkm::cont::ArrayHandleConcatenate< Concat2, Concat2>
coeffInExtended( cATemp, cDTemp );
extend2D finally works well
2016-09-05 22:53:44 +00:00
working on even length filters; need ASYM* support in Extend1D()
2016-08-15 23:54:15 +00:00
// Allocate memory for sigOut
get rid of many small memory transfers
2016-08-22 20:49:13 +00:00
sigOut.PrepareForOutput( cATempLen + cDTempLen, DeviceTag() );
fixed a memory bug, now I believe it's ready for merge
2016-08-06 07:08:29 +00:00
fix a bug, run final tests
2016-09-06 21:28:46 +00:00
vtkm::Float64 elapsedTime = 0;
working on even length filters; need ASYM* support in Extend1D()
2016-08-15 23:54:15 +00:00
if( filterLen % 2 != 0 )
{
InverseTransformOdd() keeps coeffs by itself
2016-09-08 22:28:03 +00:00
vtkm::worklet::wavelets::InverseTransformOdd<DeviceTag> inverseXformOdd
( WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
filterLen, L[0], cATempLen );
InverseTransformEven() keeps coeffs by itself
2016-09-08 22:49:20 +00:00
vtkm::worklet::DispatcherMapField<vtkm::worklet::wavelets::
InverseTransformOdd<DeviceTag>, DeviceTag> dispatcher( inverseXformOdd );
add timers to report worklet computation time
2016-08-23 23:52:40 +00:00
// use a timer
fix a cuda memory bug
2016-09-07 23:43:10 +00:00
vtkm::cont::Timer<DeviceTag> timer;
InverseTransformEven() keeps coeffs by itself
2016-09-08 22:49:20 +00:00
dispatcher.Invoke( coeffInExtended, sigOut );
add timers to report worklet computation time
2016-08-23 23:52:40 +00:00
elapsedTime = timer.GetElapsedTime();
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
}
else
{
InverseTransformEven() keeps coeffs by itself
2016-09-08 22:49:20 +00:00
vtkm::worklet::wavelets::InverseTransformEven<DeviceTag> inverseXformEven
( WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
filterLen, L[0], cATempLen, !doSymConv );
vtkm::worklet::DispatcherMapField<vtkm::worklet::wavelets::
InverseTransformEven<DeviceTag>, DeviceTag> dispatcher( inverseXformEven );
add timers to report worklet computation time
2016-08-23 23:52:40 +00:00
// use a timer
fix a cuda memory bug
2016-09-07 23:43:10 +00:00
vtkm::cont::Timer<DeviceTag> timer;
InverseTransformEven() keeps coeffs by itself
2016-09-08 22:49:20 +00:00
dispatcher.Invoke( coeffInExtended, sigOut );
add timers to report worklet computation time
2016-08-23 23:52:40 +00:00
elapsedTime = timer.GetElapsedTime();
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
}
working on even length filters; need ASYM* support in Extend1D()
2016-08-15 23:54:15 +00:00
sigOut.Shrink( L[2] );
add timers to report worklet computation time
2016-08-23 23:52:40 +00:00
return elapsedTime;
}
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
new branch on 2D wavelets
2016-08-09 23:32:13 +00:00
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
new branch on 2D wavelets
2016-08-09 23:32:13 +00:00
// Performs one level of 2D discrete wavelet transform
// It takes care of boundary conditions, etc.
// N.B.
// L[0] == L[2]
// L[1] == L[5]
// L[3] == L[7]
// L[4] == L[6]
//
// ____L[0]_______L[4]____
// | | |
// L[1] | cA | cDv | L[5]
// | (LL) | (HL) |
// | | |
// |---------------------|
// | | |
// | cDh | cDd | L[7]
// L[3] | (LH) | (HH) |
// | | |
// |__________|__________|
// L[2] L[6]
//
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
// Performs one level of 2D discrete wavelet transform on a small rectangle of input array
// The output has the same size as the small rectangle
template< typename ArrayInType, typename ArrayOutType, typename DeviceTag >
take off a typedef; some formality changes
2017-02-15 19:54:53 +00:00
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,
DeviceTag )
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
{
VTKM_ASSERT( sigDimX * sigDimY == sigIn.GetNumberOfValues() );
VTKM_ASSERT( L.size() == 10 );
L[0] = WaveletBase::GetApproxLength( sigPretendDimX ); L[2] = L[0];
L[1] = WaveletBase::GetApproxLength( sigPretendDimY ); L[5] = L[1];
L[3] = WaveletBase::GetDetailLength( sigPretendDimY ); L[7] = L[3];
L[4] = WaveletBase::GetDetailLength( sigPretendDimX ); L[6] = L[4];
L[8] = sigPretendDimX;
L[9] = sigPretendDimY;
vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
bool oddLow = true;
if( filterLen % 2 != 0 )
oddLow = false;
vtkm::Id addLen = filterLen / 2;
typedef typename ArrayInType::ValueType ValueType;
typedef vtkm::cont::ArrayHandle<ValueType> ArrayType;
typedef vtkm::worklet::wavelets::ForwardTransform2D<DeviceTag> ForwardXForm;
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
typedef typename vtkm::worklet::DispatcherMapField< ForwardXForm, DeviceTag >
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
DispatcherType;
vtkm::cont::Timer<DeviceTag> timer;
vtkm::Float64 computationTime = 0.0;
ArrayType afterX;
afterX.PrepareForOutput( sigPretendDimX * sigPretendDimY, DeviceTag() );
// First transform on rows
{
ArrayType leftExt, rightExt;
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
this->Extend2D ( sigIn,
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
sigDimX, sigDimY,
sigStartX, sigStartY,
sigPretendDimX, sigPretendDimY,
leftExt, rightExt, addLen,
WaveletBase::wmode, WaveletBase::wmode, false, false,
true, DeviceTag() ); // Extend in left-right direction
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
ForwardXForm worklet( WaveletBase::filter.GetLowDecomposeFilter(),
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
WaveletBase::filter.GetHighDecomposeFilter(),
filterLen, L[0], oddLow, true, // left-right
addLen, sigPretendDimY,
sigDimX, sigDimY,
sigStartX, sigStartY,
sigPretendDimX, sigPretendDimY,
addLen, sigPretendDimY );
DispatcherType dispatcher(worklet);
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
timer.Reset();
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
dispatcher.Invoke( leftExt, sigIn, rightExt, afterX );
computationTime += timer.GetElapsedTime();
}
// Then do transform in Y direction
{
ArrayType topExt, bottomExt;
coeffOut.PrepareForOutput( sigPretendDimX * sigPretendDimY, DeviceTag() );
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
this->Extend2D ( afterX,
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
sigPretendDimX, sigPretendDimY,
0, 0,
sigPretendDimX, sigPretendDimY,
topExt, bottomExt, addLen,
WaveletBase::wmode, WaveletBase::wmode, false, false,
false, DeviceTag() ); // Extend in top-down direction
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
ForwardXForm worklet( WaveletBase::filter.GetLowDecomposeFilter(),
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
WaveletBase::filter.GetHighDecomposeFilter(),
filterLen, L[1], oddLow, false, // top-down
sigPretendDimX, addLen,
sigPretendDimX, sigPretendDimY,
0, 0,
sigPretendDimX, sigPretendDimY,
sigPretendDimX, addLen );
DispatcherType dispatcher( worklet );
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
timer.Reset();
forward transform now avoids data copy when performing multi-level transforms
2016-10-05 22:52:07 +00:00
dispatcher.Invoke( topExt, afterX, bottomExt, coeffOut );
computationTime += timer.GetElapsedTime();
}
return computationTime;
}
sync to use gdb
2016-09-13 19:10:38 +00:00
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
// Performs one level of IDWT.
// The output array has the same dimensions as the small rectangle.
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
template< typename ArrayInType, typename ArrayOutType, typename DeviceTag >
take off a typedef; some formality changes
2017-02-15 19:54:53 +00:00
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,
DeviceTag )
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
{
VTKM_ASSERT( L.size() == 10 );
VTKM_ASSERT( inDimX * inDimY == coeffIn.GetNumberOfValues() );
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
vtkm::Id inPretendDimX = L[0] + L[4];
vtkm::Id inPretendDimY = L[1] + L[3];
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
typedef vtkm::cont::ArrayHandle<typename ArrayInType::ValueType> BasicArrayType;
typedef vtkm::worklet::wavelets::InverseTransform2D<DeviceTag> IDWT2DWorklet;
typedef vtkm::worklet::DispatcherMapField<IDWT2DWorklet, DeviceTag> Dispatcher;
sync to laptop
2016-09-22 03:25:56 +00:00
vtkm::cont::Timer<DeviceTag> timer;
vtkm::Float64 computationTime = 0.0;
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
// First inverse transform on columns
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
BasicArrayType afterY;
1 level DWT/IDWT pass all my tests
2016-09-17 02:39:42 +00:00
{
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
BasicArrayType ext1, ext2, ext3, ext4;
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
vtkm::Id extDimX = inPretendDimX;
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
vtkm::Id ext1DimY, ext2DimY, ext3DimY, ext4DimY;
style optimization
2017-02-19 04:14:20 +00:00
this->IDWTHelper2DTopDown( coeffIn,
inDimX, inDimY,
inStartX, inStartY,
inPretendDimX, inPretendDimY,
L[1], L[3],
ext1, ext2,
ext3, ext4,
ext1DimY, ext2DimY,
ext3DimY, ext4DimY,
filterLen,
wmode,
DeviceTag() );
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
afterY.PrepareForOutput( inPretendDimX * inPretendDimY, DeviceTag() );
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
IDWT2DWorklet worklet( WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
filterLen,
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
extDimX, ext1DimY, // ext1
inPretendDimX, L[1], // cA
extDimX, ext2DimY, // ext2
extDimX, ext3DimY, // ext3
inPretendDimX, L[3], // cD
extDimX, ext4DimY, // ext4
inDimX, inDimY, // coeffIn
inStartX, inStartY, // coeffIn
false ); // top-down
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
Dispatcher dispatcher( worklet );
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
timer.Reset();
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
dispatcher.Invoke( ext1, ext2, ext3, ext4, coeffIn, afterY );
sync to laptop
2016-09-22 03:25:56 +00:00
computationTime += timer.GetElapsedTime();
sync to use gdb
2016-09-13 19:10:38 +00:00
}
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
// Then inverse transform on rows
1 level DWT/IDWT pass all my tests
2016-09-17 02:39:42 +00:00
{
BasicArrayType ext1, ext2, ext3, ext4;
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
vtkm::Id extDimY = inPretendDimY;
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
vtkm::Id ext1DimX, ext2DimX, ext3DimX, ext4DimX;
style optimization
2017-02-19 04:14:20 +00:00
this->IDWTHelper2DLeftRight( afterY,
inPretendDimX, inPretendDimY,
0, 0,
inPretendDimX, inPretendDimY,
L[0], L[4],
ext1, ext2,
ext3, ext4,
ext1DimX, ext2DimX,
ext3DimX, ext4DimX,
filterLen,
wmode,
DeviceTag() );
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
sigOut.PrepareForOutput( inPretendDimX * inPretendDimY, DeviceTag() );
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
IDWT2DWorklet worklet( WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
filterLen,
ext1DimX, extDimY, // ext1
L[0], inPretendDimY, // cA
ext2DimX, extDimY, // ext2
ext3DimX, extDimY, // ext3
L[4], inPretendDimY, // cA
ext4DimX, extDimY, // ext4
inPretendDimX, inPretendDimY,
0, 0,
true ); // left-right
Dispatcher dispatcher( worklet );
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
timer.Reset();
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
dispatcher.Invoke( ext1, ext2, ext3, ext4, afterY, sigOut );
computationTime += timer.GetElapsedTime();
}
return computationTime;
}
sync to use gdb
2016-09-13 19:10:38 +00:00
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
// decides the correct extension modes for cA and cD separately,
style optimization
2017-02-19 04:14:20 +00:00
// and fill the extensions (2D matrices)
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
template< typename ArrayInType, typename ArrayOutType, typename DeviceTag >
style optimization
2017-02-19 04:14:20 +00:00
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, // output
ArrayOutType &ext3, ArrayOutType &ext4, // output
vtkm::Id &ext1DimX, vtkm::Id &ext2DimX, // output
vtkm::Id &ext3DimX, vtkm::Id &ext4DimX, // output
vtkm::Id filterLen,
DWTMode mode,
DeviceTag )
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
{
VTKM_ASSERT( inPretendDimX = cADimX + cDDimX );
// determine extension modes
DWTMode cALeft, cARight, cDLeft, cDRight;
cALeft = cARight = cDLeft = cDRight = mode;
if( mode == SYMH )
{
cDLeft = ASYMH;
if( inPretendDimX % 2 != 0 )
{
cARight = SYMW;
cDRight = ASYMW;
}
else
cDRight = ASYMH;
}
else // mode == SYMW
{
cDLeft = SYMH;
if( inPretendDimX % 2 != 0 )
{
cARight = SYMW;
cDRight = SYMH;
}
else
cARight = SYMH;
}
// determine length after extension
vtkm::Id cAExtendedDimX, cDExtendedDimX;
vtkm::Id cDPadLen = 0;
vtkm::Id addLen = filterLen / 4; // addLen == 0 for Haar kernel
if( (cADimX > cDDimX) && (mode == SYMH) )
cDPadLen = cADimX;
cAExtendedDimX = cADimX + 2 * addLen;
cDExtendedDimX = cAExtendedDimX;
// extend cA
vtkm::Id cADimY = inPretendDimY;
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
this->Extend2D ( coeffIn, inDimX, inDimY, inStartX, inStartY, cADimX, cADimY,
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
ext1, ext2, addLen, cALeft, cARight,
false, false, true, DeviceTag() );
ext1DimX = ext2DimX = addLen;
// extend cD
vtkm::Id cDDimY = inPretendDimY;
if( cDPadLen > 0 )
{
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
this->Extend2D ( coeffIn, inDimX, inDimY, inStartX + cADimX, inStartY,
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
cDDimX, cDDimY, ext3, ext4, addLen, cDLeft, cDRight,
true, false, true, DeviceTag() );
ext3DimX = addLen;
ext4DimX = addLen + 1;
}
else
{
vtkm::Id cDExtendedWouldBe = cDDimX + 2 * addLen;
if( cDExtendedWouldBe == cDExtendedDimX )
{
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
this->Extend2D ( coeffIn, inDimX, inDimY, inStartX + cADimX, inStartY,
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
cDDimX, cDDimY, ext3, ext4, addLen, cDLeft, cDRight,
false, false, true, DeviceTag());
ext3DimX = ext4DimX = addLen;
}
else if( cDExtendedWouldBe == cDExtendedDimX - 1 )
{
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
this->Extend2D ( coeffIn, inDimX, inDimY, inStartX + cADimX, inStartY,
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
cDDimX, cDDimY, ext3, ext4, addLen, cDLeft, cDRight,
false, true, true, DeviceTag());
ext3DimX = addLen;
ext4DimX = addLen + 1;
}
else
vtkm::cont::ErrorControlInternal("cDTemp Length not match!");
}
}
almost pass all tests, move to Alaska to debug one more case
2016-09-17 01:20:53 +00:00
// decides the correct extension modes for cA and cD separately,
style optimization
2017-02-19 04:14:20 +00:00
// and fill the extensions (2D matrices)
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
template< typename ArrayInType, typename ArrayOutType, typename DeviceTag >
style optimization
2017-02-19 04:14:20 +00:00
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, // output
ArrayOutType &ext3, ArrayOutType &ext4, // output
vtkm::Id &ext1DimY, vtkm::Id &ext2DimY, // output
vtkm::Id &ext3DimY, vtkm::Id &ext4DimY, // output
vtkm::Id filterLen,
DWTMode mode,
DeviceTag )
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
{
VTKM_ASSERT( inPretendDimY = cADimY + cDDimY );
// determine extension modes
debug for 3D transforms
2017-02-21 23:56:57 +00:00
DWTMode cATopMode, cADownMode, cDTopMode, cDDownMode;
cATopMode = cADownMode = cDTopMode = cDDownMode = mode;
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
if( mode == SYMH )
{
debug for 3D transforms
2017-02-21 23:56:57 +00:00
cDTopMode = ASYMH;
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
if( inPretendDimY % 2 != 0 )
{
debug for 3D transforms
2017-02-21 23:56:57 +00:00
cADownMode = SYMW;
cDDownMode = ASYMW;
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
}
else
debug for 3D transforms
2017-02-21 23:56:57 +00:00
cDDownMode = ASYMH;
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
}
else // mode == SYMW
{
debug for 3D transforms
2017-02-21 23:56:57 +00:00
cDTopMode = SYMH;
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
if( inPretendDimY % 2 != 0 )
{
debug for 3D transforms
2017-02-21 23:56:57 +00:00
cADownMode = SYMW;
cDDownMode = SYMH;
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
}
else
debug for 3D transforms
2017-02-21 23:56:57 +00:00
cADownMode = SYMH;
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
}
// determine length after extension
vtkm::Id cAExtendedDimY, cDExtendedDimY;
vtkm::Id cDPadLen = 0;
vtkm::Id addLen = filterLen / 4; // addLen == 0 for Haar kernel
if( (cADimY > cDDimY) && (mode == SYMH) )
cDPadLen = cADimY;
cAExtendedDimY = cADimY + 2 * addLen;
cDExtendedDimY = cAExtendedDimY;
// extend cA
vtkm::Id cADimX = inPretendDimX;
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
this->Extend2D ( coeffIn, inDimX, inDimY, inStartX, inStartY, cADimX, cADimY,
debug for 3D transforms
2017-02-21 23:56:57 +00:00
ext1, ext2, addLen, cATopMode, cADownMode,
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
false, false, false, DeviceTag() );
ext1DimY = ext2DimY = addLen;
// extend cD
vtkm::Id cDDimX = inPretendDimX;
if( cDPadLen > 0 )
{
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
this->Extend2D ( coeffIn, inDimX, inDimY, inStartX, inStartY + cADimY,
debug for 3D transforms
2017-02-21 23:56:57 +00:00
cDDimX, cDDimY, ext3, ext4, addLen, cDTopMode, cDDownMode,
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
true, false, false, DeviceTag() );
ext3DimY = addLen;
ext4DimY = addLen + 1;
}
else
{
vtkm::Id cDExtendedWouldBe = cDDimY + 2 * addLen;
if( cDExtendedWouldBe == cDExtendedDimY )
{
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
this->Extend2D ( coeffIn, inDimX, inDimY, inStartX, inStartY + cADimY,
debug for 3D transforms
2017-02-21 23:56:57 +00:00
cDDimX, cDDimY, ext3, ext4, addLen, cDTopMode, cDDownMode,
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
false, false, false, DeviceTag());
ext3DimY = ext4DimY = addLen;
}
else if( cDExtendedWouldBe == cDExtendedDimY - 1 )
{
2D wavelet compressor ready for Rob to review
2016-10-10 04:52:02 +00:00
this->Extend2D ( coeffIn, inDimX, inDimY, inStartX, inStartY + cADimY,
debug for 3D transforms
2017-02-21 23:56:57 +00:00
cDDimX, cDDimY, ext3, ext4, addLen, cDTopMode, cDDownMode,
new implementation of 2D transforms pass test
2016-10-10 00:07:49 +00:00
false, true, false, DeviceTag());
ext3DimY = addLen;
ext4DimY = addLen + 1;
}
else
vtkm::cont::ErrorControlInternal("cDTemp Length not match!");
}
}
debug for 3D transforms
2017-02-21 23:56:57 +00:00
// decides the correct extension modes for cA and cD separately,
// and fill the extensions (3D cubes)
template< typename ArrayInType, typename ArrayOutType, typename DeviceTag >
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, // output
ArrayOutType &ext3, ArrayOutType &ext4, // output
vtkm::Id &ext1DimX, vtkm::Id &ext2DimX, // output
vtkm::Id &ext3DimX, vtkm::Id &ext4DimX, // output
vtkm::Id filterLen,
DWTMode mode,
DeviceTag )
{
VTKM_ASSERT( inPretendDimX = cADimX + cDDimX );
// determine extension modes
DWTMode cALeftMode, cARightMode, cDLeftMode, cDRightMode;
cALeftMode = cARightMode = cDLeftMode = cDRightMode = mode;
if( mode == SYMH )
{
cDLeftMode = ASYMH;
if( inPretendDimX % 2 != 0 )
{
cARightMode = SYMW;
cDRightMode = ASYMW;
}
else
cDRightMode = ASYMH;
}
else // mode == SYMW
{
cDLeftMode = SYMH;
if( inPretendDimX % 2 != 0 )
{
cARightMode = SYMW;
cDRightMode = SYMH;
}
else
cARightMode = SYMH;
}
// determine length after extension
vtkm::Id cAExtendedDimX, cDExtendedDimX;
vtkm::Id cDPadLen = 0;
vtkm::Id addLen = filterLen / 4; // addLen == 0 for Haar kernel
if( (cADimX > cDDimX) && (mode == SYMH) )
cDPadLen = cADimX;
cAExtendedDimX = cADimX + 2 * addLen;
cDExtendedDimX = cAExtendedDimX;
// extend cA
vtkm::Id cADimY = inPretendDimY;
vtkm::Id cADimZ = inPretendDimZ;
this->Extend3DLeftRight( coeffIn,
inDimX, inDimY, inDimZ,
inStartX, inStartY, inStartZ,
cADimX, cADimY, cADimZ,
ext1, ext2,
addLen,
cALeftMode, cARightMode,
false, false,
DeviceTag() );
ext1DimX = ext2DimX = addLen;
// extend cD
vtkm::Id cDDimY = inPretendDimY;
vtkm::Id cDDimZ = inPretendDimZ;
bool pretendSigPaddedZero, padZeroAtExt2;
if( cDPadLen > 0 )
{
ext3DimX = addLen;
ext4DimX = addLen + 1;
pretendSigPaddedZero = true;
padZeroAtExt2 = false;
}
else
{
vtkm::Id cDExtendedWouldBe = cDDimX + 2 * addLen;
if( cDExtendedWouldBe == cDExtendedDimX )
{
ext3DimX = ext4DimX = addLen;
pretendSigPaddedZero = false;
padZeroAtExt2 = false;
}
else if( cDExtendedWouldBe == cDExtendedDimX - 1 )
{
ext3DimX = addLen;
ext4DimX = addLen + 1;
pretendSigPaddedZero = false;
padZeroAtExt2 = true;
}
else
vtkm::cont::ErrorControlInternal("cDTemp Length not match!");
}
this->Extend3DLeftRight( coeffIn,
inDimX, inDimY, inDimZ,
inStartX + cADimX, inStartY, inStartZ,
cDDimX, cDDimY, cDDimZ,
ext3, ext4,
addLen,
cDLeftMode, cDRightMode,
pretendSigPaddedZero, padZeroAtExt2,
DeviceTag() );
}
template< typename ArrayInType, typename ArrayOutType, typename DeviceTag >
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, // output
ArrayOutType &ext3, ArrayOutType &ext4, // output
vtkm::Id &ext1DimY, vtkm::Id &ext2DimY, // output
vtkm::Id &ext3DimY, vtkm::Id &ext4DimY, // output
vtkm::Id filterLen,
DWTMode mode,
DeviceTag )
{
VTKM_ASSERT( inPretendDimY = cADimY + cDDimY );
// determine extension modes
DWTMode cATopMode, cADownMode, cDTopMode, cDDownMode;
cATopMode = cADownMode = cDTopMode = cDDownMode = mode;
if( mode == SYMH )
{
cDTopMode = ASYMH;
if( inPretendDimY % 2 != 0 )
{
cADownMode = SYMW;
cDDownMode = ASYMW;
}
else
cDDownMode = ASYMH;
}
else // mode == SYMW
{
cDTopMode = SYMH;
if( inPretendDimY % 2 != 0 )
{
cADownMode = SYMW;
cDDownMode = SYMH;
}
else
cADownMode = SYMH;
}
// determine length after extension
vtkm::Id cAExtendedDimY, cDExtendedDimY;
vtkm::Id cDPadLen = 0;
vtkm::Id addLen = filterLen / 4; // addLen == 0 for Haar kernel
if( (cADimY > cDDimY) && (mode == SYMH) )
cDPadLen = cADimY;
cAExtendedDimY = cADimY + 2 * addLen;
cDExtendedDimY = cAExtendedDimY;
// extend cA
vtkm::Id cADimX = inPretendDimX;
vtkm::Id cADimZ = inPretendDimZ;
this->Extend3DTopDown( coeffIn,
inDimX, inDimY, inDimZ,
inStartX, inStartY, inStartZ,
cADimX, cADimY, cADimZ,
ext1, ext2,
addLen,
cATopMode, cADownMode,
false, false,
DeviceTag() );
ext1DimY = ext2DimY = addLen;
// extend cD
vtkm::Id cDDimX = inPretendDimX;
vtkm::Id cDDimZ = inPretendDimZ;
bool pretendSigPaddedZero, padZeroAtExt2;
if( cDPadLen > 0 )
{
ext3DimY = addLen;
ext4DimY = addLen + 1;
pretendSigPaddedZero = true;
padZeroAtExt2 = false;
}
else
{
vtkm::Id cDExtendedWouldBe = cDDimY + 2 * addLen;
if( cDExtendedWouldBe == cDExtendedDimY )
{
ext3DimY = ext4DimY = addLen;
pretendSigPaddedZero = false;
padZeroAtExt2 = false;
}
else if( cDExtendedWouldBe == cDExtendedDimY - 1 )
{
ext3DimY = addLen;
ext4DimY = addLen + 1;
pretendSigPaddedZero = false;
padZeroAtExt2 = true;
}
else
vtkm::cont::ErrorControlInternal("cDTemp Length not match!");
}
this->Extend3DTopDown( coeffIn,
inDimX, inDimY, inDimZ,
inStartX, inStartY + cADimY, inStartZ,
cDDimX, cDDimY, cDDimZ,
ext3, ext4,
addLen,
cDTopMode, cDDownMode,
pretendSigPaddedZero, padZeroAtExt2,
DeviceTag() );
}
template< typename ArrayInType, typename ArrayOutType, typename DeviceTag >
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, // output
ArrayOutType &ext3, ArrayOutType &ext4, // output
vtkm::Id &ext1DimZ, vtkm::Id &ext2DimZ, // output
vtkm::Id &ext3DimZ, vtkm::Id &ext4DimZ, // output
vtkm::Id filterLen,
DWTMode mode,
DeviceTag )
{
VTKM_ASSERT( inPretendDimZ = cADimZ + cDDimZ );
// determine extension modes
DWTMode cAFrontMode, cABackMode, cDFrontMode, cDBackMode;
cAFrontMode = cABackMode = cDFrontMode = cDBackMode = mode;
if( mode == SYMH )
{
cDFrontMode = ASYMH;
if( inPretendDimZ % 2 != 0 )
{
cABackMode = SYMW;
cDBackMode = ASYMW;
}
else
cDBackMode = ASYMH;
}
else // mode == SYMW
{
cDFrontMode = SYMH;
if( inPretendDimZ % 2 != 0 )
{
cABackMode = SYMW;
cDBackMode = SYMH;
}
else
cABackMode = SYMH;
}
// determine length after extension
vtkm::Id cAExtendedDimZ, cDExtendedDimZ;
vtkm::Id cDPadLen = 0;
vtkm::Id addLen = filterLen / 4; // addLen == 0 for Haar kernel
if( (cADimZ > cDDimZ) && (mode == SYMH) )
cDPadLen = cADimZ;
cAExtendedDimZ = cADimZ + 2 * addLen;
cDExtendedDimZ = cAExtendedDimZ;
// extend cA
vtkm::Id cADimX = inPretendDimX;
vtkm::Id cADimY = inPretendDimY;
this->Extend3DFrontBack( coeffIn,
inDimX, inDimY, inDimZ,
inStartX, inStartY, inStartZ,
cADimX, cADimY, cADimZ,
ext1, ext2,
addLen,
cAFrontMode, cABackMode,
false, false,
DeviceTag() );
ext1DimZ = ext2DimZ = addLen;
// extend cD
vtkm::Id cDDimX = inPretendDimX;
vtkm::Id cDDimY = inPretendDimY;
bool pretendSigPaddedZero, padZeroAtExt2;
if( cDPadLen > 0 )
{
ext3DimZ = addLen;
ext4DimZ = addLen + 1;
pretendSigPaddedZero = true;
padZeroAtExt2 = false;
}
else
{
vtkm::Id cDExtendedWouldBe = cDDimZ + 2 * addLen;
if( cDExtendedWouldBe == cDExtendedDimZ )
{
ext3DimZ = ext4DimZ = addLen;
pretendSigPaddedZero = false;
padZeroAtExt2 = false;
}
else if( cDExtendedWouldBe == cDExtendedDimZ - 1 )
{
ext3DimZ = addLen;
ext4DimZ = addLen + 1;
pretendSigPaddedZero = false;
padZeroAtExt2 = true;
}
else
vtkm::cont::ErrorControlInternal("cDTemp Length not match!");
}
this->Extend3DFrontBack( coeffIn,
inDimX, inDimY, inDimZ,
inStartX, inStartY, inStartZ + cADimZ,
cDDimX, cDDimY, cDDimZ,
ext3, ext4,
addLen,
cDFrontMode, cDBackMode,
pretendSigPaddedZero, padZeroAtExt2,
DeviceTag() );
}
};
move wavelet stuff back to worklet folder
2016-07-27 20:28:38 +00:00
} // namespace wavelets
} // namespace worklet
} // namespace vtkm
#endif
Reference in New Issue Copy Permalink