/* See {rdo_raytrace.h}. */

#define rdo_raytrace_C_COPYRIGHT "Copyright © 2008 Danillo Pereira and J. Stolfi, UNICAMP"
/* Last edited on 2008-08-16 09:14:12 by stolfi */ 

#define _GNU_SOURCE
#include <stdio.h>
#include <assert.h>

#include <frgb.h>
#include <frgb_ops.h>
#include <argparser.h>
#include <filefmt.h>
#include <nget.h>
#include <fget.h>
#include <fget_geo.h>
#include <vec.h>
#include <affirm.h>
#include <bool.h>

#include <rdo_geometry.h>
#include <rdo_solid.h>
#include <rdo_raytrace.h>

#define debug_raytrace FALSE
  /* Set to TRUE to turn on debugging printouts. */

double rdo_raytrace_pt_pt_visibility(solid_vec_t *sds, point3_t *pA, point3_t *pB)  
  {
    /* Get displacement {dir} from {pA} to {pB}: */
    vector3_t dir;
    r3_sub(pB, pA, &dir);

    /* Normalize {dir} to unit length: */
    r3_dir(&dir, &dir);
    
    /* Perturb points to avoid self-occlusion: */ 
    point3_t qA, qB;
    r3_mix(1.0, pA, +VISIBILITY_TOLERANCE, &dir, &qA);
    r3_mix(1.0, pB, -VISIBILITY_TOLERANCE, &dir, &qB);

    /* Cast a ray from {qB} to {qA}: */
    int iobR; 
    double dR;
    rdo_raytrace_first_hit(sds, &qA, &dir, &dR, &iobR);

    if ((iobR == -1) || (dR > r3_dist(&qA, &qB)))
      { /* Free line of sight: */ return 1.0; }
    else if (dR <= 0.0)
      { /* Either {pA} is inside some object, or numerical error: */ return 0.0; }
    else
      { /* The view is blocked by some solid: */ return 0.0; }
  }

double rdo_raytrace_pt_site_visibility(solid_vec_t *sds, point3_t *pA, point3_t *pB, vector3_t *nB)
  {
    /* Perturb {pB} to avoid self-occlusion: */
    point3_t pBmore;
    r3_mix(1.0, pB, RAY_ORIGIN_EPSILON, nB, &pBmore);

    /* Check whether {pA} is on the correct side of the plane of {sB}: */
    if (r3_dot(pA,nB) <= r3_dot(&pBmore,nB)) { /* lado errado: */ return 0.0; }

    /* Check occlusion by other solid objects: */
    return rdo_raytrace_pt_pt_visibility(sds, pA, &pBmore);
  }

void rdo_raytrace_first_hit(solid_vec_t *sds, point3_t *eye, vector3_t *dir, double *distP, int *iobP)
  {
    double menorDist = +INF;  /* Smallest distance to some hit. */
    int isd = -1;             /* Index of solid at that hit. */

    // !!! deveria usar algum esquema de árvore ou baldes para evitar testar todos os objs.
    int i;
    for (i = 0; (menorDist > 0) && (i < sds->ne); i++) 
      { 
        double di = rdo_raytrace_hit_distance(&(sds->e[i]), eye, dir);
        if (di < menorDist) { menorDist = di; isd = i; }        
      }
    (*distP) = menorDist;
    (*iobP) = isd;
  }

void rdo_raytrace_first_hit_site(solid_vec_t *sds, point3_t *eye, vector3_t *dir, site_t *sP, double *distP)
  {
    /* Cast ray from {eye} along {dir}: */
    double dist;
    int isd;
    rdo_raytrace_first_hit(sds, eye, dir, &dist, &isd);
      
    /* Check result: */
    if (dist < 0.0)
      { /* The {eye} is inside some solid obect: */
        (*sP) = (site_t)
          { .position = (*eye),
            .normal = (vector3_t){{0,0,0}},
            .isd = isd
          };
      }
    else if (dist == +INF) 
      { /* The ray dis not hit anything */
        (*sP) = (site_t)
          { .position = (point3_t){{0,0,0}},
            .normal = (vector3_t){{0,0,0}},
            .isd = -1
          };
      }
    else
      { /* The ray hit some object: */
        point3_t position;
        r3_mix(1.0, eye, dist, dir, &position);
        assert((isd >= 0) && (isd < sds->ne));
        rdo_site_from_point(&(sds->e[isd]), isd, &position, sP);
      }

    (*distP) = dist;
  }

  
double rdo_raytrace_hit_distance(solid_t *sd, point3_t *eye, vector3_t *dir)
  {
    if (sd == NULL) { return +INF; }
    
    if (debug_raytrace)
      { fprintf(stderr, "rdo_raytrace_hit_distance:\n"); 
        fprintf(stderr, "  solid = "); rdo_solid_write(stderr, sd); fprintf(stderr, "\n"); 
        r3_gen_print(stderr, eye, "%+8.1f", "  eye =   ( ", " ", ")\n");
        r3_gen_print(stderr, dir, "%+8.5f", "  dir =   ( ", " ", ")\n");
      }
    
    double dist;
    switch(sd->type)
      {
      case solid_type_Ovoid:
        dist = rdo_raytrace_hit_distance_ovoid(sd->corner, eye, dir); break;
        
      case solid_type_Brick:
        dist = rdo_raytrace_hit_distance_brick(sd->corner, eye, dir); break;

     default:
        demand(FALSE, "invalid object type");
        dist = NAN;
      }

    if (debug_raytrace)
      { fprintf(stderr, "  dist =  %11.5f\n", dist); }
  
    return dist;
  }

