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