如何找到两个圆的公切线？

最新推荐文章于 2021-02-27 00:50:29 发布

栉风沐雨1314

最新推荐文章于 2021-02-27 00:50:29 发布

阅读量8.9k

点赞数 3

分类专栏： python 文章标签： python 算法圆公切线 plot

本文链接：https://blog.csdn.net/gt11799/article/details/39205629

版权

python 专栏收录该内容

35 篇文章 0 订阅

订阅专栏

大家都提供了道，我来提供术。我最近在自学python，题主如果想学python，可以私信微信号给我，我们一起努力啊。
我采用的是 @vczh 的算法，期间程序有问题，我手算了一条内公切线的斜率，用的是 @程德华的方法。（我有多想念老师们出题的数值！老师们节日快乐。）
先说结论：python可以算，但是很麻烦，而且我的程序最终还是有问题。推荐使用matlab。
——————————分割线——————————————
点到直线的距离：
$d = \frac{|Ax_0+By_0+C|}{\sqrt{ A^2+B^2}}$ (还复习了Latex，真是怀念大学时光)
把分数上下除以B就可以得到k，b的等式。
$d=\frac{|kx_0-y_0+b|}{\sqrt{k^2+1}}$
圆以圆心加半径的方式传入：
$circle_1 = ((x_1, y_1), r_1),circle_2 = ((x_2,y_2), r_2)$
圆心到直线距离等于半径，我们可以得到两个等式。
$\frac{|kx_1-y_1+b|}{\sqrt{k^2+1}}=r_1\\\frac{|kx_2-y_2+b|}{\sqrt{k^2+1}}=r_2$
我们得到了一个二元二次方程组，最多的时候会有四组解。
在matlab里有solve函数，可以直接求出多元多次方程组的解，这样整个过程就结束了。
————————matlab党止步————————————————
python中有科学计算三驾马车，numpy，scipy和maplotlib。其中scipy的optimize中有fsolve函数，可以解多元多次方程组，但是，他需要给出初值，并且返回该初值附近的一组解。
是的，只有一组解。我所想到的办法就是从头开始赋初值，尽量别错过每一种可能性。这样程序就变得臃肿，而且效率低下。fsolve本来就是计算量巨大的函数，还要计算那么多次。（我多想有人过来给我一巴掌，笨蛋，＊＊模块都没听过？）
我既然是来提供术的，我会把详细的过程写下来。
核心模块，解二元二次方程组。

from scipy.optimize import fsolve
k,b = fsolve(funcForSolve, [1,1])

其中[1,1]是初值，这样就可以返回[1,1]附近的方程组的解。当然我们得要四组解（如果有的话）

from numpy import *
klist = linspace(-50, 50, 200)
blist = linspace(-100, 100, 500)
for k in klist:
    for b in blist:
        (solved_k, solved_b) = fsolve(funcForSolve, [k,b])

如何去除重复解？

if abs(k2-k1) > PRECISION and abs(b1-b2) > PRECISION:
      result.append((solved_k, solved_b))

PRICISION是个常数，我定义成了0.01。
可以想象整个程序的复杂度。但是根据我的观察，发现有几个特点可以利用：
1.四条公切线的斜率肯定不同，因此判定时只需要判定k即可。（半径不同的时候）
2.在其中一个循环中找到一组k，b，就可以结束内循环，跳到下一个k。
3.用很小的初值找到一组解，可以认为其他的解比刚求的解大。
4.当取得四组解时，整个过程就可以结束了。

def solve_func():
    result = []
    klist = linspace(-50, 50, 200)
    blist = linspace(-100, 100, 500)
    (solved_k, solved_b) = (-float('Inf'), -float('Inf'))
    for k in klist:
        for b in blist:
            if k > solved_k:
                (solved_k, solved_b) = fsolve(funcForSolve, [k,b])
                try:
                    delta_k = min([abs(solved_k - item[0]) \
                        for item in result])                        
                except(IndexError, ValueError):
                    delta_k = 10 # for first solves
                
                if delta_k > PRECISION:
                    result.append((solved_k, solved_b))
                    break
        if len(result) >= 4:
            break
    return result

