commit 4a6a75527094bdce3c38633999cd0faf304cfa4b
parent 92cfb86bfa8f37a7f8c4c81fd7afe86b9f50db7b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 19 Mar 2026 15:00:10 +0100
Mitigate numerical uncertainty in importance sampling
In order to be used as a sampling probability for a line, the value of a
node at a given wave number must be greater than or equal to the values
of its child nodes. These values are calculated from the node mesh. And
although this criterion is guaranteed by definition for the vertices of
the nodes' polylines, this is no longer true when evaluating the node’s
value elsewhere, via linear interpolation of these vertices. Hence this
commit, which proposes artificially increasing the value of the vertices
of the mesh for each parent node so that this offset compensates for
those introduced by the numerical approximations used to evaluate node
values.
Diffstat:
1 file changed, 20 insertions(+), 4 deletions(-)
diff --git a/src/sln_tree_build.c b/src/sln_tree_build.c
@@ -128,17 +128,18 @@ merge_node_polylines
FOR_EACH(ivtx, vertices_range[0], vertices_range[1]+1) {
const double nu0 = ivtx0 <= ivtx0_max ? vertices[ivtx0].wavenumber : INF;
const double nu1 = ivtx1 <= ivtx1_max ? vertices[ivtx1].wavenumber : INF;
- float nu, ka;
+ double ka;
+ float nu;
if(nu0 < nu1) {
/* The vertex comes from the node0 */
nu = vertices[ivtx0].wavenumber;
- ka = (float)(vertices[ivtx0].ka + eval_ka(tree, node1, nu));
+ ka = vertices[ivtx0].ka + eval_ka(tree, node1, nu);
++ivtx0;
} else if(nu0 > nu1) {
/* The vertex comes from the node1 */
nu = vertices[ivtx1].wavenumber;
- ka = (float)(vertices[ivtx1].ka + eval_ka(tree, node0, nu));
+ ka = vertices[ivtx1].ka + eval_ka(tree, node0, nu);
++ivtx1;
} else {
/* The vertex is shared by node0 and node1 */
@@ -149,7 +150,7 @@ merge_node_polylines
++ivtx1;
}
vertices[ivtx].wavenumber = nu;
- vertices[ivtx].ka = ka;
+ vertices[ivtx].ka = (float)ka;
}
/* Decimate the resulting polyline */
@@ -164,6 +165,21 @@ merge_node_polylines
merged_node->ivertex = vertices_range[0];
merged_node->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1);
+ /* It is necessary to ensure that the recorded vertices define a polyline
+ * along which any value (calculated by linear interpolation) is well above
+ * the sum of the corresponding values of the polylines to be merged. However,
+ * although this is guaranteed by definition for the vertices of the polyline,
+ * numerical uncertainty may nevertheless introduce errors that violate this
+ * criterion. Hence the following adjustment, which slightly increases the ka
+ * of the mesh so as to guarantee this constraint between the mesh of a node
+ * and that of its children */
+ if(tree->args.mesh_type == SLN_MESH_UPPER) {
+ FOR_EACH(ivtx, vertices_range[0], vertices_range[1]+1/*inclusive bound*/) {
+ const double ka = vertices[ivtx].ka;
+ vertices[ivtx].ka = (float)(ka + MMAX(ka*1e-2, 1e-6));
+ }
+ }
+
exit:
return res;
error:
