rendering/raytracing/ConnectivityTracer.hxx

//============================================================================
//  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.
//============================================================================
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/ErrorBadValue.h>
#include <vtkm/cont/Timer.h>
#include <vtkm/cont/TryExecute.h>
#include <vtkm/cont/internal/DeviceAdapterListHelpers.h>

#include <vtkm/rendering/raytracing/Camera.h>
#include <vtkm/rendering/raytracing/CellIntersector.h>
#include <vtkm/rendering/raytracing/CellSampler.h>
#include <vtkm/rendering/raytracing/CellTables.h>
#include <vtkm/rendering/raytracing/MeshConnectivityBase.h>
#include <vtkm/rendering/raytracing/Ray.h>
#include <vtkm/rendering/raytracing/RayOperations.h>
#include <vtkm/rendering/raytracing/RayTracingTypeDefs.h>
#include <vtkm/rendering/raytracing/Worklets.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
{
namespace detail
{
template <typename FloatType>
template <typename Device>
void RayTracking<FloatType>::Compact(vtkm::cont::ArrayHandle<FloatType>& compactedDistances,
                                     vtkm::cont::ArrayHandle<UInt8>& masks,
                                     Device)
{
  //
  // These distances are stored in the rays, and it has
  // already been compacted.
  //
  CurrentDistance = compactedDistances;

  vtkm::cont::ArrayHandleCast<vtkm::Id, vtkm::cont::ArrayHandle<vtkm::UInt8>> castedMasks(masks);

  bool distance1IsEnter = EnterDist == &Distance1;

  vtkm::cont::ArrayHandle<FloatType> compactedDistance1;
  vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(Distance1, masks, compactedDistance1);
  Distance1 = compactedDistance1;

  vtkm::cont::ArrayHandle<FloatType> compactedDistance2;
  vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(Distance2, masks, compactedDistance2);
  Distance2 = compactedDistance2;

  vtkm::cont::ArrayHandle<vtkm::Int32> compactedExitFace;
  vtkm::cont::DeviceAdapterAlgorithm<Device>::CopyIf(ExitFace, masks, compactedExitFace);
  ExitFace = compactedExitFace;

  if (distance1IsEnter)
  {
    EnterDist = &Distance1;
    ExitDist = &Distance2;
  }
  else
  {
    EnterDist = &Distance2;
    ExitDist = &Distance1;
  }
}

template <typename FloatType>
template <typename Device>
void RayTracking<FloatType>::Init(const vtkm::Id size,
                                  vtkm::cont::ArrayHandle<FloatType>& distances,
                                  Device)
{

  ExitFace.PrepareForOutput(size, Device());
  Distance1.PrepareForOutput(size, Device());
  Distance2.PrepareForOutput(size, Device());

  CurrentDistance = distances;
  //
  // Set the initial Distances
  //
  vtkm::worklet::DispatcherMapField<CopyAndOffset<FloatType>> resetDistancesDispatcher(
    CopyAndOffset<FloatType>(0.0f));
  resetDistancesDispatcher.SetDevice(Device());
  resetDistancesDispatcher.Invoke(distances, *EnterDist);

  //
  // Init the exit faces. This value is used to load the next cell
  // base on the cell and face it left
  //
  vtkm::cont::ArrayHandleConstant<vtkm::Int32> negOne(-1, size);
  vtkm::cont::Algorithm::Copy(Device(), negOne, ExitFace);

  vtkm::cont::ArrayHandleConstant<FloatType> negOnef(-1.f, size);
  vtkm::cont::Algorithm::Copy(Device(), negOnef, *ExitDist);
}

template <typename FloatType>
void RayTracking<FloatType>::Swap()
{
  vtkm::cont::ArrayHandle<FloatType>* tmpPtr;
  tmpPtr = EnterDist;
  EnterDist = ExitDist;
  ExitDist = tmpPtr;
}
} //namespace detail

template <typename FloatType, typename Device>
void ConnectivityTracer::PrintRayStatus(Ray<FloatType>& rays, Device)
{
  vtkm::Id raysExited = RayOperations::GetStatusCount(rays, RAY_EXITED_MESH, Device());
  vtkm::Id raysActive = RayOperations::GetStatusCount(rays, RAY_ACTIVE, Device());
  vtkm::Id raysAbandoned = RayOperations::GetStatusCount(rays, RAY_ABANDONED, Device());
  vtkm::Id raysExitedDom = RayOperations::GetStatusCount(rays, RAY_EXITED_DOMAIN, Device());
  std::cout << "\r Ray Status " << std::setw(10) << std::left << " Lost " << std::setw(10)
            << std::left << RaysLost << std::setw(10) << std::left << " Exited " << std::setw(10)
            << std::left << raysExited << std::setw(10) << std::left << " Active " << std::setw(10)
            << raysActive << std::setw(10) << std::left << " Abandoned " << std::setw(10)
            << raysAbandoned << " Exited Domain " << std::setw(10) << std::left << raysExitedDom
            << "\n";
}

//
//  Advance Ray
//      After a ray leaves the mesh, we need to check to see
//      of the ray re-enters the mesh within this domain. This
//      function moves the ray forward some offset to prevent
//      "shadowing" and hitting the same exit point.
//
template <typename FloatType>
class AdvanceRay : public vtkm::worklet::WorkletMapField
{
  FloatType Offset;

public:
  VTKM_CONT
  AdvanceRay(const FloatType offset = 0.00001)
    : Offset(offset)
  {
  }
  using ControlSignature = void(FieldIn<>, FieldInOut<>);
  using ExecutionSignature = void(_1, _2);

  VTKM_EXEC inline void operator()(const vtkm::UInt8& status, FloatType& distance) const
  {
    if (status == RAY_EXITED_MESH)
      distance += Offset;
  }
}; //class AdvanceRay

template <typename FloatType>
class LocateCell : public vtkm::worklet::WorkletMapField
{
private:
  const MeshConnectivityBase* MeshConn;
  CellIntersector<255> Intersector;

public:
  LocateCell(const MeshConnectivityBase* meshConn)
    : MeshConn(meshConn)
  {
  }

  using ControlSignature = void(FieldInOut<>,
                                WholeArrayIn<>,
                                FieldIn<>,
                                FieldInOut<>,
                                FieldInOut<>,
                                FieldInOut<>,
                                FieldInOut<>,
                                FieldIn<>);
  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8);

  template <typename PointPortalType>
  VTKM_EXEC inline void operator()(vtkm::Id& currentCell,
                                   PointPortalType& vertices,
                                   const vtkm::Vec<FloatType, 3>& dir,
                                   FloatType& enterDistance,
                                   FloatType& exitDistance,
                                   vtkm::Int32& enterFace,
                                   vtkm::UInt8& rayStatus,
                                   const vtkm::Vec<FloatType, 3>& origin) const
  {
    if (enterFace != -1 && rayStatus == RAY_ACTIVE)
    {
      currentCell = MeshConn->GetConnectingCell(currentCell, enterFace);
      if (currentCell == -1)
        rayStatus = RAY_EXITED_MESH;
      enterFace = -1;
    }
    //This ray is dead or exited the mesh and needs re-entry
    if (rayStatus != RAY_ACTIVE)
    {
      return;
    }
    FloatType xpoints[8];
    FloatType ypoints[8];
    FloatType zpoints[8];
    vtkm::Id cellConn[8];
    FloatType distances[6];

    const vtkm::Int32 numIndices = MeshConn->GetCellIndices(cellConn, currentCell);
    //load local cell data
    for (int i = 0; i < numIndices; ++i)
    {
      BOUNDS_CHECK(vertices, cellConn[i]);
      vtkm::Vec<FloatType, 3> point = vtkm::Vec<FloatType, 3>(vertices.Get(cellConn[i]));
      xpoints[i] = point[0];
      ypoints[i] = point[1];
      zpoints[i] = point[2];
    }
    const vtkm::UInt8 cellShape = MeshConn->GetCellShape(currentCell);
    Intersector.IntersectCell(xpoints, ypoints, zpoints, dir, origin, distances, cellShape);

    CellTables tables;
    const vtkm::Int32 numFaces = tables.FaceLookUp(tables.CellTypeLookUp(cellShape), 1);
    //vtkm::Int32 minFace = 6;
    vtkm::Int32 maxFace = -1;

    FloatType minDistance = static_cast<FloatType>(1e32);
    FloatType maxDistance = static_cast<FloatType>(-1);
    int hitCount = 0;
    for (vtkm::Int32 i = 0; i < numFaces; ++i)
    {
      FloatType dist = distances[i];

      if (dist != -1)
      {
        hitCount++;
        if (dist < minDistance)
        {
          minDistance = dist;
          //minFace = i;
        }
        if (dist > maxDistance)
        {
          maxDistance = dist;
          maxFace = i;
        }
      }
    }

    if (maxDistance <= enterDistance || minDistance == maxDistance)
    {
      rayStatus = RAY_LOST;
    }
    else
    {
      enterDistance = minDistance;
      exitDistance = maxDistance;
      enterFace = maxFace;
    }

  } //operator
};  //class LocateCell

