rnatm

commit 5b4d7990f21cf83c41e7f289a68b2a6d21b73d21
parent a53a834bd4a3be52e04130405f96bbccf381da9e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 18 Nov 2022 09:48:35 +0100

Fix numerical issue when dealing with k_ext

The way the k_ext tetrahedron is calculated differs from how the voxel
k_ext range is compuuted. This leads to numerical inconsistency where
the k_ext tetrahedron may not be strictly included in the calculated
k_ext range for the voxel. This commits corrects this problem by
slightly expanding the range of k_ext returned.

Diffstat:
M
src/rnatm.c
 | 
16
++++++++++++++--

1 file changed, 14 insertions(+), 2 deletions(-)
diff --git a/src/rnatm.c b/src/rnatm.c
@@ -369,12 +369,24 @@ rnatm_get_k_svx_voxel
     const float ks_min = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)];
     const float ks_max = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)];
 
+    /* Unlike absorption and diffusion coefficients, the range of the
+     * extinction coefficient is not stored by voxel but is recalculated from
+     * ka and ks by voxel; k_ext = ka + ks. However, the extinction coefficient
+     * of a tetrahedron is calculated by first adding ka and ks per component
+     * before adding them to obtain the k_ext of the mixture. Although the
+     * results should be identical, the reorganization of sums produces a
+     * numerical inconsistency. Therefore, the k_ext tetrahedron might not be
+     * strictly included in the range of k_ext computed for the voxel. To deal
+     * with this numerical inaccuracy, we therefore use the following epsilon
+     * to slightly increase the returned range of k_ext */
+    const float epsilon = 1e-6f;
+
     /* Empty voxel => null radiative coefficient */
     if(ka_min > ka_max) { ASSERT(ks_min > ks_max); return 0; }
 
     switch(op) {
-      case RNATM_SVX_OP_MIN: return ka_min + ks_min;
-      case RNATM_SVX_OP_MAX: return ka_max + ks_max;
+      case RNATM_SVX_OP_MIN: return (ka_min + ks_min)*(1.f-epsilon);
+      case RNATM_SVX_OP_MAX: return (ka_max + ks_max)*(1.f+epsilon);
       default: FATAL("Unreachable code\n"); break;
     }
   }