/* vide {objeto.h} */
#define objeto_C_COPYRIGHT "Copyright © 2008 Danillo Pereira e J. Stolfi, UNICAMP"
/* Last edited on 2008-07-22 20:29:33 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 "geometria_basica.h"
#include "objeto.h"

void EscreveObjeto(FILE *wr, Objeto *ob)
  { fprintf(wr, "%-12s ", ob->nome);
    int k;
    for (k = 0; k < 2; k++)
      { r3_gen_print(wr, &(ob->canto[k]), "%9.3f", " ( ", " ", " )"); }
  }

Objeto LeObjeto(FILE *rd)
  { Objeto ob;
    ob.nome = fget_string(rd);
    int k;
    for (k = 0; k < 2; k++)
      { fget_skip_and_match(rd, "(");
        ob.canto[k] = fget_r3(rd);
        fget_skip_and_match(rd, ")");
      }
    return ob;
  }

//identifica o formato atual da descrição da geometria da cena.
#define VERSAO_ARQ_LISTADEOBJETOS "2008-07-14"

void EscreveListaDeObjetos(FILE *wr, ListaDeObjetos *obs)
  {
    //escreve cabecalho:
    filefmt_write_header(wr, "ListaDeObjetos", VERSAO_ARQ_LISTADEOBJETOS);
    //aqui deveria escrever comentários fornecidos pelo cliente.

    //escreve o número de objetos:
    fprintf(wr, "num_objetos = %d\n", obs->ne);
    
    int i;
    for(i = 0; i < obs->ne; i++)
      { EscreveObjeto(wr, &(obs->e[i]));
        fprintf(wr, "\n");    
      }
    
    //escreve o rodapé e sincroniza arquivo:
    filefmt_write_footer(wr, "ListaDeObjetos");
    fflush(wr);
  }

ListaDeObjetos LeListaDeObjetos(FILE *rd)
  {
    //lê cabecalho e eventuais comentários:
    filefmt_read_header(rd, "ListaDeObjetos", VERSAO_ARQ_LISTADEOBJETOS);
    char *cmt = filefmt_read_comment(rd, '#');
    free(cmt); //por enquanto.

    //lê o número de objetos:
    int nob = nget_int(rd, "num_objetos"); fget_eol(rd);
    
    //lê os objetos:
    ListaDeObjetos obs = ListaDeObjetos_new(nob);
    int i;
    for(i = 0; i < nob; i++)
      { obs.e[i] = LeObjeto(rd);
        fget_eol(rd);    
      }
    
    //lê o rodapé:
    filefmt_read_footer(rd, "ListaDeObjetos");

    return obs;
  }

double VisibilidadePontoParaPonto(Ponto3d *pA, Ponto3d *pB, ListaDeObjetos *obs)  
  {
    //calcula vetor {dir} de {pA} para {pB}:
    Vetor3d dir;
    r3_sub(pB, pA, &dir);

    //converte {dir} em vetor unitário:
    r3_dir(&dir, &dir);
    
    //perturba ligeiramente os pontos para evitar auto-intersecção: 
    Ponto3d qA, qB;
    r3_mix(1.0, pA, +TOLERANCIA_VISIB, &dir, &qA);
    r3_mix(1.0, pB, -TOLERANCIA_VISIB, &dir, &qB);

    //lança um raio de {qB} para {qA} para verificar visibilidade:
    int iobR; 
    double dR;
    PrimeiraInterseccaoRaioObjetos(obs, &qA, &dir, &dR, &iobR);

    if ((iobR == -1) || (dR > r3_dist(&qA, &qB)))
      { /* A vista é livre: */ return 1.0; }
    else if (dR <= 0.0)
      { /* {pA} no interior de algum objeto, ou erro numérico: */ return 0.0; }
    else
      { /* Há oclusão por algum outro objeto: */ return 0.0; }
  }

void PrimeiraInterseccaoRaioObjetos(ListaDeObjetos *obs, Ponto3d *olho, Vetor3d *dir, double *distP, int *iobP)
  {
    double menorDist = +INF;  //menor distância do olho a alguma intersecção.
    int iob = -1;             //índice do objeto dessa interseção.

    //!!! deveria usar algum esquema de árvore ou baldes para evitar testar todos os objs.
    int i;
    for (i = 0; (menorDist > 0) && (i < obs->ne); i++) 
      { 
        double di = DistInterseccaoRaioObjeto(&(obs->e[i]), olho, dir);
        if (di < menorDist) { menorDist = di; iob = i; }        
      }
    (*distP) = menorDist;
    (*iobP) = iob;
  }
  
double DistInterseccaoRaioObjeto(Objeto *ob, Ponto3d *olho, Vetor3d *dir)
  {
    if (ob == NULL) { return +INF; }
    switch(ob->tipo)
      {
      case TipoDeObjeto_Bola:
        return DistInterseccaoRaioBola(ob->canto, olho, dir);
        
      case TipoDeObjeto_Cubo:
        return DistInterseccaoRaioCubo(ob->canto, olho, dir);

     default:
        demand(FALSE, "objeto de tipo inválido");
        return NAN;
      }
  }

