Implement Kapur's algorithm for binary thresholding (#696) (#699)

src/contrast.rs

@@ -102,6 +102,74 @@ pub fn otsu_level(image: &GrayImage) -> u8 {
     best_threshold
 }
 
+/// Returns the [Kapur threshold level] of an 8bpp image. This threshold
+/// maximizes the entropy of the background and foreground.
+///
+/// [Kapur threshold level]: https://doi.org/10.1016/0734-189X(85)90125-2
+pub fn kapur_level(img: &GrayImage) -> u8 {
+    // The implementation looks different to the one you can for example find in
+    // ImageMagick, because we are using the simplification of equation (18) in
+    // the original article, which allows the computation of the total entropy
+    // without having to use nested loops. The names of the variables are taken
+    // straight from the article.
+    let hist = histogram(img);
+    let histogram = &hist.channels[0];
+    const N: usize = 256;
+
+    let total_pixels = (img.width() * img.height()) as f64;
+
+    // The p_i in the article. They describe the probability of encountering
+    // gray-level i.
+    let mut p = [0.0f64; N];
+    for i in 0..N {
+        p[i] = histogram[i] as f64 / total_pixels;
+    }
+
+    // The P_s in the article, which is the probability of encountering
+    // gray-level <= s.
+    let mut cum_p = [0.0f64; N];
+    cum_p[0] = p[0];
+    for i in 1..N {
+        cum_p[i] = cum_p[i - 1] + p[i];
+    }
+
+    // The H_s in the article. These are the entropies attached to the
+    // distributions p[0],...,p[s].
+    let mut h = [0.0f64; N];
+    if p[0] > 0.0 {
+        h[0] = -p[0] * p[0].ln();
+    }
+    for s in 1..N {
+        h[s] = if p[s] > 0.0 {
+            h[s - 1] - p[s] * p[s].ln()
+        } else {
+            h[s - 1]
+        };
+    }
+
+    let mut max_entropy = f64::MIN;
+    let mut best_threshold = 0;
+
+    for s in 0..N {
+        let pq = cum_p[s] * (1.0 - cum_p[s]);
+        if pq <= 0.0 {
+            continue;
+        }
+
+        // psi_s is the sum of the total entropy of foreground and
+        // background at threshold level s. Instead of computing them
+        // separately, we use equation (18) of the original article, which
+        // simplifies it to this:
+        let psi_s = pq.ln() + h[s] / cum_p[s] + (h[255] - h[s]) / (1.0 - cum_p[s]);
+        if psi_s > max_entropy {
+            max_entropy = psi_s;
+            best_threshold = s;
+        }
+    }
+
+    best_threshold as u8
+}
+
 /// Options for how to treat the threshold value in [`threshold`] and [`threshold_mut`].
 pub enum ThresholdType {
     /// `dst(x,y) = if src(x,y) > threshold { 255 } else { 0 }`
@@ -525,6 +593,13 @@ mod tests {
         GrayImage::from_pixel(width, height, Luma([intensity]))
     }
 
+    #[test]
+    fn test_kapur_constant() {
+        assert_eq!(kapur_level(&constant_image(10, 10, 0)), 0);
+        assert_eq!(kapur_level(&constant_image(10, 10, 128)), 0);
+        assert_eq!(kapur_level(&constant_image(10, 10, 255)), 0);
+    }
+
     #[test]
     fn test_otsu_constant() {
         // Variance is 0 at any threshold, and we
@@ -667,6 +742,15 @@ mod benches {
         });
     }
 
+    #[bench]
+    fn bench_kapur_level(b: &mut Bencher) {
+        let image = gray_bench_image(200, 200);
+        b.iter(|| {
+            let level = kapur_level(&image);
+            black_box(level);
+        });
+    }
+
     #[bench]
     fn bench_stretch_contrast(b: &mut Bencher) {
         let image = gray_bench_image(200, 200);