diff --git a/README.md b/README.md index 9b52d0d..91d282e 100644 --- a/README.md +++ b/README.md @@ -7,14 +7,18 @@ documentation (library usage, speed testing, MOS correlation, etc). ## Usage -Compilation requires [Zig](https://ziglang.org/) ≥0.16.0 & a macOS, Linux, or +Compilation requires [Zig](https://ziglang.org/) ≥0.17.0 & a macOS, Linux, or Unix-like operating system. To compile, run: ```sh zig build --release=fast ``` -You may add `-Dflto=true` for FLTO, and `-Dstrip=true` to strip the binary. +You may add `-Dflto=true` for FLTO on ELF targets, and `-Dstrip=true` to strip +the binary. Zig 0.17 does not support this LTO build on macOS. + +See [the benchmark notes](tests/bench/README.md) for performance and CID22 +correlation results. `fmetrics` binary usage: diff --git a/build.zig b/build.zig index 5d3e44f..4734064 100644 --- a/build.zig +++ b/build.zig @@ -19,20 +19,13 @@ pub fn build(b: *std.Build) void { const optimize = b.standardOptimizeOption(.{}); const strip = b.option(bool, "strip", "strip symbols from the binary, defaults to false") orelse false; const flto = b.option(bool, "flto", "enable Link Time Optimization, defaults to false") orelse false; + if (flto and target.result.ofmt != .elf) { + const unsupported = b.addFail("Zig 0.17 LTO requires an ELF target; omit -Dflto=true on macOS"); + b.default_step.dependOn(&unsupported.step); + return; + } const options = b.addOptions(); - const zon = b.build_root.handle.readFileAlloc( - b.graph.io, - "build.zig.zon", - b.allocator, - .limited(1 << 20), - ) catch @panic("failed to read build.zig.zon"); - const marker = ".version = \""; - const start = std.mem.indexOf(u8, zon, marker) orelse - @panic("missing version"); - const value = zon[start + marker.len ..]; - const end = std.mem.indexOfScalar(u8, value, '"') orelse - @panic("invalid version"); - options.addOption([]const u8, "version", value[0..end]); + options.addOption([]const u8, "version", @import("build.zig.zon").version); // simpleimgio const simpleimgio_dep = b.dependency("simpleimgio", .{ @@ -83,6 +76,7 @@ pub fn build(b: *std.Build) void { translate_c.addIncludePath(b.path(".")); translate_c.addIncludePath(spng.getEmittedIncludeTree()); const c_module = translate_c.createModule(); + fmetrics_module.addImport("c", c_module); // 'libfmetrics.a' static lib const lib = b.addLibrary(.{ @@ -99,6 +93,7 @@ pub fn build(b: *std.Build) void { lib.root_module.pic = true; const lib_sources = [_][]const u8{ "src/fmetrics.c", + "src/ssimulacra2/ssimulacra2.c", "src/fmetrics_cvvdp.c", "src/iwssim/iwssim.c", "src/msssim/msssim.c", @@ -128,10 +123,10 @@ pub fn build(b: *std.Build) void { if (target.result.cpu.arch == .x86_64) { var avx_query = target.query; avx_query.cpu_features_add.addFeature( - @intFromEnum(std.Target.x86.Feature.avx2), + @backingInt(std.Target.x86.Feature.avx2), ); avx_query.cpu_features_add.addFeature( - @intFromEnum(std.Target.x86.Feature.fma), + @backingInt(std.Target.x86.Feature.fma), ); const avx = b.addObject(.{ .name = "cvvdp_avx2", diff --git a/build.zig.zon b/build.zig.zon index e2897c9..49180ab 100644 --- a/build.zig.zon +++ b/build.zig.zon @@ -3,24 +3,24 @@ .version = "0.1.0", .dependencies = .{ .spng = .{ - .url = "git+https://github.com/allyourcodebase/libspng.git#8858232bde47083739c3e25f676093423610ec5e", - .hash = "spng-0.0.1-T166cP8UAABcBGTOPKx7kUUyHEcu0eDvN8LnSIj3S_cW", + .url = "git+https://github.com/gianni-rosato/libspng.git?ref=main#bbd630d686262a56086bf8e9eddb8c4004f51c0c", + .hash = "spng-0.0.2-T166cHISAADR33WgZ6IfHBVOPj0gwjyqnxYxc-XTWF8Z", }, .lcms2 = .{ .url = "git+https://github.com/allyourcodebase/lcms2.git#66c7c0ef0a9bd29a5f78aeff797bf541d0b91624", .hash = "lcms2-2.18.0-CNsMs50PAAAbhEQvq2uI5_O2vqgvYXvcDKwFw8QbzKSk", }, .simpleimgio = .{ - .url = "git+https://github.com/gianni-rosato/simpleimgio.git#6d9b55ead31a181f82595ffaed05986476510126", - .hash = "simpleimgio-0.0.0-oHTuOdP9AADo0x_m-pBpu0Fi6GnqgVvbN_l2Ul4HbViE", + .url = "git+https://github.com/gianni-rosato/simpleimgio.git#98473d19f243eefd4ef29a2149af47bc072c2484", + .hash = "simpleimgio-0.0.1-oHTuOfMbAQBWcIe2wQBOoIIrItr8a-kBVowFbAaZkwXy", }, .fcvvdp = .{ - .url = "git+https://github.com/halidecx/fcvvdp.git#98bd21a7a876be29ba06a29325b03624047ee850", - .hash = "fcvvdp-0.4.0-DualWPL1BAD7bF7Jair50v5Se-G-VluA9Y2FElmJZe2y", + .url = "git+https://github.com/halidecx/fcvvdp.git#ab4d41df02138f5573c08e05e8e33c8ef397e613", + .hash = "fcvvdp-0.4.1-DualWIn4BAD_tcJHKx8OI1FWnSFpfeIm6dljEQHaYIFU", }, }, .fingerprint = 0xcc4003cccc882e16, - .minimum_zig_version = "0.16.0", + .minimum_zig_version = "0.17.0", .paths = .{ "build.zig", "build.zig.zon", diff --git a/pyproject.toml b/pyproject.toml index 25dfa16..72ba81a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [build-system] -requires = ["setuptools>=77", "wheel>=0.45", "ziglang==0.16.0"] +requires = ["setuptools>=77", "wheel>=0.45", "ziglang==0.17.0"] build-backend = "setuptools.build_meta" [project] diff --git a/src/fmetrics.zig b/src/fmetrics.zig index 7c9578a..dc1c8b1 100644 --- a/src/fmetrics.zig +++ b/src/fmetrics.zig @@ -14,9 +14,7 @@ // limitations under the License. const std = @import("std"); const build_opts = @import("build_opts"); -const c = @cImport({ - @cInclude("fmetrics.h"); -}); +const c = @import("c"); pub const version = build_opts.version; diff --git a/src/ssimulacra2/err_map.zig b/src/ssimulacra2/err_map.zig deleted file mode 100644 index 2a24f38..0000000 --- a/src/ssimulacra2/err_map.zig +++ /dev/null @@ -1,145 +0,0 @@ -// Copyright © 2026, Halide Compression, LLC. -// All Rights Reserved. -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at: -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. -const std = @import("std"); -const math = std.math; - -const Color = extern union { - vals: extern struct { - r: u8, - g: u8, - b: u8, - a: u8 = 0xFF, - }, - concat: u32, -}; - -fn lerp(a: u8, b: u8, t: f32) u8 { - return @trunc(@as(f32, @floatFromInt(a)) * (1.0 - t) + @as(f32, @floatFromInt(b)) * t); -} - -fn turboColor(x: f32) Color { - const kRedVec4: @Vector(4, f32) = .{ 0.13572138, 4.61539260, -42.66032258, 132.13108234 }; - const kGreenVec4: @Vector(4, f32) = .{ 0.09140261, 2.19418839, 4.84296658, -14.18503333 }; - const kBlueVec4: @Vector(4, f32) = .{ 0.10667330, 12.64194608, -60.58204836, 110.36276771 }; - const kRedVec2: @Vector(2, f32) = .{ -152.94239396, 59.28637943 }; - const kGreenVec2: @Vector(2, f32) = .{ 4.27729857, 2.82956604 }; - const kBlueVec2: @Vector(2, f32) = .{ -89.90310912, 27.34824973 }; - - const clamped_x = @max(0.0, @min(1.0, x)); - const v4: @Vector(4, f32) = .{ - 1.0, - clamped_x, - clamped_x * clamped_x, - clamped_x * clamped_x * clamped_x, - }; - const v2 = [2]f32{ - v4[2] * v4[2], - v4[3] * v4[2], - }; - const r: f32 = @max(0.0, @min(1.0, @reduce(.Add, kRedVec4 * v4) + @reduce(.Add, kRedVec2 * v2))); - const g: f32 = @max(0.0, @min(1.0, @reduce(.Add, kGreenVec4 * v4) + @reduce(.Add, kGreenVec2 * v2))); - const b: f32 = @max(0.0, @min(1.0, @reduce(.Add, kBlueVec4 * v4) + @reduce(.Add, kBlueVec2 * v2))); - - return Color{ - .vals = .{ - .r = @round(r * 255.0), - .g = @round(g * 255.0), - .b = @round(b * 255.0), - .a = 0xFF, - }, - }; -} - -pub const TURBO_MAP = blk: { - @setEvalBranchQuota(2000); - var map: [256]u32 = undefined; - for (0..256) |i| { - const x = @as(f32, @floatFromInt(i)) / 255.0; - map[i] = turboColor(x).concat; - } - break :blk map; -}; - -pub fn generateErrorMap(error_accum: []const f32, error_map: []u32, stride: u32, w: u32, h: u32) void { - for (0..h) |y| { - const row_offset = y * stride; - const out_row_offset = y * w; - for (0..w) |x| { - var ssim = error_accum[row_offset + x]; - - ssim = ssim * 0.9562382616834844; - ssim = 2.326765642916932 * ssim - 0.020884521182843837 * ssim * ssim + 6.248496625763138e-05 * ssim * ssim * ssim; - - if (ssim > 0.0) - ssim = 100.0 - 10.0 * math.pow(f32, ssim, 0.6276336467831387) - else - ssim = 100.0; - - ssim = 1.0 - ssim / 100.0; - - const value: i32 = @trunc(255.0 * @max(0.0, @min(ssim, 1.0))); - error_map[out_row_offset + x] = TURBO_MAP[@intCast(value)]; - } - } -} - -inline fn bilinearSample(src: []const f32, fx: f32, fy: f32, stride: u32, w: u32, h: u32) f32 { - const x_scaled = fx * @as(f32, @floatFromInt(w)); - const y_scaled = fy * @as(f32, @floatFromInt(h)); - - const x_floor = @floor(x_scaled); - const y_floor = @floor(y_scaled); - - const fx_frac = x_scaled - x_floor; - const fy_frac = y_scaled - y_floor; - - const ix0: i32 = @trunc(@max(0.0, @min(x_floor, @as(f32, @floatFromInt(w - 1))))); - const iy0: i32 = @trunc(@max(0.0, @min(y_floor, @as(f32, @floatFromInt(h - 1))))); - const ix1: i32 = @trunc(@max(0.0, @min(x_floor + 1.0, @as(f32, @floatFromInt(w - 1))))); - const iy1: i32 = @trunc(@max(0.0, @min(y_floor + 1.0, @as(f32, @floatFromInt(h - 1))))); - - const f00 = src[@as(usize, @intCast(iy0)) * stride + @as(usize, @intCast(ix0))]; - const f10 = src[@as(usize, @intCast(iy0)) * stride + @as(usize, @intCast(ix1))]; - const f01 = src[@as(usize, @intCast(iy1)) * stride + @as(usize, @intCast(ix0))]; - const f11 = src[@as(usize, @intCast(iy1)) * stride + @as(usize, @intCast(ix1))]; - - const lerp0 = f00 * (1.0 - fx_frac) + f10 * fx_frac; - const lerp1 = f01 * (1.0 - fx_frac) + f11 * fx_frac; - - return lerp0 * (1.0 - fy_frac) + lerp1 * fy_frac; -} - -pub fn upscaleAndAccumulate( - src: []const f32, - dst: []f32, - src_stride: u32, - src_w: u32, - src_h: u32, - dst_stride: u32, - dst_w: u32, - dst_h: u32, -) void { - const h_scale: f32 = 1.0 / @as(f32, @floatFromInt(dst_h)); - const w_scale: f32 = 1.0 / @as(f32, @floatFromInt(dst_w)); - for (0..dst_h) |y| { - const fy = @as(f32, @floatFromInt(y)) * h_scale; - const row_offset = y * dst_stride; - for (0..dst_w) |x| { - const fx = @as(f32, @floatFromInt(x)) * w_scale; - const value = bilinearSample(src, fx, fy, src_stride, src_w, src_h); - dst[row_offset + x] += value; - } - } -} diff --git a/src/ssimulacra2/linearlight.zig b/src/ssimulacra2/linearlight.zig deleted file mode 100644 index a67deda..0000000 --- a/src/ssimulacra2/linearlight.zig +++ /dev/null @@ -1,52 +0,0 @@ -// Copyright © 2026, Halide Compression, LLC. -// All Rights Reserved. -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at: -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. -const std = @import("std"); -const math = std.math; - -// sRGB to Linear Light scaling for proper metric computation - -fn makeSRGBToLinearLUT() [256]f16 { - @setEvalBranchQuota(std.math.maxInt(u32)); - var lut: [256]f16 = undefined; - for (0..256) |i| { - const c: f32 = @as(f32, @floatFromInt(i)) / 255.0; - const linear = if (c <= 0.04045) c / 12.92 else math.pow(f32, (c + 0.055) / 1.055, 2.4); - lut[i] = @floatCast(linear); - } - return lut; -} - -const SRGB_LUT: [256]f16 = makeSRGBToLinearLUT(); - -pub fn sRGBInterleavedToPlanarLinear( - src: []const u8, - dst_planes: [3][]f16, - width: u32, - height: u32, - channels: u32, -) void { - const w = @as(usize, width); - const h = @as(usize, height); - for (0..h) |y| - for (0..w) |x| { - const idx = (y * w + x) * @as(usize, channels); - const r = SRGB_LUT[src[idx + 0]]; - const g = SRGB_LUT[src[idx + 1]]; - const b = SRGB_LUT[src[idx + 2]]; - dst_planes[0][y * w + x] = r; - dst_planes[1][y * w + x] = g; - dst_planes[2][y * w + x] = b; - }; -} diff --git a/src/ssimulacra2/ssimulacra2.c b/src/ssimulacra2/ssimulacra2.c new file mode 100644 index 0000000..e0bfc3b --- /dev/null +++ b/src/ssimulacra2/ssimulacra2.c @@ -0,0 +1,468 @@ +#include "../common/mem.h" +#include "../common/color.h" +#include "../common/util.h" +#include + +static const double weight[108] = { + 0.0, + 0.0007376606707406586, + 0.0, + 0.0, + 0.0007793481682867309, + 0.0, + 0.0, + 0.0004371155730107379, + 0.0, + 1.1041726426657346, + 0.00066284834129271, + 0.00015231632783718752, + 0.0, + 0.0016406437456599754, + 0.0, + 1.8422455520539298, + 11.441172603757666, + 0.0, + 0.0007989109436015163, + 0.000176816438078653, + 0.0, + 1.8787594979546387, + 10.94906990605142, + 0.0, + 0.0007289346991508072, + 0.9677937080626833, + 0.0, + 0.00014003424285435884, + 0.9981766977854967, + 0.00031949755934435053, + 0.0004550992113792063, + 0.0, + 0.0, + 0.0013648766163243398, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 7.466890328078848, + 0.0, + 17.445833984131262, + 0.0006235601634041466, + 0.0, + 0.0, + 6.683678146179332, + 0.00037724407979611296, + 1.027889937768264, + 225.20515300849274, + 0.0, + 0.0, + 19.213238186143016, + 0.0011401524586618361, + 0.001237755635509985, + 176.39317598450694, + 0.0, + 0.0, + 24.43300999870476, + 0.28520802612117757, + 0.0004485436923833408, + 0.0, + 0.0, + 0.0, + 34.77906344483772, + 44.835625328877896, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0008680556573291698, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0005313191874358747, + 0.0, + 0.00016533814161379112, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0004179171803251336, + 0.0017290828234722833, + 0.0, + 0.0020827005846636437, + 0.0, + 0.0, + 8.826982764996862, + 23.19243343998926, + 0.0, + 95.1080498811086, + 0.9863978034400682, + 0.9834382792465353, + 0.0012286405048278493, + 171.2667255897307, + 0.9807858872435379, + 0.0, + 0.0, + 0.0, + 0.0005130064588990679, + 0.0, + 0.00010854057858411537, +}; + +typedef struct Plane { + bool ssim, edge; + float *ring[6]; + double sum[6]; +} Plane; + +typedef struct Scale { + int w, h, y, level; + float *row, *pair, *tmp, *error; + Plane plane[3]; + struct Scale *next; +} Scale; + +static inline float root(float x) { + x = fmaxf(x, 1e-12f); + uint32_t bits; + memcpy(&bits, &x, sizeof(bits)); + bits = 0x54a21d2au - bits / 3; + float r; + memcpy(&r, &bits, sizeof(r)); + const float a = x / 3.0f; + for (int i = 0; i < 2; i++) { + const float r2 = r * r; + r = 1.3333333333333333f * r - a * (r2 * r2); + } + const float r2 = r * r; + r += (r - x * r2 * r2) / 3.0f; + return x * r * r - 0.15595420055f; +} + +static inline int reflect(int x, int n) { + if (x < 0) return -x; + if (x >= n) return 2 * n - 2 - x; + return x; +} + +static void horizontal(const float *restrict src, float *restrict dst, + int w) +{ + for (int x = 0; x < 4; x++) { + dst[x] = 0.2665599883f * src[x]; + const float k[4] = {0.2134445459f, 0.1095860824f, + 0.0360749699f, 0.0076144189f}; + for (int i = 1; i <= 4; i++) + dst[x] += k[i - 1] * + (src[reflect(x - i, w)] + src[reflect(x + i, w)]); + } + for (int x = 4; x < w - 4; x++) + dst[x] = 0.2665599883f * src[x] + + 0.2134445459f * (src[x - 1] + src[x + 1]) + + 0.1095860824f * (src[x - 2] + src[x + 2]) + + 0.0360749699f * (src[x - 3] + src[x + 3]) + + 0.0076144189f * (src[x - 4] + src[x + 4]); + for (int x = w - 4; x < w; x++) { + dst[x] = 0.2665599883f * src[x]; + const float k[4] = {0.2134445459f, 0.1095860824f, + 0.0360749699f, 0.0076144189f}; + for (int i = 1; i <= 4; i++) + dst[x] += k[i - 1] * + (src[reflect(x - i, w)] + src[reflect(x + i, w)]); + } +} + +static void vertical(const float *src, float *restrict dst, int w, + const int rows[9]) +{ + const float *r0 = src + (size_t)rows[0] * w; + const float *r1 = src + (size_t)rows[1] * w; + const float *r2 = src + (size_t)rows[2] * w; + const float *r3 = src + (size_t)rows[3] * w; + const float *r4 = src + (size_t)rows[4] * w; + const float *r5 = src + (size_t)rows[5] * w; + const float *r6 = src + (size_t)rows[6] * w; + const float *r7 = src + (size_t)rows[7] * w; + const float *r8 = src + (size_t)rows[8] * w; + for (int x = 0; x < w; x++) + dst[x] = 0.2665599883f * r4[x] + + 0.2134445459f * (r3[x] + r5[x]) + + 0.1095860824f * (r2[x] + r6[x]) + + 0.0360749699f * (r1[x] + r7[x]) + + 0.0076144189f * (r0[x] + r8[x]); +} + +static void measure(Scale *s, int y) { + const int w = s->w; + int rows[9]; + for (int k = 0; k < 9; k++) + rows[k] = reflect(y + k - 4, s->h) % 9; + float *m1 = s->tmp, *m2 = m1 + w; + float *cross = m2 + w, *square = cross + w; + for (int c = 0; c < 3; c++) { + Plane *p = &s->plane[c]; + if (!p->ssim && !p->edge) continue; + vertical(p->ring[0], m1, w, rows); + vertical(p->ring[1], m2, w, rows); + if (p->ssim) { + vertical(p->ring[2], cross, w, rows); + vertical(p->ring[3], square, w, rows); + } + const float *a = p->ring[4] + (size_t)(y % 9) * w; + const float *b = p->ring[5] + (size_t)(y % 9) * w; + const double *wt = weight + c * 36 + s->level * 6; + float dsum = 0, d4sum = 0, asum = 0, a4sum = 0; + float lsum = 0, l4sum = 0; +#if defined(__clang__) +#pragma clang loop vectorize(enable) +#endif + for (int x = 0; x < w; x++) { + float d = 0, artifact = 0, lost = 0; + if (p->ssim) { + const float diff = m1[x] - m2[x]; + const float covariance = 2.0f * + (cross[x] - m1[x] * m2[x]) + 0.0009f; + const float variance = square[x] - diff * diff; + d = fmaxf((variance + diff * diff * covariance) / + (covariance + variance), 0.0f); + const float d2 = d * d; + dsum += d; + d4sum += d2 * d2; + } + if (p->edge) { + const float n1 = fabsf(a[x] - m1[x]); + const float n2 = fabsf(b[x] - m2[x]); + const float delta = (n2 - n1) / (1.0f + n1); + artifact = fmaxf(delta, 0.0f); + lost = fmaxf(-delta, 0.0f); + const float a2 = artifact * artifact; + const float l2 = lost * lost; + asum += artifact; + a4sum += a2 * a2; + lsum += lost; + l4sum += l2 * l2; + } + if (s->error) s->error[(size_t)y * w + x] += + (float)((wt[0] + wt[3]) * d + + (wt[1] + wt[4]) * artifact + + (wt[2] + wt[5]) * lost); + } + p->sum[0] += dsum; + p->sum[1] += asum; + p->sum[2] += lsum; + p->sum[3] += d4sum; + p->sum[4] += a4sum; + p->sum[5] += l4sum; + } +} + +#define XYB_ROW(X, B) \ + for (int x = 0; x < w; x++) { \ + const float cx = root(fmaxf(0.30f * r[x] + 0.622f * g[x] + \ + 0.078f * b[x] + 0.0037930734f, 0)); \ + const float cy = root(fmaxf(0.23f * r[x] + 0.692f * g[x] + \ + 0.078f * b[x] + 0.0037930734f, 0)); \ + const float yy = (cx + cy) * 0.5f; \ + if (X) ox[x] = (cx - cy) * 7.0f + 0.42f; \ + oy[x] = yy + 0.01f; \ + if (B) ob[x] = root(fmaxf(0.24342269f * r[x] + \ + 0.20476745f * g[x] + 0.55180986f * b[x] + \ + 0.0037930734f, 0)) - yy + 0.55f; \ + } + +static void push(Scale *s) { + const int w = s->w, y = s->y; + if (s->next) { + if (!(y & 1)) memcpy(s->pair, s->row, (size_t)w * 6 * 4); + if ((y & 1) || y == s->h - 1) { + const int nw = s->next->w; + for (int c = 0; c < 6; c++) { + const float *a = s->pair + (size_t)c * w; + const float *b = s->row + (size_t)c * w; + float *out = s->next->row + (size_t)c * nw; + for (int x = 0; x < w / 2; x++) + out[x] = (a[2 * x] + a[2 * x + 1] + + b[2 * x] + b[2 * x + 1]) * 0.25f; + if (w & 1) out[nw - 1] = (a[w - 1] + b[w - 1]) * 0.5f; + } + push(s->next); + } + } + for (int img = 0; img < 2; img++) { + const float *r = s->row + (size_t)img * 3 * w; + const float *g = r + w, *b = g + w; + float *restrict ox = s->plane[0].ring[4 + img]; + float *restrict oy = s->plane[1].ring[4 + img]; + float *restrict ob = s->plane[2].ring[4 + img]; + const size_t off = (size_t)(y % 9) * w; + if (ox) ox += off; + if (oy) oy += off; + if (ob) ob += off; + if (ob) { XYB_ROW(1, 1) } + else if (ox) { XYB_ROW(1, 0) } + else { XYB_ROW(0, 0) } + } + for (int c = 0; c < 3; c++) { + Plane *p = &s->plane[c]; + if (!p->ssim && !p->edge) continue; + const size_t off = (size_t)(y % 9) * w; + float *a = p->ring[4] + off, *b = p->ring[5] + off; + horizontal(a, p->ring[0] + off, w); + horizontal(b, p->ring[1] + off, w); + if (p->ssim) { + for (int x = 0; x < w; x++) s->tmp[x] = a[x] * b[x]; + horizontal(s->tmp, p->ring[2] + off, w); + for (int x = 0; x < w; x++) { + const float v = a[x] - b[x]; + s->tmp[x] = v * v; + } + horizontal(s->tmp, p->ring[3] + off, w); + } + } + if (y >= 4) measure(s, y - 4); + if (++s->y == s->h) + for (int i = s->h - 4; i < s->h; i++) measure(s, i); +} + +static double score(double v) { + v *= 0.9562382616834844; + v = (6.248496625763138e-5 * v * v) * v + + 2.326765642916932 * v - 0.020884521182843837 * v * v; + return v > 0 ? 100.0 - 10.0 * pow(v, 0.6276336467831387) : 100.0; +} + +static float sample(const Scale *s, int x, int y, int w, int h) { + const float fx = (float)x * s->w / w, fy = (float)y * s->h / h; + const int x0 = (int)fx, y0 = (int)fy; + const int x1 = x0 + (x0 + 1 < s->w); + const int y1 = y0 + (y0 + 1 < s->h); + const float dx = fx - x0, dy = fy - y0; + const float a = s->error[(size_t)y0 * s->w + x0] * (1 - dx) + + s->error[(size_t)y0 * s->w + x1] * dx; + const float b = s->error[(size_t)y1 * s->w + x0] * (1 - dx) + + s->error[(size_t)y1 * s->w + x1] * dx; + return a * (1 - dy) + b * dy; +} + +static void render_map(Scale *s, int count, uint32_t *out, int w, int h) { + for (int y = 0; y < h; y++) for (int x = 0; x < w; x++) { + float v = 0; + for (int i = 0; i < count && i < 3; i++) + v += sample(&s[i], x, y, w, h); + const int index = (int)(255 * fmin(fmax(1 - score(v) / 100, 0), 1)); + out[(size_t)y * w + x] = TURBO_COLORMAP[index]; + } +} + +static FmetricsErr compute(FmetricsWorkspace *ws, const FmetricsImg *ref, + const FmetricsImg *dist, double *result, + uint32_t *map) +{ + if (!ws || !ref || !dist || !result || !ref->data || !dist->data) + return FMETRICS_ERR_INVALID_ARGUMENT; + if (!ref->width || !ref->height || ref->width > INT_MAX - 16 || + ref->height > INT_MAX - 16) return FMETRICS_ERR_INVALID_ARGUMENT; + if (ref->width != dist->width || ref->height != dist->height) + return FMETRICS_ERR_DIMENSION_MISMATCH; + FmetricsErr err = fmetrics_validate_rgb(ref); + if (err != FMETRICS_OK) return err; + err = fmetrics_validate_rgb(dist); + if (err != FMETRICS_OK) return err; + const int w = (int)ref->width, h = (int)ref->height; + const size_t pixels = (size_t)w * h; + if ((size_t)w > (SIZE_MAX - 4096) / 1024 || + (map && pixels > (SIZE_MAX - (size_t)w * 1024 - 4096) / 8)) + return FMETRICS_ERR_OUT_OF_MEMORY; + const size_t size = (size_t)w * 1024 + 4096 + (map ? pixels * 8 : 0); + if (!workspace_reserve(ws, size)) return FMETRICS_ERR_OUT_OF_MEMORY; + Scale scales[5] = {0}; + int count = 0, sw = w, sh = h; + for (int level = 0; level <= 4; level++) { + if (level) { + sw = (sw + 1) / 2; + sh = (sh + 1) / 2; + } + if (sw < 8 || sh < 8) break; + Scale *s = &scales[count++]; + s->w = sw; + s->h = sh; + s->level = level; + s->row = scratch_alloc(&ws->scratch, (size_t)sw * 6 * 4); + s->pair = scratch_alloc(&ws->scratch, (size_t)sw * 6 * 4); + s->tmp = scratch_alloc(&ws->scratch, (size_t)sw * 4 * 4); + if (map && level < 3) { + s->error = scratch_alloc(&ws->scratch, (size_t)sw * sh * 4); + memset(s->error, 0, (size_t)sw * sh * 4); + } + for (int c = 0; c < 3; c++) { + Plane *p = &s->plane[c]; + const double *wt = weight + c * 36 + level * 6; + p->ssim = wt[0] > 0.01 || wt[3] > 0.01; + p->edge = wt[1] > 0.01 || wt[4] > 0.01 || + wt[2] > 0.01 || wt[5] > 0.01; + if (!p->ssim && !p->edge) continue; + for (int i = 0; i < 6; i++) { + if (!p->ssim && (i == 2 || i == 3)) continue; + p->ring[i] = scratch_alloc(&ws->scratch, (size_t)sw * 9 * 4); + } + } + if (count > 1) scales[count - 2].next = s; + } + if (count) { + float lut[256]; + for (int i = 0; i < 256; i++) + lut[i] = fmetrics_srgb_linear(i / 255.0f); + const FmetricsImg *images[2] = {ref, dist}; + for (int y = 0; y < h; y++) { + for (int img = 0; img < 2; img++) { + const FmetricsImg *in = images[img]; + const uint8_t *row = (const uint8_t *)in->data + + (size_t)y * in->stride; + for (int c = 0; c < 3; c++) { + float *out = scales[0].row + (size_t)(img * 3 + c) * w; + if (fmetrics_is_linear(in)) { + for (int x = 0; x < w; x++) + memcpy(out + x, row + ((size_t)x * 3 + c) * 4, 4); + } else for (int x = 0; x < w; x++) + out[x] = lut[row[(size_t)x * 3 + c]]; + } + } + push(&scales[0]); + } + } + + double total = 0; + for (int i = 0; i < count; i++) for (int c = 0; c < 3; c++) { + const Scale *s = &scales[i]; + const double n = (double)s->w * s->h; + for (int j = 0; j < 6; j++) { + double v = s->plane[c].sum[j] / n; + if (j >= 3) v = sqrt(sqrt(v)); + total += weight[c * 36 + s->level * 6 + j] * v; + } + } + *result = score(total); + if (map) render_map(scales, count, map, w, h); + return FMETRICS_OK; +} + +FmetricsErr fmetrics_ssimu2_cmp(FmetricsWorkspace *workspace, + const FmetricsImg *reference, + const FmetricsImg *distorted, double *result) +{ + return compute(workspace, reference, distorted, result, NULL); +} + +FmetricsErr fmetrics_ssimu2_cmp_map(FmetricsWorkspace *workspace, + const FmetricsImg *reference, + const FmetricsImg *distorted, double *result, + uint32_t *error_map) +{ + if (!error_map) return FMETRICS_ERR_INVALID_ARGUMENT; + return compute(workspace, reference, distorted, result, error_map); +} diff --git a/src/ssimulacra2/ssimulacra2.zig b/src/ssimulacra2/ssimulacra2.zig deleted file mode 100644 index 2f5b096..0000000 --- a/src/ssimulacra2/ssimulacra2.zig +++ /dev/null @@ -1,1150 +0,0 @@ -// Copyright © 2026, Halide Compression, LLC. -// All Rights Reserved. -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at: -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. -const std = @import("std"); -const math = std.math; -const scl = @import("linearlight.zig"); -const err_map = @import("err_map.zig"); - -pub const Ssimu2Error = error{ - InvalidChannelCount, - OutOfMemory, -}; - -const K_SIZE = 9; -const RADIUS = 4; -const VEC_LEN = std.simd.suggestVectorLength(f32) orelse 8; -const F32Vec = @Vector(VEC_LEN, f32); -const F64Vec = @Vector(VEC_LEN, f64); - -pub const Workspace = struct { - allocator: std.mem.Allocator, - width: u32, - height: u32, - stride: u32, - plane_size: usize, - planes: []f16, - temp: []f32, - scratch: []f32, - - pub fn init(allocator: std.mem.Allocator, width: u32, height: u32) !Workspace { - var ws: Workspace = .{ - .allocator = allocator, - .width = 0, - .height = 0, - .stride = 0, - .plane_size = 0, - .planes = &[_]f16{}, - .temp = &[_]f32{}, - .scratch = &[_]f32{}, - }; - try ws.ensureCapacity(width, height); - return ws; - } - - pub fn deinit(self: *Workspace) void { - if (self.planes.len > 0) self.allocator.free(self.planes); - if (self.temp.len > 0) self.allocator.free(self.temp); - if (self.scratch.len > 0) self.allocator.free(self.scratch); - self.* = undefined; - } - - pub fn ensureCapacity(self: *Workspace, width: u32, height: u32) !void { - if (self.width == width and self.height == height and self.planes.len > 0) return; - - if (self.planes.len > 0) self.allocator.free(self.planes); - if (self.temp.len > 0) self.allocator.free(self.temp); - if (self.scratch.len > 0) self.allocator.free(self.scratch); - - const pixels = @as(usize, width) * @as(usize, height); - const total_floats: usize = pixels * 3 * 2; - const plane_alloc = try self.allocator.alignedAlloc(f16, .of(f16), total_floats); - - const wh: usize = @as(usize, width) * @as(usize, height); - const temp_alloc = try self.allocator.alignedAlloc( - f32, - .of(f32), - wh * 18, - ); - const scratch_alloc = try self.allocator.alignedAlloc( - f32, - .of(f32), - width + 64, - ); - - self.width = width; - self.height = height; - self.stride = width; - self.plane_size = pixels; - self.planes = plane_alloc; - self.temp = temp_alloc; - self.scratch = scratch_alloc; - } -}; - -pub fn computeSsimu2WithWorkspace( - workspace: *Workspace, - reference: []const u8, - distorted: []const u8, - width: u32, - height: u32, - channels: u32, - error_map: ?[]u32, -) Ssimu2Error!f64 { - try workspace.ensureCapacity(width, height); - return computeSsimu2WithBuffers( - reference, - distorted, - width, - height, - channels, - error_map, - workspace.planes, - workspace.temp, - workspace.scratch, - ); -} - -pub fn computeSsimu2WithBuffers( - reference: []const u8, - distorted: []const u8, - width: u32, - height: u32, - channels: u32, - error_map: ?[]u32, - planes: []f16, - temp: []f32, - scratch: []f32, -) Ssimu2Error!f64 { - if (channels != 3 and channels != 4) return Ssimu2Error.InvalidChannelCount; - - const pixels = @as(usize, width) * @as(usize, height); - const expected_len = pixels * @as(usize, channels); - - std.debug.assert(reference.len >= expected_len); - std.debug.assert(distorted.len >= expected_len); - std.debug.assert(planes.len >= pixels * 6); - std.debug.assert(temp.len >= pixels * 18); - std.debug.assert(scratch.len >= width + 64); - - var ref_planes: [3][]f16 = undefined; - var dist_planes: [3][]f16 = undefined; - { - var off: usize = 0; - inline for (0..3) |i| { - ref_planes[i] = planes[off .. off + pixels]; - off += pixels; - } - inline for (0..3) |i| { - dist_planes[i] = planes[off .. off + pixels]; - off += pixels; - } - } - - scl.sRGBInterleavedToPlanarLinear(reference, ref_planes, width, height, channels); - scl.sRGBInterleavedToPlanarLinear(distorted, dist_planes, width, height, channels); - - const ref_const: [3][]const f16 = .{ - ref_planes[0], ref_planes[1], ref_planes[2], - }; - const dist_const: [3][]const f16 = .{ - dist_planes[0], dist_planes[1], dist_planes[2], - }; - - return processWithScratch( - ref_const, - dist_const, - width, - width, - height, - error_map, - temp, - scratch, - ); -} - -pub fn computeSsimu2( - allocator: std.mem.Allocator, - reference: []const u8, - distorted: []const u8, - width: u32, - height: u32, - channels: u32, - error_map: ?[]u32, -) Ssimu2Error!f64 { - var workspace = try Workspace.init(allocator, width, height); - defer workspace.deinit(); - return computeSsimu2WithWorkspace(&workspace, reference, distorted, width, height, channels, error_map); -} - -inline fn loadVec(src: []const f32, index: usize) F32Vec { - return @as(*align(1) const F32Vec, @ptrCast(src.ptr + index)).*; -} - -inline fn storeVec(dst: []f32, index: usize, value: F32Vec) void { - @as(*align(1) F32Vec, @ptrCast(dst.ptr + index)).* = value; -} - -fn deinterleaveMask(comptime offset: usize) [VEC_LEN]i32 { - var mask: [VEC_LEN]i32 = undefined; - for (0..VEC_LEN) |i| { - const index = i * 2 + offset; - if (index < VEC_LEN) { - mask[i] = index; - continue; - } - mask[i] = ~@as(i32, @intCast(index - VEC_LEN)); - } - return mask; -} - -inline fn multiply(src1: []const f32, src2: []const f32, dst: []f32, stride: u32, w: u32, h: u32) void { - var y: u32 = 0; - while (y < h) : (y += 1) { - const row: usize = @as(usize, y) * stride; - var x: usize = 0; - while (x + VEC_LEN <= w) : (x += VEC_LEN) { - storeVec(dst, row + x, loadVec(src1, row + x) * - loadVec(src2, row + x)); - } - while (x < w) : (x += 1) { - dst[row + x] = src1[row + x] * src2[row + x]; - } - } -} - -inline fn addSquare( - src1: []const f32, - src2: []const f32, - dst: []f32, - stride: u32, - w: u32, - h: u32, -) void { - var y: u32 = 0; - while (y < h) : (y += 1) { - const row: usize = @as(usize, y) * stride; - var x: usize = 0; - while (x + VEC_LEN <= w) : (x += VEC_LEN) { - const v = loadVec(src1, row + x) + loadVec(src2, row + x); - storeVec(dst, row + x, v * v); - } - while (x < w) : (x += 1) { - const v = src1[row + x] + src2[row + x]; - dst[row + x] = v * v; - } - } -} - -fn blurH(srcp: []const f32, dstp: []f32, kernel: [K_SIZE]f32, w: i32) void { - var j: i32 = 0; - while (j < @min(w, RADIUS)) : (j += 1) { - const dist_from_right: i32 = w - 1 - j; - var sum: f32 = 0.0; - var k: i32 = 0; - while (k < RADIUS) : (k += 1) { - const idx: i32 = if (j < RADIUS - k) @min(RADIUS - k - j, w - 1) else (j - RADIUS + k); - sum += kernel[@intCast(k)] * srcp[@intCast(idx)]; - } - k = RADIUS; - while (k < K_SIZE) : (k += 1) { - const idx: i32 = if (dist_from_right < k - RADIUS) (j - @min(k - RADIUS - dist_from_right, j)) else (j - RADIUS + k); - sum += kernel[@intCast(k)] * srcp[@intCast(idx)]; - } - dstp[@intCast(j)] = sum; - } - - const interior_end: usize = @intCast(w - @min(w, RADIUS)); - var ju: usize = RADIUS; - while (ju + VEC_LEN <= interior_end) : (ju += VEC_LEN) { - var sum: F32Vec = @splat(0.0); - inline for (0..K_SIZE) |k| - sum += @as(F32Vec, @splat(kernel[k])) * - loadVec(srcp, ju - RADIUS + k); - storeVec(dstp, ju, sum); - } - - j = @intCast(ju); - while (j < w - @min(w, RADIUS)) : (j += 1) { - var sum: f32 = kernel[RADIUS] * srcp[@intCast(j)]; - var k: i32 = 1; - while (k <= RADIUS) : (k += 1) { - const left = srcp[@intCast(j - k)]; - const right = srcp[@intCast(j + k)]; - sum += kernel[@intCast(RADIUS + k)] * (left + right); - } - dstp[@intCast(j)] = sum; - } - - j = @max(RADIUS, w - @min(w, RADIUS)); - while (j < w) : (j += 1) { - const dist_from_right: i32 = w - 1 - j; - var sum: f32 = 0.0; - var k: i32 = 0; - while (k < RADIUS) : (k += 1) { - const idx: i32 = if (j < RADIUS - k) @min(RADIUS - k - j, w - 1) else (j - RADIUS + k); - sum += kernel[@intCast(k)] * srcp[@intCast(idx)]; - } - k = RADIUS; - while (k < K_SIZE) : (k += 1) { - const idx: i32 = if (dist_from_right < k - RADIUS) (j - @min(k - RADIUS - dist_from_right, j)) else (j - RADIUS + k); - sum += kernel[@intCast(k)] * srcp[@intCast(idx)]; - } - dstp[@intCast(j)] = sum; - } -} - -inline fn blurV(src: anytype, dstp: []f32, kernel: [K_SIZE]f32, w: u32) void { - var j: usize = 0; - while (j + VEC_LEN <= w) : (j += VEC_LEN) { - var accum: F32Vec = @splat(0.0); - inline for (0..K_SIZE) |k| - accum = @mulAdd( - F32Vec, - @splat(kernel[k]), - loadVec(src[k], j), - accum, - ); - storeVec(dstp, j, accum); - } - while (j < w) : (j += 1) { - var accum: f32 = kernel[RADIUS] * - @as(f32, @floatCast(src[RADIUS][j])); - var k: u32 = 1; - while (k <= RADIUS) : (k += 1) { - const up = @as(f32, @floatCast(src[RADIUS - k][j])); - const down = @as(f32, @floatCast(src[RADIUS + k][j])); - accum += kernel[RADIUS + k] * (up + down); - } - dstp[j] = accum; - } -} - -inline fn blur(src: []const f32, dst: []f32, stride: u32, w: u32, h: u32, tmp_row: []f32) void { - const kernel = [K_SIZE]f32{ - 0.0076144188642501831054687500, - 0.0360749699175357818603515625, - 0.1095860823988914489746093750, - 0.2134445458650588989257812500, - 0.2665599882602691650390625000, - 0.2134445458650588989257812500, - 0.1095860823988914489746093750, - 0.0360749699175357818603515625, - 0.0076144188642501831054687500, - }; - var i: i32 = 0; - const ih: i32 = @bitCast(h); - while (i < ih) : (i += 1) { - const ui: u32 = @bitCast(i); - var srcp_rows: [K_SIZE][]const f32 = undefined; - const dstp_row: []f32 = dst[(ui * stride)..]; - const dist_from_bottom: i32 = ih - 1 - i; - - var k: i32 = 0; - while (k < RADIUS) : (k += 1) { - const row: i32 = if (i < RADIUS - k) (@min(RADIUS - k - i, ih - 1)) else (i - RADIUS + k); - const urow: u32 = @bitCast(row); - srcp_rows[@intCast(k)] = src[(urow * stride)..]; - } - k = RADIUS; - while (k < K_SIZE) : (k += 1) { - const row: i32 = if (dist_from_bottom < k - RADIUS) - (i - @min(k - RADIUS - dist_from_bottom, i)) - else - (i - RADIUS + k); - const urow: u32 = @bitCast(row); - srcp_rows[@intCast(k)] = src[(urow * stride)..]; - } - - blurV(srcp_rows, tmp_row, kernel, w); - blurH(tmp_row, dstp_row, kernel, @intCast(w)); - } -} - -const K_D0: f32 = 0.0037930734; -const K_D1: f32 = std.math.lossyCast(f32, math.cbrt(K_D0)); - -const V00: f32 = 0.0; -const V05: f32 = 0.5; -const V10: f32 = 1.0; -const V001: f32 = 0.01; -const V055: f32 = 0.55; -const V042: f32 = 0.42; -const V140: f32 = 14.0; - -const weight = [108]f64{ - 0.0, - 0.0007376606707406586, - 0.0, - 0.0, - 0.0007793481682867309, - 0.0, - 0.0, - 0.0004371155730107379, - 0.0, - 1.1041726426657346, - 0.00066284834129271, - 0.00015231632783718752, - 0.0, - 0.0016406437456599754, - 0.0, - 1.8422455520539298, - 11.441172603757666, - 0.0, - 0.0007989109436015163, - 0.000176816438078653, - 0.0, - 1.8787594979546387, - 10.94906990605142, - 0.0, - 0.0007289346991508072, - 0.9677937080626833, - 0.0, - 0.00014003424285435884, - 0.9981766977854967, - 0.00031949755934435053, - 0.0004550992113792063, - 0.0, - 0.0, - 0.0013648766163243398, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 7.466890328078848, - 0.0, - 17.445833984131262, - 0.0006235601634041466, - 0.0, - 0.0, - 6.683678146179332, - 0.00037724407979611296, - 1.027889937768264, - 225.20515300849274, - 0.0, - 0.0, - 19.213238186143016, - 0.0011401524586618361, - 0.001237755635509985, - 176.39317598450694, - 0.0, - 0.0, - 24.43300999870476, - 0.28520802612117757, - 0.0004485436923833408, - 0.0, - 0.0, - 0.0, - 34.77906344483772, - 44.835625328877896, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0008680556573291698, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0005313191874358747, - 0.0, - 0.00016533814161379112, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0004179171803251336, - 0.0017290828234722833, - 0.0, - 0.0020827005846636437, - 0.0, - 0.0, - 8.826982764996862, - 23.19243343998926, - 0.0, - 95.1080498811086, - 0.9863978034400682, - 0.9834382792465353, - 0.0012286405048278493, - 171.2667255897307, - 0.9807858872435379, - 0.0, - 0.0, - 0.0, - 0.0005130064588990679, - 0.0, - 0.00010854057858411537, -}; - -inline fn getWeight(c: i32, scale: i32, map: i32, norm: i32) f64 { - std.debug.assert(c >= 0 and c < 3); - std.debug.assert(scale >= 0 and scale < 6); - std.debug.assert(map >= 0 and map < 3); - std.debug.assert(norm >= 0 and norm < 2); - const idx = 36 * c + 6 * scale + 3 * norm + map; - std.debug.assert(idx >= 0 and idx < 108); - return weight[@intCast(idx)]; -} - -const SkipInfo = struct { - ssim: bool, - artifact: bool, - detailloss: bool, - - inline fn all(self: SkipInfo) bool { - return self.ssim and self.artifact and self.detailloss; - } -}; - -const skip_table: [3][6]SkipInfo = blk: { - var table: [3][6]SkipInfo = undefined; - for (0..3) |plane| { - for (0..6) |scale| { - const base = plane * 36 + scale * 6; - table[plane][scale] = .{ - .ssim = weight[base] <= 0.01 and weight[base + 3] <= 0.01, - .artifact = weight[base + 1] <= 0.01 and - weight[base + 4] <= 0.01, - .detailloss = weight[base + 2] <= 0.01 and - weight[base + 5] <= 0.01, - }; - } - } - break :blk table; -}; - -const K_M02 = 0.078; -const K_M00 = 0.30; -const K_M01 = V10 - K_M02 - K_M00; -const K_M12 = 0.078; -const K_M10 = 0.23; -const K_M11 = V10 - K_M12 - K_M10; -const K_M20 = 0.24342269; -const K_M21 = 0.20476745; -const K_M22 = V10 - K_M20 - K_M21; - -const OPSIN_ABSORBANCE_MATRIX = [_]f32{ - K_M00, K_M01, K_M02, - K_M10, K_M11, K_M12, - K_M20, K_M21, K_M22, -}; -const OPSIN_ABSORBANCE_BIAS: f32 = K_D0; -const ABSORBANCE_BIAS: f32 = -K_D1; - -inline fn cubeRootAndAdd(x: f32, add: f32) f32 { - const k_exp_bias: i32 = 0x54800000; - const k_exp_mul: i32 = 0x002AAAAA; - const k1_3: f32 = 1.0 / 3.0; - const k4_3: f32 = 4.0 / 3.0; - const xa_3 = k1_3 * x; - const m1: i32 = @bitCast(x); - const m2 = if (m1 == 0) 0 else k_exp_bias - (m1 >> 23) * k_exp_mul; - var r: f32 = @bitCast(m2); - var i: u32 = 0; - while (i < 3) : (i += 1) { - const r2 = r * r; - r = -xa_3 * (r2 * r2) + k4_3 * r; - } - var r2 = r * r; - r = k1_3 * (-x * r2 * r2 + r) + r; - r2 = r * r; - return r2 * x + add; -} - -inline fn cubeRootVec(x: F32Vec) F32Vec { - const k1_3: F32Vec = @splat(1.0 / 3.0); - const k4_3: F32Vec = @splat(4.0 / 3.0); - const xa_3 = k1_3 * x; - const m1: @Vector(VEC_LEN, i32) = @bitCast(x); - const zero: @Vector(VEC_LEN, i32) = @splat(0); - const m2 = @select( - i32, - m1 == zero, - zero, - @as(@Vector(VEC_LEN, i32), @splat(0x54800000)) - - (m1 >> @splat(23)) * - @as(@Vector(VEC_LEN, i32), @splat(0x002AAAAA)), - ); - var r: F32Vec = @bitCast(m2); - var i: u32 = 0; - while (i < 3) : (i += 1) { - const r2 = r * r; - r = -xa_3 * (r2 * r2) + k4_3 * r; - } - var r2 = r * r; - r = k1_3 * (-x * r2 * r2 + r) + r; - r2 = r * r; - return r2 * x - @as(F32Vec, @splat(K_D1)); -} - -inline fn opsinVec( - r: F32Vec, - g: F32Vec, - b: F32Vec, - m0: f32, - m1: f32, - m2: f32, - bias: F32Vec, -) F32Vec { - const mixed_b = @mulAdd(F32Vec, @splat(m2), b, bias); - const mixed_g = @mulAdd(F32Vec, @splat(m1), g, mixed_b); - return @mulAdd(F32Vec, @splat(m0), r, mixed_g); -} - -inline fn opsinAbsorbance(rgb: [3]f32) [3]f32 { - var out: [3]f32 = undefined; - const r = rgb[0]; - const g = rgb[1]; - const b = rgb[2]; - out[0] = @mulAdd(f32, OPSIN_ABSORBANCE_MATRIX[0], r, @mulAdd(f32, OPSIN_ABSORBANCE_MATRIX[1], g, @mulAdd(f32, OPSIN_ABSORBANCE_MATRIX[2], b, OPSIN_ABSORBANCE_BIAS))); - out[1] = @mulAdd(f32, OPSIN_ABSORBANCE_MATRIX[3], r, @mulAdd(f32, OPSIN_ABSORBANCE_MATRIX[4], g, @mulAdd(f32, OPSIN_ABSORBANCE_MATRIX[5], b, OPSIN_ABSORBANCE_BIAS))); - out[2] = @mulAdd(f32, OPSIN_ABSORBANCE_MATRIX[6], r, @mulAdd(f32, OPSIN_ABSORBANCE_MATRIX[7], g, @mulAdd(f32, OPSIN_ABSORBANCE_MATRIX[8], b, OPSIN_ABSORBANCE_BIAS))); - return out; -} - -inline fn mixedToXYB(mixed: [3]f32) [3]f32 { - return .{ - V05 * (mixed[0] - mixed[1]), - V05 * (mixed[0] + mixed[1]), - mixed[2], - }; -} - -inline fn linearRGBtoXYB(input: [3]f32) [3]f32 { - var mixed = opsinAbsorbance(input); - var i: u32 = 0; - while (i < 3) : (i += 1) { - if (mixed[i] < V00) mixed[i] = V00; - mixed[i] = cubeRootAndAdd(mixed[i], ABSORBANCE_BIAS); - } - return mixedToXYB(mixed); -} - -inline fn makePositiveXYB(xyb_in: *[3]f32) void { - // Safely dereference the pointer into a local array, mutate, then write back. - var arr = xyb_in.*; - arr[2] = (arr[2] - arr[1]) + V055; - arr[0] = arr[0] * V140 + V042; - arr[1] = arr[1] + V001; - xyb_in.* = arr; -} - -inline fn xyb(src: [3][]const f32, dst: [3][]f32, idx: usize) void { - const rgb = [3]f32{ - src[0][idx], - src[1][idx], - src[2][idx], - }; - var out = linearRGBtoXYB(rgb); - makePositiveXYB(&out); - inline for (0..3) |i| { - dst[i][idx] = out[i]; - } -} - -inline fn toXYB(srcp: [3][]const f32, dstp: [3][]f32, stride: u32, w: u32, h: u32) void { - const zero: F32Vec = @splat(0.0); - const bias: F32Vec = @splat(OPSIN_ABSORBANCE_BIAS); - var src = srcp; - var dst = dstp; - var y: u32 = 0; - while (y < h) : (y += 1) { - if (w >= VEC_LEN) { - var x: usize = 0; - while (true) { - if (x + VEC_LEN > w) x = w - VEC_LEN; - const r = loadVec(src[0], x); - const g = loadVec(src[1], x); - const b = loadVec(src[2], x); - const ox = opsinVec(r, g, b, K_M00, K_M01, K_M02, bias); - const oy = opsinVec(r, g, b, K_M10, K_M11, K_M12, bias); - const oz = opsinVec(r, g, b, K_M20, K_M21, K_M22, bias); - const cx = cubeRootVec(@max(ox, zero)); - const cy = cubeRootVec(@max(oy, zero)); - const cz = cubeRootVec(@max(oz, zero)); - const xv = @as(F32Vec, @splat(0.5)) * (cx - cy); - const yv = @as(F32Vec, @splat(0.5)) * (cx + cy); - storeVec(dst[0], x, xv * @as(F32Vec, @splat(14.0)) + - @as(F32Vec, @splat(0.42))); - storeVec(dst[1], x, yv + @as(F32Vec, @splat(0.01))); - storeVec(dst[2], x, cz - yv + - @as(F32Vec, @splat(0.55))); - if (x + VEC_LEN >= w) break; - x += VEC_LEN; - } - } else { - var x: usize = 0; - while (x < w) : (x += 1) xyb(src, dst, x); - } - inline for (0..3) |i| { - src[i] = src[i][stride..]; - dst[i] = dst[i][stride..]; - } - } -} - -inline fn tothe4th(y: f64) f64 { - const x = y * y; - return x * x; -} - -inline fn ssimMap( - sumsquared: []f32, - s12: []f32, - mu1: []f32, - mu2: []f32, - stride: u32, - w: u32, - h: u32, - plane: u32, - one_per_pixels: f64, - plane_avg_ssim: []f64, - error_map: ?[]f32, - scale: u32, -) void { - var sum1 = [2]f64{ 0.0, 0.0 }; - var y: u32 = 0; - while (y < h) : (y += 1) { - const row = y * stride; - var x: usize = 0; - if (error_map == null) { - var sum_d: F64Vec = @splat(0.0); - var sum_d4: F64Vec = @splat(0.0); - while (x + VEC_LEN <= w) : (x += VEC_LEN) { - const m1 = loadVec(mu1, row + x); - const m2 = loadVec(mu2, row + x); - const sq = loadVec(sumsquared, row + x); - const s12v = loadVec(s12, row + x); - const m11 = m1 * m1; - const m22 = m2 * m2; - const m12 = m1 * m2; - const m_diff = m1 - m2; - const num_m: F64Vec = @floatCast(@mulAdd( - F32Vec, - m_diff, - -m_diff, - @as(F32Vec, @splat(1.0)), - )); - const num_s: F64Vec = @floatCast(@mulAdd( - F32Vec, - s12v - m12, - @as(F32Vec, @splat(2.0)), - @as(F32Vec, @splat(0.0009)), - )); - const denom_s: F64Vec = @floatCast(sq - - @as(F32Vec, @splat(2.0)) * s12v - m11 - m22 + - @as(F32Vec, @splat(0.0009))); - const d = @max(@as(F64Vec, @splat(1.0)) - - ((num_m * num_s) / denom_s), @as(F64Vec, @splat(0.0))); - const d2 = d * d; - sum_d += d; - sum_d4 += d2 * d2; - } - sum1[0] += @reduce(.Add, sum_d); - sum1[1] += @reduce(.Add, sum_d4); - } - while (x < w) : (x += 1) { - const m1: f32 = @floatCast(mu1[row + x]); - const m2: f32 = @floatCast(mu2[row + x]); - const m11 = m1 * m1; - const m22 = m2 * m2; - const m12 = m1 * m2; - const m_diff = m1 - m2; - const num_m: f64 = @mulAdd(f32, m_diff, -m_diff, 1.0); - const num_s: f64 = @mulAdd(f32, (@as(f32, @floatCast(s12[row + x])) - m12), 2.0, 0.0009); - const denom_s: f64 = sumsquared[row + x] - - 2.0 * s12[row + x] - m11 - m22 + 0.0009; - const d1: f64 = @max(1.0 - ((num_m * num_s) / denom_s), 0.0); - sum1[0] += d1; - sum1[1] += tothe4th(d1); - - if (error_map) |emap| { - const d: f32 = @floatCast(d1); - const weight1: f32 = @floatCast(getWeight(@intCast(plane), @intCast(scale), 0, 0)); - const weight2: f32 = @floatCast(getWeight(@intCast(plane), @intCast(scale), 0, 1)); - emap[row + x] += weight1 * d; - emap[row + x] += weight2 * d; - } - } - } - plane_avg_ssim[plane * 2] = one_per_pixels * sum1[0]; - plane_avg_ssim[plane * 2 + 1] = @sqrt(@sqrt(one_per_pixels * sum1[1])); -} - -inline fn edgeMap( - im1: []f32, - im2: []f32, - mu1: []f32, - mu2: []f32, - stride: u32, - w: u32, - h: u32, - plane: u32, - one_per_pixels: f64, - plane_avg_edge: []f64, - error_map: ?[]f32, - scale: u32, -) void { - var sum2 = [4]f64{ 0.0, 0.0, 0.0, 0.0 }; - var y: u32 = 0; - while (y < h) : (y += 1) { - const row = y * stride; - var x: usize = 0; - if (error_map == null) { - var sum_artifact: F64Vec = @splat(0.0); - var sum_artifact4: F64Vec = @splat(0.0); - var sum_lost: F64Vec = @splat(0.0); - var sum_lost4: F64Vec = @splat(0.0); - while (x + VEC_LEN <= w) : (x += VEC_LEN) { - const n2: F64Vec = @floatCast(@abs( - loadVec(im2, row + x) - loadVec(mu2, row + x), - )); - const n1: F64Vec = @floatCast(@abs( - loadVec(im1, row + x) - loadVec(mu1, row + x), - )); - const one: F64Vec = @splat(1.0); - const zero: F64Vec = @splat(0.0); - const d = (one + n2) / (one + n1) - one; - const artifact = @max(d, zero); - const detail_lost = @max(-d, zero); - const artifact2 = artifact * artifact; - const lost2 = detail_lost * detail_lost; - sum_artifact += artifact; - sum_artifact4 += artifact2 * artifact2; - sum_lost += detail_lost; - sum_lost4 += lost2 * lost2; - } - sum2[0] += @reduce(.Add, sum_artifact); - sum2[1] += @reduce(.Add, sum_artifact4); - sum2[2] += @reduce(.Add, sum_lost); - sum2[3] += @reduce(.Add, sum_lost4); - } - while (x < w) : (x += 1) { - const d1: f64 = (1.0 + - @as(f64, @abs(im2[row + x] - mu2[row + x]))) / - (1.0 + @as(f64, @abs(im1[row + x] - mu1[row + x]))) - - 1.0; - const artifact = @max(d1, 0.0); - sum2[0] += artifact; - sum2[1] += artifact * artifact * artifact * artifact; - const detail_lost = @max(-d1, 0.0); - sum2[2] += detail_lost; - sum2[3] += detail_lost * detail_lost * detail_lost * detail_lost; - - if (error_map) |emap| { - const weight1: f32 = @floatCast(getWeight( - @intCast(plane), - @intCast(scale), - 1, - 0, - )); - const weight2: f32 = @floatCast(getWeight( - @intCast(plane), - @intCast(scale), - 1, - 1, - )); - const weight3: f32 = @floatCast(getWeight( - @intCast(plane), - @intCast(scale), - 2, - 0, - )); - const weight4: f32 = @floatCast(getWeight( - @intCast(plane), - @intCast(scale), - 2, - 1, - )); - emap[row + x] += weight1 * @as(f32, @floatCast(@abs(artifact))); - emap[row + x] += weight2 * @as(f32, @floatCast(@abs(artifact))); - emap[row + x] += weight3 * - @as(f32, @floatCast(@abs(detail_lost))); - emap[row + x] += weight4 * - @as(f32, @floatCast(@abs(detail_lost))); - } - } - } - plane_avg_edge[plane * 4] = one_per_pixels * sum2[0]; - plane_avg_edge[plane * 4 + 1] = @sqrt(@sqrt(one_per_pixels * sum2[1])); - plane_avg_edge[plane * 4 + 2] = one_per_pixels * sum2[2]; - plane_avg_edge[plane * 4 + 3] = @sqrt(@sqrt(one_per_pixels * sum2[3])); -} - -inline fn score(plane_avg_ssim: [6][6]f64, plane_avg_edge: [6][12]f64) f64 { - var ssim_accum: f64 = 0.0; - var idx: usize = 0; - - for (0..3) |plane| { - for (0..6) |scale| { - for (0..2) |n| { - ssim_accum = @mulAdd(f64, weight[idx], @abs(plane_avg_ssim[scale][plane * 2 + n]), ssim_accum); - idx += 1; - ssim_accum = @mulAdd(f64, weight[idx], @abs(plane_avg_edge[scale][plane * 4 + n]), ssim_accum); - idx += 1; - ssim_accum = @mulAdd(f64, weight[idx], @abs(plane_avg_edge[scale][plane * 4 + n + 2]), ssim_accum); - idx += 1; - } - } - } - - ssim_accum *= 0.9562382616834844; - ssim_accum = (6.248496625763138e-5 * ssim_accum * ssim_accum) * ssim_accum + - 2.326765642916932 * ssim_accum - - 0.020884521182843837 * ssim_accum * ssim_accum; - - if (ssim_accum > 0.0) { - ssim_accum = math.pow(f64, ssim_accum, 0.6276336467831387) * -10.0 + 100.0; - } else { - ssim_accum = 100.0; - } - return ssim_accum; -} - -inline fn downscale(src: [3][]f32, dst: [3][]f32, src_stride: u32, in_w: u32, in_h: u32) void { - const fscale: f32 = 2.0; - const uscale: u32 = 2; - const out_w = @divTrunc((in_w + uscale - 1), uscale); - const out_h = @divTrunc((in_h + uscale - 1), uscale); - const dst_stride = @divTrunc((src_stride + uscale - 1), uscale); - const normalize: f32 = 1.0 / (fscale * fscale); - const norm: F32Vec = @splat(normalize); - const full_w = in_w / 2; - - var plane: u32 = 0; - while (plane < 3) : (plane += 1) { - const srcp = src[plane]; - var dstp = dst[plane]; - var oy: u32 = 0; - while (oy < out_h) : (oy += 1) { - var ox: u32 = 0; - if (oy * 2 + 1 < in_h) { - const row0 = srcp[oy * 2 * src_stride ..]; - const row1 = srcp[(oy * 2 + 1) * src_stride ..]; - while (ox + VEC_LEN <= full_w) : (ox += VEC_LEN) { - const a0 = loadVec(row0, ox * 2); - const b0 = loadVec(row0, ox * 2 + VEC_LEN); - const a1 = loadVec(row1, ox * 2); - const b1 = loadVec(row1, ox * 2 + VEC_LEN); - const even_mask = comptime deinterleaveMask(0); - const odd_mask = comptime deinterleaveMask(1); - const e0 = @shuffle(f32, a0, b0, even_mask); - const o0 = @shuffle(f32, a0, b0, odd_mask); - const e1 = @shuffle(f32, a1, b1, even_mask); - const o1 = @shuffle(f32, a1, b1, odd_mask); - storeVec(dstp, ox, ((e0 + o0) + e1 + o1) * norm); - } - } - while (ox < out_w) : (ox += 1) { - var sum: f32 = 0.0; - var iy: u32 = 0; - while (iy < uscale) : (iy += 1) { - var ix: u32 = 0; - while (ix < uscale) : (ix += 1) { - const x: u32 = @min((ox * uscale + ix), (in_w - 1)); - const y: u32 = @min((oy * uscale + iy), (in_h - 1)); - sum += srcp[y * src_stride + x]; - } - } - dstp[ox] = sum * normalize; - } - dstp = dstp[dst_stride..]; - } - } -} - -fn process( - allocator: std.mem.Allocator, - srcp1: [3][]const f16, - srcp2: [3][]const f16, - stride: u32, - w: u32, - h: u32, - error_map: ?[]u32, -) f64 { - const wh: u32 = stride * h; - const error_map_count: u32 = if (error_map != null) 2 else 0; // error_accum and error_scale - const temp_alloc = allocator.alignedAlloc( - f32, - .of(f32), - wh * (16 + error_map_count), - ) catch unreachable; - defer allocator.free(temp_alloc); - const scratch = allocator.alignedAlloc( - f32, - .of(f32), - stride + 64, - ) catch unreachable; - defer allocator.free(scratch); - - return processWithScratch(srcp1, srcp2, stride, w, h, error_map, temp_alloc, scratch); -} - -pub fn processWithScratch( - srcp1: anytype, - srcp2: anytype, - stride: u32, - w: u32, - h: u32, - error_map: ?[]u32, - temp: []f32, - scratch: []f32, -) f64 { - const wh: u32 = stride * h; - const error_map_count: u32 = if (error_map != null) 2 else 0; // error_accum and error_scale - const required: usize = @as(usize, wh) * - @as(usize, 16 + error_map_count); - std.debug.assert(temp.len >= required); - std.debug.assert(scratch.len >= stride + 64); - var temp_use = temp[0..required]; - - var temp4x3: [4][3][]f32 = undefined; - var x: u32 = 0; - for (0..4) |i| { - for (0..3) |ii| { - temp4x3[i][ii] = temp_use[x..(x + wh)]; - x += wh; - } - } - - const tmpp3 = temp_use[x..(x + wh)]; - x += wh; - const tmpsq = temp_use[x..(x + wh)]; - x += wh; - const tmpps12 = temp_use[x..(x + wh)]; - x += wh; - const tmppmu1 = temp_use[x..(x + wh)]; - x += wh; - - var error_accum: []f32 = undefined; - var error_scale: []f32 = undefined; - if (error_map != null) { - error_accum = temp_use[x..(x + wh)]; - x += wh; - error_scale = temp_use[x..(x + wh)]; - x += wh; - @memset(error_scale, 0.0); - } - - const srcp1b = temp4x3[0]; - const srcp2b = temp4x3[1]; - const tmpp1 = temp4x3[2]; - const tmpp2 = temp4x3[3]; - - inline for (0..3) |i| { - for (srcp1[i], 0..) |v, j| - srcp1b[i][j] = @floatCast(v); - for (srcp2[i], 0..) |v, j| - srcp2b[i][j] = @floatCast(v); - } - - var plane_avg_ssim: [6][6]f64 = std.mem.zeroes([6][6]f64); - var plane_avg_edge: [6][12]f64 = std.mem.zeroes([6][12]f64); - var stride2 = stride; - var w2 = w; - var h2 = h; - - var scale: u32 = 0; - while (scale < 6) : (scale += 1) { - if (w2 < 8 or h2 < 8) break; - - if (scale > 0) { - downscale(srcp1b, srcp1b, stride2, w2, h2); - downscale(srcp2b, srcp2b, stride2, w2, h2); - stride2 = @divTrunc((stride2 + 1), 2); - w2 = @divTrunc((w2 + 1), 2); - h2 = @divTrunc((h2 + 1), 2); - if (w2 < 8 or h2 < 8) break; - } - - if (error_map != null) { - const current_wh = stride2 * h2; - @memset(error_scale[0..current_wh], 0.0); - } - - const one_per_pixels: f64 = 1.0 / @as(f64, @floatFromInt(w2 * h2)); - toXYB(srcp1b, tmpp1, stride2, w2, h2); - toXYB(srcp2b, tmpp2, stride2, w2, h2); - - for (0..3) |plane| { - const skip = skip_table[plane][scale]; - if (skip.all()) continue; - const ssim_on = !skip.ssim; - const edge_on = !skip.artifact or !skip.detailloss; - - if (ssim_on) { - multiply(tmpp1[plane], tmpp2[plane], tmpp3, stride2, w2, h2); - blur(tmpp3, tmpps12, stride2, w2, h2, scratch); - addSquare(tmpp1[plane], tmpp2[plane], tmpp3, stride2, w2, h2); - blur(tmpp3, tmpsq, stride2, w2, h2, scratch); - } - - if (ssim_on or edge_on) { - blur(tmpp1[plane], tmppmu1, stride2, w2, h2, scratch); - blur(tmpp2[plane], tmpp3, stride2, w2, h2, scratch); - } - - if (ssim_on) - ssimMap( - tmpsq, - tmpps12, - tmppmu1, - tmpp3, - stride2, - w2, - h2, - @intCast(plane), - one_per_pixels, - &plane_avg_ssim[scale], - if (error_map != null) error_scale else null, - scale, - ); - - if (edge_on) - edgeMap( - tmpp1[plane], - tmpp2[plane], - tmppmu1, - tmpp3, - stride2, - w2, - h2, - @intCast(plane), - one_per_pixels, - &plane_avg_edge[scale], - if (error_map != null) error_scale else null, - @intCast(scale), - ); - } - - if (error_map != null) { - if (scale == 0) - @memcpy(error_accum, error_scale) - else if (scale < 3) - err_map.upscaleAndAccumulate(error_scale, error_accum, stride2, w2, h2, stride, w, h); - } - } - - if (error_map) |emap| - err_map.generateErrorMap(error_accum, emap, stride, w, h); - - return score(plane_avg_ssim, plane_avg_edge); -} diff --git a/src/zig_to_c.zig b/src/zig_to_c.zig index 92dc491..3778692 100644 --- a/src/zig_to_c.zig +++ b/src/zig_to_c.zig @@ -12,181 +12,9 @@ // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // See the License for the specific language governing permissions and // limitations under the License. -const std = @import("std"); const build_opts = @import("build_opts"); -const ssimu2 = @import("ssimulacra2/ssimulacra2.zig"); -const c = @cImport({ - @cInclude("src/fmetrics.h"); -}); - -extern fn fmetrics_workspace_reserve_data( - workspace: ?*c.FmetricsWorkspace, - size: usize, -) ?*anyopaque; export fn fmetrics_version_str() callconv(.c) [*:0]const u8 { const version = build_opts.version ++ "\x00"; return version.ptr; } - -fn validate( - reference: ?*const c.FmetricsImg, - distorted: ?*const c.FmetricsImg, -) c.FmetricsErr { - const ref = reference orelse return c.FMETRICS_ERR_INVALID_ARGUMENT; - const dist = distorted orelse return c.FMETRICS_ERR_INVALID_ARGUMENT; - if (ref.data == null or dist.data == null) - return c.FMETRICS_ERR_INVALID_ARGUMENT; - for ([_]*const c.FmetricsImg{ ref, dist }) |img| { - if (!((img.format == c.FMETRICS_PIX_FMT_RGB_UINT8 and - img.colorspace == c.FMETRICS_COLORSPACE_SRGB) or - (img.format == c.FMETRICS_PIX_FMT_RGB_FLOAT and - img.colorspace == c.FMETRICS_COLORSPACE_LINEAR_SRGB))) - return c.FMETRICS_ERR_UNSUPPORTED_FORMAT; - const row_bytes = @as(usize, img.width) * - @as(usize, if (img.format == c.FMETRICS_PIX_FMT_RGB_FLOAT) 12 else 3); - if (img.stride < row_bytes) return c.FMETRICS_ERR_INVALID_ARGUMENT; - } - if (ref.width != dist.width or ref.height != dist.height) - return c.FMETRICS_ERR_DIMENSION_MISMATCH; - if (ref.width == 0 or ref.height == 0) - return c.FMETRICS_ERR_INVALID_ARGUMENT; - return c.FMETRICS_OK; -} - -fn packedImage(img: *const c.FmetricsImg, buffer: []u8) []const u8 { - const row_bytes = @as(usize, img.width) * 3; - const pixels = @as(usize, img.width) * @as(usize, img.height); - const expected_len = pixels * 3; - const data: [*]const u8 = @ptrCast(img.data.?); - const src = data[0..(@as(usize, img.stride) * img.height)]; - if (img.stride == row_bytes) return src[0..expected_len]; - - if (img.stride == row_bytes) return src[0..expected_len]; - var y: usize = 0; - while (y < img.height) : (y += 1) { - const src_off = y * @as(usize, img.stride); - const dst_off = y * row_bytes; - @memcpy(buffer[dst_off..][0..row_bytes], src[src_off..][0..row_bytes]); - } - return buffer[0..expected_len]; -} - -fn toFmetricsErr(err: ssimu2.Ssimu2Error) c.FmetricsErr { - return switch (err) { - error.InvalidChannelCount => c.FMETRICS_ERR_INVALID_ARGUMENT, - error.OutOfMemory => c.FMETRICS_ERR_OUT_OF_MEMORY, - }; -} - -fn compute( - workspace: ?*c.FmetricsWorkspace, - reference: ?*const c.FmetricsImg, - distorted: ?*const c.FmetricsImg, - result: ?*f64, - error_map: ?[*]u32, -) c.FmetricsErr { - const out = result orelse return c.FMETRICS_ERR_INVALID_ARGUMENT; - const ws = workspace orelse return c.FMETRICS_ERR_INVALID_ARGUMENT; - const valid = validate(reference, distorted); - if (valid != c.FMETRICS_OK) return valid; - - const ref = reference.?; - const dist = distorted.?; - const pixels = @as(usize, ref.width) * @as(usize, ref.height); - const packed_len = std.math.mul(usize, pixels, 3) catch - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const planes_len = std.math.mul(usize, pixels, 6) catch - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const temp_len = std.math.mul(usize, pixels, 18) catch - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const scratch_len = std.math.add(usize, ref.width, 64) catch - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const planes_bytes = std.math.mul(usize, planes_len, @sizeOf(f16)) catch - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const temp_off = std.mem.alignForward(usize, planes_bytes, @alignOf(f32)); - const temp_bytes = std.math.mul(usize, temp_len, @sizeOf(f32)) catch - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const scratch_off = std.math.add(usize, temp_off, temp_bytes) catch - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const scratch_bytes = std.math.mul(usize, scratch_len, @sizeOf(f32)) catch - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const ref_off = std.math.add(usize, scratch_off, scratch_bytes) catch - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const dist_off = std.math.add(usize, ref_off, packed_len) catch - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const total = std.math.add(usize, dist_off, packed_len) catch - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const data = fmetrics_workspace_reserve_data(ws, total) orelse - return c.FMETRICS_ERR_OUT_OF_MEMORY; - const bytes: [*]u8 = @ptrCast(data); - const planes = @as([*]f16, @ptrCast(@alignCast(bytes)))[0..planes_len]; - const temp = @as([*]f32, @ptrCast(@alignCast(bytes + temp_off)))[0..temp_len]; - const scratch = @as([*]f32, @ptrCast(@alignCast(bytes + scratch_off)))[0..scratch_len]; - if (ref.format == c.FMETRICS_PIX_FMT_RGB_FLOAT or - dist.format == c.FMETRICS_PIX_FMT_RGB_FLOAT) - { - const float_planes = temp[0 .. pixels * 6]; - for ([_]*const c.FmetricsImg{ ref, dist }, 0..) |img, image_index| { - const src: [*]const u8 = @ptrCast(img.data.?); - for (0..img.height) |y| { - for (0..img.width) |x| { - for (0..3) |ch| { - const value = if (img.format == c.FMETRICS_PIX_FMT_RGB_FLOAT) blk: { - const offset = y * img.stride + (x * 3 + ch) * 4; - const v: f32 = @bitCast(std.mem.readInt(u32, src[offset..][0..4], @import("builtin").cpu.arch.endian())); - if (!std.math.isFinite(v)) return c.FMETRICS_ERR_INVALID_ARGUMENT; - break :blk v; - } else blk: { - const v: f32 = @as(f32, @floatFromInt(src[y * img.stride + x * 3 + ch])) / 255.0; - break :blk if (v <= 0.04045) v / 12.92 else std.math.pow(f32, (v + 0.055) / 1.055, 2.4); - }; - float_planes[(image_index * 3 + ch) * pixels + y * img.width + x] = value; - } - } - } - } - const rp: [3][]const f32 = .{ float_planes[0..pixels], float_planes[pixels .. pixels * 2], float_planes[pixels * 2 .. pixels * 3] }; - const dp: [3][]const f32 = .{ float_planes[pixels * 3 .. pixels * 4], float_planes[pixels * 4 .. pixels * 5], float_planes[pixels * 5 .. pixels * 6] }; - out.* = ssimu2.processWithScratch(rp, dp, ref.width, ref.width, ref.height, if (error_map) |map| map[0..pixels] else null, temp, scratch); - return c.FMETRICS_OK; - } - const ref_slice = packedImage(ref, (bytes + ref_off)[0..packed_len]); - const dist_slice = packedImage(dist, (bytes + dist_off)[0..packed_len]); - const map_slice = if (error_map) |map| map[0..pixels] else null; - - const score = ssimu2.computeSsimu2WithBuffers( - ref_slice, - dist_slice, - ref.width, - ref.height, - 3, - map_slice, - planes, - temp, - scratch, - ) catch |err| return toFmetricsErr(err); - - out.* = score; - return c.FMETRICS_OK; -} - -export fn fmetrics_ssimu2_cmp( - workspace: ?*c.FmetricsWorkspace, - reference: ?*const c.FmetricsImg, - distorted: ?*const c.FmetricsImg, - result: ?*f64, -) callconv(.c) c.FmetricsErr { - return compute(workspace, reference, distorted, result, null); -} - -export fn fmetrics_ssimu2_cmp_map( - workspace: ?*c.FmetricsWorkspace, - reference: ?*const c.FmetricsImg, - distorted: ?*const c.FmetricsImg, - result: ?*f64, - error_map: ?[*]u32, -) callconv(.c) c.FmetricsErr { - if (error_map == null) return c.FMETRICS_ERR_INVALID_ARGUMENT; - return compute(workspace, reference, distorted, result, error_map); -} diff --git a/tests/bench/README.md b/tests/bench/README.md new file mode 100644 index 0000000..7598ff5 --- /dev/null +++ b/tests/bench/README.md @@ -0,0 +1,142 @@ +# SSIMULACRA2 validation + +The production implementation is C with compiler vectorization. + +```sh +zig build --release=fast +``` + +Each active scale consumes linear RGB rows, sends 2x2 averages to the next +scale, converts only the required XYB channels, and retains nine filtered +rows. No full-frame linear RGB or XYB planes are allocated. Score-only scratch +space is bounded by `1024 * width + 4096` bytes, independent of height. Error +maps additionally retain error planes for the first three scales. + +SSIM uses the blurred squared pixel difference to compute the variance of the +difference directly. This avoids subtracting nearly equal second moments; +identical inputs score exactly 100. Fourth-power norms use float row reductions +and double accumulation across rows. sRGB input is decoded through a float LUT +rather than quantized to half precision. + +## Upstream comparison + +`ssimulacra2.c` compares the installed library with the unmodified +[vapoursynth-zip 23.0.0 filter](https://github.com/dnjulek/vapoursynth-zip/blob/23.0.0/src/filters/ssimulacra2.zig). +Both consume the same deterministic linear RGB samples: interleaved for the +fmetrics API and planar for upstream. Inputs are shared read-only across +workers; each worker has its own reusable scratch allocation. The timed region +excludes input generation, allocation, and three warm-up frames. fmetrics input +validation remains included. This compares filter implementations directly, +without VapourSynth scheduling, decoding, or CLI I/O. + +To reproduce on macOS, from the repository root: + +```sh +zig build --release=fast +work=$(mktemp -d) +mkdir "$work/filters" +for name in ssimulacra2 simd srgb; do + gh api "repos/dnjulek/vapoursynth-zip/contents/src/filters/$name.zig?ref=23.0.0" \ + -H 'Accept: application/vnd.github.raw+json' > "$work/filters/$name.zig" +done +gh api 'repos/dnjulek/vapoursynth-zip/contents/src/vcl.zig?ref=23.0.0' \ + -H 'Accept: application/vnd.github.raw+json' > "$work/vcl.zig" +cp tests/bench/vszip.zig "$work/bench.zig" +zig build-lib -dynamic -OReleaseFast -lc "$work/bench.zig" \ + -femit-bin="$work/libvszip.dylib" +cc -O3 -ffast-math -I . tests/bench/ssimulacra2.c \ + zig-out/lib/libfmetrics.a "$work/libvszip.dylib" -lpthread -lm \ + -o "$work/bench" +"$work/bench" 1920 1080 1 60 +"$work/bench" 1920 1080 8 20 +"$work/bench" 1920 1080 16 12 +"$work/bench" 3840 2160 1 15 +``` + +On Linux, change `.dylib` to `.so`. Output columns are implementation, width, +height, workers, frames/second, scratch bytes/worker, and score. + +Measured on an 18-core Apple M5 Pro, native Zig 0.17 ReleaseFast, macOS arm64: + +| Size | Workers | fmetrics fps | upstream fps | Speedup | +| --- | ---: | ---: | ---: | ---: | +| 1920x1080 | 1 | 47.07 | 40.37 | 1.17x | +| 1920x1080 | 8 | 322.86 | 273.37 | 1.18x | +| 1920x1080 | 16 | 555.80 | 463.74 | 1.20x | +| 3840x2160 | 1 | 11.75 | 9.97 | 1.18x | + +Scratch at 1080p is 1,970,176 bytes versus 13,909,760 upstream; at 4K it is +3,936,256 bytes versus 52,655,360. These results do not establish x86 throughput. + +The existing 3840x2160 PAM CLI pair takes 92.1 ms versus 173.1 ms before the +change (10 runs, two warm-ups); process peak memory falls from 714.3 MiB to +116.0 MiB. Scores are 70.5839 versus 70.5898. + +## Subjective correlation + +All 4,341 pairs in `../mos-correlation/assets/CID22_validation_set.csv` were +scored with the before/after CLIs against the same decoded PNGs. Correlations +use the CSV's `MCOS` values and SciPy's Kendall, Spearman, and Pearson functions, +matching `../mos-correlation/mos.py`. Every pair produced a score, including +negative scores. The baseline was rebuilt from `e8bf3cf` with only the build +and C-import adaptations required by Zig 0.17. + +| Implementation | KRCC | SRCC | PCC | +| --- | ---: | ---: | ---: | +| Previous fmetrics | 0.708176 | 0.891111 | 0.873170 | +| Streaming C | 0.708025 | 0.890985 | 0.873013 | + +Reproduce with SciPy installed, saving every prediction as well as correlations: + +```sh +python3 tests/bench/correlation.py zig-out/bin/fmetrics \ + --output /tmp/fmetrics-cid22.json +``` + +Run the same command with the baseline executable to compare paired predictions. + +Each coefficient changes by less than 0.0002. These measurements cover CID22; +correlation on other subjective datasets has not been measured. + +## Video clip comparison + +The executable script uses uv inline metadata to install its Python +dependencies, including vszip 23.0.0 and fmetrics built from this checkout: + +```sh +./tests/bench/compare_video.py reference.mkv distorted.mkv \ + --threads 16 --runs 3 --csv scores.csv +``` + +Install uv, the VapourSynth system runtime, and a source plugin separately. +FFMS2 is the default decoder; its VapourSynth plugin must be installed. +`--source bestsource` uses `core.bs.VideoSource` instead. Plugins outside the +autoload path can be supplied with repeated `--load-plugin /path/to/plugin` +arguments. The script also finds the plugin in an installed `vapoursynth-vszip` +wheel when autoloading is unavailable. It prints the loaded plugin version and +path. + +Both metrics receive identical float linear RGB frames. YUV is converted with +VapourSynth's Bicubic resizer, using frame color metadata and `--matrix 709` as +the matrix fallback. RGB samples are clamped to [0, 1] before linearization. +This comparison targets SDR video. `--range limited` or `--range full` overrides +the input range; the default uses source metadata. No frame resizing occurs. + +The CSV records both scores and their signed delta for every frame. The summary +includes means, mean absolute delta, largest absolute delta, and median +throughput. Timing includes decoding, shared RGB preprocessing, and the Python +packing needed by the fmetrics bindings. It is a video-pipeline comparison; +use the C harness above for direct filter throughput. Metric output caches are +disabled, each pass warms up three frames, and pass order alternates between +runs. Decoder and operating-system caches can still affect timing. + +To sample 100 frames starting at frame 50, every fifth frame: + +```sh +./tests/bench/compare_video.py reference.mkv distorted.mkv \ + --start 50 --frames 100 --step 5 --csv sample.csv +``` + +`--frames` counts sampled frames. `--dist-start` selects a different starting +frame in the distorted clip when alignment requires it. Dimensions and frame +rates must match. Different lengths require an explicit common `--frames` span. diff --git a/tests/bench/compare_video.py b/tests/bench/compare_video.py new file mode 100755 index 0000000..b8f13ad --- /dev/null +++ b/tests/bench/compare_video.py @@ -0,0 +1,199 @@ +#!/usr/bin/env -S uv run --script +# /// script +# requires-python = ">=3.14" +# dependencies = [ +# "image-fmetrics", +# "numpy", +# "vapoursynth>=77", +# "vapoursynth-vszip==23.0.0", +# ] +# [tool.uv.sources] +# image-fmetrics = { path = "../.." } +# /// + +import argparse +import csv +from importlib import metadata +import math +import os +import queue +from pathlib import Path +import statistics +import sys +import time + + +def arguments(): + parser = argparse.ArgumentParser( + description="Compare fmetrics and vszip SSIMULACRA2 on SDR video.", + epilog="Requires image-fmetrics, numpy, VapourSynth, vszip, and a source plugin.", + ) + parser.add_argument("reference", type=Path) + parser.add_argument("distorted", type=Path) + parser.add_argument("--csv", type=Path, default=Path("ssimu2-comparison.csv")) + parser.add_argument("--source", choices=["ffms2", "bestsource"], default="ffms2") + parser.add_argument("--load-plugin", action="append", default=[], type=Path) + parser.add_argument("--threads", type=int, default=os.cpu_count() or 1) + parser.add_argument("--runs", type=int, default=3) + parser.add_argument("--start", type=int, default=0) + parser.add_argument("--dist-start", type=int) + parser.add_argument("--frames", type=int) + parser.add_argument("--step", type=int, default=1) + parser.add_argument("--matrix", default="709", help="YUV matrix fallback (default: 709)") + parser.add_argument("--range", choices=["auto", "limited", "full"], default="auto") + args = parser.parse_args() + if min(args.threads, args.runs, args.step) < 1: + parser.error("--threads, --runs, and --step must be positive") + if args.start < 0 or (args.dist_start is not None and args.dist_start < 0): + parser.error("frame offsets cannot be negative") + if args.frames is not None and args.frames < 1: + parser.error("--frames must be positive") + for path in [args.reference, args.distorted, *args.load_plugin]: + if not path.is_file(): + parser.error(f"file does not exist: {path}") + if args.csv.resolve() in {args.reference.resolve(), args.distorted.resolve()}: + parser.error("CSV output must differ from the input files") + return args + + +def main(): + args = arguments() + try: + import fmetrics + import numpy as np + import vapoursynth as vs + except ImportError as error: + raise RuntimeError( + f"{error}. Run this script with uv and install the VapourSynth " + "runtime plus a source plugin." + ) from error + + core = vs.core + core.num_threads = args.threads + for path in args.load_plugin: + core.std.LoadPlugin(str(path.resolve())) + if not hasattr(core, "vszip"): + try: + package = metadata.distribution("vapoursynth-vszip") + for path in package.files or []: + if path.name in {"libvszip.dylib", "libvszip.so", "vszip.dll"}: + core.std.LoadPlugin(str(package.locate_file(path))) + break + except metadata.PackageNotFoundError: + pass + if not hasattr(core, "vszip"): + raise RuntimeError("vszip is missing; install it or pass --load-plugin PATH") + namespace = "ffms2" if args.source == "ffms2" else "bs" + if not hasattr(core, namespace): + raise RuntimeError(f"{namespace} is missing; install it or pass --load-plugin PATH") + + def source(path): + if args.source == "ffms2": + clip = core.ffms2.Source(str(path.resolve()), cache=False) + else: + clip = core.bs.VideoSource(str(path.resolve())) + core.std.SetVideoCache(clip, mode=0) + return clip + + ref, dist = source(args.reference), source(args.distorted) + if not ref.format or not dist.format: + raise ValueError("variable-format video is not supported") + if (ref.width, ref.height) != (dist.width, dist.height): + raise ValueError("clips must have the same dimensions") + if ref.fps_num * dist.fps_den != dist.fps_num * ref.fps_den: + raise ValueError("clips must have the same frame rate") + dist_start = args.start if args.dist_start is None else args.dist_start + lengths = [ref.num_frames - args.start, dist.num_frames - dist_start] + if args.frames is None and lengths[0] != lengths[1]: + raise ValueError("clip lengths differ; select a common span with --frames") + span = min(lengths) if args.frames is None else (args.frames - 1) * args.step + 1 + if span < 1 or span > min(lengths): + raise ValueError("requested frame selection exceeds the clips") + ref = ref[args.start:args.start + span:args.step] + dist = dist[dist_start:dist_start + span:args.step] + + def linear_rgb(clip): + if args.range != "auto": + clip = core.std.SetFrameProps(clip, _ColorRange=int(args.range == "limited")) + options = {"matrix_in_s": args.matrix} if clip.format.color_family == vs.YUV else {} + encoded = core.resize.Bicubic(clip, format=vs.RGBS, **options) + rgb = core.std.Expr(encoded, expr="x 0 max 1 min") + linear = core.resize.Bicubic(rgb, transfer_in_s="srgb", transfer_s="linear") + for node in [clip, encoded, rgb, linear]: + core.std.SetVideoCache(node, mode=0) + return linear + + ref, dist = linear_rgb(ref), linear_rgb(dist) + workspaces = queue.SimpleQueue() + + def compare(n, f): + try: + workspace, rgb = workspaces.get_nowait() + except queue.Empty: + workspace = fmetrics.Workspace() + rgb = [np.empty((ref.height, ref.width, 3), np.float32) for _ in range(2)] + try: + for image, frame in zip(rgb, f): + for channel in range(3): + image[:, :, channel] = np.asarray(frame[channel]) + score = workspace.ssimu2(*rgb) + out = f[0].copy() + out.props["FMETRICS_SSIMU2"] = score + return out + finally: + workspaces.put((workspace, rgb)) + + clips = { + "fmetrics": (core.std.ModifyFrame(ref, clips=[ref, dist], selector=compare), + "FMETRICS_SSIMU2"), + "vszip": (core.vszip.SSIMULACRA2(ref, dist), "SSIMULACRA2"), + } + for clip, _ in clips.values(): + core.std.SetVideoCache(clip, mode=0) + print(f"fmetrics {fmetrics.__version__}; vszip {core.vszip.version}") + print(f"vszip plugin: {core.vszip.plugin_path}") + print(f"{ref.width}x{ref.height}, {ref.num_frames} frames, {args.threads} threads") + print("Shared linear RGB; timing includes decoding, conversion, and Python packing.") + elapsed = {name: [] for name in clips} + scores = {} + for run in range(args.runs): + order = list(clips) if run % 2 == 0 else list(reversed(clips)) + for name in order: + clip, prop = clips[name] + for frame in clip[:min(3, ref.num_frames)].frames(close=True): + pass + start = time.perf_counter() + values = [float(frame.props[prop]) for frame in + clip.frames(prefetch=args.threads, close=True)] + duration = time.perf_counter() - start + if not all(map(math.isfinite, values)): + raise ValueError(f"{name} returned a non-finite score") + if name in scores and values != scores[name]: + raise RuntimeError(f"{name} scores changed between runs") + scores[name] = values + elapsed[name].append(duration) + print(f"Run {run + 1}: {name:8} {duration:.3f} s, " + f"{len(values) / duration:.2f} fps, mean {statistics.mean(values):.6f}", + flush=True) + + deltas = [a - b for a, b in zip(scores["fmetrics"], scores["vszip"])] + with args.csv.open("w", newline="") as stream: + writer = csv.writer(stream) + writer.writerow(["reference_frame", "distorted_frame", "fmetrics", "vszip", "delta"]) + for i, (a, b, delta) in enumerate(zip(scores["fmetrics"], scores["vszip"], deltas)): + writer.writerow([args.start + i * args.step, dist_start + i * args.step, a, b, delta]) + print(f"Delta (fmetrics - vszip): mean {statistics.mean(deltas):+.6f}, " + f"MAE {statistics.mean(map(abs, deltas)):.6f}, max abs {max(map(abs, deltas)):.6f}") + for name in clips: + duration = statistics.median(elapsed[name]) + print(f"{name:8} median: {duration:.3f} s, {ref.num_frames / duration:.2f} fps") + ratio = statistics.median(elapsed["vszip"]) / statistics.median(elapsed["fmetrics"]) + print(f"Throughput ratio fmetrics/vszip: {ratio:.3f}x") + print(f"Per-frame scores: {args.csv.resolve()}") + + +if __name__ == "__main__": + try: + main() + except (RuntimeError, ValueError, OSError) as error: + sys.exit(f"Error: {error}") diff --git a/tests/bench/correlation.py b/tests/bench/correlation.py new file mode 100644 index 0000000..1d5c067 --- /dev/null +++ b/tests/bench/correlation.py @@ -0,0 +1,49 @@ +import argparse +from concurrent.futures import ThreadPoolExecutor +import csv +import json +import math +from pathlib import Path +import re +import subprocess + +from scipy import stats + + +parser = argparse.ArgumentParser() +parser.add_argument("binary", type=Path) +parser.add_argument("--assets", type=Path, + default=Path("../mos-correlation/assets")) +parser.add_argument("--workers", type=int, default=6) +parser.add_argument("--output", type=Path, required=True) +args = parser.parse_args() +rows = list(csv.DictReader( + (args.assets / "CID22_validation_set.csv").open() +)) + + +def score(row): + paths = [args.assets / "CID22_validation_set" / row[key] + for key in ("reference_img", "distorted_img")] + result = subprocess.run( + [str(args.binary.resolve()), "ssimu2", *map(str, paths)], + capture_output=True, text=True, check=True, + ) + values = re.findall(r"^\s*(-?\d+(?:\.\d+)?)\s*$", + result.stdout + "\n" + result.stderr, re.MULTILINE) + if len(values) != 1 or not math.isfinite(float(values[0])): + raise ValueError(f"Invalid score for {paths}: {result}") + return float(values[0]) + + +with ThreadPoolExecutor(max_workers=args.workers) as pool: + scores = list(pool.map(score, rows)) +truth = [float(row["MCOS"]) for row in rows] +result = { + "pairs": len(rows), + "KRCC": float(stats.kendalltau(truth, scores).statistic), + "SRCC": float(stats.spearmanr(truth, scores).statistic), + "PCC": float(stats.pearsonr(truth, scores).statistic), +} +print(json.dumps(result, indent=2)) +args.output.write_text(json.dumps({**result, "scores": scores}, indent=2)) diff --git a/tests/bench/ssimulacra2.c b/tests/bench/ssimulacra2.c new file mode 100644 index 0000000..2cabf4f --- /dev/null +++ b/tests/bench/ssimulacra2.c @@ -0,0 +1,101 @@ +#include "src/fmetrics.h" +#include "src/common/mem.h" +#include +#include +#include +#include + +extern size_t upstream_size(unsigned, unsigned); +extern double upstream_cmp(const float *, const float *, unsigned, unsigned, + float *); + +static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER; +static pthread_cond_t cond = PTHREAD_COND_INITIALIZER; +static int ready, go, mode, width, height, iterations; +static float *planar[2]; +static FmetricsImg images[2]; + +typedef struct Result { + double score; + size_t bytes; +} Result; + +static double seconds(void) { + struct timespec ts; + clock_gettime(CLOCK_MONOTONIC, &ts); + return ts.tv_sec + ts.tv_nsec * 1e-9; +} + +static void *worker(void *arg) { + Result *r = arg; + FmetricsWorkspace *ws = fmetrics_workspace_create(); + float *scratch = mode ? malloc(upstream_size(width, height) * 4) : NULL; + for (int i = -3; i < iterations; i++) { + if (!i) { + pthread_mutex_lock(&mutex); + ready++; + pthread_cond_broadcast(&cond); + while (!go) pthread_cond_wait(&cond, &mutex); + pthread_mutex_unlock(&mutex); + } + if (mode) r->score = upstream_cmp(planar[0], planar[1], width, + height, scratch); + else if (fmetrics_ssimu2_cmp(ws, &images[0], &images[1], &r->score)) + abort(); + } + r->bytes = mode ? upstream_size(width, height) * 4 : ws->scratch.size; + free(scratch); + fmetrics_workspace_destroy(ws); + return NULL; +} + +int main(int argc, char **argv) { + width = argc > 1 ? atoi(argv[1]) : 1920; + height = argc > 2 ? atoi(argv[2]) : 1080; + const int threads = argc > 3 ? atoi(argv[3]) : 1; + iterations = argc > 4 ? atoi(argv[4]) : 30; + if (width < 8 || height < 8 || threads < 1 || iterations < 1) return 1; + const size_t n = (size_t)width * height; + for (int img = 0; img < 2; img++) { + float *rgb = malloc(n * 12); + planar[img] = malloc(n * 12); + if (!rgb || !planar[img]) return 1; + images[img] = (FmetricsImg){rgb, width, height, width * 12, + FMETRICS_PIX_FMT_RGB_FLOAT, FMETRICS_COLORSPACE_LINEAR_SRGB, false}; + for (int y = 0; y < height; y++) for (int x = 0; x < width; x++) + for (int c = 0; c < 3; c++) { + const size_t i = (size_t)y * width + x; + const float v = 0.4 + 0.2 * sin(x * 0.03 + c) + + 0.15 * cos(y * 0.04 + c) + 0.05 * sin(x + y * 0.8) + + img * 0.01 * sin(x * 0.4 + y * 0.3 + c); + rgb[i * 3 + c] = planar[img][c * n + i] = v; + } + } + pthread_t *handles = malloc((size_t)threads * sizeof(*handles)); + Result *results = calloc((size_t)threads, sizeof(*results)); + if (!handles || !results) return 1; + for (mode = 0; mode < 2; mode++) { + ready = go = 0; + for (int t = 0; t < threads; t++) + if (pthread_create(handles + t, NULL, worker, results + t)) abort(); + pthread_mutex_lock(&mutex); + while (ready != threads) pthread_cond_wait(&cond, &mutex); + const double start = seconds(); + go = 1; + pthread_cond_broadcast(&cond); + pthread_mutex_unlock(&mutex); + for (int t = 0; t < threads; t++) pthread_join(handles[t], NULL); + const double elapsed = seconds() - start; + printf("%s,%d,%d,%d,%.3f,%zu,%.9f\n", + mode ? "vszip23" : "fmetrics", width, height, threads, + threads * iterations / elapsed, results[0].bytes, + results[0].score); + fflush(stdout); + } + for (int img = 0; img < 2; img++) { + free((void *)images[img].data); + free(planar[img]); + } + free(handles); + free(results); +} diff --git a/tests/bench/vszip.zig b/tests/bench/vszip.zig new file mode 100644 index 0000000..e99437e --- /dev/null +++ b/tests/bench/vszip.zig @@ -0,0 +1,10 @@ +const f = @import("filters/ssimulacra2.zig"); +export fn upstream_size(w: u32, h: u32) usize { + return f.scratchSize(w, h); +} +export fn upstream_cmp(a: [*]const f32, b: [*]const f32, w: u32, h: u32, scratch: [*]f32) f64 { + const n: usize = @as(usize, w) * h; + const r = f.Input{ .kind = .f32, .planes = .{ @ptrCast(a), @ptrCast(a + n), @ptrCast(a + 2 * n) }, .stride = @as(usize, w) * 4, .srgb = false }; + const d = f.Input{ .kind = .f32, .planes = .{ @ptrCast(b), @ptrCast(b + n), @ptrCast(b + 2 * n) }, .stride = @as(usize, w) * 4, .srgb = false }; + return f.process(r, d, w, h, scratch[0..f.scratchSize(w, h)]); +}