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

  Program:   Visualization Toolkit
  Module:    vtkVolumeRayCastMapper.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  Copyright (c) 1993-2002 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 "vtkVolumeRayCastMapper.h"

#include "vtkCamera.h"
#include "vtkEncodedGradientEstimator.h"
#include "vtkEncodedGradientShader.h"
#include "vtkFiniteDifferenceGradientEstimator.h"
#include "vtkGraphicsFactory.h"
#include "vtkImageData.h"
#include "vtkMath.h"
#include "vtkMultiThreader.h"
#include "vtkObjectFactory.h"
#include "vtkPlaneCollection.h"
#include "vtkRenderWindow.h"
#include "vtkRenderer.h"
#include "vtkTimerLog.h"
#include "vtkTransform.h"
#include "vtkVolumeProperty.h"
#include "vtkVolumeRayCastFunction.h"

#include <math.h>

vtkCxxRevisionMacro(vtkVolumeRayCastMapper, "1.95");

//----------------------------------------------------------------------------
// Needed when we don't use the vtkStandardNewMacro.
vtkInstantiatorNewMacro(vtkVolumeRayCastMapper);
//----------------------------------------------------------------------------

vtkCxxSetObjectMacro(vtkVolumeRayCastMapper,VolumeRayCastFunction,
                     vtkVolumeRayCastFunction );

#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 vtkVRCMultiplyViewPointMacro( 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[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]; }

#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]

vtkVolumeRayCastMapper* vtkVolumeRayCastMapper::New()
{
  // First try to create the object from the vtkObjectFactory
  vtkObject* ret=vtkGraphicsFactory::CreateInstance("vtkVolumeRayCastMapper");
  return (vtkVolumeRayCastMapper*)ret;
}

// Construct a new vtkVolumeRayCastMapper with default values
vtkVolumeRayCastMapper::vtkVolumeRayCastMapper()
{
  this->SampleDistance             =  1.0;
  this->ImageSampleDistance        =  1.0;
  this->MinimumImageSampleDistance =  1.0;
  this->MaximumImageSampleDistance = 10.0;
  this->AutoAdjustSampleDistances  =  1;
  this->VolumeRayCastFunction      = NULL;
  
  this->GradientEstimator  = vtkFiniteDifferenceGradientEstimator::New();
  this->GradientShader     = vtkEncodedGradientShader::New();
  
  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->ImageMemorySize[0]     = 0;
  this->ImageMemorySize[1]     = 0;

  this->Threader               = vtkMultiThreader::New();
  this->NumberOfThreads        = this->Threader->GetNumberOfThreads();

  this->Image                  = NULL;
  this->RowBounds              = NULL;
  this->OldRowBounds           = NULL;

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

  this->ZBuffer                = NULL;
  this->ZBufferSize[0]         = 0;
  this->ZBufferSize[1]         = 0;
  this->ZBufferOrigin[0]       = 0;
  this->ZBufferOrigin[1]       = 0;
  
  this->IntermixIntersectingGeometry = 1;
  
}

// Destruct a vtkVolumeRayCastMapper - clean up any memory used
vtkVolumeRayCastMapper::~vtkVolumeRayCastMapper()
{
  if ( this->GradientEstimator )
    {
    this->GradientEstimator->UnRegister(this);
    this->GradientEstimator = NULL;
    }

  this->GradientShader->Delete();

  this->SetVolumeRayCastFunction(NULL);
  
  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();
  
  if ( this->Image )
    {
    delete [] this->Image;
    }
  
  if ( this->RenderTableSize )
    {
    delete [] this->RenderTimeTable;
    delete [] this->RenderVolumeTable;
    delete [] this->RenderRendererTable;
    }
  
  if ( this->RowBounds )
    {
    delete [] this->RowBounds;
    delete [] this->OldRowBounds;
    }
}

float vtkVolumeRayCastMapper::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;
}

void vtkVolumeRayCastMapper::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 vtkVolumeRayCastMapper::SetNumberOfThreads( int num )
{
  this->Threader->SetNumberOfThreads( num );
}

void vtkVolumeRayCastMapper::SetGradientEstimator( 
                                       vtkEncodedGradientEstimator *gradest )
{

  // If we are setting it to its current value, don't do anything
  if ( this->GradientEstimator == gradest )
    {
    return;
    }

  // If we already have a gradient estimator, unregister it.
  if ( this->GradientEstimator )
    {
    this->GradientEstimator->UnRegister(this);
    this->GradientEstimator = NULL;
    }

  // If we are passing in a non-NULL estimator, register it
  if ( gradest )
    {
    gradest->Register( this );
    }

  // Actually set the estimator, and consider the object Modified
  this->GradientEstimator = gradest;
  this->Modified();
}

float vtkVolumeRayCastMapper::GetGradientMagnitudeScale()
{
  if ( !this->GradientEstimator ) 
    {
    vtkErrorMacro( "You must have a gradient estimator set to get the scale" );
    return 1.0;
    }
  
  return this->GradientEstimator->GetGradientMagnitudeScale();  
}

float vtkVolumeRayCastMapper::GetGradientMagnitudeBias()
{
  if ( !this->GradientEstimator ) 
    {
    vtkErrorMacro( "You must have a gradient estimator set to get the bias" );
    return 1.0;
    }
  
  return this->GradientEstimator->GetGradientMagnitudeBias();  
}

void vtkVolumeRayCastMapper::ReleaseGraphicsResources(vtkWindow *)
{
}

