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
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269

use ndarray::{Array2, Axis, s};

use imgal::parameter::omega;
use imgal::phasor::calibration::{
    calibrate_coords, calibrate_gs_image, calibrate_gs_image_mut, modulation_and_phase,
};
use imgal::phasor::plot::{gs_mask, gs_modulation, gs_phase, monoexponential_coords};
use imgal::phasor::time_domain::{gs_image, imaginary_coord, real_coord};
use imgal::prelude::*;
use imgal::simulation::decay::{gaussian_exponential_decay_3d, ideal_exponential_decay_1d};
use imgal::simulation::noise::poisson_noise_mut;

const TOLERANCE: f64 = 1e-10;
const SAMPLES: usize = 256;
const PERIOD: f64 = 12.5;
const TAUS: [f64; 2] = [1.0, 3.0];
const FRACTIONS: [f64; 2] = [0.7, 0.3];
const TOTAL_COUNTS: f64 = 5000.0;
const IRF_CENTER: f64 = 3.0;
const IRF_WIDTH: f64 = 0.5;
const SHAPE: (usize, usize) = (10, 10);
const MODULATION: f64 = 0.7;
const PHASE: f64 = -0.981;
const THREADS: Option<usize> = Some(0);

fn approx_equal(a: f64, b: f64) -> bool {
    (a - b).abs() < TOLERANCE
}

fn get_circle_mask(shape: (usize, usize), center: (isize, isize), radius: isize) -> Array2<bool> {
    let (row, col) = shape;
    let (cx, cy) = center;
    let r2 = radius * radius;
    let y_min = (cy - radius).max(0);
    let y_max = (cy + radius).min(row as isize - 1);
    let x_min = (cx - radius).max(0);
    let x_max = (cx + radius).min(col as isize - 1);
    let mut mask = Array2::<bool>::default(shape);
    for y in y_min..=y_max {
        for x in x_min..=x_max {
            let dx = cx - x;
            let dy = cy - y;
            // use the squared distance formula for a quick circle mask
            if dx * dx + dy * dy <= r2 {
                mask[[y as usize, x as usize]] = true;
            }
        }
    }
    mask
}

/// Tests that `calibrate_coords` returns the expected calibrated G and S
/// values.
#[test]
fn calibration_calibrate_coords_expected_results() {
    let g = -0.37;
    let s = 0.68;
    let coords_cal = calibrate_coords(g, s, MODULATION, PHASE);
    assert!(approx_equal(coords_cal.0, 0.2515280246));
    assert!(approx_equal(coords_cal.1, 0.4799902632));
}

/// Tests that `calibrate_gs_image` returns the expected calibrated G/S values
/// in a new image array.
#[test]
fn calibration_calibrate_gs_image_expected_results() -> Result<(), ImgalError> {
    let data = gaussian_exponential_decay_3d(
        SAMPLES,
        PERIOD,
        &TAUS,
        &FRACTIONS,
        TOTAL_COUNTS,
        IRF_CENTER,
        IRF_WIDTH,
        SHAPE,
        None,
    )?;
    let gs_arr = gs_image(data.view(), PERIOD, None, None, None, None)?;
    let cal_gs_arr_par = calibrate_gs_image(gs_arr.view(), MODULATION, PHASE, None, THREADS);
    let cal_gs_arr_seq = calibrate_gs_image(gs_arr.view(), MODULATION, PHASE, None, None);
    let g_mean_par = cal_gs_arr_par.index_axis(Axis(2), 0).mean().unwrap();
    let g_mean_seq = cal_gs_arr_seq.index_axis(Axis(2), 0).mean().unwrap();
    let s_mean_par = cal_gs_arr_par.index_axis(Axis(2), 1).mean().unwrap();
    let s_mean_seq = cal_gs_arr_seq.index_axis(Axis(2), 1).mean().unwrap();
    assert!(approx_equal(cal_gs_arr_par[[5, 5, 0]], 0.2536762376));
    assert!(approx_equal(cal_gs_arr_seq[[5, 5, 0]], 0.2536762376));
    assert!(approx_equal(cal_gs_arr_par[[5, 5, 1]], 0.4819949555));
    assert!(approx_equal(cal_gs_arr_seq[[5, 5, 1]], 0.4819949555));
    assert!(approx_equal(g_mean_par, 0.2536762376));
    assert!(approx_equal(g_mean_seq, 0.2536762376));
    assert!(approx_equal(s_mean_par, 0.4819949555));
    assert!(approx_equal(s_mean_seq, 0.4819949555));
    Ok(())
}

