From 292dcac9c9cdf2d10e4c75881ef3a74a0bf69b2d Mon Sep 17 00:00:00 2001 From: Rafael Sachetto Date: Thu, 27 Jun 2024 11:59:17 -0300 Subject: [PATCH] More progress in the Eikonal solver --- bsbash/build_functions.sh | 5 +- example_configs/rabbit_mesh_eikonal.ini | 19 ++++++++ src/config/config_parser.h | 1 + src/eikonal/eikonal_solver.c | 45 +++++++----------- src/main_eikonal.c | 62 ++++++++++++++++++------- 5 files changed, 81 insertions(+), 51 deletions(-) create mode 100644 example_configs/rabbit_mesh_eikonal.ini diff --git a/bsbash/build_functions.sh b/bsbash/build_functions.sh index 42644a16..8b0d390a 100755 --- a/bsbash/build_functions.sh +++ b/bsbash/build_functions.sh @@ -21,7 +21,7 @@ LIBRARY_OUTPUT_DIRECTORY=$ROOT_DIR GLOBAL_FORCE_COMPILATION="" QUIET='' WAIT_ENTER='' -WRITE_COMPILE_COMMANDS='y' +WRITE_COMPILE_COMMANDS='' DEFAULT_BUILD_DIR="build_" COMPILE_COMMANDS_FILE="${ROOT_DIR}/compile_commands.json" @@ -397,9 +397,6 @@ COMPILE_OBJECT () { ANY_COMPILED="y" fi - - - } CHECK_HEADERS() { diff --git a/example_configs/rabbit_mesh_eikonal.ini b/example_configs/rabbit_mesh_eikonal.ini new file mode 100644 index 00000000..6df592a2 --- /dev/null +++ b/example_configs/rabbit_mesh_eikonal.ini @@ -0,0 +1,19 @@ +[main] +simulation_time = 100 +dt = 1 + +[domain] +;These values are mandatory +name=UCLA Rabbit Mesh +;this mesh always start at 250.0 +maximum_discretization = 250.0 +main_function=initialize_grid_with_rabbit_mesh +;These can be optional depending on the domain main_function +mesh_file=meshes/rabheart.alg + +[save_result] +;/////mandatory///////// +output_dir=./outputs/rabbit_eikonal +main_function=save_as_text_or_binary +;////////////////// +file_prefix=T diff --git a/src/config/config_parser.h b/src/config/config_parser.h index e3ea72c4..13f5d267 100755 --- a/src/config/config_parser.h +++ b/src/config/config_parser.h @@ -246,6 +246,7 @@ void free_batch_options(struct batch_options * options); void free_visualization_options(struct visualization_options * options); void free_conversion_options(struct conversion_options *options); void free_fibers_conversion_options(struct fibers_conversion_options *options); +void free_eikonal_options(struct eikonal_options *options); void set_or_overwrite_common_data(struct config* config, const char *key, const char *value, const char *section, const char *config_file); diff --git a/src/eikonal/eikonal_solver.c b/src/eikonal/eikonal_solver.c index e4eed189..63045ec6 100644 --- a/src/eikonal/eikonal_solver.c +++ b/src/eikonal/eikonal_solver.c @@ -1,3 +1,7 @@ +// +// Created by sachetto on 06/06/24. +// Based on https://github.com/SCIInstitute/StructuredEikonal + #include #include #include "eikonal_solver.h" @@ -28,7 +32,6 @@ struct eikonal_solver * new_eikonal_solver(bool verbose) { void free_eikonal_solver(struct eikonal_solver *solver) { - //TODO: free 3D arrays FREE_3D_ARRAY(solver->speeds, solver->width, solver->height); FREE_3D_ARRAY(solver->answer, solver->width, solver->height); FREE_3D_ARRAY(solver->mask, solver->width, solver->height); @@ -251,19 +254,26 @@ static void get_solution(struct eikonal_solver *solver) { for(int k = 0; k < BLOCK_LENGTH; k++) { for(int j = 0; j < BLOCK_LENGTH; j++) { for(int i = 0; i < BLOCK_LENGTH; i++) { - double d = solver->memory_struct.h_sol[baseAddr + - k * BLOCK_LENGTH * BLOCK_LENGTH + - j * BLOCK_LENGTH + i]; + double d = solver->memory_struct.h_sol[baseAddr + k * BLOCK_LENGTH * BLOCK_LENGTH + j * BLOCK_LENGTH + i]; if ((i + bx * BLOCK_LENGTH) < solver->width && (j + by * BLOCK_LENGTH) < solver->height && (k + bz * BLOCK_LENGTH) < solver->depth) { - solver->answer[(i + bx * BLOCK_LENGTH)][(j + - by * BLOCK_LENGTH)][k + bz * BLOCK_LENGTH] = d; + solver->answer[(i + bx * BLOCK_LENGTH)][(j + by * BLOCK_LENGTH)][k + bz * BLOCK_LENGTH] = d; } } } } } + + for(int i = 0 ; i < solver->num_active_cells; i++) { + solver->active_cells[i]->v = solver->answer[(int)solver->active_cells[i]->center.x][(int)solver->active_cells[i]->center.y][(int)solver->active_cells[i]->center.z]; + + //Translating back to original space + solver->active_cells[i]->center.x = solver->active_cells[i]->center.x * solver->active_cells[i]->discretization.x + solver->active_cells[i]->discretization.x/2.0; + solver->active_cells[i]->center.y = solver->active_cells[i]->center.y * solver->active_cells[i]->discretization.y + solver->active_cells[i]->discretization.y/2.0; + solver->active_cells[i]->center.z = solver->active_cells[i]->center.z * solver->active_cells[i]->discretization.z + solver->active_cells[i]->discretization.z/2.0; + + } } void map_generator(struct eikonal_solver *solver) { @@ -389,29 +399,6 @@ void use_seeds(struct eikonal_solver *solver) { check_cuda_error( cudaMemset(solver->memory_struct.d_con, 1, volSize*sizeof(bool)) ); } -void write_alg(struct eikonal_solver *solver, char *filename) { - - FILE *out = fopen(filename, "w"); - - if (out == NULL) { - log_error_and_exit("Error opening file: %s\n", filename); - } - - for (int k = 0; k < solver->depth; k++) { - for (int j = 0; j < solver->height; j++) { - for (int i = 0; i < solver->width; i++) { - if(solver->mask[i][j][k] == true) { - real d = solver->answer[i][j][k]; - fprintf(out, "%lf, %lf, %lf, 0.5, 0.5, 0.5, %lf\n", i + 0.5, j + 0.5, k + 0.5, d); - } - } - } - } - - fclose(out); - -} - void solve_eikonal(struct eikonal_solver *solver) { if (solver->speeds == NULL) { diff --git a/src/main_eikonal.c b/src/main_eikonal.c index 76dc7f7f..61d7418c 100644 --- a/src/main_eikonal.c +++ b/src/main_eikonal.c @@ -1,13 +1,13 @@ -#include +#include "3dparty/ini_parser/ini.h" +#include "3dparty/stb_ds.h" #include "alg/cell/cell.h" #include "alg/grid/grid.h" -#include "3dparty/ini_parser/ini.h" +#include "config/config_parser.h" #include "config/domain_config.h" #include "eikonal/eikonal_solver.h" -#include "3dparty/stb_ds.h" #include "logger/logger.h" -#include "config/config_parser.h" - +#include "utils/file_utils.h" +#include static int compare_coordinates(const void *a, const void *b) { @@ -15,15 +15,15 @@ static int compare_coordinates(const void *a, const void *b) { struct cell_node *coord2 = *(struct cell_node **) b; if (coord1->center.y != coord2->center.y) { - return coord1->center.y - coord2->center.y; + return (int) (coord1->center.y - coord2->center.y); } if (coord1->center.x != coord2->center.x) { - return coord1->center.x - coord2->center.x; + return (int) (coord1->center.x - coord2->center.x); } - return coord1->center.z - coord2->center.z; + return (int) (coord1->center.z - coord2->center.z); } @@ -59,10 +59,15 @@ int main(int argc, char *argv[]) { } + //Translating to 0 for(int i = 0 ; i < grid->num_active_cells; i++) { - int new_pos_x = (grid->active_cells[i]->center.x-grid->active_cells[i]->discretization.x/2)/grid->active_cells[i]->discretization.x; - int new_pos_y = (grid->active_cells[i]->center.y-grid->active_cells[i]->discretization.y/2)/grid->active_cells[i]->discretization.y; - int new_pos_z = (grid->active_cells[i]->center.z-grid->active_cells[i]->discretization.z/2)/grid->active_cells[i]->discretization.z; + double new_pos_x = (grid->active_cells[i]->center.x-grid->active_cells[i]->discretization.x/2.0)/grid->active_cells[i]->discretization.x; + double new_pos_y = (grid->active_cells[i]->center.y-grid->active_cells[i]->discretization.y/2.0)/grid->active_cells[i]->discretization.y; + double new_pos_z = (grid->active_cells[i]->center.z-grid->active_cells[i]->discretization.z/2.0)/grid->active_cells[i]->discretization.z; + + if(new_pos_x != floor(new_pos_x) || new_pos_y != floor(new_pos_y) || new_pos_z != floor(new_pos_z)) { + log_error_and_exit("The current version only accepts integer coordinates\n"); + } grid->active_cells[i]->center.x = new_pos_x; grid->active_cells[i]->center.y = new_pos_y; @@ -72,25 +77,46 @@ int main(int argc, char *argv[]) { size_t itersPerBlock = 10, type = 1; struct eikonal_solver *solver = new_eikonal_solver(false); - solver->width = grid->cube_side_length.x/grid->active_cells[0]->discretization.x; - solver->height = grid->cube_side_length.y/grid->active_cells[0]->discretization.y; - solver->depth = grid->cube_side_length.z/grid->active_cells[0]->discretization.z; + solver->width = (int) (grid->cube_side_length.x/grid->active_cells[0]->discretization.x); + solver->height = (int) (grid->cube_side_length.y/grid->active_cells[0]->discretization.y); + solver->depth = (int) (grid->cube_side_length.z/grid->active_cells[0]->discretization.z); solver->solver_type = type; solver->iters_per_block = itersPerBlock; size_t *initial_seed = calloc(sizeof(size_t), 3); - initial_seed[0] = grid->active_cells[0]->center.x; - initial_seed[1] = grid->active_cells[0]->center.y; - initial_seed[2] = grid->active_cells[0]->center.z; + initial_seed[0] = 128; + initial_seed[1] = 48; + initial_seed[2] = 63; arrput(solver->seeds, initial_seed); solver->active_cells = grid->active_cells; solver->num_active_cells = grid->num_active_cells; solve_eikonal(solver); - write_alg(solver, "test.alg"); + struct config * save_result_config = eikonal_options->save_mesh_config; + + if(save_result_config) { + + init_config_functions(save_result_config, "./shared_libs/libdefault_save_mesh.so", "save_result"); + + print_save_mesh_config_values(save_result_config); + log_msg(LOG_LINE_SEPARATOR); + + char *out_dir_name = NULL; + GET_PARAMETER_STRING_VALUE_OR_USE_DEFAULT(out_dir_name, save_result_config, "output_dir"); + if(out_dir_name != NULL) { + create_dir(out_dir_name); + struct time_info ti = ZERO_TIME_INFO; + ((save_mesh_fn *)save_result_config->main_function)(&ti, save_result_config, grid, NULL, NULL); + } else { + log_warn("Not output dir provided. The result will not be saved!"); + } + } + + free(initial_seed); free_eikonal_solver(solver); + free_eikonal_options(eikonal_options); }