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

  Program:   Visualization Toolkit
  Module:    vtkUnstructuredGridLinearRayIntegrator.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.

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

/*
 * Copyright 2004 Sandia Corporation.
 * Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
 * license for use of this work by or on behalf of the
 * U.S. Government. Redistribution and use in source and binary forms, with
 * or without modification, are permitted provided that this Notice and any
 * statement of authorship are reproduced on all copies.
 */

#include "vtkUnstructuredGridLinearRayIntegrator.h"

#include "vtkObjectFactory.h"
#include "vtkVolumeProperty.h"
#include "vtkVolume.h"
#include "vtkDoubleArray.h"
#include "vtkPiecewiseFunction.h"
#include "vtkColorTransferFunction.h"
#include "vtkMath.h"

#include <vector>
#include <set>
#include <algorithm>

#include <math.h>

#ifndef M_SQRTPI
#define M_SQRTPI        1.77245385090551602792981
#endif
#ifndef M_2_SQRTPI
#define M_2_SQRTPI      1.12837916709551257390
#endif
#ifndef M_1_SQRTPI
#define M_1_SQRTPI      (0.5*M_2_SQRTPI)
#endif

//-----------------------------------------------------------------------------

// VTK's native classes for defining transfer functions is actually slow to
// access, so we have to cache it somehow.  This class is straightforward
// copy of the transfer function.
class vtkLinearRayIntegratorTransferFunction
{
public:
  vtkLinearRayIntegratorTransferFunction();
  ~vtkLinearRayIntegratorTransferFunction();

  void GetTransferFunction(vtkColorTransferFunction *color,
                           vtkPiecewiseFunction *opacity,
                           double unit_distance,
                           double scalar_range[2]);
  void GetTransferFunction(vtkPiecewiseFunction *intensity,
                           vtkPiecewiseFunction *opacity,
                           double unit_distance,
                           double scalar_range[2]);

  inline void GetColor(double x, double c[4]);

  struct acolor {
    double c[4];
  };
  double *ControlPoints;
  int NumControlPoints;
  acolor *Colors;

private:
  vtkLinearRayIntegratorTransferFunction(const vtkLinearRayIntegratorTransferFunction&);  // Not implemented.
  void operator=(const vtkLinearRayIntegratorTransferFunction &);  // Not implemented.
};

vtkLinearRayIntegratorTransferFunction::vtkLinearRayIntegratorTransferFunction()
{
  this->ControlPoints = NULL;
  this->Colors = NULL;

  this->NumControlPoints = 0;
}

vtkLinearRayIntegratorTransferFunction::~vtkLinearRayIntegratorTransferFunction()
{
  if (this->ControlPoints) delete[] this->ControlPoints;
  if (this->Colors) delete[] this->Colors;
}

static const double huebends[6] = {
  1.0/6.0, 1.0/3.0, 0.5, 2.0/3.0, 5.0/6.0, 1.0
};

