/**************************************************************************
  Copyright (C) 1992 Guy Moreillon
  All rights reserved.

  This software may be freely copied, modified, and redistributed
  provided that this copyright notice is preserved on all copies.

  You may not distribute this software, in whole or in part, as part of
  any commercial product without the express consent of the authors.

  There is no warranty or other guarantee of fitness of this software
  for any purpose.  It is provided solely "as is".
**************************************************************************/
/**************************************************************************

				RAD
		       (Interactive Radiosity)
		    	    Version 1.0
			       1992		    	    
		    	    
		  
		    	  Guy Moreillon
		    	  
**************************************************************************/

/**************************************************************************
  Fichier	: ray.c
  Description	: fonctions de calcul des ombres portees (raytracing)
**************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <malloc.h>
#include "types.h"
#include "macros.h"
#include "rad.h"
#include "radio.h"
#include "ray.h"

#define INIT_FAR	10000.0
#define ENL_PRCNT	0.01

BOOLEAN check_inter(VOXELP	voxel,
		    TRIANGLEP	tri)
/**************************************************************************
  But	: Determine si le triangle tri intersecte le voxel
  Entree: voxel :   Le voxel
	  tri :	    Le triangle
  Sortie: TRUE si tri intersecte voxel, FALSE sinon
**************************************************************************/
{
  VECTOR    sc;
  INTERVAL  Vu, Vs;
  VOXEL	    vox;

  Vox_sub_v(vox, (*voxel), tri->tri.poly.vertex[0].point);
  Vox_dot_v(Vu, vox, tri->tri.poly.normal);

  if ((Vu.inf > 0.0) || (Vu.sup < 0.0))
    return FALSE;

  Vox_dot_v(Vs, vox, tri->tri.poly.sa);
  if ((Vs.sup < 0.0) || (Vs.inf > 1.0))
    return FALSE;

  Vox_dot_v(Vs, vox, tri->tri.poly.sb);
  if ((Vs.sup < 0.0) || (Vs.inf > 1.0))
    return FALSE;

  V_add(sc, tri->tri.poly.sa, tri->tri.poly.sb);
  Vox_dot_v(Vs, vox, sc);
  if ((Vs.sup < 0.0) || (Vs.inf > 1.0))
    return FALSE;

  return TRUE;
}

#define SP sizeof(POLY_LP)
#define	ACCESS(I, J, K)	(grid->nb_vox[2]*(grid->nb_vox[1]*I + J) + K)

void digit_tri(GRIDP	    grid,
	       TRIANGLEP    tri)
/**************************************************************************
  But	: Met le triangle tri dans la liste de tous les voxels qu'il
	  intersecte dans la grille grid
  Entree: grid :    La grille
	  tri :	    Le triangle
  Sortie: neant
**************************************************************************/
{
  int	    imin[3], imax[3];
  VECTOR    bboxi, bboxa;

  register  int	    i, j, k;
  register  VERTEXP vertex = tri->tri.poly.vertex;

  V_init(bboxi, HUGE, HUGE, HUGE);
  V_init(bboxa, -HUGE, -HUGE, -HUGE);
  for(i = 0; i < 3; i++)
    {
      V_min(bboxi, bboxi, vertex[i].point);
      V_max(bboxa, bboxa, vertex[i].point);
    };
  
  for(i = 0; i < 3; i++)
    {
      imin[i] = (int)ffloor((bboxi[i]-grid->org[i])/grid->size[i]);
      imax[i] = (int)ffloor((bboxa[i]-grid->org[i])/grid->size[i]);
    };

  for(i = imin[0]; i <= imax[0]; i++)
    for(j = imin[1]; j <= imax[1]; j++)
      for(k = imin[2]; k <= imax[2]; k++)
	{
	  VOXEL vox;

	  vox.x.inf = grid->org[0] + (float )i*grid->size[0];
	  vox.x.sup = vox.x.inf + grid->size[0];
	  vox.y.inf = grid->org[1] + (float )j*grid->size[1];
	  vox.y.sup = vox.y.inf + grid->size[1];
	  vox.z.inf = grid->org[2] + (float )k*grid->size[2];
	  vox.z.sup = vox.z.inf + grid->size[2];

	  if (check_inter(&vox, tri))
	    {
	      POLY_LP	tmp;
	      int	vox = ACCESS(i, j, k);

	      tmp = (POLY_LP)malloc(sizeof(POLY_L));
	      assert(tmp);
	      tmp->tri  = tri;
	      tmp->next = grid->list[vox];
	      grid->list[vox] = tmp;
	    };
	};
}

