Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 21:33:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

updates to the connectivity tracer

Browse Source
This commit is contained in:
Matt Larsen 2018-10-30 07:46:42 -07:00
parent 8237db209e
commit 5a1685d8aa
28 changed files with 1996 additions and 2690 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
vtkm/rendering/CMakeLists.txt
View File

@ -178,6 +178,7 @@ set(device_sources
raytracing/SphereExtractor.cxx
raytracing/SphereIntersector.cxx
raytracing/TriangleExtractor.cxx
raytracing/TriangleIntersector.cxx
raytracing/VolumeRendererStructured.cxx
)

108
vtkm/rendering/ConnectivityProxy.cxx
View File

@ -68,6 +68,9 @@ public:
~InternalsType() {}
VTKM_CONT
void SetUnitScalar(vtkm::Float32 unitScalar) { Tracer.SetUnitScalar(unitScalar); }
void SetSampleDistance(const vtkm::Float32& distance)
{
if (Mode != VOLUME_MODE)
@ -119,7 +122,12 @@ public:
vtkm::Bounds GetSpatialBounds() const { return SpatialBounds; }
VTKM_CONT
vtkm::Range GetScalarRange() const { return ScalarRange; }
vtkm::Range GetScalarFieldRange()
{
const vtkm::cont::ArrayHandle<vtkm::Range> range = this->ScalarField.GetRange();
ScalarRange = range.GetPortalConstControl().Get(0);
return ScalarRange;
}
VTKM_CONT
void SetScalarRange(const vtkm::Range& range) { ScalarRange = range; }
@ -141,7 +149,7 @@ public:
this->EmissionField);
}
Tracer.Trace(rays);
Tracer.FullTrace(rays);
}
VTKM_CONT
@ -160,7 +168,46 @@ public:
this->EmissionField);
}
Tracer.Trace(rays);
Tracer.FullTrace(rays);
}
VTKM_CONT
PartialVector64 PartialTrace(vtkm::rendering::raytracing::Ray<vtkm::Float64>& rays)
{
if (Mode == VOLUME_MODE)
{
Tracer.SetVolumeData(this->ScalarField, this->ScalarRange, this->Cells, this->Coords);
}
else
{
Tracer.SetEnergyData(this->ScalarField,
rays.Buffers.at(0).GetNumChannels(),
this->Cells,
this->Coords,
this->EmissionField);
}
return Tracer.PartialTrace(rays);
}
VTKM_CONT
PartialVector32 PartialTrace(vtkm::rendering::raytracing::Ray<vtkm::Float32>& rays)
{
if (Mode == VOLUME_MODE)
{
Tracer.SetVolumeData(this->ScalarField, this->ScalarRange, this->Cells, this->Coords);
}
else
{
Tracer.SetEnergyData(this->ScalarField,
rays.Buffers.at(0).GetNumChannels(),
this->Cells,
this->Coords,
this->EmissionField);
}
return Tracer.PartialTrace(rays);
}
VTKM_CONT
@ -188,7 +235,7 @@ public:
throw vtkm::cont::ErrorBadValue("ENERGY MODE Not implemented for this use case\n");
}
Tracer.Trace(rays);
Tracer.FullTrace(rays);
canvas->WriteToCanvas(rays, rays.Buffers.at(0).Buffer, camera);
if (CompositeBackground)
@ -266,9 +313,9 @@ vtkm::Bounds ConnectivityProxy::GetSpatialBounds()
}
VTKM_CONT
vtkm::Range ConnectivityProxy::GetScalarRange()
vtkm::Range ConnectivityProxy::GetScalarFieldRange()
{
return Internals->GetScalarRange();
return Internals->GetScalarFieldRange();
}
VTKM_CONT
@ -301,6 +348,27 @@ void ConnectivityProxy::Trace(vtkm::rendering::raytracing::Ray<vtkm::Float64>& r
logger->CloseLogEntry(-1.0);
}
VTKM_CONT
PartialVector32 ConnectivityProxy::PartialTrace(
vtkm::rendering::raytracing::Ray<vtkm::Float32>& rays)
{
raytracing::Logger* logger = raytracing::Logger::GetInstance();
logger->OpenLogEntry("connectivity_trace_32");
if (Internals->GetRenderMode() == VOLUME_MODE)
{
logger->AddLogData("volume_mode", "true");
}
else
{
logger->AddLogData("volume_mode", "false");
}
PartialVector32 res = Internals->PartialTrace(rays);
logger->CloseLogEntry(-1.0);
return res;
}
VTKM_CONT
void ConnectivityProxy::Trace(vtkm::rendering::raytracing::Ray<vtkm::Float32>& rays)
{
@ -320,6 +388,27 @@ void ConnectivityProxy::Trace(vtkm::rendering::raytracing::Ray<vtkm::Float32>& r
logger->CloseLogEntry(-1.0);
}
VTKM_CONT
PartialVector64 ConnectivityProxy::PartialTrace(
vtkm::rendering::raytracing::Ray<vtkm::Float64>& rays)
{
raytracing::Logger* logger = raytracing::Logger::GetInstance();
logger->OpenLogEntry("connectivity_trace_64");
if (Internals->GetRenderMode() == VOLUME_MODE)
{
logger->AddLogData("volume_mode", "true");
}
else
{
logger->AddLogData("volume_mode", "false");
}
PartialVector64 res = Internals->PartialTrace(rays);
logger->CloseLogEntry(-1.0);
return res;
}
VTKM_CONT
void ConnectivityProxy::Trace(const vtkm::rendering::Camera& camera,
vtkm::rendering::CanvasRayTracer* canvas)
@ -332,10 +421,17 @@ void ConnectivityProxy::Trace(const vtkm::rendering::Camera& camera,
logger->CloseLogEntry(-1.0);
}
VTKM_CONT
void ConnectivityProxy::SetDebugPrints(bool on)
{
Internals->SetDebugPrints(on);
}
VTKM_CONT
void ConnectivityProxy::SetUnitScalar(vtkm::Float32 unitScalar)
{
Internals->SetUnitScalar(unitScalar);
}
}
} // namespace vtkm::rendering

10
vtkm/rendering/ConnectivityProxy.h
View File

