vtk-m/vtkm/worklet/WaveletCompressor.h

649 lines
21 KiB
C
Raw Normal View History

//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
2019-04-15 23:24:21 +00:00
//
// 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.
//============================================================================
#ifndef vtk_m_worklet_waveletcompressor_h
#define vtk_m_worklet_waveletcompressor_h
#include <vtkm/worklet/wavelets/WaveletDWT.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayGetValues.h>
2017-05-18 14:29:41 +00:00
namespace vtkm
{
namespace worklet
{
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
template <typename SignalArrayType, typename CoeffArrayType>
2017-05-18 14:29:41 +00:00
VTKM_CONT vtkm::Id WaveDecompose(const SignalArrayType& sigIn, // Input
vtkm::Id nLevels, // n levels of DWT
CoeffArrayType& coeffOut,
std::vector<vtkm::Id>& L)
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
{
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
{
vtkm::cont::ArrayCopy(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;
std::vector<vtkm::Id> L1d(3, 0);
2016-08-02 20:09:29 +00:00
// Use an intermediate array
2018-02-22 13:29:13 +00:00
using OutputValueType = typename CoeffArrayType::ValueType;
using InterArrayType = vtkm::cont::ArrayHandle<OutputValueType>;
2016-08-02 20:09:29 +00:00
// Define a few more types
2018-02-22 13:29:13 +00:00
using IdArrayType = vtkm::cont::ArrayHandleCounting<vtkm::Id>;
using PermutArrayType = vtkm::cont::ArrayHandlePermutation<IdArrayType, CoeffArrayType>;
2016-08-02 20:09:29 +00:00
vtkm::cont::ArrayCopy(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
WaveletDWT::DWT1D(input, output, L1d);
2016-08-02 20:09:29 +00:00
2016-08-03 15:34:04 +00:00
// move intermediate results to final array
WaveletBase::DeviceCopyStartX(output, coeffOut, cptr);
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;
}
return 0;
}
// Multi-level 1D wavelet reconstruction
template <typename CoeffArrayType, typename SignalArrayType>
2017-05-18 14:29:41 +00:00
VTKM_CONT vtkm::Id WaveReconstruct(const CoeffArrayType& coeffIn, // Input
vtkm::Id nLevels, // n levels of DWT
std::vector<vtkm::Id>& L,
SignalArrayType& sigOut)
{
2017-05-18 14:29:41 +00:00
VTKM_ASSERT(nLevels > 0);
vtkm::Id LLength = nLevels + 2;
2017-05-18 14:29:41 +00:00
VTKM_ASSERT(vtkm::Id(L.size()) == LLength);
2017-05-18 14:29:41 +00:00
std::vector<vtkm::Id> L1d(3, 0); // three elements
L1d[0] = L[0];
L1d[1] = L[1];
2018-02-22 13:29:13 +00:00
using OutValueType = typename SignalArrayType::ValueType;
using OutArrayBasic = vtkm::cont::ArrayHandle<OutValueType>;
using IdArrayType = vtkm::cont::ArrayHandleCounting<vtkm::Id>;
using PermutArrayType = vtkm::cont::ArrayHandlePermutation<IdArrayType, SignalArrayType>;
2016-08-02 20:35:00 +00:00
vtkm::cont::ArrayCopy(coeffIn, sigOut);
2017-05-18 14:29:41 +00:00
for (vtkm::Id i = 1; i <= nLevels; i++)
{
2017-05-18 14:29:41 +00:00
L1d[2] = this->GetApproxLengthLevN(L[size_t(LLength - 1)], nLevels - i);
// 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
// Make an output array
2016-08-02 20:35:00 +00:00
OutArrayBasic output;
2017-05-22 20:47:13 +00:00
WaveletDWT::IDWT1D(input, L1d, output);
2017-05-18 14:29:41 +00:00
VTKM_ASSERT(output.GetNumberOfValues() == L1d[2]);
// Move output to intermediate array
WaveletBase::DeviceCopyStartX(output, sigOut, 0);
L1d[0] = L1d[2];
2017-05-18 14:29:41 +00:00
L1d[1] = L[size_t(i + 1)];
}
return 0;
}
2017-02-21 23:56:57 +00:00
// Multi-level 3D wavelet decomposition
template <typename InArrayType, typename OutArrayType>
2017-05-18 14:29:41 +00:00
VTKM_CONT vtkm::Float64 WaveDecompose3D(InArrayType& sigIn, // Input
vtkm::Id nLevels, // n levels of DWT
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?
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
{
vtkm::cont::ArrayCopy(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
2018-02-22 13:29:13 +00:00
using OutValueType = typename OutArrayType::ValueType;
using OutBasicArray = vtkm::cont::ArrayHandle<OutValueType>;
2017-02-21 23:56:57 +00:00
// 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);
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-21 23:56:57 +00:00
OutBasicArray tempOutput;
computationTime += WaveletDWT::DWT3D(
coeffOut, inX, inY, inZ, 0, 0, 0, currentLenX, currentLenY, currentLenZ, tempOutput, false);
2017-02-21 23:56:57 +00:00
// copy results to coeffOut
WaveletBase::DeviceCubeCopyTo(
tempOutput, currentLenX, currentLenY, currentLenZ, coeffOut, inX, inY, inZ, 0, 0, 0);
}
2017-02-21 23:56:57 +00:00
return computationTime;
}
// Multi-level 3D wavelet reconstruction
template <typename InArrayType, typename OutArrayType>
2017-05-18 14:29:41 +00:00
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?
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
}
2018-02-22 13:29:13 +00:00
using OutValueType = typename OutArrayType::ValueType;
using OutBasicArray = vtkm::cont::ArrayHandle<OutValueType>;
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
{
vtkm::cont::ArrayCopy(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
{
vtkm::cont::ArrayCopy(arrIn, outBuffer);
2017-05-18 17:51:35 +00:00
}
2017-02-21 23:56:57 +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
// 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
{
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);
2017-02-21 23:56:57 +00:00
// copy back reconstructed block
WaveletBase::DeviceCubeCopyTo(
tempOutput, L3d[24], L3d[25], L3d[26], outBuffer, inX, inY, inZ, 0, 0, 0);
2017-05-22 20:47:13 +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
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];
computationTime += WaveletDWT::IDWT3D(outBuffer, inX, inY, inZ, 0, 0, 0, L3d, arrOut, true);
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
template <typename InArrayType, typename OutArrayType>
2017-05-18 14:29:41 +00:00
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)
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
{
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
{
vtkm::cont::ArrayCopy(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;
2018-02-22 13:29:13 +00:00
using OutValueType = typename OutArrayType::ValueType;
using OutBasicArray = vtkm::cont::ArrayHandle<OutValueType>;
2016-08-12 00:49:48 +00:00
// First level transform operates writes to the output array
computationTime += WaveletDWT::DWT2D(
sigIn, currentLenX, currentLenY, 0, 0, currentLenX, currentLenY, coeffOut, L2d);
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
// 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
{
OutBasicArray tempOutput;
computationTime +=
WaveletDWT::DWT2D(coeffOut, inX, inY, 0, 0, currentLenX, currentLenY, tempOutput, L2d);
// copy results to coeffOut
WaveletBase::DeviceRectangleCopyTo(
tempOutput, currentLenX, currentLenY, coeffOut, inX, inY, 0, 0);
// update currentLen
2017-05-18 14:29:41 +00:00
currentLenX = WaveletBase::GetApproxLength(currentLenX);
currentLenY = WaveletBase::GetApproxLength(currentLenY);
}
return computationTime;
}
2016-08-12 00:49:48 +00:00
// Multi-level 2D wavelet reconstruction
template <typename InArrayType, typename OutArrayType>
2017-05-18 14:29:41 +00:00
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)
{
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))
{
throw vtkm::cont::ErrorBadValue("Number of levels of transform is not supported! ");
2016-08-12 00:49:48 +00:00
}
2018-02-22 13:29:13 +00:00
using OutValueType = typename OutArrayType::ValueType;
using OutBasicArray = vtkm::cont::ArrayHandle<OutValueType>;
2016-09-22 03:25:56 +00:00
vtkm::Float64 computationTime = 0.0;
2017-05-22 20:47:13 +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
{
vtkm::cont::ArrayCopy(arrIn, arrOut);
2017-05-22 20:47:13 +00:00
return 0;
}
2016-09-22 03:25:56 +00:00
else
2017-05-18 17:51:35 +00:00
{
vtkm::cont::ArrayCopy(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);
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++)
{
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)
2017-05-18 14:29:41 +00:00
OutBasicArray tempOutput;
// IDWT
computationTime += WaveletDWT::IDWT2D(outBuffer, inX, inY, 0, 0, L2d, tempOutput);
// copy back reconstructed block
WaveletBase::DeviceRectangleCopyTo(tempOutput, L2d[8], L2d[9], outBuffer, inX, inY, 0, 0);
2017-05-22 20:47:13 +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];
}
// 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];
computationTime += WaveletDWT::IDWT2D(outBuffer, inX, inY, 0, 0, L2d, arrOut);
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
}
// Squash coefficients smaller than a threshold
template <typename CoeffArrayType>
vtkm::Id SquashCoefficients(CoeffArrayType& coeffIn, vtkm::Float64 ratio)
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();
2018-02-22 13:29:13 +00:00
using ValueType = typename CoeffArrayType::ValueType;
using CoeffArrayBasic = vtkm::cont::ArrayHandle<ValueType>;
2016-07-29 01:12:57 +00:00
CoeffArrayBasic sortedArray;
vtkm::cont::ArrayCopy(coeffIn, sortedArray);
2017-03-09 03:29:27 +00:00
WaveletBase::DeviceSort(sortedArray);
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>(vtkm::cont::ArrayGetValue(n, sortedArray));
2017-05-18 14:29:41 +00:00
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
}
2018-02-22 13:29:13 +00:00
using ThresholdType = vtkm::worklet::wavelets::ThresholdWorklet;
2017-05-18 14:29:41 +00:00
ThresholdType thresholdWorklet(nthVal);
vtkm::worklet::DispatcherMapField<ThresholdType> dispatcher(thresholdWorklet);
2017-05-18 14:29:41 +00:00
dispatcher.Invoke(coeffIn);
2017-05-22 20:47:13 +00:00
}
2016-07-28 22:40:37 +00:00
return 0;
}
2016-07-29 03:07:05 +00:00
// Report statistics on reconstructed array
template <typename ArrayType>
vtkm::Id EvaluateReconstruction(const ArrayType& original, const ArrayType& reconstruct)
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);
2016-07-29 03:07:05 +00:00
2018-02-22 13:29:13 +00:00
using ValueType = typename ArrayType::ValueType;
using ArrayBasic = vtkm::cont::ArrayHandle<ValueType>;
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
2018-02-22 13:29:13 +00:00
using DifferencerWorklet = vtkm::worklet::wavelets::Differencer;
2016-07-29 03:07:05 +00:00
DifferencerWorklet dw;
vtkm::worklet::DispatcherMapField<DifferencerWorklet> dwDispatcher(dw);
2017-05-18 14:29:41 +00:00
dwDispatcher.Invoke(original, reconstruct, errorArray);
2016-07-29 18:27:37 +00:00
2018-02-22 13:29:13 +00:00
using SquareWorklet = vtkm::worklet::wavelets::SquareWorklet;
2016-07-29 18:27:37 +00:00
SquareWorklet sw;
vtkm::worklet::DispatcherMapField<SquareWorklet> swDispatcher(sw);
2017-05-18 14:29:41 +00:00
swDispatcher.Invoke(errorArray, errorSquare);
2016-07-29 03:07:05 +00:00
VAL varErr = WaveletBase::DeviceCalculateVariance(errorArray);
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
VAL origMax = WaveletBase::DeviceMax(original);
VAL origMin = WaveletBase::DeviceMin(original);
VAL errorMax = WaveletBase::DeviceMaxAbs(errorArray);
2017-05-18 14:29:41 +00:00
VAL range = origMax - origMin;
2016-07-29 03:07:05 +00:00
VAL squareSum = WaveletBase::DeviceSum(errorSquare);
2017-05-18 14:29:41 +00:00
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)
{
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--)
{
2017-05-18 14:29:41 +00:00
L[i - 1] = WaveletBase::GetApproxLength(L[i]);
L[i] = WaveletBase::GetDetailLength(L[i]);
}
}
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)
{
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
{
sum += L[i];
2017-05-18 17:51:35 +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;
}
// 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)
{
vtkm::Id cALen = sigInLen;
2017-05-18 14:29:41 +00:00
for (vtkm::Id i = 0; i < levN; i++)
{
2017-05-18 14:29:41 +00:00
cALen = WaveletBase::GetApproxLength(cALen);
if (cALen == 0)
2017-05-18 17:51:35 +00:00
{
return cALen;
2017-05-18 17:51:35 +00:00
}
}
return cALen;
}
2017-05-18 14:29:41 +00:00
}; // class WaveletCompressor
2017-05-18 14:29:41 +00:00
} // namespace worklet
} // namespace vtkm
2017-05-18 14:29:41 +00:00
#endif