void create_grid(OBJECT_OP  obj)
/**************************************************************************
  But	: Cree la grille de voxels de l'objet composite passe en parametre,
	  et la rempli
  Entree: obj :	    l'objet
  Sortie: neant
**************************************************************************/
{
  float     smax;
  GRID	    *grid;
  VECTOR    dum, size;
  OBJECT_OP pure;

  register  int	    i, si;

  grid = (GRID *)malloc(sizeof(GRID));
  assert(grid);

  V_sub(size, obj->bbox[1],
	      obj->bbox[0]);
  V_scale(dum, ENL_PRCNT, size);
  V_add_self(size, dum);

  if (size[0] > size[1])
    if (size[0] > size[2])
      si = 0;
    else
      si = 2;
  else
    if (size[1] > size[2])
      si = 1;
    else
      si = 2;
  smax = size[si] / (float )GRID_SIZE;

  for(i = 0; i < 3; i++)
    if (i == si)
      grid->nb_vox[i] = GRID_SIZE;
    else
      grid->nb_vox[i] = (int)ffloor(size[i]/smax) + 1;

  V_init(grid->size, smax, smax, smax);

  V_scale(dum, 0.5, dum);
  V_sub(grid->org, obj->bbox[0], dum);

  grid->list = (POLY_LP *)calloc(grid->nb_vox[0]*grid->nb_vox[1]*grid->nb_vox[2], SP);
  assert(grid->list);
  
  if (obj->type == composite_type)
    {
      for(pure = obj->object.composite.sons; pure != NULL; pure = pure->next)
        for(i = 0; i < pure->object.pure.nb_polys; i++)
          digit_tri(grid, &(pure->object.pure.polys[i]));
      obj->object.composite.grid = grid;
    }
  else
    {
      for(i = 0; i < obj->object.pure.nb_polys; i++)
        digit_tri(grid, &(obj->object.pure.polys[i]));
      obj->object.pure.grid = grid;
    };
}

void init_grids(OBJECT_OP   object)
/**************************************************************************
  But	: Initialise la grille de voxels de chaque objet composite ou pure,
	  si celui-ci ne fait pas partie d'un composite
  Entree: object :    l'objet
  Sortie: neant
**************************************************************************/
{
  register  OBJECT_OP obj;

  if (object == NULL)
    return;

  switch (object->type) {
    case assembly_type :
      for(obj = object->object.assembly.sons; obj != NULL; obj = obj->next)
        init_grids(obj);
      break;
    case composite_type :
    case pure_type :
      create_grid(object);
      break;
    default :
      Bye_bye("Major bug: wrong hierarchy\n");
  }
}

BOOLEAN ray_hit_bbox(RAYP   ray,
		     VECTOR bbox[2],
		     float  dis,
		     float  *ti,
		     float  *ta)
/**************************************************************************
  But	: Determine si le rayon ray intersecte la bounding box bbox, et
	  pour quelles valeurs de t (ray = org + t*dir)
  Entree: ray :	    le rayon
	  bbox :    la bounding box
	  dis :	    la distance maximale
	  ti, ta :  les valeurs min et max de t
  Sortie: TRUE s'il touche, FALSE sinon
**************************************************************************/
{
  float     t, bi, ba;

  register  int	    i;

  nb_bbox_trials++;

  *ta = dis;
  *ti = 0.0;

  for(i = 0; i < 3; i++)
    {
      if (ray->dir[i] >= 0.0)
	bi = bbox[0][i], ba = bbox[1][i];
      else
	bi = bbox[1][i], ba = bbox[0][i];

      if ((t = (ba - ray->org[i])/ray->dir[i]) < 0.0)
	return FALSE;

      if (t < *ta)
	{
	  if (t < *ti)
	    return FALSE;
	  *ta = t;
	};

      if ((t = (bi - ray->org[i])/ray->dir[i]) >= 0.0)
	if (t >= *ti)
	  {
	    if (t > *ta)
	      return FALSE;
	    *ti = t;
	  };
    };
  
  nb_bbox_hits++;

  return TRUE;
}

