【问题标题】:Plotting implicit equations in 3d在 3d 中绘制隐式方程
【发布时间】:2020-09-08 14:28:35
【问题描述】:

我想在 3D 中绘制隐式方程 F(x,y,z) = 0。在 Matplotlib 中可以吗?

【问题讨论】:

  • 您可以在以下位置找到一些示例:matplotlib.sourceforge.net/examples/mplot3d/index.html
  • 你需要用matplotlib做吗?如果没有,您可能想看看3d contour plots in Mayavi
  • @Sven Marnach:谢谢,但不幸的是我必须使用 Matplotlib。
  • 看起来不错,只是您没有将 z 轮廓的位置移动 z 的值。它们都以 0.0 绘制!另外,您需要手动定义绘制的限制,因为 z 等高线解决方案将超出您所需的等高线间隔。
  • 我用更简洁的代码更新了我的帖子,其中包括沿其他轴绘制的等高线间隔(切片)。

标签: python numpy matplotlib plot


【解决方案1】:

您可以欺骗 matplotlib 以 3D 形式绘制隐式方程。只需为所需范围内的每个 z 值制作方程的一级等高线图。您也可以沿 y 轴和 z 轴重复该过程,以获得更立体的形状。

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np

def plot_implicit(fn, bbox=(-2.5,2.5)):
    ''' create a plot of an implicit function
    fn  ...implicit function (plot where fn==0)
    bbox ..the x,y,and z limits of plotted interval'''
    xmin, xmax, ymin, ymax, zmin, zmax = bbox*3
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    A = np.linspace(xmin, xmax, 100) # resolution of the contour
    B = np.linspace(xmin, xmax, 15) # number of slices
    A1,A2 = np.meshgrid(A,A) # grid on which the contour is plotted

    for z in B: # plot contours in the XY plane
        X,Y = A1,A2
        Z = fn(X,Y,z)
        cset = ax.contour(X, Y, Z+z, [z], zdir='z')
        # [z] defines the only level to plot for this contour for this value of z

    for y in B: # plot contours in the XZ plane
        X,Z = A1,A2
        Y = fn(X,y,Z)
        cset = ax.contour(X, Y+y, Z, [y], zdir='y')

    for x in B: # plot contours in the YZ plane
        Y,Z = A1,A2
        X = fn(x,Y,Z)
        cset = ax.contour(X+x, Y, Z, [x], zdir='x')

    # must set plot limits because the contour will likely extend
    # way beyond the displayed level.  Otherwise matplotlib extends the plot limits
    # to encompass all values in the contour.
    ax.set_zlim3d(zmin,zmax)
    ax.set_xlim3d(xmin,xmax)
    ax.set_ylim3d(ymin,ymax)

    plt.show()

Goursat Tangle 的情节如下:

def goursat_tangle(x,y,z):
    a,b,c = 0.0,-5.0,11.8
    return x**4+y**4+z**4+a*(x**2+y**2+z**2)**2+b*(x**2+y**2+z**2)+c

plot_implicit(goursat_tangle)

您可以通过使用创意色彩映射添加深度提示来更轻松地进行可视化:

OP 的情节如下所示:

def hyp_part1(x,y,z):
    return -(x**2) - (y**2) + (z**2) - 1

plot_implicit(hyp_part1, bbox=(-100.,100.))

奖励:您可以使用 python 在功能上组合这些隐式函数:

def sphere(x,y,z):
    return x**2 + y**2 + z**2 - 2.0**2

def translate(fn,x,y,z):
    return lambda a,b,c: fn(x-a,y-b,z-c)

def union(*fns):
    return lambda x,y,z: np.min(
        [fn(x,y,z) for fn in fns], 0)

def intersect(*fns):
    return lambda x,y,z: np.max(
        [fn(x,y,z) for fn in fns], 0)

def subtract(fn1, fn2):
    return intersect(fn1, lambda *args:-fn2(*args))

plot_implicit(union(sphere,translate(sphere, 1.,1.,1.)), (-2.,3.))

【讨论】:

  • @bpowah:根据 matplotlib 参考投影 = '3d' 无效。也许投影='直线'?
  • 我记得在升级到 matplotlib 1.0.1 之前遇到了这个错误。我忘记了如何解决以前的版本。
  • @bpowah:哦,您使用的是较新的 matplotlib 版本。无论如何,我会尝试解决这个问题。谢谢))
  • @bpowah:我正在尝试查看您的代码,但我不明白 Z = sphere(X,Y,z) 的工作原理。
  • @Paul 您使用的创意色彩映射是什么?
【解决方案2】:

更新:我终于找到了一种使用matplotlibscikit-image 渲染3D 隐式表面的简单方法,请参阅我的other answer。我把这个留给了对绘制参数化 3D 表面感兴趣的人。

动机

