使用Cython处理变量大小的列表

Question

我想把Sutherland-Hogman algorithm的python implementation符号化。该算法根据非常简单的规则（边的内部或外部等）更新顶点列表，但细节并不重要。这是python版本，它接受顺时针方向的多边形顶点列表。例如：

sP=[(50, 150),  (200, 50),  (350, 150), (350, 300), (250, 300), (200, 250), (150, 350),(100, 250), (100, 200)]
cP=[(100, 100), (300, 100), (300, 300), (100, 300)]

然后计算它们的交集：

^{pr2}$

这是在rosettacode上找到的代码，稍微修改一下，如果没有交集，则返回一个空列表。在

def clip(subjectPolygon, clipPolygon):
   def inside(p):
      return(cp2[0]-cp1[0])*(p[1]-cp1[1]) > (cp2[1]-cp1[1])*(p[0]-cp1[0])

   def computeIntersection():
      dc = [ cp1[0] - cp2[0], cp1[1] - cp2[1] ]
      dp = [ s[0] - e[0], s[1] - e[1] ]
      n1 = cp1[0] * cp2[1] - cp1[1] * cp2[0]
      n2 = s[0] * e[1] - s[1] * e[0] 
      n3 = 1.0 / (dc[0] * dp[1] - dc[1] * dp[0])
      return [(n1*dp[0] - n2*dc[0]) * n3, (n1*dp[1] - n2*dc[1]) * n3]

   outputList = subjectPolygon
   cp1 = clipPolygon[-1]

   for clipVertex in clipPolygon:
      cp2 = clipVertex
      inputList = outputList
      outputList = []
      s = inputList[-1]

      for subjectVertex in inputList:
         e = subjectVertex
         if inside(e):
            if not inside(s):
               outputList.append(computeIntersection())
            outputList.append(e)
         elif inside(s):
            outputList.append(computeIntersection())
         s = e
      if len(outputList)<1:
          return []
      cp1 = cp2
   return(outputList)

这个函数对于我的应用程序来说非常慢，所以我试图用numpy来实现它的cythozing。这是我的赛顿版本。我必须定义两个输入端的错误信息。在

天鹅1

cimport cython
import numpy as np
cimport numpy as np

def clip(np.ndarray[np.float32_t, ndim=2] subjectPolygon,np.ndarray[np.float32_t, ndim=2] clipPolygon):


    outputList = list(subjectPolygon)
    cdef np.ndarray[np.float32_t, ndim=1] cp1 = clipPolygon[-1,:]
    cdef np.ndarray[np.float32_t, ndim=1] cp2 

    for i in xrange(clipPolygon.shape[0]):
       cp2 = clipPolygon[i]
       inputList = outputList
       outputList = []
       s = inputList[-1]

       for subjectVertex in inputList:
          e = subjectVertex
          if inside(e, cp1, cp2):
             if not inside(s, cp1, cp2):
                outputList.append(computeIntersection(cp1, cp2, e, s))
             outputList.append(e)
          elif inside(s, cp1, cp2):
             outputList.append(computeIntersection(cp1, cp2, e, s))
          s = e
       if len(outputList)<1:
         return []
       cp1 = cp2

    return(outputList)


def computeIntersection(np.ndarray[np.float32_t, ndim=1] cp1, np.ndarray[np.float32_t, ndim=1] cp2, np.ndarray[np.float32_t, ndim=1] e, np.ndarray[np.float32_t, ndim=1] s):
    cdef np.ndarray[np.float32_t, ndim=1] dc = cp1-cp2
    cdef np.ndarray[np.float32_t, ndim=1] dp = s-e
    cdef np.float32_t n1 = cp1[0] * cp2[1] - cp1[1] * cp2[0]
    cdef np.float32_t n2 = s[0] * e[1] - s[1] * e[0] 
    cdef np.float32_t n3 = 1.0 / (dc[0] * dp[1] - dc[1] * dp[0])
    cdef np.ndarray[np.float32_t, ndim=1] res=np.array([(n1*dp[0] - n2*dc[0]) * n3, (n1*dp[1] - n2*dc[1]) * n3], dtype=np.float32)
    return res

def inside(np.ndarray[np.float32_t, ndim=1] p, np.ndarray[np.float32_t, ndim=1] cp1, np.ndarray[np.float32_t, ndim=1] cp2):
    cdef bint b=(cp2[0]-cp1[0])*(p[1]-cp1[1]) > (cp2[1]-cp1[1])*(p[0]-cp1[0])
    return b 

