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