Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-10-08 03:18:58 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Add new ContinuousScatterPlot filter

Browse Source 
This new filtered designed for bi-variate analysis builds the continuous scatterplot of a 3D tetrahedralized mesh for two given scalar point fields. The continuous scatterplot is an extension of the discrete scatterplot for continuous bi-variate analysis.
This commit is contained in:
Louis Gombert 2023-03-21 09:43:35 +01:00
parent 90706d53f8
commit 4ba15a993d
9 changed files with 1317 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
docs/changelog/continuous-scatterplot.md Normal file
View File

@ -0,0 +1,7 @@
# Continuous Scatterplot filter
This new filter, designed for bi-variate analysis, computes the continuous scatter-plot of a 3D mesh for two given scalar point fields.
The continuous scatter-plot is an extension of the discrete scatter-plot for continuous bi-variate analysis. The output dataset consists of triangle-shaped cells, whose coordinates on the 2D plane represent respectively the values of both scalar fields. The point centered scalar field generated on the triangular mesh quantifies the density of values in the data domain.
This VTK-m implementation is based on the algorithm presented in the paper "Continuous Scatterplots" by S. Bachthaler and D. Weiskopf.

2
vtkm/filter/density_estimate/CMakeLists.txt
View File

@ -8,6 +8,7 @@
## PURPOSE. See the above copyright notice for more information.
##============================================================================
set(density_estimate_headers
ContinuousScatterPlot.h
Entropy.h
Histogram.h
NDEntropy.h
@ -19,6 +20,7 @@ set(density_estimate_headers
)
set(density_estimate_sources_device
ContinuousScatterPlot.cxx
Entropy.cxx
Histogram.cxx
NDEntropy.cxx

82
vtkm/filter/density_estimate/ContinuousScatterPlot.cxx Normal file
View File

@ -0,0 +1,82 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/CoordinateSystem.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/density_estimate/ContinuousScatterPlot.h>
#include <vtkm/filter/density_estimate/worklet/ContinuousScatterPlot.h>
#include <vtkm/filter/geometry_refinement/Tetrahedralize.h>
namespace vtkm
{
namespace filter
{
namespace density_estimate
{
VTKM_CONT vtkm::cont::DataSet ContinuousScatterPlot::DoExecute(const vtkm::cont::DataSet& input)
{
// This algorithm only operate on tetra cells, we need to apply the tetrahedralize filter first.
auto tetrahedralizeFilter = vtkm::filter::geometry_refinement::Tetrahedralize();
auto tetraInput = tetrahedralizeFilter.Execute(input);
vtkm::cont::CellSetSingleType<> tetraCellSet;
tetraInput.GetCellSet().AsCellSet(tetraCellSet);
vtkm::cont::Field scalarField1 = input.GetField(GetActiveFieldName(0));
vtkm::cont::Field scalarField2 = input.GetField(GetActiveFieldName(1));
if (!(scalarField1.IsPointField() && scalarField2.IsPointField()))
{
throw vtkm::cont::ErrorFilterExecution("Point fields expected.");
}
const auto& coords = input.GetCoordinateSystem().GetDataAsMultiplexer();
vtkm::cont::CoordinateSystem activeCoordSystem = input.GetCoordinateSystem();
vtkm::cont::DataSet scatterplotDataSet;
vtkm::worklet::ContinuousScatterPlot worklet;
auto resolveFieldType = [&](const auto& resolvedScalar) {
using FieldType = typename std::decay_t<decltype(resolvedScalar)>::ValueType;
vtkm::cont::CellSetSingleType<> scatterplotCellSet;
vtkm::cont::ArrayHandle<FieldType> density;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> newCoords;
// Both fields need to resolve to the same type in order to perform calculations
vtkm::cont::ArrayHandle<FieldType> resolvedScalar2;
vtkm::cont::ArrayCopyShallowIfPossible(scalarField2.GetData(), resolvedScalar2);
worklet.Run(tetraCellSet,
coords,
newCoords,
density,
resolvedScalar,
resolvedScalar2,
scatterplotCellSet);
// Populate the new dataset representing the continuous scatterplot
// Using the density field and coordinates calculated by the worklet
activeCoordSystem = vtkm::cont::CoordinateSystem(activeCoordSystem.GetName(), newCoords);
scatterplotDataSet.AddCoordinateSystem(activeCoordSystem);
scatterplotDataSet.SetCellSet(scatterplotCellSet);
scatterplotDataSet.AddPointField(this->GetOutputFieldName(), density);
};
this->CastAndCallScalarField(scalarField1.GetData(), resolveFieldType);
return scatterplotDataSet;
}
} // namespace density_estimate
} // namespace filter
} // namespace vtkm

60
vtkm/filter/density_estimate/ContinuousScatterPlot.h Normal file
View File

@ -0,0 +1,60 @@
//============================================================================
// 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.
//============================================================================
#ifndef vtk_m_filter_density_estimate_ContinuousScatterPlot_h
#define vtk_m_filter_density_estimate_ContinuousScatterPlot_h
#include <vtkm/filter/FilterField.h>
#include <vtkm/filter/density_estimate/vtkm_filter_density_estimate_export.h>
namespace vtkm
{
namespace filter
{
namespace density_estimate
{
/// \brief Constructs the continuous scatterplot for two given scalar point fields of a mesh.
///
/// The continuous scatterplot is an extension of the discrete scatterplot for continuous bi-variate analysis.
/// This filter outputs an ExplicitDataSet of triangle-shaped cells, whose coordinates on the 2D plane represent respectively
/// the values of both scalar fields. Triangles' points are associated with a scalar field, representing the
/// density of values in the data domain. The filter tetrahedralizes the input dataset before operating.
///
/// If both fields provided don't have the same floating point precision, the output will
/// have the precision of the first one of the pair.
///
/// This implementation is based on the algorithm presented in the publication :
///
/// S. Bachthaler and D. Weiskopf, "Continuous Scatterplots"
/// in IEEE Transactions on Visualization and Computer Graphics,
/// vol. 14, no. 6, pp. 1428-1435, Nov.-Dec. 2008
/// doi: 10.1109/TVCG.2008.119.
class VTKM_FILTER_DENSITY_ESTIMATE_EXPORT ContinuousScatterPlot : public vtkm::filter::FilterField
{
public:
VTKM_CONT ContinuousScatterPlot() { this->SetOutputFieldName("density"); }
/// Select both point fields to use when running the filter.
/// Replaces setting each one individually using `SetActiveField` on indices 0 and 1.
VTKM_CONT
void SetActiveFieldsPair(const std::string& fieldName1, const std::string& fieldName2)
{
SetActiveField(0, fieldName1, vtkm::cont::Field::Association::Points);
SetActiveField(1, fieldName2, vtkm::cont::Field::Association::Points);
};
private:
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override;
};
} // namespace density_estimate
} // namespace filter
} // namespace vtm
#endif //vtk_m_filter_density_estimate_ContinuousScatterPlot_h

1
vtkm/filter/density_estimate/testing/CMakeLists.txt
View File

@ -9,6 +9,7 @@
##============================================================================
set(unit_tests
UnitTestContinuousScatterPlot.cxx
UnitTestEntropyFilter.cxx
UnitTestHistogramFilter.cxx
UnitTestNDEntropyFilter.cxx

690
vtkm/filter/density_estimate/testing/UnitTestContinuousScatterPlot.cxx Normal file
View File

@ -0,0 +1,690 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderExplicit.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/density_estimate/ContinuousScatterPlot.h>
namespace
{
std::vector<vtkm::Vec3f> tetraCoords{ vtkm::Vec3f(0.0f, 0.0f, 0.0f),
vtkm::Vec3f(2.0f, 0.0f, 0.0f),
vtkm::Vec3f(2.0f, 2.0f, 0.0f),
vtkm::Vec3f(1.0f, 0.0f, 2.0f) };
std::vector<vtkm::UInt8> tetraShape{ vtkm::CELL_SHAPE_TETRA };
std::vector<vtkm::IdComponent> tetraIndex{ 4 };
std::vector<vtkm::Id> tetraConnectivity{ 0, 1, 2, 3 };
std::vector<vtkm::Vec3f> multiCoords{
vtkm::Vec3f(0.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 2.0f, 0.0f),
vtkm::Vec3f(1.0f, 0.0f, 2.0f), vtkm::Vec3f(0.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 0.0f, 0.0f),
vtkm::Vec3f(2.0f, 2.0f, 0.0f), vtkm::Vec3f(1.0f, 0.0f, 2.0f), vtkm::Vec3f(0.0f, 0.0f, 0.0f),
vtkm::Vec3f(2.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 2.0f, 0.0f), vtkm::Vec3f(1.0f, 0.0f, 2.0f)
};
std::vector<vtkm::UInt8> multiShapes{ vtkm::CELL_SHAPE_TETRA,
vtkm::CELL_SHAPE_TETRA,
vtkm::CELL_SHAPE_TETRA };
std::vector<vtkm::IdComponent> multiIndices{ 4, 4, 4 };
std::vector<vtkm::Id> multiConnectivity{ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
template <typename FieldType1, typename FieldType2>
vtkm::cont::DataSet makeDataSet(const std::vector<vtkm::Vec3f>& ds_coords,
const std::vector<vtkm::UInt8>& ds_shapes,
const std::vector<vtkm::IdComponent>& ds_indices,
const std::vector<vtkm::Id>& ds_connectivity,
const FieldType1& scalar1,
const FieldType2& scalar2)
{
vtkm::cont::DataSetBuilderExplicit dsb;
vtkm::cont::DataSet ds = dsb.Create(ds_coords, ds_shapes, ds_indices, ds_connectivity);
ds.AddPointField("scalar1", scalar1, 4);
ds.AddPointField("scalar2", scalar2, 4);
return ds;
}
vtkm::cont::DataSet executeFilter(const vtkm::cont::DataSet ds)
{
auto continuousSCP = vtkm::filter::density_estimate::ContinuousScatterPlot();
continuousSCP.SetActiveFieldsPair("scalar1", "scalar2");
return continuousSCP.Execute(ds);
}
template <typename PositionsPortalType, typename FieldType1, typename FieldType2>
void testCoordinates(const PositionsPortalType& positionsP,
const FieldType1& scalar1,
const FieldType2& scalar2,
const vtkm::IdComponent numberOfPoints)
{
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
{
VTKM_TEST_ASSERT(test_equal(positionsP.Get(i)[0], scalar1[i]), "Wrong point coordinates");
VTKM_TEST_ASSERT(test_equal(positionsP.Get(i)[1], scalar2[i]), "Wrong point coordinates");
VTKM_TEST_ASSERT(test_equal(positionsP.Get(i)[2], 0.0f),
"Z coordinate value in the scatterplot should always be null");
}
}
template <typename DensityArrayType, typename FieldType>
void testDensity(const DensityArrayType& density,
const vtkm::IdComponent& centerId,
const FieldType& centerDensity)
{
for (vtkm::IdComponent i = 0; i < density.GetNumberOfValues(); i++)
{
if (i == centerId)
{
VTKM_TEST_ASSERT(test_equal(density.Get(i), centerDensity),
"Wrong density in the middle point of the cell");
}
else
{
VTKM_TEST_ASSERT(test_equal(density.Get(i), 0.0f),
"Density on the edge of the tetrahedron should be null");
}
}
}
template <typename CellSetType>
void testShapes(const CellSetType& cellSet)
{
for (vtkm::IdComponent i = 0; i < cellSet.GetNumberOfCells(); i++)
{
VTKM_TEST_ASSERT(cellSet.GetCellShape(i) == vtkm::CELL_SHAPE_TRIANGLE);
}
}
template <typename CellSetType>
void testConnectivity(const CellSetType& cellSet,
const vtkm::cont::ArrayHandle<vtkm::IdComponent>& expectedConnectivityArray)
{
VTKM_TEST_ASSERT(
test_equal_ArrayHandles(
cellSet.GetConnectivityArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint()),
expectedConnectivityArray),
"Wrong connectivity");
}
void TestSingleTetraProjectionQuadConvex()
{
constexpr vtkm::FloatDefault scalar1[4] = {
0.0f,
1.0f,
0.0f,
-2.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
-1.0f,
0.0f,
2.0f,
0.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 4),
"Wrong number of projected triangles in the continuous scatter plot");
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfPoints(), 5),
"Wrong number of projected points in the continuous scatter plot");
// Test point positions
auto positions = scatterPlot.GetCoordinateSystem().GetDataAsMultiplexer();
auto positionsP = positions.ReadPortal();
testCoordinates(positionsP, scalar1, scalar2, 4);
// Test intersection point coordinates
VTKM_TEST_ASSERT(test_equal(positionsP.Get(4), vtkm::Vec3f{ 0.0f, 0.0f, 0.0f }),
"Wrong intersection point coordinates");
// Test for triangle shapes
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
testShapes(cellSet);
// Test connectivity
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 1, 4, 1, 2, 4, 2, 3, 4, 3, 0, 4 });
testConnectivity(cellSet, expectedConnectivityArray);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
testDensity(density, 4, 0.888889f);
}
void TestSingleTetraProjectionQuadSelfIntersect()
{
constexpr vtkm::FloatDefault scalar1[4] = {
0.0f,
0.0f,
1.0f,
-2.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
-1.0f,
2.0f,
0.0f,
0.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 4),
"Wrong number of projected triangles in the continuous scatter plot");
// Test point positions
auto positions = scatterPlot.GetCoordinateSystem().GetDataAsMultiplexer();
auto positionsP = positions.ReadPortal();
testCoordinates(positionsP, scalar1, scalar2, 4);
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 2, 4, 2, 1, 4, 1, 3, 4, 3, 0, 4 });
testConnectivity(cellSet, expectedConnectivityArray);
}
void TestSingleTetraProjectionQuadInverseOrder()
{
constexpr vtkm::FloatDefault scalar1[4] = {
-2.0f,
0.0f,
1.0f,
0.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
0.0f,
2.0f,
0.0f,
-1.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 4),
"Wrong number of projected triangles in the continuous scatter plot");
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray = vtkm::cont::make_ArrayHandle<vtkm::IdComponent>(
{ 0,
1,
4,
1,
2,
4,
2,
3,
4,
3,
0,
4 }); // Inverting the order of points should not change connectivity
testConnectivity(cellSet, expectedConnectivityArray);
}
void TestSingleTetraProjectionQuadSelfIntersectSecond()
{
constexpr vtkm::FloatDefault scalar1[4] = {
0.0f,
1.0f,
-2.0f,
0.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
-1.0f,
0.0f,
0.0f,
2.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 4),
"Wrong number of projected triangles in the continuous scatter plot");
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 2, 4, 2, 3, 4, 3, 1, 4, 1, 0, 4 });
testConnectivity(cellSet, expectedConnectivityArray);
}
void TestSingleTetra_projection_triangle_point0Inside()
{
constexpr vtkm::FloatDefault scalar1[4] = {
3.0f,
3.0f,
4.0f,
1.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
1.0f,
0.0f,
2.0f,
2.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 3),
"Wrong number of projected triangles in the continuous scatter plot");
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfPoints(), 4),
"Wrong number of projected points in the continuous scatter plot");
// Test point positions
auto positions = scatterPlot.GetCoordinateSystem().GetDataAsMultiplexer();
auto positionsP = positions.ReadPortal();
testCoordinates(positionsP, scalar1, scalar2, 3);
// Test for triangle shapes
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
testShapes(cellSet);
// Test connectivity
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 1, 2, 0, 2, 3, 0, 3, 1, 0 });
testConnectivity(cellSet, expectedConnectivityArray);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
testDensity(density, 0, 1.33333f);
}
void TestSingleTetra_projection_triangle_point1Inside()
{
constexpr vtkm::FloatDefault scalar1[4] = {
3.0f,
3.0f,
4.0f,
1.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
0.0f,
1.0f,
2.0f,
2.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 2, 1, 2, 3, 1, 3, 0, 1 });
testConnectivity(cellSet, expectedConnectivityArray);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
testDensity(density, 1, 1.33333f);
}
void TestSingleTetra_projection_triangle_point2Inside()
{
constexpr vtkm::FloatDefault scalar1[4] = {
3.0f,
4.0f,
3.0f,
1.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
0.0f,
2.0f,
1.0f,
2.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 1, 2, 1, 3, 2, 3, 0, 2 });
testConnectivity(cellSet, expectedConnectivityArray);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
testDensity(density, 2, 1.33333f);
}
void TestSingleTetra_projection_triangle_point3Inside()
{
constexpr vtkm::FloatDefault scalar1[4] = {
3.0f,
4.0f,
1.0f,
3.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
0.0f,
2.0f,
2.0f,
1.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 1, 3, 1, 2, 3, 2, 0, 3 });
testConnectivity(cellSet, expectedConnectivityArray);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
testDensity(density, 3, 1.33333f);
}
void TestNullSpatialVolume()
{
constexpr vtkm::FloatDefault scalar1[4] = {
3.0f,
3.0f,
4.0f,
1.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
1.0f,
0.0f,
2.0f,
2.0f,
};
std::vector<vtkm::Vec3f> nullCoords{ vtkm::Vec3f(1.0f, 1.0f, 1.0f),
vtkm::Vec3f(1.0f, 1.0f, 1.0f),
vtkm::Vec3f(1.0f, 1.0f, 1.0f),
vtkm::Vec3f(1.0f, 1.0f, 1.0f) };
vtkm::cont::DataSet ds =
makeDataSet(nullCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
VTKM_TEST_ASSERT(vtkm::IsInf(density.Get(0)), "Density should be infinite for a null volume");
}
void TestNullDataSurface()
{
constexpr vtkm::FloatDefault scalar1[4] = {
0.0f,
1.0f,
3.0f,
2.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
0.0f,
1.0f,
3.0f,
2.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
VTKM_TEST_ASSERT(vtkm::IsInf(density.Get(4)), "Density should be infinite for a null volume");
}
void TestMultipleTetra()
{
constexpr vtkm::FloatDefault multiscalar1[12] = { 3.0f, 3.0f, 4.0f, 1.0f, 0.0f, 1.0f,
0.0f, -2.0f, 3.0f, 3.0f, 4.0f, 1.0f };
constexpr vtkm::FloatDefault multiscalar2[12] = { 1.0f, 0.0f, 2.0f, 2.0f, -1.0f, 0.0f,
2.0f, 0.0f, 1.0f, 0.0f, 2.0f, 2.0f };
vtkm::cont::DataSetBuilderExplicit dsb;
vtkm::cont::DataSet ds = dsb.Create(multiCoords, multiShapes, multiIndices, multiConnectivity);
ds.AddPointField("scalar1", multiscalar1, 12);
ds.AddPointField("scalar2", multiscalar2, 12);
// Filtering
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 10),
"Wrong number of projected triangles in the continuous scatter plot");
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfPoints(), 13),
"Wrong number of projected points in the continuous scatter plot");
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
testShapes(cellSet);
}
void TestNonTetra()
{
std::vector<vtkm::Vec3f> wedgeCoords{
vtkm::Vec3f(0.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 4.0f, 0.0f),
vtkm::Vec3f(0.0f, 4.0f, 0.0f), vtkm::Vec3f(1.0f, 0.0f, 3.0f), vtkm::Vec3f(1.0f, 4.0f, 3.0f),
};
constexpr vtkm::FloatDefault scalar1[6] = { 0.0f, 3.0f, 3.0f, 2.0f, 2.0f, 1.0f };
constexpr vtkm::FloatDefault scalar2[6] = { 0.0f, 1.0f, 3.0f, 2.0f, 0.0f, 1.0f };
std::vector<vtkm::UInt8> w_shape{ vtkm::CELL_SHAPE_WEDGE };
std::vector<vtkm::IdComponent> w_indices{ 6 };
std::vector<vtkm::Id> w_connectivity{ 0, 1, 2, 3, 4, 5 };
vtkm::cont::DataSetBuilderExplicit dsb;
vtkm::cont::DataSet ds = dsb.Create(wedgeCoords, w_shape, w_indices, w_connectivity);
ds.AddPointField("scalar1", scalar1, 6);
ds.AddPointField("scalar2", scalar2, 6);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 12),
"Wrong number of projected triangles in the continuous scatter plot");
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfPoints(), 15),
"Wrong number of projected points in the continuous scatter plot");
}
void TestNonPointFields()
{
constexpr vtkm::FloatDefault cellField1[1] = { 1.0f };
constexpr vtkm::FloatDefault cellField2[1] = {
0.0f,
};
vtkm::cont::DataSetBuilderExplicit dsb;
vtkm::cont::DataSet ds = dsb.Create(tetraCoords, tetraShape, tetraIndex, tetraConnectivity);
ds.AddCellField("scalar1", cellField1, 1);
ds.AddCellField("scalar2", cellField2, 1);
auto continuousSCP = vtkm::filter::density_estimate::ContinuousScatterPlot();
continuousSCP.SetActiveField(0, "scalar1", vtkm::cont::Field::Association::Cells);
continuousSCP.SetActiveField(1, "scalar2", vtkm::cont::Field::Association::Cells);
try
{
auto scatterPlot = continuousSCP.Execute(ds);
VTKM_TEST_FAIL(
"Filter execution was not aborted after providing active fields not associated with points");
}
catch (const vtkm::cont::ErrorFilterExecution&)
{
std::cout << "Execution successfully aborted" << std::endl;
}
}
void TestDataTypes()
{
constexpr vtkm::Float32 scalar1_f32[4] = {
-2.0f,
0.0f,
1.0f,
0.0f,
};
constexpr vtkm::Float32 scalar2_f32[4] = {
0.0f,
2.0f,
0.0f,
-1.0f,
};
constexpr vtkm::Float64 scalar1_f64[4] = {
-2.0,
0.0,
1.0,
0.0,
};
constexpr vtkm::Float64 scalar2_f64[4] = {
0.0,
2.0,
0.0,
-1.0,
};
// The filter should run whatever float precision we use for both fields
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1_f32, scalar2_f32);
auto scatterPlot = executeFilter(ds);
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>()
.ReadPortal();
testDensity(density, 4, 0.888889f);
vtkm::cont::DataSet ds2 =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1_f32, scalar2_f64);
auto scatterPlot2 = executeFilter(ds);
auto density2 = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>()
.ReadPortal();
testDensity(density2, 4, 0.888889f);
vtkm::cont::DataSet ds3 =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1_f64, scalar2_f32);
auto scatterPlot3 = executeFilter(ds);
auto density3 = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>()
.ReadPortal();
testDensity(density3, 4, 0.888889f);
vtkm::cont::DataSet ds4 =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1_f64, scalar2_f64);
auto scatterPlot4 = executeFilter(ds);
auto density4 = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>()
.ReadPortal();
testDensity(density4, 4, 0.888889f);
}
void TestContinuousScatterPlot()
{
// Projection forms 4 triangles in the data domain
TestSingleTetraProjectionQuadConvex();
TestSingleTetraProjectionQuadSelfIntersect();
TestSingleTetraProjectionQuadInverseOrder();
TestSingleTetraProjectionQuadSelfIntersectSecond();
// 3 triangles in the data domain
TestSingleTetra_projection_triangle_point0Inside();
TestSingleTetra_projection_triangle_point1Inside();
TestSingleTetra_projection_triangle_point2Inside();
TestSingleTetra_projection_triangle_point3Inside();
// // Edge cases
TestNullSpatialVolume();
TestNullDataSurface();
TestMultipleTetra();
TestNonTetra();
TestNonPointFields();
TestDataTypes();
}
}
int UnitTestContinuousScatterPlot(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestContinuousScatterPlot, argc, argv);
}

