/* See {rdo_synthesis_basis.h}. */

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <frgb.h>
#include <frgb_ops.h>
#include <argparser.h>
#include <float_image_paint.h>
#include <dspmat.h>
#include <dspmat_linsys_GS.h>
#include <dspmat_linsys_ALT.h>

#include <rdo_geometry.h>
#include <rdo_solid.h>
#include <rdo_finish.h>
#include <rdo_site.h>
#include <rdo_lighting.h>
#include <rdo_radiosity.h>
#include <rdo_sampling.h>
#include <rdo_synthesis.h>
#include <rdo_lighting.h>

#include <rdo_synthesis_basis.h>

#define rdo_synthesis_basis_SOLVE_METHOD 0
  /* Which method to use when solving the radiosit equation: 
    0 = Gauss-Seidel, 1 = Alernative. */

void rdo_synthesis_basis_coeffs_from_site_radiances
  ( double B[],
    int n_smp, 
    dspmat_t *M,
    double L[],
    int n_bas
  )
  { 
    fprintf(stderr, "converting %d site radiances {B[]} into %d basis coeffs {L[]} ...\n", n_smp, n_bas);
    demand(M->rows == n_smp, "row count mismatch");
    demand(M->cols == n_bas, "col count mismatch");

    /* Parameters for iterative system solving: */
    int max_iter = LINSOL_MAX_ITER;  /* Max Gauss-Seidel iterations. */
    double omega = LINSOL_OMEGA;     /* Relaxation coefficient. */
    double abs_tol = LINSOL_ABS_TOL; /* Absolute stopping criterion. */
    double rel_tol = LINSOL_REL_TOL; /* Relative stopping criterion. */
    switch(rdo_synthesis_basis_SOLVE_METHOD)
      { case 0: 
          dspmat_linsys_GS_solve
            ( B, n_smp, M, L, n_bas, max_iter, omega, abs_tol, rel_tol );
          break;
        case 1:
          dspmat_linsys_ALT_solve
            ( B, n_smp, M, L, n_bas, max_iter, omega, abs_tol, rel_tol );
          break;
        default:
          demand(FALSE, "invalid linsys solution method");
      }
  }

