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