731 lines
26 KiB
C++
731 lines
26 KiB
C++
//============================================================================
|
|
// 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_worklet_waveletcompressor_h
|
|
#define vtk_m_worklet_waveletcompressor_h
|
|
|
|
#include <vtkm/worklet/wavelets/WaveletDWT.h>
|
|
|
|
namespace vtkm {
|
|
namespace worklet {
|
|
|
|
class WaveletCompressor : public vtkm::worklet::wavelets::WaveletDWT
|
|
{
|
|
public:
|
|
|
|
// Constructor
|
|
WaveletCompressor( wavelets::WaveletName name ) : WaveletDWT( name ) {}
|
|
|
|
|
|
// Multi-level 1D wavelet decomposition
|
|
template< typename SignalArrayType, typename CoeffArrayType, typename DeviceTag >
|
|
VTKM_CONT
|
|
vtkm::Id WaveDecompose( const SignalArrayType &sigIn, // Input
|
|
vtkm::Id nLevels, // n levels of DWT
|
|
CoeffArrayType &coeffOut,
|
|
std::vector<vtkm::Id> &L,
|
|
DeviceTag )
|
|
{
|
|
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
|
|
if( nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel( sigInLen ) )
|
|
{
|
|
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
|
|
}
|
|
if( nLevels == 0 ) // 0 levels means no transform
|
|
{
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
|
|
return 0;
|
|
}
|
|
|
|
this->ComputeL( sigInLen, nLevels, L ); // memory for L is allocated by ComputeL().
|
|
vtkm::Id CLength = this->ComputeCoeffLength( L, nLevels );
|
|
VTKM_ASSERT( CLength == sigInLen );
|
|
|
|
vtkm::Id sigInPtr = 0; // pseudo pointer for the beginning of input array
|
|
vtkm::Id len = sigInLen;
|
|
vtkm::Id cALen = WaveletBase::GetApproxLength( len );
|
|
vtkm::Id cptr; // pseudo pointer for the beginning of output array
|
|
vtkm::Id tlen = 0;
|
|
std::vector<vtkm::Id> L1d(3, 0);
|
|
|
|
// Use an intermediate array
|
|
typedef typename CoeffArrayType::ValueType OutputValueType;
|
|
typedef vtkm::cont::ArrayHandle< OutputValueType > InterArrayType;
|
|
|
|
// Define a few more types
|
|
typedef vtkm::cont::ArrayHandleCounting< vtkm::Id > IdArrayType;
|
|
typedef vtkm::cont::ArrayHandlePermutation< IdArrayType, CoeffArrayType >
|
|
PermutArrayType;
|
|
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
|
|
|
|
for( vtkm::Id i = nLevels; i > 0; i-- )
|
|
{
|
|
tlen += L[ size_t(i) ];
|
|
cptr = 0 + CLength - tlen - cALen;
|
|
|
|
// make input array (permutation array)
|
|
IdArrayType inputIndices( sigInPtr, 1, len );
|
|
PermutArrayType input( inputIndices, coeffOut );
|
|
// make output array
|
|
InterArrayType output;
|
|
|
|
WaveletDWT::DWT1D( input, output, L1d, DeviceTag() );
|
|
|
|
// move intermediate results to final array
|
|
WaveletBase::DeviceCopyStartX( output, coeffOut, cptr, DeviceTag() );
|
|
|
|
// update pseudo pointers
|
|
len = cALen;
|
|
cALen = WaveletBase::GetApproxLength( cALen );
|
|
sigInPtr = cptr;
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
|
|
// Multi-level 1D wavelet reconstruction
|
|
template< typename CoeffArrayType, typename SignalArrayType, typename DeviceTag >
|
|
VTKM_CONT
|
|
vtkm::Id WaveReconstruct( const CoeffArrayType &coeffIn, // Input
|
|
vtkm::Id nLevels, // n levels of DWT
|
|
std::vector<vtkm::Id> &L,
|
|
SignalArrayType &sigOut,
|
|
DeviceTag )
|
|
{
|
|
VTKM_ASSERT( nLevels > 0 );
|
|
vtkm::Id LLength = nLevels + 2;
|
|
VTKM_ASSERT( vtkm::Id(L.size()) == LLength );
|
|
|
|
std::vector<vtkm::Id> L1d(3, 0); // three elements
|
|
L1d[0] = L[0];
|
|
L1d[1] = L[1];
|
|
|
|
typedef typename SignalArrayType::ValueType OutValueType;
|
|
typedef vtkm::cont::ArrayHandle< OutValueType > OutArrayBasic;
|
|
typedef vtkm::cont::ArrayHandleCounting< vtkm::Id > IdArrayType;
|
|
typedef vtkm::cont::ArrayHandlePermutation< IdArrayType, SignalArrayType >
|
|
PermutArrayType;
|
|
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( coeffIn, sigOut );
|
|
|
|
for( vtkm::Id i = 1; i <= nLevels; i++ )
|
|
{
|
|
L1d[2] = this->GetApproxLengthLevN( L[ size_t(LLength-1) ], nLevels-i );
|
|
|
|
// Make an input array
|
|
IdArrayType inputIndices( 0, 1, L1d[2] );
|
|
PermutArrayType input( inputIndices, sigOut );
|
|
|
|
// Make an output array
|
|
OutArrayBasic output;
|
|
|
|
WaveletDWT::IDWT1D( input, L1d, output, DeviceTag() );
|
|
VTKM_ASSERT( output.GetNumberOfValues() == L1d[2] );
|
|
|
|
// Move output to intermediate array
|
|
WaveletBase::DeviceCopyStartX( output, sigOut, 0, DeviceTag() );
|
|
|
|
L1d[0] = L1d[2];
|
|
L1d[1] = L[ size_t(i+1) ];
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
|
|
// Multi-level 3D wavelet decomposition
|
|
template< typename InArrayType, typename OutArrayType, typename DeviceTag>
|
|
VTKM_CONT
|
|
vtkm::Float64 WaveDecompose3D(
|
|
InArrayType &sigIn, // Input
|
|
vtkm::Id nLevels, // n levels of DWT
|
|
vtkm::Id inX, vtkm::Id inY, vtkm::Id inZ,
|
|
OutArrayType &coeffOut,
|
|
bool discardSigIn, // can we discard sigIn for more memory?
|
|
DeviceTag )
|
|
{
|
|
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
|
|
VTKM_ASSERT( inX * inY * inZ == sigInLen );
|
|
if( nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel( inX ) ||
|
|
nLevels > WaveletBase::GetWaveletMaxLevel( inY ) ||
|
|
nLevels > WaveletBase::GetWaveletMaxLevel( inZ ) )
|
|
{
|
|
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
|
|
}
|
|
if( nLevels == 0 ) // 0 levels means no transform
|
|
{
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
|
|
return 0;
|
|
}
|
|
|
|
vtkm::Id currentLenX = inX;
|
|
vtkm::Id currentLenY = inY;
|
|
vtkm::Id currentLenZ = inZ;
|
|
std::vector<vtkm::Id> L3d(27, 0);
|
|
|
|
typedef typename OutArrayType::ValueType OutValueType;
|
|
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
|
|
|
|
// First level transform writes to the output array
|
|
vtkm::Float64 computationTime = WaveletDWT::DWT3D(
|
|
sigIn,
|
|
inX, inY, inZ,
|
|
0, 0, 0,
|
|
currentLenX, currentLenY, currentLenZ,
|
|
coeffOut,
|
|
discardSigIn,
|
|
DeviceTag() );
|
|
|
|
// Successor transforms writes to a temporary array
|
|
for( vtkm::Id i = nLevels-1; i > 0; i-- )
|
|
{
|
|
currentLenX = WaveletBase::GetApproxLength( currentLenX );
|
|
currentLenY = WaveletBase::GetApproxLength( currentLenY );
|
|
currentLenZ = WaveletBase::GetApproxLength( currentLenZ );
|
|
|
|
OutBasicArray tempOutput;
|
|
|
|
computationTime += WaveletDWT::DWT3D(
|
|
coeffOut,
|
|
inX, inY, inZ,
|
|
0, 0, 0,
|
|
currentLenX, currentLenY, currentLenZ,
|
|
tempOutput,
|
|
false,
|
|
DeviceTag() );
|
|
|
|
// copy results to coeffOut
|
|
WaveletBase::DeviceCubeCopyTo( tempOutput,
|
|
currentLenX, currentLenY, currentLenZ,
|
|
coeffOut,
|
|
inX, inY, inZ,
|
|
0, 0, 0,
|
|
DeviceTag() );
|
|
}
|
|
|
|
return computationTime;
|
|
}
|
|
|
|
|
|
|
|
// Multi-level 3D wavelet reconstruction
|
|
template< typename InArrayType, typename OutArrayType, typename DeviceTag>
|
|
VTKM_CONT
|
|
vtkm::Float64 WaveReconstruct3D(
|
|
InArrayType &arrIn, // Input
|
|
vtkm::Id nLevels, // n levels of DWT
|
|
vtkm::Id inX, vtkm::Id inY, vtkm::Id inZ,
|
|
OutArrayType &arrOut,
|
|
bool discardArrIn, // can we discard input for more memory?
|
|
DeviceTag )
|
|
{
|
|
vtkm::Id arrInLen = arrIn.GetNumberOfValues();
|
|
VTKM_ASSERT( inX * inY * inZ == arrInLen );
|
|
if( nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel( inX ) ||
|
|
nLevels > WaveletBase::GetWaveletMaxLevel( inY ) ||
|
|
nLevels > WaveletBase::GetWaveletMaxLevel( inZ ) )
|
|
{
|
|
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
|
|
}
|
|
typedef typename OutArrayType::ValueType OutValueType;
|
|
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
|
|
vtkm::Float64 computationTime = 0.0;
|
|
|
|
OutBasicArray outBuffer;
|
|
if( nLevels == 0 ) // 0 levels means no transform
|
|
{
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, arrOut );
|
|
return 0;
|
|
}
|
|
else if ( discardArrIn )
|
|
outBuffer = arrIn;
|
|
else
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, outBuffer );
|
|
|
|
std::vector<vtkm::Id> L;
|
|
this->ComputeL3( inX, inY, inZ, nLevels, L );
|
|
std::vector<vtkm::Id> L3d(27, 0);
|
|
|
|
// All transforms but the last level operate on temporary arrays
|
|
for( size_t i = 0; i < 24; i++ )
|
|
L3d[i] = L[i];
|
|
for( size_t i = 1; i < static_cast<size_t>(nLevels); i++ )
|
|
{
|
|
L3d[24] = L3d[0] + L3d[12]; // Total X dim; this is always true for Biorthogonal wavelets
|
|
L3d[25] = L3d[1] + L3d[7]; // Total Y dim
|
|
L3d[26] = L3d[2] + L3d[5]; // Total Z dim
|
|
|
|
OutBasicArray tempOutput;
|
|
|
|
// IDWT
|
|
computationTime += WaveletDWT::IDWT3D( outBuffer,
|
|
inX, inY, inZ,
|
|
0, 0, 0,
|
|
L3d,
|
|
tempOutput,
|
|
false,
|
|
DeviceTag() );
|
|
|
|
// copy back reconstructed block
|
|
WaveletBase::DeviceCubeCopyTo( tempOutput,
|
|
L3d[24], L3d[25], L3d[26],
|
|
outBuffer,
|
|
inX, inY, inZ,
|
|
0, 0, 0,
|
|
DeviceTag() );
|
|
|
|
// update L3d array
|
|
L3d[0] = L3d[24];
|
|
L3d[1] = L3d[25];
|
|
L3d[2] = L3d[26];
|
|
for( size_t j = 3; j < 24; j++ )
|
|
L3d[j] = L[ 21 * i + j ];
|
|
}
|
|
|
|
// The last transform outputs to the final output
|
|
L3d[24] = L3d[0] + L3d[12];
|
|
L3d[25] = L3d[1] + L3d[7];
|
|
L3d[26] = L3d[2] + L3d[5];
|
|
computationTime += WaveletDWT::IDWT3D( outBuffer,
|
|
inX, inY, inZ,
|
|
0, 0, 0,
|
|
L3d,
|
|
arrOut,
|
|
true,
|
|
DeviceTag() );
|
|
|
|
return computationTime;
|
|
}
|
|
|
|
|
|
|
|
// Multi-level 2D wavelet decomposition
|
|
template< typename InArrayType, typename OutArrayType, typename DeviceTag>
|
|
VTKM_CONT
|
|
vtkm::Float64 WaveDecompose2D( const InArrayType &sigIn, // Input
|
|
vtkm::Id nLevels, // n levels of DWT
|
|
vtkm::Id inX, // Input X dim
|
|
vtkm::Id inY, // Input Y dim
|
|
OutArrayType &coeffOut,
|
|
std::vector<vtkm::Id> &L,
|
|
DeviceTag )
|
|
{
|
|
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
|
|
VTKM_ASSERT( inX * inY == sigInLen );
|
|
if( nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel( inX ) ||
|
|
nLevels > WaveletBase::GetWaveletMaxLevel( inY ) )
|
|
{
|
|
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
|
|
}
|
|
if( nLevels == 0 ) // 0 levels means no transform
|
|
{
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
|
|
return 0;
|
|
}
|
|
|
|
this->ComputeL2( inX, inY, nLevels, L );
|
|
vtkm::Id CLength = this->ComputeCoeffLength2( L, nLevels );
|
|
VTKM_ASSERT( CLength == sigInLen );
|
|
|
|
vtkm::Id currentLenX = inX;
|
|
vtkm::Id currentLenY = inY;
|
|
std::vector<vtkm::Id> L2d(10, 0);
|
|
vtkm::Float64 computationTime = 0.0;
|
|
|
|
typedef typename OutArrayType::ValueType OutValueType;
|
|
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
|
|
|
|
// First level transform operates writes to the output array
|
|
computationTime += WaveletDWT::DWT2D ( sigIn,
|
|
currentLenX, currentLenY,
|
|
0, 0,
|
|
currentLenX, currentLenY,
|
|
coeffOut, L2d, DeviceTag() );
|
|
VTKM_ASSERT( coeffOut.GetNumberOfValues() == currentLenX * currentLenY );
|
|
currentLenX = WaveletBase::GetApproxLength( currentLenX );
|
|
currentLenY = WaveletBase::GetApproxLength( currentLenY );
|
|
|
|
// Successor transforms writes to a temporary array
|
|
for( vtkm::Id i = nLevels-1; i > 0; i-- )
|
|
{
|
|
OutBasicArray tempOutput;
|
|
|
|
computationTime +=
|
|
WaveletDWT::DWT2D ( coeffOut,
|
|
inX, inY,
|
|
0, 0,
|
|
currentLenX, currentLenY,
|
|
tempOutput, L2d, DeviceTag() );
|
|
|
|
// copy results to coeffOut
|
|
WaveletBase::DeviceRectangleCopyTo( tempOutput, currentLenX, currentLenY,
|
|
coeffOut, inX, inY, 0, 0, DeviceTag() );
|
|
|
|
// update currentLen
|
|
currentLenX = WaveletBase::GetApproxLength( currentLenX );
|
|
currentLenY = WaveletBase::GetApproxLength( currentLenY );
|
|
}
|
|
|
|
return computationTime;
|
|
}
|
|
|
|
|
|
|
|
// Multi-level 2D wavelet reconstruction
|
|
template< typename InArrayType, typename OutArrayType, typename DeviceTag>
|
|
VTKM_CONT
|
|
vtkm::Float64 WaveReconstruct2D( const InArrayType &arrIn, // Input
|
|
vtkm::Id nLevels, // n levels of DWT
|
|
vtkm::Id inX, // Input X dim
|
|
vtkm::Id inY, // Input Y dim
|
|
OutArrayType &arrOut,
|
|
std::vector<vtkm::Id> &L,
|
|
DeviceTag )
|
|
{
|
|
vtkm::Id arrInLen = arrIn.GetNumberOfValues();
|
|
VTKM_ASSERT( inX * inY == arrInLen );
|
|
if( nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel( inX ) ||
|
|
nLevels > WaveletBase::GetWaveletMaxLevel( inY ) )
|
|
{
|
|
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
|
|
}
|
|
typedef typename OutArrayType::ValueType OutValueType;
|
|
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
|
|
vtkm::Float64 computationTime = 0.0;
|
|
|
|
OutBasicArray outBuffer;
|
|
if( nLevels == 0 ) // 0 levels means no transform
|
|
{
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, arrOut );
|
|
return 0;
|
|
}
|
|
else
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, outBuffer );
|
|
|
|
VTKM_ASSERT( vtkm::Id(L.size()) == 6 * nLevels + 4 );
|
|
|
|
std::vector<vtkm::Id> L2d(10, 0);
|
|
L2d[0] = L[0];
|
|
L2d[1] = L[1];
|
|
L2d[2] = L[2];
|
|
L2d[3] = L[3];
|
|
L2d[4] = L[4];
|
|
L2d[5] = L[5];
|
|
L2d[6] = L[6];
|
|
L2d[7] = L[7];
|
|
|
|
// All transforms but the last operate on temporary arrays
|
|
for( size_t i = 1; i < static_cast<size_t>(nLevels); i++ )
|
|
{
|
|
L2d[8] = L2d[0] + L2d[4]; // This is always true for Biorthogonal wavelets
|
|
L2d[9] = L2d[1] + L2d[3]; // (same above)
|
|
|
|
OutBasicArray tempOutput;
|
|
|
|
// IDWT
|
|
computationTime +=
|
|
WaveletDWT::IDWT2D ( outBuffer, inX, inY, 0, 0, L2d, tempOutput, DeviceTag() );
|
|
|
|
// copy back reconstructed block
|
|
WaveletBase::DeviceRectangleCopyTo( tempOutput, L2d[8], L2d[9],
|
|
outBuffer, inX, inY, 0, 0, DeviceTag() );
|
|
|
|
// update L2d array
|
|
L2d[0] = L2d[8];
|
|
L2d[1] = L2d[9];
|
|
L2d[2] = L[6*i+2];
|
|
L2d[3] = L[6*i+3];
|
|
L2d[4] = L[6*i+4];
|
|
L2d[5] = L[6*i+5];
|
|
L2d[6] = L[6*i+6];
|
|
L2d[7] = L[6*i+7];
|
|
|
|
}
|
|
|
|
// The last transform outputs to the final output
|
|
L2d[8] = L2d[0] + L2d[4];
|
|
L2d[9] = L2d[1] + L2d[3];
|
|
computationTime +=
|
|
WaveletDWT::IDWT2D ( outBuffer, inX, inY, 0, 0, L2d, arrOut, DeviceTag() );
|
|
|
|
return computationTime;
|
|
}
|
|
|
|
|
|
|
|
// Squash coefficients smaller than a threshold
|
|
template< typename CoeffArrayType, typename DeviceTag >
|
|
vtkm::Id SquashCoefficients( CoeffArrayType &coeffIn,
|
|
vtkm::Float64 ratio,
|
|
DeviceTag )
|
|
{
|
|
if( ratio > 1.0 )
|
|
{
|
|
vtkm::Id coeffLen = coeffIn.GetNumberOfValues();
|
|
typedef typename CoeffArrayType::ValueType ValueType;
|
|
typedef vtkm::cont::ArrayHandle< ValueType > CoeffArrayBasic;
|
|
CoeffArrayBasic 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 );
|
|
vtkm::Float64 nthVal = static_cast<vtkm::Float64>
|
|
(sortedArray.GetPortalConstControl().Get(n));
|
|
if( nthVal < 0.0 )
|
|
nthVal *= -1.0;
|
|
typedef vtkm::worklet::wavelets::ThresholdWorklet ThresholdType;
|
|
ThresholdType thresholdWorklet( nthVal );
|
|
vtkm::worklet::DispatcherMapField< ThresholdType, DeviceTag > dispatcher( thresholdWorklet );
|
|
dispatcher.Invoke( coeffIn );
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
|
|
// Report statistics on reconstructed array
|
|
template< typename ArrayType, typename DeviceTag >
|
|
vtkm::Id EvaluateReconstruction( const ArrayType &original,
|
|
const ArrayType &reconstruct,
|
|
DeviceTag )
|
|
{
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
VAL VarOrig = WaveletBase::DeviceCalculateVariance( original, DeviceTag() );
|
|
|
|
typedef typename ArrayType::ValueType ValueType;
|
|
typedef vtkm::cont::ArrayHandle< ValueType > ArrayBasic;
|
|
ArrayBasic errorArray, errorSquare;
|
|
|
|
// Use a worklet to calculate point-wise error, and its square
|
|
typedef vtkm::worklet::wavelets::Differencer DifferencerWorklet;
|
|
DifferencerWorklet dw;
|
|
vtkm::worklet::DispatcherMapField< DifferencerWorklet > dwDispatcher( dw );
|
|
dwDispatcher.Invoke( original, reconstruct, errorArray );
|
|
|
|
typedef vtkm::worklet::wavelets::SquareWorklet SquareWorklet;
|
|
SquareWorklet sw;
|
|
vtkm::worklet::DispatcherMapField< SquareWorklet > swDispatcher( sw );
|
|
swDispatcher.Invoke( errorArray, errorSquare );
|
|
|
|
VAL varErr = WaveletBase::DeviceCalculateVariance( errorArray, DeviceTag() );
|
|
VAL snr, decibels;
|
|
if( varErr != 0.0 )
|
|
{
|
|
snr = VarOrig / varErr;
|
|
decibels = 10 * vtkm::Log10( snr );
|
|
}
|
|
else
|
|
{
|
|
snr = vtkm::Infinity64();
|
|
decibels = vtkm::Infinity64();
|
|
}
|
|
|
|
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, DeviceTag() );
|
|
VAL rmse = vtkm::Sqrt( MAKEVAL(squareSum) / MAKEVAL(errorArray.GetNumberOfValues()) );
|
|
|
|
std::cout << "Data range = " << range << std::endl;
|
|
std::cout << "SNR = " << snr << std::endl;
|
|
std::cout << "SNR in decibels = " << decibels << std::endl;
|
|
std::cout << "L-infy norm = " << errorMax
|
|
<< ", after normalization = " << errorMax / range << std::endl;
|
|
std::cout << "RMSE = " << rmse
|
|
<< ", after normalization = " << rmse / range << std::endl;
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
|
|
// Compute the book keeping array L for 1D DWT
|
|
void ComputeL( vtkm::Id sigInLen,
|
|
vtkm::Id nLev,
|
|
std::vector<vtkm::Id> &L )
|
|
{
|
|
size_t nLevels = static_cast<size_t>( nLev ); // cast once
|
|
L.resize( nLevels + 2 );
|
|
L[ nLevels+1 ] = sigInLen;
|
|
L[ nLevels ] = sigInLen;
|
|
for( size_t i = nLevels; i > 0; i-- )
|
|
{
|
|
L[i-1] = WaveletBase::GetApproxLength( L[i] );
|
|
L[i] = WaveletBase::GetDetailLength( L[i] );
|
|
}
|
|
}
|
|
|
|
|
|
|
|
// Compute the book keeping array L for 2D DWT
|
|
void ComputeL2( vtkm::Id inX,
|
|
vtkm::Id inY,
|
|
vtkm::Id nLev,
|
|
std::vector<vtkm::Id> &L )
|
|
{
|
|
size_t nLevels = static_cast<size_t>( nLev );
|
|
L.resize( nLevels*6 + 4 );
|
|
L[ nLevels*6 ] = inX;
|
|
L[ nLevels*6 + 1 ] = inY;
|
|
L[ nLevels*6 + 2 ] = inX;
|
|
L[ nLevels*6 + 3 ] = inY;
|
|
|
|
for( size_t i = nLevels; i > 0; i-- )
|
|
{
|
|
// cA
|
|
L[ i*6 - 6 ] = WaveletBase::GetApproxLength( L[ i*6 + 0 ]);
|
|
L[ i*6 - 5 ] = WaveletBase::GetApproxLength( L[ i*6 + 1 ]);
|
|
|
|
// cDh
|
|
L[ i*6 - 4 ] = WaveletBase::GetApproxLength( L[ i*6 + 0 ]);
|
|
L[ i*6 - 3 ] = WaveletBase::GetDetailLength( L[ i*6 + 1 ]);
|
|
|
|
// cDv
|
|
L[ i*6 - 2 ] = WaveletBase::GetDetailLength( L[ i*6 + 0 ]);
|
|
L[ i*6 - 1 ] = WaveletBase::GetApproxLength( L[ i*6 + 1 ]);
|
|
|
|
// cDv - overwrites previous value!
|
|
L[ i*6 - 0 ] = WaveletBase::GetDetailLength( L[ i*6 + 0 ]);
|
|
L[ i*6 + 1 ] = WaveletBase::GetDetailLength( L[ i*6 + 1 ]);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
// Compute the bookkeeping array L for 3D DWT
|
|
void ComputeL3( vtkm::Id inX,
|
|
vtkm::Id inY,
|
|
vtkm::Id inZ,
|
|
vtkm::Id nLev,
|
|
std::vector<vtkm::Id> &L )
|
|
{
|
|
size_t n = static_cast<size_t>( nLev );
|
|
L.resize( n * 21 + 6 );
|
|
L[ n * 21 + 0 ] = inX;
|
|
L[ n * 21 + 1 ] = inY;
|
|
L[ n * 21 + 2 ] = inZ;
|
|
L[ n * 21 + 3 ] = inX;
|
|
L[ n * 21 + 4 ] = inY;
|
|
L[ n * 21 + 5 ] = inZ;
|
|
|
|
for( size_t i = n; i > 0; i-- )
|
|
{
|
|
// cLLL
|
|
L[ i * 21 - 21 ] = WaveletBase::GetApproxLength( L[ i * 21 + 0 ] );
|
|
L[ i * 21 - 20 ] = WaveletBase::GetApproxLength( L[ i * 21 + 1 ] );
|
|
L[ i * 21 - 19 ] = WaveletBase::GetApproxLength( L[ i * 21 + 2 ] );
|
|
|
|
// cLLH
|
|
L[ i * 21 - 18 ] = L[ i * 21 - 21 ];
|
|
L[ i * 21 - 17 ] = L[ i * 21 - 20 ];
|
|
L[ i * 21 - 16 ] = WaveletBase::GetDetailLength( L[ i * 21 + 2 ] );
|
|
|
|
// cLHL
|
|
L[ i * 21 - 15 ] = L[ i * 21 - 21 ];
|
|
L[ i * 21 - 14 ] = WaveletBase::GetDetailLength( L[ i * 21 + 1 ] );
|
|
L[ i * 21 - 13 ] = L[ i * 21 - 19 ];
|
|
|
|
// cLHH
|
|
L[ i * 21 - 12 ] = L[ i * 21 - 21 ];
|
|
L[ i * 21 - 11 ] = L[ i * 21 - 14 ];
|
|
L[ i * 21 - 10 ] = L[ i * 21 - 16 ];
|
|
|
|
// cHLL
|
|
L[ i * 21 - 9 ] = WaveletBase::GetDetailLength( L[ i * 21 + 0 ] );
|
|
L[ i * 21 - 8 ] = L[ i * 21 - 20 ];
|
|
L[ i * 21 - 7 ] = L[ i * 21 - 19 ];
|
|
|
|
// cHLH
|
|
L[ i * 21 - 6 ] = L[ i * 21 - 9 ];
|
|
L[ i * 21 - 5 ] = L[ i * 21 - 20 ];
|
|
L[ i * 21 - 3 ] = L[ i * 21 - 16 ];
|
|
|
|
// cHHL
|
|
L[ i * 21 - 3 ] = L[ i * 21 - 9 ];
|
|
L[ i * 21 - 2 ] = L[ i * 21 - 14 ];
|
|
L[ i * 21 - 1 ] = L[ i * 21 - 19 ];
|
|
|
|
// cHHH - overwrites previous value
|
|
L[ i * 21 + 0 ] = L[ i * 21 - 9 ];
|
|
L[ i * 21 + 1 ] = L[ i * 21 - 14 ];
|
|
L[ i * 21 + 2 ] = L[ i * 21 - 16 ];
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// Compute the length of coefficients for 1D transforms
|
|
vtkm::Id ComputeCoeffLength( std::vector<vtkm::Id> &L,
|
|
vtkm::Id nLevels )
|
|
{
|
|
vtkm::Id sum = L[0]; // 1st level cA
|
|
for( size_t i = 1; i <= size_t(nLevels); i++ )
|
|
sum += L[i];
|
|
return sum;
|
|
}
|
|
// Compute the length of coefficients for 2D transforms
|
|
vtkm::Id ComputeCoeffLength2( std::vector<vtkm::Id> &L,
|
|
vtkm::Id nLevels )
|
|
{
|
|
vtkm::Id sum = (L[0] * L[1]); // 1st level cA
|
|
for( size_t i = 1; i <= size_t(nLevels); i++ )
|
|
{
|
|
sum += L[ i*6 - 4 ] * L[ i*6 - 3 ]; // cDh
|
|
sum += L[ i*6 - 2 ] * L[ i*6 - 1 ]; // cDv
|
|
sum += L[ i*6 - 0 ] * L[ i*6 + 1 ]; // cDd
|
|
}
|
|
return sum;
|
|
}
|
|
|
|
// Compute approximate coefficient length at a specific level
|
|
vtkm::Id 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;
|
|
}
|
|
|
|
|
|
}; // class WaveletCompressor
|
|
|
|
} // namespace worklet
|
|
} // namespace vtkm
|
|
|
|
#endif
|