Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-16 17:22:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge topic 'support_3DWavelet_compression'

Browse Source 
fecc2e97 code style improvements
7782ff47 Take off a few debug statements; replaced tabs with spaces
3e4095f4 WaveletDWT.h
c4853885 eliminate another warning
d4deced5 get rid of some asserts used for debugging
98cf7f8d got rid of unused variable warnings
33de20da fix print type warnings
7d123455 remove my own note
...

Acked-by: Kitware Robot <kwrobot@kitware.com>
Acked-by: Robert Maynard <robert.maynard@kitware.com>
Merge-request: !765
This commit is contained in:
Robert Maynard 2017-05-22 18:00:50 +00:00 committed by Kitware Robot
commit 206556916d
6 changed files with 3849 additions and 478 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

296
vtkm/worklet/WaveletCompressor.h
View File

@ -34,7 +34,6 @@ public:
WaveletCompressor( wavelets::WaveletName name ) : WaveletDWT( name ) {}
// Multi-level 1D wavelet decomposition
template< typename SignalArrayType, typename CoeffArrayType, typename DeviceTag >
VTKM_CONT
@ -155,16 +154,191 @@ public:
// Multi-level 3D wavelet decomposition
template< typename InArrayType, typename OutArrayType, typename DeviceTag>
VTKM_CONT
vtkm::Float64 WaveDecompose3D( InArrayType &sigIn, // Input
vtkm::Id nLevels, // n levels of DWT
vtkm::Id inX,
vtkm::Id inY,
vtkm::Id inZ,
OutArrayType &coeffOut,
bool discardSigIn, // can we discard sigIn on devices?
DeviceTag )
{
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
VTKM_ASSERT( inX * inY * inZ == sigInLen );
if( nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel( inX ) ||
nLevels > WaveletBase::GetWaveletMaxLevel( inY ) ||
nLevels > WaveletBase::GetWaveletMaxLevel( inZ ) )
{
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
}
if( nLevels == 0 ) // 0 levels means no transform
{
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
return 0;
}
vtkm::Id currentLenX = inX;
vtkm::Id currentLenY = inY;
vtkm::Id currentLenZ = inZ;
std::vector<vtkm::Id> L3d(27, 0);
typedef typename OutArrayType::ValueType OutValueType;
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
// First level transform writes to the output array
vtkm::Float64 computationTime = WaveletDWT::DWT3D(
sigIn,
inX, inY, inZ,
0, 0, 0,
currentLenX, currentLenY, currentLenZ,
coeffOut,
discardSigIn,
DeviceTag() );
// Successor transforms writes to a temporary array
for( vtkm::Id i = nLevels-1; i > 0; i-- )
{
currentLenX = WaveletBase::GetApproxLength( currentLenX );
currentLenY = WaveletBase::GetApproxLength( currentLenY );
currentLenZ = WaveletBase::GetApproxLength( currentLenZ );
OutBasicArray tempOutput;
computationTime += WaveletDWT::DWT3D(
coeffOut,
inX, inY, inZ,
0, 0, 0,
currentLenX, currentLenY, currentLenZ,
tempOutput,
false,
DeviceTag() );
// copy results to coeffOut
WaveletBase::DeviceCubeCopyTo( tempOutput,
currentLenX, currentLenY, currentLenZ,
coeffOut,
inX, inY, inZ,
0, 0, 0,
DeviceTag() );
}
return computationTime;
}
// Multi-level 3D wavelet reconstruction
template< typename InArrayType, typename OutArrayType, typename DeviceTag>
VTKM_CONT
vtkm::Float64 WaveReconstruct3D(
InArrayType &arrIn, // Input
vtkm::Id nLevels, // n levels of DWT
vtkm::Id inX, vtkm::Id inY, vtkm::Id inZ,
OutArrayType &arrOut,
bool discardArrIn, // can we discard input for more memory?
DeviceTag )
{
vtkm::Id arrInLen = arrIn.GetNumberOfValues();
VTKM_ASSERT( inX * inY * inZ == arrInLen );
if( nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel( inX ) ||
nLevels > WaveletBase::GetWaveletMaxLevel( inY ) ||
nLevels > WaveletBase::GetWaveletMaxLevel( inZ ) )
{
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
}
typedef typename OutArrayType::ValueType OutValueType;
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
vtkm::Float64 computationTime = 0.0;
OutBasicArray outBuffer;
if( nLevels == 0 ) // 0 levels means no transform
{
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, arrOut );
return 0;
}
else if ( discardArrIn )
{
outBuffer = arrIn;
}
else
{
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, outBuffer );
}
std::vector<vtkm::Id> L;
this->ComputeL3( inX, inY, inZ, nLevels, L );
std::vector<vtkm::Id> L3d(27, 0);
// All transforms but the last level operate on temporary arrays
for( size_t i = 0; i < 24; i++ )
{
L3d[i] = L[i];
}
for( size_t i = 1; i < static_cast<size_t>(nLevels); i++ )
{
L3d[24] = L3d[0] + L3d[12]; // Total X dim; this is always true for Biorthogonal wavelets
L3d[25] = L3d[1] + L3d[7]; // Total Y dim
L3d[26] = L3d[2] + L3d[5]; // Total Z dim
OutBasicArray tempOutput;
// IDWT
computationTime += WaveletDWT::IDWT3D( outBuffer,
inX, inY, inZ,
0, 0, 0,
L3d,
tempOutput,
false,
DeviceTag() );
// copy back reconstructed block
WaveletBase::DeviceCubeCopyTo( tempOutput,
L3d[24], L3d[25], L3d[26],
outBuffer,
inX, inY, inZ,
0, 0, 0,
DeviceTag() );
// update L3d array
L3d[0] = L3d[24];
L3d[1] = L3d[25];
L3d[2] = L3d[26];
for( size_t j = 3; j < 24; j++ )
{
L3d[j] = L[ 21 * i + j ];
}
}
// The last transform outputs to the final output
L3d[24] = L3d[0] + L3d[12];
L3d[25] = L3d[1] + L3d[7];
L3d[26] = L3d[2] + L3d[5];
computationTime += WaveletDWT::IDWT3D( outBuffer,
inX, inY, inZ,
0, 0, 0,
L3d,
arrOut,
true,
DeviceTag() );
return computationTime;
}
// Multi-level 2D wavelet decomposition
template< typename InArrayType, typename OutArrayType, typename DeviceTag>
VTKM_CONT
FLOAT_64 WaveDecompose2D( const InArrayType &sigIn, // Input
vtkm::Id nLevels, // n levels of DWT
vtkm::Id inX, // Input X dim
vtkm::Id inY, // Input Y dim
OutArrayType &coeffOut,
std::vector<vtkm::Id> &L,
DeviceTag )
vtkm::Float64 WaveDecompose2D( const InArrayType &sigIn, // Input
vtkm::Id nLevels, // n levels of DWT
vtkm::Id inX, // Input X dim
vtkm::Id inY, // Input Y dim
OutArrayType &coeffOut,
std::vector<vtkm::Id> &L,
DeviceTag )
{
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
VTKM_ASSERT( inX * inY == sigInLen );
@ -230,13 +404,13 @@ public:
// Multi-level 2D wavelet reconstruction
template< typename InArrayType, typename OutArrayType, typename DeviceTag>
VTKM_CONT
FLOAT_64 WaveReconstruct2D( const InArrayType &arrIn, // Input
vtkm::Id nLevels, // n levels of DWT
vtkm::Id inX, // Input X dim
vtkm::Id inY, // Input Y dim
OutArrayType &arrOut,
std::vector<vtkm::Id> &L,
DeviceTag )
vtkm::Float64 WaveReconstruct2D( const InArrayType &arrIn, // Input
vtkm::Id nLevels, // n levels of DWT
vtkm::Id inX, // Input X dim
vtkm::Id inY, // Input Y dim
OutArrayType &arrOut,
std::vector<vtkm::Id> &L,
DeviceTag )
{
vtkm::Id arrInLen = arrIn.GetNumberOfValues();
VTKM_ASSERT( inX * inY == arrInLen );
@ -256,7 +430,9 @@ public:
return 0;
}
else
{
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, outBuffer );
}
VTKM_ASSERT( vtkm::Id(L.size()) == 6 * nLevels + 4 );
@ -315,24 +491,27 @@ public:
vtkm::Float64 ratio,
DeviceTag )
{
if( ratio > 1 )
if( ratio > 1.0 )
{
vtkm::Id coeffLen = coeffIn.GetNumberOfValues();
typedef typename CoeffArrayType::ValueType ValueType;
typedef vtkm::cont::ArrayHandle< ValueType > CoeffArrayBasic;
CoeffArrayBasic sortedArray;
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( coeffIn, sortedArray );
WaveletBase::DeviceSort( sortedArray, DeviceTag() );
vtkm::Id n = coeffLen -
static_cast<vtkm::Id>( static_cast<vtkm::Float64>(coeffLen)/ratio );
vtkm::Float64 nthVal = static_cast<vtkm::Float64>
(sortedArray.GetPortalConstControl().Get(n));
if( nthVal < 0.0 )
{
nthVal *= -1.0;
}
typedef vtkm::worklet::wavelets::ThresholdWorklet ThresholdType;
ThresholdType tw( nthVal );
vtkm::worklet::DispatcherMapField< ThresholdType, DeviceTag > dispatcher( tw );
ThresholdType thresholdWorklet( nthVal );
vtkm::worklet::DispatcherMapField< ThresholdType, DeviceTag > dispatcher( thresholdWorklet );
dispatcher.Invoke( coeffIn );
}
@ -402,7 +581,7 @@ public:
// Compute the book keeping array L for 1D wavelet decomposition
// Compute the book keeping array L for 1D DWT
void ComputeL( vtkm::Id sigInLen,
vtkm::Id nLev,
std::vector<vtkm::Id> &L )
@ -420,7 +599,7 @@ public:
// Compute the book keeping array L for 2D wavelet decomposition
// Compute the book keeping array L for 2D DWT
void ComputeL2( vtkm::Id inX,
vtkm::Id inY,
vtkm::Id nLev,
@ -428,10 +607,10 @@ public:
{
size_t nLevels = static_cast<size_t>( nLev );
L.resize( nLevels*6 + 4 );
L[ nLevels*6 + 0 ] = inX;
L[ nLevels*6 + 1 ] = inY;
L[ nLevels*6 + 2 ] = inX;
L[ nLevels*6 + 3 ] = inY;
L[ nLevels*6 ] = inX;
L[ nLevels*6 + 1 ] = inY;
L[ nLevels*6 + 2 ] = inX;
L[ nLevels*6 + 3 ] = inY;
for( size_t i = nLevels; i > 0; i-- )
{
@ -455,13 +634,78 @@ public:
// Compute the bookkeeping array L for 3D DWT
void ComputeL3( vtkm::Id inX,
vtkm::Id inY,
vtkm::Id inZ,
vtkm::Id nLev,
std::vector<vtkm::Id> &L )
{
size_t n = static_cast<size_t>( nLev );
L.resize( n * 21 + 6 );
L[ n * 21 + 0 ] = inX;
L[ n * 21 + 1 ] = inY;
L[ n * 21 + 2 ] = inZ;
L[ n * 21 + 3 ] = inX;
L[ n * 21 + 4 ] = inY;
L[ n * 21 + 5 ] = inZ;
for( size_t i = n; i > 0; i-- )
{
// cLLL
L[ i * 21 - 21 ] = WaveletBase::GetApproxLength( L[ i * 21 + 0 ] );
L[ i * 21 - 20 ] = WaveletBase::GetApproxLength( L[ i * 21 + 1 ] );
L[ i * 21 - 19 ] = WaveletBase::GetApproxLength( L[ i * 21 + 2 ] );
// cLLH
L[ i * 21 - 18 ] = L[ i * 21 - 21 ];
L[ i * 21 - 17 ] = L[ i * 21 - 20 ];
L[ i * 21 - 16 ] = WaveletBase::GetDetailLength( L[ i * 21 + 2 ] );
// cLHL
L[ i * 21 - 15 ] = L[ i * 21 - 21 ];
L[ i * 21 - 14 ] = WaveletBase::GetDetailLength( L[ i * 21 + 1 ] );
L[ i * 21 - 13 ] = L[ i * 21 - 19 ];
// cLHH
L[ i * 21 - 12 ] = L[ i * 21 - 21 ];
L[ i * 21 - 11 ] = L[ i * 21 - 14 ];
L[ i * 21 - 10 ] = L[ i * 21 - 16 ];
// cHLL
L[ i * 21 - 9 ] = WaveletBase::GetDetailLength( L[ i * 21 + 0 ] );
L[ i * 21 - 8 ] = L[ i * 21 - 20 ];
L[ i * 21 - 7 ] = L[ i * 21 - 19 ];
// cHLH
L[ i * 21 - 6 ] = L[ i * 21 - 9 ];
L[ i * 21 - 5 ] = L[ i * 21 - 20 ];
L[ i * 21 - 3 ] = L[ i * 21 - 16 ];
// cHHL
L[ i * 21 - 3 ] = L[ i * 21 - 9 ];
L[ i * 21 - 2 ] = L[ i * 21 - 14 ];
L[ i * 21 - 1 ] = L[ i * 21 - 19 ];
// cHHH - overwrites previous value
L[ i * 21 + 0 ] = L[ i * 21 - 9 ];
L[ i * 21 + 1 ] = L[ i * 21 - 14 ];
L[ i * 21 + 2 ] = L[ i * 21 - 16 ];
}
}
// Compute the length of coefficients for 1D transforms
vtkm::Id ComputeCoeffLength( std::vector<vtkm::Id> &L,
vtkm::Id nLevels )
{
vtkm::Id sum = L[0]; // 1st level cA
for( size_t i = 1; i <= size_t(nLevels); i++ )
{
sum += L[i];
}
return sum;
}
// Compute the length of coefficients for 2D transforms
@ -486,7 +730,9 @@ public:
{
cALen = WaveletBase::GetApproxLength( cALen );
if( cALen == 0 )
{
return cALen;
}
}
return cALen;

198
vtkm/worklet/testing/UnitTestWaveletCompressor.cxx
View File

@ -34,20 +34,6 @@ namespace worklet
namespace wavelets
{
class SineWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldInOut<>);
typedef void ExecutionSignature(_1, WorkIndex);
template<typename T>
VTKM_EXEC
void operator()(T& x, const vtkm::Id& workIdx) const
{
x = vtkm::Sin(vtkm::Float64(workIdx) / 100.0) * 100.0;
}
};
class GaussianWorklet2D : public vtkm::worklet::WorkletMapField
{
public:
@ -97,17 +83,59 @@ private: // see wikipedia page
vtkm::Float64 sigmaX2, sigmaY2; // 2 * sigma * sigma
};
}
}
}
template< typename ArrayType >
void FillArray( ArrayType& array )
template< typename T>
class GaussianWorklet3D : public vtkm::worklet::WorkletMapField
{
typedef vtkm::worklet::wavelets::SineWorklet SineWorklet;
SineWorklet worklet;
vtkm::worklet::DispatcherMapField< SineWorklet > dispatcher( worklet );
dispatcher.Invoke( array );
public:
typedef void ControlSignature(FieldInOut<>);
typedef void ExecutionSignature(_1, WorkIndex);
VTKM_EXEC
GaussianWorklet3D( vtkm::Id dx, vtkm::Id dy, vtkm::Id dz )
: dimX( dx ), dimY( dy ), dimZ( dz )
{
amp = (T)20.0;
sigmaX = (T)dimX / (T)4.0; sigmaX2 = sigmaX * sigmaX * (T)2.0;
sigmaY = (T)dimY / (T)4.0; sigmaY2 = sigmaY * sigmaY * (T)2.0;
sigmaZ = (T)dimZ / (T)4.0; sigmaZ2 = sigmaZ * sigmaZ * (T)2.0;
}
VTKM_EXEC
void Sig1Dto3D( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y, vtkm::Id &z ) const
{
z = idx / (dimX * dimY);
y = (idx - z * dimX * dimY) / dimX;
x = idx % dimX;
}
VTKM_EXEC
T GetGaussian( T x, T y, T z ) const
{
x -= (T)dimX / (T)2.0; // translate to center at (0, 0, 0)
y -= (T)dimY / (T)2.0;
z -= (T)dimZ / (T)2.0;
T power = x*x/sigmaX2 + y*y/sigmaY2 + z*z/sigmaZ2;
return vtkm::Exp( power * (T)-1.0 ) * amp;
}
VTKM_EXEC
void operator()(T& val, const vtkm::Id& workIdx) const
{
vtkm::Id x, y, z;
Sig1Dto3D( workIdx, x, y, z );
val = GetGaussian( (T)x, (T)y, (T)z );
}
private:
const vtkm::Id dimX, dimY, dimZ; // extent
T amp; // amplitude
T sigmaX, sigmaY, sigmaZ; // spread
T sigmaX2, sigmaY2, sigmaZ2; // sigma * sigma * 2
};
}
}
}
template< typename ArrayType >
@ -122,12 +150,98 @@ void FillArray2D( ArrayType& array, vtkm::Id dimX, vtkm::Id dimY )
vtkm::worklet::DispatcherMapField< WorkletType > dispatcher( worklet );
dispatcher.Invoke( array );
}
void TestDecomposeReconstruct2D()
template< typename ArrayType >
void FillArray3D( ArrayType& array, vtkm::Id dimX, vtkm::Id dimY, vtkm::Id dimZ )
{
std::cout << "Testing 2D wavelet compressor on a 1000x1000 square: " << std::endl;
typedef vtkm::worklet::wavelets::GaussianWorklet3D< typename ArrayType::ValueType> WorkletType;
WorkletType worklet( dimX, dimY, dimZ );
vtkm::worklet::DispatcherMapField< WorkletType > dispatcher( worklet );
dispatcher.Invoke( array );
}
void TestDecomposeReconstruct3D( vtkm::Float64 cratio )
{
vtkm::Id sigX = 99;
vtkm::Id sigY = 99;
vtkm::Id sigZ = 99;
vtkm::Id sigLen = sigX * sigY * sigZ;
std::cout << "Testing 3D wavelet compressor on a (99x99x99) cube..." << std::endl;
// make input data array handle
vtkm::cont::ArrayHandle<vtkm::Float32> inputArray;
inputArray.PrepareForOutput( sigLen, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
FillArray3D( inputArray, sigX, sigY, sigZ );
vtkm::cont::ArrayHandle<vtkm::Float32> outputArray;
// Use a WaveletCompressor
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::BIOR4_4;
if( wname == vtkm::worklet::wavelets::BIOR1_1 )
std::cout << "Using wavelet kernel = Bior1.1 (HAAR)" << std::endl;
else if( wname == vtkm::worklet::wavelets::BIOR2_2 )
std::cout << "Using wavelet kernel = Bior2.2 (CDF 5/3)" << std::endl;
else if( wname == vtkm::worklet::wavelets::BIOR3_3 )
std::cout << "Using wavelet kernel = Bior3.3 (CDF 8/4)" << std::endl;
else if( wname == vtkm::worklet::wavelets::BIOR4_4 )
std::cout << "Using wavelet kernel = Bior4.4 (CDF 9/7)" << std::endl;
vtkm::worklet::WaveletCompressor compressor( wname );
vtkm::Id XMaxLevel = compressor.GetWaveletMaxLevel( sigX );
vtkm::Id YMaxLevel = compressor.GetWaveletMaxLevel( sigY );
vtkm::Id ZMaxLevel = compressor.GetWaveletMaxLevel( sigZ );
vtkm::Id nLevels = vtkm::Min( vtkm::Min(XMaxLevel, YMaxLevel), ZMaxLevel );
std::cout << "Decomposition levels = " << nLevels << std::endl;
vtkm::Float64 computationTime = 0.0;
vtkm::Float64 elapsedTime1, elapsedTime2, elapsedTime3;
// Decompose
vtkm::cont::Timer<> timer;
computationTime =
compressor.WaveDecompose3D( inputArray, nLevels, sigX, sigY, sigZ, outputArray,
false, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
elapsedTime1 = timer.GetElapsedTime();
std::cout << "Decompose time = " << elapsedTime1 << std::endl;
std::cout << " ->computation time = " << computationTime << std::endl;
// Squash small coefficients
timer.Reset();
compressor.SquashCoefficients( outputArray, cratio, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
elapsedTime2 = timer.GetElapsedTime();
std::cout << "Squash time = " << elapsedTime2 << std::endl;
// Reconstruct
vtkm::cont::ArrayHandle<vtkm::Float32> reconstructArray;
timer.Reset();
computationTime =
compressor.WaveReconstruct3D( outputArray, nLevels, sigX, sigY, sigZ, reconstructArray,
false, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
elapsedTime3 = timer.GetElapsedTime();
std::cout << "Reconstruction time = " << elapsedTime3 << std::endl;
std::cout << " ->computation time = " << computationTime << std::endl;
std::cout << "Total time = "
<< (elapsedTime1 + elapsedTime2 + elapsedTime3) << std::endl;
outputArray.ReleaseResources();
compressor.EvaluateReconstruction( inputArray, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
timer.Reset();
for( vtkm::Id i = 0; i < reconstructArray.GetNumberOfValues(); i++ )
{
VTKM_TEST_ASSERT( test_equal( reconstructArray.GetPortalConstControl().Get(i),
inputArray.GetPortalConstControl().Get(i) ),
"WaveletCompressor 3D failed..." );
}
elapsedTime1 = timer.GetElapsedTime();
std::cout << "Verification time = " << elapsedTime1 << std::endl;
}
void TestDecomposeReconstruct2D( vtkm::Float64 cratio )
{
std::cout << "Testing 2D wavelet compressor on a (1000x1000) square... " << std::endl;
vtkm::Id sigX = 1000;
vtkm::Id sigY = 1000;
vtkm::Id sigLen = sigX * sigY;
@ -163,7 +277,6 @@ void TestDecomposeReconstruct2D()
// Squash small coefficients
timer.Reset();
vtkm::Float64 cratio = 1.0; // X:1 compression, where X >= 1
compressor.SquashCoefficients( outputArray, cratio, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
elapsedTime2 = timer.GetElapsedTime();
std::cout << "Squash time = " << elapsedTime2 << std::endl;
@ -196,11 +309,10 @@ void TestDecomposeReconstruct2D()
}
void TestDecomposeReconstruct1D()
void TestDecomposeReconstruct1D( vtkm::Float64 cratio )
{
std::cout << "Testing 1D wavelet compressor on a 1 million sized array " << std::endl;
vtkm::Id million = 1000000;
vtkm::Id sigLen = million * 1;
std::cout << "Testing 1D wavelet compressor on a 1 million sized array... " << std::endl;
vtkm::Id sigLen = 1000000;
// make input data array handle
std::vector<vtkm::Float64> tmpVector;
@ -218,7 +330,7 @@ void TestDecomposeReconstruct1D()
std::cout << "Wavelet kernel = CDF 9/7" << std::endl;
vtkm::worklet::WaveletCompressor compressor( wname );
// User maximum decompose levels, and no compression
// User maximum decompose levels
vtkm::Id maxLevel = compressor.GetWaveletMaxLevel( sigLen );
vtkm::Id nLevels = maxLevel;
std::cout << "Decomposition levels = " << nLevels << std::endl;
@ -231,6 +343,12 @@ void TestDecomposeReconstruct1D()
vtkm::Float64 elapsedTime = timer.GetElapsedTime();
std::cout << "Decompose time = " << elapsedTime << std::endl;
// Squash small coefficients
timer.Reset();
compressor.SquashCoefficients( outputArray, cratio, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
elapsedTime = timer.GetElapsedTime();
std::cout << "Squash time = " << elapsedTime << std::endl;
// Reconstruct
vtkm::cont::ArrayHandle<vtkm::Float64> reconstructArray;
@ -254,9 +372,15 @@ void TestDecomposeReconstruct1D()
void TestWaveletCompressor()
{
TestDecomposeReconstruct1D();
vtkm::Float64 cratio = 2.0; // X:1 compression, where X >= 1
std::cout << "Compression ratio = " << cratio << ":1 ";
std::cout << "(Reconstruction using higher compression ratios may result in failure in verification)" << std::endl;
TestDecomposeReconstruct1D( cratio );
std::cout << std::endl;
TestDecomposeReconstruct2D();
TestDecomposeReconstruct2D( cratio );
std::cout << std::endl;
TestDecomposeReconstruct3D( cratio );
}
int UnitTestWaveletCompressor(int, char *[])

94
vtkm/worklet/wavelets/WaveletBase.h
View File

@ -60,9 +60,13 @@ public:
vtkm::Id GetApproxLength( vtkm::Id sigInLen )
{
if (sigInLen % 2 != 0)
{
return((sigInLen+1) / 2);
}
else
{
return((sigInLen) / 2);
}
}
@ -70,9 +74,13 @@ public:
vtkm::Id GetDetailLength( vtkm::Id sigInLen )
{
if (sigInLen % 2 != 0)
{
return((sigInLen-1) / 2);
}
else
{
return((sigInLen) / 2);
}
}
@ -158,6 +166,51 @@ public:
// Assign zeros to a plane that's perpendicular to the X axis (Left-Right direction)
template< typename ArrayType, typename DeviceTag >
void DeviceAssignZero3DPlaneX( ArrayType &array, // input array
vtkm::Id dimX, vtkm::Id dimY, vtkm::Id dimZ, // dims of input
vtkm::Id zeroX, // X idx to set zero
DeviceTag )
{
typedef vtkm::worklet::wavelets::AssignZero3DWorklet AssignZero3DType;
AssignZero3DType zeroWorklet( dimX, dimY, dimZ, zeroX, -1, -1 );
vtkm::worklet::DispatcherMapField< AssignZero3DType, DeviceTag > dispatcher( zeroWorklet );
dispatcher.Invoke( array );
}
// Assign zeros to a plane that's perpendicular to the Y axis (Top-Down direction)
template< typename ArrayType, typename DeviceTag >
void DeviceAssignZero3DPlaneY( ArrayType &array, // input array
vtkm::Id dimX, vtkm::Id dimY, vtkm::Id dimZ, // dims of input
vtkm::Id zeroY, // Y idx to set zero
DeviceTag )
{
typedef vtkm::worklet::wavelets::AssignZero3DWorklet AssignZero3DType;
AssignZero3DType zeroWorklet( dimX, dimY, dimZ, -1, zeroY, -1 );
vtkm::worklet::DispatcherMapField< AssignZero3DType, DeviceTag > dispatcher( zeroWorklet );
dispatcher.Invoke( array );
}
// Assign zeros to a plane that's perpendicular to the Z axis (Front-Back direction)
template< typename ArrayType, typename DeviceTag >
void DeviceAssignZero3DPlaneZ( ArrayType &array, // input array
vtkm::Id dimX, vtkm::Id dimY, vtkm::Id dimZ, // dims of input
vtkm::Id zeroZ, // Y idx to set zero
DeviceTag )
{
typedef vtkm::worklet::wavelets::AssignZero3DWorklet AssignZero3DType;
AssignZero3DType zeroWorklet( dimX, dimY, dimZ, -1, -1, zeroZ );
vtkm::worklet::DispatcherMapField< AssignZero3DType, DeviceTag > dispatcher( zeroWorklet );
dispatcher.Invoke( array );
}
// Sort by the absolute value on device
struct SortLessAbsFunctor
{
@ -182,7 +235,7 @@ public:
typename ArrayType::ValueType DeviceSum( const ArrayType &array, DeviceTag )
{
return vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Reduce
( array, 0.0 );
( array, static_cast<typename ArrayType::ValueType>(0.0) );
}
@ -287,27 +340,24 @@ public:
// Fill a small rectangle from a portion of a big rectangle
template< typename SmallArrayType, typename BigArrayType, typename DeviceTag >
void DeviceRectangleCopyFrom( SmallArrayType &smallRect,
vtkm::Id smallX,
vtkm::Id smallY,
const BigArrayType &bigRect,
vtkm::Id bigX,
vtkm::Id bigY,
vtkm::Id startX,
vtkm::Id startY,
DeviceTag )
// Copy a small cube to a big cube
template< typename SmallArrayType, typename BigArrayType, typename DeviceTag>
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,
DeviceTag )
{
smallRect.PrepareForOutput( smallX*smallY, DeviceTag() );
typedef vtkm::worklet::wavelets::RectangleCopyFrom CopyFromWorklet;
CopyFromWorklet cpFrom( smallX, smallY, bigX, bigY, startX, startY );
vtkm::worklet::DispatcherMapField< CopyFromWorklet, DeviceTag > dispatcherFrom( cpFrom );
dispatcherFrom.Invoke( smallRect, bigRect );
typedef vtkm::worklet::wavelets::CubeCopyTo CopyToWorklet;
CopyToWorklet cp( smallX, smallY, smallZ, bigX, bigY, bigZ, startX, startY, startZ );
vtkm::worklet::DispatcherMapField< CopyToWorklet, DeviceTag > dispatcher( cp );
dispatcher.Invoke(smallCube, bigCube);
}
template< typename ArrayType >
void Print2DArray( const std::string &str, const ArrayType &arr, vtkm::Id dimX )
{
@ -316,7 +366,9 @@ public:
{
std::cerr << arr.GetPortalConstControl().Get(i) << " ";
if( i % dimX == dimX - 1 )
{
std::cerr << std::endl;
}
}
}
@ -330,11 +382,15 @@ protected:
void WaveLengthValidate( vtkm::Id sigInLen, vtkm::Id filterLength, vtkm::Id &level)
{
if( sigInLen < filterLength )
{
level = 0;
}
else
level = static_cast<vtkm::Id>( vtkm::Floor(
{
level = static_cast<vtkm::Id>( vtkm::Floor( 1.0 +
vtkm::Log2( static_cast<vtkm::Float64>(sigInLen) /
static_cast<vtkm::Float64>(filterLength) ) + 1.0 ) );
static_cast<vtkm::Float64>(filterLength) ) ) );
}
}
}; // class WaveletBase.

1388
vtkm/worklet/wavelets/WaveletDWT.h
View File

File diff suppressed because it is too large Load Diff

21
vtkm/worklet/wavelets/WaveletFilter.h
View File

@ -40,7 +40,7 @@ enum WaveletName {
BIOR4_4, // the same as CDF9_7
BIOR3_3, // the same as CDF8_4
BIOR2_2, // the same as CDF5_3
BIOR1_1 // the same as HAAE
BIOR1_1 // the same as HAAR
};
// Wavelet filter class;
@ -106,8 +106,15 @@ public:
}
}
vtkm::Id GetFilterLength() { return this->FilterLength; }
bool isSymmetric() { return this->Symmetricity; }
vtkm::Id GetFilterLength()
{
return this->FilterLength;
}
bool isSymmetric()
{
return this->Symmetricity;
}
typedef vtkm::cont::ArrayHandle<vtkm::Float64> FilterType;
@ -160,7 +167,9 @@ private:
void wrev( const vtkm::Float64* arrIn, vtkm::Float64* arrOut, vtkm::Id length )
{
for( vtkm::Id count = 0; count < length; count++)
{
arrOut[count] = arrIn[length - count - 1];
}
}
// Quadrature mirror filtering operation: helper function to initialize a filter.
@ -172,7 +181,9 @@ private:
{
arrOut[count] = arrIn[length - count - 1];
if (count % 2 != 0)
{
arrOut[count] = -1.0 * arrOut[count];
}
}
}
else
@ -181,7 +192,9 @@ private:
{
arrOut[count] = arrIn[length - count - 1];
if (count % 2 == 0)
{
arrOut[count] = -1.0 * arrOut[count];
}
}
}
}
@ -203,7 +216,9 @@ private:
void verbatim_copy ( const vtkm::Float64* arrIn, vtkm::Float64* arrOut, vtkm::Id length )
{
for (vtkm::Id count = 0; count < length; count++)
{
arrOut[count] = arrIn[count];
}
}
}; // class WaveletFilter.

2330
vtkm/worklet/wavelets/WaveletTransforms.h
View File

File diff suppressed because it is too large Load Diff