// ===================================================================
// $Id$
//
// functions.h
//    Different useful functions
//
//
// Initial coding by Polina Kondratieva (PK), January 2004.
// ===================================================================
#ifndef __FUNCTIONS_H__
#define __FUNCTIONS_H__

// BRDFLIB library
#include "configh.h"

// BRDFLIB headers for Database-useful functions
#include "color.h"
#include "vector3.h"
// Library for hearest neighbouring computation on n-dimensional datasets
#include <ANN/ANN.h>
#include <iostream>
#include <fstream>
#include <cstdio>
#include <vector>

__BEGIN_BRDFLIB_HEADER

using namespace std;

class CMatrix4;
class CRay;


/*
  Reference:
  ------------
  1. Philip Dutre, Global Illumination Compendium. Program of Computer
     Graphics, Cornell University.
  2. Arvo J., Stratified Sampling of Spherical Triangles, SIGGRAPH 1995.
  3. BTF Database (Dana, etc.):
     http://www1.cs.columbia.edu/CAVE/curet/.index.html
*/


// Make local coordinate system from the two given vectors (normal and
// direction if anisotropy). The normal vectors is assumed to be Z-Axis in the
// local CS and direction of anisotropy is assumed to be X-Axis.
//=============================================================================
static void MakeLocalCS(CVector3D &baseZ, CVector3D &baseX, CVector3D &locAxisX,
                        CVector3D &locAxisY, CVector3D &locAxisZ)
{
  locAxisZ = baseZ;
  locAxisX = baseX;
  locAxisY = CrossProd(locAxisZ, locAxisX);
  locAxisX = CrossProd(locAxisY, locAxisZ);
  locAxisX.Normalize();
  locAxisY.Normalize();
  locAxisZ.Normalize();
}
//=============================================================================

// Make arbitrary local coordinate system by the given vector (normal). The
// given vector is assumed to be Z-Axis in the local CS.
//=============================================================================
static void MakeArbitraryLCS(CVector3D &baseZ,    CVector3D &locAxisX,
                            CVector3D &locAxisY, CVector3D &locAxisZ)
{
  CVector3D Xaxis(1,0,0), Yaxis(0,1,0), Zaxis(0,0,1);
  locAxisX = CVector3D(1,0,0);
  locAxisY = CVector3D(0,1,0);
  locAxisZ = CVector3D(0,0,1);

  float dotX, dotY, dotZ, max = 0;
  int maxi = 0;
  dotX = fabs(DotProd(Xaxis, baseZ));
  dotY = fabs(DotProd(Yaxis, baseZ));
  dotZ = fabs(DotProd(Zaxis, baseZ));

  // 'base' is given in the World coordinate system. Choose the nearest to
  // 'base' axis 'i' in the World coordinate system and set axis 'i' in
  // local coordinate system to 'base'
  if (dotZ >= dotX) {
     max = dotZ;
     maxi = 2;        // base vector tends to Z-axis
  }
  else {
     max = dotX;
     maxi = 0;        // base vector tends to X-axis
  }

  if (max < dotY) {
     max = dotY;
     maxi = 1;        // base vector tends to Y-axis
  }

  CVector3D tmp(0,0,1);
  locAxisZ = baseZ;

  // Find other axes around the given Z-axis (locAxisZ = baseZ)
  switch(maxi)
  {
    case 1:
    case 2:
       tmp = CrossProd(locAxisZ, Xaxis);
       locAxisX  = CrossProd(tmp, locAxisZ);
       locAxisY  = CrossProd(locAxisZ, locAxisX);
       locAxisY.Normalize();
       locAxisX.Normalize();
       break;

    default:
       tmp = CrossProd(locAxisZ, Yaxis);
       locAxisY  = CrossProd(tmp, locAxisZ);
       locAxisX  = CrossProd(locAxisY, locAxisZ);
       locAxisY.Normalize();
       locAxisX.Normalize();
       break;
  }
}
//=============================================================================


// Translate Cartesian coordinates into Polar coordinates. Normal should be
// specified to set the orientation of the sphere which we need to construct
// to get polar coordinates. In this way normal is the vertical axis for
// this sphere.
//=============================================================================
static void CartesianToPolar(CVector3D &normal, CVector3D &xyz, float &phi,
                             float &theta, float &radius)
{
  CVector3D locX(1,0,0), locY(0,1,0), locZ(0,0,1);
  MakeArbitraryLCS(normal, locX, locY, locZ);

  CVector3D loc_xyz(0,0,1);
  loc_xyz[0] = DotProd(xyz, locX);
  loc_xyz[1] = DotProd(xyz, locY);
  loc_xyz[2] = DotProd(xyz, locZ);

  radius = Magnitude(loc_xyz);
  assert(radius > 0);
  theta = acos(loc_xyz[2]/radius);
  CVector3D tmp = loc_xyz;
  tmp[2] = 0;
  Normalize(tmp);

  if(tmp[1] >= 0 && tmp[0] >= 0)
  {
    if(loc_xyz[1] < 1e-6 && loc_xyz[0] < 1e-6)
    {
      phi = 0;
    }
    else
    {
      if(tmp[0] < 1e-6)
        phi = M_PI/2;
      else
        phi = atan(loc_xyz[1]/loc_xyz[0]);
    }
  }

  if(tmp[0] < 0)
  {
    if(fabs(loc_xyz[1]) < 1e-6 && fabs(loc_xyz[0]) < 1e-6)
    {
      if(tmp[0] >= 1.0 - 1e-6) phi = 2.0*M_PI;
      else                     phi = 0;
    }
    else
    {
      phi = atan(loc_xyz[1]/loc_xyz[0]) + M_PI;
    }
  }

  if(tmp[0] >= 0 && tmp[1] < 0)
  {
    if(loc_xyz[0] < 1e-6 && fabs(loc_xyz[1]) < 1e-6)
    {
      phi = 2.0*M_PI;
    }
    else
    {
      if(tmp[0] < 1e-6)
        phi = 3.0*M_PI/2;
      else
        phi = atan(loc_xyz[1]/loc_xyz[0]) + 2.0* M_PI;
    }
  }
}

