// ===================================================================
// $Id: phong.cpp,v 1.22 2003/11/04 15:27:18 havran Exp $
//
// phong.cpp
//     Implementation of Phong shading model.
//
// Licence and usage conditions in November 2003:
// This file can be used and modified only with written permission of
// Vlastimil Havran, Pohodli 27, 57001 Litomysl, the Czech Republic.
// This library is granted to be used by Polina Kondratieva within
// her MPI HIWI job working on BRDFLIB library. The usage condition
// may change in the next version of the library. At the moment,
// the library development is supported by Max-Planck-Institut for
// Informatik, Stuhlsatzenhausweg 85, 66123 Saarbruecken Germany
// See WWW: http://www.mpi-sb.mpg.de
// REPLACEMENT_STRING
//
// Initial coding by Nicholas Wilt,
// rewritten by Vlastimil Havran & Jaroslav Krivanek, 20000.
//
// Given to BRDFLIB project for the Polina Kondratieva's
// HIWI JOB assignment in November 2003.


// BRDFLIB headers
#include "config.h"
#include "db_brdf.h"
#include "rnd.h"
#include "texture.h"
#include "functions.h"
#include <iostream>

__BEGIN_BRDFLIB_SOURCE

// This variable is used to store intermediate result of random number
// generation if Selection Roulette is used.
double CDBDana_BRDFAnisotropic::aux_rnd = 0.0;
typedef vector<double> rnd_array;
rnd_array CDBDana_BRDFAnisotropic::diffuse_rnd;
rnd_array CDBDana_BRDFAnisotropic::specular_rnd;
rnd_array CDBDana_BRDFAnisotropic::refraction_rnd;

/**********************************************************************
 *                                                                    *
 *                 CDBDana_BRDFAnisotropic methods                    *
 *                                                                    *
 **********************************************************************/

// describes the surface to a given stream
void
CDBDana_BRDFAnisotropic::Describe(ostream &app, int ind) const
{
  app << "\n";
  app << "*********************************************";
  app << "*************************\n";
  indent(app, ind);
  app << "CDBDana_BRDFAnisotropic:\n";
  app << "*********************************************";
  app << "*************************\n\n";

  if (_db_color) {
    indent(app, ind+2);
    app << "  db_color = \t";
    _db_color->Describe(app, ind); app << "\n";
  }
  /*if(_brdf_i.size() > 0) {
    //indent(app, ind+2);
    app << "BRDF data for correspondent incoming/outgoing\n";
    app << "directions (" << _brdf_i.size() << " samples):\n";
    app << "-----------------------------------------------\n";
    for(int i = 0; i < (int)_brdf_i.size(); i++)
    {
      app << "[" << i << "] ";
      for(int j = 0; j < (sizeof(_brdf_i[i].p)/sizeof((_brdf_i[i].p)[0])); j++)
        app << (_brdf_i[i].p)[j];
      app << " -> " << _brdf_i[i].brdf << "\n";
    }
  }*/
  PrintDB_SearchTree(cout);
  PrintVAngles(cout);
  PrintHAngles(cout);
  PrintGridAngles(cout);
  PrintBRDF_DB(cout);
  PrintBRDF_avg_val(cout);
  PrintBRDF_rgb_val(cout);

  app << flush;

  return;
}
//========================================================================
float
CDBDana_BRDFAnisotropic::GetAverageDiffuseness(void) const
{
  cout << "Warning: for Database BRDF it is impossile to define\n";
  cout << "average diffuseness as well as diffuseness itself\n";
  return 0;
}

float
CDBDana_BRDFAnisotropic::GetDiffuseness(const CHitPointInfo &hitinfo,
			       EDistributionType distType) const
{
  cout << "Warning: for Database BRDF it is impossile to define\n";
  cout << "diffuseness\n";
  return 0;
}

bool
CDBDana_BRDFAnisotropic::GetTransparency(CValues &transparency,
			       const CHitPointInfo &hitinfo) const
{
  cout << "Warning: for Database BRDF it is impossile to define\n";
  cout << "transparency\n";
  return false;
}


//========================================================================
bool
CDBDana_BRDFAnisotropic::GetDefaultColor(CValues &defaultColor,
			       const CHitPointInfo &hitInfo) const
{
  // there is the only color => it is the default one
  if (_db_color) {
    hitInfo.GetTextureValue(defaultColor, *_db_color, CMappingBase::EE_KD);
    return false;
  }

  defaultColor.Set(0.0);
  return true;
}
//========================================================================


// Should be checked !!!!!!!!!!!!!!!!!!!
//========================================================================
CSurface *
CDBDana_BRDFAnisotropic::Dup() const
{
  CDBDana_BRDFAnisotropic *ret = CreateNewObject();

  if (_db_color)
    ret->SetDB_color(_db_color);
  if (_brdf_i.size() > 0)
    ret->_brdf_i = _brdf_i;
  if(_db_s_tree)
    ret->_db_s_tree = _db_s_tree;
  if(_v_count)
    ret->_v_count = _v_count;
  if(_h_count)
    ret->_h_count = _h_count;
  if(_v_angles)
    ret->_v_angles = _v_angles;
  if(_h_angles)
    ret->_h_angles = _h_angles;
  if(_gridAngles.size() > 0)
    ret->_gridAngles = _gridAngles;
  if(_BRDF_avg_val.size() > 0)
    ret->_BRDF_avg_val = _BRDF_avg_val;
  if(_BRDF_rgb_val.size() > 0)
    ret->_BRDF_rgb_val = _BRDF_rgb_val;
  if(_hemigons.size() > 0)
    ret->_hemigons = _hemigons;

  //ret->_supportedType = _supportedType;
  return ret;
}
//========================================================================

// For a now is incorrect
//========================================================================
bool
CDBDana_BRDFAnisotropic::GetAlbedo(
  CValues &albedo,
  const CHitPointInfo &hitInfo,
  EDistributionType   distType) const
{
  CValues sum(CValues::EE_Weigth3); // result - zero by default
  CValues d(CValues::EE_Weigth3);
  CValues c(0,0,0);
  c.SetRepresentation(CValues::EE_Weigth3);
  bool retval = false;

  //assert(_search_tree != NULL);
  assert(_brdf_i.size() > 0);
  if(_db_s_tree == NULL)
  {
    cout << "Please, assign the search tree\n";
    return retval;
  }

  if((_v_angles == NULL) || (_h_angles == NULL))
  {
    cout << "Please, specify the horizontal and vertical angles";
    cout << "for a grid on the hemisphere.\n";
    return retval;
  }

  CVector3D N = hitInfo.GetNormal();
  if (hitInfo.GetDotProd() > 0.0)
    N = -N;
  CVector3D U = ArbitraryNormal(N);     // U is perpendicular to N
  CVector3D V = CrossProd(N, U);        // orthonormal base with N and U
  CVector3D dirInit = -((hitInfo.GetRay())->GetDir());

  //CValues kd_weighted(kd_local * ( 1.0 / maxRays));

  // number of rays created
  int num = 0;
  int samples_valid = 0;

  // Incoming direction is given in Cartesian coordinates. So, to compute
  // the query pointm they should be represented into the polar form.
  //   theta = acos(Z/R),      R = 1
  //   phi = arctan(Y/X)
  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  float incomingTheta = acos(dirInit.z);
  float incomingPhi = 0.0;
  float tann = 0.0;
  if(fabs(dirInit.x) < 1e-10)
    tann = MAXFLOAT;
  else
    tann = dirInit.y/dirInit.x;

  if(dirInit.x >=0 && dirInit.y >=0)
     incomingPhi = atan(tann);

  if((dirInit.x < 0 && dirInit.y >=0) || (dirInit.x <=0 && dirInit.y <= 0))
     incomingPhi = atan(tann) + M_PI;

  if(dirInit.x > 0 && dirInit.y < 0)
     incomingPhi = atan(tann) + 2.0*M_PI;


//#define _DEBUG_PHYSICAL_PLAUSIBILITY


#ifdef _DEBUG_PHYSICAL_PLAUSIBILITY
  // test if the albedo is smaller than 1.0 -- physically plausible
  if (c.Max() > 1.0) {
    if (_diffuse) {
      hitInfo.GetTextureValue(d, *_diffuse, CMappingBase::EE_KD);
      DEBUG << "diffuse = " << d << "\n";
    }
    if (_specular) {
      hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
      DEBUG << "specular = " << d << "\n";
    }
    if (_refract) {
      hitInfo.GetTextureValue(d, *_refract, CMappingBase::EE_KR);
      DEBUG << "refr = " << d << endl;
    }
  }
#endif // _DEBUG_PHYSICAL_PLAUSIBILITY

  if (retval) {
    albedo = c;
    return true;
  }

  // unable to compute the value for some reason
#ifdef _DEBUG_PHYSICAL_PLAUSIBILITY
  if ((!_diffuse) && (!_specular) && (!_refract)) {
    WARNING << "This is the black body - albedo cannot be computed" << endl;
  }
#endif // _DEBUG

  return false;
}

