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.
|
|
|
|
//
|
2017-09-20 21:33:44 +00:00
|
|
|
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
|
2016-07-27 20:28:38 +00:00
|
|
|
// Copyright 2014 UT-Battelle, LLC.
|
|
|
|
// Copyright 2014 Los Alamos National Security.
|
|
|
|
//
|
2017-09-20 21:33:44 +00:00
|
|
|
// Under the terms of Contract DE-NA0003525 with NTESS,
|
2016-07-27 20:28:38 +00:00
|
|
|
// 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>
|
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
namespace vtkm
|
|
|
|
{
|
|
|
|
namespace worklet
|
|
|
|
{
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
class WaveletCompressor : public vtkm::worklet::wavelets::WaveletDWT
|
|
|
|
{
|
|
|
|
public:
|
|
|
|
// Constructor
|
2017-05-18 14:29:41 +00:00
|
|
|
WaveletCompressor(wavelets::WaveletName name)
|
|
|
|
: WaveletDWT(name)
|
|
|
|
{
|
|
|
|
}
|
2016-08-12 00:49:48 +00:00
|
|
|
|
2016-08-02 20:09:29 +00:00
|
|
|
// Multi-level 1D wavelet decomposition
|
2017-05-18 14:29:41 +00:00
|
|
|
template <typename SignalArrayType, typename CoeffArrayType, typename DeviceTag>
|
|
|
|
VTKM_CONT vtkm::Id WaveDecompose(const SignalArrayType& sigIn, // Input
|
|
|
|
vtkm::Id nLevels, // n levels of DWT
|
2017-05-26 17:53:28 +00:00
|
|
|
CoeffArrayType& coeffOut,
|
|
|
|
std::vector<vtkm::Id>& L,
|
|
|
|
DeviceTag)
|
2016-08-02 20:09:29 +00:00
|
|
|
{
|
|
|
|
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
|
2017-05-18 14:29:41 +00:00
|
|
|
if (nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel(sigInLen))
|
2016-08-02 20:09:29 +00:00
|
|
|
{
|
2017-01-09 21:15:32 +00:00
|
|
|
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
|
2016-08-02 20:09:29 +00:00
|
|
|
}
|
2017-05-18 14:29:41 +00:00
|
|
|
if (nLevels == 0) // 0 levels means no transform
|
2016-08-02 20:09:29 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm<DeviceTag>::Copy(sigIn, coeffOut);
|
2016-08-02 20:09:29 +00:00
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
this->ComputeL(sigInLen, nLevels, L); // memory for L is allocated by ComputeL().
|
|
|
|
vtkm::Id CLength = this->ComputeCoeffLength(L, nLevels);
|
|
|
|
VTKM_ASSERT(CLength == sigInLen);
|
2016-08-02 20:09:29 +00:00
|
|
|
|
2017-05-18 14:29:41 +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;
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::Id cALen = WaveletBase::GetApproxLength(len);
|
|
|
|
vtkm::Id cptr; // pseudo pointer for the beginning of output array
|
2016-08-02 20:09:29 +00:00
|
|
|
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
|
2017-05-18 14:29:41 +00:00
|
|
|
typedef typename CoeffArrayType::ValueType OutputValueType;
|
|
|
|
typedef vtkm::cont::ArrayHandle<OutputValueType> InterArrayType;
|
2016-08-02 20:09:29 +00:00
|
|
|
|
|
|
|
// Define a few more types
|
2017-05-18 14:29:41 +00:00
|
|
|
typedef vtkm::cont::ArrayHandleCounting<vtkm::Id> IdArrayType;
|
|
|
|
typedef vtkm::cont::ArrayHandlePermutation<IdArrayType, CoeffArrayType> PermutArrayType;
|
2016-08-02 20:09:29 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm<DeviceTag>::Copy(sigIn, coeffOut);
|
2016-08-02 20:09:29 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
for (vtkm::Id i = nLevels; i > 0; i--)
|
2016-08-02 20:09:29 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
tlen += L[size_t(i)];
|
2016-08-02 20:09:29 +00:00
|
|
|
cptr = 0 + CLength - tlen - cALen;
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2016-08-02 20:09:29 +00:00
|
|
|
// make input array (permutation array)
|
2017-05-18 14:29:41 +00:00
|
|
|
IdArrayType inputIndices(sigInPtr, 1, len);
|
|
|
|
PermutArrayType input(inputIndices, coeffOut);
|
2017-05-22 20:47:13 +00:00
|
|
|
// make output array
|
2017-05-18 14:29:41 +00:00
|
|
|
InterArrayType output;
|
2016-08-02 20:09:29 +00:00
|
|
|
|
2017-05-18 14:29:41 +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
|
2017-05-18 14:29:41 +00:00
|
|
|
WaveletBase::DeviceCopyStartX(output, coeffOut, cptr, DeviceTag());
|
2016-08-02 20:09:29 +00:00
|
|
|
|
|
|
|
// update pseudo pointers
|
|
|
|
len = cALen;
|
2017-05-18 14:29:41 +00:00
|
|
|
cALen = WaveletBase::GetApproxLength(cALen);
|
2016-08-02 20:09:29 +00:00
|
|
|
sigInPtr = cptr;
|
|
|
|
}
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
|
|
|
// Multi-level 1D wavelet reconstruction
|
2017-05-18 14:29:41 +00:00
|
|
|
template <typename CoeffArrayType, typename SignalArrayType, typename DeviceTag>
|
|
|
|
VTKM_CONT vtkm::Id WaveReconstruct(const CoeffArrayType& coeffIn, // Input
|
|
|
|
vtkm::Id nLevels, // n levels of DWT
|
2017-05-26 17:53:28 +00:00
|
|
|
std::vector<vtkm::Id>& L,
|
|
|
|
SignalArrayType& sigOut,
|
|
|
|
DeviceTag)
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
VTKM_ASSERT(nLevels > 0);
|
2016-07-27 20:28:38 +00:00
|
|
|
vtkm::Id LLength = nLevels + 2;
|
2017-05-18 14:29:41 +00:00
|
|
|
VTKM_ASSERT(vtkm::Id(L.size()) == LLength);
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
std::vector<vtkm::Id> L1d(3, 0); // three elements
|
2016-08-08 23:59:15 +00:00
|
|
|
L1d[0] = L[0];
|
|
|
|
L1d[1] = L[1];
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2017-05-18 14:29:41 +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-02 20:35:00 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm<DeviceTag>::Copy(coeffIn, sigOut);
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
for (vtkm::Id i = 1; i <= nLevels; i++)
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +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
|
2017-05-18 14:29:41 +00:00
|
|
|
IdArrayType inputIndices(0, 1, L1d[2]);
|
|
|
|
PermutArrayType input(inputIndices, sigOut);
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2016-07-27 20:28:38 +00:00
|
|
|
// Make an output array
|
2016-08-02 20:35:00 +00:00
|
|
|
OutArrayBasic output;
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
WaveletDWT::IDWT1D(input, L1d, output, DeviceTag());
|
|
|
|
VTKM_ASSERT(output.GetNumberOfValues() == L1d[2]);
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
// Move output to intermediate array
|
2017-05-18 14:29:41 +00:00
|
|
|
WaveletBase::DeviceCopyStartX(output, sigOut, 0, DeviceTag());
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
L1d[0] = L1d[2];
|
2017-05-18 14:29:41 +00:00
|
|
|
L1d[1] = L[size_t(i + 1)];
|
2016-07-27 20:28:38 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
2017-02-21 23:56:57 +00:00
|
|
|
// Multi-level 3D wavelet decomposition
|
2017-05-18 14:29:41 +00:00
|
|
|
template <typename InArrayType, typename OutArrayType, typename DeviceTag>
|
|
|
|
VTKM_CONT vtkm::Float64 WaveDecompose3D(InArrayType& sigIn, // Input
|
|
|
|
vtkm::Id nLevels, // n levels of DWT
|
2017-05-26 17:53:28 +00:00
|
|
|
vtkm::Id inX,
|
|
|
|
vtkm::Id inY,
|
|
|
|
vtkm::Id inZ,
|
2017-05-18 14:29:41 +00:00
|
|
|
OutArrayType& coeffOut,
|
|
|
|
bool discardSigIn, // can we discard sigIn on devices?
|
|
|
|
DeviceTag)
|
2017-02-21 23:56:57 +00:00
|
|
|
{
|
|
|
|
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
|
2017-05-18 14:29:41 +00:00
|
|
|
VTKM_ASSERT(inX * inY * inZ == sigInLen);
|
|
|
|
if (nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel(inX) ||
|
|
|
|
nLevels > WaveletBase::GetWaveletMaxLevel(inY) ||
|
|
|
|
nLevels > WaveletBase::GetWaveletMaxLevel(inZ))
|
2017-02-21 23:56:57 +00:00
|
|
|
{
|
2017-04-14 20:14:09 +00:00
|
|
|
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
|
2017-02-21 23:56:57 +00:00
|
|
|
}
|
2017-05-18 14:29:41 +00:00
|
|
|
if (nLevels == 0) // 0 levels means no transform
|
2017-02-21 23:56:57 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm<DeviceTag>::Copy(sigIn, coeffOut);
|
2017-02-21 23:56:57 +00:00
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::Id currentLenX = inX;
|
|
|
|
vtkm::Id currentLenY = inY;
|
|
|
|
vtkm::Id currentLenZ = inZ;
|
2017-02-21 23:56:57 +00:00
|
|
|
std::vector<vtkm::Id> L3d(27, 0);
|
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
typedef typename OutArrayType::ValueType OutValueType;
|
|
|
|
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
|
2017-02-21 23:56:57 +00:00
|
|
|
|
2017-02-24 18:56:24 +00:00
|
|
|
// First level transform writes to the output array
|
2017-05-26 17:53:28 +00:00
|
|
|
vtkm::Float64 computationTime = WaveletDWT::DWT3D(sigIn,
|
|
|
|
inX,
|
|
|
|
inY,
|
|
|
|
inZ,
|
|
|
|
0,
|
|
|
|
0,
|
|
|
|
0,
|
|
|
|
currentLenX,
|
|
|
|
currentLenY,
|
|
|
|
currentLenZ,
|
|
|
|
coeffOut,
|
|
|
|
discardSigIn,
|
|
|
|
DeviceTag());
|
2017-02-24 18:56:24 +00:00
|
|
|
|
2017-02-21 23:56:57 +00:00
|
|
|
// Successor transforms writes to a temporary array
|
2017-05-18 14:29:41 +00:00
|
|
|
for (vtkm::Id i = nLevels - 1; i > 0; i--)
|
2017-02-21 23:56:57 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
currentLenX = WaveletBase::GetApproxLength(currentLenX);
|
|
|
|
currentLenY = WaveletBase::GetApproxLength(currentLenY);
|
|
|
|
currentLenZ = WaveletBase::GetApproxLength(currentLenZ);
|
2017-02-25 18:55:00 +00:00
|
|
|
|
2017-02-21 23:56:57 +00:00
|
|
|
OutBasicArray tempOutput;
|
|
|
|
|
2017-05-26 17:53:28 +00:00
|
|
|
computationTime += WaveletDWT::DWT3D(coeffOut,
|
|
|
|
inX,
|
|
|
|
inY,
|
|
|
|
inZ,
|
|
|
|
0,
|
|
|
|
0,
|
|
|
|
0,
|
|
|
|
currentLenX,
|
|
|
|
currentLenY,
|
|
|
|
currentLenZ,
|
|
|
|
tempOutput,
|
|
|
|
false,
|
|
|
|
DeviceTag());
|
2017-02-21 23:56:57 +00:00
|
|
|
|
|
|
|
// copy results to coeffOut
|
2017-05-26 17:53:28 +00:00
|
|
|
WaveletBase::DeviceCubeCopyTo(tempOutput,
|
|
|
|
currentLenX,
|
|
|
|
currentLenY,
|
|
|
|
currentLenZ,
|
|
|
|
coeffOut,
|
|
|
|
inX,
|
|
|
|
inY,
|
|
|
|
inZ,
|
|
|
|
0,
|
|
|
|
0,
|
|
|
|
0,
|
|
|
|
DeviceTag());
|
2017-02-25 18:55:00 +00:00
|
|
|
}
|
2017-02-21 23:56:57 +00:00
|
|
|
|
|
|
|
return computationTime;
|
|
|
|
}
|
|
|
|
|
|
|
|
// Multi-level 3D wavelet reconstruction
|
2017-05-18 14:29:41 +00:00
|
|
|
template <typename InArrayType, typename OutArrayType, typename DeviceTag>
|
|
|
|
VTKM_CONT vtkm::Float64 WaveReconstruct3D(
|
|
|
|
InArrayType& arrIn, // Input
|
|
|
|
vtkm::Id nLevels, // n levels of DWT
|
2017-05-26 17:53:28 +00:00
|
|
|
vtkm::Id inX,
|
|
|
|
vtkm::Id inY,
|
|
|
|
vtkm::Id inZ,
|
|
|
|
OutArrayType& arrOut,
|
2017-05-18 14:29:41 +00:00
|
|
|
bool discardArrIn, // can we discard input for more memory?
|
|
|
|
DeviceTag)
|
2017-02-21 23:56:57 +00:00
|
|
|
{
|
|
|
|
vtkm::Id arrInLen = arrIn.GetNumberOfValues();
|
2017-05-18 14:29:41 +00:00
|
|
|
VTKM_ASSERT(inX * inY * inZ == arrInLen);
|
|
|
|
if (nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel(inX) ||
|
|
|
|
nLevels > WaveletBase::GetWaveletMaxLevel(inY) ||
|
|
|
|
nLevels > WaveletBase::GetWaveletMaxLevel(inZ))
|
2017-02-21 23:56:57 +00:00
|
|
|
{
|
2017-04-14 20:14:09 +00:00
|
|
|
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
|
2017-02-21 23:56:57 +00:00
|
|
|
}
|
2017-05-18 14:29:41 +00:00
|
|
|
typedef typename OutArrayType::ValueType OutValueType;
|
|
|
|
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
|
2017-02-26 02:49:12 +00:00
|
|
|
vtkm::Float64 computationTime = 0.0;
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2017-02-21 23:56:57 +00:00
|
|
|
OutBasicArray outBuffer;
|
2017-05-18 14:29:41 +00:00
|
|
|
if (nLevels == 0) // 0 levels means no transform
|
2017-05-22 20:47:13 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm<DeviceTag>::Copy(arrIn, arrOut);
|
2017-05-22 20:47:13 +00:00
|
|
|
return 0;
|
2017-02-21 23:56:57 +00:00
|
|
|
}
|
2017-05-18 14:29:41 +00:00
|
|
|
else if (discardArrIn)
|
2017-05-18 17:51:35 +00:00
|
|
|
{
|
2017-05-22 20:47:13 +00:00
|
|
|
outBuffer = arrIn;
|
2017-05-18 17:51:35 +00:00
|
|
|
}
|
2017-05-22 20:47:13 +00:00
|
|
|
else
|
2017-05-18 17:51:35 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm<DeviceTag>::Copy(arrIn, outBuffer);
|
2017-05-18 17:51:35 +00:00
|
|
|
}
|
2017-02-21 23:56:57 +00:00
|
|
|
|
2017-02-25 18:55:00 +00:00
|
|
|
std::vector<vtkm::Id> L;
|
2017-05-18 14:29:41 +00:00
|
|
|
this->ComputeL3(inX, inY, inZ, nLevels, L);
|
2017-02-21 23:56:57 +00:00
|
|
|
std::vector<vtkm::Id> L3d(27, 0);
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2017-02-25 18:55:00 +00:00
|
|
|
// All transforms but the last level operate on temporary arrays
|
2017-05-18 14:29:41 +00:00
|
|
|
for (size_t i = 0; i < 24; i++)
|
2017-05-18 17:51:35 +00:00
|
|
|
{
|
2017-02-25 18:55:00 +00:00
|
|
|
L3d[i] = L[i];
|
2017-05-18 17:51:35 +00:00
|
|
|
}
|
2017-05-18 14:29:41 +00:00
|
|
|
for (size_t i = 1; i < static_cast<size_t>(nLevels); i++)
|
2017-02-21 23:56:57 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
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
|
2017-02-21 23:56:57 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
OutBasicArray tempOutput;
|
2017-02-21 23:56:57 +00:00
|
|
|
|
|
|
|
// IDWT
|
2017-05-18 14:29:41 +00:00
|
|
|
computationTime +=
|
|
|
|
WaveletDWT::IDWT3D(outBuffer, inX, inY, inZ, 0, 0, 0, L3d, tempOutput, false, DeviceTag());
|
2017-02-21 23:56:57 +00:00
|
|
|
|
|
|
|
// copy back reconstructed block
|
2017-05-26 17:53:28 +00:00
|
|
|
WaveletBase::DeviceCubeCopyTo(
|
|
|
|
tempOutput, L3d[24], L3d[25], L3d[26], outBuffer, inX, inY, inZ, 0, 0, 0, DeviceTag());
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2017-02-25 18:55:00 +00:00
|
|
|
// update L3d array
|
2017-05-18 14:29:41 +00:00
|
|
|
L3d[0] = L3d[24];
|
|
|
|
L3d[1] = L3d[25];
|
|
|
|
L3d[2] = L3d[26];
|
|
|
|
for (size_t j = 3; j < 24; j++)
|
2017-05-18 17:51:35 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
L3d[j] = L[21 * i + j];
|
2017-05-18 17:51:35 +00:00
|
|
|
}
|
2017-02-21 23:56:57 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
// The last transform outputs to the final output
|
2017-02-25 18:55:00 +00:00
|
|
|
L3d[24] = L3d[0] + L3d[12];
|
|
|
|
L3d[25] = L3d[1] + L3d[7];
|
2017-05-22 20:47:13 +00:00
|
|
|
L3d[26] = L3d[2] + L3d[5];
|
2017-05-18 14:29:41 +00:00
|
|
|
computationTime +=
|
|
|
|
WaveletDWT::IDWT3D(outBuffer, inX, inY, inZ, 0, 0, 0, L3d, arrOut, true, DeviceTag());
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2017-02-26 02:49:12 +00:00
|
|
|
return computationTime;
|
2017-02-21 23:56:57 +00:00
|
|
|
}
|
|
|
|
|
2016-08-12 00:49:48 +00:00
|
|
|
// Multi-level 2D wavelet decomposition
|
2017-05-18 14:29:41 +00:00
|
|
|
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
|
2017-05-26 17:53:28 +00:00
|
|
|
OutArrayType& coeffOut,
|
|
|
|
std::vector<vtkm::Id>& L,
|
2017-05-18 14:29:41 +00:00
|
|
|
DeviceTag)
|
2016-08-12 00:49:48 +00:00
|
|
|
{
|
|
|
|
vtkm::Id sigInLen = sigIn.GetNumberOfValues();
|
2017-05-18 14:29:41 +00:00
|
|
|
VTKM_ASSERT(inX * inY == sigInLen);
|
|
|
|
if (nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel(inX) ||
|
|
|
|
nLevels > WaveletBase::GetWaveletMaxLevel(inY))
|
2016-08-12 00:49:48 +00:00
|
|
|
{
|
2017-01-09 21:15:32 +00:00
|
|
|
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
|
2016-08-12 00:49:48 +00:00
|
|
|
}
|
2017-05-18 14:29:41 +00:00
|
|
|
if (nLevels == 0) // 0 levels means no transform
|
2016-08-12 00:49:48 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm<DeviceTag>::Copy(sigIn, coeffOut);
|
2016-08-12 00:49:48 +00:00
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
this->ComputeL2(inX, inY, nLevels, L);
|
|
|
|
vtkm::Id CLength = this->ComputeCoeffLength2(L, nLevels);
|
|
|
|
VTKM_ASSERT(CLength == sigInLen);
|
2016-08-12 00:49:48 +00:00
|
|
|
|
2017-05-18 14:29:41 +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
|
|
|
|
2017-05-18 14:29:41 +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
|
2017-05-26 17:53:28 +00:00
|
|
|
computationTime += WaveletDWT::DWT2D(
|
|
|
|
sigIn, currentLenX, currentLenY, 0, 0, currentLenX, currentLenY, coeffOut, L2d, DeviceTag());
|
2017-05-18 14:29:41 +00:00
|
|
|
VTKM_ASSERT(coeffOut.GetNumberOfValues() == currentLenX * currentLenY);
|
|
|
|
currentLenX = WaveletBase::GetApproxLength(currentLenX);
|
|
|
|
currentLenY = WaveletBase::GetApproxLength(currentLenY);
|
2016-09-22 03:25:56 +00:00
|
|
|
|
2016-10-05 22:52:07 +00:00
|
|
|
// Successor transforms writes to a temporary array
|
2017-05-18 14:29:41 +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;
|
|
|
|
|
2017-05-26 17:53:28 +00:00
|
|
|
computationTime += WaveletDWT::DWT2D(
|
|
|
|
coeffOut, inX, inY, 0, 0, currentLenX, currentLenY, tempOutput, L2d, DeviceTag());
|
2016-08-14 22:29:27 +00:00
|
|
|
|
|
|
|
// copy results to coeffOut
|
2017-05-26 17:53:28 +00:00
|
|
|
WaveletBase::DeviceRectangleCopyTo(
|
|
|
|
tempOutput, currentLenX, currentLenY, coeffOut, inX, inY, 0, 0, DeviceTag());
|
2016-08-14 22:29:27 +00:00
|
|
|
|
|
|
|
// update currentLen
|
2017-05-18 14:29:41 +00:00
|
|
|
currentLenX = WaveletBase::GetApproxLength(currentLenX);
|
|
|
|
currentLenY = WaveletBase::GetApproxLength(currentLenY);
|
2016-08-14 22:29:27 +00:00
|
|
|
}
|
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
|
|
|
// Multi-level 2D wavelet reconstruction
|
2017-05-18 14:29:41 +00:00
|
|
|
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
|
2017-05-26 17:53:28 +00:00
|
|
|
OutArrayType& arrOut,
|
|
|
|
std::vector<vtkm::Id>& L,
|
2017-05-18 14:29:41 +00:00
|
|
|
DeviceTag)
|
2016-08-14 22:29:27 +00:00
|
|
|
{
|
|
|
|
vtkm::Id arrInLen = arrIn.GetNumberOfValues();
|
2017-05-18 14:29:41 +00:00
|
|
|
VTKM_ASSERT(inX * inY == arrInLen);
|
|
|
|
if (nLevels < 0 || nLevels > WaveletBase::GetWaveletMaxLevel(inX) ||
|
|
|
|
nLevels > WaveletBase::GetWaveletMaxLevel(inY))
|
2016-08-14 22:29:27 +00:00
|
|
|
{
|
2017-01-09 21:15:32 +00:00
|
|
|
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
|
2016-08-12 00:49:48 +00:00
|
|
|
}
|
2017-05-18 14:29:41 +00:00
|
|
|
typedef typename OutArrayType::ValueType OutValueType;
|
|
|
|
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
|
2016-09-22 03:25:56 +00:00
|
|
|
vtkm::Float64 computationTime = 0.0;
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2016-10-10 00:07:49 +00:00
|
|
|
OutBasicArray outBuffer;
|
2017-05-18 14:29:41 +00:00
|
|
|
if (nLevels == 0) // 0 levels means no transform
|
2017-05-22 20:47:13 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm<DeviceTag>::Copy(arrIn, arrOut);
|
2017-05-22 20:47:13 +00:00
|
|
|
return 0;
|
2016-08-14 22:29:27 +00:00
|
|
|
}
|
2016-09-22 03:25:56 +00:00
|
|
|
else
|
2017-05-18 17:51:35 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm<DeviceTag>::Copy(arrIn, outBuffer);
|
2017-05-18 17:51:35 +00:00
|
|
|
}
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
VTKM_ASSERT(vtkm::Id(L.size()) == 6 * nLevels + 4);
|
2016-08-14 22:29:27 +00:00
|
|
|
|
|
|
|
std::vector<vtkm::Id> L2d(10, 0);
|
2017-05-18 14:29:41 +00:00
|
|
|
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];
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2016-09-22 03:25:56 +00:00
|
|
|
// All transforms but the last operate on temporary arrays
|
2017-05-18 14:29:41 +00:00
|
|
|
for (size_t i = 1; i < static_cast<size_t>(nLevels); i++)
|
2016-08-14 22:29:27 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
L2d[8] = L2d[0] + L2d[4]; // This is always true for Biorthogonal wavelets
|
|
|
|
L2d[9] = L2d[1] + L2d[3]; // (same above)
|
2016-08-14 22:29:27 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
OutBasicArray tempOutput;
|
2016-08-14 22:29:27 +00:00
|
|
|
|
|
|
|
// IDWT
|
2016-08-23 23:52:40 +00:00
|
|
|
computationTime +=
|
2017-05-18 14:29:41 +00:00
|
|
|
WaveletDWT::IDWT2D(outBuffer, inX, inY, 0, 0, L2d, tempOutput, DeviceTag());
|
2016-08-14 22:29:27 +00:00
|
|
|
|
|
|
|
// copy back reconstructed block
|
2017-05-26 17:53:28 +00:00
|
|
|
WaveletBase::DeviceRectangleCopyTo(
|
|
|
|
tempOutput, L2d[8], L2d[9], outBuffer, inX, inY, 0, 0, DeviceTag());
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2016-08-14 22:29:27 +00:00
|
|
|
// update L2d array
|
2017-05-18 14:29:41 +00:00
|
|
|
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-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];
|
2017-05-18 14:29:41 +00:00
|
|
|
computationTime += WaveletDWT::IDWT2D(outBuffer, inX, inY, 0, 0, L2d, arrOut, DeviceTag());
|
2016-09-22 03:25:56 +00:00
|
|
|
|
2017-05-22 20:47:13 +00:00
|
|
|
return computationTime;
|
2016-08-12 00:49:48 +00:00
|
|
|
}
|
|
|
|
|
2016-07-27 20:28:38 +00:00
|
|
|
// Squash coefficients smaller than a threshold
|
2017-05-18 14:29:41 +00:00
|
|
|
template <typename CoeffArrayType, typename DeviceTag>
|
|
|
|
vtkm::Id SquashCoefficients(CoeffArrayType& coeffIn, vtkm::Float64 ratio, DeviceTag)
|
2016-07-28 22:40:37 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
if (ratio > 1.0)
|
2016-07-29 01:12:57 +00:00
|
|
|
{
|
|
|
|
vtkm::Id coeffLen = coeffIn.GetNumberOfValues();
|
|
|
|
typedef typename CoeffArrayType::ValueType ValueType;
|
2017-05-18 14:29:41 +00:00
|
|
|
typedef vtkm::cont::ArrayHandle<ValueType> CoeffArrayBasic;
|
2016-07-29 01:12:57 +00:00
|
|
|
CoeffArrayBasic sortedArray;
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm<DeviceTag>::Copy(coeffIn, sortedArray);
|
2017-03-09 03:29:27 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
WaveletBase::DeviceSort(sortedArray, DeviceTag());
|
2017-03-09 03:29:27 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
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)
|
2017-05-18 17:51:35 +00:00
|
|
|
{
|
2016-08-23 15:56:03 +00:00
|
|
|
nthVal *= -1.0;
|
2017-05-18 17:51:35 +00:00
|
|
|
}
|
2017-05-30 16:39:46 +00:00
|
|
|
|
2016-08-03 22:23:34 +00:00
|
|
|
typedef vtkm::worklet::wavelets::ThresholdWorklet ThresholdType;
|
2017-05-18 14:29:41 +00:00
|
|
|
ThresholdType thresholdWorklet(nthVal);
|
|
|
|
vtkm::worklet::DispatcherMapField<ThresholdType, DeviceTag> dispatcher(thresholdWorklet);
|
|
|
|
dispatcher.Invoke(coeffIn);
|
2017-05-22 20:47:13 +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
|
|
|
// Report statistics on reconstructed array
|
2017-05-18 14:29:41 +00:00
|
|
|
template <typename ArrayType, typename DeviceTag>
|
2017-05-26 17:53:28 +00:00
|
|
|
vtkm::Id EvaluateReconstruction(const ArrayType& original,
|
|
|
|
const ArrayType& reconstruct,
|
2017-05-18 14:29:41 +00:00
|
|
|
DeviceTag)
|
2016-07-29 03:07:05 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
#define VAL vtkm::Float64
|
|
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
|
|
VAL VarOrig = WaveletBase::DeviceCalculateVariance(original, DeviceTag());
|
2016-07-29 03:07:05 +00:00
|
|
|
|
|
|
|
typedef typename ArrayType::ValueType ValueType;
|
2017-05-18 14:29:41 +00:00
|
|
|
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;
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::worklet::DispatcherMapField<DifferencerWorklet> dwDispatcher(dw);
|
|
|
|
dwDispatcher.Invoke(original, reconstruct, errorArray);
|
2016-07-29 18:27:37 +00:00
|
|
|
|
|
|
|
typedef vtkm::worklet::wavelets::SquareWorklet SquareWorklet;
|
|
|
|
SquareWorklet sw;
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::worklet::DispatcherMapField<SquareWorklet> swDispatcher(sw);
|
|
|
|
swDispatcher.Invoke(errorArray, errorSquare);
|
2016-07-29 03:07:05 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
VAL varErr = WaveletBase::DeviceCalculateVariance(errorArray, DeviceTag());
|
2016-07-29 21:16:10 +00:00
|
|
|
VAL snr, decibels;
|
2017-05-18 14:29:41 +00:00
|
|
|
if (varErr != 0.0)
|
2016-07-29 21:16:10 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
snr = VarOrig / varErr;
|
|
|
|
decibels = 10 * vtkm::Log10(snr);
|
2016-07-29 21:16:10 +00:00
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
snr = vtkm::Infinity64();
|
|
|
|
decibels = vtkm::Infinity64();
|
2016-07-29 21:16:10 +00:00
|
|
|
}
|
2016-07-29 03:07:05 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
VAL origMax = WaveletBase::DeviceMax(original, DeviceTag());
|
|
|
|
VAL origMin = WaveletBase::DeviceMin(original, DeviceTag());
|
|
|
|
VAL errorMax = WaveletBase::DeviceMaxAbs(errorArray, DeviceTag());
|
|
|
|
VAL range = origMax - origMin;
|
2016-07-29 03:07:05 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
VAL squareSum = WaveletBase::DeviceSum(errorSquare, DeviceTag());
|
|
|
|
VAL rmse = vtkm::Sqrt(MAKEVAL(squareSum) / MAKEVAL(errorArray.GetNumberOfValues()));
|
2016-07-29 04:33:47 +00:00
|
|
|
|
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;
|
2017-05-22 20:47:13 +00:00
|
|
|
std::cout << "L-infy norm = " << errorMax
|
2016-07-29 18:27:37 +00:00
|
|
|
<< ", after normalization = " << errorMax / range << std::endl;
|
2017-05-18 14:29:41 +00:00
|
|
|
std::cout << "RMSE = " << rmse << ", after normalization = " << rmse / range
|
|
|
|
<< std::endl;
|
|
|
|
#undef MAKEVAL
|
|
|
|
#undef VAL
|
2016-07-29 03:07:05 +00:00
|
|
|
|
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
2017-02-21 23:56:57 +00:00
|
|
|
// Compute the book keeping array L for 1D DWT
|
2017-05-18 14:29:41 +00:00
|
|
|
void ComputeL(vtkm::Id sigInLen, vtkm::Id nLev, std::vector<vtkm::Id>& L)
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +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
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
L[i - 1] = WaveletBase::GetApproxLength(L[i]);
|
|
|
|
L[i] = WaveletBase::GetDetailLength(L[i]);
|
2016-07-27 20:28:38 +00:00
|
|
|
}
|
|
|
|
}
|
2016-10-10 04:52:02 +00:00
|
|
|
|
2017-02-21 23:56:57 +00:00
|
|
|
// Compute the book keeping array L for 2D DWT
|
2017-05-18 14:29:41 +00:00
|
|
|
void ComputeL2(vtkm::Id inX, vtkm::Id inY, vtkm::Id nLev, std::vector<vtkm::Id>& L)
|
2016-08-12 00:49:48 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
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--)
|
2016-08-12 00:49:48 +00:00
|
|
|
{
|
|
|
|
// cA
|
2017-05-18 14:29:41 +00:00
|
|
|
L[i * 6 - 6] = WaveletBase::GetApproxLength(L[i * 6 + 0]);
|
|
|
|
L[i * 6 - 5] = WaveletBase::GetApproxLength(L[i * 6 + 1]);
|
2016-08-12 00:49:48 +00:00
|
|
|
|
|
|
|
// cDh
|
2017-05-18 14:29:41 +00:00
|
|
|
L[i * 6 - 4] = WaveletBase::GetApproxLength(L[i * 6 + 0]);
|
|
|
|
L[i * 6 - 3] = WaveletBase::GetDetailLength(L[i * 6 + 1]);
|
2016-08-12 00:49:48 +00:00
|
|
|
|
|
|
|
// cDv
|
2017-05-18 14:29:41 +00:00
|
|
|
L[i * 6 - 2] = WaveletBase::GetDetailLength(L[i * 6 + 0]);
|
|
|
|
L[i * 6 - 1] = WaveletBase::GetApproxLength(L[i * 6 + 1]);
|
2016-08-12 00:49:48 +00:00
|
|
|
|
|
|
|
// cDv - overwrites previous value!
|
2017-05-18 14:29:41 +00:00
|
|
|
L[i * 6 - 0] = WaveletBase::GetDetailLength(L[i * 6 + 0]);
|
|
|
|
L[i * 6 + 1] = WaveletBase::GetDetailLength(L[i * 6 + 1]);
|
2016-08-12 00:49:48 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2017-02-21 23:56:57 +00:00
|
|
|
// Compute the bookkeeping array L for 3D DWT
|
2017-05-18 14:29:41 +00:00
|
|
|
void ComputeL3(vtkm::Id inX, vtkm::Id inY, vtkm::Id inZ, vtkm::Id nLev, std::vector<vtkm::Id>& L)
|
2017-02-21 23:56:57 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
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--)
|
2017-02-21 23:56:57 +00:00
|
|
|
{
|
|
|
|
// cLLL
|
2017-05-18 14:29:41 +00:00
|
|
|
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]);
|
2017-02-21 23:56:57 +00:00
|
|
|
|
|
|
|
// cLLH
|
2017-05-18 14:29:41 +00:00
|
|
|
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]);
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2017-02-21 23:56:57 +00:00
|
|
|
// cLHL
|
2017-05-18 14:29:41 +00:00
|
|
|
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];
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2017-02-21 23:56:57 +00:00
|
|
|
// cLHH
|
2017-05-18 14:29:41 +00:00
|
|
|
L[i * 21 - 12] = L[i * 21 - 21];
|
|
|
|
L[i * 21 - 11] = L[i * 21 - 14];
|
|
|
|
L[i * 21 - 10] = L[i * 21 - 16];
|
2017-05-22 20:47:13 +00:00
|
|
|
|
2017-02-21 23:56:57 +00:00
|
|
|
// cHLL
|
2017-05-18 14:29:41 +00:00
|
|
|
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];
|
2017-02-21 23:56:57 +00:00
|
|
|
|
|
|
|
// cHLH
|
2017-05-18 14:29:41 +00:00
|
|
|
L[i * 21 - 6] = L[i * 21 - 9];
|
|
|
|
L[i * 21 - 5] = L[i * 21 - 20];
|
|
|
|
L[i * 21 - 3] = L[i * 21 - 16];
|
2017-02-21 23:56:57 +00:00
|
|
|
|
|
|
|
// cHHL
|
2017-05-18 14:29:41 +00:00
|
|
|
L[i * 21 - 3] = L[i * 21 - 9];
|
|
|
|
L[i * 21 - 2] = L[i * 21 - 14];
|
|
|
|
L[i * 21 - 1] = L[i * 21 - 19];
|
2017-02-21 23:56:57 +00:00
|
|
|
|
|
|
|
// cHHH - overwrites previous value
|
2017-05-18 14:29:41 +00:00
|
|
|
L[i * 21 + 0] = L[i * 21 - 9];
|
|
|
|
L[i * 21 + 1] = L[i * 21 - 14];
|
|
|
|
L[i * 21 + 2] = L[i * 21 - 16];
|
2017-02-21 23:56:57 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2016-08-12 00:49:48 +00:00
|
|
|
// Compute the length of coefficients for 1D transforms
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::Id ComputeCoeffLength(std::vector<vtkm::Id>& L, vtkm::Id nLevels)
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::Id sum = L[0]; // 1st level cA
|
|
|
|
for (size_t i = 1; i <= size_t(nLevels); i++)
|
2017-05-18 17:51:35 +00:00
|
|
|
{
|
2016-07-27 20:28:38 +00:00
|
|
|
sum += L[i];
|
2017-05-18 17:51:35 +00:00
|
|
|
}
|
2016-07-27 20:28:38 +00:00
|
|
|
return sum;
|
|
|
|
}
|
2016-08-12 00:49:48 +00:00
|
|
|
// Compute the length of coefficients for 2D transforms
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::Id ComputeCoeffLength2(std::vector<vtkm::Id>& L, vtkm::Id nLevels)
|
2016-08-12 00:49:48 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::Id sum = (L[0] * L[1]); // 1st level cA
|
|
|
|
for (size_t i = 1; i <= size_t(nLevels); i++)
|
2016-08-12 00:49:48 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
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
|
2016-08-12 00:49:48 +00:00
|
|
|
}
|
|
|
|
return sum;
|
|
|
|
}
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
// Compute approximate coefficient length at a specific level
|
2017-05-18 14:29:41 +00:00
|
|
|
vtkm::Id GetApproxLengthLevN(vtkm::Id sigInLen, vtkm::Id levN)
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
|
|
|
vtkm::Id cALen = sigInLen;
|
2017-05-18 14:29:41 +00:00
|
|
|
for (vtkm::Id i = 0; i < levN; i++)
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2017-05-18 14:29:41 +00:00
|
|
|
cALen = WaveletBase::GetApproxLength(cALen);
|
|
|
|
if (cALen == 0)
|
2017-05-18 17:51:35 +00:00
|
|
|
{
|
2016-07-27 20:28:38 +00:00
|
|
|
return cALen;
|
2017-05-18 17:51:35 +00:00
|
|
|
}
|
2016-07-27 20:28:38 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
return cALen;
|
|
|
|
}
|
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
}; // class WaveletCompressor
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
} // namespace worklet
|
|
|
|
} // namespace vtkm
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2017-05-18 14:29:41 +00:00
|
|
|
#endif
|