More progress in the Eikonal solver

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() {

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

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);
 

src/eikonal/eikonal_solver.c

@@ -1,3 +1,7 @@
+//
+// Created by sachetto on 06/06/24.
+// Based on https://github.com/SCIInstitute/StructuredEikonal
+
 #include <assert.h>
 #include <cuda_runtime.h>
 #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) {

src/main_eikonal.c

@@ -1,29 +1,29 @@
-#include <stdlib.h>
+#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 <stdlib.h>
 
 static int compare_coordinates(const void *a, const void *b) {
 
     struct cell_node *coord1 = *(struct  cell_node **) a;
     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);
 
 }