// Maximum albedo
//========================================================================
bool
CDBDana_BRDFAnisotropic::GetAlbedoMax(
  CValues &albedo,
  const CHitPointInfo &hitInfo,
  EDistributionType   distType) const
{
  //CValues c(CValues::EE_Weigth3); // result - zero by default
  //CValues d(CValues::EE_Weigth3);
  //bool retval = false;

  return false;
}

// Memory allocation
//========================================================================
void CDBDana_BRDFAnisotropic::AllocArrays()
{
  assert((_v_count > 0) && (_h_count > 0));

  if(_v_angles != NULL)
    delete [] _v_angles;
  _v_angles = new double[_v_count];

  if(_h_angles != NULL)
    delete [] _h_angles;
  _h_angles = new double[_h_count];

  if(_gridAngles.size() > 0)
    _gridAngles.clear();

  S_polarCoords pt;
  pt.theta = 0;    pt.phi = 0;
  _gridAngles = T_vec_polars(_v_count*_h_count, pt);

  CValues black(0,0,0);
  black.SetRepresentation(CValues::EE_Weigth3);
  if(_BRDF_avg_val.size() > 0)
    _BRDF_avg_val.clear();
  _BRDF_avg_val = T_vec_vec_d(_v_count*_h_count,
                              T_vec_d(_v_count*_h_count, 0));

  if(_BRDF_rgb_val.size() > 0)
    _BRDF_rgb_val.clear();
  _BRDF_rgb_val = T_vec_vec_cv(_v_count*_h_count,
                               T_vec_cv(_v_count*_h_count, black));

  if(_hemigons.size() > 0)
    _hemigons.clear();
  _hemigons = T_vec_hemigons(_v_count*_h_count, CHemisphereGonio());

  _validHemigons = false;
}


// Calculate vertical and horizontal angles by a given numbers
// v_count and h_count
//========================================================================
void CDBDana_BRDFAnisotropic::InitGridAngles()
{
  int thetaIndex, phiIndex;

  // set horizontal angles in uniform way from 0 to 2.0 * M_PI
  for (phiIndex = 0; phiIndex < _h_count; phiIndex++)
  {
    _h_angles[phiIndex] = phiIndex * 2.0 * M_PI / (double)(_h_count-1);
  }

  // set vertical angles in uniform way, from 0.0 to M_PI/2
  for (thetaIndex = 0; thetaIndex < _v_count; thetaIndex++)
  {
    _v_angles[thetaIndex] = thetaIndex*M_PI_2 / (double)(_v_count-1);
  }

  ComputeGrid();
  _validHemigons = false;

  return;
}

// Create a grid (get coordinates for each node of the grid and
// put them into the array). _v_angles, _h_angles should be initialized
//========================================================================
void CDBDana_BRDFAnisotropic::ComputeGrid()
{
  assert((_v_count > 0) && (_v_angles != NULL));
  assert((_h_count > 0) && (_h_angles != NULL));

  int thetaIndex, phiIndex;

  for(phiIndex = 0; phiIndex < _h_count; phiIndex++)
  {
    for (thetaIndex = 0; thetaIndex < _v_count; thetaIndex++)
    {
       int index = GetIndex(_v_count, thetaIndex, phiIndex);
       (_gridAngles[index]).theta = _v_angles[thetaIndex];
       (_gridAngles[index]).phi = _h_angles[phiIndex];
    }
  }
  _validHemigons = false;
  return;
}

// Read the measur_num ammount of data corresponding to the sample_num
// ID number of sample from the specified files and put it into the
// special data storage _brdf_i.
//========================================================================
void
CDBDana_BRDFAnisotropic::ReadDataDB(int &sample_num, int &measur_num,
                                    int &rdim, vector<const char*> &names_db)
{
  ReadData(sample_num, measur_num, rdim, names_db, _brdf_i);
  _validHemigons = false;
}

// Read the measur_num ammount of data corresponding to the sample_num
// ID number of sample from the specified files and put it into the
// special data storage _brdf_i and build the searching tree trough
// the 3D coordinates.
//========================================================================
void
CDBDana_BRDFAnisotropic::ReadDataBuildSTree(int &sample_num, int &measur_num,
                                    int &rdim, vector<const char*> &names_db)
{
  _db_s_tree = ReadDataBT(sample_num, measur_num, rdim, names_db, _brdf_i);
  _validHemigons = false;
}

