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
|
|
|
|
{
|
|
|
|
class SineWorklet : public vtkm::worklet::WorkletMapField
|
|
|
|
{
|
|
|
|
public:
|
|
|
|
typedef void ControlSignature(FieldInOut<>);
|
|
|
|
typedef void ExecutionSignature(_1, WorkIndex);
|
|
|
|
|
|
|
|
template<typename T>
|
|
|
|
VTKM_EXEC_EXPORT
|
|
|
|
void operator()(T& x, const vtkm::Id& workIdx) const
|
|
|
|
{
|
|
|
|
x = vtkm::Sin(vtkm::Float64(workIdx) / 100.0) * 100.0;
|
|
|
|
}
|
|
|
|
};
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
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-08-26 16:54:47 +00:00
|
|
|
void Debug2DExtend()
|
|
|
|
{
|
2016-08-31 22:03:25 +00:00
|
|
|
vtkm::Id NX = 10;
|
|
|
|
vtkm::Id NY = 11;
|
2016-09-02 00:47:50 +00:00
|
|
|
typedef vtkm::cont::ArrayHandle< vtkm::Float64 > ArrayType;
|
|
|
|
ArrayType left, center, right;
|
|
|
|
|
2016-08-31 22:03:25 +00:00
|
|
|
center.PrepareForOutput( NX * NY, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
|
|
|
for( vtkm::Id i = 0; i < NX*NY; i++ )
|
2016-08-26 20:51:23 +00:00
|
|
|
center.GetPortalControl().Set(i, i);
|
2016-08-26 16:54:47 +00:00
|
|
|
|
2016-09-02 00:47:50 +00:00
|
|
|
ArrayType output1, output2;
|
2016-08-31 22:03:25 +00:00
|
|
|
std::vector<vtkm::Id> L(10, 0);
|
2016-08-26 20:51:23 +00:00
|
|
|
|
|
|
|
vtkm::worklet::wavelets::WaveletDWT dwt( vtkm::worklet::wavelets::CDF9_7 );
|
2016-09-02 00:47:50 +00:00
|
|
|
|
2016-08-31 22:03:25 +00:00
|
|
|
dwt.DWT2D(center, NX, NY, output1, L, VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
|
2016-09-02 00:47:50 +00:00
|
|
|
|
|
|
|
dwt.DWT2Dv2( center, NX, NY, output2, L, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
|
|
|
|
|
|
|
for( vtkm::Id i = 0; i < output1.GetNumberOfValues(); i++ )
|
|
|
|
{
|
|
|
|
VTKM_TEST_ASSERT( test_equal( output1.GetPortalConstControl().Get(i),
|
|
|
|
output2.GetPortalConstControl().Get(i)),
|
|
|
|
"WaveletCompressor worklet failed..." );
|
|
|
|
}
|
|
|
|
|
|
|
|
std::cout << "true results:" << std::endl;
|
|
|
|
for( vtkm::Id i = 0; i < output1.GetNumberOfValues(); i++ )
|
|
|
|
{
|
|
|
|
std::cout << output1.GetPortalConstControl().Get(i) << " ";
|
|
|
|
if( i % NX == NX - 1 )
|
|
|
|
std::cout << std::endl;
|
|
|
|
}
|
2016-08-31 22:03:25 +00:00
|
|
|
|
|
|
|
std::cout << "test results:" << std::endl;
|
2016-09-02 00:47:50 +00:00
|
|
|
for( vtkm::Id i = 0; i < output2.GetNumberOfValues(); i++ )
|
2016-08-26 16:54:47 +00:00
|
|
|
{
|
2016-09-02 00:47:50 +00:00
|
|
|
std::cout << output2.GetPortalConstControl().Get(i) << " ";
|
|
|
|
if( i % NX == NX - 1 )
|
|
|
|
std::cout << std::endl;
|
2016-08-26 16:54:47 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2016-07-29 21:16:10 +00:00
|
|
|
void DebugDWTIDWT1D()
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2016-08-16 21:20:19 +00:00
|
|
|
vtkm::Id sigLen = 21;
|
2016-07-27 20:28:38 +00:00
|
|
|
std::cout << "Testing Wavelets Worklet" << std::endl;
|
2016-07-28 23:55:38 +00:00
|
|
|
std::cout << "Input a size to test." << std::endl;
|
2016-08-16 21:20:19 +00:00
|
|
|
std::cin >> sigLen;
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
// make input data array handle
|
|
|
|
std::vector<vtkm::Float64> tmpVector;
|
|
|
|
for( vtkm::Id i = 0; i < sigLen; i++ )
|
|
|
|
tmpVector.push_back( static_cast<vtkm::Float64>( i ) );
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> inputArray =
|
|
|
|
vtkm::cont::make_ArrayHandle(tmpVector);
|
|
|
|
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> coeffOut;
|
2016-08-08 23:59:15 +00:00
|
|
|
std::vector<vtkm::Id> L(3, 0);
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
// Forward Transform
|
2016-08-16 21:20:19 +00:00
|
|
|
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF8_4;
|
2016-08-10 23:24:22 +00:00
|
|
|
vtkm::worklet::wavelets::WaveletDWT waveletdwt( wname );
|
2016-08-22 20:49:13 +00:00
|
|
|
waveletdwt.DWT1D( inputArray, coeffOut, L, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
std::cout << "Forward Wavelet Transform: result coeff length = " <<
|
|
|
|
coeffOut.GetNumberOfValues() << std::endl;
|
2016-07-27 22:42:29 +00:00
|
|
|
|
2016-07-28 23:55:38 +00:00
|
|
|
for( vtkm::Id i = 0; i < coeffOut.GetNumberOfValues(); i++ )
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
|
|
|
if( i == 0 )
|
|
|
|
std::cout << " <-- cA --> " << std::endl;
|
|
|
|
else if( i == L[0] )
|
|
|
|
std::cout << " <-- cD --> " << std::endl;
|
|
|
|
std::cout << coeffOut.GetPortalConstControl().Get(i) << std::endl;
|
|
|
|
}
|
|
|
|
|
|
|
|
// Inverse Transform
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> reconstructArray;
|
2016-08-22 20:49:13 +00:00
|
|
|
waveletdwt.IDWT1D( coeffOut, L, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-07-27 20:28:38 +00:00
|
|
|
std::cout << "Inverse Wavelet Transform: result signal length = " <<
|
|
|
|
reconstructArray.GetNumberOfValues() << std::endl;
|
2016-07-28 23:55:38 +00:00
|
|
|
for( vtkm::Id i = 0; i < reconstructArray.GetNumberOfValues(); i++ )
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
|
|
|
std::cout << reconstructArray.GetPortalConstControl().Get(i) << std::endl;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2016-08-06 07:08:29 +00:00
|
|
|
|
2016-08-14 06:28:29 +00:00
|
|
|
void DebugRectangleCopy()
|
|
|
|
{
|
2016-08-14 22:29:27 +00:00
|
|
|
vtkm::Id sigX = 12;
|
|
|
|
vtkm::Id sigY = 12;
|
2016-08-14 06:28:29 +00:00
|
|
|
vtkm::Id sigLen = sigX * sigY;
|
|
|
|
|
|
|
|
// make input data array handle
|
|
|
|
std::vector<vtkm::Float64> tmpVector;
|
|
|
|
for( vtkm::Id i = 0; i < sigLen; i++ )
|
|
|
|
tmpVector.push_back( static_cast<vtkm::Float64>( i ) );
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> inputArray =
|
|
|
|
vtkm::cont::make_ArrayHandle(tmpVector);
|
|
|
|
|
2016-08-14 22:29:27 +00:00
|
|
|
// make output array
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> outputArray;
|
2016-08-14 06:28:29 +00:00
|
|
|
|
2016-08-14 22:29:27 +00:00
|
|
|
// make bookkeeping array
|
|
|
|
std::vector<vtkm::Id> L;
|
2016-08-14 06:28:29 +00:00
|
|
|
|
2016-08-22 05:10:05 +00:00
|
|
|
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF5_3;
|
|
|
|
vtkm::worklet::WaveletCompressor wavelet( wname );
|
|
|
|
wavelet.WaveDecompose2D( inputArray, 2, sigX, sigY, outputArray, L,
|
|
|
|
VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-08-14 06:28:29 +00:00
|
|
|
|
2016-08-14 22:29:27 +00:00
|
|
|
for( vtkm::Id i = 0; i < outputArray.GetNumberOfValues(); i++ )
|
2016-08-14 06:28:29 +00:00
|
|
|
{
|
2016-08-14 22:29:27 +00:00
|
|
|
std::cout << std::setw( 10 );
|
|
|
|
std::cout << outputArray.GetPortalConstControl().Get(i) << "\t";
|
2016-08-14 06:28:29 +00:00
|
|
|
if( i % sigX == sigX-1 )
|
|
|
|
std::cout << std::endl;
|
|
|
|
}
|
|
|
|
std::cout << std::endl;
|
|
|
|
}
|
|
|
|
|
2016-08-22 16:13:45 +00:00
|
|
|
|
|
|
|
void TestDecomposeReconstruct2D()
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2016-08-22 16:13:45 +00:00
|
|
|
vtkm::Id sigX = 1000;
|
|
|
|
vtkm::Id sigY = 2000;
|
2016-08-23 15:56:03 +00:00
|
|
|
std::cout << "Please input X to test a X^2 square: " << std::endl;
|
|
|
|
std::cin >> sigX;
|
|
|
|
sigY = sigX;
|
|
|
|
//std::cout << "Testing wavelet compressor on 1000x2000 rectangle" << std::endl;
|
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() );
|
|
|
|
FillArray( inputArray );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> outputArray;
|
|
|
|
|
|
|
|
// Use a WaveletCompressor
|
2016-08-23 15:56:03 +00:00
|
|
|
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF9_7;
|
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-08-23 23:52:40 +00:00
|
|
|
//nLevels = 1;
|
2016-08-08 23:59:15 +00:00
|
|
|
std::vector<vtkm::Id> L;
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-07-29 18:27:37 +00:00
|
|
|
// Decompose
|
2016-07-27 20:28:38 +00:00
|
|
|
vtkm::cont::Timer<> timer;
|
2016-08-23 23:52:40 +00:00
|
|
|
vtkm::Float64 computationTime =
|
2016-08-22 16:13:45 +00:00
|
|
|
compressor.WaveDecompose2D( inputArray, nLevels, sigX, sigY, outputArray, L,
|
|
|
|
VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-07-27 20:28:38 +00:00
|
|
|
vtkm::Float64 elapsedTime = timer.GetElapsedTime();
|
2016-07-29 18:27:37 +00:00
|
|
|
std::cout << "Decompose time = " << elapsedTime << 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
|
|
|
|
timer.Reset();
|
2016-08-23 23:52:40 +00:00
|
|
|
vtkm::Float64 cratio = 1.0;
|
2016-08-23 15:56:03 +00:00
|
|
|
compressor.SquashCoefficients( outputArray, cratio, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
|
|
|
elapsedTime = timer.GetElapsedTime();
|
|
|
|
std::cout << "Squash time = " << elapsedTime << std::endl;
|
|
|
|
|
2016-07-27 20:28:38 +00:00
|
|
|
// Reconstruct
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> reconstructArray;
|
|
|
|
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-07-27 20:28:38 +00:00
|
|
|
elapsedTime = timer.GetElapsedTime();
|
2016-07-29 18:27:37 +00:00
|
|
|
std::cout << "Reconstruction time = " << elapsedTime << 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-10 20:33:14 +00:00
|
|
|
compressor.EvaluateReconstruction( inputArray, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
|
2016-07-29 03:07:05 +00:00
|
|
|
|
2016-07-27 20:28:38 +00:00
|
|
|
timer.Reset();
|
|
|
|
for( vtkm::Id i = 0; i < reconstructArray.GetNumberOfValues(); i++ )
|
|
|
|
{
|
|
|
|
VTKM_TEST_ASSERT( test_equal( reconstructArray.GetPortalConstControl().Get(i),
|
2016-08-23 15:56:03 +00:00
|
|
|
inputArray.GetPortalConstControl().Get(i) ),
|
2016-07-27 20:28:38 +00:00
|
|
|
"output value not the same..." );
|
|
|
|
}
|
|
|
|
elapsedTime = timer.GetElapsedTime();
|
2016-07-29 21:16:10 +00:00
|
|
|
std::cout << "Verification time = " << elapsedTime << std::endl;
|
|
|
|
}
|
|
|
|
|
|
|
|
|
2016-08-23 15:56:03 +00:00
|
|
|
void TestDecomposeReconstruct1D()
|
2016-07-29 21:16:10 +00:00
|
|
|
{
|
|
|
|
std::cout << "Testing WaveletCompressor on a 2 million sized array " << std::endl;
|
|
|
|
vtkm::Id million = 1000000;
|
|
|
|
vtkm::Id sigLen = million * 2;
|
|
|
|
|
|
|
|
// 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-08-16 21:20:19 +00:00
|
|
|
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF8_4;
|
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-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++ )
|
|
|
|
{
|
|
|
|
VTKM_TEST_ASSERT( test_equal( reconstructArray.GetPortalConstControl().Get(i),
|
|
|
|
100.0 * vtkm::Sin( static_cast<vtkm::Float64>(i)/100.0 )),
|
|
|
|
"WaveletCompressor worklet failed..." );
|
|
|
|
}
|
|
|
|
elapsedTime = timer.GetElapsedTime();
|
|
|
|
std::cout << "Verification time = " << elapsedTime << std::endl;
|
2016-07-27 20:28:38 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
void TestWaveletCompressor()
|
|
|
|
{
|
2016-07-29 21:16:10 +00:00
|
|
|
//DebugDWTIDWT1D();
|
2016-08-14 22:29:27 +00:00
|
|
|
//DebugRectangleCopy();
|
2016-08-23 15:56:03 +00:00
|
|
|
//TestDecomposeReconstruct1D();
|
2016-08-26 16:54:47 +00:00
|
|
|
//TestDecomposeReconstruct2D();
|
|
|
|
Debug2DExtend();
|
2016-07-27 20:28:38 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
int UnitTestWaveletCompressor(int, char *[])
|
|
|
|
{
|
|
|
|
return vtkm::cont::testing::Testing::Run(TestWaveletCompressor);
|
|
|
|
}
|