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