schuwi/dwd_cloud_cover_rs
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Improve rendering: barycentric interpolation, Gaussian blur, white-cloud palette

Browse Source 
Replace nearest-neighbour BFS flood fill with Delaunay triangulation
(spade) + barycentric interpolation for smooth cloud boundaries.
Add NaN-aware separable Gaussian blur (σ=8 px, rayon-parallelised) to
remove triangulation facets. Switch colour scheme from blue bands to a
continuous green-to-white opacity blend matching the DWD app style.
Pixel interpolation loop is also parallelised with rayon.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Schuwi
2026-03-06 22:19:46 +01:00
parent 4ee33a882d
commit ee9b039d56
3 changed files with 253 additions and 69 deletions
Download Patch File Download Diff File Expand all files Collapse all files

220
src/render.rs
View File

@@ -2,7 +2,8 @@ use crate::projection::OrthoProjection;
use anyhow::{Context, Result};
use font8x8::{UnicodeFonts, BASIC_FONTS};
use plotters::prelude::*;
use std::collections::VecDeque;
use rayon::prelude::*;
use spade::{DelaunayTriangulation, HasPosition, Point2, Triangulation};
use std::path::Path;
const IMG_WIDTH: u32 = 900;
@@ -14,20 +15,8 @@ const CENTER_LON: f64 = 13.08;
const HALF_W_DEG: f64 = 0.8; // horizontal half-extent in degrees of arc
const HALF_H_DEG: f64 = 0.4; // vertical half-extent
/// Cloud-cover thresholds and their associated RGB colours.
/// Values below the first threshold are drawn as the background colour.
const LEVELS: [(f32, (u8, u8, u8)); 7] = [
(1.0, (208, 232, 248)), // very thin
(2.0, (160, 200, 232)),
(5.0, (112, 168, 216)),
(10.0, ( 64, 128, 192)),
(20.0, ( 32, 80, 160)),
(50.0, ( 16, 32, 96)),
(100.0, ( 5, 10, 48)),
];
const BACKGROUND: RGBColor = RGBColor(240, 248, 255); // pale sky-blue
const LAND_BG: RGBColor = RGBColor(232, 236, 220); // muted green for land
const LAND_BG: RGBColor = RGBColor(142, 178, 35); // strong green for land
const BLACK: RGBColor = RGBColor( 0, 0, 0);
const RED: RGBColor = RGBColor(220, 30, 30);
@@ -42,16 +31,28 @@ const LOCATIONS: &[(&str, f64, f64)] = &[
("Berlin (Mitte)", 13.39, 52.52),
];
#[derive(Clone, Copy)]
struct CloudVertex {
position: Point2<f64>,
cloud: f32,
}
impl HasPosition for CloudVertex {
type Scalar = f64;
fn position(&self) -> Point2<f64> {
self.position
}
}
/// Map cloud cover (0–100 %) to a colour by blending the land background toward
/// white. Returns `None` for near-zero cover so the land background shows through.
fn cloud_color(cover: f32) -> Option<RGBColor> {
if cover < 1.0 {
return None; // clear sky — use background
return None;
}
for &(threshold, (r, g, b)) in LEVELS.iter().rev() {
if cover >= threshold {
return Some(RGBColor(r, g, b));
}
}
Some(RGBColor(LEVELS[0].1 .0, LEVELS[0].1 .1, LEVELS[0].1 .2))
let t = (cover / 100.0).clamp(0.0, 1.0);
let blend = |base: u8| -> u8 { (base as f32 + t * (255.0 - base as f32)).round() as u8 };
Some(RGBColor(blend(LAND_BG.0), blend(LAND_BG.1), blend(LAND_BG.2)))
}
/// Draw ASCII text using the built-in 8×8 pixel font.
@@ -91,6 +92,57 @@ fn draw_pixel_text(
}
}
/// Separable Gaussian blur on a flat f32 grid.
///
/// NaN pixels are treated as transparent: they contribute zero weight and the
/// result is renormalised over non-NaN neighbours only. Pixels with no non-NaN
/// neighbours remain NaN.
fn gaussian_blur(grid: &[f32], w: usize, h: usize, sigma: f32) -> Vec<f32> {
let radius = (3.0 * sigma).ceil() as usize;
let kernel: Vec<f32> = (0..=2 * radius)
.map(|i| {
let x = i as f32 - radius as f32;
(-x * x / (2.0 * sigma * sigma)).exp()
})
.collect();
// Horizontal pass
let temp: Vec<f32> = (0..h * w)
.into_par_iter()
.map(|idx| {
let row = idx / w;
let col = idx % w;
let mut val = 0.0f32;
let mut wsum = 0.0f32;
for (ki, &k) in kernel.iter().enumerate() {
let c = col as i64 + ki as i64 - radius as i64;
if c < 0 || c >= w as i64 { continue; }
let v = grid[row * w + c as usize];
if !v.is_nan() { val += k * v; wsum += k; }
}
if wsum > 0.0 { val / wsum } else { f32::NAN }
})
.collect();
// Vertical pass
(0..h * w)
.into_par_iter()
.map(|idx| {
let row = idx / w;
let col = idx % w;
let mut val = 0.0f32;
let mut wsum = 0.0f32;
for (ki, &k) in kernel.iter().enumerate() {
let r = row as i64 + ki as i64 - radius as i64;
if r < 0 || r >= h as i64 { continue; }
let v = temp[r as usize * w + col];
if !v.is_nan() { val += k * v; wsum += k; }
}
if wsum > 0.0 { val / wsum } else { f32::NAN }
})
.collect()
}
/// Render one forecast frame as a PNG.
///
/// - `lats`, `lons`, `cloud_cover` must all have the same length.
@@ -119,54 +171,89 @@ pub fn render_frame(
map_area.fill(&LAND_BG).context("fill map area")?;
// Nearest-neighbour fill via multi-source BFS.
//
// Each data point seeds one pixel in a flat grid. BFS then propagates each
// value outward to all 4-connected NaN neighbours, so every reachable pixel
// ends up holding the value of the closest source point (by Manhattan distance).
// Complexity: O(map_w × map_h) — ~410 k array ops, no hashing.
let grid_w = map_w as usize;
let grid_h = map_h as usize;
let mut grid = vec![f32::NAN; grid_w * grid_h];
let mut queue: VecDeque<(i32, i32)> = VecDeque::new();
// Seed the grid with projected data points.
// Build a Delaunay triangulation in continuous pixel-space coordinates.
// Include all grid points that project within a small margin outside the image
// bounds so that boundary pixels interpolate cleanly without edge artefacts.
const TRI_MARGIN: f64 = 10.0;
let mut tri = DelaunayTriangulation::<CloudVertex>::new();
for i in 0..lats.len() {
let cover = cloud_cover[i];
if cover.is_nan() {
continue;
}
let Some((px, py)) = proj.project(lats[i] as f64, lons[i] as f64) else { continue };
let Some((col, row)) = proj.to_pixel(px, py, map_w, map_h) else { continue };
let idx = row as usize * grid_w + col as usize;
if grid[idx].is_nan() {
grid[idx] = cover;
queue.push_back((col, row));
// Continuous pixel coordinates — no rounding.
let col_f = (px + HALF_W_DEG) / (2.0 * HALF_W_DEG) * map_w as f64;
let row_f = (-py + HALF_H_DEG) / (2.0 * HALF_H_DEG) * map_h as f64;
if col_f < -TRI_MARGIN
|| col_f > map_w as f64 + TRI_MARGIN
|| row_f < -TRI_MARGIN
|| row_f > map_h as f64 + TRI_MARGIN
{
continue;
}
let _ = tri.insert(CloudVertex { position: Point2::new(col_f, row_f), cloud: cover });
}
// BFS flood fill.
const DIRS: [(i32, i32); 4] = [(-1, 0), (1, 0), (0, -1), (0, 1)];
while let Some((col, row)) = queue.pop_front() {
let cover = grid[row as usize * grid_w + col as usize];
for (dc, dr) in DIRS {
let nc = col + dc;
let nr = row + dr;
if nc < 0 || nc >= map_w as i32 || nr < 0 || nr >= map_h as i32 {
continue;
// For each output pixel, locate the containing triangle and barycentric-interpolate
// the cloud cover from its three vertices. The triangulation is immutable at this
// point, so the work is embarrassingly parallel via rayon.
let grid_w = map_w as usize;
let grid_h = map_h as usize;
let cover_grid: Vec<f32> = (0..grid_h * grid_w)
.into_par_iter()
.map(|i| {
let col = (i % grid_w) as f64 + 0.5;
let row = (i / grid_w) as f64 + 0.5;
let p = Point2::new(col, row);
use spade::PositionInTriangulation::*;
match tri.locate(p) {
OnVertex(v) => tri.vertex(v).data().cloud,
OnEdge(e) => {
// Linear interpolation along the edge.
let edge = tri.directed_edge(e);
let a = edge.from();
let b = edge.to();
let (ax, ay) = (a.position().x, a.position().y);
let (bx, by) = (b.position().x, b.position().y);
let ab2 = (bx - ax).powi(2) + (by - ay).powi(2);
let t = if ab2 == 0.0 {
0.0
} else {
(((p.x - ax) * (bx - ax) + (p.y - ay) * (by - ay)) / ab2)
.clamp(0.0, 1.0)
};
a.data().cloud * (1.0 - t) as f32 + b.data().cloud * t as f32
}
OnFace(f) => {
let face = tri.face(f);
let [v0, v1, v2] = face.vertices();
let (ax, ay) = (v0.position().x, v0.position().y);
let (bx, by) = (v1.position().x, v1.position().y);
let (cx, cy) = (v2.position().x, v2.position().y);
let denom = (by - cy) * (ax - cx) + (cx - bx) * (ay - cy);
if denom == 0.0 {
return v0.data().cloud;
}
let l0 = ((by - cy) * (p.x - cx) + (cx - bx) * (p.y - cy)) / denom;
let l1 = ((cy - ay) * (p.x - cx) + (ax - cx) * (p.y - cy)) / denom;
let l2 = 1.0 - l0 - l1;
v0.data().cloud * l0 as f32
+ v1.data().cloud * l1 as f32
+ v2.data().cloud * l2 as f32
}
OutsideOfConvexHull(_) | NoTriangulation => f32::NAN,
}
let nidx = nr as usize * grid_w + nc as usize;
if grid[nidx].is_nan() {
grid[nidx] = cover;
queue.push_back((nc, nr));
}
}
}
})
.collect();
// Paint the filled grid.
const BLUR_SIGMA: f32 = 8.0;
let cover_grid = gaussian_blur(&cover_grid, grid_w, grid_h, BLUR_SIGMA);
// Paint the interpolated grid.
for row in 0..map_h as i32 {
for col in 0..map_w as i32 {
let cover = grid[row as usize * grid_w + col as usize];
let cover = cover_grid[row as usize * grid_w + col as usize];
if let Some(color) = cloud_color(cover) {
map_area.draw_pixel((col, row), &color).ok();
}
@@ -231,22 +318,17 @@ fn draw_legend(
);
};
// "< 1%" row (clear sky)
draw_entry("< 1%", LAND_BG, y as i32 + 22);
// Remaining threshold rows — use ASCII hyphen instead of en-dash
let entries = [
(1.0f32, "1-2%"),
(2.0, "2-5%"),
(5.0, "5-10%"),
(10.0, "10-20%"),
(20.0, "20-50%"),
(50.0, "50-100%"),
(0.0f32, "< 1%"),
(25.0, "~25%"),
(50.0, "~50%"),
(75.0, "~75%"),
(100.0, "100%"),
];
for (i, &(lo, label)) in entries.iter().enumerate() {
let row_y = y as i32 + 22 + (i as i32 + 1) * (box_h + 2);
let color = cloud_color(lo).unwrap_or(BACKGROUND);
for (i, &(cover, label)) in entries.iter().enumerate() {
let row_y = y as i32 + 22 + i as i32 * (box_h + 2);
let color = cloud_color(cover).unwrap_or(LAND_BG);
draw_entry(label, color, row_y);
}