Panda3D
fftCompressor.cxx
Go to the documentation of this file.
1 /**
2  * PANDA 3D SOFTWARE
3  * Copyright (c) Carnegie Mellon University. All rights reserved.
4  *
5  * All use of this software is subject to the terms of the revised BSD
6  * license. You should have received a copy of this license along
7  * with this source code in a file named "LICENSE."
8  *
9  * @file fftCompressor.cxx
10  * @author drose
11  * @date 2000-12-11
12  */
13 
14 #include "fftCompressor.h"
15 #include "config_mathutil.h"
16 #include "config_linmath.h"
17 
18 #include "datagram.h"
19 #include "datagramIterator.h"
20 #include "compose_matrix.h"
21 #include "pmap.h"
22 #include <math.h>
23 
24 #ifdef HAVE_FFTW
25 
26 // hack..... this is a hack to help interrogate sort out a macro in the system
27 // poll and select definitions
28 #ifdef howmany
29 #undef howmany
30 #endif
31 
32 #include <fftw3.h>
33 
34 // These FFTW support objects can only be defined if we actually have the FFTW
35 // library available.
36 static fftw_plan get_real_compress_plan(int length);
37 static fftw_plan get_real_decompress_plan(int length);
38 
39 typedef pmap<int, fftw_plan> RealPlans;
40 static RealPlans _real_compress_plans;
41 static RealPlans _real_decompress_plans;
42 
43 #endif
44 
45 /**
46  * Constructs a new compressor object with default parameters.
47  */
49 FFTCompressor() {
50  _bam_minor_version = 0;
51  set_quality(-1);
52  _use_error_threshold = false;
53  _transpose_quats = false;
54 }
55 
56 /**
57  * Returns true if the FFTW library is compiled in, so that this class is
58  * actually capable of doing useful compression/decompression work. Returns
59  * false otherwise, in which case any attempt to write a compressed stream
60  * will actually write an uncompressed stream, and any attempt to read a
61  * compressed stream will fail.
62  */
65 #ifndef HAVE_FFTW
66  return false;
67 #else
68  return true;
69 #endif
70 }
71 
72 /**
73  * Sets the quality factor for the compression. This is an integer in the
74  * range 0 - 100 that roughly controls how aggressively the reals are
75  * compressed; lower numbers mean smaller output, and more data loss.
76  *
77  * There are a few special cases. Quality -1 means to use whatever individual
78  * parameters are set in the user's Configrc file, rather than the single
79  * quality dial. Quality 101 or higher means to generate lossless output
80  * (this is the default if libfftw is not available).
81  *
82  * Quality 102 writes all four components of quaternions to the output file,
83  * rather than just three, quality 103 converts hpr to matrix (instead of
84  * quat) and writes a 9-component matrix, and quality 104 just writes out hpr
85  * directly. Quality levels 102 and greater are strictly for debugging
86  * purposes, and are only available if NDEBUG is not defined.
87  */
89 set_quality(int quality) {
90 #ifndef HAVE_FFTW
91  // If we don't actually have FFTW, we can't really compress anything.
92  if (_quality <= 100) {
93  mathutil_cat.warning()
94  << "FFTW library is not available; generating uncompressed output.\n";
95  }
96  _quality = 101;
97 
98 #else
99  _quality = quality;
100 
101  if (_quality < 0) {
102  // A negative quality indicates we should read the various parameters from
103  // individual config variables.
104  _fft_offset = fft_offset;
105  _fft_factor = fft_factor;
106  _fft_exponent = fft_exponent;
107 
108  } else if (_quality < 40) {
109  // 0 - 40 : fft-offset 1.0 - 0.001 fft-factor 1.0 fft-exponent 4.0
110 
111  double t = (double)_quality / 40.0;
112  _fft_offset = interpolate(t, 1.0, 0.001);
113  _fft_factor = 1.0;
114  _fft_exponent = 4.0;
115 
116  } else if (_quality < 95) {
117  // 40 - 95: fft-offset 0.001 fft-factor 1.0 - 0.1 fft-exponent 4.0
118 
119  double t = (double)(_quality - 40) / 55.0;
120  _fft_offset = 0.001;
121  _fft_factor = interpolate(t, 1.0, 0.1);
122  _fft_exponent = 4.0;
123 
124  } else {
125  // 95 - 100: fft-offset 0.001 fft-factor 0.1 - 0.0 fft-exponent 4.0
126 
127  double t = (double)(_quality - 95) / 5.0;
128  _fft_offset = 0.001;
129  _fft_factor = interpolate(t, 0.1, 0.0);
130  _fft_exponent = 4.0;
131  }
132 #endif
133 }
134 
135 /**
136  * Returns the quality number that was previously set via set_quality().
137  */
139 get_quality() const {
140  return _quality;
141 }
142 
143 /**
144  * Enables or disables the use of the error threshold measurement to put a cap
145  * on the amount of damage done by lossy compression. When this is enabled,
146  * the potential results of the compression are analyzed before the data is
147  * written; if it is determined that the compression will damage a particular
148  * string of reals too much, that particular string of reals is written
149  * uncompressed.
150  */
152 set_use_error_threshold(bool use_error_threshold) {
153  _use_error_threshold = use_error_threshold;
154 }
155 
156 /**
157  * Returns whether the error threshold measurement is enabled. See
158  * set_use_error_threshold().
159  */
161 get_use_error_threshold() const {
162  return _use_error_threshold;
163 }
164 
165 /**
166  * Sets the transpose_quats flag. This is provided mainly for backward
167  * compatibility with old bam files that were written out with the quaternions
168  * inadvertently transposed.
169  */
171 set_transpose_quats(bool flag) {
172  _transpose_quats = flag;
173 }
174 
175 /**
176  * Returns the transpose_quats flag. See set_transpose_quats().
177  */
179 get_transpose_quats() const {
180  return _transpose_quats;
181 }
182 
183 /**
184  * Writes the compression parameters to the indicated datagram. It is
185  * necessary to call this before writing anything else to the datagram, since
186  * these parameters will be necessary to correctly decompress the data later.
187  */
189 write_header(Datagram &datagram) {
190  datagram.add_int8(_quality);
191  if (_quality < 0) {
192  datagram.add_float64(_fft_offset);
193  datagram.add_float64(_fft_factor);
194  datagram.add_float64(_fft_exponent);
195  }
196 }
197 
198 /**
199  * Writes an array of floating-point numbers to the indicated datagram.
200  */
202 write_reals(Datagram &datagram, const PN_stdfloat *array, int length) {
203  datagram.add_int32(length);
204 
205  if (_quality > 100) {
206  // Special case: lossless output.
207  for (int i = 0; i < length; i++) {
208  datagram.add_stdfloat(array[i]);
209  }
210  return;
211  }
212 
213 #ifndef HAVE_FFTW
214  // If we don't have FFTW, we shouldn't get here.
215  nassertv(false);
216 
217 #else
218 
219  if (length == 0) {
220  // Special case: do nothing.
221  return;
222  }
223 
224  if (length == 1) {
225  // Special case: just write out the one number.
226  datagram.add_stdfloat(array[0]);
227  return;
228  }
229 
230  // Normal case: FFT the array, and write that out.
231 
232  // First, check the compressability.
233  bool reject_compression = false;
234 
235  // This logic needs a closer examination. Not sure it's useful as-is.
236  /*
237  if (_use_error_threshold) {
238  // Don't encode the data if it moves too erratically.
239  PN_stdfloat error = get_compressability(array, length);
240  if (error > fft_error_threshold) {
241  // No good: the data probably won't compress well. Just write out
242  // lossless data.
243  reject_compression = true;
244  }
245  }
246  */
247 
248  datagram.add_bool(reject_compression);
249  if (reject_compression) {
250  if (mathutil_cat.is_debug()) {
251  mathutil_cat.debug()
252  << "Writing stream of " << length << " numbers uncompressed.\n";
253  }
254  for (int i = 0; i < length; i++) {
255  datagram.add_stdfloat(array[i]);
256  }
257  return;
258  }
259 
260  // Now generate the Fourier transform.
261  int fft_length = length / 2 + 1;
262  fftw_complex *fft_bins = (fftw_complex *)alloca(fft_length * sizeof(fftw_complex));
263 
264  // This is for an in-place transform. It doesn't violate strict aliasing
265  // rules because &fft_bins[0][0] is still a double pointer. This saves on
266  // precious stack space.
267  double *data = &fft_bins[0][0];
268  int i;
269  for (i = 0; i < length; i++) {
270  data[i] = array[i];
271  }
272 
273  // Note: This is an in-place DFT. `data` and `fft_bins` are aliases.
274  fftw_plan plan = get_real_compress_plan(length);
275  fftw_execute_dft_r2c(plan, data, fft_bins);
276 
277 
278  // Now encode the numbers, run-length encoded by size, so we only write out
279  // the number of bits we need for each number.
280  // Note that Panda3D has conventionally always used FFTW2's halfcomplex
281  // format for serializing the bins. In short, this means that for an n-length
282  // FFT, it stores:
283  // 1) The real components for bins 0 through floor(n/2), followed by...
284  // 2) The imaginary components for bins floor((n+1)/2)-1 through 1.
285  // (Imaginary component for bin 0 is never stored, as that's always zero.)
286 
287  vector_double run;
288  RunWidth run_width = RW_invalid;
289  int num_written = 0;
290 
291  for (i = 0; i < length; i++) {
292  static const double max_range_32 = 2147483647.0;
293  static const double max_range_16 = 32767.0;
294  static const double max_range_8 = 127.0;
295 
296  int bin; // which FFT bin we're storing
297  int j; // 0=real; 1=imag
298  if (i < fft_length) {
299  bin = i;
300  j = 0;
301  } else {
302  bin = length - i;
303  j = 1;
304  }
305 
306  double scale_factor = get_scale_factor(bin, fft_length);
307  double num = cfloor(fft_bins[bin][j] / scale_factor + 0.5);
308 
309  // How many bits do we need to encode this integer?
310  double a = fabs(num);
311  RunWidth num_width;
312 
313  if (a == 0.0) {
314  num_width = RW_0;
315 
316  } else if (a <= max_range_8) {
317  num_width = RW_8;
318 
319  } else if (a <= max_range_16) {
320  num_width = RW_16;
321 
322  } else if (a <= max_range_32) {
323  num_width = RW_32;
324 
325  } else {
326  num_width = RW_double;
327  }
328 
329  // A special case: if we're writing a string of one-byters and we come
330  // across a single intervening zero, don't interrupt the run just for
331  // that.
332  if (run_width == RW_8 && num_width == RW_0) {
333  if (run.back() != 0) {
334  num_width = RW_8;
335  }
336  }
337 
338  if (num_width != run_width) {
339  // Now we need to flush the last run.
340 
341  // First, however, take care of the special case above: if we're
342  // switching from RW_8 to RW_0, there could be a zero at the end, which
343  // should be reclaimed into the RW_0 run.
344  bool reclaimed_zero = (run_width == RW_8 && num_width == RW_0 &&
345  run.back() == 0);
346  if (reclaimed_zero) {
347  run.pop_back();
348  }
349 
350  num_written += write_run(datagram, run_width, run);
351  run.clear();
352  run_width = num_width;
353 
354  if (reclaimed_zero) {
355  run.push_back(0);
356  }
357  }
358 
359  run.push_back(num);
360  }
361 
362  num_written += write_run(datagram, run_width, run);
363  nassertv(num_written == length);
364 #endif
365 }
366 
367 /**
368  * Writes an array of HPR angles to the indicated datagram.
369  */
371 write_hprs(Datagram &datagram, const LVecBase3 *array, int length) {
372 #ifndef NDEBUG
373  if (_quality >= 104) {
374  // If quality level is at least 104, we don't even convert hpr at all.
375  // This is just for debugging.
376  vector_stdfloat h, p, r;
377 
378  h.reserve(length);
379  p.reserve(length);
380  r.reserve(length);
381 
382  for (int i = 0; i < length; i++) {
383  h.push_back(array[i][0]);
384  p.push_back(array[i][1]);
385  r.push_back(array[i][2]);
386  }
387 
388  if (length == 0) {
389  write_reals(datagram, nullptr, length);
390  write_reals(datagram, nullptr, length);
391  write_reals(datagram, nullptr, length);
392  } else {
393  write_reals(datagram, &h[0], length);
394  write_reals(datagram, &p[0], length);
395  write_reals(datagram, &r[0], length);
396  }
397  return;
398  }
399  if (_quality >= 103) {
400  // If quality level is 103, we convert hpr to a table of matrices. This
401  // is just for debugging.
402  vector_stdfloat
403  m00, m01, m02,
404  m10, m11, m12,
405  m20, m21, m22;
406  for (int i = 0; i < length; i++) {
407  LMatrix3 mat;
408  compose_matrix(mat, LVecBase3(1.0, 1.0, 1.0), LVecBase3(0.0, 0.0, 0.0),
409  array[i]);
410  m00.push_back(mat(0, 0));
411  m01.push_back(mat(0, 1));
412  m02.push_back(mat(0, 2));
413  m10.push_back(mat(1, 0));
414  m11.push_back(mat(1, 1));
415  m12.push_back(mat(1, 2));
416  m20.push_back(mat(2, 0));
417  m21.push_back(mat(2, 1));
418  m22.push_back(mat(2, 2));
419  }
420 
421  if (length == 0) {
422  write_reals(datagram, nullptr, length);
423  write_reals(datagram, nullptr, length);
424  write_reals(datagram, nullptr, length);
425  write_reals(datagram, nullptr, length);
426  write_reals(datagram, nullptr, length);
427  write_reals(datagram, nullptr, length);
428  write_reals(datagram, nullptr, length);
429  write_reals(datagram, nullptr, length);
430  write_reals(datagram, nullptr, length);
431  } else {
432  write_reals(datagram, &m00[0], length);
433  write_reals(datagram, &m01[0], length);
434  write_reals(datagram, &m02[0], length);
435  write_reals(datagram, &m10[0], length);
436  write_reals(datagram, &m11[0], length);
437  write_reals(datagram, &m12[0], length);
438  write_reals(datagram, &m20[0], length);
439  write_reals(datagram, &m21[0], length);
440  write_reals(datagram, &m22[0], length);
441  }
442  return;
443  }
444 #endif
445 
446  // First, convert the HPR's to quats. We expect quats to have better FFT
447  // consistency, and therefore compress better, even though they have an
448  // extra component.
449 
450  // However, because the quaternion will be normalized, we don't even have to
451  // write out all three components; any three can be used to determine the
452  // fourth (provided we ensure consistency of sign).
453 
454  vector_stdfloat qr, qi, qj, qk;
455 
456  qr.reserve(length);
457  qi.reserve(length);
458  qj.reserve(length);
459  qk.reserve(length);
460 
461  for (int i = 0; i < length; i++) {
462  LMatrix3 mat;
463  compose_matrix(mat, LVecBase3(1.0, 1.0, 1.0), LVecBase3(0.0, 0.0, 0.0),
464  array[i]);
465  if (_transpose_quats) {
466  mat.transpose_in_place();
467  }
468 
469  LOrientation rot(mat);
470  rot.normalize(); // This may not be necessary, but let's not take chances.
471 
472  if (rot.get_r() < 0) {
473  // Since rot == -rot, we can flip the quarternion if need be to keep the
474  // r component positive. This has two advantages. One, it makes it
475  // possible to infer r completely given i, j, and k (since we know it
476  // must be >= 0), and two, it helps protect against poor continuity
477  // caused by inadvertent flipping of the quarternion's sign between
478  // frames.
479 
480  // The choice of leaving r implicit rather than any of the other three
481  // seems to work the best in terms of guaranteeing continuity.
482  rot.set(-rot.get_r(), -rot.get_i(), -rot.get_j(), -rot.get_k());
483  }
484 
485 #ifdef NOTIFY_DEBUG
486  if (mathutil_cat.is_warning()) {
487  LMatrix3 mat2;
488  rot.extract_to_matrix(mat2);
489  if (!mat.almost_equal(mat2, 0.0001)) {
490  LVecBase3 hpr1, hpr2;
491  LVecBase3 scale, shear;
492  decompose_matrix(mat, scale, shear, hpr1);
493  decompose_matrix(mat2, scale, shear, hpr2);
494  mathutil_cat.warning()
495  << "Converted hpr to quaternion incorrectly!\n"
496  << " Source hpr: " << array[i] << ", or " << hpr1 << "\n";
497  mathutil_cat.warning(false)
498  << " Quaternion: " << rot << "\n"
499  << " Which represents: hpr " << hpr2 << " scale "
500  << scale << "\n";
501  }
502  }
503 #endif
504 
505  qr.push_back(rot.get_r());
506  qi.push_back(rot.get_i());
507  qj.push_back(rot.get_j());
508  qk.push_back(rot.get_k());
509  }
510 
511  // If quality is at least 102, we write all four quat components, instead of
512  // just the three. This is just for debugging.
513 #ifndef NDEBUG
514  if (_quality >= 102) {
515  if (length == 0) {
516  write_reals(datagram, nullptr, length);
517  } else {
518  write_reals(datagram, &qr[0], length);
519  }
520  }
521 #endif
522  if (length == 0) {
523  write_reals(datagram, nullptr, length);
524  write_reals(datagram, nullptr, length);
525  write_reals(datagram, nullptr, length);
526  } else {
527  write_reals(datagram, &qi[0], length);
528  write_reals(datagram, &qj[0], length);
529  write_reals(datagram, &qk[0], length);
530  }
531 }
532 
533 /**
534  * Reads the compression header that was written previously. This fills in
535  * the compression parameters necessary to correctly decompress the following
536  * data.
537  *
538  * Returns true if the header is read successfully, false otherwise.
539  */
541 read_header(DatagramIterator &di, int bam_minor_version) {
542  _bam_minor_version = bam_minor_version;
543  _quality = di.get_int8();
544 
545  if (mathutil_cat.is_debug()) {
546  mathutil_cat.debug()
547  << "Found compressed data at quality level " << _quality << "\n";
548  }
549 
550 #ifndef HAVE_FFTW
551  if (_quality <= 100) {
552  mathutil_cat.error()
553  << "FFTW library is not available; cannot read compressed data.\n";
554  return false;
555  }
556 #endif
557 
558  set_quality(_quality);
559 
560  if (_quality < 0) {
561  _fft_offset = di.get_float64();
562  _fft_factor = di.get_float64();
563  _fft_exponent = di.get_float64();
564  }
565 
566  return true;
567 }
568 
569 /**
570  * Reads an array of floating-point numbers. The result is pushed onto the
571  * end of the indicated vector, which is not cleared first; it is the user's
572  * responsibility to ensure that the array is initially empty. Returns true
573  * if the data is read correctly, false if there is an error.
574  */
576 read_reals(DatagramIterator &di, vector_stdfloat &array) {
577  int length = di.get_int32();
578 
579  if (_quality > 100) {
580  array.reserve(array.size() + length);
581 
582  // Special case: lossless output.
583  for (int i = 0; i < length; i++) {
584  array.push_back(di.get_stdfloat());
585  }
586  return true;
587  }
588 
589 #ifndef HAVE_FFTW
590  // If we don't have FFTW, we shouldn't get here.
591  return false;
592 
593 #else
594 
595  if (length == 0) {
596  // Special case: do nothing.
597  return true;
598  }
599 
600  if (length == 1) {
601  // Special case: just read in the one number.
602  array.push_back(di.get_stdfloat());
603  return true;
604  }
605 
606  // Normal case: read in the FFT array, and convert it back to (nearly) the
607  // original numbers.
608 
609  // First, check the reject_compression flag. If it's set, we decided to
610  // just write out the stream uncompressed.
611  bool reject_compression = di.get_bool();
612  if (reject_compression) {
613  array.reserve(array.size() + length);
614  for (int i = 0; i < length; i++) {
615  array.push_back(di.get_stdfloat());
616  }
617  return true;
618  }
619 
620  vector_double half_complex;
621  half_complex.reserve(length);
622  int num_read = 0;
623  while (num_read < length) {
624  num_read += read_run(di, half_complex);
625  }
626  nassertr(num_read == length, false);
627  nassertr((int)half_complex.size() == length, false);
628 
629  int fft_length = length / 2 + 1;
630  fftw_complex *fft_bins = (fftw_complex *)alloca(fft_length * sizeof(fftw_complex));
631 
632  int i;
633  for (i = 0; i < fft_length; i++) {
634  double scale_factor = get_scale_factor(i, fft_length);
635 
636  // For an explanation of this, see the compression code's comment about the
637  // halfcomplex format.
638 
639  fft_bins[i][0] = half_complex[i] * scale_factor;
640  if (i == 0) {
641  // First bin doesn't store imaginary component
642  fft_bins[i][1] = 0.0;
643  } else if ((i == fft_length - 1) && !(length & 1)) {
644  // Last bin doesn't store imaginary component with even lengths
645  fft_bins[i][1] = 0.0;
646  } else {
647  fft_bins[i][1] = half_complex[length - i] * scale_factor;
648  }
649  }
650 
651  // This is for an in-place transform. It doesn't violate strict aliasing
652  // rules because &fft_bins[0][0] is still a double pointer. This saves on
653  // precious stack space.
654  double *data = &fft_bins[0][0];
655 
656  // Note: This is an in-place DFT. `data` and `fft_bins` are aliases.
657  fftw_plan plan = get_real_decompress_plan(length);
658  fftw_execute_dft_c2r(plan, fft_bins, data);
659 
660  double scale = 1.0 / (double)length;
661  array.reserve(array.size() + length);
662  for (i = 0; i < length; i++) {
663  array.push_back(data[i] * scale);
664  }
665 
666  return true;
667 #endif
668 }
669 
670 /**
671  * Reads an array of HPR angles. The result is pushed onto the end of the
672  * indicated vector, which is not cleared first; it is the user's
673  * responsibility to ensure that the array is initially empty.
674  *
675  * new_hpr is a temporary, transitional parameter. If it is set false, the
676  * hprs are decompressed according to the old, broken hpr calculation; if
677  * true, the hprs are decompressed according to the new, correct hpr
678  * calculation.
679  */
681 read_hprs(DatagramIterator &di, pvector<LVecBase3> &array, bool new_hpr) {
682 #ifndef NDEBUG
683  if (_quality >= 104) {
684  // If quality level is at least 104, we don't even convert hpr to quat.
685  // This is just for debugging.
686  vector_stdfloat h, p, r;
687  bool okflag = true;
688  okflag =
689  read_reals(di, h) &&
690  read_reals(di, p) &&
691  read_reals(di, r);
692 
693  if (okflag) {
694  nassertr(h.size() == p.size() && p.size() == r.size(), false);
695  for (int i = 0; i < (int)h.size(); i++) {
696  array.push_back(LVecBase3(h[i], p[i], r[i]));
697  }
698  }
699 
700  return okflag;
701  }
702  if (_quality >= 103) {
703  // If quality level is 103, we read in a table of 3x3 rotation matrices.
704  // This is just for debugging.
705  vector_stdfloat
706  m00, m01, m02,
707  m10, m11, m12,
708  m20, m21, m22;
709  bool okflag = true;
710  okflag =
711  read_reals(di, m00) &&
712  read_reals(di, m01) &&
713  read_reals(di, m02) &&
714  read_reals(di, m10) &&
715  read_reals(di, m11) &&
716  read_reals(di, m12) &&
717  read_reals(di, m20) &&
718  read_reals(di, m21) &&
719  read_reals(di, m22);
720 
721  if (okflag) {
722  for (int i = 0; i < (int)m00.size(); i++) {
723  LMatrix3 mat(m00[i], m01[i], m02[i],
724  m10[i], m11[i], m12[i],
725  m20[i], m21[i], m22[i]);
726  LVecBase3 scale, shear, hpr;
727  if (new_hpr) {
728  decompose_matrix(mat, scale, shear, hpr);
729  } else {
730  decompose_matrix_old_hpr(mat, scale, shear, hpr);
731  }
732  array.push_back(hpr);
733  }
734  }
735 
736  return okflag;
737  }
738 #endif
739 
740  vector_stdfloat qr, qi, qj, qk;
741 
742  bool okflag = true;
743 
744 #ifndef NDEBUG
745  if (_quality >= 102) {
746  okflag = read_reals(di, qr);
747  }
748 #endif
749 
750  okflag =
751  okflag &&
752  read_reals(di, qi) &&
753  read_reals(di, qj) &&
754  read_reals(di, qk);
755 
756  if (okflag) {
757  nassertr(qi.size() == qj.size() && qj.size() == qk.size(), false);
758 
759  array.reserve(array.size() + qi.size());
760  for (int i = 0; i < (int)qi.size(); i++) {
761  LOrientation rot;
762 
763  // Infer the r component from the remaining three.
764  PN_stdfloat qr2 = 1.0 - (qi[i] * qi[i] + qj[i] * qj[i] + qk[i] * qk[i]);
765  PN_stdfloat qr1 = qr2 < 0.0 ? 0.0 : sqrtf(qr2);
766 
767  rot.set(qr1, qi[i], qj[i], qk[i]);
768 
769 #ifndef NDEBUG
770  if (_quality >= 102) {
771  // If we have written out all four components, use them.
772  rot[0] = qr[i];
773 
774  if (!IS_THRESHOLD_EQUAL(qr[i], qr1, 0.001)) {
775  mathutil_cat.warning()
776  << "qr[" << i << "] = " << qr[i] << ", qr1 = " << qr1
777  << ", diff is " << qr1 - qr[i] << "\n";
778  }
779  } else
780 #endif
781 
782  rot.normalize(); // Just for good measure.
783 
784  LMatrix3 mat;
785  rot.extract_to_matrix(mat);
786  if (_transpose_quats) {
787  mat.transpose_in_place();
788  }
789  LVecBase3 scale, shear, hpr;
790  if (new_hpr) {
791  decompose_matrix(mat, scale, shear, hpr);
792  } else {
793  decompose_matrix_old_hpr(mat, scale, shear, hpr);
794  }
795  array.push_back(hpr);
796  }
797  }
798 
799  return okflag;
800 }
801 
802 /**
803  * Reads an array of HPR angles. The result is pushed onto the end of the
804  * indicated vector, which is not cleared first; it is the user's
805  * responsibility to ensure that the array is initially empty.
806  */
809  return read_hprs(di, array, true);
810 }
811 
812 
813 /**
814  * Frees memory that has been allocated during past runs of the FFTCompressor.
815  * This is an optional call, but it may be made from time to time to empty the
816  * global cache that the compressor objects keep to facilitate fast
817  * compression/decompression.
818  */
820 free_storage() {
821 #ifdef HAVE_FFTW
822  RealPlans::iterator pi;
823  for (pi = _real_compress_plans.begin();
824  pi != _real_compress_plans.end();
825  ++pi) {
826  fftw_destroy_plan((*pi).second);
827  }
828  _real_compress_plans.clear();
829 
830  for (pi = _real_decompress_plans.begin();
831  pi != _real_decompress_plans.end();
832  ++pi) {
833  fftw_destroy_plan((*pi).second);
834  }
835  _real_decompress_plans.clear();
836 #endif
837 }
838 
839 /**
840  * Writes a sequence of integers that all require the same number of bits.
841  * Returns the number of integers written, i.e. run.size().
842  */
843 int FFTCompressor::
844 write_run(Datagram &datagram, FFTCompressor::RunWidth run_width,
845  const vector_double &run) {
846  if (run.empty()) {
847  return 0;
848  }
849  nassertr(run_width != RW_invalid, 0);
850 
851  if (run_width != RW_double) {
852  // If the width is anything other than RW_double, we write a single byte
853  // indicating the width and length of the upcoming run.
854 
855  if (run.size() <= RW_length_mask &&
856  ((int)run_width | run.size()) != RW_double) {
857  // If there are enough bits remaining in the byte, use them to indicate
858  // the length of the run. We have to be a little careful, however, not
859  // to accidentally write a byte that looks like an RW_double flag.
860  datagram.add_uint8((int)run_width | run.size());
861 
862  } else {
863  // Otherwise, write zero as the length, to indicate that we'll write the
864  // actual length in the following 16-bit word.
865  datagram.add_uint8(run_width);
866 
867  // Assuming, of course, that the length fits within 16 bits.
868  nassertr(run.size() < 65536, 0);
869  nassertr(run.size() != 0, 0);
870 
871  datagram.add_uint16(run.size());
872  }
873  }
874 
875  // Now write the data itself.
876  vector_double::const_iterator ri;
877  switch (run_width) {
878  case RW_0:
879  // If it's a string of zeroes, we're done!
880  break;
881 
882  case RW_8:
883  for (ri = run.begin(); ri != run.end(); ++ri) {
884  datagram.add_int8((int)*ri);
885  }
886  break;
887 
888  case RW_16:
889  for (ri = run.begin(); ri != run.end(); ++ri) {
890  datagram.add_int16((int)*ri);
891  }
892  break;
893 
894  case RW_32:
895  for (ri = run.begin(); ri != run.end(); ++ri) {
896  datagram.add_int32((int)*ri);
897  }
898  break;
899 
900  case RW_double:
901  for (ri = run.begin(); ri != run.end(); ++ri) {
902  // In the case of RW_double, we only write the numbers one at a time,
903  // each time preceded by the RW_double flag. Hopefully this will happen
904  // only rarely.
905  datagram.add_int8((int8_t)RW_double);
906  datagram.add_float64(*ri);
907  }
908  break;
909 
910  default:
911  break;
912  }
913 
914  return run.size();
915 }
916 
917 /**
918  * Reads a sequence of integers that all require the same number of bits.
919  * Returns the number of integers read. It is the responsibility of the user
920  * to clear the vector before calling this function, or the numbers read will
921  * be appended to the end.
922  */
923 int FFTCompressor::
924 read_run(DatagramIterator &di, vector_double &run) {
925  uint8_t start = di.get_uint8();
926  RunWidth run_width;
927  int length;
928 
929  if ((start & 0xff) == RW_double) {
930  // RW_double is a special case, and requires the whole byte. In this
931  // case, we don't encode a length, but assume it's only one.
932  run_width = RW_double;
933  length = 1;
934 
935  } else {
936  run_width = (RunWidth)(start & RW_width_mask);
937  length = start & RW_length_mask;
938  }
939 
940  if (length == 0) {
941  // If the length was zero, it means the actual length follows as a 16-bit
942  // word.
943  length = di.get_uint16();
944  }
945  nassertr(length != 0, 0);
946 
947  run.reserve(run.size() + length);
948 
949  int i;
950  switch (run_width) {
951  case RW_0:
952  for (i = 0; i < length; i++) {
953  run.push_back(0.0);
954  }
955  break;
956 
957  case RW_8:
958  for (i = 0; i < length; i++) {
959  run.push_back((double)(int)di.get_int8());
960  }
961  break;
962 
963  case RW_16:
964  for (i = 0; i < length; i++) {
965  run.push_back((double)(int)di.get_int16());
966  }
967  break;
968 
969  case RW_32:
970  for (i = 0; i < length; i++) {
971  run.push_back((double)(int)di.get_int32());
972  }
973  break;
974 
975  case RW_double:
976  for (i = 0; i < length; i++) {
977  run.push_back(di.get_float64());
978  }
979  break;
980 
981  default:
982  break;
983  }
984 
985  return length;
986 }
987 
988 /**
989  * Returns the appropriate scaling for the given bin in the FFT output.
990  *
991  * The scale factor is the value of one integer in the quantized data. As such,
992  * greater bins (higher, more noticeable frequencies) have *lower* scaling
993  * factors, which means greater precision.
994  */
995 double FFTCompressor::
996 get_scale_factor(int i, int length) const {
997  nassertr(i < length, 1.0);
998 
999  return _fft_offset +
1000  _fft_factor * pow((double)(length - i) / (double)(length), _fft_exponent);
1001 }
1002 
1003 /**
1004  * Returns a number between a and b, inclusive, according to the value of t
1005  * between 0 and 1, inclusive.
1006  */
1007 double FFTCompressor::
1008 interpolate(double t, double a, double b) {
1009  return a + t * (b - a);
1010 }
1011 
1012 /**
1013  * Returns a factor that indicates how erratically the values are changing.
1014  * The lower the result, the calmer the numbers, and the greater its
1015  * likelihood of being successfully compressed.
1016  */
1017 PN_stdfloat FFTCompressor::
1018 get_compressability(const PN_stdfloat *data, int length) const {
1019  // The result returned is actually the standard deviation of the table of
1020  // deltas between consecutive frames. This number is larger if the frames
1021  // have wildly different values.
1022 
1023  if (length <= 2) {
1024  return 0.0;
1025  }
1026 
1027  PN_stdfloat sum = 0.0;
1028  PN_stdfloat sum2 = 0.0;
1029  for (int i = 1; i < length; i++) {
1030  PN_stdfloat delta = data[i] - data[i - 1];
1031 
1032  sum += delta;
1033  sum2 += delta * delta;
1034  }
1035  PN_stdfloat variance = (sum2 - (sum * sum) / (length - 1)) / (length - 2);
1036  if (variance < 0.0) {
1037  // This can only happen due to tiny roundoff error.
1038  return 0.0;
1039  }
1040 
1041  PN_stdfloat std_deviation = csqrt(variance);
1042 
1043  return std_deviation;
1044 }
1045 
1046 
1047 
1048 #ifdef HAVE_FFTW
1049 
1050 /**
1051  * Returns a FFTW plan suitable for compressing a float array of the indicated
1052  * length.
1053  */
1054 static fftw_plan
1055 get_real_compress_plan(int length) {
1056  RealPlans::iterator pi;
1057  pi = _real_compress_plans.find(length);
1058  if (pi != _real_compress_plans.end()) {
1059  return (*pi).second;
1060  }
1061 
1062  fftw_plan plan;
1063  plan = fftw_plan_dft_r2c_1d(length, nullptr, nullptr, FFTW_ESTIMATE);
1064  _real_compress_plans.insert(RealPlans::value_type(length, plan));
1065 
1066  return plan;
1067 }
1068 
1069 /**
1070  * Returns a FFTW plan suitable for decompressing a float array of the
1071  * indicated length.
1072  */
1073 static fftw_plan
1074 get_real_decompress_plan(int length) {
1075  RealPlans::iterator pi;
1076  pi = _real_decompress_plans.find(length);
1077  if (pi != _real_decompress_plans.end()) {
1078  return (*pi).second;
1079  }
1080 
1081  fftw_plan plan;
1082  plan = fftw_plan_dft_c2r_1d(length, nullptr, nullptr, FFTW_ESTIMATE);
1083  _real_decompress_plans.insert(RealPlans::value_type(length, plan));
1084 
1085  return plan;
1086 }
1087 
1088 #endif
A class to retrieve the individual data elements previously stored in a Datagram.
uint8_t get_uint8()
Extracts an unsigned 8-bit integer.
int16_t get_int16()
Extracts a signed 16-bit integer.
PN_stdfloat get_stdfloat()
Extracts either a 32-bit or a 64-bit floating-point number, according to Datagram::set_stdfloat_doubl...
uint16_t get_uint16()
Extracts an unsigned 16-bit integer.
int8_t get_int8()
Extracts a signed 8-bit integer.
PN_float64 get_float64()
Extracts a 64-bit floating-point number.
bool get_bool()
Extracts a boolean value.
int32_t get_int32()
Extracts a signed 32-bit integer.
An ordered list of data elements, formatted in memory for transmission over a socket or writing to a ...
Definition: datagram.h:38
void add_int16(int16_t value)
Adds a signed 16-bit integer to the datagram.
Definition: datagram.I:58
void add_int32(int32_t value)
Adds a signed 32-bit integer to the datagram.
Definition: datagram.I:67
void add_uint8(uint8_t value)
Adds an unsigned 8-bit integer to the datagram.
Definition: datagram.I:50
void add_stdfloat(PN_stdfloat value)
Adds either a 32-bit or a 64-bit floating-point number, according to set_stdfloat_double().
Definition: datagram.I:133
void add_bool(bool value)
Adds a boolean value to the datagram.
Definition: datagram.I:34
void add_uint16(uint16_t value)
Adds an unsigned 16-bit integer to the datagram.
Definition: datagram.I:85
void add_float64(PN_float64 value)
Adds a 64-bit floating-point number to the datagram.
Definition: datagram.I:123
void add_int8(int8_t value)
Adds a signed 8-bit integer to the datagram.
Definition: datagram.I:42
static void free_storage()
Frees memory that has been allocated during past runs of the FFTCompressor.
bool read_reals(DatagramIterator &di, vector_stdfloat &array)
Reads an array of floating-point numbers.
void set_use_error_threshold(bool use_error_threshold)
Enables or disables the use of the error threshold measurement to put a cap on the amount of damage d...
void set_transpose_quats(bool flag)
Sets the transpose_quats flag.
FFTCompressor()
Constructs a new compressor object with default parameters.
bool read_hprs(DatagramIterator &di, pvector< LVecBase3 > &array, bool new_hpr)
Reads an array of HPR angles.
bool get_use_error_threshold() const
Returns whether the error threshold measurement is enabled.
void set_quality(int quality)
Sets the quality factor for the compression.
bool read_header(DatagramIterator &di, int bam_minor_version)
Reads the compression header that was written previously.
int get_quality() const
Returns the quality number that was previously set via set_quality().
static bool is_compression_available()
Returns true if the FFTW library is compiled in, so that this class is actually capable of doing usef...
void write_header(Datagram &datagram)
Writes the compression parameters to the indicated datagram.
void write_reals(Datagram &datagram, const PN_stdfloat *array, int length)
Writes an array of floating-point numbers to the indicated datagram.
void write_hprs(Datagram &datagram, const LVecBase3 *array, int length)
Writes an array of HPR angles to the indicated datagram.
bool get_transpose_quats() const
Returns the transpose_quats flag.
This is our own Panda specialization on the default STL map.
Definition: pmap.h:49
This is our own Panda specialization on the default STL vector.
Definition: pvector.h:42
PANDA 3D SOFTWARE Copyright (c) Carnegie Mellon University.
PANDA 3D SOFTWARE Copyright (c) Carnegie Mellon University.
PANDA 3D SOFTWARE Copyright (c) Carnegie Mellon University.
PANDA 3D SOFTWARE Copyright (c) Carnegie Mellon University.
PANDA 3D SOFTWARE Copyright (c) Carnegie Mellon University.
PANDA 3D SOFTWARE Copyright (c) Carnegie Mellon University.
PANDA 3D SOFTWARE Copyright (c) Carnegie Mellon University.