lib/floor1.c

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