// Translate Polar coordinates into Cartesian coordinates. Normal should be
// specified to set the orientation of the sphere which we need to construct
// to get polar coordinates. In this way normal - is the vertical axis for
// this sphere.
//=============================================================================
static void PolarToCartesian(CVector3D &normal, float &phi, float &theta,
                             float &radius, CVector3D &xyz)
{
  assert(phi <= (2.0*M_PI + 1e-6) && phi >= (0 - 1e-6));
  assert(theta <= (M_PI + 1e-6) && theta >= (0 - 1e-6));
  assert(radius > 0);

  CVector3D tmp(0,0,1);
  tmp[0] = radius * sin(theta) * cos(phi);
  tmp[1] = radius * sin(theta) * sin(phi);
  tmp[2] = radius * cos(theta);

  CVector3D U = ArbitraryNormal(normal);     // U is perpendicular to N
  CVector3D V = CrossProd(normal, U);        // orthonormal base with N and U

  xyz = tmp[0] * U + tmp[1] * V + tmp[2] * normal;
}


// Function for evaluation of integral Int[(cos(Q,M))^n, Q] over digon Q2
// where Q is outgoing direction, M is Cosine lobe axis.
// Cosine lobe is nonzero only on digon Q2.
// (See Compendium, p.17, initially presented by Arvo)
//=============================================================================
static float IntegralCosN(const CVector3D &axis2,
                          const CVector3D &normal, float &nnn)
{
   CVector3D locX(1,0,0), locY(0,1,0), locZ(0,0,1);
   CVector3D norm = normal;
   // Compute local coordinate system
   MakeArbitraryLCS(norm, locX, locY, locZ);

   // Translate normal and axis2 into the local coordinate system
   CVector3D localInc, localNorm;
   localInc[0] = DotProd(locX, axis2);
   localInc[1] = DotProd(locY, axis2);
   localInc[2] = DotProd(locZ, axis2);

   localNorm[0] = DotProd(locX, norm);
   localNorm[1] = DotProd(locY, norm);
   localNorm[2] = DotProd(locZ, norm);

   float Fi, sinus, cosinus, sum = 0;

   int v_nnn = (int)nnn;
   int iii;

   // See Compendium...
   if(v_nnn % 2 == 0)
   {
      Fi = 2.0;
      cosinus = DotProd(localInc, localNorm);

      if(fabs(cosinus) >= (1.0 - 1e-06))
      {
        sinus = 0;
        cosinus /= fabs(cosinus);
      }
      else
        sinus = sqrt(1.0f - cosinus*cosinus);

      sum += Fi*sinus;

      for(iii = 3; iii <= nnn-1; iii+=2)
      {
         Fi *= (iii-1.0f+(nnn-v_nnn))/((float)iii+(nnn-v_nnn));
	 sum += Fi * pow(sinus, ((float)iii+(nnn-v_nnn)));
      }
      sum *= cosinus;
      sum += 2.0*(M_PI - acos(cosinus));
      sum /= (nnn+1.0);
   }
   else
   {
      Fi = M_PI;
      cosinus = DotProd(localInc, localNorm);

      if(fabs(cosinus) >= (1.0 - 1e-06))
      {
        sinus = 0;
        cosinus /= fabs(cosinus);
      }
      else
        sinus = sqrt(1.0f - cosinus*cosinus);
      sum += Fi;

      for(iii = 2; iii <= nnn-1.0; iii+=2)
      {
         Fi *= (iii-1.0f+(nnn-v_nnn))/((float)iii+(nnn-v_nnn));
	 sum += Fi * pow(sinus, ((float)iii+(nnn-v_nnn)));
      }
      sum *= cosinus;
      sum += M_PI;
      sum /= (nnn+1.0);
   }
   return sum;
}
//=============================================================================


// Function for evaluation of integral Int[(cos(Q,M))^n*Cos(N,Q), Q] over digon
// Q2 where Q is outgoing direction, M is Cosine lobe axis.
// Cosine lobe is nonzero only on digon Q2.
// (See Compendium, p.17, initially presented by Arvo)
// (Cosine lobe axis (axis2) is vector m (see p.16 Compendium))
//=============================================================================
static float IntegralCosNCos(const CVector3D  &axis2,
                             const CVector3D &normal, float &nnn)
{
   CVector3D locX(1,0,0), locY(0,1,0), locZ(0,0,1);
   CVector3D norm = normal;
   // Compute local coordinate system
   MakeArbitraryLCS(norm, locX, locY, locZ);

   // Translate normal and axis2 into the local coordinate system
   CVector3D localInc, localNorm;
   localInc[0] = DotProd(locX, axis2);
   localInc[1] = DotProd(locY, axis2);
   localInc[2] = DotProd(locZ, axis2);

   localNorm[0] = DotProd(locX, norm);
   localNorm[1] = DotProd(locY, norm);
   localNorm[2] = DotProd(locZ, norm);

   float deg, angle, term, sinus, cosinus, sum = 0, result;
   int v_nnn;

   v_nnn = (int)nnn;
   cosinus = DotProd(localInc, localNorm);

   if(fabs(cosinus) >= (1.0 - 1e-06))
   {
     sinus = 0;
     cosinus /= fabs(cosinus);
   }
   else
     sinus = sqrt(1.0f - cosinus*cosinus);

   if(v_nnn % 2 == 0)
   {
	  term = M_PI_2;
      angle = M_PI_2;
      deg = 0.0f + (nnn - v_nnn);
   }
   else
   {
      term = sinus;
      angle = M_PI - acos(cosinus);
      deg = 1.0f + (nnn - v_nnn);
   }

   while(deg <= nnn-2.0f)
   {
      sum +=  term;
      term *= sinus * sinus * (deg+1.0f)/(deg+2.0f);
      deg += 2.0f;
   }
   result = 2.0f * (term + cosinus*angle + cosinus*cosinus*sum)/(nnn+2.0f);
   return result;
}
//=============================================================================



