Finer earth tiles for when you need really crisp pictures
This commit is contained in:
connor committed
1 parent
5702671960
commit
7d91c87fe0
21 files changed
+2644
-122
No files matched your search
Generated
+18
@@ -1478,6 +1478,12 @@ version = "0.1.30"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "d55d956fa96f5ec02be2e13af0e20391a5aa83d6a074e3ad368959d0fab299ea"
|
||||
|
||||
[[package]]
|
||||
name = "quick-error"
|
||||
version = "2.0.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "a993555f31e5a609f617c12db6250dedcac1b0a85076912c436e6fc9b2c8e6a3"
|
||||
|
||||
[[package]]
|
||||
name = "quote"
|
||||
version = "1.0.47"
|
||||
@@ -1615,6 +1621,7 @@ dependencies = [
|
||||
"intel_tex_2",
|
||||
"resolved",
|
||||
"serde",
|
||||
"tiff",
|
||||
"toml",
|
||||
]
|
||||
|
||||
@@ -1907,6 +1914,17 @@ dependencies = [
|
||||
"syn 3.0.9",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "tiff"
|
||||
version = "0.11.3"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "b63feaf3343d35b6ca4d50483f94843803b0f51634937cc2ec519fc32232bc52"
|
||||
dependencies = [
|
||||
"flate2",
|
||||
"half",
|
||||
"quick-error",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "tinytemplate"
|
||||
version = "1.2.1"
|
||||
|
||||
@@ -31,6 +31,8 @@ intel_tex_2 = "0.5"
|
||||
half = { version = "2.7", features = ["bytemuck"] }
|
||||
bytemuck = "1.0"
|
||||
hifitime = { version = "4.3", optional = true }
|
||||
# ESA WorldCover's GeoTIFFs for `earth-regions` (deflate-compressed blocks).
|
||||
tiff = { version = "0.11", default-features = false, features = ["deflate"] }
|
||||
|
||||
[dev-dependencies]
|
||||
# JPEG too, so the tests can write a small stand-in for the Blue Marble images.
|
||||
|
||||
@@ -0,0 +1,743 @@
|
||||
//! `resolved earth-regions`: 1 x 1 degree cells of ESA WorldCover's Sentinel-2 composite
|
||||
//! (2021, version 200: the median of the year's cloud-free L2A surface reflectance at
|
||||
//! 10 m), as ESA publishes them (12000 x 12000 GeoTIFFs of red, green, blue and
|
||||
//! near-infrared, 1e-4 reflectance per unit, 0 where there is no data), each made into
|
||||
//! one tile of [`DetailLayout::WORLDCOVER`] and compressed to BC6H, once, ahead of
|
||||
//! rendering.
|
||||
//!
|
||||
//! A cell's border holds its neighbors' texels where their sources are on hand when it is
|
||||
//! built, and repeats its own edge elsewhere; its header says which. A cell is rebuilt
|
||||
//! when the neighbors on hand change, so a region grown later joins up seamlessly.
|
||||
|
||||
use std::io::Write;
|
||||
use std::path::{Path, PathBuf};
|
||||
use std::time::Instant;
|
||||
|
||||
use resolved::earth::{
|
||||
Cell, DetailFormat, DetailLayout, EarthDetail, REGION_COARSE, REGION_COARSE_PADDED,
|
||||
};
|
||||
use tiff::decoder::{Decoder, DecodingResult, Limits};
|
||||
use tiff::tags::Tag;
|
||||
|
||||
use crate::earth_tiles::compress;
|
||||
|
||||
/// Texels across a cell's interior.
|
||||
const SIZE: usize = 12000;
|
||||
|
||||
/// Reflectance per unit of the source's values.
|
||||
const SCALE: f32 = 1e-4;
|
||||
|
||||
/// The source file for `cell`, as ESA names it.
|
||||
pub fn source_name(cell: Cell) -> String {
|
||||
format!("ESA_WorldCover_10m_2021_v200_{cell}_S2RGBNIR.tif")
|
||||
}
|
||||
|
||||
/// The tiled file for `cell`.
|
||||
pub fn tiles_name(cell: Cell) -> String {
|
||||
format!("{cell}.tiles")
|
||||
}
|
||||
|
||||
fn open(path: &Path) -> Result<Decoder<std::io::BufReader<std::fs::File>>, String> {
|
||||
let fail = |e: &dyn std::fmt::Display| format!("{}: {e}", path.display());
|
||||
let file = std::fs::File::open(path).map_err(|e| fail(&e))?;
|
||||
let mut decoder = Decoder::new(std::io::BufReader::new(file))
|
||||
.map_err(|e| fail(&e))?
|
||||
.with_limits(Limits::unlimited());
|
||||
if decoder.dimensions().map_err(|e| fail(&e))? != (SIZE as u32, SIZE as u32) {
|
||||
return Err(fail(&format!("not {SIZE} x {SIZE}")));
|
||||
}
|
||||
// Four 16-bit bands: red, green, blue, near-infrared.
|
||||
let kind = decoder.colortype().map_err(|e| fail(&e))?;
|
||||
if !matches!(
|
||||
kind,
|
||||
tiff::ColorType::Multiband {
|
||||
bit_depth: 16,
|
||||
num_samples: 4
|
||||
} | tiff::ColorType::RGBA(16)
|
||||
| tiff::ColorType::CMYK(16)
|
||||
) {
|
||||
return Err(fail(&format!("{kind:?}, not four 16-bit bands")));
|
||||
}
|
||||
Ok(decoder)
|
||||
}
|
||||
|
||||
/// Checks that the file at `path` covers `cell`: its north-west corner at the cell's,
|
||||
/// and 1 / 12000 degree a texel.
|
||||
fn check_place(
|
||||
decoder: &mut Decoder<std::io::BufReader<std::fs::File>>,
|
||||
path: &Path,
|
||||
cell: Cell,
|
||||
) -> Result<(), String> {
|
||||
let fail = |e: &dyn std::fmt::Display| format!("{}: {e}", path.display());
|
||||
let scale = decoder
|
||||
.get_tag(Tag::Unknown(33550))
|
||||
.and_then(|v| v.into_f64_vec())
|
||||
.map_err(|e| fail(&format!("pixel scale: {e}")))?;
|
||||
let tie = decoder
|
||||
.get_tag(Tag::Unknown(33922))
|
||||
.and_then(|v| v.into_f64_vec())
|
||||
.map_err(|e| fail(&format!("tie point: {e}")))?;
|
||||
let texel = 1.0 / SIZE as f64;
|
||||
let ok = scale.len() >= 2
|
||||
&& tie.len() >= 6
|
||||
&& (scale[0] - texel).abs() < 1e-9
|
||||
&& (scale[1] - texel).abs() < 1e-9
|
||||
&& tie[..2] == [0.0, 0.0]
|
||||
&& (tie[3] - f64::from(cell.west)).abs() < 1e-9
|
||||
&& (tie[4] - f64::from(cell.south + 1)).abs() < 1e-9;
|
||||
if ok {
|
||||
Ok(())
|
||||
} else {
|
||||
Err(fail(&format!(
|
||||
"placed at {tie:?} with texels {scale:?}, not {cell}'s"
|
||||
)))
|
||||
}
|
||||
}
|
||||
|
||||
/// The red, green and blue of every texel in rows from the north-west, as the source's
|
||||
/// values, whether each holds data, and whether it is water: green brighter than near
|
||||
/// infrared (McFeeters' water index above 0).
|
||||
struct Bands {
|
||||
size: usize,
|
||||
rgb: Vec<[u16; 3]>,
|
||||
valid: Vec<bool>,
|
||||
water: Vec<bool>,
|
||||
}
|
||||
|
||||
fn decode(path: &Path, cell: Cell) -> Result<Bands, String> {
|
||||
let fail = |e: &dyn std::fmt::Display| format!("{}: {e}", path.display());
|
||||
let mut decoder = open(path)?;
|
||||
check_place(&mut decoder, path, cell)?;
|
||||
let values = match decoder.read_image().map_err(|e| fail(&e))? {
|
||||
DecodingResult::U16(v) => v,
|
||||
_ => return Err(fail(&"not 16-bit")),
|
||||
};
|
||||
if values.len() != 4 * SIZE * SIZE {
|
||||
return Err(fail(&format!("{} values, not 4 per texel", values.len())));
|
||||
}
|
||||
let rgb: Vec<[u16; 3]> = values.chunks_exact(4).map(|p| [p[0], p[1], p[2]]).collect();
|
||||
let valid = rgb.iter().map(|p| p.iter().all(|&v| v > 0)).collect();
|
||||
let water = values.chunks_exact(4).map(|p| p[1] > p[3]).collect();
|
||||
Ok(Bands {
|
||||
size: SIZE,
|
||||
rgb,
|
||||
valid,
|
||||
water,
|
||||
})
|
||||
}
|
||||
|
||||
/// The texels of the strip of the source at `path` (`cell`) that a neighbor's border
|
||||
/// needs: columns `x0..x1` and rows `y0..y1`, read only from the 1024 x 1024 blocks it
|
||||
/// crosses, as reflectance in rows; `None` for each texel without data.
|
||||
fn read_strip(
|
||||
path: &Path,
|
||||
cell: Cell,
|
||||
[x0, x1, y0, y1]: [usize; 4],
|
||||
) -> Result<Vec<Option<[f32; 3]>>, String> {
|
||||
let fail = |e: &dyn std::fmt::Display| format!("{}: {e}", path.display());
|
||||
let mut decoder = open(path)?;
|
||||
check_place(&mut decoder, path, cell)?;
|
||||
let (cw, ch) = decoder.chunk_dimensions();
|
||||
let (cw, ch) = (cw as usize, ch as usize);
|
||||
let across = SIZE.div_ceil(cw);
|
||||
let (w, h) = (x1 - x0, y1 - y0);
|
||||
let mut out = vec![None; w * h];
|
||||
for by in y0 / ch..=(y1 - 1) / ch {
|
||||
for bx in x0 / cw..=(x1 - 1) / cw {
|
||||
let index = (by * across + bx) as u32;
|
||||
let (dw, _) = decoder.chunk_data_dimensions(index);
|
||||
let values = match decoder.read_chunk(index).map_err(|e| fail(&e))? {
|
||||
DecodingResult::U16(v) => v,
|
||||
_ => return Err(fail(&"not 16-bit")),
|
||||
};
|
||||
for y in (by * ch).max(y0)..((by + 1) * ch).min(y1) {
|
||||
for x in (bx * cw).max(x0)..((bx + 1) * cw).min(x1) {
|
||||
let i = 4 * ((y - by * ch) * dw as usize + (x - bx * cw));
|
||||
let p = [values[i], values[i + 1], values[i + 2]];
|
||||
if p.iter().all(|&v| v > 0) {
|
||||
out[(y - y0) * w + (x - x0)] = Some(p.map(|v| f32::from(v) * SCALE));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
Ok(out)
|
||||
}
|
||||
|
||||
/// Texels across a cell at the level where holes are filled across a group of cells
|
||||
/// ([`GroupFill`]): about 300 m.
|
||||
const JOINT: usize = 375;
|
||||
|
||||
/// Per texel of a `size` x `size` grid, the sums of red, green and blue over what it
|
||||
/// covers and how much of it counts (0 to 1), filled where nothing counts from what
|
||||
/// does around it (pull-push): each level of a pyramid holds the sums of 2 x 2 texels
|
||||
/// of the one above, and from the coarsest down, each level's means are blended with
|
||||
/// the coarser level's filled values, read bilinearly, by the share of it that counts.
|
||||
/// With `seed`, a level of that size takes its filled values instead of those of the
|
||||
/// levels above it. Returns the mean of every texel, filled.
|
||||
fn pull_push(
|
||||
size: usize,
|
||||
first: Vec<[f32; 4]>,
|
||||
seed: Option<(usize, &[[f32; 3]])>,
|
||||
) -> Vec<[f32; 3]> {
|
||||
let mut sums = vec![(size, first)];
|
||||
while sums.last().expect("a level").0 > 1 {
|
||||
let (n, above) = sums.last().expect("a level");
|
||||
let (n, m) = (*n, n.div_ceil(2));
|
||||
let mut below = vec![[0.0f32; 4]; m * m];
|
||||
resolved::earth::parallel_rows(&mut below, m, |y, row| {
|
||||
for (x, t) in row.iter_mut().enumerate() {
|
||||
for (dx, dy) in [(0, 0), (1, 0), (0, 1), (1, 1)] {
|
||||
let (sx, sy) = ((2 * x + dx).min(n - 1), (2 * y + dy).min(n - 1));
|
||||
let p = above[sy * n + sx];
|
||||
for c in 0..4 {
|
||||
t[c] += p[c];
|
||||
}
|
||||
}
|
||||
}
|
||||
});
|
||||
sums.push((m, below));
|
||||
}
|
||||
// The level the push starts from, filled: the seed's, or the coarsest's mean.
|
||||
let start =
|
||||
seed.and_then(|(m, values)| sums.iter().position(|(n, _)| *n == m).map(|k| (k, values)));
|
||||
let (first_k, mut filled, mut width) = match start {
|
||||
Some((k, values)) => (k, values.to_vec(), sums[k].0),
|
||||
None => {
|
||||
let top = sums.last().expect("a level").1[0];
|
||||
let mean = if top[3] > 0.0 {
|
||||
[top[0] / top[3], top[1] / top[3], top[2] / top[3]]
|
||||
} else {
|
||||
[0.0; 3]
|
||||
};
|
||||
(sums.len() - 1, vec![mean], 1)
|
||||
}
|
||||
};
|
||||
for k in (0..first_k).rev() {
|
||||
let (n, level) = &sums[k];
|
||||
let n = *n;
|
||||
let coarse = &filled;
|
||||
let m = width;
|
||||
let full = 4f32.powi(k as i32);
|
||||
let at = |x: f32, y: f32| -> [f32; 3] {
|
||||
let (x, y) = (x.clamp(0.0, (m - 1) as f32), y.clamp(0.0, (m - 1) as f32));
|
||||
let (x0, y0) = (x.floor() as usize, y.floor() as usize);
|
||||
let (x1, y1) = ((x0 + 1).min(m - 1), (y0 + 1).min(m - 1));
|
||||
let (fx, fy) = (x - x0 as f32, y - y0 as f32);
|
||||
[0, 1, 2].map(|c| {
|
||||
let v = |x: usize, y: usize| coarse[y * m + x][c];
|
||||
(v(x0, y0) * (1.0 - fx) + v(x1, y0) * fx) * (1.0 - fy)
|
||||
+ (v(x0, y1) * (1.0 - fx) + v(x1, y1) * fx) * fy
|
||||
})
|
||||
};
|
||||
let mut next = vec![[0.0f32; 3]; n * n];
|
||||
resolved::earth::parallel_rows(&mut next, n, |y, row| {
|
||||
for (x, t) in row.iter_mut().enumerate() {
|
||||
let s = level[y * n + x];
|
||||
let share = (s[3] / full).min(1.0);
|
||||
let up = if share < 1.0 {
|
||||
at((x as f32 + 0.5) / 2.0 - 0.5, (y as f32 + 0.5) / 2.0 - 0.5)
|
||||
} else {
|
||||
[0.0; 3]
|
||||
};
|
||||
for c in 0..3 {
|
||||
let mean = if s[3] > 0.0 { s[c] / s[3] } else { 0.0 };
|
||||
t[c] = mean * share + up[c] * (1.0 - share);
|
||||
}
|
||||
}
|
||||
});
|
||||
filled = next;
|
||||
width = n;
|
||||
}
|
||||
filled
|
||||
}
|
||||
|
||||
/// Fills the texels without data from the water around them (`pull_push` over the
|
||||
/// water, green brighter than near infrared), below `seed`'s level when given: the
|
||||
/// fill across its group ([`GroupFill`]). Holes in the composite are blocks of
|
||||
/// open water about 10 km across; filled from all the data around them, they took on
|
||||
/// the land across a gulf narrower than that. A cell without water fills from all its
|
||||
/// data.
|
||||
fn fill(bands: &mut Bands, seed: Option<(usize, &[[f32; 3]])>) {
|
||||
if bands.valid.iter().all(|&v| v) {
|
||||
return;
|
||||
}
|
||||
let any_water = bands.valid.iter().zip(&bands.water).any(|(&v, &w)| v && w);
|
||||
let first: Vec<[f32; 4]> = bands
|
||||
.rgb
|
||||
.iter()
|
||||
.zip(&bands.valid)
|
||||
.zip(&bands.water)
|
||||
.map(|((p, &ok), &water)| {
|
||||
if ok && (water || !any_water) {
|
||||
[f32::from(p[0]), f32::from(p[1]), f32::from(p[2]), 1.0]
|
||||
} else {
|
||||
[0.0; 4]
|
||||
}
|
||||
})
|
||||
.collect();
|
||||
let filled = pull_push(bands.size, first, seed);
|
||||
for (i, ok) in bands.valid.iter().enumerate() {
|
||||
if !ok {
|
||||
bands.rgb[i] = filled[i].map(|v| v.round().clamp(1.0, 65535.0) as u16);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// The water of the source at `path` (`cell`) at JOINT texels a side, from its 750 x 750
|
||||
/// overview: per texel, the sums of red, green and blue over its water and its share of
|
||||
/// water, then the same over all its data.
|
||||
fn coarse_bands(path: &Path, cell: Cell) -> Result<[Vec<[f32; 4]>; 2], String> {
|
||||
let fail = |e: &dyn std::fmt::Display| format!("{}: {e}", path.display());
|
||||
let mut decoder = open(path)?;
|
||||
check_place(&mut decoder, path, cell)?;
|
||||
decoder.seek_to_image(4).map_err(|e| fail(&e))?;
|
||||
let size = 2 * JOINT;
|
||||
if decoder.dimensions().map_err(|e| fail(&e))? != (size as u32, size as u32) {
|
||||
return Err(fail(&format!("its fourth overview is not {size} x {size}")));
|
||||
}
|
||||
let values = match decoder.read_image().map_err(|e| fail(&e))? {
|
||||
DecodingResult::U16(v) => v,
|
||||
_ => return Err(fail(&"not 16-bit")),
|
||||
};
|
||||
let mut water = vec![[0.0f32; 4]; JOINT * JOINT];
|
||||
let mut data = vec![[0.0f32; 4]; JOINT * JOINT];
|
||||
for (i, p) in values.chunks_exact(4).enumerate() {
|
||||
if p[..3].contains(&0) {
|
||||
continue;
|
||||
}
|
||||
let at = (i / size / 2) * JOINT + (i % size) / 2;
|
||||
let rgb = [f32::from(p[0]), f32::from(p[1]), f32::from(p[2])];
|
||||
for (sums, counts) in [(&mut data, true), (&mut water, p[1] > p[3])] {
|
||||
if counts {
|
||||
for (sum, v) in sums[at].iter_mut().zip(rgb) {
|
||||
*sum += v / 4.0;
|
||||
}
|
||||
sums[at][3] += 0.25;
|
||||
}
|
||||
}
|
||||
}
|
||||
Ok([water, data])
|
||||
}
|
||||
|
||||
/// The cells joined to `cell` through neighbors whose sources are in `source`, itself
|
||||
/// included, each with its place from `cell` in cells east and north (counted the short
|
||||
/// way round across the 180 degree meridian, from `cell` outward).
|
||||
fn group(source: &Path, cell: Cell) -> Vec<(Cell, (i32, i32))> {
|
||||
let mut found = vec![(cell, (0, 0))];
|
||||
let mut seen = std::collections::HashSet::from([cell]);
|
||||
let mut next = 0;
|
||||
while next < found.len() {
|
||||
let (c, (x, y)) = found[next];
|
||||
next += 1;
|
||||
for &(e, n) in &Cell::NEIGHBORS {
|
||||
if let Some(other) = c.neighbor(e, n) {
|
||||
if source.join(source_name(other)).is_file() && seen.insert(other) {
|
||||
found.push((other, (x + e, y + n)));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
found.sort_by_key(|(c, _)| *c);
|
||||
found
|
||||
}
|
||||
|
||||
/// A key for a group of cells: FNV-1a over their names, in order.
|
||||
fn group_key(cells: &[(Cell, (i32, i32))]) -> u32 {
|
||||
cells
|
||||
.iter()
|
||||
.flat_map(|(c, _)| c.name().into_bytes())
|
||||
.fold(0x811c_9dc5u32, |h, b| {
|
||||
(h ^ u32::from(b)).wrapping_mul(0x0100_0193)
|
||||
})
|
||||
}
|
||||
|
||||
/// A group's water at JOINT texels a cell, filled across the whole group at once.
|
||||
/// Filled one cell at a time, holes took each cell's own water, which differs from cell
|
||||
/// to cell (N35W005's is half as bright as N36W005's), and the change showed along the
|
||||
/// cells' edge across open sea; filled over each cell's 3 x 3 neighbors, the windows of
|
||||
/// neighboring cells still differed, and left a fainter line.
|
||||
struct GroupFill {
|
||||
/// The group's cells with their places, and the westmost and northmost of these.
|
||||
cells: Vec<(Cell, (i32, i32))>,
|
||||
corner: (i32, i32),
|
||||
/// Cells across the grid, and its filled texels in rows from the north-west.
|
||||
width: usize,
|
||||
filled: Vec<[f32; 3]>,
|
||||
}
|
||||
|
||||
impl GroupFill {
|
||||
fn new(source: &Path, cells: Vec<(Cell, (i32, i32))>) -> Result<Self, String> {
|
||||
let west = cells.iter().map(|(_, (x, _))| *x).min().expect("one cell");
|
||||
let east = cells.iter().map(|(_, (x, _))| *x).max().expect("one cell");
|
||||
let south = cells.iter().map(|(_, (_, y))| *y).min().expect("one cell");
|
||||
let north = cells.iter().map(|(_, (_, y))| *y).max().expect("one cell");
|
||||
let (width, height) = ((east - west + 1) as usize, (north - south + 1) as usize);
|
||||
// Square, as `pull_push` takes it; rows past the group's hold nothing.
|
||||
let size = width.max(height) * JOINT;
|
||||
let mut mosaic = [
|
||||
vec![[0.0f32; 4]; size * size],
|
||||
vec![[0.0f32; 4]; size * size],
|
||||
];
|
||||
for &(c, (x, y)) in &cells {
|
||||
let bands = coarse_bands(&source.join(source_name(c)), c)?;
|
||||
let (bx, by) = ((x - west) as usize * JOINT, (north - y) as usize * JOINT);
|
||||
for (into, from) in mosaic.iter_mut().zip(&bands) {
|
||||
for row in 0..JOINT {
|
||||
let to = (by + row) * size + bx;
|
||||
into[to..to + JOINT].copy_from_slice(&from[row * JOINT..(row + 1) * JOINT]);
|
||||
}
|
||||
}
|
||||
}
|
||||
let [water, data] = mosaic;
|
||||
let first = if water.iter().any(|t| t[3] > 0.0) {
|
||||
water
|
||||
} else {
|
||||
data
|
||||
};
|
||||
Ok(Self {
|
||||
cells,
|
||||
corner: (west, north),
|
||||
width: size / JOINT,
|
||||
filled: pull_push(size, first, None),
|
||||
})
|
||||
}
|
||||
|
||||
/// The filled texels over `cell`, JOINT a side.
|
||||
fn seed(&self, cell: Cell) -> Vec<[f32; 3]> {
|
||||
let (_, (x, y)) = *self
|
||||
.cells
|
||||
.iter()
|
||||
.find(|(c, _)| *c == cell)
|
||||
.expect("in the group");
|
||||
let size = self.width * JOINT;
|
||||
let (bx, by) = (
|
||||
(x - self.corner.0) as usize * JOINT,
|
||||
(self.corner.1 - y) as usize * JOINT,
|
||||
);
|
||||
(0..JOINT)
|
||||
.flat_map(|row| {
|
||||
let from = (by + row) * size + bx;
|
||||
self.filled[from..from + JOINT].to_vec()
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
}
|
||||
|
||||
/// A cell's colors in REGION_COARSE blocks a side, as reflectance, from its 750 x 750
|
||||
/// overview, with the texels it has no data for taken from `seed`, its holes as its
|
||||
/// group fills them at JOINT texels a side.
|
||||
fn overview_coarse(path: &Path, cell: Cell, seed: &[[f32; 3]]) -> Result<Vec<[f32; 3]>, String> {
|
||||
let fail = |e: &dyn std::fmt::Display| format!("{}: {e}", path.display());
|
||||
let mut decoder = open(path)?;
|
||||
check_place(&mut decoder, path, cell)?;
|
||||
decoder.seek_to_image(4).map_err(|e| fail(&e))?;
|
||||
let size = 2 * JOINT;
|
||||
let values = match decoder.read_image().map_err(|e| fail(&e))? {
|
||||
DecodingResult::U16(v) => v,
|
||||
_ => return Err(fail(&"not 16-bit")),
|
||||
};
|
||||
let block = size / REGION_COARSE as usize;
|
||||
let n = REGION_COARSE as usize;
|
||||
let mut sums = vec![[0.0f64; 3]; n * n];
|
||||
for (i, p) in values.chunks_exact(4).enumerate() {
|
||||
let (x, y) = (i % size, i / size);
|
||||
let rgb = if p[..3].contains(&0) {
|
||||
seed[(y / 2) * JOINT + x / 2]
|
||||
} else {
|
||||
[f32::from(p[0]), f32::from(p[1]), f32::from(p[2])]
|
||||
};
|
||||
let at = (y / block) * n + x / block;
|
||||
for (sum, v) in sums[at].iter_mut().zip(rgb) {
|
||||
*sum += f64::from(v);
|
||||
}
|
||||
}
|
||||
let area = (block * block) as f64;
|
||||
Ok(sums
|
||||
.iter()
|
||||
.map(|s| s.map(|v| (v / area) as f32 * SCALE))
|
||||
.collect())
|
||||
}
|
||||
|
||||
/// `cell`'s coarse map ([`EarthDetail::coarse`]): its filled texels' means in
|
||||
/// REGION_COARSE blocks a side, and a border of one block each side from the neighbors
|
||||
/// in `neighbors` ([`overview_coarse`]), or its own edge where none is on hand.
|
||||
fn coarse_map(
|
||||
source: &Path,
|
||||
cell: Cell,
|
||||
neighbors: u8,
|
||||
bands: &Bands,
|
||||
group: &GroupFill,
|
||||
) -> Result<Vec<[f32; 3]>, String> {
|
||||
let n = REGION_COARSE as usize;
|
||||
let block = SIZE / n;
|
||||
let mut interior = vec![[0.0f32; 3]; n * n];
|
||||
resolved::earth::parallel_rows(&mut interior, n, |by, row| {
|
||||
for (bx, t) in row.iter_mut().enumerate() {
|
||||
let mut sum = [0.0f64; 3];
|
||||
for y in by * block..(by + 1) * block {
|
||||
for p in &bands.rgb[y * SIZE + bx * block..y * SIZE + (bx + 1) * block] {
|
||||
for c in 0..3 {
|
||||
sum[c] += f64::from(p[c]);
|
||||
}
|
||||
}
|
||||
}
|
||||
*t = sum.map(|s| (s / (block * block) as f64) as f32 * SCALE);
|
||||
}
|
||||
});
|
||||
let mut around = std::collections::HashMap::new();
|
||||
for (i, &(e, nn)) in Cell::NEIGHBORS.iter().enumerate() {
|
||||
if neighbors & (1 << i) != 0 {
|
||||
let other = cell.neighbor(e, nn).expect("on hand");
|
||||
let path = source.join(source_name(other));
|
||||
around.insert((e, nn), overview_coarse(&path, other, &group.seed(other))?);
|
||||
}
|
||||
}
|
||||
let padded = REGION_COARSE_PADDED as usize;
|
||||
let mut out = vec![[0.0f32; 3]; padded * padded];
|
||||
for py in 0..padded {
|
||||
for px in 0..padded {
|
||||
let (x, y) = (px as i64 - 1, py as i64 - 1);
|
||||
let step = |v: i64| match v {
|
||||
-1 => -1,
|
||||
v if v == n as i64 => 1,
|
||||
_ => 0,
|
||||
};
|
||||
let (e, s) = (step(x), step(y));
|
||||
out[py * padded + px] = match around.get(&(e, -s)) {
|
||||
Some(map) if (e, s) != (0, 0) => {
|
||||
let (ox, oy) = (x - i64::from(e) * n as i64, y - i64::from(s) * n as i64);
|
||||
map[oy as usize * n + ox as usize]
|
||||
}
|
||||
_ => {
|
||||
let (cx, cy) = (x.clamp(0, n as i64 - 1), y.clamp(0, n as i64 - 1));
|
||||
interior[cy as usize * n + cx as usize]
|
||||
}
|
||||
};
|
||||
}
|
||||
}
|
||||
Ok(out)
|
||||
}
|
||||
|
||||
/// Which neighbors of `cell` have their sources in `source`, a bit each in the order of
|
||||
/// [`Cell::NEIGHBORS`].
|
||||
pub fn neighbors_on_hand(source: &Path, cell: Cell) -> u8 {
|
||||
Cell::NEIGHBORS
|
||||
.iter()
|
||||
.enumerate()
|
||||
.filter(|(_, (e, n))| {
|
||||
cell.neighbor(*e, *n)
|
||||
.is_some_and(|c| source.join(source_name(c)).is_file())
|
||||
})
|
||||
.fold(0, |bits, (i, _)| bits | 1 << i)
|
||||
}
|
||||
|
||||
/// Whether `out` holds `cell` built with `neighbors`' texels in its border and its
|
||||
/// holes filled with the group `fill` stands for.
|
||||
fn up_to_date(out: &Path, cell: Cell, neighbors: u8, fill: u32) -> bool {
|
||||
EarthDetail::open(out).is_ok_and(|d| {
|
||||
d.cell() == Some(cell)
|
||||
&& d.neighbors() == neighbors
|
||||
&& d.fill_key() == fill
|
||||
&& d.coarse().is_some()
|
||||
&& d.layout() == DetailLayout::WORLDCOVER
|
||||
&& d.format() == DetailFormat::Bc6h
|
||||
})
|
||||
}
|
||||
|
||||
/// Builds each of `cells` from its source in `source` into `out`, unless it is there
|
||||
/// already with the neighbors and the group of cells on hand now: a cell joining a
|
||||
/// group remakes every cell of it, whose holes are filled across it. Cells without a
|
||||
/// source are skipped.
|
||||
pub fn build(source: &Path, out: &Path, cells: &[Cell]) -> Result<(), String> {
|
||||
std::fs::create_dir_all(out).map_err(|e| format!("{}: {e}", out.display()))?;
|
||||
let mut fills = std::collections::HashMap::new();
|
||||
for &cell in cells {
|
||||
let from = source.join(source_name(cell));
|
||||
if !from.is_file() {
|
||||
eprintln!("{cell}: no source in {}, skipped", source.display());
|
||||
continue;
|
||||
}
|
||||
let neighbors = neighbors_on_hand(source, cell);
|
||||
let members = group(source, cell);
|
||||
let key = group_key(&members);
|
||||
let to = out.join(tiles_name(cell));
|
||||
if up_to_date(&to, cell, neighbors, key) {
|
||||
continue;
|
||||
}
|
||||
if let std::collections::hash_map::Entry::Vacant(slot) = fills.entry(key) {
|
||||
eprintln!("filling holes across {} cells", members.len());
|
||||
slot.insert(GroupFill::new(source, members)?);
|
||||
}
|
||||
build_cell(source, cell, neighbors, (key, &fills[&key]), &to)?;
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// A neighbor's texels along this cell's edge: columns `span[0]..span[1]` and rows
|
||||
/// `span[2]..span[3]` in this cell's own, `None` where it has no data.
|
||||
struct Strip {
|
||||
from: Cell,
|
||||
span: [i64; 4],
|
||||
texels: Vec<Option<[f32; 3]>>,
|
||||
}
|
||||
|
||||
fn build_cell(
|
||||
source: &Path,
|
||||
cell: Cell,
|
||||
neighbors: u8,
|
||||
(key, group): (u32, &GroupFill),
|
||||
out: &PathBuf,
|
||||
) -> Result<(), String> {
|
||||
let started = Instant::now();
|
||||
let layout = DetailLayout::WORLDCOVER;
|
||||
let g = layout.gutter as usize;
|
||||
eprintln!("{cell}: reading");
|
||||
let mut bands = decode(&source.join(source_name(cell)), cell)?;
|
||||
let holes = bands.valid.iter().filter(|&&v| !v).count();
|
||||
if holes > 0 {
|
||||
let seed = group.seed(cell);
|
||||
fill(&mut bands, Some((JOINT, &seed)));
|
||||
}
|
||||
let coarse = coarse_map(source, cell, neighbors, &bands, group)?;
|
||||
|
||||
// The neighbors' strips along this cell's edges, as each border needs them.
|
||||
let mut strips = Vec::new();
|
||||
for (i, &(e, n)) in Cell::NEIGHBORS.iter().enumerate() {
|
||||
if neighbors & (1 << i) == 0 {
|
||||
continue;
|
||||
}
|
||||
let other = cell.neighbor(e, n).expect("on hand, so it exists");
|
||||
// The strip in the neighbor's own columns and rows (rows from its north).
|
||||
let span = |step: i32| match step {
|
||||
-1 => (SIZE - g, SIZE),
|
||||
0 => (0, SIZE),
|
||||
_ => (0, g),
|
||||
};
|
||||
let (x0, x1) = span(e);
|
||||
let (y0, y1) = span(-n);
|
||||
let texels = read_strip(&source.join(source_name(other)), other, [x0, x1, y0, y1])?;
|
||||
// The same strip in this cell's columns and rows.
|
||||
let shift = |step: i32| i64::from(step) * SIZE as i64;
|
||||
let (sx, sy) = (shift(e), shift(-n));
|
||||
strips.push(Strip {
|
||||
from: other,
|
||||
span: [
|
||||
x0 as i64 + sx,
|
||||
x1 as i64 + sx,
|
||||
y0 as i64 + sy,
|
||||
y1 as i64 + sy,
|
||||
],
|
||||
texels,
|
||||
});
|
||||
}
|
||||
|
||||
let own = |x: i64, y: i64| -> [f32; 3] {
|
||||
let (x, y) = (
|
||||
x.clamp(0, SIZE as i64 - 1) as usize,
|
||||
y.clamp(0, SIZE as i64 - 1) as usize,
|
||||
);
|
||||
bands.rgb[y * SIZE + x].map(|v| f32::from(v) * SCALE)
|
||||
};
|
||||
let texel = |x: i64, y: i64| -> [f32; 3] {
|
||||
let inside = |v: i64| (0..SIZE as i64).contains(&v);
|
||||
if inside(x) && inside(y) {
|
||||
return own(x, y);
|
||||
}
|
||||
for strip in &strips {
|
||||
let [x0, x1, y0, y1] = strip.span;
|
||||
if (x0..x1).contains(&x) && (y0..y1).contains(&y) {
|
||||
let w = (x1 - x0) as usize;
|
||||
let i = (y - y0) as usize * w + (x - x0) as usize;
|
||||
if let Some(rgb) = strip.texels[i] {
|
||||
return rgb;
|
||||
}
|
||||
}
|
||||
}
|
||||
own(x, y)
|
||||
};
|
||||
eprintln!("{cell}: tiling");
|
||||
let levels = EarthDetail::build_levels(&layout, texel);
|
||||
let fail = |e: std::io::Error| format!("{}: {e}", out.display());
|
||||
let part = out.with_extension("part");
|
||||
let mut file = std::fs::File::create(&part).map_err(fail)?;
|
||||
EarthDetail::write_region_header(&mut file, layout, DetailFormat::Bc6h, cell, neighbors, key)
|
||||
.map_err(fail)?;
|
||||
for (l, level) in levels.iter().enumerate() {
|
||||
file.write_all(&compress(level, layout.padded(l as u32)))
|
||||
.map_err(fail)?;
|
||||
}
|
||||
EarthDetail::write_coarse(&mut file, &coarse).map_err(fail)?;
|
||||
file.sync_all().map_err(fail)?;
|
||||
drop(file);
|
||||
std::fs::rename(&part, out).map_err(fail)?;
|
||||
EarthDetail::open(out).map_err(|e| e.to_string())?;
|
||||
let joined: Vec<String> = strips.iter().map(|s| s.from.to_string()).collect();
|
||||
eprintln!(
|
||||
"{cell}: wrote {} in {:.0} s; {:.2} % of texels filled{}",
|
||||
out.display(),
|
||||
started.elapsed().as_secs_f64(),
|
||||
100.0 * holes as f64 / (SIZE * SIZE) as f64,
|
||||
if joined.is_empty() {
|
||||
String::new()
|
||||
} else {
|
||||
format!("; border from {}", joined.join(", "))
|
||||
}
|
||||
);
|
||||
Ok(())
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn names_follow_esa() {
|
||||
let cell = Cell::parse("N29E032").unwrap();
|
||||
assert_eq!(
|
||||
source_name(cell),
|
||||
"ESA_WorldCover_10m_2021_v200_N29E032_S2RGBNIR.tif"
|
||||
);
|
||||
assert_eq!(tiles_name(cell), "N29E032.tiles");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn holes_fill_with_the_water_around_them() {
|
||||
// Land (bright, near infrared above green) in the west quarter, water elsewhere,
|
||||
// and a hole in the water 8 texels across: it fills with water, untouched data
|
||||
// stays as it was, and land does not bleed in.
|
||||
let size = 32;
|
||||
let land = [3000, 2500, 2000];
|
||||
let water = [300, 500, 600];
|
||||
let mut rgb = Vec::new();
|
||||
let mut valid = Vec::new();
|
||||
let mut is_water = Vec::new();
|
||||
for y in 0..size {
|
||||
for x in 0..size {
|
||||
let wet = x >= size / 4;
|
||||
rgb.push(if wet { water } else { land });
|
||||
let hole = (12..20).contains(&x) && (12..20).contains(&y);
|
||||
valid.push(!hole);
|
||||
is_water.push(wet);
|
||||
}
|
||||
}
|
||||
let before = rgb.clone();
|
||||
let mut bands = Bands {
|
||||
size,
|
||||
rgb,
|
||||
valid,
|
||||
water: is_water,
|
||||
};
|
||||
fill(&mut bands, None);
|
||||
for (i, (now, was)) in bands.rgb.iter().zip(&before).enumerate() {
|
||||
if bands.valid[i] {
|
||||
assert_eq!(now, was, "texel {i} with data changed");
|
||||
} else {
|
||||
for c in 0..3 {
|
||||
let off = (f64::from(now[c]) / f64::from(water[c]) - 1.0).abs();
|
||||
assert!(off < 0.01, "texel {i} channel {c}: {now:?}");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -90,7 +90,7 @@ impl Edges {
|
||||
|
||||
/// BC6H blocks of a `size` x `size` RGBA f16 image, its edges repeated out to a
|
||||
/// multiple of 4, compressed in bands on every core.
|
||||
fn compress(texels: &[f16], size: u32) -> Vec<u8> {
|
||||
pub fn compress(texels: &[f16], size: u32) -> Vec<u8> {
|
||||
let padded = size.div_ceil(4) * 4;
|
||||
let mut image = vec![f16::ZERO; 4 * (padded * padded) as usize];
|
||||
for y in 0..padded {
|
||||
|
||||
@@ -1,5 +1,6 @@
|
||||
//! The `resolved` command: render a scene file to PNG without writing a host program.
|
||||
|
||||
mod earth_regions;
|
||||
mod earth_tiles;
|
||||
mod moon_tiles;
|
||||
mod pose;
|
||||
@@ -14,8 +15,8 @@ use clap::{Args, Parser, Subcommand};
|
||||
use resolved::gpu::AdapterChoice;
|
||||
use resolved::terrain::MoonTerrain;
|
||||
use resolved::{
|
||||
CameraTier, CloudMap, EarthDetail, EarthTexture, Fidelity, Frame, LightingTier, Look,
|
||||
MoonTexture, Preview, Renderer, RendererOptions, StarCatalog, StarPhotometry,
|
||||
CameraTier, CloudMap, EarthDetail, EarthRegions, EarthTexture, Fidelity, Frame, LightingTier,
|
||||
Look, MoonTexture, Preview, Renderer, RendererOptions, StarCatalog, StarPhotometry,
|
||||
};
|
||||
|
||||
use crate::scene_file::SceneFile;
|
||||
@@ -49,6 +50,21 @@ enum Command {
|
||||
#[arg(long)]
|
||||
out: PathBuf,
|
||||
},
|
||||
/// Make 1 x 1 degree cells of Earth's surface at 10 m from ESA WorldCover's
|
||||
/// Sentinel-2 composite (scripts/fetch-earth-regions.py runs this): 194 MB of BC6H
|
||||
/// a cell. Cells already made with the neighbors and group of cells on hand now are
|
||||
/// left alone.
|
||||
EarthRegions {
|
||||
/// Directory with ESA_WorldCover_10m_2021_v200_<cell>_S2RGBNIR.tif.
|
||||
#[arg(long)]
|
||||
source: PathBuf,
|
||||
/// Directory to write <cell>.tiles into.
|
||||
#[arg(long)]
|
||||
out: PathBuf,
|
||||
/// The cells, as N29E032.
|
||||
#[arg(required = true)]
|
||||
cells: Vec<String>,
|
||||
},
|
||||
/// Make the Moon's tiled terrain from LOLA's 128 pixel-per-degree height map
|
||||
/// (scripts/fetch-moon.py --terrain runs this): about 3.3 GB.
|
||||
MoonTiles {
|
||||
@@ -176,6 +192,11 @@ fn main() -> ExitCode {
|
||||
}
|
||||
Command::Render(args) => render(&args),
|
||||
Command::EarthTiles { source, month, out } => earth_tiles::build(&source, month, &out),
|
||||
Command::EarthRegions { source, out, cells } => cells
|
||||
.iter()
|
||||
.map(|c| resolved::earth::Cell::parse(c).map_err(|e| e.to_string()))
|
||||
.collect::<Result<Vec<_>, _>>()
|
||||
.and_then(|cells| earth_regions::build(&source, &out, &cells)),
|
||||
Command::MoonTiles { source, out } => moon_tiles::build(&source, &out),
|
||||
Command::Pose {
|
||||
scene,
|
||||
@@ -243,6 +264,46 @@ fn data_file(data: Option<&Path>, name: &str, script: &str, tier: &str) -> Resul
|
||||
})
|
||||
}
|
||||
|
||||
/// The fine regions `file` names, as far as their cells have been made
|
||||
/// (scripts/fetch-earth-regions.py), from the first data directory that has each. Cells
|
||||
/// not made are drawn from the sharp map, and named on stderr; cells over open water
|
||||
/// have no data, which the script marks with a `.none` file.
|
||||
fn earth_regions(
|
||||
file: &SceneFile,
|
||||
scene: &Path,
|
||||
data: Option<&Path>,
|
||||
) -> Result<Option<Arc<EarthRegions>>, String> {
|
||||
let cells = file.earth_cells()?;
|
||||
let find = |name: String| {
|
||||
data_directories(data)
|
||||
.into_iter()
|
||||
.map(|dir| dir.join("earth/regions").join(&name))
|
||||
.find(|p| p.is_file())
|
||||
};
|
||||
let mut maps = Vec::new();
|
||||
let mut missing = Vec::new();
|
||||
for cell in cells {
|
||||
if let Some(path) = find(earth_regions::tiles_name(cell)) {
|
||||
maps.push(EarthDetail::open(&path).map_err(|e| e.to_string())?);
|
||||
} else if find(format!("{cell}.none")).is_none() {
|
||||
missing.push(cell.name());
|
||||
}
|
||||
}
|
||||
if !missing.is_empty() {
|
||||
eprintln!(
|
||||
"Earth's fine regions: {} not fetched, drawn from the sharp map (scripts/fetch-earth-regions.py --scene {})",
|
||||
missing.join(", "),
|
||||
scene.display()
|
||||
);
|
||||
}
|
||||
if maps.is_empty() {
|
||||
return Ok(None);
|
||||
}
|
||||
Ok(Some(Arc::new(
|
||||
EarthRegions::new(maps).map_err(|e| e.to_string())?,
|
||||
)))
|
||||
}
|
||||
|
||||
fn render(args: &RenderArgs) -> Result<(), String> {
|
||||
let (out, preview, timings) = (args.out.as_deref(), args.preview.as_deref(), args.timings);
|
||||
if out.is_none() && preview.is_none() && timings.is_none() {
|
||||
@@ -283,6 +344,11 @@ fn render(args: &RenderArgs) -> Result<(), String> {
|
||||
} else {
|
||||
None
|
||||
};
|
||||
let earth_regions = if args.lighting >= LightingTier::L1 {
|
||||
earth_regions(&file, path, args.data.as_deref())?
|
||||
} else {
|
||||
None
|
||||
};
|
||||
// The Moon's terrain, if made (scripts/fetch-moon.py --terrain): optional.
|
||||
let moon_terrain = if args.lighting >= LightingTier::L2 {
|
||||
match data_directories(args.data.as_deref())
|
||||
@@ -381,6 +447,7 @@ fn render(args: &RenderArgs) -> Result<(), String> {
|
||||
},
|
||||
earth_texture,
|
||||
earth_detail: earth_detail.clone(),
|
||||
earth_regions: earth_regions.clone(),
|
||||
moon_terrain: moon_terrain.clone(),
|
||||
moon_texture,
|
||||
clouds,
|
||||
@@ -430,11 +497,25 @@ fn render(args: &RenderArgs) -> Result<(), String> {
|
||||
println!(
|
||||
"maps {:.1} ms (Earth and Moon image decode){}",
|
||||
texturing.as_secs_f64() * 1e3,
|
||||
match (earth_detail.is_some(), moon_terrain.is_some()) {
|
||||
(true, true) => "; Earth's sharp map and the Moon's terrain in use",
|
||||
(true, false) => "; Earth's sharp tiled map in use",
|
||||
(false, true) => "; the Moon's terrain in use",
|
||||
(false, false) => "",
|
||||
{
|
||||
let mut used = Vec::new();
|
||||
if earth_detail.is_some() {
|
||||
used.push("Earth's sharp map".to_owned());
|
||||
}
|
||||
if let Some(regions) = &earth_regions {
|
||||
used.push(format!(
|
||||
"{} cells of Earth's fine regions",
|
||||
regions.maps().len()
|
||||
));
|
||||
}
|
||||
if moon_terrain.is_some() {
|
||||
used.push("the Moon's terrain".to_owned());
|
||||
}
|
||||
if used.is_empty() {
|
||||
String::new()
|
||||
} else {
|
||||
format!("; {} in use", used.join(", "))
|
||||
}
|
||||
}
|
||||
);
|
||||
println!(
|
||||
|
||||
@@ -9,6 +9,7 @@ use std::path::{Path, PathBuf};
|
||||
|
||||
use resolved::bodies::Appearance;
|
||||
use resolved::brdf::{Brdf, DIELECTRIC_F0};
|
||||
use resolved::earth::Cell;
|
||||
use resolved::library::spectrum_from_csv;
|
||||
use resolved::model::rgb_spectrum;
|
||||
use resolved::nalgebra::{Quaternion, UnitQuaternion, Vector2, Vector3};
|
||||
@@ -204,6 +205,11 @@ pub struct AppearanceFile {
|
||||
/// The day (YYYY-MM-DD, UTC) whose clouds lighting L2 draws, for scenes without an
|
||||
/// epoch; with one, the epoch's day is the default.
|
||||
pub cloud_date: Option<String>,
|
||||
/// Boxes of ground, each [south, west, north, east] in degrees, where lighting L1
|
||||
/// draws Earth's surface from ESA WorldCover's 10 m cells, once they are fetched
|
||||
/// (scripts/fetch-earth-regions.py --scene).
|
||||
#[serde(default)]
|
||||
pub earth_regions: Vec<[f64; 4]>,
|
||||
}
|
||||
|
||||
/// One spacecraft.
|
||||
@@ -497,6 +503,20 @@ impl SceneFile {
|
||||
}
|
||||
}
|
||||
|
||||
/// The 1-degree cells that cover `earth_regions`, each once.
|
||||
pub fn earth_cells(&self) -> Result<Vec<Cell>> {
|
||||
let mut cells = Vec::new();
|
||||
for &[south, west, north, east] in &self.appearance.earth_regions {
|
||||
cells.extend(
|
||||
Cell::covering(south, west, north, east)
|
||||
.map_err(|e| format!("earth_regions: {e}"))?,
|
||||
);
|
||||
}
|
||||
cells.sort();
|
||||
cells.dedup();
|
||||
Ok(cells)
|
||||
}
|
||||
|
||||
/// The day (YYYY-MM-DD) whose clouds lighting L2 draws: `cloud_date`, else the
|
||||
/// epoch's day in UTC.
|
||||
pub fn cloud_date(&self) -> Result<Option<String>> {
|
||||
@@ -1040,6 +1060,32 @@ mod tests {
|
||||
std::fs::remove_dir_all(&dir).unwrap();
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn earth_regions_name_their_cells() {
|
||||
let with = |extra: &str| -> SceneFile {
|
||||
toml::from_str(&format!("{MINIMAL}\n[appearance]\n{extra}\n")).unwrap()
|
||||
};
|
||||
let names = |file: SceneFile| -> Vec<String> {
|
||||
file.earth_cells().unwrap().iter().map(Cell::name).collect()
|
||||
};
|
||||
assert!(names(with("")).is_empty());
|
||||
// Overlapping boxes name each cell once; edges on whole degrees take no more.
|
||||
assert_eq!(
|
||||
names(with(
|
||||
"earth_regions = [[29, 32, 30.5, 34], [29.2, 33.1, 29.4, 33.5]]"
|
||||
)),
|
||||
["N29E032", "N29E033", "N30E032", "N30E033"]
|
||||
);
|
||||
// Across the 180 degree meridian, and south of the equator.
|
||||
assert_eq!(
|
||||
names(with("earth_regions = [[-1.5, 179.5, -0.5, -179.5]]")),
|
||||
["S02W180", "S02E179", "S01W180", "S01E179"]
|
||||
);
|
||||
assert!(with("earth_regions = [[30, 32, 29, 34]]")
|
||||
.earth_cells()
|
||||
.is_err());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn earth_month_comes_from_the_setting() {
|
||||
let with = |extra: &str| -> SceneFile {
|
||||
|
||||
+501
-14
@@ -117,11 +117,7 @@ impl EarthTexture {
|
||||
|
||||
/// Runs `f(row_index, row)` over the rows of `data` (rows of `row_len` values) on all
|
||||
/// cores.
|
||||
pub(crate) fn parallel_rows<T: Send>(
|
||||
data: &mut [T],
|
||||
row_len: usize,
|
||||
f: impl Fn(usize, &mut [T]) + Sync,
|
||||
) {
|
||||
pub fn parallel_rows<T: Send>(data: &mut [T], row_len: usize, f: impl Fn(usize, &mut [T]) + Sync) {
|
||||
let rows = data.len() / row_len;
|
||||
let threads = std::thread::available_parallelism()
|
||||
.map_or(1, |n| n.get())
|
||||
@@ -171,6 +167,17 @@ impl DetailLayout {
|
||||
levels: 3,
|
||||
};
|
||||
|
||||
/// One 1 x 1 degree cell of ESA WorldCover's Sentinel-2 composite (an
|
||||
/// [`EarthRegions`] cell): 12000 x 12000 (0.3 arcseconds, about 9.3 m at the
|
||||
/// equator), with levels from 9.3 m to 296 m: coarser than that, Blue Marble's sharp
|
||||
/// map is as sharp. Cell edges fall on texel edges at every level.
|
||||
pub const WORLDCOVER: Self = Self {
|
||||
tiles: [1, 1],
|
||||
tile: 12000,
|
||||
gutter: 32,
|
||||
levels: 6,
|
||||
};
|
||||
|
||||
/// Texels across and down the whole map at full size.
|
||||
pub fn size(&self) -> [u32; 2] {
|
||||
[self.tiles[0] * self.tile, self.tiles[1] * self.tile]
|
||||
@@ -252,8 +259,23 @@ pub struct EarthDetail {
|
||||
layout: DetailLayout,
|
||||
format: DetailFormat,
|
||||
data: DetailData,
|
||||
/// For one cell of [`EarthRegions`]: the cell, which of its neighbors' texels fill
|
||||
/// its border ([`Cell::NEIGHBORS`], a bit each), and a key for the cells its holes
|
||||
/// were filled together with.
|
||||
region: Option<(Cell, u8, u32)>,
|
||||
/// For one cell of [`EarthRegions`], its colors in [`REGION_COARSE`] blocks a side
|
||||
/// ([`EarthDetail::coarse`]).
|
||||
coarse: Option<std::sync::Arc<Vec<[f32; 3]>>>,
|
||||
}
|
||||
|
||||
/// Blocks across a cell of [`EarthRegions`] in its coarse map: 4.4 km for ESA
|
||||
/// WorldCover's cells. The renderer brings each cell's colors at this scale to the
|
||||
/// month's map's, and keeps the cell's finer detail.
|
||||
pub const REGION_COARSE: u32 = 25;
|
||||
|
||||
/// Texels across a cell's coarse map: [`REGION_COARSE`] and a border of one each side.
|
||||
pub const REGION_COARSE_PADDED: u32 = REGION_COARSE + 2;
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
enum DetailData {
|
||||
Memory(std::sync::Arc<Vec<u8>>),
|
||||
@@ -292,18 +314,53 @@ impl EarthDetail {
|
||||
levels: word(6),
|
||||
};
|
||||
layout.validate(format)?;
|
||||
let detail = Self {
|
||||
let region = match word(7) {
|
||||
0 => None,
|
||||
1 => Some((
|
||||
Cell::new(word(8) as i32, word(9) as i32)
|
||||
.map_err(|e| fail(format!("header: {e}")))?,
|
||||
word(10) as u8,
|
||||
word(11),
|
||||
)),
|
||||
n => return Err(fail(format!("unknown map kind {n}"))),
|
||||
};
|
||||
let coarse_bytes = match (region.is_some(), word(12)) {
|
||||
(_, 0) => 0,
|
||||
(true, n) if n == REGION_COARSE_PADDED => 12 * u64::from(n * n),
|
||||
(_, n) => return Err(fail(format!("a coarse map of {n} texels across"))),
|
||||
};
|
||||
let mut detail = Self {
|
||||
layout,
|
||||
format,
|
||||
data: DetailData::File(path.to_owned()),
|
||||
region,
|
||||
coarse: None,
|
||||
};
|
||||
let want = DETAIL_HEADER + detail.tile_bytes() * u64::from(layout.layers());
|
||||
let levels = DETAIL_HEADER + detail.tile_bytes() * u64::from(layout.layers());
|
||||
let want = levels + coarse_bytes;
|
||||
let have = file.metadata().map_err(|e| fail(e.to_string()))?.len();
|
||||
if have != want {
|
||||
return Err(fail(format!(
|
||||
"{have} bytes, not the {want} its header implies"
|
||||
)));
|
||||
}
|
||||
if coarse_bytes > 0 {
|
||||
use std::io::{Seek, SeekFrom};
|
||||
file.seek(SeekFrom::Start(levels))
|
||||
.map_err(|e| fail(e.to_string()))?;
|
||||
let mut bytes = vec![0u8; coarse_bytes as usize];
|
||||
file.read_exact(&mut bytes)
|
||||
.map_err(|e| fail(format!("coarse map: {e}")))?;
|
||||
let values: Vec<[f32; 3]> = bytes
|
||||
.chunks_exact(12)
|
||||
.map(|c| {
|
||||
std::array::from_fn(|i| {
|
||||
f32::from_le_bytes(c[4 * i..4 * i + 4].try_into().unwrap())
|
||||
})
|
||||
})
|
||||
.collect();
|
||||
detail.coarse = Some(std::sync::Arc::new(values));
|
||||
}
|
||||
Ok(detail)
|
||||
}
|
||||
|
||||
@@ -318,6 +375,8 @@ impl EarthDetail {
|
||||
layout,
|
||||
format,
|
||||
data: DetailData::Memory(std::sync::Arc::new(bytes)),
|
||||
region: None,
|
||||
coarse: None,
|
||||
};
|
||||
let want = detail.tile_bytes() * u64::from(layout.layers());
|
||||
let have = match &detail.data {
|
||||
@@ -338,9 +397,45 @@ impl EarthDetail {
|
||||
out: &mut impl std::io::Write,
|
||||
layout: DetailLayout,
|
||||
format: DetailFormat,
|
||||
) -> std::io::Result<()> {
|
||||
Self::write_any_header(out, layout, format, None)
|
||||
}
|
||||
|
||||
/// Writes the header of a file holding one cell of [`EarthRegions`], whose border
|
||||
/// holds the texels of the neighbors in `neighbors` ([`Cell::NEIGHBORS`], a bit
|
||||
/// each) and repeats its own edge toward the others, and whose holes were filled
|
||||
/// together with the cells `fill` stands for (a key its maker chooses). Its tile
|
||||
/// follows, then its coarse map ([`EarthDetail::write_coarse`]).
|
||||
pub fn write_region_header(
|
||||
out: &mut impl std::io::Write,
|
||||
layout: DetailLayout,
|
||||
format: DetailFormat,
|
||||
cell: Cell,
|
||||
neighbors: u8,
|
||||
fill: u32,
|
||||
) -> std::io::Result<()> {
|
||||
Self::write_any_header(out, layout, format, Some((cell, neighbors, fill)))
|
||||
}
|
||||
|
||||
fn write_any_header(
|
||||
out: &mut impl std::io::Write,
|
||||
layout: DetailLayout,
|
||||
format: DetailFormat,
|
||||
region: Option<(Cell, u8, u32)>,
|
||||
) -> std::io::Result<()> {
|
||||
let mut header = [0u8; DETAIL_HEADER as usize];
|
||||
header[..8].copy_from_slice(DETAIL_MAGIC);
|
||||
let (kind, south, west, neighbors, fill, coarse) = match region {
|
||||
Some((cell, n, f)) => (
|
||||
1,
|
||||
cell.south as u32,
|
||||
cell.west as u32,
|
||||
u32::from(n),
|
||||
f,
|
||||
REGION_COARSE_PADDED,
|
||||
),
|
||||
None => (0, 0, 0, 0, 0, 0),
|
||||
};
|
||||
let words = [
|
||||
1,
|
||||
format.code(),
|
||||
@@ -349,6 +444,12 @@ impl EarthDetail {
|
||||
layout.tile,
|
||||
layout.gutter,
|
||||
layout.levels,
|
||||
kind,
|
||||
south,
|
||||
west,
|
||||
neighbors,
|
||||
fill,
|
||||
coarse,
|
||||
];
|
||||
for (i, w) in words.iter().enumerate() {
|
||||
header[8 + 4 * i..12 + 4 * i].copy_from_slice(&w.to_le_bytes());
|
||||
@@ -366,6 +467,90 @@ impl EarthDetail {
|
||||
self.format
|
||||
}
|
||||
|
||||
/// The cell this map covers, if it is one cell of [`EarthRegions`].
|
||||
pub fn cell(&self) -> Option<Cell> {
|
||||
self.region.map(|(cell, _, _)| cell)
|
||||
}
|
||||
|
||||
/// For one cell of [`EarthRegions`], the neighbors whose texels fill its border
|
||||
/// ([`Cell::NEIGHBORS`], a bit each); 0 for a whole map.
|
||||
pub fn neighbors(&self) -> u8 {
|
||||
self.region.map_or(0, |(_, n, _)| n)
|
||||
}
|
||||
|
||||
/// For one cell of [`EarthRegions`], the key of the cells its holes were filled
|
||||
/// together with; 0 for a whole map.
|
||||
pub fn fill_key(&self) -> u32 {
|
||||
self.region.map_or(0, |(_, _, f)| f)
|
||||
}
|
||||
|
||||
/// The same map as one cell of [`EarthRegions`]: `cell`, with its border holding
|
||||
/// the texels of the neighbors in `neighbors`, and its coarse map `coarse`
|
||||
/// ([`EarthDetail::coarse`]). For maps held in memory; a cell's file says so in its
|
||||
/// header ([`EarthDetail::write_region_header`]).
|
||||
pub fn at_cell(self, cell: Cell, neighbors: u8, coarse: Vec<[f32; 3]>) -> Result<Self, Error> {
|
||||
let n = REGION_COARSE_PADDED as usize;
|
||||
if coarse.len() != n * n {
|
||||
return Err(Error::Texture(format!(
|
||||
"a coarse map of {} texels, not {n} x {n}",
|
||||
coarse.len()
|
||||
)));
|
||||
}
|
||||
Ok(Self {
|
||||
region: Some((cell, neighbors, 0)),
|
||||
coarse: Some(std::sync::Arc::new(coarse)),
|
||||
..self
|
||||
})
|
||||
}
|
||||
|
||||
/// For one cell of [`EarthRegions`], its mean linear red, green and blue in
|
||||
/// [`REGION_COARSE`] x [`REGION_COARSE`] blocks, with a border of one block each
|
||||
/// side from its neighbors: [`REGION_COARSE_PADDED`] a side, in rows from the
|
||||
/// north-west.
|
||||
pub fn coarse(&self) -> Option<&[[f32; 3]]> {
|
||||
self.coarse.as_deref().map(Vec::as_slice)
|
||||
}
|
||||
|
||||
/// Writes a cell's coarse map, after its tile, in a file begun by
|
||||
/// [`EarthDetail::write_region_header`]: f32 values, little-endian.
|
||||
pub fn write_coarse(out: &mut impl std::io::Write, coarse: &[[f32; 3]]) -> std::io::Result<()> {
|
||||
let bytes: Vec<u8> = coarse
|
||||
.iter()
|
||||
.flatten()
|
||||
.flat_map(|v| v.to_le_bytes())
|
||||
.collect();
|
||||
out.write_all(&bytes)
|
||||
}
|
||||
|
||||
/// A cell's coarse map ([`EarthDetail::coarse`]) from `texel(column, row)`, counted
|
||||
/// as in [`EarthDetail::build_levels`] from the interior's north-west corner, over
|
||||
/// the interior and a block beyond it each side. `layout.tile` must be a multiple
|
||||
/// of [`REGION_COARSE`].
|
||||
pub fn coarse_map(
|
||||
layout: &DetailLayout,
|
||||
texel: impl Fn(i64, i64) -> [f32; 3] + Sync,
|
||||
) -> Vec<[f32; 3]> {
|
||||
let block = i64::from(layout.tile / REGION_COARSE);
|
||||
let n = REGION_COARSE_PADDED as usize;
|
||||
let mut out = vec![[0.0f32; 3]; n * n];
|
||||
parallel_rows(&mut out, n, |by, row| {
|
||||
for (bx, t) in row.iter_mut().enumerate() {
|
||||
let (x0, y0) = ((bx as i64 - 1) * block, (by as i64 - 1) * block);
|
||||
let mut sum = [0.0f64; 3];
|
||||
for y in y0..y0 + block {
|
||||
for x in x0..x0 + block {
|
||||
let v = texel(x, y);
|
||||
for c in 0..3 {
|
||||
sum[c] += f64::from(v[c]);
|
||||
}
|
||||
}
|
||||
}
|
||||
*t = sum.map(|s| (s / (block * block) as f64) as f32);
|
||||
}
|
||||
});
|
||||
out
|
||||
}
|
||||
|
||||
/// Bytes of one tile, all its levels.
|
||||
pub fn tile_bytes(&self) -> u64 {
|
||||
(0..self.layout.levels)
|
||||
@@ -412,18 +597,28 @@ impl EarthDetail {
|
||||
texel: impl Fn(u32, u32) -> [f32; 3] + Sync,
|
||||
) -> Vec<Vec<f16>> {
|
||||
let [width, height] = layout.size();
|
||||
let (x0, y0) = (i64::from(x * layout.tile), i64::from(y * layout.tile));
|
||||
Self::build_levels(layout, |col, row| {
|
||||
let gx = (x0 + col).rem_euclid(i64::from(width)) as u32;
|
||||
let gy = (y0 + row).clamp(0, i64::from(height) - 1) as u32;
|
||||
texel(gx, gy)
|
||||
})
|
||||
}
|
||||
|
||||
/// One tile of `layout` at each mip level, as [`EarthDetail::build_tile`] gives
|
||||
/// them, from `texel(column, row)` counted from the tile's interior's north-west
|
||||
/// corner: from `-gutter` to `tile + gutter - 1` each way.
|
||||
pub fn build_levels(
|
||||
layout: &DetailLayout,
|
||||
texel: impl Fn(i64, i64) -> [f32; 3] + Sync,
|
||||
) -> Vec<Vec<f16>> {
|
||||
let size = layout.padded(0) as usize;
|
||||
let one = f16::from_f32(1.0);
|
||||
let mut first = vec![one; 4 * size * size];
|
||||
let (x0, y0) = (
|
||||
i64::from(x * layout.tile) - i64::from(layout.gutter),
|
||||
i64::from(y * layout.tile) - i64::from(layout.gutter),
|
||||
);
|
||||
let gutter = i64::from(layout.gutter);
|
||||
parallel_rows(&mut first, 4 * size, |row, out| {
|
||||
let gy = (y0 + row as i64).clamp(0, i64::from(height) - 1) as u32;
|
||||
for (col, t) in out.chunks_exact_mut(4).enumerate() {
|
||||
let gx = (x0 + col as i64).rem_euclid(i64::from(width)) as u32;
|
||||
let rgb = texel(gx, gy);
|
||||
let rgb = texel(col as i64 - gutter, row as i64 - gutter);
|
||||
for c in 0..3 {
|
||||
t[c] = f16::from_f32(rgb[c]);
|
||||
}
|
||||
@@ -453,6 +648,197 @@ impl EarthDetail {
|
||||
}
|
||||
}
|
||||
|
||||
/// A 1 x 1 degree cell of Earth's surface between whole degrees of geodetic latitude
|
||||
/// and longitude, named by its south-west corner as ESA WorldCover names its tiles:
|
||||
/// N29E032 runs from 29 to 30 degrees north and 32 to 33 east.
|
||||
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
|
||||
pub struct Cell {
|
||||
/// The southern edge's latitude, -90 to 89.
|
||||
pub south: i32,
|
||||
/// The western edge's longitude, -180 to 179.
|
||||
pub west: i32,
|
||||
}
|
||||
|
||||
impl Cell {
|
||||
/// The eight neighbors in the order of the bits of [`EarthDetail::neighbors`], as
|
||||
/// steps east and north: north-west, north, north-east, west, east, south-west,
|
||||
/// south, south-east.
|
||||
pub const NEIGHBORS: [(i32, i32); 8] = [
|
||||
(-1, 1),
|
||||
(0, 1),
|
||||
(1, 1),
|
||||
(-1, 0),
|
||||
(1, 0),
|
||||
(-1, -1),
|
||||
(0, -1),
|
||||
(1, -1),
|
||||
];
|
||||
|
||||
/// The cell whose south-west corner is at `south`, `west` degrees.
|
||||
pub fn new(south: i32, west: i32) -> Result<Self, Error> {
|
||||
if !(-90..90).contains(&south) || !(-180..180).contains(&west) {
|
||||
return Err(Error::Texture(format!(
|
||||
"no 1-degree cell has its south-west corner at {south}, {west}"
|
||||
)));
|
||||
}
|
||||
Ok(Self { south, west })
|
||||
}
|
||||
|
||||
/// Its name, as N29E032 or S01W075.
|
||||
pub fn name(&self) -> String {
|
||||
format!(
|
||||
"{}{:02}{}{:03}",
|
||||
if self.south < 0 { 'S' } else { 'N' },
|
||||
self.south.abs(),
|
||||
if self.west < 0 { 'W' } else { 'E' },
|
||||
self.west.abs()
|
||||
)
|
||||
}
|
||||
|
||||
/// The cell named `name`, as N29E032.
|
||||
pub fn parse(name: &str) -> Result<Self, Error> {
|
||||
let bad = || Error::Texture(format!("{name:?} does not name a cell, as N29E032 does"));
|
||||
let b = name.as_bytes();
|
||||
if b.len() != 7 || !b[1..3].iter().chain(&b[4..]).all(u8::is_ascii_digit) {
|
||||
return Err(bad());
|
||||
}
|
||||
let number = |s: &str| s.parse::<i32>().map_err(|_| bad());
|
||||
let south = match b[0] {
|
||||
b'N' => number(&name[1..3])?,
|
||||
b'S' => -number(&name[1..3])?,
|
||||
_ => return Err(bad()),
|
||||
};
|
||||
let west = match b[3] {
|
||||
b'E' => number(&name[4..])?,
|
||||
b'W' => -number(&name[4..])?,
|
||||
_ => return Err(bad()),
|
||||
};
|
||||
Self::new(south, west)
|
||||
}
|
||||
|
||||
/// The cells that cover the box from `south` to `north` and from `west` east to
|
||||
/// `east`, in degrees. A box whose `west` is east of its `east` crosses the 180
|
||||
/// degree meridian.
|
||||
pub fn covering(south: f64, west: f64, north: f64, east: f64) -> Result<Vec<Self>, Error> {
|
||||
let valid = (-90.0..=90.0).contains(&south)
|
||||
&& (-90.0..=90.0).contains(&north)
|
||||
&& south < north
|
||||
&& (-180.0..=180.0).contains(&west)
|
||||
&& (-180.0..=180.0).contains(&east)
|
||||
&& west != east;
|
||||
if !valid {
|
||||
return Err(Error::Texture(format!(
|
||||
"{south}, {west} to {north}, {east} is not a box of latitude and longitude"
|
||||
)));
|
||||
}
|
||||
let rows = (south.floor() as i32).max(-90)..(north.ceil() as i32).min(90);
|
||||
let first = west.floor() as i32;
|
||||
let mut last = east.ceil() as i32;
|
||||
if last <= first {
|
||||
last += 360;
|
||||
}
|
||||
let mut cells = Vec::new();
|
||||
for south in rows {
|
||||
for west in first..last {
|
||||
let west = (west + 180).rem_euclid(360) - 180;
|
||||
cells.push(Self { south, west });
|
||||
}
|
||||
}
|
||||
cells.sort();
|
||||
cells.dedup();
|
||||
Ok(cells)
|
||||
}
|
||||
|
||||
/// Its neighbor `east` and `north` cells away, wrapping across the 180 degree
|
||||
/// meridian; none past a pole.
|
||||
pub fn neighbor(&self, east: i32, north: i32) -> Option<Self> {
|
||||
let south = self.south + north;
|
||||
(-90..90).contains(&south).then(|| Self {
|
||||
south,
|
||||
west: (self.west + east + 180).rem_euclid(360) - 180,
|
||||
})
|
||||
}
|
||||
|
||||
/// Its index among all 360 x 180 cells, in rows from the north-west, as the shader
|
||||
/// finds it.
|
||||
pub fn index(&self) -> usize {
|
||||
(89 - self.south) as usize * 360 + (self.west + 180) as usize
|
||||
}
|
||||
}
|
||||
|
||||
impl std::fmt::Display for Cell {
|
||||
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
|
||||
f.write_str(&self.name())
|
||||
}
|
||||
}
|
||||
|
||||
/// Fine maps of Earth's surface over chosen 1 x 1 degree [`Cell`]s, drawn where a view
|
||||
/// resolves ground finer than the sharp map does, and blended into it at their edges.
|
||||
/// Each cell is a one-tile [`EarthDetail`] of one layout, made from ESA WorldCover's
|
||||
/// Sentinel-2 composite by `resolved earth-regions`
|
||||
/// (`scripts/fetch-earth-regions.py`), whose files are read as they are uploaded.
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct EarthRegions {
|
||||
cells: Vec<EarthDetail>,
|
||||
}
|
||||
|
||||
impl EarthRegions {
|
||||
/// The regions made of `cells`: one-tile maps of one layout and format, each of a
|
||||
/// different cell.
|
||||
pub fn new(cells: Vec<EarthDetail>) -> Result<Self, Error> {
|
||||
let Some(first) = cells.first() else {
|
||||
return Err(Error::Texture("regions need at least one cell".into()));
|
||||
};
|
||||
let (layout, format) = (first.layout(), first.format());
|
||||
let mut seen = std::collections::HashSet::new();
|
||||
for map in &cells {
|
||||
let Some(cell) = map.cell() else {
|
||||
return Err(Error::Texture(
|
||||
"a whole map is not a cell of the regions".into(),
|
||||
));
|
||||
};
|
||||
if map.coarse().is_none() {
|
||||
return Err(Error::Texture(format!("{cell} has no coarse map")));
|
||||
}
|
||||
if map.layout() != layout
|
||||
|| map.format() != format
|
||||
|| layout.tiles != [1, 1]
|
||||
|| layout.tile % REGION_COARSE != 0
|
||||
{
|
||||
return Err(Error::Texture(format!(
|
||||
"{cell} is {:?} in {:?}, not one tile of {layout:?} in {format:?} as the first",
|
||||
map.layout(),
|
||||
map.format()
|
||||
)));
|
||||
}
|
||||
if !seen.insert(cell) {
|
||||
return Err(Error::Texture(format!("{cell} is given twice")));
|
||||
}
|
||||
}
|
||||
Ok(Self { cells })
|
||||
}
|
||||
|
||||
/// The layout each cell's map has.
|
||||
pub fn layout(&self) -> DetailLayout {
|
||||
self.cells[0].layout()
|
||||
}
|
||||
|
||||
/// The texel format of every cell.
|
||||
pub fn format(&self) -> DetailFormat {
|
||||
self.cells[0].format()
|
||||
}
|
||||
|
||||
/// The cells' maps.
|
||||
pub fn maps(&self) -> &[EarthDetail] {
|
||||
&self.cells
|
||||
}
|
||||
|
||||
/// Bytes of device memory they take.
|
||||
pub fn bytes(&self) -> u64 {
|
||||
self.cells[0].tile_bytes() * self.cells.len() as u64
|
||||
}
|
||||
}
|
||||
|
||||
/// Linear light from an sRGB-encoded value in [0, 1].
|
||||
pub fn srgb_to_linear(v: f32) -> f32 {
|
||||
if v <= 0.040_45 {
|
||||
@@ -574,6 +960,107 @@ mod tests {
|
||||
.is_ok());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn cells_are_named_as_esa_names_them() {
|
||||
for (name, south, west) in [
|
||||
("N29E032", 29, 32),
|
||||
("S01W075", -1, -75),
|
||||
("N00W001", 0, -1),
|
||||
("S90W180", -90, -180),
|
||||
("N89E179", 89, 179),
|
||||
] {
|
||||
let cell = Cell::parse(name).unwrap();
|
||||
assert_eq!((cell.south, cell.west), (south, west));
|
||||
assert_eq!(cell.name(), name);
|
||||
}
|
||||
for bad in ["N90E000", "N00E180", "X29E032", "N29E32", "n29e032", ""] {
|
||||
assert!(Cell::parse(bad).is_err(), "{bad}");
|
||||
}
|
||||
// Rows from the north-west, as the shader counts them.
|
||||
assert_eq!(Cell::parse("N89W180").unwrap().index(), 0);
|
||||
assert_eq!(Cell::parse("S90E179").unwrap().index(), 360 * 180 - 1);
|
||||
let cell = Cell::parse("N00E179").unwrap();
|
||||
assert_eq!(cell.neighbor(1, 0), Cell::parse("N00W180").ok());
|
||||
assert_eq!(Cell::parse("N89E000").unwrap().neighbor(0, 1), None);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn boxes_are_covered_by_whole_cells() {
|
||||
let names = |s, w, n, e| -> Vec<String> {
|
||||
Cell::covering(s, w, n, e)
|
||||
.unwrap()
|
||||
.iter()
|
||||
.map(Cell::name)
|
||||
.collect()
|
||||
};
|
||||
assert_eq!(names(29.0, 32.0, 30.0, 33.0), ["N29E032"]);
|
||||
assert_eq!(
|
||||
names(29.5, 32.5, 30.5, 33.5),
|
||||
["N29E032", "N29E033", "N30E032", "N30E033"]
|
||||
);
|
||||
assert_eq!(names(-0.5, 179.5, 0.5, -179.5).len(), 4);
|
||||
assert_eq!(names(89.5, -180.0, 90.0, 180.0).len(), 360);
|
||||
assert!(Cell::covering(30.0, 32.0, 29.0, 33.0).is_err());
|
||||
assert!(Cell::covering(29.0, 32.0, 30.0, 32.0).is_err());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn region_cells_say_where_they_are() {
|
||||
let layout = DetailLayout {
|
||||
tiles: [1, 1],
|
||||
tile: 50,
|
||||
gutter: 2,
|
||||
levels: 2,
|
||||
};
|
||||
let format = DetailFormat::Rgba16Float;
|
||||
let texel = |c: i64, r: i64| [c as f32, r as f32, 0.0];
|
||||
let mut bytes = Vec::new();
|
||||
for level in EarthDetail::build_levels(&layout, texel) {
|
||||
bytes.extend_from_slice(bytemuck::cast_slice(&level));
|
||||
}
|
||||
// Blocks of 2 x 2 texels: the first, beyond the corner, averages columns and rows
|
||||
// -2 and -1; the last, 48 and 49 inside the far border.
|
||||
let coarse = EarthDetail::coarse_map(&layout, texel);
|
||||
assert_eq!(coarse.len(), 27 * 27);
|
||||
assert_eq!(coarse[0], [-1.5, -1.5, 0.0]);
|
||||
assert_eq!(coarse[27 * 27 - 1], [50.5, 50.5, 0.0]);
|
||||
let cell = Cell::parse("S12W034").unwrap();
|
||||
let path = std::env::temp_dir().join(format!("resolved-cell-{}.bin", std::process::id()));
|
||||
let mut file = std::fs::File::create(&path).unwrap();
|
||||
EarthDetail::write_region_header(&mut file, layout, format, cell, 0b1010_0101, 77).unwrap();
|
||||
std::io::Write::write_all(&mut file, &bytes).unwrap();
|
||||
EarthDetail::write_coarse(&mut file, &coarse).unwrap();
|
||||
drop(file);
|
||||
let opened = EarthDetail::open(&path).unwrap();
|
||||
assert_eq!(opened.cell(), Some(cell));
|
||||
assert_eq!(opened.coarse(), Some(coarse.as_slice()));
|
||||
assert_eq!(opened.neighbors(), 0b1010_0101);
|
||||
assert_eq!(opened.fill_key(), 77);
|
||||
// The border counts from the interior's corner: -2 at the first texel.
|
||||
let first = opened.tile_levels(0).unwrap()[0][..8].to_vec();
|
||||
std::fs::remove_file(&path).unwrap();
|
||||
let texel: &[f16] = bytemuck::cast_slice(&first);
|
||||
assert_eq!((texel[0].to_f32(), texel[1].to_f32()), (-2.0, -2.0));
|
||||
// A whole map is no cell, and regions take only cells of one layout, once each.
|
||||
let whole = EarthDetail::from_bytes(layout, format, bytes.clone()).unwrap();
|
||||
assert_eq!(whole.cell(), None);
|
||||
assert!(EarthRegions::new(vec![whole.clone()]).is_err());
|
||||
let at = |name| {
|
||||
whole
|
||||
.clone()
|
||||
.at_cell(Cell::parse(name).unwrap(), 0, coarse.clone())
|
||||
.unwrap()
|
||||
};
|
||||
assert!(whole.clone().at_cell(cell, 0, vec![]).is_err());
|
||||
assert!(EarthRegions::new(vec![at("N00E000"), at("N00E000")]).is_err());
|
||||
let regions = EarthRegions::new(vec![at("N00E000"), at("N00E001")]).unwrap();
|
||||
assert_eq!(regions.bytes(), 2 * whole.tile_bytes());
|
||||
assert!(EarthRegions::new(vec![]).is_err());
|
||||
assert!(DetailLayout::WORLDCOVER
|
||||
.validate(DetailFormat::Bc6h)
|
||||
.is_ok());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn wrong_sizes_are_rejected() {
|
||||
assert!(EarthTexture::from_rgba8(4, 2, vec![0; 31]).is_err());
|
||||
|
||||
@@ -19,7 +19,10 @@ use crate::atmosphere::{
|
||||
use crate::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS};
|
||||
use crate::camera::Intrinsics;
|
||||
use crate::clouds::{CloudMap, CloudRelief};
|
||||
use crate::earth::{DetailFormat, DetailLayout, EarthDetail, EarthTexture};
|
||||
use crate::earth::{
|
||||
DetailFormat, DetailLayout, EarthDetail, EarthRegions, EarthTexture, REGION_COARSE,
|
||||
REGION_COARSE_PADDED,
|
||||
};
|
||||
use crate::model::rgb_spectrum;
|
||||
use crate::moon::{MoonTexture, NORMALIZED_INCIDENCE, PHASE_SCALE};
|
||||
use crate::radiometry::{Radiometry, AU, SUN_RADIUS};
|
||||
@@ -33,6 +36,14 @@ use crate::Error;
|
||||
/// is 4.9 GB, so each array holds an eighth to a quarter of it.
|
||||
pub const DETAIL_GROUPS: usize = 4;
|
||||
|
||||
/// The most bytes one texture may take on lavapipe.
|
||||
const MAX_TEXTURE_BYTES: u64 = 1 << 31;
|
||||
|
||||
/// How far into a region its cells fade into the sharp map where no fetched cell lies
|
||||
/// beyond, in degrees of latitude (about 11 km), so the change in sharpness does not
|
||||
/// show as a line.
|
||||
pub const REGION_FADE: f32 = 0.1;
|
||||
|
||||
/// Channels one radiance target holds.
|
||||
pub const MAX_CHANNELS: usize = 4;
|
||||
|
||||
@@ -96,6 +107,9 @@ pub struct BodyRadiometry {
|
||||
pub earth_textured: bool,
|
||||
/// The sharp tiled map's layout, once one is set.
|
||||
pub earth_detail: Option<DetailLayout>,
|
||||
/// The fine regions' layout and how many cells each of their texture arrays holds,
|
||||
/// once they are set.
|
||||
pub earth_regions: Option<(DetailLayout, u32)>,
|
||||
/// Whether Earth and the Moon shadow each other: lunar and solar eclipses.
|
||||
pub bodies_eclipsed: bool,
|
||||
/// Per channel, Earth's albedo from `Appearance::earth_albedo` (no texture).
|
||||
@@ -155,6 +169,7 @@ impl BodyRadiometry {
|
||||
earth_rgb: [basis(0), basis(1), basis(2)],
|
||||
atmosphere: false,
|
||||
earth_detail: None,
|
||||
earth_regions: None,
|
||||
earth_textured: false,
|
||||
bodies_eclipsed: false,
|
||||
earth_albedo: albedo(&appearance.earth_albedo),
|
||||
@@ -212,6 +227,12 @@ pub struct BodiesUniforms {
|
||||
projection: [f32; 4],
|
||||
// A Kannala-Brandt fisheye's coefficients.
|
||||
fisheye: [f32; 4],
|
||||
// The fine regions: x 1 if set, y mip levels, z a cell's interior over its stored
|
||||
// size, w its border over its stored size.
|
||||
earth_regions: [f32; 4],
|
||||
// x: cells per texture array, y, z: texels across and down the whole Earth at full
|
||||
// size, w: how far cells fade in at a region's edge, in degrees (REGION_FADE).
|
||||
earth_regions_size: [f32; 4],
|
||||
}
|
||||
|
||||
impl BodiesUniforms {
|
||||
@@ -318,6 +339,23 @@ impl BodiesUniforms {
|
||||
]
|
||||
}),
|
||||
clouds: [0.0; 4],
|
||||
earth_regions: radiometry.earth_regions.map_or([0.0; 4], |(l, _)| {
|
||||
let padded = l.padded(0) as f32;
|
||||
[
|
||||
1.0,
|
||||
l.levels as f32,
|
||||
l.tile as f32 / padded,
|
||||
l.gutter as f32 / padded,
|
||||
]
|
||||
}),
|
||||
earth_regions_size: radiometry.earth_regions.map_or([0.0; 4], |(l, per)| {
|
||||
[
|
||||
per as f32,
|
||||
(360 * l.tile) as f32,
|
||||
(180 * l.tile) as f32,
|
||||
REGION_FADE,
|
||||
]
|
||||
}),
|
||||
projection: projection[0],
|
||||
fisheye: projection[1],
|
||||
beta: scene.velocities.relative().map(|b| {
|
||||
@@ -380,6 +418,9 @@ pub struct BodiesPass {
|
||||
atmosphere: AtmosphereTextures,
|
||||
detail: [wgpu::TextureView; DETAIL_GROUPS],
|
||||
detail_sampler: wgpu::Sampler,
|
||||
/// The fine regions once they are set, which compile `regions.wgsl` into the
|
||||
/// pipelines.
|
||||
regions: Option<RegionTextures>,
|
||||
/// Rays through the atmosphere (`Grazing`), as `bodies.wgsl` declares them: the
|
||||
/// transmittance per entry, then the bending.
|
||||
grazing: wgpu::Buffer,
|
||||
@@ -424,7 +465,16 @@ impl BodiesPass {
|
||||
/// buffer.
|
||||
pub fn new(ctx: &GpuContext, samples: u32, area: bool) -> Self {
|
||||
let device = &ctx.device;
|
||||
let pipeline = compile(ctx, samples, area, false, Clouds::None, ANISOTROPY_L1, None);
|
||||
let pipeline = compile(
|
||||
ctx,
|
||||
samples,
|
||||
area,
|
||||
false,
|
||||
Clouds::None,
|
||||
ANISOTROPY_L1,
|
||||
None,
|
||||
false,
|
||||
);
|
||||
let uniforms = device.create_buffer(&wgpu::BufferDescriptor {
|
||||
label: Some("bodies uniforms"),
|
||||
size: std::mem::size_of::<BodiesUniforms>() as u64,
|
||||
@@ -485,6 +535,7 @@ impl BodiesPass {
|
||||
None,
|
||||
None,
|
||||
None,
|
||||
None,
|
||||
);
|
||||
let clouds = upload_clouds(ctx, &CloudMap::from_values(1, 1, vec![0]).expect("1 x 1"));
|
||||
// Unread unless the uniforms say the Moon has Earthshine.
|
||||
@@ -510,6 +561,7 @@ impl BodiesPass {
|
||||
atmosphere,
|
||||
detail,
|
||||
detail_sampler,
|
||||
regions: None,
|
||||
grazing,
|
||||
earth_shadow,
|
||||
clouds,
|
||||
@@ -557,6 +609,7 @@ impl BodiesPass {
|
||||
self.cloud_kind(),
|
||||
self.anisotropy,
|
||||
self.terrain.as_ref().map(|t| &t.shape),
|
||||
self.regions.is_some(),
|
||||
)
|
||||
}
|
||||
|
||||
@@ -760,98 +813,125 @@ impl BodiesPass {
|
||||
detail: &EarthDetail,
|
||||
) -> Result<(), Error> {
|
||||
let layout = detail.layout();
|
||||
let format = match detail.format() {
|
||||
DetailFormat::Rgba16Float => wgpu::TextureFormat::Rgba16Float,
|
||||
DetailFormat::Bc6h => {
|
||||
if !ctx
|
||||
.device
|
||||
.features()
|
||||
.contains(wgpu::Features::TEXTURE_COMPRESSION_BC)
|
||||
{
|
||||
return Err(Error::Unsupported(format!(
|
||||
"{} has no BC texture compression, which the sharp Earth map needs",
|
||||
ctx.info.name
|
||||
)));
|
||||
}
|
||||
wgpu::TextureFormat::Bc6hRgbUfloat
|
||||
}
|
||||
};
|
||||
let size = layout.padded(0);
|
||||
let per_group = layout.layers().div_ceil(DETAIL_GROUPS as u32);
|
||||
let mut views = Vec::new();
|
||||
for group in 0..DETAIL_GROUPS as u32 {
|
||||
let first = group * per_group;
|
||||
let count = per_group.min(layout.layers().saturating_sub(first)).max(1);
|
||||
let texture = ctx.device.create_texture(&wgpu::TextureDescriptor {
|
||||
label: Some("earth detail"),
|
||||
size: wgpu::Extent3d {
|
||||
width: size,
|
||||
height: size,
|
||||
depth_or_array_layers: count,
|
||||
},
|
||||
mip_level_count: layout.levels,
|
||||
sample_count: 1,
|
||||
dimension: wgpu::TextureDimension::D2,
|
||||
format,
|
||||
usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::COPY_DST,
|
||||
view_formats: &[],
|
||||
});
|
||||
for local in 0..count {
|
||||
let layer = first + local;
|
||||
if layer >= layout.layers() {
|
||||
break;
|
||||
}
|
||||
for (level, bytes) in detail.tile_levels(layer)?.iter().enumerate() {
|
||||
let s = layout.padded(level as u32);
|
||||
// A compressed level is written in whole blocks, past its edge if
|
||||
// need be.
|
||||
let (width, rows, row_bytes) = match detail.format() {
|
||||
DetailFormat::Rgba16Float => (s, s, 8 * s),
|
||||
DetailFormat::Bc6h => {
|
||||
(s.div_ceil(4) * 4, s.div_ceil(4), 16 * s.div_ceil(4))
|
||||
}
|
||||
};
|
||||
ctx.queue.write_texture(
|
||||
wgpu::TexelCopyTextureInfo {
|
||||
texture: &texture,
|
||||
mip_level: level as u32,
|
||||
origin: wgpu::Origin3d {
|
||||
x: 0,
|
||||
y: 0,
|
||||
z: local,
|
||||
},
|
||||
aspect: wgpu::TextureAspect::All,
|
||||
},
|
||||
bytes,
|
||||
wgpu::TexelCopyBufferLayout {
|
||||
offset: 0,
|
||||
bytes_per_row: Some(row_bytes),
|
||||
rows_per_image: Some(rows),
|
||||
},
|
||||
wgpu::Extent3d {
|
||||
width,
|
||||
height: width,
|
||||
depth_or_array_layers: 1,
|
||||
},
|
||||
);
|
||||
}
|
||||
// Hand each tile to the device before reading the next, so the staging
|
||||
// copies do not pile up: the whole map is 4.9 GB.
|
||||
ctx.queue.submit([]);
|
||||
ctx.device
|
||||
.poll(wgpu::PollType::wait_indefinitely())
|
||||
.map_err(|e| Error::Gpu(e.to_string()))?;
|
||||
}
|
||||
views.push(texture.create_view(&wgpu::TextureViewDescriptor {
|
||||
dimension: Some(wgpu::TextureViewDimension::D2Array),
|
||||
..Default::default()
|
||||
}));
|
||||
}
|
||||
let views = upload_tiles(
|
||||
ctx,
|
||||
layout,
|
||||
detail.format(),
|
||||
layout.layers(),
|
||||
per_group,
|
||||
|layer| detail.tile_levels(layer),
|
||||
)?;
|
||||
self.detail = views.try_into().expect("one view per group");
|
||||
self.rebind(ctx);
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Uploads the fine regions of Earth's surface, a cell at a time, and binds them for
|
||||
/// later frames, in place of any set before. Returns how many cells each texture
|
||||
/// array holds. A BC6H map needs a device with BC compression; more cells than
|
||||
/// DETAIL_GROUPS arrays of lavapipe's largest texture hold are refused.
|
||||
pub fn set_earth_regions(
|
||||
&mut self,
|
||||
ctx: &GpuContext,
|
||||
regions: &EarthRegions,
|
||||
) -> Result<u32, Error> {
|
||||
let maps = regions.maps();
|
||||
let cells = maps.len() as u32;
|
||||
let most = (MAX_TEXTURE_BYTES / maps[0].tile_bytes()) as u32;
|
||||
if cells > most * DETAIL_GROUPS as u32 {
|
||||
return Err(Error::Unsupported(format!(
|
||||
"{cells} cells of Earth's fine regions, {:.1} GB; at most {} fit",
|
||||
regions.bytes() as f64 / 1e9,
|
||||
most * DETAIL_GROUPS as u32
|
||||
)));
|
||||
}
|
||||
let per_group = cells.div_ceil(DETAIL_GROUPS as u32);
|
||||
let views = upload_tiles(
|
||||
ctx,
|
||||
regions.layout(),
|
||||
regions.format(),
|
||||
cells,
|
||||
per_group,
|
||||
|layer| maps[layer as usize].tile_levels(0),
|
||||
)?;
|
||||
let mut slots = vec![-1i32; 360 * 180];
|
||||
for (layer, map) in maps.iter().enumerate() {
|
||||
slots[map.cell().expect("a cell").index()] = layer as i32;
|
||||
}
|
||||
let slots = {
|
||||
use wgpu::util::DeviceExt;
|
||||
ctx.device
|
||||
.create_buffer_init(&wgpu::util::BufferInitDescriptor {
|
||||
label: Some("earth region slots"),
|
||||
contents: bytemuck::cast_slice(&slots),
|
||||
usage: wgpu::BufferUsages::STORAGE,
|
||||
})
|
||||
};
|
||||
// The coarse maps as RGBA f16, a layer per cell.
|
||||
let n = REGION_COARSE_PADDED;
|
||||
let coarse = ctx.device.create_texture(&wgpu::TextureDescriptor {
|
||||
label: Some("earth region coarse maps"),
|
||||
size: wgpu::Extent3d {
|
||||
width: n,
|
||||
height: n,
|
||||
depth_or_array_layers: cells,
|
||||
},
|
||||
mip_level_count: 1,
|
||||
sample_count: 1,
|
||||
dimension: wgpu::TextureDimension::D2,
|
||||
format: wgpu::TextureFormat::Rgba16Float,
|
||||
usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::COPY_DST,
|
||||
view_formats: &[],
|
||||
});
|
||||
for (layer, map) in maps.iter().enumerate() {
|
||||
let texels: Vec<half::f16> = map
|
||||
.coarse()
|
||||
.expect("regions have coarse maps")
|
||||
.iter()
|
||||
.flat_map(|rgb| [rgb[0], rgb[1], rgb[2], 1.0].map(half::f16::from_f32))
|
||||
.collect();
|
||||
ctx.queue.write_texture(
|
||||
wgpu::TexelCopyTextureInfo {
|
||||
texture: &coarse,
|
||||
mip_level: 0,
|
||||
origin: wgpu::Origin3d {
|
||||
x: 0,
|
||||
y: 0,
|
||||
z: layer as u32,
|
||||
},
|
||||
aspect: wgpu::TextureAspect::All,
|
||||
},
|
||||
bytemuck::cast_slice(&texels),
|
||||
wgpu::TexelCopyBufferLayout {
|
||||
offset: 0,
|
||||
bytes_per_row: Some(8 * n),
|
||||
rows_per_image: Some(n),
|
||||
},
|
||||
wgpu::Extent3d {
|
||||
width: n,
|
||||
height: n,
|
||||
depth_or_array_layers: 1,
|
||||
},
|
||||
);
|
||||
}
|
||||
let coarse = coarse.create_view(&wgpu::TextureViewDescriptor {
|
||||
dimension: Some(wgpu::TextureViewDimension::D2Array),
|
||||
..Default::default()
|
||||
});
|
||||
let compiled = self.regions.is_some();
|
||||
self.regions = Some(RegionTextures {
|
||||
cells: views.try_into().expect("one view per group"),
|
||||
slots,
|
||||
coarse,
|
||||
});
|
||||
if !compiled {
|
||||
self.recompile(ctx);
|
||||
}
|
||||
self.rebind(ctx);
|
||||
Ok(per_group)
|
||||
}
|
||||
|
||||
/// Uploads Earth's atmosphere's tables and binds them for later frames; the
|
||||
/// uniforms' atmosphere flag turns them on.
|
||||
pub fn set_atmosphere(&mut self, ctx: &GpuContext, tables: &AtmosphereTables) {
|
||||
@@ -920,6 +1000,7 @@ impl BodiesPass {
|
||||
&self.moon,
|
||||
&self.atmosphere,
|
||||
(&self.detail, &self.detail_sampler),
|
||||
self.regions.as_ref(),
|
||||
bending.then_some((&self.grazing, &self.earth_shadow)),
|
||||
self.clouded.then_some(&self.clouds),
|
||||
self.terrain.as_ref().map(|t| (&t.views, &t.slots)),
|
||||
@@ -1053,6 +1134,129 @@ fn placeholder_detail(ctx: &GpuContext) -> wgpu::TextureView {
|
||||
})
|
||||
}
|
||||
|
||||
/// The fine regions on the GPU: their cells, per 1-degree cell (`Cell::index`) its
|
||||
/// layer among them or -1, and their coarse maps, a layer each.
|
||||
#[derive(Debug)]
|
||||
struct RegionTextures {
|
||||
cells: [wgpu::TextureView; DETAIL_GROUPS],
|
||||
slots: wgpu::Buffer,
|
||||
coarse: wgpu::TextureView,
|
||||
}
|
||||
|
||||
/// The texture format of tiles in `format`; BC6H needs a device with BC compression.
|
||||
fn tile_format(ctx: &GpuContext, format: DetailFormat) -> Result<wgpu::TextureFormat, Error> {
|
||||
match format {
|
||||
DetailFormat::Rgba16Float => Ok(wgpu::TextureFormat::Rgba16Float),
|
||||
DetailFormat::Bc6h => {
|
||||
if !ctx
|
||||
.device
|
||||
.features()
|
||||
.contains(wgpu::Features::TEXTURE_COMPRESSION_BC)
|
||||
{
|
||||
return Err(Error::Unsupported(format!(
|
||||
"{} has no BC texture compression, which Earth's tiled maps need",
|
||||
ctx.info.name
|
||||
)));
|
||||
}
|
||||
Ok(wgpu::TextureFormat::Bc6hRgbUfloat)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Uploads `layers` tiles of `layout` in `format`, with their levels from
|
||||
/// `levels(layer)`, a tile at a time, into DETAIL_GROUPS texture arrays of
|
||||
/// `per_group` tiles in order (unused arrays hold one empty layer), and returns their
|
||||
/// views. Running out of device memory is an error, not a panic.
|
||||
fn upload_tiles(
|
||||
ctx: &GpuContext,
|
||||
layout: DetailLayout,
|
||||
format: DetailFormat,
|
||||
layers: u32,
|
||||
per_group: u32,
|
||||
levels: impl Fn(u32) -> Result<Vec<Vec<u8>>, Error>,
|
||||
) -> Result<Vec<wgpu::TextureView>, Error> {
|
||||
let texture_format = tile_format(ctx, format)?;
|
||||
let size = layout.padded(0);
|
||||
let mut views = Vec::new();
|
||||
for group in 0..DETAIL_GROUPS as u32 {
|
||||
let first = group * per_group;
|
||||
let count = per_group.min(layers.saturating_sub(first)).max(1);
|
||||
let scope = ctx.device.push_error_scope(wgpu::ErrorFilter::OutOfMemory);
|
||||
let texture = ctx.device.create_texture(&wgpu::TextureDescriptor {
|
||||
label: Some("earth tiles"),
|
||||
size: wgpu::Extent3d {
|
||||
width: size,
|
||||
height: size,
|
||||
depth_or_array_layers: count,
|
||||
},
|
||||
mip_level_count: layout.levels,
|
||||
sample_count: 1,
|
||||
dimension: wgpu::TextureDimension::D2,
|
||||
format: texture_format,
|
||||
usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::COPY_DST,
|
||||
view_formats: &[],
|
||||
});
|
||||
if let Some(e) = pollster::block_on(scope.pop()) {
|
||||
return Err(Error::Gpu(format!(
|
||||
"{layers} tiles of {:.0} MB need more memory than {} has: {e}",
|
||||
(0..layout.levels)
|
||||
.map(|l| format.bytes(layout.padded(l)))
|
||||
.sum::<u64>() as f64
|
||||
/ 1e6,
|
||||
ctx.info.name
|
||||
)));
|
||||
}
|
||||
for local in 0..count {
|
||||
let layer = first + local;
|
||||
if layer >= layers {
|
||||
break;
|
||||
}
|
||||
for (level, bytes) in levels(layer)?.iter().enumerate() {
|
||||
let s = layout.padded(level as u32);
|
||||
// A compressed level is written in whole blocks, past its edge if need be.
|
||||
let (width, rows, row_bytes) = match format {
|
||||
DetailFormat::Rgba16Float => (s, s, 8 * s),
|
||||
DetailFormat::Bc6h => (s.div_ceil(4) * 4, s.div_ceil(4), 16 * s.div_ceil(4)),
|
||||
};
|
||||
ctx.queue.write_texture(
|
||||
wgpu::TexelCopyTextureInfo {
|
||||
texture: &texture,
|
||||
mip_level: level as u32,
|
||||
origin: wgpu::Origin3d {
|
||||
x: 0,
|
||||
y: 0,
|
||||
z: local,
|
||||
},
|
||||
aspect: wgpu::TextureAspect::All,
|
||||
},
|
||||
bytes,
|
||||
wgpu::TexelCopyBufferLayout {
|
||||
offset: 0,
|
||||
bytes_per_row: Some(row_bytes),
|
||||
rows_per_image: Some(rows),
|
||||
},
|
||||
wgpu::Extent3d {
|
||||
width,
|
||||
height: width,
|
||||
depth_or_array_layers: 1,
|
||||
},
|
||||
);
|
||||
}
|
||||
// Hand each tile to the device before reading the next, so the staging
|
||||
// copies do not pile up: Blue Marble's sharp map is 4.9 GB.
|
||||
ctx.queue.submit([]);
|
||||
ctx.device
|
||||
.poll(wgpu::PollType::wait_indefinitely())
|
||||
.map_err(|e| Error::Gpu(e.to_string()))?;
|
||||
}
|
||||
views.push(texture.create_view(&wgpu::TextureViewDescriptor {
|
||||
dimension: Some(wgpu::TextureViewDimension::D2Array),
|
||||
..Default::default()
|
||||
}));
|
||||
}
|
||||
Ok(views)
|
||||
}
|
||||
|
||||
/// Uploads Earth's texture with its mip chain as linear RGBA f16.
|
||||
fn upload_earth(ctx: &GpuContext, texture: &EarthTexture) -> wgpu::TextureView {
|
||||
let levels = texture.mip_chain();
|
||||
@@ -1175,6 +1379,7 @@ fn compile(
|
||||
clouds: Clouds,
|
||||
anisotropy: u32,
|
||||
terrain: Option<&TerrainShape>,
|
||||
regions: bool,
|
||||
) -> wgpu::RenderPipeline {
|
||||
let device = &ctx.device;
|
||||
let (radiance, depth) = if samples > 1 {
|
||||
@@ -1216,7 +1421,13 @@ fn compile(
|
||||
crate::clouds::wgsl() + include_str!("../shaders/relief.wgsl")
|
||||
}
|
||||
}
|
||||
+ &terrain.map_or_else(|| SPHERE.to_owned(), TerrainShape::wgsl))
|
||||
+ &terrain.map_or_else(|| SPHERE.to_owned(), TerrainShape::wgsl)
|
||||
+ &if regions {
|
||||
format!("const REGION_COARSE: f32 = {REGION_COARSE}.0;\n")
|
||||
+ include_str!("../shaders/regions.wgsl")
|
||||
} else {
|
||||
COARSE.to_owned()
|
||||
})
|
||||
.into(),
|
||||
),
|
||||
});
|
||||
@@ -1452,6 +1663,13 @@ fn moon_band(c: vec3<f32>, d: vec3<f32>) -> vec2<f32> {
|
||||
}
|
||||
";
|
||||
|
||||
/// `read_earth` without the fine regions.
|
||||
const COARSE: &str = "
|
||||
fn read_earth(at: MapLookup) -> vec3<f32> {
|
||||
return read_coarse(at);
|
||||
}
|
||||
";
|
||||
|
||||
/// `through_air` and `earth_shadow` where no body is seen or lit through the
|
||||
/// atmosphere: every ray straight from the camera and whole, and Earth's shadow its hard
|
||||
/// edge's.
|
||||
@@ -1477,6 +1695,7 @@ fn bind(
|
||||
moon: &wgpu::TextureView,
|
||||
atmosphere: &AtmosphereTextures,
|
||||
detail: (&[wgpu::TextureView; DETAIL_GROUPS], &wgpu::Sampler),
|
||||
regions: Option<&RegionTextures>,
|
||||
air: Option<(&wgpu::Buffer, &wgpu::Buffer)>,
|
||||
clouds: Option<&wgpu::TextureView>,
|
||||
terrain: Option<(&[wgpu::TextureView; 3], &wgpu::Buffer)>,
|
||||
@@ -1520,6 +1739,25 @@ fn bind(
|
||||
binding: 10 + DETAIL_GROUPS as u32,
|
||||
resource: wgpu::BindingResource::Sampler(detail.1),
|
||||
}])
|
||||
.chain(regions.into_iter().flat_map(|r| {
|
||||
r.cells
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(i, view)| wgpu::BindGroupEntry {
|
||||
binding: 22 + i as u32,
|
||||
resource: wgpu::BindingResource::TextureView(view),
|
||||
})
|
||||
.chain([
|
||||
wgpu::BindGroupEntry {
|
||||
binding: 22 + DETAIL_GROUPS as u32,
|
||||
resource: r.slots.as_entire_binding(),
|
||||
},
|
||||
wgpu::BindGroupEntry {
|
||||
binding: 23 + DETAIL_GROUPS as u32,
|
||||
resource: wgpu::BindingResource::TextureView(&r.coarse),
|
||||
},
|
||||
])
|
||||
}))
|
||||
.chain(air.into_iter().flat_map(|(grazing, shadow)| {
|
||||
[
|
||||
wgpu::BindGroupEntry {
|
||||
|
||||
@@ -50,7 +50,7 @@ pub use camera::{
|
||||
Blooming, BrownConrady, CameraConfig, Channel, Detector, Fisheye, Intrinsics, Psf, PsfImage,
|
||||
};
|
||||
pub use clouds::CloudMap;
|
||||
pub use earth::{EarthDetail, EarthTexture};
|
||||
pub use earth::{EarthDetail, EarthRegions, EarthTexture};
|
||||
pub use error::Error;
|
||||
pub use fidelity::{CameraTier, Fidelity, LightingTier};
|
||||
pub use fixed_pattern::FixedPattern;
|
||||
|
||||
@@ -17,7 +17,7 @@ use std::sync::Arc;
|
||||
|
||||
use crate::atmosphere::{AtmosphereTables, EarthShadow};
|
||||
use crate::clouds::{CloudMap, CloudRelief};
|
||||
use crate::earth::{EarthDetail, EarthTexture};
|
||||
use crate::earth::{EarthDetail, EarthRegions, EarthTexture};
|
||||
use crate::eclipse::{
|
||||
moon_near_shadow, sun_visible_fraction, sun_visible_light, sun_visible_through_air,
|
||||
};
|
||||
@@ -60,6 +60,10 @@ pub struct RendererOptions {
|
||||
/// resolves ground finer than `earth_texture` does; optional. Uploading Blue
|
||||
/// Marble's takes about 4.9 GB of device memory.
|
||||
pub earth_detail: Option<Arc<EarthDetail>>,
|
||||
/// Fine maps of Earth's surface over chosen 1-degree cells, drawn from lighting tier
|
||||
/// L1 where a view resolves ground finer than the sharp map does; optional. Each of
|
||||
/// ESA WorldCover's cells takes 194 MB of device memory.
|
||||
pub earth_regions: Option<Arc<EarthRegions>>,
|
||||
/// The Moon's surface, needed from lighting tier L1 up; ignored below it.
|
||||
pub moon_texture: Option<Arc<MoonTexture>>,
|
||||
/// The Moon's terrain, drawn from lighting tier L2 in place of a smooth sphere;
|
||||
@@ -503,6 +507,9 @@ impl Renderer {
|
||||
if let Some(detail) = &options.earth_detail {
|
||||
renderer.set_earth_detail(detail)?;
|
||||
}
|
||||
if let Some(regions) = &options.earth_regions {
|
||||
renderer.set_earth_regions(regions)?;
|
||||
}
|
||||
if let Some(texture) = &options.moon_texture {
|
||||
renderer.set_moon_texture(texture);
|
||||
}
|
||||
@@ -562,6 +569,25 @@ impl Renderer {
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Draws Earth's surface from the fine maps `regions` from now on where a view
|
||||
/// resolves ground finer than the sharp map does, replacing any given in
|
||||
/// [`RendererOptions::earth_regions`]. Uploads them a cell at a time from their
|
||||
/// files. A BC6H map needs a device with BC compression.
|
||||
///
|
||||
/// Below lighting tier L1 there is no Earth texture, and this does nothing.
|
||||
pub fn set_earth_regions(&mut self, regions: &EarthRegions) -> Result<(), Error> {
|
||||
if self.fidelity.lighting < LightingTier::L1 {
|
||||
log::warn!(
|
||||
"lighting {} draws no Earth texture; ignoring the fine regions",
|
||||
self.fidelity.lighting
|
||||
);
|
||||
return Ok(());
|
||||
}
|
||||
let per_group = self.bodies.set_earth_regions(&self.ctx, regions)?;
|
||||
self.body_radiometry.earth_regions = Some((regions.layout(), per_group));
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Draws the Moon with `terrain`'s shape from now on, in place of a smooth sphere:
|
||||
/// its limb, its slopes and its shadows. Uploads it a tile at a time from its file,
|
||||
/// the full-size level later, for the tiles close views need.
|
||||
|
||||
@@ -12,7 +12,8 @@
|
||||
// targets, MeshRadiance and MeshDepth (multisampled when SAMPLES > 1), and
|
||||
// SAMPLE_POSITIONS, where in a pixel each sample lies, and ANISOTROPY, the most reads
|
||||
// `read_map` and `read_detail` take along a footprint's long axis, and adds
|
||||
// `through_air`, `earth_shadow` and `through_clouds` (bodies.rs), and the Moon's shape:
|
||||
// `through_air`, `earth_shadow` and `through_clouds` (bodies.rs), `read_earth`
|
||||
// (regions.wgsl, or the coarser maps alone), and the Moon's shape:
|
||||
// `moon_distance`, `moon_surface`, `moon_shade` and `moon_band`, for a sphere or for
|
||||
// its terrain (terrain.wgsl).
|
||||
|
||||
@@ -78,6 +79,12 @@ struct Uniforms {
|
||||
projection: vec4<f32>,
|
||||
// A Kannala-Brandt fisheye's coefficients.
|
||||
fisheye: vec4<f32>,
|
||||
// The fine regions (earth.rs, `EarthRegions`): x 1 if set, y mip levels, z a cell's
|
||||
// interior over its stored size, w its border over its stored size.
|
||||
earth_regions: vec4<f32>,
|
||||
// x: cells per texture array, y, z: texels across and down the whole Earth at full
|
||||
// size, w: how far cells fade in at a region's edge, in degrees.
|
||||
earth_regions_size: vec4<f32>,
|
||||
}
|
||||
|
||||
@group(0) @binding(0) var<uniform> u: Uniforms;
|
||||
@@ -98,6 +105,7 @@ const DETAIL_GROUPS: u32 = 4u;
|
||||
@group(0) @binding(13) var earth_detail_3: texture_2d_array<f32>;
|
||||
@group(0) @binding(14) var detail_sampler: sampler;
|
||||
|
||||
|
||||
// The mesh pass's targets, SAMPLES per pixel.
|
||||
@group(1) @binding(0) var mesh_radiance: MeshRadiance;
|
||||
@group(1) @binding(1) var mesh_depth: MeshDepth;
|
||||
@@ -166,6 +174,9 @@ struct MapLookup {
|
||||
dx: vec2<f32>,
|
||||
dy: vec2<f32>,
|
||||
u_stretch: f32,
|
||||
// The body-fixed unit normal, for the fine regions (regions.wgsl), which find
|
||||
// their place from it more finely than from `uv`.
|
||||
normal: vec3<f32>,
|
||||
}
|
||||
|
||||
// East and north at body-fixed unit normal `n`, and cos(latitude); at a pole, east is
|
||||
@@ -197,6 +208,7 @@ fn map_lookup(n: vec3<f32>) -> MapLookup {
|
||||
ground_gradient(dpdx(n), frame),
|
||||
ground_gradient(dpdy(n), frame),
|
||||
1.0 / frame.cos_lat,
|
||||
n,
|
||||
);
|
||||
}
|
||||
|
||||
@@ -222,6 +234,7 @@ fn map_lookup_at(n: vec3<f32>, d: vec3<f32>, t: f32, spacing: f32, radius: f32)
|
||||
ground_gradient(ground * (along / radius), frame),
|
||||
ground_gradient(side * (across / radius), frame),
|
||||
1.0 / frame.cos_lat,
|
||||
n,
|
||||
);
|
||||
}
|
||||
|
||||
@@ -313,10 +326,11 @@ fn read_tile(st: vec2<f32>, group: u32, local: i32, dx: vec2<f32>, dy: vec2<f32>
|
||||
}
|
||||
}
|
||||
|
||||
// Earth's surface reflectance at `at`. With the sharp map, it is read where a pixel
|
||||
// covers less ground than its coarsest level's texels resolve, the global map beyond
|
||||
// that, and the two are blended across half a level so the change does not show.
|
||||
fn read_earth(at: MapLookup) -> vec3<f32> {
|
||||
// Earth's surface reflectance at `at` without the fine regions. With the sharp map, it
|
||||
// is read where a pixel covers less ground than its coarsest level's texels resolve,
|
||||
// the global map beyond that, and the two are blended across half a level so the
|
||||
// change does not show.
|
||||
fn read_coarse(at: MapLookup) -> vec3<f32> {
|
||||
if u.earth_detail.x < 0.5 {
|
||||
return read_map(earth_texture, at).rgb;
|
||||
}
|
||||
|
||||
@@ -0,0 +1,189 @@
|
||||
// Earth's fine regions (earth.rs, `EarthRegions`; bodies.rs compiles this in only when
|
||||
// they are set, and otherwise has `read_earth` read `read_coarse` alone): 1 x 1 degree
|
||||
// cells of a finer map, read where a view resolves them and blended into the map
|
||||
// below.
|
||||
|
||||
// The fine regions' cells, one per layer, split in order among DETAIL_GROUPS arrays
|
||||
// as the sharp map's tiles are, read with its sampler; and for each 1 x 1 degree cell of
|
||||
// Earth, in rows from the north-west, its layer or -1.
|
||||
@group(0) @binding(22) var earth_region_0: texture_2d_array<f32>;
|
||||
@group(0) @binding(23) var earth_region_1: texture_2d_array<f32>;
|
||||
@group(0) @binding(24) var earth_region_2: texture_2d_array<f32>;
|
||||
@group(0) @binding(25) var earth_region_3: texture_2d_array<f32>;
|
||||
@group(0) @binding(26) var<storage, read> region_slots: array<i32>;
|
||||
// Each cell's mean colors in REGION_COARSE x REGION_COARSE blocks, with a border of one
|
||||
// block from its neighbors (earth.rs, `EarthDetail::coarse`), a layer per cell.
|
||||
@group(0) @binding(27) var region_coarse: texture_2d_array<f32>;
|
||||
|
||||
// Longitude and latitude in degrees at body-fixed unit normal `n`: held to 0.05 of the
|
||||
// regions' 9 m texels at 33 degrees east and 0.18 at 180, against about 0.26 for the
|
||||
// map coordinates `uv` the other maps use.
|
||||
fn lon_lat_degrees(n: vec3<f32>) -> vec2<f32> {
|
||||
return vec2<f32>(atan2(n.y, n.x), asin(clamp(n.z, -1.0, 1.0))) * (180.0 / PI);
|
||||
}
|
||||
|
||||
// One read of cell `local` of group `group` of the fine regions at `st`, with
|
||||
// gradients `dx` and `dy`.
|
||||
fn read_region_tile(st: vec2<f32>, group: u32, local: i32, dx: vec2<f32>, dy: vec2<f32>) -> vec3<f32> {
|
||||
switch group {
|
||||
case 0u: {
|
||||
return textureSampleGrad(earth_region_0, detail_sampler, st, local, dx, dy).rgb;
|
||||
}
|
||||
case 1u: {
|
||||
return textureSampleGrad(earth_region_1, detail_sampler, st, local, dx, dy).rgb;
|
||||
}
|
||||
case 2u: {
|
||||
return textureSampleGrad(earth_region_2, detail_sampler, st, local, dx, dy).rgb;
|
||||
}
|
||||
default: {
|
||||
return textureSampleGrad(earth_region_3, detail_sampler, st, local, dx, dy).rgb;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// The fine regions' layer for the cell `column` east of -180 degrees and `row` south of
|
||||
// +90, wrapping in longitude; -1 where none was fetched, and past the poles.
|
||||
fn region_slot(column: i32, row: i32) -> i32 {
|
||||
if row < 0 || row > 179 {
|
||||
return -1;
|
||||
}
|
||||
return region_slots[u32(row) * 360u + u32((column + 360) % 360)];
|
||||
}
|
||||
|
||||
// The cell at longitude and latitude `degrees`, as a column east of -180 degrees and a
|
||||
// row south of +90, and the place in it, from its north-west corner, in cells.
|
||||
struct InCell {
|
||||
cell: vec2<i32>,
|
||||
local: vec2<f32>,
|
||||
}
|
||||
|
||||
fn in_cell(degrees: vec2<f32>) -> InCell {
|
||||
let lat = clamp(degrees.y, -89.99999, 89.99999);
|
||||
let corner = floor(vec2<f32>(degrees.x, lat));
|
||||
let cell = vec2<i32>((i32(corner.x) + 540) % 360, 89 - i32(corner.y));
|
||||
return InCell(cell, vec2<f32>(degrees.x - corner.x, 1.0 - (lat - corner.y)));
|
||||
}
|
||||
|
||||
// How much of the fine regions to draw at longitude and latitude `degrees`, by place
|
||||
// alone: 1 inside them, fading to 0 within REGION_FADE (u.earth_regions_size.w) degrees
|
||||
// of ground of an edge with no fetched cell beyond, and 0 outside. `cos_lat` shortens
|
||||
// degrees of longitude to ground.
|
||||
fn region_cover(degrees: vec2<f32>, cos_lat: f32) -> f32 {
|
||||
let at = in_cell(degrees);
|
||||
let cell = at.cell;
|
||||
if region_slot(cell.x, cell.y) < 0 {
|
||||
return 0.0;
|
||||
}
|
||||
// Ground distance in degrees to each side of the cell: west, east, north, south.
|
||||
let west = at.local.x * cos_lat;
|
||||
let east = (1.0 - at.local.x) * cos_lat;
|
||||
let north = at.local.y;
|
||||
let south = 1.0 - at.local.y;
|
||||
var nearest = 1e9;
|
||||
for (var dy = -1; dy <= 1; dy++) {
|
||||
for (var dx = -1; dx <= 1; dx++) {
|
||||
if (dx == 0 && dy == 0) || region_slot(cell.x + dx, cell.y + dy) >= 0 {
|
||||
continue;
|
||||
}
|
||||
let across = select(select(0.0, east, dx > 0), west, dx < 0);
|
||||
let down = select(select(0.0, south, dy > 0), north, dy < 0);
|
||||
nearest = min(nearest, length(vec2<f32>(across, down)));
|
||||
}
|
||||
}
|
||||
return smoothstep(0.0, u.earth_regions_size.w, nearest);
|
||||
}
|
||||
|
||||
// Earth's surface at `at` from the fine regions, filtered as `read_detail` does, with
|
||||
// footprint `f` in their texels: each read from the cell it falls in, or where that
|
||||
// cell was not fetched, from the cell at `at` itself, read past its edge (its border
|
||||
// and then the clamped sampler repeat its edge).
|
||||
fn read_region(at: MapLookup, degrees: vec2<f32>, f: Footprint) -> vec3<f32> {
|
||||
let cells = vec2<f32>(360.0, 180.0);
|
||||
let scale = cells * u.earth_regions.z;
|
||||
let step = f.major / f32(f.taps) * scale;
|
||||
let minor = f.minor * scale;
|
||||
let per_group = u32(u.earth_regions_size.x);
|
||||
let home = in_cell(degrees).cell;
|
||||
// The long axis in degrees east and north.
|
||||
let major = f.major * vec2<f32>(360.0 * at.u_stretch, -180.0);
|
||||
var sum = vec3<f32>(0.0);
|
||||
for (var i = 0u; i < f.taps; i++) {
|
||||
let offset = (f32(i) + 0.5) / f32(f.taps) - 0.5;
|
||||
var p = in_cell(degrees + offset * major);
|
||||
var slot = region_slot(p.cell.x, p.cell.y);
|
||||
if slot < 0 {
|
||||
// Toward the cell at `at`, the short way round in longitude.
|
||||
var dx = f32(p.cell.x - home.x);
|
||||
dx -= 360.0 * round(dx / 360.0);
|
||||
p.local += vec2<f32>(dx, f32(p.cell.y - home.y));
|
||||
slot = region_slot(home.x, home.y);
|
||||
}
|
||||
let st = u.earth_regions.w + p.local * u.earth_regions.z;
|
||||
let layer = u32(max(slot, 0));
|
||||
sum += read_region_tile(st, layer / per_group, i32(layer % per_group), step, minor);
|
||||
}
|
||||
return sum / f32(f.taps);
|
||||
}
|
||||
|
||||
// The month's map's colors and the cell's own coarse map at `at`, both over blocks of
|
||||
// REGION_COARSE a cell (4.4 km).
|
||||
struct MonthMatch {
|
||||
month: vec3<f32>,
|
||||
own: vec3<f32>,
|
||||
}
|
||||
|
||||
const REGION_LUMA = vec3<f32>(0.2126, 0.7152, 0.0722);
|
||||
|
||||
fn month_match(at: MapLookup, degrees: vec2<f32>) -> MonthMatch {
|
||||
let p = in_cell(degrees);
|
||||
let layer = max(region_slot(p.cell.x, p.cell.y), 0);
|
||||
let st = (1.0 + p.local * REGION_COARSE) / (REGION_COARSE + 2.0);
|
||||
let own = textureSampleLevel(region_coarse, detail_sampler, st, layer, 0.0).rgb;
|
||||
let block = vec2<f32>(1.0 / (360.0 * REGION_COARSE), 1.0 / (180.0 * REGION_COARSE));
|
||||
let wide = MapLookup(at.uv, vec2<f32>(block.x, 0.0), vec2<f32>(0.0, block.y), at.u_stretch, at.normal);
|
||||
return MonthMatch(read_coarse(wide), own);
|
||||
}
|
||||
|
||||
// The fine regions' texel `rgb` brought to the month's map's colors by `m`: scaled per
|
||||
// channel by the month's map over the cell's coarse map. The cells are one year's
|
||||
// median, greener, grayer and bluer than a dry July (around Gibraltar, July has 1.2
|
||||
// times their red and half their blue, relative to green); brought to the month's map
|
||||
// at 4.4 km, they keep their finer detail and the seasons, and their edges match the
|
||||
// map around them. The scale is held to 0.4 to 2.5 per channel.
|
||||
//
|
||||
// A texel under a sixth as bright as the month's map over its block is water the map
|
||||
// does not have (Lake Timsah, in a 2.45 km texel of desert, or a coast it draws
|
||||
// elsewhere), and keeps its own colors, fully from a third: scaled like the land
|
||||
// around it, the lake's shallows turned orange. Open sea is as dark in both, and
|
||||
// matched.
|
||||
fn month_matched(rgb: vec3<f32>, m: MonthMatch) -> vec3<f32> {
|
||||
let ratio = clamp(m.month / max(m.own, vec3<f32>(1e-4)), vec3<f32>(0.4), vec3<f32>(2.5));
|
||||
let like = smoothstep(1.0 / 6.0, 1.0 / 3.0, dot(rgb, REGION_LUMA) / max(dot(m.month, REGION_LUMA), 1e-4));
|
||||
return rgb * mix(vec3<f32>(1.0), ratio, like);
|
||||
}
|
||||
|
||||
// Earth's surface reflectance at `at`. The fine regions are read where they were
|
||||
// fetched and a pixel covers less ground than their coarsest level's texels resolve,
|
||||
// blended into what lies below across half a level and toward a region's open edges
|
||||
// (`region_cover`), with their colors brought to the month's map's (`month_matched`).
|
||||
fn read_earth(at: MapLookup) -> vec3<f32> {
|
||||
let texels = u.earth_regions_size.yz;
|
||||
let f = footprint(at.dx, at.dy, texels);
|
||||
let read = max(length(f.minor * texels), length(f.major * texels) / f32(f.taps));
|
||||
let level = log2(max(read, 1e-6));
|
||||
let coarsest = u.earth_regions.y - 1.0;
|
||||
var fine = 1.0 - smoothstep(coarsest - 0.5, coarsest, level);
|
||||
var degrees = vec2<f32>(0.0);
|
||||
if fine > 0.0 {
|
||||
degrees = lon_lat_degrees(at.normal);
|
||||
fine *= region_cover(degrees, 1.0 / at.u_stretch);
|
||||
}
|
||||
var rgb = vec3<f32>(0.0);
|
||||
if fine > 0.0 {
|
||||
rgb = month_matched(read_region(at, degrees, f), month_match(at, degrees)) * fine;
|
||||
}
|
||||
if fine < 1.0 {
|
||||
rgb += read_coarse(at) * (1.0 - fine);
|
||||
}
|
||||
return rgb;
|
||||
}
|
||||
@@ -0,0 +1,276 @@
|
||||
//! Earth's fine regions (`EarthRegions`): cells read as one map across their edges, at
|
||||
//! magnifications and at coarser levels; read only where a view resolves them and only
|
||||
//! over cells that are there; brought to the month's map's colors at 4.4 km; and faded
|
||||
//! into the map below at a region's open edges. On the GPU and on the CPU (lavapipe).
|
||||
|
||||
use crate::common;
|
||||
|
||||
use resolved::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS};
|
||||
use resolved::earth::{Cell, DetailFormat, DetailLayout, EarthDetail, EarthRegions};
|
||||
use resolved::gpu::AdapterChoice;
|
||||
use resolved::nalgebra::{UnitQuaternion, Vector3};
|
||||
use resolved::radiometry::AU;
|
||||
use resolved::{
|
||||
CameraConfig, EarthTexture, Fidelity, LightingTier, Pose, Renderer, RendererOptions, Scene,
|
||||
};
|
||||
use std::sync::Arc;
|
||||
|
||||
const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu];
|
||||
|
||||
/// Cells of 1000 texels, 111 m at the equator, with levels to 445 m, and coarse maps of
|
||||
/// 40-texel blocks.
|
||||
const LAYOUT: DetailLayout = DetailLayout {
|
||||
tiles: [1, 1],
|
||||
tile: 1000,
|
||||
gutter: 4,
|
||||
levels: 3,
|
||||
};
|
||||
|
||||
/// Linear reflectance at `x` texels east of 0 degrees and `y` south of the equator:
|
||||
/// dark like water and bright like desert in bands 300 texels wide, with smooth detail
|
||||
/// at scales from 3 to 60 texels.
|
||||
fn pattern(x: f64, y: f64) -> [f32; 3] {
|
||||
let detail = (x / 3.1).sin() * (y / 4.3).cos() * 0.3
|
||||
+ (x / 11.0 + y / 7.0).sin() * 0.4
|
||||
+ (y / 60.0).cos() * (x / 45.0).sin() * 0.3;
|
||||
let bright = (x / 300.0).floor().rem_euclid(2.0) == 1.0;
|
||||
let level = |base: f64, swing: f64| (base + swing * detail).max(0.005) as f32;
|
||||
if bright {
|
||||
[level(0.4, 0.1), level(0.3, 0.08), level(0.2, 0.05)]
|
||||
} else {
|
||||
[level(0.02, 0.01), level(0.04, 0.015), level(0.06, 0.02)]
|
||||
}
|
||||
}
|
||||
|
||||
/// The cells with corners from `south`, `west` to `north`, `east`, each a map of
|
||||
/// `texel(x, y)` at its texels' centers (texels east of 0 degrees and south of the
|
||||
/// equator), in f16, with its neighbors' texels in its border and its coarse map.
|
||||
fn cells(
|
||||
[south, west, north, east]: [i32; 4],
|
||||
texel: impl Fn(f64, f64) -> [f32; 3] + Sync + Copy,
|
||||
) -> EarthRegions {
|
||||
let n = f64::from(LAYOUT.tile);
|
||||
let mut maps = Vec::new();
|
||||
for s in south..north {
|
||||
for w in west..east {
|
||||
let (x0, y0) = (f64::from(w) * n, -f64::from(s + 1) * n);
|
||||
let at = |col: i64, row: i64| texel(x0 + col as f64 + 0.5, y0 + row as f64 + 0.5);
|
||||
let levels = EarthDetail::build_levels(&LAYOUT, at);
|
||||
let coarse = EarthDetail::coarse_map(&LAYOUT, at);
|
||||
let bytes: Vec<u8> = levels
|
||||
.iter()
|
||||
.flat_map(|l| bytemuck::cast_slice::<_, u8>(l).to_vec())
|
||||
.collect();
|
||||
let map = EarthDetail::from_bytes(LAYOUT, DetailFormat::Rgba16Float, bytes).unwrap();
|
||||
maps.push(map.at_cell(Cell::new(s, w).unwrap(), 0xff, coarse).unwrap());
|
||||
}
|
||||
}
|
||||
EarthRegions::new(maps).unwrap()
|
||||
}
|
||||
|
||||
/// A camera of the example's optics, 320 x 240 pixels.
|
||||
fn camera() -> CameraConfig {
|
||||
CameraConfig {
|
||||
width: 320,
|
||||
height: 240,
|
||||
..CameraConfig::example()
|
||||
}
|
||||
}
|
||||
|
||||
/// Lighting L1 over a uniform global map of linear reflectance `ground` (sRGB-encoded,
|
||||
/// so close to it), with `regions` if given, and the bare surface (no atmosphere).
|
||||
fn renderer(adapter: &AdapterChoice, ground: u8, regions: Option<EarthRegions>) -> Renderer {
|
||||
let options = RendererOptions {
|
||||
adapter: adapter.clone(),
|
||||
fidelity: Fidelity {
|
||||
lighting: LightingTier::L1,
|
||||
..Default::default()
|
||||
},
|
||||
earth_texture: Some(Arc::new(
|
||||
EarthTexture::from_rgba8(64, 32, [ground, ground, ground, 255].repeat(64 * 32))
|
||||
.unwrap(),
|
||||
)),
|
||||
earth_regions: regions.map(Arc::new),
|
||||
moon_texture: Some(common::uniform_moon()),
|
||||
appearance: Appearance {
|
||||
earth_atmosphere: false,
|
||||
..Default::default()
|
||||
},
|
||||
read_radiance: true,
|
||||
..Default::default()
|
||||
};
|
||||
Renderer::new(camera(), options).unwrap()
|
||||
}
|
||||
|
||||
/// Looking straight down from `altitude` over geodetic `lat`, `lon` in degrees, north
|
||||
/// up, the Sun behind the camera: so views over different places are the same view.
|
||||
fn view(lat: f64, lon: f64, altitude: f64) -> (Scene, Pose) {
|
||||
let (la, lo) = (lat.to_radians(), lon.to_radians());
|
||||
let normal = Vector3::new(la.cos() * lo.cos(), la.cos() * lo.sin(), la.sin());
|
||||
let (a, b) = (EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS);
|
||||
let e2 = 1.0 - (b * b) / (a * a);
|
||||
let n = a / (1.0 - e2 * la.sin().powi(2)).sqrt();
|
||||
let ground = Vector3::new(
|
||||
n * la.cos() * lo.cos(),
|
||||
n * la.cos() * lo.sin(),
|
||||
n * (1.0 - e2) * la.sin(),
|
||||
);
|
||||
let scene = Scene {
|
||||
earth_attitude: UnitQuaternion::identity(),
|
||||
sun: normal * AU,
|
||||
moon: Vector3::new(0.0, 0.0, -3.844e8),
|
||||
moon_attitude: UnitQuaternion::identity(),
|
||||
targets: vec![],
|
||||
velocities: Default::default(),
|
||||
};
|
||||
let camera = Pose::looking_at(ground + normal * altitude, ground, Vector3::z());
|
||||
(scene, camera.unwrap())
|
||||
}
|
||||
|
||||
fn radiance(renderer: &mut Renderer, view: &(Scene, Pose)) -> Vec<[f32; 4]> {
|
||||
let frame = renderer.submit(&view.0, &view.1).unwrap().wait().unwrap();
|
||||
frame.radiance.unwrap()
|
||||
}
|
||||
|
||||
/// The largest and the mean difference between two images, over the brightest value
|
||||
/// in them.
|
||||
fn differences(a: &[[f32; 4]], b: &[[f32; 4]]) -> (f64, f64) {
|
||||
let peak = a
|
||||
.iter()
|
||||
.flat_map(|p| &p[..3])
|
||||
.fold(0.0f32, |m, &v| m.max(v));
|
||||
let all: Vec<f64> = a
|
||||
.iter()
|
||||
.zip(b)
|
||||
.flat_map(|(p, q)| (0..3).map(move |c| f64::from((p[c] - q[c]).abs() / peak)))
|
||||
.collect();
|
||||
let worst = all.iter().fold(0.0, |m: f64, &v| m.max(v));
|
||||
(worst, all.iter().sum::<f64>() / all.len() as f64)
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn cells_read_as_one_map() {
|
||||
// The same ground seen twice: once where four cells meet, at 0, 0, and once 0.52
|
||||
// degrees north-east, inside one cell, with the pattern moved there too. 520 texels
|
||||
// is a whole number of texels at every level and of the coarse maps' blocks, so the
|
||||
// two read the same texels, and a border that did not hold the neighbor's texels
|
||||
// would show along the cells' edges, by up to the contrast between neighbors. From
|
||||
// 500 km a pixel covers 55 m, half a texel (magnified); from 2500 km, 270 m, between
|
||||
// the second and third levels. The pattern's bands change from dark to bright on the
|
||||
// meridian, so a border repeating a cell's own edge there is off by up to 24 % of the
|
||||
// peak (measured). The two views differ by up to 0.38 % (500 km) and 1.14 % (2500 km)
|
||||
// on the RTX 3060 and 0.98 % and 1.47 % on lavapipe, at the bands' edges: texture
|
||||
// filtering blends neighbors in steps of 1/256 of the contrast between the bands,
|
||||
// and the two views' texel coordinates round differently. Brought to the map's gray,
|
||||
// the dark bands are scaled up 5 times and the bright ones down; the mean
|
||||
// difference is up to 9e-5 on the RTX 3060 and 1.3e-4 on lavapipe.
|
||||
let _devices = common::device_lock();
|
||||
let corner = cells([-1, -1, 1, 1], pattern);
|
||||
let shift = 520.0;
|
||||
let middle = cells([-1, -1, 2, 2], move |x, y| pattern(x - shift, y + shift));
|
||||
for adapter in &ADAPTERS {
|
||||
let mut four = renderer(adapter, 128, Some(corner.clone()));
|
||||
let mut one = renderer(adapter, 128, Some(middle.clone()));
|
||||
let name = four.adapter_info().name.clone();
|
||||
for altitude in [5e5, 2.5e6] {
|
||||
let a = radiance(&mut four, &view(0.0, 0.0, altitude));
|
||||
let b = radiance(&mut one, &view(0.52, 0.52, altitude));
|
||||
let (worst, mean) = differences(&a, &b);
|
||||
assert!(worst < 0.02, "{name} from {altitude} m: {worst}");
|
||||
assert!(mean < 2e-4, "{name} from {altitude} m: mean {mean}");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn regions_are_read_only_where_they_resolve_the_ground() {
|
||||
// Far views, past the cells' coarsest level, and views of ground with no cell, draw
|
||||
// what they would without the regions, to the bit.
|
||||
let _devices = common::device_lock();
|
||||
let region = cells([0, 0, 1, 1], pattern);
|
||||
for adapter in &ADAPTERS {
|
||||
let mut plain = renderer(adapter, 128, None);
|
||||
let mut fine = renderer(adapter, 128, Some(region.clone()));
|
||||
let name = plain.adapter_info().name.clone();
|
||||
// 2e7 m up a pixel covers 2.2 km, 20 texels; 2 degrees away from the cell, none.
|
||||
for (lat, lon, altitude) in [(0.5, 0.5, 2e7), (0.5, 3.0, 5e5)] {
|
||||
let v = view(lat, lon, altitude);
|
||||
let (a, b) = (radiance(&mut plain, &v), radiance(&mut fine, &v));
|
||||
let differing = a.iter().zip(&b).filter(|(p, q)| p != q).count();
|
||||
assert_eq!(differing, 0, "{name} at {lat}, {lon} from {altitude} m");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn cells_take_the_months_colors() {
|
||||
// A uniform cell twice as bright as the month's map below it draws the map's
|
||||
// brightness: its colors are brought to the map's at the coarse maps' scale.
|
||||
let _devices = common::device_lock();
|
||||
let bright = cells([0, 0, 1, 1], |_, _| [0.2; 3]);
|
||||
for adapter in &ADAPTERS {
|
||||
let mut plain = renderer(adapter, 0x59, None);
|
||||
let mut fine = renderer(adapter, 0x59, Some(bright.clone()));
|
||||
let name = plain.adapter_info().name.clone();
|
||||
let v = view(0.5, 0.5, 1.5e6);
|
||||
let (a, b) = (radiance(&mut plain, &v), radiance(&mut fine, &v));
|
||||
let (worst, _) = differences(&a, &b);
|
||||
assert!(worst < 0.01, "{name}: {worst}");
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn open_edges_fade_into_the_map_below() {
|
||||
// A cell of stripes 1.1 km wide, 0.1 and 0.3, over a uniform map of 0.1, seen
|
||||
// across its east edge from 1500 km (165 m a pixel). Brought to the map's colors,
|
||||
// the stripes swing about the map's brightness inside the cell; within REGION_FADE
|
||||
// (0.1 degrees, 11 km) of the edge they fade out, and past it the map is drawn as
|
||||
// it is. The brightness averaged over the stripes stays the map's throughout.
|
||||
let _devices = common::device_lock();
|
||||
let stripes = cells([0, 0, 1, 1], |x, _| {
|
||||
[if (x / 10.0).floor().rem_euclid(2.0) == 0.0 {
|
||||
0.1
|
||||
} else {
|
||||
0.3
|
||||
}; 3]
|
||||
});
|
||||
for adapter in &ADAPTERS {
|
||||
let mut plain = renderer(adapter, 0x59, None);
|
||||
let mut fine = renderer(adapter, 0x59, Some(stripes.clone()));
|
||||
let name = plain.adapter_info().name.clone();
|
||||
let v = view(0.5, 1.0, 1.5e6);
|
||||
let (a, b) = (radiance(&mut plain, &v), radiance(&mut fine, &v));
|
||||
let width = camera().width as usize;
|
||||
// The middle row, west to east: the camera looks down with north up.
|
||||
let row = 120 * width;
|
||||
let ratio: Vec<f64> = (0..width)
|
||||
.map(|x| f64::from(b[row + x][1]) / f64::from(a[row + x][1]))
|
||||
.collect();
|
||||
// The edge at the middle column, the fade's start 67 px west.
|
||||
let (edge, fade) = (width / 2, 111_320.0f64 * 0.1 / 165.0);
|
||||
let start = edge - fade.ceil() as usize - 2;
|
||||
// Two stripes' periods, 27 px.
|
||||
let window = |from: usize| {
|
||||
let w = &ratio[from..from + 27];
|
||||
let mean = w.iter().sum::<f64>() / w.len() as f64;
|
||||
let swing =
|
||||
w.iter().fold(0.0f64, |m, &r| m.max(r)) - w.iter().fold(9.0f64, |m, &r| m.min(r));
|
||||
(mean, swing)
|
||||
};
|
||||
let (_, inside) = window(start - 30);
|
||||
assert!(inside > 0.5, "{name}: the stripes swing by {inside} inside");
|
||||
let (_, near) = window(edge - 27);
|
||||
// The last 27 px are the last 40 % of the fade, where the cell weighs at most 0.35.
|
||||
assert!(
|
||||
near < 0.4 * inside,
|
||||
"{name}: they swing by {near} at the edge"
|
||||
);
|
||||
for from in (0..edge - 27).step_by(9) {
|
||||
let (mean, _) = window(from);
|
||||
assert!((mean - 1.0).abs() < 0.03, "{name}: mean {mean} from {from}");
|
||||
}
|
||||
for (x, r) in ratio.iter().enumerate().skip(edge + 2) {
|
||||
assert!((r - 1.0).abs() < 1e-6, "{name}: {r} at {x}, past the edge");
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -13,6 +13,7 @@ mod clouds;
|
||||
mod convolve;
|
||||
mod distortion;
|
||||
mod earth_detail;
|
||||
mod earth_regions;
|
||||
mod earth_texture;
|
||||
mod earthshine;
|
||||
mod eclipses;
|
||||
|
||||
+136
-2
@@ -1233,6 +1233,139 @@ average at lighting L1 and camera C1 (0.8 DN with the global map alone,
|
||||
1.3 DN with the sharp map; 2.7 DN before the atmosphere, which lowers the
|
||||
contrast). The tests of the maps compare each device with itself.
|
||||
|
||||
**Fine regions.** Over ground a scene names, the surface can come from a 10 m map
|
||||
instead (`EarthRegions` in `earth.rs`, `shaders/regions.wgsl`; made by
|
||||
`resolved earth-regions`, fetched by `scripts/fetch-earth-regions.py`). The source is
|
||||
ESA WorldCover's Sentinel-2 composite for 2021: for each texel, the median of the
|
||||
year's cloud-free Sentinel-2 L2A surface reflectance, with bands B04, B03 and B02 as
|
||||
red, green and blue (B08, near infrared, only finds water). Unlike Blue Marble it is a
|
||||
calibrated reflectance product, but it has no seasons. Its texels are 0.3 arcseconds
|
||||
(9.3 m at the equator), in 1 x 1 degree cells of 12000 x 12000 named by their
|
||||
south-west corner (N29E032 is 29 to 30 N and 32 to 33 E). The bucket has a file only
|
||||
for cells with land in them, about 650 MB each.
|
||||
|
||||
A scene names boxes of latitude and longitude (`[appearance] earth_regions`), and the
|
||||
script fetches the cells that cover them, ahead of time: nothing else is fetched, and
|
||||
the renderer never fetches anything. Each file is checked against the bucket's ETag
|
||||
(the MD5 of its 8 MiB upload parts). A cell with no file is all water and is marked so
|
||||
the CLI does not ask for it again. A cell named but not fetched is drawn from the
|
||||
sharp map, and the CLI names it once.
|
||||
|
||||
Each cell is made into one tile: 6 levels from 9.3 m to 296 m, and a border of 32
|
||||
texels holding its neighbors' texels where their files are on hand, else its own edge
|
||||
repeated. The tile's header says which neighbors it holds, and a cell is remade when
|
||||
one arrives, so a region grown later joins up. As BC6H, a cell takes 194 MB of device
|
||||
memory, made in about 10 s.
|
||||
|
||||
The composite has holes, blocks of open water about 10 km across with no data: 13 %
|
||||
of the texels in N29E032, the north end of the Gulf of Suez, and 56 to 87 % of those
|
||||
over the sea around Gibraltar. They are filled from the water around them: a pyramid
|
||||
of the means of the water texels (those with green brighter than near infrared,
|
||||
McFeeters' water index above 0), blended from the coarsest level down into each
|
||||
level's holes (pull-push). Down to 300 m the pyramid spans the whole group of fetched
|
||||
cells joined to the cell through neighbors, read from the files' 750 x 750 overviews,
|
||||
so neighboring cells fill alike; below that, each cell's own texels. A tile records
|
||||
which group it was filled with and is remade when the group changes. Three ways of
|
||||
filling were tried before:
|
||||
|
||||
- From all the data around a hole: the holes took the color of the desert across the
|
||||
Gulf of Suez, which is narrower than a hole is wide, and showed as brown squares in
|
||||
the sea.
|
||||
- From each cell's own water: water differs from cell to cell (N35W005's is half as
|
||||
bright as N36W005's), and a cell's whole sea changed brightness at its edge, as a
|
||||
straight-edged rectangle east of Gibraltar.
|
||||
- From each cell and its eight neighbors: neighboring cells' windows differ, and left
|
||||
a fainter line along the edge between them.
|
||||
|
||||
What stays is in the data: where the composite has water next to a hole, its water
|
||||
can be darker than the filled sea beside it, in patches with the holes' straight
|
||||
edges (off the Moroccan coast east of Ceuta, in the Gibraltar reference).
|
||||
|
||||
The shader finds a read's cell from its latitude and longitude in degrees, through a
|
||||
table of the 360 x 180 cells. It reads the regions where a pixel covers less ground
|
||||
than their coarsest level's texels and blends into the sharp map across half a
|
||||
level, as the sharp map does into the global one. A read that falls in a cell not
|
||||
fetched takes the nearest fetched cell's edge. The map coordinates the other maps
|
||||
use hold a position to 6e-8, about 0.26 of these texels; in degrees, the position
|
||||
is held to 0.05 texel at Sinai's longitude and 0.18 at 180 degrees. With the map
|
||||
coordinates, two views of the same ground differed by 1.4 % of the peak at a band
|
||||
edge on the meridian.
|
||||
|
||||
The two maps differ in color. Over N29E032 at 460 m, Blue Marble's July is 28 %
|
||||
brighter than WorldCover in red, 25 % in green and 3 % in blue, and its January 9 %,
|
||||
6 % and -20 % (desert, sea and coast together). Around Gibraltar, July has 1.2 times
|
||||
the cells' red and half their blue over land, relative to green: the cells are one
|
||||
year's median, greener, grayer and bluer than a dry July. Drawn as they are, they
|
||||
left Gibraltar's land less red against its sea than the render from Blue Marble's
|
||||
July (1.18 and 1.37 times the sea's red against green, at two places, against 1.27
|
||||
and 1.56; the photo has 2.04 and 2.37), and the matching of the render to the photo
|
||||
could not fix land and sea at once.
|
||||
|
||||
So each cell is brought to the scene's month's colors at 4.4 km, and keeps its own
|
||||
finer detail (`month_matched` in `shaders/regions.wgsl`). Each cell carries a coarse
|
||||
map, its mean color in 25 x 25 blocks with a border of one block from its neighbors
|
||||
(`EarthDetail::coarse`, made with the tile). The shader scales each texel, per
|
||||
channel, by the month's global map over the same block (read with a 4.4 km
|
||||
footprint) over the coarse map there, held to 0.4 to 2.5. Gibraltar then has Blue
|
||||
Marble's July balance with the cells' detail (1.26 and 1.55), and the straight-edged
|
||||
patches in the sea are gone, being darker than the block they sit in. The cost is
|
||||
the cells' calibration: their broad colors are Blue Marble's.
|
||||
|
||||
A texel under a sixth as bright as the month's map over its block keeps its own
|
||||
colors, and is fully matched from a third: water that the 2.45 km map draws as land
|
||||
(a coast it draws elsewhere). Two other ways of telling such places apart were tried
|
||||
on Lake Timsah, which the month's map has as desert: fading the match out where the
|
||||
two maps' blocks differ in brightness by more than 1.5 times, and where they differ
|
||||
in color by more than 1.4 times between channels. The first left the lake's shallows
|
||||
orange, since the lake's block is about as bright in both maps; the second also
|
||||
switched the match off around Gibraltar, where July and the year differ by 2.3 to
|
||||
2.75 times between red and blue. What remains: the shallows of Lake Timsah, brighter
|
||||
than the gate, are tinted orange, and farms newer than Blue Marble's 2004 take its
|
||||
desert's tone.
|
||||
|
||||
Where a region ends with no fetched cell beyond, its cells fade into the map below
|
||||
over 0.1 degrees of ground (11 km, `REGION_FADE` in `gpu/bodies.rs`), measured from
|
||||
the edge, so the change in sharpness does not show as a line; with the colors
|
||||
matched, the brightness does not change there.
|
||||
|
||||
Tests (`tests/gpu/earth_regions.rs`), on the GPU and on lavapipe:
|
||||
|
||||
- The same ground seen where four cells meet and inside one, with the pattern moved
|
||||
0.52 degrees, a whole number of texels at every level and of coarse blocks. A
|
||||
pattern with a dark-to-bright edge on the seam differs by up to 0.38 % and 1.14 % of
|
||||
the peak on the GPU (500 and 2500 km) and 0.98 % and 1.47 % on lavapipe. With
|
||||
borders repeating each cell's own edge, it differs by 24 %.
|
||||
- Views past the coarsest level, and views of ground with no cell, draw what they
|
||||
would without the regions, to the bit.
|
||||
- A uniform cell twice as bright as the month's map draws the map's brightness, to
|
||||
1 % of the peak.
|
||||
- A cell of stripes over a uniform map, seen across its edge: inside, the stripes
|
||||
swing about the map's brightness; within 0.1 degrees of the edge they fade out; past
|
||||
it the map is drawn as it is; averaged over two stripes, the brightness stays the
|
||||
map's to 3 % throughout.
|
||||
|
||||
Cost at 1920 x 1080 with camera C1, the Suez scene at lighting L1 with and without
|
||||
its four cells (medians of two runs each):
|
||||
|
||||
| | Bodies pass | Frame | Setup |
|
||||
|---|---|---|---|
|
||||
| RTX 3060, without | 2.83 to 2.86 ms | 6.37 to 6.45 ms | |
|
||||
| RTX 3060, with | 3.99 to 4.00 ms | 7.30 to 8.09 ms | 2.4 s more |
|
||||
| lavapipe, without | 153 ms | 216 ms | |
|
||||
| lavapipe, with | 257 ms | 316 ms | 1.7 s more |
|
||||
|
||||
The colors' match is 0.3 ms of that on the RTX 3060 (3.69 to 4.00 ms).
|
||||
|
||||
The shader's code for the regions is compiled in only when they are set. Compiled
|
||||
into every pipeline, it cost 0.12 to 0.16 ms in the bodies pass of scenes with no
|
||||
regions (Landsat 8 and ICESat-2), from the larger shader alone; compiled in only when
|
||||
set, those scenes are as fast as before, within 0.05 ms on the RTX 3060 and 10 ms on
|
||||
lavapipe.
|
||||
Memory caps how many cells a run can hold: lavapipe allocates no more than 2 GB for
|
||||
one texture, so the cells are split among 4 arrays of up to 11, 44 cells (8.5 GB) in
|
||||
all, and a run that names more is refused at setup. By their sizes, about 25 fit
|
||||
on a 12 GB device beside the sharp map; that is not measured.
|
||||
|
||||
Two details matter for brightness:
|
||||
|
||||
- Averaging has to happen on reflectance, not on the encoded values image files
|
||||
@@ -2944,7 +3077,8 @@ is a run of the script with `--gaia 16`, at twice the size.
|
||||
- The full Moon's edge as bright as its center (a lunar reflectance model)
|
||||
- Clouds as volumes: soft edges, light through a cloud's side, and gaps below it
|
||||
(lighting L2 draws their tops as relief)
|
||||
- Earth's surface finer than Blue Marble's 460 m: mosaics of tens of meters
|
||||
(Landsat, Sentinel-2) need tiles streamed in as a view needs them
|
||||
- Earth's surface finer than Blue Marble's 460 m beyond the regions a scene names
|
||||
ahead of time (a view that moves on, as in a long simulation, would need cells
|
||||
streamed in), and the fine regions' seasons (one year's median)
|
||||
- The atmosphere seen from inside it (a camera below 100 km)
|
||||
- Motion blur, rolling shutter, stray light
|
||||
@@ -16,7 +16,7 @@ camera_tier = "C1"
|
||||
differences = [
|
||||
"The camera's timestamp is about 17 s off: the landmarks fit with the ISS 133 km further along its track than Horizons puts it at the timestamp.",
|
||||
"Clouds: the day's MODIS map puts a bank of cloud over the western strait and Tangier, where the photo has only a few clouds near Cape Spartel (#76). It is left out of the matching.",
|
||||
"The render is much softer than the photo: the photo resolves about 20 m a pixel at this range, and Blue Marble's sharp map has 460 m texels.",
|
||||
"The ground comes from ESA WorldCover's 10 m Sentinel-2 composite (2021), brought to Blue Marble's July colors at 4.4 km: its land is less red against the sea than the photo's (about 1.3 and 1.55 times the sea's red against green, against 2.0 and 2.4), as it was with Blue Marble alone. The render's pixels cover about 55 m here; the photo resolves about 20 m.",
|
||||
"A small dark disk shows inside the render's cloud bank (L2 relief; see #73).",
|
||||
]
|
||||
|
||||
|
||||
@@ -28,3 +28,8 @@ velocity = [-2442.16, -5810.64, 4363.89]
|
||||
# RMS in the photo's pixels. The first guess looked at the strait (35.95 N, 5.55 W)
|
||||
# with image up to the west.
|
||||
attitude = [0.595618782400, 0.625218378205, -0.347056536633, 0.365912566042]
|
||||
|
||||
# The 10 m cells over the strait and the coasts in view, with room for their fade
|
||||
# (scripts/fetch-earth-regions.py --scene examples/references/gibraltar/scene.toml).
|
||||
[appearance]
|
||||
earth_regions = [[35.1, -6.4, 36.9, -4.7]]
|
||||
@@ -0,0 +1,25 @@
|
||||
# Suez: Earth straight down from 500 km over the head of the Gulf of Suez, north up,
|
||||
# with the Sun 60 degrees up in the south-east. A pixel covers 55 m of ground. At
|
||||
# lighting L1 the ground comes from ESA WorldCover's 10 m Sentinel-2 composite, where
|
||||
# its cells are fetched (scripts/fetch-earth-regions.py --scene
|
||||
# examples/scenes/suez.toml, four cells, 2.6 GB), and from Blue Marble's 460 m map
|
||||
# otherwise. The four cells meet at the image's center. Generated by
|
||||
# scripts/make-example-scenes.py.
|
||||
#
|
||||
# resolved render --scene examples/scenes/suez.toml --lighting L1 --camera C1 --pretty --preview suez.png
|
||||
|
||||
sun = [87679155672, 119715976017, 18972971961]
|
||||
moon = [0, 0, -3.844e8]
|
||||
|
||||
[observer]
|
||||
position = [5005200.16, 3238009.20, 3420373.74]
|
||||
look_at = [4641634.09, 3002807.76, 3170373.74]
|
||||
up = [-0.419809932, -0.271587225, 0.866025404]
|
||||
|
||||
# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh), and the
|
||||
# 10 m cells over these boxes.
|
||||
[appearance]
|
||||
earth_month = 7
|
||||
earth_regions = [[29, 32, 31, 34]]
|
||||
# At lighting L2, the clouds of 2024 July 15 (from scripts/fetch-dated-clouds.py).
|
||||
cloud_date = "2024-07-15"
|
||||
@@ -130,7 +130,13 @@ At lighting L1, Earth's surface comes from NASA's Blue Marble, which
|
||||
from `[appearance] earth_month` in a scene without one. `scripts/fetch-earth.sh
|
||||
--sharp 07` also makes July's sharp map (460 m per texel instead of 2.45 km; 330 MB
|
||||
to download and about 3 minutes to make 4.9 GB), which the CLI and the movie use
|
||||
when it is there. Clouds over Earth come from the Blue Marble's cloud map (Terra
|
||||
when it is there. Finer still, a scene can name boxes of ground in `[appearance]
|
||||
earth_regions` to draw from ESA WorldCover's 10 m Sentinel-2 composite;
|
||||
`scripts/fetch-earth-regions.py --scene FILE` downloads the 1-degree cells that
|
||||
cover them (about 650 MB each) and makes each into `data/earth/regions/` (194 MB
|
||||
of GPU memory a cell). Only the cells named are fetched, and only ahead of time:
|
||||
the CLI draws cells that are not there from the sharp map and names them.
|
||||
`suez.toml` looks straight down from 500 km where four such cells meet. Clouds over Earth come from the Blue Marble's cloud map (Terra
|
||||
MODIS, 2001 July 29), which `scripts/fetch-clouds.sh` downloads (420 MB) and
|
||||
resizes into `data/clouds/blue-marble.png`; `[appearance] clouds = false` leaves
|
||||
them out. Lighting L2 draws the clouds of the scene's day instead (its epoch's, or
|
||||
|
||||
Executable
+180
@@ -0,0 +1,180 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Downloads the 1 x 1 degree cells of Earth's surface at 10 m that scenes name, and
|
||||
makes the renderer's fine regions from them.
|
||||
|
||||
scripts/fetch-earth-regions.py [--scene FILE] [--box S W N E] [CELL ...]
|
||||
|
||||
Cells come from each scene's `[appearance] earth_regions` boxes, from each --box
|
||||
(south, west, north, east in degrees; west east of east crosses the 180 degree
|
||||
meridian), and from names such as N29E032 (29 to 30 N, 32 to 33 E). Only those cells
|
||||
are fetched; nothing else is, and the renderer never fetches anything itself.
|
||||
|
||||
Source: ESA WorldCover 10 m 2021 v200, Sentinel-2 RGBNIR composite: the median of the
|
||||
year's cloud-free Sentinel-2 L2A surface reflectance in bands B04, B03, B02 and B08,
|
||||
as 12000 x 12000 GeoTIFFs of 16-bit values (1e-4 reflectance per unit, 0 where there
|
||||
is no data), one per cell with land in it, from the public bucket
|
||||
s3://esa-worldcover-s2 (Zanaga et al. 2022, doi:10.5281/zenodo.7254221). License:
|
||||
CC BY 4.0; credit "ESA WorldCover project 2021 / Contains modified Copernicus Sentinel
|
||||
data (2021) processed by ESA WorldCover consortium".
|
||||
|
||||
Each file (about 650 MB) goes to data/earth/regions/source/ and is checked against
|
||||
the bucket's ETag, the MD5 of its 8 MiB parts. A cell the bucket has no file for is
|
||||
all water, and gets an empty data/earth/regions/<cell>.none so the renderer does not
|
||||
ask for it again. Then `resolved earth-regions` makes data/earth/regions/<cell>.tiles
|
||||
(194 MB of BC6H, about 10 s a cell; it builds the CLI first), and remakes cells made
|
||||
before whose neighbors have arrived since, so their borders join up and their holes
|
||||
are filled alike.
|
||||
|
||||
Needs only the Python standard library.
|
||||
"""
|
||||
|
||||
import argparse
|
||||
import hashlib
|
||||
import math
|
||||
import subprocess
|
||||
import sys
|
||||
import tomllib
|
||||
import urllib.error
|
||||
import urllib.request
|
||||
from pathlib import Path
|
||||
|
||||
BUCKET = "https://esa-worldcover-s2.s3.eu-central-1.amazonaws.com/rgbnir/2021"
|
||||
PART = 8 * 1024 * 1024 # the bucket's multipart upload part size
|
||||
ROOT = Path(__file__).resolve().parent.parent
|
||||
OUT = ROOT / "data" / "earth" / "regions"
|
||||
SOURCE = OUT / "source"
|
||||
|
||||
|
||||
def name(south, west):
|
||||
"""The cell's name, as N29E032 (as Cell::name in earth.rs)."""
|
||||
ns = "S" if south < 0 else "N"
|
||||
ew = "W" if west < 0 else "E"
|
||||
return f"{ns}{abs(south):02}{ew}{abs(west):03}"
|
||||
|
||||
|
||||
def parse(cell):
|
||||
"""South and west of a name such as N29E032."""
|
||||
ok = len(cell) == 7 and cell[0] in "NS" and cell[3] in "EW"
|
||||
ok = ok and cell[1:3].isdigit() and cell[4:].isdigit()
|
||||
if not ok:
|
||||
sys.exit(f"{cell!r} does not name a cell, as N29E032 does")
|
||||
south = int(cell[1:3]) * (-1 if cell[0] == "S" else 1)
|
||||
west = int(cell[4:]) * (-1 if cell[3] == "W" else 1)
|
||||
if not (-90 <= south < 90 and -180 <= west < 180):
|
||||
sys.exit(f"{cell} is not on Earth")
|
||||
return south, west
|
||||
|
||||
|
||||
def covering(south, west, north, east):
|
||||
"""The cells over a box, as Cell::covering in earth.rs."""
|
||||
valid = -90 <= south < north <= 90 and -180 <= west <= 180
|
||||
valid = valid and -180 <= east <= 180 and west != east
|
||||
if not valid:
|
||||
sys.exit(f"{south}, {west} to {north}, {east} is not a box of latitude and longitude")
|
||||
first, last = math.floor(west), math.ceil(east)
|
||||
if last <= first:
|
||||
last += 360
|
||||
cells = set()
|
||||
for s in range(max(math.floor(south), -90), min(math.ceil(north), 90)):
|
||||
for w in range(first, last):
|
||||
cells.add((s, (w + 180) % 360 - 180))
|
||||
return cells
|
||||
|
||||
|
||||
def file_name(cell):
|
||||
return f"ESA_WorldCover_10m_2021_v200_{cell}_S2RGBNIR.tif"
|
||||
|
||||
|
||||
def url(cell):
|
||||
return f"{BUCKET}/{cell[:3]}/{file_name(cell)}"
|
||||
|
||||
|
||||
def etag(cell):
|
||||
"""The bucket's ETag for the cell, or None if it has no file for it."""
|
||||
request = urllib.request.Request(url(cell), method="HEAD")
|
||||
try:
|
||||
with urllib.request.urlopen(request, timeout=60) as response:
|
||||
return response.headers["ETag"].strip('"')
|
||||
except urllib.error.HTTPError as e:
|
||||
if e.code == 404:
|
||||
return None
|
||||
raise
|
||||
|
||||
|
||||
def local_etag(path, parts):
|
||||
"""The ETag S3 gives a file uploaded in `parts` parts of PART bytes, or whole."""
|
||||
if parts == 0:
|
||||
return hashlib.md5(path.read_bytes()).hexdigest()
|
||||
digests = []
|
||||
with open(path, "rb") as f:
|
||||
while chunk := f.read(PART):
|
||||
digests.append(hashlib.md5(chunk).digest())
|
||||
return f"{hashlib.md5(b''.join(digests)).hexdigest()}-{len(digests)}"
|
||||
|
||||
|
||||
def matches(path, tag):
|
||||
parts = int(tag.split("-")[1]) if "-" in tag else 0
|
||||
return local_etag(path, parts) == tag
|
||||
|
||||
|
||||
def fetch(cell):
|
||||
"""Downloads the cell's file if it is not here and checked; False if the bucket
|
||||
has none."""
|
||||
tag = etag(cell)
|
||||
if tag is None:
|
||||
(OUT / f"{cell}.none").touch()
|
||||
print(f"{cell}: all water, no file")
|
||||
return False
|
||||
path = SOURCE / file_name(cell)
|
||||
if path.is_file() and matches(path, tag):
|
||||
return True
|
||||
print(f"{cell}: fetching")
|
||||
part = path.with_suffix(".part")
|
||||
with urllib.request.urlopen(url(cell), timeout=60) as response, open(part, "wb") as f:
|
||||
while chunk := response.read(1 << 20):
|
||||
f.write(chunk)
|
||||
if not matches(part, tag):
|
||||
part.unlink()
|
||||
sys.exit(f"{cell}: download does not match the bucket's ETag {tag}")
|
||||
part.rename(path)
|
||||
return True
|
||||
|
||||
|
||||
def main():
|
||||
parser = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
|
||||
parser.add_argument("--scene", action="append", default=[], type=Path)
|
||||
parser.add_argument("--box", action="append", default=[], nargs=4, type=float)
|
||||
parser.add_argument("cells", nargs="*")
|
||||
args = parser.parse_args()
|
||||
|
||||
wanted = {parse(c) for c in args.cells}
|
||||
for box in args.box:
|
||||
wanted |= covering(*box)
|
||||
for scene in args.scene:
|
||||
with open(scene, "rb") as f:
|
||||
boxes = tomllib.load(f).get("appearance", {}).get("earth_regions", [])
|
||||
if not boxes:
|
||||
print(f"{scene} names no earth_regions")
|
||||
for box in boxes:
|
||||
wanted |= covering(*box)
|
||||
if not wanted:
|
||||
sys.exit("no cells named: give --scene, --box or cell names")
|
||||
|
||||
OUT.mkdir(parents=True, exist_ok=True)
|
||||
SOURCE.mkdir(parents=True, exist_ok=True)
|
||||
names = sorted(name(s, w) for s, w in wanted)
|
||||
fetched = [c for c in names if not (OUT / f"{c}.none").is_file() and fetch(c)]
|
||||
# Cells made before are passed too: those next to new ones are remade.
|
||||
made = {p.stem for p in OUT.glob("*.tiles")}
|
||||
build = sorted(set(fetched) | made)
|
||||
if build:
|
||||
subprocess.run(
|
||||
["cargo", "run", "--quiet", "--release", "-p", "resolved-cli", "--",
|
||||
"earth-regions", "--source", str(SOURCE), "--out", str(OUT), *build],
|
||||
cwd=ROOT, check=True,
|
||||
)
|
||||
print(f"{len(fetched)} cells in {OUT}")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -675,6 +675,54 @@ attitude = {fmt(quaternion(attitude), 9)}
|
||||
"""
|
||||
|
||||
|
||||
# Suez: Earth straight down from 500 km over the head of the Gulf of Suez, where the
|
||||
# four 1-degree cells of the fine regions from 29 to 31 N and 32 to 34 E meet, with the
|
||||
# Sun 60 degrees up in the south-east.
|
||||
SUEZ = (30.0, 32.9) # geodetic latitude and longitude, degrees
|
||||
SUEZ_SUN = (135.0, 60.0) # the Sun's azimuth from north and elevation, degrees
|
||||
|
||||
|
||||
def suez():
|
||||
lat, lon = (math.radians(x) for x in SUEZ)
|
||||
e2 = 1.0 - (6356752.314245 / RE) ** 2
|
||||
n = RE / math.sqrt(1.0 - e2 * math.sin(lat) ** 2)
|
||||
up = [math.cos(lat) * math.cos(lon), math.cos(lat) * math.sin(lon), math.sin(lat)]
|
||||
ground = [n * up[0], n * up[1], n * (1.0 - e2) * up[2]]
|
||||
east = norm([-math.sin(lon), math.cos(lon), 0.0])
|
||||
north = cross(up, east)
|
||||
az, el = (math.radians(x) for x in SUEZ_SUN)
|
||||
sun = [math.cos(el) * (math.cos(az) * north[i] + math.sin(az) * east[i]) + math.sin(el) * up[i]
|
||||
for i in range(3)]
|
||||
camera = [ground[i] + 500e3 * up[i] for i in range(3)]
|
||||
fmt = lambda v, d=2: "[" + ", ".join(f"{x:.{d}f}" for x in v) + "]"
|
||||
return f"""# Suez: Earth straight down from 500 km over the head of the Gulf of Suez, north up,
|
||||
# with the Sun 60 degrees up in the south-east. A pixel covers 55 m of ground. At
|
||||
# lighting L1 the ground comes from ESA WorldCover's 10 m Sentinel-2 composite, where
|
||||
# its cells are fetched (scripts/fetch-earth-regions.py --scene
|
||||
# examples/scenes/suez.toml, four cells, 2.6 GB), and from Blue Marble's 460 m map
|
||||
# otherwise. The four cells meet at the image's center. Generated by
|
||||
# scripts/make-example-scenes.py.
|
||||
#
|
||||
# resolved render --scene examples/scenes/suez.toml --lighting L1 --camera C1 --pretty --preview suez.png
|
||||
|
||||
sun = {fmt([AU * x for x in sun], 0)}
|
||||
moon = [0, 0, -3.844e8]
|
||||
|
||||
[observer]
|
||||
position = {fmt(camera)}
|
||||
look_at = {fmt(ground)}
|
||||
up = {fmt(north, 9)}
|
||||
|
||||
# At lighting L1, Earth's surface shows July (from scripts/fetch-earth.sh), and the
|
||||
# 10 m cells over these boxes.
|
||||
[appearance]
|
||||
earth_month = 7
|
||||
earth_regions = [[29, 32, 31, 34]]
|
||||
# At lighting L2, the clouds of 2024 July 15 (from scripts/fetch-dated-clouds.py).
|
||||
cloud_date = "2024-07-15"
|
||||
"""
|
||||
|
||||
|
||||
# For the long-exposure cameras, which leave the example camera's 1 in 10,000 hot.
|
||||
CALIBRATED = """# Calibrated: hot pixels mapped out, all but 1 value in a million.
|
||||
hot_pixels = [1e-6, 50]
|
||||
@@ -813,6 +861,9 @@ print(f"wrote {path}")
|
||||
path = out / "sunset.toml"
|
||||
path.write_text(with_presets(sunset()))
|
||||
print(f"wrote {path}")
|
||||
path = out / "suez.toml"
|
||||
path.write_text(suez())
|
||||
print(f"wrote {path}")
|
||||
path = out / "stars-pleiades.toml"
|
||||
path.write_text(stars_pleiades())
|
||||
print(f"wrote {path}")
|
||||
|
||||
Reference in new issue
Block a user