2016-07-27 20:28:38 +00:00
|
|
|
//============================================================================
|
|
|
|
// 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_waveletbase_h
|
|
|
|
#define vtk_m_worklet_wavelets_waveletbase_h
|
|
|
|
|
|
|
|
|
|
|
|
#include <vtkm/worklet/wavelets/WaveletFilter.h>
|
2016-07-29 04:33:47 +00:00
|
|
|
#include <vtkm/worklet/wavelets/WaveletTransforms.h>
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
#include <vtkm/Math.h>
|
2016-07-27 23:23:26 +00:00
|
|
|
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
namespace vtkm {
|
|
|
|
namespace worklet {
|
|
|
|
|
|
|
|
namespace wavelets {
|
|
|
|
|
|
|
|
// Functionalities are similar to MatWaveBase in VAPoR.
|
|
|
|
class WaveletBase
|
|
|
|
{
|
|
|
|
public:
|
|
|
|
|
|
|
|
// Constructor
|
2016-08-11 16:39:20 +00:00
|
|
|
WaveletBase( WaveletName name ) : wname ( name ),
|
|
|
|
filter( name )
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2016-08-15 19:46:35 +00:00
|
|
|
if( wname == CDF9_7 || wname == BIOR4_4 ||
|
|
|
|
wname == CDF5_3 || wname == BIOR2_2 )
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
|
|
|
this->wmode = SYMW; // Default extension mode, see MatWaveBase.cpp
|
|
|
|
}
|
2016-08-16 21:20:19 +00:00
|
|
|
else if( wname == HAAR || wname == BIOR1_1 ||
|
|
|
|
wname == CDF8_4 || wname == BIOR3_3 )
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2016-08-15 19:46:35 +00:00
|
|
|
this->wmode = SYMH;
|
2016-07-27 20:28:38 +00:00
|
|
|
}
|
2016-07-29 04:33:47 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
// Returns length of approximation coefficients from a decompostition pass.
|
|
|
|
vtkm::Id GetApproxLength( vtkm::Id sigInLen )
|
|
|
|
{
|
2016-08-10 23:24:22 +00:00
|
|
|
vtkm::Id filterLen = this->filter.GetFilterLength();
|
2016-07-29 04:33:47 +00:00
|
|
|
|
2016-08-15 19:46:35 +00:00
|
|
|
if (this->filter.isSymmetric())
|
2016-07-29 04:33:47 +00:00
|
|
|
{
|
|
|
|
if ( (this->wmode == SYMW && (filterLen % 2 != 0)) ||
|
|
|
|
(this->wmode == SYMH && (filterLen % 2 == 0)) )
|
|
|
|
{
|
|
|
|
if (sigInLen % 2 != 0)
|
|
|
|
return((sigInLen+1) / 2);
|
|
|
|
else
|
|
|
|
return((sigInLen) / 2);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2016-08-22 16:13:45 +00:00
|
|
|
return ((sigInLen + filterLen - 1) / 2);
|
2016-07-29 04:33:47 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
// Returns length of detail coefficients from a decompostition pass
|
|
|
|
vtkm::Id GetDetailLength( vtkm::Id sigInLen )
|
|
|
|
{
|
2016-08-10 23:24:22 +00:00
|
|
|
vtkm::Id filterLen = this->filter.GetFilterLength();
|
2016-07-29 04:33:47 +00:00
|
|
|
|
2016-08-15 19:46:35 +00:00
|
|
|
if (this->filter.isSymmetric())
|
2016-07-29 04:33:47 +00:00
|
|
|
{
|
|
|
|
if ( (this->wmode == SYMW && (filterLen % 2 != 0)) ||
|
|
|
|
(this->wmode == SYMH && (filterLen % 2 == 0)) )
|
|
|
|
{
|
|
|
|
if (sigInLen % 2 != 0)
|
|
|
|
return((sigInLen-1) / 2);
|
|
|
|
else
|
|
|
|
return((sigInLen) / 2);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
return static_cast<vtkm::Id>( vtkm::Floor(
|
|
|
|
static_cast<vtkm::Float64>(sigInLen + filterLen - 1) / 2.0 ) );
|
|
|
|
}
|
|
|
|
|
|
|
|
// Returns length of coefficients generated in a decompostition pass
|
|
|
|
vtkm::Id GetCoeffLength( vtkm::Id sigInLen )
|
|
|
|
{
|
|
|
|
return( GetApproxLength( sigInLen ) + GetDetailLength( sigInLen ) );
|
|
|
|
}
|
|
|
|
vtkm::Id GetCoeffLength2( vtkm::Id sigInX, vtkm::Id sigInY )
|
|
|
|
{
|
|
|
|
return( GetCoeffLength( sigInX) * GetCoeffLength( sigInY ) );
|
|
|
|
}
|
|
|
|
vtkm::Id GetCoeffLength3( vtkm::Id sigInX, vtkm::Id sigInY, vtkm::Id sigInZ)
|
|
|
|
{
|
|
|
|
return( GetCoeffLength( sigInX) * GetCoeffLength( sigInY ) * GetCoeffLength( sigInZ ) );
|
|
|
|
}
|
|
|
|
|
|
|
|
// Returns maximum wavelet decompostion level
|
|
|
|
vtkm::Id GetWaveletMaxLevel( vtkm::Id sigInLen )
|
|
|
|
{
|
2016-08-10 23:24:22 +00:00
|
|
|
vtkm::Id filterLen = this->filter.GetFilterLength();
|
|
|
|
vtkm::Id level;
|
|
|
|
this->WaveLengthValidate( sigInLen, filterLen, level );
|
|
|
|
return level;
|
2016-07-27 23:23:26 +00:00
|
|
|
}
|
2016-07-29 01:12:57 +00:00
|
|
|
|
2016-08-03 15:34:04 +00:00
|
|
|
// perform a device copy. The whole 1st array to a certain start location of the 2nd array
|
2016-08-22 20:49:13 +00:00
|
|
|
template< typename ArrayType1, typename ArrayType2, typename DeviceTag >
|
2016-08-03 15:34:04 +00:00
|
|
|
void DeviceCopyStartX( const ArrayType1 &srcArray,
|
|
|
|
ArrayType2 &dstArray,
|
2016-08-22 20:49:13 +00:00
|
|
|
vtkm::Id startIdx,
|
|
|
|
DeviceTag )
|
2016-08-03 15:34:04 +00:00
|
|
|
{
|
|
|
|
typedef vtkm::worklet::wavelets::CopyWorklet CopyType;
|
|
|
|
CopyType cp( startIdx );
|
2016-08-22 20:49:13 +00:00
|
|
|
vtkm::worklet::DispatcherMapField< CopyType, DeviceTag > dispatcher( cp );
|
2016-08-03 15:34:04 +00:00
|
|
|
dispatcher.Invoke( srcArray, dstArray );
|
|
|
|
}
|
|
|
|
|
2016-08-15 19:46:35 +00:00
|
|
|
// Assign zero value to a certain location of an array
|
2016-08-22 20:49:13 +00:00
|
|
|
template< typename ArrayType, typename DeviceTag >
|
|
|
|
void DeviceAssignZero( ArrayType &array, vtkm::Id index, DeviceTag )
|
2016-08-15 19:46:35 +00:00
|
|
|
{
|
|
|
|
typedef vtkm::worklet::wavelets::AssignZeroWorklet ZeroWorklet;
|
|
|
|
ZeroWorklet worklet( index );
|
2016-08-22 20:49:13 +00:00
|
|
|
vtkm::worklet::DispatcherMapField< ZeroWorklet, DeviceTag > dispatcher( worklet );
|
2016-08-15 19:46:35 +00:00
|
|
|
dispatcher.Invoke( array );
|
|
|
|
}
|
|
|
|
|
2016-08-30 00:20:49 +00:00
|
|
|
// Assign zeros to a certain column to a matrix
|
|
|
|
template< typename ArrayType, typename DeviceTag >
|
|
|
|
void DeviceAssignZero2DColumn( ArrayType &array, vtkm::Id dimX, vtkm::Id dimY, // input
|
|
|
|
vtkm::Id colIdx, DeviceTag )
|
|
|
|
{
|
|
|
|
typedef vtkm::worklet::wavelets::AssignZero2DColumnWorklet AssignZero2DType;
|
|
|
|
AssignZero2DType zeroWorklet( dimX, dimY, colIdx );
|
|
|
|
vtkm::worklet::DispatcherMapField< AssignZero2DType, DeviceTag >
|
|
|
|
dispatcher( zeroWorklet );
|
|
|
|
dispatcher.Invoke( array );
|
|
|
|
}
|
|
|
|
|
2016-07-29 01:12:57 +00:00
|
|
|
// Sort by the absolute value on device
|
2016-07-29 03:07:05 +00:00
|
|
|
struct SortLessAbsFunctor
|
2016-07-29 01:12:57 +00:00
|
|
|
{
|
|
|
|
template< typename T >
|
2016-08-15 19:46:35 +00:00
|
|
|
VTKM_EXEC_EXPORT
|
2016-07-29 01:12:57 +00:00
|
|
|
bool operator()(const T& x, const T& y) const
|
|
|
|
{
|
|
|
|
return vtkm::Abs(x) < vtkm::Abs(y);
|
|
|
|
}
|
|
|
|
};
|
2016-08-10 20:33:14 +00:00
|
|
|
template< typename ArrayType, typename DeviceTag >
|
|
|
|
void DeviceSort( ArrayType &array, DeviceTag )
|
2016-07-28 22:40:37 +00:00
|
|
|
{
|
2016-08-10 20:33:14 +00:00
|
|
|
vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Sort
|
2016-07-29 03:07:05 +00:00
|
|
|
( array, SortLessAbsFunctor() );
|
2016-07-29 01:12:57 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
// Reduce to the sum of all values on device
|
2016-08-10 20:33:14 +00:00
|
|
|
template< typename ArrayType, typename DeviceTag >
|
|
|
|
typename ArrayType::ValueType DeviceSum( const ArrayType &array, DeviceTag )
|
2016-07-29 01:12:57 +00:00
|
|
|
{
|
2016-08-10 20:33:14 +00:00
|
|
|
return vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Reduce
|
2016-07-29 03:07:05 +00:00
|
|
|
( array, 0.0 );
|
|
|
|
}
|
|
|
|
|
|
|
|
// Find the max and min of an array
|
|
|
|
struct minFunctor
|
|
|
|
{
|
|
|
|
template< typename FieldType >
|
2016-08-15 19:46:35 +00:00
|
|
|
VTKM_EXEC_EXPORT
|
2016-07-29 03:07:05 +00:00
|
|
|
FieldType operator()(const FieldType &x, const FieldType &y) const {
|
|
|
|
return Min(x, y);
|
|
|
|
}
|
|
|
|
};
|
|
|
|
struct maxFunctor
|
|
|
|
{
|
|
|
|
template< typename FieldType >
|
2016-08-15 19:46:35 +00:00
|
|
|
VTKM_EXEC_EXPORT
|
2016-07-29 03:07:05 +00:00
|
|
|
FieldType operator()(const FieldType& x, const FieldType& y) const {
|
2016-08-22 20:49:13 +00:00
|
|
|
return vtkm::Max(x, y);
|
2016-07-29 03:07:05 +00:00
|
|
|
}
|
|
|
|
};
|
2016-08-10 20:33:14 +00:00
|
|
|
template< typename ArrayType, typename DeviceTag >
|
|
|
|
typename ArrayType::ValueType DeviceMax( const ArrayType &array, DeviceTag )
|
2016-07-29 03:07:05 +00:00
|
|
|
{
|
|
|
|
typename ArrayType::ValueType initVal = array.GetPortalConstControl().Get(0);
|
2016-08-10 20:33:14 +00:00
|
|
|
return vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Reduce
|
2016-07-29 03:07:05 +00:00
|
|
|
( array, initVal, maxFunctor() );
|
|
|
|
}
|
2016-08-10 20:33:14 +00:00
|
|
|
template< typename ArrayType, typename DeviceTag >
|
|
|
|
typename ArrayType::ValueType DeviceMin( const ArrayType &array, DeviceTag )
|
2016-07-29 03:07:05 +00:00
|
|
|
{
|
|
|
|
typename ArrayType::ValueType initVal = array.GetPortalConstControl().Get(0);
|
2016-08-10 20:33:14 +00:00
|
|
|
return vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Reduce
|
2016-07-29 03:07:05 +00:00
|
|
|
( array, initVal, minFunctor() );
|
|
|
|
}
|
|
|
|
|
2016-07-29 04:33:47 +00:00
|
|
|
// Max absolute value of an array
|
|
|
|
struct maxAbsFunctor
|
|
|
|
{
|
|
|
|
template< typename FieldType >
|
2016-08-15 19:46:35 +00:00
|
|
|
VTKM_EXEC_EXPORT
|
2016-07-29 04:33:47 +00:00
|
|
|
FieldType operator()(const FieldType& x, const FieldType& y) const {
|
2016-08-22 20:49:13 +00:00
|
|
|
return vtkm::Max( vtkm::Abs(x), vtkm::Abs(y) );
|
2016-07-29 04:33:47 +00:00
|
|
|
}
|
|
|
|
};
|
2016-08-10 20:33:14 +00:00
|
|
|
template< typename ArrayType, typename DeviceTag >
|
|
|
|
typename ArrayType::ValueType DeviceMaxAbs( const ArrayType &array, DeviceTag )
|
2016-07-29 04:33:47 +00:00
|
|
|
{
|
|
|
|
typename ArrayType::ValueType initVal = array.GetPortalConstControl().Get(0);
|
2016-08-10 20:33:14 +00:00
|
|
|
return vtkm::cont::DeviceAdapterAlgorithm< DeviceTag >::Reduce
|
2016-07-29 04:33:47 +00:00
|
|
|
( array, initVal, maxAbsFunctor() );
|
|
|
|
}
|
2016-07-27 23:23:26 +00:00
|
|
|
|
2016-07-29 04:33:47 +00:00
|
|
|
// Calculate variance of an array
|
2016-08-10 20:33:14 +00:00
|
|
|
template< typename ArrayType, typename DeviceTag >
|
|
|
|
vtkm::Float64 DeviceCalculateVariance( ArrayType &array, DeviceTag )
|
2016-07-27 20:28:38 +00:00
|
|
|
{
|
2016-08-10 20:33:14 +00:00
|
|
|
vtkm::Float64 mean = static_cast<vtkm::Float64>(this->DeviceSum( array, DeviceTag() )) /
|
2016-07-29 04:33:47 +00:00
|
|
|
static_cast<vtkm::Float64>(array.GetNumberOfValues());
|
|
|
|
|
|
|
|
vtkm::cont::ArrayHandle< vtkm::Float64 > squaredDeviation;
|
|
|
|
|
|
|
|
// Use a worklet
|
|
|
|
typedef vtkm::worklet::wavelets::SquaredDeviation SDWorklet;
|
|
|
|
SDWorklet sdw( mean );
|
2016-08-22 20:49:13 +00:00
|
|
|
vtkm::worklet::DispatcherMapField< SDWorklet, DeviceTag > dispatcher( sdw );
|
2016-07-29 04:33:47 +00:00
|
|
|
dispatcher.Invoke( array, squaredDeviation );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-08-10 20:33:14 +00:00
|
|
|
vtkm::Float64 sdMean = this->DeviceSum( squaredDeviation, DeviceTag() ) /
|
2016-07-29 04:33:47 +00:00
|
|
|
static_cast<vtkm::Float64>( squaredDeviation.GetNumberOfValues() );
|
2016-07-27 20:28:38 +00:00
|
|
|
|
2016-07-29 04:33:47 +00:00
|
|
|
return sdMean;
|
2016-07-27 20:28:38 +00:00
|
|
|
}
|
|
|
|
|
2016-08-09 23:02:10 +00:00
|
|
|
// Transpose a matrix in an array
|
2016-08-22 20:49:13 +00:00
|
|
|
template< typename InputArrayType, typename OutputArrayType, typename DeviceTag >
|
2016-08-09 23:02:10 +00:00
|
|
|
void DeviceTranspose( const InputArrayType &inputArray,
|
|
|
|
OutputArrayType &outputArray,
|
|
|
|
vtkm::Id inputX,
|
2016-08-22 20:49:13 +00:00
|
|
|
vtkm::Id inputY,
|
|
|
|
DeviceTag )
|
2016-08-09 23:02:10 +00:00
|
|
|
{
|
|
|
|
// use a worklet
|
|
|
|
typedef vtkm::worklet::wavelets::TransposeWorklet TransposeType;
|
|
|
|
TransposeType tw ( inputX, inputY );
|
2016-08-22 20:49:13 +00:00
|
|
|
vtkm::worklet::DispatcherMapField< TransposeType, DeviceTag > dispatcher( tw );
|
2016-08-09 23:02:10 +00:00
|
|
|
dispatcher.Invoke( inputArray, outputArray );
|
|
|
|
}
|
|
|
|
|
2016-08-14 22:29:27 +00:00
|
|
|
// Copy a small rectangle to a big rectangle
|
2016-08-22 20:49:13 +00:00
|
|
|
template< typename SmallArrayType, typename BigArrayType, typename DeviceTag>
|
2016-08-14 22:29:27 +00:00
|
|
|
void DeviceRectangleCopyTo( const SmallArrayType &smallRect,
|
|
|
|
vtkm::Id smallX,
|
|
|
|
vtkm::Id smallY,
|
|
|
|
BigArrayType &bigRect,
|
|
|
|
vtkm::Id bigX,
|
|
|
|
vtkm::Id bigY,
|
|
|
|
vtkm::Id startX,
|
2016-08-22 20:49:13 +00:00
|
|
|
vtkm::Id startY,
|
|
|
|
DeviceTag )
|
2016-08-14 22:29:27 +00:00
|
|
|
{
|
|
|
|
typedef vtkm::worklet::wavelets::RectangleCopyTo CopyToWorklet;
|
|
|
|
CopyToWorklet cp( smallX, smallY, bigX, bigY, startX, startY );
|
2016-08-22 20:49:13 +00:00
|
|
|
vtkm::worklet::DispatcherMapField< CopyToWorklet, DeviceTag > dispatcher( cp );
|
2016-08-14 22:29:27 +00:00
|
|
|
dispatcher.Invoke(smallRect, bigRect);
|
|
|
|
}
|
|
|
|
|
|
|
|
// Fill a small rectangle from a portion of a big rectangle
|
2016-08-22 20:49:13 +00:00
|
|
|
template< typename SmallArrayType, typename BigArrayType, typename DeviceTag >
|
2016-08-14 22:29:27 +00:00
|
|
|
void DeviceRectangleCopyFrom( SmallArrayType &smallRect,
|
|
|
|
vtkm::Id smallX,
|
|
|
|
vtkm::Id smallY,
|
|
|
|
const BigArrayType &bigRect,
|
|
|
|
vtkm::Id bigX,
|
|
|
|
vtkm::Id bigY,
|
|
|
|
vtkm::Id startX,
|
2016-08-22 20:49:13 +00:00
|
|
|
vtkm::Id startY,
|
|
|
|
DeviceTag )
|
2016-08-14 22:29:27 +00:00
|
|
|
{
|
2016-08-22 20:49:13 +00:00
|
|
|
smallRect.PrepareForOutput( smallX*smallY, DeviceTag() );
|
2016-08-14 22:29:27 +00:00
|
|
|
typedef vtkm::worklet::wavelets::RectangleCopyFrom CopyFromWorklet;
|
|
|
|
CopyFromWorklet cpFrom( smallX, smallY, bigX, bigY, startX, startY );
|
2016-08-22 20:49:13 +00:00
|
|
|
vtkm::worklet::DispatcherMapField< CopyFromWorklet, DeviceTag > dispatcherFrom( cpFrom );
|
2016-08-14 22:29:27 +00:00
|
|
|
dispatcherFrom.Invoke( smallRect, bigRect );
|
|
|
|
}
|
|
|
|
|
|
|
|
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
protected:
|
2016-08-10 23:24:22 +00:00
|
|
|
WaveletName wname;
|
|
|
|
DWTMode wmode;
|
|
|
|
WaveletFilter filter;
|
2016-07-27 20:28:38 +00:00
|
|
|
|
|
|
|
void WaveLengthValidate( vtkm::Id sigInLen, vtkm::Id filterLength, vtkm::Id &level)
|
|
|
|
{
|
|
|
|
if( sigInLen < filterLength )
|
|
|
|
level = 0;
|
|
|
|
else
|
|
|
|
level = static_cast<vtkm::Id>( vtkm::Floor(
|
|
|
|
vtkm::Log2( static_cast<vtkm::Float64>(sigInLen) /
|
|
|
|
static_cast<vtkm::Float64>(filterLength) ) + 1.0 ) );
|
|
|
|
}
|
|
|
|
}; // class WaveletBase.
|
|
|
|
|
|
|
|
|
|
|
|
} // namespace wavelets
|
|
|
|
|
|
|
|
} // namespace worklet
|
|
|
|
} // namespace vtkm
|
|
|
|
|
|
|
|
#endif
|