/* See {rdo_synthesis.h}. */

#define rdo_synthesis_C_COPYRIGHT "Copyright © 2008 Danillo Pereira and J. Stolfi, UNICAMP"
/* Last edited on 2024-12-21 11:51:11 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 <rdo_geometry.h>
#include <rdo_solid.h>
#include <rdo_finish.h>
#include <rdo_site.h>
#include <rdo_lighting.h>
#include <rdo_camera.h>
#include <rdo_radiosity.h>
#include <rdo_sampling.h>
#include <rdo_synthesis_plain.h>
#include <rdo_synthesis_basis.h>

#include <rdo_synthesis.h>

#define debug_synthesis (FALSE)
  /* Define this as {TRUE} to get some debugging printouts. */
  
void rdo_synthesis_require_sampling_sites
  ( solid_vec_t *sds, 
    point3_t *eye,
    dist_type_t tpDist,
    sampling_options_t *opSmp,
    site_vec_t *smp
  )
  { assert(sds->ne > 0);
    if (smp->ne == 0)
      { /* Choose the sampling sites: */
        (*smp) = rdo_sampling_choose_approx_basis_centers(sds, tpDist, opSmp, eye); 
      }
    fprintf(stderr, "sampling site list {smp} has %d sites\n", smp->ne);
  }
  
void rdo_synthesis_require_approx_basis
  ( solid_vec_t *sds,
    site_vec_t *smp, 
    dist_type_t tpDist,
    basis_elem_type_t tpElem, 
    approx_basis_options_t *opBas,
    approx_basis_t *bas
  )
  { assert(sds->ne > 0);
    assert(smp->ne > 0); 
    if (bas->ne == 0) 
      { double rd[smp->ne];
        int k;
        for (k = 0; k < smp->ne; k++)
          { rd[k] = rdo_approx_basis_pick_elem_radius(tpDist, smp, k, opBas->numNeighbors);
            demand(isfinite(rd[k]), "coud not determine the element's nominal radius");
          }
        (*bas) = rdo_approx_basis_create(tpDist, tpElem, smp->e, rd, smp->ne);
      }
    fprintf(stderr, "approx basis {bas} has %d elements\n", bas->ne);
    assert(bas->ne == smp->ne);
  }
  
void rdo_synthesis_require_basis_matrix
  ( site_vec_t *smp,
    approx_basis_t *bas,
    int avgElemsPerColumn,
    dspmat_t *M
  )
  { assert(smp->ne > 0); 
    assert(bas->ne > 0); 
    if (M->ents == 0) 
      { fprintf(stderr, "computing the basis matrix {M}...\n");
        (*M) = rdo_approx_basis_build_matrix(smp, bas, avgElemsPerColumn);
        dspmat_sort_entries(M, +2, 0);
        /* !!! Should check it for diagonal dominance !!! */
      }
    double frow = ((double)M->ents)/((double)M->rows);
    fprintf(stderr, "basis matrix {M} (%d × %d)", M->rows, M->cols);
    fprintf(stderr, " has %d nonzero elements (%.1f per row)\n", M->ents, frow);
    assert(M->rows == smp->ne);
    assert(M->cols == bas->ne);
  }

void rdo_synthesis_require_site_transfer_matrix
  ( solid_vec_t *sds,
    site_vec_t *smp, 
    dspmat_t *T
  )
  { assert(sds->ne > 0);
    assert(smp->ne > 0); 
    if (T->ents == 0) 
      { fprintf(stderr, "computing the site-to-site radiance transfer matrix {T}...\n");
        int ns = smp->ne;
        double *ar = rdo_site_estimate_all_effective_areas(smp);
        (*T) = rdo_lighting_compute_site_transfer_matrix(sds, smp->e, ar, ns, smp->e, ns);
        dspmat_sort_entries(T, +2, 0);
        free(ar);
      }
    double frow = ((double)T->ents)/((double)T->rows);
    fprintf(stderr, "site-to-site radiance transfer matrix {T} (%d × %d)", T->rows, T->cols);
    fprintf(stderr, " has %d nonzero elements (%.1f per row)\n", T->ents, frow);
    assert(T->rows == smp->ne);
    assert(T->cols == smp->ne);
  }

