SDF optimizations

[pca006132]

It seems that OpenSCAD folks are interested in adding SDF support (via libfive), so I was curious about our SDF performance comparing with libfive and curv (libfive + GPU). It turns out the speed of libfive is roughly similar to us, and curv is a lot faster. So I looked into the code a bit and was thinking if we can do something about the performance.

1. Can we do something like search refinement? E.g. When |sdf(p)| >= r, we skip evaluating points in a radius of r. Probably not radius but something like an octree. This can probably reduce most of our queries.
2. Do we really need multithreaded hashtable in our case? I did a simple profiling and found that the atomic operations for hashtable insertion are taking the majority of time (like >70% of total SDF time, including query time). I am thinking about using thread local hashtables, and each thread will handle a block of morton codes. The block assignment is recorded into a global hashtable (smaller and less access), which can be accessed later. This should remove those atomic operations.
3. Provide expression tree API for python. Currently our Python API is very flexible but it is not parallelizable due to GIL. We can probably provide a similar API as libfive that allows users to build a sdf, and we can interpret the expression in multiple threads.
4. JIT. Given an expression tree, it should be pretty trivial to use other libraries like MIR to do JIT.
5. Compute shader. Similar to JIT, it should be possible to generate compute shader, and considering our library is originally designed to be CUDA accelerated, this should hopefully not be too hard? From the performance of curv (10x faster than us), GPU acceleration for SDF should be worthwhile.

[pca006132]

@gsohler in case you have any opinion about this. I don't have much idea about sdf.

[elalish]

So, 1 assumes you have a proper distance function (rather than just something that changes sign), which is assuming a lot for an SDF. It would disallow too many nice functions, including my heart example, I'm pretty sure.

3 sounds cool, but how complex would it be? Also, if curv already does such a good job, why not just use it? Does it not guarantee manifoldness (a common marching cubes problem)?

Regarding 2, 4, and 5, I trust your judgement when it comes to speed. Certainly SDF is a highly parallelizable problem. I'm also not sure which approach is faster: evaluating the SDF values twice (once to set up the data structures, and again to fill in their values), or just once like I do now where you have to find a way to index into a disorganized array (hence the hash table).

[pca006132]

There is a discussion about mesh quality here: libfive/libfive#284

[elalish]

Honestly, until we have a bunch of users of SDF complaining about our speed, I'd prefer to have a small, simple code base. Speed is great, but robustness is much more important, and big part of that is code clarity and simplicity. I think your other parallelization ideas sound simpler than trying to figure out a distance bound.

[pca006132]

Indeed, I think the simple stuff can get us pretty far away. The major problem is how to get more users, OpenSCAD SDF integration is still in its early stages, and it feels like most people are not using curv/libfive, maybe sdf is a bit too complicated for normal users.
I wanted to push for this for OpenSCAD as an alternative to libfive, by showing that our performance is great as well and have better guarantee.

[hrgdavor]

is curv still actively developed? I was following it on discord, and it went silent some time ago. And it was far from complete at that time.

[pca006132]

I wish to know as well, I feel like I am navigating abandoned cities when I check various sdf projects: Nice claims, seems working well and quite a lot of stars, but somehow there is not much recent activities and no one is talking about them.

[gsohler]

Hi Emmett, i pass whatever distance function to libfive the user provides. libfive can handle both proper distance functions and just values with a sign, but proper distance values work better. Yes, i know curv and i used it a lot in the past and its very good at evaluating SDF's However, My vision is to have a tool which can handle SDF and BREP objects in the same language/program. This is why i have taken the effort to use libfive in openscad. Now having the power of python i can even generate SDF equations or CSG trees programmatically. I am sure there will arise applications for that. Some time back I learned about the marching cubes algorithm and its a perfect fit for evaluating SDF's. its bullet proof and very usable for parallelism. However, on the other hand it runs in O(n^3) , it sounds like shooting birds with cannons which is not good enough for me. My Vision is to establish an Algorithm to evaluate SDF's in O(n^2) only. My vision is that I enclose the object with a big cube first. Then I iteratively add points to the mesh and attach them to the object surface similar like vacuum deep drawing process of plastic parts(mainly packaging). I am aware that there are many problems with will arise(like separate objects or holes inside objects) i hope that i find the right algorithm for that soon)

…

On Thu, Sep 21, 2023 at 5:39 PM Emmett Lalish ***@***.***> wrote: So, 1 assumes you have a proper distance function (rather than just something that changes sign), which is assuming a lot for an SDF. It would disallow too many nice functions, including my heart example, I'm pretty sure. 3 sounds cool, but how complex would it be? Also, if curv already does such a good job, why not just use it? Does it not guarantee manifoldness (a common marching cubes problem)? Regarding 2, 4, and 5, I trust your judgement when it comes to speed. Certainly SDF is a highly parallelizable problem. I'm also not sure which approach is faster: evaluating the SDF values twice (once to set up the data structures, and again to fill in their values), or just once like I do now where you have to find a way to index into a disorganized array (hence the hash table). — Reply to this email directly, view it on GitHub <#561 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCO4MTF2KW4GEEBU7QPHU3X3RNTTANCNFSM6AAAAAA5BXFX7I> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[elalish]

Agreed. Also, the output in number of triangles is worst-case O(d^3), and commonly O(d^2). I like to think of the number of triangles as n, so really we're talking about O(n^(3/2)) vs O(n) for some common cases. That makes me feel like it's not so ridiculous, especially if it's well-parallelized. And I think the reliability is very much worthwhile, as the vacuum drawing approach sounds very difficult to guarantee.

[elalish]

Also, I'd say for this, memory scaling is a bigger problem than compute scaling. Our current algorithm is O(d^3) in compute, but only O(n) in memory.

[ofmetal]

This is why i have taken the effort to use libfive in openscad.

Do you mean pythonscad.org ?
I tried that one couple of month ago, it was pretty buggy.

[ofmetal]

Our current algorithm is O(d^3) in compute, but only O(n) in memory.

Digging into the code, seems it is actually O(d^3) in memory.
Take a look, sdf.cpp:

Manifold Manifold::LevelSet(std::function<double(vec3)> sdf, Box bounds,
                            double edgeLength, double level, double tolerance,
                            bool canParallel) {
...
  const vec3 dim = bounds.Size();
  const ivec3 gridSize(dim / edgeLength + 1.0);
...

  const ivec3 gridPow(la::log2(gridSize + 2) + 1);
  const Uint64 maxIndex = EncodeIndex(ivec4(gridSize + 2, 1), gridPow);
...
  Vec<double> voxels(maxIndex);

So, If I set the bounds to (-12, -12, -12) .. (+12,+12,+12).
And a edgeLength to 0.1.
This will set grid size to 241x241x241 // 241 = (12+12)/0.1 + 1
Log2 of which is (8, 8, 8) - the gridPow.
That constructs maxIndex equal to (241+2) * (2^16 + 2^8 + 2^0) * 2 + 1 = 31975399
Which is a 256 (power-of-2 upper bound of 243) squared and then multiplied by 243.
So it's lower bound is ~2*d^3.

We see there is
Vec of maxIndex size.

It's constructor allocates sizeof(double) * maxIndex = 8*31975399 bytes = 255 803 192 bytes ~= 2.28 * 241^3 doubles.

Vec(size_t size) {
    reserve(size);
    this->size_ = size;
  }
void reserve(size_t n, bool seq = false) {
    if (n > capacity_) {
      T *newBuffer = reinterpret_cast<T *>(malloc(n * sizeof(T)));
  

[elalish]

Yeah, I actually recently changed this - it is now On^3 in memory, but with 7x few function evaluations, which picked up the speed considerably. The next optimization that I haven't quite gotten to is batching: it shouldn't be too hard to make this block of memory a fixed size and then iterate it over the sample space to keep the memory from exploding. If you'd like to take that on, I'd very much appreciate it!

[gsohler]

Hi Emmett, I am excited to read, that processing time goes commonly to O(d^2) and i dont understand why this is achieved. in marching cube algorithm you really have to march all cubes to find out where is actually void. this results in O(d^3) where is the hidden shortcut ?

…

On Thu, Sep 21, 2023 at 8:26 PM Emmett Lalish ***@***.***> wrote: Agreed. Also, the output in number of triangles is worst-case O(d^3), and commonly O(d^2). I like to think of the number of triangles as n, so really we're talking about O(n^(3/2)) vs O(n) for some common cases. That makes me feel like it's not *so* ridiculous, especially if it's well-parallelized. And I think the reliability is very much worthwhile, as the vacuum drawing approach sounds very difficult to guarantee. — Reply to this email directly, view it on GitHub <#561 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCO4MSIQV6N7CJHOOILDK3X3SBFZANCNFSM6AAAAAA5BXFX7I> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[elalish]

I think you misread what I said - can you quote the part you're referring to?

[gsohler]

Hi Emmett, yes you are right you stated: "Agreed. Also, the output in number of triangles is worst-case O(d^3), and commonly O(d^2). " correct, but I fear you need to find out in O(d^3) time that the amount of triangles is O(d^2) if you know of any O(d^2) algorithm, even if they probably cannot handle things like gyroids, I'd like to learn about them ... Guenther

…

On Thu, Sep 21, 2023 at 9:09 PM Emmett Lalish ***@***.***> wrote: I think you misread what I said - can you quote the part you're referring to? — Reply to this email directly, view it on GitHub <#561 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCO4MXD2WWISE62GMS5W6LX3SGFRANCNFSM6AAAAAA5BXFX7I> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[pca006132]

I managed to reduce the number of atomic operations required, thus minimizing contention and improved the performance of sdf quite a bit. (about 30% improvement in some cases? the code is quite ugly so I am not posting now)
Can probably improve it more but the interface for HashTable will need some overhaul, so will take some time.

However, I found that importing meshes into manifold is pretty slow. For example, when I use size = 60 for GyroidModule (default is 20), ~30% of the time is spent on evaluating sdf, and the remaining 70% is spent on converting Mesh into Manifold. In particular, SplitPinchedVerts and CreateFaces (finding the connected component + other per face operations) took the majority of time in single threaded mode and cannot be parallelized. This is probably a general enough problem that we should open a separate issue for, because it will also impact other mesh import.

[pca006132]

In particular I should put this next to my monitor so I can see it whenever I think atomic is free:

[elalish]

Nice! Yes, that doesn't surprise me too much - I didn't put that much effort into optimizing our constructor. I figured generally there would be more time spent on further operations than on the input, but if it's exceeding the SDF time, that's not great. I bet SplitPinchedVerts could be improved a lot - it's in general a no-op, I believe provably so for an SDF. Finding problem verts should be parallelizable at least.

Connected components is definitely slow - I looked for quite a while for a decent parallel solution, but couldn't find anything for general graphs. The trouble is they mostly scale with the diameter of the graph, which for us is ~O(n^1/2), which is pretty bad. I wonder if we could create faces by hashing their normal vector and origin offset into precision-sized buckets? It would relate disconnected parts of the same plane, but I don't think that's a problem (that's what Boolean ops will do anyway).

[pca006132]

Yeah I don't think we can speed up connected components a lot. We can perhaps open an issue to track performance work on constructor later.
At least it would be nice if we can improve the performance for SplitPinchedVerts (or skip it all together with sdf).

[pca006132]

I just had an idea: Why bother with Lipschitz constants when what we are caring is just whether or not the function can cross zero at some neighboring points? We can just do a interval value analysis. Say the original SDF is $f: \mathbb{R}^3 \to \mathbb{R}$, we can reinterpret in another domain to make it $f': \mathbb{R}^3 \times \mathbb{R}^3 \to \mathbb{R}\times\mathbb{R}$, in the sense that the input is the bounding box of the input space, and the result is the min and max (approximation) of $f$ inside this input. Every operator is interpreted as interval arithmetic. If the resulting interval does not contain 0, we will just skip the whole bounding box. Supporting loops may be tricky because it involves fixed-point calculation, we can just don't provide loop construct :)

[gsohler]

Not sure if I get your idea. Of course we can shortcut marching Cube algorithm by using intervals.but if sdf function is below zero in all 8 corners of the interval cube it does not mean that there is nothing inside.we still might have missed a rising and a falling edge inside.not sure if we can use intervals without being sure not to miss anything. If your idea is different, please explain your idea for dummies 😅 I don't understand these mathbb expressions unfortunately.

…

On Tuesday, October 31, 2023, Chun Kit LAM ***@***.***> wrote: I just had an idea: Why bother with Lipschitz constants when what we are caring is just whether or not the function can cross zero at some neighboring points? We can just do a interval value analysis. Say the original SDF is $f: \mathbb{R}^3 \to \mathbb{R}$, we can reinterpret in another domain to make it $f': \mathbb{R}^3 \times \mathbb{R}^3 \to \mathbb{R}\times\mathbb{R}$, in the sense that the input is the bounding box of the input space, and the result is the min and max (approximation) of $f$ inside this input. Every operator is interpreted as interval arithmetic. If the resulting interval does not contain 0, we will just skip the whole bounding box. Supporting loops may be tricky because it involves fixed-point calculation, we can just don't provide loop construct :) — Reply to this email directly, view it on GitHub <#561 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCO4MXYMNAVK4Q2VVRYLJLYCEAABAVCNFSM6AAAAAA5BXFX7KVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TIMZVGY3TO> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- Null

[pca006132]

The interval overaproximates the set of values the function can take. If it crosses zero somewhere within a cube, the interval should include 0.

[elalish]

If you know a way to quickly bound an arbitrary function, then sure, but that sounds hard? Also, what do you mean by loops?

[pca006132]

We don't need to bound arbitrary functions, we just need to bound those that are commonly used, like addition, multiplication, trigonometric functions. This is a common technique called abstract interpretation, I can write a demo later.

By loops I mean if the sdf contains some for/while loops. These loops make abstract interpretation harder and slower, and it seems that they are rarely used in sdfs.

[gsohler]

Just to understand your idea correctly:for your idea to work the equation must be clearly visible. A hidden equation just providing values for given inputs will not work, right?

…

On Tue, Oct 31, 2023, 16:59 Chun Kit LAM ***@***.***> wrote: The interval overaproximates the set of values the function can take. If it crosses zero somewhere within a cube, the interval should include 0. — Reply to this email directly, view it on GitHub <#561 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCO4MRLVW3K7XDNUSGQXS3YCEN4ZAVCNFSM6AAAAAA5BXFX7KVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TIMZXGAZTS> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[pca006132]

Yes. Let me give a 1D example:

f(x) = sin(x) * x

When this function is evaluated against [pi/2, pi], the range of sin(x) is [0, 1], and the range of x is [pi/2, pi], so the multiplication result range is [0, pi].

[gsohler]

Yes I understand now. For 1d you basically do zero value analysis and it can easily be done symbolically with the equation. However if you do with 2d the zero values on x depend on the actual value of y and it starts to be complicated very soon.can you help me to understand for 2d ?

…

On Tuesday, October 31, 2023, Chun Kit LAM ***@***.***> wrote: Yes. Let me give a 1D example: f(x) = sin(x) * x When this function is evaluated against [pi/2, pi], the range of sin(x) is [0, 1], and the range of x is [pi/2, pi], so the multiplication result range is [0, pi]. — Reply to this email directly, view it on GitHub <#561 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCO4MTDZ75BZDBTY5VZ3YLYCEQ7HAVCNFSM6AAAAAA5BXFX7KVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TIMZXGMYDC> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- Null

[pca006132]

No, we only work with the intervals without considering their relations, so 2D is similar. Just change the last x in the above example to y, the analysis just works.

[gsohler]

Yes understood. Using Intervals you can reduce the sdf analysis outer bounding box to the touching bound box of meshed object. This might be an advantage of the size is the meshed object is not known bit usually you set the size of the outer bounding box quite precise. But the equation sub-term topic is interesting, I keep it in mind.

…

On Tue, Oct 31, 2023, 18:09 Chun Kit LAM ***@***.***> wrote: No, we only work with the intervals without considering their relations, so 2D is similar. Just change the last x in the above example to y, the analysis just works. — Reply to this email directly, view it on GitHub <#561 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCO4MUA5BBD36YDHN4WLL3YCEWEFAVCNFSM6AAAAAA5BXFX7KVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TIMZXG42TK> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[pca006132]

No, setting the outer bounding box is actually not the main advantage of this. The most important thing is that we can skip a lot of evaluations by iteratively refine the search space, like an octree. If there are zero crossings in the current box, we divide the box into 8 and evaluate in parallel, until the boxes are small enough. It is hard to deduce the range of xyz, but given xyz we can usually get good interval of the output (for simple functions).

[gsohler]

now i understand! the equation is not evaluated for a single value, but always for intervals and there are rules how to transfer the interval from the inputs to the output. if the interval spans over 0, it's split into 8 octants, if not, the interval can be skipped. So the O(n^3) can really almost turn into o(n^2) and the result will initially look like Minecraft items, but this can resolve after smoothening. I believe this is what's actually implemented in libfive. Everything has intervals there but I did not understand until now.

…

On Tue, Oct 31, 2023 at 6:33 PM Chun Kit LAM ***@***.***> wrote: No, setting the outer bounding box is actually not the main advantage of this. The most important thing is that we can skip a lot of evaluations by iteratively refine the search space, like an octree. If there are zero crossings in the current box, we divide the box into 8 and evaluate in parallel, until the boxes are small enough. It is hard to deduce the range of xyz, but given xyz we can usually get good interval of the output (for simple functions). — Reply to this email directly, view it on GitHub <#561 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCO4MXXCQUGYKAIQEZCW4TYCEY7PAVCNFSM6AAAAAA5BXFX7KVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TIMZXHE4TC> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[gsohler]

In contrast to marching cubes, my SDF result (sphere here) does not follow an orthogonal grid, but is amorph ... [image: image.png] On Tue, Oct 31, 2023 at 8:21 PM Guenther Sohler ***@***.***> wrote:

…

now i understand! the equation is not evaluated for a single value, but always for intervals and there are rules how to transfer the interval from the inputs to the output. if the interval spans over 0, it's split into 8 octants, if not, the interval can be skipped. So the O(n^3) can really almost turn into o(n^2) and the result will initially look like Minecraft items, but this can resolve after smoothening. I believe this is what's actually implemented in libfive. Everything has intervals there but I did not understand until now. On Tue, Oct 31, 2023 at 6:33 PM Chun Kit LAM ***@***.***> wrote: > No, setting the outer bounding box is actually not the main advantage of > this. The most important thing is that we can skip a lot of evaluations by > iteratively refine the search space, like an octree. If there are zero > crossings in the current box, we divide the box into 8 and evaluate in > parallel, until the boxes are small enough. It is hard to deduce the range > of xyz, but given xyz we can usually get good interval of the output (for > simple functions). > > — > Reply to this email directly, view it on GitHub > <#561 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/ACCO4MXXCQUGYKAIQEZCW4TYCEY7PAVCNFSM6AAAAAA5BXFX7KVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TIMZXHE4TC> > . > You are receiving this because you were mentioned.Message ID: > ***@***.***> >

[pca006132]

we are unable to see your image, we can just see [image: image.png]

[gsohler]

Sorry, i have now published on imgpile The picture is visible here: https://imgpile.com/i/GUUJBg

…

On Thu, Nov 2, 2023 at 7:50 AM Chun Kit LAM ***@***.***> wrote: we are unable to see your image, we can just see [image: image.png] — Reply to this email directly, view it on GitHub <#561 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCO4MXXUMKLKQKQDGW54VDYCM7CDAVCNFSM6AAAAAA5BXFX7KVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TINJTGA3TS> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[ofmetal]

So, looks like Signed-Distance-Function here is more like Signed-Sign-Function. Distances are completely ignored.
If we create a simple sphere:

struct sdf_Sph {
	double r;
	inline static long long cnt = 0;
	sdf_Sph(double r) : r(r) {};
	double operator()(const vec3& p) const {
		double d = sqrt(p.x*p.x + p.y*p.y + p.z*p.z);
		cnt++;
		if (cnt % 1000000 <= 5) fprintf(ff,"%9lld. %.3f %.3f %.3f -> %.3f\n", cnt, p.x, p.y, p.z, r-d);
		return r - d;
	}
	string descr() const { char buf[128]; sprintf(buf, "%lld x sph(r=%.3f)", cnt, r); return string(buf); };
};

And render it like this:

	auto sph = sdf_Sph(10);
	auto m = Manifold::LevelSet(sph, { vec3(-20),vec3(20) }, 0.25);

This would call distance function 8 424 809 times, creating a sphere with 54 398 vertices.
Doubling precision like this:
auto m = Manifold::LevelSet(sph, { vec3(-20),vec3(20) }, 0.125);
calls distance function 66 462 409 times (8 times more), creating a sphere with 218 390 vertices (4x).

If we look into the sdf (well, ssf) calls we see these:

 25000000. -5.000 17.500 -13.875 -> -12.886
 25000001. -5.063 17.438 -13.813 -> -12.814
 25000002. -5.000 17.500 -13.750 -> -12.810
 25000003. -5.063 17.438 -13.688 -> -12.739
 25000004. -5.000 17.500 -13.625 -> -12.735
 25000005. -5.063 17.438 -13.563 -> -12.664
...
 65000005. 19.250 -2.875 -3.625 -> -9.798
 66000000. 19.875 -9.375 -15.750 -> -17.036
 66000001. 19.813 -9.438 -15.688 -> -16.976
 66000002. 19.875 -9.375 -15.625 -> -16.964
 66000003. 19.813 -9.438 -15.563 -> -16.903
 66000004. 19.875 -9.375 -15.500 -> -16.892
 66000005. 19.813 -9.438 -15.438 -> -16.831

At the call number 25 000 000 we see the distance result of -12.886 meaning that there are no border within this distance.
So we may safe advance 12.886/0.125 = 103 steps of the grid at once in any direction, knowing there is no border there. But the call 25 000 002 is only 1 precision step (0.125) away.

I understand, that current implementation works with any weird functions in SSF mode.
But the whole sense of SDF is Distance! Why not use this useful information for dramatically faster processing?

For the case above we'll have 66462409 / 66462409^(2/3) ~= 405 times faster processing time (0.07s vs 29s).
Which is way more faster then utilizing multiple CPU cores.

It's like using quicksort instead of bubblesort!.

[pca006132]

This is on our radar, but SDF is not receiving that much attention for now, so optimization on it is not our priority (yet).

[elalish]

If you're interesting on contributing here, PRs are most welcome! Most level set functions are not actually true distance functions, and I didn't want to get into a situation where if your function wasn't quite well-behaved enough you'd end up with a non-manifold result. We could perhaps assume the function to be a conservative bound on the distance, then do a manifoldness check and fall back to searching the space more rigorously if it fails.

[gsohler]

I just tested latest versions The SDF Example from the homepage( and probably many more) works as expected II am happy to receive bug reports at https://github.com/gsohler/openscad/ Unless I learn about the bugs you experience, I got no chance to fix them, because here in my workspace, everything "appears" to work like a charm. Latest challenge is importing STEP files into PythonSCAD (and i even intend to store META-Data next to the Faces and Indices so i can even export arcs and circles round again ) So, please feel free to announce the bugs found

…

On Wed, Nov 27, 2024 at 12:52 PM ofmetal ***@***.***> wrote: This is why i have taken the effort to use libfive in openscad. Do you mean pythonscad.org ? I tried that one couple of month ago, it was pretty buggy. — Reply to this email directly, view it on GitHub <#561 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCO4MQO6NOLRWH44ZUMBLD2CWW67AVCNFSM6AAAAABSSTSX4SVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCMZZGM4TAOA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[gsohler]

sorry, forgot my pic from what i see https://imgur.com/a/Bb0pFug On Wed, Nov 27, 2024 at 4:27 PM Guenther Sohler ***@***.***> wrote:

…

I just tested latest versions The SDF Example from the homepage( and probably many more) works as expected II am happy to receive bug reports at https://github.com/gsohler/openscad/ Unless I learn about the bugs you experience, I got no chance to fix them, because here in my workspace, everything "appears" to work like a charm. Latest challenge is importing STEP files into PythonSCAD (and i even intend to store META-Data next to the Faces and Indices so i can even export arcs and circles round again ) So, please feel free to announce the bugs found On Wed, Nov 27, 2024 at 12:52 PM ofmetal ***@***.***> wrote: > This is why i have taken the effort to use libfive in openscad. > > Do you mean pythonscad.org ? > I tried that one couple of month ago, it was pretty buggy. > > — > Reply to this email directly, view it on GitHub > <#561 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/ACCO4MQO6NOLRWH44ZUMBLD2CWW67AVCNFSM6AAAAABSSTSX4SVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCMZZGM4TAOA> > . > You are receiving this because you were mentioned.Message ID: > ***@***.***> >

[gsohler]

Hi Emmett, I get your point. Right now for PythonSCAD SDF i use libfive because it yields acceptable results. In any case I am happy to move to Manifold SDF when it becomes compatible in functionality, because there would be one dependency less (manifold is integrated already) Also libfive has the problem, that it crashes in a static constructor in windows build, so with libfive compiled in, you cannot even start up. For SDF itself do not have enough knowledge even, when i tried once to do myself. I think you set the success threshold lower for SDF, because anything which displays appears just fine. This is especially true with the Gryoid SDF which is infinitely thin by nature and only becomes viable by adding an artificial width.

…

On Wed, Nov 27, 2024 at 5:11 PM Emmett Lalish ***@***.***> wrote: If you're interesting on contributing here, PRs are most welcome! Most level set functions are not actually true distance functions, and I didn't want to get into a situation where if your function wasn't quite well-behaved enough you'd end up with a non-manifold result. We could perhaps assume the function to be a conservative bound on the distance, then do a manifoldness check and fall back to searching the space more rigorously if it fails. — Reply to this email directly, view it on GitHub <#561 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCO4MT3DZRP2QJL7LYBR3T2CXVLZAVCNFSM6AAAAABSSTSX4SVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCMZZGY3DOOA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[pca006132]

OK I looked at various sources, the most confusing thing is that SDFs can be used for both solid modeling and graphics rendering, and the requirement for them are quite different.

People write shaders for SDFs for rendering, and require the function to be a lower bound of the actual distance (e.g. with Lipschitz constant to be in (0, 1]) to make ray marching safe. People also make custom languages/libraries (libfive, fidget) to express SDFs, and use program analysis techniques to reinterpret the function in the interval domain, using subdivision to skip unnecessary work and do optimizations on-the-fly (similar to JIT, but it is just branch pruning). They also do optimization to use SIMD to evaluate the function, to speed up rasterization.

There is something orthogonal to evaluation, however: The precision of the function affects how advanced operations behave on the implicit function. For example, things like smooth union, offset, morphing, do not work for general implicit functions, and are not precise if the function is not exactly a distance function (offset may work if the internal is exact?). We can provide slightly more complicated version of union and intersection that are exact, but for things like smooth union it seems that it would be impossible to do so. This will probably affect how we should design our APIs and documentations.

---

For manifold, we only need interval evaluation and the on-the-fly branch pruning to speed up our level set construction using octree. This interval evaluation can be made safe without assuming anything about the function (continuity, smoothness etc.), unlike ray marching which requires the function to be a lower bound of the actual distance. Also, if we want gradient, we can also integrate automatic differentiation into the IR, it is a well-known technique now thanks to deep learning research.

---

Problem now: How expressive should our IR be?

In mpr for example, the opcode is quite limited, but it seems to be good enough for some simple SDFs that can be expressed in closed-form. If we want more complicated IR for things like control flow and data manipulation, things will probably get a bit more complicated (and maybe slow).

One way would be to do it in two stages, get a simple set of IR first, and then expand it as needed. Users using the simple API will have to perform function inling by themselves, and avoid iteration and using data structures in the SDF.

[hrgdavor]

that would be lovely :) . I agree with focusing on SDF strengths :) . Yes, examples :)

[pca006132]

By API, I really mean how to write the SDF, and we don't have that now. We can only compare to C++/JavaScript code, but that doesn't provide things like interval evaluation and gradient automatically. Instead, you can think of it like designing a language, because it is indeed a language, but the API will be about the low level representation instead of the surface language. Will come up with some examples tmr, probably from libfive or curv.

For our current example, I think we only have 3 relatively simple SDFs: Gyroid, CubeVoid and Layers. Gyroid and CubeVoid are simple to represent in any sensible SDF language. Layers is different: it requires conditional, but not all languages for SDFs have that, e.g. opcode for mpr doesn't have that.

[elalish]

Agreed, and I consider Layers to be much more an edge case than a compelling example. We also have spherical gyroid, which I think will fit your system better, particularly in avoiding the calculus necessary to compute the gradient currently. There's also SineSurface. I converted the Heart to a warping function, but it was originally an SDF - we should probably make an example of that too.

[elalish]

Oh, and Blobs! That's maybe the most important one - metaballs have been used quite a lot in industry and it shows a clear way to scale SDF complexity.

[pca006132]

Note that there are two things here: expressiveness and performance. For expressiveness, if you can write it as closed-form without conditionals (min/max is fine), it can be expressed with the simple set of APIs that I plan to implement first. The only issue is probably how our API will look like, but it should not affect end users anyway if they are not using C++ APIs directly. For performance, things are harder and we need an actual implementation to check that.