@ -26,6 +26,7 @@
#include <vtkm/rendering/Mapper.h>
#include <vtkm/rendering/View.h>
#include <vtkm/rendering/raytracing/Camera.h>
#include <vtkm/rendering/raytracing/PartialComposite.h>
#include <vtkm/rendering/raytracing/Ray.h>
namespace vtkm
@ -33,6 +34,9 @@ namespace vtkm
namespace rendering
{
using PartialVector64 = std::vector<vtkm::rendering::raytracing::PartialComposite<vtkm::Float64>>;
using PartialVector32 = std::vector<vtkm::rendering::raytracing::PartialComposite<vtkm::Float32>>;
class VTKM_RENDERING_EXPORT ConnectivityProxy
{
public:
@ -57,14 +61,18 @@ public:
void SetColorMap(vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>>& colormap);
void SetCompositeBackground(bool on);
void SetDebugPrints(bool on);
void SetUnitScalar(vtkm::Float32 unitScalar);
vtkm::Bounds GetSpatialBounds();
vtkm::Range GetScalarRange();
vtkm::Range GetScalarFieldRange();
void Trace(const vtkm::rendering::Camera& camera, vtkm::rendering::CanvasRayTracer* canvas);
void Trace(vtkm::rendering::raytracing::Ray<vtkm::Float64>& rays);
void Trace(vtkm::rendering::raytracing::Ray<vtkm::Float32>& rays);
PartialVector64 PartialTrace(vtkm::rendering::raytracing::Ray<vtkm::Float64>& rays);
PartialVector32 PartialTrace(vtkm::rendering::raytracing::Ray<vtkm::Float32>& rays);
protected:
struct InternalsType;
struct BoundsFunctor;

4
vtkm/rendering/MapperCylinder.cxx
View File

@ -29,12 +29,8 @@
#include <vtkm/rendering/raytracing/CylinderExtractor.h>
#include <vtkm/rendering/raytracing/CylinderIntersector.h>
#include <vtkm/rendering/raytracing/Logger.h>
#include <vtkm/rendering/raytracing/QuadExtractor.h>
#include <vtkm/rendering/raytracing/QuadIntersector.h>
#include <vtkm/rendering/raytracing/RayOperations.h>
#include <vtkm/rendering/raytracing/RayTracer.h>
#include <vtkm/rendering/raytracing/SphereExtractor.h>
#include <vtkm/rendering/raytracing/SphereIntersector.h>
#include <vtkm/rendering/raytracing/Worklets.h>
namespace vtkm

2
vtkm/rendering/MapperPoint.cxx
View File

@ -21,7 +21,6 @@
#include <vtkm/rendering/MapperPoint.h>
#include <vtkm/cont/Timer.h>
#include <vtkm/cont/TryExecute.h>
#include <vtkm/rendering/CanvasRayTracer.h>
#include <vtkm/rendering/internal/RunTriangulator.h>
@ -31,7 +30,6 @@
#include <vtkm/rendering/raytracing/RayTracer.h>
#include <vtkm/rendering/raytracing/SphereExtractor.h>
#include <vtkm/rendering/raytracing/SphereIntersector.h>
#include <vtkm/rendering/raytracing/TriangleExtractor.h>
namespace vtkm
{

4
vtkm/rendering/MapperQuad.cxx
View File

@ -26,15 +26,11 @@
#include <vtkm/rendering/CanvasRayTracer.h>
#include <vtkm/rendering/Cylinderizer.h>
#include <vtkm/rendering/raytracing/Camera.h>
#include <vtkm/rendering/raytracing/CylinderExtractor.h>
#include <vtkm/rendering/raytracing/CylinderIntersector.h>
#include <vtkm/rendering/raytracing/Logger.h>
#include <vtkm/rendering/raytracing/QuadExtractor.h>
#include <vtkm/rendering/raytracing/QuadIntersector.h>
#include <vtkm/rendering/raytracing/RayOperations.h>
#include <vtkm/rendering/raytracing/RayTracer.h>
#include <vtkm/rendering/raytracing/SphereExtractor.h>
#include <vtkm/rendering/raytracing/SphereIntersector.h>
namespace vtkm
{

62
vtkm/rendering/raytracing/BVHTraverser.h
View File

@ -78,28 +78,13 @@ VTKM_EXEC inline bool IntersectAABB(const BVHPortalType& bvh,
return (min0 > min1);
}
template <template <typename> class LeafIntersectorType>
class BVHTraverser
{
public:
template <typename Device>
class Intersector : public vtkm::worklet::WorkletMapField
{
public:
typedef typename vtkm::cont::ArrayHandle<Vec<vtkm::Float32, 4>> Float4ArrayHandle;
typedef typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Int32, 2>> Int2Handle;
typedef typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> Id4Handle;
typedef typename vtkm::cont::ArrayHandle<vtkm::Id> IdHandle;
typedef typename Float4ArrayHandle::ExecutionTypes<Device>::PortalConst Float4ArrayPortal;
typedef typename Int2Handle::ExecutionTypes<Device>::PortalConst Int2ArrayPortal;
typedef typename IdHandle::ExecutionTypes<Device>::PortalConst IdArrayPortal;
typedef typename Id4Handle::ExecutionTypes<Device>::PortalConst Id4ArrayPortal;
private:
bool Occlusion;
Float4ArrayPortal FlatBVH;
IdArrayPortal Leafs;
LeafIntersectorType<Device> LeafIntersector;
VTKM_EXEC
inline vtkm::Float32 rcp(vtkm::Float32 f) const { return 1.0f / f; }
@ -118,11 +103,8 @@ public:
public:
VTKM_CONT
Intersector(bool occlusion, LinearBVH& bvh, LeafIntersectorType<Device>& leafIntersector)
Intersector(bool occlusion)
: Occlusion(occlusion)
, FlatBVH(bvh.FlatBVH.PrepareForInput(Device()))
, Leafs(bvh.Leafs.PrepareForInput(Device()))
, LeafIntersector(leafIntersector)
{
}
using ControlSignature = void(FieldIn<>,
@ -133,11 +115,18 @@ public:
FieldOut<>,
FieldOut<>,
FieldOut<>,
WholeArrayIn<Vec3RenderingTypes>);
using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, _9);
WholeArrayIn<Vec3RenderingTypes>,
ExecObject leafIntersector,
WholeArrayIn<>,
WholeArrayIn<>);
using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12);
template <typename PointPortalType, typename Precision>
template <typename PointPortalType,
typename Precision,
typename LeafType,
typename InnerNodePortalType,
typename LeafPortalType>
VTKM_EXEC void operator()(const vtkm::Vec<Precision, 3>& dir,
const vtkm::Vec<Precision, 3>& origin,
Precision& distance,
@ -146,7 +135,10 @@ public:
Precision& minU,
Precision& minV,
vtkm::Id& hitIndex,
const PointPortalType& points) const
const PointPortalType& points,
LeafType& leafIntersector,
const InnerNodePortalType& flatBVH,
const LeafPortalType& leafs) const
{
Precision closestDistance = maxDistance;
distance = maxDistance;
@ -174,7 +166,7 @@ public:
bool hitLeftChild, hitRightChild;
bool rightCloser = IntersectAABB(FlatBVH,
bool rightCloser = IntersectAABB(flatBVH,
currentNode,
originDir,
invDir,
@ -191,7 +183,7 @@ public:
else
{
vtkm::Vec<vtkm::Float32, 4> children =
FlatBVH.Get(currentNode + 3); //Children.Get(currentNode);
flatBVH.Get(currentNode + 3); //Children.Get(currentNode);
vtkm::Int32 leftChild;
memcpy(&leftChild, &children[0], 4);
vtkm::Int32 rightChild;
@ -217,7 +209,7 @@ public:
if (currentNode < 0 && currentNode != barrier) //check register usage
{
currentNode = -currentNode - 1; //swap the neg address
LeafIntersector.IntersectLeaf(currentNode,
leafIntersector.IntersectLeaf(currentNode,
origin,
dir,
points,
@ -225,7 +217,7 @@ public:
closestDistance,
minU,
minV,
Leafs,
leafs,
minDistance);
currentNode = todo[stackptr];
stackptr--;
@ -239,16 +231,13 @@ public:
};
template <typename Precision, typename Device>
template <typename Precision, typename LeafIntersectorType>
VTKM_CONT void IntersectRays(Ray<Precision>& rays,
LinearBVH& bvh,
LeafIntersectorType<Device> leafIntersector,
vtkm::cont::CoordinateSystem& coordsHandle,
Device vtkmNotUsed(Device))
LeafIntersectorType& leafIntersector,
vtkm::cont::CoordinateSystem& coordsHandle)
{
vtkm::worklet::DispatcherMapField<Intersector<Device>> intersectDispatch(
Intersector<Device>(false, bvh, leafIntersector));
intersectDispatch.SetDevice(Device());
vtkm::worklet::DispatcherMapField<Intersector> intersectDispatch(Intersector(false));
intersectDispatch.Invoke(rays.Dir,
rays.Origin,
rays.Distance,
@ -257,7 +246,10 @@ public:
rays.U,
rays.V,
rays.HitIdx,
coordsHandle);
coordsHandle,
leafIntersector,
bvh.FlatBVH,
bvh.Leafs);
}
}; // BVHTraverser
#undef END_FLAG

2
vtkm/rendering/raytracing/CMakeLists.txt
View File

@ -35,6 +35,7 @@ set(headers
MeshConnectivityContainers.h
MeshConnectivityBase.h
MortonCodes.h
PartialComposite.h
QuadExtractor.h
QuadIntersector.h
Ray.h
@ -46,6 +47,7 @@ set(headers
SphereExtractor.h
SphereIntersector.h
TriangleExtractor.h
TriangleIntersections.h
TriangleIntersector.h
VolumeRendererStructured.h
Worklets.h

3
vtkm/rendering/raytracing/CellIntersector.h
View File

@ -20,8 +20,9 @@
#ifndef vtk_m_rendering_raytracing_Cell_Intersector_h
#define vtk_m_rendering_raytracing_Cell_Intersector_h
#include <vtkm/CellShape.h>
#include <vtkm/rendering/raytracing/CellTables.h>
#include <vtkm/rendering/raytracing/TriangleIntersector.h>
#include <vtkm/rendering/raytracing/TriangleIntersections.h>
namespace vtkm
{

9
vtkm/rendering/raytracing/ChannelBuffer.cxx
View File

@ -511,6 +511,15 @@ void ChannelBuffer<Precision>::SetNumChannels(const vtkm::Int32 numChannels)
ResizeChannelFunctor<Precision> functor(this, numChannels);
vtkm::cont::TryExecute(functor);
}
template <typename Precision>
ChannelBuffer<Precision> ChannelBuffer<Precision>::Copy()
{
ChannelBuffer res(NumChannels, Size);
vtkm::cont::Algorithm::Copy(this->Buffer, res.Buffer);
return res;
}
// Instantiate supported types
template class ChannelBuffer<vtkm::Float32>;
template class ChannelBuffer<vtkm::Float64>;

2
vtkm/rendering/raytracing/ChannelBuffer.h
View File

@ -83,6 +83,8 @@ public:
const vtkm::Id outputSize,
Precision initValue = 1.f);
ChannelBuffer<Precision> Copy();
void InitConst(const Precision value);
void InitChannels(const vtkm::cont::ArrayHandle<Precision>& signature);
void Normalize(bool invert);

12
vtkm/rendering/raytracing/ChannelBufferOperations.h
View File

@ -103,23 +103,21 @@ public:
class ChannelBufferOperations
{
public:
template <typename Device, typename Precision>
template <typename Precision>
static void Compact(ChannelBuffer<Precision>& buffer,
vtkm::cont::ArrayHandle<UInt8>& masks,
const vtkm::Id& newSize,
Device)
const vtkm::Id& newSize)
{
vtkm::cont::ArrayHandle<vtkm::Id> offsets;
offsets.PrepareForOutput(buffer.Size, Device());
offsets.Allocate(buffer.Size);
vtkm::cont::ArrayHandleCast<vtkm::Id, vtkm::cont::ArrayHandle<vtkm::UInt8>> castedMasks(masks);
vtkm::cont::DeviceAdapterAlgorithm<Device>::ScanExclusive(castedMasks, offsets);
vtkm::cont::Algorithm::ScanExclusive(castedMasks, offsets);
vtkm::cont::ArrayHandle<Precision> compactedBuffer;
compactedBuffer.PrepareForOutput(newSize * buffer.NumChannels, Device());
compactedBuffer.Allocate(newSize * buffer.NumChannels);
vtkm::worklet::DispatcherMapField<detail::CompactBuffer> dispatcher(
detail::CompactBuffer(buffer.NumChannels));
dispatcher.SetDevice(Device());
dispatcher.Invoke(masks, buffer.Buffer, offsets, compactedBuffer);
buffer.Buffer = compactedBuffer;
buffer.Size = newSize;

1301
vtkm/rendering/raytracing/ConnectivityTracer.cxx
View File

File diff suppressed because it is too large Load Diff

108
vtkm/rendering/raytracing/ConnectivityTracer.h
View File

@ -22,18 +22,11 @@
#include <vtkm/rendering/vtkm_rendering_export.h>
#include <iomanip>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/rendering/raytracing/MeshConnectivityContainers.h>
#include <vtkm/rendering/raytracing/PartialComposite.h>
#ifndef CELL_SHAPE_ZOO
#define CELL_SHAPE_ZOO 255
#endif
#ifndef CELL_SHAPE_STRUCTURED
#define CELL_SHAPE_STRUCTURED 254
#endif
namespace vtkm
{
@ -68,26 +61,29 @@ public:
ExitDist = &Distance2;
}
template <typename Device>
void Compact(vtkm::cont::ArrayHandle<FloatType>& compactedDistances,
vtkm::cont::ArrayHandle<UInt8>& masks,
Device);
vtkm::cont::ArrayHandle<UInt8>& masks);
template <typename Device>
void Init(const vtkm::Id size, vtkm::cont::ArrayHandle<FloatType>& distances, Device);
void Init(const vtkm::Id size, vtkm::cont::ArrayHandle<FloatType>& distances);
void Swap();
};
} //namespace detail
class ConnectivityTracer
/**
* \brief ConnectivityTracer is volumetric ray tracer for unstructured
* grids. Capabilities include volume rendering and integrating
* absorption and emission of N energy groups for simulated
* radiograhy.
*/
class VTKM_RENDERING_EXPORT ConnectivityTracer
{
public:
ConnectivityTracer()
: MeshContainer(nullptr)
, CountRayStatus(false)
, UnitScalar(1.f)
{
}
@ -120,57 +116,75 @@ public:
void SetSampleDistance(const vtkm::Float32& distance);
void SetColorMap(const vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>>& colorMap);
void Trace(Ray<vtkm::Float32>& rays);
void Trace(Ray<vtkm::Float64>& rays);
MeshConnContainer* GetMeshContainer() { return MeshContainer; }
void Init();
void SetDebugOn(bool on) { CountRayStatus = on; }
void SetUnitScalar(const vtkm::Float32 unitScalar) { UnitScalar = unitScalar; }
vtkm::Id GetNumberOfMeshCells() const;
void ResetTimers();
void LogTimers();
///
/// Traces rays fully through the mesh. Rays can exit and re-enter
/// multiple times before leaving the domain. This is fast path for
/// structured meshs or meshes that are not interlocking.
/// Note: rays will be compacted
///
template <typename FloatType>
void FullTrace(Ray<FloatType>& rays);
///
/// Integrates rays through the mesh. If rays leave the mesh and
/// re-enter, then those become two separate partial composites.
/// This is need to support domain decompositions that are like
/// puzzle pieces. Note: rays will be compacted
///
template <typename FloatType>
std::vector<PartialComposite<FloatType>> PartialTrace(Ray<FloatType>& rays);
///
/// Integrates the active rays though the mesh until all rays
/// have exited.
/// Precondition: rays.HitIdx is set to a valid mesh cell
///
template <typename FloatType>
void IntegrateMeshSegment(Ray<FloatType>& rays);
///
/// Find the entry point in the mesh
///
template <typename FloatType>
void FindMeshEntry(Ray<FloatType>& rays);
private:
friend struct detail::RenderFunctor;
template <typename FloatType, typename Device>
void IntersectCell(Ray<FloatType>& rays,
detail::RayTracking<FloatType>& tracker,
const MeshConnectivityBase* meshConn,
Device);
template <typename FloatType>
void IntersectCell(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker);
template <typename FloatType, typename Device>
void AccumulatePathLengths(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker, Device);
template <typename FloatType>
void AccumulatePathLengths(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker);
template <typename FloatType, typename Device>
void FindLostRays(Ray<FloatType>& rays,
detail::RayTracking<FloatType>& tracker,
const MeshConnectivityBase* meshConn,
Device);
template <typename FloatType>
void FindLostRays(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker);
template <typename FloatType, typename Device>
void SampleCells(Ray<FloatType>& rays,
detail::RayTracking<FloatType>& tracker,
const MeshConnectivityBase* meshConn,
Device);
template <typename FloatType>
void SampleCells(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker);
template <typename FloatType, typename Device>
void IntegrateCells(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker, Device);
template <typename FloatType>
void IntegrateCells(Ray<FloatType>& rays, detail::RayTracking<FloatType>& tracker);
template <typename FloatType, typename Device>
void OffsetMinDistances(Ray<FloatType>& rays, Device);
template <typename FloatType>
void OffsetMinDistances(Ray<FloatType>& rays);
template <typename FloatType, typename Device>
void PrintRayStatus(Ray<FloatType>& rays, Device);
template <typename FloatType>
void PrintRayStatus(Ray<FloatType>& rays);
protected:
template <typename Device, typename FloatType>
void RenderOnDevice(Ray<FloatType>& rays, Device);
// Data set info
vtkm::cont::Field ScalarField;
vtkm::cont::Field EmissionField;
@ -180,7 +194,6 @@ protected:
vtkm::Float32 BoundingBox[6];
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>> ColorMap;
vtkm::cont::ArrayHandle<vtkm::Id> PreviousCellIds;
vtkm::Vec<vtkm::Float32, 4> BackgroundColor;
vtkm::Float32 SampleDistance;
@ -204,6 +217,7 @@ protected:
vtkm::Float64 SampleTime;
vtkm::Float64 LostRayTime;
vtkm::Float64 MeshEntryTime;
vtkm::Float32 UnitScalar;
}; // class ConnectivityTracer<CellType,ConnectivityType>
}

1276
vtkm/rendering/raytracing/ConnectivityTracer.hxx
View File

File diff suppressed because it is too large Load Diff

186
vtkm/rendering/raytracing/CylinderIntersector.cxx
View File

@ -34,8 +34,126 @@ namespace raytracing
namespace detail
{
class FindCylinderAABBs : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
FindCylinderAABBs() {}
typedef void ControlSignature(FieldIn<>,
FieldIn<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
WholeArrayIn<Vec3RenderingTypes>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8, _9);
template <typename PointPortalType>
VTKM_EXEC void operator()(const vtkm::Id3 cylId,
const vtkm::Float32& radius,
vtkm::Float32& xmin,
vtkm::Float32& ymin,
vtkm::Float32& zmin,
vtkm::Float32& xmax,
vtkm::Float32& ymax,
vtkm::Float32& zmax,
const PointPortalType& points) const
{
// cast to Float32
vtkm::Vec<vtkm::Float32, 3> point1, point2;
vtkm::Vec<vtkm::Float32, 3> temp;
point1 = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(cylId[1]));
point2 = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(cylId[2]));
temp[0] = radius;
temp[1] = 0.0f;
temp[2] = 0.0f;
xmin = ymin = zmin = vtkm::Infinity32();
xmax = ymax = zmax = vtkm::NegativeInfinity32();
//set first point to max and min
Bounds(point1, radius, xmin, ymin, zmin, xmax, ymax, zmax);
Bounds(point2, radius, xmin, ymin, zmin, xmax, ymax, zmax);
}
VTKM_EXEC void Bounds(const vtkm::Vec<vtkm::Float32, 3>& point,
const vtkm::Float32& radius,
vtkm::Float32& xmin,
vtkm::Float32& ymin,
vtkm::Float32& zmin,
vtkm::Float32& xmax,
vtkm::Float32& ymax,
vtkm::Float32& zmax) const
{
vtkm::Vec<vtkm::Float32, 3> temp, p;
temp[0] = radius;
temp[1] = 0.0f;
temp[2] = 0.0f;
p = point + temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
p = point - temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
temp[0] = 0.f;
temp[1] = radius;
temp[2] = 0.f;
p = point + temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
p = point - temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
temp[0] = 0.f;
temp[1] = 0.f;
temp[2] = radius;
p = point + temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
p = point - temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
}
}; //class FindAABBs
template <typename Device>
class CylinderLeafIntersector : public vtkm::cont::ExecutionObjectBase
class CylinderLeafIntersector
{
public:
using IdHandle = vtkm::cont::ArrayHandle<vtkm::Id3>;
@ -45,6 +163,8 @@ public:
IdArrayPortal CylIds;
FloatPortal Radii;
CylinderLeafIntersector() {}
CylinderLeafIntersector(const IdHandle& cylIds, const FloatHandle& radii)
: CylIds(cylIds.PrepareForInput(Device()))
, Radii(radii.PrepareForInput(Device()))
@ -183,17 +303,25 @@ public:
}
};
struct IntersectFunctor
class CylinderLeafWrapper : public vtkm::cont::ExecutionObjectBase
{
template <typename Device, typename Precision>
VTKM_CONT bool operator()(Device,
CylinderIntersector* self,
Ray<Precision>& rays,
bool returnCellIndex)
protected:
using IdHandle = vtkm::cont::ArrayHandle<vtkm::Id3>;
using FloatHandle = vtkm::cont::ArrayHandle<vtkm::Float32>;
IdHandle CylIds;
FloatHandle Radii;
public:
CylinderLeafWrapper(IdHandle& cylIds, FloatHandle radii)
: CylIds(cylIds)
, Radii(radii)
{
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
self->IntersectRaysImp(Device(), rays, returnCellIndex);
return true;
}
template <typename Device>
VTKM_CONT CylinderLeafIntersector<Device> PrepareForExecution(Device) const
{
return CylinderLeafIntersector<Device>(CylIds, Radii);
}
};
@ -302,27 +430,47 @@ CylinderIntersector::~CylinderIntersector()
{
}
void CylinderIntersector::SetData(const vtkm::cont::CoordinateSystem& coords,
vtkm::cont::ArrayHandle<vtkm::Id3> cylIds,
vtkm::cont::ArrayHandle<vtkm::Float32> radii)
{
this->Radii = radii;
this->CylIds = cylIds;
this->CoordsHandle = coords;
AABBs AABB;
vtkm::worklet::DispatcherMapField<detail::FindCylinderAABBs>(detail::FindCylinderAABBs())
.Invoke(this->CylIds,
this->Radii,
AABB.xmins,
AABB.ymins,
AABB.zmins,
AABB.xmaxs,
AABB.ymaxs,
AABB.zmaxs,
CoordsHandle);
this->SetAABBs(AABB);
}
void CylinderIntersector::IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex)
{
vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
IntersectRaysImp(rays, returnCellIndex);
}
void CylinderIntersector::IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex)
{
vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
IntersectRaysImp(rays, returnCellIndex);
}
template <typename Device, typename Precision>
void CylinderIntersector::IntersectRaysImp(Device,
Ray<Precision>& rays,
bool vtkmNotUsed(returnCellIndex))
template <typename Precision>
void CylinderIntersector::IntersectRaysImp(Ray<Precision>& rays, bool vtkmNotUsed(returnCellIndex))
{
detail::CylinderLeafIntersector<Device> leafIntersector(this->CylIds, Radii);
detail::CylinderLeafWrapper leafIntersector(this->CylIds, Radii);
BVHTraverser<detail::CylinderLeafIntersector> traverser;
traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle, Device());
BVHTraverser traverser;
traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle);
RayOperations::UpdateRayStatus(rays);
}

