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.11 2000/11/08 13:16:27 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       break;
 277     }
 278
 279     return(r);
 280   }
 281   return(NULL);
 282 }
 283
 284 void vorbis_staticbook_clear(static_codebook *b){
 285   if(b->allocedp){
 286     if(b->quantlist)free(b->quantlist);
 287     if(b->lengthlist)free(b->lengthlist);
 288     if(b->nearest_tree){
 289       free(b->nearest_tree->ptr0);
 290       free(b->nearest_tree->ptr1);
 291       free(b->nearest_tree->p);
 292       free(b->nearest_tree->q);
 293       memset(b->nearest_tree,0,sizeof(encode_aux_nearestmatch));
 294       free(b->nearest_tree);
 295     }
 296     if(b->thresh_tree){
 297       free(b->thresh_tree->quantthresh);
 298       free(b->thresh_tree->quantmap);
 299       memset(b->thresh_tree,0,sizeof(encode_aux_threshmatch));
 300       free(b->thresh_tree);
 301     }
 302
 303     memset(b,0,sizeof(static_codebook));
 304   }
 305 }
 306
 307 void vorbis_staticbook_destroy(static_codebook *b){
 308   if(b->allocedp){
 309     vorbis_staticbook_clear(b);
 310     free(b);
 311   }
 312 }
 313
 314 void vorbis_book_clear(codebook *b){
 315   /* static book is not cleared; we're likely called on the lookup and
 316      the static codebook belongs to the info struct */
 317   if(b->decode_tree){
 318     free(b->decode_tree->ptr0);
 319     free(b->decode_tree->ptr1);
 320     memset(b->decode_tree,0,sizeof(decode_aux));
 321     free(b->decode_tree);
 322   }
 323   if(b->valuelist)free(b->valuelist);
 324   if(b->codelist)free(b->codelist);
 325   memset(b,0,sizeof(codebook));
 326 }
 327
 328 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
 329   long j,k;
 330   memset(c,0,sizeof(codebook));
 331   c->c=s;
 332   c->entries=s->entries;
 333   c->dim=s->dim;
 334   c->codelist=_make_words(s->lengthlist,s->entries);
 335   c->valuelist=_book_unquantize(s);
 336
 337   /* set the 'zero entry' */
 338   c->zeroentry=-1;
 339   if(c->valuelist){
 340     for(j=0;j<s->entries;j++){
 341       int flag=1;
 342       for(k=0;k<s->dim;k++){
 343         if(fabs(c->valuelist[j*s->dim+k])>1e-12){
 344           flag=0;
 345           break;
 346         }
 347       }
 348       if(flag)
 349         c->zeroentry=j;
 350     }
 351   }
 352
 353   return(0);
 354 }
 355
 356 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
 357   memset(c,0,sizeof(codebook));
 358   c->c=s;
 359   c->entries=s->entries;
 360   c->dim=s->dim;
 361   c->valuelist=_book_unquantize(s);
 362   c->decode_tree=_make_decode_tree(c);
 363   if(c->decode_tree==NULL)goto err_out;
 364   return(0);
 365  err_out:
 366   vorbis_book_clear(c);
 367   return(-1);
 368 }
 369
 370 static float _dist(int el,float *ref, float *b,int step){
 371   int i;
 372   float acc=0.;
 373   for(i=0;i<el;i++){
 374     float val=(ref[i]-b[i*step]);
 375     acc+=val*val;
 376   }
 377   return(acc);
 378 }
 379
 380 #include <stdio.h>
 381 int _best(codebook *book, float *a, int step){
 382   encode_aux_nearestmatch *nt=book->c->nearest_tree;
 383   encode_aux_threshmatch *tt=book->c->thresh_tree;
 384   encode_aux_pigeonhole *pt=book->c->pigeon_tree;
 385   int dim=book->dim;
 386   int ptr=0,k,o;
 387   /*int savebest=-1;
 388     float saverr;*/
 389
 390   /* do we have a threshhold encode hint? */
 391   if(tt){
 392     int index=0;
 393     /* find the quant val of each scalar */
 394     for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
 395       int i;
 396       /* linear search the quant list for now; it's small and although
 397          with > 8 entries, it would be faster to bisect, this would be
 398          a misplaced optimization for now */
 399       for(i=0;i<tt->threshvals-1;i++)
 400         if(a[o]<tt->quantthresh[i])break;
 401
 402       index=(index*tt->quantvals)+tt->quantmap[i];
 403     }
 404     /* regular lattices are easy :-) */
 405     if(book->c->lengthlist[index]>0) /* is this unused?  If so, we'll
 406                                         use a decision tree after all
 407                                         and fall through*/
 408       return(index);
 409   }
 410
 411   /* do we have a pigeonhole encode hint? */
 412   if(pt){
 413     const static_codebook *c=book->c;
 414     int i,besti=-1;
 415     float best;
 416     int entry=0;
 417
 418     /* dealing with sequentialness is a pain in the ass */
 419     if(c->q_sequencep){
 420       int pv;
 421       long mul=1;
 422       float qlast=0;
 423       for(k=0,o=0;k<dim;k++,o+=step){
 424         pv=(int)((a[o]-qlast-pt->min)/pt->del);
 425         if(pv<0 || pv>=pt->mapentries)break;
 426         entry+=pt->pigeonmap[pv]*mul;
 427         mul*=pt->quantvals;
 428         qlast+=pv*pt->del+pt->min;
 429       }
 430     }else{
 431       for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
 432         int pv=(int)((a[o]-pt->min)/pt->del);
 433         if(pv<0 || pv>=pt->mapentries)break;
 434         entry=entry*pt->quantvals+pt->pigeonmap[pv];
 435       }
 436     }
 437
 438     /* must be within the pigeonholable range; if we quant outside (or
 439        in an entry that we define no list for), brute force it */
 440     if(k==dim && pt->fitlength[entry]){
 441       /* search the abbreviated list */
 442       long *list=pt->fitlist+pt->fitmap[entry];
 443       for(i=0;i<pt->fitlength[entry];i++){
 444         float this=_dist(dim,book->valuelist+list[i]*dim,a,step);
 445         if(besti==-1 || this<best){
 446           best=this;
 447           besti=list[i];
 448         }
 449       }
 450
 451       return(besti);
 452     }
 453   }
 454
 455   if(nt){
 456     /* optimized using the decision tree */
 457     while(1){
 458       float c=0.;
 459       float *p=book->valuelist+nt->p[ptr];
 460       float *q=book->valuelist+nt->q[ptr];
 461
 462       for(k=0,o=0;k<dim;k++,o+=step)
 463         c+=(p[k]-q[k])*(a[o]-(p[k]+q[k])*.5);
 464
 465       if(c>0.) /* in A */
 466         ptr= -nt->ptr0[ptr];
 467       else     /* in B */
 468         ptr= -nt->ptr1[ptr];
 469       if(ptr<=0)break;
 470     }
 471     return(-ptr);
 472   }
 473
 474   /* brute force it! */
 475   {
 476     const static_codebook *c=book->c;
 477     int i,besti=-1;
 478     float best;
 479     float *e=book->valuelist;
 480     for(i=0;i<book->entries;i++){
 481       if(c->lengthlist[i]>0){
 482         float this=_dist(dim,e,a,step);
 483         if(besti==-1 || this<best){
 484           best=this;
 485           besti=i;
 486         }
 487       }
 488       e+=dim;
 489     }
 490
 491     /*if(savebest!=-1 && savebest!=besti){
 492       fprintf(stderr,"brute force/pigeonhole disagreement:\n"
 493               "original:");
 494       for(i=0;i<dim*step;i+=step)fprintf(stderr,"%g,",a[i]);
 495       fprintf(stderr,"\n"
 496               "pigeonhole (entry %d, err %g):",savebest,saverr);
 497       for(i=0;i<dim;i++)fprintf(stderr,"%g,",
 498                                 (book->valuelist+savebest*dim)[i]);
 499       fprintf(stderr,"\n"
 500               "bruteforce (entry %d, err %g):",besti,best);
 501       for(i=0;i<dim;i++)fprintf(stderr,"%g,",
 502                                 (book->valuelist+besti*dim)[i]);
 503       fprintf(stderr,"\n");
 504       }*/
 505     return(besti);
 506   }
 507 }
 508
 509 /* returns the entry number and *modifies a* to the remainder value ********/
 510 int vorbis_book_besterror(codebook *book,float *a,int step,int addmul){
 511   int dim=book->dim,i,o;
 512   int best=_best(book,a,step);
 513   switch(addmul){
 514   case 0:
 515     for(i=0,o=0;i<dim;i++,o+=step)
 516       a[o]-=(book->valuelist+best*dim)[i];
 517     break;
 518   case 1:
 519     for(i=0,o=0;i<dim;i++,o+=step){
 520       float val=(book->valuelist+best*dim)[i];
 521       if(val==0){
 522         a[o]=0;
 523       }else{
 524         a[o]/=val;
 525       }
 526     }
 527     break;
 528   }
 529   return(best);
 530 }
 531
 532 long vorbis_book_codeword(codebook *book,int entry){
 533   return book->codelist[entry];
 534 }
 535
 536 long vorbis_book_codelen(codebook *book,int entry){
 537   return book->c->lengthlist[entry];
 538 }
 539
 540 #ifdef _V_SELFTEST
 541
 542 /* Unit tests of the dequantizer; this stuff will be OK
 543    cross-platform, I simply want to be sure that special mapping cases
 544    actually work properly; a bug could go unnoticed for a while */
 545
 546 #include <stdio.h>
 547
 548 /* cases:
 549
 550    no mapping
 551    full, explicit mapping
 552    algorithmic mapping
 553
 554    nonsequential
 555    sequential
 556 */
 557
 558 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
 559 static long partial_quantlist1[]={0,7,2};
 560
 561 /* no mapping */
 562 static_codebook test1={
 563   4,16,
 564   NULL,
 565   0,
 566   0,0,0,0,
 567   NULL,
 568   NULL,NULL
 569 };
 570 static float *test1_result=NULL;
 571
 572 /* linear, full mapping, nonsequential */
 573 static_codebook test2={
 574   4,3,
 575   NULL,
 576   2,
 577   -533200896,1611661312,4,0,
 578   full_quantlist1,
 579   NULL,NULL
 580 };
 581 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
 582
 583 /* linear, full mapping, sequential */
 584 static_codebook test3={
 585   4,3,
 586   NULL,
 587   2,
 588   -533200896,1611661312,4,1,
 589   full_quantlist1,
 590   NULL,NULL
 591 };
 592 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
 593
 594 /* linear, algorithmic mapping, nonsequential */
 595 static_codebook test4={
 596   3,27,
 597   NULL,
 598   1,
 599   -533200896,1611661312,4,0,
 600   partial_quantlist1,
 601   NULL,NULL
 602 };
 603 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
 604                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
 605                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
 606                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
 607                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
 608                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
 609                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
 610                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
 611                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
 612
 613 /* linear, algorithmic mapping, sequential */
 614 static_codebook test5={
 615   3,27,
 616   NULL,
 617   1,
 618   -533200896,1611661312,4,1,
 619   partial_quantlist1,
 620   NULL,NULL
 621 };
 622 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
 623                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
 624                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
 625                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
 626                               -3, 1, 5, 4, 8,12, -1, 3, 7,
 627                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
 628                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
 629                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
 630                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
 631
 632 void run_test(static_codebook *b,float *comp){
 633   float *out=_book_unquantize(b);
 634   int i;
 635
 636   if(comp){
 637     if(!out){
 638       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
 639       exit(1);
 640     }
 641
 642     for(i=0;i<b->entries*b->dim;i++)
 643       if(fabs(out[i]-comp[i])>.0001){
 644         fprintf(stderr,"disagreement in unquantized and reference data:\n"
 645                 "position %d, %g != %g\n",i,out[i],comp[i]);
 646         exit(1);
 647       }
 648
 649   }else{
 650     if(out){
 651       fprintf(stderr,"_book_unquantize returned a value array: \n"
 652               " correct result should have been NULL\n");
 653       exit(1);
 654     }
 655   }
 656 }
 657
 658 int main(){
 659   /* run the nine dequant tests, and compare to the hand-rolled results */
 660   fprintf(stderr,"Dequant test 1... ");
 661   run_test(&test1,test1_result);
 662   fprintf(stderr,"OK\nDequant test 2... ");
 663   run_test(&test2,test2_result);
 664   fprintf(stderr,"OK\nDequant test 3... ");
 665   run_test(&test3,test3_result);
 666   fprintf(stderr,"OK\nDequant test 4... ");
 667   run_test(&test4,test4_result);
 668   fprintf(stderr,"OK\nDequant test 5... ");
 669   run_test(&test5,test5_result);
 670   fprintf(stderr,"OK\n\n");
 671
 672   return(0);
 673 }
 674
 675 #endif