1 /*
2  *  Copyright (c) 2010 The WebM project authors. All Rights Reserved.
3  *
4  *  Use of this source code is governed by a BSD-style license
5  *  that can be found in the LICENSE file in the root of the source
6  *  tree. An additional intellectual property rights grant can be found
7  *  in the file PATENTS.  All contributing project authors may
8  *  be found in the AUTHORS file in the root of the source tree.
9  */
10 
11 
12 #include "onyx_int.h"
13 
vp8_ssim_parms_16x16_c(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)14 void vp8_ssim_parms_16x16_c
15 (
16     unsigned char *s,
17     int sp,
18     unsigned char *r,
19     int rp,
20     unsigned long *sum_s,
21     unsigned long *sum_r,
22     unsigned long *sum_sq_s,
23     unsigned long *sum_sq_r,
24     unsigned long *sum_sxr
25 )
26 {
27     int i,j;
28     for(i=0;i<16;i++,s+=sp,r+=rp)
29      {
30          for(j=0;j<16;j++)
31          {
32              *sum_s += s[j];
33              *sum_r += r[j];
34              *sum_sq_s += s[j] * s[j];
35              *sum_sq_r += r[j] * r[j];
36              *sum_sxr += s[j] * r[j];
37          }
38      }
39 }
vp8_ssim_parms_8x8_c(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)40 void vp8_ssim_parms_8x8_c
41 (
42     unsigned char *s,
43     int sp,
44     unsigned char *r,
45     int rp,
46     unsigned long *sum_s,
47     unsigned long *sum_r,
48     unsigned long *sum_sq_s,
49     unsigned long *sum_sq_r,
50     unsigned long *sum_sxr
51 )
52 {
53     int i,j;
54     for(i=0;i<8;i++,s+=sp,r+=rp)
55      {
56          for(j=0;j<8;j++)
57          {
58              *sum_s += s[j];
59              *sum_r += r[j];
60              *sum_sq_s += s[j] * s[j];
61              *sum_sq_r += r[j] * r[j];
62              *sum_sxr += s[j] * r[j];
63          }
64      }
65 }
66 
67 const static int64_t cc1 =  26634; // (64^2*(.01*255)^2
68 const static int64_t cc2 = 239708; // (64^2*(.03*255)^2
69 
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)70 static double similarity
71 (
72     unsigned long sum_s,
73     unsigned long sum_r,
74     unsigned long sum_sq_s,
75     unsigned long sum_sq_r,
76     unsigned long sum_sxr,
77     int count
78 )
79 {
80     int64_t ssim_n, ssim_d;
81     int64_t c1, c2;
82 
83     //scale the constants by number of pixels
84     c1 = (cc1*count*count)>>12;
85     c2 = (cc2*count*count)>>12;
86 
87     ssim_n = (2*sum_s*sum_r+ c1)*((int64_t) 2*count*sum_sxr-
88           (int64_t) 2*sum_s*sum_r+c2);
89 
90     ssim_d = (sum_s*sum_s +sum_r*sum_r+c1)*
91         ((int64_t)count*sum_sq_s-(int64_t)sum_s*sum_s +
92         (int64_t)count*sum_sq_r-(int64_t) sum_r*sum_r +c2) ;
93 
94     return ssim_n * 1.0 / ssim_d;
95 }
96 
ssim_16x16(unsigned char * s,int sp,unsigned char * r,int rp)97 static double ssim_16x16(unsigned char *s,int sp, unsigned char *r,int rp)
98 {
99     unsigned long sum_s=0,sum_r=0,sum_sq_s=0,sum_sq_r=0,sum_sxr=0;
100     vp8_ssim_parms_16x16(s, sp, r, rp, &sum_s, &sum_r, &sum_sq_s, &sum_sq_r, &sum_sxr);
101     return similarity(sum_s, sum_r, sum_sq_s, sum_sq_r, sum_sxr, 256);
102 }
ssim_8x8(unsigned char * s,int sp,unsigned char * r,int rp)103 static double ssim_8x8(unsigned char *s,int sp, unsigned char *r,int rp)
104 {
105     unsigned long sum_s=0,sum_r=0,sum_sq_s=0,sum_sq_r=0,sum_sxr=0;
106     vp8_ssim_parms_8x8(s, sp, r, rp, &sum_s, &sum_r, &sum_sq_s, &sum_sq_r, &sum_sxr);
107     return similarity(sum_s, sum_r, sum_sq_s, sum_sq_r, sum_sxr, 64);
108 }
109 
110 // TODO: (jbb) tried to scale this function such that we may be able to use it
111 // for distortion metric in mode selection code ( provided we do a reconstruction)
dssim(unsigned char * s,int sp,unsigned char * r,int rp)112 long dssim(unsigned char *s,int sp, unsigned char *r,int rp)
113 {
114     unsigned long sum_s=0,sum_r=0,sum_sq_s=0,sum_sq_r=0,sum_sxr=0;
115     int64_t ssim3;
116     int64_t ssim_n1,ssim_n2;
117     int64_t ssim_d1,ssim_d2;
118     int64_t ssim_t1,ssim_t2;
119     int64_t c1, c2;
120 
121     // normalize by 256/64
122     c1 = cc1*16;
123     c2 = cc2*16;
124 
125     vp8_ssim_parms_16x16(s, sp, r, rp, &sum_s, &sum_r, &sum_sq_s, &sum_sq_r, &sum_sxr);
126     ssim_n1 = (2*sum_s*sum_r+ c1);
127 
128     ssim_n2 =((int64_t) 2*256*sum_sxr-(int64_t) 2*sum_s*sum_r+c2);
129 
130     ssim_d1 =((int64_t)sum_s*sum_s +(int64_t)sum_r*sum_r+c1);
131 
132     ssim_d2 = (256 * (int64_t) sum_sq_s-(int64_t) sum_s*sum_s +
133                     (int64_t) 256*sum_sq_r-(int64_t) sum_r*sum_r +c2) ;
134 
135     ssim_t1 = 256 - 256 * ssim_n1 / ssim_d1;
136     ssim_t2 = 256 - 256 * ssim_n2 / ssim_d2;
137 
138     ssim3 = 256 *ssim_t1 * ssim_t2;
139     if(ssim3 <0 )
140         ssim3=0;
141     return (long)( ssim3  );
142 }
143 
144 // We are using a 8x8 moving window with starting location of each 8x8 window
145 // on the 4x4 pixel grid. Such arrangement allows the windows to overlap
146 // block boundaries to penalize blocking artifacts.
vp8_ssim2(unsigned char * img1,unsigned char * img2,int stride_img1,int stride_img2,int width,int height)147 double vp8_ssim2
148 (
149     unsigned char *img1,
150     unsigned char *img2,
151     int stride_img1,
152     int stride_img2,
153     int width,
154     int height
155 )
156 {
157     int i,j;
158     int samples =0;
159     double ssim_total=0;
160 
161     // sample point start with each 4x4 location
162     for(i=0; i < height-8; i+=4, img1 += stride_img1*4, img2 += stride_img2*4)
163     {
164         for(j=0; j < width-8; j+=4 )
165         {
166             double v = ssim_8x8(img1+j, stride_img1, img2+j, stride_img2);
167             ssim_total += v;
168             samples++;
169         }
170     }
171     ssim_total /= samples;
172     return ssim_total;
173 }
vp8_calc_ssim(YV12_BUFFER_CONFIG * source,YV12_BUFFER_CONFIG * dest,int lumamask,double * weight)174 double vp8_calc_ssim
175 (
176     YV12_BUFFER_CONFIG *source,
177     YV12_BUFFER_CONFIG *dest,
178     int lumamask,
179     double *weight
180 )
181 {
182     double a, b, c;
183     double ssimv;
184 
185     a = vp8_ssim2(source->y_buffer, dest->y_buffer,
186                  source->y_stride, dest->y_stride, source->y_width,
187                  source->y_height);
188 
189     b = vp8_ssim2(source->u_buffer, dest->u_buffer,
190                  source->uv_stride, dest->uv_stride, source->uv_width,
191                  source->uv_height);
192 
193     c = vp8_ssim2(source->v_buffer, dest->v_buffer,
194                  source->uv_stride, dest->uv_stride, source->uv_width,
195                  source->uv_height);
196 
197     ssimv = a * .8 + .1 * (b + c);
198 
199     *weight = 1;
200 
201     return ssimv;
202 }
203 
vp8_calc_ssimg(YV12_BUFFER_CONFIG * source,YV12_BUFFER_CONFIG * dest,double * ssim_y,double * ssim_u,double * ssim_v)204 double vp8_calc_ssimg
205 (
206     YV12_BUFFER_CONFIG *source,
207     YV12_BUFFER_CONFIG *dest,
208     double *ssim_y,
209     double *ssim_u,
210     double *ssim_v
211 )
212 {
213     double ssim_all = 0;
214     double a, b, c;
215 
216     a = vp8_ssim2(source->y_buffer, dest->y_buffer,
217                  source->y_stride, dest->y_stride, source->y_width,
218                  source->y_height);
219 
220     b = vp8_ssim2(source->u_buffer, dest->u_buffer,
221                  source->uv_stride, dest->uv_stride, source->uv_width,
222                  source->uv_height);
223 
224     c = vp8_ssim2(source->v_buffer, dest->v_buffer,
225                  source->uv_stride, dest->uv_stride, source->uv_width,
226                  source->uv_height);
227     *ssim_y = a;
228     *ssim_u = b;
229     *ssim_v = c;
230     ssim_all = (a * 4 + b + c) /6;
231 
232     return ssim_all;
233 }
234