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.
|
|
|
|
//============================================================================
|
|
|
|
|
|
|
|
#include <vtkm/worklet/WaveletCompressor.h>
|
|
|
|
|
|
|
|
#include <vtkm/cont/testing/Testing.h>
|
|
|
|
#include <vtkm/cont/ArrayHandlePermutation.h>
|
|
|
|
#include <vtkm/cont/Timer.h>
|
|
|
|
|
|
|
|
#include <vector>
|
2016-08-10 23:28:59 +00:00
|
|
|
#include <iomanip>
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-08-23 15:56:03 +00:00
|
|
|
namespace vtkm
|
|
|
|
{
|
|
|
|
namespace worklet
|
|
|
|
{
|
|
|
|
namespace wavelets
|
|
|
|
{
|
2016-09-22 03:25:56 +00:00
|
|
|
|
2016-08-23 15:56:03 +00:00
|
|
|
class SineWorklet : public vtkm::worklet::WorkletMapField
|
|
|
|
{
|
|
|
|
public:
|
|
|
|
typedef void ControlSignature(FieldInOut<>);
|
|
|
|
typedef void ExecutionSignature(_1, WorkIndex);
|
|
|
|
|
|
|
|
template<typename T>
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_EXEC
|
2016-08-23 15:56:03 +00:00
|
|
|
void operator()(T& x, const vtkm::Id& workIdx) const
|
|
|
|
{
|
|
|
|
x = vtkm::Sin(vtkm::Float64(workIdx) / 100.0) * 100.0;
|
|
|
|
}
|
|
|
|
};
|
2016-09-22 03:25:56 +00:00
|
|
|
|
|
|
|
class GaussianWorklet2D : public vtkm::worklet::WorkletMapField
|
|
|
|
{
|
|
|
|
public:
|
|
|
|
typedef void ControlSignature(FieldInOut<>);
|
|
|
|
typedef void ExecutionSignature(_1, WorkIndex);
|
|
|
|
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_EXEC
|
2016-09-22 03:25:56 +00:00
|
|
|
GaussianWorklet2D( vtkm::Id dx, vtkm::Id dy, vtkm::Float64 a,
|
|
|
|
vtkm::Float64 x, vtkm::Float64 y,
|
|
|
|
vtkm::Float64 sx, vtkm::Float64 xy )
|
|
|
|
: dimX( dx ), dimY( dy ), amp (a),
|
|
|
|
x0( x ), y0( y ),
|
|
|
|
sigmaX( sx ), sigmaY( xy )
|
|
|
|
{
|
|
|
|
sigmaX2 = 2 * sigmaX * sigmaX;
|
|
|
|
sigmaY2 = 2 * sigmaY * sigmaY;
|
|
|
|
}
|
|
|
|
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_EXEC
|
2016-09-22 03:25:56 +00:00
|
|
|
void Sig1Dto2D( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y ) const
|
|
|
|
{
|
|
|
|
x = idx % dimX;
|
|
|
|
y = idx / dimX;
|
|
|
|
}
|
|
|
|
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_EXEC
|
2016-09-22 03:25:56 +00:00
|
|
|
vtkm::Float64 GetGaussian( vtkm::Float64 x, vtkm::Float64 y ) const
|
|
|
|
{
|
|
|
|
vtkm::Float64 power = (x-x0) * (x-x0) / sigmaX2 + (y-y0) * (y-y0) / sigmaY2;
|
|
|
|
return vtkm::Exp( power * -1.0 ) * amp;
|
|
|
|
}
|
|
|
|
|
|
|
|
template<typename T>
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_EXEC
|
2016-09-22 03:25:56 +00:00
|
|
|
void operator()(T& val, const vtkm::Id& workIdx) const
|
|
|
|
{
|
|
|
|
vtkm::Id x, y;
|
|
|
|
Sig1Dto2D( workIdx, x, y );
|
|
|
|
val = GetGaussian( static_cast<vtkm::Float64>(x), static_cast<vtkm::Float64>(y) );
|
|
|
|
}
|
|
|
|
|
|
|
|
private: // see wikipedia page
|
|
|
|
const vtkm::Id dimX, dimY; // 2D extent
|
|
|
|
const vtkm::Float64 amp; // amplitude
|
|
|
|
const vtkm::Float64 x0, y0; // center
|
|
|
|
const vtkm::Float64 sigmaX, sigmaY; // spread
|
|
|
|
vtkm::Float64 sigmaX2, sigmaY2; // 2 * sigma * sigma
|
|
|
|
};
|
|
|
|
|
2017-02-24 23:40:14 +00:00
|
|
|
template< typename T>
|
2017-02-21 23:56:57 +00:00
|
|
|
class GaussianWorklet3D : public vtkm::worklet::WorkletMapField
|
|
|
|
{
|
|
|
|
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 )
|
2017-02-26 00:19:33 +00:00
|
|
|
{
|
|
|
|
amp = 10.0;
|
|
|
|
sigmaX = (T)dimX / 4.0; sigmaX2 = sigmaX * sigmaX * 2.0;
|
|
|
|
sigmaY = (T)dimY / 4.0; sigmaY2 = sigmaY * sigmaY * 2.0;
|
|
|
|
sigmaZ = (T)dimZ / 4.0; sigmaZ2 = sigmaZ * sigmaZ * 2.0;
|
|
|
|
}
|
2017-02-21 23:56:57 +00:00
|
|
|
|
|
|
|
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
|
2017-02-24 23:40:14 +00:00
|
|
|
T GetGaussian( T x, T y, T z ) const
|
2017-02-21 23:56:57 +00:00
|
|
|
{
|
2017-02-26 00:19:33 +00:00
|
|
|
x -= (T)dimX / 2.0; // translate to center at (0, 0, 0)
|
|
|
|
y -= (T)dimY / 2.0;
|
|
|
|
z -= (T)dimZ / 2.0;
|
|
|
|
T power = x*x/sigmaX2 + y*y/sigmaY2 + z*z/sigmaZ2;
|
|
|
|
|
|
|
|
return vtkm::Exp( power * -1.0 ) * amp;
|
2017-02-21 23:56:57 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
VTKM_EXEC
|
|
|
|
void operator()(T& val, const vtkm::Id& workIdx) const
|
|
|
|
{
|
|
|
|
vtkm::Id x, y, z;
|
|
|
|
Sig1Dto3D( workIdx, x, y, z );
|
2017-02-26 00:19:33 +00:00
|
|
|
val = GetGaussian( (T)x, (T)y, (T)z );
|
2017-02-21 23:56:57 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
private:
|
2017-02-26 00:19:33 +00:00
|
|
|
const vtkm::Id dimX, dimY, dimZ; // extent
|
|
|
|
T amp; // amplitude
|
|
|
|
T sigmaX, sigmaY, sigmaZ; // spread
|
|
|
|
T sigmaX2, sigmaY2, sigmaZ2; // sigma * sigma * 2
|
2017-02-21 23:56:57 +00:00
|
|
|
};
|
|
|
|
|
2016-08-23 15:56:03 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
template< typename ArrayType >
|
|
|
|
void FillArray( ArrayType& array )
|
|
|
|
{
|
|
|
|
typedef vtkm::worklet::wavelets::SineWorklet SineWorklet;
|
|
|
|
SineWorklet worklet;
|
|
|
|
vtkm::worklet::DispatcherMapField< SineWorklet > dispatcher( worklet );
|
|
|
|
dispatcher.Invoke( array );
|
|
|
|
}
|
|
|
|
|
2016-09-22 03:25:56 +00:00
|
|
|
template< typename ArrayType >
|
|
|
|
void FillArray2D( ArrayType& array, vtkm::Id dimX, vtkm::Id dimY )
|
2016-08-26 16:54:47 +00:00
|
|
|
{
|
2016-09-22 03:25:56 +00:00
|
|
|
typedef vtkm::worklet::wavelets::GaussianWorklet2D WorkletType;
|
2016-09-22 23:13:50 +00:00
|
|
|
WorkletType worklet( dimX, dimY, 100.0,
|
2016-09-22 03:25:56 +00:00
|
|
|
static_cast<vtkm::Float64>(dimX)/2.0, // center
|
|
|
|
static_cast<vtkm::Float64>(dimY)/2.0, // center
|
|
|
|
static_cast<vtkm::Float64>(dimX)/4.0, // spread
|
|
|
|
static_cast<vtkm::Float64>(dimY)/4.0);// spread
|
|
|
|
vtkm::worklet::DispatcherMapField< WorkletType > dispatcher( worklet );
|
|
|
|
dispatcher.Invoke( array );
|
2016-09-11 22:14:01 +00:00
|
|
|
}
|
|
|
|
|
2017-02-21 23:56:57 +00:00
|
|
|
template< typename ArrayType >
|
|
|
|
void FillArray3D( ArrayType& array, vtkm::Id dimX, vtkm::Id dimY, vtkm::Id dimZ )
|
|
|
|
{
|
2017-02-24 23:40:14 +00:00
|
|
|
typedef vtkm::worklet::wavelets::GaussianWorklet3D< typename ArrayType::ValueType> WorkletType;
|
2017-02-21 23:56:57 +00:00
|
|
|
WorkletType worklet( dimX, dimY, dimZ );
|
|
|
|
vtkm::worklet::DispatcherMapField< WorkletType > dispatcher( worklet );
|
|
|
|
dispatcher.Invoke( array );
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
void TestDecomposeReconstruct3D()
|
|
|
|
{
|
|
|
|
std::cout << "Testing 3D wavelet compressor on a 10x10x10 cube: " << std::endl;
|
2017-02-25 20:23:40 +00:00
|
|
|
vtkm::Id sigX = 20;
|
|
|
|
vtkm::Id sigY = 20;
|
|
|
|
vtkm::Id sigZ = 20;
|
2017-02-21 23:56:57 +00:00
|
|
|
vtkm::Id sigLen = sigX * sigY * sigZ;
|
|
|
|
|
|
|
|
// make input data array handle
|
2017-02-24 23:40:14 +00:00
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float32> inputArray;
|
2017-02-21 23:56:57 +00:00
|
|
|
inputArray.PrepareForOutput( sigLen, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
|
|
|
FillArray3D( inputArray, sigX, sigY, sigZ );
|
|
|
|
|
2017-02-24 23:40:14 +00:00
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float32> outputArray;
|
2017-02-21 23:56:57 +00:00
|
|
|
|
|
|
|
// Use a WaveletCompressor
|
|
|
|
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF9_7;
|
|
|
|
std::cout << "Wavelet kernel = 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 );
|
2017-02-26 00:19:33 +00:00
|
|
|
nLevels = 2;
|
2017-02-21 23:56:57 +00:00
|
|
|
std::cout << "Decomposition levels = " << nLevels << std::endl;
|
|
|
|
vtkm::Float64 computationTime = 0.0;
|
|
|
|
vtkm::Float64 elapsedTime1, elapsedTime2, elapsedTime3;
|
|
|
|
|
|
|
|
// Decompose
|
|
|
|
//vtkm::cont::Timer<> timer;
|
|
|
|
//computationTime =
|
2017-02-25 20:23:40 +00:00
|
|
|
compressor.WaveDecompose3D( inputArray, nLevels, sigX, sigY, sigZ, outputArray,
|
2017-02-21 23:56:57 +00:00
|
|
|
VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
|
|
|
//elapsedTime1 = timer.GetElapsedTime();
|
|
|
|
//std::cout << "Decompose time = " << elapsedTime1 << std::endl;
|
|
|
|
//std::cout << " ->computation time = " << computationTime << std::endl;
|
|
|
|
|
|
|
|
std::cout << "finish decomposition" << std::endl;
|
|
|
|
|
|
|
|
// Squash small coefficients
|
|
|
|
//timer.Reset();
|
2017-02-26 00:19:33 +00:00
|
|
|
vtkm::Float64 cratio = 4.0; // X:1 compression, where X >= 1
|
|
|
|
compressor.SquashCoefficients( outputArray, cratio, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2017-02-21 23:56:57 +00:00
|
|
|
//elapsedTime2 = timer.GetElapsedTime();
|
|
|
|
//std::cout << "Squash time = " << elapsedTime2 << std::endl;
|
|
|
|
|
|
|
|
// Reconstruct
|
2017-02-24 23:40:14 +00:00
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float32> reconstructArray;
|
2017-02-21 23:56:57 +00:00
|
|
|
//timer.Reset();
|
|
|
|
//computationTime =
|
2017-02-25 20:23:40 +00:00
|
|
|
compressor.WaveReconstruct3D( outputArray, nLevels, sigX, sigY, sigZ, reconstructArray,
|
2017-02-21 23:56:57 +00:00
|
|
|
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;
|
|
|
|
|
|
|
|
std::cout << "finish reconstruction" << 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;
|
|
|
|
}
|
|
|
|
|
2016-09-11 22:14:01 +00:00
|
|
|
|
2016-08-06 07:08:29 +00:00
|
|
|
|
2016-08-22 16:13:45 +00:00
|
|
|
void TestDecomposeReconstruct2D()
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2016-10-10 04:52:02 +00:00
|
|
|
std::cout << "Testing 2D wavelet compressor on a 1000x1000 square: " << std::endl;
|
2016-10-10 00:07:49 +00:00
|
|
|
vtkm::Id sigX = 1000;
|
|
|
|
vtkm::Id sigY = 1000;
|
2016-08-22 16:13:45 +00:00
|
|
|
vtkm::Id sigLen = sigX * sigY;
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
// make input data array handle
|
2016-08-23 15:56:03 +00:00
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> inputArray;
|
|
|
|
inputArray.PrepareForOutput( sigLen, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-09-22 03:25:56 +00:00
|
|
|
FillArray2D( inputArray, sigX, sigY );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> outputArray;
|
|
|
|
|
|
|
|
// Use a WaveletCompressor
|
2016-10-10 04:52:02 +00:00
|
|
|
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF9_7;
|
|
|
|
std::cout << "Wavelet kernel = CDF 9/7" << std::endl;
|
2016-08-10 23:24:22 +00:00
|
|
|
vtkm::worklet::WaveletCompressor compressor( wname );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-08-22 16:13:45 +00:00
|
|
|
vtkm::Id XMaxLevel = compressor.GetWaveletMaxLevel( sigX );
|
|
|
|
vtkm::Id YMaxLevel = compressor.GetWaveletMaxLevel( sigY );
|
|
|
|
vtkm::Id nLevels = vtkm::Min( XMaxLevel, YMaxLevel );
|
2016-09-06 21:28:46 +00:00
|
|
|
std::cout << "Decomposition levels = " << nLevels << std::endl;
|
2016-08-08 23:59:15 +00:00
|
|
|
std::vector<vtkm::Id> L;
|
2016-09-06 22:54:12 +00:00
|
|
|
vtkm::Float64 computationTime = 0.0;
|
2016-09-22 03:25:56 +00:00
|
|
|
vtkm::Float64 elapsedTime1, elapsedTime2, elapsedTime3;
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-07-29 18:27:37 +00:00
|
|
|
// Decompose
|
2016-09-11 22:14:01 +00:00
|
|
|
vtkm::cont::Timer<> timer;
|
2016-09-06 22:54:12 +00:00
|
|
|
computationTime =
|
2016-08-22 16:13:45 +00:00
|
|
|
compressor.WaveDecompose2D( inputArray, nLevels, sigX, sigY, outputArray, L,
|
|
|
|
VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-09-22 03:25:56 +00:00
|
|
|
elapsedTime1 = timer.GetElapsedTime();
|
|
|
|
std::cout << "Decompose time = " << elapsedTime1 << std::endl;
|
2016-08-23 23:52:40 +00:00
|
|
|
std::cout << " ->computation time = " << computationTime << std::endl;
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-08-23 15:56:03 +00:00
|
|
|
// Squash small coefficients
|
2016-10-10 05:34:32 +00:00
|
|
|
timer.Reset();
|
2016-10-10 04:52:02 +00:00
|
|
|
vtkm::Float64 cratio = 1.0; // X:1 compression, where X >= 1
|
2016-08-23 15:56:03 +00:00
|
|
|
compressor.SquashCoefficients( outputArray, cratio, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-09-22 03:25:56 +00:00
|
|
|
elapsedTime2 = timer.GetElapsedTime();
|
2016-10-10 05:34:32 +00:00
|
|
|
std::cout << "Squash time = " << elapsedTime2 << std::endl;
|
2016-08-23 15:56:03 +00:00
|
|
|
|
2016-07-27 20:28:38 +00:00
|
|
|
// Reconstruct
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> reconstructArray;
|
2016-09-11 22:14:01 +00:00
|
|
|
timer.Reset();
|
2016-08-23 23:52:40 +00:00
|
|
|
computationTime =
|
2016-08-22 16:13:45 +00:00
|
|
|
compressor.WaveReconstruct2D( outputArray, nLevels, sigX, sigY, reconstructArray, L,
|
|
|
|
VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-09-22 03:25:56 +00:00
|
|
|
elapsedTime3 = timer.GetElapsedTime();
|
|
|
|
std::cout << "Reconstruction time = " << elapsedTime3 << std::endl;
|
2016-08-23 23:52:40 +00:00
|
|
|
std::cout << " ->computation time = " << computationTime << std::endl;
|
2016-09-22 03:25:56 +00:00
|
|
|
std::cout << "Total time = "
|
|
|
|
<< (elapsedTime1 + elapsedTime2 + elapsedTime3) << std::endl;
|
|
|
|
|
|
|
|
outputArray.ReleaseResources();
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-08-10 20:33:14 +00:00
|
|
|
compressor.EvaluateReconstruction( inputArray, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-07-29 03:07:05 +00:00
|
|
|
|
2016-09-11 22:14:01 +00:00
|
|
|
timer.Reset();
|
2016-07-27 20:28:38 +00:00
|
|
|
for( vtkm::Id i = 0; i < reconstructArray.GetNumberOfValues(); i++ )
|
|
|
|
{
|
2016-09-06 21:28:46 +00:00
|
|
|
VTKM_TEST_ASSERT( test_equal( reconstructArray.GetPortalConstControl().Get(i),
|
|
|
|
inputArray.GetPortalConstControl().Get(i) ),
|
2016-10-11 14:26:14 +00:00
|
|
|
"WaveletCompressor 2D failed..." );
|
2016-07-27 20:28:38 +00:00
|
|
|
}
|
2016-09-22 03:25:56 +00:00
|
|
|
elapsedTime1 = timer.GetElapsedTime();
|
|
|
|
std::cout << "Verification time = " << elapsedTime1 << std::endl;
|
2016-07-29 21:16:10 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
|
2016-08-23 15:56:03 +00:00
|
|
|
void TestDecomposeReconstruct1D()
|
2016-07-29 21:16:10 +00:00
|
|
|
{
|
2016-10-10 04:52:02 +00:00
|
|
|
std::cout << "Testing 1D wavelet compressor on a 1 million sized array " << std::endl;
|
2016-07-29 21:16:10 +00:00
|
|
|
vtkm::Id million = 1000000;
|
2016-10-10 04:52:02 +00:00
|
|
|
vtkm::Id sigLen = million * 1;
|
2016-07-29 21:16:10 +00:00
|
|
|
|
|
|
|
// make input data array handle
|
|
|
|
std::vector<vtkm::Float64> tmpVector;
|
|
|
|
for( vtkm::Id i = 0; i < sigLen; i++ )
|
2016-08-16 21:20:19 +00:00
|
|
|
{
|
2016-07-29 21:16:10 +00:00
|
|
|
tmpVector.push_back( 100.0 * vtkm::Sin(static_cast<vtkm::Float64>(i)/100.0 ));
|
2016-08-16 21:20:19 +00:00
|
|
|
}
|
2016-07-29 21:16:10 +00:00
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> inputArray =
|
|
|
|
vtkm::cont::make_ArrayHandle(tmpVector);
|
|
|
|
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> outputArray;
|
|
|
|
|
|
|
|
// Use a WaveletCompressor
|
2016-10-10 04:52:02 +00:00
|
|
|
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF9_7;
|
|
|
|
std::cout << "Wavelet kernel = CDF 9/7" << std::endl;
|
2016-08-10 23:24:22 +00:00
|
|
|
vtkm::worklet::WaveletCompressor compressor( wname );
|
2016-07-29 21:16:10 +00:00
|
|
|
|
|
|
|
// User maximum decompose levels, and no compression
|
|
|
|
vtkm::Id maxLevel = compressor.GetWaveletMaxLevel( sigLen );
|
|
|
|
vtkm::Id nLevels = maxLevel;
|
2016-10-10 04:52:02 +00:00
|
|
|
std::cout << "Decomposition levels = " << nLevels << std::endl;
|
2016-07-29 21:16:10 +00:00
|
|
|
|
2016-08-08 23:59:15 +00:00
|
|
|
std::vector<vtkm::Id> L;
|
2016-07-29 21:16:10 +00:00
|
|
|
|
|
|
|
// Decompose
|
|
|
|
vtkm::cont::Timer<> timer;
|
2016-08-10 20:33:14 +00:00
|
|
|
compressor.WaveDecompose( inputArray, nLevels, outputArray, L, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-07-29 21:16:10 +00:00
|
|
|
|
|
|
|
vtkm::Float64 elapsedTime = timer.GetElapsedTime();
|
|
|
|
std::cout << "Decompose time = " << elapsedTime << std::endl;
|
|
|
|
|
|
|
|
// Reconstruct
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> reconstructArray;
|
|
|
|
timer.Reset();
|
2016-08-10 20:33:14 +00:00
|
|
|
compressor.WaveReconstruct( outputArray, nLevels, L, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-07-29 21:16:10 +00:00
|
|
|
elapsedTime = timer.GetElapsedTime();
|
|
|
|
std::cout << "Reconstruction time = " << elapsedTime << std::endl;
|
|
|
|
|
2016-08-10 20:33:14 +00:00
|
|
|
compressor.EvaluateReconstruction( inputArray, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-07-29 21:16:10 +00:00
|
|
|
|
|
|
|
timer.Reset();
|
|
|
|
for( vtkm::Id i = 0; i < reconstructArray.GetNumberOfValues(); i++ )
|
|
|
|
{
|
2016-09-22 23:13:50 +00:00
|
|
|
VTKM_TEST_ASSERT( test_equal( reconstructArray.GetPortalConstControl().Get(i),
|
|
|
|
inputArray.GetPortalConstControl().Get(i)),
|
2016-10-11 14:26:14 +00:00
|
|
|
"WaveletCompressor 1D failed..." );
|
2016-07-29 21:16:10 +00:00
|
|
|
}
|
|
|
|
elapsedTime = timer.GetElapsedTime();
|
|
|
|
std::cout << "Verification time = " << elapsedTime << std::endl;
|
2016-07-27 20:28:38 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
void TestWaveletCompressor()
|
|
|
|
{
|
2017-02-21 23:56:57 +00:00
|
|
|
//TestDecomposeReconstruct1D();
|
|
|
|
//std::cout << std::endl;
|
|
|
|
//TestDecomposeReconstruct2D();
|
|
|
|
TestDecomposeReconstruct3D();
|
2016-07-27 20:28:38 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
int UnitTestWaveletCompressor(int, char *[])
|
|
|
|
{
|
|
|
|
return vtkm::cont::testing::Testing::Run(TestWaveletCompressor);
|
|
|
|
}
|