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

  Program:   Visualization Toolkit
  Module:    vtkTessellatedBoxSource.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 "vtkTessellatedBoxSource.h"
#include "vtkCellArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkStreamingDemandDrivenPipeline.h"

#include <cassert>

vtkStandardNewMacro(vtkTessellatedBoxSource);

// ----------------------------------------------------------------------------
vtkTessellatedBoxSource::vtkTessellatedBoxSource()
{
  this->Bounds[0] = -0.5;
  this->Bounds[1] = 0.5;
  this->Bounds[2] = -0.5;
  this->Bounds[3] = 0.5;
  this->Bounds[4] = -0.5;
  this->Bounds[5] = 0.5;

  this->Level = 0;
  this->DuplicateSharedPoints = 0;
  this->Quads = 0;
  this->OutputPointsPrecision = SINGLE_PRECISION;

  this->SetNumberOfInputPorts(0); // this is a source.
}

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

// Duplicate point method.
// Each face of the box is defined by the 3 points: an origin, a point along
// a first axis, a point along a second axis. The axes are oriented
// counterclockwise.
// Point id use voxel numbering, not hexahedron numbering.
// voxelPointId=boundingBoxQuads[faceId][facePointId]
// facePointId=0 : origin
// facePointId=1 : point along the first axis
// facePointId=2 : point along the second axis
static int boundingBoxQuads[6][3] = { { 0, 4, 2 }, // -x face
  { 5, 1, 7 },                                     // +x face
  { 0, 1, 4 },                                     // -y face
  { 6, 7, 2 },                                     // +y face
  { 1, 0, 3 },                                     // -z face
  { 4, 5, 6 } };                                   // +z face

// Minimal number of points method.
// Each edge of the box is defined by two vertices in increasing id order.
// vertexId=edges[edge][lowestId=0, highestId=1]
static int edges[12][2] = { { 0, 1 }, // 0
  { 0, 2 },                           // 1
  { 0, 4 },                           // 2
  { 1, 3 },                           // 3
  { 1, 5 },                           // 4
  { 2, 3 },                           // 5
  { 2, 6 },                           // 6
  { 3, 7 },                           // 7
  { 4, 5 },                           // 8
  { 4, 6 },                           // 9
  { 5, 7 },                           // 10
  { 6, 7 } };                         // 11

// Minimal number of points method.
// signed (edge id+1)=faces[face][edge]
// +1 because we cannot encode -0....
static int faces[6][4] = { { 3, 10, -7, -2 }, // 0: -x face
  { -5, 4, 8, -11 },                          // 1: +x face
  { 1, 5, -9, -3 },                           // 2: -y face
  { 12, -8, -6, 7 },                          // 3: +y face
  { -1, 2, 6, -4 },                           // 4: -z face
  { 9, 11, -12, -10 } };                      // 5: +z face

// ----------------------------------------------------------------------------
// Description:
// Called by the superclass. Actual creation of the points and cells
// happens here.
int vtkTessellatedBoxSource::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
{
  // get the info object
  vtkInformation* outInfo = outputVector->GetInformationObject(0);

  // get the output
  vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

  // get the bounds
  double* bounds = this->Bounds;

  vtkPoints* points = vtkPoints::New();

  // Set the desired precision for the points in the output.
  if (this->OutputPointsPrecision == vtkAlgorithm::DOUBLE_PRECISION)
  {
    points->SetDataType(VTK_DOUBLE);
  }
  else
  {
    points->SetDataType(VTK_FLOAT);
  }

  output->SetPoints(points);
  points->Delete();

  // Always create a new vtkCellArray, otherwise it uses the this->Dummy of
  // vtkPolyData...
  vtkCellArray* polys = vtkCellArray::New();
  output->SetPolys(polys);
  polys->Delete();

  if (this->DuplicateSharedPoints)
  {
    this->DuplicateSharedPointsMethod(bounds, points, polys);
  }
  else
  {
    this->MinimalPointsMethod(bounds, points, polys);
  }
  return 1;
}

void vtkTessellatedBoxSource::DuplicateSharedPointsMethod(
  double* bounds, vtkPoints* points, vtkCellArray* polys)
{
  int numberOfPoints = (this->Level + 2) * (this->Level + 2) * 6;
  int changed = points->GetNumberOfPoints() != numberOfPoints;

  if (changed)
  {
    // Topology changed.
    points->SetNumberOfPoints(numberOfPoints);
    polys->Initialize();
  }

  // Iterate over the 6 faces.
  double facePoints[3][3]; // 3 points of 3 coordinates

  int face = 0;
  vtkIdType firstPointId = 0;
  while (face < 6)
  {
    int i = 0;
    while (i < 3)
    {
      int pointId = boundingBoxQuads[face][i];
      facePoints[i][0] = bounds[pointId & 1];
      facePoints[i][1] = bounds[2 + ((pointId >> 1) & 1)];
      facePoints[i][2] = bounds[4 + ((pointId >> 2) & 1)];
      ++i;
    }
    this->BuildFace(points, polys, firstPointId, facePoints, changed);
    firstPointId += (this->Level + 2) * (this->Level + 2);
    ++face;
  }
}

