Roofs, Bevels, and Skeletons: Introducing the Raybevel Package

Data Visualization

R

CGAL

Mesh Processing

GIS

Rayshader

Rayrender

Rayvertex

Straight Skeletons

Package Development

Author

Tyler Morgan-Wall

Intro

It all started with a tweet that caught my eye - a slick 3D rendering of beveled polygons by Yi Shen @pissang1.

A couple things nerd-sniped me about this interaction: first, an algorithm I’ve never heard of? Check. One that can be used to create cool 3D meshes? Check again. And (after a few minutes of research) one that had no existing implementation in R? Check++.

The more I read, the more I realized that this would be a great new capability to add to R–and for more than just my rather niche interest in the 3D meshing and fancy beveling aspect. But first: what is a straight skeleton, and why is it useful?

What is a straight skeleton?

#Load all libraries used in this plotlibrary(spData)library(ggplot2)library(sf)library(av)library(patchwork)library(rnaturalearth)#All the rayverse packageslibrary(rayrender)library(rayvertex)library(raybevel)#Helper function: Center mesh in x/z coordinatescenter_mesh_xz =function(x) { mesh_bbox =get_mesh_center(x) mesh_bbox[2] =0translate_mesh(x,-mesh_bbox)}pennsylvania = spData::us_states[spData::us_states$NAME =="Pennsylvania",]pa_skeleton =skeletonize(pennsylvania)plot_skeleton(pa_skeleton)

1

Load the non-rayverse packages

2

Load the rayverse packages

3

Create a function used throughout the function to center meshes

4

Extract polygon of Pennsylvania from the us_states dataset

5

Skeletonize the polygon

6

Plot the skeleton (using raybevel, the package introduced in this post)

The straight skeleton of a polygon is the geometric object formed by tracing the vertices of a polygon as the edges propagate inwards. You can visualize this as a interior wavefront formed by all the edges of the polygon: vertices occur at places two edge’s wavefronts combine or split.

Get the maximum internal distance out of the straight skeleton structure

2

Create offsets for these, except for the first and last values

3

Create frames of animations

4

Create animation mp4 file

Straight Skeleton Applications

It’s useful for both 3D and 2D applications: one useful thing this algorithm can do create inset polygons that respect the topological and geometric properties of the original shape. These interior polygons preserve sharp corners, in contrast to st_buffer() which generates curved internal polygons (see Figure 4).

This is a much better option if you’ve ever resorted to negative buffers for manipulating polygons for aesthetic reasons: this package allows you to generate smaller polygons without introducing curves (for example, you can use this to create small gaps between directly-adjacent polygons without resorting to hacks with the border color).

This polygon shrinking process can also be represented continuously in 3D. Here’s where one of the really cool applications comes into play: 3D rooftops!

You can also use the straight skeleton to obtain internal polygon distance data and generate arbitrary bevels (the original goal of all this effort!).

There have also been other creative uses for the straight skeleton in GIS, such as road network simplification and estimating water drainage paths from height contours.

Developing the {raybevel} package

By itself, generating 3D rooftops was a strong enough motivator for me to add this capability to R and the rayverse: plain flat extruded polygons never had enough “character” to represent something as human and familiar as a house. It also serves as a data visualization cue: it’s nice to have a visual indicator (in this case, a sloping rooftop) that helps differentiate objects representing abstract 3D data from 3D buildings. But primarily, it would help provide a bit of humanity to visualizations of neighborhoods, which are plentiful in mapping.

Weighing the costs of complex dependencies

When I approach building a new package, I always start by doing some basic research: did anyone else already solve this in a header library I can wrap? Or is there an existing implementation in another language I can port to R? Straight skeleton generation is indeed available in the CGAL library, which has headers available via the RcppCGAL package. The Computational Geometry Algorithms Library (CGAL) is a robust, open-source software project that provides easy access to efficient and reliable geometric algorithms in the form of a C++ library. It encompasses a wide range of tools for computational geometry tasks, from basic geometric primitives to complex operations like mesh processing and geometry analysis. That power comes at a cost: the library is huge and complex, and due to its size requires a non-CRAN distribution method for the actual header files. And one thing I’ve learned over years of R package development is the danger of depending on another package for compiled code, particularly if it’s a non-standard installation process. The CRAN policies are fairly static, but the CRAN’s tooling is not: packages that are fine for years can suddenly have a two week deadline for removal looming over them because the latest version of GCC the CRAN uses is now throwing warnings.

Note

This just happened recently with the CRAN throwing -Wformat errors for the progress package, which propagated a 2-week removal warning to one of my packages–and by the time the progress package was uploaded and fixed, I only had a few days to issue my change! Additionally, the RcppCGAL package was taken down for several months this past year, and I worked with the maintainer to get it back up on the CRAN.

One of my choices as a solo R developer creating a universe of packages while also raising a kid, maintaining a house (e.g. surprise! your fridge is dead! hope you had no other weekend plans!), and working full time is to prioritize limiting surprise workloads with a deadline. It’s becoming more rare that I can drop everything and figure out a potentially complex workaround for a critical software dependency, so I try to only depend on my own universe of packages.

First attempts at a homegrown implementation

So I wanted to try to forgo the heavy CGAL dependency first, and luckily my research indicated that there were alternatives. Several people had seemed to have successfully created implementations of the Felkel and Obdržálek 1998 straight skeleton algorithm (polyskel and ladybug geometry). In particular, I saw that someone had adapted it for a plug-in for blender for creating 3D bpypolyskel and had done some testing showing robust performance (successful on 99.99% of OpenStreetMap building footprints). This algorithm was an attractive choice due to it’s relative simplicity–it only required implementing a circular linked-list data structure with some processing functions, and otherwise was a fairly small library. I decided to implement it in C++ for both speed and ease of translation, as working with circular reference-based data structures isn’t very natural in R. And after brushing up on my python and diving in, I finally had a working algorithm!

Well, sort of. My data suggested the 99.99% performance number was a bit optimistic . My own numbers indicated that the algorithm failed significantly more frequently.

the nice results from the package testing on OSM buildings was likely due to most buildings having a fairly simple and similar structure

Here, the “correctness” of the straight skeleton structure depends on minute differences between nodes, and even with double precision I ran into issues.

Note

Part of the difficulty comes from the fact that the straight skeleton is inherently a non-local structure: tiny changes to the inputs can lead to wildly different skeleton configurations. Ensuring the resulting nodes are purely hierarchical (e.g. no loops or crossing links) is critical to the meshing process.

Poor results and lessons learned

It turns out the algorithm failed more frequently than I’d liked. I knew it wasn’t perfect going in, but the frequency of the failures (a rate closer to 1%-0.1% than 0.01% in my experience, and made worse as polygon complexity increased) made it pretty unusable for large-scale geographic visualizations. But it was not all wasted time–first of all, I got to brush up on some rusty python skills by translating an entire package with relatively complex data structures to C++! But more importantly, I had a much more intimate understanding of the straight skeleton structure, which was helpful when implementing 3D mesh/bevel generation. And with that understanding came the hard-earned knowledge that the straight skeleton structure was inherently computationally difficult to generate: it needs exact arithmetic for robustness, and thus the CGAL dependency is warranted. 🤷

Creating a CGAL-based straight skeleton package

Adopting CGAL wasn’t an insubstantial amount of work: significant amounts of infrastructure was needed to import R spatial/geometry data, calculate the skeleton, export out what I needed to R, and then actually do something interesting with the skeleton. CGAL outputs the straight skeleton structure as a half-edge data structure, which we will reduce to a simple directed graph. We’re then left with a simple directed graph with heights for each node. Additionally, the straight skeleton was only half the battle: I then needed to transform that data into 2D polygons and 3D meshes suitable for rendering in R. And not just simple 3D rooftops and simple 3D bevels: I wanted to support complex curved bevels, and thus needed a custom mesh generation algorithm.

Note

CGAL actually added the ability to generate meshes and simple 3D bevels halfway through 2023, about 4 months after I started working on this project: strange how these gaps in capability often seem to be noticed and solved in parallel! Thankfully, my implementation was more ambitious so my work was not wasted, but if that’s all I was after I could have just stopped here and written a lightweight wrapper around CGAL!

3D Meshing Algorithms

Roof generation was the easy part: let’s say we have a straight skeleton of a polygon and we want to generate a roof model. The roof height is already present by virtue of the straight skeleton providing a distance to each node from the edge: all we have to do is turn all the closed loops in the straight skeleton into polygons, and use the distance as the height of each vertex. This will give us a sloped roof (which you can scale to get different angles). The algorithm for this is relatively simple:

Call an internal {raybevel} function that converts a straight skeleton structure to polygons indices, pointing to vertices in the skeleton node structure

2

Generate the straight skeleton ggplot2 layers

3

Generate the straight skeleton polygon layers (one for each polygon)

4

Loop generating frames of the animation.

3D Roof Meshes

Select a link not directly on the edge. Mark it as visited. Loop around the skeleton, taking the left-most turn at each node and recording the indices of the nodes you visit. When you return to the original vertex, those vertices form an interior polygon. Repeat for all non-edge links. When you’ve exhausted all links, de-duplicate the polygons by (in a copy of the list of indices) sorting each polygon’s indices and hashing the resulting vector. Earcut the polygons to triangulate the mesh. See an animated version of this algorithm above in Figure 11.

The straight skeleton structure defines the height of each vertex as the distance from the nearest edge, which means this algorithm works great for constant-slope rooftop generation. Let’s apply this algorithm to a few thousand houses and see how it looks in rayshader, using the new render_buildings() feature:

Specify the bounding box for the region we are pulling data from.

3

Get building polygon data road line data from OpenStreetMap

4

Get the actual extent the loaded building polygon data and use that to load elevation data using the elevatr package.

5

Crop DEM, but note that the cropped DEM will have an extent slightly different than what’s specified in e. Save that new extent to new_e.

6

Crop the DEM.

7

Visualize areas less than one meter as water (approximate tidal range)

8

Render the 3D rayshader scene

9

Generate and render the 3D building meshes from the polygon data

10

Render the scene with rayrender’s interactive pathtracer

3D Arbitrary Bevel Meshes

However, if we want to generate an arbitrary bevel instead of a simple roof, it gets a lot more complex: we need to insert new nodes at each distance defined in the bevel and link those horizontally nodes so they form constant height contours around the polygon. We then split the existing links at those bevel points and traverse the straight skeleton to link the new nodes together to form our contours. This latter step is trickier: as the interior polygon “shrinks” at higher interior distances, regions of the polygon can become disconnected from each other.

Get the maximum internal distance out of the straight skeleton structure

4

Create offsets for these, except for the first and last values

5

Loop generating different offsets, showing the polygon breaking into several

6

Animate using {av}

This complicates the meshing process significantly. There can be multiple polygons for a single interior polygon, and you need to be careful when traveling through the straight skeleton structure to only access areas that aren’t cut off yet. We ensure we do by scanning the entire straight skeleton network for local maxima: the peaks on the roof correspond to interior polygons that were topologically cut-off at one point.

We then always start our search from these maxima (when our contour is below that height) and only search the areas of the skeleton reachable from that point by always turning around when we insert a new node. This ensures we won’t connect nodes to other nodes that aren’t actually in the same interior polygon.

After modifying our straight skeleton structure with new links, we just apply the same polygonization algorithm as in the simple rooftop case! To apply our custom bevel, we just take each node’s distance from the edge and interpolate its height to the height specified in the input bevel. This transforms the interior of our polygon into any 3D profile we want! Here’s some that are build into {raybevel}:

raybevel also includes helper functions to help you create complex bevels, although you can also just manually pass a list with x/y coordinates with your own bevel.

Generate the bevel manually using mathematical functions

2

Render 3D beveled California mesh using rayrender

Or crazy things like make SVG fonts into double-sided bubble letters by applying a circular edge contour:

I’ve packaged up all these functions into the new package raybevel: simply input polygons or sf objects and transform them to 3D, or create 2D inset polygons! I’ve used it throughout this post, I’ve also included interfaces built-in to rayshader (render_buildings()) and render_beveled_polygons()).

Creating the {raybevel} hexlogo

Finally, it’s my personal philosophy that any R package focused on making cool and interesting visualizations should use itself to generate it’s package hexlogo. So I pulled up Inkscape and started working on a logo…

I decided I wanted to have bolts added to the edge plates on the center lines of the polygons, but there was no way in Inkscape to do that easily. If only there was some way to get the mid-line of the polygon! Oh wait, didn’t I just write a package for generating straight skeletons? One round trip to R/raybevel and back to SVG later…

With bevels applied to each polygon along with a metallic material and knurling added as a bump map, we then end up with this sparkling new pathtraced hexlogo:

Awesome! The package is complete. Now R users can make 2D inset polygons, beautiful beveled 3D shapes, and slick 3D rooftops. Happy New Year!