8

我有一堆带有两组数据的交叉图,并且一直在寻找一种 matploltib 方法来用平滑的多边形轮廓突出显示它们的绘制区域。

目前我只是使用 Adob​​e Illustrator 并修改保存的绘图,但这并不理想。例子:在此处输入图像描述

对于任何指向示例的指针/链接,我将不胜感激。

干杯

4

2 回答 2

18

在这里,你有一个例子。我写了主要思想,但显然,你可以做得更好。

一个简短的解释:

1)您需要计算凸包(http://en.wikipedia.org/wiki/Convex_hull

2) 使用船体,您可以对其进行缩放以将所有数据保存在里面。

3)您必须对结果曲线进行插值。

第一部分是在http://wiki.scipy.org/Cookbook/Finding_Convex_Hull中完成的。第二个是微不足道的。第三个很笼统,你可以执行任何方法,有很多不同的方法可以做到这一点。我采用了@Jaime 的方法(任意轮廓的平滑样条表示,f(length) --> x,y),我认为这是一个非常好的方法。

我希望它可以帮助你...

#Taken from http://wiki.scipy.org/Cookbook/Finding_Convex_Hull

import numpy as n, pylab as p, time

def _angle_to_point(point, centre):
    '''calculate angle in 2-D between points and x axis'''
    delta = point - centre
    res = n.arctan(delta[1] / delta[0])
    if delta[0] < 0:
        res += n.pi
    return res

def _draw_triangle(p1, p2, p3, **kwargs):
    tmp = n.vstack((p1,p2,p3))
    x,y = [x[0] for x in zip(tmp.transpose())]
    p.fill(x,y, **kwargs)

def area_of_triangle(p1, p2, p3):
    '''calculate area of any triangle given co-ordinates of the corners'''
    return n.linalg.norm(n.cross((p2 - p1), (p3 - p1)))/2.


def convex_hull(points, graphic=False, smidgen=0.0075):
    '''
    Calculate subset of points that make a convex hull around points
    Recursively eliminates points that lie inside two neighbouring points until only convex hull is remaining.

    :Parameters:
    points : ndarray (2 x m)
    array of points for which to find hull
    graphic : bool
    use pylab to show progress?
    smidgen : float
    offset for graphic number labels - useful values depend on your data range

    :Returns:
    hull_points : ndarray (2 x n)
    convex hull surrounding points
    '''

    if graphic:
        p.clf()
        p.plot(points[0], points[1], 'ro')
    n_pts = points.shape[1]
    assert(n_pts > 5)
    centre = points.mean(1)
    if graphic: p.plot((centre[0],),(centre[1],),'bo')
    angles = n.apply_along_axis(_angle_to_point, 0, points, centre)
    pts_ord = points[:,angles.argsort()]
    if graphic:
        for i in xrange(n_pts):
            p.text(pts_ord[0,i] + smidgen, pts_ord[1,i] + smidgen, \
                   '%d' % i)
    pts = [x[0] for x in zip(pts_ord.transpose())]
    prev_pts = len(pts) + 1
    k = 0
    while prev_pts > n_pts:
        prev_pts = n_pts
        n_pts = len(pts)
        if graphic: p.gca().patches = []
        i = -2
        while i < (n_pts - 2):
            Aij = area_of_triangle(centre, pts[i],     pts[(i + 1) % n_pts])
            Ajk = area_of_triangle(centre, pts[(i + 1) % n_pts], \
                                   pts[(i + 2) % n_pts])
            Aik = area_of_triangle(centre, pts[i],     pts[(i + 2) % n_pts])
            if graphic:
                _draw_triangle(centre, pts[i], pts[(i + 1) % n_pts], \
                               facecolor='blue', alpha = 0.2)
                _draw_triangle(centre, pts[(i + 1) % n_pts], \
                               pts[(i + 2) % n_pts], \
                               facecolor='green', alpha = 0.2)
                _draw_triangle(centre, pts[i], pts[(i + 2) % n_pts], \
                               facecolor='red', alpha = 0.2)
            if Aij + Ajk < Aik:
                if graphic: p.plot((pts[i + 1][0],),(pts[i + 1][1],),'go')
                del pts[i+1]
            i += 1
            n_pts = len(pts)
        k += 1
    return n.asarray(pts)

if __name__ == "__main__":

    import scipy.interpolate as interpolate

#    fig = p.figure(figsize=(10,10))

    theta = 2*n.pi*n.random.rand(1000)
    r = n.random.rand(1000)**0.5
    x,y = r*p.cos(theta),r*p.sin(theta)

    points = n.ndarray((2,len(x)))
    points[0,:],points[1,:] = x,y

    scale = 1.03
    hull_pts = scale*convex_hull(points)

    p.plot(x,y,'ko')

    x,y = [],[]
    convex = scale*hull_pts
    for point in convex:
        x.append(point[0])
        y.append(point[1])
    x.append(convex[0][0])
    y.append(convex[0][1])

    x,y = n.array(x),n.array(y)

