Improved to be about 5x faster

.cargo/config.toml

@@ -0,0 +1,2 @@
+[build]
+rustflags = ["-C", "target-cpu=native"]

.github/workflows/rust.yml

@@ -17,12 +17,12 @@ jobs:
     steps:
       - uses: actions/checkout@v4
       # rust-version from Cargo.toml
-      - name: Install Rust 1.65.0
-        uses: dtolnay/rust-toolchain@1.65.0
+      - name: Install Rust 1.82.0
+        uses: dtolnay/rust-toolchain@1.82.0
       - name: Use predefined lockfile
         run: mv Cargo.lock.MSRV Cargo.lock
       - name: Build
-        run: cargo check --locked
+        run: cargo check
 
   build:
     runs-on: ubuntu-latest

.vscode/settings.json

@@ -0,0 +1,6 @@
+{
+    "rust-analyzer.cargo.features": [
+        "libblur",
+        "rayon"
+    ]
+}

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "ssimulacra2"
-version = "0.5.1"
+version = "0.6.0"
 edition = "2021"
 description = "Rust implementation of the SSIMULACRA2 metric"
 repository = "https://github.com/rust-av/ssimulacra2"
@@ -10,24 +10,25 @@ exclude = ["test_data"]
 license = "BSD-2-Clause"
 
 # When changing MSRV: Also update README and .github/workflows/rust.yml
-rust-version = "1.65.0"
+rust-version = "1.82.0"
 
 [features]
 default = ["rayon"]
 
 [dependencies]
-num-traits = "0.2.15"
-rayon = { version = "1.5.3", optional = true }
-thiserror = "2.0.9"
-yuvxyb = "0.4.1"
+libblur = { version = "0.17", optional = true }
+num-traits = "0.2.19"
+rayon = { version = "1.10", optional = true }
+thiserror = "2.0.12"
+yuvxyb = "0.4.2"
 
 [build-dependencies.yuvxyb-math]
 version = "0.1"
 
 [dev-dependencies]
 criterion = "0.5.0"
 image = "0.25.0"
-rand = "0.8.5"
+rand = "0.9.1"
 
 [[bench]]
 name = "benches"

benches/benches.rs

@@ -2,9 +2,10 @@ use criterion::{black_box, criterion_group, criterion_main, Criterion};
 use num_traits::clamp;
 use rand::Rng;
 use ssimulacra2::{
-    compute_frame_ssimulacra2, Blur, ColorPrimaries, Frame, MatrixCoefficients, Plane,
-    TransferCharacteristic, Yuv, YuvConfig,
+    compute_frame_ssimulacra2, xyb_to_planar, Blur, BlurOperator, ColorPrimaries, Frame,
+    MatrixCoefficients, Plane, TransferCharacteristic, Yuv, YuvConfig,
 };
+use yuvxyb::{LinearRgb, Xyb};
 
 fn make_yuv(
     ss: (u8, u8),
@@ -36,10 +37,10 @@ fn make_yuv(
             ),
         ],
     };