143
vtkm/rendering/raytracing/CylinderIntersector.h
View File

@ -31,125 +31,6 @@ namespace raytracing
{
namespace detail
{
class FindCylinderAABBs : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
FindCylinderAABBs() {}
typedef void ControlSignature(FieldIn<>,
FieldIn<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
WholeArrayIn<Vec3RenderingTypes>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8, _9);
template <typename PointPortalType>
VTKM_EXEC void operator()(const vtkm::Id3 cylId,
const vtkm::Float32& radius,
vtkm::Float32& xmin,
vtkm::Float32& ymin,
vtkm::Float32& zmin,
vtkm::Float32& xmax,
vtkm::Float32& ymax,
vtkm::Float32& zmax,
const PointPortalType& points) const
{
// cast to Float32
vtkm::Vec<vtkm::Float32, 3> point1, point2;
vtkm::Vec<vtkm::Float32, 3> temp;
point1 = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(cylId[1]));
point2 = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(cylId[2]));
temp[0] = radius;
temp[1] = 0.0f;
temp[2] = 0.0f;
xmin = ymin = zmin = vtkm::Infinity32();
xmax = ymax = zmax = vtkm::NegativeInfinity32();
//set first point to max and min
Bounds(point1, radius, xmin, ymin, zmin, xmax, ymax, zmax);
Bounds(point2, radius, xmin, ymin, zmin, xmax, ymax, zmax);
}
VTKM_EXEC void Bounds(const vtkm::Vec<vtkm::Float32, 3>& point,
const vtkm::Float32& radius,
vtkm::Float32& xmin,
vtkm::Float32& ymin,
vtkm::Float32& zmin,
vtkm::Float32& xmax,
vtkm::Float32& ymax,
vtkm::Float32& zmax) const
{
vtkm::Vec<vtkm::Float32, 3> temp, p;
temp[0] = radius;
temp[1] = 0.0f;
temp[2] = 0.0f;
p = point + temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
p = point - temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
temp[0] = 0.f;
temp[1] = radius;
temp[2] = 0.f;
p = point + temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
p = point - temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
temp[0] = 0.f;
temp[1] = 0.f;
temp[2] = radius;
p = point + temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
p = point - temp;
xmin = vtkm::Min(xmin, p[0]);
xmax = vtkm::Max(xmax, p[0]);
ymin = vtkm::Min(ymin, p[1]);
ymax = vtkm::Max(ymax, p[1]);
zmin = vtkm::Min(zmin, p[2]);
zmax = vtkm::Max(zmax, p[2]);
}
}; //class FindAABBs
}
class CylinderIntersector : public ShapeIntersector
{
@ -164,33 +45,15 @@ public:
void SetData(const vtkm::cont::CoordinateSystem& coords,
vtkm::cont::ArrayHandle<vtkm::Id3> cylIds,
vtkm::cont::ArrayHandle<vtkm::Float32> radii)
{
this->Radii = radii;
this->CylIds = cylIds;
this->CoordsHandle = coords;
AABBs AABB;
vtkm::cont::ArrayHandle<vtkm::Float32> radii);
vtkm::worklet::DispatcherMapField<detail::FindCylinderAABBs>(detail::FindCylinderAABBs())
.Invoke(this->CylIds,
this->Radii,
AABB.xmins,
AABB.ymins,
AABB.zmins,
AABB.xmaxs,
AABB.ymaxs,
AABB.zmaxs,
CoordsHandle);
this->SetAABBs(AABB);
}
void IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex = false) override;
void IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex = false) override;
template <typename Device, typename Precision>
void IntersectRaysImp(Device, Ray<Precision>& rays, bool returnCellIndex);
template <typename Precision>
void IntersectRaysImp(Ray<Precision>& rays, bool returnCellIndex);
template <typename Precision>

34
vtkm/rendering/raytracing/MeshConnectivityBase.h
View File

