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-2009             *
   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];
 859         for(k=0;k<csub;k++){
 860           int booknum=info->class_subbook[class][k];
 861           if(booknum<0){
 862             maxval[k]=1;
 863           }else{
 864             maxval[k]=sbooks[info->class_subbook[class][k]]->entries;
 865           }
 866         }
 867         for(k=0;k<cdim;k++){
 868           for(l=0;l<csub;l++){
 869             int val=out[j+k];
 870             if(val<maxval[l]){
 871               bookas[k]=l;
 872               break;
 873             }
 874           }
 875           cval|= bookas[k]<<cshift;
 876           cshift+=csubbits;
 877         }
 878         /* write it */
 879         look->phrasebits+=
 880           vorbis_book_encode(books+info->class_book[class],cval,opb);
 881
 882 #ifdef TRAIN_FLOOR1
 883         {
 884           FILE *of;
 885           char buffer[80];
 886           sprintf(buffer,"line_%dx%ld_class%d.vqd",
 887                   vb->pcmend/2,posts-2,class);
 888           of=fopen(buffer,"a");
 889           fprintf(of,"%d\n",cval);
 890           fclose(of);
 891         }
 892 #endif
 893       }
 894
 895       /* write post values */
 896       for(k=0;k<cdim;k++){
 897         int book=info->class_subbook[class][bookas[k]];
 898         if(book>=0){
 899           /* hack to allow training with 'bad' books */
 900           if(out[j+k]<(books+book)->entries)
 901             look->postbits+=vorbis_book_encode(books+book,
 902                                                out[j+k],opb);
 903           /*else
 904             fprintf(stderr,"+!");*/
 905
 906 #ifdef TRAIN_FLOOR1
 907           {
 908             FILE *of;
 909             char buffer[80];
 910             sprintf(buffer,"line_%dx%ld_%dsub%d.vqd",
 911                     vb->pcmend/2,posts-2,class,bookas[k]);
 912             of=fopen(buffer,"a");
 913             fprintf(of,"%d\n",out[j+k]);
 914             fclose(of);
 915           }
 916 #endif
 917         }
 918       }
 919       j+=cdim;
 920     }
 921
 922     {
 923       /* generate quantized floor equivalent to what we'd unpack in decode */
 924       /* render the lines */
 925       int hx=0;
 926       int lx=0;
 927       int ly=post[0]*info->mult;
 928       int n=ci->blocksizes[vb->W]/2;
 929
 930       for(j=1;j<look->posts;j++){
 931         int current=look->forward_index[j];
 932         int hy=post[current]&0x7fff;
 933         if(hy==post[current]){
 934
 935           hy*=info->mult;
 936           hx=info->postlist[current];
 937
 938           render_line0(n,lx,hx,ly,hy,ilogmask);
 939
 940           lx=hx;
 941           ly=hy;
 942         }
 943       }
 944       for(j=hx;j<vb->pcmend/2;j++)ilogmask[j]=ly; /* be certain */
 945       return(1);
 946     }
 947   }else{
 948     oggpack_write(opb,0,1);
 949     memset(ilogmask,0,vb->pcmend/2*sizeof(*ilogmask));
 950     return(0);
 951   }
 952 }
 953
 954 static void *floor1_inverse1(vorbis_block *vb,vorbis_look_floor *in){
 955   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
 956   vorbis_info_floor1 *info=look->vi;
 957   codec_setup_info   *ci=vb->vd->vi->codec_setup;
 958
 959   int i,j,k;
 960   codebook *books=ci->fullbooks;
 961
 962   /* unpack wrapped/predicted values from stream */
 963   if(oggpack_read(&vb->opb,1)==1){
 964     int *fit_value=_vorbis_block_alloc(vb,(look->posts)*sizeof(*fit_value));
 965
 966     fit_value[0]=oggpack_read(&vb->opb,ov_ilog(look->quant_q-1));
 967     fit_value[1]=oggpack_read(&vb->opb,ov_ilog(look->quant_q-1));
 968
 969     /* partition by partition */
 970     for(i=0,j=2;i<info->partitions;i++){
 971       int class=info->partitionclass[i];
 972       int cdim=info->class_dim[class];
 973       int csubbits=info->class_subs[class];
 974       int csub=1<<csubbits;
 975       int cval=0;
 976
 977       /* decode the partition's first stage cascade value */
 978       if(csubbits){
 979         cval=vorbis_book_decode(books+info->class_book[class],&vb->opb);
 980
 981         if(cval==-1)goto eop;
 982       }
 983
 984       for(k=0;k<cdim;k++){
 985         int book=info->class_subbook[class][cval&(csub-1)];
 986         cval>>=csubbits;
 987         if(book>=0){
 988           if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1)
 989             goto eop;
 990         }else{
 991           fit_value[j+k]=0;
 992         }
 993       }
 994       j+=cdim;
 995     }
 996
 997     /* unwrap positive values and reconsitute via linear interpolation */
 998     for(i=2;i<look->posts;i++){
 999       int predicted=render_point(info->postlist[look->loneighbor[i-2]],
1000                                  info->postlist[look->hineighbor[i-2]],
1001                                  fit_value[look->loneighbor[i-2]],
1002                                  fit_value[look->hineighbor[i-2]],
1003                                  info->postlist[i]);
1004       int hiroom=look->quant_q-predicted;
1005       int loroom=predicted;
1006       int room=(hiroom<loroom?hiroom:loroom)<<1;
1007       int val=fit_value[i];
1008
1009       if(val){
1010         if(val>=room){
1011           if(hiroom>loroom){
1012             val = val-loroom;
1013           }else{
1014             val = -1-(val-hiroom);
1015           }
1016         }else{
1017           if(val&1){
1018             val= -((val+1)>>1);
1019           }else{
1020             val>>=1;
1021           }
1022         }
1023
1024         fit_value[i]=(val+predicted)&0x7fff;
1025         fit_value[look->loneighbor[i-2]]&=0x7fff;
1026         fit_value[look->hineighbor[i-2]]&=0x7fff;
1027
1028       }else{
1029         fit_value[i]=predicted|0x8000;
1030       }
1031
1032     }
1033
1034     return(fit_value);
1035   }
1036  eop:
1037   return(NULL);
1038 }
1039
1040 static int floor1_inverse2(vorbis_block *vb,vorbis_look_floor *in,void *memo,
1041                           float *out){
1042   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1043   vorbis_info_floor1 *info=look->vi;
1044
1045   codec_setup_info   *ci=vb->vd->vi->codec_setup;
1046   int                  n=ci->blocksizes[vb->W]/2;
1047   int j;
1048
1049   if(memo){
1050     /* render the lines */
1051     int *fit_value=(int *)memo;
1052     int hx=0;
1053     int lx=0;
1054     int ly=fit_value[0]*info->mult;
1055     /* guard lookup against out-of-range values */
1056     ly=(ly<0?0:ly>255?255:ly);
1057
1058     for(j=1;j<look->posts;j++){
1059       int current=look->forward_index[j];
1060       int hy=fit_value[current]&0x7fff;
1061       if(hy==fit_value[current]){
1062
1063         hx=info->postlist[current];
1064         hy*=info->mult;
1065         /* guard lookup against out-of-range values */
1066         hy=(hy<0?0:hy>255?255:hy);
1067
1068         render_line(n,lx,hx,ly,hy,out);
1069
1070         lx=hx;
1071         ly=hy;
1072       }
1073     }
1074     for(j=hx;j<n;j++)out[j]*=FLOOR1_fromdB_LOOKUP[ly]; /* be certain */
1075     return(1);
1076   }
1077   memset(out,0,sizeof(*out)*n);
1078   return(0);
1079 }
1080
1081 /* export hooks */
1082 const vorbis_func_floor floor1_exportbundle={
1083   &floor1_pack,&floor1_unpack,&floor1_look,&floor1_free_info,
1084   &floor1_free_look,&floor1_inverse1,&floor1_inverse2
1085 };