获取由两点和切线定义的圆，Lua

Question

我编写了一个代码来获取给定两点和切线的圆（或两个圆）。

我正在制定算法，用给定的两点和该圆的切线林德来制作一个圆。

function tangentCircle (cx, cy, cr, px, py) -- circle, point
    local dx, dy = px-cx, py-cy
    local dist = math.sqrt(dx*dx+dy*dy)
    if dist <= cr then return end
    local a, b = math.atan2 (dy, dx), math.acos (cr/dist)   
    local x1 = cx + cr * math.cos (a-b)
    local y1 = cy + cr * math.sin (a-b)
    local x2 = cx + cr * math.cos (a+b)
    local y2 = cy + cr * math.sin (a+b)
    return x1, y1, x2, y2
end

function linesIntersect(x1, y1, x2, y2, x3, y3, x4, y4)
    local dx1, dy1 = x2 - x1, y2 - y1
    local dx2, dy2 = x4 - x3, y4 - y3
    local det = dx2 * dy1 - dx1 * dy2
    if det == 0 then return end
    local c1 = x1 * y2 - y1 * x2
    local c2 = x3 * y4 - y3 * x4
    local x = (c1 * dx2 - c2 * dx1) / det
    local y = (c1 * dy2 - c2 * dy1) / det
    return x, y
end

function circleWithTwoPoints (x1, y1, x2, y2)
    local cx = (x1 + x2) / 2
    local cy = (y1 + y2) / 2
    local cr = math.sqrt((cx - x1)^2 + (cy - y1)^2)
    return cx, cy, cr
end


function intersectCircleCenterLine(cx, cy, cr, tx1, ty1, tx2, ty2)
    local angle = math.atan2 (ty2-ty1, tx2-tx1)
    local ex1, ey1 = cx + cr * math.cos (angle), cy + cr * math.sin (angle)
    local ex2, ey2 = cx + cr * math.cos (angle+math.pi), cy + cr * math.sin (angle+math.pi)
    return ex1, ey1, ex2, ey2
end

function nearestPointOnLine(tx1, ty1, tx2, ty2, cx, cy)
    local dx = tx2 - tx1
    local dy = ty2 - ty1
    local t = ((cx - tx1) * dx + (cy - ty1) * dy) / (dx * dx + dy * dy)
    local nx = tx1 + t * dx
    local ny = ty1 + t * dy
    return nx, ny
end

function circleFromThreePoints(x1, y1, x2, y2, x3, y3)
    local function distance(x1, y1, x2, y2)
        return math.sqrt((x1-x2)^2 + (y1-y2)^2)
    end

    local a = distance(x2, y2, x3, y3)
    local b = distance(x1, y1, x3, y3)
    local c = distance(x1, y1, x2, y2)
    local s = (a + b + c) / 2
    local gr = (a * b * c) / (4 * math.sqrt(s * (s - a) * (s - b) * (s - c)))
    local A = x1*(y2-y3) - y1*(x2-x3) + x2*y3 - x3*y2
    local B = (x1*x1 + y1*y1)*(y3-y2) + (x2*x2 + y2*y2)*(y1-y3) + (x3*x3 + y3*y3)*(y2-y1)
    local C = (x1*x1 + y1*y1)*(x2-x3) + (x2*x2 + y2*y2)*(x3-x1) + (x3*x3 + y3*y3)*(x1-x2)
    local D = 2*((x1-x2)*(y3-y2) - (y1-y2)*(x3-x2))

    local gx = B / D
    local gy = C / D
    return gx, gy, gr
end


function circlesFromTwoPointsAndTangent(ax, ay, bx, by, tx1, ty1, tx2, ty2)
    local px, py = linesIntersect(ax, ay, bx, by, tx1, ty1, tx2, ty2)
    print ('px, py', px, py) -- must be -10, 0, ok

    if px then
        -- two solutions
        local cx, cy, cr = circleWithTwoPoints (ax, ay, bx, by)
        print ('cx, cy, cr', cx, cy, cr) -- must be 65, 75, 49.497475, ok

        local dx1, dy1, dx2, dy2 = tangentCircle (cx, cy, cr, px, py)
        print ('dx1, dy1', dx1, dy1) -- must be 17.7115, 89.6218, ok
        print ('dx2, dy2', dx2, dy2) -- must be 79.6218, 27.7115, ok


