← Back to branch summary

~ubuntu-branches/ubuntu/wily/libde265/wily

« back to all changes in this revision

Viewing changes to tools/bjoentegaard.cc

  • Committer: Package Import Robot
  • Author(s): Joachim Bauch
  • Date: 2015-07-16 11:07:46 UTC
  • mfrom: (2.1.2 sid)
  • Revision ID: package-import@ubuntu.com-20150716110746-76vsv24j3yux7tnu
Tags: 1.0.2-1
* Imported Upstream version 1.0.2
* Added new files to copyright information.
* Only export decoder API and update symbols for new version.

Show diffs side-by-side

added added

removed removed

Lines of Context:
tools/bjoentegaard.cc
 
1
/*
 
2
 * H.265 video codec.
 
3
 * Copyright (c) 2013-2014 struktur AG, Dirk Farin <farin@struktur.de>
 
4
 *
 
5
 * This file is part of libde265.
 
6
 *
 
7
 * libde265 is free software: you can redistribute it and/or modify
 
8
 * it under the terms of the GNU General Public License as published by
 
9
 * the Free Software Foundation, either version 3 of the License, or
 
10
 * (at your option) any later version.
 
11
 *
 
12
 * libde265 is distributed in the hope that it will be useful,
 
13
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
14
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
15
 * GNU General Public License for more details.
 
16
 *
 
17
 * You should have received a copy of the GNU General Public License
 
18
 * along with libde265.  If not, see <http://www.gnu.org/licenses/>.
 
19
 */
 
20
 
 
21
#include <stdio.h>
 
22
#include <stdlib.h>
 
23
#include <fstream>
 
24
#include <sstream>
 
25
#include <string>
 
26
#include <vector>
 
27
#include <math.h>
 
28
#include <unistd.h>
 
29
 
 
30
 
 
31
const bool D = false;
 
32
 
 
33
 
 
34
/* There are numerical stability problems in the matrix inverse.
 
35
   Switching to long double seems to provide enough accuracy.
 
36
   TODO: in the long term, use a better regression algorithm.
 
37
 */
 
38
typedef long double FP;
 
39
 
 
40
 
 
41
struct datapoint
 
42
{
 
43
  double rate;
 
44
  double distortion;
 
45
};
 
46
 
 
47
struct BjoentegaardParams
 
48
{
 
49
  // a*log^3 R + b*log^2 R + c*log R + d
 
50
  double a,b,c,d;
 
51
 
 
52
  double minRate, maxRate;
 
53
};
 
54
 
 
55
std::vector<datapoint> curveA,curveB;
 
56
BjoentegaardParams paramsA,paramsB;
 
57
 
 
58
#define RATE_NORMALIZATION_FACTOR 1 //(1/1000.0)
 
59
 
 
60
 
 
61
 
 
62
FP invf(int i,int j,const FP* m)
 
63
{
 
64
  int o = 2+(j-i);
 
65
 
 
66
  i += 4+o;
 
67
  j += 4-o;
 
68
 
 
69
#define e(a,b) m[ ((j+b)%4)*4 + ((i+a)%4) ]
 
70
 
 
71
    FP inv =
 
72
      + e(+1,-1)*e(+0,+0)*e(-1,+1)
 
73
      + e(+1,+1)*e(+0,-1)*e(-1,+0)
 
74
      + e(-1,-1)*e(+1,+0)*e(+0,+1)
 
75
      - e(-1,-1)*e(+0,+0)*e(+1,+1)
 
76
      - e(-1,+1)*e(+0,-1)*e(+1,+0)
 
77
      - e(+1,-1)*e(-1,+0)*e(+0,+1);
 
78
 
 
79
    return (o%2)?inv : -inv;
 
80
 
 
81
    #undef e
 
82
 
 
83
}
 
84
 
 
85
bool inverseMatrix4x4(const FP *m, FP *out)
 
86
{
 
87
  FP inv[16];
 
88
 
 
89
  for(int i=0;i<4;i++)
 
90
    for(int j=0;j<4;j++)
 
91
      inv[j*4+i] = invf(i,j,m);
 
92
 
 
93
  FP D = 0;
 
94
 
 
95
  for(int k=0;k<4;k++) D += m[k] * inv[k*4];
 
96
 
 
97
  if (D == 0) return false;
 
98
 
 
99
  D = 1.0 / D;
 
100
 
 
101
  for (int i = 0; i < 16; i++)
 
102
    out[i] = inv[i] * D;
 
103
 
 
104
  return true;
 
105
 
 
106
}
 