void rdo_synthesis_require_basis_transfer_matrix
  ( solid_vec_t *sds,
    finish_vec_t *fns,
    site_vec_t *smp, 
    int channel,
    dspmat_t *T,
    dspmat_t *M,
    dspmat_t *R
  )
  { assert(sds->ne > 0);
    assert(fns->ne > 0); 
    assert(smp->ne > 0); 
    assert(T->ents > 0); 
    assert(T->rows == smp->ne); 
    assert(T->cols == smp->ne); 
    assert(M->ents > 0); 
    assert(M->rows == smp->ne); 
    assert(M->cols == smp->ne); 
    if (R->ents == 0) 
      { fprintf(stderr, "computing the basis-to-basis radiance transfer matrix {R}...\n");
        dspmat_t S = dspmat_new(0,0,0);
        dspmat_mul(T, M, &S);
        rdo_lighting_apply_diffusion_coeffs_to_matrix(&S, smp, sds, fns, channel);
        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. */
        dspmat_inv_mul_GS(M, &S, R, max_iter, omega, abs_tol, rel_tol);
        dspmat_sort_entries(R, +1, 0);
        dspmat_trim(&S, 0);
      }
    double frow = ((double)R->ents)/((double)R->rows);
    fprintf(stderr, "basis-to-basis radiance transfer matrix {R} (%d × %d)", R->rows, R->cols);
    fprintf(stderr, " has %d nonzero elements (%.1f per row)\n", R->ents, frow);
    assert(R->rows == smp->ne);
    assert(R->cols == smp->ne);
  }

void rdo_synthesis_require_rendering_matrix
  ( site_vec_t *pix,
    approx_basis_t *bas,
    int avgElemsPerColumn,
    dspmat_t *Q
  )
  { assert(pix->ne > 0); 
    if (Q->ents == 0) 
      { fprintf(stderr, "computing the basis-to-pixel rendering matrix {Q}...\n");
        (*Q) = rdo_approx_basis_build_matrix(pix, bas, avgElemsPerColumn);
        dspmat_sort_entries(Q, +2, 0);
      }
    double frow = ((double)Q->ents)/((double)Q->rows);
    double fcol = ((double)Q->ents)/((double)Q->cols);
    fprintf(stderr, "rendering matrix {Q} (%d × %d)", Q->rows, Q->cols);
    fprintf(stderr, " has %d nonzero elements (%.1f per row, %.1f per col)\n", Q->ents, frow, fcol);
    assert(Q->rows == pix->ne);
    assert(Q->cols == bas->ne);
  }

site_vec_t rdo_synthesis_compute_pixel_sites
  ( solid_vec_t *sds,
    camera_t *cam, 
    int nh,
    int nv
  )
  {
    int n_pix = nh*nv;
    
    fprintf(stderr, "generating %d × %d = %d pixel sites ...\n ", nh, nv, n_pix);
    
    site_vec_t pix = site_vec_new(n_pix);

    /* Scan the image pixels: */
    int ntot = 0;  /* Counts the number of image pixels generated so far. */
    int ih, iv;    /* Pixel indices (from top left). */
    for(iv = 0; iv < nv; iv++)
      for(ih = 0; ih < nh; ih++)
        {
          /* Compute the scene coordinates {p} of the pixel's center: */
          double h, v; /* Image coords of pixel center. */
          rdo_image_compute_pixel_center_coords(ih, iv, &h, &v);
          double x, y; /* Projection coordinates of the same. */
          rdo_image_coords_to_projection_coords(h, v, nh, nv, &x, &y);
          point3_t p = cam->ctr;
          r3_mix_in(x * cam->radius, &cam->r, &p);
          r3_mix_in(y * cam->radius, &cam->s, &p);

          /* Determine the ray's direction: */
          vector3_t dir = rdo_geometry_pt_pt_dir(&(cam->eye), &p);
          
          /* Trace the ray and save the result (whatever it is): */
          site_t sR;
          double dR;
          rdo_raytrace_first_hit_site(sds, &(cam->eye), &dir, &sR, &dR);
          assert(ntot < pix.ne);
          assert(ntot%nh == ih);
          assert(ntot/nh == iv);
          pix.e[ntot] = sR;
          ntot++; 
        }
    assert(ntot == pix.ne);
    printf("gerados %d sites de pixels\n", ntot);

    return pix;
  }

