imgal 0.3.1

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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

use std::cmp::Ordering;

use ndarray::{
    Array, ArrayBase, ArrayD, ArrayView1, ArrayViewMut1, AsArray, Axis, Dimension, IxDyn, ViewRepr,
};
use rayon::prelude::*;

use crate::copy::copy_into_flat;
use crate::prelude::*;

/// Compute the linear percentile over an n-dimensional image.
///
/// # Description
///
/// Calculates percentiles using linear interpolation between data points. The
/// computation can be performed either on the entire array (flattened) or along
/// a specified axis. The linear percentile is computed using:
///
/// ```text
/// h = (n - 1) × p
/// j = ⌊h⌋
/// γ = h - j
/// percentile = (1 - γ) × v[j] + γ × v[j+1]
/// ```
///
/// Where:
/// - `n` is the array length.
/// - `p` is the percentile in range `0` to `100`.
/// - `v[j]` is the value at index `j`.
/// - `⌊h⌋` is the floor function of `h`.
///
/// When `γ` is close to zero (within `epsilon`), the result is simply `v[j]`,
/// avoiding unnecessary interpolation.
///
/// # Arguments
///
/// * `data`: An n-dimensional image.
/// * `percentile`: The percentile value in thae range `0.0` to `100.0`. Values
///   out side this range will be clamped.
/// * `axis`: The axis to compute percentiles along. If `None`, the input `data`
///   is flattened and a single percentile value is returned.
/// * `epsilon`: The tolerance value used to decide the if the fractional index
///   is an integer. If `None`, then `epsilon = 1e-12`.
/// * `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
///
/// * `Ok(ArrayD<f64>)`: The linear percentile of the input data. If `axis` is
///   `None`, the result shape is `(1,)` and contains a single percentile value
///   of the flattened input `data`. If `axis` is a valid axis value, the
///   result has the same shape as `data` with `axis` removed and contains the
///   percentiles calculated along `axis`.
/// * `Err(ImgalError)`: If `axis >= data.ndim()`.
#[inline]
pub fn linear_percentile<'a, T, A, D>(
    data: A,
    percentile: f64,
    axis: Option<usize>,
    epsilon: Option<f64>,
    threads: Option<usize>,
) -> Result<ArrayD<f64>, 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" });
    }
    let per_arr = match axis {
        Some(ax) => {
            if ax >= data.ndim() {
                return Err(ImgalError::InvalidAxis {
                    axis_idx: ax,
                    dim_len: data.ndim(),
                });
            }
            let mut shape = data.shape().to_vec();
            shape.remove(ax);
            let mut arr = ArrayD::<f64>::zeros(IxDyn(&shape));
            // compute the percentile for each 1D lane along "axis"
            let lanes = data.lanes(Axis(ax));
            let lin_per_calc = |ln: ArrayView1<T>, pr: &mut f64| {
                let mut ln = Array::from_vec(ln.to_vec());
                *pr = linear_percentile_1d(ln.view_mut(), percentile, epsilon);
            };
            par!(threads,
                seq_exp: lanes.into_iter().zip(arr.iter_mut())
                    .for_each(|(ln, pr)| lin_per_calc(ln, pr)),
                par_exp: lanes.into_iter().zip(arr.iter_mut()).par_bridge()
                    .for_each(|(ln, pr)| lin_per_calc(ln, pr)));
            arr
        }
        None => {
            let mut arr = copy_into_flat(&data, threads);
            let per = linear_percentile_1d(arr.view_mut(), percentile, epsilon);
            Array::from_vec(vec![per]).into_dyn()
        }
    };
    Ok(per_arr)
}

/// 1D linear percentile.
///
/// The input data must be contiguous. If it is not then the input data is *not*
/// 1D and `0.0` is returned. It is up to the caller to ensure the input data is
/// contiguous
fn linear_percentile_1d<T>(mut data: ArrayViewMut1<T>, percentile: f64, epsilon: Option<f64>) -> f64
where
    T: AsNumeric,
{
    let epsilon = epsilon.unwrap_or(1e-12);
    // compute the percentile value using linear interpolation instead of
    // sorting the value array, get the "j" element via unstable selection if
    // "h" is an integer with epsilon value, return the percentile value
    let p_clamp = percentile.clamp(0.0, 100.0);
    match data.as_slice_mut() {
        Some(val_arr) => {
            let p = p_clamp / 100.0;
            let h = (val_arr.len() as f64 - 1.0) * p;
            let j = h.floor() as usize;
            let gamma = h - j as f64;
            if gamma.abs() < epsilon {
                val_arr
                    .select_nth_unstable_by(j, |a, b| a.partial_cmp(b).unwrap_or(Ordering::Less));
                return val_arr[j].to_f64();
            }
            val_arr.select_nth_unstable_by(j, |a, b| a.partial_cmp(b).unwrap_or(Ordering::Less));
            let v_j = val_arr[j].to_f64();
            val_arr
                .select_nth_unstable_by(j + 1, |a, b| a.partial_cmp(b).unwrap_or(Ordering::Less));
            let v_j1 = val_arr[j + 1].to_f64();
            (1.0 - gamma) * v_j + gamma * v_j1
        }
        None => 0.0,
    }
}