template <typename FloatType, typename Device>
class RayBumper : public vtkm::worklet::WorkletMapField
{
private:
  using FloatTypeHandle = typename vtkm::cont::ArrayHandle<FloatType>;
  using FloatTypePortal = typename FloatTypeHandle::template ExecutionTypes<Device>::Portal;
  FloatTypePortal DirectionsX;
  FloatTypePortal DirectionsY;
  FloatTypePortal DirectionsZ;
  const MeshConnectivityBase* MeshConn;
  CellIntersector<255> Intersector;
  const vtkm::UInt8 FailureStatus; // the status to assign ray if we fail to find the intersection
public:
  RayBumper(FloatTypeHandle dirsx,
            FloatTypeHandle dirsy,
            FloatTypeHandle dirsz,
            const MeshConnectivityBase* meshConn,
            vtkm::UInt8 failureStatus = RAY_ABANDONED)
    : DirectionsX(dirsx.PrepareForInPlace(Device()))
    , DirectionsY(dirsy.PrepareForInPlace(Device()))
    , DirectionsZ(dirsz.PrepareForInPlace(Device()))
    , MeshConn(meshConn)
    , FailureStatus(failureStatus)
  {
  }


  using ControlSignature = void(FieldInOut<>,
                                WholeArrayIn<>,
                                FieldInOut<>,
                                FieldInOut<>,
                                FieldInOut<>,
                                FieldInOut<>,
                                FieldIn<>);
  using ExecutionSignature = void(_1, _2, _3, _4, _5, WorkIndex, _6, _7);

