Metaballs: blobs of goo

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Post Reply
angros47
Posts: 2321
Joined: Jun 21, 2005 19:04

Metaballs: blobs of goo

Post by angros47 »

A metaball, or blob, is a spherical object that can be merged with other blobs; it can be useful to represent atoms and molecules, or cells (mitosis of a cell can be represented by separating two metaballs).

Here is an example, in 2d (there are also examples in 3d)

http://www.youtube.com/watch?v=3SLsB32a_FI

To draw metaballs, the trick is: first of all, imagine a "force field" that is inversely proportional to the squared distance to the center of the metaball: when it drop under a threshold, the surface of the metaball is drawn; if two metaballs are too close, between them the sum of their force fields will be higher than the threshold, and so they will merge.

The simplest way to draw metaballs is to check, for every point on the screen, if its force field (the sum of the field of all metaballs) is above the threshold or not:

Code: Select all

Dim metaballX(5) as integer
Dim metaballY(5) as integer


randomize timer

screenres 800,600,32

dim as ubyte ptr target

for i as integer=1 to 5
	metaballX(i)=800*rnd(1)
	metaballY(i)=600*rnd(1)
next


do
	target = screenPtr 
	for i as integer=1 to 5
		metaballX(i)+=10*rnd(1)-5
		metaballY(i)+=10*rnd(1)-5
	next

	screenLock 
	for y as integer = 1 to 600
		for x as integer = 1 to 800

			dim d as single
			d=0

			for i as integer=1 to 5
				d+=500.0/((x-metaballX(i))^2+(y-metaballY(i))^2)
			next


			
			if d>.5 then
				*target = 255    '  Blue level.
				*(target+1) = 255    ' Green level.
				*(target+2) = 255    ' Red level.
			else
				*target = 0    '  Blue level.
				*(target+1) = 0    ' Green level.
				*(target+2) = 0    ' Red level.
			end if
			target += 4
		next 
	next
	screenUnlock
	sleep 1

loop until inkey<>""

sleep
As you can see, it works, but it's awfully slow.

