// ---------------------------------------------------------------------------
// Deterministic sampling of goniometric light sources using
// density stretching scheme in 2 dimensions. Works for full sphere.
// Avoids rejecting any input sample.
// ---------------------------------------------------------------------------
//
// Written by Vlastimil Havran, August, 2002.
// e-mail: havran@mpi-sb.mpg.de
//
// Max-Planck-Institut fur Informatik, Germany
// ---------------------------------------------------------------------------

#ifndef __FULLGON_H__
#define __FULLGON_H__

// BRDFLIB headers
#include "configh.h"
#include "hemigon.h"

__BEGIN_BRDFLIB_HEADER

// -----------------------------------------------------------------------
// GONIMETRIC SAMPLING
// -----------------------------------------------------------------------
// Generial Notes:
//
// Let our sphere is parametrized this way:
//   phi in range <0,2*M_PI> ... horizontal angle
// theta in range <0,M_PI>   ... vertical angle
// x = r * cos(phi) * sin(theta)
// y = r * sin(phi) * sin(theta)
// z = r * cos(theta)
//
//
// -----------------------------------------------------------------------
// Class CSphereGonio
// -----------------------------------------------------------------------
class CSphereGonio
{
private:
  // the multiplier used to stretch the sample in 2D (radius)
  double MULT;

  // the treshold to recognize between top and bottom hemisphere
  double RADIUS_THRESHOLD;
  
  // the total power computed as the sum over all the patches of the
  // goniometric web diagram.
  double totalPower;
  // the total power for top hemisphere
  double totalPowerTop;
  // the total power for the bottom hemisphere
  double totalPowerBottom;

  // the two hemispheres used
  CHemisphereGonio *topHemisphere, *bottomHemisphere;
  
  // hemisphere gonio light
  void CleanUp();
public:
  // constructor and destructor
  CSphereGonio();
  ~CSphereGonio() { CleanUp();}

  enum EType {
    // The continuity is not preserved at all, but the pattern on 2D image
    // is preserved on both hemispheres: (Denoted M_A^1 in the paper)
    EE_TwoSquaresNoContinuity,
    
    // In this case we have mapping where the continuity is preserved on
    // the one edge of the square, for the three edges it is not continuous.
    // (Denoted M_A^2 in the paper.) For this mapping remapping must
    // be set in advance and 'topHemisphereFlag' for forward mapping must
    // be set correctly as well!
    EE_TwoSquaresContinuousEdge,

    // In this case we map 2D square to 3D that we have a single
    // point of discontinuity at the pole of the bottom hemisphere (vecZ = -1.0).
    // So the edges of the square in 2D meet at one point. The mapping does
    // NOT preserve the continuity along the equator in general.
    // (This mapping is denoted M_A^3 in the paper.)
    EE_SingleDiscontinuityPoint,

    // In this case we map 2D square to 3D that we have a single
    // point of discontinuity at the pole of the bottom hemisphere (vecZ = -1.0).
    // So the edges of the square in 2D meet at one point. The mapping preserves
    // continuity of mapping along the equator between top and bottom hemisphere in general.
    // (This mapping is denoted M_A^4 in the paper.)
    EE_FullEquatorContinuity,    
  };
  
  // 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,
		       EType type = EE_SingleDiscontinuityPoint);

  // -----------------------------------------------------------------------
  // 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>
  //   topHemisphereFlag is set to 'true' when the sample should be used
  //   for top hemisphere, otherwise ('false') bottom hemisphere is used.
  bool SampleOne(double &rad, double &phi,
		 bool &topHemisphereFlag, 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
  //   topHemisphereFlag specifies whether we should use top or bottom hemisphere.
  //
  // 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,
			bool topHemisphereFlag, 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);

  // 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).
  // For top hemisphere, vecZ >= 0.0, for bottom hemisphere vecZ < 0.0
  inline void
  Map2Dto3D(double &vecX, double &vecY, double &vecZ,
	    double rad, double phi, bool topHemisphere);

  // Given the vector, 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>. If vecZ >= 0.0, it sets top hemisphere flag
  // and assumes that top hemisphere is used.
  inline void
  Map3Dto2D(double &rad, double &phi, bool &topHemisphere,
	    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;}

protected:
  // the mapping type
  EType mappingType;  
};

// ------------------------------------------------------------
// 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).
// For top hemisphere, vecZ >= 0.0, for bottom hemisphere vecZ < 0.0
inline void
CSphereGonio::Map2Dto3D(double &vecX, double &vecY, double &vecZ,
			double rad, double phi, bool topHemisphere)
{
  // 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;
  if (!topHemisphere)
    vecZ = -o1; // vecZ < 0.0
  return;
}

// Given the vector, 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>. If vecZ >= 0.0, it sets top hemisphere flag
// and assumes that top hemisphere is used.
inline void
CSphereGonio::Map3Dto2D(double &rad, double &phi, bool &topHemisphere,
			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;
  if (vecZ < 0.0) {
    vecZ = -vecZ;
    topHemisphere = false;
  }
  else
    topHemisphere = true;
  
  // 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 // __FULLGON_H__