1
vtkm/filter/density_estimate/vtkm.module
View File

@ -6,6 +6,7 @@ DEPENDS
vtkm_filter_core
PRIVATE_DEPENDS
vtkm_worklet
vtkm_filter_geometry_refinement
OPTIONAL_DEPENDS
MPI::MPI_CXX
TEST_DEPENDS

1
vtkm/filter/density_estimate/worklet/CMakeLists.txt
View File

@ -9,6 +9,7 @@
##============================================================================
set(headers
ContinuousScatterPlot.h
FieldEntropy.h
FieldHistogram.h
NDimsEntropy.h

473
vtkm/filter/density_estimate/worklet/ContinuousScatterPlot.h Normal file
View File

@ -0,0 +1,473 @@
//============================================================================
// 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.
//============================================================================
#ifndef vtk_m_worklet_ContinuousScatterPlot_h
#define vtk_m_worklet_ContinuousScatterPlot_h
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayHandleGroupVec.h>
#include <vtkm/cont/ArrayHandleGroupVecVariable.h>
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/worklet/ScatterCounting.h>
#include <vtkm/worklet/WorkletMapTopology.h>
namespace vtkm
{
namespace worklet
{
class ContinuousScatterPlot
{
public:
template <typename FieldType>
struct ClassifyTetra : vtkm::worklet::WorkletVisitCellsWithPoints
{
using ControlSignature = void(CellSetIn,
FieldInPoint scalars,
FieldOutCell pointsOrder,
FieldOutCell numberOfPoints,
FieldOutCell numberOfTris);
using ExecutionSignature = void(_2 scalar, _3 pointsOrder, _4 numberOfPoints, _5 numberOfTris);
using InputDomain = _1;
VTKM_EXEC FieldType ZCrossProduct(const vtkm::Pair<FieldType, FieldType>& point1,
const vtkm::Pair<FieldType, FieldType>& point2,
const vtkm::Pair<FieldType, FieldType>& point3) const
{
return vtkm::DifferenceOfProducts(point2.first - point1.first,
point3.second - point1.second,
point2.second - point1.second,
point3.first - point1.first);
}
VTKM_EXEC bool DifferentSign(const FieldType& value1, const FieldType& value2) const
{
return (vtkm::IsNegative(value1) != vtkm::IsNegative(value2));
}
VTKM_EXEC bool AllSameSign(const FieldType& value1,
const FieldType& value2,
const FieldType& value3) const
{
return (!(DifferentSign(value1, value2) || DifferentSign(value2, value3)));
}
template <typename ScalarField, typename PointsOutOrder>
VTKM_EXEC void operator()(const ScalarField& scalar,
PointsOutOrder& pointsOrder,
vtkm::IdComponent& numberOfPoints,
vtkm::IdComponent& numberOfTris) const
{
// To classify our tetras following their projection in the 2D scalar domain,
// We consider them as quads, with their coordinates being their respective scalar values.
// To identify the projection, we want to know if the polygon formed by the 4 points of the quad is convex.
// For this, we compute the Z component of the cross product of the vectors of the polygon's edges.
vtkm::Vec<FieldType, 4> scalarCrossProduct{ ZCrossProduct(scalar[0], scalar[1], scalar[2]),
ZCrossProduct(scalar[1], scalar[2], scalar[3]),
ZCrossProduct(scalar[2], scalar[3], scalar[0]),
ZCrossProduct(scalar[3], scalar[0], scalar[1]) };
// If every cross product of consecutive edges of the quad is the same sign, it means that it is convex
// In the case 2 of them are negative and 2 positive, the quad is self-intersecting.
// If one or 3 of them are negative, we have found a non-convex quad, projecting 3 triangles.
if ((DifferentSign(scalarCrossProduct[0], scalarCrossProduct[1]) !=
DifferentSign(scalarCrossProduct[2], scalarCrossProduct[3])))
{
numberOfPoints = 4;
numberOfTris = 3;
// Here, one of the 4 points is in the triangle formed by the 3 other.
// We can guess which one it is by analyzing which cross product has a different sign than the other 3.
// Assign this specific point's id to element 3 of our order array.
if (AllSameSign(scalarCrossProduct[1], scalarCrossProduct[2], scalarCrossProduct[3]))
{
pointsOrder[0] = 0;
pointsOrder[1] = 2;
pointsOrder[2] = 3;
pointsOrder[3] = 1;
}
else if (AllSameSign(scalarCrossProduct[0], scalarCrossProduct[2], scalarCrossProduct[3]))
{
pointsOrder[0] = 0;
pointsOrder[1] = 1;
pointsOrder[2] = 3;
pointsOrder[3] = 2;
}
else if (AllSameSign(scalarCrossProduct[0], scalarCrossProduct[1], scalarCrossProduct[3]))
{
pointsOrder[0] = 0;
pointsOrder[1] = 1;
pointsOrder[2] = 2;
pointsOrder[3] = 3;
}
else
{
pointsOrder[0] = 1;
pointsOrder[1] = 2;
pointsOrder[2] = 3;
pointsOrder[3] = 0;
}
}
else
{
// This tetra projection yields 4 triangles,
// And forms a convex quad in the data plane
numberOfPoints = 5;
numberOfTris = 4;
// Find an order of points which forms a non self-intersecting quad,
// Still using the Z cross-product signs
if (DifferentSign(scalarCrossProduct[0], scalarCrossProduct[1]))
{
// Polygon Self-intersects on the second diagonal
pointsOrder[0] = 0;
pointsOrder[1] = 2;
pointsOrder[2] = 3;
pointsOrder[3] = 1;
}
else if (DifferentSign(scalarCrossProduct[1], scalarCrossProduct[2]))
{
// Self-intersect on the first diagonal
pointsOrder[0] = 0;
pointsOrder[1] = 2;
pointsOrder[2] = 1;
pointsOrder[3] = 3;
}
else
{
// Convex
pointsOrder[0] = 0;
pointsOrder[1] = 1;
pointsOrder[2] = 2;
pointsOrder[3] = 3;
}
}
}
};
template <typename FieldType>
struct VolumeMeasure : vtkm::worklet::WorkletVisitCellsWithPoints
{
using ControlSignature = void(CellSetIn,
FieldInPoint scalars,
FieldInPoint coords,
FieldInCell numberOfTris,
FieldInCell ord,
FieldOutCell newCoords,
FieldOutCell density);
using ExecutionSignature =
void(_2 scalars, _3 coords, _4 numberOfTris, _5 ord, _6 newCoords, _7 density);
using InputDomain = _1;
using Vec3t = vtkm::Vec<FieldType, 3>;
template <typename VecType>
VTKM_EXEC Vec3t DiagonalIntersection(const VecType& point1,
const VecType& point2,
const VecType& point3,
const VecType& point4) const
{
FieldType denominator = vtkm::DifferenceOfProducts(
point1[0] - point2[0], point3[1] - point4[1], point1[1] - point2[1], point3[0] - point4[0]);
// In case the points are aligned, return arbitrarily the first point as intersection.
// These points represent the diagonals of a convex polygon,
// so they are either intersecting or lying on the same line
// The surface area of the quad in the data domain will be null in this case anyways.
if (denominator == 0)
{
return point1;
}
// Otherwise, compute the intersection point of the quad's diagonals, defined by 4 points.
// This vector is the point we get when equating the line equations (point1,point2) and (point3,point4)
return Vec3t{ vtkm::DifferenceOfProducts(
point3[0] - point4[0],
vtkm::DifferenceOfProducts(point1[0], point2[1], point1[1], point2[0]),
point1[0] - point2[0],
vtkm::DifferenceOfProducts(point3[0], point4[1], point3[1], point4[0])) /
denominator,
vtkm::DifferenceOfProducts(
point3[1] - point4[1],
vtkm::DifferenceOfProducts(point1[0], point2[1], point1[1], point2[0]),
point1[1] - point2[1],
vtkm::DifferenceOfProducts(point3[0], point4[1], point3[1], point4[0])) /
denominator,
0.0f };
}
VTKM_EXEC FieldType TriangleArea(Vec3t point1, Vec3t point2, Vec3t point3) const
{
return 0.5f * vtkm::Magnitude(vtkm::Cross(point2 - point1, point3 - point1));
}
VTKM_EXEC FieldType Distance(Vec3t point1, Vec3t point2) const
{
return vtkm::Magnitude(point1 - point2);
}
template <typename ScalarField,
typename CoordsType,
typename NewCoordsType,
typename PointsOutOrder,
typename DensityField>
VTKM_EXEC void operator()(const ScalarField& scalar,
const CoordsType& coords,
const vtkm::IdComponent& numberOfTris,
PointsOutOrder& ord,
NewCoordsType& newCoords,
DensityField& density) const
{
// Write points coordinates in the 2D plane based on provided scalar values
for (int i = 0; i < 4; i++)
{
newCoords[i] = Vec3t{ scalar[i].first, scalar[i].second, 0.0f };
}
if (numberOfTris == 4)
{
// If the projection has 4 triangles, the fifth imaginary point is defined as the intersection of the quad's diagonals.
newCoords[4] = DiagonalIntersection(
newCoords[ord[0]], newCoords[ord[2]], newCoords[ord[1]], newCoords[ord[3]]);
}
// Compute densities
// The density on the borders of the data domain is always null.
// For each tetra projection the only density > 0 is associated either to the point
// located inside the triangle formed by the others (for 3 triangles projection),
// Or to the imaginary point calculated above.
density[ord[0]] = 0.0f;
density[ord[1]] = 0.0f;
density[ord[2]] = 0.0f;
density[ord[3]] = 0.0f;
// Pre-compute some vectors to reuse later
vtkm::Vec<Vec3t, 3> spatialVector = { coords[1] - coords[0],
coords[2] - coords[0],
coords[3] - coords[0] };
vtkm::Vec<Vec3t, 3> spatialCrossProducts = { vtkm::Cross(spatialVector[1], spatialVector[0]),
vtkm::Cross(spatialVector[0], spatialVector[2]),
vtkm::Cross(spatialVector[2],
spatialVector[1]) };
vtkm::Vec<Vec3t, 2> scalarArray = { Vec3t{ scalar[1].first - scalar[0].first,
scalar[2].first - scalar[0].first,
scalar[3].first - scalar[0].first },
Vec3t{ Vec3t{ scalar[1].second - scalar[0].second,
scalar[2].second - scalar[0].second,
scalar[3].second - scalar[0].second } } };
// We need to calculate the determinant in the spatial domain, using the triple product formula
FieldType determinant = vtkm::Dot(spatialVector[2], spatialCrossProducts[0]);
// Calculate the spatial gradient for both scalar fields in the tetrahedron
vtkm::Vec<Vec3t, 2> scalar_gradient;
// The determinant is null, therefore the volume is null
if (determinant == 0.0f)
{
scalar_gradient = { Vec3t{ 0.0f, 0.0f, 0.0f }, Vec3t{ 0.0f, 0.0f, 0.0f } };
}
else
{
// This gradient formulation is derived from the calculation of the inverse Jacobian matrix,
// dividing the adjugate matrix of the transformation by the determinant.
// Each column of the matrix is then multiplied by the scalar difference
// with the scalar value for point with index 0 and summed to get the gradient, for each scalar field.
FieldType inv_determinant = 1.0f / determinant;
scalar_gradient = vtkm::Vec<Vec3t, 2>{ inv_determinant *
(spatialCrossProducts[0] * scalarArray[0][2] +
spatialCrossProducts[1] * scalarArray[0][1] +
spatialCrossProducts[2] * scalarArray[0][0]),
inv_determinant *
(spatialCrossProducts[0] * scalarArray[1][2] +
spatialCrossProducts[1] * scalarArray[1][1] +
spatialCrossProducts[2] * scalarArray[1][0]) };
}
// We get the volume measure, defined as the magnitude of the cross product of the gradient
// See formula (10) in the "continuous scatterplots" paper for the demonstration
FieldType volume = vtkm::Magnitude(vtkm::Cross(scalar_gradient[0], scalar_gradient[1]));
if (numberOfTris == 3)
{
// Calculate the area of the triangle on the backface of the projected tetra
FieldType fullArea =
this->TriangleArea(newCoords[ord[0]], newCoords[ord[1]], newCoords[ord[2]]);
if (volume == 0.0f || fullArea == 0.0f)
{
// For a tetrahedra of null volume, the density is infinite
density[ord[3]] = vtkm::Infinity<FieldType>();
return;
}
// The density for the central point is the distance to the backface divided by the volume
// We interpolate the position of point [3] using the scalar values of the other points
Vec3t contribs{
this->TriangleArea(newCoords[ord[1]], newCoords[ord[2]], newCoords[ord[3]]) / fullArea,
this->TriangleArea(newCoords[ord[0]], newCoords[ord[2]], newCoords[ord[3]]) / fullArea,
this->TriangleArea(newCoords[ord[0]], newCoords[ord[1]], newCoords[ord[3]]) / fullArea
};
Vec3t backfaceProjectionInterpolated = contribs[0] * coords[ord[0]] +
contribs[1] * coords[ord[1]] + contribs[2] * coords[ord[2]];
density[ord[3]] = this->Distance(coords[ord[3]], backfaceProjectionInterpolated) / volume;
}
else
{
// 4 triangles projection
FieldType distanceToIntersection1 = this->Distance(newCoords[4], newCoords[ord[0]]);
FieldType diagonalLength1 = this->Distance(newCoords[ord[2]], newCoords[ord[0]]);
FieldType distanceToIntersection2 = this->Distance(newCoords[4], newCoords[ord[1]]);
FieldType diagonalLength2 = this->Distance(newCoords[ord[1]], newCoords[ord[3]]);
// Spatial volume is null, or data domain surface is null
if (volume == 0.0f || diagonalLength1 == 0.0f || diagonalLength2 == 0.0f)
{
density[4] = vtkm::Infinity<FieldType>();
return;
}
// Interpolate the intersection of diagonals to get its scalar values
FieldType interpolateRatio1 = distanceToIntersection1 / diagonalLength1;
FieldType interpolateRatio2 = distanceToIntersection2 / diagonalLength2;
Vec3t interpolatedPos1 =
coords[ord[0]] + (coords[ord[2]] - coords[ord[0]]) * interpolateRatio1;
Vec3t interpolatedPos2 =
coords[ord[1]] + (coords[ord[3]] - coords[ord[1]]) * interpolateRatio2;
// Eventually, the density is calculated by dividing the scalar mass density by the volume of the cell
density[4] = this->Distance(interpolatedPos1, interpolatedPos2) / volume;
}
}
};
struct ComputeTriangles : vtkm::worklet::WorkletVisitCellsWithPoints
{
using ControlSignature = void(CellSetIn,
FieldInCell pointsOrder,
FieldInCell numberOfTris,
FieldInCell offsets,
FieldOutCell connectivity);
using ExecutionSignature =
void(_2 pointsOrder, _3 numberOfTris, _4 offsets, _5 connectivity, VisitIndex, InputIndex);
using InputDomain = _1;
using ScatterType = vtkm::worklet::ScatterCounting;
template <typename OrderType, typename ConnectivityType>
VTKM_EXEC void operator()(const OrderType& order,
const vtkm::IdComponent& numberOfTris,
const vtkm::Id& offsets,
ConnectivityType& connectivity,
const vtkm::Id& visitIndex,
const vtkm::Id& cellId) const
{
vtkm::Id secondPoint, thirdPoint;
if (numberOfTris == 3)
{
secondPoint = order[(visitIndex + 1) % 3]; // 1,2,0
// The one point in the triangle formed by the 3 other is the common vertex to all 3 visible triangles,
// And has index 3 in the order array
thirdPoint = order[3];
}
else
{
secondPoint = order[(visitIndex + 1) % 4]; // 1,2,3,0
thirdPoint =
4; // Imaginary point at the intersection of diagonals, connected to every triangle
}
connectivity[0] = cellId + offsets + order[static_cast<vtkm::IdComponent>(visitIndex)];
connectivity[1] = cellId + offsets + secondPoint;
connectivity[2] = cellId + offsets + thirdPoint;
}
};
template <typename CoordsComType,
typename CoordsComTypeOut,
typename CoordsInStorageType,
typename OutputCellSetType,
typename CoordsOutStorageType,
typename FieldType>
void Run(const vtkm::cont::CellSetSingleType<>& inputCellSet,
const vtkm::cont::ArrayHandle<vtkm::Vec<CoordsComType, 3>, CoordsInStorageType>& coords,
vtkm::cont::ArrayHandle<vtkm::Vec<CoordsComTypeOut, 3>, CoordsOutStorageType>& newCoords,
vtkm::cont::ArrayHandle<FieldType>& density,
const vtkm::cont::ArrayHandle<FieldType>& field1,
const vtkm::cont::ArrayHandle<FieldType>& field2,
OutputCellSetType& outputCellset)
{
vtkm::cont::Invoker invoke;
// Use zip to pass both scalar fields to worklets as a single argument
auto scalars = vtkm::cont::make_ArrayHandleZip(field1, field2);
// We want to project every tetrahedron in the 2-dimensional data domain using its scalar values,
// following the tetra projection algorithm
// (see "A polygonal approximation to direct scalar volume rendering" by Shirley and Tuchman)
// Minus degenerate cases, this projection makes 3 or 4 triangles in the 2D plane
// This first worklet generates the number of points and triangles needed to project a tetrahedron,
// And the order in which to take them to build the cells
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::IdComponent, 4>> pointsOrder;
vtkm::cont::ArrayHandle<vtkm::IdComponent> numberOfPoints;
vtkm::cont::ArrayHandle<vtkm::IdComponent> numberOfTris;
invoke(
ClassifyTetra<FieldType>{}, inputCellSet, scalars, pointsOrder, numberOfPoints, numberOfTris);
vtkm::Id totalPoints;
vtkm::cont::ArrayHandle<vtkm::Id> offsets =
vtkm::cont::ConvertNumComponentsToOffsets(numberOfPoints, totalPoints);
// Now, compute the tetra's coordinates in the data plane,
// and the density of each projected point
vtkm::cont::ArrayHandle<FieldType> volumeMeasure;
newCoords.Allocate(totalPoints);
density.Allocate(totalPoints);
invoke(VolumeMeasure<FieldType>{},
inputCellSet,
scalars,
coords,
numberOfTris,
pointsOrder,
vtkm::cont::make_ArrayHandleGroupVecVariable(newCoords, offsets),
vtkm::cont::make_ArrayHandleGroupVecVariable(density, offsets));
// Finally, write triangle connectivity in the data domain
vtkm::worklet::ScatterCounting scatter(numberOfTris, true);
vtkm::cont::ArrayHandle<vtkm::Id> offsets_connectivity = scatter.GetInputToOutputMap();
vtkm::cont::ArrayHandle<vtkm::Id> outConnectivity;
invoke(ComputeTriangles{},
scatter,
inputCellSet,
pointsOrder,
numberOfTris,
offsets_connectivity,
vtkm::cont::make_ArrayHandleGroupVec<3>(outConnectivity));
// Create the new dataset
outputCellset.Fill(totalPoints, vtkm::CellShapeTagTriangle::Id, 3, outConnectivity);
}
};
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_ContinuousScatterPlot_h