Surface moments

[devonmpowell]

No description provided.

[devonmpowell]

@devonmpowell @raovgarimella @shevitz
I opened this thread for the discussion on 2D moments

[MackKenamond]

(copied post from elsewhere)
We are finding that we'll need to know moments of the clip face. Given a clip plane and polytope, clip the polytope with the plane and return the moments of integration over the clip face, maybe in addition to the remainder polytope and its moments. There are several models that could use those area moments. For now, my best workaround solution is to clip the polytope, shift the clip plane in its normal direction by epsilon and clip the remainder polytope with the shifted plane. Then use that remainder polytope (thin cheese slice) to deduce a "good enough" set of moments for the clip plane (divide the slice volume by epsilon to get the clip face area, divide the int_V{xdV} by epsilon to get the area centroid. Etc. It's not ideal, though.

[raovgarimella]

(copied from PR thread)
So if the algorithm returns the clip polygon (I am sure it knows this), you could rotate that polygon into the XY plane and compute the moments (either by converting it to an R2D poly and doing r2d_reduce or independently using divergence theorem).

[MackKenamond]

Or that. For SSD closure, we need the area, but for some other remapping work, we could use the higher moments to compute things at higher order accuracy (in contour integrals).

[MackKenamond]

And yes, if r3d returned a 2D polygon (in 3D coordinates), then we could transform it to 2D and use r2d to compute moments. That works for me. Similar analogy for clip segment for polygon clipping - return the line segment in 2D and use r1d to compute moments...oh wait! there's no r1d ;)

[MackKenamond]

One issue for 2D clipping is that you'd somehow need to return multiple segments for non-convex input polygons. In 2D, there can be disjoint polygons connected by a line segment. All the moments are correct. Not sure how to connect two disjoint line segments in one usable data structure.

[devonmpowell]

I think it's worth taking a look at Koehl 2012 (https://ieeexplore.ieee.org/document/6127882). This is the recursive algorithm used in r3d_reduce.

An intermediate step in the derivation is to obtain the surface moments (as a sum of triangles). Then it is extended to 3d by adding the origin as another vertex to make tetrahedra. So it's probably possible to get the 2D moments (embedded in 3D) with a few cleverly-placed lines of code. Then we can just use the face normal to project the moments onto the face itself, without any intermediate r2d step.

[raovgarimella]

The issue I see is that r3d_poly has no notion of a face unlike r3d_brep. So it is impossible to tell r3d_reduce which face to get moments for. On the other hand, we do have that info when we are clipping or splitting and we throw it away.

[MackKenamond]

Maybe this ends up being tied to the non-trivalent r3d option and the suggestion to use half-edge data structures. Each half-edge is associated with two faces connected by that half-edge: a previous face and next face. Face ordering could be explicitly defined as part of the half-edge or it might be defined implicitly based on a greedy approach and the half-edge ordering. Half-edge 0 implicitly defines face indexes 0 and 1. Or maybe it only defines face index 0 (it's previous face but not its next).
A similar concept would still apply for the tri-valent representation. Each vertex bounds three faces: between its vertex 0 and 1, between 1 and 2 and between 2 and 0. However, if degenerate edges and replicated vertexes are introduced for non-tri-valent cases, there can be a degenerate face thrown into the mix.
Whatever the approach, if an edge is clipped, the two faces connected by that edge are also clipped but can maintain their old face id.
The new faces generated from the clip pose a problem, however. Although there is one planar "clip face" introduced with associated moments, it can be disjoint and topologically represents multiple faces. Perhaps (probably) there is another implicit numbering scheme for these multiple new faces that the caller could rely on to identify these faces, but things seem to be getting uglier at this point.
A problem I see for implicit face id's is that the remainder polytope from a clipping operation will likely produce face id's that are inconsistent with the input polytope. Maybe there's some clever implicit ordering scheme that preserves the face id's. Explicit ordering avoids this issue, but it adds to the memory footprint.

[MackKenamond]

Another thought on defining face ids: The input must include one vertex or half-edge associated with each face. For tri-valent polytopes, you have to give one vertex that has the face between its 0th and 1st nodes. For half-edge, you have to give one half-edge that has the face as its "previous" face. Less input data required that way. If you had to define faces for every vertex or half-edge, there would be a lot of redundant data.

[devonmpowell]

I think the half-edge is worth looking into.
I am tinkering with a "son of r3d" that uses the half-edge structures, partly for my own intellectual satisfaction and partly because I think it will make a lot of operations easier for you guys.
I should be able to finish a quick version in the next couple of days, then maybe you all can take a look and see what you think?

[raovgarimella]

Is there some advantage that to Half-edge data structures I am missing entirely? Maybe we are talking of entirely different data structures. I knew half-edge data structures from the place they originated (Weiler and Shephard). Each half-edge or edge use only knows about its first vertex use, second vertex use, its previous edge, its next edge and its opposite half edge. Since edge uses are defined by vertex uses, you would have as many vertex uses for a vertex as there are faces coming together at the vertex. Each face is described by a set of ordered half edges. Each solid is described by a set of faces (for non-manifold solids you can actually have half-faces or face uses). The structure was really meant for modeling non-manifold geometric models rather than mesh elements.

So, you lose the trivalent vertex structure you have now and you make it a lot more verbose than before. You will be jumping around in memory a lot more and apps will once again spend time converting from BRep to the half-edge data structure and back (same or more cost as init_r3dpoly right now) - this is nearly half of our cost of usage of .R3D.

My strong suggestion if you are making the effort to rework the code, work on a straight up BRep data structure. We do it all the time in Portage, Tangram. I use it in MSTK and Jali. It is the most compact representation of a polyhedron that has a notion of faces and vertices. (I have reach far back into my brain, but perhaps FLAG MPygs are also the same?). I would've loved to do this but can't right now.

The simple BRep structure that's already in the R3D code is given below. It uses dynamic arrays for now but if you want, you can convert it to use static arrays. If you need to, you can devise one that explicitly represents edges.

typedef struct r3d_brep
{
  r3d_rvec3 *vertices;
  r3d_int numvertices;
  r3d_int **faceinds;
  r3d_int *numvertsperface;
  r3d_int numfaces;
} r3d_brep;

The static version can look like

typedef struct r3d_brep
{
  r3d_rvec3 vertices[R3D_MAX_VERTS];
  r3d_int numvertices;
  r3d_int faceverts[R3D_MAX_FACES][R3D_MAX_VERTS_PER_FACE];
  r3d_int numvertsperface[R3D_MAX_FACES];
  r3d_int numfaces;
} r3d_brep;

R3D_MAX_FACES can be same as R3D_MAX_VERTS (since you won't be proliferating vertices by enforcing the trivalent condition)

[raovgarimella]

As for the original topic of surface moments, I interpreted Mack's original query as wanting the moments of the clip face or the face between the two halves. My suggestion there is to hew close to what I floated earlier. Let the R3D clip and split routines just give you a R2D poly that you can call R2D_reduce on.

[devonmpowell]

Hmm, this is definitely true about jumping around in memory. The half-edge is more of a linked list than anything else.
I see two difficulties with the brep approach. One is that it would take some extra looping to identify "edges" (pairs of vertex indices on neighboring faces) with one another, or even neighboring faces for that matter. The other is the question of, for a non-convex poly, do we treat disjoint faces created by the same clip operation as the same face?

One other questions about projecting faces into 2D is, does it matter what the x-y orientation is on the original face?

[MackKenamond]

I vote for Devon trying out the half-edge approach and maybe comparing performance to the current version. It sounds like he's been thinking about it and is ready to go for it. There are issues with the current version and there are at least two proposed solutions, each certainly with pro's and con's. Until/unless someone tries both out or does a methodical on-paper evaluation of pro's and con's, we won't know which solution is superior. My guess is that both solutions will be slower and/or use more memory in certain cases and a "winner" will be subjective (for example, I don't care about an increase to R3D_MAX_VERTS).
A solution that enables face moments is going to be a requirement for us in the not-too-distant future (~6 months) and it seems like a logical extension of the current r3d. Having a slower solution is infinitely better than not having a solution. We'll probably use the thin sliver option until a proper solution is available. It requires twice as many clipping operations, but it provides the solution we'll need.
As for the clip face and moments, I think that having a "single" (but potentially disjoint) face (or moments for that face) would work for us in the short term, but being able to extract face moments for all polytope faces would be the more generally-useful option. There obviously needs to be some way for the caller to identify individual faces, but I think both brep and half-edge solutions can offer this. If an input polytope includes face indexes, any surviving face ID can remain unchanged and any clip faces (plural) could have IDs greater than the total number of input faces, so the caller can still identify all faces after a clip.
I agree that converting to a 2D polygon poses challenges and I don't even think it's necessary. If you have the 3D moments, the caller can do whatever is needed with them. We'd require a new function that returns face moments given an input polytope, and (again), there has to be some way for the caller to identify each face.

[MackKenamond]

I'd also like to ask a philosophical question about r3d. On one hand, every option could be implemented via some function. R3d could support trivalent, half-edge and brep. It could offer dynamic and static allocation options. On the other hand, one "compromise" set of functions could be provided that is sufficiently general. The first option is harder to maintain has a more complicated API but is completely flexible and probably offers the best possible performance in the most use cases. The other option is easier to maintain, is simpler to use, but will probably be the least performant. I'm curious what your opinions are on this.

[MackKenamond]

FYI, we've ported r2d/r3d to modern fortran for our project (so that we can easily pass around polytopes on the fortran side of things in our fortran host code and get things to perform on GPUs) and have done a lot of extending. We added a non-trivalent type (which is 2.5x faster on our polytopes) in addition to tri-valent (can use either). We've added the ability to optionally return the clip face that results from clipping (r3d_clip optionally returns an r3d_polygon type and r2d_clip returns an r2d_line type). We've added the ability to compute moments on these degenerate types (moments of r3d_polygon and r2d_line types). We will probably also add the ability to linearly interpolate polytope vertex fields (interpolate edge endpoint values to the clip vertex as each edge is clipped). There was a bit of debate about how to implement all of this. I'm guessing it would be a similar situation for C or C++.

[devonmpowell]

Very nice, Mack.
For the field interpolation, it may be useful to just add a generic interpolation field to the polytope struct. Something like a set of barycentric coordinates for tetrahedral cells, which are automatically interpolated as you clip. For hexahedral cells, you could use some field of 8 numbers (which range from 0 to 1 inside the cell) and then trilinearly interpolate
Then you can easily map physical quantities onto it after the geometric operations are done.
Let me know if that makes sense...

[MackKenamond]

Because all new vertexes are added on the perimeter of the polytope and always at edges, I believe that edge interpolation would produce the same result as tetrahedralization and linear barycentric interpolation. As for hexes, you could certainly do something more exotic/accurate. Currently, I'm just thinking about using this for visualization purposes where we have clipped-up a multi-material cell into pure sub-pieces and want to provide smooth-looking solutions near material interfaces, so it's more qualitative. Any actual remapping we do relies on exact volume integrals via moments.
Just got done writing up r1d just so our 1D/2D/3D code can use consistent algorithms in all D.

[devonmpowell]

Cool! It sounds totally sufficient for your purposes.