// For a given incIndex direction compute the BRDF values (RGB/average) in
// the corresponding grid nodes by interpolation of raw data _brdf_i which
// were read from the Database (see ReadDataBuildSTree/ReadData).
//========================================================================
void
CDBDana_BRDFAnisotropic::GetBRDFforInc_i(int &rdim, double &incomingTheta,
                double &incomingPhi, int &incIndex, int &neigh, double &eps)
{
  assert(_db_s_tree != NULL);
  assert(_v_angles != NULL);
  assert(_h_angles != NULL);
  assert((neigh < (int)_brdf_i.size()) && (neigh > 0));
  assert((incIndex >= 0) && (incIndex < _v_count*_h_count));
  assert((eps >= 0) && (eps <= 1));

  // loop variables
  int thetaIndex, phiIndex;

  //int dim = 4;
  ANNpoint      query_pt;       // query point
  ANNidxArray   nn_idx;         // near neighbor indices
  ANNdistArray  dists;          // near neighbor distances
  query_pt = annAllocPt(rdim);   // allocate query point
  nn_idx = new ANNidx[neigh];   // allocate near neigh indices
  dists = new ANNdist[neigh];   // allocate near neighbor dists

  double *diagram = new double[_v_count*_h_count]; // Average BRDF values

  // Make interpolation between 'neigh' neighours for each node of
  // of the grid
  //----------------------------------------------------------------------
  for (thetaIndex = 0; thetaIndex < _v_count; thetaIndex++) {
    for (phiIndex = 0; phiIndex < _h_count; phiIndex++) {
      query_pt[0] = _v_angles[thetaIndex];   // thetta viewing
      query_pt[1] = _h_angles[phiIndex];     // phi viewing
      query_pt[2] = incomingTheta;           // thetta incoming
      if(rdim == 4)
        query_pt[3] = incomingPhi;           // thetta incoming

      //cout << "Query point: ";             // echo query point
      //printPt(cout, query_pt, dim);

      // search through the tree (raw data from the database) for 'neigh'
      // nearest to the query point neighbours.
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      _db_s_tree->annkSearch(query_pt, // query point
                             neigh,    // number of near neighbors
                             nn_idx,   // nearest neighbors (returned)
                             dists,    // distance (returned)
                             eps);     // error bound

      // We got i nearest neighbours for a query point. Now we can
      // get albedo for this query point by interpolation between
      // the albedos of these neighbours.
      //
      // Bilinear interpolation is unappropriate here, because the form
      // of the quadrilateral made of neighbours is unknown and it is
      // possible that query point lies outside this quadrilateral
      //
      // We now the distances between the neighbours and the query point in
      // 4D. Assume that the weight of the neighbour is proportiolnal to
      // to the distance between this neighbour and the query point. =>
      //   albedo_q = (alb1*w1 + alb2*w2 + ... + albn*wn)/(w1+w2+...+wn)
      // =>
      //   albedo_q = (alb1/dist1 + alb2/dist2 + ... + albn/distn)/
      //              (1/dist1 + 1/dist2 + ... + 1/distn)
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      //cout << "NN:\tIndex\tDistance\n";
      CValues sum(0,0,0);
      float sum_avg = 0;
      sum.SetRepresentation(CValues::EE_Weigth3);
      float total_weight = 0;
      int ind_coincide = -1;
      bool coincide = false;
      int r_idx_c = 0;

      // Get number of points in the tree. It should be equal three times
      // more than the length of vector _brdf_i
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ANNkdStats stats;
      _db_s_tree->getStats(stats);
      int repl_len = stats.n_pts/3;
      assert( repl_len == (int)_brdf_i.size());

      // Initial 4D coordinates (incoming/outgoing spherical coordinates)
      // in searching kdTree were read from the DB and then replicated
      // as shown in 'dbfunctions.h', functions ReadDataBT,
      // ReadData. In this way the number of points in the search
      // tree (n_pts) is three times more than the number of points in the
      // DB (N = _v_count*_h_count). However, the number of entries in the
      // vector _brdf_i (4D-coordinates with cooresponding BRDF values in
      // RGB format and average BRDF values) is N.
      // Therefore, to get an appropriate index for _brdf_i array we should
      // get index = found_index % N.
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      for (int i = 0; i < neigh; i++) {
        dists[i] = sqrt(dists[i]);      // unsquare distance
	int real_idx = nn_idx[i] % _brdf_i.size();
        // print summary
	/*cout << i << "\t" << real_idx << "\t" << dists[i] << "\tdata (";
        for(int yy = 0; yy < dim; yy++)
        {
           cout << ((_brdf_i[real_idx]).p)[yy] << ", ";
        }
        cout << ") -> (";
        for(int yy = 0; yy < dim; yy++)
        {
           cout << RadToDegree(((_brdf_i[real_idx]).p)[yy]) << ", ";
        }
        cout << ") brdf = " << ((_brdf_i[real_idx]).brdf) << ", ";
	cout << "avg = " << ((_brdf_i[real_idx]).avg) << "\n";*/

        // Sometimes it can happens that query point coincide with
	// one/more of the neighbours. In this case its BRDF value should
	// be equal to that of this neighbour.
	//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	if(dists[i] < 1e-3)
	{
	  coincide = true;
	  ind_coincide = i;
	  r_idx_c = nn_idx[ind_coincide] % _brdf_i.size();
	  //cout << "coincide!\n";
	}

	if(!coincide)
	{
	  sum += ((_brdf_i[real_idx]).brdf) * (1.0/dists[i]);
	  sum_avg += ((_brdf_i[real_idx]).avg) * (1.0/dists[i]);
	  total_weight += 1.0/dists[i];
	}
      }
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      if(!coincide)
      {
        sum *= 1.0/total_weight;
        sum_avg /= total_weight;
      }
      else
      {
        sum = (_brdf_i[r_idx_c]).brdf;
        sum_avg = (_brdf_i[r_idx_c]).avg;
      }

      // BRDF values should be in the range [0..1] !!!!
      if(sum_avg < 0) sum_avg = 0;
      if(sum_avg > 1) sum_avg = 1;
      sum.Clip(0, 1);
      //cout << "so, interpolated weight = " << sum;
      //cout << ", avg_sum = " << sum_avg << "\n";

      int index = GetIndex(_v_count, thetaIndex, phiIndex);
      _BRDF_rgb_val[incIndex][index] = sum;     // set RGB BRDF values
      _BRDF_avg_val[incIndex][index] = sum_avg; // set average BRDF values
      diagram[index] = sum_avg;
    } // for phiIndex
  } // for thetaIndex
  //----------------------------------------------------------------------

  // we have to set the intensity for the top pole the same
  //----------------------------------------------------------------------
  int tmp_ind = GetIndex(_v_count, 0, 0);
  CValues tmp_brdf = _BRDF_rgb_val[incIndex][tmp_ind];
  tmp_brdf.SetRepresentation(CValues::EE_Weigth3);
  double tmp_b_a = _BRDF_avg_val[incIndex][tmp_ind];

  for (phiIndex = 0; phiIndex < _h_count; phiIndex++)
  {
    int index = GetIndex(_v_count, 0, phiIndex);
    _BRDF_rgb_val[incIndex][index] = tmp_brdf;
    _BRDF_avg_val[incIndex][index] = tmp_b_a;
    diagram[index] = tmp_b_a;
  } // for phiIndex

  /*
  for (thetaIndex = 0; thetaIndex < _v_count; thetaIndex++) {
    for (phiIndex = 0; phiIndex < _h_count; phiIndex++) {
      cout << "diag[" << thetaIndex << "][" << phiIndex << "] = diag(";
      cout << GetIndex(_v_count, thetaIndex, phiIndex) << ") = ";
      cout << diagram[GetIndex(_v_count, thetaIndex, phiIndex)] << "\n";
    } // for phiIndex
  } // for thetaIndex
  */

  // For each incoming direction there is a CHemisphereGonio object
  // (goniometric web). The following is an initialization of the average
  // BRDF values in keypoints of goniometric web. It returns 'false' on
  // success, 'true' on failure.
  //----------------------------------------------------------------------
  _hemigons[incIndex].InitIntensities(_h_count, _h_angles, _v_count,
                                      _v_angles, diagram);

  _validHemigons = false;
  delete [] diagram;

  return;
}

// Compute the BRDF values for all the grids. The number of grids in
// isotropic case (this case) is equal to the number of vertical_angles
//========================================================================
void
CDBDana_BRDFAnisotropic::GetBRDFforAllInc(int &neigh, double &eps)
{
  assert(_db_s_tree != NULL);
  assert((_v_angles != NULL) && (_v_count > 0));
  assert((_h_angles != NULL) && (_h_count > 0));
  assert((int)_BRDF_avg_val.size() == _v_count*_h_count);
  assert((int)_BRDF_rgb_val.size() == _v_count*_h_count);
  assert((int)_hemigons.size() == _v_count*_h_count);

  int thetaIndex, phiIndex;
  int dim = 4;

  for(phiIndex = 0; phiIndex < _h_count; phiIndex++)
  {
    //double incPhi = _h_angles[phiIndex];
    for (thetaIndex = 0; thetaIndex < _v_count; thetaIndex++)
    {
       int index = GetIndex(_v_count, thetaIndex, phiIndex);
       //double incTheta = _v_angles[thetaIndex];
       //GetBRDFforInc_i(incTheta, incPhi, thetaIndex, neigh, eps);
       GetBRDFforInc_i(dim, _v_angles[thetaIndex], _h_angles[phiIndex],
                       index, neigh, eps);
    //_hemigons[thetaIndex].PrintAll();
    }
  }
  _validHemigons = true;

  return;
}


// Data processing from the scratch. No precomputations are provided.
// Only the minimum number of parameters are given.
// Function provides: vertical, horizontal angles, grid coordinates,
//   reading and interpolation of the BRDF data for each incoming direction
//   (each grid), precomputations for further mapping (InitIntensities
//   for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFAnisotropic::InitData(int &sample_num, int &measur_num,
                       vector<const char*> &names_db, int &neigh, double &eps,
		       int &v_count, int &h_count )
{
  _v_count = v_count;  _h_count = h_count;
  int dim = 4;  // dimension of the tree

  AllocArrays();
  InitGridAngles();
  ReadDataBuildSTree(sample_num, measur_num, dim, names_db);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _h_count*_v_count) _validHemigons = true;
  else _validHemigons = false;
}

// Data processing: vertical and horizontal angles are precomputed.
// Function provides: grid coordinates, reading and interpolation of the BRDF
//   data for each incoming direction (each grid), precomputations for further
//   mapping (InitIntensities for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFAnisotropic::InitData(int &sample_num, int &measur_num,
                       vector<const char*> &names_db, int &neigh, double &eps,
		       int &v_count, int &h_count,
		       double *v_angles, double *h_angles)
{
  _v_count = v_count;  _h_count = h_count;
  int dim = 4;  // dimension of the tree

  AllocArrays();
  SetVHAngles(v_count, v_angles, h_count, h_angles);
  ComputeGrid();
  ReadDataBuildSTree(sample_num, measur_num, dim, names_db);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _h_count*_v_count) _validHemigons = true;
  else _validHemigons = false;
}

// Data processing: vertical, horizontal angles, and grid coordinates
//   are precomputed.
// Function provides: reading and interpolation of the BRDF data for each
//   incoming direction (each grid), precomputations for further mapping
//   (InitIntensities for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFAnisotropic::InitData(int &sample_num, int &measur_num,
                       vector<const char*> &names_db,
                       int &neigh, double &eps, int &v_count, int &h_count,
                       double *v_angles, double *h_angles,
                       T_vec_polars &grid_angles)
{
  _v_count = v_count;  _h_count = h_count;
  int dim = 4;  // dimension of the tree

  AllocArrays();
  SetGrid(v_count, v_angles, h_count, h_angles, grid_angles);
  ReadDataBuildSTree(sample_num, measur_num, dim, names_db);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _h_count*_v_count) _validHemigons = true;
  else _validHemigons = false;
}

// Data processing: search tree for coordinates is built.
// Function provides: vertical, horizontal angles, grid coordinates, reading
//   and interpolation of the BRDF data for each incoming direction (each
//   grid), precomputations for further mapping (InitIntensities for
//   goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFAnisotropic::InitData(int &sample_num, int &measur_num,
                       vector<const char*> &names_db, int &neigh, double &eps,
		       int &v_count, int &h_count, ANNkd_tree *tree)
{
  _v_count = v_count;  _h_count = h_count;
  assert(tree != NULL);
  int dim = 4;  // dimension of the tree

  AllocArrays();
  InitGridAngles();
  ReadDataDB(sample_num, measur_num, dim, names_db);

  // the number of points after the replication
  int num_pts = measur_num*3;

  // Has the 'tree' an appropriate number of points and dimension.
  ANNkdStats st;
  tree->getStats(st);

  if((st.dim != dim) || (st.n_pts != num_pts))
  {
    cout << "Please, set an appropriate search tree. It's n_pts should\n";
    cout << "be equal to the measur_num*3, and dim should be equal to 4\n";
    exit(3);
  }

  SetDBSearchTree(num_pts, dim, tree);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _h_count*_v_count) _validHemigons = true;
  else _validHemigons = false;
}