// ----------------------------------------------------------------------------
void vtkTessellatedBoxSource::MinimalPointsMethod(
  double* bounds, vtkPoints* points, vtkCellArray* polys)
{
  int numberOfInternalPointsPerEdge = this->Level;
  int numberOfInternalPointsPerFace = this->Level * this->Level;
  int numberOfPoints = 8 + 12 * numberOfInternalPointsPerEdge + 6 * numberOfInternalPointsPerFace;

  int changed = points->GetNumberOfPoints() != numberOfPoints;
  if (changed)
  {
    // Topology changed.
    points->SetNumberOfPoints(numberOfPoints);
    polys->Initialize();
  }

  // Compute point coordinates.

  // First the 8 vertices
  // Then the 12 * (level) internal edge points (static list)
  // Then the 6 * (level*level) internal face points (static list)

  // The 8 vertices, voxel numbering: id=4*k+2*j+i
  int i = 0;
  double p[3];
  while (i < 8)
  {
    p[0] = bounds[i & 1];
    p[1] = bounds[2 + ((i >> 1) & 1)];
    p[2] = bounds[4 + ((i >> 2) & 1)];
    points->SetPoint(i, p);
    ++i;
  }

  // The 12*(level) internal edge points
  // edges are describe by there lowest vertex Id and highest vertex Id.
  // Numbering start from edges with vertex 0 for lowest Id, in increasing
  // order.
  int currentPointId = i;
  double edgeDirection[3];
  if (this->Level > 0)
  {
    int e = 0;
    while (e < 12)
    {
      assert("check: valid_currentPointId" && (currentPointId == (8 + e * this->Level)));
      double firstPoint[3];
      double lastPoint[3];

      points->GetPoint(edges[e][0], firstPoint);
      points->GetPoint(edges[e][1], lastPoint);
      i = 0;
      while (i < 3)
      {
        edgeDirection[i] = (lastPoint[i] - firstPoint[i]) / (this->Level + 1);
        p[i] = firstPoint[i];
        ++i;
      }

      int j = 1;
      while (j <= this->Level)
      {
        i = 0;
        while (i < 3)
        {
          p[i] += edgeDirection[i];
          ++i;
        }
        points->SetPoint(currentPointId, p);
        ++currentPointId;
        ++j;
      }
      ++e;
    }

    assert("check: valid_currentPointId" && (currentPointId == (8 + 12 * this->Level)));

    int f = 0;
    while (f < 6)
    {
      assert("check: valid_currentPointId" &&
        (currentPointId == (8 + 12 * this->Level + f * this->Level * this->Level)));
      double faceDirections[2][3];
      int facePointId[3];
      e = faces[f][0];
      if (e < 0)
      {
        e = -e;
        facePointId[0] = edges[e - 1][1];
        facePointId[1] = edges[e - 1][0];
      }
      else
      {
        assert("check: not_null_edge" && e > 0);
        facePointId[0] = edges[e - 1][0];
        facePointId[1] = edges[e - 1][1];
      }
      e = faces[f][3];
      if (e < 0)
      {
        e = -e;
        facePointId[2] = edges[e - 1][1];
      }
      else
      {
        facePointId[2] = edges[e - 1][0];
      }

      double facePoints[3][3];

      i = 0;
      while (i < 3)
      {
        int pointId = facePointId[i];
        facePoints[i][0] = bounds[pointId & 1];
        facePoints[i][1] = bounds[2 + ((pointId >> 1) & 1)];
        facePoints[i][2] = bounds[4 + ((pointId >> 2) & 1)];
        ++i;
      }

      int j = 0;
      while (j < 2)
      {
        i = 0;
        while (i < 3)
        {
          faceDirections[j][i] = (facePoints[j + 1][i] - facePoints[0][i]) / (this->Level + 1);
          ++i;
        }
        ++j;
      }

      // internal face points
      j = 1;
      while (j < (this->Level + 1))
      {
        i = 1;
        while (i < (this->Level + 1))
        {
          int comp = 0;
          while (comp < 3)
          {
            p[comp] =
              facePoints[0][comp] + i * faceDirections[0][comp] + j * faceDirections[1][comp];
            ++comp;
          }
          points->SetPoint(currentPointId, p);
          ++currentPointId;
          ++i;
        }
        ++j;
      }

      ++f;
    }
  }
  assert("check: valid_currentPointId" &&
    (currentPointId == (8 + 12 * this->Level + 6 * this->Level * this->Level)));

  // Faces
  vtkIdType poly[4]; // fit a triangle or a quad depending on the Quads flag.
  int f = 0;
  while (f < 6)
  {
    // Create the cells. Two triangle per subquad.
    int j = 0;
    while (j < (this->Level + 1))
    {
      i = 0;
      while (i < (this->Level + 1))
      {
        if (this->Quads)
        {
          // (i,j)
          poly[0] = this->LocalFacePointCoordinatesToPointId(f, i, j);
          // (i+1,j)
          poly[1] = this->LocalFacePointCoordinatesToPointId(f, i + 1, j);
          // (i+1,j+1)
          poly[2] = this->LocalFacePointCoordinatesToPointId(f, i + 1, j + 1);
          // (i,j+1)
          poly[3] = this->LocalFacePointCoordinatesToPointId(f, i, j + 1);
          polys->InsertNextCell(4, poly);
        }
        else
        {
          // First triangle
          // (i,j)
          poly[0] = this->LocalFacePointCoordinatesToPointId(f, i, j);
          // (i+1,j)
          poly[1] = this->LocalFacePointCoordinatesToPointId(f, i + 1, j);
          // (i+1,j+1)
          poly[2] = this->LocalFacePointCoordinatesToPointId(f, i + 1, j + 1);
          polys->InsertNextCell(3, poly);

          // Second triangle
          // poly[0]=(i,j), same than for the first triangle
          // (i+1,j+1)
          poly[1] = poly[2];
          // (i,j+1)
          poly[2] = this->LocalFacePointCoordinatesToPointId(f, i, j + 1);
          polys->InsertNextCell(3, poly);
        }
        ++i;
      }
      ++j;
    }
    ++f;
  }
}

