1 
2 #ifdef HAVE_CONFIG_H
3 #include "config.h"
4 #endif
5 #include <schroedinger/schro.h>
6 #include <schroedinger/schrohistogram.h>
7 #include <string.h>
8 #include <math.h>
9 
10 //#define DUMP_SUBBAND_CURVES
11 
12 void schro_encoder_choose_quantisers_simple (SchroEncoderFrame * frame);
13 void schro_encoder_choose_quantisers_rdo_bit_allocation (SchroEncoderFrame *
14     frame);
15 void schro_encoder_choose_quantisers_lossless (SchroEncoderFrame * frame);
16 void schro_encoder_choose_quantisers_lowdelay (SchroEncoderFrame * frame);
17 void schro_encoder_choose_quantisers_rdo_lambda (SchroEncoderFrame * frame);
18 void schro_encoder_choose_quantisers_rdo_cbr (SchroEncoderFrame * frame);
19 void schro_encoder_choose_quantisers_constant_error (SchroEncoderFrame * frame);
20 
21 double schro_encoder_entropy_to_lambda (SchroEncoderFrame * frame,
22     double entropy);
23 static double schro_encoder_lambda_to_entropy (SchroEncoderFrame * frame,
24     double lambda);
25 static void schro_encoder_generate_subband_histograms (SchroEncoderFrame *
26     frame);
27 
28 
29 static int schro_subband_pick_quant (SchroEncoderFrame * frame,
30     int component, int i, double lambda);
31 
32 #define CURVE_SIZE 128
33 
34 double
schro_encoder_perceptual_weight_moo(double cpd)35 schro_encoder_perceptual_weight_moo (double cpd)
36 {
37   /* I pretty much pulled this out of my ass (ppd = pixels per degree) */
38   if (cpd < 4)
39     return 1;
40   return 0.68 * cpd * exp (-0.25 * cpd);
41 }
42 
43 double
schro_encoder_perceptual_weight_constant(double cpd)44 schro_encoder_perceptual_weight_constant (double cpd)
45 {
46   return 1;
47 }
48 
49 double
schro_encoder_perceptual_weight_ccir959(double cpd)50 schro_encoder_perceptual_weight_ccir959 (double cpd)
51 {
52   double w;
53   w = 0.255 * pow (1 + 0.2561 * cpd * cpd, -0.75);
54 
55   /* return normalized value */
56   return w / 0.255;
57 }
58 
59 double
schro_encoder_perceptual_weight_manos_sakrison(double cpd)60 schro_encoder_perceptual_weight_manos_sakrison (double cpd)
61 {
62   double w;
63 
64   if (cpd < 4)
65     return 1;
66   w = 2.6 * (0.0192 + 0.114 * cpd) * exp (-pow (0.114 * cpd, 1.1));
67 
68   /* return normalized value */
69   return w / 0.980779694777866;
70 }
71 
72 static double
weighted_sum(const float * h1,const float * v1,double * weight)73 weighted_sum (const float *h1, const float *v1, double *weight)
74 {
75   int i, j;
76   double sum;
77   double rowsum;
78 
79   sum = 0;
80   for (j = 0; j < CURVE_SIZE; j++) {
81     rowsum = 0;
82     for (i = 0; i < CURVE_SIZE; i++) {
83       rowsum += h1[i] * v1[j] * weight[CURVE_SIZE * j + i];
84     }
85     sum += rowsum;
86   }
87   return sum;
88 }
89 
90 static double
dot_product(const float * h1,const float * v1,const float * h2,const float * v2,double * weight)91 dot_product (const float *h1, const float *v1, const float *h2,
92     const float *v2, double *weight)
93 {
94   int i, j;
95   double sum;
96   double rowsum;
97 
98   sum = 0;
99   for (j = 0; j < CURVE_SIZE; j++) {
100     rowsum = 0;
101     for (i = 0; i < CURVE_SIZE; i++) {
102       rowsum +=
103           h1[i] * v1[j] * h2[i] * v2[j] * weight[CURVE_SIZE * j +
104           i] * weight[CURVE_SIZE * j + i];
105     }
106     sum += rowsum;
107   }
108   return sum;
109 }
110 
111 static void
solve(double * matrix,double * column,int n)112 solve (double *matrix, double *column, int n)
113 {
114   int i;
115   int j;
116   int k;
117   double x;
118 
119   for (i = 0; i < n; i++) {
120     x = 1 / matrix[i * n + i];
121     for (k = i; k < n; k++) {
122       matrix[i * n + k] *= x;
123     }
124     column[i] *= x;
125 
126     for (j = i + 1; j < n; j++) {
127       x = matrix[j * n + i];
128       for (k = i; k < n; k++) {
129         matrix[j * n + k] -= matrix[i * n + k] * x;
130       }
131       column[j] -= column[i] * x;
132     }
133   }
134 
135   for (i = n - 1; i > 0; i--) {
136     for (j = i - 1; j >= 0; j--) {
137       column[j] -= matrix[j * n + i] * column[i];
138       matrix[j * n + i] = 0;
139     }
140   }
141 }
142 
143 #ifdef unused
144 void
schro_encoder_set_default_subband_weights(SchroEncoder * encoder)145 schro_encoder_set_default_subband_weights (SchroEncoder * encoder)
146 {
147   schro_encoder_calculate_subband_weights (encoder,
148       schro_encoder_perceptual_weight_constant);
149 }
150 #endif
151 
152 void
schro_encoder_calculate_subband_weights(SchroEncoder * encoder,double (* perceptual_weight)(double))153 schro_encoder_calculate_subband_weights (SchroEncoder * encoder,
154     double (*perceptual_weight) (double))
155 {
156   int wavelet;
157   int n_levels;
158   double *matrix_intra;
159   double *matrix_inter;
160   int n;
161   int i, j;
162   double column_intra[SCHRO_LIMIT_SUBBANDS];
163   double column_inter[SCHRO_LIMIT_SUBBANDS];
164   double *weight_intra;
165   double *weight_inter;
166 
167   matrix_intra =
168       schro_malloc (sizeof (double) * SCHRO_LIMIT_SUBBANDS *
169       SCHRO_LIMIT_SUBBANDS);
170   matrix_inter =
171       schro_malloc (sizeof (double) * SCHRO_LIMIT_SUBBANDS *
172       SCHRO_LIMIT_SUBBANDS);
173   weight_intra = schro_malloc (sizeof (double) * CURVE_SIZE * CURVE_SIZE);
174   weight_inter = schro_malloc (sizeof (double) * CURVE_SIZE * CURVE_SIZE);
175 
176   for (j = 0; j < CURVE_SIZE; j++) {
177     for (i = 0; i < CURVE_SIZE; i++) {
178       double fv_intra =
179           j * encoder->cycles_per_degree_vert * (1.0 / CURVE_SIZE);
180       double fh_intra =
181           i * encoder->cycles_per_degree_horiz * (1.0 / CURVE_SIZE);
182 
183       double fv_inter =
184           j * encoder->cycles_per_degree_vert * (1.0 / CURVE_SIZE) *
185           encoder->magic_inter_cpd_scale;
186       double fh_inter =
187           i * encoder->cycles_per_degree_horiz * (1.0 / CURVE_SIZE) *
188           encoder->magic_inter_cpd_scale;
189 
190       weight_intra[j * CURVE_SIZE + i] =
191           perceptual_weight (sqrt (fv_intra * fv_intra + fh_intra * fh_intra));
192       weight_inter[j * CURVE_SIZE + i] =
193           perceptual_weight (sqrt (fv_intra * fv_inter + fh_inter * fh_inter));
194     }
195   }
196 
197   for (wavelet = 0; wavelet < SCHRO_N_WAVELETS; wavelet++) {
198     for (n_levels = 1; n_levels <= SCHRO_LIMIT_ENCODER_TRANSFORM_DEPTH;
199         n_levels++) {
200       const float *h_curve[SCHRO_LIMIT_SUBBANDS];
201       const float *v_curve[SCHRO_LIMIT_SUBBANDS];
202       int hi[SCHRO_LIMIT_SUBBANDS];
203       int vi[SCHRO_LIMIT_SUBBANDS];
204 
205       n = 3 * n_levels + 1;
206 
207       for (i = 0; i < n; i++) {
208         int position = schro_subband_get_position (i);
209         int n_transforms;
210 
211         n_transforms = n_levels - SCHRO_SUBBAND_SHIFT (position);
212         if (position & 1) {
213           hi[i] = (n_transforms - 1) * 2;
214         } else {
215           hi[i] = (n_transforms - 1) * 2 + 1;
216         }
217         if (position & 2) {
218           vi[i] = (n_transforms - 1) * 2;
219         } else {
220           vi[i] = (n_transforms - 1) * 2 + 1;
221         }
222         h_curve[i] = schro_tables_wavelet_noise_curve[wavelet][hi[i]];
223         v_curve[i] = schro_tables_wavelet_noise_curve[wavelet][vi[i]];
224       }
225 
226       if (0) {
227         for (i = 0; i < n; i++) {
228           column_intra[i] = weighted_sum (h_curve[i], v_curve[i], weight_intra);
229           column_inter[i] = weighted_sum (h_curve[i], v_curve[i], weight_inter);
230           matrix_intra[i * n + i] = dot_product (h_curve[i], v_curve[i],
231               h_curve[i], v_curve[i], weight_intra);
232           matrix_inter[i * n + i] = dot_product (h_curve[i], v_curve[i],
233               h_curve[i], v_curve[i], weight_inter);
234           for (j = i + 1; j < n; j++) {
235             matrix_intra[i * n + j] = dot_product (h_curve[i], v_curve[i],
236                 h_curve[j], v_curve[j], weight_intra);
237             matrix_intra[j * n + i] = matrix_intra[i * n + j];
238             matrix_inter[i * n + j] = dot_product (h_curve[i], v_curve[i],
239                 h_curve[j], v_curve[j], weight_inter);
240             matrix_inter[j * n + i] = matrix_inter[i * n + j];
241           }
242         }
243 
244         solve (matrix_intra, column_intra, n);
245         solve (matrix_inter, column_inter, n);
246 
247         for (i = 0; i < n; i++) {
248           if (column_intra[i] < 0 || column_inter[i] < 0) {
249             SCHRO_ERROR ("BROKEN wavelet %d n_levels %d", wavelet, n_levels);
250             break;
251           }
252         }
253 
254         SCHRO_DEBUG ("wavelet %d n_levels %d", wavelet, n_levels);
255         for (i = 0; i < n; i++) {
256           SCHRO_DEBUG ("%g", 1.0 / sqrt (column_intra[i]));
257           SCHRO_DEBUG ("%g", 1.0 / sqrt (column_inter[i]));
258           encoder->intra_subband_weights[wavelet][n_levels - 1][i] =
259               sqrt (column_intra[i]);
260           encoder->inter_subband_weights[wavelet][n_levels - 1][i] =
261               sqrt (column_inter[i]);
262         }
263       } else {
264         for (i = 0; i < n; i++) {
265           int position = schro_subband_get_position (i);
266           int n_transforms;
267           double size;
268 
269           n_transforms = n_levels - SCHRO_SUBBAND_SHIFT (position);
270           size = (1.0 / CURVE_SIZE) * (1 << n_transforms);
271           encoder->intra_subband_weights[wavelet][n_levels - 1][i] =
272               1.0 / (size * sqrt (weighted_sum (h_curve[i], v_curve[i],
273                       weight_intra)));
274           encoder->inter_subband_weights[wavelet][n_levels - 1][i] =
275               1.0 / (size * sqrt (weighted_sum (h_curve[i], v_curve[i],
276                       weight_inter)));
277         }
278       }
279     }
280   }
281 
282   schro_free (weight_intra);
283   schro_free (matrix_intra);
284   schro_free (weight_inter);
285   schro_free (matrix_inter);
286 }
287 
288 void
schro_encoder_choose_quantisers(SchroEncoderFrame * frame)289 schro_encoder_choose_quantisers (SchroEncoderFrame * frame)
290 {
291   switch (frame->encoder->quantiser_engine) {
292     case SCHRO_QUANTISER_ENGINE_SIMPLE:
293       schro_encoder_choose_quantisers_simple (frame);
294       break;
295     case SCHRO_QUANTISER_ENGINE_RDO_BIT_ALLOCATION:
296       schro_encoder_choose_quantisers_rdo_bit_allocation (frame);
297       break;
298     case SCHRO_QUANTISER_ENGINE_CBR:
299       schro_encoder_choose_quantisers_rdo_cbr (frame);
300       break;
301     case SCHRO_QUANTISER_ENGINE_LOSSLESS:
302       schro_encoder_choose_quantisers_lossless (frame);
303       break;
304     case SCHRO_QUANTISER_ENGINE_LOWDELAY:
305       schro_encoder_choose_quantisers_lowdelay (frame);
306       break;
307     case SCHRO_QUANTISER_ENGINE_RDO_LAMBDA:
308       schro_encoder_choose_quantisers_rdo_lambda (frame);
309       break;
310     case SCHRO_QUANTISER_ENGINE_CONSTANT_ERROR:
311       schro_encoder_choose_quantisers_constant_error (frame);
312       break;
313     default:
314       SCHRO_ASSERT (0);
315   }
316 }
317 
318 void
schro_encoder_choose_quantisers_lossless(SchroEncoderFrame * frame)319 schro_encoder_choose_quantisers_lossless (SchroEncoderFrame * frame)
320 {
321   int i;
322   int component;
323 
324   for (component = 0; component < 3; component++) {
325     for (i = 0; i < 1 + 3 * frame->params.transform_depth; i++) {
326       schro_encoder_frame_set_quant_index (frame, component, i, -1, -1, 0);
327     }
328   }
329 }
330 
331 void
schro_encoder_choose_quantisers_simple(SchroEncoderFrame * frame)332 schro_encoder_choose_quantisers_simple (SchroEncoderFrame * frame)
333 {
334   SchroParams *params = &frame->params;
335   int i;
336   int component;
337   double noise_amplitude;
338   double a;
339   double max;
340   double *table;
341   double range;
342 
343   /* FIXME this should really be based directly on excursion */
344   range = (1<<frame->encoder->bit_depth) - 1.0;
345 
346   noise_amplitude = range * pow (0.1, frame->encoder->noise_threshold * 0.05);
347   SCHRO_DEBUG ("noise %g", noise_amplitude);
348 
349   if (frame->num_refs == 0) {
350     table = frame->encoder->intra_subband_weights[params->wavelet_filter_index]
351         [MAX (0, params->transform_depth - 1)];
352   } else {
353     table = frame->encoder->inter_subband_weights[params->wavelet_filter_index]
354         [MAX (0, params->transform_depth - 1)];
355   }
356 
357   for (component = 0; component < 3; component++) {
358     for (i = 0; i < 1 + 3 * params->transform_depth; i++) {
359       a = noise_amplitude * table[i];
360 
361       schro_encoder_frame_set_quant_index (frame, component, i, -1, -1,
362           schro_utils_multiplier_to_quant_index (a));
363     }
364   }
365 
366   max = 1.0;
367 
368   for (i = 0; i < 1 + 3 * params->transform_depth; i++) {
369     params->quant_matrix[i] =
370         schro_utils_multiplier_to_quant_index (max / table[i]);
371     SCHRO_DEBUG ("%g %g %d", table[i], max / table[i], params->quant_matrix[i]);
372   }
373 }
374 
375 void
schro_encoder_choose_quantisers_lowdelay(SchroEncoderFrame * frame)376 schro_encoder_choose_quantisers_lowdelay (SchroEncoderFrame * frame)
377 {
378   SchroParams *params = &frame->params;
379   int i;
380   int component;
381   int base;
382   const int *table;
383 
384   /* completely made up */
385   base = 12 + (30 - frame->encoder->noise_threshold) / 2;
386 
387   table = schro_tables_lowdelay_quants[params->wavelet_filter_index]
388       [MAX (0, params->transform_depth - 1)];
389 
390   for (component = 0; component < 3; component++) {
391     schro_encoder_frame_set_quant_index (frame, component, 0, -1, -1,
392         base - table[0]);
393 
394     for (i = 0; i < params->transform_depth; i++) {
395       schro_encoder_frame_set_quant_index (frame, component, 1 + 3 * i + 0, -1,
396           -1, base - table[1 + 2 * i + 0]);
397       schro_encoder_frame_set_quant_index (frame, component, 1 + 3 * i + 1, -1,
398           -1, base - table[1 + 2 * i + 0]);
399       schro_encoder_frame_set_quant_index (frame, component, 1 + 3 * i + 2, -1,
400           -1, base - table[1 + 2 * i + 1]);
401     }
402   }
403 
404 }
405 
406 #ifdef DUMP_SUBBAND_CURVES
407 static double
pow2(double x)408 pow2 (double x)
409 {
410   return x * x;
411 }
412 #endif
413 
414 #ifdef DUMP_SUBBAND_CURVES
415 static double
measure_error_subband(SchroEncoderFrame * frame,int component,int index,int quant_index)416 measure_error_subband (SchroEncoderFrame * frame, int component, int index,
417     int quant_index)
418 {
419   int i;
420   int j;
421   int16_t *line;
422   int skip = 1;
423   double error = 0;
424   int q;
425   int quant_factor;
426   int quant_offset;
427   int value;
428   int position;
429   SchroFrameData fd;
430 
431   position = schro_subband_get_position (index);
432   schro_subband_get_frame_data (&fd, frame->iwt_frame, component, position,
433       &frame->params);
434 
435   quant_factor = schro_table_quant[quant_index];
436   if (frame->params.num_refs > 0) {
437     quant_offset = schro_table_offset_3_8[quant_index];
438   } else {
439     quant_offset = schro_table_offset_1_2[quant_index];
440   }
441 
442   error = 0;
443   if (index == 0) {
444     for (j = 0; j < fd.height; j += skip) {
445       line = SCHRO_FRAME_DATA_GET_LINE (&fd, j);
446       for (i = 1; i < fd.width; i += skip) {
447         q = schro_quantise (abs (line[i] - line[i - 1]), quant_factor,
448             quant_offset);
449         value = schro_dequantise (q, quant_factor, quant_offset);
450         error += pow2 (value - abs (line[i] - line[i - 1]));
451       }
452     }
453   } else {
454     for (j = 0; j < fd.height; j += skip) {
455       line = SCHRO_FRAME_DATA_GET_LINE (&fd, j);
456       for (i = 0; i < fd.width; i += skip) {
457         q = schro_quantise (line[i], quant_factor, quant_offset);
458         value = schro_dequantise (q, quant_factor, quant_offset);
459         error += pow2 (value - line[i]);
460       }
461     }
462   }
463   error *= skip * skip;
464 
465   return error;
466 }
467 #endif
468 
469 typedef struct _ErrorFuncInfo ErrorFuncInfo;
470 struct _ErrorFuncInfo
471 {
472   int quant_factor;
473   int quant_offset;
474   double power;
475 };
476 
477 static double
error_pow(int x,void * priv)478 error_pow (int x, void *priv)
479 {
480   ErrorFuncInfo *efi = priv;
481   int q;
482   int value;
483   int y;
484 
485   q = schro_quantise (x, efi->quant_factor, efi->quant_offset);
486   value = schro_dequantise (q, efi->quant_factor, efi->quant_offset);
487 
488   y = abs (value - x);
489 
490   return pow (y, efi->power);
491 }
492 
493 void
schro_encoder_init_error_tables(SchroEncoder * encoder)494 schro_encoder_init_error_tables (SchroEncoder * encoder)
495 {
496   int i;
497 
498   for (i = 0; i < 60; i++) {
499     ErrorFuncInfo efi;
500 
501     efi.quant_factor = schro_table_quant[i];
502     efi.quant_offset = schro_table_offset_1_2[i];
503     efi.power = encoder->magic_error_power;
504 
505     schro_histogram_table_generate (encoder->intra_hist_tables + i,
506         error_pow, &efi);
507   }
508 
509 }
510 
511 #ifdef unused
512 static double
schro_histogram_estimate_error(SchroHistogram * hist,int quant_index,int num_refs)513 schro_histogram_estimate_error (SchroHistogram * hist, int quant_index,
514     int num_refs)
515 {
516   SchroHistogramTable *table;
517 
518   if (num_refs == 0) {
519     table =
520         (SchroHistogramTable
521         *) (schro_table_error_hist_shift3_1_2[quant_index]);
522   } else {
523     /* FIXME the 3/8 table doesn't exist yet */
524     //table = (SchroHistogramTable *)(schro_table_error_hist_shift3_3_8[quant_index]);
525     table =
526         (SchroHistogramTable
527         *) (schro_table_error_hist_shift3_1_2[quant_index]);
528   }
529   return schro_histogram_apply_table (hist, table);
530 }
531 #endif
532 
533 #ifdef DUMP_SUBBAND_CURVES
534 static double
schro_encoder_estimate_subband_arith(SchroEncoderFrame * frame,int component,int index,int quant_index)535 schro_encoder_estimate_subband_arith (SchroEncoderFrame * frame, int component,
536     int index, int quant_index)
537 {
538   int i;
539   int j;
540   int16_t *line;
541   int q;
542   int quant_factor;
543   int quant_offset;
544   int estimated_entropy;
545   SchroArith *arith;
546   int position;
547   SchroFrameData fd;
548 
549   arith = schro_arith_new ();
550   schro_arith_estimate_init (arith);
551 
552   position = schro_subband_get_position (index);
553   schro_subband_get_frame_data (&fd, frame->iwt_frame, component, position,
554       &frame->params);
555 
556   quant_factor = schro_table_quant[quant_index];
557   quant_offset = schro_table_offset_1_2[quant_index];
558 
559   if (index == 0) {
560     for (j = 0; j < fd.height; j++) {
561       line = SCHRO_FRAME_DATA_GET_LINE (&fd, j);
562       schro_arith_estimate_sint (arith,
563           SCHRO_CTX_ZPZN_F1, SCHRO_CTX_COEFF_DATA, SCHRO_CTX_SIGN_ZERO, 0);
564       for (i = 1; i < fd.width; i++) {
565         q = schro_quantise (line[i] - line[i - 1], quant_factor, quant_offset);
566         schro_arith_estimate_sint (arith,
567             SCHRO_CTX_ZPZN_F1, SCHRO_CTX_COEFF_DATA, SCHRO_CTX_SIGN_ZERO, q);
568       }
569     }
570   } else {
571     for (j = 0; j < fd.height; j++) {
572       line = SCHRO_FRAME_DATA_GET_LINE (&fd, j);
573       for (i = 0; i < fd.width; i++) {
574         q = schro_quantise (line[i], quant_factor, quant_offset);
575         schro_arith_estimate_sint (arith,
576             SCHRO_CTX_ZPZN_F1, SCHRO_CTX_COEFF_DATA, SCHRO_CTX_SIGN_ZERO, q);
577       }
578     }
579   }
580 
581   estimated_entropy = 0;
582 
583   estimated_entropy += arith->contexts[SCHRO_CTX_ZPZN_F1].n_bits;
584   estimated_entropy += arith->contexts[SCHRO_CTX_ZP_F2].n_bits;
585   estimated_entropy += arith->contexts[SCHRO_CTX_ZP_F3].n_bits;
586   estimated_entropy += arith->contexts[SCHRO_CTX_ZP_F4].n_bits;
587   estimated_entropy += arith->contexts[SCHRO_CTX_ZP_F5].n_bits;
588   estimated_entropy += arith->contexts[SCHRO_CTX_ZP_F6p].n_bits;
589 
590   estimated_entropy += arith->contexts[SCHRO_CTX_COEFF_DATA].n_bits;
591 
592   estimated_entropy += arith->contexts[SCHRO_CTX_SIGN_ZERO].n_bits;
593 
594   schro_arith_free (arith);
595 
596   return estimated_entropy;
597 }
598 #endif
599 
600 static void
schro_encoder_generate_subband_histogram(SchroEncoderFrame * frame,int component,int index,SchroHistogram * hist,int skip)601 schro_encoder_generate_subband_histogram (SchroEncoderFrame * frame,
602     int component, int index, SchroHistogram * hist, int skip)
603 {
604   int position;
605   SchroFrameData fd;
606 
607   position = schro_subband_get_position (index);
608   schro_subband_get_frame_data (&fd, frame->iwt_frame, component, position,
609       &frame->params);
610 
611   if (index == 0 && frame->num_refs == 0) {
612     schro_frame_data_generate_histogram_dc_predict (&fd, hist, skip, 0, 0);
613   } else {
614     schro_frame_data_generate_histogram (&fd, hist, skip);
615   }
616 }
617 
618 static void
schro_encoder_generate_subband_histograms(SchroEncoderFrame * frame)619 schro_encoder_generate_subband_histograms (SchroEncoderFrame * frame)
620 {
621   SchroParams *params = &frame->params;
622   int i;
623   int component;
624   int pos;
625   int skip;
626 
627   for (component = 0; component < 3; component++) {
628     for (i = 0; i < 1 + 3 * params->transform_depth; i++) {
629       pos = schro_subband_get_position (i);
630       skip = 1 << MAX (0, SCHRO_SUBBAND_SHIFT (pos) - 1);
631       schro_encoder_generate_subband_histogram (frame, component, i,
632           &frame->subband_hists[component][i], skip);
633     }
634   }
635 
636   frame->have_histograms = TRUE;
637 }
638 
639 #ifdef DUMP_SUBBAND_CURVES
640 static void
schro_encoder_dump_subband_curves(SchroEncoderFrame * frame)641 schro_encoder_dump_subband_curves (SchroEncoderFrame * frame)
642 {
643   SchroParams *params = &frame->params;
644   int i;
645   int component;
646   int position;
647 
648   SCHRO_ASSERT (frame->have_histograms);
649 
650   for (component = 0; component < 3; component++) {
651     for (i = 0; i < 1 + 3 * params->transform_depth; i++) {
652       int vol;
653       SchroFrameData fd;
654       int j;
655       SchroHistogram *hist;
656 
657       position = schro_subband_get_position (i);
658       schro_subband_get_frame_data (&fd, frame->iwt_frame, component,
659           position, &frame->params);
660       vol = fd.width * fd.height;
661       hist = &frame->subband_hists[component][i];
662 
663       for (j = 0; j < 60; j++) {
664         double est_entropy;
665         double error;
666         double est_error;
667         double arith_entropy;
668 
669         error = measure_error_subband (frame, component, i, j);
670         est_entropy =
671             schro_histogram_estimate_entropy (&frame->subband_hists[component]
672             [i], j, params->is_noarith);
673         est_error =
674             schro_histogram_apply_table (hist,
675             &frame->encoder->intra_hist_tables[j]);
676         arith_entropy =
677             schro_encoder_estimate_subband_arith (frame, component, i, j);
678 
679         schro_dump (SCHRO_DUMP_SUBBAND_CURVE, "%d %d %d %g %g %g %g\n",
680             component, i, j, est_entropy / vol, arith_entropy / vol,
681             est_error / vol, error / vol);
682       }
683     }
684   }
685 }
686 #endif
687 
688 static void
schro_encoder_calc_estimates(SchroEncoderFrame * frame)689 schro_encoder_calc_estimates (SchroEncoderFrame * frame)
690 {
691   SchroParams *params = &frame->params;
692   int i;
693   int j;
694   int component;
695   double *arith_context_ratios;
696 
697   SCHRO_ASSERT (frame->have_histograms);
698 
699 #ifdef DUMP_SUBBAND_CURVES
700   schro_encoder_dump_subband_curves (frame);
701 #endif
702 
703   for (component = 0; component < 3; component++) {
704     if (frame->num_refs == 0) {
705       arith_context_ratios =
706           frame->encoder->average_arith_context_ratios_intra[component];
707     } else {
708       arith_context_ratios =
709           frame->encoder->average_arith_context_ratios_inter[component];
710     }
711 
712     for (i = 0; i < 1 + 3 * params->transform_depth; i++) {
713       for (j = 0; j < 60; j++) {
714         int position;
715         SchroHistogram *hist;
716         SchroFrameData fd;
717 
718         position = schro_subband_get_position (i);
719         schro_subband_get_frame_data (&fd, frame->iwt_frame, component,
720             position, &frame->params);
721 
722         hist = &frame->subband_hists[component][i];
723         frame->est_entropy[component][i][j] =
724             schro_histogram_estimate_entropy (hist, j, params->is_noarith);
725         frame->est_entropy[component][i][j] *= arith_context_ratios[i];
726         frame->est_error[component][i][j] = schro_histogram_apply_table (hist,
727             &frame->encoder->intra_hist_tables[j]);
728       }
729     }
730   }
731   frame->have_estimate_tables = TRUE;
732 }
733 
734 /*
735  * Quantiser engine which picks the best RDO quantisers to fit in
736  * a frame bit allocation.
737  */
738 void
schro_encoder_choose_quantisers_rdo_bit_allocation(SchroEncoderFrame * frame)739 schro_encoder_choose_quantisers_rdo_bit_allocation (SchroEncoderFrame * frame)
740 {
741   double frame_lambda;
742   int bits;
743 
744   schro_encoder_generate_subband_histograms (frame);
745   schro_encoder_calc_estimates (frame);
746 
747   SCHRO_ASSERT (frame->have_estimate_tables);
748 
749   bits = frame->allocated_residual_bits;
750 
751   frame_lambda = schro_encoder_entropy_to_lambda (frame, bits);
752   frame->frame_lambda = frame_lambda;
753   SCHRO_DEBUG ("LAMBDA: %d %g %d", frame->frame_number, frame_lambda, bits);
754 
755   schro_encoder_lambda_to_entropy (frame, frame_lambda);
756 }
757 
758 void
schro_encoder_choose_quantisers_rdo_lambda(SchroEncoderFrame * frame)759 schro_encoder_choose_quantisers_rdo_lambda (SchroEncoderFrame * frame)
760 {
761   SCHRO_DEBUG ("Using rdo_lambda quant selection on frame %d with lambda %g",
762       frame->frame_number, frame->frame_lambda);
763 
764   schro_encoder_generate_subband_histograms (frame);
765   schro_encoder_calc_estimates (frame);
766 
767   SCHRO_ASSERT (frame->have_estimate_tables);
768 
769   schro_encoder_lambda_to_entropy (frame, frame->frame_lambda);
770 }
771 
772 
773 void
schro_encoder_choose_quantisers_rdo_cbr(SchroEncoderFrame * frame)774 schro_encoder_choose_quantisers_rdo_cbr (SchroEncoderFrame * frame)
775 {
776   schro_encoder_generate_subband_histograms (frame);
777   schro_encoder_calc_estimates (frame);
778 
779   SCHRO_ASSERT (frame->have_estimate_tables);
780 
781   schro_encoder_lambda_to_entropy (frame, frame->frame_lambda);
782 }
783 
784 void
schro_encoder_estimate_entropy(SchroEncoderFrame * frame)785 schro_encoder_estimate_entropy (SchroEncoderFrame * frame)
786 {
787   SchroParams *params = &frame->params;
788   int i;
789   int component;
790   int n = 0;
791 
792   for (component = 0; component < 3; component++) {
793     for (i = 0; i < 1 + 3 * params->transform_depth; i++) {
794       n += frame->
795           est_entropy[component][i][frame->quant_indices[component][i][0]];
796     }
797   }
798   frame->estimated_residual_bits = n;
799 
800   if (frame->allocated_residual_bits > 0 &&
801       frame->estimated_residual_bits >
802       2 * frame->encoder->bits_per_picture + frame->allocated_residual_bits) {
803     SCHRO_WARNING ("%d: estimated entropy too big (%d vs %d)",
804         frame->frame_number,
805         frame->estimated_residual_bits, frame->allocated_residual_bits);
806   }
807 }
808 
809 static int
schro_subband_pick_quant(SchroEncoderFrame * frame,int component,int i,double lambda)810 schro_subband_pick_quant (SchroEncoderFrame * frame, int component, int i,
811     double lambda)
812 {
813   double x;
814   double min;
815   int j;
816   int j_min;
817   double entropy;
818   double error;
819 
820   SCHRO_ASSERT (frame->have_estimate_tables);
821 
822   j_min = -1;
823   min = 0;
824   for (j = 0; j < 60; j++) {
825     entropy = frame->est_entropy[component][i][j];
826     error = frame->est_error[component][i][j];
827 
828     x = entropy + lambda * error;
829     if (j == 0 || x < min) {
830       j_min = j;
831       min = x;
832     }
833   }
834 
835   return j_min;
836 }
837 
838 static double
schro_encoder_lambda_to_entropy(SchroEncoderFrame * frame,double frame_lambda)839 schro_encoder_lambda_to_entropy (SchroEncoderFrame * frame, double frame_lambda)
840 {
841   SchroParams *params = &frame->params;
842   int i;
843   int component;
844   double entropy = 0;
845   double *table;
846 
847   if (frame->num_refs == 0) {
848     table = frame->encoder->intra_subband_weights[params->wavelet_filter_index]
849         [MAX (0, params->transform_depth - 1)];
850   } else {
851     table = frame->encoder->inter_subband_weights[params->wavelet_filter_index]
852         [MAX (0, params->transform_depth - 1)];
853   }
854 
855   for (component = 0; component < 3; component++) {
856     for (i = 0; i < 1 + 3 * params->transform_depth; i++) {
857       double lambda;
858       double weight;
859       int quant_index;
860       int position = schro_subband_get_position (i);
861 
862       lambda = frame_lambda;
863 
864       if (i == 0) {
865         lambda *= frame->encoder->magic_subband0_lambda_scale;
866       }
867       if (component > 0) {
868         lambda *= frame->encoder->magic_chroma_lambda_scale;
869       }
870       if (SCHRO_SUBBAND_IS_DIAGONALLY_ORIENTED (position)) {
871         lambda *= frame->encoder->magic_diagonal_lambda_scale;
872       }
873 
874       weight = table[i];
875       lambda /= weight * weight;
876 
877       quant_index = schro_subband_pick_quant (frame, component, i, lambda);
878       entropy += frame->est_entropy[component][i][quant_index];
879       schro_encoder_frame_set_quant_index (frame, component, i, -1, -1,
880           quant_index);
881     }
882   }
883 
884   return entropy;
885 }
886 
887 double
schro_encoder_entropy_to_lambda(SchroEncoderFrame * frame,double entropy)888 schro_encoder_entropy_to_lambda (SchroEncoderFrame * frame, double entropy)
889 {
890   int j;
891   double lambda_hi, lambda_lo, lambda_mid;
892   double entropy_hi, entropy_lo, entropy_mid;
893 
894   lambda_hi = 1;
895   entropy_hi = schro_encoder_lambda_to_entropy (frame, lambda_hi);
896   SCHRO_DEBUG ("start target=%g lambda=%g entropy=%g",
897       entropy, lambda_hi, entropy_hi);
898 
899   if (entropy_hi < entropy) {
900     entropy_lo = entropy_hi;
901     lambda_lo = lambda_hi;
902 
903     for (j = 0; j < 5; j++) {
904       lambda_hi = lambda_lo * 100;
905       entropy_hi = schro_encoder_lambda_to_entropy (frame, lambda_hi);
906 
907       SCHRO_DEBUG ("have: lambda=[%g,%g] entropy=[%g,%g] target=%g",
908           lambda_lo, lambda_hi, entropy_lo, entropy_hi, entropy);
909       if (entropy_hi > entropy)
910         break;
911 
912       SCHRO_DEBUG ("--> step up");
913 
914       entropy_lo = entropy_hi;
915       lambda_lo = lambda_hi;
916     }
917     SCHRO_DEBUG ("--> stopping");
918   } else {
919     for (j = 0; j < 5; j++) {
920       lambda_lo = lambda_hi * 0.01;
921       entropy_lo = schro_encoder_lambda_to_entropy (frame, lambda_lo);
922 
923       SCHRO_DEBUG ("have: lambda=[%g,%g] entropy=[%g,%g] target=%g",
924           lambda_lo, lambda_hi, entropy_lo, entropy_hi, entropy);
925 
926       SCHRO_DEBUG ("--> step down");
927       if (entropy_lo < entropy)
928         break;
929 
930       entropy_hi = entropy_lo;
931       lambda_hi = lambda_lo;
932     }
933     SCHRO_DEBUG ("--> stopping");
934   }
935   if (entropy_lo == entropy_hi) {
936     return sqrt (lambda_lo * lambda_hi);
937   }
938 
939   if (entropy_lo > entropy || entropy_hi < entropy) {
940     SCHRO_ERROR ("entropy not bracketed");
941   }
942 
943   for (j = 0; j < 7; j++) {
944     if (entropy_hi == entropy_lo)
945       break;
946 
947     SCHRO_DEBUG ("have: lambda=[%g,%g] entropy=[%g,%g] target=%g",
948         lambda_lo, lambda_hi, entropy_lo, entropy_hi, entropy);
949 
950     lambda_mid = sqrt (lambda_lo * lambda_hi);
951     entropy_mid = schro_encoder_lambda_to_entropy (frame, lambda_mid);
952 
953     SCHRO_DEBUG ("picking lambda_mid=%g entropy=%g", lambda_mid, entropy_mid);
954 
955     if (entropy_mid > entropy) {
956       lambda_hi = lambda_mid;
957       entropy_hi = entropy_mid;
958       SCHRO_DEBUG ("--> focus up");
959     } else {
960       lambda_lo = lambda_mid;
961       entropy_lo = entropy_mid;
962       SCHRO_DEBUG ("--> focus down");
963     }
964   }
965 
966   lambda_mid = sqrt (lambda_hi * lambda_lo);
967   SCHRO_DEBUG ("done %g", lambda_mid);
968   return lambda_mid;
969 }
970 
971 static double
schro_encoder_lambda_to_error(SchroEncoderFrame * frame,double frame_lambda)972 schro_encoder_lambda_to_error (SchroEncoderFrame * frame, double frame_lambda)
973 {
974   SchroParams *params = &frame->params;
975   int i;
976   int component;
977   double error = 0;
978   double *table;
979 
980   if (frame->num_refs == 0) {
981     table = frame->encoder->intra_subband_weights[params->wavelet_filter_index]
982         [MAX (0, params->transform_depth - 1)];
983   } else {
984     table = frame->encoder->inter_subband_weights[params->wavelet_filter_index]
985         [MAX (0, params->transform_depth - 1)];
986   }
987 
988   for (component = 0; component < 3; component++) {
989     for (i = 0; i < 1 + 3 * params->transform_depth; i++) {
990       double lambda;
991       double weight;
992       int quant_index;
993 
994       lambda = frame_lambda;
995 
996       if (i == 0) {
997         lambda *= frame->encoder->magic_subband0_lambda_scale;
998       }
999       if (component > 0) {
1000         lambda *= frame->encoder->magic_chroma_lambda_scale;
1001       }
1002 
1003       weight = table[i];
1004       lambda /= weight * weight;
1005 
1006       quant_index = schro_subband_pick_quant (frame, component, i, lambda);
1007       error += frame->est_error[component][i][quant_index];
1008       schro_encoder_frame_set_quant_index (frame, component, i, -1, -1,
1009           quant_index);
1010     }
1011   }
1012 
1013   return error;
1014 }
1015 
1016 static double
schro_encoder_error_to_lambda(SchroEncoderFrame * frame,double error)1017 schro_encoder_error_to_lambda (SchroEncoderFrame * frame, double error)
1018 {
1019   int j;
1020   double lambda_hi, lambda_lo, lambda_mid;
1021   double error_hi, error_lo, error_mid;
1022 
1023   lambda_lo = 1;
1024   error_lo = schro_encoder_lambda_to_error (frame, lambda_lo);
1025   SCHRO_DEBUG ("start target=%g lambda=%g error=%g",
1026       error, lambda_lo, error_lo, lambda_lo, error);
1027 
1028   if (error < error_lo) {
1029     error_hi = error_lo;
1030     lambda_hi = lambda_lo;
1031 
1032     for (j = 0; j < 5; j++) {
1033       lambda_lo = lambda_hi * 100;
1034       error_lo = schro_encoder_lambda_to_error (frame, lambda_lo);
1035 
1036       SCHRO_DEBUG ("have: lambda=[%g,%g] error=[%g,%g] target=%g",
1037           lambda_lo, lambda_hi, error_lo, error_hi, error);
1038       if (error > error_lo)
1039         break;
1040 
1041       SCHRO_DEBUG ("--> step up");
1042 
1043       error_hi = error_lo;
1044       lambda_hi = lambda_lo;
1045     }
1046     SCHRO_DEBUG ("--> stopping");
1047   } else {
1048     for (j = 0; j < 5; j++) {
1049       lambda_hi = lambda_lo * 0.01;
1050       error_hi = schro_encoder_lambda_to_error (frame, lambda_hi);
1051 
1052       SCHRO_DEBUG ("have: lambda=[%g,%g] error=[%g,%g] target=%g",
1053           lambda_lo, lambda_hi, error_lo, error_hi, error);
1054 
1055       SCHRO_DEBUG ("--> step down");
1056       if (error < error_hi)
1057         break;
1058 
1059       error_lo = error_hi;
1060       lambda_lo = lambda_hi;
1061     }
1062     SCHRO_DEBUG ("--> stopping");
1063   }
1064   if (error_lo == error_hi) {
1065     return sqrt (lambda_lo * lambda_hi);
1066   }
1067 
1068   if (error_lo > error || error_hi < error) {
1069     SCHRO_DEBUG ("error not bracketed");
1070   }
1071 
1072   for (j = 0; j < 14; j++) {
1073     if (error_hi == error_lo)
1074       break;
1075 
1076     SCHRO_DEBUG ("have: lambda=[%g,%g] error=[%g,%g] target=%g",
1077         lambda_lo, lambda_hi, error_lo, error_hi, error);
1078 
1079     lambda_mid = sqrt (lambda_lo * lambda_hi);
1080     error_mid = schro_encoder_lambda_to_error (frame, lambda_mid);
1081 
1082     SCHRO_DEBUG ("picking lambda_mid=%g error=%g", lambda_mid, error_mid);
1083 
1084     if (error_mid > error) {
1085       lambda_hi = lambda_mid;
1086       error_hi = error_mid;
1087       SCHRO_DEBUG ("--> focus up");
1088     } else {
1089       lambda_lo = lambda_mid;
1090       error_lo = error_mid;
1091       SCHRO_DEBUG ("--> focus down");
1092     }
1093   }
1094 
1095   lambda_mid = sqrt (lambda_hi * lambda_lo);
1096   SCHRO_DEBUG ("done %g", lambda_mid);
1097   return lambda_mid;
1098 }
1099 
1100 void
schro_encoder_choose_quantisers_constant_error(SchroEncoderFrame * frame)1101 schro_encoder_choose_quantisers_constant_error (SchroEncoderFrame * frame)
1102 {
1103   double frame_lambda;
1104   double error;
1105 
1106   schro_encoder_generate_subband_histograms (frame);
1107   schro_encoder_calc_estimates (frame);
1108 
1109   SCHRO_ASSERT (frame->have_estimate_tables);
1110 
1111   error = 255.0 * pow (0.1, frame->encoder->noise_threshold * 0.05);
1112   error *= frame->params.video_format->width *
1113       frame->params.video_format->height;
1114 
1115   frame_lambda = schro_encoder_error_to_lambda (frame, error);
1116 
1117   frame->frame_lambda = frame_lambda;
1118   SCHRO_DEBUG ("LAMBDA: %d %g", frame->frame_number, frame_lambda);
1119 }
1120