// Data processing: data are taken from the DB, search tree for coordinates
//   is also built.
// Function provides: vertical, horizontal angles, grid coordinates,
//   interpolation of the BRDF data for each incoming direction (each grid),
//   precomputations for further mapping (InitIntensities for goniometric
//   web) for each grid
//========================================================================
void
CDBDana_BRDFAnisotropic::InitData(int &neigh, double &eps, int &v_count,
                         int &h_count, T_dataDB &brdf_i, ANNkd_tree *tree)
{
  _v_count = v_count;  _h_count = h_count;
  assert(tree != NULL);
  int dim = 4;  // dimension of the tree

  if((int)brdf_i.size() <= 0)
  {
    cout << "Please, set an appropriate brdf_i vector. It's size should\n";
    cout << "be equal to the measur_num.\n";
    exit(3);
  }

  AllocArrays();
  InitGridAngles();

  // the number of points after the replication
  int num_pts = (int)brdf_i.size()*3;

  // Has the 'tree' an appropriate number of points and dimension.
  ANNkdStats st;
  tree->getStats(st);

  if((st.dim != dim) || (st.n_pts != num_pts))
  {
    cout << "Please, set an appropriate search tree. It's n_pts should\n";
    cout << "be equal to the measur_num*3, and dim should be equal to 4\n";
    exit(3);
  }
  SetAllDBData(brdf_i, num_pts, dim, tree);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _h_count*_v_count) _validHemigons = true;
  else _validHemigons = false;
}



// Data processing: vertical and horizontal angles are precomputed, search
//   tree for coordinates is built.
// Function provides: grid coordinates, reading and interpolation of the
//   BRDF data for each incoming direction (each grid), precomputations
//   for further mapping (InitIntensities for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFAnisotropic::InitData(int &sample_num, int &measur_num,
                       vector<const char*> &names_db, int &neigh, double &eps,
		       int &v_count, int &h_count, double *v_angles,
		       double *h_angles, ANNkd_tree *tree)
{
  _v_count = v_count;  _h_count = h_count;
  assert(tree != NULL);
  int dim = 4;  // dimension of the tree

  AllocArrays();
  SetVHAngles(v_count, v_angles, h_count, h_angles);
  ComputeGrid();
  ReadDataDB(sample_num, measur_num, dim, names_db);

  // the number of points after the replication
  int num_pts = measur_num*3;

  // Has the 'tree' an appropriate number of points and dimension.
  ANNkdStats st;
  tree->getStats(st);

  if((st.dim != dim) || (st.n_pts != num_pts))
  {
    cout << "Please, set an appropriate search tree. It's n_pts should\n";
    cout << "be equal to the measur_num*3, and dim should be equal to 4\n";
    exit(3);
  }

  SetDBSearchTree(num_pts, dim, tree);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _h_count*_v_count) _validHemigons = true;
  else _validHemigons = false;
}


// Data processing: vertical and horizontal angles are precomputed, data are
//   taken from the DB, search tree for coordinates is also built,
// Function provides: grid coordinates, interpolation of the BRDF data for
//   each incoming direction (each grid), precomputations for further mapping
//   (InitIntensities for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFAnisotropic::InitData(int &neigh, double &eps, int &v_count,
		       int &h_count, double *v_angles, double *h_angles,
                       T_dataDB &brdf_i, ANNkd_tree *tree)
{
  _v_count = v_count;  _h_count = h_count;
  assert(tree != NULL);

  if((int)brdf_i.size() <= 0)
  {
    cout << "Please, set an appropriate brdf_i vector. It's size should\n";
    cout << "be equal to the measur_num.\n";
    exit(3);
  }

  AllocArrays();
  SetVHAngles(v_count, v_angles, h_count, h_angles);
  ComputeGrid();

  // the number of points after the replication
  int num_pts = (int)brdf_i.size()*3;
  // dimension of the tree
  int dim = 4;

  // Has the 'tree' an appropriate number of points and dimension.
  ANNkdStats st;
  tree->getStats(st);

  if((st.dim != dim) || (st.n_pts != num_pts))
  {
    cout << "Please, set an appropriate search tree. It's n_pts should\n";
    cout << "be equal to the measur_num*3, and dim should be equal to 4\n";
    exit(3);
  }
  SetAllDBData(brdf_i, num_pts, dim, tree);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _h_count*_v_count) _validHemigons = true;
  else _validHemigons = false;
}


// Data processing: vertical, horizontal angles, and grid coordinates are
//   precomputed, data are taken from the DB, search tree for coordinates is
//   also built,
// Function provides: interpolation of the BRDF data for each incoming
//   direction (each grid), precomputations for further mapping
//   (InitIntensities for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFAnisotropic::InitData(int &neigh, double &eps, int &v_count,
                       int &h_count, double *v_angles, double *h_angles,
                       T_vec_polars &grid_angles, T_dataDB &brdf_i,
                       ANNkd_tree *tree)
{
  _v_count = v_count;  _h_count = h_count;
  assert(tree != NULL);

  if((int)brdf_i.size() <= 0)
  {
    cout << "Please, set an appropriate brdf_i vector. It's size should\n";
    cout << "be equal to the measur_num.\n";
    exit(3);
  }
  AllocArrays();
  SetGrid(v_count, v_angles, h_count, h_angles, grid_angles);

  // the number of points after the replication
  int num_pts = (int)brdf_i.size()*3;
  // dimension of the tree
  int dim = 4;

  // Has the 'tree' an appropriate number of points and dimension.
  ANNkdStats st;
  tree->getStats(st);

  if((st.dim != dim) || (st.n_pts != num_pts))
  {
    cout << "Please, set an appropriate search tree. It's n_pts should\n";
    cout << "be equal to the measur_num*3, and dim should be equal to 4\n";
    exit(3);
  }
  SetAllDBData(brdf_i, num_pts, dim, tree);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _h_count*_v_count) _validHemigons = true;
  else _validHemigons = false;
}


//========================================================================
int
CDBDana_BRDFAnisotropic::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  return maxRays;
}
//========================================================================


//========================================================================
int
CDBDana_BRDFAnisotropic::SpawnStochasticRays(
  CContributionArray   &carray,
  const CHitPointInfo  &hitInfo,
  CSamplingInfo        &samplingInfo,
  EDistributionType     distType) const
{
  return 1;
}
//========================================================================


//========================================================================
bool
CDBDana_BRDFAnisotropic::EvaluateBSDF(CValues &bsdf,
		           const CHitPointInfo &hitInfo,
			   const CVector3D &dirOut,
			   EDistributionType distType) const
{
  return false;
}
//========================================================================


//************************************************************************
// Functions to print computed/set class members
//========================================================================

//========================================================================
void CDBDana_BRDFAnisotropic::PrintDB_SearchTree(ostream &app) const
{
  if(_db_s_tree == NULL)
  {
    app << "Data is uninithalized! Please, apply 'InitData' method\n";
    app << "and then try again\n";
    return;
  }
  app << "Search tree for 4D coordinates (outgoing/incoming theta and\n";
  app << "phi angles which where read from the database):\n";
  app << "------------------------------------------------------------\n";
  _db_s_tree->Print(ANNtrue, app);
}

//========================================================================
void CDBDana_BRDFAnisotropic::PrintVAngles(ostream &app) const
{
  if(_v_count > 0)
  {
    app << "Vertical angles (specified " << _v_count << " angles):\n";
    app << "------------------------------------------------------------\n";

    assert(_v_angles != NULL);

    for(int i = 0; i < _v_count; i++)
    {
      app << "a[" << i << "]: RAD(" << _v_angles[i] << "), ";
      app << "DEGREES(" << _v_angles[i]*180.0/M_PI << ")\n";
    }
  }
  return;
}

//========================================================================
void CDBDana_BRDFAnisotropic::PrintHAngles(ostream &app) const
{
  if(_h_count > 0)
  {
    app << "Horizontal angles (specified " << _h_count << " angles):\n";
    app << "------------------------------------------------------------\n";

    assert(_h_angles != NULL);

    for(int i = 0; i < _h_count; i++)
    {
      app << "a[" << i << "]: RAD(" << _h_angles[i] << "), ";
      app << "DEGREES(" << _h_angles[i]*180.0/M_PI << ")\n";
    }
  }
  return;
}


