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
 755   /* quantize values to multiplier spec */
 756   if(post){
 757     for(i=0;i<posts;i++){
 758       int val=post[i]&0x7fff;
 759       switch(info->mult){
 760       case 1: /* 1024 -> 256 */
 761         val>>=2;
 762         break;
 763       case 2: /* 1024 -> 128 */
 764         val>>=3;
 765         break;
 766       case 3: /* 1024 -> 86 */
 767         val/=12;
 768         break;
 769       case 4: /* 1024 -> 64 */
 770         val>>=4;
 771         break;
 772       }
 773       post[i]=val | (post[i]&0x8000);
 774     }
 775
 776     out[0]=post[0];
 777     out[1]=post[1];
 778
 779     /* find prediction values for each post and subtract them */
 780     for(i=2;i<posts;i++){
 781       int ln=look->loneighbor[i-2];
 782       int hn=look->hineighbor[i-2];
 783       int x0=info->postlist[ln];
 784       int x1=info->postlist[hn];
 785       int y0=post[ln];
 786       int y1=post[hn];
 787
 788       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
 789
 790       if((post[i]&0x8000) || (predicted==post[i])){
 791         post[i]=predicted|0x8000; /* in case there was roundoff jitter
 792                                      in interpolation */
 793         out[i]=0;
 794       }else{
 795         int headroom=(look->quant_q-predicted<predicted?
 796                       look->quant_q-predicted:predicted);
 797
 798         int val=post[i]-predicted;
 799
 800         /* at this point the 'deviation' value is in the range +/- max
 801            range, but the real, unique range can always be mapped to
 802            only [0-maxrange).  So we want to wrap the deviation into
 803            this limited range, but do it in the way that least screws
 804            an essentially gaussian probability distribution. */
 805
 806         if(val<0)
 807           if(val<-headroom)
 808             val=headroom-val-1;
 809           else
 810             val=-1-(val<<1);
 811         else
 812           if(val>=headroom)
 813             val= val+headroom;
 814           else
 815             val<<=1;
 816
 817         out[i]=val;
 818         post[ln]&=0x7fff;
 819         post[hn]&=0x7fff;
 820       }
 821     }
 822
 823     /* we have everything we need. pack it out */
 824     /* mark nontrivial floor */
 825     oggpack_write(opb,1,1);
 826
 827     /* beginning/end post */
 828     look->frames++;
 829     look->postbits+=ilog(look->quant_q-1)*2;
 830     oggpack_write(opb,out[0],ilog(look->quant_q-1));
 831     oggpack_write(opb,out[1],ilog(look->quant_q-1));
 832
 833
 834     /* partition by partition */
 835     for(i=0,j=2;i<info->partitions;i++){
 836       int class=info->partitionclass[i];
 837       int cdim=info->class_dim[class];
 838       int csubbits=info->class_subs[class];
 839       int csub=1<<csubbits;
 840       int bookas[8]={0,0,0,0,0,0,0,0};
 841       int cval=0;
 842       int cshift=0;
 843       int k,l;
 844
 845       /* generate the partition's first stage cascade value */
 846       if(csubbits){
 847         int maxval[8];
 848         for(k=0;k<csub;k++){
 849           int booknum=info->class_subbook[class][k];
 850           if(booknum<0){
 851             maxval[k]=1;
 852           }else{
 853             maxval[k]=sbooks[info->class_subbook[class][k]]->entries;
 854           }
 855         }
 856         for(k=0;k<cdim;k++){
 857           for(l=0;l<csub;l++){
 858             int val=out[j+k];
 859             if(val<maxval[l]){
 860               bookas[k]=l;
 861               break;
 862             }
 863           }
 864           cval|= bookas[k]<<cshift;
 865           cshift+=csubbits;
 866         }
 867         /* write it */
 868         look->phrasebits+=
 869           vorbis_book_encode(books+info->class_book[class],cval,opb);
 870
 871 #ifdef TRAIN_FLOOR1
 872         {
 873           FILE *of;
 874           char buffer[80];
 875           sprintf(buffer,"line_%dx%ld_class%d.vqd",
 876                   vb->pcmend/2,posts-2,class);
 877           of=fopen(buffer,"a");
 878           fprintf(of,"%d\n",cval);
 879           fclose(of);
 880         }
 881 #endif
 882       }
 883
 884       /* write post values */
 885       for(k=0;k<cdim;k++){
 886         int book=info->class_subbook[class][bookas[k]];
 887         if(book>=0){
 888           /* hack to allow training with 'bad' books */
 889           if(out[j+k]<(books+book)->entries)
 890             look->postbits+=vorbis_book_encode(books+book,
 891                                                out[j+k],opb);
 892           /*else
 893             fprintf(stderr,"+!");*/
 894
 895 #ifdef TRAIN_FLOOR1
 896           {
 897             FILE *of;
 898             char buffer[80];
 899             sprintf(buffer,"line_%dx%ld_%dsub%d.vqd",
 900                     vb->pcmend/2,posts-2,class,bookas[k]);
 901             of=fopen(buffer,"a");
 902             fprintf(of,"%d\n",out[j+k]);
 903             fclose(of);
 904           }
 905 #endif
 906         }
 907       }
 908       j+=cdim;
 909     }
 910
 911     {
 912       /* generate quantized floor equivalent to what we'd unpack in decode */
 913       /* render the lines */
 914       int hx=0;
 915       int lx=0;
 916       int ly=post[0]*info->mult;
 917       for(j=1;j<look->posts;j++){
 918         int current=look->forward_index[j];
 919         int hy=post[current]&0x7fff;
 920         if(hy==post[current]){
 921
 922           hy*=info->mult;
 923           hx=info->postlist[current];
 924
 925           render_line0(lx,hx,ly,hy,ilogmask);
 926
 927           lx=hx;
 928           ly=hy;
 929         }
 930       }
 931       for(j=hx;j<vb->pcmend/2;j++)ilogmask[j]=ly; /* be certain */
 932       return(1);
 933     }
 934   }else{
 935     oggpack_write(opb,0,1);
 936     memset(ilogmask,0,vb->pcmend/2*sizeof(*ilogmask));
 937     return(0);
 938   }
 939 }
 940
 941 static void *floor1_inverse1(vorbis_block *vb,vorbis_look_floor *in){
 942   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
 943   vorbis_info_floor1 *info=look->vi;
 944   codec_setup_info   *ci=vb->vd->vi->codec_setup;
 945
 946   int i,j,k;
 947   codebook *books=ci->fullbooks;
 948
 949   /* unpack wrapped/predicted values from stream */
 950   if(oggpack_read(&vb->opb,1)==1){
 951     int *fit_value=_vorbis_block_alloc(vb,(look->posts)*sizeof(*fit_value));
 952
 953     fit_value[0]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
 954     fit_value[1]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
 955
 956     /* partition by partition */
 957     for(i=0,j=2;i<info->partitions;i++){
 958       int class=info->partitionclass[i];
 959       int cdim=info->class_dim[class];
 960       int csubbits=info->class_subs[class];
 961       int csub=1<<csubbits;
 962       int cval=0;
 963
 964       /* decode the partition's first stage cascade value */
 965       if(csubbits){
 966         cval=vorbis_book_decode(books+info->class_book[class],&vb->opb);
 967
 968         if(cval==-1)goto eop;
 969       }
 970
 971       for(k=0;k<cdim;k++){
 972         int book=info->class_subbook[class][cval&(csub-1)];
 973         cval>>=csubbits;
 974         if(book>=0){
 975           if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1)
 976             goto eop;
 977         }else{
 978           fit_value[j+k]=0;
 979         }
 980       }
 981       j+=cdim;
 982     }
 983
 984     /* unwrap positive values and reconsitute via linear interpolation */
 985     for(i=2;i<look->posts;i++){
 986       int predicted=render_point(info->postlist[look->loneighbor[i-2]],
 987                                  info->postlist[look->hineighbor[i-2]],
 988                                  fit_value[look->loneighbor[i-2]],
 989                                  fit_value[look->hineighbor[i-2]],
 990                                  info->postlist[i]);
 991       int hiroom=look->quant_q-predicted;
 992       int loroom=predicted;
 993       int room=(hiroom<loroom?hiroom:loroom)<<1;
 994       int val=fit_value[i];
 995
 996       if(val){
 997         if(val>=room){
 998           if(hiroom>loroom){
 999             val = val-loroom;
1000           }else{
1001             val = -1-(val-hiroom);
1002           }
1003         }else{
1004           if(val&1){
1005             val= -((val+1)>>1);
1006           }else{
1007             val>>=1;
1008           }
1009         }
1010
1011         fit_value[i]=val+predicted;
1012         fit_value[look->loneighbor[i-2]]&=0x7fff;
1013         fit_value[look->hineighbor[i-2]]&=0x7fff;
1014
1015       }else{
1016         fit_value[i]=predicted|0x8000;
1017       }
1018
1019     }
1020
1021     return(fit_value);
1022   }
1023  eop:
1024   return(NULL);
1025 }
1026
1027 static int floor1_inverse2(vorbis_block *vb,vorbis_look_floor *in,void *memo,
1028                           float *out){
1029   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1030   vorbis_info_floor1 *info=look->vi;
1031
1032   codec_setup_info   *ci=vb->vd->vi->codec_setup;
1033   int                  n=ci->blocksizes[vb->W]/2;
1034   int j;
1035
1036   if(memo){
1037     /* render the lines */
1038     int *fit_value=(int *)memo;
1039     int hx=0;
1040     int lx=0;
1041     int ly=fit_value[0]*info->mult;
1042     for(j=1;j<look->posts;j++){
1043       int current=look->forward_index[j];
1044       int hy=fit_value[current]&0x7fff;
1045       if(hy==fit_value[current]){
1046
1047         hy*=info->mult;
1048         hx=info->postlist[current];
1049
1050         render_line(n,lx,hx,ly,hy,out);
1051
1052         lx=hx;
1053         ly=hy;
1054       }
1055     }
1056     for(j=hx;j<n;j++)out[j]*=FLOOR1_fromdB_LOOKUP[ly]; /* be certain */
1057     return(1);
1058   }
1059   memset(out,0,sizeof(*out)*n);
1060   return(0);
1061 }
1062
1063 /* export hooks */
1064 vorbis_func_floor floor1_exportbundle={
1065   &floor1_pack,&floor1_unpack,&floor1_look,&floor1_free_info,
1066   &floor1_free_look,&floor1_inverse1,&floor1_inverse2
1067 };