spatialrust_vision/
analysis.rs

1//! Thresholding, histograms, contrast enhancement, and integral images.
2
3use spatialrust_image::{Image, ImageView};
4
5use crate::border::fetch;
6use crate::{BorderMode, PixelComponent, VisionError, VisionResult};
7
8/// Point-wise threshold transformation.
9#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
10pub enum ThresholdType {
11    /// `max_value` above the threshold, zero otherwise.
12    Binary,
13    /// Zero above the threshold, `max_value` otherwise.
14    BinaryInv,
15    /// Clamp values above the threshold to the threshold.
16    Truncate,
17    /// Preserve values above the threshold, zero otherwise.
18    ToZero,
19    /// Zero values above the threshold, preserve the others.
20    ToZeroInv,
21}
22
23/// Adaptive neighborhood statistic.
24#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
25pub enum AdaptiveThresholdMethod {
26    /// Uniform neighborhood mean.
27    Mean,
28    /// Gaussian-weighted neighborhood mean.
29    Gaussian,
30}
31
32/// Applies a fixed threshold independently to every component.
33pub fn threshold<T: PixelComponent, const CHANNELS: usize>(
34    input: ImageView<'_, T, CHANNELS>,
35    threshold: f64,
36    max_value: f64,
37    threshold_type: ThresholdType,
38) -> VisionResult<Image<T, CHANNELS>> {
39    if !threshold.is_finite() || !max_value.is_finite() {
40        return Err(VisionError::InvalidParameter("threshold and max_value must be finite".into()));
41    }
42    let mut output = Vec::with_capacity(input.width() * input.height() * CHANNELS);
43    for y in 0..input.height() {
44        for x in 0..input.width() {
45            let pixel = input.get(x, y).expect("input coordinate in bounds");
46            output.extend(std::array::from_fn::<_, CHANNELS, _>(|channel| {
47                T::from_f64(apply_threshold(
48                    pixel[channel].to_f64(),
49                    threshold,
50                    max_value,
51                    threshold_type,
52                ))
53            }));
54        }
55    }
56    Ok(Image::try_new_with_metadata(input.width(), input.height(), output, input.metadata())?)
57}
58
59/// Selects an Otsu threshold for an 8-bit grayscale image and applies it.
60pub fn otsu_threshold_u8(
61    input: ImageView<'_, u8, 1>,
62    max_value: u8,
63    threshold_type: ThresholdType,
64) -> VisionResult<(u8, Image<u8, 1>)> {
65    let histogram = histogram_u8(input);
66    let selected = otsu_from_histogram(&histogram, input.width() * input.height()) as u8;
67    let output = threshold(input, f64::from(selected), f64::from(max_value), threshold_type)?;
68    Ok((selected, output))
69}
70
71/// Selects an Otsu threshold for a 16-bit grayscale image and applies it.
72pub fn otsu_threshold_u16(
73    input: ImageView<'_, u16, 1>,
74    max_value: u16,
75    threshold_type: ThresholdType,
76) -> VisionResult<(u16, Image<u16, 1>)> {
77    let mut histogram = vec![0_u64; 65_536];
78    for y in 0..input.height() {
79        for x in 0..input.width() {
80            histogram[input.get(x, y).expect("coordinate in bounds")[0] as usize] += 1;
81        }
82    }
83    let selected = otsu_from_histogram(&histogram, input.width() * input.height()) as u16;
84    let output = threshold(input, f64::from(selected), f64::from(max_value), threshold_type)?;
85    Ok((selected, output))
86}
87
88/// Applies adaptive binary thresholding to an 8-bit grayscale image.
89pub fn adaptive_threshold(
90    input: ImageView<'_, u8, 1>,
91    max_value: u8,
92    method: AdaptiveThresholdMethod,
93    threshold_type: ThresholdType,
94    block_size: usize,
95    c: f64,
96    border: BorderMode<u8, 1>,
97) -> VisionResult<Image<u8, 1>> {
98    if !matches!(threshold_type, ThresholdType::Binary | ThresholdType::BinaryInv) {
99        return Err(VisionError::InvalidParameter(
100            "adaptive threshold supports only binary and binary-inverse output".into(),
101        ));
102    }
103    if block_size <= 1 || block_size % 2 == 0 {
104        return Err(VisionError::InvalidParameter(
105            "adaptive block size must be odd and greater than one".into(),
106        ));
107    }
108    if !c.is_finite() {
109        return Err(VisionError::InvalidParameter("adaptive constant must be finite".into()));
110    }
111    let weights = match method {
112        AdaptiveThresholdMethod::Mean => vec![1.0 / (block_size * block_size) as f64; block_size],
113        AdaptiveThresholdMethod::Gaussian => gaussian_weights(block_size),
114    };
115    let radius = (block_size / 2) as isize;
116    let mut output = Vec::with_capacity(input.width() * input.height());
117    for y in 0..input.height() {
118        for x in 0..input.width() {
119            let mut local = 0.0;
120            for ky in 0..block_size {
121                for kx in 0..block_size {
122                    let pixel = fetch(
123                        input,
124                        x as isize + kx as isize - radius,
125                        y as isize + ky as isize - radius,
126                        border,
127                    )[0];
128                    let weight = match method {
129                        AdaptiveThresholdMethod::Mean => weights[kx],
130                        AdaptiveThresholdMethod::Gaussian => weights[kx] * weights[ky],
131                    };
132                    local += f64::from(pixel) * weight;
133                }
134            }
135            let source = i64::from(input.get(x, y).expect("coordinate in bounds")[0]);
136            let local = cv_round(local).clamp(0, 255);
137            let delta = match threshold_type {
138                ThresholdType::Binary => c.ceil() as i64,
139                ThresholdType::BinaryInv => c.floor() as i64,
140                _ => unreachable!("adaptive threshold type checked above"),
141            };
142            let selected = match threshold_type {
143                ThresholdType::Binary => source - local > -delta,
144                ThresholdType::BinaryInv => source - local <= -delta,
145                _ => unreachable!("adaptive threshold type checked above"),
146            };
147            output.push(if selected { max_value } else { 0 });
148        }
149    }
150    Ok(Image::try_new_with_metadata(input.width(), input.height(), output, input.metadata())?)
151}
152
153/// Computes a configurable one-channel histogram, optionally through a mask.
154pub fn histogram<T: PixelComponent, const CHANNELS: usize>(
155    input: ImageView<'_, T, CHANNELS>,
156    channel: usize,
157    bins: usize,
158    range: (f64, f64),
159    mask: Option<ImageView<'_, u8, 1>>,
160) -> VisionResult<Vec<u64>> {
161    if channel >= CHANNELS {
162        return Err(VisionError::InvalidParameter(format!(
163            "histogram channel {channel} is outside {CHANNELS} channels"
164        )));
165    }
166    if bins == 0 {
167        return Err(VisionError::InvalidParameter("histogram bins must be non-zero".into()));
168    }
169    if !range.0.is_finite() || !range.1.is_finite() || range.1 <= range.0 {
170        return Err(VisionError::InvalidParameter(
171            "histogram range must be finite and increasing".into(),
172        ));
173    }
174    if let Some(mask) = mask {
175        if mask.width() != input.width() || mask.height() != input.height() {
176            return Err(VisionError::ShapeMismatch(
177                "histogram mask dimensions must match the image".into(),
178            ));
179        }
180    }
181    let mut result = vec![0_u64; bins];
182    let scale = bins as f64 / (range.1 - range.0);
183    for y in 0..input.height() {
184        for x in 0..input.width() {
185            if mask.is_some_and(|mask| mask.get(x, y).expect("mask coordinate in bounds")[0] == 0) {
186                continue;
187            }
188            let value = input.get(x, y).expect("coordinate in bounds")[channel].to_f64();
189            if value >= range.0 && value < range.1 {
190                let bin = ((value - range.0) * scale).floor() as usize;
191                result[bin.min(bins - 1)] += 1;
192            }
193        }
194    }
195    Ok(result)
196}
197
198/// Computes the exact 256-bin histogram of an 8-bit grayscale image.
199#[must_use]
200pub fn histogram_u8(input: ImageView<'_, u8, 1>) -> Vec<u64> {
201    let mut result = vec![0_u64; 256];
202    for y in 0..input.height() {
203        for x in 0..input.width() {
204            result[input.get(x, y).expect("coordinate in bounds")[0] as usize] += 1;
205        }
206    }
207    result
208}
209
210/// Equalizes an 8-bit grayscale histogram.
211pub fn equalize_histogram(input: ImageView<'_, u8, 1>) -> VisionResult<Image<u8, 1>> {
212    let histogram = histogram_u8(input);
213    let total = input.width() * input.height();
214    if total == 0 {
215        return Ok(Image::try_new_with_metadata(0, 0, Vec::new(), input.metadata())?);
216    }
217    let first = histogram.iter().position(|&count| count != 0).expect("non-empty image");
218    if histogram[first] as usize == total {
219        return Ok(Image::try_new_with_metadata(
220            input.width(),
221            input.height(),
222            vec![first as u8; total],
223            input.metadata(),
224        )?);
225    }
226    let scale = 255.0 / (total as f64 - histogram[first] as f64);
227    let mut cumulative = 0_u64;
228    let mut lookup = [0_u8; 256];
229    for (value, &count) in histogram.iter().enumerate() {
230        cumulative += count;
231        lookup[value] = if value <= first {
232            0
233        } else {
234            cv_round((cumulative - histogram[first]) as f64 * scale).clamp(0, 255) as u8
235        };
236    }
237    map_u8(input, &lookup)
238}
239
240/// Applies contrast-limited adaptive histogram equalization to grayscale u8.
241pub fn clahe(
242    input: ImageView<'_, u8, 1>,
243    clip_limit: f64,
244    tiles_x: usize,
245    tiles_y: usize,
246) -> VisionResult<Image<u8, 1>> {
247    if !clip_limit.is_finite() || clip_limit < 0.0 {
248        return Err(VisionError::InvalidParameter(
249            "CLAHE clip limit must be finite and non-negative".into(),
250        ));
251    }
252    if tiles_x == 0 || tiles_y == 0 {
253        return Err(VisionError::InvalidParameter(
254            "CLAHE tile grid dimensions must be non-zero".into(),
255        ));
256    }
257    if input.width() == 0 || input.height() == 0 {
258        return Ok(Image::try_new_with_metadata(0, 0, Vec::new(), input.metadata())?);
259    }
260    let tile_width = input.width().div_ceil(tiles_x);
261    let tile_height = input.height().div_ceil(tiles_y);
262    let tile_area = tile_width * tile_height;
263    let clip_count = if clip_limit > 0.0 {
264        ((clip_limit * tile_area as f64 / 256.0) as usize).max(1)
265    } else {
266        0
267    };
268    let mut lookups = vec![[0_u8; 256]; tiles_x * tiles_y];
269    for tile_y in 0..tiles_y {
270        for tile_x in 0..tiles_x {
271            let mut histogram = [0_usize; 256];
272            for local_y in 0..tile_height {
273                for local_x in 0..tile_width {
274                    let x = tile_x * tile_width + local_x;
275                    let y = tile_y * tile_height + local_y;
276                    let value = fetch(input, x as isize, y as isize, BorderMode::Reflect101)[0];
277                    histogram[value as usize] += 1;
278                }
279            }
280            if clip_count > 0 {
281                clip_histogram(&mut histogram, clip_count);
282            }
283            let scale = 255.0 / tile_area as f64;
284            let mut cumulative = 0_usize;
285            for (value, &count) in histogram.iter().enumerate() {
286                cumulative += count;
287                lookups[tile_y * tiles_x + tile_x][value] =
288                    cv_round(cumulative as f64 * scale).clamp(0, 255) as u8;
289            }
290        }
291    }
292    interpolate_clahe(input, &lookups, tiles_x, tiles_y, tile_width, tile_height)
293}
294
295/// Summed-area table with a zero top row and left column.
296#[derive(Clone, Debug, PartialEq)]
297pub struct IntegralImage {
298    source_width: usize,
299    source_height: usize,
300    data: Vec<f64>,
301}
302
303impl IntegralImage {
304    /// Source image width, excluding the zero border.
305    #[must_use]
306    pub const fn source_width(&self) -> usize {
307        self.source_width
308    }
309
310    /// Source image height, excluding the zero border.
311    #[must_use]
312    pub const fn source_height(&self) -> usize {
313        self.source_height
314    }
315
316    /// Integral table width (`source_width + 1`).
317    #[must_use]
318    pub const fn width(&self) -> usize {
319        self.source_width + 1
320    }
321
322    /// Integral table height (`source_height + 1`).
323    #[must_use]
324    pub const fn height(&self) -> usize {
325        self.source_height + 1
326    }
327
328    /// Row-major integral values.
329    #[must_use]
330    pub fn as_slice(&self) -> &[f64] {
331        &self.data
332    }
333
334    /// Returns the sum over half-open source rectangle `[x0,x1) × [y0,y1)`.
335    pub fn sum_region(&self, x0: usize, y0: usize, x1: usize, y1: usize) -> VisionResult<f64> {
336        if x0 > x1 || y0 > y1 || x1 > self.source_width || y1 > self.source_height {
337            return Err(VisionError::InvalidParameter(
338                "integral region is outside source dimensions".into(),
339            ));
340        }
341        let stride = self.width();
342        Ok(self.data[y1 * stride + x1] - self.data[y0 * stride + x1] - self.data[y1 * stride + x0]
343            + self.data[y0 * stride + x0])
344    }
345}
346
347/// Computes a summed-area table for one selected channel.
348pub fn integral_image<T: PixelComponent, const CHANNELS: usize>(
349    input: ImageView<'_, T, CHANNELS>,
350    channel: usize,
351) -> VisionResult<IntegralImage> {
352    if channel >= CHANNELS {
353        return Err(VisionError::InvalidParameter(format!(
354            "integral channel {channel} is outside {CHANNELS} channels"
355        )));
356    }
357    let width = input
358        .width()
359        .checked_add(1)
360        .ok_or_else(|| VisionError::InvalidDimensions("integral width overflows".into()))?;
361    let height = input
362        .height()
363        .checked_add(1)
364        .ok_or_else(|| VisionError::InvalidDimensions("integral height overflows".into()))?;
365    let len = width
366        .checked_mul(height)
367        .ok_or_else(|| VisionError::InvalidDimensions("integral allocation overflows".into()))?;
368    let mut data = vec![0.0; len];
369    for y in 0..input.height() {
370        let mut row_sum = 0.0;
371        for x in 0..input.width() {
372            row_sum += input.get(x, y).expect("coordinate in bounds")[channel].to_f64();
373            data[(y + 1) * width + x + 1] = data[y * width + x + 1] + row_sum;
374        }
375    }
376    Ok(IntegralImage { source_width: input.width(), source_height: input.height(), data })
377}
378
379fn apply_threshold(
380    value: f64,
381    threshold: f64,
382    max_value: f64,
383    threshold_type: ThresholdType,
384) -> f64 {
385    match threshold_type {
386        ThresholdType::Binary => {
387            if value > threshold {
388                max_value
389            } else {
390                0.0
391            }
392        }
393        ThresholdType::BinaryInv => {
394            if value > threshold {
395                0.0
396            } else {
397                max_value
398            }
399        }
400        ThresholdType::Truncate => value.min(threshold),
401        ThresholdType::ToZero => {
402            if value > threshold {
403                value
404            } else {
405                0.0
406            }
407        }
408        ThresholdType::ToZeroInv => {
409            if value > threshold {
410                0.0
411            } else {
412                value
413            }
414        }
415    }
416}
417
418fn otsu_from_histogram(histogram: &[u64], total: usize) -> usize {
419    if total == 0 {
420        return 0;
421    }
422    let global_mean = histogram
423        .iter()
424        .enumerate()
425        .map(|(value, &count)| value as f64 * count as f64)
426        .sum::<f64>();
427    let mut background_count = 0_u64;
428    let mut background_sum = 0.0;
429    let mut best_variance = -1.0;
430    let mut best = 0;
431    for (value, &count) in histogram.iter().enumerate() {
432        background_count += count;
433        if background_count == 0 {
434            continue;
435        }
436        let foreground_count = total as u64 - background_count;
437        if foreground_count == 0 {
438            break;
439        }
440        background_sum += value as f64 * count as f64;
441        let background_mean = background_sum / background_count as f64;
442        let foreground_mean = (global_mean - background_sum) / foreground_count as f64;
443        let difference = background_mean - foreground_mean;
444        let variance = background_count as f64 * foreground_count as f64 * difference * difference;
445        if variance > best_variance {
446            best_variance = variance;
447            best = value;
448        }
449    }
450    best
451}
452
453fn gaussian_weights(size: usize) -> Vec<f64> {
454    let exact: Option<&[f64]> = match size {
455        1 => Some(&[1.0]),
456        3 => Some(&[0.25, 0.5, 0.25]),
457        5 => Some(&[0.0625, 0.25, 0.375, 0.25, 0.0625]),
458        7 => Some(&[0.03125, 0.109375, 0.21875, 0.28125, 0.21875, 0.109375, 0.03125]),
459        9 => Some(&[
460            0.015625, 0.05078125, 0.1171875, 0.19921875, 0.234375, 0.19921875, 0.1171875,
461            0.05078125, 0.015625,
462        ]),
463        _ => None,
464    };
465    if let Some(exact) = exact {
466        return exact.to_vec();
467    }
468    let sigma = 0.3 * ((size as f64 - 1.0) * 0.5 - 1.0) + 0.8;
469    let center = (size / 2) as f64;
470    let mut weights = (0..size)
471        .map(|index| {
472            let offset = index as f64 - center;
473            (-(offset * offset) / (2.0 * sigma * sigma)).exp()
474        })
475        .collect::<Vec<_>>();
476    let sum = weights.iter().sum::<f64>();
477    weights.iter_mut().for_each(|weight| *weight /= sum);
478    weights
479}
480
481fn map_u8(input: ImageView<'_, u8, 1>, lookup: &[u8; 256]) -> VisionResult<Image<u8, 1>> {
482    let mut output = Vec::with_capacity(input.width() * input.height());
483    for y in 0..input.height() {
484        for x in 0..input.width() {
485            output.push(lookup[input.get(x, y).expect("coordinate in bounds")[0] as usize]);
486        }
487    }
488    Ok(Image::try_new_with_metadata(input.width(), input.height(), output, input.metadata())?)
489}
490
491fn clip_histogram(histogram: &mut [usize; 256], clip_limit: usize) {
492    let mut clipped = 0_usize;
493    for count in histogram.iter_mut() {
494        if *count > clip_limit {
495            clipped += *count - clip_limit;
496            *count = clip_limit;
497        }
498    }
499    let batch = clipped / 256;
500    if batch > 0 {
501        histogram.iter_mut().for_each(|count| *count += batch);
502    }
503    let residual = clipped - batch * 256;
504    let step = 256_usize.checked_div(residual).unwrap_or(1).max(1);
505    for index in (0..256).step_by(step).take(residual) {
506        histogram[index] += 1;
507    }
508}
509
510fn interpolate_clahe(
511    input: ImageView<'_, u8, 1>,
512    lookups: &[[u8; 256]],
513    tiles_x: usize,
514    tiles_y: usize,
515    tile_width: usize,
516    tile_height: usize,
517) -> VisionResult<Image<u8, 1>> {
518    let mut output = Vec::with_capacity(input.width() * input.height());
519    for y in 0..input.height() {
520        let ty = y as f64 / tile_height as f64 - 0.5;
521        let y0_raw = ty.floor() as isize;
522        let ya = ty - ty.floor();
523        let y0 = y0_raw.clamp(0, tiles_y as isize - 1) as usize;
524        let y1 = (y0_raw + 1).clamp(0, tiles_y as isize - 1) as usize;
525        for x in 0..input.width() {
526            let tx = x as f64 / tile_width as f64 - 0.5;
527            let x0_raw = tx.floor() as isize;
528            let xa = tx - tx.floor();
529            let x0 = x0_raw.clamp(0, tiles_x as isize - 1) as usize;
530            let x1 = (x0_raw + 1).clamp(0, tiles_x as isize - 1) as usize;
531            let value = input.get(x, y).expect("coordinate in bounds")[0] as usize;
532            let top = f64::from(lookups[y0 * tiles_x + x0][value]) * (1.0 - xa)
533                + f64::from(lookups[y0 * tiles_x + x1][value]) * xa;
534            let bottom = f64::from(lookups[y1 * tiles_x + x0][value]) * (1.0 - xa)
535                + f64::from(lookups[y1 * tiles_x + x1][value]) * xa;
536            output.push(cv_round(top * (1.0 - ya) + bottom * ya).clamp(0, 255) as u8);
537        }
538    }
539    Ok(Image::try_new_with_metadata(input.width(), input.height(), output, input.metadata())?)
540}
541
542fn cv_round(value: f64) -> i64 {
543    let floor = value.floor();
544    let fraction = value - floor;
545    if fraction < 0.5 {
546        floor as i64
547    } else if fraction > 0.5 {
548        floor as i64 + 1
549    } else {
550        let base = floor as i64;
551        if base % 2 == 0 {
552            base
553        } else {
554            base + 1
555        }
556    }
557}
558
559#[cfg(test)]
560mod tests {
561    use super::{
562        adaptive_threshold, clahe, equalize_histogram, histogram, integral_image,
563        otsu_threshold_u8, threshold, AdaptiveThresholdMethod, ThresholdType,
564    };
565    use crate::BorderMode;
566    use spatialrust_image::{Image, ImageRegion};
567
568    #[test]
569    fn fixed_threshold_types_match_definitions_on_roi() {
570        let image = Image::<u16, 1>::try_new(5, 1, vec![0, 10, 20, 30, 40]).unwrap();
571        let roi = image.view().subview(ImageRegion::new(1, 0, 3, 1)).unwrap();
572        assert_eq!(
573            threshold(roi, 20.0, 100.0, ThresholdType::Binary).unwrap().as_slice(),
574            &[0, 0, 100]
575        );
576        assert_eq!(
577            threshold(roi, 20.0, 100.0, ThresholdType::Truncate).unwrap().as_slice(),
578            &[10, 20, 20]
579        );
580    }
581
582    #[test]
583    fn otsu_separates_bimodal_values() {
584        let image = Image::<u8, 1>::try_new(6, 1, vec![10, 10, 10, 200, 200, 200]).unwrap();
585        let (selected, output) =
586            otsu_threshold_u8(image.view(), 255, ThresholdType::Binary).unwrap();
587        assert_eq!(selected, 10);
588        assert_eq!(output.as_slice(), &[0, 0, 0, 255, 255, 255]);
589    }
590
591    #[test]
592    fn adaptive_threshold_rejects_even_blocks() {
593        let image = Image::<u8, 1>::from_pixel(3, 3, [10]).unwrap();
594        assert!(adaptive_threshold(
595            image.view(),
596            255,
597            AdaptiveThresholdMethod::Mean,
598            ThresholdType::Binary,
599            2,
600            0.0,
601            BorderMode::Replicate
602        )
603        .is_err());
604    }
605
606    #[test]
607    fn histogram_honors_channel_range_and_mask() {
608        let image =
609            Image::<f32, 2>::try_new(2, 2, vec![0.0, 0.2, 0.5, 0.7, 0.9, 1.1, 1.0, 0.4]).unwrap();
610        let mask = Image::<u8, 1>::try_new(2, 2, vec![1, 0, 1, 1]).unwrap();
611        assert_eq!(
612            histogram(image.view(), 0, 2, (0.0, 1.0), Some(mask.view())).unwrap(),
613            vec![1, 1]
614        );
615    }
616
617    #[test]
618    fn equalization_spreads_low_contrast_values() {
619        let image = Image::<u8, 1>::try_new(4, 1, vec![10, 10, 20, 30]).unwrap();
620        assert_eq!(equalize_histogram(image.view()).unwrap().as_slice(), &[0, 0, 128, 255]);
621    }
622
623    #[test]
624    fn clahe_preserves_dimensions_and_constant_input() {
625        let image = Image::<u8, 1>::from_pixel(7, 5, [42]).unwrap();
626        let output = clahe(image.view(), 2.0, 3, 2).unwrap();
627        assert_eq!((output.width(), output.height()), (7, 5));
628    }
629
630    #[test]
631    fn integral_region_matches_direct_sum() {
632        let image = Image::<f32, 1>::try_new(3, 2, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
633        let integral = integral_image(image.view(), 0).unwrap();
634        assert_eq!((integral.width(), integral.height()), (4, 3));
635        assert_eq!(integral.sum_region(1, 0, 3, 2).unwrap(), 16.0);
636    }
637}
spatialrust_vision/analysis.rs

spatialrust_vision/
analysis.rs