void rdo_synthesis_identify_highlighted_sites
  ( synthesis_options_t *opSyn,
    solid_vec_t *sds, 
    site_vec_t *smp,
    camera_t *cam,
    image_t *img,
    site_vec_t *pix
  )
  {
    int nh = (int)img->sz[1];
    int nv = (int)img->sz[2];
    int n_des = opSyn->highlight_n;  /* Count of seleted sites. */
    int *des_smp = &(opSyn->highlight_smp[0]);
    int *des_pix = &(opSyn->highlight_pix[0]);
    int jd;
    for (jd = 0; jd < n_des; jd++)
      { double h = opSyn->highlight_hv[jd].c[0];
        double v = opSyn->highlight_hv[jd].c[1];
        /* Project the selected point onto the scene's surface: */
        double x, y; /* Projection plane coordinates of selected image point. */
        rdo_image_coords_to_projection_coords(h, v, nh, nv, &x, &y);
        vector3_t dir;
        rdo_camera_compute_ray_direction(x, y, cam, &dir);
        des_smp[jd] = rdo_raytrace_find_nearest_visible(sds, &cam->eye, &dir, smp);
        des_pix[jd] = rdo_raytrace_find_nearest_visible(sds, &cam->eye, &dir, pix);
        fprintf(stderr, "near image point (%.2f,%.2f):\n", h, v);
        fprintf(stderr, "  sampling site %d", des_smp[jd]);
        if (des_smp[jd] >= 0) { fprintf(stderr, " on solid %d", smp->e[des_smp[jd]].isd); }
        fprintf(stderr, "\n");
        fprintf(stderr, "  pixel site %d = col %d row %d\n", des_pix[jd], des_pix[jd]%nh, des_pix[jd]/nh); 
      }
  }
  
image_type_t rdo_synthesis_options_image_type_parse(argparser_t *pp)
  {
    char *str = argparser_get_next(pp);
    image_type_t tpImage;
    if (strcmp(str, "Diffusion") == 0) 
      { tpImage = image_type_Diffusion; }
    else if (strcmp(str, "Emission") == 0) 
      { tpImage = image_type_Emission; }
    else if (strcmp(str, "Normal") == 0) 
      { tpImage = image_type_Normal; }
    else if (strcmp(str, "Quadrants") == 0) 
      { tpImage = image_type_Quadrants; }
    else if (strcmp(str, "Depth") == 0) 
      { tpImage = image_type_Depth; }
    else if (strcmp(str, "SiteDist") == 0) 
      { tpImage = image_type_SiteDist; }
    else if (strcmp(str, "BasisElems") == 0) 
      { tpImage = image_type_BasisElems; }
    else if (strcmp(str, "QuadrantsApprox") == 0) 
      { tpImage = image_type_QuadrantsApprox; }
    else if (strcmp(str, "Interpolators") == 0) 
      { tpImage = image_type_Interpolators; }
    else if (strcmp(str, "ElemTransfer") == 0) 
      { tpImage = image_type_ElemTransfer; }
    else if (strcmp(str, "Realistic") == 0) 
      { tpImage = image_type_Realistic; }
    else
      { argparser_error(pp, "unknown image type");
        tpImage = image_type_Diffusion;
      }
    return tpImage;
  }

