sln_tree_build.c (11394B)
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 size_t nnodes_total = 0; 180 size_t nnodes_processed = 0; 181 int progress = 0; 182 res_T res = RES_OK; 183 ASSERT(tree); 184 185 #define LOG_MSG "Meshing: %3d%%\n" 186 #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id)) 187 #define IS_LEAF(Id) (NODE(Id)->offset == 0) 188 189 nnodes_total = darray_node_size_get(&tree->nodes); 190 191 /* Print that nothing has been done yet */ 192 INFO(tree->sln, LOG_MSG, progress); 193 194 /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */ 195 stack[istack++] = SIZE_MAX; 196 197 inode = 0; /* Root node */ 198 while(inode != SIZE_MAX) { 199 const size_t istack_saved = istack; 200 201 if(IS_LEAF(inode)) { 202 res = build_leaf_polyline(tree, NODE(inode)); 203 if(res != RES_OK) goto error; 204 205 inode = stack[--istack]; /* Pop the next node */ 206 207 } else { 208 const size_t ichild0 = inode + NODE(inode)->offset + 0; 209 const size_t ichild1 = inode + NODE(inode)->offset + 1; 210 211 /* Child nodes have their polyline created */ 212 if(NODE(ichild0)->nvertices) { 213 ASSERT(NODE(ichild1)->nvertices != 0); 214 215 /* Build the polyline of the current node by merging the polylines of 216 * its children */ 217 res = merge_node_polylines 218 (tree, NODE(ichild0), NODE(ichild1), NODE(inode)); 219 if(res != RES_OK) goto error; 220 inode = stack[--istack]; 221 222 } else { 223 ASSERT(NODE(ichild1)->nvertices == 0); 224 225 ASSERT(istack < STACK_SIZE - 2); 226 stack[istack++] = inode; /* Push the current node */ 227 stack[istack++] = ichild1; /* Push child1 */ 228 229 /* Recursively build the polyline of child0 */ 230 inode = ichild0; 231 } 232 } 233 234 /* Handle progression bar */ 235 if(istack < istack_saved) { 236 int pcent = 0; 237 238 nnodes_processed += istack_saved - istack; 239 ASSERT(nnodes_processed <= nnodes_total); 240 241 /* Print progress update */ 242 pcent = (int)((double)nnodes_processed*100.0/(double)nnodes_total+0.5); 243 if(pcent/10 > progress/10) { 244 progress = pcent; 245 INFO(tree->sln, LOG_MSG, progress); 246 } 247 } 248 } 249 ASSERT(nnodes_processed == nnodes_total); 250 251 #undef NODE 252 #undef IS_LEAF 253 #undef LOG_MSG 254 255 exit: 256 return res; 257 error: 258 goto exit; 259 } 260 261 static res_T 262 partition_lines(struct sln_tree* tree) 263 { 264 size_t stack[STACK_SIZE]; 265 size_t istack = 0; 266 size_t inode = 0; 267 size_t nlines = 0; 268 size_t nlines_total = 0; 269 size_t nlines_processed = 0; 270 int progress = 0; 271 res_T res = RES_OK; 272 273 ASSERT(tree); 274 SHTR(line_list_get_size(tree->args.lines, &nlines)); 275 276 nlines_total = nlines; 277 nlines_processed = 0; 278 279 #define LOG_MSG "Partitioning: %3d%%\n" 280 #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id)) 281 #define CREATE_NODE { \ 282 res = darray_node_push_back(&tree->nodes, &SLN_NODE_NULL); \ 283 if(res != RES_OK) goto error; \ 284 } (void)0 285 286 CREATE_NODE; /* Alloc the root node */ 287 288 /* Setup the root node */ 289 NODE(0)->range[0] = 0; 290 NODE(0)->range[1] = nlines - 1; 291 292 /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */ 293 stack[istack++] = SIZE_MAX; 294 295 /* Print that although nothing has been done yet, 296 * the calculation is nevertheless in progress */ 297 INFO(tree->sln, LOG_MSG, progress); 298 299 inode = 0; /* Root node */ 300 while(inode != SIZE_MAX) { 301 /* #lines into the node */ 302 nlines = NODE(inode)->range[1] - NODE(inode)->range[0] + 1; 303 304 /* Make a leaf */ 305 if(nlines <= MAX_NLINES_PER_LEAF) { 306 int pcent = 0; 307 308 NODE(inode)->offset = 0; 309 inode = stack[--istack]; /* Pop the next node */ 310 311 ASSERT(nlines_processed + nlines <= nlines_total); 312 nlines_processed += nlines; 313 314 /* Print progress update */ 315 pcent = (int)((double)nlines_processed*100.0/(double)nlines_total+0.5); 316 if(pcent/10 > progress/10) { 317 progress = pcent; 318 INFO(tree->sln, LOG_MSG, progress); 319 } 320 321 /* Split the node */ 322 } else { 323 /* Median split */ 324 const size_t split = NODE(inode)->range[0] + (nlines + 1/*ceil*/)/2; 325 326 /* Compute the index toward the 2 children (first is the left child, 327 * followed by the right child) */ 328 size_t ichildren = darray_node_size_get(&tree->nodes); 329 330 ASSERT(NODE(inode)->range[0] <= split); 331 ASSERT(NODE(inode)->range[1] >= split); 332 ASSERT(ichildren > inode); 333 334 /* Define the offset from the current node to its children */ 335 NODE(inode)->offset = ui64_to_ui32((uint64_t)(ichildren - inode)); 336 337 CREATE_NODE; /* Alloc left child */ 338 CREATE_NODE; /* Alloc right child */ 339 340 /* Setup the left child */ 341 NODE(ichildren+0)->range[0] = NODE(inode)->range[0]; 342 NODE(ichildren+0)->range[1] = split-1; 343 /* Setup the right child */ 344 NODE(ichildren+1)->range[0] = split; 345 NODE(ichildren+1)->range[1] = NODE(inode)->range[1]; 346 347 inode = ichildren + 0; /* Make the left child current node */ 348 349 ASSERT(istack < STACK_SIZE - 1); 350 stack[istack++] = ichildren + 1; /* Push the right child */ 351 } 352 } 353 ASSERT(nlines_processed == nlines_total); 354 355 #undef NODE 356 #undef CREATE_NODE 357 358 exit: 359 return res; 360 error: 361 darray_node_clear(&tree->nodes); 362 goto exit; 363 } 364 365 /******************************************************************************* 366 * Local functions 367 ******************************************************************************/ 368 res_T 369 tree_build(struct sln_tree* tree) 370 { 371 res_T res = RES_OK; 372 373 if((res = partition_lines(tree)) != RES_OK) goto error; 374 if((res = build_polylines(tree)) != RES_OK) goto error; 375 376 exit: 377 return res; 378 error: 379 goto exit; 380 }