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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171

use ndarray::{Array1, ArrayBase, AsArray, Ix1, ViewRepr, Zip};
use rayon::prelude::*;
use rustfft::{FftPlanner, num_complex::Complex, num_traits::Zero};

use crate::prelude::*;

/// Convolve two 1D signals using the Fast Fourier Transform (FFT).
///
/// # Description
///
/// Computes the convolution of two discrete signals (`data_a` and `data_b`) by
/// transforming them into the frequency domain, multiplying them, and then
/// transforming the result back into a signal. This function uses "same-length"
/// trimming with the first parameter `data_a`. This means that the returned
/// convolution's array length will have the same length as `data_a`.
///
/// # Arguments
///
/// * `data_a`: The first input signal to FFT convolve. Returned convolution
///   arrays will be "same-length" trimmed to `data_a`'s length.
/// * `data_b`: The second input signal to FFT convolve.
/// * `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
///
/// * `Array1<f64>`: The FFT convolved result of the same length as input signal
///   `data_a`.
#[inline]
pub fn fft_convolve_1d<'a, T, A>(data_a: A, data_b: A, threads: Option<usize>) -> Array1<f64>
where
    A: AsArray<'a, T, Ix1>,
    T: 'a + AsNumeric,
{
    let data_a: ArrayBase<ViewRepr<&'a T>, Ix1> = data_a.into();
    let data_b: ArrayBase<ViewRepr<&'a T>, Ix1> = data_b.into();
    let n_a = data_a.len();
    let n_b = data_b.len();
    let n_fft = n_a + n_b - 1;
    let fft_size = n_fft.next_power_of_two();
    let mut a_fft_buf = vec![Complex::zero(); fft_size];
    let mut b_fft_buf = vec![Complex::zero(); fft_size];
    let load_buffers = |a_buf: &mut Complex<f64>, b_buf: &mut Complex<f64>, a: &T, b: &T| {
        *a_buf = Complex::new(a.to_f64(), 0.0);
        *b_buf = Complex::new(b.to_f64(), 0.0);
    };
    par!(threads,
        seq_exp: Zip::from(&mut a_fft_buf[..n_a]).and(&mut b_fft_buf[..n_b])
            .and(data_a.view())
            .and(data_b.view())
            .for_each(&load_buffers),
        par_exp: Zip::from(&mut a_fft_buf[..n_a]).and(&mut b_fft_buf[..n_b])
            .and(data_a.view())
            .and(data_b.view())
            .par_for_each(&load_buffers));
    let mut planner = FftPlanner::new();
    let fft = planner.plan_fft_forward(fft_size);
    let ifft = planner.plan_fft_inverse(fft_size);
    fft.process(&mut a_fft_buf);
    fft.process(&mut b_fft_buf);
    // multiply in the frequency domain and extract the real component (scaled
    // and input length trimmed)
    let mul_calc = |a: &mut Complex<f64>, b: &Complex<f64>| {
        *a *= b;
    };
    par!(threads,
        seq_exp: a_fft_buf.iter_mut().zip(b_fft_buf.iter())
            .for_each(|(a, b)| mul_calc(a, b)),
        par_exp: a_fft_buf.par_iter_mut().zip(b_fft_buf.par_iter())
            .for_each(|(a, b)| mul_calc(a, b)));
    ifft.process(&mut a_fft_buf);
    let scale = 1.0 / fft_size as f64;
    par!(threads,
        seq_exp: (0..n_a).zip(a_fft_buf.iter()).map(|(_, v)| v.re * scale)
            .collect::<Array1<f64>>(),
        par_exp: Array1::from_vec((0..n_a).into_par_iter()
            .zip(a_fft_buf.par_iter())
            .map(|(_, v)| v.re * scale)
            .collect::<Vec<f64>>()))
}

/// Deconvolve two 1D signals using the Fast Fourier Transform (FFT).
///
/// # Description
///
/// Computes the deconvolution of two discrete signals (`data_a` and `data_b`)
/// by transforming them into the frequency domain, dividing them, and then
/// transforming the result back into a signal. This function uses "same-length"
/// trimming with the first parameter `data_a`. This means that the returned
/// deconvolution's array length will have the same length as `data_a`.
///
/// # Arguments
///
/// * `data_a`: The first input signal to FFT deconvolve. Returned deconvolution
///   arrays will be "same-length" trimmed to `data_a`'s length.
/// * `data_b`: The second input singal to FFT deconvolve.
/// * `epsilon`: An epsilon value to prevent division by zero errors (default =
///   `1e-8`).
/// * `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
///
/// * `Array1<f64>`: The FFT deconvolved result of the same length as input
///   signal `data_a`.
#[inline]
pub fn fft_deconvolve_1d<'a, T, A>(
    data_a: A,
    data_b: A,
    epsilon: Option<f64>,
    threads: Option<usize>,
) -> Array1<f64>
where
    A: AsArray<'a, T, Ix1>,
    T: 'a + AsNumeric,
{
    let data_a: ArrayBase<ViewRepr<&'a T>, Ix1> = data_a.into();
    let data_b: ArrayBase<ViewRepr<&'a T>, Ix1> = data_b.into();
    let epsilon = epsilon.unwrap_or(1e-8);
    let n_a = data_a.len();
    let n_b = data_b.len();
    let n_fft = n_a + n_b - 1;
    let fft_size = n_fft.next_power_of_two();
    let mut a_fft_buf = vec![Complex::zero(); fft_size];
    let mut b_fft_buf = vec![Complex::zero(); fft_size];
    let load_buffers = |a_buf: &mut Complex<f64>, b_buf: &mut Complex<f64>, a: &T, b: &T| {
        *a_buf = Complex::new(a.to_f64(), 0.0);
        *b_buf = Complex::new(b.to_f64(), 0.0);
    };
    par!(threads,
        seq_exp: Zip::from(&mut a_fft_buf[..n_a]).and(&mut b_fft_buf[..n_b])
            .and(data_a.view())
            .and(data_b.view())
            .for_each(&load_buffers),
        par_exp: Zip::from(&mut a_fft_buf[..n_a]).and(&mut b_fft_buf[..n_b])
            .and(data_a.view())
            .and(data_b.view())
            .par_for_each(&load_buffers));
    let mut planner = FftPlanner::new();
    let fft = planner.plan_fft_forward(fft_size);
    let ifft = planner.plan_fft_inverse(fft_size);
    fft.process(&mut a_fft_buf);
    fft.process(&mut b_fft_buf);
    // divide in the frequency domain with epsilon value and extract the real
    // component (scaled and input length trimmed)
    let div_calc = |a: &mut Complex<f64>, b: &Complex<f64>| {
        if a.norm_sqr() > epsilon {
            *a /= b;
        } else {
            *a = Complex::zero();
        }
    };
    par!(threads,
        seq_exp: a_fft_buf.iter_mut().zip(b_fft_buf.iter())
            .for_each(|(a, b)| div_calc(a, b)),
        par_exp: a_fft_buf.par_iter_mut().zip(b_fft_buf.par_iter())
            .for_each(|(a, b)| div_calc(a, b)));
    ifft.process(&mut a_fft_buf);
    let scale = 1.0 / fft_size as f64;
    par!(threads,
        seq_exp: (0..n_a).zip(a_fft_buf.iter()).map(|(_, v)| v.re * scale)
            .collect::<Array1<f64>>(),
        par_exp: Array1::from_vec((0..n_a).into_par_iter()
            .zip(a_fft_buf.par_iter())
            .map(|(_, v)| v.re * scale)
            .collect::<Vec<f64>>()))
}