1pub const COLOC_BINS: usize = 1024;
2pub const COLOC_SCATTER: usize = 128;
3
4#[derive(Clone, Copy, Debug, PartialEq, Eq)]
5pub struct ColocPair {
6 a: usize,
7 b: usize,
8}
9
10impl ColocPair {
11 pub fn new(a: usize, b: usize) -> Self {
12 assert!(a != b, "a colocalization pair names two distinct channels");
13 Self { a, b }
14 }
15
16 pub fn a(self) -> usize {
17 self.a
18 }
19
20 pub fn b(self) -> usize {
21 self.b
22 }
23}
24
25#[derive(Clone, Copy, Debug, PartialEq)]
26pub struct ColocThresholds {
27 a: f32,
28 b: f32,
29}
30
31impl ColocThresholds {
32 pub fn new(a: f32, b: f32) -> Self {
33 assert!(
34 (0.0..=1.0).contains(&a) && (0.0..=1.0).contains(&b),
35 "a colocalization threshold is a unit fraction"
36 );
37 Self { a, b }
38 }
39
40 pub fn a(self) -> f32 {
41 self.a
42 }
43
44 pub fn b(self) -> f32 {
45 self.b
46 }
47
48 fn bins(self) -> [usize; 2] {
49 [self.a, self.b]
50 .map(|threshold| ((threshold * COLOC_BINS as f32) as usize).min(COLOC_BINS - 1))
51 }
52}
53
54#[derive(Clone, Copy, Debug, Default, PartialEq)]
55pub struct ColocMoments {
56 voxels: u64,
57 sum_a: f64,
58 sum_b: f64,
59 sum_aa: f64,
60 sum_bb: f64,
61 sum_ab: f64,
62}
63
64impl ColocMoments {
65 pub fn new(voxels: u64, sums: [f64; 5]) -> Self {
66 let [sum_a, sum_b, sum_aa, sum_bb, sum_ab] = sums;
67 assert!(
68 sums.iter().all(|value| value.is_finite() && *value >= 0.0),
69 "colocalization moments are finite and not negative"
70 );
71 Self {
72 voxels,
73 sum_a,
74 sum_b,
75 sum_aa,
76 sum_bb,
77 sum_ab,
78 }
79 }
80
81 pub fn voxels(self) -> u64 {
82 self.voxels
83 }
84
85 pub fn sums(self) -> [f64; 5] {
86 [
87 self.sum_a,
88 self.sum_b,
89 self.sum_aa,
90 self.sum_bb,
91 self.sum_ab,
92 ]
93 }
94
95 pub fn push(&mut self, a: f64, b: f64) {
96 self.voxels += 1;
97 self.sum_a += a;
98 self.sum_b += b;
99 self.sum_aa += a * a;
100 self.sum_bb += b * b;
101 self.sum_ab += a * b;
102 }
103
104 pub fn merge(&mut self, other: &Self) {
105 self.voxels += other.voxels;
106 self.sum_a += other.sum_a;
107 self.sum_b += other.sum_b;
108 self.sum_aa += other.sum_aa;
109 self.sum_bb += other.sum_bb;
110 self.sum_ab += other.sum_ab;
111 }
112
113 fn pearson(self) -> Option<f64> {
114 let count = self.voxels as f64;
115 let left = count * self.sum_aa - self.sum_a * self.sum_a;
116 let right = count * self.sum_bb - self.sum_b * self.sum_b;
117 if left <= 0.0 || right <= 0.0 {
118 return None;
119 }
120 let numerator = count * self.sum_ab - self.sum_a * self.sum_b;
121 Some((numerator / (left * right).sqrt()).clamp(-1.0, 1.0))
122 }
123
124 fn overlap(self) -> Option<f64> {
125 let mass = self.sum_aa * self.sum_bb;
126 if mass <= 0.0 {
127 return None;
128 }
129 Some((self.sum_ab / mass.sqrt()).clamp(0.0, 1.0))
130 }
131}
132
133#[derive(Clone, Copy, Debug, PartialEq)]
134pub struct ColocStatistics {
135 voxels: u64,
136 pearson: Option<f64>,
137 overlap: Option<f64>,
138 manders: [Option<f64>; 2],
139}
140
141impl ColocStatistics {
142 pub fn voxels(self) -> u64 {
143 self.voxels
144 }
145
146 pub fn pearson(self) -> Option<f64> {
147 self.pearson
148 }
149
150 pub fn overlap(self) -> Option<f64> {
151 self.overlap
152 }
153
154 pub fn manders(self) -> [Option<f64>; 2] {
155 self.manders
156 }
157}
158
159#[derive(Clone, Debug, PartialEq)]
160pub struct ColocFrame {
161 moments: ColocMoments,
162 a_by_b: Vec<f64>,
163 b_by_a: Vec<f64>,
164}
165
166impl ColocFrame {
167 pub fn new() -> Self {
168 Self {
169 moments: ColocMoments::default(),
170 a_by_b: vec![0.0; COLOC_BINS],
171 b_by_a: vec![0.0; COLOC_BINS],
172 }
173 }
174
175 pub fn from_parts(moments: ColocMoments, a_by_b: Vec<f64>, b_by_a: Vec<f64>) -> Self {
176 assert!(
177 a_by_b.len() == COLOC_BINS && b_by_a.len() == COLOC_BINS,
178 "a colocalization frame carries one sum per intensity bin"
179 );
180 Self {
181 moments,
182 a_by_b,
183 b_by_a,
184 }
185 }
186
187 pub fn moments(&self) -> ColocMoments {
188 self.moments
189 }
190
191 pub fn a_by_b(&self) -> &[f64] {
192 &self.a_by_b
193 }
194
195 pub fn b_by_a(&self) -> &[f64] {
196 &self.b_by_a
197 }
198
199 pub fn push(&mut self, a: f32, b: f32) {
200 let bin_a = bin(a);
201 let bin_b = bin(b);
202 self.moments.push(f64::from(a), f64::from(b));
203 self.a_by_b[bin_b] += f64::from(a);
204 self.b_by_a[bin_a] += f64::from(b);
205 }
206
207 pub fn merge(&mut self, other: &Self) {
208 self.moments.merge(&other.moments);
209 for (target, source) in self.a_by_b.iter_mut().zip(&other.a_by_b) {
210 *target += source;
211 }
212 for (target, source) in self.b_by_a.iter_mut().zip(&other.b_by_a) {
213 *target += source;
214 }
215 }
216
217 pub fn statistics(&self, thresholds: ColocThresholds) -> ColocStatistics {
218 let [cut_a, cut_b] = thresholds.bins();
219 let a_mass = self.moments.sum_a;
220 let b_mass = self.moments.sum_b;
221 let first = (a_mass > 0.0).then(|| (above(&self.a_by_b, cut_b) / a_mass).clamp(0.0, 1.0));
222 let second = (b_mass > 0.0).then(|| (above(&self.b_by_a, cut_a) / b_mass).clamp(0.0, 1.0));
223 ColocStatistics {
224 voxels: self.moments.voxels,
225 pearson: self.moments.pearson(),
226 overlap: self.moments.overlap(),
227 manders: [first, second],
228 }
229 }
230}
231
232impl Default for ColocFrame {
233 fn default() -> Self {
234 Self::new()
235 }
236}
237
238#[derive(Clone, Debug, PartialEq)]
239pub struct ColocSeries {
240 pair: ColocPair,
241 frames: Vec<ColocFrame>,
242 counts_a: Vec<u64>,
243 counts_b: Vec<u64>,
244 scatter: Vec<u64>,
245 total: ColocFrame,
246}
247
248impl ColocSeries {
249 pub fn new(
250 pair: ColocPair,
251 frames: Vec<ColocFrame>,
252 counts_a: Vec<u64>,
253 counts_b: Vec<u64>,
254 scatter: Vec<u64>,
255 ) -> Self {
256 let mut total = ColocFrame::new();
257 for frame in &frames {
258 total.merge(frame);
259 }
260 Self::from_parts(pair, frames, counts_a, counts_b, scatter, total)
261 }
262
263 pub fn from_parts(
264 pair: ColocPair,
265 frames: Vec<ColocFrame>,
266 counts_a: Vec<u64>,
267 counts_b: Vec<u64>,
268 scatter: Vec<u64>,
269 total: ColocFrame,
270 ) -> Self {
271 assert!(
272 counts_a.len() == COLOC_BINS && counts_b.len() == COLOC_BINS,
273 "a colocalization series counts every intensity bin"
274 );
275 assert!(
276 scatter.len() == COLOC_SCATTER * COLOC_SCATTER,
277 "a colocalization series maps the square scatter grid"
278 );
279 Self {
280 pair,
281 frames,
282 counts_a,
283 counts_b,
284 scatter,
285 total,
286 }
287 }
288
289 pub fn pair(&self) -> ColocPair {
290 self.pair
291 }
292
293 pub fn frames(&self) -> &[ColocFrame] {
294 &self.frames
295 }
296
297 pub fn frame(&self, time: u32) -> Option<&ColocFrame> {
298 self.frames.get(time as usize)
299 }
300
301 pub fn counts_a(&self) -> &[u64] {
302 &self.counts_a
303 }
304
305 pub fn counts_b(&self) -> &[u64] {
306 &self.counts_b
307 }
308
309 pub fn scatter(&self) -> &[u64] {
310 &self.scatter
311 }
312
313 pub fn total(&self) -> &ColocFrame {
314 &self.total
315 }
316
317 pub fn scatter_at(&self, a: usize, b: usize) -> u64 {
318 assert!(
319 a < COLOC_SCATTER && b < COLOC_SCATTER,
320 "the scatter cell lies on the grid"
321 );
322 self.scatter[a + COLOC_SCATTER * b]
323 }
324
325 pub fn statistics(&self, thresholds: ColocThresholds) -> ColocStatistics {
326 self.total.statistics(thresholds)
327 }
328
329 pub fn frame_statistics(
330 &self,
331 time: u32,
332 thresholds: ColocThresholds,
333 ) -> Option<ColocStatistics> {
334 self.frame(time).map(|frame| frame.statistics(thresholds))
335 }
336
337 pub fn otsu(&self) -> ColocThresholds {
338 ColocThresholds::new(otsu(&self.counts_a), otsu(&self.counts_b))
339 }
340
341 pub fn csv(&self, thresholds: ColocThresholds) -> String {
342 let mut text = String::from("scope,frame,voxels,pearson,overlap,manders_a,manders_b\n");
343 let mut rows = vec![(
344 "series".to_string(),
345 String::new(),
346 self.statistics(thresholds),
347 )];
348 for (time, frame) in self.frames.iter().enumerate() {
349 rows.push((
350 "frame".to_string(),
351 time.to_string(),
352 frame.statistics(thresholds),
353 ));
354 }
355 for (scope, frame, statistics) in rows {
356 text.push_str(&format!(
357 "{scope},{frame},{},{},{},{},{}\n",
358 statistics.voxels(),
359 option(statistics.pearson()),
360 option(statistics.overlap()),
361 option(statistics.manders()[0]),
362 option(statistics.manders()[1])
363 ));
364 }
365 text
366 }
367}
368
369fn bin(unit: f32) -> usize {
370 ((unit.clamp(0.0, 1.0) * COLOC_BINS as f32) as usize).min(COLOC_BINS - 1)
371}
372
373fn above(sums: &[f64], cut: usize) -> f64 {
374 sums.get(cut + 1..).map_or(0.0, |tail| tail.iter().sum())
375}
376
377fn otsu(counts: &[u64]) -> f32 {
378 let total: f64 = counts.iter().map(|count| *count as f64).sum();
379 if total <= 0.0 {
380 return 0.5;
381 }
382 let mut sum_all = 0.0;
383 for (index, count) in counts.iter().enumerate() {
384 sum_all += (index as f64 + 0.5) * (*count as f64);
385 }
386 let mut weight = 0.0;
387 let mut sum = 0.0;
388 let mut best = 0.0;
389 let mut split = 0;
390 for (index, count) in counts.iter().enumerate() {
391 weight += *count as f64;
392 sum += (index as f64 + 0.5) * (*count as f64);
393 let remaining = total - weight;
394 if weight <= 0.0 || remaining <= 0.0 {
395 continue;
396 }
397 let between = weight * remaining * (sum / weight - (sum_all - sum) / remaining).powi(2);
398 if between >= best {
399 best = between;
400 split = index;
401 }
402 }
403 split as f32 / counts.len() as f32
404}
405
406fn option(value: Option<f64>) -> String {
407 value.map_or_else(String::new, |value| format!("{value:.6}"))
408}