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

  Program:   Visualization Toolkit
  Module:    vtkMassProperties.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 "vtkMassProperties.h"

#include "vtkCell.h"
#include "vtkDataObject.h"
#include "vtkIdList.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkSmartPointer.h"

vtkStandardNewMacro(vtkMassProperties);

//----------------------------------------------------------------------------
// TODO: replace with std::cbrt() when C++11 is minimum requirement.
static inline double vtkCubeRoot(double x)
{
  if (x > 0.0)
  {
    return pow(x, 1.0 / 3.0);
  }
  else if (x < 0.0)
  {
    return -pow(-x, 1.0 / 3.0);
  }
  else
  {
    return x;
  }
}

//----------------------------------------------------------------------------
// Constructs with initial 0 values.
vtkMassProperties::vtkMassProperties()
{
  this->SurfaceArea = 0.0;
  this->MinCellArea = 0.0;
  this->MaxCellArea = 0.0;
  this->Volume = 0.0;
  this->VolumeProjected = 0.0;
  this->VolumeX = 0.0;
  this->VolumeY = 0.0;
  this->VolumeZ = 0.0;
  this->Kx = 0.0;
  this->Ky = 0.0;
  this->Kz = 0.0;
  this->NormalizedShapeIndex = 0.0;

  this->SetNumberOfOutputPorts(0);
}

//----------------------------------------------------------------------------
// Destroy any allocated memory.
vtkMassProperties::~vtkMassProperties() = default;

