Panda3D
Loading...
Searching...
No Matches
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.
36static fftw_plan get_real_compress_plan(int length);
37static fftw_plan get_real_decompress_plan(int length);
38
39typedef pmap<int, fftw_plan> RealPlans;
40static RealPlans _real_compress_plans;
41static RealPlans _real_decompress_plans;
42
43#endif
44
45/**
46 * Constructs a new compressor object with default parameters.
47 */
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 */
89set_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 */
139get_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 */
152set_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 */
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 */
171set_transpose_quats(bool flag) {
172 _transpose_quats = flag;
173}
174
175/**
176 * Returns the transpose_quats flag. See set_transpose_quats().
177 */
179get_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 */
189write_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 */
202write_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 */
371write_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 */
541read_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 */
576read_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 */
681read_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 */
820free_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 */
843int FFTCompressor::
844write_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 */
923int FFTCompressor::
924read_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 */
995double FFTCompressor::
996get_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 */
1007double FFTCompressor::
1008interpolate(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 */
1017PN_stdfloat FFTCompressor::
1018get_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 */
1054static fftw_plan
1055get_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 */
1073static fftw_plan
1074get_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.