# Metaballs script

• Thank you. By the way, I can't speak highly enough of Wolfram-Alpha; 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) = (a-r)(b-r)/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)

f(r) = (r^4 - r^2 + 0.25)

also

f(r) = 1-(6r^5 - 15r^4 + 10r^3)

• Wikipedia has

f(r) = (1-r^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 re-creation, 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! 8-D

• 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 thorus-like 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):

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.

• What a pain in the rear is SeoGeometry() and Python callbacks!!! The thing will crash randomly on you with no information of what's wrong, and the manual provides minimum insights on any of the callbacks. For example, it has this very useless description of a typical callback:

SetVertexCallback
Explanation
Set a per-update callback to process actor (vertices) while in local deformed
space. The user-defined function will be called back on each full update of an
actor and allows for manipulation of the vertices (or other information) at the
point in the display pipeline of Poser which processes vertex level deformations.
Arguments
The callback function should take an actor object and make desired
modifications to it.

It refers to a "display pipeline" but gives no information of how that pipeline is, so I have no clue where this actually fits in the flow. It tells I can make "desired modifications" on "other information" but doesn't say what I can change - can I change textures? geometry? only vertex coordinates? Also, what's a "full update"? If the frame changes but the actor is static will it update?

The Python API is incredibly vague and limited about geometry changes! What a shame! It seems that whoever wrote it thought paper was expensive, so he/she limited the description to the smallest amount possible.

So because someone didnt' want to spend 15 minutes describing an API function fully, I (and probably other customers) had to spend 15 hours of my life trying to identify what I can and cannot do inside a callback.

I'm giving up on this sleazy function. Mission now is to dynamically change the geometry of a prop without using that loser of a function which is SetGeometry().

• @fbs7 would you like me to see whether the Mac version is less prone to the crashes you are experiencing? I'm quite happy to bulk the think up with debugging statements, and I have unbuffered logging routines that write to an external file which can be externally monitored in real time, so Poser crashing won't prevent any of the previously issued logging being swallowed by the output buffering.
I could, alternatively, send the logging routines to you, if you prefer.

• @fbs7 have you made progress on transferring material trees? I have developed Python routines which can save and load pose, character and scene files to and from the library, bypassing Poser's own load and save routines, but have yet to implement the material tree creation. It's been a while since I last looked at it, but IIRC I found something missing from the Python materials API which meant I couldn't just create all the nodes and trees. I'll have another look, but the simplest solution for transferring materials may end up being just save from the field prop and reload them from the library onto the metaball prop(s).

• @anomalaus said in Metaballs script:

@fbs7 would you like me to see whether the Mac version is less prone to the crashes you are experiencing? I'm quite happy to bulk the think up with debugging statements, and I have unbuffered logging routines that write to an external file which can be externally monitored in real time, so Poser crashing won't prevent any of the previously issued logging being swallowed by the output buffering.

That's a good idea. Try these two scripts, if you would. This one is "Create Metaball.py":

``````import poser
from numpy import oldnumeric

#
# This script creates a cube to limit an isosurface, and a handler for a metaball
# fbs 12/7/2017
#

# Check if 'IsoField' already exists, if not create one
scene = poser.Scene()
try:
isoField = scene.Actor('IsoField')
except:
# Vertices of a cube
c = 1.0         # Cube size
verts = oldnumeric.array( [[0,0,0],[0,0,c],[0,c,0],[0,c,c],[c,0,0],[c,0,c],[c,c,0],[c,c,c]] )

# Sets = vertex indexes
sets = numpy.array( [0,1,3,2,  0,4,5,1,  0,2,6,4,  7,5,4,6,  7,3,1,5,  7,6,2,3] )

# Polys = faces, index to sets
polys = oldnumeric.array( [[0,4],[4,4],[8,4],[12,4],[16,4],[20,4]] )

# Create the geometry
meshGeom = poser.NewGeometry()

# Create the prop
isoField = scene.CreatePropFromGeom(meshGeom,'IsoField')

# Set draw style to wireframe
isoField.SetDisplayStyle( poser.kDisplayCodeWIREFRAME )

# End Frame
parmEndFrame = isoField.CreateValueParameter('End Frame')
parmEndFrame.SetValue(30)
parmEndFrame.SetMinValue(1)
parmEndFrame.SetSensitivity(1)
parmEndFrame.SetForceLimits(1)

# Start Frame
parmStartFrame = isoField.CreateValueParameter('Start Frame')
parmStartFrame.SetValue(1)
parmStartFrame.SetMinValue(1)
parmStartFrame.SetSensitivity(1)
parmStartFrame.SetForceLimits(1)

# Resolution
parmResolution = isoField.CreateValueParameter('Resolution')
parmResolution.SetValue(10)
parmResolution.SetMinValue(1)
parmResolution.SetMaxValue(20)
parmResolution.SetSensitivity(1)
parmResolution.SetForceLimits(1)

# Threshold
parmThreshold = isoField.CreateValueParameter('Threshold')
parmThreshold.SetValue(1)
scene.DrawAll()

# Find the smallest N that 'MetaBallN' does not already exist as children of the iso field
childList = isoField.Children()
found = 0
for n in xrange(1,101):
metaBallName = u'MetaBall' + str(n)	# Poser needs Unicode
found = 0
for child in childList:
if child.Name() == metaBallName:
found = 1
break
if found == 0:
break

# Now let's create a metaball handler; it will be a little cube; code is the same as above
c = 0.01
verts = oldnumeric.array( [[0,0,0],[0,0,c],[0,c,0],[0,c,c],[c,0,0],[c,0,c],[c,c,0],[c,c,c]] ) - c/2
sets = oldnumeric.array( [0,1,3,2,  0,4,5,1,  0,2,6,4,  7,5,4,6,  7,3,1,5,  7,6,2,3] )
polys = oldnumeric.array( [[0,4],[4,4],[8,4],[12,4],[16,4],[20,4]] )

# Create the geometry
meshGeom = poser.NewGeometry()

# Create the prop and make it a child of isoField
metaBall = scene.CreatePropFromGeom(meshGeom,metaBallName )
parmStickness = metaBall.CreateValueParameter('Stickness')
parmStickness.SetValue(2.0)
metaBall.SetParent(isoField, 0, 0)
parmStrength = metaBall.CreateValueParameter('Strength')
parmStrength.SetValue(5.0)

# Position the prop in the middle of the iso field
geom = isoField.Geometry()
vertex7 = geom.Vertex(7)
sX = vertex7.X()
sY = vertex7.Y()
sZ = vertex7.Z()
metaBall.SetParameter('xTran',sX/2)
metaBall.SetParameter('yTran',sY/2)
metaBall.SetParameter('zTran',sZ/2)

# All done; now we have a the cube for the iso field, and a smaller cube to center metaballs in that field
``````