当我计算两个版本的时间时，我只得到了2倍的加速，我需要至少10倍（或100倍！）。有什么事要做吗？ 如何与赛顿处理列表？在

编辑1:我听从了@DavidW的建议，我分配了numpy数组并对其进行了修剪，而不是使用list，我现在使用的是cdef函数，它应该可以带来10倍的速度，但不幸的是，我没有看到任何加速！在

半胱氨酸

cimport cython
import numpy as np
cimport numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
def clip(np.ndarray[np.float32_t, ndim=2] subjectPolygon,np.ndarray[np.float32_t, ndim=2] clipPolygon):
    return clip_in_c(subjectPolygon, clipPolygon)

@cython.boundscheck(False)
@cython.wraparound(False)
cdef np.ndarray[np.float32_t, ndim=2] clip_in_c(np.ndarray[np.float32_t, ndim=2] subjectPolygon,np.ndarray[np.float32_t, ndim=2] clipPolygon):





    cdef int cp_size=clipPolygon.shape[0]
    cdef int outputList_effective_size=subjectPolygon.shape[0]
    cdef int inputList_effective_size=outputList_effective_size
    #We allocate a fixed size array of size 
    cdef int max_size_inter=outputList_effective_size*cp_size
    cdef int k=-1

    cdef np.ndarray[np.float32_t, ndim=2] outputList=np.empty((max_size_inter,2), dtype=np.float32)
    cdef np.ndarray[np.float32_t, ndim=2] inputList=np.empty((max_size_inter,2), dtype=np.float32)
    cdef np.ndarray[np.float32_t, ndim=1] cp1 = clipPolygon[cp_size-1,:]
    cdef np.ndarray[np.float32_t, ndim=1] cp2=np.empty((2,), dtype=np.float32)

    outputList[:outputList_effective_size]=subjectPolygon
    for i in xrange(cp_size):

        cp2 = clipPolygon[i]
        inputList[:outputList_effective_size] = outputList[:outputList_effective_size]
        inputList_effective_size=outputList_effective_size
        outputList_effective_size=0
        s = inputList[inputList_effective_size-1]

        for j in xrange(inputList_effective_size):
            e = inputList[j]
            if inside(e, cp1, cp2):
                if not inside(s, cp1, cp2):                    
                    k+=1
                    outputList[k]=computeIntersection(cp1, cp2, e, s)

                k+=1
                outputList[k]=e


            elif inside(s, cp1, cp2):
                k+=1
                outputList[k]=computeIntersection(cp1, cp2, e, s)

            s = e

        if k<0:
            return np.empty((0,0),dtype=np.float32)





        outputList_effective_size=k+1

        cp1 = cp2
        k=-1

    return outputList[:outputList_effective_size]


cdef np.ndarray[np.float32_t, ndim=1] computeIntersection(np.ndarray[np.float32_t, ndim=1] cp1, np.ndarray[np.float32_t, ndim=1] cp2, np.ndarray[np.float32_t, ndim=1] e, np.ndarray[np.float32_t, ndim=1] s):
    cdef np.ndarray[np.float32_t, ndim=1] dc = cp1-cp2
    cdef np.ndarray[np.float32_t, ndim=1] dp = s-e
    cdef np.float32_t n1 = cp1[0] * cp2[1] - cp1[1] * cp2[0]
    cdef np.float32_t n2 = s[0] * e[1] - s[1] * e[0] 
    cdef np.float32_t n3 = 1.0 / (dc[0] * dp[1] - dc[1] * dp[0])
    return np.array([(n1*dp[0] - n2*dc[0]) * n3, (n1*dp[1] - n2*dc[1]) * n3], dtype=np.float32)


cdef bint inside(np.ndarray[np.float32_t, ndim=1] p, np.ndarray[np.float32_t, ndim=1] cp1, np.ndarray[np.float32_t, ndim=1] cp2):
    return (cp2[0]-cp1[0])*(p[1]-cp1[1]) > (cp2[1]-cp1[1])*(p[0]-cp1[0])

以下是基准：

import numpy as np
from cython1 import clip_cython1
from cython2 import clip_cython2
import time


sp=np.array([[50, 150],[200,50],[350,150],[250,300],[200,250],[150,350],[100,250],[100,200]],dtype=np.float32)
cp=np.array([[100,100],[300,100],[300,300],[100,300]],dtype=np.float32)

t1=time.time()
for i in xrange(120000):
    a=clip_cython1(sp, cp)
t2=time.time()
print (t2-t1)

t1=time.time()
for i in xrange(120000):
    a=clip_cython2(sp, cp)