void vtkLinearRayIntegratorTransferFunction::GetTransferFunction(
                                              vtkColorTransferFunction *color,
                                              vtkPiecewiseFunction *opacity,
                                              double unit_distance,
                                              double scalar_range[2])
{
  std::set<double> cpset;

  double *function_range = color->GetRange();
  double *function = color->GetDataPointer();
  while (1)
    {
    cpset.insert(function[0]);
    if (function[0] == function_range[1]) break;
    function += 4;
    }

  if (color->GetColorSpace() != VTK_CTF_RGB)
    {
    // If we are in an HSV color space, we must insert control points
    // in places where the RGB bends.
    double rgb[3], hsv[3];
    double hue1, hue2;
    double x1, x2;
    std::set<double>::iterator i = cpset.begin();
    x1 = *i;
    color->GetColor(x1, rgb);
    vtkMath::RGBToHSV(rgb, hsv);
    hue1 = hsv[0];
    for (i++; i != cpset.end(); i++)
      {
      x2 = *i;
      color->GetColor(x2, rgb);
      vtkMath::RGBToHSV(rgb, hsv);
      hue2 = hsv[0];

      // Are we crossing the 0/1 boundary?
      if (   (color->GetColorSpace() == VTK_CTF_HSV && color->GetHSVWrap() )
          && ((hue1 - hue2 > 0.5) || (hue2 - hue1 > 0.5)) )
        {
        // Yes, we are crossing the boundary.
        if (hue1 > hue2)
          {
          int j;
          for (j = 0; huebends[j] <= hue2; j++)
            {
            double interp = (1-hue1+huebends[j])/(1-hue1+hue2);
            cpset.insert((x2-x1)*interp + x1);
            }
          while (huebends[j] < hue1) j++;
          for ( ; j < 6; j++)
            {
            double interp = (huebends[j]-hue1)/(1-hue1+hue2);
            cpset.insert((x2-x1)*interp + x1);
            }
          }
        else
          {
          int j;
          for (j = 0; huebends[j] <= hue1; j++)
            {
            double interp = (hue1-huebends[j])/(1-hue2+hue1);
            cpset.insert((x2-x1)*interp + x1);
            }
          while (huebends[j] < hue2) j++;
          for ( ; j < 6; j++)
            {
            double interp = (1-huebends[j]+hue1)/(1-hue2+hue1);
            cpset.insert((x2-x1)*interp + x1);
            }
          }
        }
      else
        {
        // No, we are not crossing the boundary.
        int j = 0;
        double minh, maxh;
        if (hue1 < hue2)
          {
          minh = hue1;  maxh = hue2;
          }
        else
          {
          minh = hue2;  maxh = hue1;
          }
        while (huebends[j] < minh) j++;
        for (j = 0; huebends[j] < maxh; j++)
          {
          double interp = (huebends[j]-hue1)/(hue2-hue1);
          cpset.insert((x2-x1)*interp + x1);
          }
        }

      x1 = x2;
      hue1 = hue2;
      }
    }

  function_range = opacity->GetRange();
  function = opacity->GetDataPointer();
  while (1)
    {
    cpset.insert(function[0]);
    if (function[0] == function_range[0]) break;
    function += 2;
    }

  // Add the scalar at the beginning and end of the range so the interpolation
  // is correct there.
  cpset.insert(scalar_range[0]);
  cpset.insert(scalar_range[1]);
  // Make extra sure there are at least two entries in cpset.
  if (cpset.size() < 2)
    {
    cpset.insert(0.0);
    cpset.insert(1.0);
    }

  // Now record control points and colors.
  if (this->ControlPoints) delete[] this->ControlPoints;
  if (this->Colors) delete[] this->Colors;
  this->NumControlPoints = static_cast<int>(cpset.size());
  this->ControlPoints = new double[this->NumControlPoints];
  this->Colors = new acolor[this->NumControlPoints];

  std::copy(cpset.begin(), cpset.end(), this->ControlPoints);
  for (int i = 0; i < this->NumControlPoints; i++)
    {
    color->GetColor(this->ControlPoints[i], this->Colors[i].c);
    this->Colors[i].c[3] = (  opacity->GetValue(this->ControlPoints[i])
                            / unit_distance);
    }
}

void vtkLinearRayIntegratorTransferFunction::GetTransferFunction(
                                              vtkPiecewiseFunction *intensity,
                                              vtkPiecewiseFunction *opacity,
                                              double unit_distance,
                                              double scalar_range[2])
{
  std::set<double> cpset;

  double *function_range = intensity->GetRange();
  double *function = intensity->GetDataPointer();
  while (1)
    {
    cpset.insert(function[0]);
    if (function[0] == function_range[1]) break;
    function += 2;
    }

  function_range = opacity->GetRange();
  function = opacity->GetDataPointer();
  while (1)
    {
    cpset.insert(function[0]);
    if (function[0] == function_range[0]) break;
    function += 2;
    }

  // Add the scalar at the beginning and end of the range so the interpolation
  // is correct there.
  cpset.insert(scalar_range[0]);
  cpset.insert(scalar_range[1]);
  // Make extra sure there are at least two entries in cpset.
  if (cpset.size() < 2)
    {
    cpset.insert(0.0);
    cpset.insert(1.0);
    }

  // Now record control points and colors.
  if (this->ControlPoints) delete[] this->ControlPoints;
  if (this->Colors) delete[] this->Colors;
  this->NumControlPoints = static_cast<int>(cpset.size());
  this->ControlPoints = new double[this->NumControlPoints];
  this->Colors = new acolor[this->NumControlPoints];

  std::copy(cpset.begin(), cpset.end(), this->ControlPoints);
  for (int i = 0; i < this->NumControlPoints; i++)
    {
    // Is setting all the colors to the same value the right thing to do?
    this->Colors[i].c[0] = this->Colors[i].c[1] = this->Colors[i].c[2]
      = intensity->GetValue(this->ControlPoints[i]);
    this->Colors[i].c[3] = (  opacity->GetValue(this->ControlPoints[i])
                            / unit_distance);
    }
}

