<p>as promised, here my answer to your question.</p>
<p>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.</p>
<p><a href="https://i.sstatic.net/atula.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/atula.png" alt="enter image description here" /></a></p>
<p>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 <em>how much part of the blob is within a grid block</em>. 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 &quot;blob&quot; 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).</p>
<p><a href="https://i.sstatic.net/FtXMR.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/FtXMR.png" alt="enter image description here" /></a></p>
<p>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 &quot;weights&quot; more.</p>
<p>More mathematically spoken:
Important to note here is that the <code>exLerp</code> function is returning the relative value for the <code>y</code> value of <code>Q</code>. This is always expected to be a value between <code>0.0</code> to <code>1.0</code>. This property is also guranteed as long as you don't mess up with supplying only pairs of values to your <code>exLerp</code> function with one of the values being below <code>1.0</code> (our treshold value) and the other value being above <code>1.0</code>.</p>
<p>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 <code>0.0</code> and <code>1.0</code>. 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).</p>
<p>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).</p>
<p>For your algorithm this means you need to find out at the code position when you call <code>interpolate</code> whether the return value regarding a specific point is meant to replace the <code>x</code> or <code>y</code> value of the point you're returning. I guess one variant could be here to check which of the values is <code>0.5</code> (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.</p>
<p>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.</p>
<p>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 ;) ).</p>
<pre class="lang-fsharp prettyprint-override"><code>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  |&gt; Array.sumBy(fun metaball -&gt; 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), [| |]
        |] |&gt; Array.mapi(fun i (_case, _) -&gt; (_case, firstHalf.[7 - i] |&gt; (fun (_, dirs: (Direction * Direction) []) -&gt; dirs)))
    Map.ofArray (firstHalf |&gt; Array.append secondHalf)
    