/******************************************************************************
  Functions useful for:
  1) implementation of BRDF reconstruction from the goniometric web
     diagram or from the database.
  2) importance sampling of BTFs (BRDF written as a texture and
     available from some databases)
******************************************************************************/
// ========================================================================
// $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
//

// For debugging purposes
//#define _DEBUG

//==========================================================================
//==========================================================================
inline double
sqr(double a) { return a * a; }

// -----------------------------------------------------------------------
// Generate random variable in the range <a, b>
inline double
RandomValue(double a, double b) {
  double range = fabs(a - b);
  return ((double)rand() / RAND_MAX) * range + ((a < b) ? a : b);
}

// -----------------------------------------------------------------------
// N-th power of base in integer
inline int
intpow(int base, int exponent) {
  int res = 1;
  for(int i=0;i<exponent;i++)
    res *= base;
  return res;
}

// -----------------------------------------------------------------------
/**
   Generate random point in unit disk in 2D using concentric map
   with density p(x) = 1/A.

   @param r,phi   the generated point in 2D spherical coordinates
   (r,phi), radius 'r' in <0,1>, angle 'phi' in <0,2*PI> (output)

   @param wx,wy   two random numbers between 0 and 1 (input)

   Taken from JGT, paper by Peter Shirley, Kenneth Chiu:
   "A Low Distortion Map Between Disk and Square"
*/
inline void
to_unit_disk(double &r, double &phi,
	     const double &wx, const double &wy)
{
  assert(wx>=0.0);
  assert(wx<=1.0);
  assert(wy>=0.0);
  assert(wy<=1.0);

  double vx = -1.0 + 2.0 * wx;
  double vy = -1.0 + 2.0 * wy;

  if (vx > -vy) {
    if (vx > vy) {
      r = vx; phi = M_PI_4 * (vy / vx);
      // This modification is required to get range <0;2*M_PI>
      if (phi < 0.0)
	phi += 2.0 * M_PI;
    }
    else {
      r = vy; phi = M_PI_4 * (2.0 - (vx / vy) );
    }
  }
  else {
    // It holds: vx <= -vy
    if (vx < vy) {
      r = -vx; phi = M_PI_4 * (4.0 + (vy/vx) );
    }
    else {
      r = -vy;
      if (vy != 0.0)
	phi = M_PI_4 * (6.0 - (vx/vy) );
      else
	phi = 0.0;
    }
  }

  assert(r <= 1.00001);
  assert(phi >= -0.000001);
  assert(phi <= 2.0*M_PI+0.000001);
  return;
};

// -----------------------------------------------------------------------
/**
   The inverse function to the above - takes 'radius' and 'phi' angle
   and computes

   @param wx,wy --  the position in the square [0-1]x[0-1] (output)

   @param r,phi --  the generated point in 2D spherical coordinates
   (r,phi), radius 'r' in <0,1>, angle 'phi' in <0,2*PI> (input)

   Taken from JGT, paper by Peter Shirley, Kenneth Chiu:
   "A Low Distortion Map Between Disk and Square"
*/
inline void
from_unit_disk(double &x, double &y,
	       double r, double phi)
{
  double a, b;
  assert(phi >= 0.0);
  assert(phi <= 2.0*M_PI);

  if (phi > M_PI * 7 / 4)
    phi -= 2*M_PI;
  // phi is now in range [-pi/4,7pi/4]

  if ( phi < 3*M_PI_4 ) {
    if ( phi < M_PI_4) { // region 1
      a = r;
      b = phi * a / (M_PI_4);
    }
    else {     // region 2
      b = r;
      a = -(phi - M_PI / 2.0) * b / (M_PI_4);
    }
  }
  else {
    if ( phi < 5 * M_PI_4 ) { // region 3
      a = -r;
      b = (phi - M_PI) * a / (M_PI_4);
    }
    else { // region 4
      b = -r;
      a = -(phi - 3*M_PI / 2.0) * b / (M_PI_4);
    }
  }

  x = (a + 1.0) / 2.0;
  y = (b + 1.0) / 2.0;
  return;
}
//==========================================================================
//==========================================================================


/******************************************************************************
  Functions useful for:
  1) reading of the BTF/BRDF databases.
  2) construction of specific data structures for BTF/BRDF data
******************************************************************************/

// Database-data functions
//-------------------------------------------------------------------------
// Contains 4D coordinates and correspondent BRDF values (average and RGB)
//========================================================================
struct S_entryDB
{
   ANNpoint p;
   CValues  brdf;
   float    avg;
};

inline ostream& operator<< (ostream &s, const S_entryDB &entry_db)
{
  s << "(";
  assert((int)sizeof(entry_db.p) <= 4);
  for(int i = 0; i < ((int)sizeof(entry_db.p)); i++)
  {
     s << (entry_db.p)[i];
     if(i < (((int)sizeof(entry_db.p))-1)) s << ", ";
  }
  return s << ") -> " << entry_db.brdf << ", avg: " << entry_db.avg;
}

