Skip to main content

iris/histogram/
mod.rs

1use crate::error::{IrisError, Result};
2use crate::image::Image;
3use burn::tensor::{Tensor, TensorData, backend::Backend};
4
5impl<B: Backend> Image<B> {
6    /// Computes the 256-bin histogram for each channel.
7    /// Returns a vector of vectors (one per channel), each containing 256 counts.
8    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    /// Performs histogram equalization on a grayscale image to enhance contrast.
32    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        // 1. Compute histogram
44        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        // 2. Compute cumulative distribution function (CDF)
51        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        // 3. Find CDF minimum non-zero value
59        let cdf_min = cdf.iter().find(|&&x| x > 0).copied().unwrap_or(0) as f32;
60        let total = (h * w) as f32;
61
62        // 4. Equalize mapping
63        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    /// Performs histogram equalization on a color image by converting to YCrCb,
81    /// equalizing the Y channel, and converting back.
82    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        // Convert RGB to YCrCb
91        let ycrcb = self.rgb_to_ycrcb()?;
92
93        // Extract Y channel as single-channel image
94        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        // Equalize histogram on Y channel
98        let y_equalized = y_img.equalize_hist()?;
99
100        // Reconstruct YCrCb with equalized Y
101        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        // Convert back to RGB
107        ycrcb_equalized.ycrcb_to_rgb()
108    }
109
110    /// Contrast Limited Adaptive Histogram Equalization (CLAHE).
111    /// Divides the image into `grid_size x grid_size` tiles and applies
112    /// histogram equalization to each tile independently with clip limit.
113    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                // Compute histogram for this tile
147                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                // Apply clip limit: redistribute excess bins
156                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                    // Redistribute excess evenly
166                    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                // Compute CDF
177                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                // Apply LUT to tile
196                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    /// Applies a lookup table (LUT) to map pixel values.
211    /// `lut` must have 256 entries mapping each possible pixel value (0..255) to a new value.
212    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    /// Compares two histograms using the specified method.
234    /// Both histograms must have the same length.
235    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    /// Compares two color histograms using the specified method.
297    /// Returns per-channel comparison results.
298    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    /// Computes a 2D histogram over two channels of a multi-channel image.
319    ///
320    /// `channel_x` and `channel_y` select which channels to histogram (0-indexed).
321    /// The result is a 2D tensor of shape `[bins, bins]` with counts normalized
322    /// so that the maximum bin value equals 1.0.
323    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        // Normalize to [0, 1] range
365        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    /// Adaptive histogram equalization (non-CLAHE version).
375    ///
376    /// Divides the image into `grid_size x grid_size` tiles, equalizes each tile,
377    /// and blends neighboring tiles using bilinear interpolation to avoid artifacts
378    /// at tile boundaries.
379    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        // Compute LUT for each tile
403        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                // Compute histogram for this tile
415                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                // Apply clip limit
424                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                // Compute CDF and LUT
444                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        // Apply LUTs with bilinear interpolation at tile boundaries
467        let mut out_vals = vec![0.0f32; h * w];
468
469        for y in 0..h {
470            for x in 0..w {
471                // Determine which tile center this pixel is nearest to
472                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                // Bilinear interpolation
496                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; // Invert
593        }
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); // 0.0 -> 1.0
598        assert!((vals[2] - 0.0).abs() < 1e-5); // 1.0 -> 0.0
599    }
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); // Identical histograms => correlation = 1
607
608        let chi = Image::<TestBackend>::compare_hist(&hist_a, &hist_b, "chi_square").unwrap();
609        assert!((chi).abs() < 1e-5); // Identical => chi_square = 0
610
611        let inter = Image::<TestBackend>::compare_hist(&hist_a, &hist_b, "intersection").unwrap();
612        assert!((inter - 10.0).abs() < 1e-5); // sum of min(a, b) = 1+2+3+4 = 10
613
614        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        // Create a 2-channel image of shape [2, 4, 4]
629        let mut flat_data = Vec::new();
630        // Channel 0: values from 0 to 1
631        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        // Channel 1: values from 1 to 0 (reversed)
637        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        // All bins should be >= 0 and at least one should be > 0
651        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        // All output values should be in [0, 1]
666        assert!(vals.iter().all(|&v| (0.0..=1.0).contains(&v)));
667    }
668}