107
 
 
108
 
 
109
 
 
110
BjoentegaardParams fitParams(const std::vector<datapoint>& curve)
 
111
{
 
112
  // build regression matrix
 
113
 
 
114
  int n = curve.size();
 
115
 
 
116
  FP X[4*n];  // regression matrix
 
117
  FP XT[n*4]; // transpose of X
 
118
 
 
119
  for (int i=0;i<n;i++) {
 
120
    FP x = log(curve[i].rate) * RATE_NORMALIZATION_FACTOR;
 
121
 
 
122
    X[4*i + 0] = 1;
 
123
    X[4*i + 1] = x;
 
124
    X[4*i + 2] = x*x;
 
125
    X[4*i + 3] = x*x*x;
 
126
 
 
127
    if (D) printf("%f %f %f %f ;\n",1.0,(double)x,(double)(x*x),(double)(x*x*x));
 
128
 
 
129
    XT[i+0*n] = 1;
 
130
    XT[i+1*n] = x;
 
131
    XT[i+2*n] = x*x;
 
132
    XT[i+3*n] = x*x*x;
 
133
  }
 
134
 
 
135
  if (D) {
 
136
    printf("rate: ");
 
137
    for (int i=0;i<n;i++) {
 
138
      printf("%f ; ",curve[i].rate);
 
139
    }
 
140
    printf("\n");
 
141
 
 
142
    printf("distortion: ");
 
143
    for (int i=0;i<n;i++) {
 
144
      printf("%f ; ",curve[i].distortion);
 
145
    }
 
146
    printf("\n");
 
147
  }
 
148
 
 
149
  // calc X^T * X
 
150
 
 
151
  FP XTX[4*4];
 
152
  for (int y=0;y<4;y++)
 
153
    for (int x=0;x<4;x++) {
 
154
      FP sum=0;
 
155
 
 
156
      for (int i=0;i<n;i++)
 
157
        {
 
158
          sum += XT[y*n + i] * X[x + i*4];
 
159
        }
 
160
 
 
161
      XTX[y*4+x] = sum;
 
162
    }
 
163
 
 
164
  FP XTXinv[4*4];
 
165
 
 
166
  inverseMatrix4x4(XTX, XTXinv);
 
167
 
 
168
  if (D) {
 
169
    for (int y=0;y<4;y++) {
 
170
      for (int x=0;x<4;x++) {
 
171
        printf("%f ",(double)XTXinv[y*4+x]);
 
172
      }
 
173
      printf("\n");
 
174
    }
 
175
  }
 
176
 
 
177
  // calculate pseudo-inverse XP = (X^T * X)^-1 * X^T
 
178
 
 
179
  FP XP[n*4];
 
180
 
 
181
  for (int y=0;y<4;y++) {
 
182
    for (int x=0;x<n;x++) {
 
183
      FP sum=0;
 
184
 
 
185
      for (int i=0;i<4;i++)
 
186
        {
 
187
          sum += XTXinv[y*4 + i] * XT[x + i*n];
 
188
        }
 
189
 
 
190
      XP[y*n+x] = sum;
 
191
    }
 
192
  }
 
193
 
 
194
  // calculate regression parameters
 
195
 
 
196
  FP p[4];
 
197
 
 
198
  for (int k=0;k<4;k++)
 
199
    {
 
200
      FP sum=0;
 
201
 
 
202
      for (int i=0;i<n;i++) {
 
203
        sum += XP[k*n + i] * curve[i].distortion;
 
204
      }
 
205
 
 
206
      p[k]=sum;
 
207
    }
 
208
 
 
209
 
 
210
  BjoentegaardParams param;
 
211
  param.d = p[0];
 
212
  param.c = p[1];
 
213
  param.b = p[2];
 
214
  param.a = p[3];
 
215
 
 
216
 
 
217
  // find min and max rate
 
218
 
 
219
  param.minRate = curve[0].rate;
 
220
  param.maxRate = curve[0].rate;
 
221
 
 
222
  for (int i=1;i<n;i++) {
 
223
    param.minRate = std::min(param.minRate, curve[i].rate);
 
224
    param.maxRate = std::max(param.maxRate, curve[i].rate);
 
225
  }
 
226
 
 
227
  return param;
 
228
}
 
229
 
 
230
FP evalIntegralAt(const BjoentegaardParams& p, double x)
 