// Polar coordinates (Phi and Theta)
//========================================================================
struct S_polarCoords
{
   float theta;
   float phi;
};

inline ostream& operator<< (ostream &s, const S_polarCoords &polarCoords)
{
  return s << "(" << polarCoords.theta << ", " << polarCoords.phi << ")";
}

typedef vector<S_entryDB>  T_dataDB;

//========================================================================
inline float DegreeToRad(float in_degree)
{
   return in_degree*M_PI/180.0;
}

inline float RadToDegree(float in_rad)
{
   return in_rad*180.0/M_PI;
}


// Some functions with usage of ANN library http://www.cs.umd.edu/~mount/ANN
//==========================================================================
//==========================================================================

// Read point in ANN format (false on EOF)
//========================================================================
inline ANNbool readPt(ifstream &in, ANNpoint p, int &rdim)
{
    int num;
    if(!(in >> num)) return ANNfalse;

    for (int i = 0; i < rdim; i++) {
	if(!(in >> p[i])) return ANNfalse;
    }
    return ANNtrue;
}

// Print ANN point
inline void printPt(ostream &out, ANNpoint p, int &rdim)
{
    out << "(" << p[0];
    for (int i = 1; i < rdim; i++) {
	out << ", " << p[i];
    }
    out << ")";

    out << " -> (" << RadToDegree(p[0]);
    for (int i = 1; i < rdim; i++) {
	out << ", " << RadToDegree(p[i]);
    }
    out << ")\n";
}

// Lookup for the data taken from a table and stored into an array
//========================================================================
inline int
GetIndex(int maxTheta, int theta, int phi) {
  return theta + maxTheta * phi;
}


// Helpful functions for reading BRDF databases
//==========================================================================
//==========================================================================

// Read Dana`s database [3]
//========================================================================
// The difference between the isotropic and anisotropic data: isotropic
// data do not depend on Phi angle but on Theta angle. In this way, only
// #theta_angles grids with BRDF values in the nodes is enough to get
// BRDF value for any incoming and outgoing directions.
//
// Input:
// ---------------------------------------------------------------------
// sample_num - number of sample in the database (for example, there 61
//              samples in Dana's database).
// measur_num - number of measurments for each sample, number of
//              combinations of incoming and outgoing directions (there
//              are 205 combinations in Dana's database)
// names_db   - contains the names of the files in DB to be read
//   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//   names_bd[0] should be the name of the file containing polar coordinates
//               for incoming/outgoing directions
//   grey-scale case:
//   ----------------
//      names_bd[1] should be the name of the file containing BRDF values
//                  in grey-scale. It is clear, that in this case average
//                  BRDF value is the same as the corresponding grey-scale
//                  single value
//   RGB-color case:
//   ----------------
//      names_bd[1] should be the name of the file containing BRDF values
//                  for red color
//      names_bd[2] should be the name of the file containing BRDF values
//                  for green color
//      names_bd[3] should be the name of the file containing BRDF values
//                  for blue color
//      names_bd[4] should be the name of the file containing average
//                  BRDF values (for example, for albedo computation)
//   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// Output:
// ---------------------------------------------------------------------
// BRDF_i     - contains the entries consisted of 4D-coordinates (incoming
//              and outgoing directions), corresponding BRDF values in
//              RBG format, and average values of BRDFs (read from DB).
// s_treeDB   - search tree to find the query 4D-coordinates in 4D data
//              read from the DB.
//========================================================================

