Use RTINs

[kylebarron]

The idea is

terrain rgb -> martini -> viewshed algorithm based on RTINs

This would probably be really fast, and it could allow for the possibility of more accuracy closer to the viewer and less accuracy further away.

"Algorithms for visibility computation on terrains: a survey" has a great general overview of different viewshed modeling techniques and algorithms. Of special interest is chapter 6, "Visibility on multiresolution terrain models". Martini lets me create an RTIN, a type of multiresolution terrain model, really fast. Section 6.1 talks about ways to create these multiresolution models. Section 6.2 discusses ways to run visibility algorithms on these models:

Algorithms have been presented that compute the visibility map on a hierarchical TIN by navigating the tree-like structure of the model and thus do not need the explicit construction of a TIN at the given resolution (De Floriani and Magillo, 1997). Triangles forming the DEM are extracted one by one from the multiresolution model, in front-to-back order. The key point is that resolution decreases when moving from front to back.

This guarantees that triangles that are extracted first may result in the refinement of triangles that are extracted later, but the opposite cannot happen. Thus, back-tracking to already-processed triangles is never needed. The method was originally designed for nested hierarchies of TINs, but it can be applied, with no need for change, to hierarchies of right triangles, which recently have become popular terrain representations for gridded data.

De Floriani and Magillo, 1997 is here.

[kylebarron]

Martini example:

var require = require("esm")(module/*, options*/)
Martini = require('@mapbox/martini').default
var getPixels = require("get-pixels")

// Tuolumne meadows
x = 690
y = 1581
z = 12

var path = `/Users/kyle/github/mapping/nst-guide/hillshade/data/terrain_png/${z}/${x}/${y}.png`
let data
getPixels(path, function(err, pixels) {
  if(err) {
    console.log("Bad image path")
    return
  }
  data = pixels
  console.log("got pixels", pixels.shape.slice())
})

tileSize = data.shape[0]
gridSize = tileSize + 1

var terrain = new Float32Array(gridSize * gridSize);

let R, G, B;
for (let x = 0; x < tileSize; x++) {
  for (let y = 0; y < tileSize; y++) {
    R = data.get(x, y, 0)
    G = data.get(x, y, 1)
    B = data.get(x, y, 2)
    height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
    terrain[y * gridSize + x] = height;
  }
}

// backfill right and bottom borders
for (let x = 0; x < gridSize - 1; x++) {
  terrain[gridSize * (gridSize - 1) + x] = terrain[gridSize * (gridSize - 2) + x];
}
for (let y = 0; y < gridSize; y++) {
  terrain[gridSize * y + gridSize - 1] = terrain[gridSize * y + gridSize - 2];
}

// set up mesh generator for a certain 2^k+1 grid size
var martini = new Martini(gridSize);

// generate RTIN hierarchy from terrain data (an array of size^2 length)
var tile = martini.createTile(terrain);

// get a mesh (vertices and triangles indices) for a 10m error
var mesh = tile.getMesh(10);

Now mesh has two arrays:

{
  vertices: Uint16Array [
    288,  96, 320, 128, 320,  96, 320,  64, 304,  80, 288, 128,
    272, 112, 256, 128, 336, 112, 352, 128, 352, 112, 352,  96,
    344, 104, 368, 112, 368, 128, 384, 128, 376, 120, 360, 104,
    336,  80, 336,  96, 344,  88, 280,  56, 288,  48, 280,  48,
    276,  52, 276,  48, 272,  48, 274,  50, 288,  64, 288,  32,
    280,  40, 276,  44, 268,  60, 272,  56, 268,  56, 264,  56,
    272,  64, 272,  60, 268,  52, 272,  52, 256,  64, 264,  64,
    280,  64, 276,  60, 276,  64, 274,  62, 274,  64, 276,  56,
    274,  54, 274,  56,
    ... 4008 more items
  ],
  triangles: Uint16Array [
     0,  1,  2,  2,  3,  4,  0,  2,  4,  5,  0,  6,
     7,  5,  6,  1,  0,  5,  8,  9, 10, 10, 11, 12,
     8, 10, 12,  1,  9,  8, 13,  9, 14, 14, 15, 16,
    13, 14, 16, 10, 13, 17, 11, 10, 17,  9, 13, 10,
    18,  2, 19, 19, 11, 20, 18, 19, 20,  3,  2, 18,
     2,  1,  8, 19,  8, 12, 11, 19, 12,  2,  8, 19,
    21, 22, 23, 24, 23, 25, 25, 26, 27, 24, 25, 27,
    21, 23, 24, 28, 22, 21, 22, 29, 30, 23, 30, 31,
    31, 26, 25, 23,
    ... 11750 more items
  ]
}

From inspection, it looks like vertices is a flattened list of x, y coordinates. So that the vertices are (288, 96), (320, 128) and so on. Then triangles references the order of these coordinates, in its own flattened list of triples. So 0, 1, 2 is a triangle, which corresponds to the first three coordinates of vertices, or (288, 96), (320, 128), (320, 96). You can tell because these are right triangles. So two of the x or y coordinates must be the same, here 320 and 96.

[kylebarron]

So process, at least for a single sightline could be:

• load png data (what's the fastest way to get the raw image data if it's already been downloaded for hillshading?)
• Stitch multiple hillshade files together. Note that martini needs a square array... So either run martini on each tile and then stitch together, or stitch together tiles ahead of time, even if it's more than what's in the viewport, if the viewport isn't square.
• Calculate terrain height values array
• Pass to Martini
• Make kdbush index of triangle vertices (or rbush index of triangle edges?)
• Find coordinate of viewer in the image data's coordinate system
• Search index

For more than just a single sightline, also:

• Convert pixel coordinates to lat/lon
• Create GeoJSON

[kylebarron]

Note that since you're not creating a 3d terrain, running this per tile might be the best. The issue with the tile-based method, instead of stitching all the elevation values together from the outset, is that you could get artifacts where the tiles meet, where the mesh doesn't line up exactly. But this would only matter in my case if right there at the tile boundary there should exist but is missing an elevation value that blocks visibility. But it's likely that neighboring triangles would block anyways. See Martini #7 for an example of the problem with 3D imagery.

[kylebarron]

https://www.researchgate.net/profile/Leila_De_Floriani2/publication/220387766_Visibility_Computations_on_Hierarchical_Triangulated_Terrain_Models/links/55c3b29608aeb97567401b18/Visibility-Computations-on-Hierarchical-Triangulated-Terrain-Models.pdf https://www.cs.ubc.ca/~will/papers/rtin.pdf