@ -37,7 +37,10 @@ namespace rendering
{
namespace raytracing
{
//
// Base class for different types of face-to-connecting-cell
// and other mesh information
//
class VTKM_ALWAYS_EXPORT MeshConnectivityBase : public VirtualObjectBase
{
public:
@ -51,6 +54,35 @@ public:
virtual vtkm::UInt8 GetCellShape(const vtkm::Id& cellId) const = 0;
};
// A simple concrete type to wrap MeshConnectivityBase so we can
// pass an ExeObject to worklets.
class MeshWrapper
{
private:
MeshConnectivityBase* MeshConn;
public:
MeshWrapper() {}
MeshWrapper(MeshConnectivityBase* meshConn)
: MeshConn(meshConn){};
VTKM_EXEC_CONT
vtkm::Id GetConnectingCell(const vtkm::Id& cellId, const vtkm::Id& face) const
{
return MeshConn->GetConnectingCell(cellId, face);
}
VTKM_EXEC_CONT
vtkm::Int32 GetCellIndices(vtkm::Id cellIndices[8], const vtkm::Id& cellId) const
{
return MeshConn->GetCellIndices(cellIndices, cellId);
}
VTKM_EXEC_CONT
vtkm::UInt8 GetCellShape(const vtkm::Id& cellId) const { return MeshConn->GetCellShape(cellId); }
};
class VTKM_ALWAYS_EXPORT MeshConnStructured : public MeshConnectivityBase
{
protected:

81
vtkm/rendering/raytracing/MeshConnectivityBuilder.cxx
View File

@ -40,10 +40,10 @@ namespace rendering
namespace raytracing
{
struct IsExternal
struct IsUnique
{
VTKM_EXEC_CONT
inline bool operator()(const vtkm::Id& x) const { return (x < 0); }
inline bool operator()(const vtkm::Int32& x) const { return (x < 0); }
}; //struct IsExternal
class CountFaces : public vtkm::worklet::WorkletMapField
@ -93,8 +93,9 @@ public:
WholeArrayIn<>,
WholeArrayIn<>,
WholeArrayIn<>,
WholeArrayOut<>);
using ExecutionSignature = void(_1, _2, WorkIndex, _3, _4, _5, _6);
WholeArrayOut<>,
WholeArrayInOut<>);
using ExecutionSignature = void(_1, _2, WorkIndex, _3, _4, _5, _6, _7);
VTKM_EXEC
inline vtkm::Int32 GetShapeOffset(const vtkm::UInt8& shapeType) const
@ -145,14 +146,16 @@ public:
typename ConnPortalType,
typename ShapePortalType,
typename OffsetPortalType,
typename ExternalFaceFlagType>
typename ExternalFaceFlagType,
typename UniqueFacesType>
VTKM_EXEC inline void operator()(const MortonPortalType& mortonCodes,
FaceIdPairsPortalType& faceIdPairs,
const vtkm::Id& index,
const ConnPortalType& connectivity,
const ShapePortalType& shapes,
const OffsetPortalType& offsets,
ExternalFaceFlagType& flags) const
ExternalFaceFlagType& flags,
UniqueFacesType& uniqueFaces) const
{
if (index == 0)
{
@ -252,6 +255,10 @@ public:
flags.Set(currentIndex, myCell);
BOUNDS_CHECK(flags, index);
flags.Set(index, connectedCell);
// for unstructured, we want all unique faces to intersect with
// just choose one and mark as unique so the other gets culled.
uniqueFaces.Set(index, 1);
}
}
}; //class Neighbor
@ -522,7 +529,8 @@ VTKM_CONT void GenerateFaceConnnectivity(
vtkm::cont::ArrayHandle<vtkm::Id>& faceConnectivity,
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 3>>& cellFaceId,
vtkm::Float32 BoundingBox[6],
vtkm::cont::ArrayHandle<vtkm::Id>& faceOffsets)
vtkm::cont::ArrayHandle<vtkm::Id>& faceOffsets,
vtkm::cont::ArrayHandle<vtkm::Int32>& uniqueFaces)
{
vtkm::cont::Timer<vtkm::cont::DeviceAdapterTagSerial> timer;
@ -578,8 +586,8 @@ VTKM_CONT void GenerateFaceConnnectivity(
vtkm::cont::ArrayHandle<vtkm::UInt32> faceMortonCodes;
cellFaceId.Allocate(totalFaces);
faceMortonCodes.Allocate(totalFaces);
//cellFaceId.PrepareForOutput(totalFaces, DeviceAdapter());
//faceMortonCodes.PrepareForOutput(totalFaces, DeviceAdapter());
uniqueFaces.Allocate(totalFaces);
vtkm::worklet::DispatcherMapTopology<MortonCodeFace>(MortonCodeFace(inverseExtent, minPoint))
.Invoke(cellSet, coordinates, cellOffsets, faceMortonCodes, cellFaceId);
@ -594,19 +602,21 @@ VTKM_CONT void GenerateFaceConnnectivity(
vtkm::cont::ArrayHandleConstant<vtkm::Id> negOne(-1, totalFaces);
vtkm::cont::Algorithm::Copy(negOne, faceConnectivity);
vtkm::cont::ArrayHandleConstant<vtkm::Int32> negOne32(-1, totalFaces);
vtkm::cont::Algorithm::Copy(negOne32, uniqueFaces);
vtkm::worklet::DispatcherMapField<MortonNeighbor>(MortonNeighbor())
.Invoke(faceMortonCodes, cellFaceId, conn, shapes, shapeOffsets, faceConnectivity);
.Invoke(faceMortonCodes, cellFaceId, conn, shapes, shapeOffsets, faceConnectivity, uniqueFaces);
vtkm::Float64 time = timer.GetElapsedTime();
Logger::GetInstance()->AddLogData("gen_face_conn", time);
}
template <typename ShapeHandleType, typename OffsetsHandleType, typename ConnHandleType>
VTKM_CONT vtkm::cont::ArrayHandle<vtkm::Vec<Id, 4>> ExtractExternalFaces(
VTKM_CONT vtkm::cont::ArrayHandle<vtkm::Vec<Id, 4>> ExtractFaces(
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 3>>
cellFaceId, // Map of cell, face, and connecting cell
vtkm::cont::ArrayHandle<vtkm::Id>
faceConnectivity, // -1 if the face does not connect to any other face
cellFaceId, // Map of cell, face, and connecting cell
vtkm::cont::ArrayHandle<vtkm::Int32> uniqueFaces, // -1 if the face is unique
const ShapeHandleType& shapes,
const ConnHandleType& conn,
const OffsetsHandleType& shapeOffsets)
@ -614,7 +624,7 @@ VTKM_CONT vtkm::cont::ArrayHandle<vtkm::Vec<Id, 4>> ExtractExternalFaces(
vtkm::cont::Timer<vtkm::cont::DeviceAdapterTagSerial> timer;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 3>> externalFacePairs;
vtkm::cont::Algorithm::CopyIf(cellFaceId, faceConnectivity, externalFacePairs, IsExternal());
vtkm::cont::Algorithm::CopyIf(cellFaceId, uniqueFaces, externalFacePairs, IsUnique());
// We need to count the number of triangle per external face
// If it is a single cell type and it is a tet or hex, this is a special case
@ -690,6 +700,7 @@ void MeshConnectivityBuilder::BuildConnectivity(
vtkm::cont::ArrayHandle<vtkm::Id> faceConnectivity;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 3>> cellFaceId;
vtkm::cont::ArrayHandle<vtkm::Int32> uniqueFaces;
GenerateFaceConnnectivity(cellSetUnstructured,
shapes,
@ -699,13 +710,13 @@ void MeshConnectivityBuilder::BuildConnectivity(
faceConnectivity,
cellFaceId,
BoundingBox,
FaceOffsets);
FaceOffsets,
uniqueFaces);
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> externalTriangles;
//External Faces
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> triangles;
//Faces
externalTriangles =
ExtractExternalFaces(cellFaceId, faceConnectivity, shapes, conn, shapeOffsets);
triangles = ExtractFaces(cellFaceId, uniqueFaces, shapes, conn, shapeOffsets);
// scatter the coonectivity into the original order
@ -713,7 +724,7 @@ void MeshConnectivityBuilder::BuildConnectivity(
.Invoke(cellFaceId, this->FaceOffsets, faceConnectivity);
FaceConnectivity = faceConnectivity;
OutsideTriangles = externalTriangles;
Triangles = triangles;
vtkm::Float64 time = timer.GetElapsedTime();
logger->CloseLogEntry(time);
@ -748,6 +759,7 @@ void MeshConnectivityBuilder::BuildConnectivity(
vtkm::cont::ArrayHandle<vtkm::Id> faceConnectivity;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 3>> cellFaceId;
vtkm::cont::ArrayHandle<vtkm::Int32> uniqueFaces;
GenerateFaceConnnectivity(cellSetUnstructured,
shapes,
@ -757,20 +769,20 @@ void MeshConnectivityBuilder::BuildConnectivity(
faceConnectivity,
cellFaceId,
BoundingBox,
FaceOffsets);
FaceOffsets,
uniqueFaces);
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> externalTriangles;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> triangles;
//
//External Faces
externalTriangles =
ExtractExternalFaces(cellFaceId, faceConnectivity, shapes, conn, shapeOffsets);
//Faces
triangles = ExtractFaces(cellFaceId, uniqueFaces, shapes, conn, shapeOffsets);
// scatter the coonectivity into the original order
vtkm::worklet::DispatcherMapField<WriteFaceConn>(WriteFaceConn())
.Invoke(cellFaceId, this->FaceOffsets, faceConnectivity);
FaceConnectivity = faceConnectivity;
OutsideTriangles = externalTriangles;
Triangles = triangles;
vtkm::Float64 time = timer.GetElapsedTime();
logger->CloseLogEntry(time);
@ -829,9 +841,9 @@ vtkm::cont::ArrayHandle<vtkm::Id> MeshConnectivityBuilder::GetFaceOffsets()
return FaceOffsets;
}
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> MeshConnectivityBuilder::GetExternalTriangles()
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> MeshConnectivityBuilder::GetTriangles()
{
return OutsideTriangles;
return Triangles;
}
VTKM_CONT
@ -894,21 +906,20 @@ MeshConnContainer* MeshConnectivityBuilder::BuildConnectivity(
{
vtkm::cont::CellSetExplicit<> cells = cellset.Cast<vtkm::cont::CellSetExplicit<>>();
this->BuildConnectivity(cells, coordinates.GetData(), coordBounds);
meshConn = new UnstructuredContainer(
cells, coordinates, FaceConnectivity, FaceOffsets, OutsideTriangles);
meshConn =
new UnstructuredContainer(cells, coordinates, FaceConnectivity, FaceOffsets, Triangles);
}
else if (type == UnstructuredSingle)
{
vtkm::cont::CellSetSingleType<> cells = cellset.Cast<vtkm::cont::CellSetSingleType<>>();
this->BuildConnectivity(cells, coordinates.GetData(), coordBounds);
meshConn =
new UnstructuredSingleContainer(cells, coordinates, FaceConnectivity, OutsideTriangles);
meshConn = new UnstructuredSingleContainer(cells, coordinates, FaceConnectivity, Triangles);
}
else if (type == Structured)
{
vtkm::cont::CellSetStructured<3> cells = cellset.Cast<vtkm::cont::CellSetStructured<3>>();
OutsideTriangles = this->ExternalTrianglesStructured(cells);
meshConn = new StructuredContainer(cells, coordinates, OutsideTriangles);
Triangles = this->ExternalTrianglesStructured(cells);
meshConn = new StructuredContainer(cells, coordinates, Triangles);
}
else
{

4
vtkm/rendering/raytracing/MeshConnectivityBuilder.h
View File

@ -48,7 +48,7 @@ public:
vtkm::cont::ArrayHandle<vtkm::Id> GetFaceOffsets();
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> GetExternalTriangles();
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> GetTriangles();
protected:
VTKM_CONT
@ -63,7 +63,7 @@ protected:
vtkm::cont::ArrayHandle<vtkm::Id> FaceConnectivity;
vtkm::cont::ArrayHandle<vtkm::Id> FaceOffsets;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> OutsideTriangles;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> Triangles;
};
}
}

84
vtkm/rendering/raytracing/MeshConnectivityContainers.cxx
View File

@ -39,41 +39,28 @@ MeshConnContainer::MeshConnContainer(){};
MeshConnContainer::~MeshConnContainer(){};
template <typename T>
VTKM_CONT void MeshConnContainer::FindEntryImpl(Ray<T>& rays,
const vtkm::cont::DeviceAdapterId deviceId)
VTKM_CONT void MeshConnContainer::FindEntryImpl(Ray<T>& rays)
{
bool getCellIndex = true;
Intersector.SetUseWaterTight(true);
switch (deviceId.GetValue())
{
#ifdef VTKM_ENABLE_TBB
case VTKM_DEVICE_ADAPTER_TBB:
Intersector.IntersectRays(rays, vtkm::cont::DeviceAdapterTagTBB(), getCellIndex);
break;
#endif
#ifdef VTKM_ENABLE_CUDA
case VTKM_DEVICE_ADAPTER_CUDA:
Intersector.IntersectRays(rays, vtkm::cont::DeviceAdapterTagCuda(), getCellIndex);
break;
#endif
default:
Intersector.IntersectRays(rays, vtkm::cont::DeviceAdapterTagSerial(), getCellIndex);
break;
}
Intersector.IntersectRays(rays, getCellIndex);
}
void MeshConnContainer::FindEntry(Ray<vtkm::Float32>& rays,
const vtkm::cont::DeviceAdapterId deviceId)
MeshWrapper MeshConnContainer::PrepareForExecution(const vtkm::cont::DeviceAdapterId deviceId)
{
this->FindEntryImpl(rays, deviceId);
return MeshWrapper(const_cast<MeshConnectivityBase*>(this->Construct(deviceId)));
}
void MeshConnContainer::FindEntry(Ray<vtkm::Float64>& rays,
const vtkm::cont::DeviceAdapterId deviceId)
void MeshConnContainer::FindEntry(Ray<vtkm::Float32>& rays)
{
this->FindEntryImpl(rays, deviceId);
this->FindEntryImpl(rays);
}
void MeshConnContainer::FindEntry(Ray<vtkm::Float64>& rays)
{
this->FindEntryImpl(rays);
}
VTKM_CONT
@ -81,13 +68,13 @@ UnstructuredContainer::UnstructuredContainer(const vtkm::cont::CellSetExplicit<>
const vtkm::cont::CoordinateSystem& coords,
IdHandle& faceConn,
IdHandle& faceOffsets,
Id4Handle& externalTriangles)
Id4Handle& triangles)
: FaceConnectivity(faceConn)
, FaceOffsets(faceOffsets)
, Cellset(cellset)
, Coords(coords)
{
this->ExternalTriangles = externalTriangles;
this->Triangles = triangles;
//
// Grab the cell arrays
//
@ -97,7 +84,7 @@ UnstructuredContainer::UnstructuredContainer(const vtkm::cont::CellSetExplicit<>
Cellset.GetIndexOffsetArray(vtkm::TopologyElementTagPoint(), vtkm::TopologyElementTagCell());
Shapes = Cellset.GetShapesArray(vtkm::TopologyElementTagPoint(), vtkm::TopologyElementTagCell());
Intersector.SetData(Coords, ExternalTriangles);
Intersector.SetData(Coords, Triangles);
}
UnstructuredContainer::~UnstructuredContainer(){};
@ -107,6 +94,19 @@ const MeshConnectivityBase* UnstructuredContainer::Construct(
{
switch (deviceId.GetValue())
{
#ifdef VTKM_ENABLE_OPENMP
case VTKM_DEVICE_ADAPTER_OPENMP:
using OMP = vtkm::cont::DeviceAdapterTagOpenMP;
{
MeshConnUnstructured<OMP> conn(this->FaceConnectivity,
this->FaceOffsets,
this->CellConn,
this->CellOffsets,
this->Shapes);
Handle = make_MeshConnHandle(conn);
}
return Handle.PrepareForExecution(OMP());
#endif
#ifdef VTKM_ENABLE_TBB
case VTKM_DEVICE_ADAPTER_TBB:
using TBB = vtkm::cont::DeviceAdapterTagTBB;
@ -157,13 +157,13 @@ UnstructuredSingleContainer::UnstructuredSingleContainer(
const vtkm::cont::CellSetSingleType<>& cellset,
const vtkm::cont::CoordinateSystem& coords,
IdHandle& faceConn,
Id4Handle& externalTriangles)
Id4Handle& triangles)
: FaceConnectivity(faceConn)
, Coords(coords)
, Cellset(cellset)
{
this->ExternalTriangles = externalTriangles;
this->Triangles = triangles;
this->Intersector.SetUseWaterTight(true);
@ -192,7 +192,7 @@ UnstructuredSingleContainer::UnstructuredSingleContainer(
logger->OpenLogEntry("mesh_conn_construction");
vtkm::cont::Timer<cont::DeviceAdapterTagSerial> timer;
Intersector.SetData(Coords, ExternalTriangles);
Intersector.SetData(Coords, Triangles);
}
const MeshConnectivityBase* UnstructuredSingleContainer::Construct(
@ -200,6 +200,20 @@ const MeshConnectivityBase* UnstructuredSingleContainer::Construct(
{
switch (deviceId.GetValue())
{
#ifdef VTKM_ENABLE_OPENMP
case VTKM_DEVICE_ADAPTER_OPENMP:
using OMP = vtkm::cont::DeviceAdapterTagOpenMP;
{
MeshConnSingleType<OMP> conn(this->FaceConnectivity,
this->CellConnectivity,
this->CellOffsets,
this->ShapeId,
this->NumIndices,
this->NumFaces);
Handle = make_MeshConnHandle(conn);
}
return Handle.PrepareForExecution(OMP());
#endif
#ifdef VTKM_ENABLE_TBB
case VTKM_DEVICE_ADAPTER_TBB:
using TBB = vtkm::cont::DeviceAdapterTagTBB;
@ -245,18 +259,18 @@ const MeshConnectivityBase* UnstructuredSingleContainer::Construct(
StructuredContainer::StructuredContainer(const vtkm::cont::CellSetStructured<3>& cellset,
const vtkm::cont::CoordinateSystem& coords,
Id4Handle& externalTriangles)
Id4Handle& triangles)
: Coords(coords)
, Cellset(cellset)
{
ExternalTriangles = externalTriangles;
Triangles = triangles;
Intersector.SetUseWaterTight(true);
PointDims = Cellset.GetPointDimensions();
CellDims = Cellset.GetCellDimensions();
this->Intersector.SetData(Coords, ExternalTriangles);
this->Intersector.SetData(Coords, Triangles);
}
const MeshConnectivityBase* StructuredContainer::Construct(
@ -268,6 +282,10 @@ const MeshConnectivityBase* StructuredContainer::Construct(
switch (deviceId.GetValue())
{
#ifdef VTKM_ENABLE_OPENMP
case VTKM_DEVICE_ADAPTER_OPENMP:
return Handle.PrepareForExecution(vtkm::cont::DeviceAdapterTagOpenMP());
#endif
#ifdef VTKM_ENABLE_TBB
case VTKM_DEVICE_ADAPTER_TBB:
return Handle.PrepareForExecution(vtkm::cont::DeviceAdapterTagTBB());

16
vtkm/rendering/raytracing/MeshConnectivityContainers.h
View File

@ -31,7 +31,7 @@ namespace rendering
namespace raytracing
{
class MeshConnContainer
class MeshConnContainer : vtkm::cont::ExecutionObjectBase
{
public:
MeshConnContainer();
@ -39,17 +39,19 @@ public:
virtual const MeshConnectivityBase* Construct(const vtkm::cont::DeviceAdapterId deviceId) = 0;
MeshWrapper PrepareForExecution(const vtkm::cont::DeviceAdapterId deviceId);
template <typename T>
VTKM_CONT void FindEntryImpl(Ray<T>& rays, const vtkm::cont::DeviceAdapterId deviceId);
VTKM_CONT void FindEntryImpl(Ray<T>& rays);
void FindEntry(Ray<vtkm::Float32>& rays, const vtkm::cont::DeviceAdapterId deviceId);
void FindEntry(Ray<vtkm::Float32>& rays);
void FindEntry(Ray<vtkm::Float64>& rays, const vtkm::cont::DeviceAdapterId deviceId);
void FindEntry(Ray<vtkm::Float64>& rays);
protected:
using Id4Handle = typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>>;
// Mesh Boundary
Id4Handle ExternalTriangles;
Id4Handle Triangles;
TriangleIntersector Intersector;
MeshConnHandle Handle;
};
@ -83,7 +85,7 @@ public:
const vtkm::cont::CoordinateSystem& coords,
IdHandle& faceConn,
IdHandle& faceOffsets,
Id4Handle& externalTriangles);
Id4Handle& triangles);
virtual ~UnstructuredContainer();
@ -108,7 +110,7 @@ public:
VTKM_CONT
StructuredContainer(const vtkm::cont::CellSetStructured<3>& cellset,
const vtkm::cont::CoordinateSystem& coords,
Id4Handle& externalTriangles);
Id4Handle& triangles);
const MeshConnectivityBase* Construct(const vtkm::cont::DeviceAdapterId deviceId) override;

141
vtkm/rendering/raytracing/QuadIntersector.cxx
View File

@ -34,17 +34,91 @@ namespace raytracing
namespace detail
{
#define QUAD_AABB_EPSILON 1.0e-4f
class FindQuadAABBs : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
FindQuadAABBs() {}
typedef void ControlSignature(FieldIn<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
WholeArrayIn<Vec3RenderingTypes>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8);
template <typename PointPortalType>
VTKM_EXEC void operator()(const vtkm::Vec<vtkm::Id, 5> quadId,
vtkm::Float32& xmin,
vtkm::Float32& ymin,
vtkm::Float32& zmin,
vtkm::Float32& xmax,
vtkm::Float32& ymax,
vtkm::Float32& zmax,
const PointPortalType& points) const
{
// cast to Float32
vtkm::Vec<vtkm::Float32, 3> q, r, s, t;
q = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[1]));
r = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[2]));
s = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[3]));
t = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[4]));
xmin = q[0];
ymin = q[1];
zmin = q[2];
xmax = xmin;
ymax = ymin;
zmax = zmin;
xmin = vtkm::Min(xmin, r[0]);
ymin = vtkm::Min(ymin, r[1]);
zmin = vtkm::Min(zmin, r[2]);
xmax = vtkm::Max(xmax, r[0]);
ymax = vtkm::Max(ymax, r[1]);
zmax = vtkm::Max(zmax, r[2]);
xmin = vtkm::Min(xmin, s[0]);
ymin = vtkm::Min(ymin, s[1]);
zmin = vtkm::Min(zmin, s[2]);
xmax = vtkm::Max(xmax, s[0]);
ymax = vtkm::Max(ymax, s[1]);
zmax = vtkm::Max(zmax, s[2]);
xmin = vtkm::Min(xmin, t[0]);
ymin = vtkm::Min(ymin, t[1]);
zmin = vtkm::Min(zmin, t[2]);
xmax = vtkm::Max(xmax, t[0]);
ymax = vtkm::Max(ymax, t[1]);
zmax = vtkm::Max(zmax, t[2]);
vtkm::Float32 xEpsilon, yEpsilon, zEpsilon;
const vtkm::Float32 minEpsilon = 1e-6f;
xEpsilon = vtkm::Max(minEpsilon, QUAD_AABB_EPSILON * (xmax - xmin));
yEpsilon = vtkm::Max(minEpsilon, QUAD_AABB_EPSILON * (ymax - ymin));
zEpsilon = vtkm::Max(minEpsilon, QUAD_AABB_EPSILON * (zmax - zmin));
xmin -= xEpsilon;
ymin -= yEpsilon;
zmin -= zEpsilon;
xmax += xEpsilon;
ymax += yEpsilon;
zmax += zEpsilon;
}
}; //class FindAABBs
template <typename Device>
class QuadLeafIntersector : public vtkm::cont::ExecutionObjectBase
class QuadLeafIntersector
{
public:
using IdType = vtkm::Vec<vtkm::Id, 5>;
using IdHandle = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 5>>;
using IdArrayPortal = typename IdHandle::ExecutionTypes<Device>::PortalConst;
using FloatHandle = vtkm::cont::ArrayHandle<vtkm::Float32>;
using FloatPortal = typename FloatHandle::ExecutionTypes<Device>::PortalConst;
IdArrayPortal QuadIds;
QuadLeafIntersector() {}
QuadLeafIntersector(const IdHandle& quadIds)
: QuadIds(quadIds.PrepareForInput(Device()))
{
@ -228,17 +302,23 @@ public:
}
};
struct IntersectFunctor
class QuadExecWrapper : public vtkm::cont::ExecutionObjectBase
{
template <typename Device, typename Precision>
VTKM_CONT bool operator()(Device,
QuadIntersector* self,
Ray<Precision>& rays,
bool returnCellIndex)
protected:
using IdType = vtkm::Vec<vtkm::Id, 5>;
using IdHandle = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 5>>;
IdHandle QuadIds;
public:
QuadExecWrapper(IdHandle& quadIds)
: QuadIds(quadIds)
{
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
self->IntersectRaysImp(Device(), rays, returnCellIndex);
return true;
}
template <typename Device>
VTKM_CONT QuadLeafIntersector<Device> PrepareForExecution(Device) const
{
return QuadLeafIntersector<Device>(QuadIds);
}
};
@ -343,24 +423,22 @@ QuadIntersector::~QuadIntersector()
void QuadIntersector::IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex)
{
vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
IntersectRaysImp(rays, returnCellIndex);
}
void QuadIntersector::IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex)
{
vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
IntersectRaysImp(rays, returnCellIndex);
}
template <typename Device, typename Precision>
void QuadIntersector::IntersectRaysImp(Device,
Ray<Precision>& rays,
bool vtkmNotUsed(returnCellIndex))
template <typename Precision>
void QuadIntersector::IntersectRaysImp(Ray<Precision>& rays, bool vtkmNotUsed(returnCellIndex))
{
detail::QuadLeafIntersector<Device> leafIntersector(this->QuadIds);
detail::QuadExecWrapper leafIntersector(this->QuadIds);
BVHTraverser<detail::QuadLeafIntersector> traverser;
traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle, Device());
BVHTraverser traverser;
traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle);
RayOperations::UpdateRayStatus(rays);
}
@ -401,6 +479,27 @@ void QuadIntersector::IntersectionData(Ray<vtkm::Float64>& rays,
IntersectionDataImp(rays, scalarField, scalarRange);
}
void QuadIntersector::SetData(const vtkm::cont::CoordinateSystem& coords,
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 5>> quadIds)
{
this->QuadIds = quadIds;
this->CoordsHandle = coords;
AABBs AABB;
vtkm::worklet::DispatcherMapField<detail::FindQuadAABBs> faabbsInvoker;
faabbsInvoker.Invoke(this->QuadIds,
AABB.xmins,
AABB.ymins,
AABB.zmins,
AABB.xmaxs,
AABB.ymaxs,
AABB.zmaxs,
CoordsHandle);
this->SetAABBs(AABB);
}
vtkm::Id QuadIntersector::GetNumberOfShapes() const
{
return QuadIds.GetNumberOfValues();

107
vtkm/rendering/raytracing/QuadIntersector.h
View File

@ -31,80 +31,6 @@ namespace raytracing
{
namespace detail
{
#define AABB_EPSILON 1.0e-4f
class FindQuadAABBs : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
FindQuadAABBs() {}
typedef void ControlSignature(FieldIn<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
WholeArrayIn<Vec3RenderingTypes>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8);
template <typename PointPortalType>
VTKM_EXEC void operator()(const vtkm::Vec<vtkm::Id, 5> quadId,
vtkm::Float32& xmin,
vtkm::Float32& ymin,
vtkm::Float32& zmin,
vtkm::Float32& xmax,
vtkm::Float32& ymax,
vtkm::Float32& zmax,
const PointPortalType& points) const
{
// cast to Float32
vtkm::Vec<vtkm::Float32, 3> q, r, s, t;
q = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[1]));
r = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[2]));
s = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[3]));
t = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(quadId[4]));
xmin = q[0];
ymin = q[1];
zmin = q[2];
xmax = xmin;
ymax = ymin;
zmax = zmin;
xmin = vtkm::Min(xmin, r[0]);
ymin = vtkm::Min(ymin, r[1]);
zmin = vtkm::Min(zmin, r[2]);
xmax = vtkm::Max(xmax, r[0]);
ymax = vtkm::Max(ymax, r[1]);
zmax = vtkm::Max(zmax, r[2]);
xmin = vtkm::Min(xmin, s[0]);
ymin = vtkm::Min(ymin, s[1]);
zmin = vtkm::Min(zmin, s[2]);
xmax = vtkm::Max(xmax, s[0]);
ymax = vtkm::Max(ymax, s[1]);
zmax = vtkm::Max(zmax, s[2]);
xmin = vtkm::Min(xmin, t[0]);
ymin = vtkm::Min(ymin, t[1]);
zmin = vtkm::Min(zmin, t[2]);
xmax = vtkm::Max(xmax, t[0]);
ymax = vtkm::Max(ymax, t[1]);
zmax = vtkm::Max(zmax, t[2]);
vtkm::Float32 xEpsilon, yEpsilon, zEpsilon;
const vtkm::Float32 minEpsilon = 1e-6f;
xEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (xmax - xmin));
yEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (ymax - ymin));
zEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (zmax - zmin));
xmin -= xEpsilon;
ymin -= yEpsilon;
zmin -= zEpsilon;
xmax += xEpsilon;
ymax += yEpsilon;
zmax += zEpsilon;
}
}; //class FindAABBs
}
class QuadIntersector : public ShapeIntersector
{
@ -117,42 +43,15 @@ public:
void SetData(const vtkm::cont::CoordinateSystem& coords,
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 5>> quadIds)
{
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 5>> quadIds);
this->QuadIds = quadIds;
this->CoordsHandle = coords;
AABBs AABB;
vtkm::worklet::DispatcherMapField<detail::FindQuadAABBs> faabbsInvoker;
faabbsInvoker.Invoke(this->QuadIds,
AABB.xmins,
AABB.ymins,
AABB.zmins,
AABB.xmaxs,
AABB.ymaxs,
AABB.zmaxs,
CoordsHandle);
//vtkm::worklet::DispatcherMapField<detail::FindQuadAABBs> faabbsInvoker(detail::FindQuadAABBs())
// .Invoke(this->QuadIds,
// AABB.xmins,
// AABB.ymins,
// AABB.zmins,
// AABB.xmaxs,
// AABB.ymaxs,
// AABB.zmaxs,
// CoordsHandle);
this->SetAABBs(AABB);
}
void IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex = false) override;
void IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex = false) override;
template <typename Device, typename Precision>
void IntersectRaysImp(Device, Ray<Precision>& rays, bool returnCellIndex);
template <typename Precision>
void IntersectRaysImp(Ray<Precision>& rays, bool returnCellIndex);
template <typename Precision>