double rdo_raytrace_hit_distance_ovoid(point3_t corner[], point3_t *eye, vector3_t *dir)
  { 
    /* Get the ovoid's center {ctr}: */
    point3_t ctr;
    r3_mix(+0.5, &(corner[0]), +0.5, &(corner[1]), &ctr);
    
    /* Get the semidiameters {rad[0..2]} of the ovoid: */
    vector3_t rad;
    r3_mix(-0.5, &(corner[0]), +0.5, &(corner[1]), &rad);
    
    /* Get the displacement {dif} from {ctr} to {position}: */
    vector3_t dif;        
    r3_sub(eye, &ctr, &dif);

    /* Assemble the characteristic equation from {t^2 + B*t + C = 0}: */
    double A = 0.0, B = 0.0, C = -1.0;
    int k;
    for (k = 0; k < 3; k++)
      { double ak = dir->c[k]/rad.c[k];
        double bk = (eye->c[k] - ctr.c[k])/rad.c[k];
        A += ak*ak;
        B += 2*ak*bk;
        C += bk*bk;
      }
    assert(A != 0.0);
    B /= A;
    C /= A;
     
    /* Find the roots of the characteristic equation: */
    double discr = B*B - 4*C;
    if (discr < 0) 
      { /* The straight line misses the ovoid: */
        return +INF;
      }
      
    /* Find the two real solutions: */
    double srd = sqrt(discr);
    double tmin = (-B - srd)/2; /* Smallest root. */
    double tmax = (-B + srd)/2; /* Largest root. */
    if ( tmax < 0 )
      { /* The line meets the ovoid, but the ray points away from it; */
        return +INF;
      }
    else
      { /* Return the parameter of the first hit. */
        /* Note that the parameter is negative iff {eye} is inside the ovoid. */
        return tmin;
      }
  }

double rdo_raytrace_hit_distance_brick(point3_t corner[], point3_t *eye, vector3_t *dir)
  {
    /* A brick solid is the intersection of three orthogonal slabs. */
    /* The intersection of the ray's line with each slab is an interval of {t}s. */
    /* So, we just compute the intersection of those three intervals. */
    
    double tmin = -INF;  /* Lower endpoint of interval. */
    double tmax = +INF;  /* Upper endpoint of interval. */
   
    int i;
    for (i = 0; i < 3; i++)
      { /* Grab the coordinates of {eue} and {dir} along axis {i}: */
        double oi = eye->c[i];
        double di = dir->c[i];
        
        /* Grab the brick's coordinates along axis {i}: */
        double c0 = corner[0].c[i];
        double c1 = corner[1].c[i];
        
        /* Make sure that the brick coordinates are sorted: */
        double cmin = fmin(c0, c1);
        double cmax = fmax(c0, c1);
      
        if (((oi < cmin) && (di <= 0)) || ((oi > cmax) && (di >= 0))) 
          { /* The {eye} is outside the brick and points away from it: */
            return +INF;
          }
        else if (di != 0)
          { /* Compute the ray's parameter range {[t0 _t1]} inside the slab: */
            double t0, t1; 
            if (di > 0)
              { t0 = (cmin - oi)/di; t1 = (cmax - oi)/di; }
            else
              { t0 = (cmax - oi)/di; t1 = (cmin - oi)/di; }
            assert(t0 <= t1);
            
            /* Intersect {[tmin _tmax]} with {[t0 _t1]}: */
            tmin = fmax(tmin, t0);
            tmax = fmin(tmax, t1);
            
            if (tmin >= tmax) 
              { /* The ray does not meet the brick. */
                return +INF;
              }
          }
        else
          { /* The ray is parallel to the brick -- nothing to do on this axis. */
          }
      }
      
    /* If we got this far, the ray hits the brick. */
    /* Return the first root (negative is {eye} is inside the brick): */
    return tmin;
  }


int rdo_raytrace_find_nearest_visible(solid_vec_t *sds, point3_t *eye, vector3_t *dir, site_vec_t *sts)
  {
    site_t sRef; /* Scene site visible along the given ray. */
    double dRef;
    rdo_raytrace_first_hit_site(sds, eye, dir, &sRef, &dRef);

    int isd_best = -1;

    if(sRef.isd >= 0)
      { /* Reference site {sRef} is valid, find nearest visible in {sts}: */
        double d_best = +INF;
        int i;
        for(i = 0; i < sts->ne; i++)
          { site_t *si = &(sts->e[i]);
            double di = r3_dist(&(sRef.position), &(si->position));
            if(di < d_best) 
              { double vis = rdo_raytrace_pt_site_visibility
                  (sds, eye, &(si->position), &(si->normal));
                if (vis > 0.0) { d_best = di; isd_best = i; }
              }
          }
      }

    return isd_best;       
  }