迟到的答案,我只需要做同样的事情,我在某种程度上找到了另一种方法。所以我分享这个另一种观点。

本帖不回答:(1)如何绘制任何隐函数F(x,y,z)=0?但确实回答了:(2) 如何使用带有matplotlib 的网格绘制参数曲面(不是所有隐式函数,而是其中一些函数)?

@Paul 的方法具有非参数化的优点,因此我们可以在每个轴上使用轮廓方法绘制几乎任何我们想要的东西,它完全解决了 (1)。但是matplotlib 无法通过这种方法轻松构建网格,因此我们无法直接从中获取曲面,而是获取各个方向的平面曲线。这就是我回答的动机,我想解决 (2)。

渲染网格

如果我们能够对我们想要绘制的表面进行参数化(这可能很难或不可能),最多使用 2 个参数,那么我们可以使用matplotlib.plot_trisurf 方法绘制它。

也就是说,从一个隐式方程F(x,y,z)=0,如果我们能够得到一个参数系统S={x=f(u,v), y=g(u,v), z=h(u,v)},那么我们可以用matplotlib轻松绘制它,而不必求助于contour

那么,渲染这样一个 3D 表面归结为:

# Render:
ax = plt.axes(projection='3d')
ax.plot_trisurf(x, y, z, triangles=tri.triangles, cmap='jet', antialiased=True) 

其中(x, y, z) 是向量(不是meshgrid,请参见ravel),从参数(u, v)triangles 函数计算得出,参数是从(u,v) 参数派生的三角剖分,以承担网格构造。

进口

必需的导入是:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
from matplotlib.tri import Triangulation

一些表面

让我们参数化一些表面...

领域
# Parameters:
theta = np.linspace(0, 2*np.pi, 20)
phi = np.linspace(0, np.pi, 20)
theta, phi = np.meshgrid(theta, phi)
rho = 1

# Parametrization:
x = np.ravel(rho*np.cos(theta)*np.sin(phi))
y = np.ravel(rho*np.sin(theta)*np.sin(phi))
z = np.ravel(rho*np.cos(phi))

# Triangulation:
tri = Triangulation(np.ravel(theta), np.ravel(phi))

锥体
theta = np.linspace(0, 2*np.pi, 20)
rho = np.linspace(-2, 2, 20)
theta, rho = np.meshgrid(theta, rho)

x = np.ravel(rho*np.cos(theta))
y = np.ravel(rho*np.sin(theta))
z = np.ravel(rho)

tri = Triangulation(np.ravel(theta), np.ravel(rho))

环面
a, c = 1, 4
u = np.linspace(0, 2*np.pi, 20)
v = u.copy()
u, v = np.meshgrid(u, v)

x = np.ravel((c + a*np.cos(v))*np.cos(u))
y = np.ravel((c + a*np.cos(v))*np.sin(u))
z = np.ravel(a*np.sin(v))

tri = Triangulation(np.ravel(u), np.ravel(v))

Möbius Strip
u = np.linspace(0, 2*np.pi, 20)
v = np.linspace(-1, 1, 20)
u, v = np.meshgrid(u, v)

x = np.ravel((2 + (v/2)*np.cos(u/2))*np.cos(u))
y = np.ravel((2 + (v/2)*np.cos(u/2))*np.sin(u))
z = np.ravel(v/2*np.sin(u/2))

tri = Triangulation(np.ravel(u), np.ravel(v))

限制

大多数时候,Triangulation 是必需的,以便协调plot_trisurf 方法的网格构造,并且该对象只接受两个参数,因此我们仅限于二维参数曲面。我们不太可能用这种方法来表示 Goursat Tangle。

【讨论】:

  • @Paul,我找到了一种使用tri_surfTriangulation 将2d 参数化3d 函数表示为表面的方法,你知道吗?我想知道是否有办法从你的轮廓构建一个有效的网格。
【解决方案3】:

Matplotlib 需要一系列点;如果你能弄清楚如何渲染你的方程,它会做绘图。

参考Is it possible to plot implicit equations using Matplotlib? Mike Graham 的回答建议使用 scipy.optimize 对隐式函数进行数值探索。

http://xrt.wikidot.com/gallery:implicit 有一个有趣的画廊,展示了各种光线追踪隐式函数 - 如果您的方程与其中一个匹配,它可能会让您更好地了解您正在查看的内容。

如果你愿意分享实际的方程式,也许有人可以提出更简单的方法。