  template <typename PointPortalType>
  VTKM_EXEC inline void operator()(vtkm::Id& currentCell,
                                   PointPortalType& vertices,
                                   FloatType& enterDistance,
                                   FloatType& exitDistance,
                                   vtkm::Int32& enterFace,
                                   const vtkm::Id& pixelIndex,
                                   vtkm::UInt8& rayStatus,
                                   const vtkm::Vec<FloatType, 3>& origin) const
  {
    // We only process lost rays
    if (rayStatus != RAY_LOST)
    {
      return;
    }

    FloatType xpoints[8];
    FloatType ypoints[8];
    FloatType zpoints[8];
    vtkm::Id cellConn[8];
    FloatType distances[6];

    vtkm::Vec<FloatType, 3> centroid(0., 0., 0.);

    const vtkm::Int32 numIndices = MeshConn->GetCellIndices(cellConn, currentCell);
    //load local cell data
    for (int i = 0; i < numIndices; ++i)
    {
      BOUNDS_CHECK(vertices, cellConn[i]);
      vtkm::Vec<FloatType, 3> point = vtkm::Vec<FloatType, 3>(vertices.Get(cellConn[i]));
      centroid = centroid + point;
      xpoints[i] = point[0];
      ypoints[i] = point[1];
      zpoints[i] = point[2];
    }

    FloatType invNumIndices = static_cast<FloatType>(1.) / static_cast<FloatType>(numIndices);
    centroid[0] = centroid[0] * invNumIndices;
    centroid[1] = centroid[1] * invNumIndices;
    centroid[2] = centroid[2] * invNumIndices;

    vtkm::Vec<FloatType, 3> toCentroid = centroid - origin;
    vtkm::Normalize(toCentroid);

    vtkm::Vec<FloatType, 3> dir(
      DirectionsX.Get(pixelIndex), DirectionsY.Get(pixelIndex), DirectionsZ.Get(pixelIndex));
    vtkm::Vec<FloatType, 3> bump = toCentroid - dir;
    dir = dir + RAY_TUG_EPSILON * bump;

    vtkm::Normalize(dir);

    DirectionsX.Set(pixelIndex, dir[0]);
    DirectionsY.Set(pixelIndex, dir[1]);
    DirectionsZ.Set(pixelIndex, dir[2]);

    const vtkm::UInt8 cellShape = MeshConn->GetCellShape(currentCell);
    Intersector.IntersectCell(xpoints, ypoints, zpoints, dir, origin, distances, cellShape);

    CellTables tables;
    const vtkm::Int32 numFaces = tables.FaceLookUp(tables.CellTypeLookUp(cellShape), 1);

    //vtkm::Int32 minFace = 6;
    vtkm::Int32 maxFace = -1;
    FloatType minDistance = static_cast<FloatType>(1e32);
    FloatType maxDistance = static_cast<FloatType>(-1);
    int hitCount = 0;
    for (int i = 0; i < numFaces; ++i)
    {
      FloatType dist = distances[i];

      if (dist != -1)
      {
        hitCount++;
        if (dist < minDistance)
        {
          minDistance = dist;
          //minFace = i;
        }
        if (dist >= maxDistance)
        {
          maxDistance = dist;
          maxFace = i;
        }
      }
    }
    if (minDistance >= maxDistance)
    {
      rayStatus = FailureStatus;
    }
    else
    {
      enterDistance = minDistance;
      exitDistance = maxDistance;
      enterFace = maxFace;
      rayStatus = RAY_ACTIVE; //re-activate ray
    }

  } //operator
};  //class RayBumper

template <typename FloatType>
class AddPathLengths : public vtkm::worklet::WorkletMapField
{
public:
  VTKM_CONT
  AddPathLengths() {}

  using ControlSignature = void(FieldIn<RayStatusType>,            // ray status
                                FieldIn<ScalarRenderingTypes>,     // cell enter distance
                                FieldIn<ScalarRenderingTypes>,     // cell exit distance
                                FieldInOut<ScalarRenderingTypes>); // ray absorption data

  using ExecutionSignature = void(_1, _2, _3, _4);

  VTKM_EXEC inline void operator()(const vtkm::UInt8& rayStatus,
                                   const FloatType& enterDistance,
                                   const FloatType& exitDistance,
                                   FloatType& distance) const
  {
    if (rayStatus != RAY_ACTIVE)
    {
      return;
    }

    if (exitDistance <= enterDistance)
    {
      return;
    }

    FloatType segmentLength = exitDistance - enterDistance;
    distance += segmentLength;
  }
};

template <typename FloatType>
class Integrate : public vtkm::worklet::WorkletMapField
{
private:
  const vtkm::Int32 NumBins;

public:
  VTKM_CONT
  Integrate(const vtkm::Int32 numBins)
    : NumBins(numBins)
  {
  }

  using ControlSignature = void(FieldIn<RayStatusType>,                // ray status
                                FieldIn<ScalarRenderingTypes>,         // cell enter distance
                                FieldIn<ScalarRenderingTypes>,         // cell exit distance
                                FieldInOut<ScalarRenderingTypes>,      // current distance
                                WholeArrayIn<ScalarRenderingTypes>,    // cell absorption data array
                                WholeArrayInOut<ScalarRenderingTypes>, // ray absorption data
                                FieldIn<IdType>);                      // current cell

  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, WorkIndex);

  template <typename CellDataPortalType, typename RayDataPortalType>
  VTKM_EXEC inline void operator()(const vtkm::UInt8& rayStatus,
                                   const FloatType& enterDistance,
                                   const FloatType& exitDistance,
                                   FloatType& currentDistance,
                                   const CellDataPortalType& cellData,
                                   RayDataPortalType& energyBins,
                                   const vtkm::Id& currentCell,
                                   const vtkm::Id& rayIndex) const
  {
    if (rayStatus != RAY_ACTIVE)
    {
      return;
    }
    if (exitDistance <= enterDistance)
    {
      return;
    }

    FloatType segmentLength = exitDistance - enterDistance;

    vtkm::Id rayOffset = NumBins * rayIndex;
    vtkm::Id cellOffset = NumBins * currentCell;
    for (vtkm::Int32 i = 0; i < NumBins; ++i)
    {
      BOUNDS_CHECK(cellData, cellOffset + i);
      FloatType absorb = static_cast<FloatType>(cellData.Get(cellOffset + i));

      absorb = vtkm::Exp(-absorb * segmentLength);
      BOUNDS_CHECK(energyBins, rayOffset + i);
      FloatType intensity = static_cast<FloatType>(energyBins.Get(rayOffset + i));
      energyBins.Set(rayOffset + i, intensity * absorb);
    }
    currentDistance = exitDistance;
  }
};

