// ========================================================================
// $Id: lightgon.h,v 1.1.1.1 2002/10/21 13:20:03 damez Exp $
//
// lightgon.h
//     Header file for CLightGonio.
//
// Class: CRestrictedGonioLight, CGeneralGonioLight
//
// Licence: the use and distribution of this file is severely limited, please
// see the file 'doc/Licence.txt'. Any non-authorized use can be prosecuted
// under International Law. For further questions, please, e-mail to
// VHavran@seznam.cz or mail to Vlastimil Havran,
// Pohodli 27, 57001 Litomysl, the Czech Republic.
// REPLACEMENT_STRING
//
// Initial coding by Vlastimil Havran, 2002
//

#ifndef __HEMIGON_H__
#define __HEMIGON_H__

// standard headers
#include <stdexcept>
#include <complex>
#include <cmath>
#include <iomanip>

// BRDFLIB headers
#include "configh.h"
#include "functions.h"
//#include "auxfunc.h"

__BEGIN_BRDFLIB_HEADER

// -----------------------------------------------------------------------
// GONIMETRIC SAMPLING
// -----------------------------------------------------------------------
// Generial Notes:
//
// Let our hemisphere is parametrized this way:
//   phi in range <0,2*M_PI> ... horizontal angle
// theta in range <0,M_PI/2> ... vertical angle
// x = r * cos(phi) * sin(theta)
// y = r * sin(phi) * sin(theta)
// z = r * cos(theta)
//
//
// -----------------------------------------------------------------------
// Class CHemisphereGonio
// -----------------------------------------------------------------------
class CHemisphereGonio
{
private:
  //-----------------------------------------------------------------------
  // Local variables, they are initialized during InitIntensities()

  // the number of vertical angles (theta)
  int VERTICAL_COUNT;
  // the number of horizontal angles (phi)
  int HORIZONTAL_COUNT;

  // ----------------------------------------------------------------------
  // description of angles

  // the angles in rads in vertical direction (theta)
  // Note that VERTICAL_ANGLES[0] must be always 0.0 to represent
  // pole of the hemisphere !!
  // Maximum angle can be at most M_PI/2 !!!
  double *VERTICAL_ANGLES;

  // the angles in rads in horizontal direction (phi)
  double *HORIZONTAL_ANGLES;

  // Access to the intensity via index.
  inline int GetIndex(int theta, int phi) {
    assert(theta>=0);
    assert(theta<VERTICAL_COUNT);
    assert(phi>=0);
    assert(phi<HORIZONTAL_COUNT);
    //---------------------------
    return theta + VERTICAL_COUNT * phi;
  }

  // --------------------------------------------------------------------
  // this is the description of the intensities in goniometric diagram
  double *INTENSITY;

  // --------------------------------------------------------------------
  // Auxiliary variables used by the algorithm:

  // sin of vertical angles (theta), after remapping from 3D to 2D plane
  double *SIN_VERTICAL_ANGLES;
  // square of the value above
  double *SIN_VERTICAL_ANGLES_SQR;
  // the multiplier used to stretch the sample in 2D (radius)
  double MULT;
  // index of phi angle that has as a sum maximum energy
  int MULT_index;
  // precomputed array of values required to get fast search
  // in the patches changing the vertical angle
  double *GG;
  // for rescaling in phi angle
  double *ALPHA;
  // multiplier to stretch phi angle
  double _maxPhiRescale;

  // for purposes of visualization
  int last_band, last_sector;

  // -------------------------------------------------------------------
  // Precomputed aggregate values

  // the maximum intensity found in the array
  double maximumIntensity;
  // the total power computed as the sum over all the patches of the
  // goniometric web diagram.
  double totalPower;

  // ------------------------------------------------------------------
  // Allocate arrayes to represent hemisphere intensities
  // and also auxiliary variables.
  void AllocArrays();

  // Initialization of variables required for sampling in theta direction
  void InitTheta();

  // Initialization of variables required for sampling in phi direction
  void InitPhi();

  // Computes total power (flux) over the hemisphere using the intesities
  // description via the goniometric web
  void ComputeTotalPower();

  // Solve a special case of equation: R^3 + R^2 * B + C = 0
  // Returns always real root found (root is always >= 0.0)
  double SolveCubicRoot(bool ed, double B, double C);

  // Searches a digon with remapping of yy variable
  void LocateDigonStretch(double &yy, double &xx,
			  int &left, int &right,
			  double &alpha, double &effectiveDensity);

  // --------------------------------------------------------------------------
  // Searches a digon with remapping of yy variable for inverse setting (going
  // from the hemisphere direction to the 2D rectangle)
  void LocateDigonStretchInverse(
    double &phi, double &radius,
    int &left, int &right,
    double &alpha, double &effectiveDensity);

