Extruded Polygon Side Face Normals and Other Issues

[brodieG]

Thanks!

First, thanks a bunch for implementing this feature. It really adds a new
dimension (!) to rayrender. It came in quite handy as I had it on my backburner
to implement something similar manually.

More generally, thanks for getting path tracing into R. I might never have
played this closely with it otherwise, and it's been a great value-add to my life.

Now, in using it I discovered a few possible issues. All these issues
ended up being closely related and as I tried understanding what was going on I
fixed them and they all kind of rely on each other so it is more work than I
think it's worth to split them up into separate issues (though some could be).

Also, I realize that it's bad form of me to submit a patch without first
submitting the bug report, particularly one as extensive as this one. I'm
loathe to submit bug reports without at least doing some work to figure what's
going on and one thing led to another and here we are. I realize that there is
a good chance you're already fixed much of this in your internal versions, and
that you might balk at the magnitude of this patch, but I guess that's the risk
I took by getting hooked on wanting to fix this.

I've tried to be reasonably careful about what I'm doing, and tried to test as
much as possible, but I'll admit finding flipped face normals is tricky, and
this is not my domain of expertise.

All the code I used to generate the below is in this gist.

Base Polygons are Assumed to be Closed

AFAICT this is what SF does, but many polygon applications in R will self close
polygons, including decido, graphics::polygon, and others. It might be
worth either documenting and checking that inputs comply, or alternatively as in
this patch auto-close any open polygons.

Here we compare an extruded cube (left) with a cube (right). Notice how the
back is open:

Auto-closing polygons as in the patch resolves this:

Another option is to document that polygons are required to be closed, but it's
easy enough to do this automatically.

Polygon Path Direction Affect Side Triangle Normals

With the current logic depending on which direction the polygon input in the
side polygons will end up with flipped normals. The colors indicate that the side
face normals are flipped.

Note we are now viewing the scene from the opposite direction as the previous one,
so this time the cube cube is on the left, and the extruded one on the right.

This may only be an issue with manually specified polygons, but it's easy and
cheap to check the winding with a signed area calculation so we can specify the
correct side as this patch does:

The effect is particularly noticeable in holes, showing on the left a version
hacked together with cubes, and on the right an extrusion with a hole:

After patch:

And the same thing now only extruded versions rotated so we can better see them:

After patch:

Decido Earcut Hole Index Is First Vertex of Hole

Currently, examples use the last index of the outside polygon as the starting
vertex for the holes. This happens to work when polygons are explicitly closed,
but I believe it adds an extra vertex to the hole, you end up with:

intended?        actual?
+-----+          +-----+
|     |          |\    |   <<- diagonal face is part of the hole
| +-+ |          | +-+ |
| | | |          | | | |
| +-+ |          | +-+ |
|     |          |     |
+-----+          +-----+

With closed polygons you can't see this, but if that edge ends up
crossing another hole or a convexity in the outer polygon it will become
visible. Additionally, it seems worth changing if only to match what decido
does.

I'm going off the decido vignette which I read to indicate we want the
first index in the hole polygon..

Only one Hole is Supported

Logic assumes that only one hole is present, but decido documents (and
supports) multiple holes, as does SF. The patch adds support for multiple
holes. Attempting two holes with manually specified polygons fails:

Error in data.frame(x = x, y = y, hole = holes) : 
  arguments imply differing number of rows: 15, 2

With the patch applied it works:

This will appear to work with sf polygons with holes in some cases, but it is
because the hole is rendered like so:

+----------+
|          |
| +-+--+-+ |  <<- One hole masquerading as two
| | |  | | |
| +-+  +-+ |
|          |
+----------+

Here is an image comparing after the patch (left) and before the patch (right) with
the code modified to omit the top/bottom faces so you can see inside. Notice
how the pre-patch (right) has a face connecting the two holes. This is the
same polygons used for the sf tests near the end of this PR.

This is probably fine in many cases, but could cause problems with non-convex
polygons or joining edges that show up inside holes.

Material Ids Were Implicitly Required to Be Ordered

The C++ logic assumes that material ids will show up in incrementing order. If
this is not the case as happens if material ids are passed to render_scene in
non-increasing order, materials are incorrectly re-used. Here you can see the
same color re-used:

The patch changes the logic for generating material ids to ensure they are
always ordered:

Harmonizing SF and Base Logic

The logic for flipping vertical/horizontal, earcutting, and centering is
essentially the same for either SF or base polygons. I normally wouldn't have
bothered harmonizing them because it works fine, but in this case the flipping
part was only applied to the SF derived polygons. You can see that as we
attempt to render one polygon with all the flip permutations, and all
four polygons are rendered right on top of each other:

Post-patch you see all four correctly flipped horizontally/vertically:

Related, the following was inside the poly loop when generating the faces.
Since decido always returns CCW faces this should only be needed for the
side polygons, but because this was inside the loop it will affect the face
polygons after the first loop iteration.

if(!(flip_horizontal && flip_vertical) && ((flip_horizontal || flip_vertical))) {
  reversed = !reversed
}

Because the patch adds logic to auto-detect the desired winding for the side
polygons this reversing is no longer needed so I removed it.

Multi-polygons in SF

sf objects may contain multi-polygons, where a single feature defines many
polygons. The internal logic did not seem to explicitly handle this. I didn't
dig too much into as I had to change so much to handle this as well as handle
multiple holes that it was easier to re-write. With the multi-polygon polygons
in North Carolina this is what it looks like:

With patch applied.

All of NC:

With patch:

I also built a bunch of variations of multipolygon/polygon/base to compare them.
With these the main issues seemed to be alignment (in addition to side face normals).

Various examples follow. Two simple polygons viewed from two angles, original:

With patch:

Same polygons, but this time as a single MULTIPOLYGON with two polygons:

With patch (look at the shadows in the holes):

The two POLYGONS, as well as the MULTIPOLYGON together all in one sf object:

With patch:

Other than the alignment, you can see how the face normals of the hole insides is wrong prior to patch.

Examples

I ran some of the examples to make sure those still look good with the patch
applied:

I could not run all the states together as I got a C stack error (will report
separately).

[tylermorganwall]

Hi Brodie, this looks great. It's going to take me a while to go through it due to some personal issues ongoing now, but thank you for the great work and the detailed pull request.

[brodieG]

No rush on my account. I have my fork working for what I need it to for now. Sorry you've got stuff going on and hope it resolves quickly and well.

[mdsumner]

I didn't realize you wrote this code @brodieG amazing analysis here so much to explore 👍 and @tylermorganwall I appreciate better now how you are going about this. I'm putting a note here in the hope that we can standardize approaches somewhat, hope it's of interest but no probs.

With @dcooley we've been discussing the best ways to apply this triangulation (the best way is to follow Dave's lead for having C++ converters for these things IMO), in my context it's coming from silicate where I support more formats than sp/sf but here's a barebones sf-to-triangles approach.

dcooley/geometries#4 (comment)

You can ignore the idxTRI0 class construction at the end - the list of triangles is straightforward and in feature_triangles with row indexes into coords (here coords is not de-duplicated here in the usual silicate way, not that it makes a difference given the row index).