template <typename FloatType>
class IntegrateEmission : public vtkm::worklet::WorkletMapField
{
private:
  const vtkm::Int32 NumBins;
  bool DivideEmisByAbsorb;

public:
  VTKM_CONT
  IntegrateEmission(const vtkm::Int32 numBins, const bool divideEmisByAbsorb)
    : NumBins(numBins)
    , DivideEmisByAbsorb(divideEmisByAbsorb)
  {
  }

  using ControlSignature = void(FieldIn<>,                          // ray status
                                FieldIn<>,                          // cell enter distance
                                FieldIn<>,                          // cell exit distance
                                FieldInOut<>,                       // current distance
                                WholeArrayIn<ScalarRenderingTypes>, // cell absorption data array
                                WholeArrayIn<ScalarRenderingTypes>, // cell emission data array
                                WholeArrayInOut<>,                  // ray absorption data
                                WholeArrayInOut<>,                  // ray emission data
                                FieldIn<>);                         // current cell

  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, _9, WorkIndex);

  template <typename CellAbsPortalType, typename CellEmisPortalType, typename RayDataPortalType>
  VTKM_EXEC inline void operator()(const vtkm::UInt8& rayStatus,
                                   const FloatType& enterDistance,
                                   const FloatType& exitDistance,
                                   FloatType& currentDistance,
                                   const CellAbsPortalType& absorptionData,
                                   const CellEmisPortalType& emissionData,
                                   RayDataPortalType& absorptionBins,
                                   RayDataPortalType& emissionBins,
                                   const vtkm::Id& currentCell,
                                   const vtkm::Id& rayIndex) const
  {
    if (rayStatus != RAY_ACTIVE)
    {
      return;
    }
    if (exitDistance <= enterDistance)
    {
      return;
    }

    FloatType segmentLength = exitDistance - enterDistance;

    vtkm::Id rayOffset = NumBins * rayIndex;
    vtkm::Id cellOffset = NumBins * currentCell;
    for (vtkm::Int32 i = 0; i < NumBins; ++i)
    {
      BOUNDS_CHECK(absorptionData, cellOffset + i);
      FloatType absorb = static_cast<FloatType>(absorptionData.Get(cellOffset + i));
      BOUNDS_CHECK(emissionData, cellOffset + i);
      FloatType emission = static_cast<FloatType>(emissionData.Get(cellOffset + i));

      if (DivideEmisByAbsorb)
      {
        emission /= absorb;
      }

      FloatType tmp = vtkm::Exp(-absorb * segmentLength);
      BOUNDS_CHECK(absorptionBins, rayOffset + i);

      //
      // Traditionally, we would only keep track of a single intensity value per ray
      // per bin and we would integrate from the beginning to end of the ray. In a
      // distributed memory setting, we would move cell data around so that the
      // entire ray could be traced, but in situ, moving that much cell data around
      // could blow memory. Here we are keeping track of two values. Total absorption
      // through this contiguous segment of the mesh, and the amount of emitted energy
      // that makes it out of this mesh segment. If this is really run on a single node,
      // we can get the final energy value by multiplying the background intensity by
      // the total absorption of the mesh segment and add in the amount of emitted
      // energy that escapes.
      //
      FloatType absorbIntensity = static_cast<FloatType>(absorptionBins.Get(rayOffset + i));
      FloatType emissionIntensity = static_cast<FloatType>(emissionBins.Get(rayOffset + i));

      absorptionBins.Set(rayOffset + i, absorbIntensity * tmp);

      emissionIntensity = emissionIntensity * tmp + emission * (1.f - tmp);

      BOUNDS_CHECK(emissionBins, rayOffset + i);
      emissionBins.Set(rayOffset + i, emissionIntensity);
    }
    currentDistance = exitDistance;
  }
};
//
//  IdentifyMissedRay is a debugging routine that detects
//  rays that fail to have any value because of a external
//  intersection and cell intersection mismatch
//
//
class IdentifyMissedRay : public vtkm::worklet::WorkletMapField
{
public:
  vtkm::Id Width;
  vtkm::Id Height;
  vtkm::Vec<vtkm::Float32, 4> BGColor;
  IdentifyMissedRay(const vtkm::Id width,
                    const vtkm::Id height,
                    vtkm::Vec<vtkm::Float32, 4> bgcolor)
    : Width(width)
    , Height(height)
    , BGColor(bgcolor)
  {
  }
  using ControlSignature = void(FieldIn<>, WholeArrayIn<>);
  using ExecutionSignature = void(_1, _2);


  VTKM_EXEC inline bool IsBGColor(const vtkm::Vec<vtkm::Float32, 4> color) const
  {
    bool isBG = false;

    if (color[0] == BGColor[0] && color[1] == BGColor[1] && color[2] == BGColor[2] &&
        color[3] == BGColor[3])
      isBG = true;
    return isBG;
  }

  template <typename ColorBufferType>
  VTKM_EXEC inline void operator()(const vtkm::Id& pixelId, ColorBufferType& buffer) const
  {
    vtkm::Id x = pixelId % Width;
    vtkm::Id y = pixelId / Width;

    // Conservative check, we only want to check pixels in the middle
    if (x <= 0 || y <= 0)
      return;
    if (x >= Width - 1 || y >= Height - 1)
      return;
    vtkm::Vec<vtkm::Float32, 4> pixel;
    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 0));
    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 1));
    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 2));
    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(pixelId * 4 + 3));
    if (!IsBGColor(pixel))
      return;
    vtkm::Id p0 = (y)*Width + (x + 1);
    vtkm::Id p1 = (y)*Width + (x - 1);
    vtkm::Id p2 = (y + 1) * Width + (x);
    vtkm::Id p3 = (y - 1) * Width + (x);
    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 0));
    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 1));
    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 2));
    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p0 * 4 + 3));
    if (IsBGColor(pixel))
      return;
    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 0));
    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 1));
    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 2));
    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p1 * 4 + 3));
    if (IsBGColor(pixel))
      return;
    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 0));
    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 1));
    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 2));
    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p2 * 4 + 3));
    if (IsBGColor(pixel))
      return;
    pixel[0] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 0));
    pixel[1] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 1));
    pixel[2] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 2));
    pixel[3] = static_cast<vtkm::Float32>(buffer.Get(p3 * 4 + 3));
    if (IsBGColor(pixel))
      return;

    printf("Possible error ray missed ray %d\n", (int)pixelId);
  }
};