void vtkVolumeRayCastMapper::Render( vtkRenderer *ren, vtkVolume *vol )
{
  // make sure that we have scalar input and update the scalar input
  if ( this->GetInput() == NULL ) 
    {
    vtkErrorMacro(<< "No Input!");
    return;
    }
  else
    {
    this->GetInput()->UpdateInformation();
    this->GetInput()->SetUpdateExtentToWholeExtent();
    this->GetInput()->Update();
    } 


  int scalarType = this->GetInput()->GetPointData()->GetScalars()->GetDataType();
  if (scalarType != VTK_UNSIGNED_SHORT && scalarType != VTK_UNSIGNED_CHAR)
    {
    vtkErrorMacro ("Cannot volume render data of type " 
                   << vtkImageScalarTypeNameMacro(scalarType) 
                   << ", only unsigned char or unsigned short.");
    return;
    }
  // Start timing now. We didn't want to capture the update of the
  // input data in the times
  this->Timer->StartTimer();
  
  this->ConvertCroppingRegionPlanesToVoxels();
  
  this->UpdateShadingTables( ren, vol );

  // This is the input of this mapper
  vtkImageData *input = this->GetInput();
  
  // 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();
  float *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->GetPerspectiveTransformMatrix(aspect[0]/aspect[1], 
                                                        0.0, 1.0 ));
  this->PerspectiveTransform->Concatenate(cam->GetViewTransformMatrix());
  this->PerspectiveMatrix->DeepCopy(this->PerspectiveTransform->GetMatrix());

  // Compute some matrices from voxels to view and vice versa based 
  // on the whole input
  this->ComputeMatrices( input, vol );


  // How big is the viewport in pixels?
  float *viewport   =  ren->GetViewport();
  int *renWinSize   =  ren->GetRenderWindow()->GetSize();

  // Save this so that we can restore it if the image is cancelled
  float oldImageSampleDistance = this->ImageSampleDistance;
  
  // 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 )
    {
    float oldTime = this->RetrieveRenderTime( ren, vol );
    float newTime = vol->GetAllocatedRenderTime();
    this->ImageSampleDistance *= sqrt(oldTime / newTime);
    this->ImageSampleDistance = 
      (this->ImageSampleDistance>this->MaximumImageSampleDistance)?
      (this->MaximumImageSampleDistance):(this->ImageSampleDistance);
    this->ImageSampleDistance = 
      (this->ImageSampleDistance<this->MinimumImageSampleDistance)?
      (this->MinimumImageSampleDistance):(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->ImageViewportSize[0] = 
    static_cast<int>(width/this->ImageSampleDistance);
  this->ImageViewportSize[1] = 
    static_cast<int>(height/this->ImageSampleDistance);
  
  // 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 ( this->ComputeRowBounds( vol, ren ) )
    {
    vtkVolumeRayCastStaticInfo *staticInfo = new vtkVolumeRayCastStaticInfo;
    staticInfo->ClippingPlane = NULL;
    staticInfo->Volume = vol;
    staticInfo->Renderer = ren;
    staticInfo->ScalarDataPointer = 
      this->GetInput()->GetPointData()->GetScalars()->GetVoidPointer(0);
    staticInfo->ScalarDataType = 
      this->GetInput()->GetPointData()->GetScalars()->GetDataType();

    // 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 this->ImageOrigin into (x1,y1) in window (not viewport!)
      // coordinates. 
      x1 = static_cast<int> (
        viewport[0] * static_cast<float>(renWinSize[0]) +
        static_cast<float>(this->ImageOrigin[0]) * this->ImageSampleDistance );
      y1 = static_cast<int> (
        viewport[1] * static_cast<float>(renWinSize[1]) +
        static_cast<float>(this->ImageOrigin[1]) * this->ImageSampleDistance);
      
      // compute z buffer size
      this->ZBufferSize[0] = static_cast<int>(
        static_cast<float>(this->ImageInUseSize[0]) * this->ImageSampleDistance);
      this->ZBufferSize[1] = static_cast<int>(
        static_cast<float>(this->ImageInUseSize[1]) * this->ImageSampleDistance);
      
      // Use the size to compute (x2,y2) in window coordinates
      x2 = x1 + this->ZBufferSize[0] - 1;
      y2 = y1 + this->ZBufferSize[1] - 1;
      
      // This is the z buffer origin (in viewport coordinates)
      this->ZBufferOrigin[0] = static_cast<int>(
        static_cast<float>(this->ImageOrigin[0]) * this->ImageSampleDistance);
      this->ZBufferOrigin[1] = static_cast<int>(
        static_cast<float>(this->ImageOrigin[1]) * this->ImageSampleDistance);
      
      // Capture the z buffer
      this->ZBuffer = ren->GetRenderWindow()->GetZbufferData(x1,y1,x2,y2);
      }
    
    // This must be done before FunctionInitialize since FunctionInitialize
    // depends on the gradient opacity constant (computed in here) to 
    // determine whether to save the gradient magnitudes
    vol->UpdateTransferFunctions( ren );
    
    // Requires UpdateTransferFunctions to have been called first
    this->VolumeRayCastFunction->FunctionInitialize( ren, vol, 
                                                     staticInfo );
    
    vol->UpdateScalarOpacityforSampleSize( ren, this->SampleDistance );

    staticInfo->CameraThickness = 
      static_cast<float>(ren->GetActiveCamera()->GetThickness());
    
    // Copy the viewToVoxels matrix to 16 floats
    int i, j;
    for ( j = 0; j < 4; j++ )
      {
      for ( i = 0; i < 4; i++ )
        {
        staticInfo->ViewToVoxelsMatrix[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++ )
        {
        staticInfo->WorldToVoxelsMatrix[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++ )
        {
        staticInfo->VoxelsToWorldMatrix[j*4+i] = 
          static_cast<float>(this->VoxelsToWorldMatrix->GetElement(j,i));
        }
      }

    if ( this->ClippingPlanes )
      {
      this->InitializeClippingPlanes( staticInfo, this->ClippingPlanes );
      }
    else
      {
      staticInfo->NumberOfClippingPlanes = 0;
      }
    

    // Copy in the image info
    staticInfo->ImageInUseSize[0]    = this->ImageInUseSize[0];
    staticInfo->ImageInUseSize[1]    = this->ImageInUseSize[1];
    staticInfo->ImageMemorySize[0]   = this->ImageMemorySize[0];
    staticInfo->ImageMemorySize[1]   = this->ImageMemorySize[1];
    staticInfo->ImageViewportSize[0] = this->ImageViewportSize[0];
    staticInfo->ImageViewportSize[1] = this->ImageViewportSize[1];
    
    staticInfo->ImageOrigin[0] = this->ImageOrigin[0];
    staticInfo->ImageOrigin[1] = this->ImageOrigin[1];
    
    staticInfo->Image     = this->Image;
    staticInfo->RowBounds = this->RowBounds;
    
    // Set the number of threads to use for ray casting,
    // then set the execution method and do it.
    this->Threader->SetNumberOfThreads( this->NumberOfThreads );
    this->Threader->SetSingleMethod( VolumeRayCastMapper_CastRays, 
                                     (void *)staticInfo);
    this->Threader->SingleMethodExecute();

    if ( !ren->GetRenderWindow()->GetAbortRender() )
      {
      this->RenderTexture( vol, ren );
      this->Timer->StopTimer();
      this->TimeToDraw = this->Timer->GetElapsedTime();
      this->StoreRenderTime( ren, vol, this->TimeToDraw );
      }
    // Restore the image sample distance so that automatic adjustment
    // will work correctly
    else
      {
      this->ImageSampleDistance = oldImageSampleDistance;
      }
    
    if ( staticInfo->ClippingPlane )
      {
      delete staticInfo->ClippingPlane;
      }
    delete staticInfo;

    if ( this->ZBuffer )
      {
      delete [] this->ZBuffer;
      this->ZBuffer = NULL;
      }
    
    }  
}
VTK_THREAD_RETURN_TYPE VolumeRayCastMapper_CastRays( void *arg )
{
  // Get the info out of the input structure
  int threadID    = ((vtkMultiThreader::ThreadInfo *)(arg))->ThreadID;
  int threadCount = ((vtkMultiThreader::ThreadInfo *)(arg))->NumberOfThreads;
  vtkVolumeRayCastStaticInfo *staticInfo  = 
    (vtkVolumeRayCastStaticInfo *)((vtkMultiThreader::ThreadInfo *)arg)->UserData;
  
  int i, j, k;
  unsigned char *ucptr;

  vtkVolumeRayCastMapper *me = 
    vtkVolumeRayCastMapper::SafeDownCast( staticInfo->Volume->GetMapper() );
  
  if ( !me )
    {
    vtkGenericWarningMacro("The volume does not have a ray cast mapper!");
    return VTK_THREAD_RETURN_VALUE;
    }
  
  vtkVolumeRayCastDynamicInfo *dynamicInfo = new vtkVolumeRayCastDynamicInfo;

  // Initialize this to avoid purify problems
  dynamicInfo->ScalarValue = 0;
  
  float *rayStart     = dynamicInfo->TransformedStart;
  float *rayEnd       = dynamicInfo->TransformedEnd;
  float *rayDirection = dynamicInfo->TransformedDirection;
  float *rayStep      = dynamicInfo->TransformedIncrement;
  
  float norm;
  float viewRay[3];
  float rayCenter[3];
  float absStep[3];
  float voxelPoint[4];
 
  // We need to know what the center ray is (in voxel coordinates) to
  // compute sampling distance later on. Save it in an instance variable.
  
  // This is tbe near end of the center ray in view coordinates
  // Convert it to voxel coordinates
  viewRay[0] = viewRay[1] = viewRay[2] = 0.0;  
  vtkVRCMultiplyPointMacro( viewRay, rayStart,
                            staticInfo->ViewToVoxelsMatrix );
  
  // This is the far end of the center ray in view coordinates
  // Convert it to voxel coordiantes
  viewRay[2] = 1.0;  
  vtkVRCMultiplyPointMacro( viewRay, voxelPoint,
                            staticInfo->ViewToVoxelsMatrix );

  // Turn these two points into a vector
  rayCenter[0] = voxelPoint[0] - rayStart[0];
  rayCenter[1] = voxelPoint[1] - rayStart[1];
  rayCenter[2] = voxelPoint[2] - rayStart[2];
  
  // normalize the vector based on world coordinate distance
  // This way we can scale by sample distance and it will work
  // out even though we are in voxel coordinates
  rayCenter[0] /= staticInfo->CameraThickness;
  rayCenter[1] /= staticInfo->CameraThickness;
  rayCenter[2] /= staticInfo->CameraThickness;

  float centerScale = sqrt( (rayCenter[0] * rayCenter[0]) +
                            (rayCenter[1] * rayCenter[1]) +
                            (rayCenter[2] * rayCenter[2]) );
  
  rayCenter[0] /= centerScale;
  rayCenter[1] /= centerScale;
  rayCenter[2] /= centerScale;
  
  float bounds[6];
  int dim[3];
  
  me->GetInput()->GetDimensions(dim);
  bounds[0] = bounds[2] = bounds[4] = 0.0;
  bounds[1] = dim[0]-1;
  bounds[3] = dim[1]-1;
  bounds[5] = dim[2]-1;
               
  // If we have a simple crop box then we can tighten the bounds
  if ( me->Cropping && me->CroppingRegionFlags == 0x2000 )
    {
    bounds[0] = me->VoxelCroppingRegionPlanes[0];
    bounds[1] = me->VoxelCroppingRegionPlanes[1];
    bounds[2] = me->VoxelCroppingRegionPlanes[2];
    bounds[3] = me->VoxelCroppingRegionPlanes[3];
    bounds[4] = me->VoxelCroppingRegionPlanes[4];
    bounds[5] = me->VoxelCroppingRegionPlanes[5];
    }

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

  int *imageInUseSize     = staticInfo->ImageInUseSize;
  int *imageMemorySize    = staticInfo->ImageMemorySize;
  int *imageViewportSize  = staticInfo->ImageViewportSize;
  int *imageOrigin        = staticInfo->ImageOrigin;
  int *rowBounds          = staticInfo->RowBounds;
  
  unsigned char *imagePtr = staticInfo->Image;
  
  float sampleDistance    = me->GetSampleDistance();

  float val;
  
  vtkRenderWindow *renWin = staticInfo->Renderer->GetRenderWindow();

  // Compute the offset valuex for viewing rays - this is the 1 / fullSize
  // value to add to the computed location so that they falls between
  // -1 + 1/fullSize and 1 - 1/fullSize and are each 2/fullSize apart.
  // fullSize is the viewport size along the corresponding direction (in 
  // pixels)
  float offsetX = 1.0 / static_cast<float>(imageViewportSize[0]);
  float offsetY = 1.0 / static_cast<float>(imageViewportSize[1]);

  // Some variables needed for non-subvolume cropping
  float fullRayStart[3];
  float fullRayEnd[3];
  float fullRayDirection[3];
  int bitLoop, bitFlag;
  float rgbaArray[40], distanceArray[10], scalarArray[10];
  float tmp, tmpArray[4];
  int arrayCount;
  
  for ( j = 0; j < imageInUseSize[1]; j++ )
    {
    if ( j%threadCount != threadID )
      {
      continue;
      }

    if ( !threadID )
      {
      if ( renWin->CheckAbortStatus() )
        {
        break;
        }
      }
    else if ( renWin->GetAbortRender() )
      {
      break;
      }
    
    ucptr = imagePtr + 4*(j*imageMemorySize[0] + 
                             rowBounds[j*2]);

    // 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 ass offsetX (which is 1/fullSize) to 
    // center it.
    viewRay[1] = ((static_cast<float>(j) + static_cast<float>(imageOrigin[1])) /
                  imageViewportSize[1]) * 2.0 - 1.0 + offsetY;

    for ( i = rowBounds[j*2]; i <= rowBounds[j*2+1]; i++ )
      {
      // Initialize for the cases where the ray doesn't intersect anything
      ucptr[0] = 0;
      ucptr[1] = 0;
      ucptr[2] = 0;
      ucptr[3] = 0;
      
      // Convert the view coordinates for the start and end of the
      // ray into voxel coordinates, then compute the origin and direction.

      // 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 ass offsetX (which is 1/fullSize) to 
      // center it.
      viewRay[0] = ((static_cast<float>(i) + 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,
                                staticInfo->ViewToVoxelsMatrix );

      viewRay[2] = 
        (me->ZBuffer)?(me->GetZBufferValue(i,j)):(1.0);
      vtkVRCMultiplyPointMacro( viewRay, rayEnd,
                                staticInfo->ViewToVoxelsMatrix );

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

      // If cropping is off, or we are just doing a subvolume, we can
      // do the easy thing here
      if ( !me->Cropping || me->CroppingRegionFlags == 0x2000 )
        {
        if ( me->ClipRayAgainstVolume( dynamicInfo, bounds ) &&
             ( staticInfo->NumberOfClippingPlanes == 0 ||
               me->ClipRayAgainstClippingPlanes( dynamicInfo, staticInfo ) ) )
          { 
          rayDirection[0] = rayEnd[0] - rayStart[0];
          rayDirection[1] = rayEnd[1] - rayStart[1];
          rayDirection[2] = rayEnd[2] - rayStart[2];
          
          // Find the length of the input ray. It is not normalized
          norm = sqrt( rayDirection[0] * rayDirection[0] +
                       rayDirection[1] * rayDirection[1] +
                       rayDirection[2] * rayDirection[2] );
          
          // Normalize this ray into rayStep
          rayStep[0] = rayDirection[0] / norm;
          rayStep[1] = rayDirection[1] / norm;
          rayStep[2] = rayDirection[2] / norm;
          
          // Correct for perspective in the sample distance. 1.0 over the
          // dot product between this ray and the center ray is the
          // correction factor to allow samples to be taken on parallel
          // planes rather than concentric hemispheres. This factor will
          // compute to be 1.0 for parallel.
          val =  ( rayStep[0] * rayCenter[0] +
                   rayStep[1] * rayCenter[1] +
                   rayStep[2] * rayCenter[2] );
          norm = (val != 0)?(1.0/val):(1.0);
          
          // Now multiple the normalized step by the sample distance and this
          // correction factor to find the actual step
          rayStep[0] *= norm * sampleDistance * centerScale;
          rayStep[1] *= norm * sampleDistance * centerScale;
          rayStep[2] *= norm * sampleDistance * centerScale;
          
          absStep[0] = ( rayStep[0] < 0.0 )?(-rayStep[0]):(rayStep[0]);
          absStep[1] = ( rayStep[1] < 0.0 )?(-rayStep[1]):(rayStep[1]);
          absStep[2] = ( rayStep[2] < 0.0 )?(-rayStep[2]):(rayStep[2]);
          
          if ( absStep[0] >= absStep[1] && absStep[0] >= absStep[2] )
            {
            dynamicInfo->NumberOfStepsToTake = static_cast<int>
              ((rayEnd[0]-rayStart[0]) / rayStep[0]);
            }
          else if ( absStep[1] >= absStep[2] && absStep[1] >= absStep[0] )
            {
            dynamicInfo->NumberOfStepsToTake = static_cast<int>
              ((rayEnd[1]-rayStart[1]) / rayStep[1]);
            }
          else
            {
            dynamicInfo->NumberOfStepsToTake = static_cast<int>
              ((rayEnd[2]-rayStart[2]) / rayStep[2]);
            }
          
          me->VolumeRayCastFunction->CastRay( dynamicInfo, staticInfo );
          if ( dynamicInfo->Color[3] > 0.0 )
            {
            val = (dynamicInfo->Color[0]/dynamicInfo->Color[3])*255.0;
            val = (val > 255.0)?(255.0):(val);
            val = (val <   0.0)?(  0.0):(val);
            ucptr[0] = static_cast<unsigned char>(val);
            
            val = (dynamicInfo->Color[1]/dynamicInfo->Color[3])*255.0;
            val = (val > 255.0)?(255.0):(val);
            val = (val <   0.0)?(  0.0):(val);
            ucptr[1] = static_cast<unsigned char>(val);
            
            val = (dynamicInfo->Color[2]/dynamicInfo->Color[3])*255.0;
            val = (val > 255.0)?(255.0):(val);
            val = (val <   0.0)?(  0.0):(val);
            ucptr[2] = static_cast<unsigned char>(val);
            
            val = dynamicInfo->Color[3]*255.0;
            val = (val > 255.0)?(255.0):(val);
            val = (val <   0.0)?(  0.0):(val);
            ucptr[3] = static_cast<unsigned char>(val);
            }
          }
        }
      // Otherwise, cropping is on and we don't have a simple subvolume.
      // We'll have to cast a ray for each of the 27 regions that is on 
      // and composite the results.
      else
        {
        // We'll keep an array of regions that we intersect, the arrayCount
        // variable will count how many of them we have
        arrayCount = 0;

        // Save the ray start, end, and direction. We will modify this
        // during each iteration of the loop for the current cropping region.
        fullRayStart[0] = rayStart[0];
        fullRayStart[1] = rayStart[1];
        fullRayStart[2] = rayStart[2];
          
        fullRayEnd[0] = rayEnd[0];
        fullRayEnd[1] = rayEnd[1];
        fullRayEnd[2] = rayEnd[2];

        fullRayDirection[0] = rayDirection[0];
        fullRayDirection[1] = rayDirection[1];
        fullRayDirection[2] = rayDirection[2];

        // Loop through the twenty seven cropping regions
        bitFlag = 1;
        for ( bitLoop = 0; bitLoop < 27; bitLoop++ )
          {
          // Check if this cropping region is on
          if ( !(me->CroppingRegionFlags & bitFlag) )
            {
            bitFlag = bitFlag << 1;
            continue;
            }

          // Restore the ray information
          rayStart[0] = fullRayStart[0];
          rayStart[1] = fullRayStart[1];
          rayStart[2] = fullRayStart[2];
          
          rayEnd[0] = fullRayEnd[0];
          rayEnd[1] = fullRayEnd[1];
          rayEnd[2] = fullRayEnd[2];

          rayDirection[0] = fullRayDirection[0];
          rayDirection[1] = fullRayDirection[1];
          rayDirection[2] = fullRayDirection[2];  

          // Figure out the bounds of the cropping region
          // along the X axis
          switch ( bitLoop % 3 )
            {
            case 0:
              bounds[0] = 0;
              bounds[1] = me->VoxelCroppingRegionPlanes[0];
              break;
            case 1:
              bounds[0] = me->VoxelCroppingRegionPlanes[0];
              bounds[1] = me->VoxelCroppingRegionPlanes[1];
              break;
            case 2:
              bounds[0] = me->VoxelCroppingRegionPlanes[1];
              bounds[1] = staticInfo->DataSize[0] - 1;
              break;
            }
          
          // Figure out the bounds of the cropping region
          // along the Y axis
          switch ( (bitLoop % 9) / 3 )
            {
            case 0:
              bounds[2] = 0;
              bounds[3] = me->VoxelCroppingRegionPlanes[2];
              break;
            case 1:
              bounds[2] = me->VoxelCroppingRegionPlanes[2];
              bounds[3] = me->VoxelCroppingRegionPlanes[3];
              break;
            case 2:
              bounds[2] = me->VoxelCroppingRegionPlanes[3];
              bounds[3] = staticInfo->DataSize[1] - 1;
              break;
            }
          
          // Figure out the bounds of the cropping region
          // along the Z axis
          switch ( bitLoop / 9 )
            {
            case 0:
              bounds[4] = 0;
              bounds[5] = me->VoxelCroppingRegionPlanes[4];
              break;
            case 1:
              bounds[4] = me->VoxelCroppingRegionPlanes[4];
              bounds[5] = me->VoxelCroppingRegionPlanes[5];
              break;
            case 2:
              bounds[4] = me->VoxelCroppingRegionPlanes[5];
              bounds[5] = staticInfo->DataSize[2] - 1;
              break;
            }
          
          // Check against the bounds of the volume
          for ( k = 0; k < 3; k++ )
            {
            if ( bounds[2*k] < 0 )
              {
              bounds[2*k] = 0;
              }
            if ( bounds[2*k + 1] > (staticInfo->DataSize[k]-1) )
              {
              bounds[2*k + 1] = (staticInfo->DataSize[k] - 1);
              }
            }

          // Clip against the volume and the clipping planes
          if ( me->ClipRayAgainstVolume( dynamicInfo, bounds ) &&
               ( staticInfo->NumberOfClippingPlanes == 0 ||
                 me->ClipRayAgainstClippingPlanes( dynamicInfo, 
                                                   staticInfo ) ) )
            {
            // The ray start and end may have changed - recompute
            // the direction
            rayDirection[0] = rayEnd[0] - rayStart[0];
            rayDirection[1] = rayEnd[1] - rayStart[1];
            rayDirection[2] = rayEnd[2] - rayStart[2];

            // Find the length of the ray. It is not normalized yet
            norm = sqrt( rayDirection[0] * rayDirection[0] +
                         rayDirection[1] * rayDirection[1] +
                         rayDirection[2] * rayDirection[2] );
          
            // Normalize this ray into rayStep
            rayStep[0] = rayDirection[0] / norm;
            rayStep[1] = rayDirection[1] / norm;
            rayStep[2] = rayDirection[2] / norm;
          
            // Correct for perspective in the sample distance. 1.0 over the
            // dot product between this ray and the center ray is the
            // correction factor to allow samples to be taken on parallel
            // planes rather than concentric hemispheres. This factor will
            // compute to be 1.0 for parallel.
            val =  ( rayStep[0] * rayCenter[0] +
                     rayStep[1] * rayCenter[1] +
                     rayStep[2] * rayCenter[2] );
            norm = (val != 0)?(1.0/val):(1.0);
            
            // Now multiple the normalized step by the sample distance and this
            // correction factor to find the actual step
            rayStep[0] *= norm * sampleDistance * centerScale;
            rayStep[1] *= norm * sampleDistance * centerScale;
            rayStep[2] *= norm * sampleDistance * centerScale;
          
            // Find the major direction to determine the number of
            // steps to take
            absStep[0] = ( rayStep[0] < 0.0 )?(-rayStep[0]):(rayStep[0]);
            absStep[1] = ( rayStep[1] < 0.0 )?(-rayStep[1]):(rayStep[1]);
            absStep[2] = ( rayStep[2] < 0.0 )?(-rayStep[2]):(rayStep[2]);
            if ( absStep[0] >= absStep[1] && absStep[0] >= absStep[2] )
              {
              dynamicInfo->NumberOfStepsToTake = static_cast<int>
                ((rayEnd[0]-rayStart[0]) / rayStep[0]);
              }
            else if ( absStep[1] >= absStep[2] && absStep[1] >= absStep[0] )
              {
              dynamicInfo->NumberOfStepsToTake = static_cast<int>
                ((rayEnd[1]-rayStart[1]) / rayStep[1]);
              }
            else
              {
              dynamicInfo->NumberOfStepsToTake = static_cast<int>
                ((rayEnd[2]-rayStart[2]) / rayStep[2]);
              }
            
            // Cast the ray
            me->VolumeRayCastFunction->CastRay( dynamicInfo, staticInfo );
            
            // If the ray returns a non-transparent color, store this
            // in our arrays of distances and colors
            if ( dynamicInfo->Color[3] > 0.0 )
              {
              // Figure out the distance from this ray start to the full ray
              // start and use this to sort the ray segments
              for ( k = 0; k < 3; k++ )
                {
                if ( absStep[k] >= absStep[(k+1)%3] && 
                     absStep[k] >= absStep[(k+2)%3] )
                  {
                  distanceArray[arrayCount] = 
                    (rayStart[k] - fullRayStart[k]) / rayStep[k];
                  break;
                  }
                }
              
              // Store the ray color
              rgbaArray[4*arrayCount  ] = dynamicInfo->Color[0];
              rgbaArray[4*arrayCount+1] = dynamicInfo->Color[1];
              rgbaArray[4*arrayCount+2] = dynamicInfo->Color[2];
              rgbaArray[4*arrayCount+3] = dynamicInfo->Color[3];
              scalarArray[arrayCount] = dynamicInfo->ScalarValue;

              if ( !staticInfo->MIPFunction )
                {
                // Do a sort pass (one iteration of bubble sort each time an
                // element is added. The array stores element from farthest
                // to closest
                for ( k = arrayCount; 
                      k > 0 && distanceArray[k] > distanceArray[k-1]; k-- )
                  {
                  tmp = distanceArray[k];
                  distanceArray[k] = distanceArray[k-1];
                  distanceArray[k-1] = tmp;
                  
                  tmpArray[0] = rgbaArray[4*k  ];
                  tmpArray[1] = rgbaArray[4*k+1];
                  tmpArray[2] = rgbaArray[4*k+2];
                  tmpArray[3] = rgbaArray[4*k+3];
                  
                  rgbaArray[4*k  ] = rgbaArray[4*(k-1)  ];
                  rgbaArray[4*k+1] = rgbaArray[4*(k-1)+1];
                  rgbaArray[4*k+2] = rgbaArray[4*(k-1)+2];
                  rgbaArray[4*k+3] = rgbaArray[4*(k-1)+3];
                  
                  rgbaArray[4*(k-1)  ] = tmpArray[0];
                  rgbaArray[4*(k-1)+1] = tmpArray[1];
                  rgbaArray[4*(k-1)+2] = tmpArray[2];
                  rgbaArray[4*(k-1)+3] = tmpArray[3];
                  }
                }
              
              arrayCount++;
              }
            }
          
          // Move the bit over by one
          bitFlag = bitFlag << 1;
          }

        // We have encountered something in at least one crop region - 
        // merge all results into one RGBA value
        if ( arrayCount )
          {
          // Is this MIP compositing? We need to treat this differently
          if ( staticInfo->MIPFunction )
            {
            dynamicInfo->Color[0] = 0.0;
            dynamicInfo->Color[1] = 0.0;
            dynamicInfo->Color[2] = 0.0;
            dynamicInfo->Color[3] = 0.0;
            dynamicInfo->ScalarValue = 0.0;
            
            // If we are maximizing the opacity, find the max Color[3]
            if ( staticInfo->MaximizeOpacity )
              {
              for ( k = 0; k < arrayCount; k++ )
                {
                if ( rgbaArray[k*4+3] > dynamicInfo->Color[3] )
                  {
                  dynamicInfo->Color[0] = rgbaArray[k*4  ];
                  dynamicInfo->Color[1] = rgbaArray[k*4+1];
                  dynamicInfo->Color[2] = rgbaArray[k*4+2];
                  dynamicInfo->Color[3] = rgbaArray[k*4+3];                  
                  }                
                }
              }
            // Otherwise we are maximizing scalar value
            else
              {
              for ( k = 0; k < arrayCount; k++ )
                {
                if ( scalarArray[k] > dynamicInfo->ScalarValue )
                  {
                  dynamicInfo->Color[0] = rgbaArray[k*4  ];
                  dynamicInfo->Color[1] = rgbaArray[k*4+1];
                  dynamicInfo->Color[2] = rgbaArray[k*4+2];
                  dynamicInfo->Color[3] = rgbaArray[k*4+3];         
                  dynamicInfo->ScalarValue = scalarArray[k];
                  }                
                }
              }
            }
          else
            {
            // Now we have the sorted distances / colors, put them together
            // in a back-to-front order. First, initialize the color to black
            // and the remaining opacity (color[3]) to 1.0
            dynamicInfo->Color[0] = 0.0;
            dynamicInfo->Color[1] = 0.0;
            dynamicInfo->Color[2] = 0.0;
            dynamicInfo->Color[3] = 1.0;
          
           // Now do alpha blending, keeping remaining opacity in color[3]
            for ( k = 0; k < arrayCount; k++ )
              {
              dynamicInfo->Color[0] = dynamicInfo->Color[0] * 
                (1.0 - rgbaArray[k*4 + 3]) + rgbaArray[k*4 + 0];
              dynamicInfo->Color[1] = dynamicInfo->Color[1] * 
                (1.0 - rgbaArray[k*4 + 3]) + rgbaArray[k*4 + 1];
              dynamicInfo->Color[2] = dynamicInfo->Color[2] * 
                (1.0 - rgbaArray[k*4 + 3]) + rgbaArray[k*4 + 2];
              dynamicInfo->Color[3] *= 1.0 - rgbaArray[k*4 + 3];
              }
          
            // Take 1.0 - color[3] to convert from remaining opacity to alpha.
            dynamicInfo->Color[3] = 1.0 - dynamicInfo->Color[3];
            }
          
          val = (dynamicInfo->Color[0]/dynamicInfo->Color[3])*255.0;
          val = (val > 255.0)?(255.0):(val);
          val = (val <   0.0)?(  0.0):(val);
          ucptr[0] = static_cast<unsigned char>(val);
          
          val = (dynamicInfo->Color[1]/dynamicInfo->Color[3])*255.0;
          val = (val > 255.0)?(255.0):(val);
          val = (val <   0.0)?(  0.0):(val);
          ucptr[1] = static_cast<unsigned char>(val);
          
          val = (dynamicInfo->Color[2]/dynamicInfo->Color[3])*255.0;
          val = (val > 255.0)?(255.0):(val);
          val = (val <   0.0)?(  0.0):(val);
          ucptr[2] = static_cast<unsigned char>(val);
          
          val = dynamicInfo->Color[3]*255.0;
          val = (val > 255.0)?(255.0):(val);
          val = (val <   0.0)?(  0.0):(val);
          ucptr[3] = static_cast<unsigned char>(val);
          }
        }
      
      // Increment the image pointer
      ucptr+=4;
      }
    }

  delete dynamicInfo;
  
  return VTK_THREAD_RETURN_VALUE;
}

float vtkVolumeRayCastMapper::GetZBufferValue(int x, int y)
{
  int xPos, yPos;
  
  xPos = static_cast<int>(static_cast<float>(x) * this->ImageSampleDistance);
  yPos = static_cast<int>(static_cast<float>(y) * this->ImageSampleDistance);
  
  xPos = (xPos >= this->ZBufferSize[0])?(this->ZBufferSize[0]-1):(xPos);
  yPos = (yPos >= this->ZBufferSize[1])?(this->ZBufferSize[1]-1):(yPos);
  
  return *(this->ZBuffer + yPos*this->ZBufferSize[0] + xPos);
}

int vtkVolumeRayCastMapper::ComputeRowBounds(vtkVolume   *vol,
                                             vtkRenderer *ren)
{
  float voxelPoint[3];
  float viewPoint[8][4];
  int i, j, k;
  unsigned char *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];
  
  this->GetInput()->GetDimensions(dim);
  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);
  
  float camPos[3];
  float worldBounds[6];
  vol->GetBounds( worldBounds );
  int insideFlag = 0;
  ren->GetActiveCamera()->GetPosition( camPos );
  if ( camPos[0] >= worldBounds[0] &&
       camPos[0] <= worldBounds[1] &&
       camPos[1] >= worldBounds[2] &&
       camPos[1] <= worldBounds[3] &&
       camPos[2] >= worldBounds[4] &&
       camPos[2] <= worldBounds[5] )
    {
    insideFlag = 1;
    }
  
  
  // 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];
    }
  
  
  // 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));
  
  // 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 * this->ImageViewportSize[0] - 2; 
  minY = ( minY + 1.0 ) * 0.5 * this->ImageViewportSize[1] - 2; 
  maxX = ( maxX + 1.0 ) * 0.5 * this->ImageViewportSize[0] + 2; 
  maxY = ( maxY + 1.0 ) * 0.5 * this->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 > this->ImageViewportSize[0]-1 &&
         maxX > this->ImageViewportSize[0]-1 ) ||
       ( minX > this->ImageViewportSize[0]-1 &&
         maxX > this->ImageViewportSize[0]-1 ) )
    {
    return 0;
    }

  int oldImageMemorySize[2];
  oldImageMemorySize[0] = this->ImageMemorySize[0];
  oldImageMemorySize[1] = this->ImageMemorySize[1];
  
  // Swap the row bounds
  int *tmpptr;
  tmpptr = this->RowBounds;
  this->RowBounds = this->OldRowBounds;
  this->OldRowBounds = tmpptr;
  

  // 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>this->ImageViewportSize[0]-1)?
    (this->ImageViewportSize[0]-1):(maxX);
  maxY = (maxY>this->ImageViewportSize[1]-1)?
    (this->ImageViewportSize[1]-1):(maxY);

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

  // What is a power of 2 size big enough to fit this image?
  this->ImageMemorySize[0] = 32;
  this->ImageMemorySize[1] = 32;
  while ( this->ImageMemorySize[0] < this->ImageInUseSize[0] )
    {
    this->ImageMemorySize[0] *= 2;
    }
  while ( this->ImageMemorySize[1] < this->ImageInUseSize[1] )
    {
    this->ImageMemorySize[1] *= 2;
    }
  
  this->ImageOrigin[0] = static_cast<int>(minX);
  this->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] > 2*this->ImageMemorySize[0] ||
       oldImageMemorySize[1] > 2*this->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] >= this->ImageMemorySize[0] &&
       oldImageMemorySize[1] >= this->ImageMemorySize[1] )
    {
    this->ImageMemorySize[0] = oldImageMemorySize[0];
    this->ImageMemorySize[1] = oldImageMemorySize[1];
    }
  
  // Do we already have a texture big enough? If not, create a new one and
  // clear it.
  if ( !this->Image ||
       this->ImageMemorySize[0] > oldImageMemorySize[0] ||
       this->ImageMemorySize[1] > oldImageMemorySize[1] )
    {
    // If there is an image there must be row bounds
    if ( this->Image )
      {
      delete [] this->Image;
      delete [] this->RowBounds;
      delete [] this->OldRowBounds;
      }
    
    this->Image = new unsigned char[(this->ImageMemorySize[0] *
                                     this->ImageMemorySize[1] * 4)];
    
    // Create the row bounds array. This will store the start / stop pixel
    // for each row. This helps eleminate 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*this->ImageMemorySize[1]];
    this->OldRowBounds = new int [2*this->ImageMemorySize[1]];

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

  // 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 < this->ImageInUseSize[1]; j++ )
      {
      this->RowBounds[j*2] = 0;
      this->RowBounds[j*2+1] = this->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*this->ImageViewportSize[0] - this->ImageOrigin[0];
      
      y1 = (viewPoint[lineIndex[i][0]][1]+1.0) * 
        0.5*this->ImageViewportSize[1] - this->ImageOrigin[1];
      
      x2 = (viewPoint[lineIndex[i][1]][0]+1.0) * 
        0.5*this->ImageViewportSize[0] - this->ImageOrigin[0];
      
      y2 = (viewPoint[lineIndex[i][1]][1]+1.0) * 
        0.5*this->ImageViewportSize[1] - this->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 < this->ImageInUseSize[1]; j++ )
      {
      this->RowBounds[j*2] = this->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>this->ImageInUseSize[0]-1)?
            (this->ImageInUseSize[0]-1):(xlow);

          xhigh = (xhigh<0)?(0):(xhigh);
          xhigh = (xhigh>this->ImageInUseSize[0]-1)?(
            this->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] = this->ImageMemorySize[0];
        this->RowBounds[j*2+1] = -1;
        }
      }
    }
  
  for ( j = this->ImageInUseSize[1]; j < this->ImageMemorySize[1]; j++ )
    {
    this->RowBounds[j*2] = this->ImageMemorySize[0];
    this->RowBounds[j*2+1] = -1;
    }

  for ( j = 0; j < this->ImageMemorySize[1]; j++ )
    {
    // 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 = this->Image + 4*( j*this->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 = this->Image + 4*( j*this->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 = this->Image + 4*( j*this->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 vtkVolumeRayCastMapper::ComputeMatrices( vtkImageData *data, 
                                                 vtkVolume *vol )
{
  // Get the data spacing. This scaling is not accounted for in
  // the volume's matrix, so we must add it in.
  float volumeSpacing[3];
  data->GetSpacing( volumeSpacing );
  
  // Get the origin of the data.  This translation is not accounted for in
  // the volume's matrix, so we must add it in.
  float volumeOrigin[3];
  data->GetOrigin( volumeOrigin );

  // Get the dimensions of the data.
  int volumeDimensions[3];
  data->GetDimensions( volumeDimensions );

  vtkTransform *voxelsTransform = this->VoxelsTransform;
  vtkTransform *voxelsToViewTransform = this->VoxelsToViewTransform;
  
  // 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.
  this->VolumeMatrix->DeepCopy( vol->GetMatrix() );
  voxelsToViewTransform->SetMatrix( VolumeMatrix );

  // Create a transform that will account for the scaling and translation of
  // the scalar data. The is the volume to voxels matrix.
  voxelsTransform->Identity();
  voxelsTransform->Translate(volumeOrigin[0], 
                             volumeOrigin[1], 
                             volumeOrigin[2] );
  
  voxelsTransform->Scale( volumeSpacing[0],
                          volumeSpacing[1],
                          volumeSpacing[2] );
  
  // Now concatenate the volume's matrix with this scalar data matrix
  voxelsToViewTransform->PreMultiply();
  voxelsToViewTransform->Concatenate( voxelsTransform->GetMatrix() );

  // Now we actually have the world to voxels matrix - copy it out
  this->WorldToVoxelsMatrix->DeepCopy( voxelsToViewTransform->GetMatrix() );
  this->WorldToVoxelsMatrix->Invert();
  
  // We also want to invert this to get voxels to world
  this->VoxelsToWorldMatrix->DeepCopy( voxelsToViewTransform->GetMatrix() );
  
  // Compute the voxels to view transform by concatenating the
  // voxels to world matrix with the projection matrix (world to view)
  voxelsToViewTransform->PostMultiply();
  voxelsToViewTransform->Concatenate( this->PerspectiveMatrix );
  
  this->VoxelsToViewMatrix->DeepCopy( voxelsToViewTransform->GetMatrix() );
  
  this->ViewToVoxelsMatrix->DeepCopy( this->VoxelsToViewMatrix );
  this->ViewToVoxelsMatrix->Invert();  
}

