Allowing some non-manifold geometry

[ramcdona]

Thanks for all your work on this library and for making it available.

My program has a ~20 year old home-grown CSG engine that is a bit long in the tooth. When it was pulled together, the developer didn't know about modern ideas of robust geometry etc. I am looking to replace our code with something external and modular that I can rely on.

However, our program has several good use cases for non-manifold geometry. I am curious if this library will tolerate them well enough for my purposes.

Specifically, we sometimes create surfaces that are not intended to represent a solid. Usually these are planar. Sometimes they extend beyond the solid's bounding box and are used as a cutting plane. Other times, they have finite dimension (say a disk).

We often intersect solids with the membrane. All of the normal triangle mesh-mesh intersection stuff should work as normal.

Then, we make the decision of which triangles to keep. Since the membrane has no interior or exterior, it can never contribute to the decision as to whether to keep a triangle associated with a triangle. However, we can define rules for whether to keep a membrane triangle depending on whether it is inside or outside a solid.

Sometimes we use the plane to construct a cross section cut -- say for slicing in 3D printing applications. In that case, we want to keep the triangles inside the solid. We also use this rule to construct ribs and spars and other structures inside an aircraft's wing. These are eventually fed to FEA software using shell elements.

Other times, we want to keep the part of the membrane that is outside the solid -- say when the membrane represents an actuator disk (a model of a propeller) or a wake shed by a wing.

The final challenge representing these meshes is that a single edge can have three triangles -- two for the solid and one for an attached membrane.

However, if we have two meshes...

MeshA; // Manifold solid
MeshB; // Non-watertight surface, has boundary edges, but no 'bad' triangles.

And the ability to compute the intersections (i.e. add the intersection edges and re-triangulate all cut tris)

intersect(MeshA, MeshB);

And then the ability to classify the meshes...

classify( MeshB, MeshA, MeshBinA, MeshBoutA );

Normally, one would also...

classify( MeshA, MeshB, MeshAinB, MeshAoutB );

But since B is not a solid, this has no meaning.

Normally, one would construct a Boolean operation by merging appropriate results...

merge( MeshBoutA, MeshAoutB); // Boolean union

Which will result in a watertight manifold solid.

In my case, I would want to

merge(MeshA, MeshBinA)

or

merge(MeshA, MeshBoutA)

Which would result in non-manifold meshes.

However, if I could successfully do 'intersect' and 'classify', using your library, I could handle the rest myself. I think you will find that the intersection and one-way classification do not require the manifold assumption and might be permissible with this library.

Thanks for any suggestions.

[elalish]

No, this library fundamentally can't handle non-manifold geometry, as all the intersections are based on winding numbers, which are not defined for the kinds of meshes you describe. If you want to do major surgery on your app, I'd recommend looking at our SplitByPlane method, which will allow you to e.g. cut ribs into a wing, but it will make the wing into two objects whose rib triangles face each other, thus retaining manifoldness. You can use meshID and faceID to relate the verts on each side to each other. Likewise, you can always represent a plane as a manifold by simply making a zero-thickness box (though I haven't tested that much - file an issue if you find a bug).

[ramcdona]

OpenVSP's surfaces stop short of full NURBS. We use piecewise Bezier surfaces.

What is the difference? We ditch the R, which has some tradeoffs. Most of the literature takes the R for granted, but there are some good reasons to not use rational functions.

Second, we don't bookeeep our surfaces as splines with (sometimes repeated) knots. Instead, we work with piecewise patches. This is mostly a data structure and construction artifact as the same surfaces can be represented in either form.

The bigger departure from traditional CAD is that we don't use a traditional BREP. Instead, we use the full extents of a parametric surface (we don't use trimming curves). For most objects, that surface is required to fold over in a way that it encloses a watertight volume.

While Bezier surfaces are the core in-memory representation, we go in lots of directions depending on the use case. We have a surface-surface intersection algorithm that operates on Bezier surfaces. That can generate a more traditional CAD BREP. We also have our own home-grown isotropic surface mesh generator that will generate a mesh on the trimmed Bezier surfaces.

We also have a handful of other algorithms that start with a simple wireframe evaluated from the Bezier surface and treat the geometry as discrete elements from there on.

To start with, this is how we visualize the geometry. We evaluate a simple structured wireframe mesh on each Bezier surface. That wireframe is displayed to OpenGL in either wireframe, hidden, shaded, or other view modes.

If we want to know the wetted areas of all the components in the model, we start with the wireframe (quad) mesh. I drop in diagonals to turn each object into a self-watertight triangle mesh. I perform the Boolean Union of all components in the model, discarding all interior triangles. Since all triangles are tagged with their origin component, I can then sweep through and calculated the wetted areas of the trimmed components.

I have a similar algorithm that calculates the mass properties. This algorithm uses slicing planes to track the contributions of components carried inside the aircraft -- say a fuel tank or an engine inside a wing or a fuselage. Slicing planes are constructed, exterior triangles on the slicing plane are discarded. Then, the triangles are assigned a density according to which components they slice through. That slice is stretched to a thickness, and mass properties are summed.

I have another algorithm that uses cutting planes to calculate the cross sectional area of the aircraft. When performing a linearized wave drag calculation for supersonic aircraft, you consider a cylindrical control volume far from the aircraft. When integrating over every point on this far-field cylinder, you consider the cone of influence on each point (at the Mach angle). By the time the cone of influence reaches the aircraft, it can be approximated by a flat plane. To integrate over the cylinder, we construct cutting planes (at the Mach angle) from nose-to-tail down the aircraft and also 360deg around the cylinder. This allows us to construct an equivalent body of revolution for the linearlized wave drag calculation.

Most importantly, OpenVSP is accompanied by a potential flow solver VSPAERO. Historically, VSPAERO has two modes -- a panel code mode (with thick surfaces) and a VLM mode (like AVL) where lifting surfaces are replaced with a thin mean camber surface.

Recently, we have been working on a hybrid of the two -- where fuselages and bodies would typically be represented as thick bodies and wings would be represented by their mean camber surface. We will intersect the two and trim the wing at the juncture. This generates a non-manifold geometry, but it works really well for the flow solver.

OpenVSP does all of these algorithms (and more) with a home-grown CSG engine that is over 20 years old. It was not written with modern ideas of robust geometry in mind. It uses different tolerances and heuristics all over the place. It generally works pretty well, but there are cases where it falls on its face.

One of those is when a planar surface is intersected with another object where the meshes perfectly align. This used to be a rare occurrence, but with the new mixed thick/thin VSPAERO, it is a typical case for a vertical tail intersecting a fuselage.

I am faced with the choice of adding more heuristics and here and there fixes to our own CSG algorithm, or possibly switching to a stand-alone library that has been independently designed, tested, and verified.

[hrgdavor]

@ramcdona: is there a robust alg in OpenVSP or elsewhere that handles bezier intersections in 2D that results in beziers .... or even booleans of curves enclosed by beziers that does not force them into triangulation ?

I am just curious, thought you might know, I have not found nice resources on this topic/problem.

[ramcdona]

I do not know of one. All of OpenVSP's Bezier stuff is 3D.

I do not know of much (any) literature on Robust Bezier curve/surface algorithms (where I use Robust in the formal Robust geometry sense). I will use robust in the sense of an algorithm that works pretty well in practice.

OpenVSP's surface-surface intersection is done with an AABB quadtree search. If two AABB overlap, the surfaces they contain are split in U and V and four new AABB are created. We then recurse the comparison. I think we alternate which surface is split, so we always recurse 1 into 4 rather than 1 into 16.

Once we get to a point where the two intersecting surfaces are equivalent to a plane (to some tolerance), we stop recursion, construct two triangles from the planes and intersect those tris. We calculate the XYZ and Ua,Va and Ub,Vb of the intersection line segments.

None of these algorithms are currently done in a Robust manner, but they often work. They also often fail.

I think you will find the above is a pretty standard approach to surface-surface intersection of parametric surfaces.

Perhaps more interesting is how OpenVSP handles the 3D point-surface projection problem (and several other similar algorithms).

For point-projection, the standard algorithm is a Newton's method search. Minimize the distance from the surface S(u,v) to the point P by varying (u,v). Formulate this as a distance, estimate the derivatives, run a Newton's method solver to drive the derivatives to zero. This algorithm does OK with an initial guess sufficiently close to the solution. However, the problem is plagued with local minima, and if you don't have a good initial guess, you're in a tough spot. You typically end up doing a grid search to find an initial guess to start from. If your surfaces are piecewise, you can do an AABB search to find which patches you need to start searching first.

Since OpenVSP uses only Beziers (not Rational Beziers), we gain a bunch of advantages. Primarily, Beziers are closed under addition, subtraction, multiplication by a scalar, multiplication by another Bezier, differentiation, integration, and composition. They can be split and promoted (in order) exactly. Importantly, all of these operations can be done by just operating on the control points (you don't have to resort to repeated evaluation of the curves).

So, for something like a curve-point projection calculation, we do what I call a combined geometric/numerical algorithm (where the classical approach is purely numerical).

Start with curve XYZ(u) and point Pxyz. Lets say that the curve is 3rd order.

Subtract Pxyz from the control points of XYZ(u). This gives a curve that is the displacement from P to XYZ for every u. The displacement curve is still 3rd order.

Dxyz(u) = XYZ(u)-Pxyz;

Multiply the X, Y, and Z components of the displacement curve by themselves. This gives a curve whose components are the components of the displacement squared. These curves are now 6th order.

D2xyz(u) = Dxyz(u) .* Dxyz(u); % using a Matlab-like component-by-component notation here.

Sum the components to get a 6th order scalar curve of the distance squared between the point and the curve.

D2(u) = sum_i=1^3 D2xyz(u)

Take the derivative of this curve. The derivative of a 6th order Bezier is a 5th order Bezier.

dD2(u) = d/du D2(u)

This is an analytical Bezier curve whose zeros are local minima of the point-projection problem.

Next, we use the variation diminishing property of Beziers to count how many zeros there can be in a given curve segment. We can then bisect the curve (by splitting) and recurse until we isolate curve segments that contain only one root.

We now have a relatively narrow bound (we can recurse as far as we want if it needs to be narrower) with a single root. This means that we can start Newton's method with a good guess, or we can use some other root finding algorithm that capitalizes on knowing the narrow bounds.

Up until the final root polishing step, everything has been done analytically for the whole curve -- using operations only on the control points.

OpenVSP's implementations of these kinds of algorithms are not Robust. However, in practice they are much more robust (and faster) than the traditional purely numerical approach.

It would be very interesting for someone to research the possibility of Robust algorithms like these. If you know of someone looking for a Ph.D. topic in computational geometry, this would be a good one.

It is not practical to take the same approach with Rational Beziers (or NURBS). Rational functions are closed under a similar (not quite the same) set of operations. However, I am not sure if the restriction of NURBS weights summing to 1.0 changes which operations NURBS are closed under. Even if they are similarly closed for a useful set of operations, these operations cause a much more severe explosion in order of the curves for Rationals. Every time you need to find a compatible denominator, you go to a higher order. Even operations like differentiation (which reduces the order of a Bezier) increases the order of a Rational Bezier. Consequently, this kind of approach is impractical for NURBS curves and surfaces.

[elalish]

Thanks for the detailed explanations - I certainly learned a few things. When you say you "use the full extents of a parametric surface", does that mean you are basically operating on a CSG tree, vaguely akin to OpenSCAD (except using bezier surfaces)? That certainly seems amenable to the use of this library, once those beziers are triangulated. As for cross-sections, I wonder if #426 would help? It would be helpful to have a power-user involved to make sure the API design for that is good, since I imagine you may need relations.

[ramcdona]

By "use the full extents of a parametric surface" I mean that the domain of a given surface is U in [0,1] and V in [0,1]. OpenVSP uses that entire domain. The boundaries of a given surface will always be U=0, U=1, V=0, and V=1.

This is contrasted with the standard BREP approach that most CAD programs take. In a CAD-BREP approach, there will likely be several curves defined on a given surface. These curves are used to trim the surface and consequently define the boundaries of a face in the BREP. I.e. the curves define a subdomain. There are regions of the surface outside these curves that therefore go "un-used".

OpenVSP requires each of our components to be individually watertight (with some exceptions of course). This typically means that the V=0 and V=1 edges of a surface must coincide -- and that the U=0 edge will either fold over itself or will collapse to a point. The U=1 edge will behave similarly.

Multiple (individually watertight) components are arranged (usually overlapping). This implies a CSG union of the components. We do not visualize or compute this CSG Union on the fly.

When needed, we compute the CSG Union of the model either using our built-in triangle soup CSG engine or our built-in Bezier surface intersection code.

We do not currently support a full CSG operation tree. For a long time, we only supported the Union of all components. More recently (still a long time ago) we added support for 'negative components'. These imply a CSG tree with one branch constructing the union of all normal components -- the other branch with the union of all negative components -- and a subtraction operation at to combine these two branches together.

So far, our users and sponsors have not prioritized implementation of a comprehensive CSG tree -- all the operations, but more importantly, the ability to specify them in arbitrary orders and combinations. Our geometry engines are pretty much ready for this, but integrating the complexity of specifying a CSG tree (when our model is already assembled with a different tree construct) presents a UI challenge.

I do really like the idea of investigating a dedicated planar slice computing capability. It makes sense that it could be faster and more robust than using a general purpose CSG engine for the task. However, I will likely keep that on the back burner as a potential down the road optimization. For now, I need to make sure I can handle normal CSG and also our somewhat odd thin-surface CSG (either via manifold thin surfaces or some other approach).

[ramcdona]

I am working on converting my non-manifold thin surfaces into manifolds (duplicate mesh, swap triangle node order, swap normal vector direction, then merge boundary points and connect front/back mesh together).

I am creating a MeshGL and am supplying unique faceID for every triangle (this helped). I am also assigning a unique (integer) property to every vertex (this didn't help).

When I create a manifold from my MeshGL, it moves two vertices to apparently lose some triangles. However, because the faceID's are unique, the tris are still there, they are simply collapsed in space.

The blue line is the mesh I start with. The gold is what I get out of the manifold immediately after creation.

How can I keep manifold from collapsing those tris?

The blue is a thin manifold surface. I.e. all of the perimeter vertices for the front and back copy are shared, but the interior six vertices on the front/back sides are independent.

So, the two triangles that vanish have three shared vertices and two shared edges. The third edges should be unique, but they actually connect identical verts, so the only way to distinguish the top/bottom edges of those tris is by which faceID's they pair up to.

The MeshGL data structure has mergeFromVert and mergeToVert, which makes me think I could use this to make the thin surface manifold (instead of my DIY approach). I.e. I could send the non-manifold mesh twice (with verts and tris in identical order), and then tell it to merge the verts around the perimeter to make a manifold.

That said, I think this would put me in the same position.

[elalish]

Hmm, that's surprising. You can go ahead an file an issue on it with a simple repro. I'd guess something funny is going on in SimplifyTopology() if you want to take a look yourself. We certainly don't have a test case on this yet.

[ramcdona]

What file format is preferred? I haven't built assimp yet -- looks like I don't have a choice now.

I see the GLTF writer in native JavaScript -- is there a lossless file I/O capability native to the C++ library?

[pca006132]

No, we don't have any import/export functions. If you want to avoid using any additional libraries, https://paulbourke.net/dataformats/ply/ should be a very simple format for dumping out the data.

[ramcdona]

OK, thanks. I've got some simple ones myself, but I was looking for something that would carry forward all the properties and faceID etc. to be sure my problem isn't someplace non-obvious.

[ramcdona]

Here is a sample *.ply file of the information I'm putting into a MeshGL to construct a Manifold.

This loses the two 'corner' triangles as discussed before.
surface0.ply.zip