BOOLEAN ray_hit_tri(RAYP	ray,
		    TRIANGLEP	tri,
		    float 	ti,
		    float 	ta,
		    float 	*t)
/**************************************************************************
  But	: Determine si le rayon ray touche le polygone tri, et rend la
	  la valeur de t correspondante. t doit etre compris entre ti et ta
  Entree: ray :	    le rayon
	  tri :	    le polygone
	  ti, ta :  les valeurs min et max pour t
	  t :	    la valeur (rendue) de t au point d'intersection
  Sortie: TRUE s'il touche, FALSE sinon
**************************************************************************/
{
  VECTOR    p, sc;
  float     al, be, ga;

  nb_tests++;

  *t = (V_dot(tri->tri.poly.normal, tri->tri.poly.vertex[0].point) -
	V_dot(tri->tri.poly.normal, ray->org)) /
       V_dot(tri->tri.poly.normal, ray->dir);
  if ((*t < ti) || (*t > ta))
    return FALSE;

  V_add_mul(p, ray->org, *t, ray->dir);
  V_sub(p, p, tri->tri.poly.vertex[0].point);

  al = V_dot(p, tri->tri.poly.sa);
  if ((al < 0.0) || (al > 1.0))
    return FALSE;

  be = V_dot(p, tri->tri.poly.sb);
  if ((be < 0.0) || (be > 1.0))
    return FALSE;

  V_add(sc, tri->tri.poly.sa, tri->tri.poly.sb);
  ga = V_dot(p, sc);
  if ((ga < 0.0) || (ga > 1.0))
    return FALSE;

  nb_hits++;

  return TRUE;
}

BOOLEAN first_p_hit(RAYP    ray,
		    POLY_LP list,
		    float   ti,
		    float   ta)
/**************************************************************************
  But	: Determine si le rayon ray touche un polygone de la liste list
	  avant la distance maximale ta et apres la distance min ti
  Entree: ray :	    le rayon
	  list :    la liste des polygones
	  ti, ta :  les distances min et max
  Sortie: TRUE s'il touche, FALSE sinon
**************************************************************************/
{
  POLY_LP   l;
  float     t;

  for(l = list; l != NULL; l = l->next)
    if (ray_hit_tri(ray, l->tri, ti, ta, &t))
      return TRUE;

  return FALSE;
}

#define ind_of_min(RES, T1, T2, T3) \
    {\
      float  ft[3];\
      ft[0] = _ABS(T1);\
      ft[1] = _ABS(T2);\
      ft[2] = _ABS(T3);\
      if ((ft[0] < ft[1]) && (ft[0] < ft[2]))\
	RES = 0;\
      else if (ft[1] < ft[2])\
	RES = 1;\
      else\
	RES = 2;\
    }

BOOLEAN first_hit(RAYP	    ray,
		  float     ti,
		  float     ta,
		  GRID	    *grid)
/**************************************************************************
  But	: Determine si le rayon ray touche un polygone dans la grille grid
	  avant la distance maximale ta et apres la distance min ti
  Entree: ray :	    le rayon
	  ti, ta :  les distances min et max
	  grid :    la grille
  Sortie: TRUE s'il touche, FALSE sinon
**************************************************************************/
{
  POLY_LP   list;
  VECTOR    p, t, d, dum;
  float     t0, t1;
  int	    g[3];

  register  int	    i, ind;

  t0 = ti;
  V_add_mul(p, ray->org, t0, ray->dir);
  V_sub(dum, p, grid->org);
  for(i = 0; i < 3; i++)
    {
      g[i] = (int)floor((double)dum[i]/(double)grid->size[i]);
      if ((g[i] < 0) || (g[i] >= grid->nb_vox[i]))
	return FALSE;
    };
  
  for(i = 0; i < 3; i++)
    {
      t[i] = (float)((double)(grid->org[i] + g[i]*grid->size[i] - ray->org[i])
		    /(double)ray->dir[i]);
      d[i] = (float)((double)grid->size[i]/(double)_ABS(ray->dir[i]));
      if (ray->dir[i] >= 0.0)
	t[i] += d[i];
    };

  while ((t0 < ta) &&
	 (g[0] < grid->nb_vox[0]) && (g[0] >= 0) &&
	 (g[1] < grid->nb_vox[1]) && (g[1] >= 0) &&
	 (g[2] < grid->nb_vox[2]) && (g[2] >= 0))
    {
      ind_of_min(ind, t[0], t[1], t[2]);
      t1 = t[ind];

      if ((list = grid->list[ACCESS(g[0], g[1], g[2])]) != NULL)
	if (first_p_hit(ray, list, ti, ta))
	  return TRUE;

      t0 = t1;
      t[ind] += d[ind];

      if (ray->dir[ind] > 0.0)
	g[ind]++;
      else if (ray->dir[ind] < 0.0)
	g[ind]--;
    };

  return FALSE;
}

