1 /********************************************************************
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. *
8 * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999 *
9 * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company *
10 * http://www.xiph.org/ *
12 ********************************************************************
14 function: build a VQ codebook
15 author: Monty <xiphmont@mit.edu>
16 modifications by: Monty
17 last modification date: Dec 10 1999
19 ********************************************************************/
21 /* This code is *not* part of libvorbis. It is used to generate
22 trained codebooks offline and then spit the results into a
23 pregenerated codebook that is compiled into libvorbis. It is an
24 expensive (but good) algorithm. Run it on big iron. */
26 /* There are so many optimizations to explore in *both* stages that
27 considering the undertaking is almost withering. For now, we brute
37 /* Codebook generation happens in two steps:
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.
45 2) Generate the codebook and auxiliary data from the trained data set
48 /* Building a codebook from trained set **********************************
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 */
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 */
63 int iascsort(const void *a,const void *b){
64 long av=isortvals[*((long *)a)];
65 long bv=isortvals[*((long *)b)];
69 extern double _dist_sq(vqgen *v,double *a, double *b);
71 /* goes through the split, but just counts it and returns a metric*/
72 void vqsp_count(vqgen *v,long *membership,
73 long *entryindex,long entries,
74 long *pointindex,long points,
75 long *entryA,long *entryB,
77 long *entriesA,long *entriesB,long *entriesC){
81 memset(entryA,0,sizeof(long)*entries);
82 memset(entryB,0,sizeof(long)*entries);
84 /* Do the points belonging to this cell occur on sideA, sideB or
87 for(i=0;i<points;i++){
88 double *ppt=_point(v,pointindex[i]);
89 long firstentry=membership[i];
92 if(!entryA[firstentry] || !entryB[firstentry]){
93 for(k=0;k<v->elements;k++)
94 position+=ppt[k]*n[k];
102 /* The entry splitting isn't total, so that storage has to be
103 allocated for recursion. Reuse the entryA/entryB vectors */
104 /* keep the entries in ascending order (relative to the original
105 list); we rely on that stability when ordering p/q choice */
106 for(j=0;j<entries;j++){
107 if(entryA[j] && entryB[j])C++;
108 if(entryA[j])entryA[A++]=entryindex[j];
109 if(entryB[j])entryB[B++]=entryindex[j];
116 void pq_in_out(vqgen *v,double *n,double *c,double *p,double *q){
119 for(k=0;k<v->elements;k++){
120 double center=(p[k]+q[k])/2.;
121 n[k]=(center-q[k])*2.;
126 static void spinnit(void){
128 static long lasttime=0;
130 struct timeval thistime;
132 gettimeofday(&thistime,NULL);
133 test=thistime.tv_sec*10+thistime.tv_usec/100000;
140 fprintf(stderr,"|\b");
143 fprintf(stderr,"/\b");
146 fprintf(stderr,"-\b");
149 fprintf(stderr,"\\\b");
156 int lp_split(vqgen *v,vqbook *b,
157 long *entryindex,long entries,
158 long *pointindex,long points,
161 /* The encoder, regardless of book, will be using a straight
162 euclidian distance-to-point metric to determine closest point.
163 Thus we split the cells using the same (we've already trained the
164 codebook set spacing and distribution using special metrics and
165 even a midpoint division won't disturb the basic properties) */
168 double *n=alloca(sizeof(double)*v->elements);
170 long *entryA=calloc(entries,sizeof(long));
171 long *entryB=calloc(entries,sizeof(long));
178 long *membership=malloc(sizeof(long)*points);
183 /* which cells do points belong to? Do this before n^2 best pair chooser. */
185 for(i=0;i<points;i++){
186 double *ppt=_point(v,pointindex[i]);
188 double firstmetric=_dist_sq(v,_now(v,entryindex[0]),ppt);
190 if(points*entries>64*1024)spinnit();
192 for(j=1;j<entries;j++){
193 double thismetric=_dist_sq(v,_now(v,entryindex[j]),ppt);
194 if(thismetric<firstmetric){
195 firstmetric=thismetric;
200 membership[i]=firstentry;
203 /* We need to find the dividing hyperplane. find the median of each
204 axis as the centerpoint and the normal facing farthest point */
206 /* more than one way to do this part. For small sets, we can brute
209 if(entries<8 || points*entries*entries<128*1024*1024){
210 /* try every pair possibility */
213 for(i=0;i<entries-1;i++){
214 for(j=i+1;j<entries;j++){
216 pq_in_out(v,n,&c,_now(v,entryindex[i]),_now(v,entryindex[j]));
217 vqsp_count(v,membership,
222 &entriesA,&entriesB,&entriesC);
223 this=(entriesA-entriesC)*(entriesB-entriesC);
225 /* when choosing best, we also want some form of stability to
226 make sure more branches are pared later; secondary
227 weighting isn;t needed as the entry lists are in ascending
228 order, and we always try p/q in the same sequence */
241 double *p=alloca(v->elements*sizeof(double));
242 double *q=alloca(v->elements*sizeof(double));
245 /* try COG/normal and furthest pairs */
247 /* eventually, we want to select the closest entry and figure n/c
248 from p/q (because storing n/c is too large */
249 for(k=0;k<v->elements;k++){
253 for(j=0;j<entries;j++)
254 p[k]+=v->entrylist[entryindex[j]*v->elements+k];
259 /* we go through the entries one by one, looking for the entry on
260 the other side closest to the point of reflection through the
263 for(i=0;i<entries;i++){
264 double *ppi=_now(v,entryindex[i]);
270 for(k=0;k<v->elements;k++)
273 for(j=0;j<entries;j++){
275 double this=_dist_sq(v,q,_now(v,entryindex[j]));
276 if(ref_j==-1 || this<ref_best){
283 pq_in_out(v,n,&c,ppi,_now(v,ref_j));
284 vqsp_count(v,membership,
289 &entriesA,&entriesB,&entriesC);
290 this=(entriesA-entriesC)*(entriesB-entriesC);
292 /* when choosing best, we also want some form of stability to
293 make sure more branches are pared later; secondary
294 weighting isn;t needed as the entry lists are in ascending
295 order, and we always try p/q in the same sequence */
314 /* find cells enclosing points */
315 /* count A/B points */
317 pq_in_out(v,n,&c,_now(v,besti),_now(v,bestj));
318 vqsp_count(v,membership,
323 &entriesA,&entriesB,&entriesC);
327 /* the point index is split, so we do an Order n rearrangement into
328 A first/B last and just pass it on */
335 for(k=0;k<v->elements;k++)
336 position+=_point(v,pointindex[Aptr])[k]*n[k];
337 if(position<=0.)break; /* not in A */
342 for(k=0;k<v->elements;k++)
343 position+=_point(v,pointindex[Bptr])[k]*n[k];
344 if(position>0.)break; /* not in B */
348 long temp=pointindex[Aptr];
349 pointindex[Aptr]=pointindex[Bptr];
350 pointindex[Bptr]=temp;
356 fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
357 entries,depth,entriesA-entriesC,entriesC,entriesB-entriesC);
359 long thisaux=b->aux++;
360 if(b->aux>=b->alloc){
362 b->ptr0=realloc(b->ptr0,sizeof(long)*b->alloc);
363 b->ptr1=realloc(b->ptr1,sizeof(long)*b->alloc);
364 b->p=realloc(b->p,sizeof(long)*b->alloc);
365 b->q=realloc(b->q,sizeof(long)*b->alloc);
373 b->ptr0[thisaux]=entryA[0];
375 b->ptr0[thisaux]= -b->aux;
376 ret=lp_split(v,b,entryA,entriesA,pointindex,pointsA,depth+1);
380 b->ptr1[thisaux]=entryB[0];
382 b->ptr1[thisaux]= -b->aux;
383 ret+=lp_split(v,b,entryB,entriesB,pointindex+pointsA,points-pointsA,
392 static int _node_eq(vqbook *v, long a, long b){
393 long Aptr0=v->ptr0[a];
394 long Aptr1=v->ptr1[a];
397 long Bptr0=v->ptr0[b];
398 long Bptr1=v->ptr1[b];
403 /* the possibility of choosing the same p and q, but switched, can;t
404 happen because we always look for the best p/q in the same search
405 order and the search is stable */
407 if(Aptr0==Bptr0 && Aptr1==Bptr1 && Ap==Bp && Aq==Bq)
413 void vqsp_book(vqgen *v, vqbook *b, long *quantlist){
414 long *entryindex=malloc(sizeof(long)*v->entries);
415 long *pointindex=malloc(sizeof(long)*v->points);
416 long *membership=malloc(sizeof(long)*v->points);
419 memset(b,0,sizeof(vqbook));
420 for(i=0;i<v->entries;i++)entryindex[i]=i;
421 for(i=0;i<v->points;i++)pointindex[i]=i;
424 b->entries=v->entries;
427 b->ptr0=malloc(sizeof(long)*b->alloc);
428 b->ptr1=malloc(sizeof(long)*b->alloc);
429 b->p=malloc(sizeof(long)*b->alloc);
430 b->q=malloc(sizeof(long)*b->alloc);
431 b->lengthlist=calloc(b->entries,sizeof(long));
433 /* which cells do points belong to? Only do this once. */
435 for(i=0;i<v->points;i++){
436 double *ppt=_point(v,i);
438 double firstmetric=_dist_sq(v,_now(v,0),ppt);
440 for(j=1;j<v->entries;j++){
441 double thismetric=_dist_sq(v,_now(v,j),ppt);
442 if(thismetric<firstmetric){
443 firstmetric=thismetric;
448 membership[i]=firstentry;
452 /* first, generate the encoding decision heirarchy */
453 fprintf(stderr,"Total leaves: %d\n",
454 lp_split(v,b,entryindex,v->entries,
455 pointindex,v->points,0));
460 fprintf(stderr,"Paring and rerouting redundant branches:\n");
461 /* The tree is likely big and redundant. Pare and reroute branches */
464 long *unique_entries=alloca(b->aux*sizeof(long));
470 fprintf(stderr,"\t...");
471 /* span the tree node by node; list unique decision nodes and
472 short circuit redundant branches */
477 /* check list of unique decisions */
479 if(_node_eq(b,i,unique_entries[j]))break;
483 unique_entries[nodes++]=i++;
485 /* a redundant entry; find all higher nodes referencing it and
486 short circuit them to the previously noted unique entry */
488 for(k=0;k<b->aux;k++){
489 if(b->ptr0[k]==-i)b->ptr0[k]=-unique_entries[j];
490 if(b->ptr1[k]==-i)b->ptr1[k]=-unique_entries[j];
493 /* Now, we need to fill in the hole from this redundant
494 entry in the listing. Insert the last entry in the list.
495 Fix the forward pointers to that last entry */
497 b->ptr0[i]=b->ptr0[b->aux];
498 b->ptr1[i]=b->ptr1[b->aux];
499 b->p[i]=b->p[b->aux];
500 b->q[i]=b->q[b->aux];
501 for(k=0;k<b->aux;k++){
502 if(b->ptr0[k]==-b->aux)b->ptr0[k]=-i;
503 if(b->ptr1[k]==-b->aux)b->ptr1[k]=-i;
510 fprintf(stderr," %ld remaining\n",b->aux);
514 /* run all training points through the decision tree to get a final
517 long *probability=malloc(b->entries*sizeof(long));
518 long *membership=malloc(b->entries*sizeof(long));
520 for(i=0;i<b->entries;i++)probability[i]=1; /* trivial guard */
522 b->valuelist=v->entrylist; /* temporary for vqenc_entry; replaced later */
523 for(i=0;i<v->points;i++){
524 int ret=vqenc_entry(b,v->pointlist+i*v->elements);
527 for(i=0;i<v->entries;i++)membership[i]=i;
529 /* find codeword lengths */
530 /* much more elegant means exist. Brute force n^2, minimum thought */
531 for(i=v->entries;i>1;i--){
532 int first=-1,second=-1;
533 long least=v->points+1;
535 /* find the two nodes to join */
536 for(j=0;j<v->entries;j++)
537 if(probability[j]<least){
538 least=probability[j];
542 for(j=0;j<v->entries;j++)
543 if(probability[j]<least && membership[j]!=first){
544 least=probability[j];
545 second=membership[j];
547 if(first==-1 || second==-1){
548 fprintf(stderr,"huffman fault; no free branch\n");
553 least=probability[first]+probability[second];
554 for(j=0;j<v->entries;j++)
555 if(membership[j]==first || membership[j]==second){
557 probability[j]=least;
561 for(i=0;i<v->entries-1;i++)
562 if(membership[i]!=membership[i+1]){
563 fprintf(stderr,"huffman fault; failed to build single tree\n");
567 /* unneccessary metric */
568 memset(probability,0,sizeof(long)*v->entries);
569 for(i=0;i<v->points;i++){
570 int ret=vqenc_entry(b,v->pointlist+i*v->elements);
578 /* Sort the entries by codeword length, short to long (eases
579 assignment and packing to do it now) */
581 long *wordlen=b->lengthlist;
582 long *index=malloc(b->entries*sizeof(long));
583 long *revindex=malloc(b->entries*sizeof(long));
585 for(i=0;i<b->entries;i++)index[i]=i;
586 isortvals=b->lengthlist;
587 qsort(index,b->entries,sizeof(long),iascsort);
589 /* rearrange storage; ptr0/1 first as it needs a reverse index */
590 /* n and c stay unchanged */
591 for(i=0;i<b->entries;i++)revindex[index[i]]=i;
592 for(i=0;i<b->aux;i++){
593 if(b->ptr0[i]>=0)b->ptr0[i]=revindex[b->ptr0[i]];
594 if(b->ptr1[i]>=0)b->ptr1[i]=revindex[b->ptr1[i]];
595 b->p[i]=revindex[b->p[i]];
596 b->q[i]=revindex[b->q[i]];
600 /* map lengthlist and vallist with index */
601 b->lengthlist=calloc(b->entries,sizeof(long));
602 b->valuelist=malloc(sizeof(double)*b->entries*b->dim);
603 b->quantlist=malloc(sizeof(long)*b->entries*b->dim);
604 for(i=0;i<b->entries;i++){
606 for(k=0;k<b->dim;k++){
607 b->valuelist[i*b->dim+k]=v->entrylist[e*b->dim+k];
608 b->quantlist[i*b->dim+k]=quantlist[e*b->dim+k];
610 b->lengthlist[i]=wordlen[e];
616 /* generate the codewords (short to long) */
620 b->codelist=malloc(sizeof(long)*b->entries);
621 for(i=0;i<b->entries;i++){
622 if(length != b->lengthlist[i]){
623 current<<=(b->lengthlist[i]-length);
624 length=b->lengthlist[i];
626 b->codelist[i]=current;
631 /* sanity check the codewords */
632 for(i=0;i<b->entries;i++){
633 for(j=i+1;j<b->entries;j++){
634 if(b->codelist[i]==b->codelist[j]){
635 fprintf(stderr,"Error; codewords for %ld and %ld both equal %lx\n",
636 i, j, b->codelist[i]);
641 fprintf(stderr,"Done.\n\n");
644 /* slow version for use here that does not use a preexpanded n/c. */
645 int vqenc_entry(vqbook *b,double *val){
647 double *n=alloca(b->dim*sizeof(double));
651 double *p=b->valuelist+b->p[ptr]*b->dim;
652 double *q=b->valuelist+b->q[ptr]*b->dim;
654 for(k=0;k<b->dim;k++){
660 for(k=0;k<b->dim;k++)