// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause
#include "vtkImageResliceMapper.h"

#include "vtkAbstractImageInterpolator.h"
#include "vtkCamera.h"
#include "vtkGarbageCollector.h"
#include "vtkImageChangeInformation.h"
#include "vtkImageData.h"
#include "vtkImageProperty.h"
#include "vtkImageResliceToColors.h"
#include "vtkImageSlice.h"
#include "vtkImageSliceMapper.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkLinearTransform.h"
#include "vtkLookupTable.h"
#include "vtkMath.h"
#include "vtkMatrix3x3.h"
#include "vtkMatrix4x4.h"
#include "vtkObjectFactory.h"
#include "vtkPlane.h"
#include "vtkPoints.h"
#include "vtkRenderer.h"
#include "vtkStreamingDemandDrivenPipeline.h"

// A tolerance to compensate for roundoff errors
#define VTK_RESLICE_MAPPER_VOXEL_TOL 7.62939453125e-06

VTK_ABI_NAMESPACE_BEGIN
vtkStandardNewMacro(vtkImageResliceMapper);

//------------------------------------------------------------------------------
vtkImageResliceMapper::vtkImageResliceMapper()
{
  this->ChangeInformation = vtkImageChangeInformation::New();
  this->SliceMapper = vtkImageSliceMapper::New();
  this->ImageReslice = vtkImageResliceToColors::New();
  this->ResliceMatrix = vtkMatrix4x4::New();
  this->WorldToDataMatrix = vtkMatrix4x4::New();
  this->DataToSliceMatrix = vtkMatrix4x4::New();
  this->SliceToWorldMatrix = vtkMatrix4x4::New();

  this->JumpToNearestSlice = 0;
  this->AutoAdjustImageQuality = 1;
  this->SeparateWindowLevelOperation = 1;
  this->SlabType = VTK_IMAGE_SLAB_MEAN;
  this->SlabThickness = 0.0;
  this->SlabSampleFactor = 2;
  this->ImageSampleFactor = 1;
  this->ResampleToScreenPixels = 1;
  this->InternalResampleToScreenPixels = 0;
  this->ResliceNeedUpdate = 0;

  // streaming requires an output port
  this->SetNumberOfOutputPorts(1);
}

//------------------------------------------------------------------------------
vtkImageResliceMapper::~vtkImageResliceMapper()
{
  if (this->ChangeInformation)
  {
    this->ChangeInformation->Delete();
  }
  if (this->SliceMapper)
  {
    this->SliceMapper->Delete();
  }
  if (this->ImageReslice)
  {
    this->ImageReslice->Delete();
  }
  if (this->ResliceMatrix)
  {
    this->ResliceMatrix->Delete();
  }
  if (this->WorldToDataMatrix)
  {
    this->WorldToDataMatrix->Delete();
  }
  if (this->DataToSliceMatrix)
  {
    this->DataToSliceMatrix->Delete();
  }
  if (this->SliceToWorldMatrix)
  {
    this->SliceToWorldMatrix->Delete();
  }
}

//------------------------------------------------------------------------------
void vtkImageResliceMapper::SetSlicePlane(vtkPlane* plane)
{
  if (this->SlicePlane == plane)
  {
    return;
  }
  if (this->SlicePlane)
  {
    this->SlicePlane->Delete();
  }
  if (!plane)
  {
    this->SlicePlane = vtkPlane::New();
  }
  else
  {
    this->SlicePlane = plane;
    plane->Register(this);
  }

  this->Modified();
}

//------------------------------------------------------------------------------
void vtkImageResliceMapper::SetInterpolator(vtkAbstractImageInterpolator* interpolator)
{
  vtkMTimeType mtime = this->ImageReslice->GetMTime();

  this->ImageReslice->SetInterpolator(interpolator);

  if (this->ImageReslice->GetMTime() > mtime)
  {
    this->Modified();
  }
}

//------------------------------------------------------------------------------
vtkAbstractImageInterpolator* vtkImageResliceMapper::GetInterpolator()
{
  return this->ImageReslice->GetInterpolator();
}

//------------------------------------------------------------------------------
void vtkImageResliceMapper::ReleaseGraphicsResources(vtkWindow* win)
{
  this->SliceMapper->ReleaseGraphicsResources(win);
}

//------------------------------------------------------------------------------
void vtkImageResliceMapper::Render(vtkRenderer* ren, vtkImageSlice* prop)
{
  if (this->ResliceNeedUpdate)
  {
    this->ImageReslice->SetInputConnection(this->GetInputConnection(0, 0));
    this->ImageReslice->UpdateWholeExtent();
    this->ResliceNeedUpdate = 0;

    // Adjust the reslice output to put it into a coordinate system that is
    // aligned with the SlicePlane.  Since this coordinate system depends
    // only on the SlicePlane and nothing else, it is identical for all
    // instances of vtkImageResliceMapper that use the same slice plane,
    // regardless of whether the prop orientation or image orientation
    // differs between mappers.  Hence we expect the mappers to paint their
    // pixels at the same depth in the depth buffer, reducing z-fighting.
    // (Note that vtkImageChangeInformation passes the input data arrays
    // unchanged to the output, so there is no duplication/copying of data).
    this->ChangeInformation->SetInputData(this->ImageReslice->GetOutput());
    double direction[9] = { 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 };
    this->ChangeInformation->SetOutputDirection(direction);
    double origin[4] = { 0.0, 0.0, 0.0, 1.0 };
    this->ImageReslice->GetOutputOrigin(origin);
    this->DataToSliceMatrix->MultiplyPoint(origin, origin);
    this->ChangeInformation->SetOutputOrigin(origin);
    this->ChangeInformation->UpdateWholeExtent();
  }

  // apply checkerboard pattern (should have timestamps)
  vtkImageProperty* property = prop->GetProperty();
  if (property && property->GetCheckerboard() && this->InternalResampleToScreenPixels &&
    !this->SeparateWindowLevelOperation && this->SliceFacesCamera)
  {
    this->CheckerboardImage(this->ImageReslice->GetOutput(), ren->GetActiveCamera(), property);
  }

  // delegate to vtkImageSliceMapper
  this->SliceMapper->SetInputData(this->ChangeInformation->GetOutput());
  // the SliceMapper uses a coordinate system aligned with the SlicePlane
  this->SliceMapper->GetDataToWorldMatrix()->DeepCopy(this->SliceToWorldMatrix);
  // the mapper uses SliceFacesCamera to decide whether to use a polygon
  // for the texture versus using a quad the size of the window
  this->SliceMapper->SetSliceFacesCamera(
    (this->SliceFacesCamera && !this->SeparateWindowLevelOperation));
  this->SliceMapper->SetExactPixelMatch(this->InternalResampleToScreenPixels);
  this->SliceMapper->SetBorder((this->Border || this->InternalResampleToScreenPixels));
  this->SliceMapper->SetBackground((this->Background &&
    !(this->SliceFacesCamera && this->InternalResampleToScreenPixels &&
      !this->SeparateWindowLevelOperation)));
  this->SliceMapper->SetPassColorData(!this->SeparateWindowLevelOperation);
  this->SliceMapper->SetDisplayExtent(this->ImageReslice->GetOutputExtent());

  // render pass info for members of vtkImageStack
  this->SliceMapper->MatteEnable = this->MatteEnable;
  this->SliceMapper->ColorEnable = this->ColorEnable;
  this->SliceMapper->DepthEnable = this->DepthEnable;

  // let vtkImageSliceMapper do the rest of the work
  this->SliceMapper->SetNumberOfThreads(this->NumberOfThreads);
  this->SliceMapper->SetClippingPlanes(this->ClippingPlanes);
  this->SliceMapper->Render(ren, prop);
}

