1 
2 /*!
3  *************************************************************************************
4  * \file img_dist_ssim.c
5  *
6  * \brief
7  *    Compute structural similarity (SSIM) index using the encoded image and the reference image
8  *
9  * \author
10  *    Main contributors (see contributors.h for copyright, address and affiliation details)
11  *     - Woo-Shik Kim                    <wooshik.kim@usc.edu>
12  *     - Zhen Li                         <zli@dolby.com>
13  *     - Alexis Michael Tourapis         <alexismt@ieee.org>
14  *************************************************************************************
15  */
16 #include "contributors.h"
17 #include "global.h"
18 #include "img_distortion.h"
19 #include "enc_statistics.h"
20 
21 //#define UNBIASED_VARIANCE // unbiased estimation of the variance
22 
compute_ssim(VideoParameters * p_Vid,InputParameters * p_Inp,imgpel ** refImg,imgpel ** encImg,int height,int width,int win_height,int win_width,int comp)23 float compute_ssim (VideoParameters *p_Vid, InputParameters *p_Inp, imgpel **refImg, imgpel **encImg, int height, int width, int win_height, int win_width, int comp)
24 {
25 
26   static const float K1 = 0.01f, K2 = 0.03f;
27   float max_pix_value_sqd;
28   float C1, C2;
29   float win_pixels = (float) (win_width * win_height);
30 #ifdef UNBIASED_VARIANCE
31   float win_pixels_bias = win_pixels - 1;
32 #else
33   float win_pixels_bias = win_pixels;
34 #endif
35   float mb_ssim, meanOrg, meanEnc;
36   float varOrg, varEnc, covOrgEnc;
37   int imeanOrg, imeanEnc, ivarOrg, ivarEnc, icovOrgEnc;
38   float cur_distortion = 0.0;
39   int i, j, n, m, win_cnt = 0;
40   int overlapSize = p_Inp->SSIMOverlapSize;
41 
42   max_pix_value_sqd = (float) (p_Vid->max_pel_value_comp[comp] * p_Vid->max_pel_value_comp[comp]);
43   C1 = K1 * K1 * max_pix_value_sqd;
44   C2 = K2 * K2 * max_pix_value_sqd;
45 
46   for (j = 0; j <= height - win_height; j += overlapSize)
47   {
48     for (i = 0; i <= width - win_width; i += overlapSize)
49     {
50       imeanOrg = 0;
51       imeanEnc = 0;
52       ivarOrg  = 0;
53       ivarEnc  = 0;
54       icovOrgEnc = 0;
55 
56       for ( n = j;n < j + win_height;n ++)
57       {
58         for (m = i;m < i + win_width;m ++)
59         {
60           imeanOrg   += refImg[n][m];
61           imeanEnc   += encImg[n][m];
62           ivarOrg    += refImg[n][m] * refImg[n][m];
63           ivarEnc    += encImg[n][m] * encImg[n][m];
64           icovOrgEnc += refImg[n][m] * encImg[n][m];
65         }
66       }
67 
68       meanOrg = (float) imeanOrg / win_pixels;
69       meanEnc = (float) imeanEnc / win_pixels;
70 
71       varOrg    = ((float) ivarOrg - ((float) imeanOrg) * meanOrg) / win_pixels_bias;
72       varEnc    = ((float) ivarEnc - ((float) imeanEnc) * meanEnc) / win_pixels_bias;
73       covOrgEnc = ((float) icovOrgEnc - ((float) imeanOrg) * meanEnc) / win_pixels_bias;
74 
75       mb_ssim  = (float) ((2.0 * meanOrg * meanEnc + C1) * (2.0 * covOrgEnc + C2));
76       mb_ssim /= (float) (meanOrg * meanOrg + meanEnc * meanEnc + C1) * (varOrg + varEnc + C2);
77 
78       cur_distortion += mb_ssim;
79       win_cnt++;
80     }
81   }
82 
83   cur_distortion /= (float) win_cnt;
84 
85   if (cur_distortion >= 1.0 && cur_distortion < 1.01) // avoid float accuracy problem at very low QP(e.g.2)
86     cur_distortion = 1.0;
87 
88   return cur_distortion;
89 }
90 
91 /*!
92  ************************************************************************
93  * \brief
94  *    Find SSIM for all three components
95  ************************************************************************
96  */
find_ssim(VideoParameters * p_Vid,InputParameters * p_Inp,ImageStructure * ref,ImageStructure * src,DistMetric * metricSSIM)97 void find_ssim (VideoParameters *p_Vid, InputParameters *p_Inp, ImageStructure *ref, ImageStructure *src, DistMetric *metricSSIM)
98 {
99   DistortionParams *p_Dist = p_Vid->p_Dist;
100   FrameFormat *format = &ref->format;
101 
102   metricSSIM->value[0] = compute_ssim (p_Vid, p_Inp, ref->data[0], src->data[0], format->height[0], format->width[0], BLOCK_SIZE_8x8, BLOCK_SIZE_8x8, 0);
103   // Chroma.
104   if (format->yuv_format != YUV400)
105   {
106     metricSSIM->value[1]  = compute_ssim (p_Vid, p_Inp, ref->data[1], src->data[1], format->height[1], format->width[1], p_Vid->mb_cr_size_y, p_Vid->mb_cr_size_x, 1);
107     metricSSIM->value[2]  = compute_ssim (p_Vid, p_Inp, ref->data[2], src->data[2], format->height[1], format->width[1], p_Vid->mb_cr_size_y, p_Vid->mb_cr_size_x, 2);
108   }
109 
110   {
111     accumulate_average(metricSSIM,  p_Dist->frame_ctr);
112     accumulate_avslice(metricSSIM,  p_Vid->type, p_Vid->p_Stats->frame_ctr[p_Vid->type]);
113   }
114 }
115 
116