t2=time.time()
print (t2-t1)

39.45秒

44.12款

第二个更糟！在

编辑2CodeReview的@Peter Taylor给出的最佳答案是，每次在U s内部计算时，都是多余的，因为s=e，而您已经在_e内部计算（并将dc和n1分解到函数之外，但这没有多大帮助）。在

cimport cython
import numpy as np
cimport numpy as np


def clip(np.ndarray[np.float32_t, ndim=2] subjectPolygon,np.ndarray[np.float32_t, ndim=2] clipPolygon):


    outputList = list(subjectPolygon)
    cdef np.ndarray[np.float32_t, ndim=1] cp1 = clipPolygon[-1,:]
    cdef np.ndarray[np.float32_t, ndim=1] cp2 
    cdef bint inside_e, inside_s
    cdef np.float32_t n1
    cdef np.ndarray[np.float32_t, ndim=1] dc 
    cdef int i



    for i in range(clipPolygon.shape[0]):

       cp2 = clipPolygon[i]
       #intermediate
       n1 = cp1[0] * cp2[1] - cp1[1] * cp2[0]
       dc=cp1-cp2
       inputList = outputList
       outputList = []
       s = inputList[-1]

       inside_s=inside(s, cp1, dc)
       for index, subjectVertex in enumerate(inputList):

           e = subjectVertex
           inside_e=inside(e, cp1, dc)
           if inside_e:

               if not inside_s:
                   outputList.append(computeIntersection(dc, n1, e, s))
               outputList.append(e)

           elif inside_s:
               outputList.append(computeIntersection(dc, n1, e, s))

           s = e
           inside_s=inside_e

       if len(outputList)<1:
           return []
       cp1 = cp2

    return(outputList)


cdef np.ndarray[np.float32_t, ndim=1] computeIntersection(np.ndarray[np.float32_t, ndim=1] dc, np.float32_t n1, np.ndarray[np.float32_t, ndim=1] e, np.ndarray[np.float32_t, ndim=1] s):
    cdef np.ndarray[np.float32_t, ndim=1] dp = s-e
    cdef np.float32_t n2 = s[0] * e[1] - s[1] * e[0] 
    cdef np.float32_t n3 = 1.0 / (dc[0] * dp[1] - dc[1] * dp[0])
    return np.array([(n1*dp[0] - n2*dc[0]) * n3, (n1*dp[1] - n2*dc[1]) * n3], dtype=np.float32)


cdef bint inside(np.ndarray[np.float32_t, ndim=1] p, np.ndarray[np.float32_t, ndim=1] cp1, np.ndarray[np.float32_t, ndim=1] dc):
    return (-dc[0])*(p[1]-cp1[1]) > (-dc[1])*(p[0]-cp1[0])

混合了两个版本（只有numpy数组和@Peter Taylor的技巧效果稍差）。不知道为什么？可能是因为我们要分配一个很长的列表特殊形状[0]*cp.形状[0]？在

Answer 1

我的速度提高了15倍：

In [12]: timeit clippy.clip(clippy.sP, clippy.cP)
10000 loops, best of 3: 126 µs per loop
In [13]: timeit clippy.clip1(clippy.sP, clippy.cP)
10000 loops, best of 3: 75.9 µs per loop
In [14]: timeit myclip.clip(clippy.sP, clippy.cP)
10000 loops, best of 3: 47.1 µs per loop
In [15]: timeit myclip.clip1(clippy.sP, clippy.cP)
100000 loops, best of 3: 8.2 µs per loop

clippy.clip是您的原始函数。在

clippy.clip1也是Python，但是用元组解包代替了大多数列表索引。在

^{pr2}$

myclip.clip是原始的cythonized；仍在处理列表，而不是数组。在

myclip.clip1是第二个cythonized：

cdef computeIntersection1(double cp10, double cp11, double cp20, double cp21,double s0, double s1,double e0, double e1):
  dc0, dc1 = cp10 - cp20, cp11 - cp21
  dp0, dp1 =  s0 - e0, s1 - e1
  n1 = cp10 * cp21 - cp11 * cp20
  n2 = s0 * e1 - s1 * e0
  n3 = 1.0 / (dc0 * dp1 - dc1 * dp0)
  return (n1*dp0 - n2*dc0) * n3, (n1*dp1 - n2*dc1) * n3