BOOLEAN object_hit(RAYP		ray,
		   OBJECT_OP	object,
		   float 	dis)
/**************************************************************************
  But	: Determine si le rayon ray touche l'objet obj avant la distance
	  maximale dis
  Entree: ray :	    le rayon
	  object :  l'objet
	  dis :	    la distance maximale
  Sortie: TRUE s'il touche, FALSE sinon
**************************************************************************/
{
  OBJECT_OP obj;
  BOOLEAN   hit;
  float     ti, ta;

  if (ray_hit_bbox(ray, object->bbox, dis, &ti, &ta))
    switch (object->type) {
      case assembly_type :
        for(obj = object->object.assembly.sons; obj != NULL; obj = obj->next)
          if (hit = object_hit(ray, obj, dis))
            break;
        return hit;
      case composite_type :
        return (first_hit(ray, ti, ta, object->object.composite.grid));
      case pure_type :
        if (!object->object.pure.shadows)
	  {
	    nb_shad_test_av++;  
            return FALSE;   
	  }
        return (first_hit(ray, ti, ta, object->object.pure.grid));
      default :
	Bye_bye("Wrong object type !!\n");
    }
  else
    return FALSE;
}

BOOLEAN shadow(RAY	*ray,
	       SCENE_OP	scene,
	       float 	dis)
/**************************************************************************
  But	: Determine si le rayon ray touche un objet avant la distance max
  Entree: ray	    : le rayon
	  scene	    : la scene
	  dis	    : la distance maximale
  Sortie: TRUE si le point est a l'ombre, FALSE sinon
**************************************************************************/
{
  nb_rays++;

  return (object_hit(ray, scene->scene, dis));
}

BOOLEAN closest_p_hit(SELECTIONP    select,
		      RAYP	    ray,
		      POLY_LP	    list,
		      float	    ti,
		      float	    ta)
/**************************************************************************
  But	: Determine le polygone de la liste list le plus proche que le
	  rayon ray touche avant la distance maximale ta et apres la
	  distance min ti
  Entree: select :  la struct de selection partiellement completee (rendue)
	  ray :	    le rayon
	  list :    la liste des polygones
	  ti, ta :  les distances min et max
  Sortie: TRUE si le rayon touche un polygone, FALSE sinon
**************************************************************************/
{
  POLY_LP   l;
  float     t;

  select->tri  = NULL;
  select->dist = INIT_FAR;

  for(l = list; l != NULL; l = l->next)
    if (ray_hit_tri(ray, l->tri, ti, ta, &t))
      if (t < select->dist)
	{
	  select->dist = t;
	  select->tri  = l->tri;
	}
  
  if (select->tri != NULL)
    return TRUE;
  else
    return FALSE;
}

BOOLEAN closest_hit(SELECTIONP	select,
		    RAYP	ray,
		    float	ti,
		    float	ta,
		    GRIDP	grid)
