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.3 2000/05/17 22:34:22 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdlib.h>
  20 #include <math.h>
  21 #include <string.h>
  22 #include "vorbis/codec.h"
  23 #include "vorbis/codebook.h"
  24 #include "bitwise.h"
  25 #include "scales.h"
  26 #include "sharedbook.h"
  27 #include "os.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 int _best(codebook *book, double *a, int step){
 323   encode_aux_nearestmatch *nt=book->c->nearest_tree;
 324   encode_aux_threshmatch *tt=book->c->thresh_tree;
 325   int dim=book->dim;
 326   int ptr=0,k,o;
 327
 328   /* we assume for now that a thresh tree is the only other possibility */
 329   if(tt){
 330     int index=0;
 331     /* find the quant val of each scalar */
 332     for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
 333       int i;
 334       /* linear search the quant list for now; it's small and although
 335          with > 8 entries, it would be faster to bisect, this would be
 336          a misplaced optimization for now */
 337       for(i=0;i<tt->threshvals-1;i++)
 338         if(a[o]<tt->quantthresh[i])break;
 339
 340       index=(index*tt->quantvals)+tt->quantmap[i];
 341     }
 342     /* regular lattices are easy :-) */
 343     if(book->c->lengthlist[index]>0) /* is this unused?  If so, we'll
 344                                         use a decision tree after all
 345                                         and fall through*/
 346       return(index);
 347   }
 348
 349   if(nt){
 350     /* optimized using the decision tree */
 351     while(1){
 352       double c=0.;
 353       double *p=book->valuelist+nt->p[ptr];
 354       double *q=book->valuelist+nt->q[ptr];
 355
 356       for(k=0,o=0;k<dim;k++,o+=step)
 357         c+=(p[k]-q[k])*(a[o]-(p[k]+q[k])*.5);
 358
 359       if(c>0.) /* in A */
 360         ptr= -nt->ptr0[ptr];
 361       else     /* in B */
 362         ptr= -nt->ptr1[ptr];
 363       if(ptr<=0)break;
 364     }
 365     return(-ptr);
 366   }
 367
 368   return(-1);
 369 }
 370
 371 static double _dist(int el,double *a, double *b){
 372   int i;
 373   double acc=0.;
 374   for(i=0;i<el;i++){
 375     double val=(a[i]-b[i]);
 376     acc+=val*val;
 377   }
 378   return(acc);
 379 }
 380
 381 /* returns the entry number and *modifies a* to the remainder value ********/
 382 int vorbis_book_besterror(codebook *book,double *a,int step,int addmul){
 383   int dim=book->dim,i,o;
 384   int best=_best(book,a,step);
 385   switch(addmul){
 386   case 0:
 387     for(i=0,o=0;i<dim;i++,o+=step)
 388       a[o]-=(book->valuelist+best*dim)[i];
 389     break;
 390   case 1:
 391     for(i=0,o=0;i<dim;i++,o+=step){
 392       double val=(book->valuelist+best*dim)[i];
 393       if(val==0){
 394         a[o]=0;
 395       }else{
 396         a[o]/=val;
 397       }
 398     }
 399     break;
 400   }
 401   return(best);
 402 }
 403
 404 long vorbis_book_codeword(codebook *book,int entry){
 405   return book->codelist[entry];
 406 }
 407
 408 long vorbis_book_codelen(codebook *book,int entry){
 409   return book->c->lengthlist[entry];
 410 }
 411
 412 #ifdef _V_SELFTEST
 413
 414 /* Unit tests of the dequantizer; this stuff will be OK
 415    cross-platform, I simply want to be sure that special mapping cases
 416    actually work properly; a bug could go unnoticed for a while */
 417
 418 #include <stdio.h>
 419
 420 /* cases:
 421
 422    no mapping
 423    full, explicit mapping
 424    algorithmic mapping
 425
 426    nonsequential
 427    sequential
 428 */
 429
 430 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
 431 static long partial_quantlist1[]={0,7,2};
 432
 433 /* no mapping */
 434 static_codebook test1={
 435   4,16,
 436   NULL,
 437   0,
 438   0,0,0,0,
 439   NULL,
 440   NULL,NULL
 441 };
 442 static double *test1_result=NULL;
 443
 444 /* linear, full mapping, nonsequential */
 445 static_codebook test2={
 446   4,3,
 447   NULL,
 448   2,
 449   -533200896,1611661312,4,0,
 450   full_quantlist1,
 451   NULL,NULL
 452 };
 453 static double test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
 454
 455 /* linear, full mapping, sequential */
 456 static_codebook test3={
 457   4,3,
 458   NULL,
 459   2,
 460   -533200896,1611661312,4,1,
 461   full_quantlist1,
 462   NULL,NULL
 463 };
 464 static double test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
 465
 466 /* linear, algorithmic mapping, nonsequential */
 467 static_codebook test4={
 468   3,27,
 469   NULL,
 470   1,
 471   -533200896,1611661312,4,0,
 472   partial_quantlist1,
 473   NULL,NULL
 474 };
 475 static double test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
 476                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
 477                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
 478                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
 479                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
 480                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
 481                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
 482                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
 483                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
 484
 485 /* linear, algorithmic mapping, sequential */
 486 static_codebook test5={
 487   3,27,
 488   NULL,
 489   1,
 490   -533200896,1611661312,4,1,
 491   partial_quantlist1,
 492   NULL,NULL
 493 };
 494 static double test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
 495                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
 496                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
 497                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
 498                               -3, 1, 5, 4, 8,12, -1, 3, 7,
 499                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
 500                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
 501                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
 502                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
 503
 504 void run_test(static_codebook *b,double *comp){
 505   double *out=_book_unquantize(b);
 506   int i;
 507
 508   if(comp){
 509     if(!out){
 510       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
 511       exit(1);
 512     }
 513
 514     for(i=0;i<b->entries*b->dim;i++)
 515       if(fabs(out[i]-comp[i])>.0001){
 516         fprintf(stderr,"disagreement in unquantized and reference data:\n"
 517                 "position %d, %g != %g\n",i,out[i],comp[i]);
 518         exit(1);
 519       }
 520
 521   }else{
 522     if(out){
 523       fprintf(stderr,"_book_unquantize returned a value array: \n"
 524               " correct result should have been NULL\n");
 525       exit(1);
 526     }
 527   }
 528 }
 529
 530 int main(){
 531   /* run the nine dequant tests, and compare to the hand-rolled results */
 532   fprintf(stderr,"Dequant test 1... ");
 533   run_test(&test1,test1_result);
 534   fprintf(stderr,"OK\nDequant test 2... ");
 535   run_test(&test2,test2_result);
 536   fprintf(stderr,"OK\nDequant test 3... ");
 537   run_test(&test3,test3_result);
 538   fprintf(stderr,"OK\nDequant test 4... ");
 539   run_test(&test4,test4_result);
 540   fprintf(stderr,"OK\nDequant test 5... ");
 541   run_test(&test5,test5_result);
 542   fprintf(stderr,"OK\n\n");
 543
 544   return(0);
 545 }
 546
 547 #endif