Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

sync to laptop

Browse Source
This commit is contained in:
Samuel Li 2016-09-21 20:25:56 -07:00
parent 9a3c2ea2a0
commit d19bbfac81
4 changed files with 261 additions and 147 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
vtkm/worklet/WaveletCompressor.h
View File

@ -190,14 +190,22 @@ public:
vtkm::Id currentLenX = inX;
vtkm::Id currentLenY = inY;
std::vector<vtkm::Id> L2d(10, 0);
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
vtkm::Float64 computationTime = 0.0;
typedef typename OutArrayType::ValueType OutValueType;
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
vtkm::Float64 computationTime = 0.0;
for( vtkm::Id i = nLevels; i > 0; i-- )
//vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( sigIn, coeffOut );
// First level transform operates on the input array
computationTime +=
WaveletDWT::DWT2Dv3( sigIn, currentLenX, currentLenY, coeffOut, L2d, DeviceTag());
VTKM_ASSERT( coeffOut.GetNumberOfValues() == currentLenX * currentLenY );
currentLenX = WaveletBase::GetApproxLength( currentLenX );
currentLenY = WaveletBase::GetApproxLength( currentLenY );
// Successor transforms operate on a temporary array
for( vtkm::Id i = nLevels-1; i > 0; i-- )
{
// make temporary input array
OutBasicArray tempInput;
@ -207,7 +215,7 @@ public:
OutBasicArray tempOutput;
computationTime +=
WaveletDWT::DWT2Dv2( tempInput, currentLenX, currentLenY, tempOutput, L2d, DeviceTag());
WaveletDWT::DWT2Dv3( tempInput, currentLenX, currentLenY, tempOutput, L2d, DeviceTag());
// copy results to coeffOut
WaveletBase::DeviceRectangleCopyTo( tempOutput, currentLenX, currentLenY,
@ -240,17 +248,21 @@ public:
{
throw vtkm::cont::ErrorControlBadValue("Number of levels of transform is not supported! ");
}
// fill the output array
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, arrOut );
if( nLevels == 0 ) // 0 levels means no transform
{
return 0;
}
VTKM_ASSERT( vtkm::Id(L.size()) == 6 * nLevels + 4 );
typedef typename OutArrayType::ValueType OutValueType;
typedef vtkm::cont::ArrayHandle<OutValueType> OutBasicArray;
OutBasicArray outBuffer;
vtkm::Float64 computationTime = 0.0;
if( nLevels == 0 ) // 0 levels means no transform
{
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, arrOut );
return 0;
}
else
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Copy( arrIn, outBuffer );
VTKM_ASSERT( vtkm::Id(L.size()) == 6 * nLevels + 4 );
std::vector<vtkm::Id> L2d(10, 0);
L2d[0] = L[0];
L2d[1] = L[1];
@ -261,8 +273,8 @@ public:
L2d[6] = L[6];
L2d[7] = L[7];
vtkm::Float64 computationTime = 0.0;
for( size_t i = 1; i <= static_cast<size_t>(nLevels); i++ )
// All transforms but the last operate on temporary arrays
for( size_t i = 1; i < static_cast<size_t>(nLevels); i++ )
{
L2d[8] = L2d[0] + L2d[4]; // This is always true for Biorthogonal wavelets
L2d[9] = L2d[1] + L2d[3]; // (same above)
@ -270,15 +282,15 @@ public:
// make input, output array
OutBasicArray tempInput, tempOutput;
WaveletBase::DeviceRectangleCopyFrom( tempInput, L2d[8], L2d[9],
arrOut, inX, inY, 0, 0, DeviceTag() );
outBuffer, inX, inY, 0, 0, DeviceTag() );
// IDWT
computationTime +=
WaveletDWT::IDWT2Dv2( tempInput, L2d, tempOutput, DeviceTag() );
WaveletDWT::IDWT2Dv3( tempInput, L2d, tempOutput, DeviceTag() );
// copy back reconstructed block
WaveletBase::DeviceRectangleCopyTo( tempOutput, L2d[8], L2d[9],
arrOut, inX, inY, 0, 0, DeviceTag() );
outBuffer, inX, inY, 0, 0, DeviceTag() );
// update L2d array
L2d[0] = L2d[8];
@ -291,6 +303,11 @@ public:
L2d[7] = L[6*i+7];
}
// The last transform takes place in place
L2d[8] = L2d[0] + L2d[4];
L2d[9] = L2d[1] + L2d[3];
computationTime += WaveletDWT::IDWT2Dv3( outBuffer, L2d, arrOut, DeviceTag() );
return computationTime;
}

