/*
 *  Ray++ - Object-oriented ray tracing library
 *  Copyright (C) 1998-2001 Martin Reinecke and others.
 *  See the AUTHORS file for more information.
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Library General Public
 *  License as published by the Free Software Foundation; either
 *  version 2 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Library General Public License for more details.
 *
 *  You should have received a copy of the GNU Library General Public
 *  License along with this library; if not, write to the Free
 *  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 *  See the README file for more information.
 */

/********************************************************************
This file contains the program for rendering gaseous volumes whose
density is given by a density function.
REFERENCE: Ray Tracing Volume Densities" by James T. Kajiya and 
           Brian P. Von Herzen (SIGGRAPH '84)
AUTHOR:    Neeharika Adabala

*********************************************************************/

#include "volumes/gas.h"

namespace RAYPP {

extern HANDLE<RENDERER> Renderer;

GAS::GAS ()
  : Opacity (COLOUR(float4(0.4),float4(0.4),float4(0.4))), phase(RAYLEIGH),
    min_step(float4(0.01))
  {}

GAS::GAS (const COLOUR &opac, phase_func ph,
  const HANDLE<DENSITY_FUNCTOR> &Dens)
  : Opacity (opac), phase(ph), Densfunc (Dens),
    min_step(float4(0.01))
  {}

float8 GAS::Integrate_Density_Function(const VECTOR &Lower,
  const VECTOR &Higher) const
  {
  const uint4 min_levels = 3, max_levels = 12;
  const float8 tolerance = 1.0/256;

  VECTOR diff = Higher - Lower;
  float8 step = diff.Length();

  if (step < Small_float4) return 0;

  float8 sum = 0.5* (Densfunc->Get_Dens(Higher)+Densfunc->Get_Dens(Lower));

  float8 result = step * sum;
  float8 last_result = result;

  uint4 total_intervals=1;
  float8 ptdist=1.0;

  for (uint4 iter=0; iter<max_levels; ++iter)
    {
    for (float8 offset=0.5*ptdist; offset<1.0; offset+=ptdist)
      {
      sum += Densfunc->Get_Dens(Lower+offset*diff);
      ++total_intervals;
      }
    last_result = result;
    result = sum/total_intervals*step;
    if ((iter>=min_levels) &&
      ((abs(result-last_result)<tolerance)
    || (ptdist < min_step)))
      return result;

    ptdist *= 0.5;
    }

  Log_Stream << "no convergence in density: "
    << last_result << " " << result << endl;
  return result;
  }

//virtual
void GAS::Init ()
  {
  if (initialized) return;
  if (!Densfunc)
    error ("GAS: no density function given!");
  initialized = true;
  }

//virtual
void GAS::Transform (const TRANSFORM &trans)
  {
  cni();
  Trans.Add_Transform (trans);
  }

//virtual 
COLOUR GAS::Calc_new_Importance (const RAY &Ray) const
  {
  ci();

  VECTOR Near = Trans.InvTransPoint (Ray.evalmin());
  VECTOR Far  = Trans.InvTransPoint (Ray.evalmax());

  return (Ray.Importance
    * (-Opacity*Integrate_Density_Function(Near,Far)).exp());
  }

//private
float4 GAS::Phase_Function(const float4 cosangle) const
  {
  switch (phase)
    {
    case LAMBERT:   //Lambert Surface Phase Function
      {      
      const float8 a = acos(cosangle);
      return ( 8/(3*Pi) * (sin(a) + (Pi - a) * cosangle));
      }

    case ANISOTROPIC: //Anisotropic
      {
      const float4 x = float4(0.8);
      return(1 + x * cosangle);
      }

    case RAYLEIGH: //Rayleigh Scattering
      return ( 0.75*(1+cosangle*cosangle));

    default:
      return 1.0;
    }
  }

COLOUR GAS::Get_Contrib (const VECTOR &Entry, const VECTOR &Interest,
  const VECTOR &dir) const
  {
  COLOUR term1;
  COLOUR term2(0,0,0);
  float8 term3;
  INCIDENT_ARRAY Arr;

  term1 = (-Opacity * Integrate_Density_Function(Entry, Interest)).exp();

  Renderer->Calc_Illumination (Interest, COLOUR(1,1,1), Arr);
  for(uint4 i=0; i<Arr.size(); ++i)
    {
    float8 cosangle = Dot(Arr[i].Direction,dir);
    term2 += Arr[i].Intensity * Phase_Function(cosangle);
    }
  term2 += Arr.Ambient;

  term3 = Densfunc->Get_Dens(Interest);

  return  term2 * term1 * term3;
  }

//private
COLOUR GAS::Integrate_Brightness(const VECTOR &Lower, const VECTOR &Higher,
  const VECTOR &dir) const
  {
  const uint4 min_levels = 3, max_levels = 9;
  const float8 tolerance = 1.0/256;

  VECTOR diff = Higher - Lower;
  float8 step = diff.Length();

  if (step < Small_float4) return COLOUR(0,0,0);

  COLOUR sum = 0.5* (Get_Contrib(Lower, Lower , dir)
                    +Get_Contrib(Lower, Higher, dir));

  COLOUR result = step * sum;
  COLOUR last_result = result;

  uint4 total_intervals=1;
  float8 ptdist=1.0;

  COLOUR cdiff;

  for (uint4 iter=0; iter<max_levels; ++iter)
    {
    for (float8 offset=0.5*ptdist; offset<1.0; offset+=ptdist)
      {
      sum += Get_Contrib(Lower,Lower+offset*diff, dir);
      ++total_intervals;
      }
    last_result = result;
    result = sum/total_intervals*step;

    cdiff=result-last_result;
    if ((iter>=min_levels) &&
       ((cdiff.AbsoluteSum() < tolerance)
       || (ptdist < min_step)))
      return result;

    ptdist *= 0.5;
    }

  Log_Stream << "no convergence in brightness: "
    << last_result << " " << result << endl;
  return result;
  }

//virtual
COLOUR GAS::Calc_Modified_Colour (const RAY &Ray, const COLOUR &Ingoing) const
  {
  ci();

  VECTOR Near = Trans.InvTransPoint (Ray.evalmin());
  VECTOR Far  = Trans.InvTransPoint (Ray.evalmax());

  return   Integrate_Brightness(Near, Far, Ray.dir)
         + Ingoing * (-Opacity*Integrate_Density_Function(Near,Far)).exp();
  }

//virtual 
COLOUR GAS::Get_Attenuated_Light (const RAY &Ray, const COLOUR &Ingoing) const
  {
  ci ();

  VECTOR Near = Trans.InvTransPoint (Ray.evalmin());
  VECTOR Far  = Trans.InvTransPoint (Ray.evalmax());

  return(Ingoing * (-Opacity*Integrate_Density_Function(Near,Far)).exp());
  }

void GAS::Set_Opacity  (const COLOUR &new_opacity)
  {
  cni();
  Opacity = new_opacity;
  }

void GAS::Set_Phase (phase_func pf)
  {
  cni();
  phase = pf;
  }

void GAS::Set_Density  (const HANDLE<DENSITY_FUNCTOR> &Dens)
  {
  cni();
  Densfunc = Dens;
  }

} // namespace RAYPP