2016-07-27 20:28:38 +00:00
|
|
|
//============================================================================
|
|
|
|
// 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
|
2016-08-11 16:39:20 +00:00
|
|
|
WaveletCompressor( wavelets::WaveletName name ) : WaveletDWT( name ) {}
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-08-12 00:49:48 +00:00
|
|
|
|
2016-10-10 04:52:02 +00:00
|
|
|
|
2016-08-02 20:09:29 +00:00
|
|
|
// Multi-level 1D wavelet decomposition
|
2016-08-10 20:33:14 +00:00
|
|
|
template< typename SignalArrayType, typename CoeffArrayType, typename DeviceTag >
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_CONT
|
2016-08-03 15:34:04 +00:00
|
|
|
vtkm::Id WaveDecompose( const SignalArrayType &sigIn, // Input
|
2016-08-02 20:09:29 +00:00
|
|
|
vtkm::Id nLevels, // n levels of DWT
|
2016-08-03 15:34:04 +00:00
|
|
|
CoeffArrayType &coeffOut,
|
2016-08-10 20:33:14 +00:00
|
|
|
std::vector<vtkm::Id> &L,
|
|
|
|
DeviceTag )
|
2016-08-02 20:09:29 +00:00
|
|
|
{
|
|
|
|
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
|
2016-08-08 23:59:15 +00:00
|
|
|
if( nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel( sigInLen ) )
|
2016-08-02 20:09:29 +00:00
|
|
|
{
|
2016-08-05 23:08:25 +00:00
|
|
|
throw vtkm::cont::ErrorControlBadValue("Number of levels of transform is not supported! ");
|
2016-08-02 20:09:29 +00:00
|
|
|
}
|
|
|
|
if( nLevels == 0 ) // 0 levels means no transform
|
|
|
|
{
|
2016-08-10 20:33:14 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
|
2016-08-02 20:09:29 +00:00
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
2016-08-12 00:49:48 +00:00
|
|
|
this->ComputeL( sigInLen, nLevels, L ); // memory for L is allocated by ComputeL().
|
2016-08-02 20:09:29 +00:00
|
|
|
vtkm::Id CLength = this->ComputeCoeffLength( L, nLevels );
|
2016-08-03 15:34:04 +00:00
|
|
|
VTKM_ASSERT( CLength == sigInLen );
|
2016-08-02 20:09:29 +00:00
|
|
|
|
|
|
|
vtkm::Id sigInPtr = 0; // pseudo pointer for the beginning of input array
|
2016-08-03 15:34:04 +00:00
|
|
|
vtkm::Id len = sigInLen;
|
2016-08-02 20:09:29 +00:00
|
|
|
vtkm::Id cALen = WaveletBase::GetApproxLength( len );
|
|
|
|
vtkm::Id cptr; // pseudo pointer for the beginning of output array
|
|
|
|
vtkm::Id tlen = 0;
|
2016-08-08 23:59:15 +00:00
|
|
|
std::vector<vtkm::Id> L1d(3, 0);
|
2016-08-02 20:09:29 +00:00
|
|
|
|
|
|
|
// 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;
|
|
|
|
|
2016-08-10 20:33:14 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
|
2016-08-02 20:09:29 +00:00
|
|
|
|
|
|
|
for( vtkm::Id i = nLevels; i > 0; i-- )
|
|
|
|
{
|
2016-08-08 23:59:15 +00:00
|
|
|
tlen += L[ size_t(i) ];
|
2016-08-02 20:09:29 +00:00
|
|
|
cptr = 0 + CLength - tlen - cALen;
|
|
|
|
|
|
|
|
// make input array (permutation array)
|
2016-08-03 15:34:04 +00:00
|
|
|
IdArrayType inputIndices( sigInPtr, 1, len );
|
|
|
|
PermutArrayType input( inputIndices, coeffOut );
|
2016-08-02 20:09:29 +00:00
|
|
|
// make output array
|
2016-08-03 15:34:04 +00:00
|
|
|
InterArrayType output;
|
2016-08-02 20:09:29 +00:00
|
|
|
|
2016-08-22 20:49:13 +00:00
|
|
|
WaveletDWT::DWT1D( input, output, L1d, DeviceTag() );
|
2016-08-02 20:09:29 +00:00
|
|
|
|
2016-08-03 15:34:04 +00:00
|
|
|
// move intermediate results to final array
|
2016-08-22 20:49:13 +00:00
|
|
|
WaveletBase::DeviceCopyStartX( output, coeffOut, cptr, DeviceTag() );
|
2016-08-02 20:09:29 +00:00
|
|
|
|
|
|
|
// update pseudo pointers
|
|
|
|
len = cALen;
|
|
|
|
cALen = WaveletBase::GetApproxLength( cALen );
|
|
|
|
sigInPtr = cptr;
|
|
|
|
}
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
2016-08-12 00:49:48 +00:00
|
|
|
|
2016-10-10 04:52:02 +00:00
|
|
|
|
2016-07-27 20:28:38 +00:00
|
|
|
// Multi-level 1D wavelet reconstruction
|
2016-08-10 20:33:14 +00:00
|
|
|
template< typename CoeffArrayType, typename SignalArrayType, typename DeviceTag >
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_CONT
|
2016-07-27 20:28:38 +00:00
|
|
|
vtkm::Id WaveReconstruct( const CoeffArrayType &coeffIn, // Input
|
|
|
|
vtkm::Id nLevels, // n levels of DWT
|
2016-08-08 23:59:15 +00:00
|
|
|
std::vector<vtkm::Id> &L,
|
2016-08-10 20:33:14 +00:00
|
|
|
SignalArrayType &sigOut,
|
|
|
|
DeviceTag )
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
|
|
|
VTKM_ASSERT( nLevels > 0 );
|
|
|
|
vtkm::Id LLength = nLevels + 2;
|
2016-08-08 23:59:15 +00:00
|
|
|
VTKM_ASSERT( vtkm::Id(L.size()) == LLength );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-08-08 23:59:15 +00:00
|
|
|
std::vector<vtkm::Id> L1d(3, 0); // three elements
|
|
|
|
L1d[0] = L[0];
|
|
|
|
L1d[1] = L[1];
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-08-02 20:35:00 +00:00
|
|
|
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;
|
|
|
|
|
2016-08-10 20:33:14 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( coeffIn, sigOut );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
for( vtkm::Id i = 1; i <= nLevels; i++ )
|
|
|
|
{
|
2016-08-08 23:59:15 +00:00
|
|
|
L1d[2] = this->GetApproxLengthLevN( L[ size_t(LLength-1) ], nLevels-i );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
// Make an input array
|
|
|
|
IdArrayType inputIndices( 0, 1, L1d[2] );
|
2016-08-02 20:35:00 +00:00
|
|
|
PermutArrayType input( inputIndices, sigOut );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
// Make an output array
|
2016-08-02 20:35:00 +00:00
|
|
|
OutArrayBasic output;
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-10-10 04:52:02 +00:00
|
|
|
WaveletDWT::IDWT1D( input, L1d, output, DeviceTag() );
|
2016-08-06 07:08:29 +00:00
|
|
|
VTKM_ASSERT( output.GetNumberOfValues() == L1d[2] );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
// Move output to intermediate array
|
2016-08-22 20:49:13 +00:00
|
|
|
WaveletBase::DeviceCopyStartX( output, sigOut, 0, DeviceTag() );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
L1d[0] = L1d[2];
|
2016-08-08 23:59:15 +00:00
|
|
|
L1d[1] = L[ size_t(i+1) ];
|
2016-07-27 20:28:38 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
2016-08-12 00:49:48 +00:00
|
|
|
|
2016-10-10 04:52:02 +00:00
|
|
|
|
2016-08-12 00:49:48 +00:00
|
|
|
// Multi-level 2D wavelet decomposition
|
2016-08-22 05:10:05 +00:00
|
|
|
template< typename InArrayType, typename OutArrayType, typename DeviceTag>
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_CONT
|
2016-08-23 23:52:40 +00:00
|
|
|
FLOAT_64 WaveDecompose2D( const InArrayType &sigIn, // Input
|
2016-08-12 00:49:48 +00:00
|
|
|
vtkm::Id nLevels, // n levels of DWT
|
|
|
|
vtkm::Id inX, // Input X dim
|
|
|
|
vtkm::Id inY, // Input Y dim
|
|
|
|
OutArrayType &coeffOut,
|
2016-08-22 05:10:05 +00:00
|
|
|
std::vector<vtkm::Id> &L,
|
|
|
|
DeviceTag )
|
2016-08-12 00:49:48 +00:00
|
|
|
{
|
|
|
|
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
|
|
|
|
VTKM_ASSERT( inX * inY == sigInLen );
|
|
|
|
if( nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel( inX ) ||
|
|
|
|
nLevels > WaveletBase::GetWaveletMaxLevel( inY ) )
|
|
|
|
{
|
|
|
|
throw vtkm::cont::ErrorControlBadValue("Number of levels of transform is not supported! ");
|
|
|
|
}
|
|
|
|
if( nLevels == 0 ) // 0 levels means no transform
|
|
|
|
{
|
2016-08-22 05:10:05 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
|
2016-08-12 00:49:48 +00:00
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
|
|
|
this->ComputeL2( inX, inY, nLevels, L );
|
|
|
|
vtkm::Id CLength = this->ComputeCoeffLength2( L, nLevels );
|
|
|
|
VTKM_ASSERT( CLength == sigInLen );
|
|
|
|
|
2016-08-14 22:29:27 +00:00
|
|
|
vtkm::Id currentLenX = inX;
|
|
|
|
vtkm::Id currentLenY = inY;
|
2016-08-12 00:49:48 +00:00
|
|
|
std::vector<vtkm::Id> L2d(10, 0);
|
2016-09-22 03:25:56 +00:00
|
|
|
vtkm::Float64 computationTime = 0.0;
|
2016-08-14 22:29:27 +00:00
|
|
|
|
|
|
|
typedef typename OutArrayType::ValueType OutValueType;
|
|
|
|
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
|
2016-08-12 00:49:48 +00:00
|
|
|
|
2016-10-05 22:52:07 +00:00
|
|
|
// First level transform operates writes to the output array
|
2016-10-10 04:52:02 +00:00
|
|
|
computationTime += WaveletDWT::DWT2D ( sigIn,
|
2016-10-05 22:52:07 +00:00
|
|
|
currentLenX, currentLenY,
|
|
|
|
0, 0,
|
|
|
|
currentLenX, currentLenY,
|
|
|
|
coeffOut, L2d, DeviceTag() );
|
2016-09-22 03:25:56 +00:00
|
|
|
VTKM_ASSERT( coeffOut.GetNumberOfValues() == currentLenX * currentLenY );
|
|
|
|
currentLenX = WaveletBase::GetApproxLength( currentLenX );
|
|
|
|
currentLenY = WaveletBase::GetApproxLength( currentLenY );
|
|
|
|
|
2016-10-05 22:52:07 +00:00
|
|
|
// Successor transforms writes to a temporary array
|
2016-09-22 03:25:56 +00:00
|
|
|
for( vtkm::Id i = nLevels-1; i > 0; i-- )
|
2016-08-12 00:49:48 +00:00
|
|
|
{
|
2016-08-14 22:29:27 +00:00
|
|
|
OutBasicArray tempOutput;
|
|
|
|
|
2016-08-23 23:52:40 +00:00
|
|
|
computationTime +=
|
2016-10-10 04:52:02 +00:00
|
|
|
WaveletDWT::DWT2D ( coeffOut,
|
2016-10-05 22:52:07 +00:00
|
|
|
inX, inY,
|
|
|
|
0, 0,
|
|
|
|
currentLenX, currentLenY,
|
|
|
|
tempOutput, L2d, DeviceTag() );
|
2016-08-14 22:29:27 +00:00
|
|
|
|
|
|
|
// copy results to coeffOut
|
|
|
|
WaveletBase::DeviceRectangleCopyTo( tempOutput, currentLenX, currentLenY,
|
2016-08-22 20:49:13 +00:00
|
|
|
coeffOut, inX, inY, 0, 0, DeviceTag() );
|
2016-08-14 22:29:27 +00:00
|
|
|
|
|
|
|
// update currentLen
|
|
|
|
currentLenX = WaveletBase::GetApproxLength( currentLenX );
|
|
|
|
currentLenY = WaveletBase::GetApproxLength( currentLenY );
|
|
|
|
}
|
2016-08-23 23:52:40 +00:00
|
|
|
|
|
|
|
return computationTime;
|
2016-08-14 22:29:27 +00:00
|
|
|
}
|
2016-08-12 00:49:48 +00:00
|
|
|
|
2016-08-14 22:29:27 +00:00
|
|
|
|
2016-10-10 04:52:02 +00:00
|
|
|
|
2016-08-14 22:29:27 +00:00
|
|
|
// Multi-level 2D wavelet reconstruction
|
2016-08-22 05:10:05 +00:00
|
|
|
template< typename InArrayType, typename OutArrayType, typename DeviceTag>
|
2016-10-19 22:42:58 +00:00
|
|
|
VTKM_CONT
|
2016-08-23 23:52:40 +00:00
|
|
|
FLOAT_64 WaveReconstruct2D( const InArrayType &arrIn, // Input
|
2016-08-14 22:29:27 +00:00
|
|
|
vtkm::Id nLevels, // n levels of DWT
|
|
|
|
vtkm::Id inX, // Input X dim
|
|
|
|
vtkm::Id inY, // Input Y dim
|
|
|
|
OutArrayType &arrOut,
|
2016-08-22 05:10:05 +00:00
|
|
|
std::vector<vtkm::Id> &L,
|
|
|
|
DeviceTag )
|
2016-08-14 22:29:27 +00:00
|
|
|
{
|
|
|
|
vtkm::Id arrInLen = arrIn.GetNumberOfValues();
|
|
|
|
VTKM_ASSERT( inX * inY == arrInLen );
|
|
|
|
if( nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel( inX ) ||
|
|
|
|
nLevels > WaveletBase::GetWaveletMaxLevel( inY ) )
|
|
|
|
{
|
|
|
|
throw vtkm::cont::ErrorControlBadValue("Number of levels of transform is not supported! ");
|
2016-08-12 00:49:48 +00:00
|
|
|
}
|
2016-09-22 03:25:56 +00:00
|
|
|
typedef typename OutArrayType::ValueType OutValueType;
|
|
|
|
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
|
|
|
|
vtkm::Float64 computationTime = 0.0;
|
|
|
|
|
2016-10-10 00:07:49 +00:00
|
|
|
OutBasicArray outBuffer;
|
2016-08-14 22:29:27 +00:00
|
|
|
if( nLevels == 0 ) // 0 levels means no transform
|
2016-10-10 00:07:49 +00:00
|
|
|
{
|
2016-09-22 03:25:56 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, arrOut );
|
2016-10-10 00:07:49 +00:00
|
|
|
return 0;
|
2016-08-14 22:29:27 +00:00
|
|
|
}
|
2016-09-22 03:25:56 +00:00
|
|
|
else
|
|
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, outBuffer );
|
|
|
|
|
2016-08-14 22:29:27 +00:00
|
|
|
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];
|
2016-08-12 00:49:48 +00:00
|
|
|
|
2016-09-22 03:25:56 +00:00
|
|
|
// All transforms but the last operate on temporary arrays
|
|
|
|
for( size_t i = 1; i < static_cast<size_t>(nLevels); i++ )
|
2016-08-14 22:29:27 +00:00
|
|
|
{
|
|
|
|
L2d[8] = L2d[0] + L2d[4]; // This is always true for Biorthogonal wavelets
|
|
|
|
L2d[9] = L2d[1] + L2d[3]; // (same above)
|
|
|
|
|
2016-10-10 00:07:49 +00:00
|
|
|
OutBasicArray tempOutput;
|
2016-08-14 22:29:27 +00:00
|
|
|
|
|
|
|
// IDWT
|
2016-08-23 23:52:40 +00:00
|
|
|
computationTime +=
|
2016-10-10 04:52:02 +00:00
|
|
|
WaveletDWT::IDWT2D ( outBuffer, inX, inY, 0, 0, L2d, tempOutput, DeviceTag() );
|
2016-08-14 22:29:27 +00:00
|
|
|
|
|
|
|
// copy back reconstructed block
|
|
|
|
WaveletBase::DeviceRectangleCopyTo( tempOutput, L2d[8], L2d[9],
|
2016-09-22 03:25:56 +00:00
|
|
|
outBuffer, inX, inY, 0, 0, DeviceTag() );
|
2016-08-14 22:29:27 +00:00
|
|
|
|
|
|
|
// 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];
|
2016-10-10 00:07:49 +00:00
|
|
|
|
2016-08-14 22:29:27 +00:00
|
|
|
}
|
|
|
|
|
2016-10-10 04:52:02 +00:00
|
|
|
// The last transform outputs to the final output
|
2016-09-22 03:25:56 +00:00
|
|
|
L2d[8] = L2d[0] + L2d[4];
|
|
|
|
L2d[9] = L2d[1] + L2d[3];
|
2016-10-10 00:07:49 +00:00
|
|
|
computationTime +=
|
2016-10-10 04:52:02 +00:00
|
|
|
WaveletDWT::IDWT2D ( outBuffer, inX, inY, 0, 0, L2d, arrOut, DeviceTag() );
|
2016-09-22 03:25:56 +00:00
|
|
|
|
2016-08-23 23:52:40 +00:00
|
|
|
return computationTime;
|
2016-08-12 00:49:48 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
|
2016-10-10 04:52:02 +00:00
|
|
|
|
2016-07-27 20:28:38 +00:00
|
|
|
// Squash coefficients smaller than a threshold
|
2016-08-10 20:33:14 +00:00
|
|
|
template< typename CoeffArrayType, typename DeviceTag >
|
2016-07-28 22:40:37 +00:00
|
|
|
vtkm::Id SquashCoefficients( CoeffArrayType &coeffIn,
|
2016-08-10 20:33:14 +00:00
|
|
|
vtkm::Float64 ratio,
|
2016-08-22 05:10:05 +00:00
|
|
|
DeviceTag )
|
2016-07-28 22:40:37 +00:00
|
|
|
{
|
2016-07-29 01:12:57 +00:00
|
|
|
if( ratio > 1 )
|
|
|
|
{
|
|
|
|
vtkm::Id coeffLen = coeffIn.GetNumberOfValues();
|
|
|
|
typedef typename CoeffArrayType::ValueType ValueType;
|
|
|
|
typedef vtkm::cont::ArrayHandle< ValueType > CoeffArrayBasic;
|
|
|
|
CoeffArrayBasic sortedArray;
|
2016-08-10 20:33:14 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( coeffIn, sortedArray );
|
|
|
|
WaveletBase::DeviceSort( sortedArray, DeviceTag() );
|
2016-07-29 01:12:57 +00:00
|
|
|
|
2016-08-06 02:45:30 +00:00
|
|
|
vtkm::Id n = coeffLen -
|
|
|
|
static_cast<vtkm::Id>( static_cast<vtkm::Float64>(coeffLen)/ratio );
|
2016-08-23 15:56:03 +00:00
|
|
|
vtkm::Float64 nthVal = static_cast<vtkm::Float64>
|
|
|
|
(sortedArray.GetPortalConstControl().Get(n));
|
|
|
|
if( nthVal < 0.0 )
|
|
|
|
nthVal *= -1.0;
|
2016-08-03 22:23:34 +00:00
|
|
|
typedef vtkm::worklet::wavelets::ThresholdWorklet ThresholdType;
|
2016-08-23 15:56:03 +00:00
|
|
|
ThresholdType tw( nthVal );
|
|
|
|
vtkm::worklet::DispatcherMapField< ThresholdType, DeviceTag > dispatcher( tw );
|
|
|
|
dispatcher.Invoke( coeffIn );
|
2016-07-29 01:12:57 +00:00
|
|
|
}
|
2016-07-28 22:40:37 +00:00
|
|
|
|
|
|
|
return 0;
|
|
|
|
}
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-07-29 03:07:05 +00:00
|
|
|
|
2016-10-10 04:52:02 +00:00
|
|
|
|
2016-07-29 03:07:05 +00:00
|
|
|
// Report statistics on reconstructed array
|
2016-08-10 20:33:14 +00:00
|
|
|
template< typename ArrayType, typename DeviceTag >
|
2016-07-29 03:07:05 +00:00
|
|
|
vtkm::Id EvaluateReconstruction( const ArrayType &original,
|
2016-08-10 20:33:14 +00:00
|
|
|
const ArrayType &reconstruct,
|
|
|
|
DeviceTag )
|
2016-07-29 03:07:05 +00:00
|
|
|
{
|
|
|
|
#define VAL vtkm::Float64
|
|
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
2016-08-10 20:33:14 +00:00
|
|
|
VAL VarOrig = WaveletBase::DeviceCalculateVariance( original, DeviceTag() );
|
2016-07-29 03:07:05 +00:00
|
|
|
|
|
|
|
typedef typename ArrayType::ValueType ValueType;
|
|
|
|
typedef vtkm::cont::ArrayHandle< ValueType > ArrayBasic;
|
2016-07-29 18:27:37 +00:00
|
|
|
ArrayBasic errorArray, errorSquare;
|
2016-07-29 03:07:05 +00:00
|
|
|
|
2016-07-29 18:27:37 +00:00
|
|
|
// Use a worklet to calculate point-wise error, and its square
|
2016-07-29 03:07:05 +00:00
|
|
|
typedef vtkm::worklet::wavelets::Differencer DifferencerWorklet;
|
|
|
|
DifferencerWorklet dw;
|
2016-07-29 18:27:37 +00:00
|
|
|
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 );
|
2016-07-29 03:07:05 +00:00
|
|
|
|
2016-08-10 20:33:14 +00:00
|
|
|
VAL varErr = WaveletBase::DeviceCalculateVariance( errorArray, DeviceTag() );
|
2016-07-29 21:16:10 +00:00
|
|
|
VAL snr, decibels;
|
|
|
|
if( varErr != 0.0 )
|
|
|
|
{
|
|
|
|
snr = VarOrig / varErr;
|
|
|
|
decibels = 10 * vtkm::Log10( snr );
|
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
|
|
|
snr = vtkm::Infinity64();
|
|
|
|
decibels = vtkm::Infinity64();
|
|
|
|
}
|
2016-07-29 03:07:05 +00:00
|
|
|
|
2016-08-10 20:33:14 +00:00
|
|
|
VAL origMax = WaveletBase::DeviceMax( original, DeviceTag() );
|
|
|
|
VAL origMin = WaveletBase::DeviceMin( original, DeviceTag() );
|
|
|
|
VAL errorMax = WaveletBase::DeviceMaxAbs( errorArray, DeviceTag() );
|
2016-07-29 03:07:05 +00:00
|
|
|
VAL range = origMax - origMin;
|
|
|
|
|
2016-08-10 20:33:14 +00:00
|
|
|
VAL squareSum = WaveletBase::DeviceSum( errorSquare, DeviceTag() );
|
2016-07-29 04:33:47 +00:00
|
|
|
VAL rmse = vtkm::Sqrt( MAKEVAL(squareSum) / MAKEVAL(errorArray.GetNumberOfValues()) );
|
|
|
|
|
2016-07-29 18:27:37 +00:00
|
|
|
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;
|
2016-07-29 03:07:05 +00:00
|
|
|
#undef MAKEVAL
|
|
|
|
#undef VAL
|
|
|
|
|
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
2016-10-10 04:52:02 +00:00
|
|
|
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
// Compute the book keeping array L for 1D wavelet decomposition
|
2016-08-08 23:59:15 +00:00
|
|
|
void ComputeL( vtkm::Id sigInLen,
|
2016-08-12 00:49:48 +00:00
|
|
|
vtkm::Id nLev,
|
2016-08-08 23:59:15 +00:00
|
|
|
std::vector<vtkm::Id> &L )
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2016-08-12 00:49:48 +00:00
|
|
|
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-- )
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
|
|
|
L[i-1] = WaveletBase::GetApproxLength( L[i] );
|
|
|
|
L[i] = WaveletBase::GetDetailLength( L[i] );
|
|
|
|
}
|
|
|
|
}
|
2016-10-10 04:52:02 +00:00
|
|
|
|
|
|
|
|
|
|
|
|
2016-08-12 00:49:48 +00:00
|
|
|
// Compute the book keeping array L for 2D wavelet decomposition
|
|
|
|
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 + 0 ] = 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 ]);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-10-10 04:52:02 +00:00
|
|
|
|
2016-08-12 00:49:48 +00:00
|
|
|
// Compute the length of coefficients for 1D transforms
|
2016-08-08 23:59:15 +00:00
|
|
|
vtkm::Id ComputeCoeffLength( std::vector<vtkm::Id> &L,
|
|
|
|
vtkm::Id nLevels )
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2016-08-08 23:59:15 +00:00
|
|
|
vtkm::Id sum = L[0]; // 1st level cA
|
|
|
|
for( size_t i = 1; i <= size_t(nLevels); i++ )
|
2016-07-27 20:28:38 +00:00
|
|
|
sum += L[i];
|
|
|
|
return sum;
|
|
|
|
}
|
2016-08-12 00:49:48 +00:00
|
|
|
// 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;
|
|
|
|
}
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
// 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
|