Overhauled API, rewrote docs and examples

[amatulic]

Consider this a draft for review. All examples tested and working with new API.

Slight speed improvement using matrix operations to calculate metaball contributions to the bounding volume.

I kept the axis parameter on the torus because I found it convenient to have the torus start aligned a particular way before applying rotation. axis parameter has been removed.

I also retained the ellipsoid because there's an ambiguity when scaling and rotating a sphere (you get different results depending on the order, and while that order may always be the same, the user may not know that). I found it clearer to have an ellipsoid with defined axis, which I then rotate.

[adrianVmariano]

Here's a list of various outstanding issues and topics (some maybe unresolved) from our discussion:

1. Move computation of suppress and 1/(r/coef) outside of all the mb functions into _metaball_fieldfunc so that the individual metaball functions only compute and return r/coef.
2. Consider stiffness/influence measure, maybe varying exponent on (coef/r)? How is this concept related to or independent of the existing cutoff idea?
3. Assess/revise suppression, perhaps to have negligible effect until near the suppression point and then smooth drop off. Logistic function?
4. Examples of the cutoff in action with simple case like 2 balls. What goes wrong. How do you fix it with cutoff. What goes wrong if you take the cutoff to an extreme
5. Correct roundcube so that it has r=(dv.x^e+dv.y^e+dv.z^e)^(1/e). Remove exponent from coef. Have special case check for the perfectly square case when squareness is 1.
6. After the above, roundcube should be able to replace cube
7. Have function interfaces parallel standard openscad, e.g. r/d for sphere, h/r for cylinder, size for cube
8. remove ellipsoid(?) Does it do something different than scaled sphere
9. Make isovalue 1 the default and make isovalue a named parameter (?) (Related to 7)
10. Add new cylinder function with rounded caps
11. Replace E field language with discussion of a more direct CAD idea that we are making objects that meld together when they get close.
12. Example of when/how you would adjust isovalue for a goal. What is problem arises with the default isovalue and how do you solve it, and what goes wrong if you take it too far
13. Would it make sense to change isovalue to threshold? (Maybe more understandable to users? I'll admit I find isovalue easier to type...)
14. Fix tetrahedron example to be math-free
15. Consider changing order so metaballs appears first because they are probably more accessible/useful to users
16. Merge isosurface with isosurface_array
17. Simple basic example showing 2 objects separated, showing them closer together starting the blend, and then showing them very close, fully combined into one object.
18. Pass through spin, anchor, orient, cp (I think) and atype to vnf_polyhedron. Also I think convexity.
19. Remove origin parameter from isosurface_array
20. Add 2d version?

That's all I could find going back and searching through the discussion.

[amatulic]

Providing updates to the points above as I progress. This will take a while.

1. Moving 1/r^n and suppress outside individual functions not done. Not only is _metaball_fieldfunc no longer in use, but I also disagree that one value should apply to all metaballs. I deliberately designed this to allow for flexibility of setting different values for different balls.
2. Done, added mb_shaping() function that is called at the end of each _mb* function. Revised documentation to account for this. There is a speed penalty however, about 1.5 sec for Example 1.
3. Logistic function doesn't hit zero. This one is nice $y=\frac{1}{2}\cos\left(180\left(\frac{x}{c}\right)^{n}\right)+\frac{1}{2}$ - https://www.desmos.com/calculator/b8x99jlpmp - with $n=4$ the falloff starts halfway to the cutoff $c$. Small but noticeable impact on execution time. Already impelemented in # 2.
4. Done in first 5 examples.
5. Done, mb_roundcube corrected.
6. mb_cube removed, now handled by roundcube with squareness=1. Note that with surfaces generated by marching cubes, one can never get sharp edges; by design, a surface always cuts across corners of a voxel cube. This is why the octahedron is so sharp.
7. Done, metaball function arguments now conform to their OpenSCAD or BOSL2 counterparts.
8. Done, ellipsoid removed.
9. Done, although I disagree that isovalue should be 1 or any other default, given examples I have found online, but I suppose it makes sense for the way we're designing the field functions.
10. Cylinder function added.
11. Done, removed most references to e-field, rewrote things, needs review and copyediting.
12. I think I covered this with the first 5 new examples.
13. I explained in the documentation that for consistency we use "isovalue" throughout but it means "threshold" in the context of metaballs.
14. Tetrahedron example is now math free, but doesn't use loops because it got messy.
15. Documentation revised so metaballs appear first.
16. isosurface_array() is no more, it is now merged into isosurface().
17. How about the animation Richard made?
18. Done, all anchoring args and convexity now passed through to vnf_polyhedron.
19. Origin parameter removed, at the cost of complicating metaballs() because it has an array and a bounding box, but you can't pass both to isosurface(), and isosurface() now returns whatever internal bounding box it calculated centered on the origin. It's much cleaner to have an origin argument, or have bounding_box do double duty for that.
20. Considering 2D version for later, out of scope for this PR.

[adrianVmariano]

Note that the order of items above is arbitrary and not necessarily the natural order for implementation.

I think if we are going to add exponential influence adjustment then isovalue default DEFINITELY should be 1, because that's the only way that the shapes are invariant and you change the influence value.

Having a default doesn't mean you can't change it! It has become apparent that the isovalue is not a very important parameter, so the current ordering doesn't necessarily make sense. In other words, changing the current ordering may be the right thing to do.

[amatulic]

By the way, here's the integral for the cone-shaped metaball, with the "wire" of length L going from $z=-L/2$ to $z=L/2$. If $r_1$ and $r_2$ are the radii of the cone ends, the radius interpolation along the length of the wire would be

$r(z) = \frac 1 2 (r_1+r_2) + \frac z L (r_2-r_1)$

The field contribution at point $(x_0,z_0)$, with an influence of 1/p, is

$\int_{-L/2}^{+L/2} r(z)\left[x_{0}^{2}+\left(z_{0}-z\right)^{2}\right]^{-p/2} dz$

Wolfram Alpha doesn't want to solve it.

[adrianVmariano]

I wonder how dumbed down that is compared to the commercial Mathematica. In any case, the solution if it exists probably depends on special functions that are hard to compute, or is a very long complicated expression, or both.

[amatulic]

It would be nice if I could do anything with it. If there's a closed form for specific values of p, like p=1 or p=2, that would be helpful.
I tried fooling Wolfram Alpha by substituting constants, which sometimes works:
$\int_{-a}^a \frac{b + c x}{p + (q-x)^2} dx$
assuming the exponent is -1... looks simple but it still refused to solve it.

[adrianVmariano]

Solution is not horrible:

[adrianVmariano]

Regarding 19, why can't you just apply move(origin,vnf) to the result of isosurface? Part of the argument for eliminating origin is that it's very easy to do the same thing with move, so there shouldn't be a significant complication.

[amatulic]

The solutions aren't horrible individually, but there doesn't seem to be a straightforward way to have a single influence exponent like the other metaball types.

Regarding 19, why can't you just apply move(origin,vnf) to the result of isosurface?

Because that doesn't work.

Say my isosurface bounding box is [[20,20,20], [40,40,40]], and the voxel size is 3. Then internally, because the box size isn't divisible by the voxel size, the actual bounds are [[20,20,20], [41,41,41]], which fits a 7x7x7 array of voxels. Then it returns the isosurface with those new bounds centered on the origin, so the bounds are now [[-10.5,-10.5,-10.5], [10.5,10.5,10.5]]. But my original origin was supposed to be [20,20,20]. If I then perform move(origin, vnf), the new bottom of the bounding box is at [9.5,9.5,9.5], nowhere near where I wanted it to be.

That is why the origin parameter was useful.

However, I now have metaballs() calculate the bounding box just like isosurface() does so that it knows what's being returned and adjusts for it. It's a kludge, inelegant to do the same thing twice, I don't like it, but it meets the requirement of removing origin.

[adrianVmariano]

If adding a function call is adding 1/2 second I'll bet the run time for those integral solutions would be substantially longer. And I also suspect it'll create an object that behaves differently than the existing balls types. Regarding influence you could assume that the output is 1/r and apply influence with an exponential.

Regarding the bounds, I hadn't realized this issue at all. My first thinking was that voxels should be adjusted to fit the bounding box, but that creates non-square voxels, which I'm not sure is a good idea. I don't like the idea of the output having an unknown size, but I'm also not sure I see a solution I like. Like specify the bounding box in voxels instead of units?

It does seem like if you're going to alter the size of the bounding box, the new box should be centered on the requested box. That is, if you request [20,40], you should not be computing [20,41] but rather [19.5,40.5], so that the error is symmetric.

You could calculate the number of voxels on a range, figure out the adjusted start, and then use count() to generate the data points.

But rounding error aside for voxel fit, I don't understand how you can end up "nowhere" near the right place. In your example the bounding box requested is [20,40], so you calculate the origin to be 30. You compute and get a solution that you think is [-10,10] but it's actually [-10.5,10.5] due to rounding errors. You correct the origin by adding 30 to this and now you think you have [20,40] but you actually have [19.5,40.5], which the answer incorrectly shifted by 1/2 unit. But if you centered the bounding box it would be correct. I have no clue where you pulled the 9.5 bound from in your example above. And it seems like the origin parameter is mainly needed to compensate for not centering the bounding box on the requested bounding box.

Note that it's also ok to have internal parameters between the module and function or other internal functions. But if I have overlooked something and there's a compelling reason why users need the origin parameter, we should understand that.

[adrianVmariano]

Would it make sense to allow boundingbox = [xsize, ysize, zsize] to produced a centered box? And maybe even bounding_box=5 to make a centered cubic box.

[amatulic]

The simplest and most elegant solution is to pass _origin as a hidden parameter and use it if it's defined. isosurface() needs to know an origin regarless, so it uses either the one for the bounding box supplied, or calculates it from the center of the box implied by the array. So it already has a notion of "origin". If a private _origin argument is defined, it would just use that.

Each voxel that contains the surface is a data structure that includes an x,y,z coordinate among other things, and that coordinate is determined by the origin of the bounding box.

In the case of metaballs, a bounding box must be specified, but that can't be passed into isosurface() when passing an array, so it would be cleaner to pass the origin of the bounding box as a private argument. isosurface() needs an origin anyway to work, whether it calculates one or receives one.

I've just made this change, which eliminated a handful of lines of duplicate code, and also eliminated the need to dereference vnf to move the coordinates in it.

[adrianVmariano]

I still don't understand why you need the origin parameter. Did you fix the off-centering of the bounding box?

It seems like if isosurface produces a centered output then you will always know what to pass to move() to fix the output if necessary. If you don't for some reason, then users might have this problem too, so we need to understand that.

[amatulic]

Earlier timing tests were using the wrong version of code without the matrix operation to build the fieldarray. Now using the correct version.

Timing tests (approx average of 5 runs each case), using Example 1.

Case 1. This is my baseline, and takes 4.4 seconds to execute.

function mb_shaping(r, coeff, cutoff, influence) =
    sign(influence) * 0.5*(cos(180*(r/cutoff)^4)+1) * (coeff/r)^(1/abs(influence));

function _mb_sphere(dv, coeff, cutoff, influence) =
    let(r=norm(dv)) r>=cutoff ? 0
    : influence != 1 ? mb_shaping(r, coeff, cutoff, influence) : coeff/r;

Case 2. Moving the r>=cutoff and influence==1 checks into mb_shaping(), execution time increased to 5.0 seconds.

Case 3. Removing the check for influence==1 and just letting it calculate, execution time increased to 5.6 seconds.

Case 4. Moving the r>=cutoff check back into mb_sphere and letting mb_shaping() calculate: 5.5 seconds (one trial was 5.6 but average was 5.5).

Case 5. Bypassing mb_shaping() altogether and putting everything inline in mb_sphere: slightly under 4.4 seconds, negligible difference from "baseline".

This is all with influence=1. If I set influence=1.1 for all 5 balls in Example 1, the execution time for the baseline case becomes about 5.7 seconds, comparable to case 3 and 4 (because the cos x^4 function is being calculated every time).

[adrianVmariano]

So if I am understanding this correctly, simply adding a function call is increasing time from 4.4s to 5s, which is a 14% increase. That suggests that it might be necessary to inline everything. There is one possibility worth exploring if you want to maximize run time, which is to change the function literal. You aren't showing everything in the above quote because you don't give the definition of mb_sphere. It could look something like

function mb_sphere(r,influence, cutoff) = 
    influence==1 && cutoff==INF ? function (v) _mb_sphere_basic(v,r)
    cutoff==INF ? function (v) _mb_sphere_influenced(v,r,sign(influence),1/abs(influence))
                           : function(v) _mb_sphere_full(v,r,sign(influence),1/abs(influence),cutoff);

Something along those lines. So then there is a version of mb_sphere that just computes r/norm(v), another version that computes sign_inf*(r/norm(v))^inv_influence, and a third version that does the works. This is obviously more of a nuisance to write and maintain., But presumably it will deliver the best overall response since there is no extra function call and in baseline cases, no unnecessary tests or computation at all.

[adrianVmariano]

Note also, which parameter do you think is more important to users, the sphere influence or the cutoff? I imagined influence to be more important....but I don't know....which ever one is more important should go first so it can be used positionally without having to set the other one.

(Oh, also the tests I show above ==INF I think don't work. Need to look up how to check that. Or maybe just have no default and check if the parameter is set.)

[adrianVmariano]

Test for infinity is is_finite()

[amatulic]

Inlining all of that would actually have four separate outcomes, for each combination of influence =1 or !=1 and cutoff>r or <r. I think what I have now is a good compromise between maintainability and speed.

which parameter do you think is more important to users, the sphere influence or the cutoff?

As I play around with it, I find I use cutoff more often. Influence is nice, however, to cause all metaballs except one to grow or shrink by changing the influence of that one metaball.

It looks like the remaining open tasks in the list above are mostly examples, and having arguments like r and d like in OpenSCAD. You know I am not in favor of that; I can imagine the complaints, from those who don't read documentation carefully, about not getting the diameter asked for. I prefer coeff or some name that clearly means the coefficient is analogous to a radius or size, but not actually radius or size. Has anybody else weighed in on this?

Keep in mind also that every metaball function would need logic checks and conversions between d and r. I think this is unnecessary overhead in time-critical functions for no actual benefit.

[adrianVmariano]

By writing the stuff in here we kind of hide it from everybody else.

There will not be complaints from people about balls not having the "right" size. The other thing is that we should write this so that it's the easiest to use and easiest to remember how to use. I would say it's more likely we'd get complaints from users that mb_sphere(d=2) doesn't work. The problem with "coeff" is that it means nothing and doesn't convey what the argument does. So actually it's more hostile to users who don't read the documentation. And a user reading the cheatsheet, or just looking at the usage message won't have guidance about what the argument means. It will only get worse if you have multiple arguments to a ball. It's also more limiting because you can't for example have the option of giving diameter or radius to a sphere, or specifying cubes by size in analogy with how cube() works.

In terms of what other people think, we know that lightwave and blender both have you specify the radius of your ball. So it appears standard to do that.

[adrianVmariano]

I was looking at the docs for blender metaballs and realized that they split negative influence off as a totally separate parameter. We might consider doing this. As noted, right now the sign is stripped from influence anyway, so it is a separated calculation. Setting a ball to negative influence appears to be rare and also a very confusing mistake to make. If we required influence to be positive we could specify negative balls with negative=true which would be very visible and consistent with the already-existing separation of the sign and numerical influence.

[adrianVmariano]

Does there need to be mb_custom(f,cutoff,influence,negative) that takes a user-supplied function that returns dist?

[amatulic]

Does there need to be mb_custom(f,cutoff,influence,negative) that takes a user-supplied function that returns dist?

I don't see the need. Besides the fact that the user function isn't invoked the same way, the only requirements for the user function are that it takes the distance vector dv argument and returns a single value. Anything else is up to the user, including whether to have influence, cutoff, or negative as arguments.

[adrianVmariano]

Ok. I was just thinking that providing mb_custom would automatically provide the proper scaling, influence and cutoff. I don't really think custom functions are likely to get much use. Consider how much trouble we're having coming up with a cone function.

[adrianVmariano]

You might consider using list_shape in isosurface to get nx, ny and nz. That will check that the matrix given is actually a cube. (It will return undef for variable dimensions.

The code that does the bounding box is kinda hard to follow. Just reading the code, is the bounding box centered on the requested bounding box now? I think we need to understand why you think you need the origin parameter, because if there's a real need there that I don't understand, then users might actually need it. And if there isn't a real need, your code could be cleaner because you shouldn't need it either.

I would have done the bounding box calculation without all those half-voxel additions and roundings, something like

origin = mean(boundingbox),   // center of bounding box
box_counts = v_ceil( (boundingbox[1]-boundbox[0]) / voxel),
box_width = box_counts*voxel,
point_sets = [for(i=[0:2]) lerpn(origin-box_width[i]/2,origin+box_width[i]/2, box_counts[i]+1)]

which I guess is longer code? The point sets could also be separated into xset, yset and zset instead and instead of calculating x, y and z in the loops you just pull them from the list, which will be (slightly) faster. (In the case of calculating the function you can loop directly on x, y and z, which is a little bit faster as well.)

I can't quite tell what's going on with _origin. If the point is just to handle positioning of the bounding box for arrays I don't think it's necessary assuming that you center the box. After calculation you just apply move(origin) with the above definition of origin. You can compute only origin in isosurface() and the rest of the bounding box stuff can happen once in the private function instead of being repeated.

Am I missing something?

[adrianVmariano]

I think the API is looking good and just a few details remain to be cleaned up:

1. Resolve issue with need for origin argument and centering of the box
2. Make sure to use is_finite() to test parameters where infinity is not a meaningful input (every parameter except cutoff)
3. For parallelism change mb_roundcube to mb_cube. Might consider bumping default squareness to 0.75. And also accept an array size argument. I know it seems crazy since I asked for ellipsoid to be removed because scale could do it all...and it can...but cube() takes an array size argument. So for parallelism it should be allowed to have a size as a scalar like currently done, or a 3-vector. And if you're really missing ellipsoid you can use cube with squareness=0! I'm not certain but possibly the fastest way to apply this 3-vector would be to pass the inverse size as input either 1/size or [1/size.x, 1/size.y, 1/size.z] and start the code with let(dv=dv*size) which will do the right thing for both the scalar and vector cases. Remove the reference to size that is there already, and that should be it. Possibly using v_abs is slower than necessary because it does error checking that its input is a numeric vector. Not sure how big a difference that makes.
4. For parallelism with cylinder, the arg order needs to be length,r first of all, for consistency with cyl and cylinder. Also d should be supported, and for height, also l, h and height. You can handle the latter with length = one_defined([h,l,length,height],"h,l,length,height") before you do any error checking on h. (We got into this annoying situation with the arguments because of varying use across the library and also inconsistent use by OpenSCAD itself and it seemed like supporting all four of these was the only robust consistent solution, so we always support all four of these arguments together. Note that you used "length" but native OpenSCAD calls this argument "h", presumably for height.)
5. For parallelism with torus at least change rmaj and rmin to r_maj and r_min. There's a bunch of other options on torus...kind of a lot. Ideally you would handle them so mb_torus handles all the options that torus does. But it's a lot, so if you want to skip them I won't complain. Actually you could probably just cut the arg processing code out of torus() and drop it into your code, which wouldn't be so horrible to do. It looks like when the dust settles it has defined the major and minor radii from its 8 input parameters.
6. Is it necessary to call vnf_merge_points? How bad is the time penalty?

[amatulic]

1. Resolved in chat, also made metaball bounding box centered around the requested bounding box in a less clunky way, also did same for bounding box passed with a function into isosurface().
2. Done, is_finite checks added.
3. mb_roundcube is now mb_cuboid, mb_cylinder is mb_cyl. Left squareness as 0.5 because this is the default in the only other place squareness is used, in squircle() (and in fact mb_cuboid calls a private squircle function to get the exponent).
4. Done, supporting h|l|height|length in mb_cyl
5. I don't understand. mb_torus already uses r_maj/d_maj and r_min/d_min as arguments. What "bunch of options"? (Edit: Done, now supporting torus() args as described in comment below)
6. We should absolutely not call vnf_merge_points. The time penalty is immense. I just tried this on the final example (Neovus) in isosurface(). Without vnf_merge_points it finished in 11 seconds. With vnf_merge_points it took nearly 6 minutes. It may be more feasible when we get to 2D contours, in fact I think we'll need it to generate a path from a soup of two-point segments.

Now that's all done, I still need to come up with missing examples from the previous to-do list.

[adrianVmariano]

5. The version of mb_torus I looked at used rmaj and rmin instead of r_maj and r_min. The regular torus() command that you are parallel to has:

//   r_maj = major radius of torus ring. (use with 'r_min', or 'd_min')
//   r_min = minor radius of torus ring. (use with 'r_maj', or 'd_maj')
//   ---
//   d_maj  = major diameter of torus ring. (use with 'r_min', or 'd_min')
//   d_min = minor diameter of torus ring. (use with 'r_maj', or 'd_maj')
//   or = outer radius of the torus. (use with 'ir', or 'id')
//   ir = inside radius of the torus. (use with 'or', or 'od')
//   od = outer diameter of the torus. (use with 'ir' or 'id')
//   id = inside diameter of the torus. (use with 'or' or 'od')

so 8 different args in total that can be combined in different ways.