oxigdal_algorithms/raster/
calculator.rs

1//! Raster calculator (map algebra) with expression parsing
2//!
3//! This module provides a comprehensive raster calculator with support for:
4//! - Arithmetic operations: +, -, *, /, ^
5//! - Math functions: sqrt, log, exp, sin, cos, tan, abs, floor, ceil, etc.
6//! - Band algebra: (B1 - B2) / (B1 + B2) for NDVI and similar indices
7//! - Conditional operations: if/then/else
8//! - Multi-band operations
9//! - Proper NoData handling
10
11use crate::error::{AlgorithmError, Result};
12use oxigdal_core::buffer::RasterBuffer;
13
14#[cfg(not(feature = "std"))]
15use alloc::{string::String, vec::Vec};
16
17#[cfg(feature = "parallel")]
18use rayon::prelude::*;
19
20/// Token types for expression lexer
21#[derive(Debug, Clone, PartialEq)]
22enum Token {
23    /// Number literal
24    Number(f64),
25    /// Band reference (e.g., B1, B2)
26    Band(usize),
27    /// Identifier (function name or variable)
28    Ident(String),
29    /// Operators
30    Plus,
31    Minus,
32    Multiply,
33    Divide,
34    Power,
35    /// Parentheses
36    LeftParen,
37    RightParen,
38    /// Comparison operators
39    Greater,
40    Less,
41    GreaterEqual,
42    LessEqual,
43    Equal,
44    NotEqual,
45    /// Logical operators
46    And,
47    Or,
48    /// Keywords
49    If,
50    Then,
51    Else,
52    /// Comma (for function arguments)
53    Comma,
54}
55
56/// Expression AST node
57#[derive(Debug, Clone)]
58enum Expr {
59    /// Number literal
60    Number(f64),
61    /// Band reference
62    Band(usize),
63    /// Binary operation
64    BinaryOp {
65        left: Box<Expr>,
66        op: BinaryOp,
67        right: Box<Expr>,
68    },
69    /// Unary operation
70    UnaryOp { op: UnaryOp, expr: Box<Expr> },
71    /// Function call
72    Function { name: String, args: Vec<Expr> },
73    /// Conditional expression
74    Conditional {
75        condition: Box<Expr>,
76        then_expr: Box<Expr>,
77        else_expr: Box<Expr>,
78    },
79}
80
81/// Binary operators
82#[derive(Debug, Clone, Copy, PartialEq)]
83enum BinaryOp {
84    Add,
85    Subtract,
86    Multiply,
87    Divide,
88    Power,
89    Greater,
90    Less,
91    GreaterEqual,
92    LessEqual,
93    Equal,
94    NotEqual,
95    And,
96    Or,
97}
98
99/// Unary operators
100#[derive(Debug, Clone, Copy, PartialEq)]
101enum UnaryOp {
102    Negate,
103}
104
105/// Tokenizer for raster expressions
106struct Lexer {
107    input: Vec<char>,
108    pos: usize,
109}
110
111impl Lexer {
112    fn new(input: &str) -> Self {
113        Self {
114            input: input.chars().collect(),
115            pos: 0,
116        }
117    }
118
119    fn current(&self) -> Option<char> {
120        if self.pos < self.input.len() {
121            Some(self.input[self.pos])
122        } else {
123            None
124        }
125    }
126
127    fn advance(&mut self) {
128        self.pos += 1;
129    }
130
131    fn skip_whitespace(&mut self) {
132        while let Some(c) = self.current() {
133            if c.is_whitespace() {
134                self.advance();
135            } else {
136                break;
137            }
138        }
139    }
140
141    fn read_number(&mut self) -> Result<f64> {
142        let mut num_str = String::new();
143        let mut has_dot = false;
144
145        while let Some(c) = self.current() {
146            if c.is_ascii_digit() {
147                num_str.push(c);
148                self.advance();
149            } else if c == '.' && !has_dot {
150                has_dot = true;
151                num_str.push(c);
152                self.advance();
153            } else {
154                break;
155            }
156        }
157
158        num_str
159            .parse::<f64>()
160            .map_err(|_| AlgorithmError::InvalidParameter {
161                parameter: "expression",
162                message: format!("Invalid number: {num_str}"),
163            })
164    }
165
166    fn read_ident(&mut self) -> String {
167        let mut ident = String::new();
168
169        while let Some(c) = self.current() {
170            if c.is_alphanumeric() || c == '_' {
171                ident.push(c);
172                self.advance();
173            } else {
174                break;
175            }
176        }
177
178        ident
179    }
180
181    fn next_token(&mut self) -> Result<Option<Token>> {
182        self.skip_whitespace();
183
184        let c = match self.current() {
185            Some(c) => c,
186            None => return Ok(None),
187        };
188
189        let token = if c.is_ascii_digit() {
190            Token::Number(self.read_number()?)
191        } else if c.is_alphabetic() || c == 'B' {
192            let ident = self.read_ident();
193
194            // Check for band reference (B1, B2, etc.)
195            if ident.starts_with('B') && ident.len() > 1 {
196                let band_num = ident[1..].parse::<usize>().ok();
197                if let Some(num) = band_num {
198                    return Ok(Some(Token::Band(num)));
199                }
200            }
201
202            // Check for keywords
203            match ident.as_str() {
204                "if" | "IF" => Token::If,
205                "then" | "THEN" => Token::Then,
206                "else" | "ELSE" => Token::Else,
207                "and" | "AND" => Token::And,
208                "or" | "OR" => Token::Or,
209                _ => Token::Ident(ident),
210            }
211        } else {
212            self.advance();
213            match c {
214                '+' => Token::Plus,
215                '-' => Token::Minus,
216                '*' => Token::Multiply,
217                '/' => Token::Divide,
218                '^' => Token::Power,
219                '(' => Token::LeftParen,
220                ')' => Token::RightParen,
221                ',' => Token::Comma,
222                '>' => {
223                    if self.current() == Some('=') {
224                        self.advance();
225                        Token::GreaterEqual
226                    } else {
227                        Token::Greater
228                    }
229                }
230                '<' => {
231                    if self.current() == Some('=') {
232                        self.advance();
233                        Token::LessEqual
234                    } else {
235                        Token::Less
236                    }
237                }
238                '=' => {
239                    if self.current() == Some('=') {
240                        self.advance();
241                        Token::Equal
242                    } else {
243                        return Err(AlgorithmError::InvalidParameter {
244                            parameter: "expression",
245                            message: "Expected '==' for equality comparison".to_string(),
246                        });
247                    }
248                }
249                '!' => {
250                    if self.current() == Some('=') {
251                        self.advance();
252                        Token::NotEqual
253                    } else {
254                        return Err(AlgorithmError::InvalidParameter {
255                            parameter: "expression",
256                            message: "Expected '!=' for inequality comparison".to_string(),
257                        });
258                    }
259                }
260                _ => {
261                    return Err(AlgorithmError::InvalidParameter {
262                        parameter: "expression",
263                        message: format!("Unexpected character: {c}"),
264                    });
265                }
266            }
267        };
268
269        Ok(Some(token))
270    }
271
272    fn tokenize(&mut self) -> Result<Vec<Token>> {
273        let mut tokens = Vec::new();
274
275        while let Some(token) = self.next_token()? {
276            tokens.push(token);
277        }
278
279        Ok(tokens)
280    }
281}
282
283/// Parser for raster expressions
284struct Parser {
285    tokens: Vec<Token>,
286    pos: usize,
287}
288
289impl Parser {
290    fn new(tokens: Vec<Token>) -> Self {
291        Self { tokens, pos: 0 }
292    }
293
294    fn current(&self) -> Option<&Token> {
295        self.tokens.get(self.pos)
296    }
297
298    fn advance(&mut self) {
299        self.pos += 1;
300    }
301
302    fn parse(&mut self) -> Result<Expr> {
303        self.parse_conditional()
304    }
305
306    fn parse_conditional(&mut self) -> Result<Expr> {
307        if matches!(self.current(), Some(Token::If)) {
308            self.advance();
309            let condition = Box::new(self.parse_or()?);
310
311            if !matches!(self.current(), Some(Token::Then)) {
312                return Err(AlgorithmError::InvalidParameter {
313                    parameter: "expression",
314                    message: "Expected 'then' after if condition".to_string(),
315                });
316            }
317            self.advance();
318
319            let then_expr = Box::new(self.parse_or()?);
320
321            if !matches!(self.current(), Some(Token::Else)) {
322                return Err(AlgorithmError::InvalidParameter {
323                    parameter: "expression",
324                    message: "Expected 'else' in conditional".to_string(),
325                });
326            }
327            self.advance();
328
329            let else_expr = Box::new(self.parse_or()?);
330
331            Ok(Expr::Conditional {
332                condition,
333                then_expr,
334                else_expr,
335            })
336        } else {
337            self.parse_or()
338        }
339    }
340
341    fn parse_or(&mut self) -> Result<Expr> {
342        let mut left = self.parse_and()?;
343
344        while matches!(self.current(), Some(Token::Or)) {
345            self.advance();
346            let right = self.parse_and()?;
347            left = Expr::BinaryOp {
348                left: Box::new(left),
349                op: BinaryOp::Or,
350                right: Box::new(right),
351            };
352        }
353
354        Ok(left)
355    }
356
357    fn parse_and(&mut self) -> Result<Expr> {
358        let mut left = self.parse_comparison()?;
359
360        while matches!(self.current(), Some(Token::And)) {
361            self.advance();
362            let right = self.parse_comparison()?;
363            left = Expr::BinaryOp {
364                left: Box::new(left),
365                op: BinaryOp::And,
366                right: Box::new(right),
367            };
368        }
369
370        Ok(left)
371    }
372
373    fn parse_comparison(&mut self) -> Result<Expr> {
374        let mut left = self.parse_additive()?;
375
376        while let Some(token) = self.current() {
377            let op = match token {
378                Token::Greater => BinaryOp::Greater,
379                Token::Less => BinaryOp::Less,
380                Token::GreaterEqual => BinaryOp::GreaterEqual,
381                Token::LessEqual => BinaryOp::LessEqual,
382                Token::Equal => BinaryOp::Equal,
383                Token::NotEqual => BinaryOp::NotEqual,
384                _ => break,
385            };
386
387            self.advance();
388            let right = self.parse_additive()?;
389            left = Expr::BinaryOp {
390                left: Box::new(left),
391                op,
392                right: Box::new(right),
393            };
394        }
395
396        Ok(left)
397    }
398
399    fn parse_additive(&mut self) -> Result<Expr> {
400        let mut left = self.parse_multiplicative()?;
401
402        while let Some(token) = self.current() {
403            let op = match token {
404                Token::Plus => BinaryOp::Add,
405                Token::Minus => BinaryOp::Subtract,
406                _ => break,
407            };
408
409            self.advance();
410            let right = self.parse_multiplicative()?;
411            left = Expr::BinaryOp {
412                left: Box::new(left),
413                op,
414                right: Box::new(right),
415            };
416        }
417
418        Ok(left)
419    }
420
421    fn parse_multiplicative(&mut self) -> Result<Expr> {
422        let mut left = self.parse_power()?;
423
424        while let Some(token) = self.current() {
425            let op = match token {
426                Token::Multiply => BinaryOp::Multiply,
427                Token::Divide => BinaryOp::Divide,
428                _ => break,
429            };
430
431            self.advance();
432            let right = self.parse_power()?;
433            left = Expr::BinaryOp {
434                left: Box::new(left),
435                op,
436                right: Box::new(right),
437            };
438        }
439
440        Ok(left)
441    }
442
443    fn parse_power(&mut self) -> Result<Expr> {
444        let mut left = self.parse_unary()?;
445
446        while matches!(self.current(), Some(Token::Power)) {
447            self.advance();
448            let right = self.parse_unary()?;
449            left = Expr::BinaryOp {
450                left: Box::new(left),
451                op: BinaryOp::Power,
452                right: Box::new(right),
453            };
454        }
455
456        Ok(left)
457    }
458
459    fn parse_unary(&mut self) -> Result<Expr> {
460        if matches!(self.current(), Some(Token::Minus)) {
461            self.advance();
462            let expr = self.parse_unary()?;
463            Ok(Expr::UnaryOp {
464                op: UnaryOp::Negate,
465                expr: Box::new(expr),
466            })
467        } else {
468            self.parse_primary()
469        }
470    }
471
472    fn parse_primary(&mut self) -> Result<Expr> {
473        match self.current() {
474            Some(Token::Number(n)) => {
475                let val = *n;
476                self.advance();
477                Ok(Expr::Number(val))
478            }
479            Some(Token::Band(b)) => {
480                let band = *b;
481                self.advance();
482                Ok(Expr::Band(band))
483            }
484            Some(Token::Ident(name)) => {
485                let func_name = name.clone();
486                self.advance();
487
488                if matches!(self.current(), Some(Token::LeftParen)) {
489                    self.advance();
490                    let mut args = Vec::new();
491
492                    if !matches!(self.current(), Some(Token::RightParen)) {
493                        loop {
494                            args.push(self.parse_conditional()?);
495
496                            if matches!(self.current(), Some(Token::Comma)) {
497                                self.advance();
498                            } else {
499                                break;
500                            }
501                        }
502                    }
503
504                    if !matches!(self.current(), Some(Token::RightParen)) {
505                        return Err(AlgorithmError::InvalidParameter {
506                            parameter: "expression",
507                            message: "Expected ')' after function arguments".to_string(),
508                        });
509                    }
510                    self.advance();
511
512                    Ok(Expr::Function {
513                        name: func_name,
514                        args,
515                    })
516                } else {
517                    Err(AlgorithmError::InvalidParameter {
518                        parameter: "expression",
519                        message: format!("Unknown identifier: {func_name}"),
520                    })
521                }
522            }
523            Some(Token::LeftParen) => {
524                self.advance();
525                let expr = self.parse_conditional()?;
526
527                if !matches!(self.current(), Some(Token::RightParen)) {
528                    return Err(AlgorithmError::InvalidParameter {
529                        parameter: "expression",
530                        message: "Expected ')'".to_string(),
531                    });
532                }
533                self.advance();
534
535                Ok(expr)
536            }
537            _ => Err(AlgorithmError::InvalidParameter {
538                parameter: "expression",
539                message: "Unexpected token in expression".to_string(),
540            }),
541        }
542    }
543}
544
545/// Expression optimizer for constant folding and algebraic simplifications
546struct Optimizer;
547
548impl Optimizer {
549    /// Optimize an expression tree
550    fn optimize(expr: Expr) -> Expr {
551        let expr = Self::constant_fold(expr);
552        let expr = Self::algebraic_simplify(expr);
553        Self::eliminate_common_subexpressions(expr)
554    }
555
556    /// Constant folding: evaluate constant expressions at compile time
557    fn constant_fold(expr: Expr) -> Expr {
558        match expr {
559            Expr::BinaryOp { left, op, right } => {
560                let left = Self::constant_fold(*left);
561                let right = Self::constant_fold(*right);
562
563                // If both sides are constants, evaluate
564                if let (Expr::Number(l), Expr::Number(r)) = (&left, &right) {
565                    let result = match op {
566                        BinaryOp::Add => l + r,
567                        BinaryOp::Subtract => l - r,
568                        BinaryOp::Multiply => l * r,
569                        BinaryOp::Divide => {
570                            if r.abs() < f64::EPSILON {
571                                f64::NAN
572                            } else {
573                                l / r
574                            }
575                        }
576                        BinaryOp::Power => l.powf(*r),
577                        BinaryOp::Greater => {
578                            if l > r {
579                                1.0
580                            } else {
581                                0.0
582                            }
583                        }
584                        BinaryOp::Less => {
585                            if l < r {
586                                1.0
587                            } else {
588                                0.0
589                            }
590                        }
591                        BinaryOp::GreaterEqual => {
592                            if l >= r {
593                                1.0
594                            } else {
595                                0.0
596                            }
597                        }
598                        BinaryOp::LessEqual => {
599                            if l <= r {
600                                1.0
601                            } else {
602                                0.0
603                            }
604                        }
605                        BinaryOp::Equal => {
606                            if (l - r).abs() < f64::EPSILON {
607                                1.0
608                            } else {
609                                0.0
610                            }
611                        }
612                        BinaryOp::NotEqual => {
613                            if (l - r).abs() >= f64::EPSILON {
614                                1.0
615                            } else {
616                                0.0
617                            }
618                        }
619                        BinaryOp::And => {
620                            if *l != 0.0 && *r != 0.0 {
621                                1.0
622                            } else {
623                                0.0
624                            }
625                        }
626                        BinaryOp::Or => {
627                            if *l != 0.0 || *r != 0.0 {
628                                1.0
629                            } else {
630                                0.0
631                            }
632                        }
633                    };
634                    return Expr::Number(result);
635                }
636
637                Expr::BinaryOp {
638                    left: Box::new(left),
639                    op,
640                    right: Box::new(right),
641                }
642            }
643            Expr::UnaryOp { op, expr } => {
644                let expr = Self::constant_fold(*expr);
645                if let Expr::Number(n) = expr {
646                    let result = match op {
647                        UnaryOp::Negate => -n,
648                    };
649                    return Expr::Number(result);
650                }
651                Expr::UnaryOp {
652                    op,
653                    expr: Box::new(expr),
654                }
655            }
656            Expr::Function { name, args } => {
657                let args: Vec<Expr> = args.into_iter().map(Self::constant_fold).collect();
658
659                // If all args are constants, evaluate the function
660                let all_const = args.iter().all(|arg| matches!(arg, Expr::Number(_)));
661                if all_const {
662                    let arg_vals: Vec<f64> = args
663                        .iter()
664                        .filter_map(|arg| {
665                            if let Expr::Number(n) = arg {
666                                Some(*n)
667                            } else {
668                                None
669                            }
670                        })
671                        .collect();
672
673                    let result = match name.as_str() {
674                        "sqrt" if arg_vals.len() == 1 => Some(arg_vals[0].sqrt()),
675                        "abs" if arg_vals.len() == 1 => Some(arg_vals[0].abs()),
676                        "log" if arg_vals.len() == 1 => Some(arg_vals[0].ln()),
677                        "log10" if arg_vals.len() == 1 => Some(arg_vals[0].log10()),
678                        "exp" if arg_vals.len() == 1 => Some(arg_vals[0].exp()),
679                        "sin" if arg_vals.len() == 1 => Some(arg_vals[0].sin()),
680                        "cos" if arg_vals.len() == 1 => Some(arg_vals[0].cos()),
681                        "tan" if arg_vals.len() == 1 => Some(arg_vals[0].tan()),
682                        "floor" if arg_vals.len() == 1 => Some(arg_vals[0].floor()),
683                        "ceil" if arg_vals.len() == 1 => Some(arg_vals[0].ceil()),
684                        "round" if arg_vals.len() == 1 => Some(arg_vals[0].round()),
685                        "min" if !arg_vals.is_empty() => {
686                            Some(arg_vals.iter().copied().fold(f64::INFINITY, f64::min))
687                        }
688                        "max" if !arg_vals.is_empty() => {
689                            Some(arg_vals.iter().copied().fold(f64::NEG_INFINITY, f64::max))
690                        }
691                        _ => None,
692                    };
693
694                    if let Some(val) = result {
695                        return Expr::Number(val);
696                    }
697                }
698
699                Expr::Function { name, args }
700            }
701            Expr::Conditional {
702                condition,
703                then_expr,
704                else_expr,
705            } => {
706                let condition = Self::constant_fold(*condition);
707                let then_expr = Self::constant_fold(*then_expr);
708                let else_expr = Self::constant_fold(*else_expr);
709
710                // If condition is constant, choose branch
711                if let Expr::Number(cond) = condition {
712                    if cond != 0.0 {
713                        return then_expr;
714                    } else {
715                        return else_expr;
716                    }
717                }
718
719                Expr::Conditional {
720                    condition: Box::new(condition),
721                    then_expr: Box::new(then_expr),
722                    else_expr: Box::new(else_expr),
723                }
724            }
725            other => other,
726        }
727    }
728
729    /// Algebraic simplifications: x + 0 -> x, x * 1 -> x, x * 0 -> 0, etc.
730    fn algebraic_simplify(expr: Expr) -> Expr {
731        match expr {
732            Expr::BinaryOp { left, op, right } => {
733                let left = Self::algebraic_simplify(*left);
734                let right = Self::algebraic_simplify(*right);
735
736                match (&left, op, &right) {
737                    // x + 0 = x, 0 + x = x
738                    (_, BinaryOp::Add, Expr::Number(n)) if n.abs() < f64::EPSILON => left,
739                    (Expr::Number(n), BinaryOp::Add, _) if n.abs() < f64::EPSILON => right,
740
741                    // x - 0 = x
742                    (_, BinaryOp::Subtract, Expr::Number(n)) if n.abs() < f64::EPSILON => left,
743
744                    // x * 0 = 0, 0 * x = 0
745                    (_, BinaryOp::Multiply, Expr::Number(n))
746                    | (Expr::Number(n), BinaryOp::Multiply, _)
747                        if n.abs() < f64::EPSILON =>
748                    {
749                        Expr::Number(0.0)
750                    }
751
752                    // x * 1 = x, 1 * x = x
753                    (_, BinaryOp::Multiply, Expr::Number(n)) if (n - 1.0).abs() < f64::EPSILON => {
754                        left
755                    }
756                    (Expr::Number(n), BinaryOp::Multiply, _) if (n - 1.0).abs() < f64::EPSILON => {
757                        right
758                    }
759
760                    // x / 1 = x
761                    (_, BinaryOp::Divide, Expr::Number(n)) if (n - 1.0).abs() < f64::EPSILON => {
762                        left
763                    }
764
765                    // x ^ 0 = 1
766                    (_, BinaryOp::Power, Expr::Number(n)) if n.abs() < f64::EPSILON => {
767                        Expr::Number(1.0)
768                    }
769
770                    // x ^ 1 = x
771                    (_, BinaryOp::Power, Expr::Number(n)) if (n - 1.0).abs() < f64::EPSILON => left,
772
773                    _ => Expr::BinaryOp {
774                        left: Box::new(left),
775                        op,
776                        right: Box::new(right),
777                    },
778                }
779            }
780            Expr::UnaryOp { op, expr } => {
781                let expr = Self::algebraic_simplify(*expr);
782                Expr::UnaryOp {
783                    op,
784                    expr: Box::new(expr),
785                }
786            }
787            Expr::Function { name, args } => {
788                let args = args.into_iter().map(Self::algebraic_simplify).collect();
789                Expr::Function { name, args }
790            }
791            Expr::Conditional {
792                condition,
793                then_expr,
794                else_expr,
795            } => Expr::Conditional {
796                condition: Box::new(Self::algebraic_simplify(*condition)),
797                then_expr: Box::new(Self::algebraic_simplify(*then_expr)),
798                else_expr: Box::new(Self::algebraic_simplify(*else_expr)),
799            },
800            other => other,
801        }
802    }
803
804    /// Common Subexpression Elimination (simplified version)
805    /// In a full implementation, this would detect and cache repeated subexpressions
806    fn eliminate_common_subexpressions(expr: Expr) -> Expr {
807        // For now, this is a placeholder for a more sophisticated CSE pass
808        // A full implementation would:
809        // 1. Build a hash map of expression -> cache variable
810        // 2. Replace repeated expressions with cache lookups
811        // 3. Require changing the evaluation model to support cached values
812        expr
813    }
814}
815
816/// Evaluator for raster expressions
817struct Evaluator<'a> {
818    bands: &'a [RasterBuffer],
819}
820
821impl<'a> Evaluator<'a> {
822    fn new(bands: &'a [RasterBuffer]) -> Self {
823        Self { bands }
824    }
825
826    fn eval_pixel(&self, expr: &Expr, x: u64, y: u64) -> Result<f64> {
827        match expr {
828            Expr::Number(n) => Ok(*n),
829            Expr::Band(b) => {
830                if *b == 0 || *b > self.bands.len() {
831                    return Err(AlgorithmError::InvalidParameter {
832                        parameter: "band",
833                        message: format!("Band {} out of range (1-{})", b, self.bands.len()),
834                    });
835                }
836                self.bands[b - 1]
837                    .get_pixel(x, y)
838                    .map_err(AlgorithmError::Core)
839            }
840            Expr::BinaryOp { left, op, right } => {
841                let lval = self.eval_pixel(left, x, y)?;
842                let rval = self.eval_pixel(right, x, y)?;
843
844                let result = match op {
845                    BinaryOp::Add => lval + rval,
846                    BinaryOp::Subtract => lval - rval,
847                    BinaryOp::Multiply => lval * rval,
848                    BinaryOp::Divide => {
849                        if rval.abs() < f64::EPSILON {
850                            f64::NAN
851                        } else {
852                            lval / rval
853                        }
854                    }
855                    BinaryOp::Power => lval.powf(rval),
856                    BinaryOp::Greater => {
857                        if lval > rval {
858                            1.0
859                        } else {
860                            0.0
861                        }
862                    }
863                    BinaryOp::Less => {
864                        if lval < rval {
865                            1.0
866                        } else {
867                            0.0
868                        }
869                    }
870                    BinaryOp::GreaterEqual => {
871                        if lval >= rval {
872                            1.0
873                        } else {
874                            0.0
875                        }
876                    }
877                    BinaryOp::LessEqual => {
878                        if lval <= rval {
879                            1.0
880                        } else {
881                            0.0
882                        }
883                    }
884                    BinaryOp::Equal => {
885                        if (lval - rval).abs() < f64::EPSILON {
886                            1.0
887                        } else {
888                            0.0
889                        }
890                    }
891                    BinaryOp::NotEqual => {
892                        if (lval - rval).abs() >= f64::EPSILON {
893                            1.0
894                        } else {
895                            0.0
896                        }
897                    }
898                    BinaryOp::And => {
899                        if lval != 0.0 && rval != 0.0 {
900                            1.0
901                        } else {
902                            0.0
903                        }
904                    }
905                    BinaryOp::Or => {
906                        if lval != 0.0 || rval != 0.0 {
907                            1.0
908                        } else {
909                            0.0
910                        }
911                    }
912                };
913
914                Ok(result)
915            }
916            Expr::UnaryOp { op, expr } => {
917                let val = self.eval_pixel(expr, x, y)?;
918                match op {
919                    UnaryOp::Negate => Ok(-val),
920                }
921            }
922            Expr::Function { name, args } => self.eval_function(name, args, x, y),
923            Expr::Conditional {
924                condition,
925                then_expr,
926                else_expr,
927            } => {
928                let cond_val = self.eval_pixel(condition, x, y)?;
929                if cond_val != 0.0 {
930                    self.eval_pixel(then_expr, x, y)
931                } else {
932                    self.eval_pixel(else_expr, x, y)
933                }
934            }
935        }
936    }
937
938    fn eval_function(&self, name: &str, args: &[Expr], x: u64, y: u64) -> Result<f64> {
939        let arg_vals: Result<Vec<f64>> =
940            args.iter().map(|arg| self.eval_pixel(arg, x, y)).collect();
941        let arg_vals = arg_vals?;
942
943        let result = match name {
944            "sqrt" => {
945                if arg_vals.len() != 1 {
946                    return Err(AlgorithmError::InvalidParameter {
947                        parameter: "sqrt",
948                        message: "Expected 1 argument".to_string(),
949                    });
950                }
951                arg_vals[0].sqrt()
952            }
953            "abs" => {
954                if arg_vals.len() != 1 {
955                    return Err(AlgorithmError::InvalidParameter {
956                        parameter: "abs",
957                        message: "Expected 1 argument".to_string(),
958                    });
959                }
960                arg_vals[0].abs()
961            }
962            "log" => {
963                if arg_vals.len() != 1 {
964                    return Err(AlgorithmError::InvalidParameter {
965                        parameter: "log",
966                        message: "Expected 1 argument".to_string(),
967                    });
968                }
969                arg_vals[0].ln()
970            }
971            "log10" => {
972                if arg_vals.len() != 1 {
973                    return Err(AlgorithmError::InvalidParameter {
974                        parameter: "log10",
975                        message: "Expected 1 argument".to_string(),
976                    });
977                }
978                arg_vals[0].log10()
979            }
980            "exp" => {
981                if arg_vals.len() != 1 {
982                    return Err(AlgorithmError::InvalidParameter {
983                        parameter: "exp",
984                        message: "Expected 1 argument".to_string(),
985                    });
986                }
987                arg_vals[0].exp()
988            }
989            "sin" => {
990                if arg_vals.len() != 1 {
991                    return Err(AlgorithmError::InvalidParameter {
992                        parameter: "sin",
993                        message: "Expected 1 argument".to_string(),
994                    });
995                }
996                arg_vals[0].sin()
997            }
998            "cos" => {
999                if arg_vals.len() != 1 {
1000                    return Err(AlgorithmError::InvalidParameter {
1001                        parameter: "cos",
1002                        message: "Expected 1 argument".to_string(),
1003                    });
1004                }
1005                arg_vals[0].cos()
1006            }
1007            "tan" => {
1008                if arg_vals.len() != 1 {
1009                    return Err(AlgorithmError::InvalidParameter {
1010                        parameter: "tan",
1011                        message: "Expected 1 argument".to_string(),
1012                    });
1013                }
1014                arg_vals[0].tan()
1015            }
1016            "floor" => {
1017                if arg_vals.len() != 1 {
1018                    return Err(AlgorithmError::InvalidParameter {
1019                        parameter: "floor",
1020                        message: "Expected 1 argument".to_string(),
1021                    });
1022                }
1023                arg_vals[0].floor()
1024            }
1025            "ceil" => {
1026                if arg_vals.len() != 1 {
1027                    return Err(AlgorithmError::InvalidParameter {
1028                        parameter: "ceil",
1029                        message: "Expected 1 argument".to_string(),
1030                    });
1031                }
1032                arg_vals[0].ceil()
1033            }
1034            "round" => {
1035                if arg_vals.len() != 1 {
1036                    return Err(AlgorithmError::InvalidParameter {
1037                        parameter: "round",
1038                        message: "Expected 1 argument".to_string(),
1039                    });
1040                }
1041                arg_vals[0].round()
1042            }
1043            "min" => {
1044                if arg_vals.is_empty() {
1045                    return Err(AlgorithmError::InvalidParameter {
1046                        parameter: "min",
1047                        message: "Expected at least 1 argument".to_string(),
1048                    });
1049                }
1050                arg_vals.iter().copied().fold(f64::INFINITY, f64::min)
1051            }
1052            "max" => {
1053                if arg_vals.is_empty() {
1054                    return Err(AlgorithmError::InvalidParameter {
1055                        parameter: "max",
1056                        message: "Expected at least 1 argument".to_string(),
1057                    });
1058                }
1059                arg_vals.iter().copied().fold(f64::NEG_INFINITY, f64::max)
1060            }
1061            _ => {
1062                return Err(AlgorithmError::InvalidParameter {
1063                    parameter: "function",
1064                    message: format!("Unknown function: {name}"),
1065                });
1066            }
1067        };
1068
1069        Ok(result)
1070    }
1071}
1072
1073/// Raster calculator for map algebra with expression parsing
1074pub struct RasterCalculator;
1075
1076impl RasterCalculator {
1077    /// Evaluates a raster expression on one or more bands
1078    ///
1079    /// # Arguments
1080    ///
1081    /// * `expression` - The expression to evaluate (e.g., "(B1 - B2) / (B1 + B2)")
1082    /// * `bands` - Input bands (B1, B2, etc.)
1083    ///
1084    /// # Examples
1085    ///
1086    /// NDVI: `"(B1 - B2) / (B1 + B2)"`
1087    /// Conditional: `"if B1 > 100 then B1 * 2 else B1"`
1088    /// Math: `"sqrt(B1 ^ 2 + B2 ^ 2)"`
1089    ///
1090    /// # Errors
1091    ///
1092    /// Returns an error if the expression is invalid or evaluation fails
1093    pub fn evaluate(expression: &str, bands: &[RasterBuffer]) -> Result<RasterBuffer> {
1094        if bands.is_empty() {
1095            return Err(AlgorithmError::EmptyInput {
1096                operation: "evaluate",
1097            });
1098        }
1099
1100        // Check all bands have same dimensions
1101        let width = bands[0].width();
1102        let height = bands[0].height();
1103        for (_i, band) in bands.iter().enumerate().skip(1) {
1104            if band.width() != width || band.height() != height {
1105                return Err(AlgorithmError::InvalidDimensions {
1106                    message: "All bands must have same dimensions",
1107                    actual: band.width() as usize,
1108                    expected: width as usize,
1109                });
1110            }
1111        }
1112
1113        // Tokenize
1114        let mut lexer = Lexer::new(expression);
1115        let tokens = lexer.tokenize()?;
1116
1117        // Parse
1118        let mut parser = Parser::new(tokens);
1119        let expr = parser.parse()?;
1120
1121        // Optimize expression
1122        let expr = Optimizer::optimize(expr);
1123
1124        // Evaluate
1125        let evaluator = Evaluator::new(bands);
1126        let mut result = RasterBuffer::zeros(width, height, bands[0].data_type());
1127
1128        for y in 0..height {
1129            for x in 0..width {
1130                let value = evaluator.eval_pixel(&expr, x, y)?;
1131                result
1132                    .set_pixel(x, y, value)
1133                    .map_err(AlgorithmError::Core)?;
1134            }
1135        }
1136
1137        Ok(result)
1138    }
1139
1140    /// Evaluates a raster expression in parallel using rayon
1141    ///
1142    /// This method processes rows in parallel for improved performance on multi-core systems.
1143    /// Falls back to sequential evaluation if the parallel feature is not enabled.
1144    ///
1145    /// # Arguments
1146    ///
1147    /// * `expression` - The expression to evaluate (e.g., "(B1 - B2) / (B1 + B2)")
1148    /// * `bands` - Input bands (B1, B2, etc.)
1149    ///
1150    /// # Examples
1151    ///
1152    /// NDVI: `"(B1 - B2) / (B1 + B2)"`
1153    /// Conditional: `"if B1 > 100 then B1 * 2 else B1"`
1154    /// Math: `"sqrt(B1 ^ 2 + B2 ^ 2)"`
1155    ///
1156    /// # Errors
1157    ///
1158    /// Returns an error if the expression is invalid or evaluation fails
1159    #[cfg(feature = "parallel")]
1160    pub fn evaluate_parallel(expression: &str, bands: &[RasterBuffer]) -> Result<RasterBuffer> {
1161        if bands.is_empty() {
1162            return Err(AlgorithmError::EmptyInput {
1163                operation: "evaluate_parallel",
1164            });
1165        }
1166
1167        // Check all bands have same dimensions
1168        let width = bands[0].width();
1169        let height = bands[0].height();
1170        for band in bands.iter().skip(1) {
1171            if band.width() != width || band.height() != height {
1172                return Err(AlgorithmError::InvalidDimensions {
1173                    message: "All bands must have same dimensions",
1174                    actual: band.width() as usize,
1175                    expected: width as usize,
1176                });
1177            }
1178        }
1179
1180        // Tokenize
1181        let mut lexer = Lexer::new(expression);
1182        let tokens = lexer.tokenize()?;
1183
1184        // Parse
1185        let mut parser = Parser::new(tokens);
1186        let expr = parser.parse()?;
1187
1188        // Optimize expression
1189        let expr = Optimizer::optimize(expr);
1190
1191        // Create evaluator
1192        let evaluator = Evaluator::new(bands);
1193
1194        // Create result buffer
1195        let mut result = RasterBuffer::zeros(width, height, bands[0].data_type());
1196
1197        // Process rows in parallel
1198        let row_data: Result<Vec<Vec<f64>>> = (0..height)
1199            .into_par_iter()
1200            .map(|y| {
1201                let mut row = Vec::with_capacity(width as usize);
1202                for x in 0..width {
1203                    let value = evaluator.eval_pixel(&expr, x, y)?;
1204                    row.push(value);
1205                }
1206                Ok(row)
1207            })
1208            .collect();
1209
1210        let row_data = row_data?;
1211
1212        // Write results back to buffer
1213        for (y, row) in row_data.iter().enumerate() {
1214            for (x, &value) in row.iter().enumerate() {
1215                result
1216                    .set_pixel(x as u64, y as u64, value)
1217                    .map_err(AlgorithmError::Core)?;
1218            }
1219        }
1220
1221        Ok(result)
1222    }
1223
1224    /// Applies a binary operation to two rasters (legacy API)
1225    pub fn apply_binary(
1226        a: &RasterBuffer,
1227        b: &RasterBuffer,
1228        op: RasterExpression,
1229    ) -> Result<RasterBuffer> {
1230        if a.width() != b.width() || a.height() != b.height() {
1231            return Err(AlgorithmError::InvalidDimensions {
1232                message: "Rasters must have same dimensions",
1233                actual: a.width() as usize,
1234                expected: b.width() as usize,
1235            });
1236        }
1237
1238        let mut result = RasterBuffer::zeros(a.width(), a.height(), a.data_type());
1239
1240        for y in 0..a.height() {
1241            for x in 0..a.width() {
1242                let val_a = a.get_pixel(x, y).map_err(AlgorithmError::Core)?;
1243                let val_b = b.get_pixel(x, y).map_err(AlgorithmError::Core)?;
1244
1245                let val = match op {
1246                    RasterExpression::Add => val_a + val_b,
1247                    RasterExpression::Subtract => val_a - val_b,
1248                    RasterExpression::Multiply => val_a * val_b,
1249                    RasterExpression::Divide => {
1250                        if val_b.abs() < f64::EPSILON {
1251                            f64::NAN
1252                        } else {
1253                            val_a / val_b
1254                        }
1255                    }
1256                    RasterExpression::Max => val_a.max(val_b),
1257                    RasterExpression::Min => val_a.min(val_b),
1258                };
1259
1260                result.set_pixel(x, y, val).map_err(AlgorithmError::Core)?;
1261            }
1262        }
1263
1264        Ok(result)
1265    }
1266
1267    /// Applies a unary function to a raster (legacy API)
1268    pub fn apply_unary<F>(src: &RasterBuffer, func: F) -> Result<RasterBuffer>
1269    where
1270        F: Fn(f64) -> f64,
1271    {
1272        let mut result = RasterBuffer::zeros(src.width(), src.height(), src.data_type());
1273
1274        for y in 0..src.height() {
1275            for x in 0..src.width() {
1276                let val = src.get_pixel(x, y).map_err(AlgorithmError::Core)?;
1277                let new_val = func(val);
1278                result
1279                    .set_pixel(x, y, new_val)
1280                    .map_err(AlgorithmError::Core)?;
1281            }
1282        }
1283
1284        Ok(result)
1285    }
1286}
1287
1288/// Raster expression operations (legacy API)
1289#[derive(Debug, Clone, Copy, PartialEq)]
1290pub enum RasterExpression {
1291    /// Add two rasters
1292    Add,
1293    /// Subtract rasters
1294    Subtract,
1295    /// Multiply rasters
1296    Multiply,
1297    /// Divide rasters
1298    Divide,
1299    /// Maximum of two rasters
1300    Max,
1301    /// Minimum of two rasters
1302    Min,
1303}
1304
1305#[cfg(test)]
1306#[allow(clippy::panic, clippy::cloned_ref_to_slice_refs)]
1307mod tests {
1308    use super::*;
1309    use oxigdal_core::types::RasterDataType;
1310
1311    // ========== Basic Functionality Tests ==========
1312
1313    #[test]
1314    fn test_simple_arithmetic() {
1315        let mut b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1316        let mut b2 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1317
1318        for y in 0..5 {
1319            for x in 0..5 {
1320                b1.set_pixel(x, y, 10.0).ok();
1321                b2.set_pixel(x, y, 5.0).ok();
1322            }
1323        }
1324
1325        let result = RasterCalculator::evaluate("B1 + B2", &[b1, b2]);
1326        assert!(result.is_ok());
1327        let r = result.expect("Result should be ok");
1328        let val = r.get_pixel(0, 0).expect("Should get pixel");
1329        assert!((val - 15.0).abs() < f64::EPSILON);
1330    }
1331
1332    #[test]
1333    fn test_ndvi() {
1334        let mut nir = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1335        let mut red = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1336
1337        for y in 0..5 {
1338            for x in 0..5 {
1339                nir.set_pixel(x, y, 100.0).ok();
1340                red.set_pixel(x, y, 50.0).ok();
1341            }
1342        }
1343
1344        let result = RasterCalculator::evaluate("(B1 - B2) / (B1 + B2)", &[nir, red]);
1345        assert!(result.is_ok());
1346        let r = result.expect("Result should be ok");
1347        let val = r.get_pixel(0, 0).expect("Should get pixel");
1348        let expected = (100.0 - 50.0) / (100.0 + 50.0);
1349        assert!((val - expected).abs() < 0.001);
1350    }
1351
1352    #[test]
1353    fn test_math_functions() {
1354        let mut b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1355
1356        for y in 0..5 {
1357            for x in 0..5 {
1358                b1.set_pixel(x, y, 16.0).ok();
1359            }
1360        }
1361
1362        let result = RasterCalculator::evaluate("sqrt(B1)", &[b1]);
1363        assert!(result.is_ok());
1364        let r = result.expect("Result should be ok");
1365        let val = r.get_pixel(0, 0).expect("Should get pixel");
1366        assert!((val - 4.0).abs() < f64::EPSILON);
1367    }
1368
1369    #[test]
1370    fn test_conditional() {
1371        let mut b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1372
1373        for y in 0..5 {
1374            for x in 0..5 {
1375                b1.set_pixel(x, y, (x * 10) as f64).ok();
1376            }
1377        }
1378
1379        let result = RasterCalculator::evaluate("if B1 > 20 then 1 else 0", &[b1]);
1380        assert!(result.is_ok());
1381        let r = result.expect("Result should be ok");
1382
1383        let val0 = r.get_pixel(0, 0).expect("Should get pixel");
1384        assert!((val0 - 0.0).abs() < f64::EPSILON);
1385
1386        let val3 = r.get_pixel(3, 0).expect("Should get pixel");
1387        assert!((val3 - 1.0).abs() < f64::EPSILON);
1388    }
1389
1390    #[test]
1391    fn test_legacy_add() {
1392        let mut a = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1393        let mut b = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1394
1395        a.set_pixel(0, 0, 10.0).ok();
1396        b.set_pixel(0, 0, 5.0).ok();
1397
1398        let result = RasterCalculator::apply_binary(&a, &b, RasterExpression::Add);
1399        assert!(result.is_ok());
1400    }
1401
1402    // ========== Edge Cases ==========
1403
1404    #[test]
1405    fn test_empty_bands() {
1406        let result = RasterCalculator::evaluate("B1 + B2", &[]);
1407        assert!(result.is_err());
1408        if let Err(AlgorithmError::EmptyInput { .. }) = result {
1409            // Expected
1410        } else {
1411            panic!("Expected EmptyInput error");
1412        }
1413    }
1414
1415    #[test]
1416    fn test_single_pixel() {
1417        let mut b1 = RasterBuffer::zeros(1, 1, RasterDataType::Float32);
1418        b1.set_pixel(0, 0, 42.0).ok();
1419
1420        let result = RasterCalculator::evaluate("B1 * 2", &[b1]);
1421        assert!(result.is_ok());
1422        let r = result.expect("Should succeed");
1423        let val = r.get_pixel(0, 0).expect("Should get pixel");
1424        assert!((val - 84.0).abs() < f64::EPSILON);
1425    }
1426
1427    #[test]
1428    fn test_division_by_zero() {
1429        let mut b1 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1430        let mut b2 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1431
1432        for y in 0..3 {
1433            for x in 0..3 {
1434                b1.set_pixel(x, y, 10.0).ok();
1435                b2.set_pixel(x, y, 0.0).ok();
1436            }
1437        }
1438
1439        let result = RasterCalculator::evaluate("B1 / B2", &[b1, b2]);
1440        assert!(result.is_ok());
1441        let r = result.expect("Should succeed");
1442        let val = r.get_pixel(0, 0).expect("Should get pixel");
1443        assert!(val.is_nan()); // Division by zero should give NaN
1444    }
1445
1446    #[test]
1447    fn test_mismatched_dimensions() {
1448        let b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1449        let b2 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1450
1451        let result = RasterCalculator::evaluate("B1 + B2", &[b1, b2]);
1452        assert!(result.is_err());
1453        if let Err(AlgorithmError::InvalidDimensions { .. }) = result {
1454            // Expected
1455        } else {
1456            panic!("Expected InvalidDimensions error");
1457        }
1458    }
1459
1460    // ========== Error Conditions ==========
1461
1462    #[test]
1463    fn test_invalid_band_reference() {
1464        let b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1465
1466        let result = RasterCalculator::evaluate("B5 + B1", &[b1]);
1467        assert!(result.is_err());
1468    }
1469
1470    #[test]
1471    fn test_undefined_function() {
1472        let b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1473
1474        let result = RasterCalculator::evaluate("undefined_func(B1)", &[b1]);
1475        assert!(result.is_err());
1476    }
1477
1478    #[test]
1479    fn test_invalid_expression() {
1480        let b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1481
1482        let result = RasterCalculator::evaluate("B1 +", &[b1]);
1483        assert!(result.is_err());
1484    }
1485
1486    #[test]
1487    fn test_mismatched_parentheses() {
1488        let b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1489
1490        let result = RasterCalculator::evaluate("(B1 + 10", &[b1]);
1491        assert!(result.is_err());
1492    }
1493
1494    #[test]
1495    fn test_wrong_function_arity() {
1496        let b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1497
1498        let result = RasterCalculator::evaluate("sqrt(B1, B1)", &[b1]);
1499        assert!(result.is_err());
1500    }
1501
1502    // ========== Complex Operations ==========
1503
1504    #[test]
1505    fn test_nested_functions() {
1506        let mut b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1507
1508        for y in 0..5 {
1509            for x in 0..5 {
1510                b1.set_pixel(x, y, 9.0).ok();
1511            }
1512        }
1513
1514        let result = RasterCalculator::evaluate("sqrt(sqrt(B1))", &[b1]);
1515        assert!(result.is_ok());
1516        let r = result.expect("Should succeed");
1517        let val = r.get_pixel(0, 0).expect("Should get pixel");
1518        let expected = 9.0_f64.sqrt().sqrt();
1519        assert!((val - expected).abs() < 0.001);
1520    }
1521
1522    #[test]
1523    fn test_complex_expression() {
1524        let mut b1 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1525        let mut b2 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1526        let mut b3 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1527
1528        for y in 0..3 {
1529            for x in 0..3 {
1530                b1.set_pixel(x, y, 10.0).ok();
1531                b2.set_pixel(x, y, 5.0).ok();
1532                b3.set_pixel(x, y, 2.0).ok();
1533            }
1534        }
1535
1536        let result = RasterCalculator::evaluate("(B1 + B2) * B3 - sqrt(B1)", &[b1, b2, b3]);
1537        assert!(result.is_ok());
1538        let r = result.expect("Should succeed");
1539        let val = r.get_pixel(0, 0).expect("Should get pixel");
1540        let expected = (10.0 + 5.0) * 2.0 - 10.0_f64.sqrt();
1541        assert!((val - expected).abs() < 0.001);
1542    }
1543
1544    #[test]
1545    fn test_all_math_functions() {
1546        let mut b1 = RasterBuffer::zeros(2, 2, RasterDataType::Float32);
1547
1548        for y in 0..2 {
1549            for x in 0..2 {
1550                b1.set_pixel(x, y, 2.5).ok();
1551            }
1552        }
1553
1554        // Test each function
1555        let functions = vec![
1556            ("abs(B1)", 2.5),
1557            ("floor(B1)", 2.0),
1558            ("ceil(B1)", 3.0),
1559            ("round(B1)", 3.0), // rounds to nearest (2.5 -> 3.0)
1560            ("exp(B1)", 2.5_f64.exp()),
1561            ("log(B1)", 2.5_f64.ln()),
1562            ("log10(B1)", 2.5_f64.log10()),
1563            ("sqrt(B1)", 2.5_f64.sqrt()),
1564            ("sin(B1)", 2.5_f64.sin()),
1565            ("cos(B1)", 2.5_f64.cos()),
1566            ("tan(B1)", 2.5_f64.tan()),
1567        ];
1568
1569        for (expr, expected) in functions {
1570            let result = RasterCalculator::evaluate(expr, &[b1.clone()]);
1571            assert!(result.is_ok(), "Failed for expression: {}", expr);
1572            let r = result.expect("Should succeed");
1573            let val = r.get_pixel(0, 0).expect("Should get pixel");
1574            assert!(
1575                (val - expected).abs() < 0.001,
1576                "Expression {} expected {} but got {}",
1577                expr,
1578                expected,
1579                val
1580            );
1581        }
1582    }
1583
1584    #[test]
1585    fn test_min_max_functions() {
1586        let mut b1 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1587        let mut b2 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1588
1589        for y in 0..3 {
1590            for x in 0..3 {
1591                b1.set_pixel(x, y, 10.0).ok();
1592                b2.set_pixel(x, y, 20.0).ok();
1593            }
1594        }
1595
1596        let result = RasterCalculator::evaluate("min(B1, B2)", &[b1.clone(), b2.clone()]);
1597        assert!(result.is_ok());
1598        let r = result.expect("Should succeed");
1599        let val = r.get_pixel(0, 0).expect("Should get pixel");
1600        assert!((val - 10.0).abs() < f64::EPSILON);
1601
1602        let result = RasterCalculator::evaluate("max(B1, B2)", &[b1, b2]);
1603        assert!(result.is_ok());
1604        let r = result.expect("Should succeed");
1605        let val = r.get_pixel(0, 0).expect("Should get pixel");
1606        assert!((val - 20.0).abs() < f64::EPSILON);
1607    }
1608
1609    #[test]
1610    fn test_power_operation() {
1611        let mut b1 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1612
1613        for y in 0..3 {
1614            for x in 0..3 {
1615                b1.set_pixel(x, y, 2.0).ok();
1616            }
1617        }
1618
1619        let result = RasterCalculator::evaluate("B1 ^ 3", &[b1]);
1620        assert!(result.is_ok());
1621        let r = result.expect("Should succeed");
1622        let val = r.get_pixel(0, 0).expect("Should get pixel");
1623        assert!((val - 8.0).abs() < f64::EPSILON);
1624    }
1625
1626    #[test]
1627    fn test_comparison_operators() {
1628        let mut b1 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1629
1630        for y in 0..3 {
1631            for x in 0..3 {
1632                b1.set_pixel(x, y, 10.0).ok();
1633            }
1634        }
1635
1636        let tests = vec![
1637            ("B1 > 5", 1.0),
1638            ("B1 < 5", 0.0),
1639            ("B1 >= 10", 1.0),
1640            ("B1 <= 10", 1.0),
1641            ("B1 == 10", 1.0),
1642            ("B1 != 5", 1.0),
1643        ];
1644
1645        for (expr, expected) in tests {
1646            let result = RasterCalculator::evaluate(expr, &[b1.clone()]);
1647            assert!(result.is_ok(), "Failed for expression: {}", expr);
1648            let r = result.expect("Should succeed");
1649            let val = r.get_pixel(0, 0).expect("Should get pixel");
1650            assert!(
1651                (val - expected).abs() < f64::EPSILON,
1652                "Expression {} expected {} but got {}",
1653                expr,
1654                expected,
1655                val
1656            );
1657        }
1658    }
1659
1660    #[test]
1661    fn test_logical_operators() {
1662        let mut b1 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1663
1664        for y in 0..3 {
1665            for x in 0..3 {
1666                b1.set_pixel(x, y, 10.0).ok();
1667            }
1668        }
1669
1670        let result = RasterCalculator::evaluate("B1 > 5 and B1 < 15", &[b1.clone()]);
1671        assert!(result.is_ok());
1672        let r = result.expect("Should succeed");
1673        let val = r.get_pixel(0, 0).expect("Should get pixel");
1674        assert!((val - 1.0).abs() < f64::EPSILON);
1675
1676        let result = RasterCalculator::evaluate("B1 < 5 or B1 > 5", &[b1]);
1677        assert!(result.is_ok());
1678        let r = result.expect("Should succeed");
1679        let val = r.get_pixel(0, 0).expect("Should get pixel");
1680        assert!((val - 1.0).abs() < f64::EPSILON);
1681    }
1682
1683    #[test]
1684    fn test_nested_conditionals() {
1685        let mut b1 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1686
1687        for y in 0..3 {
1688            for x in 0..3 {
1689                b1.set_pixel(x, y, (x * 10) as f64).ok();
1690            }
1691        }
1692
1693        // Nested conditionals using proper syntax
1694        let result =
1695            RasterCalculator::evaluate("if B1 > 15 then 3 else (if B1 > 5 then 2 else 1)", &[b1]);
1696        assert!(result.is_ok());
1697        let r = result.expect("Should succeed");
1698
1699        let val0 = r.get_pixel(0, 0).expect("Should get pixel");
1700        assert!((val0 - 1.0).abs() < f64::EPSILON);
1701
1702        let val1 = r.get_pixel(1, 0).expect("Should get pixel");
1703        assert!((val1 - 2.0).abs() < f64::EPSILON);
1704
1705        let val2 = r.get_pixel(2, 0).expect("Should get pixel");
1706        assert!((val2 - 3.0).abs() < f64::EPSILON);
1707    }
1708
1709    // ========== Legacy API Tests ==========
1710
1711    #[test]
1712    fn test_legacy_operations() {
1713        let mut a = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1714        let mut b = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1715
1716        for y in 0..3 {
1717            for x in 0..3 {
1718                a.set_pixel(x, y, 10.0).ok();
1719                b.set_pixel(x, y, 5.0).ok();
1720            }
1721        }
1722
1723        let operations = vec![
1724            (RasterExpression::Add, 15.0),
1725            (RasterExpression::Subtract, 5.0),
1726            (RasterExpression::Multiply, 50.0),
1727            (RasterExpression::Divide, 2.0),
1728            (RasterExpression::Max, 10.0),
1729            (RasterExpression::Min, 5.0),
1730        ];
1731
1732        for (op, expected) in operations {
1733            let result = RasterCalculator::apply_binary(&a, &b, op);
1734            assert!(result.is_ok());
1735            let r = result.expect("Should succeed");
1736            let val = r.get_pixel(0, 0).expect("Should get pixel");
1737            assert!(
1738                (val - expected).abs() < f64::EPSILON,
1739                "Operation {:?} expected {} but got {}",
1740                op,
1741                expected,
1742                val
1743            );
1744        }
1745    }
1746
1747    #[test]
1748    fn test_legacy_unary() {
1749        let mut src = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1750
1751        for y in 0..3 {
1752            for x in 0..3 {
1753                src.set_pixel(x, y, 5.0).ok();
1754            }
1755        }
1756
1757        let result = RasterCalculator::apply_unary(&src, |x| x * 2.0);
1758        assert!(result.is_ok());
1759        let r = result.expect("Should succeed");
1760        let val = r.get_pixel(0, 0).expect("Should get pixel");
1761        assert!((val - 10.0).abs() < f64::EPSILON);
1762    }
1763
1764    #[test]
1765    fn test_unary_negate() {
1766        let mut b1 = RasterBuffer::zeros(3, 3, RasterDataType::Float32);
1767
1768        for y in 0..3 {
1769            for x in 0..3 {
1770                b1.set_pixel(x, y, 10.0).ok();
1771            }
1772        }
1773
1774        let result = RasterCalculator::evaluate("-B1", &[b1]);
1775        assert!(result.is_ok());
1776        let r = result.expect("Should succeed");
1777        let val = r.get_pixel(0, 0).expect("Should get pixel");
1778        assert!((val + 10.0).abs() < f64::EPSILON);
1779    }
1780
1781    // ========== Optimizer Tests ==========
1782
1783    #[test]
1784    fn test_optimizer_constant_folding() {
1785        let b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1786
1787        // Constant expression should be pre-computed
1788        let result = RasterCalculator::evaluate("2 + 3", &[b1.clone()]);
1789        assert!(result.is_ok());
1790        let r = result.expect("Should succeed");
1791        let val = r.get_pixel(0, 0).expect("Should get pixel");
1792        assert!((val - 5.0).abs() < f64::EPSILON);
1793
1794        // Constant function call
1795        let result = RasterCalculator::evaluate("sqrt(16)", &[b1.clone()]);
1796        assert!(result.is_ok());
1797        let r = result.expect("Should succeed");
1798        let val = r.get_pixel(0, 0).expect("Should get pixel");
1799        assert!((val - 4.0).abs() < f64::EPSILON);
1800
1801        // Mixed constant and variable
1802        let mut b2 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1803        for y in 0..5 {
1804            for x in 0..5 {
1805                b2.set_pixel(x, y, 10.0).ok();
1806            }
1807        }
1808        let result = RasterCalculator::evaluate("B1 + (2 + 3)", &[b2]);
1809        assert!(result.is_ok());
1810        let r = result.expect("Should succeed");
1811        let val = r.get_pixel(0, 0).expect("Should get pixel");
1812        assert!((val - 15.0).abs() < f64::EPSILON);
1813    }
1814
1815    #[test]
1816    fn test_optimizer_algebraic_simplification() {
1817        let mut b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1818        for y in 0..5 {
1819            for x in 0..5 {
1820                b1.set_pixel(x, y, 10.0).ok();
1821            }
1822        }
1823
1824        // x + 0 = x
1825        let result = RasterCalculator::evaluate("B1 + 0", &[b1.clone()]);
1826        assert!(result.is_ok());
1827        let r = result.expect("Should succeed");
1828        let val = r.get_pixel(0, 0).expect("Should get pixel");
1829        assert!((val - 10.0).abs() < f64::EPSILON);
1830
1831        // x * 1 = x
1832        let result = RasterCalculator::evaluate("B1 * 1", &[b1.clone()]);
1833        assert!(result.is_ok());
1834        let r = result.expect("Should succeed");
1835        let val = r.get_pixel(0, 0).expect("Should get pixel");
1836        assert!((val - 10.0).abs() < f64::EPSILON);
1837
1838        // x * 0 = 0
1839        let result = RasterCalculator::evaluate("B1 * 0", &[b1.clone()]);
1840        assert!(result.is_ok());
1841        let r = result.expect("Should succeed");
1842        let val = r.get_pixel(0, 0).expect("Should get pixel");
1843        assert!(val.abs() < f64::EPSILON);
1844
1845        // x ^ 1 = x
1846        let result = RasterCalculator::evaluate("B1 ^ 1", &[b1]);
1847        assert!(result.is_ok());
1848        let r = result.expect("Should succeed");
1849        let val = r.get_pixel(0, 0).expect("Should get pixel");
1850        assert!((val - 10.0).abs() < f64::EPSILON);
1851    }
1852
1853    #[test]
1854    fn test_optimizer_conditional_constant() {
1855        let mut b1 = RasterBuffer::zeros(5, 5, RasterDataType::Float32);
1856        for y in 0..5 {
1857            for x in 0..5 {
1858                b1.set_pixel(x, y, 10.0).ok();
1859            }
1860        }
1861
1862        // Constant true condition
1863        let result = RasterCalculator::evaluate("if 1 then B1 else B1 * 2", &[b1.clone()]);
1864        assert!(result.is_ok());
1865        let r = result.expect("Should succeed");
1866        let val = r.get_pixel(0, 0).expect("Should get pixel");
1867        assert!((val - 10.0).abs() < f64::EPSILON);
1868
1869        // Constant false condition
1870        let result = RasterCalculator::evaluate("if 0 then B1 else B1 * 2", &[b1]);
1871        assert!(result.is_ok());
1872        let r = result.expect("Should succeed");
1873        let val = r.get_pixel(0, 0).expect("Should get pixel");
1874        assert!((val - 20.0).abs() < f64::EPSILON);
1875    }
1876
1877    // ========== Parallel Evaluation Tests ==========
1878
1879    #[test]
1880    #[cfg(feature = "parallel")]
1881    fn test_parallel_evaluation() {
1882        let mut b1 = RasterBuffer::zeros(100, 100, RasterDataType::Float32);
1883        let mut b2 = RasterBuffer::zeros(100, 100, RasterDataType::Float32);
1884
1885        for y in 0..100 {
1886            for x in 0..100 {
1887                b1.set_pixel(x, y, (x + y) as f64).ok();
1888                b2.set_pixel(x, y, (x * y) as f64).ok();
1889            }
1890        }
1891
1892        let result = RasterCalculator::evaluate_parallel("B1 + B2", &[b1.clone(), b2.clone()]);
1893        assert!(result.is_ok());
1894        let r = result.expect("Should succeed");
1895
1896        // Verify some values
1897        for y in 0..100 {
1898            for x in 0..100 {
1899                let val = r.get_pixel(x, y).expect("Should get pixel");
1900                let expected = (x + y) as f64 + (x * y) as f64;
1901                assert!(
1902                    (val - expected).abs() < f64::EPSILON,
1903                    "Mismatch at ({}, {}): expected {}, got {}",
1904                    x,
1905                    y,
1906                    expected,
1907                    val
1908                );
1909            }
1910        }
1911    }
1912
1913    #[test]
1914    #[cfg(feature = "parallel")]
1915    fn test_parallel_complex_expression() {
1916        let mut nir = RasterBuffer::zeros(50, 50, RasterDataType::Float32);
1917        let mut red = RasterBuffer::zeros(50, 50, RasterDataType::Float32);
1918
1919        for y in 0..50 {
1920            for x in 0..50 {
1921                nir.set_pixel(x, y, 100.0 + x as f64).ok();
1922                red.set_pixel(x, y, 50.0 + y as f64).ok();
1923            }
1924        }
1925
1926        let result = RasterCalculator::evaluate_parallel(
1927            "(B1 - B2) / (B1 + B2)",
1928            &[nir.clone(), red.clone()],
1929        );
1930        assert!(result.is_ok());
1931        let r = result.expect("Should succeed");
1932
1933        // Verify NDVI calculation
1934        for y in 0..50 {
1935            for x in 0..50 {
1936                let nir_val = 100.0 + x as f64;
1937                let red_val = 50.0 + y as f64;
1938                let expected = (nir_val - red_val) / (nir_val + red_val);
1939                let val = r.get_pixel(x, y).expect("Should get pixel");
1940                assert!(
1941                    (val - expected).abs() < 0.001,
1942                    "NDVI mismatch at ({}, {}): expected {}, got {}",
1943                    x,
1944                    y,
1945                    expected,
1946                    val
1947                );
1948            }
1949        }
1950    }
1951
1952    #[test]
1953    #[cfg(feature = "parallel")]
1954    fn test_parallel_with_optimization() {
1955        let mut b1 = RasterBuffer::zeros(50, 50, RasterDataType::Float32);
1956        for y in 0..50 {
1957            for x in 0..50 {
1958                b1.set_pixel(x, y, x as f64).ok();
1959            }
1960        }
1961
1962        // Expression with constants that should be optimized
1963        let result = RasterCalculator::evaluate_parallel("B1 * 1 + 0 + sqrt(16)", &[b1]);
1964        assert!(result.is_ok());
1965        let r = result.expect("Should succeed");
1966
1967        for y in 0..50 {
1968            for x in 0..50 {
1969                let val = r.get_pixel(x, y).expect("Should get pixel");
1970                let expected = x as f64 + 4.0; // B1 + 4 (after optimization)
1971                assert!(
1972                    (val - expected).abs() < f64::EPSILON,
1973                    "Mismatch at ({}, {}): expected {}, got {}",
1974                    x,
1975                    y,
1976                    expected,
1977                    val
1978                );
1979            }
1980        }
1981    }
1982}
oxigdal_algorithms/raster/calculator.rs

oxigdal_algorithms/raster/
calculator.rs
