Geek asks and answers

Using Cython Algorithm to Accelerate Connected Components

Firstly, I use the python [2.7.2], numpy [1.6.2rc1], cython [0.16], gcc [MinGW] compiler on a Windows XP machine.

I needed an algorithm for three-dimensional connected components to process three-dimensional binary data (i.e. 1 s and 0 s) stored in numpy arrays. Unfortunately, I could not find any existing code, so I adapted the code found here to work with 3D arrays. Everything works fine, but speed is desirable for processing huge data sets. As a result, I came across cython and decided to give it a try.

So far, cython has improved speed: Cython: 0.339 s Python: 0.635 s

Using cProfile, my time length is in a pure python version:

new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax,yMin:yMax,zMin:zMax].ravel()))

Question: What is the correct way to "cythonize" strings:

new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax,yMin:yMax,zMin:zMax].ravel()))
for x,y,z in zip(ind[0],ind[1],ind[2]):

Any help would be appreciated and hopefully this work will help others.

Pure version of python [* .py]:

import numpy as np

def find_regions_3D(Array):
    x_dim=np.size(Array,0)
    y_dim=np.size(Array,1)
    z_dim=np.size(Array,2)
    regions = {}
    array_region = np.zeros((x_dim,y_dim,z_dim),)
    equivalences = {}
    n_regions = 0
    #first pass. find regions.
    ind=np.where(Array==1)
    for x,y,z in zip(ind[0],ind[1],ind[2]):

        # get the region number from all surrounding cells including diagnols (27) or create new region                        
        xMin=max(x-1,0)
        xMax=min(x+1,x_dim-1)
        yMin=max(y-1,0)
        yMax=min(y+1,y_dim-1)
        zMin=max(z-1,0)
        zMax=min(z+1,z_dim-1)

        max_region=array_region[xMin:xMax+1,yMin:yMax+1,zMin:zMax+1].max()

        if max_region > 0:
            #a neighbour already has a region, new region is the smallest > 0
            new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax+1,yMin:yMax+1,zMin:zMax+1].ravel()))
            #update equivalences
            if max_region > new_region:
                if max_region in equivalences:
                    equivalences[max_region].add(new_region)
                else:
                    equivalences[max_region] = set((new_region, ))
        else:
            n_regions += 1
            new_region = n_regions

        array_region[x,y,z] = new_region


    #Scan Array again, assigning all equivalent regions the same region value.
    for x,y,z in zip(ind[0],ind[1],ind[2]):
        r = array_region[x,y,z]
        while r in equivalences:
            r= min(equivalences[r])
        array_region[x,y,z]=r

    #return list(regions.itervalues())
    return array_region

Net python accelerations:

#Original line:
new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax+1,yMin:yMax+1,zMin:zMax+1].ravel()))

#ver A:
new_region = array_region[xMin:xMax+1,yMin:yMax+1,zMin:zMax+1]
min(new_region[new_region>0])

#ver B:
new_region = min( i for i in array_region[xMin:xMax,yMin:yMax,zMin:zMax].ravel() if i>0)

#ver C:
sub=array_region[xMin:xMax,yMin:yMax,zMin:zMax]
nlist=np.where(sub>0)
minList=[]
for x,y,z in zip(nlist[0],nlist[1],nlist[2]):
    minList.append(sub[x,y,z])
new_region=min(minList)

Time Results:
O: 0.0220445
: 0.0002161
B: 0.0173195
C: 0.0002560

Cython Version [* .pyx]:

import numpy as np
cimport numpy as np

DTYPE = np.int
ctypedef np.int_t DTYPE_t

cdef inline int int_max(int a, int b): return a if a >= b else b
cdef inline int int_min(int a, int b): return a if a <= b else b

