Commit 3b2fb5c8 authored by Roxana Bujack's avatar Roxana Bujack Committed by Kitware Robot

Merge topic 'realRadiusFix'

9c41daa2 change input to spacing and radius to make it more intuitive
b991d3c6 merge style differences
7d0ee16f add the real double radius for scaling
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Merge-request: !1825
parents 4faf1bf0 9c41daa2
......@@ -28,12 +28,15 @@ public:
const vtkm::filter::FieldMetadata& fieldMetadata,
const vtkm::filter::PolicyBase<DerivedPolicy>&);
VTKM_CONT void SetRadius(vtkm::Vec<vtkm::Int32, 3> _radius) { this->Radius = _radius; }
VTKM_CONT void SetRadius(double _radius) { this->Radius = _radius; }
VTKM_CONT void SetSpacing(vtkm::Vec<double, 3> _spacing) { this->Spacing = _spacing; }
VTKM_CONT void SetOrder(vtkm::Int32 _order) { this->Order = _order; }
private:
vtkm::Vec<vtkm::Int32, 3> Radius = { 1, 1, 1 };
double Radius = 1;
vtkm::Vec<double, 3> Spacing = { 1, 1, 1 };
vtkm::Int32 Order = 0;
};
}
......
......@@ -39,7 +39,7 @@ inline VTKM_CONT vtkm::cont::DataSet ComputeMoments::DoExecute(
vtkm::cont::DataSet output;
output.CopyStructure(input);
auto worklet = vtkm::worklet::moments::ComputeMoments(this->Radius);
auto worklet = vtkm::worklet::moments::ComputeMoments(this->Spacing, this->Radius);
worklet.Run(input.GetCellSet(), field, this->Order, output);
......
......@@ -32,13 +32,15 @@ namespace moments
struct ComputeMoments2D : public vtkm::worklet::WorkletPointNeighborhood
{
public:
ComputeMoments2D(const vtkm::Vec<vtkm::Int32, 3>& _radius, int _p, int _q)
: Radius(_radius)
ComputeMoments2D(const vtkm::Vec<double, 3>& _spacing, const double _radius, int _p, int _q)
: Spacing(_spacing)
, Radius(_radius)
, p(_p)
, q(_q)
{
assert(_radius[0] >= 1);
assert(_radius[1] >= 1);
assert(_spacing[0] > 1e-10);
assert(_spacing[1] > 1e-10);
assert(_spacing[2] > 1e-10);
assert(_p >= 0);
assert(_q >= 0);
......@@ -55,28 +57,31 @@ public:
{
// TODO: type safety and numerical precision
auto sum = vtkm::TypeTraits<T>::ZeroInitialization();
vtkm::Vec<vtkm::Float64, 2> recp{ 1.0 / Radius[0], 1.0 / Radius[1] };
// vtkm::Vec<vtkm::Float64, 2> recp{ 1.0 / RadiusReal[0], 1.0 / RadiusReal[1] };
vtkm::Vec<vtkm::Int32, 3> RadiusDiscrete = { this->Radius / (this->Spacing[0] - 1e-10),
this->Radius / (this->Spacing[1] - 1e-10),
this->Radius / (this->Spacing[2] - 1e-10) };
// Clamp the radius to the dataset bounds (discard out-of-bounds points).
const auto minRadius = boundary.ClampNeighborIndex(-this->Radius);
const auto maxRadius = boundary.ClampNeighborIndex(this->Radius);
const auto minRadius = boundary.ClampNeighborIndex(-RadiusDiscrete);
const auto maxRadius = boundary.ClampNeighborIndex(RadiusDiscrete);
for (vtkm::IdComponent j = minRadius[1]; j <= maxRadius[1]; ++j)
{
if (j > -this->Radius[1] && boundary.IJK[1] + j == 0)
if (j > -RadiusDiscrete[1] && boundary.IJK[1] + j == 0)
{ // Don't double count samples that exist on other nodes:
continue;
}
for (vtkm::IdComponent i = minRadius[0]; i <= maxRadius[0]; ++i)
{
if (i > -this->Radius[0] && boundary.IJK[0] + i == 0)
if (i > -RadiusDiscrete[0] && boundary.IJK[0] + i == 0)
{ // Don't double count samples that exist on other nodes:
continue;
}
const vtkm::Float64 r0 = i * recp[0];
const vtkm::Float64 r1 = j * recp[1];
const vtkm::Float64 r0 = i * 1. / RadiusDiscrete[0];
const vtkm::Float64 r1 = j * 1. / RadiusDiscrete[1];
if (r0 * r0 + r1 * r1 <= 1)
{
......@@ -85,11 +90,13 @@ public:
}
}
moment = T(sum * recp[0] * recp[1]);
moment = T(sum * Spacing[0] * Spacing[1]);
}
private:
const vtkm::Vec<vtkm::Int32, 3> Radius;
vtkm::Vec<vtkm::Int32, 3> RadiusDiscrete;
const double Radius;
const vtkm::Vec<double, 3>& Spacing;
const int p;
const int q;
};
......@@ -97,15 +104,20 @@ private:
struct ComputeMoments3D : public vtkm::worklet::WorkletPointNeighborhood
{
public:
ComputeMoments3D(const vtkm::Vec<vtkm::Int32, 3>& _radius, int _p, int _q, int _r)
: Radius(_radius)
ComputeMoments3D(const vtkm::Vec<double, 3>& _spacing,
const double _radius,
int _p,
int _q,
int _r)
: Spacing(_spacing)
, Radius(_radius)
, p(_p)
, q(_q)
, r(_r)
{
assert(_radius[0] >= 1);
assert(_radius[1] >= 1);
assert(_radius[2] >= 1);
assert(_spacing[0] > 1e-10);
assert(_spacing[1] > 1e-10);
assert(_spacing[2] > 1e-10);
assert(_p >= 0);
assert(_q >= 0);
......@@ -123,38 +135,41 @@ public:
{
// TODO: type safety and numerical precision
auto sum = vtkm::TypeTraits<T>::ZeroInitialization();
const vtkm::Vec<vtkm::Float64, 3> recp{ 1.0 / this->Radius[0],
1.0 / this->Radius[1],
1.0 / this->Radius[2] };
// const vtkm::Vec<vtkm::Float64, 3> recp{ 1.0 / this->RadiusReal[0],
// 1.0 / this->RadiusReal[1],
// 1.0 / this->RadiusReal[2] };
vtkm::Vec<vtkm::Int32, 3> RadiusDiscrete = { this->Radius / (this->Spacing[0] - 1e-10),
this->Radius / (this->Spacing[1] - 1e-10),
this->Radius / (this->Spacing[2] - 1e-10) };
// Clamp the radius to the dataset bounds (discard out-of-bounds points).
const auto minRadius = boundary.ClampNeighborIndex(-this->Radius);
const auto maxRadius = boundary.ClampNeighborIndex(this->Radius);
const auto minRadius = boundary.ClampNeighborIndex(-RadiusDiscrete);
const auto maxRadius = boundary.ClampNeighborIndex(RadiusDiscrete);
for (vtkm::IdComponent k = minRadius[2]; k <= maxRadius[2]; ++k)
{
if (k > -this->Radius[2] && boundary.IJK[2] + k == 0)
if (k > -RadiusDiscrete[2] && boundary.IJK[2] + k == 0)
{ // Don't double count samples that exist on other nodes:
continue;
}
for (vtkm::IdComponent j = minRadius[1]; j <= maxRadius[1]; ++j)
{
if (j > -this->Radius[1] && boundary.IJK[1] + j == 0)
if (j > -RadiusDiscrete[1] && boundary.IJK[1] + j == 0)
{ // Don't double count samples that exist on other nodes:
continue;
}
for (vtkm::IdComponent i = minRadius[0]; i <= maxRadius[0]; ++i)
{
if (i > -this->Radius[0] && boundary.IJK[0] + i == 0)
if (i > -RadiusDiscrete[0] && boundary.IJK[0] + i == 0)
{ // Don't double count samples that exist on other nodes:
continue;
}
const vtkm::Float64 r0 = i * recp[0];
const vtkm::Float64 r1 = j * recp[1];
const vtkm::Float64 r2 = k * recp[2];
const vtkm::Float64 r0 = i * 1. / RadiusDiscrete[0];
const vtkm::Float64 r1 = j * 1. / RadiusDiscrete[1];
const vtkm::Float64 r2 = k * 1. / RadiusDiscrete[2];
if (r0 * r0 + r1 * r1 + r2 * r2 <= 1)
{
......@@ -164,11 +179,12 @@ public:
}
}
moment = T(sum * recp[0] * recp[1] * recp[2]);
moment = T(sum * Spacing[0] * Spacing[1] * Spacing[2]);
}
private:
const vtkm::Vec<vtkm::Int32, 3>& Radius;
const double Radius;
const vtkm::Vec<double, 3>& Spacing;
const int p;
const int q;
const int r;
......@@ -177,8 +193,9 @@ private:
class ComputeMoments
{
public:
ComputeMoments(vtkm::Vec<vtkm::Int32, 3> _radius)
: Radius(_radius)
ComputeMoments(vtkm::Vec<double, 3>& _spacing, const double _radius)
: Spacing(_spacing)
, Radius(_radius)
{
}
......@@ -188,7 +205,8 @@ public:
template <typename T, typename S>
void operator()(const vtkm::cont::CellSetStructured<2>& input,
const vtkm::cont::ArrayHandle<T, S>& pixels,
vtkm::Vec<vtkm::Int32, 3> Radius,
vtkm::Vec<double, 3> Spacing,
double Radius,
int maxOrder,
vtkm::cont::DataSet& output) const
{
......@@ -203,7 +221,7 @@ public:
vtkm::cont::ArrayHandle<T> moments;
DispatcherType dispatcher(WorkletType{ Radius, p, q });
DispatcherType dispatcher(WorkletType{ Spacing, Radius, p, q });
dispatcher.Invoke(input, pixels, moments);
std::string fieldName = std::string("index") + std::string(p, '0') + std::string(q, '1');
......@@ -218,7 +236,8 @@ public:
template <typename T, typename S>
void operator()(const vtkm::cont::CellSetStructured<3>& input,
const vtkm::cont::ArrayHandle<T, S>& pixels,
vtkm::Vec<vtkm::Int32, 3> Radius,
vtkm::Vec<double, 3> Spacing,
double Radius,
int maxOrder,
vtkm::cont::DataSet& output) const
{
......@@ -236,7 +255,7 @@ public:
vtkm::cont::ArrayHandle<T> moments;
DispatcherType dispatcher(WorkletType{ Radius, p, q, r });
DispatcherType dispatcher(WorkletType{ Spacing, Radius, p, q, r });
dispatcher.Invoke(input, pixels, moments);
std::string fieldName = std::string("index") + std::string(p, '0') +
......@@ -258,11 +277,12 @@ public:
vtkm::cont::DataSet& output) const
{
input.ResetCellSetList(vtkm::cont::CellSetListTagStructured())
.CastAndCall(ResolveDynamicCellSet(), pixels, this->Radius, maxOrder, output);
.CastAndCall(ResolveDynamicCellSet(), pixels, this->Spacing, this->Radius, maxOrder, output);
}
private:
const vtkm::Vec<vtkm::Int32, 3> Radius = { 1, 1, 1 };
const double Radius = 1;
const vtkm::Vec<double, 3> Spacing = { 1, 1, 1 };
};
}
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment