Metaballs script

I'm moving a discussion on a metaball script from the another thread with a generic name to here.
This is today's test of the thingie. First part has the metaballs movig and changing strength while the field keeps its resolution and threshold; then the middle part is the field changing it's resolution (I made that animated too, so that parts that require more resolution can be set with a finer resolution); last part has the field changing threshold (which is just a global scale factor for the metaball strengths):

So, what has been achieved up so far with this thingie:
 Basic field calculations for f(r) = str/r^2 (which is a bad formula, but works for now)
 Field and mesh calcs are based on volume crawling instead of multiple loops
 Meshes are now welded and have smooth shading
 tranX/Y/Z, strength, threshold and resolution are animated
 Baking for a range is set, and the meshes are saved in memory
 Playback runs at full speed, as the meshes are already baked
What I'm planning to do next is in these directions, and the are easy to do:
 Add a frame frame as a dial
 Add scalingX/Y/Z and rotationX/Y/Z for metaellipsis instead of metaballs
 Add metacubes with scaling and rotation
I'm also planning this, but somewhat fuzzy on how to do that yet:
 Add a metacone
 Add a metacylinder
 I'm just deleting and recreating the prop on every frame; that goes surprisingly fast, but I lose the materials on every frame. Gotta find a better way of changing the geometry during playback.
Finally, this is last thing, but I have no idea on how to do it, as the API for materials seems to be quite out of whack in the Python API:
 Copy the render nodes from the isofield prop to the isosurface during playback
How the heck there is a GetShaderTree() but there isn't a SetShaderTree() function in the API? Looks like SM wants people to go the looooooooooooong way on copying materials around. This will kill my plaback speed!

It's snowing outside, so I'm stuck inside, therefore back to Metaballs. I don't think these metaballs are very sticky. I think the problem is the field function. I'm using
f(r) = str/r^2
so the isosurface of a single metaball is defined by
str/r^2 = thld, or r = sqrt(str/thld)
so if I put two metaballs of the same str apart by d, then they will stick if
str/(d/2)^2+str(d/2)^2 = thld, or d = sqrt(8*str/thld); so d/r = sqrt(8) = 2.83
but each metaball has a radius of r/r=1.00, so each metaball willl only stick at 40% of its own radius, or 20% of its diameter. I think that's why the balls don't look very sticky. Problem is indeed the field function.
Say I want the metaballs to stick at 2x their own radius, and at the same time I need to use r^2 to avoid calculating a square root. And then I can get rid of that division too if I invert it. So, say u=r^2/(str*thld)^2; then I need a function f(u) such that f(0) > 1, f(1) = 1, f(4) = 0.5 and f(8) = 0.
That way I can set threshold = 1.0 and it becomes just a global field scale, and then the metaball will have a radius = strength; can just rename str to r0 to make it more clear, and will be consistent with the tranX/Y/Z units. The radius will be r0, and it will stick at 2*r0.
Now, say I added another parameter "Stickness", call is K, that controls that. So I need a function
f(0) > 1, f(1) = 1, f(K) = 0.5, f(2*K) = 0
3 points define a parabola, so
f(u) = (ua)(ub)/c
=> (1x)(1y)/z = 1, (Kx)(Ky)/z = 0.5, (2Kx)(2Ky)/z = 0Save time with wolframalpha for linear algebra (haha), and get
a = 2 K
b = 2 K^2  2 K + 1
c = 2 (2 K^3  3 K^2 + K)so for K=2 then a=4, b=5, c=12, or f(r) = (1/12)(r4)(r5), which indeed goes f(1) = 1, f(2) = 1/2, f(4) = 0. Easy peasy.
So let's see how this guy goes:
u = r^2/(strthld)^2
f(r) = (1/c)(ra)(rb) if r < (strthrld)(4K^2)
f(r) = 0 if r < (strthrld)*(4K^2)

This is the result of the field function
f(u,K) = (ua)(ub)/c
It's so and so. That test has two balls with stickness = 2.0 in the first part, which looks good until then  they do get sticky further apart than before; but then the 2nd part has stickness increasing in both balls from 2.0 to 5.0  that makes them stick even if they don't move and don't change their strengths. That worked ok, but the link between them gets too thick too fast; may be a range thing, but I wanted a smoother effect.
Also the radius seems to be getting affected by K; I messed up something to cause that. Gotta check that field function again...


@fbs7 said in Metaballs script:
This is the result of the field function
f(u,K) = (ua)(ub)/c
It's so and so. That test has two balls with stickness = 2.0 in the first part, which looks good until then  they do get sticky further apart than before; but then the 2nd part has stickness increasing in both balls from 2.0 to 5.0  that makes them stick even if they don't move and don't change their strengths. That worked ok, but the link between them gets too thick too fast; may be a range thing, but I wanted a smoother effect.
Also the radius seems to be getting affected by K; I messed up something to cause that. Gotta check that field function again...
This looks pretty good!
Looking forward how you progress with this!

Thank you. By the way, I can't speak highly enough of WolframAlpha; if anyone has any kind of algebraic problem, the thing is really a miraculous helper.
I've found other simple solutions for a "sticky" field equation. One idea is this:
 Still quadratic form f(r,K) = (ar)(br)/c such that f(1,K) = 1, f(K,K) = 0.5, a=2K so that field goes to 0 at 2K:
Solution:
a = 2 K, b = 2 K^2  2 K + 1, c = 2 (2 K^3  3 K^2 + K) This is from http://www.geisswerks.com/ryan/BLOBS/blobs.html:
f(r) = (r^4  r^2 + 0.25)
also
f(r) = 1(6r^5  15r^4 + 10r^3)
 Wikipedia has
f(r) = (1r^2)^2

Houdini lets you choose sum/max/min form for f(r) (nice!!!), and also the field function

Pixar's RenderMan uses a quite popular form:
f(r) = 1  3R^2 + 3R^4  R^6
 Blinn (first guy that proposed metaballs) used:
f(r) = a.exp(br^2)
...yikies, not touching that Maya uses a cubic form, although I didn't write it down

Hmm... I see now the problem with that stickness thing. It's not due to the field equation. Everybody uses a similar polynomial form. The problem is that the field is symmetric; so if I want a small decrease in the field f(r,theta) for a small increase in r, then another point with quite a different theta will have a similar field, so if another metaball comes close the thing will bulge on several different angles. All symmetric field equations will give the same result.
So I could add a field linking all metaball centers, but that's exactly the same thing as the user adding a metacylinder between the balls. So I'm leaving that alone, at least for now.

So, moving on. Added scaling in the code. Scale is pretty trivial; it just follows the scale of the handlers. These are the handlers before the simulation:
And this is the simulation:

Alright, now rotations. I was never good at understanding rotation matrices in college, so any help will be helpful. I have a field prop, say f, which may be parented to something and transformed per the user's wishes. That guy will have a world matrix Wf.
Now I have a metaball, say m, which is parented to f and transformed, and it has a local matrix Lf.
Now I have a coordinate XYZ; that's the coordinate in the field prop; that's where I make the calculations, so let's call it XYZf. So in order to get this coordinate relative to m, say XYZm, then will I use this?
XYZm = inv(Lf)*XYZf
Is that right? Or should I use the direct transformation XYZm = (Lf)*XYZf?
Now that I have the field potential at XYZm, I just assign that as the field potential at XYZf, as both refer to the same point.
Then I build the isosurface as a series of coordinates XYZi. These are of course relative to "f". But I'm building the prop with those coordinates in world space (CreatePropFromGeometry only works in world space, I think). Now I want to set my isosurface as a child of "f". The moment I do that, I think Poser will transform my coordinates, because I see that happening all the time in the GUI. So Poser things that XYZ is a world coordinate. I think it wil do this:
XYZf = (inv Wf) * XYZi
But that will break my surface, as it's already in "f" coordinates. So I think I gotta do build the isosurface prop as
XYZii = Wf * XYZi
and give that to Poser to CreatePropFromGeometry() and SetParent(), so that when Poser transforms it will come back to local coordinates.
Or should I just SetParent() with realign = 0? Will that stop Poser from applying (inv Wf) to my prop?

Hmm... rewriting to correct a couple of typos:
Alright, now rotations. I was never good at understanding rotation matrices in college, so any help will be helpful. I have a field prop, say f, which may be parented to something and transformed per the user's wishes. That guy will have a world matrix Wf.
Now I have a metaball, say m, which is parented to f and transformed, and it has a local matrix Lm.
Now I have a coordinate XYZ; that's the coordinate in the field prop; that's where I make the calculations, so let's call it XYZf. So in order to get this coordinate relative to m, say XYZm, then will I use this?
XYZm = inv(Lm)*XYZf
Is that right? Or should I use the direct transformation XYZm = (Lf)*XYZf?
Now that I have the field potential at XYZm, I just assign that as the field potential at XYZf, as both refer to the same point.
Then I build the isosurface as a series of coordinates XYZi. These are relative to "f". But I'm building the prop with those coordinates in world space (CreatePropFromGeometry only works in world space, I think). Now I want to set my isosurface as a child of "f". The moment I do that, I think Poser will transform my coordinates, because I see that happening all the time in the GUI. So Poser things that XYZ is a world coordinate. I think it wil do this:
XYZf = (inv Wf) * XYZi
But that will break my surface, as it's already in "f" coordinates. So I think I gotta do build the isosurface prop as
XYZii = Wf * XYZi
and give that to Poser to CreatePropFromGeometry() and SetParent(), so that when Poser transforms it will come back to local coordinates.
Or should I just SetParent() with realign = 0? Will that stop Poser from transforming my prop upon parenting?

