Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-19 18:45:43 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge topic 'polyLinePathGeom'

Browse Source 
fc69b9c1b Merge branch 'master' of https://gitlab.kitware.com/vtk/vtk-m into polyLinePathGeom
b866b31f4 Fix cuda compile error.
023e6bb62 Remove redundant code for computing w.
e74a0800b Move helper worklets from the detail namespace into the Tube class.
fc0e1b703 Merge branch 'polyLinePathGeom' of gitlab.kitware.com:dpugmire/vtk-m into polyLinePathGeom
2dad0301e Handle polylines with only 1 point. Add testing for linear polylines.
ca34f7319 Merge branch 'polyLineExample' into 'polyLinePathGeom'
d8faed848 Polyline examples: More intelligible comments.
...

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !1710
This commit is contained in:
Dave Pugmire 2019-07-03 16:47:44 +00:00 committed by Kitware Robot
commit 5a80d7042e
8 changed files with 935 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
examples/CMakeLists.txt
View File

@ -26,6 +26,7 @@ if(VTKm_ENABLE_EXAMPLES)
add_subdirectory(multi_backend)
add_subdirectory(oscillator)
add_subdirectory(particle_advection)
add_subdirectory(polyline_archimedean_helix)
add_subdirectory(redistribute_points)
add_subdirectory(rendering)
add_subdirectory(streamline)

18
examples/polyline_archimedean_helix/CMakeLists.txt Normal file
View File

@ -0,0 +1,18 @@
##============================================================================
## 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.
##============================================================================
cmake_minimum_required(VERSION 3.8...3.14 FATAL_ERROR)
project(PolyLineArchimedeanHelix CXX)
find_package(VTKm REQUIRED QUIET)
add_executable(PolyLineArchimedeanHelix PolyLineArchimedeanHelix.cxx)
if(TARGET vtkm::cuda)
set_source_files_properties(PolyLineArchimedeanHelix.cxx PROPERTIES LANGUAGE "CUDA")
endif()
target_link_libraries(PolyLineArchimedeanHelix PRIVATE vtkm_filter vtkm_cont vtkm_rendering)

146
examples/polyline_archimedean_helix/PolylineArchimedeanHelix.cxx Normal file
View File