let calculatePointsViaLinearInterpolation (points: Map&lt;int * int, float&gt;) 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 |&gt; Array.map(fun (x0, y0) -&gt; (x0 + x, y0 + y))
                let edgeValues = edgePositions |&gt; Array.map(fun key -&gt; points.[key])
                let [| va; vb; vd; vc |] = edgeValues
                let edgeConfiguration =
                    let [| a; b; d; c |] = edgeValues |&gt; Array.map(fun value -&gt; value &gt;= 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 -&gt; mkPoint 0.5 0.0
                    | Bottom -&gt; mkPoint 0.5 1.0
                    | Left -&gt; mkPoint 0.0 0.5
                    | Right -&gt; mkPoint 1.0 0.5
                let pointValuesForDirection = function
                    | Top -&gt; va, vb
                    | Left -&gt; va, vc
                    | Right -&gt; vb, vd
                    | Bottom -&gt; vc, vd
    
                yield!
                    edgeConfiguration
                    |&gt; Array.map(fun (dirA, dirB) -&gt;
                        let applyInterpolation direction =
                            let (x, y) = pointValuesForDirection direction
                            let pos = getPosForDirection direction
                            let factor =
                                if dontInterpolate then
                                    match direction with
                                    | Left | Right -&gt; pos.y

                                    | Top | Bottom -&gt; pos.x
                                else
                                    match direction with
                                    | Right | Left -&gt;
                                        0.0 + ((1.0 - x) / (y - x)) * (1.0 - 0.0)
                                    | Top | Bottom -&gt;
                                        0.0 + ((1.0 - x) / (y - x)) * (1.0 - 0.0)
                            match direction with
                            | Top | Bottom -&gt; factor + xF, pos.y + yF
                            | Left | Right -&gt; pos.x + xF, factor + yF
                            |&gt; 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

</code></pre>
<p>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.</p>
<p><a href="https://i.sstatic.net/mlQB1.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/mlQB1.png" alt="non interpolated" /></a>
<a href="https://i.sstatic.net/NN8NS.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/NN8NS.png" alt="interpolated" /></a></p>
<p>EDIT:</p>
<p>So visually you can better see what happens if you show the values with interpolated colors for each tile.</p>
<p><a href="https://i.sstatic.net/lbUvY.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/lbUvY.png" alt="enter image description here" /></a></p>
<p>This one is using linear interpolation for the color values with every iso surface value above 1.0 being just blue.</p>
<p>Below the same picture without line.</p>
<p><a href="https://i.sstatic.net/10nbk.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/10nbk.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2021-06-06T05:26:36.923Z
Appreciate for your answer. I just posted my answer bellow with working code. I still have one question...
Arsenius
CC BY-SA 4.0
2021-06-03T23:01:47.253Z
2021-06-06T11:27:54.490Z
2021-06-06T11:27:54.490Z
realvictorprm
<p>Thanks @realvictorprm I figured out how to properly adjust the code:</p>
<p><a href="https://i.sstatic.net/gpM9s.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/gpM9s.png" alt="enter image description here" /></a></p>
<p><a href="https://i.sstatic.net/jIJoM.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/jIJoM.png" alt="enter image description here" /></a></p>
<pre><code>extends Node2D


class Blob:
    var pos_x
    var pos_y
    var radius
    var velocity
    func _init(x, y, r, v):
        pos_x = x
        pos_y = y
        radius = r
        velocity = v

const cell_size = 16
const blobs_count = 10
const blob_size = [20, 40]
var screen_size
var blobs
var allowUpdate = true
#var font = DynamicFont.new()

const drawMap = {
    0: null,
    1: [-0.5, 0, 0, -0.5], 
    2: [-0.5, -1, 0, -0.5], 
    3: [-0.5, 0, -0.5, -1], 
    4: [-1, -0.5, -0.5, -1], 
    5: [-1, -0.5, -0.5, 0, -0.5, -1, 0, -0.5],
    6: [-1, -0.5, 0, -0.5],
    7: [-1, -0.5, -0.5, 0],
    8: [-1, -0.5, -0.5, 0],
    9: [-1, -0.5, 0, -0.5],
    10: [-1, -0.5, -0.5, -1, -0.5, 0, 0, -0.5],
    11: [-1, -0.5, -0.5, -1],
    12: [-0.5, -1, -0.5, 0],
    13: [-0.5, -1, 0, -0.5],
    14: [-0.5, 0, 0, -0.5],
    15: null
}


func calcIsoSurface(x1, x2, y1, y2, r):
    var dx = abs(x1 - x2)
    var dy = abs(y1 - y2)
    var sd = dx*dx + dy*dy
    var res = float(r*r) / float(sd)
    return res

# Called when the node enters the scene tree for the first time.
func _ready():
    #font.font_data = load(&quot;res://font.ttf&quot;)
    #font.size = 7
    screen_size = get_viewport().size
    blobs = Array()
    var rng = RandomNumberGenerator.new()
    var r = rng.randi_range(blob_size[0], blob_size[1])
    var x = rng.randi_range(r, screen_size.x - r)
    var y = rng.randi_range(r, screen_size.y - r)
    for n in range(blobs_count):
        blobs.push_back(
            Blob.new(
                x, 
                y, 
                r, 
                Vector2(rng.randf_range(-3, 3), rng.randf_range(-3, 3))
            )
        )
    print(screen_size)


func formDrawIndex(x, y, sum, vertexes):
    var drawIndex = 0
    var corners = []
    if x &gt; 0 &amp;&amp; y &gt; 0:
        var vertex = vertexes.pop_front()
        
        if sum &gt;= 1:
            drawIndex |= 1

        if vertexes.back() &gt;= 1:
            drawIndex |= 2
            
        if vertex &gt;= 1:
            drawIndex |= 4  

        if vertexes.front() &gt;= 1:
            drawIndex |= 8
        
        corners.push_back(sum)
        corners.push_back(vertexes.back())
        corners.push_back(vertex)
        corners.push_back(vertexes.front())
    return  {&quot;draw_index&quot;: drawIndex, &quot;corners&quot;: corners}   


func exLerp(oneSum, zeroSum):
    if oneSum == zeroSum:
        return null
    return -(1-((1 - oneSum) / (zeroSum - oneSum)))


func interpolateLines(lines, corners):
    if lines == null:
        return lines
    for i in range(0, lines.size(), 2):
        var x = lines[i]
        var y = lines[i+1]
        
        if (x == 0 || x == -1) &amp;&amp; (y == 0 || y == -1):
            continue
        if (x == 0 || x == -1):
            lines[i+1] = exLerp(corners[1], corners[0]) if x == 0 else exLerp(corners[2], corners[3])
        if (y == 0 || y == -1):
            lines[i] = exLerp(corners[3], corners[0]) if y == 0 else exLerp(corners[2], corners[1])
    
    return lines


func drawLines(x, y, lines):
    if lines != null &amp;&amp; lines.size() &gt;= 4:
        draw_line(
            Vector2(x + (cell_size*lines[0]), y + (cell_size*lines[1])), 
            Vector2(x + (cell_size*lines[2]), y + (cell_size*lines[3])), 
            Color.green
        )

    if lines != null &amp;&amp; lines.size() == 8:
        draw_line(
            Vector2(x + (cell_size*lines[4]), y + (cell_size*lines[5])), 
            Vector2(x + (cell_size*lines[6]), y + (cell_size*lines[7])), 
            Color.green
        )


# Called after update() in the _process()
func _draw():
    var vertexes = []
    for x in range(0, screen_size.x, cell_size):
        for y in range(0, screen_size.y, cell_size):
            var sum = 0
            for blob in blobs:
                sum += calcIsoSurface(x, blob.pos_x, y, blob.pos_y, blob.radius)
                #var c = Color.blue
                #c.a = 0.0001
                #draw_circle(Vector2(blob.pos_x,blob.pos_y), blob.radius, c)
            
            #if sum &gt;= 1:   
            #   draw_rect(Rect2(x, y, 1, 1), Color.red)
            #else:
            #   draw_rect(Rect2(x, y, 1, 1), Color.black)
        
            #draw_string(font, Vector2(x, y), &quot;%s&quot; % sum, Color.black)

            var indexies = formDrawIndex(x, y, sum, vertexes)
            var lines = drawMap[indexies[&quot;draw_index&quot;]]
            var corners = indexies[&quot;corners&quot;]
            lines = interpolateLines(lines, corners)
            
            drawLines(x, y, lines)
            vertexes.push_back(sum)

        if x &gt; 0:
            vertexes.pop_front()



func _input(event):
    if event is InputEventMouseButton &amp;&amp; event.is_pressed():
        allowUpdate = !allowUpdate
        print(allowUpdate)


# Called every frame. 'delta' is the elapsed time since the previous frame.

func _process(delta):
    if !allowUpdate:        
        return      
    update()

    for blob in blobs:
        blob.pos_x += blob.velocity.x
        blob.pos_y += blob.velocity.y
        if blob.pos_x &gt; screen_size.x || blob.pos_x &lt; 0:
            blob.velocity.x *= -1
        if blob.pos_y &gt; screen_size.y || blob.pos_y &lt; 0:
            blob.velocity.y *= -1
</code></pre>
<p>First calculation of corners weight has to be float not integer:</p>
<pre><code>    func calcIsoSurface(x1, x2, y1, y2, r):
        var dx = abs(x1 - x2)
        var dy = abs(y1 - y2)
        var sd = dx*dx + dy*dy
        var res = float(r*r) / float(sd)
        return res 
</code></pre>
<p>The second lerp has to use inversion (since my coordinate system from right to left and from bottom to top):</p>
<pre><code>func exLerp(oneSum, zeroSum):
    if oneSum == zeroSum:
        return null
    return -(1-((1 - oneSum) / (zeroSum - oneSum)))
</code></pre>
<p>And finally interpolation itself:</p>
<pre><code>func interpolateLines(lines, corners):
    if lines == null:
        return lines
    for i in range(0, lines.size(), 2):
        var x = lines[i]
        var y = lines[i+1]

        if (x == 0 || x == -1) &amp;&amp; (y == 0 || y == -1):
            continue
        if (x == 0 || x == -1):
            lines[i+1] = exLerp(corners[1], corners[0]) if x == 0 else exLerp(corners[2], corners[3])
        if (y == 0 || y == -1):
            lines[i] = exLerp(corners[3], corners[0]) if y == 0 else exLerp(corners[2], corners[1])

    return lines
</code></pre>
<p>But I still not fully understand how it works. Why ratio of metaball radius and distance from centre of metaball to corner of marching square is correctly affect x or y coordinate for interpolation when we apply lerp? I can't fully understand it visually. Another word how <code>calcIsoSurface</code> relate to <code>exLerp</code> in order to figure out correct <code>x</code> or <code>y</code> for interpolation?</p>

CC BY-SA 4.0
2021-06-06T10:15:59.387Z
I will edit my answer to adress your last question
realvictorprm
CC BY-SA 4.0
2021-06-06T10:57:57.857Z
I'm not exactly sure what I else I can tell you or show you visually. Does the last image I added to my answer help you somehow @Arsenius?
realvictorprm
CC BY-SA 4.0
2021-06-11T07:49:54.547Z
@realvictorprm thanks for edit. Visually I mean show square and isosurface and how corner weights (sum of ratio between metaball radius and distance to square corner from centre of metaball) affect line curving when this interpolation formula is applying   (1 - cornerWeight1) / (cornerWeight2 - cornerWeight1) 
Good visual example shows here https://www.youtube.com/watch?v=wnl0O_xkgzs why a^2+b^2 = c^2 this useful when we calculate distance in 2d from one point to in another in the space. If u have time few slides to show all this relations to visually understand how it all works.
Arsenius
CC BY-SA 4.0
2021-06-06T05:00:32.473Z
2021-06-06T05:19:23.523Z
2021-06-06T05:19:23.523Z
Arsenius
<p>I struggle do understand how linear interpolation works in the marching square rendering algorithm context.</p>

<p>I created simple example in GDScript (Godot) of random floating metaballs to demonstrate the problem:</p>

<pre><code>extends Node2D

class Blob:
    var pos_x
    var pos_y
    var radius
    var velocity
    func _init(x, y, r, v):
        pos_x = x
        pos_y = y
        radius = r
        velocity = v

# Declare member variables here. Examples:
# const a = 2
# var b = "text"

const cell_size = 16
const blobs_count = 10
const blob_size = [20, 40]
const max_sum = 1
var screen_size
var blobs
var allowUpdate = true

const drawMap = {
    0: null,
    1: [-0.5, 0, 0, -0.5],
    2: [-0.5, -1, 0, -0.5],
    3: [-0.5, 0, -0.5, -1],
    4: [-1, -0.5, -0.5, -1],
    5: [-1, -0.5, -0.5, 0, -0.5, -1, 0, -0.5],
    6: [-1, -0.5, 0, -0.5],
    7: [-1, -0.5, -0.5, 0],
    8: [-1, -0.5, -0.5, 0],
    9: [-1, -0.5, 0, -0.5],
    10: [-1, -0.5, -0.5, -1, -0.5, 0, 0, -0.5],
    11: [-1, -0.5, -0.5, -1],
    12: [-0.5, -1, -0.5, 0],
    13: [-0.5, -1, 0, -0.5],
    14: [-0.5, 0, 0, -0.5],
    15: null
}


func calcIsoSurface(x1, x2, y1, y2, r):
    var dx = x1 - x2
    var dy = y1 - y2
    var sd = dx*dx + dy*dy
    var res = r*r / sd
    return res

# Called when the node enters the scene tree for the first time.
func _ready():
    screen_size = get_viewport().size
    blobs = Array()
    var rng = RandomNumberGenerator.new()
    var r = rng.randi_range(blob_size[0], blob_size[1])
    var x = rng.randi_range(r, screen_size.x - r)
    var y = rng.randi_range(r, screen_size.y - r)
    for n in range(blobs_count):
        blobs.push_back(
            Blob.new(
                x, 
                y, 
                r, 
                Vector2(rng.randf_range(-1, 1), rng.randf_range(-1, 1))
            )
        )
    print(screen_size)


func formDrawIndex(x, y, sum, vertexes):
    var drawIndex = 0
    var corners = []
    if x &gt; 0 &amp;&amp; y &gt; 0:
        if sum &gt;= 1:
            drawIndex |= 1

        corners.push_back(sum)

        if vertexes.back() &gt;= 1:
            drawIndex |= 2

        corners.push_back(vertexes.back())
        corners.push_back(vertexes.pop_front())

        if corners.back() &gt;= 1:
            drawIndex |= 4  
        if vertexes.front() &gt;= 1:
            drawIndex |= 8

        corners.push_back(vertexes.front())
    return  {"draw_index": drawIndex, "corners": corners}   


func exLerp(oneSum, zeroSum):
    if oneSum == zeroSum:
        return null
    return (1 - oneSum) / (zeroSum - oneSum)


func interpolateLines(lines, corners):
    if lines == null:
        return lines
    for i in range(0, lines.size(), 2):
        var x = lines[i]
        var y = lines[i+1]
        #somehow implement correct interpolation here
    return lines


func drawLines(x, y, lines):
    if lines != null &amp;&amp; lines.size() &gt;= 4:
        draw_line(
            Vector2(x + (cell_size*lines[0]), y + (cell_size*lines[1])), 
            Vector2(x + (cell_size*lines[2]), y + (cell_size*lines[3])), 
            Color.green
        )

    if lines != null &amp;&amp; lines.size() == 8:
        draw_line(
            Vector2(x + (cell_size*lines[4]), y + (cell_size*lines[5])), 
            Vector2(x + (cell_size*lines[6]), y + (cell_size*lines[7])), 
            Color.green
        )

# Called after update() in the _process()
func _draw():
    var vertexes = []
    for x in range(0, screen_size.x, cell_size):
        for y in range(0, screen_size.y, cell_size):
            var sum = 0
            for blob in blobs:
                sum += calcIsoSurface(x, blob.pos_x, y, blob.pos_y, blob.radius)
                #draw_circle(Vector2(blob.pos_x,blob.pos_y), blob.radius, Color.darkgray)

            #if sum &gt;= 1:
                #draw_rect(Rect2(x, y, 1, 1), Color.red)
            #else:  draw_rect(Rect2(x, y, 1, 1), Color.black)

            var indexies = formDrawIndex(x, y, sum, vertexes)
            var lines = drawMap[indexies["draw_index"]]
            lines = interpolateLines(lines, indexies["corners"])            
            drawLines(x, y, lines)

            vertexes.push_back(sum)

        if x &gt; 0:
            vertexes.pop_front()



func _input(event):
    if event is InputEventMouseButton &amp;&amp; event.is_pressed():
        allowUpdate = !allowUpdate
        print(allowUpdate)


# Called every frame. 'delta' is the elapsed time since the previous frame.
func _process(delta):
    if !allowUpdate:        
        return      
    update()

    for blob in blobs:
        blob.pos_x += blob.velocity.x
        blob.pos_y += blob.velocity.y
        if blob.pos_x &gt; screen_size.x || blob.pos_x &lt; 0:
            blob.velocity.x *= -1
        if blob.pos_y &gt; screen_size.y || blob.pos_y &lt; 0:
            blob.velocity.y *= -1
</code></pre>

<p>The outcome looks something like this:
<a href="https://i.sstatic.net/vCcqM.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/vCcqM.png" alt="enter image description here"></a></p>

<p>Now I would like to apply linear interpolation to make my meatballs smoother. This is where I stuck. The desire outcome is transform rendering of this:
<a href="https://i.sstatic.net/t1Nve.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/t1Nve.png" alt="enter image description here"></a></p>

<p>To this:
<a href="https://i.sstatic.net/v6V2m.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/v6V2m.png" alt="enter image description here"></a></p>

<p>I already created <code>interpolateLines</code> and <code>exLerp</code> functions, but don't understand how exactly in this context linear interpolation works and why?
I used this <a href="http://jamie-wong.com/2014/08/19/metaballs-and-marching-squares/" rel="nofollow noreferrer">material</a> as a theory for code implementation, but last part is still blurry for me.</p>

<p>Anyone can provide working code example and explain theory in dump language with more deep math explanation and visual demonstration?
So I can finally breakthrough this problem.   </p>

<rendering><algorithm><2d-graphics><physics><fluid-sim>
2D metaballs with marching squares and linear interpolation
Arsenius
# 2D metaballs with marching squares and linear interpolation

(Asked by: username_0 on 2020-05-29)

I struggle do understand how linear interpolation works in the marching square rendering algorithm context.

I created simple example in GDScript (Godot) of random floating metaballs to demonstrate the problem:

```
extends Node2D

class Blob:
    var pos_x
    var pos_y
    var radius
    var velocity
    func _init(x, y, r, v):
        pos_x = x
        pos_y = y
        radius = r
        velocity = v

# Declare member variables here. Examples:
# const a = 2
# var b = "text"

const cell_size = 16
const blobs_count = 10
const blob_size = [20, 40]
const max_sum = 1
var screen_size
var blobs
var allowUpdate = true

const drawMap = {
    0: null,
    1: [-0.5, 0, 0, -0.5],
    2: [-0.5, -1, 0, -0.5],
    3: [-0.5, 0, -0.5, -1],
    4: [-1, -0.5, -0.5, -1],
    5: [-1, -0.5, -0.5, 0, -0.5, -1, 0, -0.5],
    6: [-1, -0.5, 0, -0.5],
    7: [-1, -0.5, -0.5, 0],
    8: [-1, -0.5, -0.5, 0],
    9: [-1, -0.5, 0, -0.5],
    10: [-1, -0.5, -0.5, -1, -0.5, 0, 0, -0.5],
    11: [-1, -0.5, -0.5, -1],
    12: [-0.5, -1, -0.5, 0],
    13: [-0.5, -1, 0, -0.5],
    14: [-0.5, 0, 0, -0.5],
    15: null
}

func calcIsoSurface(x1, x2, y1, y2, r):
    var dx = x1 - x2
    var dy = y1 - y2
    var sd = dx*dx + dy*dy
    var res = r*r / sd
    return res

# Called when the node enters the scene tree for the first time.
func _ready():
    screen_size = get_viewport().size
    blobs = Array()
    var rng = RandomNumberGenerator.new()
    var r = rng.randi_range(blob_size[0], blob_size[1])
    var x = rng.randi_range(r, screen_size.x - r)
    var y = rng.randi_range(r, screen_size.y - r)
    for n in range(blobs_count):
        blobs.push_back(
            Blob.new(
                x, 
                y, 
                r, 
                Vector2(rng.randf_range(-1, 1), rng.randf_range(-1, 1))
            )
        )
    print(screen_size)

func formDrawIndex(x, y, sum, vertexes):
    var drawIndex = 0
    var corners = []
    if x > 0 && y > 0:
        if sum >= 1:
            drawIndex |= 1

        corners.push_back(sum)

        if vertexes.back() >= 1:
            drawIndex |= 2

        corners.push_back(vertexes.back())
        corners.push_back(vertexes.pop_front())

        if corners.back() >= 1:
            drawIndex |= 4  
        if vertexes.front() >= 1:
            drawIndex |= 8

        corners.push_back(vertexes.front())
    return  {"draw_index": drawIndex, "corners": corners}   

func exLerp(oneSum, zeroSum):
    if oneSum == zeroSum:
        return null
    return (1 - oneSum) / (zeroSum - oneSum)

func interpolateLines(lines, corners):
    if lines == null:
        return lines
    for i in range(0, lines.size(), 2):
        var x = lines[i]
        var y = lines[i+1]
        #somehow implement correct interpolation here
    return lines

func drawLines(x, y, lines):
    if lines != null && lines.size() >= 4:
        draw_line(
            Vector2(x + (cell_size*lines[0]), y + (cell_size*lines[1])), 
            Vector2(x + (cell_size*lines[2]), y + (cell_size*lines[3])), 
            Color.green
        )

    if lines != null && lines.size() == 8:
        draw_line(
            Vector2(x + (cell_size*lines[4]), y + (cell_size*lines[5])), 
            Vector2(x + (cell_size*lines[6]), y + (cell_size*lines[7])), 
            Color.green
        )

# Called after update() in the _process()
func _draw():
    var vertexes = []
    for x in range(0, screen_size.x, cell_size):
        for y in range(0, screen_size.y, cell_size):
            var sum = 0
            for blob in blobs:
                sum += calcIsoSurface(x, blob.pos_x, y, blob.pos_y, blob.radius)
                #draw_circle(Vector2(blob.pos_x,blob.pos_y), blob.radius, Color.darkgray)

            #if sum >= 1:
                #draw_rect(Rect2(x, y, 1, 1), Color.red)
            #else:  draw_rect(Rect2(x, y, 1, 1), Color.black)

            var indexies = formDrawIndex(x, y, sum, vertexes)
            var lines = drawMap[indexies["draw_index"]]
            lines = interpolateLines(lines, indexies["corners"])            
            drawLines(x, y, lines)

            vertexes.push_back(sum)

        if x > 0:
            vertexes.pop_front()

func _input(event):
    if event is InputEventMouseButton && event.is_pressed():
        allowUpdate = !allowUpdate
        print(allowUpdate)

# Called every frame. 'delta' is the elapsed time since the previous frame.
func _process(delta):
    if !allowUpdate:        
        return      
    update()

    for blob in blobs:
        blob.pos_x += blob.velocity.x
        blob.pos_y += blob.velocity.y
        if blob.pos_x > screen_size.x || blob.pos_x < 0:
            blob.velocity.x *= -1
        if blob.pos_y > screen_size.y || blob.pos_y < 0:
            blob.velocity.y *= -1

```

The outcome looks something like this:
[![enter image description here](https://i.sstatic.net/vCcqM.png)](https://i.sstatic.net/vCcqM.png)

Now I would like to apply linear interpolation to make my meatballs smoother. This is where I stuck. The desire outcome is transform rendering of this:
[![enter image description here](https://i.sstatic.net/t1Nve.png)](https://i.sstatic.net/t1Nve.png)

To this:
[![enter image description here](https://i.sstatic.net/v6V2m.png)](https://i.sstatic.net/v6V2m.png)

I already created `interpolateLines` and `exLerp` functions, but don't understand how exactly in this context linear interpolation works and why?
I used this [material](http://jamie-wong.com/2014/08/19/metaballs-and-marching-squares/) as a theory for code implementation, but last part is still blurry for me.

Anyone can provide working code example and explain theory in dump language with more deep math explanation and visual demonstration?
So I can finally breakthrough this problem.

**Tags:** rendering, algorithm, 2d-graphics, physics, fluid-sim


---

## 2 answers

### Answer by username_1 on 2021-06-03. Score: 3

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.

[![enter image description here](https://i.sstatic.net/atula.png)](https://i.sstatic.net/atula.png)

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

[![enter image description here](https://i.sstatic.net/FtXMR.png)](https://i.sstatic.net/FtXMR.png)

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

```fsharp
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.

[![non interpolated](https://i.sstatic.net/mlQB1.png)](https://i.sstatic.net/mlQB1.png)
[![interpolated](https://i.sstatic.net/NN8NS.png)](https://i.sstatic.net/NN8NS.png)

EDIT:

So visually you can better see what happens if you show the values with interpolated colors for each tile.

[![enter image description here](https://i.sstatic.net/lbUvY.png)](https://i.sstatic.net/lbUvY.png)

This one is using linear interpolation for the color values with every iso surface value above 1.0 being just blue.

Below the same picture without line.

[![enter image description here](https://i.sstatic.net/10nbk.png)](https://i.sstatic.net/10nbk.png)

> **Comments:**
> **username_0**: Appreciate for your answer. I just posted my answer bellow with working code. I still have one question...  

---

### Answer by username_0 on 2021-06-06. Score: 0

Thanks @username_1 I figured out how to properly adjust the code:

[![enter image description here](https://i.sstatic.net/gpM9s.png)](https://i.sstatic.net/gpM9s.png)

[![enter image description here](https://i.sstatic.net/jIJoM.png)](https://i.sstatic.net/jIJoM.png)

```
extends Node2D

class Blob:
    var pos_x
    var pos_y
    var radius
    var velocity
    func _init(x, y, r, v):
        pos_x = x
        pos_y = y
        radius = r
        velocity = v

const cell_size = 16
const blobs_count = 10
const blob_size = [20, 40]
var screen_size
var blobs
var allowUpdate = true
#var font = DynamicFont.new()

const drawMap = {
    0: null,
    1: [-0.5, 0, 0, -0.5], 
    2: [-0.5, -1, 0, -0.5], 
    3: [-0.5, 0, -0.5, -1], 
    4: [-1, -0.5, -0.5, -1], 
    5: [-1, -0.5, -0.5, 0, -0.5, -1, 0, -0.5],
    6: [-1, -0.5, 0, -0.5],
    7: [-1, -0.5, -0.5, 0],
    8: [-1, -0.5, -0.5, 0],
    9: [-1, -0.5, 0, -0.5],
    10: [-1, -0.5, -0.5, -1, -0.5, 0, 0, -0.5],
    11: [-1, -0.5, -0.5, -1],
    12: [-0.5, -1, -0.5, 0],
    13: [-0.5, -1, 0, -0.5],
    14: [-0.5, 0, 0, -0.5],
    15: null
}

func calcIsoSurface(x1, x2, y1, y2, r):
    var dx = abs(x1 - x2)
    var dy = abs(y1 - y2)
    var sd = dx*dx + dy*dy
    var res = float(r*r) / float(sd)
    return res

# Called when the node enters the scene tree for the first time.
func _ready():
    #font.font_data = load("res://font.ttf")
    #font.size = 7
    screen_size = get_viewport().size
    blobs = Array()
    var rng = RandomNumberGenerator.new()
    var r = rng.randi_range(blob_size[0], blob_size[1])
    var x = rng.randi_range(r, screen_size.x - r)
    var y = rng.randi_range(r, screen_size.y - r)
    for n in range(blobs_count):
        blobs.push_back(
            Blob.new(
                x, 
                y, 
                r, 
                Vector2(rng.randf_range(-3, 3), rng.randf_range(-3, 3))
            )
        )
    print(screen_size)

func formDrawIndex(x, y, sum, vertexes):
    var drawIndex = 0
    var corners = []
    if x > 0 && y > 0:
        var vertex = vertexes.pop_front()
        
        if sum >= 1:
            drawIndex |= 1

        if vertexes.back() >= 1:
            drawIndex |= 2
            
        if vertex >= 1:
            drawIndex |= 4  

        if vertexes.front() >= 1:
            drawIndex |= 8
        
        corners.push_back(sum)
        corners.push_back(vertexes.back())
        corners.push_back(vertex)
        corners.push_back(vertexes.front())
    return  {"draw_index": drawIndex, "corners": corners}   

func exLerp(oneSum, zeroSum):
    if oneSum == zeroSum:
        return null
    return -(1-((1 - oneSum) / (zeroSum - oneSum)))

func interpolateLines(lines, corners):
    if lines == null:
        return lines
    for i in range(0, lines.size(), 2):
        var x = lines[i]
        var y = lines[i+1]
        
        if (x == 0 || x == -1) && (y == 0 || y == -1):
            continue
        if (x == 0 || x == -1):
            lines[i+1] = exLerp(corners[1], corners[0]) if x == 0 else exLerp(corners[2], corners[3])
        if (y == 0 || y == -1):
            lines[i] = exLerp(corners[3], corners[0]) if y == 0 else exLerp(corners[2], corners[1])
    
    return lines

func drawLines(x, y, lines):
    if lines != null && lines.size() >= 4:
        draw_line(
            Vector2(x + (cell_size*lines[0]), y + (cell_size*lines[1])), 
            Vector2(x + (cell_size*lines[2]), y + (cell_size*lines[3])), 
            Color.green
        )

    if lines != null && lines.size() == 8:
        draw_line(
            Vector2(x + (cell_size*lines[4]), y + (cell_size*lines[5])), 
            Vector2(x + (cell_size*lines[6]), y + (cell_size*lines[7])), 
            Color.green
        )

# Called after update() in the _process()
func _draw():
    var vertexes = []
    for x in range(0, screen_size.x, cell_size):
        for y in range(0, screen_size.y, cell_size):
            var sum = 0
            for blob in blobs:
                sum += calcIsoSurface(x, blob.pos_x, y, blob.pos_y, blob.radius)
                #var c = Color.blue
                #c.a = 0.0001
                #draw_circle(Vector2(blob.pos_x,blob.pos_y), blob.radius, c)
            
            #if sum >= 1:   
            #   draw_rect(Rect2(x, y, 1, 1), Color.red)
            #else:
            #   draw_rect(Rect2(x, y, 1, 1), Color.black)
        
            #draw_string(font, Vector2(x, y), "%s" % sum, Color.black)

            var indexies = formDrawIndex(x, y, sum, vertexes)
            var lines = drawMap[indexies["draw_index"]]
            var corners = indexies["corners"]
            lines = interpolateLines(lines, corners)
            
            drawLines(x, y, lines)
            vertexes.push_back(sum)

        if x > 0:
            vertexes.pop_front()

func _input(event):
    if event is InputEventMouseButton && event.is_pressed():
        allowUpdate = !allowUpdate
        print(allowUpdate)

# Called every frame. 'delta' is the elapsed time since the previous frame.

func _process(delta):
    if !allowUpdate:        
        return      
    update()

    for blob in blobs:
        blob.pos_x += blob.velocity.x
        blob.pos_y += blob.velocity.y
        if blob.pos_x > screen_size.x || blob.pos_x < 0:
            blob.velocity.x *= -1
        if blob.pos_y > screen_size.y || blob.pos_y < 0:
            blob.velocity.y *= -1

```

First calculation of corners weight has to be float not integer:

```
    func calcIsoSurface(x1, x2, y1, y2, r):
        var dx = abs(x1 - x2)
        var dy = abs(y1 - y2)
        var sd = dx*dx + dy*dy
        var res = float(r*r) / float(sd)
        return res 

```

The second lerp has to use inversion (since my coordinate system from right to left and from bottom to top):

```
func exLerp(oneSum, zeroSum):
    if oneSum == zeroSum:
        return null
    return -(1-((1 - oneSum) / (zeroSum - oneSum)))

```

And finally interpolation itself:

```
func interpolateLines(lines, corners):
    if lines == null:
        return lines
    for i in range(0, lines.size(), 2):
        var x = lines[i]
        var y = lines[i+1]

        if (x == 0 || x == -1) && (y == 0 || y == -1):
            continue
        if (x == 0 || x == -1):
            lines[i+1] = exLerp(corners[1], corners[0]) if x == 0 else exLerp(corners[2], corners[3])
        if (y == 0 || y == -1):
            lines[i] = exLerp(corners[3], corners[0]) if y == 0 else exLerp(corners[2], corners[1])

    return lines

```

But I still not fully understand how it works. Why ratio of metaball radius and distance from centre of metaball to corner of marching square is correctly affect x or y coordinate for interpolation when we apply lerp? I can't fully understand it visually. Another word how `calcIsoSurface` relate to `exLerp` in order to figure out correct `x` or `y` for interpolation?

> **Comments:**
> **username_1**: I will edit my answer to adress your last question  
> **username_1**: I'm not exactly sure what I else I can tell you or show you visually. Does the last image I added to my answer help you somehow @username_0?  
> **username_0**: @username_1 thanks for edit. Visually I mean show square and isosurface and how corner weights (sum of ratio between metaball radius and distance to square corner from centre of metaball) affect line curving when this interpolation formula is applying   (1 - cornerWeight1) / (cornerWeight2 - cornerWeight1)  Good visual example shows here https://www.youtube.com/watch?v=wnl0O_xkgzs why a^2+b^2 = c^2 this useful when we calculate distance in 2d from one point to in another in the space. If u have time few slides to show all this relations to visually understand how it all works.  

---

