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