diff --git a/vtkm/filter/CMakeLists.txt b/vtkm/filter/CMakeLists.txt index c19739ff3..33a9b63b5 100644 --- a/vtkm/filter/CMakeLists.txt +++ b/vtkm/filter/CMakeLists.txt @@ -36,7 +36,6 @@ set(headers ResultField.h Threshold.h VertexClustering.h - WaveletCompressor.h ) set(header_template_sources @@ -50,7 +49,6 @@ set(header_template_sources PointElevation.hxx Threshold.hxx VertexClustering.hxx - WaveletCompressor.hxx ) vtkm_declare_headers(${headers}) diff --git a/vtkm/filter/internal/CMakeLists.txt b/vtkm/filter/internal/CMakeLists.txt index 621a8c1a4..00465f863 100644 --- a/vtkm/filter/internal/CMakeLists.txt +++ b/vtkm/filter/internal/CMakeLists.txt @@ -22,10 +22,6 @@ set(headers ResolveFieldTypeAndExecute.h ResolveFieldTypeAndMap.h RuntimeDeviceTracker.h - FilterBanks.h - WaveletFilter.h - WaveletBase.h - WaveletDWT.h ) vtkm_declare_headers(${headers}) diff --git a/vtkm/filter/testing/CMakeLists.txt b/vtkm/filter/testing/CMakeLists.txt index 1d31cd864..f8674bc6e 100644 --- a/vtkm/filter/testing/CMakeLists.txt +++ b/vtkm/filter/testing/CMakeLists.txt @@ -26,7 +26,6 @@ set(unit_tests UnitTestMarchingCubesFilter.cxx UnitTestThresholdFilter.cxx UnitTestVertexClusteringFilter.cxx - UnitTestWaveletCompressorFilter.cxx ) vtkm_unit_tests(SOURCES ${unit_tests}) diff --git a/vtkm/worklet/CMakeLists.txt b/vtkm/worklet/CMakeLists.txt index de0d70fb4..d731fb5d4 100644 --- a/vtkm/worklet/CMakeLists.txt +++ b/vtkm/worklet/CMakeLists.txt @@ -44,12 +44,13 @@ set(headers TriangulateUniformGrid.h Threshold.h VertexClustering.h - WaveletTransforms.h + WaveletCompressor.h ) #----------------------------------------------------------------------------- add_subdirectory(internal) add_subdirectory(splatkernels) +add_subdirectory(wavelets) vtkm_declare_headers(${headers}) diff --git a/vtkm/worklet/WaveletCompressor.h b/vtkm/worklet/WaveletCompressor.h new file mode 100644 index 000000000..11272706d --- /dev/null +++ b/vtkm/worklet/WaveletCompressor.h @@ -0,0 +1,238 @@ +//============================================================================ +// 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. +//============================================================================ + +#ifndef vtk_m_worklet_waveletcompressor_h +#define vtk_m_worklet_waveletcompressor_h + +#include + +//#include + +#include +#include +#include + +#include + +namespace vtkm { +namespace worklet { + +//template< typename DeviceAdapter > +class WaveletCompressor : public vtkm::worklet::wavelets::WaveletDWT +{ +public: + + // Constructor + WaveletCompressor( const std::string &w_name ) : WaveletDWT( w_name ) {} + + // Multi-level 1D wavelet decomposition + template< typename SignalArrayType, typename CoeffArrayType> + VTKM_EXEC_CONT_EXPORT + vtkm::Id WaveDecompose( const SignalArrayType &sigIn, // Input + vtkm::Id nLevels, // n levels of DWT + CoeffArrayType &coeffOut, + vtkm::Id* L ) + { + vtkm::Id sigInLen = sigIn.GetNumberOfValues(); + if( nLevels < 1 || nLevels > WaveletBase::GetWaveletMaxLevel( sigInLen ) ) + { + std::cerr << "nLevel is not supported: " << nLevels << std::endl; + // throw an error + } + /* 0 levels means no transform + if( nLevels == 0 ) + { + vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Copy + (sigIn, C ); + return 0; + } + */ + + this->ComputeL( sigInLen, nLevels, L ); + vtkm::Id CLength = this->ComputeCoeffLength( L, nLevels ); + + // Use 64bit floats for intermediate calculation + #define VAL vtkm::Float64 + + vtkm::Id sigInPtr = 0; // pseudo pointer for the beginning of input array + vtkm::Id len = sigIn.GetNumberOfValues(); + vtkm::Id cALen = WaveletBase::GetApproxLength( len ); + vtkm::Id cptr; // pseudo pointer for the beginning of output array + vtkm::Id tlen = 0; + vtkm::Id L1d[3]; + + // Use an intermediate array + typedef vtkm::cont::ArrayHandle< VAL > InterArrayType; + typedef typename InterArrayType::PortalControl InterPortalType; + InterArrayType interArray; + interArray.Allocate( CLength ); + vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Copy + (sigIn, interArray ); + + // Define a few more types + typedef vtkm::cont::ArrayHandleCounting< vtkm::Id > IdArrayType; + typedef vtkm::cont::ArrayHandlePermutation< IdArrayType, InterArrayType > PermutArrayType; + + for( vtkm::Id i = nLevels; i > 0; i-- ) + { + tlen += L[i]; + cptr = 0 + CLength - tlen - cALen; + + // make input array (permutation array) + IdArrayType inputIndices( sigInPtr, 1, len ); + PermutArrayType input( inputIndices, interArray ); + // make output array + InterArrayType output; + + WaveletDWT::DWT1D( input, output, L1d ); + + // update interArray + vtkm::cont::ArrayPortalToIterators< InterPortalType > + outputIter( output.GetPortalControl() ); + vtkm::cont::ArrayPortalToIterators< InterPortalType > + interArrayIter( interArray.GetPortalControl() ); + std::copy( outputIter.GetBegin(), outputIter.GetEnd(), interArrayIter.GetBegin() + cptr ); + + // update pseudo pointers + len = cALen; + cALen = WaveletBase::GetApproxLength( cALen ); + sigInPtr = cptr; + } + + vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Copy + (interArray, coeffOut ); + + #undef VAL + + return 0; + } + + // Multi-level 1D wavelet reconstruction + template< typename CoeffArrayType, typename SignalArrayType > + VTKM_EXEC_CONT_EXPORT + vtkm::Id WaveReconstruct( const CoeffArrayType &coeffIn, // Input + vtkm::Id nLevels, // n levels of DWT + vtkm::Id* L, + SignalArrayType &sigOut ) + { + VTKM_ASSERT( nLevels > 0 ); + + vtkm::Id LLength = nLevels + 2; + + vtkm::Id L1d[3] = {L[0], L[1], 0}; + + // Use 64bit floats for intermediate calculation + #define VAL vtkm::Float64 + + // Use intermediate arrays + typedef vtkm::cont::ArrayHandle< VAL > InterArrayType; + typedef typename InterArrayType::PortalControl InterPortalType; + typedef vtkm::cont::ArrayHandleCounting< vtkm::Id > IdArrayType; + typedef vtkm::cont::ArrayHandlePermutation< IdArrayType, InterArrayType > PermutArrayType; + + InterArrayType interArray; + vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Copy + (coeffIn, interArray ); + + for( vtkm::Id i = 1; i <= nLevels; i++ ) + { + L1d[2] = this->GetApproxLengthLevN( L[ LLength-1 ], nLevels-i ); + + // Make an input array + IdArrayType inputIndices( 0, 1, L1d[2] ); + PermutArrayType input( inputIndices, interArray ); + + // Make an output array + InterArrayType output; + + WaveletDWT::IDWT1D( input, L1d, output ); + + // Move output to intermediate array + vtkm::cont::ArrayPortalToIterators< InterPortalType > + outputIter( output.GetPortalControl() ); + vtkm::cont::ArrayPortalToIterators< InterPortalType > + interArrayIter( interArray.GetPortalControl() ); + std::copy( outputIter.GetBegin(), outputIter.GetEnd(), interArrayIter.GetBegin() ); + + L1d[0] = L1d[2]; + L1d[1] = L[i+1]; + } + + vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Copy + ( interArray, sigOut ); + + #undef VAL + + + return 0; + } + + // Squash coefficients smaller than a threshold + /* + template< typename CoeffArrayType > + VTKM_EXEC_CONT_EXPORT + vtkm::Id SquashCoefficients( CoeffArrayType, &coeffIn, + vtkm::Id ratio ); + */ + + + // Compute the book keeping array L for 1D wavelet decomposition + void ComputeL( vtkm::Id sigInLen, vtkm::Id nLevels, vtkm::Id* L ) + { + L[nLevels+1] = sigInLen; + L[nLevels] = sigInLen; + for( vtkm::Id i = nLevels; i > 0; i-- ) + { + L[i-1] = WaveletBase::GetApproxLength( L[i] ); + L[i] = WaveletBase::GetDetailLength( L[i] ); + } + } + + // Compute the length of coefficients + vtkm::Id ComputeCoeffLength( const vtkm::Id* L, vtkm::Id nLevels ) + { + vtkm::Id sum = L[0]; // 1st level cA + for( vtkm::Id i = 1; i <= nLevels; i++ ) + sum += L[i]; + return sum; + } + + // Compute approximate coefficient length at a specific level + vtkm::Id GetApproxLengthLevN( vtkm::Id sigInLen, vtkm::Id levN ) + { + vtkm::Id cALen = sigInLen; + for( vtkm::Id i = 0; i < levN; i++ ) + { + cALen = WaveletBase::GetApproxLength( cALen ); + if( cALen == 0 ) + return cALen; + } + + return cALen; + } + + + +}; // class WaveletCompressor + +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/testing/CMakeLists.txt b/vtkm/worklet/testing/CMakeLists.txt index 886c0b051..9bdfea4f7 100644 --- a/vtkm/worklet/testing/CMakeLists.txt +++ b/vtkm/worklet/testing/CMakeLists.txt @@ -38,7 +38,7 @@ set(unit_tests UnitTestWorkletMapTopologyExplicit.cxx UnitTestWorkletMapTopologyUniform.cxx UnitTestVertexClustering.cxx - UnitTestWaveletTransforms.cxx + UnitTestWaveletCompressor.cxx ) vtkm_save_worklet_unit_tests(${unit_tests}) diff --git a/vtkm/worklet/testing/UnitTestWaveletCompressor.cxx b/vtkm/worklet/testing/UnitTestWaveletCompressor.cxx new file mode 100644 index 000000000..b866a52dd --- /dev/null +++ b/vtkm/worklet/testing/UnitTestWaveletCompressor.cxx @@ -0,0 +1,208 @@ +//============================================================================ +// 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 + +#include +#include +#include + +#include + + + +void TestExtend1D() +{ + // make input data array handle + typedef vtkm::Float64 T; + typedef vtkm::cont::ArrayHandle ArrayType; + vtkm::Id sigLen = 20; + std::vector tmpVector; + for( vtkm::Id i = 0; i < sigLen; i++ ) + tmpVector.push_back( static_cast(i) ); + + vtkm::cont::ArrayHandle inputArray = + vtkm::cont::make_ArrayHandle(tmpVector); + ArrayType outputArray; + + vtkm::worklet::wavelets::WaveletDWT w("CDF9/7"); + /* + typedef vtkm::cont::ArrayHandleConcatenate< ArrayType, ArrayType> + ArrayConcat; + typedef vtkm::cont::ArrayHandleConcatenate< ArrayConcat, ArrayType > ArrayConcat2; + ArrayConcat2 outputArray; + */ + w.Extend1D( inputArray, outputArray, 4, + vtkm::worklet::wavelets::SYMW, vtkm::worklet::wavelets::SYMW ); + + std::cout << "Start testing Extend1D" << std::endl; + for (vtkm::Id i = 0; i < outputArray.GetNumberOfValues(); ++i) + std::cout << outputArray.GetPortalConstControl().Get(i) << std::endl; + std::cout << "\nFinish testing Extend1D" << std::endl; +} + +VTKM_EXEC_CONT_EXPORT +void TestDWTIDWT1D() +{ + vtkm::Id sigLen = 20; + std::cout << "Testing Wavelets Worklet" << std::endl; + std::cout << "Default test size is 20. " << std::endl; + std::cout << "Input a new size to test." << std::endl; + std::cout << "Input 0 to stick with 20." << std::endl; + vtkm::Id tmpIn; + vtkm::Id million = 1;//1000000; + std::cin >> tmpIn; + if( tmpIn != 0 ) + sigLen = tmpIn * million; + + // make input data array handle + std::vector tmpVector; + for( vtkm::Id i = 0; i < sigLen; i++ ) + tmpVector.push_back( static_cast( i ) ); + vtkm::cont::ArrayHandle inputArray = + vtkm::cont::make_ArrayHandle(tmpVector); + + vtkm::cont::ArrayHandle coeffOut; + vtkm::Id L[3]; + + // Forward Transform + vtkm::worklet::wavelets::WaveletDWT waveletdwt( "CDF9/7" ); + waveletdwt.DWT1D( inputArray, coeffOut, L ); + + std::cout << "Forward Wavelet Transform: result coeff length = " << + coeffOut.GetNumberOfValues() << std::endl; + /* + printf("L[0] = %lld, L[1] = %lld, L[2] = %lld\n", L[0], L[1], L[2] ); + */ + for( vtkm::Id i; i < coeffOut.GetNumberOfValues(); i++ ) + { + 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 reconstructArray; + waveletdwt.IDWT1D( coeffOut, L, reconstructArray ); + std::cout << "Inverse Wavelet Transform: result signal length = " << + reconstructArray.GetNumberOfValues() << std::endl; + for( vtkm::Id i; i < reconstructArray.GetNumberOfValues(); i++ ) + { + std::cout << reconstructArray.GetPortalConstControl().Get(i) << std::endl; + } +} + +VTKM_EXEC_CONT_EXPORT +void TestWaveDecomposeReconstruct() +{ + vtkm::Id sigLen = 20; + std::cout << "Testing Wavelets Worklet" << std::endl; + std::cout << "Default test size is 20. " << std::endl; + std::cout << "Input a new size to test." << std::endl; + std::cout << "Input 0 to stick with 20." << std::endl; + vtkm::Id tmpIn; + vtkm::Id million = 1000000; + std::cin >> tmpIn; + if( tmpIn != 0 ) + sigLen = tmpIn * million; + + // make input data array handle + std::vector tmpVector; + for( vtkm::Id i = 0; i < sigLen; i++ ) + tmpVector.push_back( vtkm::Sin(static_cast( i ) )); + vtkm::cont::ArrayHandle inputArray = + vtkm::cont::make_ArrayHandle(tmpVector); + + vtkm::cont::ArrayHandle outputArray; + + // Use a WaveletCompressor + vtkm::Id nLevels = 2; + vtkm::worklet::WaveletCompressor compressor("CDF9/7"); + + // User input of decompose levels + vtkm::Id maxLevel = compressor.GetWaveletMaxLevel( sigLen ); + std::cout << "Please input how many wavelet transform levels to perform, between 1 and " + << maxLevel << std::endl; + vtkm::Id levTemp; + std::cin >> levTemp; + if( levTemp > 0 && levTemp <= maxLevel ) + nLevels = levTemp; + else + { + std::cerr << "not valid levels of transforms" << std::endl; + exit(1); + } + vtkm::Id L[ nLevels + 2 ]; + + // Use a timer and decompose + vtkm::cont::Timer<> timer; + compressor.WaveDecompose( inputArray, nLevels, outputArray, L ); + vtkm::Float64 elapsedTime = timer.GetElapsedTime(); + std::cout << "Decompose takes time: " << elapsedTime << std::endl; + + #if 0 + std::cout << "Output coefficients has length = " << + outputArray.GetNumberOfValues() << std::endl; + for( vtkm::Id i = 0; i < outputArray.GetNumberOfValues(); i++ ) + { + std::cout << outputArray.GetPortalConstControl().Get(i) << std::endl; + } + #endif + + // Sort all coefficients + + + // Reconstruct + vtkm::cont::ArrayHandle reconstructArray; + timer.Reset(); + compressor.WaveReconstruct( outputArray, nLevels, L, reconstructArray ); + elapsedTime = timer.GetElapsedTime(); + std::cout << "Reconstruction takes time: " << elapsedTime << std::endl; + + //std::cout << "Reconstruct array has length = " << + // reconstructArray.GetNumberOfValues() << std::endl; + timer.Reset(); + for( vtkm::Id i = 0; i < reconstructArray.GetNumberOfValues(); i++ ) + { + //std::cout << reconstructArray.GetPortalConstControl().Get(i) << std::endl; + + VTKM_TEST_ASSERT( test_equal( reconstructArray.GetPortalConstControl().Get(i), + vtkm::Sin( static_cast(i) )), + "output value not the same..." ); + } + elapsedTime = timer.GetElapsedTime(); + std::cout << "Verification takes time: " << elapsedTime << std::endl; + +} + +void TestWaveletCompressor() +{ + std::cout << "Welcome to WaveletCompressor test program :) " << std::endl; + //TestExtend1D(); + //TestDWTIDWT1D(); + TestWaveDecomposeReconstruct(); +} + +int UnitTestWaveletCompressor(int, char *[]) +{ + return vtkm::cont::testing::Testing::Run(TestWaveletCompressor); +} diff --git a/vtkm/worklet/wavelets/CMakeLists.txt b/vtkm/worklet/wavelets/CMakeLists.txt new file mode 100644 index 000000000..d1758b5df --- /dev/null +++ b/vtkm/worklet/wavelets/CMakeLists.txt @@ -0,0 +1,30 @@ +##============================================================================ +## 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 2016 Sandia Corporation. +## Copyright 2016 UT-Battelle, LLC. +## Copyright 2016 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. +##============================================================================ + +set(headers + FilterBanks.h + WaveletFilter.h + WaveletBase.h + WaveletTransforms.h + WaveletDWT.h + ) + +#----------------------------------------------------------------------------- +vtkm_declare_headers(${headers}) diff --git a/vtkm/worklet/wavelets/FilterBanks.h b/vtkm/worklet/wavelets/FilterBanks.h new file mode 100644 index 000000000..4814174d9 --- /dev/null +++ b/vtkm/worklet/wavelets/FilterBanks.h @@ -0,0 +1,80 @@ +//============================================================================= +// +// 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 2015 Sandia Corporation. +// Copyright 2015 UT-Battelle, LLC. +// Copyright 2015 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. +// +//============================================================================= + +#ifndef vtk_m_worklet_wavelets_filterbanks_h +#define vtk_m_worklet_wavelets_filterbanks_h + + +#include + +namespace vtkm { +namespace worklet { + +namespace wavelets { + + const vtkm::Float64 hm4_44[9] = { + 0.037828455507264, + -0.023849465019557, + -0.110624404418437, + 0.377402855612831, + 0.852698679008894, + 0.377402855612831, + -0.110624404418437, + -0.023849465019557, + 0.037828455507264 + }; + + const vtkm::Float64 h4[9] = { + 0.0, + -0.064538882628697, + -0.040689417609164, + 0.418092273221617, + 0.788485616405583, + 0.418092273221617, + -0.0406894176091641, + -0.0645388826286971, + 0.0 + }; + + const double hm2_22[6] = { + -0.1767766952966368811002110905262, + 0.3535533905932737622004221810524, + 1.0606601717798212866012665431573, + 0.3535533905932737622004221810524, + -0.1767766952966368811002110905262 + }; + + const double h2[18] = { + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.3535533905932737622004221810524, + 0.7071067811865475244008443621048, + 0.3535533905932737622004221810524, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0 + }; + +}; + +} +} + +#endif diff --git a/vtkm/worklet/wavelets/WaveletBase.h b/vtkm/worklet/wavelets/WaveletBase.h new file mode 100644 index 000000000..f1f057c89 --- /dev/null +++ b/vtkm/worklet/wavelets/WaveletBase.h @@ -0,0 +1,183 @@ +//============================================================================ +// 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. +//============================================================================ + +#ifndef vtk_m_worklet_wavelets_waveletbase_h +#define vtk_m_worklet_wavelets_waveletbase_h + + +#include + +#include + +namespace vtkm { +namespace worklet { + +namespace wavelets { + +enum DWTMode { // boundary extension modes + INVALID = -1, + ZPD, + SYMH, + SYMW, + ASYMH, ASYMW, SP0, SP1, PPD, PER +}; + +// Functionalities are similar to MatWaveBase in VAPoR. +class WaveletBase +{ +public: + + // Constructor + WaveletBase( const std::string &w_name ) + { + this->filter = NULL; + this->wmode = PER; + this->wname = w_name; + if( wname.compare("CDF9/7") == 0 ) + { + this->wmode = SYMW; // Default extension mode, see MatWaveBase.cpp + this->filter = new vtkm::worklet::wavelets::WaveletFilter( wname ); + } + else if( wname.compare("CDF5/3") == 0 ) + { + this->wmode = SYMW; + this->filter = new vtkm::worklet::wavelets::WaveletFilter( wname ); + } + else + { + std::cerr << "This wavelet kernel is not supported: " << wname << std::endl; + // throw an error + } + } + + // Destructor + virtual ~WaveletBase() + { + if( filter ) + delete filter; + filter = NULL; + } + + // Get the wavelet filter + const WaveletFilter* GetWaveletFilter() + { + if( this->filter == NULL ) + { + // throw an error + } + return filter; + } + + // Returns length of approximation coefficients from a decompostition pass. + vtkm::Id GetApproxLength( vtkm::Id sigInLen ) + { + vtkm::Id filterLen = this->filter->GetFilterLength(); + + if (this->wmode == PER) + return static_cast(vtkm::Ceil( (static_cast(sigInLen)) / 2.0 )); + else if (this->filter->isSymmetric()) + { + if ( (this->wmode == SYMW && (filterLen % 2 != 0)) || + (this->wmode == SYMH && (filterLen % 2 == 0)) ) + { + if (sigInLen % 2 != 0) + return((sigInLen+1) / 2); + else + return((sigInLen) / 2); + } + } + + return static_cast( vtkm::Floor( + static_cast(sigInLen + filterLen - 1) / 2.0 ) ); + } + + // Returns length of detail coefficients from a decompostition pass + vtkm::Id GetDetailLength( vtkm::Id sigInLen ) + { + vtkm::Id filterLen = this->filter->GetFilterLength(); + + if (this->wmode == PER) + return static_cast(vtkm::Ceil( (static_cast(sigInLen)) / 2.0 )); + else if (this->filter->isSymmetric()) + { + if ( (this->wmode == SYMW && (filterLen % 2 != 0)) || + (this->wmode == SYMH && (filterLen % 2 == 0)) ) + { + if (sigInLen % 2 != 0) + return((sigInLen-1) / 2); + else + return((sigInLen) / 2); + } + } + + return static_cast( vtkm::Floor( + static_cast(sigInLen + filterLen - 1) / 2.0 ) ); + } + + // Returns length of coefficients generated in a decompostition pass + vtkm::Id GetCoeffLength( vtkm::Id sigInLen ) + { + return( GetApproxLength( sigInLen ) + GetDetailLength( sigInLen ) ); + } + vtkm::Id GetCoeffLength2( vtkm::Id sigInX, vtkm::Id sigInY ) + { + return( GetCoeffLength( sigInX) * GetCoeffLength( sigInY ) ); + } + vtkm::Id GetCoeffLength3( vtkm::Id sigInX, vtkm::Id sigInY, vtkm::Id sigInZ) + { + return( GetCoeffLength( sigInX) * GetCoeffLength( sigInY ) * GetCoeffLength( sigInZ ) ); + } + + // Returns maximum wavelet decompostion level + vtkm::Id GetWaveletMaxLevel( vtkm::Id sigInLen ) + { + if( ! this->filter ) + return 0; + else { + vtkm::Id filterLen = this->filter->GetFilterLength(); + vtkm::Id level; + this->WaveLengthValidate( sigInLen, filterLen, level ); + return level; + } + } + +protected: + vtkm::worklet::wavelets::DWTMode wmode; + WaveletFilter* filter; + std::string wname; + + void WaveLengthValidate( vtkm::Id sigInLen, vtkm::Id filterLength, vtkm::Id &level) + { + if( sigInLen < filterLength ) + level = 0; + else + level = static_cast( vtkm::Floor( + vtkm::Log2( static_cast(sigInLen) / + static_cast(filterLength) ) + 1.0 ) ); + } +}; // class WaveletBase. + + +} // namespace wavelets + +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/wavelets/WaveletDWT.h b/vtkm/worklet/wavelets/WaveletDWT.h new file mode 100644 index 000000000..12f837e73 --- /dev/null +++ b/vtkm/worklet/wavelets/WaveletDWT.h @@ -0,0 +1,407 @@ +//============================================================================ +// 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. +//============================================================================ + +#ifndef vtk_m_worklet_wavelets_waveletdwt_h +#define vtk_m_worklet_wavelets_waveletdwt_h + +#include + +#include + +#include +#include +#include + +#include + +namespace vtkm { +namespace worklet { +namespace wavelets { + +class WaveletDWT : public WaveletBase +{ +public: + + // Constructor + WaveletDWT( const std::string &w_name ) : WaveletBase( w_name ) {} + + + // Func: Extend 1D signal + template< typename SigInArrayType, typename SigExtendedArrayType > + vtkm::Id Extend1D( const SigInArrayType &sigIn, // Input + SigExtendedArrayType &sigOut, // Output + vtkm::Id addLen, + vtkm::worklet::wavelets::DWTMode leftExtMethod, + vtkm::worklet::wavelets::DWTMode rightExtMethod ) + { + typedef typename SigInArrayType::ValueType ValueType; + typedef vtkm::cont::ArrayHandle< ValueType > ExtensionArrayType; + ExtensionArrayType leftExtend, rightExtend; + leftExtend.Allocate( addLen ); + rightExtend.Allocate( addLen ); + + typedef typename ExtensionArrayType::PortalControl PortalType; + typedef typename SigInArrayType::PortalConstControl SigInPortalConstType; + + typedef vtkm::cont::ArrayHandleConcatenate< ExtensionArrayType, SigInArrayType> + ArrayConcat; + typedef vtkm::cont::ArrayHandleConcatenate< ArrayConcat, ExtensionArrayType > + ArrayConcat2; + + PortalType leftExtendPortal = leftExtend.GetPortalControl(); + PortalType rightExtendPortal = rightExtend.GetPortalControl(); + SigInPortalConstType sigInPortal = sigIn.GetPortalConstControl(); + vtkm::Id sigInLen = sigIn.GetNumberOfValues(); + + switch( leftExtMethod ) + { + case vtkm::worklet::wavelets::SYMH: + { + for( vtkm::Id count = 0; count < addLen; count++ ) + leftExtendPortal.Set( count, sigInPortal.Get( addLen - count - 1) ); + break; + } + case vtkm::worklet::wavelets::SYMW: + { + for( vtkm::Id count = 0; count < addLen; count++ ) + leftExtendPortal.Set( count, sigInPortal.Get( addLen - count ) ); + break; + } + default: + { + // throw out an error + return 1; + } + } + + switch( rightExtMethod ) + { + case vtkm::worklet::wavelets::SYMH: + { + for( vtkm::Id count = 0; count < addLen; count++ ) + rightExtendPortal.Set( count, sigInPortal.Get( sigInLen - count - 1 ) ); + break; + } + case SYMW: + { + for( vtkm::Id count = 0; count < addLen; count++ ) + rightExtendPortal.Set( count, sigInPortal.Get( sigInLen - count - 2 ) ); + break; + } + default: + { + // throw out an error + return 1; + } + } + + ArrayConcat leftOn( leftExtend, sigIn ); + ArrayConcat2 rightOn(leftOn, rightExtend ); + vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Copy + ( rightOn, sigOut ); + + return 0; + } + + + + // Func: + // Performs one level of 1D discrete wavelet transform + // It takes care of boundary conditions, etc. + template< typename SignalArrayType, typename CoeffArrayType> + VTKM_EXEC_CONT_EXPORT + vtkm::Id DWT1D( const SignalArrayType &sigIn, // Input + CoeffArrayType &coeffOut, // Output: cA followed by cD + vtkm::Id L[3] ) // Output: how many cA and cD. + { + + vtkm::Id sigInLen = sigIn.GetNumberOfValues(); + if( GetWaveletMaxLevel( sigInLen ) < 1 ) + { + // throw an error + std::cerr << "Cannot transform signal of length " << sigInLen << std::endl; + return -1; + } + + L[0] = WaveletBase::GetApproxLength( sigInLen ); + L[1] = WaveletBase::GetDetailLength( sigInLen ); + L[2] = sigInLen; + + VTKM_ASSERT( L[0] + L[1] == L[2] ); + + vtkm::Id filterLen = WaveletBase::filter->GetFilterLength(); + + bool doSymConv = false; + if( WaveletBase::filter->isSymmetric() ) + { + 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; + + typedef typename SignalArrayType::ValueType SigInValueType; + typedef vtkm::cont::ArrayHandle SignalArrayTypeBasic; + + SignalArrayTypeBasic sigInExtended; + + this->Extend1D( sigIn, sigInExtended, addLen, WaveletBase::wmode, WaveletBase::wmode ); + VTKM_ASSERT( sigInExtended.GetNumberOfValues() == sigExtendedLen ); + + // Coefficients in coeffOutTmp are interleaving, + // e.g. cA are at 0, 2, 4...; cD are at 1, 3, 5... + typedef typename CoeffArrayType::ValueType CoeffOutValueType; + typedef vtkm::cont::ArrayHandle< CoeffOutValueType > CoeffOutArrayBasic; + CoeffOutArrayBasic coeffOutTmp; + + + // initialize a worklet + vtkm::worklet::wavelets::ForwardTransform forwardTransform; + forwardTransform.SetFilterLength( filterLen ); + forwardTransform.SetCoeffLength( L[0], L[1] ); + forwardTransform.SetOddness( oddLow, oddHigh ); + + vtkm::worklet::DispatcherMapField + dispatcher(forwardTransform); + dispatcher.Invoke( sigInExtended, + WaveletBase::filter->GetLowDecomposeFilter(), + WaveletBase::filter->GetHighDecomposeFilter(), + coeffOutTmp ); + + // Separate cA and cD. + typedef vtkm::cont::ArrayHandleCounting< vtkm::Id > IdArrayType; + typedef vtkm::cont::ArrayHandlePermutation< IdArrayType, CoeffArrayType > PermutArrayType; + + IdArrayType approxIndices( 0, 2, L[0] ); + IdArrayType detailIndices( 1, 2, L[1] ); + PermutArrayType cATmp( approxIndices, coeffOutTmp ); + PermutArrayType cDTmp( detailIndices, coeffOutTmp ); + + typedef vtkm::cont::ArrayHandleConcatenate< PermutArrayType, PermutArrayType > + PermutArrayConcatenated; + PermutArrayConcatenated coeffOutConcat( cATmp, cDTmp ); + + vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Copy( + coeffOutConcat, coeffOut ); + + return 0; + } + + + + // Func: + // Performs one level of inverse wavelet transform + // It takes care of boundary conditions, etc. + template< typename CoeffArrayType, typename SignalArrayType> + VTKM_EXEC_CONT_EXPORT + vtkm::Id IDWT1D( const CoeffArrayType &coeffIn, // Input, cA followed by cD + const vtkm::Id L[3], // Input, how many cA and cD + SignalArrayType &sigOut ) // Output + { + #if 0 + std::cerr << "coeffIn len = " << coeffIn.GetNumberOfValues() << std::endl; + std::cerr << L[0] << ", " << L[1] << ", " << L[2] << std::endl; + #endif + + VTKM_ASSERT( coeffIn.GetNumberOfValues() == L[2] ); + + vtkm::Id filterLen = WaveletBase::filter->GetFilterLength(); + 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; + + if( WaveletBase::filter->isSymmetric() ) + { + if(( WaveletBase::wmode == SYMW && (filterLen % 2 != 0) ) || + ( WaveletBase::wmode == SYMH && (filterLen % 2 == 0) ) ) + { + doSymConv = true; + + 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; + } + } + } + + vtkm::Id cATempLen, cDTempLen, reconTempLen; + vtkm::Id addLen = 0; + vtkm::Id cDPadLen = 0; + if( doSymConv ) // extend cA and cD + { + addLen = filterLen / 4; + if( (L[0] > L[1]) && (WaveletBase::wmode == SYMH) ) + cDPadLen = L[0]; + cATempLen = L[0] + 2 * addLen; + cDTempLen = cATempLen; // same length + } + else // not extend cA and cD + { + cATempLen = L[0]; + cDTempLen = L[1]; + } + reconTempLen = L[2]; + if( reconTempLen % 2 != 0 ) + reconTempLen++; + + typedef vtkm::cont::ArrayHandleCounting< vtkm::Id > IdArrayType; + typedef vtkm::cont::ArrayHandlePermutation< IdArrayType, CoeffArrayType > PermutArrayType; + + // 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 > CoeffArrayBasic; + typedef vtkm::cont::ArrayHandle< CoeffValueType > ExtensionArrayType; + + CoeffArrayBasic cATemp, cDTemp; + + + if( doSymConv ) // Actually extend cA and cD + { + this->Extend1D( cA, cATemp, addLen, cALeftMode, cARightMode ); + if( cDPadLen > 0 ) + { + // Add back the missing final cD: 0.0 + ExtensionArrayType singleValArray; + singleValArray.Allocate(1); + singleValArray.GetPortalControl().Set(0, 0.0); + vtkm::cont::ArrayHandleConcatenate< PermutArrayType, ExtensionArrayType > + cDPad( cD, singleValArray ); + + this->Extend1D( cDPad, cDTemp, addLen, cDLeftMode, cDRightMode ); + } + else + { + this->Extend1D( cD, cDTemp, addLen, cDLeftMode, cDRightMode ); + + // Attached an zero if cDTemp is shorter than cDTempLen + if( cDTemp.GetNumberOfValues() != cDTempLen ) + { + VTKM_ASSERT( cDTemp.GetNumberOfValues() == cDTempLen - 1 ); + CoeffArrayBasic singleValArray; + singleValArray.Allocate(1); + singleValArray.GetPortalControl().Set(0, 0.0); + vtkm::cont::ArrayHandleConcatenate< CoeffArrayBasic, CoeffArrayBasic> + concat1( cDTemp, singleValArray ); + CoeffArrayBasic cDStorage; + vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG >::Copy + ( concat1, cDStorage ); + cDTemp = cDStorage; + } + } + } // end if( doSymConv ) + else // Make cATemp and cDTemp from cA and cD + { + vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG >::Copy + (cA, cATemp ); + vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG >::Copy + (cD, cDTemp ); + } + + #if 0 + std::cerr << "cATemp has length: " << cATemp.GetNumberOfValues() << std::endl; + for( vtkm::Id i = 0; i < cATemp.GetNumberOfValues(); i++ ) + std::cout << cATemp.GetPortalConstControl().Get(i) << std::endl; + std::cerr << "cDTemp has length: " << cDTemp.GetNumberOfValues() << std::endl; + for( vtkm::Id i = 0; i < cDTemp.GetNumberOfValues(); i++ ) + std::cout << cDTemp.GetPortalConstControl().Get(i) << std::endl; + #endif + + if( filterLen % 2 != 0 ) + { + vtkm::cont::ArrayHandleConcatenate< CoeffArrayBasic, CoeffArrayBasic> + coeffInExtended( cATemp, cDTemp ); + + // Initialize a worklet + vtkm::worklet::wavelets::InverseTransformOdd inverseXformOdd; + inverseXformOdd.SetFilterLength( filterLen ); + inverseXformOdd.SetCALength( L[0], cATempLen ); + vtkm::worklet::DispatcherMapField< vtkm::worklet::wavelets::InverseTransformOdd > + dispatcher( inverseXformOdd ); + dispatcher.Invoke( coeffInExtended, + WaveletBase::filter->GetLowReconstructFilter(), + WaveletBase::filter->GetHighReconstructFilter(), + sigOut ); + #if 0 + std::cerr << "coeffInExtended len = " << coeffInExtended.GetNumberOfValues() << std::endl; + std::cerr << sigOut.GetNumberOfValues() << ", " << L[2] << std::endl; + #endif + + VTKM_ASSERT( sigOut.GetNumberOfValues() >= L[2] ); + sigOut.Shrink( L[2] ); + } + else + { + // need to implement the even filter length worklet first + } + + return 0; + + } // function IDWT1D + + +}; // class WaveletDWT + +} // namespace wavelets +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/wavelets/WaveletFilter.h b/vtkm/worklet/wavelets/WaveletFilter.h new file mode 100644 index 000000000..32bea1fd4 --- /dev/null +++ b/vtkm/worklet/wavelets/WaveletFilter.h @@ -0,0 +1,182 @@ +//============================================================================ +// 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. +//============================================================================ + +#ifndef vtk_m_worklet_wavelets_waveletfilter_h +#define vtk_m_worklet_wavelets_waveletfilter_h + +#include + +#include + +#include + +namespace vtkm { +namespace worklet { + +namespace wavelets { + +// Wavelet filter class; +// functionally equivalent to WaveFiltBase and its subclasses in VAPoR. +class WaveletFilter +{ +public: + // constructor + WaveletFilter( const std::string &wname ) + { + lowDecomposeFilter = highDecomposeFilter = + lowReconstructFilter = highReconstructFilter = NULL; + this->filterLength = 0; + if( wname.compare("CDF9/7") == 0 ) + { + this->symmetricity= true; + this->filterLength = 9; + this->AllocateFilterMemory(); + wrev( vtkm::worklet::wavelets::hm4_44, lowDecomposeFilter, filterLength ); + qmf_wrev( vtkm::worklet::wavelets::h4, highDecomposeFilter, filterLength ); + verbatim_copy( vtkm::worklet::wavelets::h4, lowReconstructFilter, filterLength ); + qmf_even( vtkm::worklet::wavelets::hm4_44, highReconstructFilter, filterLength ); + } + else if( wname.compare("CDF5/3") == 0 ) + { + this->symmetricity = true; + this->filterLength = 5; + this->AllocateFilterMemory(); + wrev( vtkm::worklet::wavelets::hm2_22, lowDecomposeFilter, filterLength ); + qmf_wrev( vtkm::worklet::wavelets::h2+6, highDecomposeFilter, filterLength ); + verbatim_copy( vtkm::worklet::wavelets::h2+6, lowReconstructFilter, filterLength ); + qmf_even( vtkm::worklet::wavelets::hm2_22, highReconstructFilter, filterLength ); + } + else + { + std::cerr << "Not supported wavelet kernel: " << wname << std::endl; + // throw an error here + } + } + + // destructor + virtual ~WaveletFilter() + { + if( lowDecomposeFilter ) delete[] lowDecomposeFilter; + if( highDecomposeFilter ) delete[] highDecomposeFilter; + if( lowReconstructFilter ) delete[] lowReconstructFilter; + if( highReconstructFilter ) delete[] highReconstructFilter; + } + + vtkm::Id GetFilterLength() { return this->filterLength; } + bool isSymmetric() { return this->symmetricity; } + + typedef vtkm::cont::ArrayHandle FilterType; + FilterType GetLowDecomposeFilter() const + { + return vtkm::cont::make_ArrayHandle( lowDecomposeFilter, filterLength ); + } + FilterType GetHighDecomposeFilter() const + { + return vtkm::cont::make_ArrayHandle( highDecomposeFilter, filterLength ); + } + FilterType GetLowReconstructFilter() const + { + return vtkm::cont::make_ArrayHandle( lowReconstructFilter, filterLength); + } + FilterType GetHighReconstructFilter() const + { + return vtkm::cont::make_ArrayHandle( highReconstructFilter, filterLength ); + } + +protected: + bool symmetricity; + vtkm::Id filterLength; + vtkm::Float64* lowDecomposeFilter; + vtkm::Float64* highDecomposeFilter; + vtkm::Float64* lowReconstructFilter; + vtkm::Float64* highReconstructFilter; + + void AllocateFilterMemory() + { + lowDecomposeFilter = new vtkm::Float64[ this->filterLength ]; + highDecomposeFilter = new vtkm::Float64[ this->filterLength ]; + lowReconstructFilter = new vtkm::Float64[ this->filterLength ]; + highReconstructFilter = new vtkm::Float64[ this->filterLength ]; + } + + // Flipping operation; helper function to initialize a filter. + void wrev( const vtkm::Float64* sigIn, vtkm::Float64* sigOut, vtkm::Id sigLength ) + { + for( vtkm::Id count = 0; count < sigLength; count++) + sigOut[count] = sigIn[sigLength - count - 1]; + } + + // Quadrature mirror filtering operation: helper function to initialize a filter. + void qmf_even ( const vtkm::Float64* sigIn, vtkm::Float64* sigOut, vtkm::Id sigLength ) + { + for (vtkm::Id count = 0; count < sigLength; count++) + { + sigOut[count] = sigIn[sigLength - count - 1]; + + if (sigLength % 2 == 0) { + if (count % 2 != 0) + sigOut[count] = -1.0 * sigOut[count]; + } + else { + if (count % 2 == 0) + sigOut[count] = -1.0 * sigOut[count]; + } + } + } + + // Flipping and QMF at the same time: helper function to initialize a filter. + void qmf_wrev ( const vtkm::Float64* sigIn, vtkm::Float64* sigOut, vtkm::Id sigLength ) + { + for (vtkm::Id count = 0; count < sigLength; count++) { + sigOut[count] = sigIn[sigLength - count - 1]; + + if (sigLength % 2 == 0) { + if (count % 2 != 0) + sigOut[count] = -1 * sigOut[count]; + } + else { + if (count % 2 == 0) + sigOut[count] = -1 * sigOut[count]; + } + } + + vtkm::Float64 tmp; + for (vtkm::Id count = 0; count < sigLength/2; count++) { + tmp = sigOut[count]; + sigOut[count] = sigOut[sigLength - count - 1]; + sigOut[sigLength - count - 1] = tmp; + } + } + + // Verbatim Copying: helper function to initialize a filter. + void verbatim_copy ( const vtkm::Float64* sigIn, vtkm::Float64* sigOut, vtkm::Id sigLength ) + { + for (vtkm::Id count = 0; count < sigLength; count++) + sigOut[count] = sigIn[count]; + } +}; // class WaveletFilter. + +} // namespace wavelets. + + +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/worklet/WaveletTransforms.h b/vtkm/worklet/wavelets/WaveletTransforms.h similarity index 98% rename from vtkm/worklet/WaveletTransforms.h rename to vtkm/worklet/wavelets/WaveletTransforms.h index 2a2c53e97..df77e8d77 100644 --- a/vtkm/worklet/WaveletTransforms.h +++ b/vtkm/worklet/wavelets/WaveletTransforms.h @@ -28,6 +28,7 @@ namespace vtkm { namespace worklet { +namespace wavelets { // Worklet: perform a simple forward transform class ForwardTransform: public vtkm::worklet::WorkletMapField @@ -227,9 +228,10 @@ private: vtkm::Id cALen; // Number of actual cAs vtkm::Id cALenExtended; // Number of extended cA at the beginning of input array -}; // Finish class ForwardTransform +}; // class ForwardTransform -} // Finish namespace worlet -} // Finish namespace vtkm +} // namespace wavelets +} // namespace worlet +} // namespace vtkm #endif // vtk_m_worklet_Wavelets_h