sln_tree_build.c (12223B)
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 double ka; 132 float nu; 133 134 if(nu0 < nu1) { 135 /* The vertex comes from the node0 */ 136 nu = vertices[ivtx0].wavenumber; 137 ka = vertices[ivtx0].ka + eval_ka(tree, node1, nu); 138 ++ivtx0; 139 } else if(nu0 > nu1) { 140 /* The vertex comes from the node1 */ 141 nu = vertices[ivtx1].wavenumber; 142 ka = vertices[ivtx1].ka + eval_ka(tree, node0, nu); 143 ++ivtx1; 144 } else { 145 /* The vertex is shared by node0 and node1 */ 146 nu = vertices[ivtx0].wavenumber; 147 ka = vertices[ivtx0].ka + vertices[ivtx1].ka; 148 --vertices_range[1]; /* Remove duplicate */ 149 ++ivtx0; 150 ++ivtx1; 151 } 152 vertices[ivtx].wavenumber = nu; 153 vertices[ivtx].ka = (float)ka; 154 } 155 156 /* Decimate the resulting polyline */ 157 res = polyline_decimate(tree->sln, darray_vertex_data_get(&tree->vertices), 158 vertices_range, (float)tree->args.mesh_decimation_err, tree->args.mesh_type); 159 if(res != RES_OK) goto error; 160 161 /* Shrink the size of the vertices */ 162 darray_vertex_resize(&tree->vertices, vertices_range[1] + 1); 163 164 /* Setup the node polyline */ 165 merged_node->ivertex = vertices_range[0]; 166 merged_node->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1); 167 168 /* It is necessary to ensure that the recorded vertices define a polyline 169 * along which any value (calculated by linear interpolation) is well above 170 * the sum of the corresponding values of the polylines to be merged. However, 171 * although this is guaranteed by definition for the vertices of the polyline, 172 * numerical uncertainty may nevertheless introduce errors that violate this 173 * criterion. Hence the following adjustment, which slightly increases the ka 174 * of the mesh so as to guarantee this constraint between the mesh of a node 175 * and that of its children */ 176 if(tree->args.mesh_type == SLN_MESH_UPPER) { 177 FOR_EACH(ivtx, vertices_range[0], vertices_range[1]+1/*inclusive bound*/) { 178 const double ka = vertices[ivtx].ka; 179 vertices[ivtx].ka = (float)(ka + MMAX(ka*1e-2, 1e-6)); 180 } 181 } 182 183 exit: 184 return res; 185 error: 186 goto exit; 187 } 188 189 static res_T 190 build_polylines(struct sln_tree* tree) 191 { 192 size_t stack[STACK_SIZE*2]; 193 size_t istack = 0; 194 size_t inode = 0; 195 size_t nnodes_total = 0; 196 size_t nnodes_processed = 0; 197 int progress = 0; 198 res_T res = RES_OK; 199 ASSERT(tree); 200 201 #define LOG_MSG "Meshing: %3d%%\n" 202 #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id)) 203 #define IS_LEAF(Id) (NODE(Id)->offset == 0) 204 205 nnodes_total = darray_node_size_get(&tree->nodes); 206 207 /* Print that nothing has been done yet */ 208 INFO(tree->sln, LOG_MSG, progress); 209 210 /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */ 211 stack[istack++] = SIZE_MAX; 212 213 inode = 0; /* Root node */ 214 while(inode != SIZE_MAX) { 215 const size_t istack_saved = istack; 216 217 if(IS_LEAF(inode)) { 218 res = build_leaf_polyline(tree, NODE(inode)); 219 if(res != RES_OK) goto error; 220 221 inode = stack[--istack]; /* Pop the next node */ 222 223 } else { 224 const size_t ichild0 = inode + NODE(inode)->offset + 0; 225 const size_t ichild1 = inode + NODE(inode)->offset + 1; 226 227 /* Child nodes have their polyline created */ 228 if(NODE(ichild0)->nvertices) { 229 ASSERT(NODE(ichild1)->nvertices != 0); 230 231 /* Build the polyline of the current node by merging the polylines of 232 * its children */ 233 res = merge_node_polylines 234 (tree, NODE(ichild0), NODE(ichild1), NODE(inode)); 235 if(res != RES_OK) goto error; 236 inode = stack[--istack]; 237 238 } else { 239 ASSERT(NODE(ichild1)->nvertices == 0); 240 241 ASSERT(istack < STACK_SIZE - 2); 242 stack[istack++] = inode; /* Push the current node */ 243 stack[istack++] = ichild1; /* Push child1 */ 244 245 /* Recursively build the polyline of child0 */ 246 inode = ichild0; 247 } 248 } 249 250 /* Handle progression bar */ 251 if(istack < istack_saved) { 252 int pcent = 0; 253 254 nnodes_processed += istack_saved - istack; 255 ASSERT(nnodes_processed <= nnodes_total); 256 257 /* Print progress update */ 258 pcent = (int)((double)nnodes_processed*100.0/(double)nnodes_total+0.5); 259 if(pcent/10 > progress/10) { 260 progress = pcent; 261 INFO(tree->sln, LOG_MSG, progress); 262 } 263 } 264 } 265 ASSERT(nnodes_processed == nnodes_total); 266 267 #undef NODE 268 #undef IS_LEAF 269 #undef LOG_MSG 270 271 exit: 272 return res; 273 error: 274 goto exit; 275 } 276 277 static res_T 278 partition_lines(struct sln_tree* tree) 279 { 280 size_t stack[STACK_SIZE]; 281 size_t istack = 0; 282 size_t inode = 0; 283 size_t nlines = 0; 284 size_t nlines_total = 0; 285 size_t nlines_processed = 0; 286 int progress = 0; 287 res_T res = RES_OK; 288 289 ASSERT(tree); 290 SHTR(line_list_get_size(tree->args.lines, &nlines)); 291 292 nlines_total = nlines; 293 nlines_processed = 0; 294 295 #define LOG_MSG "Partitioning: %3d%%\n" 296 #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id)) 297 #define CREATE_NODE { \ 298 res = darray_node_push_back(&tree->nodes, &SLN_NODE_NULL); \ 299 if(res != RES_OK) goto error; \ 300 } (void)0 301 302 CREATE_NODE; /* Alloc the root node */ 303 304 /* Setup the root node */ 305 NODE(0)->range[0] = 0; 306 NODE(0)->range[1] = nlines - 1; 307 308 /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */ 309 stack[istack++] = SIZE_MAX; 310 311 /* Print that although nothing has been done yet, 312 * the calculation is nevertheless in progress */ 313 INFO(tree->sln, LOG_MSG, progress); 314 315 inode = 0; /* Root node */ 316 while(inode != SIZE_MAX) { 317 /* #lines into the node */ 318 nlines = NODE(inode)->range[1] - NODE(inode)->range[0] + 1; 319 320 /* Make a leaf */ 321 if(nlines <= MAX_NLINES_PER_LEAF) { 322 int pcent = 0; 323 324 NODE(inode)->offset = 0; 325 inode = stack[--istack]; /* Pop the next node */ 326 327 ASSERT(nlines_processed + nlines <= nlines_total); 328 nlines_processed += nlines; 329 330 /* Print progress update */ 331 pcent = (int)((double)nlines_processed*100.0/(double)nlines_total+0.5); 332 if(pcent/10 > progress/10) { 333 progress = pcent; 334 INFO(tree->sln, LOG_MSG, progress); 335 } 336 337 /* Split the node */ 338 } else { 339 /* Median split */ 340 const size_t split = NODE(inode)->range[0] + (nlines + 1/*ceil*/)/2; 341 342 /* Compute the index toward the 2 children (first is the left child, 343 * followed by the right child) */ 344 size_t ichildren = darray_node_size_get(&tree->nodes); 345 346 ASSERT(NODE(inode)->range[0] <= split); 347 ASSERT(NODE(inode)->range[1] >= split); 348 ASSERT(ichildren > inode); 349 350 /* Define the offset from the current node to its children */ 351 NODE(inode)->offset = ui64_to_ui32((uint64_t)(ichildren - inode)); 352 353 CREATE_NODE; /* Alloc left child */ 354 CREATE_NODE; /* Alloc right child */ 355 356 /* Setup the left child */ 357 NODE(ichildren+0)->range[0] = NODE(inode)->range[0]; 358 NODE(ichildren+0)->range[1] = split-1; 359 /* Setup the right child */ 360 NODE(ichildren+1)->range[0] = split; 361 NODE(ichildren+1)->range[1] = NODE(inode)->range[1]; 362 363 inode = ichildren + 0; /* Make the left child current node */ 364 365 ASSERT(istack < STACK_SIZE - 1); 366 stack[istack++] = ichildren + 1; /* Push the right child */ 367 } 368 } 369 ASSERT(nlines_processed == nlines_total); 370 371 #undef NODE 372 #undef CREATE_NODE 373 374 exit: 375 return res; 376 error: 377 darray_node_clear(&tree->nodes); 378 goto exit; 379 } 380 381 /******************************************************************************* 382 * Local functions 383 ******************************************************************************/ 384 res_T 385 tree_build(struct sln_tree* tree) 386 { 387 res_T res = RES_OK; 388 389 if((res = partition_lines(tree)) != RES_OK) goto error; 390 if((res = build_polylines(tree)) != RES_OK) goto error; 391 392 exit: 393 return res; 394 error: 395 goto exit; 396 }