//------------------------------------------------------------------------------
void vtkImageResliceMapper::Update(int port)
{
  // I don't like to override Update, or call Modified() in Update,
  // but this allows updates to be forced where MTimes can't be used
  bool resampleToScreenPixels = (this->ResampleToScreenPixels != 0);
  vtkRenderer* ren = nullptr;

  if (this->AutoAdjustImageQuality && resampleToScreenPixels)
  {
    // only use image-size texture if image is smaller than render window,
    // since otherwise there is far less advantage in doing so
    vtkImageSlice* prop = this->GetCurrentProp();
    ren = this->GetCurrentRenderer();
    if (ren && prop)
    {
      const int* rsize = ren->GetSize();
      int maxrsize = (rsize[0] > rsize[1] ? rsize[0] : rsize[1]);
      int* isize = this->GetInput()->GetDimensions();
      int maxisize = (isize[0] > isize[1] ? isize[0] : isize[1]);
      maxisize = (isize[2] > maxisize ? isize[2] : maxisize);
      if (maxisize <= maxrsize && maxisize <= 1024)
      {
        resampleToScreenPixels = (prop->GetAllocatedRenderTime() >= 1.0);
      }
    }
  }

  if (resampleToScreenPixels)
  {
    // force update if quality has increased to "ResampleToScreenPixels"
    if (!this->InternalResampleToScreenPixels)
    {
      this->Modified();
    }
    else
    {
      // force update if renderer size has changes, since the texture
      // size is equal to the renderer size for "ResampleToScreenPixels"
      if (!ren)
      {
        ren = this->GetCurrentRenderer();
      }
      if (ren)
      {
        int* extent = this->ImageReslice->GetOutputExtent();
        const int* size = ren->GetSize();
        if (size[0] != (extent[1] - extent[0] + 1) || size[1] != (extent[3] - extent[2] + 1))
        {
          this->Modified();
        }
      }
    }
  }
  else if (this->InternalResampleToScreenPixels)
  {
    // if execution reaches this point in the code, then the
    // rendering has just switched to interactive quality, and it is
    // necessary to force update if modified since the last update
    if (this->GetMTime() > this->UpdateTime.GetMTime())
    {
      this->Modified();
    }
    else
    {
      // don't switch yet: wait until the camera changes position,
      // which will cause the MTime to change
      resampleToScreenPixels = true;
    }
  }

  this->InternalResampleToScreenPixels = resampleToScreenPixels;

  // Always update if something else caused the input to update
  vtkImageData* input = this->GetInput();
  if (input && input->GetUpdateTime() > this->UpdateTime.GetMTime())
  {
    this->Modified();
  }

  this->Superclass::Update(port);
  this->UpdateTime.Modified();
}

//------------------------------------------------------------------------------
void vtkImageResliceMapper::Update()
{
  this->Superclass::Update();
}

//------------------------------------------------------------------------------
vtkTypeBool vtkImageResliceMapper::Update(int port, vtkInformationVector*)
{
  // One can't really make requests of a mapper so default to regular
  // update.
  this->Update(port);
  return 1;
}

//------------------------------------------------------------------------------
vtkTypeBool vtkImageResliceMapper::Update(vtkInformation*)
{
  // One can't really make requests of a mapper so default to regular
  // update.
  this->Update();
  return 1;
}

//------------------------------------------------------------------------------
vtkTypeBool vtkImageResliceMapper::ProcessRequest(
  vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  if (request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_INFORMATION()))
  {
    // use superclass method to update some info
    this->Superclass::ProcessRequest(request, inputVector, outputVector);

    // need the prop and renderer
    vtkImageSlice* prop = this->GetCurrentProp();
    vtkRenderer* ren = this->GetCurrentRenderer();

    if (ren && prop)
    {
      vtkImageProperty* property = prop->GetProperty();

      // Get point/normal from camera
      if (this->SliceFacesCamera || this->SliceAtFocalPoint)
      {
        vtkCamera* camera = ren->GetActiveCamera();

        if (this->SliceFacesCamera)
        {
          double normal[3];
          camera->GetDirectionOfProjection(normal);
          normal[0] = -normal[0];
          normal[1] = -normal[1];
          normal[2] = -normal[2];
          this->SlicePlane->SetNormal(normal);
        }

        if (this->SliceAtFocalPoint)
        {
          double point[4];
          camera->GetFocalPoint(point);

          if (this->JumpToNearestSlice)
          {
            double normal[4];
            this->SlicePlane->GetNormal(normal);
            normal[3] = -vtkMath::Dot(point, normal);
            point[3] = 1.0;

            // convert normal from world to image coordinates
            double worldToData[16];
            vtkMatrix4x4* dataToWorld = this->GetDataToWorldMatrix();
            vtkMatrix4x4::Transpose(dataToWorld->GetData(), worldToData);
            vtkMatrix4x4::MultiplyPoint(worldToData, normal, normal);
            double dataToImage[9];
            vtkMatrix3x3::Transpose(this->DataDirection, dataToImage);
            vtkMatrix3x3::MultiplyPoint(dataToImage, normal, normal);

            // find the slice orientation from the normal
            int k = 0;
            double maxsq = 0;
            double sumsq = 0;
            for (int i = 0; i < 3; i++)
            {
              double tmpsq = normal[i] * normal[i];
              sumsq += tmpsq;
              if (tmpsq > maxsq)
              {
                maxsq = tmpsq;
                k = i;
              }
            }

            // make sure slice is aligned with image, i.e. not oblique
            if ((1.0 - maxsq / sumsq) < 1e-12)
            {
              // get the point in image coordinates
              vtkMatrix4x4::Invert(dataToWorld->GetData(), worldToData);
              vtkMatrix4x4::MultiplyPoint(worldToData, point, point);
              vtkMath::Subtract(point, this->DataOrigin, point);
              vtkMatrix3x3::Invert(this->DataDirection, dataToImage);
              vtkMatrix3x3::MultiplyPoint(dataToImage, point, point);

              // set the point to lie exactly on a slice
              double z = point[k] / this->DataSpacing[k];
              if (z > VTK_INT_MIN && z < VTK_INT_MAX)
              {
                int j = vtkMath::Floor(z + 0.5);
                point[k] = j * this->DataSpacing[k];
              }

              // convert back to world coordinates
              vtkMatrix3x3::MultiplyPoint(this->DataDirection, point, point);
              vtkMath::Add(point, this->DataOrigin, point);
              dataToWorld->MultiplyPoint(point, point);
            }
          }

          this->SlicePlane->SetOrigin(point);
        }

      } // end of "Get point/normal from camera"

      // set the matrices
      this->UpdateResliceMatrix(ren, prop);

      // update the coords for the polygon to be textured
      this->UpdatePolygonCoords(ren);

      // set the reslice spacing/origin/extent/axes
      this->UpdateResliceInformation(ren);

      // set the reslice bits related to the property
      this->UpdateResliceInterpolation(property);

      // update anything related to the image coloring
      this->UpdateColorInformation(property);
    }

    // set the number of threads to use when executing
    this->ImageReslice->SetNumberOfThreads(this->NumberOfThreads);

    // delegate request to vtkImageReslice (generally not a good thing to
    // do, but I'm familiar with the vtkImageReslice code that gets called).
    return this->ImageReslice->ProcessRequest(request, inputVector, outputVector);
  }

  if (request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT()))
  {
    if (this->Streaming)
    {
      // delegate request to vtkImageReslice (generally not a good thing to
      // do, but I'm familiar with the vtkImageReslice code that gets called).
      return this->ImageReslice->ProcessRequest(request, inputVector, outputVector);
    }
    else
    {
      vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
      int ext[6];
      inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), ext);
      inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), ext, 6);
    }
    return 1;
  }

  if (request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_DATA()))
  {
    vtkInformation* outInfo = outputVector->GetInformationObject(0);
    vtkImageData* output = vtkImageData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

    // set output extent to avoid re-execution
    output->GetInformation()->Set(vtkDataObject::DATA_EXTENT(),
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT()), 6);

    // do an update of Reslice on the next render
    this->ResliceNeedUpdate = 1;

    return 1;
  }

  return this->Superclass::ProcessRequest(request, inputVector, outputVector);
}

