使用Cython加速连接组件算法

Question

首先,我在windows xp机器上使用python [2.7.2],numpy [1.6.2rc1],cython [0.16],gcc [MinGW]编译器.

我需要一个3D连通分量算法来处理存储在numpy数组中的一些3D二进制数据(即1和0).不幸的是,我找不到任何现有的代码,所以我改编了这里的代码来处理3D数组.一切都很好,但是处理大量数据集的速度是可取的.结果我偶然发现了cython,并决定尝试一下.

到目前为止,cython已经提高了速度:Cython:0.339 s Python:0.635 s

使用cProfile,我在纯python版本中的耗时行是:

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

---

问题: "cythonize"线路的正确方法是什么:

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

任何帮助将不胜感激,希望这项工作将有助于其他人.

---

纯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

纯蟒蛇加速:

#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)

时间结果:
O:0.0220445
A:0.0002161
B:0.0173195
C:0.0002560

---

Cython版本[*.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加速:

使用:

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])

时间:0.170,原价:0.339秒

---

结果

在考虑了所提供的许多有用的评论和答案后,我目前的算法运行在:
Cython:0.0219
Python:0.4309

与纯蟒蛇相比,Cython的速度提高了20倍.

目前的Cython代码:

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

设置文件[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++",
                             )]
)

构建命令:

python setup.py build_ext --inplace

Answer 1

正如@gotgenes指出的那样,你绝对应该使用cython -a <file>,并试图减少你看到的黄色量.黄色对应的情况越来越严重.

我发现减少黄色量的事情:

1. 这看起来像是一种永远不会有任何越界数组访问的情况,只要输入Array有3个维度,所以可以关闭边界检查:

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

2. 给编译器以获取更多信息有效的索引,即只要你cdef的ndarray给予尽可能多的信息,您可以:

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

3. 为编译器提供有关正/负的更多信息.即如果你知道某个变量总是积极的,cdef那就unsigned int不是int,因为这意味着Cython可以消除任何负索引检查.

4. 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]
   

5. 避免使用该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]
   

6. 避免做花式索引/切片.特别是避免两次这样做!并避免使用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在此代码.)

7. 您不需要存储所有等价项,因为您对集合执行的唯一操作是查找最小元素.所以,如果你不是有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]
   

8. 毕竟,最后要做的就是不使用任何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.(请注意,它将需要一个C++编译器.)