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 expr = parser.parse()?;
57
58        // Optimize expression
59        let expr = Optimizer::optimize(expr);
60
61        // Evaluate
62        let evaluator = Evaluator::new(bands);
63        let mut result = RasterBuffer::zeros(width, height, bands[0].data_type());
64
65        for y in 0..height {
66            for x in 0..width {
67                let value = evaluator.eval_pixel(&expr, x, y)?;
68                result
69                    .set_pixel(x, y, value)
70                    .map_err(AlgorithmError::Core)?;
71            }
72        }
73
74        Ok(result)
75    }
76
77    /// Evaluates a raster expression in parallel using rayon
78    ///
79    /// This method processes rows in parallel for improved performance on multi-core systems.
80    /// Falls back to sequential evaluation if the parallel feature is not enabled.
81    ///
82    /// # Arguments
83    ///
84    /// * `expression` - The expression to evaluate (e.g., "(B1 - B2) / (B1 + B2)")
85    /// * `bands` - Input bands (B1, B2, etc.)
86    ///
87    /// # Examples
88    ///
89    /// NDVI: `"(B1 - B2) / (B1 + B2)"`
90    /// Conditional: `"if B1 > 100 then B1 * 2 else B1"`
91    /// Math: `"sqrt(B1 ^ 2 + B2 ^ 2)"`
92    ///
93    /// # Errors
94    ///
95    /// Returns an error if the expression is invalid or evaluation fails
96    #[cfg(feature = "parallel")]
97    pub fn evaluate_parallel(expression: &str, bands: &[RasterBuffer]) -> Result<RasterBuffer> {
98        if bands.is_empty() {
99            return Err(AlgorithmError::EmptyInput {
100                operation: "evaluate_parallel",
101            });
102        }
103
104        // Check all bands have same dimensions
105        let width = bands[0].width();
106        let height = bands[0].height();
107        for band in bands.iter().skip(1) {
108            if band.width() != width || band.height() != height {
109                return Err(AlgorithmError::InvalidDimensions {
110                    message: "All bands must have same dimensions",
111                    actual: band.width() as usize,
112                    expected: width as usize,
113                });
114            }
115        }
116
117        // Tokenize
118        let mut lexer = Lexer::new(expression);
119        let tokens = lexer.tokenize()?;
120
121        // Parse
122        let mut parser = Parser::new(tokens);
123        let expr = parser.parse()?;
124
125        // Optimize expression
126        let expr = Optimizer::optimize(expr);
127
128        // Create evaluator
129        let evaluator = Evaluator::new(bands);
130
131        // Create result buffer
132        let mut result = RasterBuffer::zeros(width, height, bands[0].data_type());
133
134        // Process rows in parallel
135        let row_data: Result<Vec<Vec<f64>>> = (0..height)
136            .into_par_iter()
137            .map(|y| {
138                let mut row = Vec::with_capacity(width as usize);
139                for x in 0..width {
140                    let value = evaluator.eval_pixel(&expr, x, y)?;
141                    row.push(value);
142                }
143                Ok(row)
144            })
145            .collect();
146
147        let row_data = row_data?;
148
149        // Write results back to buffer
150        for (y, row) in row_data.iter().enumerate() {
151            for (x, &value) in row.iter().enumerate() {
152                result
153                    .set_pixel(x as u64, y as u64, value)
154                    .map_err(AlgorithmError::Core)?;
155            }
156        }
157
158        Ok(result)
159    }
160
161    /// Applies a binary operation to two rasters (legacy API)
162    pub fn apply_binary(
163        a: &RasterBuffer,
164        b: &RasterBuffer,
165        op: RasterExpression,
166    ) -> Result<RasterBuffer> {
167        if a.width() != b.width() || a.height() != b.height() {
168            return Err(AlgorithmError::InvalidDimensions {
169                message: "Rasters must have same dimensions",
170                actual: a.width() as usize,
171                expected: b.width() as usize,
172            });
173        }
174
175        let mut result = RasterBuffer::zeros(a.width(), a.height(), a.data_type());
176
177        for y in 0..a.height() {
178            for x in 0..a.width() {
179                let val_a = a.get_pixel(x, y).map_err(AlgorithmError::Core)?;
180                let val_b = b.get_pixel(x, y).map_err(AlgorithmError::Core)?;
181
182                let val = match op {
183                    RasterExpression::Add => val_a + val_b,
184                    RasterExpression::Subtract => val_a - val_b,
185                    RasterExpression::Multiply => val_a * val_b,
186                    RasterExpression::Divide => {
187                        if val_b.abs() < f64::EPSILON {
188                            f64::NAN
189                        } else {
190                            val_a / val_b
191                        }
192                    }
193                    RasterExpression::Max => val_a.max(val_b),
194                    RasterExpression::Min => val_a.min(val_b),
195                };
196
197                result.set_pixel(x, y, val).map_err(AlgorithmError::Core)?;
198            }
199        }
200
201        Ok(result)
202    }
203
204    /// Applies a unary function to a raster (legacy API)
205    pub fn apply_unary<F>(src: &RasterBuffer, func: F) -> Result<RasterBuffer>
206    where
207        F: Fn(f64) -> f64,
208    {
209        let mut result = RasterBuffer::zeros(src.width(), src.height(), src.data_type());
210
211        for y in 0..src.height() {
212            for x in 0..src.width() {
213                let val = src.get_pixel(x, y).map_err(AlgorithmError::Core)?;
214                let new_val = func(val);
215                result
216                    .set_pixel(x, y, new_val)
217                    .map_err(AlgorithmError::Core)?;
218            }
219        }
220
221        Ok(result)
222    }
223}
224
225/// Raster expression operations (legacy API)
226#[derive(Debug, Clone, Copy, PartialEq)]
227pub enum RasterExpression {
228    /// Add two rasters
229    Add,
230    /// Subtract rasters
231    Subtract,
232    /// Multiply rasters
233    Multiply,
234    /// Divide rasters
235    Divide,
236    /// Maximum of two rasters
237    Max,
238    /// Minimum of two rasters
239    Min,
240}
oxigdal_algorithms/raster/calculator/ops.rs

oxigdal_algorithms/raster/calculator/
ops.rs