void vtkVolumeRayCastMapper::InitializeClippingPlanes( 
                                           vtkVolumeRayCastStaticInfo *staticInfo,
                                           vtkPlaneCollection *planes )
{
  vtkPlane *onePlane;
  float    worldNormal[3], worldOrigin[3];
  float    volumeOrigin[4];
  int      i;
  float    *worldToVoxelsMatrix;
  float    *voxelsToWorldMatrix;
  int      count;
  float    *clippingPlane;
  float    t;
  
  count = planes->GetNumberOfItems();
  staticInfo->NumberOfClippingPlanes = count;

  if ( !count )
    {
    return;
    }
  
  worldToVoxelsMatrix = staticInfo->WorldToVoxelsMatrix;
  voxelsToWorldMatrix = staticInfo->VoxelsToWorldMatrix;

  staticInfo->ClippingPlane = new float [4*count];
  
  // loop through all the clipping planes
  for ( i = 0; i < count; i++ )
    {
    onePlane = (vtkPlane *)planes->GetItemAsObject(i);
    onePlane->GetNormal(worldNormal);
    onePlane->GetOrigin(worldOrigin);
    clippingPlane = staticInfo->ClippingPlane + 4*i;
    vtkVRCMultiplyNormalMacro( worldNormal, 
                               clippingPlane,
                               voxelsToWorldMatrix );
    vtkVRCMultiplyPointMacro( worldOrigin, volumeOrigin,
                              worldToVoxelsMatrix );

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

    clippingPlane[3] = -(clippingPlane[0]*volumeOrigin[0] + 
                         clippingPlane[1]*volumeOrigin[1] + 
                         clippingPlane[2]*volumeOrigin[2]);
    }
}