47
vtkm/rendering/raytracing/RayOperations.h
View File

@ -113,11 +113,11 @@ public:
class RayOperations
{
public:
template <typename Device, typename T>
static void ResetStatus(Ray<T>& rays, vtkm::UInt8 status, Device)
template <typename T>
static void ResetStatus(Ray<T>& rays, vtkm::UInt8 status)
{
vtkm::cont::ArrayHandleConstant<vtkm::UInt8> statusHandle(status, rays.NumRays);
vtkm::cont::Algorithm::Copy(Device(), statusHandle, rays.Status);
vtkm::cont::Algorithm::Copy(statusHandle, rays.Status);
}
//
@ -146,8 +146,8 @@ public:
const vtkm::rendering::Camera& camera,
const vtkm::rendering::CanvasRayTracer& canvas);
template <typename Device, typename T>
static vtkm::Id RaysInMesh(Ray<T>& rays, Device)
template <typename T>
static vtkm::Id RaysInMesh(Ray<T>& rays)
{
vtkm::Vec<UInt8, 2> maskValues;
maskValues[0] = RAY_ACTIVE;
@ -157,17 +157,16 @@ public:
vtkm::worklet::DispatcherMapField<ManyMask<vtkm::UInt8, 2>> dispatcher{ (
ManyMask<vtkm::UInt8, 2>{ maskValues }) };
dispatcher.SetDevice(Device());
dispatcher.Invoke(rays.Status, masks);
vtkm::cont::ArrayHandleCast<vtkm::Id, vtkm::cont::ArrayHandle<vtkm::UInt8>> castedMasks(masks);
const vtkm::Id initVal = 0;
vtkm::Id count = vtkm::cont::DeviceAdapterAlgorithm<Device>::Reduce(castedMasks, initVal);
vtkm::Id count = vtkm::cont::Algorithm::Reduce(castedMasks, initVal);
return count;
}
template <typename Device, typename T>
static vtkm::Id GetStatusCount(Ray<T>& rays, vtkm::Id status, Device)
template <typename T>
static vtkm::Id GetStatusCount(Ray<T>& rays, vtkm::Id status)
{
vtkm::UInt8 statusUInt8;
if (status < 0 || status > 255)
@ -180,17 +179,16 @@ public:
vtkm::worklet::DispatcherMapField<Mask<vtkm::UInt8>> dispatcher{ (
Mask<vtkm::UInt8>{ statusUInt8 }) };
dispatcher.SetDevice(Device());
dispatcher.Invoke(rays.Status, masks);
vtkm::cont::ArrayHandleCast<vtkm::Id, vtkm::cont::ArrayHandle<vtkm::UInt8>> castedMasks(masks);
const vtkm::Id initVal = 0;
vtkm::Id count = vtkm::cont::DeviceAdapterAlgorithm<Device>::Reduce(castedMasks, initVal);
vtkm::Id count = vtkm::cont::Algorithm::Reduce(castedMasks, initVal);
return count;
}
template <typename Device, typename T>
static vtkm::Id RaysProcessed(Ray<T>& rays, Device)
template <typename T>
static vtkm::Id RaysProcessed(Ray<T>& rays)
{
vtkm::Vec<UInt8, 3> maskValues;
maskValues[0] = RAY_TERMINATED;
@ -201,17 +199,16 @@ public:
vtkm::worklet::DispatcherMapField<ManyMask<vtkm::UInt8, 3>> dispatcher{ (
ManyMask<vtkm::UInt8, 3>{ maskValues }) };
dispatcher.SetDevice(Device());
dispatcher.Invoke(rays.Status, masks);
vtkm::cont::ArrayHandleCast<vtkm::Id, vtkm::cont::ArrayHandle<vtkm::UInt8>> castedMasks(masks);
const vtkm::Id initVal = 0;
vtkm::Id count = vtkm::cont::DeviceAdapterAlgorithm<Device>::Reduce(castedMasks, initVal);
vtkm::Id count = vtkm::cont::Algorithm::Reduce(castedMasks, initVal);
return count;
}
template <typename Device, typename T>
static vtkm::cont::ArrayHandle<vtkm::UInt8> CompactActiveRays(Ray<T>& rays, Device)
template <typename T>
static vtkm::cont::ArrayHandle<vtkm::UInt8> CompactActiveRays(Ray<T>& rays)
{
vtkm::Vec<UInt8, 1> maskValues;
maskValues[0] = RAY_ACTIVE;
@ -220,7 +217,6 @@ public:
vtkm::worklet::DispatcherMapField<Mask<vtkm::UInt8>> dispatcher{ (
Mask<vtkm::UInt8>{ statusUInt8 }) };
dispatcher.SetDevice(Device());
dispatcher.Invoke(rays.Status, masks);
vtkm::cont::ArrayHandle<T> emptyHandle;
@ -261,7 +257,7 @@ public:
break;
}
vtkm::cont::ArrayHandle<T> compacted;
vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(*floatArrayPointers[i], masks, compacted);
vtkm::cont::Algorithm::CopyIf(*floatArrayPointers[i], masks, compacted);
*floatArrayPointers[i] = compacted;
}
@ -275,15 +271,15 @@ public:
rays.Dir = vtkm::cont::make_ArrayHandleCompositeVector(rays.DirX, rays.DirY, rays.DirZ);
vtkm::cont::ArrayHandle<vtkm::Id> compactedHits;
vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(rays.HitIdx, masks, compactedHits);
vtkm::cont::Algorithm::CopyIf(rays.HitIdx, masks, compactedHits);
rays.HitIdx = compactedHits;
vtkm::cont::ArrayHandle<vtkm::Id> compactedPixels;
vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(rays.PixelIdx, masks, compactedPixels);
vtkm::cont::Algorithm::CopyIf(rays.PixelIdx, masks, compactedPixels);
rays.PixelIdx = compactedPixels;
vtkm::cont::ArrayHandle<vtkm::UInt8> compactedStatus;
vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(rays.Status, masks, compactedStatus);
vtkm::cont::Algorithm::CopyIf(rays.Status, masks, compactedStatus);
rays.Status = compactedStatus;
rays.NumRays = rays.Status.GetPortalConstControl().GetNumberOfValues();
@ -291,7 +287,7 @@ public:
const size_t bufferCount = static_cast<size_t>(rays.Buffers.size());
for (size_t i = 0; i < bufferCount; ++i)
{
ChannelBufferOperations::Compact(rays.Buffers[i], masks, rays.NumRays, Device());
ChannelBufferOperations::Compact(rays.Buffers[i], masks, rays.NumRays);
}
return masks;
}
@ -340,12 +336,11 @@ public:
}
}
template <typename Device, typename T>
static void CopyDistancesToMin(Ray<T> rays, Device, const T offset = 0.f)
template <typename T>
static void CopyDistancesToMin(Ray<T> rays, const T offset = 0.f)
{
vtkm::worklet::DispatcherMapField<CopyAndOffsetMask<T>> dispatcher{ (
CopyAndOffsetMask<T>{ offset, RAY_EXITED_MESH }) };
dispatcher.SetDevice(Device());
dispatcher.Invoke(rays.Distance, rays.MinDistance, rays.Status);
}
};

