Find neighbored segments/labels in Matrix

Hi all!
I have applied a segmentation on an array and want to find the neighbors of each segment.
For example:

Input:
4×6 Matrix{Int64}:
1  1  1  2  2  2
1  1  5  5  4  4
7  7  3  3  3  3
9  9  9  9  8  8

Output:
Dict{Int64, Matrix{Int64}} with 8 entries:
 5 => [1 2 … 7 3]
 4 => [2 5 3]
 7 => [1 5 3 9]
 2 => [1 5 4]
 9 => [7 3 8]
 8 => [3 9]
 3 => [5 4 … 8 9]
 1 => [2 5 7 3]

I’ve tried region_adjacent_graph function of the ImageSegmentation.jl package which does what I need but it takes 20 minutes for a 10000x10000 Matrix with ~15000 segments.
Any ideas how to get the neighbors faster?
Bjoern

Is each segment always contained in a row? I know nothing about this field, but my guess is that it would be easy to hand-craft a faster algorithm for your particular use-case.

No, segments could also be over several rows.

And are the values in the matrix always numbers between 1 and n when there are n segments?

yes

Here’s a simple implementation, let me know how it works :slight_smile:

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 = CartesianIndex(-1, -1):CartesianIndex(1, 1)
    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
2 Likes

Thanks a lot @sijo. Your function works pretty well and fast, just 2 seconds on the big array.

1 Like

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 :slight_smile: )

2 Likes