numpy 在三维中绘制隐式方程

9udxz4iz  于 2023-04-30  发布在  其他
关注(0)|答案(8)|浏览(165)

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

svgewumm

svgewumm1#

您可以通过欺骗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)

您可以通过使用创意色彩Map添加深度线索来简化可视化:

下面是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.))

nxowjjhe

nxowjjhe2#

更新:我终于找到了一个简单的方法来渲染3D隐式曲面,看看我的other answer。我离开了这一个谁是有兴趣在绘制参数化三维曲面。

动机

迟来的回答,我只是需要做同样的事情,我找到了另一种方式来做它在某种程度上。所以我分享另一个视角。
这篇文章没有回答:(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)是根据参数(u, v)功能计算的矢量(不是meshgrid,参见ravel),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缠结。

6rqinv9w

6rqinv9w3#

实际上,有一种简单的方法可以使用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)

它返回:

gab6jxml

gab6jxml4#

Matplotlib期望一系列点;如果你能弄清楚如何渲染你的方程,它就会做绘图。
参考Is it possible to plot implicit equations using Matplotlib?,Mike Graham的回答建议使用scipy。优化以数值地探索隐函数。
http://xrt.wikidot.com/gallery:implicit上有一个有趣的图库,展示了各种光线跟踪隐式函数-如果您的方程与其中一个相匹配,它可能会给予您更好地了解您正在查看的内容。
如果做不到这一点,如果你愿意分享实际的等式,也许有人可以提出一个更简单的方法。

dzhpxtsq

dzhpxtsq5#

据我所知,这是不可能的。你得自己用数值方法解这个方程。使用scipy。优化是个好主意。最简单的情况是,你知道你想要绘制的曲面的范围,只需要在x和y上做一个规则的网格,然后尝试求解方程F(xi,yi,z)=0,给出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)
wkftcu5l

wkftcu5l6#

你在matplotlib上看过mplot3d吗?

7lrncoxx

7lrncoxx7#

最后,我做到了(我更新了我的matplotlib到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)

结果如下:

谢谢你,保罗!

7vux5j2d

7vux5j2d8#

MathGL(GPL绘图库)可以很容易地绘图。只需创建具有函数值f[i,j,k]的数据网格,并使用Surf3()函数在值f[i,j,k]=0处生成等值面。看看这个样本。

相关问题