By the way, with just one line changed the function will work for arrays of any dimension:
julia> function find_neighbors(A; nsegments=maximum(A))
# These two should be compared to see which one performs better:
flags = zeros(Bool, nsegments, nsegments)
#flags = falses(nsegments, nsegments) # Uses a BitArray instead of an array of Bool
indices = CartesianIndices(A)
neighbors = -oneunit(indices[1]):oneunit(indices[1]) # ← this line changed
for I in indices
for N in neighbors
if I+N in indices
flags[A[I+N], A[I]] = true
end
end
end
# Unset diagonal elements since a segment is not neighbor of itself
flags[diagind(flags)] .= false
return Dict(seg => findall(col) for (seg, col) in enumerate(eachcol(flags)))
end
find_neighbors (generic function with 1 method)
julia> A = [1 2 3
1 1 3
4 4 4;;;
1 2 3
1 3 3
1 5 5]
3×3×2 Array{Int64, 3}:
[:, :, 1] =
1 2 3
1 1 3
4 4 4
[:, :, 2] =
1 2 3
1 3 3
1 5 5
julia> find_neighbors(A)
Dict{Int64, Vector{Int64}} with 5 entries:
5 => [1, 3, 4]
4 => [1, 3, 5]
2 => [1, 3]
3 => [1, 2, 4, 5]
1 => [2, 3, 4, 5]
Isn’t this cool? (yeah I know it’s just a rehash of that one 
)