imgal 0.3.0

use std::collections::HashSet;

use ndarray::{Array2, ArrayBase, ArrayView1, AsArray, Axis, Ix1, Ix3, ViewRepr, Zip};

use crate::prelude::*;

/// Map G and S coordinates back to the input phasor array as a boolean mask.
///
/// # Description
///
/// Maps the G and S coordinates back to the input G/S phasor array and returns
/// a 2D boolean mask where `true` indicates G and S coordiantes representing
/// the `g_coords` and `s_coords` arrays.
///
/// # Arguments
///
/// * `data`: The G/S 3D array.
/// * `g_coords`: A 1D array of `g` coordinates in the `data` array. The
///   `g_coords` and `s_coords` array lengths must match.
/// * `s_coords`: A 1D array of `s` coordiantes in the `data` array. The
///   `s_coords` and `g_coords` array lengths must match.
/// * `axis`: The channel axis. If `None`, then `axis = 2`.
/// * `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(Array2<bool>)`: A 2-dimensional boolean mask where `true` pixels
///   represent values found in the `g_coords` and `s_coords` arrays.
/// * `Err(ImgalError)`: If `g_coords.len() != s_coords.len()`.
#[inline]
pub fn gs_mask<'a, T, A, B>(
    data: A,
    g_coords: B,
    s_coords: B,
    axis: Option<usize>,
    threads: Option<usize>,
) -> Result<Array2<bool>, ImgalError>
where
    A: AsArray<'a, T, Ix3>,
    B: AsArray<'a, f64, Ix1>,
    T: 'a + AsNumeric,
{
    let g_coords: ArrayBase<ViewRepr<&'a f64>, Ix1> = g_coords.into();
    let s_coords: ArrayBase<ViewRepr<&'a f64>, Ix1> = s_coords.into();
    let gl = g_coords.len();
    let sl = s_coords.len();
    if gl != sl {
        return Err(ImgalError::MismatchedArrayLengths {
            a_arr_name: "g_coords",
            a_arr_len: gl,
            b_arr_name: "s_coords",
            b_arr_len: sl,
        });
    }
    let a = axis.unwrap_or(2);
    if a >= 3 {
        return Err(ImgalError::InvalidAxis {
            axis_idx: a,
            dim_len: 3,
        });
    }
    // create a HashSet of G/S coordinates and check if a given G/S value pair
    // is within the set
    let data: ArrayBase<ViewRepr<&'a T>, Ix3> = data.into();
    let mut coords_set: HashSet<(u64, u64)> = HashSet::with_capacity(gl);
    g_coords.iter().zip(s_coords.iter()).for_each(|(g, s)| {
        coords_set.insert((g.to_bits(), s.to_bits()));
    });
    let mut shape = data.shape().to_vec();
    shape.remove(a);
    let mut map_arr = Array2::<bool>::default((shape[0], shape[1]));
    let lanes = data.lanes(Axis(a));
    let gs_mask_calc = |ln: ArrayView1<T>, p: &mut bool| {
        let dg = ln[0].to_f64();
        let ds = ln[1].to_f64();
        if (!dg.is_nan() || !ds.is_nan() || dg != 0.0 && ds != 0.0)
            && coords_set.contains(&(dg.to_bits(), ds.to_bits()))
        {
            *p = true;
        }
    };
    par!(threads,
        seq_exp: Zip::from(lanes).and(map_arr.view_mut())
            .for_each(&gs_mask_calc),
        par_exp: Zip::from(lanes).and(map_arr.view_mut())
            .par_for_each(&gs_mask_calc));
    Ok(map_arr)
}

/// Compute the modulation of phasor G and S coordinates.
///
/// # Description
///
/// Computes the modulation (M) of phasor G and S coordinates using the
/// pythagorean theorem to find the hypotenuse (*i.e.* the modulation):
///
/// ```text
/// M = √(G² + S²)
/// ````
///
/// # Arguments
///
/// * `g`: The real component, G.
/// * `s`: The imaginary component, S.
///
/// # Returns
///
/// * `f64`: The modulation (M) of the (G, S) phasor coordinates.
#[inline]
pub fn gs_modulation(g: f64, s: f64) -> f64 {
    let g_sqr: f64 = g * g;
    let s_sqr: f64 = s * s;
    (g_sqr + s_sqr).sqrt()
}

/// Compute the phase angle of phasor G and S coordinates.
///
/// # Description
///
/// Computes the phase angle or phi (φ) of phasor G and S coordinates using:
///
/// ```text
/// φ = tan⁻¹(S / G)
/// ````
///
/// This implementation uses atan2 and computes the four quadrant arctangent of
/// the phasor coordinates.
///
/// # Arguments
///
/// * `g`: The real component, G.
/// * `s`: The imaginary component, S.
///
/// # Returns
///
/// * `f64`: The phase (phi, φ)  of the (G, S) phasor coordinates.
#[inline]
pub fn gs_phase(g: f64, s: f64) -> f64 {
    s.atan2(g)
}

/// Compute the G and S coordinates for a monoexponential decay.
///
/// # Description
///
/// Computes the G and S coordinates for a monoexponential decay given as:
///
/// ```text
/// G = 1 / 1 + (ωτ)²
/// S = ωτ / 1 + (ωτ)²
/// ```
///
/// # Arguments
///
/// * `tau`: The lifetime of a monoexponential decay.
/// * `omega`: The angular frequency.
///
/// # Returns
///
/// * `(f64, f64)`: The monoexponential decay coordinates, (G, S).
///
/// # Reference
///
/// <https://doi.org/10.1117/1.JBO.25.7.071203>
#[inline]
pub fn monoexponential_coords(tau: f64, omega: f64) -> (f64, f64) {
    let ot = omega * tau;
    let denom = 1.0 + (ot * ot);
    let g = 1.0 / denom;
    let s = (omega * tau) / denom;
    (g, s)
}