linesimplify.py 文件源码

def reduce(seq, maxpoints=20, straighten=False):
    if min(seq) == max(seq):
        return np.array(seq)
    x = map(float, range(0, len(seq)))
    # plt.plot(seq)
    VWpts = simplify(x, seq, maxpoints)
    xn, yn = map(list, zip(*VWpts))
    # plt.plot(xn, yn, linewidth=2, color='r')
    simple = np.interp(x, xn, yn)
    seq2 = np.copy(seq)
    points = []
    for xi, yi in enumerate(seq2):
        point = ogr.Geometry(ogr.wkbPoint)
        point.AddPoint(xi, yi)
        points.append(point)
    points = np.array(points)
    while True:
        poly = createPoly(xn, yn, len(seq2), max(seq2))
        line = ogr.Geometry(ogr.wkbLineString)
        for xi, yi in zip(xn, yn):
            line.AddPoint(xi, yi)
        dists = np.array([line.Distance(point) for point in points])
        contain = np.array([point.Within(poly) for point in points])
        dists[~contain] = 0
        points = points[(dists > 0)]
        dists = dists[(dists > 0)]
        if len(dists) == 0:
            break
        candidate = points[np.argmax(dists)]
        cp = candidate.GetPoint()
        index = np.argmin(np.array(xn) < cp[0])
        xn.insert(index, cp[0])
        yn.insert(index, cp[1])
    if straighten:
        indices = [i for i in range(0, len(xn)) if (xn[i], yn[i]) in VWpts]
        for i, j in enumerate(indices):
            if i < (len(indices)-1):
                if indices[i+1] > j+1:
                    dx = abs(xn[j] - xn[indices[i + 1]])
                    dy = abs(yn[j] - yn[indices[i + 1]])
                    if dx > dy:
                        seg_y = yn[j:indices[i + 1]+1]
                        for k in range(j, indices[i + 1]+1):
                            yn[k] = min(seg_y)
    # plt.plot(xn, yn, linewidth=2, color='g')
    # plt.show()
    return np.interp(x, xn, yn).astype(int)