double DistInterseccaoRaioBola(Ponto3d canto[], Ponto3d *olho, Vetor3d *dir)
  { 
    //calcula o centro {ctr} da bola:
    Ponto3d ctr;
    r3_mix(+0.5, &(canto[0]), +0.5, &(canto[1]), &ctr);
    
    //calcula o semidiametro da bola {rad} em cada eixo:
    Vetor3d rad;
    r3_mix(-0.5, &(canto[0]), +0.5, &(canto[1]), &rad);
    
    //calcula o deslocamento {dif} do centro para o {olho}
    Vetor3d dif;        
    r3_sub(olho, &ctr, &dif);

    //determina a equação característica {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 = (olho->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;
     
    //determina as raízes da equação característica:
    double discr = B*B - 4*C;
    if (discr < 0) 
      { //a reta não encontra a bola:
        return +INF;
      }
      
    //determina as duas soluções reais:
    double srd = sqrt(discr);
    double tmin = (-B - srd)/2; //menor raiz.
    double tmax = (-B + srd)/2; //maior raiz.
    if ( tmax < 0 )
      { //a reta encontra a bola, mas o raio aponta para longe da mesma:
        return +INF;
      }
    else
      { //retorna parâmetro da primeira intersecção da reta.
        //note que será negativo se e somente se o olho estiver dentro da bola.
        return tmin;
      }
  }

double DistInterseccaoRaioCubo(Ponto3d canto[], Ponto3d *olho, Vetor3d *dir)
  {
    //um {Cubo} é a intersecção de três placas ortogonais.
    //a intesecção do raio com cada placa dá um intervalo de {t}'s.
    //deve portanto retornar a intersecção desses intervalos.
    
    double tmin = -INF;  //limite inferior para a intersecção.
    double tmax = +INF;  //limite superior para a intersecção.
   
    int i;
    for (i = 0; i < 3; i++)
      { //pega as coordenadas do olho e da direção no eixo {i}:
        double oi = olho->c[i];
        double di = dir->c[i];
        
        //pega as coordenadas das faces do cubo perpendiculares ao eixo {i}:
        double c0 = canto[0].c[i];
        double c1 = canto[1].c[i];
        
        //ordena as coordenadas:
        //!!! deve exigir que entrem já ordenadas? !!!
        double cmin = fmin(c0, c1);
        double cmax = fmax(c0, c1);
      
        if (((oi < cmin) && (di <= 0)) || ((oi > cmax) && (di >= 0))) 
          { //o {olho} está fora do cubo e o raio aponta para longe:
            return +INF;
          }
        else if (di != 0)
          { //calcula o intervalo {[t0 _t1]} do parâmetro que está dentro da placa:
            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);
            
            //intersecta {[tmin _tmax]} com {[t0 _t1]}:
            tmin = fmax(tmin, t0);
            tmax = fmin(tmax, t1);
            
            if (tmin >= tmax) 
              { //o raio não intercepta o cubo:
                return +INF;
              }
          }
        else
          { // raio é paralelo à placa; nada a fazer neste eixo.
          }
      }
      
    //se chegamos até aqui, o raio intercepta o cubo.
    //retorna a primeira raiz (que pode ser negativa se o olho está dentro):
    return tmin;
  }
    
void DeterminaNormalDeObjeto(Objeto *ob, Ponto3d *posicao, Vetor3d *normal)
  {
    switch (ob->tipo)
      { 
        case TipoDeObjeto_Bola: 
          DeterminaNormalDeBola(ob->canto, posicao, normal);
          break;
        case TipoDeObjeto_Cubo:
          DeterminaNormalDeCubo(ob->canto, posicao, normal);
          break;
        default: //não deve haver outro tipo de intersecção:
          assert(FALSE);
      }
  }

void DeterminaNormalDeBola(Ponto3d canto[], Ponto3d *posicao, Vetor3d *normal)
  {
    //calcula o centro {ctr} da bola:
    Ponto3d ctr;
    r3_mix(+0.5, &(canto[0]), +0.5, &(canto[1]), &ctr);
    
    //calcula o semidiametro da bola {rad} em cada eixo:
    Vetor3d rad;
    r3_mix(-0.5, &(canto[0]), +0.5, &(canto[1]), &rad);
    
    //calcula o deslocamento a partir do centro:
    r3_sub(posicao, &ctr, normal);
    
    //calcula o gradiente:
    int k;
    for (k = 0; k < 3; k++) { normal->c[k] /= (rad.c[k]*rad.c[k]); }
    
    //normaliza o gradiente:
    r3_dir(normal, normal);
  }

void DeterminaNormalDeCubo(Ponto3d canto[], Ponto3d *posicao, Vetor3d *normal)
  {
    //determina o centro do cubo:
    Ponto3d centro;
    r3_mix(0.5, &(canto[0]), 0.5, &(canto[1]), &centro);

    //determina o vetor do centro ao ponto dado:
    Vetor3d dif;
    r3_sub(posicao, &centro, &dif);

    //determina o eixo {imax} perpendicular à face one o ponto está:
    int imax = -1;    
    double rmax = 0;  //distância relativa com maior valor absoluto.
    int i;
    for (i = 0; i < 3; i++)
      { //pega dados references às projeções do cubo e do ponto no eixo {i}:
        double c0 = canto[0].c[i];   //uma ponta do cubo.
        double c1 = canto[1].c[i];   //a outra ponta.
        double centro = (c0 + c1)/2;      //centro do cubo.
        double raio = fabs(c0 - c1)/2;    //semi-largura do cubo.
        assert(raio > 0);
        double dist = posicao->c[i] - centro; //distância (com sinal) do centro ao ponto.
        double r = dist/raio;                 //distância relativa ao raio.
        if (fabs(r) > fabs(rmax)) { rmax = r; imax = i; }
      }
    assert(rmax != 0); //senão o ponto dado é o centro do cubo.
    assert(imax >= 0); //idem.

    //a normal é a direção do eixo {imax} com o sinal de {rmax}:
    (*normal) = (Vetor3d){{ 0,0,0 }};
    normal->c[imax] = copysign(1.0, rmax);
  }