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

  Program:   Visualization Toolkit
  Module:    vtkVolumePicker.cxx

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

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

=========================================================================*/
#include "vtkVolumePicker.h"
#include "vtkObjectFactory.h"

#include "vtkBox.h"
#include "vtkImageData.h"
#include "vtkVolume.h"
#include "vtkVolumeMapper.h"

vtkStandardNewMacro(vtkVolumePicker);

//------------------------------------------------------------------------------
vtkVolumePicker::vtkVolumePicker()
{
  this->PickCroppingPlanes = 0;
  this->CroppingPlaneId = -1;
}

//------------------------------------------------------------------------------
vtkVolumePicker::~vtkVolumePicker() = default;

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

  os << indent << "PickCroppingPlanes: " << (this->PickCroppingPlanes ? "On" : "Off") << "\n";

  os << indent << "CroppingPlaneId: " << this->CroppingPlaneId << "\n";
}

//------------------------------------------------------------------------------
void vtkVolumePicker::ResetPickInfo()
{
  this->Superclass::ResetPickInfo();

  this->CroppingPlaneId = -1;
}

//------------------------------------------------------------------------------
// Intersect a vtkVolume with a line by ray casting.  Compared to the
// same method in the superclass, this method will look for cropping planes.

double vtkVolumePicker::IntersectVolumeWithLine(const double p1[3], const double p2[3], double t1,
  double t2, vtkProp3D* prop, vtkAbstractVolumeMapper* mapper)
{
  double tMin = VTK_DOUBLE_MAX;

  vtkImageData* data = vtkImageData::SafeDownCast(mapper->GetDataSetInput());
  vtkVolumeMapper* vmapper = vtkVolumeMapper::SafeDownCast(mapper);

  if (data == nullptr)
  {
    // This picker only works with image inputs
    return VTK_DOUBLE_MAX;
  }

  // Convert ray to structured coordinates
  double spacing[3], origin[3];
  int extent[6];
  data->GetSpacing(spacing);
  data->GetOrigin(origin);
  data->GetExtent(extent);

  double x1[3], x2[3];
  for (int i = 0; i < 3; i++)
  {
    x1[i] = (p1[i] - origin[i]) / spacing[i];
    x2[i] = (p2[i] - origin[i]) / spacing[i];
  }

  // These are set to the plane that the ray enters through
  int planeId = -1;
  int extentPlaneId = -1;

  // There might be multiple regions, depending on cropping flags
  int numSegments = 1;
  double t1List[16], t2List[16], s1List[16];
  int planeIdList[16];
  t1List[0] = t1;
  t2List[0] = t2;
  // s1 is the cropping plane intersection, initialize to large value
  double s1 = s1List[0] = VTK_DOUBLE_MAX;
  planeIdList[0] = -1;

  // Find the cropping bounds in structured coordinates
  double bounds[6];
  for (int j = 0; j < 6; j++)
  {
    bounds[j] = extent[j];
  }

  if (vmapper && vmapper->GetCropping())
  {
    vmapper->GetCroppingRegionPlanes(bounds);
    for (int j = 0; j < 3; j++)
    {
      double b1 = (bounds[2 * j] - origin[j]) / spacing[j];
      double b2 = (bounds[2 * j + 1] - origin[j]) / spacing[j];
      bounds[2 * j] = (b1 < b2 ? b1 : b2);
      bounds[2 * j + 1] = (b1 < b2 ? b2 : b1);
      if (bounds[2 * j] < extent[2 * j])
      {
        bounds[2 * j] = extent[2 * j];
      }
      if (bounds[2 * j + 1] > extent[2 * j + 1])
      {
        bounds[2 * j + 1] = extent[2 * j + 1];
      }
      if (bounds[2 * j] > bounds[2 * j + 1])
      {
        return VTK_DOUBLE_MAX;
      }
    }

    // Get all of the line segments that intersect the visible blocks
    int flags = vmapper->GetCroppingRegionFlags();
    if (!vtkVolumePicker::ClipLineWithCroppingRegion(bounds, extent, flags, x1, x2, t1, t2,
          extentPlaneId, numSegments, t1List, t2List, s1List, planeIdList))
    {
      return VTK_DOUBLE_MAX;
    }
  }
  else
  {
    // If no cropping, then use volume bounds
    double s2;
    if (!vtkVolumePicker::ClipLineWithExtent(extent, x1, x2, s1, s2, extentPlaneId))
    {
      return VTK_DOUBLE_MAX;
    }
    s1List[0] = s1;
    t1List[0] = ((s1 > t1) ? s1 : t1);
    t2List[0] = ((s2 < t2) ? s2 : t2);
  }

  if (this->PickCroppingPlanes && vmapper && vmapper->GetCropping())
  {
    // Only require information about the first intersection
    s1 = s1List[0];
    if (s1 > t1)
    {
      planeId = planeIdList[0];
    }

    // Set data values at the intersected cropping or clipping plane
    if ((tMin = t1List[0]) < this->GlobalTMin)
    {
      this->ResetPickInfo();
      this->DataSet = data;
      this->Mapper = vmapper;

      double x[3];
      for (int j = 0; j < 3; j++)
      {
        x[j] = x1[j] * (1.0 - tMin) + x2[j] * tMin;
        if (planeId >= 0 && j == planeId / 2)
        {
          x[j] = bounds[planeId];
        }
        else if (planeId < 0 && extentPlaneId >= 0 && j == extentPlaneId / 2)
        {
          x[j] = extent[extentPlaneId];
        }
        this->MapperPosition[j] = x[j] * spacing[j] + origin[j];
      }

      this->SetImageDataPickInfo(x, extent);
    }
  }
  else
  {
    // Go through the segments in order, until a hit occurs
    for (int segment = 0; segment < numSegments; segment++)
    {
      if ((tMin = this->Superclass::IntersectVolumeWithLine(
             p1, p2, t1List[segment], t2List[segment], prop, mapper)) < VTK_DOUBLE_MAX)
      {
        s1 = s1List[segment];
        // Keep the first planeId that was set at the first intersection
        // that occurred after t1
        if (planeId < 0 && s1 > t1)
        {
          planeId = planeIdList[segment];
        }
        break;
      }
    }
  }

  if (tMin < this->GlobalTMin)
  {
    this->CroppingPlaneId = planeId;
    // If t1 is at a cropping or extent plane, use the plane normal
    if (planeId < 0)
    {
      planeId = extentPlaneId;
    }
    if (planeId >= 0 && tMin == s1)
    {
      this->MapperNormal[0] = 0.0;
      this->MapperNormal[1] = 0.0;
      this->MapperNormal[2] = 0.0;
      this->MapperNormal[planeId / 2] = 2.0 * (planeId % 2) - 1.0;
      if (spacing[planeId / 2] < 0)
      {
        this->MapperNormal[planeId / 2] = -this->MapperNormal[planeId / 2];
      }
    }
  }

  return tMin;
}

