ci: fix ci

Cargo.toml

@@ -2,7 +2,7 @@
 name = "imageproc"
 version = "0.26.0"
 authors = ["theotherphil"]
-rust-version = "1.82.0"
+rust-version = "1.85.1"
 edition = "2021"
 license = "MIT"
 description = "Image processing operations"
@@ -46,6 +46,10 @@ assert_approx_eq = "1.1.0"
 proptest = "1.4.0"
 wasm-bindgen-test = "0.3.38"
 
+[build-dependencies]
+cargo_metadata = "0.23"
+semver = "1.0.27"
+
 [package.metadata.docs.rs]
 # See https://github.com/image-rs/imageproc/issues/358
 # all-features = true

build.rs

@@ -0,0 +1,22 @@
+use cargo_metadata::MetadataCommand;
+use semver::Version;
+
+fn main() {
+    println!("cargo::rustc-check-cfg=cfg(lagacy_image)");
+    println!("cargo:rerun-if-changed=Cargo.lock");
+
+    let metadata = MetadataCommand::new()
+        .exec()
+        .expect("Failed to fetch metadata");
+
+    let image_version = &metadata
+        .packages
+        .iter()
+        .find(|p| p.name == "image")
+        .expect("image crate not found in dependencies")
+        .version;
+
+    if image_version < &Version::parse("0.25.8").unwrap() {
+        println!("cargo:rustc-cfg=lagacy_image");
+    }
+}

src/contrast.rs

@@ -46,55 +46,94 @@ pub fn adaptive_threshold(image: &GrayImage, block_radius: u32, delta: i32) -> G
     out
 }
 
-/// Returns the [Otsu threshold level] of an 8bpp image.
+/// Returns the Otsu threshold level for an 8bpp grayscale image.
 ///
