4
\$\begingroup\$

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[0]
    cdef int ymax = arr.shape[1]
    cdef int xmax = arr.shape[2]
    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[0], reverse=True):
        print(item)
        labels[labels == item[0]] = item[1]
    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.
    Link: http://slideplayer.com/slide/8645709/
    """
    cdef np.ndarray[np.int_t, ndim=2] neighbors = np.array([[index[0], index[1]-1, index[2]],
                                                           [index[0], index[1]+1, index[2]],
                                                           [index[0], index[1], index[2]-1], 
                                                           [index[0], index[1], index[2]+1],
                                                           [index[0]-1, index[1], index[2]], 
                                                           [index[0]+1, index[1], index[2]]], 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.
    Link: http://slideplayer.com/slide/8645709/

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

    zz,yy,xx = np.mgrid[(index[0]-1):(index[0]+2) , (index[1]-1):(index[1]+2), (index[2]-1):(index[2]+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.

\$\endgroup\$
2
\$\begingroup\$

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[0]):
    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.

\$\endgroup\$
2
\$\begingroup\$

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.

\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.