注：1.斜率和截距超出范围外就无能为力了。
2.如果两个圆是分离的，提高斜率和截距的上下限不会影响太多性能。
其中try是为了得到第一组值的时候，result是空列表，遍历会出错。
定义一个函数funcForSolve，输入两个圆，以列表的形式返回方程组即可。
——————————回答完毕————————————————
我感觉题主应该是python新手，以下是更详细的过程。
思考一个程序的结构，我的习惯是从核心模块开始，往上和往下去思考。
——————
往上：
核心模块需要什么？一个方程组表达式的列表
如何得到这个列表？把圆的坐标带入公式
要是有圆内含怎么办？需要根据圆心距离判断
需要一个计算点到点距离的模块。
往下：
核心模块输出了什么？一个列表，里面k,b组成的tuple
用k,b组成直线的方程式输出
绘图。
绘图需要什么？
圆的坐标，直线的方程，坐标范围。
————————
我喜欢混沌的思维，总是尽可能把程序砸碎，看看有没有可以重用的模块。根据k,b返回直线方程就是可以重用的模块。

其中黑色箭头是程序引用的关系，红色是主函数中的顺序。我一直试图把输入控制在主函数中，后来发现函数的调用变得难以理解。

首先是怎么得到方程组。

def circle_to_func(circle_1, circle_2):
    (x1, y1) = circle_1[0]
    r1 = circle_1[1]
    (x2, y2) = circle_2[0]
    r2 = circle_2[1]
    
    dist = point_to_point(circle_1[0], circle_2[0])
    assert dist > circle_1[1] and dist > circle_2[1], \
    "Containing circle has no common tangent."
    
    return '(k*%f - %f + b) ** 2 - %f**2 * (k**2 + 1)' %(x1,y1,r1),\
    '(k*%f - %f + b) ** 2 - %f**2 * (k**2 + 1)' %(x2,y2,r2)
    
def funcForSolve(x):
    'convert the string of functions to the polynomial'
    k,b = x.tolist()
    func = circle_to_func(CIRCLE_1, CIRCLE_2)
    return [eval(func[0]), eval(func[1])]

python中不支持表达式，只能把要计算的等式当作字符串存储。circle_to_func就是把圆的坐标带入公式中，而funcForSolve则是要把字符串转换成方程组。而给fsolve的参数函数最好是只有一个输入，方便输入初值。我把内含圆的判断也放到了这里。
然后就是得到k,b后，转换成方程式输出。

def argu_to_line(arguments):
    result = []
    for (k,b) in arguments:
        result.append("%.4f*x%+.4f" %(k,b))
    return result

注意这里的方程式是要用来计算的，因此表达式得是符合python运算法则的。
然后画圆

fig = plt.figure(dpi=120)
ax = fig.add_subplot(111)

circle_plot_1 = plt.Circle(circle_1[0], circle_1[1], fill=False)
circle_plot_2 = plt.Circle(circle_2[0], circle_2[1], fill=False)
ax.add_artist(circle_plot_1)
ax.add_artist(circle_plot_2)

plt是pyplot的简称，是matplotlib中的模块
画线

limit = axes_limits(circle_1, circle_2)
x = linspace(limit[0][0], limit[0][1], 100)
for func in line_func:
    plt.plot(x, eval(func), 'r')

我根据圆心和半径建立了一个函数，返回坐标的范围。也就是上面的axes_limits()
然后就是简单的设置，显示或者保存图像。

ax.set_xlim(limit[0])
ax.set_ylim(limit[1])
ax.set_xlabel('x')
ax.set_ylabel('y')
    
plt.show()
fig.savefig('common_tangent.png', transparent=True, dpi=300)

当输入是：

PRECISION = 0.01
CIRCLE_1 = ((10,10),1)
CIRCLE_2 = ((5,6),2)

输出是

y=0.1893*x+7.0888
y=1.0969*x-2.4528
y=0.5698*x+5.4528
y=2.3107*x-10.5888

以及：

图比较简陋，实在不好意思呐。

源代码：gt11799/Common_tangent · GitHub

原载于知乎

——————————————

github主页：https://github.com/gt11799

E-mail：gting405@163.com

栉风沐雨1314

关注

3
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
如何找到两个圆的公切线？

大家都提供了道，我来提供术。我最近在自学python，题主如果想学python，可以私信微信号给我，我们一起努力啊。我采用的是 @vczh的算法，期间程序有问题，我手算了一条内公切线的斜率，用的是 @程德华的方法。（我有多想念老师们出题的数值！老师们节日快乐。）先说结论：python可以算，但是很麻烦，而且我的程序最终还是有问题。推荐使用matlab。——————————分割线————
复制链接

扫一扫