#Taken from https://stackoverflow.com/questions/14344099/numpy-scipy-smooth-spline-representation-of-an-arbitrary-contour-flength
    nt = n.linspace(0, 1, 100)
    t = n.zeros(x.shape)
    t[1:] = n.sqrt((x[1:] - x[:-1])**2 + (y[1:] - y[:-1])**2)
    t = n.cumsum(t)
    t /= t[-1]
    x2 = interpolate.spline(t, x, nt)
    y2 = interpolate.spline(t, y, nt)
    p.plot(x2, y2,'r--',linewidth=2)

    p.show()

在此处输入图像描述

有一些有用的论文,例如:

http://repositorium.sdum.uminho.pt/bitstream/1822/6429/1/ConcaveHull_ACM_MYS.pdf

此外,您可以尝试: http: //resources.arcgis.com/en/help/main/10.1/index.html#//007000000013000000

我对arcgis一无所知,但它看起来不错。

于 2013-07-09T20:43:50.837 回答
0

我遇到了这个并实现了易于使用的功能以及一些替代方案/改进。

改进:

  • 使用确保平滑的周期性插值
  • 使用二次插值
  • 现在也只适用于正点
  • scipy.interpolate.spline使用已弃用功能的替代方案

备择方案:

  • 许多不同且可配置的插值方案
  • 圆角凸包版本

希望这对一路上的人有所帮助。

import sklearn.preprocessing
import sklearn.pipeline
import scipy.spatial
import numpy as np


def calculate_hull(
        X, 
        scale=1.1, 
        padding="scale", 
        n_interpolate=100, 
        interpolation="quadratic_periodic", 
        return_hull_points=False):
    """
    Calculates a "smooth" hull around given points in `X`.
    The different settings have different drawbacks but the given defaults work reasonably well.
    Parameters
    ----------
    X : np.ndarray
        2d-array with 2 columns and `n` rows
    scale : float, optional
        padding strength, by default 1.1
    padding : str, optional
        padding mode, by default "scale"
    n_interpolate : int, optional
        number of interpolation points, by default 100
    interpolation : str or callable(ix,iy,x), optional
        interpolation mode, by default "quadratic_periodic"

    Inspired by: https://stackoverflow.com/a/17557853/991496
    """
    
    if padding == "scale":

        # scaling based padding
        scaler = sklearn.pipeline.make_pipeline(
            sklearn.preprocessing.StandardScaler(with_std=False),
            sklearn.preprocessing.MinMaxScaler(feature_range=(-1,1)))
        points_scaled = scaler.fit_transform(X) * scale
        hull_scaled = scipy.spatial.ConvexHull(points_scaled, incremental=True)
        hull_points_scaled = points_scaled[hull_scaled.vertices]
        hull_points = scaler.inverse_transform(hull_points_scaled)
        hull_points = np.concatenate([hull_points, hull_points[:1]])
    
    elif padding == "extend" or isinstance(padding, (float, int)):
        # extension based padding
        # TODO: remove?
        if padding == "extend":
            add = (scale - 1) * np.max([
                X[:,0].max() - X[:,0].min(), 
                X[:,1].max() - X[:,1].min()])
        else:
            add = padding
        points_added = np.concatenate([
            X + [0,add], 
            X - [0,add], 
            X + [add, 0], 
            X - [add, 0]])
        hull = scipy.spatial.ConvexHull(points_added)
        hull_points = points_added[hull.vertices]
        hull_points = np.concatenate([hull_points, hull_points[:1]])
    else:
        raise ValueError(f"Unknown padding mode: {padding}")
    
    # number of interpolated points
    nt = np.linspace(0, 1, n_interpolate)
    
    x, y = hull_points[:,0], hull_points[:,1]
    
    # ensures the same spacing of points between all hull points
    t = np.zeros(x.shape)
    t[1:] = np.sqrt((x[1:] - x[:-1])**2 + (y[1:] - y[:-1])**2)
    t = np.cumsum(t)
    t /= t[-1]

    # interpolation types
    if interpolation is None or interpolation == "linear":
        x2 = scipy.interpolate.interp1d(t, x, kind="linear")(nt)
        y2 = scipy.interpolate.interp1d(t, y, kind="linear")(nt)
    elif interpolation == "quadratic":
        x2 = scipy.interpolate.interp1d(t, x, kind="quadratic")(nt)
        y2 = scipy.interpolate.interp1d(t, y, kind="quadratic")(nt)

    elif interpolation == "quadratic_periodic":
        x2 = scipy.interpolate.splev(nt, scipy.interpolate.splrep(t, x, per=True, k=4))
        y2 = scipy.interpolate.splev(nt, scipy.interpolate.splrep(t, y, per=True, k=4))
    
    elif interpolation == "cubic":
        x2 = scipy.interpolate.CubicSpline(t, x, bc_type="periodic")(nt)
        y2 = scipy.interpolate.CubicSpline(t, y, bc_type="periodic")(nt)
    else:
        x2 = interpolation(t, x, nt)
        y2 = interpolation(t, y, nt)
    
    X_hull = np.concatenate([x2.reshape(-1,1), y2.reshape(-1,1)], axis=1)
    if return_hull_points:
        return X_hull, hull_points
    else:
        return X_hull


