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   int quantized=vorbis_dBquant(flr+x0);
 403
 404   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;
 405
 406   memset(a,0,sizeof(*a));
 407   a->x0=x0;
 408   a->x1=x1;
 409   if(x1>=n)x1=n-1;
 410
 411   for(i=x0;i<=x1;i++){
 412     int quantized=vorbis_dBquant(flr+i);
 413     if(quantized){
 414       if(mdct[i]+info->twofitatten>=flr[i]){
 415         xa  += i;
 416         ya  += quantized;
 417         x2a += i*i;
 418         y2a += quantized*quantized;
 419         xya += i*quantized;
 420         na++;
 421       }else{
 422         xb  += i;
 423         yb  += quantized;
 424         x2b += i*i;
 425         y2b += quantized*quantized;
 426         xyb += i*quantized;
 427         nb++;
 428       }
 429     }
 430   }
 431
 432   xb+=xa;
 433   yb+=ya;
 434   x2b+=x2a;
 435   y2b+=y2a;
 436   xyb+=xya;
 437   nb+=na;
 438
 439   /* weight toward the actually used frequencies if we meet the threshhold */
 440   {
 441     int weight=nb*info->twofitweight/(na+1);
 442
 443     a->xa=xa*weight+xb;
 444     a->ya=ya*weight+yb;
 445     a->x2a=x2a*weight+x2b;
 446     a->y2a=y2a*weight+y2b;
 447     a->xya=xya*weight+xyb;
 448     a->an=na*weight+nb;
 449   }
 450
 451   return(na);
 452 }
 453
 454 static void fit_line(lsfit_acc *a,int fits,int *y0,int *y1){
 455   long x=0,y=0,x2=0,y2=0,xy=0,an=0,i;
 456   long x0=a[0].x0;
 457   long x1=a[fits-1].x1;
 458
 459   for(i=0;i<fits;i++){
 460     x+=a[i].xa;
 461     y+=a[i].ya;
 462     x2+=a[i].x2a;
 463     y2+=a[i].y2a;
 464     xy+=a[i].xya;
 465     an+=a[i].an;
 466   }
 467
 468   if(*y0>=0){
 469     x+=   x0;
 470     y+=  *y0;
 471     x2+=  x0 *  x0;
 472     y2+= *y0 * *y0;
 473     xy+= *y0 *  x0;
 474     an++;
 475   }
 476
 477   if(*y1>=0){
 478     x+=   x1;
 479     y+=  *y1;
 480     x2+=  x1 *  x1;
 481     y2+= *y1 * *y1;
 482     xy+= *y1 *  x1;
 483     an++;
 484   }
 485
 486   if(an){
 487     /* need 64 bit multiplies, which C doesn't give portably as int */
 488     double fx=x;
 489     double fy=y;
 490     double fx2=x2;
 491     double fxy=xy;
 492     double denom=1./(an*fx2-fx*fx);
 493     double a=(fy*fx2-fxy*fx)*denom;
 494     double b=(an*fxy-fx*fy)*denom;
 495     *y0=rint(a+b*x0);
 496     *y1=rint(a+b*x1);
 497
 498     /* limit to our range! */
 499     if(*y0>1023)*y0=1023;
 500     if(*y1>1023)*y1=1023;
 501     if(*y0<0)*y0=0;
 502     if(*y1<0)*y1=0;
 503
 504   }else{
 505     *y0=0;
 506     *y1=0;
 507   }
 508 }
 509
 510 /*static void fit_line_point(lsfit_acc *a,int fits,int *y0,int *y1){
 511   long y=0;
 512   int i;
 513
 514   for(i=0;i<fits && y==0;i++)
 515     y+=a[i].ya;
 516
 517   *y0=*y1=y;
 518   }*/
 519
 520 static int inspect_error(int x0,int x1,int y0,int y1,const float *mask,
 521                          const float *mdct,
 522                          vorbis_info_floor1 *info){
 523   int dy=y1-y0;
 524   int adx=x1-x0;
 525   int ady=abs(dy);
 526   int base=dy/adx;
 527   int sy=(dy<0?base-1:base+1);
 528   int x=x0;
 529   int y=y0;
 530   int err=0;
 531   int val=vorbis_dBquant(mask+x);
 532   int mse=0;
 533   int n=0;
 534
 535   ady-=abs(base*adx);
 536
 537   mse=(y-val);
 538   mse*=mse;
 539   n++;
 540   if(mdct[x]+info->twofitatten>=mask[x]){
 541     if(y+info->maxover<val)return(1);
 542     if(y-info->maxunder>val)return(1);
 543   }
 544
 545   while(++x<x1){
 546     err=err+ady;
 547     if(err>=adx){
 548       err-=adx;
 549       y+=sy;
 550     }else{
 551       y+=base;
 552     }
 553
 554     val=vorbis_dBquant(mask+x);
 555     mse+=((y-val)*(y-val));
 556     n++;
 557     if(mdct[x]+info->twofitatten>=mask[x]){
 558       if(val){
 559         if(y+info->maxover<val)return(1);
 560         if(y-info->maxunder>val)return(1);
 561       }
 562     }
 563   }
 564
 565   if(info->maxover*info->maxover/n>info->maxerr)return(0);
 566   if(info->maxunder*info->maxunder/n>info->maxerr)return(0);
 567   if(mse/n>info->maxerr)return(1);
 568   return(0);
 569 }
 570
 571 static int post_Y(int *A,int *B,int pos){
 572   if(A[pos]<0)
 573     return B[pos];
 574   if(B[pos]<0)
 575     return A[pos];
 576
 577   return (A[pos]+B[pos])>>1;
 578 }
 579
 580 static int seq=0;
 581
 582 int *floor1_fit(vorbis_block *vb,vorbis_look_floor1 *look,
 583                           const float *logmdct,   /* in */
 584                           const float *logmask){
 585   long i,j;
 586   vorbis_info_floor1 *info=look->vi;
 587   long n=look->n;
 588   long posts=look->posts;
 589   long nonzero=0;
 590   lsfit_acc fits[VIF_POSIT+1];
 591   int fit_valueA[VIF_POSIT+2]; /* index by range list position */
 592   int fit_valueB[VIF_POSIT+2]; /* index by range list position */
 593
 594   int loneighbor[VIF_POSIT+2]; /* sorted index of range list position (+2) */
 595   int hineighbor[VIF_POSIT+2];
 596   int *output=NULL;
 597   int memo[VIF_POSIT+2];
 598
 599   for(i=0;i<posts;i++)fit_valueA[i]=-200; /* mark all unused */
 600   for(i=0;i<posts;i++)fit_valueB[i]=-200; /* mark all unused */
 601   for(i=0;i<posts;i++)loneighbor[i]=0; /* 0 for the implicit 0 post */
 602   for(i=0;i<posts;i++)hineighbor[i]=1; /* 1 for the implicit post at n */
 603   for(i=0;i<posts;i++)memo[i]=-1;      /* no neighbor yet */
 604
 605   /* quantize the relevant floor points and collect them into line fit
 606      structures (one per minimal division) at the same time */
 607   if(posts==0){
 608     nonzero+=accumulate_fit(logmask,logmdct,0,n,fits,n,info);
 609   }else{
 610     for(i=0;i<posts-1;i++)
 611       nonzero+=accumulate_fit(logmask,logmdct,look->sorted_index[i],
 612                               look->sorted_index[i+1],fits+i,
 613                               n,info);
 614   }
 615
 616   if(nonzero){
 617     /* start by fitting the implicit base case.... */
 618     int y0=-200;
 619     int y1=-200;
 620     fit_line(fits,posts-1,&y0,&y1);
 621
 622     fit_valueA[0]=y0;
 623     fit_valueB[0]=y0;
 624     fit_valueB[1]=y1;
 625     fit_valueA[1]=y1;
 626
 627     /* Non degenerate case */
 628     /* start progressive splitting.  This is a greedy, non-optimal
 629        algorithm, but simple and close enough to the best
 630        answer. */
 631     for(i=2;i<posts;i++){
 632       int sortpos=look->reverse_index[i];
 633       int ln=loneighbor[sortpos];
 634       int hn=hineighbor[sortpos];
 635
 636       /* eliminate repeat searches of a particular range with a memo */
 637       if(memo[ln]!=hn){
 638         /* haven't performed this error search yet */
 639         int lsortpos=look->reverse_index[ln];
 640         int hsortpos=look->reverse_index[hn];
 641         memo[ln]=hn;
 642
 643         {
 644           /* A note: we want to bound/minimize *local*, not global, error */
 645           int lx=info->postlist[ln];
 646           int hx=info->postlist[hn];
 647           int ly=post_Y(fit_valueA,fit_valueB,ln);
 648           int hy=post_Y(fit_valueA,fit_valueB,hn);
 649
 650           if(ly==-1 || hy==-1){
 651             exit(1);
 652           }
 653
 654           if(inspect_error(lx,hx,ly,hy,logmask,logmdct,info)){
 655             /* outside error bounds/begin search area.  Split it. */
 656             int ly0=-200;
 657             int ly1=-200;
 658             int hy0=-200;
 659             int hy1=-200;
 660             fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1);
 661             fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1);
 662
 663             /* store new edge values */
 664             fit_valueB[ln]=ly0;
 665             if(ln==0)fit_valueA[ln]=ly0;
 666             fit_valueA[i]=ly1;
 667             fit_valueB[i]=hy0;
 668             fit_valueA[hn]=hy1;
 669             if(hn==1)fit_valueB[hn]=hy1;
 670
 671             if(ly1>=0 || hy0>=0){
 672               /* store new neighbor values */
 673               for(j=sortpos-1;j>=0;j--)
 674                 if(hineighbor[j]==hn)
 675                   hineighbor[j]=i;
 676                 else
 677                   break;
 678               for(j=sortpos+1;j<posts;j++)
 679                 if(loneighbor[j]==ln)
 680                   loneighbor[j]=i;
 681                 else
 682                   break;
 683
 684             }
 685           }else{
 686
 687             fit_valueA[i]=-200;
 688             fit_valueB[i]=-200;
 689           }
 690         }
 691       }
 692     }
 693
 694     output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
 695
 696     output[0]=post_Y(fit_valueA,fit_valueB,0);
 697     output[1]=post_Y(fit_valueA,fit_valueB,1);
 698
 699     /* fill in posts marked as not using a fit; we will zero
 700        back out to 'unused' when encoding them so long as curve
 701        interpolation doesn't force them into use */
 702     for(i=2;i<posts;i++){
 703       int ln=look->loneighbor[i-2];
 704       int hn=look->hineighbor[i-2];
 705       int x0=info->postlist[ln];
 706       int x1=info->postlist[hn];
 707       int y0=output[ln];
 708       int y1=output[hn];
 709
 710       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
 711       int vx=post_Y(fit_valueA,fit_valueB,i);
 712
 713       if(vx>=0 && predicted!=vx){
 714         output[i]=vx;
 715       }else{
 716         output[i]= predicted|0x8000;
 717       }
 718     }
 719   }
 720
 721   return(output);
 722
 723 }
 724
 725 int *floor1_interpolate_fit(vorbis_block *vb,vorbis_look_floor1 *look,
 726                           int *A,int *B,
 727                           int del){
 728
 729   long i;
 730   long posts=look->posts;
 731   int *output=NULL;
 732
 733   if(A && B){
 734     output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
 735
 736     for(i=0;i<posts;i++){
 737       output[i]=((65536-del)*(A[i]&0x7fff)+del*(B[i]&0x7fff)+32768)>>16;
 738       if(A[i]&0x8000 && B[i]&0x8000)output[i]|=0x8000;
 739     }
 740   }
 741
 742   return(output);
 743 }
 744
 745
 746 int floor1_encode(oggpack_buffer *opb,vorbis_block *vb,
 747                   vorbis_look_floor1 *look,
 748                   int *post,int *ilogmask){
 749
 750   long i,j;
 751   vorbis_info_floor1 *info=look->vi;
 752   long n=look->n;
 753   long posts=look->posts;
 754   codec_setup_info *ci=vb->vd->vi->codec_setup;
 755   int out[VIF_POSIT+2];
 756   static_codebook **sbooks=ci->book_param;
 757   codebook *books=ci->fullbooks;
 758   static long seq=0;
 759
 760   /* quantize values to multiplier spec */
 761   if(post){
 762     for(i=0;i<posts;i++){
 763       int val=post[i]&0x7fff;
 764       switch(info->mult){
 765       case 1: /* 1024 -> 256 */
 766         val>>=2;
 767         break;
 768       case 2: /* 1024 -> 128 */
 769         val>>=3;
 770         break;
 771       case 3: /* 1024 -> 86 */
 772         val/=12;
 773         break;
 774       case 4: /* 1024 -> 64 */
 775         val>>=4;
 776         break;
 777       }
 778       post[i]=val | (post[i]&0x8000);
 779     }
 780
 781     out[0]=post[0];
 782     out[1]=post[1];
 783
 784     /* find prediction values for each post and subtract them */
 785     for(i=2;i<posts;i++){
 786       int ln=look->loneighbor[i-2];
 787       int hn=look->hineighbor[i-2];
 788       int x0=info->postlist[ln];
 789       int x1=info->postlist[hn];
 790       int y0=post[ln];
 791       int y1=post[hn];
 792
 793       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
 794
 795       if((post[i]&0x8000) || (predicted==post[i])){
 796         post[i]=predicted|0x8000; /* in case there was roundoff jitter
 797                                      in interpolation */
 798         out[i]=0;
 799       }else{
 800         int headroom=(look->quant_q-predicted<predicted?
 801                       look->quant_q-predicted:predicted);
 802
 803         int val=post[i]-predicted;
 804
 805         /* at this point the 'deviation' value is in the range +/- max
 806            range, but the real, unique range can always be mapped to
 807            only [0-maxrange).  So we want to wrap the deviation into
 808            this limited range, but do it in the way that least screws
 809            an essentially gaussian probability distribution. */
 810
 811         if(val<0)
 812           if(val<-headroom)
 813             val=headroom-val-1;
 814           else
 815             val=-1-(val<<1);
 816         else
 817           if(val>=headroom)
 818             val= val+headroom;
 819           else
 820             val<<=1;
 821
 822         out[i]=val;
 823         post[ln]&=0x7fff;
 824         post[hn]&=0x7fff;
 825       }
 826     }
 827
 828     /* we have everything we need. pack it out */
 829     /* mark nontrivial floor */
 830     oggpack_write(opb,1,1);
 831
 832     /* beginning/end post */
 833     look->frames++;
 834     look->postbits+=ilog(look->quant_q-1)*2;
 835     oggpack_write(opb,out[0],ilog(look->quant_q-1));
 836     oggpack_write(opb,out[1],ilog(look->quant_q-1));
 837
 838
 839     /* partition by partition */
 840     for(i=0,j=2;i<info->partitions;i++){
 841       int class=info->partitionclass[i];
 842       int cdim=info->class_dim[class];
 843       int csubbits=info->class_subs[class];
 844       int csub=1<<csubbits;
 845       int bookas[8]={0,0,0,0,0,0,0,0};
 846       int cval=0;
 847       int cshift=0;
 848       int k,l;
 849
 850       /* generate the partition's first stage cascade value */
 851       if(csubbits){
 852         int maxval[8];
 853         for(k=0;k<csub;k++){
 854           int booknum=info->class_subbook[class][k];
 855           if(booknum<0){
 856             maxval[k]=1;
 857           }else{
 858             maxval[k]=sbooks[info->class_subbook[class][k]]->entries;
 859           }
 860         }
 861         for(k=0;k<cdim;k++){
 862           for(l=0;l<csub;l++){
 863             int val=out[j+k];
 864             if(val<maxval[l]){
 865               bookas[k]=l;
 866               break;
 867             }
 868           }
 869           cval|= bookas[k]<<cshift;
 870           cshift+=csubbits;
 871         }
 872         /* write it */
 873         look->phrasebits+=
 874           vorbis_book_encode(books+info->class_book[class],cval,opb);
 875
 876 #ifdef TRAIN_FLOOR1
 877         {
 878           FILE *of;
 879           char buffer[80];
 880           sprintf(buffer,"line_%dx%ld_class%d.vqd",
 881                   vb->pcmend/2,posts-2,class);
 882           of=fopen(buffer,"a");
 883           fprintf(of,"%d\n",cval);
 884           fclose(of);
 885         }
 886 #endif
 887       }
 888
 889       /* write post values */
 890       for(k=0;k<cdim;k++){
 891         int book=info->class_subbook[class][bookas[k]];
 892         if(book>=0){
 893           /* hack to allow training with 'bad' books */
 894           if(out[j+k]<(books+book)->entries)
 895             look->postbits+=vorbis_book_encode(books+book,
 896                                                out[j+k],opb);
 897           /*else
 898             fprintf(stderr,"+!");*/
 899
 900 #ifdef TRAIN_FLOOR1
 901           {
 902             FILE *of;
 903             char buffer[80];
 904             sprintf(buffer,"line_%dx%ld_%dsub%d.vqd",
 905                     vb->pcmend/2,posts-2,class,bookas[k]);
 906             of=fopen(buffer,"a");
 907             fprintf(of,"%d\n",out[j+k]);
 908             fclose(of);
 909           }
 910 #endif
 911         }
 912       }
 913       j+=cdim;
 914     }
 915
 916     {
 917       /* generate quantized floor equivalent to what we'd unpack in decode */
 918       /* render the lines */
 919       int hx=0;
 920       int lx=0;
 921       int ly=post[0]*info->mult;
 922       for(j=1;j<look->posts;j++){
 923         int current=look->forward_index[j];
 924         int hy=post[current]&0x7fff;
 925         if(hy==post[current]){
 926
 927           hy*=info->mult;
 928           hx=info->postlist[current];
 929
 930           render_line0(lx,hx,ly,hy,ilogmask);
 931
 932           lx=hx;
 933           ly=hy;
 934         }
 935       }
 936       for(j=hx;j<vb->pcmend/2;j++)ilogmask[j]=ly; /* be certain */
 937       seq++;
 938       return(1);
 939     }
 940   }else{
 941     oggpack_write(opb,0,1);
 942     memset(ilogmask,0,vb->pcmend/2*sizeof(*ilogmask));
 943     seq++;
 944     return(0);
 945   }
 946 }
 947
 948 static void *floor1_inverse1(vorbis_block *vb,vorbis_look_floor *in){
 949   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
 950   vorbis_info_floor1 *info=look->vi;
 951   codec_setup_info   *ci=vb->vd->vi->codec_setup;
 952
 953   int i,j,k;
 954   codebook *books=ci->fullbooks;
 955
 956   /* unpack wrapped/predicted values from stream */
 957   if(oggpack_read(&vb->opb,1)==1){
 958     int *fit_value=_vorbis_block_alloc(vb,(look->posts)*sizeof(*fit_value));
 959
 960     fit_value[0]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
 961     fit_value[1]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
 962
 963     /* partition by partition */
 964     for(i=0,j=2;i<info->partitions;i++){
 965       int class=info->partitionclass[i];
 966       int cdim=info->class_dim[class];
 967       int csubbits=info->class_subs[class];
 968       int csub=1<<csubbits;
 969       int cval=0;
 970
 971       /* decode the partition's first stage cascade value */
 972       if(csubbits){
 973         cval=vorbis_book_decode(books+info->class_book[class],&vb->opb);
 974
 975         if(cval==-1)goto eop;
 976       }
 977
 978       for(k=0;k<cdim;k++){
 979         int book=info->class_subbook[class][cval&(csub-1)];
 980         cval>>=csubbits;
 981         if(book>=0){
 982           if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1)
 983             goto eop;
 984         }else{
 985           fit_value[j+k]=0;
 986         }
 987       }
 988       j+=cdim;
 989     }
 990
 991     /* unwrap positive values and reconsitute via linear interpolation */
 992     for(i=2;i<look->posts;i++){
 993       int predicted=render_point(info->postlist[look->loneighbor[i-2]],
 994                                  info->postlist[look->hineighbor[i-2]],
 995                                  fit_value[look->loneighbor[i-2]],
 996                                  fit_value[look->hineighbor[i-2]],
 997                                  info->postlist[i]);
 998       int hiroom=look->quant_q-predicted;
 999       int loroom=predicted;
1000       int room=(hiroom<loroom?hiroom:loroom)<<1;
1001       int val=fit_value[i];
1002
1003       if(val){
1004         if(val>=room){
1005           if(hiroom>loroom){
1006             val = val-loroom;
1007           }else{
1008             val = -1-(val-hiroom);
1009           }
1010         }else{
1011           if(val&1){
1012             val= -((val+1)>>1);
1013           }else{
1014             val>>=1;
1015           }
1016         }
1017
1018         fit_value[i]=val+predicted;
1019         fit_value[look->loneighbor[i-2]]&=0x7fff;
1020         fit_value[look->hineighbor[i-2]]&=0x7fff;
1021
1022       }else{
1023         fit_value[i]=predicted|0x8000;
1024       }
1025
1026     }
1027
1028     return(fit_value);
1029   }
1030  eop:
1031   return(NULL);
1032 }
1033
1034 static int floor1_inverse2(vorbis_block *vb,vorbis_look_floor *in,void *memo,
1035                           float *out){
1036   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1037   vorbis_info_floor1 *info=look->vi;
1038
1039   codec_setup_info   *ci=vb->vd->vi->codec_setup;
1040   int                  n=ci->blocksizes[vb->W]/2;
1041   int j;
1042
1043   if(memo){
1044     /* render the lines */
1045     int *fit_value=(int *)memo;
1046     int hx=0;
1047     int lx=0;
1048     int ly=fit_value[0]*info->mult;
1049     for(j=1;j<look->posts;j++){
1050       int current=look->forward_index[j];
1051       int hy=fit_value[current]&0x7fff;
1052       if(hy==fit_value[current]){
1053
1054         hy*=info->mult;
1055         hx=info->postlist[current];
1056
1057         render_line(n,lx,hx,ly,hy,out);
1058
1059         lx=hx;
1060         ly=hy;
1061       }
1062     }
1063     for(j=hx;j<n;j++)out[j]*=FLOOR1_fromdB_LOOKUP[ly]; /* be certain */
1064     return(1);
1065   }
1066   memset(out,0,sizeof(*out)*n);
1067   return(0);
1068 }
1069
1070 /* export hooks */
1071 vorbis_func_floor floor1_exportbundle={
1072   &floor1_pack,&floor1_unpack,&floor1_look,&floor1_free_info,
1073   &floor1_free_look,&floor1_inverse1,&floor1_inverse2
1074 };
1075
1076