计算两个任意椭圆的交点

作者 : RY    标签: sympy

计算两个任意椭圆的交点需要求解一个一元四次方程,而通过椭圆的参数计算出这个方程的系数也是很麻烦的事情。本文通过Sympy进行符号运算,并自动生成计算交点的程序。

In [20]:
%pylab inline
#%load_ext sympyprinting
%load_ext sympy.interactive.ipythonprinting
%reload_ext sympy.interactive.ipythonprinting
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.
The sympy.interactive.ipythonprinting extension is already loaded. To reload it, use:
  %reload_ext sympy.interactive.ipythonprinting

首先载入我们需要的库,由于SymPy的一个BUG,为了让cse()能正常工作,需要运行下面两行代码:

oo.__class__.q = oo
oo.__class__.p = 1
In [45]:
import numpy as np
import pylab as pl
from sympy import *
oo.__class__.q = oo
oo.__class__.p = 1

椭圆可以用如下的参数方程表示,其中$a$为长轴的长度,$b$为短轴长度,$(x_c, y_c)$为椭圆的中心坐标,$\theta$为椭圆长轴和X轴的夹角,$t$为参变量,$(x(t),y(t))$为椭圆上各点的坐标:

$x(t)=x_c + a\,\cos t\,\cos \theta - b\,\sin t\,\sin\theta$

$y(t)=y_c + a\,\cos t\,\sin \theta + b\,\sin t\,\cos\theta$

下面是用参数方程计算椭圆上各点,并绘图的函数:

In [3]:
def ellipse(e, t):
    from numpy import cos, sin
    a, b, x_c, y_c, theta = e
    ct, st, cth, sth = cos(t), sin(t), cos(theta), sin(theta)
    x = a*ct*cth - b*st*sth + x_c
    y = a*sth*ct + b*st*cth + y_c
    return x, y

def plot_ellipse(e):
    t = np.linspace(0, np.pi*2, 100)
    x, y = ellipse(e, t)
    pl.plot(x, y)    
    
e1 = (5.0, 3.0, 1.0, 2.0, 0.3)
plot_ellipse(e1)
x, y = ellipse(e1, [0, np.pi])
pl.plot(x, y, "--")
pl.plot(x.mean(), y.mean(), "o")
pl.axis("equal");

标准椭圆的隐函数方程为:

$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$

为了计算任意椭圆的的隐函数方程,需要进行旋转和平移两种坐标变换。假设$(x', y')$为任意椭圆上的一点,它和标准椭圆上的点$(x, y)$之间有如下关系:

