WIP: adding eikonal solver

build.sh

@@ -87,6 +87,7 @@ for i in "${BUILD_ARGS[@]}"; do
             COMPILE_POSTPROCESSOR='y'
             COMPILE_EXPAND='y'
             COMPILE_CLIP='y'
+            COMPILE_EIKONAL='y'
             ;;
         simulator)
             COMPILE_SIMULATOR='y'
@@ -294,6 +295,11 @@ if [ -n "$COMPILE_CLIP" ]; then
     COMPILE_EXECUTABLE "MonoAlg3D_clip_mesh" "src/main_clip_mesh.c" "" "$STATIC_DEPS" "$DYNAMIC_DEPS" "$EXECUTABLES_LIBRARY_PATH $EXTRA_LIB_PATH"
 fi
 
+if [ -n "$COMPILE_EIKONAL" ]; then
+    COMPILE_EXECUTABLE "MonoAlg3D_eikonal_solver" "src/main_eikonal.c" "" "$STATIC_DEPS" "$DYNAMIC_DEPS" "$EXECUTABLES_LIBRARY_PATH $EXTRA_LIB_PATH"
+fi
+
+
 FIND_CRITERION
 
 if [ -n "$CRITERION_FOUND" ]; then

src/gui/gui.c

@@ -1151,14 +1151,15 @@ void init_and_open_gui_window(struct gui_shared_info *gui_config) {
             Vector2 info_pos;
 
             if(gui_config->draw_type == DRAW_SIMULATION) {
-                text = TextFormat("Mouse is on Volume: %.2lf, %.2lf, %.2lf with grid position %i", gui_state->current_mouse_over_volume.position_draw.x,
+                text = TextFormat("Mouse is on Volume: %.2lf, %.2lf, %.2lf with grid position %i and value %.2lf", gui_state->current_mouse_over_volume.position_draw.x,
                         gui_state->current_mouse_over_volume.position_draw.y, gui_state->current_mouse_over_volume.position_draw.z,
-                        gui_state->current_mouse_over_volume.matrix_position + 1);
+                        gui_state->current_mouse_over_volume.matrix_position + gui_state->current_mouse_over_volume.v);
 
             } else {
 
-                text = TextFormat("Mouse is on Volume: %.2lf, %.2lf, %.2lf", gui_state->current_mouse_over_volume.position_draw.x,
-                        gui_state->current_mouse_over_volume.position_draw.y, gui_state->current_mouse_over_volume.position_draw.z);
+                text = TextFormat("Mouse is on Volume: %.2lf, %.2lf, %.2lf with value %.2lf", gui_state->current_mouse_over_volume.position_draw.x,
+                        gui_state->current_mouse_over_volume.position_draw.y, gui_state->current_mouse_over_volume.position_draw.z, 
+                        gui_state->current_mouse_over_volume.v);
             }
 
             text_size = MeasureTextEx(gui_state->font, text, gui_state->font_size_big, gui_state->font_spacing_big);

src/gui/gui_mesh_helpers.c

@@ -168,6 +168,7 @@ static void check_mouse_over_volume(struct voxel *voxel, struct gui_state *gui_s
             gui_state->current_mouse_over_volume.position_draw = (Vector3){p_mesh.x, p_mesh.y, p_mesh.z};
             gui_state->current_mouse_over_volume.matrix_position = voxel->matrix_position;
             gui_state->ray_mouse_over_hit_distance = collision_mouse_over.distance;
+            gui_state->current_mouse_over_volume.v = voxel->v;
         }
     }
 

src/main_eikonal.c