//========================================================================
static ANNkd_tree *ReadDataBT(int &sample_num, int &measur_num, int &rdim,
                              vector<const char*> &names_db, T_dataDB &BRDF_i)
{
  assert((names_db.size() == 2) || (names_db.size() == 5));
  if((rdim != 3) && (rdim != 4))
  {
    cout << "Note: dimension of search tree can be only 3 (for isotropic\n";
    cout << "case) or 4 (for anisotropic case)\n";
    abort();
  }

  // In Dana et al.'s database there are measurments for 61 samples. There
  // are 205 measurments for each sample (205 incoming/outgoing angles
  // combinations).
  //----------------------------------------------------------------------
  assert(((sample_num > 0) && (sample_num < 62)) &&
         ((measur_num > 0) && (measur_num < 206)));

  char dummy_name[] = "dummyDBData.txt";
  if(names_db.size() == 2)
  {
    names_db.push_back(dummy_name);
    names_db.push_back(dummy_name);
    names_db.push_back(dummy_name);
  }

  std::ifstream in_tb( names_db[0]);     // incoming/outgoing angles
  std::ifstream in_red( names_db[1] );   // red component of BRDF
  std::ifstream in_green( names_db[2] ); // green component of BRDF
  std::ifstream in_blue( names_db[3] );  // blue component of BRDF
  std::ifstream in_avg( names_db[4] );   // average of BRDF components

  // If one/more specified files do not exist - exit
  //----------------------------------------------------------------------
  if( !in_tb || !in_red || (!in_green && (names_db.size() == 5)) ||
      (!in_blue && (names_db.size() == 5)) ||
      (!in_avg && (names_db.size() == 5)))
  {
    cout << "No such file or directory ";
    if(!in_tb) cout << "\"" << names_db[0] << "\" ";
    if(!in_red) cout << "\"" << names_db[1] << "\" ";
    if(!in_green && (names_db.size() == 5))
      cout << "\"" << names_db[2] << "\" ";
    if(!in_blue && (names_db.size() == 5))
      cout << "\"" << names_db[3] << "\" ";
    if(!in_avg && (names_db.size() == 5))
      cout << "\"" << names_db[4] << "\" ";

    cout << "\n";
    abort();
  }

  // Read table of angles and put them into the 4D kd-tree
  //----------------------------------------------------------------------
  // 'Data points' are polar coordinates here. In this case we get 4D
  // structure: 2 incoming coordinates and 2 outgoing coordinates.
  // Put these coordinates into the 1D array according to the
  // requirements of ANN-library

  // Maximum number of data points: 205 combinations of incoming/outgoing
  // direction with 2 polar coordinates for each direction. Of course,
  // it is possible to take less than 205 combinations.

  // To get correct distances on the boundaries of the data set we should
  // replicate the data for Phi angles: 0 and 2Pi - is the same => We should
  // put Phi angles on both left and right sides of original range of Phi
  // angles as following:
  //
  // Theta1  |------------|------------|------------|
  //         |            |  Original  |            |
  // Theta0  |------------|------------|------------|
  //         -2Pi         0           2Pi          4Pi
  int m_pts = measur_num;                 // measured number of points
  int m_pts_r = m_pts*3;                  // number of points after replication
  int dim_tmp = 4;                        // dimension of point to be read

  int           n_pts;                    // actual number of data points
  ANNpointArray data_pts;                 // data points
  ANNpoint point_tmp;                     // temporary points
  data_pts = annAllocPts(m_pts_r, rdim);   // allocate data points
  point_tmp = annAllocPt(dim_tmp);
  n_pts = 0;

  //std::cout << "Data Points:\n";

  // Read the table of incoming/outgoing polar angles from the table
  //----------------------------------------------------------------------
  while (n_pts < m_pts && readPt(in_tb, point_tmp, dim_tmp)) {
     //printPt(cout, point_tmp, dim_tmp);
     // Replicate the data
     for(int tt = 0; tt < 3; tt++)
     {
       data_pts[n_pts + m_pts*tt][0] = point_tmp[0];
       data_pts[n_pts + m_pts*tt][1] = point_tmp[1] + (tt-1)*2.0*M_PI;
       data_pts[n_pts + m_pts*tt][2] = point_tmp[2];

       if(rdim == 4)
         data_pts[n_pts + m_pts*tt][3] = point_tmp[3] + (tt-1)*2.0*M_PI;
     }
     n_pts++;
  }

  assert(n_pts == m_pts);
  //----------------------------------------------------------------------

  char tmp_array[2000];
  int j;

  // Read the table of BRDF values for sample number sample_num
  //----------------------------------------------------------------------
  // Grey BRDF values
  if(names_db.size() == 2)
  {
    // Read first sample_num-1 lines from the file to reach required sample
    for(j = 0; j < sample_num-1; j++)
      in_red.getline(tmp_array, sizeof(tmp_array), '\n');

    int num_r;
    in_red >> num_r;

    for(j = 0; j < measur_num; j++)
    {
      float tmp_value = 0;
      in_red >> tmp_value;
      CValues color_v(CValues::EE_Weigth1);
      color_v.Set(0, tmp_value);
      S_entryDB entry_db;
      entry_db.brdf = color_v;
      entry_db.p = data_pts[j+m_pts];
      BRDF_i.push_back(entry_db);  // Fill the table of BRDF values
    }
  }
  else  // RGB BRDF values
  {
    // Read first sample_num-1 lines from the file to reach required sample
    for(j = 0; j < sample_num-1; j++)
    {
       in_red.getline(tmp_array, sizeof(tmp_array), '\n');
       in_green.getline(tmp_array, sizeof(tmp_array), '\n');
       in_blue.getline(tmp_array, sizeof(tmp_array), '\n');
       in_avg.getline(tmp_array, sizeof(tmp_array), '\n');
    }
    int num_r, num_green, num_blue, num_avg;

    // Read sample number (it should be required sample number)
    in_red >> num_r;   in_green >> num_green;   in_blue >> num_blue;
    in_avg >> num_avg;
    assert((num_r == num_green) && (num_green == num_blue)
           && (num_r == num_avg));

    for(j = 0; j < measur_num; j++)
    {
      float tmp_valueR = 0, tmp_valueG = 0, tmp_valueB = 0, tmp_avg = 0;
      in_red >> tmp_valueR;
      in_green >> tmp_valueG;
      in_blue >> tmp_valueB;
      in_avg  >> tmp_avg;
      CValues color_v(tmp_valueR, tmp_valueG, tmp_valueB);
      color_v.SetRepresentation(CValues::EE_Weigth3);
      S_entryDB entry_db;
      entry_db.brdf = color_v;
      entry_db.p = data_pts[j+m_pts];
      entry_db.avg = tmp_avg;
      BRDF_i.push_back(entry_db);  // Fill the table of BRDF values
    }
  }

  // build search structure for the 4D coordinates
  ANNkd_tree *the_tree = new ANNkd_tree( data_pts, // the data points
                                         m_pts_r,  // number of points
                                         rdim);    // dimension of space
  return the_tree;
}


