sln_tree_build.c (9900B)
1 /* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2026 Université de Lorraine 3 * Copyright (C) 2022 Centre National de la Recherche Scientifique 4 * Copyright (C) 2022 Université Paul Sabatier 5 * 6 * This program is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 18 19 #include "sln.h" 20 #include "sln_device_c.h" 21 #include "sln_polyline.h" 22 #include "sln_tree_c.h" 23 24 #include <star/shtr.h> 25 26 #include <rsys/cstr.h> 27 28 /* Maximum number of lines per leaf */ 29 #define MAX_NLINES_PER_LEAF 1 30 31 /* STACK_SIZE is set to the maximum depth of the partitioning tree generated by 32 * the tree_build function. This is actually the maximum stack size where tree 33 * nodes are pushed by the recursive build process */ 34 #define STACK_SIZE 64 35 36 /******************************************************************************* 37 * Helper functions 38 ******************************************************************************/ 39 static FINLINE uint32_t 40 ui64_to_ui32(const uint64_t ui64) 41 { 42 if(ui64 > UINT32_MAX) 43 VFATAL("%s: overflow %lu.\n", ARG2(FUNC_NAME, ((unsigned long)ui64))); 44 return (uint32_t)ui64; 45 } 46 47 static INLINE res_T 48 build_leaf_polyline(struct sln_tree* tree, struct sln_node* leaf) 49 { 50 size_t vertices_range[2]; 51 res_T res = RES_OK; 52 53 /* Currently assume that we have only one line per leaf */ 54 ASSERT(leaf->range[0] == leaf->range[1]); 55 56 /* Line meshing */ 57 res = line_mesh(tree, leaf->range[0], tree->args.nvertices_hint, vertices_range); 58 if(res != RES_OK) goto error; 59 60 /* Decimate the line mesh */ 61 res = polyline_decimate(tree->sln, darray_vertex_data_get(&tree->vertices), 62 vertices_range, (float)tree->args.mesh_decimation_err, tree->args.mesh_type); 63 if(res != RES_OK) goto error; 64 65 /* Shrink the size of the vertices */ 66 darray_vertex_resize(&tree->vertices, vertices_range[1] + 1); 67 68 /* Setup the leaf polyline */ 69 leaf->ivertex = vertices_range[0]; 70 leaf->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1); 71 72 exit: 73 return res; 74 error: 75 goto exit; 76 } 77 78 static INLINE double 79 eval_ka 80 (const struct sln_tree* tree, 81 const struct sln_node* node, 82 const double wavenumber) 83 { 84 struct sln_mesh mesh = SLN_MESH_NULL; 85 double ka = 0; 86 ASSERT(tree && node); 87 88 /* Whether the mesh to be constructed corresponds to the spectrum or its upper 89 * limit, use the node mesh to calculate the value of ka at a given wave 90 * number. Calculating the value from the node lines would take far too long*/ 91 SLN(node_get_mesh(tree, node, &mesh)); 92 ka = sln_mesh_eval(&mesh, wavenumber); 93 94 return ka; 95 } 96 97 static res_T 98 merge_node_polylines 99 (struct sln_tree* tree, 100 const struct sln_node* node0, 101 const struct sln_node* node1, 102 struct sln_node* merged_node) 103 { 104 struct sln_vertex* vertices = NULL; 105 size_t vertices_range[2]; 106 size_t ivtx; 107 size_t ivtx0, ivtx0_max; 108 size_t ivtx1, ivtx1_max; 109 res_T res = RES_OK; 110 ASSERT(tree && node0 && node1 && merged_node); 111 112 /* Define the vertices range of the merged polyline */ 113 vertices_range[0] = darray_vertex_size_get(&tree->vertices); 114 vertices_range[1] = vertices_range[0] + node0->nvertices + node1->nvertices - 1; 115 116 /* Allocate the memory space to store the new polyline */ 117 res = darray_vertex_resize(&tree->vertices, vertices_range[1] + 1); 118 if(res != RES_OK) { 119 ERROR(tree->sln, "Error in merging polylines -- %s.\n", res_to_cstr(res)); 120 goto error; 121 } 122 vertices = darray_vertex_data_get(&tree->vertices); 123 124 ivtx0 = node0->ivertex; 125 ivtx1 = node1->ivertex; 126 ivtx0_max = ivtx0 + node0->nvertices - 1; 127 ivtx1_max = ivtx1 + node1->nvertices - 1; 128 FOR_EACH(ivtx, vertices_range[0], vertices_range[1]+1) { 129 const double nu0 = ivtx0 <= ivtx0_max ? vertices[ivtx0].wavenumber : INF; 130 const double nu1 = ivtx1 <= ivtx1_max ? vertices[ivtx1].wavenumber : INF; 131 float nu, ka; 132 133 if(nu0 < nu1) { 134 /* The vertex comes from the node0 */ 135 nu = vertices[ivtx0].wavenumber; 136 ka = (float)(vertices[ivtx0].ka + eval_ka(tree, node1, nu)); 137 ++ivtx0; 138 } else if(nu0 > nu1) { 139 /* The vertex comes from the node1 */ 140 nu = vertices[ivtx1].wavenumber; 141 ka = (float)(vertices[ivtx1].ka + eval_ka(tree, node0, nu)); 142 ++ivtx1; 143 } else { 144 /* The vertex is shared by node0 and node1 */ 145 nu = vertices[ivtx0].wavenumber; 146 ka = vertices[ivtx0].ka + vertices[ivtx1].ka; 147 --vertices_range[1]; /* Remove duplicate */ 148 ++ivtx0; 149 ++ivtx1; 150 } 151 vertices[ivtx].wavenumber = nu; 152 vertices[ivtx].ka = ka; 153 } 154 155 /* Decimate the resulting polyline */ 156 res = polyline_decimate(tree->sln, darray_vertex_data_get(&tree->vertices), 157 vertices_range, (float)tree->args.mesh_decimation_err, tree->args.mesh_type); 158 if(res != RES_OK) goto error; 159 160 /* Shrink the size of the vertices */ 161 darray_vertex_resize(&tree->vertices, vertices_range[1] + 1); 162 163 /* Setup the node polyline */ 164 merged_node->ivertex = vertices_range[0]; 165 merged_node->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1); 166 167 exit: 168 return res; 169 error: 170 goto exit; 171 } 172 173 static res_T 174 build_polylines(struct sln_tree* tree) 175 { 176 size_t stack[STACK_SIZE*2]; 177 size_t istack = 0; 178 size_t inode = 0; 179 res_T res = RES_OK; 180 ASSERT(tree); 181 182 #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id)) 183 #define IS_LEAF(Id) (NODE(Id)->offset == 0) 184 185 /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */ 186 stack[istack++] = SIZE_MAX; 187 188 inode = 0; /* Root node */ 189 while(inode != SIZE_MAX) { 190 191 if(IS_LEAF(inode)) { 192 res = build_leaf_polyline(tree, NODE(inode)); 193 if(res != RES_OK) goto error; 194 195 inode = stack[--istack]; /* Pop the next node */ 196 197 } else { 198 const size_t ichild0 = inode + NODE(inode)->offset + 0; 199 const size_t ichild1 = inode + NODE(inode)->offset + 1; 200 201 /* Child nodes have their polyline created */ 202 if(NODE(ichild0)->nvertices) { 203 ASSERT(NODE(ichild1)->nvertices != 0); 204 205 /* Build the polyline of the current node by merging the polylines of 206 * its children */ 207 res = merge_node_polylines 208 (tree, NODE(ichild0), NODE(ichild1), NODE(inode)); 209 if(res != RES_OK) goto error; 210 inode = stack[--istack]; 211 212 } else { 213 ASSERT(NODE(ichild1)->nvertices == 0); 214 215 ASSERT(istack < STACK_SIZE - 2); 216 stack[istack++] = inode; /* Push the current node */ 217 stack[istack++] = ichild1; /* Push child1 */ 218 219 /* Recursively build the polyline of child0 */ 220 inode = ichild0; 221 } 222 } 223 } 224 225 #undef NODE 226 227 exit: 228 return res; 229 error: 230 goto exit; 231 } 232 233 static res_T 234 partition_lines(struct sln_tree* tree) 235 { 236 size_t stack[STACK_SIZE]; 237 size_t istack = 0; 238 size_t inode; 239 size_t nlines; 240 res_T res = RES_OK; 241 242 ASSERT(tree); 243 SHTR(line_list_get_size(tree->args.lines, &nlines)); 244 245 #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id)) 246 #define CREATE_NODE { \ 247 res = darray_node_push_back(&tree->nodes, &SLN_NODE_NULL); \ 248 if(res != RES_OK) goto error; \ 249 } (void)0 250 251 CREATE_NODE; /* Alloc the root node */ 252 253 /* Setup the root node */ 254 NODE(0)->range[0] = 0; 255 NODE(0)->range[1] = nlines - 1; 256 257 /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */ 258 stack[istack++] = SIZE_MAX; 259 260 inode = 0; /* Root node */ 261 while(inode != SIZE_MAX) { 262 /* #lines into the node */ 263 nlines = NODE(inode)->range[1] - NODE(inode)->range[0] + 1; 264 265 /* Make a leaf */ 266 if(nlines <= MAX_NLINES_PER_LEAF) { 267 NODE(inode)->offset = 0; 268 inode = stack[--istack]; /* Pop the next node */ 269 270 /* Split the node */ 271 } else { 272 /* Median split */ 273 const size_t split = NODE(inode)->range[0] + (nlines + 1/*ceil*/)/2; 274 275 /* Compute the index toward the 2 children (first is the left child, 276 * followed by the right child) */ 277 size_t ichildren = darray_node_size_get(&tree->nodes); 278 279 ASSERT(NODE(inode)->range[0] <= split); 280 ASSERT(NODE(inode)->range[1] >= split); 281 ASSERT(ichildren > inode); 282 283 /* Define the offset from the current node to its children */ 284 NODE(inode)->offset = ui64_to_ui32((uint64_t)(ichildren - inode)); 285 286 CREATE_NODE; /* Alloc left child */ 287 CREATE_NODE; /* Alloc right child */ 288 289 /* Setup the left child */ 290 NODE(ichildren+0)->range[0] = NODE(inode)->range[0]; 291 NODE(ichildren+0)->range[1] = split-1; 292 /* Setup the right child */ 293 NODE(ichildren+1)->range[0] = split; 294 NODE(ichildren+1)->range[1] = NODE(inode)->range[1]; 295 296 inode = ichildren + 0; /* Make the left child current node */ 297 298 ASSERT(istack < STACK_SIZE - 1); 299 stack[istack++] = ichildren + 1; /* Push the right child */ 300 } 301 } 302 #undef NODE 303 #undef CREATE_NODE 304 305 exit: 306 return res; 307 error: 308 darray_node_clear(&tree->nodes); 309 goto exit; 310 } 311 312 /******************************************************************************* 313 * Local functions 314 ******************************************************************************/ 315 res_T 316 tree_build(struct sln_tree* tree) 317 { 318 res_T res = RES_OK; 319 320 if((res = partition_lines(tree)) != RES_OK) goto error; 321 if((res = build_polylines(tree)) != RES_OK) goto error; 322 323 exit: 324 return res; 325 error: 326 goto exit; 327 }