/* See {rdo_lighting.h}. */

#define rdo_lighting_C_COPYRIGHT "Copyright © 2008 Danillo Pereira and J. Stolfi, UNICAMP"
/* Last edited on 2024-12-21 11:52:16 by stolfi */ 

#include <stdio.h>

#include <dspmat.h>

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

#include <rdo_lighting.h>

finish_t *rdo_lighting_get_site_finish(site_t *st, solid_vec_t *sds, finish_vec_t *fns)
  { if (st == NULL) { return NULL; }
    int isd = st->isd;
    if (isd < 0) { return NULL; }
    affirm(isd < sds->ne, "site with invalid solid object index");
    int ifn = sds->e[isd].ifn;
    if (ifn < 0) { return NULL; }
    affirm(ifn < fns->ne, "solid objects with invalid finish index");
    return &(fns->e[ifn]);
  }

dspmat_t rdo_lighting_compute_site_transfer_matrix
  ( solid_vec_t *sds,
    site_t so[], 
    double ao[],
    int no, 
    site_t sd[], 
    int nd
  )
  { dspmat_t T = dspmat_new(nd, no, 10*no);
    int i;
    dspmat_pos_t posT = 0;
    for(i = 0; i < nd; i++)
      { /* Get hold of site {sd[i]}: */
        point3_t *pdi = &(sd[i].position);
        vector3_t *ndi = &(sd[i].normal);
        int j;
        for(j = 0; j < no; j++)
          { /* Get hold of site {so[j]} and its area {ao[j]: */
            point3_t *poj = &(so[j].position);
            vector3_t *noj = &(so[j].normal);
            double Tij = ao[j] * rdo_lighting_compute_radiance_coeff_between_sites(sds, poj, noj, pdi, ndi);
            if (Tij != 0.0) 
              { posT = dspmat_add_element(&T, posT, i, j, Tij); }
          } 
      }
    dspmat_trim(&T, posT);
    return T;
  }

void rdo_lighting_apply_diffusion_coeffs_to_matrix
  ( dspmat_t *S,
    site_vec_t *sts, 
    solid_vec_t *sds, 
    finish_vec_t *fns, 
    int channel
  )
  { 
    assert(sds->ne > 0);
    assert(fns->ne > 0); 
    assert(sts->ne > 0); 
    assert(S->ents > 0); 
    assert(S->rows == sts->ne); 
    
    int nB = sts->ne;
    double *B = NOTNULL(malloc(nB * sizeof(double)));
    rdo_lighting_get_site_vec_diffusion_coeffs(sts, sds, fns, channel, B, nB);
    dspmat_pos_t posS;
    for (posS = 0; posS < S->ents; posS++)
      { dspmat_entry_t *eP = &(S->e[posS]);
        assert(eP->col < S->cols);
        assert(eP->row < S->rows);
        eP->val *= B[eP->row];
      }
    free(B);
  }

double rdo_lighting_compute_radiance_coeff_between_sites
  ( solid_vec_t *sds, 
    point3_t *pO, 
    vector3_t *nO, 
    point3_t *pD, 
    vector3_t *nD
  )
  { /* Get the vector from {pO} to {pD}, and the squared distance {dOD2} between them: */
    vector3_t rOD; r3_sub(pD, pO, &rOD);
    double dOD2 = r3_norm_sqr(&rOD);
    if (dOD2 < VISIBILITY_TOLERANCE * VISIBILITY_TOLERANCE) 
      { /* The sites seem to be the same site: */
        return 0.0; /* No transfer to self. */
      }
      
    /* Get the projections {cO,cD} of {rOD} on the normal directions: */
    double cO = +r3_dot(&rOD, nO); if (cO <= 0.0) { return 0.0; }
    double cD = -r3_dot(&rOD, nD); if (cD <= 0.0) { return 0.0; }
    /* Perturb slightly the two points away from the surface: */
    point3_t qO, qD; /* Perturbed points. */
    r3_mix(RAY_ORIGIN_EPSILON, nO, 1.0, pO, &qO);
    r3_mix(RAY_ORIGIN_EPSILON, nD, 1.0, pD, &qD);
    /* Compute the visibility factor between the two sites: */
    double fvis = rdo_raytrace_pt_pt_visibility(sds, &qO, &qD);
    /* Combine the factors: */
    return fvis*cO*cD/(dOD2*dOD2);
  }