inline void vtkLinearRayIntegratorTransferFunction::GetColor(double x,
                                                             double c[4])
{
  int i = 1;
  while ((this->ControlPoints[i] < x) && (i < this->NumControlPoints-1)) i++;

  double before = this->ControlPoints[i-1];
  double after = this->ControlPoints[i];

  double interp = (x-before)/(after-before);

  double *beforec = this->Colors[i-1].c;
  double *afterc = this->Colors[i].c;
  c[0] = (1-interp)*beforec[0] + interp*afterc[0];
  c[1] = (1-interp)*beforec[1] + interp*afterc[1];
  c[2] = (1-interp)*beforec[2] + interp*afterc[2];
  c[3] = (1-interp)*beforec[3] + interp*afterc[3];
}

//-----------------------------------------------------------------------------

vtkStandardNewMacro(vtkUnstructuredGridLinearRayIntegrator);

vtkUnstructuredGridLinearRayIntegrator::vtkUnstructuredGridLinearRayIntegrator()
{
  this->Property = NULL;
  this->TransferFunctions = NULL;
  this->NumIndependentComponents = 0;
}

//-----------------------------------------------------------------------------

vtkUnstructuredGridLinearRayIntegrator::~vtkUnstructuredGridLinearRayIntegrator()
{
  delete[] this->TransferFunctions;
}

//-----------------------------------------------------------------------------

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

//-----------------------------------------------------------------------------

void vtkUnstructuredGridLinearRayIntegrator::Initialize(
                                                    vtkVolume *volume,
                                                    vtkDataArray *scalars)
{
  vtkVolumeProperty *property = volume->GetProperty();

  if (   (property == this->Property)
      && (this->TransferFunctionsModified > property->GetMTime()) )
    {
    // Nothing has changed from the last time Initialize was run.
    return;
    }

  int numcomponents = scalars->GetNumberOfComponents();

  this->Property = property;
  this->TransferFunctionsModified.Modified();

  if (!property->GetIndependentComponents())
    {
    // The scalars actually hold material properties.
    if ((numcomponents != 4) && (numcomponents != 2) )
      {
      vtkErrorMacro("Only 2-tuples and 4-tuples allowed for dependent components.");
      }
    return;
    }

  delete[] this->TransferFunctions;

  this->NumIndependentComponents = numcomponents;
  this->TransferFunctions
    = new vtkLinearRayIntegratorTransferFunction[numcomponents];

  for (int component = 0; component < numcomponents; component++)
    {
    if (property->GetColorChannels(component) == 1)
      {
      this->TransferFunctions[component]
        .GetTransferFunction(property->GetGrayTransferFunction(component),
                             property->GetScalarOpacity(component),
                             property->GetScalarOpacityUnitDistance(component),
                             scalars->GetRange(component));
      }
    else
      {
      this->TransferFunctions[component]
        .GetTransferFunction(property->GetRGBTransferFunction(component),
                             property->GetScalarOpacity(component),
                             property->GetScalarOpacityUnitDistance(component),
                             scalars->GetRange(component));
      }
    }
}

//-----------------------------------------------------------------------------

