// ---------------------------------------------------------------------------
// 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
// ---------------------------------------------------------------------------

// standard headers
#include <iostream>
#include <cstdio>
using namespace std;

// headers
#include "fullgon.h"


CSphereGonio::CSphereGonio():
  MULT(0.0), RADIUS_THRESHOLD(0.0),
  totalPower(0.0), totalPowerTop(0.0),
  totalPowerBottom(0.0),
  topHemisphere(0),
  bottomHemisphere(0)
{

}

// hemisphere gonio light
void
CSphereGonio::CleanUp()
{
  delete topHemisphere;
  topHemisphere = 0;
  delete bottomHemisphere;
  bottomHemisphere = 0;

  MULT = 0.0;
  RADIUS_THRESHOLD = 0.0;
  totalPower = 0.0;
  totalPowerTop = 0.0;
  totalPowerBottom = 0.0;

  return;
}

// 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
CSphereGonio::InitIntensities(int phi_num, double *phiAngles,
			      int theta_num, double *thetaAngles,
			      double *diagram,
			      EType type)
{
  double minThetaMPI2 = M_PI;
  int    minThetaMPI2index = 0;
  double maxThetaMPI2 = M_PI;
  int    maxThetaMPI2index = theta_num-1;  

  // one of the two mappings;
  mappingType = type;
  
  const float epsilon = 1e-6;  
  // find out the theta angle for which we have M_PI/2
  for (int i = 0; i < theta_num; i++) {
    // diff
    double diff = thetaAngles[i] - M_PI/2;
    // minimum angle lower than M_PI/2
    if (diff <= epsilon) {
      if (fabs(diff) < fabs(minThetaMPI2)) {
	minThetaMPI2 = diff;
	minThetaMPI2index = i;
      }
    }
    // minimum angle greater than M_PI/2
    if (diff >= -epsilon) {
      if (fabs(diff) < fabs(maxThetaMPI2)) {
	maxThetaMPI2 = diff;
	maxThetaMPI2index = i;
      }
    }	
  } // for i

  if (minThetaMPI2index == maxThetaMPI2index) {
    // we have to compute the intensities for top and bottom
    // hemisphere separately
    int thetaTop = minThetaMPI2index+1;
    double *intenAux = new double [theta_num * phi_num];
    double *thetaAux = new double [theta_num];
    for (int theta = 0; theta < thetaTop; theta++) {
      thetaAux[theta] = thetaAngles[theta];
      for (int phi = 0; phi < phi_num; phi++) {
	intenAux[theta + thetaTop * phi] = diagram[theta + theta_num * phi];
      } // for phi
#ifdef _DEBUG      
      if (theta == 0) {
	if (thetaAux[0] != 0.0) {
	  cout << "Top hemisphere - angle is not correctly specified" << endl;
	  exit(3);
	}
	double value = intenAux[0];
	const float epsilon = 1e-5;
	for (int i = 0; i < phi_num; i++) {
	  if (fabs(intenAux[0 + thetaTop * i] - value) > epsilon) {
	    cout << "Top hemisphere - the density/intensity for the pole must be the same!"
		 << endl;
	    exit(3);
	  }
	} // for i
      }
#endif // _DEBUG
    } // for theta

    // allocate and initialize the top hemisphere
    topHemisphere = new CHemisphereGonio();
    bool ok = topHemisphere->InitIntensities(phi_num, phiAngles,
					     thetaTop, thetaAux,
					     intenAux);
    if (ok) {
      cout << "Initialization for the top hemisphere failed" << endl;
      exit(3);
    }

    // ---------------------------------------------------------------------
    // now the intensities for the bottom hemisphere
    int thetaBottom = theta_num - maxThetaMPI2index;
    for (int theta = 0; theta < thetaTop; theta++) {
      thetaAux[theta] = M_PI - thetaAngles[theta_num - theta -1];
      for (int phi = 0; phi < phi_num; phi++) {
	intenAux[theta + thetaBottom * phi] =
	  diagram[theta_num - theta - 1 + theta_num * phi];
      } // for phi
#ifdef _DEBUG      
      if (theta == 0) {
	if (thetaAux[0] != 0.0) {
	  cout << "Bottom hemisphere - angle is not correctly specified" << endl;
	  exit(3);
	}
	double value = intenAux[0];
	const float epsilon = 1e-5;
	for (int i = 0; i < phi_num; i++) {
	  if (fabs(intenAux[0 + thetaBottom * i] - value) > epsilon) {
	    cout << "Bottom hemisphere - the density/intensity for the pole must be the same!"
		 << endl;
	    exit(3);
	  }
	} // for i
      }
#endif // _DEBUG
    } // for theta    

    // allocate and initialize the top hemisphere
    bottomHemisphere = new CHemisphereGonio();
    ok = bottomHemisphere->InitIntensities(phi_num, phiAngles,
					   thetaBottom, thetaAux,
					   intenAux);
    if (ok) {
      cout << "Initialization for the bottom hemisphere failed" << endl;
      exit(3);
    }

    // now compute the total power
    totalPowerTop = topHemisphere->TotalPower();
    totalPowerBottom = bottomHemisphere->TotalPower();
    totalPower = totalPowerTop + totalPowerBottom;

    assert(totalPowerTop > 0.0);
    MULT = sqrt(totalPowerBottom/totalPowerTop + 1.0);
    
    return false; // ok
  }

  // we need to compute the interpolated intensities
  // for equation (theta = M_PI/2);
  cout << "NOT YET IMPLEMENTED" << endl;
  exit(3);
  
  return false;
}


