star-4v_s

An invariant property of diffuse random walks
git clone git://git.meso-star.fr/star-4v_s.git
Log | Files | Refs | README | LICENSE

s4vs_realization.c (3924B)


      1 /* Copyright (C) 2015-2018, 2021, 2024 |Méso|Star> (contact@meso-star.com)
      2  *
      3  * This program is free software: you can redistribute it and/or modify
      4  * it under the terms of the GNU General Public License as published by
      5  * the Free Software Foundation, either version 3 of the License, or
      6  * (at your option) any later version.
      7  *
      8  * This program is distributed in the hope that it will be useful,
      9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     11  * GNU General Public License for more details.
     12  *
     13  * You should have received a copy of the GNU General Public License
     14  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     15 
     16 #include <rsys/float3.h>
     17 
     18 #include <star/s3d.h>
     19 #include <star/ssp.h>
     20 #include <star/smc.h>
     21 
     22 #include "s4vs_realization.h"
     23 
     24 /*******************************************************************************
     25  * Helper function
     26  ******************************************************************************/
     27 int
     28 s4vs_discard_self_hit
     29   (const struct s3d_hit* hit,
     30    const float ray_org[3],
     31    const float ray_dir[3],
     32    const float ray_range[2],
     33    void* ray_data,
     34    void* filter_data)
     35 {
     36   const struct s3d_primitive* prim_from = ray_data;
     37 
     38   /* Avoid unused variable warn */
     39   (void)ray_org, (void)ray_dir, (void)ray_range, (void)filter_data;
     40   return prim_from ? S3D_PRIMITIVE_EQ(prim_from, &hit->prim) : 0;
     41 }
     42 
     43 /*******************************************************************************
     44  * 4V/S integrand
     45  ******************************************************************************/
     46 res_T
     47 s4vs_realization
     48   (void* out_length,
     49    struct ssp_rng* rng,
     50    const unsigned ithread,
     51    const uint64_t irealisation,
     52    void* context)
     53 {
     54   struct s4vs_context* ctx = (struct s4vs_context*)context;
     55   struct s3d_attrib attrib;
     56   struct s3d_primitive prim;
     57   double sample[3];
     58   double normal[3];
     59   float u[3], x[3], st[2];
     60   const float range[2] = {0.f, FLT_MAX};
     61   struct s3d_hit hit;
     62   double w = 0;
     63   double sigma = 0;
     64   int keep_running = 0;
     65   float r0, r1, r2;
     66   (void)ithread, (void)irealisation; /* Avoid "unused variable" warning */
     67 
     68   /* Sample a surface location, i.e. primitive ID and parametric coordinates */
     69   r0 = ssp_rng_canonical_float(rng);
     70   r1 = ssp_rng_canonical_float(rng);
     71   r2 = ssp_rng_canonical_float(rng);
     72   S3D(scene_view_sample(ctx->view, r0, r1, r2, &prim, st));
     73 
     74   /* retrieve the sampled geometric normal and position */
     75   S3D(primitive_get_attrib(&prim, S3D_GEOMETRY_NORMAL, st, &attrib));
     76   d3_normalize(normal, d3_set_f3(normal, attrib.value));
     77   S3D(primitive_get_attrib(&prim, S3D_POSITION, st, &attrib));
     78   f3_set(x, attrib.value);
     79 
     80   /* Cosine weighted sampling of the hemisphere around the sampled normal */
     81   ssp_ran_hemisphere_cos(rng, normal, sample, NULL);
     82   f3_set_d3(u, sample);
     83 
     84   /* Find the 1st hit from the sampled location along the sampled direction */
     85   S3D(scene_view_trace_ray(ctx->view, x, u, range, &prim, &hit));
     86 
     87   /* No intersection <=> numerical imprecision or geometry leakage */
     88   if(S3D_HIT_NONE(&hit)) return RES_UNKNOWN_ERR;
     89 
     90   keep_running = 1;
     91   while(keep_running) { /* Here we go for the diffuse random walk */
     92 
     93     /* Sample a length according to ks */
     94     sigma = ssp_ran_exp(rng, ctx->ks);
     95 
     96     if(sigma < hit.distance) {
     97       int i;
     98       FOR_EACH(i, 0, 3) x[i] = x[i] + (float)sigma*u[i];
     99       d3_normalize(sample, d3_set_f3(sample, u));
    100       f3_set_d3(u, ssp_ran_sphere_hg(rng, sample, ctx->g, sample, NULL));
    101 
    102       /* sample a new direction */
    103       S3D(scene_view_trace_ray(ctx->view, x, u, range, NULL, &hit));
    104 
    105       w = w + sigma;
    106 
    107       /* No intersection <=> numerical imprecision or geometry leakage */
    108       if(S3D_HIT_NONE(&hit)) return RES_UNKNOWN_ERR;
    109 
    110     } else { /* Stop the random walk */
    111       w = w + hit.distance;
    112       keep_running = 0;
    113     }
    114   }
    115 
    116   SMC_DOUBLE(out_length) = w;
    117   return RES_OK;
    118 }