void vtkUnstructuredGridLinearRayIntegrator::Integrate(
                                            vtkDoubleArray *intersectionLengths,
                                            vtkDataArray *nearIntersections,
                                            vtkDataArray *farIntersections,
                                            float color[4])
{
  int numintersections = intersectionLengths->GetNumberOfTuples();
  if (this->Property->GetIndependentComponents())
    {
    int numscalars = nearIntersections->GetNumberOfComponents();
    double *nearScalars = new double[numscalars];
    double *farScalars = new double[numscalars];
    std::set<double> segments;
    for (vtkIdType i = 0; i < numintersections; i++)
      {
      double total_length = intersectionLengths->GetValue(i);
      nearIntersections->GetTuple(i, nearScalars);
      farIntersections->GetTuple(i, farScalars);

      // Split up segment on control points, because it is nonlinear in
      // these regions.
      segments.erase(segments.begin(), segments.end());
      segments.insert(0.0);
      segments.insert(1.0);
      for (int j = 0; j < numscalars; j++)
        {
        double *cp = this->TransferFunctions[j].ControlPoints;
        vtkIdType numcp = this->TransferFunctions[j].NumControlPoints;
        double minscalar, maxscalar;
        if (nearScalars[j] < farScalars[j])
          {
          minscalar = nearScalars[j];  maxscalar = farScalars[j];
          }
        else
          {
          minscalar = farScalars[j];  maxscalar = nearScalars[j];
          }
        for (int k = 0; k < numcp; k++)
          {
          if (cp[k] <= minscalar) continue;
          if (cp[k] >= maxscalar) break;
          // If we are here, we need to break the segment at the given scalar.
          // Find the fraction between the near and far segment points.
          segments.insert(  (cp[k]-nearScalars[j])
                          / (farScalars[j]-nearScalars[j]));
          }
        }

      // Iterate over all the segment pieces (from front to back) and
      // integrate each piece.
      std::set<double>::iterator segi = segments.begin();
      double nearInterpolant = *segi;
      for (segi++; segi != segments.end(); segi++)
        {
        double farInterpolant = *segi;
        double nearcolor[4] = {0.0, 0.0, 0.0, 0.0};
        double farcolor[4] = {0.0, 0.0, 0.0, 0.0};
        double length = total_length*(farInterpolant-nearInterpolant);
        // Here we handle the mixing of material properties.  This never
        // seems to be defined very clearly.  I handle this by assuming
        // that each scalar represents a cloud of particles of a certian
        // color and a certain density.  We mix the scalars in the same way
        // as mixing these particles together.  By necessity, the density
        // becomes greater.  The "opacity" parameter is really interpreted
        // as the attenuation coefficient (which is proportional to
        // density) and can therefore easily be greater than one.  The
        // opacity of the resulting color will, however, always be scaled
        // between 0 and 1.
        for (int j = 0; j < numscalars; j++)
          {
          double scalar
            = (farScalars[j]-nearScalars[j])*nearInterpolant + nearScalars[j];
          if (j == 0)
            {
            this->TransferFunctions[j].GetColor(scalar, nearcolor);
            }
          else
            {
            double c[4];
            this->TransferFunctions[j].GetColor(scalar, c);
            if (c[3] + nearcolor[3] > 1.0e-8f)
              {
              nearcolor[0] *= nearcolor[3]/(c[3] + nearcolor[3]);
              nearcolor[1] *= nearcolor[3]/(c[3] + nearcolor[3]);
              nearcolor[2] *= nearcolor[3]/(c[3] + nearcolor[3]);
              nearcolor[0] += c[0]*c[3]/(c[3] + nearcolor[3]);
              nearcolor[1] += c[1]*c[3]/(c[3] + nearcolor[3]);
              nearcolor[2] += c[2]*c[3]/(c[3] + nearcolor[3]);
              nearcolor[3] += c[3];
              }
            }

          scalar
            = (farScalars[j]-nearScalars[j])*farInterpolant + nearScalars[j];
          if (j == 0)
            {
            this->TransferFunctions[j].GetColor(scalar, farcolor);
            }
          else
            {
            double c[4];
            this->TransferFunctions[j].GetColor(scalar, c);
            if (c[3] + farcolor[3] > 1.0e-8f)
              {
              farcolor[0] *= farcolor[3]/(c[3] + farcolor[3]);
              farcolor[1] *= farcolor[3]/(c[3] + farcolor[3]);
              farcolor[2] *= farcolor[3]/(c[3] + farcolor[3]);
              farcolor[0] += c[0]*c[3]/(c[3] + farcolor[3]);
              farcolor[1] += c[1]*c[3]/(c[3] + farcolor[3]);
              farcolor[2] += c[2]*c[3]/(c[3] + farcolor[3]);
              farcolor[3] += c[3];
              }
            }
          }
        this->IntegrateRay(length, nearcolor, nearcolor[3],
                           farcolor, farcolor[3], color);

        nearInterpolant = farInterpolant;
        }
      }
      delete[] nearScalars;
      delete[] farScalars;
    }
  else
    {
    double unitdistance = this->Property->GetScalarOpacityUnitDistance();
    if (nearIntersections->GetNumberOfComponents() == 4)
      {
      for (vtkIdType i = 0; i < numintersections; i++)
        {
        double length = intersectionLengths->GetValue(i);
        double *nearcolor = nearIntersections->GetTuple(i);
        double *farcolor = farIntersections->GetTuple(i);
        this->IntegrateRay(length, nearcolor, nearcolor[3]/unitdistance,
                           farcolor, farcolor[3]/unitdistance, color);
        }
      }
    else  // Two components.
      {
      for (vtkIdType i = 0; i < numintersections; i++)
        {
        double length = intersectionLengths->GetValue(i);
        double *nearcolor = nearIntersections->GetTuple(i);
        double *farcolor = farIntersections->GetTuple(i);
        this->IntegrateRay(length, nearcolor[0], nearcolor[1]/unitdistance,
                           farcolor[0], farcolor[1]/unitdistance, color);
        }
      }
    }
}

