给定非平滑基础数据,如何在 matplotlib 中获得平滑轮廓线
How to get a smoothed contour line in matplotlib given non-smooth underlying data
我在二维网格上有布尔数据,想使用 matplotlib
在数据为 True
和 False
的区域之间绘制一个轮廓图。
但是,这些区域之间的分离在实际数据中并不平滑。给定此数据,我如何计算平滑的计数?
这是一个最小的例子:
import numpy as np
import matplotlib.pyplot as plt
# generate some non-smooth example data
MESHSIZE = 10
REFINEMENT = 4*MESHSIZE
x = np.linspace(-MESHSIZE, MESHSIZE, REFINEMENT)
xv, yv = np.meshgrid(x, x)
xvf = xv.reshape(-1)
yvf = yv.reshape(-1)
def choppy_circle(x, y):
inner = x.astype(int)**2+y.astype(int)**2 < 10.0
return inner
# consider this the *actual* data given to me as-is
my_x = xvf
my_y = yvf
my_z = choppy_circle(xvf, yvf)
# need to visualize the contour that separates areas where
# my_z is True/False
plt.tricontour(my_x, my_y, my_z, levels=np.array([1.0-1e-3]))
plt.scatter(xv, yv, s=0.1)
plt.show()
这产生了以下情节,它忠实于数据,但不是我要找的:
如何使用 my_x
、my_y
和 my_z
中给出的数据在 my_z
为 True
的域周围构建平滑轮廓?
像这样:
您可以将样条曲线拟合到轮廓。并通过选择样条曲线的平滑参数使其尽可能平滑。
首先获取边界点
import functools
import itertools
mask = my_z.reshape(40,40)
mask &= functools.reduce(np.logical_or,[~np.roll(np.roll(mask, shift_x, 0),shift_y,1)
for shift_x,shift_y in itertools.product((-1,0,1),repeat=2)])
x,y = my_x[mask.reshape(-1)],my_y[mask.reshape(-1)]
plt.scatter(x,y)
现在我们根据相应复数的自变量对您的分数进行排序。如果您不明白我的意思,那就是点与原点和点 (1,0) 所成的角度。并为其拟合一条样条曲线。
import scipy.interpolate as interpolate
import matplotlib.pyplot as plt
arr = np.array(sorted(zip(x,y), key=lambda x: cmath.phase(x[0]+1j*x[1])))
s=1
tck, u = interpolate.splprep([arr[:,0],arr[:,1]],per=1, s=s)
x_i, y_i = interpolate.splev(np.linspace(0, 1, 10**4), tck)
ax = plt.gca()
ax.plot(x_i, y_i)
ax.scatter(arr[:,0],arr[:,1])
ax.set_title(f"{s=}")
ax.set_aspect('equal')
根据 s
,结果看起来会有所不同。我为你绘制了一些:
您可以使用shapely
获取任意形状的质心和边界框,然后绘制一个圆:
# […] same as previously
# get points
cs = plt.tricontour(my_x, my_y, my_z, levels=np.array([1.0-1e-3]))
v = cs.collections[0].get_paths()[0].vertices
from shapely.geometry import Polygon
# find centroid coordinates and bounding box
p = Polygon(v)
x,y =p.centroid.coords[0]
minx, miny, maxx, maxy = p.bounds
# plot circle
# depending on the data, one could also plot an ellipse or rectangle
r = max((maxx-minx)/2, (maxy-miny)/2)
circle = plt.Circle((x, y), r, color='r', fill=False)
plt.gca().add_patch(circle)
输出:
像this答案中提出的那样提取轮廓数据,并按照@user2640045提出的样条插值允许对任意轮廓进行操作:
# my_x, my_y, my_z as above...
# get contour data
cs = plt.tricontour(my_x, my_y, my_z, levels=np.array([1.0-1e-3]))
print(type(cs))
# process each contour
for contour in cs.collections[0].get_paths():
vert = contour.vertices
vert_x = vert[:, 0]
vert_y = vert[:, 1]
# plot contour
plt.plot(vert_x, vert_y)
# interpolate contour points
s = 20
tck, u = interpolate.splprep([vert_x, vert_y], per=1, s=s)
x_interp, y_interp = interpolate.splev(np.linspace(0, 1, 10**3), tck)
# plot interpolated contour
plt.plot(x_interp, y_interp)
# plot grid
plt.scatter(xv, yv, s=0.1)
# display plot
plt.show()
重要的一点是循环头
for contour in cs.collections[0].get_paths():
其中获取每条等高线的x-y数据
我在二维网格上有布尔数据,想使用 matplotlib
在数据为 True
和 False
的区域之间绘制一个轮廓图。
但是,这些区域之间的分离在实际数据中并不平滑。给定此数据,我如何计算平滑的计数?
这是一个最小的例子:
import numpy as np
import matplotlib.pyplot as plt
# generate some non-smooth example data
MESHSIZE = 10
REFINEMENT = 4*MESHSIZE
x = np.linspace(-MESHSIZE, MESHSIZE, REFINEMENT)
xv, yv = np.meshgrid(x, x)
xvf = xv.reshape(-1)
yvf = yv.reshape(-1)
def choppy_circle(x, y):
inner = x.astype(int)**2+y.astype(int)**2 < 10.0
return inner
# consider this the *actual* data given to me as-is
my_x = xvf
my_y = yvf
my_z = choppy_circle(xvf, yvf)
# need to visualize the contour that separates areas where
# my_z is True/False
plt.tricontour(my_x, my_y, my_z, levels=np.array([1.0-1e-3]))
plt.scatter(xv, yv, s=0.1)
plt.show()
这产生了以下情节,它忠实于数据,但不是我要找的:
如何使用 my_x
、my_y
和 my_z
中给出的数据在 my_z
为 True
的域周围构建平滑轮廓?
像这样:
您可以将样条曲线拟合到轮廓。并通过选择样条曲线的平滑参数使其尽可能平滑。
首先获取边界点
import functools
import itertools
mask = my_z.reshape(40,40)
mask &= functools.reduce(np.logical_or,[~np.roll(np.roll(mask, shift_x, 0),shift_y,1)
for shift_x,shift_y in itertools.product((-1,0,1),repeat=2)])
x,y = my_x[mask.reshape(-1)],my_y[mask.reshape(-1)]
plt.scatter(x,y)
现在我们根据相应复数的自变量对您的分数进行排序。如果您不明白我的意思,那就是点与原点和点 (1,0) 所成的角度。并为其拟合一条样条曲线。
import scipy.interpolate as interpolate
import matplotlib.pyplot as plt
arr = np.array(sorted(zip(x,y), key=lambda x: cmath.phase(x[0]+1j*x[1])))
s=1
tck, u = interpolate.splprep([arr[:,0],arr[:,1]],per=1, s=s)
x_i, y_i = interpolate.splev(np.linspace(0, 1, 10**4), tck)
ax = plt.gca()
ax.plot(x_i, y_i)
ax.scatter(arr[:,0],arr[:,1])
ax.set_title(f"{s=}")
ax.set_aspect('equal')
根据 s
,结果看起来会有所不同。我为你绘制了一些:
您可以使用shapely
获取任意形状的质心和边界框,然后绘制一个圆:
# […] same as previously
# get points
cs = plt.tricontour(my_x, my_y, my_z, levels=np.array([1.0-1e-3]))
v = cs.collections[0].get_paths()[0].vertices
from shapely.geometry import Polygon
# find centroid coordinates and bounding box
p = Polygon(v)
x,y =p.centroid.coords[0]
minx, miny, maxx, maxy = p.bounds
# plot circle
# depending on the data, one could also plot an ellipse or rectangle
r = max((maxx-minx)/2, (maxy-miny)/2)
circle = plt.Circle((x, y), r, color='r', fill=False)
plt.gca().add_patch(circle)
输出:
像this答案中提出的那样提取轮廓数据,并按照@user2640045提出的样条插值允许对任意轮廓进行操作:
# my_x, my_y, my_z as above...
# get contour data
cs = plt.tricontour(my_x, my_y, my_z, levels=np.array([1.0-1e-3]))
print(type(cs))
# process each contour
for contour in cs.collections[0].get_paths():
vert = contour.vertices
vert_x = vert[:, 0]
vert_y = vert[:, 1]
# plot contour
plt.plot(vert_x, vert_y)
# interpolate contour points
s = 20
tck, u = interpolate.splprep([vert_x, vert_y], per=1, s=s)
x_interp, y_interp = interpolate.splev(np.linspace(0, 1, 10**3), tck)
# plot interpolated contour
plt.plot(x_interp, y_interp)
# plot grid
plt.scatter(xv, yv, s=0.1)
# display plot
plt.show()
重要的一点是循环头
for contour in cs.collections[0].get_paths():
其中获取每条等高线的x-y数据