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
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

use ndarray::{Array, ArrayBase, AsArray, Dimension, ViewRepr};

use crate::image::{histogram, histogram_bin_midpoint};
use crate::prelude::*;
use crate::statistics::min_max;
use crate::threshold::manual::manual_mask;

/// Create a boolean mask using Otsu's method.
///
/// # Description
///
/// Creates a boolean mask using Nobuyuki Otsu's automatic threshold method. The
/// Otsu threshold value used to create the mask is calculated by maximizing the
/// between-class variance of the assumed bimodal distribution in the image
/// histogram.
///
/// # Arguments
///
/// * `data`: The input n-dimensional image.
/// * `bins`: The number of bins to use to construct the image histogram for
///   Otsu's method. If `None`, then `bins = 256`.
/// * `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(Array<bool, D>)`: A boolean image of the same shape as the input image
///   with pixels that are greater than the computed Otsu threshold value set
///   as `true` and pixels that are below the Otsu threshold value set as
///   `false`.
/// * `Err(ImgalError)`: If `data.is_empty() == true` or `bins == 0`.
///
/// # Reference
///
/// <https://doi.org/10.1109/TSMC.1979.4310076>
#[inline]
pub fn otsu_mask<'a, T, A, D>(
    data: A,
    bins: Option<usize>,
    threads: Option<usize>,
) -> Result<Array<bool, D>, ImgalError>
where
    A: AsArray<'a, T, D>,
    D: Dimension,
    T: 'a + AsNumeric,
{
    let data: ArrayBase<ViewRepr<&'a T>, D> = data.into();
    let threshold = otsu_value(&data, bins, threads)?;
    Ok(manual_mask(data, threshold, threads))
}

/// Compute an image threshold with Otsu's method.
///
/// # Description
///
/// Calculates an image threshold value using Nobuyuki Otsu's automatic image
/// threshold method. The Otsu threshold value is calculated by maximizing the
/// between-class variance of the assumed bimodal distribution in the image
/// histogram.
///
/// # Arguments
///
/// * `data`: The input n-dimensional image.
/// * `bins`: The number of bins to use to construct the image histogram for
///   Otsu's method. If `None`, the `bins = 256`.
/// * `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(T)`: The Otsu threshold value.
/// * `Err(ImgalError)`: If `data.is_empty() == true` or `bins == 0`.
///
/// # Reference
///
/// <https://doi.org/10.1109/TSMC.1979.4310076>
#[inline]
pub fn otsu_value<'a, T, A, D>(
    data: A,
    bins: Option<usize>,
    threads: Option<usize>,
) -> Result<T, ImgalError>
where
    A: AsArray<'a, T, D>,
    D: Dimension,
    T: 'a + AsNumeric,
{
    let data: ArrayBase<ViewRepr<&'a T>, D> = data.into();
    let hist = histogram(&data, bins, threads)?;
    let dl = hist.len();
    let (min, max) = min_max(data, threads)?;
    let mut bcv: f64 = 0.0;
    let mut bcv_max: f64 = 0.0;
    let mut hist_sum: f64 = 0.0;
    let mut hist_inten: f64 = 0.0;
    let mut inten_k: f64 = 0.0;
    let mut k_star: usize = 0;
    let mut n_k: f64 = 0.0;
    hist.iter().enumerate().for_each(|(i, &v)| {
        let v = v as f64;
        hist_sum += v;
        hist_inten += i as f64 * v;
    });
    // compute threshold, here "k" is the current threshold at index "i"
    hist.iter().take(dl - 1).enumerate().for_each(|(i, &v)| {
        let v = v as f64;
        inten_k += i as f64 * v;
        n_k += v;
        let denom = n_k * (hist_sum - n_k);
        if denom != 0.0 {
            let num = (n_k / hist_sum) * hist_inten - inten_k;
            bcv = (num * num) / denom;
        } else {
            bcv = 0.0;
        }
        if bcv >= bcv_max {
            bcv_max = bcv;
            k_star = i;
        }
    });
    histogram_bin_midpoint(k_star, min, max, bins.unwrap_or(256))
}