$ \begin{bmatrix} x' \\ y' \\ \end{bmatrix} = \begin{bmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \\ \end{bmatrix}\begin{bmatrix} x \\ y \\ \end{bmatrix} + \begin{bmatrix} x_c \\ y_c \\ \end{bmatrix}$

其中$\begin{bmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \\ \end{bmatrix}$为旋转矩阵,它将点$(x, y)$绕坐标原点旋转角度$\theta$。将上面的公式进行变换,得到:

$ \begin{bmatrix} x \\ y \\ \end{bmatrix} = \begin{bmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \\ \end{bmatrix}\begin{bmatrix} x'-x_c \\ y'-y_c \\ \end{bmatrix}$

将$(x, y)$带入到标准隐函数方程就可以得到$(x', y')$的隐函数方程。下面我们用SymPy来完成这个计算。

In [4]:
x, y, theta, xc, yc, a, b = symbols("x y theta x_c y_c a b", real=True)
M = Matrix([[cos(theta),  sin(theta)],
            [-sin(theta), cos(theta)]])
P = Matrix([[x-xc],[y-yc]])
P2 = M * P
eq = (x/a)**2 + (y/b)**2 - 1
eq2 = eq.subs({x:P2[0], y:P2[1]}, simultaneous=True)
print eq2
-1 + (-(x - x_c)*sin(theta) + (y - y_c)*cos(theta))**2/b**2 + ((x - x_c)*cos(theta) + (y - y_c)*sin(theta))**2/a**2

隐函数方程的曲线可以通过等值线函数contour()绘制:

In [9]:
def implicit_ellipse(e, x, y):
    from numpy import cos, sin
    a, b, x_c, y_c, theta = e    
    return (-1 + (-(x - x_c)*sin(theta) + (y - y_c)*cos(theta))**2/b**2 + 
            ((x - x_c)*cos(theta) + (y - y_c)*sin(theta))**2/a**2)

yg, xg = np.mgrid[-2:6:100j, -5:7:100j]
zg = implicit_ellipse(e1, xg, yg)

pl.contour(xg, yg, zg, levels=[0]);

两个任意椭圆的交点可以通过如下步骤进行计算:

  1. 将椭圆1用隐函数方程表示。由于椭圆是一个二阶曲线,它的隐函数方程的一般形式为:$A x^2 + B x y + C y^2 + D x + E y + F = 0$

  2. 将椭圆2的参数方程形式:$x = f(t), y = g(t)$代入到上面的隐函数方程中,得到一个关于$sin(t), cos(t)$的二阶方程。

  3. 将上述方程中的$\cos{t}$替换为$c_t$,$\sin{t}$替换为$\sqrt{1-c_t^2}$,并通过一些变换将其转换成一个关于$c_t$的一元四次方程。

  4. 求解上述方程,对$c_t$在-1到1之间的实数解,求反余弦:$t = \pm \arccos{c_t}$。

  5. 将$t$的值代入到参数方程中得到椭圆2上的一点$(x_t, y_t)$,并验证$(x_t, y_t)$是否满足椭圆1的隐函数方程。满足的点就为交点。

由上面的步骤可以看出,其中涉及到许多符号表达式变换。通过SymPy的符号运算功能,可以帮助我们简化程序的编写。

首先让SymPy帮我们计算隐函数方程的一般形式中的那些系数:$A, B, C, D, E, F$. 通过as_poly()方法可以将一个表达式转换成多项式对象,然后通过多项式对象的coeffs()方法可以获得其各项的系数。多元多项式的系数排列可以有很多种方式,这里使用"grlex"排列,它是按照阶数、未知数的字母顺序排列,其结果和步骤1中的一般形式的顺序相同。

In [31]:
p = eq2.as_poly(x, y)
for coef in p.coeffs(order="grlex"):
    display(coef)
$$\frac{a^{2} \sin^{2}{\left (\theta \right )} + b^{2} \cos^{2}{\left (\theta \right )}}{a^{2} b^{2}}$$
$$\frac{- 2 a^{2} \sin{\left (\theta \right )} \cos{\left (\theta \right )} + 2 b^{2} \sin{\left (\theta \right )} \cos{\left (\theta \right )}}{a^{2} b^{2}}$$
$$\frac{a^{2} \cos^{2}{\left (\theta \right )} + b^{2} \sin^{2}{\left (\theta \right )}}{a^{2} b^{2}}$$
$$\frac{- 2 a^{2} x_{c} \sin^{2}{\left (\theta \right )} + 2 a^{2} y_{c} \sin{\left (\theta \right )} \cos{\left (\theta \right )} - 2 b^{2} x_{c} \cos^{2}{\left (\theta \right )} - 2 b^{2} y_{c} \sin{\left (\theta \right )} \cos{\left (\theta \right )}}{a^{2} b^{2}}$$
$$\frac{2 a^{2} x_{c} \sin{\left (\theta \right )} \cos{\left (\theta \right )} - 2 a^{2} y_{c} \cos^{2}{\left (\theta \right )} - 2 b^{2} x_{c} \sin{\left (\theta \right )} \cos{\left (\theta \right )} - 2 b^{2} y_{c} \sin^{2}{\left (\theta \right )}}{a^{2} b^{2}}$$
$$\frac{- a^{2} b^{2} + a^{2} x_{c}^{2} \sin^{2}{\left (\theta \right )} - 2 a^{2} x_{c} y_{c} \sin{\left (\theta \right )} \cos{\left (\theta \right )} + a^{2} y_{c}^{2} \cos^{2}{\left (\theta \right )} + b^{2} x_{c}^{2} \cos^{2}{\left (\theta \right )} + 2 b^{2} x_{c} y_{c} \sin{\left (\theta \right )} \cos{\left (\theta \right )} + b^{2} y_{c}^{2} \sin^{2}{\left (\theta \right )}}{a^{2} b^{2}}$$

我们可以直接将上述结果复制到一个Python函数中,不过仔细观察各项系数,可以发现它们有许多相同的子项,这样会造成许多不必要的重复计算。我们可以使用cse()将一个复杂的表达式分成多步进行计算。下面先看一个简单的例子。

对于为了计算P2中的两个表达式,可以先计算四个临时变量$x_0, x_1, x_2, x_3$,然后P2的就可以简化为$x_0 x_3 + x_1 x_2, -x_0 x_1 + x_2 x_3$。这样就避免了$x-x_c$等子项的重复运算。

cse(expression)expression中的公共子项提取出来用临时变量表示,返回一个二元元组,其中第一个元素为一个列表,列表中的每个元素对应一个临时变量的表达式;元组的第二个元素为采用各个临时变量比表示的原表达式。

In [33]:
display(P2)
cse(P2)
Out[33]:
$$\left[\begin{smallmatrix}\left(x - x_{c}\right) \cos{\left (\theta \right )} + \left(y - y_{c}\right) \sin{\left (\theta \right )}\\- \left(x - x_{c}\right) \sin{\left (\theta \right )} + \left(y - y_{c}\right) \cos{\left (\theta \right )}\end{smallmatrix}\right]$$
$$\begin{pmatrix}\begin{bmatrix}\begin{pmatrix}x_{0}, & x - x_{c}\end{pmatrix}, & \begin{pmatrix}x_{1}, & \sin{\left (\theta \right )}\end{pmatrix}, & \begin{pmatrix}x_{2}, & y - y_{c}\end{pmatrix}, & \begin{pmatrix}x_{3}, & \cos{\left (\theta \right )}\end{pmatrix}\end{bmatrix}, & \begin{bmatrix}\left[\begin{smallmatrix}x_{0} x_{3} + x_{1} x_{2}\\- x_{0} x_{1} + x_{2} x_{3}\end{smallmatrix}\right]\end{bmatrix}\end{pmatrix}$$

下面的两个函数可以将cse()得到结果转换成一个Python函数,用于数值运算。这段程序要解决几个问题:

  1. 将整数用浮点数表示:由于Python的整数除法,为了让诸如1/3的算式能得到正确的结果,需要将它表示为1.0/3.0

  2. 幂为分数,底数为负数:当幂运算符的底数为负数、幂为分数时,Python会抛出ValueError: negative number cannot be raised to a fractional power异常,这时需要采用复数进行运算。因此凡是遇到幂为分数的项,都将底数用complex()转换为复数。

In [37]:
print complex(-0.3)**0.5
print (-0.3)**0.5
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-37-a22dcb2942d6> in <module>()
      1 print complex(-0.3)**0.5
----> 2 print (-0.3)**0.5

ValueError: negative number cannot be raised to a fractional power
(3.35383338435e-17+0.547722557505j)
In [39]:
from sympy.printing import StrPrinter

class FloatExprPrinter(StrPrinter):
    def _print_Integer(self, expr):
        return str(expr)+".0"
    
    def _print_Rational(self, expr):
        return '%s.0/%s' % (expr.p, expr.q)
    
    def _print_Pow(self, expr):
        exp = expr.exp
        if exp.is_Rational and exp.q != 1 and exp is not S.Half:
            return 'complex(%s)**(%s)' % (self._print(expr.base), self._print(exp))
        else:
            return super(FloatExprPrinter, self)._print_Pow(expr)

cse_funcs = {}
def cse2func(funcname, precodes, seq, printer_class=FloatExprPrinter):
    import textwrap
    printer = printer_class()
    codes = ["def %s:" % funcname]
    if isinstance(precodes, basestring):
        precodes = [precodes]
    for line in precodes:
        codes.append("    %s" % line)
    for variable, value in seq[0]:
        codes.append("    %s = %s" % (variable, printer._print(value)))
    returns = "    return (%s)" % ", ".join([printer._print(value) for value in seq[1]])
    codes.append("\n".join(textwrap.wrap(returns, 80)))
    code = "\n".join(codes)
    cse_funcs[funcname] = code
    return code

下面我们用cse2func()得到计算椭圆的隐函数方程系数的函数ellipse_equation():

In [41]:
seq = cse(p.coeffs(order="grlex"))
exec cse2func("ellipse_equation(a, b, x_c, y_c, theta)", "from math import sin, cos", seq)
print cse_funcs["ellipse_equation(a, b, x_c, y_c, theta)"]
def ellipse_equation(a, b, x_c, y_c, theta):
    from math import sin, cos
    x0 = a**2.0
    x1 = b**2.0
    x2 = sin(theta)
    x3 = cos(theta)
    x4 = x_c**2.0
    x5 = y_c**2.0
    x6 = -2.0
    x7 = x3**(-x6)
    x8 = x2**(-x6)
    x9 = x0*y_c
    x10 = -x0
    x11 = x1*y_c
    x12 = x2*x3
    x13 = x0*x7
    x14 = x12*x_c
    x15 = 1.0/(x0*x1)
    x16 = x1*x8
    x17 = x0*x8
    x18 = x1*x7
    return (x15*(x17 + x18), x12*x15*x6*(x0 - x1), x15*(x13 + x16),
x15*x6*(x11*x12 - x12*x9 + x17*x_c + x18*x_c), x15*x6*(x1*x14 + x10*x14 + x11*x8
+ x7*x9), x15*(x1*x10 - x11*x14*x6 + x13*x5 + x14*x6*x9 + x16*x5 + x17*x4 +
x18*x4))

下面是前面的椭圆e1的隐函数方程中各项的系数:

In [43]:
print ellipse_equation(*e1)
(0.046210289247655884, -0.04015235366364695, 0.10490082186345523, -0.01211587116801786, -0.37945093379017397, -0.6144911306258171)

首先参数方程:

In [47]:
theta, xc, yc, a, b, t = symbols("theta x_c y_c a b t", real=True)
x = a*cos(t)
y = b*sin(t)
M = Matrix([[cos(theta), -sin(theta)],[sin(theta), cos(theta)]])
P = Matrix([[x], [y]])
xt, yt = M*P +Matrix([[xc], [yc]])
(xt, yt)
Out[47]:
$$\begin{pmatrix}a \cos{\left (t \right )} \cos{\left (\theta \right )} - b \sin{\left (t \right )} \sin{\left (\theta \right )} + x_{c}, & a \sin{\left (\theta \right )} \cos{\left (t \right )} + b \sin{\left (t \right )} \cos{\left (\theta \right )} + y_{c}\end{pmatrix}$$

将参数方程代入到隐函数方程中:

In [50]:
A, B, C, D, E, F, x, y = symbols("A B C D E F x y")
eq = A*x**2 + B*x*y + C*y**2 + D*x + E*y + F
eq = expand(eq.subs({x:xt, y:yt}))

将$\cos{t}$和$\sin{t}$替换为$t_c$和$\sqrt{1-t_c^2}$:

In [52]:
tc = symbols("t_c")
sqrt_term = sqrt(1-tc**2)
eq1 = eq.subs({cos(t):tc, sin(t):sqrt_term})

leq为含有根号的项,req为不含根号的项。两边平方之后应该相等:

In [53]:
leq = sqrt_term * eq1.coeff(sqrt_term)
req = eq1 - expand(leq)
eq_square = expand(leq**2 - req**2)

eq_square转换为$t_c$的多项式,并对其系数进行cse()进行处理,最后将其转换为Python函数:

In [55]:
p = eq_square.as_poly(tc)
seq = cse(p.coeffs())
exec cse2func(
    "ellipses_intersection_equation(A, B, C, D, E, F, a, b, x_c, y_c, theta)", 
    "from math import sin, cos", seq)

下面通过cse()得到一元四次方程$a x^4 + b x^3 + c x^2 + d x + e = 0$的解。注意这里需要判断最高项的系数是否接近0,如果接近0的话,就变为求解低一次的方程了。一元N次方程的解可以通过roots(Poly(expression, x))得到。

In [56]:
a,b,c,d,e,x = symbols("a,b,c,d,e,x", real=True)
r = roots(Poly(a*x**4 + b*x**3 + c*x**2 + d*x + e), x).keys()
seq = cse(r)
exec cse2func("roots4(a,b,c,d,e)", ["from cmath import sqrt", "if abs(a)<1e-10: return roots3(b,c,d,e)"], seq)

r = roots(Poly(b*x**3 + c*x**2 + d*x + e), x).keys()
seq = cse(r)
exec cse2func("roots3(b,c,d,e)", ["from cmath import sqrt", "if abs(b)<1e-10: return roots2(c,d,e)"], seq)

r = roots(Poly(c*x**2 + d*x + e), x).keys()
seq = cse(r)
exec cse2func("roots2(c,d,e)", ["from cmath import sqrt", "if abs(c)<1e-10: return roots1(d,e)"], seq)

r = roots(Poly(d*x + e), x).keys()
seq = cse(r)
exec cse2func("roots1(d,e)", ["from cmath import sqrt", "if abs(d)<1e-10: return []"], seq)

使用前面自动生成的各个函数,我们可以很方便地编写计算椭圆交点的函数。

ellipse(e, t):计算参数方程系数为e的椭圆在t处的坐标,如果t为单个数值,则采用math下的函数,否则采用numpy下的函数。

ellipse_implicit(e, x, y):计算隐函数系数为e的椭圆在(x, y)处的值,若值为0则表示(x, y)在此椭圆之上。

ellipse_intersections(e1, e2):计算两个用参数方程系数表示的椭圆的交点坐标。

In [57]:
def ellipse(e, t):
    if np.isscalar(t):
        from math import cos, sin
    else:
        from numpy import cos, sin
    a, b, x_c, y_c, theta = e
    ct, st, cth, sth = cos(t), sin(t), cos(theta), sin(theta)
    x = a*ct*cth - b*st*sth + x_c
    y = a*sth*ct + b*st*cth + y_c
    return x, y

def ellipse_implicit(e, x, y):
    a, b, c, d, e, f = e
    return a*x**2 + b*x*y + c*y**2 + d*x + e*y + f

def ellipse_intersections(e1, e2, eps=1.0e-6):
    from math import acos, sqrt
    e1imp = ellipse_equation(*e1)
    p =  ellipses_intersection_equation(*(tuple(e1imp) + e2))
    roots = [r.real for r in roots4(*p) if abs(r.imag) < eps]
    
    points = []
    for root in roots:
        t = acos(float(root))
        for t2 in [t, -t]:
            xp, yp = ellipse(e2, t2)
            if abs(ellipse_implicit(e1imp, xp, yp)) < eps:
                if not points:
                    points.append((xp, yp))
                else:
                    mindist = min([sqrt((po[0]-xp)**2+(po[1]-yp)**2) for po in points])
                    if mindist > eps:
                        points.append((xp, yp))
    return points
In [58]:
e1 = (5.0, 3.0, 1.0, 2.0, 0.2)
e2 = (7.0, 4.0, 0, 0, -0.3)
e3 = (5.0, 3.0, 1.0, 2.0, 0.2 + np.pi/2)

def ellipse_plot(e1, e2, draw_e1=True, draw_e2=True):
    t = np.linspace(0, np.pi*2, 100)
    x1, y1 = ellipse(e1, t)
    x2, y2 = ellipse(e2, t)
    if draw_e1:
        pl.plot(x1, y1)
    if draw_e2:
        pl.plot(x2, y2)
    
    points = ellipse_intersections(e1, e2)
    xp = [p[0] for p in points]
    yp = [p[1] for p in points]
    pl.plot(xp, yp, "o");
    
ellipse_plot(e1, e2)
ellipse_plot(e1, e3, False, True)
ellipse_plot(e2, e3, False, False)
In [59]:
e1 = (5.0, 3.0, 1.0, 2.0, 0.2)
e2 = (5.0, 3.0, -1.0, 1.2, 0.2)
ellipse_plot(e1, e2)
In [60]:
ellipse_plot((5.0, 3.0, 1.0, 2.0, 0.0), (5.0, 3.0, 0.0, 0.0, 0.0))
In [61]:
ellipse_plot((5.0, 3.0, 1.0, 2.0, 0.0), (5.0, 3.0, 1.0, 2.0, 0.01))
In [62]:
ellipse_plot((5.0, 3.0, 1.0, 2.0, 0.0), (5.0, 3.0, 1.0, 2.1, 0.0))

loading...