ref: eaa8043a31be5304946cd33d8f3f9f72547bff4a
parent: 4bc6ae434207c854150c0d9dd014fb527da591c2
parent: c7b7011b9bb959db69b9c4a65d99c8057cebb6cb
author: Aℓex Converse <[email protected]>
date: Fri Aug 7 12:43:28 EDT 2015
Merge "Move VP9 SSIM metrics to vpx_dsp."
--- a/test/consistency_test.cc
+++ b/test/consistency_test.cc
@@ -23,11 +23,11 @@
#include "test/clear_system_state.h"
#include "test/register_state_check.h"
#include "test/util.h"
-#include "vp9/encoder/vp9_ssim.h"
+#include "vpx_dsp/ssim.h"
#include "vpx_mem/vpx_mem.h"
extern "C"
-double vp9_get_ssim_metrics(uint8_t *img1, int img1_pitch,
+double vpx_get_ssim_metrics(uint8_t *img1, int img1_pitch,
uint8_t *img2, int img2_pitch,
int width, int height,
Ssimv *sv2, Metrics *m,
@@ -144,7 +144,7 @@
double CheckConsistency(int frame) {EXPECT_LT(frame, 2)<< "Frame to check has to be less than 2.";
return
- vp9_get_ssim_metrics(source_data_[frame], source_stride_,
+ vpx_get_ssim_metrics(source_data_[frame], source_stride_,
reference_data_[frame], reference_stride_,
width_, height_, ssim_array_, &metrics_, 1);
}
--- a/vp9/common/vp9_rtcd_defs.pl
+++ b/vp9/common/vp9_rtcd_defs.pl
@@ -261,17 +261,6 @@
specialize qw/vp9_fdct8x8_quant sse2 ssse3 neon/;
}
-#
-# Structured Similarity (SSIM)
-#
-if (vpx_config("CONFIG_INTERNAL_STATS") eq "yes") {- add_proto qw/void vp9_ssim_parms_8x8/, "uint8_t *s, int sp, uint8_t *r, int rp, unsigned long *sum_s, unsigned long *sum_r, unsigned long *sum_sq_s, unsigned long *sum_sq_r, unsigned long *sum_sxr";
- specialize qw/vp9_ssim_parms_8x8/, "$sse2_x86_64";
-
- add_proto qw/void vp9_ssim_parms_16x16/, "uint8_t *s, int sp, uint8_t *r, int rp, unsigned long *sum_s, unsigned long *sum_r, unsigned long *sum_sq_s, unsigned long *sum_sq_r, unsigned long *sum_sxr";
- specialize qw/vp9_ssim_parms_16x16/, "$sse2_x86_64";
-}
-
# fdct functions
if (vpx_config("CONFIG_VP9_HIGHBITDEPTH") eq "yes") {@@ -329,14 +318,6 @@
add_proto qw/void vp9_highbd_quantize_fp_32x32/, "const tran_low_t *coeff_ptr, intptr_t n_coeffs, int skip_block, const int16_t *zbin_ptr, const int16_t *round_ptr, const int16_t *quant_ptr, const int16_t *quant_shift_ptr, tran_low_t *qcoeff_ptr, tran_low_t *dqcoeff_ptr, const int16_t *dequant_ptr, uint16_t *eob_ptr, const int16_t *scan, const int16_t *iscan";
specialize qw/vp9_highbd_quantize_fp_32x32/;
-
- #
- # Structured Similarity (SSIM)
- #
- if (vpx_config("CONFIG_INTERNAL_STATS") eq "yes") {- add_proto qw/void vp9_highbd_ssim_parms_8x8/, "uint16_t *s, int sp, uint16_t *r, int rp, uint32_t *sum_s, uint32_t *sum_r, uint32_t *sum_sq_s, uint32_t *sum_sq_r, uint32_t *sum_sxr";
- specialize qw/vp9_highbd_ssim_parms_8x8/;
- }
# fdct functions
add_proto qw/void vp9_highbd_fht4x4/, "const int16_t *input, tran_low_t *output, int stride, int tx_type";
--- a/vp9/encoder/vp9_encoder.c
+++ b/vp9/encoder/vp9_encoder.c
@@ -18,6 +18,9 @@
#include "./vpx_scale_rtcd.h"
#include "vpx/internal/vpx_psnr.h"
#include "vpx_dsp/vpx_filter.h"
+#if CONFIG_INTERNAL_STATS
+#include "vpx_dsp/ssim.h"
+#endif
#include "vpx_ports/mem.h"
#include "vpx_ports/vpx_timer.h"
#include "vpx_scale/vpx_scale.h"
@@ -51,9 +54,6 @@
#include "vp9/encoder/vp9_segmentation.h"
#include "vp9/encoder/vp9_skin_detection.h"
#include "vp9/encoder/vp9_speed_features.h"
-#if CONFIG_INTERNAL_STATS
-#include "vp9/encoder/vp9_ssim.h"
-#endif
#include "vp9/encoder/vp9_svc_layercontext.h"
#include "vp9/encoder/vp9_temporal_filter.h"
@@ -4416,13 +4416,13 @@
#if CONFIG_VP9_HIGHBITDEPTH
if (cm->use_highbitdepth) {- frame_ssim2 = vp9_highbd_calc_ssim(orig, recon, &weight,
+ frame_ssim2 = vpx_highbd_calc_ssim(orig, recon, &weight,
(int)cm->bit_depth);
} else {- frame_ssim2 = vp9_calc_ssim(orig, recon, &weight);
+ frame_ssim2 = vpx_calc_ssim(orig, recon, &weight);
}
#else
- frame_ssim2 = vp9_calc_ssim(orig, recon, &weight);
+ frame_ssim2 = vpx_calc_ssim(orig, recon, &weight);
#endif // CONFIG_VP9_HIGHBITDEPTH
cpi->worst_ssim= MIN(cpi->worst_ssim, frame_ssim2);
@@ -4431,13 +4431,13 @@
#if CONFIG_VP9_HIGHBITDEPTH
if (cm->use_highbitdepth) {- frame_ssim2 = vp9_highbd_calc_ssim(
+ frame_ssim2 = vpx_highbd_calc_ssim(
orig, &cm->post_proc_buffer, &weight, (int)cm->bit_depth);
} else {- frame_ssim2 = vp9_calc_ssim(orig, &cm->post_proc_buffer, &weight);
+ frame_ssim2 = vpx_calc_ssim(orig, &cm->post_proc_buffer, &weight);
}
#else
- frame_ssim2 = vp9_calc_ssim(orig, &cm->post_proc_buffer, &weight);
+ frame_ssim2 = vpx_calc_ssim(orig, &cm->post_proc_buffer, &weight);
#endif // CONFIG_VP9_HIGHBITDEPTH
cpi->summedp_quality += frame_ssim2 * weight;
@@ -4472,7 +4472,7 @@
if (!cm->use_highbitdepth)
#endif
{- double this_inconsistency = vp9_get_ssim_metrics(
+ double this_inconsistency = vpx_get_ssim_metrics(
cpi->Source->y_buffer, cpi->Source->y_stride,
cm->frame_to_show->y_buffer, cm->frame_to_show->y_stride,
cpi->Source->y_width, cpi->Source->y_height, cpi->ssim_vars,
@@ -4492,14 +4492,14 @@
double y, u, v, frame_all;
#if CONFIG_VP9_HIGHBITDEPTH
if (cm->use_highbitdepth) {- frame_all = vp9_highbd_calc_ssimg(cpi->Source, cm->frame_to_show, &y,
+ frame_all = vpx_highbd_calc_ssimg(cpi->Source, cm->frame_to_show, &y,
&u, &v, (int)cm->bit_depth);
} else {- frame_all = vp9_calc_ssimg(cpi->Source, cm->frame_to_show, &y, &u,
+ frame_all = vpx_calc_ssimg(cpi->Source, cm->frame_to_show, &y, &u,
&v);
}
#else
- frame_all = vp9_calc_ssimg(cpi->Source, cm->frame_to_show, &y, &u, &v);
+ frame_all = vpx_calc_ssimg(cpi->Source, cm->frame_to_show, &y, &u, &v);
#endif // CONFIG_VP9_HIGHBITDEPTH
adjust_image_stat(y, u, v, frame_all, &cpi->ssimg);
}
@@ -4508,7 +4508,7 @@
#endif
{double y, u, v, frame_all;
- frame_all = vp9_calc_fastssim(cpi->Source, cm->frame_to_show, &y, &u,
+ frame_all = vpx_calc_fastssim(cpi->Source, cm->frame_to_show, &y, &u,
&v);
adjust_image_stat(y, u, v, frame_all, &cpi->fastssim);
/* TODO(JBB): add 10/12 bit support */
@@ -4518,7 +4518,7 @@
#endif
{double y, u, v, frame_all;
- frame_all = vp9_psnrhvs(cpi->Source, cm->frame_to_show, &y, &u, &v);
+ frame_all = vpx_psnrhvs(cpi->Source, cm->frame_to_show, &y, &u, &v);
adjust_image_stat(y, u, v, frame_all, &cpi->psnrhvs);
}
}
--- a/vp9/encoder/vp9_encoder.h
+++ b/vp9/encoder/vp9_encoder.h
@@ -16,6 +16,10 @@
#include "./vpx_config.h"
#include "vpx/internal/vpx_codec_internal.h"
#include "vpx/vp8cx.h"
+#if CONFIG_INTERNAL_STATS
+#include "vpx_dsp/ssim.h"
+#endif
+#include "vpx_dsp/variance.h"
#include "vpx_util/vpx_thread.h"
#include "vp9/common/vp9_alloccommon.h"
@@ -34,13 +38,9 @@
#include "vp9/encoder/vp9_quantize.h"
#include "vp9/encoder/vp9_ratectrl.h"
#include "vp9/encoder/vp9_rd.h"
-#if CONFIG_INTERNAL_STATS
-#include "vp9/encoder/vp9_ssim.h"
-#endif
#include "vp9/encoder/vp9_speed_features.h"
#include "vp9/encoder/vp9_svc_layercontext.h"
#include "vp9/encoder/vp9_tokenize.h"
-#include "vpx_dsp/variance.h"
#if CONFIG_VP9_TEMPORAL_DENOISING
#include "vp9/encoder/vp9_denoiser.h"
--- a/vp9/encoder/vp9_fastssim.c
+++ /dev/null
@@ -1,465 +1,0 @@
-/*
- * Copyright (c) 2010 The WebM project authors. All Rights Reserved.
- *
- * Use of this source code is governed by a BSD-style license
- * that can be found in the LICENSE file in the root of the source
- * tree. An additional intellectual property rights grant can be found
- * in the file PATENTS. All contributing project authors may
- * be found in the AUTHORS file in the root of the source tree.
- *
- * This code was originally written by: Nathan E. Egge, at the Daala
- * project.
- */
-#include <math.h>
-#include <string.h>
-#include "./vpx_config.h"
-#include "./vp9_rtcd.h"
-#include "vp9/encoder/vp9_ssim.h"
-/* TODO(jbb): High bit depth version of this code needed */
-typedef struct fs_level fs_level;
-typedef struct fs_ctx fs_ctx;
-
-#define SSIM_C1 (255 * 255 * 0.01 * 0.01)
-#define SSIM_C2 (255 * 255 * 0.03 * 0.03)
-
-#define FS_MINI(_a, _b) ((_a) < (_b) ? (_a) : (_b))
-#define FS_MAXI(_a, _b) ((_a) > (_b) ? (_a) : (_b))
-
-struct fs_level {- uint16_t *im1;
- uint16_t *im2;
- double *ssim;
- int w;
- int h;
-};
-
-struct fs_ctx {- fs_level *level;
- int nlevels;
- unsigned *col_buf;
-};
-
-static void fs_ctx_init(fs_ctx *_ctx, int _w, int _h, int _nlevels) {- unsigned char *data;
- size_t data_size;
- int lw;
- int lh;
- int l;
- lw = (_w + 1) >> 1;
- lh = (_h + 1) >> 1;
- data_size = _nlevels * sizeof(fs_level)
- + 2 * (lw + 8) * 8 * sizeof(*_ctx->col_buf);
- for (l = 0; l < _nlevels; l++) {- size_t im_size;
- size_t level_size;
- im_size = lw * (size_t) lh;
- level_size = 2 * im_size * sizeof(*_ctx->level[l].im1);
- level_size += sizeof(*_ctx->level[l].ssim) - 1;
- level_size /= sizeof(*_ctx->level[l].ssim);
- level_size += im_size;
- level_size *= sizeof(*_ctx->level[l].ssim);
- data_size += level_size;
- lw = (lw + 1) >> 1;
- lh = (lh + 1) >> 1;
- }
- data = (unsigned char *) malloc(data_size);
- _ctx->level = (fs_level *) data;
- _ctx->nlevels = _nlevels;
- data += _nlevels * sizeof(*_ctx->level);
- lw = (_w + 1) >> 1;
- lh = (_h + 1) >> 1;
- for (l = 0; l < _nlevels; l++) {- size_t im_size;
- size_t level_size;
- _ctx->level[l].w = lw;
- _ctx->level[l].h = lh;
- im_size = lw * (size_t) lh;
- level_size = 2 * im_size * sizeof(*_ctx->level[l].im1);
- level_size += sizeof(*_ctx->level[l].ssim) - 1;
- level_size /= sizeof(*_ctx->level[l].ssim);
- level_size *= sizeof(*_ctx->level[l].ssim);
- _ctx->level[l].im1 = (uint16_t *) data;
- _ctx->level[l].im2 = _ctx->level[l].im1 + im_size;
- data += level_size;
- _ctx->level[l].ssim = (double *) data;
- data += im_size * sizeof(*_ctx->level[l].ssim);
- lw = (lw + 1) >> 1;
- lh = (lh + 1) >> 1;
- }
- _ctx->col_buf = (unsigned *) data;
-}
-
-static void fs_ctx_clear(fs_ctx *_ctx) {- free(_ctx->level);
-}
-
-static void fs_downsample_level(fs_ctx *_ctx, int _l) {- const uint16_t *src1;
- const uint16_t *src2;
- uint16_t *dst1;
- uint16_t *dst2;
- int w2;
- int h2;
- int w;
- int h;
- int i;
- int j;
- w = _ctx->level[_l].w;
- h = _ctx->level[_l].h;
- dst1 = _ctx->level[_l].im1;
- dst2 = _ctx->level[_l].im2;
- w2 = _ctx->level[_l - 1].w;
- h2 = _ctx->level[_l - 1].h;
- src1 = _ctx->level[_l - 1].im1;
- src2 = _ctx->level[_l - 1].im2;
- for (j = 0; j < h; j++) {- int j0offs;
- int j1offs;
- j0offs = 2 * j * w2;
- j1offs = FS_MINI(2 * j + 1, h2) * w2;
- for (i = 0; i < w; i++) {- int i0;
- int i1;
- i0 = 2 * i;
- i1 = FS_MINI(i0 + 1, w2);
- dst1[j * w + i] = src1[j0offs + i0] + src1[j0offs + i1]
- + src1[j1offs + i0] + src1[j1offs + i1];
- dst2[j * w + i] = src2[j0offs + i0] + src2[j0offs + i1]
- + src2[j1offs + i0] + src2[j1offs + i1];
- }
- }
-}
-
-static void fs_downsample_level0(fs_ctx *_ctx, const unsigned char *_src1,
- int _s1ystride, const unsigned char *_src2,
- int _s2ystride, int _w, int _h) {- uint16_t *dst1;
- uint16_t *dst2;
- int w;
- int h;
- int i;
- int j;
- w = _ctx->level[0].w;
- h = _ctx->level[0].h;
- dst1 = _ctx->level[0].im1;
- dst2 = _ctx->level[0].im2;
- for (j = 0; j < h; j++) {- int j0;
- int j1;
- j0 = 2 * j;
- j1 = FS_MINI(j0 + 1, _h);
- for (i = 0; i < w; i++) {- int i0;
- int i1;
- i0 = 2 * i;
- i1 = FS_MINI(i0 + 1, _w);
- dst1[j * w + i] = _src1[j0 * _s1ystride + i0]
- + _src1[j0 * _s1ystride + i1] + _src1[j1 * _s1ystride + i0]
- + _src1[j1 * _s1ystride + i1];
- dst2[j * w + i] = _src2[j0 * _s2ystride + i0]
- + _src2[j0 * _s2ystride + i1] + _src2[j1 * _s2ystride + i0]
- + _src2[j1 * _s2ystride + i1];
- }
- }
-}
-
-static void fs_apply_luminance(fs_ctx *_ctx, int _l) {- unsigned *col_sums_x;
- unsigned *col_sums_y;
- uint16_t *im1;
- uint16_t *im2;
- double *ssim;
- double c1;
- int w;
- int h;
- int j0offs;
- int j1offs;
- int i;
- int j;
- w = _ctx->level[_l].w;
- h = _ctx->level[_l].h;
- col_sums_x = _ctx->col_buf;
- col_sums_y = col_sums_x + w;
- im1 = _ctx->level[_l].im1;
- im2 = _ctx->level[_l].im2;
- for (i = 0; i < w; i++)
- col_sums_x[i] = 5 * im1[i];
- for (i = 0; i < w; i++)
- col_sums_y[i] = 5 * im2[i];
- for (j = 1; j < 4; j++) {- j1offs = FS_MINI(j, h - 1) * w;
- for (i = 0; i < w; i++)
- col_sums_x[i] += im1[j1offs + i];
- for (i = 0; i < w; i++)
- col_sums_y[i] += im2[j1offs + i];
- }
- ssim = _ctx->level[_l].ssim;
- c1 = (double) (SSIM_C1 * 4096 * (1 << 4 * _l));
- for (j = 0; j < h; j++) {- unsigned mux;
- unsigned muy;
- int i0;
- int i1;
- mux = 5 * col_sums_x[0];
- muy = 5 * col_sums_y[0];
- for (i = 1; i < 4; i++) {- i1 = FS_MINI(i, w - 1);
- mux += col_sums_x[i1];
- muy += col_sums_y[i1];
- }
- for (i = 0; i < w; i++) {- ssim[j * w + i] *= (2 * mux * (double) muy + c1)
- / (mux * (double) mux + muy * (double) muy + c1);
- if (i + 1 < w) {- i0 = FS_MAXI(0, i - 4);
- i1 = FS_MINI(i + 4, w - 1);
- mux += col_sums_x[i1] - col_sums_x[i0];
- muy += col_sums_x[i1] - col_sums_x[i0];
- }
- }
- if (j + 1 < h) {- j0offs = FS_MAXI(0, j - 4) * w;
- for (i = 0; i < w; i++)
- col_sums_x[i] -= im1[j0offs + i];
- for (i = 0; i < w; i++)
- col_sums_y[i] -= im2[j0offs + i];
- j1offs = FS_MINI(j + 4, h - 1) * w;
- for (i = 0; i < w; i++)
- col_sums_x[i] += im1[j1offs + i];
- for (i = 0; i < w; i++)
- col_sums_y[i] += im2[j1offs + i];
- }
- }
-}
-
-#define FS_COL_SET(_col, _joffs, _ioffs) \
- do { \- unsigned gx; \
- unsigned gy; \
- gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
- gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
- col_sums_gx2[(_col)] = gx * (double)gx; \
- col_sums_gy2[(_col)] = gy * (double)gy; \
- col_sums_gxgy[(_col)] = gx * (double)gy; \
- } \
- while (0)
-
-#define FS_COL_ADD(_col, _joffs, _ioffs) \
- do { \- unsigned gx; \
- unsigned gy; \
- gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
- gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
- col_sums_gx2[(_col)] += gx * (double)gx; \
- col_sums_gy2[(_col)] += gy * (double)gy; \
- col_sums_gxgy[(_col)] += gx * (double)gy; \
- } \
- while (0)
-
-#define FS_COL_SUB(_col, _joffs, _ioffs) \
- do { \- unsigned gx; \
- unsigned gy; \
- gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
- gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
- col_sums_gx2[(_col)] -= gx * (double)gx; \
- col_sums_gy2[(_col)] -= gy * (double)gy; \
- col_sums_gxgy[(_col)] -= gx * (double)gy; \
- } \
- while (0)
-
-#define FS_COL_COPY(_col1, _col2) \
- do { \- col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)]; \
- col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)]; \
- col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)]; \
- } \
- while (0)
-
-#define FS_COL_HALVE(_col1, _col2) \
- do { \- col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)] * 0.5; \
- col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)] * 0.5; \
- col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)] * 0.5; \
- } \
- while (0)
-
-#define FS_COL_DOUBLE(_col1, _col2) \
- do { \- col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)] * 2; \
- col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)] * 2; \
- col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)] * 2; \
- } \
- while (0)
-
-static void fs_calc_structure(fs_ctx *_ctx, int _l) {- uint16_t *im1;
- uint16_t *im2;
- unsigned *gx_buf;
- unsigned *gy_buf;
- double *ssim;
- double col_sums_gx2[8];
- double col_sums_gy2[8];
- double col_sums_gxgy[8];
- double c2;
- int stride;
- int w;
- int h;
- int i;
- int j;
- w = _ctx->level[_l].w;
- h = _ctx->level[_l].h;
- im1 = _ctx->level[_l].im1;
- im2 = _ctx->level[_l].im2;
- ssim = _ctx->level[_l].ssim;
- gx_buf = _ctx->col_buf;
- stride = w + 8;
- gy_buf = gx_buf + 8 * stride;
- memset(gx_buf, 0, 2 * 8 * stride * sizeof(*gx_buf));
- c2 = SSIM_C2 * (1 << 4 * _l) * 16 * 104;
- for (j = 0; j < h + 4; j++) {- if (j < h - 1) {- for (i = 0; i < w - 1; i++) {- unsigned g1;
- unsigned g2;
- unsigned gx;
- unsigned gy;
- g1 = abs(im1[(j + 1) * w + i + 1] - im1[j * w + i]);
- g2 = abs(im1[(j + 1) * w + i] - im1[j * w + i + 1]);
- gx = 4 * FS_MAXI(g1, g2) + FS_MINI(g1, g2);
- g1 = abs(im2[(j + 1) * w + i + 1] - im2[j * w + i]);
- g2 = abs(im2[(j + 1) * w + i] - im2[j * w + i + 1]);
- gy = 4 * FS_MAXI(g1, g2) + FS_MINI(g1, g2);
- gx_buf[(j & 7) * stride + i + 4] = gx;
- gy_buf[(j & 7) * stride + i + 4] = gy;
- }
- } else {- memset(gx_buf + (j & 7) * stride, 0, stride * sizeof(*gx_buf));
- memset(gy_buf + (j & 7) * stride, 0, stride * sizeof(*gy_buf));
- }
- if (j >= 4) {- int k;
- col_sums_gx2[3] = col_sums_gx2[2] = col_sums_gx2[1] = col_sums_gx2[0] = 0;
- col_sums_gy2[3] = col_sums_gy2[2] = col_sums_gy2[1] = col_sums_gy2[0] = 0;
- col_sums_gxgy[3] = col_sums_gxgy[2] = col_sums_gxgy[1] =
- col_sums_gxgy[0] = 0;
- for (i = 4; i < 8; i++) {- FS_COL_SET(i, -1, 0);
- FS_COL_ADD(i, 0, 0);
- for (k = 1; k < 8 - i; k++) {- FS_COL_DOUBLE(i, i);
- FS_COL_ADD(i, -k - 1, 0);
- FS_COL_ADD(i, k, 0);
- }
- }
- for (i = 0; i < w; i++) {- double mugx2;
- double mugy2;
- double mugxgy;
- mugx2 = col_sums_gx2[0];
- for (k = 1; k < 8; k++)
- mugx2 += col_sums_gx2[k];
- mugy2 = col_sums_gy2[0];
- for (k = 1; k < 8; k++)
- mugy2 += col_sums_gy2[k];
- mugxgy = col_sums_gxgy[0];
- for (k = 1; k < 8; k++)
- mugxgy += col_sums_gxgy[k];
- ssim[(j - 4) * w + i] = (2 * mugxgy + c2) / (mugx2 + mugy2 + c2);
- if (i + 1 < w) {- FS_COL_SET(0, -1, 1);
- FS_COL_ADD(0, 0, 1);
- FS_COL_SUB(2, -3, 2);
- FS_COL_SUB(2, 2, 2);
- FS_COL_HALVE(1, 2);
- FS_COL_SUB(3, -4, 3);
- FS_COL_SUB(3, 3, 3);
- FS_COL_HALVE(2, 3);
- FS_COL_COPY(3, 4);
- FS_COL_DOUBLE(4, 5);
- FS_COL_ADD(4, -4, 5);
- FS_COL_ADD(4, 3, 5);
- FS_COL_DOUBLE(5, 6);
- FS_COL_ADD(5, -3, 6);
- FS_COL_ADD(5, 2, 6);
- FS_COL_DOUBLE(6, 7);
- FS_COL_ADD(6, -2, 7);
- FS_COL_ADD(6, 1, 7);
- FS_COL_SET(7, -1, 8);
- FS_COL_ADD(7, 0, 8);
- }
- }
- }
- }
-}
-
-#define FS_NLEVELS (4)
-
-/*These weights were derived from the default weights found in Wang's original
- Matlab implementation: {0.0448, 0.2856, 0.2363, 0.1333}.- We drop the finest scale and renormalize the rest to sum to 1.*/
-
-static const double FS_WEIGHTS[FS_NLEVELS] = {0.2989654541015625,- 0.3141326904296875, 0.2473602294921875, 0.1395416259765625};
-
-static double fs_average(fs_ctx *_ctx, int _l) {- double *ssim;
- double ret;
- int w;
- int h;
- int i;
- int j;
- w = _ctx->level[_l].w;
- h = _ctx->level[_l].h;
- ssim = _ctx->level[_l].ssim;
- ret = 0;
- for (j = 0; j < h; j++)
- for (i = 0; i < w; i++)
- ret += ssim[j * w + i];
- return pow(ret / (w * h), FS_WEIGHTS[_l]);
-}
-
-static double calc_ssim(const unsigned char *_src, int _systride,
- const unsigned char *_dst, int _dystride, int _w, int _h) {- fs_ctx ctx;
- double ret;
- int l;
- ret = 1;
- fs_ctx_init(&ctx, _w, _h, FS_NLEVELS);
- fs_downsample_level0(&ctx, _src, _systride, _dst, _dystride, _w, _h);
- for (l = 0; l < FS_NLEVELS - 1; l++) {- fs_calc_structure(&ctx, l);
- ret *= fs_average(&ctx, l);
- fs_downsample_level(&ctx, l + 1);
- }
- fs_calc_structure(&ctx, l);
- fs_apply_luminance(&ctx, l);
- ret *= fs_average(&ctx, l);
- fs_ctx_clear(&ctx);
- return ret;
-}
-
-static double convert_ssim_db(double _ssim, double _weight) {- return 10 * (log10(_weight) - log10(_weight - _ssim));
-}
-
-double vp9_calc_fastssim(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
- double *ssim_y, double *ssim_u, double *ssim_v) {- double ssimv;
- vp9_clear_system_state();
-
- *ssim_y = calc_ssim(source->y_buffer, source->y_stride, dest->y_buffer,
- dest->y_stride, source->y_crop_width,
- source->y_crop_height);
-
- *ssim_u = calc_ssim(source->u_buffer, source->uv_stride, dest->u_buffer,
- dest->uv_stride, source->uv_crop_width,
- source->uv_crop_height);
-
- *ssim_v = calc_ssim(source->v_buffer, source->uv_stride, dest->v_buffer,
- dest->uv_stride, source->uv_crop_width,
- source->uv_crop_height);
- ssimv = (*ssim_y) * .8 + .1 * ((*ssim_u) + (*ssim_v));
-
- return convert_ssim_db(ssimv, 1.0);
-}
--- a/vp9/encoder/vp9_psnrhvs.c
+++ /dev/null
@@ -1,224 +1,0 @@
-/*
- * Copyright (c) 2010 The WebM project authors. All Rights Reserved.
- *
- * Use of this source code is governed by a BSD-style license
- * that can be found in the LICENSE file in the root of the source
- * tree. An additional intellectual property rights grant can be found
- * in the file PATENTS. All contributing project authors may
- * be found in the AUTHORS file in the root of the source tree.
- *
- * This code was originally written by: Gregory Maxwell, at the Daala
- * project.
- */
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-
-#include "./vpx_config.h"
-#include "./vp9_rtcd.h"
-#include "./vpx_dsp_rtcd.h"
-#include "vp9/encoder/vp9_ssim.h"
-
-#if !defined(M_PI)
-# define M_PI (3.141592653589793238462643)
-#endif
-#include <string.h>
-
-void od_bin_fdct8x8(tran_low_t *y, int ystride, const int16_t *x, int xstride) {- (void) xstride;
- vpx_fdct8x8(x, y, ystride);
-}
-
-/* Normalized inverse quantization matrix for 8x8 DCT at the point of
- * transparency. This is not the JPEG based matrix from the paper,
- this one gives a slightly higher MOS agreement.*/
-float csf_y[8][8] = {{1.6193873005, 2.2901594831, 2.08509755623, 1.48366094411,- 1.00227514334, 0.678296995242, 0.466224900598, 0.3265091542}, {2.2901594831,- 1.94321815382, 2.04793073064, 1.68731108984, 1.2305666963, 0.868920337363,
- 0.61280991668, 0.436405793551}, {2.08509755623, 2.04793073064,- 1.34329019223, 1.09205635862, 0.875748795257, 0.670882927016,
- 0.501731932449, 0.372504254596}, {1.48366094411, 1.68731108984,- 1.09205635862, 0.772819797575, 0.605636379554, 0.48309405692,
- 0.380429446972, 0.295774038565}, {1.00227514334, 1.2305666963,- 0.875748795257, 0.605636379554, 0.448996256676, 0.352889268808,
- 0.283006984131, 0.226951348204}, {0.678296995242, 0.868920337363,- 0.670882927016, 0.48309405692, 0.352889268808, 0.27032073436,
- 0.215017739696, 0.17408067321}, {0.466224900598, 0.61280991668,- 0.501731932449, 0.380429446972, 0.283006984131, 0.215017739696,
- 0.168869545842, 0.136153931001}, {0.3265091542, 0.436405793551,- 0.372504254596, 0.295774038565, 0.226951348204, 0.17408067321,
- 0.136153931001, 0.109083846276}};
-float csf_cb420[8][8] = {- {1.91113096927, 2.46074210438, 1.18284184739, 1.14982565193, 1.05017074788,- 0.898018824055, 0.74725392039, 0.615105596242}, {2.46074210438,- 1.58529308355, 1.21363250036, 1.38190029285, 1.33100189972,
- 1.17428548929, 0.996404342439, 0.830890433625}, {1.18284184739,- 1.21363250036, 0.978712413627, 1.02624506078, 1.03145147362,
- 0.960060382087, 0.849823426169, 0.731221236837}, {1.14982565193,- 1.38190029285, 1.02624506078, 0.861317501629, 0.801821139099,
- 0.751437590932, 0.685398513368, 0.608694761374}, {1.05017074788,- 1.33100189972, 1.03145147362, 0.801821139099, 0.676555426187,
- 0.605503172737, 0.55002013668, 0.495804539034}, {0.898018824055,- 1.17428548929, 0.960060382087, 0.751437590932, 0.605503172737,
- 0.514674450957, 0.454353482512, 0.407050308965}, {0.74725392039,- 0.996404342439, 0.849823426169, 0.685398513368, 0.55002013668,
- 0.454353482512, 0.389234902883, 0.342353999733}, {0.615105596242,- 0.830890433625, 0.731221236837, 0.608694761374, 0.495804539034,
- 0.407050308965, 0.342353999733, 0.295530605237}};
-float csf_cr420[8][8] = {- {2.03871978502, 2.62502345193, 1.26180942886, 1.11019789803, 1.01397751469,- 0.867069376285, 0.721500455585, 0.593906509971}, {2.62502345193,- 1.69112867013, 1.17180569821, 1.3342742857, 1.28513006198,
- 1.13381474809, 0.962064122248, 0.802254508198}, {1.26180942886,- 1.17180569821, 0.944981930573, 0.990876405848, 0.995903384143,
- 0.926972725286, 0.820534991409, 0.706020324706}, {1.11019789803,- 1.3342742857, 0.990876405848, 0.831632933426, 0.77418706195,
- 0.725539939514, 0.661776842059, 0.587716619023}, {1.01397751469,- 1.28513006198, 0.995903384143, 0.77418706195, 0.653238524286,
- 0.584635025748, 0.531064164893, 0.478717061273}, {0.867069376285,- 1.13381474809, 0.926972725286, 0.725539939514, 0.584635025748,
- 0.496936637883, 0.438694579826, 0.393021669543}, {0.721500455585,- 0.962064122248, 0.820534991409, 0.661776842059, 0.531064164893,
- 0.438694579826, 0.375820256136, 0.330555063063}, {0.593906509971,- 0.802254508198, 0.706020324706, 0.587716619023, 0.478717061273,
- 0.393021669543, 0.330555063063, 0.285345396658}};
-
-static double convert_score_db(double _score, double _weight) {- return 10 * (log10(255 * 255) - log10(_weight * _score));
-}
-
-static double calc_psnrhvs(const unsigned char *_src, int _systride,
- const unsigned char *_dst, int _dystride,
- double _par, int _w, int _h, int _step,
- float _csf[8][8]) {- float ret;
- int16_t dct_s[8 * 8], dct_d[8 * 8];
- tran_low_t dct_s_coef[8 * 8], dct_d_coef[8 * 8];
- float mask[8][8];
- int pixels;
- int x;
- int y;
- (void) _par;
- ret = pixels = 0;
- /*In the PSNR-HVS-M paper[1] the authors describe the construction of
- their masking table as "we have used the quantization table for the
- color component Y of JPEG [6] that has been also obtained on the
- basis of CSF. Note that the values in quantization table JPEG have
- been normalized and then squared." Their CSF matrix (from PSNR-HVS)
- was also constructed from the JPEG matrices. I can not find any obvious
- scheme of normalizing to produce their table, but if I multiply their
- CSF by 0.38857 and square the result I get their masking table.
- I have no idea where this constant comes from, but deviating from it
- too greatly hurts MOS agreement.
-
- [1] Nikolay Ponomarenko, Flavia Silvestri, Karen Egiazarian, Marco Carli,
- Jaakko Astola, Vladimir Lukin, "On between-coefficient contrast masking
- of DCT basis functions", CD-ROM Proceedings of the Third
- International Workshop on Video Processing and Quality Metrics for Consumer
- Electronics VPQM-07, Scottsdale, Arizona, USA, 25-26 January, 2007, 4 p.*/
- for (x = 0; x < 8; x++)
- for (y = 0; y < 8; y++)
- mask[x][y] = (_csf[x][y] * 0.3885746225901003)
- * (_csf[x][y] * 0.3885746225901003);
- for (y = 0; y < _h - 7; y += _step) {- for (x = 0; x < _w - 7; x += _step) {- int i;
- int j;
- float s_means[4];
- float d_means[4];
- float s_vars[4];
- float d_vars[4];
- float s_gmean = 0;
- float d_gmean = 0;
- float s_gvar = 0;
- float d_gvar = 0;
- float s_mask = 0;
- float d_mask = 0;
- for (i = 0; i < 4; i++)
- s_means[i] = d_means[i] = s_vars[i] = d_vars[i] = 0;
- for (i = 0; i < 8; i++) {- for (j = 0; j < 8; j++) {- int sub = ((i & 12) >> 2) + ((j & 12) >> 1);
- dct_s[i * 8 + j] = _src[(y + i) * _systride + (j + x)];
- dct_d[i * 8 + j] = _dst[(y + i) * _dystride + (j + x)];
- s_gmean += dct_s[i * 8 + j];
- d_gmean += dct_d[i * 8 + j];
- s_means[sub] += dct_s[i * 8 + j];
- d_means[sub] += dct_d[i * 8 + j];
- }
- }
- s_gmean /= 64.f;
- d_gmean /= 64.f;
- for (i = 0; i < 4; i++)
- s_means[i] /= 16.f;
- for (i = 0; i < 4; i++)
- d_means[i] /= 16.f;
- for (i = 0; i < 8; i++) {- for (j = 0; j < 8; j++) {- int sub = ((i & 12) >> 2) + ((j & 12) >> 1);
- s_gvar += (dct_s[i * 8 + j] - s_gmean) * (dct_s[i * 8 + j] - s_gmean);
- d_gvar += (dct_d[i * 8 + j] - d_gmean) * (dct_d[i * 8 + j] - d_gmean);
- s_vars[sub] += (dct_s[i * 8 + j] - s_means[sub])
- * (dct_s[i * 8 + j] - s_means[sub]);
- d_vars[sub] += (dct_d[i * 8 + j] - d_means[sub])
- * (dct_d[i * 8 + j] - d_means[sub]);
- }
- }
- s_gvar *= 1 / 63.f * 64;
- d_gvar *= 1 / 63.f * 64;
- for (i = 0; i < 4; i++)
- s_vars[i] *= 1 / 15.f * 16;
- for (i = 0; i < 4; i++)
- d_vars[i] *= 1 / 15.f * 16;
- if (s_gvar > 0)
- s_gvar = (s_vars[0] + s_vars[1] + s_vars[2] + s_vars[3]) / s_gvar;
- if (d_gvar > 0)
- d_gvar = (d_vars[0] + d_vars[1] + d_vars[2] + d_vars[3]) / d_gvar;
- od_bin_fdct8x8(dct_s_coef, 8, dct_s, 8);
- od_bin_fdct8x8(dct_d_coef, 8, dct_d, 8);
- for (i = 0; i < 8; i++)
- for (j = (i == 0); j < 8; j++)
- s_mask += dct_s_coef[i * 8 + j] * dct_s_coef[i * 8 + j] * mask[i][j];
- for (i = 0; i < 8; i++)
- for (j = (i == 0); j < 8; j++)
- d_mask += dct_d_coef[i * 8 + j] * dct_d_coef[i * 8 + j] * mask[i][j];
- s_mask = sqrt(s_mask * s_gvar) / 32.f;
- d_mask = sqrt(d_mask * d_gvar) / 32.f;
- if (d_mask > s_mask)
- s_mask = d_mask;
- for (i = 0; i < 8; i++) {- for (j = 0; j < 8; j++) {- float err;
- err = fabs(dct_s_coef[i * 8 + j] - dct_d_coef[i * 8 + j]);
- if (i != 0 || j != 0)
- err = err < s_mask / mask[i][j] ? 0 : err - s_mask / mask[i][j];
- ret += (err * _csf[i][j]) * (err * _csf[i][j]);
- pixels++;
- }
- }
- }
- }
- ret /= pixels;
- return ret;
-}
-double vp9_psnrhvs(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
- double *y_psnrhvs, double *u_psnrhvs, double *v_psnrhvs) {- double psnrhvs;
- double par = 1.0;
- int step = 7;
- vp9_clear_system_state();
- *y_psnrhvs = calc_psnrhvs(source->y_buffer, source->y_stride, dest->y_buffer,
- dest->y_stride, par, source->y_crop_width,
- source->y_crop_height, step, csf_y);
-
- *u_psnrhvs = calc_psnrhvs(source->u_buffer, source->uv_stride, dest->u_buffer,
- dest->uv_stride, par, source->uv_crop_width,
- source->uv_crop_height, step, csf_cb420);
-
- *v_psnrhvs = calc_psnrhvs(source->v_buffer, source->uv_stride, dest->v_buffer,
- dest->uv_stride, par, source->uv_crop_width,
- source->uv_crop_height, step, csf_cr420);
- psnrhvs = (*y_psnrhvs) * .8 + .1 * ((*u_psnrhvs) + (*v_psnrhvs));
-
- return convert_score_db(psnrhvs, 1.0);
-}
--- a/vp9/encoder/vp9_ssim.c
+++ /dev/null
@@ -1,500 +1,0 @@
-/*
- * Copyright (c) 2010 The WebM project authors. All Rights Reserved.
- *
- * Use of this source code is governed by a BSD-style license
- * that can be found in the LICENSE file in the root of the source
- * tree. An additional intellectual property rights grant can be found
- * in the file PATENTS. All contributing project authors may
- * be found in the AUTHORS file in the root of the source tree.
- */
-
-#include <math.h>
-#include "./vp9_rtcd.h"
-#include "vpx_ports/mem.h"
-#include "vp9/encoder/vp9_ssim.h"
-
-void vp9_ssim_parms_16x16_c(uint8_t *s, int sp, uint8_t *r,
- int rp, unsigned long *sum_s, unsigned long *sum_r,
- unsigned long *sum_sq_s, unsigned long *sum_sq_r,
- unsigned long *sum_sxr) {- int i, j;
- for (i = 0; i < 16; i++, s += sp, r += rp) {- for (j = 0; j < 16; j++) {- *sum_s += s[j];
- *sum_r += r[j];
- *sum_sq_s += s[j] * s[j];
- *sum_sq_r += r[j] * r[j];
- *sum_sxr += s[j] * r[j];
- }
- }
-}
-void vp9_ssim_parms_8x8_c(uint8_t *s, int sp, uint8_t *r, int rp,
- unsigned long *sum_s, unsigned long *sum_r,
- unsigned long *sum_sq_s, unsigned long *sum_sq_r,
- unsigned long *sum_sxr) {- int i, j;
- for (i = 0; i < 8; i++, s += sp, r += rp) {- for (j = 0; j < 8; j++) {- *sum_s += s[j];
- *sum_r += r[j];
- *sum_sq_s += s[j] * s[j];
- *sum_sq_r += r[j] * r[j];
- *sum_sxr += s[j] * r[j];
- }
- }
-}
-
-#if CONFIG_VP9_HIGHBITDEPTH
-void vp9_highbd_ssim_parms_8x8_c(uint16_t *s, int sp, uint16_t *r, int rp,
- uint32_t *sum_s, uint32_t *sum_r,
- uint32_t *sum_sq_s, uint32_t *sum_sq_r,
- uint32_t *sum_sxr) {- int i, j;
- for (i = 0; i < 8; i++, s += sp, r += rp) {- for (j = 0; j < 8; j++) {- *sum_s += s[j];
- *sum_r += r[j];
- *sum_sq_s += s[j] * s[j];
- *sum_sq_r += r[j] * r[j];
- *sum_sxr += s[j] * r[j];
- }
- }
-}
-#endif // CONFIG_VP9_HIGHBITDEPTH
-
-static const int64_t cc1 = 26634; // (64^2*(.01*255)^2
-static const int64_t cc2 = 239708; // (64^2*(.03*255)^2
-
-static double similarity(unsigned long sum_s, unsigned long sum_r,
- unsigned long sum_sq_s, unsigned long sum_sq_r,
- unsigned long sum_sxr, int count) {- int64_t ssim_n, ssim_d;
- int64_t c1, c2;
-
- // scale the constants by number of pixels
- c1 = (cc1 * count * count) >> 12;
- c2 = (cc2 * count * count) >> 12;
-
- ssim_n = (2 * sum_s * sum_r + c1) * ((int64_t) 2 * count * sum_sxr -
- (int64_t) 2 * sum_s * sum_r + c2);
-
- ssim_d = (sum_s * sum_s + sum_r * sum_r + c1) *
- ((int64_t)count * sum_sq_s - (int64_t)sum_s * sum_s +
- (int64_t)count * sum_sq_r - (int64_t) sum_r * sum_r + c2);
-
- return ssim_n * 1.0 / ssim_d;
-}
-
-static double ssim_8x8(uint8_t *s, int sp, uint8_t *r, int rp) {- unsigned long sum_s = 0, sum_r = 0, sum_sq_s = 0, sum_sq_r = 0, sum_sxr = 0;
- vp9_ssim_parms_8x8(s, sp, r, rp, &sum_s, &sum_r, &sum_sq_s, &sum_sq_r,
- &sum_sxr);
- return similarity(sum_s, sum_r, sum_sq_s, sum_sq_r, sum_sxr, 64);
-}
-
-#if CONFIG_VP9_HIGHBITDEPTH
-static double highbd_ssim_8x8(uint16_t *s, int sp, uint16_t *r, int rp,
- unsigned int bd) {- uint32_t sum_s = 0, sum_r = 0, sum_sq_s = 0, sum_sq_r = 0, sum_sxr = 0;
- const int oshift = bd - 8;
- vp9_highbd_ssim_parms_8x8(s, sp, r, rp, &sum_s, &sum_r, &sum_sq_s, &sum_sq_r,
- &sum_sxr);
- return similarity(sum_s >> oshift,
- sum_r >> oshift,
- sum_sq_s >> (2 * oshift),
- sum_sq_r >> (2 * oshift),
- sum_sxr >> (2 * oshift),
- 64);
-}
-#endif // CONFIG_VP9_HIGHBITDEPTH
-
-// We are using a 8x8 moving window with starting location of each 8x8 window
-// on the 4x4 pixel grid. Such arrangement allows the windows to overlap
-// block boundaries to penalize blocking artifacts.
-double vp9_ssim2(uint8_t *img1, uint8_t *img2, int stride_img1,
- int stride_img2, int width, int height) {- int i, j;
- int samples = 0;
- double ssim_total = 0;
-
- // sample point start with each 4x4 location
- for (i = 0; i <= height - 8;
- i += 4, img1 += stride_img1 * 4, img2 += stride_img2 * 4) {- for (j = 0; j <= width - 8; j += 4) {- double v = ssim_8x8(img1 + j, stride_img1, img2 + j, stride_img2);
- ssim_total += v;
- samples++;
- }
- }
- ssim_total /= samples;
- return ssim_total;
-}
-
-#if CONFIG_VP9_HIGHBITDEPTH
-double vp9_highbd_ssim2(uint8_t *img1, uint8_t *img2, int stride_img1,
- int stride_img2, int width, int height,
- unsigned int bd) {- int i, j;
- int samples = 0;
- double ssim_total = 0;
-
- // sample point start with each 4x4 location
- for (i = 0; i <= height - 8;
- i += 4, img1 += stride_img1 * 4, img2 += stride_img2 * 4) {- for (j = 0; j <= width - 8; j += 4) {- double v = highbd_ssim_8x8(CONVERT_TO_SHORTPTR(img1 + j), stride_img1,
- CONVERT_TO_SHORTPTR(img2 + j), stride_img2,
- bd);
- ssim_total += v;
- samples++;
- }
- }
- ssim_total /= samples;
- return ssim_total;
-}
-#endif // CONFIG_VP9_HIGHBITDEPTH
-
-double vp9_calc_ssim(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
- double *weight) {- double a, b, c;
- double ssimv;
-
- a = vp9_ssim2(source->y_buffer, dest->y_buffer,
- source->y_stride, dest->y_stride,
- source->y_crop_width, source->y_crop_height);
-
- b = vp9_ssim2(source->u_buffer, dest->u_buffer,
- source->uv_stride, dest->uv_stride,
- source->uv_crop_width, source->uv_crop_height);
-
- c = vp9_ssim2(source->v_buffer, dest->v_buffer,
- source->uv_stride, dest->uv_stride,
- source->uv_crop_width, source->uv_crop_height);
-
- ssimv = a * .8 + .1 * (b + c);
-
- *weight = 1;
-
- return ssimv;
-}
-
-double vp9_calc_ssimg(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
- double *ssim_y, double *ssim_u, double *ssim_v) {- double ssim_all = 0;
- double a, b, c;
-
- a = vp9_ssim2(source->y_buffer, dest->y_buffer,
- source->y_stride, dest->y_stride,
- source->y_crop_width, source->y_crop_height);
-
- b = vp9_ssim2(source->u_buffer, dest->u_buffer,
- source->uv_stride, dest->uv_stride,
- source->uv_crop_width, source->uv_crop_height);
-
- c = vp9_ssim2(source->v_buffer, dest->v_buffer,
- source->uv_stride, dest->uv_stride,
- source->uv_crop_width, source->uv_crop_height);
- *ssim_y = a;
- *ssim_u = b;
- *ssim_v = c;
- ssim_all = (a * 4 + b + c) / 6;
-
- return ssim_all;
-}
-
-// traditional ssim as per: http://en.wikipedia.org/wiki/Structural_similarity
-//
-// Re working out the math ->
-//
-// ssim(x,y) = (2*mean(x)*mean(y) + c1)*(2*cov(x,y)+c2) /
-// ((mean(x)^2+mean(y)^2+c1)*(var(x)+var(y)+c2))
-//
-// mean(x) = sum(x) / n
-//
-// cov(x,y) = (n*sum(xi*yi)-sum(x)*sum(y))/(n*n)
-//
-// var(x) = (n*sum(xi*xi)-sum(xi)*sum(xi))/(n*n)
-//
-// ssim(x,y) =
-// (2*sum(x)*sum(y)/(n*n) + c1)*(2*(n*sum(xi*yi)-sum(x)*sum(y))/(n*n)+c2) /
-// (((sum(x)*sum(x)+sum(y)*sum(y))/(n*n) +c1) *
-// ((n*sum(xi*xi) - sum(xi)*sum(xi))/(n*n)+
-// (n*sum(yi*yi) - sum(yi)*sum(yi))/(n*n)+c2)))
-//
-// factoring out n*n
-//
-// ssim(x,y) =
-// (2*sum(x)*sum(y) + n*n*c1)*(2*(n*sum(xi*yi)-sum(x)*sum(y))+n*n*c2) /
-// (((sum(x)*sum(x)+sum(y)*sum(y)) + n*n*c1) *
-// (n*sum(xi*xi)-sum(xi)*sum(xi)+n*sum(yi*yi)-sum(yi)*sum(yi)+n*n*c2))
-//
-// Replace c1 with n*n * c1 for the final step that leads to this code:
-// The final step scales by 12 bits so we don't lose precision in the constants.
-
-double ssimv_similarity(Ssimv *sv, int64_t n) {- // Scale the constants by number of pixels.
- const int64_t c1 = (cc1 * n * n) >> 12;
- const int64_t c2 = (cc2 * n * n) >> 12;
-
- const double l = 1.0 * (2 * sv->sum_s * sv->sum_r + c1) /
- (sv->sum_s * sv->sum_s + sv->sum_r * sv->sum_r + c1);
-
- // Since these variables are unsigned sums, convert to double so
- // math is done in double arithmetic.
- const double v = (2.0 * n * sv->sum_sxr - 2 * sv->sum_s * sv->sum_r + c2)
- / (n * sv->sum_sq_s - sv->sum_s * sv->sum_s + n * sv->sum_sq_r
- - sv->sum_r * sv->sum_r + c2);
-
- return l * v;
-}
-
-// The first term of the ssim metric is a luminance factor.
-//
-// (2*mean(x)*mean(y) + c1)/ (mean(x)^2+mean(y)^2+c1)
-//
-// This luminance factor is super sensitive to the dark side of luminance
-// values and completely insensitive on the white side. check out 2 sets
-// (1,3) and (250,252) the term gives ( 2*1*3/(1+9) = .60
-// 2*250*252/ (250^2+252^2) => .99999997
-//
-// As a result in this tweaked version of the calculation in which the
-// luminance is taken as percentage off from peak possible.
-//
-// 255 * 255 - (sum_s - sum_r) / count * (sum_s - sum_r) / count
-//
-double ssimv_similarity2(Ssimv *sv, int64_t n) {- // Scale the constants by number of pixels.
- const int64_t c1 = (cc1 * n * n) >> 12;
- const int64_t c2 = (cc2 * n * n) >> 12;
-
- const double mean_diff = (1.0 * sv->sum_s - sv->sum_r) / n;
- const double l = (255 * 255 - mean_diff * mean_diff + c1) / (255 * 255 + c1);
-
- // Since these variables are unsigned, sums convert to double so
- // math is done in double arithmetic.
- const double v = (2.0 * n * sv->sum_sxr - 2 * sv->sum_s * sv->sum_r + c2)
- / (n * sv->sum_sq_s - sv->sum_s * sv->sum_s +
- n * sv->sum_sq_r - sv->sum_r * sv->sum_r + c2);
-
- return l * v;
-}
-void ssimv_parms(uint8_t *img1, int img1_pitch, uint8_t *img2, int img2_pitch,
- Ssimv *sv) {- vp9_ssim_parms_8x8(img1, img1_pitch, img2, img2_pitch,
- &sv->sum_s, &sv->sum_r, &sv->sum_sq_s, &sv->sum_sq_r,
- &sv->sum_sxr);
-}
-
-double vp9_get_ssim_metrics(uint8_t *img1, int img1_pitch,
- uint8_t *img2, int img2_pitch,
- int width, int height,
- Ssimv *sv2, Metrics *m,
- int do_inconsistency) {- double dssim_total = 0;
- double ssim_total = 0;
- double ssim2_total = 0;
- double inconsistency_total = 0;
- int i, j;
- int c = 0;
- double norm;
- double old_ssim_total = 0;
- vp9_clear_system_state();
- // We can sample points as frequently as we like start with 1 per 4x4.
- for (i = 0; i < height; i += 4,
- img1 += img1_pitch * 4, img2 += img2_pitch * 4) {- for (j = 0; j < width; j += 4, ++c) {- Ssimv sv = {0};- double ssim;
- double ssim2;
- double dssim;
- uint32_t var_new;
- uint32_t var_old;
- uint32_t mean_new;
- uint32_t mean_old;
- double ssim_new;
- double ssim_old;
-
- // Not sure there's a great way to handle the edge pixels
- // in ssim when using a window. Seems biased against edge pixels
- // however you handle this. This uses only samples that are
- // fully in the frame.
- if (j + 8 <= width && i + 8 <= height) {- ssimv_parms(img1 + j, img1_pitch, img2 + j, img2_pitch, &sv);
- }
-
- ssim = ssimv_similarity(&sv, 64);
- ssim2 = ssimv_similarity2(&sv, 64);
-
- sv.ssim = ssim2;
-
- // dssim is calculated to use as an actual error metric and
- // is scaled up to the same range as sum square error.
- // Since we are subsampling every 16th point maybe this should be
- // *16 ?
- dssim = 255 * 255 * (1 - ssim2) / 2;
-
- // Here I introduce a new error metric: consistency-weighted
- // SSIM-inconsistency. This metric isolates frames where the
- // SSIM 'suddenly' changes, e.g. if one frame in every 8 is much
- // sharper or blurrier than the others. Higher values indicate a
- // temporally inconsistent SSIM. There are two ideas at work:
- //
- // 1) 'SSIM-inconsistency': the total inconsistency value
- // reflects how much SSIM values are changing between this
- // source / reference frame pair and the previous pair.
- //
- // 2) 'consistency-weighted': weights de-emphasize areas in the
- // frame where the scene content has changed. Changes in scene
- // content are detected via changes in local variance and local
- // mean.
- //
- // Thus the overall measure reflects how inconsistent the SSIM
- // values are, over consistent regions of the frame.
- //
- // The metric has three terms:
- //
- // term 1 -> uses change in scene Variance to weight error score
- // 2 * var(Fi)*var(Fi-1) / (var(Fi)^2+var(Fi-1)^2)
- // larger changes from one frame to the next mean we care
- // less about consistency.
- //
- // term 2 -> uses change in local scene luminance to weight error
- // 2 * avg(Fi)*avg(Fi-1) / (avg(Fi)^2+avg(Fi-1)^2)
- // larger changes from one frame to the next mean we care
- // less about consistency.
- //
- // term3 -> measures inconsistency in ssim scores between frames
- // 1 - ( 2 * ssim(Fi)*ssim(Fi-1)/(ssim(Fi)^2+sssim(Fi-1)^2).
- //
- // This term compares the ssim score for the same location in 2
- // subsequent frames.
- var_new = sv.sum_sq_s - sv.sum_s * sv.sum_s / 64;
- var_old = sv2[c].sum_sq_s - sv2[c].sum_s * sv2[c].sum_s / 64;
- mean_new = sv.sum_s;
- mean_old = sv2[c].sum_s;
- ssim_new = sv.ssim;
- ssim_old = sv2[c].ssim;
-
- if (do_inconsistency) {- // We do the metric once for every 4x4 block in the image. Since
- // we are scaling the error to SSE for use in a psnr calculation
- // 1.0 = 4x4x255x255 the worst error we can possibly have.
- static const double kScaling = 4. * 4 * 255 * 255;
-
- // The constants have to be non 0 to avoid potential divide by 0
- // issues other than that they affect kind of a weighting between
- // the terms. No testing of what the right terms should be has been
- // done.
- static const double c1 = 1, c2 = 1, c3 = 1;
-
- // This measures how much consistent variance is in two consecutive
- // source frames. 1.0 means they have exactly the same variance.
- const double variance_term = (2.0 * var_old * var_new + c1) /
- (1.0 * var_old * var_old + 1.0 * var_new * var_new + c1);
-
- // This measures how consistent the local mean are between two
- // consecutive frames. 1.0 means they have exactly the same mean.
- const double mean_term = (2.0 * mean_old * mean_new + c2) /
- (1.0 * mean_old * mean_old + 1.0 * mean_new * mean_new + c2);
-
- // This measures how consistent the ssims of two
- // consecutive frames is. 1.0 means they are exactly the same.
- double ssim_term = pow((2.0 * ssim_old * ssim_new + c3) /
- (ssim_old * ssim_old + ssim_new * ssim_new + c3),
- 5);
-
- double this_inconsistency;
-
- // Floating point math sometimes makes this > 1 by a tiny bit.
- // We want the metric to scale between 0 and 1.0 so we can convert
- // it to an snr scaled value.
- if (ssim_term > 1)
- ssim_term = 1;
-
- // This converts the consistency metric to an inconsistency metric
- // ( so we can scale it like psnr to something like sum square error.
- // The reason for the variance and mean terms is the assumption that
- // if there are big changes in the source we shouldn't penalize
- // inconsistency in ssim scores a bit less as it will be less visible
- // to the user.
- this_inconsistency = (1 - ssim_term) * variance_term * mean_term;
-
- this_inconsistency *= kScaling;
- inconsistency_total += this_inconsistency;
- }
- sv2[c] = sv;
- ssim_total += ssim;
- ssim2_total += ssim2;
- dssim_total += dssim;
-
- old_ssim_total += ssim_old;
- }
- old_ssim_total += 0;
- }
-
- norm = 1. / (width / 4) / (height / 4);
- ssim_total *= norm;
- ssim2_total *= norm;
- m->ssim2 = ssim2_total;
- m->ssim = ssim_total;
- if (old_ssim_total == 0)
- inconsistency_total = 0;
-
- m->ssimc = inconsistency_total;
-
- m->dssim = dssim_total;
- return inconsistency_total;
-}
-
-
-#if CONFIG_VP9_HIGHBITDEPTH
-double vp9_highbd_calc_ssim(YV12_BUFFER_CONFIG *source,
- YV12_BUFFER_CONFIG *dest,
- double *weight, unsigned int bd) {- double a, b, c;
- double ssimv;
-
- a = vp9_highbd_ssim2(source->y_buffer, dest->y_buffer,
- source->y_stride, dest->y_stride,
- source->y_crop_width, source->y_crop_height, bd);
-
- b = vp9_highbd_ssim2(source->u_buffer, dest->u_buffer,
- source->uv_stride, dest->uv_stride,
- source->uv_crop_width, source->uv_crop_height, bd);
-
- c = vp9_highbd_ssim2(source->v_buffer, dest->v_buffer,
- source->uv_stride, dest->uv_stride,
- source->uv_crop_width, source->uv_crop_height, bd);
-
- ssimv = a * .8 + .1 * (b + c);
-
- *weight = 1;
-
- return ssimv;
-}
-
-double vp9_highbd_calc_ssimg(YV12_BUFFER_CONFIG *source,
- YV12_BUFFER_CONFIG *dest, double *ssim_y,
- double *ssim_u, double *ssim_v, unsigned int bd) {- double ssim_all = 0;
- double a, b, c;
-
- a = vp9_highbd_ssim2(source->y_buffer, dest->y_buffer,
- source->y_stride, dest->y_stride,
- source->y_crop_width, source->y_crop_height, bd);
-
- b = vp9_highbd_ssim2(source->u_buffer, dest->u_buffer,
- source->uv_stride, dest->uv_stride,
- source->uv_crop_width, source->uv_crop_height, bd);
-
- c = vp9_highbd_ssim2(source->v_buffer, dest->v_buffer,
- source->uv_stride, dest->uv_stride,
- source->uv_crop_width, source->uv_crop_height, bd);
- *ssim_y = a;
- *ssim_u = b;
- *ssim_v = c;
- ssim_all = (a * 4 + b + c) / 6;
-
- return ssim_all;
-}
-#endif // CONFIG_VP9_HIGHBITDEPTH
--- a/vp9/encoder/vp9_ssim.h
+++ /dev/null
@@ -1,96 +1,0 @@
-/*
- * Copyright (c) 2014 The WebM project authors. All Rights Reserved.
- *
- * Use of this source code is governed by a BSD-style license
- * that can be found in the LICENSE file in the root of the source
- * tree. An additional intellectual property rights grant can be found
- * in the file PATENTS. All contributing project authors may
- * be found in the AUTHORS file in the root of the source tree.
- */
-
-#ifndef VP9_ENCODER_VP9_SSIM_H_
-#define VP9_ENCODER_VP9_SSIM_H_
-
-#ifdef __cplusplus
-extern "C" {-#endif
-
-#include "vpx_scale/yv12config.h"
-
-// metrics used for calculating ssim, ssim2, dssim, and ssimc
-typedef struct {- // source sum ( over 8x8 region )
- uint64_t sum_s;
-
- // reference sum (over 8x8 region )
- uint64_t sum_r;
-
- // source sum squared ( over 8x8 region )
- uint64_t sum_sq_s;
-
- // reference sum squared (over 8x8 region )
- uint64_t sum_sq_r;
-
- // sum of source times reference (over 8x8 region)
- uint64_t sum_sxr;
-
- // calculated ssim score between source and reference
- double ssim;
-} Ssimv;
-
-// metrics collected on a frame basis
-typedef struct {- // ssim consistency error metric ( see code for explanation )
- double ssimc;
-
- // standard ssim
- double ssim;
-
- // revised ssim ( see code for explanation)
- double ssim2;
-
- // ssim restated as an error metric like sse
- double dssim;
-
- // dssim converted to decibels
- double dssimd;
-
- // ssimc converted to decibels
- double ssimcd;
-} Metrics;
-
-double vp9_get_ssim_metrics(uint8_t *img1, int img1_pitch, uint8_t *img2,
- int img2_pitch, int width, int height, Ssimv *sv2,
- Metrics *m, int do_inconsistency);
-
-double vp9_calc_ssim(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
- double *weight);
-
-double vp9_calc_ssimg(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
- double *ssim_y, double *ssim_u, double *ssim_v);
-
-double vp9_calc_fastssim(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
- double *ssim_y, double *ssim_u, double *ssim_v);
-
-double vp9_psnrhvs(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
- double *ssim_y, double *ssim_u, double *ssim_v);
-
-#if CONFIG_VP9_HIGHBITDEPTH
-double vp9_highbd_calc_ssim(YV12_BUFFER_CONFIG *source,
- YV12_BUFFER_CONFIG *dest,
- double *weight,
- unsigned int bd);
-
-double vp9_highbd_calc_ssimg(YV12_BUFFER_CONFIG *source,
- YV12_BUFFER_CONFIG *dest,
- double *ssim_y,
- double *ssim_u,
- double *ssim_v,
- unsigned int bd);
-#endif // CONFIG_VP9_HIGHBITDEPTH
-
-#ifdef __cplusplus
-} // extern "C"
-#endif
-
-#endif // VP9_ENCODER_VP9_SSIM_H_
--- a/vp9/encoder/x86/vp9_ssim_opt_x86_64.asm
+++ /dev/null
@@ -1,216 +1,0 @@
-;
-; Copyright (c) 2010 The WebM project authors. All Rights Reserved.
-;
-; Use of this source code is governed by a BSD-style license
-; that can be found in the LICENSE file in the root of the source
-; tree. An additional intellectual property rights grant can be found
-; in the file PATENTS. All contributing project authors may
-; be found in the AUTHORS file in the root of the source tree.
-;
-
-%include "vpx_ports/x86_abi_support.asm"
-
-; tabulate_ssim - sums sum_s,sum_r,sum_sq_s,sum_sq_r, sum_sxr
-%macro TABULATE_SSIM 0
- paddusw xmm15, xmm3 ; sum_s
- paddusw xmm14, xmm4 ; sum_r
- movdqa xmm1, xmm3
- pmaddwd xmm1, xmm1
- paddd xmm13, xmm1 ; sum_sq_s
- movdqa xmm2, xmm4
- pmaddwd xmm2, xmm2
- paddd xmm12, xmm2 ; sum_sq_r
- pmaddwd xmm3, xmm4
- paddd xmm11, xmm3 ; sum_sxr
-%endmacro
-
-; Sum across the register %1 starting with q words
-%macro SUM_ACROSS_Q 1
- movdqa xmm2,%1
- punpckldq %1,xmm0
- punpckhdq xmm2,xmm0
- paddq %1,xmm2
- movdqa xmm2,%1
- punpcklqdq %1,xmm0
- punpckhqdq xmm2,xmm0
- paddq %1,xmm2
-%endmacro
-
-; Sum across the register %1 starting with q words
-%macro SUM_ACROSS_W 1
- movdqa xmm1, %1
- punpcklwd %1,xmm0
- punpckhwd xmm1,xmm0
- paddd %1, xmm1
- SUM_ACROSS_Q %1
-%endmacro
-;void ssim_parms_sse2(
-; unsigned char *s,
-; int sp,
-; unsigned char *r,
-; int rp
-; unsigned long *sum_s,
-; unsigned long *sum_r,
-; unsigned long *sum_sq_s,
-; unsigned long *sum_sq_r,
-; unsigned long *sum_sxr);
-;
-; TODO: Use parm passing through structure, probably don't need the pxors
-; ( calling app will initialize to 0 ) could easily fit everything in sse2
-; without too much hastle, and can probably do better estimates with psadw
-; or pavgb At this point this is just meant to be first pass for calculating
-; all the parms needed for 16x16 ssim so we can play with dssim as distortion
-; in mode selection code.
-global sym(vp9_ssim_parms_16x16_sse2) PRIVATE
-sym(vp9_ssim_parms_16x16_sse2):
- push rbp
- mov rbp, rsp
- SHADOW_ARGS_TO_STACK 9
- SAVE_XMM 15
- push rsi
- push rdi
- ; end prolog
-
- mov rsi, arg(0) ;s
- mov rcx, arg(1) ;sp
- mov rdi, arg(2) ;r
- mov rax, arg(3) ;rp
-
- pxor xmm0, xmm0
- pxor xmm15,xmm15 ;sum_s
- pxor xmm14,xmm14 ;sum_r
- pxor xmm13,xmm13 ;sum_sq_s
- pxor xmm12,xmm12 ;sum_sq_r
- pxor xmm11,xmm11 ;sum_sxr
-
- mov rdx, 16 ;row counter
-.NextRow:
-
- ;grab source and reference pixels
- movdqu xmm5, [rsi]
- movdqu xmm6, [rdi]
- movdqa xmm3, xmm5
- movdqa xmm4, xmm6
- punpckhbw xmm3, xmm0 ; high_s
- punpckhbw xmm4, xmm0 ; high_r
-
- TABULATE_SSIM
-
- movdqa xmm3, xmm5
- movdqa xmm4, xmm6
- punpcklbw xmm3, xmm0 ; low_s
- punpcklbw xmm4, xmm0 ; low_r
-
- TABULATE_SSIM
-
- add rsi, rcx ; next s row
- add rdi, rax ; next r row
-
- dec rdx ; counter
- jnz .NextRow
-
- SUM_ACROSS_W xmm15
- SUM_ACROSS_W xmm14
- SUM_ACROSS_Q xmm13
- SUM_ACROSS_Q xmm12
- SUM_ACROSS_Q xmm11
-
- mov rdi,arg(4)
- movd [rdi], xmm15;
- mov rdi,arg(5)
- movd [rdi], xmm14;
- mov rdi,arg(6)
- movd [rdi], xmm13;
- mov rdi,arg(7)
- movd [rdi], xmm12;
- mov rdi,arg(8)
- movd [rdi], xmm11;
-
- ; begin epilog
- pop rdi
- pop rsi
- RESTORE_XMM
- UNSHADOW_ARGS
- pop rbp
- ret
-
-;void ssim_parms_sse2(
-; unsigned char *s,
-; int sp,
-; unsigned char *r,
-; int rp
-; unsigned long *sum_s,
-; unsigned long *sum_r,
-; unsigned long *sum_sq_s,
-; unsigned long *sum_sq_r,
-; unsigned long *sum_sxr);
-;
-; TODO: Use parm passing through structure, probably don't need the pxors
-; ( calling app will initialize to 0 ) could easily fit everything in sse2
-; without too much hastle, and can probably do better estimates with psadw
-; or pavgb At this point this is just meant to be first pass for calculating
-; all the parms needed for 16x16 ssim so we can play with dssim as distortion
-; in mode selection code.
-global sym(vp9_ssim_parms_8x8_sse2) PRIVATE
-sym(vp9_ssim_parms_8x8_sse2):
- push rbp
- mov rbp, rsp
- SHADOW_ARGS_TO_STACK 9
- SAVE_XMM 15
- push rsi
- push rdi
- ; end prolog
-
- mov rsi, arg(0) ;s
- mov rcx, arg(1) ;sp
- mov rdi, arg(2) ;r
- mov rax, arg(3) ;rp
-
- pxor xmm0, xmm0
- pxor xmm15,xmm15 ;sum_s
- pxor xmm14,xmm14 ;sum_r
- pxor xmm13,xmm13 ;sum_sq_s
- pxor xmm12,xmm12 ;sum_sq_r
- pxor xmm11,xmm11 ;sum_sxr
-
- mov rdx, 8 ;row counter
-.NextRow:
-
- ;grab source and reference pixels
- movq xmm3, [rsi]
- movq xmm4, [rdi]
- punpcklbw xmm3, xmm0 ; low_s
- punpcklbw xmm4, xmm0 ; low_r
-
- TABULATE_SSIM
-
- add rsi, rcx ; next s row
- add rdi, rax ; next r row
-
- dec rdx ; counter
- jnz .NextRow
-
- SUM_ACROSS_W xmm15
- SUM_ACROSS_W xmm14
- SUM_ACROSS_Q xmm13
- SUM_ACROSS_Q xmm12
- SUM_ACROSS_Q xmm11
-
- mov rdi,arg(4)
- movd [rdi], xmm15;
- mov rdi,arg(5)
- movd [rdi], xmm14;
- mov rdi,arg(6)
- movd [rdi], xmm13;
- mov rdi,arg(7)
- movd [rdi], xmm12;
- mov rdi,arg(8)
- movd [rdi], xmm11;
-
- ; begin epilog
- pop rdi
- pop rsi
- RESTORE_XMM
- UNSHADOW_ARGS
- pop rbp
- ret
--- a/vp9/vp9cx.mk
+++ b/vp9/vp9cx.mk
@@ -33,7 +33,6 @@
VP9_CX_SRCS-yes += encoder/vp9_ethread.h
VP9_CX_SRCS-yes += encoder/vp9_ethread.c
VP9_CX_SRCS-yes += encoder/vp9_extend.c
-VP9_CX_SRCS-$(CONFIG_INTERNAL_STATS) += encoder/vp9_fastssim.c
VP9_CX_SRCS-yes += encoder/vp9_firstpass.c
VP9_CX_SRCS-yes += encoder/vp9_block.h
VP9_CX_SRCS-yes += encoder/vp9_bitstream.h
@@ -57,7 +56,6 @@
VP9_CX_SRCS-yes += encoder/vp9_encoder.c
VP9_CX_SRCS-yes += encoder/vp9_picklpf.c
VP9_CX_SRCS-yes += encoder/vp9_picklpf.h
-VP9_CX_SRCS-$(CONFIG_INTERNAL_STATS) += encoder/vp9_psnrhvs.c
VP9_CX_SRCS-yes += encoder/vp9_quantize.c
VP9_CX_SRCS-yes += encoder/vp9_ratectrl.c
VP9_CX_SRCS-yes += encoder/vp9_rd.c
@@ -72,8 +70,6 @@
VP9_CX_SRCS-yes += encoder/vp9_svc_layercontext.c
VP9_CX_SRCS-yes += encoder/vp9_resize.c
VP9_CX_SRCS-yes += encoder/vp9_resize.h
-VP9_CX_SRCS-$(CONFIG_INTERNAL_STATS) += encoder/vp9_ssim.c
-VP9_CX_SRCS-$(CONFIG_INTERNAL_STATS) += encoder/vp9_ssim.h
VP9_CX_SRCS-$(CONFIG_INTERNAL_STATS) += encoder/vp9_blockiness.c
VP9_CX_SRCS-yes += encoder/vp9_tokenize.c
@@ -113,7 +109,6 @@
VP9_CX_SRCS-$(HAVE_SSSE3) += encoder/x86/vp9_dct_ssse3_x86_64.asm
endif
endif
-VP9_CX_SRCS-$(ARCH_X86_64) += encoder/x86/vp9_ssim_opt_x86_64.asm
VP9_CX_SRCS-$(HAVE_SSE2) += encoder/x86/vp9_dct_sse2.c
VP9_CX_SRCS-$(HAVE_SSSE3) += encoder/x86/vp9_dct_ssse3.c
--- /dev/null
+++ b/vpx_dsp/fastssim.c
@@ -1,0 +1,465 @@
+/*
+ * Copyright (c) 2010 The WebM project authors. All Rights Reserved.
+ *
+ * Use of this source code is governed by a BSD-style license
+ * that can be found in the LICENSE file in the root of the source
+ * tree. An additional intellectual property rights grant can be found
+ * in the file PATENTS. All contributing project authors may
+ * be found in the AUTHORS file in the root of the source tree.
+ *
+ * This code was originally written by: Nathan E. Egge, at the Daala
+ * project.
+ */
+#include <math.h>
+#include <string.h>
+#include "./vpx_config.h"
+#include "./vpx_dsp_rtcd.h"
+#include "vpx_dsp/ssim.h"
+/* TODO(jbb): High bit depth version of this code needed */
+typedef struct fs_level fs_level;
+typedef struct fs_ctx fs_ctx;
+
+#define SSIM_C1 (255 * 255 * 0.01 * 0.01)
+#define SSIM_C2 (255 * 255 * 0.03 * 0.03)
+
+#define FS_MINI(_a, _b) ((_a) < (_b) ? (_a) : (_b))
+#define FS_MAXI(_a, _b) ((_a) > (_b) ? (_a) : (_b))
+
+struct fs_level {+ uint16_t *im1;
+ uint16_t *im2;
+ double *ssim;
+ int w;
+ int h;
+};
+
+struct fs_ctx {+ fs_level *level;
+ int nlevels;
+ unsigned *col_buf;
+};
+
+static void fs_ctx_init(fs_ctx *_ctx, int _w, int _h, int _nlevels) {+ unsigned char *data;
+ size_t data_size;
+ int lw;
+ int lh;
+ int l;
+ lw = (_w + 1) >> 1;
+ lh = (_h + 1) >> 1;
+ data_size = _nlevels * sizeof(fs_level)
+ + 2 * (lw + 8) * 8 * sizeof(*_ctx->col_buf);
+ for (l = 0; l < _nlevels; l++) {+ size_t im_size;
+ size_t level_size;
+ im_size = lw * (size_t) lh;
+ level_size = 2 * im_size * sizeof(*_ctx->level[l].im1);
+ level_size += sizeof(*_ctx->level[l].ssim) - 1;
+ level_size /= sizeof(*_ctx->level[l].ssim);
+ level_size += im_size;
+ level_size *= sizeof(*_ctx->level[l].ssim);
+ data_size += level_size;
+ lw = (lw + 1) >> 1;
+ lh = (lh + 1) >> 1;
+ }
+ data = (unsigned char *) malloc(data_size);
+ _ctx->level = (fs_level *) data;
+ _ctx->nlevels = _nlevels;
+ data += _nlevels * sizeof(*_ctx->level);
+ lw = (_w + 1) >> 1;
+ lh = (_h + 1) >> 1;
+ for (l = 0; l < _nlevels; l++) {+ size_t im_size;
+ size_t level_size;
+ _ctx->level[l].w = lw;
+ _ctx->level[l].h = lh;
+ im_size = lw * (size_t) lh;
+ level_size = 2 * im_size * sizeof(*_ctx->level[l].im1);
+ level_size += sizeof(*_ctx->level[l].ssim) - 1;
+ level_size /= sizeof(*_ctx->level[l].ssim);
+ level_size *= sizeof(*_ctx->level[l].ssim);
+ _ctx->level[l].im1 = (uint16_t *) data;
+ _ctx->level[l].im2 = _ctx->level[l].im1 + im_size;
+ data += level_size;
+ _ctx->level[l].ssim = (double *) data;
+ data += im_size * sizeof(*_ctx->level[l].ssim);
+ lw = (lw + 1) >> 1;
+ lh = (lh + 1) >> 1;
+ }
+ _ctx->col_buf = (unsigned *) data;
+}
+
+static void fs_ctx_clear(fs_ctx *_ctx) {+ free(_ctx->level);
+}
+
+static void fs_downsample_level(fs_ctx *_ctx, int _l) {+ const uint16_t *src1;
+ const uint16_t *src2;
+ uint16_t *dst1;
+ uint16_t *dst2;
+ int w2;
+ int h2;
+ int w;
+ int h;
+ int i;
+ int j;
+ w = _ctx->level[_l].w;
+ h = _ctx->level[_l].h;
+ dst1 = _ctx->level[_l].im1;
+ dst2 = _ctx->level[_l].im2;
+ w2 = _ctx->level[_l - 1].w;
+ h2 = _ctx->level[_l - 1].h;
+ src1 = _ctx->level[_l - 1].im1;
+ src2 = _ctx->level[_l - 1].im2;
+ for (j = 0; j < h; j++) {+ int j0offs;
+ int j1offs;
+ j0offs = 2 * j * w2;
+ j1offs = FS_MINI(2 * j + 1, h2) * w2;
+ for (i = 0; i < w; i++) {+ int i0;
+ int i1;
+ i0 = 2 * i;
+ i1 = FS_MINI(i0 + 1, w2);
+ dst1[j * w + i] = src1[j0offs + i0] + src1[j0offs + i1]
+ + src1[j1offs + i0] + src1[j1offs + i1];
+ dst2[j * w + i] = src2[j0offs + i0] + src2[j0offs + i1]
+ + src2[j1offs + i0] + src2[j1offs + i1];
+ }
+ }
+}
+
+static void fs_downsample_level0(fs_ctx *_ctx, const unsigned char *_src1,
+ int _s1ystride, const unsigned char *_src2,
+ int _s2ystride, int _w, int _h) {+ uint16_t *dst1;
+ uint16_t *dst2;
+ int w;
+ int h;
+ int i;
+ int j;
+ w = _ctx->level[0].w;
+ h = _ctx->level[0].h;
+ dst1 = _ctx->level[0].im1;
+ dst2 = _ctx->level[0].im2;
+ for (j = 0; j < h; j++) {+ int j0;
+ int j1;
+ j0 = 2 * j;
+ j1 = FS_MINI(j0 + 1, _h);
+ for (i = 0; i < w; i++) {+ int i0;
+ int i1;
+ i0 = 2 * i;
+ i1 = FS_MINI(i0 + 1, _w);
+ dst1[j * w + i] = _src1[j0 * _s1ystride + i0]
+ + _src1[j0 * _s1ystride + i1] + _src1[j1 * _s1ystride + i0]
+ + _src1[j1 * _s1ystride + i1];
+ dst2[j * w + i] = _src2[j0 * _s2ystride + i0]
+ + _src2[j0 * _s2ystride + i1] + _src2[j1 * _s2ystride + i0]
+ + _src2[j1 * _s2ystride + i1];
+ }
+ }
+}
+
+static void fs_apply_luminance(fs_ctx *_ctx, int _l) {+ unsigned *col_sums_x;
+ unsigned *col_sums_y;
+ uint16_t *im1;
+ uint16_t *im2;
+ double *ssim;
+ double c1;
+ int w;
+ int h;
+ int j0offs;
+ int j1offs;
+ int i;
+ int j;
+ w = _ctx->level[_l].w;
+ h = _ctx->level[_l].h;
+ col_sums_x = _ctx->col_buf;
+ col_sums_y = col_sums_x + w;
+ im1 = _ctx->level[_l].im1;
+ im2 = _ctx->level[_l].im2;
+ for (i = 0; i < w; i++)
+ col_sums_x[i] = 5 * im1[i];
+ for (i = 0; i < w; i++)
+ col_sums_y[i] = 5 * im2[i];
+ for (j = 1; j < 4; j++) {+ j1offs = FS_MINI(j, h - 1) * w;
+ for (i = 0; i < w; i++)
+ col_sums_x[i] += im1[j1offs + i];
+ for (i = 0; i < w; i++)
+ col_sums_y[i] += im2[j1offs + i];
+ }
+ ssim = _ctx->level[_l].ssim;
+ c1 = (double) (SSIM_C1 * 4096 * (1 << 4 * _l));
+ for (j = 0; j < h; j++) {+ unsigned mux;
+ unsigned muy;
+ int i0;
+ int i1;
+ mux = 5 * col_sums_x[0];
+ muy = 5 * col_sums_y[0];
+ for (i = 1; i < 4; i++) {+ i1 = FS_MINI(i, w - 1);
+ mux += col_sums_x[i1];
+ muy += col_sums_y[i1];
+ }
+ for (i = 0; i < w; i++) {+ ssim[j * w + i] *= (2 * mux * (double) muy + c1)
+ / (mux * (double) mux + muy * (double) muy + c1);
+ if (i + 1 < w) {+ i0 = FS_MAXI(0, i - 4);
+ i1 = FS_MINI(i + 4, w - 1);
+ mux += col_sums_x[i1] - col_sums_x[i0];
+ muy += col_sums_x[i1] - col_sums_x[i0];
+ }
+ }
+ if (j + 1 < h) {+ j0offs = FS_MAXI(0, j - 4) * w;
+ for (i = 0; i < w; i++)
+ col_sums_x[i] -= im1[j0offs + i];
+ for (i = 0; i < w; i++)
+ col_sums_y[i] -= im2[j0offs + i];
+ j1offs = FS_MINI(j + 4, h - 1) * w;
+ for (i = 0; i < w; i++)
+ col_sums_x[i] += im1[j1offs + i];
+ for (i = 0; i < w; i++)
+ col_sums_y[i] += im2[j1offs + i];
+ }
+ }
+}
+
+#define FS_COL_SET(_col, _joffs, _ioffs) \
+ do { \+ unsigned gx; \
+ unsigned gy; \
+ gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
+ gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
+ col_sums_gx2[(_col)] = gx * (double)gx; \
+ col_sums_gy2[(_col)] = gy * (double)gy; \
+ col_sums_gxgy[(_col)] = gx * (double)gy; \
+ } \
+ while (0)
+
+#define FS_COL_ADD(_col, _joffs, _ioffs) \
+ do { \+ unsigned gx; \
+ unsigned gy; \
+ gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
+ gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
+ col_sums_gx2[(_col)] += gx * (double)gx; \
+ col_sums_gy2[(_col)] += gy * (double)gy; \
+ col_sums_gxgy[(_col)] += gx * (double)gy; \
+ } \
+ while (0)
+
+#define FS_COL_SUB(_col, _joffs, _ioffs) \
+ do { \+ unsigned gx; \
+ unsigned gy; \
+ gx = gx_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
+ gy = gy_buf[((j + (_joffs)) & 7) * stride + i + (_ioffs)]; \
+ col_sums_gx2[(_col)] -= gx * (double)gx; \
+ col_sums_gy2[(_col)] -= gy * (double)gy; \
+ col_sums_gxgy[(_col)] -= gx * (double)gy; \
+ } \
+ while (0)
+
+#define FS_COL_COPY(_col1, _col2) \
+ do { \+ col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)]; \
+ col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)]; \
+ col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)]; \
+ } \
+ while (0)
+
+#define FS_COL_HALVE(_col1, _col2) \
+ do { \+ col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)] * 0.5; \
+ col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)] * 0.5; \
+ col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)] * 0.5; \
+ } \
+ while (0)
+
+#define FS_COL_DOUBLE(_col1, _col2) \
+ do { \+ col_sums_gx2[(_col1)] = col_sums_gx2[(_col2)] * 2; \
+ col_sums_gy2[(_col1)] = col_sums_gy2[(_col2)] * 2; \
+ col_sums_gxgy[(_col1)] = col_sums_gxgy[(_col2)] * 2; \
+ } \
+ while (0)
+
+static void fs_calc_structure(fs_ctx *_ctx, int _l) {+ uint16_t *im1;
+ uint16_t *im2;
+ unsigned *gx_buf;
+ unsigned *gy_buf;
+ double *ssim;
+ double col_sums_gx2[8];
+ double col_sums_gy2[8];
+ double col_sums_gxgy[8];
+ double c2;
+ int stride;
+ int w;
+ int h;
+ int i;
+ int j;
+ w = _ctx->level[_l].w;
+ h = _ctx->level[_l].h;
+ im1 = _ctx->level[_l].im1;
+ im2 = _ctx->level[_l].im2;
+ ssim = _ctx->level[_l].ssim;
+ gx_buf = _ctx->col_buf;
+ stride = w + 8;
+ gy_buf = gx_buf + 8 * stride;
+ memset(gx_buf, 0, 2 * 8 * stride * sizeof(*gx_buf));
+ c2 = SSIM_C2 * (1 << 4 * _l) * 16 * 104;
+ for (j = 0; j < h + 4; j++) {+ if (j < h - 1) {+ for (i = 0; i < w - 1; i++) {+ unsigned g1;
+ unsigned g2;
+ unsigned gx;
+ unsigned gy;
+ g1 = abs(im1[(j + 1) * w + i + 1] - im1[j * w + i]);
+ g2 = abs(im1[(j + 1) * w + i] - im1[j * w + i + 1]);
+ gx = 4 * FS_MAXI(g1, g2) + FS_MINI(g1, g2);
+ g1 = abs(im2[(j + 1) * w + i + 1] - im2[j * w + i]);
+ g2 = abs(im2[(j + 1) * w + i] - im2[j * w + i + 1]);
+ gy = 4 * FS_MAXI(g1, g2) + FS_MINI(g1, g2);
+ gx_buf[(j & 7) * stride + i + 4] = gx;
+ gy_buf[(j & 7) * stride + i + 4] = gy;
+ }
+ } else {+ memset(gx_buf + (j & 7) * stride, 0, stride * sizeof(*gx_buf));
+ memset(gy_buf + (j & 7) * stride, 0, stride * sizeof(*gy_buf));
+ }
+ if (j >= 4) {+ int k;
+ col_sums_gx2[3] = col_sums_gx2[2] = col_sums_gx2[1] = col_sums_gx2[0] = 0;
+ col_sums_gy2[3] = col_sums_gy2[2] = col_sums_gy2[1] = col_sums_gy2[0] = 0;
+ col_sums_gxgy[3] = col_sums_gxgy[2] = col_sums_gxgy[1] =
+ col_sums_gxgy[0] = 0;
+ for (i = 4; i < 8; i++) {+ FS_COL_SET(i, -1, 0);
+ FS_COL_ADD(i, 0, 0);
+ for (k = 1; k < 8 - i; k++) {+ FS_COL_DOUBLE(i, i);
+ FS_COL_ADD(i, -k - 1, 0);
+ FS_COL_ADD(i, k, 0);
+ }
+ }
+ for (i = 0; i < w; i++) {+ double mugx2;
+ double mugy2;
+ double mugxgy;
+ mugx2 = col_sums_gx2[0];
+ for (k = 1; k < 8; k++)
+ mugx2 += col_sums_gx2[k];
+ mugy2 = col_sums_gy2[0];
+ for (k = 1; k < 8; k++)
+ mugy2 += col_sums_gy2[k];
+ mugxgy = col_sums_gxgy[0];
+ for (k = 1; k < 8; k++)
+ mugxgy += col_sums_gxgy[k];
+ ssim[(j - 4) * w + i] = (2 * mugxgy + c2) / (mugx2 + mugy2 + c2);
+ if (i + 1 < w) {+ FS_COL_SET(0, -1, 1);
+ FS_COL_ADD(0, 0, 1);
+ FS_COL_SUB(2, -3, 2);
+ FS_COL_SUB(2, 2, 2);
+ FS_COL_HALVE(1, 2);
+ FS_COL_SUB(3, -4, 3);
+ FS_COL_SUB(3, 3, 3);
+ FS_COL_HALVE(2, 3);
+ FS_COL_COPY(3, 4);
+ FS_COL_DOUBLE(4, 5);
+ FS_COL_ADD(4, -4, 5);
+ FS_COL_ADD(4, 3, 5);
+ FS_COL_DOUBLE(5, 6);
+ FS_COL_ADD(5, -3, 6);
+ FS_COL_ADD(5, 2, 6);
+ FS_COL_DOUBLE(6, 7);
+ FS_COL_ADD(6, -2, 7);
+ FS_COL_ADD(6, 1, 7);
+ FS_COL_SET(7, -1, 8);
+ FS_COL_ADD(7, 0, 8);
+ }
+ }
+ }
+ }
+}
+
+#define FS_NLEVELS (4)
+
+/*These weights were derived from the default weights found in Wang's original
+ Matlab implementation: {0.0448, 0.2856, 0.2363, 0.1333}.+ We drop the finest scale and renormalize the rest to sum to 1.*/
+
+static const double FS_WEIGHTS[FS_NLEVELS] = {0.2989654541015625,+ 0.3141326904296875, 0.2473602294921875, 0.1395416259765625};
+
+static double fs_average(fs_ctx *_ctx, int _l) {+ double *ssim;
+ double ret;
+ int w;
+ int h;
+ int i;
+ int j;
+ w = _ctx->level[_l].w;
+ h = _ctx->level[_l].h;
+ ssim = _ctx->level[_l].ssim;
+ ret = 0;
+ for (j = 0; j < h; j++)
+ for (i = 0; i < w; i++)
+ ret += ssim[j * w + i];
+ return pow(ret / (w * h), FS_WEIGHTS[_l]);
+}
+
+static double calc_ssim(const unsigned char *_src, int _systride,
+ const unsigned char *_dst, int _dystride, int _w, int _h) {+ fs_ctx ctx;
+ double ret;
+ int l;
+ ret = 1;
+ fs_ctx_init(&ctx, _w, _h, FS_NLEVELS);
+ fs_downsample_level0(&ctx, _src, _systride, _dst, _dystride, _w, _h);
+ for (l = 0; l < FS_NLEVELS - 1; l++) {+ fs_calc_structure(&ctx, l);
+ ret *= fs_average(&ctx, l);
+ fs_downsample_level(&ctx, l + 1);
+ }
+ fs_calc_structure(&ctx, l);
+ fs_apply_luminance(&ctx, l);
+ ret *= fs_average(&ctx, l);
+ fs_ctx_clear(&ctx);
+ return ret;
+}
+
+static double convert_ssim_db(double _ssim, double _weight) {+ return 10 * (log10(_weight) - log10(_weight - _ssim));
+}
+
+double vpx_calc_fastssim(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
+ double *ssim_y, double *ssim_u, double *ssim_v) {+ double ssimv;
+ vpx_clear_system_state();
+
+ *ssim_y = calc_ssim(source->y_buffer, source->y_stride, dest->y_buffer,
+ dest->y_stride, source->y_crop_width,
+ source->y_crop_height);
+
+ *ssim_u = calc_ssim(source->u_buffer, source->uv_stride, dest->u_buffer,
+ dest->uv_stride, source->uv_crop_width,
+ source->uv_crop_height);
+
+ *ssim_v = calc_ssim(source->v_buffer, source->uv_stride, dest->v_buffer,
+ dest->uv_stride, source->uv_crop_width,
+ source->uv_crop_height);
+ ssimv = (*ssim_y) * .8 + .1 * ((*ssim_u) + (*ssim_v));
+
+ return convert_ssim_db(ssimv, 1.0);
+}
--- /dev/null
+++ b/vpx_dsp/psnrhvs.c
@@ -1,0 +1,223 @@
+/*
+ * Copyright (c) 2010 The WebM project authors. All Rights Reserved.
+ *
+ * Use of this source code is governed by a BSD-style license
+ * that can be found in the LICENSE file in the root of the source
+ * tree. An additional intellectual property rights grant can be found
+ * in the file PATENTS. All contributing project authors may
+ * be found in the AUTHORS file in the root of the source tree.
+ *
+ * This code was originally written by: Gregory Maxwell, at the Daala
+ * project.
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "./vpx_config.h"
+#include "./vpx_dsp_rtcd.h"
+#include "vpx_dsp/ssim.h"
+
+#if !defined(M_PI)
+# define M_PI (3.141592653589793238462643)
+#endif
+#include <string.h>
+
+void od_bin_fdct8x8(tran_low_t *y, int ystride, const int16_t *x, int xstride) {+ (void) xstride;
+ vpx_fdct8x8(x, y, ystride);
+}
+
+/* Normalized inverse quantization matrix for 8x8 DCT at the point of
+ * transparency. This is not the JPEG based matrix from the paper,
+ this one gives a slightly higher MOS agreement.*/
+float csf_y[8][8] = {{1.6193873005, 2.2901594831, 2.08509755623, 1.48366094411,+ 1.00227514334, 0.678296995242, 0.466224900598, 0.3265091542}, {2.2901594831,+ 1.94321815382, 2.04793073064, 1.68731108984, 1.2305666963, 0.868920337363,
+ 0.61280991668, 0.436405793551}, {2.08509755623, 2.04793073064,+ 1.34329019223, 1.09205635862, 0.875748795257, 0.670882927016,
+ 0.501731932449, 0.372504254596}, {1.48366094411, 1.68731108984,+ 1.09205635862, 0.772819797575, 0.605636379554, 0.48309405692,
+ 0.380429446972, 0.295774038565}, {1.00227514334, 1.2305666963,+ 0.875748795257, 0.605636379554, 0.448996256676, 0.352889268808,
+ 0.283006984131, 0.226951348204}, {0.678296995242, 0.868920337363,+ 0.670882927016, 0.48309405692, 0.352889268808, 0.27032073436,
+ 0.215017739696, 0.17408067321}, {0.466224900598, 0.61280991668,+ 0.501731932449, 0.380429446972, 0.283006984131, 0.215017739696,
+ 0.168869545842, 0.136153931001}, {0.3265091542, 0.436405793551,+ 0.372504254596, 0.295774038565, 0.226951348204, 0.17408067321,
+ 0.136153931001, 0.109083846276}};
+float csf_cb420[8][8] = {+ {1.91113096927, 2.46074210438, 1.18284184739, 1.14982565193, 1.05017074788,+ 0.898018824055, 0.74725392039, 0.615105596242}, {2.46074210438,+ 1.58529308355, 1.21363250036, 1.38190029285, 1.33100189972,
+ 1.17428548929, 0.996404342439, 0.830890433625}, {1.18284184739,+ 1.21363250036, 0.978712413627, 1.02624506078, 1.03145147362,
+ 0.960060382087, 0.849823426169, 0.731221236837}, {1.14982565193,+ 1.38190029285, 1.02624506078, 0.861317501629, 0.801821139099,
+ 0.751437590932, 0.685398513368, 0.608694761374}, {1.05017074788,+ 1.33100189972, 1.03145147362, 0.801821139099, 0.676555426187,
+ 0.605503172737, 0.55002013668, 0.495804539034}, {0.898018824055,+ 1.17428548929, 0.960060382087, 0.751437590932, 0.605503172737,
+ 0.514674450957, 0.454353482512, 0.407050308965}, {0.74725392039,+ 0.996404342439, 0.849823426169, 0.685398513368, 0.55002013668,
+ 0.454353482512, 0.389234902883, 0.342353999733}, {0.615105596242,+ 0.830890433625, 0.731221236837, 0.608694761374, 0.495804539034,
+ 0.407050308965, 0.342353999733, 0.295530605237}};
+float csf_cr420[8][8] = {+ {2.03871978502, 2.62502345193, 1.26180942886, 1.11019789803, 1.01397751469,+ 0.867069376285, 0.721500455585, 0.593906509971}, {2.62502345193,+ 1.69112867013, 1.17180569821, 1.3342742857, 1.28513006198,
+ 1.13381474809, 0.962064122248, 0.802254508198}, {1.26180942886,+ 1.17180569821, 0.944981930573, 0.990876405848, 0.995903384143,
+ 0.926972725286, 0.820534991409, 0.706020324706}, {1.11019789803,+ 1.3342742857, 0.990876405848, 0.831632933426, 0.77418706195,
+ 0.725539939514, 0.661776842059, 0.587716619023}, {1.01397751469,+ 1.28513006198, 0.995903384143, 0.77418706195, 0.653238524286,
+ 0.584635025748, 0.531064164893, 0.478717061273}, {0.867069376285,+ 1.13381474809, 0.926972725286, 0.725539939514, 0.584635025748,
+ 0.496936637883, 0.438694579826, 0.393021669543}, {0.721500455585,+ 0.962064122248, 0.820534991409, 0.661776842059, 0.531064164893,
+ 0.438694579826, 0.375820256136, 0.330555063063}, {0.593906509971,+ 0.802254508198, 0.706020324706, 0.587716619023, 0.478717061273,
+ 0.393021669543, 0.330555063063, 0.285345396658}};
+
+static double convert_score_db(double _score, double _weight) {+ return 10 * (log10(255 * 255) - log10(_weight * _score));
+}
+
+static double calc_psnrhvs(const unsigned char *_src, int _systride,
+ const unsigned char *_dst, int _dystride,
+ double _par, int _w, int _h, int _step,
+ float _csf[8][8]) {+ float ret;
+ int16_t dct_s[8 * 8], dct_d[8 * 8];
+ tran_low_t dct_s_coef[8 * 8], dct_d_coef[8 * 8];
+ float mask[8][8];
+ int pixels;
+ int x;
+ int y;
+ (void) _par;
+ ret = pixels = 0;
+ /*In the PSNR-HVS-M paper[1] the authors describe the construction of
+ their masking table as "we have used the quantization table for the
+ color component Y of JPEG [6] that has been also obtained on the
+ basis of CSF. Note that the values in quantization table JPEG have
+ been normalized and then squared." Their CSF matrix (from PSNR-HVS)
+ was also constructed from the JPEG matrices. I can not find any obvious
+ scheme of normalizing to produce their table, but if I multiply their
+ CSF by 0.38857 and square the result I get their masking table.
+ I have no idea where this constant comes from, but deviating from it
+ too greatly hurts MOS agreement.
+
+ [1] Nikolay Ponomarenko, Flavia Silvestri, Karen Egiazarian, Marco Carli,
+ Jaakko Astola, Vladimir Lukin, "On between-coefficient contrast masking
+ of DCT basis functions", CD-ROM Proceedings of the Third
+ International Workshop on Video Processing and Quality Metrics for Consumer
+ Electronics VPQM-07, Scottsdale, Arizona, USA, 25-26 January, 2007, 4 p.*/
+ for (x = 0; x < 8; x++)
+ for (y = 0; y < 8; y++)
+ mask[x][y] = (_csf[x][y] * 0.3885746225901003)
+ * (_csf[x][y] * 0.3885746225901003);
+ for (y = 0; y < _h - 7; y += _step) {+ for (x = 0; x < _w - 7; x += _step) {+ int i;
+ int j;
+ float s_means[4];
+ float d_means[4];
+ float s_vars[4];
+ float d_vars[4];
+ float s_gmean = 0;
+ float d_gmean = 0;
+ float s_gvar = 0;
+ float d_gvar = 0;
+ float s_mask = 0;
+ float d_mask = 0;
+ for (i = 0; i < 4; i++)
+ s_means[i] = d_means[i] = s_vars[i] = d_vars[i] = 0;
+ for (i = 0; i < 8; i++) {+ for (j = 0; j < 8; j++) {+ int sub = ((i & 12) >> 2) + ((j & 12) >> 1);
+ dct_s[i * 8 + j] = _src[(y + i) * _systride + (j + x)];
+ dct_d[i * 8 + j] = _dst[(y + i) * _dystride + (j + x)];
+ s_gmean += dct_s[i * 8 + j];
+ d_gmean += dct_d[i * 8 + j];
+ s_means[sub] += dct_s[i * 8 + j];
+ d_means[sub] += dct_d[i * 8 + j];
+ }
+ }
+ s_gmean /= 64.f;
+ d_gmean /= 64.f;
+ for (i = 0; i < 4; i++)
+ s_means[i] /= 16.f;
+ for (i = 0; i < 4; i++)
+ d_means[i] /= 16.f;
+ for (i = 0; i < 8; i++) {+ for (j = 0; j < 8; j++) {+ int sub = ((i & 12) >> 2) + ((j & 12) >> 1);
+ s_gvar += (dct_s[i * 8 + j] - s_gmean) * (dct_s[i * 8 + j] - s_gmean);
+ d_gvar += (dct_d[i * 8 + j] - d_gmean) * (dct_d[i * 8 + j] - d_gmean);
+ s_vars[sub] += (dct_s[i * 8 + j] - s_means[sub])
+ * (dct_s[i * 8 + j] - s_means[sub]);
+ d_vars[sub] += (dct_d[i * 8 + j] - d_means[sub])
+ * (dct_d[i * 8 + j] - d_means[sub]);
+ }
+ }
+ s_gvar *= 1 / 63.f * 64;
+ d_gvar *= 1 / 63.f * 64;
+ for (i = 0; i < 4; i++)
+ s_vars[i] *= 1 / 15.f * 16;
+ for (i = 0; i < 4; i++)
+ d_vars[i] *= 1 / 15.f * 16;
+ if (s_gvar > 0)
+ s_gvar = (s_vars[0] + s_vars[1] + s_vars[2] + s_vars[3]) / s_gvar;
+ if (d_gvar > 0)
+ d_gvar = (d_vars[0] + d_vars[1] + d_vars[2] + d_vars[3]) / d_gvar;
+ od_bin_fdct8x8(dct_s_coef, 8, dct_s, 8);
+ od_bin_fdct8x8(dct_d_coef, 8, dct_d, 8);
+ for (i = 0; i < 8; i++)
+ for (j = (i == 0); j < 8; j++)
+ s_mask += dct_s_coef[i * 8 + j] * dct_s_coef[i * 8 + j] * mask[i][j];
+ for (i = 0; i < 8; i++)
+ for (j = (i == 0); j < 8; j++)
+ d_mask += dct_d_coef[i * 8 + j] * dct_d_coef[i * 8 + j] * mask[i][j];
+ s_mask = sqrt(s_mask * s_gvar) / 32.f;
+ d_mask = sqrt(d_mask * d_gvar) / 32.f;
+ if (d_mask > s_mask)
+ s_mask = d_mask;
+ for (i = 0; i < 8; i++) {+ for (j = 0; j < 8; j++) {+ float err;
+ err = fabs(dct_s_coef[i * 8 + j] - dct_d_coef[i * 8 + j]);
+ if (i != 0 || j != 0)
+ err = err < s_mask / mask[i][j] ? 0 : err - s_mask / mask[i][j];
+ ret += (err * _csf[i][j]) * (err * _csf[i][j]);
+ pixels++;
+ }
+ }
+ }
+ }
+ ret /= pixels;
+ return ret;
+}
+double vpx_psnrhvs(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
+ double *y_psnrhvs, double *u_psnrhvs, double *v_psnrhvs) {+ double psnrhvs;
+ double par = 1.0;
+ int step = 7;
+ vpx_clear_system_state();
+ *y_psnrhvs = calc_psnrhvs(source->y_buffer, source->y_stride, dest->y_buffer,
+ dest->y_stride, par, source->y_crop_width,
+ source->y_crop_height, step, csf_y);
+
+ *u_psnrhvs = calc_psnrhvs(source->u_buffer, source->uv_stride, dest->u_buffer,
+ dest->uv_stride, par, source->uv_crop_width,
+ source->uv_crop_height, step, csf_cb420);
+
+ *v_psnrhvs = calc_psnrhvs(source->v_buffer, source->uv_stride, dest->v_buffer,
+ dest->uv_stride, par, source->uv_crop_width,
+ source->uv_crop_height, step, csf_cr420);
+ psnrhvs = (*y_psnrhvs) * .8 + .1 * ((*u_psnrhvs) + (*v_psnrhvs));
+
+ return convert_score_db(psnrhvs, 1.0);
+}
--- /dev/null
+++ b/vpx_dsp/ssim.c
@@ -1,0 +1,500 @@
+/*
+ * Copyright (c) 2010 The WebM project authors. All Rights Reserved.
+ *
+ * Use of this source code is governed by a BSD-style license
+ * that can be found in the LICENSE file in the root of the source
+ * tree. An additional intellectual property rights grant can be found
+ * in the file PATENTS. All contributing project authors may
+ * be found in the AUTHORS file in the root of the source tree.
+ */
+
+#include <math.h>
+#include "./vpx_dsp_rtcd.h"
+#include "vpx_dsp/ssim.h"
+#include "vpx_ports/mem.h"
+
+void vpx_ssim_parms_16x16_c(uint8_t *s, int sp, uint8_t *r,
+ int rp, unsigned long *sum_s, unsigned long *sum_r,
+ unsigned long *sum_sq_s, unsigned long *sum_sq_r,
+ unsigned long *sum_sxr) {+ int i, j;
+ for (i = 0; i < 16; i++, s += sp, r += rp) {+ for (j = 0; j < 16; j++) {+ *sum_s += s[j];
+ *sum_r += r[j];
+ *sum_sq_s += s[j] * s[j];
+ *sum_sq_r += r[j] * r[j];
+ *sum_sxr += s[j] * r[j];
+ }
+ }
+}
+void vpx_ssim_parms_8x8_c(uint8_t *s, int sp, uint8_t *r, int rp,
+ unsigned long *sum_s, unsigned long *sum_r,
+ unsigned long *sum_sq_s, unsigned long *sum_sq_r,
+ unsigned long *sum_sxr) {+ int i, j;
+ for (i = 0; i < 8; i++, s += sp, r += rp) {+ for (j = 0; j < 8; j++) {+ *sum_s += s[j];
+ *sum_r += r[j];
+ *sum_sq_s += s[j] * s[j];
+ *sum_sq_r += r[j] * r[j];
+ *sum_sxr += s[j] * r[j];
+ }
+ }
+}
+
+#if CONFIG_VP9_HIGHBITDEPTH
+void vpx_highbd_ssim_parms_8x8_c(uint16_t *s, int sp, uint16_t *r, int rp,
+ uint32_t *sum_s, uint32_t *sum_r,
+ uint32_t *sum_sq_s, uint32_t *sum_sq_r,
+ uint32_t *sum_sxr) {+ int i, j;
+ for (i = 0; i < 8; i++, s += sp, r += rp) {+ for (j = 0; j < 8; j++) {+ *sum_s += s[j];
+ *sum_r += r[j];
+ *sum_sq_s += s[j] * s[j];
+ *sum_sq_r += r[j] * r[j];
+ *sum_sxr += s[j] * r[j];
+ }
+ }
+}
+#endif // CONFIG_VP9_HIGHBITDEPTH
+
+static const int64_t cc1 = 26634; // (64^2*(.01*255)^2
+static const int64_t cc2 = 239708; // (64^2*(.03*255)^2
+
+static double similarity(unsigned long sum_s, unsigned long sum_r,
+ unsigned long sum_sq_s, unsigned long sum_sq_r,
+ unsigned long sum_sxr, int count) {+ int64_t ssim_n, ssim_d;
+ int64_t c1, c2;
+
+ // scale the constants by number of pixels
+ c1 = (cc1 * count * count) >> 12;
+ c2 = (cc2 * count * count) >> 12;
+
+ ssim_n = (2 * sum_s * sum_r + c1) * ((int64_t) 2 * count * sum_sxr -
+ (int64_t) 2 * sum_s * sum_r + c2);
+
+ ssim_d = (sum_s * sum_s + sum_r * sum_r + c1) *
+ ((int64_t)count * sum_sq_s - (int64_t)sum_s * sum_s +
+ (int64_t)count * sum_sq_r - (int64_t) sum_r * sum_r + c2);
+
+ return ssim_n * 1.0 / ssim_d;
+}
+
+static double ssim_8x8(uint8_t *s, int sp, uint8_t *r, int rp) {+ unsigned long sum_s = 0, sum_r = 0, sum_sq_s = 0, sum_sq_r = 0, sum_sxr = 0;
+ vpx_ssim_parms_8x8(s, sp, r, rp, &sum_s, &sum_r, &sum_sq_s, &sum_sq_r,
+ &sum_sxr);
+ return similarity(sum_s, sum_r, sum_sq_s, sum_sq_r, sum_sxr, 64);
+}
+
+#if CONFIG_VP9_HIGHBITDEPTH
+static double highbd_ssim_8x8(uint16_t *s, int sp, uint16_t *r, int rp,
+ unsigned int bd) {+ uint32_t sum_s = 0, sum_r = 0, sum_sq_s = 0, sum_sq_r = 0, sum_sxr = 0;
+ const int oshift = bd - 8;
+ vpx_highbd_ssim_parms_8x8(s, sp, r, rp, &sum_s, &sum_r, &sum_sq_s, &sum_sq_r,
+ &sum_sxr);
+ return similarity(sum_s >> oshift,
+ sum_r >> oshift,
+ sum_sq_s >> (2 * oshift),
+ sum_sq_r >> (2 * oshift),
+ sum_sxr >> (2 * oshift),
+ 64);
+}
+#endif // CONFIG_VP9_HIGHBITDEPTH
+
+// We are using a 8x8 moving window with starting location of each 8x8 window
+// on the 4x4 pixel grid. Such arrangement allows the windows to overlap
+// block boundaries to penalize blocking artifacts.
+double vpx_ssim2(uint8_t *img1, uint8_t *img2, int stride_img1,
+ int stride_img2, int width, int height) {+ int i, j;
+ int samples = 0;
+ double ssim_total = 0;
+
+ // sample point start with each 4x4 location
+ for (i = 0; i <= height - 8;
+ i += 4, img1 += stride_img1 * 4, img2 += stride_img2 * 4) {+ for (j = 0; j <= width - 8; j += 4) {+ double v = ssim_8x8(img1 + j, stride_img1, img2 + j, stride_img2);
+ ssim_total += v;
+ samples++;
+ }
+ }
+ ssim_total /= samples;
+ return ssim_total;
+}
+
+#if CONFIG_VP9_HIGHBITDEPTH
+double vpx_highbd_ssim2(uint8_t *img1, uint8_t *img2, int stride_img1,
+ int stride_img2, int width, int height,
+ unsigned int bd) {+ int i, j;
+ int samples = 0;
+ double ssim_total = 0;
+
+ // sample point start with each 4x4 location
+ for (i = 0; i <= height - 8;
+ i += 4, img1 += stride_img1 * 4, img2 += stride_img2 * 4) {+ for (j = 0; j <= width - 8; j += 4) {+ double v = highbd_ssim_8x8(CONVERT_TO_SHORTPTR(img1 + j), stride_img1,
+ CONVERT_TO_SHORTPTR(img2 + j), stride_img2,
+ bd);
+ ssim_total += v;
+ samples++;
+ }
+ }
+ ssim_total /= samples;
+ return ssim_total;
+}
+#endif // CONFIG_VP9_HIGHBITDEPTH
+
+double vpx_calc_ssim(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
+ double *weight) {+ double a, b, c;
+ double ssimv;
+
+ a = vpx_ssim2(source->y_buffer, dest->y_buffer,
+ source->y_stride, dest->y_stride,
+ source->y_crop_width, source->y_crop_height);
+
+ b = vpx_ssim2(source->u_buffer, dest->u_buffer,
+ source->uv_stride, dest->uv_stride,
+ source->uv_crop_width, source->uv_crop_height);
+
+ c = vpx_ssim2(source->v_buffer, dest->v_buffer,
+ source->uv_stride, dest->uv_stride,
+ source->uv_crop_width, source->uv_crop_height);
+
+ ssimv = a * .8 + .1 * (b + c);
+
+ *weight = 1;
+
+ return ssimv;
+}
+
+double vpx_calc_ssimg(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
+ double *ssim_y, double *ssim_u, double *ssim_v) {+ double ssim_all = 0;
+ double a, b, c;
+
+ a = vpx_ssim2(source->y_buffer, dest->y_buffer,
+ source->y_stride, dest->y_stride,
+ source->y_crop_width, source->y_crop_height);
+
+ b = vpx_ssim2(source->u_buffer, dest->u_buffer,
+ source->uv_stride, dest->uv_stride,
+ source->uv_crop_width, source->uv_crop_height);
+
+ c = vpx_ssim2(source->v_buffer, dest->v_buffer,
+ source->uv_stride, dest->uv_stride,
+ source->uv_crop_width, source->uv_crop_height);
+ *ssim_y = a;
+ *ssim_u = b;
+ *ssim_v = c;
+ ssim_all = (a * 4 + b + c) / 6;
+
+ return ssim_all;
+}
+
+// traditional ssim as per: http://en.wikipedia.org/wiki/Structural_similarity
+//
+// Re working out the math ->
+//
+// ssim(x,y) = (2*mean(x)*mean(y) + c1)*(2*cov(x,y)+c2) /
+// ((mean(x)^2+mean(y)^2+c1)*(var(x)+var(y)+c2))
+//
+// mean(x) = sum(x) / n
+//
+// cov(x,y) = (n*sum(xi*yi)-sum(x)*sum(y))/(n*n)
+//
+// var(x) = (n*sum(xi*xi)-sum(xi)*sum(xi))/(n*n)
+//
+// ssim(x,y) =
+// (2*sum(x)*sum(y)/(n*n) + c1)*(2*(n*sum(xi*yi)-sum(x)*sum(y))/(n*n)+c2) /
+// (((sum(x)*sum(x)+sum(y)*sum(y))/(n*n) +c1) *
+// ((n*sum(xi*xi) - sum(xi)*sum(xi))/(n*n)+
+// (n*sum(yi*yi) - sum(yi)*sum(yi))/(n*n)+c2)))
+//
+// factoring out n*n
+//
+// ssim(x,y) =
+// (2*sum(x)*sum(y) + n*n*c1)*(2*(n*sum(xi*yi)-sum(x)*sum(y))+n*n*c2) /
+// (((sum(x)*sum(x)+sum(y)*sum(y)) + n*n*c1) *
+// (n*sum(xi*xi)-sum(xi)*sum(xi)+n*sum(yi*yi)-sum(yi)*sum(yi)+n*n*c2))
+//
+// Replace c1 with n*n * c1 for the final step that leads to this code:
+// The final step scales by 12 bits so we don't lose precision in the constants.
+
+double ssimv_similarity(Ssimv *sv, int64_t n) {+ // Scale the constants by number of pixels.
+ const int64_t c1 = (cc1 * n * n) >> 12;
+ const int64_t c2 = (cc2 * n * n) >> 12;
+
+ const double l = 1.0 * (2 * sv->sum_s * sv->sum_r + c1) /
+ (sv->sum_s * sv->sum_s + sv->sum_r * sv->sum_r + c1);
+
+ // Since these variables are unsigned sums, convert to double so
+ // math is done in double arithmetic.
+ const double v = (2.0 * n * sv->sum_sxr - 2 * sv->sum_s * sv->sum_r + c2)
+ / (n * sv->sum_sq_s - sv->sum_s * sv->sum_s + n * sv->sum_sq_r
+ - sv->sum_r * sv->sum_r + c2);
+
+ return l * v;
+}
+
+// The first term of the ssim metric is a luminance factor.
+//
+// (2*mean(x)*mean(y) + c1)/ (mean(x)^2+mean(y)^2+c1)
+//
+// This luminance factor is super sensitive to the dark side of luminance
+// values and completely insensitive on the white side. check out 2 sets
+// (1,3) and (250,252) the term gives ( 2*1*3/(1+9) = .60
+// 2*250*252/ (250^2+252^2) => .99999997
+//
+// As a result in this tweaked version of the calculation in which the
+// luminance is taken as percentage off from peak possible.
+//
+// 255 * 255 - (sum_s - sum_r) / count * (sum_s - sum_r) / count
+//
+double ssimv_similarity2(Ssimv *sv, int64_t n) {+ // Scale the constants by number of pixels.
+ const int64_t c1 = (cc1 * n * n) >> 12;
+ const int64_t c2 = (cc2 * n * n) >> 12;
+
+ const double mean_diff = (1.0 * sv->sum_s - sv->sum_r) / n;
+ const double l = (255 * 255 - mean_diff * mean_diff + c1) / (255 * 255 + c1);
+
+ // Since these variables are unsigned, sums convert to double so
+ // math is done in double arithmetic.
+ const double v = (2.0 * n * sv->sum_sxr - 2 * sv->sum_s * sv->sum_r + c2)
+ / (n * sv->sum_sq_s - sv->sum_s * sv->sum_s +
+ n * sv->sum_sq_r - sv->sum_r * sv->sum_r + c2);
+
+ return l * v;
+}
+void ssimv_parms(uint8_t *img1, int img1_pitch, uint8_t *img2, int img2_pitch,
+ Ssimv *sv) {+ vpx_ssim_parms_8x8(img1, img1_pitch, img2, img2_pitch,
+ &sv->sum_s, &sv->sum_r, &sv->sum_sq_s, &sv->sum_sq_r,
+ &sv->sum_sxr);
+}
+
+double vpx_get_ssim_metrics(uint8_t *img1, int img1_pitch,
+ uint8_t *img2, int img2_pitch,
+ int width, int height,
+ Ssimv *sv2, Metrics *m,
+ int do_inconsistency) {+ double dssim_total = 0;
+ double ssim_total = 0;
+ double ssim2_total = 0;
+ double inconsistency_total = 0;
+ int i, j;
+ int c = 0;
+ double norm;
+ double old_ssim_total = 0;
+ vpx_clear_system_state();
+ // We can sample points as frequently as we like start with 1 per 4x4.
+ for (i = 0; i < height; i += 4,
+ img1 += img1_pitch * 4, img2 += img2_pitch * 4) {+ for (j = 0; j < width; j += 4, ++c) {+ Ssimv sv = {0};+ double ssim;
+ double ssim2;
+ double dssim;
+ uint32_t var_new;
+ uint32_t var_old;
+ uint32_t mean_new;
+ uint32_t mean_old;
+ double ssim_new;
+ double ssim_old;
+
+ // Not sure there's a great way to handle the edge pixels
+ // in ssim when using a window. Seems biased against edge pixels
+ // however you handle this. This uses only samples that are
+ // fully in the frame.
+ if (j + 8 <= width && i + 8 <= height) {+ ssimv_parms(img1 + j, img1_pitch, img2 + j, img2_pitch, &sv);
+ }
+
+ ssim = ssimv_similarity(&sv, 64);
+ ssim2 = ssimv_similarity2(&sv, 64);
+
+ sv.ssim = ssim2;
+
+ // dssim is calculated to use as an actual error metric and
+ // is scaled up to the same range as sum square error.
+ // Since we are subsampling every 16th point maybe this should be
+ // *16 ?
+ dssim = 255 * 255 * (1 - ssim2) / 2;
+
+ // Here I introduce a new error metric: consistency-weighted
+ // SSIM-inconsistency. This metric isolates frames where the
+ // SSIM 'suddenly' changes, e.g. if one frame in every 8 is much
+ // sharper or blurrier than the others. Higher values indicate a
+ // temporally inconsistent SSIM. There are two ideas at work:
+ //
+ // 1) 'SSIM-inconsistency': the total inconsistency value
+ // reflects how much SSIM values are changing between this
+ // source / reference frame pair and the previous pair.
+ //
+ // 2) 'consistency-weighted': weights de-emphasize areas in the
+ // frame where the scene content has changed. Changes in scene
+ // content are detected via changes in local variance and local
+ // mean.
+ //
+ // Thus the overall measure reflects how inconsistent the SSIM
+ // values are, over consistent regions of the frame.
+ //
+ // The metric has three terms:
+ //
+ // term 1 -> uses change in scene Variance to weight error score
+ // 2 * var(Fi)*var(Fi-1) / (var(Fi)^2+var(Fi-1)^2)
+ // larger changes from one frame to the next mean we care
+ // less about consistency.
+ //
+ // term 2 -> uses change in local scene luminance to weight error
+ // 2 * avg(Fi)*avg(Fi-1) / (avg(Fi)^2+avg(Fi-1)^2)
+ // larger changes from one frame to the next mean we care
+ // less about consistency.
+ //
+ // term3 -> measures inconsistency in ssim scores between frames
+ // 1 - ( 2 * ssim(Fi)*ssim(Fi-1)/(ssim(Fi)^2+sssim(Fi-1)^2).
+ //
+ // This term compares the ssim score for the same location in 2
+ // subsequent frames.
+ var_new = sv.sum_sq_s - sv.sum_s * sv.sum_s / 64;
+ var_old = sv2[c].sum_sq_s - sv2[c].sum_s * sv2[c].sum_s / 64;
+ mean_new = sv.sum_s;
+ mean_old = sv2[c].sum_s;
+ ssim_new = sv.ssim;
+ ssim_old = sv2[c].ssim;
+
+ if (do_inconsistency) {+ // We do the metric once for every 4x4 block in the image. Since
+ // we are scaling the error to SSE for use in a psnr calculation
+ // 1.0 = 4x4x255x255 the worst error we can possibly have.
+ static const double kScaling = 4. * 4 * 255 * 255;
+
+ // The constants have to be non 0 to avoid potential divide by 0
+ // issues other than that they affect kind of a weighting between
+ // the terms. No testing of what the right terms should be has been
+ // done.
+ static const double c1 = 1, c2 = 1, c3 = 1;
+
+ // This measures how much consistent variance is in two consecutive
+ // source frames. 1.0 means they have exactly the same variance.
+ const double variance_term = (2.0 * var_old * var_new + c1) /
+ (1.0 * var_old * var_old + 1.0 * var_new * var_new + c1);
+
+ // This measures how consistent the local mean are between two
+ // consecutive frames. 1.0 means they have exactly the same mean.
+ const double mean_term = (2.0 * mean_old * mean_new + c2) /
+ (1.0 * mean_old * mean_old + 1.0 * mean_new * mean_new + c2);
+
+ // This measures how consistent the ssims of two
+ // consecutive frames is. 1.0 means they are exactly the same.
+ double ssim_term = pow((2.0 * ssim_old * ssim_new + c3) /
+ (ssim_old * ssim_old + ssim_new * ssim_new + c3),
+ 5);
+
+ double this_inconsistency;
+
+ // Floating point math sometimes makes this > 1 by a tiny bit.
+ // We want the metric to scale between 0 and 1.0 so we can convert
+ // it to an snr scaled value.
+ if (ssim_term > 1)
+ ssim_term = 1;
+
+ // This converts the consistency metric to an inconsistency metric
+ // ( so we can scale it like psnr to something like sum square error.
+ // The reason for the variance and mean terms is the assumption that
+ // if there are big changes in the source we shouldn't penalize
+ // inconsistency in ssim scores a bit less as it will be less visible
+ // to the user.
+ this_inconsistency = (1 - ssim_term) * variance_term * mean_term;
+
+ this_inconsistency *= kScaling;
+ inconsistency_total += this_inconsistency;
+ }
+ sv2[c] = sv;
+ ssim_total += ssim;
+ ssim2_total += ssim2;
+ dssim_total += dssim;
+
+ old_ssim_total += ssim_old;
+ }
+ old_ssim_total += 0;
+ }
+
+ norm = 1. / (width / 4) / (height / 4);
+ ssim_total *= norm;
+ ssim2_total *= norm;
+ m->ssim2 = ssim2_total;
+ m->ssim = ssim_total;
+ if (old_ssim_total == 0)
+ inconsistency_total = 0;
+
+ m->ssimc = inconsistency_total;
+
+ m->dssim = dssim_total;
+ return inconsistency_total;
+}
+
+
+#if CONFIG_VP9_HIGHBITDEPTH
+double vpx_highbd_calc_ssim(YV12_BUFFER_CONFIG *source,
+ YV12_BUFFER_CONFIG *dest,
+ double *weight, unsigned int bd) {+ double a, b, c;
+ double ssimv;
+
+ a = vpx_highbd_ssim2(source->y_buffer, dest->y_buffer,
+ source->y_stride, dest->y_stride,
+ source->y_crop_width, source->y_crop_height, bd);
+
+ b = vpx_highbd_ssim2(source->u_buffer, dest->u_buffer,
+ source->uv_stride, dest->uv_stride,
+ source->uv_crop_width, source->uv_crop_height, bd);
+
+ c = vpx_highbd_ssim2(source->v_buffer, dest->v_buffer,
+ source->uv_stride, dest->uv_stride,
+ source->uv_crop_width, source->uv_crop_height, bd);
+
+ ssimv = a * .8 + .1 * (b + c);
+
+ *weight = 1;
+
+ return ssimv;
+}
+
+double vpx_highbd_calc_ssimg(YV12_BUFFER_CONFIG *source,
+ YV12_BUFFER_CONFIG *dest, double *ssim_y,
+ double *ssim_u, double *ssim_v, unsigned int bd) {+ double ssim_all = 0;
+ double a, b, c;
+
+ a = vpx_highbd_ssim2(source->y_buffer, dest->y_buffer,
+ source->y_stride, dest->y_stride,
+ source->y_crop_width, source->y_crop_height, bd);
+
+ b = vpx_highbd_ssim2(source->u_buffer, dest->u_buffer,
+ source->uv_stride, dest->uv_stride,
+ source->uv_crop_width, source->uv_crop_height, bd);
+
+ c = vpx_highbd_ssim2(source->v_buffer, dest->v_buffer,
+ source->uv_stride, dest->uv_stride,
+ source->uv_crop_width, source->uv_crop_height, bd);
+ *ssim_y = a;
+ *ssim_u = b;
+ *ssim_v = c;
+ ssim_all = (a * 4 + b + c) / 6;
+
+ return ssim_all;
+}
+#endif // CONFIG_VP9_HIGHBITDEPTH
--- /dev/null
+++ b/vpx_dsp/ssim.h
@@ -1,0 +1,105 @@
+/*
+ * Copyright (c) 2014 The WebM project authors. All Rights Reserved.
+ *
+ * Use of this source code is governed by a BSD-style license
+ * that can be found in the LICENSE file in the root of the source
+ * tree. An additional intellectual property rights grant can be found
+ * in the file PATENTS. All contributing project authors may
+ * be found in the AUTHORS file in the root of the source tree.
+ */
+
+#ifndef VPX_ENCODER_VP9_SSIM_H_
+#define VPX_ENCODER_VP9_SSIM_H_
+
+#ifdef __cplusplus
+extern "C" {+#endif
+
+#include "./vpx_config.h"
+#include "vpx_scale/yv12config.h"
+
+// TODO(aconverse): Unify vp8/vp9_clear_system_state
+#if ARCH_X86 || ARCH_X86_64
+void vpx_reset_mmx_state(void);
+#define vpx_clear_system_state() vpx_reset_mmx_state()
+#else
+#define vpx_clear_system_state()
+#endif
+
+// metrics used for calculating ssim, ssim2, dssim, and ssimc
+typedef struct {+ // source sum ( over 8x8 region )
+ uint64_t sum_s;
+
+ // reference sum (over 8x8 region )
+ uint64_t sum_r;
+
+ // source sum squared ( over 8x8 region )
+ uint64_t sum_sq_s;
+
+ // reference sum squared (over 8x8 region )
+ uint64_t sum_sq_r;
+
+ // sum of source times reference (over 8x8 region)
+ uint64_t sum_sxr;
+
+ // calculated ssim score between source and reference
+ double ssim;
+} Ssimv;
+
+// metrics collected on a frame basis
+typedef struct {+ // ssim consistency error metric ( see code for explanation )
+ double ssimc;
+
+ // standard ssim
+ double ssim;
+
+ // revised ssim ( see code for explanation)
+ double ssim2;
+
+ // ssim restated as an error metric like sse
+ double dssim;
+
+ // dssim converted to decibels
+ double dssimd;
+
+ // ssimc converted to decibels
+ double ssimcd;
+} Metrics;
+
+double vpx_get_ssim_metrics(uint8_t *img1, int img1_pitch, uint8_t *img2,
+ int img2_pitch, int width, int height, Ssimv *sv2,
+ Metrics *m, int do_inconsistency);
+
+double vpx_calc_ssim(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
+ double *weight);
+
+double vpx_calc_ssimg(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
+ double *ssim_y, double *ssim_u, double *ssim_v);
+
+double vpx_calc_fastssim(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
+ double *ssim_y, double *ssim_u, double *ssim_v);
+
+double vpx_psnrhvs(YV12_BUFFER_CONFIG *source, YV12_BUFFER_CONFIG *dest,
+ double *ssim_y, double *ssim_u, double *ssim_v);
+
+#if CONFIG_VP9_HIGHBITDEPTH
+double vpx_highbd_calc_ssim(YV12_BUFFER_CONFIG *source,
+ YV12_BUFFER_CONFIG *dest,
+ double *weight,
+ unsigned int bd);
+
+double vpx_highbd_calc_ssimg(YV12_BUFFER_CONFIG *source,
+ YV12_BUFFER_CONFIG *dest,
+ double *ssim_y,
+ double *ssim_u,
+ double *ssim_v,
+ unsigned int bd);
+#endif // CONFIG_VP9_HIGHBITDEPTH
+
+#ifdef __cplusplus
+} // extern "C"
+#endif
+
+#endif // VPX_ENCODER_VP9_SSIM_H_
--- a/vpx_dsp/vpx_dsp.mk
+++ b/vpx_dsp/vpx_dsp.mk
@@ -22,6 +22,10 @@
DSP_SRCS-yes += bitwriter.c
DSP_SRCS-yes += bitwriter_buffer.c
DSP_SRCS-yes += bitwriter_buffer.h
+DSP_SRCS-$(CONFIG_INTERNAL_STATS) += ssim.c
+DSP_SRCS-$(CONFIG_INTERNAL_STATS) += ssim.h
+DSP_SRCS-$(CONFIG_INTERNAL_STATS) += psnrhvs.c
+DSP_SRCS-$(CONFIG_INTERNAL_STATS) += fastssim.c
endif
ifeq ($(CONFIG_DECODERS),yes)
@@ -294,6 +298,10 @@
DSP_SRCS-$(HAVE_SSE2) += x86/variance_sse2.c # Contains SSE2 and SSSE3
DSP_SRCS-$(HAVE_AVX2) += x86/variance_avx2.c
DSP_SRCS-$(HAVE_AVX2) += x86/variance_impl_avx2.c
+
+ifeq ($(ARCH_X86_64),yes)
+DSP_SRCS-$(HAVE_SSE2) += x86/ssim_opt_x86_64.asm
+endif # ARCH_X86_64
ifeq ($(CONFIG_USE_X86INC),yes)
DSP_SRCS-$(HAVE_SSE2) += x86/subpel_variance_sse2.asm # Contains SSE2 and SSSE3
--- a/vpx_dsp/vpx_dsp_rtcd_defs.pl
+++ b/vpx_dsp/vpx_dsp_rtcd_defs.pl
@@ -990,6 +990,17 @@
add_proto qw/void vpx_sad4x4x4d/, "const uint8_t *src_ptr, int src_stride, const uint8_t * const ref_ptr[], int ref_stride, uint32_t *sad_array";
specialize qw/vpx_sad4x4x4d msa/, "$sse_x86inc";
+#
+# Structured Similarity (SSIM)
+#
+if (vpx_config("CONFIG_INTERNAL_STATS") eq "yes") {+ add_proto qw/void vpx_ssim_parms_8x8/, "uint8_t *s, int sp, uint8_t *r, int rp, unsigned long *sum_s, unsigned long *sum_r, unsigned long *sum_sq_s, unsigned long *sum_sq_r, unsigned long *sum_sxr";
+ specialize qw/vpx_ssim_parms_8x8/, "$sse2_x86_64";
+
+ add_proto qw/void vpx_ssim_parms_16x16/, "uint8_t *s, int sp, uint8_t *r, int rp, unsigned long *sum_s, unsigned long *sum_r, unsigned long *sum_sq_s, unsigned long *sum_sq_r, unsigned long *sum_sxr";
+ specialize qw/vpx_ssim_parms_16x16/, "$sse2_x86_64";
+}
+
if (vpx_config("CONFIG_VP9_HIGHBITDEPTH") eq "yes") {#
# Block subtraction
@@ -1176,6 +1187,13 @@
add_proto qw/void vpx_highbd_sad4x4x4d/, "const uint8_t *src_ptr, int src_stride, const uint8_t* const ref_ptr[], int ref_stride, uint32_t *sad_array";
specialize qw/vpx_highbd_sad4x4x4d/, "$sse2_x86inc";
+ #
+ # Structured Similarity (SSIM)
+ #
+ if (vpx_config("CONFIG_INTERNAL_STATS") eq "yes") {+ add_proto qw/void vpx_highbd_ssim_parms_8x8/, "uint16_t *s, int sp, uint16_t *r, int rp, uint32_t *sum_s, uint32_t *sum_r, uint32_t *sum_sq_s, uint32_t *sum_sq_r, uint32_t *sum_sxr";
+ specialize qw/vpx_highbd_ssim_parms_8x8/;
+ }
} # CONFIG_VP9_HIGHBITDEPTH
} # CONFIG_ENCODERS
--- /dev/null
+++ b/vpx_dsp/x86/ssim_opt_x86_64.asm
@@ -1,0 +1,216 @@
+;
+; Copyright (c) 2010 The WebM project authors. All Rights Reserved.
+;
+; Use of this source code is governed by a BSD-style license
+; that can be found in the LICENSE file in the root of the source
+; tree. An additional intellectual property rights grant can be found
+; in the file PATENTS. All contributing project authors may
+; be found in the AUTHORS file in the root of the source tree.
+;
+
+%include "vpx_ports/x86_abi_support.asm"
+
+; tabulate_ssim - sums sum_s,sum_r,sum_sq_s,sum_sq_r, sum_sxr
+%macro TABULATE_SSIM 0
+ paddusw xmm15, xmm3 ; sum_s
+ paddusw xmm14, xmm4 ; sum_r
+ movdqa xmm1, xmm3
+ pmaddwd xmm1, xmm1
+ paddd xmm13, xmm1 ; sum_sq_s
+ movdqa xmm2, xmm4
+ pmaddwd xmm2, xmm2
+ paddd xmm12, xmm2 ; sum_sq_r
+ pmaddwd xmm3, xmm4
+ paddd xmm11, xmm3 ; sum_sxr
+%endmacro
+
+; Sum across the register %1 starting with q words
+%macro SUM_ACROSS_Q 1
+ movdqa xmm2,%1
+ punpckldq %1,xmm0
+ punpckhdq xmm2,xmm0
+ paddq %1,xmm2
+ movdqa xmm2,%1
+ punpcklqdq %1,xmm0
+ punpckhqdq xmm2,xmm0
+ paddq %1,xmm2
+%endmacro
+
+; Sum across the register %1 starting with q words
+%macro SUM_ACROSS_W 1
+ movdqa xmm1, %1
+ punpcklwd %1,xmm0
+ punpckhwd xmm1,xmm0
+ paddd %1, xmm1
+ SUM_ACROSS_Q %1
+%endmacro
+;void ssim_parms_sse2(
+; unsigned char *s,
+; int sp,
+; unsigned char *r,
+; int rp
+; unsigned long *sum_s,
+; unsigned long *sum_r,
+; unsigned long *sum_sq_s,
+; unsigned long *sum_sq_r,
+; unsigned long *sum_sxr);
+;
+; TODO: Use parm passing through structure, probably don't need the pxors
+; ( calling app will initialize to 0 ) could easily fit everything in sse2
+; without too much hastle, and can probably do better estimates with psadw
+; or pavgb At this point this is just meant to be first pass for calculating
+; all the parms needed for 16x16 ssim so we can play with dssim as distortion
+; in mode selection code.
+global sym(vpx_ssim_parms_16x16_sse2) PRIVATE
+sym(vpx_ssim_parms_16x16_sse2):
+ push rbp
+ mov rbp, rsp
+ SHADOW_ARGS_TO_STACK 9
+ SAVE_XMM 15
+ push rsi
+ push rdi
+ ; end prolog
+
+ mov rsi, arg(0) ;s
+ mov rcx, arg(1) ;sp
+ mov rdi, arg(2) ;r
+ mov rax, arg(3) ;rp
+
+ pxor xmm0, xmm0
+ pxor xmm15,xmm15 ;sum_s
+ pxor xmm14,xmm14 ;sum_r
+ pxor xmm13,xmm13 ;sum_sq_s
+ pxor xmm12,xmm12 ;sum_sq_r
+ pxor xmm11,xmm11 ;sum_sxr
+
+ mov rdx, 16 ;row counter
+.NextRow:
+
+ ;grab source and reference pixels
+ movdqu xmm5, [rsi]
+ movdqu xmm6, [rdi]
+ movdqa xmm3, xmm5
+ movdqa xmm4, xmm6
+ punpckhbw xmm3, xmm0 ; high_s
+ punpckhbw xmm4, xmm0 ; high_r
+
+ TABULATE_SSIM
+
+ movdqa xmm3, xmm5
+ movdqa xmm4, xmm6
+ punpcklbw xmm3, xmm0 ; low_s
+ punpcklbw xmm4, xmm0 ; low_r
+
+ TABULATE_SSIM
+
+ add rsi, rcx ; next s row
+ add rdi, rax ; next r row
+
+ dec rdx ; counter
+ jnz .NextRow
+
+ SUM_ACROSS_W xmm15
+ SUM_ACROSS_W xmm14
+ SUM_ACROSS_Q xmm13
+ SUM_ACROSS_Q xmm12
+ SUM_ACROSS_Q xmm11
+
+ mov rdi,arg(4)
+ movd [rdi], xmm15;
+ mov rdi,arg(5)
+ movd [rdi], xmm14;
+ mov rdi,arg(6)
+ movd [rdi], xmm13;
+ mov rdi,arg(7)
+ movd [rdi], xmm12;
+ mov rdi,arg(8)
+ movd [rdi], xmm11;
+
+ ; begin epilog
+ pop rdi
+ pop rsi
+ RESTORE_XMM
+ UNSHADOW_ARGS
+ pop rbp
+ ret
+
+;void ssim_parms_sse2(
+; unsigned char *s,
+; int sp,
+; unsigned char *r,
+; int rp
+; unsigned long *sum_s,
+; unsigned long *sum_r,
+; unsigned long *sum_sq_s,
+; unsigned long *sum_sq_r,
+; unsigned long *sum_sxr);
+;
+; TODO: Use parm passing through structure, probably don't need the pxors
+; ( calling app will initialize to 0 ) could easily fit everything in sse2
+; without too much hastle, and can probably do better estimates with psadw
+; or pavgb At this point this is just meant to be first pass for calculating
+; all the parms needed for 16x16 ssim so we can play with dssim as distortion
+; in mode selection code.
+global sym(vpx_ssim_parms_8x8_sse2) PRIVATE
+sym(vpx_ssim_parms_8x8_sse2):
+ push rbp
+ mov rbp, rsp
+ SHADOW_ARGS_TO_STACK 9
+ SAVE_XMM 15
+ push rsi
+ push rdi
+ ; end prolog
+
+ mov rsi, arg(0) ;s
+ mov rcx, arg(1) ;sp
+ mov rdi, arg(2) ;r
+ mov rax, arg(3) ;rp
+
+ pxor xmm0, xmm0
+ pxor xmm15,xmm15 ;sum_s
+ pxor xmm14,xmm14 ;sum_r
+ pxor xmm13,xmm13 ;sum_sq_s
+ pxor xmm12,xmm12 ;sum_sq_r
+ pxor xmm11,xmm11 ;sum_sxr
+
+ mov rdx, 8 ;row counter
+.NextRow:
+
+ ;grab source and reference pixels
+ movq xmm3, [rsi]
+ movq xmm4, [rdi]
+ punpcklbw xmm3, xmm0 ; low_s
+ punpcklbw xmm4, xmm0 ; low_r
+
+ TABULATE_SSIM
+
+ add rsi, rcx ; next s row
+ add rdi, rax ; next r row
+
+ dec rdx ; counter
+ jnz .NextRow
+
+ SUM_ACROSS_W xmm15
+ SUM_ACROSS_W xmm14
+ SUM_ACROSS_Q xmm13
+ SUM_ACROSS_Q xmm12
+ SUM_ACROSS_Q xmm11
+
+ mov rdi,arg(4)
+ movd [rdi], xmm15;
+ mov rdi,arg(5)
+ movd [rdi], xmm14;
+ mov rdi,arg(6)
+ movd [rdi], xmm13;
+ mov rdi,arg(7)
+ movd [rdi], xmm12;
+ mov rdi,arg(8)
+ movd [rdi], xmm11;
+
+ ; begin epilog
+ pop rdi
+ pop rsi
+ RESTORE_XMM
+ UNSHADOW_ARGS
+ pop rbp
+ ret