# 3D Connected Component in Cython

This is my implementation of 3D connected component algorithm, which I use for a 255x512x512 binary matrix. Although it is written in Cython, it still takes quite some time both for 6-neighbors and 26-neighbors. Do you see any inefficiency in the algorithm? I would really appreciate if you could give me any kind of feedback regarding efficiency, readability and maintainability.

connectedComponent3D.pyx

import numpy as np
cimport numpy as np
cimport cython
from cpython cimport array
import array

@cython.wraparound(False)
@cython.nonecheck(False)
@cython.boundscheck(False)
def main(np.ndarray[np.int_t, ndim=3] arr):
"""
Input arguments: arr -> The array which connected component algo. is applied
type: 3D numpy array dtype=np.int
Returns: numberOfComponents -> Number of distinct components in the input array -> int
labels -> It has the same dims with input array. It contains the labels of each
element inside the input array(label 0 is background)

equivalencydict -> Contains which component label is connected to which component label -> dict
"""
cdef np.ndarray[np.int_t, ndim=3] labels = np.zeros_like(arr, dtype=np.int)
cdef dict equivalencydict = {}
cdef int numberOfComponents

labels, equivalencydict = getConnectedComponents(arr, equivalencydict)
labels = applyEquivalencyDictToLabels(labels, equivalencydict)
numberOfComponents = len(np.unique(labels)) - 1
return [numberOfComponents, labels]

@cython.wraparound(False)
@cython.nonecheck(False)
@cython.boundscheck(False)
cdef getConnectedComponents(np.ndarray[np.int_t, ndim=3] arr, dict equivalencydict):
"""
Input arguments: arr -> The array which connected component is applied -> 3D numpy array
dtype=np.int
equivalencydict -> Contains which component label is connected to which
component label type: dict
Returns: labels -> It has the same dims with input array. It contains the labels of each
element inside the input array(label 0 is background)
equivalencydict

For each element in the input array, if the element is True on the input array, the function
checks the labeled neighbors of that element. If there is no labeled neighbor, the element
gets a new label buy incrementing currentComp by 1. If there are labeled
neighbors, then the minimum of the neighbor labels is assigned to the
element.
"""
cdef np.ndarray[np.int_t, ndim=3] labels = np.zeros_like(arr, dtype=np.int)
cdef int currentComp = 0
cdef int z
cdef int y
cdef int x
cdef int zmax = arr.shape
cdef int ymax = arr.shape
cdef int xmax = arr.shape
cdef array.array neighborLabels = array.array('i', [])

#Iterate over each element of arr
for z in range(1, zmax - 1):
print(z)
for y in range(ymax):
for x in range(xmax):
if arr[z, y, x] == 0:
continue
else:
neighborLabels = getNeighborLabels(arr, labels, (z, y, x))
labels, equivalencydict, currentComp = labeling(labels, neighborLabels, currentComp,
equivalencydict, z, y, x)
return labels, equivalencydict

@cython.wraparound(False)
@cython.nonecheck(False)
@cython.boundscheck(False)
cdef getNeighborLabels(np.ndarray[np.int_t, ndim=3] arr, np.ndarray[np.int_t, ndim=3] labels,
tuple centerIndex):
"""
Input Arguments: centerIndex -> The index of the element which its neighbors are going to be searched
Returns: neighborLabels -> It contains the labels of suitable neighbors
(Neighbors which are true on the input array and are not labeled before during runtime are considered
as suitable)

neighbors: This is 2D numpy array which contains all the 26 neighbors around the element which is
located at centerIndex(z, y, x)
neighborsLabels: This is array which contains the labels of neighbors if there are any suitable
neighbors (Due to performance considerations, array is used instead of list)

If a neighbor is True on the input array and if it has been labeled before during runtime
(Otherwise it would be 0), then add the label of this neighbor to neighborLabels array.

"""
cdef np.ndarray[np.int_t, ndim=2] neighbors = get26Neighbors(centerIndex)
cdef array.array neighborLabels = array.array('i', [])
cdef np.ndarray[np.int_t, ndim=1] neighbor_index

