f82ecb46a05b96f3fcf353bf8591f8f841aeac93
[theora.git] / lib / collect.c
1 /********************************************************************
2  *                                                                  *
3  * THIS FILE IS PART OF THE OggTheora SOFTWARE CODEC SOURCE CODE.   *
4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
7  *                                                                  *
8  * THE Theora SOURCE CODE IS COPYRIGHT (C) 2002-2009                *
9  * by the Xiph.Org Foundation http://www.xiph.org/                  *
10  *                                                                  *
11  ********************************************************************
12
13   function: mode selection code
14   last mod: $Id$
15
16  ********************************************************************/
17 #include <stdio.h>
18 #include <limits.h>
19 #include <math.h>
20 #include <string.h>
21 #include "collect.h"
22
23 #if defined(OC_COLLECT_METRICS)
24
25 int              OC_HAS_MODE_METRICS;
26 double           OC_MODE_RD_WEIGHT[OC_LOGQ_BINS][3][2][OC_SAD_BINS];
27 oc_mode_metrics  OC_MODE_METRICS[OC_LOGQ_BINS-1][3][2][OC_SAD_BINS];
28 const char      *OC_MODE_METRICS_FILENAME="modedec.stats";
29
30 void oc_mode_metrics_add(oc_mode_metrics *_metrics,
31  double _w,int _s,int _q,int _r,double _d){
32   if(_metrics->w>0){
33     double ds;
34     double dq;
35     double dr;
36     double dd;
37     double ds2;
38     double dq2;
39     double s2;
40     double sq;
41     double q2;
42     double sr;
43     double qr;
44     double sd;
45     double qd;
46     double s2q;
47     double sq2;
48     double w;
49     double wa;
50     double rwa;
51     double rwa2;
52     double rwb;
53     double rwb2;
54     double rw2;
55     double rw3;
56     double rw4;
57     wa=_metrics->w;
58     ds=_s-_metrics->s/wa;
59     dq=_q-_metrics->q/wa;
60     dr=_r-_metrics->r/wa;
61     dd=_d-_metrics->d/wa;
62     ds2=ds*ds;
63     dq2=dq*dq;
64     s2=_metrics->s2;
65     sq=_metrics->sq;
66     q2=_metrics->q2;
67     sr=_metrics->sr;
68     qr=_metrics->qr;
69     sd=_metrics->sd;
70     qd=_metrics->qd;
71     s2q=_metrics->s2q;
72     sq2=_metrics->sq2;
73     w=wa+_w;
74     rwa=wa/w;
75     rwb=_w/w;
76     rwa2=rwa*rwa;
77     rwb2=rwb*rwb;
78     rw2=wa*rwb;
79     rw3=rw2*(rwa2-rwb2);
80     rw4=_w*rwa2*rwa2+wa*rwb2*rwb2;
81     _metrics->s2q2+=-2*(ds*sq2+dq*s2q)*rwb
82      +(ds2*q2+4*ds*dq*sq+dq2*s2)*rwb2+ds2*dq2*rw4;
83     _metrics->s2q+=(-2*ds*sq-dq*s2)*rwb+ds2*dq*rw3;
84     _metrics->sq2+=(-ds*q2-2*dq*sq)*rwb+ds*dq2*rw3;
85     _metrics->sqr+=(-ds*qr-dq*sr-dr*sq)*rwb+ds*dq*dr*rw3;
86     _metrics->sqd+=(-ds*qd-dq*sd-dd*sq)*rwb+ds*dq*dd*rw3;
87     _metrics->s2+=ds2*rw2;
88     _metrics->sq+=ds*dq*rw2;
89     _metrics->q2+=dq2*rw2;
90     _metrics->sr+=ds*dr*rw2;
91     _metrics->qr+=dq*dr*rw2;
92     _metrics->r2+=dr*dr*rw2;
93     _metrics->sd+=ds*dd*rw2;
94     _metrics->qd+=dq*dd*rw2;
95     _metrics->d2+=dd*dd*rw2;
96   }
97   _metrics->w+=_w;
98   _metrics->s+=_s*_w;
99   _metrics->q+=_q*_w;
100   _metrics->r+=_r*_w;
101   _metrics->d+=_d*_w;
102 }
103
104 void oc_mode_metrics_merge(oc_mode_metrics *_dst,
105  const oc_mode_metrics *_src,int _n){
106   int i;
107   /*Find a non-empty set of metrics.*/
108   for(i=0;i<_n&&_src[i].w==0;i++);
109   if(i>=_n){
110     memset(_dst,0,sizeof(*_dst));
111     return;
112   }
113   memcpy(_dst,_src+i,sizeof(*_dst));
114   /*And iterate over the remaining non-empty sets of metrics.*/
115   for(i++;i<_n;i++)if(_src[i].w!=0){
116     double ds;
117     double dq;
118     double dr;
119     double dd;
120     double ds2;
121     double dq2;
122     double s2a;
123     double s2b;
124     double sqa;
125     double sqb;
126     double q2a;
127     double q2b;
128     double sra;
129     double srb;
130     double qra;
131     double qrb;
132     double sda;
133     double sdb;
134     double qda;
135     double qdb;
136     double s2qa;
137     double s2qb;
138     double sq2a;
139     double sq2b;
140     double w;
141     double wa;
142     double wb;
143     double rwa;
144     double rwb;
145     double rwa2;
146     double rwb2;
147     double rw2;
148     double rw3;
149     double rw4;
150     wa=_dst->w;
151     wb=_src[i].w;
152     ds=_src[i].s/wb-_dst->s/wa;
153     dq=_src[i].q/wb-_dst->q/wa;
154     dr=_src[i].r/wb-_dst->r/wa;
155     dd=_src[i].d/wb-_dst->d/wa;
156     ds2=ds*ds;
157     dq2=dq*dq;
158     s2a=_dst->s2;
159     sqa=_dst->sq;
160     q2a=_dst->q2;
161     sra=_dst->sr;
162     qra=_dst->qr;
163     sda=_dst->sd;
164     qda=_dst->qd;
165     s2qa=_dst->s2q;
166     sq2a=_dst->sq2;
167     s2b=_src[i].s2;
168     sqb=_src[i].sq;
169     q2b=_src[i].q2;
170     srb=_src[i].sr;
171     qrb=_src[i].qr;
172     sdb=_src[i].sd;
173     qdb=_src[i].qd;
174     s2qb=_src[i].s2q;
175     sq2b=_src[i].sq2;
176     w=wa+wb;
177     if(w==0)rwa=rwb=0;
178     else{
179       rwa=wa/w;
180       rwb=wb/w;
181     }
182     rwa2=rwa*rwa;
183     rwb2=rwb*rwb;
184     rw2=wa*rwb;
185     rw3=rw2*(rwa2-rwb2);
186     rw4=wb*rwa2*rwa2+wa*rwb2*rwb2;
187     /*
188     (1,1,1) ->
189      (0,0,0)#
190      (1,0,0) C(1,1)*C(1,0)*C(1,0)->  d^{1,0,0}*(rwa*B_{0,1,1}-rwb*A_{0,1,1})
191      (0,1,0) C(1,0)*C(1,1)*C(1,0)->  d^{0,1,0}*(rwa*B_{1,0,1}-rwb*A_{1,0,1})
192      (0,0,1) C(1,0)*C(1,0)*C(1,1)->  d^{0,0,1}*(rwa*B_{1,1,0}-rwb*A_{1,1,0})
193      (1,1,0)*
194      (1,0,1)*
195      (0,1,1)*
196      (1,1,1) C(1,1)*C(1,1)*C(1,1)->  d^{1,1,1}*(rwa^3*wb-rwb^3*wa)
197     (2,1) ->
198      (0,0)#
199      (1,0) C(2,1)*C(1,1)->2*d^{1,0}*(rwa*B_{1,1}-rwb*A_{1,1})
200      (0,1) C(2,0)*C(1,1)->  d^{0,1}*(rwa*B_{2,0}-rwb*A_{2,0})
201      (2,0)*
202      (1,1)*
203      (2,1) C(2,2)*C(1,1)->  d^{2,1}*(rwa^3*wb-rwb^3*wa)
204     (2,2) ->
205      (0,0)#
206      (1,0) C(2,1)*C(2,0)->2*d^{1,0}*(rwa*B_{1,2}-rwb*A_{1,2})
207      (0,1) C(2,0)*C(2,1)->2*d^{0,1}*(rwa*B_{2,1}-rwb*A_{2,1})
208      (2,0) C(2,2)*C(2,0)->  d^{2,0}*(rwa^2*B_{0,2}+rwb^2*A_{0,2})
209      (1,1) C(2,1)*C(2,1)->4*d^{1,1}*(rwa^2*B_{1,1}+rwb^2*A_{1,1})
210      (0,2) C(2,0)*C(2,2)->  d^{0,2}*(rwa^2*B_{2,0}+rwb^2*A_{2,0})
211      (1,2)*
212      (2,1)*
213      (2,2) C(2,2)*C(2,2)*d^{2,2}*(rwa^4*wb+rwb^4*wa)
214     */
215     _dst->s2q2+=_src[i].s2q2+2*(ds*(rwa*sq2b-rwb*sq2a)+dq*(rwa*s2qb-rwb*s2qa))
216      +ds2*(rwa2*q2b+rwb2*q2a)+4*ds*dq*(rwa2*sqb+rwb2*sqa)
217      +dq2*(rwa2*s2b+rwb2*s2a)+ds2*dq2*rw4;
218     _dst->s2q+=_src[i].s2q+2*ds*(rwa*sqb-rwb*sqa)
219      +dq*(rwa*s2b-rwb*s2a)+ds2*dq*rw3;
220     _dst->sq2+=_src[i].sq2+ds*(rwa*q2b-rwb*q2a)
221      +2*dq*(rwa*sqb-rwb*sqa)+ds*dq2*rw3;
222     _dst->sqr+=_src[i].sqr+ds*(rwa*qrb-rwb*qra)+dq*(rwa*srb-rwb*sra)
223      +dr*(rwa*sqb-rwb*sqa)+ds*dq*dr*rw3;
224     _dst->sqd+=_src[i].sqd+ds*(rwa*qdb-rwb*qda)+dq*(rwa*sdb-rwb*sda)
225      +dd*(rwa*sqb-rwb*sqa)+ds*dq*dd*rw3;
226     _dst->s2+=_src[i].s2+ds2*rw2;
227     _dst->sq+=_src[i].sq+ds*dq*rw2;
228     _dst->q2+=_src[i].q2+dq2*rw2;
229     _dst->sr+=_src[i].sr+ds*dr*rw2;
230     _dst->qr+=_src[i].qr+dq*dr*rw2;
231     _dst->r2+=_src[i].r2+dr*dr*rw2;
232     _dst->sd+=_src[i].sd+ds*dd*rw2;
233     _dst->qd+=_src[i].qd+dq*dd*rw2;
234     _dst->d2+=_src[i].d2+dd*dd*rw2;
235     _dst->w+=_src[i].w;
236     _dst->s+=_src[i].s;
237     _dst->q+=_src[i].q;
238     _dst->r+=_src[i].r;
239     _dst->d+=_src[i].d;
240   }
241 }
242
243 /*Adjust a single corner of a set of metric bins to minimize the squared
244    prediction error of R and D.
245   Each bin is assumed to cover a quad like so:
246     (s0,q0)    (s1,q0)
247        A----------B
248        |          |
249        |          |
250        |          |
251        |          |
252        C----------Z
253     (s0,q1)    (s1,q1)
254   The values A, B, and C are fixed, and Z is the free parameter.
255   Then, for example, R_i is predicted via bilinear interpolation as
256     x_i=(s_i-s0)/(s1-s0)
257     y_i=(q_i-q0)/(q1-q0)
258     dRds1_i=A+(B-A)*x_i
259     dRds2_i=C+(Z-C)*x_i
260     R_i=dRds1_i+(dRds2_i-dRds1_i)*y_i
261   To find the Z that minimizes the squared prediction error over i, this can
262    be rewritten as
263     R_i-(A+(B-A)*x_i+(C-A)*y_i+(A-B-C)*x_i*y_i)=x_i*y_i*Z
264   Letting X={...,x_i*y_i,...}^T and
265    Y={...,R_i-(A+(B-A)*x_i+(C-A)*y_i+(A-B-C)*x_i*y_i),...}^T,
266    the optimal Z is given by Z=(X^T.Y)/(X^T.X).
267   Now, we need to compute these dot products without actually storing data for
268    each sample.
269   Starting with X^T.X, we have
270    X^T.X = sum(x_i^2*y_i^2) = sum((s_i-s0)^2*(q_i-q0)^2)/((s1-s0)^2*(q1-q0)^2).
271   Expanding the interior of the sum in a monomial basis of s_i and q_i gives
272     s0^2*q0^2  *(1)
273      -2*s0*q0^2*(s_i)
274      -2*s0^2*q0*(q_i)
275      +q0^2     *(s_i^2)
276      +4*s0*q0  *(s_i*q_i)
277      +s0^2     *(q_i^2)
278      -2*q0     *(s_i^2*q_i)
279      -2*s0     *(s_i*q_i^2)
280      +1        *(s_i^2*q_i^2).
281   However, computing things directly in this basis leads to gross numerical
282    errors, as most of the terms will have similar size and destructive
283    cancellation results.
284   A much better basis is the central (co-)moment basis:
285     {1,s_i-sbar,q_i-qbar,(s_i-sbar)^2,(s_i-sbar)*(q_i-qbar),(q_i-qbar)^2,
286      (s_i-sbar)^2*(q_i-qbar),(s_i-sbar)*(q_i-qbar)^2,(s_i-sbar)^2*(q_i-qbar)^2},
287    where sbar and qbar are the average s and q values over the bin,
288    respectively.
289   In that basis, letting ds=sbar-s0 and dq=qbar-q0, (s_i-s0)^2*(q_i-q0)^2 is
290     ds^2*dq^2*(1)
291      +dq^2   *((s_i-sbar)^2)
292      +4*ds*dq*((s_i-sbar)*(q_i-qbar))
293      +ds^2   *((q_i-qbar)^2)
294      +2*dq   *((s_i-sbar)^2*(q_i-qbar))
295      +2*ds   *((s_i-sbar)*(q_i-qbar)^2)
296      +1      *((s_i-sbar)^2*(q_i-qbar)^2).
297   With these expressions in the central (co-)moment bases, all we need to do
298    is compute sums over the (co-)moment terms, which can be done
299    incrementally (see oc_mode_metrics_add() and oc_mode_metrics_merge()),
300    with no need to store the individual samples.
301   Now, for X^T.Y, we have
302     X^T.Y = sum((R_i-A-((B-A)/(s1-s0))*(s_i-s0)-((C-A)/(q1-q0))*(q_i-q0)
303      -((A-B-C)/((s1-s0)*(q1-q0)))*(s_i-s0)*(q_i-q0))*(s_i-s0)*(q_i-q0))/
304      ((s1-s0)*(q1-q0)),
305    or, rewriting the constants to simplify notation,
306     X^T.Y = sum((C0+C1*(s_i-s0)+C2*(q_i-q0)
307      +C3*(s_i-s0)*(q_i-q0)+R_i)*(s_i-s0)*(q_i-q0))/((s1-s0)*(q1-q0)).
308   Again, converting to the central (co-)moment basis, the interior of the
309    above sum is
310     ds*dq*(rbar+C0+C1*ds+C2*dq+C3*ds*dq)  *(1)
311      +(C1*dq+C3*dq^2)                     *((s_i-sbar)^2)
312      +(rbar+C0+2*C1*ds+2*C2*dq+4*C3*ds*dq)*((s_i-sbar)*(q_i-qbar))
313      +(C2*ds+C3*ds^2)                     *((q_i-qbar)^2)
314      +dq                                  *((s_i-sbar)*(r_i-rbar))
315      +ds                                  *((q_i-qbar)*(r_i-rbar))
316      +(C1+2*C3*dq)                        *((s_i-sbar)^2*(q_i-qbar))
317      +(C2+2*C3*ds)                        *((s_i-sbar)*(q_i-qbar)^2)
318      +1                                   *((s_i-sbar)*(q_i-qbar)*(r_i-rbar))
319      +C3                                  *((s_i-sbar)^2*(q_i-qbar)^2).
320   You might think it would be easier (if perhaps slightly less robust) to
321    accumulate terms directly around s0 and q0.
322   However, we update each corner of the bins in turn, so we would have to
323    change basis to move the sums from corner to corner anyway.*/
324 double oc_mode_metrics_solve(double *_r,double *_d,
325  const oc_mode_metrics *_metrics,const int *_s0,const int *_s1,
326  const int *_q0,const int *_q1,
327  const double *_ra,const double *_rb,const double *_rc,
328  const double *_da,const double *_db,const double *_dc,int _n){
329   double xx;
330   double rxy;
331   double dxy;
332   double wt;
333   int i;
334   xx=rxy=dxy=wt=0;
335   for(i=0;i<_n;i++)if(_metrics[i].w>0){
336     double s10;
337     double q10;
338     double sq10;
339     double ds;
340     double dq;
341     double ds2;
342     double dq2;
343     double r;
344     double d;
345     double s2;
346     double sq;
347     double q2;
348     double sr;
349     double qr;
350     double sd;
351     double qd;
352     double s2q;
353     double sq2;
354     double sqr;
355     double sqd;
356     double s2q2;
357     double c0;
358     double c1;
359     double c2;
360     double c3;
361     double w;
362     w=_metrics[i].w;
363     wt+=w;
364     s10=_s1[i]-_s0[i];
365     q10=_q1[i]-_q0[i];
366     sq10=s10*q10;
367     ds=_metrics[i].s/w-_s0[i];
368     dq=_metrics[i].q/w-_q0[i];
369     ds2=ds*ds;
370     dq2=dq*dq;
371     s2=_metrics[i].s2;
372     sq=_metrics[i].sq;
373     q2=_metrics[i].q2;
374     s2q=_metrics[i].s2q;
375     sq2=_metrics[i].sq2;
376     s2q2=_metrics[i].s2q2;
377     xx+=(dq2*(ds2*w+s2)+4*ds*dq*sq+ds2*q2+2*(dq*s2q+ds*sq2)+s2q2)/(sq10*sq10);
378     r=_metrics[i].r/w;
379     sr=_metrics[i].sr;
380     qr=_metrics[i].qr;
381     sqr=_metrics[i].sqr;
382     c0=-_ra[i];
383     c1=-(_rb[i]-_ra[i])/s10;
384     c2=-(_rc[i]-_ra[i])/q10;
385     c3=-(_ra[i]-_rb[i]-_rc[i])/sq10;
386     rxy+=(ds*dq*(r+c0+c1*ds+c2*dq+c3*ds*dq)*w+(c1*dq+c3*dq2)*s2
387      +(r+c0+2*(c1*ds+(c2+2*c3*ds)*dq))*sq+(c2*ds+c3*ds2)*q2+dq*sr+ds*qr
388      +(c1+2*c3*dq)*s2q+(c2+2*c3*ds)*sq2+sqr+c3*s2q2)/sq10;
389     d=_metrics[i].d/w;
390     sd=_metrics[i].sd;
391     qd=_metrics[i].qd;
392     sqd=_metrics[i].sqd;
393     c0=-_da[i];
394     c1=-(_db[i]-_da[i])/s10;
395     c2=-(_dc[i]-_da[i])/q10;
396     c3=-(_da[i]-_db[i]-_dc[i])/sq10;
397     dxy+=(ds*dq*(d+c0+c1*ds+c2*dq+c3*ds*dq)*w+(c1*dq+c3*dq2)*s2
398      +(d+c0+2*(c1*ds+(c2+2*c3*ds)*dq))*sq+(c2*ds+c3*ds2)*q2+dq*sd+ds*qd
399      +(c1+2*c3*dq)*s2q+(c2+2*c3*ds)*sq2+sqd+c3*s2q2)/sq10;
400   }
401   if(xx>1E-3){
402     *_r=rxy/xx;
403     *_d=dxy/xx;
404   }
405   else{
406     *_r=0;
407     *_d=0;
408   }
409   return wt;
410 }
411
412 /*Compile collected SATD/logq/rate/RMSE metrics into a form that's immediately
413    useful for mode decision.*/
414 void oc_mode_metrics_update(int _niters_min,int _reweight){
415   int niters;
416   int prevdr;
417   int prevdd;
418   int dr;
419   int dd;
420   int pli;
421   int qti;
422   int qi;
423   int si;
424   dd=dr=INT_MAX;
425   niters=0;
426   /*The encoder interpolates rate and RMSE terms bilinearly from an
427      OC_LOGQ_BINS by OC_SAD_BINS grid of sample points in OC_MODE_RD.
428     To find the sample values at the grid points that minimize the total
429      squared prediction error actually requires solving a relatively sparse
430      linear system with a number of variables equal to the number of grid
431      points.
432     Instead of writing a general sparse linear system solver, we just use
433      Gauss-Seidel iteration, i.e., we update one grid point at time until
434      they stop changing.*/
435   do{
436     prevdr=dr;
437     prevdd=dd;
438     dd=dr=0;
439     for(pli=0;pli<3;pli++){
440       for(qti=0;qti<2;qti++){
441         for(qi=0;qi<OC_LOGQ_BINS;qi++){
442           for(si=0;si<OC_SAD_BINS;si++){
443             oc_mode_metrics m[4];
444             int             s0[4];
445             int             s1[4];
446             int             q0[4];
447             int             q1[4];
448             double          ra[4];
449             double          rb[4];
450             double          rc[4];
451             double          da[4];
452             double          db[4];
453             double          dc[4];
454             double          r;
455             double          d;
456             int             rate;
457             int             rmse;
458             int             ds;
459             int             n;
460             n=0;
461             /*Collect the statistics for the (up to) four bins grid point
462                (si,qi) touches.*/
463             if(qi>0&&si>0){
464               q0[n]=OC_MODE_LOGQ[qi-1][pli][qti];
465               q1[n]=OC_MODE_LOGQ[qi][pli][qti];
466               s0[n]=si-1<<OC_SAD_SHIFT;
467               s1[n]=si<<OC_SAD_SHIFT;
468               ra[n]=ldexp(OC_MODE_RD[qi-1][pli][qti][si-1].rate,-OC_BIT_SCALE);
469               da[n]=ldexp(OC_MODE_RD[qi-1][pli][qti][si-1].rmse,-OC_RMSE_SCALE);
470               rb[n]=ldexp(OC_MODE_RD[qi-1][pli][qti][si].rate,-OC_BIT_SCALE);
471               db[n]=ldexp(OC_MODE_RD[qi-1][pli][qti][si].rmse,-OC_RMSE_SCALE);
472               rc[n]=ldexp(OC_MODE_RD[qi][pli][qti][si-1].rate,-OC_BIT_SCALE);
473               dc[n]=ldexp(OC_MODE_RD[qi][pli][qti][si-1].rmse,-OC_RMSE_SCALE);
474               *(m+n++)=*(OC_MODE_METRICS[qi-1][pli][qti]+si-1);
475             }
476             if(qi>0){
477               ds=si+1<OC_SAD_BINS?1:-1;
478               q0[n]=OC_MODE_LOGQ[qi-1][pli][qti];
479               q1[n]=OC_MODE_LOGQ[qi][pli][qti];
480               s0[n]=si+ds<<OC_SAD_SHIFT;
481               s1[n]=si<<OC_SAD_SHIFT;
482               ra[n]=ldexp(OC_MODE_RD[qi-1][pli][qti][si+ds].rate,-OC_BIT_SCALE);
483               da[n]=
484                ldexp(OC_MODE_RD[qi-1][pli][qti][si+ds].rmse,-OC_RMSE_SCALE);
485               rb[n]=ldexp(OC_MODE_RD[qi-1][pli][qti][si].rate,-OC_BIT_SCALE);
486               db[n]=ldexp(OC_MODE_RD[qi-1][pli][qti][si].rmse,-OC_RMSE_SCALE);
487               rc[n]=ldexp(OC_MODE_RD[qi][pli][qti][si+ds].rate,-OC_BIT_SCALE);
488               dc[n]=ldexp(OC_MODE_RD[qi][pli][qti][si+ds].rmse,-OC_RMSE_SCALE);
489               *(m+n++)=*(OC_MODE_METRICS[qi-1][pli][qti]+si);
490             }
491             if(qi+1<OC_LOGQ_BINS&&si>0){
492               q0[n]=OC_MODE_LOGQ[qi+1][pli][qti];
493               q1[n]=OC_MODE_LOGQ[qi][pli][qti];
494               s0[n]=si-1<<OC_SAD_SHIFT;
495               s1[n]=si<<OC_SAD_SHIFT;
496               ra[n]=ldexp(OC_MODE_RD[qi+1][pli][qti][si-1].rate,-OC_BIT_SCALE);
497               da[n]=ldexp(OC_MODE_RD[qi+1][pli][qti][si-1].rmse,-OC_RMSE_SCALE);
498               rb[n]=ldexp(OC_MODE_RD[qi+1][pli][qti][si].rate,-OC_BIT_SCALE);
499               db[n]=ldexp(OC_MODE_RD[qi+1][pli][qti][si].rmse,-OC_RMSE_SCALE);
500               rc[n]=ldexp(OC_MODE_RD[qi][pli][qti][si-1].rate,-OC_BIT_SCALE);
501               dc[n]=ldexp(OC_MODE_RD[qi][pli][qti][si-1].rmse,-OC_RMSE_SCALE);
502               *(m+n++)=*(OC_MODE_METRICS[qi][pli][qti]+si-1);
503             }
504             if(qi+1<OC_LOGQ_BINS){
505               ds=si+1<OC_SAD_BINS?1:-1;
506               q0[n]=OC_MODE_LOGQ[qi+1][pli][qti];
507               q1[n]=OC_MODE_LOGQ[qi][pli][qti];
508               s0[n]=si+ds<<OC_SAD_SHIFT;
509               s1[n]=si<<OC_SAD_SHIFT;
510               ra[n]=ldexp(OC_MODE_RD[qi+1][pli][qti][si+ds].rate,-OC_BIT_SCALE);
511               da[n]=
512                ldexp(OC_MODE_RD[qi+1][pli][qti][si+ds].rmse,-OC_RMSE_SCALE);
513               rb[n]=ldexp(OC_MODE_RD[qi+1][pli][qti][si].rate,-OC_BIT_SCALE);
514               db[n]=ldexp(OC_MODE_RD[qi+1][pli][qti][si].rmse,-OC_RMSE_SCALE);
515               rc[n]=ldexp(OC_MODE_RD[qi][pli][qti][si+ds].rate,-OC_BIT_SCALE);
516               dc[n]=ldexp(OC_MODE_RD[qi][pli][qti][si+ds].rmse,-OC_RMSE_SCALE);
517               *(m+n++)=*(OC_MODE_METRICS[qi][pli][qti]+si);
518             }
519             /*On the first pass, initialize with a simple weighted average of
520                the neighboring bins.*/
521             if(!OC_HAS_MODE_METRICS&&niters==0){
522               double w;
523               w=r=d=0;
524               while(n-->0){
525                 w+=m[n].w;
526                 r+=m[n].r;
527                 d+=m[n].d;
528               }
529               r=w>1E-3?r/w:0;
530               d=w>1E-3?d/w:0;
531               OC_MODE_RD_WEIGHT[qi][pli][qti][si]=w;
532             }
533             else{
534               /*Update the grid point and save the weight for later.*/
535               OC_MODE_RD_WEIGHT[qi][pli][qti][si]=
536                oc_mode_metrics_solve(&r,&d,m,s0,s1,q0,q1,ra,rb,rc,da,db,dc,n);
537             }
538             rate=OC_CLAMPI(-32768,(int)(ldexp(r,OC_BIT_SCALE)+0.5),32767);
539             rmse=OC_CLAMPI(-32768,(int)(ldexp(d,OC_RMSE_SCALE)+0.5),32767);
540             dr+=abs(rate-OC_MODE_RD[qi][pli][qti][si].rate);
541             dd+=abs(rmse-OC_MODE_RD[qi][pli][qti][si].rmse);
542             OC_MODE_RD[qi][pli][qti][si].rate=(ogg_int16_t)rate;
543             OC_MODE_RD[qi][pli][qti][si].rmse=(ogg_int16_t)rmse;
544           }
545         }
546       }
547     }
548   }
549   /*After a fixed number of initial iterations, only iterate so long as the
550      total change is decreasing.
551     This ensures we don't oscillate forever, which is a danger, as all of our
552      results are rounded fairly coarsely.*/
553   while((dr>0||dd>0)&&(niters++<_niters_min||(dr<prevdr&&dd<prevdd)));
554   if(_reweight){
555     /*Now, reduce the values of the optimal solution until we get enough
556        samples in each bin to overcome the constant OC_ZWEIGHT factor.
557       This encourages sampling under-populated bins and prevents a single large
558        sample early on from discouraging coding in that bin ever again.*/
559     for(pli=0;pli<3;pli++){ 
560       for(qti=0;qti<2;qti++){
561         for(qi=0;qi<OC_LOGQ_BINS;qi++){
562           for(si=0;si<OC_SAD_BINS;si++){
563             double wt;
564             wt=OC_MODE_RD_WEIGHT[qi][pli][qti][si];
565             wt/=OC_ZWEIGHT+wt;
566             OC_MODE_RD[qi][pli][qti][si].rate=(ogg_int16_t)
567              (OC_MODE_RD[qi][pli][qti][si].rate*wt+0.5);
568             OC_MODE_RD[qi][pli][qti][si].rmse=(ogg_int16_t)
569              (OC_MODE_RD[qi][pli][qti][si].rmse*wt+0.5);
570           }
571         }
572       }
573     }
574   }
575 }
576
577 void oc_mode_metrics_dump(void){
578   FILE *fmetrics;
579   fmetrics=fopen(OC_MODE_METRICS_FILENAME,"wb");
580   if(fmetrics!=NULL){
581     (void)fwrite(OC_MODE_METRICS,sizeof(OC_MODE_METRICS),1,fmetrics);
582     (void)fwrite(OC_MODE_LOGQ,sizeof(OC_MODE_LOGQ),1,fmetrics);
583     fclose(fmetrics);
584   }
585 }
586
587 void oc_mode_metrics_print(FILE *_fout){
588   int qii;
589   fprintf(_fout,
590    "/*File generated by libtheora with OC_COLLECT_METRICS"
591    " defined at compile time.*/\n"
592    "#if !defined(_modedec_H)\n"
593    "# define _modedec_H (1)\n"
594    "# include \"encint.h\"\n"
595    "\n"
596    "\n"
597    "\n"
598    "/*The log of the average quantizer for each of the OC_MODE_RD table rows\n"
599    "   (e.g., for the represented qi's, and each pli and qti), in Q10 format.\n"
600    "  The actual statistics used by the encoder will be interpolated from\n"
601    "   that table based on log_plq for the actual quantization matrix used.*/\n"
602    "# if !defined(OC_COLLECT_METRICS)\n"
603    "static const\n"
604    "# endif\n"
605    "ogg_int16_t OC_MODE_LOGQ[OC_LOGQ_BINS][3][2]={\n");
606   for(qii=0;qii<OC_LOGQ_BINS;qii++){
607     fprintf(_fout,"  { {0x%04X,0x%04X},{0x%04X,0x%04X},{0x%04X,0x%04X} }%s\n",
608      OC_MODE_LOGQ[qii][0][0],OC_MODE_LOGQ[qii][0][1],OC_MODE_LOGQ[qii][1][0],
609      OC_MODE_LOGQ[qii][1][1],OC_MODE_LOGQ[qii][2][0],OC_MODE_LOGQ[qii][2][1],
610      qii+1<OC_LOGQ_BINS?",":"");
611   }
612   fprintf(_fout,
613    "};\n"
614    "\n"
615    "# if !defined(OC_COLLECT_METRICS)\n"
616    "static const\n"
617    "# endif\n"
618    "oc_mode_rd OC_MODE_RD[OC_LOGQ_BINS][3][2][OC_SAD_BINS]={\n");
619   for(qii=0;qii<OC_LOGQ_BINS;qii++){
620     int pli;
621     fprintf(_fout,"  {\n");
622     for(pli=0;pli<3;pli++){
623       int qti;
624       fprintf(_fout,"    {\n");
625       for(qti=0;qti<2;qti++){
626         int bin;
627         int qi;
628         static const char *pl_names[3]={"Y'","Cb","Cr"};
629         static const char *qti_names[2]={"INTRA","INTER"};
630         qi=(63*qii+(OC_LOGQ_BINS-1>>1))/(OC_LOGQ_BINS-1);
631         fprintf(_fout,"      /*%s  qi=%i  %s*/\n",
632          pl_names[pli],qi,qti_names[qti]);
633         fprintf(_fout,"      {\n");
634         fprintf(_fout,"        ");
635         for(bin=0;bin<OC_SAD_BINS;bin++){
636           if(bin&&!(bin&0x3))fprintf(_fout,"\n        ");
637           fprintf(_fout,"{%5i,%5i}",
638            OC_MODE_RD[qii][pli][qti][bin].rate,
639            OC_MODE_RD[qii][pli][qti][bin].rmse);
640           if(bin+1<OC_SAD_BINS)fprintf(_fout,",");
641         }
642         fprintf(_fout,"\n      }");
643         if(qti<1)fprintf(_fout,",");
644         fprintf(_fout,"\n");
645       }
646       fprintf(_fout,"    }");
647       if(pli<2)fprintf(_fout,",");
648       fprintf(_fout,"\n");
649     }
650     fprintf(_fout,"  }");
651     if(qii+1<OC_LOGQ_BINS)fprintf(_fout,",");
652     fprintf(_fout,"\n");
653   }
654   fprintf(_fout,
655    "};\n"
656    "\n"
657    "#endif\n");
658 }
659
660
661 # if !defined(OC_COLLECT_NO_ENC_FUNCS)
662 void oc_enc_mode_metrics_load(oc_enc_ctx *_enc){
663   oc_restore_fpu(&_enc->state);
664   /*Load any existing mode metrics if we haven't already.*/
665   if(!OC_HAS_MODE_METRICS){
666     FILE *fmetrics;
667     memset(OC_MODE_METRICS,0,sizeof(OC_MODE_METRICS));
668     fmetrics=fopen(OC_MODE_METRICS_FILENAME,"rb");
669     if(fmetrics!=NULL){
670       (void)fread(OC_MODE_METRICS,sizeof(OC_MODE_METRICS),1,fmetrics);
671       (void)fread(OC_MODE_LOGQ,sizeof(OC_MODE_LOGQ),1,fmetrics);
672       fclose(fmetrics);
673     }
674     else{
675       int qii;
676       int qi;
677       int pli;
678       int qti;
679       for(qii=0;qii<OC_LOGQ_BINS;qii++){
680         qi=(63*qii+(OC_LOGQ_BINS-1>>1))/(OC_LOGQ_BINS-1);
681         for(pli=0;pli<3;pli++)for(qti=0;qti<2;qti++){
682           OC_MODE_LOGQ[qii][pli][qti]=_enc->log_plq[qi][pli][qti];
683         }
684       }
685     }
686     oc_mode_metrics_update(100,1);
687     OC_HAS_MODE_METRICS=1;
688   }
689 }
690
691 /*The following token skipping code used to also be used in the decoder (and
692    even at one point other places in the encoder).
693   However, it was obsoleted by other optimizations, and is now only used here.
694   It has been moved here to avoid generating the code when it's not needed.*/
695
696 /*Determines the number of blocks or coefficients to be skipped for a given
697    token value.
698   _token:      The token value to skip.
699   _extra_bits: The extra bits attached to this token.
700   Return: A positive value indicates that number of coefficients are to be
701            skipped in the current block.
702           Otherwise, the negative of the return value indicates that number of
703            blocks are to be ended.*/
704 typedef ptrdiff_t (*oc_token_skip_func)(int _token,int _extra_bits);
705
706 /*Handles the simple end of block tokens.*/
707 static ptrdiff_t oc_token_skip_eob(int _token,int _extra_bits){
708   int nblocks_adjust;
709   nblocks_adjust=OC_UNIBBLE_TABLE32(0,1,2,3,7,15,0,0,_token)+1;
710   return -_extra_bits-nblocks_adjust;
711 }
712
713 /*The last EOB token has a special case, where an EOB run of size zero ends all
714    the remaining blocks in the frame.*/
715 static ptrdiff_t oc_token_skip_eob6(int _token,int _extra_bits){
716   /*Note: We want to return -PTRDIFF_MAX, but that requires C99, which is not
717      yet available everywhere; this should be equivalent.*/
718   if(!_extra_bits)return -(~(size_t)0>>1);
719   return -_extra_bits;
720 }
721
722 /*Handles the pure zero run tokens.*/
723 static ptrdiff_t oc_token_skip_zrl(int _token,int _extra_bits){
724   return _extra_bits+1;
725 }
726
727 /*Handles a normal coefficient value token.*/
728 static ptrdiff_t oc_token_skip_val(void){
729   return 1;
730 }
731
732 /*Handles a category 1A zero run/coefficient value combo token.*/
733 static ptrdiff_t oc_token_skip_run_cat1a(int _token){
734   return _token-OC_DCT_RUN_CAT1A+2;
735 }
736
737 /*Handles category 1b, 1c, 2a, and 2b zero run/coefficient value combo tokens.*/
738 static ptrdiff_t oc_token_skip_run(int _token,int _extra_bits){
739   int run_cati;
740   int ncoeffs_mask;
741   int ncoeffs_adjust;
742   run_cati=_token-OC_DCT_RUN_CAT1B;
743   ncoeffs_mask=OC_BYTE_TABLE32(3,7,0,1,run_cati);
744   ncoeffs_adjust=OC_BYTE_TABLE32(7,11,2,3,run_cati);
745   return (_extra_bits&ncoeffs_mask)+ncoeffs_adjust;
746 }
747
748 /*A jump table for computing the number of coefficients or blocks to skip for
749    a given token value.
750   This reduces all the conditional branches, etc., needed to parse these token
751    values down to one indirect jump.*/
752 static const oc_token_skip_func OC_TOKEN_SKIP_TABLE[TH_NDCT_TOKENS]={
753   oc_token_skip_eob,
754   oc_token_skip_eob,
755   oc_token_skip_eob,
756   oc_token_skip_eob,
757   oc_token_skip_eob,
758   oc_token_skip_eob,
759   oc_token_skip_eob6,
760   oc_token_skip_zrl,
761   oc_token_skip_zrl,
762   (oc_token_skip_func)oc_token_skip_val,
763   (oc_token_skip_func)oc_token_skip_val,
764   (oc_token_skip_func)oc_token_skip_val,
765   (oc_token_skip_func)oc_token_skip_val,
766   (oc_token_skip_func)oc_token_skip_val,
767   (oc_token_skip_func)oc_token_skip_val,
768   (oc_token_skip_func)oc_token_skip_val,
769   (oc_token_skip_func)oc_token_skip_val,
770   (oc_token_skip_func)oc_token_skip_val,
771   (oc_token_skip_func)oc_token_skip_val,
772   (oc_token_skip_func)oc_token_skip_val,
773   (oc_token_skip_func)oc_token_skip_val,
774   (oc_token_skip_func)oc_token_skip_val,
775   (oc_token_skip_func)oc_token_skip_val,
776   (oc_token_skip_func)oc_token_skip_run_cat1a,
777   (oc_token_skip_func)oc_token_skip_run_cat1a,
778   (oc_token_skip_func)oc_token_skip_run_cat1a,
779   (oc_token_skip_func)oc_token_skip_run_cat1a,
780   (oc_token_skip_func)oc_token_skip_run_cat1a,
781   oc_token_skip_run,
782   oc_token_skip_run,
783   oc_token_skip_run,
784   oc_token_skip_run
785 };
786
787 /*Determines the number of blocks or coefficients to be skipped for a given
788    token value.
789   _token:      The token value to skip.
790   _extra_bits: The extra bits attached to this token.
791   Return: A positive value indicates that number of coefficients are to be
792            skipped in the current block.
793           Otherwise, the negative of the return value indicates that number of
794            blocks are to be ended.
795           0 will never be returned, so that at least one coefficient in one
796            block will always be decoded for every token.*/
797 static ptrdiff_t oc_dct_token_skip(int _token,int _extra_bits){
798   return (*OC_TOKEN_SKIP_TABLE[_token])(_token,_extra_bits);
799 }
800
801
802 void oc_enc_mode_metrics_collect(oc_enc_ctx *_enc){
803   static const unsigned char OC_ZZI_HUFF_OFFSET[64]={
804      0,16,16,16,16,16,32,32,
805     32,32,32,32,32,32,32,48,
806     48,48,48,48,48,48,48,48,
807     48,48,48,48,64,64,64,64,
808     64,64,64,64,64,64,64,64,
809     64,64,64,64,64,64,64,64,
810     64,64,64,64,64,64,64,64
811   };
812   const oc_fragment *frags;
813   const unsigned    *frag_satd;
814   const unsigned    *frag_ssd;
815   const ptrdiff_t   *coded_fragis;
816   ptrdiff_t          ncoded_fragis;
817   ptrdiff_t          fragii;
818   double             fragw;
819   int                modelines[3][3][2];
820   int                qti;
821   int                qii;
822   int                qi;
823   int                pli;
824   int                zzi;
825   int                token;
826   int                eb;
827   oc_restore_fpu(&_enc->state);
828   /*Figure out which metric bins to use for this frame's quantizers.*/
829   for(qii=0;qii<_enc->state.nqis;qii++){
830     for(pli=0;pli<3;pli++){
831       for(qti=0;qti<2;qti++){
832         int log_plq;
833         int modeline;
834         log_plq=_enc->log_plq[_enc->state.qis[qii]][pli][qti];
835         for(modeline=0;modeline<OC_LOGQ_BINS-1&&
836          OC_MODE_LOGQ[modeline+1][pli][qti]>log_plq;modeline++);
837         modelines[qii][pli][qti]=modeline;
838       }
839     }
840   }
841   qti=_enc->state.frame_type;
842   frags=_enc->state.frags;
843   frag_satd=_enc->frag_satd;
844   frag_ssd=_enc->frag_ssd;
845   coded_fragis=_enc->state.coded_fragis;
846   ncoded_fragis=fragii=0;
847   /*Weight the fragments by the inverse frame size; this prevents HD content
848      from dominating the statistics.*/
849   fragw=1.0/_enc->state.nfrags;
850   for(pli=0;pli<3;pli++){
851     ptrdiff_t ti[64];
852     int       eob_token[64];
853     int       eob_run[64];
854     /*Set up token indices and eob run counts.
855       We don't bother trying to figure out the real cost of the runs that span
856        coefficients; instead we use the costs that were available when R-D
857        token optimization was done.*/
858     for(zzi=0;zzi<64;zzi++){
859       ti[zzi]=_enc->dct_token_offs[pli][zzi];
860       if(ti[zzi]>0){
861         token=_enc->dct_tokens[pli][zzi][0];
862         eb=_enc->extra_bits[pli][zzi][0];
863         eob_token[zzi]=token;
864         eob_run[zzi]=-oc_dct_token_skip(token,eb);
865       }
866       else{
867         eob_token[zzi]=OC_NDCT_EOB_TOKEN_MAX;
868         eob_run[zzi]=0;
869       }
870     }
871     /*Scan the list of coded fragments for this plane.*/
872     ncoded_fragis+=_enc->state.ncoded_fragis[pli];
873     for(;fragii<ncoded_fragis;fragii++){
874       ptrdiff_t fragi;
875       int       frag_bits;
876       int       huffi;
877       int       skip;
878       int       mb_mode;
879       unsigned  satd;
880       int       bin;
881       int       qtj;
882       fragi=coded_fragis[fragii];
883       frag_bits=0;
884       for(zzi=0;zzi<64;){
885         if(eob_run[zzi]>0){
886           /*We've reached the end of the block.*/
887           eob_run[zzi]--;
888           break;
889         }
890         huffi=_enc->huff_idxs[qti][zzi>0][pli+1>>1]
891          +OC_ZZI_HUFF_OFFSET[zzi];
892         if(eob_token[zzi]<OC_NDCT_EOB_TOKEN_MAX){
893           /*This token caused an EOB run to be flushed.
894             Therefore it gets the bits associated with it.*/
895           frag_bits+=_enc->huff_codes[huffi][eob_token[zzi]].nbits
896            +OC_DCT_TOKEN_EXTRA_BITS[eob_token[zzi]];
897           eob_token[zzi]=OC_NDCT_EOB_TOKEN_MAX;
898         }
899         token=_enc->dct_tokens[pli][zzi][ti[zzi]];
900         eb=_enc->extra_bits[pli][zzi][ti[zzi]];
901         ti[zzi]++;
902         skip=oc_dct_token_skip(token,eb);
903         if(skip<0){
904           eob_token[zzi]=token;
905           eob_run[zzi]=-skip;
906         }
907         else{
908           /*A regular DCT value token; accumulate the bits for it.*/
909           frag_bits+=_enc->huff_codes[huffi][token].nbits
910            +OC_DCT_TOKEN_EXTRA_BITS[token];
911           zzi+=skip;
912         }
913       }
914       mb_mode=frags[fragi].mb_mode;
915       qii=frags[fragi].qii;
916       qi=_enc->state.qis[qii];
917       satd=frag_satd[fragi]<<(pli+1&2);
918       bin=OC_MINI(satd>>OC_SAD_SHIFT,OC_SAD_BINS-1);
919       qtj=mb_mode!=OC_MODE_INTRA;
920       /*Accumulate statistics.
921         The rate (frag_bits) and RMSE (sqrt(frag_ssd)) are not scaled by
922          OC_BIT_SCALE and OC_RMSE_SCALE; this lets us change the scale factor
923          yet still use old data.*/
924       oc_mode_metrics_add(
925        OC_MODE_METRICS[modelines[qii][pli][qtj]][pli][qtj]+bin,
926        fragw,satd,_enc->log_plq[qi][pli][qtj],frag_bits,sqrt(frag_ssd[fragi]));
927     }
928   }
929   /*Update global SATD/logq/rate/RMSE estimation matrix.*/
930   oc_mode_metrics_update(4,1);
931 }
932 # endif
933
934 #endif