--  radius from P to D or from P to E
        local pr = math.sqrt ((dx1-px)^2+(dy1-py)^2)
        print ('pr', pr) -- must be 93.808315

        local ex1, ey1, ex2, ey2 = intersectCircleCenterLine(px, py, pr, tx1, ty1, tx2, ty2)
        local gx1, gy1, gr1 = circleFromThreePoints(ax, ay, bx, by, ex1, ey1)
        print ('gx1, gy1, gr1', gx1, gy1, gr1)
        -- must be 83.8083, 56.1917, 56.191685, ok

        local gx2, gy2, gr2 = circleFromThreePoints(ax, ay, bx, by, ex2, ey2)
        print ('gx2, gy2, gr2', gx2, gy2, gr2)
        -- must be -103.80831519647 243.80831519647 243.80831519647 ok

        return gx1, gy1, gr1, gx2, gy2, gr2
    else -- parallel
        local cx, cy, cr = circleWithTwoPoints (ax, ay, bx, by)

        local ex, ey = nearestPointOnLine(tx1, ty1, tx2, ty2, cx, cy)

        local gx, gy, gr = circleFromThreePoints(ax, ay, bx, by, ex, ey)

        return gx, gy, gr
    end
end

示例：

local x1, y1 = 30, 40
local x2, y2 = 100, 110
local tx1, ty1 = -100, 0
local tx2, ty2 = 100, 0 -- horizontal line

local gx1, gy1, gr1, gx2, gy2, gr2 = circlesFromTwoPointsAndTangent(x1, y1, x2, y2, tx1, ty1, tx2, ty2)

print ('circle 1', gx1, gy1, gr1) -- must be 83.8083    56.1917 56.191685
print ('circle 2', gx2, gy2, gr2) -- must be -103.808   243.808 243.808315

如何查看是否无解？或者输入有限制，例如垂直或水平线有问题？它有什么限制？

Answer 1

我已经用Python实现了解决方案。

它是如何工作的：为了简化公式，我执行坐标移动将第一个点放置在原点，旋转将第二个点放置在 OX 轴上，然后向左移动一半的距离。现在圆心将位于OY轴上，它的坐标是

(0,y)

。

我对线点进行变换，然后将线方程转换为正规方程形式(

px+qy+r=0

)。

然后求解二次方程：半径的平方应等于到直线的距离的平方（所选直线方程形式的简单公式）。

二次方程可能没有解（例如，点之间的线）、一个解（例如，穿过第二个点并垂直于第一第二向量的线）或两个解。

然后将坐标变换回来。为我的图片和您的数据提供了示例。

import math
def circle2ptstangenettoline(ax, ay, bx, by, tx1, ty1, tx2, ty2):
    dist = math.hypot(ax-bx, ay-by)
    h = dist/2
    ang = math.atan2(by-ay, bx-ax)
    ux1 = tx1 - ax
    uy1 = ty1 - ay
    ux2 = tx2 - ax
    uy2 = ty2 - ay
    vx1 = ux1 * math.cos(-ang) - uy1 * math.sin(-ang) - h
    vy1 = ux1 * math.sin(-ang) + uy1 * math.cos(-ang)
    vx2 = ux2 * math.cos(-ang) - uy2 * math.sin(-ang) - h
    vy2 = ux2 * math.sin(-ang) + uy2 * math.cos(-ang)
    vd = math.hypot(vx2-vx1, vy2-vy1)
    p = (vy2-vy1) / vd
    q = (vx1-vx2) / vd
    r = (vx2*vy1-vx1*vy2) / vd
    a = q*q-1
    b = 2*q*r
    c = r*r - h*h
    if a==0:
        if b==0:
            return None
        else:
            y = -c/b
            cr = math.hypot(y, h)
            cx = h * math.cos(ang) - y * math.sin(ang) + ax
            cy = h * math.sin(ang) + y * math.cos(ang) + ay
            return (cx, cy, cr)
    Dis = b*b-4*a*c
    if Dis < 0:
        return None
    if Dis == 0:
        y = -0.5*b/a
        cr = math.hypot(y, h)
        cx = h * math.cos(ang) - y * math.sin(ang) + ax
        cy = h * math.sin(ang) + y * math.cos(ang) + ay
        return (cx, cy, cr)
    y1 = 0.5*(-b - math.sqrt(Dis))/a
    cr1 = math.hypot(y1, h)
    cx1 = h * math.cos(ang) - y1 * math.sin(ang) + ax
    cy1 = h * math.sin(ang) + y1 * math.cos(ang) + ay

    y2 = 0.5*(-b + math.sqrt(Dis))/a
    cr2 = math.hypot(y2, h)
    cx2 = h * math.cos(ang) - y2 * math.sin(ang) + ax
    cy2 = h * math.sin(ang) + y2 * math.cos(ang) + ay
    return ((cx1,cy1,cr1), (cx2,cy2,cr2))


print(circle2ptstangenettoline(-3, 0, 3, 0, 5, 4, -7, -11))
print(circle2ptstangenettoline(-3, 0, 3, 0, 5, 4, 7, -11))
print(circle2ptstangenettoline(30, 40, 100, 110, -100, 0, 100, 0))

None
((0.0, 3.9528611794604025, 4.962369545296388), (0.0, -5.428416735015959, 6.202234133681292))
((-103.8083151964687, 243.8083151964687, 243.80831519646873), (83.80831519646861, 56.1916848035314, 56.191684803531416))