def find_regions_3D(np.ndarray Array not None):
    cdef int x_dim=np.size(Array,0)
    cdef int y_dim=np.size(Array,1)
    cdef int z_dim=np.size(Array,2)
    regions = {}
    cdef np.ndarray array_region = np.zeros((x_dim,y_dim,z_dim),dtype=DTYPE)
    equivalences = {}
    cdef int n_regions = 0
    #first pass. find regions.
    ind=np.where(Array==1)
    cdef int xMin, xMax, yMin, yMax, zMin, zMax, max_region, new_region, x, y, z
    for x,y,z in zip(ind[0],ind[1],ind[2]):

        # get the region number from all surrounding cells including diagnols (27) or create new region                        
        xMin=int_max(x-1,0)
        xMax=int_min(x+1,x_dim-1)+1
        yMin=int_max(y-1,0)
        yMax=int_min(y+1,y_dim-1)+1
        zMin=int_max(z-1,0)
        zMax=int_min(z+1,z_dim-1)+1

        max_region=array_region[xMin:xMax,yMin:yMax,zMin:zMax].max()

        if max_region > 0:
            #a neighbour already has a region, new region is the smallest > 0
            new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax,yMin:yMax,zMin:zMax].ravel()))
            #update equivalences
            if max_region > new_region:
                if max_region in equivalences:
                    equivalences[max_region].add(new_region)
                else:
                    equivalences[max_region] = set((new_region, ))
        else:
            n_regions += 1
            new_region = n_regions

        array_region[x,y,z] = new_region


    #Scan Array again, assigning all equivalent regions the same region value.
    cdef int r
    for x,y,z in zip(ind[0],ind[1],ind[2]):
        r = array_region[x,y,z]
        while r in equivalences:
            r= min(equivalences[r])
        array_region[x,y,z]=r

    #return list(regions.itervalues())
    return array_region

Cython Accelerations:

Using:

cdef np.ndarray region = np.zeros((3,3,3),dtype=DTYPE)
...
        region=array_region[xMin:xMax,yMin:yMax,zMin:zMax]
        new_region=np.min(region[region>0])

Time: 0.170, original: 0.339 s

results

After reviewing many useful comments and answers, my current algorithms work on:
Cython: 0.0219
Python: 0.4309

Cython provides a 20-fold increase in speed compared to pure python.

Current Cython Code:

import numpy as np
import cython
cimport numpy as np
cimport cython

from libcpp.map cimport map

DTYPE = np.int
ctypedef np.int_t DTYPE_t

cdef inline int int_max(int a, int b): return a if a >= b else b
cdef inline int int_min(int a, int b): return a if a <= b else b

@cython.boundscheck(False)
def find_regions_3D(np.ndarray[DTYPE_t,ndim=3] Array not None):
    cdef unsigned int x_dim=np.size(Array,0),y_dim=np.size(Array,1),z_dim=np.size(Array,2)
    regions = {}
    cdef np.ndarray[DTYPE_t,ndim=3] array_region = np.zeros((x_dim,y_dim,z_dim),dtype=DTYPE)
    cdef np.ndarray region = np.zeros((3,3,3),dtype=DTYPE)
    cdef map[int,int] equivalences
    cdef unsigned int n_regions = 0

    #first pass. find regions.
    ind=np.where(Array==1)
    cdef np.ndarray[DTYPE_t,ndim=1] ind_x = ind[0], ind_y = ind[1], ind_z = ind[2]
    cells=range(len(ind_x))
    cdef unsigned int xMin, xMax, yMin, yMax, zMin, zMax, max_region, new_region, x, y, z, i, xi, yi, zi, val
    for i in cells:

        x=ind_x[i]
        y=ind_y[i]
        z=ind_z[i]

        # get the region number from all surrounding cells including diagnols (27) or create new region                        
        xMin=int_max(x-1,0)
        xMax=int_min(x+1,x_dim-1)+1
        yMin=int_max(y-1,0)
        yMax=int_min(y+1,y_dim-1)+1
        zMin=int_max(z-1,0)
        zMax=int_min(z+1,z_dim-1)+1

        max_region = 0
        new_region = 2000000000 # huge number
        for xi in range(xMin, xMax):
            for yi in range(yMin, yMax):
                for zi in range(zMin, zMax):
                    val = array_region[xi,yi,zi]
                    if val > max_region: # val is the new maximum
                        max_region = val

                    if 0 < val < new_region: # val is the new minimum
                        new_region = val

        if max_region > 0:
           if max_region > new_region:
                if equivalences.count(max_region) == 0 or new_region < equivalences[max_region]:
                    equivalences[max_region] = new_region
        else:
           n_regions += 1
           new_region = n_regions

        array_region[x,y,z] = new_region


    #Scan Array again, assigning all equivalent regions the same region value.
    cdef int r
    for i in cells:
        x=ind_x[i]
        y=ind_y[i]
        z=ind_z[i]

        r = array_region[x,y,z]
        while equivalences.count(r) > 0:
            r= equivalences[r]
        array_region[x,y,z]=r

    return array_region