44
vtkm/rendering/raytracing/SphereIntersector.cxx
View File

@ -141,6 +141,8 @@ public:
IdArrayPortal PointIds;
FloatPortal Radii;
SphereLeafIntersector() {}
SphereLeafIntersector(const IdHandle& pointIds, const FloatHandle& radii)
: PointIds(pointIds.PrepareForInput(Device()))
, Radii(radii.PrepareForInput(Device()))
@ -193,17 +195,25 @@ public:
}
};
struct IntersectFunctor
class SphereLeafWrapper : public vtkm::cont::ExecutionObjectBase
{
template <typename Device, typename Precision>
VTKM_CONT bool operator()(Device,
SphereIntersector* self,
Ray<Precision>& rays,
bool returnCellIndex)
protected:
using IdHandle = vtkm::cont::ArrayHandle<vtkm::Id>;
using FloatHandle = vtkm::cont::ArrayHandle<vtkm::Float32>;
IdHandle PointIds;
FloatHandle Radii;
public:
SphereLeafWrapper(IdHandle& pointIds, FloatHandle radii)
: PointIds(pointIds)
, Radii(radii)
{
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
self->IntersectRaysImp(Device(), rays, returnCellIndex);
return true;
}
template <typename Device>
VTKM_CONT SphereLeafIntersector<Device> PrepareForExecution(Device) const
{
return SphereLeafIntersector<Device>(PointIds, Radii);
}
};
@ -321,24 +331,22 @@ void SphereIntersector::SetData(const vtkm::cont::CoordinateSystem& coords,
void SphereIntersector::IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex)
{
vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
IntersectRaysImp(rays, returnCellIndex);
}
void SphereIntersector::IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex)
{
vtkm::cont::TryExecute(detail::IntersectFunctor(), this, rays, returnCellIndex);
IntersectRaysImp(rays, returnCellIndex);
}
template <typename Device, typename Precision>
void SphereIntersector::IntersectRaysImp(Device,
Ray<Precision>& rays,
bool vtkmNotUsed(returnCellIndex))
template <typename Precision>
void SphereIntersector::IntersectRaysImp(Ray<Precision>& rays, bool vtkmNotUsed(returnCellIndex))
{
detail::SphereLeafIntersector<Device> leafIntersector(this->PointIds, Radii);
detail::SphereLeafWrapper leafIntersector(this->PointIds, Radii);
BVHTraverser<detail::SphereLeafIntersector> traverser;
traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle, Device());
BVHTraverser traverser;
traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle);
RayOperations::UpdateRayStatus(rays);
}

4
vtkm/rendering/raytracing/SphereIntersector.h
View File

@ -48,8 +48,8 @@ public:
void IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex = false) override;
template <typename Device, typename Precision>
void IntersectRaysImp(Device, Ray<Precision>& rays, bool returnCellIndex);
template <typename Precision>
void IntersectRaysImp(Ray<Precision>& rays, bool returnCellIndex);
template <typename Precision>

891
vtkm/rendering/raytracing/TriangleIntersector.h
View File