@@ -0,0 +1,225 @@
+#include <stdlib.h>
+#include <math.h>
+#include "3dparty/stb_ds.h"
+#include "alg/grid/grid.h"
+#include "config/domain_config.h"
+#include "config/save_mesh_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;
+}
+
+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;
+        if(heap->arr[i].grid_cell->center.x == cell->center.x && heap->arr[i].grid_cell->center.y == cell->center.y && 
+            heap->arr[i].grid_cell->center.z == cell->center.z) {
+            return i;
+        }
+    }
+    return -1;
+}
+
+static real euclidean_distance(struct cell_node *cell1, struct cell_node *cell2) {
+
+    real dx = cell1->center.x - cell2->center.x;
+    real dy = cell1->center.y - cell2->center.y;
+    real dz = cell1->center.z - cell2->center.z;
+
+    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 *neighbour_grid_cell = neighbour;
+    struct cell_node *grid_cell = x.grid_cell;
+    bool has_found;
+
+    struct transition_node *white_neighbor_cell;
+
+    uint16_t neighbour_grid_cell_level = ((struct basic_cell_data *)(neighbour_grid_cell))->level;
+    enum cell_type neighbour_grid_cell_type = ((struct basic_cell_data *)(neighbour_grid_cell))->type;
+
+    if(neighbour_grid_cell_level > grid_cell->cell_data.level) {
+        if(neighbour_grid_cell_type == TRANSITION_NODE) {
+            has_found = false;
+            while(!has_found) {
+                if(neighbour_grid_cell_type == TRANSITION_NODE) {
+                    white_neighbor_cell = (struct transition_node *)neighbour_grid_cell;
+                    if(white_neighbor_cell->single_connector == NULL) {
+                        has_found = true;
+                    } else {
+                        neighbour_grid_cell = white_neighbor_cell->quadruple_connector1;
+                        neighbour_grid_cell_type = ((struct basic_cell_data *)(neighbour_grid_cell))->type;
+                    }
+                } else {
+                    break;
+                }
+            }
+        }
+    } else {
+        if(neighbour_grid_cell_level <= grid_cell->cell_data.level && (neighbour_grid_cell_type == TRANSITION_NODE)) {
+            has_found = false;
+            while(!has_found) {
+                if(neighbour_grid_cell_type == TRANSITION_NODE) {
+                    white_neighbor_cell = (struct transition_node *)(neighbour_grid_cell);
+                    if(white_neighbor_cell->single_connector == 0) {
+                        has_found = true;
+                    } else {
+                        neighbour_grid_cell = white_neighbor_cell->single_connector;
+                        neighbour_grid_cell_type = ((struct basic_cell_data *)(neighbour_grid_cell))->type;
+                    }
+                } else {
+                    break;
+                }
+            }
+        }
+    }
+
+    if(neighbour_grid_cell_type == CELL_NODE) {
+
+        struct cell_node *neighbour_cell = (struct cell_node *)(neighbour_grid_cell);                
+
+        if(!neighbour_cell->active || in_array(known, neighbour_cell) != -1) {
+            return;
+        }
+
+        real tentative_distance = x.distance + wave_speed * euclidean_distance(grid_cell, neighbour_grid_cell);
+        real neighbour_distance = hmget(*distances, neighbour_cell->center);
+
+        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);
+
+
+        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);
+
+            if(index == -1) {
+                xi.distance = tentative_distance;                
+                heap_push(heap, xi);
+            } else {
+                heap->arr[index].distance = tentative_distance;
+            }
+        }
+    }
+}
+
+
+int main(int argc, char **argv) {
+
+    struct grid *grid = new_grid();   
+
+    char *discretization = "200.0";
+    char *side_length = "20000.0";
+    char *num_layers = "1";
+    real wave_speed = 1.0;       
+
+    struct config *domain_config = domain_config = alloc_and_init_config_data();
+
+    domain_config->main_function_name = strdup("initialize_grid_with_square_mesh");
+    shput_dup_value(domain_config->config_data, "name", "Test custom mesh");
+
+    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");
+
+    int success = ((set_spatial_domain_fn*)domain_config->main_function)(domain_config, grid);
+
+    if(!success) {
+        clean_and_free_grid(grid);
+        free_config_data(domain_config);
+        log_error_and_exit("Error loading the domain");
+    }
+
+    order_grid_cells(grid);
+
+    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 point_hash_entry *times = NULL;
+    hmdefault(times, INFINITY);
+
+    struct heap_point *known = NULL;
+
+    //TODO: this has to come from the stimulus
+    struct heap_point first_node;
+
+    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);
+        }
+
+    }
+
+    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 j = 0; j <= LEFT; j++) {
+            compute_t(x, x.grid_cell->neighbours[j], known, wave_speed, trial, &distances, &times); 
+        }        
+    }
+
+    while(trial->size > 0) {
+        struct heap_point x = heap_pop(trial);
+        arrpush(known, x);
+        for (int i = 0; i <= LEFT; i++) {
+            compute_t(x, x.grid_cell->neighbours[i], known, wave_speed, trial, &distances, &times); 
+        }
+    }
+
+    // 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;
+
+    ((save_mesh_fn *)save_mesh_config->main_function)(&ti, save_mesh_config, grid, NULL, NULL);
+
+    free_config_data(save_mesh_config);
+
+}

