1
2 /*!
3 *************************************************************************************
4 * \file img_dist_ms_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 * - Peshala Pahalawatta <ppaha@dolby.com>
13 * - Zhen Li <zli@dolby.com>
14 * - Alexis Michael Tourapis <alexismt@ieee.org>
15 *************************************************************************************
16 */
17 #include "contributors.h"
18 #include "global.h"
19 #include "img_distortion.h"
20 #include "enc_statistics.h"
21 #include "memalloc.h"
22 #include "math.h"
23
24 #ifndef UNBIASED_VARIANCE
25 #define UNBIASED_VARIANCE // unbiased estimation of the variance
26 #endif
27
28 #define MS_SSIM_PAD 6
29 #define MS_SSIM_PAD2 3
30
31 #define MS_SSIM_BETA0 0.0448
32 #define MS_SSIM_BETA1 0.2856
33 #define MS_SSIM_BETA2 0.3001
34 #define MS_SSIM_BETA3 0.2363
35 #define MS_SSIM_BETA4 0.1333
36
37 #define MAX_SSIM_LEVELS 5
38
39 //Computes the product of the contrast and structure componenents of the structural similarity metric.
compute_structural_components(VideoParameters * p_Vid,InputParameters * p_Inp,imgpel ** refImg,imgpel ** encImg,int height,int width,int win_height,int win_width,int comp)40 float compute_structural_components (VideoParameters *p_Vid, InputParameters *p_Inp, imgpel **refImg, imgpel **encImg, int height, int width, int win_height, int win_width, int comp)
41 {
42 static const float K2 = 0.03f;
43 float max_pix_value_sqd;
44 float C2;
45 float win_pixels = (float) (win_width * win_height);
46 #ifdef UNBIASED_VARIANCE
47 float win_pixels_bias = win_pixels - 1;
48 #else
49 float win_pixels_bias = win_pixels;
50 #endif
51 float mb_ssim, meanOrg, meanEnc;
52 float varOrg, varEnc, covOrgEnc;
53 int imeanOrg, imeanEnc, ivarOrg, ivarEnc, icovOrgEnc;
54 float cur_distortion = 0.0;
55 int i, j, n, m, win_cnt = 0;
56 int overlapSize = p_Inp->SSIMOverlapSize;
57
58 max_pix_value_sqd = (float) (p_Vid->max_pel_value_comp[comp] * p_Vid->max_pel_value_comp[comp]);
59 C2 = K2 * K2 * max_pix_value_sqd;
60
61 for (j = 0; j <= height - win_height; j += overlapSize)
62 {
63 for (i = 0; i <= width - win_width; i += overlapSize)
64 {
65 imeanOrg = 0;
66 imeanEnc = 0;
67 ivarOrg = 0;
68 ivarEnc = 0;
69 icovOrgEnc = 0;
70
71 for ( n = j;n < j + win_height;n ++)
72 {
73 for (m = i;m < i + win_width;m ++)
74 {
75 imeanOrg += refImg[n][m];
76 imeanEnc += encImg[n][m];
77 ivarOrg += refImg[n][m] * refImg[n][m];
78 ivarEnc += encImg[n][m] * encImg[n][m];
79 icovOrgEnc += refImg[n][m] * encImg[n][m];
80 }
81 }
82
83 meanOrg = (float) imeanOrg / win_pixels;
84 meanEnc = (float) imeanEnc / win_pixels;
85
86 varOrg = ((float) ivarOrg - ((float) imeanOrg) * meanOrg) / win_pixels_bias;
87 varEnc = ((float) ivarEnc - ((float) imeanEnc) * meanEnc) / win_pixels_bias;
88 covOrgEnc = ((float) icovOrgEnc - ((float) imeanOrg) * meanEnc) / win_pixels_bias;
89
90 mb_ssim = (float) (2.0 * covOrgEnc + C2);
91 mb_ssim /= (float) (varOrg + varEnc + C2);
92
93 cur_distortion += mb_ssim;
94 win_cnt++;
95 }
96 }
97
98 cur_distortion /= (float) win_cnt;
99
100 if (cur_distortion >= 1.0 && cur_distortion < 1.01) // avoid float accuracy problem at very low QP(e.g.2)
101 cur_distortion = 1.0;
102
103 return cur_distortion;
104 }
105
compute_luminance_component(VideoParameters * p_Vid,InputParameters * p_Inp,imgpel ** refImg,imgpel ** encImg,int height,int width,int win_height,int win_width,int comp)106 float compute_luminance_component (VideoParameters *p_Vid, InputParameters *p_Inp, imgpel **refImg, imgpel **encImg, int height, int width, int win_height, int win_width, int comp)
107 {
108 static const float K1 = 0.01f;
109 float max_pix_value_sqd;
110 float C1;
111 float win_pixels = (float) (win_width * win_height);
112 float mb_ssim, meanOrg, meanEnc;
113 int imeanOrg, imeanEnc;
114 float cur_distortion = 0.0;
115 int i, j, n, m, win_cnt = 0;
116 int overlapSize = p_Inp->SSIMOverlapSize;
117
118 max_pix_value_sqd = (float) (p_Vid->max_pel_value_comp[comp] * p_Vid->max_pel_value_comp[comp]);
119 C1 = K1 * K1 * max_pix_value_sqd;
120
121 for (j = 0; j <= height - win_height; j += overlapSize)
122 {
123 for (i = 0; i <= width - win_width; i += overlapSize)
124 {
125 imeanOrg = 0;
126 imeanEnc = 0;
127
128 for ( n = j;n < j + win_height;n ++)
129 {
130 for (m = i;m < i + win_width;m ++)
131 {
132 imeanOrg += refImg[n][m];
133 imeanEnc += encImg[n][m];
134 }
135 }
136
137 meanOrg = (float) imeanOrg / win_pixels;
138 meanEnc = (float) imeanEnc / win_pixels;
139
140 mb_ssim = (float) (2.0 * meanOrg * meanEnc + C1);
141 mb_ssim /= (float) (meanOrg * meanOrg + meanEnc * meanEnc + C1);
142
143 cur_distortion += mb_ssim;
144 win_cnt++;
145 }
146 }
147
148 cur_distortion /= (float) win_cnt;
149
150 if (cur_distortion >= 1.0 && cur_distortion < 1.01) // avoid float accuracy problem at very low QP(e.g.2)
151 cur_distortion = 1.0;
152
153 return cur_distortion;
154 }
155
horizontal_symmetric_extension(int ** buffer,int width,int height)156 void horizontal_symmetric_extension(int **buffer, int width, int height )
157 {
158 int j;
159 int* buf;
160
161 int height_plus_pad2 = height + MS_SSIM_PAD2;
162 int width_plus_pad2_minus_one = width + MS_SSIM_PAD2 - 1;
163
164 // horizontal PSE
165 for (j = MS_SSIM_PAD2; j < height_plus_pad2; j++ ) {
166 // left column
167 buf = &buffer[j][MS_SSIM_PAD2];
168 buf[-1] = buf[1];
169 buf[-2] = buf[2];
170
171 // right column
172 buf = &buffer[j][width_plus_pad2_minus_one];
173 buf[1] = buf[-1];
174 buf[2] = buf[-2];
175 buf[3] = buf[-3];
176 }
177 }
178
vertical_symmetric_extension(int ** buffer,int width,int height)179 void vertical_symmetric_extension(int **buffer, int width, int height)
180 {
181 int i;
182
183 int *bufminus1 = &buffer[MS_SSIM_PAD2-1][MS_SSIM_PAD2];
184 int *bufminus2 = &buffer[MS_SSIM_PAD2-2][MS_SSIM_PAD2];
185 int *bufplus1 = &buffer[MS_SSIM_PAD2+1][MS_SSIM_PAD2];
186 int *bufplus2 = &buffer[MS_SSIM_PAD2+2][MS_SSIM_PAD2];
187
188 int *bufhminus1 = &buffer[height + MS_SSIM_PAD2-2][MS_SSIM_PAD2];
189 int *bufhminus2 = &buffer[height + MS_SSIM_PAD2-3][MS_SSIM_PAD2];
190 int *bufhminus3 = &buffer[height + MS_SSIM_PAD2-4][MS_SSIM_PAD2];
191 int *bufhplus1 = &buffer[height + MS_SSIM_PAD2][MS_SSIM_PAD2];
192 int *bufhplus2 = &buffer[height + MS_SSIM_PAD2+1][MS_SSIM_PAD2];
193 int *bufhplus3 = &buffer[height + MS_SSIM_PAD2+2][MS_SSIM_PAD2];
194
195 for (i = 0; i < width; i++)
196 {
197 //Top Rows
198 bufminus1[i] = bufplus1[i];
199 bufminus2[i] = bufplus2[i];
200 //Bottom Rows
201 bufhplus1[i] = bufhminus1[i];
202 bufhplus2[i] = bufhminus2[i];
203 bufhplus3[i] = bufhminus3[i];
204 }
205 }
206
imgpel_to_padded_int(imgpel ** src,int ** buffer,int width,int height)207 static void imgpel_to_padded_int(imgpel** src, int **buffer, int width, int height)
208 {
209 int i, j;
210 int* tmpDst;
211 imgpel* tmpSrc;
212
213 tmpDst = buffer[MS_SSIM_PAD2];
214 for (j = 0; j < height; j++)
215 {
216 tmpSrc = src[j];
217 tmpDst = &buffer[j + MS_SSIM_PAD2][MS_SSIM_PAD2];
218 for (i = 0; i < width; i++)
219 {
220 tmpDst[i] = (int)tmpSrc[i];
221 }
222 }
223 }
224
downsample(imgpel ** src,imgpel ** out,int height,int width)225 void downsample(imgpel** src, imgpel** out, int height, int width)
226 {
227 int height2 = height >> 1;
228 int width2 = width >> 1;
229 int i, j;
230 int ii, jj;
231 int iDst;
232 int tmp, tmp1, tmp2;
233 int* tmpDst;
234 int* tmpSrc;
235
236 int** itemp;
237 int** dest;
238
239 get_mem2Dint(&itemp, height + MS_SSIM_PAD, width + MS_SSIM_PAD);
240 get_mem2Dint(&dest, height + MS_SSIM_PAD, width2 + MS_SSIM_PAD);
241
242 imgpel_to_padded_int(src, itemp, width, height);
243 horizontal_symmetric_extension(itemp, width, height);
244
245 for (j = 0; j < height; j++)
246 {
247 tmpDst = dest[j+MS_SSIM_PAD2];
248 tmpSrc = itemp[j+MS_SSIM_PAD2];
249 iDst = MS_SSIM_PAD2;
250 for (i = 0; i < width2; i++, iDst++)
251 {
252 ii = (i << 1) + MS_SSIM_PAD2;
253 tmp1 = tmpSrc[ii-1] + tmpSrc[ii+2];
254 tmp2 = tmpSrc[ii] + tmpSrc[ii+1];
255 tmpDst[iDst] = tmpSrc[ii-2] + tmpSrc[ii+3] + (tmp1 << 1) + tmp1 + (tmp2 << 5) - (tmp2 << 2);
256 tmpDst[iDst] >>= 6;
257 }
258 }
259
260 //Periodic extension
261 vertical_symmetric_extension(dest, width2, height);
262
263 for (i = 0; i < width2; i++)
264 {
265 ii = i + MS_SSIM_PAD2;
266 for (j = 0; j < height2; j++)
267 {
268 jj = (j << 1) + MS_SSIM_PAD2;
269 tmp1 = dest[jj-1][ii] + dest[jj+2][ii];
270 tmp2 = dest[jj][ii] + dest[jj+1][ii];
271 tmp = dest[jj-2][ii] + dest[jj+3][ii] + (tmp1 << 1) + tmp1 + (tmp2 << 5) - (tmp2 << 2);
272 out[j][i] = (unsigned char) (tmp >> 6); //Note: Should change for different bit depths
273 }
274 }
275 free_mem2Dint(itemp);
276 free_mem2Dint(dest);
277 }
278
compute_ms_ssim(VideoParameters * p_Vid,InputParameters * p_Inp,imgpel ** refImg,imgpel ** encImg,int height,int width,int win_height,int win_width,int comp)279 float compute_ms_ssim(VideoParameters *p_Vid, InputParameters *p_Inp, imgpel **refImg, imgpel **encImg, int height, int width, int win_height, int win_width, int comp)
280 {
281 float structural[MAX_SSIM_LEVELS];
282 float cur_distortion;
283 float luminance;
284 imgpel** dsRef;
285 imgpel** dsEnc;
286 int m;
287 static const int max_ssim_levels_minus_one = MAX_SSIM_LEVELS - 1;
288 static const float exponent[5] = {(float)MS_SSIM_BETA0, (float)MS_SSIM_BETA1, (float)MS_SSIM_BETA2, (float)MS_SSIM_BETA3, (float)MS_SSIM_BETA4};
289
290 dsRef = NULL;
291 dsEnc = NULL;
292 get_mem2Dpel(&dsRef, height>>1, width>>1);
293 get_mem2Dpel(&dsEnc, height>>1, width>>1);
294
295 structural[0] = compute_structural_components(p_Vid, p_Inp, refImg, encImg, height, width, win_height, win_width, comp);
296 cur_distortion = (float)pow(structural[0], exponent[0]);
297
298 downsample(refImg, dsRef, height, width);
299 downsample(encImg, dsEnc, height, width);
300
301 for (m = 1; m < MAX_SSIM_LEVELS; m++)
302 {
303 height >>= 1;
304 width >>= 1;
305 structural[m] = compute_structural_components(p_Vid, p_Inp, dsRef, dsEnc, height, width, imin(win_height,height), imin(win_width,width), comp);
306 cur_distortion *= (float)pow(structural[m], exponent[m]);
307 if (m < max_ssim_levels_minus_one)
308 {
309 downsample(dsRef, dsRef, height, width);
310 downsample(dsEnc, dsEnc, height, width);
311 }
312 else
313 {
314 luminance = compute_luminance_component(p_Vid, p_Inp, dsRef, dsEnc, height, width, imin(win_height,height), imin(win_width,width), comp);
315 cur_distortion *= (float)pow(luminance, exponent[m]);
316 }
317 }
318 free_mem2Dpel(dsRef);
319 free_mem2Dpel(dsEnc);
320 dsRef = NULL;
321 dsEnc = NULL;
322
323 return cur_distortion;
324 }
325
326 /*!
327 ************************************************************************
328 * \brief
329 * Find MS-SSIM for all three components
330 ************************************************************************
331 */
find_ms_ssim(VideoParameters * p_Vid,InputParameters * p_Inp,ImageStructure * ref,ImageStructure * src,DistMetric * metricSSIM)332 void find_ms_ssim (VideoParameters *p_Vid, InputParameters *p_Inp, ImageStructure *ref, ImageStructure *src, DistMetric *metricSSIM)
333 {
334 DistortionParams *p_Dist = p_Vid->p_Dist;
335 FrameFormat *format = &ref->format;
336
337 metricSSIM->value[0] = compute_ms_ssim (p_Vid, p_Inp, ref->data[0], src->data[0], format->height[0], format->width[0], BLOCK_SIZE_8x8, BLOCK_SIZE_8x8, 0);
338 // Chroma.
339 if (format->yuv_format != YUV400)
340 {
341 metricSSIM->value[1] = compute_ms_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);
342 metricSSIM->value[2] = compute_ms_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);
343 }
344
345 {
346 accumulate_average(metricSSIM, p_Dist->frame_ctr);
347 accumulate_avslice(metricSSIM, p_Vid->type, p_Vid->p_Stats->frame_ctr[p_Vid->type]);
348 }
349 }
350