lib/sharedbook.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
   5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
   6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
   7  *                                                                  *
   8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
   9  * by the Xiph.Org Foundation http://www.xiph.org/                  *
  10  *                                                                  *
  11  ********************************************************************
  12
  13  function: basic shared codebook operations
  14  last mod: $Id$
  15
  16  ********************************************************************/
  17
  18 #include <stdlib.h>
  19 #include <math.h>
  20 #include <string.h>
  21 #include <ogg/ogg.h>
  22 #include "os.h"
  23 #include "misc.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.f)+.001); //+epsilon
  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   double mant=val&0x1fffff;
  64   int    sign=val&0x80000000;
  65   long   exp =(val&0x7fe00000L)>>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 ogg_uint32_t *_make_words(long *l,long n,long sparsecount){
  74   long i,j,count=0;
  75   ogg_uint32_t marker[33];
  76   ogg_uint32_t *r=_ogg_malloc((sparsecount?sparsecount:n)*sizeof(*r));
  77   memset(marker,0,sizeof(marker));
  78
  79   for(i=0;i<n;i++){
  80     long length=l[i];
  81     if(length>0){
  82       ogg_uint32_t 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         _ogg_free(r);
  93         return(NULL);
  94       }
  95       r[count++]=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     }else
 125       if(sparsecount==0)count++;
 126   }
 127
 128   /* sanity check the huffman tree; an underpopulated tree must be
 129      rejected. The only exception is the one-node pseudo-nil tree,
 130      which appears to be underpopulated because the tree doesn't
 131      really exist; there's only one possible 'codeword' or zero bits,
 132      but the above tree-gen code doesn't mark that. */
 133   if(sparsecount != 1){
 134     for(i=1;i<33;i++)
 135       if(marker[i] & (0xffffffffUL>>(32-i))){
 136         _ogg_free(r);
 137         return(NULL);
 138       }
 139   }
 140
 141   /* bitreverse the words because our bitwise packer/unpacker is LSb
 142      endian */
 143   for(i=0,count=0;i<n;i++){
 144     ogg_uint32_t temp=0;
 145     for(j=0;j<l[i];j++){
 146       temp<<=1;
 147       temp|=(r[count]>>j)&1;
 148     }
 149
 150     if(sparsecount){
 151       if(l[i])
 152         r[count++]=temp;
 153     }else
 154       r[count++]=temp;
 155   }
 156
 157   return(r);
 158 }
 159
 160 /* there might be a straightforward one-line way to do the below
 161    that's portable and totally safe against roundoff, but I haven't
 162    thought of it.  Therefore, we opt on the side of caution */
 163 long _book_maptype1_quantvals(const static_codebook *b){
 164   long vals=floor(pow((float)b->entries,1.f/b->dim));
 165
 166   /* the above *should* be reliable, but we'll not assume that FP is
 167      ever reliable when bitstream sync is at stake; verify via integer
 168      means that vals really is the greatest value of dim for which
 169      vals^b->bim <= b->entries */
 170   /* treat the above as an initial guess */
 171   while(1){
 172     long acc=1;
 173     long acc1=1;
 174     int i;
 175     for(i=0;i<b->dim;i++){
 176       acc*=vals;
 177       acc1*=vals+1;
 178     }
 179     if(acc<=b->entries && acc1>b->entries){
 180       return(vals);
 181     }else{
 182       if(acc>b->entries){
 183         vals--;
 184       }else{
 185         vals++;
 186       }
 187     }
 188   }
 189 }
 190
 191 /* unpack the quantized list of values for encode/decode ***********/
 192 /* we need to deal with two map types: in map type 1, the values are
 193    generated algorithmically (each column of the vector counts through
 194    the values in the quant vector). in map type 2, all the values came
 195    in in an explicit list.  Both value lists must be unpacked */
 196 float *_book_unquantize(const static_codebook *b,int n,int *sparsemap){
 197   long j,k,count=0;
 198   if(b->maptype==1 || b->maptype==2){
 199     int quantvals;
 200     float mindel=_float32_unpack(b->q_min);
 201     float delta=_float32_unpack(b->q_delta);
 202     float *r=_ogg_calloc(n*b->dim,sizeof(*r));
 203
 204     /* maptype 1 and 2 both use a quantized value vector, but
 205        different sizes */
 206     switch(b->maptype){
 207     case 1:
 208       /* most of the time, entries%dimensions == 0, but we need to be
 209          well defined.  We define that the possible vales at each
 210          scalar is values == entries/dim.  If entries%dim != 0, we'll
 211          have 'too few' values (values*dim<entries), which means that
 212          we'll have 'left over' entries; left over entries use zeroed
 213          values (and are wasted).  So don't generate codebooks like
 214          that */
 215       quantvals=_book_maptype1_quantvals(b);
 216       for(j=0;j<b->entries;j++){
 217         if((sparsemap && b->lengthlist[j]) || !sparsemap){
 218           float last=0.f;
 219           int indexdiv=1;
 220           for(k=0;k<b->dim;k++){
 221             int index= (j/indexdiv)%quantvals;
 222             float val=b->quantlist[index];
 223             val=fabs(val)*delta+mindel+last;
 224             if(b->q_sequencep)last=val;
 225             if(sparsemap)
 226               r[sparsemap[count]*b->dim+k]=val;
 227             else
 228               r[count*b->dim+k]=val;
 229             indexdiv*=quantvals;
 230           }
 231           count++;
 232         }
 233
 234       }
 235       break;
 236     case 2:
 237       for(j=0;j<b->entries;j++){
 238         if((sparsemap && b->lengthlist[j]) || !sparsemap){
 239           float last=0.f;
 240
 241           for(k=0;k<b->dim;k++){
 242             float val=b->quantlist[j*b->dim+k];
 243             val=fabs(val)*delta+mindel+last;
 244             if(b->q_sequencep)last=val;
 245             if(sparsemap)
 246               r[sparsemap[count]*b->dim+k]=val;
 247             else
 248               r[count*b->dim+k]=val;
 249           }
 250           count++;
 251         }
 252       }
 253       break;
 254     }
 255
 256     return(r);
 257   }
 258   return(NULL);
 259 }
 260
 261 void vorbis_staticbook_destroy(static_codebook *b){
 262   if(b->allocedp){
 263     if(b->quantlist)_ogg_free(b->quantlist);
 264     if(b->lengthlist)_ogg_free(b->lengthlist);
 265     memset(b,0,sizeof(*b));
 266     _ogg_free(b);
 267   } /* otherwise, it is in static memory */
 268 }
 269
 270 void vorbis_book_clear(codebook *b){
 271   /* static book is not cleared; we're likely called on the lookup and
 272      the static codebook belongs to the info struct */
 273   if(b->valuelist)_ogg_free(b->valuelist);
 274   if(b->codelist)_ogg_free(b->codelist);
 275
 276   if(b->dec_index)_ogg_free(b->dec_index);
 277   if(b->dec_codelengths)_ogg_free(b->dec_codelengths);
 278   if(b->dec_firsttable)_ogg_free(b->dec_firsttable);
 279
 280   memset(b,0,sizeof(*b));
 281 }
 282
 283 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
 284
 285   memset(c,0,sizeof(*c));
 286   c->c=s;
 287   c->entries=s->entries;
 288   c->used_entries=s->entries;
 289   c->dim=s->dim;
 290   c->codelist=_make_words(s->lengthlist,s->entries,0);
 291   //c->valuelist=_book_unquantize(s,s->entries,NULL);
 292   c->quantvals=_book_maptype1_quantvals(s);
 293   c->minval=(int)rint(_float32_unpack(s->q_min));
 294   c->delta=(int)rint(_float32_unpack(s->q_delta));
 295
 296   return(0);
 297 }
 298
 299 static ogg_uint32_t bitreverse(ogg_uint32_t x){
 300   x=    ((x>>16)&0x0000ffffUL) | ((x<<16)&0xffff0000UL);
 301   x=    ((x>> 8)&0x00ff00ffUL) | ((x<< 8)&0xff00ff00UL);
 302   x=    ((x>> 4)&0x0f0f0f0fUL) | ((x<< 4)&0xf0f0f0f0UL);
 303   x=    ((x>> 2)&0x33333333UL) | ((x<< 2)&0xccccccccUL);
 304   return((x>> 1)&0x55555555UL) | ((x<< 1)&0xaaaaaaaaUL);
 305 }
 306
 307 static int sort32a(const void *a,const void *b){
 308   return ( **(ogg_uint32_t **)a>**(ogg_uint32_t **)b)-
 309     ( **(ogg_uint32_t **)a<**(ogg_uint32_t **)b);
 310 }
 311
 312 /* decode codebook arrangement is more heavily optimized than encode */
 313 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
 314   int i,j,n=0,tabn;
 315   int *sortindex;
 316   memset(c,0,sizeof(*c));
 317
 318   /* count actually used entries */
 319   for(i=0;i<s->entries;i++)
 320     if(s->lengthlist[i]>0)
 321       n++;
 322
 323   c->entries=s->entries;
 324   c->used_entries=n;
 325   c->dim=s->dim;
 326
 327   if(n>0){
 328
 329     /* two different remappings go on here.
 330
 331     First, we collapse the likely sparse codebook down only to
 332     actually represented values/words.  This collapsing needs to be
 333     indexed as map-valueless books are used to encode original entry
 334     positions as integers.
 335
 336     Second, we reorder all vectors, including the entry index above,
 337     by sorted bitreversed codeword to allow treeless decode. */
 338
 339     /* perform sort */
 340     ogg_uint32_t *codes=_make_words(s->lengthlist,s->entries,c->used_entries);
 341     ogg_uint32_t **codep=alloca(sizeof(*codep)*n);
 342
 343     if(codes==NULL)goto err_out;
 344
 345     for(i=0;i<n;i++){
 346       codes[i]=bitreverse(codes[i]);
 347       codep[i]=codes+i;
 348     }
 349
 350     qsort(codep,n,sizeof(*codep),sort32a);
 351
 352     sortindex=alloca(n*sizeof(*sortindex));
 353     c->codelist=_ogg_malloc(n*sizeof(*c->codelist));
 354     /* the index is a reverse index */
 355     for(i=0;i<n;i++){
 356       int position=codep[i]-codes;
 357       sortindex[position]=i;
 358     }
 359
 360     for(i=0;i<n;i++)
 361       c->codelist[sortindex[i]]=codes[i];
 362     _ogg_free(codes);
 363
 364
 365     c->valuelist=_book_unquantize(s,n,sortindex);
 366     c->dec_index=_ogg_malloc(n*sizeof(*c->dec_index));
 367
 368     for(n=0,i=0;i<s->entries;i++)
 369       if(s->lengthlist[i]>0)
 370         c->dec_index[sortindex[n++]]=i;
 371
 372     c->dec_codelengths=_ogg_malloc(n*sizeof(*c->dec_codelengths));
 373     for(n=0,i=0;i<s->entries;i++)
 374       if(s->lengthlist[i]>0)
 375         c->dec_codelengths[sortindex[n++]]=s->lengthlist[i];
 376
 377     c->dec_firsttablen=_ilog(c->used_entries)-4; /* this is magic */
 378     if(c->dec_firsttablen<5)c->dec_firsttablen=5;
 379     if(c->dec_firsttablen>8)c->dec_firsttablen=8;
 380
 381     tabn=1<<c->dec_firsttablen;
 382     c->dec_firsttable=_ogg_calloc(tabn,sizeof(*c->dec_firsttable));
 383     c->dec_maxlength=0;
 384
 385     for(i=0;i<n;i++){
 386       if(c->dec_maxlength<c->dec_codelengths[i])
 387         c->dec_maxlength=c->dec_codelengths[i];
 388       if(c->dec_codelengths[i]<=c->dec_firsttablen){
 389         ogg_uint32_t orig=bitreverse(c->codelist[i]);
 390         for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
 391           c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
 392       }
 393     }
 394
 395     /* now fill in 'unused' entries in the firsttable with hi/lo search
 396        hints for the non-direct-hits */
 397     {
 398       ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
 399       long lo=0,hi=0;
 400
 401       for(i=0;i<tabn;i++){
 402         ogg_uint32_t word=i<<(32-c->dec_firsttablen);
 403         if(c->dec_firsttable[bitreverse(word)]==0){
 404           while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
 405           while(    hi<n && word>=(c->codelist[hi]&mask))hi++;
 406
 407           /* we only actually have 15 bits per hint to play with here.
 408              In order to overflow gracefully (nothing breaks, efficiency
 409              just drops), encode as the difference from the extremes. */
 410           {
 411             unsigned long loval=lo;
 412             unsigned long hival=n-hi;
 413
 414             if(loval>0x7fff)loval=0x7fff;
 415             if(hival>0x7fff)hival=0x7fff;
 416             c->dec_firsttable[bitreverse(word)]=
 417               0x80000000UL | (loval<<15) | hival;
 418           }
 419         }
 420       }
 421     }
 422   }
 423
 424   return(0);
 425  err_out:
 426   vorbis_book_clear(c);
 427   return(-1);
 428 }
 429
 430 long vorbis_book_codeword(codebook *book,int entry){
 431   if(book->c) /* only use with encode; decode optimizations are
 432                  allowed to break this */
 433     return book->codelist[entry];
 434   return -1;
 435 }
 436
 437 long vorbis_book_codelen(codebook *book,int entry){
 438   if(book->c) /* only use with encode; decode optimizations are
 439                  allowed to break this */
 440     return book->c->lengthlist[entry];
 441   return -1;
 442 }
 443
 444 #ifdef _V_SELFTEST
 445
 446 /* Unit tests of the dequantizer; this stuff will be OK
 447    cross-platform, I simply want to be sure that special mapping cases
 448    actually work properly; a bug could go unnoticed for a while */
 449
 450 #include <stdio.h>
 451
 452 /* cases:
 453
 454    no mapping
 455    full, explicit mapping
 456    algorithmic mapping
 457
 458    nonsequential
 459    sequential
 460 */
 461
 462 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
 463 static long partial_quantlist1[]={0,7,2};
 464
 465 /* no mapping */
 466 static_codebook test1={
 467   4,16,
 468   NULL,
 469   0,
 470   0,0,0,0,
 471   NULL,
 472   0
 473 };
 474 static float *test1_result=NULL;
 475
 476 /* linear, full mapping, nonsequential */
 477 static_codebook test2={
 478   4,3,
 479   NULL,
 480   2,
 481   -533200896,1611661312,4,0,
 482   full_quantlist1,
 483   0
 484 };
 485 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
 486
 487 /* linear, full mapping, sequential */
 488 static_codebook test3={
 489   4,3,
 490   NULL,
 491   2,
 492   -533200896,1611661312,4,1,
 493   full_quantlist1,
 494   0
 495 };
 496 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
 497
 498 /* linear, algorithmic mapping, nonsequential */
 499 static_codebook test4={
 500   3,27,
 501   NULL,
 502   1,
 503   -533200896,1611661312,4,0,
 504   partial_quantlist1,
 505   0
 506 };
 507 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
 508                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
 509                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
 510                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
 511                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
 512                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
 513                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
 514                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
 515                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
 516
 517 /* linear, algorithmic mapping, sequential */
 518 static_codebook test5={
 519   3,27,
 520   NULL,
 521   1,
 522   -533200896,1611661312,4,1,
 523   partial_quantlist1,
 524   0
 525 };
 526 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
 527                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
 528                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
 529                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
 530                               -3, 1, 5, 4, 8,12, -1, 3, 7,
 531                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
 532                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
 533                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
 534                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
 535
 536 void run_test(static_codebook *b,float *comp){
 537   float *out=_book_unquantize(b,b->entries,NULL);
 538   int i;
 539
 540   if(comp){
 541     if(!out){
 542       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
 543       exit(1);
 544     }
 545
 546     for(i=0;i<b->entries*b->dim;i++)
 547       if(fabs(out[i]-comp[i])>.0001){
 548         fprintf(stderr,"disagreement in unquantized and reference data:\n"
 549                 "position %d, %g != %g\n",i,out[i],comp[i]);
 550         exit(1);
 551       }
 552
 553   }else{
 554     if(out){
 555       fprintf(stderr,"_book_unquantize returned a value array: \n"
 556               " correct result should have been NULL\n");
 557       exit(1);
 558     }
 559   }
 560 }
 561
 562 int main(){
 563   /* run the nine dequant tests, and compare to the hand-rolled results */
 564   fprintf(stderr,"Dequant test 1... ");
 565   run_test(&test1,test1_result);
 566   fprintf(stderr,"OK\nDequant test 2... ");
 567   run_test(&test2,test2_result);
 568   fprintf(stderr,"OK\nDequant test 3... ");
 569   run_test(&test3,test3_result);
 570   fprintf(stderr,"OK\nDequant test 4... ");
 571   run_test(&test4,test4_result);
 572   fprintf(stderr,"OK\nDequant test 5... ");
 573   run_test(&test5,test5_result);
 574   fprintf(stderr,"OK\n\n");
 575
 576   return(0);
 577 }
 578
 579 #endif