synthesis_options_t *rdo_synthesis_options_parse(argparser_t *pp)
{
  synthesis_options_t *opSyn = NOTNULL(malloc(sizeof(synthesis_options_t)));
  /* ---------------------------------------------------------------------- */
  if (argparser_keyword_present(pp, "-imageType"))
    { opSyn->imageType = rdo_synthesis_options_image_type_parse(pp); }
  else
    { argparser_error(pp, "falta parametro \"-imageType\""); }
  /* ---------------------------------------------------------------------- */
  if (argparser_keyword_present(pp, "-maxBounces"))
    { opSyn->maxBounces = (int)argparser_get_next_int(pp, 0, 1000); }
  else
    { argparser_error(pp, "falta parametro \"-maxBounces\""); }
  /* ---------------------------------------------------------------------- */
  if (argparser_keyword_present(pp, "-crossSize"))
    { opSyn->crossSize = argparser_get_next_double(pp, 0.0, 1000.0); }
  else
    { opSyn->crossSize = 1.0; }
  /* ---------------------------------------------------------------------- */
  if (argparser_keyword_present(pp, "-showSamplingSites"))
    { opSyn->showSamplingSites = TRUE; }
  else
    { opSyn->showSamplingSites = FALSE; }
  /* ---------------------------------------------------------------------- */
  if (argparser_keyword_present(pp, "-signedValues"))
    { opSyn->signedValues = TRUE; }
  else
    { opSyn->signedValues = FALSE; }
  /* ---------------------------------------------------------------------- */
  int n_des = 0;
  while (argparser_keyword_present(pp, "-highlight"))
    { r2_t hv_des;      /* {H,V}oordinates of highlighted pixel. */
      frgb_t color_des; /* Color to use for highlighted pixel. */
      hv_des.c[0] = argparser_get_next_double(pp, -1e200, +1e200);
      hv_des.c[1] = argparser_get_next_double(pp, -1e200, +1e200);
      if (argparser_keyword_present_next(pp, "-color"))
        { color_des = frgb_parse_color(pp); }
      else 
        { color_des = (frgb_t){{ 0.300f, 0.300f, 0.300f }};
          int r = 1 + (n_des % 7);
          if ((r & 1) != 0) { color_des.c[0] = 1.000; }
          if ((r & 2) != 0) { color_des.c[1] = 1.000; }
          if ((r & 4) != 0) { color_des.c[2] = 1.000; }
        }
      if (n_des >= MAX_HIGHLIGHTS) 
        { argparser_error(pp, "too many \"-highlight\" arguments"); break; }
      opSyn->highlight_hv[n_des] = hv_des;
      opSyn->highlight_color[n_des] = color_des;
      n_des++;
    }
  opSyn->highlight_n = n_des;
  /* ---------------------------------------------------------------------- */
  /* !!! completar !!!  */
  return opSyn;
}

void rdo_synthesis_show_site_vec
  ( site_vec_t sts[],
    solid_vec_t *sds,
    camera_t *cam, 
    image_t *img,
    bool_t vis_only, 
    frgb_t *color0, 
    frgb_t *color1, 
    double rad, 
    double hwd,
    bool_t diagonal
  )
  { int nc = (int)img->sz[0];
    int i; 
    for(i = 0; i < sts->ne; i++)
      { site_t *si = &(sts->e[i]);
        int channel;
        for (channel = 0; channel < nc; channel++)
          { float value0 = color0->c[channel]; /* Value if invisible. */
            float value1 = color1->c[channel]; /* Value if visible. */
            rdo_synthesis_show_site(si, sds, cam, img, channel, vis_only, value0, value1, rad, hwd, diagonal);
          }
      }
  } 

