Hi,

I’ve been trying to implement a version of the famous marching cubes algorithm and have arrived at the following abomination:

```
function marching_cubes(values::Array{T,3},points::Array{Point3D{T},3},iso,cube_size=1) where T
vertices = Array{Point3f0,1}()
indices = Array{Tuple{Int,Int,Int},1}()
vertex_list = Array{Point3f0,1}(12)
ind_pushed = false
cube_indices = [(0,0,0),(cube_size,0,0),(cube_size,cube_size,0),(0,cube_size,0),(0,0,cube_size),(cube_size,0,cube_size),(cube_size,cube_size,cube_size),(0,cube_size,cube_size)]
edge_table,tri_table = get_edge_tri_table()
i_max,j_max,k_max = size(points)
for k=1:cube_size:k_max-1
for j=1:cube_size:j_max-1
for i=1:cube_size:i_max-1
cube_index = zero(Int32)
for n = 1:8
if values[i+cube_indices[n][1],j+cube_indices[n][2],k+cube_indices[n][3]]<iso
cube_index |= Int32(2^(n-1))
end
end
cube_index += Int32(1)
edge_entry = edge_table[cube_index]
if edge_entry == 0
continue
end
if edge_entry & 1 == 1
vertex_list[1] = vertex_interp(iso,points[i+cube_indices[1][1],j+cube_indices[1][2],k+cube_indices[1][3]],points[i+cube_indices[2][1],j+cube_indices[2][2],k+cube_indices[2][3]],values[i+cube_indices[1][1],j+cube_indices[1][2],k+cube_indices[1][3]],values[i+cube_indices[2][1],j+cube_indices[2][2],k+cube_indices[2][3]])
end
if edge_entry & 2 == 2
vertex_list[2] = vertex_interp(iso,points[i+cube_indices[2][1],j+cube_indices[2][2],k+cube_indices[2][3]],points[i+cube_indices[3][1],j+cube_indices[3][2],k+cube_indices[3][3]],values[i+cube_indices[2][1],j+cube_indices[2][2],k+cube_indices[2][3]],values[i+cube_indices[3][1],j+cube_indices[3][2],k+cube_indices[3][3]])
end
if edge_entry & 4 == 4
vertex_list[3] = vertex_interp(iso,points[i+cube_indices[3][1],j+cube_indices[3][2],k+cube_indices[3][3]],points[i+cube_indices[4][1],j+cube_indices[4][2],k+cube_indices[4][3]],values[i+cube_indices[3][1],j+cube_indices[3][2],k+cube_indices[3][3]],values[i+cube_indices[4][1],j+cube_indices[4][2],k+cube_indices[4][3]])
end
if edge_entry & 8 == 8
vertex_list[4] = vertex_interp(iso,points[i+cube_indices[4][1],j+cube_indices[4][2],k+cube_indices[4][3]],points[i+cube_indices[1][1],j+cube_indices[1][2],k+cube_indices[1][3]],values[i+cube_indices[4][1],j+cube_indices[4][2],k+cube_indices[4][3]],values[i+cube_indices[1][1],j+cube_indices[1][2],k+cube_indices[1][3]])
end
if edge_entry & 16 == 16
vertex_list[5] = vertex_interp(iso,points[i+cube_indices[5][1],j+cube_indices[5][2],k+cube_indices[5][3]],points[i+cube_indices[6][1],j+cube_indices[6][2],k+cube_indices[6][3]],values[i+cube_indices[5][1],j+cube_indices[5][2],k+cube_indices[5][3]],values[i+cube_indices[6][1],j+cube_indices[6][2],k+cube_indices[6][3]])
end
if edge_entry & 32 == 32
vertex_list[6] = vertex_interp(iso,points[i+cube_indices[6][1],j+cube_indices[6][2],k+cube_indices[6][3]],points[i+cube_indices[7][1],j+cube_indices[7][2],k+cube_indices[7][3]],values[i+cube_indices[6][1],j+cube_indices[6][2],k+cube_indices[6][3]],values[i+cube_indices[7][1],j+cube_indices[7][2],k+cube_indices[7][3]])
end
if edge_entry & 64 == 64
vertex_list[7] = vertex_interp(iso,points[i+cube_indices[7][1],j+cube_indices[7][2],k+cube_indices[7][3]],points[i+cube_indices[8][1],j+cube_indices[8][2],k+cube_indices[8][3]],values[i+cube_indices[7][1],j+cube_indices[7][2],k+cube_indices[7][3]],values[i+cube_indices[8][1],j+cube_indices[8][2],k+cube_indices[8][3]])
end
if edge_entry & 128 == 128
vertex_list[8] = vertex_interp(iso,points[i+cube_indices[8][1],j+cube_indices[8][2],k+cube_indices[8][3]],points[i+cube_indices[5][1],j+cube_indices[5][2],k+cube_indices[5][3]],values[i+cube_indices[8][1],j+cube_indices[8][2],k+cube_indices[8][3]],values[i+cube_indices[5][1],j+cube_indices[5][2],k+cube_indices[5][3]])
end
if edge_entry & 256 == 256
vertex_list[9] = vertex_interp(iso,points[i+cube_indices[1][1],j+cube_indices[1][2],k+cube_indices[1][3]],points[i+cube_indices[5][1],j+cube_indices[5][2],k+cube_indices[5][3]],values[i+cube_indices[1][1],j+cube_indices[1][2],k+cube_indices[1][3]],values[i+cube_indices[5][1],j+cube_indices[5][2],k+cube_indices[5][3]])
end
if edge_entry & 512 == 512
vertex_list[10] = vertex_interp(iso,points[i+cube_indices[2][1],j+cube_indices[2][2],k+cube_indices[2][3]],points[i+cube_indices[6][1],j+cube_indices[6][2],k+cube_indices[6][3]],values[i+cube_indices[2][1],j+cube_indices[2][2],k+cube_indices[2][3]],values[i+cube_indices[6][1],j+cube_indices[6][2],k+cube_indices[6][3]])
end
if edge_entry & 1024 == 1024
vertex_list[11] = vertex_interp(iso,points[i+cube_indices[3][1],j+cube_indices[3][2],k+cube_indices[3][3]],points[i+cube_indices[7][1],j+cube_indices[7][2],k+cube_indices[7][3]],values[i+cube_indices[3][1],j+cube_indices[3][2],k+cube_indices[3][3]],values[i+cube_indices[7][1],j+cube_indices[7][2],k+cube_indices[7][3]])
end
if edge_entry & 2048 == 2048
vertex_list[12] = vertex_interp(iso,points[i+cube_indices[4][1],j+cube_indices[4][2],k+cube_indices[4][3]],points[i+cube_indices[8][1],j+cube_indices[8][2],k+cube_indices[8][3]],values[i+cube_indices[4][1],j+cube_indices[4][2],k+cube_indices[4][3]],values[i+cube_indices[8][1],j+cube_indices[8][2],k+cube_indices[8][3]])
end
table_entry = tri_table[cube_index]
n = 1
while table_entry[n]!=-1
vert = vertex_list[table_entry[n]+1]
push!(vertices,vert)
if !ind_pushed
push!(indices,(i,j,k))
ind_pushed = true
end
n+=1
end
ind_pushed = false
end
end
end
return vertices,indices
end
```

This results in no real unwanted allocations through each loop, having @time output of

`0.026461 seconds (720 allocations: 88.844 KiB)`

However as you can see this code looks awful and I was wondering how to do this more cleanly.

I’ve tried a lot of things which didn’t work. Defining arrays t_val, t_p through selection of the correct elements resulted in allocations of the order of 300mb, defining array views resulted in allocations of 100mb and slow execution, defining the correct values of indices in each loop resulted in 200mb allocations.

I guess what I’m trying to understand is how to not allocate any extra arrays or values inside the loop, but be able to define what part of the array to use inside the loop.

If there is a totally different way of doing this that leads to more readable code and similar performance please let me know i’m very keen to learn how to utilize all the tools in julia!

cheers

EDIT: Upon request, a working code via pastebin:

https://pastebin.com/0sKpb643

Only thing that changed is Point3D -> Point3f0 (Point3f0 is basically static array of 3 Float32’s) . Both have same functionality but one is custom the other isn’t.

Minimal requirement: GeometryTypes

Visualization requirement: GLVisualize,ColorTypes