for neighbor_index in neighbors:
if arr[tuple(neighbor_index)] == 1 and labels[tuple(neighbor_index)] != 0:
neighborLabels.append(labels[tuple(neighbor_index)])
return neighborLabels

@cython.wraparound(False)
@cython.nonecheck(False)
@cython.boundscheck(False)
cdef applyEquivalencyDictToLabels(np.ndarray[np.int_t, ndim=3] labels, dict equivalencydict):
"""
The equivalencydict is a sorted dictionary. By iterating over it in reverse direction,
the equivalencydict is applied to labels array.
"""
cdef tuple item

for item in sorted(list(equivalencydict.items()), key=lambda x:x, reverse=True):
print(item)
labels[labels == item] = item
return labels

@cython.wraparound(False)
@cython.nonecheck(False)
@cython.boundscheck(False)
cdef get6Neighbors(tuple index):
"""
Input arguments: index -> This is the index(z, y, x) of element whose neighbors are need to be
calculated. type: tuple
Returns: neighbors -> indices of 6-neighbors

This function calculates all 6 neighbors of an element in 3D space.
In order to see what a 6-neighbors is check the 29/38 slide in below link. Left figure is 6-n and
right one is 26-neighbors.
"""
cdef np.ndarray[np.int_t, ndim=2] neighbors = np.array([[index, index-1, index],
[index, index+1, index],
[index, index, index-1],
[index, index, index+1],
[index-1, index, index],
[index+1, index, index]], dtype=np.int)
return np.resize(neighbors, (6,3))

@cython.wraparound(False)
@cython.nonecheck(False)
@cython.boundscheck(False)
cdef labeling(np.ndarray[np.int_t, ndim=3] labels, array.array neighborLabels, int currentComp,
dict equivalencydict, int z, int y, int x):
"""
This function assigns the appropriate label to element(with indices z, y, x which are passed in as
input arguments).
If there is no suitable neighbor around the element(z, y, x) then a new component is created and
assigned to the element. Else, the minimum of the neigbors' labels around the element is
assigned to the element.
"""
if len(neighborLabels) == 0:
currentComp += 1
labels[z, y, x] = currentComp
else:
labels[z, y, x] = np.amin(neighborLabels)
equivalencydict = addLabelsToEquivalencyDict(labels, neighborLabels,
equivalencydict, z, y, x)
return labels, equivalencydict, currentComp

@cython.wraparound(False)
@cython.nonecheck(False)
@cython.boundscheck(False)
cdef addLabelsToEquivalencyDict(np.ndarray[np.int_t, ndim=3] labels, array.array neighborLabels,
dict equivalencydict, int z, int y, int x):
"""
This function creates the spatial relationship between the neighbors of an element(z, y, x).
The neighbors of an element are also connected to each other but they may not have the same
component due to the complexity of the objects in the input array. Equivalencydict dictionary
saves the information of which label is connected to which label. This equivalencydict is then
applied to labels at the end of the main function.
"""
cdef int label

for label in neighborLabels:
if label != labels[z, y, x]:
equivalencydict[label] = labels[z, y, x]
return equivalencydict

@cython.wraparound(False)
@cython.nonecheck(False)
@cython.boundscheck(False)
cdef get26Neighbors(tuple index):
"""
Input arguments: index -> This is the index(z, y, x) of element whose neighbors are need to be
calculated. type: tuple
Returns: neighbors -> indices of 26-neighbors

This function calculates all 16 neighbors of an element in 3D space.
In order to see what a 26-neighbors is check the 29/38 slide in below link. Left figure is 6-n and
right one is 26-neighbors.

"""
cdef np.ndarray zz
cdef np.ndarray yy
cdef np.ndarray xx

zz,yy,xx = np.mgrid[(index-1):(index+2) , (index-1):(index+2), (index-1):(index+2)]
cdef np.ndarray[np.int_t, ndim=2] neighbors = np.vstack((zz.flatten(), yy.flatten(), xx.flatten())).T.astype(np.int)
#Delete the center which is not a neighbor but the element itself and resize the neighbors
np.delete(neighbors, [36,37,38])
return np.resize(neighbors, (26,3))


