htsky

commit e218bbb4a3238af68c950167b4f07a1319d948bf
parent bbf347fa822717ae0e40b680f2bbf262dc772176
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu,  9 Nov 2023 09:53:21 +0100

Fix invalid memory read for clouds with irregular Z dimensions

For clouds with irregular Z dimensions, a look-up table is used to speed
up the calculation of the voxel index along the Z dimension
corresponding to a given position. This queried position can lie on a
bound of a LUT's cell. In such case, one should define whether the bound
is inclusive or exclusive, i.e. if the corresponding cell must be
considered or not. Actually, the boundary is inclusive on the lower
bound of the cell and exclusive on its upper bound. And there was an
issue on this : an exclusive boundary was considered inclusive. Although
this bug had no impact on the accuracy of calculations, it resulted in
an invalid memory read when accessing the cell considered erroneously.

Diffstat:
M
src/htsky.c
 | 
4
++++
M
src/htsky_cloud.c
 | 
47
++++++++++++++++++++++++++++++++++++++---------

2 files changed, 42 insertions(+), 9 deletions(-)
diff --git a/src/htsky.c b/src/htsky.c
@@ -829,8 +829,10 @@ htsky_fetch_raw_property
        * position in the svx2htcp_z Look Up Table and use the LUT to define the
        * voxel index into the HTCP descriptor */
       const struct split* splits = darray_split_cdata_get(&sky->svx2htcp_z);
+      const size_t nsplits = darray_split_size_get(&sky->svx2htcp_z);
       const size_t ilut = (size_t)
         ((pos_cs[2] - cloud_desc->lower[2]) / sky->lut_cell_sz);
+      if(ilut >= nsplits) FATAL("LUT index out of range\n");
       ivox[2] = splits[ilut].index + (pos_cs[2] > splits[ilut].height);
     }
 
@@ -935,7 +937,9 @@ htsky_fetch_temperature(const struct htsky* sky, const double pos[3])
        * position in the svx2htcp_z Look Up Table and use the LUT to define the
        * voxel index into the HTCP descripptor */
       const struct split* splits = darray_split_cdata_get(&sky->svx2htcp_z);
+      const size_t nsplits = darray_split_size_get(&sky->svx2htcp_z);
       const size_t ilut = (size_t)((pos_cs[2] - desc->lower[2]) / sky->lut_cell_sz);
+      if(ilut >= nsplits) FATAL("LUT index out of range\n");
       ivox[2] = splits[ilut].index + (pos_cs[2] > splits[ilut].height);
     }
 
diff --git a/src/htsky_cloud.c b/src/htsky_cloud.c
@@ -214,6 +214,7 @@ compute_k_bounds_irregular_z
   double Ca_avg;
   double Cs_avg;
   double pos_z;
+  double lut_low, lut_upp;
   size_t ilut_low, ilut_upp;
   size_t ilut;
   size_t ivox_z_prev = SIZE_MAX;
@@ -241,8 +242,20 @@ compute_k_bounds_irregular_z
 
   /* Compute the *inclusive* bounds of the indices of the LUT cells
    * overlapped by the SVX voxel */
-  ilut_low = (size_t)((vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz);
-  ilut_upp = (size_t)((vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz);
+  lut_low = (vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz;
+  lut_upp = (vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz;
+  ilut_low = (size_t)lut_low;
+  ilut_upp = (size_t)lut_upp;
+
+  /* A LUT coordinate equal to its search index means that this coordinate lies
+   * strictly on the boundary of a LUT cell. If this coordinate is the upper
+   * limit of a coordinate range, we subtract one from the LUT index to ensure
+   * that the corresponding cell is indeed included in the submitted range
+   * coordinates */
+  if(lut_upp == ilut_upp) {
+    --ilut_upp;
+  }
+
   ASSERT(ilut_low <= ilut_upp);
 
   /* Prepare the iteration over the LES voxels overlapped by the SVX voxel */
@@ -250,8 +263,10 @@ compute_k_bounds_irregular_z
   ka_bounds[1] = ks_bounds[1] = kext_bounds[1] =-DBL_MAX;
   ivox_z_prev = SIZE_MAX;
   pos_z = vox_svx_low[2];
-  ASSERT(pos_z >= (double)ilut_low * sky->lut_cell_sz);
-  ASSERT(pos_z <= (double)ilut_upp * sky->lut_cell_sz);
+  ASSERT(vox_svx_low[2] >=
+    sky->htcp_desc.lower[2] + sky->lut_cell_sz*(double)(ilut_low));
+  ASSERT(vox_svx_upp[2] <=
+    sky->htcp_desc.lower[2] + sky->lut_cell_sz*(double)(ilut_upp+1));
 
   /* Iterate over the LUT cells overlapped by the voxel */
   FOR_EACH(ilut, ilut_low, ilut_upp+1) {
@@ -414,6 +429,7 @@ compute_xh2o_range_irregular_z
   size_t ilut_low, ilut_upp;
   size_t ilut;
   size_t ivox_z_prev;
+  double lut_low, lut_upp;
   double low[2], upp[2];
   double pos_z;
   double ipart;
@@ -435,17 +451,30 @@ compute_xh2o_range_irregular_z
 
   /* Compute the *inclusive* bounds of the indices of the LUT cells
    * overlapped by the SVX voxel */
-  ilut_low = (size_t)((vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz);
-  ilut_upp = (size_t)((vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz);
+  lut_low = (vox_svx_low[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz;
+  lut_upp = (vox_svx_upp[2]-sky->htcp_desc.lower[2])/sky->lut_cell_sz;
+  ilut_low = (size_t)lut_low;
+  ilut_upp = (size_t)lut_upp;
+
+  /* A LUT coordinate equal to its search index means that this coordinate lies
+   * strictly on the boundary of a LUT cell. If this coordinate is the upper
+   * limit of a coordinate range, we subtract one from the LUT index to ensure
+   * that the corresponding cell is indeed included in the submitted range
+   * coordinates */
+  if(lut_upp == ilut_upp) {
+    --ilut_upp;
+  }
+
   ASSERT(ilut_low <= ilut_upp);
 
   /* Prepare the iteration over the LES voxels overlapped by the SVX voxel */
   x_h2o_range[0] = DBL_MAX;
   x_h2o_range[1] =-DBL_MAX;
   ivox_z_prev = SIZE_MAX;
-  pos_z = vox_svx_low[2];
-  ASSERT(pos_z >= (double)ilut_low * sky->lut_cell_sz);
-  ASSERT(pos_z <= (double)ilut_upp * sky->lut_cell_sz);
+  ASSERT(vox_svx_low[2] >=
+    sky->htcp_desc.lower[2] + sky->lut_cell_sz*(double)(ilut_low));
+  ASSERT(vox_svx_upp[2] <=
+    sky->htcp_desc.lower[2] + sky->lut_cell_sz*(double)(ilut_upp+1));
 
   /* Iterate over the LUT cells overlapped by the voxel */
   FOR_EACH(ilut, ilut_low, ilut_upp+1) {