8859636672
VTK-m has been updated to replace old per device worklet testing executables with a device dependent shared library so that it's able to accept a device adapter at runtime. Meanwhile, it updates the testing infrastructure APIs. vtkm::cont::testing::Run function would call ForceDevice when needed and if users need the device adapter info at runtime, RunOnDevice function would pass the adapter into the functor. Optional Parser is bumped from 1.3 to 1.7.
658 lines
22 KiB
C++
658 lines
22 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
|
|
// Copyright 2014 UT-Battelle, LLC.
|
|
// Copyright 2014 Los Alamos National Security.
|
|
//
|
|
// Under the terms of Contract DE-NA0003525 with NTESS,
|
|
// 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>
|
|
|
|
#include <vtkm/cont/ArrayCopy.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>
|
|
VTKM_CONT vtkm::Id WaveDecompose(const SignalArrayType& sigIn, // Input
|
|
vtkm::Id nLevels, // n levels of DWT
|
|
CoeffArrayType& coeffOut,
|
|
std::vector<vtkm::Id>& L)
|
|
{
|
|
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::ArrayCopy(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
|
|
using OutputValueType = typename CoeffArrayType::ValueType;
|
|
using InterArrayType = vtkm::cont::ArrayHandle<OutputValueType>;
|
|
|
|
// Define a few more types
|
|
using IdArrayType = vtkm::cont::ArrayHandleCounting<vtkm::Id>;
|
|
using PermutArrayType = vtkm::cont::ArrayHandlePermutation<IdArrayType, CoeffArrayType>;
|
|
|
|
vtkm::cont::ArrayCopy(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);
|
|
|
|
// move intermediate results to final array
|
|
WaveletBase::DeviceCopyStartX(output, coeffOut, cptr);
|
|
|
|
// update pseudo pointers
|
|
len = cALen;
|
|
cALen = WaveletBase::GetApproxLength(cALen);
|
|
sigInPtr = cptr;
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
// Multi-level 1D wavelet reconstruction
|
|
template <typename CoeffArrayType, typename SignalArrayType>
|
|
VTKM_CONT vtkm::Id WaveReconstruct(const CoeffArrayType& coeffIn, // Input
|
|
vtkm::Id nLevels, // n levels of DWT
|
|
std::vector<vtkm::Id>& L,
|
|
SignalArrayType& sigOut)
|
|
{
|
|
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];
|
|
|
|
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>;
|
|
|
|
vtkm::cont::ArrayCopy(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);
|
|
VTKM_ASSERT(output.GetNumberOfValues() == L1d[2]);
|
|
|
|
// Move output to intermediate array
|
|
WaveletBase::DeviceCopyStartX(output, sigOut, 0);
|
|
|
|
L1d[0] = L1d[2];
|
|
L1d[1] = L[size_t(i + 1)];
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
// Multi-level 3D wavelet decomposition
|
|
template <typename InArrayType, typename OutArrayType>
|
|
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 on devices?
|
|
{
|
|
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::ArrayCopy(sigIn, coeffOut);
|
|
return 0;
|
|
}
|
|
|
|
vtkm::Id currentLenX = inX;
|
|
vtkm::Id currentLenY = inY;
|
|
vtkm::Id currentLenZ = inZ;
|
|
|
|
using OutValueType = typename OutArrayType::ValueType;
|
|
using OutBasicArray = vtkm::cont::ArrayHandle<OutValueType>;
|
|
|
|
// 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);
|
|
|
|
// 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);
|
|
|
|
// copy results to coeffOut
|
|
WaveletBase::DeviceCubeCopyTo(
|
|
tempOutput, currentLenX, currentLenY, currentLenZ, coeffOut, inX, inY, inZ, 0, 0, 0);
|
|
}
|
|
|
|
return computationTime;
|
|
}
|
|
|
|
// Multi-level 3D wavelet reconstruction
|
|
template <typename InArrayType, typename OutArrayType>
|
|
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?
|
|
{
|
|
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! ");
|
|
}
|
|
using OutValueType = typename OutArrayType::ValueType;
|
|
using OutBasicArray = vtkm::cont::ArrayHandle<OutValueType>;
|
|
vtkm::Float64 computationTime = 0.0;
|
|
|
|
OutBasicArray outBuffer;
|
|
if (nLevels == 0) // 0 levels means no transform
|
|
{
|
|
vtkm::cont::ArrayCopy(arrIn, arrOut);
|
|
return 0;
|
|
}
|
|
else if (discardArrIn)
|
|
{
|
|
outBuffer = arrIn;
|
|
}
|
|
else
|
|
{
|
|
vtkm::cont::ArrayCopy(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);
|
|
|
|
// copy back reconstructed block
|
|
WaveletBase::DeviceCubeCopyTo(
|
|
tempOutput, L3d[24], L3d[25], L3d[26], outBuffer, inX, inY, inZ, 0, 0, 0);
|
|
|
|
// 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);
|
|
|
|
return computationTime;
|
|
}
|
|
|
|
// Multi-level 2D wavelet decomposition
|
|
template <typename InArrayType, typename OutArrayType>
|
|
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)
|
|
{
|
|
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::ArrayCopy(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;
|
|
|
|
using OutValueType = typename OutArrayType::ValueType;
|
|
using OutBasicArray = vtkm::cont::ArrayHandle<OutValueType>;
|
|
|
|
// First level transform operates writes to the output array
|
|
computationTime += WaveletDWT::DWT2D(
|
|
sigIn, currentLenX, currentLenY, 0, 0, currentLenX, currentLenY, coeffOut, L2d);
|
|
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);
|
|
|
|
// copy results to coeffOut
|
|
WaveletBase::DeviceRectangleCopyTo(
|
|
tempOutput, currentLenX, currentLenY, coeffOut, inX, inY, 0, 0);
|
|
|
|
// update currentLen
|
|
currentLenX = WaveletBase::GetApproxLength(currentLenX);
|
|
currentLenY = WaveletBase::GetApproxLength(currentLenY);
|
|
}
|
|
|
|
return computationTime;
|
|
}
|
|
|
|
// Multi-level 2D wavelet reconstruction
|
|
template <typename InArrayType, typename OutArrayType>
|
|
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();
|
|
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! ");
|
|
}
|
|
using OutValueType = typename OutArrayType::ValueType;
|
|
using OutBasicArray = vtkm::cont::ArrayHandle<OutValueType>;
|
|
vtkm::Float64 computationTime = 0.0;
|
|
|
|
OutBasicArray outBuffer;
|
|
if (nLevels == 0) // 0 levels means no transform
|
|
{
|
|
vtkm::cont::ArrayCopy(arrIn, arrOut);
|
|
return 0;
|
|
}
|
|
else
|
|
{
|
|
vtkm::cont::ArrayCopy(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);
|
|
|
|
// copy back reconstructed block
|
|
WaveletBase::DeviceRectangleCopyTo(tempOutput, L2d[8], L2d[9], outBuffer, inX, inY, 0, 0);
|
|
|
|
// 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);
|
|
|
|
return computationTime;
|
|
}
|
|
|
|
// Squash coefficients smaller than a threshold
|
|
template <typename CoeffArrayType>
|
|
vtkm::Id SquashCoefficients(CoeffArrayType& coeffIn, vtkm::Float64 ratio)
|
|
{
|
|
if (ratio > 1.0)
|
|
{
|
|
vtkm::Id coeffLen = coeffIn.GetNumberOfValues();
|
|
using ValueType = typename CoeffArrayType::ValueType;
|
|
using CoeffArrayBasic = vtkm::cont::ArrayHandle<ValueType>;
|
|
CoeffArrayBasic sortedArray;
|
|
vtkm::cont::ArrayCopy(coeffIn, sortedArray);
|
|
|
|
WaveletBase::DeviceSort(sortedArray);
|
|
|
|
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;
|
|
}
|
|
|
|
using ThresholdType = vtkm::worklet::wavelets::ThresholdWorklet;
|
|
ThresholdType thresholdWorklet(nthVal);
|
|
vtkm::worklet::DispatcherMapField<ThresholdType> dispatcher(thresholdWorklet);
|
|
dispatcher.Invoke(coeffIn);
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
// Report statistics on reconstructed array
|
|
template <typename ArrayType>
|
|
vtkm::Id EvaluateReconstruction(const ArrayType& original, const ArrayType& reconstruct)
|
|
{
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
VAL VarOrig = WaveletBase::DeviceCalculateVariance(original);
|
|
|
|
using ValueType = typename ArrayType::ValueType;
|
|
using ArrayBasic = vtkm::cont::ArrayHandle<ValueType>;
|
|
ArrayBasic errorArray, errorSquare;
|
|
|
|
// Use a worklet to calculate point-wise error, and its square
|
|
using DifferencerWorklet = vtkm::worklet::wavelets::Differencer;
|
|
DifferencerWorklet dw;
|
|
vtkm::worklet::DispatcherMapField<DifferencerWorklet> dwDispatcher(dw);
|
|
dwDispatcher.Invoke(original, reconstruct, errorArray);
|
|
|
|
using SquareWorklet = vtkm::worklet::wavelets::SquareWorklet;
|
|
SquareWorklet sw;
|
|
vtkm::worklet::DispatcherMapField<SquareWorklet> swDispatcher(sw);
|
|
swDispatcher.Invoke(errorArray, errorSquare);
|
|
|
|
VAL varErr = WaveletBase::DeviceCalculateVariance(errorArray);
|
|
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);
|
|
VAL origMin = WaveletBase::DeviceMin(original);
|
|
VAL errorMax = WaveletBase::DeviceMaxAbs(errorArray);
|
|
VAL range = origMax - origMin;
|
|
|
|
VAL squareSum = WaveletBase::DeviceSum(errorSquare);
|
|
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
|