Search code examples
mathluageometrytrigonometry

Get third control point quadratic Bezier curve for parabola with given fucus and directrix, Lua


I am trying to make the quadatic Bezier curve that is equivalent to the given parabola, that was defined as focus (fx, fy) and directrix dirY. (algorithm to make Voronoi diagram with sweeping line, the Fortune's algorithm)

The curve is limited by points A (ax, ay) and B (bx, by).

I need to find the third control point C (cx, cy) as here: https://en.m.wikipedia.org/wiki/File:Quadratic_Bezier_parabola_equivalence.svg or https://en.wikipedia.org/wiki/File:Relationship_between_parabola_and_quadratic_Bezier.svg

function getBezierControlPoint_focus_directrix(fx, fy, ax, ay, bx, by)
    -- (x-h)^2=4*p*(y-k)
    if (ay == by) then
        -- exception: horizontal AB
        local k = (fy + dirY) / 2
        local h = fx
        local cx = h
        local cy = ay + 2*(k-ay)
        return cx, cy -- it works perfect
    else
--  [Axis direction](https://en.wikipedia.org/wiki/Parabola#Axis_direction)
        local h = fx
        local k = (fy + dirY) / 2
        local cx = (ax+bx)/2
        local f = (k-dirY)/2
        
        -- vertex
        -- h = -b/(2*a)
        -- k = (4*a*c-b*b)/(4*a)
        
        local a = -1/(dirY-ay)
        local b = -2 * a * h
        local c = h * h * a + k
        print ('h', h, -b/(2*a)) -- ok, same
        print ('k', k, (4*a*c-b*b)/(4*a)) -- ok, same
        
        -- derivative value of parabola in point A (ax):
        local day = a*2*ax + b
        local cy = ay + day * (bx-ax)/2
        
        return cx, cy -- not right, here is problem
    end
end

usage: three control points to draw parabola:

controlPoints = {ax, ay, cx, cy, bx, by}

Parabola as beach line, drawn with bezier-curve

Parabola as beach line, drawn with bezier-curve Update: https://www.desmos.com/calculator/joeugogi8b

Ready function:

function getBezierControlPoint (fx, fy, ax, bx)
    local function f(x) return (x*x-2*fx*x+fx*fx+fy*fy-dirY*dirY) / (2*(fy-dirY)) end
    local function df(x) return (x-fx) / (fy-dirY) end
    if (fy == dirY) then return end -- not parabola
    local ay, by = f(ax), f(bx)
    local ad, dx = df(ax), (bx-ax)/2
    local cx, cy = ax+dx, ay+ad*dx
    return cx, cy
end

Solution

  • If you're working with a focus x=a,y=b and directrix y=c, then we're already axis-aligned, so that saves a bunch of work.

    A focus/directrix definition of the parabola is based on distance equivalences, so let's see which distances we can find for the "tip" x,y of the parabola:

    1. dist((x,y), (a,b)) = sqrt((x-a)² + (y-b)²), and
    2. dist((x,y), (x,c)) = abs(y - c)

    Easy enough, we can solve the equality sqrt((x-a)² + (y-b)²) = abs(y - c) (either by hand or by using something like wolfram alpha) and get:

            x² - 2ax + a² + b² -c²
    f(x) = -----------------------
                   2(b-c)
    

    which in lua is:

    function f(x, a, b, c)
        local n = x*x - 2*a*x + a*a + b*b - c*c
        local d = 2*(b-c)
        return n / d
    end
    

    (and we're ignoring what happens when b==c, because we'll never call this function if that is the case, because if the focus lies on the directrix, that's not a parabola)

    Which means we now have a trivial-to-implement function for getting y=f(x) given any values a, b, and c, as long as b and c are not the same value.

    We then want to find the equivalent Bezier expression for the segment with end points A=(xa,f(xa)) and B=(xb,f(xb)), which means we're "done" already: for any quadratic Bezier curve, the control point lies at the intersection of the tangent line at the start, and the tangent line at the end of the segment.

    Also note that we do not get a choice in what ay and by are: you only get to pick ax and bx, because f(x) defines what the corresponding ay and by must be.

    Now, we (implicitly) know what the tangent lines are, because we know the derivative of f(x) (again, either because we worked that out by hand, or we just told a computer to do it for us):

            x-a
    f'(x) = ---
            b-c
    

    which in lua is:

    function df(x, a, b, c)
        return (x-a) / (b-c)
    end
    

    And so finding that control point is a matter of calculating a standard line/line intersection, because you know A and the tangent at A, and we know B and the tangent at B. Lots of posts on Stackoverflow and pages on the web that explain how to calculate a line/line intersection, but let's use this one:

    function lli(x1, y1, x2, y2, x3, y3, x4, y4)
        -- compact this as much as you feel the need to: lots of repeating values.
        local nx = (x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * y4 - y3 * x4)
        local ny = (x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * y4 - y3 * x4)
        local d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4)
    
        if (d == 0) then
            error("these lines don't intersect")
        end
    
        return { nx / d, ny / d}
    end
    

    And with that, we have everything we need:

    function fdParabolaToBezier(a, b, c, ax, bx)
        if (b == c) then
            error("this is not a parabola")
        end
    
        -- y coordinate for A:
        local ay = f(ax, a, b, c)
    
        -- tangent slope for A:
        local ad = df(ax, a, b, c)
    
        -- y coordinate for B:
        local by = f(bx, a, b, c)
    
        -- tangent slope for B:
        local bd = df(bx, a, b, c)
    
        -- and then the tangent intersection:
        local C = lli(
            ax, ay,          -- A, which by definition lies on A's tangent
            ax + 1, ay + ad, -- The point on A's tangent at x=ax+1
            bx, by,          -- B, which by definition lies on B's tangent
            bx + a, by + bd  -- The point on B's tangent at x=ab+1
        )
    
        return {ax, ay, C[1], C[2], bx, by}
    end
    

    And then an example call for getting the Bezier equivalent of the curve segment on a parabola with focus (2,5) and directrix y=3, starting at x=-5 and ending at x=4:

    local coords = fdParabolaToBezier(2, 5, 3, -5, 4)
    print(table.concat(coords, ', '))
    

    Which yields -5, 16.25, 6.9117647058824, 7.9117647058824, 4, 5.0