Metaballs script

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.

Hi! I'm really Interested in this.
Will you make it compatible with octane render?

Hello Spartan; I'm not sure how Octane works, so it's difficult to predict.
For rendering a single frame I'd say most probably yes, because that field is just a regular prop with regular texture nodes, which is regenerated on each frame.
For rendering a movie that will have to be seen, as the prop has variable geometry, so it has to be fully regenerated on each frame. That's done through a Python callback. So if Octane refreshes everything for each frame, then it should work, but if it retains information from the previous frame then it might not work.
One thing I'm thinking, though, to avoid the trouble with the callback, is to build the field as a series of props, and just hide/make them visible in the right frames. That way it should be compatible with any renderers.