为什么SymPy计算平面交点时会出错？

Question

我有一个奇怪的问题，SymPy中平面的交点可以处理简单的例子，但在坐标更复杂的情况下会失败。我发布了一个可行的简单示例和一个失败的示例。如Povray图像所示，我有三个平面通过多面体的顶点，并且垂直于通过相应顶点和中心的线。我想计算这些平面相交的点，但是SymPy给出的相交线是错误的。在图像中，可以看到正确的交点是短线（使用CSG相交创建）。与它们平行的长线是SymPy计算出的线。

我做错了什么，还是这是SymPy的一个bug？

更多图像在这里：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)]

- Watchduck

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()])

但是它在绘图中并不容易使用。

好消息

当前开发版本1.1.2.dev已修复了linsolve方法中的错误。从GitHub repo获取它。

早期版本的解决方法

用浮点数替换根式:

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))

这并不能完全解决问题，但它可以减小bug的影响，让你能够得到合理的图表交集。

- user6655984

Answer 2

这个bug在正在开发的版本中没有完全解决。特别是涉及到非简单坐标线和圆之间的交点。您可以使用shapely代替：

[在Python中找到线和圆的交点最有效的方法是什么？]

- Ella