//------------------------------------------------------------------------------
// Update the WorldToData transformation matrix, which is just the
// inverse of the vtkProp3D matrix.
void vtkImageResliceMapper::UpdateWorldToDataMatrix(vtkImageSlice* prop)
{
  // copy the matrix, but only if it has changed (we do this to
  // preserve the modified time of the matrix)
  double tmpmat[16] = { 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 };
  if (!prop->GetIsIdentity())
  {
    vtkMatrix4x4::Invert(*prop->GetMatrix()->Element, tmpmat);
  }
  double* mat = *this->WorldToDataMatrix->Element;
  for (int i = 0; i < 16; i++)
  {
    if (mat[i] != tmpmat[i])
    {
      this->WorldToDataMatrix->DeepCopy(tmpmat);
      break;
    }
  }
}

//------------------------------------------------------------------------------
// Update the SliceToWorld transformation matrix
void vtkImageResliceMapper::UpdateSliceToWorldMatrix(vtkCamera* camera)
{
  // Get slice plane in world coords by passing null as the prop matrix
  double plane[4];
  this->GetSlicePlaneInDataCoords(nullptr, plane);

  // Make sure normal is facing towards camera
  vtkMatrix4x4* viewMatrix = camera->GetViewTransformMatrix();
  double* ndop = viewMatrix->Element[2];
  if (vtkMath::Dot(ndop, plane) < 0)
  {
    plane[0] = -plane[0];
    plane[1] = -plane[1];
    plane[2] = -plane[2];
    plane[3] = -plane[3];
  }

  // The normal is the first three elements
  double* normal = plane;

  // The last element is -dot(normal, origin)
  double dp = -plane[3];

  // Compute rotation angle between camera axis and slice normal
  double vec[3];
  vtkMath::Cross(ndop, normal, vec);
  double costheta = vtkMath::Dot(ndop, normal);
  double sintheta = vtkMath::Norm(vec);
  double theta = atan2(sintheta, costheta);
  if (sintheta != 0)
  {
    vec[0] /= sintheta;
    vec[1] /= sintheta;
    vec[2] /= sintheta;
  }
  // convert to quaternion
  costheta = cos(0.5 * theta);
  sintheta = sin(0.5 * theta);
  double quat[4];
  quat[0] = costheta;
  quat[1] = vec[0] * sintheta;
  quat[2] = vec[1] * sintheta;
  quat[3] = vec[2] * sintheta;
  // convert to matrix
  double mat[3][3];
  vtkMath::QuaternionToMatrix3x3(quat, mat);

  // Create a slice-to-world transform matrix
  // The columns are v1, v2, normal
  vtkMatrix4x4* sliceToWorld = this->SliceToWorldMatrix;

  double v1[3], v2[3];
  vtkMath::Multiply3x3(mat, viewMatrix->Element[0], v1);
  vtkMath::Multiply3x3(mat, viewMatrix->Element[1], v2);

  sliceToWorld->Element[0][0] = v1[0];
  sliceToWorld->Element[1][0] = v1[1];
  sliceToWorld->Element[2][0] = v1[2];
  sliceToWorld->Element[3][0] = 0.0;

  sliceToWorld->Element[0][1] = v2[0];
  sliceToWorld->Element[1][1] = v2[1];
  sliceToWorld->Element[2][1] = v2[2];
  sliceToWorld->Element[3][1] = 0.0;

  sliceToWorld->Element[0][2] = normal[0];
  sliceToWorld->Element[1][2] = normal[1];
  sliceToWorld->Element[2][2] = normal[2];
  sliceToWorld->Element[3][2] = 0.0;

  sliceToWorld->Element[0][3] = -dp * normal[0];
  sliceToWorld->Element[1][3] = -dp * normal[1];
  sliceToWorld->Element[2][3] = dp - dp * normal[2];
  sliceToWorld->Element[3][3] = 1.0;
}

