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

remove unnecessary files

Browse Source
This commit is contained in:
Samuel Li 2016-07-29 13:05:53 -07:00
parent ec46fd58a0
commit 3cc6385426
7 changed files with 0 additions and 1272 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

88
vtkm/filter/WaveletCompressor.h
View File

@ -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 <vtkm/filter/internal/WaveletDWT.h>
//#include <vtkm/worklet/WaveletTransforms.h>
#include <vtkm/cont/ArrayHandleConcatenate.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/ArrayHandlePermutation.h>
#include <vtkm/Math.h>
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 <vtkm/filter/WaveletCompressor.hxx>
#endif

225
vtkm/filter/WaveletCompressor.hxx
View File

@ -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 <vtkm/filter/WaveletCompressor.h>
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

80
vtkm/filter/internal/FilterBanks.h
View File

@ -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 <vtkm/Types.h>
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

183
vtkm/filter/internal/WaveletBase.h
View File

@ -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 <vtkm/filter/internal/WaveletFilter.h>
#include <vtkm/Math.h>
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::Id>(vtkm::Ceil( (static_cast<vtkm::Float64>(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::Id>( vtkm::Floor(
static_cast<vtkm::Float64>(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::Id>(vtkm::Ceil( (static_cast<vtkm::Float64>(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::Id>( vtkm::Floor(
static_cast<vtkm::Float64>(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::Id>( vtkm::Floor(
vtkm::Log2( static_cast<vtkm::Float64>(sigInLen) /
static_cast<vtkm::Float64>(filterLength) ) + 1.0 ) );
}
}; // Finish class WaveletBase.
} // Finish namespace internal
} // Finish namespace filter
} // Finish namespace vtkm
#endif

410
vtkm/filter/internal/WaveletDWT.h
View File

@ -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 <vtkm/worklet/WorkletMapField.h>
//#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/filter/internal/WaveletBase.h>
#include <vtkm/worklet/WaveletTransforms.h>
#include <vtkm/cont/ArrayHandleConcatenate.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/ArrayHandlePermutation.h>
#include <vtkm/Math.h>
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<SigInValueType> 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<vtkm::worklet::ForwardTransform>
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

182
vtkm/filter/internal/WaveletFilter.h
View File

@ -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 <vtkm/cont/ArrayHandle.h>
#include <vtkm/filter/internal/FilterBanks.h>
#include <vtkm/Math.h>
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<vtkm::Float64> 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

104
vtkm/worklet/testing/UnitTestWaveletTransforms.cxx
View File

@ -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 <vtkm/worklet/WaveletTransforms.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/filter/internal/FilterBanks.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/cont/Timer.h>
#include <vector>
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<vtkm::Float64> tmpVector;
for( vtkm::Id i = 0; i < sigLen + 8; i++ )
tmpVector.push_back( static_cast<vtkm::Float64>(i%100+1) );
vtkm::cont::ArrayHandle<vtkm::Float64> input1DArray =
vtkm::cont::make_ArrayHandle(tmpVector);
// output array handle
vtkm::cont::ArrayHandle<vtkm::Float64> outputArray1;
// make two filter array handles
vtkm::cont::ArrayHandle<vtkm::Float64> lowFilter =
vtkm::cont::make_ArrayHandle(vtkm::filter::internal::hm4_44, 9);
vtkm::cont::ArrayHandle<vtkm::Float64> 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<vtkm::worklet::ForwardTransform>
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;
}