def draw_hull(
        X, 
        scale=1.1, 
        padding="scale", 
        n_interpolate=100, 
        interpolation="quadratic_periodic",
        plot_kwargs=None, 
        ax=None):
    """Uses `calculate_hull` to draw a hull around given points.

    Parameters
    ----------
    X : np.ndarray
        2d-array with 2 columns and `n` rows
    scale : float, optional
        padding strength, by default 1.1
    padding : str, optional
        padding mode, by default "scale"
    n_interpolate : int, optional
        number of interpolation points, by default 100
    interpolation : str or callable(ix,iy,x), optional
        interpolation mode, by default "quadratic_periodic"
    plot_kwargs : dict, optional
        `matplotlib.pyplot.plot` kwargs, by default None
    ax : `matplotlib.axes.Axes`, optional
        [description], by default None
    """

    if plot_kwargs is None:
        plot_kwargs = {}

    X_hull = calculate_hull(
        X, scale=scale, padding=padding, n_interpolate=n_interpolate, interpolation=interpolation)
    if ax is None:
        ax= plt.gca()
    plt.plot(X_hull[:,0], X_hull[:,1], **plot_kwargs)


def draw_rounded_hull(X, padding=0.1, line_kwargs=None, ax=None):
    """Plots a convex hull around points with rounded corners and a given padding.

    Parameters
    ----------
    X : np.array
        2d array with two columns and n rows
    padding : float, optional
        padding between hull and points, by default 0.1
    line_kwargs : dict, optional
        line kwargs (used for `matplotlib.pyplot.plot` and `matplotlib.patches.Arc`), by default None
    ax : matplotlib.axes.Axes, optional
        axes to plat on, by default None
    """

    default_line_kwargs = dict(
        color="black",
        linewidth=1
    )
    if line_kwargs is None:
        line_kwargs = default_line_kwargs
    else:
        line_kwargs = {**default_line_kwargs, **line_kwargs}

    if ax is None:
        ax = plt.gca()

    hull = scipy.spatial.ConvexHull(X)
    hull_points = X[hull.vertices]

    hull_points = np.concatenate([hull_points[[-1]], hull_points, hull_points[[0]]])

    diameter = padding * 2
    for i in range(1, hull_points.shape[0] - 1):

        # line
        
        # source: https://stackoverflow.com/a/1243676/991496
        
        norm_next = np.flip(hull_points[i] - hull_points[i + 1]) * [-1, 1]
        norm_next /= np.linalg.norm(norm_next)

        norm_prev = np.flip(hull_points[i - 1] - hull_points[i]) * [-1, 1]
        norm_prev /= np.linalg.norm(norm_prev)

        # plot line
        line = hull_points[i:i+2] + norm_next * diameter / 2
        ax.plot(line[:,0], line[:,1], **line_kwargs) 

        # arc

        angle_next = np.rad2deg(np.arccos(np.dot(norm_next, [1,0])))
        if norm_next[1] < 0:
            angle_next = 360 - angle_next

        angle_prev = np.rad2deg(np.arccos(np.dot(norm_prev, [1,0])))
        if norm_prev[1] < 0:
            angle_prev = 360 - angle_prev

        arc = patches.Arc(
            hull_points[i], 
            diameter, diameter,
            angle=0, fill=False, theta1=angle_prev, theta2=angle_next,
            **line_kwargs)

        ax.add_patch(arc)


if __name__ == '__main__':

    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib import patches

    # np.random.seed(42)
    X = np.random.random((20,2))

    fig, ax = plt.subplots(1,1, figsize=(10,10))
    ax.scatter(X[:,0], X[:,1])
    draw_rounded_hull(X, padding=0.1)
    draw_hull(X)
    
    ax.set(xlim=[-1,2], ylim= [-1,2])
    fig.savefig("_out/test.png")

在此处输入图像描述

于 2022-03-05T08:58:58.677 回答