htsky

commit 2db0b59ceb6a0f3f395438e80efaeccd6009fb97
parent bbc29781f3007b7927f854b1a898ee8e73371a41
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu, 16 Apr 2020 21:45:00 +0200

Handle a possible numerical issue in htsky_fetch_raw_property

The computed voxel coordinates might be greater than the grid resolution.

Diffstat:
M
src/htsky.c
 | 
13
+++++++++++++

1 file changed, 13 insertions(+), 0 deletions(-)
diff --git a/src/htsky.c b/src/htsky.c
@@ -772,6 +772,8 @@ htsky_fetch_raw_property
     /* Not in atmopshere => No gas */
     if(!in_atmosphere) comp_mask &= ~HTSKY_CPNT_FLAG_GAS;
   } else {
+    const double* upp;
+    const size_t* def;
     world_to_cloud(sky, pos, pos_cs);
 
     /* Compute the index of the voxel to fetch */
@@ -789,6 +791,17 @@ htsky_fetch_raw_property
         ((pos_cs[2] - cloud_desc->lower[2]) / sky->lut_cell_sz);
       ivox[2] = splits[ilut].index + (pos_cs[2] > splits[ilut].height);
     }
+
+    /* Handle numerical issues that may lead to a position lying onto the cloud
+     * upper boundaries */
+    def = sky->htcp_desc.spatial_definition;
+    upp = cloud_desc->upper;
+    if(ivox[0] == def[0] && eq_eps(pos_cs[0], upp[0], 1.e-6)) ivox[0] = def[0];
+    if(ivox[1] == def[1] && eq_eps(pos_cs[1], upp[1], 1.e-6)) ivox[1] = def[1];
+    if(ivox[2] == def[2] && eq_eps(pos_cs[2], upp[2], 1.e-6)) ivox[2] = def[2];
+    if(ivox[0] >= def[0]) FATAL("Out of bound X voxel coordinate\n");
+    if(ivox[1] >= def[1]) FATAL("Out of bound Y voxel coordinate\n");
+    if(ivox[2] >= def[2]) FATAL("Out of bound Z voxel coordinate\n");
   }
 
   if(comp_mask & HTSKY_CPNT_FLAG_PARTICLES) {