[cvs] / xvidcore / src / plugins / plugin_psnrhvsm.c Repository:
ViewVC logotype

Annotation of /xvidcore/src/plugins/plugin_psnrhvsm.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (view) (download)

1 : Isibaar 1.1 /*****************************************************************************
2 :     *
3 :     * XVID MPEG-4 VIDEO CODEC
4 :     * - PSNR-HVS-M plugin: computes the PSNR-HVS-M metric -
5 :     *
6 :     * Copyright(C) 2010 Michael Militzer <michael@xvid.org>
7 :     *
8 :     * This program is free software ; you can redistribute it and/or modify
9 :     * it under the terms of the GNU General Public License as published by
10 :     * the Free Software Foundation ; either version 2 of the License, or
11 :     * (at your option) any later version.
12 :     *
13 :     * This program is distributed in the hope that it will be useful,
14 :     * but WITHOUT ANY WARRANTY ; without even the implied warranty of
15 :     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 :     * GNU General Public License for more details.
17 :     *
18 :     * You should have received a copy of the GNU General Public License
19 :     * along with this program ; if not, write to the Free Software
20 :     * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 :     *
22 : Isibaar 1.2 * $Id: plugin_psnrhvsm.c,v 1.1 2010/10/10 19:19:46 Isibaar Exp $
23 : Isibaar 1.1 *
24 :     ****************************************************************************/
25 :    
26 :     /*****************************************************************************
27 :     *
28 :     * The PSNR-HVS-M metric is described in the following paper:
29 :     *
30 :     * "On between-coefficient contrast masking of DCT basis functions", by
31 :     * N. Ponomarenko, F. Silvestri, K. Egiazarian, M. Carli, J. Astola, V. Lukin,
32 :     * in Proceedings of the Third International Workshop on Video Processing and
33 :     * Quality Metrics for Consumer Electronics VPQM-07, January, 2007, 4 p.
34 :     *
35 :     * http://www.ponomarenko.info/psnrhvsm.htm
36 :     *
37 :     ****************************************************************************/
38 :    
39 :     #include <stdlib.h>
40 :     #include <stdio.h>
41 :     #include <math.h>
42 :     #include "../portab.h"
43 :     #include "../xvid.h"
44 :     #include "../dct/fdct.h"
45 :     #include "../image/image.h"
46 :     #include "../utils/mem_transfer.h"
47 :     #include "../utils/emms.h"
48 :    
49 :     typedef struct {
50 :    
51 :     uint64_t mse_sum; /* for avrg psnrhvsm */
52 :     long frame_cnt;
53 :    
54 :     } psnrhvsm_data_t; /* internal plugin data */
55 :    
56 :     static const float CSF_Coeff[64] =
57 :     { 1.608443f, 2.339554f, 2.573509f, 1.608443f, 1.072295f, 0.643377f, 0.504610f, 0.421887f,
58 :     2.144591f, 2.144591f, 1.838221f, 1.354478f, 0.989811f, 0.443708f, 0.428918f, 0.467911f,
59 :     1.838221f, 1.979622f, 1.608443f, 1.072295f, 0.643377f, 0.451493f, 0.372972f, 0.459555f,
60 :     1.838221f, 1.513829f, 1.169777f, 0.887417f, 0.504610f, 0.295806f, 0.321689f, 0.415082f,
61 :     1.429727f, 1.169777f, 0.695543f, 0.459555f, 0.378457f, 0.236102f, 0.249855f, 0.334222f,
62 :     1.072295f, 0.735288f, 0.467911f, 0.402111f, 0.317717f, 0.247453f, 0.227744f, 0.279729f,
63 :     0.525206f, 0.402111f, 0.329937f, 0.295806f, 0.249855f, 0.212687f, 0.214459f, 0.254803f,
64 :     0.357432f, 0.279729f, 0.270896f, 0.262603f, 0.229778f, 0.257351f, 0.249855f, 0.259950f
65 :     };
66 :    
67 :     static const float Mask_Coeff[64] =
68 :     { 0.000000f, 0.826446f, 1.000000f, 0.390625f, 0.173611f, 0.062500f, 0.038447f, 0.026874f,
69 :     0.694444f, 0.694444f, 0.510204f, 0.277008f, 0.147929f, 0.029727f, 0.027778f, 0.033058f,
70 :     0.510204f, 0.591716f, 0.390625f, 0.173611f, 0.062500f, 0.030779f, 0.021004f, 0.031888f,
71 :     0.510204f, 0.346021f, 0.206612f, 0.118906f, 0.038447f, 0.013212f, 0.015625f, 0.026015f,
72 :     0.308642f, 0.206612f, 0.073046f, 0.031888f, 0.021626f, 0.008417f, 0.009426f, 0.016866f,
73 :     0.173611f, 0.081633f, 0.033058f, 0.024414f, 0.015242f, 0.009246f, 0.007831f, 0.011815f,
74 :     0.041649f, 0.024414f, 0.016437f, 0.013212f, 0.009426f, 0.006830f, 0.006944f, 0.009803f,
75 :     0.019290f, 0.011815f, 0.011080f, 0.010412f, 0.007972f, 0.010000f, 0.009426f, 0.010203f
76 :     };
77 :    
78 : Isibaar 1.2 #if 0 /* Floating-point implementation */
79 : Isibaar 1.1
80 :     static uint32_t Calc_MSE_H(int16_t *DCT_A, int16_t *DCT_B, uint8_t *IMG_A, uint8_t *IMG_B, int stride)
81 :     {
82 :     int x, y, i, j;
83 :     uint32_t Global_A, Global_B, Sum_A = 0, Sum_B = 0;
84 :     uint32_t Local[8] = {0, 0, 0, 0, 0, 0, 0, 0};
85 :     uint32_t Local_Square[8] = {0, 0, 0, 0, 0, 0, 0, 0};
86 :     float MASK_A = 0.f, MASK_B = 0.f;
87 : Isibaar 1.2 float Mult1 = 1.f, Mult2 = 1.f;
88 : Isibaar 1.1 uint32_t MSE_H = 0;
89 :    
90 :     /* Step 1: Calculate CSF weighted energy of DCT coefficients */
91 :     for (y = 0; y < 8; y++) {
92 :     for (x = 0; x < 8; x++) {
93 :     MASK_A += (float)(DCT_A[y*8 + x]*DCT_A[y*8 + x])*Mask_Coeff[y*8 + x];
94 :     MASK_B += (float)(DCT_B[y*8 + x]*DCT_B[y*8 + x])*Mask_Coeff[y*8 + x];
95 :     }
96 :     }
97 :    
98 :     /* Step 2: Determine local variances compared to entire block variance */
99 :     for (y = 0; y < 2; y++) {
100 :     for (x = 0; x < 2; x++) {
101 : Isibaar 1.2 for (j = 0; j < 4; j++) {
102 :     for (i = 0; i < 4; i++) {
103 :     uint8_t A = IMG_A[(y*4+j)*stride + 4*x + i];
104 :     uint8_t B = IMG_B[(y*4+j)*stride + 4*x + i];
105 :    
106 :     Local[y*2 + x] += A;
107 :     Local[y*2 + x + 4] += B;
108 :     Local_Square[y*2 + x] += A*A;
109 :     Local_Square[y*2 + x + 4] += B*B;
110 :     }
111 : Isibaar 1.1 }
112 :     }
113 :     }
114 :    
115 :     Global_A = Local[0] + Local[1] + Local[2] + Local[3];
116 :     Global_B = Local[4] + Local[5] + Local[6] + Local[7];
117 :    
118 :     for (i = 0; i < 8; i++)
119 :     Local[i] = (Local_Square[i]<<4) - (Local[i]*Local[i]); /* 16*Var(Di) */
120 :    
121 :     Local_Square[0] += (Local_Square[1] + Local_Square[2] + Local_Square[3]);
122 :     Local_Square[4] += (Local_Square[5] + Local_Square[6] + Local_Square[7]);
123 :    
124 :     Global_A = (Local_Square[0]<<6) - Global_A*Global_A; /* 64*Var(D) */
125 :     Global_B = (Local_Square[4]<<6) - Global_B*Global_B; /* 64*Var(D) */
126 :    
127 :     /* Step 3: Calculate contrast masking threshold */
128 : Isibaar 1.2 if (Global_A)
129 :     Mult1 = (float)(Local[0]+Local[1]+Local[2]+Local[3])/((float)(Global_A)/4.f);
130 :    
131 :     if (Global_B)
132 :     Mult2 = (float)(Local[4]+Local[5]+Local[6]+Local[7])/((float)(Global_B)/4.f);
133 :    
134 :     MASK_A = (float)sqrt(MASK_A * Mult1) / 32.f;
135 :     MASK_B = (float)sqrt(MASK_B * Mult2) / 32.f;
136 : Isibaar 1.1
137 :     if (MASK_B > MASK_A) MASK_A = MASK_B; /* MAX(MASK_A, MASK_B) */
138 :    
139 :     /* Step 4: Calculate MSE of DCT coeffs reduced by masking effect */
140 :     for (j = 0; j < 8; j++) {
141 :     for (i = 0; i < 8; i++) {
142 :     float u = (float)abs(DCT_A[j*8 + i] - DCT_B[j*8 + i]);
143 :    
144 :     if ((i|j)>0) {
145 :     if (u < (MASK_A / Mask_Coeff[j*8 + i]))
146 :     u = 0; /* The error is not perceivable */
147 :     else
148 :     u -= (MASK_A / Mask_Coeff[j*8 + i]);
149 :     }
150 :    
151 :     MSE_H += (uint32_t) ((256.f*(u * CSF_Coeff[j*8 + i])*(u * CSF_Coeff[j*8 + i])) + 0.5f);
152 :     }
153 :     }
154 :     return MSE_H; /* Fixed-point value right-shifted by eight */
155 :     }
156 :    
157 :     #else /* First draft of a fixed-point implementation.
158 :     Might serve as a template for MMX/SSE code */
159 :    
160 :     static const uint16_t iMask_Coeff[64] =
161 :     { 0, 59577, 65535, 40959, 27306, 16384, 12850, 10743,
162 :     54612, 54612, 46811, 34492, 25206, 11299, 10923, 11915,
163 :     46811, 50412, 40959, 27306, 16384, 11497, 9498, 11703,
164 :     46811, 38550, 29789, 22598, 12850, 7533, 8192, 10570,
165 :     36408, 29789, 17712, 11703, 9637, 6012, 6363, 8511,
166 :     27306, 18724, 11915, 10240, 8091, 6302, 5799, 7123,
167 :     13374, 10240, 8402, 7533, 6363, 5416, 5461, 6489,
168 :     9102, 7123, 6898, 6687, 5851, 6553, 6363, 6620
169 :     };
170 :    
171 :     static const uint16_t Inv_iMask_Coeff[64] =
172 :     { 0, 310, 256, 655, 1475, 4096, 6659, 9526,
173 :     369, 369, 502, 924, 1731, 8612, 9216, 7744,
174 :     502, 433, 655, 1475, 4096, 8317, 12188, 8028,
175 :     502, 740, 1239, 2153, 6659, 19376, 16384, 9840,
176 :     829, 1239, 3505, 8028, 11838, 30415, 27159, 15178,
177 :     1475, 3136, 7744, 10486, 16796, 27688, 32691, 21667,
178 :     6147, 10486, 15575, 19376, 27159, 37482, 36866, 26114,
179 :     13271, 21667, 23105, 24587, 32112, 25600, 27159, 25091
180 :     };
181 :    
182 :     static const uint16_t iCSF_Coeff[64] =
183 :     { 1647, 2396, 2635, 1647, 1098, 659, 517, 432,
184 :     2196, 2196, 1882, 1387, 1014, 454, 439, 479,
185 :     1882, 2027, 1647, 1098, 659, 462, 382, 471,
186 :     1882, 1550, 1198, 909, 517, 303, 329, 425,
187 :     1464, 1198, 712, 471, 388, 242, 256, 342,
188 :     1098, 753, 479, 412, 325, 253, 233, 286,
189 :     538, 412, 338, 303, 256, 218, 220, 261,
190 :     366, 286, 277, 269, 235, 264, 256, 266
191 :     };
192 :    
193 : Isibaar 1.2 static __inline uint32_t isqrt(unsigned long n)
194 :     {
195 :     uint32_t c = 0x8000;
196 :     uint32_t g = 0x8000;
197 :    
198 :     for(;;) {
199 :     if(g*g > n)
200 :     g ^= c;
201 :     c >>= 1;
202 :     if(c == 0)
203 :     return g;
204 :     g |= c;
205 :     }
206 :     }
207 :    
208 : Isibaar 1.1 static uint32_t Calc_MSE_H(int16_t *DCT_A, int16_t *DCT_B, uint8_t *IMG_A, uint8_t *IMG_B, int stride)
209 :     {
210 :     int x, y, i, j;
211 :     uint32_t Global_A, Global_B, Sum_A = 0, Sum_B = 0;
212 :     uint32_t Local[8] = {0, 0, 0, 0, 0, 0, 0, 0};
213 :     uint32_t Local_Square[8] = {0, 0, 0, 0, 0, 0, 0, 0};
214 :     uint32_t MASK;
215 :     uint32_t MSE_H = 0;
216 :    
217 :     /* Step 1: Calculate CSF weighted energy of DCT coefficients */
218 :     for (y = 0; y < 8; y++) {
219 :     for (x = 0; x < 8; x++) {
220 :     uint16_t A = (abs(DCT_A[y*8 + x]) * iMask_Coeff[y*8 + x] + 4096) >> 13;
221 :     uint16_t B = (abs(DCT_B[y*8 + x]) * iMask_Coeff[y*8 + x] + 4096) >> 13;
222 :    
223 :     Sum_A += ((A*A) >> 3); /* PMADDWD */
224 :     Sum_B += ((B*B) >> 3);
225 :     }
226 :     }
227 :    
228 :     /* Step 2: Determine local variances compared to entire block variance */
229 :     for (y = 0; y < 2; y++) {
230 :     for (x = 0; x < 2; x++) {
231 : Isibaar 1.2 for (j = 0; j < 4; j++) {
232 :     for (i = 0; i < 4; i++) {
233 :     uint8_t A = IMG_A[(y*4+j)*stride + 4*x + i];
234 :     uint8_t B = IMG_B[(y*4+j)*stride + 4*x + i];
235 :    
236 :     Local[y*2 + x] += A;
237 :     Local[y*2 + x + 4] += B;
238 :     Local_Square[y*2 + x] += A*A;
239 :     Local_Square[y*2 + x + 4] += B*B;
240 :     }
241 : Isibaar 1.1 }
242 :     }
243 :     }
244 :    
245 :     Global_A = Local[0] + Local[1] + Local[2] + Local[3];
246 :     Global_B = Local[4] + Local[5] + Local[6] + Local[7];
247 :    
248 :     for (i = 0; i < 8; i++)
249 :     Local[i] = (Local_Square[i]<<4) - (Local[i]*Local[i]); /* 16*Var(Di) */
250 :    
251 :     Local_Square[0] += (Local_Square[1] + Local_Square[2] + Local_Square[3]);
252 :     Local_Square[4] += (Local_Square[5] + Local_Square[6] + Local_Square[7]);
253 :    
254 :     Global_A = (Local_Square[0]<<6) - Global_A*Global_A; /* 64*Var(D) */
255 :     Global_B = (Local_Square[4]<<6) - Global_B*Global_B; /* 64*Var(D) */
256 :    
257 :     /* Step 3: Calculate contrast masking threshold */
258 :     {
259 : Isibaar 1.2 uint32_t MASK_A, MASK_B;
260 :     uint32_t Mult1 = 64, Mult2 = 64;
261 :    
262 :     if (Global_A)
263 :     Mult1 = ((Local[0]+Local[1]+Local[2]+Local[3])<<8) / Global_A;
264 : Isibaar 1.1
265 : Isibaar 1.2 if (Global_B)
266 :     Mult2 = ((Local[4]+Local[5]+Local[6]+Local[7])<<8) / Global_B;
267 : Isibaar 1.1
268 : Isibaar 1.2 MASK_A = isqrt(2*Sum_A*Mult1) + 16;
269 :     MASK_B = isqrt(2*Sum_B*Mult2) + 16;
270 : Isibaar 1.1
271 :     if (MASK_B > MASK_A) /* MAX(MASK_A, MASK_B) */
272 : Isibaar 1.2 MASK = (uint32_t) MASK_B;
273 : Isibaar 1.1 else
274 : Isibaar 1.2 MASK = (uint32_t) MASK_A;
275 : Isibaar 1.1 }
276 :    
277 :     /* Step 4: Calculate MSE of DCT coeffs reduced by masking effect */
278 :     for (j = 0; j < 8; j++) {
279 :     for (i = 0; i < 8; i++) {
280 :     uint32_t Thresh = (MASK * Inv_iMask_Coeff[j*8 + i] + 128) >> 8;
281 :     uint32_t u = abs(DCT_A[j*8 + i] - DCT_B[j*8 + i]) << 10;
282 :    
283 :     if ((i|j)>0) {
284 :     if (u < Thresh)
285 :     u = 0; /* The error is not perceivable */
286 :     else
287 :     u -= Thresh;
288 :     }
289 :    
290 :     {
291 :     uint64_t tmp = (u * iCSF_Coeff[j*8 + i] + 512) >> 10;
292 :     MSE_H += (uint32_t) ((tmp * tmp) >> 12);
293 :     }
294 :     }
295 :     }
296 :     return MSE_H;
297 :     }
298 :    
299 :     #endif
300 :    
301 :     static void psnrhvsm_after(xvid_plg_data_t *data, psnrhvsm_data_t *psnrhvsm)
302 :     {
303 :     DECLARE_ALIGNED_MATRIX(DCT, 2, 64, int16_t, CACHE_LINE);
304 : Isibaar 1.2 uint32_t x, y, stride = data->original.stride[0];
305 : Isibaar 1.1 int16_t *DCT_A = &DCT[0], *DCT_B = &DCT[64];
306 :     uint8_t *IMG_A = (uint8_t *) data->original.plane[0];
307 :     uint8_t *IMG_B = (uint8_t *) data->current.plane[0];
308 : Isibaar 1.2 uint64_t MSE_H = 0;
309 : Isibaar 1.1
310 :     for (y = 0; y < data->height; y += 8) {
311 :     for (x = 0; x < data->width; x += 8) {
312 :     int offset = y*stride + x;
313 :    
314 :     emms();
315 :    
316 :     /* Transfer data */
317 :     transfer_8to16copy(DCT_A, IMG_A + offset, stride);
318 :     transfer_8to16copy(DCT_B, IMG_B + offset, stride);
319 :    
320 :     /* Perform DCT */
321 :     fdct(DCT_A);
322 :     fdct(DCT_B);
323 :    
324 :     emms();
325 :    
326 :     /* Calculate MSE_H reduced by contrast masking effect */
327 :     MSE_H += Calc_MSE_H(DCT_A, DCT_B, IMG_A + offset, IMG_B + offset, stride);
328 :     }
329 :     }
330 :    
331 : Isibaar 1.2 x = 4*MSE_H / (data->width * data->height);
332 :     psnrhvsm->mse_sum += x;
333 : Isibaar 1.1 psnrhvsm->frame_cnt++;
334 :    
335 : Isibaar 1.2 printf(" psnrhvsm: %2.2f\n", sse_to_PSNR(x, 1024));
336 : Isibaar 1.1 }
337 :    
338 :     static int psnrhvsm_create(xvid_plg_create_t *create, void **handle)
339 :     {
340 :     psnrhvsm_data_t *psnrhvsm;
341 :     psnrhvsm = (psnrhvsm_data_t *) malloc(sizeof(psnrhvsm_data_t));
342 :    
343 :     psnrhvsm->mse_sum = 0;
344 :     psnrhvsm->frame_cnt = 0;
345 :    
346 :     *(handle) = (void*) psnrhvsm;
347 :     return 0;
348 :     }
349 :    
350 :     int xvid_plugin_psnrhvsm(void *handle, int opt, void *param1, void *param2)
351 :     {
352 :     switch(opt) {
353 :     case(XVID_PLG_INFO):
354 :     ((xvid_plg_info_t *)param1)->flags = XVID_REQORIGINAL;
355 :     break;
356 :     case(XVID_PLG_CREATE):
357 :     psnrhvsm_create((xvid_plg_create_t *)param1,(void **)param2);
358 :     break;
359 :     case(XVID_PLG_BEFORE):
360 :     case(XVID_PLG_FRAME):
361 :     break;
362 :     case(XVID_PLG_AFTER):
363 :     psnrhvsm_after((xvid_plg_data_t *)param1, (psnrhvsm_data_t *)handle);
364 :     break;
365 :     case(XVID_PLG_DESTROY):
366 :     {
367 :     uint32_t MSE_H;
368 :     psnrhvsm_data_t *psnrhvsm = (psnrhvsm_data_t *)handle;
369 :    
370 :     if (psnrhvsm) {
371 :     MSE_H = (uint32_t) (psnrhvsm->mse_sum / psnrhvsm->frame_cnt);
372 :    
373 :     emms();
374 : Isibaar 1.2 printf("Average psnrhvsm: %2.2f\n", sse_to_PSNR(MSE_H, 1024));
375 : Isibaar 1.1 free(psnrhvsm);
376 :     }
377 :     }
378 :     break;
379 :     default:
380 :     break;
381 :     }
382 :     return 0;
383 :     };

No admin address has been configured
ViewVC Help
Powered by ViewVC 1.0.4