// ----------------------------------------------------------------------------
// Description:
// Compute the pointId of point (i,j) of face f.
// Used by MinimalPointsMethod().
// \pre valid_face: f>=0 && f<6
// \pre valid_i: i>=0 && i<=(this->Level+1)
// \pre valid_j: j>=0 && j<=(this->Level+1)
vtkIdType vtkTessellatedBoxSource::LocalFacePointCoordinatesToPointId(int f, int i, int j)
{
  assert("pre: valid_face" && f >= 0 && f < 6);
  assert("pre: valid_i" && i >= 0 && i <= (this->Level + 1));
  assert("pre: valid_j" && j >= 0 && j <= (this->Level + 1));

  vtkIdType result;

  int e;

  // vertex point?
  if (i == 0 && j == 0)
  {
    e = faces[f][0];
    if (e < 0)
    {
      e = -e;
      result = edges[e - 1][1];
    }
    else
    {
      result = edges[e - 1][0];
    }
  }
  else
  {
    if (i == this->Level + 1 && j == 0)
    {
      e = faces[f][0];
      if (e < 0)
      {
        e = -e;
        result = edges[e - 1][0];
      }
      else
      {
        result = edges[e - 1][1];
      }
    }
    else
    {
      if (i == this->Level + 1 && j == this->Level + 1)
      {
        e = faces[f][1];
        if (e < 0)
        {
          e = -e;
          result = edges[e - 1][0];
        }
        else
        {
          result = edges[e - 1][1];
        }
      }
      else
      {
        if (i == 0 && j == this->Level + 1)
        {
          e = faces[f][3];
          if (e < 0)
          {
            e = -e;
            result = edges[e - 1][1];
          }
          else
          {
            result = edges[e - 1][0];
          }
        }
        else
        {
          // internal edge point?
          if (i == 0)
          {
            e = faces[f][3];
            if (e < 0)
            {
              e = -e;
              result = 8 + (e - 1) * this->Level + j - 1;
            }
            else
            {
              result = 8 + (e - 1) * this->Level + this->Level - j;
            }
          }
          else
          {
            if (i == this->Level + 1)
            {
              e = faces[f][1];
              if (e < 0)
              {
                e = -e;
                result = 8 + (e - 1) * this->Level + this->Level - j;
              }
              else
              {
                result = 8 + (e - 1) * this->Level + j - 1;
              }
            }
            else
            {
              if (j == 0)
              {
                e = faces[f][0];
                if (e < 0)
                {
                  e = -e;
                  result = 8 + (e - 1) * this->Level + this->Level - i;
                }
                else
                {
                  result = 8 + (e - 1) * this->Level + i - 1;
                }
              }
              else
              {
                if (j == this->Level + 1)
                {
                  e = faces[f][2];
                  if (e < 0)
                  {
                    e = -e;
                    result = 8 + (e - 1) * this->Level + i - 1;
                  }
                  else
                  {
                    result = 8 + (e - 1) * this->Level + this->Level - i;
                  }
                }
                else
                {
                  // internal face point.
                  result =
                    (8 + (12 + f * this->Level) * this->Level) + (j - 1) * this->Level + (i - 1);
                }
              }
            }
          }
        }
      }
    }
  }
  assert("post: valid_result" && result >= 0 &&
    result < (8 + 12 * this->Level + 6 * this->Level * this->Level));
  return result;
}

