/* See {rdo_site.h} */
#define rdo_site_C_COPYRIGHT "Copyright © 2008 Danillo Pereira and J. Stolfi, UNICAMP"
/* Last edited on 2024-12-21 11:51:30 by stolfi */ 

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

#include <frgb.h>
#include <frgb_ops.h>
#include <filefmt.h>
#include <jsrandom.h>
#include <nget.h>
#include <fget.h>
#include <fget_geo.h>
#include <vec.h>
#include <affirm.h>
#include <bool.h>

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

vec_typeimpl(site_vec_t, site_vec, site_t);

void rdo_site_from_point(solid_t *sd, int isd, point3_t *position, site_t *sP)
  { sP->position = (*position);
    rdo_solid_get_normal(sd, position, &(sP->normal));
    sP->isd = isd;
  }

void rdo_site_vec_scramble(site_vec_t *sts)
  {
    site_t *s = sts->e;
    int n = sts->ne;
    int i;
    for(i = 0; i < n; i++)
      { /* Swaps {s[i]} with {s[j]} where {j} is random in {i..n-1}: */
        int j = int32_abrandom(i, n - 1);
        if (j != i) { site_t t = s[j]; s[j] = s[i]; s[i] = t; }
      } 
  }

void rdo_site_find_closest_sites
  ( site_t *sr,
    site_vec_t *sts, 
    int max_nc, 
    dist_type_t tpDist, 
    bool_t same_solid,
    int ind[], 
    double dst[], 
    int *ncP
  )
  {
    int nc = 0; /* Number of candidates in {ind[],dst[]}. */
    
    auto void insert_in_queue(int ind_c, double dst_c);
    /* Inserts {sts->e[ind_c]} with distance {dst_c} in the list of the
      {max_nc} best candidates, namely {ind[0..nc-1],dst[0..nc-1]}, in
      order of increasing distance. If {nc < max_nc}, increments {nc}; else
      discards {ind_c} or the most distant of the current candidates,
      preserving {nc == max_nc}. */

    void insert_in_queue(int ind_c, double dst_c)
      { /* Look for the position {j} in the queue where to insert the candidate: */
        int j = nc;
        while ((j > 0) && (dst_c < dst[j-1]))
          { /* Push candidate {j-1} to position {j}; */
            /* If position {j} does not exist, discard the candidate {j-1}: */
            if (j < max_nc)
              { ind[j] = ind[j-1]; dst[j] = dst[j-1]; }
            /* The empty slot now is {j-1}: */
            j--;
          }
        if (j < max_nc) 
          { /* Insert the new candidate {ind_c,dst_c} in position {j}:  */
            ind[j] = ind_c; dst[j] = dst_c;
          }
        /* If the queue is not full yet, we got one more candidate: */
        if (nc < max_nc) { nc++; }
      }

    /* Scan all sites and push candidates into the queue: */
    int i;
    for(i = 0; i < sts->ne; i++)
      { site_t *sti = &(sts->e[i]);
        if ((! same_solid) || ((sr->isd >= 0) && (sr->isd == sti->isd)))
          { double d = rdo_site_distance(sr, sti, tpDist);
            if ((d < +INF) && (d != 0.0)) { insert_in_queue(i, d); }
          }
      }
    
    (*ncP) = nc;
  }

double rdo_site_estimate_effective_area(site_t *sr, site_vec_t *sts)
  {
    int max_clos = 7; /* Number of closest sites to look for. */
    
    int n_clos = 0; /* Number of closest sites actually found ({max_clos} or less). */
    
    /* The {n_clos} sites closest to {st[ind]}, excluding itself:  */
    int ind_clos[max_clos];    /* Indices of those sites. */
    double dst_clos[max_clos]; /* Their distances to {st[ind]} in increasing order. */

    dist_type_t tpDist = dist_type_Euclidean; /* Ignore normals (!!! confirm) */
    
    if (sr->isd >= 0) 
      { /* Look first for sites on the same solid object: */
        rdo_site_find_closest_sites(sr, sts, max_clos, tpDist, TRUE, ind_clos, dst_clos, &n_clos);
      }

    if (n_clos == 0)
      { /* No other sites on the same solid object; try again with all sites: */
        rdo_site_find_closest_sites(sr, sts, max_clos, tpDist, FALSE, ind_clos, dst_clos, &n_clos);
      }
    
    if (n_clos == 0)
      { /* No clue about {sr}'s effective area: */
        return NAN;
      }
    else 
      { /* Estimate the effective area from the distance to the {n_clos}th neighbor: */
        double R = dst_clos[n_clos-1];
        double aR = M_PI*R*R;
        /* Assume that sites are uniformly distributed on a mostly flat surface: */
        return aR/(n_clos - 0.5);
      }
  }  

