diff --git a/libvmaf/src/feature/ciede.c b/libvmaf/src/feature/ciede.c index 7bbfa8f60..bd9104f8a 100644 --- a/libvmaf/src/feature/ciede.c +++ b/libvmaf/src/feature/ciede.c @@ -68,9 +68,9 @@ static void scale_chroma_planes_hbd(VmafPicture *in, VmafPicture *out) uint16_t *out_buf = out->data[p]; for (unsigned i = 0; i < out->h[p]; i++) { for (unsigned j = 0; j < out->w[p]; j++) { - out_buf[j] = in_buf[(j / ((p && ss_ver) ? 2 : 1))]; + out_buf[j] = in_buf[(j / ((p && ss_hor) ? 2 : 1))]; } - in_buf += (((p && ss_hor) ? i % 2 : 1) * in->stride[p]) / 2; + in_buf += (((p && ss_ver) ? i % 2 : 1) * in->stride[p]) / 2; out_buf += out->stride[p] / 2; } } @@ -86,9 +86,9 @@ static void scale_chroma_planes(VmafPicture *in, VmafPicture *out) uint8_t *out_buf = out->data[p]; for (unsigned i = 0; i < out->h[p]; i++) { for (unsigned j = 0; j < out->w[p]; j++) { - out_buf[j] = in_buf[(j / ((p && ss_ver) ? 2 : 1))]; + out_buf[j] = in_buf[(j / ((p && ss_hor) ? 2 : 1))]; } - in_buf += ((p && ss_hor) ? i % 2 : 1) * in->stride[p]; + in_buf += ((p && ss_ver) ? i % 2 : 1) * in->stride[p]; out_buf += out->stride[p]; } } diff --git a/libvmaf/src/feature/cuda/ciede/ciede.cu b/libvmaf/src/feature/cuda/ciede/ciede.cu new file mode 100644 index 000000000..6c80f2150 --- /dev/null +++ b/libvmaf/src/feature/cuda/ciede/ciede.cu @@ -0,0 +1,331 @@ +/** + * + * Copyright 2026 Bardie Høgh Joensen + * + * Licensed under the BSD+Patent License (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://opensource.org/licenses/BSDplusPatent + * + * 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. + * + */ + +// CUDA port of feature/ciede.c, which is in large part a port of the +// ciede2000 implementation from av-metrics +// (https://github.com/rust-av/av-metrics) with the following license: + +/* +The MIT License (MIT) +Copyright (c) 2019 Joshua Holmer + +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in +the Software without restriction, including without limitation the rights to +use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies +of the Software, and to permit persons to whom the Software is furnished to do +so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +/* + * Unlike the psnr_cuda/ssim_cuda kernels, this port is NOT bit-exact with the + * CPU extractor: ciede is dominated by libm transcendentals (pow/atan2/sin/ + * cos/exp), which differ between glibc and CUDA in the low bits regardless of + * precision, and double-precision throughput is 1/64 rate on consumer GPUs. + * The device math therefore runs in float32 — the CPU reference truncates + * every intermediate to float anyway — and parity with the CPU extractor is + * validated to a per-frame score tolerance (|delta| <= 1e-3) instead of + * bit-exactness. The per-block partial sums and the final pooling stay in + * double, summed in a deterministic order. + */ + +#include "cuda_helper.cuh" + +#include "common.h" + +#define M_PI_F 3.14159265358979323846f + +typedef struct LABColorf { + float l; + float a; + float b; +} LABColorf; + +__device__ __forceinline__ float get_h_prime(const float x, const float y) +{ + if ((x == 0.0f) && (y == 0.0f)) + return 0.0f; + float hue_angle = atan2f(x, y); + if (hue_angle < 0.0f) + hue_angle += 2.f * M_PI_F; + return hue_angle; +} + +__device__ __forceinline__ float get_delta_h_prime(const float c1, + const float c2, const float h_prime_1, const float h_prime_2) +{ + if ((c1 == 0.0f) || (c2 == 0.0f)) + return 0.0f; + if (fabsf(h_prime_1 - h_prime_2) <= M_PI_F) + return h_prime_2 - h_prime_1; + if (h_prime_2 <= h_prime_1) + return h_prime_2 - h_prime_1 + 2.f * M_PI_F; + else + return h_prime_2 - h_prime_1 - 2.f * M_PI_F; +} + +__device__ __forceinline__ float get_upcase_h_bar_prime(const float h_prime_1, + const float h_prime_2) +{ + return fabsf(h_prime_1 - h_prime_2) > M_PI_F ? + (h_prime_1 + h_prime_2 + 2.0f * M_PI_F) / 2.0f : + (h_prime_1 + h_prime_2) / 2.0f; +} + +__device__ __forceinline__ float get_upcase_t(const float upcase_h_bar_prime) +{ + return 1.0f - + 0.17f * cosf(upcase_h_bar_prime - M_PI_F / 6.0f) + + 0.24f * cosf(2.0f * upcase_h_bar_prime) + + 0.32f * cosf(3.0f * upcase_h_bar_prime + M_PI_F / 30.0f) - + 0.20f * cosf(4.0f * upcase_h_bar_prime - 7.0f * M_PI_F / 20.0f); +} + +__device__ __forceinline__ float get_r_sub_t(const float c_bar_prime, + const float upcase_h_bar_prime) +{ + const float degrees = + (upcase_h_bar_prime * (180.0f / M_PI_F) - 275.0f) * (1.0f / 25.0f); + const float c7 = powf(c_bar_prime, 7); + + return -2.0f * + sqrtf(c7 / (c7 + powf(25.f, 7))) * + sinf((60.0f * expf(-(degrees * degrees))) * (M_PI_F / 180.0f)); +} + +__device__ float ciede2000(LABColorf color_1, LABColorf color_2) +{ + // default ksub from feature/ciede.c: l = 0.65, c = 1.0, h = 4.0 + const float ksub_l = 0.65f, ksub_c = 1.0f, ksub_h = 4.0f; + + const float delta_l_prime = color_2.l - color_1.l; + const float l_bar = (color_1.l + color_2.l) / 2; + const float c1 = sqrtf(color_1.a * color_1.a + color_1.b * color_1.b); + const float c2 = sqrtf(color_2.a * color_2.a + color_2.b * color_2.b); + const float c_bar = (c1 + c2) / 2; + const float c_bar7 = powf(c_bar, 7); + const float chroma_shift = 1 - sqrtf(c_bar7 / (c_bar7 + powf(25.f, 7))); + const float a_prime_1 = color_1.a + (color_1.a / 2) * chroma_shift; + const float a_prime_2 = color_2.a + (color_2.a / 2) * chroma_shift; + const float c_prime_1 = + sqrtf(a_prime_1 * a_prime_1 + color_1.b * color_1.b); + const float c_prime_2 = + sqrtf(a_prime_2 * a_prime_2 + color_2.b * color_2.b); + const float c_bar_prime = (c_prime_1 + c_prime_2) / 2; + const float delta_c_prime = c_prime_2 - c_prime_1; + const float s_sub_l = 1.f + ((0.015f * ((l_bar - 50) * (l_bar - 50))) / + sqrtf(20 + ((l_bar - 50) * (l_bar - 50)))); + const float s_sub_c = 1.f + 0.045f * c_bar_prime; + const float h_prime_1 = get_h_prime(color_1.b, a_prime_1); + const float h_prime_2 = get_h_prime(color_2.b, a_prime_2); + const float delta_h_prime = get_delta_h_prime(c1, c2, h_prime_1, h_prime_2); + const float delta_upcase_h_prime = + 2.0f * sqrtf(c_prime_1 * c_prime_2) * sinf(delta_h_prime / 2.0f); + const float upcase_h_bar_prime = + get_upcase_h_bar_prime(h_prime_1, h_prime_2); + const float upcase_t = get_upcase_t(upcase_h_bar_prime); + const float s_sub_upcase_h = 1.0f + 0.015f * c_bar_prime * upcase_t; + const float r_sub_t = get_r_sub_t(c_bar_prime, upcase_h_bar_prime); + const float lightness = delta_l_prime / (ksub_l * s_sub_l); + const float chroma = delta_c_prime / (ksub_c * s_sub_c); + const float hue = delta_upcase_h_prime / (ksub_h * s_sub_upcase_h); + + return sqrtf(lightness * lightness + chroma * chroma + + hue * hue + r_sub_t * chroma * hue); +} + +__device__ __forceinline__ float rgb_to_xyz_map(float c) +{ + if (c > 10.f / 255.f) { + const float A = 0.055f; + const float D = 1.0f / 1.055f; + return powf((c + A) * D, 2.4f); + } else { + const float D = 1.0f / 12.92f; + return (c * D); + } +} + +__device__ __forceinline__ float xyz_to_lab_map(float c) +{ + const float KAPPA = 24389.0f / 27.0f; + const float EPSILON = 216.0f / 24389.0f; + + if (c > EPSILON) { + return powf(c, 1.0f / 3.0f); + } else { + return (KAPPA * c + 16.0f) * (1.0f / 116.0f); + } +} + +__device__ LABColorf get_lab_color(float y, float u, float v, unsigned bpc) +{ + const float scale = 1 << (bpc - 8); + + y = (y - 16.f * scale) * (1.f / (219.f * scale)); + u = (u - 128.f * scale) * (1.f / (224.f * scale)); + v = (v - 128.f * scale) * (1.f / (224.f * scale)); + + // Assumes BT.709 + float r = y + 1.28033f * v; + float g = y - 0.21482f * u - 0.38059f * v; + float b = y + 2.12798f * u; + + r = rgb_to_xyz_map(r); + g = rgb_to_xyz_map(g); + b = rgb_to_xyz_map(b); + + float x = r * 0.4124564390896921f + g * 0.357576077643909f + + b * 0.18043748326639894f; + y = r * 0.21267285140562248f + g * 0.715152155287818f + + b * 0.07217499330655958f; + float z = r * 0.019333895582329317f + g * 0.119192025881303f + + b * 0.9503040785363677f; + + x = xyz_to_lab_map(x * (1.0f / 0.95047f)); + y = xyz_to_lab_map(y); + z = xyz_to_lab_map(z * (1.0f / 1.08883f)); + + LABColorf lab_color = { + (116.0f * y) - 16.0f, + 500.0f * (x - y), + 200.0f * (y - z), + }; + + return lab_color; +} + +// Block-reduce per-thread ΔE (double partials) and write one partial per +// block; the host sums the partials sequentially so pooling is deterministic +__device__ __forceinline__ void block_reduce_de(double de, + VmafCudaBuffer partials) +{ + __shared__ double sh[256]; + + const int t = threadIdx.y * blockDim.x + threadIdx.x; + sh[t] = de; + __syncthreads(); + for (int stride = 128; stride > 0; stride >>= 1) { + if (t < stride) + sh[t] += sh[t + stride]; + __syncthreads(); + } + if (t == 0) + reinterpret_cast(partials.data)[ + blockIdx.y * gridDim.x + blockIdx.x] = sh[0]; +} + +// scale_chroma_planes in feature/ciede.c halves the column index for +// horizontal subsampling and advances the source row every second output row +// for vertical subsampling. +__device__ __forceinline__ int chroma_col(int j, int ss_hor) +{ + return ss_hor ? j / 2 : j; +} + +__device__ __forceinline__ int chroma_row(int i, int ss_ver) +{ + return ss_ver ? i / 2 : i; +} + +extern "C" { + +__global__ void ciede_kernel_8bpc(const VmafPicture ref, const VmafPicture dis, + VmafCudaBuffer partials, unsigned width, unsigned height, + int ss_hor, int ss_ver) +{ + const int j = blockIdx.x * blockDim.x + threadIdx.x; + const int i = blockIdx.y * blockDim.y + threadIdx.y; + + double de00 = 0.0; + if (j < width && i < height) { + const int cj = chroma_col(j, ss_hor); + const int ci = chroma_row(i, ss_ver); + + const float r_y = (reinterpret_cast(ref.data[0]) + + (size_t)i * ref.stride[0])[j]; + const float r_u = (reinterpret_cast(ref.data[1]) + + (size_t)ci * ref.stride[1])[cj]; + const float r_v = (reinterpret_cast(ref.data[2]) + + (size_t)ci * ref.stride[2])[cj]; + const float d_y = (reinterpret_cast(dis.data[0]) + + (size_t)i * dis.stride[0])[j]; + const float d_u = (reinterpret_cast(dis.data[1]) + + (size_t)ci * dis.stride[1])[cj]; + const float d_v = (reinterpret_cast(dis.data[2]) + + (size_t)ci * dis.stride[2])[cj]; + + const LABColorf color_1 = get_lab_color(r_y, r_u, r_v, 8); + const LABColorf color_2 = get_lab_color(d_y, d_u, d_v, 8); + de00 = ciede2000(color_1, color_2); + } + + block_reduce_de(de00, partials); +} + +__global__ void ciede_kernel_16bpc(const VmafPicture ref, const VmafPicture dis, + VmafCudaBuffer partials, unsigned width, unsigned height, + int ss_hor, int ss_ver) +{ + const int j = blockIdx.x * blockDim.x + threadIdx.x; + const int i = blockIdx.y * blockDim.y + threadIdx.y; + + double de00 = 0.0; + if (j < width && i < height) { + const int cj = chroma_col(j, ss_hor); + const int ci = chroma_row(i, ss_ver); + + const float r_y = reinterpret_cast( + reinterpret_cast(ref.data[0]) + + (size_t)i * ref.stride[0])[j]; + const float r_u = reinterpret_cast( + reinterpret_cast(ref.data[1]) + + (size_t)ci * ref.stride[1])[cj]; + const float r_v = reinterpret_cast( + reinterpret_cast(ref.data[2]) + + (size_t)ci * ref.stride[2])[cj]; + const float d_y = reinterpret_cast( + reinterpret_cast(dis.data[0]) + + (size_t)i * dis.stride[0])[j]; + const float d_u = reinterpret_cast( + reinterpret_cast(dis.data[1]) + + (size_t)ci * dis.stride[1])[cj]; + const float d_v = reinterpret_cast( + reinterpret_cast(dis.data[2]) + + (size_t)ci * dis.stride[2])[cj]; + + const LABColorf color_1 = get_lab_color(r_y, r_u, r_v, ref.bpc); + const LABColorf color_2 = get_lab_color(d_y, d_u, d_v, dis.bpc); + de00 = ciede2000(color_1, color_2); + } + + block_reduce_de(de00, partials); +} + +} diff --git a/libvmaf/src/feature/cuda/ciede_cuda.c b/libvmaf/src/feature/cuda/ciede_cuda.c new file mode 100644 index 000000000..c63d1555b --- /dev/null +++ b/libvmaf/src/feature/cuda/ciede_cuda.c @@ -0,0 +1,257 @@ +/** + * + * Copyright 2026 Bardie Høgh Joensen + * + * Licensed under the BSD+Patent License (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://opensource.org/licenses/BSDplusPatent + * + * 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. + * + */ + +#include +#include +#include +#include +#include + +#include "common.h" +#include "feature_collector.h" +#include "feature_extractor.h" +#include "cuda/ciede_cuda.h" +#include "picture.h" +#include "picture_cuda.h" +#include "cuda_helper.cuh" + +typedef struct CiedeStateCuda { + CUevent finished, consumed; + CUevent slot_done[2]; + CUfunction funcbpc8, funcbpc16; + CUstream str, host_stream; + VmafCudaBuffer *partials; + double *partials_host; + void *write_score_parameters; + unsigned w, h; + unsigned bpc; + unsigned n_partials; + int ss_hor, ss_ver; +} CiedeStateCuda; + +typedef struct write_score_parameters_ciede { + VmafFeatureCollector *feature_collector; + CiedeStateCuda *s; + const double *partials; + unsigned index; +} write_score_parameters_ciede; + +static int init_fex_cuda(VmafFeatureExtractor *fex, enum VmafPixelFormat pix_fmt, + unsigned bpc, unsigned w, unsigned h) +{ + CiedeStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + + if (pix_fmt == VMAF_PIX_FMT_YUV400P) + return -EINVAL; + switch (bpc) { + case 8: + case 10: + case 12: + case 16: + break; + default: + return -EINVAL; + } + + CHECK_CUDA(cu_f, cuCtxPushCurrent(fex->cu_state->ctx)); + // the work stream is deliberately legacy-blocking: producers like the + // ffmpeg libvmaf_cuda filter fill device pictures with synchronous-API + // copies that are queued on the legacy NULL stream (device-to-device + // memcpy does not block the host), and only blocking-flavor streams are + // implicitly ordered after NULL-stream work. Other extractors' created + // streams are unaffected, so kernel overlap with them is preserved. + CHECK_CUDA(cu_f, cuStreamCreateWithPriority(&s->str, CU_STREAM_DEFAULT, 0)); + CHECK_CUDA(cu_f, cuStreamCreateWithPriority(&s->host_stream, CU_STREAM_NON_BLOCKING, 0)); + CHECK_CUDA(cu_f, cuEventCreate(&s->finished, CU_EVENT_DEFAULT)); + CHECK_CUDA(cu_f, cuEventCreate(&s->consumed, CU_EVENT_DEFAULT)); + CHECK_CUDA(cu_f, cuEventCreate(&s->slot_done[0], CU_EVENT_DEFAULT)); + CHECK_CUDA(cu_f, cuEventCreate(&s->slot_done[1], CU_EVENT_DEFAULT)); + + CUmodule module; + CHECK_CUDA(cu_f, cuModuleLoadData(&module, ciede_ptx)); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->funcbpc8, module, "ciede_kernel_8bpc")); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->funcbpc16, module, "ciede_kernel_16bpc")); + + CHECK_CUDA(cu_f, cuCtxPopCurrent(NULL)); + + s->w = w; + s->h = h; + s->bpc = bpc; + s->ss_hor = pix_fmt != VMAF_PIX_FMT_YUV444P; + s->ss_ver = pix_fmt == VMAF_PIX_FMT_YUV420P; + s->n_partials = DIV_ROUND_UP(w, 16) * DIV_ROUND_UP(h, 16); + + int ret = 0; + + // two write_score slots + two pinned readback slots so frame i+1 never + // has to wait for frame i's host callback (see slot_done in extract) + s->write_score_parameters = malloc(sizeof(write_score_parameters_ciede) * 2); + if (!s->write_score_parameters) return -ENOMEM; + for (unsigned i = 0; i < 2; i++) + ((write_score_parameters_ciede*)s->write_score_parameters)[i].s = s; + + ret |= vmaf_cuda_buffer_alloc(fex->cu_state, &s->partials, + sizeof(double) * s->n_partials); + ret |= vmaf_cuda_buffer_host_alloc(fex->cu_state, (void**)&s->partials_host, + sizeof(double) * s->n_partials * 2); + if (ret) return -ENOMEM; + + return 0; +} + +static int write_scores(write_score_parameters_ciede *params) +{ + CiedeStateCuda *s = params->s; + VmafFeatureCollector *feature_collector = params->feature_collector; + + // sequential sum over block partials keeps the result deterministic + double de00_sum = 0.; + for (unsigned b = 0; b < s->n_partials; b++) + de00_sum += params->partials[b]; + + // identical frames give de00_sum == 0 and score +inf, like the CPU fex + const double score = 45. - 20. * + log10(de00_sum / ((double)s->w * s->h)); + return vmaf_feature_collector_append(feature_collector, "ciede2000", score, + params->index); +} + +static int extract_fex_cuda(VmafFeatureExtractor *fex, VmafPicture *ref_pic, + VmafPicture *ref_pic_90, VmafPicture *dist_pic, + VmafPicture *dist_pic_90, unsigned index, + VmafFeatureCollector *feature_collector) +{ + CiedeStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + + (void) ref_pic_90; + (void) dist_pic_90; + + // two slots: wait for frame index-2's host callback (effectively always + // complete) instead of stalling on the whole previous frame's work + const unsigned slot = index & 1; + CHECK_CUDA(cu_f, cuEventSynchronize(s->slot_done[slot])); + + // kernels run on the extractor's own stream so they overlap with other + // extractors' work on the picture streams; wait for both uploads first + CHECK_CUDA(cu_f, cuStreamWaitEvent(s->str, + vmaf_cuda_picture_get_ready_event(ref_pic), + CU_EVENT_WAIT_DEFAULT)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(s->str, + vmaf_cuda_picture_get_ready_event(dist_pic), + CU_EVENT_WAIT_DEFAULT)); + + { + unsigned width = s->w, height = s->h; + int ss_hor = s->ss_hor, ss_ver = s->ss_ver; + void *kernel_params[] = { + (void*) ref_pic, (void*) dist_pic, (void*) s->partials, + &width, &height, &ss_hor, &ss_ver, + }; + const CUfunction func = + (ref_pic->bpc == 8) ? s->funcbpc8 : s->funcbpc16; + CHECK_CUDA(cu_f, cuLaunchKernel(func, + DIV_ROUND_UP(width, 16), DIV_ROUND_UP(height, 16), 1, + 16, 16, 1, 0, s->str, kernel_params, NULL)); + } + + // lifetime handshake: the pool recycles a picture once the `finished` + // event its own stream records (after the fex loop) has completed, so + // make both picture streams wait for our reads + CHECK_CUDA(cu_f, cuEventRecord(s->consumed, s->str)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(vmaf_cuda_picture_get_stream(ref_pic), + s->consumed, CU_EVENT_WAIT_DEFAULT)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(vmaf_cuda_picture_get_stream(dist_pic), + s->consumed, CU_EVENT_WAIT_DEFAULT)); + + // Download block partials into this slot's readback segment + double *partials_host = s->partials_host + slot * s->n_partials; + CHECK_CUDA(cu_f, cuMemcpyDtoHAsync(partials_host, s->partials->data, + sizeof(double) * s->n_partials, s->str)); + CHECK_CUDA(cu_f, cuEventRecord(s->finished, s->str)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(s->host_stream, s->finished, + CU_EVENT_WAIT_DEFAULT)); + + write_score_parameters_ciede *params = + &((write_score_parameters_ciede*)s->write_score_parameters)[slot]; + params->feature_collector = feature_collector; + params->partials = partials_host; + params->index = index; + CHECK_CUDA(cu_f, cuLaunchHostFunc(s->host_stream, (CUhostFn*)write_scores, + params)); + CHECK_CUDA(cu_f, cuEventRecord(s->slot_done[slot], s->host_stream)); + + return 0; +} + +static int flush_fex_cuda(VmafFeatureExtractor *fex, + VmafFeatureCollector *feature_collector) +{ + (void)feature_collector; + CiedeStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + + // drain the pending write_scores host callback so the final frame's + // score is in the collector before anything reads it + CHECK_CUDA(cu_f, cuStreamSynchronize(s->str)); + CHECK_CUDA(cu_f, cuStreamSynchronize(s->host_stream)); + return 1; +} + +static int close_fex_cuda(VmafFeatureExtractor *fex) +{ + CiedeStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + CHECK_CUDA(cu_f, cuStreamSynchronize(s->str)); + CHECK_CUDA(cu_f, cuStreamSynchronize(s->host_stream)); + CHECK_CUDA(cu_f, cuEventDestroy(s->finished)); + CHECK_CUDA(cu_f, cuEventDestroy(s->consumed)); + CHECK_CUDA(cu_f, cuEventDestroy(s->slot_done[0])); + CHECK_CUDA(cu_f, cuEventDestroy(s->slot_done[1])); + CHECK_CUDA(cu_f, cuStreamDestroy(s->str)); + CHECK_CUDA(cu_f, cuStreamDestroy(s->host_stream)); + + int ret = 0; + if (s->partials) { + ret |= vmaf_cuda_buffer_free(fex->cu_state, s->partials); + free(s->partials); + } + if (s->partials_host) + ret |= vmaf_cuda_buffer_host_free(fex->cu_state, s->partials_host); + if (s->write_score_parameters) + free(s->write_score_parameters); + + return ret; +} + +static const char *provided_features[] = { + "ciede2000", + NULL +}; + +VmafFeatureExtractor vmaf_fex_ciede_cuda = { + .name = "ciede_cuda", + .init = init_fex_cuda, + .extract = extract_fex_cuda, + .flush = flush_fex_cuda, + .close = close_fex_cuda, + .priv_size = sizeof(CiedeStateCuda), + .provided_features = provided_features, + .flags = VMAF_FEATURE_EXTRACTOR_CUDA | VMAF_FEATURE_EXTRACTOR_CUDA_CHROMA, +}; diff --git a/libvmaf/src/feature/cuda/ciede_cuda.h b/libvmaf/src/feature/cuda/ciede_cuda.h new file mode 100644 index 000000000..6710a8310 --- /dev/null +++ b/libvmaf/src/feature/cuda/ciede_cuda.h @@ -0,0 +1,26 @@ +/** + * + * Copyright 2026 Bardie Høgh Joensen + * + * Licensed under the BSD+Patent License (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://opensource.org/licenses/BSDplusPatent + * + * 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. + * + */ + +#ifndef FEATURE_CIEDE_CUDA_H_ +#define FEATURE_CIEDE_CUDA_H_ + +#include +#include "common.h" + +extern const unsigned char ciede_ptx[]; +#endif /* _FEATURE_CIEDE_CUDA_H_ */ diff --git a/libvmaf/src/feature/cuda/float_ssim/ssim.cu b/libvmaf/src/feature/cuda/float_ssim/ssim.cu new file mode 100644 index 000000000..c1960717e --- /dev/null +++ b/libvmaf/src/feature/cuda/float_ssim/ssim.cu @@ -0,0 +1,287 @@ +/** + * + * Copyright 2026 Bardie Høgh Joensen + * + * Licensed under the BSD+Patent License (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://opensource.org/licenses/BSDplusPatent + * + * 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. + * + */ + +/* + * GPU port of the float_ssim feature extractor (feature/ssim.c + feature/iqa). + * + * The kernels below replicate the iqa reference operation-for-operation so + * per-pixel intermediates match the CPU extractor to float precision: + * - decimation samples at (x*factor, y*factor) with a box kernel and + * symmetric boundary mirroring (iqa/decimate.c, _iqa_filter_pixel) + * - the 11-tap separable Gaussian runs in valid mode (output shrinks by + * kernel-1) with double accumulation per 1-D pass and a float round + * in between, exactly like _iqa_convolve with IQA_CONVOLVE_1D + * - per-pixel l/c/s uses the same float/double mixing as _iqa_ssim; + * float-typed reference operations use explicit __f*_rn intrinsics so + * nvcc cannot contract them into fma with different rounding + * The frame mean is a per-block tree reduction; the block partials are + * summed sequentially on the host, so results are deterministic run-to-run. + */ + +#include "cuda_helper.cuh" + +#include "common.h" + +__constant__ float gaussian_k[11] = { + 0.001028f, 0.007599f, 0.036001f, 0.109361f, 0.213006f, 0.266012f, + 0.213006f, 0.109361f, 0.036001f, 0.007599f, 0.001028f, +}; + +// Mirrors an idx along its valid [0, sup) range like iqa's KBND_SYMMETRIC +__device__ __forceinline__ int mirror_sym(int idx, const int sup) +{ + if (idx < 0) idx = -1 - idx; + else if (idx >= sup) idx = (sup - (idx - sup)) - 1; + return idx; +} + +extern "C" { + +__global__ void ssim_normalize_8bpc(const VmafPicture pic, VmafCudaBuffer out, + unsigned w, unsigned h) +{ + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + if (x >= w || y >= h) return; + + const uint8_t v = (reinterpret_cast(pic.data[0]) + + y * pic.stride[0])[x]; + reinterpret_cast(out.data)[y * w + x] = static_cast(v); +} + +__global__ void ssim_normalize_16bpc(const VmafPicture pic, VmafCudaBuffer out, + unsigned w, unsigned h, float scaler) +{ + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + if (x >= w || y >= h) return; + + const uint16_t v = reinterpret_cast( + reinterpret_cast(pic.data[0]) + + y * pic.stride[0])[x]; + // picture_copy: (float)data / scaler; scaler is a power of two, exact + reinterpret_cast(out.data)[y * w + x] = + __fdiv_rn(static_cast(v), scaler); +} + +// iqa/decimate.c: out[y*sw+x] = box(in, x*factor, y*factor), symmetric +// bounds — fused with the float conversion of picture_copy so the full-res +// float image is never materialized. The conversion is exact per sample +// (power-of-two divide), so the box filter sees bit-identical floats. +__global__ void ssim_decimate_8bpc(const VmafPicture pic, VmafCudaBuffer out, + int w, int h, int sw, int sh, int factor) +{ + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + if (x >= sw || y >= sh) return; + + const uint8_t *base = reinterpret_cast(pic.data[0]); + const int cx = x * factor; + const int cy = y * factor; + const int uc = factor / 2; + const int k_even = (factor & 1) ? 0 : 1; + const float k_val = 1.0f / (factor * factor); + + double sum = 0.0; + if (cx >= uc && cy >= uc && cx < w - uc && cy < h - uc) { + for (int v = -uc; v <= uc - k_even; v++) { + const uint8_t *row = base + (cy + v) * pic.stride[0] + cx; + for (int u = -uc; u <= uc - k_even; u++) { + const float px = static_cast(row[u]); + sum += px * k_val; + } + } + } else { + for (int v = -uc; v <= uc - k_even; v++) { + const int yy = mirror_sym(cy + v, h); + const uint8_t *row = base + yy * pic.stride[0]; + for (int u = -uc; u <= uc - k_even; u++) { + const float px = static_cast(row[mirror_sym(cx + u, w)]); + sum += px * k_val; + } + } + } + // _iqa_filter_pixel returns (float)(sum * kscale) with kscale == 1.0 + reinterpret_cast(out.data)[y * sw + x] = static_cast(sum); +} + +__global__ void ssim_decimate_16bpc(const VmafPicture pic, VmafCudaBuffer out, + int w, int h, int sw, int sh, int factor, float scaler) +{ + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + if (x >= sw || y >= sh) return; + + const uint8_t *base = reinterpret_cast(pic.data[0]); + const int cx = x * factor; + const int cy = y * factor; + const int uc = factor / 2; + const int k_even = (factor & 1) ? 0 : 1; + const float k_val = 1.0f / (factor * factor); + + double sum = 0.0; + if (cx >= uc && cy >= uc && cx < w - uc && cy < h - uc) { + for (int v = -uc; v <= uc - k_even; v++) { + const uint16_t *row = reinterpret_cast( + base + (cy + v) * pic.stride[0]) + cx; + for (int u = -uc; u <= uc - k_even; u++) { + const float px = __fdiv_rn(static_cast(row[u]), scaler); + sum += px * k_val; + } + } + } else { + for (int v = -uc; v <= uc - k_even; v++) { + const int yy = mirror_sym(cy + v, h); + const uint16_t *row = reinterpret_cast( + base + yy * pic.stride[0]); + for (int u = -uc; u <= uc - k_even; u++) { + const float px = __fdiv_rn(static_cast( + row[mirror_sym(cx + u, w)]), scaler); + sum += px * k_val; + } + } + } + reinterpret_cast(out.data)[y * sw + x] = static_cast(sum); +} + +__global__ void ssim_products(const VmafCudaBuffer ref, const VmafCudaBuffer cmp, + VmafCudaBuffer ref2, VmafCudaBuffer cmp2, VmafCudaBuffer both, int n) +{ + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n) return; + + const float r = reinterpret_cast(ref.data)[i]; + const float c = reinterpret_cast(cmp.data)[i]; + reinterpret_cast(ref2.data)[i] = __fmul_rn(r, r); + reinterpret_cast(cmp2.data)[i] = __fmul_rn(c, c); + reinterpret_cast(both.data)[i] = __fmul_rn(r, c); +} + +// _iqa_convolve horizontal pass (valid in x): rows span the full height so +// the vertical pass has its apron, exactly like the CPU img_cache +__global__ void ssim_conv_h(const VmafCudaBuffer in, VmafCudaBuffer cache, + int w, int h, int dst_w) +{ + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int ky = blockIdx.y * blockDim.y + threadIdx.y; + if (x >= dst_w || ky >= h) return; + + const float *img = reinterpret_cast(in.data); + const int kx = x + 5; + const int img_offset = ky * w + kx; + + double sum = 0.0; +#pragma unroll + for (int u = -5; u <= 5; u++) + sum += img[img_offset + u] * gaussian_k[u + 5]; + reinterpret_cast(cache.data)[img_offset] = static_cast(sum); +} + +__global__ void ssim_conv_v(const VmafCudaBuffer cache, VmafCudaBuffer out, + int w, int dst_w, int dst_h) +{ + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + if (x >= dst_w || y >= dst_h) return; + + const float *c = reinterpret_cast(cache.data); + const int img_offset = (y + 5) * w + x + 5; + + double sum = 0.0; +#pragma unroll + for (int v = -5; v <= 5; v++) + sum += c[img_offset + v * w] * gaussian_k[v + 5]; + reinterpret_cast(out.data)[y * dst_w + x] = static_cast(sum); +} + +// _iqa_ssim per-pixel l/c/s and frame sums; one partial per block, laid out +// as partials[block][4] = {ssim, l, c, s} +__global__ void ssim_map_reduce(const VmafCudaBuffer mu1_buf, + const VmafCudaBuffer mu2_buf, const VmafCudaBuffer cref2_buf, + const VmafCudaBuffer ccmp2_buf, const VmafCudaBuffer cboth_buf, + VmafCudaBuffer partials, int n, float c1, float c2, float c3) +{ + __shared__ double sh[256][4]; + + const int i = blockIdx.x * blockDim.x + threadIdx.x; + const int t = threadIdx.x; + + double lcs = 0.0, l = 0.0, c = 0.0, s = 0.0; + if (i < n) { + const float mu1 = reinterpret_cast(mu1_buf.data)[i]; + const float mu2 = reinterpret_cast(mu2_buf.data)[i]; + + // sigma buffers are float and computed as conv(x^2) - mu^2 in float + float ref_sigma_sqd = __fsub_rn( + reinterpret_cast(cref2_buf.data)[i], + __fmul_rn(mu1, mu1)); + float cmp_sigma_sqd = __fsub_rn( + reinterpret_cast(ccmp2_buf.data)[i], + __fmul_rn(mu2, mu2)); + ref_sigma_sqd = MAX(0.0f, ref_sigma_sqd); + cmp_sigma_sqd = MAX(0.0f, cmp_sigma_sqd); + const float sigma_both = __fsub_rn( + reinterpret_cast(cboth_buf.data)[i], + __fmul_rn(mu1, mu2)); + + // float sqrt over a double-promoted float product, as in _iqa_ssim + const float sigma_ref_sigma_cmp = static_cast(sqrt( + static_cast(__fmul_rn(ref_sigma_sqd, cmp_sigma_sqd)))); + + // l and c divide a double numerator by a float-summed denominator + l = (2.0 * mu1 * mu2 + c1) / + static_cast(__fadd_rn(__fadd_rn(__fmul_rn(mu1, mu1), + __fmul_rn(mu2, mu2)), c1)); + c = (2.0 * sigma_ref_sigma_cmp + c2) / + static_cast(__fadd_rn(__fadd_rn(ref_sigma_sqd, + cmp_sigma_sqd), c2)); + + const float clamped_sigma_both = + (sigma_both < 0.0f && sigma_ref_sigma_cmp <= 0.0f) ? + 0.0f : sigma_both; + // s is a float division in the reference, promoted on assignment + s = static_cast(__fdiv_rn(__fadd_rn(clamped_sigma_both, c3), + __fadd_rn(sigma_ref_sigma_cmp, c3))); + + lcs = l * c * s; + } + + sh[t][0] = lcs; + sh[t][1] = l; + sh[t][2] = c; + sh[t][3] = s; + __syncthreads(); + + for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) { + if (t < stride) { +#pragma unroll + for (int j = 0; j < 4; j++) + sh[t][j] += sh[t + stride][j]; + } + __syncthreads(); + } + + if (t == 0) { +#pragma unroll + for (int j = 0; j < 4; j++) + reinterpret_cast(partials.data)[blockIdx.x * 4 + j] = + sh[0][j]; + } +} + +} diff --git a/libvmaf/src/feature/cuda/float_ssim_cuda.c b/libvmaf/src/feature/cuda/float_ssim_cuda.c new file mode 100644 index 000000000..8da95ef1b --- /dev/null +++ b/libvmaf/src/feature/cuda/float_ssim_cuda.c @@ -0,0 +1,491 @@ +/** + * + * Copyright 2026 Bardie Høgh Joensen + * + * Licensed under the BSD+Patent License (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://opensource.org/licenses/BSDplusPatent + * + * 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. + * + */ + +#include +#include +#include +#include +#include + +#include "common.h" +#include "feature_collector.h" +#include "feature_extractor.h" +#include "cuda/float_ssim_cuda.h" +#include "opt.h" +#include "picture.h" +#include "picture_cuda.h" +#include "cuda_helper.cuh" + +#define GAUSSIAN_LEN 11 +#define REDUCE_BLOCK 256 + +typedef struct SsimStateCuda { + CUevent finished, consumed; + CUevent slot_done[2]; + CUfunction f_norm8, f_norm16, f_dec8, f_dec16, f_products; + CUfunction f_conv_h, f_conv_v, f_map_reduce; + CUstream str, host_stream; + VmafCudaBuffer *ref_f, *cmp_f; + VmafCudaBuffer *refd, *cmpd; + VmafCudaBuffer *ref2, *cmp2, *both; + VmafCudaBuffer *cache; + VmafCudaBuffer *mu1, *mu2, *cref2, *ccmp2, *cboth; + VmafCudaBuffer *partials; + double *partials_host; + void *write_score_parameters; + unsigned w, h, sw, sh, cw, ch; + unsigned n_blocks; + unsigned bpc; + int factor; + bool enable_lcs; + bool enable_db; + bool clip_db; + double max_db; + int scale; +} SsimStateCuda; + +static const VmafOption options[] = { + { + .name = "enable_lcs", + .help = "enable luminance, contrast and structure intermediate output", + .offset = offsetof(SsimStateCuda, enable_lcs), + .type = VMAF_OPT_TYPE_BOOL, + .default_val.b = false, + }, + { + .name = "enable_db", + .help = "write SSIM values as dB", + .offset = offsetof(SsimStateCuda, enable_db), + .type = VMAF_OPT_TYPE_BOOL, + .default_val.b = false, + }, + { + .name = "clip_db", + .help = "clip dB scores", + .offset = offsetof(SsimStateCuda, clip_db), + .type = VMAF_OPT_TYPE_BOOL, + .default_val.b = false, + }, + { + .name = "scale", + .help = "decimation scale factor (0=auto, 1=no downscaling, 2-10=explicit)", + .offset = offsetof(SsimStateCuda, scale), + .type = VMAF_OPT_TYPE_INT, + .default_val.i = 0, + .min = 0, + .max = 10, + }, + { 0 } +}; + +typedef struct write_score_parameters_ssim { + VmafFeatureCollector *feature_collector; + SsimStateCuda *s; + const double *partials; + unsigned index; +} write_score_parameters_ssim; + +// iqa/math_utils.c _round: round half away from zero +static int iqa_round(float a) +{ + int sign_a = a > 0.0f ? 1 : -1; + return a - (int)a >= 0.5 ? (int)a + sign_a : (int)a; +} + +static int alloc_buf(VmafFeatureExtractor *fex, VmafCudaBuffer **buf, + size_t size) +{ + return vmaf_cuda_buffer_alloc(fex->cu_state, buf, size); +} + +static int init_fex_cuda(VmafFeatureExtractor *fex, enum VmafPixelFormat pix_fmt, + unsigned bpc, unsigned w, unsigned h) +{ + SsimStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + + (void) pix_fmt; + + CHECK_CUDA(cu_f, cuCtxPushCurrent(fex->cu_state->ctx)); + // the work stream is deliberately legacy-blocking: producers like the + // ffmpeg libvmaf_cuda filter fill device pictures with synchronous-API + // copies that are queued on the legacy NULL stream (device-to-device + // memcpy does not block the host), and only blocking-flavor streams are + // implicitly ordered after NULL-stream work. Other extractors' created + // streams are unaffected, so kernel overlap with them is preserved. + CHECK_CUDA(cu_f, cuStreamCreateWithPriority(&s->str, CU_STREAM_DEFAULT, 0)); + CHECK_CUDA(cu_f, cuStreamCreateWithPriority(&s->host_stream, CU_STREAM_NON_BLOCKING, 0)); + CHECK_CUDA(cu_f, cuEventCreate(&s->finished, CU_EVENT_DEFAULT)); + CHECK_CUDA(cu_f, cuEventCreate(&s->consumed, CU_EVENT_DEFAULT)); + CHECK_CUDA(cu_f, cuEventCreate(&s->slot_done[0], CU_EVENT_DEFAULT)); + CHECK_CUDA(cu_f, cuEventCreate(&s->slot_done[1], CU_EVENT_DEFAULT)); + + CUmodule module; + CHECK_CUDA(cu_f, cuModuleLoadData(&module, ssim_ptx)); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->f_norm8, module, "ssim_normalize_8bpc")); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->f_norm16, module, "ssim_normalize_16bpc")); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->f_dec8, module, "ssim_decimate_8bpc")); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->f_dec16, module, "ssim_decimate_16bpc")); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->f_products, module, "ssim_products")); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->f_conv_h, module, "ssim_conv_h")); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->f_conv_v, module, "ssim_conv_v")); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->f_map_reduce, module, "ssim_map_reduce")); + + CHECK_CUDA(cu_f, cuCtxPopCurrent(NULL)); + + s->w = w; + s->h = h; + s->bpc = bpc; + + // compute_ssim: scale = max(1, round(min(w,h) / 256.0)), or the override + const unsigned min_wh = w < h ? w : h; + s->factor = s->scale > 0 ? + s->scale : (iqa_round((float)min_wh / 256.0f) < 1 ? + 1 : iqa_round((float)min_wh / 256.0f)); + + if (s->factor > 1) { + // _iqa_decimate: sw = w/factor + (w&1) + s->sw = w / s->factor + (w & 1); + s->sh = h / s->factor + (h & 1); + } else { + s->sw = w; + s->sh = h; + } + if (s->sw < GAUSSIAN_LEN || s->sh < GAUSSIAN_LEN) + return -EINVAL; + s->cw = s->sw - GAUSSIAN_LEN + 1; + s->ch = s->sh - GAUSSIAN_LEN + 1; + s->n_blocks = DIV_ROUND_UP(s->cw * s->ch, REDUCE_BLOCK); + + const unsigned peak = (1 << bpc) - 1; + if (s->clip_db) { + const double mse = 0.5 / (w * h); + s->max_db = ceil(10. * log10(peak * peak / mse)); + } else { + s->max_db = INFINITY; + } + + int ret = 0; + + // two write_score slots + two pinned readback slots so frame i+1 never + // has to wait for frame i's host callback (see slot_done in extract) + s->write_score_parameters = malloc(sizeof(write_score_parameters_ssim) * 2); + if (!s->write_score_parameters) return -ENOMEM; + for (unsigned i = 0; i < 2; i++) + ((write_score_parameters_ssim*)s->write_score_parameters)[i].s = s; + + const size_t full = sizeof(float) * w * h; + const size_t dec = sizeof(float) * s->sw * s->sh; + const size_t conv = sizeof(float) * s->cw * s->ch; + + // factor>1 reads pictures directly in the fused decimate kernels, so + // the full-resolution float images are only needed at factor==1 + if (s->factor > 1) { + ret |= alloc_buf(fex, &s->refd, dec); + ret |= alloc_buf(fex, &s->cmpd, dec); + } else { + ret |= alloc_buf(fex, &s->ref_f, full); + ret |= alloc_buf(fex, &s->cmp_f, full); + } + ret |= alloc_buf(fex, &s->ref2, dec); + ret |= alloc_buf(fex, &s->cmp2, dec); + ret |= alloc_buf(fex, &s->both, dec); + ret |= alloc_buf(fex, &s->cache, dec); + ret |= alloc_buf(fex, &s->mu1, conv); + ret |= alloc_buf(fex, &s->mu2, conv); + ret |= alloc_buf(fex, &s->cref2, conv); + ret |= alloc_buf(fex, &s->ccmp2, conv); + ret |= alloc_buf(fex, &s->cboth, conv); + ret |= alloc_buf(fex, &s->partials, sizeof(double) * 4 * s->n_blocks); + ret |= vmaf_cuda_buffer_host_alloc(fex->cu_state, (void**)&s->partials_host, + sizeof(double) * 4 * s->n_blocks * 2); + if (ret) return -ENOMEM; + + return 0; +} + +#define MIN(x, y) (((x) < (y)) ? (x) : (y)) + +static double convert_to_db(double score, double max_db) +{ + return MIN(-10. * log10(1 - score), max_db); +} + +static int write_scores(write_score_parameters_ssim *params) +{ + SsimStateCuda *s = params->s; + VmafFeatureCollector *feature_collector = params->feature_collector; + + // sequential sum over block partials keeps the result deterministic + double ssim_sum = 0., l_sum = 0., c_sum = 0., s_sum = 0.; + for (unsigned b = 0; b < s->n_blocks; b++) { + ssim_sum += params->partials[b * 4 + 0]; + l_sum += params->partials[b * 4 + 1]; + c_sum += params->partials[b * 4 + 2]; + s_sum += params->partials[b * 4 + 3]; + } + + // _iqa_ssim returns float means; compute_ssim widens them to double + const double n = (double)(s->cw * s->ch); + double score = (double)(float)(ssim_sum / n); + const double l_score = (double)(float)(l_sum / n); + const double c_score = (double)(float)(c_sum / n); + const double s_score = (double)(float)(s_sum / n); + + if (s->enable_db) + score = convert_to_db(score, s->max_db); + + int err = vmaf_feature_collector_append(feature_collector, "float_ssim", + score, params->index); + if (s->enable_lcs) { + err |= vmaf_feature_collector_append(feature_collector, "float_ssim_l", + l_score, params->index); + err |= vmaf_feature_collector_append(feature_collector, "float_ssim_c", + c_score, params->index); + err |= vmaf_feature_collector_append(feature_collector, "float_ssim_s", + s_score, params->index); + } + + return err; +} + +static void launch_conv(SsimStateCuda *s, CudaFunctions *cu_f, CUstream stream, + VmafCudaBuffer *in, VmafCudaBuffer *out) +{ + int w = s->sw, h = s->sh, dst_w = s->cw, dst_h = s->ch; + { + void *args[] = { (void*)in, (void*)s->cache, &w, &h, &dst_w }; + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_conv_h, + DIV_ROUND_UP(dst_w, 16), DIV_ROUND_UP(h, 16), 1, + 16, 16, 1, 0, stream, args, NULL)); + } + { + void *args[] = { (void*)s->cache, (void*)out, &w, &dst_w, &dst_h }; + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_conv_v, + DIV_ROUND_UP(dst_w, 16), DIV_ROUND_UP(dst_h, 16), 1, + 16, 16, 1, 0, stream, args, NULL)); + } +} + +static int extract_fex_cuda(VmafFeatureExtractor *fex, VmafPicture *ref_pic, + VmafPicture *ref_pic_90, VmafPicture *dist_pic, + VmafPicture *dist_pic_90, unsigned index, + VmafFeatureCollector *feature_collector) +{ + SsimStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + + (void) ref_pic_90; + (void) dist_pic_90; + + // two slots: wait for frame index-2's host callback (effectively always + // complete) instead of stalling on the whole previous frame's work + const unsigned slot = index & 1; + CHECK_CUDA(cu_f, cuEventSynchronize(s->slot_done[slot])); + + // kernels run on the extractor's own stream so they overlap with other + // extractors' work on the picture streams; wait for both uploads first + CHECK_CUDA(cu_f, cuStreamWaitEvent(s->str, + vmaf_cuda_picture_get_ready_event(ref_pic), + CU_EVENT_WAIT_DEFAULT)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(s->str, + vmaf_cuda_picture_get_ready_event(dist_pic), + CU_EVENT_WAIT_DEFAULT)); + + unsigned w = s->w, h = s->h; + float scaler = 4.0f; + if (s->bpc == 12) scaler = 16.0f; + if (s->bpc == 16) scaler = 256.0f; + + VmafCudaBuffer *ref_in, *cmp_in; + if (s->factor > 1) { + // fused normalize+decimate reads the pictures directly + int iw = w, ih = h, sw = s->sw, sh = s->sh, factor = s->factor; + if (s->bpc == 8) { + void *a1[] = { (void*)ref_pic, (void*)s->refd, &iw, &ih, &sw, &sh, &factor }; + void *a2[] = { (void*)dist_pic, (void*)s->cmpd, &iw, &ih, &sw, &sh, &factor }; + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_dec8, DIV_ROUND_UP(sw, 16), + DIV_ROUND_UP(sh, 16), 1, 16, 16, 1, 0, s->str, a1, NULL)); + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_dec8, DIV_ROUND_UP(sw, 16), + DIV_ROUND_UP(sh, 16), 1, 16, 16, 1, 0, s->str, a2, NULL)); + } else { + void *a1[] = { (void*)ref_pic, (void*)s->refd, &iw, &ih, &sw, &sh, &factor, &scaler }; + void *a2[] = { (void*)dist_pic, (void*)s->cmpd, &iw, &ih, &sw, &sh, &factor, &scaler }; + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_dec16, DIV_ROUND_UP(sw, 16), + DIV_ROUND_UP(sh, 16), 1, 16, 16, 1, 0, s->str, a1, NULL)); + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_dec16, DIV_ROUND_UP(sw, 16), + DIV_ROUND_UP(sh, 16), 1, 16, 16, 1, 0, s->str, a2, NULL)); + } + ref_in = s->refd; + cmp_in = s->cmpd; + } else { + if (s->bpc == 8) { + void *a1[] = { (void*)ref_pic, (void*)s->ref_f, &w, &h }; + void *a2[] = { (void*)dist_pic, (void*)s->cmp_f, &w, &h }; + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_norm8, DIV_ROUND_UP(w, 16), + DIV_ROUND_UP(h, 16), 1, 16, 16, 1, 0, s->str, a1, NULL)); + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_norm8, DIV_ROUND_UP(w, 16), + DIV_ROUND_UP(h, 16), 1, 16, 16, 1, 0, s->str, a2, NULL)); + } else { + void *a1[] = { (void*)ref_pic, (void*)s->ref_f, &w, &h, &scaler }; + void *a2[] = { (void*)dist_pic, (void*)s->cmp_f, &w, &h, &scaler }; + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_norm16, DIV_ROUND_UP(w, 16), + DIV_ROUND_UP(h, 16), 1, 16, 16, 1, 0, s->str, a1, NULL)); + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_norm16, DIV_ROUND_UP(w, 16), + DIV_ROUND_UP(h, 16), 1, 16, 16, 1, 0, s->str, a2, NULL)); + } + ref_in = s->ref_f; + cmp_in = s->cmp_f; + } + + // lifetime handshake right after the only kernels that read picture + // memory: the pool recycles a picture once the `finished` event its own + // stream records (after the fex loop) has completed, so make both + // picture streams wait for those reads — everything below works on our + // own buffers and can outlive the pictures + CHECK_CUDA(cu_f, cuEventRecord(s->consumed, s->str)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(vmaf_cuda_picture_get_stream(ref_pic), + s->consumed, CU_EVENT_WAIT_DEFAULT)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(vmaf_cuda_picture_get_stream(dist_pic), + s->consumed, CU_EVENT_WAIT_DEFAULT)); + + { + int n = s->sw * s->sh; + void *args[] = { (void*)ref_in, (void*)cmp_in, (void*)s->ref2, + (void*)s->cmp2, (void*)s->both, &n }; + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_products, + DIV_ROUND_UP(n, REDUCE_BLOCK), 1, 1, + REDUCE_BLOCK, 1, 1, 0, s->str, args, NULL)); + } + + launch_conv(s, cu_f, s->str, ref_in, s->mu1); + launch_conv(s, cu_f, s->str, cmp_in, s->mu2); + launch_conv(s, cu_f, s->str, s->ref2, s->cref2); + launch_conv(s, cu_f, s->str, s->cmp2, s->ccmp2); + launch_conv(s, cu_f, s->str, s->both, s->cboth); + + { + // _iqa_ssim: C1 = (K1*L)^2, C2 = (K2*L)^2, C3 = C2/2, L = 255 + float c1 = (0.01f * 255) * (0.01f * 255); + float c2 = (0.03f * 255) * (0.03f * 255); + float c3 = c2 / 2.0f; + int n = s->cw * s->ch; + void *args[] = { (void*)s->mu1, (void*)s->mu2, (void*)s->cref2, + (void*)s->ccmp2, (void*)s->cboth, (void*)s->partials, + &n, &c1, &c2, &c3 }; + CHECK_CUDA(cu_f, cuLaunchKernel(s->f_map_reduce, + s->n_blocks, 1, 1, REDUCE_BLOCK, 1, 1, 0, + s->str, args, NULL)); + } + + // Download block partials into this slot's readback segment + double *partials_host = s->partials_host + slot * 4 * s->n_blocks; + CHECK_CUDA(cu_f, cuMemcpyDtoHAsync(partials_host, s->partials->data, + sizeof(double) * 4 * s->n_blocks, s->str)); + CHECK_CUDA(cu_f, cuEventRecord(s->finished, s->str)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(s->host_stream, s->finished, + CU_EVENT_WAIT_DEFAULT)); + + write_score_parameters_ssim *params = + &((write_score_parameters_ssim*)s->write_score_parameters)[slot]; + params->feature_collector = feature_collector; + params->partials = partials_host; + params->index = index; + CHECK_CUDA(cu_f, cuLaunchHostFunc(s->host_stream, (CUhostFn*)write_scores, + params)); + CHECK_CUDA(cu_f, cuEventRecord(s->slot_done[slot], s->host_stream)); + + return 0; +} + +static int flush_fex_cuda(VmafFeatureExtractor *fex, + VmafFeatureCollector *feature_collector) +{ + (void)feature_collector; + SsimStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + + // drain the pending write_scores host callback so the final frame's + // score is in the collector before anything reads it + CHECK_CUDA(cu_f, cuStreamSynchronize(s->str)); + CHECK_CUDA(cu_f, cuStreamSynchronize(s->host_stream)); + return 1; +} + +static int free_buf(VmafFeatureExtractor *fex, VmafCudaBuffer *buf) +{ + int ret = 0; + if (buf) { + ret = vmaf_cuda_buffer_free(fex->cu_state, buf); + free(buf); + } + return ret; +} + +static int close_fex_cuda(VmafFeatureExtractor *fex) +{ + SsimStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + CHECK_CUDA(cu_f, cuStreamSynchronize(s->str)); + CHECK_CUDA(cu_f, cuStreamSynchronize(s->host_stream)); + CHECK_CUDA(cu_f, cuEventDestroy(s->finished)); + CHECK_CUDA(cu_f, cuEventDestroy(s->consumed)); + CHECK_CUDA(cu_f, cuEventDestroy(s->slot_done[0])); + CHECK_CUDA(cu_f, cuEventDestroy(s->slot_done[1])); + CHECK_CUDA(cu_f, cuStreamDestroy(s->str)); + CHECK_CUDA(cu_f, cuStreamDestroy(s->host_stream)); + + int ret = 0; + ret |= free_buf(fex, s->ref_f); + ret |= free_buf(fex, s->cmp_f); + ret |= free_buf(fex, s->refd); + ret |= free_buf(fex, s->cmpd); + ret |= free_buf(fex, s->ref2); + ret |= free_buf(fex, s->cmp2); + ret |= free_buf(fex, s->both); + ret |= free_buf(fex, s->cache); + ret |= free_buf(fex, s->mu1); + ret |= free_buf(fex, s->mu2); + ret |= free_buf(fex, s->cref2); + ret |= free_buf(fex, s->ccmp2); + ret |= free_buf(fex, s->cboth); + ret |= free_buf(fex, s->partials); + if (s->partials_host) + ret |= vmaf_cuda_buffer_host_free(fex->cu_state, s->partials_host); + if (s->write_score_parameters) + free(s->write_score_parameters); + + return ret; +} + +static const char *provided_features[] = { + "float_ssim", + NULL +}; + +VmafFeatureExtractor vmaf_fex_float_ssim_cuda = { + .name = "ssim_cuda", + .options = options, + .init = init_fex_cuda, + .extract = extract_fex_cuda, + .flush = flush_fex_cuda, + .close = close_fex_cuda, + .priv_size = sizeof(SsimStateCuda), + .provided_features = provided_features, + .flags = VMAF_FEATURE_EXTRACTOR_CUDA, +}; diff --git a/libvmaf/src/feature/cuda/float_ssim_cuda.h b/libvmaf/src/feature/cuda/float_ssim_cuda.h new file mode 100644 index 000000000..e71373377 --- /dev/null +++ b/libvmaf/src/feature/cuda/float_ssim_cuda.h @@ -0,0 +1,26 @@ +/** + * + * Copyright 2026 Bardie Høgh Joensen + * + * Licensed under the BSD+Patent License (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://opensource.org/licenses/BSDplusPatent + * + * 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. + * + */ + +#ifndef FEATURE_SSIM_CUDA_H_ +#define FEATURE_SSIM_CUDA_H_ + +#include +#include "common.h" + +extern const unsigned char ssim_ptx[]; +#endif /* _FEATURE_SSIM_CUDA_H_ */ diff --git a/libvmaf/src/feature/cuda/integer_psnr/psnr.cu b/libvmaf/src/feature/cuda/integer_psnr/psnr.cu new file mode 100644 index 000000000..25fd9a3c7 --- /dev/null +++ b/libvmaf/src/feature/cuda/integer_psnr/psnr.cu @@ -0,0 +1,135 @@ +/** + * + * Copyright 2026 Bardie Høgh Joensen + * + * Licensed under the BSD+Patent License (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://opensource.org/licenses/BSDplusPatent + * + * 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. + * + */ + +#include "cuda_helper.cuh" + +#include "common.h" + +// Block-reduce a per-thread value and let thread 0 issue ONE atomicAdd per +// block: one atomic per warp serializes on the single global accumulator, +// one per block is ~32x fewer +__device__ __forceinline__ void block_reduce_add(uint64_t v, + unsigned long long *accum) +{ + __shared__ uint64_t warp_sums[8]; // 256 threads = 8 warps + + const int t = threadIdx.y * blockDim.x + threadIdx.x; +#pragma unroll + for (int i = 16; i > 0; i >>= 1) { + v += uint64_t(__shfl_down_sync(0xffffffff, uint32_t(v), i)) | + (uint64_t(__shfl_down_sync(0xffffffff, uint32_t(v >> 32), i)) << 32); + } + if ((t % 32) == 0) + warp_sums[t / 32] = v; + __syncthreads(); + if (t == 0) { + uint64_t sum = 0; +#pragma unroll + for (int i = 0; i < 8; i++) + sum += warp_sums[i]; + if (sum) + atomicAdd(accum, static_cast(sum)); + } +} + +extern "C" { + +// Grid-stride over uchar4 vectors so each warp reads full 128-byte segments +// (cuMemAllocPitch aligns every row start); the <=3 tail pixels of each row +// when the width isn't a multiple of 4 are covered by a scalar per-row loop. +// Integer sums are order-independent, so the score is unchanged. +__global__ void psnr_kernel_8bpc(const VmafPicture ref, const VmafPicture dis, + VmafCudaBuffer sse, unsigned plane, unsigned width, unsigned height) +{ + const uint8_t *rbase = reinterpret_cast(ref.data[plane]); + const uint8_t *dbase = reinterpret_cast(dis.data[plane]); + const long vecs = width / 4; + const long total = vecs * height; + const long tid = blockIdx.x * (long)blockDim.x + threadIdx.x; + const long step = (long)gridDim.x * blockDim.x; + + uint64_t sq = 0; + for (long i = tid; i < total; i += step) { + const long y = i / vecs; + const long x = (i - y * vecs) * 4; + const uchar4 r = *reinterpret_cast( + rbase + y * ref.stride[plane] + x); + const uchar4 d = *reinterpret_cast( + dbase + y * dis.stride[plane] + x); + int e = r.x - d.x; sq += e * e; + e = r.y - d.y; sq += e * e; + e = r.z - d.z; sq += e * e; + e = r.w - d.w; sq += e * e; + } + + const unsigned tail = width & 3; + if (tail) { + for (long y = tid; y < height; y += step) { + const uint8_t *r = rbase + y * ref.stride[plane]; + const uint8_t *d = dbase + y * dis.stride[plane]; + for (unsigned x = width - tail; x < width; x++) { + const int e = r[x] - d[x]; + sq += e * e; + } + } + } + + block_reduce_add(sq, + reinterpret_cast(sse.data) + plane); +} + +__global__ void psnr_kernel_16bpc(const VmafPicture ref, const VmafPicture dis, + VmafCudaBuffer sse, unsigned plane, unsigned width, unsigned height) +{ + const uint8_t *rbase = reinterpret_cast(ref.data[plane]); + const uint8_t *dbase = reinterpret_cast(dis.data[plane]); + const long vecs = width / 2; + const long total = vecs * height; + const long tid = blockIdx.x * (long)blockDim.x + threadIdx.x; + const long step = (long)gridDim.x * blockDim.x; + + uint64_t sq = 0; + for (long i = tid; i < total; i += step) { + const long y = i / vecs; + const long x = (i - y * vecs) * 4; // byte offset of the ushort2 + const ushort2 r = *reinterpret_cast( + rbase + y * ref.stride[plane] + x); + const ushort2 d = *reinterpret_cast( + dbase + y * dis.stride[plane] + x); + int e = r.x - d.x; + sq += static_cast(static_cast(e) * e); + e = r.y - d.y; + sq += static_cast(static_cast(e) * e); + } + + if (width & 1) { + const unsigned x = width - 1; + for (long y = tid; y < height; y += step) { + const int e = reinterpret_cast( + rbase + y * ref.stride[plane])[x] - + reinterpret_cast( + dbase + y * dis.stride[plane])[x]; + sq += static_cast(static_cast(e) * e); + } + } + + block_reduce_add(sq, + reinterpret_cast(sse.data) + plane); +} + +} diff --git a/libvmaf/src/feature/cuda/integer_psnr_cuda.c b/libvmaf/src/feature/cuda/integer_psnr_cuda.c new file mode 100644 index 000000000..ce9565b02 --- /dev/null +++ b/libvmaf/src/feature/cuda/integer_psnr_cuda.c @@ -0,0 +1,378 @@ +/** + * + * Copyright 2026 Bardie Høgh Joensen + * + * Licensed under the BSD+Patent License (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://opensource.org/licenses/BSDplusPatent + * + * 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. + * + */ + +#include +#include +#include +#include +#include +#include + +#include "common.h" +#include "feature_collector.h" +#include "feature_extractor.h" +#include "cuda/integer_psnr_cuda.h" +#include "opt.h" +#include "picture.h" +#include "picture_cuda.h" +#include "cuda_helper.cuh" + +typedef struct PsnrStateCuda { + CUevent finished, consumed; + CUevent slot_done[2]; + CUfunction funcbpc8, funcbpc16; + CUstream str, host_stream; + VmafCudaBuffer *sse; + uint64_t *sse_host; + void *write_score_parameters; + unsigned bpc; + bool enable_chroma; + bool enable_mse; + bool enable_apsnr; + bool reduced_hbd_peak; + uint32_t peak; + double psnr_max[3]; + double min_sse; + struct { + uint64_t sse[3]; + uint64_t n_pixels[3]; + } apsnr; +} PsnrStateCuda; + +static const VmafOption options[] = { + { + .name = "enable_chroma", + .help = "enable calculation for chroma channels", + .offset = offsetof(PsnrStateCuda, enable_chroma), + .type = VMAF_OPT_TYPE_BOOL, + .default_val.b = true, + }, + { + .name = "enable_mse", + .help = "enable MSE calculation", + .offset = offsetof(PsnrStateCuda, enable_mse), + .type = VMAF_OPT_TYPE_BOOL, + .default_val.b = false, + }, + { + .name = "enable_apsnr", + .help = "enable APSNR calculation", + .offset = offsetof(PsnrStateCuda, enable_apsnr), + .type = VMAF_OPT_TYPE_BOOL, + .default_val.b = false, + }, + { + .name = "reduced_hbd_peak", + .help = "reduce hbd peak value to align with scaled 8-bit content", + .offset = offsetof(PsnrStateCuda, reduced_hbd_peak), + .type = VMAF_OPT_TYPE_BOOL, + .default_val.b = false, + }, + { + .name = "min_sse", + .help = "constrain the minimum possible sse", + .offset = offsetof(PsnrStateCuda, min_sse), + .type = VMAF_OPT_TYPE_DOUBLE, + .default_val.d = 0.0, + .min = 0.0, + .max = DBL_MAX, + }, + { 0 } +}; + +typedef struct write_score_parameters_psnr { + VmafFeatureCollector *feature_collector; + PsnrStateCuda *s; + const uint64_t *sse; + unsigned w[3], h[3]; + unsigned index; +} write_score_parameters_psnr; + +static int init_fex_cuda(VmafFeatureExtractor *fex, enum VmafPixelFormat pix_fmt, + unsigned bpc, unsigned w, unsigned h) +{ + PsnrStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + + CHECK_CUDA(cu_f, cuCtxPushCurrent(fex->cu_state->ctx)); + // the work stream is deliberately legacy-blocking: producers like the + // ffmpeg libvmaf_cuda filter fill device pictures with synchronous-API + // copies that are queued on the legacy NULL stream (device-to-device + // memcpy does not block the host), and only blocking-flavor streams are + // implicitly ordered after NULL-stream work. Other extractors' created + // streams are unaffected, so kernel overlap with them is preserved. + CHECK_CUDA(cu_f, cuStreamCreateWithPriority(&s->str, CU_STREAM_DEFAULT, 0)); + CHECK_CUDA(cu_f, cuStreamCreateWithPriority(&s->host_stream, CU_STREAM_NON_BLOCKING, 0)); + CHECK_CUDA(cu_f, cuEventCreate(&s->finished, CU_EVENT_DEFAULT)); + CHECK_CUDA(cu_f, cuEventCreate(&s->consumed, CU_EVENT_DEFAULT)); + CHECK_CUDA(cu_f, cuEventCreate(&s->slot_done[0], CU_EVENT_DEFAULT)); + CHECK_CUDA(cu_f, cuEventCreate(&s->slot_done[1], CU_EVENT_DEFAULT)); + + CUmodule module; + CHECK_CUDA(cu_f, cuModuleLoadData(&module, psnr_ptx)); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->funcbpc8, module, "psnr_kernel_8bpc")); + CHECK_CUDA(cu_f, cuModuleGetFunction(&s->funcbpc16, module, "psnr_kernel_16bpc")); + + CHECK_CUDA(cu_f, cuCtxPopCurrent(NULL)); + + s->bpc = bpc; + s->peak = s->reduced_hbd_peak ? 255 * 1 << (bpc - 8) : (1 << bpc) - 1; + + if (pix_fmt == VMAF_PIX_FMT_YUV400P) + s->enable_chroma = false; + + for (unsigned i = 0; i < 3; i++) { + if (s->min_sse != 0.0) { + const int ss_hor = pix_fmt != VMAF_PIX_FMT_YUV444P; + const int ss_ver = pix_fmt == VMAF_PIX_FMT_YUV420P; + const double mse = s->min_sse / + (((i && ss_hor) ? w / 2 : w) * ((i && ss_ver) ? h / 2 : h)); + s->psnr_max[i] = ceil(10. * log10(s->peak * s->peak / mse)); + } else { + s->psnr_max[i] = (6 * bpc) + 12; + } + } + + int ret = 0; + + // two write_score slots + two pinned readback slots so frame i+1 never + // has to wait for frame i's host callback (see slot_done in extract) + s->write_score_parameters = malloc(sizeof(write_score_parameters_psnr) * 2); + if (!s->write_score_parameters) goto free_buf; + for (unsigned i = 0; i < 2; i++) + ((write_score_parameters_psnr*)s->write_score_parameters)[i].s = s; + + ret |= vmaf_cuda_buffer_alloc(fex->cu_state, &s->sse, sizeof(uint64_t) * 3); + if (ret) goto free_buf; + ret |= vmaf_cuda_buffer_host_alloc(fex->cu_state, (void**)&s->sse_host, + sizeof(uint64_t) * 3 * 2); + if (ret) goto free_buf; + + return 0; + +free_buf: + if (s->sse) { + ret |= vmaf_cuda_buffer_free(fex->cu_state, s->sse); + free(s->sse); + } + if (s->write_score_parameters) + free(s->write_score_parameters); + + return -ENOMEM; +} + +#define MAX(x, y) (((x) > (y)) ? (x) : (y)) + +static char *mse_name[3] = { "mse_y", "mse_cb", "mse_cr" }; +static char *psnr_name[3] = { "psnr_y", "psnr_cb", "psnr_cr" }; + +static int write_scores(write_score_parameters_psnr *params) +{ + PsnrStateCuda *s = params->s; + VmafFeatureCollector *feature_collector = params->feature_collector; + + const double peak = (s->bpc == 8) ? 255. : (double) s->peak; + const unsigned n = s->enable_chroma ? 3 : 1; + + int err = 0; + for (unsigned p = 0; p < n; p++) { + const uint64_t sse = params->sse[p]; + + if (s->enable_apsnr) { + s->apsnr.sse[p] += sse; + s->apsnr.n_pixels[p] += (uint64_t)params->w[p] * params->h[p]; + } + + const double mse = + ((double) sse) / ((double)params->w[p] * params->h[p]); + const double psnr = + MIN(10. * log10(peak * peak / MAX(mse, 1e-16)), s->psnr_max[p]); + + err |= vmaf_feature_collector_append(feature_collector, psnr_name[p], + psnr, params->index); + if (s->enable_mse) { + err |= vmaf_feature_collector_append(feature_collector, mse_name[p], + mse, params->index); + } + } + + return err; +} + +static int extract_fex_cuda(VmafFeatureExtractor *fex, VmafPicture *ref_pic, + VmafPicture *ref_pic_90, VmafPicture *dist_pic, + VmafPicture *dist_pic_90, unsigned index, + VmafFeatureCollector *feature_collector) +{ + PsnrStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + + (void) ref_pic_90; + (void) dist_pic_90; + + // two slots: wait for frame index-2's host callback (effectively always + // complete) instead of stalling on the whole previous frame's work + const unsigned slot = index & 1; + CHECK_CUDA(cu_f, cuEventSynchronize(s->slot_done[slot])); + + const unsigned n = s->enable_chroma ? 3 : 1; + + // kernels run on the extractor's own stream so they overlap with other + // extractors' work on the picture streams; wait for both uploads first + CHECK_CUDA(cu_f, cuStreamWaitEvent(s->str, + vmaf_cuda_picture_get_ready_event(ref_pic), + CU_EVENT_WAIT_DEFAULT)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(s->str, + vmaf_cuda_picture_get_ready_event(dist_pic), + CU_EVENT_WAIT_DEFAULT)); + CHECK_CUDA(cu_f, cuMemsetD8Async(s->sse->data, 0, sizeof(uint64_t) * 3, + s->str)); + + const CUfunction func = (ref_pic->bpc == 8) ? s->funcbpc8 : s->funcbpc16; + const unsigned vec = (ref_pic->bpc == 8) ? 4 : 2; + for (unsigned p = 0; p < n; p++) { + unsigned plane = p; + unsigned width = ref_pic->w[p]; + unsigned height = ref_pic->h[p]; + // 1-D grid-stride kernel over vectorized loads; cap the grid so + // tail-of-wave blocks stay busy + unsigned n_blocks = DIV_ROUND_UP(width / vec * height, 256); + if (n_blocks > 2048) n_blocks = 2048; + if (!n_blocks) n_blocks = 1; + void *kernel_params[] = { + (void*) ref_pic, (void*) dist_pic, (void*) s->sse, + &plane, &width, &height, + }; + CHECK_CUDA(cu_f, cuLaunchKernel(func, + n_blocks, 1, 1, 256, 1, 1, 0, + s->str, kernel_params, NULL)); + } + + // lifetime handshake: the pool recycles a picture once the `finished` + // event its own stream records (after the fex loop) has completed, so + // make both picture streams wait for our reads + CHECK_CUDA(cu_f, cuEventRecord(s->consumed, s->str)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(vmaf_cuda_picture_get_stream(ref_pic), + s->consumed, CU_EVENT_WAIT_DEFAULT)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(vmaf_cuda_picture_get_stream(dist_pic), + s->consumed, CU_EVENT_WAIT_DEFAULT)); + + // Download sse into this slot's readback segment + uint64_t *sse_host = s->sse_host + slot * 3; + CHECK_CUDA(cu_f, cuMemcpyDtoHAsync(sse_host, s->sse->data, + sizeof(uint64_t) * 3, s->str)); + CHECK_CUDA(cu_f, cuEventRecord(s->finished, s->str)); + CHECK_CUDA(cu_f, cuStreamWaitEvent(s->host_stream, s->finished, + CU_EVENT_WAIT_DEFAULT)); + + write_score_parameters_psnr *params = + &((write_score_parameters_psnr*)s->write_score_parameters)[slot]; + params->feature_collector = feature_collector; + params->sse = sse_host; + for (unsigned p = 0; p < n; p++) { + params->w[p] = ref_pic->w[p]; + params->h[p] = ref_pic->h[p]; + } + params->index = index; + CHECK_CUDA(cu_f, cuLaunchHostFunc(s->host_stream, (CUhostFn*)write_scores, + params)); + CHECK_CUDA(cu_f, cuEventRecord(s->slot_done[slot], s->host_stream)); + + return 0; +} + +static int flush_fex_cuda(VmafFeatureExtractor *fex, + VmafFeatureCollector *feature_collector) +{ + PsnrStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + const char *apsnr_name[3] = { "apsnr_y", "apsnr_cb", "apsnr_cr" }; + + CHECK_CUDA(cu_f, cuStreamSynchronize(s->str)); + CHECK_CUDA(cu_f, cuStreamSynchronize(s->host_stream)); + + // aggregates only: set_aggregate is idempotent, so a second flush from + // the threaded + CUDA flush paths cannot double-append + int err = 0; + if (s->enable_apsnr) { + for (unsigned i = 0; i < 3; i++) { + + double apsnr = 10 * (log10(s->peak * s->peak) + + log10(s->apsnr.n_pixels[i]) - + log10(s->apsnr.sse[i])); + + double max_apsnr = + ceil(10 * log10(s->peak * s->peak * + s->apsnr.n_pixels[i] * + 2)); + + err |= + vmaf_feature_collector_set_aggregate(feature_collector, + apsnr_name[i], + MIN(apsnr, max_apsnr)); + } + } + + return (err < 0) ? err : !err; +} + +static int close_fex_cuda(VmafFeatureExtractor *fex) +{ + PsnrStateCuda *s = fex->priv; + CudaFunctions *cu_f = fex->cu_state->f; + CHECK_CUDA(cu_f, cuStreamSynchronize(s->str)); + CHECK_CUDA(cu_f, cuStreamSynchronize(s->host_stream)); + CHECK_CUDA(cu_f, cuEventDestroy(s->finished)); + CHECK_CUDA(cu_f, cuEventDestroy(s->consumed)); + CHECK_CUDA(cu_f, cuEventDestroy(s->slot_done[0])); + CHECK_CUDA(cu_f, cuEventDestroy(s->slot_done[1])); + CHECK_CUDA(cu_f, cuStreamDestroy(s->str)); + CHECK_CUDA(cu_f, cuStreamDestroy(s->host_stream)); + + int ret = 0; + + if (s->sse) { + ret |= vmaf_cuda_buffer_free(fex->cu_state, s->sse); + free(s->sse); + } + if (s->sse_host) + ret |= vmaf_cuda_buffer_host_free(fex->cu_state, s->sse_host); + + if (s->write_score_parameters) + free(s->write_score_parameters); + + return ret; +} + +static const char *provided_features[] = { + "psnr_y", "psnr_cb", "psnr_cr", + NULL +}; + +VmafFeatureExtractor vmaf_fex_integer_psnr_cuda = { + .name = "psnr_cuda", + .options = options, + .init = init_fex_cuda, + .extract = extract_fex_cuda, + .flush = flush_fex_cuda, + .close = close_fex_cuda, + .priv_size = sizeof(PsnrStateCuda), + .provided_features = provided_features, + .flags = VMAF_FEATURE_EXTRACTOR_TEMPORAL | VMAF_FEATURE_EXTRACTOR_CUDA | + VMAF_FEATURE_EXTRACTOR_CUDA_CHROMA, +}; diff --git a/libvmaf/src/feature/cuda/integer_psnr_cuda.h b/libvmaf/src/feature/cuda/integer_psnr_cuda.h new file mode 100644 index 000000000..07a7db7c7 --- /dev/null +++ b/libvmaf/src/feature/cuda/integer_psnr_cuda.h @@ -0,0 +1,26 @@ +/** + * + * Copyright 2026 Bardie Høgh Joensen + * + * Licensed under the BSD+Patent License (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://opensource.org/licenses/BSDplusPatent + * + * 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. + * + */ + +#ifndef FEATURE_PSNR_CUDA_H_ +#define FEATURE_PSNR_CUDA_H_ + +#include +#include "common.h" + +extern const unsigned char psnr_ptx[]; +#endif /* _FEATURE_PSNR_CUDA_H_ */ diff --git a/libvmaf/src/feature/feature_extractor.c b/libvmaf/src/feature/feature_extractor.c index 1df405940..7cf7de9fe 100644 --- a/libvmaf/src/feature/feature_extractor.c +++ b/libvmaf/src/feature/feature_extractor.c @@ -54,6 +54,9 @@ extern VmafFeatureExtractor vmaf_fex_cambi; extern VmafFeatureExtractor vmaf_fex_integer_adm_cuda; extern VmafFeatureExtractor vmaf_fex_integer_vif_cuda; extern VmafFeatureExtractor vmaf_fex_integer_motion_cuda; +extern VmafFeatureExtractor vmaf_fex_integer_psnr_cuda; +extern VmafFeatureExtractor vmaf_fex_float_ssim_cuda; +extern VmafFeatureExtractor vmaf_fex_ciede_cuda; #endif extern VmafFeatureExtractor vmaf_fex_null; @@ -81,6 +84,9 @@ static VmafFeatureExtractor *feature_extractor_list[] = { &vmaf_fex_integer_adm_cuda, &vmaf_fex_integer_vif_cuda, &vmaf_fex_integer_motion_cuda, + &vmaf_fex_integer_psnr_cuda, + &vmaf_fex_float_ssim_cuda, + &vmaf_fex_ciede_cuda, #endif &vmaf_fex_null, NULL diff --git a/libvmaf/src/feature/feature_extractor.h b/libvmaf/src/feature/feature_extractor.h index ecd25f1a6..fb6f69bc2 100644 --- a/libvmaf/src/feature/feature_extractor.h +++ b/libvmaf/src/feature/feature_extractor.h @@ -39,6 +39,7 @@ enum VmafFeatureExtractorFlags { VMAF_FEATURE_EXTRACTOR_CUDA = 1 << 1, VMAF_FEATURE_FRAME_SYNC = 1 << 2, VMAF_FEATURE_EXTRACTOR_PREV_REF = 1 << 3, + VMAF_FEATURE_EXTRACTOR_CUDA_CHROMA = 1 << 4, }; typedef struct VmafFeatureExtractor { diff --git a/libvmaf/src/libvmaf.c b/libvmaf/src/libvmaf.c index f914ddaa2..ea24779af 100644 --- a/libvmaf/src/libvmaf.c +++ b/libvmaf/src/libvmaf.c @@ -649,6 +649,12 @@ static int flush_context_threaded(VmafContext *vmaf) for (unsigned i = 0; i < rfe.cnt; i++) { if (!(rfe.fex_ctx[i]->fex->flags & VMAF_FEATURE_EXTRACTOR_TEMPORAL)) continue; + /* CUDA feature extractors are flushed by the HAVE_CUDA block in + * flush_context(); flushing them here as well double-appends their + * final scores. */ + if (rfe.fex_ctx[i]->fex->flags & VMAF_FEATURE_EXTRACTOR_CUDA) + continue; + err |= vmaf_feature_extractor_context_flush(rfe.fex_ctx[i], vmaf->feature_collector); } @@ -734,6 +740,21 @@ enum { HW_FLAG_DEVICE = 1 << 1, }; +static uint8_t rfe_cuda_plane_mask(RegisteredFeatureExtractors *rfe) +{ + // existing CUDA feature extractors are luma-only, so only upload the + // chroma planes when a registered extractor declares it reads them + uint8_t mask = 0x1; + for (unsigned i = 0; i < rfe->cnt; i++) { + const uint64_t flags = rfe->fex_ctx[i]->fex->flags; + if ((flags & VMAF_FEATURE_EXTRACTOR_CUDA) && + (flags & VMAF_FEATURE_EXTRACTOR_CUDA_CHROMA)) + mask |= 0x6; + } + + return mask; +} + static int translate_picture_host(VmafContext *vmaf, VmafPicture *pic, VmafPicture *pic_device, unsigned hw_flags) { @@ -748,7 +769,8 @@ static int translate_picture_host(VmafContext *vmaf, VmafPicture *pic, if (!vmaf->cuda.state.ctx) return -EINVAL; err |= vmaf_ring_buffer_fetch_next_picture(vmaf->cuda.ring_buffer, pic_device); - err |= vmaf_cuda_picture_upload_async(pic_device, pic, 0x1); + err |= vmaf_cuda_picture_upload_async(pic_device, pic, + rfe_cuda_plane_mask(&vmaf->registered_feature_extractors)); if (err) { vmaf_log(VMAF_LOG_LEVEL_ERROR, "problem moving host pic into cuda device buffer\n"); @@ -778,7 +800,9 @@ static int translate_picture_device(VmafContext *vmaf, VmafPicture *pic, return err; } - err = vmaf_cuda_picture_download_async(pic, pic_host, 0x1); + // host pictures always carry every plane, so CPU feature extractors + // that read chroma (psnr, ciede, ...) expect all of them here + err = vmaf_cuda_picture_download_async(pic, pic_host, 0x7); if (err) { vmaf_log(VMAF_LOG_LEVEL_ERROR, "problem moving cuda pic into host buffer\n"); diff --git a/libvmaf/src/meson.build b/libvmaf/src/meson.build index 6900c572d..dc92680d7 100644 --- a/libvmaf/src/meson.build +++ b/libvmaf/src/meson.build @@ -354,6 +354,9 @@ if is_cuda_enabled 'adm_decouple' : [feature_src_dir + 'cuda/integer_adm/adm_decouple.cu'], 'filter1d' : [feature_src_dir + 'cuda/integer_vif/filter1d.cu'], 'motion_score' : [feature_src_dir + 'cuda/integer_motion/motion_score.cu'], + 'psnr' : [feature_src_dir + 'cuda/integer_psnr/psnr.cu'], + 'ssim' : [feature_src_dir + 'cuda/float_ssim/ssim.cu'], + 'ciede' : [feature_src_dir + 'cuda/ciede/ciede.cu'], } message(cuda_cu_sources) cuda_sources = [ @@ -541,6 +544,9 @@ if is_cuda_enabled feature_src_dir + 'cuda/integer_adm_cuda.c', feature_src_dir + 'cuda/integer_vif_cuda.c', feature_src_dir + 'cuda/integer_motion_cuda.c', + feature_src_dir + 'cuda/integer_psnr_cuda.c', + feature_src_dir + 'cuda/float_ssim_cuda.c', + feature_src_dir + 'cuda/ciede_cuda.c', ] endif diff --git a/libvmaf/test/meson.build b/libvmaf/test/meson.build index 17f1c7e3f..d1879c8e5 100644 --- a/libvmaf/test/meson.build +++ b/libvmaf/test/meson.build @@ -195,8 +195,26 @@ test_cuda_pic_preallocation = executable('test_cuda_pic_preallocation', c_args: ['-DHAVE_CUDA=1'] ) +test_cuda_psnr_ssim_parity = executable('test_cuda_psnr_ssim_parity', + ['test.c', 'test_cuda_psnr_ssim_parity.c'], + include_directories : [libvmaf_inc, test_inc], + link_with : get_option('default_library') == 'both' ? libvmaf.get_static_lib() : libvmaf, + dependencies: cuda_dependency, + c_args: ['-DHAVE_CUDA=1'] +) + +test_cuda_ciede_parity = executable('test_cuda_ciede_parity', + ['test.c', 'test_cuda_ciede_parity.c'], + include_directories : [libvmaf_inc, test_inc], + link_with : get_option('default_library') == 'both' ? libvmaf.get_static_lib() : libvmaf, + dependencies: cuda_dependency, + c_args: ['-DHAVE_CUDA=1'] +) + test('test_ring_buffer', test_ring_buffer) test('test_cuda_pic_preallocation', test_cuda_pic_preallocation) +test('test_cuda_psnr_ssim_parity', test_cuda_psnr_ssim_parity) +test('test_cuda_ciede_parity', test_cuda_ciede_parity) endif test_pic_preallocation = executable('test_pic_preallocation', diff --git a/libvmaf/test/test_cuda_ciede_parity.c b/libvmaf/test/test_cuda_ciede_parity.c new file mode 100644 index 000000000..038aa0882 --- /dev/null +++ b/libvmaf/test/test_cuda_ciede_parity.c @@ -0,0 +1,200 @@ +/** + * + * Copyright 2026 Bardie Høgh Joensen + * + * Licensed under the BSD+Patent License (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://opensource.org/licenses/BSDplusPatent + * + * 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. + * + */ + +/* + * CPU vs CUDA parity test for the ciede_cuda feature extractor. Unlike the + * psnr_cuda/ssim_cuda parity test this asserts a small tolerance rather than + * bit-exactness: ciede is dominated by libm transcendentals, which differ + * between glibc and CUDA in the low bits, and the device math runs in + * float32 (the CPU reference truncates every intermediate to float anyway). + * + * Uses YUV420P and YUV422P input so both axes of the fused nearest-neighbor + * chroma upsampling are exercised independently against the CPU's + * scale_chroma_planes. Also checks the identical-frame case, where both + * implementations must return +inf. + * + * Exits with meson's SKIP code (77) when no CUDA device is available so CI + * without a GPU reports the test as skipped. + */ + +#include +#include +#include +#include + +#include "test.h" + +#include "libvmaf/libvmaf.h" +#include "libvmaf/libvmaf_cuda.h" +#include "libvmaf/picture.h" + +#define N_FRAMES 4 +#define CIEDE_EPS 1e-3 + +static uint32_t lcg_state; + +static uint32_t lcg_next(void) +{ + lcg_state = lcg_state * 1664525u + 1013904223u; + return lcg_state >> 16; +} + +// frame N_FRAMES-1 is generated identical (ref == dist) to exercise the +// de00_sum == 0 -> +inf path +static void fill_pictures(VmafPicture *ref, VmafPicture *dist, unsigned bpc, + unsigned index) +{ + const unsigned peak = (1 << bpc) - 1; + const int identical = index == N_FRAMES - 1; + lcg_state = 54321u + index * 7919u; + + for (unsigned p = 0; p < 3; p++) { + for (unsigned i = 0; i < ref->h[p]; i++) { + for (unsigned j = 0; j < ref->w[p]; j++) { + const int v = (i + j + lcg_next()) % (peak + 1); + const int noise = identical ? 0 : (int)(lcg_next() % 31) - 15; + int vd = v + noise; + if (vd < 0) vd = 0; + if (vd > (int)peak) vd = peak; + if (bpc == 8) { + ((uint8_t*)ref->data[p])[i * ref->stride[p] + j] = v; + ((uint8_t*)dist->data[p])[i * dist->stride[p] + j] = vd; + } else { + ((uint16_t*)ref->data[p])[i * (ref->stride[p] / 2) + j] = v; + ((uint16_t*)dist->data[p])[i * (dist->stride[p] / 2) + j] = vd; + } + } + } + } +} + +// returns 0 on success, 1 when CUDA is unavailable (caller should skip) +static int run_pass(int use_cuda, enum VmafPixelFormat pix_fmt, unsigned bpc, + unsigned w, unsigned h, double scores[N_FRAMES], char **fail) +{ + int err = 0; + *fail = NULL; + + VmafConfiguration cfg = { + .log_level = VMAF_LOG_LEVEL_ERROR, + }; + + VmafContext *vmaf; + err = vmaf_init(&vmaf, cfg); + if (err) { *fail = "problem during vmaf_init"; return 0; } + + if (use_cuda) { + VmafCudaState *cu_state; + VmafCudaConfiguration cuda_cfg = { 0 }; + err = vmaf_cuda_state_init(&cu_state, cuda_cfg); + if (err) { + vmaf_close(vmaf); + return 1; // no CUDA device, skip + } + err = vmaf_cuda_import_state(vmaf, cu_state); + if (err) { *fail = "problem during vmaf_cuda_import_state"; return 0; } + } + + err = vmaf_use_feature(vmaf, use_cuda ? "ciede_cuda" : "ciede", NULL); + if (err) { *fail = "problem during vmaf_use_feature"; return 0; } + + for (unsigned i = 0; i < N_FRAMES; i++) { + VmafPicture ref, dist; + err = vmaf_picture_alloc(&ref, pix_fmt, bpc, w, h); + err |= vmaf_picture_alloc(&dist, pix_fmt, bpc, w, h); + if (err) { *fail = "problem during vmaf_picture_alloc"; return 0; } + fill_pictures(&ref, &dist, bpc, i); + err = vmaf_read_pictures(vmaf, &ref, &dist, i); + if (err) { *fail = "problem during vmaf_read_pictures"; return 0; } + } + + err = vmaf_read_pictures(vmaf, NULL, NULL, 0); + if (err) { *fail = "problem during vmaf_read_pictures flush"; return 0; } + + for (unsigned i = 0; i < N_FRAMES; i++) { + err = vmaf_feature_score_at_index(vmaf, "ciede2000", &scores[i], i); + if (err) { *fail = "problem during vmaf_feature_score_at_index"; return 0; } + } + + err = vmaf_close(vmaf); + if (err) { *fail = "problem during vmaf_close"; return 0; } + + return 0; +} + +static char *parity(enum VmafPixelFormat pix_fmt, unsigned bpc, + unsigned w, unsigned h) +{ + double cpu[N_FRAMES], gpu[N_FRAMES]; + char *fail = NULL; + + run_pass(0, pix_fmt, bpc, w, h, cpu, &fail); + if (fail) return fail; + + if (run_pass(1, pix_fmt, bpc, w, h, gpu, &fail)) { + fprintf(stderr, "no CUDA device available, skipping\n"); + exit(77); + } + if (fail) return fail; + + // the identical last frame must be +inf on both sides + mu_assert("cpu identical-frame score must be +inf", + isinf(cpu[N_FRAMES - 1]) && cpu[N_FRAMES - 1] > 0); + mu_assert("cuda identical-frame score must be +inf", + isinf(gpu[N_FRAMES - 1]) && gpu[N_FRAMES - 1] > 0); + + for (unsigned i = 0; i < N_FRAMES - 1; i++) { + if (fabs(cpu[i] - gpu[i]) > CIEDE_EPS) { + fprintf(stderr, "mismatch format %d, %u bpc, %ux%u, frame %u: " + "cpu=%.9f gpu=%.9f\n", pix_fmt, bpc, w, h, i, + cpu[i], gpu[i]); + return "cpu/cuda ciede2000 score mismatch"; + } + } + + return NULL; +} + +static char *test_ciede_cuda_parity_8bpc(void) +{ + return parity(VMAF_PIX_FMT_YUV420P, 8, 768, 432); +} + +static char *test_ciede_cuda_parity_10bpc(void) +{ + return parity(VMAF_PIX_FMT_YUV420P, 10, 768, 432); +} + +static char *test_ciede_cuda_parity_422_8bpc(void) +{ + return parity(VMAF_PIX_FMT_YUV422P, 8, 64, 48); +} + +static char *test_ciede_cuda_parity_422_10bpc(void) +{ + return parity(VMAF_PIX_FMT_YUV422P, 10, 64, 48); +} + +char *run_tests() +{ + mu_run_test(test_ciede_cuda_parity_8bpc); + mu_run_test(test_ciede_cuda_parity_10bpc); + mu_run_test(test_ciede_cuda_parity_422_8bpc); + mu_run_test(test_ciede_cuda_parity_422_10bpc); + return NULL; +} diff --git a/libvmaf/test/test_cuda_psnr_ssim_parity.c b/libvmaf/test/test_cuda_psnr_ssim_parity.c new file mode 100644 index 000000000..36f9476ac --- /dev/null +++ b/libvmaf/test/test_cuda_psnr_ssim_parity.c @@ -0,0 +1,349 @@ +/** + * + * Copyright 2026 Bardie Høgh Joensen + * + * Licensed under the BSD+Patent License (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://opensource.org/licenses/BSDplusPatent + * + * 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. + * + */ + +/* + * CPU vs CUDA parity test for the psnr_cuda and ssim_cuda feature + * extractors. Runs the same deterministic synthetic frames through the CPU + * extractors (psnr, float_ssim) and the CUDA extractors (psnr_cuda, + * ssim_cuda), then asserts bit-exact per-frame scores and matching serialized + * APSNR aggregates. + * + * Exits with meson's SKIP code (77) when no CUDA device is available so CI + * without a GPU reports the test as skipped. + */ + +#include +#include +#include +#include + +#include "test.h" + +#include "libvmaf/libvmaf.h" +#include "libvmaf/libvmaf_cuda.h" +#include "libvmaf/picture.h" + +#define N_FRAMES 5 + +static const char *score_keys[] = { + "psnr_y", "psnr_cb", "psnr_cr", + "mse_y", "mse_cb", "mse_cr", + "float_ssim", "float_ssim_l", "float_ssim_c", "float_ssim_s", +}; +#define N_KEYS (sizeof(score_keys) / sizeof(score_keys[0])) + +static const char *apsnr_keys[] = { "apsnr_y", "apsnr_cb", "apsnr_cr" }; +#define N_APSNR (sizeof(apsnr_keys) / sizeof(apsnr_keys[0])) + +typedef struct ParityCase { + const char *name; + enum VmafPixelFormat pix_fmt; + unsigned bpc, w, h; + int scale; + bool enable_db, clip_db; + bool reduced_hbd_peak; + double min_sse; +} ParityCase; + +static uint32_t lcg_state; + +static uint32_t lcg_next(void) +{ + lcg_state = lcg_state * 1664525u + 1013904223u; + return lcg_state >> 16; +} + +static void fill_pictures(VmafPicture *ref, VmafPicture *dist, unsigned bpc, + unsigned index) +{ + const unsigned peak = (1 << bpc) - 1; + lcg_state = 12345u + index * 7919u; + + for (unsigned p = 0; p < 3; p++) { + if (bpc == 8) { + uint8_t *r = ref->data[p]; + uint8_t *d = dist->data[p]; + for (unsigned i = 0; i < ref->h[p]; i++) { + for (unsigned j = 0; j < ref->w[p]; j++) { + const int v = (i + j + lcg_next()) % (peak + 1); + const int noise = (int)(lcg_next() % 15) - 7; + int vd = v + noise; + if (vd < 0) vd = 0; + if (vd > (int)peak) vd = peak; + r[j] = v; + d[j] = vd; + } + r += ref->stride[p]; + d += dist->stride[p]; + } + } else { + uint16_t *r = ref->data[p]; + uint16_t *d = dist->data[p]; + for (unsigned i = 0; i < ref->h[p]; i++) { + for (unsigned j = 0; j < ref->w[p]; j++) { + const int v = (i + j + lcg_next()) % (peak + 1); + const int noise = (int)(lcg_next() % 61) - 30; + int vd = v + noise; + if (vd < 0) vd = 0; + if (vd > (int)peak) vd = peak; + r[j] = v; + d[j] = vd; + } + r += ref->stride[p] / 2; + d += dist->stride[p] / 2; + } + } + } +} + +static int read_apsnr(VmafContext *vmaf, double scores[N_APSNR]) +{ + char path[128]; + const int len = snprintf(path, sizeof(path), + "test_cuda_psnr_ssim_parity_%p.json", + (void *)vmaf); + if (len < 0 || (size_t)len >= sizeof(path)) + return -1; + + if (vmaf_write_output(vmaf, path, VMAF_OUTPUT_FORMAT_JSON)) + return -1; + + FILE *file = fopen(path, "r"); + if (!file) { + remove(path); + return -1; + } + + bool found[N_APSNR] = { false }; + char line[256]; + while (fgets(line, sizeof(line), file)) { + for (unsigned k = 0; k < N_APSNR; k++) { + if (!strstr(line, apsnr_keys[k])) + continue; + + char *value = strchr(line, ':'); + char *end; + if (!value) + continue; + scores[k] = strtod(value + 1, &end); + found[k] = end != value + 1; + } + } + + const int close_err = fclose(file); + const int remove_err = remove(path); + if (close_err || remove_err) + return -1; + + for (unsigned k = 0; k < N_APSNR; k++) { + if (!found[k]) + return -1; + } + + return 0; +} + +// returns 0 on success, 1 when CUDA is unavailable (caller should skip) +static int run_pass(int use_cuda, const ParityCase *test_case, + double scores[N_FRAMES][N_KEYS], + double apsnr[N_APSNR], char **fail) +{ + int err = 0; + *fail = NULL; + + VmafConfiguration cfg = { + .log_level = VMAF_LOG_LEVEL_ERROR, + .n_threads = use_cuda ? 2 : 0, // threads + CUDA: the double-flush path + }; + + VmafContext *vmaf; + err = vmaf_init(&vmaf, cfg); + if (err) { *fail = "problem during vmaf_init"; return 0; } + + if (use_cuda) { + VmafCudaState *cu_state; + VmafCudaConfiguration cuda_cfg = { 0 }; + err = vmaf_cuda_state_init(&cu_state, cuda_cfg); + if (err) { + vmaf_close(vmaf); + return 1; // no CUDA device, skip + } + err = vmaf_cuda_import_state(vmaf, cu_state); + if (err) { *fail = "problem during vmaf_cuda_import_state"; return 0; } + } + + // APSNR combined with n_threads regression-tests flush idempotency. MSE, + // reduced peak and min_sse exercise the remaining PSNR scoring options. + VmafFeatureDictionary *psnr_dict = NULL; + err = vmaf_feature_dictionary_set(&psnr_dict, "enable_mse", "true"); + err |= vmaf_feature_dictionary_set(&psnr_dict, "enable_apsnr", "true"); + if (test_case->reduced_hbd_peak) + err |= vmaf_feature_dictionary_set(&psnr_dict, "reduced_hbd_peak", "true"); + if (test_case->min_sse > 0.0) { + char min_sse[32]; + snprintf(min_sse, sizeof(min_sse), "%.17g", test_case->min_sse); + err |= vmaf_feature_dictionary_set(&psnr_dict, "min_sse", min_sse); + } + if (err) { *fail = "problem configuring psnr options"; return 0; } + + err = vmaf_use_feature(vmaf, use_cuda ? "psnr_cuda" : "psnr", psnr_dict); + if (err) { *fail = "problem during vmaf_use_feature psnr"; return 0; } + + VmafFeatureDictionary *ssim_dict = NULL; + char scale[16]; + snprintf(scale, sizeof(scale), "%d", test_case->scale); + err = vmaf_feature_dictionary_set(&ssim_dict, "enable_lcs", "true"); + err |= vmaf_feature_dictionary_set(&ssim_dict, "scale", scale); + if (test_case->enable_db) + err |= vmaf_feature_dictionary_set(&ssim_dict, "enable_db", "true"); + if (test_case->clip_db) + err |= vmaf_feature_dictionary_set(&ssim_dict, "clip_db", "true"); + if (err) { *fail = "problem configuring ssim options"; return 0; } + + err = vmaf_use_feature(vmaf, use_cuda ? "ssim_cuda" : "float_ssim", + ssim_dict); + if (err) { *fail = "problem during vmaf_use_feature ssim"; return 0; } + + for (unsigned i = 0; i < N_FRAMES; i++) { + VmafPicture ref, dist; + err = vmaf_picture_alloc(&ref, test_case->pix_fmt, test_case->bpc, + test_case->w, test_case->h); + err |= vmaf_picture_alloc(&dist, test_case->pix_fmt, test_case->bpc, + test_case->w, test_case->h); + if (err) { *fail = "problem during vmaf_picture_alloc"; return 0; } + fill_pictures(&ref, &dist, test_case->bpc, i); + err = vmaf_read_pictures(vmaf, &ref, &dist, i); + if (err) { *fail = "problem during vmaf_read_pictures"; return 0; } + } + + err = vmaf_read_pictures(vmaf, NULL, NULL, 0); + if (err) { *fail = "problem during vmaf_read_pictures flush"; return 0; } + + for (unsigned i = 0; i < N_FRAMES; i++) { + for (unsigned k = 0; k < N_KEYS; k++) { + err = vmaf_feature_score_at_index(vmaf, score_keys[k], + &scores[i][k], i); + if (err) { *fail = "problem during vmaf_feature_score_at_index"; return 0; } + } + } + + if (read_apsnr(vmaf, apsnr)) { + *fail = "problem reading APSNR aggregates"; + return 0; + } + + err = vmaf_close(vmaf); + if (err) { *fail = "problem during vmaf_close"; return 0; } + + return 0; +} + +static char *parity(const ParityCase *test_case) +{ + double cpu[N_FRAMES][N_KEYS], gpu[N_FRAMES][N_KEYS]; + double cpu_apsnr[N_APSNR], gpu_apsnr[N_APSNR]; + char *fail = NULL; + + run_pass(0, test_case, cpu, cpu_apsnr, &fail); + if (fail) return fail; + + if (run_pass(1, test_case, gpu, gpu_apsnr, &fail)) { + // meson exitcode protocol: 77 = SKIP, so CI without a GPU reports + // this as skipped rather than silently passing + fprintf(stderr, "no CUDA device available, skipping\n"); + exit(77); + } + if (fail) return fail; + + for (unsigned i = 0; i < N_FRAMES; i++) { + for (unsigned k = 0; k < N_KEYS; k++) { + if (memcmp(&cpu[i][k], &gpu[i][k], sizeof(cpu[i][k]))) { + fprintf(stderr, "mismatch %s, frame %u, %s: " + "cpu=%a gpu=%a\n", test_case->name, i, + score_keys[k], cpu[i][k], gpu[i][k]); + return "cpu/cuda score mismatch"; + } + } + } + + for (unsigned k = 0; k < N_APSNR; k++) { + if (memcmp(&cpu_apsnr[k], &gpu_apsnr[k], sizeof(cpu_apsnr[k]))) { + fprintf(stderr, "aggregate mismatch %s, %s: cpu=%a gpu=%a\n", + test_case->name, apsnr_keys[k], cpu_apsnr[k], + gpu_apsnr[k]); + return "cpu/cuda APSNR aggregate mismatch"; + } + } + + return NULL; +} + +static char *test_psnr_ssim_cuda_parity_420_8bpc_auto_scale(void) +{ + const ParityCase test_case = { + .name = "yuv420p 8 bpc 768x432 auto scale", + .pix_fmt = VMAF_PIX_FMT_YUV420P, + .bpc = 8, .w = 768, .h = 432, + }; + return parity(&test_case); +} + +static char *test_psnr_ssim_cuda_parity_420_10bpc_auto_scale(void) +{ + const ParityCase test_case = { + .name = "yuv420p 10 bpc 768x432 auto scale", + .pix_fmt = VMAF_PIX_FMT_YUV420P, + .bpc = 10, .w = 768, .h = 432, + }; + return parity(&test_case); +} + +static char *test_psnr_ssim_cuda_parity_444_12bpc_odd_scale_3(void) +{ + const ParityCase test_case = { + .name = "yuv444p 12 bpc 65x49 scale 3", + .pix_fmt = VMAF_PIX_FMT_YUV444P, + .bpc = 12, .w = 65, .h = 49, + .scale = 3, + }; + return parity(&test_case); +} + +static char *test_psnr_ssim_cuda_parity_422_16bpc_scale_1_options(void) +{ + const ParityCase test_case = { + .name = "yuv422p 16 bpc 64x48 scale 1 with options", + .pix_fmt = VMAF_PIX_FMT_YUV422P, + .bpc = 16, .w = 64, .h = 48, + .scale = 1, + .enable_db = true, + .clip_db = true, + .reduced_hbd_peak = true, + .min_sse = 1.0, + }; + return parity(&test_case); +} + +char *run_tests() +{ + mu_run_test(test_psnr_ssim_cuda_parity_420_8bpc_auto_scale); + mu_run_test(test_psnr_ssim_cuda_parity_420_10bpc_auto_scale); + mu_run_test(test_psnr_ssim_cuda_parity_444_12bpc_odd_scale_3); + mu_run_test(test_psnr_ssim_cuda_parity_422_16bpc_scale_1_options); + return NULL; +} diff --git a/libvmaf/test/test_feature_extractor.c b/libvmaf/test/test_feature_extractor.c index aac8a067a..b8b0de1ea 100644 --- a/libvmaf/test/test_feature_extractor.c +++ b/libvmaf/test/test_feature_extractor.c @@ -47,6 +47,18 @@ static char *test_get_feature_extractor_by_name_and_feature_name() "VMAF_integer_feature_adm2_score", flags); mu_assert("problem during vmaf_get_feature_extractor_by_feature_name", fex && !strcmp(fex->name, "adm_cuda")); + + fex = vmaf_get_feature_extractor_by_feature_name("psnr_y", flags); + mu_assert("problem during vmaf_get_feature_extractor_by_feature_name", + fex && !strcmp(fex->name, "psnr_cuda")); + + fex = vmaf_get_feature_extractor_by_feature_name("float_ssim", flags); + mu_assert("problem during vmaf_get_feature_extractor_by_feature_name", + fex && !strcmp(fex->name, "ssim_cuda")); + + fex = vmaf_get_feature_extractor_by_feature_name("ciede2000", flags); + mu_assert("problem during vmaf_get_feature_extractor_by_feature_name", + fex && !strcmp(fex->name, "ciede_cuda")); #endif return NULL;