cdef cclip1(subjectPolygon, clipPolygon):
   cdef double cp10, cp11, cp20, cp21
   cdef double s0, s1, e0, e1
   cdef double s_in, e_in

   outputList = subjectPolygon
   cp10, cp11 = clipPolygon[-1]

   for cp20, cp21 in clipPolygon:

      inputList = outputList
      #print(inputList)
      outputList = []
      s0,s1 = inputList[-1]
      #s_in = inside(s0, s1)
      s_in = (cp20-cp10)*(s1-cp11) - (cp21-cp11)*(s0-cp10)
      for e0, e1  in inputList:
         #e_in = inside(e0, e1)
         e_in = (cp20-cp10)*(e1-cp11) - (cp21-cp11)*(e0-cp10)
         if e_in>0:
            if s_in<=0:
               outputList.append(computeIntersection1(cp10,cp11,cp20,cp21,s0,s1,e0,e1))
            outputList.append((e0, e1))
         elif s_in>0:
            outputList.append(computeIntersection1(cp10,cp11,cp20,cp21,s0,s1,e0,e1))
         s0,s1,s_in = e0,e1,e_in
      if len(outputList)<1:
          return []
      cp10, cp11 = cp20, cp21
   return outputList

def clip1(subjectPolygon, clipPolygon):
    return cclip1(subjectPolygon, clipPolygon)

带-a注释的html仍然显示了相当多的黄色，但是大多数计算不需要Python。在compute函数中，有一个Python检查0除数，以及Python调用来构建返回元组。元组解包仍然调用Python。所以还有改进的余地。在

在Python代码中，使用numpy没有任何好处。列表很小，访问列表元素的速度更快。但是在cython数组可能是typed memoryviews和纯C代码的垫脚石。在

---

其他时间。在

第二次编辑：

In [24]: timeit edit2.clip(np.array(clippy.sP,np.float32), np.array(clippy.cP,np
    ...: .float32))
1000 loops, best of 3: 228 µs per loop

@Matt's边框

In [25]: timeit clippy.polygon_clip(clippy.rp,clippy.cp,100,100,300,300)
1000 loops, best of 3: 208 µs per loop

扩展类

我通过定义一个扩展类来清理代码

cdef class Point:
    cdef public double x, y
    def __init__(self, x, y):
        self.x = x
        self.y = y

这样我就可以写下：

  s = inputList[-1]
  s_in =  insideP(s, cp1, cp2)

“cover”函数必须将元组列表转换为点和v.v的列表

sP = [Point(*x) for x in subjectPolygon]

这有轻微的速度损失。在

Answer 2

在搞乱了你的Cython代码之后，我发现你的库已经在其他地方实现了会更容易些，所以请查看scikit映像版本，它只是几行Numpy代码，以及你从matplotlib中寻找的算法：

import numpy as np
from matplotlib import path, transforms

def polygon_clip(rp, cp, r0, c0, r1, c1):
"""Clip a polygon to the given bounding box.

Parameters
     
rp, cp : (N,) ndarray of double
    Row and column coordinates of the polygon.
(r0, c0), (r1, c1) : double
    Top-left and bottom-right coordinates of the bounding box.

Returns
   -
r_clipped, c_clipped : (M,) ndarray of double
    Coordinates of clipped polygon.

Notes
  -
This makes use of Sutherland-Hodgman clipping as implemented in
AGG 2.4 and exposed in Matplotlib.

"""
    poly = path.Path(np.vstack((rp, cp)).T, closed=True)
    clip_rect = transforms.Bbox([[r0, c0], [r1, c1]])
    poly_clipped = poly.clip_to_bbox(clip_rect).to_polygons()[0]

    # This should be fixed in matplotlib >1.5
    if np.all(poly_clipped[-1] == poly_clipped[-2]):
        poly_clipped = poly_clipped[:-1]

    return poly_clipped[:, 0], poly_clipped[:, 1]

如果没有其他任何东西，以上应该更容易转换成Cython。在

[更新] 从另一个Cython回答剖析中，尝试已经从C++实现多边形裁剪的这个包，称为https://pypi.python.org/pypi/pyclipper用法：

导入pyclipper

主题=( （（180，200），（260，200），（260，150），（180，150））， （（215160），（230，190），（200，190）） )在

夹子=（（190，210），（240，210），（240，130），（190，130））

pc=皮夹。皮夹() 电脑地址路径（夹子，pyclipper.PT_剪辑，正确） 个人电脑地址路径（主题：，pyclipper.PT_主题，正确）

解决方案=pc.执行(pyclipper.CT_交叉口, Pycliepper.PFT_EVENODD公司, Pycliepper.PFT_EVENODD公司)在

上面的速度和我糟糕的AMD工作电脑BTW9US上的快速Cython代码答案差不多。在