//----------------------------------------------------------------------------
// Description:
// This method measures volume, surface area, and normalized shape index.
// Currently, the input is a ploydata which consists of triangles.
int vtkMassProperties::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* vtkNotUsed(outputVector))
{
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);

  // call ExecuteData
  vtkPolyData* input = vtkPolyData::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkIdType cellId, numCells, numPts, numIds;
  double p[3];

  numCells = input->GetNumberOfCells();
  numPts = input->GetNumberOfPoints();
  if (numCells < 1 || numPts < 1)
  {
    vtkErrorMacro(<< "No data to measure...!");
    return 1;
  }

  vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
  ptIds->Allocate(VTK_CELL_SIZE);

  // Traverse all cells, obtaining node coordinates.
  //
  double vol[3], kxyz[3];
  double xp[3]; // to compute volumeproj
  double munc[3], wxyz, wxy, wxz, wyz;
  double area, surfacearea;
  double volumeproj;
  double mincellarea, maxcellarea;
  double a, b, c, s;
  double x[3], y[3], z[3];
  double i[3], j[3], k[3], u[3], absu[3], length;
  double ii[3], jj[3], kk[3];
  double xavg, yavg, zavg;
  vtkIdType idx;

  // Initialize variables ...
  //
  surfacearea = 0.0;
  volumeproj = 0.0;
  mincellarea = VTK_DOUBLE_MAX;
  maxcellarea = 0.0;
  wxyz = 0;
  wxy = 0.0;
  wxz = 0.0;
  wyz = 0.0;
  for (idx = 0; idx < 3; idx++)
  {
    munc[idx] = 0.0;
    vol[idx] = 0.0;
    kxyz[idx] = 0.0;
  }

  for (cellId = 0; cellId < numCells; cellId++)
  {
    if (input->GetCellType(cellId) != VTK_TRIANGLE)
    {
      vtkWarningMacro(<< "Input data type must be VTK_TRIANGLE not " << input->GetCellType(cellId));
      continue;
    }
    input->GetCellPoints(cellId, ptIds);
    numIds = ptIds->GetNumberOfIds();
    assert(numIds == 3);

    // store current vertex (x,y,z) coordinates ...
    //
    for (idx = 0; idx < numIds; idx++)
    {
      input->GetPoint(ptIds->GetId(idx), p);
      x[idx] = p[0];
      y[idx] = p[1];
      z[idx] = p[2];
    }

    // get i j k vectors ...
    //
    i[0] = (x[1] - x[0]);
    j[0] = (y[1] - y[0]);
    k[0] = (z[1] - z[0]);
    i[1] = (x[2] - x[0]);
    j[1] = (y[2] - y[0]);
    k[1] = (z[2] - z[0]);
    i[2] = (x[2] - x[1]);
    j[2] = (y[2] - y[1]);
    k[2] = (z[2] - z[1]);

    // cross product between two vectors, to determine normal vector
    //
    u[0] = (j[0] * k[1] - k[0] * j[1]);
    u[1] = (k[0] * i[1] - i[0] * k[1]);
    u[2] = (i[0] * j[1] - j[0] * i[1]);

    // normalize normal vector to 1
    //
    length = sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]);
    if (length != 0.0)
    {
      u[0] /= length;
      u[1] /= length;
      u[2] /= length;
    }
    else
    {
      u[0] = u[1] = u[2] = 0.0;
    }

    // determine max unit normal component...
    //
    absu[0] = fabs(u[0]);
    absu[1] = fabs(u[1]);
    absu[2] = fabs(u[2]);

    if ((absu[0] > absu[1]) && (absu[0] > absu[2]))
    {
      munc[0]++;
    }
    else if ((absu[1] > absu[0]) && (absu[1] > absu[2]))
    {
      munc[1]++;
    }
    else if ((absu[2] > absu[0]) && (absu[2] > absu[1]))
    {
      munc[2]++;
    }
    else if ((absu[0] == absu[1]) && (absu[0] == absu[2]))
    {
      wxyz++;
    }
    else if ((absu[0] == absu[1]) && (absu[0] > absu[2]))
    {
      wxy++;
    }
    else if ((absu[0] == absu[2]) && (absu[0] > absu[1]))
    {
      wxz++;
    }
    else if ((absu[1] == absu[2]) && (absu[0] < absu[2]))
    {
      wyz++;
    }
    else
    {
      vtkErrorMacro(<< "Unpredicted situation...!");
      return 1;
    }

    // This is reduced to ...
    //
    ii[0] = i[0] * i[0];
    ii[1] = i[1] * i[1];
    ii[2] = i[2] * i[2];
    jj[0] = j[0] * j[0];
    jj[1] = j[1] * j[1];
    jj[2] = j[2] * j[2];
    kk[0] = k[0] * k[0];
    kk[1] = k[1] * k[1];
    kk[2] = k[2] * k[2];

    // area of a triangle...
    //
    a = sqrt(ii[1] + jj[1] + kk[1]);
    b = sqrt(ii[0] + jj[0] + kk[0]);
    c = sqrt(ii[2] + jj[2] + kk[2]);
    s = 0.5 * (a + b + c);
    area = sqrt(fabs(s * (s - a) * (s - b) * (s - c)));
    surfacearea += area;
    if (area < mincellarea)
    {
      mincellarea = area;
    }
    if (area > maxcellarea)
    {
      maxcellarea = area;
    }

    // volume elements ...
    //
    zavg = (z[0] + z[1] + z[2]) / 3.0;
    yavg = (y[0] + y[1] + y[2]) / 3.0;
    xavg = (x[0] + x[1] + x[2]) / 3.0;

    vol[2] += (area * u[2] * zavg);
    vol[1] += (area * u[1] * yavg);
    vol[0] += (area * u[0] * xavg);

    // V  =  (z1+z2+z3)(x1y2-x2y1+x2y3-x3y2+x3y1-x1y3)/6
    // Volume under triangle is projected area of the triangle times
    // the average of the three z values
    vtkMath::Cross(x, y, xp);
    volumeproj += zavg * (xp[0] + xp[1] + xp[2]) / 2;
  }

  // Surface Area ...
  //
  this->SurfaceArea = surfacearea;
  this->MinCellArea = mincellarea;
  this->MaxCellArea = maxcellarea;

  // Weighting factors in Discrete Divergence theorem for volume calculation.
  //
  kxyz[0] = (munc[0] + (wxyz / 3.0) + ((wxy + wxz) / 2.0)) / numCells;
  kxyz[1] = (munc[1] + (wxyz / 3.0) + ((wxy + wyz) / 2.0)) / numCells;
  kxyz[2] = (munc[2] + (wxyz / 3.0) + ((wxz + wyz) / 2.0)) / numCells;
  this->VolumeX = vol[0];
  this->VolumeY = vol[1];
  this->VolumeZ = vol[2];
  this->Kx = kxyz[0];
  this->Ky = kxyz[1];
  this->Kz = kxyz[2];
  this->Volume = (kxyz[0] * vol[0] + kxyz[1] * vol[1] + kxyz[2] * vol[2]);
  this->Volume = fabs(this->Volume);
  this->VolumeProjected = volumeproj;
  this->NormalizedShapeIndex = (sqrt(surfacearea) / vtkCubeRoot(this->Volume)) / 2.199085233;

  return 1;
}

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

  vtkPolyData* input = vtkPolyData::SafeDownCast(this->GetInput(0));
  if (!input)
  {
    return;
  }
  os << indent << "VolumeX: " << this->GetVolumeX() << "\n";
  os << indent << "VolumeY: " << this->GetVolumeY() << "\n";
  os << indent << "VolumeZ: " << this->GetVolumeZ() << "\n";
  os << indent << "Kx: " << this->GetKx() << "\n";
  os << indent << "Ky: " << this->GetKy() << "\n";
  os << indent << "Kz: " << this->GetKz() << "\n";
  os << indent << "Volume:  " << this->GetVolume() << "\n";
  // os << indent << "Volume Projected:  " << this->GetVolumeProjected  () << "\n";
  // os << indent << "Volume Error:  " <<
  //  fabs(this->GetVolume() - this->GetVolumeProjected())   << "\n";
  os << indent << "Surface Area: " << this->GetSurfaceArea() << "\n";
  os << indent << "Min Cell Area: " << this->GetMinCellArea() << "\n";
  os << indent << "Max Cell Area: " << this->GetMaxCellArea() << "\n";
  os << indent << "Normalized Shape Index: " << this->GetNormalizedShapeIndex() << "\n";
}