/// Tests that `calibrate_gs_image_mut` mutates the input data with the expected
/// G/S calibrated values.
#[test]
fn calibration_calibrate_gs_image_mut_expected_results() -> Result<(), ImgalError> {
    let data = gaussian_exponential_decay_3d(
        SAMPLES,
        PERIOD,
        &TAUS,
        &FRACTIONS,
        TOTAL_COUNTS,
        IRF_CENTER,
        IRF_WIDTH,
        SHAPE,
        None,
    )?;
    let mut gs_arr_par = gs_image(data.view(), PERIOD, None, None, None, None)?;
    let mut gs_arr_seq = gs_arr_par.clone();
    calibrate_gs_image_mut(gs_arr_par.view_mut(), MODULATION, PHASE, None, THREADS);
    calibrate_gs_image_mut(gs_arr_seq.view_mut(), MODULATION, PHASE, None, None);
    let g_mean_par = gs_arr_par.index_axis(Axis(2), 0).mean().unwrap();
    let g_mean_seq = gs_arr_seq.index_axis(Axis(2), 0).mean().unwrap();
    let s_mean_par = gs_arr_par.index_axis(Axis(2), 1).mean().unwrap();
    let s_mean_seq = gs_arr_seq.index_axis(Axis(2), 1).mean().unwrap();
    assert!(approx_equal(gs_arr_par[[5, 5, 0]], 0.2536762376));
    assert!(approx_equal(gs_arr_seq[[5, 5, 0]], 0.2536762376));
    assert!(approx_equal(gs_arr_par[[5, 5, 1]], 0.4819949555));
    assert!(approx_equal(gs_arr_seq[[5, 5, 1]], 0.4819949555));
    assert!(approx_equal(g_mean_par, 0.2536762376));
    assert!(approx_equal(g_mean_seq, 0.2536762376));
    assert!(approx_equal(s_mean_par, 0.4819949555));
    assert!(approx_equal(s_mean_seq, 0.4819949555));
    Ok(())
}

/// Tests that `modulation_and_phase` returns the expected modulation and phase
/// values for the given parameters.
#[test]
fn calibration_modulation_and_phase_expected_results() {
    let w = omega(PERIOD);
    let mod_phs = modulation_and_phase(-0.055, 0.59, 1.1, w);
    assert!(approx_equal(mod_phs.0, 1.4768757234));
    assert!(approx_equal(mod_phs.1, -1.1586655116));
}

/// Tests that `gs_mask` maps G and S coordinates back to the original input
/// image as a boolean mask.
#[test]
fn plot_gs_mask_expected_results() -> Result<(), ImgalError> {
    let mut data = gaussian_exponential_decay_3d(
        SAMPLES,
        PERIOD,
        &TAUS,
        &FRACTIONS,
        TOTAL_COUNTS,
        IRF_CENTER,
        IRF_WIDTH,
        (50, 50),
        None,
    )?;
    poisson_noise_mut(data.view_mut().into_dyn(), 0.3, None, None);
    let gs_arr = gs_image(data.view(), PERIOD, None, None, None, None)?;
    let g_coords = gs_arr.slice(s![25..30, 25..30, 0]).flatten().to_vec();
    let s_coords = gs_arr.slice(s![25..30, 25..30, 1]).flatten().to_vec();
    let mask_par = gs_mask(gs_arr.view(), &g_coords, &s_coords, None, THREADS)?;
    let mask_seq = gs_mask(gs_arr.view(), &g_coords, &s_coords, None, None)?;
    assert_eq!(mask_par[[28, 28]], true);
    assert_eq!(mask_par[[5, 5]], false);
    assert_eq!(mask_seq[[28, 28]], true);
    assert_eq!(mask_seq[[5, 5]], false);
    Ok(())
}

/// Tests that `gs_modulation` returns the expected modulation for a G and S
/// pair.
#[test]
fn plot_gs_modulation_expected_results() {
    let m = gs_modulation(0.71, 0.43);
    assert!(approx_equal(m, 0.8300602387));
}

/// Tests that `gs_phase` returns the expected phase for a G and S pair.
#[test]
fn plot_gs_phase_expected_results() {
    let p = gs_phase(0.71, 0.43);
    assert!(approx_equal(p, 0.5445517081));
}

/// Tests that `monoexponential_coords` returns the expected G and S values for
/// the given tau and omega values.
#[test]
fn plot_monoexponential_coords_expected_results() {
    let w = omega(PERIOD);
    let coords = monoexponential_coords(1.1, w);
    assert!(approx_equal(coords.0, 0.765860473));
    assert!(approx_equal(coords.1, 0.4234598078));
}