231
{
 
232
  FP sum = 0;
 
233
 
 
234
  // integral of: d
 
235
 
 
236
  sum += p.d * x;
 
237
 
 
238
  // integral of: c*log(x)
 
239
 
 
240
  sum += p.c * x* (log(x)-1);
 
241
 
 
242
  // integral of: b*log(x)^2
 
243
 
 
244
  sum += p.b * x * ((log(x)-2)*log(x)+2);
 
245
 
 
246
  // integral of: a*log(x)^3
 
247
 
 
248
  sum += p.a * x * (log(x)*((log(x)-3)*log(x)+6)-6);
 
249
 
 
250
  return sum;
 
251
}
 
252
 
 
253
 
 
254
double calcBjoentegaard(const BjoentegaardParams& paramsA,
 
255
                        const BjoentegaardParams& paramsB,
 
256
                        double min_rate, double max_rate)
 
257
{
 
258
  double mini = std::max(paramsA.minRate, paramsB.minRate);
 
259
  double maxi = std::min(paramsA.maxRate, paramsB.maxRate);
 
260
 
 
261
  if (min_rate >= 0) mini = std::max(mini, min_rate);
 
262
  if (max_rate >= 0) maxi = std::min(maxi, max_rate);
 
263
 
 
264
  if (D) printf("range: %f %f\n",mini,maxi);
 
265
 
 
266
  FP intA = evalIntegralAt(paramsA, maxi) - evalIntegralAt(paramsA, mini);
 
267
  FP intB = evalIntegralAt(paramsB, maxi) - evalIntegralAt(paramsB, mini);
 
268
 
 
269
  if (D) printf("int1:%f int2:%f\n",(double)intA,(double)intB);
 
270
 
 
271
  return (intA-intB)/(maxi-mini);
 
272
}
 
273
 
 
274
 
 
275
std::vector<datapoint> readRDFile(const char* filename, float min_rate, float max_rate)
 
276
{
 
277
  std::vector<datapoint> curve;
 
278
  std::ifstream istr(filename);
 
279
 
 
280
  for (;;)
 
281
    {
 
282
      std::string line;
 
283
      getline(istr, line);
 
284
      if (istr.eof())
 
285
        break;
 
286
 
 
287
      if (line[0]=='#') continue;
 
288
 
 
289
      std::stringstream sstr(line);
 
290
      datapoint p;
 
291
      sstr >> p.rate >> p.distortion;
 
292
 
 
293
      if (min_rate>=0 && p.rate < min_rate) continue;
 
294
      if (max_rate>=0 && p.rate > max_rate) continue;
 
295
 
 
296
      curve.push_back(p);
 
297
    }
 
298
 
 
299
  return curve;
 
300
}
 
301
 
 
302
 
 
303
int main(int argc, char** argv)
 
304
{
 
305
  float min_rate = -1;
 
306
  float max_rate = -1;
 
307
 
 
308
  int c;
 
309
  while ((c=getopt(argc,argv, "l:h:")) != -1) {
 
310
    switch (c) {
 
311
    case 'l': min_rate = atof(optarg); break;
 
312
    case 'h': max_rate = atof(optarg); break;
 
313
    }
 
314
  }
 
315
 
 
316
  curveA = readRDFile(argv[optind], min_rate, max_rate);
 
317
  paramsA = fitParams(curveA);
 
318
 
 
319
  printf("params A: %f %f %f %f\n",paramsA.a,paramsA.b,paramsA.c,paramsA.d);
 
320
 
 
321
  printf("gnuplot: %f*log(x)**3+%f*log(x)**2+%f*log(x)+%f\n",paramsA.a,paramsA.b,paramsA.c,paramsA.d);
 
322
 
 
323
  if (optind+1<argc) {
 
324
    curveB = readRDFile(argv[optind+1], min_rate, max_rate);
 
325
    paramsB = fitParams(curveB);
 
326
 
 
327
    printf("params B: %f %f %f %f\n",paramsB.a,paramsB.b,paramsB.c,paramsB.d);
 
328
 
 
329
    printf("gnuplot: %f*log(x)**3+%f*log(x)**2+%f*log(x)+%f\n",paramsB.a,paramsB.b,paramsB.c,paramsB.d);
 
330
 
 
331
    double delta = calcBjoentegaard(paramsA,paramsB, min_rate,max_rate);
 
332
 
 
333
    printf("Bjoentegaard delta: %f dB   (A-B -> >0 -> first (A) is better)\n",delta);
 
334
 
 
335
    if (delta>=0) {
 
336
      printf("-> first is better by %f dB\n",delta);
 
337
    }
 
338
    else {
 
339
      printf("-> second is better by %f dB\n",-delta);
 
340
    }
 
341
  }
 
342
 
 
343
  return 0;
 
344
}