connectedComponent3D_Setup.py file

from distutils.core import setup
from Cython.Build import cythonize
import numpy

setup(
ext_modules = cythonize("connectedComponent3D.pyx"),
include_dirs=[numpy.get_include()]
)


If you want to run the function first cd into a folder where the above two functions are at and run the following code on your command window which compiles the code:

python connectedComponent3D_Setup.py build_ext --inplace


Then you can import connectedComponent3D and call the connectedComponent3D.main(arr) function.

There's quite a bit of code here so I'm not going to be comprehensive, but here are some observations.

## Profiling

The first thing you should do is to try to profile your code and work out where time is actually being spent. Concentrate your optimization in these functions.

## get26Neighbors

I doubt if this function is limiting your speed. However, np.delete is often a little slow since it involves moving elements about in an already allocated array. I'd be tempted to allocate an empty array and fill it instead:

cdef np.ndarray[np.int_t, ndim=2] neighbors = np.empty((26,3),dtype=np.int)
neighbors[:12,0] = xx[:12]
neighbors[12:,0] = xx[13:]
neighbors[:12,1] = yy[:12]
# etc...
# (check my maths on which element to skip)


## getNeighborLabels

Indexing using arr[tuple(neighbor_index)] will be slow. In order to take advantage of Cython's ability to index arrays quickly you need to index using only integers and slices, not general Python objects like tuple. It's better to explicitly spell out the indices (event if this is a little more code):

for i in range(neighbors.shape):
xi = neighbors[i,0]
yi = neighbors[i,1]
zi = neighbors[i,2]
if arr[xi,yi,zi] == 1 and labels[xi,yi,zi] != 0:
neighborLabels.append(labels[xi,yi,zi])


I've also changed the outer loop to an integer indexed loop since this is much quicker in Cython than using the Python iterator protocol (for x in neighbours) - I know this is generally considered worse style in Python, but it is faster here. There's other places you could apply this.

I'd also be tempted to make neighborLabels a list rather than a numpy array, since Python lists are designed to have quick appends while numpy arrays aren't.

## General comments on typing

You've specified explicit types for (almost) everything. You don't generally get much benefit in Cython from specifying basic Python types such as dict or tuple. You also don't benefit from specifying the types of Numpy arrays unless you're indexing them directly. You may also lost a little speed due to unnecessary type checking. Therefore, for something like main (and other functions here...) where you don't index your arrays directly there is really very little point in specifying the types.

You also never get a lot out of np.ndarray (without specifying the type stored in the array).

## Compilation decorators

You applied

@cython.wraparound(False)
@cython.nonecheck(False)
@cython.boundscheck(False)


to every function. These only make a difference for Cython's fast indexing into Numpy arrays (with ints and slices) so for most functions here they will make no difference at all. It is better to think about where you actually need these (and where bounds checking might be helpful, for example) rather than applying these blindly in the hope they'll do something.

The is so useful, thank you for making this.

May I suggest one improvement: adding a small zero padding / zero stripping utility.

Because currently, if a connected component is at the edge, looking for it's out of bound neighbors will result an out of bound index error. Instead of checking for index, maybe simply padding the 3d volume with a 'Null' value will be a more versatile solution.

def zero_pad_3d(np_3d):
return np.pad(np_3d, ((1,1), (1,1), (1,1)), 'constant')

def zero_strip_3d(np_3d):
return np_3d[1:-1, 1:-1, 1:-1]


Then after calculations, strip away the paddings to return back the original size.

You can try adding annotate=True to your cythonize call, like so:

from distutils.core import setup
from Cython.Build import cythonize
import numpy

setup(
ext_modules = cythonize("connectedComponent3D.pyx",annotate=True),
include_dirs=[numpy.get_include()]
)


This will produce an HTML report next to the compilation results highlighting which lines of code still use the Python API as opposed to pure C with native types. This information, along with the profiling that DavidW suggested, will help you remove bottlenecks.