-    let mut rng = rand::thread_rng();
+    let mut rng = rand::rng();
     for (i, plane) in data.planes.iter_mut().enumerate() {
         for val in plane.data_origin_mut().iter_mut() {
-            *val = rng.gen_range(if full_range {
+            *val = rng.random_range(if full_range {
                 0..=255
             } else if i == 0 {
                 16..=235
@@ -64,15 +65,15 @@ fn make_yuv(
 }
 
 fn distort_yuv(input: &Yuv<u8>) -> Yuv<u8> {
-    let mut rng = rand::thread_rng();
+    let mut rng = rand::rng();
     let mut planes = [
         input.data()[0].clone(),
         input.data()[1].clone(),
         input.data()[2].clone(),
     ];
     for plane in &mut planes {
         for pix in plane.data_origin_mut() {
-            *pix = clamp(i16::from(*pix) + rng.gen_range(-16..=16), 0, 255) as u8;
+            *pix = clamp(i16::from(*pix) + rng.random_range(-16..=16), 0, 255) as u8;
         }
     }
     let data: Frame<u8> = Frame { planes };
@@ -121,10 +122,169 @@ fn bench_blur(c: &mut Criterion) {
 
         // Blur the image
         let mut blur = Blur::new(width, height);
+        let mut out = std::array::from_fn(|_| vec![0.0f32; width * height]);
 
-        b.iter(|| blur.blur(black_box(&image)))
+        b.iter(|| blur.blur(black_box(&image), black_box(&mut out)));
     });
 }
 
-criterion_group!(benches, bench_ssimulacra2, bench_blur);
+fn create_linearrgb_from_image(path: &str) -> LinearRgb {
+    let img = image::open(path).unwrap();
+    let img = img.to_rgb8();
+    let (width, height) = img.dimensions();
+
+    let mut data = vec![[0.0f32; 3]; (width * height) as usize];
+    for (i, pixel) in img.pixels().enumerate() {
+        data[i] = [
+            pixel[0] as f32 / 255.0,
+            pixel[1] as f32 / 255.0,
+            pixel[2] as f32 / 255.0,
+        ];
+    }
+
+    LinearRgb::new(data, width as usize, height as usize).expect("이미지 변환 실패")
+}
+
+// downscale_by_2
+fn bench_downscale_by_2(c: &mut Criterion) {
+    // let mut group = c.benchmark_group("downscale");
+    // group.measurement_time(std::time::Duration::from_secs_f32(9.2)); //Unable to complete 100 samples in 9.0s, so increased to 9.2s
+
+    c.bench_function("downscale_by_2", |b| {
+        // load the image
+        let image = create_linearrgb_from_image("test_data/tank_source.png");
+
+        b.iter(|| ssimulacra2::downscale_by_2(black_box(&image)))
+    });
+    // c.finish();
+}
+
+fn bench_image_multiply(c: &mut Criterion) {
+    c.bench_function("image_multiply", |b| {
+        let image = create_linearrgb_from_image("test_data/tank_source.png");
+        let width = image.width();
+        let height = image.height();
+
+        let img1 = xyb_to_planar(&Xyb::from(image.clone()));
+        let img2 = img1.clone();
+
+        let mut out = [
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+        ];
+
+        b.iter(|| {
+            ssimulacra2::image_multiply(black_box(&img1), black_box(&img2), black_box(&mut out))
+        })
+    });
+}
+
+fn bench_ssim_map(c: &mut Criterion) {
+    let mut group = c.benchmark_group("ssim_map");
+    group.measurement_time(std::time::Duration::from_secs(9));
+
+    group.bench_function("standard", |b| {
+        let source_image = create_linearrgb_from_image("test_data/tank_source.png");
+        let distorted_image = create_linearrgb_from_image("test_data/tank_distorted.png");
+        let width = source_image.width();
+        let height = source_image.height();
+
+        let mut img1_xyb = Xyb::from(source_image);
+        let mut img2_xyb = Xyb::from(distorted_image);
+
+        for pix in img1_xyb.data_mut().iter_mut() {
+            pix[2] = (pix[2] - pix[1]) + 0.55;
+            pix[0] = (pix[0]).mul_add(14.0, 0.42);
+            pix[1] += 0.01;
+        }
+        for pix in img2_xyb.data_mut().iter_mut() {
+            pix[2] = (pix[2] - pix[1]) + 0.55;
+            pix[0] = (pix[0]).mul_add(14.0, 0.42);
+            pix[1] += 0.01;
+        }
+
+        let img1_planar = xyb_to_planar(&img1_xyb);
+        let img2_planar = xyb_to_planar(&img2_xyb);
+
+        let mut blur = Blur::new(width, height);
+
+        let mut mul = [
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+        ];
+        let mut sigma1_sq = [
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+        ];
+        let mut sigma2_sq = [
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+        ];
+        let mut sigma12 = [
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+        ];
+        let mut mu1 = [
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+        ];
+        let mut mu2 = [
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+            vec![0.0f32; width * height],
+        ];
+
+        ssimulacra2::image_multiply(&img1_planar, &img1_planar, &mut mul);
+        blur.blur(&mul, &mut sigma1_sq).expect("blur failed");
+
+        ssimulacra2::image_multiply(&img2_planar, &img2_planar, &mut mul);
+        blur.blur(&mul, &mut sigma2_sq).expect("blur failed");
+
+        ssimulacra2::image_multiply(&img1_planar, &img2_planar, &mut mul);
+        blur.blur(&mul, &mut sigma12).expect("blur failed");
+
+        blur.blur(&img1_planar, &mut mu1).expect("blur failed");
+        blur.blur(&img2_planar, &mut mu2).expect("blur failed");
+
+        b.iter(|| {
+            ssimulacra2::ssim_map(
+                black_box(width),
+                black_box(height),
+                black_box(&mu1),
+                black_box(&mu2),
+                black_box(&sigma1_sq),
+                black_box(&sigma2_sq),
+                black_box(&sigma12),
+            )
+        });
+    });
+
+    group.finish();
+}
+
+fn bench_xyb_to_planar(c: &mut Criterion) {
+    c.bench_function("xyb_to_planar", |b| {
+        let image = create_linearrgb_from_image("test_data/tank_source.png");
+        let xyb = Xyb::from(image.clone());
+
+        b.iter(|| xyb_to_planar(black_box(&xyb)))
+    });
+}
+
+// criterion_group!(benches, bench_ssimulacra2, bench_blur);
+criterion_group!(
+    benches,
+    bench_ssimulacra2,
+    bench_blur,
+    bench_downscale_by_2,
+    bench_image_multiply,
+    bench_ssim_map,
+    bench_xyb_to_planar
+);
 criterion_main!(benches);

build.rs

@@ -121,13 +121,13 @@ fn init_recursive_gaussian(out_path: &str) -> io::Result<()> {
 
     write_const_usize(&mut out_file, "RADIUS", radius as usize)?;
 
-    write_const_f32(&mut out_file, "VERT_MUL_IN_1", n2[0] as f32)?;
-    write_const_f32(&mut out_file, "VERT_MUL_IN_3", n2[1] as f32)?;
-    write_const_f32(&mut out_file, "VERT_MUL_IN_5", n2[2] as f32)?;
+    // write_const_f32(&mut out_file, "VERT_MUL_IN_1", n2[0] as f32)?;
+    // write_const_f32(&mut out_file, "VERT_MUL_IN_3", n2[1] as f32)?;
+    // write_const_f32(&mut out_file, "VERT_MUL_IN_5", n2[2] as f32)?;
 
-    write_const_f32(&mut out_file, "VERT_MUL_PREV_1", d1[0] as f32)?;
-    write_const_f32(&mut out_file, "VERT_MUL_PREV_3", d1[1] as f32)?;
-    write_const_f32(&mut out_file, "VERT_MUL_PREV_5", d1[2] as f32)?;
+    // write_const_f32(&mut out_file, "VERT_MUL_PREV_1", d1[0] as f32)?;
+    // write_const_f32(&mut out_file, "VERT_MUL_PREV_3", d1[1] as f32)?;
+    // write_const_f32(&mut out_file, "VERT_MUL_PREV_5", d1[2] as f32)?;
 
     write_const_f32(&mut out_file, "MUL_IN_1", mul_in[0])?;
     write_const_f32(&mut out_file, "MUL_IN_3", mul_in[4])?;