void rdo_synthesis_show_selected_sites
  ( site_vec_t *sts,
    int index[],
    frgb_t color[],
    int n,
    solid_vec_t *sds,
    camera_t *cam,
    image_t *img,
    double rad,
    double hwd,
    bool_t diagonal
  )
  {
    int nc = (int)img->sz[0];
    double threshold_Y = 0.250; /* Luminosity threshold for contrast fudging. */
    int jd;
    for (jd = 0; jd < n; jd++)
      { int i = index[jd];
        if ((i >= 0) && (i < sts->ne))
          { site_t *s = &(sts->e[i]);
            frgb_t crj = color[jd];
            /* Fudge the color {crj} so that the cross contrasts with the element: */
            if (frgb_get_Y(&crj) < threshold_Y)
              { /* Mix {crj} with white: */
                frgb_t white = frgb_White;
                crj = frgb_mix(0.5, &crj, 0.5, &white);
              }
            else
              { /* Mix {crj} with black: */
                crj = frgb_scale(0.5, &crj);
              }
            int channel;
            for (channel = 0; channel < nc; channel++)
              { float value1 = crj.c[channel];  /* Value if visible. */
                float value0 = (float)(threshold_Y + crj.c[channel])/2;
                rdo_synthesis_show_site(s, sds, cam, img, channel, FALSE, value0, value1, rad, hwd, diagonal);
              }
          }
      }
  }

void rdo_synthesis_show_site
  ( site_t *s, 
    solid_vec_t *sds,
    camera_t *cam, 
    image_t *img, 
    int channel,
    bool_t vis_only, 
    float value0, 
    float value1, 
    double rad, 
    double hwd,
    bool_t diagonal
  )
  { if (s == NULL) { return; }
    int nh = (int)img->sz[1];
    int nv = (int)img->sz[2];
    /* Get the pixel indices {ih,iv} on the image: */
    double x, y; /* Projection plane coordinates of point's projection. */
    rdo_camera_project_point(&s->position, cam, &x, &y);
    double h, v; /* Image coords of point's projection. */
    rdo_image_coords_from_projection_coords(x, y, nh, nv, &h, &v);
    
    /* Determine the site's visibility from the camera: */
    double fvis = rdo_raytrace_pt_site_visibility(sds, &cam->eye, &s->position, &s->normal);
    
    if ((! vis_only) || (fvis > 0))
      { /* Paint the site: */
        float val = (float)((1 - fvis) * value0 + fvis * value1);
        (void)float_image_paint_cross(img, channel, h, nv - v, rad, hwd, diagonal, val, 3);
      }
  } 

void rdo_synthesis_compute_pixel_radiances_from_coeffs
  ( double L[],
    int nL,
    dspmat_t *Q, 
    double F[],
    int nF
  )
  {
    dspmat_map_col(Q, L, nL, F, nF);
  }

void rdo_synthesis_paint_pixel_radiances
  ( image_t *img, 
    int channel, 
    double F[],
    int nF
  )
  { int nh = (int)img->sz[1];
    int nv = (int)img->sz[2];
    int n_pix = nh*nv;
    assert(nF == n_pix);

    int ih, iv;
    for(iv = 0; iv < nv; iv++)
      { for(ih = 0; ih < nh; ih++)
          { int k = iv*nh + ih;
            if ((debug_synthesis) && ((k < 20) || (k % 100 == 0)))
              { fprintf(stderr, " F[%d] = %9.5f\n", k, F[k]); }
            float_image_set_sample(img, channel, ih, nv - 1 - iv, (float)F[k]);
          }
      }
  }