template <typename FloatType, typename Device>
class SampleCellAssocCells : public vtkm::worklet::WorkletMapField
{
private:
  using ColorHandle = typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>>;
  using ColorBuffer = typename vtkm::cont::ArrayHandle<FloatType>;
  using ColorConstPortal = typename ColorHandle::ExecutionTypes<Device>::PortalConst;
  using ColorPortal = typename ColorBuffer::template ExecutionTypes<Device>::Portal;

  CellSampler<255> Sampler;
  FloatType SampleDistance;
  FloatType MinScalar;
  FloatType InvDeltaScalar;
  ColorPortal FrameBuffer;
  ColorConstPortal ColorMap;
  const MeshConnectivityBase* MeshConn;
  vtkm::Int32 ColorMapSize;

public:
  SampleCellAssocCells(const FloatType& sampleDistance,
                       const FloatType& minScalar,
                       const FloatType& maxScalar,
                       ColorHandle& colorMap,
                       ColorBuffer& frameBuffer,
                       const MeshConnectivityBase* meshConn)
    : SampleDistance(sampleDistance)
    , MinScalar(minScalar)
    , ColorMap(colorMap.PrepareForInput(Device()))
    , MeshConn(meshConn)
  {
    InvDeltaScalar = (minScalar == maxScalar) ? 1.f : 1.f / (maxScalar - minScalar);
    ColorMapSize = static_cast<vtkm::Int32>(ColorMap.GetNumberOfValues());
    this->FrameBuffer = frameBuffer.PrepareForOutput(frameBuffer.GetNumberOfValues(), Device());
  }


  using ControlSignature = void(FieldIn<>,
                                WholeArrayIn<ScalarRenderingTypes>,
                                FieldIn<>,
                                FieldIn<>,
                                FieldInOut<>,
                                FieldInOut<>);
  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, WorkIndex);

  template <typename ScalarPortalType>
  VTKM_EXEC inline void operator()(const vtkm::Id& currentCell,
                                   ScalarPortalType& scalarPortal,
                                   const FloatType& enterDistance,
                                   const FloatType& exitDistance,
                                   FloatType& currentDistance,
                                   vtkm::UInt8& rayStatus,
                                   const vtkm::Id& pixelIndex) const
  {

    if (rayStatus != RAY_ACTIVE)
      return;

    vtkm::Vec<vtkm::Float32, 4> color;
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 0);
    color[0] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 0));
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 1);
    color[1] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 1));
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 2);
    color[2] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 2));
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 3);
    color[3] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 3));

    vtkm::Float32 scalar;
    BOUNDS_CHECK(scalarPortal, currentCell);
    scalar = vtkm::Float32(scalarPortal.Get(currentCell));
    //
    // There can be mismatches in the initial enter distance and the current distance
    // due to lost rays at cell borders. For now,
    // we will just advance the current position to the enter distance, since otherwise,
    // the pixel would never be sampled.
    //
    if (currentDistance < enterDistance)
      currentDistance = enterDistance;

    vtkm::Float32 lerpedScalar;
    lerpedScalar = static_cast<vtkm::Float32>((scalar - MinScalar) * InvDeltaScalar);
    vtkm::Id colorIndex = vtkm::Id(lerpedScalar * vtkm::Float32(ColorMapSize));
    if (colorIndex < 0)
      colorIndex = 0;
    if (colorIndex >= ColorMapSize)
      colorIndex = ColorMapSize - 1;
    BOUNDS_CHECK(ColorMap, colorIndex);
    vtkm::Vec<vtkm::Float32, 4> sampleColor = ColorMap.Get(colorIndex);

    while (enterDistance <= currentDistance && currentDistance <= exitDistance)
    {
      //composite
      vtkm::Float32 alpha = sampleColor[3] * (1.f - color[3]);
      color[0] = color[0] + sampleColor[0] * alpha;
      color[1] = color[1] + sampleColor[1] * alpha;
      color[2] = color[2] + sampleColor[2] * alpha;
      color[3] = alpha + color[3];

      if (color[3] > 1.)
      {
        rayStatus = RAY_TERMINATED;
        break;
      }
      currentDistance += SampleDistance;
    }

    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 0);
    FrameBuffer.Set(pixelIndex * 4 + 0, color[0]);
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 1);
    FrameBuffer.Set(pixelIndex * 4 + 1, color[1]);
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 2);
    FrameBuffer.Set(pixelIndex * 4 + 2, color[2]);
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 3);
    FrameBuffer.Set(pixelIndex * 4 + 3, color[3]);
  }
}; //class Sample cell

template <typename FloatType, typename Device>
class SampleCellAssocPoints : public vtkm::worklet::WorkletMapField
{
private:
  using ColorHandle = typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 4>>;
  using ColorBuffer = typename vtkm::cont::ArrayHandle<FloatType>;
  using ColorConstPortal = typename ColorHandle::ExecutionTypes<Device>::PortalConst;
  using ColorPortal = typename ColorBuffer::template ExecutionTypes<Device>::Portal;

  CellSampler<255> Sampler;
  FloatType SampleDistance;
  const MeshConnectivityBase* MeshConn;
  FloatType MinScalar;
  FloatType InvDeltaScalar;
  ColorPortal FrameBuffer;
  ColorConstPortal ColorMap;
  vtkm::Id ColorMapSize;

public:
  SampleCellAssocPoints(const FloatType& sampleDistance,
                        const FloatType& minScalar,
                        const FloatType& maxScalar,
                        ColorHandle& colorMap,
                        ColorBuffer& frameBuffer,
                        const MeshConnectivityBase* meshConn)
    : SampleDistance(sampleDistance)
    , MeshConn(meshConn)
    , MinScalar(minScalar)
    , ColorMap(colorMap.PrepareForInput(Device()))
  {
    InvDeltaScalar = (minScalar == maxScalar) ? 1.f : 1.f / (maxScalar - minScalar);
    ColorMapSize = ColorMap.GetNumberOfValues();
    this->FrameBuffer = frameBuffer.PrepareForOutput(frameBuffer.GetNumberOfValues(), Device());
  }


  using ControlSignature = void(FieldIn<>,
                                WholeArrayIn<Vec3>,
                                WholeArrayIn<ScalarRenderingTypes>,
                                FieldIn<>,
                                FieldIn<>,
                                FieldInOut<>,
                                FieldIn<>,
                                FieldInOut<>,
                                FieldIn<>);
  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7, _8, WorkIndex, _9);

  template <typename PointPortalType, typename ScalarPortalType>
  VTKM_EXEC inline void operator()(const vtkm::Id& currentCell,
                                   PointPortalType& vertices,
                                   ScalarPortalType& scalarPortal,
                                   const FloatType& enterDistance,
                                   const FloatType& exitDistance,
                                   FloatType& currentDistance,
                                   const vtkm::Vec<vtkm::Float32, 3>& dir,
                                   vtkm::UInt8& rayStatus,
                                   const vtkm::Id& pixelIndex,
                                   const vtkm::Vec<FloatType, 3>& origin) const
  {

    if (rayStatus != RAY_ACTIVE)
      return;

    vtkm::Vec<vtkm::Float32, 4> color;
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 0);
    color[0] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 0));
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 1);
    color[1] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 1));
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 2);
    color[2] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 2));
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 3);
    color[3] = static_cast<vtkm::Float32>(FrameBuffer.Get(pixelIndex * 4 + 3));

    if (color[3] >= 1.f)
    {
      rayStatus = RAY_TERMINATED;
      return;
    }
    vtkm::Vec<vtkm::Float32, 8> scalars;
    vtkm::Vec<vtkm::Vec<FloatType, 3>, 8> points;
    // silence "may" be uninitialized warning
    for (vtkm::Int32 i = 0; i < 8; ++i)
    {
      scalars[i] = 0.f;
      points[i] = vtkm::Vec<FloatType, 3>(0.f, 0.f, 0.f);
    }
    //load local scalar cell data
    vtkm::Id cellConn[8];
    const vtkm::Int32 numIndices = MeshConn->GetCellIndices(cellConn, currentCell);
    for (int i = 0; i < numIndices; ++i)
    {
      BOUNDS_CHECK(scalarPortal, cellConn[i]);
      scalars[i] = static_cast<vtkm::Float32>(scalarPortal.Get(cellConn[i]));
      BOUNDS_CHECK(vertices, cellConn[i]);
      points[i] = vtkm::Vec<FloatType, 3>(vertices.Get(cellConn[i]));
    }
    //
    // There can be mismatches in the initial enter distance and the current distance
    // due to lost rays at cell borders. For now,
    // we will just advance the current position to the enter distance, since otherwise,
    // the pixel would never be sampled.
    //
    if (currentDistance < enterDistance)
    {
      currentDistance = enterDistance;
    }

    const vtkm::Int32 cellShape = MeshConn->GetCellShape(currentCell);
    while (enterDistance <= currentDistance && currentDistance <= exitDistance)
    {
      vtkm::Vec<FloatType, 3> sampleLoc = origin + currentDistance * dir;
      vtkm::Float32 lerpedScalar;
      bool validSample =
        Sampler.SampleCell(points, scalars, sampleLoc, lerpedScalar, *this, cellShape);
      if (!validSample)
      {
        //
        // There is a slight mismatch between intersections and parametric coordinates
        // which results in a invalid sample very close to the cell edge. Just throw
        // this sample away, and move to the next sample.
        //

        //There should be a sample here, so offset and try again.

        currentDistance += 0.00001f;
        continue;
      }
      lerpedScalar = static_cast<vtkm::Float32>((lerpedScalar - MinScalar) * InvDeltaScalar);
      vtkm::Id colorIndex = vtkm::Id(lerpedScalar * vtkm::Float32(ColorMapSize));

      colorIndex = vtkm::Min(vtkm::Max(colorIndex, vtkm::Id(0)), ColorMapSize - 1);
      BOUNDS_CHECK(ColorMap, colorIndex);
      vtkm::Vec<vtkm::Float32, 4> sampleColor = ColorMap.Get(colorIndex);
      //composite
      sampleColor[3] *= (1.f - color[3]);
      color[0] = color[0] + sampleColor[0] * sampleColor[3];
      color[1] = color[1] + sampleColor[1] * sampleColor[3];
      color[2] = color[2] + sampleColor[2] * sampleColor[3];
      color[3] = sampleColor[3] + color[3];

      if (color[3] >= 1.0)
      {
        rayStatus = RAY_TERMINATED;
        break;
      }
      currentDistance += SampleDistance;
    }

    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 0);
    FrameBuffer.Set(pixelIndex * 4 + 0, color[0]);
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 1);
    FrameBuffer.Set(pixelIndex * 4 + 1, color[1]);
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 2);
    FrameBuffer.Set(pixelIndex * 4 + 2, color[2]);
    BOUNDS_CHECK(FrameBuffer, pixelIndex * 4 + 3);
    FrameBuffer.Set(pixelIndex * 4 + 3, color[3]);
  }
}; //class Sample cell