void rdo_synthesis_basis_compute_image
  ( scene_t *scn,
    point3_t *eye, 
    dist_type_t tpDist,
    basis_elem_type_t tpElem,
    sampling_options_t *opSmp, 
    approx_basis_options_t *opBas, 
    dspmat_t *M,
    dspmat_t *T,
    synthesis_options_t *opSyn,
    int channel,
    site_vec_t *pix,
    dspmat_t *Q, 
    double F[]
  )
  { fprintf(stderr, "entering {%s} ...\n", __FUNCTION__);
    
    int n_pix = pix->ne; /* Count of image pixel sites. */
    
    /* Get hold of the scene objects and finishes: */
    solid_vec_t *sds = &(scn->sds); assert(sds->ne > 0);
    finish_vec_t *fns = &(scn->fns); assert(fns->ne > 0);

    /* Make sure that we have the sampling sites {smp}: */
    site_vec_t *smp = &(scn->smp);
    rdo_synthesis_require_sampling_sites(sds, eye, tpDist, opSmp, smp);
    int n_smp = smp->ne;
    
    /* Make sure that we have the approximation basis {bas}: */
    approx_basis_t *bas = &(scn->bas);
    rdo_synthesis_require_approx_basis(sds, smp, tpDist, tpElem, opBas, bas);
    int n_bas = bas->ne;

    demand(n_bas == n_smp, "approx basis does not match sampling sites");
    
    /* Get hold of the highlighted pixel sites and sampling sites: */
    int des_n = opSyn->highlight_n;  /* Count of seleted sites. */
    int *des_smp = &(opSyn->highlight_smp[0]);
    /*  int *des_pix = &(opSyn->highlight_pix[0]); */
    frgb_t *des_color = &(opSyn->highlight_color[0]);

    /* Vector {B} holds the surface radiance at the sampling sites: */
    double *B = NOTNULL(malloc(n_bas * sizeof(double)));

    /* Vector {L} holds the coeffs of the surface radiance in terms of the basis: */
    double *L = NOTNULL(malloc(n_bas * sizeof(double)));

    /* Vector {E} holds the coeffs of the spontaneus emission function: */
    double *E = NOTNULL(malloc(n_bas * sizeof(double)));

    /* Estimate the average number of pixels in the support of each element: */
    int avgPixPerElem = opBas->numNeighbors*n_pix/n_bas;

    /* Matrix {Q} is the rendering matrix, that maps {L} to {F}:  */
    rdo_synthesis_require_rendering_matrix(pix, bas, avgPixPerElem, Q);
    
    /* {R} is the basis radiance transfer matrix (channel-dependent; we may not need it). */
    dspmat_t R = dspmat_new(0,0,0);
    
    /* Extract the sample values {des_val[0..des_n-1]} for highlighting the selected items: */
    double des_val[des_n];
    int jd;
    for (jd = 0; jd < des_n; jd++) { des_val[jd] = des_color[jd].c[channel]; }

    /* Compute the basis coeffs {L} of the scene radiance: */
    switch (opSyn->imageType)
      {
        case image_type_BasisElems:
          { /* Turns on only the selected basis elements: */
            rdo_linalg_set_some_elements(L, n_bas, des_smp, des_val, des_n);
            break;
          }

        case image_type_Interpolators:
          { /* Sets nonzero radiance {B[i]} only at the selected sampling sites: */
            rdo_linalg_set_some_elements(B, n_smp, des_smp, des_val, des_n);
            /* Finds the coefficients of {L} that interpolate {B[i]} at the sampling sites:  */
            rdo_synthesis_require_basis_matrix(smp, bas, opBas->numNeighbors, M);
            rdo_synthesis_basis_coeffs_from_site_radiances(B, n_smp, M, L, n_bas);
            break;
          }

        case image_type_QuadrantsApprox:
          { /* Sets {B[i]} to the simple texture radiance: */
            rdo_lighting_compute_site_vec_radiances_quadrants(smp, sds, channel, B, n_smp);
            /* Finds the coefficients of {L} that interpolate {B[i]} at the sampling sites:  */
            rdo_synthesis_require_basis_matrix(smp, bas, opBas->numNeighbors, M);
            rdo_synthesis_basis_coeffs_from_site_radiances(B, n_smp, M, L, n_bas);
            break;
          }

        case image_type_ElemTransfer:
          { /* Sets selected emittance basis coeffs to 1.0: */
            rdo_linalg_set_some_elements(E, n_bas, des_smp, des_val, des_n);

            /* Initialize the total radiance {L} to {E}: */
            rdo_radiosity_initialize_basis_radiance(E, L, n_bas);
            /* Perform {opSyn->maxBounces} iteration of the radiosity algorithm: */
            int bc; /* Photon bounce counter. */
            for (bc = 0; bc < opSyn->maxBounces; bc++)
              { /* Perform a single-bounce basis radiance transfer:  */
                rdo_synthesis_require_basis_matrix(smp, bas, opBas->numNeighbors + 1, M);
                rdo_synthesis_require_site_transfer_matrix(sds, smp, T);
                rdo_synthesis_require_basis_transfer_matrix(sds, fns, smp, channel, T, M, &R);
                rdo_radiosity_basis_transfer_iteration(E, &R, L, n_bas);
              }
            break;
          }

        case image_type_Realistic:
          { /* Sets the radiance at each sampling site to be its emission coeff: */
            rdo_lighting_get_site_vec_emission_coeffs(smp, sds, fns, channel, B, n_bas);
            /* Convert the site radiances {B} to basis coeffs {E}: */
            rdo_synthesis_require_basis_matrix(smp, bas, opBas->numNeighbors, &(scn->M));
            rdo_synthesis_basis_coeffs_from_site_radiances(B, n_smp, &(scn->M), E, n_bas);
            /* Initialize the total radiance {L} to {E}: */
            rdo_radiosity_initialize_basis_radiance(E, L, n_bas);
            /* Perform {opSyn->maxBounces} iteration of the radiosity algorithm: */
            int bc; /* Photon bounce counter. */
            for (bc = 0; bc < opSyn->maxBounces; bc++)
              { /* Perform a single-bounce basis radiance transfer:  */
                rdo_synthesis_require_site_transfer_matrix(sds, smp, T);
                rdo_synthesis_require_basis_transfer_matrix(sds, fns, smp, channel, T, M, &R);
                rdo_radiosity_basis_transfer_iteration(E, &R, L, n_bas);
              }
            break;
          }

        default:
          assert(FALSE); /* Unknown option. */
      }

    rdo_synthesis_compute_pixel_radiances_from_coeffs(L, n_bas, Q, F, n_pix); 

    /* Free the temporary storage: */
    free(B);
    free(L);
    dspmat_trim(&R, 0);
  }