Smooth union, diff. and inter. only handles two shapes

[TLC123]

Can this be solved by recursively smin a list of shapes? I'm thinking Smin_recursive equals smin the lists first element with the Smin_recursive of the lists tail.

[doug-moen]

If you try to combine 3 shapes using smooth union, the order in which you combine them matters, and the results are not symmetrical, which looks bad. That's the problem I want to solve. I want an n-ary smooth union operation which is commutative and associative: you get the same results regardless of how the shapes are ordered.

Here's what the current situation looks like (smooth union of 3 circles). Try the following program:

let
blend3 n (s1,s2,s3) = smooth n .union (s1, smooth n .union (s2, s3));
circles = [for (a in 0..<tau by tau/3) circle 1 >> move [1,0,0] >> rotate a];

in
union (
    blend3 1.27 circles,
    union circles >> colour black
)

The two circles on the left are blended together first, and then the result is blended with the right circle. Notice that the output is not radially symmetrical.

[bluecube]

@doug-moen Do you have an idea of what the commutative and associative version it should look like when applied in this example? I tried to think about it a bit and I don't :-)

[doug-moen]

The underlying problem is that the smooth_min function is not distributive. If we define smin3 like this:

let
blend3 k (s1,s2,s3) =
    make_shape {
        dist p : smin3(s1.dist p, s2.dist p, s3.dist p, k),
        bbox : [ min(s1.bbox[MIN], s2.bbox[MIN], s3.bbox[MIN]),
                 max(s1.bbox[MAX], s2.bbox[MAX], s3.bbox[MAX]) ],
        is_2d : s1.is_2d && s2.is_2d && s3.is_2d,
        is_3d : s1.is_3d && s3.is_3d && s3.is_3d,
    };
smin3(a,b,c,k) = smooth_min(a,smooth_min(b,c,k),k);
circles = [for (a in 0..<tau by tau/3) circle 1 >> move [1,0,0] >> rotate a];

in
union (
    blend3 1.27 circles,
    union circles >> colour black
)

then we get the same result as above. To solve the problem, we need to start by finding a clever implementation of smin3 that gives the same result, regardless of the order the a, b, c arguments.

[bluecube]

I have an idea about how to bend Inigo Quilez's exponential smin to support more than two arguments... Oh well, now I have to install curv to try it :-)

[doug-moen]

We'll start by visualizing the min function.

union(
    everything >> colour white,
    make_shape {
        dist(a,b,z,t) = min(a,b);
        is_2d = true;
    } >> colour black)
>> show_axes

The min function is <= 0 whenever either of its arguments is <= 0. If you interpret this as a distance function, it's an infinite shape that covers all of 2D space, except for the positive X,Y quadrant.

[doug-moen]

Now let's visualize the smooth_min function: smooth_min(x,y,r). This is just like the min function, except that a fillet of radius r has been added to smooth out the sharp corner where the X and Y axes meet.

let
r = 1;

in
union(
    everything >> colour white,
    make_shape {
        dist(x,y,z,t) = smooth_min(x,y,r);
        is_2d = true;
    } >> colour black)

[doug-moen]

In order to generalize this to 3 dimensions, first consider the 3D min function: min(x,y,z). This is a distance function that covers all of 3 dimensional space, except for the positive X,Y,Z octant.

Next, consider the 3D smooth min function: smin3(x,y,z,r). This is just like the 3D min function. except that in the positive X,Y,Z octant, we have applied fillets of radius r all along the positive X, Y and Z axes. At the origin, all 3 fillets should meet to form an eighth-of-a-sphere shape (the analog of the quarter circle shape seen in the 2D case).

[bluecube]

... i can't install it due to some compiler problems.

What I meant was to dismantle the exponential smin into its two parts: float res = exp( -k*a ) + exp( -k*b ); and return -log( res )/k; and make use of the commutative and associative addition in the first half like this: float res = sum(exp(-k*a) for a in inputs); (if you forgive the python syntax :-) )

If this didn't work, my second try would be float res = sum(exp(-k*a) for a in inputs) / len(inputs);.

[TLC123]

Case 1: one or none of a, b, c... n is smaller than r. - > classic union. Case 2: two or more of a, b, c... n is less than r. In the n-dimentional corner. Here there may be a need to special-case every combination of two corners, three corners up to n corners. But my guess is that simply truncate all values larger than r to r and translate to the origo .

From there it's easy to think that the distance would be something like r minus mag(a, b, c... n) but if i remember radi in hyper dimensions there are some caveats.

On Mon, 7 May 2018, 15:50 Doug Moen, notifications@github.com wrote:

In order to generalize this to 3 dimensions, first consider the 3D min function: min(x,y,z). This is a distance function that covers all of 3 dimensional space, except for the positive X,Y,Z octant.

Next, consider the 3D smooth min function: smin3(x,y,z,r). This is just like the 3D min function. except that in the positive X,Y,Z octant, we have applied fillets of radius r all along the positive X, Y and Z axes. At the origin, all 3 fillets should meet to form an eighth-of-a-sphere shape (the analog of the quarter circle shape seen in the 2D case).

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/doug-moen/curv/issues/15#issuecomment-387071021, or mute the thread https://github.com/notifications/unsubscribe-auth/AKcAaXVOBTen_oFIEvhlhsmrKuVXK2B_ks5twFE_gaJpZM4T0osw .

[doug-moen]

One more special property of smooth_min that I should mention. If you visualize smooth_min(x,y,r) as a distance field,

let
r = 1;

in
make_shape {
    dist(x,y,z,t) = smooth_min(x,y,r);
    is_2d = true;
}
>> show_dist

Then you see that it is not a Euclidean distance field. The fillet of radius r is repeated endlessly without changing its radius. This is important. If a Euclidean version is used, then when 2 shapes are blended, the result is not guaranteed to be be Lipschitz(1), and sphere-tracing does not render the blend correctly. That is why "something like r minus mag(a, b, c... n)" will not give the desired results.

[TLC123]

I had some too. warnings was treated as errors??? Maybe you got something else.

On Mon, 7 May 2018, 16:00 Kuba Marek, notifications@github.com wrote:

... i can't install it due to some compiler problems.

What I meant was to dismantle the exponential smin into its two parts: float res = exp( -ka ) + exp( -kb ); and return -log( res )/k; and make use of the commutative and associative addition in the first half like this: float res = sum(exp(-k*a) for a in inputs); (if you forgive the python syntax :-) )

If this didn't work, my second try would be float res = sum(exp(-k*a) for a in inputs) / len(inputs);.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/doug-moen/curv/issues/15#issuecomment-387073805, or mute the thread https://github.com/notifications/unsubscribe-auth/AKcAaZtFk8rwbes1iG3gE5eixHbJQBLjks5twFNfgaJpZM4T0osw .

[TLC123]

let
r = 4;
sminr( x,y,r) =        if (x< r  && y<r)          r-mag(x-r,y-r)
         else       min(x,y) ;
in
make_shape {
    dist(x,y,z,t) = sminr(x,y,r);
    is_2d = true;
}
>> show_dist

that gets me

[TLC123]

for the black circles example i get this: looks like nested sminr produce expected results but I'm just not sure how to make smin3(a,b,c,k) recurse over a longer list smin3(list,k)

let 
sminr( x,y,r) =        if ( (x< r  && y<r )   )       r-mag(x-r,y-r) 
    else  min(x,y) ;
blend3 k (s1,s2,s3) =
    make_shape {
        dist p : smin3(s1.dist p, s2.dist p, s3.dist p, k),
        bbox : [ min(s1.bbox[MIN], s2.bbox[MIN], s3.bbox[MIN]),
                 max(s1.bbox[MAX], s2.bbox[MAX], s3.bbox[MAX]) ],
        is_2d : s1.is_2d && s2.is_2d && s3.is_2d,
        is_3d : s1.is_3d && s3.is_3d && s3.is_3d,
    };
smin3(a,b,c,k) = sminr(a,sminr(b,c,k),k);
circles = [for (a in 0..<tau by tau/3) circle 1 >> move [1,0,0] >> rotate a];

in
union (
    blend3 1.27 circles,
    union circles >> colour black
)

[doug-moen]

@TLC123: There is a problem with your sminr function, when it is used for blending shapes. The problem is that the blended distance field overestimates the distance to the surface of the shape. This is most evident when unioning two shapes that meet at an angle > 90 degrees. And this, in turn, causes the 3D previewer (which uses sphere-tracing) to not render correctly. There are visual artifacts: the shape seems to deform as you rotate it.

The distance field debugger displays purple or red in areas where the distance field is bad. (If you substitute smooth_min for sminr in the definition of blend below, then the problem goes away.)

let
sminr(x,y,r) = if (x<r && y<r) r-mag(x-r,y-r) else min(x,y);
blend k (s1,s2) =
    make_shape {
        dist p : sminr(s1.dist p, s2.dist p, k),
        bbox : [ min(s1.bbox[MIN], s2.bbox[MIN]),
                 max(s1.bbox[MAX], s2.bbox[MAX]) ],
        is_2d : s1.is_2d && s2.is_2d,
        is_3d : s1.is_3d && s2.is_3d,
    };

in
blend .5 (
    rect(1,5),
    rect(1,5) >> rotate(10*deg),
)
>> show_dist

[doug-moen]

For reference, here is Inigo Quilez's polynomial smooth min, which I currently use as smooth_min:

smooth_min(a,b,k) =
    let h = clamp( 0.5+0.5*(b-a)/k, 0, 1 );
    in lerp( b, a, h ) - k*h*(1.0-h);

And here is Dave Smith's soft min, which is exactly the same function as smooth min, but implemented using a different algorithm, which I think might be faster.

// credit: Dave Smith @ Media Molecule
// http://media.lolrus.mediamolecule.com/AlexEvans_SIGGRAPH-2015.pdf
// Same distance field as IQ's smooth polynomial min.
let
soft_min(a, b, r) =
    let e = max(r - abs(a - b), 0);
    in min(a, b) - e*e*0.25/r;

I find soft_min easier to work with. I used it as the basis for implementing chamfer.

So, I was planning to spend time contemplating these two functions, and figure out how to generalize one of them to 3D. Once the 3D case is solved, then generalize to N-D.

[TLC123]

Thanks for shedding light on these challenges. Intuitively i first thought that knowing the gradient normal would be very helpful for finding the correct rounding, reducing the rounding if gradients where co-planar and so on.

On Mon, May 7, 2018 at 6:53 PM Doug Moen notifications@github.com wrote:

For reference, here is Inigo Quilez's polynomial smooth min, which I currently use as smooth_min:

smooth_min(a,b,k) = let h = clamp( 0.5+0.5(b-a)/k, 0, 1 ); in lerp( b, a, h ) - kh*(1.0-h);

And here is Dave Smith's soft min, which is exactly the same function as smooth min, but implemented using a different algorithm, which I think might be faster.

// credit: Dave Smith @ Media Molecule // http://media.lolrus.mediamolecule.com/AlexEvans_SIGGRAPH-2015.pdf // Same distance field as IQ's smooth polynomial min. let soft_min(a, b, r) = let e = max(r - abs(a - b), 0); in min(a, b) - ee0.25/r;

I find soft_min easier to work with. I used it as the basis for implementing chamfer.

So, I was planning to spend time contemplating these two functions, and figure out how to generalize one of them to 3D. Once the 3D case is solved, then generalize to N-D.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/doug-moen/curv/issues/15#issuecomment-387129752, or mute the thread https://github.com/notifications/unsubscribe-auth/AKcAaZtAywGzGMt7MVAm6zIMKjV7V_3Oks5twHwdgaJpZM4T0osw .

[TLC123]

"Can this be solved by recursively smin a list of shapes?" Yes! Smin_recursive equals ->sort the list of distances in ascending or descending order. smin the lists first element with the Smin_recursive of the lists tail.

Works fine in my unfinished openscad polygonizer. Images here:

Parts of the code i used:

function evalv(p)=
max( 
abs(p.z)-0.7, 
let(r=6,hr=r/2)
minRlist([
tube(p+[sin(0)*r,cos(0)*r,0],hr),
tube(p+[sin(120)*r,cos(120)*r,0],hr),
tube(p+[sin(240)*r,cos(240)*r,0],hr),
tube(p+[sin(60)*r,cos(60)*r,0],hr),
tube(p+[sin(180)*r,cos(180)*r,0],hr),
tube(p+[sin(300)*r,cos(300)*r,0],hr)
]
,1)
);

