Implement L2 Norm for morphology transforms (#572)

src/distance_transform.rs

@@ -22,6 +22,12 @@ pub enum Norm {
     /// Defines d((x1, y1), (x2, y2)) to be abs(x1 - x2) + abs(y1 - y2).
     /// Also known as the Manhattan or city block norm.
     L1,
+    /// Defines d((x1, y1), (x2, y2)) to be sqrt((x1 - x2)^2 + (y1 - y2)^2).
+    /// Also known as the Euclidean norm.
+    /// when working with integer distances, we take the smallest
+    /// greater integer value, also know as 'ceiling'
+    /// example : d((0,0),(1,2)) = ceil(sqrt(1+2^2)) = ceil(2.236...) = 3
+    L2,
     /// Defines d((x1, y1), (x2, y2)) to be max(abs(x1 - x2), abs(y1 - y2)).
     /// Also known as the chessboard norm.
     LInf,
@@ -60,6 +66,17 @@ pub enum Norm {
 ///
 /// assert_pixels_eq!(distance_transform(&image, Norm::L1), l1_distances);
 ///
+/// // L2 norm
+/// let l2_distances = gray_image!(
+///     3,   3,   2,   3,   3;
+///     3,   2,   1,   2,   3;
+///     2,   1,   0,   1,   2;
+///     3,   2,   1,   2,   3;
+///     3,   3,   2,   3,   3
+/// );
+///
+/// assert_pixels_eq!(distance_transform(&image, Norm::L2), l2_distances);
+///
 /// // LInf norm
 /// let linf_distances = gray_image!(
 ///     2,   2,   2,   2,   2;
@@ -70,6 +87,7 @@ pub enum Norm {
 /// );
 ///
 /// assert_pixels_eq!(distance_transform(&image, Norm::LInf), linf_distances);
+///
 /// # }
 /// ```
 pub fn distance_transform(image: &GrayImage, norm: Norm) -> GrayImage {
@@ -95,8 +113,38 @@ pub(crate) enum DistanceFrom {
 }
 
 pub(crate) fn distance_transform_impl(image: &mut GrayImage, norm: Norm, from: DistanceFrom) {
+    match norm {
+        Norm::LInf => distance_transform_impl_linf_or_l1::<true>(image, from),
+        Norm::L1 => distance_transform_impl_linf_or_l1::<false>(image, from),
+        Norm::L2 => {
+            match from {
+                DistanceFrom::Foreground => (),
+                DistanceFrom::Background => image
+                    .iter_mut()
+                    .for_each(|p| *p = if *p == 0 { 1 } else { 0 }),
+            }
+            let float_dist: ImageBuffer<Luma<f64>, Vec<f64>> =
+                euclidean_squared_distance_transform(image);
+            image
+                .iter_mut()
+                .zip(float_dist.iter())
+                .for_each(|(u, v)| *u = v.sqrt().clamp(0.0, 255.0).ceil() as u8);
+        }
+    }
+}
+
+fn distance_transform_impl_linf_or_l1<const IS_LINF: bool>(
+    image: &mut GrayImage,
+    from: DistanceFrom,
+) {
     let max_distance = Luma([min(image.width() + image.height(), 255u32) as u8]);
 
+    // We use an unsafe code block for optimisation purposes here
+    // We use the 'unsafe_get_pixel' and 'check' unsafe functions,
+    // which are faster than safe functions,
+    // and we garantee that they are used safely
+    // by making sure we are always within the bounds of the image
+
     unsafe {
         // Top-left to bottom-right
         for y in 0..image.height() {
@@ -120,7 +168,7 @@ pub(crate) fn distance_transform_impl(image: &mut GrayImage, norm: Norm, from: D
                 if y > 0 {
                     check(image, x, y, x, y - 1);
 
-                    if norm == Norm::LInf {
+                    if IS_LINF {
                         if x > 0 {
                             check(image, x, y, x - 1, y - 1);
                         }
@@ -142,7 +190,7 @@ pub(crate) fn distance_transform_impl(image: &mut GrayImage, norm: Norm, from: D
                 if y < image.height() - 1 {
                     check(image, x, y, x, y + 1);
 
-                    if norm == Norm::LInf {
+                    if IS_LINF {
                         if x < image.width() - 1 {
                             check(image, x, y, x + 1, y + 1);
                         }
@@ -615,6 +663,9 @@ mod tests {
     bench_distance_transform!(bench_distance_transform_l1_10, Norm::L1, side: 10);
     bench_distance_transform!(bench_distance_transform_l1_100, Norm::L1, side: 100);
     bench_distance_transform!(bench_distance_transform_l1_200, Norm::L1, side: 200);
+    bench_distance_transform!(bench_distance_transform_l2_10, Norm::L2, side: 10);
+    bench_distance_transform!(bench_distance_transform_l2_100, Norm::L2, side: 100);
+    bench_distance_transform!(bench_distance_transform_l2_200, Norm::L2, side: 200);
     bench_distance_transform!(bench_distance_transform_linf_10, Norm::LInf, side: 10);
     bench_distance_transform!(bench_distance_transform_linf_100, Norm::LInf, side: 100);
     bench_distance_transform!(bench_distance_transform_linf_200, Norm::LInf, side: 200);