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

merge updates from the 1D wavelets branch

Browse Source
This commit is contained in:
Samuel Li 2016-08-18 22:42:29 -06:00
commit 372cd6729b
9 changed files with 660 additions and 326 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
vtkm/cont/ArrayHandleConcatenate.h
View File

@ -142,10 +142,9 @@ public:
VTKM_CONT_EXPORT
void Allocate( vtkm::Id numberOfValues )
{
VTKM_ASSERT( numberOfValues >= 0 );
VTKM_ASSERT( this->valid );
(void)numberOfValues; // dummy statement to avoid a warning
throw vtkm::cont::ErrorControlInternal(
"ArrayHandleConcatenate should not be allocated explicitly.");
"ArrayHandleConcatenate should not be allocated explicitly. " );
}
VTKM_CONT_EXPORT
@ -242,9 +241,9 @@ public:
VTKM_CONT_EXPORT
PortalExecution PrepareForOutput( vtkm::Id numberOfValues )
{
VTKM_ASSERT( numberOfValues > 0 );
throw vtkm::cont::ErrorControlBadValue(
"ArrayHandleConcatenate is derived and read-only.");
(void)numberOfValues; // dummy statement to avoid a warning
throw vtkm::cont::ErrorControlInternal(
"ArrayHandleConcatenate is derived and read-only. " );
}
VTKM_CONT_EXPORT

2
vtkm/cont/testing/UnitTestArrayHandleConcatenate.cxx
View File

