libvorbis-1.0.1/vq/latticehint.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-2001             *
   9  * by the XIPHOPHORUS Company http://www.xiph.org/                  *
  10  *                                                                  *
  11  ********************************************************************
  12
  13  function: utility main for building thresh/pigeonhole encode hints
  14  last mod: $Id: latticehint.c,v 1.12 2001/12/20 01:00:39 segher Exp $
  15
  16  ********************************************************************/
  17
  18 #include <stdlib.h>
  19 #include <stdio.h>
  20 #include <math.h>
  21 #include <string.h>
  22 #include <errno.h>
  23 #include "../lib/scales.h"
  24 #include "bookutil.h"
  25 #include "vqgen.h"
  26 #include "vqsplit.h"
  27
  28 /* The purpose of this util is to build encode hints for lattice
  29    codebooks so that brute forcing each codebook entry isn't needed.
  30    Threshhold hints are for books in which each scalar in the vector
  31    is independant (eg, residue) and pigeonhole lookups provide a
  32    minimum error fit for words where the scalars are interdependant
  33    (each affecting the fit of the next in sequence) as in an LSP
  34    sequential book (or can be used along with a sparse threshhold map,
  35    like a splitting tree that need not be trained)
  36
  37    If the input book is non-sequential, a threshhold hint is built.
  38    If the input book is sequential, a pigeonholing hist is built.
  39    If the book is sparse, a pigeonholing hint is built, possibly in addition
  40      to the threshhold hint
  41
  42    command line:
  43    latticehint book.vqh [threshlist]
  44
  45    latticehint produces book.vqh on stdout */
  46
  47 static int longsort(const void *a, const void *b){
  48   return(**((long **)a)-**((long **)b));
  49 }
  50
  51 static int addtosearch(int entry,long **tempstack,long *tempcount,int add){
  52   long *ptr=tempstack[entry];
  53   long i=tempcount[entry];
  54
  55   if(ptr){
  56     while(i--)
  57       if(*ptr++==add)return(0);
  58     tempstack[entry]=_ogg_realloc(tempstack[entry],
  59                              (tempcount[entry]+1)*sizeof(long));
  60   }else{
  61     tempstack[entry]=_ogg_malloc(sizeof(long));
  62   }
  63
  64   tempstack[entry][tempcount[entry]++]=add;
  65   return(1);
  66 }
  67
  68 static void setvals(int dim,encode_aux_pigeonhole *p,
  69                     long *temptrack,float *tempmin,float *tempmax,
  70                     int seqp){
  71   int i;
  72   float last=0.f;
  73   for(i=0;i<dim;i++){
  74     tempmin[i]=(temptrack[i])*p->del+p->min+last;
  75     tempmax[i]=tempmin[i]+p->del;
  76     if(seqp)last=tempmin[i];
  77   }
  78 }
  79
  80 /* note that things are currently set up such that input fits that
  81    quantize outside the pigeonmap are dropped and brute-forced.  So we
  82    can ignore the <0 and >=n boundary cases in min/max error */
  83
  84 static float minerror(int dim,float *a,encode_aux_pigeonhole *p,
  85                        long *temptrack,float *tempmin,float *tempmax){
  86   int i;
  87   float err=0.f;
  88   for(i=0;i<dim;i++){
  89     float eval=0.f;
  90     if(a[i]<tempmin[i]){
  91       eval=tempmin[i]-a[i];
  92     }else if(a[i]>tempmax[i]){
  93       eval=a[i]-tempmax[i];
  94     }
  95     err+=eval*eval;
  96   }
  97   return(err);
  98 }
  99
 100 static float maxerror(int dim,float *a,encode_aux_pigeonhole *p,
 101                        long *temptrack,float *tempmin,float *tempmax){
 102   int i;
 103   float err=0.f,eval;
 104   for(i=0;i<dim;i++){
 105     if(a[i]<tempmin[i]){
 106       eval=tempmax[i]-a[i];
 107     }else if(a[i]>tempmax[i]){
 108       eval=a[i]-tempmin[i];
 109     }else{
 110       float t1=a[i]-tempmin[i];
 111       eval=tempmax[i]-a[i];
 112       if(t1>eval)eval=t1;
 113     }
 114     err+=eval*eval;
 115   }
 116   return(err);
 117 }
 118
 119 int main(int argc,char *argv[]){
 120   codebook *b;
 121   static_codebook *c;
 122   int entries=-1,dim=-1;
 123   float min,del;
 124   char *name;
 125   long i,j;
 126   float *suggestions;
 127   int suggcount=0;
 128
 129   if(argv[1]==NULL){
 130     fprintf(stderr,"Need a lattice book on the command line.\n");
 131     exit(1);
 132   }
 133
 134   {
 135     char *ptr;
 136     char *filename=strdup(argv[1]);
 137
 138     b=codebook_load(filename);
 139     c=(static_codebook *)(b->c);
 140
 141     ptr=strrchr(filename,'.');
 142     if(ptr){
 143       *ptr='\0';
 144       name=strdup(filename);
 145     }else{
 146       name=strdup(filename);
 147     }
 148   }
 149
 150   if(c->maptype!=1){
 151     fprintf(stderr,"Provided book is not a latticebook.\n");
 152     exit(1);
 153   }
 154
 155   entries=b->entries;
 156   dim=b->dim;
 157   min=_float32_unpack(c->q_min);
 158   del=_float32_unpack(c->q_delta);
 159
 160   /* Do we want to gen a threshold hint? */
 161   if(c->q_sequencep==0){
 162     /* yes. Discard any preexisting threshhold hint */
 163     long quantvals=_book_maptype1_quantvals(c);
 164     long **quantsort=alloca(quantvals*sizeof(long *));
 165     encode_aux_threshmatch *t=_ogg_calloc(1,sizeof(encode_aux_threshmatch));
 166     c->thresh_tree=t;
 167
 168     fprintf(stderr,"Adding threshold hint to %s...\n",name);
 169
 170     /* partial/complete suggestions */
 171     if(argv[2]){
 172       char *ptr=strdup(argv[2]);
 173       suggestions=alloca(sizeof(float)*quantvals);
 174
 175       for(suggcount=0;ptr && suggcount<quantvals;suggcount++){
 176         char *ptr2=strchr(ptr,',');
 177         if(ptr2)*ptr2++='\0';
 178         suggestions[suggcount]=atof(ptr);
 179         ptr=ptr2;
 180       }
 181     }
 182
 183     /* simplest possible threshold hint only */
 184     t->quantthresh=_ogg_calloc(quantvals-1,sizeof(float));
 185     t->quantmap=_ogg_calloc(quantvals,sizeof(int));
 186     t->threshvals=quantvals;
 187     t->quantvals=quantvals;
 188
 189     /* the quantvals may not be in order; sort em first */
 190     for(i=0;i<quantvals;i++)quantsort[i]=c->quantlist+i;
 191     qsort(quantsort,quantvals,sizeof(long *),longsort);
 192
 193     /* ok, gen the map and thresholds */
 194     for(i=0;i<quantvals;i++)t->quantmap[i]=quantsort[i]-c->quantlist;
 195     for(i=0;i<quantvals-1;i++){
 196       float v1=*(quantsort[i])*del+min;
 197       float v2=*(quantsort[i+1])*del+min;
 198
 199       for(j=0;j<suggcount;j++)
 200         if(v1<suggestions[j] && suggestions[j]<v2){
 201           t->quantthresh[i]=suggestions[j];
 202           break;
 203         }
 204
 205       if(j==suggcount){
 206         t->quantthresh[i]=(v1+v2)*.5;
 207       }
 208     }
 209   }
 210
 211   /* Do we want to gen a pigeonhole hint? */
 212 #if 0
 213   for(i=0;i<entries;i++)if(c->lengthlist[i]==0)break;
 214   if(c->q_sequencep || i<entries){
 215     long **tempstack;
 216     long *tempcount;
 217     long *temptrack;
 218     float *tempmin;
 219     float *tempmax;
 220     long totalstack=0;
 221     long pigeons;
 222     long subpigeons;
 223     long quantvals=_book_maptype1_quantvals(c);
 224     int changep=1,factor;
 225
 226     encode_aux_pigeonhole *p=_ogg_calloc(1,sizeof(encode_aux_pigeonhole));
 227     c->pigeon_tree=p;
 228
 229     fprintf(stderr,"Adding pigeonhole hint to %s...\n",name);
 230
 231     /* the idea is that we quantize uniformly, even in a nonuniform
 232        lattice, so that quantization of one scalar has a predictable
 233        result on the next sequential scalar in a greedy matching
 234        algorithm.  We generate a lookup based on the quantization of
 235        the vector (pigeonmap groups quantized entries together) and
 236        list the entries that could possible be the best fit for any
 237        given member of that pigeonhole.  The encode process then has a
 238        much smaller list to brute force */
 239
 240     /* find our pigeonhole-specific quantization values, fill in the
 241        quant value->pigeonhole map */
 242     factor=3;
 243     p->del=del;
 244     p->min=min;
 245     p->quantvals=quantvals;
 246     {
 247       int max=0;
 248       for(i=0;i<quantvals;i++)if(max<c->quantlist[i])max=c->quantlist[i];
 249       p->mapentries=max;
 250     }
 251     p->pigeonmap=_ogg_malloc(p->mapentries*sizeof(long));
 252     p->quantvals=(quantvals+factor-1)/factor;
 253
 254     /* pigeonhole roughly on the boundaries of the quantvals; the
 255        exact pigeonhole grouping is an optimization issue, not a
 256        correctness issue */
 257     for(i=0;i<p->mapentries;i++){
 258       float thisval=del*i+min; /* middle of the quant zone */
 259       int quant=0;
 260       float err=fabs(c->quantlist[0]*del+min-thisval);
 261       for(j=1;j<quantvals;j++){
 262         float thiserr=fabs(c->quantlist[j]*del+min-thisval);
 263         if(thiserr<err){
 264           quant=j/factor;
 265           err=thiserr;
 266         }
 267       }
 268       p->pigeonmap[i]=quant;
 269     }
 270
 271     /* pigeonmap complete.  Now do the grungy business of finding the
 272     entries that could possibly be the best fit for a value appearing
 273     in the pigeonhole. The trick that allows the below to work is the
 274     uniform quantization; even though the scalars may be 'sequential'
 275     (each a delta from the last), the uniform quantization means that
 276     the error variance is *not* dependant.  Given a pigeonhole and an
 277     entry, we can find the minimum and maximum possible errors
 278     (relative to the entry) for any point that could appear in the
 279     pigeonhole */
 280
 281     /* must iterate over both pigeonholes and entries */
 282     /* temporarily (in order to avoid thinking hard), we grow each
 283        pigeonhole seperately, the build a stack of 'em later */
 284     pigeons=1;
 285     subpigeons=1;
 286     for(i=0;i<dim;i++)subpigeons*=p->mapentries;
 287     for(i=0;i<dim;i++)pigeons*=p->quantvals;
 288     temptrack=_ogg_calloc(dim,sizeof(long));
 289     tempmin=_ogg_calloc(dim,sizeof(float));
 290     tempmax=_ogg_calloc(dim,sizeof(float));
 291     tempstack=_ogg_calloc(pigeons,sizeof(long *));
 292     tempcount=_ogg_calloc(pigeons,sizeof(long));
 293
 294     while(1){
 295       float errorpost=-1;
 296       char buffer[80];
 297
 298       /* map our current pigeonhole to a 'big pigeonhole' so we know
 299          what list we're after */
 300       int entry=0;
 301       for(i=dim-1;i>=0;i--)entry=entry*p->quantvals+p->pigeonmap[temptrack[i]];
 302       setvals(dim,p,temptrack,tempmin,tempmax,c->q_sequencep);
 303       sprintf(buffer,"Building pigeonhole search list [%ld]...",totalstack);
 304
 305
 306       /* Search all entries to find the one with the minimum possible
 307          maximum error.  Record that error */
 308       for(i=0;i<entries;i++){
 309         if(c->lengthlist[i]>0){
 310           float this=maxerror(dim,b->valuelist+i*dim,p,
 311                                temptrack,tempmin,tempmax);
 312           if(errorpost==-1 || this<errorpost)errorpost=this;
 313           spinnit(buffer,subpigeons);
 314         }
 315       }
 316
 317       /* Our search list will contain all entries with a minimum
 318          possible error <= our errorpost */
 319       for(i=0;i<entries;i++)
 320         if(c->lengthlist[i]>0){
 321           spinnit(buffer,subpigeons);
 322           if(minerror(dim,b->valuelist+i*dim,p,
 323                       temptrack,tempmin,tempmax)<errorpost)
 324             totalstack+=addtosearch(entry,tempstack,tempcount,i);
 325         }
 326
 327       for(i=0;i<dim;i++){
 328         temptrack[i]++;
 329         if(temptrack[i]<p->mapentries)break;
 330         temptrack[i]=0;
 331       }
 332       if(i==dim)break;
 333       subpigeons--;
 334     }
 335
 336     fprintf(stderr,"\r                                                     "
 337             "\rTotal search list size (all entries): %ld\n",totalstack);
 338
 339     /* pare the index of lists for improbable quantizations (where
 340        improbable is determined by c->lengthlist; we assume that
 341        pigeonholing is in sync with the codeword cells, which it is */
 342     /*for(i=0;i<entries;i++){
 343       float probability= 1.f/(1<<c->lengthlist[i]);
 344       if(c->lengthlist[i]==0 || probability*entries<cutoff){
 345         totalstack-=tempcount[i];
 346         tempcount[i]=0;
 347       }
 348       }*/
 349
 350     /* pare the list of shortlists; merge contained and similar lists
 351        together */
 352     p->fitmap=_ogg_malloc(pigeons*sizeof(long));
 353     for(i=0;i<pigeons;i++)p->fitmap[i]=-1;
 354     while(changep){
 355       char buffer[80];
 356       changep=0;
 357
 358       for(i=0;i<pigeons;i++){
 359         if(p->fitmap[i]<0 && tempcount[i]){
 360           for(j=i+1;j<pigeons;j++){
 361             if(p->fitmap[j]<0 && tempcount[j]){
 362               /* is one list a superset, or are they sufficiently similar? */
 363               int amiss=0,bmiss=0,ii,jj;
 364               for(ii=0;ii<tempcount[i];ii++){
 365                 for(jj=0;jj<tempcount[j];jj++)
 366                   if(tempstack[i][ii]==tempstack[j][jj])break;
 367                 if(jj==tempcount[j])amiss++;
 368               }
 369               for(jj=0;jj<tempcount[j];jj++){
 370                 for(ii=0;ii<tempcount[i];ii++)
 371                   if(tempstack[i][ii]==tempstack[j][jj])break;
 372                 if(ii==tempcount[i])bmiss++;
 373               }
 374               if(amiss==0 ||
 375                  bmiss==0 ||
 376                  (amiss*2<tempcount[i] && bmiss*2<tempcount[j] &&
 377                  tempcount[i]+bmiss<entries/30)){
 378
 379                 /*superset/similar  Add all of one to the other. */
 380                 for(jj=0;jj<tempcount[j];jj++)
 381                   totalstack+=addtosearch(i,tempstack,tempcount,
 382                                           tempstack[j][jj]);
 383                 totalstack-=tempcount[j];
 384                 p->fitmap[j]=i;
 385                 changep=1;
 386               }
 387             }
 388           }
 389           sprintf(buffer,"Consolidating [%ld total, %s]... ",totalstack,
 390                   changep?"reit":"nochange");
 391           spinnit(buffer,pigeons-i);
 392         }
 393       }
 394     }
 395
 396     /* repack the temp stack in final form */
 397     fprintf(stderr,"\r                                                       "
 398             "\rFinal total list size: %ld\n",totalstack);
 399
 400
 401     p->fittotal=totalstack;
 402     p->fitlist=_ogg_malloc((totalstack+1)*sizeof(long));
 403     p->fitlength=_ogg_malloc(pigeons*sizeof(long));
 404     {
 405       long usage=0;
 406       for(i=0;i<pigeons;i++){
 407         if(p->fitmap[i]==-1){
 408           if(tempcount[i])
 409             memcpy(p->fitlist+usage,tempstack[i],tempcount[i]*sizeof(long));
 410           p->fitmap[i]=usage;
 411           p->fitlength[i]=tempcount[i];
 412           usage+=tempcount[i];
 413           if(usage>totalstack){
 414             fprintf(stderr,"Internal error; usage>totalstack\n");
 415             exit(1);
 416           }
 417         }else{
 418           p->fitlength[i]=p->fitlength[p->fitmap[i]];
 419           p->fitmap[i]=p->fitmap[p->fitmap[i]];
 420         }
 421       }
 422     }
 423   }
 424 #endif
 425
 426   write_codebook(stdout,name,c);
 427   fprintf(stderr,"\r                                                     "
 428           "\nDone.\n");
 429   exit(0);
 430 }