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