int vtkVolumeRayCastMapper::ClipRayAgainstClippingPlanes( 
                                           vtkVolumeRayCastDynamicInfo *dynamicInfo, 
                                           vtkVolumeRayCastStaticInfo *staticInfo )
{
  float    *clippingPlane;
  int      i;
  float    rayDir[3];
  float    t, point[3], dp;
  float    *rayStart, *rayEnd;

  rayStart = dynamicInfo->TransformedStart;
  rayEnd = dynamicInfo->TransformedEnd;

  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 < staticInfo->NumberOfClippingPlanes; i++ )
    {
    clippingPlane = staticInfo->ClippingPlane + 4*i;

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

    if ( dp != 0.0 )
      {
      t = 
        -( clippingPlane[0]*rayStart[0] + 
           clippingPlane[1]*rayStart[1] + 
           clippingPlane[2]*rayStart[2] + clippingPlane[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;
          }
        }
      }
    }

  return 1;
}

int vtkVolumeRayCastMapper::ClipRayAgainstVolume( 
                                            vtkVolumeRayCastDynamicInfo *dynamicInfo, 
                                            float bounds[6] )
{
  int    loop;
  float  diff;
  float  t;
  float  *rayStart, *rayEnd, *rayDirection;

  rayStart     = dynamicInfo->TransformedStart;
  rayEnd       = dynamicInfo->TransformedEnd;
  rayDirection = dynamicInfo->TransformedDirection;
  
  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;
    }
    
  return 1;
}