/**************************************************************************
  But	: Determine le polygone le plus proche que le rayon ray touche
	  dans la grille grid avant la distance maximale ta et apres la
	  distance min ti
  Entree: select :  la struct de selection partiellement completee
	  ray :	    le rayon
	  ti, ta :  les distances min et max
	  grid :    la grille
  Sortie: TRUE si le rayon touche un polygone, FALSE sinon
**************************************************************************/
{
  POLY_LP   list;
  VECTOR    p, t, d, dum;
  float     t0, t1;
  BOOLEAN   no_hit = TRUE;
  int	    g[3];

  register  int	    i, ind;

  t0 = ti;
  V_add_mul(p, ray->org, t0, ray->dir);
  V_sub(dum, p, grid->org);
  for(i = 0; i < 3; i++)
    {
      g[i] = (int)floor((double)dum[i]/(double)grid->size[i]);
      if ((g[i] < 0) || (g[i] >= grid->nb_vox[i]))
	return FALSE;
    };
  
  for(i = 0; i < 3; i++)
    {
      t[i] = (float)((double)(grid->org[i] + g[i]*grid->size[i] - ray->org[i])
		    /(double)ray->dir[i]);
      d[i] = (float)((double)grid->size[i]/(double)_ABS(ray->dir[i]));
      if (ray->dir[i] >= 0.0)
	t[i] += d[i];
    };

  while ((t0 < ta) &&
	 (g[0] < grid->nb_vox[0]) && (g[0] >= 0) &&
	 (g[1] < grid->nb_vox[1]) && (g[1] >= 0) &&
	 (g[2] < grid->nb_vox[2]) && (g[2] >= 0) &&
	 no_hit)
    {
      ind_of_min(ind, t[0], t[1], t[2]);
      t1 = t[ind];

      if ((list = grid->list[ACCESS(g[0], g[1], g[2])]) != NULL)
	if (closest_p_hit(select, ray, list, ti, t1))
	  no_hit = FALSE;

      t0 = t1;
      t[ind] += d[ind];

      if (ray->dir[ind] > 0.0)
	g[ind]++;
      else if (ray->dir[ind] < 0.0)
	g[ind]--;
    };

  return (!no_hit);
}

BOOLEAN closest_object_hit(SELECTIONP	select,
			   RAYP		ray,
			   OBJECT_OP	object,
			   float 	dis)
/**************************************************************************
  But	: Determine si le rayon ray touche l'objet obj avant la distance
	  maximale dis, et si cet objet est un assemblage ou un composite,
	  il determine quel objet pur il touche en premier
  Entree: select :  la struct de selection partiellement completee
	  ray :	    le rayon
	  object :  l'objet
	  dis :	    la distance maximale
  Sortie: TRUE s'il touche, FALSE sinon
**************************************************************************/
{
  OBJECT_OP obj;
  SELECTION s;
  BOOLEAN   hit = FALSE;
  float     ti, ta;

  select->dist = INIT_FAR;

  if (ray_hit_bbox(ray, object->bbox, dis, &ti, &ta))
    switch (object->type) {
      case assembly_type :
        for(obj = object->object.assembly.sons; obj != NULL; obj = obj->next)
          if (closest_object_hit(&s, ray, obj, dis))
	    if (s.dist < select->dist)
	      {
	        *select	      = s;
		hit	      = TRUE;
	      }
        break;
      case composite_type :
        if (hit = closest_hit(select, ray, ti, ta, object->object.composite.grid))
	  select->pure = select->tri->dad;
        break;
      case pure_type :
        if (hit = closest_hit(select, ray, ti, ta, object->object.pure.grid))
	  select->pure = object;
	break;
      default :
	Bye_bye("Wrong object type !!\n");
    }
  else
    select->pure = NULL;

  return hit;
}

BOOLEAN object_in_sight(SELECTIONP  select,
			RAYP	    ray,
			SCENE_OP    scene,
			float	    dis)
/**************************************************************************
  But	: Determine s'il y un objet dans la ligne de visee
  Entree: select    : la structure de selection completee
	  ray	    : le rayon
	  scene	    : la scene
	  dis	    : la distance maximale
  Sortie: TRUE s'il existe un objet dans la ligne de visee, FALSE sinon
**************************************************************************/
{
  SELECTION s;
  BOOLEAN   hit = FALSE;

  select->dist = INIT_FAR;

  if (closest_object_hit(&s, ray, scene->scene, dis))
    if (s.dist < select->dist)
      {
     	*select = s;
	hit     = TRUE;
      }

  if (hit)
    {
      if (select->tri->tri.poly.sons != NULL)
        {
          SAMPLE sample;

          V_add_mul(sample.point, ray->org, select->dist, ray->dir);
          find_subdiv(select->tri->tri.poly.sons, &sample);
          select->elem = sample.el;
        }
      else
        select->elem = NULL;
    };

  return hit;
}

#undef	ind_of_min
#undef	ACCESS