@ -0,0 +1,146 @@
//============================================================================
// 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 <complex>
#include <vtkm/cont/DataSetBuilderExplicit.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/io/writer/VTKDataSetWriter.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/Tube.h>
#include <vtkm/cont/ColorTable.h>
#include <vtkm/cont/CoordinateSystem.h>
#include <vtkm/cont/DataSetFieldAdd.h>
#include <vtkm/rendering/Camera.h>
#include <vtkm/rendering/CanvasRayTracer.h>
#include <vtkm/rendering/MapperRayTracer.h>
#include <vtkm/rendering/Scene.h>
#include <vtkm/rendering/View3D.h>
vtkm::Vec<vtkm::FloatDefault, 3> ArchimedeanSpiralToCartesian(
vtkm::Vec<vtkm::FloatDefault, 3> const& p)
{
// p[0] = r, p[1] = theta, p[2] = z:
vtkm::Vec<vtkm::FloatDefault, 3> xyz;
auto c = std::polar(p[0], p[1]);
xyz[0] = c.real();
xyz[1] = c.imag();
xyz[2] = p[2];
return xyz;
}
void TubeThatSpiral(vtkm::FloatDefault radius, vtkm::Id numLineSegments, vtkm::Id numSides)
{
vtkm::cont::DataSetBuilderExplicitIterative dsb;
std::vector<vtkm::Id> ids;
// The Archimedian spiral is defined by the equation r = a + b*theta.
// To extend to a 3D curve, use z = t, theta = t, r = a + b t.
vtkm::FloatDefault a = vtkm::FloatDefault(0.2);
vtkm::FloatDefault b = vtkm::FloatDefault(0.8);
for (vtkm::Id i = 0; i < numLineSegments; ++i)
{
vtkm::FloatDefault t = 4 * vtkm::FloatDefault(3.1415926) * (i + 1) /
numLineSegments; // roughly two spins around. Doesn't need to be perfect.
vtkm::FloatDefault r = a + b * t;
vtkm::FloatDefault theta = t;
vtkm::Vec<vtkm::FloatDefault, 3> cylindricalCoordinate{ r, theta, t };
vtkm::Vec<vtkm::FloatDefault, 3> spiralSample =
ArchimedeanSpiralToCartesian(cylindricalCoordinate);
vtkm::Id pid = dsb.AddPoint(spiralSample);
ids.push_back(pid);
}
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
vtkm::cont::DataSet ds = dsb.Create();
vtkm::worklet::Tube tubeWorklet(
/*capEnds = */ true,
/* how smooth the cylinder is; infinitely smooth as n->infty */ numSides,
radius);
// You added lines, but you need to extend it to a tube.
// This generates a new pointset, and new cell set.
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> tubePoints;
vtkm::cont::CellSetSingleType<> tubeCells;
tubeWorklet.Run(ds.GetCoordinateSystem(0), ds.GetCellSet(0), tubePoints, tubeCells);
vtkm::cont::DataSet tubeDataset;
tubeDataset.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coords", tubePoints));
tubeDataset.AddCellSet(tubeCells);
vtkm::Bounds coordsBounds = tubeDataset.GetCoordinateSystem().GetBounds();
vtkm::Vec<vtkm::Float64, 3> totalExtent(
coordsBounds.X.Length(), coordsBounds.Y.Length(), coordsBounds.Z.Length());
vtkm::Float64 mag = vtkm::Magnitude(totalExtent);
vtkm::Normalize(totalExtent);
// setup a camera and point it to towards the center of the input data
vtkm::rendering::Camera camera;
camera.ResetToBounds(coordsBounds);
camera.SetLookAt(totalExtent * (mag * .5f));
camera.SetViewUp(vtkm::make_Vec(0.f, 1.f, 0.f));
camera.SetClippingRange(1.f, 100.f);
camera.SetFieldOfView(60.f);
camera.SetPosition(totalExtent * (mag * 2.f));
vtkm::cont::ColorTable colorTable("inferno");
vtkm::rendering::Scene scene;
vtkm::rendering::MapperRayTracer mapper;
vtkm::rendering::CanvasRayTracer canvas(2048, 2048);
vtkm::rendering::Color bg(0.2f, 0.2f, 0.2f, 1.0f);
vtkm::cont::DataSetFieldAdd dsfa;
std::vector<vtkm::FloatDefault> v(tubePoints.GetNumberOfValues());
// The first value is a cap:
v[0] = 0;
for (vtkm::Id i = 1; i < vtkm::Id(v.size()); i += numSides)
{
vtkm::FloatDefault t = 4 * vtkm::FloatDefault(3.1415926) * (i + 1) / numSides;
vtkm::FloatDefault r = a + b * t;
for (vtkm::Id j = i; j < i + numSides && j < vtkm::Id(v.size()); ++j)
{
v[j] = r;
}
}
// Point at the end cap should be the same color as the surroundings:
v[v.size() - 1] = v[v.size() - 2];
dsfa.AddPointField(tubeDataset, "Spiral Radius", v);
scene.AddActor(vtkm::rendering::Actor(tubeDataset.GetCellSet(),
tubeDataset.GetCoordinateSystem(),
tubeDataset.GetField("Spiral Radius"),
colorTable));
vtkm::rendering::View3D view(scene, mapper, canvas, camera, bg);
view.Initialize();
view.Paint();
std::string output_filename = "tube_output_" + std::to_string(numSides) + "_sides.pnm";
view.SaveAs(output_filename);
}
int main()
{
// Radius of the tube:
vtkm::FloatDefault radius = 0.5f;
// How many segments is the tube decomposed into?
vtkm::Id numLineSegments = 100;
// As numSides->infty, the tubes becomes perfectly cylindrical:
vtkm::Id numSides = 50;
TubeThatSpiral(radius, numLineSegments, numSides);
// Setting numSides = 4 makes a square around the polyline:
numSides = 4;
TubeThatSpiral(radius, numLineSegments, numSides);
return 0;
}

6
vtkm/cont/DataSetBuilderExplicit.cxx
View File