//------------------------------------------------------------------------------
// Update the reslice matrix, which is the slice-to-data matrix.
void vtkImageResliceMapper::UpdateResliceMatrix(vtkRenderer* ren, vtkImageSlice* prop)
{
  // Save the old matrix
  double* matrixElements = *this->ResliceMatrix->Element;
  double oldMatrixElements[16];
  vtkMatrix4x4::DeepCopy(oldMatrixElements, matrixElements);

  // Get world-to-data matrix from the prop matrix
  this->UpdateWorldToDataMatrix(prop);

  // Check if prop matrix is orthonormal
  bool propMatrixIsOrthonormal = false;
  vtkMatrix4x4* propMatrix = nullptr;
  if (!this->InternalResampleToScreenPixels)
  {
    static double tol = 1e-12;
    propMatrix = prop->GetMatrix();
    double* row0 = propMatrix->Element[0];
    double* row1 = propMatrix->Element[1];
    double* row2 = propMatrix->Element[2];
    propMatrixIsOrthonormal =
      (fabs(vtkMath::Dot(row0, row0) - 1.0) < tol && fabs(vtkMath::Dot(row1, row1) - 1.0) < tol &&
        fabs(vtkMath::Dot(row2, row2) - 1.0) < tol && fabs(vtkMath::Dot(row0, row1)) < tol &&
        fabs(vtkMath::Dot(row0, row2)) < tol && fabs(vtkMath::Dot(row1, row2)) < tol);
  }

  // Compute SliceToWorld matrix from camera if prop matrix is not
  // orthonormal or if InternalResampleToScreenPixels is set
  if (this->InternalResampleToScreenPixels || !propMatrixIsOrthonormal)
  {
    this->UpdateSliceToWorldMatrix(ren->GetActiveCamera());
    vtkMatrix4x4::Multiply4x4(
      this->WorldToDataMatrix, this->SliceToWorldMatrix, this->ResliceMatrix);
    vtkMatrix4x4::Invert(this->ResliceMatrix, this->DataToSliceMatrix);
  }
  else
  {
    // Get the matrices used to compute the reslice matrix
    vtkMatrix4x4* resliceMatrix = this->ResliceMatrix;
    vtkMatrix4x4* viewMatrix = ren->GetActiveCamera()->GetViewTransformMatrix();

    // Get slice plane in world coords by passing null as the matrix
    double wplane[4];
    this->GetSlicePlaneInDataCoords(nullptr, wplane);

    // Check whether normal is facing towards camera, the "ndop" is
    // the negative of the direction of projection for the camera
    double* ndop = viewMatrix->Element[2];
    double dotprod = vtkMath::Dot(ndop, wplane);

    // Get slice plane in data coords by passing the prop matrix, flip
    // normal to face the camera
    double plane[4];
    this->GetSlicePlaneInDataCoords(propMatrix, plane);
    if (dotprod < 0)
    {
      plane[0] = -plane[0];
      plane[1] = -plane[1];
      plane[2] = -plane[2];
      plane[3] = -plane[3];

      wplane[0] = -wplane[0];
      wplane[1] = -wplane[1];
      wplane[2] = -wplane[2];
      wplane[3] = -wplane[3];
    }

    // The normal is the first three elements of the plane
    double* normal = plane;

    // Convert normal to a coordinate system aligned with the image grid,
    // using the transposed inverse of the data-to-image transform matrix.
    // (DataDirection is already the inverse of the data-to-image matrix,
    // so we just have to transpose it to get the matrix we need)
    double transposed[9];
    vtkMatrix3x3::Transpose(this->DataDirection, transposed);
    double inormal[3];
    vtkMatrix3x3::MultiplyPoint(transposed, normal, inormal);

    // Find the largest component of the normal
    int maxi = 0;
    double maxv = 0.0;
    for (int i = 0; i < 3; i++)
    {
      double tmp = inormal[i] * inormal[i];
      if (tmp > maxv)
      {
        maxi = i;
        maxv = tmp;
      }
    }

    // Create the corresponding axis
    double axis[3];
    axis[0] = 0.0;
    axis[1] = 0.0;
    axis[2] = 0.0;
    axis[maxi] = ((inormal[maxi] < 0.0) ? -1.0 : 1.0);

    // Create two orthogonal axes
    double saxis[3], taxis[3];
    taxis[0] = 0.0;
    taxis[1] = 1.0;
    taxis[2] = 0.0;
    if (maxi == 1)
    {
      taxis[1] = 0.0;
      taxis[2] = 1.0;
    }
    vtkMath::Cross(taxis, axis, saxis);

    // The last element is -dot(normal, origin)
    double dp = (-plane[3] + wplane[0] * propMatrix->Element[0][3] +
      wplane[1] * propMatrix->Element[1][3] + wplane[2] * propMatrix->Element[2][3]);

    // Compute the rotation angle between the axis and the normal
    double vec[3];
    vtkMath::Cross(axis, inormal, vec);
    double costheta = vtkMath::Dot(axis, inormal);
    double sintheta = vtkMath::Norm(vec);
    double theta = atan2(sintheta, costheta);
    if (sintheta != 0)
    {
      vec[0] /= sintheta;
      vec[1] /= sintheta;
      vec[2] /= sintheta;
    }
    // convert to quaternion
    costheta = cos(0.5 * theta);
    sintheta = sin(0.5 * theta);
    double quat[4];
    quat[0] = costheta;
    quat[1] = vec[0] * sintheta;
    quat[2] = vec[1] * sintheta;
    quat[3] = vec[2] * sintheta;
    // convert to matrix
    double mat[3][3];
    vtkMath::QuaternionToMatrix3x3(quat, mat);
    // move the rotation from image coords to data coords
    vtkMatrix3x3::Multiply3x3(this->DataDirection, *mat, *mat);

    // Create a slice-to-data transform matrix
    // The columns are v1, v2, normal
    double v1[3], v2[3];
    vtkMath::Multiply3x3(mat, saxis, v1);
    vtkMath::Multiply3x3(mat, taxis, v2);

    resliceMatrix->Element[0][0] = v1[0];
    resliceMatrix->Element[1][0] = v1[1];
    resliceMatrix->Element[2][0] = v1[2];
    resliceMatrix->Element[3][0] = 0.0;

    resliceMatrix->Element[0][1] = v2[0];
    resliceMatrix->Element[1][1] = v2[1];
    resliceMatrix->Element[2][1] = v2[2];
    resliceMatrix->Element[3][1] = 0.0;

    resliceMatrix->Element[0][2] = normal[0];
    resliceMatrix->Element[1][2] = normal[1];
    resliceMatrix->Element[2][2] = normal[2];
    resliceMatrix->Element[3][2] = 0.0;

    resliceMatrix->Element[0][3] = dp * (propMatrix->Element[2][0] - normal[0]) -
      (propMatrix->Element[0][3] * propMatrix->Element[0][0] +
        propMatrix->Element[1][3] * propMatrix->Element[1][0] +
        propMatrix->Element[2][3] * propMatrix->Element[2][0]);
    resliceMatrix->Element[1][3] = dp * (propMatrix->Element[2][1] - normal[1]) -
      (propMatrix->Element[0][3] * propMatrix->Element[0][1] +
        propMatrix->Element[1][3] * propMatrix->Element[1][1] +
        propMatrix->Element[2][3] * propMatrix->Element[2][1]);
    resliceMatrix->Element[2][3] = dp * (propMatrix->Element[2][2] - normal[2]) -
      (propMatrix->Element[0][3] * propMatrix->Element[0][2] +
        propMatrix->Element[1][3] * propMatrix->Element[1][2] +
        propMatrix->Element[2][3] * propMatrix->Element[2][2]);
    resliceMatrix->Element[3][3] = 1.0;

    // Compute the DataToSliceMatrix
    vtkMatrix4x4::Invert(resliceMatrix, this->DataToSliceMatrix);

    // Compute the SliceToWorldMatrix
    vtkMatrix4x4::Multiply4x4(propMatrix, resliceMatrix, this->SliceToWorldMatrix);
  }

  // If matrix changed, mark as modified so that Reslice will update
  int matrixChanged = 0;
  for (int j = 0; j < 16; j++)
  {
    matrixChanged |= (matrixElements[j] != oldMatrixElements[j]);
  }
  if (matrixChanged)
  {
    this->ResliceMatrix->Modified();
  }
}