double *rdo_site_estimate_all_effective_areas(site_vec_t *sts)
  {
    int n = sts->ne;
    double *ar = NOTNULL(malloc(n * sizeof(double)));
    int i;
    for (i = 0; i < n; i++) 
      { /* Find the radius {ri} that contains {nin} other sites: */
        ar[i] = rdo_site_estimate_effective_area(&(sts->e[i]), sts);
        demand(isfinite(ar[i]), "could not determine the site's effective area");
      }
    return ar;
  }  

void rdo_site_remove_closest(site_vec_t *sts, double minDist, dist_type_t td)
  {
    /* Greedy randomized algorithm. */
    /* !!! custo quadrático! deveria usar algum type de bucketing... !!! */

    int nin = sts->ne;  /* Count of original sites. */
    site_t *s = sts->e; /* Original sites. */
    
    int nsb = 0; /* Count of surviving sites. */

    /* Scan all sites: */
    int i;
    for(i = 0; i < nin; i++)
      { /* At this moment, the chosen sites are {s[0..nsb-1]}. */
        /* At this moment, sites {s[i..nin-1]} are still to be considered. */

        /* Checks whether site {i} can be selected: */
        bool_t i_OK = TRUE; /* In principle... */
        int j;
        for(j = 0; i_OK && (j < nsb); j++)
          { double dist = rdo_site_distance(&(s[i]), &(s[j]), td);
            if (dist < minDist) { i_OK = FALSE; }
          }

        if (i_OK)
          { /* Site {i} lies far enough from other sites, pick it: */
            s[nsb] = s[i]; 
            nsb++;
          }
      }

    /* keep only the selected sites: */
    site_vec_trim(sts, nsb);
    fprintf(stderr, "removing close sites:  %d input sites  %d output sites\n", nin, sts->ne);
  }

double rdo_site_distance(site_t *sA, site_t *sB, dist_type_t td)
  {
    if (td == dist_type_Euclidean)
      { 
        double disitioAB = r3_dist(&(sA->position), &(sB->position));
        return disitioAB;
      }
    else if (td == dist_type_Directional)
      { 
        double cosAB = r3_dot(&(sA->normal), &(sB->normal));   
        if (cosAB <= 0.0)
          { return +INF; }
        else
          { double distAB = r3_dist(&(sA->position), &(sB->position));
            //!!! se {distAB == 0} mas {cosAB < 1}, deveria ser {> 0} !!!
            return distAB/cosAB;
          }
      }
    else
      { assert(FALSE); return NAN; }
  }


#define site_vec_FILE_VERSION "2008-07-14"
  /* Identifies the format of a site list file. */

void rdo_site_vec_write(FILE *wr, site_vec_t *sts)
  {
    /* Write the header line: */
    filefmt_write_header(wr, "site_vec_t", site_vec_FILE_VERSION);
    /* Here whould write comments provided by client. */

    /* Write the site count: */
    fprintf(wr, "num_sitios = %d\n", sts->ne);
    
    /* Write the sites: */
    int i;
    for(i = 0; i < sts->ne; i++)
      { rdo_site_write(wr, &(sts->e[i]));
        fprintf(wr, "\n");    
      }
    
    /* Write the footer line, and synchronize: */
    filefmt_write_footer(wr, "site_vec_t");
    fflush(wr);
  }

site_vec_t rdo_site_vec_read(FILE *rd)
  {
    /* Read the header line and comment lines (if any): */
    filefmt_read_header(rd, "site_vec_t", site_vec_FILE_VERSION);
    char *cmt = filefmt_read_comment(rd, '#');
    free(cmt); /* For now. */

    /* Read the number of sites */
    int nob = nget_int32(rd, "num_sitios"); fget_eol(rd);
    
    /* Read the sites: */
    site_vec_t sts = site_vec_new(nob);
    int i;
    for(i = 0; i < nob; i++)
      { sts.e[i] = rdo_site_read(rd);
        fget_eol(rd);    
      }
    
    /* Read the footer line: */
    filefmt_read_footer(rd, "site_vec_t");

    return sts;
  }

void rdo_site_write(FILE *wr, site_t *st)
  { r3_gen_print(wr, &(st->position), "%9.3f", " ( ", " ", " )");
    r3_gen_print(wr, &(st->normal), "%9.3f", " ( ", " ", " )");
    fprintf(wr, " %5d", st->isd);
  }

site_t rdo_site_read(FILE *rd)
  { site_t st;
    
    fget_skip_spaces_and_match(rd, "(");
    st.position = fget_r3(rd);
    fget_skip_spaces_and_match(rd, ")");
    
    fget_skip_spaces_and_match(rd, "(");
    st.normal = fget_r3(rd);
    fget_skip_spaces_and_match(rd, ")");
    
    st.isd = fget_int32(rd);
    
    return st;
  }