/// Tests that `gs_image` returns the expected G/S phasor image by checking
/// points inside the image (with and without a mask) and the mean of each
/// channel.
#[test]
fn time_domain_gs_image_expected_results() -> Result<(), ImgalError> {
    let data = gaussian_exponential_decay_3d(
        SAMPLES,
        PERIOD,
        &TAUS,
        &FRACTIONS,
        TOTAL_COUNTS,
        IRF_CENTER,
        IRF_WIDTH,
        (100, 100),
        None,
    )?;
    let mask = get_circle_mask((100, 100), (50, 50), 8);
    let gs_no_mask_par = gs_image(data.view(), PERIOD, None, None, None, THREADS)?;
    let gs_no_mask_seq = gs_image(data.view(), PERIOD, None, None, None, None)?;
    let gs_with_mask_par = gs_image(data.view(), PERIOD, Some(mask.view()), None, None, THREADS)?;
    let gs_with_mask_seq = gs_image(data.view(), PERIOD, Some(mask.view()), None, None, None)?;
    let g_no_mask_view_par = gs_no_mask_par.index_axis(Axis(2), 0);
    let g_no_mask_view_seq = gs_no_mask_seq.index_axis(Axis(2), 0);
    let s_no_mask_view_par = gs_no_mask_par.index_axis(Axis(2), 1);
    let s_no_mask_view_seq = gs_no_mask_seq.index_axis(Axis(2), 1);
    let g_with_mask_view_par = gs_with_mask_par.index_axis(Axis(2), 0);
    let g_with_mask_view_seq = gs_with_mask_seq.index_axis(Axis(2), 0);
    let s_with_mask_view_par = gs_with_mask_par.index_axis(Axis(2), 1);
    let s_with_mask_view_seq = gs_with_mask_seq.index_axis(Axis(2), 1);
    assert!(approx_equal(
        g_no_mask_view_par.mean().unwrap(),
        -0.3706731273
    ));
    assert!(approx_equal(
        g_no_mask_view_seq.mean().unwrap(),
        -0.3706731273
    ));
    assert!(approx_equal(
        s_no_mask_view_par.mean().unwrap(),
        0.6841432489
    ));
    assert!(approx_equal(
        s_no_mask_view_seq.mean().unwrap(),
        0.6841432489
    ));
    assert!(approx_equal(g_with_mask_view_par[[45, 52]], -0.3706731273));
    assert!(approx_equal(g_with_mask_view_seq[[45, 52]], -0.3706731273));
    assert!(approx_equal(s_with_mask_view_par[[45, 52]], 0.6841432489));
    assert!(approx_equal(s_with_mask_view_seq[[45, 52]], 0.6841432489));
    assert_eq!(g_with_mask_view_par[[5, 8]], 0.0);
    assert_eq!(g_with_mask_view_seq[[5, 8]], 0.0);
    assert_eq!(s_with_mask_view_par[[5, 8]], 0.0);
    assert_eq!(s_with_mask_view_seq[[5, 8]], 0.0);
    Ok(())
}

/// Tests that `imaginary_coord` returns the expected imaginary (S) coordinate.
#[test]
fn time_domain_imaginary_coord_expected_results() -> Result<(), ImgalError> {
    let data = ideal_exponential_decay_1d(SAMPLES, PERIOD, &TAUS, &FRACTIONS, TOTAL_COUNTS, None)?;
    let s_coord_par = imaginary_coord(&data, PERIOD, None, THREADS);
    let s_coord_seq = imaginary_coord(&data, PERIOD, None, None);
    assert!(approx_equal(s_coord_par, 0.410217863));
    assert!(approx_equal(s_coord_seq, 0.410217863));
    Ok(())
}

/// Tests that `real_coord` returns the expected real (G) coordinate.
#[test]
fn time_domain_real_coord_expected_results() -> Result<(), ImgalError> {
    let data = ideal_exponential_decay_1d(SAMPLES, PERIOD, &TAUS, &FRACTIONS, TOTAL_COUNTS, None)?;
    let g_coord_par = real_coord(&data, PERIOD, None, THREADS);
    let g_coord_seq = real_coord(&data, PERIOD, None, None);
    assert!(approx_equal(g_coord_par, 0.660137605));
    assert!(approx_equal(g_coord_seq, 0.660137605));
    Ok(())
}