/* See {rdo_sampling.h} */

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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#include <frgb.h>
#include <argparser.h>
#include <jsrandom.h>

#include <rdo_geometry.h>
#include <rdo_solid.h>
#include <rdo_raytrace.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>
  
site_vec_t rdo_sampling_choose_approx_basis_centers
  ( solid_vec_t *sds,
    dist_type_t tpDist,
    sampling_options_t *opSmp, 
    point3_t *eye 
  )
  {
    fprintf(stderr, "generating sampling sites:");
    fprintf(stderr, "  initial sites = %d\n", opSmp->numInitialSites);
    fprintf(stderr, "  max sites = %d\n", opSmp->maxNumSites);
    
    site_vec_t sts = site_vec_new(0);
    
    /* Get the scene's bounding box: */
    point3_t corner[2];
    rdo_solid_vec_bounding_box(sds, corner);
    
    int n_ini = opSmp->numInitialSites; /* Total number of initial samples. */
    
    /* Divide the initial samples by the three methods: */
    double w_eye = (opSmp->sampleFromObserver ? 1.0 : 0.0);
    double w_asd = (opSmp->sampleFromAllSides ? 4.0 : 0.0);
    double w_lts = (opSmp->sampleFromLights ? 5.0 : 0.0);
    double w_tot = w_eye + w_asd + w_lts;
    assert(w_tot > 0);
    int n_eye = (int)floor(n_ini * w_eye/w_tot);
    int n_asd = (int)floor(n_ini * w_asd/w_tot);
    int n_lts = n_ini - n_eye - n_asd;

    if (opSmp->sampleFromObserver) 
      { rdo_sampling_generate_sites_by_projection(sds, corner, eye, n_eye, &sts); }

    if (opSmp->sampleFromAllSides)
      { 
        /* Get the center {mid} of the scene's box, and its radius vector {rad}: */
        point3_t mid; r3_mix(+0.5, &(corner[0]), +0.5, &(corner[1]), &mid);
        point3_t rad; r3_mix(-0.5, &(corner[0]), +0.5, &(corner[1]), &rad);

        /* !!! Should generate random inward rays from the surface of the bbox !!! */
        /* Pick four points {f1,f2,f3,f4} outside and above the scene: */
        double dist = 4 * r3_norm(&rad);
        point3_t f1 = (point3_t){{ mid.c[0] - dist, mid.c[1] - dist, mid.c[2] + dist }};
        point3_t f2 = (point3_t){{ mid.c[0] - dist, mid.c[1] + dist, mid.c[2] + dist }};
        point3_t f3 = (point3_t){{ mid.c[0] + dist, mid.c[1] - dist, mid.c[2] + dist }};
        point3_t f4 = (point3_t){{ mid.c[0] + dist, mid.c[1] + dist, mid.c[2] + dist }};

        rdo_sampling_generate_sites_by_projection(sds, corner, &f1, n_asd/4, &sts);
        rdo_sampling_generate_sites_by_projection(sds, corner, &f2, n_asd/4, &sts);
        rdo_sampling_generate_sites_by_projection(sds, corner, &f3, n_asd/4, &sts);
        rdo_sampling_generate_sites_by_projection(sds, corner, &f4, n_asd/4, &sts); 
      }

    if (opSmp->sampleFromLights)
      { 
        fprintf(stderr, "option \"-sampleFromLights\" (%d) not implemented yet\n", n_lts); 
      }
    
    rdo_site_remove_closest(&sts, opSmp->minDist, tpDist);
    
    /* Generate some secondary sites: */
    int j;
    for(j = 0; j < opSmp->secondarySamples_numGens; j++) 
      { int old_ne = sts.ne;
        rdo_sampling_generate_secondary_sites(sds, &sts, opSmp->secondarySamples_numRays);
        rdo_site_remove_closest(&sts, opSmp->minDist, tpDist);
        fprintf(stderr, "generation %d:", j);
        fprintf(stderr, "  %d additional sites", sts.ne - old_ne);
        fprintf(stderr, "  %d total sites\n", sts.ne);
      } 
      
    return sts;
  }

void rdo_sampling_generate_sites_by_projection
  ( solid_vec_t *sds, 
    point3_t corner[],
    point3_t *eye, 
    int n,
    site_vec_t *sts
  )
  {
    fprintf(stderr, "generating %d additional sites by projection ...\n", n);
    
    int nold = sts->ne;  /* Number of original sites in {sts} */
    int ntot = nold;     /* Currente number of sites in {sts}. */
    int k;
    for(k = 0; k < n; k++)
      { /* Pick a point {aim} on the scene bounding box: */
        point3_t aim;
        int jface = int32_abrandom(0,2);  /* Axis perpendicular to the chosen face. */
        int j;
        for (j = 0; j < 3; j++)
          { if (j == jface)
              { /* Pick one of the two faces: */
                aim.c[j] = corner[int32_abrandom(0,1)].c[j];
              }
            else
              { /* Pick a coordinate at random in the bounding range: */
                double s = drandom(), t = 1 - s;
                aim.c[j] = s*corner[0].c[j] + t*corner[1].c[j];
              }
          }

        /* Shoot the ray: */
        vector3_t dir = rdo_geometry_pt_pt_dir(eye, &aim);
        site_t sR;
        double dR;
        rdo_raytrace_first_hit_site(sds, eye, &dir, &sR, &dR);
        if (isfinite(dR) && (sR.isd >= 0))
          { site_vec_expand(sts, ntot);
            sts->e[ntot] = sR;
            ntot++; 
          } 
      }
    fprintf(stderr, "generated %d new sites (%d total sites)\n", ntot - nold, ntot);
    site_vec_trim(sts, ntot);
  }