-/// [Otsu threshold level]: https://en.wikipedia.org/wiki/Otsu%27s_method
+/// This implementation uses a numerically stable algorithm to ensure consistent
+/// results across different floating-point environments, such as `cargo test` vs `cargo miri test`.
+///
+/// See: [Otsu's method on Wikipedia](https://en.wikipedia.org/wiki/Otsu%27s_method)
+///
+/// # Note on Numerical Stability
+///
+/// The standard formula for inter-class variance is:
+///
+/// $$ \sigma_B^2(t) = w_b(t) w_f(t)[\mu_b(t) - \mu_f(t)]^2 $$
+///
+/// where $w$ are the weights (pixel counts) and $\mu$ are the means of the background ($b$)
+/// and foreground ($f$) classes. This formula is susceptible to "[catastrophic cancellation](https://en.wikipedia.org/wiki/Catastrophic_cancellation)"
+/// when the means $\mu_b$ and $\mu_f$ are very close, leading to a significant loss of precision.
+///
+/// This implementation uses an algebraically equivalent, but more stable formula.
+///
+/// ## Derivation
+///
+/// Let:
+/// - $W_T$ be the total number of pixels (total weight).
+/// - $S_T$ be the sum of all pixel intensity levels.
+/// - $w_b(t)$ be the background weight at threshold $t$.
+/// - $S_b(t)$ be the background sum of intensities at threshold $t$.
+///
+/// We can express the means as $\mu_b(t) = S_b(t) / w_b(t)$ and $\mu_f(t) = S_f(t) / w_f(t)$.
+/// By substituting $w_f = W_T - w_b$ and $S_f = S_T - S_b$, the term inside the square of the
+/// original formula can be rewritten:
+///
+/// $$
+/// \begin{aligned}
+/// \mu_b - \mu_f &= \frac{S_b}{w_b} - \frac{S_f}{w_f} \\
+///              &= \frac{S_b(W_T - w_b) - (S_T - S_b)w_b}{w_b w_f} \\
+///              &= \frac{S_b W_T - S_b w_b - S_T w_b + S_b w_b}{w_b w_f} \\
+///              &= \frac{S_b W_T - S_T w_b}{w_b w_f}
+/// \end{aligned}
+/// $$
+///
+/// Substituting this back into the variance formula yields:
+///
+/// $$ \sigma_B^2(t) = w_b w_f \left( \frac{S_b W_T - S_T w_b}{w_b w_f} \right)^2 = \frac{[S_b(t) W_T - S_T w_b(t)]^2}{w_b(t) w_f(t)} $$
+///
+/// This final formula avoids the direct subtraction of means and is therefore much more robust
+/// against floating-point rounding errors.
+///
+#[cfg_attr(feature = "katexit", katexit::katexit)]
 pub fn otsu_level(image: &GrayImage) -> u8 {
     let hist = histogram(image);
     let (width, height) = image.dimensions();
     let total_weight = width * height;
 
-    // Sum of all pixel intensities, to use when calculating means.
     let total_pixel_sum = hist.channels[0]
         .iter()
         .enumerate()
         .fold(0f64, |sum, (t, h)| sum + (t as u32 * h) as f64);
 
-    // Sum of all pixel intensities in the background class.
-    let mut background_pixel_sum = 0f64;
-
-    // The weight of a class (background or foreground) is
-    // the number of pixels which belong to that class at
-    // the current threshold.
+    let mut background_pixel_sum = 0u64;
     let mut background_weight = 0u32;
-    let mut foreground_weight;
 
     let mut largest_variance = 0f64;
     let mut best_threshold = 0u8;
 
     for (threshold, hist_count) in hist.channels[0].iter().enumerate() {
-        background_weight += hist_count;
-        if background_weight == 0 {
+        let hist_count = *hist_count;
+        if hist_count == 0 {
             continue;
-        };
+        }
 
-        foreground_weight = total_weight - background_weight;
+        background_weight += hist_count;
+        background_pixel_sum += (threshold as u64) * (hist_count as u64);
+
+        let foreground_weight = total_weight - background_weight;
         if foreground_weight == 0 {
             break;
-        };
-
-        background_pixel_sum += (threshold as u32 * hist_count) as f64;
-        let foreground_pixel_sum = total_pixel_sum - background_pixel_sum;
+        }
 
-        let background_mean = background_pixel_sum / (background_weight as f64);
-        let foreground_mean = foreground_pixel_sum / (foreground_weight as f64);
+        let background_weight_f = background_weight as f64;
 
-        let mean_diff_squared = (background_mean - foreground_mean).powi(2);
-        let intra_class_variance =
-            (background_weight as f64) * (foreground_weight as f64) * mean_diff_squared;
+        let inter_class_variance = (background_pixel_sum as f64 * total_weight as f64
+            - total_pixel_sum * background_weight_f)
+            .powi(2)
+            / background_weight_f
+            / foreground_weight as f64;
 
-        if intra_class_variance > largest_variance {
-            largest_variance = intra_class_variance;
+        if inter_class_variance > largest_variance {
+            largest_variance = inter_class_variance;
             best_threshold = threshold as u8;
         }
     }

src/corners.rs

@@ -615,7 +615,7 @@ mod tests {
             00, 00, 00, 00, 00, 10, 00;
             00, 00, 00, 10, 10, 00, 00);
 
-        assert_eq!(
+        assert_approx_eq!(
             intensity_centroid(&image, 3, 3, 3),
             -std::f32::consts::FRAC_PI_4
         );

src/filter/bilateral.rs

@@ -11,20 +11,20 @@ use crate::definitions::Image;
 /// This trait is used with the `bilateral_filter()` function.
 pub trait ColorDistance<P> {
     /// Returns a distance measure between the two pixels based on their colors
-    fn color_distance(&self, pixel1: &P, pixel2: &P) -> f32;
+    fn color_distance(&self, pixel1: &P, pixel2: &P) -> f64;
 }
 
 /// A gaussian function of the euclidean distance between two pixel's colors.
 ///
 /// This implements [`ColorDistance`].
 pub struct GaussianEuclideanColorDistance {
-    sigma_squared: f32,
+    sigma_squared: f64,
 }
 impl GaussianEuclideanColorDistance {
     /// Creates a new [`GaussianEuclideanColorDistance`] using a given sigma value.
     ///
     /// Internally, this is stored as sigma squared for performance.
-    pub fn new(sigma: f32) -> Self {
+    pub fn new(sigma: f64) -> Self {
         GaussianEuclideanColorDistance {
             sigma_squared: sigma.powi(2),
         }
@@ -36,13 +36,13 @@ where
     P: Pixel,
     f32: From<P::Subpixel>,
 {
-    fn color_distance(&self, pixel1: &P, pixel2: &P) -> f32 {
+    fn color_distance(&self, pixel1: &P, pixel2: &P) -> f64 {
         let euclidean_distance_squared = pixel1
             .channels()
             .iter()
             .zip(pixel2.channels().iter())
-            .map(|(c1, c2)| (f32::from(*c1) - f32::from(*c2)).powi(2))
-            .sum::<f32>();
+            .map(|(c1, c2)| (f32::from(*c1) as f64 - f32::from(*c2) as f64).powi(2))
+            .sum::<f64>();
 
         gaussian_weight(euclidean_distance_squared, self.sigma_squared)
     }
@@ -94,15 +94,16 @@ where
 pub fn bilateral_filter<I, P, C>(
     image: &I,
     radius: u8,
-    spatial_sigma: f32,
+    spatial_sigma: f64, // changed f32 to f64
     color_distance: C,
 ) -> Image<P>
 where
     I: GenericImage<Pixel = P>,
     P: Pixel,
-    C: ColorDistance<P>,
+    C: ColorDistance<P>, // Assume ColorDistance also returns f64
     <P as image::Pixel>::Subpixel: 'static,
-    f32: From<P::Subpixel> + AsPrimitive<P::Subpixel>,
+    // changed generic bounds to use f64
+    f64: From<P::Subpixel> + AsPrimitive<P::Subpixel>,
 {
     assert!(!image.width() > i32::MAX as u32);
     assert!(!image.height() > i32::MAX as u32);
@@ -118,7 +119,8 @@ where
         .map(|(w_y, w_x)| {
             //The gaussian of the euclidean spatial distance
             gaussian_weight(
-                <f32 as From<i16>>::from(w_x).powi(2) + <f32 as From<i16>>::from(w_y).powi(2),
+                // All calculations now use f64
+                <f64 as From<i16>>::from(w_x).powi(2) + <f64 as From<i16>>::from(w_y).powi(2),
                 spatial_sigma.powi(2),
             )
         })
@@ -135,23 +137,21 @@ where
             .clone()
             .cartesian_product(window_range.clone())
             .map(|(w_y, w_x)| {
-                // these casts will always be correct due to asserts made at the beginning of the
-                // function about the image width/height
-                //
-                // the subtraction will also never overflow due to the `is_empty()` assert
                 let window_y = (i32::from(w_y) + (y as i32)).clamp(0, (height as i32) - 1);
                 let window_x = (i32::from(w_x) + (x as i32)).clamp(0, (width as i32) - 1);
 
                 let (window_y, window_x) = (window_y as u32, window_x as u32);
 
                 debug_assert!(image.in_bounds(window_x, window_y));
-                // Safety: we clamped `window_x` and `window_y` to be in bounds.
                 let window_pixel = unsafe { image.unsafe_get_pixel(window_x, window_y) };
 
                 let spatial_distance = spatial_distance_lookup
                     [(window_len * (w_y + radius) + (w_x + radius)) as usize];
-                let color_distance = color_distance.color_distance(&center_pixel, &window_pixel);
-                let weight = spatial_distance * color_distance;
+
+                // Assuming color_distance now returns f64. You may need to adjust its implementation.
+                let color_distance_val =
+                    color_distance.color_distance(&center_pixel, &window_pixel);
+                let weight = spatial_distance * color_distance_val;
 
                 (weight, window_pixel)
             });
@@ -162,17 +162,19 @@ where
     Image::from_fn(width, height, bilateral_pixel_filter)
 }
 
-fn weighted_average<P>(weights_and_values: impl Iterator<Item = (f32, P)>) -> P
+fn weighted_average<P>(weights_and_values: impl Iterator<Item = (f64, P)>) -> P
+// changed f32 to f64
 where
     P: Pixel,
     <P as image::Pixel>::Subpixel: 'static,
-    f32: From<P::Subpixel> + AsPrimitive<P::Subpixel>,
+    // changed generic bounds to use f64
+    f64: From<P::Subpixel> + AsPrimitive<P::Subpixel>,
 {
     let (weights_sum, weighted_channel_sums) = weights_and_values
         .map(|(w, v)| {
             (
                 w,
-                v.channels().iter().map(|s| w * f32::from(*s)).collect_vec(),
+                v.channels().iter().map(|s| w * f64::from(*s)).collect_vec(), // changed f32 to f64
             )
         })
         .reduce(|(w1, channels1), (w2, channels2)| {
@@ -187,17 +189,22 @@ where
         })
         .expect("cannot find a weighted average given no weights and values");
 
+    // Handle the case where all weights are zero to avoid division by zero
+    // if weights_sum == 0.0 {
+    //     return *P::from_slice(&[0.as_(); P::CHANNEL_COUNT]);
+    // }
+
     let channel_averages = weighted_channel_sums.iter().map(|x| x / weights_sum);
 
     *P::from_slice(
         &channel_averages
-            .map(<f32 as AsPrimitive<P::Subpixel>>::as_)
+            .map(<f64 as AsPrimitive<P::Subpixel>>::as_) // changed f32 to f64
             .collect_vec(),
     )
 }
 
 /// Un-normalized Gaussian Weight
-fn gaussian_weight(x_squared: f32, sigma_squared: f32) -> f32 {
+fn gaussian_weight(x_squared: f64, sigma_squared: f64) -> f64 {
     (-0.5 * x_squared / sigma_squared).exp()
 }
 
@@ -236,8 +243,8 @@ mod proptests {
         fn proptest_bilateral_filter_greyscale(
             img in arbitrary_image::<Luma<u8>>(1..40, 1..40),
             radius in 0..5u8,
-            color_sigma in any::<f32>(),
-            spatial_sigma in any::<f32>(),
+            color_sigma in any::<f64>(),
+            spatial_sigma in any::<f64>(),
         ) {
             let out = bilateral_filter(&img, radius, spatial_sigma, GaussianEuclideanColorDistance::new(color_sigma));
             assert_eq!(out.dimensions(), img.dimensions());
@@ -247,8 +254,8 @@ mod proptests {
         fn proptest_bilateral_filter_rgb(
             img in arbitrary_image::<Rgb<u8>>(1..40, 1..40),
             radius in 0..5u8,
-            color_sigma in any::<f32>(),
-            spatial_sigma in any::<f32>(),
+            color_sigma in any::<f64>(),
+            spatial_sigma in any::<f64>(),
         ) {
             let out = bilateral_filter(&img, radius, spatial_sigma, GaussianEuclideanColorDistance::new(color_sigma));
             assert_eq!(out.dimensions(), img.dimensions());

src/geometry.rs

@@ -282,7 +282,7 @@ mod tests {
             arc_length(&[Point::new(1.0, 1.0), Point::new(4.0, 5.0)], false),
             5.0
         );
-        assert_eq!(
+        assert_approx_eq!(
             arc_length(
                 &[
                     Point::new(1.0, 1.0),