mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-16 17:22:55 +00:00
e28309f09b
If a global static array is declared with VTKM_EXEC_CONSTANT and the code is compiled by nvcc (for multibackend code) then the array is only accesible on the GPU. If for some reason a worklet fails on the cuda backend and it is re-executed on any of the CPU backends, it will continue to fail. We couldn't find a simple way to declare the array once and have it available on both CPU and GPU. The approach we are using here is to declare the arrays as static inside some "Get" function which is marked as VTKM_EXEC_CONT.
298 lines
10 KiB
C++
298 lines
10 KiB
C++
//============================================================================
|
|
// 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.
|
|
//
|
|
// Copyright 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
|
|
// Copyright 2015 UT-Battelle, LLC.
|
|
// Copyright 2015 Los Alamos National Security.
|
|
//
|
|
// Under the terms of Contract DE-NA0003525 with NTESS,
|
|
// the U.S. Government retains certain rights in this software.
|
|
//
|
|
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
|
|
// Laboratory (LANL), the U.S. Government retains certain rights in
|
|
// this software.
|
|
//============================================================================
|
|
#ifndef vtk_m_rendering_raytracing_MortonCodes_h
|
|
#define vtk_m_rendering_raytracing_MortonCodes_h
|
|
|
|
#include <vtkm/VectorAnalysis.h>
|
|
|
|
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
|
|
|
|
#include <vtkm/rendering/raytracing/CellTables.h>
|
|
|
|
#include <vtkm/worklet/DispatcherMapField.h>
|
|
#include <vtkm/worklet/DispatcherMapTopology.h>
|
|
#include <vtkm/worklet/WorkletMapField.h>
|
|
#include <vtkm/worklet/WorkletMapTopology.h>
|
|
|
|
namespace vtkm
|
|
{
|
|
namespace rendering
|
|
{
|
|
namespace raytracing
|
|
{
|
|
|
|
|
|
//Note: if this takes a long time. we could use a lookup table
|
|
//expands 10-bit unsigned int into 30 bits
|
|
VTKM_EXEC inline vtkm::UInt32 ExpandBits32(vtkm::UInt32 x32)
|
|
{
|
|
x32 = (x32 | (x32 << 16)) & 0x030000FF;
|
|
x32 = (x32 | (x32 << 8)) & 0x0300F00F;
|
|
x32 = (x32 | (x32 << 4)) & 0x030C30C3;
|
|
x32 = (x32 | (x32 << 2)) & 0x09249249;
|
|
return x32;
|
|
}
|
|
|
|
VTKM_EXEC inline vtkm::UInt64 ExpandBits64(vtkm::UInt32 x)
|
|
{
|
|
vtkm::UInt64 x64 = x & 0x1FFFFF;
|
|
x64 = (x64 | x64 << 32) & 0x1F00000000FFFF;
|
|
x64 = (x64 | x64 << 16) & 0x1F0000FF0000FF;
|
|
x64 = (x64 | x64 << 8) & 0x100F00F00F00F00F;
|
|
x64 = (x64 | x64 << 4) & 0x10c30c30c30c30c3;
|
|
x64 = (x64 | x64 << 2) & 0x1249249249249249;
|
|
|
|
return x64;
|
|
}
|
|
|
|
//Returns 30 bit morton code for coordinates for
|
|
//coordinates in the unit cude
|
|
VTKM_EXEC inline vtkm::UInt32 Morton3D(vtkm::Float32& x, vtkm::Float32& y, vtkm::Float32& z)
|
|
{
|
|
//take the first 10 bits
|
|
x = vtkm::Min(vtkm::Max(x * 1024.0f, 0.0f), 1023.0f);
|
|
y = vtkm::Min(vtkm::Max(y * 1024.0f, 0.0f), 1023.0f);
|
|
z = vtkm::Min(vtkm::Max(z * 1024.0f, 0.0f), 1023.0f);
|
|
//expand the 10 bits to 30
|
|
vtkm::UInt32 xx = ExpandBits32((vtkm::UInt32)x);
|
|
vtkm::UInt32 yy = ExpandBits32((vtkm::UInt32)y);
|
|
vtkm::UInt32 zz = ExpandBits32((vtkm::UInt32)z);
|
|
//interleave coordinates
|
|
return (zz << 2 | yy << 1 | xx);
|
|
}
|
|
|
|
//Returns 30 bit morton code for coordinates for
|
|
//coordinates in the unit cude
|
|
VTKM_EXEC inline vtkm::UInt64 Morton3D64(vtkm::Float32& x, vtkm::Float32& y, vtkm::Float32& z)
|
|
{
|
|
//take the first 21 bits
|
|
x = vtkm::Min(vtkm::Max(x * 2097152.0f, 0.0f), 2097151.0f);
|
|
y = vtkm::Min(vtkm::Max(y * 2097152.0f, 0.0f), 2097151.0f);
|
|
z = vtkm::Min(vtkm::Max(z * 2097152.0f, 0.0f), 2097151.0f);
|
|
//expand the 10 bits to 30
|
|
vtkm::UInt64 xx = ExpandBits64((vtkm::UInt32)x);
|
|
vtkm::UInt64 yy = ExpandBits64((vtkm::UInt32)y);
|
|
vtkm::UInt64 zz = ExpandBits64((vtkm::UInt32)z);
|
|
//interleave coordinates
|
|
return (zz << 2 | yy << 1 | xx);
|
|
}
|
|
|
|
class MortonCodeFace : public vtkm::worklet::WorkletMapPointToCell
|
|
{
|
|
private:
|
|
// (1.f / dx),(1.f / dy), (1.f, / dz)
|
|
vtkm::Vec<vtkm::Float32, 3> InverseExtent;
|
|
vtkm::Vec<vtkm::Float32, 3> MinCoordinate;
|
|
|
|
VTKM_EXEC inline void Normalize(vtkm::Vec<vtkm::Float32, 3>& point) const
|
|
{
|
|
point = (point - MinCoordinate) * InverseExtent;
|
|
}
|
|
|
|
VTKM_EXEC inline void Sort4(vtkm::Vec<vtkm::Id, 4>& indices) const
|
|
{
|
|
if (indices[0] < indices[1])
|
|
{
|
|
vtkm::Id temp = indices[1];
|
|
indices[1] = indices[0];
|
|
indices[0] = temp;
|
|
}
|
|
if (indices[2] < indices[3])
|
|
{
|
|
vtkm::Id temp = indices[3];
|
|
indices[3] = indices[2];
|
|
indices[2] = temp;
|
|
}
|
|
if (indices[0] < indices[2])
|
|
{
|
|
vtkm::Id temp = indices[2];
|
|
indices[2] = indices[0];
|
|
indices[0] = temp;
|
|
}
|
|
if (indices[1] < indices[3])
|
|
{
|
|
vtkm::Id temp = indices[3];
|
|
indices[3] = indices[1];
|
|
indices[1] = temp;
|
|
}
|
|
if (indices[1] < indices[2])
|
|
{
|
|
vtkm::Id temp = indices[2];
|
|
indices[2] = indices[1];
|
|
indices[1] = temp;
|
|
}
|
|
}
|
|
|
|
public:
|
|
VTKM_CONT
|
|
MortonCodeFace(const vtkm::Vec<vtkm::Float32, 3>& inverseExtent,
|
|
const vtkm::Vec<vtkm::Float32, 3>& minCoordinate)
|
|
: InverseExtent(inverseExtent)
|
|
, MinCoordinate(minCoordinate)
|
|
{
|
|
}
|
|
|
|
typedef void ControlSignature(CellSetIn cellset,
|
|
WholeArrayIn<>,
|
|
FieldInTo<>,
|
|
WholeArrayOut<>,
|
|
WholeArrayOut<>);
|
|
|
|
typedef void ExecutionSignature(CellShape, FromIndices, WorkIndex, _2, _3, _4, _5);
|
|
|
|
template <typename CellShape,
|
|
typename CellNodeVecType,
|
|
typename PointPortalType,
|
|
typename MortonPortalType,
|
|
typename CellFaceIdsPortalType>
|
|
VTKM_EXEC inline void operator()(const CellShape& cellShape,
|
|
const CellNodeVecType& cellIndices,
|
|
const vtkm::Id& cellId,
|
|
const PointPortalType& points,
|
|
const vtkm::Id& offset,
|
|
MortonPortalType& mortonCodes,
|
|
CellFaceIdsPortalType& cellFaceIds) const
|
|
{
|
|
vtkm::Int32 faceCount;
|
|
vtkm::Int32 tableOffset;
|
|
|
|
if (cellShape.Id == vtkm::CELL_SHAPE_TETRA)
|
|
{
|
|
faceCount = CellTables::Get().FaceLookUp[1][1];
|
|
tableOffset = CellTables::Get().FaceLookUp[1][0];
|
|
}
|
|
else if (cellShape.Id == vtkm::CELL_SHAPE_HEXAHEDRON)
|
|
{
|
|
faceCount = CellTables::Get().FaceLookUp[0][1];
|
|
tableOffset = CellTables::Get().FaceLookUp[0][0];
|
|
}
|
|
else if (cellShape.Id == vtkm::CELL_SHAPE_WEDGE)
|
|
{
|
|
faceCount = CellTables::Get().FaceLookUp[2][1];
|
|
tableOffset = CellTables::Get().FaceLookUp[2][0];
|
|
}
|
|
else if (cellShape.Id == vtkm::CELL_SHAPE_PYRAMID)
|
|
{
|
|
faceCount = CellTables::Get().FaceLookUp[3][1];
|
|
tableOffset = CellTables::Get().FaceLookUp[3][0];
|
|
}
|
|
else
|
|
{
|
|
printf("Unknown shape type %d\n", (int)cellShape.Id);
|
|
return;
|
|
}
|
|
|
|
//calc the morton code at the center of each face
|
|
for (vtkm::Int32 i = 0; i < faceCount; ++i)
|
|
{
|
|
vtkm::Vec<vtkm::Float32, 3> center;
|
|
vtkm::UInt32 code;
|
|
vtkm::Vec<vtkm::Id, 3> cellFace;
|
|
cellFace[0] = cellId;
|
|
|
|
// We must be sure that this calulation is the same for all faces. If we didn't
|
|
// then it is possible for the same face to end up in multiple morton "buckets" due to
|
|
// the wonders of floating point math. This is bad. If we calculate in the same order
|
|
// for all faces, then at worst, two different faces can enter the same bucket, which
|
|
// we currently check for.
|
|
vtkm::Vec<vtkm::Id, 4> faceIndices;
|
|
faceIndices[3] = -1;
|
|
//Number of indices this face has
|
|
const vtkm::Id indiceCount = CellTables::Get().ShapesFaceList[tableOffset + i][0];
|
|
for (vtkm::Int32 j = 1; j <= indiceCount; j++)
|
|
{
|
|
faceIndices[j - 1] = cellIndices[CellTables::Get().ShapesFaceList[tableOffset + i][j]];
|
|
}
|
|
//sort the indices in descending order
|
|
Sort4(faceIndices);
|
|
|
|
vtkm::Int32 count = 1;
|
|
BOUNDS_CHECK(points, faceIndices[0]);
|
|
center = points.Get(faceIndices[0]);
|
|
for (int idx = 1; idx < indiceCount; ++idx)
|
|
{
|
|
BOUNDS_CHECK(points, faceIndices[idx]);
|
|
center = center + points.Get(faceIndices[idx]);
|
|
count++;
|
|
}
|
|
//TODO: we could make this a recipical, but this is not a bottleneck.
|
|
center[0] = center[0] / vtkm::Float32(count);
|
|
center[1] = center[1] / vtkm::Float32(count);
|
|
center[2] = center[2] / vtkm::Float32(count);
|
|
Normalize(center);
|
|
code = Morton3D(center[0], center[1], center[2]);
|
|
BOUNDS_CHECK(mortonCodes, offset + i);
|
|
mortonCodes.Set(offset + i, code);
|
|
cellFace[1] = i;
|
|
cellFace[2] = -1; //Need to initialize this for the next step
|
|
BOUNDS_CHECK(cellFaceIds, offset + i);
|
|
cellFaceIds.Set(offset + i, cellFace);
|
|
}
|
|
}
|
|
}; // class MortonCodeFace
|
|
|
|
class MortonCodeAABB : public vtkm::worklet::WorkletMapField
|
|
{
|
|
private:
|
|
// (1.f / dx),(1.f / dy), (1.f, / dz)
|
|
vtkm::Vec<vtkm::Float32, 3> InverseExtent;
|
|
vtkm::Vec<vtkm::Float32, 3> MinCoordinate;
|
|
|
|
public:
|
|
VTKM_CONT
|
|
MortonCodeAABB(const vtkm::Vec<vtkm::Float32, 3>& inverseExtent,
|
|
const vtkm::Vec<vtkm::Float32, 3>& minCoordinate)
|
|
: InverseExtent(inverseExtent)
|
|
, MinCoordinate(minCoordinate)
|
|
{
|
|
}
|
|
|
|
typedef void
|
|
ControlSignature(FieldIn<>, FieldIn<>, FieldIn<>, FieldIn<>, FieldIn<>, FieldIn<>, FieldOut<>);
|
|
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6, _7);
|
|
typedef _7 InputDomain;
|
|
|
|
VTKM_EXEC
|
|
void operator()(const vtkm::Float32& xmin,
|
|
const vtkm::Float32& ymin,
|
|
const vtkm::Float32& zmin,
|
|
const vtkm::Float32& xmax,
|
|
const vtkm::Float32& ymax,
|
|
const vtkm::Float32& zmax,
|
|
vtkm::UInt32& mortonCode) const
|
|
{
|
|
vtkm::Vec<vtkm::Float32, 3> direction(xmax - xmin, ymax - ymin, zmax - zmin);
|
|
vtkm::Float32 halfDistance = sqrtf(vtkm::dot(direction, direction)) * 0.5f;
|
|
vtkm::Normalize(direction);
|
|
vtkm::Float32 centroidx = xmin + halfDistance * direction[0] - MinCoordinate[0];
|
|
vtkm::Float32 centroidy = ymin + halfDistance * direction[1] - MinCoordinate[1];
|
|
vtkm::Float32 centroidz = zmin + halfDistance * direction[2] - MinCoordinate[2];
|
|
//normalize the centroid tp 10 bits
|
|
centroidx *= InverseExtent[0];
|
|
centroidy *= InverseExtent[1];
|
|
centroidz *= InverseExtent[2];
|
|
mortonCode = Morton3D(centroidx, centroidy, centroidz);
|
|
}
|
|
}; // class MortonCodeAABB
|
|
}
|
|
}
|
|
} //namespace vtkm::rendering::raytracing
|
|
#endif
|