vq/bookutil.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: utility functions for loading .vqh and .vqd files
  15  last mod: $Id: bookutil.c,v 1.11 2000/02/21 01:12:53 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdlib.h>
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <string.h>
  23 #include <errno.h>
  24 #include "vorbis/codebook.h"
  25 #include "bookutil.h"
  26
  27 void codebook_unquantize(codebook *b){
  28   long j,k;
  29   const static_codebook *c=b->c;
  30   double mindel=float24_unpack(c->q_min);
  31   double delta=float24_unpack(c->q_delta);
  32   if(!b->valuelist)b->valuelist=malloc(sizeof(double)*c->entries*c->dim);
  33
  34   for(j=0;j<c->entries;j++){
  35     double last=0.;
  36     for(k=0;k<c->dim;k++){
  37       double val=c->quantlist[j*c->dim+k]*delta+last+mindel;
  38       b->valuelist[j*c->dim+k]=val;
  39       if(c->q_sequencep)last=val;
  40
  41     }
  42   }
  43 }
  44
  45 /* A few little utils for reading files */
  46 /* read a line.  Use global, persistent buffering */
  47 static char *linebuffer=NULL;
  48 static int  lbufsize=0;
  49 char *get_line(FILE *in){
  50   long sofar=0;
  51   if(feof(in))return NULL;
  52
  53   while(1){
  54     int gotline=0;
  55
  56     while(!gotline){
  57       if(sofar+1>=lbufsize){
  58         if(!lbufsize){
  59           lbufsize=1024;
  60           linebuffer=malloc(lbufsize);
  61         }else{
  62           lbufsize*=2;
  63           linebuffer=realloc(linebuffer,lbufsize);
  64         }
  65       }
  66       {
  67         long c=fgetc(in);
  68         switch(c){
  69         case EOF:
  70           if(sofar==0)return(NULL);
  71           /* fallthrough correct */
  72         case '\n':
  73           linebuffer[sofar]='\0';
  74           gotline=1;
  75           break;
  76         default:
  77           linebuffer[sofar++]=c;
  78           linebuffer[sofar]='\0';
  79           break;
  80         }
  81       }
  82     }
  83
  84     if(linebuffer[0]=='#'){
  85       sofar=0;
  86     }else{
  87       return(linebuffer);
  88     }
  89   }
  90 }
  91
  92 /* read the next numerical value from the given file */
  93 static char *value_line_buff=NULL;
  94
  95 int get_line_value(FILE *in,double *value){
  96   char *next;
  97
  98   if(!value_line_buff)return(-1);
  99
 100   *value=strtod(value_line_buff, &next);
 101   if(next==value_line_buff){
 102     value_line_buff=NULL;
 103     return(-1);
 104   }else{
 105     value_line_buff=next;
 106     while(*value_line_buff>44)value_line_buff++;
 107     if(*value_line_buff==44)value_line_buff++;
 108     return(0);
 109   }
 110 }
 111
 112 int get_next_value(FILE *in,double *value){
 113   while(1){
 114     if(get_line_value(in,value)){
 115       value_line_buff=get_line(in);
 116       if(!value_line_buff)return(-1);
 117     }else{
 118       return(0);
 119     }
 120   }
 121 }
 122
 123 int get_next_ivalue(FILE *in,long *ivalue){
 124   double value;
 125   int ret=get_next_value(in,&value);
 126   *ivalue=value;
 127   return(ret);
 128 }
 129
 130 static double sequence_base=0.;
 131 static int v_sofar=0;
 132 void reset_next_value(void){
 133   value_line_buff=NULL;
 134   sequence_base=0.;
 135   v_sofar=0;
 136 }
 137
 138 int get_vector(codebook *b,FILE *in,int start, int n,double *a){
 139   int i;
 140   const static_codebook *c=b->c;
 141
 142   while(1){
 143
 144     if(v_sofar==n || get_line_value(in,a)){
 145       reset_next_value();
 146       if(get_next_value(in,a))
 147         break;
 148       for(i=0;i<start;i++){
 149         sequence_base=*a;
 150         get_line_value(in,a);
 151       }
 152     }
 153
 154     for(i=1;i<c->dim;i++)
 155       if(get_line_value(in,a+i))
 156         break;
 157
 158     if(i==c->dim){
 159       double temp=a[c->dim-1];
 160       for(i=0;i<c->dim;i++)a[i]-=sequence_base;
 161       if(c->q_sequencep)sequence_base=temp;
 162       v_sofar++;
 163       return(0);
 164     }
 165     sequence_base=0.;
 166   }
 167
 168   return(-1);
 169 }
 170
 171 /* read lines fromt he beginning until we find one containing the
 172    specified string */
 173 char *find_seek_to(FILE *in,char *s){
 174   rewind(in);
 175   while(1){
 176     char *line=get_line(in);
 177     if(line){
 178       if(strstr(line,s))
 179         return(line);
 180     }else
 181       return(NULL);
 182   }
 183 }
 184
 185
 186 /* this reads the format as written by vqbuild; innocent (legal)
 187    tweaking of the file that would not affect its valid header-ness
 188    will break this routine */
 189
 190 codebook *codebook_load(char *filename){
 191   codebook *b=calloc(1,sizeof(codebook));
 192   static_codebook *c=(static_codebook *)(b->c=calloc(1,sizeof(static_codebook)));
 193   encode_aux *a=calloc(1,sizeof(encode_aux));
 194   FILE *in=fopen(filename,"r");
 195   char *line;
 196   long i;
 197
 198   c->encode_tree=a;
 199
 200   if(in==NULL){
 201     fprintf(stderr,"Couldn't open codebook %s\n",filename);
 202     exit(1);
 203   }
 204
 205   /* find the codebook struct */
 206   find_seek_to(in,"static static_codebook _vq_book_");
 207
 208   /* get the major important values */
 209   line=get_line(in);
 210   if(sscanf(line,"%ld, %ld, %ld, %ld, %d, %d",
 211             &(c->dim),&(c->entries),&(c->q_min),&(c->q_delta),&(c->q_quant),
 212             &(c->q_sequencep))!=6){
 213     fprintf(stderr,"1: syntax in %s in line:\t %s",filename,line);
 214     exit(1);
 215   }
 216
 217   /* find the auxiliary encode struct (if any) */
 218   find_seek_to(in,"static encode_aux _vq_aux_");
 219   /* how big? */
 220   line=get_line(in);
 221   line=get_line(in);
 222   line=get_line(in);
 223   line=get_line(in);
 224   line=get_line(in);
 225   if(sscanf(line,"%ld, %ld",&(a->aux),&(a->alloc))!=2){
 226     fprintf(stderr,"2: syntax in %s in line:\t %s",filename,line);
 227     exit(1);
 228   }
 229
 230   /* load the quantized entries */
 231   find_seek_to(in,"static long _vq_quantlist_");
 232   reset_next_value();
 233   c->quantlist=malloc(sizeof(long)*c->entries*c->dim);
 234   for(i=0;i<c->entries*c->dim;i++)
 235     if(get_next_ivalue(in,c->quantlist+i)){
 236       fprintf(stderr,"out of data while reading codebook %s\n",filename);
 237       exit(1);
 238     }
 239
 240   /* load the lengthlist */
 241   find_seek_to(in,"static long _vq_lengthlist");
 242   reset_next_value();
 243   c->lengthlist=malloc(sizeof(long)*c->entries);
 244   for(i=0;i<c->entries;i++)
 245     if(get_next_ivalue(in,c->lengthlist+i)){
 246       fprintf(stderr,"out of data while reading codebook %s\n",filename);
 247       exit(1);
 248     }
 249
 250   /* load ptr0 */
 251   find_seek_to(in,"static long _vq_ptr0");
 252   reset_next_value();
 253   a->ptr0=malloc(sizeof(long)*a->aux);
 254   for(i=0;i<a->aux;i++)
 255     if(get_next_ivalue(in,a->ptr0+i)){
 256       fprintf(stderr,"out of data while reading codebook %s\n",filename);
 257       exit(1);
 258     }
 259
 260   /* load ptr1 */
 261   find_seek_to(in,"static long _vq_ptr1");
 262   reset_next_value();
 263   a->ptr1=malloc(sizeof(long)*a->aux);
 264   for(i=0;i<a->aux;i++)
 265     if(get_next_ivalue(in,a->ptr1+i)){
 266       fprintf(stderr,"out of data while reading codebook %s\n",filename);
 267       exit(1);
 268     }
 269
 270
 271   /* load p */
 272   find_seek_to(in,"static long _vq_p_");
 273   reset_next_value();
 274   a->p=malloc(sizeof(long)*a->aux);
 275   for(i=0;i<a->aux;i++)
 276     if(get_next_ivalue(in,a->p+i)){
 277       fprintf(stderr,"out of data while reading codebook %s\n",filename);
 278       exit(1);
 279     }
 280
 281   /* load q */
 282   find_seek_to(in,"static long _vq_q_");
 283   reset_next_value();
 284   a->q=malloc(sizeof(long)*a->aux);
 285   for(i=0;i<a->aux;i++)
 286     if(get_next_ivalue(in,a->q+i)){
 287       fprintf(stderr,"out of data while reading codebook %s\n",filename);
 288       exit(1);
 289     }
 290
 291   /* got it all */
 292   fclose(in);
 293
 294   /* might as well unquantize the entries while we're at it */
 295   codebook_unquantize(b);
 296
 297   /* don't need n and c */
 298   return(b);
 299 }
 300
 301
 302 /* the new version!  we have ptr0[0] distinct trees in the auxiliary
 303    encoding structure.  Find the best match in each one, choosing the
 304    best global match */
 305
 306 int codebook_entry(codebook *b,double *val){
 307   const static_codebook *c=b->c;
 308   encode_aux *t=c->encode_tree;
 309   int trees=t->ptr0[0];
 310
 311   /*{
 312     brute force
 313     double this,best=_dist(c->dim,val,b->valuelist);
 314     int i;
 315     for(i=1;i<c->entries;i++){
 316       this=_dist(c->dim,val,b->valuelist+i*c->dim);
 317       if(this<best){
 318         ptr=-i;
 319         best=this;
 320       }
 321     }
 322   }*/
 323
 324   double *n=alloca(c->dim*sizeof(double));
 325   double bestmetric;
 326   double best=-1;
 327   while(trees>0){
 328     int ptr=t->ptr0[--trees],k;
 329
 330     while(1){
 331       double C=0.;
 332       double *p=b->valuelist+t->p[ptr];
 333       double *q=b->valuelist+t->q[ptr];
 334
 335       for(k=0;k<c->dim;k++){
 336         n[k]=p[k]-q[k];
 337         C-=(p[k]+q[k])*n[k];
 338       }
 339       C/=2.;
 340
 341       for(k=0;k<c->dim;k++)
 342         C+=n[k]*val[k];
 343
 344       if(C>0.) /* in A */
 345         ptr= -t->ptr0[ptr];
 346       else     /* in B */
 347         ptr= -t->ptr1[ptr];
 348       if(ptr<=0)break;
 349     }
 350
 351     {
 352       double dist=0.;
 353       for(k=0;k<c->dim;k++){
 354         double one=val[k]-(b->valuelist-ptr*c->dim)[k];
 355         dist+=one*one;
 356       }
 357       if(best==-1 || dist<bestmetric){
 358         best=-ptr;
 359         bestmetric=dist;
 360       }
 361     }
 362   }
 363
 364   return(best);
 365 }
 366
 367 /* 24 bit float (not IEEE; nonnormalized mantissa +
 368    biased exponent ): neeeeemm mmmmmmmm mmmmmmmm */
 369
 370 long float24_pack(double val){
 371   int sign=0;
 372   long exp;
 373   long mant;
 374   if(val<0){
 375     sign=0x800000;
 376     val= -val;
 377   }
 378   exp= floor(log(val)/log(2));
 379   mant=rint(ldexp(val,17-exp));
 380   exp=(exp+VQ_FEXP_BIAS)<<18;
 381
 382   return(sign|exp|mant);
 383 }
 384
 385 double float24_unpack(long val){
 386   double mant=val&0x3ffff;
 387   double sign=val&0x800000;
 388   double exp =(val&0x7c0000)>>18;
 389   if(sign)mant= -mant;
 390   return(ldexp(mant,exp-17-VQ_FEXP_BIAS));
 391 }
 392
 393 void spinnit(char *s,int n){
 394   static int p=0;
 395   static long lasttime=0;
 396   long test;
 397   struct timeval thistime;
 398
 399   gettimeofday(&thistime,NULL);
 400   test=thistime.tv_sec*10+thistime.tv_usec/100000;
 401   if(lasttime!=test){
 402     lasttime=test;
 403
 404     fprintf(stderr,"%s%d ",s,n);
 405
 406     p++;if(p>3)p=0;
 407     switch(p){
 408     case 0:
 409       fprintf(stderr,"|    \r");
 410       break;
 411     case 1:
 412       fprintf(stderr,"/    \r");
 413       break;
 414     case 2:
 415       fprintf(stderr,"-    \r");
 416       break;
 417     case 3:
 418       fprintf(stderr,"\\    \r");
 419       break;
 420     }
 421     fflush(stderr);
 422   }
 423 }
 424
 425 void build_tree_from_lengths(int vals, long *hist, long *lengths){
 426   int i,j;
 427   long *membership=malloc(vals*sizeof(long));
 428
 429   for(i=0;i<vals;i++)membership[i]=i;
 430
 431   /* find codeword lengths */
 432   /* much more elegant means exist.  Brute force n^2, minimum thought */
 433   for(i=vals;i>1;i--){
 434     int first=-1,second=-1;
 435     long least=-1;
 436
 437     spinnit("building... ",i);
 438
 439     /* find the two nodes to join */
 440     for(j=0;j<vals;j++)
 441       if(least==-1 || hist[j]<least){
 442         least=hist[j];
 443         first=membership[j];
 444       }
 445     least=-1;
 446     for(j=0;j<vals;j++)
 447       if((least==-1 || hist[j]<least) && membership[j]!=first){
 448         least=hist[j];
 449         second=membership[j];
 450       }
 451     if(first==-1 || second==-1){
 452       fprintf(stderr,"huffman fault; no free branch\n");
 453       exit(1);
 454     }
 455
 456     /* join them */
 457     least=hist[first]+hist[second];
 458     for(j=0;j<vals;j++)
 459       if(membership[j]==first || membership[j]==second){
 460         membership[j]=first;
 461         hist[j]=least;
 462         lengths[j]++;
 463       }
 464   }
 465   for(i=0;i<vals-1;i++)
 466     if(membership[i]!=membership[i+1]){
 467       fprintf(stderr,"huffman fault; failed to build single tree\n");
 468       exit(1);
 469     }
 470
 471   free(membership);
 472 }