constant kernel crash in Delaunay triangulation code (Cython)

EDIT:

After much testing the problem has been found to be in the following loop. I put the loop in its own function. The use of it is to accept a pointer to an array of Vertex structs and remove all repeated Vertices. This is done by checking their id and not comparing their coordinates. Because each Vertex is given a unique id at the very start, it makes it more straightforward like this.

The code to fix would then simplify to this

The basic structures:

ctypedef struct Vertex:
    double x
    double y
    int vid_nr

    


cdef Vertex create_vertex(double x_coord, double y_coord, int vid_number):
    cdef Vertex v
    v.x = x_coord
    v.y = y_coord
    v.vid_nr = vid_number
    return v


cdef Vertex* alloc_varr(int size, double var):
    cdef Vertex* varray = <Vertex*>malloc(size * sizeof(Vertex))
    cdef size_t i
    for i in range(size):
        varray[i] = create_vertex(var, var, i)
    return varray

The basic deallocation method. Tested and functions as expected

cdef Vertex* dealloc_from_varray(Vertex* original_array, int current_size, int remove_index):
    #careful remove_index is actual position - 1

    cdef Vertex* new_array = <Vertex*>malloc((current_size-1) * sizeof(Vertex))
    cdef size_t i
    cdef size_t j
    cdef int upper_limit = remove_index
    for i in range(upper_limit):
        new_array[i] = original_array[i]
    for j in range(upper_limit+1, current_size):
        new_array[j-1] = original_array[j]
        
    free(original_array)
    return new_array

And now the method that deallocates repeated Vertices by vid_nr

cdef void dealloc_non_uniques(Vertex* free_varray, int size):
    cdef int k = 0 # Initialize k
    cdef int l = 1 #initialized at first iteration but to be safe
    cdef int vsize = size
    #dealloc all non-unique verrtices from the free arr
    while k < vsize - 1:
        l = k + 1
        while l < vsize:
            if free_varray[l].vid_nr == free_varray[k].vid_nr:
                free_varray = dealloc_from_varray(free_varray, vsize, l)
                vsize = vsize - 1
            else:
                l = l + 1
        k = k + 1
    size = vsize

Below a simple test

cdef int varr_size = 7
cdef Vertex* og_array = alloc_varr(varr_size, 0)

og_array[0].x = 0.2
og_array[0].y = 2
og_array[1].x = -1.8
og_array[1].y = 0.15
og_array[2].x = -0.1
og_array[2].y = -2.3
og_array[3].x = 2
og_array[3].y = -0.3
og_array[4].x = -1.8
og_array[4].y = 6
og_array[5].x = -3.9
og_array[5].y = -2
og_array[6].x = 1.7
og_array[6].y = -3.4
og_array[6].vid_nr = og_array[5].vid_nr
og_array[2].vid_nr = og_array[4].vid_nr
og_array[3].vid_nr = og_array[4].vid_nr




print(varr_size)
for i in range(varr_size):
    print(og_array[i])


#print("\ntest dealloc_non_uniques:\nwe expect varr_size to be 4")

#dealloc_non_uniques(og_array, varr_size)

will crash the Kernel once i remove the comment # and run the dealloc_non_uniques method. Why ?

After a long time i finally found what the problems were. First of all, declared variables are immutable. So when i declare the size of the whole mesh as a global, to call it later for plotting reasons, i cannot update it on the fly with the triangulate function. This was fairly easy to get around. I can just allocate the appropriate size in memory like <int*>global_mesh_size = .... and then modify it inside the triangulate function by equating global_mesh_size[0] to the temporary mesh_size int i create inside the function itself.

The second error in the first version i posted here was in how the loops were written. Because the final mesh size is unknown, when i iterate through its elements with the intention to remove i have to make sure the loop limits are adjusted dynamically. This could not be done in a for range() loop because range() takes the initial input and creates the literal range that will remain static until the loop ends. Therefore it was necessary to update the range manually in the function and use a while loop with a simple break condition. This allowed me to deallocate bad Triangles from the final_mesh without reaching out of bounds elements.

The third error was a very simple syntax mistake in the sort_by_angle() function where i was accidentally using the .x instead of the .y coordinate to calculate the temporary vector magnitude before doing the norm and the projection on the x-axis. This lead to the Vertexes not being correctly sorted most of the time. This was the most time-consuming error to spot because there were many lines of code and i just ended up adding a lot of print() statements to see if the algebra was being done correctly.

Overall for the following test unit

cdef int varr_ssize = 4
cdef Vertex* og_array = alloc_varr(varr_ssize, 0)
cdef int* global_mesh_size = <int*>malloc(sizeof(int))
global_mesh_size[0] = 0 #dynamically changed through triangulate fx

og_array[0].x = 0.2
og_array[0].y = 2
og_array[1].x = -1.8
og_array[1].y = 0.15
og_array[2].x = -0.1
og_array[2].y = -2.3
og_array[3].x = 2
og_array[3].y = -0.3


cdef Face* testmesh = triangulate(og_array, varr_ssize, global_mesh_size)

print("\nglobal mesh size is:", global_mesh_size[0])
for f in range(global_mesh_size[0]):
    print(testmesh[f])

free(testmesh)
free(og_array)
free(global_mesh_size)

gives the correct triangulation of

global mesh size is: 2

{'ver1': {'x': 2.0, 'y': -0.3, 'vid_nr': 3}, 'ver2': {'x': 0.2, 'y': 2.0, 'vid_nr': 0}, 'ver3': {'x': -1.8, 'y': 0.15, 'vid_nr': 1}, 'fid_nr': 0}
{'ver1': {'x': 2.0, 'y': -0.3, 'vid_nr': 3}, 'ver2': {'x': -1.8, 'y': 0.15, 'vid_nr': 1}, 'ver3': {'x': -0.1, 'y': -2.3, 'vid_nr': 2}, 'fid_nr': 1}

To avoid clutter in this thread i will only link the Github page to the code because it is too long to read comfortably in the forum layout. As usual only the first few Cells in the notebook are relevant. Anything below the first output is just for testing purposes

Delaunay triangulation through Bowyer-Watson algorithm

Thanks to everyone who helped <3

EDIT: For a more complex example