void vtkVolumeRayCastMapper::UpdateShadingTables( vtkRenderer *ren, 
                                                  vtkVolume *vol )
{
  int                   shading;
  vtkVolumeProperty     *volume_property;

  volume_property = vol->GetProperty();

  shading = volume_property->GetShade();

  this->GradientEstimator->SetInput( this->GetInput() );

  if ( shading )
    {
    this->GradientShader->UpdateShadingTable( ren, vol, 
                                              this->GradientEstimator );
    }
}

float vtkVolumeRayCastMapper::GetZeroOpacityThreshold( vtkVolume *vol )
{
  return( this->VolumeRayCastFunction->GetZeroOpacityThreshold( vol ) );
}

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

  os << indent << "Sample Distance: " << this->SampleDistance << "\n";
  os << indent << "Image Sample Distance: " 
     << this->ImageSampleDistance << "\n";
  os << indent << "Minimum Image Sample Distance: " 
     << this->MinimumImageSampleDistance << "\n";
  os << indent << "Maximum Image Sample Distance: " 
     << this->MaximumImageSampleDistance << "\n";
  os << indent << "Auto Adjust Sample Distances: " 
     << this->AutoAdjustSampleDistances << "\n";
  os << indent << "Intermix Intersecting Geometry: "
    << (this->IntermixIntersectingGeometry ? "On\n" : "Off\n");
  
  os << indent << "Number Of Threads: " << this->NumberOfThreads << "\n";

  if ( this->VolumeRayCastFunction )
    {
    os << indent << "Ray Cast Function: " << this->VolumeRayCastFunction<<"\n";
    }
  else
    {
    os << indent << "Ray Cast Function: (none)\n";
    }

  if ( this->GradientEstimator )
    {
      os << indent << "Gradient Estimator: " << (this->GradientEstimator) <<
        endl;
    }
  else
    {
      os << indent << "Gradient Estimator: (none)" << endl;
    }

  if ( this->GradientShader )
    {
      os << indent << "Gradient Shader: " << (this->GradientShader) << endl;
    }
  else
    {
      os << indent << "Gradient Shader: (none)" << endl;
    }

}

