Add optional normalization of SSE template matching

examples/template_matching.rs

@@ -8,11 +8,46 @@ use std::fs;
 use std::f32;
 use image::{open, GenericImage, GrayImage, Luma, Rgb, RgbImage};
 use imageproc::definitions::Image;
-use imageproc::template_matching::match_template;
+use imageproc::template_matching::{match_template, MatchTemplateMethod};
 use imageproc::map::map_colors;
 use imageproc::rect::Rect;
 use imageproc::drawing::draw_hollow_rect_mut;
 
+struct TemplateMatchingArgs {
+    input_path: PathBuf,
+    output_dir: PathBuf,
+    template_x: u32,
+    template_y: u32,
+    template_w: u32,
+    template_h: u32
+}
+
+impl TemplateMatchingArgs {
+    fn parse(args: Vec<String>) -> TemplateMatchingArgs {
+        if args.len() != 7 {
+            panic!(r#"
+Usage:
+
+     cargo run --example template_matching input_path output_dir template_x template_y template_w template_h
+
+Loads the image at input_path and extracts a region with the given location and size to use as the matching
+template. Calls match_template on the input image and this template, and saves the results to output_dir.
+"#);
+        }
+
+        let input_path = PathBuf::from(&args[1]);
+        let output_dir = PathBuf::from(&args[2]);
+        let template_x = args[3].parse().unwrap();
+        let template_y = args[4].parse().unwrap();
+        let template_w = args[5].parse().unwrap();
+        let template_h = args[6].parse().unwrap();
+
+        TemplateMatchingArgs {
+            input_path, output_dir, template_x, template_y, template_w, template_h
+        }
+    }
+}
+
 /// Convert an f32-valued image to a 8 bit depth, covering the whole
 /// available intensity range.
 fn convert_to_gray_image(image: &Image<Luma<f32>>) -> GrayImage {
@@ -50,39 +85,24 @@ fn draw_green_rect(image: &GrayImage, rect: Rect) -> RgbImage {
     color_image
 }
 
-struct TemplateMatchingArgs {
-    input_path: PathBuf,
-    output_dir: PathBuf,
-    template_x: u32,
-    template_y: u32,
-    template_w: u32,
-    template_h: u32
-}
-
-impl TemplateMatchingArgs {
-    fn parse(args: Vec<String>) -> TemplateMatchingArgs {
-        if args.len() != 7 {
-            panic!(r#"
-Usage:
-
-     cargo run --example template_matching input_path output_dir template_x template_y template_w template_h
+fn run_match_template(
+    args: &TemplateMatchingArgs,
+    image: &GrayImage,
+    template: &GrayImage,
+    method: MatchTemplateMethod) -> RgbImage {
+    // Match the template and convert to u8 depth to display
+    let result = match_template(&image, &template, method);
+    let result_scaled = convert_to_gray_image(&result);
 
-Loads the image at input_path and extracts a region with the given location and size to use as the matching
-template. Calls match_template on the input image and this template, and saves the results to output_dir.
-"#);
-        }
+    // Pad the result to the same size as the input image, to make them easier to compare
+    let mut result_padded = GrayImage::new(image.width(), image.height());
+    result_padded.copy_from(&result_scaled, args.template_w / 2, args.template_h / 2);
 
-        let input_path = PathBuf::from(&args[1]);
-        let output_dir = PathBuf::from(&args[2]);
-        let template_x = args[3].parse().unwrap();
-        let template_y = args[4].parse().unwrap();
-        let template_w = args[5].parse().unwrap();
-        let template_h = args[6].parse().unwrap();
+    // Show location the template was extracted from
+    let roi = Rect::at(args.template_x as i32, args.template_y as i32)
+        .of_size(args.template_w, args.template_h);
 
-        TemplateMatchingArgs {
-            input_path, output_dir, template_x, template_y, template_w, template_h
-        }
-    }
+    draw_green_rect(&result_padded, roi)
 }
 
 fn main() {
@@ -107,26 +127,23 @@ fn main() {
     // Extract the requested image sub-region to use as the template
     let template = copy_sub_image(&image, args.template_x, args.template_y, args.template_w, args.template_h);
 
-    // Match the template and convert to u8 depth to display
-    let result = match_template(&image, &template);
-    let result_scaled = convert_to_gray_image(&result);
-
-    // Pad the result to the same size as the input image, to make them easier to compare
-    let mut result_padded = GrayImage::new(image.width(), image.height());
-    result_padded.copy_from(&result_scaled, args.template_w / 2, args.template_h / 2);
+    // Match using all available match methods
+    let sse = run_match_template(&args, &image, &template, MatchTemplateMethod::SumOfSquaredErrors);
+    let sse_norm = run_match_template(&args, &image, &template, MatchTemplateMethod::SumOfSquaredErrorsNormalized);
 
     // Show location the template was extracted from
     let roi = Rect::at(args.template_x as i32, args.template_y as i32)
         .of_size(args.template_w, args.template_h);
 
     let image_with_roi = draw_green_rect(&image, roi);
-    let result_with_roi = draw_green_rect(&result_padded, roi);
-
+
     // Save images to output_dir
     let template_path = output_dir.join("template.png");
     template.save(&template_path).unwrap();
     let source_path = output_dir.join("image.png");
     image_with_roi.save(&source_path).unwrap();
-    let result_path = output_dir.join("result.png");
-    result_with_roi.save(&result_path).unwrap();
+    let sse_path = output_dir.join("result_sse.png");
+    sse.save(&sse_path).unwrap();
+    let sse_path = output_dir.join("result_sse_norm.png");
+    sse_norm.save(&sse_path).unwrap();
 }

src/template_matching.rs

@@ -1,9 +1,21 @@
 //! Functions for performing template matching.
 use definitions::Image;
+use rect::Rect;
+use integral_image::{integral_squared_image, sum_image_pixels};
 use image::{GenericImage, GrayImage, Luma};
 
-/// Slides a `template` over an `image` and computes the sum of squared pixel intensity
-/// differences at each point.
+/// Method used to compute the matching score between a template and an image region.
+#[derive(Copy, Clone, Debug, PartialEq, Eq)]
+pub enum MatchTemplateMethod {
+    /// Sum of the squares of the difference between image and template pixel
+    /// intensities.
+    SumOfSquaredErrors,
+    /// Divides the sum computed using `SumOfSquaredErrors` by a normalization term.
+    SumOfSquaredErrorsNormalized,
+}
+
+/// Slides a `template` over an `image` and scores the match at each point using
+/// the requested `method`.
 ///
 /// The returned image has dimensions `image.width() - template.width() + 1` by
 /// `image.height() - template.height() + 1`.
@@ -12,34 +24,64 @@ use image::{GenericImage, GrayImage, Luma};
 ///
 /// If either dimension of `template` is not strictly less than the corresponding dimension
 /// of `image`.
-pub fn match_template(image: &GrayImage, template: &GrayImage) -> Image<Luma<f32>> {
+pub fn match_template(image: &GrayImage, template: &GrayImage, method: MatchTemplateMethod) -> Image<Luma<f32>> {
     let (image_width, image_height) = image.dimensions();
     let (template_width, template_height) = template.dimensions();
 
     assert!(image_width > template_width, "image width must strictly exceed template width");
     assert!(image_height > template_height, "image height must strictly exceed template height");
 
+    let should_normalize = method == MatchTemplateMethod::SumOfSquaredErrorsNormalized;
+    let image_squared_integral = if should_normalize { Some(integral_squared_image(&image)) } else { None };
+    let template_squared_sum = if should_normalize { Some(sum_squares(&template)) } else { None };
+
     let mut result = Image::new(image_width - template_width + 1, image_height - template_height + 1);
 
     for y in 0..result.height() {
         for x in 0..result.width() {
-            let mut sse = 0f32;
+            let mut score = 0f32;
 
             for dy in 0..template_height {
                 for dx in 0..template_width {
                     let image_value = unsafe { image.unsafe_get_pixel(x + dx, y + dy)[0] as f32 };
                     let template_value = unsafe { template.unsafe_get_pixel(dx, dy)[0] as f32 };
-                    sse += (image_value - template_value).powf(2.0);
+                    score += (image_value - template_value).powf(2.0);
                 }
             }
 
-            result.put_pixel(x, y, Luma([sse]));
+            if let (&Some(ref i), &Some(t)) = (&image_squared_integral, &template_squared_sum) {
+                let region = Rect::at(x as i32, y as i32).of_size(template_width, template_height);
+                let norm = normalization_term(i, t, region);
+                score /= norm;
+            }
+
+            result.put_pixel(x, y, Luma([score]));
         }
     }
 
     result
 }
 
+fn sum_squares(template: &GrayImage) -> f32 {
+    template.iter().map(|p| *p as f32 * *p as f32).sum()
+}
+
+/// Returns the square root of the product of the sum of squares of
+/// pixel intensities in template and the provided region of image.
+fn normalization_term(
+    image_squared_integral: &Image<Luma<u32>>,
+    template_squared_sum: f32,
+    region: Rect) -> f32 {
+    let image_sum = sum_image_pixels(
+        image_squared_integral,
+        region.left() as u32,
+        region.top() as u32,
+        region.right() as u32,
+        region.bottom() as u32
+    ) as f32;
+    (image_sum * template_squared_sum).sqrt()
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -50,22 +92,22 @@ mod tests {
     #[test]
     #[should_panic]
     fn match_template_panics_if_image_width_does_not_exceed_template_width() {
-        let _ = match_template(&GrayImage::new(5, 5), &GrayImage::new(5, 4));
+        let _ = match_template(&GrayImage::new(5, 5), &GrayImage::new(5, 4), MatchTemplateMethod::SumOfSquaredErrors);
     }
 
     #[test]
     #[should_panic]
     fn match_template_panics_if_image_height_does_not_exceed_template_height() {
-        let _ = match_template(&GrayImage::new(5, 5), &GrayImage::new(4, 5));
+        let _ = match_template(&GrayImage::new(5, 5), &GrayImage::new(4, 5), MatchTemplateMethod::SumOfSquaredErrors);
     }
 
     #[test]
     fn match_template_accepts_valid_template_size() {
-        let _ = match_template(&GrayImage::new(5, 5), &GrayImage::new(4, 4));
+        let _ = match_template(&GrayImage::new(5, 5), &GrayImage::new(4, 4), MatchTemplateMethod::SumOfSquaredErrors);
     }
 
     #[test]
-    fn match_template_example() {
+    fn match_template_sum_of_squared_errors() {
         let image = gray_image!(
             1, 4, 2;
             2, 1, 3;
@@ -76,7 +118,7 @@ mod tests {
             3, 4
         );
 
-        let actual = match_template(&image, &template);
+        let actual = match_template(&image, &template, MatchTemplateMethod::SumOfSquaredErrors);
         let expected = gray_image!(type: f32,
             14.0, 14.0;
             3.0, 1.0
@@ -85,21 +127,75 @@ mod tests {
         assert_pixels_eq!(actual, expected);
     }
 
+    #[test]
+    fn match_template_sum_of_squared_errors_normalized() {
+        let image = gray_image!(
+            1, 4, 2;
+            2, 1, 3;
+            3, 3, 4
+        );
+        let template = gray_image!(
+            1, 2;
+            3, 4
+        );
+
+        let actual = match_template(&image, &template, MatchTemplateMethod::SumOfSquaredErrorsNormalized);
+        let tss = 30f32;
+        let expected = gray_image!(type: f32,
+            14.0 / (22.0 * tss).sqrt(), 14.0 / (30.0 * tss).sqrt();
+            3.0 / (23.0 * tss).sqrt(), 1.0 / (35.0 * tss).sqrt()
+        );
+
+        assert_pixels_eq!(actual, expected);
+    }
+
     macro_rules! bench_match_template {
-        ($name:ident, image_size: $s:expr, template_size: $t:expr) => {
+        ($name:ident, image_size: $s:expr, template_size: $t:expr, method: $m:expr) => {
             #[bench]
             fn $name(b: &mut Bencher) {
                 let image = gray_bench_image($s, $s);
                 let template = gray_bench_image($t, $t);
                 b.iter(|| {
-                    let result = match_template(&image, &template);
+                    let result = match_template(&image, &template, MatchTemplateMethod::SumOfSquaredErrors);
                     black_box(result);
                 })
             }
         }
     }
 
-    bench_match_template!(bench_match_template_s100_t1, image_size: 100, template_size: 1);
-    bench_match_template!(bench_match_template_s100_t4, image_size: 100, template_size: 4);
-    bench_match_template!(bench_match_template_s100_t16, image_size: 100, template_size: 16);
+    bench_match_template!(
+        bench_match_template_s100_t1_sse,
+        image_size: 100,
+        template_size: 1,
+        method: MatchTemplateMethod::SumOfSquaredErrors);
+
+    bench_match_template!(
+        bench_match_template_s100_t4_sse,
+        image_size: 100,
+        template_size: 4,
+        method: MatchTemplateMethod::SumOfSquaredErrors);
+
+    bench_match_template!(
+        bench_match_template_s100_t16_sse,
+        image_size: 100,
+        template_size: 16,
+        method: MatchTemplateMethod::SumOfSquaredErrors);
+
+        bench_match_template!(
+        bench_match_template_s100_t1_sse_norm,
+        image_size: 100,
+        template_size: 1,
+        method: MatchTemplateMethod::SumOfSquaredErrorsNormalized);
+
+    bench_match_template!(
+        bench_match_template_s100_t4_sse_norm,
+        image_size: 100,
+        template_size: 4,
+        method: MatchTemplateMethod::SumOfSquaredErrorsNormalized);
+
+    bench_match_template!(
+        bench_match_template_s100_t16_sse_norm,
+        image_size: 100,
+        template_size: 16,
+        method: MatchTemplateMethod::SumOfSquaredErrorsNormalized);
 }