【讨论】:

    【解决方案4】:

    据我所知,这是不可能的。你必须自己解决这个方程。使用 scipy.optimize 是个好主意。最简单的情况是,您知道要绘制的曲面范围,只需在 x 和 y 中制作规则网格,并尝试求解方程 F(xi,yi,z)=0 的 z,给出一个起点z 点。以下是一个非常肮脏的代码,可能会对您有所帮助

    from scipy import *
    from scipy import optimize
    
    xrange = (0,1)
    yrange = (0,1)
    density = 100
    startz = 1
    
    def F(x,y,z):
        return x**2+y**2+z**2-10
    
    x = linspace(xrange[0],xrange[1],density)
    y = linspace(yrange[0],yrange[1],density)
    
    points = []
    for xi in x:
        for yi in y:
            g = lambda z:F(xi,yi,z)
            res = optimize.fsolve(g, startz, full_output=1)
            if res[2] == 1:
                zi = res[0]
                points.append([xi,yi,zi])
    
    points = array(points)
    

    【讨论】:

    • 我认为这不是一个好主意。尝试 F(x, y, z) = x^2 + y^2 + z^2 - 1 来查看问题。它应该绘制一个球体,但您的代码只会绘制半个球体。
    • 哦,刚刚注意到你已经在使用基本相同的功能了:)
    • 确实如此。如果隐函数是多值的,那就有问题了。正如我所说,这只是一个让最简单的事情成为可能的肮脏代码。#
    【解决方案5】:

    实际上,使用 scikit-image 包可以轻松绘制隐式 3D 曲面。关键是marching_cubes 方法。

    import numpy as np
    from skimage import measure
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import axes3d
    

    然后我们在 3D 网格上计算函数,在这个例子中,我们使用在其答案中定义的 goursat_tangle 方法 @Paul

    xl = np.linspace(-3, 3, 50)
    X, Y, Z = np.meshgrid(xl, xl, xl)
    F = goursat_tangle(X, Y, Z)
    

    marching_cubes 在这里发生了奇迹:

    verts, faces, normals, values = measure.marching_cubes(F, 0, spacing=[np.diff(xl)[0]]*3)
    verts -= 3
    

    我们只需要校正以Voxel 坐标表示的顶点坐标(因此使用spacing 开关和随后的原点移位进行缩放)。

    最后只是使用tri_surface 渲染等值面:

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_trisurf(verts[:, 0], verts[:, 1], faces, verts[:, 2], cmap='jet', lw=0)
    

    返回:

    【讨论】:

    • @qutron 和@Paul,我找到了一种使用matplotlibscikit-image 绘制隐式3D 表面的简单方法。干杯
    • 必须添加 plot.show() 才能使其工作。我还建议您添加 Goursat Tangle 或其他一些功能以确保完整性,并考虑将它们全部添加到一个代码块中以便于复制粘贴。
    • 太棒了。我认为x轴和y轴可以互换。对于非对称函数ax.plot_trisurf(verts[:, 1], verts[:, 0], faces, verts[:, 2], cmap='jet', lw=0) 符合我的期望。
    【解决方案6】:

    你看过 matplotlib 上的mplot3d 吗?

    【讨论】:

    • 是的,我有。主要问题是函数是隐式的。 AFAIK,Matplotlib 不绘制方程,它绘制一系列点。我只是不知道如何计算适合我的隐式方程的所有点系列。
    【解决方案7】:

    最后,我做到了(我将 matplotlib 更新到 1.0.1)。 这是代码:

    import matplotlib.pyplot as plt
    import numpy as np
    from mpl_toolkits.mplot3d import Axes3D
    
    def hyp_part1(x,y,z):
        return -(x**2) - (y**2) + (z**2) - 1
    
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    
    x_range = np.arange(-100,100,10) 
    y_range = np.arange(-100,100,10)
    X,Y = np.meshgrid(x_range,y_range)
    A = np.linspace(-100, 100, 15)
    
    A1,A2 = np.meshgrid(A,A)    
    
    for z in A: 
        X,Y = A1, A2
        Z = hyp_part1(X,Y,z)
        ax.contour(X, Y, Z+z, [z], zdir='z')
    
    for y in A: 
        X,Z= A1, A2
        Y = hyp_part1(X,y,Z)
        ax.contour(X, Y+y, Z, [y], zdir='y')
    
    for x in A:
        Y,Z = A1, A2 
        X = hyp_part1(x,Y,Z)
        ax.contour(X+x, Y, Z, [x], zdir='x')
    
    ax.set_zlim3d(-100,100)
    ax.set_xlim3d(-100,100)
    ax.set_ylim3d(-100,100)
    

    这是结果:

    谢谢你,保罗!

    【讨论】:

      【解决方案8】:

      MathGL(GPL 绘图库)可以轻松绘制它。只需创建一个具有函数值 f[i,j,k] 的数据网格,并使用 Surf3() 函数在值 f[i,j,k]=0 处生成等值面。看到这个sample

      【讨论】:

      • @链接示例已失效。
      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 2015-05-17
      • 1970-01-01
      • 2022-11-26
      • 2022-11-10
      • 2023-02-05
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多