3133 lines
111 KiB
C++
3133 lines
111 KiB
C++
//============================================================================
|
|
// Copyright (c) Kitware, Inc.
|
|
// All rights reserved.
|
|
// See LICENSE.txt for details.
|
|
// This software is distributed WITHOUT ANY WARRANTY; without even
|
|
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|
// PURPOSE. See the above copyright notice for more information.
|
|
//
|
|
// Copyright 2014 Sandia Corporation.
|
|
// Copyright 2014 UT-Battelle, LLC.
|
|
// Copyright 2014 Los Alamos National Security.
|
|
//
|
|
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
|
|
// the U.S. Government retains certain rights in this software.
|
|
//
|
|
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
|
|
// Laboratory (LANL), the U.S. Government retains certain rights in
|
|
// this software.
|
|
//============================================================================
|
|
|
|
#ifndef vtk_m_worklet_Wavelets_h
|
|
#define vtk_m_worklet_Wavelets_h
|
|
|
|
#include <vtkm/worklet/WorkletMapField.h>
|
|
#include <vtkm/worklet/DispatcherMapField.h>
|
|
|
|
#include <vtkm/Math.h>
|
|
|
|
namespace vtkm {
|
|
namespace worklet {
|
|
namespace wavelets {
|
|
|
|
enum DWTMode { // boundary extension modes
|
|
SYMH,
|
|
SYMW,
|
|
ASYMH,
|
|
ASYMW
|
|
};
|
|
|
|
enum ExtensionDirection { // which side of a cube to extend
|
|
LEFT, // X direction
|
|
RIGHT, // X direction Y
|
|
TOP, // Y direction | Z
|
|
BOTTOM, // Y direction | /
|
|
FRONT, // Z direction | /
|
|
BACK // Z direction |/________ X
|
|
};
|
|
|
|
|
|
// Worklet for 3D signal extension
|
|
// It operates on a specified part of a big cube
|
|
class ExtensionWorklet3D : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension
|
|
WholeArrayIn < ScalarAll > ); // signal
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
ExtensionWorklet3D ( vtkm::Id extdimX, vtkm::Id extdimY, vtkm::Id extdimZ,
|
|
vtkm::Id sigdimX, vtkm::Id sigdimY, vtkm::Id sigdimZ,
|
|
vtkm::Id sigstartX, vtkm::Id sigstartY, vtkm::Id sigstartZ,
|
|
vtkm::Id sigpretendX, vtkm::Id sigpretendY, vtkm::Id sigpretendZ,
|
|
DWTMode m, // SYMH, SYMW, etc.
|
|
ExtensionDirection dir,
|
|
bool pad_zero )
|
|
:
|
|
extDimX( extdimX ), extDimY( extdimY ), extDimZ( extdimZ ),
|
|
sigDimX( sigdimX ), sigDimY( sigdimY ), sigDimZ( sigdimZ ),
|
|
sigStartX( sigstartX ), sigStartY( sigstartY ), sigStartZ( sigstartZ ),
|
|
sigPretendDimX( sigpretendX ),
|
|
sigPretendDimY( sigpretendY ),
|
|
sigPretendDimZ( sigpretendZ ),
|
|
mode(m),
|
|
direction( dir ),
|
|
padZero( pad_zero )
|
|
{}
|
|
|
|
// Index translation helper
|
|
VTKM_EXEC_CONT
|
|
void Ext1Dto3D ( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y, vtkm::Id &z ) const
|
|
{
|
|
z = idx / (extDimX * extDimY);
|
|
y = (idx - z * extDimX * extDimY) / extDimX;
|
|
x = idx % extDimX;
|
|
}
|
|
|
|
// Index translation helper
|
|
VTKM_EXEC_CONT
|
|
vtkm::Id Sig3Dto1D( vtkm::Id x, vtkm::Id y, vtkm::Id z) const
|
|
{
|
|
return z * sigDimX * sigDimY + y * sigDimX + x;
|
|
}
|
|
|
|
// Index translation helper
|
|
VTKM_EXEC_CONT
|
|
vtkm::Id SigPretend3Dto1D( vtkm::Id x, vtkm::Id y, vtkm::Id z ) const
|
|
{
|
|
return (z + sigStartZ) * sigDimX * sigDimY + (y + sigStartY) * sigDimX + x + sigStartX;
|
|
}
|
|
|
|
template< typename PortalOutType, typename PortalInType >
|
|
VTKM_EXEC
|
|
void operator()( PortalOutType &portalOut,
|
|
const PortalInType &portalIn,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
vtkm::Id extX, extY, extZ;
|
|
vtkm::Id sigPretendX, sigPretendY, sigPretendZ;
|
|
Ext1Dto3D( workIndex, extX, extY, extZ );
|
|
typename PortalOutType::ValueType sym = 1.0;
|
|
if( mode == ASYMH || mode == ASYMW )
|
|
sym = -1.0;
|
|
if( direction == LEFT )
|
|
{
|
|
sigPretendY = extY;
|
|
sigPretendZ = extZ;
|
|
if( mode == SYMH || mode == ASYMH )
|
|
sigPretendX = extDimX - extX - 1;
|
|
else // mode == SYMW || mode == ASYMW
|
|
sigPretendX = extDimX - extX;
|
|
}
|
|
else if( direction == RIGHT )
|
|
{
|
|
sigPretendY = extY;
|
|
sigPretendZ = extZ;
|
|
if( mode == SYMH || mode == ASYMH )
|
|
sigPretendX = sigPretendDimX - extX - 1;
|
|
else
|
|
sigPretendX = sigPretendDimX - extX - 2;
|
|
if( padZero )
|
|
sigPretendX++;
|
|
}
|
|
else if( direction == TOP )
|
|
{
|
|
sigPretendX = extX;
|
|
sigPretendZ = extZ;
|
|
if( mode == SYMH || mode == ASYMH )
|
|
sigPretendY = extDimY - extY - 1;
|
|
else // mode == SYMW || mode == ASYMW
|
|
sigPretendY = extDimY - extY;
|
|
}
|
|
else if( direction == BOTTOM )
|
|
{
|
|
sigPretendX = extX;
|
|
sigPretendZ = extZ;
|
|
if( mode == SYMH || mode == ASYMH )
|
|
sigPretendY = sigPretendDimY - extY - 1;
|
|
else
|
|
sigPretendY = sigPretendDimY - extY - 2;
|
|
if( padZero )
|
|
sigPretendY++;
|
|
}
|
|
else if( direction == FRONT )
|
|
{
|
|
sigPretendX = extX;
|
|
sigPretendY = extY;
|
|
if( mode == SYMH || mode == ASYMH )
|
|
sigPretendZ = extDimZ - extZ - 1;
|
|
else // mode == SYMW || mode == ASYMW
|
|
sigPretendZ = extDimZ - extZ;
|
|
}
|
|
else if( direction == BACK )
|
|
{
|
|
sigPretendX = extX;
|
|
sigPretendY = extY;
|
|
if( mode == SYMH || mode == ASYMH )
|
|
sigPretendZ = sigPretendDimZ - extZ - 1;
|
|
else
|
|
sigPretendZ = sigPretendDimZ - extZ - 2;
|
|
if( padZero )
|
|
sigPretendZ++;
|
|
}
|
|
else
|
|
{
|
|
sigPretendX = sigPretendY = sigPretendZ = 0; // so the compiler doesn't complain
|
|
vtkm::cont::ErrorControlInternal("Invalid extension mode for cubes!");
|
|
}
|
|
|
|
if( sigPretendX == sigPretendDimX || // decides to pad a zero
|
|
sigPretendY == sigPretendDimY ||
|
|
sigPretendZ == sigPretendDimZ )
|
|
portalOut.Set( workIndex, 0.0 );
|
|
else
|
|
portalOut.Set( workIndex, sym * portalIn.Get(
|
|
SigPretend3Dto1D(sigPretendX, sigPretendY, sigPretendZ) ));
|
|
}
|
|
|
|
private:
|
|
const vtkm::Id extDimX, extDimY, extDimZ, sigDimX, sigDimY, sigDimZ;
|
|
const vtkm::Id sigStartX, sigStartY, sigStartZ; // defines a small cube to work on
|
|
const vtkm::Id sigPretendDimX, sigPretendDimY, sigPretendDimZ; // small cube dims
|
|
const DWTMode mode;
|
|
const ExtensionDirection direction;
|
|
const bool padZero; // treat sigIn as having a zero at the end
|
|
};
|
|
|
|
|
|
//=============================================================================
|
|
|
|
// Y
|
|
//
|
|
// | Z
|
|
// | /
|
|
// | /
|
|
// | /
|
|
// | /
|
|
// | /
|
|
// |/------------- X
|
|
//
|
|
// The following 3 classes perform the same functionaliry in 3 directions
|
|
class IndexTranslator3CubesLeftRight
|
|
{
|
|
public:
|
|
IndexTranslator3CubesLeftRight (
|
|
vtkm::Id x_1, vtkm::Id y_1, vtkm::Id z_1,
|
|
vtkm::Id x_2, vtkm::Id y_2, vtkm::Id z_2,
|
|
vtkm::Id startx_2, vtkm::Id starty_2, vtkm::Id startz_2,
|
|
vtkm::Id pretendx_2, vtkm::Id pretendy_2, vtkm::Id pretendz_2,
|
|
vtkm::Id x_3, vtkm::Id y_3, vtkm::Id z_3 )
|
|
:
|
|
dimX1(x_1), dimY1(y_1), dimZ1(z_1),
|
|
dimX2(x_2), dimY2(y_2), dimZ2(z_2),
|
|
startX2( startx_2 ), startY2( starty_2 ), startZ2(startz_2),
|
|
pretendDimX2(pretendx_2), pretendDimY2(pretendy_2), pretendDimZ2(pretendz_2),
|
|
dimX3(x_3), dimY3(y_3), dimZ3(z_3)
|
|
{ (void)dimY2; }
|
|
|
|
VTKM_EXEC_CONT
|
|
void Translate3Dto1D( vtkm::Id inX, vtkm::Id inY, vtkm::Id inZ, // 3D indices as input
|
|
vtkm::Id &cube, vtkm::Id &idx ) const // which cube, and idx of that cube
|
|
{
|
|
if ( dimX1 <= inX && inX < (dimX1 + pretendDimX2) )
|
|
{
|
|
vtkm::Id inX_local = inX - dimX1;
|
|
cube = 2;
|
|
idx = (inZ + startZ2) * dimX2 * dimY2 + (inY + startY2) * dimX2 + (inX_local + startX2);
|
|
}
|
|
else if ( 0 <= inX && inX < dimX1 )
|
|
{
|
|
cube = 1;
|
|
idx = inZ * dimX1 * dimY1 + inY * dimX1 + inX;
|
|
}
|
|
else if ( (dimX1 + pretendDimX2) <= inX && inX < (dimX1 + pretendDimX2 + dimX3) )
|
|
{
|
|
vtkm::Id inX_local = inX - dimX1 - pretendDimX2;
|
|
cube = 3;
|
|
idx = inZ * dimX3 * dimY3 + inY * dimX3 + inX_local;
|
|
}
|
|
else
|
|
vtkm::cont::ErrorControlInternal("Invalid index!");
|
|
}
|
|
|
|
private:
|
|
const vtkm::Id dimX1, dimY1, dimZ1; // left extension
|
|
const vtkm::Id dimX2, dimY2, dimZ2; // actual signal dims
|
|
const vtkm::Id startX2, startY2, startZ2, pretendDimX2, pretendDimY2, pretendDimZ2;
|
|
const vtkm::Id dimX3, dimY3, dimZ3; // right extension
|
|
};
|
|
|
|
class IndexTranslator3CubesTopDown
|
|
{
|
|
public:
|
|
IndexTranslator3CubesTopDown (
|
|
vtkm::Id x_1, vtkm::Id y_1, vtkm::Id z_1,
|
|
vtkm::Id x_2, vtkm::Id y_2, vtkm::Id z_2,
|
|
vtkm::Id startx_2, vtkm::Id starty_2, vtkm::Id startz_2,
|
|
vtkm::Id pretendx_2, vtkm::Id pretendy_2, vtkm::Id pretendz_2,
|
|
vtkm::Id x_3, vtkm::Id y_3, vtkm::Id z_3 )
|
|
:
|
|
dimX1(x_1), dimY1(y_1), dimZ1(z_1),
|
|
dimX2(x_2), dimY2(y_2), dimZ2(z_2),
|
|
startX2( startx_2 ), startY2( starty_2 ), startZ2(startz_2),
|
|
pretendDimX2(pretendx_2), pretendDimY2(pretendy_2), pretendDimZ2(pretendz_2),
|
|
dimX3(x_3), dimY3(y_3), dimZ3(z_3)
|
|
{ }
|
|
|
|
VTKM_EXEC_CONT
|
|
void Translate3Dto1D( vtkm::Id inX, vtkm::Id inY, vtkm::Id inZ, // 3D indices as input
|
|
vtkm::Id &cube, vtkm::Id &idx ) const // which cube, and idx of that cube
|
|
{
|
|
if ( dimY1 <= inY && inY < (dimY1 + pretendDimY2) )
|
|
{
|
|
vtkm::Id inY_local = inY - dimY1;
|
|
cube = 2;
|
|
idx = (inZ + startZ2) * dimX2 * dimY2 + (inY_local + startY2) * dimX2 + inX + startX2;
|
|
}
|
|
else if ( 0 <= inY && inY < dimY1 )
|
|
{
|
|
cube = 1;
|
|
idx = inZ * dimX1 * dimY1 + inY * dimX1 + inX;
|
|
}
|
|
else if ( (dimY1 + pretendDimY2) <= inY && inY < (dimY1 + pretendDimY2 + dimY3) )
|
|
{
|
|
vtkm::Id inY_local = inY - dimY1 - pretendDimY2;
|
|
cube = 3;
|
|
idx = inZ * dimX3 * dimY3 + inY_local * dimX3 + inX;
|
|
}
|
|
else
|
|
vtkm::cont::ErrorControlInternal("Invalid index!");
|
|
}
|
|
|
|
private:
|
|
const vtkm::Id dimX1, dimY1, dimZ1; // left extension
|
|
const vtkm::Id dimX2, dimY2, dimZ2; // actual signal dims
|
|
const vtkm::Id startX2, startY2, startZ2, pretendDimX2, pretendDimY2, pretendDimZ2;
|
|
const vtkm::Id dimX3, dimY3, dimZ3; // right extension
|
|
};
|
|
|
|
class IndexTranslator3CubesFrontBack
|
|
{
|
|
public:
|
|
IndexTranslator3CubesFrontBack (
|
|
vtkm::Id x_1, vtkm::Id y_1, vtkm::Id z_1,
|
|
vtkm::Id x_2, vtkm::Id y_2, vtkm::Id z_2,
|
|
vtkm::Id startx_2, vtkm::Id starty_2, vtkm::Id startz_2,
|
|
vtkm::Id pretendx_2, vtkm::Id pretendy_2, vtkm::Id pretendz_2,
|
|
vtkm::Id x_3, vtkm::Id y_3, vtkm::Id z_3 )
|
|
:
|
|
dimX1(x_1), dimY1(y_1), dimZ1(z_1),
|
|
dimX2(x_2), dimY2(y_2), dimZ2(z_2),
|
|
startX2( startx_2 ), startY2( starty_2 ), startZ2(startz_2),
|
|
pretendDimX2(pretendx_2), pretendDimY2(pretendy_2), pretendDimZ2(pretendz_2),
|
|
dimX3(x_3), dimY3(y_3), dimZ3(z_3)
|
|
{ }
|
|
|
|
VTKM_EXEC_CONT
|
|
void Translate3Dto1D( vtkm::Id inX, vtkm::Id inY, vtkm::Id inZ, // 3D indices as input
|
|
vtkm::Id &cube, vtkm::Id &idx ) const // which cube, and idx of that cube
|
|
{
|
|
if ( dimZ1 <= inZ && inZ < (dimZ1 + pretendDimZ2) )
|
|
{
|
|
vtkm::Id inZ_local = inZ - dimZ1;
|
|
cube = 2;
|
|
idx = (inZ_local + startZ2) * dimX2 * dimY2 + (inY + startY2) * dimX2 + inX + startX2;
|
|
}
|
|
else if ( 0 <= inZ && inZ < dimZ1 )
|
|
{
|
|
cube = 1;
|
|
idx = inZ * dimX1 * dimY1 + inY * dimX1 + inX;
|
|
}
|
|
else if ( (dimZ1 + pretendDimZ2) <= inZ && inZ < (dimZ1 + pretendDimZ2 + dimZ3) )
|
|
{
|
|
vtkm::Id inZ_local = inZ - dimZ1 - pretendDimZ2;
|
|
cube = 3;
|
|
idx = inZ_local * dimX3 * dimY3 + inY * dimX3 + inX;
|
|
}
|
|
else
|
|
vtkm::cont::ErrorControlInternal("Invalid index!");
|
|
}
|
|
|
|
private:
|
|
const vtkm::Id dimX1, dimY1, dimZ1; // left extension
|
|
const vtkm::Id dimX2, dimY2, dimZ2; // actual signal dims
|
|
const vtkm::Id startX2, startY2, startZ2, pretendDimX2, pretendDimY2, pretendDimZ2;
|
|
const vtkm::Id dimX3, dimY3, dimZ3; // right extension
|
|
};
|
|
|
|
|
|
//=============================================================================
|
|
|
|
//
|
|
// ---------------------------------------------------
|
|
// | | | | | | |
|
|
// |cube1 | cube5 |cube2 |cube3 | cube5 |cube4 |
|
|
// | ext1 | cA | ext2 | ext3 | cD | ext4 |
|
|
// | (x1) | (xa) | (x2) | (x3) | (xd) | (x4) |
|
|
// | | | | | | |
|
|
// ----------------------------------------------------
|
|
// The following 3 classes perform the same functionaliry in 3 directions
|
|
class IndexTranslator6CubesLeftRight
|
|
{
|
|
public:
|
|
IndexTranslator6CubesLeftRight(
|
|
vtkm::Id x_1, vtkm::Id y_1, vtkm::Id z_1,
|
|
vtkm::Id x_2, vtkm::Id y_2, vtkm::Id z_2,
|
|
vtkm::Id x_3, vtkm::Id y_3, vtkm::Id z_3,
|
|
vtkm::Id x_4, vtkm::Id y_4, vtkm::Id z_4,
|
|
vtkm::Id x_a, vtkm::Id y_a, vtkm::Id z_a,
|
|
vtkm::Id x_d, vtkm::Id y_d, vtkm::Id z_d,
|
|
vtkm::Id x_5, vtkm::Id y_5, vtkm::Id z_5,
|
|
vtkm::Id start_x5, vtkm::Id start_y5, vtkm::Id start_z5 )
|
|
:
|
|
dimX1(x_1), dimY1(y_1), dimZ1(z_1),
|
|
dimX2(x_2), dimY2(y_2), dimZ2(z_2),
|
|
dimX3(x_3), dimY3(y_3), dimZ3(z_3),
|
|
dimX4(x_4), dimY4(y_4), dimZ4(z_4),
|
|
dimXa(x_a), dimYa(y_a), dimZa(z_a),
|
|
dimXd(x_d), dimYd(y_d), dimZd(z_d),
|
|
dimX5(x_5), dimY5(y_5), dimZ5(z_5),
|
|
startX5(start_x5), startY5(start_y5), startZ5(start_z5)
|
|
{ }
|
|
|
|
VTKM_EXEC_CONT
|
|
void Translate3Dto1D( vtkm::Id inX, vtkm::Id inY, vtkm::Id inZ, // 3D indices as input
|
|
vtkm::Id &cube, vtkm::Id &idx ) const // which cube, and idx of that cube
|
|
{
|
|
if ( dimX1 <= inX && inX < (dimX1 + dimXa) )
|
|
{
|
|
vtkm::Id inX_local = inX - dimX1;
|
|
cube = 5; // cAcD
|
|
idx = (inZ + startZ5) * dimX5 * dimY5 + (inY + startY5) * dimX5 + (inX_local + startX5 );
|
|
}
|
|
else if ( (dimX1 + dimXa + dimX2 + dimX3) <= inX && inX < (dimX1 + dimXa + dimX2 + dimX3 + dimXd) )
|
|
{
|
|
vtkm::Id inX_local = inX - dimX1 - dimX2 - dimX3; // no -dimXa since this is on the same cube
|
|
cube = 5; // cAcD
|
|
idx = (inZ + startZ5) * dimX5 * dimY5 + (inY + startY5) * dimX5 + (inX_local + startX5 );
|
|
}
|
|
else if ( 0 <= inX && inX < dimX1 )
|
|
{
|
|
cube = 1; // ext1
|
|
idx = inZ * dimX1 * dimY1 + inY * dimX1 + inX;
|
|
}
|
|
else if ( (dimX1 + dimXa) <= inX && inX < (dimX1 + dimXa + dimX2) )
|
|
{
|
|
vtkm::Id inX_local = inX - dimX1 - dimXa;
|
|
cube = 2; // ext2
|
|
idx = inZ * dimX2 * dimY2 + inY * dimX2 + inX_local;
|
|
}
|
|
else if ( (dimX1 + dimXa + dimX2) <= inX && inX < (dimX1 + dimXa + dimX2 + dimX3) )
|
|
{
|
|
vtkm::Id inX_local = inX - dimX1 - dimXa - dimX2;
|
|
cube = 3; // ext3
|
|
idx = inZ * dimX3 * dimY3 + inY * dimX3 + inX_local;
|
|
}
|
|
else if ( (dimX1 + dimXa + dimX2 + dimX3 + dimXd) <= inX &&
|
|
inX < (dimX1 + dimXa + dimX2 + dimX3 + dimXd + dimX4) )
|
|
{
|
|
vtkm::Id inX_local = inX - dimX1 - dimXa - dimX2 - dimX3 - dimXd;
|
|
cube = 4; // ext4
|
|
idx = inZ * dimX4 * dimY4 + inY * dimX4 + inX_local;
|
|
}
|
|
else
|
|
vtkm::cont::ErrorControlInternal("Invalid index!");
|
|
}
|
|
|
|
private:
|
|
const vtkm::Id dimX1, dimY1, dimZ1, dimX2, dimY2, dimZ2, dimX3, dimY3, dimZ3, dimX4, dimY4, dimZ4; // extension cube sizes
|
|
const vtkm::Id dimXa, dimYa, dimZa, dimXd, dimYd, dimZd; // signal cube sizes
|
|
const vtkm::Id dimX5, dimY5, dimZ5, startX5, startY5, startZ5; // entire cube size
|
|
};
|
|
|
|
class IndexTranslator6CubesTopDown
|
|
{
|
|
public:
|
|
IndexTranslator6CubesTopDown(
|
|
vtkm::Id x_1, vtkm::Id y_1, vtkm::Id z_1,
|
|
vtkm::Id x_2, vtkm::Id y_2, vtkm::Id z_2,
|
|
vtkm::Id x_3, vtkm::Id y_3, vtkm::Id z_3,
|
|
vtkm::Id x_4, vtkm::Id y_4, vtkm::Id z_4,
|
|
vtkm::Id x_a, vtkm::Id y_a, vtkm::Id z_a,
|
|
vtkm::Id x_d, vtkm::Id y_d, vtkm::Id z_d,
|
|
vtkm::Id x_5, vtkm::Id y_5, vtkm::Id z_5,
|
|
vtkm::Id start_x5, vtkm::Id start_y5, vtkm::Id start_z5 )
|
|
:
|
|
dimX1(x_1), dimY1(y_1), dimZ1(z_1),
|
|
dimX2(x_2), dimY2(y_2), dimZ2(z_2),
|
|
dimX3(x_3), dimY3(y_3), dimZ3(z_3),
|
|
dimX4(x_4), dimY4(y_4), dimZ4(z_4),
|
|
dimXa(x_a), dimYa(y_a), dimZa(z_a),
|
|
dimXd(x_d), dimYd(y_d), dimZd(z_d),
|
|
dimX5(x_5), dimY5(y_5), dimZ5(z_5),
|
|
startX5(start_x5), startY5(start_y5), startZ5(start_z5)
|
|
{ }
|
|
|
|
VTKM_EXEC_CONT
|
|
void Translate3Dto1D( vtkm::Id inX, vtkm::Id inY, vtkm::Id inZ, // 3D indices as input
|
|
vtkm::Id &cube, vtkm::Id &idx ) const // which cube, and idx of that cube
|
|
{
|
|
if ( dimY1 <= inY && inY < (dimY1 + dimYa) )
|
|
{
|
|
vtkm::Id inY_local = inY - dimY1;
|
|
cube = 5; // cAcD
|
|
idx = (inZ + startZ5) * dimX5 * dimY5 + (inY_local + startY5 ) * dimX5 + (inX + startX5);
|
|
}
|
|
else if ( (dimY1 + dimYa + dimY2 + dimY3) <= inY && inY < (dimY1 + dimYa + dimY2 + dimY3 + dimYd) )
|
|
{
|
|
vtkm::Id inY_local = inY - dimY1 - dimY2 - dimY3;
|
|
cube = 5; // cAcD
|
|
idx = (inZ + startZ5) * dimX5 * dimY5 + (inY_local + startY5) * dimX5 + (inX + startX5);
|
|
}
|
|
else if ( 0 <= inY && inY < dimY1 )
|
|
{
|
|
cube = 1; // ext1
|
|
idx = inZ * dimX1 * dimY1 + inY * dimX1 + inX;
|
|
}
|
|
else if ( (dimY1 + dimYa) <= inY && inY < (dimY1 + dimYa + dimY2) )
|
|
{
|
|
vtkm::Id inY_local = inY - dimY1 - dimYa;
|
|
cube = 2; // ext2
|
|
idx = inZ * dimX2 * dimY2 + inY_local * dimX2 + inX;
|
|
}
|
|
else if ( (dimY1 + dimYa + dimY2) <= inY && inY < (dimY1 + dimYa + dimY2 + dimY3) )
|
|
{
|
|
vtkm::Id inY_local = inY - dimY1 - dimYa - dimY2;
|
|
cube = 3; // ext3
|
|
idx = inZ * dimX3 * dimY3 + inY_local * dimX3 + inX;
|
|
}
|
|
else if ( (dimY1 + dimYa + dimY2 + dimY3 + dimYd) <= inY &&
|
|
inY < (dimY1 + dimYa + dimY2 + dimY3 + dimYd + dimY4) )
|
|
{
|
|
vtkm::Id inY_local = inY - dimY1 - dimYa - dimY2 - dimY3 - dimYd;
|
|
cube = 4; // ext4
|
|
idx = inZ * dimX4 * dimY4 + inY_local * dimX4 + inX;
|
|
}
|
|
else
|
|
vtkm::cont::ErrorControlInternal("Invalid index!");
|
|
}
|
|
|
|
private:
|
|
const vtkm::Id dimX1, dimY1, dimZ1, dimX2, dimY2, dimZ2, dimX3, dimY3, dimZ3, dimX4, dimY4, dimZ4; // extension cube sizes
|
|
const vtkm::Id dimXa, dimYa, dimZa, dimXd, dimYd, dimZd; // signal cube sizes
|
|
const vtkm::Id dimX5, dimY5, dimZ5, startX5, startY5, startZ5; // entire cube size
|
|
};
|
|
|
|
class IndexTranslator6CubesFrontBack
|
|
{
|
|
public:
|
|
VTKM_EXEC_CONT
|
|
IndexTranslator6CubesFrontBack(
|
|
vtkm::Id x_1, vtkm::Id y_1, vtkm::Id z_1,
|
|
vtkm::Id x_2, vtkm::Id y_2, vtkm::Id z_2,
|
|
vtkm::Id x_3, vtkm::Id y_3, vtkm::Id z_3,
|
|
vtkm::Id x_4, vtkm::Id y_4, vtkm::Id z_4,
|
|
vtkm::Id x_a, vtkm::Id y_a, vtkm::Id z_a,
|
|
vtkm::Id x_d, vtkm::Id y_d, vtkm::Id z_d,
|
|
vtkm::Id x_5, vtkm::Id y_5, vtkm::Id z_5,
|
|
vtkm::Id start_x5, vtkm::Id start_y5, vtkm::Id start_z5 )
|
|
:
|
|
dimX1(x_1), dimY1(y_1), dimZ1(z_1),
|
|
dimX2(x_2), dimY2(y_2), dimZ2(z_2),
|
|
dimX3(x_3), dimY3(y_3), dimZ3(z_3),
|
|
dimX4(x_4), dimY4(y_4), dimZ4(z_4),
|
|
dimXa(x_a), dimYa(y_a), dimZa(z_a),
|
|
dimXd(x_d), dimYd(y_d), dimZd(z_d),
|
|
dimX5(x_5), dimY5(y_5), dimZ5(z_5),
|
|
startX5(start_x5), startY5(start_y5), startZ5(start_z5)
|
|
{
|
|
/* printf("IndexTranslator6CubesFrontBack: \n" );
|
|
printf(" cube1 dims: (%lld, %lld, %lld)\n", dimX1, dimY1, dimZ1 );
|
|
printf(" cube2 dims: (%lld, %lld, %lld)\n", dimX2, dimY2, dimZ2 );
|
|
printf(" cube3 dims: (%lld, %lld, %lld)\n", dimX3, dimY3, dimZ3 );
|
|
printf(" cube4 dims: (%lld, %lld, %lld)\n", dimX4, dimY4, dimZ4 );
|
|
printf(" cube5 dims: (%lld, %lld, %lld)\n", dimX5, dimY5, dimZ5 );
|
|
printf(" cA dims: (%lld, %lld, %lld)\n", dimXa, dimYa, dimZa );
|
|
printf(" cD dims: (%lld, %lld, %lld)\n", dimXd, dimYd, dimZd );
|
|
printf(" start idx : (%lld, %lld, %lld)\n", startX5, startY5, startZ5 ); */
|
|
}
|
|
|
|
VTKM_EXEC_CONT
|
|
void Translate3Dto1D( vtkm::Id inX, vtkm::Id inY, vtkm::Id inZ, // 3D indices as input
|
|
vtkm::Id &cube, vtkm::Id &idx ) const // which cube, and idx of that cube
|
|
{
|
|
if ( dimZ1 <= inZ && inZ < (dimZ1 + dimZa) )
|
|
{
|
|
vtkm::Id inZ_local = inZ - dimZ1;
|
|
cube = 5; // cAcD
|
|
idx = (inZ_local + startZ5) * dimX5 * dimY5 + (inY + startY5 ) * dimX5 + (inX + startX5);
|
|
}
|
|
else if ( (dimZ1 + dimZa + dimZ2 + dimZ3) <= inZ && inZ < (dimZ1 + dimZa + dimZ2 + dimZ3 + dimZd) )
|
|
{
|
|
vtkm::Id inZ_local = inZ - dimZ1 - dimZ2 - dimZ3;
|
|
cube = 5; // cAcD
|
|
idx = (inZ_local + startZ5) * dimX5 * dimY5 + (inY + startY5) * dimX5 + (inX + startX5);
|
|
}
|
|
else if ( 0 <= inZ && inZ < dimZ1 )
|
|
{
|
|
cube = 1; // ext1
|
|
idx = inZ * dimX1 * dimY1 + inY * dimX1 + inX;
|
|
}
|
|
else if ( (dimZ1 + dimZa) <= inZ && inZ < (dimZ1 + dimZa + dimZ2) )
|
|
{
|
|
vtkm::Id inZ_local = inZ - dimZ1 - dimZa;
|
|
cube = 2; // ext2
|
|
idx = inZ_local * dimX2 * dimY2 + inY * dimX2 + inX;
|
|
}
|
|
else if ( (dimZ1 + dimZa + dimZ2) <= inZ && inZ < (dimZ1 + dimZa + dimZ2 + dimZ3) )
|
|
{
|
|
vtkm::Id inZ_local = inZ - dimZ1 - dimZa - dimZ2;
|
|
cube = 3; // ext3
|
|
idx = inZ_local * dimX3 * dimY3 + inY * dimX3 + inX;
|
|
}
|
|
else if ( (dimZ1 + dimZa + dimZ2 + dimZ3 + dimZd) <= inZ &&
|
|
inZ < (dimZ1 + dimZa + dimZ2 + dimZ3 + dimZd + dimZ4) )
|
|
{
|
|
vtkm::Id inZ_local = inZ - dimZ1 - dimZa - dimZ2 - dimZ3 - dimZd;
|
|
cube = 4; // ext4
|
|
idx = inZ_local * dimX4 * dimY4 + inY * dimX4 + inX;
|
|
}
|
|
else
|
|
vtkm::cont::ErrorControlInternal("Invalid index!");
|
|
}
|
|
|
|
private:
|
|
const vtkm::Id dimX1, dimY1, dimZ1, dimX2, dimY2, dimZ2, dimX3, dimY3, dimZ3, dimX4, dimY4, dimZ4; // extension cube sizes
|
|
const vtkm::Id dimXa, dimYa, dimZa, dimXd, dimYd, dimZd; // signal cube sizes
|
|
const vtkm::Id dimX5, dimY5, dimZ5, startX5, startY5, startZ5; // entire cube size
|
|
};
|
|
|
|
|
|
|
|
//=============================================================================
|
|
|
|
// Worklet: 3D forward transform along X (left-right)
|
|
template< typename DeviceTag >
|
|
class ForwardTransform3DLeftRight: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature(WholeArrayIn<ScalarAll>, // left extension
|
|
WholeArrayIn<ScalarAll>, // signal
|
|
WholeArrayIn<ScalarAll>, // right extension
|
|
WholeArrayOut<ScalarAll>); // cA followed by cD
|
|
typedef void ExecutionSignature(_1, _2, _3, _4, WorkIndex);
|
|
typedef _4 InputDomain;
|
|
|
|
VTKM_EXEC_CONT
|
|
ForwardTransform3DLeftRight (
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &loFilter,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &hiFilter,
|
|
vtkm::Id filter_len,
|
|
vtkm::Id approx_len,
|
|
bool odd_low,
|
|
vtkm::Id dimX1, vtkm::Id dimY1, vtkm::Id dimZ1,
|
|
vtkm::Id dimX2, vtkm::Id dimY2, vtkm::Id dimZ2,
|
|
vtkm::Id startX2, vtkm::Id startY2, vtkm::Id startZ2,
|
|
vtkm::Id pretendX2, vtkm::Id pretendY2, vtkm::Id pretendZ2,
|
|
vtkm::Id dimX3, vtkm::Id dimY3, vtkm::Id dimZ3 )
|
|
:
|
|
lowFilter( loFilter.PrepareForInput( DeviceTag() ) ),
|
|
highFilter( hiFilter.PrepareForInput( DeviceTag() ) ),
|
|
filterLen( filter_len ),
|
|
approxLen( approx_len ),
|
|
outDimX( pretendX2 ), outDimY( pretendY2 ), outDimZ( pretendZ2 ),
|
|
translator( dimX1, dimY1, dimZ1,
|
|
dimX2, dimY2, dimZ2,
|
|
startX2, startY2, startZ2,
|
|
pretendX2, pretendY2, pretendZ2,
|
|
dimX3, dimY3, dimZ3 )
|
|
{
|
|
this->lstart = odd_low ? 1 : 0;
|
|
this->hstart = 1;
|
|
}
|
|
|
|
VTKM_EXEC_CONT
|
|
void Output1Dto3D( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y, vtkm::Id &z ) const
|
|
{
|
|
z = idx / (outDimX * outDimY);
|
|
y = (idx - z * outDimX * outDimY) / outDimX;
|
|
x = idx % outDimX;
|
|
}
|
|
VTKM_EXEC_CONT
|
|
vtkm::Id Output3Dto1D( vtkm::Id x, vtkm::Id y, vtkm::Id z ) const
|
|
{
|
|
return z * outDimX * outDimY + y * outDimX + x;
|
|
}
|
|
|
|
// Use 64-bit float for convolution calculation
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
template <typename InPortalType1, typename InPortalType2, typename InPortalType3 >
|
|
VTKM_EXEC_CONT
|
|
VAL GetVal( const InPortalType1 &portal1, const InPortalType2 &portal2,
|
|
const InPortalType3 &portal3, vtkm::Id inCube, vtkm::Id inIdx ) const
|
|
{
|
|
if( inCube == 2 )
|
|
return MAKEVAL( portal2.Get(inIdx) );
|
|
else if( inCube == 1 )
|
|
return MAKEVAL( portal1.Get(inIdx) );
|
|
else if( inCube == 3 )
|
|
return MAKEVAL( portal3.Get(inIdx) );
|
|
else
|
|
{
|
|
vtkm::cont::ErrorControlInternal("Invalid cube index!");
|
|
return -1;
|
|
}
|
|
}
|
|
|
|
template <typename InPortalType1, typename InPortalType2,
|
|
typename InPortalType3, typename OutputPortalType>
|
|
VTKM_EXEC_CONT
|
|
void operator()(const InPortalType1 &inPortal1, // left extension
|
|
const InPortalType2 &inPortal2, // signal
|
|
const InPortalType3 &inPortal3, // right extension
|
|
OutputPortalType &coeffOut,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
vtkm::Id workX, workY, workZ, output1D;
|
|
Output1Dto3D( workIndex, workX, workY, workZ );
|
|
vtkm::Id inputCube, inputIdx;
|
|
typedef typename OutputPortalType::ValueType OutputValueType;
|
|
|
|
if( workX % 2 == 0 ) // calculate cA
|
|
{
|
|
vtkm::Id xl = lstart + workX;
|
|
VAL sum = MAKEVAL(0.0);
|
|
for( vtkm::Id k = filterLen - 1; k > -1; k-- )
|
|
{
|
|
translator.Translate3Dto1D( xl, workY, workZ, inputCube, inputIdx );
|
|
sum += lowFilter.Get(k) *
|
|
GetVal( inPortal1, inPortal2, inPortal3, inputCube, inputIdx );
|
|
xl++;
|
|
}
|
|
output1D = Output3Dto1D( workX/2, workY, workZ );
|
|
coeffOut.Set( output1D, static_cast<OutputValueType>(sum) );
|
|
}
|
|
else // calculate cD
|
|
{
|
|
vtkm::Id xh = hstart + workX - 1;
|
|
VAL sum=MAKEVAL(0.0);
|
|
for( vtkm::Id k = filterLen - 1; k > -1; k-- )
|
|
{
|
|
translator.Translate3Dto1D( xh, workY, workZ, inputCube, inputIdx );
|
|
sum += highFilter.Get(k) *
|
|
GetVal( inPortal1, inPortal2, inPortal3, inputCube, inputIdx );
|
|
xh++;
|
|
}
|
|
output1D = Output3Dto1D( (workX-1)/2 + approxLen, workY, workZ );
|
|
coeffOut.Set( output1D, static_cast<OutputValueType>(sum) );
|
|
}
|
|
}
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
private:
|
|
const typename vtkm::cont::ArrayHandle<vtkm::Float64>::ExecutionTypes<DeviceTag>::PortalConst
|
|
lowFilter, highFilter;
|
|
const vtkm::Id filterLen, approxLen;
|
|
const vtkm::Id outDimX, outDimY, outDimZ;
|
|
const IndexTranslator3CubesLeftRight translator;
|
|
vtkm::Id lstart, hstart;
|
|
};
|
|
|
|
template< typename DeviceTag >
|
|
class ForwardTransform3DTopDown: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature(WholeArrayIn<ScalarAll>, // left extension
|
|
WholeArrayIn<ScalarAll>, // signal
|
|
WholeArrayIn<ScalarAll>, // right extension
|
|
WholeArrayOut<ScalarAll>); // cA followed by cD
|
|
typedef void ExecutionSignature(_1, _2, _3, _4, WorkIndex);
|
|
typedef _4 InputDomain;
|
|
|
|
VTKM_EXEC_CONT
|
|
ForwardTransform3DTopDown (
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &loFilter,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &hiFilter,
|
|
vtkm::Id filter_len,
|
|
vtkm::Id approx_len,
|
|
bool odd_low,
|
|
vtkm::Id dimX1, vtkm::Id dimY1, vtkm::Id dimZ1,
|
|
vtkm::Id dimX2, vtkm::Id dimY2, vtkm::Id dimZ2,
|
|
vtkm::Id startX2, vtkm::Id startY2, vtkm::Id startZ2,
|
|
vtkm::Id pretendX2, vtkm::Id pretendY2, vtkm::Id pretendZ2,
|
|
vtkm::Id dimX3, vtkm::Id dimY3, vtkm::Id dimZ3 )
|
|
:
|
|
lowFilter( loFilter.PrepareForInput( DeviceTag() ) ),
|
|
highFilter( hiFilter.PrepareForInput( DeviceTag() ) ),
|
|
filterLen( filter_len ),
|
|
approxLen( approx_len ),
|
|
outDimX( pretendX2 ), outDimY( pretendY2 ), outDimZ( pretendZ2 ),
|
|
translator( dimX1, dimY1, dimZ1,
|
|
dimX2, dimY2, dimZ2,
|
|
startX2, startY2, startZ2,
|
|
pretendX2, pretendY2, pretendZ2,
|
|
dimX3, dimY3, dimZ3 )
|
|
{
|
|
this->lstart = odd_low ? 1 : 0;
|
|
this->hstart = 1;
|
|
}
|
|
|
|
VTKM_EXEC_CONT
|
|
void Output1Dto3D( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y, vtkm::Id &z ) const
|
|
{
|
|
z = idx / (outDimX * outDimY);
|
|
y = (idx - z * outDimX * outDimY) / outDimX;
|
|
x = idx % outDimX;
|
|
}
|
|
VTKM_EXEC_CONT
|
|
vtkm::Id Output3Dto1D( vtkm::Id x, vtkm::Id y, vtkm::Id z ) const
|
|
{
|
|
return z * outDimX * outDimY + y * outDimX + x;
|
|
}
|
|
|
|
// Use 64-bit float for convolution calculation
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
template <typename InPortalType1, typename InPortalType2, typename InPortalType3 >
|
|
VTKM_EXEC_CONT
|
|
VAL GetVal( const InPortalType1 &portal1, const InPortalType2 &portal2,
|
|
const InPortalType3 &portal3, vtkm::Id inCube, vtkm::Id inIdx ) const
|
|
{
|
|
if( inCube == 2 )
|
|
return MAKEVAL( portal2.Get(inIdx) );
|
|
else if( inCube == 1 )
|
|
return MAKEVAL( portal1.Get(inIdx) );
|
|
else if( inCube == 3 )
|
|
return MAKEVAL( portal3.Get(inIdx) );
|
|
else
|
|
{
|
|
vtkm::cont::ErrorControlInternal("Invalid cube index!");
|
|
return -1;
|
|
}
|
|
}
|
|
|
|
template <typename InPortalType1, typename InPortalType2,
|
|
typename InPortalType3, typename OutputPortalType>
|
|
VTKM_EXEC_CONT
|
|
void operator()(const InPortalType1 &inPortal1, // top extension
|
|
const InPortalType2 &inPortal2, // signal
|
|
const InPortalType3 &inPortal3, // down extension
|
|
OutputPortalType &coeffOut,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
vtkm::Id workX, workY, workZ, output1D;
|
|
Output1Dto3D( workIndex, workX, workY, workZ );
|
|
vtkm::Id inputCube, inputIdx;
|
|
typedef typename OutputPortalType::ValueType OutputValueType;
|
|
|
|
if( workY % 2 == 0 ) // calculate cA
|
|
{
|
|
vtkm::Id yl = lstart + workY;
|
|
VAL sum = MAKEVAL(0.0);
|
|
for( vtkm::Id k = filterLen - 1; k > -1; k-- )
|
|
{
|
|
translator.Translate3Dto1D( workX, yl, workZ, inputCube, inputIdx );
|
|
sum += lowFilter.Get(k) *
|
|
GetVal( inPortal1, inPortal2, inPortal3, inputCube, inputIdx );
|
|
yl++;
|
|
}
|
|
output1D = Output3Dto1D( workX, workY/2, workZ );
|
|
coeffOut.Set( output1D, static_cast<OutputValueType>(sum) );
|
|
}
|
|
else // calculate cD
|
|
{
|
|
vtkm::Id yh = hstart + workY - 1;
|
|
VAL sum=MAKEVAL(0.0);
|
|
for( vtkm::Id k = filterLen - 1; k > -1; k-- )
|
|
{
|
|
translator.Translate3Dto1D( workX, yh, workZ, inputCube, inputIdx );
|
|
sum += highFilter.Get(k) *
|
|
GetVal( inPortal1, inPortal2, inPortal3, inputCube, inputIdx );
|
|
yh++;
|
|
}
|
|
output1D = Output3Dto1D( workX, (workY-1)/2 + approxLen, workZ );
|
|
coeffOut.Set( output1D, static_cast<OutputValueType>(sum) );
|
|
}
|
|
}
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
private:
|
|
const typename vtkm::cont::ArrayHandle<vtkm::Float64>::ExecutionTypes<DeviceTag>::PortalConst
|
|
lowFilter, highFilter;
|
|
const vtkm::Id filterLen, approxLen;
|
|
const vtkm::Id outDimX, outDimY, outDimZ;
|
|
const IndexTranslator3CubesTopDown translator;
|
|
vtkm::Id lstart, hstart;
|
|
};
|
|
|
|
template< typename DeviceTag >
|
|
class ForwardTransform3DFrontBack: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature(WholeArrayIn<ScalarAll>, // left extension
|
|
WholeArrayIn<ScalarAll>, // signal
|
|
WholeArrayIn<ScalarAll>, // right extension
|
|
WholeArrayOut<ScalarAll>); // cA followed by cD
|
|
typedef void ExecutionSignature(_1, _2, _3, _4, WorkIndex);
|
|
typedef _4 InputDomain;
|
|
|
|
VTKM_EXEC_CONT
|
|
ForwardTransform3DFrontBack (
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &loFilter,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &hiFilter,
|
|
vtkm::Id filter_len,
|
|
vtkm::Id approx_len,
|
|
bool odd_low,
|
|
vtkm::Id dimX1, vtkm::Id dimY1, vtkm::Id dimZ1,
|
|
vtkm::Id dimX2, vtkm::Id dimY2, vtkm::Id dimZ2,
|
|
vtkm::Id startX2, vtkm::Id startY2, vtkm::Id startZ2,
|
|
vtkm::Id pretendX2, vtkm::Id pretendY2, vtkm::Id pretendZ2,
|
|
vtkm::Id dimX3, vtkm::Id dimY3, vtkm::Id dimZ3 )
|
|
:
|
|
lowFilter( loFilter.PrepareForInput( DeviceTag() ) ),
|
|
highFilter( hiFilter.PrepareForInput( DeviceTag() ) ),
|
|
filterLen( filter_len ),
|
|
approxLen( approx_len ),
|
|
outDimX( pretendX2 ), outDimY( pretendY2 ), outDimZ( pretendZ2 ),
|
|
translator( dimX1, dimY1, dimZ1,
|
|
dimX2, dimY2, dimZ2,
|
|
startX2, startY2, startZ2,
|
|
pretendX2, pretendY2, pretendZ2,
|
|
dimX3, dimY3, dimZ3 )
|
|
{
|
|
this->lstart = odd_low ? 1 : 0;
|
|
this->hstart = 1;
|
|
}
|
|
|
|
VTKM_EXEC_CONT
|
|
void Output1Dto3D( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y, vtkm::Id &z ) const
|
|
{
|
|
z = idx / (outDimX * outDimY);
|
|
y = (idx - z * outDimX * outDimY) / outDimX;
|
|
x = idx % outDimX;
|
|
}
|
|
VTKM_EXEC_CONT
|
|
vtkm::Id Output3Dto1D( vtkm::Id x, vtkm::Id y, vtkm::Id z ) const
|
|
{
|
|
return z * outDimX * outDimY + y * outDimX + x;
|
|
}
|
|
|
|
// Use 64-bit float for convolution calculation
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
template <typename InPortalType1, typename InPortalType2, typename InPortalType3 >
|
|
VTKM_EXEC_CONT
|
|
VAL GetVal( const InPortalType1 &portal1, const InPortalType2 &portal2,
|
|
const InPortalType3 &portal3, vtkm::Id inCube, vtkm::Id inIdx ) const
|
|
{
|
|
if( inCube == 2 )
|
|
return MAKEVAL( portal2.Get(inIdx) );
|
|
else if( inCube == 1 )
|
|
return MAKEVAL( portal1.Get(inIdx) );
|
|
else if( inCube == 3 )
|
|
return MAKEVAL( portal3.Get(inIdx) );
|
|
else
|
|
{
|
|
vtkm::cont::ErrorControlInternal("Invalid cube index!");
|
|
return -1;
|
|
}
|
|
}
|
|
|
|
template <typename InPortalType1, typename InPortalType2,
|
|
typename InPortalType3, typename OutputPortalType>
|
|
VTKM_EXEC_CONT
|
|
void operator()(const InPortalType1 &inPortal1, // top extension
|
|
const InPortalType2 &inPortal2, // signal
|
|
const InPortalType3 &inPortal3, // down extension
|
|
OutputPortalType &coeffOut,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
vtkm::Id workX, workY, workZ, output1D;
|
|
Output1Dto3D( workIndex, workX, workY, workZ );
|
|
vtkm::Id inputCube, inputIdx;
|
|
typedef typename OutputPortalType::ValueType OutputValueType;
|
|
|
|
if( workZ % 2 == 0 ) // calculate cA
|
|
{
|
|
vtkm::Id zl = lstart + workZ;
|
|
VAL sum = MAKEVAL(0.0);
|
|
for( vtkm::Id k = filterLen - 1; k > -1; k-- )
|
|
{
|
|
translator.Translate3Dto1D( workX, workY, zl, inputCube, inputIdx );
|
|
sum += lowFilter.Get(k) *
|
|
GetVal( inPortal1, inPortal2, inPortal3, inputCube, inputIdx );
|
|
zl++;
|
|
}
|
|
output1D = Output3Dto1D( workX, workY, workZ/2 );
|
|
coeffOut.Set( output1D, static_cast<OutputValueType>(sum) );
|
|
}
|
|
else // calculate cD
|
|
{
|
|
vtkm::Id zh = hstart + workZ - 1;
|
|
VAL sum=MAKEVAL(0.0);
|
|
for( vtkm::Id k = filterLen - 1; k > -1; k-- )
|
|
{
|
|
translator.Translate3Dto1D( workX, workY, zh, inputCube, inputIdx );
|
|
sum += highFilter.Get(k) *
|
|
GetVal( inPortal1, inPortal2, inPortal3, inputCube, inputIdx );
|
|
zh++;
|
|
}
|
|
output1D = Output3Dto1D( workX, workY, (workZ-1)/2 + approxLen );
|
|
coeffOut.Set( output1D, static_cast<OutputValueType>(sum) );
|
|
}
|
|
}
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
private:
|
|
const typename vtkm::cont::ArrayHandle<vtkm::Float64>::ExecutionTypes<DeviceTag>::PortalConst
|
|
lowFilter, highFilter;
|
|
const vtkm::Id filterLen, approxLen;
|
|
const vtkm::Id outDimX, outDimY, outDimZ;
|
|
const IndexTranslator3CubesFrontBack translator;
|
|
vtkm::Id lstart, hstart;
|
|
};
|
|
|
|
|
|
|
|
//=============================================================================
|
|
|
|
//
|
|
// ---------------------------------------------------
|
|
// | | | | | | |
|
|
// |cube1 | cube5 |cube2 |cube3 | cube5 |cube4 |
|
|
// | ext1 | cA | ext2 | ext3 | cD | ext4 |
|
|
// | (x1) | (xa) | (x2) | (x3) | (xd) | (x4) |
|
|
// | | | | | | |
|
|
// ----------------------------------------------------
|
|
// The following 3 classes perform the same functionaliry in 3 directions
|
|
//
|
|
// Worklet: perform a simple 3D inverse transform along X direction (Left-right)
|
|
template< typename DeviceTag >
|
|
class InverseTransform3DLeftRight: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayIn< ScalarAll >, // ext1
|
|
WholeArrayIn< ScalarAll >, // ext2
|
|
WholeArrayIn< ScalarAll >, // ext3
|
|
WholeArrayIn< ScalarAll >, // ext4
|
|
WholeArrayIn< ScalarAll >, // cA+cD (signal)
|
|
FieldOut< ScalarAll> ); // outptu coefficients
|
|
typedef void ExecutionSignature( _1, _2, _3, _4, _5, _6, WorkIndex );
|
|
typedef _6 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
InverseTransform3DLeftRight(
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &lo_fil,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &hi_fil,
|
|
vtkm::Id fil_len,
|
|
vtkm::Id x_1, vtkm::Id y_1, vtkm::Id z_1, // ext1
|
|
vtkm::Id x_2, vtkm::Id y_2, vtkm::Id z_2, // ext2
|
|
vtkm::Id x_3, vtkm::Id y_3, vtkm::Id z_3, // ext3
|
|
vtkm::Id x_4, vtkm::Id y_4, vtkm::Id z_4, // ext4
|
|
vtkm::Id x_a, vtkm::Id y_a, vtkm::Id z_a, // cA
|
|
vtkm::Id x_d, vtkm::Id y_d, vtkm::Id z_d, // cD
|
|
vtkm::Id x_5, vtkm::Id y_5, vtkm::Id z_5, // signal, actual dims
|
|
vtkm::Id startX5, vtkm::Id startY5, vtkm::Id startZ5 )
|
|
:
|
|
lowFilter( lo_fil.PrepareForInput( DeviceTag() ) ),
|
|
highFilter( hi_fil.PrepareForInput( DeviceTag() ) ),
|
|
filterLen( fil_len ),
|
|
outDimX(x_a + x_d), outDimY( y_a ), outDimZ( z_a ),
|
|
cALenExtended( x_1 + x_a + x_2 ),
|
|
translator(x_1, y_1, z_1,
|
|
x_2, y_2, z_2,
|
|
x_3, y_3, z_3,
|
|
x_4, y_4, z_4,
|
|
x_a, y_a, z_a,
|
|
x_d, y_d, z_d,
|
|
x_5, y_5, z_5,
|
|
startX5, startY5, startZ5 )
|
|
{ }
|
|
|
|
VTKM_EXEC_CONT
|
|
void Output1Dto3D( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y, vtkm::Id &z ) const
|
|
{
|
|
z = idx / (outDimX * outDimY);
|
|
y = (idx - z * outDimX * outDimY) / outDimX;
|
|
x = idx % outDimX;
|
|
}
|
|
|
|
// Use 64-bit float for convolution calculation
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
template <typename InPortalType1, typename InPortalType2, typename InPortalType3,
|
|
typename InPortalType4, typename InPortalType5 >
|
|
VTKM_EXEC_CONT
|
|
VAL GetVal( const InPortalType1 &ext1,
|
|
const InPortalType2 &ext2,
|
|
const InPortalType3 &ext3,
|
|
const InPortalType4 &ext4,
|
|
const InPortalType5 &sig5,
|
|
vtkm::Id inCube,
|
|
vtkm::Id inIdx ) const
|
|
{
|
|
if( inCube == 2 )
|
|
return MAKEVAL( ext2.Get(inIdx) );
|
|
else if( inCube == 4 )
|
|
return MAKEVAL( ext4.Get(inIdx) );
|
|
else if( inCube == 1 )
|
|
return MAKEVAL( ext1.Get(inIdx) );
|
|
else if( inCube == 3 )
|
|
return MAKEVAL( ext3.Get(inIdx) );
|
|
else if( inCube == 5 )
|
|
return MAKEVAL( sig5.Get(inIdx) );
|
|
else
|
|
{
|
|
vtkm::cont::ErrorControlInternal("Invalid matrix index!");
|
|
return -1;
|
|
}
|
|
}
|
|
|
|
template< typename InPortalType1, typename InPortalType2, typename InPortalType3,
|
|
typename InPortalType4, typename InPortalType5, typename OutputValueType >
|
|
VTKM_EXEC
|
|
void operator() (const InPortalType1 &portal1,
|
|
const InPortalType2 &portal2,
|
|
const InPortalType3 &portal3,
|
|
const InPortalType4 &portal4,
|
|
const InPortalType5 &portal5,
|
|
OutputValueType &coeffOut,
|
|
const vtkm::Id &workIdx ) const
|
|
{
|
|
vtkm::Id workX, workY, workZ;
|
|
vtkm::Id k1, k2, xi;
|
|
vtkm::Id inputCube, inputIdx;
|
|
Output1Dto3D( workIdx, workX, workY, workZ );
|
|
|
|
if( filterLen % 2 != 0 ) // odd filter
|
|
{
|
|
if( workX % 2 != 0 )
|
|
{
|
|
k1 = filterLen - 2; k2 = filterLen - 1;
|
|
}
|
|
else
|
|
{
|
|
k1 = filterLen - 1; k2 = filterLen - 2;
|
|
}
|
|
|
|
VAL sum = 0.0;
|
|
xi = (workX + 1) / 2;
|
|
while( k1 > -1 )
|
|
{
|
|
translator.Translate3Dto1D( xi, workY, workZ, inputCube, inputIdx );
|
|
sum += lowFilter.Get(k1) * GetVal( portal1, portal2, portal3, portal4, portal5,
|
|
inputCube, inputIdx );
|
|
xi++;
|
|
k1 -= 2;
|
|
}
|
|
xi = workX / 2;
|
|
while( k2 > -1 )
|
|
{
|
|
translator.Translate3Dto1D( xi + cALenExtended, workY, workZ, inputCube, inputIdx );
|
|
sum += highFilter.Get(k2) * GetVal( portal1, portal2, portal3, portal4, portal5,
|
|
inputCube, inputIdx );
|
|
xi++;
|
|
k2 -= 2;
|
|
}
|
|
coeffOut = static_cast< OutputValueType> (sum);
|
|
}
|
|
else // even filter
|
|
{
|
|
if( (filterLen/2) % 2 != 0 ) // odd length half filter
|
|
{
|
|
xi = workX / 2;
|
|
if( workX % 2 != 0 )
|
|
k1 = filterLen - 1;
|
|
else
|
|
k1 = filterLen - 2;
|
|
}
|
|
else // even length half filter
|
|
{
|
|
xi = (workX + 1) / 2;
|
|
if( workX % 2 != 0 )
|
|
k1 = filterLen - 2;
|
|
else
|
|
k1 = filterLen - 1;
|
|
}
|
|
VAL cA, cD;
|
|
VAL sum = 0.0;
|
|
while( k1 > -1 )
|
|
{
|
|
translator.Translate3Dto1D( xi, workY, workZ, inputCube, inputIdx );
|
|
cA = GetVal( portal1, portal2, portal3, portal4, portal5, inputCube, inputIdx );
|
|
translator.Translate3Dto1D( xi + cALenExtended, workY, workZ, inputCube, inputIdx );
|
|
cD = GetVal( portal1, portal2, portal3, portal4, portal5, inputCube, inputIdx );
|
|
sum += lowFilter.Get(k1) * cA + highFilter.Get(k1) * cD;
|
|
xi++;
|
|
k1 -= 2;
|
|
}
|
|
coeffOut = static_cast< OutputValueType >(sum);
|
|
}
|
|
}
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
private:
|
|
const typename vtkm::cont::ArrayHandle<vtkm::Float64>::ExecutionTypes<DeviceTag>::PortalConst
|
|
lowFilter, highFilter;
|
|
const vtkm::Id filterLen;
|
|
vtkm::Id outDimX, outDimY, outDimZ;
|
|
vtkm::Id cALenExtended; // Number of cA at the beginning of input, followed by cD
|
|
const IndexTranslator6CubesLeftRight translator;
|
|
};
|
|
|
|
template< typename DeviceTag >
|
|
class InverseTransform3DTopDown: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayIn< ScalarAll >, // ext1
|
|
WholeArrayIn< ScalarAll >, // ext2
|
|
WholeArrayIn< ScalarAll >, // ext3
|
|
WholeArrayIn< ScalarAll >, // ext4
|
|
WholeArrayIn< ScalarAll >, // cA+cD (signal)
|
|
FieldOut< ScalarAll> ); // outptu coefficients
|
|
typedef void ExecutionSignature( _1, _2, _3, _4, _5, _6, WorkIndex );
|
|
typedef _6 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
InverseTransform3DTopDown(
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &lo_fil,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &hi_fil,
|
|
vtkm::Id fil_len,
|
|
vtkm::Id x_1, vtkm::Id y_1, vtkm::Id z_1, // ext1
|
|
vtkm::Id x_2, vtkm::Id y_2, vtkm::Id z_2, // ext2
|
|
vtkm::Id x_3, vtkm::Id y_3, vtkm::Id z_3, // ext3
|
|
vtkm::Id x_4, vtkm::Id y_4, vtkm::Id z_4, // ext4
|
|
vtkm::Id x_a, vtkm::Id y_a, vtkm::Id z_a, // cA
|
|
vtkm::Id x_d, vtkm::Id y_d, vtkm::Id z_d, // cD
|
|
vtkm::Id x_5, vtkm::Id y_5, vtkm::Id z_5, // signal, actual dims
|
|
vtkm::Id startX5, vtkm::Id startY5, vtkm::Id startZ5 )
|
|
:
|
|
lowFilter( lo_fil.PrepareForInput( DeviceTag() ) ),
|
|
highFilter( hi_fil.PrepareForInput( DeviceTag() ) ),
|
|
filterLen( fil_len ),
|
|
outDimX( x_a ), outDimY(y_a + y_d), outDimZ( z_a ),
|
|
cALenExtended( y_1 + y_a + y_2 ),
|
|
translator(x_1, y_1, z_1,
|
|
x_2, y_2, z_2,
|
|
x_3, y_3, z_3,
|
|
x_4, y_4, z_4,
|
|
x_a, y_a, z_a,
|
|
x_d, y_d, z_d,
|
|
x_5, y_5, z_5,
|
|
startX5, startY5, startZ5 )
|
|
{ }
|
|
|
|
VTKM_EXEC_CONT
|
|
void Output1Dto3D( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y, vtkm::Id &z ) const
|
|
{
|
|
z = idx / (outDimX * outDimY);
|
|
y = (idx - z * outDimX * outDimY) / outDimX;
|
|
x = idx % outDimX;
|
|
}
|
|
|
|
// Use 64-bit float for convolution calculation
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
template <typename InPortalType1, typename InPortalType2, typename InPortalType3,
|
|
typename InPortalType4, typename InPortalType5 >
|
|
VTKM_EXEC_CONT
|
|
VAL GetVal( const InPortalType1 &ext1,
|
|
const InPortalType2 &ext2,
|
|
const InPortalType3 &ext3,
|
|
const InPortalType4 &ext4,
|
|
const InPortalType5 &sig5,
|
|
vtkm::Id inCube,
|
|
vtkm::Id inIdx ) const
|
|
{
|
|
if( inCube == 2 )
|
|
return MAKEVAL( ext2.Get(inIdx) );
|
|
else if( inCube == 4 )
|
|
return MAKEVAL( ext4.Get(inIdx) );
|
|
else if( inCube == 1 )
|
|
return MAKEVAL( ext1.Get(inIdx) );
|
|
else if( inCube == 3 )
|
|
return MAKEVAL( ext3.Get(inIdx) );
|
|
else if( inCube == 5 )
|
|
return MAKEVAL( sig5.Get(inIdx) );
|
|
else
|
|
{
|
|
vtkm::cont::ErrorControlInternal("Invalid matrix index!");
|
|
return -1;
|
|
}
|
|
}
|
|
|
|
template< typename InPortalType1, typename InPortalType2, typename InPortalType3,
|
|
typename InPortalType4, typename InPortalType5, typename OutputValueType >
|
|
VTKM_EXEC
|
|
void operator() (const InPortalType1 &portal1,
|
|
const InPortalType2 &portal2,
|
|
const InPortalType3 &portal3,
|
|
const InPortalType4 &portal4,
|
|
const InPortalType5 &portal5,
|
|
OutputValueType &coeffOut,
|
|
const vtkm::Id &workIdx ) const
|
|
{
|
|
vtkm::Id workX, workY, workZ;
|
|
vtkm::Id k1, k2, yi;
|
|
vtkm::Id inputCube, inputIdx;
|
|
Output1Dto3D( workIdx, workX, workY, workZ );
|
|
|
|
if( filterLen % 2 != 0 ) // odd filter
|
|
{
|
|
if( workY % 2 != 0 )
|
|
{
|
|
k1 = filterLen - 2; k2 = filterLen - 1;
|
|
}
|
|
else
|
|
{
|
|
k1 = filterLen - 1; k2 = filterLen - 2;
|
|
}
|
|
|
|
VAL sum = 0.0;
|
|
yi = (workY + 1) / 2;
|
|
while( k1 > -1 )
|
|
{
|
|
translator.Translate3Dto1D( workX, yi, workZ, inputCube, inputIdx );
|
|
sum += lowFilter.Get(k1) * GetVal( portal1, portal2, portal3, portal4, portal5,
|
|
inputCube, inputIdx );
|
|
yi++;
|
|
k1 -= 2;
|
|
}
|
|
yi = workY / 2;
|
|
while( k2 > -1 )
|
|
{
|
|
translator.Translate3Dto1D( workX, yi + cALenExtended, workZ, inputCube, inputIdx );
|
|
sum += highFilter.Get(k2) * GetVal( portal1, portal2, portal3, portal4, portal5,
|
|
inputCube, inputIdx );
|
|
yi++;
|
|
k2 -= 2;
|
|
}
|
|
coeffOut = static_cast< OutputValueType >(sum);
|
|
}
|
|
else // even filter
|
|
{
|
|
if( (filterLen/2) % 2 != 0 )
|
|
{
|
|
yi = workY / 2;
|
|
if( workY % 2 != 0 )
|
|
k1 = filterLen - 1;
|
|
else
|
|
k1 = filterLen - 2;
|
|
}
|
|
else
|
|
{
|
|
yi = (workY + 1) / 2;
|
|
if( workY % 2 != 0 )
|
|
k1 = filterLen - 2;
|
|
else
|
|
k1 = filterLen - 1;
|
|
}
|
|
VAL cA, cD;
|
|
VAL sum = 0.0;
|
|
while( k1 > -1 )
|
|
{
|
|
translator.Translate3Dto1D( workX, yi, workZ, inputCube, inputIdx );
|
|
cA = GetVal( portal1, portal2, portal3, portal4, portal5, inputCube, inputIdx );
|
|
translator.Translate3Dto1D( workX, yi + cALenExtended, workZ, inputCube, inputIdx );
|
|
cD = GetVal( portal1, portal2, portal3, portal4, portal5, inputCube, inputIdx );
|
|
sum += lowFilter.Get(k1) * cA + highFilter.Get(k1) * cD;
|
|
yi++;
|
|
k1 -= 2;
|
|
}
|
|
coeffOut = static_cast< OutputValueType >(sum);
|
|
}
|
|
}
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
private:
|
|
const typename vtkm::cont::ArrayHandle<vtkm::Float64>::ExecutionTypes<DeviceTag>::PortalConst
|
|
lowFilter, highFilter;
|
|
const vtkm::Id filterLen;
|
|
vtkm::Id outDimX, outDimY, outDimZ;
|
|
vtkm::Id cALenExtended; // Number of cA at the beginning of input, followed by cD
|
|
const IndexTranslator6CubesTopDown translator;
|
|
};
|
|
|
|
template< typename DeviceTag >
|
|
class InverseTransform3DFrontBack: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayIn< ScalarAll >, // ext1
|
|
WholeArrayIn< ScalarAll >, // ext2
|
|
WholeArrayIn< ScalarAll >, // ext3
|
|
WholeArrayIn< ScalarAll >, // ext4
|
|
WholeArrayIn< ScalarAll >, // cA+cD (signal)
|
|
FieldOut< ScalarAll> ); // outptu coefficients
|
|
typedef void ExecutionSignature( _1, _2, _3, _4, _5, _6, WorkIndex );
|
|
typedef _6 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
InverseTransform3DFrontBack(
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &lo_fil,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &hi_fil,
|
|
vtkm::Id fil_len,
|
|
vtkm::Id x_1, vtkm::Id y_1, vtkm::Id z_1, // ext1
|
|
vtkm::Id x_2, vtkm::Id y_2, vtkm::Id z_2, // ext2
|
|
vtkm::Id x_3, vtkm::Id y_3, vtkm::Id z_3, // ext3
|
|
vtkm::Id x_4, vtkm::Id y_4, vtkm::Id z_4, // ext4
|
|
vtkm::Id x_a, vtkm::Id y_a, vtkm::Id z_a, // cA
|
|
vtkm::Id x_d, vtkm::Id y_d, vtkm::Id z_d, // cD
|
|
vtkm::Id x_5, vtkm::Id y_5, vtkm::Id z_5, // signal, actual dims
|
|
vtkm::Id startX5, vtkm::Id startY5, vtkm::Id startZ5 )
|
|
:
|
|
lowFilter( lo_fil.PrepareForInput( DeviceTag() ) ),
|
|
highFilter( hi_fil.PrepareForInput( DeviceTag() ) ),
|
|
filterLen( fil_len ),
|
|
outDimX( x_a ), outDimY( y_a ), outDimZ(z_a + z_d),
|
|
cALenExtended( z_1 + z_a + z_2 ),
|
|
translator(x_1, y_1, z_1,
|
|
x_2, y_2, z_2,
|
|
x_3, y_3, z_3,
|
|
x_4, y_4, z_4,
|
|
x_a, y_a, z_a,
|
|
x_d, y_d, z_d,
|
|
x_5, y_5, z_5,
|
|
startX5, startY5, startZ5 )
|
|
{
|
|
/* printf("InverseTransform3DFrontBack: \n");
|
|
printf(" output dims: (%lld, %lld, %lld)\n", outDimX, outDimY, outDimZ ); */
|
|
}
|
|
|
|
VTKM_EXEC_CONT
|
|
void Output1Dto3D( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y, vtkm::Id &z ) const
|
|
{
|
|
z = idx / (outDimX * outDimY);
|
|
y = (idx - z * outDimX * outDimY) / outDimX;
|
|
x = idx % outDimX;
|
|
}
|
|
|
|
// Use 64-bit float for convolution calculation
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
template <typename InPortalType1, typename InPortalType2, typename InPortalType3,
|
|
typename InPortalType4, typename InPortalType5 >
|
|
VTKM_EXEC_CONT
|
|
VAL GetVal( const InPortalType1 &ext1,
|
|
const InPortalType2 &ext2,
|
|
const InPortalType3 &ext3,
|
|
const InPortalType4 &ext4,
|
|
const InPortalType5 &sig5,
|
|
vtkm::Id inCube,
|
|
vtkm::Id inIdx ) const
|
|
{
|
|
if( inCube == 2 )
|
|
return MAKEVAL( ext2.Get(inIdx) );
|
|
else if( inCube == 4 )
|
|
return MAKEVAL( ext4.Get(inIdx) );
|
|
else if( inCube == 1 )
|
|
return MAKEVAL( ext1.Get(inIdx) );
|
|
else if( inCube == 3 )
|
|
return MAKEVAL( ext3.Get(inIdx) );
|
|
else if( inCube == 5 )
|
|
return MAKEVAL( sig5.Get(inIdx) );
|
|
else
|
|
{
|
|
vtkm::cont::ErrorControlInternal("Invalid matrix index!");
|
|
return -1;
|
|
}
|
|
}
|
|
|
|
template< typename InPortalType1, typename InPortalType2, typename InPortalType3,
|
|
typename InPortalType4, typename InPortalType5, typename OutputValueType >
|
|
VTKM_EXEC
|
|
void operator() (const InPortalType1 &portal1,
|
|
const InPortalType2 &portal2,
|
|
const InPortalType3 &portal3,
|
|
const InPortalType4 &portal4,
|
|
const InPortalType5 &portal5,
|
|
OutputValueType &coeffOut,
|
|
const vtkm::Id &workIdx ) const
|
|
{
|
|
vtkm::Id workX, workY, workZ;
|
|
vtkm::Id k1, k2, zi;
|
|
vtkm::Id inputCube, inputIdx;
|
|
Output1Dto3D( workIdx, workX, workY, workZ );
|
|
//printf("A new iteration: \n");
|
|
//printf(" workIdx = %lld, workX = %lld, workY = %lld, workZ = %lld\n", workIdx, workX, workY, workZ );
|
|
|
|
if( filterLen % 2 != 0 ) // odd filter
|
|
{
|
|
if( workZ % 2 != 0 )
|
|
{
|
|
k1 = filterLen - 2; k2 = filterLen - 1;
|
|
}
|
|
else
|
|
{
|
|
k1 = filterLen - 1; k2 = filterLen - 2;
|
|
}
|
|
|
|
VAL sum = 0.0;
|
|
zi = (workZ + 1) / 2;
|
|
while( k1 > -1 )
|
|
{
|
|
translator.Translate3Dto1D( workX, workY, zi, inputCube, inputIdx );
|
|
sum += lowFilter.Get(k1) * GetVal( portal1, portal2, portal3, portal4, portal5,
|
|
inputCube, inputIdx );
|
|
zi++;
|
|
k1 -= 2;
|
|
}
|
|
zi = workZ / 2;
|
|
while( k2 > -1 )
|
|
{
|
|
translator.Translate3Dto1D( workX, workY, zi + cALenExtended, inputCube, inputIdx );
|
|
sum += highFilter.Get(k2) * GetVal( portal1, portal2, portal3, portal4, portal5,
|
|
inputCube, inputIdx );
|
|
zi++;
|
|
k2 -= 2;
|
|
}
|
|
coeffOut = static_cast< OutputValueType >(sum);
|
|
}
|
|
else // even filter
|
|
{
|
|
if( (filterLen/2) % 2 != 0 )
|
|
{
|
|
zi = workZ / 2;
|
|
if( workZ % 2 != 0 )
|
|
k1 = filterLen - 1;
|
|
else
|
|
k1 = filterLen - 2;
|
|
}
|
|
else
|
|
{
|
|
zi = (workZ + 1) / 2;
|
|
if( workZ % 2 != 0 )
|
|
k1 = filterLen - 2;
|
|
else
|
|
k1 = filterLen - 1;
|
|
}
|
|
VAL cA, cD;
|
|
VAL sum = 0.0;
|
|
while( k1 > -1 )
|
|
{
|
|
translator.Translate3Dto1D( workX, workY, zi, inputCube, inputIdx );
|
|
cA = GetVal( portal1, portal2, portal3, portal4, portal5, inputCube, inputIdx );
|
|
translator.Translate3Dto1D( workX, workY, zi + cALenExtended, inputCube, inputIdx );
|
|
cD = GetVal( portal1, portal2, portal3, portal4, portal5, inputCube, inputIdx );
|
|
sum += lowFilter.Get(k1) * cA + highFilter.Get(k1) * cD;
|
|
zi++;
|
|
k1 -= 2;
|
|
}
|
|
coeffOut = static_cast< OutputValueType >(sum);
|
|
}
|
|
}
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
private:
|
|
const typename vtkm::cont::ArrayHandle<vtkm::Float64>::ExecutionTypes<DeviceTag>::PortalConst
|
|
lowFilter, highFilter;
|
|
const vtkm::Id filterLen;
|
|
vtkm::Id outDimX, outDimY, outDimZ;
|
|
vtkm::Id cALenExtended; // Number of cA at the beginning of input, followed by cD
|
|
const IndexTranslator6CubesFrontBack translator;
|
|
};
|
|
|
|
|
|
|
|
|
|
//=============================================================================
|
|
|
|
// ---------------------------------------------------
|
|
// | | | | | | |
|
|
// | | | | | | |
|
|
// | ext1 | cA | ext2 | ext3 | cD | ext4 |
|
|
// | (x1) | (xa) | (x2) | (x3) | (xd) | (x4) |
|
|
// | | | | | | |
|
|
// ----------------------------------------------------
|
|
// matrix1: ext1
|
|
// matrix2: ext2
|
|
// matrix3: ext3
|
|
// matrix4: ext4
|
|
// matrix5: cA + cD
|
|
class IndexTranslator6Matrices
|
|
{
|
|
public:
|
|
IndexTranslator6Matrices( vtkm::Id x_1, vtkm::Id y_1,
|
|
vtkm::Id x_a, vtkm::Id y_a,
|
|
vtkm::Id x_2, vtkm::Id y_2,
|
|
vtkm::Id x_3, vtkm::Id y_3,
|
|
vtkm::Id x_d, vtkm::Id y_d,
|
|
vtkm::Id x_4, vtkm::Id y_4,
|
|
vtkm::Id x_5, vtkm::Id y_5, // actual size of matrix5
|
|
vtkm::Id start_x5, vtkm::Id start_y5,// start indices of pretend matrix
|
|
bool mode )
|
|
: x1(x_1), y1(y_1),
|
|
xa(x_a), ya(y_a),
|
|
x2(x_2), y2(y_2),
|
|
x3(x_3), y3(y_3),
|
|
xd(x_d), yd(y_d),
|
|
x4(x_4), y4(y_4),
|
|
x5(x_5), y5(y_5),
|
|
startX5(start_x5), startY5(start_y5),
|
|
modeLR (mode)
|
|
{
|
|
// Get pretend matrix dims
|
|
if( modeLR )
|
|
{
|
|
pretendX5 = xa + xd;
|
|
pretendY5 = y1;
|
|
}
|
|
else
|
|
{
|
|
pretendX5 = x1;
|
|
pretendY5 = ya + yd;
|
|
}
|
|
(void)y5;
|
|
}
|
|
|
|
VTKM_EXEC_CONT
|
|
void Translate2Dto1D( vtkm::Id inX, vtkm::Id inY, // 2D indices as input
|
|
vtkm::Id &mat, vtkm::Id &idx ) const // which matrix, and idx of that matrix
|
|
{
|
|
if( modeLR ) // left-right mode
|
|
{
|
|
if ( 0 <= inX && inX < x1 )
|
|
{
|
|
mat = 1; // ext1
|
|
idx = inY * x1 + inX;
|
|
}
|
|
else if ( x1 <= inX && inX < (x1 + xa) )
|
|
{
|
|
mat = 5; // cAcD
|
|
idx = (inY + startY5) * x5 + (inX - x1 + startX5 );
|
|
}
|
|
else if ( (x1 + xa) <= inX && inX < (x1 + xa + x2) )
|
|
{
|
|
mat = 2; // ext2
|
|
idx = inY * x2 + (inX - x1 - xa);
|
|
}
|
|
else if ( (x1 + xa + x2) <= inX && inX < (x1 + xa + x2 + x3) )
|
|
{
|
|
mat = 3; // ext3
|
|
idx = inY * x3 + (inX - x1 - xa - x2);
|
|
}
|
|
else if ( (x1 + xa + x2 + x3) <= inX && inX < (x1 + xa + x2 + x3 + xd) )
|
|
{
|
|
mat = 5; // cAcD
|
|
idx = (inY + startY5) * x5 + (inX - x1 - x2 - x3 + startX5 );
|
|
}
|
|
else if ( (x1 + xa + x2 + x3 + xd) <= inX && inX < (x1 + xa + x2 + x3 + xd + x4) )
|
|
{
|
|
mat = 4; // ext4
|
|
idx = inY * x4 + (inX - x1 - xa - x2 - x3 - xd);
|
|
}
|
|
else
|
|
vtkm::cont::ErrorControlInternal("Invalid index!");
|
|
}
|
|
else // top-down mode
|
|
{
|
|
if ( 0 <= inY && inY < y1 )
|
|
{
|
|
mat = 1; // ext1
|
|
idx = inY * x1 + inX;
|
|
}
|
|
else if ( y1 <= inY && inY < (y1 + ya) )
|
|
{
|
|
mat = 5; // cAcD
|
|
idx = (inY - y1 + startY5 ) * x5 + inX + startX5;
|
|
}
|
|
else if ( (y1 + ya) <= inY && inY < (y1 + ya + y2) )
|
|
{
|
|
mat = 2; // ext2
|
|
idx = (inY - y1 - ya) * x2 + inX;
|
|
}
|
|
else if ( (y1 + ya + y2) <= inY && inY < (y1 + ya + y2 + y3) )
|
|
{
|
|
mat = 3; // ext3
|
|
idx = (inY - y1 - ya - y2) * x3 + inX;
|
|
}
|
|
else if ( (y1 + ya + y2 + y3) <= inY && inY < (y1 + ya + y2 + y3 + yd) )
|
|
{
|
|
mat = 5; // cAcD
|
|
idx = (inY - y1 - y2 - y3 + startY5 ) * x5 + inX + startX5;
|
|
}
|
|
else if ( (y1 + ya + y2 + y3 + yd) <= inY && inY < (y1 + ya + y2 + y3 + yd + y4) )
|
|
{
|
|
mat = 4; // ext4
|
|
idx = (inY - y1 - ya - y2 - y3 - yd) * x4 + inX;
|
|
}
|
|
else
|
|
vtkm::cont::ErrorControlInternal("Invalid index!");
|
|
}
|
|
}
|
|
|
|
private:
|
|
const vtkm::Id x1, y1, xa, ya, x2, y2, x3, y3, xd, yd, x4, y4;
|
|
vtkm::Id x5, y5, startX5, startY5, pretendX5, pretendY5;
|
|
const bool modeLR ; // true = left-right mode; false = top-down mode.
|
|
};
|
|
|
|
|
|
|
|
// ................
|
|
// . .
|
|
// -----.--------------.-----
|
|
// | . | | . |
|
|
// | . | | . |
|
|
// | ext1 | mat2 | ext2 |
|
|
// | (x1) | (x2) | (x3) |
|
|
// | . | | . |
|
|
// -----.--------------.-----
|
|
// ................
|
|
class IndexTranslator3Matrices
|
|
{
|
|
public:
|
|
IndexTranslator3Matrices( vtkm::Id x_1, vtkm::Id y_1,
|
|
vtkm::Id x_2, vtkm::Id y_2, // actual dims of mat2
|
|
vtkm::Id startx_2, vtkm::Id starty_2, // start idx of pretend
|
|
vtkm::Id pretendx_2, vtkm::Id pretendy_2,// pretend dims
|
|
vtkm::Id x_3, vtkm::Id y_3,
|
|
bool mode )
|
|
:
|
|
dimX1(x_1), dimY1(y_1),
|
|
dimX2(x_2), dimY2(y_2),
|
|
startX2( startx_2 ), startY2( starty_2 ),
|
|
pretendDimX2( pretendx_2 ), pretendDimY2( pretendy_2 ),
|
|
dimX3(x_3), dimY3(y_3),
|
|
mode_lr(mode)
|
|
{ (void)dimY2; }
|
|
|
|
VTKM_EXEC_CONT
|
|
void Translate2Dto1D( vtkm::Id inX, vtkm::Id inY, // 2D indices as input
|
|
vtkm::Id &mat, vtkm::Id &idx ) const // which matrix, and idx of that matrix
|
|
{
|
|
if( mode_lr ) // left-right mode
|
|
{
|
|
if ( 0 <= inX && inX < dimX1 )
|
|
{
|
|
mat = 1;
|
|
idx = inY * dimX1 + inX;
|
|
}
|
|
else if ( dimX1 <= inX && inX < (dimX1 + pretendDimX2) )
|
|
{
|
|
mat = 2;
|
|
idx = (inY + startY2) * dimX2 + (inX + startX2 - dimX1);
|
|
}
|
|
else if ( (dimX1 + pretendDimX2) <= inX && inX < (dimX1 + pretendDimX2 + dimX3) )
|
|
{
|
|
mat = 3;
|
|
idx = inY * dimX3 + (inX - dimX1 - pretendDimX2);
|
|
}
|
|
else
|
|
vtkm::cont::ErrorControlInternal("Invalid index!");
|
|
}
|
|
else // top-down mode
|
|
{
|
|
if ( 0 <= inY && inY < dimY1 )
|
|
{
|
|
mat = 1;
|
|
idx = inY * dimX1 + inX;
|
|
}
|
|
else if ( dimY1 <= inY && inY < (dimY1 + pretendDimY2) )
|
|
{
|
|
mat = 2;
|
|
idx = (inY + startY2 - dimY1) * dimX2 + inX + startX2;
|
|
}
|
|
else if ( (dimY1 + pretendDimY2) <= inY && inY < (dimY1 + pretendDimY2 + dimY3) )
|
|
{
|
|
mat = 3;
|
|
idx = (inY - dimY1 - pretendDimY2) * dimX3 + inX;
|
|
}
|
|
else
|
|
vtkm::cont::ErrorControlInternal("Invalid index!");
|
|
}
|
|
}
|
|
|
|
private:
|
|
const vtkm::Id dimX1, dimY1;
|
|
const vtkm::Id dimX2, dimY2, startX2, startY2, pretendDimX2, pretendDimY2;
|
|
const vtkm::Id dimX3, dimY3;
|
|
const bool mode_lr; // true: left right mode; false: top down mode.
|
|
};
|
|
|
|
|
|
|
|
// Worklet for 2D signal extension
|
|
// This implementation operates on a specified part of a big rectangle
|
|
class ExtensionWorklet2D : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
|
|
WholeArrayIn < ScalarAll > ); // signal part
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
ExtensionWorklet2D ( vtkm::Id extdimX, vtkm::Id extdimY,
|
|
vtkm::Id sigdimX, vtkm::Id sigdimY,
|
|
vtkm::Id sigstartX, vtkm::Id sigstartY,
|
|
vtkm::Id sigpretendX, vtkm::Id sigpretendY,
|
|
DWTMode m,
|
|
ExtensionDirection dir,
|
|
bool pad_zero)
|
|
:
|
|
extDimX( extdimX ), extDimY( extdimY ),
|
|
sigDimX( sigdimX ), sigDimY( sigdimY ),
|
|
sigStartX( sigstartX ), sigStartY( sigstartY ),
|
|
sigPretendDimX( sigpretendX ), sigPretendDimY( sigpretendY ),
|
|
mode( m ),
|
|
direction( dir ),
|
|
padZero( pad_zero )
|
|
{ (void)sigDimY; }
|
|
|
|
// Index translation helper
|
|
VTKM_EXEC_CONT
|
|
void Ext1Dto2D ( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y ) const
|
|
{
|
|
x = idx % extDimX;
|
|
y = idx / extDimX;
|
|
}
|
|
|
|
// Index translation helper
|
|
VTKM_EXEC_CONT
|
|
vtkm::Id Sig2Dto1D( vtkm::Id x, vtkm::Id y ) const
|
|
{
|
|
return y * sigDimX + x;
|
|
}
|
|
|
|
// Index translation helper
|
|
VTKM_EXEC_CONT
|
|
vtkm::Id SigPretend2Dto1D( vtkm::Id x, vtkm::Id y ) const
|
|
{
|
|
return (y + sigStartY) * sigDimX + x + sigStartX;
|
|
}
|
|
|
|
template< typename PortalOutType, typename PortalInType >
|
|
VTKM_EXEC
|
|
void operator()( PortalOutType &portalOut,
|
|
const PortalInType &portalIn,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
vtkm::Id extX, extY, sigPretendX, sigPretendY;
|
|
sigPretendX = sigPretendY = 0;
|
|
Ext1Dto2D( workIndex, extX, extY );
|
|
typename PortalOutType::ValueType sym = 1.0;
|
|
if( mode == ASYMH || mode == ASYMW )
|
|
sym = -1.0;
|
|
if( direction == LEFT )
|
|
{
|
|
sigPretendY = extY;
|
|
if( mode == SYMH || mode == ASYMH )
|
|
sigPretendX = extDimX - extX - 1;
|
|
else // mode == SYMW || mode == ASYMW
|
|
sigPretendX = extDimX - extX;
|
|
}
|
|
else if( direction == TOP )
|
|
{
|
|
sigPretendX = extX;
|
|
if( mode == SYMH || mode == ASYMH )
|
|
sigPretendY = extDimY - extY - 1;
|
|
else // mode == SYMW || mode == ASYMW
|
|
sigPretendY = extDimY - extY;
|
|
}
|
|
else if( direction == RIGHT )
|
|
{
|
|
sigPretendY = extY;
|
|
if( mode == SYMH || mode == ASYMH )
|
|
sigPretendX = sigPretendDimX - extX - 1;
|
|
else
|
|
sigPretendX = sigPretendDimX - extX - 2;
|
|
if( padZero )
|
|
sigPretendX++;
|
|
}
|
|
else if( direction == BOTTOM )
|
|
{
|
|
sigPretendX = extX;
|
|
if( mode == SYMH || mode == ASYMH )
|
|
sigPretendY = sigPretendDimY - extY - 1;
|
|
else
|
|
sigPretendY = sigPretendDimY - extY - 2;
|
|
if( padZero )
|
|
sigPretendY++;
|
|
}
|
|
else
|
|
vtkm::cont::ErrorControlInternal("Invalid extension mode for matrices!");
|
|
|
|
if( sigPretendX == sigPretendDimX || sigPretendY == sigPretendDimY )
|
|
portalOut.Set( workIndex, 0.0 );
|
|
else
|
|
portalOut.Set( workIndex, sym *
|
|
portalIn.Get( SigPretend2Dto1D(sigPretendX, sigPretendY) ));
|
|
}
|
|
|
|
private:
|
|
const vtkm::Id extDimX, extDimY, sigDimX, sigDimY;
|
|
const vtkm::Id sigStartX, sigStartY, sigPretendDimX, sigPretendDimY;
|
|
const DWTMode mode;
|
|
const ExtensionDirection direction;
|
|
const bool padZero; // treat sigIn as having a column/row zeros
|
|
};
|
|
|
|
|
|
|
|
// Worklet: perform a simple 2D forward transform
|
|
template< typename DeviceTag >
|
|
class ForwardTransform2D: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature(WholeArrayIn<ScalarAll>, // left/top extension
|
|
WholeArrayIn<ScalarAll>, // sigIn
|
|
WholeArrayIn<ScalarAll>, // right/bottom extension
|
|
WholeArrayOut<ScalarAll>); // cA followed by cD
|
|
typedef void ExecutionSignature(_1, _2, _3, _4, WorkIndex);
|
|
typedef _4 InputDomain;
|
|
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
ForwardTransform2D ( const vtkm::cont::ArrayHandle<vtkm::Float64> &loFilter,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &hiFilter,
|
|
vtkm::Id filter_len, vtkm::Id approx_len,
|
|
bool odd_low, bool mode_lr,
|
|
vtkm::Id x1, vtkm::Id y1, // dims of left/top extension
|
|
vtkm::Id x2, vtkm::Id y2, // dims of signal
|
|
vtkm::Id startx2, vtkm::Id starty2, // start idx of signal
|
|
vtkm::Id pretendx2, vtkm::Id pretendy2, // pretend dims of signal
|
|
vtkm::Id x3, vtkm::Id y3 ) // dims of right/bottom extension
|
|
:
|
|
lowFilter( loFilter.PrepareForInput( DeviceTag() ) ),
|
|
highFilter( hiFilter.PrepareForInput( DeviceTag() ) ),
|
|
filterLen( filter_len ), approxLen( approx_len ),
|
|
outDimX( pretendx2 ), outDimY( pretendy2 ),
|
|
oddlow( odd_low ), modeLR( mode_lr ),
|
|
translator( x1, y1,
|
|
x2, y2,
|
|
startx2, starty2,
|
|
pretendx2, pretendy2,
|
|
x3, y3,
|
|
mode_lr )
|
|
{ this->SetStartPosition(); }
|
|
|
|
VTKM_EXEC_CONT
|
|
void Output1Dto2D( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y ) const
|
|
{
|
|
x = idx % outDimX;
|
|
y = idx / outDimX;
|
|
}
|
|
VTKM_EXEC_CONT
|
|
vtkm::Id Output2Dto1D( vtkm::Id x, vtkm::Id y ) const
|
|
{
|
|
return y * outDimX + x;
|
|
}
|
|
|
|
// Use 64-bit float for convolution calculation
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
template <typename InPortalType1, typename InPortalType2, typename InPortalType3 >
|
|
VTKM_EXEC_CONT
|
|
VAL GetVal( const InPortalType1 &portal1, const InPortalType2 &portal2,
|
|
const InPortalType3 &portal3, vtkm::Id inMatrix, vtkm::Id inIdx ) const
|
|
{
|
|
if( inMatrix == 1 )
|
|
return MAKEVAL( portal1.Get(inIdx) );
|
|
else if( inMatrix == 2 )
|
|
return MAKEVAL( portal2.Get(inIdx) );
|
|
else if( inMatrix == 3 )
|
|
return MAKEVAL( portal3.Get(inIdx) );
|
|
else
|
|
{
|
|
vtkm::cont::ErrorControlInternal("Invalid matrix index!");
|
|
return -1;
|
|
}
|
|
}
|
|
|
|
template <typename InPortalType1, typename InPortalType2,
|
|
typename InPortalType3, typename OutputPortalType>
|
|
VTKM_EXEC_CONT
|
|
void operator()(const InPortalType1 &inPortal1, // left/top extension
|
|
const InPortalType2 &inPortal2, // signal
|
|
const InPortalType3 &inPortal3, // right/bottom extension
|
|
OutputPortalType &coeffOut,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
vtkm::Id workX, workY, output1D;
|
|
Output1Dto2D( workIndex, workX, workY );
|
|
vtkm::Id inputMatrix, inputIdx;
|
|
typedef typename OutputPortalType::ValueType OutputValueType;
|
|
|
|
if( modeLR )
|
|
{
|
|
if( workX % 2 == 0 ) // calculate cA
|
|
{
|
|
vtkm::Id xl = lstart + workX;
|
|
VAL sum = MAKEVAL(0.0);
|
|
for( vtkm::Id k = filterLen - 1; k > -1; k-- )
|
|
{
|
|
translator.Translate2Dto1D( xl, workY, inputMatrix, inputIdx );
|
|
sum += lowFilter.Get(k) *
|
|
GetVal( inPortal1, inPortal2, inPortal3, inputMatrix, inputIdx );
|
|
xl++;
|
|
}
|
|
output1D = Output2Dto1D( workX/2, workY );
|
|
coeffOut.Set( output1D, static_cast<OutputValueType>(sum) );
|
|
}
|
|
else // calculate cD
|
|
{
|
|
vtkm::Id xh = hstart + workX - 1;
|
|
VAL sum=MAKEVAL(0.0);
|
|
for( vtkm::Id k = filterLen - 1; k > -1; k-- )
|
|
{
|
|
translator.Translate2Dto1D( xh, workY, inputMatrix, inputIdx );
|
|
sum += highFilter.Get(k) *
|
|
GetVal( inPortal1, inPortal2, inPortal3, inputMatrix, inputIdx );
|
|
xh++;
|
|
}
|
|
output1D = Output2Dto1D( (workX-1)/2 + approxLen, workY );
|
|
coeffOut.Set( output1D, static_cast<OutputValueType>(sum) );
|
|
}
|
|
}
|
|
else // top-down order
|
|
{
|
|
if( workY % 2 == 0 ) // calculate cA
|
|
{
|
|
vtkm::Id yl = lstart + workY;
|
|
VAL sum = MAKEVAL(0.0);
|
|
for( vtkm::Id k = filterLen - 1; k > -1; k-- )
|
|
{
|
|
translator.Translate2Dto1D( workX, yl, inputMatrix, inputIdx );
|
|
sum += lowFilter.Get(k) *
|
|
GetVal( inPortal1, inPortal2, inPortal3, inputMatrix, inputIdx );
|
|
yl++;
|
|
}
|
|
output1D = Output2Dto1D( workX, workY/2 );
|
|
coeffOut.Set( output1D, static_cast<OutputValueType>(sum) );
|
|
}
|
|
else // calculate cD
|
|
{
|
|
vtkm::Id yh = hstart + workY - 1;
|
|
VAL sum=MAKEVAL(0.0);
|
|
for( vtkm::Id k = filterLen - 1; k > -1; k-- )
|
|
{
|
|
translator.Translate2Dto1D( workX, yh, inputMatrix, inputIdx );
|
|
sum += highFilter.Get(k) *
|
|
GetVal( inPortal1, inPortal2, inPortal3, inputMatrix, inputIdx );
|
|
yh++;
|
|
}
|
|
output1D = Output2Dto1D( workX, (workY-1)/2 + approxLen );
|
|
coeffOut.Set( output1D, static_cast<OutputValueType>(sum) );
|
|
}
|
|
}
|
|
}
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
private:
|
|
const typename vtkm::cont::ArrayHandle<vtkm::Float64>::ExecutionTypes<DeviceTag>::
|
|
PortalConst lowFilter, highFilter;
|
|
const vtkm::Id filterLen, approxLen;
|
|
const vtkm::Id outDimX, outDimY;
|
|
bool oddlow;
|
|
bool modeLR; // true = left right; false = top down.
|
|
const IndexTranslator3Matrices translator;
|
|
vtkm::Id lstart, hstart;
|
|
|
|
VTKM_EXEC_CONT
|
|
void SetStartPosition()
|
|
{
|
|
this->lstart = this->oddlow ? 1 : 0;
|
|
this->hstart = 1;
|
|
}
|
|
};
|
|
|
|
|
|
|
|
// ---------------------------------------------------
|
|
// | | | | | | |
|
|
// | | | | | | |
|
|
// | ext1 | cA | ext2 | ext3 | cD | ext4 |
|
|
// | (x1) | (xa) | (x2) | (x3) | (xd) | (x4) |
|
|
// | | | | | | |
|
|
// ----------------------------------------------------
|
|
// portal1: ext1
|
|
// portal2: ext2
|
|
// portal3: ext3
|
|
// portal4: ext4
|
|
// portal5: cA + cD
|
|
// Worklet: perform a simple 2D inverse transform
|
|
template< typename DeviceTag >
|
|
class InverseTransform2D: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayIn< ScalarAll >, // ext1
|
|
WholeArrayIn< ScalarAll >, // ext2
|
|
WholeArrayIn< ScalarAll >, // ext3
|
|
WholeArrayIn< ScalarAll >, // ext4
|
|
WholeArrayIn< ScalarAll >, // cA+cD (signal)
|
|
FieldOut< ScalarAll> ); // outptu coeffs
|
|
typedef void ExecutionSignature( _1, _2, _3, _4, _5, _6, WorkIndex );
|
|
typedef _6 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
InverseTransform2D( const vtkm::cont::ArrayHandle<vtkm::Float64> &lo_fil,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &hi_fil,
|
|
vtkm::Id fil_len,
|
|
vtkm::Id x_1, vtkm::Id y_1, // ext1
|
|
vtkm::Id x_a, vtkm::Id y_a, // cA
|
|
vtkm::Id x_2, vtkm::Id y_2, // ext2
|
|
vtkm::Id x_3, vtkm::Id y_3, // ext3
|
|
vtkm::Id x_d, vtkm::Id y_d, // cD
|
|
vtkm::Id x_4, vtkm::Id y_4, // ext4
|
|
vtkm::Id x_5, vtkm::Id y_5,
|
|
vtkm::Id startX5, vtkm::Id startY5,
|
|
bool mode_lr )
|
|
: lowFilter( lo_fil.PrepareForInput( DeviceTag() ) ),
|
|
highFilter( hi_fil.PrepareForInput( DeviceTag() ) ),
|
|
filterLen( fil_len ),
|
|
translator(x_1, y_1, x_a, y_a, x_2, y_2,
|
|
x_3, y_3, x_d, y_d, x_4, y_4,
|
|
x_5, y_5, startX5, startY5, mode_lr ),
|
|
modeLR( mode_lr )
|
|
{
|
|
if( modeLR )
|
|
{
|
|
outputDimX = x_a + x_d;
|
|
outputDimY = y_1;
|
|
cALenExtended = x_1 + x_a + x_2;
|
|
}
|
|
else
|
|
{
|
|
outputDimX = x_1;
|
|
outputDimY = y_a + y_d;
|
|
cALenExtended = y_1 + y_a + y_2;
|
|
}
|
|
}
|
|
|
|
VTKM_EXEC_CONT
|
|
void Output1Dto2D( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y ) const
|
|
{
|
|
x = idx % outputDimX;
|
|
y = idx / outputDimX;
|
|
}
|
|
|
|
// Use 64-bit float for convolution calculation
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
template <typename InPortalType1, typename InPortalType2, typename InPortalType3,
|
|
typename InPortalType4, typename InPortalTypecAcD >
|
|
VTKM_EXEC_CONT
|
|
VAL GetVal( const InPortalType1 &ext1,
|
|
const InPortalType2 &ext2,
|
|
const InPortalType3 &ext3,
|
|
const InPortalType4 &ext4,
|
|
const InPortalTypecAcD &cAcD,
|
|
vtkm::Id inMatrix, vtkm::Id inIdx ) const
|
|
{
|
|
if( inMatrix == 1 )
|
|
return MAKEVAL( ext1.Get(inIdx) );
|
|
else if( inMatrix == 2 )
|
|
return MAKEVAL( ext2.Get(inIdx) );
|
|
else if( inMatrix == 3 )
|
|
return MAKEVAL( ext3.Get(inIdx) );
|
|
else if( inMatrix == 4 )
|
|
return MAKEVAL( ext4.Get(inIdx) );
|
|
else if( inMatrix == 5 )
|
|
return MAKEVAL( cAcD.Get(inIdx) );
|
|
else
|
|
{
|
|
vtkm::cont::ErrorControlInternal("Invalid matrix index!");
|
|
return -1;
|
|
}
|
|
}
|
|
|
|
template< typename InPortalType1, typename InPortalType2, typename InPortalType3,
|
|
typename InPortalType4, typename InPortalTypecAcD,
|
|
typename OutputValueType >
|
|
VTKM_EXEC
|
|
void operator() (const InPortalType1 &portal1,
|
|
const InPortalType2 &portal2,
|
|
const InPortalType3 &portal3,
|
|
const InPortalType4 &portal4,
|
|
const InPortalTypecAcD &portalcAcD,
|
|
OutputValueType &coeffOut,
|
|
const vtkm::Id &workIdx ) const
|
|
{
|
|
vtkm::Id workX, workY;
|
|
vtkm::Id k1, k2, xi, yi, inputMatrix, inputIdx;
|
|
Output1Dto2D( workIdx, workX, workY );
|
|
|
|
// left-right, odd filter
|
|
if( modeLR && (filterLen % 2 != 0) )
|
|
{
|
|
if( workX % 2 != 0 )
|
|
{
|
|
k1 = filterLen - 2; k2 = filterLen - 1;
|
|
}
|
|
else
|
|
{
|
|
k1 = filterLen - 1; k2 = filterLen - 2;
|
|
}
|
|
|
|
VAL sum = 0.0;
|
|
xi = (workX + 1) / 2;
|
|
while( k1 > -1 )
|
|
{
|
|
translator.Translate2Dto1D( xi, workY, inputMatrix, inputIdx );
|
|
sum += lowFilter.Get(k1) * GetVal( portal1, portal2, portal3, portal4,
|
|
portalcAcD, inputMatrix, inputIdx );
|
|
xi++;
|
|
k1 -= 2;
|
|
}
|
|
xi = workX / 2;
|
|
while( k2 > -1 )
|
|
{
|
|
translator.Translate2Dto1D( xi + cALenExtended, workY, inputMatrix, inputIdx );
|
|
sum += highFilter.Get(k2) * GetVal( portal1, portal2, portal3, portal4,
|
|
portalcAcD, inputMatrix, inputIdx );
|
|
xi++;
|
|
k2 -= 2;
|
|
}
|
|
coeffOut = static_cast< OutputValueType> (sum);
|
|
}
|
|
|
|
// top-down, odd filter
|
|
else if ( !modeLR && (filterLen % 2 != 0) )
|
|
{
|
|
if( workY % 2 != 0 )
|
|
{
|
|
k1 = filterLen - 2; k2 = filterLen - 1;
|
|
}
|
|
else
|
|
{
|
|
k1 = filterLen - 1; k2 = filterLen - 2;
|
|
}
|
|
|
|
VAL sum = 0.0;
|
|
yi = (workY + 1) / 2;
|
|
while( k1 > -1 )
|
|
{
|
|
translator.Translate2Dto1D( workX, yi, inputMatrix, inputIdx );
|
|
VAL cA = GetVal( portal1, portal2, portal3, portal4, portalcAcD, inputMatrix, inputIdx );
|
|
sum += lowFilter.Get(k1) * cA;
|
|
yi++;
|
|
k1 -= 2;
|
|
}
|
|
yi = workY / 2;
|
|
while( k2 > -1 )
|
|
{
|
|
translator.Translate2Dto1D( workX, yi + cALenExtended, inputMatrix, inputIdx );
|
|
VAL cD = GetVal( portal1, portal2, portal3, portal4, portalcAcD, inputMatrix, inputIdx );
|
|
sum += highFilter.Get(k2) * cD;
|
|
yi++;
|
|
k2 -= 2;
|
|
}
|
|
coeffOut = static_cast< OutputValueType >(sum);
|
|
}
|
|
|
|
// left-right, even filter
|
|
else if( modeLR && (filterLen % 2 == 0) )
|
|
{
|
|
if( (filterLen/2) % 2 != 0 ) // odd length half filter
|
|
{
|
|
xi = workX / 2;
|
|
if( workX % 2 != 0 )
|
|
k1 = filterLen - 1;
|
|
else
|
|
k1 = filterLen - 2;
|
|
}
|
|
else // even length half filter
|
|
{
|
|
xi = (workX + 1) / 2;
|
|
if( workX % 2 != 0 )
|
|
k1 = filterLen - 2;
|
|
else
|
|
k1 = filterLen - 1;
|
|
}
|
|
VAL cA, cD;
|
|
VAL sum = 0.0;
|
|
while( k1 > -1 )
|
|
{
|
|
translator.Translate2Dto1D( xi, workY, inputMatrix, inputIdx );
|
|
cA = GetVal( portal1, portal2, portal3, portal4, portalcAcD,
|
|
inputMatrix, inputIdx );
|
|
translator.Translate2Dto1D( xi + cALenExtended, workY, inputMatrix, inputIdx );
|
|
cD = GetVal( portal1, portal2, portal3, portal4, portalcAcD,
|
|
inputMatrix, inputIdx );
|
|
sum += lowFilter.Get(k1) * cA + highFilter.Get(k1) * cD;
|
|
xi++;
|
|
k1 -= 2;
|
|
}
|
|
coeffOut = static_cast< OutputValueType >(sum);
|
|
}
|
|
|
|
// top-down, even filter
|
|
else
|
|
{
|
|
if( (filterLen/2) % 2 != 0 )
|
|
{
|
|
yi = workY / 2;
|
|
if( workY % 2 != 0 )
|
|
k1 = filterLen - 1;
|
|
else
|
|
k1 = filterLen - 2;
|
|
}
|
|
else
|
|
{
|
|
yi = (workY + 1) / 2;
|
|
if( workY % 2 != 0 )
|
|
k1 = filterLen - 2;
|
|
else
|
|
k1 = filterLen - 1;
|
|
}
|
|
VAL cA, cD;
|
|
VAL sum = 0.0;
|
|
while( k1 > -1 )
|
|
{
|
|
translator.Translate2Dto1D( workX, yi, inputMatrix, inputIdx );
|
|
cA = GetVal( portal1, portal2, portal3, portal4, portalcAcD,
|
|
inputMatrix, inputIdx );
|
|
translator.Translate2Dto1D( workX, yi + cALenExtended, inputMatrix, inputIdx );
|
|
cD = GetVal( portal1, portal2, portal3, portal4, portalcAcD,
|
|
inputMatrix, inputIdx );
|
|
sum += lowFilter.Get(k1) * cA + highFilter.Get(k1) * cD;
|
|
yi++;
|
|
k1 -= 2;
|
|
}
|
|
coeffOut = static_cast< OutputValueType >(sum);
|
|
}
|
|
}
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
private:
|
|
const typename vtkm::cont::ArrayHandle<vtkm::Float64>::
|
|
ExecutionTypes<DeviceTag>::PortalConst lowFilter, highFilter;
|
|
const vtkm::Id filterLen;
|
|
vtkm::Id outputDimX, outputDimY;
|
|
vtkm::Id cALenExtended; // Number of cA at the beginning of input, followed by cD
|
|
const IndexTranslator6Matrices translator;
|
|
const bool modeLR;
|
|
};
|
|
|
|
|
|
|
|
// Worklet: perform a simple 1D forward transform
|
|
template< typename DeviceTag >
|
|
class ForwardTransform: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature(WholeArrayIn<ScalarAll>, // sigIn
|
|
WholeArrayOut<ScalarAll>); // cA followed by cD
|
|
typedef void ExecutionSignature(_1, _2, WorkIndex);
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
ForwardTransform( const vtkm::cont::ArrayHandle<vtkm::Float64> &loFilter,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &hiFilter,
|
|
vtkm::Id filLen, vtkm::Id approx_len, vtkm::Id detail_len,
|
|
bool odd_low, bool odd_high ) :
|
|
lowFilter( loFilter.PrepareForInput(DeviceTag()) ),
|
|
highFilter( hiFilter.PrepareForInput(DeviceTag()) ),
|
|
filterLen( filLen ),
|
|
approxLen( approx_len ),
|
|
detailLen( detail_len ),
|
|
oddlow ( odd_low ),
|
|
oddhigh ( odd_high )
|
|
{ this->SetStartPosition(); }
|
|
|
|
// Use 64-bit float for convolution calculation
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
template <typename InputPortalType, typename OutputPortalType>
|
|
VTKM_EXEC
|
|
void operator()(const InputPortalType &signalIn,
|
|
OutputPortalType &coeffOut,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
typedef typename OutputPortalType::ValueType OutputValueType;
|
|
if( workIndex < approxLen + detailLen )
|
|
{
|
|
if( workIndex % 2 == 0 ) // calculate cA
|
|
{
|
|
vtkm::Id xl = xlstart + workIndex;
|
|
VAL sum=MAKEVAL(0.0);
|
|
for( vtkm::Id k = filterLen - 1; k >= 0; k-- )
|
|
sum += lowFilter.Get(k) * MAKEVAL( signalIn.Get(xl++) );
|
|
vtkm::Id outputIdx = workIndex / 2; // put cA at the beginning
|
|
coeffOut.Set( outputIdx, static_cast<OutputValueType>(sum) );
|
|
}
|
|
else // calculate cD
|
|
{
|
|
VAL sum=MAKEVAL(0.0);
|
|
vtkm::Id xh = xhstart + workIndex - 1;
|
|
for( vtkm::Id k = filterLen - 1; k >= 0; k-- )
|
|
sum += highFilter.Get(k) * MAKEVAL( signalIn.Get(xh++) );
|
|
vtkm::Id outputIdx = approxLen + (workIndex-1) / 2; // put cD after cA
|
|
coeffOut.Set( outputIdx, static_cast<OutputValueType>(sum) );
|
|
}
|
|
}
|
|
}
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
private:
|
|
const typename vtkm::cont::ArrayHandle<vtkm::Float64>::
|
|
ExecutionTypes<DeviceTag>::PortalConst lowFilter, highFilter;
|
|
const vtkm::Id filterLen, approxLen, detailLen; // filter and outcome coeff length.
|
|
bool oddlow, oddhigh;
|
|
vtkm::Id xlstart, xhstart;
|
|
|
|
VTKM_EXEC_CONT
|
|
void SetStartPosition()
|
|
{
|
|
this->xlstart = this->oddlow ? 1 : 0;
|
|
this->xhstart = this->oddhigh ? 1 : 0;
|
|
}
|
|
};
|
|
|
|
|
|
|
|
// Worklet: perform an 1D inverse transform for odd length, symmetric filters.
|
|
template< typename DeviceTag >
|
|
class InverseTransformOdd: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature(WholeArrayIn<ScalarAll>, // Input: coeffs,
|
|
// cA followed by cD
|
|
WholeArrayOut<ScalarAll>); // output
|
|
typedef void ExecutionSignature(_1, _2, WorkIndex);
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
InverseTransformOdd( const vtkm::cont::ArrayHandle<vtkm::Float64> &loFilter,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &hiFilter,
|
|
vtkm::Id filLen, vtkm::Id ca_len, vtkm::Id ext_len )
|
|
: lowFilter( loFilter.PrepareForInput(DeviceTag()) ),
|
|
highFilter( hiFilter.PrepareForInput(DeviceTag()) ),
|
|
filterLen( filLen ), cALen( ca_len ),
|
|
cALen2( ca_len * 2 ), cALenExtended( ext_len ) {}
|
|
|
|
// Use 64-bit float for convolution calculation
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
template <typename InputPortalType, typename OutputPortalType>
|
|
VTKM_EXEC
|
|
void operator()(const InputPortalType &coeffs,
|
|
OutputPortalType &sigOut,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
if( workIndex < cALen2 ) // valid calculation region
|
|
{
|
|
vtkm::Id xi; // coeff indices
|
|
vtkm::Id k1, k2; // indices for low and high filter
|
|
VAL sum = 0.0;
|
|
|
|
if( workIndex % 2 != 0 )
|
|
{
|
|
k1 = this->filterLen - 2;
|
|
k2 = this->filterLen - 1;
|
|
}
|
|
else
|
|
{
|
|
k1 = this->filterLen - 1;
|
|
k2 = this->filterLen - 2;
|
|
}
|
|
|
|
xi = (workIndex+1) / 2;
|
|
while( k1 > -1 ) // k1 >= 0
|
|
{
|
|
sum += lowFilter.Get(k1) * MAKEVAL( coeffs.Get(xi++) );
|
|
k1 -= 2;
|
|
}
|
|
|
|
xi = workIndex / 2;
|
|
while( k2 > -1 ) // k2 >= 0
|
|
{
|
|
sum += highFilter.Get(k2) * MAKEVAL( coeffs.Get( this->cALenExtended + xi++ ) );
|
|
k2 -= 2;
|
|
}
|
|
|
|
sigOut.Set(workIndex, static_cast<typename OutputPortalType::ValueType>( sum ) );
|
|
}
|
|
|
|
}
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
private:
|
|
const typename vtkm::cont::ArrayHandle<vtkm::Float64>::ExecutionTypes<DeviceTag>::
|
|
PortalConst lowFilter, highFilter;
|
|
const vtkm::Id filterLen; // filter length.
|
|
const vtkm::Id cALen; // Number of actual cAs
|
|
const vtkm::Id cALen2; // = cALen * 2
|
|
const vtkm::Id cALenExtended; // Number of cA at the beginning of input, followed by cD
|
|
};
|
|
|
|
|
|
|
|
// Worklet: perform an 1D inverse transform for even length, symmetric filters.
|
|
template< typename DeviceTag >
|
|
class InverseTransformEven: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature(WholeArrayIn<ScalarAll>, // Input: coeffs,
|
|
// cA followed by cD
|
|
WholeArrayOut<ScalarAll>); // output
|
|
typedef void ExecutionSignature(_1, _2, WorkIndex);
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
InverseTransformEven( const vtkm::cont::ArrayHandle<vtkm::Float64> &loFilter,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64> &hiFilter,
|
|
vtkm::Id filtL, vtkm::Id cAL, vtkm::Id cALExt, bool m )
|
|
: lowFilter( loFilter.PrepareForInput(DeviceTag()) ),
|
|
highFilter( hiFilter.PrepareForInput(DeviceTag()) ),
|
|
filterLen( filtL ), cALen( cAL ), cALen2( cAL * 2 ),
|
|
cALenExtended( cALExt ), matlab( m ) {}
|
|
|
|
// Use 64-bit float for convolution calculation
|
|
#define VAL vtkm::Float64
|
|
#define MAKEVAL(a) (static_cast<VAL>(a))
|
|
template <typename InputPortalType, typename OutputPortalType>
|
|
VTKM_EXEC
|
|
void operator()(const InputPortalType &coeffs,
|
|
OutputPortalType &sigOut,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
if( workIndex < cALen2 ) // valid calculation region
|
|
{
|
|
vtkm::Id xi; // coeff indices
|
|
vtkm::Id k; // indices for low and high filter
|
|
VAL sum = 0.0;
|
|
|
|
if( matlab || (filterLen/2) % 2 != 0 ) // odd length half filter
|
|
{
|
|
xi = workIndex / 2;
|
|
if( workIndex % 2 != 0 )
|
|
k = filterLen - 1;
|
|
else
|
|
k = filterLen - 2;
|
|
}
|
|
else
|
|
{
|
|
xi = (workIndex + 1) / 2;
|
|
if( workIndex % 2 != 0 )
|
|
k = filterLen - 2;
|
|
else
|
|
k = filterLen - 1;
|
|
}
|
|
|
|
while( k > -1 ) // k >= 0
|
|
{
|
|
sum += lowFilter.Get(k) * MAKEVAL( coeffs.Get( xi ) ) + // cA
|
|
highFilter.Get(k) * MAKEVAL( coeffs.Get( xi + cALenExtended) ); // cD
|
|
xi++;
|
|
k -= 2;
|
|
}
|
|
|
|
sigOut.Set(workIndex, static_cast<typename OutputPortalType::ValueType>( sum ) );
|
|
}
|
|
}
|
|
#undef MAKEVAL
|
|
#undef VAL
|
|
|
|
private:
|
|
const typename vtkm::cont::ArrayHandle<vtkm::Float64>::ExecutionTypes<DeviceTag>::
|
|
PortalConst lowFilter, highFilter;
|
|
const vtkm::Id filterLen; // filter length.
|
|
const vtkm::Id cALen; // Number of actual cAs
|
|
const vtkm::Id cALen2; // = cALen * 2
|
|
const vtkm::Id cALenExtended; // Number of cA at the beginning of input, followed by cD
|
|
bool matlab; // followed the naming convention from VAPOR
|
|
// It's always false for the 1st 4 filters.
|
|
};
|
|
|
|
|
|
|
|
class ThresholdWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( FieldInOut <ScalarAll> ); // Thresholding in-place
|
|
typedef void ExecutionSignature( _1 );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
ThresholdWorklet( vtkm::Float64 t ) : threshold( t ), // must pass in a positive val
|
|
neg_threshold( t*-1.0 ) {}
|
|
|
|
template <typename ValueType >
|
|
VTKM_EXEC
|
|
void operator()( ValueType &coeffVal ) const
|
|
{
|
|
if( neg_threshold < coeffVal && coeffVal < threshold )
|
|
coeffVal = 0.0;
|
|
}
|
|
|
|
private:
|
|
vtkm::Float64 threshold; // positive
|
|
vtkm::Float64 neg_threshold; // negative
|
|
};
|
|
|
|
|
|
|
|
class SquaredDeviation: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature(FieldIn<ScalarAll>,
|
|
FieldOut<ScalarAll>);
|
|
typedef _2 ExecutionSignature( _1 );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
template <typename ValueType>
|
|
VTKM_EXEC_CONT
|
|
SquaredDeviation( ValueType t )
|
|
{
|
|
this->mean = static_cast<vtkm::Float64>(t);
|
|
}
|
|
|
|
template <typename ValueType>
|
|
VTKM_EXEC
|
|
ValueType operator()( const ValueType &num ) const
|
|
{
|
|
vtkm::Float64 num64 = static_cast<vtkm::Float64>( num );
|
|
vtkm::Float64 diff = this->mean - num64;
|
|
return static_cast<ValueType>( diff * diff );
|
|
}
|
|
|
|
private:
|
|
vtkm::Float64 mean;
|
|
};
|
|
|
|
|
|
|
|
class Differencer: public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature(FieldIn<ScalarAll>,
|
|
FieldIn<ScalarAll>,
|
|
FieldOut<ScalarAll>);
|
|
typedef _3 ExecutionSignature( _1, _2 );
|
|
typedef _1 InputDomain;
|
|
|
|
template <typename ValueType1, typename ValueType2 >
|
|
VTKM_EXEC
|
|
ValueType1 operator()( const ValueType1 &v1, const ValueType2 &v2 ) const
|
|
{
|
|
return v1 - static_cast<ValueType1>(v2);
|
|
}
|
|
};
|
|
|
|
|
|
|
|
class SquareWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature(FieldIn< ScalarAll>,
|
|
FieldOut<ScalarAll>);
|
|
typedef _2 ExecutionSignature( _1 );
|
|
typedef _1 InputDomain;
|
|
|
|
template <typename ValueType>
|
|
VTKM_EXEC
|
|
ValueType operator()( const ValueType &v ) const
|
|
{
|
|
return (v * v);
|
|
}
|
|
};
|
|
|
|
|
|
|
|
class CopyWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayIn< ScalarAll >,
|
|
WholeArrayOut< ScalarAll > );
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
CopyWorklet( vtkm::Id idx )
|
|
{
|
|
this->startIdx = idx;
|
|
}
|
|
|
|
template< typename PortalInType, typename PortalOutType >
|
|
VTKM_EXEC
|
|
void operator()( const PortalInType &portalIn,
|
|
PortalOutType &portalOut,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
portalOut.Set( (startIdx + workIndex), portalIn.Get(workIndex) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id startIdx;
|
|
};
|
|
|
|
|
|
|
|
// Worklet for 1D signal extension no. 1
|
|
class LeftSYMHExtentionWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
|
|
WholeArrayIn < ScalarAll > ); // signal part
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
LeftSYMHExtentionWorklet( vtkm::Id len ) : addLen( len ) {}
|
|
|
|
template< typename PortalOutType, typename PortalInType >
|
|
VTKM_EXEC
|
|
void operator()( PortalOutType &portalOut,
|
|
const PortalInType &portalIn,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
portalOut.Set( workIndex, portalIn.Get(this->addLen - workIndex - 1) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id addLen;
|
|
};
|
|
|
|
|
|
|
|
// Worklet for 1D signal extension no. 2
|
|
class LeftSYMWExtentionWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
|
|
WholeArrayIn < ScalarAll > ); // signal part
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
LeftSYMWExtentionWorklet( vtkm::Id len ) : addLen( len ) {}
|
|
|
|
template< typename PortalOutType, typename PortalInType >
|
|
VTKM_EXEC
|
|
void operator()( PortalOutType &portalOut,
|
|
const PortalInType &portalIn,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
portalOut.Set( workIndex, portalIn.Get(this->addLen - workIndex) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id addLen;
|
|
};
|
|
|
|
|
|
|
|
// Worklet for 1D signal extension no. 3
|
|
class LeftASYMHExtentionWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
|
|
WholeArrayIn < ScalarAll > ); // signal part
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
LeftASYMHExtentionWorklet( vtkm::Id len ) : addLen (len) {}
|
|
|
|
template< typename PortalOutType, typename PortalInType >
|
|
VTKM_EXEC_CONT
|
|
void operator()( PortalOutType &portalOut,
|
|
const PortalInType &portalIn,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
portalOut.Set( workIndex, portalIn.Get( addLen - workIndex - 1) * (-1.0) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id addLen;
|
|
};
|
|
|
|
|
|
|
|
// Worklet for 1D signal extension no. 4
|
|
class LeftASYMWExtentionWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
|
|
WholeArrayIn < ScalarAll > ); // signal part
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
LeftASYMWExtentionWorklet( vtkm::Id len ) : addLen (len) {}
|
|
|
|
template< typename PortalOutType, typename PortalInType >
|
|
VTKM_EXEC_CONT
|
|
void operator()( PortalOutType &portalOut,
|
|
const PortalInType &portalIn,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
portalOut.Set( workIndex, portalIn.Get( addLen - workIndex ) * (-1.0) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id addLen;
|
|
};
|
|
|
|
|
|
|
|
// Worklet for 1D signal extension no. 5
|
|
class RightSYMHExtentionWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
|
|
WholeArrayIn < ScalarAll > ); // signal part
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
RightSYMHExtentionWorklet ( vtkm::Id sigInl ) : sigInLen( sigInl ) {}
|
|
|
|
template< typename PortalOutType, typename PortalInType >
|
|
VTKM_EXEC
|
|
void operator()( PortalOutType &portalOut,
|
|
const PortalInType &portalIn,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
portalOut.Set( workIndex, portalIn.Get(this->sigInLen - workIndex - 1) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id sigInLen;
|
|
};
|
|
|
|
|
|
|
|
// Worklet for 1D signal extension no. 6
|
|
class RightSYMWExtentionWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
|
|
WholeArrayIn < ScalarAll > ); // signal part
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
RightSYMWExtentionWorklet ( vtkm::Id sigInl ) : sigInLen( sigInl ) {}
|
|
|
|
template< typename PortalOutType, typename PortalInType >
|
|
VTKM_EXEC
|
|
void operator()( PortalOutType &portalOut,
|
|
const PortalInType &portalIn,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
portalOut.Set( workIndex, portalIn.Get(this->sigInLen - workIndex - 2) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id sigInLen;
|
|
};
|
|
|
|
|
|
|
|
// Worklet for 1D signal extension no. 7
|
|
class RightASYMHExtentionWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
|
|
WholeArrayIn < ScalarAll > ); // signal part
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
RightASYMHExtentionWorklet ( vtkm::Id sigInl ) : sigInLen( sigInl ) {}
|
|
|
|
template< typename PortalOutType, typename PortalInType >
|
|
VTKM_EXEC_CONT
|
|
void operator()( PortalOutType &portalOut,
|
|
const PortalInType &portalIn,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
portalOut.Set( workIndex, portalIn.Get( sigInLen - workIndex - 1) * (-1.0) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id sigInLen;
|
|
};
|
|
|
|
|
|
|
|
// Worklet for 1D signal extension no. 8
|
|
class RightASYMWExtentionWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayOut < ScalarAll >, // extension part
|
|
WholeArrayIn < ScalarAll > ); // signal part
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
typedef _1 InputDomain;
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
RightASYMWExtentionWorklet ( vtkm::Id sigInl ) : sigInLen( sigInl ) {}
|
|
|
|
template< typename PortalOutType, typename PortalInType >
|
|
VTKM_EXEC_CONT
|
|
void operator()( PortalOutType &portalOut,
|
|
const PortalInType &portalIn,
|
|
const vtkm::Id &workIndex) const
|
|
{
|
|
portalOut.Set( workIndex, portalIn.Get( sigInLen - workIndex - 2) * (-1.0) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id sigInLen;
|
|
};
|
|
|
|
|
|
|
|
// Assign zero to a single index
|
|
class AssignZeroWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayInOut< ScalarAll > );
|
|
typedef void ExecutionSignature( _1, WorkIndex );
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
AssignZeroWorklet( vtkm::Id idx ) : zeroIdx( idx ) { }
|
|
|
|
template< typename PortalType >
|
|
VTKM_EXEC
|
|
void operator()( PortalType &array,
|
|
const vtkm::Id &workIdx ) const
|
|
{
|
|
if( workIdx == this->zeroIdx )
|
|
array.Set( workIdx, static_cast<typename PortalType::ValueType>(0.0) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id zeroIdx;
|
|
};
|
|
|
|
|
|
|
|
// Assign zero to a row or a column in a 2D array.
|
|
// Change row or column is controlled by negative indices.
|
|
class AssignZero2DWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayInOut< ScalarAll > );
|
|
typedef void ExecutionSignature( _1, WorkIndex );
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
AssignZero2DWorklet( vtkm::Id x, vtkm::Id y, vtkm::Id zero_x, vtkm::Id zero_y )
|
|
: dimX( x ), dimY( y ), zeroX( zero_x ), zeroY( zero_y )
|
|
{ (void)dimY; }
|
|
|
|
// Index translation helper
|
|
VTKM_EXEC_CONT
|
|
void GetLogicalDim( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y ) const
|
|
{
|
|
x = idx % dimX;
|
|
y = idx / dimX;
|
|
}
|
|
|
|
template< typename PortalType >
|
|
VTKM_EXEC
|
|
void operator()( PortalType &array,
|
|
const vtkm::Id &workIdx ) const
|
|
{
|
|
vtkm::Id x, y;
|
|
GetLogicalDim( workIdx, x, y );
|
|
if( zeroY < 0 && x == zeroX ) // assign zero to a column
|
|
array.Set( workIdx, static_cast<typename PortalType::ValueType>(0.0) );
|
|
else if( zeroX < 0 && y == zeroY ) // assign zero to a row
|
|
array.Set( workIdx, static_cast<typename PortalType::ValueType>(0.0) );
|
|
// else if( x == zeroX && y == zeroY ) // assign zero to an element
|
|
// array.Set( workIdx, static_cast<typename PortalType::ValueType>(0.0) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id dimX, dimY;
|
|
vtkm::Id zeroX, zeroY; // element at (zeroX, zeroY) will be assigned zero.
|
|
// each becomes a wild card if negative
|
|
};
|
|
|
|
|
|
|
|
// Assign zero to a plane (2D) in a 3D cube.
|
|
// Which plane to assign zero is controlled by negative indices.
|
|
class AssignZero3DWorklet : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( WholeArrayInOut< ScalarAll > );
|
|
typedef void ExecutionSignature( _1, WorkIndex );
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
AssignZero3DWorklet( vtkm::Id x, vtkm::Id y, vtkm::Id z,
|
|
vtkm::Id zero_x, vtkm::Id zero_y, vtkm::Id zero_z )
|
|
: dimX( x ), dimY( y ), dimZ( z ),
|
|
zeroX( zero_x ), zeroY( zero_y ), zeroZ( zero_z )
|
|
{ }
|
|
|
|
// Index translation helper
|
|
VTKM_EXEC_CONT
|
|
void GetLogicalDim( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y, vtkm::Id &z ) const
|
|
{
|
|
z = idx / (dimX * dimY);
|
|
y = (idx - z * dimX * dimY) / dimX;
|
|
x = idx % dimX;
|
|
}
|
|
|
|
template< typename PortalType >
|
|
VTKM_EXEC
|
|
void operator()( PortalType &array,
|
|
const vtkm::Id &workIdx ) const
|
|
{
|
|
vtkm::Id x, y, z;
|
|
GetLogicalDim( workIdx, x, y, z );
|
|
if( zeroZ < 0 && zeroY < 0 && x == zeroX ) // plane perpendicular to X axis
|
|
array.Set( workIdx, static_cast<typename PortalType::ValueType>(0.0) );
|
|
else if( zeroZ < 0 && zeroX < 0 && y == zeroY ) // plane perpendicular to Y axis
|
|
array.Set( workIdx, static_cast<typename PortalType::ValueType>(0.0) );
|
|
else if( zeroY < 0 && zeroX < 0 && z == zeroZ ) // plane perpendicular to Z axis
|
|
array.Set( workIdx, static_cast<typename PortalType::ValueType>(0.0) );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id dimX, dimY, dimZ;
|
|
vtkm::Id zeroX, zeroY, zeroZ; // element at (zeroX, zeroY, zeroZ) will be assigned zero.
|
|
// each becomes a wild card if negative
|
|
};
|
|
|
|
|
|
|
|
// Worklet: Copys a small rectangle to become a part of a big rectangle
|
|
class RectangleCopyTo : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( FieldIn< ScalarAll >, // Input, small rectangle
|
|
WholeArrayOut< ScalarAll > ); // Output, big rectangle
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
RectangleCopyTo( vtkm::Id inx, vtkm::Id iny,
|
|
vtkm::Id outx, vtkm::Id outy,
|
|
vtkm::Id xStart, vtkm::Id yStart )
|
|
:
|
|
inXLen (inx), inYLen ( iny ),
|
|
outXLen (outx), outYLen ( outy ),
|
|
outXStart(xStart), outYStart( yStart )
|
|
{ }
|
|
|
|
VTKM_EXEC_CONT
|
|
void GetLogicalDimOfInputRect( const vtkm::Id &idx,
|
|
vtkm::Id &x,
|
|
vtkm::Id &y ) const
|
|
{
|
|
x = idx % inXLen;
|
|
y = idx / inXLen;
|
|
}
|
|
|
|
VTKM_EXEC_CONT
|
|
vtkm::Id Get1DIdxOfOutputRect( vtkm::Id x,
|
|
vtkm::Id y ) const
|
|
{
|
|
return y * outXLen + x;
|
|
}
|
|
|
|
template< typename ValueInType, typename PortalOutType >
|
|
VTKM_EXEC
|
|
void operator()( const ValueInType &valueIn,
|
|
PortalOutType &arrayOut,
|
|
const vtkm::Id &workIdx ) const
|
|
{
|
|
vtkm::Id xOfIn, yOfIn;
|
|
GetLogicalDimOfInputRect( workIdx, xOfIn, yOfIn );
|
|
vtkm::Id outputIdx = Get1DIdxOfOutputRect( xOfIn+outXStart, yOfIn+outYStart );
|
|
arrayOut.Set( outputIdx, valueIn );
|
|
}
|
|
|
|
private:
|
|
vtkm::Id inXLen, inYLen;
|
|
vtkm::Id outXLen, outYLen;
|
|
vtkm::Id outXStart, outYStart;
|
|
};
|
|
|
|
|
|
|
|
// Worklet: Copys a small cube to become a part of a big cube
|
|
class CubeCopyTo : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature( FieldIn< ScalarAll >, // Input, small cube
|
|
WholeArrayOut< ScalarAll > ); // Output, big cube
|
|
typedef void ExecutionSignature( _1, _2, WorkIndex );
|
|
|
|
// Constructor
|
|
VTKM_EXEC_CONT
|
|
CubeCopyTo( vtkm::Id inx, vtkm::Id iny, vtkm::Id inz,
|
|
vtkm::Id outx, vtkm::Id outy, vtkm::Id outz,
|
|
vtkm::Id xStart, vtkm::Id yStart, vtkm::Id zStart )
|
|
:
|
|
inDimX( inx ), inDimY( iny ), inDimZ( inz ),
|
|
outDimX( outx ), outDimY( outy ), outDimZ( outz ),
|
|
outStartX( xStart), outStartY( yStart ), outStartZ( zStart )
|
|
{ }
|
|
|
|
VTKM_EXEC_CONT
|
|
void GetLogicalDimOfInputCube( vtkm::Id idx,
|
|
vtkm::Id &x, vtkm::Id &y, vtkm::Id &z ) const
|
|
{
|
|
z = idx / (inDimX * inDimY);
|
|
y = (idx - z * inDimX * inDimY) / inDimX;
|
|
x = idx % inDimX;
|
|
}
|
|
|
|
VTKM_EXEC_CONT
|
|
vtkm::Id Get1DIdxOfOutputCube( vtkm::Id x, vtkm::Id y, vtkm::Id z ) const
|
|
{
|
|
return z * outDimX * outDimZ + y * outDimX + x;
|
|
}
|
|
|
|
template< typename ValueInType, typename PortalOutType >
|
|
VTKM_EXEC
|
|
void operator()( const ValueInType &valueIn,
|
|
PortalOutType &arrayOut,
|
|
const vtkm::Id &workIdx ) const
|
|
{
|
|
vtkm::Id inX, inY, inZ;
|
|
GetLogicalDimOfInputCube( workIdx, inX, inY, inZ );
|
|
vtkm::Id outputIdx = Get1DIdxOfOutputCube( inX+outStartX, inY+outStartY, inZ+outStartZ );
|
|
arrayOut.Set( outputIdx, valueIn );
|
|
}
|
|
|
|
private:
|
|
const vtkm::Id inDimX, inDimY, inDimZ; // input small cube
|
|
const vtkm::Id outDimX, outDimY, outDimZ; // output big cube
|
|
const vtkm::Id outStartX, outStartY, outStartZ; // where to put
|
|
};
|
|
|
|
|
|
} // namespace wavelets
|
|
} // namespace worlet
|
|
} // namespace vtkm
|
|
|
|
#endif // vtk_m_worklet_Wavelets_h
|