Is there a faster solution? Of course: you can evaluate just some points, and interpolate the rest, with an algorithm called "marching squares" (http://en.wikipedia.org/wiki/Marching_squares):

Here is the result:

Code: Select all

Dim shared metaballX(5) as integer
Dim shared metaballY(5) as integer


randomize timer

screenres 800,600,32,2
ScreenSet 1,0


for i as integer=1 to 5
	metaballX(i)=800*rnd(1)
	metaballY(i)=600*rnd(1)
next


function scalarfield (x as integer, y as integer) as single
	dim d as single

	for i as integer=1 to 5
		d+=500.0/((x-metaballX(i))^2+(y-metaballY(i))^2)
	next

	return d
end function



do
	cls

	for i as integer=1 to 5
		metaballX(i)+=2*rnd(1)-1
		metaballY(i)+=2*rnd(1)-1
	next

	for y as integer = 1 to 600 step 20
		for x as integer = 1 to 800 step 20


			dim as single A,B,C,D
			A=scalarfield(x,y)
			B=scalarfield(x+20,y)
			C=scalarfield(x,y+20)
			D=scalarfield(x+20,y+20)

			dim as integer side
			if A<.5 then side or=1
			if B<.5 then side or=2
			if C<.5 then side or=4
			if D<.5 then side or=8

			dim as integer x1,y1,x2,y2
			select case side
			case 1, 9, 14:
				x1=x+20*(.5-A)/(B-A)
				y1=y+20*(.5-A)/(C-A)
				line(x,y1)-(x1,y)
			case 2, 6, 13:
				x1=x-20*(.5-B)/(A-B)+20
				y1=y+20*(.5-B)/(D-B)
				line(x1,y)-(x+20,y1)
			case 3, 12:
				y1=y+20*(.5-A)/(C-A)
				y2=y+20*(.5-B)/(D-B)
				line(x,y1)-(x+20,y2)
			case 4,11:
				x1=x+20*(.5-C)/(D-C)
				y1=y-20*(.5-C)/(A-C)+20
				line(x,y1)-(x1,y+20)
			case 5,10:
				x1=x+20*(.5-A)/(B-A)
				x2=x+20*(.5-C)/(D-C)
				line(x1,y)-(x2,y+20)
			case 7,8:
				x1=x-20*(.5-D)/(C-D)+20
				y1=y-20*(.5-D)/(B-D)+20
				line(x+20,y1)-(x1,y+20)
				
			end select
			
		next 
	next
	flip 1,0
	sleep 1

loop until inkey<>""

sleep
badidea
Posts: 2586
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Metaballs: blobs of goo

Post by badidea »

A bit dusty code, but exactly what I was looking for (marching squares). There is however a bug in it as can be better seen it this slightly modified version:

Code: Select all

const N_BALL = 20
Dim shared metaballX(1 to N_BALL) as integer
Dim shared metaballY(1 to N_BALL) as integer

'randomize timer
randomize 27

screenres 800,600,32',2
width 800\8,600\16
'ScreenSet 1,0

for i as integer = 1 to N_BALL
	metaballX(i) = 800*rnd(1)
	metaballY(i) = 600*rnd(1)
next

function scalarfield (x as integer, y as integer) as single
	dim d as single
	for i as integer = 1 to N_BALL
		d += 500.0 / ((x-metaballX(i))^2 + (y-metaballY(i))^2)
	next
	return d
end function

dim as long key
do
	for i as integer = 1 to N_BALL
		metaballX(i) += 2*rnd(1)-1
		metaballY(i) += 2*rnd(1)-1
	next

	screenlock
	cls
	for y as integer = 1 to 600 step 20
		for x as integer = 1 to 800 step 20
			pset (x,y), rgb(125,125,125) 'grey

			dim as single A,B,C,D
			A = scalarfield(x,y)
			B = scalarfield(x+20,y)
			C = scalarfield(x,y+20)
			D = scalarfield(x+20,y+20)

			dim as integer side
			if A < .5 then side or= 1
			if B < .5 then side or= 2
			if C < .5 then side or= 4
			if D < .5 then side or= 8

			dim as integer x1,y1,x2,y2
			select case side
			case 1, 9, 14: '0001, 1001, 1110
				x1 = x+20*(.5-A)/(B-A)
				y1 = y+20*(.5-A)/(C-A)
				line(x,y1)-(x1,y), rgb(255,0,0) 'red
			case 2, 6, 13: '0010, 0110, 1101
				x1 = x-20*(.5-B)/(A-B)+20
				y1 = y+20*(.5-B)/(D-B)
				line(x1,y)-(x+20,y1), rgb(255,255,0) 'yellow
			case 3, 12: '0011, 1100
				y1 = y+20*(.5-A)/(C-A)
				y2 = y+20*(.5-B)/(D-B)
				line(x,y1)-(x+20,y2), rgb(0,255,0) 'green
			case 4, 11: '0100, 1011
				x1 = x+20*(.5-C)/(D-C)
				y1 = y-20*(.5-C)/(A-C)+20
				line(x,y1)-(x1,y+20), rgb(0,255,255) 'cyan
			case 5, 10: '0101, 1010
				x1 = x+20*(.5-A)/(B-A)
				x2 = x+20*(.5-C)/(D-C)
				line(x1,y)-(x2,y+20), rgb(0,0,255) 'blue
			case 7, 8: '0111, 1000
				x1 = x-20*(.5-D)/(C-D)+20
				y1 = y-20*(.5-D)/(B-D)+20
				line(x+20,y1)-(x1,y+20), rgb(255,255,255) 'white
			end select
		next 
	next
	draw string (5,5), "Print <escape> to exit, <any key> to change positions"
	screenunlock

	'flip 1,0
	key = getkey

loop while key <> 27
print "End."
See how the 2 balls at the left connect at the start.
Could be the 'saddle point ambiguity' that the Wikipedia article mentioned. I'll look further into this tomorrow.
UEZ
Posts: 972
Joined: May 05, 2017 19:59
Location: Germany

Re: Metaballs: blobs of goo

Post by UEZ »

Nice.
angros47
Posts: 2321
Joined: Jun 21, 2005 19:04

Re: Metaballs: blobs of goo

Post by angros47 »

I confirm. The cause of the bug is that every square can contain only one line (only one case is executed per iteration, and each case can draw only one line). If you look at the squares, you can see how the one affected by the bug should contain two lines, to render the balls correctly.

To fix the bug, a special routine should be added to handle that specific case, but it would likely be too time consuming. An alternative might perhaps be replacing the marching squares with marching triangles, since no ambiguity would be possible with triangles
UEZ
Posts: 972
Joined: May 05, 2017 19:59
Location: Germany

Re: Metaballs: blobs of goo

Post by UEZ »

Maybe you know it: you can use the trick with blur and contrast to achieve similar effect.

Code: Select all

'Coded by UEZ

#Include "fbgfx.bi"
Using FB

'https://gist.github.com/mattdesl/4383372
'Copyright (c) 2007, Romain Guy All rights reserved.
Type tImage 
    As Integer width, height, pitch
    As Long Ptr pixels
End Type

Function BlurPass(img As tImage, iRadius As Ubyte, iW As Ulong, iH As Ulong) As Any Ptr
    Dim As Ulong iW1 = iW - 1, iH1 = iH - 1
   
	Dim pImage_blurred As Any Ptr = Imagecreate(img.height, img.width, 0, 32)

	Dim As tImage img_b
	Imageinfo(pImage_blurred, img_b.Width, img_b.height, , img_b.pitch, img_b.pixels)
	
	Dim As Long previousPixelIndex, sumAlpha, sumRed, sumGreen, sumBlue, i
	
	Dim As Ulong windowSize = iRadius * 2 + 1, radiusPlusOne = iRadius + 1, _
                 srcIndex = 0, sumLookupTable(256 * windowSize), indexLookupTable(radiusPlusOne), dstIndex, x, y, nextPixelIndex
                
	Union _Color
		As Ulong argb
		Type
			As Ubyte b, g, r, a
		End Type
	End Union
	Dim As _Color pixel, nextPixel, previousPixel
	
	Dim As Integer pitch_img = img.pitch Shr 2, pitch_b = img_b.pitch Shr 2

	For i = 0 To Ubound(sumLookupTable) - 1
		sumLookupTable(i) = i \ windowSize
	Next
	
	If iRadius < iW Then
		For i = 0 To Ubound(indexLookupTable) - 1
			indexLookupTable(i) = i
		Next
	Else
		For i = 0 To iW - 1
			indexLookupTable(i) = i
		Next
		For i = iW To Ubound(indexLookupTable) - 1
			indexLookupTable(i) = iW1
		Next
	Endif
	
	For y = 0 To iH1
		sumAlpha = 0: sumRed = 0: sumGreen = 0: sumBlue = 0
		dstIndex = y
		
		pixel.argb = img.pixels[srcIndex]
		sumAlpha += radiusPlusOne * pixel.a
		sumRed   += radiusPlusOne * pixel.r
		sumGreen += radiusPlusOne * pixel.g
		sumBlue  += radiusPlusOne * pixel.b
		
		For i = 1 To iRadius
			pixel.argb = img.pixels[srcIndex + indexLookupTable(i)]
			sumAlpha += pixel.a			
			sumRed   += pixel.r
			sumGreen += pixel.g
			sumBlue  += pixel.b
		Next
		
		For x = 0 To iW1
			img_b.pixels[dstIndex] = sumLookupTable(sumAlpha) Shl 24 Or _
									 sumLookupTable(sumRed)   Shl 16 Or _
									 sumLookupTable(sumGreen) Shl  8 Or _
									 sumLookupTable(sumBlue)
			
			'img_b.pixels[dstIndex] = Rgba(sumLookupTable(sumRed), sumLookupTable(sumGreen), sumLookupTable(sumBlue), sumLookupTable(sumAlpha))
			
			dstIndex += pitch_b
			
			nextPixelIndex = x + radiusPlusOne
			If nextPixelIndex >= iW Then nextPixelIndex = iW1
			
			previousPixelIndex = x - iRadius
			If previousPixelIndex < 0 Then previousPixelIndex = 0
			
			nextPixel.argb = 	img.pixels[srcIndex + nextPixelIndex]
			previousPixel.argb = img.pixels[srcIndex + previousPixelIndex]
			
			sumAlpha += nextPixel.a
			sumAlpha -= previousPixel.a
			
			sumRed += nextPixel.r
			sumRed -= previousPixel.r
			
			sumGreen += nextPixel.g
			sumGreen -= previousPixel.g
			
			sumBlue += nextPixel.b
			sumBlue -= previousPixel.b
		Next
	
		srcIndex += pitch_img
	Next
	
	Return pImage_Blurred
End Function

Function FastBlur(img As tImage, iRadius As Ubyte) As Any Ptr
	'iRadius = Iif(iRadius < 0, 0, iRadius)

	Dim pImgPassH As Any Ptr = BlurPass(img, iRadius, img.width, img.height) 'horizontal pass

	Dim As tImage img2
	Imageinfo(pImgPassH, img2.Width, img2.height, , img2.pitch, img2.pixels)

	Dim pImgPassW As Any Ptr = BlurPass(img2, iRadius, img.height, img.Width) 'vertical pass
	ImageDestroy(pImgPassH)

	Return pImgPassW
End Function

Sub ImageContrast(pImage as any Pointer, contrast as Integer, brightness as Integer = 0, iW As Ushort, iH As Ushort)
   #Define Red(colour) ((colors Shr 16) And 255)
   #Define Green(colour) ((colors Shr 8) And 255)
   #Define Blue(colour) (colors And 255)
   #Define Truncate(colour) (Iif(colour < 0, 0, Iif(colour > 255, 255, CUbyte(colour))))
   Dim as Ulong colors
   Dim As Integer w, h, pitch
   Dim As Any Pointer pixdata
   Imageinfo(pImage, w, h, , pitch, pixdata)
   
   Dim as Single factor, contrastLevel = (((100.0 + contrast) / 100.0) * ((100.0 + contrast) / 100.0))

   For y as UShort = 0 to iH - 1
      For x as Ushort = 0 to iW - 1
         colors = *CPtr(ULong ptr, pixdata + y * pitch + x Shl 2)
 
         *CPtr(ULong ptr, pixdata + y * pitch + x Shl 2) = RGB(Truncate(brightness + (((Red(colour) / 255 - 0.5) * contrastLevel) + 0.5) * 255.0), _
                                                               Truncate(brightness + (((Green(colour) / 255 - 0.5) * contrastLevel) + 0.5) * 255.0), _
                                                               Truncate(brightness + (((Blue(colour) / 255 - 0.5) * contrastLevel) + 0.5) * 255.0))
      Next
   Next
End Sub

Const w = 1920 Shr 1
Const h = 1080 Shr 1
Const _t = 1 / 60

Dim As Long _s = 1, iw = w, ih = Int(iw * 9 / 16)

Screenres iw * _s, ih * _s, 32, 2, GFX_ALPHA_PRIMITIVES Or GFX_NO_SWITCH 'Or GFX_NO_FRAME 'Or GFX_ALWAYS_ON_TOP 'Or GFX_FULLSCREEN
Screenset 1, 0
Color &hFF, &hFFFFFFFF
Cls

Dim As Any Ptr pImage, pImage_blurred
pImage = Imagecreate(iw, ih, 0, 32)
Dim As tImage img
Imageinfo(pImage, img.Width, img.height, , img.pitch, img.pixels)

Randomize , 2

Dim As Ulong iFPS, cfps = 0
Dim As Double t = 0, fTimer = Timer
Dim As Long i
Type v2
	As Single r, x, y, vx, vy
End Type
Dim As v2 balls(19)
For i = 0 To Ubound(balls)
	balls(i).r = 50
	balls(i).x = balls(i).r + Rnd() * (iw - 2 * balls(i).r)
	balls(i).y = balls(i).r + Rnd() * (ih - 2 * balls(i).r)
	balls(i).vx = Rnd() * 2 - 1
	balls(i).vy = Rnd() * 2 - 1
Next

Do
	Line pImage, (0, 0) - (iw - 1, ih - 1), &hFFFFFFFF, BF
	For i = 0 To Ubound(balls)
		Circle pImage, (balls(i).x, balls(i).y), balls(i).r, &hF0000000,,,, F
		balls(i).x += balls(i).vx
		balls(i).y += balls(i).vy
		If balls(i).x < balls(i).r Or balls(i).x > iw - balls(i).r Then balls(i).vx = -balls(i).vx
		If balls(i).y < balls(i).r Or balls(i).y > ih - balls(i).r Then balls(i).vy = -balls(i).vy	
	Next
	pImage_blurred = FastBlur(img, 30)
	ImageContrast(pImage_blurred, 500, 0, iw, ih)
	Put (0, 0), pImage_blurred, Pset
	Imagedestroy(pImage_blurred)
	
	t += _t

	Draw String(4, 4), iFPS & " fps", &hFFFF0000
	
	Flip
	cfps += 1
	If Timer - fTimer > 0.99 Then
		iFPS = cfps
		cfps = 0
		fTimer = Timer
	End If
	Sleep (1)
Loop Until Len(Inkey())

Imagedestroy(pImage)
badidea
Posts: 2586
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Metaballs: blobs of goo

Post by badidea »

It think I fixed it. If you keep a button pressed, then at loopCount 260 a ball (bottom-right) starts separating.
The smaller dots are the grid centers. When white, then part of the blob.
Note: I changed and '<' to '>'. Same result, but now positive logic, easier for my mind.

Code: Select all

const N_BALL = 20
Dim shared metaballX(1 to N_BALL) as integer
Dim shared metaballY(1 to N_BALL) as integer

'randomize timer
randomize 29 '28 '27

screenres 800,600,32',2
width 800\8,600\16
'ScreenSet 1,0

for i as integer = 1 to N_BALL
	metaballX(i) = 800*rnd(1)
	metaballY(i) = 600*rnd(1)
next

function scalarfield (x as integer, y as integer) as single
	dim d as single
	for i as integer = 1 to N_BALL
		d += 500.0 / ((x-metaballX(i))^2 + (y-metaballY(i))^2)
	next
	return d
end function

const GRID_SZ = 40 '20 '10

dim as long key
dim as ulong dotColor, loopCount
do
	for i as integer = 1 to N_BALL
		metaballX(i) += 2*rnd(1)-1
		metaballY(i) += 2*rnd(1)-1
	next

	screenlock
	cls
	for y as integer = 1 to 600 step GRID_SZ
		for x as integer = 1 to 800 step GRID_SZ

			dim as single A,B,C,D, E
			A = scalarfield(x, y)
			B = scalarfield(x + GRID_SZ, y)
			C = scalarfield(x, y + GRID_SZ)
			D = scalarfield(x + GRID_SZ, y + GRID_SZ)
			E = scalarfield(x + GRID_SZ \ 2, y + GRID_SZ \ 2) 'center value
			'A--B
			'|  |
			'C--D
			dim as integer side
			if A > .5 then side or= 1
			if B > .5 then side or= 2
			if C > .5 then side or= 4
			if D > .5 then side or= 8

			dotColor = iif(A > .5, rgb(255,255,255), rgb(100,100,100))
			circle (x,y), 3, dotColor 'grey or white
			dotColor = iif(E > .5, rgb(255,255,255), rgb(100,100,100))
			circle (x + GRID_SZ \ 2, y + GRID_SZ \ 2), 1, dotColor 'grey or white

			dim as integer x1,y1,x2,y2
			select case side
			case 1, 14: '0001, 1110 (A is different)
				x1 = x + GRID_SZ *(.5-A)/(B-A)
				y1 = y + GRID_SZ *(.5-A)/(C-A)
				line(x,y1)-(x1,y), rgb(255,0,0) 'red
			case 9: '1001
				if E > .5 then
					x1 = x - GRID_SZ *(.5-B)/(A-B) + GRID_SZ 
					y1 = y + GRID_SZ *(.5-B)/(D-B)
					line(x1,y)-(x + GRID_SZ ,y1), rgb(255,255,255) 'white
					x1 = x + GRID_SZ *(.5-C)/(D-C)
					y1 = y - GRID_SZ *(.5-C)/(A-C) + GRID_SZ 
					line(x,y1)-(x1,y + GRID_SZ ), rgb(255,255,255) 'white
				else
					x1 = x + GRID_SZ *(.5-A)/(B-A)
					y1 = y + GRID_SZ *(.5-A)/(C-A)
					line(x,y1)-(x1,y), rgb(255,255,255) 'white
					x1 = x - GRID_SZ *(.5-D)/(C-D) + GRID_SZ 
					y1 = y - GRID_SZ *(.5-D)/(B-D) + GRID_SZ 
					line(x + GRID_SZ ,y1)-(x1,y + GRID_SZ ), rgb(255,255,255) 'white
				end if
			case 2, 13: '0010, 1101 (B is different)
				x1 = x - GRID_SZ *(.5-B)/(A-B) + GRID_SZ 
				y1 = y + GRID_SZ *(.5-B)/(D-B)
				line(x1,y)-(x + GRID_SZ ,y1), rgb(255,255,0) 'yellow
			case 6: '0110
				if E > .5 then
					x1 = x + GRID_SZ *(.5-A)/(B-A)
					y1 = y + GRID_SZ *(.5-A)/(C-A)
					line(x,y1)-(x1,y), rgb(255,255,255) 'white
					x1 = x - GRID_SZ *(.5-D)/(C-D) + GRID_SZ 
					y1 = y - GRID_SZ *(.5-D)/(B-D) + GRID_SZ 
					line(x + GRID_SZ ,y1)-(x1,y + GRID_SZ ), rgb(255,255,255) 'white
				else
					x1 = x - GRID_SZ *(.5-B)/(A-B) + GRID_SZ 
					y1 = y + GRID_SZ *(.5-B)/(D-B)
					line(x1,y)-(x + GRID_SZ ,y1), rgb(255,255,255) 'white
					x1 = x + GRID_SZ *(.5-C)/(D-C)
					y1 = y - GRID_SZ *(.5-C)/(A-C) + GRID_SZ 
					line(x,y1)-(x1,y + GRID_SZ ), rgb(255,255,255) 'white
				end if
			case 3, 12: '0011, 1100 (horizontal)
				y1 = y + GRID_SZ *(.5-A)/(C-A)
				y2 = y + GRID_SZ *(.5-B)/(D-B)
				line(x,y1)-(x + GRID_SZ ,y2), rgb(0,255,0) 'green
			case 4, 11: '0100, 1011 (C is different)
				x1 = x + GRID_SZ *(.5-C)/(D-C)
				y1 = y - GRID_SZ *(.5-C)/(A-C) + GRID_SZ 
				line(x,y1)-(x1,y + GRID_SZ ), rgb(0,255,255) 'cyan
			case 5, 10: '0101, 1010 (vertical)
				x1 = x + GRID_SZ *(.5-A)/(B-A)
				x2 = x + GRID_SZ *(.5-C)/(D-C)
				line(x1,y)-(x2,y + GRID_SZ ), rgb(0,0,255) 'blue
			case 7, 8: '0111, 1000 (D is different)
				x1 = x - GRID_SZ *(.5-D)/(C-D) + GRID_SZ 
				y1 = y - GRID_SZ *(.5-D)/(B-D) + GRID_SZ 
				line(x + GRID_SZ ,y1)-(x1,y + GRID_SZ ), rgb(255,0,255) 'pink
			end select
		next 
	next
	draw string (10,5), "Print <escape> to exit, <any key> to change positions"
	draw string (10,25), "loopCount: " & loopCount
	screenunlock

	'flip 1,0
	key = getkey
	loopCount += 1

loop while key <> 27
print "End."
Case 6 & 9 are each opposite, so more efficient rewrite to be done...

@UEZ: Also looks cool, but blur is not what I want.
UEZ
Posts: 972
Joined: May 05, 2017 19:59
Location: Germany

Re: Metaballs: blobs of goo

Post by UEZ »

badidea wrote: Apr 07, 2022 22:34 @UEZ: Also looks cool, but blur is not what I want.
I know what you want - not "fake" meta balls but real calculated ones. :wink:
badidea
Posts: 2586
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Metaballs: blobs of goo

Post by badidea »

UEZ wrote: Apr 08, 2022 6:02 I know what you want - not "fake" meta balls but real calculated ones. :wink:
Actually, I do not know what metaballs are. Wikipedia says: "metaballs are organic-looking n-dimensional isosurfaces, characterised by their ability to meld together when in close proximity to create single, contiguous objects". Also: "Not to be confused with Meatballs."
But I am anyway mostly interested in the contour lines that one can make with the marching squares. See modified code:

Code: Select all

const N_BALL = 20
Dim shared metaballX(1 to N_BALL) as integer
Dim shared metaballY(1 to N_BALL) as integer

'randomize timer
randomize 29 '28 '27

screenres 800, 600, 32',2
width 800\8, 600\16
'ScreenSet 1,0

for i as integer = 1 to N_BALL
	metaballX(i) = 800 * rnd(1)
	metaballY(i) = 600 * rnd(1)
next

function scalarfield (x as integer, y as integer) as single
	dim d as single
	for i as integer = 1 to N_BALL
		d += 500.0 / ((x-metaballX(i))^2 + (y-metaballY(i))^2)
	next
	return d
end function

const GRID_SZ = 10

dim as ulong dotColor
dim as single cl1 = 0.1 'contour line
do
	cl1 *= 1.2
	if cl1 > 4 then exit do
	for y as integer = 1 to 600 step GRID_SZ
		for x as integer = 1 to 800 step GRID_SZ

			dim as single A, B, C, D, E
			A = scalarfield(x, y)
			B = scalarfield(x + GRID_SZ, y)
			C = scalarfield(x, y + GRID_SZ)
			D = scalarfield(x + GRID_SZ, y + GRID_SZ)
			E = scalarfield(x + GRID_SZ \ 2, y + GRID_SZ \ 2) 'center value
			'A--B
			'|  |
			'C--D
			dim as integer side
			if A > cl1 then side or= 1
			if B > cl1 then side or= 2
			if C > cl1 then side or= 4
			if D > cl1 then side or= 8

			'draw grid
			dotColor = iif(A > cl1, rgb(255,255,255), rgb(100,100,100))
			pset(x,y), dotColor 'grey or white
			dotColor = iif(E > cl1, rgb(255,255,255), rgb(100,100,100))
			pset(x + GRID_SZ \ 2, y + GRID_SZ \ 2), dotColor 'grey or white

			dim as integer x1,y1,x2,y2
			select case side
			case 1, 14: '0001, 1110 (A is different)
				x1 = x + GRID_SZ * (cl1-A) / (B-A)
				y1 = y + GRID_SZ * (cl1-A) / (C-A)
				line(x, y1)-(x1, y), rgb(255,0,0) 'red
			case 2, 13: '0010, 1101 (B is different)
				x1 = x + GRID_SZ * (cl1-A) / (B-A)
				y1 = y + GRID_SZ * (cl1-B) / (D-B)
				line(x1, y)-(x + GRID_SZ, y1), rgb(255,255,0) 'yellow
			case 4, 11: '0100, 1011 (C is different)
				x1 = x + GRID_SZ * (cl1-C) / (D-C)
				y1 = y + GRID_SZ * (cl1-A) / (C-A)
				line(x, y1)-(x1, y + GRID_SZ), rgb(0,255,255) 'cyan
			case 7, 8: '0111, 1000 (D is different)
				x1 = x + GRID_SZ * (cl1-C) / (D-C)
				y1 = y + GRID_SZ * (cl1-B) / (D-B)
				line(x + GRID_SZ, y1)-(x1, y + GRID_SZ), rgb(255,0,255) 'pink
			case 3, 12: '0011, 1100 (horizontal)
				y1 = y + GRID_SZ * (cl1-A) / (C-A)
				y2 = y + GRID_SZ * (cl1-B) / (D-B)
				line(x, y1)-(x + GRID_SZ, y2), rgb(0,255,0) 'green
			case 5, 10: '0101, 1010 (vertical)
				x1 = x + GRID_SZ * (cl1-A) / (B-A)
				x2 = x + GRID_SZ * (cl1-C) / (D-C)
				line(x1, y)-(x2, y + GRID_SZ), rgb(0,0,255) 'blue
			case 6, 9: '0110, 1001
				if (side = 6 and E > cl1) or (side = 9 and E < cl1) then
					x1 = x + GRID_SZ * (cl1-A) / (B-A)
					y1 = y + GRID_SZ * (cl1-A) / (C-A)
					line(x, y1)-(x1, y), rgb(255,255,255) 'white
					x1 = x + GRID_SZ * (cl1-C) / (D-C)
					y1 = y + GRID_SZ * (cl1-B) / (D-B)
					line(x + GRID_SZ, y1)-(x1, y + GRID_SZ ), rgb(255,255,255) 'white
				else
					x1 = x + GRID_SZ * (cl1-A) / (B-A)
					y1 = y + GRID_SZ * (cl1-B) / (D-B)
					line(x1,y)-(x + GRID_SZ ,y1), rgb(255,255,255) 'white
					x1 = x + GRID_SZ * (cl1-C) / (D-C)
					y1 = y + GRID_SZ * (cl1-A) / (C-A)
					line(x,y1)-(x1,y + GRID_SZ ), rgb(255,255,255) 'white
				end if
			end select
		next 
	next
	sleep 400
loop
print "End."
getkey
The metaballs from angros47 are more like infinity pillars.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Metaballs: blobs of goo

Post by dodicat »

Some non metallic blobs.

Code: Select all

'#cmdline "-exx" 'checked
Type pt
    As Single x,y
    As Long bound
End Type

Type vpt Extends pt
    As Single dx,dy
    Declare Sub set(xx As Single,yy As Single,dxx As Single,dyy As Single)
End Type

Sub vpt.set(xx As Single,yy As Single,dxx As Single,dyy As Single)
    x=xx:y=yy:dx=dxx:dy=dyy
End Sub

Type blob
    As pt a(Any)
    As Long number
    As Long size
End Type

Sub GetClosest(a() As pt, ans As blob,v As pt)
    #define incircle(cx,cy,radius,x,y) (cx-x)*(cx-x) +(cy-y)*(cy-y)<= radius*radius
    Dim As Single r=ans.size,ctr,rad=ans.size
    Do
        r+=.5
        ctr=0
        For n As Long=Lbound(a) To Ubound(a)
            If incircle(v.x,v.y,r,a(n).x,a(n).y) And incircle(v.x,v.y,rad,a(n).x,a(n).y)=0 Then
                ctr+=1
                Redim Preserve ans.a(1 To ctr)
                ans.a(ctr)=a(n)
            End If
        Next n
    Loop Until Ubound(ans.a)>=ans.number
End Sub

Sub circulate(p As blob)
    #macro Circlesort() 
    For p1 As Long  = Lbound(p.a) To Ubound(p.a)-1
        For p2 As Long  = p1 + 1 To Ubound(p.a)
            If Atan2(p.a(p1).y-c.y,p.a(p1).x-c.x)< Atan2(p.a(p2).y-c.y,p.a(p2).x-c.x) Then
                Swap p.a(p1),p.a(p2)
            End If
        Next p2
    Next p1
    #endmacro
    Dim As pt C '--centroid of points
    Dim As Long counter
    For n As Long=Lbound(p.a) To Ubound(p.a)
        counter+=1
        c.x+=p.a(n).x
        c.y+=p.a(n).y
    Next n
    c.x=c.x/counter
    c.y=c.y/counter
    CircleSort()
End Sub

Sub move(v As vpt)
    Dim As Long xres,yres
    Screeninfo xres,yres
    v.x+=v.dx
    v.y+=v.dy
    If v.x<50 Then v.x=50:v.dx=-v.dx
    If  v.x>xres-50 Then v.x=xres-50:v.dx=-v.dx
    If v.y<50 Then v.y=50:v.dy=-v.dy
    If v.y>yres-50 Then v.y=yres-50:v.dy=-v.dy
End Sub

Function inpolygon(p1() As pt, p2 As pt) As Long
    #define Winder(L1,L2,p) ((L1.x-L2.x)*(p.y-L2.y)-(p.x-L2.x)*(L1.y-L2.y))
    Dim As Long index,nextindex,k=Ubound(p1)+1,wn
    For n As Long=1 To Ubound(p1)
        index=n Mod k:nextindex=(n+1) Mod k
        If nextindex=0 Then nextindex=1
        If p1(index).y<=p2.y Then
            If p1(nextindex).y>p2.y Andalso  Winder(p1(index),p1(nextindex),p2)>0 Then wn+=1
        Else
            If p1(nextindex).y<=p2.y Andalso Winder(p1(index),p1(nextindex),p2)<0 Then wn-=1
        End If
    Next n
    Return wn
End Function

Sub bline(x1 As Long,y1 As Long,x2 As Long,y2 As Long,col As Ulong, p() As pt,a() As pt,Byref count As Long)
    Var dx=Abs(x2-x1),dy=Abs(y2-y1),sx=Sgn(x2-x1),sy=Sgn(y2-y1)
    Dim As Long e
    If dx<dy Then  e=dx\2 Else e=dy\2
    Dim As pt tmp
    Do
        count+=1
        a(count).x=x1
        a(count).y=y1
        If inpolygon(p(),a(count)) Then a(count).bound=1
        If x1 = x2 Then If y1 = y2 Then Exit Do
        If dx > dy Then
            x1 += sx : e -= dy : If e < 0 Then e += dx : y1 += sy
        Else
            y1 += sy : e -= dx : If e < 0 Then e += dy : x1 += sx
        End If
    Loop
End Sub

Function GetRims(c() As blob,p() As pt,a() As pt) As Long
    Dim As Long n,count
    For n1 As Long=Lbound(c) To Ubound(c)-1
        For n2 As Long=n1+1 To Ubound(c)
            For n =Lbound(c(n1).a) To Ubound(c(n1).a)-1
                bline(c(n1).a(n).x,c(n1).a(n).y,c(n1).a(n+1).x,c(n1).a(n+1).y,Rgb(100,0,0),c(n2).a(),a(),count)
            Next n
            bline(c(n1).a(n).x,c(n1).a(n).y,c(n1).a(1).x,c(n1).a(1).y,Rgb(100,0,0),c(n2).a(),a(),count)
            For n =Lbound(c(n2).a) To Ubound(c(n2).a)-1
                bline(c(n2).a(n).x,c(n2).a(n).y,c(n2).a(n+1).x,c(n2).a(n+1).y,Rgb(100,0,0),c(n1).a(),a(),count)
            Next n
            bline(c(n2).a(n).x,c(n2).a(n).y,c(n2).a(1).x,c(n2).a(1).y,Rgb(100,0,0),c(n1).a(),a(),count) 
        Next n2
    Next n1
    Redim Preserve a(1 To count)
    Return 0
End Function

Sub init(p() As pt,st As Long=15)
    Dim As Long xres,yres
    Screeninfo xres,yres
    Dim As Long count
    Var k=-20  '+300
    For x As Long=-k To xres+k Step st
        For y As Long=-k To yres+k Step st
            count+=1
            Next: Next
            Redim p(1 To count)
            count=0
    For x As Long=-k To xres+k Step st
        For y As Long=-k To yres+k Step st
    count+=1
    p(count)=Type(x,y)
    Next: Next  
End Sub 

Function Regulate(Byval MyFps As Long,Byref fps As Long) As Long
    Static As Double timervalue,lastsleeptime,t3,frames
    frames+=1
    If (Timer-t3)>=1 Then t3=Timer:fps=frames:frames=0
    Var sleeptime=lastsleeptime+((1/myfps)-Timer+timervalue)*1000
    If sleeptime<1 Then sleeptime=1
    lastsleeptime=sleeptime
    timervalue=Timer
    Return sleeptime
End Function
                
    'set up the arrays
    #define lim 5
    Redim As pt p()
    Redim As blob shape(1 To lim)
    shape(1).number=15:shape(1).size=60
    shape(2).number=30:shape(2).size=120
    shape(3).number=25:shape(3).size=100
    shape(4).number=25:shape(4).size=100
    shape(5).number=30:shape(5).size=120
                
    Dim As vpt v(1 To lim)
    v(1).set(500,300,1.15,1.12)
    v(2).set(200,300,1.15,-1.12)
    v(3).set(400,300,1.1,1.2)
    v(4).set(10,10,1.2,1.2)
    v(5).set(700,100,1.2,1.2)
    
     Screen 20,32
     Dim As Ulong background=Rgb(255,255,255)
     Color ,background
Dim As Any Ptr i=Imagecreate(1024,768)

    init(p())
    
     Line i,(0,0)-(1023,767),Rgb(0,100,255),b           
    For n As Long=1 To Ubound(p)        'get grid points to image
        Circle i,(p(n).x,p(n).y),1,Rgb(100,100,100),,,,f
    Next n
    
    Dim As Long fps            
Do
       Redim As pt a(1 To 100000)
                    
        For n As Long=1 To Ubound(v)
            move(v(n))
        Next
                    
        For n As Long=1 To Ubound(v)
            GetClosest(p(),shape(n),v(n))
        Next
                    
        For n As Long=1 To Ubound(shape)
            circulate(shape(n))
        Next n
                    
        GetRims(shape(),p(),a())
        
        Screenlock
        Cls
        For n As Long=1 To Ubound(a)
            Circle(a(n).x,a(n).y),1,Rgb(0,0,0),,,,f
        Next
                    
        For n As Long=1 To Ubound(a)
            If a(n).bound=1  Then Circle (a(n).x,a(n).y),1,background,,,,f
        Next
                    
        Put(0,0),i,trans
        Draw String (50,50),"Framerate = " &fps,Rgb(0,0,0)
        Draw String (50,70),"Press <esc> to end ",Rgb(0,0,0)
        Screenunlock
        Sleep regulate(40,fps)                    
Loop Until Inkey=Chr(27)
                
        Sleep
        Imagedestroy i
                 
Post Reply