Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

half done even length filters implementation

Browse Source
This commit is contained in:
Samuel Li 2016-08-15 13:46:35 -06:00
parent ee32ea4cf9
commit a6efad0448
4 changed files with 115 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -72,6 +72,18 @@ namespace wavelets {
0.0
};
const double hm1_11[2] = {
0.70710678118654752440084436210,
0.70710678118654752440084436210
};
const double h1[10] = {
0.0, 0.0, 0.0, 0.0,
0.70710678118654752440084436210,
0.70710678118654752440084436210,
0.0, 0.0, 0.0, 0.0
} ;
};
}

54
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.
@ -55,13 +49,14 @@ public:
WaveletBase( WaveletName name ) : wname ( name ),
filter( name )
{
if( this->wname == CDF9_7 )
if( wname == CDF9_7 || wname == BIOR4_4 ||
wname == CDF5_3 || wname == BIOR2_2 )
{
this->wmode = SYMW; // Default extension mode, see MatWaveBase.cpp
}
else if( this->wname == CDF5_3 )
else if( wname == HAAR || wname == BIOR1_1 )
{
this->wmode = SYMW;
this->wmode = SYMH;
}
}
@ -70,9 +65,7 @@ public:
{
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)) )
@ -93,9 +86,7 @@ public:
{
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)) )
@ -136,7 +127,7 @@ public:
// 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)
@ -147,18 +138,29 @@ 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, typename DeviceTag >
VTKM_EXEC_CONT_EXPORT
VTKM_EXEC_EXPORT
void DeviceSort( ArrayType &array, DeviceTag )
{
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Sort
@ -167,7 +169,7 @@ public:
// Reduce to the sum of all values on device
template< typename ArrayType, typename DeviceTag >
VTKM_EXEC_CONT_EXPORT
VTKM_EXEC_EXPORT
typename ArrayType::ValueType DeviceSum( const ArrayType &array, DeviceTag )
{
return vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Reduce
@ -178,7 +180,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);
}
@ -186,13 +188,13 @@ 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, typename DeviceTag >
VTKM_EXEC_CONT_EXPORT
VTKM_EXEC_EXPORT
typename ArrayType::ValueType DeviceMax( const ArrayType &array, DeviceTag )
{
typename ArrayType::ValueType initVal = array.GetPortalConstControl().Get(0);
@ -200,7 +202,7 @@ public:
( array, initVal, maxFunctor() );
}
template< typename ArrayType, typename DeviceTag >
VTKM_EXEC_CONT_EXPORT
VTKM_EXEC_EXPORT
typename ArrayType::ValueType DeviceMin( const ArrayType &array, DeviceTag )
{
typename ArrayType::ValueType initVal = array.GetPortalConstControl().Get(0);
@ -212,13 +214,13 @@ 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, typename DeviceTag >
VTKM_EXEC_CONT_EXPORT
VTKM_EXEC_EXPORT
typename ArrayType::ValueType DeviceMaxAbs( const ArrayType &array, DeviceTag )
{
typename ArrayType::ValueType initVal = array.GetPortalConstControl().Get(0);

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

@ -52,8 +52,11 @@ 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;
@ -90,9 +93,9 @@ public:
}
ExtensionArrayType rightExtend;
if( attachRightZero )
if( attachZeroRightRight )
rightExtend.Allocate( addLen + 1 );
else
else // attachZeroRightLeft mode will attach zero later.
rightExtend.Allocate( addLen );
switch( rightExtMethod )
@ -118,12 +121,15 @@ public:
}
}
if( attachRightZero )
if( attachZeroRightRight )
WaveletBase::DeviceAssignZero( rightExtend, addLen );
else if( attachZeroRightLeft )
{
typedef vtkm::worklet::wavelets::AssignZeroWorklet ZeroWorklet;
ZeroWorklet worklet( addLen );
vtkm::worklet::DispatcherMapField< ZeroWorklet > dispatcher( worklet );
dispatcher.Invoke( rightExtend );
ExtensionArrayType rightExtendPlusOne;
rightExtendPlusOne.Allocate( addLen + 1 );
WaveletBase::DeviceCopyStartX( rightExtend, rightExtendPlusOne, 1 );
WaveletBase::DeviceAssignZero( rightExtendPlusOne, 0 );
rightExtend = rightExtendPlusOne ;
}
typedef vtkm::cont::ArrayHandleConcatenate< ExtensionArrayType, SigInArrayType>
@ -197,10 +203,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] );
@ -231,7 +237,7 @@ 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();
bool doSymConv = false;
@ -279,15 +285,16 @@ public:
{
addLen = filterLen / 4;
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++;
@ -314,27 +321,20 @@ 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 );
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, at the beginning of right 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,11 +342,19 @@ 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 );
}
if( filterLen % 2 != 0 )

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

@ -34,7 +34,11 @@ namespace wavelets {
enum WaveletName {
CDF9_7,
CDF5_3
CDF5_3,
HAAR,
BIOR4_4, // the same as CDF9_7
BIOR2_2, // the same as CDF5_3
BIOR1_1 // the same as HAAE
};
// Flipping operation; helper function to initialize a filter.
@ -101,37 +105,48 @@ public:
lowReconstructFilter(NULL),
highReconstructFilter(NULL)
{
if( wtype == CDF9_7 )
if( wtype == CDF9_7 || wtype == BIOR4_4 )
{
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 );
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( wtype == CDF5_3 )
else if( wtype == CDF5_3 || wtype == BIOR2_2 )
{
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 );
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 if( wtype == HAAR || wtype == BIOR1_1 )
{
this->filterLength = 2;
this->AllocateFilterMemory();
wrev( vtkm::worklet::wavelets::hm1_11, lowDecomposeFilter, filterLength );
qmf_wrev( vtkm::worklet::wavelets::h1+4, highDecomposeFilter, filterLength );
verbatim_copy( vtkm::worklet::wavelets::h1+4, lowReconstructFilter, filterLength );
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
@ -161,10 +176,10 @@ private:
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;
}
}; // class WaveletFilter.
} // namespace wavelets.