@ -19,16 +19,11 @@
//============================================================================
#ifndef vtk_m_rendering_raytracing_TriagnleIntersector_h
#define vtk_m_rendering_raytracing_TriagnleIntersector_h
#include <cstring>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCompositeVector.h>
#include <vtkm/rendering/internal/RunTriangulator.h>
#include <vtkm/rendering/raytracing/BVHTraverser.h>
#include <vtkm/rendering/raytracing/BoundingVolumeHierarchy.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/rendering/raytracing/Ray.h>
#include <vtkm/rendering/raytracing/RayOperations.h>
#include <vtkm/rendering/raytracing/ShapeIntersector.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/rendering/vtkm_rendering_export.h>
namespace vtkm
{
@ -37,896 +32,44 @@ namespace rendering
namespace raytracing
{
class Moller
{
public:
template <typename Precision>
VTKM_EXEC void IntersectTri(const vtkm::Vec<Precision, 3>& a,
const vtkm::Vec<Precision, 3>& b,
const vtkm::Vec<Precision, 3>& c,
const vtkm::Vec<Precision, 3>& dir,
Precision& distance,
Precision& u,
Precision& v,
const vtkm::Vec<Precision, 3>& origin) const
{
const vtkm::Float32 EPSILON2 = 0.0001f;
vtkm::Vec<Precision, 3> e1 = b - a;
vtkm::Vec<Precision, 3> e2 = c - a;
vtkm::Vec<Precision, 3> p;
p[0] = dir[1] * e2[2] - dir[2] * e2[1];
p[1] = dir[2] * e2[0] - dir[0] * e2[2];
p[2] = dir[0] * e2[1] - dir[1] * e2[0];
Precision dot = e1[0] * p[0] + e1[1] * p[1] + e1[2] * p[2];
if (dot != 0.f)
{
dot = 1.f / dot;
vtkm::Vec<Precision, 3> t;
t = origin - a;
u = (t[0] * p[0] + t[1] * p[1] + t[2] * p[2]) * dot;
if (u >= (0.f - EPSILON2) && u <= (1.f + EPSILON2))
{
vtkm::Vec<Precision, 3> q; // = t % e1;
q[0] = t[1] * e1[2] - t[2] * e1[1];
q[1] = t[2] * e1[0] - t[0] * e1[2];
q[2] = t[0] * e1[1] - t[1] * e1[0];
v = (dir[0] * q[0] + dir[1] * q[1] + dir[2] * q[2]) * dot;
if (v >= (0.f - EPSILON2) && v <= (1.f + EPSILON2) && !(u + v > 1.f))
{
distance = (e2[0] * q[0] + e2[1] * q[1] + e2[2] * q[2]) * dot;
}
}
}
}
}; //Moller
// TODO: optimization for sorting ray dims before this call.
// This is called multiple times and kz,kx, and ky are
// constant for the ray
class WaterTight
{
public:
template <typename Precision>
VTKM_EXEC inline void FindDir(const vtkm::Vec<Precision, 3>& dir,
vtkm::Vec<Precision, 3>& s,
vtkm::Vec<Int32, 3>& k) const
{
//Find max ray direction
k[2] = 0;
if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[1]))
{
if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[2]))
k[2] = 0;
else
k[2] = 2;
}
else
{
if (vtkm::Abs(dir[1]) > vtkm::Abs(dir[2]))
k[2] = 1;
else
k[2] = 2;
}
k[0] = k[2] + 1;
if (k[0] == 3)
k[0] = 0;
k[1] = k[0] + 1;
if (k[1] == 3)
k[1] = 0;
if (dir[k[2]] < 0.f)
{
vtkm::Int32 temp = k[1];
k[1] = k[0];
k[0] = temp;
}
s[0] = dir[k[0]] / dir[k[2]];
s[1] = dir[k[1]] / dir[k[2]];
s[2] = 1.f / dir[k[2]];
}
template <typename Precision>
VTKM_EXEC_CONT inline void IntersectTri(const vtkm::Vec<Precision, 3>& a,
const vtkm::Vec<Precision, 3>& b,
const vtkm::Vec<Precision, 3>& c,
const vtkm::Vec<Precision, 3>& dir,
Precision& distance,
Precision& u,
Precision& v,
const vtkm::Vec<Precision, 3>& origin) const
{
vtkm::Vec<Int32, 3> k;
vtkm::Vec<Precision, 3> s;
//Find max ray direction
k[2] = 0;
if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[1]))
{
if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[2]))
k[2] = 0;
else
k[2] = 2;
}
else
{
if (vtkm::Abs(dir[1]) > vtkm::Abs(dir[2]))
k[2] = 1;
else
k[2] = 2;
}
k[0] = k[2] + 1;
if (k[0] == 3)
k[0] = 0;
k[1] = k[0] + 1;
if (k[1] == 3)
k[1] = 0;
if (dir[k[2]] < 0.f)
{
vtkm::Int32 temp = k[1];
k[1] = k[0];
k[0] = temp;
}
s[0] = dir[k[0]] / dir[k[2]];
s[1] = dir[k[1]] / dir[k[2]];
s[2] = 1.f / dir[k[2]];
vtkm::Vec<Precision, 3> A, B, C;
A = a - origin;
B = b - origin;
C = c - origin;
const Precision Ax = A[k[0]] - s[0] * A[k[2]];
const Precision Ay = A[k[1]] - s[1] * A[k[2]];
const Precision Bx = B[k[0]] - s[0] * B[k[2]];
const Precision By = B[k[1]] - s[1] * B[k[2]];
const Precision Cx = C[k[0]] - s[0] * C[k[2]];
const Precision Cy = C[k[1]] - s[1] * C[k[2]];
//scaled barycentric coords
u = Cx * By - Cy * Bx;
v = Ax * Cy - Ay * Cx;
Precision w = Bx * Ay - By * Ax;
if (u == 0.f || v == 0.f || w == 0.f)
{
vtkm::Float64 CxBy = vtkm::Float64(Cx) * vtkm::Float64(By);
vtkm::Float64 CyBx = vtkm::Float64(Cy) * vtkm::Float64(Bx);
u = vtkm::Float32(CxBy - CyBx);
vtkm::Float64 AxCy = vtkm::Float64(Ax) * vtkm::Float64(Cy);
vtkm::Float64 AyCx = vtkm::Float64(Ay) * vtkm::Float64(Cx);
v = vtkm::Float32(AxCy - AyCx);
vtkm::Float64 BxAy = vtkm::Float64(Bx) * vtkm::Float64(Ay);
vtkm::Float64 ByAx = vtkm::Float64(By) * vtkm::Float64(Ax);
w = vtkm::Float32(BxAy - ByAx);
}
Precision low = vtkm::Min(u, vtkm::Min(v, w));
Precision high = vtkm::Max(u, vtkm::Max(v, w));
bool invalid = (low < 0.) && (high > 0.);
Precision det = u + v + w;
if (det == 0.)
invalid = true;
const Precision Az = s[2] * A[k[2]];
const Precision Bz = s[2] * B[k[2]];
const Precision Cz = s[2] * C[k[2]];
det = 1.f / det;
u = u * det;
v = v * det;
distance = (u * Az + v * Bz + w * det * Cz);
u = v;
v = w * det;
if (invalid)
distance = -1.;
}
template <typename Precision>
VTKM_EXEC inline void IntersectTriSn(const vtkm::Vec<Precision, 3>& a,
const vtkm::Vec<Precision, 3>& b,
const vtkm::Vec<Precision, 3>& c,
const vtkm::Vec<Precision, 3>& s,
const vtkm::Vec<Int32, 3>& k,
Precision& distance,
Precision& u,
Precision& v,
const vtkm::Vec<Precision, 3>& origin) const
{
vtkm::Vec<Precision, 3> A, B, C;
A = a - origin;
B = b - origin;
C = c - origin;
const Precision Ax = A[k[0]] - s[0] * A[k[2]];
const Precision Ay = A[k[1]] - s[1] * A[k[2]];
const Precision Bx = B[k[0]] - s[0] * B[k[2]];
const Precision By = B[k[1]] - s[1] * B[k[2]];
const Precision Cx = C[k[0]] - s[0] * C[k[2]];
const Precision Cy = C[k[1]] - s[1] * C[k[2]];
//scaled barycentric coords
u = Cx * By - Cy * Bx;
v = Ax * Cy - Ay * Cx;
Precision w = Bx * Ay - By * Ax;
if (u == 0.f || v == 0.f || w == 0.f)
{
vtkm::Float64 CxBy = vtkm::Float64(Cx) * vtkm::Float64(By);
vtkm::Float64 CyBx = vtkm::Float64(Cy) * vtkm::Float64(Bx);
u = vtkm::Float32(CxBy - CyBx);
vtkm::Float64 AxCy = vtkm::Float64(Ax) * vtkm::Float64(Cy);
vtkm::Float64 AyCx = vtkm::Float64(Ay) * vtkm::Float64(Cx);
v = vtkm::Float32(AxCy - AyCx);
vtkm::Float64 BxAy = vtkm::Float64(Bx) * vtkm::Float64(Ay);
vtkm::Float64 ByAx = vtkm::Float64(By) * vtkm::Float64(Ax);
w = vtkm::Float32(BxAy - ByAx);
}
Precision low = vtkm::Min(u, vtkm::Min(v, w));
Precision high = vtkm::Max(u, vtkm::Max(v, w));
bool invalid = (low < 0.) && (high > 0.);
Precision det = u + v + w;
if (det == 0.)
invalid = true;
const Precision Az = s[2] * A[k[2]];
const Precision Bz = s[2] * B[k[2]];
const Precision Cz = s[2] * C[k[2]];
det = 1.f / det;
u = u * det;
v = v * det;
distance = (u * Az + v * Bz + w * det * Cz);
u = v;
v = w * det;
if (invalid)
distance = -1.;
}
}; //WaterTight
template <>
VTKM_EXEC inline void WaterTight::IntersectTri<vtkm::Float64>(
const vtkm::Vec<vtkm::Float64, 3>& a,
const vtkm::Vec<vtkm::Float64, 3>& b,
const vtkm::Vec<vtkm::Float64, 3>& c,
const vtkm::Vec<vtkm::Float64, 3>& dir,
vtkm::Float64& distance,
vtkm::Float64& u,
vtkm::Float64& v,
const vtkm::Vec<vtkm::Float64, 3>& origin) const
{
//Find max ray direction
int kz = 0;
if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[1]))
{
if (vtkm::Abs(dir[0]) > vtkm::Abs(dir[2]))
kz = 0;
else
kz = 2;
}
else
{
if (vtkm::Abs(dir[1]) > vtkm::Abs(dir[2]))
kz = 1;
else
kz = 2;
}
vtkm::Int32 kx = kz + 1;
if (kx == 3)
kx = 0;
vtkm::Int32 ky = kx + 1;
if (ky == 3)
ky = 0;
if (dir[kz] < 0.f)
{
vtkm::Int32 temp = ky;
ky = kx;
kx = temp;
}
vtkm::Float64 Sx = dir[kx] / dir[kz];
vtkm::Float64 Sy = dir[ky] / dir[kz];
vtkm::Float64 Sz = 1. / dir[kz];
vtkm::Vec<vtkm::Float64, 3> A, B, C;
A = a - origin;
B = b - origin;
C = c - origin;
const vtkm::Float64 Ax = A[kx] - Sx * A[kz];
const vtkm::Float64 Ay = A[ky] - Sy * A[kz];
const vtkm::Float64 Bx = B[kx] - Sx * B[kz];
const vtkm::Float64 By = B[ky] - Sy * B[kz];
const vtkm::Float64 Cx = C[kx] - Sx * C[kz];
const vtkm::Float64 Cy = C[ky] - Sy * C[kz];
//scaled barycentric coords
u = Cx * By - Cy * Bx;
v = Ax * Cy - Ay * Cx;
vtkm::Float64 w = Bx * Ay - By * Ax;
vtkm::Float64 low = vtkm::Min(u, vtkm::Min(v, w));
vtkm::Float64 high = vtkm::Max(u, vtkm::Max(v, w));
bool invalid = (low < 0.) && (high > 0.);
vtkm::Float64 det = u + v + w;
if (det == 0.)
invalid = true;
const vtkm::Float64 Az = Sz * A[kz];
const vtkm::Float64 Bz = Sz * B[kz];
const vtkm::Float64 Cz = Sz * C[kz];
det = 1. / det;
u = u * det;
v = v * det;
distance = (u * Az + v * Bz + w * det * Cz);
u = v;
v = w * det;
if (invalid)
distance = -1.;
}
template <typename Device>
//class WaterTightTriLeafIntersector : public vtkm::exec::ExecutionObjectBase
class WaterTightTriLeafIntersector
{
public:
using Id4Handle = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>>;
using Id4ArrayPortal = typename Id4Handle::ExecutionTypes<Device>::PortalConst;
Id4ArrayPortal Triangles;
public:
WaterTightTriLeafIntersector(const Id4Handle& triangles)
: Triangles(triangles.PrepareForInput(Device()))
{
}
WaterTightTriLeafIntersector(const WaterTightTriLeafIntersector<Device>& other)
: Triangles(other.Triangles)
{
}
template <typename PointPortalType, typename LeafPortalType, typename Precision>
VTKM_EXEC inline void IntersectLeaf(const vtkm::Int32& currentNode,
const vtkm::Vec<Precision, 3>& origin,
const vtkm::Vec<Precision, 3>& dir,
const PointPortalType& points,
vtkm::Id& hitIndex,
Precision& closestDistance,
Precision& minU,
Precision& minV,
LeafPortalType leafs,
const Precision& minDistance) const
{
const vtkm::Id triangleCount = leafs.Get(currentNode);
WaterTight intersector;
for (vtkm::Id i = 1; i <= triangleCount; ++i)
{
const vtkm::Id triIndex = leafs.Get(currentNode + i);
vtkm::Vec<Id, 4> triangle = Triangles.Get(triIndex);
vtkm::Vec<Precision, 3> a = vtkm::Vec<Precision, 3>(points.Get(triangle[1]));
vtkm::Vec<Precision, 3> b = vtkm::Vec<Precision, 3>(points.Get(triangle[2]));
vtkm::Vec<Precision, 3> c = vtkm::Vec<Precision, 3>(points.Get(triangle[3]));
Precision distance = -1.;
Precision u, v;
intersector.IntersectTri(a, b, c, dir, distance, u, v, origin);
if (distance != -1. && distance < closestDistance && distance > minDistance)
{
closestDistance = distance;
minU = u;
minV = v;
hitIndex = triIndex;
}
} // for
}
};
template <typename Device>
//class MollerTriLeafIntersector : public vtkm::exec::ExecutionObjectBase
class MollerTriLeafIntersector
{
//protected:
public:
using Id4Handle = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>>;
using Id4ArrayPortal = typename Id4Handle::ExecutionTypes<Device>::PortalConst;
Id4ArrayPortal Triangles;
public:
MollerTriLeafIntersector(const Id4Handle& triangles)
: Triangles(triangles.PrepareForInput(Device()))
{
}
template <typename PointPortalType, typename LeafPortalType, typename Precision>
VTKM_EXEC inline void IntersectLeaf(const vtkm::Int32& currentNode,
const vtkm::Vec<Precision, 3>& origin,
const vtkm::Vec<Precision, 3>& dir,
const PointPortalType& points,
vtkm::Id& hitIndex,
Precision& closestDistance,
Precision& minU,
Precision& minV,
LeafPortalType leafs,
const Precision& minDistance) const
{
const vtkm::Id triangleCount = leafs.Get(currentNode);
Moller intersector;
for (vtkm::Id i = 1; i <= triangleCount; ++i)
{
const vtkm::Id triIndex = leafs.Get(currentNode + i);
vtkm::Vec<Id, 4> triangle = Triangles.Get(triIndex);
vtkm::Vec<Precision, 3> a = vtkm::Vec<Precision, 3>(points.Get(triangle[1]));
vtkm::Vec<Precision, 3> b = vtkm::Vec<Precision, 3>(points.Get(triangle[2]));
vtkm::Vec<Precision, 3> c = vtkm::Vec<Precision, 3>(points.Get(triangle[3]));
Precision distance = -1.;
Precision u, v;
intersector.IntersectTri(a, b, c, dir, distance, u, v, origin);
if (distance != -1. && distance < closestDistance && distance > minDistance)
{
closestDistance = distance;
minU = u;
minV = v;
hitIndex = triIndex;
}
} // for
}
};
template <typename T>
VTKM_EXEC inline void swap(T& a, T& b)
{
T tmp = a;
a = b;
b = tmp;
}
VTKM_EXEC
inline vtkm::Float32 up(const vtkm::Float32& a)
{
return (a > 0.f) ? a * (1.f + vtkm::Float32(2e-23)) : a * (1.f - vtkm::Float32(2e-23));
}
VTKM_EXEC
inline vtkm::Float32 down(const vtkm::Float32& a)
{
return (a > 0.f) ? a * (1.f - vtkm::Float32(2e-23)) : a * (1.f + vtkm::Float32(2e-23));
}
VTKM_EXEC
inline vtkm::Float32 upFast(const vtkm::Float32& a)
{
return a * (1.f + vtkm::Float32(2e-23));
}
VTKM_EXEC
inline vtkm::Float32 downFast(const vtkm::Float32& a)
{
return a * (1.f - vtkm::Float32(2e-23));
}
namespace detail
{
class TriangleIntersectionData
{
public:
// Worklet to calutate the normals of a triagle if
// none are stored in the data set
class CalculateNormals : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
CalculateNormals() {}
typedef void ControlSignature(FieldIn<>,
FieldIn<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
WholeArrayIn<Vec3RenderingTypes>,
WholeArrayIn<>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7);
template <typename Precision, typename PointPortalType, typename IndicesPortalType>
VTKM_EXEC inline void operator()(const vtkm::Id& hitIndex,
const vtkm::Vec<Precision, 3>& rayDir,
Precision& normalX,
Precision& normalY,
Precision& normalZ,
const PointPortalType& points,
const IndicesPortalType& indicesPortal) const
{
if (hitIndex < 0)
return;
vtkm::Vec<Id, 4> indices = indicesPortal.Get(hitIndex);
vtkm::Vec<Precision, 3> a = points.Get(indices[1]);
vtkm::Vec<Precision, 3> b = points.Get(indices[2]);
vtkm::Vec<Precision, 3> c = points.Get(indices[3]);
vtkm::Vec<Precision, 3> normal = vtkm::TriangleNormal(a, b, c);
vtkm::Normalize(normal);
//flip the normal if its pointing the wrong way
if (vtkm::dot(normal, rayDir) > 0.f)
normal = -normal;
normalX = normal[0];
normalY = normal[1];
normalZ = normal[2];
}
}; //class CalculateNormals
template <typename Precision>
class LerpScalar : public vtkm::worklet::WorkletMapField
{
private:
Precision MinScalar;
Precision invDeltaScalar;
public:
VTKM_CONT
LerpScalar(const vtkm::Float32& minScalar, const vtkm::Float32& maxScalar)
: MinScalar(minScalar)
{
//Make sure the we don't divide by zero on
//something like an iso-surface
if (maxScalar - MinScalar != 0.f)
invDeltaScalar = 1.f / (maxScalar - MinScalar);
else
invDeltaScalar = 0.f;
}
typedef void ControlSignature(FieldIn<>,
FieldIn<>,
FieldIn<>,
FieldInOut<>,
WholeArrayIn<ScalarRenderingTypes>,
WholeArrayIn<>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6);
template <typename ScalarPortalType, typename IndicesPortalType>
VTKM_EXEC void operator()(const vtkm::Id& hitIndex,
const Precision& u,
const Precision& v,
Precision& lerpedScalar,
const ScalarPortalType& scalars,
const IndicesPortalType& indicesPortal) const
{
if (hitIndex < 0)
return;
vtkm::Vec<Id, 4> indices = indicesPortal.Get(hitIndex);
Precision n = 1.f - u - v;
Precision aScalar = Precision(scalars.Get(indices[1]));
Precision bScalar = Precision(scalars.Get(indices[2]));
Precision cScalar = Precision(scalars.Get(indices[3]));
lerpedScalar = aScalar * n + bScalar * u + cScalar * v;
//normalize
lerpedScalar = (lerpedScalar - MinScalar) * invDeltaScalar;
}
}; //class LerpScalar
template <typename Precision>
class NodalScalar : public vtkm::worklet::WorkletMapField
{
private:
Precision MinScalar;
Precision invDeltaScalar;
public:
VTKM_CONT
NodalScalar(const vtkm::Float32& minScalar, const vtkm::Float32& maxScalar)
: MinScalar(minScalar)
{
//Make sure the we don't divide by zero on
//something like an iso-surface
if (maxScalar - MinScalar != 0.f)
invDeltaScalar = 1.f / (maxScalar - MinScalar);
else
invDeltaScalar = 1.f / minScalar;
}
typedef void ControlSignature(FieldIn<>,
FieldOut<>,
WholeArrayIn<ScalarRenderingTypes>,
WholeArrayIn<>);
typedef void ExecutionSignature(_1, _2, _3, _4);
template <typename ScalarPortalType, typename IndicesPortalType>
VTKM_EXEC void operator()(const vtkm::Id& hitIndex,
Precision& scalar,
const ScalarPortalType& scalars,
const IndicesPortalType& indicesPortal) const
{
if (hitIndex < 0)
return;
vtkm::Vec<Id, 4> indices = indicesPortal.Get(hitIndex);
//Todo: one normalization
scalar = Precision(scalars.Get(indices[0]));
//normalize
scalar = (scalar - MinScalar) * invDeltaScalar;
}
}; //class LerpScalar
template <typename Precision>
VTKM_CONT void Run(Ray<Precision>& rays,
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> triangles,
vtkm::cont::CoordinateSystem coordsHandle,
const vtkm::cont::Field* scalarField,
const vtkm::Range& scalarRange)
{
bool isSupportedField =
(scalarField->GetAssociation() == vtkm::cont::Field::Association::POINTS ||
scalarField->GetAssociation() == vtkm::cont::Field::Association::CELL_SET);
if (!isSupportedField)
throw vtkm::cont::ErrorBadValue("Field not accociated with cell set or points");
bool isAssocPoints = scalarField->GetAssociation() == vtkm::cont::Field::Association::POINTS;
// Find the triangle normal
vtkm::worklet::DispatcherMapField<CalculateNormals>(CalculateNormals())
.Invoke(
rays.HitIdx, rays.Dir, rays.NormalX, rays.NormalY, rays.NormalZ, coordsHandle, triangles);
// Calculate scalar value at intersection point
if (isAssocPoints)
{
vtkm::worklet::DispatcherMapField<LerpScalar<Precision>>(
LerpScalar<Precision>(vtkm::Float32(scalarRange.Min), vtkm::Float32(scalarRange.Max)))
.Invoke(rays.HitIdx, rays.U, rays.V, rays.Scalar, *scalarField, triangles);
}
else
{
vtkm::worklet::DispatcherMapField<NodalScalar<Precision>>(
NodalScalar<Precision>(vtkm::Float32(scalarRange.Min), vtkm::Float32(scalarRange.Max)))
.Invoke(rays.HitIdx, rays.Scalar, *scalarField, triangles);
}
} // Run
}; // Class IntersectionData
#define AABB_EPSILON 0.00001f
class FindTriangleAABBs : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
FindTriangleAABBs() {}
typedef void ControlSignature(FieldIn<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
FieldOut<>,
WholeArrayIn<Vec3RenderingTypes>);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7, _8);
template <typename PointPortalType>
VTKM_EXEC void operator()(const vtkm::Vec<vtkm::Id, 4> indices,
vtkm::Float32& xmin,
vtkm::Float32& ymin,
vtkm::Float32& zmin,
vtkm::Float32& xmax,
vtkm::Float32& ymax,
vtkm::Float32& zmax,
const PointPortalType& points) const
{
// cast to Float32
vtkm::Vec<vtkm::Float32, 3> point;
point = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(indices[1]));
xmin = point[0];
ymin = point[1];
zmin = point[2];
xmax = xmin;
ymax = ymin;
zmax = zmin;
point = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(indices[2]));
xmin = vtkm::Min(xmin, point[0]);
ymin = vtkm::Min(ymin, point[1]);
zmin = vtkm::Min(zmin, point[2]);
xmax = vtkm::Max(xmax, point[0]);
ymax = vtkm::Max(ymax, point[1]);
zmax = vtkm::Max(zmax, point[2]);
point = static_cast<vtkm::Vec<vtkm::Float32, 3>>(points.Get(indices[3]));
xmin = vtkm::Min(xmin, point[0]);
ymin = vtkm::Min(ymin, point[1]);
zmin = vtkm::Min(zmin, point[2]);
xmax = vtkm::Max(xmax, point[0]);
ymax = vtkm::Max(ymax, point[1]);
zmax = vtkm::Max(zmax, point[2]);
vtkm::Float32 xEpsilon, yEpsilon, zEpsilon;
const vtkm::Float32 minEpsilon = 1e-6f;
xEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (xmax - xmin));
yEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (ymax - ymin));
zEpsilon = vtkm::Max(minEpsilon, AABB_EPSILON * (zmax - zmin));
xmin -= xEpsilon;
ymin -= yEpsilon;
zmin -= zEpsilon;
xmax += xEpsilon;
ymax += yEpsilon;
zmax += zEpsilon;
}
}; //class FindAABBs
#undef AABB_EPSILON
} // namespace detail
class TriangleIntersector : public ShapeIntersector
class VTKM_RENDERING_EXPORT TriangleIntersector : public ShapeIntersector
{
protected:
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> Triangles;
bool UseWaterTight;
public:
TriangleIntersector()
: UseWaterTight(false)
{
}
TriangleIntersector();
void SetUseWaterTight(bool useIt) { UseWaterTight = useIt; }
void SetUseWaterTight(bool useIt);
void SetData(const vtkm::cont::CoordinateSystem& coords,
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> triangles)
{
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> triangles);
CoordsHandle = coords;
Triangles = triangles;
vtkm::rendering::raytracing::AABBs AABB;
vtkm::worklet::DispatcherMapField<detail::FindTriangleAABBs>(detail::FindTriangleAABBs())
.Invoke(Triangles,
AABB.xmins,
AABB.ymins,
AABB.zmins,
AABB.xmaxs,
AABB.ymaxs,
AABB.zmaxs,
CoordsHandle);
this->SetAABBs(AABB);
}
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> GetTriangles() { return Triangles; }
class CellIndexFilter : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
CellIndexFilter() {}
typedef void ControlSignature(FieldInOut<>, WholeArrayIn<>);
typedef void ExecutionSignature(_1, _2);
template <typename TrianglePortalType>
VTKM_EXEC void operator()(vtkm::Id& hitIndex, TrianglePortalType& triangles) const
{
vtkm::Id cellIndex = -1;
if (hitIndex != -1)
{
cellIndex = triangles.Get(hitIndex)[0];
}
hitIndex = cellIndex;
}
}; //class CellIndexFilter
struct IntersectFunctor
{
template <typename Device, typename Precision>
VTKM_CONT bool operator()(Device,
TriangleIntersector* self,
Ray<Precision>& rays,
bool returnCellIndex)
{
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
self->IntersectRays(rays, Device(), returnCellIndex);
return true;
}
};
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 4>> GetTriangles();
vtkm::Id GetNumberOfShapes() const override;
VTKM_CONT void IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex = false) override
{
vtkm::cont::TryExecute(IntersectFunctor(), this, rays, returnCellIndex);
}
VTKM_CONT void IntersectRays(Ray<vtkm::Float32>& rays, bool returnCellIndex = false) override;
VTKM_CONT void IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex = false) override;
VTKM_CONT void IntersectRays(Ray<vtkm::Float64>& rays, bool returnCellIndex = false) override
{
vtkm::cont::TryExecute(IntersectFunctor(), this, rays, returnCellIndex);
}
template <typename Precision, typename Device>
VTKM_CONT void IntersectRays(Ray<Precision>& rays,
Device vtkmNotUsed(Device),
bool returnCellIndex)
{
if (UseWaterTight)
{
WaterTightTriLeafIntersector<Device> leafIntersector(this->Triangles);
BVHTraverser<WaterTightTriLeafIntersector> traverser;
traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle, Device());
}
else
{
MollerTriLeafIntersector<Device> leafIntersector(this->Triangles);
BVHTraverser<MollerTriLeafIntersector> traverser;
traverser.IntersectRays(rays, this->BVH, leafIntersector, this->CoordsHandle, Device());
}
// Normally we return the index of the triangle hit,
// but in some cases we are only interested in the cell
if (returnCellIndex)
{
vtkm::worklet::DispatcherMapField<CellIndexFilter> cellIndexFilterDispatcher;
cellIndexFilterDispatcher.Invoke(rays.HitIdx, Triangles);
}
// Update ray status
RayOperations::UpdateRayStatus(rays);
}
VTKM_CONT void IntersectionData(Ray<vtkm::Float32>& rays,
const vtkm::cont::Field* scalarField,
const vtkm::Range& scalarRange) override
{
IntersectionDataImp(rays, scalarField, scalarRange);
}
const vtkm::Range& scalarRange) override;
VTKM_CONT void IntersectionData(Ray<vtkm::Float64>& rays,
const vtkm::cont::Field* scalarField,
const vtkm::Range& scalarRange) override
{
IntersectionDataImp(rays, scalarField, scalarRange);
}
const vtkm::Range& scalarRange) override;
template <typename Precision>
VTKM_CONT void IntersectRaysImp(Ray<Precision>& rays, bool returnCellIndex);
template <typename Precision>
VTKM_CONT void IntersectionDataImp(Ray<Precision>& rays,
const vtkm::cont::Field* scalarField,
const vtkm::Range& scalarRange)
{
ShapeIntersector::IntersectionPoint(rays);
detail::TriangleIntersectionData intData;
intData.Run(rays, this->Triangles, this->CoordsHandle, scalarField, scalarRange);
}
const vtkm::Range& scalarRange);
vtkm::Id GetNumberOfShapes() const override { return Triangles.GetNumberOfValues(); }
}; // class intersector
}
}