//------------------------------------------------------------------------------
// Do all the fancy math to set up the reslicing
void vtkImageResliceMapper::UpdateResliceInformation(vtkRenderer* ren)
{
  vtkMatrix4x4* resliceMatrix = this->ResliceMatrix;
  vtkImageResliceToColors* reslice = this->ImageReslice;

  int extent[6];
  double spacing[3];
  double origin[3];

  // Get current spacing and origin
  reslice->GetOutputSpacing(spacing);
  reslice->GetOutputOrigin(origin);
  reslice->GetOutputExtent(extent);

  // Get the view matrix
  vtkCamera* camera = ren->GetActiveCamera();
  vtkMatrix4x4* viewMatrix = camera->GetViewTransformMatrix();

  // Get slice plane in world coords by passing null as the matrix
  double plane[4];
  this->GetSlicePlaneInDataCoords(nullptr, plane);

  // Check whether normal is facing towards camera, the "ndop" is
  // the negative of the direction of projection for the camera
  double* ndop = viewMatrix->Element[2];
  if (vtkMath::Dot(ndop, plane) < 0)
  {
    plane[0] = -plane[0];
    plane[1] = -plane[1];
    plane[2] = -plane[2];
    plane[3] = -plane[3];
  }

  // Get the z position of the slice in slice coords
  // (requires plane to be normalized by GetSlicePlaneInDataCoords)
  double z = (plane[2] - 2.0) * plane[3];

  if (this->InternalResampleToScreenPixels)
  {
    // Get the projection matrix
    double aspect = ren->GetTiledAspectRatio();
    vtkMatrix4x4* projMatrix = camera->GetProjectionTransformMatrix(aspect, 0, 1);

    // Compute other useful matrices
    double worldToView[16];
    double viewToWorld[16];
    double planeWorldToView[16];
    vtkMatrix4x4::Multiply4x4(*projMatrix->Element, *viewMatrix->Element, worldToView);
    vtkMatrix4x4::Invert(worldToView, viewToWorld);
    vtkMatrix4x4::Transpose(viewToWorld, planeWorldToView);

    double worldToSlice[16];
    double viewToSlice[16];
    vtkMatrix4x4::Invert(*this->SliceToWorldMatrix->Element, worldToSlice);
    vtkMatrix4x4::Multiply4x4(worldToSlice, viewToWorld, viewToSlice);

    // Transform the plane into view coordinates, using the transpose
    // of the inverse of the world-to-view matrix
    vtkMatrix4x4::MultiplyPoint(planeWorldToView, plane, plane);

    // Compute the bounds in slice coords
    double xmin = VTK_DOUBLE_MAX;
    double xmax = -VTK_DOUBLE_MAX;
    double ymin = VTK_DOUBLE_MAX;
    double ymax = -VTK_DOUBLE_MAX;

    for (int i = 0; i < 4; i++)
    {
      // The four corners of the view
      double x = (((i & 1) == 0) ? -1.0 : 1.0);
      double y = (((i & 2) == 0) ? -1.0 : 1.0);

      double hpoint[4];
      hpoint[0] = x;
      hpoint[1] = y;
      hpoint[2] = 0.0;
      hpoint[3] = 1.0;

      if (fabs(plane[2]) < 1e-6)
      {
        // Looking at plane edge-on, just put some
        // points at front clipping plane, others at back plane
        hpoint[2] = (((i & 1) == 0) ? 0.0 : 1.0);
      }
      else
      {
        // Intersect with the slice plane
        hpoint[2] = -(x * plane[0] + y * plane[1] + plane[3]) / plane[2];

        // Clip to the front and back clipping planes
        if (hpoint[2] < 0)
        {
          hpoint[2] = 0.0;
        }
        else if (hpoint[2] > 1)
        {
          hpoint[2] = 1.0;
        }
      }

      // Transform into slice coords
      vtkMatrix4x4::MultiplyPoint(viewToSlice, hpoint, hpoint);

      x = hpoint[0] / hpoint[3];
      y = hpoint[1] / hpoint[3];

      // Find min/max in slice coords
      if (x < xmin)
      {
        xmin = x;
      }
      if (x > xmax)
      {
        xmax = x;
      }
      if (y < ymin)
      {
        ymin = y;
      }
      if (y > ymax)
      {
        ymax = y;
      }
    }

    // The ResliceExtent is always set to the renderer size,
    // this is the maximum size ever required and sticking to
    // this size avoids any memory reallocation on GPU or CPU
    const int* size = ren->GetSize();
    int xsize = ((size[0] <= 0) ? 1 : size[0]);
    int ysize = ((size[1] <= 0) ? 1 : size[1]);

    extent[0] = 0;
    extent[1] = xsize - 1;
    extent[2] = 0;
    extent[3] = ysize - 1;
    extent[4] = 0;
    extent[5] = 0;

    // Find the spacing
    spacing[0] = (xmax - xmin) / xsize;
    spacing[1] = (ymax - ymin) / ysize;

    // Corner of resliced plane, including half-pixel offset to
    // exactly match texels to pixels in the final rendering
    double newOrigin[4];
    newOrigin[0] = xmin + 0.5 * spacing[0];
    newOrigin[1] = ymin + 0.5 * spacing[1];
    newOrigin[2] = z;
    newOrigin[3] = 1.0;
    resliceMatrix->MultiplyPoint(newOrigin, newOrigin);
    origin[0] = newOrigin[0];
    origin[1] = newOrigin[1];
    origin[2] = newOrigin[2];
  }
  else
  {
    // Compute texel spacing from image spacing
    double inputSpacing[3];
    this->GetInput()->GetSpacing(inputSpacing);
    inputSpacing[0] = fabs(inputSpacing[0]);
    inputSpacing[1] = fabs(inputSpacing[1]);
    inputSpacing[2] = fabs(inputSpacing[2]);
    for (int j = 0; j < 2; j++)
    {
      double xc = resliceMatrix->Element[j][0];
      double yc = resliceMatrix->Element[j][1];
      double zc = resliceMatrix->Element[j][2];
      double s =
        (xc * xc * inputSpacing[0] + yc * yc * inputSpacing[1] + zc * zc * inputSpacing[2]) /
        sqrt(xc * xc + yc * yc + zc * zc);
      s /= this->ImageSampleFactor;
      // only modify if difference is greater than roundoff tolerance
      if (fabs((s - spacing[j]) / s) > 1e-12)
      {
        spacing[j] = s;
      }
    }

    // Find the bounds for the texture
    double xmin = VTK_DOUBLE_MAX;
    double xmax = -VTK_DOUBLE_MAX;
    double ymin = VTK_DOUBLE_MAX;
    double ymax = -VTK_DOUBLE_MAX;

    vtkPoints* points = this->SliceMapper->GetPoints();
    vtkIdType n = points->GetNumberOfPoints();
    if (n == 0)
    {
      // This occurs when the slice doesn't intersect the image.
      // Nothing will be drawn when this occurs, but the pipeline still
      // needs to update.  So we generate one pixel at the Origin
      // (convert the origin from data coords to slice coords).
      double point[4];
      this->GetInput()->GetOrigin(point);
      point[3] = 1.0;
      resliceMatrix->MultiplyPoint(point, point);
      xmin = point[0];
      xmax = point[0];
      ymin = point[1];
      ymax = point[1];
    }

    for (vtkIdType k = 0; k < n; k++)
    {
      double point[3];
      points->GetPoint(k, point);

      xmin = ((xmin < point[0]) ? xmin : point[0]);
      xmax = ((xmax > point[0]) ? xmax : point[0]);
      ymin = ((ymin < point[1]) ? ymin : point[1]);
      ymax = ((ymax > point[1]) ? ymax : point[1]);
    }

    double tol = VTK_RESLICE_MAPPER_VOXEL_TOL;
    int xsize = vtkMath::Floor((xmax - xmin) / spacing[0] + tol);
    int ysize = vtkMath::Floor((ymax - ymin) / spacing[1] + tol);
    if (this->Border == 0)
    {
      xsize += 1;
      ysize += 1;
    }
    if (xsize < 1)
    {
      xsize = 1;
    }
    if (ysize < 1)
    {
      ysize = 1;
    }

    // Keep old size if possible, to avoid memory reallocation
    if ((xsize - 1) > extent[1] || (ysize - 1) > extent[3] || (0.9 * extent[1] / xsize) > 1.0 ||
      (0.9 * extent[3] / ysize) > 1.0)
    {
      extent[1] = xsize - 1;
      extent[3] = ysize - 1;
    }
    extent[0] = 0;
    extent[2] = 0;
    extent[4] = 0;
    extent[5] = 0;

    double x0 = xmin + 0.5 * spacing[0] * (this->Border != 0);
    double y0 = ymin + 0.5 * spacing[1] * (this->Border != 0);

    double corner[4] = { x0, y0, z, 1.0 };
    double newOrigin[4];
    resliceMatrix->MultiplyPoint(corner, newOrigin);

    double dx = newOrigin[0] - origin[0];
    double dy = newOrigin[1] - origin[1];
    double dz = newOrigin[2] - origin[2];

    // only modify origin if it has changed by tolerance
    if (dx * dx + dy * dy + dz * dz > tol * tol * spacing[0] * spacing[1])
    {
      origin[0] = newOrigin[0];
      origin[1] = newOrigin[1];
      origin[2] = newOrigin[2];
    }
  }

  // Get the resliced image direction from the slice to data matrix
  const double* data = resliceMatrix->GetData();
  double direction[9] = {
    data[0], data[1], data[2], // 1st row
    data[4], data[5], data[6], // 2nd row
    data[8], data[9], data[10] // 3rd row
  };

  // Prepare for reslicing
  reslice->SetOutputExtent(extent);
  reslice->SetOutputSpacing(spacing);
  reslice->SetOutputDirection(direction);
  reslice->SetOutputOrigin(origin);

  if ((this->SliceFacesCamera && this->InternalResampleToScreenPixels &&
        !this->SeparateWindowLevelOperation) ||
    this->SlabThickness > 0)
  {
    // if slice follows camera, use reslice to set the border
    reslice->SetBorder(this->Border);
  }
  else
  {
    // tell reslice to use a double-thickness border,
    // since the polygon geometry will dictate the actual size
    reslice->SetBorder(true);
    reslice->SetBorderThickness(1.0);
  }
}

