as promised, here my answer to your question.

As I can't follow your code completely (I spend most of my time implementing my code :D ) I can only explain to you what the last step is.

So why do we do interpolation here at all? Well the article already describes well that this piece helps creating a smoother geometry. This works, because we take into account how much part of the blob is within a grid block. Looking at the first picture above, take a look into the first block and the second block. The bottom right point has on the left side a scalar value of 2 and on the right side a value of 5. Taking the explanation here into account this must mean that the right block contains more "blob" than the left block which can be seen as the line in the right block is nearly cutting the block in half (see picture below).

Now how does this interpolation work? We simply try to calculate a weight which drags the point more up or down, depending on how the scalar values between the two points are. E.g. if you have a value of 0 for a top point and a value of 2 for a bottom point we expect the resulting point on the line between these points to be exactly in the middle. If you now increase the value of the bottom point from 2 to e.g. 10, we expect the resulting point on the line to be much nearer to the bottom point as it "weights" more.

More mathematically spoken: Important to note here is that the exLerp function is returning the relative value for the y value of Q. This is always expected to be a value between 0.0 to 1.0. This property is also guranteed as long as you don't mess up with supplying only pairs of values to your exLerp function with one of the values being below 1.0 (our treshold value) and the other value being above 1.0.

If you wonder why this is guranteed, just take the formula you implemented (which is correct btw) and supply for x arbitrary values which are greater than 1.0 and for y values which are less than 1.0. The return value will always be between 0.0 and 1.0. If you violate these properties this doesn't hold anymore because the bottom part of the division becomes smaller than the upper part (thus the value will always be greater than 1).

So as soon as you have a surface transition, you need to call the exLerp function with the scalar value which is inside of the volume, e.g. in this case the scalar value of the point D (speaking of voxels) and the scalar value of the vertex which is outside of the volume (in this case point B).

For your algorithm this means you need to find out at the code position when you call interpolate whether the return value regarding a specific point is meant to replace the x or y value of the point you're returning. I guess one variant could be here to check which of the values is 0.5 (if it would be relative position values) because this is always the value you want to replace if you're coming from the non interpolated variant.

Below you can find some F# code which correctly implements the algorithm. I tried to not optimise that much and made use of more descriptive types so the code better represents what the intention is.

PS: I suspect that you will have questions as I wrote this answer pretty late (and as soon as it's late one makes more faults and explains less good but I had a promise to keep ;) ).

type Position = { x: float; y: float}
    
type Metaball = { radius: float; position: Position }
    
let calculateOccupiedPoints (metaballs: Metaball[]) width height gridDensity =
    Map [|  for x in 0..width - 1 do
                for y in 0.. height - 1 do
                    // The grid density specifies how big a point is, e.g. 1.0 means one block, 0.5 means half a block
                    let xF, yF = float x * gridDensity, float y * gridDensity
                    yield (x, y), metaballs  |> Array.sumBy(fun metaball -> metaball.radius ** 2.0 / ((metaball.position.x - xF) ** 2.0 + (metaball.position.y - yF) ** 2.0))
    |]
    
    
type Line = { Start: Position; End: Position }
    
type Direction =
    | Top
    | Bottom
    | Left
    | Right
    
// The second half of the table is the frist half of the table but only reversed.
let marchingSquaresTable =
    let firstHalf =
        [|
            (false,false,false,false), [||]
            (false,false,false,true), [| Left, Bottom |]
            (false,false,true,false), [| Bottom , Right|]
            (false,false,true,true), [| Left, Right |]
            (false,true,false,false), [| Top, Right |]
            (false,true,false,true), [| Left, Top; Bottom, Right |]
            (false,true,true,false), [| Top, Bottom |]
            (false,true,true,true), [| Left, Top |] |]
    let secondHalf =
        [|
            (true,false,false,false), [| |]
            (true,false,false,true), [| |]
            (true,false,true,false), [| |]
            (true,false,true,true), [| |]
            (true,true,false,false), [| |]
            (true,true,false,true), [| |]
            (true,true,true,false), [| |]
            (true,true,true,true), [| |]
        |] |> Array.mapi(fun i (_case, _) -> (_case, firstHalf.[7 - i] |> (fun (_, dirs: (Direction * Direction) []) -> dirs)))
    Map.ofArray (firstHalf |> Array.append secondHalf)
    
let calculatePointsViaLinearInterpolation (points: Map<int * int, float>) height width dontInterpolate =
    [|
        for x in 0..(width - 4) do
            for y in 0..(height - 4) do
                let edges =
                    [|
                        0, 0 // A
                        1, 0 // B
                        1, 1 // D
                        0, 1 // C
                    |]
                let xF, yF = float x, float y
                let edgePositions = edges |> Array.map(fun (x0, y0) -> (x0 + x, y0 + y))
                let edgeValues = edgePositions |> Array.map(fun key -> points.[key])
                let [| va; vb; vd; vc |] = edgeValues
                let edgeConfiguration =
                    let [| a; b; d; c |] = edgeValues |> Array.map(fun value -> value >= 1.0)
                    marchingSquaresTable.[(a, b, d, c)]
                let inline mkPoint x y = { x = x; y = y }
                let inline tupleToPoint (x, y) = mkPoint x y
                let getPosForDirection dir =
                    match dir with
                    | Top -> mkPoint 0.5 0.0
                    | Bottom -> mkPoint 0.5 1.0
                    | Left -> mkPoint 0.0 0.5
                    | Right -> mkPoint 1.0 0.5
                let pointValuesForDirection = function
                    | Top -> va, vb
                    | Left -> va, vc
                    | Right -> vb, vd
                    | Bottom -> vc, vd
    
                yield!
                    edgeConfiguration
                    |> Array.map(fun (dirA, dirB) ->
                        let applyInterpolation direction =
                            let (x, y) = pointValuesForDirection direction
                            let pos = getPosForDirection direction
                            let factor =
                                if dontInterpolate then
                                    match direction with
                                    | Left | Right -> pos.y

                                    | Top | Bottom -> pos.x
                                else
                                    match direction with
                                    | Right | Left ->
                                        0.0 + ((1.0 - x) / (y - x)) * (1.0 - 0.0)
                                    | Top | Bottom ->
                                        0.0 + ((1.0 - x) / (y - x)) * (1.0 - 0.0)
                            match direction with
                            | Top | Bottom -> factor + xF, pos.y + yF
                            | Left | Right -> pos.x + xF, factor + yF
                            |> tupleToPoint
    
                        (applyInterpolation dirA, applyInterpolation dirB)
                    )
    |]
    
let width, height = 64, 64
    
let metaballs =
    [| { radius = 3.4564; position = { x = 7.5; y = 8.0 } }
        { radius = 2.225433; position = { x = 5.5; y = 3.6 } }|]
    
let points = calculateOccupiedPoints metaballs width height (1./4.)
let lines = calculatePointsViaLinearInterpolation points width height false
let lines2 = calculatePointsViaLinearInterpolation points width height true

If I draw the results one get's this (ignore the slight artifacts, the underlying lib Avalonia probably has a bug). First pic shows the non interpolated variant, second the interpolated variant.