commit 12357da771597e74ca7b2f9893f536d93eff10f0
parent 8478852f51c71b9fe3cfa1d74e9cafbcedc07b78
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 12 May 2022 14:29:53 +0200
Builds the polyline of all the tree nodes
Diffstat:
8 files changed, 263 insertions(+), 46 deletions(-)
diff --git a/src/sln_line.c b/src/sln_line.c
@@ -283,13 +283,13 @@ error:
/* Append the line mesh into the vertices array */
static res_T
save_line_mesh
- (const struct sln_tree* tree,
+ (struct sln_tree* tree,
const size_t iline,
const struct molecule_params* mol_params,
const double spectral_range[2],
const struct darray_double* wavenumbers,
const struct darray_double* values,
- struct darray_vertex* vertices)
+ size_t vertices_range[2]) /* Range into which the line vertices are saved */
{
const struct line* line = NULL;
size_t nvertices;
@@ -298,12 +298,12 @@ save_line_mesh
size_t ivertex;
size_t i;
res_T res = RES_OK;
- ASSERT(tree && mol_params && wavenumbers && values && vertices);
+ ASSERT(tree && mol_params && wavenumbers && values && vertices_range);
ASSERT(darray_double_size_get(wavenumbers) == darray_double_size_get(values));
ASSERT(spectral_range && spectral_range[0] <= spectral_range[1]);
line = darray_line_cdata_get(&tree->lines) + iline;
- nvertices = darray_vertex_size_get(vertices);
+ nvertices = darray_vertex_size_get(&tree->vertices);
nwavenumbers = darray_double_size_get(wavenumbers);
/* Compute the overall number of vertices of the line */
@@ -312,7 +312,7 @@ save_line_mesh
- 1;/* Do not duplicate the line center */
/* Allocate the list of line vertices */
- res = darray_vertex_resize(vertices, nvertices + line_nvertices);
+ res = darray_vertex_resize(&tree->vertices, nvertices + line_nvertices);
if(res != RES_OK) {
log_err(tree->sln, "Error allocating the line vertices -- %s.\n",
res_to_cstr(res));
@@ -324,7 +324,7 @@ save_line_mesh
/* Emit the 1st vertex for the lower bound of the spectral range if the line
* is clipped by it */
if(line->wavenumber - mol_params->cutoff <= spectral_range[0]) {
- struct sln_vertex* vtx = darray_vertex_data_get(vertices) + i;
+ struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i;
vtx->wavenumber = (float)spectral_range[0];
vtx->ka = (float)line_value(tree, iline, vtx->wavenumber);
++i;
@@ -332,7 +332,7 @@ save_line_mesh
/* Copy the vertices of the line for its lower half */
FOR_EACH_REVERSE(ivertex, nwavenumbers-1, 0) {
- struct sln_vertex* vtx = darray_vertex_data_get(vertices) + i;
+ struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i;
const double ka = darray_double_cdata_get(values)[ivertex];
const double nu = /* Mirror the wavenumber */
2*line->wavenumber - darray_double_cdata_get(wavenumbers)[ivertex];
@@ -349,7 +349,7 @@ save_line_mesh
/* Copy the vertices of the line for its upper half */
FOR_EACH(ivertex, 0, nwavenumbers) {
- struct sln_vertex* vtx = darray_vertex_data_get(vertices) + i;
+ struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i;
const double ka = darray_double_cdata_get(values)[ivertex];
const double nu = darray_double_cdata_get(wavenumbers)[ivertex];
@@ -367,7 +367,7 @@ save_line_mesh
/* Emit the last vertex for the upper bound of the spectral range if the line
* is clipped by it */
if(line->wavenumber + mol_params->cutoff >= spectral_range[1]) {
- struct sln_vertex* vtx = darray_vertex_data_get(vertices) + i;
+ struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i;
vtx->wavenumber = (float)spectral_range[1];
vtx->ka = (float)line_value(tree, iline, vtx->wavenumber);
++i;
@@ -376,12 +376,16 @@ save_line_mesh
/* Several vertices could be skipped due to the clipping of the line by the
* spectral range. Resize the list of vertices to the effective number of
* registered vertices. */
- CHK(darray_vertex_resize(vertices, i) == RES_OK);
+ CHK(darray_vertex_resize(&tree->vertices, i) == RES_OK);
+
+ /* Setup the range of the line vertices */
+ vertices_range[0] = nvertices;
+ vertices_range[1] = i-1;
exit:
return res;
error:
- darray_vertex_resize(vertices, nvertices);
+ darray_vertex_resize(&tree->vertices, nvertices);
goto exit;
}
@@ -435,21 +439,32 @@ line_value
const double wavenumber)
{
const struct line* line = NULL;
- double profile;
+ const struct shtr_line* shtr_line = NULL;
+ const struct molecule_params* mol_params = NULL;
ASSERT(tree && iline < darray_line_size_get(&tree->lines));
+ /* Retrieve the molecular parameters of the line to be mesh */
+ SHTR(lines_view_get_line(tree->lines_view, iline, &shtr_line));
+ mol_params = tree_get_molecule_params(tree, shtr_line->molecule_id);
+
line = darray_line_cdata_get(&tree->lines) + iline;
- profile = sln_compute_voigt_profile
- (wavenumber, line->wavenumber, line->gamma_d, line->gamma_l);
- return line->profile_factor * profile;
+
+ if(wavenumber < line->wavenumber - mol_params->cutoff
+ || wavenumber > line->wavenumber + mol_params->cutoff) {
+ return 0;
+ } else {
+ const double profile = sln_compute_voigt_profile
+ (wavenumber, line->wavenumber, line->gamma_d, line->gamma_l);
+ return line->profile_factor * profile;
+ }
}
res_T
line_mesh
- (const struct sln_tree* tree,
+ (struct sln_tree* tree,
const size_t iline,
const size_t nvertices_hint,
- struct darray_vertex* vertices)
+ size_t vertices_range[2])
{
struct darray_double values; /* List of evaluated values */
struct darray_double wavenumbers; /* List of considered wavenumbers */
@@ -458,7 +473,7 @@ line_mesh
const struct line* line = NULL;
size_t nvertices_adjusted = 0; /* computed from nvertices_hint */
res_T res = RES_OK;
- ASSERT(tree && nvertices_hint && vertices);
+ ASSERT(tree && nvertices_hint);
ASSERT(iline < darray_line_size_get(&tree->lines));
line = darray_line_cdata_get(&tree->lines) + iline;
@@ -487,7 +502,7 @@ line_mesh
eval_mesh(tree, iline, &wavenumbers, &values);
res = save_line_mesh(tree, iline, mol_params, tree->wavenumber_range,
- &wavenumbers, &values, vertices);
+ &wavenumbers, &values, vertices_range);
if(res != RES_OK) goto error;
exit:
diff --git a/src/sln_line.h b/src/sln_line.h
@@ -54,12 +54,12 @@ line_value
const size_t iline,
const double wavenumber);
-/* Append the line vertices into `vertices' */
extern LOCAL_SYM res_T
line_mesh
- (const struct sln_tree* tree,
+ (struct sln_tree* tree,
const size_t iline,
const size_t nvertices_hint,
- struct darray_vertex* vertices);
+ /* Range of generated line vertices. Upper bound is exclusive */
+ size_t vertices_range[2]);
#endif /* SLN_LINE_H */
diff --git a/src/sln_polyline.c b/src/sln_polyline.c
@@ -92,6 +92,22 @@ find_falsest_vertex
*out_err = err_max;
}
+static INLINE void
+check_polyline_vertices
+ (const struct sln_vertex* vertices,
+ const size_t vertices_range[2])
+{
+#ifdef NDEBUG
+ (void)vertices, (void)vertices_range;
+#else
+ size_t i;
+ ASSERT(vertices);
+ FOR_EACH(i, vertices_range[0]+1, vertices_range[1]+1) {
+ CHK(vertices[i].wavenumber > vertices[i-1].wavenumber);
+ }
+#endif
+}
+
/*******************************************************************************
* Local function
******************************************************************************/
@@ -117,6 +133,7 @@ polyline_decimate
res_T res = RES_OK;
ASSERT(vertices && vertices_range && err >= 0);
ASSERT(vertices_range[0] < vertices_range[1]);
+ check_polyline_vertices(vertices, vertices_range);
darray_size_t_init(sln->allocator, &stack);
@@ -168,6 +185,8 @@ polyline_decimate
vertices_range[1] = ivtx - 1;
+ check_polyline_vertices(vertices, vertices_range);
+
exit:
darray_size_t_release(&stack);
return res;
diff --git a/src/sln_polyline.h b/src/sln_polyline.h
@@ -26,7 +26,8 @@ struct sln_device;
struct sln_vertex;
/* In place simplification of a polyline. Given a curve composed of line
- * segments, compute a similar curve with fewer points. */
+ * segments, compute a similar curve with fewer points. The range of this new
+ * polyline is returned in `vertices_range'. */
extern LOCAL_SYM res_T
polyline_decimate
(struct sln_device* sln,
diff --git a/src/sln_tree.c b/src/sln_tree.c
@@ -485,7 +485,7 @@ size_t
sln_node_get_lines_count(const struct sln_node* node)
{
ASSERT(node);
- return node->range[1] - node->range[0] + 1;/*Both boundaries are inclusives*/
+ return node->range[1] - node->range[0] + 1/*Both boundaries are inclusives*/;
}
res_T
@@ -509,7 +509,7 @@ sln_node_get_mesh
struct sln_mesh* mesh)
{
if(!tree || !node || !mesh) return RES_BAD_ARG;
- mesh->vertices = darray_vertex_cdata_get(&tree->vertices) + node->ivertices;
+ mesh->vertices = darray_vertex_cdata_get(&tree->vertices) + node->ivertex;
mesh->nvertices = node->nvertices;
return RES_OK;
}
diff --git a/src/sln_tree_build.c b/src/sln_tree_build.c
@@ -25,6 +25,9 @@
#include <rsys/cstr.h>
+/* STACK_SIZE is set to the maximum depth of the partitioning tree generated by
+ * the tree_build function. This is actually the maximum stack size where tree
+ * nodes are pushed by the recursive build process */
#define STACK_SIZE 64
/*******************************************************************************
@@ -39,7 +42,7 @@ ui64_to_ui32(const uint64_t ui64)
}
static INLINE res_T
-setup_leaf
+build_leaf_polyline
(struct sln_tree* tree,
const struct sln_tree_create_args* args,
struct sln_node* leaf)
@@ -51,12 +54,10 @@ setup_leaf
ASSERT(leaf->range[0] == leaf->range[1]);
/* Line meshing */
- vertices_range[0] = darray_vertex_size_get(&tree->vertices);
- res = line_mesh(tree, leaf->range[0], args->nvertices_hint, &tree->vertices);
+ res = line_mesh(tree, leaf->range[0], args->nvertices_hint, vertices_range);
if(res != RES_OK) goto error;
- vertices_range[1] = darray_vertex_size_get(&tree->vertices) - 1;
- /* Decimating the line mesh */
+ /* Decimate the line mesh */
res = polyline_decimate(tree->sln, darray_vertex_data_get(&tree->vertices),
vertices_range, args->mesh_decimation_err);
if(res != RES_OK) {
@@ -64,11 +65,12 @@ setup_leaf
res_to_cstr(res));
goto error;
}
+
+ /* Shrink the size of the vertices */
darray_vertex_resize(&tree->vertices, vertices_range[1] + 1);
- /* Setup leaf data */
- leaf->offset = 0;
- leaf->ivertices = vertices_range[0];
+ /* Setup the leaf polyline */
+ leaf->ivertex = vertices_range[0];
leaf->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1);
exit:
@@ -77,11 +79,167 @@ error:
goto exit;
}
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-tree_build
+static INLINE double
+compute_ka
+ (const struct sln_tree* tree,
+ const double wavenumber,
+ const size_t line_range[2])
+{
+ double ka = 0;
+ size_t iline = 0;
+ ASSERT(tree && wavenumber >= 0 && line_range);
+ ASSERT(line_range[0] <= line_range[1]);
+
+ FOR_EACH(iline, line_range[0], line_range[1]+1) {
+ ka += line_value(tree, iline, wavenumber);
+ }
+ return ka;
+}
+
+static res_T
+merge_node_polylines
+ (struct sln_tree* tree,
+ const struct sln_tree_create_args* args,
+ const struct sln_node* node0,
+ const struct sln_node* node1,
+ struct sln_node* merged_node)
+{
+ struct sln_vertex* vertices = NULL;
+ size_t vertices_range[2];
+ size_t ivtx;
+ size_t ivtx0, ivtx0_max;
+ size_t ivtx1, ivtx1_max;
+ res_T res = RES_OK;
+ ASSERT(tree && node0 && node1 && merged_node);
+
+ /* Define the vertices range of the merged polyline */
+ vertices_range[0] = darray_vertex_size_get(&tree->vertices);
+ vertices_range[1] = vertices_range[0] + node0->nvertices + node1->nvertices - 1;
+
+ /* Allocate the memory space to store the new polyline */
+ res = darray_vertex_resize(&tree->vertices, vertices_range[1] + 1);
+ if(res != RES_OK) {
+ log_err(tree->sln, "Error in merging polylines -- %s.\n", res_to_cstr(res));
+ goto error;
+ }
+ vertices = darray_vertex_data_get(&tree->vertices);
+
+ ivtx0 = node0->ivertex;
+ ivtx1 = node1->ivertex;
+ ivtx0_max = ivtx0 + node0->nvertices - 1;
+ ivtx1_max = ivtx1 + node1->nvertices - 1;
+ 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;
+
+ if(nu0 < nu1) {
+ /* The vertex comes from the node0 */
+ nu = vertices[ivtx0].wavenumber;
+ ka = (float)(vertices[ivtx0].ka + compute_ka(tree, nu, node1->range));
+ ++ivtx0;
+ } else if(nu0 > nu1) {
+ /* The vertex comes from the node1 */
+ nu = vertices[ivtx1].wavenumber;
+ ka = (float)(vertices[ivtx1].ka + compute_ka(tree, nu, node0->range));
+ ++ivtx1;
+ } else {
+ /* The vertex is shared by node0 and node1 */
+ nu = vertices[ivtx0].wavenumber;
+ ka = vertices[ivtx0].ka + vertices[ivtx1].ka;
+ --vertices_range[1]; /* Remove duplicate */
+ ++ivtx0;
+ ++ivtx1;
+ }
+ vertices[ivtx].wavenumber = nu;
+ vertices[ivtx].ka = ka;
+ }
+
+ /* Decimate the resulting polyline */
+ res = polyline_decimate(tree->sln, darray_vertex_data_get(&tree->vertices),
+ vertices_range, args->mesh_decimation_err);
+ if(res != RES_OK) {
+ log_err(tree->sln, "Error while decimating the line mesh -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ /* Shrink the size of the vertices */
+ darray_vertex_resize(&tree->vertices, vertices_range[1] + 1);
+
+ /* Setup the node polyline */
+ merged_node->ivertex = vertices_range[0];
+ merged_node->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+build_polylines
+ (struct sln_tree* tree,
+ const struct sln_tree_create_args* args)
+{
+ size_t stack[STACK_SIZE*2];
+ size_t istack = 0;
+ size_t inode = 0;
+ res_T res = RES_OK;
+ ASSERT(tree && args);
+
+ #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id))
+ #define IS_LEAF(Id) (NODE(Id)->offset == 0)
+
+ /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */
+ stack[istack++] = SIZE_MAX;
+
+ inode = 0; /* Root node */
+ while(inode != SIZE_MAX) {
+
+ if(IS_LEAF(inode)) {
+ res = build_leaf_polyline(tree, args, NODE(inode));
+ if(res != RES_OK) goto error;
+
+ inode = stack[--istack]; /* Pop the next node */
+
+ } else {
+ const size_t ichild0 = inode + NODE(inode)->offset + 0;
+ const size_t ichild1 = inode + NODE(inode)->offset + 1;
+
+ /* Child nodes have their polyline created */
+ if(NODE(ichild0)->nvertices) {
+ ASSERT(NODE(ichild1)->nvertices != 0);
+
+ /* Build the polyline of the current node by merging the polylines of
+ * its children */
+ res = merge_node_polylines
+ (tree, args, NODE(ichild0), NODE(ichild1), NODE(inode));
+ if(res != RES_OK) goto error;
+ inode = stack[--istack];
+
+ } else {
+ ASSERT(NODE(ichild1)->nvertices == 0);
+
+ stack[istack++] = inode; /* Push the current node */
+ stack[istack++] = ichild1; /* Push child1 */
+
+ /* Recursively build the polyline of child0 */
+ inode = ichild0;
+ }
+ }
+ }
+
+ #undef NODE
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+partition_lines
(struct sln_tree* tree,
const struct sln_tree_create_args* args)
{
@@ -105,18 +263,19 @@ tree_build
/* Setup the root node */
NODE(0)->range[0] = 0;
NODE(0)->range[1] = nlines - 1;
- inode = 0;
+ /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */
+ stack[istack++] = SIZE_MAX;
+
+ inode = 0; /* Root node */
while(inode != SIZE_MAX) {
/* #lines into the node */
nlines = NODE(inode)->range[1] - NODE(inode)->range[0] + 1;
/* Make a leaf */
if(nlines <= args->max_nlines_per_leaf) {
- res = setup_leaf(tree, args, NODE(inode));
- if(res != RES_OK) goto error;
-
- inode = istack ? stack[--istack] : SIZE_MAX; /* Pop the next node */
+ NODE(inode)->offset = 0;
+ inode = stack[--istack]; /* Pop the next node */
/* Split the node */
} else {
@@ -159,3 +318,24 @@ error:
darray_node_clear(&tree->nodes);
goto exit;
}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+tree_build
+ (struct sln_tree* tree,
+ const struct sln_tree_create_args* args)
+{
+ res_T res = RES_OK;
+
+ res = partition_lines(tree, args);
+ if(res != RES_OK) goto error;
+ res = build_polylines(tree, args);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/sln_tree_c.h b/src/sln_tree_c.h
@@ -30,10 +30,9 @@ struct sln_tree_create_args;
struct shtr_lines_view;
struct sln_node { /* 32 Bytes */
- /* Range of the line indices corresponding to the node. Both the lower and
- * upper indices are included into the node */
+ /* Range of the line indices corresponding to the node */
uint64_t range[2];
- uint64_t ivertices; /* Index toward the 1st vertex */
+ uint64_t ivertex; /* Index toward the 1st vertex */
uint32_t nvertices; /* #vertices */
uint32_t offset; /* Offset toward the node's children (left then right) */
};
diff --git a/src/test_sln_tree.c b/src/test_sln_tree.c
@@ -386,6 +386,9 @@ test_tree
CHK(sln_node_get_mesh(tree, node, NULL) == RES_BAD_ARG);
CHK(sln_node_get_mesh(tree, node, &mesh) == RES_OK);
+ CHK(node = sln_tree_get_root(tree));
+ CHK(sln_node_get_mesh(tree, node, &mesh) == RES_OK);
+
dump_line(stdout, &mesh);
check_tree(tree, lines_list);
