I arrived here by Googling this problem myself, read the linked articles, and produced a relatively compact solution which generates the common set of 47 tiles. It requires a 2x3 tileset for the autotiled material like this: 
With a single-tile variant at the top left, inner corners at the top right, and four outer corner tiles at the bottom (you may recognize this arrangement from RPG Maker).
The trick is to break each "logical" map tile into 4 half-tiles for rendering. furthermore, a half-tile in the tileset can only be in that position in a generated tile, so a top-left half tile can only be used in a top-left position.
These restrictions mean that you only need to check 3 full-tile neighbors per half-tile, instead of all 8 neighboring tiles.
I implemented this idea quickly to test it. Here's the proof-of-concept code (TypeScript):
//const dirs = { N: 1, E: 2, S: 4, W:8, NE: 16, SE: 32, SW: 64, NW: 128 };
const edges = { A: 1+8+128, B: 1+2+16, C: 4+8+64, D: 4+2+32 };
const mapA = { 0:8, 128:8, 1:16, 8:10, 9:2, 137:18, 136:10, 129:16 };
const mapB = { 0:11, 16:11, 1:19, 2:9, 3:3, 19:17, 18:9, 17:19 };
const mapC = { 0:20, 64:20, 4:12, 8:22, 12:6, 76:14, 72:22, 68:12 };
const mapD = { 0:23, 32:23, 4:15, 2:21, 6:7, 38:13, 34:21, 36:15 };
export function GenerateAutotileMap(_map: number[][], _tile: integer): number[][]
{
var result = [];
for (var y=0; y < _map.length; y++) {
const row = _map[y];
const Y = y*2;
// half-tiles
result[Y] = [];
result[Y+1] = [];
// each row
for (var x=0; x < row.length; x++) {
// get the tile
const t = row[x];
const X = x*2;
if (t != _tile) continue;
// Check nearby tile materials.
const neighbors = (North(_map, x, y) == t? 1:0)
+ (East(_map, x, y) == t? 2:0)
+ (South(_map, x, y) == t? 4:0)
+ (West(_map, x, y) == t? 8:0)
+ (NorthEast(_map, x, y) == t? 16:0)
+ (SouthEast(_map, x, y) == t? 32:0)
+ (SouthWest(_map, x, y) == t? 64:0)
+ (NorthWest(_map, x, y) == t? 128:0);
// Isolated tile
if (neighbors == 0) {
result[Y][X] = 0;
result[Y][X+1] = 1;
result[Y+1][X] = 4;
result[Y+1][X+1] = 5;
continue;
}
// Find half-tiles.
result[Y][X] = mapA[neighbors & edges.A];
result[Y][X+1] = mapB[neighbors & edges.B];
result[Y+1][X] = mapC[neighbors & edges.C];
result[Y+1][X+1] = mapD[neighbors & edges.D];
}
}
return result;
}
Explanation:
A is the top left part of the tile, B is the top right, C is the bottom left, D is the bottom right.edges holds bitmasks for each of these, so we can grab only the relevant neighbor info.map* are dictionaries mapping neighbor states to graphic indices in the tileset image (0..24).
- since each half-tile checks 3 neighbors, each one has 2^3=8 states.
_tile is the tile targeted for autotiling. - Since our logical tiles are twice as large as our rendering tiles, all the autotile coords (x,y) have to be doubled in the rendering map.
Anyway, here are the results (with only one tile, anyway):
