I have the rendering pipeline mostly working, but I need something to actually render.  At some point world gen will be required, so I thought I should just implement it now.

Most world gen uses perlin noise to create the heightmap, which is a critical feature of the terrain - perhaps the most important feature.  The problem, is that it isn't very realistic on a larger scale.  Here is a map of a rather large Minecraft world, generated with perlin noise like this, and it does not have realistic continents, islands, biomes, etc:
In real life, there are island chains, separated continents, deep trenches under the oceans, mountain chains, etc.  All of these features are created by plate tectonics.  But how can be simulate plate tectonics in a procedural generation algorithm?

Plate Tectonics in an Infinite Procedurally Generated World
When I started working on this world generator, I started by thinking that maybe I could generate a list of nearby points to a chunk, and make faultlines as lines between them, maybe modified with some noise.  The issue, is that with a procedurally generated and infinite world, keeping track of and searching these points is not really feasible, and as far as I know, there is not an algorithm to compute the nearby points based on arbitrary given coordinates.  This is important, because when we generate the world, we generate it one chunk at a time, where a chunk is a 16x16x64 tower of blocks that forms one section of the world.  We cannot generate the entire world at once, because the world may be too large, even infinite.  In essence, it isn't feasible to generate tectonic plates themselves.  Thankfully, there is a great workaround.

A procedural world generator like this is, in essence, a function that describes the characteristics of the world at point x, y given a seed value.  So, we need a function that describes the tectonic plate and/or nearby plates at the coordinate x, y.  And we must do this without having a list of tectonic plates to pull from, as the number of tectonic plates is infinite.

Bent Space
Instead of generating the tectonic plates, it works far better to start with a square grid of tectonic plates, then project the x and y coordinates onto this grid using a bunch of cool math.  Here is the same world as the one I showed in the first image, but without coordinate transformation.  Note that faultines are still simulated:As you can see, it is very blocky.  All I did to get from this to the world you saw in the first image, was apply changes to the x and y coordinates.  This makes the blocky grid disappear, and allows us to bypass that tectonic plate issue.

Steps of Coordinate Alteration
The sequence of steps I apply is as follows:

1. Apply a sine function to the coordinates, with a changing frequency
2. Perlin noise is added to the coordinates.  This seems to produce more realistic terrain than adding perlin noise to the heightmap.
3. Some perlin noise is added again, but with different parameters and math
4. Math is done to distort the corners of the grid.  This causes the corners to smooth and disappear into the terrain, as well as creates more interesting terrain features.
5. Sometimes the corner of a square of the grid is removed, and sometimes the entire square is divided and given to the nearby squares.  This is chosen and controlled by value noise, and makes the features even less blocky.  At this point the grid pattern we started with is nearly imperceptible.

For each of these features, a world config controls their strength, and the strength is then modulated with more perlin noise.  Each layer of perlin noise should have a different seed, although looking at my code, it seems I didn't do this.  Anyway, different regions of the world should have different coordinate-altering characteristics, leading to differences in the way the terrain looks and feels in different regions.

This coordinate remapping process produces a lot more interesting terrain in my opinion.

Generating Plate Characteristics

This is one of the simpler steps.  Once we have our coordinates, I run a function like this:

# vNoise means 'vector value noise' - produces a random 2d vector given input coordinates.
def vNoise(x: int, y: int, seed: int):
    # 2d value noise function I call 'gridNoise', as it makes a grid.
    # This function converts the converted coordinates into a hash for each plate.
    h = gridNoise(int(x), int(y), seed)
    r = random.Random(h)
    return r.random() * 2 - 1, r.random() * 2 - 1

 This generates a 2d vector, which represents the direction of movement of the plate.  Plates move in different directions and with different speeds.  By comparing the speeds of nearby plates, we can determine the type of faultline between the plates, and generate a realistic heightmap accordingly.

Since different plates also have different overall heights, a modified perlin noise function is run to generate height data for plates themselves:

# sx and sy are components of the transformed coordinates
plateHeight = noise.pnoise2((int(sx) + 0.3 * sx) / 2, (int(sy) * 0.3) / 2, 3)*2 - 0.5

There is some influence for the side of the plate you are on, but the plate itself dominates in creating the height.  This perlin noise function here is largely responsible for the differences in plate height seen in the world gen images.

So, we essentially have two characteristics for each plate:  height, and velocity.

Faultlines and Heightmap Generation

Each faultline has two characteristics, which I call 'pressure' and 'shift'.  Real geologists may have their own terms, but those are mine.  Pressure is high when plates are colliding, and it is low when they are dividing.  Shift describes how much plates are moving beside eachother.

The change to the hightmap is computed my lerping between the plate height and the fault pressure, according to how close the current x, y coordinate is to a fautline.  Perhaps it is better described in code:

# Compute the location of the nearest plate, so we can compare the velocity
edgeX = abs(sx - round(sx))
edgeY = abs(sy - round(sy))
xCloser = edgeX < edgeY
nearX = squareNear(sx) if xCloser else sx
nearY = squareNear(sy) if not xCloser else sy