@ -56,7 +56,7 @@ void TestConcatenateEmptyArray()
{
std::vector< vtkm::Float64 > vec;
for( vtkm::Id i = 0; i < ARRAY_SIZE; i++ )
vec.push_back( i*1.5 );
vec.push_back( vtkm::Float64(i) * 1.5 );
typedef vtkm::Float64 CoeffValueType;
typedef vtkm::cont::ArrayHandle<CoeffValueType> CoeffArrayTypeTmp;

50
vtkm/worklet/WaveletCompressor.h
View File

@ -37,16 +37,17 @@ class WaveletCompressor : public vtkm::worklet::wavelets::WaveletDWT
public:
// Constructor
WaveletCompressor( const std::string &w_name ) : WaveletDWT( w_name ) {}
WaveletCompressor( wavelets::WaveletName name ) : WaveletDWT( name ) {}
// Multi-level 1D wavelet decomposition
template< typename SignalArrayType, typename CoeffArrayType>
template< typename SignalArrayType, typename CoeffArrayType, typename DeviceTag >
VTKM_CONT_EXPORT
vtkm::Id WaveDecompose( const SignalArrayType &sigIn, // Input
vtkm::Id nLevels, // n levels of DWT
CoeffArrayType &coeffOut,
std::vector<vtkm::Id> &L)
std::vector<vtkm::Id> &L,
DeviceTag )
{
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
if( nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel( sigInLen ) )
@ -55,7 +56,8 @@ public:
}
if( nLevels == 0 ) // 0 levels means no transform
{
WaveletBase::DeviceCopy( sigIn, coeffOut );
//WaveletBase::DeviceCopy( sigIn, coeffOut, DeviceTag );
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
return 0;
}
@ -79,7 +81,8 @@ public:
typedef vtkm::cont::ArrayHandlePermutation< IdArrayType, CoeffArrayType >
PermutArrayType;
WaveletBase::DeviceCopy( sigIn, coeffOut );
//WaveletBase::DeviceCopy( sigIn, coeffOut );
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
for( vtkm::Id i = nLevels; i > 0; i-- )
{
@ -108,12 +111,13 @@ public:
// Multi-level 1D wavelet reconstruction
template< typename CoeffArrayType, typename SignalArrayType >
template< typename CoeffArrayType, typename SignalArrayType, typename DeviceTag >
VTKM_CONT_EXPORT
vtkm::Id WaveReconstruct( const CoeffArrayType &coeffIn, // Input
vtkm::Id nLevels, // n levels of DWT
std::vector<vtkm::Id> &L,
SignalArrayType &sigOut )
SignalArrayType &sigOut,
DeviceTag )
{
VTKM_ASSERT( nLevels > 0 );
vtkm::Id LLength = nLevels + 2;
@ -130,7 +134,8 @@ public:
typedef vtkm::cont::ArrayHandlePermutation< IdArrayType, SignalArrayType >
PermutArrayType;
WaveletBase::DeviceCopy( coeffIn, sigOut );
//WaveletBase::DeviceCopy( coeffIn, sigOut );
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( coeffIn, sigOut );
for( vtkm::Id i = 1; i <= nLevels; i++ )
{
@ -288,10 +293,10 @@ public:
// Squash coefficients smaller than a threshold
template< typename CoeffArrayType >
VTKM_CONT_EXPORT
template< typename CoeffArrayType, typename DeviceTag >
vtkm::Id SquashCoefficients( CoeffArrayType &coeffIn,
vtkm::Float64 ratio )
vtkm::Float64 ratio,
DeviceTag )
{
if( ratio > 1 )
{
@ -299,8 +304,9 @@ public:
typedef typename CoeffArrayType::ValueType ValueType;
typedef vtkm::cont::ArrayHandle< ValueType > CoeffArrayBasic;
CoeffArrayBasic sortedArray;
WaveletBase::DeviceCopy( coeffIn, sortedArray );
WaveletBase::DeviceSort( sortedArray );
//WaveletBase::DeviceCopy( coeffIn, sortedArray );
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( coeffIn, sortedArray );
WaveletBase::DeviceSort( sortedArray, DeviceTag() );
vtkm::Id n = coeffLen -
static_cast<vtkm::Id>( static_cast<vtkm::Float64>(coeffLen)/ratio );
@ -315,14 +321,14 @@ public:
// Report statistics on reconstructed array
template< typename ArrayType >
VTKM_CONT_EXPORT
template< typename ArrayType, typename DeviceTag >
vtkm::Id EvaluateReconstruction( const ArrayType &original,
const ArrayType &reconstruct )
const ArrayType &reconstruct,
DeviceTag )
{
#define VAL vtkm::Float64
#define MAKEVAL(a) (static_cast<VAL>(a))
VAL VarOrig = WaveletBase::CalculateVariance( original );
VAL VarOrig = WaveletBase::DeviceCalculateVariance( original, DeviceTag() );
typedef typename ArrayType::ValueType ValueType;
typedef vtkm::cont::ArrayHandle< ValueType > ArrayBasic;
@ -339,7 +345,7 @@ public:
vtkm::worklet::DispatcherMapField< SquareWorklet > swDispatcher( sw );
swDispatcher.Invoke( errorArray, errorSquare );
VAL varErr = WaveletBase::CalculateVariance( errorArray );
VAL varErr = WaveletBase::DeviceCalculateVariance( errorArray, DeviceTag() );
VAL snr, decibels;
if( varErr != 0.0 )
{
@ -352,12 +358,12 @@ public:
decibels = vtkm::Infinity64();
}
VAL origMax = WaveletBase::DeviceMax( original );
VAL origMin = WaveletBase::DeviceMin( original );
VAL errorMax = WaveletBase::DeviceMaxAbs( errorArray );
VAL origMax = WaveletBase::DeviceMax( original, DeviceTag() );
VAL origMin = WaveletBase::DeviceMin( original, DeviceTag() );
VAL errorMax = WaveletBase::DeviceMaxAbs( errorArray, DeviceTag() );
VAL range = origMax - origMin;
VAL squareSum = WaveletBase::DeviceSum( errorSquare );
VAL squareSum = WaveletBase::DeviceSum( errorSquare, DeviceTag() );
VAL rmse = vtkm::Sqrt( MAKEVAL(squareSum) / MAKEVAL(errorArray.GetNumberOfValues()) );
std::cout << "Data range = " << range << std::endl;

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

@ -29,14 +29,10 @@
void DebugDWTIDWT1D()
{
vtkm::Id sigLen = 20;
vtkm::Id sigLen = 21;
std::cout << "Testing Wavelets Worklet" << std::endl;
std::cout << "Input a size to test." << std::endl;
vtkm::Id tmpIn;
vtkm::Id million = 1000000;
std::cin >> tmpIn;
if( tmpIn != 0 )
sigLen = tmpIn * million;
std::cin >> sigLen;
// make input data array handle
std::vector<vtkm::Float64> tmpVector;
@ -49,7 +45,8 @@ void DebugDWTIDWT1D()
std::vector<vtkm::Id> L(3, 0);
// Forward Transform
vtkm::worklet::wavelets::WaveletDWT waveletdwt( "CDF9/7" );
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF8_4;
vtkm::worklet::wavelets::WaveletDWT waveletdwt( wname );
waveletdwt.DWT1D( inputArray, coeffOut, L );
std::cout << "Forward Wavelet Transform: result coeff length = " <<
@ -158,20 +155,18 @@ void DebugWaveDecomposeReconstruct()
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 a new size to test." << std::endl;
std::cout << "Input 0 to stick with 20." << std::endl;
vtkm::Id tmpIn;
vtkm::Id million = 1000000;
vtkm::Id thousand = 1000;
//std::cin >> tmpIn;
tmpIn = 100;
std::cin >> tmpIn;
if( tmpIn != 0 )
sigLen = tmpIn * million;
sigLen = tmpIn;
// make input data array handle
std::vector<vtkm::Float64> tmpVector;
for( vtkm::Id i = 0; i < sigLen; i++ )
tmpVector.push_back( 100.0 * vtkm::Sin(static_cast<vtkm::Float64>(i)/100.0 ));
tmpVector.push_back( static_cast<vtkm::Float64>(i) );
//tmpVector.push_back( 100.0 * vtkm::Sin(static_cast<vtkm::Float64>(i)/100.0 ));
vtkm::cont::ArrayHandle<vtkm::Float64> inputArray =
vtkm::cont::make_ArrayHandle(tmpVector);
@ -179,15 +174,15 @@ void DebugWaveDecomposeReconstruct()
// Use a WaveletCompressor
vtkm::Id nLevels = 2;
vtkm::worklet::WaveletCompressor compressor("CDF9/7");
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF8_4;
vtkm::worklet::WaveletCompressor compressor( wname );
// User input of decompose levels
vtkm::Id maxLevel = compressor.GetWaveletMaxLevel( sigLen );
std::cout << "Input how many wavelet transform levels to perform, between 1 and "
<< maxLevel << std::endl;
vtkm::Id levTemp;
//std::cin >> levTemp;
levTemp = 17;
std::cin >> levTemp;
if( levTemp > 0 && levTemp <= maxLevel )
nLevels = levTemp;
else
@ -198,35 +193,34 @@ void DebugWaveDecomposeReconstruct()
std::cout << "Input a compression ratio ( >=1 )to test. "
<< "1 means no compression. " << std::endl;
vtkm::Float64 cratio;
//std::cin >> cratio;
cratio = 10;
std::cin >> cratio;
VTKM_ASSERT ( cratio >= 1 );
std::vector<vtkm::Id> L;
// Decompose
vtkm::cont::Timer<> timer;
compressor.WaveDecompose( inputArray, nLevels, outputArray, L );
compressor.WaveDecompose( inputArray, nLevels, outputArray, L, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
vtkm::Float64 elapsedTime = timer.GetElapsedTime();
std::cout << "Decompose time = " << elapsedTime << std::endl;
// Squash small coefficients
timer.Reset();
compressor.SquashCoefficients( outputArray, cratio );
compressor.SquashCoefficients( outputArray, cratio, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
elapsedTime = timer.GetElapsedTime();
std::cout << "Thresholding time = " << elapsedTime << std::endl;
// Reconstruct
vtkm::cont::ArrayHandle<vtkm::Float64> reconstructArray;
timer.Reset();
compressor.WaveReconstruct( outputArray, nLevels, L, reconstructArray );
compressor.WaveReconstruct( outputArray, nLevels, L, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
elapsedTime = timer.GetElapsedTime();
std::cout << "Reconstruction time = " << elapsedTime << std::endl;
compressor.EvaluateReconstruction( inputArray, reconstructArray );
compressor.EvaluateReconstruction( inputArray, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
timer.Reset();
for( vtkm::Id i = 0; i < reconstructArray.GetNumberOfValues(); i++ )
@ -251,14 +245,17 @@ void TestWaveDecomposeReconstruct()
// make input data array handle
std::vector<vtkm::Float64> tmpVector;
for( vtkm::Id i = 0; i < sigLen; i++ )
{
tmpVector.push_back( 100.0 * vtkm::Sin(static_cast<vtkm::Float64>(i)/100.0 ));
}
vtkm::cont::ArrayHandle<vtkm::Float64> inputArray =
vtkm::cont::make_ArrayHandle(tmpVector);
vtkm::cont::ArrayHandle<vtkm::Float64> outputArray;
// Use a WaveletCompressor
vtkm::worklet::WaveletCompressor compressor("CDF9/7");
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF8_4;
vtkm::worklet::WaveletCompressor compressor( wname );
// User maximum decompose levels, and no compression
vtkm::Id maxLevel = compressor.GetWaveletMaxLevel( sigLen );
@ -268,7 +265,7 @@ void TestWaveDecomposeReconstruct()
// Decompose
vtkm::cont::Timer<> timer;
compressor.WaveDecompose( inputArray, nLevels, outputArray, L );
compressor.WaveDecompose( inputArray, nLevels, outputArray, L, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
vtkm::Float64 elapsedTime = timer.GetElapsedTime();
std::cout << "Decompose time = " << elapsedTime << std::endl;
@ -276,11 +273,11 @@ void TestWaveDecomposeReconstruct()
// Reconstruct
vtkm::cont::ArrayHandle<vtkm::Float64> reconstructArray;
timer.Reset();
compressor.WaveReconstruct( outputArray, nLevels, L, reconstructArray );
compressor.WaveReconstruct( outputArray, nLevels, L, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
elapsedTime = timer.GetElapsedTime();
std::cout << "Reconstruction time = " << elapsedTime << std::endl;
compressor.EvaluateReconstruction( inputArray, reconstructArray );
compressor.EvaluateReconstruction( inputArray, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
timer.Reset();
for( vtkm::Id i = 0; i < reconstructArray.GetNumberOfValues(); i++ )

66
vtkm/worklet/wavelets/FilterBanks.h
View File

@ -32,6 +32,7 @@ namespace worklet {
namespace wavelets {
const vtkm::Float64 hm4_44[9] = {
/* From VAPoR
0.037828455507264,
-0.023849465019557,
-0.110624404418437,
@ -40,10 +41,23 @@ namespace wavelets {
0.377402855612831,
-0.110624404418437,
-0.023849465019557,
0.037828455507264
0.037828455507264 */
/* From http://wavelets.pybytes.com/wavelet/bior4.4/ and its git repo:
* https://github.com/nigma/pywt/blob/035e1fa14c2cd70ca270da20b1523e834a7ae635/src/wavelets_coeffs.template.h */
0.03782845550726404,
-0.023849465019556843,
-0.11062440441843718,
0.37740285561283066,
0.85269867900889385,
0.37740285561283066,
-0.11062440441843718,
-0.023849465019556843,
0.03782845550726404
};
const vtkm::Float64 h4[9] = {
/* From VAPoR
0.0,
-0.064538882628697,
-0.040689417609164,
@ -52,10 +66,22 @@ namespace wavelets {
0.418092273221617,
-0.0406894176091641,
-0.0645388826286971,
0.0
0.0 */
/* From http://wavelets.pybytes.com/wavelet/bior4.4/ and its git repo:
* https://github.com/nigma/pywt/blob/035e1fa14c2cd70ca270da20b1523e834a7ae635/src/wavelets_coeffs.template.h */
0.0,
-0.064538882628697058,
-0.040689417609164058,
0.41809227322161724,
0.7884856164055829,
0.41809227322161724,
-0.040689417609164058,
-0.064538882628697058,
0.0
};
const double hm2_22[6] = {
const vtkm::Float64 hm2_22[6] = {
-0.1767766952966368811002110905262,
0.3535533905932737622004221810524,
1.0606601717798212866012665431573,
@ -63,7 +89,7 @@ namespace wavelets {
-0.1767766952966368811002110905262
};
const double h2[18] = {
const vtkm::Float64 h2[18] = {
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.3535533905932737622004221810524,
0.7071067811865475244008443621048,
@ -72,6 +98,38 @@ namespace wavelets {
0.0
};
const vtkm::Float64 hm1_11[2] = {
0.70710678118654752440084436210,
0.70710678118654752440084436210
};
const vtkm::Float64 h1[10] = {
0.0, 0.0, 0.0, 0.0,
0.70710678118654752440084436210,
0.70710678118654752440084436210,
0.0, 0.0, 0.0, 0.0
};
const vtkm::Float64 hm3_33[8] = {
0.0662912607362388304125791589473,
-0.1988737822087164912377374768420,
-0.1546796083845572709626847042104,
0.9943689110435824561886873842099,
0.9943689110435824561886873842099,
-0.1546796083845572709626847042104,
-0.1988737822087164912377374768420,
0.0662912607362388304125791589473
};
const vtkm::Float64 h3[20] = {
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.1767766952966368811002110905262,
0.5303300858899106433006332715786,
0.5303300858899106433006332715786,
0.1767766952966368811002110905262,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
};
};
}

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

@ -34,16 +34,10 @@ namespace worklet {
namespace wavelets {
enum DWTMode { // boundary extension modes
INVALID = -1,
ZPD,
SYMH,
SYMW,
ASYMH,
ASYMW,
SP0,
SP1,
PPD,
PER
ASYMW
};
// Functionalities are similar to MatWaveBase in VAPoR.
@ -52,53 +46,27 @@ class WaveletBase
public:
// Constructor
WaveletBase( const std::string &w_name )
WaveletBase( WaveletName name ) : wname ( name ),
filter( name )
{
this->filter = NULL;
this->wmode = PER;
this->wname = w_name;
if( wname.compare("CDF9/7") == 0 )
if( wname == CDF9_7 || wname == BIOR4_4 ||
wname == CDF5_3 || wname == BIOR2_2 )
{
this->wmode = SYMW; // Default extension mode, see MatWaveBase.cpp
this->filter = new vtkm::worklet::wavelets::WaveletFilter( wname );
}
else if( wname.compare("CDF5/3") == 0 )
else if( wname == HAAR || wname == BIOR1_1 ||
wname == CDF8_4 || wname == BIOR3_3 )
{
this->wmode = SYMW;
this->filter = new vtkm::worklet::wavelets::WaveletFilter( wname );
this->wmode = SYMH;
}
else
{
vtkm::cont::ErrorControlBadValue("This wavelet kernel is not supported! ");
}
}
// Destructor
virtual ~WaveletBase()
{
if( filter ) delete filter;
filter = NULL;
}
// Get the wavelet filter
const WaveletFilter* GetWaveletFilter()
{
if( this->filter == NULL )
{
vtkm::cont::ErrorControlInternal( "Filter is NULL right now!" );
}
return filter;
}
// Returns length of approximation coefficients from a decompostition pass.
vtkm::Id GetApproxLength( vtkm::Id sigInLen )
{
vtkm::Id filterLen = this->filter->GetFilterLength();
vtkm::Id filterLen = this->filter.GetFilterLength();
if (this->wmode == PER)
return static_cast<vtkm::Id>(vtkm::Ceil( (static_cast<vtkm::Float64>(sigInLen)) / 2.0 ));
else if (this->filter->isSymmetric())
if (this->filter.isSymmetric())
{
if ( (this->wmode == SYMW && (filterLen % 2 != 0)) ||
(this->wmode == SYMH && (filterLen % 2 == 0)) )
@ -117,11 +85,9 @@ public:
// Returns length of detail coefficients from a decompostition pass
vtkm::Id GetDetailLength( vtkm::Id sigInLen )
{
vtkm::Id filterLen = this->filter->GetFilterLength();
vtkm::Id filterLen = this->filter.GetFilterLength();
if (this->wmode == PER)
return static_cast<vtkm::Id>(vtkm::Ceil( (static_cast<vtkm::Float64>(sigInLen)) / 2.0 ));
else if (this->filter->isSymmetric())
if (this->filter.isSymmetric())
{
if ( (this->wmode == SYMW && (filterLen % 2 != 0)) ||
(this->wmode == SYMH && (filterLen % 2 == 0)) )
@ -154,29 +120,15 @@ public:
// 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;
}
}
// perform a device copy
template< typename ArrayType1, typename ArrayType2 >
VTKM_EXEC_CONT_EXPORT
void DeviceCopy( const ArrayType1 &srcArray,
ArrayType2 &dstArray)
{
vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Copy
( srcArray, dstArray );
vtkm::Id filterLen = this->filter.GetFilterLength();
vtkm::Id level;
this->WaveLengthValidate( sigInLen, filterLen, level );
return level;
}
// perform a device copy. The whole 1st array to a certain start location of the 2nd array
template< typename ArrayType1, typename ArrayType2 >
VTKM_EXEC_CONT_EXPORT
VTKM_EXEC_EXPORT
void DeviceCopyStartX( const ArrayType1 &srcArray,
ArrayType2 &dstArray,
vtkm::Id startIdx)
@ -187,30 +139,41 @@ public:
dispatcher.Invoke( srcArray, dstArray );
}
// Assign zero value to a certain location of an array
template< typename ArrayType >
VTKM_EXEC_EXPORT
void DeviceAssignZero( ArrayType &array, vtkm::Id index )
{
typedef vtkm::worklet::wavelets::AssignZeroWorklet ZeroWorklet;
ZeroWorklet worklet( index );
vtkm::worklet::DispatcherMapField< ZeroWorklet > dispatcher( worklet );
dispatcher.Invoke( array );
}
// Sort by the absolute value on device
struct SortLessAbsFunctor
{
template< typename T >
VTKM_EXEC_CONT_EXPORT
VTKM_EXEC_EXPORT
bool operator()(const T& x, const T& y) const
{
return vtkm::Abs(x) < vtkm::Abs(y);
}
};
template< typename ArrayType >
VTKM_EXEC_CONT_EXPORT
void DeviceSort( ArrayType &array )
template< typename ArrayType, typename DeviceTag >
VTKM_EXEC_EXPORT
void DeviceSort( ArrayType &array, DeviceTag )
{
vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Sort
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Sort
( array, SortLessAbsFunctor() );
}
// Reduce to the sum of all values on device
template< typename ArrayType >
VTKM_EXEC_CONT_EXPORT
typename ArrayType::ValueType DeviceSum( const ArrayType &array )
template< typename ArrayType, typename DeviceTag >
VTKM_EXEC_EXPORT
typename ArrayType::ValueType DeviceSum( const ArrayType &array, DeviceTag )
{
return vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Reduce
return vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Reduce
( array, 0.0 );
}
@ -218,7 +181,7 @@ public:
struct minFunctor
{
template< typename FieldType >
VTKM_EXEC_CONT_EXPORT
VTKM_EXEC_EXPORT
FieldType operator()(const FieldType &x, const FieldType &y) const {
return Min(x, y);
}
@ -226,25 +189,25 @@ public:
struct maxFunctor
{
template< typename FieldType >
VTKM_EXEC_CONT_EXPORT
VTKM_EXEC_EXPORT
FieldType operator()(const FieldType& x, const FieldType& y) const {
return Max(x, y);
}
};
template< typename ArrayType >
VTKM_EXEC_CONT_EXPORT
typename ArrayType::ValueType DeviceMax( const ArrayType &array )
template< typename ArrayType, typename DeviceTag >
VTKM_EXEC_EXPORT
typename ArrayType::ValueType DeviceMax( const ArrayType &array, DeviceTag )
{
typename ArrayType::ValueType initVal = array.GetPortalConstControl().Get(0);
return vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Reduce
return vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Reduce
( array, initVal, maxFunctor() );
}
template< typename ArrayType >
VTKM_EXEC_CONT_EXPORT
typename ArrayType::ValueType DeviceMin( const ArrayType &array )
template< typename ArrayType, typename DeviceTag >
VTKM_EXEC_EXPORT
typename ArrayType::ValueType DeviceMin( const ArrayType &array, DeviceTag )
{
typename ArrayType::ValueType initVal = array.GetPortalConstControl().Get(0);
return vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Reduce
return vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Reduce
( array, initVal, minFunctor() );
}
@ -252,26 +215,25 @@ public:
struct maxAbsFunctor
{
template< typename FieldType >
VTKM_EXEC_CONT_EXPORT
VTKM_EXEC_EXPORT
FieldType operator()(const FieldType& x, const FieldType& y) const {
return Max( vtkm::Abs(x), vtkm::Abs(y) );
}
};
template< typename ArrayType >
VTKM_EXEC_CONT_EXPORT
typename ArrayType::ValueType DeviceMaxAbs( const ArrayType &array )
template< typename ArrayType, typename DeviceTag >
VTKM_EXEC_EXPORT
typename ArrayType::ValueType DeviceMaxAbs( const ArrayType &array, DeviceTag )
{
typename ArrayType::ValueType initVal = array.GetPortalConstControl().Get(0);
return vtkm::cont::DeviceAdapterAlgorithm< VTKM_DEFAULT_DEVICE_ADAPTER_TAG>::Reduce
return vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Reduce
( array, initVal, maxAbsFunctor() );
}
// Calculate variance of an array
template< typename ArrayType >
VTKM_CONT_EXPORT
vtkm::Float64 CalculateVariance( ArrayType &array )
template< typename ArrayType, typename DeviceTag >
vtkm::Float64 DeviceCalculateVariance( ArrayType &array, DeviceTag )
{
vtkm::Float64 mean = static_cast<vtkm::Float64>(WaveletBase::DeviceSum( array )) /
vtkm::Float64 mean = static_cast<vtkm::Float64>(this->DeviceSum( array, DeviceTag() )) /
static_cast<vtkm::Float64>(array.GetNumberOfValues());
vtkm::cont::ArrayHandle< vtkm::Float64 > squaredDeviation;
@ -282,7 +244,7 @@ public:
vtkm::worklet::DispatcherMapField< SDWorklet > dispatcher( sdw );
dispatcher.Invoke( array, squaredDeviation );
vtkm::Float64 sdMean = this->DeviceSum( squaredDeviation ) /
vtkm::Float64 sdMean = this->DeviceSum( squaredDeviation, DeviceTag() ) /
static_cast<vtkm::Float64>( squaredDeviation.GetNumberOfValues() );
return sdMean;
@ -340,9 +302,9 @@ public:
protected:
vtkm::worklet::wavelets::DWTMode wmode;
WaveletFilter* filter;
std::string wname;
WaveletName wname;
DWTMode wmode;
WaveletFilter filter;
void WaveLengthValidate( vtkm::Id sigInLen, vtkm::Id filterLength, vtkm::Id &level)
{

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

@ -31,6 +31,7 @@
#include <vtkm/cont/ArrayHandlePermutation.h>
#include <vtkm/Math.h>
#include <vtkm/cont/Timer.h>
namespace vtkm {
namespace worklet {
@ -41,7 +42,7 @@ class WaveletDWT : public WaveletBase
public:
// Constructor
WaveletDWT( const std::string &w_name ) : WaveletBase( w_name ) {}
WaveletDWT( WaveletName name ) : WaveletBase( name ) {}
// Func: Extend 1D signal
@ -52,10 +53,15 @@ public:
vtkm::Id addLen,
vtkm::worklet::wavelets::DWTMode leftExtMethod,
vtkm::worklet::wavelets::DWTMode rightExtMethod,
bool attachRightZero )
bool attachZeroRightLeft,
bool attachZeroRightRight )
{
// "right extension" can be attached a zero on either end, but not both ends.
VTKM_ASSERT( !attachZeroRightRight || !attachZeroRightLeft );
typedef typename SigInArrayType::ValueType ValueType;
typedef vtkm::cont::ArrayHandle< ValueType > ExtensionArrayType;
ExtensionArrayType leftExtend;
leftExtend.Allocate( addLen );
@ -65,23 +71,41 @@ public:
typedef vtkm::worklet::wavelets::LeftSYMWExtentionWorklet LeftSYMW;
typedef vtkm::worklet::wavelets::RightSYMHExtentionWorklet RightSYMH;
typedef vtkm::worklet::wavelets::RightSYMWExtentionWorklet RightSYMW;
typedef vtkm::worklet::wavelets::LeftASYMHExtentionWorklet LeftASYMH;
typedef vtkm::worklet::wavelets::LeftASYMWExtentionWorklet LeftASYMW;
typedef vtkm::worklet::wavelets::RightASYMHExtentionWorklet RightASYMH;
typedef vtkm::worklet::wavelets::RightASYMWExtentionWorklet RightASYMW;
switch( leftExtMethod )
{
case vtkm::worklet::wavelets::SYMH:
case SYMH:
{
LeftSYMH worklet( addLen );
vtkm::worklet::DispatcherMapField< LeftSYMH > dispatcher( worklet );
dispatcher.Invoke( leftExtend, sigIn );
break;
}
case vtkm::worklet::wavelets::SYMW:
case SYMW:
{
LeftSYMW worklet( addLen );
vtkm::worklet::DispatcherMapField< LeftSYMW > dispatcher( worklet );
dispatcher.Invoke( leftExtend, sigIn );
break;
}
case ASYMH:
{
LeftASYMH worklet( addLen );
vtkm::worklet::DispatcherMapField< LeftASYMH > dispatcher( worklet );
dispatcher.Invoke( leftExtend, sigIn );
break;
}
case ASYMW:
{
LeftASYMW worklet( addLen );
vtkm::worklet::DispatcherMapField< LeftASYMW > dispatcher( worklet );
dispatcher.Invoke( leftExtend, sigIn );
break;
}
default:
{
vtkm::cont::ErrorControlInternal("Left extension mode not supported!");
@ -90,40 +114,110 @@ public:
}
ExtensionArrayType rightExtend;
if( attachRightZero )
rightExtend.Allocate( addLen + 1 );
else
if( !attachZeroRightLeft ) // no attach zero, or only attach on RightRight
{
// Allocate memory
if( attachZeroRightRight )
rightExtend.Allocate( addLen + 1 );
else
rightExtend.Allocate( addLen );
switch( rightExtMethod )
{
case SYMH:
{
RightSYMH worklet( sigInLen );
vtkm::worklet::DispatcherMapField< RightSYMH > dispatcher( worklet );
dispatcher.Invoke( rightExtend, sigIn );
break;
}
case SYMW:
{
RightSYMW worklet( sigInLen );
vtkm::worklet::DispatcherMapField< RightSYMW > dispatcher( worklet );
dispatcher.Invoke( rightExtend, sigIn );
break;
}
case ASYMH:
{
RightASYMH worklet( sigInLen );
vtkm::worklet::DispatcherMapField< RightASYMH > dispatcher( worklet );
dispatcher.Invoke( rightExtend, sigIn );
break;
}
case ASYMW:
{
RightASYMW worklet( sigInLen );
vtkm::worklet::DispatcherMapField< RightASYMW > dispatcher( worklet );
dispatcher.Invoke( rightExtend, sigIn );
break;
}
default:
{
vtkm::cont::ErrorControlInternal("Right extension mode not supported!");
return 1;
}
}
if( attachZeroRightRight )
WaveletBase::DeviceAssignZero( rightExtend, addLen );
}
else // attachZeroRightLeft mode
{
typedef vtkm::cont::ArrayHandleConcatenate<SigInArrayType, ExtensionArrayType>
ConcatArray;
// attach a zero at the end of sigIn
ExtensionArrayType singleValArray;
singleValArray.Allocate(1);
WaveletBase::DeviceAssignZero( singleValArray, 0 );
ConcatArray sigInPlusOne( sigIn, singleValArray );
// allocate memory for extension
rightExtend.Allocate( addLen );
switch( rightExtMethod )
{
case vtkm::worklet::wavelets::SYMH:
switch( rightExtMethod )
{
RightSYMH worklet( sigInLen );
vtkm::worklet::DispatcherMapField< RightSYMH > dispatcher( worklet );
dispatcher.Invoke( rightExtend, sigIn );
break;
case SYMH:
{
RightSYMH worklet( sigInLen + 1 );
vtkm::worklet::DispatcherMapField< RightSYMH > dispatcher( worklet );
dispatcher.Invoke( rightExtend, sigInPlusOne );
break;
}
case SYMW:
{
RightSYMW worklet( sigInLen + 1 );
vtkm::worklet::DispatcherMapField< RightSYMW > dispatcher( worklet );
dispatcher.Invoke( rightExtend, sigInPlusOne );
break;
}
case ASYMH:
{
RightASYMH worklet( sigInLen + 1 );
vtkm::worklet::DispatcherMapField< RightASYMH > dispatcher( worklet );
dispatcher.Invoke( rightExtend, sigInPlusOne );
break;
}
case ASYMW:
{
RightASYMW worklet( sigInLen + 1 );
vtkm::worklet::DispatcherMapField< RightASYMW > dispatcher( worklet );
dispatcher.Invoke( rightExtend, sigInPlusOne );
break;
}
default:
{
vtkm::cont::ErrorControlInternal("Right extension mode not supported!");
return 1;
}
}
case SYMW:
{
RightSYMW worklet( sigInLen );
vtkm::worklet::DispatcherMapField< RightSYMW > dispatcher( worklet );
dispatcher.Invoke( rightExtend, sigIn );
break;
}
default:
{
vtkm::cont::ErrorControlInternal("Right extension mode not supported!");
return 1;
}
}
if( attachRightZero )
{
typedef vtkm::worklet::wavelets::AssignZeroWorklet ZeroWorklet;
ZeroWorklet worklet( addLen );
vtkm::worklet::DispatcherMapField< ZeroWorklet > dispatcher( worklet );
dispatcher.Invoke( rightExtend );
// make a copy of rightExtend with a zero attached to the left
ExtensionArrayType rightExtendPlusOne;
rightExtendPlusOne.Allocate( addLen + 1 );
WaveletBase::DeviceCopyStartX( rightExtend, rightExtendPlusOne, 1 );
WaveletBase::DeviceAssignZero( rightExtendPlusOne, 0 );
rightExtend = rightExtendPlusOne ;
}
typedef vtkm::cont::ArrayHandleConcatenate< ExtensionArrayType, SigInArrayType>
@ -158,10 +252,10 @@ public:
VTKM_ASSERT( L[0] + L[1] == L[2] );
vtkm::Id filterLen = WaveletBase::filter->GetFilterLength();
vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
bool doSymConv = false;
if( WaveletBase::filter->isSymmetric() )
if( WaveletBase::filter.isSymmetric() )
{
if( ( WaveletBase::wmode == SYMW && ( filterLen % 2 != 0 ) ) ||
( WaveletBase::wmode == SYMH && ( filterLen % 2 == 0 ) ) )
@ -197,10 +291,10 @@ public:
ConcatType2 sigInExtended;
this->Extend1D( sigIn, sigInExtended, addLen,
WaveletBase::wmode, WaveletBase::wmode, false );
WaveletBase::wmode, WaveletBase::wmode, false, false );
VTKM_ASSERT( sigInExtended.GetNumberOfValues() == sigExtendedLen );
// initialize a worklet
// initialize a worklet for forward transform
vtkm::worklet::wavelets::ForwardTransform forwardTransform;
forwardTransform.SetFilterLength( filterLen );
forwardTransform.SetCoeffLength( L[0], L[1] );
@ -210,8 +304,8 @@ public:
vtkm::worklet::DispatcherMapField<vtkm::worklet::wavelets::ForwardTransform>
dispatcher(forwardTransform);
dispatcher.Invoke( sigInExtended,
WaveletBase::filter->GetLowDecomposeFilter(),
WaveletBase::filter->GetHighDecomposeFilter(),
WaveletBase::filter.GetLowDecomposeFilter(),
WaveletBase::filter.GetHighDecomposeFilter(),
coeffOut );
VTKM_ASSERT( L[0] + L[1] <= coeffOut.GetNumberOfValues() );
@ -231,16 +325,16 @@ public:
SignalArrayType &sigOut ) // Output
{
VTKM_ASSERT( L.size() == 3 );
VTKM_ASSERT( coeffIn.GetNumberOfValues() == L[2] );
VTKM_ASSERT( L[2] == coeffIn.GetNumberOfValues() );
vtkm::Id filterLen = WaveletBase::filter->GetFilterLength();
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::filter.isSymmetric() )
{
if(( WaveletBase::wmode == SYMW && (filterLen % 2 != 0) ) ||
( WaveletBase::wmode == SYMH && (filterLen % 2 == 0) ) )
@ -272,25 +366,23 @@ public:
}
}
vtkm::Id cATempLen, cDTempLen, reconTempLen;
vtkm::Id cATempLen, cDTempLen; //, reconTempLen;
vtkm::Id addLen = 0;
vtkm::Id cDPadLen = 0;
if( doSymConv ) // extend cA and cD
{
addLen = filterLen / 4;
addLen = filterLen / 4; // addLen == 0 for Haar kernel
if( (L[0] > L[1]) && (WaveletBase::wmode == SYMH) )
cDPadLen = L[0]; // SYMH is rare
cDPadLen = L[0];
cATempLen = L[0] + 2 * addLen;
cDTempLen = cATempLen; // same length
}
else // not extend cA and cD
{
{ // (biorthogonal kernels won't come into this case)
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 >
@ -314,27 +406,26 @@ public:
if( doSymConv ) // Actually extend cA and cD
{
this->Extend1D( cA, cATemp, addLen, cALeftMode, cARightMode, false );
// first extend cA to be cATemp
this->Extend1D( cA, cATemp, addLen, cALeftMode, cARightMode, false, false );
// Then extend cD based on extension needs
if( cDPadLen > 0 )
{
/* Add back the missing final cD: 0.0
* TODO when SYMH is needed.
*
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 );
*/
// Add back the missing final cD, 0.0, before doing extension
this->Extend1D( cD, cDTemp, addLen, cDLeftMode, cDRightMode, true, false );
}
else
{
vtkm::Id cDTempLenWouldBe = cD.GetNumberOfValues() + 2 * addLen;
vtkm::Id cDTempLenWouldBe = L[1] + 2 * addLen;
if( cDTempLenWouldBe == cDTempLen )
this->Extend1D( cD, cDTemp, addLen, cDLeftMode, cDRightMode, false);
{
this->Extend1D( cD, cDTemp, addLen, cDLeftMode, cDRightMode, false, false);
}
else if( cDTempLenWouldBe == cDTempLen - 1 )
this->Extend1D( cD, cDTemp, addLen, cDLeftMode, cDRightMode, true);
{
this->Extend1D( cD, cDTemp, addLen, cDLeftMode, cDRightMode, false, true );
}
else
{
vtkm::cont::ErrorControlInternal("cDTemp Length not match!");
@ -342,40 +433,56 @@ public:
}
}
}
else
{ /* TODO when SYMH is needed
WaveletBase::DeviceCopy( cA, cATemp );
WaveletBase::DeviceCopy( cD, cDTemp );
*/
else // !doSymConv (biorthogonals kernel won't come into this case)
{
// make cATemp
ExtensionArrayType dummyArray;
dummyArray.Allocate(0);
Concat1 cALeftOn( dummyArray, cA );
cATemp = vtkm::cont::make_ArrayHandleConcatenate< Concat1, ExtensionArrayType >
( cALeftOn, dummyArray );
// make cDTemp
Concat1 cDLeftOn( dummyArray, cD );
cDTemp = vtkm::cont::make_ArrayHandleConcatenate< Concat1, ExtensionArrayType >
( cDLeftOn, dummyArray );
}
// make sure signal extension went as expected
VTKM_ASSERT( cATemp.GetNumberOfValues() == cATempLen );
VTKM_ASSERT( cDTemp.GetNumberOfValues() == cDTempLen );
vtkm::cont::ArrayHandleConcatenate< Concat2, Concat2>
coeffInExtended( cATemp, cDTemp );
// Allocate memory for sigOut
sigOut.Allocate( cATempLen + cDTempLen );
if( filterLen % 2 != 0 )
{
vtkm::cont::ArrayHandleConcatenate< Concat2, Concat2>
coeffInExtended( cATemp, cDTemp );
// Allocate memory for sigOut
sigOut.Allocate( coeffInExtended.GetNumberOfValues() );
// 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(),
WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
sigOut );
VTKM_ASSERT( sigOut.GetNumberOfValues() >= L[2] );
sigOut.Shrink( L[2] );
}
else
{
// TODO: implement for even length filter
vtkm::worklet::wavelets::InverseTransformEven inverseXformEven
( filterLen, L[0], cATempLen, !doSymConv );
vtkm::worklet::DispatcherMapField<vtkm::worklet::wavelets::InverseTransformEven>
dispatcher( inverseXformEven );
dispatcher.Invoke( coeffInExtended,
WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
sigOut );
}
sigOut.Shrink( L[2] );
return 0;
} // function IDWT1D

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

@ -32,54 +32,81 @@ namespace worklet {
namespace wavelets {
enum WaveletName {
CDF9_7,
CDF5_3,
CDF8_4,
HAAR,
BIOR4_4, // the same as CDF9_7
BIOR3_3, // the same as CDF8_4
BIOR2_2, // the same as CDF5_3
BIOR1_1 // the same as HAAE
};
// Wavelet filter class;
// functionally equivalent to WaveFiltBase and its subclasses in VAPoR.
class WaveletFilter
{
public:
// constructor
WaveletFilter( const std::string &wname )
WaveletFilter( WaveletName wtype ) : symmetricity(true),
filterLength(0),
lowDecomposeFilter(NULL),
highDecomposeFilter(NULL),
lowReconstructFilter(NULL),
highReconstructFilter(NULL)
{
lowDecomposeFilter = highDecomposeFilter =
lowReconstructFilter = highReconstructFilter = NULL;
this->filterLength = 0;
if( wname.compare("CDF9/7") == 0 )
if( wtype == CDF9_7 || wtype == BIOR4_4 )
{
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 );
this->wrev( vtkm::worklet::wavelets::hm4_44, lowDecomposeFilter, filterLength );
this->qmf_wrev( vtkm::worklet::wavelets::h4, highDecomposeFilter, filterLength );
this->verbatim_copy( vtkm::worklet::wavelets::h4, lowReconstructFilter, filterLength );
this->qmf_even( vtkm::worklet::wavelets::hm4_44, highReconstructFilter, filterLength );
}
else if( wname.compare("CDF5/3") == 0 )
else if( wtype == CDF8_4 || wtype == BIOR3_3 )
{
this->filterLength = 8;
this->AllocateFilterMemory();
this->wrev( vtkm::worklet::wavelets::hm3_33, lowDecomposeFilter, filterLength );
this->qmf_wrev( vtkm::worklet::wavelets::h3+6, highDecomposeFilter, filterLength );
this->verbatim_copy( vtkm::worklet::wavelets::h3+6, lowReconstructFilter, filterLength );
this->qmf_even( vtkm::worklet::wavelets::hm3_33, highReconstructFilter, filterLength );
}
else if( wtype == CDF5_3 || wtype == BIOR2_2 )
{
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 );
this->wrev( vtkm::worklet::wavelets::hm2_22, lowDecomposeFilter, filterLength );
this->qmf_wrev( vtkm::worklet::wavelets::h2+6, highDecomposeFilter, filterLength );
this->verbatim_copy( vtkm::worklet::wavelets::h2+6, lowReconstructFilter, filterLength );
this->qmf_even( vtkm::worklet::wavelets::hm2_22, highReconstructFilter, filterLength );
}
else
else if( wtype == HAAR || wtype == BIOR1_1 )
{
vtkm::cont::ErrorControlBadValue("Not supported wavelet kernel type!");
this->filterLength = 2;
this->AllocateFilterMemory();
this->wrev( vtkm::worklet::wavelets::hm1_11, lowDecomposeFilter, filterLength );
this->qmf_wrev( vtkm::worklet::wavelets::h1+4, highDecomposeFilter, filterLength );
this->verbatim_copy( vtkm::worklet::wavelets::h1+4, lowReconstructFilter, filterLength );
this->qmf_even( vtkm::worklet::wavelets::hm1_11, highReconstructFilter, filterLength );
}
}
// destructor
virtual ~WaveletFilter()
{
if( lowDecomposeFilter ) delete[] lowDecomposeFilter;
if( highDecomposeFilter ) delete[] highDecomposeFilter;
if( lowReconstructFilter ) delete[] lowReconstructFilter;
if( highReconstructFilter ) delete[] highReconstructFilter;
if( lowDecomposeFilter )
{
delete[] lowDecomposeFilter;
lowDecomposeFilter = highDecomposeFilter =
lowReconstructFilter = highReconstructFilter = NULL ;
}
}
vtkm::Id GetFilterLength() { return this->filterLength; }
bool isSymmetric() { return this->symmetricity; }
vtkm::Id GetFilterLength() { return filterLength; }
bool isSymmetric() { return symmetricity; }
typedef vtkm::cont::ArrayHandle<vtkm::Float64> FilterType;
FilterType GetLowDecomposeFilter() const
@ -99,7 +126,7 @@ public:
return vtkm::cont::make_ArrayHandle( highReconstructFilter, filterLength );
}
protected:
private:
bool symmetricity;
vtkm::Id filterLength;
vtkm::Float64* lowDecomposeFilter;
@ -109,12 +136,12 @@ protected:
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 ];
lowDecomposeFilter = new vtkm::Float64[ filterLength * 4 ];
highDecomposeFilter = lowDecomposeFilter + filterLength;
lowReconstructFilter = highDecomposeFilter + filterLength;
highReconstructFilter = lowReconstructFilter + filterLength;
}
// Flipping operation; helper function to initialize a filter.
void wrev( const vtkm::Float64* sigIn, vtkm::Float64* sigOut, vtkm::Id sigLength )
{
@ -125,36 +152,30 @@ protected:
// 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++)
if( sigLength % 2 == 0 )
{
sigOut[count] = sigIn[sigLength - count - 1];
if (sigLength % 2 == 0) {
for (vtkm::Id count = 0; count < sigLength; count++)
{
sigOut[count] = sigIn[sigLength - count - 1];
if (count % 2 != 0)
sigOut[count] = -1.0 * sigOut[count];
}
else {
}
else
{
for (vtkm::Id count = 0; count < sigLength; count++)
{
sigOut[count] = sigIn[sigLength - count - 1];
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];
}
}
qmf_even( sigIn, sigOut, sigLength );
vtkm::Float64 tmp;
for (vtkm::Id count = 0; count < sigLength/2; count++) {
@ -170,8 +191,8 @@ protected:
for (vtkm::Id count = 0; count < sigLength; count++)
sigOut[count] = sigIn[count];
}
}; // class WaveletFilter.
}; // class WaveletFilter.
} // namespace wavelets.

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

@ -147,18 +147,12 @@ public:
// Constructor
VTKM_EXEC_CONT_EXPORT
InverseTransformOdd()
{
magicNum = 0.0;
filterLen = 0;
cALen = 0;
}
InverseTransformOdd() : filterLen(0), cALen(0), cALen2(0), cALenExtended(0) {}
// Set the filter length
VTKM_EXEC_CONT_EXPORT
void SetFilterLength( vtkm::Id len )
{
//VTKM_ASSERT( len % 2 == 1 );
this->filterLen = len;
}
@ -166,7 +160,8 @@ public:
VTKM_EXEC_CONT_EXPORT
void SetCALength( vtkm::Id len, vtkm::Id lenExt )
{
this->cALen = len;
this->cALen = len;
this->cALen2 = len * 2;
this->cALenExtended = lenExt;
}
@ -184,39 +179,38 @@ public:
OutputPortalType &sigOut,
const vtkm::Id &workIndex) const
{
vtkm::Id xi; // coeff indices
vtkm::Id k; // filter indices
if( workIndex < 2*cALen ) // valid calculation region
if( workIndex < cALen2 ) // valid calculation region
{
vtkm::Id xi; // coeff indices
vtkm::Id k1, k2; // indices for low and high filter
VAL sum = 0.0;
xi = (workIndex+1) / 2;
if( workIndex % 2 != 0 )
k = this->filterLen - 2;
else
k = this->filterLen - 1;
while( k >= 0 )
{
sum += lowFilter.Get(k) * MAKEVAL( coeffs.Get(xi) );
xi++;
k -= 2;
k1 = this->filterLen - 2;
k2 = this->filterLen - 1;
}
else
{
k1 = this->filterLen - 1;
k2 = this->filterLen - 2;
}
xi = (workIndex+1) / 2;
while( k1 > -1 ) // k1 >= 0
{
sum += lowFilter.Get(k1) * MAKEVAL( coeffs.Get(xi++) );
k1 -= 2;
}
xi = workIndex / 2;
if( workIndex % 2 != 0 )
k = this->filterLen - 1;
else
k = this->filterLen - 2;
while( k >= 0 )
while( k2 > -1 ) // k2 >= 0
{
sum += highFilter.Get(k) * MAKEVAL( coeffs.Get( xi + this->cALenExtended ) );
xi++;
k -= 2;
sum += highFilter.Get(k2) * MAKEVAL( coeffs.Get( this->cALenExtended + xi++ ) );
k2 -= 2;
}
sigOut.Set(workIndex,
static_cast<typename OutputPortalType::ValueType>( sum ) );
sigOut.Set(workIndex, static_cast<typename OutputPortalType::ValueType>( sum ) );
}
}
@ -225,12 +219,95 @@ public:
#undef VAL
private:
vtkm::Float64 magicNum;
vtkm::Id filterLen; // filter length.
vtkm::Id cALen; // Number of actual cAs
vtkm::Id cALenExtended; // Number of extended cA at the beginning of input array
vtkm::Id cALen2; // = cALen * 2
vtkm::Id cALenExtended; // Number of cA at the beginning of input, followed by cD
}; // class ForwardTransform
}; // class InverseTransformOdd
// Worklet: perform an inverse transform for even length, symmetric filters.
class InverseTransformEven: public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn<ScalarAll>, // Input: coeffs,
// cA followed by cD
WholeArrayIn<Scalar>, // lowFilter
WholeArrayIn<Scalar>, // highFilter
WholeArrayOut<ScalarAll>); // output
typedef void ExecutionSignature(_1, _2, _3, _4, WorkIndex);
typedef _1 InputDomain;
// Constructor
VTKM_EXEC_CONT_EXPORT
InverseTransformEven( vtkm::Id filtL, vtkm::Id cAL, vtkm::Id cALExt, bool m ) :
filterLen(filtL), cALen(cAL), cALenExtended(cALExt), matlab(m)
{
this->cALen2 = cALen * 2;
}
// Use 64-bit float for convolution calculation
#define VAL vtkm::Float64
#define MAKEVAL(a) (static_cast<VAL>(a))
template <typename InputPortalType,
typename FilterPortalType,
typename OutputPortalType>
VTKM_EXEC_EXPORT
void operator()(const InputPortalType &coeffs,
const FilterPortalType &lowFilter,
const FilterPortalType &highFilter,
OutputPortalType &sigOut,
const vtkm::Id &workIndex) const
{
if( workIndex < cALen2 ) // valid calculation region
{
vtkm::Id xi; // coeff indices
vtkm::Id k; // indices for low and high filter
VAL sum = 0.0;
if( matlab || (filterLen/2) % 2 != 0 ) // odd length half filter
{
xi = workIndex / 2;
if( workIndex % 2 != 0 )
k = filterLen - 1;
else
k = filterLen - 2;
}
else
{
xi = (workIndex + 1) / 2;
if( workIndex % 2 != 0 )
k = filterLen - 2;
else
k = filterLen - 1;
}
while( k > -1 ) // k >= 0
{
sum += lowFilter.Get(k) * MAKEVAL( coeffs.Get( xi ) ) + // cA
highFilter.Get(k) * MAKEVAL( coeffs.Get( xi + cALenExtended) ); // cD
xi++;
k -= 2;
}
sigOut.Set(workIndex, static_cast<typename OutputPortalType::ValueType>( sum ) );
}
}
#undef MAKEVAL
#undef VAL
private:
vtkm::Id filterLen; // filter length.
vtkm::Id cALen; // Number of actual cAs
vtkm::Id cALen2; // = cALen * 2
vtkm::Id cALenExtended; // Number of cA at the beginning of input, followed by cD
bool matlab; // followed the naming convention from VAPOR
}; // class InverseTransformEven
// Worklet:
@ -371,7 +448,7 @@ private:
};
// Worklet:
// Worklets for signal extension no. 1
class LeftSYMHExtentionWorklet : public vtkm::worklet::WorkletMapField
{
public:
@ -401,7 +478,7 @@ private:
};
// Worklet:
// Worklets for signal extension no. 2
class LeftSYMWExtentionWorklet : public vtkm::worklet::WorkletMapField
{
public:
@ -431,7 +508,61 @@ private:
};
// Worklet:
// Worklets for signal extension no. 3
class LeftASYMHExtentionWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
WholeArrayIn < ScalarAll > ); // signal part
typedef void ExecutionSignature( _1, _2, WorkIndex );
typedef _1 InputDomain;
// Constructor
VTKM_EXEC_CONT_EXPORT
LeftASYMHExtentionWorklet( vtkm::Id len ) : addLen (len) {}
template< typename PortalOutType, typename PortalInType >
VTKM_EXEC_CONT_EXPORT
void operator()( PortalOutType &portalOut,
const PortalInType &portalIn,
const vtkm::Id &workIndex) const
{
portalOut.Set( workIndex, portalIn.Get( addLen - workIndex - 1) * (-1.0) );
}
private:
vtkm::Id addLen;
};
// Worklets for signal extension no. 4
class LeftASYMWExtentionWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
WholeArrayIn < ScalarAll > ); // signal part
typedef void ExecutionSignature( _1, _2, WorkIndex );
typedef _1 InputDomain;
// Constructor
VTKM_EXEC_CONT_EXPORT
LeftASYMWExtentionWorklet( vtkm::Id len ) : addLen (len) {}
template< typename PortalOutType, typename PortalInType >
VTKM_EXEC_CONT_EXPORT
void operator()( PortalOutType &portalOut,
const PortalInType &portalIn,
const vtkm::Id &workIndex) const
{
portalOut.Set( workIndex, portalIn.Get( addLen - workIndex ) * (-1.0) );
}
private:
vtkm::Id addLen;
};
// Worklets for signal extension no. 5
class RightSYMHExtentionWorklet : public vtkm::worklet::WorkletMapField
{
public:
@ -461,7 +592,7 @@ private:
};
// Worklet:
// Worklets for signal extension no. 6
class RightSYMWExtentionWorklet : public vtkm::worklet::WorkletMapField
{
public:
@ -491,7 +622,60 @@ private:
};
// Worklet:
// Worklets for signal extension no. 7
class RightASYMHExtentionWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
WholeArrayIn < ScalarAll > ); // signal part
typedef void ExecutionSignature( _1, _2, WorkIndex );
typedef _1 InputDomain;
// Constructor
VTKM_EXEC_CONT_EXPORT
RightASYMHExtentionWorklet ( vtkm::Id sigInl ) : sigInLen( sigInl ) {}
template< typename PortalOutType, typename PortalInType >
VTKM_EXEC_CONT_EXPORT
void operator()( PortalOutType &portalOut,
const PortalInType &portalIn,
const vtkm::Id &workIndex) const
{
portalOut.Set( workIndex, portalIn.Get( sigInLen - workIndex - 1) * (-1.0) );
}
private:
vtkm::Id sigInLen;
};
// Worklets for signal extension no. 8
class RightASYMWExtentionWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
WholeArrayIn < ScalarAll > ); // signal part
typedef void ExecutionSignature( _1, _2, WorkIndex );
typedef _1 InputDomain;
// Constructor
VTKM_EXEC_CONT_EXPORT
RightASYMWExtentionWorklet ( vtkm::Id sigInl ) : sigInLen( sigInl ) {}
template< typename PortalOutType, typename PortalInType >
VTKM_EXEC_CONT_EXPORT
void operator()( PortalOutType &portalOut,
const PortalInType &portalIn,
const vtkm::Id &workIndex) const
{
portalOut.Set( workIndex, portalIn.Get( sigInLen - workIndex - 2) * (-1.0) );
}
private:
vtkm::Id sigInLen;
};
class AssignZeroWorklet : public vtkm::worklet::WorkletMapField
{
public: