/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkFixedPointVolumeRayCastMapper.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkFixedPointVolumeRayCastMapper.h"

#include "vtkCamera.h"
#include "vtkColorTransferFunction.h"
#include "vtkCommand.h"
#include "vtkDataArray.h"
#include "vtkEncodedGradientShader.h"
#include "vtkFiniteDifferenceGradientEstimator.h"
#include "vtkFixedPointRayCastImage.h"
#include "vtkFixedPointVolumeRayCastCompositeGOHelper.h"
#include "vtkFixedPointVolumeRayCastCompositeGOShadeHelper.h"
#include "vtkFixedPointVolumeRayCastCompositeHelper.h"
#include "vtkFixedPointVolumeRayCastCompositeShadeHelper.h"
#include "vtkFixedPointVolumeRayCastMIPHelper.h"
#include "vtkGraphicsFactory.h"
#include "vtkImageData.h"
#include "vtkLight.h"
#include "vtkMath.h"
#include "vtkMultiThreader.h"
#include "vtkObjectFactory.h"
#include "vtkPiecewiseFunction.h"
#include "vtkPlaneCollection.h"
#include "vtkPointData.h"
#include "vtkRayCastImageDisplayHelper.h"
#include "vtkRenderWindow.h"
#include "vtkRenderer.h"
#include "vtkSphericalDirectionEncoder.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkTimerLog.h"
#include "vtkTransform.h"
#include "vtkVolumeProperty.h"
#include "vtkVolumeRayCastSpaceLeapingImageFilter.h"

#include <cmath>
#include <exception>

vtkStandardNewMacro(vtkFixedPointVolumeRayCastMapper);
vtkCxxSetObjectMacro(vtkFixedPointVolumeRayCastMapper, RayCastImage, vtkFixedPointRayCastImage);

#define vtkVRCMultiplyPointMacro(A, B, M)                                                          \
  B[0] = A[0] * M[0] + A[1] * M[1] + A[2] * M[2] + M[3];                                           \
  B[1] = A[0] * M[4] + A[1] * M[5] + A[2] * M[6] + M[7];                                           \
  B[2] = A[0] * M[8] + A[1] * M[9] + A[2] * M[10] + M[11];                                         \
  B[3] = A[0] * M[12] + A[1] * M[13] + A[2] * M[14] + M[15];                                       \
  if (B[3] != 1.0)                                                                                 \
  {                                                                                                \
    B[0] /= B[3];                                                                                  \
    B[1] /= B[3];                                                                                  \
    B[2] /= B[3];                                                                                  \
  }

#define vtkVRCMultiplyNormalMacro(A, B, M)                                                         \
  B[0] = A[0] * M[0] + A[1] * M[4] + A[2] * M[8];                                                  \
  B[1] = A[0] * M[1] + A[1] * M[5] + A[2] * M[9];                                                  \
  B[2] = A[0] * M[2] + A[1] * M[6] + A[2] * M[10]

//------------------------------------------------------------------------------
template <class T>
void vtkFixedPointVolumeRayCastMapperComputeCS1CGradients(T* dataPtr, int dim[3], double spacing[3],
  double scalarRange[2], unsigned short** gradientNormal, unsigned char** gradientMagnitude,
  vtkDirectionEncoder* directionEncoder, int thread_id, int thread_count,
  vtkFixedPointVolumeRayCastMapper* me)
{
  int x, y, z;
  vtkIdType yinc, zinc;
  int x_start, x_limit;
  int y_start, y_limit;
  int z_start, z_limit;
  T* dptr;
  float n[3], t;
  float gvalue = 0;
  int xlow, xhigh;
  double aspect[3];
  float scale;
  unsigned short* dirPtr;
  unsigned char* magPtr;

  if (thread_id == 0)
  {
    me->InvokeEvent(vtkCommand::VolumeMapperComputeGradientsStartEvent, nullptr);
  }

  double avgSpacing = (spacing[0] + spacing[1] + spacing[2]) / 3.0;

  // adjust the aspect
  aspect[0] = spacing[0] * 2.0 / avgSpacing;
  aspect[1] = spacing[1] * 2.0 / avgSpacing;
  aspect[2] = spacing[2] * 2.0 / avgSpacing;

  // compute the increments
  yinc = static_cast<vtkIdType>(dim[0]);
  zinc = yinc * static_cast<vtkIdType>(dim[1]);

  if (scalarRange[1] - scalarRange[0])
  {
    scale = 255.0 / (0.25 * (scalarRange[1] - scalarRange[0]));
  }
  else
  {
    scale = 1.0;
  }

  x_start = 0;
  x_limit = dim[0];
  y_start = 0;
  y_limit = dim[1];
  z_start = (int)(((float)thread_id / (float)thread_count) * dim[2]);
  z_limit = (int)(((float)(thread_id + 1) / (float)thread_count) * dim[2]);

  // Do sanity checking on limits - make sure they are all within bounds
  // of the scalar input
  assert(z_start >= 0);
  assert(z_limit <= dim[2]);

  int* dxBuffer = new int[dim[0]];
  int* dyBuffer = new int[dim[0]];
  int* dzBuffer = new int[dim[0]];

  for (z = z_start; z < z_limit; z++)
  {
    unsigned short* gradientDirPtr = gradientNormal[z];
    unsigned char* gradientMagPtr = gradientMagnitude[z];

    for (y = y_start; y < y_limit; y++)
    {
      xlow = x_start;
      xhigh = x_limit;

      dirPtr = gradientDirPtr + y * yinc + xlow;
      magPtr = gradientMagPtr + y * yinc + xlow;

      // Working on dx - that is this row
      dptr = dataPtr + z * zinc + y * yinc + xlow;
      // add this into our dxBuffer
      dxBuffer[0] = *dptr;
      for (x = xlow + 1; x < xhigh; x++)
      {
        *(dxBuffer + x) = *(dptr + x - 1);
      }

      // subtract this from our dxBuffer
      for (x = xlow; x < xhigh - 1; x++)
      {
        *(dxBuffer + x) -= *(dptr + x + 1);
      }
      dxBuffer[xhigh - 1] -= dptr[xhigh - 1];

      // working on dy - need the row before and the row after.
      // first, the row before, or this row if we are at the edge
      if (y > 0)
      {
        dptr = dataPtr + z * zinc + (y - 1) * yinc + xlow;
      }
      else
      {
        dptr = dataPtr + z * zinc + y * yinc + xlow;
      }
      // add this into our dyBuffer
      for (x = xlow; x < xhigh; x++)
      {
        dyBuffer[x] = dptr[x];
      }

      // now the row after
      if (y < y_limit - 1)
      {
        dptr = dataPtr + z * zinc + (y + 1) * yinc + xlow;
      }
      else
      {
        dptr = dataPtr + z * zinc + y * yinc + xlow;
      }
      // subtract this from our dyBuffer
      for (x = xlow; x < xhigh; x++)
      {
        dyBuffer[x] -= dptr[x];
      }

      // Find the pointer for the slice before - use this if there is
      // no slice before
      if (z > 0)
      {
        dptr = dataPtr + (z - 1) * zinc + y * yinc + xlow;
      }
      else
      {
        dptr = dataPtr + z * zinc + y * yinc + xlow;
      }

      // add this into our dzBuffer
      for (x = xlow; x < xhigh; x++)
      {
        dzBuffer[x] = dptr[x];
      }

      // Find the pointer for the slice after - use this if there is
      // no slice after
      if (z < z_limit - 1)
      {
        dptr = dataPtr + (z + 1) * zinc + y * yinc + xlow;
      }
      else
      {
        dptr = dataPtr + z * zinc + y * yinc + xlow;
      }

      // add this into our dzBuffer
      for (x = xlow; x < xhigh; x++)
      {
        dzBuffer[x] -= dptr[x];
      }

      // now one more loop to generate the normals
      for (x = xlow; x < xhigh; x++)
      {
        n[0] = dxBuffer[x];
        n[1] = dyBuffer[x];
        n[2] = dzBuffer[x];

        // Take care of the aspect ratio of the data
        // Scaling in the vtkVolume is isotropic, so this is the
        // only place we have to worry about non-isotropic scaling.
        n[0] /= aspect[0];
        n[1] /= aspect[1];
        n[2] /= aspect[2];

        // Compute the gradient magnitude
        t = sqrt((double)(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]));

        // Encode this into an 8 bit value
        gvalue = t * scale;

        gvalue = (gvalue < 0.0) ? (0.0) : (gvalue);
        gvalue = (gvalue > 255.0) ? (255.0) : (gvalue);

        // Normalize the gradient direction
        if (t > 0.0)
        {
          n[0] /= t;
          n[1] /= t;
          n[2] /= t;
        }
        else
        {
          n[0] = n[1] = n[2] = 0.0;
        }

        *(magPtr++) = static_cast<unsigned char>(gvalue + 0.5);
        *(dirPtr++) = directionEncoder->GetEncodedDirection(n);
      }
    }

    if ((z / thread_count) % 8 == 7 && thread_id == 0)
    {
      double args[1];
      args[0] = static_cast<float>(z - z_start) / static_cast<float>(z_limit - z_start - 1);
      me->InvokeEvent(vtkCommand::VolumeMapperComputeGradientsProgressEvent, args);
    }
  }

  delete[] dxBuffer;
  delete[] dyBuffer;
  delete[] dzBuffer;

  if (thread_id == 0)
  {
    me->InvokeEvent(vtkCommand::VolumeMapperComputeGradientsEndEvent, nullptr);
  }
}

VTK_THREAD_RETURN_TYPE vtkFPVRCMSwitchOnDataType(void* arg)
{
  vtkFixedPointVolumeRayCastMapper* mapper;
  int thread_count;
  int thread_id;

  thread_id = ((vtkMultiThreader::ThreadInfo*)(arg))->ThreadID;
  thread_count = ((vtkMultiThreader::ThreadInfo*)(arg))->NumberOfThreads;
  mapper = (vtkFixedPointVolumeRayCastMapper*)(((vtkMultiThreader::ThreadInfo*)(arg))->UserData);

  vtkImageData* input = vtkImageData::SafeDownCast(mapper->GetInput());

  void* dataPtr = mapper->GetCurrentScalars()->GetVoidPointer(0);
  int scalarType = mapper->GetCurrentScalars()->GetDataType();

  int dim[3];
  double spacing[3];
  input->GetDimensions(dim);
  input->GetSpacing(spacing);

  // Find the scalar range
  double scalarRange[2];
  mapper->GetCurrentScalars()->GetRange(scalarRange, 0);

  if (scalarType == VTK_UNSIGNED_CHAR)
  {
    vtkFixedPointVolumeRayCastMapperComputeCS1CGradients((unsigned char*)(dataPtr), dim, spacing,
      scalarRange, mapper->GradientNormal, mapper->GradientMagnitude, mapper->DirectionEncoder,
      thread_id, thread_count, mapper);
  }
  else if (scalarType == VTK_UNSIGNED_SHORT)
  {
    vtkFixedPointVolumeRayCastMapperComputeCS1CGradients((unsigned short*)(dataPtr), dim, spacing,
      scalarRange, mapper->GradientNormal, mapper->GradientMagnitude, mapper->DirectionEncoder,
      thread_id, thread_count, mapper);
  }
  else if (scalarType == VTK_CHAR)
  {
    vtkFixedPointVolumeRayCastMapperComputeCS1CGradients((char*)(dataPtr), dim, spacing,
      scalarRange, mapper->GradientNormal, mapper->GradientMagnitude, mapper->DirectionEncoder,
      thread_id, thread_count, mapper);
  }
  else if (scalarType == VTK_SHORT)
  {
    vtkFixedPointVolumeRayCastMapperComputeCS1CGradients((short*)(dataPtr), dim, spacing,
      scalarRange, mapper->GradientNormal, mapper->GradientMagnitude, mapper->DirectionEncoder,
      thread_id, thread_count, mapper);
  }

  return VTK_THREAD_RETURN_VALUE;
}

template <class T>
void vtkFixedPointVolumeRayCastMapperComputeGradients(T* dataPtr, int dim[3], double spacing[3],
  int components, int independent, double scalarRange[4][2], unsigned short** gradientNormal,
  unsigned char** gradientMagnitude, vtkDirectionEncoder* directionEncoder,
  vtkFixedPointVolumeRayCastMapper* me)
{
  int x, y, z, c;
  vtkIdType yinc, zinc;
  int x_start, x_limit;
  int y_start, y_limit;
  int z_start, z_limit;
  T *dptr, *cdptr;
  float n[3], t;
  float gvalue = 0;
  int xlow, xhigh;
  double aspect[3];
  vtkIdType xstep, ystep, zstep;
  float scale[4];
  unsigned short *dirPtr, *cdirPtr;
  unsigned char *magPtr, *cmagPtr;

  me->InvokeEvent(vtkCommand::VolumeMapperComputeGradientsStartEvent, nullptr);

  double avgSpacing = (spacing[0] + spacing[1] + spacing[2]) / 3.0;

  // adjust the aspect
  aspect[0] = spacing[0] * 2.0 / avgSpacing;
  aspect[1] = spacing[1] * 2.0 / avgSpacing;
  aspect[2] = spacing[2] * 2.0 / avgSpacing;

  // compute the increments
  yinc = static_cast<vtkIdType>(dim[0]);
  zinc = yinc * static_cast<vtkIdType>(dim[1]);

  // Compute steps through the volume in x, y, and z
  xstep = components;
  ystep = components * yinc;
  zstep = components * zinc;

  if (!independent)
  {
    if (scalarRange[components - 1][1] - scalarRange[components - 1][0])
    {
      scale[0] = 255.0 / (0.25 * (scalarRange[components - 1][1] - scalarRange[components - 1][0]));
    }
    else
    {
      scale[0] = 0.0;
    }
  }
  else
  {
    for (c = 0; c < components; c++)
    {
      if (scalarRange[c][1] - scalarRange[c][0])
      {
        scale[c] = 255.0 / (0.25 * (scalarRange[c][1] - scalarRange[c][0]));
      }
      else
      {
        scale[c] = 1.0;
      }
    }
  }

  x_start = 0;
  x_limit = dim[0];
  y_start = 0;
  y_limit = dim[1];
  z_start = 0;
  z_limit = dim[2];

  int increment = (independent) ? (components) : (1);

  float tolerance[4];
  for (c = 0; c < components; c++)
  {
    tolerance[c] = .00001 * (scalarRange[c][1] - scalarRange[c][0]);
  }

  // Loop through all the data and compute the encoded normal and
  // gradient magnitude for each scalar location
  for (z = z_start; z < z_limit; z++)
  {
    unsigned short* gradientDirPtr = gradientNormal[z];
    unsigned char* gradientMagPtr = gradientMagnitude[z];

    for (y = y_start; y < y_limit; y++)
    {
      xlow = x_start;
      xhigh = x_limit;

      dptr = dataPtr + components * (z * zinc + y * yinc + xlow);

      dirPtr = gradientDirPtr + (y * yinc + xlow) * increment;
      magPtr = gradientMagPtr + (y * yinc + xlow) * increment;

      for (x = xlow; x < xhigh; x++)
      {
        for (c = 0; (independent && c < components) || c == 0; c++)
        {
          cdptr = dptr + ((independent) ? (c) : (components - 1));
          cdirPtr = dirPtr + ((independent) ? (c) : (0));
          cmagPtr = magPtr + ((independent) ? (c) : (0));

          // Allow up to 3 tries to find the gadient - looking out at a distance of
          // 1, 2, and 3 units.
          int foundGradient = 0;
          for (int d = 1; d <= 3 && !foundGradient; d++)
          {
            // Use a central difference method if possible,
            // otherwise use a forward or backward difference if
            // we are on the edge
            // Compute the X component
            if (x < d)
            {
              n[0] = 2.0 * ((float)*(cdptr) - (float)*(cdptr + d * xstep));
            }
            else if (x >= dim[0] - d)
            {
              n[0] = 2.0 * ((float)*(cdptr - d * xstep) - (float)*(cdptr));
            }
            else
            {
              n[0] = (float)*(cdptr - d * xstep) - (float)*(cdptr + d * xstep);
            }

            // Compute the Y component
            if (y < d)
            {
              n[1] = 2.0 * ((float)*(cdptr) - (float)*(cdptr + d * ystep));
            }
            else if (y >= dim[1] - d)
            {
              n[1] = 2.0 * ((float)*(cdptr - d * ystep) - (float)*(cdptr));
            }
            else
            {
              n[1] = (float)*(cdptr - d * ystep) - (float)*(cdptr + d * ystep);
            }

            // Compute the Z component
            if (z < d)
            {
              n[2] = 2.0 * ((float)*(cdptr) - (float)*(cdptr + d * zstep));
            }
            else if (z >= dim[2] - d)
            {
              n[2] = 2.0 * ((float)*(cdptr - d * zstep) - (float)*(cdptr));
            }
            else
            {
              n[2] = (float)*(cdptr - d * zstep) - (float)*(cdptr + d * zstep);
            }

            // Take care of the aspect ratio of the data
            // Scaling in the vtkVolume is isotropic, so this is the
            // only place we have to worry about non-isotropic scaling.
            n[0] /= d * aspect[0];
            n[1] /= d * aspect[1];
            n[2] /= d * aspect[2];

            // Compute the gradient magnitude
            t = sqrt((double)(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]));

            // Encode this into an 8 bit value
            gvalue = t * scale[c];

            if (d > 1)
            {
              gvalue = 0;
            }

            gvalue = (gvalue < 0.0) ? (0.0) : (gvalue);
            gvalue = (gvalue > 255.0) ? (255.0) : (gvalue);

            // Normalize the gradient direction
            if (t > tolerance[c])
            {
              n[0] /= t;
              n[1] /= t;
              n[2] /= t;
              foundGradient = 1;
            }
            else
            {
              n[0] = n[1] = n[2] = 0.0;
            }
          }

          *cmagPtr = static_cast<unsigned char>(gvalue + 0.5);
          *cdirPtr = directionEncoder->GetEncodedDirection(n);
        }

        dptr += components;
        dirPtr += increment;
        magPtr += increment;
      }
    }
    if (z % 8 == 7)
    {
      double args[1];
      args[0] = static_cast<float>(z - z_start) / static_cast<float>(z_limit - z_start - 1);
      me->InvokeEvent(vtkCommand::VolumeMapperComputeGradientsProgressEvent, args);
    }
  }

  me->InvokeEvent(vtkCommand::VolumeMapperComputeGradientsEndEvent, nullptr);
}

//------------------------------------------------------------------------------
// Construct a new vtkFixedPointVolumeRayCastMapper with default values
vtkFixedPointVolumeRayCastMapper::vtkFixedPointVolumeRayCastMapper()
{
  this->SampleDistance = 1.0;
  this->InteractiveSampleDistance = 2.0;
  this->ImageSampleDistance = 1.0;
  this->MinimumImageSampleDistance = 1.0;
  this->MaximumImageSampleDistance = 10.0;
  this->AutoAdjustSampleDistances = 1;
  this->LockSampleDistanceToInputSpacing = 0;

  // Should never be used without initialization, but
  // set here to avoid compiler warnings
  this->OldSampleDistance = 1.0;
  this->OldImageSampleDistance = 1.0;

  this->PerspectiveMatrix = vtkMatrix4x4::New();
  this->ViewToWorldMatrix = vtkMatrix4x4::New();
  this->ViewToVoxelsMatrix = vtkMatrix4x4::New();
  this->VoxelsToViewMatrix = vtkMatrix4x4::New();
  this->WorldToVoxelsMatrix = vtkMatrix4x4::New();
  this->VoxelsToWorldMatrix = vtkMatrix4x4::New();

  this->VolumeMatrix = vtkMatrix4x4::New();

  this->PerspectiveTransform = vtkTransform::New();
  this->VoxelsTransform = vtkTransform::New();
  this->VoxelsToViewTransform = vtkTransform::New();

  this->Threader = vtkMultiThreader::New();
  this->ThreadWarning = true;
  this->RayCastImage = vtkFixedPointRayCastImage::New();

  this->RowBounds = nullptr;
  this->OldRowBounds = nullptr;

  this->RenderTimeTable = nullptr;
  this->RenderVolumeTable = nullptr;
  this->RenderRendererTable = nullptr;
  this->RenderTableSize = 0;
  this->RenderTableEntries = 0;

  this->RenderWindow = nullptr;

  this->MIPHelper = vtkFixedPointVolumeRayCastMIPHelper::New();
  this->CompositeHelper = vtkFixedPointVolumeRayCastCompositeHelper::New();
  this->CompositeGOHelper = vtkFixedPointVolumeRayCastCompositeGOHelper::New();
  this->CompositeShadeHelper = vtkFixedPointVolumeRayCastCompositeShadeHelper::New();
  this->CompositeGOShadeHelper = vtkFixedPointVolumeRayCastCompositeGOShadeHelper::New();

  this->IntermixIntersectingGeometry = 1;

  int i;
  for (i = 0; i < 4; i++)
  {
    this->SavedRGBFunction[i] = nullptr;
    this->SavedGrayFunction[i] = nullptr;
    this->SavedScalarOpacityFunction[i] = nullptr;
    this->SavedGradientOpacityFunction[i] = nullptr;
    this->SavedColorChannels[i] = 0;
    this->SavedScalarOpacityDistance[i] = 0;
    this->TableSize[i] = 0;
  }

  this->SavedSampleDistance = 0;
  this->SavedBlendMode = -1;

  this->SavedGradientsInput = nullptr;
  this->SavedParametersInput = nullptr;

  this->NumberOfGradientSlices = 0;
  this->GradientNormal = nullptr;
  this->GradientMagnitude = nullptr;
  this->ContiguousGradientNormal = nullptr;
  this->ContiguousGradientMagnitude = nullptr;

  this->DirectionEncoder = vtkSphericalDirectionEncoder::New();
  this->GradientShader = vtkEncodedGradientShader::New();
  this->GradientEstimator = vtkFiniteDifferenceGradientEstimator::New();

  this->GradientEstimator->SetDirectionEncoder(this->DirectionEncoder);

  this->ShadingRequired = 0;
  this->GradientOpacityRequired = 0;

  this->CroppingRegionMask[0] = 1;
  for (i = 1; i < 27; i++)
  {
    this->CroppingRegionMask[i] = this->CroppingRegionMask[i - 1] * 2;
  }

  this->NumTransformedClippingPlanes = 0;
  this->TransformedClippingPlanes = nullptr;

  // Which scalar field are we rendering this time, and which
  // did we render last time (so we can check if it is changing)
  this->CurrentScalars = nullptr;
  this->PreviousScalars = nullptr;

  this->ImageDisplayHelper = vtkRayCastImageDisplayHelper::New();
  this->ImageDisplayHelper->PreMultipliedColorsOn();
  this->ImageDisplayHelper->SetPixelScale(2.0);

  // This is the min max volume used for space leaping. Each 4x4x4 cell from
  // the original input volume has three values per component - a minimum scalar
  // index, maximum scalar index, and a values used for both the maximum gradient
  // magnitude and a flag. The flag is used to indicate for the
  // current transfer function whether any non-zero opacity exists between the
  // minimum and maximum scalar values and up to the maximum gradient magnitude
  this->MinMaxVolume = nullptr;
  this->MinMaxVolumeSize[0] = 0;
  this->MinMaxVolumeSize[1] = 0;
  this->MinMaxVolumeSize[2] = 0;
  this->MinMaxVolumeSize[3] = 0;
  this->SavedMinMaxInput = nullptr;

  this->Volume = nullptr;

  this->FinalColorWindow = 1.0;
  this->FinalColorLevel = 0.5;

  this->FlipMIPComparison = 0;

  this->TableShift[0] = 0;
  this->TableShift[1] = 0;
  this->TableShift[2] = 0;
  this->TableShift[3] = 0;

  this->TableScale[0] = 1;
  this->TableScale[1] = 1;
  this->TableScale[2] = 1;
  this->TableScale[3] = 1;

  this->SpaceLeapFilter = vtkVolumeRayCastSpaceLeapingImageFilter::New();

  // Cached space leaping output. This is shared between runs. The output
  // of the last run is passed back to the SpaceLeapFilter and its reused
  // since we may not be updating every flag in this structure.
  this->MinMaxVolumeCache = vtkImageData::New();
}

//------------------------------------------------------------------------------
// Destruct a vtkFixedPointVolumeRayCastMapper - clean up any memory used
vtkFixedPointVolumeRayCastMapper::~vtkFixedPointVolumeRayCastMapper()
{
  this->SpaceLeapFilter->Delete();

  this->PerspectiveMatrix->Delete();
  this->ViewToWorldMatrix->Delete();
  this->ViewToVoxelsMatrix->Delete();
  this->VoxelsToViewMatrix->Delete();
  this->WorldToVoxelsMatrix->Delete();
  this->VoxelsToWorldMatrix->Delete();

  this->VolumeMatrix->Delete();

  this->VoxelsTransform->Delete();
  this->VoxelsToViewTransform->Delete();
  this->PerspectiveTransform->Delete();

  this->Threader->Delete();

  this->MIPHelper->Delete();
  this->CompositeHelper->Delete();
  this->CompositeGOHelper->Delete();
  this->CompositeShadeHelper->Delete();
  this->CompositeGOShadeHelper->Delete();

  if (this->RayCastImage)
  {
    this->RayCastImage->Delete();
    this->RayCastImage = nullptr;
  }

  delete[] this->RenderTimeTable;
  delete[] this->RenderVolumeTable;
  delete[] this->RenderRendererTable;

  delete[] this->RowBounds;
  delete[] this->OldRowBounds;

  int i;
  if (this->GradientNormal)
  {
    // Contiguous? Delete in one chunk otherwise delete slice by slice
    if (this->ContiguousGradientNormal)
    {
      delete[] this->ContiguousGradientNormal;
      this->ContiguousGradientNormal = nullptr;
    }
    else
    {
      for (i = 0; i < this->NumberOfGradientSlices; i++)
      {
        delete[] this->GradientNormal[i];
      }
    }
    delete[] this->GradientNormal;
    this->GradientNormal = nullptr;
  }

  if (this->GradientMagnitude)
  {
    // Contiguous? Delete in one chunk otherwise delete slice by slice
    if (this->ContiguousGradientMagnitude)
    {
      delete[] this->ContiguousGradientMagnitude;
      this->ContiguousGradientMagnitude = nullptr;
    }
    else
    {
      for (i = 0; i < this->NumberOfGradientSlices; i++)
      {
        delete[] this->GradientMagnitude[i];
      }
    }
    delete[] this->GradientMagnitude;
    this->GradientMagnitude = nullptr;
  }

  this->DirectionEncoder->Delete();
  this->GradientShader->Delete();
  this->GradientEstimator->Delete();

  delete[] this->TransformedClippingPlanes;

  this->ImageDisplayHelper->Delete();

  this->MinMaxVolumeCache->Delete();
}

float vtkFixedPointVolumeRayCastMapper::ComputeRequiredImageSampleDistance(
  float desiredTime, vtkRenderer* ren)
{
  return this->ComputeRequiredImageSampleDistance(desiredTime, ren, nullptr);
}

float vtkFixedPointVolumeRayCastMapper::ComputeRequiredImageSampleDistance(
  float desiredTime, vtkRenderer* ren, vtkVolume* vol)
{
  float result;

  float oldTime;

  if (vol)
  {
    oldTime = this->RetrieveRenderTime(ren, vol);
  }
  else
  {
    oldTime = this->RetrieveRenderTime(ren);
  }

  float newTime = desiredTime;

  if (oldTime == 0.0)
  {
    if (newTime > 10)
    {
      result = this->MinimumImageSampleDistance;
    }
    else
    {
      result = this->MaximumImageSampleDistance / 2.0;
    }
  }
  else
  {
    oldTime /= (this->ImageSampleDistance * this->ImageSampleDistance);
    result = this->ImageSampleDistance * sqrt(oldTime / newTime);
    result =
      (result > this->MaximumImageSampleDistance) ? (this->MaximumImageSampleDistance) : (result);
    result =
      (result < this->MinimumImageSampleDistance) ? (this->MinimumImageSampleDistance) : (result);
  }

  return result;
}

float vtkFixedPointVolumeRayCastMapper::RetrieveRenderTime(vtkRenderer* ren, vtkVolume* vol)
{
  int i;

  for (i = 0; i < this->RenderTableEntries; i++)
  {
    if (this->RenderVolumeTable[i] == vol && this->RenderRendererTable[i] == ren)
    {
      return this->RenderTimeTable[i];
    }
  }

  return 0.0;
}

float vtkFixedPointVolumeRayCastMapper::RetrieveRenderTime(vtkRenderer* ren)
{
  int i;

  for (i = 0; i < this->RenderTableEntries; i++)
  {
    if (this->RenderRendererTable[i] == ren)
    {
      return this->RenderTimeTable[i];
    }
  }

  return 0.0;
}

void vtkFixedPointVolumeRayCastMapper::StoreRenderTime(vtkRenderer* ren, vtkVolume* vol, float time)
{
  int i;
  for (i = 0; i < this->RenderTableEntries; i++)
  {
    if (this->RenderVolumeTable[i] == vol && this->RenderRendererTable[i] == ren)
    {
      this->RenderTimeTable[i] = time;
      return;
    }
  }

  // Need to increase size
  if (this->RenderTableEntries >= this->RenderTableSize)
  {
    if (this->RenderTableSize == 0)
    {
      this->RenderTableSize = 10;
    }
    else
    {
      this->RenderTableSize *= 2;
    }

    float* oldTimePtr = this->RenderTimeTable;
    vtkVolume** oldVolumePtr = this->RenderVolumeTable;
    vtkRenderer** oldRendererPtr = this->RenderRendererTable;

    this->RenderTimeTable = new float[this->RenderTableSize];
    this->RenderVolumeTable = new vtkVolume*[this->RenderTableSize];
    this->RenderRendererTable = new vtkRenderer*[this->RenderTableSize];

    for (i = 0; i < this->RenderTableEntries; i++)
    {
      this->RenderTimeTable[i] = oldTimePtr[i];
      this->RenderVolumeTable[i] = oldVolumePtr[i];
      this->RenderRendererTable[i] = oldRendererPtr[i];
    }

    delete[] oldTimePtr;
    delete[] oldVolumePtr;
    delete[] oldRendererPtr;
  }

  this->RenderTimeTable[this->RenderTableEntries] = time;
  this->RenderVolumeTable[this->RenderTableEntries] = vol;
  this->RenderRendererTable[this->RenderTableEntries] = ren;

  this->RenderTableEntries++;
}

void vtkFixedPointVolumeRayCastMapper::SetNumberOfThreads(int num)
{
  this->Threader->SetNumberOfThreads(num);
}

int vtkFixedPointVolumeRayCastMapper::GetNumberOfThreads()
{
  if (this->Threader)
  {
    return this->Threader->GetNumberOfThreads();
  }
  return 0;
}

//------------------------------------------------------------------------------
// This method should be called after UpdateColorTables since it
// relies on some information (shift and scale) computed in that method,
// as well as the last built time for the color tables.
void vtkFixedPointVolumeRayCastMapper::UpdateMinMaxVolume(vtkVolume* vol)
{
  // A three bit variable:
  //   first bit indicates need to update flags
  //   second bit indicates need to update scalars
  //   third bit indicates need to update gradient magnitudes
  int needToUpdate = 0;

  // Get the image data
  vtkImageData* input = vtkImageData::SafeDownCast(this->GetInput());

  // We'll need this info later
  int dim[3];
  input->GetDimensions(dim);

  // Has the data itself changed?
  if (input != this->SavedMinMaxInput ||
    input->GetMTime() > this->SpaceLeapFilter->GetLastMinMaxBuildTime() ||
    this->CurrentScalars != this->PreviousScalars)
  {
    needToUpdate |= 0x03;
  }

  // Do the gradient magnitudes need to be filled in?
  if (this->GradientOpacityRequired && // done
    (needToUpdate & 0x02 ||            // done
      this->SavedGradientsMTime.GetMTime() > this->SpaceLeapFilter->GetLastMinMaxBuildTime()))
  {
    needToUpdate |= 0x05;
  }

  // Have the parameters changed which means the flags need
  // to be recomputed. Actually, we could be checking just
  // a subset of these parameters (we don't need to recompute
  // the flags if the colors change, but unless these seems
  // like a significant performance problem, I'd rather not
  // complicate the code)
  if (!(needToUpdate & 0x01) &&
    this->SavedParametersMTime.GetMTime() > this->SpaceLeapFilter->GetLastMinMaxFlagTime())
  {
    needToUpdate |= 0x01;
  }

  if (!needToUpdate)
  {
    return;
  }

  // Set the update flags, telling the filter what to update...
  this->SpaceLeapFilter->SetInputConnection(this->GetInputConnection(0, 0));
  this->SpaceLeapFilter->SetCurrentScalars(this->CurrentScalars);
  this->SpaceLeapFilter->SetIndependentComponents(vol->GetProperty()->GetIndependentComponents());
  this->SpaceLeapFilter->SetComputeMinMax((needToUpdate & 0x02) ? 1 : 0);
  this->SpaceLeapFilter->SetComputeGradientOpacity((needToUpdate & 0x04) ? 1 : 0);
  this->SpaceLeapFilter->SetUpdateGradientOpacityFlags(
    (this->GradientOpacityRequired && (needToUpdate & 0x01)) ? 1 : 0);
  this->SpaceLeapFilter->SetGradientMagnitude(this->GradientMagnitude);
  this->SpaceLeapFilter->SetTableSize(this->TableSize);
  this->SpaceLeapFilter->SetTableShift(this->TableShift);
  this->SpaceLeapFilter->SetTableScale(this->TableScale);
  for (unsigned int compIdx = 0; compIdx < 4; ++compIdx)
  {
    this->SpaceLeapFilter->SetScalarOpacityTable(compIdx, this->ScalarOpacityTable[compIdx]);
    this->SpaceLeapFilter->SetGradientOpacityTable(compIdx, this->GradientOpacityTable[compIdx]);
  }
  this->SpaceLeapFilter->SetCache(this->MinMaxVolumeCache);
  this->SpaceLeapFilter->Update();
  this->MinMaxVolume = this->SpaceLeapFilter->GetMinMaxVolume(this->MinMaxVolumeSize);

  // Cached space leaping output. This is shared between runs. The output
  // of the last run is passed back to the SpaceLeapFilter and its reused
  // since we may not be updating every flag in this structure.
  this->MinMaxVolumeCache->ShallowCopy(this->SpaceLeapFilter->GetOutput());

  // For debugging if necessary, write the min-max-volume components.
  // vtkVolumeRayCastSpaceLeapingImageFilter::WriteMinMaxVolume(
  //  0,this->MinMaxVolume,this->MinMaxVolumeSize,
  //  "MinMaxVolumeNewComponent0.mha");

  // If the line below is commented out, we get reference counting loops
  this->SpaceLeapFilter->SetInputConnection(nullptr);

  if (needToUpdate & 0x02)
  {
    this->SavedMinMaxInput = input;
  }
}

//------------------------------------------------------------------------------
void vtkFixedPointVolumeRayCastMapper::UpdateCroppingRegions()
{
  this->ConvertCroppingRegionPlanesToVoxels();

  int i;
  for (i = 0; i < 6; i++)
  {
    this->FixedPointCroppingRegionPlanes[i] =
      this->ToFixedPointPosition(this->VoxelCroppingRegionPlanes[i]);
  }
}

// This is the initialization that should be done once per image
// The render has been broken into several parts to support AMR
// volume rendering. Basically, this is done by having the AMR
// mapper call the PerImageInitialization once, then the
// PerVolumeInitialization once for each volume in the hierarchical
// structure. Finally, the AMR mapper divides all the volumes
// into subvolumes in order to render everything in a back-to-front
// order. The PerSubVolumeInitialization is called for each subvolume,
// then the RenderSubVolume is called. Finally, the DisplayImage method
// is called to map the image onto the screen. When this class is used
// directly as the mapper, the Render method calls these initialization
// methods and the RenderSubVolumeMethod. The AMR mapper will set the
// multiRender flag to 1 indicating that the PerImageInitialization
// should fully polulate the RayCastImage class based on the
// origin, spacing, and extent passed in. This will result in computing
// some things twice - once for the "full" volume (the extent bounding
// all volumes in the hierarchy), then once for each volume in the
// hierarchy. This does not make sense when rendering just a single
// volume so the multiRender flag indicates whether to do this
// computation here or skip it for later.
int vtkFixedPointVolumeRayCastMapper::PerImageInitialization(vtkRenderer* ren, vtkVolume* vol,
  int multiRender, double inputOrigin[3], double inputSpacing[3], int inputExtent[6])
{
  // Save this so that we can restore it if the image is cancelled
  this->OldImageSampleDistance = this->ImageSampleDistance;
  this->OldSampleDistance = this->SampleDistance;

  // If we are automatically adjusting the size to achieve a desired frame
  // rate, then do that adjustment here. Base the new image sample distance
  // on the previous one and the previous render time. Don't let
  // the adjusted image sample distance be less than the minimum image sample
  // distance or more than the maximum image sample distance.
  if (this->AutoAdjustSampleDistances)
  {
    this->ImageSampleDistance =
      this->ComputeRequiredImageSampleDistance(vol->GetAllocatedRenderTime(), ren, vol);

    // If this is an interactive render (faster than 1 frame per second) then we'll
    // increase the sample distance along the ray to improve performance
    if (vol->GetAllocatedRenderTime() < 1.0)
    {
      this->SampleDistance = this->InteractiveSampleDistance;
    }
  }

  // Pass the ImageSampleDistance on the RayCastImage
  this->RayCastImage->SetImageSampleDistance(this->ImageSampleDistance);

  // The full image fills the viewport. First, compute the actual viewport
  // size, then divide by the ImageSampleDistance to find the full image
  // size in pixels
  int width, height;
  ren->GetTiledSize(&width, &height);
  this->RayCastImage->SetImageViewportSize(static_cast<int>(width / this->ImageSampleDistance),
    static_cast<int>(height / this->ImageSampleDistance));

  if (multiRender)
  {
    this->UpdateCroppingRegions();
    this->ComputeMatrices(inputOrigin, inputSpacing, inputExtent, ren, vol);

    if (!this->ComputeRowBounds(ren, 1, 0, inputExtent))
    {
      return 0;
    }
  }

  return 1;
}

// This is the initialization that should be done once per volume
void vtkFixedPointVolumeRayCastMapper::PerVolumeInitialization(vtkRenderer* ren, vtkVolume* vol)
{
  // This is the input of this mapper
  vtkImageData* input = vtkImageData::SafeDownCast(this->GetInput());
  this->PreviousScalars = this->CurrentScalars;

  // make sure that we have scalar input and update the scalar input
  if (input == nullptr)
  {
    vtkErrorMacro(<< "No Input!");
    return;
  }
  else
  {
    vtkAlgorithm* inAlg = this->GetInputAlgorithm();
    inAlg->UpdateWholeExtent();
  }

  int usingCellColors;
  this->CurrentScalars = vtkFixedPointVolumeRayCastMapper::GetScalars(input, this->ScalarMode,
    this->ArrayAccessMode, this->ArrayId, this->ArrayName, usingCellColors);

  if (usingCellColors)
  {
    vtkErrorMacro("Cell Scalars not supported");
    return;
  }

  // Compute some matrices from voxels to view and vice versa based
  // on the whole input
  double inputSpacing[3];
  double inputOrigin[3];
  int inputExtent[6];
  input->GetSpacing(inputSpacing);
  input->GetOrigin(inputOrigin);
  input->GetExtent(inputExtent);

  this->ComputeMatrices(inputOrigin, inputSpacing, inputExtent, ren, vol);

  this->RenderWindow = ren->GetRenderWindow();
  this->Volume = vol;

  // Adjust the sample spacing if necessary
  if (this->LockSampleDistanceToInputSpacing)
  {
    double const dist = this->SpacingAdjustedSampleDistance(inputSpacing, inputExtent);

    // Need to treat interactive renders differently, because if
    // AutoAdjustSampleDistances is on, then we doubled the sample
    // distance
    if (this->AutoAdjustSampleDistances && vol->GetAllocatedRenderTime() < 1.0)
    {
      if (this->SampleDistance / (dist * 2) < 0.999 || this->SampleDistance / (dist * 2) > 1.001)
      {
        this->OldSampleDistance = dist;
        this->SampleDistance = dist * 2;
        this->InteractiveSampleDistance = dist * 2;
      }
    }
    else if (this->SampleDistance / dist < 0.999 || this->SampleDistance / dist > 1.001)
    {
      this->OldSampleDistance = dist;
      this->SampleDistance = dist;
      this->InteractiveSampleDistance = dist * 2;
    }
  }

  this->UpdateColorTable(vol);
  this->UpdateGradients(vol);
  this->UpdateShadingTable(ren, vol);

  // Calls SpaceLeapFilter for a multi-threaded optimized computation of the
  // min-max structure.
  this->UpdateMinMaxVolume(vol);
}

// This is the initialization that should be done once per subvolume
void vtkFixedPointVolumeRayCastMapper::PerSubVolumeInitialization(
  vtkRenderer* ren, vtkVolume* vol, int multiRender)
{
  this->UpdateCroppingRegions();

  // Compute row bounds. This will also compute the size of the image to
  // render, allocate the space if necessary, and clear the image where
  // required. If no rays need to be cast, restore the old image sample
  // distance and return
  int inputExtent[6];
  vtkImageData* input = vtkImageData::SafeDownCast(this->GetInput());
  input->GetExtent(inputExtent);

  // If this is part of a multirender (AMR volume rendering) then
  // the image parameters have already been computed and we can skip
  // that. In all cases we need to compute the row bounds so pass in
  // a 1 for that flag
  int imageFlag = (multiRender) ? (0) : (1);
  if (!this->ComputeRowBounds(ren, imageFlag, 1, inputExtent))
  {
    this->AbortRender();
    return;
  }

  // If this is part of a multiRender, then we've already captured the z buffer,
  // otherwise we need to do it here
  if (!multiRender)
  {
    this->CaptureZBuffer(ren);
  }

  this->InitializeRayInfo(vol);
}

// This is the render method for the subvolume
void vtkFixedPointVolumeRayCastMapper::RenderSubVolume()
{
  // Set the number of threads to use for ray casting,
  // then set the execution method and do it.
  this->InvokeEvent(vtkCommand::VolumeMapperRenderStartEvent, nullptr);
  this->Threader->SetSingleMethod(FixedPointVolumeRayCastMapper_CastRays, (void*)this);
  this->Threader->SingleMethodExecute();
  this->InvokeEvent(vtkCommand::VolumeMapperRenderEndEvent, nullptr);
}

// This method displays the image that has been created
void vtkFixedPointVolumeRayCastMapper::DisplayRenderedImage(vtkRenderer* ren, vtkVolume* vol)
{
  float depth;
  if (this->IntermixIntersectingGeometry)
  {
    depth = this->MinimumViewDistance;
  }
  else
  {
    depth = -1;
  }

  if (this->FinalColorWindow != 1.0 || this->FinalColorLevel != 0.5)
  {
    this->ApplyFinalColorWindowLevel();
  }

  this->ImageDisplayHelper->RenderTexture(vol, ren, this->RayCastImage, depth);
}

void vtkFixedPointVolumeRayCastMapper::ApplyFinalColorWindowLevel()
{
  double scale = 1.0 / this->FinalColorWindow;
  double bias = 0.5 - this->FinalColorLevel / this->FinalColorWindow;

  unsigned short* image = this->RayCastImage->GetImage();
  unsigned short* iptr;

  int fullSize[2];
  this->RayCastImage->GetImageMemorySize(fullSize);

  int size[2];
  this->RayCastImage->GetImageInUseSize(size);

  int i, j;

  for (j = 0; j < fullSize[1]; j++)
  {
    iptr = image + 4 * j * fullSize[0];
    for (i = 0; i < size[0]; i++)
    {
      int tmp;

      // Red component
      tmp = (int)((float)(*iptr) * scale + bias * (float)(*(iptr + 3)));
      tmp = (tmp < 0) ? (0) : (tmp);
      tmp = (tmp > 32767) ? (32767) : (tmp);
      *iptr = tmp;

      // Green component
      iptr++;
      tmp = (int)((float)(*iptr) * scale + bias * (float)(*(iptr + 2)));
      tmp = (tmp < 0) ? (0) : (tmp);
      tmp = (tmp > 32767) ? (32767) : (tmp);
      *iptr = tmp;

      // Green component
      iptr++;
      tmp = (int)((float)(*iptr) * scale + bias * (float)(*(iptr + 1)));
      tmp = (tmp < 0) ? (0) : (tmp);
      tmp = (tmp > 32767) ? (32767) : (tmp);
      *iptr = tmp;

      // alpha - do nothing
      iptr++;
      iptr++;
    }
  }
}

// This method should be called when the render is aborted to restore previous values.
// Otherwise, the old time is still stored, with the newly computed sample distances,
// and that will cause problems on the next render.
void vtkFixedPointVolumeRayCastMapper::AbortRender()
{
  // Restore values
  this->ImageSampleDistance = this->OldImageSampleDistance;
  this->SampleDistance = this->OldSampleDistance;
}

// Capture the ZBuffer to use for intermixing with opaque geometry
// that has already been rendered
void vtkFixedPointVolumeRayCastMapper::CaptureZBuffer(vtkRenderer* ren)
{
  // How big is the viewport in pixels?
  double* viewport = ren->GetViewport();
  const int* renWinSize = ren->GetRenderWindow()->GetSize();

  // Do we need to capture the z buffer to intermix intersecting
  // geometry? If so, do it here
  if (this->IntermixIntersectingGeometry && ren->GetNumberOfPropsRendered())
  {
    int x1, x2, y1, y2;

    // turn ImageOrigin into (x1,y1) in window (not viewport!)
    // coordinates.
    int imageOrigin[2];
    int imageInUseSize[2];
    this->RayCastImage->GetImageOrigin(imageOrigin);
    this->RayCastImage->GetImageInUseSize(imageInUseSize);

    x1 = static_cast<int>(viewport[0] * static_cast<float>(renWinSize[0]) +
      static_cast<float>(imageOrigin[0]) * this->ImageSampleDistance);
    y1 = static_cast<int>(viewport[1] * static_cast<float>(renWinSize[1]) +
      static_cast<float>(imageOrigin[1]) * this->ImageSampleDistance);

    int zbufferSize[2];
    int zbufferOrigin[2];

    // compute z buffer size
    zbufferSize[0] =
      static_cast<int>(static_cast<float>(imageInUseSize[0]) * this->ImageSampleDistance);
    zbufferSize[1] =
      static_cast<int>(static_cast<float>(imageInUseSize[1]) * this->ImageSampleDistance);

    // Use the size to compute (x2,y2) in window coordinates
    x2 = x1 + zbufferSize[0] - 1;
    y2 = y1 + zbufferSize[1] - 1;

    // This is the z buffer origin (in viewport coordinates)
    zbufferOrigin[0] =
      static_cast<int>(static_cast<float>(imageOrigin[0]) * this->ImageSampleDistance);
    zbufferOrigin[1] =
      static_cast<int>(static_cast<float>(imageOrigin[1]) * this->ImageSampleDistance);

    this->RayCastImage->SetZBufferSize(zbufferSize);
    this->RayCastImage->SetZBufferOrigin(zbufferOrigin);
    this->RayCastImage->AllocateZBuffer();

    // Capture the z buffer
    ren->GetRenderWindow()->GetZbufferData(x1, y1, x2, y2, this->RayCastImage->GetZBuffer());

    this->RayCastImage->UseZBufferOn();
  }
  else
  {
    this->RayCastImage->UseZBufferOff();
  }
}

void vtkFixedPointVolumeRayCastMapper::Render(vtkRenderer* ren, vtkVolume* vol)
{
  //  if (this->Threader->GetNumberOfThreads() > 1 && this->ThreadWarning)
  //    {
  //    vtkWarningMacro("Number of threads "
  //                    << this->Threader->GetNumberOfThreads()
  //                    << " is  > 1. This class does not produce repeatable results when the number
  //                    of threads exceeds 1.");
  //    this->ThreadWarning = false;
  //    }

  if (vtkImageData::SafeDownCast(this->GetInput()) == nullptr)
  {
    vtkWarningMacro("Mapper supports only vtkImageData");
    return;
  }

  if (this->GetBlendMode() != vtkVolumeMapper::COMPOSITE_BLEND &&
    this->GetBlendMode() != vtkVolumeMapper::MAXIMUM_INTENSITY_BLEND &&
    this->GetBlendMode() != vtkVolumeMapper::MINIMUM_INTENSITY_BLEND &&
    this->GetBlendMode() != vtkVolumeMapper::ADDITIVE_BLEND)
  {
    vtkErrorMacro(<< "Selected blend mode not supported. "
                  << "Only Composite, MIP, MinIP and additive modes "
                  << "are supported by the fixed point implementation.");
    return;
  }

  this->Timer->StartTimer();

  // Since we are passing in a value of 0 for the multiRender flag
  // (this is a single render pass - not part of a multipass AMR render)
  // then we know the origin, spacing, and extent values will not
  // be used so just initialize everything to 0. No need to check
  // the return value of the PerImageInitialization method - since this
  // is not a multirender it will always return 1.
  double dummyOrigin[3] = { 0.0, 0.0, 0.0 };
  double dummySpacing[3] = { 0.0, 0.0, 0.0 };
  int dummyExtent[6] = { 0, 0, 0, 0, 0, 0 };
  this->PerImageInitialization(ren, vol, 0, dummyOrigin, dummySpacing, dummyExtent);

  this->PerVolumeInitialization(ren, vol);

  vtkRenderWindow* renWin = ren->GetRenderWindow();

  if (renWin && renWin->CheckAbortStatus())
  {
    this->AbortRender();
    return;
  }

  this->PerSubVolumeInitialization(ren, vol, 0);
  if (renWin && renWin->CheckAbortStatus())
  {
    this->AbortRender();
    return;
  }

  this->RenderSubVolume();

  if (renWin && renWin->CheckAbortStatus())
  {
    this->AbortRender();
    return;
  }

  this->DisplayRenderedImage(ren, vol);

  this->Timer->StopTimer();
  this->TimeToDraw = this->Timer->GetElapsedTime();
  // If we've increased the sample distance, account for that in the stored time. Since we
  // don't get linear performance improvement, use a factor of .66
  this->StoreRenderTime(ren, vol,
    this->TimeToDraw * this->ImageSampleDistance * this->ImageSampleDistance *
      (1.0 + 0.66 * (this->SampleDistance - this->OldSampleDistance) / this->OldSampleDistance));

  this->SampleDistance = this->OldSampleDistance;
}

VTK_THREAD_RETURN_TYPE FixedPointVolumeRayCastMapper_CastRays(void* arg)
{
  // Get the info out of the input structure
  int threadID = ((vtkMultiThreader::ThreadInfo*)(arg))->ThreadID;
  int threadCount = ((vtkMultiThreader::ThreadInfo*)(arg))->NumberOfThreads;

  vtkFixedPointVolumeRayCastMapper* me =
    (vtkFixedPointVolumeRayCastMapper*)(((vtkMultiThreader::ThreadInfo*)arg)->UserData);

  if (!me)
  {
    vtkGenericWarningMacro("Irrecoverable error: no mapper specified");
    return VTK_THREAD_RETURN_VALUE;
  }

  vtkVolume* vol = me->GetVolume();

  if (me->GetBlendMode() == vtkVolumeMapper::MAXIMUM_INTENSITY_BLEND ||
    me->GetBlendMode() == vtkVolumeMapper::MINIMUM_INTENSITY_BLEND)
  {
    me->GetMIPHelper()->GenerateImage(threadID, threadCount, vol, me);
  }
  else
  {
    if (me->GetShadingRequired() == 0)
    {
      if (me->GetGradientOpacityRequired() == 0)
      {
        me->GetCompositeHelper()->GenerateImage(threadID, threadCount, vol, me);
      }
      else
      {
        me->GetCompositeGOHelper()->GenerateImage(threadID, threadCount, vol, me);
      }
    }
    else
    {
      if (me->GetGradientOpacityRequired() == 0)
      {
        me->GetCompositeShadeHelper()->GenerateImage(threadID, threadCount, vol, me);
      }
      else
      {
        me->GetCompositeGOShadeHelper()->GenerateImage(threadID, threadCount, vol, me);
      }
    }
  }

  return VTK_THREAD_RETURN_VALUE;
}

// Create an image into the vtkImageData argument. Used generally for
// creating thumbnail images
void vtkFixedPointVolumeRayCastMapper::CreateCanonicalView(
  vtkVolume* vol, vtkImageData* image, int blend_mode, double direction[3], double viewUp[3])
{
  // Make sure we have as long as we'd like so that the
  // image sample distance will be 1.0
  vol->SetAllocatedRenderTime(VTK_DOUBLE_MAX, nullptr);

  // Create a renderer / camera with the right parameters
  // These will never be mapped to the screen - just used
  // to hold parameters such as size, view direction, aspect,
  // etc.
  vtkRenderWindow* renWin = vtkRenderWindow::New();
  vtkRenderer* ren = vtkRenderer::New();
  vtkCamera* cam = ren->GetActiveCamera();

  renWin->AddRenderer(ren);
  int dim[3];
  image->GetDimensions(dim);

  // The size of the window is the size of the image
  renWin->SetSize(dim[0], dim[1]);

  double* center = vol->GetCenter();

  double bnds[6];
  vol->GetBounds(bnds);
#if 0
  double d = sqrt((bnds[1]-bnds[0])*(bnds[1]-bnds[0]) +
                  (bnds[3]-bnds[2])*(bnds[3]-bnds[2]) +
                  (bnds[5]-bnds[4])*(bnds[5]-bnds[4]));
#endif

  // For now use x distance - need to change this
  double d = bnds[1] - bnds[0];

  // Set up the camera in parallel
  cam->SetFocalPoint(center);
  cam->ParallelProjectionOn();
  cam->SetPosition(
    center[0] - d * direction[0], center[1] - d * direction[1], center[2] - d * direction[2]);
  cam->SetViewUp(viewUp);

  cam->SetParallelScale(d / 2);

  // Add a light
  vtkLight* light = vtkLight::New();
  light->SetPosition(
    center[0] - d * direction[0], center[1] - d * direction[1], center[2] - d * direction[2]);
  light->SetFocalPoint(center);
  ren->AddLight(light);

  int savedBlendMode = this->BlendMode;
  this->BlendMode = blend_mode;
  int savedCropping = this->Cropping;
  this->Cropping = 0;

  // Do all the initialization. This is not a multipass image
  // so just pass in dummy origin, spacing, and extent here
  double dummyOrigin[3] = { 0.0, 0.0, 0.0 };
  double dummySpacing[3] = { 0.0, 0.0, 0.0 };
  int dummyExtent[6] = { 0, 0, 0, 0, 0, 0 };
  this->PerImageInitialization(ren, vol, 0, dummyOrigin, dummySpacing, dummyExtent);
  this->PerVolumeInitialization(ren, vol);
  this->PerSubVolumeInitialization(ren, vol, 0);

  // Render the image
  this->RenderSubVolume();

  // Now copy the image into the vtkImageData
  unsigned char* outPtr = (unsigned char*)image->GetScalarPointer();
  unsigned short* inPtr = this->RayCastImage->GetImage();

  int viewportSize[2];
  int inUseSize[2];
  int memorySize[2];
  int origin[2];

  this->RayCastImage->GetImageViewportSize(viewportSize);
  this->RayCastImage->GetImageInUseSize(inUseSize);
  this->RayCastImage->GetImageMemorySize(memorySize);
  this->RayCastImage->GetImageOrigin(origin);

  int i, j;
  for (j = 0; j < dim[1]; j++)
  {
    for (i = 0; i < dim[0]; i++)
    {
      if (j < origin[1] || (j - origin[1]) >= inUseSize[1] || i < origin[0] ||
        (i - origin[0]) >= inUseSize[0])
      {
        *(outPtr++) = 0;
        *(outPtr++) = 0;
        *(outPtr++) = 0;
      }
      else
      {
        unsigned short* tmp = inPtr + (j - origin[1]) * memorySize[0] * 4 + (i - origin[0]) * 4;
        *(outPtr++) = (*(tmp++)) >> 7;
        *(outPtr++) = (*(tmp++)) >> 7;
        *(outPtr++) = (*(tmp++)) >> 7;
      }
    }
  }

  // Restore

  this->SampleDistance = this->OldSampleDistance;
  this->BlendMode = savedBlendMode;
  this->Cropping = savedCropping;

  // Clean up
  renWin->RemoveRenderer(ren);
  ren->RemoveLight(light);

  renWin->Delete();
  ren->Delete();
  light->Delete();
}

void vtkFixedPointVolumeRayCastMapper::ComputeRayInfo(
  int x, int y, unsigned int pos[3], unsigned int dir[3], unsigned int* numSteps)
{
  float viewRay[3];
  float rayDirection[3];
  double rayStart[4], rayEnd[4];

  int imageViewportSize[2];
  int imageOrigin[2];
  this->RayCastImage->GetImageViewportSize(imageViewportSize);
  this->RayCastImage->GetImageOrigin(imageOrigin);

  double offsetX = 1.0 / static_cast<double>(imageViewportSize[0]);
  double offsetY = 1.0 / static_cast<double>(imageViewportSize[1]);

  // compute the view point y value for this row. Do this by
  // taking our pixel position, adding the image origin then dividing
  // by the full image size to get a number from 0 to 1-1/fullSize. Then,
  // multiply by two and subtract one to get a number from
  // -1 to 1 - 2/fullSize. Then add offsetX (which is 1/fullSize) to
  // center it.
  viewRay[1] =
    ((static_cast<float>(y) + static_cast<float>(imageOrigin[1])) / imageViewportSize[1]) * 2.0 -
    1.0 + offsetY;

  // compute the view point x value for this pixel. Do this by
  // taking our pixel position, adding the image origin then dividing
  // by the full image size to get a number from 0 to 1-1/fullSize. Then,
  // multiply by two and subtract one to get a number from
  // -1 to 1 - 2/fullSize. Then add offsetX (which is 1/fullSize) to
  // center it.
  viewRay[0] =
    ((static_cast<float>(x) + static_cast<float>(imageOrigin[0])) / imageViewportSize[0]) * 2.0 -
    1.0 + offsetX;

  // Now transform this point with a z value of 0 for the ray start, and
  // a z value of 1 for the ray end. This corresponds to the near and far
  // plane locations. If IntermixIntersectingGeometry is on, then use
  // the zbuffer value instead of 1.0
  viewRay[2] = 0.0;
  vtkVRCMultiplyPointMacro(viewRay, rayStart, this->ViewToVoxelsArray);

  viewRay[2] = this->RayCastImage->GetZBufferValue(x, y);

  vtkVRCMultiplyPointMacro(viewRay, rayEnd, this->ViewToVoxelsArray);

  rayDirection[0] = rayEnd[0] - rayStart[0];
  rayDirection[1] = rayEnd[1] - rayStart[1];
  rayDirection[2] = rayEnd[2] - rayStart[2];

  float originalRayStart[3];
  originalRayStart[0] = rayStart[0];
  originalRayStart[1] = rayStart[1];
  originalRayStart[2] = rayStart[2];

  // Initialize with 0, fill in with actual number of steps
  // if necessary
  *numSteps = 0;

  if (this->ClipRayAgainstVolume(rayStart, rayEnd, rayDirection, this->CroppingBounds) &&
    (this->NumTransformedClippingPlanes == 0 ||
      this->ClipRayAgainstClippingPlanes(
        rayStart, rayEnd, this->NumTransformedClippingPlanes, this->TransformedClippingPlanes)))
  {
    double worldRayDirection[3];
    worldRayDirection[0] = rayDirection[0] * this->SavedSpacing[0];
    worldRayDirection[1] = rayDirection[1] * this->SavedSpacing[1];
    worldRayDirection[2] = rayDirection[2] * this->SavedSpacing[2];
    double worldLength = vtkMath::Normalize(worldRayDirection) / this->SampleDistance;

    rayDirection[0] /= worldLength;
    rayDirection[1] /= worldLength;
    rayDirection[2] /= worldLength;

    float diff[3];
    diff[0] = (rayStart[0] - originalRayStart[0]) * ((rayDirection[0] < 0) ? (-1) : (1));
    diff[1] = (rayStart[1] - originalRayStart[1]) * ((rayDirection[1] < 0) ? (-1) : (1));
    diff[2] = (rayStart[2] - originalRayStart[2]) * ((rayDirection[2] < 0) ? (-1) : (1));

    int steps = -1;

    if (diff[0] >= diff[1] && diff[0] >= diff[2] && rayDirection[0])
    {
      steps = 1 +
        static_cast<int>(
          diff[0] / ((rayDirection[0] < 0) ? (-rayDirection[0]) : (rayDirection[0])));
    }

    if (diff[1] >= diff[0] && diff[1] >= diff[2] && rayDirection[1])
    {
      steps = 1 +
        static_cast<int>(
          diff[1] / ((rayDirection[1] < 0) ? (-rayDirection[1]) : (rayDirection[1])));
    }

    if (diff[2] >= diff[0] && diff[2] >= diff[1] && rayDirection[2])
    {
      steps = 1 +
        static_cast<int>(
          diff[2] / ((rayDirection[2] < 0) ? (-rayDirection[2]) : (rayDirection[2])));
    }

    if (steps > 0)
    {
      rayStart[0] = originalRayStart[0] + steps * rayDirection[0];
      rayStart[1] = originalRayStart[1] + steps * rayDirection[1];
      rayStart[2] = originalRayStart[2] + steps * rayDirection[2];
    }

    if (rayStart[0] > 0.0 && rayStart[1] > 0.0 && rayStart[2] > 0.0)
    {
      pos[0] = this->ToFixedPointPosition(rayStart[0]);
      pos[1] = this->ToFixedPointPosition(rayStart[1]);
      pos[2] = this->ToFixedPointPosition(rayStart[2]);
      dir[0] = this->ToFixedPointDirection(rayDirection[0]);
      dir[1] = this->ToFixedPointDirection(rayDirection[1]);
      dir[2] = this->ToFixedPointDirection(rayDirection[2]);

      int stepLoop;
      int stepsValid = 0;
      unsigned int currSteps;
      for (stepLoop = 0; stepLoop < 3; stepLoop++)
      {
        if (!(dir[stepLoop] & 0x7fffffff))
        {
          continue;
        }

        unsigned int endVal = this->ToFixedPointPosition(rayEnd[stepLoop]);

        if (dir[stepLoop] & 0x80000000)
        {
          if (endVal > pos[stepLoop])
          {
            currSteps = static_cast<unsigned int>(
              1 + (endVal - pos[stepLoop]) / (dir[stepLoop] & 0x7fffffff));
          }
          else
          {
            currSteps = 0;
          }
        }
        else
        {
          if (pos[stepLoop] > endVal)
          {
            currSteps = 1 + (pos[stepLoop] - endVal) / dir[stepLoop];
          }
          else
          {
            currSteps = 0;
          }
        }

        if (!stepsValid || currSteps < *numSteps)
        {
          *numSteps = currSteps;
          stepsValid = 1;
        }
      }
    }
  }
}

void vtkFixedPointVolumeRayCastMapper::InitializeRayInfo(vtkVolume* vol)
{
  if (!vol)
  {
    return;
  }

  // Copy the viewToVoxels matrix to 16 floats
  int i, j;
  for (j = 0; j < 4; j++)
  {
    for (i = 0; i < 4; i++)
    {
      this->ViewToVoxelsArray[j * 4 + i] =
        static_cast<float>(this->ViewToVoxelsMatrix->GetElement(j, i));
    }
  }

  // Copy the worldToVoxels matrix to 16 floats
  for (j = 0; j < 4; j++)
  {
    for (i = 0; i < 4; i++)
    {
      this->WorldToVoxelsArray[j * 4 + i] =
        static_cast<float>(this->WorldToVoxelsMatrix->GetElement(j, i));
    }
  }

  // Copy the voxelsToWorld matrix to 16 floats
  for (j = 0; j < 4; j++)
  {
    for (i = 0; i < 4; i++)
    {
      this->VoxelsToWorldArray[j * 4 + i] =
        static_cast<float>(this->VoxelsToWorldMatrix->GetElement(j, i));
    }
  }

  int dim[3];
  vtkImageData::SafeDownCast(this->GetInput())->GetDimensions(dim);
  this->CroppingBounds[0] = this->CroppingBounds[2] = this->CroppingBounds[4] = 0.0;
  this->CroppingBounds[1] = dim[0] - 1;
  this->CroppingBounds[3] = dim[1] - 1;
  this->CroppingBounds[5] = dim[2] - 1;

  // Do some initialization of the clipping planes
  this->NumTransformedClippingPlanes =
    (this->ClippingPlanes) ? (this->ClippingPlanes->GetNumberOfItems()) : (0);

  // Clear out old clipping planes
  delete[] this->TransformedClippingPlanes;
  this->TransformedClippingPlanes = nullptr;

  // Do we have any clipping planes
  if (this->NumTransformedClippingPlanes > 0)
  {
    // Allocate some space to store the plane equations
    this->TransformedClippingPlanes = new float[4 * this->NumTransformedClippingPlanes];

    // loop through all the clipping planes
    for (i = 0; i < this->NumTransformedClippingPlanes; i++)
    {
      // Convert plane into voxel coordinate system
      double worldNormal[3], worldOrigin[3];
      double inputOrigin[4];
      vtkPlane* onePlane = (vtkPlane*)this->ClippingPlanes->GetItemAsObject(i);
      onePlane->GetNormal(worldNormal);
      onePlane->GetOrigin(worldOrigin);
      float* planePtr = this->TransformedClippingPlanes + 4 * i;
      vtkVRCMultiplyNormalMacro(worldNormal, planePtr, this->VoxelsToWorldArray);
      vtkVRCMultiplyPointMacro(worldOrigin, inputOrigin, this->WorldToVoxelsArray);

      float t =
        sqrt(planePtr[0] * planePtr[0] + planePtr[1] * planePtr[1] + planePtr[2] * planePtr[2]);
      if (t)
      {
        planePtr[0] /= t;
        planePtr[1] /= t;
        planePtr[2] /= t;
      }

      planePtr[3] = -(
        planePtr[0] * inputOrigin[0] + planePtr[1] * inputOrigin[1] + planePtr[2] * inputOrigin[2]);
    }
  }

  // If we have a simple crop box then we can tighten the bounds
  if (this->Cropping && this->CroppingRegionFlags == 0x2000)
  {
    this->CroppingBounds[0] = this->VoxelCroppingRegionPlanes[0];
    this->CroppingBounds[1] = this->VoxelCroppingRegionPlanes[1];
    this->CroppingBounds[2] = this->VoxelCroppingRegionPlanes[2];
    this->CroppingBounds[3] = this->VoxelCroppingRegionPlanes[3];
    this->CroppingBounds[4] = this->VoxelCroppingRegionPlanes[4];
    this->CroppingBounds[5] = this->VoxelCroppingRegionPlanes[5];
  }

  this->CroppingBounds[0] = (this->CroppingBounds[0] < 0) ? (0) : (this->CroppingBounds[0]);
  this->CroppingBounds[0] =
    (this->CroppingBounds[0] > dim[0] - 1) ? (dim[0] - 1) : (this->CroppingBounds[0]);
  this->CroppingBounds[1] = (this->CroppingBounds[1] < 0) ? (0) : (this->CroppingBounds[1]);
  this->CroppingBounds[1] =
    (this->CroppingBounds[1] > dim[0] - 1) ? (dim[0] - 1) : (this->CroppingBounds[1]);
  this->CroppingBounds[2] = (this->CroppingBounds[2] < 0) ? (0) : (this->CroppingBounds[2]);
  this->CroppingBounds[2] =
    (this->CroppingBounds[2] > dim[1] - 1) ? (dim[1] - 1) : (this->CroppingBounds[2]);
  this->CroppingBounds[3] = (this->CroppingBounds[3] < 0) ? (0) : (this->CroppingBounds[3]);
  this->CroppingBounds[3] =
    (this->CroppingBounds[3] > dim[1] - 1) ? (dim[1] - 1) : (this->CroppingBounds[3]);
  this->CroppingBounds[4] = (this->CroppingBounds[4] < 0) ? (0) : (this->CroppingBounds[4]);
  this->CroppingBounds[4] =
    (this->CroppingBounds[4] > dim[2] - 1) ? (dim[2] - 1) : (this->CroppingBounds[4]);
  this->CroppingBounds[5] = (this->CroppingBounds[5] < 0) ? (0) : (this->CroppingBounds[5]);
  this->CroppingBounds[5] =
    (this->CroppingBounds[5] > dim[2] - 1) ? (dim[2] - 1) : (this->CroppingBounds[5]);

  // Save spacing because for some reason this call is really really slow!
  vtkImageData::SafeDownCast(this->GetInput())->GetSpacing(this->SavedSpacing);
}

// Return 0 if our volume is outside the view frustum, 1 if it
// is in the view frustum.
int vtkFixedPointVolumeRayCastMapper::ComputeRowBounds(
  vtkRenderer* ren, int imageFlag, int rowBoundsFlag, int inputExtent[6])
{
  float voxelPoint[3];
  double viewPoint[8][4];
  int i, j, k;
  unsigned short* ucptr;
  float minX, minY, maxX, maxY, minZ, maxZ;

  minX = 1.0;
  minY = 1.0;
  maxX = -1.0;
  maxY = -1.0;
  minZ = 1.0;
  maxZ = 0.0;

  float bounds[6];
  int dim[3];
  dim[0] = inputExtent[1] - inputExtent[0] + 1;
  dim[1] = inputExtent[3] - inputExtent[2] + 1;
  dim[2] = inputExtent[5] - inputExtent[4] + 1;

  bounds[0] = bounds[2] = bounds[4] = 0.0;
  bounds[1] = static_cast<float>(dim[0] - 1);
  bounds[3] = static_cast<float>(dim[1] - 1);
  bounds[5] = static_cast<float>(dim[2] - 1);

  int insideFlag = 0;
  double camPos[4];
  ren->GetActiveCamera()->GetPosition(camPos);
  camPos[3] = 1.0;
  this->WorldToVoxelsMatrix->MultiplyPoint(camPos, camPos);
  if (camPos[3])
  {
    camPos[0] /= camPos[3];
    camPos[1] /= camPos[3];
    camPos[2] /= camPos[3];
  }

  // If we have a simple crop box then we can tighten the bounds
  if (this->Cropping && this->CroppingRegionFlags == 0x2000)
  {
    bounds[0] = this->VoxelCroppingRegionPlanes[0];
    bounds[1] = this->VoxelCroppingRegionPlanes[1];
    bounds[2] = this->VoxelCroppingRegionPlanes[2];
    bounds[3] = this->VoxelCroppingRegionPlanes[3];
    bounds[4] = this->VoxelCroppingRegionPlanes[4];
    bounds[5] = this->VoxelCroppingRegionPlanes[5];
  }

  if (camPos[0] >= bounds[0] && camPos[0] <= bounds[1] && camPos[1] >= bounds[2] &&
    camPos[1] <= bounds[3] && camPos[2] >= bounds[4] && camPos[2] <= bounds[5])
  {
    insideFlag = 1;
  }

  // Copy the voxelsToView matrix to 16 floats
  float voxelsToViewMatrix[16];
  for (j = 0; j < 4; j++)
  {
    for (i = 0; i < 4; i++)
    {
      voxelsToViewMatrix[j * 4 + i] =
        static_cast<float>(this->VoxelsToViewMatrix->GetElement(j, i));
    }
  }

  // Convert the voxel bounds to view coordinates to find out the
  // size and location of the image we need to generate.
  int idx = 0;
  if (insideFlag)
  {
    minX = -1.0;
    maxX = 1.0;
    minY = -1.0;
    maxY = 1.0;
    minZ = 0.001;
    maxZ = 0.001;
  }
  else
  {
    for (k = 0; k < 2; k++)
    {
      voxelPoint[2] = bounds[4 + k];
      for (j = 0; j < 2; j++)
      {
        voxelPoint[1] = bounds[2 + j];
        for (i = 0; i < 2; i++)
        {
          voxelPoint[0] = bounds[i];
          vtkVRCMultiplyPointMacro(voxelPoint, viewPoint[idx], voxelsToViewMatrix);

          minX = (viewPoint[idx][0] < minX) ? (viewPoint[idx][0]) : (minX);
          minY = (viewPoint[idx][1] < minY) ? (viewPoint[idx][1]) : (minY);
          maxX = (viewPoint[idx][0] > maxX) ? (viewPoint[idx][0]) : (maxX);
          maxY = (viewPoint[idx][1] > maxY) ? (viewPoint[idx][1]) : (maxY);
          minZ = (viewPoint[idx][2] < minZ) ? (viewPoint[idx][2]) : (minZ);
          maxZ = (viewPoint[idx][2] > maxZ) ? (viewPoint[idx][2]) : (maxZ);
          idx++;
        }
      }
    }
  }

  if (minZ < 0.001 || maxZ > 0.9999)
  {
    minX = -1.0;
    maxX = 1.0;
    minY = -1.0;
    maxY = 1.0;
    insideFlag = 1;
  }

  this->MinimumViewDistance = (minZ < 0.001) ? (0.001) : ((minZ > 0.999) ? (0.999) : (minZ));

  int imageViewportSize[2];
  int imageOrigin[2];
  int imageMemorySize[2];
  int imageInUseSize[2];
  this->RayCastImage->GetImageViewportSize(imageViewportSize);
  this->RayCastImage->GetImageOrigin(imageOrigin);
  this->RayCastImage->GetImageMemorySize(imageMemorySize);

  // We have min/max values from -1.0 to 1.0 now - we want to convert
  // these to pixel locations. Give a couple of pixels of breathing room
  // on each side if possible
  minX = (minX + 1.0) * 0.5 * imageViewportSize[0] - 2;
  minY = (minY + 1.0) * 0.5 * imageViewportSize[1] - 2;
  maxX = (maxX + 1.0) * 0.5 * imageViewportSize[0] + 2;
  maxY = (maxY + 1.0) * 0.5 * imageViewportSize[1] + 2;

  // If we are outside the view frustum return 0 - there is no need
  // to render anything
  if ((minX < 0 && maxX < 0) || (minY < 0 && maxY < 0) ||
    (minX > imageViewportSize[0] - 1 && maxX > imageViewportSize[0] - 1) ||
    (minY > imageViewportSize[1] - 1 && maxY > imageViewportSize[1] - 1))
  {
    return 0;
  }

  int oldImageMemorySize[2];
  oldImageMemorySize[0] = imageMemorySize[0];
  oldImageMemorySize[1] = imageMemorySize[1];

  // Check the bounds - the volume might project outside of the
  // viewing box / frustum so clip it if necessary
  minX = (minX < 0) ? (0) : (minX);
  minY = (minY < 0) ? (0) : (minY);
  maxX = (maxX > imageViewportSize[0] - 1) ? (imageViewportSize[0] - 1) : (maxX);
  maxY = (maxY > imageViewportSize[1] - 1) ? (imageViewportSize[1] - 1) : (maxY);

  // Create the new image, and set its size and position
  imageInUseSize[0] = static_cast<int>(maxX - minX + 1.0);
  imageInUseSize[1] = static_cast<int>(maxY - minY + 1.0);

  // What is a power of 2 size big enough to fit this image?
  imageMemorySize[0] = 32;
  imageMemorySize[1] = 32;
  while (imageMemorySize[0] < imageInUseSize[0])
  {
    imageMemorySize[0] *= 2;
  }
  while (imageMemorySize[1] < imageInUseSize[1])
  {
    imageMemorySize[1] *= 2;
  }

  imageOrigin[0] = static_cast<int>(minX);
  imageOrigin[1] = static_cast<int>(minY);

  // If the old image size is much too big (more than twice in
  // either direction) then set the old width to 0 which will
  // cause the image to be recreated
  if (oldImageMemorySize[0] > 4 * imageMemorySize[0] ||
    oldImageMemorySize[1] > 4 * imageMemorySize[1])
  {
    oldImageMemorySize[0] = 0;
  }

  // If the old image is big enough (but not too big - we handled
  // that above) then we'll bump up our required size to the
  // previous one. This will keep us from thrashing.
  if (oldImageMemorySize[0] >= imageMemorySize[0] && oldImageMemorySize[1] >= imageMemorySize[1])
  {
    imageMemorySize[0] = oldImageMemorySize[0];
    imageMemorySize[1] = oldImageMemorySize[1];
  }

  if (imageFlag)
  {
    this->RayCastImage->SetImageOrigin(imageOrigin);
    this->RayCastImage->SetImageMemorySize(imageMemorySize);
    this->RayCastImage->SetImageInUseSize(imageInUseSize);

    // Do we already have a texture big enough? If not, create a new one and
    // clear it.
    if (imageMemorySize[0] > oldImageMemorySize[0] || imageMemorySize[1] > oldImageMemorySize[1])
    {
      this->RayCastImage->AllocateImage();
      delete[] this->RowBounds;
      delete[] this->OldRowBounds;

      this->RayCastImage->ClearImage();

      if (rowBoundsFlag)
      {
        // Create the row bounds array. This will store the start / stop pixel
        // for each row. This helps eliminate work in areas outside the bounding
        // hexahedron since a bounding box is not very tight. We keep the old ones
        // too to help with only clearing where required
        this->RowBounds = new int[2 * imageMemorySize[1]];
        this->OldRowBounds = new int[2 * imageMemorySize[1]];

        for (i = 0; i < imageMemorySize[1]; i++)
        {
          this->RowBounds[i * 2] = imageMemorySize[0];
          this->RowBounds[i * 2 + 1] = -1;
          this->OldRowBounds[i * 2] = imageMemorySize[0];
          this->OldRowBounds[i * 2 + 1] = -1;
        }
      }
    }
  }

  if (!rowBoundsFlag)
  {
    return 1;
  }

  // Swap the row bounds
  int* tmpptr;
  tmpptr = this->RowBounds;
  this->RowBounds = this->OldRowBounds;
  this->OldRowBounds = tmpptr;

  // If we are inside the volume our row bounds indicate every ray must be
  // cast - we don't need to intersect with the 12 lines
  if (insideFlag)
  {
    for (j = 0; j < imageInUseSize[1]; j++)
    {
      this->RowBounds[j * 2] = 0;
      this->RowBounds[j * 2 + 1] = imageInUseSize[0] - 1;
    }
  }
  else
  {
    // create an array of lines where the y value of the first vertex is less
    // than or equal to the y value of the second vertex. There are 12 lines,
    // each containing x1, y1, x2, y2 values.
    float lines[12][4];
    float x1, y1, x2, y2;
    int xlow, xhigh;
    int lineIndex[12][2] = { { 0, 1 }, { 2, 3 }, { 4, 5 }, { 6, 7 }, { 0, 2 }, { 1, 3 }, { 4, 6 },
      { 5, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 } };

    for (i = 0; i < 12; i++)
    {
      x1 = (viewPoint[lineIndex[i][0]][0] + 1.0) * 0.5 * imageViewportSize[0] - imageOrigin[0];

      y1 = (viewPoint[lineIndex[i][0]][1] + 1.0) * 0.5 * imageViewportSize[1] - imageOrigin[1];

      x2 = (viewPoint[lineIndex[i][1]][0] + 1.0) * 0.5 * imageViewportSize[0] - imageOrigin[0];

      y2 = (viewPoint[lineIndex[i][1]][1] + 1.0) * 0.5 * imageViewportSize[1] - imageOrigin[1];

      if (y1 < y2)
      {
        lines[i][0] = x1;
        lines[i][1] = y1;
        lines[i][2] = x2;
        lines[i][3] = y2;
      }
      else
      {
        lines[i][0] = x2;
        lines[i][1] = y2;
        lines[i][2] = x1;
        lines[i][3] = y1;
      }
    }

    // Now for each row in the image, find out the start / stop pixel
    // If min > max, then no intersection occurred
    for (j = 0; j < imageInUseSize[1]; j++)
    {
      this->RowBounds[j * 2] = imageMemorySize[0];
      this->RowBounds[j * 2 + 1] = -1;
      for (i = 0; i < 12; i++)
      {
        if (j >= lines[i][1] && j <= lines[i][3] && (lines[i][1] != lines[i][3]))
        {
          x1 = lines[i][0] +
            (static_cast<float>(j) - lines[i][1]) / (lines[i][3] - lines[i][1]) *
              (lines[i][2] - lines[i][0]);

          xlow = static_cast<int>(x1 + 1.5);
          xhigh = static_cast<int>(x1 - 1.0);

          xlow = (xlow < 0) ? (0) : (xlow);
          xlow = (xlow > imageInUseSize[0] - 1) ? (imageInUseSize[0] - 1) : (xlow);

          xhigh = (xhigh < 0) ? (0) : (xhigh);
          xhigh = (xhigh > imageInUseSize[0] - 1) ? (imageInUseSize[0] - 1) : (xhigh);

          if (xlow < this->RowBounds[j * 2])
          {
            this->RowBounds[j * 2] = xlow;
          }
          if (xhigh > this->RowBounds[j * 2 + 1])
          {
            this->RowBounds[j * 2 + 1] = xhigh;
          }
        }
      }
      // If they are the same this is either a point on the cube or
      // all lines were out of bounds (all on one side or the other)
      // It is safe to ignore the point (since the ray isn't likely
      // to travel through it enough to actually take a sample) and it
      // must be ignored in the case where all lines are out of range
      if (this->RowBounds[j * 2] == this->RowBounds[j * 2 + 1])
      {
        this->RowBounds[j * 2] = imageMemorySize[0];
        this->RowBounds[j * 2 + 1] = -1;
      }
    }
  }

  for (j = imageInUseSize[1]; j < imageMemorySize[1]; j++)
  {
    this->RowBounds[j * 2] = imageMemorySize[0];
    this->RowBounds[j * 2 + 1] = -1;
  }

  unsigned short* image = this->RayCastImage->GetImage();

  for (j = 0; j < imageMemorySize[1]; j++)
  {
    if (j % 64 == 1 && this->RenderWindow && this->RenderWindow->CheckAbortStatus())
    {
      return 0;
    }

    // New bounds are not overlapping with old bounds - clear between
    // old bounds only
    if (this->RowBounds[j * 2 + 1] < this->OldRowBounds[j * 2] ||
      this->RowBounds[j * 2] > this->OldRowBounds[j * 2 + 1])
    {
      ucptr = image + 4 * (j * imageMemorySize[0] + this->OldRowBounds[j * 2]);
      for (i = 0; i <= (this->OldRowBounds[j * 2 + 1] - this->OldRowBounds[j * 2]); i++)
      {
        *(ucptr++) = 0;
        *(ucptr++) = 0;
        *(ucptr++) = 0;
        *(ucptr++) = 0;
      }
    }
    // New bounds do overlap with old bounds
    else
    {
      // Clear from old min to new min
      ucptr = image + 4 * (j * imageMemorySize[0] + this->OldRowBounds[j * 2]);
      for (i = 0; i < (this->RowBounds[j * 2] - this->OldRowBounds[j * 2]); i++)
      {
        *(ucptr++) = 0;
        *(ucptr++) = 0;
        *(ucptr++) = 0;
        *(ucptr++) = 0;
      }

      // Clear from new max to old max
      ucptr = image + 4 * (j * imageMemorySize[0] + this->RowBounds[j * 2 + 1] + 1);
      for (i = 0; i < (this->OldRowBounds[j * 2 + 1] - this->RowBounds[j * 2 + 1]); i++)
      {
        *(ucptr++) = 0;
        *(ucptr++) = 0;
        *(ucptr++) = 0;
        *(ucptr++) = 0;
      }
    }
  }

  return 1;
}

void vtkFixedPointVolumeRayCastMapper::ComputeMatrices(double inputOrigin[3],
  double inputSpacing[3], int inputExtent[6], vtkRenderer* ren, vtkVolume* vol)
{
  // Get the camera from the renderer
  vtkCamera* cam = ren->GetActiveCamera();

  // Get the aspect ratio from the renderer. This is needed for the
  // computation of the perspective matrix
  ren->ComputeAspect();
  double* aspect = ren->GetAspect();

  // Keep track of the projection matrix - we'll need it in a couple of places
  // Get the projection matrix. The method is called perspective, but
  // the matrix is valid for perspective and parallel viewing transforms.
  // Don't replace this with the GetCompositePerspectiveTransformMatrix
  // because that turns off stereo rendering!!!
  this->PerspectiveTransform->Identity();
  this->PerspectiveTransform->Concatenate(
    cam->GetProjectionTransformMatrix(aspect[0] / aspect[1], 0.0, 1.0));
  this->PerspectiveTransform->Concatenate(cam->GetViewTransformMatrix());
  this->PerspectiveMatrix->DeepCopy(this->PerspectiveTransform->GetMatrix());

  // Compute the origin of the extent the volume origin is at voxel (0,0,0)
  // but we want to consider (0,0,0) in voxels to be at
  // (inputExtent[0], inputExtent[2], inputExtent[4]).
  double extentOrigin[3];
  extentOrigin[0] = inputOrigin[0] + inputExtent[0] * inputSpacing[0];
  extentOrigin[1] = inputOrigin[1] + inputExtent[2] * inputSpacing[1];
  extentOrigin[2] = inputOrigin[2] + inputExtent[4] * inputSpacing[2];

  // Get the volume matrix. This is a volume to world matrix right now.
  // We'll need to invert it, translate by the origin and scale by the
  // spacing to change it to a world to voxels matrix.
  vol->GetModelToWorldMatrix(this->VolumeMatrix);

  this->VoxelsToViewTransform->SetMatrix(this->VolumeMatrix);

  // Create a transform that will account for the scaling and translation of
  // the scalar data. The is the volume to voxels matrix.
  this->VoxelsTransform->Identity();
  this->VoxelsTransform->Translate(extentOrigin[0], extentOrigin[1], extentOrigin[2]);

  this->VoxelsTransform->Scale(inputSpacing[0], inputSpacing[1], inputSpacing[2]);

  // Now concatenate the volume's matrix with this scalar data matrix
  this->VoxelsToViewTransform->PreMultiply();
  this->VoxelsToViewTransform->Concatenate(this->VoxelsTransform->GetMatrix());

  // Now we actually have the world to voxels matrix - copy it out
  this->WorldToVoxelsMatrix->DeepCopy(this->VoxelsToViewTransform->GetMatrix());
  this->WorldToVoxelsMatrix->Invert();

  // We also want to invert this to get voxels to world
  this->VoxelsToWorldMatrix->DeepCopy(this->VoxelsToViewTransform->GetMatrix());

  // Compute the voxels to view transform by concatenating the
  // voxels to world matrix with the projection matrix (world to view)
  this->VoxelsToViewTransform->PostMultiply();
  this->VoxelsToViewTransform->Concatenate(this->PerspectiveMatrix);

  this->VoxelsToViewMatrix->DeepCopy(this->VoxelsToViewTransform->GetMatrix());

  this->ViewToVoxelsMatrix->DeepCopy(this->VoxelsToViewMatrix);
  this->ViewToVoxelsMatrix->Invert();
}

int vtkFixedPointVolumeRayCastMapper::ClipRayAgainstClippingPlanes(
  double rayStart[3], double rayEnd[3], int numClippingPlanes, float* clippingPlanes)
{

  float* planePtr;
  int i;
  float t, point[3], dp;
  float rayDir[3];

  rayDir[0] = rayEnd[0] - rayStart[0];
  rayDir[1] = rayEnd[1] - rayStart[1];
  rayDir[2] = rayEnd[2] - rayStart[2];

  // loop through all the clipping planes
  for (i = 0; i < numClippingPlanes; i++)
  {
    planePtr = clippingPlanes + 4 * i;

    dp = planePtr[0] * rayDir[0] + planePtr[1] * rayDir[1] + planePtr[2] * rayDir[2];

    if (dp != 0.0)
    {
      t = -(planePtr[0] * rayStart[0] + planePtr[1] * rayStart[1] + planePtr[2] * rayStart[2] +
            planePtr[3]) /
        dp;

      if (t > 0.0 && t < 1.0)
      {
        point[0] = rayStart[0] + t * rayDir[0];
        point[1] = rayStart[1] + t * rayDir[1];
        point[2] = rayStart[2] + t * rayDir[2];

        if (dp > 0.0)
        {
          rayStart[0] = point[0];
          rayStart[1] = point[1];
          rayStart[2] = point[2];
        }
        else
        {
          rayEnd[0] = point[0];
          rayEnd[1] = point[1];
          rayEnd[2] = point[2];
        }

        rayDir[0] = rayEnd[0] - rayStart[0];
        rayDir[1] = rayEnd[1] - rayStart[1];
        rayDir[2] = rayEnd[2] - rayStart[2];
      }
      // If the clipping plane is outside the ray segment, then
      // figure out if that means the ray segment goes to zero (if so
      // return 0) or doesn't affect it (if so do nothing)
      else
      {
        if (dp >= 0.0 && t >= 1.0)
        {
          return 0;
        }
        if (dp <= 0.0 && t <= 0.0)
        {
          return 0;
        }
      }
    }
    else
    {
      // rayDir is perpendicular to planePtr; i.e., the ray does not
      // intersect the plane. Determine which side of the plane the ray
      // is on.
      float side = planePtr[0] * rayStart[0] + planePtr[1] * rayStart[1] +
        planePtr[2] * rayStart[2] + planePtr[3];
      if (side < 0)
      {
        // clip
        return 0;
      }
    }
  }

  return 1;
}

int vtkFixedPointVolumeRayCastMapper::ClipRayAgainstVolume(
  double rayStart[3], double rayEnd[3], float rayDirection[3], double bounds[6])
{
  int loop;
  float diff;
  float t;

  if (rayStart[0] >= bounds[1] || rayStart[1] >= bounds[3] || rayStart[2] >= bounds[5] ||
    rayStart[0] < bounds[0] || rayStart[1] < bounds[2] || rayStart[2] < bounds[4])
  {
    for (loop = 0; loop < 3; loop++)
    {
      diff = 0;

      if (rayStart[loop] < (bounds[2 * loop] + 0.01))
      {
        diff = (bounds[2 * loop] + 0.01) - rayStart[loop];
      }
      else if (rayStart[loop] > (bounds[2 * loop + 1] - 0.01))
      {
        diff = (bounds[2 * loop + 1] - 0.01) - rayStart[loop];
      }

      if (diff)
      {
        if (rayDirection[loop] != 0.0)
        {
          t = diff / rayDirection[loop];
        }
        else
        {
          t = -1.0;
        }

        if (t > 0.0)
        {
          rayStart[0] += rayDirection[0] * t;
          rayStart[1] += rayDirection[1] * t;
          rayStart[2] += rayDirection[2] * t;
        }
      }
    }
  }

  // If the voxel still isn't inside the volume, then this ray
  // doesn't really intersect the volume

  if (rayStart[0] >= bounds[1] || rayStart[1] >= bounds[3] || rayStart[2] >= bounds[5] ||
    rayStart[0] < bounds[0] || rayStart[1] < bounds[2] || rayStart[2] < bounds[4])
  {
    return 0;
  }

  // The ray does intersect the volume, and we have a starting
  // position that is inside the volume
  if (rayEnd[0] >= bounds[1] || rayEnd[1] >= bounds[3] || rayEnd[2] >= bounds[5] ||
    rayEnd[0] < bounds[0] || rayEnd[1] < bounds[2] || rayEnd[2] < bounds[4])
  {
    for (loop = 0; loop < 3; loop++)
    {
      diff = 0;

      if (rayEnd[loop] < (bounds[2 * loop] + 0.01))
      {
        diff = (bounds[2 * loop] + 0.01) - rayEnd[loop];
      }
      else if (rayEnd[loop] > (bounds[2 * loop + 1] - 0.01))
      {
        diff = (bounds[2 * loop + 1] - 0.01) - rayEnd[loop];
      }

      if (diff)
      {
        if (rayDirection[loop] != 0.0)
        {
          t = diff / rayDirection[loop];
        }
        else
        {
          t = 1.0;
        }

        if (t < 0.0)
        {
          rayEnd[0] += rayDirection[0] * t;
          rayEnd[1] += rayDirection[1] * t;
          rayEnd[2] += rayDirection[2] * t;
        }
      }
    }
  }

  // To be absolutely certain our ray remains inside the volume,
  // recompute the ray direction (since it has changed - it is not
  // normalized and therefore changes when start/end change) and move
  // the start/end points in by 1/1000th of the distance.
  float offset;
  offset = (rayEnd[0] - rayStart[0]) * 0.001;
  rayStart[0] += offset;
  rayEnd[0] -= offset;

  offset = (rayEnd[1] - rayStart[1]) * 0.001;
  rayStart[1] += offset;
  rayEnd[1] -= offset;

  offset = (rayEnd[2] - rayStart[2]) * 0.001;
  rayStart[2] += offset;
  rayEnd[2] -= offset;

  if (rayEnd[0] >= bounds[1] || rayEnd[1] >= bounds[3] || rayEnd[2] >= bounds[5] ||
    rayEnd[0] < bounds[0] || rayEnd[1] < bounds[2] || rayEnd[2] < bounds[4])
  {
    return 0;
  }

  if ((rayEnd[0] - rayStart[0]) * rayDirection[0] < 0.0 ||
    (rayEnd[1] - rayStart[1]) * rayDirection[1] < 0.0 ||
    (rayEnd[2] - rayStart[2]) * rayDirection[2] < 0.0)
  {
    return 0;
  }

  return 1;
}

void vtkFixedPointVolumeRayCastMapper::ComputeGradients(vtkVolume* vol)
{
  vtkImageData* input = vtkImageData::SafeDownCast(this->GetInput());

  void* dataPtr = this->CurrentScalars->GetVoidPointer(0);

  int scalarType = this->CurrentScalars->GetDataType();
  int components = this->CurrentScalars->GetNumberOfComponents();
  int independent = vol->GetProperty()->GetIndependentComponents();

  int dim[3];
  double spacing[3];
  input->GetDimensions(dim);
  input->GetSpacing(spacing);

  // Find the scalar range
  double scalarRange[4][2];
  int c;
  for (c = 0; c < components; c++)
  {
    this->CurrentScalars->GetRange(scalarRange[c], c);
  }

  vtkIdType sliceSize = (static_cast<vtkIdType>(dim[0]) * static_cast<vtkIdType>(dim[1]) *
    ((independent) ? (components) : (1)));
  vtkIdType numSlices = dim[2];

  int i;

  // Delete the prior gradient normal information
  if (this->GradientNormal)
  {
    // Contiguous? Delete in one chunk otherwise delete slice by slice
    if (this->ContiguousGradientNormal)
    {
      delete[] this->ContiguousGradientNormal;
      this->ContiguousGradientNormal = nullptr;
    }
    else
    {
      for (i = 0; i < this->NumberOfGradientSlices; i++)
      {
        delete[] this->GradientNormal[i];
      }
    }
    delete[] this->GradientNormal;
    this->GradientNormal = nullptr;
  }

  // Delete the prior gradient magnitude information
  if (this->GradientMagnitude)
  {
    // Contiguous? Delete in one chunk otherwise delete slice by slice
    if (this->ContiguousGradientMagnitude)
    {
      delete[] this->ContiguousGradientMagnitude;
      this->ContiguousGradientMagnitude = nullptr;
    }
    else
    {
      for (i = 0; i < this->NumberOfGradientSlices; i++)
      {
        delete[] this->GradientMagnitude[i];
      }
    }
    delete[] this->GradientMagnitude;
    this->GradientMagnitude = nullptr;
  }

  this->NumberOfGradientSlices = numSlices;
  this->GradientNormal = new unsigned short*[numSlices];
  this->GradientMagnitude = new unsigned char*[numSlices];

  // first, attempt contiguous memory. If this fails, then go
  // for non-contiguous
  // NOTE: Standard behavior is to catch std::bad_alloc, but it's
  // CMemoryException if hosted by MFC.
  try
  {
    this->ContiguousGradientNormal = new unsigned short[numSlices * sliceSize];
  }
  catch (...)
  {
    this->ContiguousGradientNormal = nullptr;
  }
  try
  {
    this->ContiguousGradientMagnitude = new unsigned char[numSlices * sliceSize];
  }
  catch (...)
  {
    this->ContiguousGradientMagnitude = nullptr;
  }

  if (this->ContiguousGradientNormal)
  {
    // We were able to allocate contiguous space - we just need to set the
    // slice pointers here
    for (i = 0; i < numSlices; i++)
    {
      this->GradientNormal[i] = this->ContiguousGradientNormal + i * sliceSize;
    }
  }
  else
  {
    // We were not able to allocate contiguous space - allocate it slice by slice
    for (i = 0; i < numSlices; i++)
    {
      this->GradientNormal[i] = new unsigned short[sliceSize];
    }
  }

  if (this->ContiguousGradientMagnitude)
  {
    // We were able to allocate contiguous space - we just need to set the
    // slice pointers here
    for (i = 0; i < numSlices; i++)
    {
      this->GradientMagnitude[i] = this->ContiguousGradientMagnitude + i * sliceSize;
    }
  }
  else
  {
    // We were not able to allocate contiguous space - allocate it slice by slice
    for (i = 0; i < numSlices; i++)
    {
      this->GradientMagnitude[i] = new unsigned char[sliceSize];
    }
  }

  vtkTimerLog* timer = vtkTimerLog::New();
  timer->StartTimer();

  if (components == 1 &&
    (scalarType == VTK_UNSIGNED_CHAR || scalarType == VTK_CHAR ||
      scalarType == VTK_UNSIGNED_SHORT || scalarType == VTK_SHORT))
  {
    this->Threader->SetSingleMethod(vtkFPVRCMSwitchOnDataType, (vtkObject*)this);
    this->Threader->SingleMethodExecute();
  }

  else
  {
    switch (scalarType)
    {
      vtkTemplateMacro(vtkFixedPointVolumeRayCastMapperComputeGradients((VTK_TT*)(dataPtr), dim,
        spacing, components, independent, scalarRange, this->GradientNormal,
        this->GradientMagnitude, this->DirectionEncoder, this));
    }
  }

  timer->StopTimer();
  // cout << "Gradients computed in " << timer->GetElapsedTime() << " seconds " << endl;
  timer->Delete();
}

int vtkFixedPointVolumeRayCastMapper::UpdateShadingTable(vtkRenderer* ren, vtkVolume* vol)
{
  if (this->ShadingRequired == 0)
  {
    return 0;
  }

  // How many components?
  int components = this->CurrentScalars->GetNumberOfComponents();

  int c;
  for (c = 0; c < ((vol->GetProperty()->GetIndependentComponents()) ? (components) : (1)); c++)
  {
    this->GradientShader->SetActiveComponent(c);
    this->GradientShader->UpdateShadingTable(ren, vol, this->GradientEstimator);

    float* r = this->GradientShader->GetRedDiffuseShadingTable(vol);
    float* g = this->GradientShader->GetGreenDiffuseShadingTable(vol);
    float* b = this->GradientShader->GetBlueDiffuseShadingTable(vol);

    float* rptr = r;
    float* gptr = g;
    float* bptr = b;

    unsigned short* tablePtr = this->DiffuseShadingTable[c];

    int i;
    for (i = 0; i < this->DirectionEncoder->GetNumberOfEncodedDirections(); i++)
    {
      *(tablePtr++) = static_cast<unsigned short>((*(rptr++)) * VTKKW_FP_SCALE + 0.5);
      *(tablePtr++) = static_cast<unsigned short>((*(gptr++)) * VTKKW_FP_SCALE + 0.5);
      *(tablePtr++) = static_cast<unsigned short>((*(bptr++)) * VTKKW_FP_SCALE + 0.5);
    }

    r = this->GradientShader->GetRedSpecularShadingTable(vol);
    g = this->GradientShader->GetGreenSpecularShadingTable(vol);
    b = this->GradientShader->GetBlueSpecularShadingTable(vol);

    rptr = r;
    gptr = g;
    bptr = b;

    tablePtr = this->SpecularShadingTable[c];

    for (i = 0; i < this->DirectionEncoder->GetNumberOfEncodedDirections(); i++)
    {
      *(tablePtr++) = static_cast<unsigned short>((*(rptr++)) * VTKKW_FP_SCALE + 0.5);
      *(tablePtr++) = static_cast<unsigned short>((*(gptr++)) * VTKKW_FP_SCALE + 0.5);
      *(tablePtr++) = static_cast<unsigned short>((*(bptr++)) * VTKKW_FP_SCALE + 0.5);
    }
  }

  return 1;
}

int vtkFixedPointVolumeRayCastMapper::UpdateGradients(vtkVolume* vol)
{
  int needToUpdate = 0;

  this->GradientOpacityRequired = 0;
  this->ShadingRequired = 0;

  // Get the image data
  vtkImageData* input = vtkImageData::SafeDownCast(this->GetInput());

  if (vol->GetProperty()->GetShade())
  {
    needToUpdate = 1;
    this->ShadingRequired = 1;
  }

  for (int c = 0; c < this->CurrentScalars->GetNumberOfComponents(); c++)
  {
    vtkPiecewiseFunction* f = vol->GetProperty()->GetGradientOpacity(c);
    if (strcmp(f->GetType(), "Constant") != 0 || f->GetValue(0.0) != 1.0)
    {
      needToUpdate = 1;
      this->GradientOpacityRequired = 1;
    }
  }

  if (!needToUpdate)
  {
    return 0;
  }

  // Check if the input has changed
  if (input == this->SavedGradientsInput && this->CurrentScalars == this->PreviousScalars &&
    input->GetMTime() < this->SavedGradientsMTime.GetMTime())
  {
    return 0;
  }

  this->ComputeGradients(vol);

  // Time to save the input used to update the tabes
  this->SavedGradientsInput = vtkImageData::SafeDownCast(this->GetInput());
  this->SavedGradientsMTime.Modified();

  return 1;
}

int vtkFixedPointVolumeRayCastMapper::UpdateColorTable(vtkVolume* vol)
{
  int needToUpdate = 0;

  // Get the image data
  vtkImageData* input = vtkImageData::SafeDownCast(this->GetInput());

  // Has the data itself changed?
  if (input != this->SavedParametersInput || this->CurrentScalars != this->PreviousScalars ||
    input->GetMTime() > this->SavedParametersMTime.GetMTime())
  {
    needToUpdate = 1;
  }

  // What is the blending mode?
  int blendMode = this->GetBlendMode();
  if (blendMode != this->SavedBlendMode)
  {
    needToUpdate = 1;

    if (this->GetBlendMode() == vtkVolumeMapper::MINIMUM_INTENSITY_BLEND)
    {
      this->FlipMIPComparison = 1;
    }
    else
    {
      this->FlipMIPComparison = 0;
    }
  }

  // How many components?
  int components = this->CurrentScalars->GetNumberOfComponents();

  // Has the sample distance changed?
  if (this->SavedSampleDistance != this->SampleDistance)
  {
    needToUpdate = 1;
  }

  vtkColorTransferFunction* rgbFunc[4];
  vtkPiecewiseFunction* grayFunc[4];
  vtkPiecewiseFunction* scalarOpacityFunc[4];
  vtkPiecewiseFunction* gradientOpacityFunc[4];
  int colorChannels[4];
  float scalarOpacityDistance[4];

  int c;

  for (c = 0; c < ((vol->GetProperty()->GetIndependentComponents()) ? (components) : (1)); c++)
  {
    colorChannels[c] = vol->GetProperty()->GetColorChannels(c);
    if (colorChannels[c] == 1)
    {
      rgbFunc[c] = nullptr;
      grayFunc[c] = vol->GetProperty()->GetGrayTransferFunction(c);
    }
    else
    {
      rgbFunc[c] = vol->GetProperty()->GetRGBTransferFunction(c);
      grayFunc[c] = nullptr;
    }
    scalarOpacityFunc[c] = vol->GetProperty()->GetScalarOpacity(c);
    gradientOpacityFunc[c] = vol->GetProperty()->GetGradientOpacity(c);
    scalarOpacityDistance[c] = vol->GetProperty()->GetScalarOpacityUnitDistance(c);

    // Has the number of color channels changed?
    if (this->SavedColorChannels[c] != colorChannels[c])
    {
      needToUpdate = 1;
    }

    // Has the color transfer function changed in some way,
    // and we are using it?
    if (colorChannels[c] == 3)
    {
      if (this->SavedRGBFunction[c] != rgbFunc[c] ||
        this->SavedParametersMTime.GetMTime() < rgbFunc[c]->GetMTime())
      {
        needToUpdate = 1;
      }
    }

    // Has the gray transfer function changed in some way,
    // and we are using it?
    if (colorChannels[c] == 1)
    {
      if (this->SavedGrayFunction[c] != grayFunc[c] ||
        this->SavedParametersMTime.GetMTime() < grayFunc[c]->GetMTime())
      {
        needToUpdate = 1;
      }
    }

    // Has the scalar opacity transfer function changed in some way?
    if (this->SavedScalarOpacityFunction[c] != scalarOpacityFunc[c] ||
      this->SavedParametersMTime.GetMTime() < scalarOpacityFunc[c]->GetMTime())
    {
      needToUpdate = 1;
    }

    // Has the gradient opacity transfer function changed in some way?
    if (this->SavedGradientOpacityFunction[c] != gradientOpacityFunc[c] ||
      this->SavedParametersMTime.GetMTime() < gradientOpacityFunc[c]->GetMTime())
    {
      needToUpdate = 1;
    }

    // Has the distance over which the scalar opacity function is defined changed?
    if (this->SavedScalarOpacityDistance[c] != scalarOpacityDistance[c])
    {
      needToUpdate = 1;
    }
  }

  // If we have not found any need to update, return now
  if (!needToUpdate)
  {
    return 0;
  }

  for (c = 0; c < ((vol->GetProperty()->GetIndependentComponents()) ? (components) : (1)); c++)
  {
    this->SavedRGBFunction[c] = rgbFunc[c];
    this->SavedGrayFunction[c] = grayFunc[c];
    this->SavedScalarOpacityFunction[c] = scalarOpacityFunc[c];
    this->SavedGradientOpacityFunction[c] = gradientOpacityFunc[c];
    this->SavedColorChannels[c] = colorChannels[c];
    this->SavedScalarOpacityDistance[c] = scalarOpacityDistance[c];
  }

  this->SavedSampleDistance = this->SampleDistance;
  this->SavedBlendMode = blendMode;
  this->SavedParametersInput = input;

  this->SavedParametersMTime.Modified();

  int scalarType = this->CurrentScalars->GetDataType();

  int i;
  float tmpArray[3 * 32768];

  // Find the scalar range
  double scalarRange[4][2];
  for (c = 0; c < components; c++)
  {
    this->CurrentScalars->GetRange(scalarRange[c], c);

    // Is the difference between max and min less than 32768? If so, and if
    // the data is not of float or double type, use a simple offset mapping.
    // If the difference between max and min is 32768 or greater, or the data
    // is of type float or double, we must use an offset / scaling mapping.
    // In this case, the array size will be 32768 - we need to figure out the
    // offset and scale factor.
    float offset;
    float scale;

    int arraySizeNeeded;

    if (scalarType == VTK_FLOAT || scalarType == VTK_DOUBLE ||
      scalarRange[c][1] - scalarRange[c][0] > 32767)
    {
      arraySizeNeeded = 32768;
      offset = -scalarRange[c][0];

      if (scalarRange[c][1] - scalarRange[c][0])
      {
        scale = 32767.0 / (scalarRange[c][1] - scalarRange[c][0]);
      }
      else
      {
        scale = 1.0;
      }
    }
    else
    {
      arraySizeNeeded = (int)(scalarRange[c][1] - scalarRange[c][0] + 1);
      offset = -scalarRange[c][0];
      scale = 1.0;
    }

    this->TableSize[c] = arraySizeNeeded;
    this->TableShift[c] = offset;
    this->TableScale[c] = scale;
  }

  if (vol->GetProperty()->GetIndependentComponents())
  {
    for (c = 0; c < components; c++)
    {
      // Sample the transfer functions between the min and max.
      if (colorChannels[c] == 1)
      {
        float tmpArray2[32768];
        grayFunc[c]->GetTable(scalarRange[c][0], scalarRange[c][1], this->TableSize[c], tmpArray2);
        for (int index = 0; index < this->TableSize[c]; index++)
        {
          tmpArray[3 * index] = tmpArray2[index];
          tmpArray[3 * index + 1] = tmpArray2[index];
          tmpArray[3 * index + 2] = tmpArray2[index];
        }
      }
      else
      {
        rgbFunc[c]->GetTable(scalarRange[c][0], scalarRange[c][1], this->TableSize[c], tmpArray);
      }
      // Convert color to short format
      for (i = 0; i < this->TableSize[c]; i++)
      {
        this->ColorTable[c][3 * i] =
          static_cast<unsigned short>(tmpArray[3 * i] * VTKKW_FP_SCALE + 0.5);
        this->ColorTable[c][3 * i + 1] =
          static_cast<unsigned short>(tmpArray[3 * i + 1] * VTKKW_FP_SCALE + 0.5);
        this->ColorTable[c][3 * i + 2] =
          static_cast<unsigned short>(tmpArray[3 * i + 2] * VTKKW_FP_SCALE + 0.5);
      }

      scalarOpacityFunc[c]->GetTable(
        scalarRange[c][0], scalarRange[c][1], this->TableSize[c], tmpArray);

      // Correct the opacity array for the spacing between the planes if we are
      // using a composite blending operation
      if (this->BlendMode == vtkVolumeMapper::COMPOSITE_BLEND)
      {
        float* ptr = tmpArray;
        double factor;
        {
          factor = this->SampleDistance / vol->GetProperty()->GetScalarOpacityUnitDistance(c);
        }

        for (i = 0; i < this->TableSize[c]; i++)
        {
          if (*ptr > 0.0001)
          {
            *ptr = 1.0 - pow((double)(1.0 - (*ptr)), factor);
          }
          ptr++;
        }
      }

      // Convert tables to short format
      for (i = 0; i < this->TableSize[c]; i++)
      {
        this->ScalarOpacityTable[c][i] =
          static_cast<unsigned short>(tmpArray[i] * VTKKW_FP_SCALE + 0.5);
      }

      if (scalarRange[c][1] - scalarRange[c][0])
      {
        gradientOpacityFunc[c]->GetTable(
          0, (scalarRange[c][1] - scalarRange[c][0]) * 0.25, 256, tmpArray);

        for (i = 0; i < 256; i++)
        {
          this->GradientOpacityTable[c][i] =
            static_cast<unsigned short>(tmpArray[i] * VTKKW_FP_SCALE + 0.5);
        }
      }
      else
      {
        for (i = 0; i < 256; i++)
        {
          this->GradientOpacityTable[c][i] = 0x0000;
        }
      }
    }
  }
  else
  {
    if (components == 2)
    {
      // Sample the transfer functions between the min and max.
      if (colorChannels[0] == 1)
      {
        float tmpArray2[32768];
        grayFunc[0]->GetTable(scalarRange[0][0], scalarRange[0][1], this->TableSize[0], tmpArray2);
        for (int index = 0; index < this->TableSize[0]; index++)
        {
          tmpArray[3 * index] = tmpArray2[index];
          tmpArray[3 * index + 1] = tmpArray2[index];
          tmpArray[3 * index + 2] = tmpArray2[index];
        }
      }
      else
      {
        rgbFunc[0]->GetTable(scalarRange[0][0], scalarRange[0][1], this->TableSize[0], tmpArray);
      }

      // Convert color to short format
      for (i = 0; i < this->TableSize[0]; i++)
      {
        this->ColorTable[0][3 * i] =
          static_cast<unsigned short>(tmpArray[3 * i] * VTKKW_FP_SCALE + 0.5);
        this->ColorTable[0][3 * i + 1] =
          static_cast<unsigned short>(tmpArray[3 * i + 1] * VTKKW_FP_SCALE + 0.5);
        this->ColorTable[0][3 * i + 2] =
          static_cast<unsigned short>(tmpArray[3 * i + 2] * VTKKW_FP_SCALE + 0.5);
      }
    }

    // The opacity table is indexed with the last component
    scalarOpacityFunc[0]->GetTable(scalarRange[components - 1][0], scalarRange[components - 1][1],
      this->TableSize[components - 1], tmpArray);

    // Correct the opacity array for the spacing between the planes if we are
    // using a composite blending operation
    if (this->BlendMode == vtkVolumeMapper::COMPOSITE_BLEND)
    {
      float* ptr = tmpArray;
      double factor = this->SampleDistance / vol->GetProperty()->GetScalarOpacityUnitDistance();
      for (i = 0; i < this->TableSize[components - 1]; i++)
      {
        if (*ptr > 0.0001)
        {
          *ptr = 1.0 - pow((double)(1.0 - (*ptr)), factor);
        }
        ptr++;
      }
    }

    // Convert tables to short format
    for (i = 0; i < this->TableSize[components - 1]; i++)
    {
      this->ScalarOpacityTable[0][i] =
        static_cast<unsigned short>(tmpArray[i] * VTKKW_FP_SCALE + 0.5);
    }

    if (scalarRange[components - 1][1] - scalarRange[components - 1][0])
    {
      gradientOpacityFunc[0]->GetTable(
        0, (scalarRange[components - 1][1] - scalarRange[components - 1][0]) * 0.25, 256, tmpArray);

      for (i = 0; i < 256; i++)
      {
        this->GradientOpacityTable[0][i] =
          static_cast<unsigned short>(tmpArray[i] * VTKKW_FP_SCALE + 0.5);
      }
    }
    else
    {
      for (i = 0; i < 256; i++)
      {
        this->GradientOpacityTable[0][i] = 0x0000;
      }
    }
  }

  return 1;
}

int vtkFixedPointVolumeRayCastMapper::ShouldUseNearestNeighborInterpolation(vtkVolume* vol)
{
  //  return ( this->UseShortCuts ||
  //           vol->GetProperty()->GetInterpolationType() == VTK_NEAREST_INTERPOLATION );
  return (vol->GetProperty()->GetInterpolationType() == VTK_NEAREST_INTERPOLATION);
}

// Print method for vtkFixedPointVolumeRayCastMapper
void vtkFixedPointVolumeRayCastMapper::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "Sample Distance: " << this->SampleDistance << endl;
  os << indent << "Interactive Sample Distance: " << this->InteractiveSampleDistance << endl;
  os << indent << "Image Sample Distance: " << this->ImageSampleDistance << endl;
  os << indent << "Minimum Image Sample Distance: " << this->MinimumImageSampleDistance << endl;
  os << indent << "Maximum Image Sample Distance: " << this->MaximumImageSampleDistance << endl;
  os << indent << "Auto Adjust Sample Distances: " << this->AutoAdjustSampleDistances << endl;
  os << indent << "LockSampleDistanceToInputSpacing: "
     << (this->LockSampleDistanceToInputSpacing ? "On\n" : "Off\n");
  os << indent << "Intermix Intersecting Geometry: "
     << (this->IntermixIntersectingGeometry ? "On\n" : "Off\n");
  os << indent << "Final Color Window: " << this->FinalColorWindow << endl;
  os << indent << "Final Color Level: " << this->FinalColorLevel << endl;
  os << indent << "Space leaping filter: " << this->SpaceLeapFilter << endl;

  // These are all things that shouldn't be printed....
  // os << indent << "ShadingRequired: " << this->ShadingRequired << endl;
  // os << indent << "GradientOpacityRequired: " << this->GradientOpacityRequired
  //   << endl;

  // os << indent << "CurrentScalars: " << this->CurrentScalars << endl;
  // os << indent << "PreviousScalars: " << this->PreviousScalars << endl;

  // if ( this->RayCastImage )
  //  {
  //  os << indent << "Ray Cast Image:\n";
  //  this->RayCastImage->PrintSelf(os,indent.GetNextIndent());
  //  }
  // else
  //  {
  //  os << indent << "Ray Cast Image: (none)\n";
  //  }

  // os << indent << "RenderWindow: " << this->RenderWindow << endl;

  // os << indent << "CompositeHelper: " << this->CompositeHelper << endl;
  // os << indent << "CompositeShadeHelper: " << this->CompositeShadeHelper << endl;
  // os << indent << "CompositeGOHelper: " << this->CompositeGOHelper << endl;
  // os << indent << "CompositeGOShadeHelper: " << this->CompositeGOShadeHelper << endl;
  // os << indent << "MIPHelper: " << this->MIPHelper << endl;

  // os << indent << "TableShift: " << this->TableShift[0] << " "
  //   << this->TableShift[1] << " " << this->TableShift[2] << " "
  //   << this->TableShift[3] << endl;
  // os << indent << "TableScale: " << this->TableScale[0] << " "
  //   << this->TableScale[1] << " " << this->TableScale[2] << " "
  //   << this->TableScale[3] << endl;

  // os << indent << "Flip Mip Comparison" << this->FlipMIPComparison << end;"
}

void vtkFixedPointVolumeRayCastMapper::ReleaseGraphicsResources(vtkWindow* win)
{
  if (win && this->ImageDisplayHelper)
  {
    this->ImageDisplayHelper->ReleaseGraphicsResources(win);
  }
}