1use crate::error::{IrisError, Result};
2use crate::image::Image;
3use burn::tensor::{Tensor, TensorData, backend::Backend};
4
5impl<B: Backend> Image<B> {
6 pub fn calc_hist(&self) -> Result<Vec<Vec<u32>>> {
9 let dims = self.tensor.dims();
10 let c = dims[0];
11 let h = dims[1];
12 let w = dims[2];
13
14 let tensor_data = self.tensor.clone().into_data();
15 let flat_vals: Vec<f32> = tensor_data.iter::<f32>().collect();
16 let mut histograms = vec![vec![0u32; 256]; c];
17
18 for ch in 0..c {
19 for y in 0..h {
20 for x in 0..w {
21 let val = flat_vals[ch * h * w + y * w + x];
22 let bin = (val.clamp(0.0, 1.0) * 255.0) as usize;
23 histograms[ch][bin] += 1;
24 }
25 }
26 }
27
28 Ok(histograms)
29 }
30
31 pub fn equalize_hist(&self) -> Result<Self> {
33 let gray = self.grayscale()?;
34 let dims = gray.tensor.dims();
35 let h = dims[1];
36 let w = dims[2];
37
38 let device = gray.tensor.device();
39 let tensor_data = gray.tensor.clone().into_data();
40 let flat_vals: Vec<f32> = tensor_data.iter::<f32>().collect();
41 let mut out_vals = vec![0.0f32; h * w];
42
43 let mut hist = [0u32; 256];
45 for &val in &flat_vals {
46 let bin = (val.clamp(0.0, 1.0) * 255.0) as usize;
47 hist[bin] += 1;
48 }
49
50 let mut cdf = [0u32; 256];
52 let mut sum = 0u32;
53 for i in 0..256 {
54 sum += hist[i];
55 cdf[i] = sum;
56 }
57
58 let cdf_min = cdf.iter().find(|&&x| x > 0).copied().unwrap_or(0) as f32;
60 let total = (h * w) as f32;
61
62 let mut lut = [0.0f32; 256];
64 if total > cdf_min {
65 for i in 0..256 {
66 lut[i] = ((cdf[i] as f32 - cdf_min) / (total - cdf_min) * 255.0).round() / 255.0;
67 }
68 }
69
70 for i in 0..(h * w) {
71 let bin = (flat_vals[i].clamp(0.0, 1.0) * 255.0) as usize;
72 out_vals[i] = lut[bin];
73 }
74
75 let new_data = TensorData::new(out_vals, [1, h, w]);
76 let new_tensor = Tensor::<B, 3>::from_data(new_data, &device);
77 Ok(Image::new(new_tensor))
78 }
79
80 pub fn equalize_hist_color(&self) -> Result<Self> {
83 let dims = self.tensor.dims();
84 if dims[0] != 3 {
85 return Err(IrisError::InvalidParameter(
86 "Input must be a 3-channel RGB image".into(),
87 ));
88 }
89
90 let ycrcb = self.rgb_to_ycrcb()?;
92
93 let y_channel = ycrcb.tensor.clone().slice([0..1, 0..dims[1], 0..dims[2]]);
95 let y_img = Image::new(y_channel);
96
97 let y_equalized = y_img.equalize_hist()?;
99
100 let cr = ycrcb.tensor.clone().slice([1..2, 0..dims[1], 0..dims[2]]);
102 let cb = ycrcb.tensor.clone().slice([2..3, 0..dims[1], 0..dims[2]]);
103 let ycrcb_equalized =
104 Image::merge_channels(&[y_equalized, Image::new(cr), Image::new(cb)])?;
105
106 ycrcb_equalized.ycrcb_to_rgb()
108 }
109
110 pub fn clahe(&self, clip_limit: f32, grid_size: usize) -> Result<Self> {
114 if grid_size == 0 {
115 return Err(IrisError::InvalidParameter("grid_size must be > 0".into()));
116 }
117
118 let gray = self.grayscale()?;
119 let dims = gray.tensor.dims();
120 let h = dims[1];
121 let w = dims[2];
122
123 let device = gray.tensor.device();
124 let tensor_data = gray.tensor.clone().into_data();
125 let flat_vals: Vec<f32> = tensor_data.iter::<f32>().collect();
126 let mut out_vals = flat_vals.clone();
127
128 let tile_h = h / grid_size;
129 let tile_w = w / grid_size;
130
131 if tile_h == 0 || tile_w == 0 {
132 return Err(IrisError::InvalidParameter(
133 "Image too small for given grid_size".into(),
134 ));
135 }
136
137 for ty in 0..grid_size {
138 for tx in 0..grid_size {
139 let y0 = ty * tile_h;
140 let x0 = tx * tile_w;
141 let y1 = if ty == grid_size - 1 { h } else { y0 + tile_h };
142 let x1 = if tx == grid_size - 1 { w } else { x0 + tile_w };
143
144 let tile_pixels = (y1 - y0) * (x1 - x0);
145
146 let mut hist = [0u32; 256];
148 for y in y0..y1 {
149 for x in x0..x1 {
150 let bin = (flat_vals[y * w + x].clamp(0.0, 1.0) * 255.0) as usize;
151 hist[bin] += 1;
152 }
153 }
154
155 if clip_limit > 0.0 {
157 let limit = (clip_limit * tile_pixels as f32 / 256.0) as u32;
158 let mut excess = 0u32;
159 for bin in 0..256 {
160 if hist[bin] > limit {
161 excess += hist[bin] - limit;
162 hist[bin] = limit;
163 }
164 }
165 let avg_inc = excess / 256;
167 let rem = excess % 256;
168 for bin in 0..256 {
169 hist[bin] += avg_inc;
170 if bin < rem as usize {
171 hist[bin] += 1;
172 }
173 }
174 }
175
176 let mut cdf = [0u32; 256];
178 let mut sum = 0u32;
179 for i in 0..256 {
180 sum += hist[i];
181 cdf[i] = sum;
182 }
183
184 let cdf_min = cdf.iter().find(|&&x| x > 0).copied().unwrap_or(0) as f32;
185 let total = tile_pixels as f32;
186
187 let mut lut = [0.0f32; 256];
188 if total > cdf_min {
189 for i in 0..256 {
190 lut[i] =
191 ((cdf[i] as f32 - cdf_min) / (total - cdf_min) * 255.0).round() / 255.0;
192 }
193 }
194
195 for y in y0..y1 {
197 for x in x0..x1 {
198 let bin = (flat_vals[y * w + x].clamp(0.0, 1.0) * 255.0) as usize;
199 out_vals[y * w + x] = lut[bin];
200 }
201 }
202 }
203 }
204
205 let new_data = TensorData::new(out_vals, [1, h, w]);
206 let new_tensor = Tensor::<B, 3>::from_data(new_data, &device);
207 Ok(Image::new(new_tensor))
208 }
209
210 pub fn apply_lut(&self, lut: &[f32; 256]) -> Result<Self> {
213 let dims = self.tensor.dims();
214 let c = dims[0];
215 let h = dims[1];
216 let w = dims[2];
217
218 let device = self.tensor.device();
219 let tensor_data = self.tensor.clone().into_data();
220 let flat_vals: Vec<f32> = tensor_data.iter::<f32>().collect();
221 let mut out_vals = vec![0.0f32; c * h * w];
222
223 for i in 0..(c * h * w) {
224 let bin = (flat_vals[i].clamp(0.0, 1.0) * 255.0) as usize;
225 out_vals[i] = lut[bin];
226 }
227
228 let new_data = TensorData::new(out_vals, [c, h, w]);
229 let new_tensor = Tensor::<B, 3>::from_data(new_data, &device);
230 Ok(Image::new(new_tensor))
231 }
232
233 pub fn compare_hist(hist_a: &[f32], hist_b: &[f32], method: &str) -> Result<f64> {
236 if hist_a.len() != hist_b.len() {
237 return Err(IrisError::DimensionMismatch {
238 expected: vec![hist_a.len()],
239 actual: vec![hist_b.len()],
240 });
241 }
242
243 match method {
244 "correlation" => {
245 let n = hist_a.len() as f64;
246 let mean_a: f64 = hist_a.iter().map(|&x| x as f64).sum::<f64>() / n;
247 let mean_b: f64 = hist_b.iter().map(|&x| x as f64).sum::<f64>() / n;
248
249 let mut num = 0.0;
250 let mut den_a = 0.0;
251 let mut den_b = 0.0;
252 for i in 0..hist_a.len() {
253 let da = hist_a[i] as f64 - mean_a;
254 let db = hist_b[i] as f64 - mean_b;
255 num += da * db;
256 den_a += da * da;
257 den_b += db * db;
258 }
259 let den = (den_a * den_b).sqrt();
260 Ok(if den.abs() < 1e-10 { 0.0 } else { num / den })
261 }
262 "chi_square" => {
263 let mut sum = 0.0;
264 for i in 0..hist_a.len() {
265 let a = hist_a[i] as f64;
266 let b = hist_b[i] as f64;
267 if a + b > 0.0 {
268 sum += (a - b).powi(2) / (a + b);
269 }
270 }
271 Ok(sum)
272 }
273 "intersection" => {
274 let sum: f64 = hist_a
275 .iter()
276 .zip(hist_b.iter())
277 .map(|(&a, &b)| (a as f64).min(b as f64))
278 .sum();
279 Ok(sum)
280 }
281 "hellinger" => {
282 let mut sum = 0.0;
283 for i in 0..hist_a.len() {
284 let a = (hist_a[i] as f64).sqrt();
285 let b = (hist_b[i] as f64).sqrt();
286 sum += (a - b).powi(2);
287 }
288 Ok((sum / 2.0).sqrt())
289 }
290 _ => Err(IrisError::InvalidParameter(format!(
291 "Unknown comparison method: {method}. Use correlation, chi_square, intersection, or hellinger"
292 ))),
293 }
294 }
295
296 pub fn compare_hist_color(
299 hist_a: &[Vec<f32>],
300 hist_b: &[Vec<f32>],
301 method: &str,
302 ) -> Result<Vec<f64>> {
303 if hist_a.len() != hist_b.len() {
304 return Err(IrisError::DimensionMismatch {
305 expected: vec![hist_a.len()],
306 actual: vec![hist_b.len()],
307 });
308 }
309
310 let mut results = Vec::with_capacity(hist_a.len());
311 for (a, b) in hist_a.iter().zip(hist_b.iter()) {
312 let score = Self::compare_hist(a, b, method)?;
313 results.push(score);
314 }
315 Ok(results)
316 }
317
318 pub fn calc_hist_2d(
324 &self,
325 channel_x: usize,
326 channel_y: usize,
327 bins: usize,
328 ) -> Result<Tensor<B, 2>> {
329 if bins == 0 {
330 return Err(IrisError::InvalidParameter("bins must be > 0".into()));
331 }
332
333 let dims = self.tensor.dims();
334 let c = dims[0];
335 let h = dims[1];
336 let w = dims[2];
337
338 if channel_x >= c || channel_y >= c {
339 return Err(IrisError::DimensionMismatch {
340 expected: vec![c],
341 actual: vec![channel_x.max(channel_y) + 1],
342 });
343 }
344
345 let tensor_data = self.tensor.clone().into_data();
346 let flat_vals: Vec<f32> = tensor_data.iter::<f32>().collect();
347
348 let mut hist = vec![0u32; bins * bins];
349
350 for y in 0..h {
351 for x in 0..w {
352 let val_x = flat_vals[channel_x * h * w + y * w + x];
353 let val_y = flat_vals[channel_y * h * w + y * w + x];
354
355 let bin_x =
356 ((val_x.clamp(0.0, 1.0) * (bins as f32 - 1.0)).round() as usize).min(bins - 1);
357 let bin_y =
358 ((val_y.clamp(0.0, 1.0) * (bins as f32 - 1.0)).round() as usize).min(bins - 1);
359
360 hist[bin_y * bins + bin_x] += 1;
361 }
362 }
363
364 let max_val = hist.iter().copied().max().unwrap_or(1) as f32;
366 let hist_f32: Vec<f32> = hist.iter().map(|&v| v as f32 / max_val).collect();
367
368 let device = self.tensor.device();
369 let new_data = TensorData::new(hist_f32, [bins, bins]);
370 let new_tensor = Tensor::<B, 2>::from_data(new_data, &device);
371 Ok(new_tensor)
372 }
373
374 pub fn equalize_hist_adaptive(&self, clip_limit: f32, grid_size: usize) -> Result<Self> {
380 if grid_size == 0 {
381 return Err(IrisError::InvalidParameter("grid_size must be > 0".into()));
382 }
383
384 let gray = self.grayscale()?;
385 let dims = gray.tensor.dims();
386 let h = dims[1];
387 let w = dims[2];
388
389 let device = gray.tensor.device();
390 let tensor_data = gray.tensor.clone().into_data();
391 let flat_vals: Vec<f32> = tensor_data.iter::<f32>().collect();
392
393 let tile_h = h / grid_size;
394 let tile_w = w / grid_size;
395
396 if tile_h == 0 || tile_w == 0 {
397 return Err(IrisError::InvalidParameter(
398 "Image too small for given grid_size".into(),
399 ));
400 }
401
402 let mut tile_luts: Vec<Vec<f32>> = Vec::with_capacity(grid_size * grid_size);
404
405 for ty in 0..grid_size {
406 for tx in 0..grid_size {
407 let y0 = ty * tile_h;
408 let x0 = tx * tile_w;
409 let y1 = if ty == grid_size - 1 { h } else { y0 + tile_h };
410 let x1 = if tx == grid_size - 1 { w } else { x0 + tile_w };
411
412 let tile_pixels = (y1 - y0) * (x1 - x0);
413
414 let mut hist = [0u32; 256];
416 for y in y0..y1 {
417 for x in x0..x1 {
418 let bin = (flat_vals[y * w + x].clamp(0.0, 1.0) * 255.0) as usize;
419 hist[bin] += 1;
420 }
421 }
422
423 if clip_limit > 0.0 {
425 let limit = (clip_limit * tile_pixels as f32 / 256.0) as u32;
426 let mut excess = 0u32;
427 for bin in 0..256 {
428 if hist[bin] > limit {
429 excess += hist[bin] - limit;
430 hist[bin] = limit;
431 }
432 }
433 let avg_inc = excess / 256;
434 let rem = excess % 256;
435 for bin in 0..256 {
436 hist[bin] += avg_inc;
437 if bin < rem as usize {
438 hist[bin] += 1;
439 }
440 }
441 }
442
443 let mut cdf = [0u32; 256];
445 let mut sum = 0u32;
446 for i in 0..256 {
447 sum += hist[i];
448 cdf[i] = sum;
449 }
450
451 let cdf_min = cdf.iter().find(|&&x| x > 0).copied().unwrap_or(0) as f32;
452 let total = tile_pixels as f32;
453
454 let mut lut = [0.0f32; 256];
455 if total > cdf_min {
456 for i in 0..256 {
457 lut[i] =
458 ((cdf[i] as f32 - cdf_min) / (total - cdf_min) * 255.0).round() / 255.0;
459 }
460 }
461
462 tile_luts.push(lut.to_vec());
463 }
464 }
465
466 let mut out_vals = vec![0.0f32; h * w];
468
469 for y in 0..h {
470 for x in 0..w {
471 let tx = (x as f32 / tile_w as f32 - 0.5).clamp(0.0, (grid_size - 1) as f32);
473 let ty = (y as f32 / tile_h as f32 - 0.5).clamp(0.0, (grid_size - 1) as f32);
474
475 let tx0 = tx.floor() as usize;
476 let ty0 = ty.floor() as usize;
477 let tx1 = (tx0 + 1).min(grid_size - 1);
478 let ty1 = (ty0 + 1).min(grid_size - 1);
479
480 let fx = tx - tx0 as f32;
481 let fy = ty - ty0 as f32;
482
483 let bin = (flat_vals[y * w + x].clamp(0.0, 1.0) * 255.0) as usize;
484
485 let lut00 = &tile_luts[ty0 * grid_size + tx0];
486 let lut10 = &tile_luts[ty0 * grid_size + tx1];
487 let lut01 = &tile_luts[ty1 * grid_size + tx0];
488 let lut11 = &tile_luts[ty1 * grid_size + tx1];
489
490 let v00 = lut00[bin];
491 let v10 = lut10[bin];
492 let v01 = lut01[bin];
493 let v11 = lut11[bin];
494
495 let val = v00 * (1.0 - fx) * (1.0 - fy)
497 + v10 * fx * (1.0 - fy)
498 + v01 * (1.0 - fx) * fy
499 + v11 * fx * fy;
500
501 out_vals[y * w + x] = val;
502 }
503 }
504
505 let new_data = TensorData::new(out_vals, [1, h, w]);
506 let new_tensor = Tensor::<B, 3>::from_data(new_data, &device);
507 Ok(Image::new(new_tensor))
508 }
509}
510
511#[cfg(test)]
512mod tests {
513 use super::*;
514 use crate::test_helpers::{TestBackend, test_device};
515
516 #[test]
517 fn test_histogram_operations() {
518 let device = test_device();
519 let flat_data = vec![0.1f32, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8];
520 let tensor =
521 Tensor::<TestBackend, 3>::from_data(TensorData::new(flat_data, [1, 2, 4]), &device);
522 let img = Image::new(tensor);
523
524 let hists = img.calc_hist().unwrap();
525 assert_eq!(hists.len(), 1);
526 assert_eq!(hists[0].len(), 256);
527
528 let eq = img.equalize_hist().unwrap();
529 assert_eq!(eq.shape(), [1, 2, 4]);
530 }
531
532 #[test]
533 fn test_equalize_hist_color() {
534 let device = test_device();
535 let flat_data = vec![
536 0.2f32, 0.4, 0.6, 0.8, 0.1, 0.3, 0.5, 0.7, 0.9, 0.0, 0.2, 0.4,
537 ];
538 let tensor =
539 Tensor::<TestBackend, 3>::from_data(TensorData::new(flat_data, [3, 2, 2]), &device);
540 let img = Image::new(tensor);
541
542 let eq = img.equalize_hist_color().unwrap();
543 assert_eq!(eq.shape(), [3, 2, 2]);
544 }
545
546 #[test]
547 fn test_compare_hist_color() {
548 let hist_a = vec![
549 vec![1.0, 2.0, 3.0, 4.0],
550 vec![2.0, 3.0, 4.0, 5.0],
551 vec![3.0, 4.0, 5.0, 6.0],
552 ];
553 let hist_b = vec![
554 vec![1.0, 2.0, 3.0, 4.0],
555 vec![2.0, 3.0, 4.0, 5.0],
556 vec![3.0, 4.0, 5.0, 6.0],
557 ];
558
559 let results =
560 Image::<TestBackend>::compare_hist_color(&hist_a, &hist_b, "correlation").unwrap();
561 assert_eq!(results.len(), 3);
562 for r in results {
563 assert!((r - 1.0).abs() < 1e-5);
564 }
565
566 let chi_results =
567 Image::<TestBackend>::compare_hist_color(&hist_a, &hist_b, "chi_square").unwrap();
568 for r in chi_results {
569 assert!(r.abs() < 1e-5);
570 }
571 }
572
573 #[test]
574 fn test_clahe() {
575 let device = test_device();
576 let data: Vec<f32> = (0..64).map(|i| (i as f32) / 64.0).collect();
577 let tensor = Tensor::<TestBackend, 3>::from_data(TensorData::new(data, [1, 8, 8]), &device);
578 let img = Image::new(tensor);
579 let result = img.clahe(2.0, 4).unwrap();
580 assert_eq!(result.shape(), [1, 8, 8]);
581 }
582
583 #[test]
584 fn test_apply_lut() {
585 let device = test_device();
586 let data = vec![0.0f32, 0.5, 1.0];
587 let tensor = Tensor::<TestBackend, 3>::from_data(TensorData::new(data, [1, 1, 3]), &device);
588 let img = Image::new(tensor);
589
590 let mut lut = [0.0f32; 256];
591 for i in 0..256 {
592 lut[i] = 1.0 - (i as f32) / 255.0; }
594 let result = img.apply_lut(&lut).unwrap();
595 assert_eq!(result.shape(), [1, 1, 3]);
596 let vals: Vec<f32> = result.tensor.into_data().iter::<f32>().collect();
597 assert!((vals[0] - 1.0).abs() < 1e-5); assert!((vals[2] - 0.0).abs() < 1e-5); }
600
601 #[test]
602 fn test_compare_hist() {
603 let hist_a = vec![1.0, 2.0, 3.0, 4.0];
604 let hist_b = vec![1.0, 2.0, 3.0, 4.0];
605 let corr = Image::<TestBackend>::compare_hist(&hist_a, &hist_b, "correlation").unwrap();
606 assert!((corr - 1.0).abs() < 1e-5); let chi = Image::<TestBackend>::compare_hist(&hist_a, &hist_b, "chi_square").unwrap();
609 assert!((chi).abs() < 1e-5); let inter = Image::<TestBackend>::compare_hist(&hist_a, &hist_b, "intersection").unwrap();
612 assert!((inter - 10.0).abs() < 1e-5); let hel = Image::<TestBackend>::compare_hist(&hist_a, &hist_b, "hellinger").unwrap();
615 assert!((hel).abs() < 1e-5);
616 }
617
618 #[test]
619 fn test_compare_hist_invalid() {
620 let hist_a = vec![1.0, 2.0];
621 let hist_b = vec![1.0, 2.0];
622 assert!(Image::<TestBackend>::compare_hist(&hist_a, &hist_b, "invalid").is_err());
623 }
624
625 #[test]
626 fn test_calc_hist_2d() {
627 let device = test_device();
628 let mut flat_data = Vec::new();
630 for y in 0..4 {
632 for x in 0..4 {
633 flat_data.push((y * 4 + x) as f32 / 15.0);
634 }
635 }
636 for y in 0..4 {
638 for x in 0..4 {
639 flat_data.push(1.0 - (y * 4 + x) as f32 / 15.0);
640 }
641 }
642 let tensor =
643 Tensor::<TestBackend, 3>::from_data(TensorData::new(flat_data, [2, 4, 4]), &device);
644 let img = Image::new(tensor);
645
646 let hist_2d = img.calc_hist_2d(0, 1, 4).unwrap();
647 let dims = hist_2d.dims();
648 assert_eq!(dims, [4, 4]);
649 let vals: Vec<f32> = hist_2d.into_data().iter::<f32>().collect();
650 assert!(vals.iter().all(|&v| v >= 0.0));
652 assert!(vals.iter().any(|&v| v > 0.0));
653 }
654
655 #[test]
656 fn test_equalize_hist_adaptive() {
657 let device = test_device();
658 let data: Vec<f32> = (0..64).map(|i| (i as f32) / 64.0).collect();
659 let tensor = Tensor::<TestBackend, 3>::from_data(TensorData::new(data, [1, 8, 8]), &device);
660 let img = Image::new(tensor);
661
662 let result = img.equalize_hist_adaptive(2.0, 2).unwrap();
663 assert_eq!(result.shape(), [1, 8, 8]);
664 let vals: Vec<f32> = result.tensor.into_data().iter::<f32>().collect();
665 assert!(vals.iter().all(|&v| (0.0..=1.0).contains(&v)));
667 }
668}