Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-20 02:55:47 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

finish 2D re-write. Some tests pass and some others fail. Need to keep debugging

Browse Source
This commit is contained in:
Samuel Li 2016-09-05 19:52:35 -06:00
parent 48588c3852
commit 42a80662d8
4 changed files with 236 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
vtkm/worklet/WaveletCompressor.h
View File

@ -145,7 +145,7 @@ public:
// Make an output array
OutArrayBasic output;
WaveletDWT::IDWT1D( input, L1d, output, DeviceTag() );
WaveletDWT::IDWT1D( input, L1d, output, false, DeviceTag() );
VTKM_ASSERT( output.GetNumberOfValues() == L1d[2] );
// Move output to intermediate array
@ -207,7 +207,7 @@ public:
OutBasicArray tempOutput;
computationTime +=
WaveletDWT::DWT2D( tempInput, currentLenX, currentLenY, tempOutput, L2d, DeviceTag());
WaveletDWT::DWT2Dv2( tempInput, currentLenX, currentLenY, tempOutput, L2d, DeviceTag());
// copy results to coeffOut
WaveletBase::DeviceRectangleCopyTo( tempOutput, currentLenX, currentLenY,
@ -274,7 +274,7 @@ public:
// IDWT
computationTime +=
WaveletDWT::IDWT2D( tempInput, L2d, tempOutput, DeviceTag() );
WaveletDWT::IDWT2Dv2( tempInput, L2d, tempOutput, DeviceTag() );
// copy back reconstructed block
WaveletBase::DeviceRectangleCopyTo( tempOutput, L2d[8], L2d[9],

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

@ -59,10 +59,11 @@ void FillArray( ArrayType& array )
dispatcher.Invoke( array );
}
void Debug2DExtend()
{
vtkm::Id NX = 10;
vtkm::Id NY = 11;
vtkm::Id NY = 8;
typedef vtkm::cont::ArrayHandle< vtkm::Float64 > ArrayType;
ArrayType left, center, right;
@ -73,10 +74,12 @@ void Debug2DExtend()
ArrayType output1, output2;
std::vector<vtkm::Id> L(10, 0);
vtkm::worklet::wavelets::WaveletDWT dwt( vtkm::worklet::wavelets::CDF9_7 );
vtkm::worklet::wavelets::WaveletDWT dwt( vtkm::worklet::wavelets::CDF8_4 );
// get true results
dwt.DWT2D(center, NX, NY, output1, L, VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
// get test results
dwt.DWT2Dv2( center, NX, NY, output2, L, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
for( vtkm::Id i = 0; i < output1.GetNumberOfValues(); i++ )
@ -86,7 +89,7 @@ void Debug2DExtend()
"WaveletCompressor worklet failed..." );
}
std::cout << "true results:" << std::endl;
std::cout << "true results after DWT:" << std::endl;
for( vtkm::Id i = 0; i < output1.GetNumberOfValues(); i++ )
{
std::cout << output1.GetPortalConstControl().Get(i) << " ";
@ -94,13 +97,39 @@ std::cout << "true results:" << std::endl;
std::cout << std::endl;
}
std::cout << "test results:" << std::endl;
std::cout << "test results after DWT:" << std::endl;
for( vtkm::Id i = 0; i < output2.GetNumberOfValues(); i++ )
{
std::cout << output2.GetPortalConstControl().Get(i) << " ";
if( i % NX == NX - 1 )
std::cout << std::endl;
}
ArrayType idwt_out1, idwt_out2;
// true results go through IDWT
dwt.IDWT2D( output1, L, idwt_out1, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
// test results go through IDWT
dwt.IDWT2Dv2( output2, L, idwt_out2, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
std::cout << "true results after IDWT:" << std::endl;
for( vtkm::Id i = 0; i < idwt_out1.GetNumberOfValues(); i++ )
{
std::cout << idwt_out1.GetPortalConstControl().Get(i) << " ";
if( i % NX == NX - 1 )
std::cout << std::endl;
}
std::cout << "test results after IDWT:" << std::endl;
for( vtkm::Id i = 0; i < idwt_out2.GetNumberOfValues(); i++ )
{
std::cout << idwt_out2.GetPortalConstControl().Get(i) << " ";
if( i % NX == NX - 1 )
std::cout << std::endl;
}
}
void DebugDWTIDWT1D()
@ -139,7 +168,7 @@ void DebugDWTIDWT1D()
// Inverse Transform
vtkm::cont::ArrayHandle<vtkm::Float64> reconstructArray;
waveletdwt.IDWT1D( coeffOut, L, reconstructArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
waveletdwt.IDWT1D( coeffOut, L, reconstructArray, false, VTKM_DEFAULT_DEVICE_ADAPTER_TAG() );
std::cout << "Inverse Wavelet Transform: result signal length = " <<
reconstructArray.GetNumberOfValues() << std::endl;
for( vtkm::Id i = 0; i < reconstructArray.GetNumberOfValues(); i++ )
@ -202,12 +231,13 @@ void TestDecomposeReconstruct2D()
vtkm::cont::ArrayHandle<vtkm::Float64> outputArray;
// Use a WaveletCompressor
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF9_7;
vtkm::worklet::wavelets::WaveletName wname = vtkm::worklet::wavelets::CDF8_4;
vtkm::worklet::WaveletCompressor compressor( wname );
vtkm::Id XMaxLevel = compressor.GetWaveletMaxLevel( sigX );
vtkm::Id YMaxLevel = compressor.GetWaveletMaxLevel( sigY );
vtkm::Id nLevels = vtkm::Min( XMaxLevel, YMaxLevel );
std::cout << "nLevels = " << nLevels << std::endl;
//nLevels = 1;
std::vector<vtkm::Id> L;
@ -242,9 +272,14 @@ void TestDecomposeReconstruct2D()
timer.Reset();
for( vtkm::Id i = 0; i < reconstructArray.GetNumberOfValues(); i++ )
{
VTKM_TEST_ASSERT( test_equal( reconstructArray.GetPortalConstControl().Get(i),
inputArray.GetPortalConstControl().Get(i) ),
"output value not the same..." );
vtkm::Float64 real = inputArray.GetPortalConstControl().Get(i);
vtkm::Float64 rect = reconstructArray.GetPortalConstControl().Get(i);
vtkm::Float64 diff = real - rect;
if( diff < -0.00001 || diff > 0.00001 )
std::cout << i << ": " << real << ", \t" << rect << std::endl;
//VTKM_TEST_ASSERT( test_equal( reconstructArray.GetPortalConstControl().Get(i),
// inputArray.GetPortalConstControl().Get(i) ),
// "output value not the same..." );
}
elapsedTime = timer.GetElapsedTime();
std::cout << "Verification time = " << elapsedTime << std::endl;
@ -310,8 +345,8 @@ void TestWaveletCompressor()
//DebugDWTIDWT1D();
//DebugRectangleCopy();
//TestDecomposeReconstruct1D();
//TestDecomposeReconstruct2D();
Debug2DExtend();
TestDecomposeReconstruct2D();
//Debug2DExtend();
}
int UnitTestWaveletCompressor(int, char *[])

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

@ -573,7 +573,6 @@ if( print)
{
for( vtkm::Id i = 0; i < coeffInExtended.GetNumberOfValues(); i++ )
printf( "%.2e ",coeffInExtended.GetPortalConstControl().Get(i) );
//std::cout << coeffInExtended.GetPortalConstControl().Get(i) << " ";
std::cout << std::endl;
}
@ -891,32 +890,103 @@ if( print)
vtkm::Id inDimY = L[1] + L[3];
VTKM_ASSERT( inDimX * inDimY == coeffIn.GetNumberOfValues() );
std::vector<vtkm::Id> xL(3, 0);
vtkm::Id filterLen = WaveletBase::filter.GetFilterLength();
typedef vtkm::worklet::wavelets::InverseTransform2DOdd IdwtOddWorklet;
typedef vtkm::worklet::wavelets::InverseTransform2DEven IdwtEvenWorklet;
// First inverse transform on columns
// Transpose
ArrayInType beforeY, beforeYExtend;
beforeY.PrepareForOutput( inDimY * inDimX, DeviceTag() );
vtkm::Id beforeYDimX = inDimY;
vtkm::Id beforeYDimY = inDimX;
beforeY.PrepareForOutput( beforeYDimX * beforeYDimY, DeviceTag() );
WaveletBase::DeviceTranspose( coeffIn, inDimX, inDimY,
beforeY, inDimY, inDimX,
beforeY, beforeYDimX, beforeYDimY,
0, 0, DeviceTag() );
vtkm::Id beforeYExtendDimX;
vtkm::Id beforeYExtendDimY = inDimX;
this->IDWTHelper( beforeY, L[1], L[3], beforeYExtendDimY,
beforeYExtend, beforeYExtendDimX,
vtkm::Id cATempLen, beforeYExtendDimX;
vtkm::Id beforeYExtendDimY = beforeYDimY;
this->IDWTHelper( beforeY, L[1], L[3], beforeYDimY,
beforeYExtend, cATempLen, beforeYExtendDimX,
filterLen, wmode, DeviceTag() );
std::cout << "print v2 after extension" << std::endl;
for( vtkm::Id i = 0; i < beforeYExtend.GetNumberOfValues(); i++ )
beforeY.ReleaseResources();
// Inverse transform
ArrayInType afterY;
vtkm::Id afterYDimX = beforeYDimX;
vtkm::Id afterYDimY = beforeYDimY;
afterY.PrepareForOutput( afterYDimX * afterYDimY, DeviceTag() );
if( filterLen % 2 != 0 )
{
printf("%.2e ",beforeYExtend.GetPortalConstControl().Get(i));
//std::cout << beforeYExtend.GetPortalConstControl().Get(i) << " ";
if( i % beforeYExtendDimX == beforeYExtendDimX - 1 )
std::cout << std::endl;
IdwtOddWorklet worklet( filterLen, beforeYExtendDimX, beforeYExtendDimY,
afterYDimX, afterYDimY, cATempLen );
vtkm::worklet::DispatcherMapField<IdwtOddWorklet, DeviceTag>
dispatcher( worklet );
dispatcher.Invoke( beforeYExtend,
WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
afterY );
}
std::cout << "finish printing v2 after extension" << std::endl;
else
{
IdwtEvenWorklet worklet( filterLen, beforeYExtendDimX, beforeYExtendDimY,
afterYDimX, afterYDimY, cATempLen );
vtkm::worklet::DispatcherMapField<IdwtEvenWorklet, DeviceTag>
dispatcher( worklet );
dispatcher.Invoke( beforeYExtend,
WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
afterY );
}
beforeYExtend.ReleaseResources();
// Second inverse transform on rows
ArrayInType beforeX, beforeXExtend;
vtkm::Id beforeXDimX = afterYDimY;
vtkm::Id beforeXDimY = afterYDimX;
vtkm::Id beforeXExtendDimX;
vtkm::Id beforeXExtendDimY = beforeXDimY;
beforeX.PrepareForOutput( beforeXDimX * beforeXDimY, DeviceTag() );
WaveletBase::DeviceTranspose( afterY, afterYDimX, afterYDimY,
beforeX, beforeXDimX, beforeXDimY,
0, 0, DeviceTag() );
this->IDWTHelper( beforeX, L[0], L[4], beforeXDimY,
beforeXExtend, cATempLen, beforeXExtendDimX,
filterLen, wmode, DeviceTag() );
beforeX.ReleaseResources();
// Inverse transform
ArrayInType afterX;
vtkm::Id afterXDimX = beforeXDimX;
vtkm::Id afterXDimY = beforeXDimY;
afterX.PrepareForOutput( afterXDimX * afterXDimY, DeviceTag() );
if( filterLen % 2 != 0 )
{
IdwtOddWorklet worklet( filterLen, beforeXExtendDimX, beforeXExtendDimY,
afterXDimX, afterXDimY, cATempLen );
vtkm::worklet::DispatcherMapField<IdwtOddWorklet, DeviceTag>
dispatcher( worklet );
dispatcher.Invoke( beforeXExtend,
WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
afterX );
}
else
{
IdwtEvenWorklet worklet( filterLen, beforeXExtendDimX, beforeXExtendDimY,
afterXDimX, afterXDimY, cATempLen );
vtkm::worklet::DispatcherMapField<IdwtEvenWorklet, DeviceTag>
dispatcher( worklet );
dispatcher.Invoke( beforeXExtend,
WaveletBase::filter.GetLowReconstructFilter(),
WaveletBase::filter.GetHighReconstructFilter(),
afterX );
}
sigOut = afterX;
return 0;
}
@ -928,6 +998,7 @@ std::cout << "finish printing v2 after extension" << std::endl;
vtkm::Id cDDimX, // of codffIn
vtkm::Id inDimY, // of codffIn
ArrayOutType &coeffExtend,
vtkm::Id &cATempLen, // output
vtkm::Id &coeffExtendDimX, // output
vtkm::Id filterLen,
DWTMode mode,
@ -961,7 +1032,7 @@ std::cout << "finish printing v2 after extension" << std::endl;
cARight = SYMH;
}
// determine length after extension
vtkm::Id cATempLen, cDTempLen; // X dimension after extension
vtkm::Id cDTempLen; // cD dimension after extension. cA is passed in
vtkm::Id cDPadLen = 0;
vtkm::Id addLen = filterLen / 4; // addLen == 0 for Haar kernel
if( (cADimX > cDDimX) && (mode == SYMH) )
@ -993,7 +1064,6 @@ std::cout << "finish printing v2 after extension" << std::endl;
{
this->Extend2D( cD, cDDimX, cDDimY, cDTemp, addLen,
cDLeft, cDRight, true, false, DeviceTag() );
std::cout << "went here" << std::endl;
}
else
{

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

@ -309,12 +309,12 @@ public:
y = idx / outputDimX;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id Input2Dto1D( vtkm::Id &x, vtkm::Id &y ) const
vtkm::Id Input2Dto1D( const vtkm::Id &x, const vtkm::Id &y ) const
{
return y * inputDimX + x;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id Output2Dto1D( vtkm::Id &x, vtkm::Id &y ) const
vtkm::Id Output2Dto1D( const vtkm::Id &x, const vtkm::Id &y ) const
{
return y * outputDimX + x;
}
@ -403,10 +403,9 @@ public:
// Constructor
VTKM_EXEC_CONT_EXPORT
InverseTransform2DOdd( vtkm::Id fil_len, vtkm::Id x1, vtkm::Id y1, vtkm::Id x2,
vtkm::Id y2, vtkm::Id cA_len, vtkm::Id cA_len_ext )
vtkm::Id y2, vtkm::Id cA_len_ext )
: filterLen( fil_len ), inputDimX( x1 ), inputDimY( y1 ),
outputDimX( x2 ), outputDimY( y2 ), cALen( cA_len ),
cALen2( cA_len*2 ), cALenExtended( cA_len_ext ) {}
outputDimX( x2 ), outputDimY( y2 ), cALenExtended( cA_len_ext ) {}
VTKM_EXEC_CONT_EXPORT
void Output1Dto2D( const vtkm::Id &idx, vtkm::Id &x, vtkm::Id &y ) const
@ -415,15 +414,10 @@ public:
y = idx / outputDimX;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id Input2Dto1D( vtkm::Id &x, vtkm::Id &y ) const
vtkm::Id Input2Dto1D( vtkm::Id x, vtkm::Id y ) const
{
return y * inputDimX + x;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id Output2Dto1D( vtkm::Id &x, vtkm::Id &y ) const
{
return y * outputDimX + x;
}
// Use 64-bit float for convolution calculation
#define VAL vtkm::Float64
@ -467,11 +461,9 @@ public:
k1 -= 2;
}
xi = inX / 2;
vtkm::Id realXIdx;
while( k2 > -1 )
{
realXIdx = xi + cALenExtended;
sigIdx1D = Input2Dto1D( realXIdx, inY );
sigIdx1D = Input2Dto1D( xi + cALenExtended, inY );
sum += hiFilter.Get(k2) * MAKEVAL( sigIn.Get( sigIdx1D ) );
xi++;
k2 -= 2;
@ -486,7 +478,6 @@ public:
private:
vtkm::Id filterLen;
vtkm::Id inputDimX, inputDimY, outputDimX, outputDimY;
vtkm::Id cALen, cALen2; // cALen2 == cALen * 2
vtkm::Id cALenExtended; // Number of cA at the beginning of input, followed by cD
};
@ -584,6 +575,98 @@ private:
};
// Worklet: perform an inverse transform for even length, symmetric filters.
class InverseTransform2DEven: public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn<ScalarAll>, // Input: coeffs,
// cA followed by cD
WholeArrayIn<Scalar>, // lowFilter
WholeArrayIn<Scalar>, // highFilter
FieldOut<ScalarAll>); // output
typedef void ExecutionSignature(_1, _2, _3, _4, WorkIndex);
typedef _4 InputDomain;
// Constructor
VTKM_EXEC_CONT_EXPORT
InverseTransform2DEven( vtkm::Id filtL, vtkm::Id x1, vtkm::Id y1,
vtkm::Id x2, vtkm::Id y2, vtkm::Id cALExt ) :
filterLen(filtL), inputDimX( x1 ), inputDimY( y1 ),
outputDimX( x2 ), outputDimY( y2 ), cALenExtended(cALExt) {}
VTKM_EXEC_CONT_EXPORT
void Output1Dto2D( const vtkm::Id &idx, vtkm::Id &x, vtkm::Id &y ) const
{
x = idx % outputDimX;
y = idx / outputDimX;
}
VTKM_EXEC_CONT_EXPORT
vtkm::Id Input2Dto1D( vtkm::Id x, vtkm::Id y ) const
{
return y * inputDimX + x;
}
// Use 64-bit float for convolution calculation
#define VAL vtkm::Float64
#define MAKEVAL(a) (static_cast<VAL>(a))
template <typename InputPortalType,
typename FilterPortalType,
typename OutputValueType>
VTKM_EXEC_EXPORT
void operator()(const InputPortalType &coeffs,
const FilterPortalType &lowFilter,
const FilterPortalType &highFilter,
OutputValueType &sigOut,
const vtkm::Id &workIndex) const
{
vtkm::Id outX, outY;
Output1Dto2D( workIndex, outX, outY );
vtkm::Id inX = outX;
vtkm::Id inY = outY;
vtkm::Id xi, k;
VAL sum = 0.0;
if( (filterLen/2) % 2 != 0 ) // odd length half filter
{
xi = inX / 2;
if( inX % 2 != 0 )
k = filterLen - 1;
else
k = filterLen - 2;
}
else
{
xi = (inX + 1) / 2;
if( inX % 2 != 0 )
k = filterLen - 2;
else
k = filterLen - 1;
}
vtkm::Id cAIdx1D, cDIdx1D;
while( k > -1 ) // k >= 0
{
cAIdx1D = Input2Dto1D( xi, inY );
cDIdx1D = Input2Dto1D( xi + cALenExtended, inY );
sum += lowFilter.Get(k) * MAKEVAL( coeffs.Get( cAIdx1D ) ) + // cA
highFilter.Get(k) * MAKEVAL( coeffs.Get( cDIdx1D ) ); // cD
xi++;
k -= 2;
}
sigOut = static_cast<OutputValueType>( sum );
}
#undef MAKEVAL
#undef VAL
private:
vtkm::Id filterLen; // filter length.
vtkm::Id inputDimX, inputDimY, outputDimX, outputDimY;
vtkm::Id cALenExtended; // Number of cA at the beginning of input, followed by cD
};
// Worklet: perform an inverse transform for even length, symmetric filters.
class InverseTransformEven: public vtkm::worklet::WorkletMapField
{