imgal 0.3.0

A fast and open-source scientific image processing and algorithm library.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

use ndarray::{ArrayBase, AsArray, Dimension, ViewRepr};
use rayon::prelude::*;

use crate::prelude::*;

/// Compute the sum of an n-dimensional image using Kahan compensated summation.
///
/// # Description
///
/// Computes the Kahan sum of an n-dimensional image. The Kahan compensated
/// summation algorithm corrects for floating-point rounding errors and
/// precision loss at each step of the summation. To compensate for
/// floating-point error an error residual is subtracted from each value per
/// iteration.
///
/// # Arguments
///
/// * `data`: The input n-dimensional image.
///
/// # Returns
///
/// * `Ok(T)`: The Kahan sum.
/// * `Err(ImgalError)`: If `data.is_empty() == true`.
#[inline]
pub fn kahan_sum<'a, T, A, D>(data: A) -> Result<T, ImgalError>
where
    A: AsArray<'a, T, D>,
    D: Dimension,
    T: 'a + AsNumeric,
{
    let data: ArrayBase<ViewRepr<&'a T>, D> = data.into();
    if data.is_empty() {
        return Err(ImgalError::InvalidParameterEmptyArray { param_name: "data" });
    }
    Ok(data
        .iter()
        .fold((T::default(), T::default()), |acc, &v| {
            let adj = v - acc.1;
            let new_sum = acc.0 + adj;
            let comp = (new_sum - acc.0) - adj;
            (new_sum, comp)
        })
        .0)
}

/// Compute the sum of an n-dimensional image.
///
/// # Description
///
/// Computes the sum of numerical values in an n-dimensional image.
///
/// # Arguments
///
/// * `data`: The input n-dimensional image.
/// * `threads`: The requested number of threads to use for parallel execution.
///   If `None` or `Some(1)` sequential execution is used. If `Some(0)`, then
///   the maximum available parallelism is used. Thread counts are clamped to
///   the systems maximum.
///
/// # Returns
///
/// * `T`: The sum.
///
/// # Example
///
/// ```
/// use ndarray::Array1;
///
/// use imgal::statistics::sum;
///
/// let arr = [1.82, 3.35, 7.13, 9.25];
/// let total = sum(&arr, None);
/// assert_eq!(total, 21.55);
/// ```
#[inline]
pub fn sum<'a, T, A, D>(data: A, threads: Option<usize>) -> T
where
    A: AsArray<'a, T, D>,
    D: Dimension,
    T: 'a + AsNumeric,
{
    let data: ArrayBase<ViewRepr<&'a T>, D> = data.into();
    par!(threads,
        seq_exp: data.iter().fold(T::default(), |acc, &v| acc + v),
        par_exp: data.into_par_iter()
            .fold(|| T::default(), |acc, &v| acc + v)
            .reduce(|| T::default(), |acc, v| acc + v))
}