// ----------------------------------------------------------------------------
// Description:
// Build one of the face of the box with some level of tessellation.
// facePoints[0] is the lower-left point
// facePoints[1] is the point along the first axis
// facePoints[2] is the point along the second axis
// \pre positive_id: firstPointId>=0
// \pre points_exists: points!=0
// \pre polys_exists: polys!=0
void vtkTessellatedBoxSource::BuildFace(vtkPoints* points, vtkCellArray* polys,
  vtkIdType firstPointId, double facePoints[3][3], int changed)
{
  assert("pre: positive_id" && firstPointId >= 0);
  assert("pre: points_exists" && points != nullptr);
  assert("pre: polys_exists" && polys != nullptr);

  double direction[2][3];

  int j = 0;
  while (j < 2)
  {
    int i = 0;
    while (i < 3)
    {
      direction[j][i] = (facePoints[j + 1][i] - facePoints[0][i]) / (this->Level + 1);
      ++i;
    }
    ++j;
  }

  // Create the point positions.

  double p[3];

  j = 0;
  while (j <= (this->Level + 1))
  {
    int i = 0;
    while (i <= (this->Level + 1))
    {
      int comp = 0;
      while (comp < 3)
      {
        p[comp] = facePoints[0][comp] + i * direction[0][comp] + j * direction[1][comp];
        ++comp;
      }
      points->SetPoint(firstPointId + j * (this->Level + 2) + i, p);

      ++i;
    }
    ++j;
  }

  if (changed)
  {
    vtkIdType poly[4]; // fit a triangle or a quad depending on the Quads flag.

    // Create the cells. Two triangle per subquad.
    j = 0;
    while (j < (this->Level + 1))
    {
      int i = 0;
      while (i < (this->Level + 1))
      {
        if (this->Quads)
        {
          // (i,j)
          poly[0] = firstPointId + j * (this->Level + 2) + i;
          // (i+1,j)
          poly[1] = firstPointId + j * (this->Level + 2) + i + 1;
          // (i+1,j+1)
          poly[2] = firstPointId + (j + 1) * (this->Level + 2) + i + 1;
          // (i,j+1)
          poly[3] = firstPointId + (j + 1) * (this->Level + 2) + i;
          polys->InsertNextCell(4, poly);
        }
        else
        {
          // First triangle
          // (i,j)
          poly[0] = firstPointId + j * (this->Level + 2) + i;
          // (i+1,j)
          poly[1] = firstPointId + j * (this->Level + 2) + i + 1;
          // (i+1,j+1)
          poly[2] = firstPointId + (j + 1) * (this->Level + 2) + i + 1;
          polys->InsertNextCell(3, poly);

          // Second triangle
          // (i,j)
          poly[0] = firstPointId + j * (this->Level + 2) + i;
          // (i+1,j+1)
          poly[1] = firstPointId + (j + 1) * (this->Level + 2) + i + 1;
          // (i,j+1)
          poly[2] = firstPointId + (j + 1) * (this->Level + 2) + i;
          polys->InsertNextCell(3, poly);
        }
        ++i;
      }
      ++j;
    }
  }
}

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

  os << indent << "Bounds: (" << this->Bounds[0];
  int i = 1;
  while (i < 6)
  {
    os << ", " << this->Bounds[i];
    ++i;
  }
  os << ")" << endl;
  os << indent << "Level: " << this->Level << endl;

  os << indent << "DuplicateSharedPoints: ";
  if (this->DuplicateSharedPoints)
  {
    os << "true";
  }
  else
  {
    os << "false";
  }

  os << indent << "Quads: ";
  if (this->Quads)
  {
    os << "true";
  }
  else
  {
    os << "false";
  }
  os << endl;
  os << indent << "Output Points Precision: " << this->OutputPointsPrecision << endl;
}
