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_clear(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   }
 267 }
 268
 269 void vorbis_staticbook_destroy(static_codebook *b){
 270   if(b->allocedp){
 271     vorbis_staticbook_clear(b);
 272     _ogg_free(b);
 273   }
 274 }
 275
 276 void vorbis_book_clear(codebook *b){
 277   /* static book is not cleared; we're likely called on the lookup and
 278      the static codebook belongs to the info struct */
 279   if(b->valuelist)_ogg_free(b->valuelist);
 280   if(b->codelist)_ogg_free(b->codelist);
 281
 282   if(b->dec_index)_ogg_free(b->dec_index);
 283   if(b->dec_codelengths)_ogg_free(b->dec_codelengths);
 284   if(b->dec_firsttable)_ogg_free(b->dec_firsttable);
 285
 286   memset(b,0,sizeof(*b));
 287 }
 288
 289 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
 290
 291   memset(c,0,sizeof(*c));
 292   c->c=s;
 293   c->entries=s->entries;
 294   c->used_entries=s->entries;
 295   c->dim=s->dim;
 296   c->codelist=_make_words(s->lengthlist,s->entries,0);
 297   //c->valuelist=_book_unquantize(s,s->entries,NULL);
 298   c->quantvals=_book_maptype1_quantvals(s);
 299   c->minval=(int)rint(_float32_unpack(s->q_min));
 300   c->delta=(int)rint(_float32_unpack(s->q_delta));
 301
 302   return(0);
 303 }
 304
 305 static ogg_uint32_t bitreverse(ogg_uint32_t x){
 306   x=    ((x>>16)&0x0000ffffUL) | ((x<<16)&0xffff0000UL);
 307   x=    ((x>> 8)&0x00ff00ffUL) | ((x<< 8)&0xff00ff00UL);
 308   x=    ((x>> 4)&0x0f0f0f0fUL) | ((x<< 4)&0xf0f0f0f0UL);
 309   x=    ((x>> 2)&0x33333333UL) | ((x<< 2)&0xccccccccUL);
 310   return((x>> 1)&0x55555555UL) | ((x<< 1)&0xaaaaaaaaUL);
 311 }
 312
 313 static int sort32a(const void *a,const void *b){
 314   return ( **(ogg_uint32_t **)a>**(ogg_uint32_t **)b)-
 315     ( **(ogg_uint32_t **)a<**(ogg_uint32_t **)b);
 316 }
 317
 318 /* decode codebook arrangement is more heavily optimized than encode */
 319 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
 320   int i,j,n=0,tabn;
 321   int *sortindex;
 322   memset(c,0,sizeof(*c));
 323
 324   /* count actually used entries */
 325   for(i=0;i<s->entries;i++)
 326     if(s->lengthlist[i]>0)
 327       n++;
 328
 329   c->entries=s->entries;
 330   c->used_entries=n;
 331   c->dim=s->dim;
 332
 333   if(n>0){
 334
 335     /* two different remappings go on here.
 336
 337     First, we collapse the likely sparse codebook down only to
 338     actually represented values/words.  This collapsing needs to be
 339     indexed as map-valueless books are used to encode original entry
 340     positions as integers.
 341
 342     Second, we reorder all vectors, including the entry index above,
 343     by sorted bitreversed codeword to allow treeless decode. */
 344
 345     /* perform sort */
 346     ogg_uint32_t *codes=_make_words(s->lengthlist,s->entries,c->used_entries);
 347     ogg_uint32_t **codep=alloca(sizeof(*codep)*n);
 348
 349     if(codes==NULL)goto err_out;
 350
 351     for(i=0;i<n;i++){
 352       codes[i]=bitreverse(codes[i]);
 353       codep[i]=codes+i;
 354     }
 355
 356     qsort(codep,n,sizeof(*codep),sort32a);
 357
 358     sortindex=alloca(n*sizeof(*sortindex));
 359     c->codelist=_ogg_malloc(n*sizeof(*c->codelist));
 360     /* the index is a reverse index */
 361     for(i=0;i<n;i++){
 362       int position=codep[i]-codes;
 363       sortindex[position]=i;
 364     }
 365
 366     for(i=0;i<n;i++)
 367       c->codelist[sortindex[i]]=codes[i];
 368     _ogg_free(codes);
 369
 370
 371     c->valuelist=_book_unquantize(s,n,sortindex);
 372     c->dec_index=_ogg_malloc(n*sizeof(*c->dec_index));
 373
 374     for(n=0,i=0;i<s->entries;i++)
 375       if(s->lengthlist[i]>0)
 376         c->dec_index[sortindex[n++]]=i;
 377
 378     c->dec_codelengths=_ogg_malloc(n*sizeof(*c->dec_codelengths));
 379     for(n=0,i=0;i<s->entries;i++)
 380       if(s->lengthlist[i]>0)
 381         c->dec_codelengths[sortindex[n++]]=s->lengthlist[i];
 382
 383     c->dec_firsttablen=_ilog(c->used_entries)-4; /* this is magic */
 384     if(c->dec_firsttablen<5)c->dec_firsttablen=5;
 385     if(c->dec_firsttablen>8)c->dec_firsttablen=8;
 386
 387     tabn=1<<c->dec_firsttablen;
 388     c->dec_firsttable=_ogg_calloc(tabn,sizeof(*c->dec_firsttable));
 389     c->dec_maxlength=0;
 390
 391     for(i=0;i<n;i++){
 392       if(c->dec_maxlength<c->dec_codelengths[i])
 393         c->dec_maxlength=c->dec_codelengths[i];
 394       if(c->dec_codelengths[i]<=c->dec_firsttablen){
 395         ogg_uint32_t orig=bitreverse(c->codelist[i]);
 396         for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
 397           c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
 398       }
 399     }
 400
 401     /* now fill in 'unused' entries in the firsttable with hi/lo search
 402        hints for the non-direct-hits */
 403     {
 404       ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
 405       long lo=0,hi=0;
 406
 407       for(i=0;i<tabn;i++){
 408         ogg_uint32_t word=i<<(32-c->dec_firsttablen);
 409         if(c->dec_firsttable[bitreverse(word)]==0){
 410           while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
 411           while(    hi<n && word>=(c->codelist[hi]&mask))hi++;
 412
 413           /* we only actually have 15 bits per hint to play with here.
 414              In order to overflow gracefully (nothing breaks, efficiency
 415              just drops), encode as the difference from the extremes. */
 416           {
 417             unsigned long loval=lo;
 418             unsigned long hival=n-hi;
 419
 420             if(loval>0x7fff)loval=0x7fff;
 421             if(hival>0x7fff)hival=0x7fff;
 422             c->dec_firsttable[bitreverse(word)]=
 423               0x80000000UL | (loval<<15) | hival;
 424           }
 425         }
 426       }
 427     }
 428   }
 429
 430   return(0);
 431  err_out:
 432   vorbis_book_clear(c);
 433   return(-1);
 434 }
 435
 436 long vorbis_book_codeword(codebook *book,int entry){
 437   if(book->c) /* only use with encode; decode optimizations are
 438                  allowed to break this */
 439     return book->codelist[entry];
 440   return -1;
 441 }
 442
 443 long vorbis_book_codelen(codebook *book,int entry){
 444   if(book->c) /* only use with encode; decode optimizations are
 445                  allowed to break this */
 446     return book->c->lengthlist[entry];
 447   return -1;
 448 }
 449
 450 #ifdef _V_SELFTEST
 451
 452 /* Unit tests of the dequantizer; this stuff will be OK
 453    cross-platform, I simply want to be sure that special mapping cases
 454    actually work properly; a bug could go unnoticed for a while */
 455
 456 #include <stdio.h>
 457
 458 /* cases:
 459
 460    no mapping
 461    full, explicit mapping
 462    algorithmic mapping
 463
 464    nonsequential
 465    sequential
 466 */
 467
 468 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
 469 static long partial_quantlist1[]={0,7,2};
 470
 471 /* no mapping */
 472 static_codebook test1={
 473   4,16,
 474   NULL,
 475   0,
 476   0,0,0,0,
 477   NULL,
 478   0
 479 };
 480 static float *test1_result=NULL;
 481
 482 /* linear, full mapping, nonsequential */
 483 static_codebook test2={
 484   4,3,
 485   NULL,
 486   2,
 487   -533200896,1611661312,4,0,
 488   full_quantlist1,
 489   0
 490 };
 491 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
 492
 493 /* linear, full mapping, sequential */
 494 static_codebook test3={
 495   4,3,
 496   NULL,
 497   2,
 498   -533200896,1611661312,4,1,
 499   full_quantlist1,
 500   0
 501 };
 502 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
 503
 504 /* linear, algorithmic mapping, nonsequential */
 505 static_codebook test4={
 506   3,27,
 507   NULL,
 508   1,
 509   -533200896,1611661312,4,0,
 510   partial_quantlist1,
 511   0
 512 };
 513 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
 514                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
 515                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
 516                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
 517                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
 518                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
 519                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
 520                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
 521                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
 522
 523 /* linear, algorithmic mapping, sequential */
 524 static_codebook test5={
 525   3,27,
 526   NULL,
 527   1,
 528   -533200896,1611661312,4,1,
 529   partial_quantlist1,
 530   0
 531 };
 532 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
 533                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
 534                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
 535                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
 536                               -3, 1, 5, 4, 8,12, -1, 3, 7,
 537                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
 538                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
 539                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
 540                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
 541
 542 void run_test(static_codebook *b,float *comp){
 543   float *out=_book_unquantize(b,b->entries,NULL);
 544   int i;
 545
 546   if(comp){
 547     if(!out){
 548       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
 549       exit(1);
 550     }
 551
 552     for(i=0;i<b->entries*b->dim;i++)
 553       if(fabs(out[i]-comp[i])>.0001){
 554         fprintf(stderr,"disagreement in unquantized and reference data:\n"
 555                 "position %d, %g != %g\n",i,out[i],comp[i]);
 556         exit(1);
 557       }
 558
 559   }else{
 560     if(out){
 561       fprintf(stderr,"_book_unquantize returned a value array: \n"
 562               " correct result should have been NULL\n");
 563       exit(1);
 564     }
 565   }
 566 }
 567
 568 int main(){
 569   /* run the nine dequant tests, and compare to the hand-rolled results */
 570   fprintf(stderr,"Dequant test 1... ");
 571   run_test(&test1,test1_result);
 572   fprintf(stderr,"OK\nDequant test 2... ");
 573   run_test(&test2,test2_result);
 574   fprintf(stderr,"OK\nDequant test 3... ");
 575   run_test(&test3,test3_result);
 576   fprintf(stderr,"OK\nDequant test 4... ");
 577   run_test(&test4,test4_result);
 578   fprintf(stderr,"OK\nDequant test 5... ");
 579   run_test(&test5,test5_result);
 580   fprintf(stderr,"OK\n\n");
 581
 582   return(0);
 583 }
 584
 585 #endif