diff --git a/vtkm/filter/WaveletCompressor.h b/vtkm/filter/WaveletCompressor.h deleted file mode 100644 index 45a60c278..000000000 --- a/vtkm/filter/WaveletCompressor.h +++ /dev/null @@ -1,88 +0,0 @@ -//============================================================================ -// 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_filter_waveletcompressor_h -#define vtk_m_filter_waveletcompressor_h - -#include - -//#include - -#include -#include -#include - -#include - -namespace vtkm { -namespace filter { - -//template< typename DeviceAdapter > -class WaveletCompressor : public internal::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 ); - - // 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 ); - - // 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 ); - - // Compute the length of coefficients - vtkm::Id ComputeCoeffLength( const vtkm::Id* L, vtkm::Id nLevels ); - - // Compute approximate coefficient length at a specific level - vtkm::Id GetApproxLengthLevN( vtkm::Id sigInLen, vtkm::Id levN ); - - - -}; // Finish class WaveletCompressor - -} // Finish namespace filter -} // Finish namespace vtkm - -#include - -#endif diff --git a/vtkm/filter/WaveletCompressor.hxx b/vtkm/filter/WaveletCompressor.hxx deleted file mode 100644 index 304d72438..000000000 --- a/vtkm/filter/WaveletCompressor.hxx +++ /dev/null @@ -1,225 +0,0 @@ -//============================================================================ -// 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 - -namespace vtkm { -namespace filter { - - -template< typename SignalArrayType, typename CoeffArrayType> -VTKM_EXEC_CONT_EXPORT -vtkm::Id -WaveletCompressor::WaveDecompose( const SignalArrayType &sigIn, // Input - vtkm::Id nLevels, // n levels of DWT - CoeffArrayType &C, - vtkm::Id* L ) // bookkeeping array; - // len(L) = nLevels+2 -{ - vtkm::Id sigInLen = sigIn.GetNumberOfValues(); - if( nLevels < 1 || nLevels > WaveletBase::GetWaveletMaxLevel( sigInLen ) ) - { - std::cerr << "nLevel is not supported: " << nLevels << std::endl; - // throw an error - } - /* - if( nLevels == 0 ) // 0 levels means no transform - { - 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, C ); - - #undef VAL - - return 0; -} - - - -template< typename CoeffArrayType, typename SignalArrayType > -VTKM_EXEC_CONT_EXPORT -vtkm::Id -WaveletCompressor::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; -} - -/* -template< typename CoeffArrayType > -VTKM_EXEC_CONT_EXPORT -vtkm::Id -WaveletCompressor::SquashCoefficients( CoeffArrayType, &coeffIn, - vtkm::Id ratio ) -{ - typedef typename CoeffArrayType::ValueType ValueType; - vtkm::cont::ArrayHandle< ValueType > sortArray; - vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Copy - ( coeffIn, sortArray ); - vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Sort( sortArray ); -} -*/ - -vtkm::Id -WaveletCompressor::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; -} - -vtkm::Id -WaveletCompressor::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; -} - -void -WaveletCompressor::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] ); - } -} - - -} // Finish namespace filter -} // Finish namespace vtkm - diff --git a/vtkm/filter/internal/FilterBanks.h b/vtkm/filter/internal/FilterBanks.h deleted file mode 100644 index 0843a3b25..000000000 --- a/vtkm/filter/internal/FilterBanks.h +++ /dev/null @@ -1,80 +0,0 @@ -//============================================================================= -// -// 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_filter_internal_filterbanks_h -#define vtk_m_filter_internal_filterbanks_h - - -#include - -namespace vtkm { -namespace filter { - -namespace internal { - - 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/filter/internal/WaveletBase.h b/vtkm/filter/internal/WaveletBase.h deleted file mode 100644 index 29debeb41..000000000 --- a/vtkm/filter/internal/WaveletBase.h +++ /dev/null @@ -1,183 +0,0 @@ -//============================================================================ -// 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_filter_internal_waveletbase_h -#define vtk_m_filter_internal_waveletbase_h - - -#include - -#include - -namespace vtkm { -namespace filter { - -namespace internal { - -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::filter::internal::WaveletFilter( wname ); - } - else if( wname.compare("CDF5/3") == 0 ) - { - this->wmode = SYMW; - this->filter = new vtkm::filter::internal::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::filter::internal::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 ) ); - } -}; // Finish class WaveletBase. - - -} // Finish namespace internal - -} // Finish namespace filter -} // Finish namespace vtkm - -#endif diff --git a/vtkm/filter/internal/WaveletDWT.h b/vtkm/filter/internal/WaveletDWT.h deleted file mode 100644 index 348556ec9..000000000 --- a/vtkm/filter/internal/WaveletDWT.h +++ /dev/null @@ -1,410 +0,0 @@ -//============================================================================ -// 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_filter_internal_waveletdwt_h -#define vtk_m_filter_internal_waveletdwt_h - -//#include -//#include - -#include - -#include - -#include -#include -#include - -#include - -namespace vtkm { -namespace filter { -namespace internal { - -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::filter::internal::DWTMode leftExtMethod, - vtkm::filter::internal::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::filter::internal::SYMH: - { - for( vtkm::Id count = 0; count < addLen; count++ ) - leftExtendPortal.Set( count, sigInPortal.Get( addLen - count - 1) ); - break; - } - case vtkm::filter::internal::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::filter::internal::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::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::filter::internal::DWTMode cALeftMode = WaveletBase::wmode; - vtkm::filter::internal::DWTMode cARightMode = WaveletBase::wmode; - vtkm::filter::internal::DWTMode cDLeftMode = WaveletBase::wmode; - vtkm::filter::internal::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::InverseTransformOdd inverseXformOdd; - inverseXformOdd.SetFilterLength( filterLen ); - inverseXformOdd.SetCALength( L[0], cATempLen ); - vtkm::worklet::DispatcherMapField< vtkm::worklet::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; - - } // Finish function IDWT1D - - -}; // Finish class WaveletDWT - -} // Finish namespace internal -} // Finish namespace filter -} // Finish namespace vtkm - -#endif diff --git a/vtkm/filter/internal/WaveletFilter.h b/vtkm/filter/internal/WaveletFilter.h deleted file mode 100644 index 6d95e77fa..000000000 --- a/vtkm/filter/internal/WaveletFilter.h +++ /dev/null @@ -1,182 +0,0 @@ -//============================================================================ -// 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_filter_internal_waveletfilter_h -#define vtk_m_filter_internal_waveletfilter_h - -#include - -#include - -#include - -namespace vtkm { -namespace filter { - -namespace internal { - -// 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::filter::internal::hm4_44, lowDecomposeFilter, filterLength ); - qmf_wrev( vtkm::filter::internal::h4, highDecomposeFilter, filterLength ); - verbatim_copy( vtkm::filter::internal::h4, lowReconstructFilter, filterLength ); - qmf_even( vtkm::filter::internal::hm4_44, highReconstructFilter, filterLength ); - } - else if( wname.compare("CDF5/3") == 0 ) - { - this->symmetricity = true; - this->filterLength = 5; - this->AllocateFilterMemory(); - wrev( vtkm::filter::internal::hm2_22, lowDecomposeFilter, filterLength ); - qmf_wrev( vtkm::filter::internal::h2+6, highDecomposeFilter, filterLength ); - verbatim_copy( vtkm::filter::internal::h2+6, lowReconstructFilter, filterLength ); - qmf_even( vtkm::filter::internal::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]; - } -}; // Finish class WaveletFilter. - -} // Finish namespace internal. - - -} // Finish namespace filter -} // Finish namespace vtkm - -#endif diff --git a/vtkm/worklet/testing/UnitTestWaveletTransforms.cxx b/vtkm/worklet/testing/UnitTestWaveletTransforms.cxx deleted file mode 100644 index 75d9eef65..000000000 --- a/vtkm/worklet/testing/UnitTestWaveletTransforms.cxx +++ /dev/null @@ -1,104 +0,0 @@ -//============================================================================ -// 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 - -#include - - -void TestWaveletTransforms( ) -{ - 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 (in millions)." << 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 + 8; i++ ) - tmpVector.push_back( static_cast(i%100+1) ); - - vtkm::cont::ArrayHandle input1DArray = - vtkm::cont::make_ArrayHandle(tmpVector); - - // output array handle - vtkm::cont::ArrayHandle outputArray1; - - // make two filter array handles - vtkm::cont::ArrayHandle lowFilter = - vtkm::cont::make_ArrayHandle(vtkm::filter::internal::hm4_44, 9); - vtkm::cont::ArrayHandle highFilter = - vtkm::cont::make_ArrayHandle(vtkm::filter::internal::h4, 9); - - // initialize a worklet - vtkm::worklet::ForwardTransform forwardTransform; - forwardTransform.SetFilterLength( 9 ); - forwardTransform.SetCoeffLength( sigLen/2, sigLen/2 ); - forwardTransform.SetOddness( false, true ); - - // setup a timer - vtkm::cont::Timer<> timer; - - vtkm::worklet::DispatcherMapField - dispatcher(forwardTransform); - dispatcher.Invoke(input1DArray, - lowFilter, - highFilter, - outputArray1); - - srand ((unsigned int)time(NULL)); - vtkm::Id randNum = rand() % sigLen; - std::cout << "A random output: " - << outputArray1.GetPortalConstControl().Get(randNum) << std::endl; - - vtkm::Float64 elapsedTime = timer.GetElapsedTime(); - std::cerr << "Dealing array size " << sigLen/million << " millions takes time " - << elapsedTime << std::endl; - if( sigLen < 21 ) - for (vtkm::Id i = 0; i < outputArray1.GetNumberOfValues(); ++i) - { - std::cout << outputArray1.GetPortalConstControl().Get(i) << ", "; - if( i % 2 != 0 ) - std::cout << std::endl; - } -} - - -int UnitTestWaveletTransforms(int, char* []) -{ - // TestDWT1D(); - // TestExtend1D(); - return vtkm::cont::testing::Testing::Run(TestWaveletTransforms); - - return 0; -}