Working with lists inside distance functions

[p-e-w]

I am trying to build a library for working with polyhedra in Curv. I have written a function that computes the signed distance from a point to an arbitrary (convex or non-convex) polyhedron, described by a list of vertices and a list of faces, each of which consists of a list of vertex indices.

The function is working, in the sense that it gives the correct result when invoked separately, but I am hitting a brick wall trying to use it inside an actual shape's distance function. The problems all seem to come down to limitations when working with lists inside distance functions.

For instance, it does not seem to be possible to even index a list:

let
points = [[0, 0, 0]];
in
make_shape {
  dist (x, y, z, t) =
    mag points[0];    // <-- "Geometry Compiler: not a constant"
  bbox = [[0, 0, 0], [1, 1, 1]];
  is_3d = true;
}

Surrounding point[0] with parentheses gives a different error, which makes me suspect there is a bug here as the two expressions should be entirely equivalent:

let
points = [[0, 0, 0]];
in
make_shape {
  dist (x, y, z, t) =
    mag (points[0]);    // <-- "Geometry Compiler: [[0,0,0]] is not a function"
  bbox = [[0, 0, 0], [1, 1, 1]];
  is_3d = true;
}

Attempting to move the declaration of points inside the distance function reveals yet another problem:

make_shape {
  dist (x, y, z, t) =
    let
    points = [[0, 0, 0]];    // <-- "this list constructor does not support the Geometry Compiler"
    in
    mag points[0];
  bbox = [[0, 0, 0], [1, 1, 1]];
  is_3d = true;
}

The documentation does mention that standard iteration is not supported within distance functions and that imperative constructs are needed instead, so I was prepared to have to work around that. But from the above examples it appears that not only iteration, but lists in general are unsupported, as neither indexing nor construction appears to work. Is the Geometry Compiler really this limited?

Looking at the standard library, all existing distance functions seem to use either the built-in min/max primitives, which work with list arguments, or be constructed by repeated application of some fixed n-ary operator, allowing the handling of lists to be moved outside of the distance function. Neither of the two approaches will work for arbitrary nonconvex polyhedra with nonconvex faces.

I hope I am missing something here and that there is in fact a way to work with lists that has simply escaped my notice so far. If there isn't, it would unfortunately seem that many very interesting distance functions are currently out of reach :worried:

[doug-moen]

This is a limitation of the Geometry Compiler that needs to be fixed before Curv reaches 1.0 status. Because you are right, many interesting distance functions are currently out of reach.

Currently, the only list types supported are "vec2", "vec3", and "vec4", which are lists of 2, 3 and 4 numbers. The only escape hatch right now is that you can export the shape as a mesh, from the command line, without using the geometry compiler.

The next step forward is to add support for multidimensional, rectangular arrays of numbers.

[p-e-w]

Thanks for the suggestion, mesh export indeed works and for the first time I can see polyhedra!

Unfortunately, the quality of the mesh is very poor. There are artifacts surrounding the edges of the polyhedron that get smaller when vsize is reduced, but don't disappear completely. Even a simple cube, when rendered as a multi-megabyte mesh, looks nowhere near as crisp as objects tend to do in the geometry viewer :disappointed:

[doug-moen]

I've been thinking about implementation. It would be easier, as a first step, to add support for lists of numbers, and lists of points. General multidimensional arrays are more complicated, since GLSL 1.5 (the GPU language I am targeting) does not support these.

The default mesh generation algorithm is a variant of marching cubes, which works very well for the curved, organic looking shapes that I make. However, it works poorly for polyhedra. Right now, your best bet is to create a high resolution mesh using -Ovsize=, and then simplifying the mesh using MeshLab, using Quadric Edge Collapse Decimation and Planar Simplification.

That's just a workaround. On the 1.0 roadmap, I want to support multiple meshing algorithms, depending on the kind of model you are creating. For polyhedra, I want a meshing algorithm that detects sharp edges. I've tried Dual Contouring: it works, but the output contains self intersections, which causes problems for importing into other programs (eg, you can't import these meshes into OpenSCAD). @mkeeter's libfive had the same problem a while back, but last month Matt introduced a new meshing algorithm that does edge detection and produces valid meshes. I just haven't got around to trying it yet.

Other meshing algorithms on the roadmap: (1) Support for multi-material FDM printers, since I'm upgrading my Prusa with a multi-material unit next year. The idea is to automatically export one mesh for each colour or material. (2) Support for fractals like the Mandelbulb. This requires a specialized algorithm for good results.

[s-ol]

@doug-moen I don't remember whether I brought this up on the mailing list before, but have you considered targeting SPIR-V instead of GLSL? scopes is a language that targets x86 and SPIR-V and it seems that it can manage to compile some rather complex programs with a variety of types, although I can't say whether this is due to an advanced compiler or because SPIR-V leaves more room for that.

[doug-moen]

SPIR-V is part of Vulkan, which is the successor to OpenGL 4.6. GLSL is part of OpenGL. There's no important difference between SPIR-V and GLSL 4.6 in terms of my roadmap. At this point, I'd rather target OpenGL, because it has a larger ecosystem, and there are more open source projects that I can leverage. OpenGL is also better if I want to eventually run Curv in a web browser: Vulkan will never run in a web browser.

The GL compiler limitations are mostly my fault, not a problem with GLSL per se. The GL compiler was designed at a time when I was trying to duplicate the features of ImplicitCAD and Antimony -- at that time, I knew little about F-REP, and I concluded that I only needed to support arithmetic expressions. Then I discovered Shadertoy.com, and I added while statements, mutable variables, etc, so I could do fractals, etc. Now, the original compiler architecture is outdated, and needs to be rewritten.

The GLSL limitations come from the version of OpenGL that I choose to target. Right now, I target OpenGL 3.2, which runs on a lot of older hardware. To implement everything on my roadmap, I eventually need to target OpenGL 4.3, which requires a GPU built in 2012 or later, but is not supported by MacOS.

Thank you for telling me about scopes. There is a lot of overlap with Curv, and it's interesting to look and see what design decisions were made, in comparison to Curv. Scopes is for writing video game engines, and the target users are video game engine developers. From what I've read, Scopes pushes more complexity onto the user than what I want for Curv, so I probably can't use exactly the same approach for the compiler.

[p-e-w]

In the context of constructing polyhedra, one way or another you will need a list of faces (or vertices, or edges, depending on the representation chosen), each of which takes the form of a list of vertices (or faces, or edges). List elements can be either 2D/3D points, or integer indices into other lists.

Therefore, constructing arbitrary polyhedra requires at minimum support for lists, and lists of lists, of num, vec2, and vec3. Looking at my implementation, this would indeed already suffice. Non-rectangular lists, inhomogeneous lists and lists of arbitrary nesting depth are not necessary for this and a large variety of other interesting problems.

[TLC123]

There is an general issue with constructing sharp meshes, regardless the choice of algorithm. If the underlying distance field is not continuous and reasonably linear there is no guarantee of even a fair approximation .

[doug-moen]

@TLC123 said "There is an general issue with constructing sharp meshes, regardless the choice of algorithm. If the underlying distance field is not continuous and reasonably linear there is no guarantee of even a fair approximation."

That is correct. Some functions are too ill behaved to be considered distance fields, and cannot be rendered at all. Some distance fields can't be rendered by the forthcoming "sharp" mesh generator with edge detection, but still meet the requirements for sphere tracing: these can be rendered by the forthcoming "fractal" mesh export algorithm.

My initial plan is to support 3 separate meshing algorithms, and let the user select the best algorithm for the job. It would be better, from a usability perspective, to just have one, but I don't know if that will turn out to be feasible. I'll have more information once I have all 3 algorithms available for testing.

[doug-moen]

@p-e-w said "Therefore, constructing arbitrary polyhedra requires at minimum support for lists, and lists of lists, of num, vec2, and vec3."

Thanks for the advice. I'm working on it. So far, I have implemented lists of num, vec2, vec3, vec4. Two dimensional arrays of num, vec2, vec3, vec4 are not ready yet.

[doug-moen]

What I have implemented is good enough for small arrays, maybe up to 16K of data? It depends on the GPU. Large arrays need a different approach: I'd probably embed the array data in an image, and synthesize the image in the Geometry Compiler. Large arrays will require a separate feature request.

Here is my test case: https://github.com/doug-moen/curv/blob/master/examples/tests/gdf.curv

[p-e-w]

@doug-moen Thanks for the update, I'm experimenting with it now. It appears that list comprehensions are still not supported, but I'll try to work around that with while constructs.

Meanwhile, I have solved the mesh generation issues by writing a custom STL exporter for polyhedra in Curv. It has direct access to the geometry information and thus produces "perfect" meshes that exactly match the actual polyhedron using the theoretical minimum number of triangular facets (by triangulating the faces of the polyhedron).

[doug-moen]

Within a distance function, arrays are not first class values, they are required to be constants. You can bind a constant array value to a variable A using let, you can index the array using A[i] or A[i,j], and that's it. This means that array contents cannot depend on the (x,y,z,t) parameters of a distance function. If you are okay with this restriction, then you can use list comprehensions if you lift the list comprehension out of the distance function and bind it to a local variable using a let that encloses the distance function or shape expression. Building an array inside a distance function using a while loop is not supported right now. You can build the array any way you want if you do the work outside of the distance function.

I am planning to reimplement the Geometry Compiler using a better IR (Intermediate Representation), before the 1.0 release, and then some of these restrictions will be removed. The body of a distance function will be partially evaluated, which will allow you to use a larger subset of the Curv language.

[doug-moen]

@p-e-w I'd love to see your code if you are comfortable with sharing it.

[p-e-w]

@doug-moen Good to know. This complicates the code a bit, but can be worked around. I have one more request: Could you add support for the count function to the Geometry Compiler? This isn't a show stopper, but it would be nice not to have to pass lists of array lengths around everywhere just to be able to iterate.

I definitely intend to publish the code, I just want to experiment and polish for a while longer (and, hopefully, actually make it work :smile:...).

[doug-moen]

@p-e-w: "This complicates the code a bit, but can be worked around."

I added count(list) to the geometry compiler. What other features are missing that you need for this program?

[p-e-w]

@doug-moen It looks like I'm still one level of indexing short. I am working with a list of face polygons, each of which is a list of vec3 vertices, and for one algorithm I need to access individual coordinates of vertex points.

This doesn't work:

let
faces = [
  [
    [0,0,0],
  ],
];
in
make_shape {
  dist(x,y,z,t) =
    let
    face = faces[0];
    vertex = face[0];
    d = vertex[X];     <-- Error
    in
    d;
  is_3d = true;
}

I was expecting this as you mentioned that only two levels of indexing are currently supported, and I thought it wouldn't be a big deal because coordinate access can be simulated by multiplication with a suitable projection matrix. But dot doesn't work with vertex either, so I'm stuck.

How can I get the individual components of vertex in this case?

[doug-moen]

It doesn't work, but the error message is wrong. So I fixed the error message. It now says:

ERROR: Geometry Compiler: at field .dist: can't index a 2D array of vec3 with a single index
at file "A":
11|     face = faces[0];
                     ^  

You can't say faces[i], you need to say faces[i,j]. This limitation reflects limitations in the GLSL language. GLSL is a C-like language. In C, 'face = faces[0]' would be implemented by defining 'face' to be a pointer to an array of vec3, but GLSL does not support pointers.

[doug-moen]

You can now write faces[i,j,k] and get back a number.

Matrix multiplication is not yet supported by the Geometry Compiler. I consider this important; I just haven't got around to it yet.

[doug-moen]

GLSL doesn't have pointers, but you can simulate them with integers. A common polyhedron representation is two arrays. vertexes is a list of vertexes, and faces is a list of lists of vertex numbers (which are indexes into the vertexes array). If the faces are triangles, then faces is a list of vec3. Then, faces[i] is a vec3, which you can store in a variable.

[p-e-w]

Thank you for the update, the error messages are indeed much clearer now.

But now that I see what the actual problem was, the count issue is back, unfortunately. If I have a two-dimensional faces array and I can't index it with faces[i], how can I determine the number of vertices in a face, since count(faces[i]) will not work?

Perhaps a two-argument function icount could be added, with icount(list, [i, j, ...]) being equal to count(list[i, j, ...]), but without the need to evaluate list[i, j, ...]?

[doug-moen]

I can fix this. I am thinking of two alternatives. Once I dive into the code, I'll see which is more practical.

1. Modify the compiler so that count(faces[i]) works as a special case. Note that since Curv only supports rectangular arrays on the GPU, the value of i is irrelevant in count(faces[i]).
2. Add a new operation: dim(a) returns a list of the dimensions of multidimensional array a. In your case, you would use dim(faces)[1].

I'm curious about your code. Since Curv only supports rectangular arrays on the GPU, your data structure requires all faces in the faces array to have the same number of vertices. To represent a polyhedron like the cuboctahedron, which has both square and triangular faces, you would need a different data structure. Am I understanding this correctly?

[p-e-w]

I have worked around this issue by determining array lengths outside of the distance function and passing them around. The other issue you describe is now easily fixed by padding all face arrays to the same length.

At long last, my code compiles now. But the rendering in the preview window is utterly broken. I have reason to believe that this is a bug in Curv as the distance function for a cube constructed with my code seems to return exactly the same values as that of the cube from the standard library, when evaluated without rendering.

I will now write some more tests to confirm this, then publish my code and write a proper bug report.