// The same as above, but do not compute the searching tree
//========================================================================
static void ReadData(int &sample_num, int &measur_num, int &rdim,
                     vector<const char*> names_db, T_dataDB &BRDF_i)
{
  assert((names_db.size() == 2) || (names_db.size() == 5));

  assert(((sample_num > 0) && (sample_num < 62)) &&
         ((measur_num > 0) && (measur_num < 206)));

  if((rdim != 3) && (rdim != 4))
  {
    cout << "Note: dimension of search tree can be only 3 (for isotropic\n";
    cout << "case) or 4 (for anisotropic case)\n";
    abort();
  }

  char dummy_name[] = "dummyDBData.txt";
  if(names_db.size() == 2)
  {
    names_db.push_back(dummy_name);
    names_db.push_back(dummy_name);
    names_db.push_back(dummy_name);
  }

  std::ifstream in_tb   ( names_db[0] );
  std::ifstream in_red  ( names_db[1] );
  std::ifstream in_green( names_db[2] );
  std::ifstream in_blue ( names_db[3] );
  std::ifstream in_avg  ( names_db[4] );

  if( !in_tb || !in_red || (!in_green && (names_db.size() == 5)) ||
      (!in_blue && (names_db.size() == 5)) ||
      (!in_avg && (names_db.size() == 5)))
  {
    cout << "No such file or directory ";
    if(!in_tb) cout << "\"" << names_db[0] << "\" ";
    if(!in_red) cout << "\"" << names_db[1];
    if(!in_green && (names_db.size() == 5))
      cout << "\"" << names_db[2] << "\" ";
    if(!in_blue && (names_db.size() == 5))
      cout << "\"" << names_db[3] << "\" ";
    if(!in_avg && (names_db.size() == 5))
      cout << "\"" << names_db[4] << "\" ";

    cout << "\n";
    return;
  }

  // Read table of angles
  //----------------------------------------------------------------------
  int m_pts = measur_num;                 // measured number of points
  int m_pts_r = m_pts*3;                  // number of points after replication
  int dim_tmp = 4;                        // dimension of point to be read

  int           n_pts;                    // actual number of data points
  ANNpointArray data_pts;                 // data points
  ANNpoint point_tmp;                     // temporary points
  data_pts = annAllocPts(m_pts_r, rdim);  // allocate data points
  point_tmp = annAllocPt(dim_tmp);
  n_pts = 0;

  //std::cout << "Data Points:\n";

  // Read the table of incoming/outgoing polar angles from the table
  while (n_pts < m_pts && readPt(in_tb, point_tmp, dim_tmp)) {
     //printPt(cout, point_tmp, dim_tmp);
     // Replicate the data
     for(int tt = 0; tt < 3; tt++)
     {
       data_pts[n_pts + m_pts*tt][0] = point_tmp[0];
       data_pts[n_pts + m_pts*tt][1] = point_tmp[1] + (tt-1)*2.0*M_PI;
       data_pts[n_pts + m_pts*tt][2] = point_tmp[2];

       if(rdim == 4)
         data_pts[n_pts + m_pts*tt][3] = point_tmp[3] + (tt-1)*2.0*M_PI;
     }
     n_pts++;
  }
  assert(n_pts == m_pts);
  //----------------------------------------------------------------------

  char tmp_array[2000];
  int j;

  // Read the table of BRDF values for sample number sample_num
  //----------------------------------------------------------------------
  if(names_db.size() == 2)
  {
    // Read first sample_num-1 lines from the file
    for(j = 0; j < sample_num-1; j++)
      in_red.getline(tmp_array, sizeof(tmp_array), '\n');

    int num_r;
    in_red >> num_r;

    for(j = 0; j < measur_num; j++)
    {
      float tmp_value = 0;
      in_red >> tmp_value;
      CValues color_v(CValues::EE_Weigth1);
      color_v.Set(0, tmp_value);
      S_entryDB entry_db;
      entry_db.brdf = color_v;
      entry_db.p = data_pts[j+m_pts];
      BRDF_i.push_back(entry_db);
    }
  }
  else
  {
    // Read first sample_num-1 lines from the file
    for(j = 0; j < sample_num-1; j++)
    {
       in_red.getline(tmp_array, sizeof(tmp_array), '\n');
       in_green.getline(tmp_array, sizeof(tmp_array), '\n');
       in_blue.getline(tmp_array, sizeof(tmp_array), '\n');
       in_avg.getline(tmp_array, sizeof(tmp_array), '\n');
    }
    int num_r, num_green, num_blue, num_avg;

    in_red >> num_r;   in_green >> num_green;   in_blue >> num_blue;
    in_avg >> num_avg;
    assert((num_r == num_green) && (num_green == num_blue)
           && (num_r == num_avg));

    for(j = 0; j < measur_num; j++)
    {
      float tmp_valueR = 0, tmp_valueG = 0, tmp_valueB = 0, tmp_avg = 0;
      in_red >> tmp_valueR;
      in_green >> tmp_valueG;
      in_blue >> tmp_valueB;
      in_avg  >> tmp_avg;
      CValues color_v(tmp_valueR, tmp_valueG, tmp_valueB);
      color_v.SetRepresentation(CValues::EE_Weigth3);
      S_entryDB entry_db;
      entry_db.brdf = color_v;
      entry_db.p = data_pts[j+m_pts];
      entry_db.avg = tmp_avg;
      BRDF_i.push_back(entry_db);
    }
  }
  return;
}




