使用网格生成三角形轮廓

0 投票
1 回答
3535 浏览
提问于 2025-04-16 14:58

我正在尝试用 matplotlib 和 numpy 重现维基百科上关于 Dirichlet 分布 的图表(或者只要一个轮廓图也可以)。我在生成三角形的等高线图时遇到了一些困难。第一个问题是 meshgrid 这个函数不能返回一个三角形的点。如果我能得到一个三角形的点,contourf 能处理这种不规则的输入吗?

这是我目前的代码:

#!/usr/bin/env python
from __future__ import division
import matplotlib
matplotlib.use("TkAgg")
matplotlib.rc('text', usetex=True)
matplotlib.rcParams['text.latex.preamble']=r"""\usepackage{amsmath}
"""
import math
import scipy.special

root_three_over_two = np.sqrt(3) / 2

def new_figure():
    # 1.45
    plt.figure(figsize = [2.6, 2.6 * root_three_over_two], dpi = 1200)
    plt.axes([0.05, 0.10, 0.90, 0.90], frameon = False)

xsize = 1.0
ysize = root_three_over_two * xsize
plt.axis([0, xsize, 0, ysize])
resolution = 0.05

R = inclusive_arange(0.0, 1.0, resolution)
x, y = np.meshgrid(inclusive_arange(0.0, 1.0, resolution),
                   inclusive_arange(0.0, 1.0, resolution))
# UNFORTUNATELY x, and y include a lot of points where x+y>1
x = []
y = []
for yy in R:
    x.append(list(inclusive_arange(0.0, 1.0 - yy, resolution)))
    y.append([yy for xx in R])
print x
print y
z = 1 - x - y

# We can use these to convert to and from the equilateral triangle.
M = [[1, 0.5], [0, root_three_over_two]]
Mi = np.linalg.inv(M)

def dirichlet(x, y, z, a, b, c):
    if z < 0:
        return 0
    return x ** (a - 1) * y ** (b - 1) * z ** (c - 1) \
            * math.gamma(a + b + c) \
            / (math.gamma(a) * math.gamma(b) * math.gamma(c))

dirichlet = np.frompyfunc(dirichlet, 6, 1)

for (dirichlet_parm, filename) in [((5.0, 1.5, 2.5), "dir_small.pdf")]:
    new_figure()
    height = dirichlet(x, y, z, *dirichlet_parm)
    M = np.max(height)
    cs = plt.contourf(x, y, height, 50)
    S = sum(dirichlet_parm)
    plt.savefig(filename)

1 个回答

2

你并不需要一个规则的网格来制作等高线图。你可以先定义一个边界,然后进行德劳内三角剖分。坐标当然还是要保持规则的,看起来你可能需要进行一些转换。这个例子应该能帮助你创建一个二维等高线图。我也制作过一些边界不规则的三维表面图,这样会出现一些边缘的瑕疵,看起来不太好。使用一个非常细的网格可能会改善这些问题。

撰写回答