void rdo_sampling_generate_secondary_sites(solid_vec_t *sds, site_vec_t *sts, int num_rays)
  { 
    fprintf(stderr, "generating %d secondary sites from each of %d sites ...\n", num_rays, sts->ne);
    
    /* Try to generate {num_rays} child rays from each existing site: */
    int nold = sts->ne;  /* Original count of sites in {sts} */
    int ntot = nold;     /* Total count of sites in {sts}. */
    int i; 
    for(i = 0; i < nold; i++) 
      { int k;
        for (k = 0; k < num_rays; k++)
          { site_t sR = rdo_sampling_generate_secondary_site(sds, &(sts->e[i]));
            if (sR.isd >= 0)
              { site_vec_expand(sts, ntot);
                sts->e[ntot] = sR;
                ntot++; 
              }
          }
      }
      
    fprintf(stderr, "generated %d new derived sites (%d total sites)\n", ntot - nold, ntot);
    site_vec_trim(sts, ntot);
  }

site_t rdo_sampling_generate_secondary_site(solid_vec_t *sds, site_t *sO)
  {
    /* Perturb the origin site along the normal: */
    point3_t origem;
    r3_mix(1.0, &(sO->position), RAY_ORIGIN_EPSILON, &(sO->normal), &origem);

    /* Generate a random direction in the positive hemisphere of {sO}, */
    /* with distribution biased towards the normal of {sO}: */
    vector3_t dir;
    r3_throw_dir(&dir);
    r3_add(&(sO->normal), &dir, &dir);
    r3_mix_in(1.0e-12, &(sO->normal), &dir); /* Paranoia to avoi null {dir}. */
    r3_dir(&dir, &dir);
    
    assert(r3_dot(&(sO->normal), &dir) >= 0.0);
    
    /* Cast a ray from {sO} along {dir}:  */
    site_t sR;
    double dR;
    rdo_raytrace_first_hit_site(sds, &origem, &dir, &sR, &dR);

    if ((! isfinite(dR)) || (r3_dist(&(sR.position), &(sO->position)) <= VISIBILITY_TOLERANCE))
      { /* Site at infinity or self-intersection, discard: */
        sR.isd = -1;
        sR.position = (point3_t){{ 0,0,0 }};
        sR.normal = (vector3_t){{ 0,0,0 }};
      }
    
    return sR;
  }

  
#define MAX_MAX_SITES 100000
  /* Max number of sampling sites that the user may ask for. */

sampling_options_t *rdo_sampling_options_parse(argparser_t *pp)
  {
    sampling_options_t *opSmp = NOTNULL(malloc(sizeof(sampling_options_t)));

    argparser_get_keyword(pp, "-maxNumSites");
    opSmp->maxNumSites = (int)argparser_get_next_int(pp, 1, MAX_MAX_SITES);
    
    argparser_get_keyword(pp, "-minDist");
    opSmp->minDist = argparser_get_next_double(pp, 1.0e-10, 1.0e+10);
    
    if (argparser_keyword_present(pp, "-numInitialSites"))
      { opSmp->numInitialSites = (int)argparser_get_next_int(pp, 1, MAX_MAX_SITES); }
    else 
      { opSmp->numInitialSites = opSmp->maxNumSites; }
      
    opSmp->sampleFromObserver = argparser_keyword_present(pp, "-sampleFromObserver");
    opSmp->sampleFromAllSides = argparser_keyword_present(pp, "-sampleFromAllSides");
    opSmp->sampleFromLights = argparser_keyword_present(pp, "-sampleFromLights");
    if (! (opSmp->sampleFromObserver || opSmp->sampleFromAllSides || opSmp->sampleFromLights))
      { opSmp->sampleFromObserver = TRUE;
        opSmp->sampleFromAllSides = TRUE;
        opSmp->sampleFromLights = TRUE;
      }

    if (argparser_keyword_present(pp, "-secondarySamples"))
      { opSmp->secondarySamples_numGens = (int)argparser_get_next_int(pp, 1, MAX_MAX_SITES);
        opSmp->secondarySamples_numRays = (int)argparser_get_next_int(pp, 1, MAX_MAX_SITES);
      }
    else 
      { opSmp->secondarySamples_numGens = 0;
        opSmp->secondarySamples_numRays = 0;
      }

    return opSmp;
  }