Get the circle, that is defined by two points and tangent line, Lua

I've made a code to get the circle (or two circles) for given two points and tangential line.

I am making the algorithm to make a circle with given two points and tangential linde for this circle.

Task: For given A and B points and tangent line, defined as T points, build the circle G, that goes through A and B and tangent to the line T.

This: PPL: Apollonius' Problem with Two Points and a Line

Solution:

Find the intersection point P for AB and T1T2.
Create circle with center in point between And B with radius CB.
Find tangent line from P to circle C, the tangent point is D.
Get radius PD. The radius belongs to circle P.
Find two points E1 and E2 by crossing line T and circle P.
Get circles G1 and G2 by given three points (A, B, G) as circumcircle.
Both G circles are solution of the task.

Desmos with one simple solution: ox2zww716t

Image:

Help functions:

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

Main function:

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

Example:

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

How to check if it has no solutions? Or the input has limitations, for example vertical or horizontal lines have issues? What limitations it has?

Solution

I've implemented solution in Python.

How it works: to simpify formulas, I perform coordinate shift to place the first point in origin, then rotation to place the second point onto OX axis, then shift left by half of distance. So circle center will lie on OY axis, and it's coordinates are (0,y).

Now I apply transformations to line points, then transform line equation into normal equation form (px+qy+r=0).

Then solve quadratic equation: squared radius should be equal to squared distance to the line (simple formula for chosen line equation form).

Quadratic equation might give no solutions (for example, line between points), one solution (for example, line through second point and perpendicular to first-second vector) or two solutions.

Then transform coordinates back. Example is given for my picture and for your data.

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