145
vtkm/worklet/testing/UnitTestWaveletCompressor.cxx
View File

@ -33,6 +33,7 @@ namespace worklet
{
namespace wavelets
{
class SineWorklet : public vtkm::worklet::WorkletMapField
{
public:
@ -46,6 +47,56 @@ public:
x = vtkm::Sin(vtkm::Float64(workIdx) / 100.0) * 100.0;
}
};
class GaussianWorklet2D : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldInOut<>);
typedef void ExecutionSignature(_1, WorkIndex);
VTKM_EXEC_EXPORT
GaussianWorklet2D( vtkm::Id dx, vtkm::Id dy, vtkm::Float64 a,
vtkm::Float64 x, vtkm::Float64 y,
vtkm::Float64 sx, vtkm::Float64 xy )
: dimX( dx ), dimY( dy ), amp (a),
x0( x ), y0( y ),
sigmaX( sx ), sigmaY( xy )
{
sigmaX2 = 2 * sigmaX * sigmaX;
sigmaY2 = 2 * sigmaY * sigmaY;
}
VTKM_EXEC_EXPORT
void Sig1Dto2D( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y ) const
{
x = idx % dimX;
y = idx / dimX;
}
VTKM_EXEC_EXPORT
vtkm::Float64 GetGaussian( vtkm::Float64 x, vtkm::Float64 y ) const
{
vtkm::Float64 power = (x-x0) * (x-x0) / sigmaX2 + (y-y0) * (y-y0) / sigmaY2;
return vtkm::Exp( power * -1.0 ) * amp;
}
template<typename T>
VTKM_EXEC_EXPORT
void operator()(T& val, const vtkm::Id& workIdx) const
{
vtkm::Id x, y;
Sig1Dto2D( workIdx, x, y );
val = GetGaussian( static_cast<vtkm::Float64>(x), static_cast<vtkm::Float64>(y) );
}
private: // see wikipedia page
const vtkm::Id dimX, dimY; // 2D extent
const vtkm::Float64 amp; // amplitude
const vtkm::Float64 x0, y0; // center
const vtkm::Float64 sigmaX, sigmaY; // spread
vtkm::Float64 sigmaX2, sigmaY2; // 2 * sigma * sigma
};
}
}
}
@ -59,48 +110,24 @@ void FillArray( ArrayType& array )
dispatcher.Invoke( array );
}
void DebugExtend2D()
template< typename ArrayType >
void FillArray2D( ArrayType& array, vtkm::Id dimX, vtkm::Id dimY )
{
vtkm::Id NX = 10;
vtkm::Id NY = 10;
vtkm::Id addLen = 4;
typedef vtkm::cont::ArrayHandle< vtkm::Float64 > ArrayType;
ArrayType left1, left2, center, right1, right2;
ArrayType centerExtended;
center.PrepareForOutput( NX * NY, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
for( vtkm::Id i = 0; i < NX*NY; i++ )
center.GetPortalControl().Set(i, i);
typedef vtkm::worklet::wavelets::ExtensionWorklet2D ExtWorklet;
typedef vtkm::worklet::wavelets::LeftExtensionWorklet2D LeftExtWorklet;
typedef vtkm::worklet::wavelets::RightExtensionWorklet2D RightExtWorklet;
vtkm::worklet::wavelets::ExtensionDirection2D
extdirLeft = vtkm::worklet::wavelets::ExtensionDirection2D::LEFT;
vtkm::worklet::wavelets::ExtensionDirection2D
extdirRight = vtkm::worklet::wavelets::ExtensionDirection2D::RIGHT;
vtkm::worklet::wavelets::DWTMode mode = vtkm::worklet::wavelets::SYMW;
vtkm::worklet::wavelets::WaveletDWT dwt( vtkm::worklet::wavelets::CDF9_7 );
// compute real values
{
dwt.Extend2D( center, NX, NY, centerExtended, addLen, mode, mode, false, true,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
}
// compute test implementation
{
dwt.Extend2Dv3( center, NX, NY, left1, right1, addLen, mode, mode, false, true, true,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
}
typedef vtkm::worklet::wavelets::GaussianWorklet2D WorkletType;
WorkletType worklet( dimX, dimY, 3200.0,
static_cast<vtkm::Float64>(dimX)/2.0, // center
static_cast<vtkm::Float64>(dimY)/2.0, // center
static_cast<vtkm::Float64>(dimX)/4.0, // spread
static_cast<vtkm::Float64>(dimY)/4.0);// spread
vtkm::worklet::DispatcherMapField< WorkletType > dispatcher( worklet );
dispatcher.Invoke( array );
}
void DebugDWT2D()
{
vtkm::Id NX = 6;
vtkm::Id NY = 5;
vtkm::Id NX = 4;
vtkm::Id NY = 4;
typedef vtkm::cont::ArrayHandle< vtkm::Float64 > ArrayType;
ArrayType left, center, right;
@ -136,7 +163,6 @@ ArrayType idwt_out1, idwt_out2;
// test results go through IDWT
dwt.IDWT2Dv3( output3, L, idwt_out2, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
std::cerr << "finish IDWT2Dv3" << std::endl;
std::cout << "\ntrue results after IDWT:" << std::endl;
for( vtkm::Id i = 0; i < idwt_out1.GetNumberOfValues(); i++ )
@ -238,23 +264,22 @@ void DebugRectangleCopy()
void TestDecomposeReconstruct2D()
{
vtkm::Id sigX = 1000;
vtkm::Id sigY = 2000;
std::cout << "Please input X to test a X^2 square: " << std::endl;
std::cin >> sigX;
sigY = sigX;
//std::cout << "Testing wavelet compressor on 1000x2000 rectangle" << std::endl;
vtkm::Id sigX = 32768;
vtkm::Id sigY = 32768;
//std::cout << "Please input X to test a X^2 square: " << std::endl;
//std::cin >> sigX;
//sigY = sigX;
vtkm::Id sigLen = sigX * sigY;
// make input data array handle
vtkm::cont::ArrayHandle<vtkm::Float64> inputArray;
inputArray.PrepareForOutput( sigLen, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
FillArray( inputArray );
FillArray2D( inputArray, sigX, sigY );
vtkm::cont::ArrayHandle<vtkm::Float64> outputArray;
// Use a WaveletCompressor
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::HAAR;
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF9_7;
vtkm::worklet::WaveletCompressor compressor( wname );
vtkm::Id XMaxLevel = compressor.GetWaveletMaxLevel( sigX );
@ -264,23 +289,23 @@ void TestDecomposeReconstruct2D()
std::cout << "Decomposition levels = " << nLevels << std::endl;
std::vector<vtkm::Id> L;
vtkm::Float64 computationTime = 0.0;
vtkm::Float64 elapsedTime = 0.0;
vtkm::Float64 elapsedTime1, elapsedTime2, elapsedTime3;
// Decompose
vtkm::cont::Timer<> timer;
computationTime =
compressor.WaveDecompose2D( inputArray, nLevels, sigX, sigY, outputArray, L,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
elapsedTime = timer.GetElapsedTime();
std::cout << "Decompose time = " << elapsedTime << std::endl;
elapsedTime1 = timer.GetElapsedTime();
std::cout << "Decompose time = " << elapsedTime1 << std::endl;
std::cout << " ->computation time = " << computationTime << std::endl;
// Squash small coefficients
timer.Reset();
vtkm::Float64 cratio = 1.0;
vtkm::Float64 cratio = 100.0;
compressor.SquashCoefficients( outputArray, cratio, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
elapsedTime = timer.GetElapsedTime();
std::cout << "Squash time = " << elapsedTime << std::endl;
elapsedTime2 = timer.GetElapsedTime();
std::cout << "Squash time = " << elapsedTime2 << std::endl;
// Reconstruct
vtkm::cont::ArrayHandle<vtkm::Float64> reconstructArray;
@ -288,12 +313,17 @@ void TestDecomposeReconstruct2D()
computationTime =
compressor.WaveReconstruct2D( outputArray, nLevels, sigX, sigY, reconstructArray, L,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
elapsedTime = timer.GetElapsedTime();
std::cout << "Reconstruction time = " << elapsedTime << std::endl;
elapsedTime3 = timer.GetElapsedTime();
std::cout << "Reconstruction time = " << elapsedTime3 << std::endl;
std::cout << " ->computation time = " << computationTime << std::endl;
std::cout << "Total time = "
<< (elapsedTime1 + elapsedTime2 + elapsedTime3) << std::endl;
outputArray.ReleaseResources();
compressor.EvaluateReconstruction( inputArray, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
/*
timer.Reset();
for( vtkm::Id i = 0; i < reconstructArray.GetNumberOfValues(); i++ )
{
@ -301,8 +331,9 @@ void TestDecomposeReconstruct2D()
inputArray.GetPortalConstControl().Get(i) ),
"output value not the same..." );
}
elapsedTime = timer.GetElapsedTime();
std::cout << "Verification time = " << elapsedTime << std::endl;
elapsedTime1 = timer.GetElapsedTime();
std::cout << "Verification time = " << elapsedTime1 << std::endl;
*/
}
@ -365,8 +396,8 @@ void TestWaveletCompressor()
//DebugDWTIDWT1D();
//DebugRectangleCopy();
//TestDecomposeReconstruct1D();
//TestDecomposeReconstruct2D();
DebugDWT2D();
TestDecomposeReconstruct2D();
//DebugDWT2D();
//DebugExtend2D();
}

80
vtkm/worklet/wavelets/WaveletDWT.h
View File

@ -94,17 +94,16 @@ public:
vtkm::Id extDimX, extDimY;
vtkm::worklet::wavelets::ExtensionDirection2D dir;
// Work on left/top extension
{
{ // Work on left/top extension
if( modeLR )
{
dir = vtkm::worklet::wavelets::ExtensionDirection2D::LEFT;
dir = LEFT;
extDimX = addLen;
extDimY = sigDimY;
}
else
{
dir = vtkm::worklet::wavelets::ExtensionDirection2D::TOP;
dir = TOP;
extDimX = sigDimX;
extDimY = addLen;
}
@ -120,13 +119,13 @@ public:
{
if( modeLR )
{
dir = vtkm::worklet::wavelets::ExtensionDirection2D::RIGHT;
dir = RIGHT;
extDimX = addLen;
extDimY = sigDimY;
}
else
{
dir = vtkm::worklet::wavelets::ExtensionDirection2D::BOTTOM;
dir = BOTTOM;
extDimX = sigDimX;
extDimY = addLen;
}
@ -140,13 +139,13 @@ public:
{
if( modeLR )
{
dir = vtkm::worklet::wavelets::ExtensionDirection2D::RIGHT;
dir = RIGHT;
extDimX = addLen+1;
extDimY = sigDimY;
}
else
{
dir = vtkm::worklet::wavelets::ExtensionDirection2D::BOTTOM;
dir = BOTTOM;
extDimX = sigDimX;
extDimY = addLen+1;
}
@ -167,13 +166,13 @@ public:
ExtendArrayType ext2Temp;
if( modeLR )
{
dir = vtkm::worklet::wavelets::ExtensionDirection2D::RIGHT;
dir = RIGHT;
extDimX = addLen;
extDimY = sigDimY;
}
else
{
dir = vtkm::worklet::wavelets::ExtensionDirection2D::BOTTOM;
dir = BOTTOM;
extDimX = sigDimX;
extDimY = addLen;
}
@ -202,46 +201,6 @@ public:
0, DeviceTag() );
}
}
#if 0
{
vtkm::Id extDimXRight = extDimX;
if( attachZeroRightRight )
extDimXRight++;
rightExt.PrepareForOutput( extDimXRight * extDimY, DeviceTag() );
ExtensionWorklet worklet( extDimXRight, extDimY, sigDimX, sigDimY, rightExtMethod,
vtkm::worklet::wavelets::ExtensionDirection2D::RIGHT,
false );
vtkm::worklet::DispatcherMapField< ExtensionWorklet, DeviceTag >
dispatcher( worklet );
dispatcher.Invoke( rightExt, sigIn );
if( attachZeroRightRight )
WaveletBase::DeviceAssignZero2DColumn( rightExt, extDimXRight, extDimY,
extDimXRight-1, DeviceTag() );
}
else // attachZeroRightLeft mode.
{
ExtendArrayType rightExtendTmp;
rightExtendTmp.PrepareForOutput( extDimX * extDimY, DeviceTag() );
ExtensionWorklet worklet( extDimX, extDimY, sigDimX, sigDimY, rightExtMethod,
vtkm::worklet::wavelets::ExtensionDirection2D::RIGHT,
true ); // treat sigIn as having zeros
vtkm::worklet::DispatcherMapField< ExtensionWorklet, DeviceTag >
dispatcher( worklet );
dispatcher.Invoke( rightExt, sigIn );
rightExt.PrepareForOutput( (extDimX+1) * extDimY, DeviceTag() );
WaveletBase::DeviceRectangleCopyTo( rightExtendTmp, extDimX, extDimY,
rightExt, extDimX+1, extDimY,
1, 0, DeviceTag() );
WaveletBase::DeviceAssignZero2DColumn( rightExt, extDimX+1, extDimY,
0, DeviceTag() );
}
#endif
return 0;
}
@ -918,6 +877,8 @@ if( print)
ArrayType sigExtended;
vtkm::Id outDimX = sigDimX;
vtkm::Id outDimY = sigDimY;
vtkm::cont::Timer<DeviceTag> timer;
vtkm::Float64 computationTime = 0.0;
// First transform on rows
ArrayType afterX;
@ -927,14 +888,16 @@ if( print)
this->Extend2Dv3( sigIn, sigDimX, sigDimY, leftExt, rightExt, addLen,
WaveletBase::wmode, WaveletBase::wmode, false, false,
true, DeviceTag() ); // Extend in left-right direction
timer.Reset();
ForwardXFormv3 worklet( WaveletBase::filter.GetLowDecomposeFilter(),
WaveletBase::filter.GetHighDecomposeFilter(),
filterLen, L[0], oddLow, true, // left-right
addLen, sigDimY, sigDimX, sigDimY, addLen, sigDimY );
DispatcherType dispatcher(worklet);
dispatcher.Invoke( leftExt, sigIn, rightExt, afterX );
//WaveletBase::Print2DArray("\ntest results afrer DWT on rows:", afterX, sigDimX );
computationTime += timer.GetElapsedTime();
}
//WaveletBase::Print2DArray("\ntest results afrer DWT on rows: ", afterX, sigDimX );
// Then do transform in Y direction
coeffOut.PrepareForOutput( sigDimX * sigDimY, DeviceTag() );
@ -943,17 +906,16 @@ if( print)
this->Extend2Dv3( afterX, sigDimX, sigDimY, topExt, bottomExt, addLen,
WaveletBase::wmode, WaveletBase::wmode, false, false,
false, DeviceTag() ); // Extend in top-down direction
//WaveletBase::Print2DArray("top ext:", topExt, sigDimX );
//WaveletBase::Print2DArray("signal:", afterX, sigDimX );
//WaveletBase::Print2DArray("bottom ext:", bottomExt, sigDimX );
timer.Reset();
ForwardXFormv3 worklet( WaveletBase::filter.GetLowDecomposeFilter(),
WaveletBase::filter.GetHighDecomposeFilter(),
filterLen, L[1], oddLow, false, // top-down
sigDimX, addLen, sigDimX, sigDimY, sigDimX, addLen );
DispatcherType dispatcher( worklet );
dispatcher.Invoke( topExt, afterX, bottomExt, coeffOut );
computationTime += timer.GetElapsedTime();
}
return 0;
return computationTime;
}
@ -973,6 +935,8 @@ if( print)
typedef vtkm::worklet::wavelets::InverseTransform2D<DeviceTag> IDWT2DWorklet;
typedef vtkm::worklet::DispatcherMapField<IDWT2DWorklet, DeviceTag> Dispatcher;
BasicArrayType afterY;
vtkm::cont::Timer<DeviceTag> timer;
vtkm::Float64 computationTime = 0.0;
// First inverse transform on columns
{
@ -983,6 +947,7 @@ if( print)
ext1DimY, ext2DimY, ext3DimY, ext4DimY,
filterLen, wmode, DeviceTag() );
afterY.PrepareForOutput( inDimX * inDimY, DeviceTag() );
timer.Reset();
IDWT2DWorklet worklet( WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
filterLen,
@ -995,6 +960,7 @@ if( print)
false ); // top-down
Dispatcher dispatcher( worklet );
dispatcher.Invoke( ext1, ext2, ext3, ext4, coeffIn, afterY );
computationTime += timer.GetElapsedTime();
}
// Then inverse transform on rows
@ -1006,6 +972,7 @@ if( print)
ext1DimX, ext2DimX, ext3DimX, ext4DimX,
filterLen, wmode, DeviceTag() );
sigOut.PrepareForOutput( inDimX * inDimY, DeviceTag() );
timer.Reset();
IDWT2DWorklet worklet( WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
filterLen,
@ -1018,9 +985,10 @@ if( print)
true ); // left-right
Dispatcher dispatcher( worklet );
dispatcher.Invoke( ext1, ext2, ext3, ext4, afterY, sigOut );
computationTime += timer.GetElapsedTime();
}
return 0;
return computationTime;
}

130
vtkm/worklet/wavelets/WaveletTransforms.h
View File

@ -87,6 +87,7 @@ public:
}
}
VTKM_EXEC_CONT_EXPORT
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
{
@ -179,6 +180,7 @@ public:
x2(x_2), y2(y_2),
x3(x_3), y3(y_3), mode_lr(mode) {}
VTKM_EXEC_CONT_EXPORT
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
{
@ -295,7 +297,7 @@ public:
template <typename InPortalType1, typename InPortalType2,
typename InPortalType3, typename OutputPortalType>
VTKM_EXEC_EXPORT
VTKM_EXEC_CONT_EXPORT
void operator()(const InPortalType1 &inPortal1, // left/top extension
const InPortalType2 &inPortal2, // signal
const InPortalType3 &inPortal3, // right/bottom extension
@ -393,6 +395,117 @@ private:
};
// Worklet for 2D signal extension
// This implementation operates on a small rectangle. Use it after LANL.
#if 0
class ExtensionWorklet2Dv3 : 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_EXPORT
ExtensionWorklet2Dv3( vtkm::Id x1, vtkm::Id y1,
vtkm::Id x2, vtkm::Id y2,
vtkm::Id x3, vtkm::Id y3,
vtkm::Id x4, vtkm::Id y4,
DWTMode m, ExtensionDirection2D dir, bool pad_zero)
: extDimX( x1 ), extDimY( y1 ),
sigDimX( x2 ), sigDimY( y2 ),
sigPretendDimX( x3 ), sigPretendDimY( y3 ),
sigStartDimX( x4 ), sigStartDimY( y4 ),
mode(m), direction( dir ), padZero( pad_zero ) {}
// Index translation helper
VTKM_EXEC_CONT_EXPORT
void Ext1Dto2D ( vtkm::Id idx, vtkm::Id &x, vtkm::Id &y ) const
{
x = idx % extDimX;
y = idx / extDimX;
}
// Index translation helper
VTKM_EXEC_CONT_EXPORT
vtkm::Id Sig2Dto1D( vtkm::Id x, vtkm::Id y ) const
{
return y * sigDimX + x;
}
// Index translation helper
VTKM_EXEC_CONT_EXPORT
vtkm::Id SigPretend2Dto1D( vtkm::Id x, vtkm::Id y ) const
{
return (y + sigStartY) * sigDimX + x + sigStartX;
}
template< typename PortalOutType, typename PortalInType >
VTKM_EXEC_EXPORT
void operator()( PortalOutType &portalOut,
const PortalInType &portalIn,
const vtkm::Id &workIndex) const
{
vtkm::Id extX, extY, sigPretendX, sigPretendY;
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 // direction == BOTTOM
{
sigPretendX = extX;
if( mode == SYMH || mode == ASYMH )
sigPretendY = sigPretendDimY - extY - 1;
else
sigPretendY = sigPretendDimY - extY - 2;
if( padZero )
sigPretendY++;
}
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 sigPretendDimX, sigPretendDimY, sigStartX, sigStartY;
const DWTMode mode;
const ExtensionDirection2D direction;
const bool padZero; // treat sigIn as having a column/row zeros
};
#endif
// !!! useful code above !!!
// Worklet for 2D signal extension
class ExtensionWorklet2D : public vtkm::worklet::WorkletMapField
{
@ -1874,11 +1987,6 @@ public:
vtkm::Id outX, outY;
Output1Dto2D( workIdx, outX, outY );
bool print = false;
/*if( outX == 8 && outY == 9 )
{
print = true;
std::cout << "\ntrue convolution: " << std::endl;
}*/
vtkm::Id inX = outX;
vtkm::Id inY = outY;
@ -1897,32 +2005,22 @@ public:
vtkm::Id xi = (inX + 1) / 2;
vtkm::Id sigIdx1D;
if(print)
std::cout << "low filter convolution: " << std::endl;
while( k1 > -1 )
{
sigIdx1D = Input2Dto1D( xi, inY );
sum += lowFilter.Get(k1) * MAKEVAL( sigIn.Get( sigIdx1D ) );
xi++;
k1 -= 2;
if( print )
std::cout << sigIn.Get( sigIdx1D ) << std::endl;
}
xi = inX / 2;
if(print)
std::cout << "high filter convolution: " << std::endl;
while( k2 > -1 )
{
sigIdx1D = Input2Dto1D( xi + cALenExtended, inY );
sum += highFilter.Get(k2) * MAKEVAL( sigIn.Get( sigIdx1D ) );
xi++;
k2 -= 2;
if( print )
std::cout << sigIn.Get( sigIdx1D ) << std::endl;
}
coeffOut = static_cast< OutputValueType> (sum);
if(print)
std::cout << "write back: " << sum << std::endl;
}
#undef MAKEVAL