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.7 2000/07/19 18:10:46 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdlib.h>
  20 #include <math.h>
  21 #include <string.h>
  22 #include "os.h"
  23 #include "vorbis/codec.h"
  24 #include "vorbis/codebook.h"
  25 #include "bitwise.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(double 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 double _float32_unpack(long val){
  64   double mant=val&0x1fffff;
  65   double sign=val&0x80000000;
  66   double 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;
 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   return(t);
 177 }
 178
 179 /* there might be a straightforward one-line way to do the below
 180    that's portable and totally safe against roundoff, but I haven't
 181    thought of it.  Therefore, we opt on the side of caution */
 182 long _book_maptype1_quantvals(const static_codebook *b){
 183   long vals=floor(pow(b->entries,1./b->dim));
 184
 185   /* the above *should* be reliable, but we'll not assume that FP is
 186      ever reliable when bitstream sync is at stake; verify via integer
 187      means that vals really is the greatest value of dim for which
 188      vals^b->bim <= b->entries */
 189   /* treat the above as an initial guess */
 190   while(1){
 191     long acc=1;
 192     long acc1=1;
 193     int i;
 194     for(i=0;i<b->dim;i++){
 195       acc*=vals;
 196       acc1*=vals+1;
 197     }
 198     if(acc<=b->entries && acc1>b->entries){
 199       return(vals);
 200     }else{
 201       if(acc>b->entries){
 202         vals--;
 203       }else{
 204         vals++;
 205       }
 206     }
 207   }
 208 }
 209
 210 /* unpack the quantized list of values for encode/decode ***********/
 211 /* we need to deal with two map types: in map type 1, the values are
 212    generated algorithmically (each column of the vector counts through
 213    the values in the quant vector). in map type 2, all the values came
 214    in in an explicit list.  Both value lists must be unpacked */
 215 double *_book_unquantize(const static_codebook *b){
 216   long j,k;
 217   if(b->maptype==1 || b->maptype==2){
 218     int quantvals;
 219     double mindel=_float32_unpack(b->q_min);
 220     double delta=_float32_unpack(b->q_delta);
 221     double *r=calloc(b->entries*b->dim,sizeof(double));
 222
 223     /* maptype 1 and 2 both use a quantized value vector, but
 224        different sizes */
 225     switch(b->maptype){
 226     case 1:
 227       /* most of the time, entries%dimensions == 0, but we need to be
 228          well defined.  We define that the possible vales at each
 229          scalar is values == entries/dim.  If entries%dim != 0, we'll
 230          have 'too few' values (values*dim<entries), which means that
 231          we'll have 'left over' entries; left over entries use zeroed
 232          values (and are wasted).  So don't generate codebooks like
 233          that */
 234       quantvals=_book_maptype1_quantvals(b);
 235       for(j=0;j<b->entries;j++){
 236         double last=0.;
 237         int indexdiv=1;
 238         for(k=0;k<b->dim;k++){
 239           int index= (j/indexdiv)%quantvals;
 240           double val=b->quantlist[index];
 241           val=fabs(val)*delta+mindel+last;
 242           if(b->q_sequencep)last=val;
 243           r[j*b->dim+k]=val;
 244           indexdiv*=quantvals;
 245         }
 246       }
 247       break;
 248     case 2:
 249       for(j=0;j<b->entries;j++){
 250         double last=0.;
 251         for(k=0;k<b->dim;k++){
 252           double val=b->quantlist[j*b->dim+k];
 253           val=fabs(val)*delta+mindel+last;
 254           if(b->q_sequencep)last=val;
 255           r[j*b->dim+k]=val;
 256         }
 257       }
 258     }
 259     return(r);
 260   }
 261   return(NULL);
 262 }
 263
 264 void vorbis_staticbook_clear(static_codebook *b){
 265   if(b->quantlist)free(b->quantlist);
 266   if(b->lengthlist)free(b->lengthlist);
 267   if(b->nearest_tree){
 268     free(b->nearest_tree->ptr0);
 269     free(b->nearest_tree->ptr1);
 270     free(b->nearest_tree->p);
 271     free(b->nearest_tree->q);
 272     memset(b->nearest_tree,0,sizeof(encode_aux_nearestmatch));
 273     free(b->nearest_tree);
 274   }
 275   if(b->thresh_tree){
 276     free(b->thresh_tree->quantthresh);
 277     free(b->thresh_tree->quantmap);
 278     memset(b->thresh_tree,0,sizeof(encode_aux_threshmatch));
 279     free(b->thresh_tree);
 280   }
 281   memset(b,0,sizeof(static_codebook));
 282 }
 283
 284 void vorbis_book_clear(codebook *b){
 285   /* static book is not cleared; we're likely called on the lookup and
 286      the static codebook belongs to the info struct */
 287   if(b->decode_tree){
 288     free(b->decode_tree->ptr0);
 289     free(b->decode_tree->ptr1);
 290     memset(b->decode_tree,0,sizeof(decode_aux));
 291     free(b->decode_tree);
 292   }
 293   if(b->valuelist)free(b->valuelist);
 294   if(b->codelist)free(b->codelist);
 295   memset(b,0,sizeof(codebook));
 296 }
 297
 298 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
 299   memset(c,0,sizeof(codebook));
 300   c->c=s;
 301   c->entries=s->entries;
 302   c->dim=s->dim;
 303   c->codelist=_make_words(s->lengthlist,s->entries);
 304   c->valuelist=_book_unquantize(s);
 305   return(0);
 306 }
 307
 308 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
 309   memset(c,0,sizeof(codebook));
 310   c->c=s;
 311   c->entries=s->entries;
 312   c->dim=s->dim;
 313   c->valuelist=_book_unquantize(s);
 314   c->decode_tree=_make_decode_tree(c);
 315   if(c->decode_tree==NULL)goto err_out;
 316   return(0);
 317  err_out:
 318   vorbis_book_clear(c);
 319   return(-1);
 320 }
 321
 322 static double _dist(int el,double *ref, double *b,int step){
 323   int i;
 324   double acc=0.;
 325   for(i=0;i<el;i++){
 326     double val=(ref[i]-b[i*step]);
 327     acc+=val*val;
 328   }
 329   return(acc);
 330 }
 331
 332 #include <stdio.h>
 333 int _best(codebook *book, double *a, int step){
 334   encode_aux_nearestmatch *nt=book->c->nearest_tree;
 335   encode_aux_threshmatch *tt=book->c->thresh_tree;
 336   encode_aux_pigeonhole *pt=book->c->pigeon_tree;
 337   int dim=book->dim;
 338   int ptr=0,k,o;
 339   int savebest=-1;
 340   double saverr;
 341
 342   /* do we have a threshhold encode hint? */
 343   if(tt){
 344     int index=0;
 345     /* find the quant val of each scalar */
 346     for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
 347       int i;
 348       /* linear search the quant list for now; it's small and although
 349          with > 8 entries, it would be faster to bisect, this would be
 350          a misplaced optimization for now */
 351       for(i=0;i<tt->threshvals-1;i++)
 352         if(a[o]<tt->quantthresh[i])break;
 353
 354       index=(index*tt->quantvals)+tt->quantmap[i];
 355     }
 356     /* regular lattices are easy :-) */
 357     if(book->c->lengthlist[index]>0) /* is this unused?  If so, we'll
 358                                         use a decision tree after all
 359                                         and fall through*/
 360       return(index);
 361   }
 362
 363   /* do we have a pigeonhole encode hint? */
 364   if(pt){
 365     const static_codebook *c=book->c;
 366     int i,besti=-1;
 367     double best;
 368     int entry=0;
 369
 370     /* dealing with sequentialness is a pain in the ass */
 371     if(c->q_sequencep){
 372       int pv;
 373       long mul=1;
 374       double qlast=0;
 375       for(k=0,o=0;k<dim;k++,o+=step){
 376         pv=(int)((a[o]-qlast-pt->min)/pt->del);
 377         if(pv<0 || pv>=pt->mapentries)break;
 378         entry+=pt->pigeonmap[pv]*mul;
 379         mul*=pt->quantvals;
 380         qlast+=pv*pt->del+pt->min;
 381       }
 382     }else{
 383       for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
 384         int pv=(int)((a[o]-pt->min)/pt->del);
 385         if(pv<0 || pv>=pt->mapentries)break;
 386         entry=entry*pt->quantvals+pt->pigeonmap[pv];
 387       }
 388     }
 389
 390     /* must be within the pigeonholable range; if we quant outside (or
 391        in an entry that we define no list for), brute force it */
 392     if(k==dim && pt->fitlength[entry]){
 393       /* search the abbreviated list */
 394       long *list=pt->fitlist+pt->fitmap[entry];
 395       for(i=0;i<pt->fitlength[entry];i++){
 396         double this=_dist(dim,book->valuelist+list[i]*dim,a,step);
 397         if(besti==-1 || this<best){
 398           best=this;
 399           besti=list[i];
 400         }
 401       }
 402
 403       return(besti);
 404     }
 405   }
 406
 407   if(nt){
 408     /* optimized using the decision tree */
 409     while(1){
 410       double c=0.;
 411       double *p=book->valuelist+nt->p[ptr];
 412       double *q=book->valuelist+nt->q[ptr];
 413
 414       for(k=0,o=0;k<dim;k++,o+=step)
 415         c+=(p[k]-q[k])*(a[o]-(p[k]+q[k])*.5);
 416
 417       if(c>0.) /* in A */
 418         ptr= -nt->ptr0[ptr];
 419       else     /* in B */
 420         ptr= -nt->ptr1[ptr];
 421       if(ptr<=0)break;
 422     }
 423     return(-ptr);
 424   }
 425
 426   /* brute force it! */
 427   {
 428     const static_codebook *c=book->c;
 429     int i,besti=-1;
 430     double best;
 431     double *e=book->valuelist;
 432     for(i=0;i<book->entries;i++){
 433       if(c->lengthlist[i]>0){
 434         double this=_dist(dim,e,a,step);
 435         if(besti==-1 || this<best){
 436           best=this;
 437           besti=i;
 438         }
 439       }
 440       e+=dim;
 441     }
 442
 443     /*if(savebest!=-1 && savebest!=besti){
 444       fprintf(stderr,"brute force/pigeonhole disagreement:\n"
 445               "original:");
 446       for(i=0;i<dim*step;i+=step)fprintf(stderr,"%g,",a[i]);
 447       fprintf(stderr,"\n"
 448               "pigeonhole (entry %d, err %g):",savebest,saverr);
 449       for(i=0;i<dim;i++)fprintf(stderr,"%g,",
 450                                 (book->valuelist+savebest*dim)[i]);
 451       fprintf(stderr,"\n"
 452               "bruteforce (entry %d, err %g):",besti,best);
 453       for(i=0;i<dim;i++)fprintf(stderr,"%g,",
 454                                 (book->valuelist+besti*dim)[i]);
 455       fprintf(stderr,"\n");
 456       }*/
 457     return(besti);
 458   }
 459 }
 460
 461 /* returns the entry number and *modifies a* to the remainder value ********/
 462 int vorbis_book_besterror(codebook *book,double *a,int step,int addmul){
 463   int dim=book->dim,i,o;
 464   int best=_best(book,a,step);
 465   switch(addmul){
 466   case 0:
 467     for(i=0,o=0;i<dim;i++,o+=step)
 468       a[o]-=(book->valuelist+best*dim)[i];
 469     break;
 470   case 1:
 471     for(i=0,o=0;i<dim;i++,o+=step){
 472       double val=(book->valuelist+best*dim)[i];
 473       if(val==0){
 474         a[o]=0;
 475       }else{
 476         a[o]/=val;
 477       }
 478     }
 479     break;
 480   }
 481   return(best);
 482 }
 483
 484 long vorbis_book_codeword(codebook *book,int entry){
 485   return book->codelist[entry];
 486 }
 487
 488 long vorbis_book_codelen(codebook *book,int entry){
 489   return book->c->lengthlist[entry];
 490 }
 491
 492 #ifdef _V_SELFTEST
 493
 494 /* Unit tests of the dequantizer; this stuff will be OK
 495    cross-platform, I simply want to be sure that special mapping cases
 496    actually work properly; a bug could go unnoticed for a while */
 497
 498 #include <stdio.h>
 499
 500 /* cases:
 501
 502    no mapping
 503    full, explicit mapping
 504    algorithmic mapping
 505
 506    nonsequential
 507    sequential
 508 */
 509
 510 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
 511 static long partial_quantlist1[]={0,7,2};
 512
 513 /* no mapping */
 514 static_codebook test1={
 515   4,16,
 516   NULL,
 517   0,
 518   0,0,0,0,
 519   NULL,
 520   NULL,NULL
 521 };
 522 static double *test1_result=NULL;
 523
 524 /* linear, full mapping, nonsequential */
 525 static_codebook test2={
 526   4,3,
 527   NULL,
 528   2,
 529   -533200896,1611661312,4,0,
 530   full_quantlist1,
 531   NULL,NULL
 532 };
 533 static double test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
 534
 535 /* linear, full mapping, sequential */
 536 static_codebook test3={
 537   4,3,
 538   NULL,
 539   2,
 540   -533200896,1611661312,4,1,
 541   full_quantlist1,
 542   NULL,NULL
 543 };
 544 static double test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
 545
 546 /* linear, algorithmic mapping, nonsequential */
 547 static_codebook test4={
 548   3,27,
 549   NULL,
 550   1,
 551   -533200896,1611661312,4,0,
 552   partial_quantlist1,
 553   NULL,NULL
 554 };
 555 static double test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
 556                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
 557                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
 558                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
 559                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
 560                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
 561                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
 562                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
 563                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
 564
 565 /* linear, algorithmic mapping, sequential */
 566 static_codebook test5={
 567   3,27,
 568   NULL,
 569   1,
 570   -533200896,1611661312,4,1,
 571   partial_quantlist1,
 572   NULL,NULL
 573 };
 574 static double test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
 575                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
 576                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
 577                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
 578                               -3, 1, 5, 4, 8,12, -1, 3, 7,
 579                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
 580                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
 581                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
 582                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
 583
 584 void run_test(static_codebook *b,double *comp){
 585   double *out=_book_unquantize(b);
 586   int i;
 587
 588   if(comp){
 589     if(!out){
 590       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
 591       exit(1);
 592     }
 593
 594     for(i=0;i<b->entries*b->dim;i++)
 595       if(fabs(out[i]-comp[i])>.0001){
 596         fprintf(stderr,"disagreement in unquantized and reference data:\n"
 597                 "position %d, %g != %g\n",i,out[i],comp[i]);
 598         exit(1);
 599       }
 600
 601   }else{
 602     if(out){
 603       fprintf(stderr,"_book_unquantize returned a value array: \n"
 604               " correct result should have been NULL\n");
 605       exit(1);
 606     }
 607   }
 608 }
 609
 610 int main(){
 611   /* run the nine dequant tests, and compare to the hand-rolled results */
 612   fprintf(stderr,"Dequant test 1... ");
 613   run_test(&test1,test1_result);
 614   fprintf(stderr,"OK\nDequant test 2... ");
 615   run_test(&test2,test2_result);
 616   fprintf(stderr,"OK\nDequant test 3... ");
 617   run_test(&test3,test3_result);
 618   fprintf(stderr,"OK\nDequant test 4... ");
 619   run_test(&test4,test4_result);
 620   fprintf(stderr,"OK\nDequant test 5... ");
 621   run_test(&test5,test5_result);
 622   fprintf(stderr,"OK\n\n");
 623
 624   return(0);
 625 }
 626
 627 #endif