Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: downloads/libvorbis-1.2.0/vq/vqsplit.c @ 16

Last change on this file since 16 was 16, checked in by landauf, 17 years ago

added libvorbis

File size: 17.6 KB
Line 
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 Xiph.Org Foundation http://www.xiph.org/                  *
10 *                                                                  *
11 ********************************************************************
12
13 function: build a VQ codebook and the encoding decision 'tree'
14 last mod: $Id: vqsplit.c 13293 2007-07-24 00:09:47Z xiphmont $
15
16 ********************************************************************/
17
18/* This code is *not* part of libvorbis.  It is used to generate
19   trained codebooks offline and then spit the results into a
20   pregenerated codebook that is compiled into libvorbis.  It is an
21   expensive (but good) algorithm.  Run it on big iron. */
22
23/* There are so many optimizations to explore in *both* stages that
24   considering the undertaking is almost withering.  For now, we brute
25   force it all */
26
27#include <stdlib.h>
28#include <stdio.h>
29#include <math.h>
30#include <string.h>
31#include <sys/time.h>
32
33#include "vqgen.h"
34#include "vqsplit.h"
35#include "bookutil.h"
36
37/* Codebook generation happens in two steps:
38
39   1) Train the codebook with data collected from the encoder: We use
40   one of a few error metrics (which represent the distance between a
41   given data point and a candidate point in the training set) to
42   divide the training set up into cells representing roughly equal
43   probability of occurring.
44
45   2) Generate the codebook and auxiliary data from the trained data set
46*/
47
48/* Building a codebook from trained set **********************************
49
50   The codebook in raw form is technically finished once it's trained.
51   However, we want to finalize the representative codebook values for
52   each entry and generate auxiliary information to optimize encoding.
53   We generate the auxiliary coding tree using collected data,
54   probably the same data as in the original training */
55
56/* At each recursion, the data set is split in half.  Cells with data
57   points on side A go into set A, same with set B.  The sets may
58   overlap.  If the cell overlaps the deviding line only very slightly
59   (provided parameter), we may choose to ignore the overlap in order
60   to pare the tree down */
61
62long *isortvals;
63int iascsort(const void *a,const void *b){
64  long av=isortvals[*((long *)a)];
65  long bv=isortvals[*((long *)b)];
66  return(av-bv);
67}
68
69static float _Ndist(int el,float *a, float *b){
70  int i;
71  float acc=0.f;
72  for(i=0;i<el;i++){
73    float val=(a[i]-b[i]);
74    acc+=val*val;
75  }
76  return sqrt(acc);
77}
78
79#define _Npoint(i) (pointlist+dim*(i))
80#define _Nnow(i) (entrylist+dim*(i))
81
82
83/* goes through the split, but just counts it and returns a metric*/
84int vqsp_count(float *entrylist,float *pointlist,int dim,
85               long *membership,long *reventry,
86               long *entryindex,long entries, 
87               long *pointindex,long points,int splitp,
88               long *entryA,long *entryB,
89               long besti,long bestj,
90               long *entriesA,long *entriesB,long *entriesC){
91  long i,j;
92  long A=0,B=0,C=0;
93  long pointsA=0;
94  long pointsB=0;
95  long *temppointsA=NULL;
96  long *temppointsB=NULL;
97 
98  if(splitp){
99    temppointsA=_ogg_malloc(points*sizeof(long));
100    temppointsB=_ogg_malloc(points*sizeof(long));
101  }
102
103  memset(entryA,0,sizeof(long)*entries);
104  memset(entryB,0,sizeof(long)*entries);
105
106  /* Do the points belonging to this cell occur on sideA, sideB or
107     both? */
108
109  for(i=0;i<points;i++){
110    float *ppt=_Npoint(pointindex[i]);
111    long   firstentry=membership[pointindex[i]];
112
113    if(firstentry==besti){
114      entryA[reventry[firstentry]]=1;
115      if(splitp)temppointsA[pointsA++]=pointindex[i];
116      continue;
117    }
118    if(firstentry==bestj){
119      entryB[reventry[firstentry]]=1;
120      if(splitp)temppointsB[pointsB++]=pointindex[i];
121      continue;
122    }
123    {
124      float distA=_Ndist(dim,ppt,_Nnow(besti));
125      float distB=_Ndist(dim,ppt,_Nnow(bestj));
126      if(distA<distB){
127        entryA[reventry[firstentry]]=1;
128        if(splitp)temppointsA[pointsA++]=pointindex[i];
129      }else{
130        entryB[reventry[firstentry]]=1;
131        if(splitp)temppointsB[pointsB++]=pointindex[i];
132      }
133    }
134  }
135
136  /* The entry splitting isn't total, so that storage has to be
137     allocated for recursion.  Reuse the entryA/entryB vectors */
138  /* keep the entries in ascending order (relative to the original
139     list); we rely on that stability when ordering p/q choice */
140  for(j=0;j<entries;j++){
141    if(entryA[j] && entryB[j])C++;
142    if(entryA[j])entryA[A++]=entryindex[j];
143    if(entryB[j])entryB[B++]=entryindex[j];
144  }
145  *entriesA=A;
146  *entriesB=B;
147  *entriesC=C;
148  if(splitp){
149    memcpy(pointindex,temppointsA,sizeof(long)*pointsA);
150    memcpy(pointindex+pointsA,temppointsB,sizeof(long)*pointsB);
151    free(temppointsA);
152    free(temppointsB);
153  }
154  return(pointsA);
155}
156
157int lp_split(float *pointlist,long totalpoints,
158             codebook *b,
159             long *entryindex,long entries, 
160             long *pointindex,long points,
161             long *membership,long *reventry,
162             long depth, long *pointsofar){
163
164  encode_aux_nearestmatch *t=b->c->nearest_tree;
165
166  /* The encoder, regardless of book, will be using a straight
167     euclidian distance-to-point metric to determine closest point.
168     Thus we split the cells using the same (we've already trained the
169     codebook set spacing and distribution using special metrics and
170     even a midpoint division won't disturb the basic properties) */
171
172  int dim=b->dim;
173  float *entrylist=b->valuelist;
174  long ret;
175  long *entryA=_ogg_calloc(entries,sizeof(long));
176  long *entryB=_ogg_calloc(entries,sizeof(long));
177  long entriesA=0;
178  long entriesB=0;
179  long entriesC=0;
180  long pointsA=0;
181  long i,j,k;
182
183  long   besti=-1;
184  long   bestj=-1;
185 
186  char spinbuf[80];
187  sprintf(spinbuf,"splitting [%ld left]... ",totalpoints-*pointsofar);
188
189  /* one reverse index needed */
190  for(i=0;i<b->entries;i++)reventry[i]=-1;
191  for(i=0;i<entries;i++)reventry[entryindex[i]]=i;
192
193  /* We need to find the dividing hyperplane. find the median of each
194     axis as the centerpoint and the normal facing farthest point */
195
196  /* more than one way to do this part.  For small sets, we can brute
197     force it. */
198
199  if(entries<8 || (float)points*entries*entries<16.f*1024*1024){
200    /* try every pair possibility */
201    float best=0;
202    float this;
203    for(i=0;i<entries-1;i++){
204      for(j=i+1;j<entries;j++){
205        spinnit(spinbuf,entries-i);
206        vqsp_count(b->valuelist,pointlist,dim,
207                   membership,reventry,
208                   entryindex,entries, 
209                   pointindex,points,0,
210                   entryA,entryB,
211                   entryindex[i],entryindex[j],
212                   &entriesA,&entriesB,&entriesC);
213        this=(entriesA-entriesC)*(entriesB-entriesC);
214
215        /* when choosing best, we also want some form of stability to
216           make sure more branches are pared later; secondary
217           weighting isn;t needed as the entry lists are in ascending
218           order, and we always try p/q in the same sequence */
219       
220        if( (besti==-1) ||
221            (this>best) ){
222         
223          best=this;
224          besti=entryindex[i];
225          bestj=entryindex[j];
226
227        }
228      }
229    }
230  }else{
231    float *p=alloca(dim*sizeof(float));
232    float *q=alloca(dim*sizeof(float));
233    float best=0.f;
234   
235    /* try COG/normal and furthest pairs */
236    /* meanpoint */
237    /* eventually, we want to select the closest entry and figure n/c
238       from p/q (because storing n/c is too large */
239    for(k=0;k<dim;k++){
240      spinnit(spinbuf,entries);
241     
242      p[k]=0.f;
243      for(j=0;j<entries;j++)
244        p[k]+=b->valuelist[entryindex[j]*dim+k];
245      p[k]/=entries;
246
247    }
248
249    /* we go through the entries one by one, looking for the entry on
250       the other side closest to the point of reflection through the
251       center */
252
253    for(i=0;i<entries;i++){
254      float *ppi=_Nnow(entryindex[i]);
255      float ref_best=0.f;
256      float ref_j=-1;
257      float this;
258      spinnit(spinbuf,entries-i);
259     
260      for(k=0;k<dim;k++)
261        q[k]=2*p[k]-ppi[k];
262
263      for(j=0;j<entries;j++){
264        if(j!=i){
265          float this=_Ndist(dim,q,_Nnow(entryindex[j]));
266          if(ref_j==-1 || this<=ref_best){ /* <=, not <; very important */
267            ref_best=this;
268            ref_j=entryindex[j];
269          }
270        }
271      }
272
273      vqsp_count(b->valuelist,pointlist,dim,
274                 membership,reventry,
275                 entryindex,entries, 
276                 pointindex,points,0,
277                 entryA,entryB,
278                 entryindex[i],ref_j,
279                 &entriesA,&entriesB,&entriesC);
280      this=(entriesA-entriesC)*(entriesB-entriesC);
281
282        /* when choosing best, we also want some form of stability to
283           make sure more branches are pared later; secondary
284           weighting isn;t needed as the entry lists are in ascending
285           order, and we always try p/q in the same sequence */
286       
287      if( (besti==-1) ||
288          (this>best) ){
289       
290        best=this;
291        besti=entryindex[i];
292        bestj=ref_j;
293
294      }
295    }
296    if(besti>bestj){
297      long temp=besti;
298      besti=bestj;
299      bestj=temp;
300    }
301
302  }
303 
304  /* find cells enclosing points */
305  /* count A/B points */
306
307  pointsA=vqsp_count(b->valuelist,pointlist,dim,
308                     membership,reventry,
309                     entryindex,entries, 
310                     pointindex,points,1,
311                     entryA,entryB,
312                     besti,bestj,
313                     &entriesA,&entriesB,&entriesC);
314
315  /*  fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
316      entries,depth,entriesA-entriesC,entriesC,entriesB-entriesC);*/
317  {
318    long thisaux=t->aux++;
319    if(t->aux>=t->alloc){
320      t->alloc*=2;
321      t->ptr0=_ogg_realloc(t->ptr0,sizeof(long)*t->alloc);
322      t->ptr1=_ogg_realloc(t->ptr1,sizeof(long)*t->alloc);
323      t->p=_ogg_realloc(t->p,sizeof(long)*t->alloc);
324      t->q=_ogg_realloc(t->q,sizeof(long)*t->alloc);
325    }
326   
327    t->p[thisaux]=besti;
328    t->q[thisaux]=bestj;
329   
330    if(entriesA==1){
331      ret=1;
332      t->ptr0[thisaux]=entryA[0];
333      *pointsofar+=pointsA;
334    }else{
335      t->ptr0[thisaux]= -t->aux;
336      ret=lp_split(pointlist,totalpoints,b,entryA,entriesA,pointindex,pointsA,
337                   membership,reventry,depth+1,pointsofar); 
338    }
339    if(entriesB==1){
340      ret++;
341      t->ptr1[thisaux]=entryB[0];
342      *pointsofar+=points-pointsA;
343    }else{
344      t->ptr1[thisaux]= -t->aux;
345      ret+=lp_split(pointlist,totalpoints,b,entryB,entriesB,pointindex+pointsA,
346                    points-pointsA,membership,reventry,
347                    depth+1,pointsofar); 
348    }
349  }
350  free(entryA);
351  free(entryB);
352  return(ret);
353}
354
355static int _node_eq(encode_aux_nearestmatch *v, long a, long b){
356  long    Aptr0=v->ptr0[a];
357  long    Aptr1=v->ptr1[a];
358  long    Bptr0=v->ptr0[b];
359  long    Bptr1=v->ptr1[b];
360
361  /* the possibility of choosing the same p and q, but switched, can;t
362     happen because we always look for the best p/q in the same search
363     order and the search is stable */
364
365  if(Aptr0==Bptr0 && Aptr1==Bptr1)
366    return(1);
367
368  return(0);
369}
370
371void vqsp_book(vqgen *v, codebook *b, long *quantlist){
372  long i,j;
373  static_codebook *c=(static_codebook *)b->c;
374  encode_aux_nearestmatch *t;
375
376  memset(b,0,sizeof(codebook));
377  memset(c,0,sizeof(static_codebook));
378  b->c=c;
379  t=c->nearest_tree=_ogg_calloc(1,sizeof(encode_aux_nearestmatch));
380  c->maptype=2;
381
382  /* make sure there are no duplicate entries and that every
383     entry has points */
384
385  for(i=0;i<v->entries;){
386    /* duplicate? if so, eliminate */
387    for(j=0;j<i;j++){
388      if(_Ndist(v->elements,_now(v,i),_now(v,j))==0.f){
389        fprintf(stderr,"found a duplicate entry!  removing...\n");
390        v->entries--;
391        memcpy(_now(v,i),_now(v,v->entries),sizeof(float)*v->elements);
392        memcpy(quantlist+i*v->elements,quantlist+v->entries*v->elements,
393               sizeof(long)*v->elements);
394        break;
395      }
396    }
397    if(j==i)i++;
398  }
399
400  {
401    v->assigned=_ogg_calloc(v->entries,sizeof(long));
402    for(i=0;i<v->points;i++){
403      float *ppt=_point(v,i);
404      float firstmetric=_Ndist(v->elements,_now(v,0),ppt);
405      long   firstentry=0;
406
407      if(!(i&0xff))spinnit("checking... ",v->points-i);
408
409      for(j=0;j<v->entries;j++){
410        float thismetric=_Ndist(v->elements,_now(v,j),ppt);
411        if(thismetric<firstmetric){
412          firstmetric=thismetric;
413          firstentry=j;
414        }
415      }
416     
417      v->assigned[firstentry]++;
418    }
419
420    for(j=0;j<v->entries;){
421      if(v->assigned[j]==0){
422        fprintf(stderr,"found an unused entry!  removing...\n");
423        v->entries--;
424        memcpy(_now(v,j),_now(v,v->entries),sizeof(float)*v->elements);
425        v->assigned[j]=v->assigned[v->elements];
426        memcpy(quantlist+j*v->elements,quantlist+v->entries*v->elements,
427               sizeof(long)*v->elements);
428        continue;
429      }
430      j++;
431    }
432  }
433
434  fprintf(stderr,"Building a book with %ld unique entries...\n",v->entries);
435
436  {
437    long *entryindex=_ogg_malloc(v->entries*sizeof(long *));
438    long *pointindex=_ogg_malloc(v->points*sizeof(long));
439    long *membership=_ogg_malloc(v->points*sizeof(long));
440    long *reventry=_ogg_malloc(v->entries*sizeof(long));
441    long pointssofar=0;
442     
443    for(i=0;i<v->entries;i++)entryindex[i]=i;
444    for(i=0;i<v->points;i++)pointindex[i]=i;
445
446    t->alloc=4096;
447    t->ptr0=_ogg_malloc(sizeof(long)*t->alloc);
448    t->ptr1=_ogg_malloc(sizeof(long)*t->alloc);
449    t->p=_ogg_malloc(sizeof(long)*t->alloc);
450    t->q=_ogg_malloc(sizeof(long)*t->alloc);
451    t->aux=0;
452    c->dim=v->elements;
453    c->entries=v->entries;
454    c->lengthlist=_ogg_calloc(c->entries,sizeof(long));
455    b->valuelist=v->entrylist; /* temporary; replaced later */
456    b->dim=c->dim;
457    b->entries=c->entries;
458
459    for(i=0;i<v->points;i++)membership[i]=-1;
460    for(i=0;i<v->points;i++){
461      float *ppt=_point(v,i);
462      long   firstentry=0;
463      float firstmetric=_Ndist(v->elements,_now(v,0),ppt);
464   
465      if(!(i&0xff))spinnit("assigning... ",v->points-i);
466
467      for(j=1;j<v->entries;j++){
468        if(v->assigned[j]!=-1){
469          float thismetric=_Ndist(v->elements,_now(v,j),ppt);
470          if(thismetric<=firstmetric){
471            firstmetric=thismetric;
472            firstentry=j;
473          }
474        }
475      }
476     
477      membership[i]=firstentry;
478    }
479
480    fprintf(stderr,"Leaves added: %d              \n",
481            lp_split(v->pointlist,v->points,
482                     b,entryindex,v->entries,
483                     pointindex,v->points,
484                     membership,reventry,
485                     0,&pointssofar));
486     
487    free(pointindex);
488    free(membership);
489    free(reventry);
490   
491    fprintf(stderr,"Paring/rerouting redundant branches... ");
492   
493    /* The tree is likely big and redundant.  Pare and reroute branches */
494    {
495      int changedflag=1;
496     
497      while(changedflag){
498        changedflag=0;
499       
500        /* span the tree node by node; list unique decision nodes and
501           short circuit redundant branches */
502       
503        for(i=0;i<t->aux;){
504          int k;
505         
506          /* check list of unique decisions */
507          for(j=0;j<i;j++)
508            if(_node_eq(t,i,j))break;
509         
510          if(j<i){
511            /* a redundant entry; find all higher nodes referencing it and
512               short circuit them to the previously noted unique entry */
513            changedflag=1;
514            for(k=0;k<t->aux;k++){
515              if(t->ptr0[k]==-i)t->ptr0[k]=-j;
516              if(t->ptr1[k]==-i)t->ptr1[k]=-j;
517            }
518           
519            /* Now, we need to fill in the hole from this redundant
520               entry in the listing.  Insert the last entry in the list.
521               Fix the forward pointers to that last entry */
522            t->aux--;
523            t->ptr0[i]=t->ptr0[t->aux];
524            t->ptr1[i]=t->ptr1[t->aux];
525            t->p[i]=t->p[t->aux];
526            t->q[i]=t->q[t->aux];
527            for(k=0;k<t->aux;k++){
528              if(t->ptr0[k]==-t->aux)t->ptr0[k]=-i;
529              if(t->ptr1[k]==-t->aux)t->ptr1[k]=-i;
530            }
531            /* hole plugged */
532           
533          }else
534            i++;
535        }
536       
537        fprintf(stderr,"\rParing/rerouting redundant branches... "
538                "%ld remaining   ",t->aux);
539      }
540      fprintf(stderr,"\n");
541    }
542  }
543 
544  /* run all training points through the decision tree to get a final
545     probability count */
546  {
547    long *probability=_ogg_malloc(c->entries*sizeof(long));
548    for(i=0;i<c->entries;i++)probability[i]=1; /* trivial guard */
549    b->dim=c->dim;
550
551    /* sigh.  A necessary hack */
552    for(i=0;i<t->aux;i++)t->p[i]*=c->dim;
553    for(i=0;i<t->aux;i++)t->q[i]*=c->dim;
554   
555    for(i=0;i<v->points;i++){
556      /* we use the linear matcher regardless becuase the trainer
557         doesn't convert log to linear */
558      int ret=_best(b,v->pointlist+i*v->elements,1);
559      probability[ret]++;
560      if(!(i&0xff))spinnit("counting hits... ",v->points-i);
561    }
562    for(i=0;i<t->aux;i++)t->p[i]/=c->dim;
563    for(i=0;i<t->aux;i++)t->q[i]/=c->dim;
564
565    build_tree_from_lengths(c->entries,probability,c->lengthlist);
566   
567    free(probability);
568  }
569
570  /* Sort the entries by codeword length, short to long (eases
571     assignment and packing to do it now) */
572  {
573    long *wordlen=c->lengthlist;
574    long *index=_ogg_malloc(c->entries*sizeof(long));
575    long *revindex=_ogg_malloc(c->entries*sizeof(long));
576    int k;
577    for(i=0;i<c->entries;i++)index[i]=i;
578    isortvals=c->lengthlist;
579    qsort(index,c->entries,sizeof(long),iascsort);
580
581    /* rearrange storage; ptr0/1 first as it needs a reverse index */
582    /* n and c stay unchanged */
583    for(i=0;i<c->entries;i++)revindex[index[i]]=i;
584    for(i=0;i<t->aux;i++){
585      if(!(i&0x3f))spinnit("sorting... ",t->aux-i);
586
587      if(t->ptr0[i]>=0)t->ptr0[i]=revindex[t->ptr0[i]];
588      if(t->ptr1[i]>=0)t->ptr1[i]=revindex[t->ptr1[i]];
589      t->p[i]=revindex[t->p[i]];
590      t->q[i]=revindex[t->q[i]];
591    }
592    free(revindex);
593
594    /* map lengthlist and vallist with index */
595    c->lengthlist=_ogg_calloc(c->entries,sizeof(long));
596    b->valuelist=_ogg_malloc(sizeof(float)*c->entries*c->dim);
597    c->quantlist=_ogg_malloc(sizeof(long)*c->entries*c->dim);
598    for(i=0;i<c->entries;i++){
599      long e=index[i];
600      for(k=0;k<c->dim;k++){
601        b->valuelist[i*c->dim+k]=v->entrylist[e*c->dim+k];
602        c->quantlist[i*c->dim+k]=quantlist[e*c->dim+k];
603      }
604      c->lengthlist[i]=wordlen[e];
605    }
606
607    free(wordlen);
608  }
609
610  fprintf(stderr,"Done.                            \n\n");
611}
612
Note: See TracBrowser for help on using the repository browser.