function minRlist(l,r)=minRlistworker(  (quicksort(l)),len(l)-1,r) ; 
function minRlistworker(l,c,r)=c>0? minR(minRlistworker(l,c-1,r), l[c], r) : l[0]  ;

function quicksort(kvs) = 
  len(kvs)>0
     ? let( 
         pivot   = kvs[floor(len(kvs)/2)], 
         lesser  = [ for (y = kvs) if (y  < pivot) y ], 
         equal   = [ for (y = kvs) if (y == pivot) y ], 
         greater = [ for (y = kvs) if (y  > pivot) y ] )
          concat( quicksort(lesser), equal, quicksort(greater))
      : [];

 function minR(d1, d2, r = 1) = r > 0 ?
 let (m = min(d1, d2))(abs(d1) < abs(r) && abs(d2) < abs(r)) || (d1 <
 r && d2 < r) ? min(m, r - norm([r - d1, r - d2])) : m : let (m =
 min(d1, d2), rr = -r)(d1 < rr && d2 < rr) ? min(m, norm([d1, d2]) -
 rr) : m;

[doug-moen]

Thanks for pushing this forward. I need to try this in Curv.

A question. You said "sort the list of distances in ascending or descending order. smin the lists first element with the Smin_recursive of the lists tail." Is the result different, depending on whether you sort in ascending or descending order?

The code you wrote won't be accepted by the Curv geometry compiler--for one thing, recursive functions are illegal on the GPU, and the compiler has other limitations as well. But, fun fact: if you export to STL and you don't use the JIT compiler, then the geometry compiler won't be invoked, so the OpenSCAD code with the recursive quicksort will probably work in Curv under that scenario. So I might try to get your code working first using STL export, using MeshLab to visualize the result.

If I succeed with STL export, then getting it working on the GPU is the next step, but will require more work.

[TLC123]

Got me stumped there for a while on the pure-functional non-recursive sorting algorithm. That's a hard nut to crack. But then i remember that the do operation makes it easier to iterated sort.

ascending or descending should intuitively make some difference when more than two distances are within the smooth radious.

For a fast and loose version maybe its enougth to find the three or four smallest distances and in line the nested smin's.

D1=min(list, largerthan=0)
D2=min(list, largerthan=D1) 
D3=min(list, largerthan=D2) 
D4=min(list, largerthan=D3) 

[Digital-Monk]

As a hopeless pragmatist, I personally wouldn't mind if the order of application changed the final result. I'm just trying to get fillets on my 3D printed parts for mechanical strength, and the precise details don't really matter to me. :D I'm not sure how that would be any different from now, since we have to pick the order to pair everything up anyway. But it would be nice if instead of:

smooth r .difference (
    smooth r .difference (
        smooth r .difference (
            body;
            cut1;
        );
        cut2;
    );
    cut3;
)

I could just say: smooth r .difference ( body; cut1; cut2; cut3; )

It would be up to me to decide what order to do the cuts in, or whether I should union all the cuts into a single cutout object, or whatever.

On the other hand, I suppose I should try my hand at implementing the "all at once" fillet procedure that I saw somewhere for OpenSCAD: (Hoping I'm remembering this correctly)

1. Offset object by r.
2. Subtract the result of 1 from a "sufficiently large" bounding cube.
3. Offset the result of 2 by r.
4. Subtract the result of 3 from a slightly smaller cube.

This easily takes hours in OpenSCAD, but it doesn't look like that much work for Curv...

[edit] This is, in fact, instantaneous. Unfortunately, I don't understand offset properly, because step 3 doesn't work as this would need. In OpenSCAD, steps 1 and 3 are Minkowski sums with sphere(r)... That is what I get from offset in step 1, but not in step 3. Assuming that I ended up with a mitred SDF after step 2?. FWIW, my test code:

let
    cross = union [
        box(10,1,1);
        box(1,10,1);
        box(1,1,10);
    ];

    d = 4;

    balloon = offset d cross;

    cutter = difference [
            box(20,20,20);
            balloon;
    ];

    filleter = offset d cutter;

    filleted = difference [
            box(11,11,11);
            filleter;
    ];

    crosssection = box(inf,20,inf) >> move(0,10,0);
in
    row [
        cross;
        balloon;
        difference [
            cutter;
            crosssection;
        ];
        difference [
            filleter;
            crosssection;
        ];
        filleted;
    ]

Trying to figure out when I lost the exactness of the distance field:

• box is exact
• cross should be exact because union of exacts should be exact
• maybe lost it at balloon (a positive offset of cross)? Is offset of exact also exact, or mitred? If this lost exactness, then nothing after that will be Minkowski-like and I'm sad...

[doug-moen]

@Digital-Monk said: smooth r .difference ( body; cut1; cut2; cut3; )

You can use

reduce [nothing, smooth r .difference] [body, cut1, cut2, cut3]

if you want to iteratively apply smooth r .difference to successive pairs of shapes.

reduce is part of the standard map-reduce-filter vocabulary of functional programming. nothing is a standard shape with nothing in it, the geometric equivalent of zero, and smooth r .difference is of course a function value that you are passing as an argument.

[Digital-Monk]

Ah, thank you! I knew there was a way to "accumulate" over lists, but I forgot the term was reduce. Definitely cleaned up my code...

[doug-moen]

@Digital-Monk said Unfortunately, I don't understand offset properly, because step 3 doesn't work as this would need.

Unfortunately, the implementation of offset is a big cheat. It should be equivalent to minkowski of a sphere, but that equivalence breaks down when you use it with union/intersection/difference, or when you use non-similarity transformations, or a few other cases. Currently, offset is defined as the isosurface of a signed distance field, which means it isn't equivalent to minkowski of a sphere when the signed distance field is non-Euclidean.

Fixing offset is on my TODO list, but doing it intelligently, with good performance, requires different algorithms for different cases.

My plan for implementing offset intelligently for unions, intersections and differences involves rewriting the offset operation in terms of smooth union, diff and inter. on multiple shapes. Which means that this plan is waiting on bug #15 to be fixed.

Which means that implementing smooth union et al in terms of offset is probably not the way forward.

[doug-moen]

Um, so the problem is that the set operations don't produce exact/Euclidean distance fields.

union is exact on the outside, but not on the inside. difference is mitred, not exact, on the outside.

And... I never documented this behaviour. I should do that. [Done now.]

[Digital-Monk]

So it was still OK at balloon (offset of union of exacts is still exact), but then I subtracted that out of a bounding cube with difference, resulting in a mitred SDF, so everything after that was similarly mitred. OK. Good to know. (Oh, thanks again for the reduce pointer -- it simplified my code for a car-mounted Qi wireless phone charger/holder quite a bit)

[Digital-Monk]

I'm a doofus. I obviously haven't grokked SDF modeling yet and keep thinking normal CSG...

In step 2, the only reason for using difference is to get the negative of balloon. In CSG, you would subtract balloon from a surrounding box, but with an SDF, all I need to do is invert the sign.

Would that work?

1. boxes are exact
2. union of exacts is exact on outside
3. offset of exact is [edit] mitred, so probably no improvement... [edit: moot after prior edit] 4. negative of exact is exact?

[edit: moot after prior edit] If so, then both of the subtractions/differences above simply become SDF negations...

So, that won't let me do the trick either, but just out of curiosity, how do you negate an SDF? I feel a bit stupid asking, but since it's the dist member of a shape record, I'm not sure how to negate it. Do I need to "copy construct" a new record, and negate dist as I do it?

A breadcrumb here for any other "machinists" out there. Hopefully it might make the reduce pattern a little bit easier to use/read. mask was the best I could think of for the intersection form...

let
    weld = fillet -> objects -> reduce [nothing,smooth fillet .union] objects;
    mill = fillet -> objects -> reduce [nothing,smooth fillet .difference] objects;
    mask = fillet -> objects -> reduce [nothing,smooth fillet .intersection] objects;
in
    // example usage
    weld 0.5 [
        box(10,1,1),
        box(1,10,1),
        box(1,1,10)
    ]

(Also, I apologize if I'm off-topic or otherwise abusing the issue system or this issue in particular. Just excited with a new tool.)