1/********************************************************************
2 *                                                                  *
3 * THIS FILE IS PART OF THE OggVorbis 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 OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
9 * by the Xiph.Org Foundation http://www.xiph.org/                  *
10 *                                                                  *
11 ********************************************************************
12
13 function: floor backend 1 implementation
14 last mod: $Id: floor1.c 17079 2010-03-26 06:51:41Z xiphmont $
15
16 ********************************************************************/
17
18#include <stdlib.h>
19#include <string.h>
20#include <math.h>
21#include <ogg/ogg.h>
22#include "vorbis/codec.h"
23#include "codec_internal.h"
24#include "registry.h"
25#include "codebook.h"
26#include "misc.h"
27#include "scales.h"
28
29#include <stdio.h>
30
31#define floor1_rangedB 140 /* floor 1 fixed at -140dB to 0dB range */
32
33typedef struct lsfit_acc{
34  int x0;
35  int x1;
36
37  int xa;
38  int ya;
39  int x2a;
40  int y2a;
41  int xya;
42  int an;
43
44  int xb;
45  int yb;
46  int x2b;
47  int y2b;
48  int xyb;
49  int bn;
50} lsfit_acc;
51
52/***********************************************/
53
54static void floor1_free_info(vorbis_info_floor *i){
55  vorbis_info_floor1 *info=(vorbis_info_floor1 *)i;
56  if(info){
57    memset(info,0,sizeof(*info));
58    _ogg_free(info);
59  }
60}
61
62static void floor1_free_look(vorbis_look_floor *i){
63  vorbis_look_floor1 *look=(vorbis_look_floor1 *)i;
64  if(look){
65    /*fprintf(stderr,"floor 1 bit usage %f:%f (%f total)\n",
66            (float)look->phrasebits/look->frames,
67            (float)look->postbits/look->frames,
68            (float)(look->postbits+look->phrasebits)/look->frames);*/
69
70    memset(look,0,sizeof(*look));
71    _ogg_free(look);
72  }
73}
74
75static int ilog(unsigned int v){
76  int ret=0;
77  while(v){
78    ret++;
79    v>>=1;
80  }
81  return(ret);
82}
83
84static int ilog2(unsigned int v){
85  int ret=0;
86  if(v)--v;
87  while(v){
88    ret++;
89    v>>=1;
90  }
91  return(ret);
92}
93
94static void floor1_pack (vorbis_info_floor *i,oggpack_buffer *opb){
95  vorbis_info_floor1 *info=(vorbis_info_floor1 *)i;
96  int j,k;
97  int count=0;
98  int rangebits;
99  int maxposit=info->postlist[1];
100  int maxclass=-1;
101
102  /* save out partitions */
103  oggpack_write(opb,info->partitions,5); /* only 0 to 31 legal */
104  for(j=0;j<info->partitions;j++){
105    oggpack_write(opb,info->partitionclass[j],4); /* only 0 to 15 legal */
106    if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j];
107  }
108
109  /* save out partition classes */
110  for(j=0;j<maxclass+1;j++){
111    oggpack_write(opb,info->class_dim[j]-1,3); /* 1 to 8 */
112    oggpack_write(opb,info->class_subs[j],2); /* 0 to 3 */
113    if(info->class_subs[j])oggpack_write(opb,info->class_book[j],8);
114    for(k=0;k<(1<<info->class_subs[j]);k++)
115      oggpack_write(opb,info->class_subbook[j][k]+1,8);
116  }
117
118  /* save out the post list */
119  oggpack_write(opb,info->mult-1,2);     /* only 1,2,3,4 legal now */
120  oggpack_write(opb,ilog2(maxposit),4);
121  rangebits=ilog2(maxposit);
122
123  for(j=0,k=0;j<info->partitions;j++){
124    count+=info->class_dim[info->partitionclass[j]];
125    for(;k<count;k++)
126      oggpack_write(opb,info->postlist[k+2],rangebits);
127  }
128}
129
130static int icomp(const void *a,const void *b){
131  return(**(int **)a-**(int **)b);
132}
133
134static vorbis_info_floor *floor1_unpack (vorbis_info *vi,oggpack_buffer *opb){
135  codec_setup_info     *ci=vi->codec_setup;
136  int j,k,count=0,maxclass=-1,rangebits;
137
138  vorbis_info_floor1 *info=_ogg_calloc(1,sizeof(*info));
139  /* read partitions */
140  info->partitions=oggpack_read(opb,5); /* only 0 to 31 legal */
141  for(j=0;j<info->partitions;j++){
142    info->partitionclass[j]=oggpack_read(opb,4); /* only 0 to 15 legal */
143    if(info->partitionclass[j]<0)goto err_out;
144    if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j];
145  }
146
147  /* read partition classes */
148  for(j=0;j<maxclass+1;j++){
149    info->class_dim[j]=oggpack_read(opb,3)+1; /* 1 to 8 */
150    info->class_subs[j]=oggpack_read(opb,2); /* 0,1,2,3 bits */
151    if(info->class_subs[j]<0)
152      goto err_out;
153    if(info->class_subs[j])info->class_book[j]=oggpack_read(opb,8);
154    if(info->class_book[j]<0 || info->class_book[j]>=ci->books)
155      goto err_out;
156    for(k=0;k<(1<<info->class_subs[j]);k++){
157      info->class_subbook[j][k]=oggpack_read(opb,8)-1;
158      if(info->class_subbook[j][k]<-1 || info->class_subbook[j][k]>=ci->books)
159        goto err_out;
160    }
161  }
162
163  /* read the post list */
164  info->mult=oggpack_read(opb,2)+1;     /* only 1,2,3,4 legal now */
165  rangebits=oggpack_read(opb,4);
166  if(rangebits<0)goto err_out;
167
168  for(j=0,k=0;j<info->partitions;j++){
169    count+=info->class_dim[info->partitionclass[j]];
170    for(;k<count;k++){
171      int t=info->postlist[k+2]=oggpack_read(opb,rangebits);
172      if(t<0 || t>=(1<<rangebits))
173        goto err_out;
174    }
175  }
176  info->postlist[0]=0;
177  info->postlist[1]=1<<rangebits;
178
179  /* don't allow repeated values in post list as they'd result in
180     zero-length segments */
181  {
182    int *sortpointer[VIF_POSIT+2];
183    for(j=0;j<count+2;j++)sortpointer[j]=info->postlist+j;
184    qsort(sortpointer,count+2,sizeof(*sortpointer),icomp);
185
186    for(j=1;j<count+2;j++)
187      if(*sortpointer[j-1]==*sortpointer[j])goto err_out;
188  }
189
190  return(info);
191
192 err_out:
193  floor1_free_info(info);
194  return(NULL);
195}
196
197static vorbis_look_floor *floor1_look(vorbis_dsp_state *vd,
198                                      vorbis_info_floor *in){
199
200  int *sortpointer[VIF_POSIT+2];
201  vorbis_info_floor1 *info=(vorbis_info_floor1 *)in;
202  vorbis_look_floor1 *look=_ogg_calloc(1,sizeof(*look));
203  int i,j,n=0;
204
205  look->vi=info;
206  look->n=info->postlist[1];
207
208  /* we drop each position value in-between already decoded values,
209     and use linear interpolation to predict each new value past the
210     edges.  The positions are read in the order of the position
211     list... we precompute the bounding positions in the lookup.  Of
212     course, the neighbors can change (if a position is declined), but
213     this is an initial mapping */
214
215  for(i=0;i<info->partitions;i++)n+=info->class_dim[info->partitionclass[i]];
216  n+=2;
217  look->posts=n;
218
219  /* also store a sorted position index */
220  for(i=0;i<n;i++)sortpointer[i]=info->postlist+i;
221  qsort(sortpointer,n,sizeof(*sortpointer),icomp);
222
223  /* points from sort order back to range number */
224  for(i=0;i<n;i++)look->forward_index[i]=sortpointer[i]-info->postlist;
225  /* points from range order to sorted position */
226  for(i=0;i<n;i++)look->reverse_index[look->forward_index[i]]=i;
227  /* we actually need the post values too */
228  for(i=0;i<n;i++)look->sorted_index[i]=info->postlist[look->forward_index[i]];
229
230  /* quantize values to multiplier spec */
231  switch(info->mult){
232  case 1: /* 1024 -> 256 */
233    look->quant_q=256;
234    break;
235  case 2: /* 1024 -> 128 */
236    look->quant_q=128;
237    break;
238  case 3: /* 1024 -> 86 */
239    look->quant_q=86;
240    break;
241  case 4: /* 1024 -> 64 */
242    look->quant_q=64;
243    break;
244  }
245
246  /* discover our neighbors for decode where we don't use fit flags
247     (that would push the neighbors outward) */
248  for(i=0;i<n-2;i++){
249    int lo=0;
250    int hi=1;
251    int lx=0;
252    int hx=look->n;
253    int currentx=info->postlist[i+2];
254    for(j=0;j<i+2;j++){
255      int x=info->postlist[j];
256      if(x>lx && x<currentx){
257        lo=j;
258        lx=x;
259      }
260      if(x<hx && x>currentx){
261        hi=j;
262        hx=x;
263      }
264    }
265    look->loneighbor[i]=lo;
266    look->hineighbor[i]=hi;
267  }
268
269  return(look);
270}
271
272static int render_point(int x0,int x1,int y0,int y1,int x){
273  y0&=0x7fff; /* mask off flag */
274  y1&=0x7fff;
275
276  {
277    int dy=y1-y0;
278    int adx=x1-x0;
279    int ady=abs(dy);
280    int err=ady*(x-x0);
281
282    int off=err/adx;
283    if(dy<0)return(y0-off);
284    return(y0+off);
285  }
286}
287
288static int vorbis_dBquant(const float *x){
289  int i= *x*7.3142857f+1023.5f;
290  if(i>1023)return(1023);
291  if(i<0)return(0);
292  return i;
293}
294
295static const float FLOOR1_fromdB_LOOKUP[256]={
296  1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F,
297  1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F,
298  1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F,
299  2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F,
300  2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F,
301  3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F,
302  4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F,
303  6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F,
304  7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F,
305  1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F,
306  1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F,
307  1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F,
308  2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F,
309  2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F,
310  3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F,
311  4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F,
312  5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F,
313  7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F,
314  9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F,
315  1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F,
316  1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F,
317  2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F,
318  2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F,
319  3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F,
320  4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F,
321  5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F,
322  7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F,
323  9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F,
324  0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F,
325  0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F,
326  0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F,
327  0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F,
328  0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F,
329  0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F,
330  0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F,
331  0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F,
332  0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F,
333  0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F,
334  0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F,
335  0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F,
336  0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F,
337  0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F,
338  0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F,
339  0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F,
340  0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F,
341  0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F,
342  0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F,
343  0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F,
344  0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F,
345  0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F,
346  0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F,
347  0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F,
348  0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F,
349  0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F,
350  0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F,
351  0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F,
352  0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F,
353  0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F,
354  0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F,
355  0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F,
356  0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F,
357  0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F,
358  0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F,
359  0.82788260F, 0.88168307F, 0.9389798F, 1.F,
360};
361
362static void render_line(int n, int x0,int x1,int y0,int y1,float *d){
363  int dy=y1-y0;
364  int adx=x1-x0;
365  int ady=abs(dy);
366  int base=dy/adx;
367  int sy=(dy<0?base-1:base+1);
368  int x=x0;
369  int y=y0;
370  int err=0;
371
372  ady-=abs(base*adx);
373
374  if(n>x1)n=x1;
375
376  if(x<n)
377    d[x]*=FLOOR1_fromdB_LOOKUP[y];
378
379  while(++x<n){
380    err=err+ady;
381    if(err>=adx){
382      err-=adx;
383      y+=sy;
384    }else{
385      y+=base;
386    }
387    d[x]*=FLOOR1_fromdB_LOOKUP[y];
388  }
389}
390
391static void render_line0(int n, int x0,int x1,int y0,int y1,int *d){
392  int dy=y1-y0;
393  int adx=x1-x0;
394  int ady=abs(dy);
395  int base=dy/adx;
396  int sy=(dy<0?base-1:base+1);
397  int x=x0;
398  int y=y0;
399  int err=0;
400
401  ady-=abs(base*adx);
402
403  if(n>x1)n=x1;
404
405  if(x<n)
406    d[x]=y;
407
408  while(++x<n){
409    err=err+ady;
410    if(err>=adx){
411      err-=adx;
412      y+=sy;
413    }else{
414      y+=base;
415    }
416    d[x]=y;
417  }
418}
419
420/* the floor has already been filtered to only include relevant sections */
421static int accumulate_fit(const float *flr,const float *mdct,
422                          int x0, int x1,lsfit_acc *a,
423                          int n,vorbis_info_floor1 *info){
424  long i;
425
426  int xa=0,ya=0,x2a=0,y2a=0,xya=0,na=0, xb=0,yb=0,x2b=0,y2b=0,xyb=0,nb=0;
427
428  memset(a,0,sizeof(*a));
429  a->x0=x0;
430  a->x1=x1;
431  if(x1>=n)x1=n-1;
432
433  for(i=x0;i<=x1;i++){
434    int quantized=vorbis_dBquant(flr+i);
435    if(quantized){
436      if(mdct[i]+info->twofitatten>=flr[i]){
437        xa  += i;
438        ya  += quantized;
439        x2a += i*i;
440        y2a += quantized*quantized;
441        xya += i*quantized;
442        na++;
443      }else{
444        xb  += i;
445        yb  += quantized;
446        x2b += i*i;
447        y2b += quantized*quantized;
448        xyb += i*quantized;
449        nb++;
450      }
451    }
452  }
453
454  a->xa=xa;
455  a->ya=ya;
456  a->x2a=x2a;
457  a->y2a=y2a;
458  a->xya=xya;
459  a->an=na;
460
461  a->xb=xb;
462  a->yb=yb;
463  a->x2b=x2b;
464  a->y2b=y2b;
465  a->xyb=xyb;
466  a->bn=nb;
467
468  return(na);
469}
470
471static int fit_line(lsfit_acc *a,int fits,int *y0,int *y1,
472                    vorbis_info_floor1 *info){
473  double xb=0,yb=0,x2b=0,y2b=0,xyb=0,bn=0;
474  int i;
475  int x0=a[0].x0;
476  int x1=a[fits-1].x1;
477
478  for(i=0;i<fits;i++){
479    double weight = (a[i].bn+a[i].an)*info->twofitweight/(a[i].an+1)+1.;
480
481    xb+=a[i].xb + a[i].xa * weight;
482    yb+=a[i].yb + a[i].ya * weight;
483    x2b+=a[i].x2b + a[i].x2a * weight;
484    y2b+=a[i].y2b + a[i].y2a * weight;
485    xyb+=a[i].xyb + a[i].xya * weight;
486    bn+=a[i].bn + a[i].an * weight;
487  }
488
489  if(*y0>=0){
490    xb+=   x0;
491    yb+=  *y0;
492    x2b+=  x0 *  x0;
493    y2b+= *y0 * *y0;
494    xyb+= *y0 *  x0;
495    bn++;
496  }
497
498  if(*y1>=0){
499    xb+=   x1;
500    yb+=  *y1;
501    x2b+=  x1 *  x1;
502    y2b+= *y1 * *y1;
503    xyb+= *y1 *  x1;
504    bn++;
505  }
506
507  {
508    double denom=(bn*x2b-xb*xb);
509
510    if(denom>0.){
511      double a=(yb*x2b-xyb*xb)/denom;
512      double b=(bn*xyb-xb*yb)/denom;
513      *y0=rint(a+b*x0);
514      *y1=rint(a+b*x1);
515
516      /* limit to our range! */
517      if(*y0>1023)*y0=1023;
518      if(*y1>1023)*y1=1023;
519      if(*y0<0)*y0=0;
520      if(*y1<0)*y1=0;
521
522      return 0;
523    }else{
524      *y0=0;
525      *y1=0;
526      return 1;
527    }
528  }
529}
530
531static int inspect_error(int x0,int x1,int y0,int y1,const float *mask,
532                         const float *mdct,
533                         vorbis_info_floor1 *info){
534  int dy=y1-y0;
535  int adx=x1-x0;
536  int ady=abs(dy);
537  int base=dy/adx;
538  int sy=(dy<0?base-1:base+1);
539  int x=x0;
540  int y=y0;
541  int err=0;
542  int val=vorbis_dBquant(mask+x);
543  int mse=0;
544  int n=0;
545
546  ady-=abs(base*adx);
547
548  mse=(y-val);
549  mse*=mse;
550  n++;
551  if(mdct[x]+info->twofitatten>=mask[x]){
552    if(y+info->maxover<val)return(1);
553    if(y-info->maxunder>val)return(1);
554  }
555
556  while(++x<x1){
557    err=err+ady;
558    if(err>=adx){
559      err-=adx;
560      y+=sy;
561    }else{
562      y+=base;
563    }
564
565    val=vorbis_dBquant(mask+x);
566    mse+=((y-val)*(y-val));
567    n++;
568    if(mdct[x]+info->twofitatten>=mask[x]){
569      if(val){
570        if(y+info->maxover<val)return(1);
571        if(y-info->maxunder>val)return(1);
572      }
573    }
574  }
575
576  if(info->maxover*info->maxover/n>info->maxerr)return(0);
577  if(info->maxunder*info->maxunder/n>info->maxerr)return(0);
578  if(mse/n>info->maxerr)return(1);
579  return(0);
580}
581
582static int post_Y(int *A,int *B,int pos){
583  if(A[pos]<0)
584    return B[pos];
585  if(B[pos]<0)
586    return A[pos];
587
588  return (A[pos]+B[pos])>>1;
589}
590
591int *floor1_fit(vorbis_block *vb,vorbis_look_floor1 *look,
592                          const float *logmdct,   /* in */
593                          const float *logmask){
594  long i,j;
595  vorbis_info_floor1 *info=look->vi;
596  long n=look->n;
597  long posts=look->posts;
598  long nonzero=0;
599  lsfit_acc fits[VIF_POSIT+1];
600  int fit_valueA[VIF_POSIT+2]; /* index by range list position */
601  int fit_valueB[VIF_POSIT+2]; /* index by range list position */
602
603  int loneighbor[VIF_POSIT+2]; /* sorted index of range list position (+2) */
604  int hineighbor[VIF_POSIT+2];
605  int *output=NULL;
606  int memo[VIF_POSIT+2];
607
608  for(i=0;i<posts;i++)fit_valueA[i]=-200; /* mark all unused */
609  for(i=0;i<posts;i++)fit_valueB[i]=-200; /* mark all unused */
610  for(i=0;i<posts;i++)loneighbor[i]=0; /* 0 for the implicit 0 post */
611  for(i=0;i<posts;i++)hineighbor[i]=1; /* 1 for the implicit post at n */
612  for(i=0;i<posts;i++)memo[i]=-1;      /* no neighbor yet */
613
614  /* quantize the relevant floor points and collect them into line fit
615     structures (one per minimal division) at the same time */
616  if(posts==0){
617    nonzero+=accumulate_fit(logmask,logmdct,0,n,fits,n,info);
618  }else{
619    for(i=0;i<posts-1;i++)
620      nonzero+=accumulate_fit(logmask,logmdct,look->sorted_index[i],
621                              look->sorted_index[i+1],fits+i,
622                              n,info);
623  }
624
625  if(nonzero){
626    /* start by fitting the implicit base case.... */
627    int y0=-200;
628    int y1=-200;
629    fit_line(fits,posts-1,&y0,&y1,info);
630
631    fit_valueA[0]=y0;
632    fit_valueB[0]=y0;
633    fit_valueB[1]=y1;
634    fit_valueA[1]=y1;
635
636    /* Non degenerate case */
637    /* start progressive splitting.  This is a greedy, non-optimal
638       algorithm, but simple and close enough to the best
639       answer. */
640    for(i=2;i<posts;i++){
641      int sortpos=look->reverse_index[i];
642      int ln=loneighbor[sortpos];
643      int hn=hineighbor[sortpos];
644
645      /* eliminate repeat searches of a particular range with a memo */
646      if(memo[ln]!=hn){
647        /* haven't performed this error search yet */
648        int lsortpos=look->reverse_index[ln];
649        int hsortpos=look->reverse_index[hn];
650        memo[ln]=hn;
651
652        {
653          /* A note: we want to bound/minimize *local*, not global, error */
654          int lx=info->postlist[ln];
655          int hx=info->postlist[hn];
656          int ly=post_Y(fit_valueA,fit_valueB,ln);
657          int hy=post_Y(fit_valueA,fit_valueB,hn);
658
659          if(ly==-1 || hy==-1){
660            exit(1);
661          }
662
663          if(inspect_error(lx,hx,ly,hy,logmask,logmdct,info)){
664            /* outside error bounds/begin search area.  Split it. */
665            int ly0=-200;
666            int ly1=-200;
667            int hy0=-200;
668            int hy1=-200;
669            int ret0=fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1,info);
670            int ret1=fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1,info);
671
672            if(ret0){
673              ly0=ly;
674              ly1=hy0;
675            }
676            if(ret1){
677              hy0=ly1;
678              hy1=hy;
679            }
680
681            if(ret0 && ret1){
682              fit_valueA[i]=-200;
683              fit_valueB[i]=-200;
684            }else{
685              /* store new edge values */
686              fit_valueB[ln]=ly0;
687              if(ln==0)fit_valueA[ln]=ly0;
688              fit_valueA[i]=ly1;
689              fit_valueB[i]=hy0;
690              fit_valueA[hn]=hy1;
691              if(hn==1)fit_valueB[hn]=hy1;
692
693              if(ly1>=0 || hy0>=0){
694                /* store new neighbor values */
695                for(j=sortpos-1;j>=0;j--)
696                  if(hineighbor[j]==hn)
697                    hineighbor[j]=i;
698                  else
699                    break;
700                for(j=sortpos+1;j<posts;j++)
701                  if(loneighbor[j]==ln)
702                    loneighbor[j]=i;
703                  else
704                    break;
705              }
706            }
707          }else{
708            fit_valueA[i]=-200;
709            fit_valueB[i]=-200;
710          }
711        }
712      }
713    }
714
715    output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
716
717    output[0]=post_Y(fit_valueA,fit_valueB,0);
718    output[1]=post_Y(fit_valueA,fit_valueB,1);
719
720    /* fill in posts marked as not using a fit; we will zero
721       back out to 'unused' when encoding them so long as curve
722       interpolation doesn't force them into use */
723    for(i=2;i<posts;i++){
724      int ln=look->loneighbor[i-2];
725      int hn=look->hineighbor[i-2];
726      int x0=info->postlist[ln];
727      int x1=info->postlist[hn];
728      int y0=output[ln];
729      int y1=output[hn];
730
731      int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
732      int vx=post_Y(fit_valueA,fit_valueB,i);
733
734      if(vx>=0 && predicted!=vx){
735        output[i]=vx;
736      }else{
737        output[i]= predicted|0x8000;
738      }
739    }
740  }
741
742  return(output);
743
744}
745
746int *floor1_interpolate_fit(vorbis_block *vb,vorbis_look_floor1 *look,
747                          int *A,int *B,
748                          int del){
749
750  long i;
751  long posts=look->posts;
752  int *output=NULL;
753
754  if(A && B){
755    output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
756
757    /* overly simpleminded--- look again post 1.2 */
758    for(i=0;i<posts;i++){
759      output[i]=((65536-del)*(A[i]&0x7fff)+del*(B[i]&0x7fff)+32768)>>16;
760      if(A[i]&0x8000 && B[i]&0x8000)output[i]|=0x8000;
761    }
762  }
763
764  return(output);
765}
766
767
768int floor1_encode(oggpack_buffer *opb,vorbis_block *vb,
769                  vorbis_look_floor1 *look,
770                  int *post,int *ilogmask){
771
772  long i,j;
773  vorbis_info_floor1 *info=look->vi;
774  long posts=look->posts;
775  codec_setup_info *ci=vb->vd->vi->codec_setup;
776  int out[VIF_POSIT+2];
777  static_codebook **sbooks=ci->book_param;
778  codebook *books=ci->fullbooks;
779
780  /* quantize values to multiplier spec */
781  if(post){
782    for(i=0;i<posts;i++){
783      int val=post[i]&0x7fff;
784      switch(info->mult){
785      case 1: /* 1024 -> 256 */
786        val>>=2;
787        break;
788      case 2: /* 1024 -> 128 */
789        val>>=3;
790        break;
791      case 3: /* 1024 -> 86 */
792        val/=12;
793        break;
794      case 4: /* 1024 -> 64 */
795        val>>=4;
796        break;
797      }
798      post[i]=val | (post[i]&0x8000);
799    }
800
801    out[0]=post[0];
802    out[1]=post[1];
803
804    /* find prediction values for each post and subtract them */
805    for(i=2;i<posts;i++){
806      int ln=look->loneighbor[i-2];
807      int hn=look->hineighbor[i-2];
808      int x0=info->postlist[ln];
809      int x1=info->postlist[hn];
810      int y0=post[ln];
811      int y1=post[hn];
812
813      int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
814
815      if((post[i]&0x8000) || (predicted==post[i])){
816        post[i]=predicted|0x8000; /* in case there was roundoff jitter
817                                     in interpolation */
818        out[i]=0;
819      }else{
820        int headroom=(look->quant_q-predicted<predicted?
821                      look->quant_q-predicted:predicted);
822
823        int val=post[i]-predicted;
824
825        /* at this point the 'deviation' value is in the range +/- max
826           range, but the real, unique range can always be mapped to
827           only [0-maxrange).  So we want to wrap the deviation into
828           this limited range, but do it in the way that least screws
829           an essentially gaussian probability distribution. */
830
831        if(val<0)
832          if(val<-headroom)
833            val=headroom-val-1;
834          else
835            val=-1-(val<<1);
836        else
837          if(val>=headroom)
838            val= val+headroom;
839          else
840            val<<=1;
841
842        out[i]=val;
843        post[ln]&=0x7fff;
844        post[hn]&=0x7fff;
845      }
846    }
847
848    /* we have everything we need. pack it out */
849    /* mark nontrivial floor */
850    oggpack_write(opb,1,1);
851
852    /* beginning/end post */
853    look->frames++;
854    look->postbits+=ilog(look->quant_q-1)*2;
855    oggpack_write(opb,out[0],ilog(look->quant_q-1));
856    oggpack_write(opb,out[1],ilog(look->quant_q-1));
857
858
859    /* partition by partition */
860    for(i=0,j=2;i<info->partitions;i++){
861      int class=info->partitionclass[i];
862      int cdim=info->class_dim[class];
863      int csubbits=info->class_subs[class];
864      int csub=1<<csubbits;
865      int bookas[8]={0,0,0,0,0,0,0,0};
866      int cval=0;
867      int cshift=0;
868      int k,l;
869
870      /* generate the partition's first stage cascade value */
871      if(csubbits){
872        int maxval[8];
873        for(k=0;k<csub;k++){
874          int booknum=info->class_subbook[class][k];
875          if(booknum<0){
876            maxval[k]=1;
877          }else{
878            maxval[k]=sbooks[info->class_subbook[class][k]]->entries;
879          }
880        }
881        for(k=0;k<cdim;k++){
882          for(l=0;l<csub;l++){
883            int val=out[j+k];
884            if(val<maxval[l]){
885              bookas[k]=l;
886              break;
887            }
888          }
889          cval|= bookas[k]<<cshift;
890          cshift+=csubbits;
891        }
892        /* write it */
893        look->phrasebits+=
894          vorbis_book_encode(books+info->class_book[class],cval,opb);
895
896#ifdef TRAIN_FLOOR1
897        {
898          FILE *of;
899          char buffer[80];
900          sprintf(buffer,"line_%dx%ld_class%d.vqd",
901                  vb->pcmend/2,posts-2,class);
902          of=fopen(buffer,"a");
903          fprintf(of,"%d\n",cval);
904          fclose(of);
905        }
906#endif
907      }
908
909      /* write post values */
910      for(k=0;k<cdim;k++){
911        int book=info->class_subbook[class][bookas[k]];
912        if(book>=0){
913          /* hack to allow training with 'bad' books */
914          if(out[j+k]<(books+book)->entries)
915            look->postbits+=vorbis_book_encode(books+book,
916                                               out[j+k],opb);
917          /*else
918            fprintf(stderr,"+!");*/
919
920#ifdef TRAIN_FLOOR1
921          {
922            FILE *of;
923            char buffer[80];
924            sprintf(buffer,"line_%dx%ld_%dsub%d.vqd",
925                    vb->pcmend/2,posts-2,class,bookas[k]);
926            of=fopen(buffer,"a");
927            fprintf(of,"%d\n",out[j+k]);
928            fclose(of);
929          }
930#endif
931        }
932      }
933      j+=cdim;
934    }
935
936    {
937      /* generate quantized floor equivalent to what we'd unpack in decode */
938      /* render the lines */
939      int hx=0;
940      int lx=0;
941      int ly=post[0]*info->mult;
942      int n=ci->blocksizes[vb->W]/2;
943
944      for(j=1;j<look->posts;j++){
945        int current=look->forward_index[j];
946        int hy=post[current]&0x7fff;
947        if(hy==post[current]){
948
949          hy*=info->mult;
950          hx=info->postlist[current];
951
952          render_line0(n,lx,hx,ly,hy,ilogmask);
953
954          lx=hx;
955          ly=hy;
956        }
957      }
958      for(j=hx;j<vb->pcmend/2;j++)ilogmask[j]=ly; /* be certain */
959      return(1);
960    }
961  }else{
962    oggpack_write(opb,0,1);
963    memset(ilogmask,0,vb->pcmend/2*sizeof(*ilogmask));
964    return(0);
965  }
966}
967
968static void *floor1_inverse1(vorbis_block *vb,vorbis_look_floor *in){
969  vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
970  vorbis_info_floor1 *info=look->vi;
971  codec_setup_info   *ci=vb->vd->vi->codec_setup;
972
973  int i,j,k;
974  codebook *books=ci->fullbooks;
975
976  /* unpack wrapped/predicted values from stream */
977  if(oggpack_read(&vb->opb,1)==1){
978    int *fit_value=_vorbis_block_alloc(vb,(look->posts)*sizeof(*fit_value));
979
980    fit_value[0]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
981    fit_value[1]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
982
983    /* partition by partition */
984    for(i=0,j=2;i<info->partitions;i++){
985      int class=info->partitionclass[i];
986      int cdim=info->class_dim[class];
987      int csubbits=info->class_subs[class];
988      int csub=1<<csubbits;
989      int cval=0;
990
991      /* decode the partition's first stage cascade value */
992      if(csubbits){
993        cval=vorbis_book_decode(books+info->class_book[class],&vb->opb);
994
995        if(cval==-1)goto eop;
996      }
997
998      for(k=0;k<cdim;k++){
999        int book=info->class_subbook[class][cval&(csub-1)];
1000        cval>>=csubbits;
1001        if(book>=0){
1002          if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1)
1003            goto eop;
1004        }else{
1005          fit_value[j+k]=0;
1006        }
1007      }
1008      j+=cdim;
1009    }
1010
1011    /* unwrap positive values and reconsitute via linear interpolation */
1012    for(i=2;i<look->posts;i++){
1013      int predicted=render_point(info->postlist[look->loneighbor[i-2]],
1014                                 info->postlist[look->hineighbor[i-2]],
1015                                 fit_value[look->loneighbor[i-2]],
1016                                 fit_value[look->hineighbor[i-2]],
1017                                 info->postlist[i]);
1018      int hiroom=look->quant_q-predicted;
1019      int loroom=predicted;
1020      int room=(hiroom<loroom?hiroom:loroom)<<1;
1021      int val=fit_value[i];
1022
1023      if(val){
1024        if(val>=room){
1025          if(hiroom>loroom){
1026            val = val-loroom;
1027          }else{
1028            val = -1-(val-hiroom);
1029          }
1030        }else{
1031          if(val&1){
1032            val= -((val+1)>>1);
1033          }else{
1034            val>>=1;
1035          }
1036        }
1037
1038        fit_value[i]=val+predicted;
1039        fit_value[look->loneighbor[i-2]]&=0x7fff;
1040        fit_value[look->hineighbor[i-2]]&=0x7fff;
1041
1042      }else{
1043        fit_value[i]=predicted|0x8000;
1044      }
1045
1046    }
1047
1048    return(fit_value);
1049  }
1050 eop:
1051  return(NULL);
1052}
1053
1054static int floor1_inverse2(vorbis_block *vb,vorbis_look_floor *in,void *memo,
1055                          float *out){
1056  vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1057  vorbis_info_floor1 *info=look->vi;
1058
1059  codec_setup_info   *ci=vb->vd->vi->codec_setup;
1060  int                  n=ci->blocksizes[vb->W]/2;
1061  int j;
1062
1063  if(memo){
1064    /* render the lines */
1065    int *fit_value=(int *)memo;
1066    int hx=0;
1067    int lx=0;
1068    int ly=fit_value[0]*info->mult;
1069    for(j=1;j<look->posts;j++){
1070      int current=look->forward_index[j];
1071      int hy=fit_value[current]&0x7fff;
1072      if(hy==fit_value[current]){
1073
1074        hy*=info->mult;
1075        hx=info->postlist[current];
1076
1077        render_line(n,lx,hx,ly,hy,out);
1078
1079        lx=hx;
1080        ly=hy;
1081      }
1082    }
1083    for(j=hx;j<n;j++)out[j]*=FLOOR1_fromdB_LOOKUP[ly]; /* be certain */
1084    return(1);
1085  }
1086  memset(out,0,sizeof(*out)*n);
1087  return(0);
1088}
1089
1090/* export hooks */
1091const vorbis_func_floor floor1_exportbundle={
1092  &floor1_pack,&floor1_unpack,&floor1_look,&floor1_free_info,
1093  &floor1_free_look,&floor1_inverse1,&floor1_inverse2
1094};
1095