# Get the velocity of this plate and the nearest plate, then compute the plate height
dx, dy = vNoise(sx, sy, seed)
dxNear, dyNear = vNoise(nearX, nearY, seed)
plateHeight = noise.pnoise2((int(sx) + 0.3 * sx) / 2, (int(sy) * 0.3) / 2, 3)*2 - 0.5

# Compute the pressure and shift
pressure = (dx - dxNear) / (nearX - sx) if xCloser else (dy - dyNear) / (nearY - sy)
shift = (dy - dyNear) if xCloser else (dx - dxNear)

# Determine how much the current x and y are affected by this faultline
q = 0.003 / (0.003 + edgeX * edgeY)

# This is done to prevent long, thin islands in the ocean due to weak faultlines.
# Mountainous island chains are still generated.
largeFault = abs(pressure) > 1.2

# Here is how the heightmap is finally computed
heightmap[x][y] = (plateHeight * (1 - q) + q * pressure) if largeFault else plateHeight

I am not currently using shift, but it could be cool to help generate other terrain features.  Ideally the terrain on each side of a shifting fault would be offset in opposite directions, making it look like the terrain slid in the world's past.

To compare with the top image, here is what the heightmap looks like without computing for faultlines:The terrain still looks good, but there are no island chains, deep trenches, or fault-related mountain chains.

Geology Computations
The fault data can be used to generate volcanos, hydrothermal vents, mineral and rock information, etc.  This is useful for generating resources, but I haven't really implemented all that.  So far I only do volcano and vent generation here, for a little volcanic activity.  Later it would be cool to use the geology data to generate minerals and resources.

Smoothing
The methods above can produce a noisy and jagged map.  For this reason, I use a 5x5 average to smooth the terrain.  Here's what it looks like without smoothing:The smoothing may be too strong.  I don't know.  But it certainly needs some smoothing.

Areas for other work and improvement

Some parts of the map still have evidence of the underlying grid pattern for the tectonic plates.  You can see this in the upper, leftish-middle portion of the map, where there is this odd plus-shaped island.  This is an artifact from the grid used to generate the map.  It's harder to see, but there is another artifact below that.  It is particularly visible in oceans.  Here, I circled them and drew lines to make it more obvious:The orange lines show the top-to-bottom dividing lines on the grid.  The yellow lines show side-to-side.  Red circles regions that look unnatural due to the grid showing through.  Now that I point out the grid, can you see the pattern in the volcanic activity, coastlines, and islands?  Suddenly it looks a lot more organized.  Ideally, this would be completely indiscernible, and for most of the map, it is.

A compete world generator should also have mechanisms for generating the smaller, finer features of the world.  I have not implemented such mechanisms really at all.  This includes rivers, lakes, hills, biomes, etc.

Another thing I'd like to add to this system is a sensible and realistic climate, for generating biomes and snow.

Real quick, here's an image of the heightmap this world has:Some of the trenches in the north are very deep, and some of the faultline islands have tall mountains.  Very fun.

Porting to Arduino
After I got an algorithm I like, I needed to port it to Arduino.  This was much easier, and the worldgen runs quickly too:You can see that the world has generated, although my ability to look at it is quite limited at the moment.  Each chunk generates in about 16 to 17 milliseconds.  Since it is not guaranteed that any particular chunk will have land, or be without volcanoes, I had to make a function to generate a few chunks until it finds one with a suitable spawnpoint.

Rendering
Unfortunately, this world is not currently rendering, so now I have to figure that out.  I need to find a good rendering algorithm to do it fast enough, but I have some ideas that I think will help.  If I make the rendering function do a floodfill from the center of the camera, it can render only the blocks the player sees, and it should be much faster than iterating over every block in a chunk.  I can also stop the floodfill after a certain number of blocks render or after a certain amount of time, and since the most nearby blocks are rendered first, the game should still look good.  62s to render one chunk is just way too long, but I know I can get it down to something reasonable.

Other Work
The last thing I want to say - for some reason I never thought to google worldgen based on plate tectonics before writing this.  So I did google it, and found some interesting stuff:

https://www.gamedev.net/forums/topic/623145-terrain-generation-with-plate-tectonics/
OP here says "It seems to me like even the most naïve model of plate tectonics is able to produce more convincing heightmaps than conventional fractal based methods."  This seems accurate to me.  The reply from another user: "For extremely large scale terrain you simply can't have all of the terrain in memory. This is the main reason why fractal-based methods will continue to be very popular."  This is exactly the issue I addressed, and his point makes perfect sense.  The code OP wrote actually simulates the tectonic plates, so it isn't really a procedural generator.  Nonetheless, the terrain generated looks good at least at a quick glance, and may be worth looking into.  He has a thesis here: https://www.theseus.fi/handle/10024/40422

The other one I found is here: http://procgenesis.com/WorldGen/worldgen.html  It doesn't seem to be open source, but the terrain looks good and there's other cool generators on that site that may be worth looking at, including a language generator.