From 70e7a54f8508c8775daa93b09b76ecfdfc49ba2b Mon Sep 17 00:00:00 2001 From: Monty Date: Thu, 9 Dec 1999 10:16:13 +0000 Subject: [PATCH] New split code. Monty svn path=/trunk/vorbis/; revision=183 --- vq/vqgen.c | 1018 +++++++++++++++++++++++++++++++----------------------------- 1 file changed, 524 insertions(+), 494 deletions(-) diff --git a/vq/vqgen.c b/vq/vqgen.c index 06f4a71..21c3b00 100644 --- a/vq/vqgen.c +++ b/vq/vqgen.c @@ -23,13 +23,29 @@ pregenerated codebook that is compiled into libvorbis. It is an expensive (but good) algorithm. Run it on big iron. */ +/* There are so many optimizations to explore in *both* stages that + considering the undertaking is almost withering. For now, we brute + force it all */ + #include #include #include #include #include "minit.h" -/************************************************************************ +/* Codebook generation happens in two steps: + + 1) Train the codebook with data collected from the encoder: We use + one of a few error metrics (which represent the distance between a + given data point and a candidate point in the training set) to + divide the training set up into cells representing roughly equal + probability of occurring. + + 2) Generate the codebook and auxiliary data from the trained data set +*/ + +/* Codebook training **************************************************** + * * The basic idea here is that a VQ codebook is like an m-dimensional * foam with n bubbles. The bubbles compete for space/volume and are * 'pressurized' [biased] according to some metric. The basic alg @@ -88,6 +104,8 @@ double _dist_sq(vqgen *v,double *a, double *b){ } return acc; } + +/* A metric for LSP codes */ /* candidate,actual */ double _dist_and_pos(vqgen *v,double *b, double *a){ int i; @@ -163,7 +181,6 @@ double vqgen_iterate(vqgen *v){ double *new=malloc(sizeof(double)*v->entries*v->elements); long *nearcount=malloc(v->entries*sizeof(long)); double *nearbias=malloc(v->entries*desired*sizeof(double)); - long biascount; #ifdef NOISY char buff[80]; @@ -178,7 +195,7 @@ double vqgen_iterate(vqgen *v){ bias=fopen(buff,"w"); #endif - fprintf(stderr,"Pass #%ld... ",v->it); + fprintf(stderr,"Pass #%d... ",v->it); if(v->entries<2){ fprintf(stderr,"generation requires at least two entries\n"); @@ -190,8 +207,9 @@ double vqgen_iterate(vqgen *v){ memset(nearcount,0,sizeof(long)*v->entries); memset(v->assigned,0,sizeof(long)*v->entries); for(i=0;ipoints;i++){ - double firstmetric=v->metric_func(v,_now(v,0),_point(v,i))+v->bias[0]; - double secondmetric=v->metric_func(v,_now(v,1),_point(v,i))+v->bias[1]; + double *ppt=_point(v,i); + double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0]; + double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1]; long firstentry=0; long secondentry=1; if(firstmetric>secondmetric){ @@ -304,228 +322,260 @@ double vqgen_iterate(vqgen *v){ return(asserror); } -/* the additional fields are the hyperplanes we've already used to - split the set. The additional hyperplanes are necessary to bound - later splits. One per depth */ - -int lp_split(vqgen *v,long *index,long points, - double *additional_a, double *additional_b, long depth){ - static int frameno=0; - long A[points]; - long B[points]; - long ap=0; - long bp=0; - long cp=0; +/* Building a codebook from trained set ********************************** + + The codebook in raw form is technically finished once it's trained. + However, we want to finalize the representative codebook values for + each entry and generate auxiliary information to optimize encoding. + We generate the auxiliary coding tree using collected data, + probably the same data as in the original training */ + +/* At each recursion, the data set is split in half. Cells with data + points on side A go into set A, same with set B. The sets may + overlap. If the cell overlaps the deviding line only very slightly + (provided parameter), we may choose to ignore the overlap in order + to pare the tree down */ + +double *sortvals; +int els; +int iascsort(const void *a,const void *b){ + double av=sortvals[*((long *)a) * els]; + double bv=sortvals[*((long *)b) * els]; + if(av7){ - printf("leaf: points %ld, depth %ld\n",points,depth); - return(1); + for(i=0;ielements;k++) + position+=ppt[k]*n[k]; + if(position>0.){ + entryA[firstentry]++; + }else{ + entryB[firstentry]++; + } + } + + /* look to see if entries are in the slack zone */ + /* The entry splitting isn't total, so that storage has to be + allocated for recursion. Reuse the entryA/entryB vectors */ + for(j=0;jelements;k++){ + double center=(p[k]+q[k])/2.; + n[k]=(center-q[k])*2.; + *c+=center*n[k]; + } +} - if(points==1){ - printf("leaf: point %ld, depth %ld\n",index[0],depth); +void pq_center_out(vqgen *v,double *n,double *c,double *center,double *q){ + int k; + *c=0.; + for(k=0;kelements;k++){ + n[k]=(center[k]-q[k])*2.; + *c+=center[k]*n[k]; + } +} + +int lp_split(vqgen *v,long *entryindex,long entries, + long *pointindex,long points,long depth, + double slack){ + + /* The encoder, regardless of book, will be using a straight + euclidian distance-to-point metric to determine closest point. + Thus we split the cells using the same (we've already trained the + codebook set spacing and distribution using special metrics and + even a midpoint division won't disturb the basic properties) */ + + long ret; + double *p; + double *q; + double *n; + double c; + long *entryA=calloc(entries,sizeof(long)); + long *entryB=calloc(entries,sizeof(long)); + long entriesA=0; + long entriesB=0; + long entriesC=0; + long pointsA=0; + long i,j,k; + + p=alloca(sizeof(double)*v->elements); + q=alloca(sizeof(double)*v->elements); + n=alloca(sizeof(double)*v->elements); + memset(p,0,sizeof(double)*v->elements); + + /* depth limit */ + if(depth>20){ + printf("leaf: entries %ld, depth %ld\n",entries,depth); + return(entries); + } + + /* nothing to do */ + if(entries==1){ + printf("leaf: entry %ld, depth %ld\n",entryindex[0],depth); return(1); } - if(points==2){ - /* The result must be an even split */ - B[bp++]=index[0]; - A[ap++]=index[1]; + /* The result must be an even split */ + if(entries==2){ printf("even split: depth %ld\n",depth); return(2); - }else{ - /* We need to find the best split */ - for(i=0;ibest){ + } + + /* We need to find the dividing hyperplane. find the median of each + axis as the centerpoint and the normal facing farthest point */ + + /* more than one way to do this part. For small sets, we can brute + force it. */ + + if(entries<64){ + /* try every pair possibility */ + double best=0; + long besti=0; + long bestj=0; + double this; + for(i=0;ibest){ + best=this; besti=i; bestj=j; - best=_dist_sq(v,_now(v,index[i]),_now(v,index[j])); } } } - - /* We have our endpoints. initial divvy */ - { - double *cA=malloc(v->elements*sizeof(double)); - double *cB=malloc(v->elements*sizeof(double)); - double **a=malloc((points+depth)*sizeof(double *)); - double *b=malloc((points+depth)*sizeof(double)); - double *aa=malloc((points+depth)*v->elements*sizeof(double)); - double dA=0.,dB=0.; - minit l; - - minit_init(&l,v->elements,points+depth-1,0); - - /* divisor hyperplane */ - for(i=0;ielements;i++){ - double m=(_now(v,index[besti])[i]+_now(v,index[bestj])[i])/2; - cA[i]=_now(v,index[besti])[i]-_now(v,index[bestj])[i]; - cB[i]=-cA[i]; - dA+=cA[i]*m; - dB-=cA[i]*m; + pq_in_out(v,n,&c,_now(v,entryindex[besti]),_now(v,entryindex[bestj])); + }else{ + double best=0.; + long bestj=0; + + /* try COG/normal and furthest pairs */ + /* medianpoint */ + for(k=0;kelements;k++){ + /* just sort the index array */ + sortvals=v->pointlist+k; + els=v->elements; + qsort(pointindex,points,sizeof(long),iascsort); + if(points&0x1){ + p[k]=v->pointlist[(pointindex[points/2])*v->elements+k]; + }else{ + p[k]=(v->pointlist[(pointindex[points/2])*v->elements+k]+ + v->pointlist[(pointindex[points/2-1])*v->elements+k])/2.; } - - for(i=0;ielements; - - /* check each bubble to see if it intersects the divisor hyperplane */ - for(j=0;jelements;k++){ - double m=(_now(v,index[i])[k]+_now(v,index[j])[k])/2; - a[count][k]=_now(v,index[i])[k]-_now(v,index[j])[k]; - b[count]+=a[count][k]*m; - } - count++; - } - } - - /* additional bounding hyperplanes from previous splits */ - for(i=0;ielements;k++) - a[count][k]=additional_a[m++]; - b[count++]=additional_b[i]; - } - - /* on what side of the dividing hyperplane is the current test - point? */ - for(i=0;ielements;i++) - d1+=_now(v,index[j])[i]*cA[i]; - - - if(d1dA){cp++;B[bp++]=index[j];} - }else{ - B[bp++]=index[j]; - if(d2best){ + best=this; + bestj=j; } + } - /*{ - char buffer[80]; - FILE *o; - int i; - - sprintf(buffer,"set%d.m",frameno); - o=fopen(buffer,"w"); - for(i=0;ielements*(depth+1)); - double *newadd_b=alloca(sizeof(double)*(depth+1)); - if(additional_a){ - memcpy(newadd_a,additional_a,sizeof(double)*v->elements*depth); - memcpy(newadd_b,additional_b,sizeof(double)*depth); - } + pq_center_out(v,n,&c,p,_point(v,pointindex[bestj])); - memcpy(newadd_a+v->elements*depth,cA,sizeof(double)*v->elements); - newadd_b[depth]=dA; - leaves=lp_split(v,A,ap,newadd_a,newadd_b,depth+1); - memcpy(newadd_a+v->elements*depth,cB,sizeof(double)*v->elements); - newadd_b[depth]=dB; - leaves+=lp_split(v,B,bp,newadd_a,newadd_b,depth+1); - + } + + /* find cells enclosing points */ + /* count A/B points */ + + lp_count(v,entryindex,entries, + pointindex,points, + entryA,entryB, + n, c, slack, + &entriesA,&entriesB,&entriesC); + + /* the point index is split evenly, so we do an Order n + rearrangement into A first/B last and just pass it on */ + { + long Aptr=0; + long Bptr=points-1; + while(Aptrelements;k++) + position+=_point(v,pointindex[Aptr])[k]*n[k]; + if(position<0.)break; /* not in A */ + Aptr++; + } + while(Aptrelements;k++) + position+=_point(v,pointindex[Bptr])[k]*n[k]; + if(position>0.)break; /* not in B */ + Bptr--; + } + if(Aptr