//========================================================================
void CDBDana_BRDFAnisotropic::PrintGridAngles(ostream &app) const
{
  if((_h_count > 0) && (_v_count > 0))
  {
    assert((int)_gridAngles.size() == _v_count*_h_count);
    app << "Grid coordinates (size " << _v_count << "x";
    app << _h_count << "):\n";
    app << "------------------------------------------------------------\n";
    for(int i = 0; i < (int)_gridAngles.size(); i++)
    {
      app << "[" << i << "]" << _gridAngles[i] << ", ";
      if(!((i+1) % _v_count)) app << "\n";
    }
    app << "\n";
  }
  return;
}


//========================================================================
void CDBDana_BRDFAnisotropic::PrintBRDF_DB(ostream &app) const
{
  if((int)_brdf_i.size() > 0 && _validHemigons)
  {
    app << "Polar coordinates (thetta_out, phi_out, thetta_inc, phi_inc)\n";
    app << "and corresponding BRDF values in RGB format and average values\n";
    app << "read from the database\n";
    app << "------------------------------------------------------------\n";
    for(int i = 0; i < (int)_brdf_i.size(); i++)
    {
      app << "[" << i << "]: " << _brdf_i[i] << "\n";
    }
  }

  return;
}

//========================================================================
void CDBDana_BRDFAnisotropic::PrintBRDF_avg_val(ostream &app) const
{
  if((int)_BRDF_avg_val.size() > 0 && _validHemigons)
  {
    assert((int)_BRDF_avg_val.size() == _v_count*_h_count);
    app << "BRDF average values (isotropic case). The number of incoming \n";
    app << "light directions is equal to " << _BRDF_avg_val.size();
    app << "(the number of vertical\n";
    app << "angles multiplied the number of horizontal_angles)\n";
    app << "------------------------------------------------------------\n";

    for(int i = 0; i < (int)_BRDF_avg_val.size(); i++)
    {
      app << "Grid for [" << i << "]->" << _gridAngles[i];
      app << " incoming direction:\n";
      app << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
      assert((int)_BRDF_avg_val[i].size() == _h_count*_v_count);
      for(int j = 0; j < (int)_BRDF_avg_val[i].size(); j++)
      {
        app << "[" << j << "] " << _BRDF_avg_val[i][j] << "\t";
	if(!((j+1) % _v_count)) app << "\n";
      }
      app << "\n";
    }
    return;
  }
  app << "Requested _BRDF_avg_val array is empty\n";

  return;
}


//========================================================================
void CDBDana_BRDFAnisotropic::PrintBRDF_rgb_val(ostream &app) const
{
  if((int)_BRDF_rgb_val.size() > 0 && _validHemigons)
  {
    assert((int)_BRDF_avg_val.size() == _v_count*_h_count);
    app << "BRDF RGB values (isotropic case). The number of incoming \n";
    app << "light directions is equal to " << _BRDF_avg_val.size();
    app << " (the number of vertical\n";
    app << "angles miltiplied by the number of horizontal_angles)\n";
    app << "------------------------------------------------------------\n";

    for(int i = 0; i < (int)_BRDF_rgb_val.size(); i++)
    {
      app << "Grid for [" << i << "]->"<< _gridAngles[i];
      app << " incoming direction:\n";
      app << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
      assert((int)_BRDF_rgb_val[i].size() == _h_count*_v_count);
      for(int j = 0; j < (int)_BRDF_rgb_val[i].size(); j++)
      {
        app << "[" << j << "] " << _BRDF_rgb_val[i][j] << "\n";
	if(!((j+1) % _v_count)) app << "\n";
      }
      app << "\n";
    }
    return;
 }
 app << "Requested _BRDF_rgb_val array is empty\n";

 return;
}





/**********************************************************************
 *                                                                    *
 *                    CDBDana_BRDFIsotropic methods                   *
 *                                                                    *
 **********************************************************************/
// describes the surface to a given stream
void
CDBDana_BRDFIsotropic::Describe(ostream &app, int ind) const
{
  app << "\n";
  app << "*********************************************";
  app << "*************************\n";
  indent(app, ind);
  app << "CDBDana_BRDFIsotropic:\n";
  app << "*********************************************";
  app << "*************************\n\n";

  if (_db_color) {
    indent(app, ind+2);
    app << "  db_color = \t";
    _db_color->Describe(app, ind); app << "\n";
  }
  PrintDB_SearchTree(cout);
  PrintVAngles(cout);
  PrintHAngles(cout);
  PrintGridAngles(cout);
  PrintBRDF_DB(cout);
  PrintBRDF_avg_val(cout);
  PrintBRDF_rgb_val(cout);

  app << flush;

  return;
}
//========================================================================
float
CDBDana_BRDFIsotropic::GetAverageDiffuseness(void) const
{
  cout << "Warning: for Database BRDF it is impossile to define\n";
  cout << "average diffuseness as well as diffuseness itself\n";
  return 0;
}

float
CDBDana_BRDFIsotropic::GetDiffuseness(const CHitPointInfo &hitinfo,
			       EDistributionType distType) const
{
  cout << "Warning: for Database BRDF it is impossile to define\n";
  cout << "diffuseness\n";
  return 0;
}

bool
CDBDana_BRDFIsotropic::GetTransparency(CValues &transparency,
			       const CHitPointInfo &hitinfo) const
{
  cout << "Warning: for Database BRDF it is impossile to define\n";
  cout << "transparency\n";
  return false;
}

//========================================================================
bool
CDBDana_BRDFIsotropic::GetDefaultColor(CValues &defaultColor,
			       const CHitPointInfo &hitInfo) const
{
  // there is the only color => it is the default one
  if (_db_color) {
    hitInfo.GetTextureValue(defaultColor, *_db_color, CMappingBase::EE_KD);
    return false;
  }

  defaultColor.Set(0.0);
  return true;
}

// Should be checked !!!!!!!!!!!!!!!!!!!
//========================================================================
CSurface *
CDBDana_BRDFIsotropic::Dup() const
{
  CDBDana_BRDFIsotropic *ret = CreateNewObject();

  if (_db_color)
    ret->SetDB_color(_db_color);
  if (_brdf_i.size() > 0)
    ret->_brdf_i = _brdf_i;
  if(_db_s_tree)
    ret->_db_s_tree = _db_s_tree;
  if(_v_count)
    ret->_v_count = _v_count;
  if(_h_count)
    ret->_h_count = _h_count;
  if(_v_angles)
    ret->_v_angles = _v_angles;
  if(_h_angles)
    ret->_h_angles = _h_angles;
  if(_gridAngles.size() > 0)
    ret->_gridAngles = _gridAngles;
  if(_BRDF_avg_val.size() > 0)
    ret->_BRDF_avg_val = _BRDF_avg_val;
  if(_BRDF_rgb_val.size() > 0)
    ret->_BRDF_rgb_val = _BRDF_rgb_val;
  if(_hemigons.size() > 0)
    ret->_hemigons = _hemigons;

  //ret->_supportedType = _supportedType;
  return ret;
}
//========================================================================

// For a now is incorrect
//========================================================================
bool
CDBDana_BRDFIsotropic::GetAlbedo(
  CValues &albedo,
  const CHitPointInfo &hitInfo,
  EDistributionType   distType) const
{
  CValues sum(CValues::EE_Weigth3); // result - zero by default
  CValues d(CValues::EE_Weigth3);
  CValues c(0,0,0);
  c.SetRepresentation(CValues::EE_Weigth3);
  bool retval = false;

  //assert(_search_tree != NULL);
  assert(_brdf_i.size() > 0);
  if(_db_s_tree == NULL)
  {
    cout << "Please, assign the search tree\n";
    return retval;
  }

  if((_v_angles == NULL) || (_h_angles == NULL))
  {
    cout << "Please, specify the horizontal and vertical angles";
    cout << "for a grid on the hemisphere.\n";
    return retval;
  }

  CVector3D N = hitInfo.GetNormal();
  if (hitInfo.GetDotProd() > 0.0)
    N = -N;
  CVector3D U = ArbitraryNormal(N);     // U is perpendicular to N
  CVector3D V = CrossProd(N, U);        // orthonormal base with N and U
  CVector3D dirInit = -((hitInfo.GetRay())->GetDir());

  //CValues kd_weighted(kd_local * ( 1.0 / maxRays));

  // number of rays created
  int num = 0;
  int samples_valid = 0;

  // Incoming direction is given in Cartesian coordinates. So, to compute
  // the query pointm they should be represented into the polar form.
  //   theta = acos(Z/R),      R = 1
  //   phi = arctan(Y/X)
  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  float incomingTheta = acos(dirInit.z);
  float incomingPhi = 0.0;
  float tann = 0.0;
  if(fabs(dirInit.x) < 1e-10)
    tann = MAXFLOAT;
  else
    tann = dirInit.y/dirInit.x;

  if(dirInit.x >=0 && dirInit.y >=0)
     incomingPhi = atan(tann);

  if((dirInit.x < 0 && dirInit.y >=0) || (dirInit.x <=0 && dirInit.y <= 0))
     incomingPhi = atan(tann) + M_PI;

  if(dirInit.x > 0 && dirInit.y < 0)
     incomingPhi = atan(tann) + 2.0*M_PI;


//#define _DEBUG_PHYSICAL_PLAUSIBILITY


#ifdef _DEBUG_PHYSICAL_PLAUSIBILITY
  // test if the albedo is smaller than 1.0 -- physically plausible
  if (c.Max() > 1.0) {
    if (_diffuse) {
      hitInfo.GetTextureValue(d, *_diffuse, CMappingBase::EE_KD);
      DEBUG << "diffuse = " << d << "\n";
    }
    if (_specular) {
      hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
      DEBUG << "specular = " << d << "\n";
    }
    if (_refract) {
      hitInfo.GetTextureValue(d, *_refract, CMappingBase::EE_KR);
      DEBUG << "refr = " << d << endl;
    }
  }
#endif // _DEBUG_PHYSICAL_PLAUSIBILITY

  if (retval) {
    albedo = c;
    return true;
  }

  // unable to compute the value for some reason
#ifdef _DEBUG_PHYSICAL_PLAUSIBILITY
  if ((!_diffuse) && (!_specular) && (!_refract)) {
    WARNING << "This is the black body - albedo cannot be computed" << endl;
  }
#endif // _DEBUG

  return false;
}

// Memory allocation
//========================================================================
void CDBDana_BRDFIsotropic::AllocArrays()
{
  if(_v_angles != NULL)
    delete [] _v_angles;
  _v_angles = new double[_v_count];

  if(_h_angles != NULL)
    delete [] _h_angles;
  _h_angles = new double[_h_count];

  if(_gridAngles.size() > 0)
    _gridAngles.clear();

  S_polarCoords pt;
  pt.theta = 0;    pt.phi = 0;
  _gridAngles = T_vec_polars(_v_count*_h_count, pt);

  CValues black(0,0,0);
  black.SetRepresentation(CValues::EE_Weigth3);
  if(_BRDF_avg_val.size() > 0)
    _BRDF_avg_val.clear();
  _BRDF_avg_val = T_vec_vec_d(_v_count, T_vec_d(_v_count*_h_count, 0));

  if(_BRDF_rgb_val.size() > 0)
    _BRDF_rgb_val.clear();
  _BRDF_rgb_val = T_vec_vec_cv(_v_count, T_vec_cv(_v_count*_h_count, black));

  if(_hemigons.size() > 0)
    _hemigons.clear();
  _hemigons = T_vec_hemigons(_v_count, CHemisphereGonio());

  _validHemigons = false;
}

// Read the measur_num ammount of data corresponding to the sample_num
// ID number of sample from the specified files and put it into the
// special data storage _brdf_i.
//========================================================================
/*void
CDBDana_BRDFIsotropic::ReadDataDB(int &sample_num, int &measur_num,
                          int &rdim, vector<const char*> &names_db)
{
  ReadData(sample_num, measur_num, rdim, names_db, _brdf_i);
  _validHemigons = false;
}

// Read the measur_num ammount of data corresponding to the sample_num
// ID number of sample from the specified files and put it into the
// special data storage _brdf_i and build the searching tree trough
// the 3D coordinates.
//========================================================================
void
CDBDana_BRDFIsotropic::ReadDataBuildSTree(int &sample_num, int &measur_num,
                                  int &rdim, vector<const char*> &names_db)
{
  _db_s_tree = ReadDataBT(sample_num, measur_num, rdim, names_db, _brdf_i);
  _validHemigons = false;
}*/

// Compute the BRDF values for all the grids. The number of grids in
// isotropic case (this case) is equal to the number of vertical_angles
//========================================================================
void
CDBDana_BRDFIsotropic::GetBRDFforAllInc(int &neigh, double &eps)
{
  assert(_db_s_tree != NULL);
  assert(_v_angles != NULL);
  assert(_h_angles != NULL);

  double incPhi = 0;
  int thetaIndex;
  int dim = 3;

  for (thetaIndex = 0; thetaIndex < _v_count; thetaIndex++)
  {
    //CHemisphereGonio hg;
    double incTheta = _v_angles[thetaIndex];
    GetBRDFforInc_i(dim, incTheta, incPhi, thetaIndex, neigh, eps);
    //_hemigons[thetaIndex].PrintAll();
  }
  _validHemigons = true;

  return;
}

// Data processing from the scratch. No precomputations are provided.
// Only the minimum number of parameters are given.
// Function provides: vertical, horizontal angles, grid coordinates,
//   reading and interpolation of the BRDF data for each incoming direction
//   (each grid), precomputations for further mapping (InitIntensities
//   for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFIsotropic::InitData(int &sample_num, int &measur_num,
                       vector<const char*> &names_db, int &neigh, double &eps,
		       int &v_count, int &h_count )
{
  _v_count = v_count;  _h_count = h_count;
  int dim = 3;

  AllocArrays();
  InitGridAngles();
  ReadDataBuildSTree(sample_num, measur_num, dim, names_db);
  GetBRDFforAllInc(neigh, eps);
  cout << "BRDF_rgb_val[0].size = " << _BRDF_rgb_val[0].size() << "\n";
  if((int)_hemigons.size() == _v_count) _validHemigons = true;
  else _validHemigons = false;
}



// Data processing: vertical and horizontal angles are precomputed.
// Function provides: grid coordinates, reading and interpolation of the BRDF
//   data for each incoming direction (each grid), precomputations for further
//   mapping (InitIntensities for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFIsotropic::InitData(int &sample_num, int &measur_num,
                       vector<const char*> &names_db, int &neigh, double &eps,
		       int &v_count, int &h_count,
		       double *v_angles, double *h_angles)
{
  _v_count = v_count;  _h_count = h_count;
  int dim = 3;

  AllocArrays();
  SetVHAngles(v_count, v_angles, h_count, h_angles);
  ComputeGrid();
  ReadDataBuildSTree(sample_num, measur_num, dim, names_db);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _v_count) _validHemigons = true;
  else _validHemigons = false;
}

// Data processing: vertical, horizontal angles, and grid coordinates
//   are precomputed.
// Function provides: reading and interpolation of the BRDF data for each
//   incoming direction (each grid), precomputations for further mapping
//   (InitIntensities for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFIsotropic::InitData(int &sample_num, int &measur_num,
                       vector<const char*> &names_db,
                       int &neigh, double &eps, int &v_count, int &h_count,
                       double *v_angles, double *h_angles,
                       T_vec_polars &grid_angles)
{
  _v_count = v_count;  _h_count = h_count;
  int dim = 3;

  AllocArrays();
  SetGrid(v_count, v_angles, h_count, h_angles, grid_angles);
  ReadDataBuildSTree(sample_num, measur_num, dim, names_db);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _v_count) _validHemigons = true;
  else _validHemigons = false;
}


// Data processing: search tree for coordinates is built.
// Function provides: vertical, horizontal angles, grid coordinates, reading
//   and interpolation of the BRDF data for each incoming direction (each
//   grid), precomputations for further mapping (InitIntensities for
//   goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFIsotropic::InitData(int &sample_num, int &measur_num,
                       vector<const char*> &names_db, int &neigh, double &eps,
		       int &v_count, int &h_count, ANNkd_tree *tree)
{
  _v_count = v_count;  _h_count = h_count;
  assert(tree != NULL);
  int dim = 3;  // dimension of the tree

  AllocArrays();
  InitGridAngles();
  ReadDataDB(sample_num, measur_num, dim, names_db);

  // the number of points after the replication
  int num_pts = measur_num*3;

  // Has the 'tree' an appropriate number of points and dimension.
  ANNkdStats st;
  tree->getStats(st);

  if((st.dim != dim) || (st.n_pts != num_pts))
  {
    cout << "Please, set an appropriate search tree. It's n_pts should\n";
    cout << "be equal to the measur_num*3, and dim should be equal to 3\n";
    exit(3);
  }

  SetDBSearchTree(num_pts, dim, tree);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _v_count) _validHemigons = true;
  else _validHemigons = false;
}