template <typename FloatType, typename Device>
void ConnectivityTracer::IntersectCell(Ray<FloatType>& rays,
                                       detail::RayTracking<FloatType>& tracker,
                                       const MeshConnectivityBase* meshConn,
                                       Device)
{
  using LocateC = LocateCell<FloatType>;
  vtkm::cont::Timer<Device> timer;
  vtkm::worklet::DispatcherMapField<LocateC> locateDispatch(meshConn);
  locateDispatch.SetDevice(Device());
  locateDispatch.Invoke(rays.HitIdx,
                        this->Coords,
                        rays.Dir,
                        *(tracker.EnterDist),
                        *(tracker.ExitDist),
                        tracker.ExitFace,
                        rays.Status,
                        rays.Origin);

  if (this->CountRayStatus)
    RaysLost = RayOperations::GetStatusCount(rays, RAY_LOST, Device());
  this->IntersectTime += timer.GetElapsedTime();
}

template <typename FloatType, typename Device>
void ConnectivityTracer::AccumulatePathLengths(Ray<FloatType>& rays,
                                               detail::RayTracking<FloatType>& tracker,
                                               Device)
{
  vtkm::worklet::DispatcherMapField<AddPathLengths<FloatType>> dispatcher(
    (AddPathLengths<FloatType>()));
  dispatcher.SetDevice(Device());
  dispatcher.Invoke(
    rays.Status, *(tracker.EnterDist), *(tracker.ExitDist), rays.GetBuffer("path_lengths").Buffer);
}

template <typename FloatType, typename Device>
void ConnectivityTracer::FindLostRays(Ray<FloatType>& rays,
                                      detail::RayTracking<FloatType>& tracker,
                                      const MeshConnectivityBase* meshConn,
                                      Device)
{
  using RayB = RayBumper<FloatType, Device>;
  vtkm::cont::Timer<Device> timer;

  vtkm::worklet::DispatcherMapField<RayB> bumpDispatch(
    RayB(rays.DirX, rays.DirY, rays.DirZ, meshConn));
  bumpDispatch.SetDevice(Device());
  bumpDispatch.Invoke(rays.HitIdx,
                      this->Coords,
                      *(tracker.EnterDist),
                      *(tracker.ExitDist),
                      tracker.ExitFace,
                      rays.Status,
                      rays.Origin);

  this->LostRayTime += timer.GetElapsedTime();
}

template <typename FloatType, typename Device>
void ConnectivityTracer::SampleCells(Ray<FloatType>& rays,
                                     detail::RayTracking<FloatType>& tracker,
                                     const MeshConnectivityBase* meshConn,
                                     Device)
{
  using SampleP = SampleCellAssocPoints<FloatType, Device>;
  using SampleC = SampleCellAssocCells<FloatType, Device>;
  vtkm::cont::Timer<Device> timer;

  VTKM_ASSERT(rays.Buffers.at(0).GetNumChannels() == 4);

  if (FieldAssocPoints)
  {
    vtkm::worklet::DispatcherMapField<SampleP> dispatcher(
      SampleP(this->SampleDistance,
              vtkm::Float32(this->ScalarBounds.Min),
              vtkm::Float32(this->ScalarBounds.Max),
              this->ColorMap,
              rays.Buffers.at(0).Buffer,
              meshConn));
    dispatcher.SetDevice(Device());
    dispatcher.Invoke(rays.HitIdx,
                      this->Coords,
                      this->ScalarField.GetData(),
                      *(tracker.EnterDist),
                      *(tracker.ExitDist),
                      tracker.CurrentDistance,
                      rays.Dir,
                      rays.Status,
                      rays.Origin);
  }
  else
  {
    vtkm::worklet::DispatcherMapField<SampleC> dispatcher(
      SampleC(this->SampleDistance,
              vtkm::Float32(this->ScalarBounds.Min),
              vtkm::Float32(this->ScalarBounds.Max),
              this->ColorMap,
              rays.Buffers.at(0).Buffer,
              meshConn));
    dispatcher.SetDevice(Device());
    dispatcher.Invoke(rays.HitIdx,
                      this->ScalarField.GetData(),
                      *(tracker.EnterDist),
                      *(tracker.ExitDist),
                      tracker.CurrentDistance,
                      rays.Status);
  }

  this->SampleTime += timer.GetElapsedTime();
}

template <typename FloatType, typename Device>
void ConnectivityTracer::IntegrateCells(Ray<FloatType>& rays,
                                        detail::RayTracking<FloatType>& tracker,
                                        Device)
{
  vtkm::cont::Timer<Device> timer;
  if (HasEmission)
  {
    bool divideEmisByAbsorp = false;
    vtkm::cont::ArrayHandle<FloatType> absorp = rays.Buffers.at(0).Buffer;
    vtkm::cont::ArrayHandle<FloatType> emission = rays.GetBuffer("emission").Buffer;
    vtkm::worklet::DispatcherMapField<IntegrateEmission<FloatType>> dispatcher(
      IntegrateEmission<FloatType>(rays.Buffers.at(0).GetNumChannels(), divideEmisByAbsorp));
    dispatcher.SetDevice(Device());
    dispatcher.Invoke(rays.Status,
                      *(tracker.EnterDist),
                      *(tracker.ExitDist),
                      rays.Distance,
                      this->ScalarField.GetData(),
                      this->EmissionField.GetData(),
                      absorp,
                      emission,
                      rays.HitIdx);
  }
  else
  {
    vtkm::worklet::DispatcherMapField<Integrate<FloatType>> dispatcher(
      Integrate<FloatType>(rays.Buffers.at(0).GetNumChannels()));
    dispatcher.SetDevice(Device());
    dispatcher.Invoke(rays.Status,
                      *(tracker.EnterDist),
                      *(tracker.ExitDist),
                      rays.Distance,
                      this->ScalarField.GetData(),
                      rays.Buffers.at(0).Buffer,
                      rays.HitIdx);
  }

  IntegrateTime += timer.GetElapsedTime();
}

