]> git.jsancho.org Git - lugaru.git/blob - libvorbis-1.0.1/lib/floor1.c
Force XYZ to initialize to please Valgrind.
[lugaru.git] / libvorbis-1.0.1 / lib / floor1.c
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-2002             *
9  * by the XIPHOPHORUS Company http://www.xiph.org/                  *
10  *                                                                  *
11  ********************************************************************
12
13  function: floor backend 1 implementation
14  last mod: $Id: floor1.c,v 1.26 2003/02/15 07:10:07 xiphmont Exp $
15
16  ********************************************************************/
17
18 #include <stdlib.h>
19 #include <string.h>
20 #include <math.h>
21 #include <ogg/ogg.h>
22 #include "vorbis/codec.h"
23 #include "codec_internal.h"
24 #include "registry.h"
25 #include "codebook.h"
26 #include "misc.h"
27 #include "scales.h"
28
29 #include <stdio.h>
30
31 #define floor1_rangedB 140 /* floor 1 fixed at -140dB to 0dB range */
32
33 typedef struct {
34   int sorted_index[VIF_POSIT+2];
35   int forward_index[VIF_POSIT+2];
36   int reverse_index[VIF_POSIT+2];
37   
38   int hineighbor[VIF_POSIT];
39   int loneighbor[VIF_POSIT];
40   int posts;
41
42   int n;
43   int quant_q;
44   vorbis_info_floor1 *vi;
45
46   long phrasebits;
47   long postbits;
48   long frames;
49 } vorbis_look_floor1;
50
51 typedef struct lsfit_acc{
52   long x0;
53   long x1;
54
55   long xa;
56   long ya;
57   long x2a;
58   long y2a;
59   long xya; 
60   long an;
61 } lsfit_acc;
62
63 /***********************************************/
64  
65 static void floor1_free_info(vorbis_info_floor *i){
66   vorbis_info_floor1 *info=(vorbis_info_floor1 *)i;
67   if(info){
68     memset(info,0,sizeof(*info));
69     _ogg_free(info);
70   }
71 }
72
73 static void floor1_free_look(vorbis_look_floor *i){
74   vorbis_look_floor1 *look=(vorbis_look_floor1 *)i;
75   if(look){
76     /*fprintf(stderr,"floor 1 bit usage %f:%f (%f total)\n",
77             (float)look->phrasebits/look->frames,
78             (float)look->postbits/look->frames,
79             (float)(look->postbits+look->phrasebits)/look->frames);*/
80
81     memset(look,0,sizeof(*look));
82     _ogg_free(look);
83   }
84 }
85
86 static int ilog(unsigned int v){
87   int ret=0;
88   while(v){
89     ret++;
90     v>>=1;
91   }
92   return(ret);
93 }
94
95 static int ilog2(unsigned int v){
96   int ret=0;
97   if(v)--v;
98   while(v){
99     ret++;
100     v>>=1;
101   }
102   return(ret);
103 }
104
105 static void floor1_pack (vorbis_info_floor *i,oggpack_buffer *opb){
106   vorbis_info_floor1 *info=(vorbis_info_floor1 *)i;
107   int j,k;
108   int count=0;
109   int rangebits;
110   int maxposit=info->postlist[1];
111   int maxclass=-1;
112
113   /* save out partitions */
114   oggpack_write(opb,info->partitions,5); /* only 0 to 31 legal */
115   for(j=0;j<info->partitions;j++){
116     oggpack_write(opb,info->partitionclass[j],4); /* only 0 to 15 legal */
117     if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j];
118   }
119
120   /* save out partition classes */
121   for(j=0;j<maxclass+1;j++){
122     oggpack_write(opb,info->class_dim[j]-1,3); /* 1 to 8 */
123     oggpack_write(opb,info->class_subs[j],2); /* 0 to 3 */
124     if(info->class_subs[j])oggpack_write(opb,info->class_book[j],8);
125     for(k=0;k<(1<<info->class_subs[j]);k++)
126       oggpack_write(opb,info->class_subbook[j][k]+1,8);
127   }
128
129   /* save out the post list */
130   oggpack_write(opb,info->mult-1,2);     /* only 1,2,3,4 legal now */ 
131   oggpack_write(opb,ilog2(maxposit),4);
132   rangebits=ilog2(maxposit);
133
134   for(j=0,k=0;j<info->partitions;j++){
135     count+=info->class_dim[info->partitionclass[j]]; 
136     for(;k<count;k++)
137       oggpack_write(opb,info->postlist[k+2],rangebits);
138   }
139 }
140
141
142 static vorbis_info_floor *floor1_unpack (vorbis_info *vi,oggpack_buffer *opb){
143   codec_setup_info     *ci=vi->codec_setup;
144   int j,k,count=0,maxclass=-1,rangebits;
145
146   vorbis_info_floor1 *info=_ogg_calloc(1,sizeof(*info));
147   /* read partitions */
148   info->partitions=oggpack_read(opb,5); /* only 0 to 31 legal */
149   for(j=0;j<info->partitions;j++){
150     info->partitionclass[j]=oggpack_read(opb,4); /* only 0 to 15 legal */
151     if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j];
152   }
153
154   /* read partition classes */
155   for(j=0;j<maxclass+1;j++){
156     info->class_dim[j]=oggpack_read(opb,3)+1; /* 1 to 8 */
157     info->class_subs[j]=oggpack_read(opb,2); /* 0,1,2,3 bits */
158     if(info->class_subs[j]<0)
159       goto err_out;
160     if(info->class_subs[j])info->class_book[j]=oggpack_read(opb,8);
161     if(info->class_book[j]<0 || info->class_book[j]>=ci->books)
162       goto err_out;
163     for(k=0;k<(1<<info->class_subs[j]);k++){
164       info->class_subbook[j][k]=oggpack_read(opb,8)-1;
165       if(info->class_subbook[j][k]<-1 || info->class_subbook[j][k]>=ci->books)
166         goto err_out;
167     }
168   }
169
170   /* read the post list */
171   info->mult=oggpack_read(opb,2)+1;     /* only 1,2,3,4 legal now */ 
172   rangebits=oggpack_read(opb,4);
173
174   for(j=0,k=0;j<info->partitions;j++){
175     count+=info->class_dim[info->partitionclass[j]]; 
176     for(;k<count;k++){
177       int t=info->postlist[k+2]=oggpack_read(opb,rangebits);
178       if(t<0 || t>=(1<<rangebits))
179         goto err_out;
180     }
181   }
182   info->postlist[0]=0;
183   info->postlist[1]=1<<rangebits;
184
185   return(info);
186   
187  err_out:
188   floor1_free_info(info);
189   return(NULL);
190 }
191
192 static int icomp(const void *a,const void *b){
193   return(**(int **)a-**(int **)b);
194 }
195
196 static vorbis_look_floor *floor1_look(vorbis_dsp_state *vd,
197                                       vorbis_info_floor *in){
198
199   int *sortpointer[VIF_POSIT+2];
200   vorbis_info_floor1 *info=(vorbis_info_floor1 *)in;
201   vorbis_look_floor1 *look=_ogg_calloc(1,sizeof(*look));
202   int i,j,n=0;
203
204   look->vi=info;
205   look->n=info->postlist[1];
206  
207   /* we drop each position value in-between already decoded values,
208      and use linear interpolation to predict each new value past the
209      edges.  The positions are read in the order of the position
210      list... we precompute the bounding positions in the lookup.  Of
211      course, the neighbors can change (if a position is declined), but
212      this is an initial mapping */
213
214   for(i=0;i<info->partitions;i++)n+=info->class_dim[info->partitionclass[i]];
215   n+=2;
216   look->posts=n;
217
218   /* also store a sorted position index */
219   for(i=0;i<n;i++)sortpointer[i]=info->postlist+i;
220   qsort(sortpointer,n,sizeof(*sortpointer),icomp);
221
222   /* points from sort order back to range number */
223   for(i=0;i<n;i++)look->forward_index[i]=sortpointer[i]-info->postlist;
224   /* points from range order to sorted position */
225   for(i=0;i<n;i++)look->reverse_index[look->forward_index[i]]=i;
226   /* we actually need the post values too */
227   for(i=0;i<n;i++)look->sorted_index[i]=info->postlist[look->forward_index[i]];
228   
229   /* quantize values to multiplier spec */
230   switch(info->mult){
231   case 1: /* 1024 -> 256 */
232     look->quant_q=256;
233     break;
234   case 2: /* 1024 -> 128 */
235     look->quant_q=128;
236     break;
237   case 3: /* 1024 -> 86 */
238     look->quant_q=86;
239     break;
240   case 4: /* 1024 -> 64 */
241     look->quant_q=64;
242     break;
243   }
244
245   /* discover our neighbors for decode where we don't use fit flags
246      (that would push the neighbors outward) */
247   for(i=0;i<n-2;i++){
248     int lo=0;
249     int hi=1;
250     int lx=0;
251     int hx=look->n;
252     int currentx=info->postlist[i+2];
253     for(j=0;j<i+2;j++){
254       int x=info->postlist[j];
255       if(x>lx && x<currentx){
256         lo=j;
257         lx=x;
258       }
259       if(x<hx && x>currentx){
260         hi=j;
261         hx=x;
262       }
263     }
264     look->loneighbor[i]=lo;
265     look->hineighbor[i]=hi;
266   }
267
268   return(look);
269 }
270
271 static int render_point(int x0,int x1,int y0,int y1,int x){
272   y0&=0x7fff; /* mask off flag */
273   y1&=0x7fff;
274     
275   {
276     int dy=y1-y0;
277     int adx=x1-x0;
278     int ady=abs(dy);
279     int err=ady*(x-x0);
280     
281     int off=err/adx;
282     if(dy<0)return(y0-off);
283     return(y0+off);
284   }
285 }
286
287 static int vorbis_dBquant(const float *x){
288   int i= *x*7.3142857f+1023.5f;
289   if(i>1023)return(1023);
290   if(i<0)return(0);
291   return i;
292 }
293
294 static float FLOOR1_fromdB_LOOKUP[256]={
295   1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F, 
296   1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F, 
297   1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F, 
298   2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F, 
299   2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F, 
300   3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F, 
301   4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F, 
302   6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F, 
303   7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F, 
304   1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F, 
305   1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F, 
306   1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F, 
307   2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F, 
308   2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F, 
309   3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F, 
310   4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F, 
311   5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F, 
312   7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F, 
313   9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F, 
314   1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F, 
315   1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F, 
316   2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F, 
317   2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F, 
318   3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F, 
319   4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F, 
320   5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F, 
321   7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F, 
322   9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F, 
323   0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F, 
324   0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F, 
325   0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F, 
326   0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F, 
327   0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F, 
328   0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F, 
329   0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F, 
330   0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F, 
331   0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F, 
332   0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F, 
333   0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F, 
334   0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F, 
335   0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F, 
336   0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F, 
337   0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F, 
338   0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F, 
339   0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F, 
340   0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F, 
341   0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F, 
342   0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F, 
343   0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F, 
344   0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F, 
345   0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F, 
346   0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F, 
347   0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F, 
348   0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F, 
349   0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F, 
350   0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F, 
351   0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F, 
352   0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F, 
353   0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F, 
354   0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F, 
355   0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F, 
356   0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F, 
357   0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F, 
358   0.82788260F, 0.88168307F, 0.9389798F, 1.F, 
359 };
360
361 static void render_line(int x0,int x1,int y0,int y1,float *d){
362   int dy=y1-y0;
363   int adx=x1-x0;
364   int ady=abs(dy);
365   int base=dy/adx;
366   int sy=(dy<0?base-1:base+1);
367   int x=x0;
368   int y=y0;
369   int err=0;
370
371   ady-=abs(base*adx);
372
373   d[x]*=FLOOR1_fromdB_LOOKUP[y];
374   while(++x<x1){
375     err=err+ady;
376     if(err>=adx){
377       err-=adx;
378       y+=sy;
379     }else{
380       y+=base;
381     }
382     d[x]*=FLOOR1_fromdB_LOOKUP[y];
383   }
384 }
385
386 static void render_line0(int x0,int x1,int y0,int y1,int *d){
387   int dy=y1-y0;
388   int adx=x1-x0;
389   int ady=abs(dy);
390   int base=dy/adx;
391   int sy=(dy<0?base-1:base+1);
392   int x=x0;
393   int y=y0;
394   int err=0;
395
396   ady-=abs(base*adx);
397
398   d[x]=y;
399   while(++x<x1){
400     err=err+ady;
401     if(err>=adx){
402       err-=adx;
403       y+=sy;
404     }else{
405       y+=base;
406     }
407     d[x]=y;
408   }
409 }
410
411 /* the floor has already been filtered to only include relevant sections */
412 static int accumulate_fit(const float *flr,const float *mdct,
413                           int x0, int x1,lsfit_acc *a,
414                           int n,vorbis_info_floor1 *info){
415   long i;
416   int quantized=vorbis_dBquant(flr+x0);
417
418   long xa=0,ya=0,x2a=0,y2a=0,xya=0,na=0, xb=0,yb=0,x2b=0,y2b=0,xyb=0,nb=0;
419
420   memset(a,0,sizeof(*a));
421   a->x0=x0;
422   a->x1=x1;
423   if(x1>=n)x1=n-1;
424
425   for(i=x0;i<=x1;i++){
426     int quantized=vorbis_dBquant(flr+i);
427     if(quantized){
428       if(mdct[i]+info->twofitatten>=flr[i]){
429         xa  += i;
430         ya  += quantized;
431         x2a += i*i;
432         y2a += quantized*quantized;
433         xya += i*quantized;
434         na++;
435       }else{
436         xb  += i;
437         yb  += quantized;
438         x2b += i*i;
439         y2b += quantized*quantized;
440         xyb += i*quantized;
441         nb++;
442       }
443     }
444   }
445
446   xb+=xa;
447   yb+=ya;
448   x2b+=x2a;
449   y2b+=y2a;
450   xyb+=xya;
451   nb+=na;
452
453   /* weight toward the actually used frequencies if we meet the threshhold */
454   {
455     int weight=nb*info->twofitweight/(na+1);
456
457     a->xa=xa*weight+xb;
458     a->ya=ya*weight+yb;
459     a->x2a=x2a*weight+x2b;
460     a->y2a=y2a*weight+y2b;
461     a->xya=xya*weight+xyb;
462     a->an=na*weight+nb;
463   }
464
465   return(na);
466 }
467
468 static void fit_line(lsfit_acc *a,int fits,int *y0,int *y1){
469   long x=0,y=0,x2=0,y2=0,xy=0,an=0,i;
470   long x0=a[0].x0;
471   long x1=a[fits-1].x1;
472
473   for(i=0;i<fits;i++){
474     x+=a[i].xa;
475     y+=a[i].ya;
476     x2+=a[i].x2a;
477     y2+=a[i].y2a;
478     xy+=a[i].xya;
479     an+=a[i].an;
480   }
481
482   if(*y0>=0){
483     x+=   x0;
484     y+=  *y0;
485     x2+=  x0 *  x0;
486     y2+= *y0 * *y0;
487     xy+= *y0 *  x0;
488     an++;
489   }
490
491   if(*y1>=0){
492     x+=   x1;
493     y+=  *y1;
494     x2+=  x1 *  x1;
495     y2+= *y1 * *y1;
496     xy+= *y1 *  x1;
497     an++;
498   }
499   
500   if(an){
501     /* need 64 bit multiplies, which C doesn't give portably as int */
502     double fx=x;
503     double fy=y;
504     double fx2=x2;
505     double fxy=xy;
506     double denom=1./(an*fx2-fx*fx);
507     double a=(fy*fx2-fxy*fx)*denom;
508     double b=(an*fxy-fx*fy)*denom;
509     *y0=rint(a+b*x0);
510     *y1=rint(a+b*x1);
511     
512     /* limit to our range! */
513     if(*y0>1023)*y0=1023;
514     if(*y1>1023)*y1=1023;
515     if(*y0<0)*y0=0;
516     if(*y1<0)*y1=0;
517     
518   }else{
519     *y0=0;
520     *y1=0;
521   }
522 }
523
524 /*static void fit_line_point(lsfit_acc *a,int fits,int *y0,int *y1){
525   long y=0;
526   int i;
527
528   for(i=0;i<fits && y==0;i++)
529     y+=a[i].ya;
530   
531   *y0=*y1=y;
532   }*/
533
534 static int inspect_error(int x0,int x1,int y0,int y1,const float *mask,
535                          const float *mdct,
536                          vorbis_info_floor1 *info){
537   int dy=y1-y0;
538   int adx=x1-x0;
539   int ady=abs(dy);
540   int base=dy/adx;
541   int sy=(dy<0?base-1:base+1);
542   int x=x0;
543   int y=y0;
544   int err=0;
545   int val=vorbis_dBquant(mask+x);
546   int mse=0;
547   int n=0;
548
549   ady-=abs(base*adx);
550   
551   mse=(y-val);
552   mse*=mse;
553   n++;
554   if(mdct[x]+info->twofitatten>=mask[x]){
555     if(y+info->maxover<val)return(1);
556     if(y-info->maxunder>val)return(1);
557   }
558
559   while(++x<x1){
560     err=err+ady;
561     if(err>=adx){
562       err-=adx;
563       y+=sy;
564     }else{
565       y+=base;
566     }
567
568     val=vorbis_dBquant(mask+x);
569     mse+=((y-val)*(y-val));
570     n++;
571     if(mdct[x]+info->twofitatten>=mask[x]){
572       if(val){
573         if(y+info->maxover<val)return(1);
574         if(y-info->maxunder>val)return(1);
575       }
576     }
577   }
578   
579   if(info->maxover*info->maxover/n>info->maxerr)return(0);
580   if(info->maxunder*info->maxunder/n>info->maxerr)return(0);
581   if(mse/n>info->maxerr)return(1);
582   return(0);
583 }
584
585 static int post_Y(int *A,int *B,int pos){
586   if(A[pos]<0)
587     return B[pos];
588   if(B[pos]<0)
589     return A[pos];
590
591   return (A[pos]+B[pos])>>1;
592 }
593
594 static int seq=0;
595
596 int *floor1_fit(vorbis_block *vb,vorbis_look_floor1 *look,
597                           const float *logmdct,   /* in */
598                           const float *logmask){
599   long i,j;
600   vorbis_info_floor1 *info=look->vi;
601   long n=look->n;
602   long posts=look->posts;
603   long nonzero=0;
604   lsfit_acc fits[VIF_POSIT+1];
605   int fit_valueA[VIF_POSIT+2]; /* index by range list position */
606   int fit_valueB[VIF_POSIT+2]; /* index by range list position */
607
608   int loneighbor[VIF_POSIT+2]; /* sorted index of range list position (+2) */
609   int hineighbor[VIF_POSIT+2]; 
610   int *output=NULL;
611   int memo[VIF_POSIT+2];
612
613   for(i=0;i<posts;i++)fit_valueA[i]=-200; /* mark all unused */
614   for(i=0;i<posts;i++)fit_valueB[i]=-200; /* mark all unused */
615   for(i=0;i<posts;i++)loneighbor[i]=0; /* 0 for the implicit 0 post */
616   for(i=0;i<posts;i++)hineighbor[i]=1; /* 1 for the implicit post at n */
617   for(i=0;i<posts;i++)memo[i]=-1;      /* no neighbor yet */
618
619   /* quantize the relevant floor points and collect them into line fit
620      structures (one per minimal division) at the same time */
621   if(posts==0){
622     nonzero+=accumulate_fit(logmask,logmdct,0,n,fits,n,info);
623   }else{
624     for(i=0;i<posts-1;i++)
625       nonzero+=accumulate_fit(logmask,logmdct,look->sorted_index[i],
626                               look->sorted_index[i+1],fits+i,
627                               n,info);
628   }
629   
630   if(nonzero){
631     /* start by fitting the implicit base case.... */
632     int y0=-200;
633     int y1=-200;
634     fit_line(fits,posts-1,&y0,&y1);
635
636     fit_valueA[0]=y0;
637     fit_valueB[0]=y0;
638     fit_valueB[1]=y1;
639     fit_valueA[1]=y1;
640
641     /* Non degenerate case */
642     /* start progressive splitting.  This is a greedy, non-optimal
643        algorithm, but simple and close enough to the best
644        answer. */
645     for(i=2;i<posts;i++){
646       int sortpos=look->reverse_index[i];
647       int ln=loneighbor[sortpos];
648       int hn=hineighbor[sortpos];
649       
650       /* eliminate repeat searches of a particular range with a memo */
651       if(memo[ln]!=hn){
652         /* haven't performed this error search yet */
653         int lsortpos=look->reverse_index[ln];
654         int hsortpos=look->reverse_index[hn];
655         memo[ln]=hn;
656                 
657         {
658           /* A note: we want to bound/minimize *local*, not global, error */
659           int lx=info->postlist[ln];
660           int hx=info->postlist[hn];      
661           int ly=post_Y(fit_valueA,fit_valueB,ln);
662           int hy=post_Y(fit_valueA,fit_valueB,hn);
663           
664           if(ly==-1 || hy==-1){
665             exit(1);
666           }
667
668           if(inspect_error(lx,hx,ly,hy,logmask,logmdct,info)){
669             /* outside error bounds/begin search area.  Split it. */
670             int ly0=-200;
671             int ly1=-200;
672             int hy0=-200;
673             int hy1=-200;
674             fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1);
675             fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1);
676             
677             /* store new edge values */
678             fit_valueB[ln]=ly0;
679             if(ln==0)fit_valueA[ln]=ly0;
680             fit_valueA[i]=ly1;
681             fit_valueB[i]=hy0;
682             fit_valueA[hn]=hy1;
683             if(hn==1)fit_valueB[hn]=hy1;
684             
685             if(ly1>=0 || hy0>=0){
686               /* store new neighbor values */
687               for(j=sortpos-1;j>=0;j--)
688                 if(hineighbor[j]==hn)
689                   hineighbor[j]=i;
690                 else
691                   break;
692               for(j=sortpos+1;j<posts;j++)
693                 if(loneighbor[j]==ln)
694                   loneighbor[j]=i;
695                 else
696                   break;
697               
698             }
699           }else{
700             
701             fit_valueA[i]=-200;
702             fit_valueB[i]=-200;
703           }
704         }
705       }
706     }
707   
708     output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
709     
710     output[0]=post_Y(fit_valueA,fit_valueB,0);
711     output[1]=post_Y(fit_valueA,fit_valueB,1);
712     
713     /* fill in posts marked as not using a fit; we will zero
714        back out to 'unused' when encoding them so long as curve
715        interpolation doesn't force them into use */
716     for(i=2;i<posts;i++){
717       int ln=look->loneighbor[i-2];
718       int hn=look->hineighbor[i-2];
719       int x0=info->postlist[ln];
720       int x1=info->postlist[hn];
721       int y0=output[ln];
722       int y1=output[hn];
723       
724       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
725       int vx=post_Y(fit_valueA,fit_valueB,i);
726       
727       if(vx>=0 && predicted!=vx){ 
728         output[i]=vx;
729       }else{
730         output[i]= predicted|0x8000;
731       }
732     }
733   }
734
735   return(output);
736   
737 }
738                 
739 int *floor1_interpolate_fit(vorbis_block *vb,vorbis_look_floor1 *look,
740                           int *A,int *B,
741                           int del){
742
743   long i;
744   long posts=look->posts;
745   int *output=NULL;
746   
747   if(A && B){
748     output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
749     
750     for(i=0;i<posts;i++){
751       output[i]=((65536-del)*(A[i]&0x7fff)+del*(B[i]&0x7fff)+32768)>>16;
752       if(A[i]&0x8000 && B[i]&0x8000)output[i]|=0x8000;
753     }
754   }
755
756   return(output);
757 }
758
759
760 int floor1_encode(vorbis_block *vb,vorbis_look_floor1 *look,
761                   int *post,int *ilogmask){
762
763   long i,j;
764   vorbis_info_floor1 *info=look->vi;
765   long n=look->n;
766   long posts=look->posts;
767   codec_setup_info *ci=vb->vd->vi->codec_setup;
768   int out[VIF_POSIT+2];
769   static_codebook **sbooks=ci->book_param;
770   codebook *books=ci->fullbooks;
771   static long seq=0; 
772
773   /* quantize values to multiplier spec */
774   if(post){
775     for(i=0;i<posts;i++){
776       int val=post[i]&0x7fff;
777       switch(info->mult){
778       case 1: /* 1024 -> 256 */
779         val>>=2;
780         break;
781       case 2: /* 1024 -> 128 */
782         val>>=3;
783         break;
784       case 3: /* 1024 -> 86 */
785         val/=12;
786         break;
787       case 4: /* 1024 -> 64 */
788         val>>=4;
789         break;
790       }
791       post[i]=val | (post[i]&0x8000);
792     }
793
794     out[0]=post[0];
795     out[1]=post[1];
796
797     /* find prediction values for each post and subtract them */
798     for(i=2;i<posts;i++){
799       int ln=look->loneighbor[i-2];
800       int hn=look->hineighbor[i-2];
801       int x0=info->postlist[ln];
802       int x1=info->postlist[hn];
803       int y0=post[ln];
804       int y1=post[hn];
805       
806       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
807       
808       if((post[i]&0x8000) || (predicted==post[i])){
809         post[i]=predicted|0x8000; /* in case there was roundoff jitter
810                                      in interpolation */
811         out[i]=0;
812       }else{
813         int headroom=(look->quant_q-predicted<predicted?
814                       look->quant_q-predicted:predicted);
815         
816         int val=post[i]-predicted;
817         
818         /* at this point the 'deviation' value is in the range +/- max
819            range, but the real, unique range can always be mapped to
820            only [0-maxrange).  So we want to wrap the deviation into
821            this limited range, but do it in the way that least screws
822            an essentially gaussian probability distribution. */
823         
824         if(val<0)
825           if(val<-headroom)
826             val=headroom-val-1;
827           else
828             val=-1-(val<<1);
829         else
830           if(val>=headroom)
831             val= val+headroom;
832           else
833             val<<=1;
834         
835         out[i]=val;
836         post[ln]&=0x7fff;
837         post[hn]&=0x7fff;
838       }
839     }
840     
841     /* we have everything we need. pack it out */
842     /* mark nontrivial floor */
843     oggpack_write(&vb->opb,1,1);
844       
845     /* beginning/end post */
846     look->frames++;
847     look->postbits+=ilog(look->quant_q-1)*2;
848     oggpack_write(&vb->opb,out[0],ilog(look->quant_q-1));
849     oggpack_write(&vb->opb,out[1],ilog(look->quant_q-1));
850       
851       
852     /* partition by partition */
853     for(i=0,j=2;i<info->partitions;i++){
854       int class=info->partitionclass[i];
855       int cdim=info->class_dim[class];
856       int csubbits=info->class_subs[class];
857       int csub=1<<csubbits;
858       int bookas[8]={0,0,0,0,0,0,0,0};
859       int cval=0;
860       int cshift=0;
861       int k,l;
862
863       /* generate the partition's first stage cascade value */
864       if(csubbits){
865         int maxval[8];
866         for(k=0;k<csub;k++){
867           int booknum=info->class_subbook[class][k];
868           if(booknum<0){
869             maxval[k]=1;
870           }else{
871             maxval[k]=sbooks[info->class_subbook[class][k]]->entries;
872           }
873         }
874         for(k=0;k<cdim;k++){
875           for(l=0;l<csub;l++){
876             int val=out[j+k];
877             if(val<maxval[l]){
878               bookas[k]=l;
879               break;
880             }
881           }
882           cval|= bookas[k]<<cshift;
883           cshift+=csubbits;
884         }
885         /* write it */
886         look->phrasebits+=
887           vorbis_book_encode(books+info->class_book[class],cval,&vb->opb);
888         
889 #ifdef TRAIN_FLOOR1
890         {
891           FILE *of;
892           char buffer[80];
893           sprintf(buffer,"line_%dx%ld_class%d.vqd",
894                   vb->pcmend/2,posts-2,class);
895           of=fopen(buffer,"a");
896           fprintf(of,"%d\n",cval);
897           fclose(of);
898         }
899 #endif
900       }
901         
902       /* write post values */
903       for(k=0;k<cdim;k++){
904         int book=info->class_subbook[class][bookas[k]];
905         if(book>=0){
906           /* hack to allow training with 'bad' books */
907           if(out[j+k]<(books+book)->entries)
908             look->postbits+=vorbis_book_encode(books+book,
909                                                out[j+k],&vb->opb);
910           /*else
911             fprintf(stderr,"+!");*/
912           
913 #ifdef TRAIN_FLOOR1
914           {
915             FILE *of;
916             char buffer[80];
917             sprintf(buffer,"line_%dx%ld_%dsub%d.vqd",
918                     vb->pcmend/2,posts-2,class,bookas[k]);
919             of=fopen(buffer,"a");
920             fprintf(of,"%d\n",out[j+k]);
921             fclose(of);
922           }
923 #endif
924         }
925       }
926       j+=cdim;
927     }
928     
929     {
930       /* generate quantized floor equivalent to what we'd unpack in decode */
931       /* render the lines */
932       int hx=0;
933       int lx=0;
934       int ly=post[0]*info->mult;
935       for(j=1;j<look->posts;j++){
936         int current=look->forward_index[j];
937         int hy=post[current]&0x7fff;
938         if(hy==post[current]){
939           
940           hy*=info->mult;
941           hx=info->postlist[current];
942         
943           render_line0(lx,hx,ly,hy,ilogmask);
944         
945           lx=hx;
946           ly=hy;
947         }
948       }
949       for(j=hx;j<vb->pcmend/2;j++)ilogmask[j]=ly; /* be certain */    
950       seq++;
951       return(1);
952     }
953   }else{
954     oggpack_write(&vb->opb,0,1);
955     memset(ilogmask,0,vb->pcmend/2*sizeof(*ilogmask));
956     seq++;
957     return(0);
958   }
959 }
960
961 static void *floor1_inverse1(vorbis_block *vb,vorbis_look_floor *in){
962   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
963   vorbis_info_floor1 *info=look->vi;
964   codec_setup_info   *ci=vb->vd->vi->codec_setup;
965   
966   int i,j,k;
967   codebook *books=ci->fullbooks;   
968
969   /* unpack wrapped/predicted values from stream */
970   if(oggpack_read(&vb->opb,1)==1){
971     int *fit_value=_vorbis_block_alloc(vb,(look->posts)*sizeof(*fit_value));
972
973     fit_value[0]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
974     fit_value[1]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
975
976     /* partition by partition */
977     for(i=0,j=2;i<info->partitions;i++){
978       int class=info->partitionclass[i];
979       int cdim=info->class_dim[class];
980       int csubbits=info->class_subs[class];
981       int csub=1<<csubbits;
982       int cval=0;
983
984       /* decode the partition's first stage cascade value */
985       if(csubbits){
986         cval=vorbis_book_decode(books+info->class_book[class],&vb->opb);
987
988         if(cval==-1)goto eop;
989       }
990
991       for(k=0;k<cdim;k++){
992         int book=info->class_subbook[class][cval&(csub-1)];
993         cval>>=csubbits;
994         if(book>=0){
995           if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1)
996             goto eop;
997         }else{
998           fit_value[j+k]=0;
999         }
1000       }
1001       j+=cdim;
1002     }
1003
1004     /* unwrap positive values and reconsitute via linear interpolation */
1005     for(i=2;i<look->posts;i++){
1006       int predicted=render_point(info->postlist[look->loneighbor[i-2]],
1007                                  info->postlist[look->hineighbor[i-2]],
1008                                  fit_value[look->loneighbor[i-2]],
1009                                  fit_value[look->hineighbor[i-2]],
1010                                  info->postlist[i]);
1011       int hiroom=look->quant_q-predicted;
1012       int loroom=predicted;
1013       int room=(hiroom<loroom?hiroom:loroom)<<1;
1014       int val=fit_value[i];
1015
1016       if(val){
1017         if(val>=room){
1018           if(hiroom>loroom){
1019             val = val-loroom;
1020           }else{
1021             val = -1-(val-hiroom);
1022           }
1023         }else{
1024           if(val&1){
1025             val= -((val+1)>>1);
1026           }else{
1027             val>>=1;
1028           }
1029         }
1030
1031         fit_value[i]=val+predicted;
1032         fit_value[look->loneighbor[i-2]]&=0x7fff;
1033         fit_value[look->hineighbor[i-2]]&=0x7fff;
1034
1035       }else{
1036         fit_value[i]=predicted|0x8000;
1037       }
1038         
1039     }
1040
1041     return(fit_value);
1042   }
1043  eop:
1044   return(NULL);
1045 }
1046
1047 static int floor1_inverse2(vorbis_block *vb,vorbis_look_floor *in,void *memo,
1048                           float *out){
1049   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1050   vorbis_info_floor1 *info=look->vi;
1051
1052   codec_setup_info   *ci=vb->vd->vi->codec_setup;
1053   int                  n=ci->blocksizes[vb->W]/2;
1054   int j;
1055
1056   if(memo){
1057     /* render the lines */
1058     int *fit_value=(int *)memo;
1059     int hx=0;
1060     int lx=0;
1061     int ly=fit_value[0]*info->mult;
1062     for(j=1;j<look->posts;j++){
1063       int current=look->forward_index[j];
1064       int hy=fit_value[current]&0x7fff;
1065       if(hy==fit_value[current]){
1066         
1067         hy*=info->mult;
1068         hx=info->postlist[current];
1069         
1070         render_line(lx,hx,ly,hy,out);
1071         
1072         lx=hx;
1073         ly=hy;
1074       }
1075     }
1076     for(j=hx;j<n;j++)out[j]*=FLOOR1_fromdB_LOOKUP[ly]; /* be certain */    
1077     return(1);
1078   }
1079   memset(out,0,sizeof(*out)*n);
1080   return(0);
1081 }
1082
1083 /* export hooks */
1084 vorbis_func_floor floor1_exportbundle={
1085   &floor1_pack,&floor1_unpack,&floor1_look,&floor1_free_info,
1086   &floor1_free_look,&floor1_inverse1,&floor1_inverse2
1087 };
1088
1089