void rdo_synthesis_compute_image
  ( scene_t *scn, 
    camera_t *cam,
    dist_type_t tpDist,
    basis_elem_type_t tpElem,
    sampling_options_t *opSmp,
    approx_basis_options_t *opBas,
    synthesis_options_t *opSyn,
    image_t *img
  )
  {
    int nc = (int)img->sz[0];
    int nh = (int)img->sz[1];
    int nv = (int)img->sz[2];
    
    /* Determines the sites corresponding to the image pixels: */
    site_vec_t pix = rdo_synthesis_compute_pixel_sites(&(scn->sds), cam, nh, nv);
    int n_pix = pix.ne;    
    assert(n_pix == nh*nv); /* Assumes that every pixel has its site (possibly null). */

    /* Vector {F} holds the radiances at the image pixel sites: */
    double *F = NOTNULL(malloc(n_pix * sizeof(double)));

    /* Matrix {Q} is the rendering matrix, that maps {L} to {F} (we may not need it):  */
    dspmat_t Q = dspmat_new(0,0,0);
    
    /* Get hold of the scene's {T} matrix (may be empty, and we may not need it): */
    dspmat_t *T = &(scn->T);
    
    /* Get hold of the scene's {M} matrix (may be empty, and we may not need it): */
    dspmat_t *M = &(scn->M);

    /* Loop over the color channels: */
    int channel;
    for (channel = 0; channel < nc; channel++)
      { 

        /* Dispatch according to the desired image type: */
        switch (opSyn->imageType)
          {
            case image_type_Diffusion:
            case image_type_Emission:
            case image_type_Normal:
            case image_type_Quadrants:
            case image_type_Depth:
            case image_type_SiteDist:
              { if (opSyn->showSamplingSites)
                  { /* Make sure that the sampling sites have been generated: */
                    rdo_synthesis_require_sampling_sites(&scn->sds, &(cam->eye), tpDist, opSmp, &scn->smp);
                  }
                rdo_synthesis_identify_highlighted_sites(opSyn, &scn->sds, &scn->smp, cam, img, &pix);
                rdo_synthesis_plain_compute_image
                  ( &(scn->sds), &(scn->fns), cam, 
                    tpDist, opSyn, 
                    channel, &pix, F
                  );
                break;
              }

            case image_type_BasisElems:
            case image_type_Interpolators:
            case image_type_QuadrantsApprox:
            case image_type_ElemTransfer:
            case image_type_Realistic:
              { rdo_synthesis_require_sampling_sites(&scn->sds, &(cam->eye), tpDist, opSmp, &scn->smp);
                rdo_synthesis_identify_highlighted_sites(opSyn, &scn->sds, &scn->smp, cam, img, &pix);
                rdo_synthesis_basis_compute_image
                  ( scn, &(cam->eye), 
                    tpDist, tpElem, opSmp, opBas, 
                    M, T,
                    opSyn, channel, &pix, &Q, F
                  );
                break;
              }
          }

        rdo_synthesis_paint_pixel_radiances(img, channel, F, n_pix); 

        if (opSyn->showSamplingSites)
          { /* Paints crosses at the sampling sites: */
            frgb_t color0 = (frgb_t){{ 0.500f, 0.400f, 0.200f }};
            frgb_t color1 = (frgb_t){{ 1.000f, 0.800f, 0.500f }};
            double rad = opSyn->crossSize; /* Radius of crosses (pixels). */
            double hwd = 0.5;                 /* Half-widths of cross arms (pixels). */
            rdo_synthesis_show_site_vec
              ( &scn->smp, &scn->sds, cam, img, FALSE, &color0, &color1, rad, hwd, FALSE ); 
          }

        /* Show the highlighted sites: */
        { double rad = 2.0 * opSyn->crossSize; /* Radius of crosses (pixels). */
          double hwd = 1.0;                       /* Half-widths of cross arms (pixels). */
          rdo_synthesis_show_selected_sites
            ( &scn->smp, 
              opSyn->highlight_smp, opSyn->highlight_color, opSyn->highlight_n, 
              &scn->sds, cam, img,
              rad, hwd, TRUE
            );
        }
      }
                
    if (opSyn->signedValues) { rdo_image_map_zero_to_half(img); }

    /* Discard work areas (but keep {scn->T}, {scn->M}, possibly for output): */
    site_vec_trim(&pix,0);
    free(F);
    dspmat_trim(&Q, 0);
  }