1use std::ops::Deref;
5
6use num::{FromPrimitive, ToPrimitive};
7
8use crate::cv::features::{GLCM, glcm_multichannel, glcm_multichannel_object};
9use crate::im::PineappleViewBuffer;
10
11#[inline]
12pub fn texture_energy(glcm: &GLCM) -> f32 {
13 glcm.iter().fold(0.0, |acc, (_, _, x)| acc + x * x)
14}
15
16#[inline]
17pub fn texture_contrast(glcm: &GLCM) -> f32 {
18 let mut contrast = 0.0;
19 for (i, j, g_ij) in glcm.iter() {
20 let i_minus_j = i as f32 - j as f32;
21 contrast += i_minus_j * i_minus_j * g_ij;
22 }
23 contrast
24}
25
26#[inline]
27pub fn texture_correlation(glcm: &GLCM) -> f32 {
28 let (px, py) = glcm.margin_sums();
29
30 let ux = px
31 .iter()
32 .enumerate()
33 .fold(0.0, |acc, (i, &x)| acc + ((i as f32 + 1.0) * x));
34
35 let uy = py
36 .iter()
37 .enumerate()
38 .fold(0.0, |acc, (i, &y)| acc + ((i as f32 + 1.0) * y));
39
40 let sx = px
41 .iter()
42 .enumerate()
43 .fold(0.0, |acc, (i, &x)| {
44 acc + (i as f32 + 1.0 - ux) * (i as f32 + 1.0 - ux) * x
45 })
46 .sqrt();
47
48 let sy = py
49 .iter()
50 .enumerate()
51 .fold(0.0, |acc, (i, &y)| {
52 acc + (i as f32 + 1.0 - uy) * (i as f32 + 1.0 - uy) * y
53 })
54 .sqrt();
55
56 let mut correlation = 0.0;
57 for (i, j, g_ij) in glcm.iter() {
58 correlation += ((i as f32 + 1.0 - ux) * (j as f32 + 1.0 - uy) * g_ij) / (sx * sy);
59 }
60
61 correlation
62}
63
64#[inline]
65pub fn texture_sum_of_squares(glcm: &GLCM) -> f32 {
66 let (px, _) = glcm.margin_sums();
67
68 let u = px
69 .iter()
70 .enumerate()
71 .fold(0.0, |acc, (i, &x)| acc + ((i as f32 + 1.0) * x));
72
73 let mut sum_of_squares = 0.0;
74 for (i, _, g_ij) in glcm.iter() {
75 sum_of_squares += (i as f32 + 1.0 - u) * (i as f32 + 1.0 - u) * g_ij;
76 }
77
78 sum_of_squares
79}
80
81#[inline]
82pub fn texture_inverse_difference(glcm: &GLCM) -> f32 {
83 let mut inverse_difference_moment = 0.0;
84 for (i, j, g_ij) in glcm.iter() {
85 inverse_difference_moment +=
86 (1.0 / (1.0 + ((i as f32 - j as f32) * (i as f32 - j as f32)))) * g_ij;
87 }
88
89 inverse_difference_moment
90}
91
92#[inline]
93pub fn texture_sum_average(glcm: &GLCM) -> f32 {
94 let mut px_plus_y = vec![0.0; 2 * glcm.rows()];
95 for (i, j, g_ij) in glcm.iter() {
96 px_plus_y[i + j] += g_ij
97 }
98
99 let mut sum_average = 0.0;
100 for (i, px_plus_y_k) in px_plus_y.iter().enumerate().take(2 * glcm.rows()) {
101 sum_average += i as f32 * px_plus_y_k;
102 }
103
104 sum_average
105}
106
107#[inline]
108pub fn texture_sum_variance(glcm: &GLCM) -> f32 {
109 let mut px_plus_y = vec![0.0; 2 * glcm.rows()];
110 for (i, j, g_ij) in glcm.iter() {
111 px_plus_y[i + j] += g_ij;
112 }
113
114 let mut sum_average = 0.0;
115 let mut sum_variance = 0.0;
116 for (i, px_plus_y_k) in px_plus_y.iter().enumerate().take(2 * glcm.rows()) {
117 let k = i as f32;
118 sum_average += k * px_plus_y_k;
119 sum_variance += k * k * px_plus_y_k;
120 }
121
122 sum_variance - sum_average * sum_average
123}
124
125#[inline]
126pub fn texture_sum_entropy(glcm: &GLCM) -> f32 {
127 let mut px_plus_y = vec![0.0; 2 * glcm.rows()];
128 for (i, j, g_ij) in glcm.iter() {
129 px_plus_y[i + j] += g_ij;
130 }
131
132 let mut sum_entropy = 0.0;
133 for px_plus_y_k in px_plus_y.iter().take(2 * glcm.rows()) {
134 let buffer = if *px_plus_y_k <= f32::EPSILON {
135 1.0
136 } else {
137 0.0
138 };
139
140 sum_entropy += px_plus_y_k * (px_plus_y_k + buffer).log2();
141 }
142
143 -sum_entropy
144}
145
146#[inline]
147pub fn texture_entropy(glcm: &GLCM) -> f32 {
148 -glcm.iter().fold(0.0, |acc, (_, _, g_ij)| {
149 acc + g_ij * (g_ij + f32::EPSILON).log2()
150 })
151}
152
153#[inline]
154pub fn texture_difference_variance(glcm: &GLCM) -> f32 {
155 let mut px_minus_y = vec![0.0; glcm.rows()];
156 for (i, j, g_ij) in glcm.iter() {
157 let index = (i as i32 - j as i32).unsigned_abs() as usize;
158 px_minus_y[index] += g_ij;
159 }
160
161 let mean_x_minus_y = px_minus_y.iter().sum::<f32>() / px_minus_y.len() as f32;
162 let variance = px_minus_y.iter().enumerate().fold(0.0, |acc, (_, &x)| {
163 acc + (x - mean_x_minus_y) * (x - mean_x_minus_y)
164 });
165
166 variance / px_minus_y.len() as f32
167}
168
169#[inline]
170pub fn texture_difference_entropy(glcm: &GLCM) -> f32 {
171 let mut px_minus_y = vec![0.0; glcm.rows()];
172 for (i, j, g_ij) in glcm.iter() {
173 let index = (i as i32 - j as i32).unsigned_abs() as usize;
174 px_minus_y[index] += g_ij;
175 }
176
177 -px_minus_y
178 .iter()
179 .fold(0.0, |acc, &x| acc + x * (x + f32::EPSILON).log2())
180}
181
182#[inline]
183pub fn texture_infocorr_1(glcm: &GLCM) -> f32 {
184 let (px, py) = glcm.margin_sums();
185
186 let hx = -px
187 .iter()
188 .fold(0.0, |acc, &x| acc + x * (x + f32::EPSILON).log2());
189 let hy = -py
190 .iter()
191 .fold(0.0, |acc, &y| acc + y * (y + f32::EPSILON).log2());
192
193 let mut hxy1 = 0.0;
194 let mut hxy2 = 0.0;
195 for (i, j, g_ij) in glcm.iter() {
196 hxy1 += g_ij * (g_ij + f32::EPSILON).log2();
197 hxy2 += px[i] * py[j] * (px[i] * py[j] + f32::EPSILON).log2();
198 }
199
200 (hxy2 - hxy1) / hx.max(hy)
201}
202
203#[inline]
204pub fn texture_infocorr_2(glcm: &GLCM) -> f32 {
205 let (px, py) = glcm.margin_sums();
206
207 let mut hxy1 = 0.0;
208 let mut hxy2 = 0.0;
209 for (i, j, g_ij) in glcm.iter() {
210 hxy1 += g_ij * (g_ij + f32::EPSILON).log2();
211 hxy2 += px[i] * py[j] * (px[i] * py[j] + f32::EPSILON).log2();
212 }
213
214 (1.0 - (-2.0 * (hxy1 - hxy2)).exp()).sqrt()
215}
216
217#[inline]
218pub fn haralick_features(glcm: &GLCM) -> [f32; 13] {
219 let (px, py) = glcm.margin_sums();
220
221 let (mut ux, mut uy) = (0.0, 0.0);
222 let (mut hx, mut hy) = (0.0, 0.0);
223
224 for i in 0..px.len() {
225 ux += (i as f32 + 1.0) * px[i];
226 uy += (i as f32 + 1.0) * py[i];
227 hx -= px[i] * (px[i] + f32::EPSILON).log2();
228 hy -= py[i] * (py[i] + f32::EPSILON).log2();
229 }
230
231 let sx = px
232 .iter()
233 .enumerate()
234 .fold(0.0, |acc, (i, &x)| {
235 acc + (i as f32 + 1.0 - ux) * (i as f32 + 1.0 - ux) * x
236 })
237 .sqrt();
238
239 let sy = py
240 .iter()
241 .enumerate()
242 .fold(0.0, |acc, (i, &y)| {
243 acc + (i as f32 + 1.0 - uy) * (i as f32 + 1.0 - uy) * y
244 })
245 .sqrt();
246
247 let mut hxy1 = 0.0;
248 let mut hxy2 = 0.0;
249
250 let mut px_plus_y = vec![0.0; 2 * glcm.rows()];
251 let mut px_minus_y = vec![0.0; glcm.rows()];
252
253 let mut energy = 0.0;
254 let mut contrast = 0.0;
255 let mut correlation = 0.0;
256 let mut sum_of_squares = 0.0;
257 let mut inverse_difference_moment = 0.0;
258 let mut sum_average = 0.0;
259 let mut sum_variance = 0.0;
260 let mut sum_entropy = 0.0;
261 let mut entropy = 0.0;
262 let mut difference_variance = 0.0;
263 let mut difference_entropy = 0.0;
264
265 for (i, j, g_ij) in glcm.iter() {
266 hxy1 += g_ij * (g_ij + f32::EPSILON).log2();
267 hxy2 += px[i] * py[j] * (px[i] * py[j] + f32::EPSILON).log2();
268
269 let index = (i as i32 - j as i32).unsigned_abs() as usize;
270 px_minus_y[index] += g_ij;
271 px_plus_y[i + j] += g_ij;
272
273 let i = i as f32;
274 let j = j as f32;
275 let d = i - j;
276 let dsq = d * d;
277
278 energy += g_ij * g_ij;
279 contrast += dsq * g_ij;
280 correlation += ((i + 1.0 - ux) * (j + 1.0 - uy) * g_ij) / (sx * sy);
281 sum_of_squares += (i + 1.0 - ux) * (i + 1.0 - ux) * g_ij;
282 inverse_difference_moment += (1.0 / (1.0 + dsq)) * g_ij;
283 entropy += g_ij * (g_ij + f32::EPSILON).log2();
284 }
285
286 for (i, px_plus_y_k) in px_plus_y.iter().enumerate().take(2 * glcm.rows()) {
287 let k = i as f32;
288 sum_average += k * px_plus_y_k;
289 sum_variance += k * k * px_plus_y_k;
290
291 let buffer = if *px_plus_y_k <= f32::EPSILON {
292 1.0
293 } else {
294 0.0
295 };
296 sum_entropy += px_plus_y_k * (px_plus_y_k + buffer).log2();
297 }
298
299 sum_variance -= sum_average * sum_average;
300
301 let u_x_minus_y = px_minus_y.iter().sum::<f32>() / px_minus_y.len() as f32;
302
303 for px_minus_y_k in px_minus_y.iter().take(glcm.rows()) {
304 difference_variance += (px_minus_y_k - u_x_minus_y) * (px_minus_y_k - u_x_minus_y);
305 difference_entropy += px_minus_y_k * (px_minus_y_k + f32::EPSILON).log2();
306 }
307
308 difference_variance /= px_minus_y.len() as f32;
309
310 let information_measure_of_correlation_1 = (hxy2 - hxy1) / hx.max(hy);
311 let information_measure_of_correlation_2 = (1.0 - (-2.0 * (hxy1 - hxy2)).exp()).sqrt();
312
313 [
314 energy,
315 contrast,
316 correlation,
317 sum_of_squares,
318 inverse_difference_moment,
319 sum_average,
320 sum_variance,
321 -sum_entropy,
322 -entropy,
323 difference_variance,
324 -difference_entropy,
325 information_measure_of_correlation_1,
326 information_measure_of_correlation_2,
327 ]
328}
329
330#[inline]
331pub fn descriptors<T>(pixels: &[T], width: usize, height: usize, channels: usize) -> [f32; 13]
332where
333 T: ToPrimitive,
334{
335 let mut haralick: [f32; 13] = [0.0; 13];
336 for i in [0, 45, 90, 135].iter() {
337 for glcm in glcm_multichannel(pixels, width, height, channels, *i as f32, 1.0).iter() {
338 let features = haralick_features(glcm);
339 for j in 0..13 {
340 haralick[j] += features[j] / (4.0 * channels as f32);
341 }
342 }
343 }
344
345 haralick
346}
347
348#[inline]
349pub fn objects<T, Container>(object: &PineappleViewBuffer<T, Container>) -> [f32; 13]
350where
351 T: ToPrimitive + FromPrimitive,
352 Container: Deref<Target = [T]>,
353{
354 let mut haralick: [f32; 13] = [0.0; 13];
355 for i in [0, 45, 90, 135].iter() {
356 for glcm in glcm_multichannel_object(object, *i as f32, 1.0).iter() {
357 let features = haralick_features(glcm);
358 for j in 0..13 {
359 haralick[j] += features[j] / (4.0 * object.channels() as f32);
360 }
361 }
362 }
363
364 haralick
365}
366
367#[cfg(test)]
368mod test {
369
370 use super::*;
371
372 use crate::im::PineappleBuffer;
373
374 fn square_image() -> [u8; 4] {
375 [0, 0, 255, 255]
376 }
377
378 const EPS: f32 = 1e-6;
379
380 #[test]
381 fn test_texture_energy() {
382 let energy = texture_energy(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
383 assert_eq!(energy, 0.5);
384 }
385
386 #[test]
387 fn test_texture_contrast() {
388 let contrast = texture_contrast(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
389 assert_eq!(contrast, 0.0);
390 }
391
392 #[test]
393 fn test_texture_correlation() {
394 let correlation = texture_correlation(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
395 assert_eq!(correlation, 1.0);
396 }
397
398 #[test]
399 fn test_texture_sum_of_squares() {
400 let sum_of_squares =
401 texture_sum_of_squares(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
402 assert_eq!(sum_of_squares, 992.25);
403 }
404
405 #[test]
406 fn test_texture_inverse_difference() {
407 let inverse_difference =
408 texture_inverse_difference(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
409
410 assert_eq!(inverse_difference, 1.0);
411 }
412
413 #[test]
414 fn test_texture_sum_average() {
415 let sum_average = texture_sum_average(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
416 assert_eq!(sum_average, 63.0);
417 }
418
419 #[test]
420 fn test_texture_sum_variance() {
421 let sum_variance = texture_sum_variance(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
422 assert_eq!(sum_variance, 3969.0);
423 }
424
425 #[test]
426 fn test_texture_sum_entropy() {
427 let sum_entropy = texture_sum_entropy(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
428 assert_eq!(sum_entropy, 1.0);
429 }
430
431 #[test]
432 fn test_texture_entropy() {
433 let entropy = texture_entropy(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
434 assert!((entropy - 1.0).abs() < EPS);
435 }
436
437 #[test]
438 fn test_texture_difference_variance() {
439 let difference_variance =
440 texture_difference_variance(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
441
442 assert!((difference_variance - 0.015380859).abs() < EPS);
443 }
444
445 #[test]
446 fn test_texture_difference_entropy() {
447 let difference_entropy =
448 texture_difference_entropy(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
449
450 assert!((difference_entropy - 0.0).abs() < EPS);
451 }
452
453 #[test]
454 fn test_texture_information_measure_of_correlation_1() {
455 let imc1 = texture_infocorr_1(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
456 assert_eq!(imc1, -1.0);
457 }
458
459 #[test]
460 fn test_texture_information_measure_of_correlation_2() {
461 let imc2 = texture_infocorr_2(&GLCM::new(&square_image(), 2, 2, 0, 1, 0.0, 1.0));
462 assert!((imc2 - 0.92987347).abs() < EPS);
463 }
464
465 #[test]
466 fn test_haralick_features() {
467 let pixels = square_image();
468 let comatrix = GLCM::new(&pixels, 2, 2, 0, 1, 0.0, 1.0);
469 let features = haralick_features(&comatrix);
470
471 assert_eq!(features.len(), 13);
472
473 let energy = texture_energy(&comatrix);
474 let contrast = texture_contrast(&comatrix);
475 let correlation = texture_correlation(&comatrix);
476 let sum_of_squares = texture_sum_of_squares(&comatrix);
477 let inverse_difference = texture_inverse_difference(&comatrix);
478 let sum_average = texture_sum_average(&comatrix);
479 let sum_variance = texture_sum_variance(&comatrix);
480 let sum_entropy = texture_sum_entropy(&comatrix);
481 let entropy = texture_entropy(&comatrix);
482 let difference_variance = texture_difference_variance(&comatrix);
483 let difference_entropy = texture_difference_entropy(&comatrix);
484 let imc1 = texture_infocorr_1(&comatrix);
485 let imc2 = texture_infocorr_2(&comatrix);
486
487 assert_eq!(features[0], energy);
488 assert_eq!(features[1], contrast);
489 assert_eq!(features[2], correlation);
490 assert_eq!(features[3], sum_of_squares);
491 assert_eq!(features[4], inverse_difference);
492 assert_eq!(features[5], sum_average);
493 assert_eq!(features[6], sum_variance);
494 assert_eq!(features[7], sum_entropy);
495 assert_eq!(features[8], entropy);
496 assert_eq!(features[9], difference_variance);
497 assert_eq!(features[10], difference_entropy);
498 assert_eq!(features[11], imc1);
499 assert_eq!(features[12], imc2);
500 }
501
502 #[test]
503 fn test_object_texture() {
504 let pixels = square_image();
505 let buffer = PineappleBuffer::new(2, 2, 1, pixels.to_vec()).unwrap();
506 let object = PineappleViewBuffer::new(0, 0, 2, 2, &buffer);
507
508 let texture_array = descriptors(&pixels, 2, 2, 1);
509 let texture_object = objects(&object);
510
511 assert_eq!(texture_array, texture_object);
512 }
513}