//------------------------------------------------------------------------------
// This method does several things.  Given the volume CroppingRegionPlanes
// stored in bounds (in structured coords), and the volume extent, it
// casts a ray through the 27 "blocks" that the volume has been divided into.
// Each "block" is turned on or off by a bit in "flags".  The result
// of the ray cast is a collection of line segments: the parametric
// start and end of each segment is stored in t1List and t2List respectively.
// If the segment starts at a cropping plane, the planeIdList will store
// the Id of that plane, otherwise planeIdList will store -1 for that segment.

int vtkVolumePicker::ClipLineWithCroppingRegion(const double bounds[6], const int extent[6],
  int flags, const double x1[3], const double x2[3], double t1, double t2, int& extentPlaneId,
  int& numSegments, double* t1List, double* t2List, double* s1List, int* planeIdList)
{
  extentPlaneId = -1;
  numSegments = 0;
  double s1, s2;

  // Start by clipping the line with the volume extent
  if (!vtkVolumePicker::ClipLineWithExtent(extent, x1, x2, s1, s2, extentPlaneId))
  {
    return 0;
  }

  if (s1 >= t1)
  {
    t1 = s1;
  }
  if (s2 <= t2)
  {
    t2 = s2;
  }

  if (t2 < t1)
  {
    return 0;
  }

  // Compute the coordinates that correspond to t1
  double x[3];
  for (int i = 0; i < 3; i++)
  {
    x[i] = x1[i] * (1.0 - t1) + x2[i] * t1;
    // Watch for out-of-bounds due to numerical roundoff
    if (x[i] < extent[2 * i])
    {
      x[i] = extent[2 * i];
    }
    if (x[i] > extent[2 * i + 1])
    {
      x[i] = extent[2 * i + 1];
    }
  }
  if (t1 == s1 && extentPlaneId >= 0)
  {
    // If right on the boundary, set position exactly
    x[extentPlaneId / 2] = extent[extentPlaneId];
  }

  // Find out which block is hit first, store indices and bounds
  int xi[3];
  double blockBounds[6];
  for (int j = 0; j < 3; j++)
  {
    xi[j] = 0;
    blockBounds[2 * j] = extent[2 * j];
    blockBounds[2 * j + 1] = bounds[2 * j];
    // Be particular about the ray direction
    if (x[j] > bounds[2 * j] || (x[j] == bounds[2 * j] && x1[j] < x2[j]))
    {
      xi[j] = 1;
      blockBounds[2 * j] = bounds[2 * j];
      blockBounds[2 * j + 1] = bounds[2 * j + 1];
    }
    if (x[j] > bounds[2 * j + 1] || (x[j] == bounds[2 * j + 1] && x1[j] < x2[j]))
    {
      xi[j] = 2;
      blockBounds[2 * j] = bounds[2 * j + 1];
      blockBounds[2 * j + 1] = extent[2 * j + 1];
    }
  }

  // Loop through the blocks along the ray path
  int plane1 = -1;
  int plane2 = -1;
  for (;;)
  {
    if (!vtkBox::IntersectWithLine(blockBounds, x1, x2, s1, s2, nullptr, nullptr, plane1, plane2))
    {
      // This should never happen, but if it does, stop here
      break;
    }

    int blockId = xi[0] + xi[1] * 3 + xi[2] * 9;
    if ((flags >> blockId) & 1)
    {
      t1List[numSegments] = (t1 > s1 ? t1 : s1);
      t2List[numSegments] = (t2 < s2 ? t2 : s2);
      s1List[numSegments] = s1;
      planeIdList[numSegments] = -1;
      if (plane1 >= 0)
      {
        // Compute  plane1/2  and  plane1%2
        int k = (plane1 >> 1);
        int l = (plane1 & 1);

        // Need to know if the ray is entering the volume, i.e. whether
        // the adjacent block that the ray is coming from is "off", because
        // we can't define a clip plane unless it is off.
        static const int blockInc[3] = { 1, 3, 9 };
        int noPlane = 1;

        if (xi[k] == 1)
        {
          noPlane = (flags >> (blockId + blockInc[k] * (2 * l - 1)) & 1);
          if (!noPlane)
          {
            planeIdList[numSegments] = plane1;
          }
        }
        else if (xi[k] == 0)
        {
          noPlane = (flags >> (blockId + blockInc[k]) & 1);
          if (!noPlane && l == 1)
          {
            planeIdList[numSegments] = 2 * k;
          }
        }
        else if (xi[k] == 2)
        {
          noPlane = (flags >> (blockId - blockInc[k]) & 1);
          if (!noPlane && l == 0)
          {
            planeIdList[numSegments] = 2 * k + 1;
          }
        }
      }

      // Sanity check: allow no segments with negative length
      if (t1List[numSegments] <= t2List[numSegments])
      {
        if (numSegments > 0 && t1List[numSegments] == t2List[numSegments - 1])
        {
          // Concatenate this segment with the previous one
          t2List[numSegments - 1] = t2List[numSegments];
        }
        else
        {
          // Add this segment as a new segment
          numSegments++;
        }
      }
    }

    // If there is no exit plane, the ray terminated and the search is over
    if (plane2 < 0)
    {
      break;
    }

    // Use the exit plane to choose the next block
    int k = plane2 / 2;
    xi[k] += 2 * (plane2 - 2 * k) - 1;

    if (xi[k] == 0)
    {
      blockBounds[2 * k] = extent[2 * k];
      blockBounds[2 * k + 1] = bounds[2 * k];
    }
    else if (xi[k] == 1)
    {
      blockBounds[2 * k] = bounds[2 * k];
      blockBounds[2 * k + 1] = bounds[2 * k + 1];
    }
    else if (xi[k] == 2)
    {
      blockBounds[2 * k] = bounds[2 * k + 1];
      blockBounds[2 * k + 1] = extent[2 * k + 1];
    }
    else
    {
      // Exit, stage right
      break;
    }
  }

  return numSegments;
}