//------------------------------------------------------------------------------
// Do all the fancy math to set up the reslicing
void vtkImageResliceMapper::UpdateColorInformation(vtkImageProperty* property)
{
  vtkScalarsToColors* lookupTable = this->DefaultLookupTable;

  if (property)
  {
    double colorWindow = property->GetColorWindow();
    double colorLevel = property->GetColorLevel();
    if (property->GetLookupTable())
    {
      lookupTable = property->GetLookupTable();
      if (!property->GetUseLookupTableScalarRange())
      {
        lookupTable->SetRange(colorLevel - 0.5 * colorWindow, colorLevel + 0.5 * colorWindow);
      }
    }
    else
    {
      lookupTable->SetRange(colorLevel - 0.5 * colorWindow, colorLevel + 0.5 * colorWindow);
    }
  }
  else
  {
    lookupTable->SetRange(0, 255);
  }
  this->ImageReslice->SetBypass(this->SeparateWindowLevelOperation != 0);
  this->ImageReslice->SetLookupTable(lookupTable);
  double backgroundColor[4] = { 0.0, 0.0, 0.0, 0.0 };
  if (this->Background)
  {
    this->GetBackgroundColor(property, backgroundColor);
    backgroundColor[0] *= 255;
    backgroundColor[1] *= 255;
    backgroundColor[2] *= 255;
    backgroundColor[3] *= 255;
  }
  this->ImageReslice->SetBackgroundColor(backgroundColor);
}

//------------------------------------------------------------------------------
// Set the reslice interpolation and slab thickness
void vtkImageResliceMapper::UpdateResliceInterpolation(vtkImageProperty* property)
{
  // set the interpolation mode and border
  int interpMode = VTK_RESLICE_NEAREST;
  int slabSlices = 1;

  if (property)
  {
    switch (property->GetInterpolationType())
    {
      case VTK_NEAREST_INTERPOLATION:
        interpMode = VTK_RESLICE_NEAREST;
        break;
      case VTK_LINEAR_INTERPOLATION:
        interpMode = VTK_RESLICE_LINEAR;
        break;
      case VTK_CUBIC_INTERPOLATION:
        interpMode = VTK_RESLICE_CUBIC;
        break;
    }
  }

  // set up the slice spacing for slab views
  double spacing[3], inputSpacing[3];
  this->ImageReslice->GetOutputSpacing(spacing);
  this->GetInput()->GetSpacing(inputSpacing);
  inputSpacing[0] = fabs(inputSpacing[0]);
  inputSpacing[1] = fabs(inputSpacing[1]);
  inputSpacing[2] = fabs(inputSpacing[2]);
  double xc = this->ResliceMatrix->Element[2][0];
  double yc = this->ResliceMatrix->Element[2][1];
  double zc = this->ResliceMatrix->Element[2][2];
  spacing[2] = (xc * xc * inputSpacing[0] + yc * yc * inputSpacing[1] + zc * zc * inputSpacing[2]) /
    sqrt(xc * xc + yc * yc + zc * zc);

  // slab slice spacing is half the input slice spacing
  int n = vtkMath::Ceil(this->SlabThickness / spacing[2]);
  slabSlices = 1 + this->SlabSampleFactor * n;
  if (slabSlices > 1)
  {
    spacing[2] = this->SlabThickness / (slabSlices - 1);
  }
  this->ImageReslice->SetOutputSpacing(spacing);
  int slabMode = this->SlabType;
  double scalarScale = 1.0;
  if (slabMode == VTK_IMAGE_SLAB_SUM)
  {
    // "sum" means integrating over the path length of each ray through
    // the volume, so we need to include the sample spacing as a factor
    scalarScale = spacing[2];
  }

  this->ImageReslice->SetInterpolationMode(interpMode);
  this->ImageReslice->SetSlabMode(slabMode);
  this->ImageReslice->SetSlabNumberOfSlices(slabSlices);
  this->ImageReslice->SetScalarScale(scalarScale);
  this->ImageReslice->SlabTrapezoidIntegrationOn();
}

//------------------------------------------------------------------------------
void vtkImageResliceMapper::CheckerboardImage(
  vtkImageData* input, vtkCamera* camera, vtkImageProperty* property)
{
  // Use focal point as center of checkerboard pattern.  This guarantees
  // exactly the same checkerboard for all images in the scene, which is
  // useful when doing multiple overlays.
  double focalPoint[4];
  camera->GetFocalPoint(focalPoint);
  focalPoint[3] = 1.0;

  double worldToSlice[16];
  vtkMatrix4x4::Invert(*this->SliceToWorldMatrix->Element, worldToSlice);

  vtkMatrix4x4::MultiplyPoint(worldToSlice, focalPoint, focalPoint);
  if (focalPoint[3] != 0.0)
  {
    focalPoint[0] /= focalPoint[3];
    focalPoint[1] /= focalPoint[3];
    focalPoint[2] /= focalPoint[3];
  }

  // Get the checkerboard spacing and apply the offset fraction
  double checkSpacing[2], checkOffset[2];
  property->GetCheckerboardSpacing(checkSpacing);
  property->GetCheckerboardOffset(checkOffset);
  checkOffset[0] = checkOffset[0] * checkSpacing[0] + focalPoint[0];
  checkOffset[1] = checkOffset[1] * checkSpacing[1] + focalPoint[1];

  // Adjust according to the origin and spacing of the slice data
  double origin[3], spacing[3];
  input->GetSpacing(spacing);
  input->GetOrigin(origin);
  checkOffset[0] = (checkOffset[0] - origin[0]) / spacing[0];
  checkOffset[1] = (checkOffset[1] - origin[1]) / spacing[1];
  checkSpacing[0] /= spacing[0];
  checkSpacing[1] /= spacing[1];

  // Apply the checkerboard to the data
  int extent[6];
  input->GetExtent(extent);
  unsigned char* data = static_cast<unsigned char*>(input->GetScalarPointerForExtent(extent));

  vtkImageMapper3D::CheckerboardRGBA(data, extent[1] - extent[0] + 1, extent[3] - extent[2] + 1,
    checkOffset[0], checkOffset[1], checkSpacing[0], checkSpacing[1]);
}