// Database from Bonn (BTF database). We desided to get an average RGB color
// of the images and use it as BRDF value. Format of file which we get from
// the database is as following:
//   -------------------------------------------------------------------
//    ID   theta_v   phi_v  theta_l   phi_l   RGB_BRDF   Average_BRDF
//   -------------------------------------------------------------------
//
// where ID - is timestep embeded into the file name of the BTF Database
//       theta_v, phi_v - angles of measurments for viewing direction
//                        (outgoing direction)
//       theta_v, phi_v - angles of measurments for incoming light direction
//       RGB_BRDF - averaged color of specified image from the Database
//       average_BRDF - luminance of RGB_BRDF color
//
// File average.pl is used to:
//  read the images from the database,
//  average the color of the images,
//  writting the results (RGB color and Luminance) with corresponding
//    outgoing/incoming angles of measurments into the one file.
//========================================================================
static ANNkd_tree *ReadDataBonnBT(int &sample_num, int &measur_num, int &rdim,
                              vector<const char*> &names_db, T_dataDB &BRDF_i)
{
  sample_num = 1;
  assert((int)names_db.size() == 1);
  if((rdim != 3) && (rdim != 4))
  {
    cout << "Note: dimension of search tree can be only 3 (for isotropic\n";
    cout << "case) or 4 (for anisotropic case)\n";
    abort();
  }

  // In database created in Bonn there are measurments for some samples. There
  // are 6561 measurments for each sample (6561 combinations of incoming and
  // outgoing angles).
  //----------------------------------------------------------------------
  assert((measur_num > 0) && (measur_num < 6562));

  std::ifstream in_tb( names_db[0]);

  // If one/more specified files do not exist - exit
  //----------------------------------------------------------------------
  if( !in_tb )
  {
    cout << "No such file or directory ";
    if(!in_tb) cout << "\"" << names_db[0] << "\" ";
    cout << "\n";
    abort();
  }

  // Read table of angles and put them into the 4D kd-tree
  //----------------------------------------------------------------------
  // 'Data points' are polar coordinates here. In this case we get 4D
  // structure: 2 incoming coordinates and 2 outgoing coordinates.
  // Put these coordinates into the 1D array according to the
  // requirements of ANN-library

  // Maximum number of data points: 205 combinations of incoming/outgoing
  // direction with 2 polar coordinates for each direction. Of course,
  // it is possible to take less than 205 combinations.

  // To get correct distances on the boundaries of the data set we should
  // replicate the data for Phi angles: 0 and 2Pi - is the same => We should
  // put Phi angles on both left and right sides of original range of Phi
  // angles as following:
  //
  // Theta1  |------------|------------|------------|
  //         |            |  Original  |            |
  // Theta0  |------------|------------|------------|
  //         -2Pi         0           2Pi          4Pi
  int m_pts = measur_num;                 // measured number of points
  int m_pts_r = m_pts*3;                  // number of points after replication
  int dim_tmp = 4;                        // dimension of point to be read

  //int           n_pts;                    // actual number of data points
  ANNpointArray data_pts;                 // data points
  ANNpoint point_tmp;                     // temporary points
  data_pts = annAllocPts(m_pts_r, rdim);   // allocate data points
  point_tmp = annAllocPt(dim_tmp);

  int j;

  // Read the table of BRDF values for sample number sample_num
  //----------------------------------------------------------------------
  for(j = 0; j < measur_num; j++)
  {
    float theta_out = 0, phi_out = 0, theta_inc = 0, phi_inc = 0;
    float tmp_valueR = 0, tmp_valueG = 0, tmp_valueB = 0, tmp_avg = 0;
    int ID_num = 0;

    in_tb >> ID_num;
    in_tb >> theta_out;
    in_tb >> phi_out;
    in_tb >> theta_inc;
    in_tb >> phi_inc;

    for(int tt = 0; tt < 3; tt++)
    {
      // Angles in database are measured in degrees => we should
      // convert them into the radians.
      data_pts[j + m_pts*tt][0] = DegreeToRad(theta_out);
      data_pts[j + m_pts*tt][1] = DegreeToRad(phi_out) + (tt-1)*2.0*M_PI;
      data_pts[j + m_pts*tt][2] = DegreeToRad(theta_inc);

      if(rdim == 4)
        data_pts[j + m_pts*tt][3] = DegreeToRad(phi_inc) + (tt-1)*2.0*M_PI;
    }
    in_tb >> tmp_valueR;
    in_tb >> tmp_valueG;
    in_tb >> tmp_valueB;
    in_tb >> tmp_avg;
    CValues color_v(tmp_valueR, tmp_valueG, tmp_valueB);
    color_v.SetRepresentation(CValues::EE_Weigth3);
    S_entryDB entry_db;
    entry_db.brdf = color_v;
    entry_db.p = data_pts[j+m_pts];
    entry_db.avg = tmp_avg;
    BRDF_i.push_back(entry_db);  // Fill the table of BRDF values
  }

  // build search structure for the 4D coordinates
  ANNkd_tree *the_tree = new ANNkd_tree( data_pts, // the data points
                                         m_pts_r,  // number of points
                                         rdim);     // dimension of space
  return the_tree;
}




