lib/sharedbook.c

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