# 2D metaballs with marching squares and linear interpolation

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 velocity
func _init(x, y, r, v):
pos_x = x
pos_y = y
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.
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)

#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:

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:

To this:

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

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.

EDIT:

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

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.

• Appreciate for your answer. I just posted my answer bellow with working code. I still have one question... Commented Jun 6, 2021 at 5:26

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

extends Node2D

class Blob:
var pos_x
var pos_y
var velocity
func _init(x, y, r, v):
pos_x = x
pos_y = y
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.
#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

#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?

• I will edit my answer to adress your last question Commented Jun 6, 2021 at 10:15
• 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? Commented Jun 6, 2021 at 10:57
• @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 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. Commented Jun 11, 2021 at 7:49