star-line

sln_tree_build.c (53029B)

---

 /* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com)
  * Copyright (C) 2026 Université de Lorraine
  * Copyright (C) 2022 Centre National de la Recherche Scientifique
  * Copyright (C) 2022 Université Paul Sabatier
  *
  * This file is part of Star-Line.
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
 
 #include "sln.h"
 #include "sln_device_c.h"
 #include "sln_polyline.h"
 #include "sln_tree_c.h"
 
 #include <star/shtr.h>
 
 #include <rsys/cstr.h>
 #include <rsys/math.h>
 
 #include <omp.h>
 
 /* Structure used to track the progress of a construction phase */
 struct progress {
   /* Absolute progression */
   size_t total; /* Total number of items to process */
   size_t processed; /* Number of items already processed */
 
   /* Relative progression */
   int pcent; /* Currently displayed percentage */
   int pcent_step; /* Increment between displayed percentages */
 
   const char* msg; /* Message to add before the displayed percentage */
 };
 #define PROGRESS_DEFAULT__ {0,0,0,10,NULL}
 static const struct progress PROGRESS_DEFAULT = PROGRESS_DEFAULT__;
 
 /* Structure containing temporary variables used to create polyline on a leaf
  * These variables are not declared locally within the function responsible for
  * constructing this polyline in order to avoid having to allocate and
  * deallocate them with each call to the function. Instead, they are passed as
  * function inputs so as to share, as much as possible, the memory space
  * allocated during previous calls to the function, thereby reducing the cost of
  * allocation and deallocation. */
 struct build_leaf_polyline_scratch {
   struct darray_vertex vertices; /* Polyline vertices */
   struct darray_node nodes; /* Dummy leaf nodes */
 
   /* Temp vertex buffer used when merging polylines into a single one */
   struct darray_vertex vertices_tmp;
 };
 
 static res_T
 merge_child_polylines
   (struct sln_tree* tree,
    const size_t inode, /* Node at which child polylines are merged */
    const unsigned nchildren, /* #children to merge */
    struct darray_node* nodes, /* List of nodes */
    struct darray_vertex* vertices); /* List of polyline vertices */
 
 static res_T
 collapse_child_polylines
   (struct sln_tree* tree,
    const size_t inode, /* Node at which child polylines are merged */
    const unsigned nchildren, /* #children to collapse */
    struct darray_node* nodes, /* List of nodes */
    struct darray_vertex* scratch, /* Temporary buffer of vertices */
    struct darray_vertex* vertices); /* List of polyline vertices */
 
 /* A multiplier to be applied to the calculated ka values in order to mitigate
  * numerical issues related to ka interpolation in the case of a tree with an
  * upper bound. This interpolation must ensure that the interpolated ka value is
  * well above the actual ka value, which may not be the case due to numerical
  * inaccuracies. Hence this constant, which defines a percentage of ka to be
  * added to the calculated ka values to ensure that any interpolated value is
  * well above the actual ka value. */
 #define KA_ADJUSTEMENT 1.01
 
 /* Maximum queue size used for the breadth-first tree traversal. For a perfectly
  * balanced tree, this corresponds to the maximum number of leaves a tree can
  * have. Note that this value does not need to be too high, since the
  * breadth-first traversal is only used to partition the first few levels of the
  * tree until a sufficient number of subtrees have been identified so that they
  * can then be constructed in parallel using a depth-first algorithm */
 #define QUEUE_SIZE_MAX 256
 
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
 static FINLINE uint32_t
 ui64_to_ui32(const uint64_t ui64)
 {
   if(ui64 > UINT32_MAX)
     VFATAL("%s: overflow %lu.\n", ARG2(FUNC_NAME, ((unsigned long)ui64)));
   return (uint32_t)ui64;
 }
 
 static FINLINE unsigned
 size_t_to_unsigned(const size_t sz)
 {
   if(sz > UINT_MAX)
     VFATAL("%s: overflow %lu.\n", ARG2(FUNC_NAME, ((unsigned long)sz)));
   return (unsigned)sz;
 }
 
 static void
 progress_print(struct sln_device* dev, const struct progress* progress)
 {
   ASSERT(dev && progress);
   if(progress->msg) {
     INFO(dev, "%s%3d%%\n", progress->msg, progress->pcent);
   } else {
     INFO(dev, "%3d%%\n", progress->pcent);
   }
 }
 
 static void
 progress_update
   (struct sln_device* dev,
    struct progress* progress,
    const size_t count)
 {
   int pcent = 0;
   ASSERT(dev && progress);
 
   #pragma omp critical
   {
     progress->processed += count;
     ASSERT(progress->processed <= progress->total);
 
     pcent = (int)
       ( (double)progress->processed*100.0
       / (double)progress->total
       + 0.5 /* round */);
 
     if(pcent/progress->pcent_step > progress->pcent/progress->pcent_step) {
       progress->pcent = pcent;
       progress_print(dev, progress);
     }
   }
 }
 
 static INLINE void
 build_leaf_polyline_scratch_init
   (struct mem_allocator* allocator,
    struct build_leaf_polyline_scratch* scratch)
 {
   ASSERT(scratch);
   darray_vertex_init(allocator, &scratch->vertices);
   darray_vertex_init(allocator, &scratch->vertices_tmp);
   darray_node_init(allocator, &scratch->nodes);
 }
 
 static INLINE void
 build_leaf_polyline_scratch_release
   (struct build_leaf_polyline_scratch* scratch)
 {
   ASSERT(scratch);
   darray_vertex_release(&scratch->vertices);
   darray_vertex_release(&scratch->vertices_tmp);
   darray_node_release(&scratch->nodes);
 }
 
 static INLINE res_T
 build_leaf_polyline_from_1line
   (struct sln_tree* tree,
    struct sln_node* leaf,
    struct darray_vertex* vertices)
 {
   size_t vertices_range[2];
   res_T res = RES_OK;
 
   ASSERT(tree && leaf && vertices);
 
   /* Assume that there is only one line per leaf */
   ASSERT(leaf->range[0] == leaf->range[1]);
 
   /* Line meshing */
   res = line_mesh
     (/* in */  tree, leaf->range[0], tree->args.nvertices_hint,
      /* out */ vertices, vertices_range);
   if(res != RES_OK) goto error;
 
   /* Decimate the line mesh */
   res = polyline_decimate(tree->sln, darray_vertex_data_get(vertices),
      vertices_range, (float)tree->args.mesh_decimation_err, tree->args.mesh_type);
   if(res != RES_OK) goto error;
 
   /* Shrink the size of the vertices */
   darray_vertex_resize(vertices, vertices_range[1] + 1);
 
   /* Setup the leaf polyline  */
   leaf->ivertex = vertices_range[0];
   leaf->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1);
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 /* Build the polyline of a leaf node containing more than one line.
  *
  * Each line is drawn as if it were the only one on the leaf. This allows the
  * line meshing routine to be used, and consequently the algorithm that relies
  * on the properties of the lines (symmetry, monotonicity on one half). Next,
  * this set of polylines is merged/collapsed into a single polyline, which is
  * the leaf polyline.
  *
  * To reuse the merging/collapsing designed for internal nodes, the function
  * treats a leaf as if it were an internal node. A list of temporary nodes is
  * thus created, containing not only the leaf node but also its “children”
  * i.e., its lines.
  *
  * Once created, the leaf’s polyline is finally copied into the tree’s vertex
  * buffer */
 static INLINE res_T
 build_leaf_polyline_from_Nlines
   (struct sln_tree* tree,
    struct sln_node* leaf,
    struct build_leaf_polyline_scratch* scratch,
    struct darray_vertex* vertices)
 {
   const struct sln_vertex* src = NULL;
   struct sln_vertex* dst = NULL;
   size_t memsz = 0;
 
   size_t nnodes = 0;
   unsigned nlines = 0; /* Number of lines in the leaf */
   unsigned i = 0;
   res_T res = RES_OK;
 
   ASSERT(tree && leaf && scratch && vertices);
 
   /* Compute the number of lines of the leaf */
   nlines = size_t_to_unsigned(leaf->range[1] - leaf->range[0] + 1/*inclusive*/);
   ASSERT(nlines > 1); /* Assume there is more than one line per leaf */
 
   nnodes = nlines + 1/*leaf*/;
 
   /* Helper macro */
   #define NODE(Id) (darray_node_data_get(&scratch->nodes) + (Id))
 
   /* Allocate the temporary list of nodes, one node per leaf to which on
    * additionnal node is added for the leaf */
   res = darray_node_resize(&scratch->nodes, nnodes);
   if(res != RES_OK) goto error;
   memset(NODE(0), 0, sizeof(struct sln_node)*nnodes);
 
   /* Clean up the temporary list of vertices */
   darray_vertex_clear(&scratch->vertices);
 
   /* The leaf node will be the first one. Its "child" representing the node of
    * each lines, will be store after it in the node list */
   NODE(0)->offset = 1;
   NODE(0)->range[0] = leaf->range[0];
   NODE(0)->range[1] = leaf->range[1];
 
   FOR_EACH(i, 0, nlines) {
     size_t vertices_range[2] = {0, 0}; /* Range in the vertex buffer */
     const size_t ichild = NODE(0)->offset + i; /* Index of the child */
     const size_t iline = leaf->range[0] + i;
 
     /* Mesh the line in the temporary vertex buffer */
     res = line_mesh
       (/* in */  tree, iline, tree->args.nvertices_hint,
        /* out */ &scratch->vertices, vertices_range);
     if(res != RES_OK) goto error;
 
     /* Decimate the line mesh */
     res = polyline_decimate(tree->sln,
       darray_vertex_data_get(&scratch->vertices), vertices_range,
       (float)tree->args.mesh_decimation_err, tree->args.mesh_type);
     if(res != RES_OK) goto error;
 
     NODE(ichild)->ivertex = vertices_range[0];
     NODE(ichild)->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1);
     NODE(ichild)->range[0] = iline;
     NODE(ichild)->range[1] = iline;
   }
 
   /* Merge/collapse the polylines of each lines associated to nodes.
    * These nodes are the "children" of the leaf */
   if(tree->args.collapse_polylines) {
     res = collapse_child_polylines(tree, 0, nlines,
       &scratch->nodes, &scratch->vertices_tmp, &scratch->vertices);
   } else {
     res = merge_child_polylines(tree, 0, nlines,
       &scratch->nodes, &scratch->vertices);
   }
   if(res != RES_OK) goto error;
 
   /* Setup the leaf */
   leaf->ivertex = darray_vertex_size_get(vertices);
   leaf->nvertices = NODE(0/*leaf*/)->nvertices;
 
   /* Copy the leaf vertices from temporary buffer to the vertex buffer of the
    * tree */
   res = darray_vertex_resize(vertices, leaf->ivertex + leaf->nvertices);
   if(res != RES_OK) goto error;
   src = darray_vertex_cdata_get(&scratch->vertices) + NODE(0/*leaf*/)->ivertex;
   dst = darray_vertex_data_get(vertices) + leaf->ivertex;
   memsz = leaf->nvertices * sizeof(*src);
   memcpy(dst, src, memsz);
 
   #undef NODE
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 static INLINE res_T
 build_leaf_polyline
   (struct sln_tree* tree,
    struct sln_node* leaf,
    struct build_leaf_polyline_scratch* scratch,
    struct darray_vertex* vertices)
 {
   size_t nlines = 0; /* Number of lines in the leaf */
   ASSERT(tree && leaf);
 
   nlines = leaf->range[1] - leaf->range[0] + 1;
 
   if(nlines == 1) {
     return build_leaf_polyline_from_1line(tree, leaf, vertices);
   } else {
     return build_leaf_polyline_from_Nlines(tree, leaf, scratch, vertices);
   }
 }
 
 static INLINE double
 eval_ka
   (const struct sln_node* node,
    const struct sln_vertex* vertices,
    const double wavenumber)
 {
   struct sln_mesh mesh = SLN_MESH_NULL;
   double ka = 0;
   ASSERT(node && vertices);
 
   /* Whether the mesh to be constructed corresponds to the spectrum or its upper
    * limit, use the node mesh to calculate the value of ka at a given wave
    * number. Calculating the value from the node lines would take far too long*/
   mesh.vertices = vertices + node->ivertex;
   mesh.nvertices = node->nvertices;
   ka = sln_mesh_eval(&mesh, wavenumber);
   return ka;
 }
 
 /* Merge all child polylines in a single step. The polylines are merged before
  * being decimated */
 res_T
 merge_child_polylines
   (struct sln_tree* tree,
    const size_t inode,
    const unsigned nchildren,
    struct darray_node* nodes,
    struct darray_vertex* vertex_list)
 {
   /* Helper constant */
   static const size_t NO_MORE_VERTEX = SIZE_MAX;
 
   /* Polyline vertices */
   #define NCHILDREN_MAX \
     ( SLN_TREE_ARITY_MAX > SLN_LEAF_NLINES_MAX \
     ? SLN_TREE_ARITY_MAX : SLN_LEAF_NLINES_MAX)
   size_t children_ivtx[NCHILDREN_MAX] = {0};
 
   size_t vertices_range[2] = {0,0};
   struct sln_vertex* vertices = NULL;
   size_t ivtx = 0;
   size_t nvertices = 0;
 
   /* Miscellaneous */
   unsigned i = 0;
   res_T res = RES_OK;
 
   /* Pre-conditions */
   ASSERT(tree && inode < darray_node_size_get(nodes));
   ASSERT(nchildren >= 2 && nchildren <= NCHILDREN_MAX);
 
   #define NODE(Id) (darray_node_data_get(nodes) + (Id))
 
   /* Compute the number of vertices to be merged,
    * i.e., the sum of vertices of the children */
   nvertices = 0;
   FOR_EACH(i, 0, nchildren) {
     const size_t ichild = inode + NODE(inode)->offset + i;
     nvertices += NODE(ichild)->nvertices;
   }
 
   /* Define the vertices range of the merged polyline */
   vertices_range[0] = darray_vertex_size_get(vertex_list);
   vertices_range[1] = vertices_range[0] + nvertices - 1/*inclusive bound*/;
 
   /* Allocate the memory space to store the new polyline */
   res = darray_vertex_resize(vertex_list, vertices_range[1]+1);
   if(res != RES_OK) {
     ERROR(tree->sln, "Error in merging polylines -- %s.\n", res_to_cstr(res));
     goto error;
   }
   vertices = darray_vertex_data_get(vertex_list);
 
   /* Initialize the vertex index list. For each child, the initial value
    * corresponds to the index of its first vertex. This index will be
    * incremented as vertices are merged into the parent polyline. */
   FOR_EACH(i, 0, nchildren) {
     const size_t ichild = inode + NODE(inode)->offset + i;
     children_ivtx[i] = NODE(ichild)->ivertex;
   }
 
   FOR_EACH(ivtx, vertices_range[0], vertices_range[1]+1/*inclusive bound*/) {
     double ka = 0;
     double nu = INF;
 
     /* The number of vertices corresponding to the current wave number for which
      * the parent ka is calculated. It is at least equal to one, since this nu
      * is defined by the child vertices, but may be greater if multiple children
      * share the same vertex, i.e., a ka value calculated for the same nu */
     unsigned nvertices_merged = 0;
 
     /* Find the minimum wave number among the vertices of the child vertices
      * that are candidates for merging */
     FOR_EACH(i, 0, nchildren) {
       const size_t child_ivtx = children_ivtx[i];
       if(child_ivtx != NO_MORE_VERTEX) {
         nu = MMIN(nu, vertices[child_ivtx].wavenumber);
       }
     }
     ASSERT(nu != INF); /* At least one vertex must have been found */
 
     /* Compute the value of ka at the wave number determined above */
     FOR_EACH(i, 0, nchildren) {
       const size_t child_ivtx = children_ivtx[i];
       const struct sln_node* child = NODE(inode + NODE(inode)->offset + i);
 
       if(child_ivtx == NO_MORE_VERTEX
       || nu != vertices[child_ivtx].wavenumber) {
         /* The wave number does not correspond to a vertex in the current
          * child's mesh. Therefore, its contribution to the parent node's ka is
          * computed  */
         ka += eval_ka(child, vertices, nu);
 
       } else {
         /* The wave number is the one for which the child node stores a ka
          * value. Add it to the parent node's ka value and designate the child's
          * next vertex as a candidate for merging into the parent. The exception
          * is when all vertices of the child have already been merged. In this
          * case, report that the child no longer has any candidate vertices */
         ka += vertices[child_ivtx].ka;
         ++children_ivtx[i];
         if(children_ivtx[i] >= child->ivertex + child->nvertices) {
           children_ivtx[i] = NO_MORE_VERTEX;
         }
         ++nvertices_merged; /* Record that a vertex has been merged */
       }
     }
 
     /* Setup the parent vertex */
     vertices[ivtx].wavenumber = (float)nu;
     vertices[ivtx].ka = (float)ka;
 
     /* If multiple child vertices have been merged, then a single wave number
      * corresponds to a vertex with multiple children. The number of parent
      * vertices is therefore no longer the sum of the number of its children's
      * vertices, since some vertices are duplicated. Hence the following
      * adjustment, which removes the duplicate vertices. */
     vertices_range[1] -= (nvertices_merged-1);
   }
 
   /* Decimate the resulting polyline */
   res = polyline_decimate(tree->sln, darray_vertex_data_get(vertex_list),
      vertices_range, (float)tree->args.mesh_decimation_err, tree->args.mesh_type);
   if(res != RES_OK) goto error;
 
   /* Setup the node polyline */
   NODE(inode)->ivertex = vertices_range[0];
   NODE(inode)->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*KA_ADJUSTEMENT);
     }
   }
 
   /* Shrink the size of the vertices */
   darray_vertex_resize(vertex_list, vertices_range[1] + 1);
 
   #undef NODE
   #undef NCHILDREN_MAX
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 /* Merge child polylines by combining them in pairs (the resulting polyline is
  * then simplified), and repeat this process until only a single polyline
  * remains. This polyline becomes the node's polyline */
 static res_T
 collapse_child_polylines
   (struct sln_tree* tree,
    const size_t inode,
    const unsigned nchildren,
    struct darray_node* nodes,
    struct darray_vertex* scratch,
    struct darray_vertex* vertex_list)
 {
   /* Polylines to be collapsed, i.e., the ids of their first and last vertices */
   #define NCHILDREN_MAX \
     ( SLN_TREE_ARITY_MAX > SLN_LEAF_NLINES_MAX \
     ? SLN_TREE_ARITY_MAX : SLN_LEAF_NLINES_MAX)
   size_t poly_parts[NCHILDREN_MAX][2] = {0};
   size_t nparts;
 
   /* Indices of the first and last vertices of the resulting polyline.
    * These indices are absolute to the tree's vertex list */
   size_t vertices_range[2] = {0,0};
 
   /* Redux double buffering */
   struct sln_vertex* buf[2] = {NULL, NULL};
   int r, w; /* Index of the buffer in read/write */
 
   /* Miscellaneous */
   struct sln_vertex* vertices = NULL; /* Pointer to the tree's vertex buffer */
   size_t range_merge[2] = {0,0}; /* vertex range of a merged polyline */
   size_t nvertices = 0; /* #vertices of the resulting polyline */
   size_t ivtx = 0;
   unsigned i = 0;
   int ncollapses = 0; /* Number of collapse steps */
   res_T res = RES_OK;
 
   /* Pre-conditions */
   ASSERT(tree && inode < darray_node_size_get(nodes));
   ASSERT(nchildren >= 2 && nchildren <= NCHILDREN_MAX);
 
   #define NODE(Id) (darray_node_data_get(nodes) + (Id))
 
   /* Compute the number of vertices to be merged,
    * i.e., the sum of vertices of the children */
   nvertices = 0;
   FOR_EACH(i, 0, nchildren) {
     const size_t ichild = inode + NODE(inode)->offset + i;
     nvertices += NODE(ichild)->nvertices;
   }
 
   /* Define the vertices range of the merged polyline */
   vertices_range[0] = darray_vertex_size_get(vertex_list);
   vertices_range[1] = vertices_range[0] + nvertices - 1/*inclusive bound*/;
 
   /* Allocate the memory space required to store the new polyline and the
    * temporary buffer. This is the memory space in which the collapse
    * procedure's double buffering will take place. Each buffer will be used
    * alternately as a read/write buffer when reducing the polylines of the
    * child nodes */
   if((res = darray_vertex_resize(vertex_list, vertices_range[1]+1)) != RES_OK
   || (res = darray_vertex_resize(scratch, nvertices)) != RES_OK) {
     ERROR(tree->sln, "Error in merging polylines -- %s.\n", res_to_cstr(res));
     goto error;
   }
 
   /* Recover the memory space of the tree's vertices */
   vertices = darray_vertex_data_get(vertex_list);
 
   /* Recover the memory space to be used in the Redux process.
    *
    * In the vertex buffer, this refers to the newly allocated memory space used
    * during the collapse. The beginning of the buffer contains the vertices of
    * the registered nodes and must therefore not be modified. At the end of the
    * collapse, this space will contain the vertices of the polylines resulting
    * from the collapse of the child nodes' polylines.
    *
    * The scratch buffer is used as-is, since its sole purpose is to temporarily
    * store the vertices of the polylines to be merged. */
   buf[0] = vertices + vertices_range[0];
   buf[1] = darray_vertex_data_get(scratch);
 
   /* Initially, the number of partitions to be reduced is equal to the number of
    * children */
   nparts = nchildren;
 
   /* Calculate the number of reduction steps, i.e., the logarithm of the power
    * of 2 that is greater than or equal to the number of polylines to be
    * reduced, to which one is added to obtain the number of steps, not the index
    * of the last step */
   ncollapses = log2i((int)round_up_pow2(nparts)) + 1;
 
   /* Set the index of the write buffer so that, once the collapse process is
    * complete, the last write buffer is the one whose memory space is allocated
    * in the tree's vertex buffer. This way, no additional copies will be needed
    * to store the result of the reduction */
   w = ncollapses % 2 ? 0 : 1;
 
   /* Initialize the polylines to be reduced, i.e., copy the child vertices into
    * a collapse buffer and record their index ranges */
   ivtx = 0;
   for(i = 0; i < nchildren; ++i) {
     const size_t ichild = inode + NODE(inode)->offset + i;
     const struct sln_node* child = NODE(ichild);
     const size_t memsz = child->nvertices * sizeof(struct sln_vertex);
     const struct sln_vertex* src = vertices + child->ivertex;
     struct sln_vertex* dst = buf[w] + ivtx;
 
     memcpy(dst, src, memsz);
     poly_parts[i][0] = ivtx;
     poly_parts[i][1] = ivtx + child->nvertices - 1/*inclusive bound*/;
 
     ivtx += child->nvertices;
   }
 
   r = w; /* Index of the buffer from which to read the data */
   w =!r; /* Index of the buffer into which to write the data  */
 
 
   /* As long as the number of segments is not equal to one, there are still
    * polylines to be merged in pairs */
   while(nparts > 1) {
     size_t ipart = 0;
 
     for(ipart=0; ipart < nparts-1; ipart+=2) {
       struct sln_mesh mesh0 = SLN_MESH_NULL;
       struct sln_mesh mesh1 = SLN_MESH_NULL;
 
       /* Retrieve the two polylines to be merged */
       const size_t* part0 = poly_parts[ipart+0];
       const size_t* part1 = poly_parts[ipart+1];
 
       /* Calculate the partition index of the merged polyline, that is, the index
        * that will contain the information for the polyline resulting from the
        * merge: the first and last indices of its in the vertex array currently
        * being written */
       const size_t ipart_merge = ipart/2;
 
       /* Compute the number of vertices of the merged polyline */
       const size_t nvtx =
         ((part0[1] - part0[0]) + 1/*inclusive bound*/)
       + ((part1[1] - part1[0]) + 1/*inclusive bound*/);
 
       /* For each child, initialized its vertex index to its first vertex. This
        * index will be incremented as vertices are merged into the merged
        * polyline  */
       size_t ivtx0 = part0[0];
       size_t ivtx1 = part1[0];
 
       /* Set the vertex index range for the merged polyline in the write buffer */
       range_merge[0] = part0[0];
       range_merge[1] = range_merge[0] + nvtx-1/*inclusive bound*/;
 
       /* Setup the mesh of the two polylines to be merged */
       mesh0.vertices = buf[r] + part0[0]; mesh0.nvertices = part0[1]-part0[0]+1;
       mesh1.vertices = buf[r] + part1[0]; mesh1.nvertices = part1[1]-part1[0]+1;
 
       /* Merge the polylines */
       FOR_EACH(ivtx, range_merge[0], range_merge[1]+1/*inclusive bound*/) {
         const double nu0 = ivtx0 <= part0[1] ? buf[r][ivtx0].wavenumber : INF;
         const double nu1 = ivtx1 <= part1[1] ? buf[r][ivtx1].wavenumber : INF;
         double ka = 0;
         double nu = 0;
 
         /* Find the minimum wave number among the vertices of the child vertices
          * that are candidates for merging */
         if(nu0 < nu1) { /* The vertex comes from the child0 */
           nu = buf[r][ivtx0].wavenumber;
           ka = buf[r][ivtx0].ka + sln_mesh_eval(&mesh1, nu);
           ++ivtx0;
 
         } else if(nu0 > nu1) { /* The vertex comes from the child1 */
           nu = buf[r][ivtx1].wavenumber;
           ka = buf[r][ivtx1].ka + sln_mesh_eval(&mesh0, nu);
           ++ivtx1;
 
         } else { /* The vertex is shared by node0 and node1 */
           nu = buf[r][ivtx0].wavenumber;
           ka = buf[r][ivtx0].ka + buf[r][ivtx1].ka;
           --range_merge[1]; /* Remove duplicate */
           ++ivtx0;
           ++ivtx1;
         }
         buf[w][ivtx].wavenumber = (float)nu;
         buf[w][ivtx].ka = (float)ka;
       }
 
       /* Decimate the resulting polyline */
       res = polyline_decimate(tree->sln, buf[w], range_merge,
         (float)tree->args.mesh_decimation_err, tree->args.mesh_type);
       if(res != RES_OK) goto error;
 
       /* Setup the partition of the merge polyline */
       poly_parts[ipart_merge][0] = range_merge[0];
       poly_parts[ipart_merge][1] = range_merge[1];
     }
 
     /* If there is a polyline that has not been merged, copy its vertices to the
      * write buffer so that it can be processed in the next reduction step */
     if(nparts % 2) {
       const size_t* remain_part = poly_parts[nparts-1];
       const size_t nvtx = remain_part[1] - remain_part[0] + 1/*inclusive*/;
       const size_t memsz = nvtx * sizeof(struct sln_vertex);
       const struct sln_vertex* src = buf[r] + remain_part[0];
       struct sln_vertex* dst = buf[w] + remain_part[0];
 
       memcpy(dst, src, memsz);
       poly_parts[nparts/2][0] = remain_part[0];
       poly_parts[nparts/2][1] = remain_part[1];
     }
 
 #ifndef NDEBUG
     FOR_EACH(i, 1, (nparts+1/*ceil*/)/2) {
       ASSERT(poly_parts[i][0] > poly_parts[i-1][1]);
     }
 #endif
 
     /* Update the number of partitions to be reduced */
     nparts = (nparts + 1/*ceil*/)/2;
 
     /* Swap read/write buffers */
     r = !r;
     w = !w;
   }
 
   nvertices = (range_merge[1] - range_merge[0]) + 1/*inclusive bound*/;
   vertices_range[1] = vertices_range[0] + nvertices - 1/*inclusive bound*/;
 
   /* Assumed that double buffering was configured to ensure that the resulting
    * polyline is stored in the tree vertex buffer */
   ASSERT(buf[r] != darray_vertex_cdata_get(scratch));
 
   /* Setup the node */
   NODE(inode)->ivertex = vertices_range[0];
   NODE(inode)->nvertices = ui64_to_ui32(nvertices);
 
   /* 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*KA_ADJUSTEMENT);
     }
   }
 
   /* Shrink the size of the vertices */
   darray_vertex_resize(vertex_list, vertices_range[1] + 1);
 
   #undef NCHILDREN
   #undef NODE
 
 exit:
   return res;
 error:
   goto exit;
 }
 static res_T
 build_polylines
   (struct sln_tree* tree,
    const size_t root_index,
    const size_t nodes_count, /* Total number of nodes in the tree (for debug) */
    struct darray_vertex* vertices,
    struct progress* progress)
 {
   /* Stack */
   #define STACK_SIZE (SLN_TREE_DEPTH_MAX*SLN_TREE_ARITY_MAX)
   size_t stack[STACK_SIZE];
   size_t istack = 0;
 
   /* Progress */
   size_t nnodes_processed = 0;
 
   /* Miscellaneous */
   struct build_leaf_polyline_scratch scratch;
   size_t inode = 0;
   res_T res = RES_OK;
 
   ASSERT(tree && nodes_count != 0);
   ASSERT(root_index < darray_node_size_get(&tree->nodes));
   (void)nodes_count; /* Avoid "Unused variable" warning */
 
   build_leaf_polyline_scratch_init(tree->sln->allocator, &scratch);
 
   #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 = root_index; /* Root node */
   while(inode != SIZE_MAX) {
     const size_t istack_saved = istack;
 
     if(IS_LEAF(inode)) {
       res = build_leaf_polyline(tree, NODE(inode), &scratch, vertices);
       if(res != RES_OK) goto error;
 
       inode = stack[--istack]; /* Pop the next node */
 
     } else {
       const size_t ichild0 = inode + NODE(inode)->offset + 0;
       const struct sln_node* node = darray_node_cdata_get(&tree->nodes)+inode;
       const unsigned nchildren = node_child_count(node, tree->args.arity);
       int child_polylines_are_missing = 1;
       size_t i = 0;
 
       FOR_EACH(i, 0, nchildren) {
         child_polylines_are_missing = NODE(ichild0 + i)->nvertices == 0;
         if(child_polylines_are_missing) break;
       }
 
       /* Child nodes have their polyline created */
       if(!child_polylines_are_missing) {
 #ifndef NDEBUG
         /* Check that all children have their polylines created */
         FOR_EACH(i, 1, nchildren) {
           const size_t ichild = ichild0 + i;
           ASSERT(NODE(ichild)->nvertices != 0);
         }
 #endif
         if(tree->args.collapse_polylines) {
           res = collapse_child_polylines(tree, inode, nchildren, &tree->nodes,
             &scratch.vertices_tmp, vertices);
         } else {
           res = merge_child_polylines
             (tree, inode, nchildren, &tree->nodes, vertices);
         }
         if(res != RES_OK) goto error;
 
         inode = stack[--istack]; /* Pop the next node */
 
       /* Child nodes have NOT their polyline created */
       } else {
         ASSERT(istack + (nchildren - 1/*ichild0*/ + 1/*inode*/) <= STACK_SIZE);
         stack[istack++] = inode; /* Push the current node */
 
          /* Push the child nodes, except for those whose polyline has already
           * been constructed, as is the case when the child node was the root
           * of a separately constructed subtree */
         FOR_EACH_REVERSE(i, nchildren, 0) {
           const size_t ichild = ichild0 + i-1;
           if(NODE(ichild)->nvertices == 0) stack[istack++] = ichild;
         }
 
         /* Ensure that at least one child was pushed */
         ASSERT(stack[istack-1] != inode);
 
         inode = stack[--istack]; /* Recursively build the polyline of the 1st child */
       }
     }
 
     /* Handle progression bar */
     if(istack < istack_saved) {
       size_t nnodes = istack_saved - istack;
 
       nnodes_processed += nnodes;
       progress_update(tree->sln, progress, nnodes);
     }
   }
   ASSERT(nnodes_processed == nodes_count);
 
   #undef NODE
   #undef IS_LEAF
   #undef LOG_MSG
   #undef STACK_SIZE
 
 exit:
   build_leaf_polyline_scratch_release(&scratch);
   return res;
 error:
   goto exit;
 }
 
 static res_T
 partition_lines_depth_first
   (struct sln_tree* tree,
    const size_t root_index,
    struct darray_node* nodes,
    struct progress* progress)
 {
   /* Stack */
   #define STACK_SIZE (SLN_TREE_DEPTH_MAX*(SLN_TREE_ARITY_MAX-1))
   size_t stack[STACK_SIZE];
   size_t istack = 0;
 
   /* Progress */
   size_t nlines_total = 0;
   size_t nlines_processed = 0;
 
   /* Miscellaneous */
   size_t inode = 0;
   res_T res = RES_OK;
 
   /* Pre-condition */
   ASSERT(tree && nodes && progress);
   ASSERT(root_index < darray_node_size_get(nodes));
   (void)nlines_total; /* Avoid "Unused variable" warning */
 
   #define NODE(Id) (darray_node_data_get(nodes) + (Id))
   #define CREATE_NODE {                                                        \
     res = darray_node_push_back(nodes, &SLN_NODE_NULL);                        \
     if(res != RES_OK) goto error;                                              \
   } (void)0
 
   nlines_total = NODE(root_index)->range[1] - NODE(root_index)->range[0] + 1;
   nlines_processed = 0;
 
   /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */
   stack[istack++] = SIZE_MAX;
 
   inode = root_index; /* Root node */
   while(inode != SIZE_MAX) {
     /* #lines into the node */
     size_t nlines_node = NODE(inode)->range[1] - NODE(inode)->range[0] + 1;
 
     /* Make a leaf */
     if(nlines_node <= tree->args.leaf_nlines) {
 
       NODE(inode)->offset = 0;
       inode = stack[--istack]; /* Pop the next node */
 
       nlines_processed += nlines_node; /* For debug */
 
       progress_update(tree->sln, progress, nlines_node);
 
     /* Split the node  */
     } else {
       size_t node_range[2] = {0,0};
       size_t nlines_child_min = 0; /* Min #lines per child */
       size_t nlines_remain = 0;
       size_t nchildren = 0;
       size_t iline = 0;
       size_t i = 0;
 
       /* Calculate the index of the first child */
       size_t ichildren = darray_node_size_get(nodes);
 
       /* Compute how the number of children that node has */
       nchildren = node_child_count(NODE(inode), tree->args.arity);
 
       ASSERT(nchildren <= tree->args.arity);
       ASSERT(ichildren > inode);
 
       node_range[0] = NODE(inode)->range[0];
       node_range[1] = NODE(inode)->range[1];
 
       /* Define the offset from the current node to its children */
       NODE(inode)->offset = ui64_to_ui32((uint64_t)(ichildren - inode));
 
       nlines_child_min = nlines_node / nchildren;
       nlines_remain = nlines_node % nchildren;
 
       iline = node_range[0];
       FOR_EACH(i, 0, nchildren) {
         /* Compute the number of lines per child. Start by assigning the minimum
          * number of lines to each child, then distribute the remaining lines
          * among the first children */
         size_t nlines_child = nlines_child_min + (i < nlines_remain);
 
         CREATE_NODE;
 
         /* Set the range of lines line for the newly created child. Note that
          * the boundaries of the range are inclusive, which is why 1 is
          * subtracted to the upper bound */
         NODE(ichildren+i)->range[0] = iline;
         NODE(ichildren+i)->range[1] = iline + nlines_child - 1/*inclusive bound*/;
         iline += nlines_child;
 
         /* Check that the child's lines are a subset of the parent's lines */
         ASSERT(NODE(ichildren+i)->range[0] >= node_range[0]);
         ASSERT(NODE(ichildren+i)->range[1] <= node_range[1]);
       }
 
       inode = ichildren; /* Make the first child the current node */
 
       /* Push the other children */
       ASSERT(istack + (nchildren-1/*1st child*/) <= STACK_SIZE);
       FOR_EACH_REVERSE(i, nchildren-1, 0) stack[istack++] = ichildren + i;
     }
   }
   ASSERT(nlines_processed == nlines_total);
 
   #undef NODE
   #undef CREATE_NODE
   #undef STACK_SIZE
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 static res_T
 partition_lines_breadth_first
   (struct sln_tree* tree,
    const size_t root_index,
    const unsigned nleaves_max_hint) /* Advice on the maxium of leafs to create */
 {
   /* Static memory space of the queue */
   struct item {
     struct list_node link;
     size_t inode;
   } items[QUEUE_SIZE_MAX];
 
   /* Linked lists for managing queue data */
   struct list_node free_items;
   struct list_node queue;
   size_t nnodes = 0; /* Number of enqueud nodes */
 
   /* Miscellaneous */
   size_t nlines_total = 0;
   size_t nleaves = 0;
   size_t i = 0;
   res_T res = RES_OK;
 
   /* Pre-conditions */
   ASSERT(tree);
   ASSERT(root_index < darray_node_size_get(&tree->nodes));
 
   /* Macros to help in managing nodes */
   #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id))
   #define CREATE_NODE {                                                        \
     res = darray_node_push_back(&tree->nodes, &SLN_NODE_NULL);                 \
     if(res != RES_OK) goto error;                                              \
   } (void)0
 
   /* Helper macro for the queue data structure */
   #define ENQUEUE(INode) {                                                     \
     struct item* item__ = NULL;                                                \
     ASSERT(!is_list_empty(&free_items));                                       \
     item__ = CONTAINER_OF(list_head(&free_items), struct item, link);          \
     item__->inode = INode;                                                     \
     list_move_tail(&item__->link, &queue);                                     \
     ++nnodes;                                                                  \
   }
   #define DEQUEUE                                                              \
     (list_move_tail(list_head(&queue), &free_items),                           \
     --nnodes,                                                                  \
     CONTAINER_OF(list_tail(&free_items), struct item, link)->inode)
 
   /* Setup the linked lists */
   list_init(&free_items);
   list_init(&queue);
   FOR_EACH(i, 0, sizeof(items)/sizeof(items[0])) {
     items[i].inode = SIZE_MAX;
     list_init(&items[i].link);
     list_add_tail(&free_items, &items[i].link);
   }
 
   SHTR(line_list_get_size(tree->args.lines, &nlines_total));
 
   /* Setup the root node */
   NODE(root_index)->range[0] = 0;
   NODE(root_index)->range[1] = nlines_total - 1;
 
   /* Register the root node in the queue */
   ENQUEUE(root_index);
 
   /* Breadth-first partitioning */
   while(!is_list_empty(&queue)) {
     size_t node_range[2] = {0,0};
     size_t nlines_node = 0; /* #lines in the node */
     size_t nlines_child_min = 0; /* Min #lines per child */
     size_t nlines_remain = 0; /* Lines to be distributed among the children */
     size_t nchildren = 0; /* #children of the node */
     size_t ichildren = 0; /* Index of the node's first child */
     size_t iline = 0;
     size_t inode = 0;
 
 
     /* Check whether the number of registered leaves and currently processed
      * nodes is greater than or equal to the recommended maximum number of
      * leaves. If so, stop breadth-first partitioning and treat all registered
      * nodes as leaves. The policy is therefore to accept more leaves than
      * expected */
     if(nnodes + nleaves >= nleaves_max_hint) {
       nleaves += nnodes;
       break;
     }
 
     inode = DEQUEUE; /* Get the current node */
 
     /* Retrieve the range of lines contained within the node */
     node_range[0] = NODE(inode)->range[0];
     node_range[1] = NODE(inode)->range[1];
     nlines_node = node_range[1] - node_range[0] + 1/*inclusive bound*/;
     ASSERT(nlines_node >= tree->args.leaf_nlines);
 
     if(nlines_node <= tree->args.leaf_nlines) {
       /* Make a leaf */
       ++nleaves;
       continue;
     }
 
     /* Compute the number of children that node has */
     nchildren = node_child_count(NODE(inode), tree->args.arity);
     ASSERT(nchildren <= tree->args.arity);
 
     /* Check whether, after adding the child nodes to the queue, the total
      * number of nodes in the queue, plus the number of leaves already created,
      * does not exceed the max size of the queue. If so, the breadth-first
      * partitioning is complete and the nodes currently in the queue are all
      * leaves */
     if(nnodes + nchildren + nleaves > QUEUE_SIZE_MAX) {
       nleaves += nnodes + 1/*Do not forget the the current node*/;
       break;
     }
 
     /* Calculate the index of the first child */
     ichildren = darray_node_size_get(&tree->nodes);
     ASSERT(ichildren > inode);
 
     /* Define the offset from the current node to its children */
     NODE(inode)->offset = ui64_to_ui32((uint64_t)(ichildren - inode));
 
     /* Calculate the minimum number of lines per child and the number of
      * remaining lines to be distributed equally among them */
     nlines_child_min = nlines_node / nchildren;
     nlines_remain = nlines_node % nchildren;
 
     iline = node_range[0];
     FOR_EACH(i, 0, nchildren) {
       /* Compute the number of lines per child. Start by assigning the minimum
        * number of lines to each child, then distribute the remaining lines
        * among the first children */
       size_t nlines_child = nlines_child_min + (i < nlines_remain);
 
       CREATE_NODE;
 
       /* Set the range of lines line for the newly created child. Note that
        * the boundaries of the range are inclusive, which is why 1 is
        * subtracted to the upper bound */
       NODE(ichildren+i)->range[0] = iline;
       NODE(ichildren+i)->range[1] = iline + nlines_child - 1/*inclusive bound*/;
       iline += nlines_child;
 
       /* Check that the child's lines are a subset of the parent's lines */
       ASSERT(NODE(ichildren+i)->range[0] >= node_range[0]);
       ASSERT(NODE(ichildren+i)->range[1] <= node_range[1]);
 
       ENQUEUE(ichildren+i); /* Register the child node */
     }
   }
 
   #undef NODE
   #undef CREATE_NODE
   #undef ENQUEUE
   #undef DEQUEUE
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 static res_T
 build_subtrees
   (struct sln_tree* tree,
    const size_t subtrees[],
    unsigned nsubtrees,
    struct progress* meshing)
 {
   /* Subtree temporary data */
   struct scratch {
     struct darray_node nodes;
     struct darray_vertex vertices;
     size_t nnodes;
   }* scratches;
 
   /* Miscellaneous */
   struct progress partitioning = PROGRESS_DEFAULT;
   struct mem_allocator* allocator = NULL;
   size_t nlines = 0;
   unsigned nthreads = 0;
   int i = 0;
   ATOMIC res = RES_OK;
 
   /* Pre-conditions */
   ASSERT(tree && subtrees && nsubtrees >= 1);
   ASSERT(nsubtrees < INT_MAX);
 
   #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id))
   #define IS_LEAF(Buf, Id) (darray_node_cdata_get(Buf)[Id].offset == 0)
 
   allocator = tree->sln->allocator;
 
   /* Allocate the per subtree scratch */
   scratches = MEM_CALLOC(allocator, nsubtrees, sizeof(*scratches));
   if(!scratches) { res = RES_MEM_ERR; goto error; }
   FOR_EACH(i, 0, (int)nsubtrees) {
     darray_node_init(allocator, &scratches[i].nodes);
     darray_vertex_init(allocator, &scratches[i].vertices);
   }
 
   /* Setup the progress bar and print that although nothing has been done yet,
    * the calculation is nevertheless in progress */
   SHTR(line_list_get_size(tree->args.lines, &nlines));
   partitioning.total = nlines;
   partitioning.msg = "Partitioning: ";
   progress_print(tree->sln, &partitioning);
 
   /* Set the number of threads to use. The advice on the number of threads must
    * not exceed the number of threads available on the machine (see the
    * sln_tree_create function); the caller can only specify a lower number of
    * threads. urthermore, the number of subtrees is calculated so that each subtree
    * corresponds to at least one thread, ensuring that there are no more
    * subtrees than available threads.
    *
    * The actual number of threads to be used for constructing the subtrees in
    * parallel should therefore always be set to the number of subtrees. However,
    * to provide greater flexibility in the distribution of threads or the number
    * of subtrees, the number of threads is calculated below as the minimum
    * between the number of available threads and the number of subtrees */
   nthreads = MMIN(tree->args.nthreads_hint, nsubtrees);
   omp_set_num_threads((int)nthreads);
 
   #pragma omp parallel for
   for(i=0; i < (int)nsubtrees; ++i) {
     struct sln_node root = SLN_NODE_NULL;
     struct scratch* scratch = scratches + i;
     size_t nnodes_s = 0;
     size_t nnodes_t = 0;
     res_T res2 = RES_OK;
 
     /* Skip the remaining subtrees in case of an error */
     if(ATOMIC_GET(&res) != RES_OK) continue;
 
     /* Copy the subtree root in the temporary node buffer */
     #pragma omp critical
     {
       root = *NODE(subtrees[i]);
     }
     res2 = darray_node_push_back(&scratch->nodes, &root);
     if(res2 != RES_OK) { ATOMIC_SET(&res, res2); continue; }
 
     /* Partition the line */
     res2 = partition_lines_depth_first
       (tree, 0/*root*/, &scratch->nodes, &partitioning);
     if(res2 != RES_OK) { ATOMIC_SET(&res, res2); continue; }
 
     #pragma omp critical
     {
       struct sln_node* dst = NULL;
       const struct sln_node* src = NULL;
 
       /* Increase the node buffer to store the subtree nodes */
       nnodes_s = darray_node_size_get(&scratch->nodes);
       nnodes_t = darray_node_size_get(&tree->nodes);
       res2 = darray_node_resize(&tree->nodes, nnodes_t + nnodes_s - 1/*root*/);
       if(res2 == RES_OK) {
 
         if(!IS_LEAF(&scratch->nodes, 0)) {
           /* Set the offset to the first child of the root node of the subtree */
           NODE(subtrees[i])->offset = ui64_to_ui32(nnodes_t - subtrees[i]);
         }
 
         /* Copy the subtree nodes in the main node buffer */
         src = darray_node_cdata_get(&scratch->nodes) + 1/*discard root*/;
         dst = darray_node_data_get(&tree->nodes) + nnodes_t;
         memcpy(dst, src, (nnodes_s-1/*discard root*/)*sizeof(*src));
       }
     }
     if(res2 != RES_OK) { ATOMIC_SET(&res, res2); continue; }
 
     /* Free the temporary buffer but retain the number of nodes in the subtree.
      * This number is used later when constructing polylines (see below) */
     darray_node_purge(&scratch->nodes);
     scratch->nnodes = nnodes_s;
   }
 
   /* Setup the progress bar and print that although nothing has been done yet,
    * the calculation is nevertheless in progress */
   *meshing = PROGRESS_DEFAULT;
   meshing->total = darray_node_size_get(&tree->nodes);
   meshing->msg = "Meshing: ";
   progress_print(tree->sln, meshing);
 
   #pragma omp parallel for
   for(i=0; i < (int)nsubtrees; ++i) {
     struct scratch* scratch = scratches + i;
     size_t inode = subtrees[i];
     size_t nverts_s = 0;
     size_t nverts_t = 0;
     size_t j = 0;
     res_T res2 = RES_OK;
 
     /* Skip the remaining subtrees in case of an error */
     if(ATOMIC_GET(&res) != RES_OK) continue;
 
     res2 = build_polylines
       (tree, inode, scratch->nnodes, &scratch->vertices, meshing);
     if(res2 != RES_OK) { ATOMIC_SET(&res, res2); continue; }
 
     #pragma omp critical
     {
       struct sln_vertex* dst = NULL;
       const struct sln_vertex* src = NULL;
 
       /* Increase the vertex buffer to store the subtree polylines */
       nverts_s = darray_vertex_size_get(&scratch->vertices);
       nverts_t = darray_vertex_size_get(&tree->vertices);
       res2 = darray_vertex_resize(&tree->vertices, nverts_t + nverts_s);
       if(res2 == RES_OK) {
         /* Copy the subtree nodes in the main vertex buffer */
         src = darray_vertex_cdata_get(&scratch->vertices);
         dst = darray_vertex_data_get(&tree->vertices) + nverts_t;
         memcpy(dst, src, nverts_s*sizeof(*src));
       }
     }
     if(res2 != RES_OK) { ATOMIC_SET(&res, res2); continue; }
 
     /* Update the index of the first vertex of the polyline for the nodes in the
      * subtree based on their new locations, i.e., the main vertex buffer */
     NODE(subtrees[i])->ivertex += nverts_t;
     FOR_EACH(j, 0, scratch->nnodes-1/*The root has been just treated*/) {
       inode = subtrees[i] + NODE(subtrees[i])->offset + j;
       NODE(inode)->ivertex += nverts_t;
     }
 
     /* Free the temporary vertex buffer */
     darray_node_purge(&scratch->nodes);
   }
 
   #undef NODE
   #undef IS_LEAF
 
 exit:
   /* Clean up temporary data */
   FOR_EACH(i, 0, (int)nsubtrees) {
     darray_node_release(&scratches[i].nodes);
     darray_vertex_release(&scratches[i].vertices);
   }
   MEM_RM(allocator, scratches);
   return (res_T)res;
 error:
   goto exit;
 }
 
 static res_T
 tree_build_sequential(struct sln_tree* tree)
 {
   /* Progress statuses */
   struct progress partitioning = PROGRESS_DEFAULT;
   struct progress meshing = PROGRESS_DEFAULT;
 
   struct sln_node node = SLN_NODE_NULL;
   size_t nlines = 0;
   size_t nnodes = 0;
   size_t iroot = 0;
   res_T res = RES_OK;
 
   SHTR(line_list_get_size(tree->args.lines, &nlines));
 
   iroot = darray_node_size_get(&tree->nodes);
   ASSERT(iroot == 0); /* assume that the tree contains no data */
 
   /* Setup the root node */
   node.range[0] = 0;
   node.range[1] = nlines - 1/* inclusive bound */;
   res = darray_node_push_back(&tree->nodes, &node);
   if(res != RES_OK) goto error;
 
   /* Partition lines */
   partitioning.total = nlines;
   partitioning.msg = "Partitioning: ";
   progress_print(tree->sln, &partitioning);
   res = partition_lines_depth_first(tree, iroot, &tree->nodes, &partitioning);
   if(res != RES_OK) goto error;
 
   /* Create polylines */
   nnodes = darray_node_size_get(&tree->nodes);
   meshing.total = nnodes;
   meshing.msg = "Meshing: ";
   progress_print(tree->sln, &meshing);
   res = build_polylines(tree, iroot, nnodes, &tree->vertices, &meshing);
   if(res != RES_OK) goto error;
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 static res_T
 tree_build_parallel(struct sln_tree* tree)
 {
   /* Stack data structure */
   #define STACK_SIZE (SLN_TREE_DEPTH_MAX*(SLN_TREE_ARITY_MAX-1))
   size_t stack[STACK_SIZE];
   size_t istack = 0;
 
   /* Subtrees */
   size_t subtrees[QUEUE_SIZE_MAX]; /* List of sub-tree root indices */
   unsigned nsubtrees = 0;
 
   /* Progress message of the meshing step */
   struct progress meshing = PROGRESS_DEFAULT;
 
   /* Miscellaneous */
   struct sln_node root = SLN_NODE_NULL;
   size_t nlines = 0;
   size_t nnodes_upper_tree = 0; /* #nodes from the root to the subtrees */
   size_t iroot = 0;
   size_t inode = 0;
   size_t i = 0;
   res_T res = RES_OK;
 
   ASSERT(tree);
 
   iroot = darray_node_size_get(&tree->nodes);
   ASSERT(iroot == 0); /* assume that the tree contains no data */
 
   SHTR(line_list_get_size(tree->args.lines, &nlines));
 
   /* Setup the root node */
   root.range[0] = 0;
   root.range[1] = nlines;
   res = darray_node_push_back(&tree->nodes, &root);
   if(res != RES_OK) goto error;
 
   /* Partition the lines in bread-first fixing the maximum number of leafs to
    * the available number of threads */
   res = partition_lines_breadth_first(tree, iroot, tree->args.nthreads_hint);
   if(res != RES_OK) goto error;
 
   /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */
   stack[istack++] = SIZE_MAX;
 
   /* Retrieve the leaf indices generated by bread-first partitioning. These
    * correspond to the roots of the subtree to be built */
   inode = iroot; /* Root node */
   while(inode != SIZE_MAX) {
     const struct sln_node* node = darray_node_cdata_get(&tree->nodes) + inode;
     ASSERT(inode < darray_node_size_get(&tree->nodes));
 
     if(sln_node_is_leaf(node)) {
       ASSERT(nsubtrees < QUEUE_SIZE_MAX);
       subtrees[nsubtrees++] = inode;
       inode = stack[--istack]; /* Pop the next node */
 
     } else {
       const size_t ichild0 = inode + node->offset;
       const unsigned nchildren = node_child_count(node, tree->args.arity);
 
       /* Push the children except the 1st one */
       ASSERT(istack + (nchildren-1/*ichild0*/) <= STACK_SIZE);
       FOR_EACH_REVERSE(i, nchildren-1, 0) stack[istack++] = ichild0 + i;
 
       /* Recursively traverse the 1st child */
       inode = ichild0;
     }
   }
 
   /* Retrieves the number of registered nodes. This is the number of nodes from
    * the root to the subtrees, excluding the roots of those subtrees */
   nnodes_upper_tree = darray_node_size_get(&tree->nodes);
   nnodes_upper_tree -= nsubtrees; /* Exclude the root node of the subtrees */
 
   res = build_subtrees(tree, subtrees, nsubtrees, &meshing);
   if(res != RES_OK) goto error;
 
   /* Build the polylines of the upper level if necessary */
   if(nnodes_upper_tree != 0) {
     res = build_polylines
       (tree, 0/*root*/, nnodes_upper_tree, &tree->vertices, &meshing);
     if(res != RES_OK) goto error;
   }
 
   #undef STACK_SIZE
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 /*******************************************************************************
  * Local functions
  ******************************************************************************/
 res_T
 tree_build(struct sln_tree* tree)
 {
   res_T res = RES_OK;
   ASSERT(tree);
 
   res = tree->args.nthreads_hint == 1
     ? tree_build_sequential(tree)
     : tree_build_parallel(tree);
   if(res != RES_OK) goto error;
 
 exit:
   return res;
 error:
   darray_node_purge(&tree->nodes);
   darray_vertex_purge(&tree->vertices);
   goto exit;
 }