// template <vtkm::Int32 CellType>
// template <typename FloatType>
// void ConnectivityTracer<CellType>::PrintDebugRay(Ray<FloatType>& rays, vtkm::Id rayId)
// {
//   vtkm::Id index = -1;
//   for (vtkm::Id i = 0; i < rays.NumRays; ++i)
//   {
//     if (rays.PixelIdx.GetPortalControl().Get(i) == rayId)
//     {
//       index = i;
//       break;
//     }
//   }
//   if (index == -1)
//   {
//     return;
//   }

//   std::cout << "++++++++RAY " << rayId << "++++++++\n";
//   std::cout << "Status: " << (int)rays.Status.GetPortalControl().Get(index) << "\n";
//   std::cout << "HitIndex: " << rays.HitIdx.GetPortalControl().Get(index) << "\n";
//   std::cout << "Dist " << rays.Distance.GetPortalControl().Get(index) << "\n";
//   std::cout << "MinDist " << rays.MinDistance.GetPortalControl().Get(index) << "\n";
//   std::cout << "Origin " << rays.Origin.GetPortalConstControl().Get(index) << "\n";
//   std::cout << "Dir " << rays.Dir.GetPortalConstControl().Get(index) << "\n";
//   std::cout << "+++++++++++++++++++++++++\n";
// }

template <typename FloatType, typename Device>
void ConnectivityTracer::OffsetMinDistances(Ray<FloatType>& rays, Device)
{
  vtkm::worklet::DispatcherMapField<AdvanceRay<FloatType>> dispatcher(
    AdvanceRay<FloatType>(FloatType(0.001)));
  dispatcher.SetDevice(Device());
  dispatcher.Invoke(rays.Status, rays.MinDistance);
}

template <typename Device, typename FloatType>
void ConnectivityTracer::RenderOnDevice(Ray<FloatType>& rays, Device)
{
  this->RaysLost = 0;
  Logger* logger = Logger::GetInstance();
  logger->OpenLogEntry("conn_tracer");
  logger->AddLogData("device", GetDeviceString(Device()));
  this->ResetTimers();
  vtkm::cont::Timer<Device> renderTimer;

  vtkm::cont::DeviceAdapterId devId = Device();

  const MeshConnectivityBase* meshConn = MeshContainer->Construct(devId);

  bool hasPathLengths = rays.HasBuffer("path_lengths");

  vtkm::cont::Timer<Device> timer;
  this->Init(); // sets sample distance
  //
  // All Rays begin as exited to force intersection
  //
  RayOperations::ResetStatus(rays, RAY_EXITED_MESH, Device());

  detail::RayTracking<FloatType> rayTracker;

  rayTracker.Init(rays.NumRays, rays.Distance, Device());
  vtkm::Float64 time = timer.GetElapsedTime();
  logger->AddLogData("init", time);

  //CountRayStatus = true;
  bool cullMissedRays = true;
  bool workRemaining = true;
  if (this->CountRayStatus)
  {
    this->PrintRayStatus(rays, Device());
  }

  do
  {
    {
      vtkm::cont::Timer<Device> entryTimer;
      //
      // if ray misses the external face it will be marked RAY_EXITED_MESH
      //
      MeshContainer->FindEntry(rays, devId);
      MeshEntryTime += entryTimer.GetElapsedTime();
    }

    if (this->CountRayStatus)
    {
      this->PrintRayStatus(rays, Device());
    }
    if (cullMissedRays)
    {
      //TODO: if we always call this after intersection, then
      //      we could make a specialized version that only compacts
      //      hitIdx distance and status, resizing everything else.
      vtkm::cont::ArrayHandle<UInt8> activeRays;
      activeRays = RayOperations::CompactActiveRays(rays, Device());
      rayTracker.Compact(rays.Distance, activeRays, Device());
      cullMissedRays = false;
    }

    if (this->CountRayStatus)
    {
      this->PrintRayStatus(rays, Device());
    }
    // TODO: we should compact out exited rays once below a threshold
    while (RayOperations::RaysInMesh(rays, Device()))
    {
      //
      // Rays the leave the mesh will be marked as RAYEXITED_MESH
      this->IntersectCell(rays, rayTracker, meshConn, Device());
      //
      // If the ray was lost due to precision issues, we find it.
      // If it is marked RAY_ABANDONED, then something went wrong.
      //
      this->FindLostRays(rays, rayTracker, meshConn, Device());
      //
      // integrate along the ray
      //
      if (this->Integrator == Volume)
        this->SampleCells(rays, rayTracker, meshConn, Device());
      else
        this->IntegrateCells(rays, rayTracker, Device());

      if (hasPathLengths)
      {
        this->AccumulatePathLengths(rays, rayTracker, Device());
      }
      //swap enter and exit distances
      rayTracker.Swap();
      if (this->CountRayStatus)
        this->PrintRayStatus(rays, Device());
    } //for

    workRemaining = RayOperations::RaysProcessed(rays, Device()) != rays.NumRays;
    //
    // Ensure that we move the current distance forward some
    // epsilon so we don't re-enter the cell we just left.
    //
    if (workRemaining)
    {
      RayOperations::CopyDistancesToMin(rays, Device());
      this->OffsetMinDistances(rays, Device());
    }
  } while (workRemaining);

  if (rays.DebugWidth != -1 && this->Integrator == Volume)
  {

    vtkm::cont::ArrayHandleCounting<vtkm::Id> pCounter(0, 1, rays.NumRays);
    vtkm::worklet::DispatcherMapField<IdentifyMissedRay> dispatcher(
      IdentifyMissedRay(rays.DebugWidth, rays.DebugHeight, this->BackgroundColor));
    dispatcher.SetDevice(Device());
    dispatcher.Invoke(pCounter, rays.Buffers.at(0).Buffer);
  }
  vtkm::Float64 renderTime = renderTimer.GetElapsedTime();
  this->LogTimers();
  logger->AddLogData("active_pixels", rays.NumRays);
  logger->CloseLogEntry(renderTime);
} //Render
}
}
} // namespace vtkm::rendering::raytracing