commit ca3057d87598661d13365de8b3c6ba36e684075a
parent 3bc34f91bad49c5ba8e3a660b62069b3d6e0c4a4
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 17 Dec 2020 19:03:28 +0100
Add a CLI used to voxelize volumetric mesh
Diffstat:
2 files changed, 422 insertions(+), 0 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -26,6 +26,7 @@ option(NO_TEST "Do not build tests" OFF)
find_package(Embree 3.6 REQUIRED)
find_package(RCMake 0.4 REQUIRED)
find_package(RSys 0.10 REQUIRED)
+find_package(StarTetraHedra)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
include(rcmake)
@@ -99,6 +100,11 @@ if(NOT NO_TEST)
new_test(test_suvm_primitive_intersection)
endif()
+if(StarTetraHedra_FOUND)
+ add_executable(suvm_voxelize ${SUVM_SOURCE_DIR}/suvm_voxelize_sth.c)
+ target_link_libraries(suvm_voxelize RSys StarTetraHedra suvm)
+endif()
+
################################################################################
# Define output & install directories
################################################################################
diff --git a/src/suvm_voxelize_sth.c b/src/suvm_voxelize_sth.c
@@ -0,0 +1,416 @@
+/* Copyright (C) 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* getopt/close support */
+
+#include "suvm.h"
+#include <star/sth.h>
+#include <rsys/clock_time.h>
+#include <rsys/cstr.h>
+#include <rsys/mem_allocator.h>
+
+#include <string.h>
+#include <unistd.h> /* getopt & close functions */
+
+struct args {
+ const char* input_filename;
+ const char* output_filename;
+ unsigned def[3];
+ int quit;
+ int verbose;
+};
+
+static const struct args ARGS_DEFAULT = {
+ NULL, /* Input file */
+ NULL, /* Output file */
+ {64, 64, 64}, /* Definition */
+ 0, /* Quit */
+ 0 /* Verbose */
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+print_help(const char* cmd)
+{
+ ASSERT(cmd);
+ printf(
+"Usage: %s [options] [file]\n"
+"Voxelize the submitted Star-TetraHedra geometry and save the result\n"
+"into a VTK file.\n",
+ cmd);
+ printf("\n");
+ printf(
+" -d X:Y:Z voxelisation definition along the 3 axes.\n"
+" The default definition is %u:%u:%u.\n",
+ ARGS_DEFAULT.def[0],
+ ARGS_DEFAULT.def[1],
+ ARGS_DEFAULT.def[2]);
+ printf(
+" -h display this help and exit.\n");
+ printf(
+" -o <output> filename of the output VTK file. If not defined, write\n"
+" results to standard ouput\n");
+ printf(
+" -v make the program verobse.\n");
+ printf("\n");
+ printf(
+"Copyright (C) 2020 |Meso|Star> <contact@meso-star.com>.\n"
+"This is free software released under the GNU GPL license, version 3 or\n"
+"later. You are free to change or redistribute it under certain\n"
+"conditions <http://gnu.org.licenses/gpl.html>\n");
+}
+
+static void
+args_release(struct args* args)
+{
+ ASSERT(args);
+ *args = ARGS_DEFAULT;
+}
+
+static res_T
+args_init(struct args* args, const int argc, char** argv)
+{
+ size_t n;
+ int opt;
+ res_T res = RES_OK;
+ ASSERT(args);
+
+ *args = ARGS_DEFAULT;
+
+ while((opt = getopt(argc, argv, "d:ho:v")) != -1) {
+ switch(opt) {
+ case 'd':
+ res = cstr_to_list_uint(optarg, ':', args->def, &n, 3);
+ if(res == RES_OK && n != 3) res = RES_BAD_ARG;
+ if(!args->def[0] || !args->def[1] || !args->def[2]) res = RES_BAD_ARG;
+ break;
+ case 'h':
+ print_help(argv[0]);
+ args_release(args);
+ args->quit = 1;
+ break;
+ case 'o':
+ args->output_filename = optarg;
+ break;
+ case 'v': args->verbose = 1; break;
+ default: res = RES_BAD_ARG; break;
+ }
+ if(res != RES_OK) {
+ if(optarg) {
+ fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
+ argv[0], optarg, opt);
+ }
+ goto error;
+ }
+ }
+
+ if(optind < argc) {
+ args->input_filename = argv[optind];
+ }
+
+exit:
+ optind = 1;
+ return res;
+error:
+ args_release(args);
+ goto exit;
+}
+
+static void
+get_indices(const size_t itetra, size_t ids[4], void* ctx)
+{
+ const uint64_t* ui64;
+ struct sth_desc* desc = ctx;
+ ASSERT(desc);
+ ui64 = sth_desc_get_tetrahedron_indices(desc, itetra);
+ ids[0] = (size_t)ui64[0];
+ ids[1] = (size_t)ui64[1];
+ ids[2] = (size_t)ui64[2];
+ ids[3] = (size_t)ui64[3];
+}
+
+static void
+get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ const double* dbl;
+ struct sth_desc* desc = ctx;
+ ASSERT(desc);
+ dbl = sth_desc_get_vertex_position(desc, ivert);
+ pos[0] = dbl[0];
+ pos[1] = dbl[1];
+ pos[2] = dbl[2];
+}
+
+static void
+write_voxels
+ (FILE* stream,
+ const int* vxls,
+ const unsigned def[3],
+ const double vxl_sz[3],
+ const double org[3])
+{
+ size_t ix, iy, iz;
+ size_t ivxl;
+ size_t i;
+ ASSERT(stream && vxls && def && def[0] && def[1] && def[2]);
+
+ fprintf(stream, "# vtk DataFile Version 2.0\n");
+ fprintf(stream, "Volumetric grid\n");
+ fprintf(stream, "ASCII\n");
+
+ fprintf(stream, "DATASET RECTILINEAR_GRID\n");
+ fprintf(stream, "DIMENSIONS %u %u %u\n", def[0]+1, def[1]+1, def[2]+1);
+ fprintf(stream, "X_COORDINATES %u float\n", def[0]+1);
+ FOR_EACH(i, 0, def[0]+1) fprintf(stream, "%g ", (double)i*vxl_sz[0] + org[0]);
+ fprintf(stream, "\n");
+ fprintf(stream, "Y_COORDINATES %u float\n", def[1]+1);
+ FOR_EACH(i, 0, def[1]+1) fprintf(stream, "%g ", (double)i*vxl_sz[1] + org[1]);
+ fprintf(stream, "\n");
+ fprintf(stream, "Z_COORDINATES %u float\n", def[2]+1);
+ FOR_EACH(i, 0, def[2]+1) fprintf(stream, "%g ", (double)i*vxl_sz[2] + org[2]);
+ fprintf(stream, "\n");
+
+ fprintf(stream, "CELL_DATA %u\n", def[0]*def[1]*def[2]);
+ fprintf(stream, "SCALARS intersect_type int 1\n");
+ fprintf(stream, "LOOKUP_TABLE default\n");
+
+ ivxl = 0;
+ FOR_EACH(iz, 0, def[2]) {
+ FOR_EACH(iy, 0, def[1]) {
+ FOR_EACH(ix, 0, def[0]) {
+ fprintf(stream, "%i\n", vxls[ivxl]);
+ ++ivxl;
+ }
+ }
+ }
+}
+
+static res_T
+voxelize_suvm_volume
+ (const struct suvm_volume* vol,
+ const unsigned def[3],
+ int* vxls)
+{
+ double low[3], upp[3];
+ double vxl_sz[3];
+ size_t iprim;
+ size_t nprims;
+ res_T res = RES_OK;
+ ASSERT(vol && def && def[0] && def[1] && def[2] && vxls);
+
+ /* Compute the voxel size */
+ res = suvm_volume_get_aabb(vol, low, upp);
+ if(res != RES_OK) goto error;
+ vxl_sz[0] = (upp[0] - low[0]) / (double)def[0];
+ vxl_sz[1] = (upp[1] - low[1]) / (double)def[1];
+ vxl_sz[2] = (upp[2] - low[2]) / (double)def[2];
+
+ res = suvm_volume_get_primitives_count(vol, &nprims);
+ if(res != RES_OK) goto error;
+
+ FOR_EACH(iprim, 0, nprims) {
+ struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
+ struct suvm_polyhedron poly;
+ size_t ivxl_low[3];
+ size_t ivxl_upp[3];
+ size_t ivxl[3];
+ size_t i = 0;
+
+ CHK(suvm_volume_get_primitive(vol, iprim, &prim) == RES_OK);
+ CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK);
+
+ /* Transform the polyhedron AABB in voxel space */
+ ivxl_low[0] = (size_t)((poly.lower[0] - low[0]) / vxl_sz[0]);
+ ivxl_low[1] = (size_t)((poly.lower[1] - low[1]) / vxl_sz[1]);
+ ivxl_low[2] = (size_t)((poly.lower[2] - low[2]) / vxl_sz[2]);
+ ivxl_upp[0] = (size_t)ceil((poly.upper[0] - low[0]) / vxl_sz[0]);
+ ivxl_upp[1] = (size_t)ceil((poly.upper[1] - low[1]) / vxl_sz[1]);
+ ivxl_upp[2] = (size_t)ceil((poly.upper[2] - low[2]) / vxl_sz[2]);
+ CHK(ivxl_low[0] < def[0] && ivxl_upp[0] <= def[0]);
+ CHK(ivxl_low[1] < def[1] && ivxl_upp[1] <= def[1]);
+ CHK(ivxl_low[2] < def[2] && ivxl_upp[2] <= def[2]);
+
+ FOR_EACH(ivxl[2], ivxl_low[2], ivxl_upp[2]) {
+ float vxl_low[3];
+ float vxl_upp[3];
+ vxl_low[0] = (float)low[0];
+ vxl_low[1] = (float)low[1];
+ vxl_low[2] = (float)((double)ivxl[2] * vxl_sz[2] + low[2]);
+ vxl_upp[0] = (float)upp[0];
+ vxl_upp[1] = (float)upp[1];
+ vxl_upp[2] = vxl_low[2] + (float)vxl_sz[2];
+
+ FOR_EACH(ivxl[1], ivxl_low[1], ivxl_upp[1]) {
+ vxl_low[1] = (float)((double)ivxl[1] * vxl_sz[1] + low[1]);
+ vxl_upp[1] = vxl_low[1] + (float)vxl_sz[1];
+ FOR_EACH(ivxl[0], ivxl_low[0], ivxl_upp[0]) {
+ vxl_low[0] = (float)((double)ivxl[0] * vxl_sz[0] + low[0]);
+ vxl_upp[0] = vxl_low[0] + (float)vxl_sz[0];
+
+ i = ivxl[0] + ivxl[1]*def[0] + ivxl[2]*def[0]*def[1];
+ vxls[i] += (int)suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp);
+ }
+ }
+ }
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Main functions
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ char dump[128];
+ struct args args = ARGS_DEFAULT;
+ FILE* stream_out = stdout;
+ FILE* stream_in = stdin;
+ const char* stream_out_name = "stdout";
+ const char* stream_in_name = "stdin";
+ struct sth* sth = NULL;
+ struct sth_desc desc = STH_DESC_NULL;
+ struct suvm_device* suvm = NULL;
+ struct suvm_volume* vol = NULL;
+ struct suvm_tetrahedral_mesh_args msh_args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL;
+ struct time t0, t1;
+ int* vxls = NULL;
+ double low[3];
+ double upp[3];
+ double vxl_sz[3];
+ res_T res = RES_OK;
+ int err = 0;
+
+ res = args_init(&args, argc, argv);
+ if(res != RES_OK) goto error;
+ if(args.quit) goto exit;
+
+ /* Setup input stream */
+ if(args.input_filename) {
+ stream_in = fopen(args.input_filename, "r");
+ if(!stream_in) {
+ fprintf(stderr, "Could not open the input file '%s' -- %s.\n",
+ args.input_filename, strerror(errno));
+ goto error;
+ }
+ stream_in_name = args.input_filename;
+ }
+
+ /* Setup output stream */
+ if(args.output_filename) {
+ stream_out = fopen(args.output_filename, "w");
+ if(!stream_out) {
+ fprintf(stderr, "Could not open the output file '%s' -- %s.\n",
+ args.output_filename, strerror(errno));
+ goto error;
+ }
+ stream_out_name = args.output_filename;
+ (void)stream_out_name;
+ }
+
+
+ /* Load the submitted file */
+ res = sth_create(NULL, NULL, args.verbose, &sth);
+ if(res != RES_OK) goto error;
+ res = sth_load_stream(sth, stream_in, stream_in_name);
+ if(res != RES_OK) goto error;
+ res = sth_get_desc(sth, &desc);
+ if(res != RES_OK) goto error;
+ if(args.verbose) {
+ printf("#nodes: %lu; #tetrahedra: %lu\n", desc.nvertices, desc.ntetrahedra);
+ }
+
+ res = suvm_device_create(NULL, NULL, args.verbose, &suvm);
+ if(res != RES_OK) goto error;
+
+ msh_args.nvertices = desc.nvertices;
+ msh_args.ntetrahedra = desc.ntetrahedra;
+ msh_args.get_indices = get_indices;
+ msh_args.get_position = get_position;
+ msh_args.context = &desc;
+
+ /* Setup the volume from the loaded data */
+ if(args.verbose) {
+ printf("Partitionning the tetrahedral mesh '%s'\n", stream_in_name);
+ }
+ time_current(&t0);
+ res = suvm_tetrahedral_mesh_create(suvm, &msh_args, &vol);
+ if(res != RES_OK) goto error;
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ if(args.verbose) {
+ printf("Build acceleration structure in %s\n", dump);
+ }
+
+ /* Allocate the list of voxels */
+ vxls = mem_calloc(args.def[0]*args.def[1]*args.def[2], sizeof(*vxls));
+ if(!vxls) {
+ fprintf(stderr, "Could not allocate the list of voxels.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ /* Voxelize the volume */
+ if(args.verbose) {
+ printf("Voxelizing the volume\n");
+ }
+ time_current(&t0);
+ res = voxelize_suvm_volume(vol, args.def, vxls);
+ if(res != RES_OK) goto error;
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ if(args.verbose) {
+ printf("Volume voxelized in %s\n", dump);
+ }
+
+ /* Write the voxels */
+ if(args.verbose) {
+ printf("Writing voxels in '%s'\n", stream_out_name);
+ }
+ SUVM(volume_get_aabb(vol, low, upp));
+ vxl_sz[0] = (upp[0] - low[0]) / (double)args.def[0];
+ vxl_sz[1] = (upp[1] - low[1]) / (double)args.def[1];
+ vxl_sz[2] = (upp[2] - low[2]) / (double)args.def[2];
+ time_current(&t0);
+ write_voxels(stream_out, vxls, args.def, vxl_sz, low);
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ if(args.verbose) {
+ printf("Voxels written in %s\n", dump);
+ }
+
+exit:
+ if(stream_out && stream_out != stdout) fclose(stream_out);
+ if(stream_in && stream_in != stdin) fclose(stream_in);
+ if(suvm) SUVM(device_ref_put(suvm));
+ if(vol) SUVM(volume_ref_put(vol));
+ if(sth) STH(ref_put(sth));
+ if(vxls) mem_rm(vxls);
+ args_release(&args);
+ if(mem_allocated_size() != 0) {
+ fprintf(stderr, "Memory leaks: %lu Bytes.\n",
+ (unsigned long)mem_allocated_size());
+ }
+ return err;
+error:
+ err = -1;
+ goto exit;
+}
