From f361fcf9112efdc6113b712d8cba0d763363fd4e Mon Sep 17 00:00:00 2001 From: Rafael Sachetto Date: Tue, 23 Apr 2024 16:26:31 -0300 Subject: [PATCH] WIP: isotropic eikonal solver --- src/config/config_common.h | 34 ++--- src/config/config_parser.c | 110 +++++++++++++++- src/config/config_parser.h | 18 ++- src/main_batch.c | 4 +- src/main_eikonal.c | 252 ++++++++++++++++++++++--------------- 5 files changed, 295 insertions(+), 123 deletions(-) diff --git a/src/config/config_common.h b/src/config/config_common.h index 0dc06088..bc23a8b1 100644 --- a/src/config/config_common.h +++ b/src/config/config_common.h @@ -46,24 +46,24 @@ void init_config_functions(struct config *config, char *default_lib, char *confi do { \ if((s) == NULL) { \ log_info(tag " no configuration.\n"); \ - return; \ + } else { \ + log_info(tag " configuration:\n"); \ + log_info(tag " library = %s\n", (s)->library_file_path); \ + log_info(tag " main function = %s\n", (s)->main_function_name); \ + \ + if((s)->init_function_name) { \ + log_info(tag " init function = %s\n", (s)->init_function_name); \ + } \ + \ + if((s)->end_function_name) { \ + log_info(tag " end function = %s\n", (s)->end_function_name); \ + } \ + \ + for(int __i = 0; __i < arrlen((s)->extra_function_names); __i++) { \ + log_info(tag " extra function %d = %s\n", __i+1, (s)->extra_function_names[__i]); \ + } \ + STRING_HASH_PRINT_KEY_VALUE_LOG(tag, (s)->config_data); \ } \ - log_info(tag " configuration:\n"); \ - log_info(tag " library = %s\n", (s)->library_file_path); \ - log_info(tag " main function = %s\n", (s)->main_function_name); \ - \ - if((s)->init_function_name) { \ - log_info(tag " init function = %s\n", (s)->init_function_name); \ - } \ - \ - if((s)->end_function_name) { \ - log_info(tag " end function = %s\n", (s)->end_function_name); \ - } \ - \ - for(int __i = 0; __i < arrlen((s)->extra_function_names); __i++) { \ - log_info(tag " extra function %d = %s\n", __i+1, (s)->extra_function_names[__i]); \ - } \ - STRING_HASH_PRINT_KEY_VALUE_LOG(tag, (s)->config_data); \ } while(0) #define CALL_EXTRA_FUNCTIONS(fn, ...) \ diff --git a/src/config/config_parser.c b/src/config/config_parser.c index 3013c04a..7d787610 100644 --- a/src/config/config_parser.c +++ b/src/config/config_parser.c @@ -12,6 +12,9 @@ static const char *batch_opt_string = "c:h?"; static const struct option long_batch_options[] = {{"config_file", required_argument, NULL, 'c'}}; +#define eikonal_opt_string batch_opt_string +#define long_eikonal_options long_batch_options + static const char *conversion_opt_string = "i:o:c:h?"; static const struct option long_conversion_options[] = {{"input", required_argument, NULL, 'i'}, {"output", required_argument, NULL, 'o'}, @@ -134,6 +137,15 @@ void display_batch_usage(char **argv) { exit(EXIT_FAILURE); } +void display_eikonal_usage(char **argv) { + + printf("Usage: %s [options]\n\n", argv[0]); + printf("Options:\n"); + printf("--config_file | -c [configuration_file_path]. Eikonal solver configuration file. Default NULL.\n"); + printf("--help | -h. Shows this help and exit \n"); + exit(EXIT_FAILURE); +} + void display_conversion_usage(char **argv) { printf("Usage: %s [options]\n\n", argv[0]); @@ -168,7 +180,7 @@ void display_visualization_usage(char **argv) { exit(EXIT_FAILURE); } -void maybe_issue_overwrite_warning(const char *var, const char *section, const char *old_value, const char *new_value, const char *config_file) { +static void maybe_issue_overwrite_warning(const char *var, const char *section, const char *old_value, const char *new_value, const char *config_file) { if(strcmp(old_value, new_value) != 0) { fprintf(stderr, "WARNING: option %s in %s was set in the file %s to %s and is being overwritten " @@ -177,6 +189,26 @@ void maybe_issue_overwrite_warning(const char *var, const char *section, const c } } +struct eikonal_options *new_eikonal_options() { + struct eikonal_options *user_args = (struct eikonal_options *)calloc(1, sizeof(struct eikonal_options)); + + user_args->stim_configs = NULL; + sh_new_arena(user_args->stim_configs); + shdefault(user_args->stim_configs, NULL); + + user_args->domain_config = NULL; + user_args->save_mesh_config = NULL; + + return user_args; +} + +void free_eikonal_options(struct eikonal_options *options) { + shfree(options->stim_configs); + free_config_data(options->domain_config); + free_config_data(options->save_mesh_config); + free(options); +} + struct batch_options *new_batch_options() { struct batch_options *user_args = (struct batch_options *)calloc(1, sizeof(struct batch_options)); sh_new_arena(user_args->config_to_change); @@ -1022,7 +1054,7 @@ void parse_batch_options(int argc, char **argv, struct batch_options *user_args) while(opt != -1) { switch(opt) { case 'c': - user_args->batch_config_file = strdup(optarg); + user_args->config_file = strdup(optarg); break; case 'h': /* fall-through is intentional */ case '?': @@ -1037,6 +1069,31 @@ void parse_batch_options(int argc, char **argv, struct batch_options *user_args) } } +void parse_eikonal_options(int argc, char **argv, struct eikonal_options *user_args) { + + int opt = 0; + int option_index; + + opt = getopt_long_only(argc, argv, eikonal_opt_string, long_eikonal_options, &option_index); + + while(opt != -1) { + switch(opt) { + case 'c': + user_args->config_file = strdup(optarg); + break; + case 'h': /* fall-through is intentional */ + case '?': + display_eikonal_usage(argv); + break; + default: + /* You won't actually get here. */ + break; + } + + opt = getopt_long(argc, argv, eikonal_opt_string, long_eikonal_options, &option_index); + } +} + void parse_conversion_options(int argc, char **argv, struct conversion_options *user_args) { int opt = 0; @@ -1799,6 +1856,53 @@ int parse_config_file(void *user, const char *section, const char *name, const c return 1; } +int parse_eikonal_config_file(void *options, const char *section, const char *name, const char *value) { + + struct eikonal_options *pconfig = (struct eikonal_options *)options; + + if(SECTION_STARTS_WITH(STIM_SECTION)) { + + struct config *tmp = (struct config *)shget(pconfig->stim_configs, section); + + if(tmp == NULL) { + tmp = alloc_and_init_config_data(); + shput(pconfig->stim_configs, section, tmp); + } + + if(MATCH_NAME("name")) { + fprintf(stderr, + "name is a reserved word and should not be used inside a stimulus config section. Found in %s. " + "Exiting!\n", + section); + exit(EXIT_FAILURE); + } else { + set_common_data(tmp, name, value); + } + + } else if(MATCH_SECTION(DOMAIN_SECTION)) { + + if(pconfig->domain_config == NULL) { + pconfig->domain_config = alloc_and_init_config_data(); + } + + set_common_data(pconfig->domain_config, name, value); + } else if(MATCH_SECTION(SAVE_RESULT_SECTION)) { + + if(pconfig->save_mesh_config == NULL) { + pconfig->save_mesh_config = alloc_and_init_config_data(); + } + + set_common_data(pconfig->save_mesh_config, name, value); + + } else { + + fprintf(stderr, "\033[33;5;7mInvalid name %s in section %s on the config file!\033[0m\n", name, section); + return 0; + } + + return 1; +} + int parse_preprocessor_config(void *user, const char *section, const char *name, const char *value) { static int function_counter = 0; @@ -1832,6 +1936,8 @@ int parse_preprocessor_config(void *user, const char *section, const char *name, return 1; } + + #define WRITE_INI_SECTION(SECTION) fprintf(ini_file, "[%s]\n", SECTION) #define WRITE_NAME_VALUE_WITHOUT_CHECK(NAME, VALUE, SPECIFIER) fprintf(ini_file, "%s = %" SPECIFIER "\n", NAME, VALUE) diff --git a/src/config/config_parser.h b/src/config/config_parser.h index 81a7336f..d19b8489 100755 --- a/src/config/config_parser.h +++ b/src/config/config_parser.h @@ -163,14 +163,20 @@ struct user_options { }; struct batch_options { - char *batch_config_file; + char *config_file; char *output_folder; char *initial_config; int num_simulations; bool skip_existing_run; struct string_hash_entry *config_to_change; }; - +struct eikonal_options { + char *config_file; + char *output_folder; + struct string_voidp_hash_entry *stim_configs; + struct config *domain_config; + struct config *save_mesh_config; +}; struct visualization_options { char *input; char *files_prefix; @@ -198,19 +204,23 @@ struct fibers_conversion_options { char *output_file; }; -void display_usage( char** argv ); +void display_usage(char **argv); void display_batch_usage(char **argv); +void display_eikonal_usage(char **argv); void display_conversion_usage(char **argv); void display_fibers_conversion_usage(char **argv); +void display_visualization_usage(char **argv); struct user_options * new_user_options(); struct batch_options * new_batch_options(); +struct eikonal_options * new_eikonal_options(); struct visualization_options * new_visualization_options(); struct conversion_options * new_conversion_options(); struct fibers_conversion_options * new_fibers_conversion_options(); void parse_options(int argc, char**argv, struct user_options *user_args); void parse_batch_options(int argc, char**argv, struct batch_options *user_args); +void parse_eikonal_options(int argc, char**argv, struct eikonal_options *user_args); void parse_visualization_options(int argc, char**argv, struct visualization_options *user_args); void parse_conversion_options(int argc, char **argv, struct conversion_options *user_args); void parse_fibers_conversion_options(int argc, char **argv, struct fibers_conversion_options *user_args); @@ -218,6 +228,7 @@ void get_config_file(int argc, char**argv, struct user_options *user_args); int parse_config_file(void* user, const char* section, const char* name, const char* value); int parse_batch_config_file(void *user, const char *section, const char *name, const char *value); +int parse_eikonal_config_file(void *user, const char *section, const char *name, const char *value); int parse_preprocessor_config(void* user, const char* section, const char* name, const char* value); int parse_converter_config_file(void *user, const char *section, const char *name, const char *value); @@ -231,7 +242,6 @@ 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 maybe_issue_overwrite_warning(const char *var, const char *section, const char *old_value, const char *new_value, const char *config_file); 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/main_batch.c b/src/main_batch.c index cd8af017..e8cadab5 100644 --- a/src/main_batch.c +++ b/src/main_batch.c @@ -35,8 +35,8 @@ int main(int argc, char **argv) { parse_batch_options(argc, argv, batch_options); - if (ini_parse(batch_options->batch_config_file, parse_batch_config_file, batch_options) < 0) { - fprintf(stderr, "Error: Can't load the config file %s\n", batch_options->batch_config_file); + if (ini_parse(batch_options->config_file, parse_batch_config_file, batch_options) < 0) { + fprintf(stderr, "Error: Can't load the config file %s\n", batch_options->config_file); MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE); } diff --git a/src/main_eikonal.c b/src/main_eikonal.c index 5ea04db4..3517aee4 100644 --- a/src/main_eikonal.c +++ b/src/main_eikonal.c @@ -1,25 +1,27 @@ #include #include -#include "3dparty/stb_ds.h" #include "alg/grid/grid.h" #include "config/domain_config.h" #include "config/save_mesh_config.h" +#include "config/stim_config.h" #include "utils/heap.h" #include "logger/logger.h" - - -static int in_array(struct heap_point *arr, struct cell_node *cell) { - - int n = arrlen(arr); - for(int i = 0; i < n; i++) { - if(arr[i].grid_cell->center.x == cell->center.x && arr[i].grid_cell->center.y == cell->center.y && - arr[i].grid_cell->center.z == cell->center.z) { - return i; - } - } - return -1; +#include "3dparty/ini_parser/ini.h" +#include "utils/file_utils.h" + +struct heap_point_hash_entry { + struct point_3d key; + struct heap_point value; +}; + +static struct point_3d compute_wave_speed(struct point_3d conductivity_tensor) { + //Compute wave speed based on conductivity tensor (e.g., using eigenvalues) + //Example: For isotropic conductivity, return sqrt(1 / conductivity_tensor) + //For anisotropic conductivity, compute wave speed based on eigenvalues + return POINT3D(sqrt(conductivity_tensor.x), sqrt(conductivity_tensor.y), sqrt(conductivity_tensor.z)); } +//TODO: add an hash map in the heap to make the search faster static int in_heap(struct point_distance_heap *heap, struct heap_point x, int n) { for(int i = 0; i < n; i++) { struct cell_node *cell = x.grid_cell; @@ -40,8 +42,7 @@ static real euclidean_distance(struct cell_node *cell1, struct cell_node *cell2) return sqrt(dx * dx + dy * dy + dz * dz); } -void compute_t(struct heap_point x, void *neighbour, struct heap_point *known, real wave_speed, struct point_distance_heap *heap, struct point_hash_entry **distances, struct point_hash_entry **time) { - +void compute_t(struct heap_point x, void *neighbour, struct heap_point_hash_entry *known, struct point_3d condutivity_tensor, struct point_distance_heap *heap) { void *neighbour_grid_cell = neighbour; struct cell_node *grid_cell = x.grid_cell; @@ -92,33 +93,37 @@ void compute_t(struct heap_point x, void *neighbour, struct heap_point *known, r struct cell_node *neighbour_cell = (struct cell_node *)(neighbour_grid_cell); - if(!neighbour_cell->active || in_array(known, neighbour_cell) != -1) { + struct heap_point tmp = hmget(known, neighbour_cell->center); + + if(!neighbour_cell->active || tmp.grid_cell != NULL) { return; } - real tentative_distance = x.distance + wave_speed * euclidean_distance(grid_cell, neighbour_grid_cell); - real neighbour_distance = hmget(*distances, neighbour_cell->center); + struct heap_point xi; + xi.grid_cell = neighbour_cell; - if(tentative_distance < neighbour_distance) { - struct heap_point xi; - xi.grid_cell = neighbour_cell; - - hmput(*distances, neighbour_cell->center, tentative_distance); - hmput(*time, neighbour_cell->center, tentative_distance/wave_speed); - - int index = in_heap(heap, xi, heap->size); + int index = in_heap(heap, xi, heap->size); + real neighbour_distance = INFINITY; - printf("Computing for %lf, %lf, %lf. Tentative %lf, distance %lf, index %d\n", - neighbour_cell->center.x, neighbour_cell->center.y, neighbour_cell->center.z, - tentative_distance, neighbour_distance, index); + real wave_speed = compute_wave_speed(condutivity_tensor).x; - if(index == -1) { - xi.distance = tentative_distance; - heap_push(heap, xi); - } else { - heap->arr[index].distance = tentative_distance; - } + //TODO: we have to know the neighbour's direction to calculate the wave speed + real tentative_distance = x.distance + wave_speed * euclidean_distance(grid_cell, neighbour_grid_cell); + real time = tentative_distance/wave_speed; + + if(index != -1) { + neighbour_distance = heap->arr[index].distance; + } else { + xi.distance = tentative_distance; + xi.grid_cell->v = time; + heap_push(heap, xi); + return; + } + + if(tentative_distance < neighbour_distance) { + heap->arr[index].distance = tentative_distance; + heap->arr[index].grid_cell->v = time; } } } @@ -126,100 +131,151 @@ void compute_t(struct heap_point x, void *neighbour, struct heap_point *known, r int main(int argc, char **argv) { - struct grid *grid = new_grid(); + struct eikonal_options *eikonal_options = new_eikonal_options(); + parse_eikonal_options(argc, argv, eikonal_options); + struct grid *grid = new_grid(); + + if (ini_parse(eikonal_options->config_file, parse_eikonal_config_file, eikonal_options) < 0) { + fprintf(stderr, "Error: Can't load the config file %s\n", eikonal_options->config_file); + exit(EXIT_FAILURE); + } + + struct config *domain_config = eikonal_options->domain_config; + struct config *save_mesh_config = eikonal_options->save_mesh_config; + struct string_voidp_hash_entry *stimuli_configs = eikonal_options->stim_configs; - char *discretization = "200.0"; - char *side_length = "20000.0"; - char *num_layers = "1"; - real wave_speed = 1.0; + bool save_to_file = (save_mesh_config != NULL); + char *out_dir_name = NULL; + + if(save_to_file) { + init_config_functions(save_mesh_config, "./shared_libs/libdefault_save_mesh.so", "save_result"); + GET_PARAMETER_STRING_VALUE_OR_USE_DEFAULT(out_dir_name, save_mesh_config, "output_dir"); - struct config *domain_config = domain_config = alloc_and_init_config_data(); + if(out_dir_name == NULL) { + log_error("No output directory provided to save the results! Exiting!\n"); + exit(EXIT_FAILURE); + } + + create_dir(out_dir_name); - domain_config->main_function_name = strdup("initialize_grid_with_square_mesh"); - shput_dup_value(domain_config->config_data, "name", "Test custom mesh"); + } else { + log_error("No configuration provided to save the results! Exiting!\n"); + free(out_dir_name); + exit(EXIT_FAILURE); + } - shput(domain_config->config_data, "start_dx", strdup(discretization)); - shput(domain_config->config_data, "start_dy", strdup(discretization)); - shput(domain_config->config_data, "start_dz", strdup(discretization)); - shput(domain_config->config_data, "side_length", strdup(side_length)); - shput(domain_config->config_data, "num_layers", strdup(num_layers)); - init_config_functions(domain_config, "./shared_libs/libdefault_domains.so", "domain"); + // Configure the functions and set the mesh domain + if(domain_config) { - int success = ((set_spatial_domain_fn*)domain_config->main_function)(domain_config, grid); + init_config_functions(domain_config, "./shared_libs/libdefault_domains.so", "domain"); - if(!success) { - clean_and_free_grid(grid); - free_config_data(domain_config); - log_error_and_exit("Error loading the domain"); - } + print_domain_config_values(domain_config); + log_msg(LOG_LINE_SEPARATOR); + + bool success = ((set_spatial_domain_fn *)eikonal_options->domain_config->main_function)(domain_config, grid); - order_grid_cells(grid); + if(!success) { + log_error_and_exit("Error configuring the tissue domain!\n"); + } + + order_grid_cells(grid); + } + struct point_3d condutivity_tensor = SAME_POINT3D(0.00005336); + struct heap_point *data = (struct heap_point *) malloc(grid->num_active_cells * sizeof(struct heap_point)); struct point_distance_heap *trial = build_heap(data, 0, grid->num_active_cells); - struct point_hash_entry *distances = NULL; - hmdefault(distances, INFINITY); + struct heap_point_hash_entry *known = NULL; + struct heap_point default_entry = {NULL, -1}; + hmdefault(known, default_entry); - struct point_hash_entry *times = NULL; - hmdefault(times, INFINITY); - struct heap_point *known = NULL; + if(stimuli_configs) { - //TODO: this has to come from the stimulus - struct heap_point first_node; + size_t n = shlen(stimuli_configs); + struct time_info time_info = {0.0, 0.0, 0.0, 0}; - for(int i = 0; i < grid->num_active_cells; i++) { - struct cell_node *cell = grid->active_cells[i]; - if(cell->center.x == 100.0 && cell->center.y == 100.0 && cell->center.z == 100.0) { - first_node.grid_cell = cell; - first_node.distance = 0.0; - arrpush(known, first_node); - hmput(distances, cell->center, 0.0); - hmput(times, cell->center, 0.0); + if(n > 0) { + STIM_CONFIG_HASH_FOR_INIT_FUNCTIONS(stimuli_configs); } - + + for(int i = 0; i < n; i++) { + + struct string_voidp_hash_entry e = stimuli_configs[i]; + log_info("Stimulus name: %s\n", e.key); + struct config *tmp = (struct config *)e.value; + print_stim_config_values(tmp); + log_msg(LOG_LINE_SEPARATOR); + } + + set_spatial_stim(&time_info, stimuli_configs, grid, false); + + struct config *tmp = NULL; + + uint32_t n_active = grid->num_active_cells; + struct cell_node **ac = grid->active_cells; + + for(size_t k = 0; k < n; k++) { + + real stim_start = 0.0; + real stim_dur = 0.0; + real stim_period = 0.0; + + tmp = (struct config *)stimuli_configs[k].value; + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real, stim_start, tmp, "start"); + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real, stim_dur, tmp, "duration"); + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, stim_period, tmp, "period"); + + for(uint32_t i = 0; i < n_active; i++) { + real value = ((real *)(tmp->persistent_data))[i]; + if(value != 0.0) { + struct heap_point node; + node.grid_cell = ac[i]; + if(stim_start == 0.0) { + node.distance = 0.0; + } else { + node.distance = compute_wave_speed(condutivity_tensor).x / stim_start; + } + hmput(known, ac[i]->center, node); + } + } + + + // if(stim_period > 0.0) { + // if(time >= stim_start + stim_period) { + // stim_start = stim_start + stim_period; + // sds stim_start_char = sdscatprintf(sdsempty(), "%lf", stim_start); + // shput_dup_value(tmp->config_data, "start", stim_start_char); + // sdsfree(stim_start_char); + // } + // } + + // time = cur_time; + } + } - - printf("First cell: %f, %f, %f\n", first_node.grid_cell->center.x, first_node.grid_cell->center.y, first_node.grid_cell->center.z); - for(int i = 0; i < arrlen(known); i++) { - struct heap_point x = known[i]; + for(int i = 0; i < hmlen(known); i++) { + struct heap_point x = known[i].value; for (int j = 0; j <= LEFT; j++) { - compute_t(x, x.grid_cell->neighbours[j], known, wave_speed, trial, &distances, ×); + compute_t(x, x.grid_cell->neighbours[j], known, condutivity_tensor, trial); } } while(trial->size > 0) { struct heap_point x = heap_pop(trial); - arrpush(known, x); + hmput(known, x.grid_cell->center, x); for (int i = 0; i <= LEFT; i++) { - compute_t(x, x.grid_cell->neighbours[i], known, wave_speed, trial, &distances, ×); + compute_t(x, x.grid_cell->neighbours[i], known, condutivity_tensor, trial); } } - // Print times hash map - for (int i = 0; i < arrlen(known); i++) { - struct cell_node *cell = known[i].grid_cell; - real time = hmget(times, cell->center); - cell->v = time; - printf("Time for cell (%f, %f, %f): %f\n", cell->center.x, cell->center.y, cell->center.z, time); - } - - struct config *save_mesh_config = alloc_and_init_config_data(); - - save_mesh_config->main_function_name = strdup("save_as_text_or_binary"); - shput_dup_value(save_mesh_config->config_data, "output_dir", "./test_eikonal"); - shput_dup_value(save_mesh_config->config_data, "print_rate", "1"); - init_config_functions(save_mesh_config, "shared_libs/libdefault_save_mesh.so", "save_result"); - - shput(save_mesh_config->config_data, "file_prefix", strdup("AT")); - struct time_info ti = ZERO_TIME_INFO; - + CALL_INIT_SAVE_MESH(save_mesh_config); ((save_mesh_fn *)save_mesh_config->main_function)(&ti, save_mesh_config, grid, NULL, NULL); - + CALL_END_SAVE_MESH(save_mesh_config, grid); free_config_data(save_mesh_config); }