// Data processing: data are taken from the DB, search tree for coordinates
//   is also built.
// Function provides: vertical, horizontal angles, grid coordinates,
//   interpolation of the BRDF data for each incoming direction (each grid),
//   precomputations for further mapping (InitIntensities for goniometric
//   web) for each grid
//========================================================================
void
CDBDana_BRDFIsotropic::InitData(int &neigh, double &eps, int &v_count,
                         int &h_count, T_dataDB &brdf_i, ANNkd_tree *tree)
{
  _v_count = v_count;  _h_count = h_count;
  assert(tree != NULL);

  if((int)brdf_i.size() <= 0)
  {
    cout << "Please, set an appropriate brdf_i vector. It's size should\n";
    cout << "be equal to the measur_num.\n";
    exit(3);
  }
  AllocArrays();
  InitGridAngles();

  // the number of points after the replication
  int num_pts = (int)brdf_i.size()*3;
  // dimension of the tree
  int dim = 3;

  // Has the 'tree' an appropriate number of points and dimension.
  ANNkdStats st;
  tree->getStats(st);

  if((st.dim != dim) || (st.n_pts != num_pts))
  {
    cout << "Please, set an appropriate search tree. It's n_pts should\n";
    cout << "be equal to the measur_num*3, and dim should be equal to 3\n";
    exit(3);
  }
  SetAllDBData(brdf_i, num_pts, dim, tree);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _v_count) _validHemigons = true;
  else _validHemigons = false;
}



// Data processing: vertical and horizontal angles are precomputed, search
//   tree for coordinates is built.
// Function provides: grid coordinates, reading and interpolation of the
//   BRDF data for each incoming direction (each grid), precomputations
//   for further mapping (InitIntensities for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFIsotropic::InitData(int &sample_num, int &measur_num,
                       vector<const char*> &names_db, int &neigh, double &eps,
		       int &v_count, int &h_count, double *v_angles,
		       double *h_angles, ANNkd_tree *tree)
{
  _v_count = v_count;  _h_count = h_count;
  int dim = 3;  // dimension of the tree
  assert(tree != NULL);

  AllocArrays();
  SetVHAngles(v_count, v_angles, h_count, h_angles);
  ComputeGrid();
  ReadDataDB(sample_num, measur_num, dim, names_db);

  // the number of points after the replication
  int num_pts = measur_num*3;

  // Has the 'tree' an appropriate number of points and dimension.
  ANNkdStats st;
  tree->getStats(st);

  if((st.dim != dim) || (st.n_pts != num_pts))
  {
    cout << "Please, set an appropriate search tree. It's n_pts should\n";
    cout << "be equal to the measur_num*3, and dim should be equal to 3\n";
    exit(3);
  }

  SetDBSearchTree(num_pts, dim, tree);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _v_count) _validHemigons = true;
  else _validHemigons = false;
}


// Data processing: vertical and horizontal angles are precomputed, data are
//   taken from the DB, search tree for coordinates is also built,
// Function provides: grid coordinates, interpolation of the BRDF data for
//   each incoming direction (each grid), precomputations for further mapping
//   (InitIntensities for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFIsotropic::InitData(int &neigh, double &eps, int &v_count,
                       int &h_count, double *v_angles, double *h_angles,
                       T_dataDB &brdf_i, ANNkd_tree *tree)
{
  _v_count = v_count;  _h_count = h_count;
  assert(tree != NULL);

  if((int)brdf_i.size() <= 0)
  {
    cout << "Please, set an appropriate brdf_i vector. It's size should\n";
    cout << "be equal to the measur_num.\n";
    exit(3);
  }
  AllocArrays();
  SetVHAngles(v_count, v_angles, h_count, h_angles);
  ComputeGrid();

  // the number of points after the replication
  int num_pts = (int)brdf_i.size()*3;
  // dimension of the tree
  int dim = 3;

  // Has the 'tree' an appropriate number of points and dimension.
  ANNkdStats st;
  tree->getStats(st);

  if((st.dim != dim) || (st.n_pts != num_pts))
  {
    cout << "Please, set an appropriate search tree. It's n_pts should\n";
    cout << "be equal to the measur_num*3, and dim should be equal to 3\n";
    exit(3);
  }
  SetAllDBData(brdf_i, num_pts, dim, tree);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _v_count) _validHemigons = true;
  else _validHemigons = false;
}


// Data processing: vertical, horizontal angles, and grid coordinates are
//   precomputed, data are taken from the DB, search tree for coordinates is
//   also built,
// Function provides: interpolation of the BRDF data for each incoming
//   direction (each grid), precomputations for further mapping
//   (InitIntensities for goniometric web) for each grid
//========================================================================
void
CDBDana_BRDFIsotropic::InitData(int &neigh, double &eps, int &v_count,
                       int &h_count, double *v_angles, double *h_angles,
                       T_vec_polars &grid_angles, T_dataDB &brdf_i,
                       ANNkd_tree *tree)
{
  _v_count = v_count;  _h_count = h_count;
  assert(tree != NULL);

  if((int)brdf_i.size() <= 0)
  {
    cout << "Please, set an appropriate brdf_i vector. It's size should\n";
    cout << "be equal to the measur_num.\n";
    exit(3);
  }
  AllocArrays();
  SetGrid(v_count, v_angles, h_count, h_angles, grid_angles);

  // the number of points after the replication
  int num_pts = (int)brdf_i.size()*3;
  // dimension of the tree
  int dim = 3;

  // Has the 'tree' an appropriate number of points and dimension.
  ANNkdStats st;
  tree->getStats(st);

  if((st.dim != dim) || (st.n_pts != num_pts))
  {
    cout << "Please, set an appropriate search tree. It's n_pts should\n";
    cout << "be equal to the measur_num*3, and dim should be equal to 3\n";
    exit(3);
  }
  SetAllDBData(brdf_i, num_pts, dim, tree);
  GetBRDFforAllInc(neigh, eps);
  if((int)_hemigons.size() == _v_count) _validHemigons = true;
  else _validHemigons = false;
}


//========================================================================
int
CDBDana_BRDFIsotropic::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  return maxRays;
}
//========================================================================


//========================================================================
int
CDBDana_BRDFIsotropic::SpawnStochasticRays(
  CContributionArray   &carray,
  const CHitPointInfo  &hitInfo,
  CSamplingInfo        &samplingInfo,
  EDistributionType     distType) const
{
  return 1;
}
//========================================================================


//========================================================================
bool
CDBDana_BRDFIsotropic::EvaluateBSDF(CValues &bsdf,
		           const CHitPointInfo &hitInfo,
			   const CVector3D &dirOut,
			   EDistributionType distType) const
{
  return false;
}
//========================================================================


//************************************************************************
// Functions to print computed/set class members
//========================================================================
void CDBDana_BRDFIsotropic::PrintDB_SearchTree(ostream &app) const
{
  if(_db_s_tree == NULL)
  {
    app << "Data is uninithalized! Please, apply 'InitData' method\n";
    app << "and then try again\n";
    return;
  }
  app << "Search tree for 3D coordinates (outgoing theta/phi and \n";
  app << "incoming theta angles which where read from the database):\n";
  app << "------------------------------------------------------------\n";
  _db_s_tree->Print(ANNtrue, app);
}
//========================================================================

void CDBDana_BRDFIsotropic::PrintBRDF_avg_val(ostream &app) const
{
  if((int)_BRDF_avg_val.size() > 0 && _validHemigons)
  {
    assert((int)_BRDF_avg_val.size() == _v_count);
    app << "BRDF average values (isotropic case). The number of incoming \n";
    app << "light directions is equal to " << _BRDF_avg_val.size();
    app << " (the number of vertical\n";
    app << "angles multiplied by the number of horizontal angles)\n";
    app << "------------------------------------------------------------\n";

    for(int i = 0; i < (int)_BRDF_avg_val.size(); i++)
    {
      app << "Grid for [" << i << "]->(Theta = " << _v_angles[i];
      app << ") incoming direction:\n";
      app << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
      assert((int)_BRDF_avg_val[i].size() == _h_count*_v_count);
      for(int j = 0; j < (int)_BRDF_avg_val[i].size(); j++)
      {
        app << "[" << j << "] " << _BRDF_avg_val[i][j] << "\t";
	if(!((j+1) % _v_count)) app << "\n";
      }
      app << "\n";
    }
    return;
  }
  app << "Requested _BRDF_avg_val array is empty\n";

  return;
}


