sln_tree_build.c (10050B)
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, 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 switch(tree->args.mesh_type) { 89 case SLN_MESH_UPPER: 90 SLN(node_get_mesh(tree, node, &mesh)); 91 ka = sln_mesh_eval(&mesh, wavenumber); 92 break; 93 case SLN_MESH_FIT: 94 /* We choose not to make any approximation when calculating the usual 95 * mesh. One thus evaluates ka on the lines of the node and not on the 96 * mesh of the node. */ 97 ka = sln_node_eval(tree, node, wavenumber); 98 break; 99 default: FATAL("Unreachable code.\n"); break; 100 } 101 return ka; 102 } 103 104 static res_T 105 merge_node_polylines 106 (struct sln_tree* tree, 107 const struct sln_node* node0, 108 const struct sln_node* node1, 109 struct sln_node* merged_node) 110 { 111 struct sln_vertex* vertices = NULL; 112 size_t vertices_range[2]; 113 size_t ivtx; 114 size_t ivtx0, ivtx0_max; 115 size_t ivtx1, ivtx1_max; 116 res_T res = RES_OK; 117 ASSERT(tree && node0 && node1 && merged_node); 118 119 /* Define the vertices range of the merged polyline */ 120 vertices_range[0] = darray_vertex_size_get(&tree->vertices); 121 vertices_range[1] = vertices_range[0] + node0->nvertices + node1->nvertices - 1; 122 123 /* Allocate the memory space to store the new polyline */ 124 res = darray_vertex_resize(&tree->vertices, vertices_range[1] + 1); 125 if(res != RES_OK) { 126 ERROR(tree->sln, "Error in merging polylines -- %s.\n", res_to_cstr(res)); 127 goto error; 128 } 129 vertices = darray_vertex_data_get(&tree->vertices); 130 131 ivtx0 = node0->ivertex; 132 ivtx1 = node1->ivertex; 133 ivtx0_max = ivtx0 + node0->nvertices - 1; 134 ivtx1_max = ivtx1 + node1->nvertices - 1; 135 FOR_EACH(ivtx, vertices_range[0], vertices_range[1]+1) { 136 const double nu0 = ivtx0 <= ivtx0_max ? vertices[ivtx0].wavenumber : INF; 137 const double nu1 = ivtx1 <= ivtx1_max ? vertices[ivtx1].wavenumber : INF; 138 float nu, ka; 139 140 if(nu0 < nu1) { 141 /* The vertex comes from the node0 */ 142 nu = vertices[ivtx0].wavenumber; 143 ka = (float)(vertices[ivtx0].ka + eval_ka(tree, node1, nu)); 144 ++ivtx0; 145 } else if(nu0 > nu1) { 146 /* The vertex comes from the node1 */ 147 nu = vertices[ivtx1].wavenumber; 148 ka = (float)(vertices[ivtx1].ka + eval_ka(tree, node0, nu)); 149 ++ivtx1; 150 } else { 151 /* The vertex is shared by node0 and node1 */ 152 nu = vertices[ivtx0].wavenumber; 153 ka = vertices[ivtx0].ka + vertices[ivtx1].ka; 154 --vertices_range[1]; /* Remove duplicate */ 155 ++ivtx0; 156 ++ivtx1; 157 } 158 vertices[ivtx].wavenumber = nu; 159 vertices[ivtx].ka = ka; 160 } 161 162 /* Decimate the resulting polyline */ 163 res = polyline_decimate(tree->sln, darray_vertex_data_get(&tree->vertices), 164 vertices_range, tree->args.mesh_decimation_err, tree->args.mesh_type); 165 if(res != RES_OK) goto error; 166 167 /* Shrink the size of the vertices */ 168 darray_vertex_resize(&tree->vertices, vertices_range[1] + 1); 169 170 /* Setup the node polyline */ 171 merged_node->ivertex = vertices_range[0]; 172 merged_node->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1); 173 174 exit: 175 return res; 176 error: 177 goto exit; 178 } 179 180 static res_T 181 build_polylines(struct sln_tree* tree) 182 { 183 size_t stack[STACK_SIZE*2]; 184 size_t istack = 0; 185 size_t inode = 0; 186 res_T res = RES_OK; 187 ASSERT(tree); 188 189 #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id)) 190 #define IS_LEAF(Id) (NODE(Id)->offset == 0) 191 192 /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */ 193 stack[istack++] = SIZE_MAX; 194 195 inode = 0; /* Root node */ 196 while(inode != SIZE_MAX) { 197 198 if(IS_LEAF(inode)) { 199 res = build_leaf_polyline(tree, NODE(inode)); 200 if(res != RES_OK) goto error; 201 202 inode = stack[--istack]; /* Pop the next node */ 203 204 } else { 205 const size_t ichild0 = inode + NODE(inode)->offset + 0; 206 const size_t ichild1 = inode + NODE(inode)->offset + 1; 207 208 /* Child nodes have their polyline created */ 209 if(NODE(ichild0)->nvertices) { 210 ASSERT(NODE(ichild1)->nvertices != 0); 211 212 /* Build the polyline of the current node by merging the polylines of 213 * its children */ 214 res = merge_node_polylines 215 (tree, NODE(ichild0), NODE(ichild1), NODE(inode)); 216 if(res != RES_OK) goto error; 217 inode = stack[--istack]; 218 219 } else { 220 ASSERT(NODE(ichild1)->nvertices == 0); 221 222 ASSERT(istack < STACK_SIZE - 2); 223 stack[istack++] = inode; /* Push the current node */ 224 stack[istack++] = ichild1; /* Push child1 */ 225 226 /* Recursively build the polyline of child0 */ 227 inode = ichild0; 228 } 229 } 230 } 231 232 #undef NODE 233 234 exit: 235 return res; 236 error: 237 goto exit; 238 } 239 240 static res_T 241 partition_lines(struct sln_tree* tree) 242 { 243 size_t stack[STACK_SIZE]; 244 size_t istack = 0; 245 size_t inode; 246 size_t nlines; 247 res_T res = RES_OK; 248 249 ASSERT(tree); 250 SHTR(line_list_get_size(tree->args.lines, &nlines)); 251 252 #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id)) 253 #define CREATE_NODE { \ 254 res = darray_node_push_back(&tree->nodes, &SLN_NODE_NULL); \ 255 if(res != RES_OK) goto error; \ 256 } (void)0 257 258 CREATE_NODE; /* Alloc the root node */ 259 260 /* Setup the root node */ 261 NODE(0)->range[0] = 0; 262 NODE(0)->range[1] = nlines - 1; 263 264 /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */ 265 stack[istack++] = SIZE_MAX; 266 267 inode = 0; /* Root node */ 268 while(inode != SIZE_MAX) { 269 /* #lines into the node */ 270 nlines = NODE(inode)->range[1] - NODE(inode)->range[0] + 1; 271 272 /* Make a leaf */ 273 if(nlines <= MAX_NLINES_PER_LEAF) { 274 NODE(inode)->offset = 0; 275 inode = stack[--istack]; /* Pop the next node */ 276 277 /* Split the node */ 278 } else { 279 /* Median split */ 280 const size_t split = NODE(inode)->range[0] + (nlines + 1/*ceil*/)/2; 281 282 /* Compute the index toward the 2 children (first is the left child, 283 * followed by the right child) */ 284 size_t ichildren = darray_node_size_get(&tree->nodes); 285 286 ASSERT(NODE(inode)->range[0] <= split); 287 ASSERT(NODE(inode)->range[1] >= split); 288 ASSERT(ichildren > inode); 289 290 /* Define the offset from the current node to its children */ 291 NODE(inode)->offset = ui64_to_ui32((uint64_t)(ichildren - inode)); 292 293 CREATE_NODE; /* Alloc left child */ 294 CREATE_NODE; /* Alloc right child */ 295 296 /* Setup the left child */ 297 NODE(ichildren+0)->range[0] = NODE(inode)->range[0]; 298 NODE(ichildren+0)->range[1] = split-1; 299 /* Setup the right child */ 300 NODE(ichildren+1)->range[0] = split; 301 NODE(ichildren+1)->range[1] = NODE(inode)->range[1]; 302 303 inode = ichildren + 0; /* Make the left child current node */ 304 305 ASSERT(istack < STACK_SIZE - 1); 306 stack[istack++] = ichildren + 1; /* Push the right child */ 307 } 308 } 309 #undef NODE 310 #undef CREATE_NODE 311 312 exit: 313 return res; 314 error: 315 darray_node_clear(&tree->nodes); 316 goto exit; 317 } 318 319 /******************************************************************************* 320 * Local functions 321 ******************************************************************************/ 322 res_T 323 tree_build(struct sln_tree* tree) 324 { 325 res_T res = RES_OK; 326 327 if((res = partition_lines(tree)) != RES_OK) goto error; 328 if((res = build_polylines(tree)) != RES_OK) goto error; 329 330 exit: 331 return res; 332 error: 333 goto exit; 334 }