@ -46,7 +46,8 @@ vtkm::Id DataSetBuilderExplicitIterative::AddPoint(const vtkm::Vec<vtkm::Float32
{
points.push_back(pt);
vtkm::Id id = static_cast<vtkm::Id>(points.size());
return id;
//ID is zero-based.
return id - 1;
}
VTKM_CONT
@ -56,7 +57,8 @@ vtkm::Id DataSetBuilderExplicitIterative::AddPoint(const vtkm::Float32& x,
{
points.push_back(vtkm::make_Vec(x, y, z));
vtkm::Id id = static_cast<vtkm::Id>(points.size());
return id;
//ID is zero-based.
return id - 1;
}
//Define cells.

1
vtkm/worklet/CMakeLists.txt
View File

@ -70,6 +70,7 @@ set(headers
Threshold.h
ThresholdPoints.h
Triangulate.h
Tube.h
VertexClustering.h
WarpScalar.h
WarpVector.h

539
vtkm/worklet/Tube.h Normal file
View File

@ -0,0 +1,539 @@
//============================================================================
// 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_tube_h
#define vtk_m_worklet_tube_h
#include <typeinfo>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/WorkletMapTopology.h>
namespace vtkm
{
namespace worklet
{
class Tube
{
public:
//Helper worklet to count various things in each polyline.
class CountSegments : public vtkm::worklet::WorkletMapPointToCell
{
public:
VTKM_CONT
CountSegments(const bool& capping, const vtkm::Id& n)
: Capping(capping)
, NumSides(n)
, NumVertsPerCell(3)
{
}
using ControlSignature = void(CellSetIn,
FieldOut ptsPerPolyline,
FieldOut ptsPerTube,
FieldOut numTubeConnIds);
using ExecutionSignature = void(CellShape shapeType,
PointCount numPoints,
_2 ptsPerPolyline,
_3 ptsPerTube,
_4 numTubeConnIds);
using InputDomain = _1;
template <typename CellShapeTag>
VTKM_EXEC void operator()(const CellShapeTag& shapeType,
const vtkm::IdComponent& numPoints,
vtkm::Id& ptsPerPolyline,
vtkm::Id& ptsPerTube,
vtkm::Id& numTubeConnIds) const
{
// We only support polylines that contain 2 or more points.
if (shapeType.Id == vtkm::CELL_SHAPE_POLY_LINE && numPoints > 1)
{
ptsPerPolyline = numPoints;
ptsPerTube = this->NumSides * numPoints;
// (two tris per segment) X (numSides) X numVertsPerCell
numTubeConnIds = (numPoints - 1) * 2 * this->NumSides * this->NumVertsPerCell;
//Capping adds center vertex in middle of cap, plus NumSides triangles for cap.
if (this->Capping)
{
ptsPerTube += 2;
numTubeConnIds += (2 * this->NumSides * this->NumVertsPerCell);
}
}
else
{
ptsPerPolyline = 0;
ptsPerTube = 0;
numTubeConnIds = 0;
}
}
private:
bool Capping;
vtkm::Id NumSides;
vtkm::Id NumVertsPerCell;
};
//Helper worklet to generate normals at each point in the polyline.
class GenerateNormals : public vtkm::worklet::WorkletMapPointToCell
{
static constexpr vtkm::FloatDefault vecMagnitudeEps = static_cast<vtkm::FloatDefault>(1e-3);
public:
VTKM_CONT
GenerateNormals()
: DefaultNorm(0, 0, 1)
{
}
using ControlSignature = void(CellSetIn cellset,
WholeArrayIn pointCoords,
FieldInTo polylineOffset,
WholeArrayOut newNormals);
using ExecutionSignature = void(CellShape shapeType,
PointCount numPoints,
PointIndices ptIndices,
_2 inPts,
_3 polylineOffset,
_4 outNormals);
using InputDomain = _1;
template <typename InPointsType, typename PointIndexType>
VTKM_EXEC vtkm::IdComponent FindValidSegment(const InPointsType& inPts,
const PointIndexType& ptIndices,
const vtkm::IdComponent& numPoints,
vtkm::IdComponent start) const
{
auto ps = inPts.Get(ptIndices[start]);
vtkm::IdComponent end = start + 1;
while (end < numPoints)
{
auto pe = inPts.Get(ptIndices[end]);
if (vtkm::Magnitude(pe - ps) > 0)
return end - 1;
end++;
}
return numPoints;
}
template <typename CellShapeTag,
typename PointIndexType,
typename InPointsType,
typename OutNormalType>
VTKM_EXEC void operator()(const CellShapeTag& shapeType,
const vtkm::IdComponent& numPoints,
const PointIndexType& ptIndices,
const InPointsType& inPts,
const vtkm::Id& polylineOffset,
OutNormalType& outNormals) const
{
//Ignore non-polyline and polyline with less than 2 points.
if (shapeType.Id != vtkm::CELL_SHAPE_POLY_LINE || numPoints < 2)
return;
else
{
//The following follows the VTK implementation in:
//vtkPolyLine::GenerateSlidingNormals
vtkm::Vec<vtkm::FloatDefault, 3> sPrev, sNext, normal, p0, p1;
vtkm::IdComponent sNextId = FindValidSegment(inPts, ptIndices, numPoints, 0);
if (sNextId != numPoints) // at least one valid segment
{
p0 = inPts.Get(ptIndices[sNextId]);
p1 = inPts.Get(ptIndices[sNextId + 1]);
sPrev = vtkm::Normal(p1 - p0);
}
else // no valid segments. Set everything to the default normal.
{
for (vtkm::Id i = 0; i < numPoints; i++)
outNormals.Set(polylineOffset + i, this->DefaultNorm);
return;
}
// find the next valid, non-parallel segment
while (++sNextId < numPoints)
{
sNextId = FindValidSegment(inPts, ptIndices, numPoints, sNextId);
if (sNextId != numPoints)
{
p0 = inPts.Get(ptIndices[sNextId]);
p1 = inPts.Get(ptIndices[sNextId + 1]);
sNext = vtkm::Normal(p1 - p0);
// now the starting normal should simply be the cross product
// in the following if statement we check for the case where
// the two segments are parallel, in which case, continue searching
// for the next valid segment
auto n = vtkm::Cross(sPrev, sNext);
if (vtkm::Magnitude(n) > vecMagnitudeEps)
{
normal = n;
sPrev = sNext;
break;
}
}
}
//only one valid segment...
if (sNextId >= numPoints)
{
for (vtkm::IdComponent j = 0; j < 3; j++)
if (sPrev[j] != 0)
{
normal[(j + 2) % 3] = 0;
normal[(j + 1) % 3] = 1;
normal[j] = -sPrev[(j + 1) % 3] / sPrev[j];
break;
}
}
vtkm::Normalize(normal);
vtkm::Id lastNormalId = 0;
while (++sNextId < numPoints)
{
sNextId = FindValidSegment(inPts, ptIndices, numPoints, sNextId);
if (sNextId == numPoints)
break;
p0 = inPts.Get(ptIndices[sNextId]);
p1 = inPts.Get(ptIndices[sNextId + 1]);
sNext = vtkm::Normal(p1 - p0);
auto q = vtkm::Cross(sNext, sPrev);
if (vtkm::Magnitude(q) <= vtkm::Epsilon<vtkm::FloatDefault>()) //can't use this segment
continue;
vtkm::Normalize(q);
vtkm::FloatDefault f1 = vtkm::Dot(q, normal);
vtkm::FloatDefault f2 = 1 - (f1 * f1);
if (f2 > 0)
f2 = vtkm::Sqrt(f2);
else
f2 = 0;
auto c = vtkm::Normal(sNext + sPrev);
auto w = vtkm::Cross(c, q);
c = vtkm::Cross(sPrev, q);
if ((vtkm::Dot(normal, c) * vtkm::Dot(w, c)) < 0)
f2 = -f2;
for (vtkm::Id i = lastNormalId; i < sNextId; i++)
outNormals.Set(polylineOffset + i, normal);
lastNormalId = sNextId;
sPrev = sNext;
normal = (f1 * q) + (f2 * w);
}
for (vtkm::Id i = lastNormalId; i < numPoints; i++)
outNormals.Set(polylineOffset + i, normal);
}
}
private:
vtkm::Vec<vtkm::FloatDefault, 3> DefaultNorm;
};
//Helper worklet to generate the tube points
class GeneratePoints : public vtkm::worklet::WorkletMapPointToCell
{
public:
VTKM_CONT
GeneratePoints(const bool& capping, const vtkm::Id& n, const vtkm::FloatDefault& r)
: Capping(capping)
, NumSides(n)
, Radius(r)
, Theta(2 * static_cast<vtkm::FloatDefault>(vtkm::Pi()) / static_cast<vtkm::FloatDefault>(n))
{
}
using ControlSignature = void(CellSetIn cellset,
WholeArrayIn pointCoords,
WholeArrayIn normals,
FieldInTo tubePointOffsets,
FieldInTo polylineOffset,
WholeArrayOut newPointCoords);
using ExecutionSignature = void(CellShape shapeType,
PointCount numPoints,
PointIndices ptIndices,
_2 inPts,
_3 inNormals,
_4 tubePointOffsets,
_5 polylineOffset,
_6 outPts);
using InputDomain = _1;
template <typename CellShapeTag,
typename PointIndexType,
typename InPointsType,
typename InNormalsType,
typename OutPointsType>
VTKM_EXEC void operator()(const CellShapeTag& shapeType,
const vtkm::IdComponent& numPoints,
const PointIndexType& ptIndices,
const InPointsType& inPts,
const InNormalsType& inNormals,
const vtkm::Id& tubePointOffsets,
const vtkm::Id& polylineOffset,
OutPointsType& outPts) const
{
if (shapeType.Id != vtkm::CELL_SHAPE_POLY_LINE || numPoints < 2)
return;
else
{
vtkm::Vec<vtkm::FloatDefault, 3> n, p, pNext, sNext, sPrev;
vtkm::Id outIdx = tubePointOffsets;
for (vtkm::IdComponent j = 0; j < numPoints; j++)
{
if (j == 0) //first point
{
p = inPts.Get(ptIndices[j]);
pNext = inPts.Get(ptIndices[j + 1]);
sNext = pNext - p;
sPrev = sNext;
}
else if (j == numPoints - 1) //last point
{
sPrev = sNext;
p = pNext;
}
else
{
p = pNext;
pNext = inPts.Get(ptIndices[j + 1]);
sPrev = sNext;
sNext = pNext - p;
}
n = inNormals.Get(polylineOffset + j);
//Coincident points.
if (vtkm::Magnitude(sNext) <= vtkm::Epsilon<vtkm::FloatDefault>())
this->RaiseError("Coincident points in Tube worklet.");
vtkm::Normalize(sNext);
auto s = (sPrev + sNext) / 2.;
if (vtkm::Magnitude(s) <= vtkm::Epsilon<vtkm::FloatDefault>())
s = vtkm::Cross(sPrev, n);
vtkm::Normalize(s);
auto w = vtkm::Cross(s, n);
//Bad normal
if (vtkm::Magnitude(w) <= vtkm::Epsilon<vtkm::FloatDefault>())
this->RaiseError("Bad normal in Tube worklet.");
vtkm::Normalize(w);
//create orthogonal coordinate system.
auto nP = vtkm::Cross(w, s);
vtkm::Normalize(nP);
//Add the start cap vertex. This is just a point at the center of the tube (on the polyline).
if (this->Capping && j == 0)
{
outPts.Set(outIdx, p);
outIdx++;
}
//this only implements the 'sides share vertices' line 476
vtkm::Vec<vtkm::FloatDefault, 3> normal;
for (vtkm::IdComponent k = 0; k < this->NumSides; k++)
{
vtkm::FloatDefault angle = static_cast<vtkm::FloatDefault>(k) * this->Theta;
vtkm::FloatDefault cosValue = vtkm::Cos(angle);
vtkm::FloatDefault sinValue = vtkm::Sin(angle);
normal = w * cosValue + nP * sinValue;
auto newPt = p + this->Radius * normal;
outPts.Set(outIdx, newPt);
outIdx++;
}
//Add the end cap vertex. This is just a point at the center of the tube (on the polyline).
if (this->Capping && j == numPoints - 1)
{
outPts.Set(outIdx, p);
outIdx++;
}
}
}
}
private:
bool Capping;
vtkm::Id NumSides;
vtkm::FloatDefault Radius;
vtkm::FloatDefault Theta;
};
//Helper worklet to generate the tube cells
class GenerateCells : public vtkm::worklet::WorkletMapPointToCell
{
public:
VTKM_CONT
GenerateCells(const bool& capping, const vtkm::Id& n)
: Capping(capping)
, NumSides(n)
{
}
using ControlSignature = void(CellSetIn cellset,
FieldInTo tubePointOffsets,
FieldInTo tubeConnOffsets,
WholeArrayOut outConnectivity);
using ExecutionSignature = void(CellShape shapeType,
PointCount numPoints,
_2 tubePointOffset,
_3 tubeConnOffsets,
_4 outConn);
using InputDomain = _1;
template <typename CellShapeTag, typename OutConnType>
VTKM_EXEC void operator()(const CellShapeTag& shapeType,
const vtkm::IdComponent& numPoints,
const vtkm::Id& tubePointOffset,
const vtkm::Id& tubeConnOffset,
OutConnType& outConn) const
{
if (shapeType.Id != vtkm::CELL_SHAPE_POLY_LINE || numPoints < 2)
return;
else
{
vtkm::Id outIdx = tubeConnOffset;
vtkm::Id tubePtOffset = (this->Capping ? tubePointOffset + 1 : tubePointOffset);
for (vtkm::IdComponent i = 0; i < numPoints - 1; i++)
{
for (vtkm::Id j = 0; j < this->NumSides; j++)
{
//Triangle 1: verts 0,1,2
outConn.Set(outIdx + 0, tubePtOffset + i * this->NumSides + j);
outConn.Set(outIdx + 1, tubePtOffset + i * this->NumSides + (j + 1) % this->NumSides);
outConn.Set(outIdx + 2,
tubePtOffset + (i + 1) * this->NumSides + (j + 1) % this->NumSides);
outIdx += 3;
//Triangle 2: verts 0,2,3
outConn.Set(outIdx + 0, tubePtOffset + i * this->NumSides + j);
outConn.Set(outIdx + 1,
tubePtOffset + (i + 1) * this->NumSides + (j + 1) % this->NumSides);
outConn.Set(outIdx + 2, tubePtOffset + (i + 1) * this->NumSides + j);
outIdx += 3;
}
}
if (this->Capping)
{
//start cap triangles
vtkm::Id startCenterPt = 0 + tubePointOffset;
for (vtkm::Id j = 0; j < this->NumSides; j++)
{
outConn.Set(outIdx + 0, startCenterPt);
outConn.Set(outIdx + 1, startCenterPt + 1 + j);
outConn.Set(outIdx + 2, startCenterPt + 1 + ((j + 1) % this->NumSides));
outIdx += 3;
}
//end cap triangles
vtkm::Id endCenterPt = (tubePointOffset + 1) + (numPoints * this->NumSides);
vtkm::Id endOffsetPt = endCenterPt - this->NumSides;
for (vtkm::Id j = 0; j < this->NumSides; j++)
{
outConn.Set(outIdx + 0, endCenterPt);
outConn.Set(outIdx + 1, endOffsetPt + j);
outConn.Set(outIdx + 2, endOffsetPt + ((j + 1) % this->NumSides));
outIdx += 3;
}
}
}
}
private:
bool Capping;
vtkm::Id NumSides;
};
VTKM_CONT
Tube(const bool& capping, const vtkm::Id& n, const vtkm::FloatDefault& r)
: Capping(capping)
, NumSides(n)
, Radius(r)
{
}
VTKM_CONT
void Run(const vtkm::cont::CoordinateSystem& coords,
const vtkm::cont::DynamicCellSet& cellset,
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>>& newPoints,
vtkm::cont::CellSetSingleType<>& newCells)
{
using ExplCoordsType = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>>;
using NormalsType = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>>;
if (!(coords.GetData().IsType<ExplCoordsType>() &&
(cellset.IsSameType(vtkm::cont::CellSetExplicit<>()) ||
cellset.IsSameType(vtkm::cont::CellSetSingleType<>()))))
{
throw vtkm::cont::ErrorBadValue("Tube filter only supported for polyline data.");
}
//Count number of polyline pts, tube pts and tube cells
vtkm::cont::ArrayHandle<vtkm::Id> ptsPerPolyline, ptsPerTube, numTubeConnIds;
CountSegments countSegs(this->Capping, this->NumSides);
vtkm::worklet::DispatcherMapTopology<CountSegments> countInvoker(countSegs);
countInvoker.Invoke(cellset, ptsPerPolyline, ptsPerTube, numTubeConnIds);
vtkm::Id totalPolylinePts = vtkm::cont::Algorithm::Reduce(ptsPerPolyline, vtkm::Id(0));
if (totalPolylinePts == 0)
throw vtkm::cont::ErrorBadValue("Tube filter only supported for polyline data.");
vtkm::Id totalTubePts = vtkm::cont::Algorithm::Reduce(ptsPerTube, vtkm::Id(0));
vtkm::Id totalTubeConnIds = vtkm::cont::Algorithm::Reduce(numTubeConnIds, vtkm::Id(0));
vtkm::cont::ArrayHandle<vtkm::Id> polylineOffset, tubePointOffsets, tubeConnOffsets;
vtkm::cont::Algorithm::ScanExclusive(ptsPerPolyline, polylineOffset);
vtkm::cont::Algorithm::ScanExclusive(ptsPerTube, tubePointOffsets);
vtkm::cont::Algorithm::ScanExclusive(numTubeConnIds, tubeConnOffsets);
//Generate normals at each point on all polylines
ExplCoordsType inCoords = coords.GetData().Cast<ExplCoordsType>();
NormalsType normals;
normals.Allocate(totalPolylinePts);
vtkm::worklet::DispatcherMapTopology<GenerateNormals> genNormalsDisp;
genNormalsDisp.Invoke(cellset, inCoords, polylineOffset, normals);
//Generate the tube points
newPoints.Allocate(totalTubePts);
GeneratePoints genPts(this->Capping, this->NumSides, this->Radius);
vtkm::worklet::DispatcherMapTopology<GeneratePoints> genPtsDisp(genPts);
genPtsDisp.Invoke(cellset, inCoords, normals, tubePointOffsets, polylineOffset, newPoints);
//Generate tube cells
vtkm::cont::ArrayHandle<vtkm::Id> newConnectivity;
newConnectivity.Allocate(totalTubeConnIds);
GenerateCells genCells(this->Capping, this->NumSides);
vtkm::worklet::DispatcherMapTopology<GenerateCells> genCellsDisp(genCells);
genCellsDisp.Invoke(cellset, tubePointOffsets, tubeConnOffsets, newConnectivity);
newCells.Fill(totalTubePts, vtkm::CELL_SHAPE_TRIANGLE, 3, newConnectivity);
}
private:
bool Capping;
vtkm::Id NumSides;
vtkm::FloatDefault Radius;
};
}
}
#endif // vtk_m_worklet_tube_h

1
vtkm/worklet/testing/CMakeLists.txt
View File

@ -65,6 +65,7 @@ set(unit_tests
UnitTestThreshold.cxx
UnitTestThresholdPoints.cxx
UnitTestTriangulate.cxx
UnitTestTube.cxx
UnitTestWholeCellSetIn.cxx
UnitTestWorkletMapField.cxx
UnitTestWorkletMapFieldExecArg.cxx

225
vtkm/worklet/testing/UnitTestTube.cxx Normal file
View File

@ -0,0 +1,225 @@
//============================================================================
// 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/DataSetBuilderExplicit.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/io/writer/VTKDataSetWriter.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/Tube.h>
namespace
{
void appendPts(vtkm::cont::DataSetBuilderExplicitIterative& dsb,
const vtkm::Vec<vtkm::FloatDefault, 3>& pt,
std::vector<vtkm::Id>& ids)
{
vtkm::Id pid = dsb.AddPoint(pt);
ids.push_back(pid);
}
void createNonPoly(vtkm::cont::DataSetBuilderExplicitIterative& dsb)
{
std::vector<vtkm::Id> ids;
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(0, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(1, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(1, 1, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_TRIANGLE, ids);
}
vtkm::Id calcNumPoints(const std::size_t& numPtIds, const vtkm::Id& numSides, const bool& capEnds)
{
//there are 'numSides' points for each polyline vertex
//plus, 2 more for the center point of start and end caps.
return static_cast<vtkm::Id>(numPtIds) * numSides + (capEnds ? 2 : 0);
}
vtkm::Id calcNumCells(const std::size_t& numPtIds, const vtkm::Id& numSides, const bool& capEnds)
{
//Each line segment has numSides * 2 triangles.
//plus, numSides triangles for each cap.
return (2 * static_cast<vtkm::Id>(numPtIds - 1) * numSides) + (capEnds ? 2 * numSides : 0);
}
void TestTube(bool capEnds, vtkm::FloatDefault radius, vtkm::Id numSides, vtkm::Id insertNonPolyPos)
{
using VecType = vtkm::Vec<vtkm::FloatDefault, 3>;
vtkm::cont::DataSetBuilderExplicitIterative dsb;
std::vector<vtkm::Id> ids;
if (insertNonPolyPos == 0)
createNonPoly(dsb);
vtkm::Id reqNumPts = 0, reqNumCells = 0;
ids.clear();
appendPts(dsb, VecType(0, 0, 0), ids);
appendPts(dsb, VecType(1, 0, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
reqNumPts += calcNumPoints(ids.size(), numSides, capEnds);
reqNumCells += calcNumCells(ids.size(), numSides, capEnds);
if (insertNonPolyPos == 1)
createNonPoly(dsb);
ids.clear();
appendPts(dsb, VecType(0, 0, 0), ids);
appendPts(dsb, VecType(1, 0, 0), ids);
appendPts(dsb, VecType(2, 0, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
reqNumPts += calcNumPoints(ids.size(), numSides, capEnds);
reqNumCells += calcNumCells(ids.size(), numSides, capEnds);
if (insertNonPolyPos == 2)
createNonPoly(dsb);
ids.clear();
appendPts(dsb, VecType(0, 0, 0), ids);
appendPts(dsb, VecType(1, 0, 0), ids);
appendPts(dsb, VecType(2, 1, 0), ids);
appendPts(dsb, VecType(3, 0, 0), ids);
appendPts(dsb, VecType(4, 0, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
reqNumPts += calcNumPoints(ids.size(), numSides, capEnds);
reqNumCells += calcNumCells(ids.size(), numSides, capEnds);
if (insertNonPolyPos == 3)
createNonPoly(dsb);
//Add something a little more complicated...
ids.clear();
vtkm::FloatDefault x0 = 0;
vtkm::FloatDefault x1 = static_cast<vtkm::FloatDefault>(6.28);
vtkm::FloatDefault dx = static_cast<vtkm::FloatDefault>(0.05);
for (vtkm::FloatDefault x = x0; x < x1; x += dx)
appendPts(dsb, VecType(x, vtkm::Cos(x), vtkm::Sin(x) / 2), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
reqNumPts += calcNumPoints(ids.size(), numSides, capEnds);
reqNumCells += calcNumCells(ids.size(), numSides, capEnds);
if (insertNonPolyPos == 4)
createNonPoly(dsb);
//Finally, add a degenerate polyline: don't dance with the beast.
ids.clear();
appendPts(dsb, VecType(6, 6, 6), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
//Should NOT produce a tubed polyline, so don't increment reqNumPts and reqNumCells.
vtkm::cont::DataSet ds = dsb.Create();
vtkm::worklet::Tube tubeWorklet(capEnds, numSides, radius);
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> newPoints;
vtkm::cont::CellSetSingleType<> newCells;
tubeWorklet.Run(ds.GetCoordinateSystem(0), ds.GetCellSet(0), newPoints, newCells);
VTKM_TEST_ASSERT(newPoints.GetNumberOfValues() == reqNumPts,
"Wrong number of points in Tube worklet");
VTKM_TEST_ASSERT(newCells.GetNumberOfCells() == reqNumCells, "Wrong cell shape in Tube worklet");
VTKM_TEST_ASSERT(newCells.GetCellShape(0) == vtkm::CELL_SHAPE_TRIANGLE,
"Wrong cell shape in Tube worklet");
}
void TestLinearPolylines()
{
using VecType = vtkm::Vec<vtkm::FloatDefault, 3>;
//Create a number of linear polylines along a set of directions.
//We check that the tubes are all copacetic (proper number of cells, points),
//and that the tube points all lie in the proper plane.
//This will validate the code that computes the coordinate frame at each
//vertex in the polyline. There are numeric checks to handle co-linear segments.
std::vector<VecType> dirs;
for (int i = -1; i <= 1; i++)
for (int j = -1; j <= 1; j++)
for (int k = -1; k <= 1; k++)
{
if (!i && !j && !k)
continue;
dirs.push_back(vtkm::Normal(VecType(static_cast<vtkm::FloatDefault>(i),
static_cast<vtkm::FloatDefault>(j),
static_cast<vtkm::FloatDefault>(k))));
}
bool capEnds = false;
vtkm::Id numSides = 3;
vtkm::FloatDefault radius = 1;
for (auto& dir : dirs)
{
vtkm::cont::DataSetBuilderExplicitIterative dsb;
std::vector<vtkm::Id> ids;
VecType pt(0, 0, 0);
for (int i = 0; i < 5; i++)
{
appendPts(dsb, pt, ids);
pt = pt + dir;
}
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
vtkm::cont::DataSet ds = dsb.Create();
vtkm::Id reqNumPts = calcNumPoints(ids.size(), numSides, capEnds);
vtkm::Id reqNumCells = calcNumCells(ids.size(), numSides, capEnds);
vtkm::worklet::Tube tubeWorklet(capEnds, numSides, radius);
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> newPoints;
vtkm::cont::CellSetSingleType<> newCells;
tubeWorklet.Run(ds.GetCoordinateSystem(0), ds.GetCellSet(0), newPoints, newCells);
VTKM_TEST_ASSERT(newPoints.GetNumberOfValues() == reqNumPts,
"Wrong number of points in Tube worklet");
VTKM_TEST_ASSERT(newCells.GetNumberOfCells() == reqNumCells,
"Wrong cell shape in Tube worklet");
VTKM_TEST_ASSERT(newCells.GetCellShape(0) == vtkm::CELL_SHAPE_TRIANGLE,
"Wrong cell shape in Tube worklet");
//Each of the 3 points should be in the plane defined by dir.
auto portal = newPoints.GetPortalConstControl();
for (vtkm::Id i = 0; i < newPoints.GetNumberOfValues(); i += 3)
{
auto p0 = portal.Get(i + 0);
auto p1 = portal.Get(i + 1);
auto p2 = portal.Get(i + 2);
auto vec = vtkm::Normal(vtkm::Cross(p0 - p1, p0 - p2));
vtkm::FloatDefault dp = vtkm::Abs(vtkm::Dot(vec, dir));
VTKM_TEST_ASSERT((1 - dp) <= vtkm::Epsilon<vtkm::FloatDefault>(),
"Tube points in wrong orientation");
}
}
}
void TestTubeWorklets()
{
std::cout << "Testing Tube Worklet" << std::endl;
std::vector<vtkm::Id> testNumSides = { 3, 4, 8, 13, 20 };
std::vector<vtkm::FloatDefault> testRadii = { 0.01f, 0.05f, 0.10f };
std::vector<int> insertNonPolylinePos = { -1, 0, 1, 2, 3, 4 };
for (auto& i : insertNonPolylinePos)
for (auto& n : testNumSides)
for (auto& r : testRadii)
{
TestTube(false, r, n, i);
TestTube(true, r, n, i);
}
TestLinearPolylines();
}
}
int UnitTestTube(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestTubeWorklets, argc, argv);
}