Installation file [setup.py]

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy

setup(
    cmdclass = {'build_ext': build_ext},
    ext_modules = [Extension("ConnectComp", ["ConnectedComponents.pyx"],
                             include_dirs =[numpy.get_include()],
                             language="c++",
                             )]
)

Build Command:

python setup.py build_ext --inplace
+5
source share
3 answers

@gotgenes, cython -a <file> , . C.

, , :

  • , , Array 3 , :

    cimport cython

@cython.boundscheck(False)
def find_regions_3d(...):

  • , .. , cdef a ndarray , :

     def find_regions_3D(np.ndarray[DTYPE_t,ndim=3] Array not None):
     [...]
     cdef np.ndarray[DTYPE_t,ndim=3] array_region = ...
     [etc.]

  • / . , , cdef unsigned int, int, , Cython .

  • ind, ..

    ind = np.where(Array==1)
cdef np.ndarray[DTYPE_t,ndim=1] ind_x = ind[0], ind_y = ind[1], ind_z = ind[2]

  • for x,y,z in zip(..[0],..[1],..[2]).

    cdef int i
for i in range(len(ind_x)):
    x = ind_x[i]
    y = ind_y[i]
    z = ind_z[i]

  • /. ! filter!

    max_region=array_region[xMin:xMax,yMin:yMax,zMin:zMax].max()
if max_region > 0:
    new_region = min(filter(lambda i: i > 0, array_region[xMin:xMax,yMin:yMax,zMin:zMax].ravel()))
    if max_region > new_region:
        if max_region in equivalences:
            equivalences[max_region].add(new_region)
        else:
            equivalences[max_region] = set((new_region, ))


    max_region = 0
new_region = 2000000000 # "infinity"
for xi in range(xMin, xMax):
    for yi in range(yMin, yMax):
        for zi in range(zMin, zMax):
            val = array_region[xi,yi,zi]
            if val > max_region: # val is the new maximum
                max_region = val

            if 0 < val < new_region: # val is the new minimum
                new_region = val

if max_region > 0:
   if max_region > new_region:
       if max_region in equivalences:
           equivalences[max_region].add(new_region)
       else:
           equivalences[max_region] = set((new_region, ))
else:
   n_regions += 1
   new_region = n_regions

    , 10 C, - Python.

    (, cdef , xi, yi, zi val .)

  • , , , - . , equivalences int int,

    if max_region in equivalences:
    equivalences[max_region].add(new_region)
else:
    equivalences[max_region] = set((new_region, ))

[...]

while r in equivalences:
    r = min(equivalences[r])


    if max_region not in equivalences or new_region < equivalences[max_region]:
    equivalences[max_region] = new_region

[...]

while r in equivalences:
    r = equivalences[r]

  • , , - Python, , equivalences. , int int, from libcpp.map cimport map, cdef map[int,int] equivalences .. not in equivalences equivalences.count(..) == 0 .. in equivalences equivalences.count(..) > 0. ( , ++.)

+5

( , )

, scipy ndimage.label , ( , ). , :

from scipy import ndimage 
ndimage.label(your_data, connectivity_struct)

(, , ..)

+2

cython , C, Python, . - C , Py *. cdef python.

In your code, I would assume, for example, that the c loop zipcreates many python objects, and it would be much faster to iterate over the index int, which is then used to get the elements in ind[0], .... But look at the generated C code and see what seems to be calling unnecessarily many python functions.

0
source

More articles:


All Articles