//========================================================================
void CDBDana_BRDFIsotropic::PrintBRDF_rgb_val(ostream &app) const
{
  if((int)_BRDF_rgb_val.size() > 0 && _validHemigons)
  {
    assert((int)_BRDF_avg_val.size() == _v_count);
    app << "BRDF RGB values (isotropic case). The number of incoming\n";
    app << "light directions is equal to " << _BRDF_rgb_val.size();
    app << " (the number of vertical\n";
    app << "angles multiplied by the number of horizontal angles)\n";
    app << "------------------------------------------------------------\n";

    for(int i = 0; i < (int)_BRDF_rgb_val.size(); i++)
    {
      app << "Grid for [" << i << "]->(Theta = " << _v_angles[i];
      app << ") incoming direction:\n";
      app << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
      assert((int)_BRDF_rgb_val[i].size() == _h_count*_v_count);
      for(int j = 0; j < (int)_BRDF_rgb_val[i].size(); j++)
      {
        app << "[" << j << "] " << _BRDF_rgb_val[i][j] << "\t";
	if(!((j+1) % _v_count)) app << "\n";
      }
      app << "\n";
    }
    return;
 }
 app << "Requested _BRDF_rgb_val array is empty\n";

 return;
}




/**********************************************************************
 *                                                                    *
 *                    CDBBonn_BRDFAnisotropic methods                 *
 *                                                                    *
 **********************************************************************/
// describes the surface to a given stream
void
CDBBonn_BRDFAnisotropic::Describe(ostream &app, int ind) const
{
  app << "\n";
  app << "*********************************************";
  app << "*************************\n";
  indent(app, ind);
  app << "CDBBonn_BRDFAnisotropic:\n";
  app << "*********************************************";
  app << "*************************\n\n";

  if (_db_color) {
    indent(app, ind+2);
    app << "  db_color = \t";
    _db_color->Describe(app, ind); app << "\n";
  }
  /*if(_brdf_i.size() > 0) {
    //indent(app, ind+2);
    app << "BRDF data for correspondent incoming/outgoing\n";
    app << "directions (" << _brdf_i.size() << " samples):\n";
    app << "-----------------------------------------------\n";
    for(int i = 0; i < (int)_brdf_i.size(); i++)
    {
      app << "[" << i << "] ";
      for(int j = 0; j < (sizeof(_brdf_i[i].p)/sizeof((_brdf_i[i].p)[0])); j++)
        app << (_brdf_i[i].p)[j];
      app << " -> " << _brdf_i[i].brdf << "\n";
    }
  }*/
  PrintDB_SearchTree(cout);
  PrintVAngles(cout);
  PrintHAngles(cout);
  PrintGridAngles(cout);
  PrintBRDF_DB(cout);
  PrintBRDF_avg_val(cout);
  PrintBRDF_rgb_val(cout);

  app << flush;

  return;
}


//========================================================================
bool
CDBBonn_BRDFAnisotropic::GetDefaultColor(CValues &defaultColor,
			       const CHitPointInfo &hitInfo) const
{
  // there is the only color => it is the default one
  if (_db_color) {
    hitInfo.GetTextureValue(defaultColor, *_db_color, CMappingBase::EE_KD);
    return false;
  }

  defaultColor.Set(0.0);
  return true;
}
//========================================================================


// Should be checked !!!!!!!!!!!!!!!!!!!
//========================================================================
CSurface *
CDBBonn_BRDFAnisotropic::Dup() const
{
  CDBBonn_BRDFAnisotropic *ret = CreateNewObject();

  if (_db_color)
    ret->SetDB_color(_db_color);
  if (_brdf_i.size() > 0)
    ret->_brdf_i = _brdf_i;
  if(_db_s_tree)
    ret->_db_s_tree = _db_s_tree;
  if(_v_count)
    ret->_v_count = _v_count;
  if(_h_count)
    ret->_h_count = _h_count;
  if(_v_angles)
    ret->_v_angles = _v_angles;
  if(_h_angles)
    ret->_h_angles = _h_angles;
  if(_gridAngles.size() > 0)
    ret->_gridAngles = _gridAngles;
  if(_BRDF_avg_val.size() > 0)
    ret->_BRDF_avg_val = _BRDF_avg_val;
  if(_BRDF_rgb_val.size() > 0)
    ret->_BRDF_rgb_val = _BRDF_rgb_val;
  if(_hemigons.size() > 0)
    ret->_hemigons = _hemigons;

  //ret->_supportedType = _supportedType;
  return ret;
}
//========================================================================


// Read the measur_num ammount of data corresponding to the sample_num
// ID number of sample from the specified files and put it into the
// special data storage _brdf_i.
//========================================================================
void
CDBBonn_BRDFAnisotropic::ReadDataDB(int &sample_num, int &measur_num,
                                    int &rdim, vector<const char*> &names_db)
{
  ReadDataBonn(sample_num, measur_num, rdim, names_db, _brdf_i);
  _validHemigons = false;
}

// Read the measur_num ammount of data corresponding to the sample_num
// ID number of sample from the specified files and put it into the
// special data storage _brdf_i and build the searching tree trough
// the 3D coordinates.
//========================================================================
void
CDBBonn_BRDFAnisotropic::ReadDataBuildSTree(int &sample_num, int &measur_num,
                                    int &rdim, vector<const char*> &names_db)
{
  _db_s_tree = ReadDataBonnBT(sample_num, measur_num, rdim, names_db, _brdf_i);
  _validHemigons = false;

}






/**********************************************************************
 *                                                                    *
 *                    CDBBonn_BRDFIsotropic methods                 *
 *                                                                    *
 **********************************************************************/
// describes the surface to a given stream
void
CDBBonn_BRDFIsotropic::Describe(ostream &app, int ind) const
{
  app << "\n";
  app << "*********************************************";
  app << "*************************\n";
  indent(app, ind);
  app << "CDBBonn_BRDFIsotropic:\n";
  app << "*********************************************";
  app << "*************************\n\n";

  if (_db_color) {
    indent(app, ind+2);
    app << "  db_color = \t";
    _db_color->Describe(app, ind); app << "\n";
  }

  PrintDB_SearchTree(cout);
  PrintVAngles(cout);
  PrintHAngles(cout);
  PrintGridAngles(cout);
  PrintBRDF_DB(cout);
  PrintBRDF_avg_val(cout);
  PrintBRDF_rgb_val(cout);

  app << flush;

  return;
}


//========================================================================
bool
CDBBonn_BRDFIsotropic::GetDefaultColor(CValues &defaultColor,
			       const CHitPointInfo &hitInfo) const
{
  // there is the only color => it is the default one
  if (_db_color) {
    hitInfo.GetTextureValue(defaultColor, *_db_color, CMappingBase::EE_KD);
    return false;
  }

  defaultColor.Set(0.0);
  return true;
}
//========================================================================


// Should be checked !!!!!!!!!!!!!!!!!!!
//========================================================================
CSurface *
CDBBonn_BRDFIsotropic::Dup() const
{
  CDBBonn_BRDFIsotropic *ret = CreateNewObject();

  if (_db_color)
    ret->SetDB_color(_db_color);
  if (_brdf_i.size() > 0)
    ret->_brdf_i = _brdf_i;
  if(_db_s_tree)
    ret->_db_s_tree = _db_s_tree;
  if(_v_count)
    ret->_v_count = _v_count;
  if(_h_count)
    ret->_h_count = _h_count;
  if(_v_angles)
    ret->_v_angles = _v_angles;
  if(_h_angles)
    ret->_h_angles = _h_angles;
  if(_gridAngles.size() > 0)
    ret->_gridAngles = _gridAngles;
  if(_BRDF_avg_val.size() > 0)
    ret->_BRDF_avg_val = _BRDF_avg_val;
  if(_BRDF_rgb_val.size() > 0)
    ret->_BRDF_rgb_val = _BRDF_rgb_val;
  if(_hemigons.size() > 0)
    ret->_hemigons = _hemigons;

  //ret->_supportedType = _supportedType;
  return ret;
}
//========================================================================


// Read the measur_num ammount of data corresponding to the sample_num
// ID number of sample from the specified files and put it into the
// special data storage _brdf_i.
//========================================================================
void
CDBBonn_BRDFIsotropic::ReadDataDB(int &sample_num, int &measur_num,
                                    int &rdim, vector<const char*> &names_db)
{
  ReadDataBonn(sample_num, measur_num, rdim, names_db, _brdf_i);
  _validHemigons = false;
}

// Read the measur_num ammount of data corresponding to the sample_num
// ID number of sample from the specified files and put it into the
// special data storage _brdf_i and build the searching tree trough
// the 3D coordinates.
//========================================================================
void
CDBBonn_BRDFIsotropic::ReadDataBuildSTree(int &sample_num, int &measur_num,
                                    int &rdim, vector<const char*> &names_db)
{
  _db_s_tree = ReadDataBonnBT(sample_num, measur_num, rdim, names_db, _brdf_i);
  _validHemigons = false;

}




__END_BRDFLIB_SOURCE

