diff --git a/vtkm/worklet/wavelets/FilterBanks.h b/vtkm/worklet/wavelets/FilterBanks.h index 4814174d9..124a64e15 100644 --- a/vtkm/worklet/wavelets/FilterBanks.h +++ b/vtkm/worklet/wavelets/FilterBanks.h @@ -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 + } ; + }; } diff --git a/vtkm/worklet/wavelets/WaveletBase.h b/vtkm/worklet/wavelets/WaveletBase.h index a8c7040b7..1d5ef116f 100644 --- a/vtkm/worklet/wavelets/WaveletBase.h +++ b/vtkm/worklet/wavelets/WaveletBase.h @@ -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::Ceil( (static_cast(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::Ceil( (static_cast(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); diff --git a/vtkm/worklet/wavelets/WaveletDWT.h b/vtkm/worklet/wavelets/WaveletDWT.h index 05b7f65a7..a9e99578f 100644 --- a/vtkm/worklet/wavelets/WaveletDWT.h +++ b/vtkm/worklet/wavelets/WaveletDWT.h @@ -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 ) diff --git a/vtkm/worklet/wavelets/WaveletFilter.h b/vtkm/worklet/wavelets/WaveletFilter.h index 1abb9767f..1d3db9b28 100644 --- a/vtkm/worklet/wavelets/WaveletFilter.h +++ b/vtkm/worklet/wavelets/WaveletFilter.h @@ -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 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.