1use crate::retained::Retained;
4use crate::{
5 MapError,
6 model::{self, CrsIdentity, NodataPolicy, QueryConvention},
7 sink,
8 sloping::{self, Cell, Corner, Face, Mesh},
9 terrain::NativeTerrain,
10};
11use rayon::prelude::*;
12use serde::Deserialize;
13use serde_json::{Value, json};
14use std::{
15 collections::HashSet,
16 fs,
17 io::Write,
18 path::Path,
19 sync::{Arc, Mutex},
20 time::Instant,
21};
22
23type CachedGeometry = (Vec<Mesh>, Vec<Vec<Option<sloping::CellStencil>>>);
24type WindowKey = (
25 String,
26 String,
27 [usize; 4],
28 Vec<[usize; 2]>,
29 Vec<u8>,
30 Option<Vec<u8>>,
31);
32struct WindowData {
33 owner: Vec<u32>,
34 ground: Vec<f32>,
35 weights: Option<sloping::PixelWeights>,
36}
37
38#[derive(Default)]
41pub struct PreparedSession {
42 terrain: Option<(String, Arc<NativeTerrain>)>,
43 geometry: Option<(String, Arc<CachedGeometry>)>,
44 windows: Retained<WindowKey, WindowData>,
45 hydraulics: Mutex<Retained<String, Vec<sloping::PreparedProfile>>>,
46 pool: Option<(usize, Arc<rayon::ThreadPool>)>,
47}
48
49fn identity(path: &Path) -> Result<String, MapError> {
50 let path = path.canonicalize()?;
51 let metadata = path.metadata()?;
52 Ok(format!(
53 "{:?}:{}:{:?}",
54 path,
55 metadata.len(),
56 metadata.modified()?
57 ))
58}
59fn geometry_key(req: &Request, base: &Path) -> Result<String, MapError> {
60 let mut key = String::new();
61 for filename in &req.mesh_files {
62 let path = model::resolve(base, filename);
63 key.push_str(&identity(&path)?);
64 let document: Value = serde_json::from_reader(fs::File::open(&path)?)
65 .map_err(|e| MapError::Invalid(e.to_string()))?;
66 if let Some(arrays) = document["arrays"].as_object() {
67 for descriptor in arrays.values() {
68 if let Some(name) = descriptor["path"].as_str() {
69 key.push_str(&identity(&model::resolve(path.parent().unwrap(), name))?);
70 }
71 }
72 }
73 }
74 Ok(key)
75}
76
77pub const SCHEMA: &str = "ras-commander-sloping-request/1";
78
79#[derive(Deserialize)]
80#[serde(deny_unknown_fields)]
81pub struct Request {
82 pub schema: String,
83 #[serde(default)]
84 pub options: sloping::Options,
85 pub terrain: String,
86 pub epsg: u32,
87 pub crs_identity: CrsIdentity,
88 pub terrain_transform: [f64; 6],
89 pub window: [usize; 4],
90 pub geometry_file: String,
91 #[serde(default)]
92 pub aoi_file: Option<String>,
93 pub mesh_files: Vec<String>,
94 pub owner_cells: Vec<[usize; 2]>,
95 pub profiles: Vec<Profile>,
96 pub threads: usize,
97 #[serde(default = "one_profile")]
98 pub max_inflight_profiles: usize,
99 pub max_working_bytes: u64,
100 #[serde(default)]
101 pub stream_outputs: bool,
102 #[serde(default = "window_budget")]
104 pub window_cache_bytes: u64,
105 #[serde(default)]
107 pub hydraulic_cache_bytes: u64,
108 #[serde(default)]
110 pub geotiff_outputs: bool,
111 #[serde(default)]
115 pub velocity_normals_files: Option<Vec<String>>,
116
117 #[serde(default)]
119 pub terrain_scope: Option<crate::terrain_scope::TerrainScope>,
120}
121
122fn one_profile() -> usize {
123 1
124}
125fn window_budget() -> u64 {
126 64 << 20
127}
128
129fn admit_profiles(
130 requested: usize,
131 threads: usize,
132 profiles: usize,
133 budget: u64,
134 serial: u64,
135 additional: u64,
136) -> Result<(usize, u64), MapError> {
137 if serial > budget || additional == 0 || requested == 0 || threads == 0 || profiles == 0 {
138 return Err(MapError::Budget(
139 "sloping decoder plan exceeds budget".into(),
140 ));
141 }
142 let limit = requested.min(threads).min(profiles);
143 let extra = ((budget - serial) / additional).min((limit - 1) as u64);
144 Ok((extra as usize + 1, serial + extra * additional))
145}
146
147fn map_batch<T: Sync, R: Send>(
148 batch: &[T],
149 pool: &rayon::ThreadPool,
150 compute: impl Fn(&T) -> Result<R, MapError> + Sync + Send,
151) -> Result<Vec<R>, MapError> {
152 pool.install(|| batch.par_iter().map(compute).collect())
153}
154
155struct RenderedProfile {
156 depth: Vec<f32>,
157 wse: Vec<f32>,
158 velocity: Option<Vec<f32>>,
159 hydraulic_seconds: f64,
160 render_seconds: f64,
161 hydraulic_reused: bool,
162 finite_pixels: usize,
163 flooded_pixels: usize,
164}
165
166#[allow(clippy::too_many_arguments)] fn render(
168 profile: &Profile,
169 base: &Path,
170 meshes: &[Mesh],
171 stencils: &[Vec<Option<sloping::CellStencil>>],
172 window: &sloping::RasterWindow<'_>,
173 pool: &rayon::ThreadPool,
174 options: sloping::Options,
175 weights: Option<&sloping::PixelWeights>,
176 normals: Option<&[Vec<crate::velocity::FaceNormal>]>,
177 foreign: &[usize],
178 hydraulic_cache: Option<&Mutex<Retained<String, Vec<sloping::PreparedProfile>>>>,
179) -> Result<RenderedProfile, MapError> {
180 let tick = Instant::now();
181 let cache_key = if hydraulic_cache.is_some() && normals.is_none() {
182 let mut key =
183 serde_json::to_string(&options).map_err(|e| MapError::Invalid(e.to_string()))?;
184 for file in &profile.values_files {
185 key.push_str(&identity(&model::resolve(base, file))?);
186 }
187 Some(key)
188 } else {
189 None
190 };
191 let cached = match (&cache_key, hydraulic_cache) {
192 (Some(key), Some(cache)) => cache.lock().unwrap().get(key),
193 _ => None,
194 };
195 let hydraulic_reused = cached.is_some();
196 let mut new_prepared = Vec::new();
197 let mut cell_values = Vec::new();
198 if !hydraulic_reused {
199 for ((mesh, stencil), filename) in meshes.iter().zip(stencils).zip(&profile.values_files) {
200 let wse = values(&model::resolve(base, filename), mesh.cells.len())?;
201 new_prepared.push(sloping::prepare_profile_with_options(
202 mesh, stencil, &wse, 0.001, pool, options,
203 )?);
204 if normals.is_some() {
205 cell_values.push(wse);
206 }
207 }
208 }
209 let prepared = cached.unwrap_or_else(|| Arc::new(new_prepared));
210 if !hydraulic_reused && let (Some(key), Some(cache)) = (cache_key, hydraulic_cache) {
211 let bytes = prepared.iter().map(|p| p.retained_bytes()).sum::<u64>()
212 + (prepared.capacity() * std::mem::size_of::<sloping::PreparedProfile>()) as u64
213 + key.capacity() as u64
214 + (2 * std::mem::size_of::<(String, Arc<Vec<sloping::PreparedProfile>>, u64)>()
215 + std::mem::size_of::<Vec<sloping::PreparedProfile>>()
216 + 2 * std::mem::size_of::<usize>()) as u64;
217 cache.lock().unwrap().insert(key, prepared.clone(), bytes);
218 }
219 crate::terrain_scope::deny_wet_foreign(foreign, window.owner, window.owner_cells, |m, c| {
220 prepared[m].is_wet(c)
221 })?;
222 let hydraulic_seconds = tick.elapsed().as_secs_f64();
223 let tick = Instant::now();
224 let (mut depth, mut wse) =
225 sloping::render_profile_cached(stencils, &prepared, window, 0.001, pool, weights)?;
226 let render_seconds = tick.elapsed().as_secs_f64();
227 let (finite_pixels, flooded_pixels) = crate::finalize::apply(&mut depth, &mut wse, pool);
228 let velocity = match (normals, &profile.face_velocity_files) {
229 (Some(normals), Some(files)) => {
230 let mut face_velocity = Vec::new();
231 for (mesh, file) in meshes.iter().zip(files) {
232 face_velocity.push(crate::velocity::read_face_velocity(
233 &model::resolve(base, file),
234 mesh.faces.len(),
235 )?);
236 }
237 let inputs: Vec<_> = meshes
238 .iter()
239 .zip(stencils)
240 .zip(&cell_values)
241 .zip(normals)
242 .zip(&face_velocity)
243 .map(|((((mesh, stencil), wse), normal), velocity)| {
244 (
245 mesh,
246 stencil.as_slice(),
247 wse.as_slice(),
248 normal.as_slice(),
249 velocity.as_slice(),
250 )
251 })
252 .collect();
253 Some(crate::velocity::render(&inputs, window, &depth, pool)?)
254 }
255 _ => None,
256 };
257 Ok(RenderedProfile {
258 depth,
259 wse,
260 velocity,
261 hydraulic_seconds,
262 render_seconds,
263 hydraulic_reused,
264 finite_pixels,
265 flooded_pixels,
266 })
267}
268
269#[derive(Deserialize)]
270#[serde(deny_unknown_fields)]
271pub struct Profile {
272 pub id: String,
273 pub values_files: Vec<String>,
274 #[serde(default)]
276 pub face_velocity_files: Option<Vec<String>>,
277}
278
279#[derive(Deserialize)]
280pub struct Geometry {
281 #[serde(default)]
282 pub vertex_elevations: Option<Vec<f32>>,
283 pub cells: Vec<InputCell>,
284 pub faces: Vec<InputFace>,
285 pub vertices: Vec<InputVertex>,
286}
287#[derive(Deserialize)]
288pub struct InputCell {
289 pub point: [f64; 2],
290 pub minimum_elevation: f32,
291 #[serde(rename = "virtual")]
292 pub virtual_cell: bool,
293 pub faces: Vec<usize>,
294}
295#[derive(Deserialize)]
296pub struct InputFace {
297 pub cells: [usize; 2],
298 pub vertices: [usize; 2],
299 pub minimum_elevation: f32,
300 pub perimeter: bool,
301 pub polyline: Vec<[f64; 2]>,
302}
303#[derive(Deserialize)]
304pub struct InputVertex {
305 pub point: [f64; 2],
306 pub faces: Vec<usize>,
307 pub ccw: Vec<usize>,
308 pub cw: Vec<usize>,
309}
310
311impl Geometry {
312 pub fn prepare(&self) -> Result<(Mesh, Vec<Option<sloping::CellStencil>>), MapError> {
313 if self
314 .cells
315 .iter()
316 .any(|c| c.virtual_cell != (c.faces.len() == 1))
317 {
318 return Err(MapError::Invalid(
319 "sloping virtual-cell flag contradicts topology".into(),
320 ));
321 }
322 let faces =
323 self.faces
324 .iter()
325 .map(|f| {
326 let a = self.cells.get(f.cells[0]).ok_or_else(|| {
327 MapError::Invalid("sloping cell address out of range".into())
328 })?;
329 let b = self.cells.get(f.cells[1]).ok_or_else(|| {
330 MapError::Invalid("sloping cell address out of range".into())
331 })?;
332 Ok(Face {
333 cells: f.cells,
334 vertices: f.vertices,
335 minimum: f.minimum_elevation,
336 perimeter: f.perimeter,
337 application: sloping::face_application(&f.polyline, [a.point, b.point])?,
338 })
339 })
340 .collect::<Result<Vec<_>, MapError>>()?;
341 let mesh = Mesh {
342 cells: self
343 .cells
344 .iter()
345 .map(|c| Cell {
346 minimum: c.minimum_elevation,
347 virtual_cell: c.virtual_cell,
348 })
349 .collect(),
350 faces,
351 corners: self
352 .vertices
353 .iter()
354 .map(|v| Corner {
355 point: v.point,
356 faces: v.faces.clone(),
357 ccw: v.ccw.clone(),
358 cw: v.cw.clone(),
359 })
360 .collect(),
361 };
362 let cell_faces = self
363 .cells
364 .iter()
365 .map(|c| c.faces.clone())
366 .collect::<Vec<_>>();
367 let mut stencils = sloping::prepare_stencils(&mesh, &cell_faces)?;
368 for (cell, stencil) in stencils.iter_mut().enumerate() {
369 if let Some(stencil) = stencil {
370 for &face in &self.cells[cell].faces {
371 let f = &self.faces[face];
372 if f.cells[0] == cell {
373 stencil.physical_ring.extend(f.polyline.iter().copied());
374 } else {
375 stencil
376 .physical_ring
377 .extend(f.polyline.iter().rev().copied());
378 }
379 }
380 stencil.physical_ring.dedup();
381 }
382 }
383 if let Some(elevations) = &self.vertex_elevations {
384 sloping::set_vertex_elevations(&mut stencils, elevations, mesh.corners.len())?;
385 }
386 Ok((mesh, stencils))
387 }
388}
389
390pub fn method_record() -> Value {
391 json!({"method_id":"rasmapper-plain-sloping/1", "backend":"rust", "experimental":true,
392 "native_version_qualified":"7.0.1", "stream_outputs":true, "native_ownership":true,"aoi_mask":crate::aoi::METHOD,
393 "enhanced_options":true,"enhanced_method_id":"rasmapper-sloping-with-faces/1",
394 "depth_weighted_faces":true,
395 "native_donation_history":true,
396 "bounded_profile_concurrency":true,
397 "persistent_session":true,"multi_window_cache":true,"hydraulic_cache":true,"native_geotiff_outputs":true,"bounded_pixel_weights":true,
398 "velocity":crate::velocity::METHOD_ID,
399
400 "terrain_scope":crate::terrain_scope::METHOD_ID,
401 "options":{"CellRenderMode":"Sloping","CellStencilMode":"JustFacepoints",
402 "FacepointAdjustmentMode":"None","ShallowBehaviorMode":"None",
403 "UseDepthWeightedFaces":false,"FaceWSMode":"Adjusted"},
404 "native_plot_tolerance":0.001_f32,"extra_depth_cutoff":null,
405 "scope":"2D cells; plain cell-corner sloping; no 1D, slopingPretty or additional rendering options",
406 "reuse":"geometry, terrio ownership and terrain prepared once; bounded Rayon pixel and reconstruction workers"})
407}
408
409fn selected_method(options: sloping::Options) -> Value {
410 let mut record = method_record();
411 if options.render_mode == sloping::RenderMode::WithFaces {
412 record["method_id"] = "rasmapper-sloping-with-faces/1".into();
413 record["options"]["CellStencilMode"] = "WithFaces".into();
414 record["options"]["ShallowBehaviorMode"] = if options.reduce_shallow_to_horizontal {
415 "ReduceToHorizontal"
416 } else {
417 "None"
418 }
419 .into();
420 record["options"]["UseDepthWeightedFaces"] = options.use_depth_weighted_faces.into();
421 record["scope"] = "physical 2D cells; cell-corner and face sloping; no 1D or standalone storage-area surfaces".into();
422 }
423 record
424}
425
426fn values(path: &Path, count: usize) -> Result<Vec<f32>, MapError> {
427 let expected = count
428 .checked_mul(4)
429 .ok_or_else(|| MapError::Budget("sloping values size overflow".into()))?;
430 if fs::metadata(path)?.len() != expected as u64 {
431 return Err(MapError::Invalid("sloping profile length mismatch".into()));
432 }
433 let raw = fs::read(path)?;
434 if raw.len() != expected {
435 return Err(MapError::Invalid("sloping profile changed size".into()));
436 }
437 let values: Vec<_> = raw
438 .as_chunks::<4>()
439 .0
440 .iter()
441 .map(|v| f32::from_le_bytes(*v))
442 .collect();
443 if values.iter().any(|v| !v.is_finite()) {
444 return Err(MapError::Invalid(
445 "sloping profiles require finite cell WSE including virtual cells".into(),
446 ));
447 }
448 Ok(values)
449}
450
451pub fn run(req: &Request, base: &Path, out: &Path) -> Result<Value, MapError> {
452 run_with_session(req, base, out, None)
453}
454
455pub fn run_with_session(
456 req: &Request,
457 base: &Path,
458 out: &Path,
459 mut session: Option<&mut PreparedSession>,
460) -> Result<Value, MapError> {
461 req.options.validate()?;
462 if req.geotiff_outputs && (!req.stream_outputs || req.velocity_normals_files.is_some()) {
463 return Err(MapError::Unsupported(
464 "native GeoTIFF requires stream_outputs and no velocity".into(),
465 ));
466 }
467 if let Some(scope) = &req.terrain_scope {
468 scope.validate()?;
469 if req.options.render_mode == sloping::RenderMode::WithFaces {
470 return Err(MapError::Unsupported(
471 "terrain_scope is qualified for plain sloping only".into(),
472 ));
473 }
474 }
475 let started = Instant::now();
476 let [col, row, width, height] = req.window;
477 let count = width
478 .checked_mul(height)
479 .filter(|&n| n > 0 && n <= model::MAX_PROFILE_PIXELS as usize)
480 .ok_or_else(|| MapError::Budget("sloping window empty or too large".into()))?;
481 if req.schema != SCHEMA
482 || !(1..=64).contains(&req.threads)
483 || !(1..=64).contains(&req.max_inflight_profiles)
484 || req.epsg == 0
485 || req.mesh_files.is_empty()
486 || req.owner_cells.is_empty()
487 || req.owner_cells.len() >= u32::MAX as usize
488 || req.profiles.is_empty()
489 {
490 return Err(MapError::Invalid("invalid sloping request".into()));
491 }
492 if let Some(files) = &req.velocity_normals_files
493 && files.len() != req.mesh_files.len()
494 {
495 return Err(MapError::Invalid(
496 "velocity_normals_files needs one file per mesh".into(),
497 ));
498 }
499 if req.profiles.iter().any(|p| {
500 p.face_velocity_files.is_some() != req.velocity_normals_files.is_some()
501 || p.face_velocity_files
502 .as_ref()
503 .is_some_and(|f| f.len() != req.mesh_files.len())
504 }) {
505 return Err(MapError::Invalid(
506 "velocity_normals_files and every profile's face_velocity_files (one per mesh) go together".into(),
507 ));
508 }
509 let mut ids = HashSet::new();
510 for profile in &req.profiles {
511 model::validate_id(&profile.id)?;
512 if !ids.insert(model::id_key(&profile.id))
513 || profile.values_files.len() != req.mesh_files.len()
514 {
515 return Err(MapError::Invalid(
516 "duplicate sloping profile id or mesh value count mismatch".into(),
517 ));
518 }
519 }
520 let geometry_bytes = req
521 .mesh_files
522 .iter()
523 .chain(std::iter::once(&req.geometry_file))
524 .chain(req.aoi_file.iter())
525 .try_fold(0_u64, |sum, p| {
526 sum.checked_add(fs::metadata(model::resolve(base, p))?.len())
527 .ok_or_else(|| MapError::Budget("sloping geometry size overflow".into()))
528 })?;
529 let binary_bytes = req.mesh_files.iter().try_fold(0_u64, |sum, name| {
530 sum.checked_add(crate::sloping_binary::input_bytes(&model::resolve(
531 base, name,
532 ))?)
533 .ok_or_else(|| MapError::Budget("sloping binary geometry size overflow".into()))
534 })?;
535 let geometry_bytes = geometry_bytes
536 .checked_add(binary_bytes)
537 .ok_or_else(|| MapError::Budget("sloping geometry size overflow".into()))?;
538 let cache = 64_u64 << 20;
539 let persistent_bytes = if session.is_some() {
542 req.window_cache_bytes
543 } else {
544 0
545 };
546 let hydraulic_bytes = if session.is_some() {
547 req.hydraulic_cache_bytes
548 } else {
549 0
550 };
551 if let Some(s) = session.as_mut() {
552 s.windows.resize(persistent_bytes);
553 s.hydraulics.lock().unwrap().resize(hydraulic_bytes);
554 }
555 let aoi_scratch = crate::aoi::scratch_bytes(
556 req.aoi_file
557 .as_ref()
558 .map(|p| model::resolve(base, p))
559 .as_deref(),
560 req.threads,
561 )?;
562 let planned = (count as u64)
563 .checked_mul(
564 if req.options.render_mode == sloping::RenderMode::WithFaces {
565 56
566 } else {
567 40
568 },
569 )
570 .and_then(|v| {
571 geometry_bytes
572 .checked_mul(
573 if req.options.render_mode == sloping::RenderMode::WithFaces {
574 20
575 } else {
576 16
577 },
578 )
579 .and_then(|g| v.checked_add(g))
580 })
581 .and_then(|v| v.checked_add(cache))
582 .and_then(|v| v.checked_add(persistent_bytes))
583 .and_then(|v| v.checked_add(hydraulic_bytes))
584 .and_then(|v| v.checked_add(aoi_scratch))
585 .and_then(|v| {
587 v.checked_add(if req.terrain_scope.is_some() {
588 (width as u64 + 2) * (height as u64 + 2)
589 } else {
590 0
591 })
592 })
593 .ok_or_else(|| MapError::Budget("sloping array plan overflow".into()))?;
594 if planned > req.max_working_bytes {
595 return Err(MapError::Budget("sloping array plan exceeds budget".into()));
596 }
597 let key = if session.is_some() {
598 geometry_key(req, base)?
599 } else {
600 String::new()
601 };
602 if let Some(s) = session.as_mut()
604 && s.geometry.as_ref().is_some_and(|(old, _)| old != &key)
605 {
606 s.geometry = None;
607 s.hydraulics.lock().unwrap().clear();
608 s.windows.clear();
609 }
610 let terrain_path = model::resolve(base, &req.terrain);
611 let terrain_key = if session.is_some() {
612 format!(
613 "{}:{}:{}",
614 identity(&terrain_path)?,
615 req.epsg,
616 serde_json::to_string(&req.crs_identity)
617 .map_err(|e| MapError::Invalid(e.to_string()))?
618 )
619 } else {
620 String::new()
621 };
622 let terrain_reused = session.as_ref().is_some_and(|s| {
623 s.terrain
624 .as_ref()
625 .is_some_and(|(key, _)| key == &terrain_key)
626 });
627 if !terrain_reused && let Some(s) = session.as_mut() {
628 s.terrain = None;
629 s.windows.clear();
630 }
631 let terrain = if terrain_reused {
632 session
633 .as_ref()
634 .unwrap()
635 .terrain
636 .as_ref()
637 .unwrap()
638 .1
639 .clone()
640 } else {
641 Arc::new(NativeTerrain::open(
642 model::resolve(base, &req.terrain)
643 .to_str()
644 .ok_or_else(|| MapError::Invalid("terrain path is not UTF8".into()))?,
645 req.epsg,
646 &req.crs_identity,
647 QueryConvention::PixelCenter,
648 NodataPolicy::Dry,
649 cache,
650 )?)
651 };
652 if let Some(s) = session.as_mut() {
653 s.terrain = Some((terrain_key.clone(), terrain.clone()));
654 }
655 if terrain.nodata.is_none()
656 || terrain.gt != req.terrain_transform
657 || (terrain.gt[1] + terrain.gt[5]).abs() > terrain.gt[1].abs() * 1e-12
658 || col.checked_add(width).is_none_or(|v| v > terrain.width)
659 || row.checked_add(height).is_none_or(|v| v > terrain.height)
660 || terrain.max_block_bytes() > cache
661 {
662 return Err(MapError::Unsupported(
663 "sloping terrain identity, square grid, block or window mismatch".into(),
664 ));
665 }
666 let admitted = planned
667 .checked_add(
668 if req.options.render_mode == sloping::RenderMode::WithFaces {
669 (terrain.block.0 as u64)
670 .checked_mul(terrain.block.1 as u64)
671 .and_then(|v| v.checked_mul(req.threads as u64))
672 .and_then(|v| v.checked_mul(64))
673 .ok_or_else(|| {
674 MapError::Budget("native donation replay tile plan overflow".into())
675 })?
676 } else {
677 0
678 },
679 )
680 .ok_or_else(|| MapError::Budget("native donation replay plan overflow".into()))?
681 .checked_add(
682 terrain
683 .decoder_transient_bytes_per_load()
684 .checked_mul(req.threads as u64)
685 .ok_or_else(|| MapError::Budget("sloping parallel decoder plan overflow".into()))?,
686 )
687 .ok_or_else(|| MapError::Budget("sloping decoder plan overflow".into()))?;
688 let additional = (count as u64)
692 .checked_mul(
693 if req.options.render_mode == sloping::RenderMode::WithFaces {
694 32
695 } else {
696 16
697 },
698 )
699 .and_then(|v| {
700 geometry_bytes
701 .checked_mul(
702 if req.options.render_mode == sloping::RenderMode::WithFaces {
703 20
704 } else {
705 16
706 },
707 )
708 .and_then(|g| v.checked_add(g))
709 })
710 .ok_or_else(|| MapError::Budget("sloping concurrent profile plan overflow".into()))?;
711 let (effective, admitted) = admit_profiles(
712 req.max_inflight_profiles,
713 req.threads,
714 req.profiles.len(),
715 req.max_working_bytes,
716 admitted,
717 additional,
718 )?;
719 fs::create_dir(out)?;
720 let mut report = json!({"schema":"ras-commander-sloping-report/1","status":"in_progress", "method":selected_method(req.options),
721 "threads":req.threads,"window":req.window,"planned_bytes":admitted,
722 "requested_inflight_profiles":req.max_inflight_profiles,"effective_inflight_profiles":effective,
723 "additional_profile_planned_bytes":additional,"profiles":[]});
724 sink::write_json(&out.join("report.json"), &report)?;
725 let execution = std::panic::catch_unwind(std::panic::AssertUnwindSafe(
726 || -> Result<(), MapError> {
727 if let Some(s) = session.as_mut()
728 && s.pool
729 .as_ref()
730 .is_some_and(|(threads, _)| *threads != req.threads)
731 {
732 s.pool = None;
733 }
734 let pool = if let Some((threads, pool)) = session
735 .as_ref()
736 .and_then(|s| s.pool.as_ref())
737 .filter(|(threads, _)| *threads == req.threads)
738 {
739 let _ = threads;
740 pool.clone()
741 } else {
742 Arc::new(
743 rayon::ThreadPoolBuilder::new()
744 .num_threads(req.threads)
745 .build()
746 .map_err(|e| MapError::Internal(e.to_string()))?,
747 )
748 };
749 if let Some(s) = session.as_mut() {
750 s.pool = Some((req.threads, pool.clone()));
751 }
752 let tick = Instant::now();
753 let geometry_reused = session
754 .as_ref()
755 .is_some_and(|s| s.geometry.as_ref().is_some_and(|(old, _)| old == &key));
756 if !geometry_reused && let Some(s) = session.as_mut() {
757 s.geometry = None;
758 s.hydraulics.lock().unwrap().clear();
759 s.windows.clear();
760 }
761 let geometry = if geometry_reused {
762 session
763 .as_ref()
764 .unwrap()
765 .geometry
766 .as_ref()
767 .unwrap()
768 .1
769 .clone()
770 } else {
771 let mut meshes = Vec::new();
772 let mut stencils = Vec::new();
773 for filename in &req.mesh_files {
774 let source = crate::sloping_binary::load(&model::resolve(base, filename))?;
775 let (mesh, stencil) = source.prepare()?;
776 meshes.push(mesh);
777 stencils.push(stencil);
778 }
779 Arc::new((meshes, stencils))
780 };
781 if let Some(s) = session.as_mut() {
782 s.geometry = Some((key.clone(), geometry.clone()));
783 }
784 let (meshes, stencils) = &*geometry;
785 let velocity_normals = req
786 .velocity_normals_files
787 .as_ref()
788 .map(|files| -> Result<Vec<_>, MapError> {
789 files
790 .iter()
791 .zip(meshes)
792 .map(|(file, mesh)| {
793 crate::velocity::read_normals(
794 &model::resolve(base, file),
795 mesh.faces.len(),
796 )
797 })
798 .collect()
799 })
800 .transpose()?;
801 if velocity_normals.is_some() {
802 report["velocity"] = json!({"method": crate::velocity::method_record()});
803 }
804 report["session"] = json!({"terrain_reused":terrain_reused,"geometry_reused":geometry_reused,"retained_window_budget_bytes":persistent_bytes});
805 for &[mesh, cell] in &req.owner_cells {
806 if mesh >= meshes.len()
807 || cell >= meshes[mesh].cells.len()
808 || meshes[mesh].cells[cell].virtual_cell
809 {
810 return Err(MapError::Invalid(
811 "sloping ownership requires physical mesh cells".into(),
812 ));
813 }
814 }
815 report["geometry_prepare_seconds"] = tick.elapsed().as_secs_f64().into();
816 report["ownership_prepare_seconds"] = 0.0.into();
817 report["terrain_prepare_seconds"] = 0.0.into();
818 let transform = [
819 terrain.gt[0] + col as f64 * terrain.gt[1],
820 terrain.gt[1],
821 0.,
822 terrain.gt[3] + row as f64 * terrain.gt[5],
823 0.,
824 terrain.gt[5],
825 ];
826 let grid = terrio::ops::rasterize::Grid {
827 geo_transform: transform,
828 cols: width,
829 rows: height,
830 crs: Default::default(),
831 };
832 let tick = Instant::now();
833 let window_key = if session.is_some() {
834 (
835 terrain_key.clone(),
836 key.clone(),
837 req.window,
838 req.owner_cells.clone(),
839 fs::read(model::resolve(base, &req.geometry_file))?,
840 req.aoi_file
841 .as_ref()
842 .map(|p| fs::read(model::resolve(base, p)))
843 .transpose()?,
844 )
845 } else {
846 (
847 String::new(),
848 String::new(),
849 req.window,
850 Vec::new(),
851 Vec::new(),
852 None,
853 )
854 };
855 let cached_window = session.as_mut().and_then(|s| s.windows.get(&window_key));
856 let window_reused = cached_window.is_some();
857 let mut data = if let Some(data) = cached_window {
858 data
859 } else {
860 let cells = crate::ownership::load(
861 &model::resolve(base, &req.geometry_file),
862 req.owner_cells.len(),
863 false,
864 )?;
865 let mut owner = pool.install(|| crate::ownership::burn(cells, grid.clone()))?;
866 if let Some(filename) = &req.aoi_file {
867 let mask = pool.install(|| {
868 crate::aoi::mask(&model::resolve(base, filename), terrain.gt, req.window)
869 })?;
870 pool.install(|| {
871 owner
872 .par_iter_mut()
873 .zip(mask.par_iter())
874 .for_each(|(label, &m)| {
875 if m == 0 {
876 *label = 0;
877 }
878 })
879 });
880 }
881 report["ownership_prepare_seconds"] = tick.elapsed().as_secs_f64().into();
882 let tick = Instant::now();
883 let mut ground = vec![sloping::NODATA; count];
884 let (bw, bh) = terrain.block;
885 let stripe_pixels = width
889 .checked_mul(bh)
890 .ok_or_else(|| MapError::Budget("sloping terrain stripe overflow".into()))?;
891 pool.install(|| {
892 ground
893 .par_chunks_mut(stripe_pixels)
894 .enumerate()
895 .try_for_each(|(stripe, destination)| -> Result<(), MapError> {
896 let first = row + stripe * bh;
897 let last = first + destination.len() / width;
898 for by in first / bh..=(last - 1) / bh {
899 for bx in col / bw..=(col + width - 1) / bw {
900 let win = terrain.block_window(bx, by);
901 let r0 = first.max(win.y);
902 let r1 = last.min(win.y + win.height);
903 let c0 = col.max(win.x);
904 let c1 = (col + width).min(win.x + win.width);
905 if !(r0..r1).any(|r| {
906 owner[(r - row) * width + c0 - col
907 ..(r - row) * width + c1 - col]
908 .iter()
909 .any(|&label| label != 0)
910 }) {
911 continue;
912 }
913 let block = terrain.get_block(bx, by)?;
914 for r in r0..r1 {
915 for c in c0..c1 {
916 if owner[(r - row) * width + c - col] == 0 {
917 continue;
918 }
919 if let Some(value) = terrain.classify(
920 block[(r - win.y) * win.width + c - win.x],
921 ) {
922 if !value.is_finite() {
923 return Err(MapError::Terrain(
924 "nonfinite sloping terrain".into(),
925 ));
926 }
927 destination[(r - first) * width + c - col] =
928 value as f32;
929 }
930 }
931 }
932 }
933 }
934 Ok(())
935 })
936 })?;
937 report["terrain_prepare_seconds"] = tick.elapsed().as_secs_f64().into();
938 Arc::new(WindowData {
939 owner,
940 ground,
941 weights: None,
942 })
943 };
944 let window_base_bytes = count as u64 * 8
945 + (2 * std::mem::size_of::<(WindowKey, Arc<WindowData>, u64)>()
946 + std::mem::size_of::<WindowData>()
947 + 2 * std::mem::size_of::<usize>()) as u64
948 + window_key.0.len() as u64
949 + window_key.1.len() as u64
950 + (window_key.3.len() * std::mem::size_of::<[usize; 2]>()) as u64
951 + window_key.4.len() as u64
952 + window_key.5.as_ref().map_or(0, |v| v.len()) as u64;
953 let weight_budget = persistent_bytes
954 .saturating_sub(window_base_bytes)
955 .min(32 << 20);
956 let (bw, bh) = terrain.block;
957 if session.is_some()
958 && !window_reused
959 && weight_budget > 0
960 && req.options.render_mode != sloping::RenderMode::WithFaces
961 {
962 let window = sloping::RasterWindow {
963 source_transform: terrain.gt,
964 source_tile: [bw, bh],
965 width,
966 height,
967 transform,
968 owner: &data.owner,
969 terrain: &data.ground,
970 owner_cells: &req.owner_cells,
971 };
972 let weights = sloping::PixelWeights::prepare(
973 stencils,
974 &window,
975 weight_budget as usize,
976 &pool,
977 )?;
978 Arc::get_mut(&mut data).unwrap().weights = Some(weights);
979 }
980 let weight_bytes = data.weights.as_ref().map_or(0, |w| w.bytes());
981 report["session"]["window_reused"] = window_reused.into();
982 report["session"]["pixel_weight_cache_bytes"] = weight_bytes.into();
983 let window_bytes = window_base_bytes + weight_bytes as u64;
984 if let Some(s) = session.as_mut() {
985 s.windows.insert(window_key, data.clone(), window_bytes);
986 report["session"]["retained_window_bytes"] = s.windows.bytes().into();
987 report["session"]["retained_window_count"] = s.windows.len().into();
988 }
989 let window = sloping::RasterWindow {
990 source_transform: terrain.gt,
991 source_tile: [bw, bh],
992 width,
993 height,
994 transform,
995 owner: &data.owner,
996 terrain: &data.ground,
997 owner_cells: &req.owner_cells,
998 };
999 let foreign = match &req.terrain_scope {
1000 Some(scope) => {
1001 let scoped = crate::terrain_scope::check(
1002 scope,
1003 base,
1004 &terrain,
1005 req.window,
1006 &data.owner,
1007 )?;
1008 report["terrain_scope"] = scoped.record;
1009 scoped.foreign
1010 }
1011 None => Vec::new(),
1012 };
1013 for batch in req.profiles.chunks(effective) {
1014 let results = map_batch(batch, &pool, |profile| {
1015 render(
1016 profile,
1017 base,
1018 meshes,
1019 stencils,
1020 &window,
1021 &pool,
1022 req.options,
1023 data.weights.as_ref(),
1024 velocity_normals.as_deref(),
1025 &foreign,
1026 session
1027 .as_ref()
1028 .filter(|_| hydraulic_bytes > 0)
1029 .map(|s| &s.hydraulics),
1030 )
1031 })?;
1032 for (profile, rendered) in batch.iter().zip(results) {
1035 let RenderedProfile {
1036 depth,
1037 wse,
1038 velocity,
1039 hydraulic_seconds,
1040 render_seconds,
1041 hydraulic_reused,
1042 finite_pixels,
1043 flooded_pixels,
1044 } = rendered;
1045 let dir = out.join(&profile.id);
1046 fs::create_dir(&dir)?;
1047 if req.stream_outputs {
1048 let mut stream =
1049 std::io::BufWriter::with_capacity(65536, std::io::stdout().lock());
1050 let mut header =
1051 json!({"id":profile.id,"status":"complete","pixels":count});
1052 if velocity.is_some() {
1053 header["velocity"] = true.into();
1054 }
1055 if req.geotiff_outputs {
1056 let tick = Instant::now();
1057 header["products"] = "native-geotiff/1".into();
1058 header["depth_finite_pixels"] = finite_pixels.into();
1059 header["flooded_pixels"] = flooded_pixels.into();
1060 crate::geotiff::write(
1061 &dir.join("depth.tif"),
1062 depth,
1063 width,
1064 height,
1065 transform,
1066 req.epsg,
1067 )?;
1068 crate::geotiff::write(
1069 &dir.join("wse.tif"),
1070 wse,
1071 width,
1072 height,
1073 transform,
1074 req.epsg,
1075 )?;
1076 header["write_seconds"] = tick.elapsed().as_secs_f64().into();
1077 writeln!(stream, "{header}")?;
1078 stream.flush()?;
1079 } else {
1080 writeln!(stream, "{header}")?;
1081 crate::horizontal::write_f32(&mut stream, &depth)?;
1082 crate::horizontal::write_f32(&mut stream, &wse)?;
1083 if let Some(velocity) = &velocity {
1084 crate::horizontal::write_f32(&mut stream, velocity)?;
1085 }
1086 stream.flush()?;
1087 }
1088 } else {
1089 crate::horizontal::write_values(&dir.join("depth.bin"), &depth)?;
1090 crate::horizontal::write_values(&dir.join("wse.bin"), &wse)?;
1091 if let Some(velocity) = &velocity {
1092 crate::horizontal::write_values(&dir.join("velocity.bin"), velocity)?;
1093 }
1094 }
1095 let unit = json!({"id":profile.id,"status":"complete","hydraulic_seconds":hydraulic_seconds,"render_seconds":render_seconds,"hydraulic_reused":hydraulic_reused});
1096 sink::write_json(&dir.join("unit.json"), &unit)?;
1097 report["profiles"].as_array_mut().unwrap().push(unit);
1098 sink::write_json(&out.join("report.json"), &report)?;
1099 }
1100 }
1101 if let Some(s) = session.as_ref() {
1102 let cache = s.hydraulics.lock().unwrap();
1103 report["session"]["retained_hydraulic_bytes"] = cache.bytes().into();
1104 report["session"]["retained_hydraulic_count"] = cache.len().into();
1105 }
1106 Ok(())
1107 },
1108 ));
1109 match execution {
1110 Ok(Ok(())) => report["status"] = "complete".into(),
1111 other => {
1112 report["status"] = "aborted".into();
1113 report["error"] = match other {
1114 Ok(Err(e)) => e.to_string(),
1115 _ => "sloping internal panic".into(),
1116 }
1117 .into();
1118 for profile in &req.profiles {
1119 let dir = out.join(&profile.id);
1120 for filename in [
1121 "depth.bin",
1122 "wse.bin",
1123 "depth.bin.partial",
1124 "wse.bin.partial",
1125 "depth.tif",
1126 "wse.tif",
1127 ] {
1128 let _ = fs::remove_file(dir.join(filename));
1129 }
1130 if dir.is_dir() {
1131 sink::write_json(
1132 &dir.join("unit.json"),
1133 &json!({"id":profile.id,"status":"failed","error":"session aborted"}),
1134 )?;
1135 }
1136 }
1137 report["profiles"] = req
1138 .profiles
1139 .iter()
1140 .map(|p| json!({"id":p.id,"status":"failed","error":"session aborted"}))
1141 .collect::<Vec<_>>()
1142 .into();
1143 }
1144 }
1145 report["total_seconds"] = started.elapsed().as_secs_f64().into();
1146 report["terrain_cache"] = json!(terrain.stats());
1147 sink::write_json(&out.join("report.json"), &report)?;
1148 Ok(report)
1149}
1150
1151#[cfg(test)]
1152mod tests {
1153 use super::*;
1154 use std::sync::{
1155 Barrier,
1156 atomic::{AtomicUsize, Ordering},
1157 };
1158
1159 #[test]
1160 fn admission_caps_threads_profiles_memory_and_overflow() {
1161 assert_eq!(admit_profiles(8, 4, 6, 1000, 100, 50).unwrap(), (4, 250));
1162 assert_eq!(admit_profiles(8, 4, 2, 1000, 100, 50).unwrap(), (2, 150));
1163 assert_eq!(admit_profiles(8, 4, 6, 149, 100, 50).unwrap(), (1, 100));
1164 assert!(admit_profiles(8, 4, 6, 99, 100, 50).is_err());
1165 assert_eq!(
1166 admit_profiles(64, 64, 64, u64::MAX, u64::MAX - 1, u64::MAX).unwrap(),
1167 (1, u64::MAX - 1)
1168 );
1169 assert!(admit_profiles(0, 1, 1, 100, 50, 10).is_err());
1170 }
1171
1172 #[test]
1173 fn batches_really_overlap_without_exceeding_bound_and_keep_order() {
1174 let pool = rayon::ThreadPoolBuilder::new()
1175 .num_threads(2)
1176 .build()
1177 .unwrap();
1178 let active = AtomicUsize::new(0);
1179 let peak = AtomicUsize::new(0);
1180 let barrier = Barrier::new(2);
1181 let mut output = Vec::new();
1182 for batch in [0, 1, 2, 3].chunks(2) {
1183 output.extend(
1184 map_batch(batch, &pool, |&i| {
1185 let n = active.fetch_add(1, Ordering::SeqCst) + 1;
1186 peak.fetch_max(n, Ordering::SeqCst);
1187 barrier.wait();
1188 active.fetch_sub(1, Ordering::SeqCst);
1189 Ok(i)
1190 })
1191 .unwrap(),
1192 );
1193 }
1194 assert_eq!(peak.load(Ordering::SeqCst), 2);
1195 assert_eq!(output, vec![0, 1, 2, 3]);
1196 assert!(
1197 map_batch(&[0, 1], &pool, |_| -> Result<(), MapError> {
1198 Err(MapError::Invalid("failed profile".into()))
1199 })
1200 .is_err()
1201 );
1202 }
1203}