void rdo_lighting_get_site_vec_diffusion_coeffs
  ( site_vec_t *sts, 
    solid_vec_t *sds,
    finish_vec_t *fns,
    int channel, 
    double B[], 
    int nB
  )
  {
    demand(nB == sts->ne, "incompatible {sts->ne} and {nB}");
    /* Loop over the sites: */
    int k;
    for (k = 0; k < sts->ne; k++)
      { /* Grab the site {stk}: */
        site_t *stk = &(sts->e[k]);
        /* Get its finish {fn}: */
        finish_t *fn = rdo_lighting_get_site_finish(stk, sds, fns);
        B[k] = rdo_finish_get_diffusion_coeff(fn, channel);
      }
  }

void rdo_lighting_get_site_vec_emission_coeffs
  ( site_vec_t *sts, 
    solid_vec_t *sds, 
    finish_vec_t *fns, 
    int channel,
    double B[], 
    int nB
  )
  { demand(nB == sts->ne, "incompatible {sts->ne} and {nB}");
    /* Loop over the sites: */
    int k;
    for (k = 0; k < nB; k++)
      { /* Grab the site {sk}: */
        site_t *stk = &(sts->e[k]);
        /* Get its finish {fn}: */
        finish_t *fn = rdo_lighting_get_site_finish(stk, sds, fns);
        B[k] = rdo_finish_get_emission_coeff(fn, channel);
      }
  }

void rdo_lighting_compute_site_vec_radiances_quadrants
  ( site_vec_t *sts, 
    solid_vec_t *sds, 
    int channel, 
    double B[], 
    int nB
  )
  {
    demand(nB == sts->ne, "incompatible {sds->ne} and {nB}");

    /* Loop over the sites: */
    int k;
    for (k = 0; k < sts->ne; k++)
      { /* Grab the site {stk}: */
        site_t *stk = &(sts->e[k]);
        /* Get the solid {sdk}: */
        int isd = stk->isd;
        solid_t *sdk = (isd < 0 ? NULL : &(sds->e[isd]));
        B[k] = rdo_lighting_compute_site_radiance_quadrants(stk, sdk, channel);
      }
  }

double rdo_lighting_compute_site_radiance_quadrants(site_t *st, solid_t *sd, int channel)
  {
    frgb_t color;
    if (sd != NULL)
      { /* A valid intersection, get the normal {u}: */
        vector3_t *u = &st->normal;
        
        /* Default color coordinates: */
        color = (frgb_t){{0.150f, 0.150f, 0.150f}};
        
        /* Find the axis {k} where the normal is largest: */
        int k = 0;
        int j;
        for (j = 1; j < 3; j++) { if (fabs(u->c[j]) > fabs(u->c[k])) { k = j; } }
        
        /* !!! change the formulas to depend on solid-relative coords rather than normal. !!! */
        
        if (sd->type == solid_type_Brick)
          { /* Just paint each face with a different color: */
            color.c[k] = 1.0;
            if (u->c[k] < 0) { for (j = 0; j < 3; j++) { color.c[j] = (float)(1.0 - color.c[j]); } }
          }
        else if (sd->type == solid_type_Ovoid)
          { /* Start with a brick-like color schema: */
            color.c[k] = 1.0;
            /* Do something different near the corners of the brick: */
            double t = fabs(u->c[0]) + fabs(u->c[1]) + fabs(u->c[2]);
            if (t > 2*fabs(u->c[k]))
              { color.c[k] = 0.150f; }
            else if (u->c[k] < 0)
              { for (j = 0; j < 3; j++) { color.c[j] = (float)(1.0 - color.c[j]); } }
          }
      }
    else
      { /* Site lies on no object: */
        color = (frgb_t){{0.250, 0.250, 0.250}};
      }  
    return color.c[channel];
  }