From 4ba15a993d2947ea2e75cdfd3d640e1531cd346c Mon Sep 17 00:00:00 2001 From: Louis Gombert Date: Tue, 21 Mar 2023 09:43:35 +0100 Subject: [PATCH] Add new ContinuousScatterPlot filter 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. --- docs/changelog/continuous-scatterplot.md | 7 + vtkm/filter/density_estimate/CMakeLists.txt | 2 + .../ContinuousScatterPlot.cxx | 82 +++ .../density_estimate/ContinuousScatterPlot.h | 60 ++ .../density_estimate/testing/CMakeLists.txt | 1 + .../testing/UnitTestContinuousScatterPlot.cxx | 690 ++++++++++++++++++ vtkm/filter/density_estimate/vtkm.module | 1 + .../density_estimate/worklet/CMakeLists.txt | 1 + .../worklet/ContinuousScatterPlot.h | 473 ++++++++++++ 9 files changed, 1317 insertions(+) create mode 100644 docs/changelog/continuous-scatterplot.md create mode 100644 vtkm/filter/density_estimate/ContinuousScatterPlot.cxx create mode 100644 vtkm/filter/density_estimate/ContinuousScatterPlot.h create mode 100644 vtkm/filter/density_estimate/testing/UnitTestContinuousScatterPlot.cxx create mode 100644 vtkm/filter/density_estimate/worklet/ContinuousScatterPlot.h diff --git a/docs/changelog/continuous-scatterplot.md b/docs/changelog/continuous-scatterplot.md new file mode 100644 index 000000000..be85f46ea --- /dev/null +++ b/docs/changelog/continuous-scatterplot.md @@ -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. diff --git a/vtkm/filter/density_estimate/CMakeLists.txt b/vtkm/filter/density_estimate/CMakeLists.txt index 2fb862448..4b73cea47 100644 --- a/vtkm/filter/density_estimate/CMakeLists.txt +++ b/vtkm/filter/density_estimate/CMakeLists.txt @@ -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 diff --git a/vtkm/filter/density_estimate/ContinuousScatterPlot.cxx b/vtkm/filter/density_estimate/ContinuousScatterPlot.cxx new file mode 100644 index 000000000..0cfc7dc59 --- /dev/null +++ b/vtkm/filter/density_estimate/ContinuousScatterPlot.cxx @@ -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 +#include +#include +#include +#include +#include +#include + +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::ValueType; + + vtkm::cont::CellSetSingleType<> scatterplotCellSet; + vtkm::cont::ArrayHandle density; + vtkm::cont::ArrayHandle> newCoords; + + // Both fields need to resolve to the same type in order to perform calculations + vtkm::cont::ArrayHandle 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 diff --git a/vtkm/filter/density_estimate/ContinuousScatterPlot.h b/vtkm/filter/density_estimate/ContinuousScatterPlot.h new file mode 100644 index 000000000..f5c978d25 --- /dev/null +++ b/vtkm/filter/density_estimate/ContinuousScatterPlot.h @@ -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 +#include + +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 diff --git a/vtkm/filter/density_estimate/testing/CMakeLists.txt b/vtkm/filter/density_estimate/testing/CMakeLists.txt index 21b8ce4d9..c3f897b57 100644 --- a/vtkm/filter/density_estimate/testing/CMakeLists.txt +++ b/vtkm/filter/density_estimate/testing/CMakeLists.txt @@ -9,6 +9,7 @@ ##============================================================================ set(unit_tests + UnitTestContinuousScatterPlot.cxx UnitTestEntropyFilter.cxx UnitTestHistogramFilter.cxx UnitTestNDEntropyFilter.cxx diff --git a/vtkm/filter/density_estimate/testing/UnitTestContinuousScatterPlot.cxx b/vtkm/filter/density_estimate/testing/UnitTestContinuousScatterPlot.cxx new file mode 100644 index 000000000..94ed0c1d8 --- /dev/null +++ b/vtkm/filter/density_estimate/testing/UnitTestContinuousScatterPlot.cxx @@ -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 +#include +#include +#include +#include +#include + +namespace +{ + +std::vector 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 tetraShape{ vtkm::CELL_SHAPE_TETRA }; +std::vector tetraIndex{ 4 }; +std::vector tetraConnectivity{ 0, 1, 2, 3 }; + +std::vector 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 multiShapes{ vtkm::CELL_SHAPE_TETRA, + vtkm::CELL_SHAPE_TETRA, + vtkm::CELL_SHAPE_TETRA }; +std::vector multiIndices{ 4, 4, 4 }; +std::vector multiConnectivity{ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 }; + +template +vtkm::cont::DataSet makeDataSet(const std::vector& ds_coords, + const std::vector& ds_shapes, + const std::vector& ds_indices, + const std::vector& 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 +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 +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 +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 +void testConnectivity(const CellSetType& cellSet, + const vtkm::cont::ArrayHandle& 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({ 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>() + .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({ 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( + { 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({ 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({ 1, 2, 0, 2, 3, 0, 3, 1, 0 }); + testConnectivity(cellSet, expectedConnectivityArray); + + // Test density values + auto density = scatterPlot.GetField("density") + .GetData() + .AsArrayHandle>() + .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({ 0, 2, 1, 2, 3, 1, 3, 0, 1 }); + testConnectivity(cellSet, expectedConnectivityArray); + + // Test density values + auto density = scatterPlot.GetField("density") + .GetData() + .AsArrayHandle>() + .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({ 0, 1, 2, 1, 3, 2, 3, 0, 2 }); + testConnectivity(cellSet, expectedConnectivityArray); + + // Test density values + auto density = scatterPlot.GetField("density") + .GetData() + .AsArrayHandle>() + .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({ 0, 1, 3, 1, 2, 3, 2, 0, 3 }); + testConnectivity(cellSet, expectedConnectivityArray); + + // Test density values + auto density = scatterPlot.GetField("density") + .GetData() + .AsArrayHandle>() + .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 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>() + .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>() + .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 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 w_shape{ vtkm::CELL_SHAPE_WEDGE }; + std::vector w_indices{ 6 }; + std::vector 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>() + .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>() + .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>() + .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>() + .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); +} diff --git a/vtkm/filter/density_estimate/vtkm.module b/vtkm/filter/density_estimate/vtkm.module index c66fc406e..e01135424 100644 --- a/vtkm/filter/density_estimate/vtkm.module +++ b/vtkm/filter/density_estimate/vtkm.module @@ -6,6 +6,7 @@ DEPENDS vtkm_filter_core PRIVATE_DEPENDS vtkm_worklet + vtkm_filter_geometry_refinement OPTIONAL_DEPENDS MPI::MPI_CXX TEST_DEPENDS diff --git a/vtkm/filter/density_estimate/worklet/CMakeLists.txt b/vtkm/filter/density_estimate/worklet/CMakeLists.txt index 2a962fc5a..f56e922d5 100644 --- a/vtkm/filter/density_estimate/worklet/CMakeLists.txt +++ b/vtkm/filter/density_estimate/worklet/CMakeLists.txt @@ -9,6 +9,7 @@ ##============================================================================ set(headers + ContinuousScatterPlot.h FieldEntropy.h FieldHistogram.h NDimsEntropy.h diff --git a/vtkm/filter/density_estimate/worklet/ContinuousScatterPlot.h b/vtkm/filter/density_estimate/worklet/ContinuousScatterPlot.h new file mode 100644 index 000000000..47ed3352f --- /dev/null +++ b/vtkm/filter/density_estimate/worklet/ContinuousScatterPlot.h @@ -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 +#include +#include +#include +#include +#include +#include +#include + +namespace vtkm +{ +namespace worklet +{ +class ContinuousScatterPlot +{ +public: + template + 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& point1, + const vtkm::Pair& point2, + const vtkm::Pair& 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 + 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 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 + 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; + + template + 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 + 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 spatialVector = { coords[1] - coords[0], + coords[2] - coords[0], + coords[3] - coords[0] }; + vtkm::Vec spatialCrossProducts = { vtkm::Cross(spatialVector[1], spatialVector[0]), + vtkm::Cross(spatialVector[0], spatialVector[2]), + vtkm::Cross(spatialVector[2], + spatialVector[1]) }; + vtkm::Vec 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 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{ 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(); + 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(); + 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 + 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(visitIndex)]; + connectivity[1] = cellId + offsets + secondPoint; + connectivity[2] = cellId + offsets + thirdPoint; + } + }; + + template + void Run(const vtkm::cont::CellSetSingleType<>& inputCellSet, + const vtkm::cont::ArrayHandle, CoordsInStorageType>& coords, + vtkm::cont::ArrayHandle, CoordsOutStorageType>& newCoords, + vtkm::cont::ArrayHandle& density, + const vtkm::cont::ArrayHandle& field1, + const vtkm::cont::ArrayHandle& 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> pointsOrder; + vtkm::cont::ArrayHandle numberOfPoints; + vtkm::cont::ArrayHandle numberOfTris; + invoke( + ClassifyTetra{}, inputCellSet, scalars, pointsOrder, numberOfPoints, numberOfTris); + + vtkm::Id totalPoints; + vtkm::cont::ArrayHandle 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 volumeMeasure; + newCoords.Allocate(totalPoints); + density.Allocate(totalPoints); + invoke(VolumeMeasure{}, + 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 offsets_connectivity = scatter.GetInputToOutputMap(); + vtkm::cont::ArrayHandle 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