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-2000 *
9 * by Monty <monty@xiph.org> and The XIPHOPHORUS Company *
10 * http://www.xiph.org/ *
12 ********************************************************************
14 function: build a VQ codebook and the encoding decision 'tree'
15 last mod: $Id: vqsplit.c,v 1.20 2000/10/12 03:13:02 xiphmont Exp $
17 ********************************************************************/
19 /* This code is *not* part of libvorbis. It is used to generate
20 trained codebooks offline and then spit the results into a
21 pregenerated codebook that is compiled into libvorbis. It is an
22 expensive (but good) algorithm. Run it on big iron. */
24 /* There are so many optimizations to explore in *both* stages that
25 considering the undertaking is almost withering. For now, we brute
37 #include "../lib/sharedbook.h"
39 /* Codebook generation happens in two steps:
41 1) Train the codebook with data collected from the encoder: We use
42 one of a few error metrics (which represent the distance between a
43 given data point and a candidate point in the training set) to
44 divide the training set up into cells representing roughly equal
45 probability of occurring.
47 2) Generate the codebook and auxiliary data from the trained data set
50 /* Building a codebook from trained set **********************************
52 The codebook in raw form is technically finished once it's trained.
53 However, we want to finalize the representative codebook values for
54 each entry and generate auxiliary information to optimize encoding.
55 We generate the auxiliary coding tree using collected data,
56 probably the same data as in the original training */
58 /* At each recursion, the data set is split in half. Cells with data
59 points on side A go into set A, same with set B. The sets may
60 overlap. If the cell overlaps the deviding line only very slightly
61 (provided parameter), we may choose to ignore the overlap in order
62 to pare the tree down */
65 int iascsort(const void *a,const void *b){
66 long av=isortvals[*((long *)a)];
67 long bv=isortvals[*((long *)b)];
71 static float _Ndist(int el,float *a, float *b){
75 float val=(a[i]-b[i]);
81 #define _Npoint(i) (pointlist+dim*(i))
82 #define _Nnow(i) (entrylist+dim*(i))
85 /* goes through the split, but just counts it and returns a metric*/
86 int vqsp_count(float *entrylist,float *pointlist,int dim,
87 long *membership,long *reventry,
88 long *entryindex,long entries,
89 long *pointindex,long points,int splitp,
90 long *entryA,long *entryB,
91 long besti,long bestj,
92 long *entriesA,long *entriesB,long *entriesC){
97 long *temppointsA=NULL;
98 long *temppointsB=NULL;
101 temppointsA=malloc(points*sizeof(long));
102 temppointsB=malloc(points*sizeof(long));
105 memset(entryA,0,sizeof(long)*entries);
106 memset(entryB,0,sizeof(long)*entries);
108 /* Do the points belonging to this cell occur on sideA, sideB or
111 for(i=0;i<points;i++){
112 float *ppt=_Npoint(pointindex[i]);
113 long firstentry=membership[pointindex[i]];
115 if(firstentry==besti){
116 entryA[reventry[firstentry]]=1;
117 if(splitp)temppointsA[pointsA++]=pointindex[i];
120 if(firstentry==bestj){
121 entryB[reventry[firstentry]]=1;
122 if(splitp)temppointsB[pointsB++]=pointindex[i];
126 float distA=_Ndist(dim,ppt,_Nnow(besti));
127 float distB=_Ndist(dim,ppt,_Nnow(bestj));
129 entryA[reventry[firstentry]]=1;
130 if(splitp)temppointsA[pointsA++]=pointindex[i];
132 entryB[reventry[firstentry]]=1;
133 if(splitp)temppointsB[pointsB++]=pointindex[i];
138 /* The entry splitting isn't total, so that storage has to be
139 allocated for recursion. Reuse the entryA/entryB vectors */
140 /* keep the entries in ascending order (relative to the original
141 list); we rely on that stability when ordering p/q choice */
142 for(j=0;j<entries;j++){
143 if(entryA[j] && entryB[j])C++;
144 if(entryA[j])entryA[A++]=entryindex[j];
145 if(entryB[j])entryB[B++]=entryindex[j];
151 memcpy(pointindex,temppointsA,sizeof(long)*pointsA);
152 memcpy(pointindex+pointsA,temppointsB,sizeof(long)*pointsB);
159 int lp_split(float *pointlist,long totalpoints,
161 long *entryindex,long entries,
162 long *pointindex,long points,
163 long *membership,long *reventry,
164 long depth, long *pointsofar){
166 encode_aux_nearestmatch *t=b->c->nearest_tree;
168 /* The encoder, regardless of book, will be using a straight
169 euclidian distance-to-point metric to determine closest point.
170 Thus we split the cells using the same (we've already trained the
171 codebook set spacing and distribution using special metrics and
172 even a midpoint division won't disturb the basic properties) */
175 float *entrylist=b->valuelist;
177 long *entryA=calloc(entries,sizeof(long));
178 long *entryB=calloc(entries,sizeof(long));
189 sprintf(spinbuf,"splitting [%ld left]... ",totalpoints-*pointsofar);
191 /* one reverse index needed */
192 for(i=0;i<b->entries;i++)reventry[i]=-1;
193 for(i=0;i<entries;i++)reventry[entryindex[i]]=i;
195 /* We need to find the dividing hyperplane. find the median of each
196 axis as the centerpoint and the normal facing farthest point */
198 /* more than one way to do this part. For small sets, we can brute
201 if(entries<8 || (float)points*entries*entries<16.*1024*1024){
202 /* try every pair possibility */
205 for(i=0;i<entries-1;i++){
206 for(j=i+1;j<entries;j++){
207 spinnit(spinbuf,entries-i);
208 vqsp_count(b->valuelist,pointlist,dim,
213 entryindex[i],entryindex[j],
214 &entriesA,&entriesB,&entriesC);
215 this=(entriesA-entriesC)*(entriesB-entriesC);
217 /* when choosing best, we also want some form of stability to
218 make sure more branches are pared later; secondary
219 weighting isn;t needed as the entry lists are in ascending
220 order, and we always try p/q in the same sequence */
233 float *p=alloca(dim*sizeof(float));
234 float *q=alloca(dim*sizeof(float));
237 /* try COG/normal and furthest pairs */
239 /* eventually, we want to select the closest entry and figure n/c
240 from p/q (because storing n/c is too large */
242 spinnit(spinbuf,entries);
245 for(j=0;j<entries;j++)
246 p[k]+=b->valuelist[entryindex[j]*dim+k];
251 /* we go through the entries one by one, looking for the entry on
252 the other side closest to the point of reflection through the
255 for(i=0;i<entries;i++){
256 float *ppi=_Nnow(entryindex[i]);
260 spinnit(spinbuf,entries-i);
265 for(j=0;j<entries;j++){
267 float this=_Ndist(dim,q,_Nnow(entryindex[j]));
268 if(ref_j==-1 || this<=ref_best){ /* <=, not <; very important */
275 vqsp_count(b->valuelist,pointlist,dim,
281 &entriesA,&entriesB,&entriesC);
282 this=(entriesA-entriesC)*(entriesB-entriesC);
284 /* when choosing best, we also want some form of stability to
285 make sure more branches are pared later; secondary
286 weighting isn;t needed as the entry lists are in ascending
287 order, and we always try p/q in the same sequence */
306 /* find cells enclosing points */
307 /* count A/B points */
309 pointsA=vqsp_count(b->valuelist,pointlist,dim,
315 &entriesA,&entriesB,&entriesC);
317 /* fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
318 entries,depth,entriesA-entriesC,entriesC,entriesB-entriesC);*/
320 long thisaux=t->aux++;
321 if(t->aux>=t->alloc){
323 t->ptr0=realloc(t->ptr0,sizeof(long)*t->alloc);
324 t->ptr1=realloc(t->ptr1,sizeof(long)*t->alloc);
325 t->p=realloc(t->p,sizeof(long)*t->alloc);
326 t->q=realloc(t->q,sizeof(long)*t->alloc);
334 t->ptr0[thisaux]=entryA[0];
335 *pointsofar+=pointsA;
337 t->ptr0[thisaux]= -t->aux;
338 ret=lp_split(pointlist,totalpoints,b,entryA,entriesA,pointindex,pointsA,
339 membership,reventry,depth+1,pointsofar);
343 t->ptr1[thisaux]=entryB[0];
344 *pointsofar+=points-pointsA;
346 t->ptr1[thisaux]= -t->aux;
347 ret+=lp_split(pointlist,totalpoints,b,entryB,entriesB,pointindex+pointsA,
348 points-pointsA,membership,reventry,
357 static int _node_eq(encode_aux_nearestmatch *v, long a, long b){
358 long Aptr0=v->ptr0[a];
359 long Aptr1=v->ptr1[a];
360 long Bptr0=v->ptr0[b];
361 long Bptr1=v->ptr1[b];
363 /* the possibility of choosing the same p and q, but switched, can;t
364 happen because we always look for the best p/q in the same search
365 order and the search is stable */
367 if(Aptr0==Bptr0 && Aptr1==Bptr1)
373 void vqsp_book(vqgen *v, codebook *b, long *quantlist){
375 static_codebook *c=(static_codebook *)b->c;
376 encode_aux_nearestmatch *t;
378 memset(b,0,sizeof(codebook));
379 memset(c,0,sizeof(static_codebook));
381 t=c->nearest_tree=calloc(1,sizeof(encode_aux_nearestmatch));
384 /* make sure there are no duplicate entries and that every
387 for(i=0;i<v->entries;){
388 /* duplicate? if so, eliminate */
390 if(_Ndist(v->elements,_now(v,i),_now(v,j))==0.){
391 fprintf(stderr,"found a duplicate entry! removing...\n");
393 memcpy(_now(v,i),_now(v,v->entries),sizeof(float)*v->elements);
394 memcpy(quantlist+i*v->elements,quantlist+v->entries*v->elements,
395 sizeof(long)*v->elements);
403 v->assigned=calloc(v->entries,sizeof(long));
404 for(i=0;i<v->points;i++){
405 float *ppt=_point(v,i);
406 float firstmetric=_Ndist(v->elements,_now(v,0),ppt);
409 if(!(i&0xff))spinnit("checking... ",v->points-i);
411 for(j=0;j<v->entries;j++){
412 float thismetric=_Ndist(v->elements,_now(v,j),ppt);
413 if(thismetric<firstmetric){
414 firstmetric=thismetric;
419 v->assigned[firstentry]++;
422 for(j=0;j<v->entries;){
423 if(v->assigned[j]==0){
424 fprintf(stderr,"found an unused entry! removing...\n");
426 memcpy(_now(v,j),_now(v,v->entries),sizeof(float)*v->elements);
427 v->assigned[j]=v->assigned[v->elements];
428 memcpy(quantlist+j*v->elements,quantlist+v->entries*v->elements,
429 sizeof(long)*v->elements);
436 fprintf(stderr,"Building a book with %ld unique entries...\n",v->entries);
439 long *entryindex=malloc(v->entries*sizeof(long *));
440 long *pointindex=malloc(v->points*sizeof(long));
441 long *membership=malloc(v->points*sizeof(long));
442 long *reventry=malloc(v->entries*sizeof(long));
445 for(i=0;i<v->entries;i++)entryindex[i]=i;
446 for(i=0;i<v->points;i++)pointindex[i]=i;
449 t->ptr0=malloc(sizeof(long)*t->alloc);
450 t->ptr1=malloc(sizeof(long)*t->alloc);
451 t->p=malloc(sizeof(long)*t->alloc);
452 t->q=malloc(sizeof(long)*t->alloc);
455 c->entries=v->entries;
456 c->lengthlist=calloc(c->entries,sizeof(long));
457 b->valuelist=v->entrylist; /* temporary; replaced later */
459 b->entries=c->entries;
461 for(i=0;i<v->points;i++)membership[i]=-1;
462 for(i=0;i<v->points;i++){
463 float *ppt=_point(v,i);
465 float firstmetric=_Ndist(v->elements,_now(v,0),ppt);
467 if(!(i&0xff))spinnit("assigning... ",v->points-i);
469 for(j=1;j<v->entries;j++){
470 if(v->assigned[j]!=-1){
471 float thismetric=_Ndist(v->elements,_now(v,j),ppt);
472 if(thismetric<=firstmetric){
473 firstmetric=thismetric;
479 membership[i]=firstentry;
482 fprintf(stderr,"Leaves added: %d \n",
483 lp_split(v->pointlist,v->points,
484 b,entryindex,v->entries,
485 pointindex,v->points,
493 fprintf(stderr,"Paring/rerouting redundant branches... ");
495 /* The tree is likely big and redundant. Pare and reroute branches */
502 /* span the tree node by node; list unique decision nodes and
503 short circuit redundant branches */
508 /* check list of unique decisions */
510 if(_node_eq(t,i,j))break;
513 /* a redundant entry; find all higher nodes referencing it and
514 short circuit them to the previously noted unique entry */
516 for(k=0;k<t->aux;k++){
517 if(t->ptr0[k]==-i)t->ptr0[k]=-j;
518 if(t->ptr1[k]==-i)t->ptr1[k]=-j;
521 /* Now, we need to fill in the hole from this redundant
522 entry in the listing. Insert the last entry in the list.
523 Fix the forward pointers to that last entry */
525 t->ptr0[i]=t->ptr0[t->aux];
526 t->ptr1[i]=t->ptr1[t->aux];
527 t->p[i]=t->p[t->aux];
528 t->q[i]=t->q[t->aux];
529 for(k=0;k<t->aux;k++){
530 if(t->ptr0[k]==-t->aux)t->ptr0[k]=-i;
531 if(t->ptr1[k]==-t->aux)t->ptr1[k]=-i;
539 fprintf(stderr,"\rParing/rerouting redundant branches... "
540 "%ld remaining ",t->aux);
542 fprintf(stderr,"\n");
546 /* run all training points through the decision tree to get a final
549 long *probability=malloc(c->entries*sizeof(long));
550 for(i=0;i<c->entries;i++)probability[i]=1; /* trivial guard */
553 /* sigh. A necessary hack */
554 for(i=0;i<t->aux;i++)t->p[i]*=c->dim;
555 for(i=0;i<t->aux;i++)t->q[i]*=c->dim;
557 for(i=0;i<v->points;i++){
558 /* we use the linear matcher regardless becuase the trainer
559 doesn't convert log to linear */
560 int ret=_best(b,v->pointlist+i*v->elements,1);
562 if(!(i&0xff))spinnit("counting hits... ",v->points-i);
564 for(i=0;i<t->aux;i++)t->p[i]/=c->dim;
565 for(i=0;i<t->aux;i++)t->q[i]/=c->dim;
567 build_tree_from_lengths(c->entries,probability,c->lengthlist);
572 /* Sort the entries by codeword length, short to long (eases
573 assignment and packing to do it now) */
575 long *wordlen=c->lengthlist;
576 long *index=malloc(c->entries*sizeof(long));
577 long *revindex=malloc(c->entries*sizeof(long));
579 for(i=0;i<c->entries;i++)index[i]=i;
580 isortvals=c->lengthlist;
581 qsort(index,c->entries,sizeof(long),iascsort);
583 /* rearrange storage; ptr0/1 first as it needs a reverse index */
584 /* n and c stay unchanged */
585 for(i=0;i<c->entries;i++)revindex[index[i]]=i;
586 for(i=0;i<t->aux;i++){
587 if(!(i&0x3f))spinnit("sorting... ",t->aux-i);
589 if(t->ptr0[i]>=0)t->ptr0[i]=revindex[t->ptr0[i]];
590 if(t->ptr1[i]>=0)t->ptr1[i]=revindex[t->ptr1[i]];
591 t->p[i]=revindex[t->p[i]];
592 t->q[i]=revindex[t->q[i]];
596 /* map lengthlist and vallist with index */
597 c->lengthlist=calloc(c->entries,sizeof(long));
598 b->valuelist=malloc(sizeof(float)*c->entries*c->dim);
599 c->quantlist=malloc(sizeof(long)*c->entries*c->dim);
600 for(i=0;i<c->entries;i++){
602 for(k=0;k<c->dim;k++){
603 b->valuelist[i*c->dim+k]=v->entrylist[e*c->dim+k];
604 c->quantlist[i*c->dim+k]=quantlist[e*c->dim+k];
606 c->lengthlist[i]=wordlen[e];
612 fprintf(stderr,"Done. \n\n");