//-----------------------------------------------------------------------------

void vtkUnstructuredGridLinearRayIntegrator::IntegrateRay(
                                                       double length,
                                                       double intensity_front,
                                                       double attenuation_front,
                                                       double intensity_back,
                                                       double attenuation_back,
                                                       float color[4])
{
  float Psi = vtkUnstructuredGridLinearRayIntegrator::Psi(length,
                                                          attenuation_front,
                                                          attenuation_back);
  float zeta = (float)exp(-0.5*length*(attenuation_front+attenuation_back));
  float alpha = 1-zeta;

  float newintensity = (1-color[3])*(  intensity_front*(1-Psi)
                                     + intensity_back*(Psi-zeta) );
  // Is setting the RGB values the same the right thing to do?
  color[0] += newintensity;
  color[1] += newintensity;
  color[2] += newintensity;
  color[3] += (1-color[3])*alpha;
}

void vtkUnstructuredGridLinearRayIntegrator::IntegrateRay(
                                                    double length,
                                                    const double color_front[3],
                                                    double attenuation_front,
                                                    const double color_back[3],
                                                    double attenuation_back,
                                                    float color[4])
{
  float Psi = vtkUnstructuredGridLinearRayIntegrator::Psi(length,
                                                          attenuation_front,
                                                          attenuation_back);
  float zeta = (float)exp(-0.5*length*(attenuation_front+attenuation_back));
  float alpha = 1-zeta;

  color[0] += (1-color[3])*(color_front[0]*(1-Psi) + color_back[0]*(Psi-zeta));
  color[1] += (1-color[3])*(color_front[1]*(1-Psi) + color_back[1]*(Psi-zeta));
  color[2] += (1-color[3])*(color_front[2]*(1-Psi) + color_back[2]*(Psi-zeta));
  color[3] += (1-color[3])*alpha;
}

//-----------------------------------------------------------------------------

static inline float erf_fitting_function(float u)
{
  return
    - 1.26551223 + u*(1.00002368 + u*(0.37409196 + u*(0.09678418 + 
        u*(-0.18628806 + u*(0.27886807 + u*(-1.13520398 + u*(1.48851587 +
        u*(-0.82215223 + u*0.17087277))))))));
}