Alright, progress is slow when you don't know what you're doing. I got hold of these transformation matrices now, and also resolved the problem about changing geometry of the prop on the fly without having to rebuild it on each frame  that wasn't working with transformations involved, as Poser was flashing the prop before and after parenting during recreation, so that wouldn't work for preview.
Fix was to indeed to use SetParent( actor, 0, 0 ) to avoid transformation changes during the first parenting, then just use SetGeometry() / MarkGeomChanged() to change geometry on the fly. As the materials are no longer destroyed, they are kept during playback. This is tonight's test:

Actually I lie. The actual fix was, after 4 hours of trying, was this sequence:
oldParent = isoSurface.Parent()
isoSurface.SetParent( scene.Actor('UNIVERSE') )
isoSurface.SetGeometry( ...whatever geometry you want... )
isoSurface.MarkGeomChanged()
isoSurface.SetParent( oldParent )That allows playback at full speed, keeps the materials and allows the metaball field to be rotated/scaled/moved at will.

@fbs7 virtually dying to play with this!!! Holler if you need a Mac beta tester! 8D

Thank you; as soon as I get the local transformation in I'll post the updated versions.

Alright, after a lot of frustration with these local transformations (SM, please indicate in the manual that Actor.LocalMatrix() returns the transformation matrix as the transposed of what everybody else uses... ), finally got the metaball local rotations and scaling under control.
This is tonight's test:

Another fast test:

Okay, so where to go from here. Thinking of next things to tackle:

Add metacylinder (field equation seems straightforward), metacube (doesn't seem tough either), metacone (no idea about field equation for this guy), metathorus (another one with odd field equation)

Find a way to copy the whole material tree from the field prop (the big wireframe thingie) to the surface

Find a way to extend the surface beyond the boundaries of the field; now that the calculation is a volume crawler it should be possible, by replacing the field array lattice[nX,nY,nZ] with a dictionary lattice(nX,nY,nZ); not sure about performance with that  dictionaries are much slower than matrices. Another thought is to keep resizing the field lattice whenever the isosurface touches the boundaries... hmmm...

Change from a volume crawler to a surface crawler; this might expedite processing and avoid lengthy calculations of internal volumes for large metaballs, but it does create a problem about finding the proper seed points for the crawler

Add a GUI dialog to Add Metaball and Run Simulation. I guess this should be easy, I think.

Add a callback to recalculate the isosurface at the current frame whenever the user changes one of the metaballs. That might make it easier to see the effects of manipulation in real time.
Plus I have a random Poser crash that happens during playback. This is infuriating. The thing only blows up, and Poser doesn't even have the decency of telling why it crashed (like "Invalid geometry", or "Update failed" or whatever). How can I fix something if the stupid darn idiotic damned mongrel of a thing doesn't give me any debugging information!!


Hmm... actually I do have a problem with my crawler, that I didn't think of before. The volume crawler starts from the center of the metaballs and keeps crawling the volume while the field value is > threshold. I thought that would work with negative metaballs because the seed on the negative metaball would fail, but then the seed on the positive metaball would crawl its volume entirely.
But, here's the problem:
A positive metaball in green, and a negative metaball in red that makes the field potential at the center of the positive metaball to be < threshold. In this case the crawler will stop instead of calculating a volume in the shaded blue (which would make a thoruslike kind of thingie), because it will think it's outside the metaball.... hmm... how to fix this... hmmm....
I think that each lattice point I should separate the sum of the positive contributions and the sum of the negative contributions. Then if the sum of the positive contributions only is > threshold, then add the neighboring points to the crawler.
So, a volume crawler with that fix would be like this:
for all metaballs push center into volume crawler list while volume crawler list is not empty pop point for all metaballs f1 = sum of all positive metaball fields at that point f2 = sum of all negative metaball fields at that point lattice[point] = f1+f2 if f1 > 1.0 push all 26 neighbor points into the volume crawler list push all 8 west & south neighbor points into the triangle generation list
It's just 2 lines change from what I'm doing right now. But how to change that from a volume crawler to a surface crawler when there are negative metaballs? hmm... not easy problem...

During my research I came across the paper that proposed the marching cubes thing in 1987 (Lorensen & Cline):
http://academy.cba.mit.edu/classes/scanning_printing/MarchingCubes.pdf
It's pretty cool. Everybody just copies the same algorithm since then. Well done to have an algorithm copied to the letter for 30 years.
Now, the surface crawler thing is easy when all metaballs are positive (just start from the center of the metaball, then choose a direction like X and keep going until it flips from > threshold to < threshold; then crawl around the neighbor cubes that some vertices are > threshold and some are < threshold.
But when you have negative metaballs can't do that anymore, as the surfaces may be completely split. I can't think of a general way to handle that. The volume crawler with the separate f1 and f2 functions as above seems to be the only general way that I see, and should handle solids with holes and also complete splits.
So, not spending more time with the surface crawler thingie.