// 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
CSphereGonio::GetIntensity(double vecX, double vecY, double vecZ)
{
  float intensity;
  
  if (vecZ >= 0.0) {
    // top hemisphere
    intensity = topHemisphere->GetIntensity(vecX, vecY, vecZ);
  }
  else {
    // bottom hemisphere
    intensity = bottomHemisphere->GetIntensity(vecX, vecY, -vecZ);
  }

  return intensity;
}

// -----------------------------------------------------------------------
// 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>
//   topHemisphere is set when the sample should be used for top hemisphere,
//    otherwise bottom hemisphere is used.
bool
CSphereGonio::SampleOne(double &rad, double &phi,
			bool &topHemisphereFlag, int index)
{
  assert( (rad >= 0.0) && (rad <= 1.0) );
  assert( (phi >= 0.0) && (phi <= 2.0 * M_PI) );
  bool ok = false;

  if (mappingType == EE_SingleDiscontinuityPoint) {
    rad *= MULT;
    if (rad > 1.0) {
      // use the bottom sphere and modify the radius, here
      // we have one point of the discontinuity.
      rad = sqrt( (1.0/(MULT*MULT - 1.0))*(MULT*MULT - rad*rad) );
      // now we have radius again in the range <0.0; 1.0>
      assert(rad >= 0.0);
      assert(rad <= 1.0);
      ok = bottomHemisphere->SampleOne(rad, phi, index);
      topHemisphereFlag = false;
    }
    else {
      // use top sphere, no change of the radius
      ok = topHemisphere->SampleOne(rad, phi, index);
      topHemisphereFlag = true;
    }
    return false; // ok
  }

  if (mappingType == EE_TwoSquaresNoContinuity) {
    // random selection based on the total power
    float rndValue = RandomValue(0.0, totalPower);
    if (rndValue <= totalPowerTop) {
      // use the top hemisphere
      ok = topHemisphere->SampleOne(rad, phi, index);
      topHemisphereFlag = true;
    }
    else {
      // use the bottom hemisphere
      ok = bottomHemisphere->SampleOne(rad, phi, index);
      topHemisphereFlag = false;
    }    
  }

  if (mappingType == EE_TwoSquaresContinuousEdge) {
    // the stretching in one dimension, from a square we
    // do two squares with an adjacent edge. The decision
    // which hemisphere is accessed must be done above
    // based on the original random sample in 2D unit square.
        
    if (topHemisphereFlag) {
      // use the top hemisphere
      ok = topHemisphere->SampleOne(rad, phi, index);
      //topHemisphereFlag = true;
    }
    else {
      // use the bottom hemisphere
      ok = bottomHemisphere->SampleOne(rad, phi, index);
      //topHemisphereFlag = false;
    }

    
  }
  
  assert(!ok);
  
  return false; // ok
}

// 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
//   topHemisphere 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
CSphereGonio::SampleOneInverse(double &rad, double &phi,
			       bool topHemisphereFlag, int index)
{
  bool ok;

  if (topHemisphereFlag) {
    ok = topHemisphere->SampleOneInverse(rad, phi, index);

    if (mappingType == EE_SingleDiscontinuityPoint) {
      // we want to remap it back to the original range
      rad /= MULT;
    }

  }
  else {
    // bottom hemisphere
    ok = bottomHemisphere->SampleOneInverse(rad, phi, index);

    if (mappingType == EE_SingleDiscontinuityPoint) {
      // we want to remap it back to the original range
      rad = sqrt(MULT*MULT - rad * rad * (MULT*MULT-1.0));
      assert(rad >= 1.0);
      assert(rad <= MULT);
      // get to the range <0.0; 1.0>
      rad /= MULT;
    }
  }
  
  assert(!ok);

  return false; // ok
}