#if 0
// This function is not used directly.  It is here for reference.
static inline float erf(float x)
{
  /* Compute as described in Numerical Recipes in C++ by Press, et al. */
/*   x = abs(x);        In this application, x should always be >= 0. */
  float u = 1/(1 + 0.5*x);
  float ans = u*exp(-x*x + erf_fitting_function(u));
/*   return (x >= 0 ? 1 - ans : ans - 1);    x should always be >= 0. */
  return 1 - ans;
}
#endif

/* Compute Dawson's integral as described in Numerical Recipes in C++ by
   Press, et al. */
#define H 0.4
#define NMAX 6
static const float dawson_constant0 = 0.852144;
static const float dawson_constant1 = 0.236928;
static const float dawson_constant2 = 0.0183156;
static const float dawson_constant3 = 0.000393669;
static const float dawson_constant4 = 2.35258e-6;
static const float dawson_constant5 = 3.90894e-9;
static inline float dawson(float x)
{
  if (x > 0.2)
    {
/*  x = abs(x);       In this application, x should always be >= 0. */
    int n0 = 2*(int)((0.5/H)*x + 0.5);
    float xp = x - (float)n0*H;
    float e1 = exp((2*H)*xp);
    float e2 = e1*e1;
    float d1 = n0 + 1;
    float d2 = d1 - 2;
    float sum = 0;
    sum = dawson_constant0*(e1/d1 + 1/(d2*e1));
    d1 += 2;  d2 -= 2;  e1 *= e2;
    sum += dawson_constant1*(e1/d1 + 1/(d2*e1));
    d1 += 2;  d2 -= 2;  e1 *= e2;
    sum += dawson_constant2*(e1/d1 + 1/(d2*e1));
    d1 += 2;  d2 -= 2;  e1 *= e2;
    sum += dawson_constant3*(e1/d1 + 1/(d2*e1));
    d1 += 2;  d2 -= 2;  e1 *= e2;
    sum += dawson_constant4*(e1/d1 + 1/(d2*e1));
    d1 += 2;  d2 -= 2;  e1 *= e2;
    sum += dawson_constant5*(e1/d1 + 1/(d2*e1));
    return M_1_SQRTPI*exp(-xp*xp)*sum;
    }
  else
    {
    float x2 = x*x;
    return x*(1 - (2.0/3.0)*x2*(1 - .4*x2*(1 - (2.0/7.0)*x2)));
    }
}

#if 0
// This function is not used directly.  It is here for reference.
inline float erfi(float x)
{
  return M_2_SQRTPI*exp(x*x)*dawson(x);
}
#endif

float vtkUnstructuredGridLinearRayIntegrator::Psi(float length,
                                                  float attenuation_front,
                                                  float attenuation_back)
{
  float difftauD = length*fabs(attenuation_back - attenuation_front);

  if (difftauD < 1.0e-8f)
    {
    // Volume is homogeneous (with respect to attenuation).
    float tauD = length * attenuation_front;
    if (tauD < 1.0e-8f)
      {
      return 1;
      }
    else
      {
      return (1 - (float)exp(-tauD))/tauD;
      }
    }
  else
    {
    float invsqrt2difftauD = 1/(float)sqrt(2*difftauD);
    float frontterm = length*invsqrt2difftauD*attenuation_front;
    float backterm  = length*invsqrt2difftauD*attenuation_back;
    if (attenuation_back > attenuation_front)
      {
      float u, Y;
      u = 1/(1+0.5f*frontterm);
      Y = u*(float)exp(erf_fitting_function(u));
      u = 1/(1+0.5f*backterm);
      Y += -u*exp(  frontterm*frontterm-backterm*backterm
                          + erf_fitting_function(u));
      Y *= M_SQRTPI*invsqrt2difftauD;
      return Y;
      }
    else
      {
      float expterm = (float)exp(backterm*backterm-frontterm*frontterm);
      return 2*invsqrt2difftauD*(dawson(frontterm) - expterm*dawson(backterm));
      }
    }
}