//------------------------------------------------------------------------------
// Compute the vertices of the polygon in the slice coordinate system
#define VTK_IRM_MAX_VERTS 32
#define VTK_IRM_MAX_COORDS 96
void vtkImageResliceMapper::UpdatePolygonCoords(vtkRenderer* ren)
{
  // Get the projection matrix
  double aspect = ren->GetTiledAspectRatio();
  vtkCamera* camera = ren->GetActiveCamera();
  vtkMatrix4x4* viewMatrix = camera->GetViewTransformMatrix();
  vtkMatrix4x4* projMatrix = camera->GetProjectionTransformMatrix(aspect, 0, 1);

  // Compute other useful matrices
  double worldToView[16];
  double viewToWorld[16];
  vtkMatrix4x4::Multiply4x4(*projMatrix->Element, *viewMatrix->Element, worldToView);
  vtkMatrix4x4::Invert(worldToView, viewToWorld);

  // Get slice plane in world coords by passing null as the matrix
  double plane[4];
  this->GetSlicePlaneInDataCoords(nullptr, plane);

  // Check whether normal is facing towards camera, the "ndop" is
  // the negative of the direction of projection for the camera
  double* ndop = viewMatrix->Element[2];
  if (vtkMath::Dot(ndop, plane) < 0)
  {
    plane[0] = -plane[0];
    plane[1] = -plane[1];
    plane[2] = -plane[2];
    plane[3] = -plane[3];
  }

  // Get the z position of the slice in slice coords
  // (requires plane to be normalized by GetSlicePlaneInDataCoords)
  double z = (plane[2] - 2.0) * plane[3];

  // Generate a tolerance based on the screen pixel size
  double fpoint[4];
  camera->GetFocalPoint(fpoint);
  fpoint[3] = 1.0;
  vtkMatrix4x4::MultiplyPoint(worldToView, fpoint, fpoint);
  fpoint[0] /= fpoint[3];
  fpoint[1] /= fpoint[3];
  fpoint[2] /= fpoint[3];
  fpoint[3] = 1.0;

  double topOfScreen[4], botOfScreen[4];
  fpoint[1] -= 1.0;
  vtkMatrix4x4::MultiplyPoint(viewToWorld, fpoint, topOfScreen);
  fpoint[1] += 2.0;
  vtkMatrix4x4::MultiplyPoint(viewToWorld, fpoint, botOfScreen);

  topOfScreen[0] /= topOfScreen[3];
  topOfScreen[1] /= topOfScreen[3];
  topOfScreen[2] /= topOfScreen[3];
  topOfScreen[3] = 1.0;

  botOfScreen[0] /= botOfScreen[3];
  botOfScreen[1] /= botOfScreen[3];
  botOfScreen[2] /= botOfScreen[3];
  botOfScreen[3] = 1.0;

  // height of view in world coords at focal point
  double viewHeight = sqrt(vtkMath::Distance2BetweenPoints(topOfScreen, botOfScreen));

  // height of view in pixels
  int height = ren->GetSize()[1];

  double tol = (height == 0 ? 0.5 : viewHeight * 0.5 / height);

  // make the data bounding box (with or without border)
  double b = (this->Border ? 0.5 : VTK_RESLICE_MAPPER_VOXEL_TOL);
  double bounds[6];
  for (int ii = 0; ii < 3; ii++)
  {
    double c = b * this->DataSpacing[ii];
    int lo = this->DataWholeExtent[2 * ii];
    int hi = this->DataWholeExtent[2 * ii + 1];
    if (lo == hi && tol > c)
    { // apply tolerance to avoid degeneracy
      c = tol;
    }
    bounds[2 * ii] = lo * this->DataSpacing[ii] - c;
    bounds[2 * ii + 1] = hi * this->DataSpacing[ii] + c;
  }

  // transform the vertices to the slice coord system
  double xpoints[8], ypoints[8];
  double weights1[8], weights2[8];
  bool above[8], below[8];

  // arrays for the list of polygon points
  int n = 0;
  double newxpoints[VTK_IRM_MAX_VERTS];
  double newypoints[VTK_IRM_MAX_VERTS];
  double cx = 0.0;
  double cy = 0.0;

  for (int i = 0; i < 8; i++)
  {
    double point[4];
    point[0] = bounds[0 + ((i >> 0) & 1)];
    point[1] = bounds[2 + ((i >> 1) & 1)];
    point[2] = bounds[4 + ((i >> 2) & 1)];
    point[3] = 1.0;
    vtkMatrix3x3::MultiplyPoint(this->DataDirection, point, point);
    vtkMath::Add(point, this->DataOrigin, point);
    this->DataToSliceMatrix->MultiplyPoint(point, point);
    xpoints[i] = point[0] / point[3];
    ypoints[i] = point[1] / point[3];
    weights1[i] = point[2] / point[3] - z - 0.5 * this->SlabThickness;
    weights2[i] = weights1[i] + this->SlabThickness;
    below[i] = (weights1[i] < 0);
    above[i] = (weights2[i] >= 0);

    if (this->SlabThickness > 0 && above[i] && below[i])
    {
      newxpoints[n] = xpoints[i];
      newypoints[n] = ypoints[i];
      cx += xpoints[i];
      cy += ypoints[i];
      n++;
    }
  }

  // go through the edges and find the new points
  for (int j = 0; j < 12; j++)
  {
    // verts from edges (sorry about this..)
    int i1 = (j & 3) | (((j << 1) ^ (j << 2)) & 4);
    int i2 = (i1 ^ (1 << (j >> 2)));

    double* weights = weights2;
    bool* side = above;
    int m = 1 + (this->SlabThickness > 0);
    for (int k = 0; k < m; k++)
    {
      if (side[i1] ^ side[i2])
      {
        double w1 = weights[i2];
        double w2 = -weights[i1];
        double x = (w1 * xpoints[i1] + w2 * xpoints[i2]) / (w1 + w2);
        double y = (w1 * ypoints[i1] + w2 * ypoints[i2]) / (w1 + w2);
        newxpoints[n] = x;
        newypoints[n] = y;
        cx += x;
        cy += y;
        n++;
      }
      weights = weights1;
      side = below;
    }
  }

  double coords[VTK_IRM_MAX_COORDS];

  if (n > 0)
  {
    // centroid
    cx /= n;
    cy /= n;

    // sort the points to make a convex polygon
    double angles[VTK_IRM_MAX_VERTS];
    for (int k = 0; k < n; k++)
    {
      double x = newxpoints[k];
      double y = newypoints[k];
      double t = atan2(y - cy, x - cx);
      int kk;
      for (kk = 0; kk < k; kk++)
      {
        if (t < angles[kk])
        {
          break;
        }
      }
      for (int jj = k; jj > kk; --jj)
      {
        int jj3 = jj * 3;
        angles[jj] = angles[jj - 1];
        coords[jj3] = coords[jj3 - 3];
        coords[jj3 + 1] = coords[jj3 - 2];
        coords[jj3 + 2] = coords[jj3 - 1];
      }
      int kk3 = kk * 3;
      angles[kk] = t;
      coords[kk3] = x;
      coords[kk3 + 1] = y;
      coords[kk3 + 2] = z;
    }
  }

  // remove degenerate points
  if (n > 0)
  {
    bool found = true;
    int m = 0;
    do
    {
      m = 0;
      double xl = coords[3 * (n - 1) + 0];
      double yl = coords[3 * (n - 1) + 1];
      for (int k = 0; k < n; k++)
      {
        double x = coords[3 * k + 0];
        double y = coords[3 * k + 1];

        if (((x - xl) * (x - xl) + (y - yl) * (y - yl)) > tol * tol)
        {
          coords[3 * m + 0] = x;
          coords[3 * m + 1] = y;
          xl = x;
          yl = y;
          m++;
        }
      }
      found = (m < n);
      n = m;
    } while (found && n > 0);
  }

  // find convex hull
  if (this->SlabThickness > 0 && n > 0)
  {
    bool found = true;
    int m = 0;
    do
    {
      m = 0;
      double xl = coords[3 * (n - 1) + 0];
      double yl = coords[3 * (n - 1) + 1];
      for (int k = 0; k < n; k++)
      {
        double x = coords[3 * k + 0];
        double y = coords[3 * k + 1];
        int k1 = (k + 1) % n;
        double xn = coords[3 * k1 + 0];
        double yn = coords[3 * k1 + 1];

        if ((xn - xl) * (y - yl) - (yn - yl) * (x - xl) < tol * tol)
        {
          coords[3 * m + 0] = x;
          coords[3 * m + 1] = y;
          xl = x;
          yl = y;
          m++;
        }
      }
      found = (m < n);
      n = m;
    } while (found && n > 0);
  }

  vtkPoints* points = this->SliceMapper->GetPoints();
  if (!points)
  {
    points = vtkPoints::New();
    points->SetDataTypeToDouble();
    this->SliceMapper->SetPoints(points);
    points->Delete();
  }

  points->SetNumberOfPoints(n);
  for (int k = 0; k < n; k++)
  {
    points->SetPoint(k, &coords[3 * k]);
  }
  points->Modified();
}