// Database from Bonn (BTF database). We desided to get an average RGB color
// of the images and use it as BRDF value. Format of file which we get from
// the database is as following:
//
// theta_view  phi_view  theta_light  phi_light  RGB_BRDF  Average_BRDF
//
// File average.pl is used to:
//  read the images from the database,
//  average the color of the images,
//  writting the results (RGB color and Luminance) with corresponding
//    outgoing/incoming angles of measurments into the one file.
//========================================================================
static void ReadDataBonn(int &sample_num, int &measur_num, int &rdim,
                              vector<const char*> &names_db, T_dataDB &BRDF_i)
{
  sample_num = 1;
  assert((int)names_db.size() == 1);
  if((rdim != 3) && (rdim != 4))
  {
    cout << "Note: dimension of search tree can be only 3 (for isotropic\n";
    cout << "case) or 4 (for anisotropic case)\n";
    abort();
  }

  // In database created in Bonn there are measurments for some samples. There
  // are 6561 measurments for each sample (6561 combinations of incoming and
  // outgoing angles).
  //----------------------------------------------------------------------
  assert((measur_num > 0) && (measur_num < 6562));

   std::ifstream in_tb( names_db[0]);

  // If one/more specified files do not exist - exit
  //----------------------------------------------------------------------
  if( !in_tb )
  {
    cout << "No such file or directory ";
    if(!in_tb) cout << "\"" << names_db[0] << "\" ";
    cout << "\n";
    abort();
  }

  // Read table of angles and put them into the 4D kd-tree
  //----------------------------------------------------------------------
  // 'Data points' are polar coordinates here. In this case we get 4D
  // structure: 2 incoming coordinates and 2 outgoing coordinates.
  // Put these coordinates into the 1D array according to the
  // requirements of ANN-library

  // Maximum number of data points: 205 combinations of incoming/outgoing
  // direction with 2 polar coordinates for each direction. Of course,
  // it is possible to take less than 205 combinations.

  // To get correct distances on the boundaries of the data set we should
  // replicate the data for Phi angles: 0 and 2Pi - is the same => We should
  // put Phi angles on both left and right sides of original range of Phi
  // angles as following:
  //
  // Theta1  |------------|------------|------------|
  //         |            |  Original  |            |
  // Theta0  |------------|------------|------------|
  //         -2Pi         0           2Pi          4Pi
  int m_pts = measur_num;                 // measured number of points
  int m_pts_r = m_pts*3;                  // number of points after replication
  int dim_tmp = 4;                        // dimension of point to be read

  //int           n_pts;                    // actual number of data points
  ANNpointArray data_pts;                 // data points
  ANNpoint point_tmp;                     // temporary points
  data_pts = annAllocPts(m_pts_r, rdim);   // allocate data points
  point_tmp = annAllocPt(dim_tmp);

  int j;

  // Read the table of BRDF values for sample number sample_num
  //----------------------------------------------------------------------
  for(j = 0; j < measur_num; j++)
  {
    float theta_out = 0, phi_out = 0, theta_inc = 0, phi_inc = 0;
    float tmp_valueR = 0, tmp_valueG = 0, tmp_valueB = 0, tmp_avg = 0;
    int ID_num = 0;

    in_tb >> ID_num;
    in_tb >> theta_out;
    in_tb >> phi_out;
    in_tb >> theta_inc;
    in_tb >> phi_inc;

    for(int tt = 0; tt < 3; tt++)
    {
      // Angles in database are measured in degrees => we should
      // convert them into the radians.
      data_pts[j + m_pts*tt][0] = DegreeToRad(theta_out);
      data_pts[j + m_pts*tt][1] = DegreeToRad(phi_out) + (tt-1)*2.0*M_PI;
      data_pts[j + m_pts*tt][2] = DegreeToRad(theta_inc);

      if(rdim == 4)
        data_pts[j + m_pts*tt][3] = DegreeToRad(phi_inc) + (tt-1)*2.0*M_PI;
    }
    in_tb >> tmp_valueR;
    in_tb >> tmp_valueG;
    in_tb >> tmp_valueB;
    in_tb >> tmp_avg;
    CValues color_v(tmp_valueR, tmp_valueG, tmp_valueB);
    color_v.SetRepresentation(CValues::EE_Weigth3);
    S_entryDB entry_db;
    entry_db.brdf = color_v;
    entry_db.p = data_pts[j+m_pts];
    entry_db.avg = tmp_avg;
    BRDF_i.push_back(entry_db);  // Fill the table of BRDF values
  }

  return;
}


// Different Binary search implementations
//==========================================================================
//==========================================================================

//------------------------------------------------------------------------
// Binary search for an index of array component to which the given input
// number is equal.
//   input:
//      vector array - specified array,
//      float number - given number to be found in the specified array
//   output:
//      int ->  the index in the specified array of the array component to
//              which the input number is equal
//------------------------------------------------------------------------
inline int BinSearchNumber(vector<float> &array, float number)
{
  int N = (int)array.size();
  if(N <= 0)
  {
     cout << "There are no elements in the array!\n";
     return -1;
  }

  if(number < array[0] || number > array[N-1])
  {
     cout << "There are no such a number in the tree\n";
     return -1;
  }

  if(number == array[0])   return 0;
  if(number == array[N-1]) return N-1;

  int size = 0;
  int i;
  i = (int)(N/2.0);

  int dummy = 0;
  int lastR = N;
  int lastL = 0;

  while(size >= 0)
  {
    if(number == array[i])
    {
       return i;
    }
    else
    {
       if(number >= array[i])
       {
          lastL = i;
	  size = lastR - lastL;
	  i = (int)(size/2.0) + i;
	  if(size == 1) i++;
       }
       else
       {
          lastR = i;
	  size = lastR - lastL;
	  i = i - (int)(size/2.0)-1;
       }
    }
    dummy++;
    if(dummy > N) break;
  }
  return i;
}



//------------------------------------------------------------------------
// Binary search for an interval to which the given input number in the
// specified input array is fitted.
//   input:
//     double *array - specified array,
//     int rN        - the length of the specified array
//     double number - given number to be found in the specified array
//   output:
//     int ->  the index in the specified array of the interval to which
//             the input number is fitted
//------------------------------------------------------------------------
inline int BinSearchInterval(double *array, int rN, double number)
{
  if(rN <= 0)
  {
     return -1;
  }

  if(number < array[0]-1e-06 || number > array[rN-1]+1e-06)
  {
     return -1;
  }

  if(number < array[1])
  {
     return 0;
  }

  if(number >= array[rN-1])
  {
     return rN-1;
  }

  int size = 0;
  int i;
  i = (int)(rN/2.0);

  int dummy = 0;
  int lastR = rN;
  int lastL = 0;

  while(size >= 0)
  {
    if(number >= array[i] && number < array[i+1])
    {
       return i;
    }
    else
    {
       if(number >= array[i])
       {
          lastL = i;
	  size = lastR - lastL;
	  i = (int)(size/2.0) + i;
	  if(size == 1) i++;
       }
       else
       {
          lastR = i;
	  size = lastR - lastL;
	  i = i - (int)(size/2.0)-1;
       }
    }
    dummy++;
    if(dummy > rN) break;
  }
  return i;
}



__END_BRDFLIB_HEADER

#endif // __FUNCTIONS_H__
