Update documentation, version numbers, copyright dates
[platform/upstream/libvorbis.git] / vq / vqgen.c
index 74c83ea..90e08a2 100644 (file)
@@ -1,18 +1,17 @@
 /********************************************************************
  *                                                                  *
- * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
- * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
- * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
- * PLEASE READ THESE TERMS DISTRIBUTING.                            *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
  *                                                                  *
- * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000             *
- * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
- * http://www.xiph.org/                                             *
+ * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001             *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
  *                                                                  *
  ********************************************************************
 
  function: train a VQ codebook 
- last mod: $Id: vqgen.c,v 1.31 2000/05/08 20:49:51 xiphmont Exp $
+ last mod: $Id$
 
  ********************************************************************/
 
@@ -32,7 +31,6 @@
 
 #include "vqgen.h"
 #include "bookutil.h"
-#include "../lib/sharedbook.h"
 
 /* Codebook generation happens in two steps: 
 
 #define vN(data,i) (data+v->elements*i)
 
 /* default metric; squared 'distance' from desired value. */
-double _dist(vqgen *v,double *a, double *b){
+float _dist(vqgen *v,float *a, float *b){
   int i;
   int el=v->elements;
-  double acc=0.;
+  float acc=0.f;
   for(i=0;i<el;i++){
-    double val=(a[i]-b[i]);
+    float val=(a[i]-b[i]);
     acc+=val*val;
   }
   return sqrt(acc);
 }
 
-double *_weight_null(vqgen *v,double *a){
+float *_weight_null(vqgen *v,float *a){
   return a;
 }
 
@@ -78,24 +76,24 @@ double *_weight_null(vqgen *v,double *a){
 void _vqgen_seed(vqgen *v){
   long i;
   for(i=0;i<v->entries;i++)
-    memcpy(_now(v,i),_point(v,i),sizeof(double)*v->elements);
+    memcpy(_now(v,i),_point(v,i),sizeof(float)*v->elements);
+  v->seeded=1;
 }
 
 int directdsort(const void *a, const void *b){
-  double av=*((double *)a);
-  double bv=*((double *)b);
-  if(av>bv)return(-1);
-  return(1);
+  float av=*((float *)a);
+  float bv=*((float *)b);
+  return (av<bv)-(av>bv);
 }
 
 void vqgen_cellmetric(vqgen *v){
   int j,k;
-  double min=-1.,max=-1.,mean=0.,acc=0.;
+  float min=-1.f,max=-1.f,mean=0.f,acc=0.f;
   long dup=0,unused=0;
  #ifdef NOISY
   int i;
    char buff[80];
-   double spacings[v->entries];
+   float spacings[v->entries];
    int count=0;
    FILE *cells;
    sprintf(buff,"cellspace%d.m",v->it);
@@ -104,11 +102,11 @@ void vqgen_cellmetric(vqgen *v){
 
   /* minimum, maximum, cell spacing */
   for(j=0;j<v->entries;j++){
-    double localmin=-1.;
+    float localmin=-1.;
 
     for(k=0;k<v->entries;k++){
       if(j!=k){
-       double this=_dist(v,_now(v,j),_now(v,k));
+       float this=_dist(v,_now(v,j),_now(v,k));
        if(this>0){
          if(v->assigned[k] && (localmin==-1 || this<localmin))
            localmin=this;
@@ -141,7 +139,7 @@ void vqgen_cellmetric(vqgen *v){
          min,mean/acc,max,unused,dup);
 
 #ifdef NOISY
-  qsort(spacings,count,sizeof(double),directdsort);
+  qsort(spacings,count,sizeof(float),directdsort);
   for(i=0;i<count;i++)
     fprintf(cells,"%g\n",spacings[i]);
   fclose(cells);
@@ -166,18 +164,18 @@ void vqgen_cellmetric(vqgen *v){
 
 void vqgen_quantize(vqgen *v,quant_meta *q){
 
-  double maxdel;
-  double mindel;
+  float maxdel;
+  float mindel;
 
-  double delta;
-  double maxquant=((1<<q->quant)-1);
+  float delta;
+  float maxquant=((1<<q->quant)-1);
 
   int j,k;
 
   mindel=maxdel=_now(v,0)[0];
   
   for(j=0;j<v->entries;j++){
-    double last=0.;
+    float last=0.f;
     for(k=0;k<v->elements;k++){
       if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
       if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
@@ -190,7 +188,7 @@ void vqgen_quantize(vqgen *v,quant_meta *q){
      encoded.  Loosen the delta slightly to allow for additional error
      during sequence quantization */
 
-  delta=(maxdel-mindel)/((1<<q->quant)-1.5);
+  delta=(maxdel-mindel)/((1<<q->quant)-1.5f);
 
   q->min=_float32_pack(mindel);
   q->delta=_float32_pack(delta);
@@ -199,10 +197,10 @@ void vqgen_quantize(vqgen *v,quant_meta *q){
   delta=_float32_unpack(q->delta);
 
   for(j=0;j<v->entries;j++){
-    double last=0;
+    float last=0;
     for(k=0;k<v->elements;k++){
-      double val=_now(v,j)[k];
-      double now=rint((val-last-mindel)/delta);
+      float val=_now(v,j)[k];
+      float now=rint((val-last-mindel)/delta);
       
       _now(v,j)[k]=now;
       if(now<0){
@@ -225,13 +223,13 @@ void vqgen_quantize(vqgen *v,quant_meta *q){
    scales; we just make sure they're properly offset. */
 void vqgen_unquantize(vqgen *v,quant_meta *q){
   long j,k;
-  double mindel=_float32_unpack(q->min);
-  double delta=_float32_unpack(q->delta);
+  float mindel=_float32_unpack(q->min);
+  float delta=_float32_unpack(q->delta);
 
   for(j=0;j<v->entries;j++){
-    double last=0.;
+    float last=0.f;
     for(k=0;k<v->elements;k++){
-      double now=_now(v,j)[k];
+      float now=_now(v,j)[k];
       now=fabs(now)*delta+last+mindel;
       if(q->sequencep)last=now;
       _now(v,j)[k]=now;
@@ -239,22 +237,23 @@ void vqgen_unquantize(vqgen *v,quant_meta *q){
   }
 }
 
-void vqgen_init(vqgen *v,int elements,int aux,int entries,double mindist,
-               double  (*metric)(vqgen *,double *, double *),
-               double *(*weight)(vqgen *,double *)){
+void vqgen_init(vqgen *v,int elements,int aux,int entries,float mindist,
+               float  (*metric)(vqgen *,float *, float *),
+               float *(*weight)(vqgen *,float *),int centroid){
   memset(v,0,sizeof(vqgen));
 
+  v->centroid=centroid;
   v->elements=elements;
   v->aux=aux;
   v->mindist=mindist;
   v->allocated=32768;
-  v->pointlist=malloc(v->allocated*(v->elements+v->aux)*sizeof(double));
+  v->pointlist=_ogg_malloc(v->allocated*(v->elements+v->aux)*sizeof(float));
 
   v->entries=entries;
-  v->entrylist=malloc(v->entries*v->elements*sizeof(double));
-  v->assigned=malloc(v->entries*sizeof(long));
-  v->bias=calloc(v->entries,sizeof(double));
-  v->max=calloc(v->entries,sizeof(double));
+  v->entrylist=_ogg_malloc(v->entries*v->elements*sizeof(float));
+  v->assigned=_ogg_malloc(v->entries*sizeof(long));
+  v->bias=_ogg_calloc(v->entries,sizeof(float));
+  v->max=_ogg_calloc(v->entries,sizeof(float));
   if(metric)
     v->metric_func=metric;
   else
@@ -264,35 +263,86 @@ void vqgen_init(vqgen *v,int elements,int aux,int entries,double mindist,
   else
     v->weight_func=_weight_null;
 
+  v->asciipoints=tmpfile();
+
 }
 
-void vqgen_addpoint(vqgen *v, double *p,double *a){
+void vqgen_addpoint(vqgen *v, float *p,float *a){
+  int k;
+  for(k=0;k<v->elements;k++)
+    fprintf(v->asciipoints,"%.12g\n",p[k]);
+  for(k=0;k<v->aux;k++)
+    fprintf(v->asciipoints,"%.12g\n",a[k]);
+
   if(v->points>=v->allocated){
     v->allocated*=2;
-    v->pointlist=realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
-                        sizeof(double));
+    v->pointlist=_ogg_realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
+                        sizeof(float));
+  }
+
+  memcpy(_point(v,v->points),p,sizeof(float)*v->elements);
+  if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(float)*v->aux);
+  /* quantize to the density mesh if it's selected */
+  if(v->mindist>0.f){
+    /* quantize to the mesh */
+    for(k=0;k<v->elements+v->aux;k++)
+      _point(v,v->points)[k]=
+       rint(_point(v,v->points)[k]/v->mindist)*v->mindist;
   }
-  
-  memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
-  if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(double)*v->aux);
   v->points++;
-  if(v->points==v->entries)_vqgen_seed(v);
   if(!(v->points&0xff))spinnit("loading... ",v->points);
 }
 
-double vqgen_iterate(vqgen *v,int biasp){
+/* yes, not threadsafe.  These utils aren't */
+static int sortit=0;
+static int sortsize=0;
+static int meshcomp(const void *a,const void *b){
+  if(((sortit++)&0xfff)==0)spinnit("sorting mesh...",sortit);
+  return(memcmp(a,b,sortsize));
+}
+
+void vqgen_sortmesh(vqgen *v){
+  sortit=0;
+  if(v->mindist>0.f){
+    long i,march=1;
+
+    /* sort to make uniqueness detection trivial */
+    sortsize=(v->elements+v->aux)*sizeof(float);
+    qsort(v->pointlist,v->points,sortsize,meshcomp);
+
+    /* now march through and eliminate dupes */
+    for(i=1;i<v->points;i++){
+      if(memcmp(_point(v,i),_point(v,i-1),sortsize)){
+       /* a new, unique entry.  march it down */
+       if(i>march)memcpy(_point(v,march),_point(v,i),sortsize);
+       march++;
+      }
+      spinnit("eliminating density... ",v->points-i);
+    }
+
+    /* we're done */
+    fprintf(stderr,"\r%ld training points remining out of %ld"
+           " after density mesh (%ld%%)\n",march,v->points,march*100/v->points);
+    v->points=march;
+
+  }
+  v->sorted=1;
+}
+
+float vqgen_iterate(vqgen *v,int biasp){
   long   i,j,k;
-  long   biasable;
 
-  double fdesired=(double)v->points/v->entries;
-  long  desired=fdesired;
-  long  desired2=desired*2;
+  float fdesired;
+  long  desired;
+  long  desired2;
 
-  double asserror=0.;
-  double meterror=0.;
-  double *new=malloc(sizeof(double)*v->entries*v->elements);
-  long   *nearcount=malloc(v->entries*sizeof(long));
-  double *nearbias=malloc(v->entries*desired2*sizeof(double));
+  float asserror=0.f;
+  float meterror=0.f;
+  float *new;
+  float *new2;
+  long   *nearcount;
+  float *nearbias;
  #ifdef NOISY
    char buff[80];
    FILE *assig;
@@ -312,24 +362,33 @@ double vqgen_iterate(vqgen *v,int biasp){
     exit(1);
   }
 
+  if(!v->sorted)vqgen_sortmesh(v);
+  if(!v->seeded)_vqgen_seed(v);
+
+  fdesired=(float)v->points/v->entries;
+  desired=fdesired;
+  desired2=desired*2;
+  new=_ogg_malloc(sizeof(float)*v->entries*v->elements);
+  new2=_ogg_malloc(sizeof(float)*v->entries*v->elements);
+  nearcount=_ogg_malloc(v->entries*sizeof(long));
+  nearbias=_ogg_malloc(v->entries*desired2*sizeof(float));
+
   /* fill in nearest points for entry biasing */
-  /*memset(v->bias,0,sizeof(double)*v->entries);*/
+  /*memset(v->bias,0,sizeof(float)*v->entries);*/
   memset(nearcount,0,sizeof(long)*v->entries);
   memset(v->assigned,0,sizeof(long)*v->entries);
-  biasable=0;
   if(biasp){
     for(i=0;i<v->points;i++){
-      double *ppt=v->weight_func(v,_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];
+      float *ppt=v->weight_func(v,_point(v,i));
+      float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
+      float secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
       long   firstentry=0;
       long   secondentry=1;
-      int    biasflag=1;
       
       if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
       
       if(firstmetric>secondmetric){
-       double temp=firstmetric;
+       float temp=firstmetric;
        firstmetric=secondmetric;
        secondmetric=temp;
        firstentry=1;
@@ -337,7 +396,7 @@ double vqgen_iterate(vqgen *v,int biasp){
       }
       
       for(j=2;j<v->entries;j++){
-       double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
+       float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
        if(thismetric<secondmetric){
          if(thismetric<firstmetric){
            secondmetric=firstmetric;
@@ -354,8 +413,8 @@ double vqgen_iterate(vqgen *v,int biasp){
       j=firstentry;
       for(j=0;j<v->entries;j++){
        
-       double thismetric,localmetric;
-       double *nearbiasptr=nearbias+desired2*j;
+       float thismetric,localmetric;
+       float *nearbiasptr=nearbias+desired2*j;
        long k=nearcount[j];
        
        localmetric=v->metric_func(v,_now(v,j),ppt);
@@ -373,69 +432,56 @@ double vqgen_iterate(vqgen *v,int biasp){
           cells in a codebook to be roughly some minimum size (as with
           the low resolution residue books) */
        
-       if(localmetric>=v->mindist){
+       /* a cute two-stage delayed sorting hack */
+       if(k<desired){
+         nearbiasptr[k]=thismetric;
+         k++;
+         if(k==desired){
+           spinnit("biasing... ",v->points+v->points+v->entries-i);
+           qsort(nearbiasptr,desired,sizeof(float),directdsort);
+         }
          
-         /* a cute two-stage delayed sorting hack */
-         if(k<desired){
-           nearbiasptr[k]=thismetric;
-           k++;
-           if(k==desired){
-             spinnit("biasing... ",v->points+v->points+v->entries-i);
-             qsort(nearbiasptr,desired,sizeof(double),directdsort);
-           }
-           
-         }else if(thismetric>nearbiasptr[desired-1]){
-           nearbiasptr[k]=thismetric;
-           k++;
-           if(k==desired2){
-             spinnit("biasing... ",v->points+v->points+v->entries-i);
-             qsort(nearbiasptr,desired2,sizeof(double),directdsort);
-             k=desired;
-           }
+       }else if(thismetric>nearbiasptr[desired-1]){
+         nearbiasptr[k]=thismetric;
+         k++;
+         if(k==desired2){
+           spinnit("biasing... ",v->points+v->points+v->entries-i);
+           qsort(nearbiasptr,desired2,sizeof(float),directdsort);
+           k=desired;
          }
-         nearcount[j]=k;
-       }else
-         biasflag=0;
+       }
+       nearcount[j]=k;
       }
-      biasable+=biasflag;
     }
     
     /* inflate/deflate */
     
     for(i=0;i<v->entries;i++){
-      double *nearbiasptr=nearbias+desired2*i;
+      float *nearbiasptr=nearbias+desired2*i;
       
       spinnit("biasing... ",v->points+v->entries-i);
       
       /* due to the delayed sorting, we likely need to finish it off....*/
       if(nearcount[i]>desired)
-       qsort(nearbiasptr,nearcount[i],sizeof(double),directdsort);
-      
-      /* desired is the *maximum* bias queue size.  If we're using
-        minimum distance, we're not interested in the max size... we're
-        interested in the biasable number of points */
-      {
-       long localdesired=(double)biasable/v->entries;
-       if(localdesired)
-         v->bias[i]=nearbiasptr[localdesired-1];
-       else
-         v->bias[i]=nearbiasptr[0];
-      }
+       qsort(nearbiasptr,nearcount[i],sizeof(float),directdsort);
+
+      v->bias[i]=nearbiasptr[desired-1];
+
     }
   }else{ 
-    memset(v->bias,0,v->entries*sizeof(double));
+    memset(v->bias,0,v->entries*sizeof(float));
   }
 
   /* Now assign with new bias and find new midpoints */
   for(i=0;i<v->points;i++){
-    double *ppt=v->weight_func(v,_point(v,i));
-    double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
+    float *ppt=v->weight_func(v,_point(v,i));
+    float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
     long   firstentry=0;
 
     if(!(i&0xff))spinnit("centering... ",v->points-i);
 
     for(j=0;j<v->entries;j++){
-      double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
+      float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
       if(thismetric<firstmetric){
        firstmetric=thismetric;
        firstentry=j;
@@ -450,17 +496,35 @@ double vqgen_iterate(vqgen *v,int biasp){
           ppt[0],ppt[1]);
 #endif
 
-    firstmetric-=v->bias[firstentry];
+    firstmetric-=v->bias[j];
     meterror+=firstmetric;
-    /* set up midpoints for next iter */
-    if(v->assigned[j]++){
-      for(k=0;k<v->elements;k++)
-       vN(new,j)[k]+=ppt[k];
-      if(firstmetric>v->max[firstentry])v->max[firstentry]=firstmetric;
+
+    if(v->centroid==0){
+      /* set up midpoints for next iter */
+      if(v->assigned[j]++){
+       for(k=0;k<v->elements;k++)
+         vN(new,j)[k]+=ppt[k];
+       if(firstmetric>v->max[j])v->max[j]=firstmetric;
+      }else{
+       for(k=0;k<v->elements;k++)
+         vN(new,j)[k]=ppt[k];
+       v->max[j]=firstmetric;
+      }
     }else{
-      for(k=0;k<v->elements;k++)
-       vN(new,j)[k]=ppt[k];
-      v->max[firstentry]=firstmetric;
+      /* centroid */
+      if(v->assigned[j]++){
+       for(k=0;k<v->elements;k++){
+         if(vN(new,j)[k]>ppt[k])vN(new,j)[k]=ppt[k];
+         if(vN(new2,j)[k]<ppt[k])vN(new2,j)[k]=ppt[k];
+       }
+       if(firstmetric>v->max[firstentry])v->max[j]=firstmetric;
+      }else{
+       for(k=0;k<v->elements;k++){
+         vN(new,j)[k]=ppt[k];
+         vN(new2,j)[k]=ppt[k];
+       }
+       v->max[firstentry]=firstmetric;
+      }
     }
   }
 
@@ -472,9 +536,15 @@ double vqgen_iterate(vqgen *v,int biasp){
     fprintf(bias,"%g\n",v->bias[j]);
 #endif
     asserror+=fabs(v->assigned[j]-fdesired);
-    if(v->assigned[j])
-      for(k=0;k<v->elements;k++)
-       _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
+    if(v->assigned[j]){
+      if(v->centroid==0){
+       for(k=0;k<v->elements;k++)
+         _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
+      }else{
+       for(k=0;k<v->elements;k++)
+         _now(v,j)[k]=(vN(new,j)[k]+vN(new2,j)[k])/2.f;
+      }
+    }
   }
 
   asserror/=(v->entries*fdesired);