src/utils/build.sh

@@ -1,4 +1,4 @@
-UTILS_SOURCE_FILES="search.c stop_watch.c sort.c file_utils.c batch_utils.c"
-UTILS_HEADER_FILES="utils.h stop_watch.h file_utils.h batch_utils.c"
+UTILS_SOURCE_FILES="search.c stop_watch.c sort.c file_utils.c batch_utils.c heap.c"
+UTILS_HEADER_FILES="utils.h stop_watch.h file_utils.h batch_utils.c heap.h"
 
 COMPILE_STATIC_LIB "utils" "$UTILS_SOURCE_FILES" "$UTILS_HEADER_FILES"

src/utils/heap.c

@@ -0,0 +1,90 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "heap.h"
+
+static void heapify(struct point_distance_heap* h, int index) {
+    int left = index * 2 + 1;
+    int right = index * 2 + 2;
+    int min = index;
+
+    if (left >= h->size || left < 0)
+        left = -1;
+    if (right >= h->size || right < 0)
+        right = -1;
+
+    if (left != -1 && h->arr[left].distance < h->arr[index].distance)
+        min = left;
+    if (right != -1 && h->arr[right].distance < h->arr[min].distance)
+        min = right;
+
+    if (min != index) {
+        struct heap_point temp = h->arr[min];
+        h->arr[min] = h->arr[index];
+        h->arr[index] = temp;
+        heapify(h, min);
+    }
+}
+
+struct point_distance_heap * build_heap(struct heap_point* data, int size, int capacity) {
+
+    struct point_distance_heap* h = (struct point_distance_heap *) malloc(sizeof(struct point_distance_heap));
+    if (h == NULL) {
+        fprintf(stderr, "Error allocating memory for heap\n");
+        return NULL;
+    }
+
+    h->capacity = capacity; 
+    h->size = size;
+
+    h->arr = data;
+
+    int i = (h->size - 2) / 2;
+    while (i >= 0) {
+        heapify(h, i);
+        i--;
+    }
+    return h;
+}
+
+void insert_helper(struct point_distance_heap* h, int index) {
+
+    int parent = (index - 1) / 2;
+
+    if (h->arr[parent].distance > h->arr[index].distance) {
+        struct heap_point temp = h->arr[parent];
+        h->arr[parent] = h->arr[index];
+        h->arr[index] = temp;
+        insert_helper(h, parent);
+    }
+}
+
+struct heap_point heap_pop(struct point_distance_heap* h) {
+
+    struct heap_point delete_item = {NULL, -1};
+
+    if (h->size == 0) {
+        printf("\nHeap id empty.");
+        return delete_item;
+    }
+
+    delete_item = h->arr[0];
+
+    h->arr[0] = h->arr[h->size - 1];
+    h->size--;
+
+    heapify(h, 0);
+    return delete_item;
+}
+
+void heap_push(struct point_distance_heap* h, struct heap_point data) {
+
+    // Checking if heap is full or not
+    if (h->size < h->capacity) {
+        // Inserting data into an array
+        h->arr[h->size] = data;
+        // Calling insertHelper function
+        insert_helper(h, h->size);
+        // Incrementing size of array
+        h->size++;
+    }
+}

src/utils/heap.h

@@ -0,0 +1,17 @@
+
+#include "../common_types/common_types.h"
+
+struct heap_point {
+    struct cell_node *grid_cell;
+    real distance;
+};
+
+struct point_distance_heap {
+    struct heap_point * arr;
+    int size;
+    int capacity;
+};
+
+struct point_distance_heap * build_heap(struct heap_point* data, int size, int capacity);
+struct heap_point heap_pop(struct point_distance_heap* h); 
+void heap_push(struct point_distance_heap *h, struct heap_point data);