python - 为什么 SymPy 会计算出错误的平面交点？

Question

我有一个奇怪的问题，SymPy 中的平面相交适用于简单的示例，但对于具有更复杂坐标的示例却失败了。我发布了一个有效的简单示例和一个失败的示例。正如 Povray 图像所示，我有三个平面穿过 a polyhedron 的顶点。并且垂直于通过相应顶点和中心的线。我想计算这些平面相交的点，但是 SymPy 给出了平面对相交的直线的错误结果。在图像中，正确的交点可以看作是短线(使用 CSG 交点创建)。与它们平行的长线是SymPy计算的。

我是不是做错了什么，或者这是 SymPy 中的错误？

更多图片在这里:http://paste.watchduck.net/1712/sympy_planes/
有谁知道如何在页面上放置许多图像，而不会被阻止发布问题？ (“您的帖子似乎包含格式不正确的代码。”)

作品

代码:

from sympy import Point3D, Plane


pointR = Point3D(1/2, 0, 1/2)
pointG = Point3D(1, 0, 0)

planeR = Plane(pointR, pointR)
planeG = Plane(pointG, pointG)

print('\n######## Intersection of the planes:')
lineRG = planeR.intersection(planeG)[0]  # yellow
print(lineRG)

print('\n######## Intersection of plane and contained line returns the line:')
lineRG_again = planeR.intersection(lineRG)[0]
print(lineRG_again.equals(lineRG))

输出:

######## Intersection of the planes:
Line3D(Point3D(1, 0, 0), Point3D(1, 1/2, 0))

######## Intersection of plane and contained line returns the line:
True

失败

代码:

from sympy import sqrt, Point3D, Plane

pointR = Point3D(-1, 1 + sqrt(2), -2*sqrt(2) - 1)
pointG = Point3D(-sqrt(2) - 1, 1, -2*sqrt(2) - 1)
pointB = Point3D(-1, -sqrt(2) - 1, -2*sqrt(2) - 1)

planeR = Plane(pointR, pointR)
planeG = Plane(pointG, pointG)
planeB = Plane(pointB, pointB)

print('\n######## Intersections of the planes:')

lineRG = planeR.intersection(planeG)[0]  # yellow
lineRB = planeR.intersection(planeB)[0]  # magenta
lineGB = planeG.intersection(planeB)[0]  # cyan

print(lineRG)
print(lineRB)
print(lineGB)

print('\n######## Lines RG (yellow) and GB (cyan) intersect:')
print(lineRG.intersection(lineGB))
print('\n######## But line RB (magenta) is skew to both of them:')
print(lineRB.intersection(lineRG))
print(lineRB.intersection(lineGB))

print('\n######## Intersection of plane and contained line fails:')
lineRG_again = planeR.intersection(lineRG)

输出:

######## Intersections of the planes:
Line3D(Point3D(-1, 1, 0), Point3D((1 + sqrt(2))*(-2*sqrt(2) - 1) + 2*sqrt(2), -2*sqrt(2) + (-2*sqrt(2) - 1)*(-sqrt(2) - 1), -1 - (1 + sqrt(2))*(-sqrt(2) - 1)))
Line3D(Point3D(-1, 0, 0), Point3D((1 + sqrt(2))*(-2*sqrt(2) - 1) - (-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 1, 0, 2 + 2*sqrt(2)))
Line3D(Point3D(-1, 1, 0), Point3D(-(-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 2*sqrt(2) - 2, -(-2*sqrt(2) - 1)*(-sqrt(2) - 1) + 2 + 2*sqrt(2), 1 + (-sqrt(2) - 1)**2))

######## Lines RG (yellow) and GB (cyan) intersect:
[Point3D(-1, 1, 0)]

######## But line RB (magenta) is skew to both of them:
[]
[]

######## Intersection of plane and contained line fails:
Traceback (most recent call last):
  File "planes2.py", line 47, in <module>
    lineRG_again = planeR.intersection(lineRG)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/geometry/plane.py", line 390, in intersection
    p = a.subs(t, c[0])
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/basic.py", line 916, in subs
    rv = rv._subs(old, new, **kwargs)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/cache.py", line 93, in wrapper
    retval = cfunc(*args, **kwargs)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/basic.py", line 1030, in _subs
    rv = fallback(self, old, new)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/basic.py", line 1007, in fallback
    rv = self.func(*args)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/geometry/point.py", line 1104, in __new__
    args = Point(*args, **kwargs)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/geometry/point.py", line 159, in __new__
    raise ValueError('Imaginary coordinates are not permitted.')
ValueError: Imaginary coordinates are not permitted.

图片:

编辑:适用于 SymPy 1.1.2

安装 SymPy 的开发版本后(pip install git+https://github.com/sympy/sympy.git)我得到了正确的结果:

######## Intersections of pairs of planes:
Line3D(Point3D(-7 + sqrt(2)/2, -sqrt(2)/2 + 7, 0), Point3D((1 + sqrt(2))*(-2*sqrt(2) - 1) - 6 + 5*sqrt(2)/2, -5*sqrt(2)/2 + 6 + (-2*sqrt(2) - 1)*(-sqrt(2) - 1), -1 - (1 + sqrt(2))*(-sqrt(2) - 1)))
Line3D(Point3D(-13 - 6*sqrt(2), 0, 0), Point3D(-13 + (1 + sqrt(2))*(-2*sqrt(2) - 1) - (-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 6*sqrt(2), 0, 2 + 2*sqrt(2)))
Line3D(Point3D(-13/2 - 3*sqrt(2), -7*sqrt(2)/2 + 1/2, 0), Point3D(-(-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 15/2 - 5*sqrt(2), -(-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 3*sqrt(2)/2 + 3/2, 1 + (-sqrt(2) - 1)**2))

######## Intersection of all planes:
[Point3D(0, 0, -20*sqrt(2)/7 - 11/7)]

Answer 1

在 SymPy 1.1.1 及更早版本中，intersection 在法向量涉及部首时返回错误结果。这是一个更简单的示例:

p1 = Plane((1, 0, 0), (sqrt(2), 0, 0))
p2 = Plane((1, 1, 1), (1, 1, 1))
line = p1.intersection(p2)[0]      # this line is wrong
print(line.arbitrary_point())

返回线的参数方程为 Point3D(3, -sqrt(2)*t, sqrt(2)*t) 这是错误的，因为第一个平面有方程 sqrt(2 )*(x-1) = 0，即 x=1。

你仍然可以找到正确的交集方程

solve([p1.equation(), p2.equation()])

但这不会那么容易用于绘图。

好消息

该错误(在 linsolve 方法中)已在当前开发版本 1.1.2.dev 中修复。从GitHub repo获取.

早期版本的解决方法

用 float 替换部首:

pointR = Point3D(-1, N(1 + sqrt(2)), N(-2*sqrt(2) - 1))
pointG = Point3D(N(-sqrt(2) - 1), 1, N(-2*sqrt(2) - 1))
pointB = Point3D(-1, N(-sqrt(2) - 1), N(-2*sqrt(2) - 1))

这不会使一切变得完美，但应该会减少错误的影响，并且您可能能够为您的图表获得合理的交叉点。