suvm_voxelize.c (12295B)
1 /* Copyright (C) 2020-2023, 2025-2026 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #define _POSIX_C_SOURCE 200112L /* getopt/close support */ 17 18 #include "suvm.h" 19 #include <star/smsh.h> 20 #include <rsys/clock_time.h> 21 #include <rsys/cstr.h> 22 #include <rsys/mem_allocator.h> 23 24 #include <string.h> 25 #include <unistd.h> /* getopt & close functions */ 26 27 struct args { 28 const char* input_filename; 29 const char* output_filename; 30 unsigned def[3]; 31 int precompute_normals; 32 int quit; 33 int verbose; 34 }; 35 36 static const struct args ARGS_DEFAULT = { 37 NULL, /* Input file */ 38 NULL, /* Output file */ 39 {64, 64, 64}, /* Definition */ 40 0, /* Precompute normals */ 41 0, /* Quit */ 42 0 /* Verbose */ 43 }; 44 45 /******************************************************************************* 46 * Helper functions 47 ******************************************************************************/ 48 static INLINE void 49 usage(FILE* stream) 50 { 51 fprintf(stream, 52 "usage: suvm-voxelize [-hnv] [-d x,y,z] [-o output] [input]\n"); 53 } 54 55 static void 56 args_release(struct args* args) 57 { 58 ASSERT(args); 59 *args = ARGS_DEFAULT; 60 } 61 62 static res_T 63 args_init(struct args* args, const int argc, char** argv) 64 { 65 size_t n; 66 int opt; 67 res_T res = RES_OK; 68 ASSERT(args); 69 70 *args = ARGS_DEFAULT; 71 72 while((opt = getopt(argc, argv, "d:hno:v")) != -1) { 73 switch(opt) { 74 case 'd': 75 res = cstr_to_list_uint(optarg, ',', args->def, &n, 3); 76 if(res == RES_OK && n != 3) res = RES_BAD_ARG; 77 if(!args->def[0] || !args->def[1] || !args->def[2]) res = RES_BAD_ARG; 78 break; 79 case 'h': 80 usage(stdin); 81 args_release(args); 82 args->quit = 1; 83 break; 84 case 'n': 85 args->precompute_normals = 1; 86 break; 87 case 'o': 88 args->output_filename = optarg; 89 break; 90 case 'v': args->verbose = 1; break; 91 default: res = RES_BAD_ARG; break; 92 } 93 if(res != RES_OK) { 94 if(optarg) { 95 fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n", 96 argv[0], optarg, opt); 97 } 98 goto error; 99 } 100 } 101 102 if(optind < argc) { 103 args->input_filename = argv[optind]; 104 } 105 106 exit: 107 optind = 1; 108 return res; 109 error: 110 usage(stderr); 111 args_release(args); 112 goto exit; 113 } 114 115 static void 116 get_indices(const size_t itetra, size_t ids[4], void* ctx) 117 { 118 const uint64_t* ui64; 119 struct smsh_desc* desc = ctx; 120 ASSERT(desc && desc->dcell == 4); 121 ui64 = smsh_desc_get_cell(desc, itetra); 122 ids[0] = (size_t)ui64[0]; 123 ids[1] = (size_t)ui64[1]; 124 ids[2] = (size_t)ui64[2]; 125 ids[3] = (size_t)ui64[3]; 126 } 127 128 static void 129 get_position(const size_t ivert, double pos[3], void* ctx) 130 { 131 const double* dbl; 132 struct smsh_desc* desc = ctx; 133 ASSERT(desc && desc->dnode == 3); 134 dbl = smsh_desc_get_node(desc, ivert); 135 pos[0] = dbl[0]; 136 pos[1] = dbl[1]; 137 pos[2] = dbl[2]; 138 } 139 140 static void 141 write_voxels 142 (FILE* stream, 143 const int* vxls, 144 const unsigned def[3], 145 const double vxl_sz[3], 146 const double org[3]) 147 { 148 size_t ix, iy, iz; 149 size_t ivxl; 150 size_t i; 151 ASSERT(stream && vxls && def && def[0] && def[1] && def[2]); 152 153 fprintf(stream, "# vtk DataFile Version 2.0\n"); 154 fprintf(stream, "Volumetric grid\n"); 155 fprintf(stream, "ASCII\n"); 156 157 fprintf(stream, "DATASET RECTILINEAR_GRID\n"); 158 fprintf(stream, "DIMENSIONS %u %u %u\n", def[0]+1, def[1]+1, def[2]+1); 159 fprintf(stream, "X_COORDINATES %u float\n", def[0]+1); 160 FOR_EACH(i, 0, def[0]+1) fprintf(stream, "%g ", (double)i*vxl_sz[0] + org[0]); 161 fprintf(stream, "\n"); 162 fprintf(stream, "Y_COORDINATES %u float\n", def[1]+1); 163 FOR_EACH(i, 0, def[1]+1) fprintf(stream, "%g ", (double)i*vxl_sz[1] + org[1]); 164 fprintf(stream, "\n"); 165 fprintf(stream, "Z_COORDINATES %u float\n", def[2]+1); 166 FOR_EACH(i, 0, def[2]+1) fprintf(stream, "%g ", (double)i*vxl_sz[2] + org[2]); 167 fprintf(stream, "\n"); 168 169 fprintf(stream, "CELL_DATA %u\n", def[0]*def[1]*def[2]); 170 fprintf(stream, "SCALARS nintersections int 1\n"); 171 fprintf(stream, "LOOKUP_TABLE default\n"); 172 173 ivxl = 0; 174 FOR_EACH(iz, 0, def[2]) { 175 FOR_EACH(iy, 0, def[1]) { 176 FOR_EACH(ix, 0, def[0]) { 177 fprintf(stream, "%i\n", vxls[ivxl]); 178 ++ivxl; 179 } 180 } 181 } 182 } 183 184 static res_T 185 voxelize_suvm_volume 186 (const struct suvm_volume* vol, 187 const unsigned def[3], 188 int* vxls) 189 { 190 double low[3], upp[3]; 191 double vxl_sz[3]; 192 size_t iprim; 193 size_t nprims; 194 int progress = 0; 195 res_T res = RES_OK; 196 ASSERT(vol && def && def[0] && def[1] && def[2] && vxls); 197 198 /* Compute the voxel size */ 199 res = suvm_volume_get_aabb(vol, low, upp); 200 if(res != RES_OK) goto error; 201 vxl_sz[0] = (upp[0] - low[0]) / (double)def[0]; 202 vxl_sz[1] = (upp[1] - low[1]) / (double)def[1]; 203 vxl_sz[2] = (upp[2] - low[2]) / (double)def[2]; 204 205 res = suvm_volume_get_primitives_count(vol, &nprims); 206 if(res != RES_OK) goto error; 207 208 FOR_EACH(iprim, 0, nprims) { 209 struct suvm_primitive prim = SUVM_PRIMITIVE_NULL; 210 struct suvm_polyhedron poly; 211 size_t ivxl_low[3]; 212 size_t ivxl_upp[3]; 213 size_t ivxl[3]; 214 size_t i = 0; 215 int pcent; 216 217 CHK(suvm_volume_get_primitive(vol, iprim, &prim) == RES_OK); 218 CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK); 219 220 /* Transform the polyhedron AABB in voxel space */ 221 ivxl_low[0] = (size_t)((poly.lower[0] - low[0]) / vxl_sz[0]); 222 ivxl_low[1] = (size_t)((poly.lower[1] - low[1]) / vxl_sz[1]); 223 ivxl_low[2] = (size_t)((poly.lower[2] - low[2]) / vxl_sz[2]); 224 ivxl_upp[0] = (size_t)ceil((poly.upper[0] - low[0]) / vxl_sz[0]); 225 ivxl_upp[1] = (size_t)ceil((poly.upper[1] - low[1]) / vxl_sz[1]); 226 ivxl_upp[2] = (size_t)ceil((poly.upper[2] - low[2]) / vxl_sz[2]); 227 CHK(ivxl_low[0] < def[0] && ivxl_upp[0] <= def[0]); 228 CHK(ivxl_low[1] < def[1] && ivxl_upp[1] <= def[1]); 229 CHK(ivxl_low[2] < def[2] && ivxl_upp[2] <= def[2]); 230 231 FOR_EACH(ivxl[2], ivxl_low[2], ivxl_upp[2]) { 232 float vxl_low[3]; 233 float vxl_upp[3]; 234 vxl_low[0] = (float)low[0]; 235 vxl_low[1] = (float)low[1]; 236 vxl_low[2] = (float)((double)ivxl[2] * vxl_sz[2] + low[2]); 237 vxl_upp[0] = (float)upp[0]; 238 vxl_upp[1] = (float)upp[1]; 239 vxl_upp[2] = vxl_low[2] + (float)vxl_sz[2]; 240 241 FOR_EACH(ivxl[1], ivxl_low[1], ivxl_upp[1]) { 242 vxl_low[1] = (float)((double)ivxl[1] * vxl_sz[1] + low[1]); 243 vxl_upp[1] = vxl_low[1] + (float)vxl_sz[1]; 244 FOR_EACH(ivxl[0], ivxl_low[0], ivxl_upp[0]) { 245 vxl_low[0] = (float)((double)ivxl[0] * vxl_sz[0] + low[0]); 246 vxl_upp[0] = vxl_low[0] + (float)vxl_sz[0]; 247 248 i = ivxl[0] + ivxl[1]*def[0] + ivxl[2]*def[0]*def[1]; 249 vxls[i] += (int)suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp); 250 } 251 } 252 } 253 pcent = (int)((double)iprim * 100 / (double)nprims + 0.5/*round*/); 254 if(pcent > progress) { 255 progress = pcent; 256 fprintf(stderr, "Voxelizing: %3d%%\r", progress); 257 fflush(stderr); 258 } 259 } 260 261 fprintf(stderr, "Voxelizing: %3d%%\n", progress); 262 263 exit: 264 return res; 265 error: 266 goto exit; 267 } 268 269 /******************************************************************************* 270 * Main functions 271 ******************************************************************************/ 272 int 273 main(int argc, char** argv) 274 { 275 char dump[128]; 276 struct args args = ARGS_DEFAULT; 277 FILE* stream_out = stdout; 278 FILE* stream_in = stdin; 279 const char* stream_out_name = "stdout"; 280 const char* stream_in_name = "stdin"; 281 struct smsh* smsh = NULL; 282 struct smsh_create_args smsh_args = SMSH_CREATE_ARGS_DEFAULT; 283 struct smsh_load_stream_args load_stream_args = SMSH_LOAD_STREAM_ARGS_NULL; 284 struct smsh_desc desc = SMSH_DESC_NULL; 285 struct suvm_device* suvm = NULL; 286 struct suvm_volume* vol = NULL; 287 struct suvm_tetrahedral_mesh_args msh_args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL; 288 struct time t0, t1; 289 int* vxls = NULL; 290 double low[3]; 291 double upp[3]; 292 double vxl_sz[3]; 293 res_T res = RES_OK; 294 int err = 0; 295 296 res = args_init(&args, argc, argv); 297 if(res != RES_OK) goto error; 298 if(args.quit) goto exit; 299 300 /* Setup input stream */ 301 if(args.input_filename) { 302 stream_in = fopen(args.input_filename, "r"); 303 if(!stream_in) { 304 fprintf(stderr, "Could not open the input file '%s' -- %s.\n", 305 args.input_filename, strerror(errno)); 306 goto error; 307 } 308 stream_in_name = args.input_filename; 309 } 310 311 /* Setup output stream */ 312 if(args.output_filename) { 313 stream_out = fopen(args.output_filename, "w"); 314 if(!stream_out) { 315 fprintf(stderr, "Could not open the output file '%s' -- %s.\n", 316 args.output_filename, strerror(errno)); 317 goto error; 318 } 319 stream_out_name = args.output_filename; 320 (void)stream_out_name; 321 } 322 323 324 /* Load the submitted file */ 325 smsh_args.verbose = args.verbose; 326 res = smsh_create(&smsh_args, &smsh); 327 if(res != RES_OK) goto error; 328 load_stream_args.stream = stream_in; 329 load_stream_args.name = stream_in_name; 330 load_stream_args.memory_mapping = stream_in != stdin; 331 res = smsh_load_stream(smsh, &load_stream_args); 332 if(res != RES_OK) goto error; 333 res = smsh_get_desc(smsh, &desc); 334 if(res != RES_OK) goto error; 335 if(desc.dnode != 3 || desc.dcell != 4) { 336 fprintf(stderr, "Only tetrahedrical mesh are supported\n"); 337 res = RES_BAD_ARG; 338 goto error; 339 } 340 if(args.verbose) { 341 fprintf(stderr, "#nodes: %u; #tetrahedra: %u\n", desc.dnode, desc.dcell); 342 } 343 344 res = suvm_device_create(NULL, NULL, args.verbose, &suvm); 345 if(res != RES_OK) goto error; 346 347 msh_args.nvertices = desc.nnodes; 348 msh_args.ntetrahedra = desc.ncells; 349 msh_args.get_indices = get_indices; 350 msh_args.get_position = get_position; 351 msh_args.context = &desc; 352 msh_args.precompute_normals = args.precompute_normals; 353 354 /* Setup the volume from the loaded data */ 355 if(args.verbose) { 356 fprintf(stderr, "Partitionning the tetrahedral mesh '%s'\n", stream_in_name); 357 } 358 time_current(&t0); 359 res = suvm_tetrahedral_mesh_create(suvm, &msh_args, &vol); 360 if(res != RES_OK) goto error; 361 time_sub(&t0, time_current(&t1), &t0); 362 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 363 if(args.verbose) { 364 fprintf(stderr, "Build acceleration structure in %s\n", dump); 365 } 366 367 /* Allocate the list of voxels */ 368 vxls = mem_calloc(args.def[0]*args.def[1]*args.def[2], sizeof(*vxls)); 369 if(!vxls) { 370 fprintf(stderr, "Could not allocate the list of voxels.\n"); 371 res = RES_MEM_ERR; 372 goto error; 373 } 374 375 /* Voxelize the volume */ 376 time_current(&t0); 377 res = voxelize_suvm_volume(vol, args.def, vxls); 378 if(res != RES_OK) goto error; 379 time_sub(&t0, time_current(&t1), &t0); 380 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 381 if(args.verbose) { 382 fprintf(stderr, "Volume voxelized in %s\n", dump); 383 } 384 385 /* Write the voxels */ 386 if(args.verbose) { 387 fprintf(stderr, "Writing voxels in '%s'\n", stream_out_name); 388 } 389 SUVM(volume_get_aabb(vol, low, upp)); 390 vxl_sz[0] = (upp[0] - low[0]) / (double)args.def[0]; 391 vxl_sz[1] = (upp[1] - low[1]) / (double)args.def[1]; 392 vxl_sz[2] = (upp[2] - low[2]) / (double)args.def[2]; 393 time_current(&t0); 394 write_voxels(stream_out, vxls, args.def, vxl_sz, low); 395 time_sub(&t0, time_current(&t1), &t0); 396 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 397 if(args.verbose) { 398 fprintf(stderr, "Voxels written in %s\n", dump); 399 } 400 401 exit: 402 if(stream_out && stream_out != stdout) fclose(stream_out); 403 if(stream_in && stream_in != stdin) fclose(stream_in); 404 if(suvm) SUVM(device_ref_put(suvm)); 405 if(vol) SUVM(volume_ref_put(vol)); 406 if(smsh) SMSH(ref_put(smsh)); 407 if(vxls) mem_rm(vxls); 408 args_release(&args); 409 if(mem_allocated_size() != 0) { 410 fprintf(stderr, "Memory leaks: %lu Bytes.\n", 411 (unsigned long)mem_allocated_size()); 412 } 413 return err; 414 error: 415 err = -1; 416 goto exit; 417 }