//------------------------------------------------------------------------------
void vtkImageResliceMapper::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "JumpToNearestSlice: " << (this->JumpToNearestSlice ? "On\n" : "Off\n");
  os << indent << "AutoAdjustImageQuality: " << (this->AutoAdjustImageQuality ? "On\n" : "Off\n");
  os << indent
     << "SeparateWindowLevelOperation: " << (this->SeparateWindowLevelOperation ? "On\n" : "Off\n");
  os << indent << "ResampleToScreenPixels: " << (this->ResampleToScreenPixels ? "On\n" : "Off\n");
  os << indent << "SlabThickness: " << this->SlabThickness << "\n";
  os << indent << "SlabType: " << this->GetSlabTypeAsString() << "\n";
  os << indent << "SlabSampleFactor: " << this->SlabSampleFactor << "\n";
  os << indent << "ImageSampleFactor: " << this->ImageSampleFactor << "\n";
  os << indent << "Interpolator: " << this->GetInterpolator() << "\n";
}

//------------------------------------------------------------------------------
const char* vtkImageResliceMapper::GetSlabTypeAsString()
{
  switch (this->SlabType)
  {
    case VTK_IMAGE_SLAB_MIN:
      return "Min";
    case VTK_IMAGE_SLAB_MAX:
      return "Max";
    case VTK_IMAGE_SLAB_MEAN:
      return "Mean";
    case VTK_IMAGE_SLAB_SUM:
      return "Sum";
  }
  return "";
}

//------------------------------------------------------------------------------
vtkMTimeType vtkImageResliceMapper::GetMTime()
{
  vtkMTimeType mTime = this->Superclass::GetMTime();

  // Check whether interpolator has changed
  vtkAbstractImageInterpolator* interpolator = this->ImageReslice->GetInterpolator();
  if (interpolator)
  {
    vtkMTimeType mTime2 = interpolator->GetMTime();
    if (mTime2 > mTime)
    {
      mTime = mTime2;
    }
  }

  // Include camera in MTime so that REQUEST_INFORMATION
  // will be called if the camera changes
  if (this->SliceFacesCamera || this->SliceAtFocalPoint || this->InternalResampleToScreenPixels)
  {
    vtkRenderer* ren = this->GetCurrentRenderer();
    if (ren)
    {
      vtkCamera* camera = ren->GetActiveCamera();
      vtkMTimeType mTime2 = camera->GetMTime();
      mTime = (mTime2 > mTime ? mTime2 : mTime);
    }
  }

  if (!this->SliceFacesCamera || !this->SliceAtFocalPoint)
  {
    vtkMTimeType sTime = this->SlicePlane->GetMTime();
    mTime = (sTime > mTime ? sTime : mTime);
  }

  vtkImageSlice* prop = this->GetCurrentProp();
  if (prop != nullptr)
  {
    vtkMTimeType mTime2 = prop->GetUserTransformMatrixMTime();
    mTime = (mTime2 > mTime ? mTime2 : mTime);

    vtkImageProperty* property = prop->GetProperty();
    if (property != nullptr)
    {
      bool useMTime = true;
      if (this->SeparateWindowLevelOperation)
      {
        // only care about property if interpolation mode has changed,
        // since interpolation is the only property-related operation
        // done by vtkImageReslice if SeparateWindowLevelOperation is on
        int imode = this->ImageReslice->GetInterpolationMode();
        this->UpdateResliceInterpolation(property);
        useMTime = (imode != this->ImageReslice->GetInterpolationMode());
      }
      if (useMTime)
      {
        mTime2 = property->GetMTime();
        mTime = (mTime2 > mTime ? mTime2 : mTime);

        vtkScalarsToColors* lookupTable = property->GetLookupTable();
        if (lookupTable != nullptr)
        {
          // check the lookup table mtime
          mTime2 = lookupTable->GetMTime();
          mTime = (mTime2 > mTime ? mTime2 : mTime);
        }
      }
    }
  }

  return mTime;
}

//------------------------------------------------------------------------------
void vtkImageResliceMapper::GetIndexBounds(double extent[6])
{
  if (!this->GetInput())
  {
    return;
  }

  this->UpdateInformation();

  extent[0] = this->DataWholeExtent[0];
  extent[1] = this->DataWholeExtent[1];
  extent[2] = this->DataWholeExtent[2];
  extent[3] = this->DataWholeExtent[3];
  extent[4] = this->DataWholeExtent[4];
  extent[5] = this->DataWholeExtent[5];

  // expand by half a pixel if border is on
  double border = 0.5 * (this->Border != 0);

  extent[0] -= border;
  extent[1] += border;
  extent[2] -= border;
  extent[3] += border;
  extent[4] -= border;
  extent[5] += border;
}

//------------------------------------------------------------------------------
double* vtkImageResliceMapper::GetBounds()
{
  if (!this->GetInput())
  {
    vtkMath::UninitializeBounds(this->Bounds);
    return this->Bounds;
  }

  double extent[6];
  this->GetIndexBounds(extent);

  double* spacing = this->DataSpacing;
  double* origin = this->DataOrigin;
  double* direction = this->DataDirection;

  // compute bounds
  for (int k = 0; k < 2; ++k)
  {
    double kval = extent[k + 4];
    for (int j = 0; j < 2; ++j)
    {
      double jval = extent[j + 2];
      for (int i = 0; i < 2; ++i)
      {
        double ival = extent[i];
        double point[3];
        for (int c = 0; c < 3; ++c)
        {
          point[c] = ival * spacing[0] * direction[c * 3] +
            jval * spacing[1] * direction[c * 3 + 1] + kval * spacing[2] * direction[c * 3 + 2] +
            origin[c];
        }
        if (i + j + k == 0)
        {
          this->Bounds[0] = point[0];
          this->Bounds[1] = point[0];
          this->Bounds[2] = point[1];
          this->Bounds[3] = point[1];
          this->Bounds[4] = point[2];
          this->Bounds[5] = point[2];
        }
        else
        {
          for (int c = 0; c < 3; ++c)
          {
            this->Bounds[c * 2] = point[c] < this->Bounds[c * 2] ? point[c] : this->Bounds[c * 2];
            this->Bounds[c * 2 + 1] =
              point[c] > this->Bounds[c * 2 + 1] ? point[c] : this->Bounds[c * 2 + 1];
          }
        }
      }
    }
  }

  return this->Bounds;
}

//------------------------------------------------------------------------------
void vtkImageResliceMapper::ReportReferences(vtkGarbageCollector* collector)
{
  this->Superclass::ReportReferences(collector);
  // These filters share our input and are therefore involved in a
  // reference loop.
  vtkGarbageCollectorReport(collector, this->ImageReslice, "ImageReslice");
  vtkGarbageCollectorReport(collector, this->SliceMapper, "SliceMapper");
}
VTK_ABI_NAMESPACE_END