  // clean all the variables connected with the existence of the
  // hemisphere gonio light
  void CleanUp();
public:
  // constructor and destructor
  CHemisphereGonio();
  ~CHemisphereGonio() { CleanUp();}

  CHemisphereGonio *Dup() const;

  // Initialization of the intensities in keypoints of goniometric web
  // Returns 'false' on success, 'true' on failure. It assumes that the
  // intensities in the diagram are organized in the way:
  //
  //   diagram[theta + theta_num * phi] = intensityValue
  //
  // Further it assumes that theta[0] = 0.0 (hemisphere pole).
  // Further it assumes that phi[0] = 0.0.
  bool InitIntensities(int phi_num, double *phi,
		       int theta_num, double *theta,
		       double *diagram);

  // -----------------------------------------------------------------------
  // Performs that density stretching to fit the description of goniometric
  // light source: modifies 'radius' and 'phi' according to the goniometric
  // light source diagram
  // Returns false, when the sample is valid.
  // Returns true, when the sample cannot be taken as the input of
  // the goniometric light source diagram -- ERROR for current implementation.
  // Variable 'index' is just for debugging purposes, it is the index of
  // the sample to find out the bug
  //
  // Assumes input:
  //   rad to be in the range <0.0-1.0>  ... radius in 2D disk
  //   phi to be in the range <0.0-2*PI> ... phi angle
  //
  // Creates output according to the goniometric diagram description (so like
  // importance sampling):
  //   rad is in the range <0.0-1.0>
  //   phi is in the range <0.0-2*PI>
  bool SampleOne(double &rad, double &phi, int index = 0);

  // The inverse function to the above - the inverse density stretching - for
  // going from 3D space to 2D mapping coordinates
  //
  // Assumes input:
  //   rad to be in the range <0.0-1.0>  ... radius in 2D disk
  //   phi to be in the range <0.0-2*PI> ... phi angle
  //
  // Creates output:
  //   rad is in the range <0.0-1.0>
  //   phi is in the range <0.0-2*PI>
  bool SampleOneInverse(double &rad, double &phi, int index = 0);

  // Given the vector (vecZ >= 0.0), computes the interpolated intensity with
  // respect to the definition of goniometric diagram as the linear
  // interpolation between the four corners of a patch on the goniometric web.
  float GetIntensity(double vecX, double vecY, double vecZ);

  // Returns the maximum radius internally used
  double GetMULT() const { return MULT;}

  // Given the position on the disk by radius 'r' in the range <0.0-1.0>
  // and angle 'phi' in the range <0.0-2*PI>, compute the vector in 3D,
  // with the center of the hemisphere by vector <0.0, 0.0, 1.0> (rad = 0.0).
  inline void
  Map2Dto3D(double &vecX, double &vecY, double &vecZ,
	    double rad, double phi);

  // Given the vector assuming z-coordinate >= 0.0, computes the radius
  // and phi on the 2D disk. The radius computed is then in the range
  // <0.0-1.0> and phi in the range <0.0-2*PI>
  inline void
  Map3Dto2D(double &rad, double &phi,
	    double vecX, double vecY, double vecZ);

  // Returns the total power of the hemisphere integrating over the
  // all small patches of the goniometric web diagram
  double TotalPower() const { return totalPower;}

  int NumBands() const { return VERTICAL_COUNT;}
  int NumSectors() const { return HORIZONTAL_COUNT;}

  // For purposes of the visualization, set the index of the
  // last band and sector
  void GetLastBandAndSector(int &band, int &sector) {
    band = last_band; sector = last_sector;
  }

