suvm_voxelize.c (13146B)
1 /* Copyright (C) 2020-2023 |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 void 49 print_help(const char* cmd) 50 { 51 ASSERT(cmd); 52 printf( 53 "Usage: %s [options] [file]\n" 54 "Voxelize a smsh(5) file and save result following VTK file format.\n" 55 "With no file option, mesh is read from standard input\n", 56 cmd); 57 printf("\n"); 58 printf( 59 " -d X:Y:Z voxelisation definition along the 3 axes.\n" 60 " Default definition is %u:%u:%u.\n", 61 ARGS_DEFAULT.def[0], 62 ARGS_DEFAULT.def[1], 63 ARGS_DEFAULT.def[2]); 64 printf( 65 " -h display this help and exit.\n"); 66 printf( 67 " -n precompute the tetrahedra normals.\n"); 68 printf( 69 " -o <output> filename of the output VTK file. If not defined, write\n" 70 " results to standard ouput\n"); 71 printf( 72 " -v make the program verobse.\n"); 73 printf("\n"); 74 printf( 75 "This is free software released under the GNU GPL license, version 3 or\n" 76 "later. You are free to change or redistribute it under certain\n" 77 "conditions <http://gnu.org.licenses/gpl.html>\n"); 78 } 79 80 static void 81 args_release(struct args* args) 82 { 83 ASSERT(args); 84 *args = ARGS_DEFAULT; 85 } 86 87 static res_T 88 args_init(struct args* args, const int argc, char** argv) 89 { 90 size_t n; 91 int opt; 92 res_T res = RES_OK; 93 ASSERT(args); 94 95 *args = ARGS_DEFAULT; 96 97 while((opt = getopt(argc, argv, "d:hno:v")) != -1) { 98 switch(opt) { 99 case 'd': 100 res = cstr_to_list_uint(optarg, ':', args->def, &n, 3); 101 if(res == RES_OK && n != 3) res = RES_BAD_ARG; 102 if(!args->def[0] || !args->def[1] || !args->def[2]) res = RES_BAD_ARG; 103 break; 104 case 'h': 105 print_help(argv[0]); 106 args_release(args); 107 args->quit = 1; 108 break; 109 case 'n': 110 args->precompute_normals = 1; 111 break; 112 case 'o': 113 args->output_filename = optarg; 114 break; 115 case 'v': args->verbose = 1; break; 116 default: res = RES_BAD_ARG; break; 117 } 118 if(res != RES_OK) { 119 if(optarg) { 120 fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n", 121 argv[0], optarg, opt); 122 } 123 goto error; 124 } 125 } 126 127 if(optind < argc) { 128 args->input_filename = argv[optind]; 129 } 130 131 exit: 132 optind = 1; 133 return res; 134 error: 135 args_release(args); 136 goto exit; 137 } 138 139 static void 140 get_indices(const size_t itetra, size_t ids[4], void* ctx) 141 { 142 const uint64_t* ui64; 143 struct smsh_desc* desc = ctx; 144 ASSERT(desc && desc->dcell == 4); 145 ui64 = smsh_desc_get_cell(desc, itetra); 146 ids[0] = (size_t)ui64[0]; 147 ids[1] = (size_t)ui64[1]; 148 ids[2] = (size_t)ui64[2]; 149 ids[3] = (size_t)ui64[3]; 150 } 151 152 static void 153 get_position(const size_t ivert, double pos[3], void* ctx) 154 { 155 const double* dbl; 156 struct smsh_desc* desc = ctx; 157 ASSERT(desc && desc->dnode == 3); 158 dbl = smsh_desc_get_node(desc, ivert); 159 pos[0] = dbl[0]; 160 pos[1] = dbl[1]; 161 pos[2] = dbl[2]; 162 } 163 164 static void 165 write_voxels 166 (FILE* stream, 167 const int* vxls, 168 const unsigned def[3], 169 const double vxl_sz[3], 170 const double org[3]) 171 { 172 size_t ix, iy, iz; 173 size_t ivxl; 174 size_t i; 175 ASSERT(stream && vxls && def && def[0] && def[1] && def[2]); 176 177 fprintf(stream, "# vtk DataFile Version 2.0\n"); 178 fprintf(stream, "Volumetric grid\n"); 179 fprintf(stream, "ASCII\n"); 180 181 fprintf(stream, "DATASET RECTILINEAR_GRID\n"); 182 fprintf(stream, "DIMENSIONS %u %u %u\n", def[0]+1, def[1]+1, def[2]+1); 183 fprintf(stream, "X_COORDINATES %u float\n", def[0]+1); 184 FOR_EACH(i, 0, def[0]+1) fprintf(stream, "%g ", (double)i*vxl_sz[0] + org[0]); 185 fprintf(stream, "\n"); 186 fprintf(stream, "Y_COORDINATES %u float\n", def[1]+1); 187 FOR_EACH(i, 0, def[1]+1) fprintf(stream, "%g ", (double)i*vxl_sz[1] + org[1]); 188 fprintf(stream, "\n"); 189 fprintf(stream, "Z_COORDINATES %u float\n", def[2]+1); 190 FOR_EACH(i, 0, def[2]+1) fprintf(stream, "%g ", (double)i*vxl_sz[2] + org[2]); 191 fprintf(stream, "\n"); 192 193 fprintf(stream, "CELL_DATA %u\n", def[0]*def[1]*def[2]); 194 fprintf(stream, "SCALARS intersect_type int 1\n"); 195 fprintf(stream, "LOOKUP_TABLE default\n"); 196 197 ivxl = 0; 198 FOR_EACH(iz, 0, def[2]) { 199 FOR_EACH(iy, 0, def[1]) { 200 FOR_EACH(ix, 0, def[0]) { 201 fprintf(stream, "%i\n", vxls[ivxl]); 202 ++ivxl; 203 } 204 } 205 } 206 } 207 208 static res_T 209 voxelize_suvm_volume 210 (const struct suvm_volume* vol, 211 const unsigned def[3], 212 int* vxls) 213 { 214 double low[3], upp[3]; 215 double vxl_sz[3]; 216 size_t iprim; 217 size_t nprims; 218 int progress = 0; 219 res_T res = RES_OK; 220 ASSERT(vol && def && def[0] && def[1] && def[2] && vxls); 221 222 /* Compute the voxel size */ 223 res = suvm_volume_get_aabb(vol, low, upp); 224 if(res != RES_OK) goto error; 225 vxl_sz[0] = (upp[0] - low[0]) / (double)def[0]; 226 vxl_sz[1] = (upp[1] - low[1]) / (double)def[1]; 227 vxl_sz[2] = (upp[2] - low[2]) / (double)def[2]; 228 229 res = suvm_volume_get_primitives_count(vol, &nprims); 230 if(res != RES_OK) goto error; 231 232 FOR_EACH(iprim, 0, nprims) { 233 struct suvm_primitive prim = SUVM_PRIMITIVE_NULL; 234 struct suvm_polyhedron poly; 235 size_t ivxl_low[3]; 236 size_t ivxl_upp[3]; 237 size_t ivxl[3]; 238 size_t i = 0; 239 int pcent; 240 241 CHK(suvm_volume_get_primitive(vol, iprim, &prim) == RES_OK); 242 CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK); 243 244 /* Transform the polyhedron AABB in voxel space */ 245 ivxl_low[0] = (size_t)((poly.lower[0] - low[0]) / vxl_sz[0]); 246 ivxl_low[1] = (size_t)((poly.lower[1] - low[1]) / vxl_sz[1]); 247 ivxl_low[2] = (size_t)((poly.lower[2] - low[2]) / vxl_sz[2]); 248 ivxl_upp[0] = (size_t)ceil((poly.upper[0] - low[0]) / vxl_sz[0]); 249 ivxl_upp[1] = (size_t)ceil((poly.upper[1] - low[1]) / vxl_sz[1]); 250 ivxl_upp[2] = (size_t)ceil((poly.upper[2] - low[2]) / vxl_sz[2]); 251 CHK(ivxl_low[0] < def[0] && ivxl_upp[0] <= def[0]); 252 CHK(ivxl_low[1] < def[1] && ivxl_upp[1] <= def[1]); 253 CHK(ivxl_low[2] < def[2] && ivxl_upp[2] <= def[2]); 254 255 FOR_EACH(ivxl[2], ivxl_low[2], ivxl_upp[2]) { 256 float vxl_low[3]; 257 float vxl_upp[3]; 258 vxl_low[0] = (float)low[0]; 259 vxl_low[1] = (float)low[1]; 260 vxl_low[2] = (float)((double)ivxl[2] * vxl_sz[2] + low[2]); 261 vxl_upp[0] = (float)upp[0]; 262 vxl_upp[1] = (float)upp[1]; 263 vxl_upp[2] = vxl_low[2] + (float)vxl_sz[2]; 264 265 FOR_EACH(ivxl[1], ivxl_low[1], ivxl_upp[1]) { 266 vxl_low[1] = (float)((double)ivxl[1] * vxl_sz[1] + low[1]); 267 vxl_upp[1] = vxl_low[1] + (float)vxl_sz[1]; 268 FOR_EACH(ivxl[0], ivxl_low[0], ivxl_upp[0]) { 269 vxl_low[0] = (float)((double)ivxl[0] * vxl_sz[0] + low[0]); 270 vxl_upp[0] = vxl_low[0] + (float)vxl_sz[0]; 271 272 i = ivxl[0] + ivxl[1]*def[0] + ivxl[2]*def[0]*def[1]; 273 vxls[i] += (int)suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp); 274 } 275 } 276 } 277 pcent = (int)((double)iprim * 100 / (double)nprims + 0.5/*round*/); 278 if(pcent > progress) { 279 progress = pcent; 280 fprintf(stderr, "Voxelizing: %3d%%\r", progress); 281 fflush(stderr); 282 } 283 } 284 285 fprintf(stderr, "Voxelizing: %3d%%\n", progress); 286 287 exit: 288 return res; 289 error: 290 goto exit; 291 } 292 293 /******************************************************************************* 294 * Main functions 295 ******************************************************************************/ 296 int 297 main(int argc, char** argv) 298 { 299 char dump[128]; 300 struct args args = ARGS_DEFAULT; 301 FILE* stream_out = stdout; 302 FILE* stream_in = stdin; 303 const char* stream_out_name = "stdout"; 304 const char* stream_in_name = "stdin"; 305 struct smsh* smsh = NULL; 306 struct smsh_create_args smsh_args = SMSH_CREATE_ARGS_DEFAULT; 307 struct smsh_load_stream_args load_stream_args = SMSH_LOAD_STREAM_ARGS_NULL; 308 struct smsh_desc desc = SMSH_DESC_NULL; 309 struct suvm_device* suvm = NULL; 310 struct suvm_volume* vol = NULL; 311 struct suvm_tetrahedral_mesh_args msh_args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL; 312 struct time t0, t1; 313 int* vxls = NULL; 314 double low[3]; 315 double upp[3]; 316 double vxl_sz[3]; 317 res_T res = RES_OK; 318 int err = 0; 319 320 res = args_init(&args, argc, argv); 321 if(res != RES_OK) goto error; 322 if(args.quit) goto exit; 323 324 /* Setup input stream */ 325 if(args.input_filename) { 326 stream_in = fopen(args.input_filename, "r"); 327 if(!stream_in) { 328 fprintf(stderr, "Could not open the input file '%s' -- %s.\n", 329 args.input_filename, strerror(errno)); 330 goto error; 331 } 332 stream_in_name = args.input_filename; 333 } 334 335 /* Setup output stream */ 336 if(args.output_filename) { 337 stream_out = fopen(args.output_filename, "w"); 338 if(!stream_out) { 339 fprintf(stderr, "Could not open the output file '%s' -- %s.\n", 340 args.output_filename, strerror(errno)); 341 goto error; 342 } 343 stream_out_name = args.output_filename; 344 (void)stream_out_name; 345 } 346 347 348 /* Load the submitted file */ 349 smsh_args.verbose = args.verbose; 350 res = smsh_create(&smsh_args, &smsh); 351 if(res != RES_OK) goto error; 352 load_stream_args.stream = stream_in; 353 load_stream_args.name = stream_in_name; 354 load_stream_args.memory_mapping = stream_in != stdin; 355 res = smsh_load_stream(smsh, &load_stream_args); 356 if(res != RES_OK) goto error; 357 res = smsh_get_desc(smsh, &desc); 358 if(res != RES_OK) goto error; 359 if(desc.dnode != 3 || desc.dcell != 4) { 360 fprintf(stderr, "Only tetrahedrical mesh are supported\n"); 361 res = RES_BAD_ARG; 362 goto error; 363 } 364 if(args.verbose) { 365 fprintf(stderr, "#nodes: %u; #tetrahedra: %u\n", desc.dnode, desc.dcell); 366 } 367 368 res = suvm_device_create(NULL, NULL, args.verbose, &suvm); 369 if(res != RES_OK) goto error; 370 371 msh_args.nvertices = desc.nnodes; 372 msh_args.ntetrahedra = desc.ncells; 373 msh_args.get_indices = get_indices; 374 msh_args.get_position = get_position; 375 msh_args.context = &desc; 376 msh_args.precompute_normals = args.precompute_normals; 377 378 /* Setup the volume from the loaded data */ 379 if(args.verbose) { 380 fprintf(stderr, "Partitionning the tetrahedral mesh '%s'\n", stream_in_name); 381 } 382 time_current(&t0); 383 res = suvm_tetrahedral_mesh_create(suvm, &msh_args, &vol); 384 if(res != RES_OK) goto error; 385 time_sub(&t0, time_current(&t1), &t0); 386 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 387 if(args.verbose) { 388 fprintf(stderr, "Build acceleration structure in %s\n", dump); 389 } 390 391 /* Allocate the list of voxels */ 392 vxls = mem_calloc(args.def[0]*args.def[1]*args.def[2], sizeof(*vxls)); 393 if(!vxls) { 394 fprintf(stderr, "Could not allocate the list of voxels.\n"); 395 res = RES_MEM_ERR; 396 goto error; 397 } 398 399 /* Voxelize the volume */ 400 time_current(&t0); 401 res = voxelize_suvm_volume(vol, args.def, vxls); 402 if(res != RES_OK) goto error; 403 time_sub(&t0, time_current(&t1), &t0); 404 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 405 if(args.verbose) { 406 fprintf(stderr, "Volume voxelized in %s\n", dump); 407 } 408 409 /* Write the voxels */ 410 if(args.verbose) { 411 fprintf(stderr, "Writing voxels in '%s'\n", stream_out_name); 412 } 413 SUVM(volume_get_aabb(vol, low, upp)); 414 vxl_sz[0] = (upp[0] - low[0]) / (double)args.def[0]; 415 vxl_sz[1] = (upp[1] - low[1]) / (double)args.def[1]; 416 vxl_sz[2] = (upp[2] - low[2]) / (double)args.def[2]; 417 time_current(&t0); 418 write_voxels(stream_out, vxls, args.def, vxl_sz, low); 419 time_sub(&t0, time_current(&t1), &t0); 420 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 421 if(args.verbose) { 422 fprintf(stderr, "Voxels written in %s\n", dump); 423 } 424 425 exit: 426 if(stream_out && stream_out != stdout) fclose(stream_out); 427 if(stream_in && stream_in != stdin) fclose(stream_in); 428 if(suvm) SUVM(device_ref_put(suvm)); 429 if(vol) SUVM(volume_ref_put(vol)); 430 if(smsh) SMSH(ref_put(smsh)); 431 if(vxls) mem_rm(vxls); 432 args_release(&args); 433 if(mem_allocated_size() != 0) { 434 fprintf(stderr, "Memory leaks: %lu Bytes.\n", 435 (unsigned long)mem_allocated_size()); 436 } 437 return err; 438 error: 439 err = -1; 440 goto exit; 441 }