lib/sharedbook.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
   6  * PLEASE READ THESE TERMS DISTRIBUTING.                            *
   7  *                                                                  *
   8  * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000             *
   9  * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: basic shared codebook operations
  15  last mod: $Id: sharedbook.c,v 1.9 2000/10/12 03:12:54 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdlib.h>
  20 #include <math.h>
  21 #include <string.h>
  22 #include <ogg/ogg.h>
  23 #include "os.h"
  24 #include "vorbis/codec.h"
  25 #include "vorbis/codebook.h"
  26 #include "scales.h"
  27 #include "sharedbook.h"
  28
  29 /**** pack/unpack helpers ******************************************/
  30 int _ilog(unsigned int v){
  31   int ret=0;
  32   while(v){
  33     ret++;
  34     v>>=1;
  35   }
  36   return(ret);
  37 }
  38
  39 /* 32 bit float (not IEEE; nonnormalized mantissa +
  40    biased exponent) : neeeeeee eeemmmmm mmmmmmmm mmmmmmmm
  41    Why not IEEE?  It's just not that important here. */
  42
  43 #define VQ_FEXP 10
  44 #define VQ_FMAN 21
  45 #define VQ_FEXP_BIAS 768 /* bias toward values smaller than 1. */
  46
  47 /* doesn't currently guard under/overflow */
  48 long _float32_pack(float val){
  49   int sign=0;
  50   long exp;
  51   long mant;
  52   if(val<0){
  53     sign=0x80000000;
  54     val= -val;
  55   }
  56   exp= floor(log(val)/log(2));
  57   mant=rint(ldexp(val,(VQ_FMAN-1)-exp));
  58   exp=(exp+VQ_FEXP_BIAS)<<VQ_FMAN;
  59
  60   return(sign|exp|mant);
  61 }
  62
  63 float _float32_unpack(long val){
  64   float mant=val&0x1fffff;
  65   float sign=val&0x80000000;
  66   float exp =(val&0x7fe00000)>>VQ_FMAN;
  67   if(sign)mant= -mant;
  68   return(ldexp(mant,exp-(VQ_FMAN-1)-VQ_FEXP_BIAS));
  69 }
  70
  71 /* given a list of word lengths, generate a list of codewords.  Works
  72    for length ordered or unordered, always assigns the lowest valued
  73    codewords first.  Extended to handle unused entries (length 0) */
  74 long *_make_words(long *l,long n){
  75   long i,j;
  76   long marker[33];
  77   long *r=malloc(n*sizeof(long));
  78   memset(marker,0,sizeof(marker));
  79
  80   for(i=0;i<n;i++){
  81     long length=l[i];
  82     if(length>0){
  83       long entry=marker[length];
  84
  85       /* when we claim a node for an entry, we also claim the nodes
  86          below it (pruning off the imagined tree that may have dangled
  87          from it) as well as blocking the use of any nodes directly
  88          above for leaves */
  89
  90       /* update ourself */
  91       if(length<32 && (entry>>length)){
  92         /* error condition; the lengths must specify an overpopulated tree */
  93         free(r);
  94         return(NULL);
  95       }
  96       r[i]=entry;
  97
  98       /* Look to see if the next shorter marker points to the node
  99          above. if so, update it and repeat.  */
 100       {
 101         for(j=length;j>0;j--){
 102
 103           if(marker[j]&1){
 104             /* have to jump branches */
 105             if(j==1)
 106               marker[1]++;
 107             else
 108               marker[j]=marker[j-1]<<1;
 109             break; /* invariant says next upper marker would already
 110                       have been moved if it was on the same path */
 111           }
 112           marker[j]++;
 113         }
 114       }
 115
 116       /* prune the tree; the implicit invariant says all the longer
 117          markers were dangling from our just-taken node.  Dangle them
 118          from our *new* node. */
 119       for(j=length+1;j<33;j++)
 120         if((marker[j]>>1) == entry){
 121           entry=marker[j];
 122           marker[j]=marker[j-1]<<1;
 123         }else
 124           break;
 125     }
 126   }
 127
 128   /* bitreverse the words because our bitwise packer/unpacker is LSb
 129      endian */
 130   for(i=0;i<n;i++){
 131     long temp=0;
 132     for(j=0;j<l[i];j++){
 133       temp<<=1;
 134       temp|=(r[i]>>j)&1;
 135     }
 136     r[i]=temp;
 137   }
 138
 139   return(r);
 140 }
 141
 142 /* build the decode helper tree from the codewords */
 143 decode_aux *_make_decode_tree(codebook *c){
 144   const static_codebook *s=c->c;
 145   long top=0,i,j,n;
 146   decode_aux *t=malloc(sizeof(decode_aux));
 147   long *ptr0=t->ptr0=calloc(c->entries*2,sizeof(long));
 148   long *ptr1=t->ptr1=calloc(c->entries*2,sizeof(long));
 149   long *codelist=_make_words(s->lengthlist,s->entries);
 150
 151   if(codelist==NULL)return(NULL);
 152   t->aux=c->entries*2;
 153
 154   for(i=0;i<c->entries;i++){
 155     if(s->lengthlist[i]>0){
 156       long ptr=0;
 157       for(j=0;j<s->lengthlist[i]-1;j++){
 158         int bit=(codelist[i]>>j)&1;
 159         if(!bit){
 160           if(!ptr0[ptr])
 161             ptr0[ptr]= ++top;
 162           ptr=ptr0[ptr];
 163         }else{
 164           if(!ptr1[ptr])
 165             ptr1[ptr]= ++top;
 166           ptr=ptr1[ptr];
 167         }
 168       }
 169       if(!((codelist[i]>>j)&1))
 170         ptr0[ptr]=-i;
 171       else
 172         ptr1[ptr]=-i;
 173     }
 174   }
 175   free(codelist);
 176
 177   t->tabn = _ilog(c->entries)-4; /* this is magic */
 178   if(t->tabn<5)t->tabn=5;
 179   n = 1<<t->tabn;
 180   t->tab = malloc(n*sizeof(long));
 181   t->tabl = malloc(n*sizeof(int));
 182   for (i = 0; i < n; i++) {
 183     long p = 0;
 184     for (j = 0; j < t->tabn && (p > 0 || j == 0); j++) {
 185       if (i & (1 << j))
 186         p = ptr1[p];
 187       else
 188         p = ptr0[p];
 189     }
 190     /* now j == length, and p == -code */
 191     t->tab[i] = p;
 192     t->tabl[i] = j;
 193   }
 194
 195   return(t);
 196 }
 197
 198 /* there might be a straightforward one-line way to do the below
 199    that's portable and totally safe against roundoff, but I haven't
 200    thought of it.  Therefore, we opt on the side of caution */
 201 long _book_maptype1_quantvals(const static_codebook *b){
 202   long vals=floor(pow(b->entries,1./b->dim));
 203
 204   /* the above *should* be reliable, but we'll not assume that FP is
 205      ever reliable when bitstream sync is at stake; verify via integer
 206      means that vals really is the greatest value of dim for which
 207      vals^b->bim <= b->entries */
 208   /* treat the above as an initial guess */
 209   while(1){
 210     long acc=1;
 211     long acc1=1;
 212     int i;
 213     for(i=0;i<b->dim;i++){
 214       acc*=vals;
 215       acc1*=vals+1;
 216     }
 217     if(acc<=b->entries && acc1>b->entries){
 218       return(vals);
 219     }else{
 220       if(acc>b->entries){
 221         vals--;
 222       }else{
 223         vals++;
 224       }
 225     }
 226   }
 227 }
 228
 229 /* unpack the quantized list of values for encode/decode ***********/
 230 /* we need to deal with two map types: in map type 1, the values are
 231    generated algorithmically (each column of the vector counts through
 232    the values in the quant vector). in map type 2, all the values came
 233    in in an explicit list.  Both value lists must be unpacked */
 234 float *_book_unquantize(const static_codebook *b){
 235   long j,k;
 236   if(b->maptype==1 || b->maptype==2){
 237     int quantvals;
 238     float mindel=_float32_unpack(b->q_min);
 239     float delta=_float32_unpack(b->q_delta);
 240     float *r=calloc(b->entries*b->dim,sizeof(float));
 241
 242     /* maptype 1 and 2 both use a quantized value vector, but
 243        different sizes */
 244     switch(b->maptype){
 245     case 1:
 246       /* most of the time, entries%dimensions == 0, but we need to be
 247          well defined.  We define that the possible vales at each
 248          scalar is values == entries/dim.  If entries%dim != 0, we'll
 249          have 'too few' values (values*dim<entries), which means that
 250          we'll have 'left over' entries; left over entries use zeroed
 251          values (and are wasted).  So don't generate codebooks like
 252          that */
 253       quantvals=_book_maptype1_quantvals(b);
 254       for(j=0;j<b->entries;j++){
 255         float last=0.;
 256         int indexdiv=1;
 257         for(k=0;k<b->dim;k++){
 258           int index= (j/indexdiv)%quantvals;
 259           float val=b->quantlist[index];
 260           val=fabs(val)*delta+mindel+last;
 261           if(b->q_sequencep)last=val;
 262           r[j*b->dim+k]=val;
 263           indexdiv*=quantvals;
 264         }
 265       }
 266       break;
 267     case 2:
 268       for(j=0;j<b->entries;j++){
 269         float last=0.;
 270         for(k=0;k<b->dim;k++){
 271           float val=b->quantlist[j*b->dim+k];
 272           val=fabs(val)*delta+mindel+last;
 273           if(b->q_sequencep)last=val;
 274           r[j*b->dim+k]=val;
 275         }
 276       }
 277     }
 278     return(r);
 279   }
 280   return(NULL);
 281 }
 282
 283 void vorbis_staticbook_clear(static_codebook *b){
 284   if(b->quantlist)free(b->quantlist);
 285   if(b->lengthlist)free(b->lengthlist);
 286   if(b->nearest_tree){
 287     free(b->nearest_tree->ptr0);
 288     free(b->nearest_tree->ptr1);
 289     free(b->nearest_tree->p);
 290     free(b->nearest_tree->q);
 291     memset(b->nearest_tree,0,sizeof(encode_aux_nearestmatch));
 292     free(b->nearest_tree);
 293   }
 294   if(b->thresh_tree){
 295     free(b->thresh_tree->quantthresh);
 296     free(b->thresh_tree->quantmap);
 297     memset(b->thresh_tree,0,sizeof(encode_aux_threshmatch));
 298     free(b->thresh_tree);
 299   }
 300   memset(b,0,sizeof(static_codebook));
 301 }
 302
 303 void vorbis_book_clear(codebook *b){
 304   /* static book is not cleared; we're likely called on the lookup and
 305      the static codebook belongs to the info struct */
 306   if(b->decode_tree){
 307     free(b->decode_tree->ptr0);
 308     free(b->decode_tree->ptr1);
 309     memset(b->decode_tree,0,sizeof(decode_aux));
 310     free(b->decode_tree);
 311   }
 312   if(b->valuelist)free(b->valuelist);
 313   if(b->codelist)free(b->codelist);
 314   memset(b,0,sizeof(codebook));
 315 }
 316
 317 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
 318   memset(c,0,sizeof(codebook));
 319   c->c=s;
 320   c->entries=s->entries;
 321   c->dim=s->dim;
 322   c->codelist=_make_words(s->lengthlist,s->entries);
 323   c->valuelist=_book_unquantize(s);
 324   return(0);
 325 }
 326
 327 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
 328   memset(c,0,sizeof(codebook));
 329   c->c=s;
 330   c->entries=s->entries;
 331   c->dim=s->dim;
 332   c->valuelist=_book_unquantize(s);
 333   c->decode_tree=_make_decode_tree(c);
 334   if(c->decode_tree==NULL)goto err_out;
 335   return(0);
 336  err_out:
 337   vorbis_book_clear(c);
 338   return(-1);
 339 }
 340
 341 static float _dist(int el,float *ref, float *b,int step){
 342   int i;
 343   float acc=0.;
 344   for(i=0;i<el;i++){
 345     float val=(ref[i]-b[i*step]);
 346     acc+=val*val;
 347   }
 348   return(acc);
 349 }
 350
 351 #include <stdio.h>
 352 int _best(codebook *book, float *a, int step){
 353   encode_aux_nearestmatch *nt=book->c->nearest_tree;
 354   encode_aux_threshmatch *tt=book->c->thresh_tree;
 355   encode_aux_pigeonhole *pt=book->c->pigeon_tree;
 356   int dim=book->dim;
 357   int ptr=0,k,o;
 358   /*int savebest=-1;
 359     float saverr;*/
 360
 361   /* do we have a threshhold encode hint? */
 362   if(tt){
 363     int index=0;
 364     /* find the quant val of each scalar */
 365     for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
 366       int i;
 367       /* linear search the quant list for now; it's small and although
 368          with > 8 entries, it would be faster to bisect, this would be
 369          a misplaced optimization for now */
 370       for(i=0;i<tt->threshvals-1;i++)
 371         if(a[o]<tt->quantthresh[i])break;
 372
 373       index=(index*tt->quantvals)+tt->quantmap[i];
 374     }
 375     /* regular lattices are easy :-) */
 376     if(book->c->lengthlist[index]>0) /* is this unused?  If so, we'll
 377                                         use a decision tree after all
 378                                         and fall through*/
 379       return(index);
 380   }
 381
 382   /* do we have a pigeonhole encode hint? */
 383   if(pt){
 384     const static_codebook *c=book->c;
 385     int i,besti=-1;
 386     float best;
 387     int entry=0;
 388
 389     /* dealing with sequentialness is a pain in the ass */
 390     if(c->q_sequencep){
 391       int pv;
 392       long mul=1;
 393       float qlast=0;
 394       for(k=0,o=0;k<dim;k++,o+=step){
 395         pv=(int)((a[o]-qlast-pt->min)/pt->del);
 396         if(pv<0 || pv>=pt->mapentries)break;
 397         entry+=pt->pigeonmap[pv]*mul;
 398         mul*=pt->quantvals;
 399         qlast+=pv*pt->del+pt->min;
 400       }
 401     }else{
 402       for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
 403         int pv=(int)((a[o]-pt->min)/pt->del);
 404         if(pv<0 || pv>=pt->mapentries)break;
 405         entry=entry*pt->quantvals+pt->pigeonmap[pv];
 406       }
 407     }
 408
 409     /* must be within the pigeonholable range; if we quant outside (or
 410        in an entry that we define no list for), brute force it */
 411     if(k==dim && pt->fitlength[entry]){
 412       /* search the abbreviated list */
 413       long *list=pt->fitlist+pt->fitmap[entry];
 414       for(i=0;i<pt->fitlength[entry];i++){
 415         float this=_dist(dim,book->valuelist+list[i]*dim,a,step);
 416         if(besti==-1 || this<best){
 417           best=this;
 418           besti=list[i];
 419         }
 420       }
 421
 422       return(besti);
 423     }
 424   }
 425
 426   if(nt){
 427     /* optimized using the decision tree */
 428     while(1){
 429       float c=0.;
 430       float *p=book->valuelist+nt->p[ptr];
 431       float *q=book->valuelist+nt->q[ptr];
 432
 433       for(k=0,o=0;k<dim;k++,o+=step)
 434         c+=(p[k]-q[k])*(a[o]-(p[k]+q[k])*.5);
 435
 436       if(c>0.) /* in A */
 437         ptr= -nt->ptr0[ptr];
 438       else     /* in B */
 439         ptr= -nt->ptr1[ptr];
 440       if(ptr<=0)break;
 441     }
 442     return(-ptr);
 443   }
 444
 445   /* brute force it! */
 446   {
 447     const static_codebook *c=book->c;
 448     int i,besti=-1;
 449     float best;
 450     float *e=book->valuelist;
 451     for(i=0;i<book->entries;i++){
 452       if(c->lengthlist[i]>0){
 453         float this=_dist(dim,e,a,step);
 454         if(besti==-1 || this<best){
 455           best=this;
 456           besti=i;
 457         }
 458       }
 459       e+=dim;
 460     }
 461
 462     /*if(savebest!=-1 && savebest!=besti){
 463       fprintf(stderr,"brute force/pigeonhole disagreement:\n"
 464               "original:");
 465       for(i=0;i<dim*step;i+=step)fprintf(stderr,"%g,",a[i]);
 466       fprintf(stderr,"\n"
 467               "pigeonhole (entry %d, err %g):",savebest,saverr);
 468       for(i=0;i<dim;i++)fprintf(stderr,"%g,",
 469                                 (book->valuelist+savebest*dim)[i]);
 470       fprintf(stderr,"\n"
 471               "bruteforce (entry %d, err %g):",besti,best);
 472       for(i=0;i<dim;i++)fprintf(stderr,"%g,",
 473                                 (book->valuelist+besti*dim)[i]);
 474       fprintf(stderr,"\n");
 475       }*/
 476     return(besti);
 477   }
 478 }
 479
 480 /* returns the entry number and *modifies a* to the remainder value ********/
 481 int vorbis_book_besterror(codebook *book,float *a,int step,int addmul){
 482   int dim=book->dim,i,o;
 483   int best=_best(book,a,step);
 484   switch(addmul){
 485   case 0:
 486     for(i=0,o=0;i<dim;i++,o+=step)
 487       a[o]-=(book->valuelist+best*dim)[i];
 488     break;
 489   case 1:
 490     for(i=0,o=0;i<dim;i++,o+=step){
 491       float val=(book->valuelist+best*dim)[i];
 492       if(val==0){
 493         a[o]=0;
 494       }else{
 495         a[o]/=val;
 496       }
 497     }
 498     break;
 499   }
 500   return(best);
 501 }
 502
 503 long vorbis_book_codeword(codebook *book,int entry){
 504   return book->codelist[entry];
 505 }
 506
 507 long vorbis_book_codelen(codebook *book,int entry){
 508   return book->c->lengthlist[entry];
 509 }
 510
 511 #ifdef _V_SELFTEST
 512
 513 /* Unit tests of the dequantizer; this stuff will be OK
 514    cross-platform, I simply want to be sure that special mapping cases
 515    actually work properly; a bug could go unnoticed for a while */
 516
 517 #include <stdio.h>
 518
 519 /* cases:
 520
 521    no mapping
 522    full, explicit mapping
 523    algorithmic mapping
 524
 525    nonsequential
 526    sequential
 527 */
 528
 529 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
 530 static long partial_quantlist1[]={0,7,2};
 531
 532 /* no mapping */
 533 static_codebook test1={
 534   4,16,
 535   NULL,
 536   0,
 537   0,0,0,0,
 538   NULL,
 539   NULL,NULL
 540 };
 541 static float *test1_result=NULL;
 542
 543 /* linear, full mapping, nonsequential */
 544 static_codebook test2={
 545   4,3,
 546   NULL,
 547   2,
 548   -533200896,1611661312,4,0,
 549   full_quantlist1,
 550   NULL,NULL
 551 };
 552 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
 553
 554 /* linear, full mapping, sequential */
 555 static_codebook test3={
 556   4,3,
 557   NULL,
 558   2,
 559   -533200896,1611661312,4,1,
 560   full_quantlist1,
 561   NULL,NULL
 562 };
 563 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
 564
 565 /* linear, algorithmic mapping, nonsequential */
 566 static_codebook test4={
 567   3,27,
 568   NULL,
 569   1,
 570   -533200896,1611661312,4,0,
 571   partial_quantlist1,
 572   NULL,NULL
 573 };
 574 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
 575                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
 576                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
 577                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
 578                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
 579                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
 580                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
 581                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
 582                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
 583
 584 /* linear, algorithmic mapping, sequential */
 585 static_codebook test5={
 586   3,27,
 587   NULL,
 588   1,
 589   -533200896,1611661312,4,1,
 590   partial_quantlist1,
 591   NULL,NULL
 592 };
 593 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
 594                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
 595                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
 596                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
 597                               -3, 1, 5, 4, 8,12, -1, 3, 7,
 598                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
 599                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
 600                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
 601                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
 602
 603 void run_test(static_codebook *b,float *comp){
 604   float *out=_book_unquantize(b);
 605   int i;
 606
 607   if(comp){
 608     if(!out){
 609       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
 610       exit(1);
 611     }
 612
 613     for(i=0;i<b->entries*b->dim;i++)
 614       if(fabs(out[i]-comp[i])>.0001){
 615         fprintf(stderr,"disagreement in unquantized and reference data:\n"
 616                 "position %d, %g != %g\n",i,out[i],comp[i]);
 617         exit(1);
 618       }
 619
 620   }else{
 621     if(out){
 622       fprintf(stderr,"_book_unquantize returned a value array: \n"
 623               " correct result should have been NULL\n");
 624       exit(1);
 625     }
 626   }
 627 }
 628
 629 int main(){
 630   /* run the nine dequant tests, and compare to the hand-rolled results */
 631   fprintf(stderr,"Dequant test 1... ");
 632   run_test(&test1,test1_result);
 633   fprintf(stderr,"OK\nDequant test 2... ");
 634   run_test(&test2,test2_result);
 635   fprintf(stderr,"OK\nDequant test 3... ");
 636   run_test(&test3,test3_result);
 637   fprintf(stderr,"OK\nDequant test 4... ");
 638   run_test(&test4,test4_result);
 639   fprintf(stderr,"OK\nDequant test 5... ");
 640   run_test(&test5,test5_result);
 641   fprintf(stderr,"OK\n\n");
 642
 643   return(0);
 644 }
 645
 646 #endif