  void PrintAll()
  {
    cout << "CHemisphereGonio class members after initialization:\n";
    cout << "-------------------------------------------------------------\n";
    if(VERTICAL_COUNT)
      cout << "VERTICAL_COUNT = " << VERTICAL_COUNT << "\n\n";
    if(HORIZONTAL_COUNT)
      cout << "HORIZONTAL_COUNT = " << HORIZONTAL_COUNT << "\n\n";

  if(VERTICAL_ANGLES)
  {
    cout << "VERTICAL_ANGLES:\n";
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
    for(int i = 0; i < VERTICAL_COUNT; i++)
    {
      cout << "VERTICAL_ANGLES[" << i << "] = " << VERTICAL_ANGLES[i] << "\n";
    }
  }
  else cout << "VERTICAL_ANGLES is empty\n";
  cout << "\n";

  if(HORIZONTAL_ANGLES)
  {
    cout << "HORIZONTAL_ANGLES:\n";
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
    for(int i = 0; i <= HORIZONTAL_COUNT; i++)
    {
      cout << "HORIZONTAL_ANGLES[" << i << "] = ";
      cout << HORIZONTAL_ANGLES[i] << "\n";
    }
  }
  else cout << "HORIZONTAL_ANGLES is empty\n";
  cout << "\n";

  if(INTENSITY)
  {
    cout << "INTENSITY:\n";
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
    for(int i = 0; i < HORIZONTAL_COUNT*VERTICAL_COUNT; i++)
    {
      cout << "INTENSITY[" << i << "] = " << INTENSITY[i] << "\t";
      if(!((i+1) % VERTICAL_COUNT)) cout << "\n";
    }
  }
  else cout << "INTENSITY is empty\n";
  cout << "\n";

  if(SIN_VERTICAL_ANGLES)
  {
    cout << "SIN_VERTICAL_ANGLES:\n";
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
    for(int i = 0; i < VERTICAL_COUNT; i++)
    {
      cout << "SIN_VERTICAL_ANGLES[" << i << "] = ";
      cout << SIN_VERTICAL_ANGLES[i] << "\n";
    }
  }
  else cout << "SIN_VERTICAL_ANGLES is empty\n";
  cout << "\n";

  if(SIN_VERTICAL_ANGLES_SQR)
  {
    cout << "SIN_VERTICAL_ANGLES_SQR:\n";
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
    for(int i = 0; i < VERTICAL_COUNT; i++)
    {
      cout << "SIN_VERTICAL_ANGLES_SQR[" << i << "] = ";
      cout << SIN_VERTICAL_ANGLES_SQR[i] << "\n";
    }
  }
  else cout << "SIN_VERTICAL_ANGLES_SQR is empty\n";
  cout << "\n";

  if(MULT) cout << "MULT = " << MULT << "\n\n";
  if(MULT_index) cout << "MULT_index = " << MULT_index << "\n\n";

  if(GG)
  {
    cout << "GG:\n";
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
    for(int i = 0; i < HORIZONTAL_COUNT*VERTICAL_COUNT; i++)
    {
      cout << "GG[" << i << "] = " << GG[i] << "\t";
      if(!((i+1) % VERTICAL_COUNT)) cout << "\n";
    }
  }
  else cout << "GG is empty\n";
  cout << "\n";

  if(ALPHA)
  {
    cout << "ALPHA:\n";
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
    for(int i = 0; i <= HORIZONTAL_COUNT; i++)
    {
      cout << "ALPHA[" << i << "] = " << ALPHA[i] << "\n";
    }
  }
  else cout << "ALPHA is empty\n";
  cout << "\n";

  if(_maxPhiRescale) cout << "_maxPhiRescale = " << _maxPhiRescale << "\n\n";

  //if(last_band > 0 && last_band < ) cout << "last_band = " <<
  // last_band << "\n\n";
  //if(last_sector) cout << "last_sector = " << last_sector << "\n\n";

  if(maximumIntensity)
    cout << "maximumIntensity = " << maximumIntensity << "\n\n";

  if(totalPower)
    cout << "totalPower = " << totalPower << "\n\n";

  }
};

// ----------------------------------------------------------------------
// Inline functions

// Given the position on the disk by radius 'r' in the range <0.0-1.0>
// and angle 'phi' in the range <0.0-2*PI>, compute the vector in 3D,
// with the center of the hemisphere by vector <0.0, 0.0, 1.0> (rad = 0.0).
inline void
CHemisphereGonio::Map2Dto3D(double &vecX, double &vecY, double &vecZ,
			    double rad, double phi)
{
  // Casting rays to a three dimensional space: mapping 2D->3D preserving
  // the density from 2D to 3D
  double o1 = 1.0 - rad * rad;
  double o2 = rad * sqrt(1.0 + o1);
  // the final vector in 3D space is (X,Y,Z)
  vecX = o2 * cos(phi);
  vecY = o2 * sin(phi);
  vecZ = o1;
}

// Given the vector assuming z-coordinate >= 0.0, computes the radius
// and phi on the 2D disk. The radius computed is then in the range
// <0.0-1.0> and phi in the range <0.0-2*PI>
inline void
CHemisphereGonio::Map3Dto2D(double &rad, double &phi,
			    double vecX, double vecY, double vecZ)
{
  double magnitude = sqrt(vecX * vecX + vecY * vecY + vecZ * vecZ);
  assert(magnitude > 0.0);
  magnitude = 1.0 / magnitude;
  vecX *= magnitude;
  vecY *= magnitude;
  vecZ *= magnitude;
  assert(vecZ<=1.0);
  // now rad is in the range <0.0-1.0>
  rad = sqrt(1.0 - vecZ);
  double o2 = rad * sqrt(1.0 + vecZ);
  if (o2 > 1e-20) {
    double cosPHI = vecX / o2;
    double sinPHI = vecY / o2;
    // compute phi by arcusstangent over the two arguments
    phi = atan2(sinPHI, cosPHI);
    if (phi < 0.0)
      phi = 2.0 * M_PI + phi;
  }
  else
    phi = 0.0; // just set since we have no chance to compute it

  // now phi is in the range <0.0-2*PI>
  return;
}

__END_BRDFLIB_HEADER

#endif // __HEMIGON_H__
