oxigdal_algorithms/raster/calculator/
ops.rs

1//! High-level raster calculator API: `RasterCalculator` and `RasterExpression`
2
3use super::{evaluator::Evaluator, lexer::Lexer, optimizer::Optimizer, parser::Parser};
4use crate::error::{AlgorithmError, Result};
5use oxigdal_core::buffer::RasterBuffer;
6
7#[cfg(feature = "parallel")]
8use rayon::prelude::*;
9
10/// Raster calculator for map algebra with expression parsing
11pub struct RasterCalculator;
12
13impl RasterCalculator {
14    /// Evaluates a raster expression on one or more bands
15    ///
16    /// # Arguments
17    ///
18    /// * `expression` - The expression to evaluate (e.g., "(B1 - B2) / (B1 + B2)")
19    /// * `bands` - Input bands (B1, B2, etc.)
20    ///
21    /// # Examples
22    ///
23    /// NDVI: `"(B1 - B2) / (B1 + B2)"`
24    /// Conditional: `"if B1 > 100 then B1 * 2 else B1"`
25    /// Math: `"sqrt(B1 ^ 2 + B2 ^ 2)"`
26    ///
27    /// # Errors
28    ///
29    /// Returns an error if the expression is invalid or evaluation fails
30    pub fn evaluate(expression: &str, bands: &[RasterBuffer]) -> Result<RasterBuffer> {
31        if bands.is_empty() {
32            return Err(AlgorithmError::EmptyInput {
33                operation: "evaluate",
34            });
35        }
36
37        // Check all bands have same dimensions
38        let width = bands[0].width();
39        let height = bands[0].height();
40        for (_i, band) in bands.iter().enumerate().skip(1) {
41            if band.width() != width || band.height() != height {
42                return Err(AlgorithmError::InvalidDimensions {
43                    message: "All bands must have same dimensions",
44                    actual: band.width() as usize,
45                    expected: width as usize,
46                });
47            }
48        }
49
50        // Tokenize
51        let mut lexer = Lexer::new(expression);
52        let tokens = lexer.tokenize()?;
53
54        // Parse
55        let mut parser = Parser::new(tokens);
56        let raw_expr = parser.parse()?;
57
58        // Optimize — returns (rewritten_expr, cse_slots)
59        let (expr, cache_slots) = Optimizer::optimize(raw_expr);
60
61        // Build evaluator with CSE slot expressions
62        let evaluator = Evaluator::new(bands, &cache_slots);
63        let slot_count = cache_slots.len();
64
65        let mut result = RasterBuffer::zeros(width, height, bands[0].data_type());
66
67        for y in 0..height {
68            for x in 0..width {
69                // One fresh pixel-cache per pixel — reset to all-None
70                let mut pixel_cache: Vec<Option<f64>> = vec![None; slot_count];
71                let value = evaluator.eval_pixel(&expr, x, y, &mut pixel_cache)?;
72                result
73                    .set_pixel(x, y, value)
74                    .map_err(AlgorithmError::Core)?;
75            }
76        }
77
78        Ok(result)
79    }
80
81    /// Evaluates a raster expression in parallel using rayon
82    ///
83    /// This method processes rows in parallel for improved performance on multi-core systems.
84    /// Falls back to sequential evaluation if the parallel feature is not enabled.
85    ///
86    /// # Arguments
87    ///
88    /// * `expression` - The expression to evaluate (e.g., "(B1 - B2) / (B1 + B2)")
89    /// * `bands` - Input bands (B1, B2, etc.)
90    ///
91    /// # Examples
92    ///
93    /// NDVI: `"(B1 - B2) / (B1 + B2)"`
94    /// Conditional: `"if B1 > 100 then B1 * 2 else B1"`
95    /// Math: `"sqrt(B1 ^ 2 + B2 ^ 2)"`
96    ///
97    /// # Errors
98    ///
99    /// Returns an error if the expression is invalid or evaluation fails
100    #[cfg(feature = "parallel")]
101    pub fn evaluate_parallel(expression: &str, bands: &[RasterBuffer]) -> Result<RasterBuffer> {
102        if bands.is_empty() {
103            return Err(AlgorithmError::EmptyInput {
104                operation: "evaluate_parallel",
105            });
106        }
107
108        // Check all bands have same dimensions
109        let width = bands[0].width();
110        let height = bands[0].height();
111        for band in bands.iter().skip(1) {
112            if band.width() != width || band.height() != height {
113                return Err(AlgorithmError::InvalidDimensions {
114                    message: "All bands must have same dimensions",
115                    actual: band.width() as usize,
116                    expected: width as usize,
117                });
118            }
119        }
120
121        // Tokenize
122        let mut lexer = Lexer::new(expression);
123        let tokens = lexer.tokenize()?;
124
125        // Parse
126        let mut parser = Parser::new(tokens);
127        let raw_expr = parser.parse()?;
128
129        // Optimize — returns (rewritten_expr, cse_slots)
130        let (expr, cache_slots) = Optimizer::optimize(raw_expr);
131
132        // Create evaluator
133        let evaluator = Evaluator::new(bands, &cache_slots);
134        let slot_count = cache_slots.len();
135
136        // Create result buffer
137        let mut result = RasterBuffer::zeros(width, height, bands[0].data_type());
138
139        // Process rows in parallel — each row gets its own pixel_cache per pixel
140        let row_data: Result<Vec<Vec<f64>>> = (0..height)
141            .into_par_iter()
142            .map(|y| {
143                let mut row = Vec::with_capacity(width as usize);
144                for x in 0..width {
145                    let mut pixel_cache: Vec<Option<f64>> = vec![None; slot_count];
146                    let value = evaluator.eval_pixel(&expr, x, y, &mut pixel_cache)?;
147                    row.push(value);
148                }
149                Ok(row)
150            })
151            .collect();
152
153        let row_data = row_data?;
154
155        // Write results back to buffer
156        for (y, row) in row_data.iter().enumerate() {
157            for (x, &value) in row.iter().enumerate() {
158                result
159                    .set_pixel(x as u64, y as u64, value)
160                    .map_err(AlgorithmError::Core)?;
161            }
162        }
163
164        Ok(result)
165    }
166
167    /// Applies a binary operation to two rasters (legacy API)
168    pub fn apply_binary(
169        a: &RasterBuffer,
170        b: &RasterBuffer,
171        op: RasterExpression,
172    ) -> Result<RasterBuffer> {
173        if a.width() != b.width() || a.height() != b.height() {
174            return Err(AlgorithmError::InvalidDimensions {
175                message: "Rasters must have same dimensions",
176                actual: a.width() as usize,
177                expected: b.width() as usize,
178            });
179        }
180
181        let mut result = RasterBuffer::zeros(a.width(), a.height(), a.data_type());
182
183        for y in 0..a.height() {
184            for x in 0..a.width() {
185                let val_a = a.get_pixel(x, y).map_err(AlgorithmError::Core)?;
186                let val_b = b.get_pixel(x, y).map_err(AlgorithmError::Core)?;
187
188                let val = match op {
189                    RasterExpression::Add => val_a + val_b,
190                    RasterExpression::Subtract => val_a - val_b,
191                    RasterExpression::Multiply => val_a * val_b,
192                    RasterExpression::Divide => {
193                        if val_b.abs() < f64::EPSILON {
194                            f64::NAN
195                        } else {
196                            val_a / val_b
197                        }
198                    }
199                    RasterExpression::Max => val_a.max(val_b),
200                    RasterExpression::Min => val_a.min(val_b),
201                };
202
203                result.set_pixel(x, y, val).map_err(AlgorithmError::Core)?;
204            }
205        }
206
207        Ok(result)
208    }
209
210    /// Applies a unary function to a raster (legacy API)
211    pub fn apply_unary<F>(src: &RasterBuffer, func: F) -> Result<RasterBuffer>
212    where
213        F: Fn(f64) -> f64,
214    {
215        let mut result = RasterBuffer::zeros(src.width(), src.height(), src.data_type());
216
217        for y in 0..src.height() {
218            for x in 0..src.width() {
219                let val = src.get_pixel(x, y).map_err(AlgorithmError::Core)?;
220                let new_val = func(val);
221                result
222                    .set_pixel(x, y, new_val)
223                    .map_err(AlgorithmError::Core)?;
224            }
225        }
226
227        Ok(result)
228    }
229}
230
231/// Raster expression operations (legacy API)
232#[derive(Debug, Clone, Copy, PartialEq)]
233pub enum RasterExpression {
234    /// Add two rasters
235    Add,
236    /// Subtract rasters
237    Subtract,
238    /// Multiply rasters
239    Multiply,
240    /// Divide rasters
241    Divide,
242    /// Maximum of two rasters
243    Max,
244    /// Minimum of two rasters
245    Min,
246}
oxigdal_algorithms/raster/calculator/ops.rs

oxigdal_algorithms/raster/calculator/
ops.rs
