1 /*
2 * Copyright (c) 2005-2008, 2010 Genome Research Ltd.
3 * Author(s): James Bonfield
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above
12 * copyright notice, this list of conditions and the following
13 * disclaimer in the documentation and/or other materials provided
14 * with the distribution.
15 *
16 * 3. Neither the names Genome Research Ltd and Wellcome Trust Sanger
17 * Institute nor the names of its contributors may be used to endorse
18 * or promote products derived from this software without specific
19 * prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY GENOME RESEARCH LTD AND CONTRIBUTORS "AS
22 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
23 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
24 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL GENOME RESEARCH
25 * LTD OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 /*
35 * Author(s): James Bonfield
36 *
37 * Copyright (c) 2001 MEDICAL RESEARCH COUNCIL
38 * All rights reserved
39 *
40 * Redistribution and use in source and binary forms, with or without
41 * modification, are permitted provided that the following conditions are met:
42 *
43 * 1 Redistributions of source code must retain the above copyright notice,
44 * this list of conditions and the following disclaimer.
45 *
46 * 2 Redistributions in binary form must reproduce the above copyright
47 * notice, this list of conditions and the following disclaimer in
48 * the documentation and/or other materials provided with the
49 * distribution.
50 *
51 * 3 Neither the name of the MEDICAL RESEARCH COUNCIL, THE LABORATORY OF
52 * MOLECULAR BIOLOGY nor the names of its contributors may be used
53 * to endorse or promote products derived from this software without
54 * specific prior written permission.
55 *
56 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
57 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
58 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
59 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
60 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
61 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
62 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
63 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
64 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
65 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
66 * POSSIBILITY OF SUCH DAMAGE.
67 */
68
69 #ifdef HAVE_CONFIG_H
70 #include "io_lib_config.h"
71 #endif
72
73 #include <stdio.h>
74 #include <unistd.h>
75 #include <stdlib.h>
76 #include <string.h>
77 #include <zlib.h>
78 #include <assert.h>
79 #include <math.h>
80 #include <ctype.h>
81
82 #ifndef M_PI
83 # define M_PI 3.14159265358979323846
84 #endif
85
86 #include "io_lib/ztr.h"
87 #include "io_lib/os.h"
88 #include "io_lib/compression.h"
89 #include "io_lib/xalloc.h"
90
91 #ifndef NDEBUG
92 # define NDEBUG
93 #endif
94
95 /*
96 * ---------------------------------------------------------------------------
97 * ZTR_FORM_ZLIB
98 * ---------------------------------------------------------------------------
99 */
100
101 /*
102 * Some comments on zlib usage.
103 *
104 * - Ideally for trace data, after decorrelation, we should use Z_FILTERED.
105 * Empirical studies show that this gives the best compression ratio, but
106 * it is slow (about the same speed as normal gzip). MUCH faster is huffman
107 * only, and it doesn't give radically different compression ratios.
108 *
109 * - When compressing using Z_HUFFMAN_ONLY we used compression level
110 * '1' as this invokes the deflate_fast() algorithm. It makes no
111 * difference to the compression level, but it seems to be quicker still.
112 *
113 */
114
115 /*
116 * zlib_huff()
117 *
118 * Compresses data using huffman encoding, as implemented by zlib.
119 *
120 * Arguments:
121 * uncomp Uncompressed input data
122 * uncomp_len Length of uncomp data
123 * comp_len Output: length of compressed data
124 *
125 * Returns:
126 * Compressed data if successful
127 * NULL if not successful
128 */
zlib_huff(char * uncomp,int uncomp_len,int strategy,int * comp_len)129 char *zlib_huff(char *uncomp, int uncomp_len, int strategy, int *comp_len) {
130 z_stream zstr;
131 int err;
132 int comp_len_tmp = (int)(uncomp_len * 1.001 + 12); /* Maximum expansion size */
133 char *comp = (char *)xmalloc(comp_len_tmp+5);
134 int c_len;
135
136 /* Initialise zlib */
137 zstr.zalloc = (alloc_func)0;
138 zstr.zfree = (free_func)0;
139 zstr.opaque = (voidpf)0;
140
141 if ((err = deflateInit2(&zstr, 1, Z_DEFLATED, 15, 8, strategy)) != Z_OK) {
142 fprintf(stderr, "zlib errror in deflateInit2(): %d\n", err);
143 return NULL;
144 }
145
146 /* Set up input and output buffers */
147 zstr.next_in = (unsigned char *)uncomp;
148 zstr.avail_in = uncomp_len;
149 zstr.next_out = (unsigned char *)comp+5;
150 zstr.avail_out = comp_len_tmp;
151
152 /* Do the compression */
153 if ((err = deflate(&zstr, Z_FINISH)) != Z_STREAM_END) {
154 fprintf(stderr, "zlib errror in deflate(): %d\n", err);
155 return NULL;
156 }
157
158 /* Tidy up */
159 deflateEnd(&zstr);
160
161 c_len = zstr.total_out;
162
163 /* Return */
164 comp[0] = ZTR_FORM_ZLIB;
165 comp[1] = (uncomp_len >> 0) & 0xff;
166 comp[2] = (uncomp_len >> 8) & 0xff;
167 comp[3] = (uncomp_len >> 16) & 0xff;
168 comp[4] = (uncomp_len >> 24) & 0xff;
169
170 if (comp_len)
171 *comp_len = c_len+5;
172
173 return comp;
174 }
175
176 /*
177 * zlib_dehuff()
178 *
179 * Uncompresses data using huffman encoding, as implemented by zlib.
180 *
181 * Arguments:
182 * comp Compressed input data
183 * comp_len Length of comp data
184 * uncomp_len Output: length of uncompressed data
185 *
186 * Returns:
187 * Uncompressed data if successful
188 * NULL if not successful
189 */
zlib_dehuff(char * comp,int comp_len,int * uncomp_len)190 char *zlib_dehuff(char *comp, int comp_len, int *uncomp_len) {
191 z_stream zstr;
192 int err;
193 char *uncomp;
194 int ulen;
195
196 /* Allocate */
197 ulen =
198 ((unsigned char)comp[1] << 0) +
199 ((unsigned char)comp[2] << 8) +
200 ((unsigned char)comp[3] << 16) +
201 ((unsigned char)comp[4] << 24);
202 uncomp = (char *)xmalloc(ulen);
203
204 /* Initialise zlib */
205 zstr.zalloc = (alloc_func)0;
206 zstr.zfree = (free_func)0;
207 zstr.opaque = (voidpf)0;
208
209 if ((err = inflateInit(&zstr)) != Z_OK) {
210 fprintf(stderr, "zlib errror in inflateInit(): %d\n", err);
211 return NULL;
212 }
213
214 /* Set up input and output buffers */
215 zstr.next_in = (unsigned char *)comp+5;
216 zstr.avail_in = comp_len-5;
217 zstr.next_out = (unsigned char *)uncomp;
218 zstr.avail_out = ulen;
219
220 /* Do the decompression */
221 if ((err = inflate(&zstr, Z_FINISH)) != Z_STREAM_END) {
222 fprintf(stderr, "zlib errror in deflate(): %d\n", err);
223 return NULL;
224 }
225
226 /* Tidy up */
227 inflateEnd(&zstr);
228
229 if (uncomp_len)
230 *uncomp_len = ulen;
231
232 return uncomp;
233 }
234
235
236 /*
237 * ---------------------------------------------------------------------------
238 * ZTR_FORM_RLE
239 * ---------------------------------------------------------------------------
240 */
241
242 /*
243 * Run length encoding.
244 *
245 * Any run of 3 or more identical characters (up to 255 in a row) are replaced
246 * by a 'guard' byte followed by the number of characters followed by
247 * the character value itself.
248 * Any single guard value in the input is escaped using 'guard 0'.
249 *
250 * Specifying guard as -1 will automatically pick one of the least used
251 * characters in the input as the guard.
252 *
253 * Arguments:
254 * uncomp Input data
255 * uncomp_len Length of input data 'uncomp'
256 * guard Guard byte - used to encode "N" copies of data
257 * comp_len Output: length of compressed data
258 *
259 * Returns:
260 * Compressed data if successful
261 * NULL if not successful
262 */
rle(char * uncomp,int uncomp_len,int guard,int * comp_len)263 char *rle(char *uncomp, int uncomp_len, int guard, int *comp_len) {
264 int i, k, c_len = 0;
265 char *comp = xmalloc(2 * uncomp_len + 6);
266 char *out = comp + 6;
267
268 /* A guard of -1 implies to search for the least frequent symbol */
269 if (guard == -1) {
270 int cnt[256];
271 int bestcnt = uncomp_len + 1;
272
273 for (i = 0; i < 256; i++)
274 cnt[i] = 0;
275
276 for (i = 0; i < uncomp_len; i++) {
277 cnt[(unsigned char)uncomp[i]]++;
278 }
279
280 for (i = 0; i < 256; i++) {
281 if (cnt[i] < bestcnt) {
282 bestcnt = cnt[i];
283 guard = i;
284 }
285 }
286 }
287
288 for (i = 0; i < uncomp_len; i=k) {
289 /*
290 * Detect blocks of up identical bytes up to 255 bytes long.
291 */
292 for (k = i; k < uncomp_len && uncomp[i] == uncomp[k]; k++)
293 if (k-i == 255)
294 break;
295
296 /* 1, 2 or 3 bytes are best stored "as is" */
297 if (k-i < 4) {
298 do {
299 /*
300 * If we find 'guard' in our sequence, escape it by
301 * outputting 'guard' <zero>. (We know that we'll never
302 * write out zero copies of a token in our rle compression
303 * algorithm.)
304 */
305 if ((unsigned char)(uncomp[i]) == guard) {
306 out[c_len++] = guard;
307 out[c_len++] = 0;
308 } else {
309 out[c_len++] = uncomp[i];
310 }
311 i++;
312 } while (k >= i+1);
313 } else {
314 /* More than 3 bytes: store as ('guard', length, byte value) */
315 out[c_len++] = guard;
316 out[c_len++] = k-i;
317 out[c_len++] = uncomp[i];
318 }
319 }
320
321 /* Return */
322 comp[0] = ZTR_FORM_RLE;
323 comp[1] = (uncomp_len >> 0) & 0xff;
324 comp[2] = (uncomp_len >> 8) & 0xff;
325 comp[3] = (uncomp_len >> 16) & 0xff;
326 comp[4] = (uncomp_len >> 24) & 0xff;
327 comp[5] = guard;
328
329 if (comp_len)
330 *comp_len = c_len+6;
331
332 return comp;
333 }
334
335 /*
336 * Reverses run length encoding.
337 *
338 * Arguments:
339 * comp Compressed input data
340 * comp_len Length of comp data
341 * uncomp_len Output: length of uncompressed data
342 *
343 * Returns:
344 * Uncompressed data if successful
345 * NULL if not successful
346 */
347 /* ARGSUSED */
unrle(char * comp,int comp_len,int * uncomp_len)348 char *unrle(char *comp, int comp_len, int *uncomp_len) {
349 int in_i, out_i, i, val, count, out_len;
350 char *uncomp;
351 unsigned char *in = (unsigned char *)comp+6;
352 char *out;
353 int guard = (unsigned char)comp[5];
354
355 /* Allocate */
356 out_len =
357 ((unsigned char)comp[1] << 0) +
358 ((unsigned char)comp[2] << 8) +
359 ((unsigned char)comp[3] << 16) +
360 ((unsigned char)comp[4] << 24);
361 out = uncomp = (char *)xmalloc(out_len);
362
363 for (in_i = out_i = 0; out_i < out_len; in_i++) {
364 if (in[in_i] != guard) {
365 /* When not 'guard' it's easy - just output this token */
366 assert(out_i >= 0 && out_i < out_len);
367 out[out_i++] = in[in_i];
368
369 } else {
370 /*
371 * Found an 'guard' token. If next token is zero, then
372 * we were simply escaping a real 'guard' token in the input
373 * data, otherwise output a string of bytes.
374 */
375 count = in[++in_i];
376 if (count != 0) {
377 val = in[++in_i];
378 for (i = 0; i < count; i++) {
379 assert(out_i >= 0 && out_i < out_len);
380 out[out_i++] = val;
381 }
382 } else {
383 assert(out_i >= 0 && out_i < out_len);
384 out[out_i++] = guard;
385 }
386 }
387 }
388
389 if (uncomp_len)
390 *uncomp_len = out_len;
391
392 return uncomp;
393 }
394
395 /*
396 * ---------------------------------------------------------------------------
397 * ZTR_FORM_XRLE
398 * ---------------------------------------------------------------------------
399 */
400
401 /*
402 * Mutli-byte run length encoding.
403 *
404 * Any run of 3 or more identical characters (up to 255 in a row) are replaced
405 * by a 'guard' byte followed by the number of characters followed by
406 * the character value itself.
407 * Any single guard value in the input is escaped using 'guard 0'.
408 *
409 * Specifying guard as -1 will automatically pick one of the least used
410 * characters in the input as the guard.
411 *
412 * Arguments:
413 * uncomp Input data
414 * uncomp_len Length of input data 'uncomp'
415 * guard Guard byte - used to encode "N" copies of data
416 * rsz Size of blocks to compare for run checking.
417 * comp_len Output: length of compressed data
418 *
419 * Returns:
420 * Compressed data if successful
421 * NULL if not successful
422 */
xrle(char * uncomp,int uncomp_len,int guard,int rsz,int * comp_len)423 char *xrle(char *uncomp, int uncomp_len, int guard, int rsz, int *comp_len) {
424 char *comp = (char *)malloc(2 * uncomp_len + 3);
425 char *out = comp;
426 int i, j, k;
427
428 /* A guard of -1 implies to search for the least frequent symbol */
429 if (guard == -1) {
430 int cnt[256];
431 int bestcnt = uncomp_len + 1;
432
433 for (i = 0; i < 256; i++)
434 cnt[i] = 0;
435
436 for (i = 0; i < uncomp_len; i++) {
437 cnt[(unsigned char)uncomp[i]]++;
438 }
439
440 for (i = 0; i < 256; i++) {
441 if (cnt[i] < bestcnt) {
442 bestcnt = cnt[i];
443 guard = i;
444 }
445 }
446 }
447
448 *out++ = ZTR_FORM_XRLE;
449 *out++ = rsz;
450 *out++ = guard;
451
452 for (i = 0; i < uncomp_len; i = k) {
453 /* Count repeats */
454 k = i + rsz;
455 while (k <= uncomp_len - rsz && !memcmp(&uncomp[i], &uncomp[k], rsz)) {
456 k += rsz;
457 if ((k-i)/rsz == 255)
458 break;
459 }
460
461 if (k-i > rsz) {
462 /* Duplicates, so RLE */
463 *out++ = guard;
464 *out++ = (k-i)/rsz;
465 for (j = 0; j < rsz; j++)
466 *out++ = uncomp[i+j];
467 } else {
468 /* No dups, store as is escaping guarding as appropriate */
469 if ((unsigned char)(uncomp[i]) == guard) {
470 *out++ = guard;
471 *out++ = 0;
472 } else {
473 *out++ = uncomp[i];
474 }
475 k = i+1;
476 }
477 }
478
479 *comp_len = out-comp;
480
481 return comp;
482 }
483
484
485 /*
486 * Reverses multi-byte run length encoding.
487 *
488 * Arguments:
489 * comp Compressed input data
490 * comp_len Length of comp data
491 * uncomp_len Output: length of uncompressed data
492 *
493 * Returns:
494 * Uncompressed data if successful
495 * NULL if not successful
496 */
unxrle(char * comp,int comp_len,int * uncomp_len)497 char *unxrle(char *comp, int comp_len, int *uncomp_len) {
498 char *uncomp;
499 char *out;
500 int rsz = (unsigned char)comp[1];
501 int guard = (unsigned char)comp[2];
502 unsigned char *in;
503 int unclen = 0, cpos, len;
504
505 /* Calculate uncompressed length */
506 for (in = (unsigned char *)comp, cpos = 3; cpos < comp_len; unclen++) {
507 if (in[cpos++] == guard) {
508 if ((len = in[cpos++])) {
509 cpos += rsz;
510 unclen += len*rsz -1;
511 }
512 }
513 }
514 *uncomp_len = unclen;
515
516 /* Expand */
517 uncomp = out = (char *)malloc(unclen+1);
518 for (in = (unsigned char *)comp, cpos = 3; cpos < comp_len;) {
519 char c;
520 if ((c = in[cpos++]) != guard) {
521 *out++ = c;
522 } else {
523 int len = in[cpos++];
524 if (len) {
525 int i, j;
526 for (i = 0; i < len; i++) {
527 for (j = 0; j < rsz; j++) {
528 *out++ = in[cpos+j];
529 }
530 }
531 cpos += rsz;
532 } else {
533 *out++ = guard;
534 }
535 }
536 }
537 *out++ = 0;
538
539 return uncomp;
540 }
541
542 /*
543 * Mutli-byte run length encoding.
544 *
545 * Steps along in words of size 'rsz'. Unlike XRLE above this does run-length
546 * encoding by writing out an additional "length" word every time 2 or more
547 * words in a row are spotted. This removes the need for a guard byte.
548 *
549 * Additionally this method ensures that both input and output formats remain
550 * aligned on words of size 'rsz'.
551 *
552 * Arguments:
553 * uncomp Input data
554 * uncomp_len Length of input data 'uncomp'
555 * rsz Size of blocks to compare for run checking.
556 * comp_len Output: length of compressed data
557 *
558 * Returns:
559 * Compressed data if successful
560 * NULL if not successful
561 */
xrle2(char * uncomp,int uncomp_len,int rsz,int * comp_len)562 char *xrle2(char *uncomp, int uncomp_len, int rsz, int *comp_len) {
563 char *comp = (char *)malloc(1.4*uncomp_len + rsz);
564 char *last = uncomp;
565 char *out = comp;
566 int i, j, k, run_len = 0;
567
568 *out++ = ZTR_FORM_XRLE2;
569 *out++ = rsz;
570 for (i = 2; i < rsz; i++)
571 *out++ = -40;
572
573 /* FIXME: how to deal with uncomp_len not being a multiple of rsz */
574 for (i = 0; i < uncomp_len; i += rsz) {
575 /* FIXME: use inline #def versions of memcmp/memcpy for speed? */
576 memcpy(out, &uncomp[i], rsz);
577 out += rsz;
578
579 if (memcmp(last, &uncomp[i], rsz) == 0) {
580 run_len++;
581 } else {
582 run_len = 1;
583 last = &uncomp[i];
584 }
585
586 /* NB: >= 3 is more optimal in many cases */
587 if (run_len >= 2) {
588 /* Count remaining copies */
589 for (k = i+rsz; k < uncomp_len && run_len < 257; k += rsz) {
590 if (memcmp(last, &uncomp[k], rsz) != 0)
591 break;
592 run_len++;
593 }
594 run_len -= 2;
595
596 *out++ = run_len;
597 for (j = 1; j < rsz; j++) {
598 /* Padding with last reduces entropy compared to padding
599 * with copies of run_len
600 */
601 *out++ = last[j];
602 }
603 i = k-rsz;
604
605 last = out-rsz;
606 run_len = 0;
607 }
608 }
609
610 *comp_len = out-comp;
611
612 return comp;
613 }
614
615 /*
616 * Reverses multi-byte run length encoding (xrle_new).
617 *
618 * Arguments:
619 * comp Compressed input data
620 * comp_len Length of comp data
621 * uncomp_len Output: length of uncompressed data
622 *
623 * Returns:
624 * Uncompressed data if successful
625 * NULL if not successful
626 */
unxrle2(char * comp,int comp_len,int * uncomp_len)627 char *unxrle2(char *comp, int comp_len, int *uncomp_len) {
628 char *out, *last;
629 int out_len, out_alloc, rsz, i, j, run_len;
630
631 out_alloc = comp_len*2; /* just an estimate */
632 out_len = 0;
633 if (NULL == (out = (char *)malloc(out_alloc)))
634 return NULL;
635
636 if (*comp++ != ZTR_FORM_XRLE2)
637 return NULL;
638
639 /* Read rsz and swallow padding */
640 rsz = *comp++;
641 comp_len -= 2;
642 for (i = 2; i < rsz; i++) {
643 comp++;
644 comp_len--;
645 }
646
647 /* Uncompress */
648 run_len = 0;
649 last = comp;
650 for (i = 0; i < comp_len;) {
651 while (out_len + rsz > out_alloc) {
652 out_alloc *= 2;
653 if (NULL == (out = (char *)realloc(out, out_alloc)))
654 return NULL;
655 }
656 memcpy(&out[out_len], &comp[i], rsz);
657
658 if (memcmp(&out[out_len], last, rsz) == 0) {
659 run_len++;
660 } else {
661 run_len = 1;
662 }
663
664 i += rsz;
665 out_len += rsz;
666
667 /* NB: >= 3 is more optimal in many cases */
668 if (run_len >= 2) {
669 /* Count remaining copies */
670 run_len = (unsigned char)comp[i];
671 i += rsz;
672
673 while (out_len + run_len * rsz > out_alloc) {
674 out_alloc *= 2;
675 if (NULL == (out = (char *)realloc(out, out_alloc)))
676 return NULL;
677 }
678
679 for (j = 0; j < run_len; j++) {
680 memcpy(&out[out_len], last, rsz);
681 out_len += rsz;
682 }
683 run_len = 0;
684 }
685
686 last = &comp[i-rsz];
687 }
688
689 /* Shrink back down to avoid excessive memory usage */
690 out = realloc(out, out_len);
691 *uncomp_len = out_len;
692
693 return out;
694 }
695
696
697 /*
698 * ---------------------------------------------------------------------------
699 * ZTR_FORM_DELTA1
700 * ---------------------------------------------------------------------------
701 */
702
703 /*
704 * This replaces 'samples' with successive differences between samples.
705 * These implementations support 'level's of 1, 2 and 3.
706 *
707 * NB: This is analogous to our SCF delta_samples1 (etc) function, except that
708 * this function about 40% faster.
709 *
710 * Implementation ideas taken from Jean Thierry-Mieg's CTF code.
711 */
712
713 /*
714 * decorrelate1()
715 *
716 * Produce successive deltas from a 1-byte array.
717 *
718 * Arguments:
719 * uncomp Uncompressed data
720 * uncomp_len Length of uncompressed data
721 * level Differencing level (must be 1, 2 or 3)
722 * comp_len Return: where to store new compressed length
723 *
724 * Returns:
725 * Success: A decorrelated buffer (malloced)
726 * Failure: NULL
727 */
decorrelate1(char * x_uncomp,int uncomp_len,int level,int * comp_len)728 char *decorrelate1(char *x_uncomp,
729 int uncomp_len,
730 int level,
731 int *comp_len) {
732 int i, z;
733 int u1 = 0, u2 = 0, u3 = 0;
734 char *comp = (char *)xmalloc(uncomp_len + 2);
735 unsigned char *u_uncomp = (unsigned char *)x_uncomp;
736
737 if (!comp)
738 return NULL;
739
740 comp+=2;
741 switch (level) {
742 case 1:
743 for (i = 0; i < uncomp_len; i++) {
744 z = u1;
745 u1 = u_uncomp[i];
746 comp[i] = u_uncomp[i] - z;
747 }
748 break;
749
750 case 2:
751 for (i = 0; i < uncomp_len; i++) {
752 z = 2*u1 - u2;
753 u2 = u1;
754 u1 = u_uncomp[i];
755 comp[i] = u_uncomp[i] - z;
756 }
757 break;
758
759 case 3:
760 for (i = 0; i < uncomp_len; i++) {
761 z = 3*u1 - 3*u2 + u3;
762 u3 = u2;
763 u2 = u1;
764 u1 = u_uncomp[i];
765 comp[i] = u_uncomp[i] - z;
766 }
767 break;
768
769 default:
770 return NULL;
771 }
772 comp-=2;
773 comp[0] = ZTR_FORM_DELTA1;
774 comp[1] = level;
775
776 *comp_len = uncomp_len+2;
777
778 return comp;
779 }
780
781 #define ABS(a) ((a) > 0 ? (a) : -(a))
782
783 /* ZTR_FORM_DDELTA1 - experimental */
decorrelate1dyn(char * x_uncomp,int uncomp_len,int * comp_len)784 char *decorrelate1dyn(char *x_uncomp,
785 int uncomp_len,
786 int *comp_len) {
787 int i, j, z[4];
788 int u1 = 0, u2 = 0, u3 = 0;
789 char *comp = (char *)xmalloc(uncomp_len + 2);
790 unsigned char *u_uncomp = (unsigned char *)x_uncomp;
791 int level = 3; /* default level */
792 int last_level = level;
793 int best;
794
795 if (!comp)
796 return NULL;
797
798 comp+=2;
799 for (i = 0; i < uncomp_len; i++) {
800 z[1] = u1;
801 z[2] = 2*u1 - u2;
802 z[3] = 3*u1 - 3*u2 + u3;
803 comp[i] = u_uncomp[i] - z[last_level];
804 best = 10000;
805 for (j = 1; j < 3; j++) {
806 if (ABS(u_uncomp[i] - z[j]) < best) {
807 best = ABS(u_uncomp[i] - z[j]);
808 last_level = j;
809 }
810 }
811 u3 = u2;
812 u2 = u1;
813 u1 = u_uncomp[i];
814 }
815 comp-=2;
816 comp[0] = ZTR_FORM_DDELTA1;
817 comp[1] = level;
818
819 *comp_len = uncomp_len+2;
820
821 return comp;
822 }
823
824 /*
825 * recorrelate1()
826 *
827 * The reverse of decorrelate1()
828 *
829 * Arguments:
830 * comp Compressed input data
831 * comp_len Length of comp data
832 * uncomp_len Output: length of uncompressed data
833 *
834 * Returns:
835 * Success: uncompressed data
836 * Failure: NULL
837 */
recorrelate1(char * x_comp,int comp_len,int * uncomp_len)838 char *recorrelate1(char *x_comp,
839 int comp_len,
840 int *uncomp_len) {
841 int i, z;
842 int u1 = 0, u2 = 0, u3 = 0;
843 int level = x_comp[1];
844 char *uncomp;
845
846 uncomp = (char *)xmalloc(comp_len-2);
847 if (!uncomp)
848 return NULL;
849
850 x_comp+=2;
851 comp_len-=2;
852 *uncomp_len = comp_len;
853
854 switch (level) {
855 case 1:
856 for (i = 0; i < comp_len; i++) {
857 z = u1;
858 u1 = uncomp[i] = x_comp[i] + z;
859 }
860 break;
861
862 case 2:
863 for (i = 0; i < comp_len; i++) {
864 z = 2*u1 - u2;
865 u2 = u1;
866 u1 = uncomp[i] = x_comp[i] + z;
867 }
868 break;
869
870 case 3:
871 for (i = 0; i < comp_len; i++) {
872 z = 3*u1 - 3*u2 + u3;
873 u3 = u2;
874 u2 = u1;
875 u1 = uncomp[i] = x_comp[i] + z;
876 }
877 break;
878 }
879
880 return uncomp;
881 }
882
883 /*
884 * ---------------------------------------------------------------------------
885 * ZTR_FORM_DELTA2
886 * ---------------------------------------------------------------------------
887 */
888
889 /*
890 * decorrelate2()
891 *
892 * Produce successive deltas from a 2-byte array (big endian)
893 *
894 * Arguments:
895 * uncomp Uncompressed data
896 * uncomp_len Length of uncompressed data
897 * level Differencing level (must be 1, 2 or 3)
898 * comp_len Return: where to store new compressed length
899 *
900 * Returns:
901 * Success: A decorrelated buffer (malloced)
902 * Failure: NULL
903 */
decorrelate2(char * x_uncomp,int uncomp_len,int level,int * comp_len)904 char *decorrelate2(char *x_uncomp,
905 int uncomp_len,
906 int level,
907 int *comp_len) {
908 int i, z, delta;
909 int u1 = 0, u2 = 0, u3 = 0;
910 char *comp = (char *)xmalloc(uncomp_len + 2);
911 unsigned char *u_uncomp = (unsigned char *)x_uncomp;
912
913 if (!comp)
914 return NULL;
915
916 comp+=2;
917 switch (level) {
918 case 1:
919 for (i = 0; i < uncomp_len; i+=2) {
920 z = u1;
921 u1 = (u_uncomp[i] << 8) + u_uncomp[i+1];
922 delta = u1 - z;
923 comp[i ] = (delta >> 8) & 0xff;
924 comp[i+1] = (delta >> 0) & 0xff;
925 }
926 break;
927
928 case 2:
929 for (i = 0; i < uncomp_len; i+=2) {
930 z = 2*u1 - u2;
931 u2 = u1;
932 u1 = (u_uncomp[i] << 8) + u_uncomp[i+1];
933 delta = u1 - z;
934 comp[i ] = (delta >> 8) & 0xff;
935 comp[i+1] = (delta >> 0) & 0xff;
936 }
937 break;
938
939 case 3:
940 for (i = 0; i < uncomp_len; i+=2) {
941 z = 3*u1 - 3*u2 + u3;
942 u3 = u2;
943 u2 = u1;
944 u1 = (u_uncomp[i] << 8) + u_uncomp[i+1];
945 delta = u1 - z;
946 comp[i ] = (delta >> 8) & 0xff;
947 comp[i+1] = (delta >> 0) & 0xff;
948 }
949 break;
950
951 default:
952 return NULL;
953 }
954 comp-=2;
955 comp[0] = ZTR_FORM_DELTA2;
956 comp[1] = level;
957
958 *comp_len = uncomp_len+2;
959
960 return comp;
961 }
962
decorrelate2dyn(char * x_uncomp,int uncomp_len,int * comp_len)963 char *decorrelate2dyn(char *x_uncomp,
964 int uncomp_len,
965 int *comp_len) {
966 int i, j, z[4];
967 int u1 = 0, u2 = 0, u3 = 0;
968 char *comp = (char *)xmalloc(uncomp_len + 2);
969 unsigned char *u_uncomp = (unsigned char *)x_uncomp;
970 int level = 2; /* minimum level */
971 int last_level = level;
972 int best;
973
974 if (!comp)
975 return NULL;
976
977 comp+=2;
978 for (i = 0; i < uncomp_len; i+=2) {
979 z[0] = 0;
980 z[1] = u1;
981 z[2] = 2*u1 - u2;
982 z[3] = 3*u1 - 3*u2 + u3;
983 u3 = u2;
984 u2 = u1;
985 u1 = (u_uncomp[i]<<8) + u_uncomp[i+1];
986 comp[i ] = ((u1 - z[last_level]) >> 8) & 0xff;
987 comp[i+1] = ((u1 - z[last_level]) >> 0) & 0xff;
988 best = 10000;
989 for (j = level; j < 4; j++) {
990 if (ABS(u1 - z[j]) < best) {
991 best = ABS(u1 - z[j]);
992 last_level = j;
993 }
994 }
995 }
996 comp-=2;
997 comp[0] = ZTR_FORM_DDELTA2;
998 comp[1] = level;
999
1000 *comp_len = uncomp_len+2;
1001
1002 return comp;
1003 }
1004
1005 /*
1006 * recorrelate2()
1007 *
1008 * The reverse of decorrelate2()
1009 *
1010 * Arguments:
1011 * comp Compressed input data
1012 * comp_len Length of comp data
1013 * uncomp_len Output: length of uncompressed data
1014 *
1015 * Returns:
1016 * Success: uncompressed data
1017 * Failure: NULL
1018 */
recorrelate2(char * x_comp,int comp_len,int * uncomp_len)1019 char *recorrelate2(char *x_comp,
1020 int comp_len,
1021 int *uncomp_len) {
1022 int i, z;
1023 int u1 = 0, u2 = 0, u3 = 0;
1024 int level = x_comp[1];
1025 char *uncomp;
1026 unsigned char *u_comp = (unsigned char *)x_comp;
1027
1028 uncomp = (char *)xmalloc(comp_len-2);
1029 if (!uncomp)
1030 return NULL;
1031
1032 u_comp+=2;
1033 comp_len-=2;
1034 *uncomp_len = comp_len;
1035
1036 switch (level) {
1037 case 1:
1038 for (i = 0; i < comp_len; i+=2) {
1039 z = u1;
1040 u1 = ((u_comp[i] << 8) | u_comp[i+1]) + z;
1041 uncomp[i ] = (u1 >> 8) & 0xff;
1042 uncomp[i+1] = (u1 >> 0) & 0xff;
1043 }
1044 break;
1045
1046 case 2:
1047 for (i = 0; i < comp_len; i+=2) {
1048 z = 2*u1 - u2;
1049 u2 = u1;
1050 u1 = ((u_comp[i] << 8) | u_comp[i+1]) + z;
1051 uncomp[i ] = (u1 >> 8) & 0xff;
1052 uncomp[i+1] = (u1 >> 0) & 0xff;
1053 }
1054 break;
1055
1056 case 3:
1057 for (i = 0; i < comp_len; i+=2) {
1058 z = 3*u1 - 3*u2 + u3;
1059 u3 = u2;
1060 u2 = u1;
1061 u1 = ((u_comp[i] << 8) | u_comp[i+1]) + z;
1062 uncomp[i ] = (u1 >> 8) & 0xff;
1063 uncomp[i+1] = (u1 >> 0) & 0xff;
1064 }
1065 break;
1066 }
1067
1068 return uncomp;
1069 }
1070
1071 /*
1072 * ---------------------------------------------------------------------------
1073 * ZTR_FORM_DELTA4
1074 * ---------------------------------------------------------------------------
1075 */
1076
1077 /*
1078 * decorrelate4()
1079 *
1080 * Produce successive deltas from a 4-byte array (big endian)
1081 *
1082 * Arguments:
1083 * uncomp Uncompressed data
1084 * uncomp_len Length of uncompressed data
1085 * level Differencing level (must be 1, 2 or 3)
1086 * comp_len Return: where to store new compressed length
1087 *
1088 * Returns:
1089 * Success: A decorrelated buffer (malloced)
1090 * Failure: NULL
1091 */
decorrelate4(char * x_uncomp,int uncomp_len,int level,int * comp_len)1092 char *decorrelate4(char *x_uncomp,
1093 int uncomp_len,
1094 int level,
1095 int *comp_len) {
1096 int i, z, delta;
1097 int u1 = 0, u2 = 0, u3 = 0;
1098 char *comp = (char *)xmalloc(uncomp_len + 4);
1099 unsigned char *u_uncomp = (unsigned char *)x_uncomp;
1100
1101 if (!comp)
1102 return NULL;
1103
1104 comp+=4;
1105 switch (level) {
1106 case 1:
1107 for (i = 0; i < uncomp_len; i+=4) {
1108 z = u1;
1109 u1 =(u_uncomp[i ] << 24) +
1110 (u_uncomp[i+1] << 16) +
1111 (u_uncomp[i+2] << 8) +
1112 (u_uncomp[i+3] << 0);
1113 delta = u1 - z;
1114 comp[i ] = (delta >> 24) & 0xff;
1115 comp[i+1] = (delta >> 16) & 0xff;
1116 comp[i+2] = (delta >> 8) & 0xff;
1117 comp[i+3] = (delta >> 0) & 0xff;
1118 }
1119 break;
1120
1121 case 2:
1122 for (i = 0; i < uncomp_len; i+=4) {
1123 z = 2*u1 - u2;
1124 u2 = u1;
1125 u1 =(u_uncomp[i ] << 24) +
1126 (u_uncomp[i+1] << 16) +
1127 (u_uncomp[i+2] << 8) +
1128 (u_uncomp[i+3] << 0);
1129 delta = u1 - z;
1130 comp[i ] = (delta >> 24) & 0xff;
1131 comp[i+1] = (delta >> 16) & 0xff;
1132 comp[i+2] = (delta >> 8) & 0xff;
1133 comp[i+3] = (delta >> 0) & 0xff;
1134 }
1135 break;
1136
1137 case 3:
1138 for (i = 0; i < uncomp_len; i+=4) {
1139 z = 3*u1 - 3*u2 + u3;
1140 u3 = u2;
1141 u2 = u1;
1142 u1 =(u_uncomp[i ] << 24) +
1143 (u_uncomp[i+1] << 16) +
1144 (u_uncomp[i+2] << 8) +
1145 (u_uncomp[i+3] << 0);
1146 delta = u1 - z;
1147 comp[i ] = (delta >> 24) & 0xff;
1148 comp[i+1] = (delta >> 16) & 0xff;
1149 comp[i+2] = (delta >> 8) & 0xff;
1150 comp[i+3] = (delta >> 0) & 0xff;
1151 }
1152 break;
1153
1154 default:
1155 return NULL;
1156 }
1157 comp-=4;
1158 comp[0] = ZTR_FORM_DELTA4;
1159 comp[1] = level;
1160 comp[2] = 0; /* dummy - to align on 4-byte boundary */
1161 comp[3] = 0; /* dummy - to align on 4-byte boundary */
1162
1163 *comp_len = uncomp_len+4;
1164
1165 return comp;
1166 }
1167
1168 /*
1169 * recorrelate4()
1170 *
1171 * The reverse of decorrelate4()
1172 *
1173 * Arguments:
1174 * comp Compressed input data
1175 * comp_len Length of comp data
1176 * uncomp_len Output: length of uncompressed data
1177 *
1178 * Returns:
1179 * Success: uncompressed data
1180 * Failure: NULL
1181 */
recorrelate4(char * x_comp,int comp_len,int * uncomp_len)1182 char *recorrelate4(char *x_comp,
1183 int comp_len,
1184 int *uncomp_len) {
1185 int i, z;
1186 int u1 = 0, u2 = 0, u3 = 0;
1187 int level = x_comp[1];
1188 char *uncomp;
1189 unsigned char *u_comp = (unsigned char *)x_comp;
1190
1191 uncomp = (char *)xmalloc(comp_len-4);
1192 if (!uncomp)
1193 return NULL;
1194
1195 u_comp+=4;
1196 comp_len-=4;
1197 *uncomp_len = comp_len;
1198
1199 switch (level) {
1200 case 1:
1201 for (i = 0; i < comp_len; i+=4) {
1202 z = u1;
1203 u1 = z +
1204 ((u_comp[i ] << 24) |
1205 (u_comp[i+1] << 16) |
1206 (u_comp[i+2] << 8) |
1207 u_comp[i+3]);
1208 uncomp[i ] = (u1 >> 24) & 0xff;
1209 uncomp[i+1] = (u1 >> 16) & 0xff;
1210 uncomp[i+2] = (u1 >> 8) & 0xff;
1211 uncomp[i+3] = (u1 >> 0) & 0xff;
1212 }
1213 break;
1214
1215 case 2:
1216 for (i = 0; i < comp_len; i+=4) {
1217 z = 2*u1 - u2;
1218 u2 = u1;
1219 u1 = z +
1220 ((u_comp[i ] << 24) |
1221 (u_comp[i+1] << 16) |
1222 (u_comp[i+2] << 8) |
1223 u_comp[i+3]);
1224 uncomp[i ] = (u1 >> 24) & 0xff;
1225 uncomp[i+1] = (u1 >> 16) & 0xff;
1226 uncomp[i+2] = (u1 >> 8) & 0xff;
1227 uncomp[i+3] = (u1 >> 0) & 0xff;
1228 }
1229 break;
1230
1231 case 3:
1232 for (i = 0; i < comp_len; i+=4) {
1233 z = 3*u1 - 3*u2 + u3;
1234 u3 = u2;
1235 u2 = u1;
1236 u1 = z +
1237 ((u_comp[i ] << 24) |
1238 (u_comp[i+1] << 16) |
1239 (u_comp[i+2] << 8) |
1240 u_comp[i+3]);
1241 uncomp[i ] = (u1 >> 24) & 0xff;
1242 uncomp[i+1] = (u1 >> 16) & 0xff;
1243 uncomp[i+2] = (u1 >> 8) & 0xff;
1244 uncomp[i+3] = (u1 >> 0) & 0xff;
1245 }
1246 break;
1247 }
1248
1249 return uncomp;
1250 }
1251
1252
1253 /*
1254 * ---------------------------------------------------------------------------
1255 * ZTR_FORM_16TO8
1256 * ---------------------------------------------------------------------------
1257 */
1258
1259 /*
1260 * shrink_16to8()
1261 *
1262 * Stores an array of 16-bit (big endian) array elements in an 8-bit array.
1263 * We assume that most 16-bit elements encode numbers that fit in an 8-bit
1264 * value. When not possible, we store a marker followed by the 16-bit value
1265 * stored as multiple 8-bit values.
1266 *
1267 * uncomp Uncompressed data
1268 * uncomp_len Length of uncompressed data (in bytes)
1269 * comp_len Return: where to store new compressed length
1270 *
1271 * Returns:
1272 * Success: An 8-bit array (malloced)
1273 * Failure: NULL
1274 */
shrink_16to8(char * x_uncomp,int uncomp_len,int * comp_len)1275 char *shrink_16to8(char *x_uncomp, int uncomp_len, int *comp_len) {
1276 char *comp;
1277 int i, j, i16;
1278 signed char *s_uncomp = (signed char *)x_uncomp;
1279
1280 /* Allocation - worst case is 3 * (uncomp_len/2) + 1 */
1281 if (NULL == (comp = (char *)xmalloc(3 * (uncomp_len/2) + 1)))
1282 return NULL;
1283
1284 comp[0] = ZTR_FORM_16TO8;
1285 for (i = 0, j = 1; i < uncomp_len; i+=2) {
1286 i16 = (s_uncomp[i] << 8) | (unsigned char)s_uncomp[i+1];
1287 if (i16 >= -127 && i16 <= 127) {
1288 comp[j++] = i16;
1289 } else {
1290 comp[j++] = -128;
1291 comp[j++] = s_uncomp[i];
1292 comp[j++] = s_uncomp[i+1];
1293 }
1294 }
1295
1296 /* Reclaim unneeded memory */
1297 comp = xrealloc(comp, j);
1298
1299 *comp_len = j;
1300 return comp;
1301 }
1302
1303
1304 /*
1305 * expand_8to16()
1306 *
1307 * The opposite of the shrink_16to8() function.
1308 *
1309 * comp Compressed input data
1310 * comp_len Length of comp data (in bytes)
1311 * uncomp_len Output: length of uncompressed data (in bytes)
1312 *
1313 * Returns:
1314 * Success: Uncompressed data (char *)
1315 * Failure: NULL
1316 */
expand_8to16(char * x_comp,int comp_len,int * uncomp_len)1317 char *expand_8to16(char *x_comp, int comp_len, int *uncomp_len) {
1318 int i, j;
1319 char *uncomp;
1320 signed char *s_comp = (signed char *)x_comp;
1321
1322 /* Allocation - worst case is twice comp_len */
1323 if (NULL == (uncomp = (char *)xmalloc(comp_len*2)))
1324 return NULL;
1325
1326 #if 0
1327 for (i = 0, j = 1; j < comp_len; i+=2) {
1328 if (s_comp[j] != -128) {
1329 uncomp[i ] = s_comp[j] < 0 ? -1 : 0;
1330 uncomp[i+1] = s_comp[j++];
1331 } else {
1332 j++;
1333 uncomp[i ] = s_comp[j++];
1334 uncomp[i+1] = s_comp[j++];
1335 }
1336 }
1337 #endif
1338
1339 for (i = 0, j = 1; j < comp_len; i+=2) {
1340 if (s_comp[j] >= 0) {
1341 uncomp[i ] = 0;
1342 uncomp[i+1] = s_comp[j++];
1343 } else {
1344 if (s_comp[j] != -128) {
1345 uncomp[i+1] = s_comp[j++];
1346 uncomp[i ] = -1;
1347 } else {
1348 j++;
1349 uncomp[i ] = s_comp[j++];
1350 uncomp[i+1] = s_comp[j++];
1351 }
1352 }
1353 }
1354
1355 /* Reclaim unneeded memory */
1356 uncomp = xrealloc(uncomp, i);
1357
1358 *uncomp_len = i;
1359 return uncomp;
1360 }
1361
1362 /*
1363 * ---------------------------------------------------------------------------
1364 * ZTR_FORM_32TO8
1365 * ---------------------------------------------------------------------------
1366 */
1367
1368 /*
1369 * shrink_32to8()
1370 *
1371 * Stores an array of 32-bit (big endian) array elements in an 8-bit array.
1372 * We assume that most 32-bit elements encode numbers that fit in an 8-bit
1373 * value. When not possible, we store a marker followed by the 32-bit value
1374 * stored as multiple 8-bit values.
1375 *
1376 * uncomp Uncompressed data
1377 * uncomp_len Length of uncompressed data (in bytes)
1378 * comp_len Return: where to store new compressed length
1379 *
1380 * Returns:
1381 * Success: An 8-bit array (malloced)
1382 * Failure: NULL
1383 */
shrink_32to8(char * x_uncomp,int uncomp_len,int * comp_len)1384 char *shrink_32to8(char *x_uncomp, int uncomp_len, int *comp_len) {
1385 char *comp;
1386 int i, j, i32;
1387 signed char *s_uncomp = (signed char *)x_uncomp;
1388
1389 /* Allocation - worst case is 5 * (uncomp_len/4) + 1 */
1390 if (NULL == (comp = (char *)xmalloc(5 * (uncomp_len/4) + 1)))
1391 return NULL;
1392
1393 comp[0] = ZTR_FORM_32TO8;
1394 for (i = 0, j = 1; i < uncomp_len; i+=4) {
1395 i32 = (s_uncomp[i] << 24) |
1396 (s_uncomp[i+1] << 16) |
1397 (s_uncomp[i+2] << 8) |
1398 (unsigned char)s_uncomp[i+3];
1399 if (i32 >= -127 && i32 <= 127) {
1400 comp[j++] = i32;
1401 } else {
1402 comp[j++] = -128;
1403 comp[j++] = s_uncomp[i];
1404 comp[j++] = s_uncomp[i+1];
1405 comp[j++] = s_uncomp[i+2];
1406 comp[j++] = s_uncomp[i+3];
1407 }
1408 }
1409
1410 /* Reclaim unneeded memory */
1411 comp = xrealloc(comp, j);
1412
1413 *comp_len = j;
1414 return comp;
1415 }
1416
1417 /*
1418 * expand_8to32()
1419 *
1420 * The opposite of the shrink_32to8() function.
1421 *
1422 * comp Compressed input data
1423 * comp_len Length of comp data (in bytes)
1424 * uncomp_len Output: length of uncompressed data (in bytes)
1425 *
1426 * Returns:
1427 * Success: Uncompressed data (char *)
1428 * Failure: NULL
1429 */
expand_8to32(char * comp,int comp_len,int * uncomp_len)1430 char *expand_8to32(char *comp, int comp_len, int *uncomp_len) {
1431 int i, j;
1432 char *uncomp;
1433 signed char *s_comp = (signed char *)comp;
1434
1435 /* Allocation - worst case is four times comp_len */
1436 if (NULL == (uncomp = (char *)xmalloc(comp_len*4)))
1437 return NULL;
1438
1439 for (i = 0, j = 1; j < comp_len; i+=4) {
1440 if (s_comp[j] != -128) {
1441 uncomp[i ] = s_comp[j] < 0 ? -1 : 0;
1442 uncomp[i+1] = s_comp[j] < 0 ? -1 : 0;
1443 uncomp[i+2] = s_comp[j] < 0 ? -1 : 0;
1444 uncomp[i+3] = s_comp[j++];
1445 } else {
1446 j++;
1447 uncomp[i ] = s_comp[j++];
1448 uncomp[i+1] = s_comp[j++];
1449 uncomp[i+2] = s_comp[j++];
1450 uncomp[i+3] = s_comp[j++];
1451 }
1452 }
1453
1454 /* Reclaim unneeded memory */
1455 uncomp = xrealloc(uncomp, i);
1456
1457 *uncomp_len = i;
1458 return uncomp;
1459 }
1460
1461 /*
1462 * ---------------------------------------------------------------------------
1463 * ZTR_FORM_FOLLOW1
1464 * ---------------------------------------------------------------------------
1465 */
1466 static int follow_tab[256][256];
follow1(char * x_uncomp,int uncomp_len,int * comp_len)1467 char *follow1(char *x_uncomp,
1468 int uncomp_len,
1469 int *comp_len) {
1470 char *comp = (char *)xmalloc(uncomp_len + 256 + 1);
1471 unsigned char *u_uncomp = (unsigned char *)x_uncomp;
1472 signed char *s_uncomp = ((signed char *)x_uncomp);
1473 int i, j;
1474 char next[256];
1475 int count[256];
1476
1477 if (!comp)
1478 return NULL;
1479
1480 /* Count di-freqs */
1481 memset(follow_tab, 0, 256*256*sizeof(int));
1482 #if 0
1483 for (i = 0; i < uncomp_len-1; i++)
1484 follow_tab[u_uncomp[i]][u_uncomp[i+1]]++;
1485
1486 /* Pick the most frequent next byte from the preceeding byte */
1487 for (i = 0; i < 256; i++) {
1488 int bestval, bestind;
1489
1490 bestval = bestind = 0;
1491 for (j = 0; j < 256; j++) {
1492 if (follow_tab[i][j] > bestval) {
1493 bestval = follow_tab[i][j];
1494 bestind = j;
1495 }
1496 }
1497 next[i] = bestind;
1498 }
1499 #endif
1500 memset(next, 0, 256*sizeof(*next));
1501 memset(count, 0, 256*sizeof(*count));
1502 /* Pick the most frequent next byte from the preceeding byte */
1503 for (i = 0; i < uncomp_len-1; ) {
1504 int cur = u_uncomp[i];
1505 int nxt = u_uncomp[++i];
1506 int folcnt = ++follow_tab[cur][nxt];
1507 if (folcnt > count[cur]) {
1508 count[cur] = folcnt;
1509 next[cur] = nxt;
1510 }
1511 }
1512
1513 j = 0;
1514 comp[j++] = ZTR_FORM_FOLLOW1;
1515
1516 /* Output 'next' array */
1517 for (i = 0; i < 256; i++, j++)
1518 comp[j] = next[i];
1519
1520 /* Output new 'uncomp' as next['uncomp'] */
1521 comp[j++] = u_uncomp[0];
1522 for (i = 1; i < uncomp_len; i++, j++) {
1523 comp[j] = next[u_uncomp[i-1]] - s_uncomp[i];
1524 }
1525 *comp_len = j;
1526
1527 return comp;
1528 }
1529
unfollow1(char * x_comp,int comp_len,int * uncomp_len)1530 char *unfollow1(char *x_comp,
1531 int comp_len,
1532 int *uncomp_len) {
1533 unsigned char *u_uncomp;
1534 int i, j;
1535 char next[256];
1536 unsigned char *u_comp = (unsigned char *)x_comp;
1537 signed char *s_comp = (signed char *)x_comp;
1538
1539 u_uncomp = (unsigned char *)xmalloc(comp_len-256-1);
1540 if (!u_uncomp)
1541 return NULL;
1542
1543 /* Load next[] array */
1544 j = 1;
1545 for (i = 0; i < 256; i++, j++)
1546 next[i] = u_comp[j];
1547
1548 /* Replace comp[x] with next[comp[x-1]] - comp[x]*/
1549 u_uncomp[0] = u_comp[j++];
1550
1551 comp_len -= 257;
1552 s_comp += 257;
1553 for (i = 1; i < comp_len; i++) {
1554 u_uncomp[i] = next[u_uncomp[i-1]] - s_comp[i];
1555 }
1556
1557 *uncomp_len = i;
1558
1559 return (char *)u_uncomp;
1560 }
1561
1562
1563 /*
1564 * ---------------------------------------------------------------------------
1565 * ZTR_FORM_CHEB445
1566 * ---------------------------------------------------------------------------
1567 */
1568
1569 #if 0
1570 /*
1571 * Compresses using chebyshev polynomials to predict the next peak.
1572 * Based on around 96 modern ABI files it compresses by around 5% (varied
1573 * between 3.9 and 6.6).
1574 */
1575
1576 /*
1577 * For now this has been disabled in favour of the integer version below
1578 * as we cannot guarantee all systems to have the same floating point
1579 * roundings, especially with the final conversion to integer.
1580 * (Also, for some unknown reason, the integer version produces smaller
1581 * files.)
1582 */
1583
1584 char *cheb445comp(char *uncomp, int uncomp_len, int *data_len)
1585 {
1586 int i, k, l, z;
1587 int datap;
1588 float frac[5];
1589 float fz[25];
1590 signed short *d16 = (signed short *)uncomp;
1591 int nwords = uncomp_len / 2;
1592 short *dptr = d16;
1593 signed short *data = (signed short *)malloc((nwords+1)*sizeof(short));
1594
1595 data[0] = le_int2(ZTR_FORM_CHEB445);
1596 /* Check for boundary cases */
1597 if (nwords <= 4) {
1598 switch (nwords) {
1599 case 4:
1600 data[4] = be_int2(be_int2(d16[3])-be_int2(d16[2]));
1601 case 3:
1602 data[3] = be_int2(be_int2(d16[2])-be_int2(d16[1]));
1603 case 2:
1604 data[2] = be_int2(be_int2(d16[1])-be_int2(d16[0]));
1605 case 1:
1606 data[1] = be_int2(d16[0]);
1607 }
1608 *data_len = nwords*2;
1609 return (char *)data;
1610 }
1611
1612 /* First 4 values are just direct deltas */
1613 data[1] = be_int2(d16[0]);
1614 data[2] = be_int2(be_int2(d16[1])-be_int2(d16[0]));
1615 data[3] = be_int2(be_int2(d16[2])-be_int2(d16[1]));
1616 data[4] = be_int2(be_int2(d16[3])-be_int2(d16[2]));
1617 datap = 5;
1618
1619 /* Initialise - speeds up loop */
1620 for (k = 0; k < 5; k++) {
1621 float kx = cos(M_PI*(k+0.5)/5)*1.5;
1622 frac[k] = (kx + 1.5) - (int)(kx + 1.5);
1623 }
1624 for (z = l = 0; l < 5; l++) {
1625 for (k = 0; k < 5; k++, z++) {
1626 fz[z] = 0.4 * cos(l * M_PI*(k+0.5)/5);
1627 }
1628 }
1629
1630 /* Loop */
1631 for (i = 0; i < nwords-4; i++) {
1632 float dd, y = 10/3.0;
1633 float f[5], coef[5];
1634 signed short diff;
1635 int p;
1636
1637 f[0] = be_int2(dptr[2])*frac[4] + be_int2(dptr[3])*frac[0];
1638 f[1] = be_int2(dptr[2])*frac[3] + be_int2(dptr[3])*frac[1];
1639 f[2] = be_int2(dptr[1])*frac[2] + be_int2(dptr[2])*frac[2];
1640 f[3] = be_int2(dptr[0])*frac[1] + be_int2(dptr[1])*frac[3];
1641 f[4] = be_int2(dptr[0])*frac[0] + be_int2(dptr[1])*frac[4];
1642
1643 for (z = l = 0; l < 5; l++, z+=5)
1644 coef[l] = f[0] * fz[z+0] +
1645 f[1] * fz[z+1] +
1646 f[2] * fz[z+2] +
1647 f[3] * fz[z+3] +
1648 f[4] * fz[z+4];
1649
1650 dd = y*coef[3]+coef[2];
1651 p = 0.5 + 5/3.0*(y*dd-coef[3]+coef[1])-dd+coef[0]/2.0;
1652
1653 if (p < 0) p = 0;
1654 diff = be_int2(dptr[4]) - p;
1655 data[datap++] = be_int2(diff);
1656
1657 dptr++;
1658 }
1659
1660 *data_len = datap*2;
1661
1662 return (char *)data;
1663 }
1664
1665 char *cheb445uncomp(char *comp, int comp_len, int *uncomp_len)
1666 {
1667 int i, k, l, z;
1668 float frac[5];
1669 float fz[25];
1670 signed short *d16 = (signed short *)comp;
1671 int nwords = comp_len / 2;
1672 signed short *data = (signed short *)xmalloc(comp_len);
1673 short *dptr = data, *dptr2 = d16;
1674
1675 /* Check for boundary cases */
1676 if (nwords <= 3) {
1677 switch (nwords) {
1678 case 3:
1679 data[0] = be_int2(d16[1]);
1680 data[1] = be_int2(be_int2(d16[2])+be_int2(data[0]));
1681 data[2] = be_int2(be_int2(d16[3])+be_int2(data[1]));
1682 break;
1683 case 2:
1684 data[0] = be_int2(d16[1]);
1685 data[1] = be_int2(be_int2(d16[2])+be_int2(data[0]));
1686 break;
1687 case 1:
1688 data[0] = be_int2(d16[1]);
1689 break;
1690 }
1691 *uncomp_len = (nwords-1)*2;
1692 return (char *)data;
1693 }
1694
1695 /* First 3 values are just direct deltas */
1696 data[0] = be_int2(d16[1]);
1697 data[1] = be_int2(be_int2(d16[2])+be_int2(data[0]));
1698 data[2] = be_int2(be_int2(d16[3])+be_int2(data[1]));
1699 data[3] = be_int2(be_int2(d16[4])+be_int2(data[2]));
1700 dptr2 += 5;
1701
1702 /* Initialise - speeds up loop */
1703 for (k = 0; k < 5; k++) {
1704 float kx = cos(M_PI*(k+0.5)/5)*1.5;
1705 frac[k] = (kx + 1.5) - (int)(kx + 1.5);
1706 }
1707 for (z = l = 0; l < 5; l++) {
1708 for (k = 0; k < 5; k++, z++) {
1709 fz[z] = 0.4 * cos(l * M_PI*(k+0.5)/5);
1710 }
1711 }
1712
1713 /* Loop */
1714 for (i = 0; i < nwords-3; i++) {
1715 float dd, y = 10/3.0;
1716 float f[5], coef[5];
1717 signed short diff;
1718 int p;
1719
1720 f[0] = be_int2(dptr[2])*frac[4] + be_int2(dptr[3])*frac[0];
1721 f[1] = be_int2(dptr[2])*frac[3] + be_int2(dptr[3])*frac[1];
1722 f[2] = be_int2(dptr[1])*frac[2] + be_int2(dptr[2])*frac[2];
1723 f[3] = be_int2(dptr[0])*frac[1] + be_int2(dptr[1])*frac[3];
1724 f[4] = be_int2(dptr[0])*frac[0] + be_int2(dptr[1])*frac[4];
1725
1726 for (z = l = 0; l < 5; l++, z+=5)
1727 coef[l] = f[0] * fz[z+0] +
1728 f[1] * fz[z+1] +
1729 f[2] * fz[z+2] +
1730 f[3] * fz[z+3] +
1731 f[4] * fz[z+4];
1732
1733 dd = y*coef[3]+coef[2];
1734 p = 0.5 + 5/3.0*(y*dd-coef[3]+coef[1])-dd+coef[0]/2.0;
1735
1736 if (p < 0) p = 0;
1737 diff = be_int2(*dptr2) + p;
1738 dptr[4] = be_int2(diff);
1739
1740 dptr++;
1741 dptr2++;
1742 }
1743
1744 *uncomp_len = (nwords-1)*2;
1745 return (char *)data;
1746 }
1747 #endif
1748
1749 /*
1750 * ---------------------------------------------------------------------------
1751 * ZTR_FORM_ICHEB
1752 * ---------------------------------------------------------------------------
1753 */
1754
1755 /*
1756 * Integer versions of the chebyshev polynomial compressor. This uses
1757 * the polynomials to predict the next peak from the preceeding 3.
1758 * Tested on 100 ABI-3700, Megabace and Licor files it compressed by
1759 * around 7-8%. (Oddly this is slightly more than the floating point
1760 * version.)
1761 *
1762 * These require 32-bit integers and have code to make sure that arithmetic
1763 * does not overflow this.
1764 */
1765 #define CH1 150
1766 #define CH2 105
ichebcomp(char * uncomp,int uncomp_len,int * data_len)1767 char *ichebcomp(char *uncomp, int uncomp_len, int *data_len)
1768 {
1769 int i, l, z;
1770 int datap;
1771 int frac[5] = {139,57,75,93,11};
1772 int fz[20] = {42, 42, 42, 42, 42,
1773 39, 24, 0,-24,-39,
1774 33,-12,-42,-12, 33,
1775 24,-39, 0, 39,-24};
1776 int dfac;
1777 signed short *d16 = (signed short *)uncomp;
1778 int nwords = uncomp_len / 2;
1779 short *dptr = d16;
1780 signed short *data = (signed short *)malloc((nwords+1)*sizeof(short));
1781
1782 data[0] = le_int2(ZTR_FORM_ICHEB);
1783 /* Check for boundary cases */
1784 if (nwords <= 4) {
1785 switch (nwords) {
1786 case 4:
1787 data[4] = be_int2(be_int2(d16[3])-be_int2(d16[2]));
1788 case 3:
1789 data[3] = be_int2(be_int2(d16[2])-be_int2(d16[1]));
1790 case 2:
1791 data[2] = be_int2(be_int2(d16[1])-be_int2(d16[0]));
1792 case 1:
1793 data[1] = be_int2(d16[0]);
1794 }
1795 *data_len = nwords*2;
1796 return (char *)data;
1797 }
1798
1799 /* First 4 values are just direct deltas */
1800 data[1] = be_int2(d16[0]);
1801 data[2] = be_int2(be_int2(d16[1])-be_int2(d16[0]));
1802 data[3] = be_int2(be_int2(d16[2])-be_int2(d16[1]));
1803 data[4] = be_int2(be_int2(d16[3])-be_int2(d16[2]));
1804 datap = 5;
1805
1806 /* Loop */
1807 for (i = 4; i < nwords; i++) {
1808 int dd, f[5];
1809 signed int coef[4];
1810 signed short diff;
1811 int p;
1812
1813 /*
1814 * FIXME: As an alternative to the range checking below, if we
1815 * scale dptr[X] such that it's never higher than 2800 then
1816 * the 32-bit arithmetic will never overflow. Practically speaking,
1817 * all observed ABI and Megabace files have vales up to 1600 only.
1818 */
1819
1820 /*
1821 * frac[N] is always paired with frac[4-N], summing to 1.0 - or
1822 * 150 when scaled.
1823 * Hence f[0] has range 0 to 65536*150.
1824 */
1825 f[0] = ((unsigned short)be_int2(dptr[2]))*frac[4] +
1826 ((unsigned short)be_int2(dptr[3]))*frac[0];
1827 f[1] = ((unsigned short)be_int2(dptr[2]))*frac[3] +
1828 ((unsigned short)be_int2(dptr[3]))*frac[1];
1829 f[2] = ((unsigned short)be_int2(dptr[1]))*frac[2] +
1830 ((unsigned short)be_int2(dptr[2]))*frac[2];
1831 f[3] = ((unsigned short)be_int2(dptr[0]))*frac[1] +
1832 ((unsigned short)be_int2(dptr[1]))*frac[3];
1833 f[4] = ((unsigned short)be_int2(dptr[0]))*frac[0] +
1834 ((unsigned short)be_int2(dptr[1]))*frac[4];
1835
1836 /*
1837 * fz[z+0..5] sums to no more than 210 (5*42) and no less than 0.
1838 * Therefore coef[l] has range 0 to 65536*150*210, which (just)
1839 * fits in 31-bits, plus 1 for the sign.
1840 */
1841 for (z = l = 0; l < 4; l++, z+=5)
1842 coef[l] = (f[0] * fz[z+0] +
1843 f[1] * fz[z+1] +
1844 f[2] * fz[z+2] +
1845 f[3] * fz[z+3] +
1846 f[4] * fz[z+4]);
1847
1848 /*
1849 * computing p requires at most a temporary variable of
1850 * 24.1 * coef, but coef may be a full 32-bit integer.
1851 * If coef is sufficiently close to cause an integer overflow then
1852 * we scale it down.
1853 */
1854 {
1855 int max = 0;
1856
1857 for (l = 0; l < 4; l++) {
1858 if (max < ABS(coef[l]))
1859 max = ABS(coef[l]);
1860 }
1861
1862
1863 if (max > 1<<26) {
1864 dfac = max / (1<<26) + 1;
1865 for (l = 0; l < 4; l++)
1866 coef[l] /= dfac;
1867 } else {
1868 dfac = 1;
1869 }
1870 }
1871
1872 dd = (coef[3]/3)*10+coef[2];
1873 p = ((((dd/3)*10-coef[3]+coef[1])/3)*5-dd+coef[0]/2)/(CH1*CH2);
1874 p *= dfac;
1875
1876 if (p < 0) p = 0;
1877 diff = be_int2(dptr[4]) - p;
1878 data[datap++] = be_int2(diff);
1879
1880 dptr++;
1881 }
1882
1883 *data_len = datap*2;
1884
1885 return (char *)data;
1886 }
1887
ichebuncomp(char * comp,int comp_len,int * uncomp_len)1888 char *ichebuncomp(char *comp, int comp_len, int *uncomp_len)
1889 {
1890 int i, l, z;
1891 int frac[5] = {139,57,75,93,11};
1892 int fz[20] = {42, 42, 42, 42, 42,
1893 39, 24, 0,-24,-39,
1894 33,-12,-42,-12, 33,
1895 24,-39, 0, 39,-24};
1896 signed short *d16 = (signed short *)comp;
1897 int nwords = comp_len / 2 - 1;
1898 signed short *data = (signed short *)xmalloc(comp_len);
1899 short *dptr = data, *dptr2 = d16;
1900 int dfac;
1901
1902 /* Check for boundary cases */
1903 if (nwords <= 4) {
1904 switch (nwords) {
1905 case 4:
1906 data[0] = be_int2(d16[1]);
1907 data[1] = be_int2(be_int2(d16[2])+be_int2(data[0]));
1908 data[2] = be_int2(be_int2(d16[3])+be_int2(data[1]));
1909 data[3] = be_int2(be_int2(d16[4])+be_int2(data[2]));
1910 break;
1911 case 3:
1912 data[0] = be_int2(d16[1]);
1913 data[1] = be_int2(be_int2(d16[2])+be_int2(data[0]));
1914 data[2] = be_int2(be_int2(d16[3])+be_int2(data[1]));
1915 break;
1916 case 2:
1917 data[0] = be_int2(d16[1]);
1918 data[1] = be_int2(be_int2(d16[2])+be_int2(data[0]));
1919 break;
1920 case 1:
1921 data[0] = be_int2(d16[1]);
1922 break;
1923 }
1924 *uncomp_len = nwords*2;
1925 return (char *)data;
1926 }
1927
1928 /* First 4 values are just direct deltas */
1929 data[0] = be_int2(d16[1]);
1930 data[1] = be_int2(be_int2(d16[2])+be_int2(data[0]));
1931 data[2] = be_int2(be_int2(d16[3])+be_int2(data[1]));
1932 data[3] = be_int2(be_int2(d16[4])+be_int2(data[2]));
1933 dptr2 += 5;
1934
1935 /* Loop */
1936 for (i = 4; i < nwords; i++) {
1937 int dd, coef[5], f[5];
1938 signed short diff;
1939 int p;
1940
1941 f[0] = ((unsigned short)be_int2(dptr[2]))*frac[4] +
1942 ((unsigned short)be_int2(dptr[3]))*frac[0];
1943 f[1] = ((unsigned short)be_int2(dptr[2]))*frac[3] +
1944 ((unsigned short)be_int2(dptr[3]))*frac[1];
1945 f[2] = ((unsigned short)be_int2(dptr[1]))*frac[2] +
1946 ((unsigned short)be_int2(dptr[2]))*frac[2];
1947 f[3] = ((unsigned short)be_int2(dptr[0]))*frac[1] +
1948 ((unsigned short)be_int2(dptr[1]))*frac[3];
1949 f[4] = ((unsigned short)be_int2(dptr[0]))*frac[0] +
1950 ((unsigned short)be_int2(dptr[1]))*frac[4];
1951
1952 for (z = l = 0; l < 4; l++, z+=5)
1953 coef[l] = f[0] * fz[z+0] +
1954 f[1] * fz[z+1] +
1955 f[2] * fz[z+2] +
1956 f[3] * fz[z+3] +
1957 f[4] * fz[z+4];
1958
1959 /*
1960 * computing p requires at most a temporary variable of
1961 * 24.1 * coef, but coef may be a full 32-bit integer.
1962 * If coef is sufficiently close to cause an integer overflow then
1963 * we scale it down.
1964 */
1965 {
1966 int max = 0;
1967
1968 for (l = 0; l < 4; l++) {
1969 if (max < ABS(coef[l]))
1970 max = ABS(coef[l]);
1971 }
1972
1973
1974 if (max > 1<<26) {
1975 dfac = max / (1<<26) + 1;
1976 for (l = 0; l < 4; l++)
1977 coef[l] /= dfac;
1978 } else {
1979 dfac = 1;
1980 }
1981 }
1982
1983 dd = (coef[3]/3)*10+coef[2];
1984 p = ((((dd/3)*10-coef[3]+coef[1])/3)*5-dd+coef[0]/2)/(CH1*CH2);
1985 p *= dfac;
1986
1987 if (p < 0) p = 0;
1988 diff = be_int2(*dptr2) + p;
1989 dptr[4] = be_int2(diff);
1990
1991 dptr++;
1992 dptr2++;
1993 }
1994
1995 *uncomp_len = nwords*2;
1996 return (char *)data;
1997 }
1998
1999 /*
2000 * ---------------------------------------------------------------------------
2001 * ZTR_FORM_LOG
2002 * ---------------------------------------------------------------------------
2003 */
2004
2005 /*
2006 * This is a LOSSY compression. It replaces N with 10 * log2(N).
2007 * (Just an idea, and not a great one it seems.)
2008 */
log2_data(char * x_uncomp,int uncomp_len,int * comp_len)2009 char *log2_data(char *x_uncomp,
2010 int uncomp_len,
2011 int *comp_len) {
2012 int i, u1, l1;
2013 char *comp = (char *)xmalloc(uncomp_len + 2);
2014 unsigned char *u_uncomp = (unsigned char *)x_uncomp;
2015
2016 if (!comp)
2017 return NULL;
2018
2019 comp+=2;
2020 for (i = 0; i < uncomp_len; i+=2) {
2021 u1 =(u_uncomp[i ] << 8) +
2022 (u_uncomp[i+1] << 0);
2023 l1 = (int)(10 * log(u1+1) / log(2));
2024 comp[i ] = (l1 >> 8) & 0xff;
2025 comp[i+1] = (l1 >> 0) & 0xff;
2026 }
2027
2028 comp-=2;
2029 comp[0] = ZTR_FORM_LOG2;
2030 comp[1] = 0; /* dummy - to align on 2-byte boundary */
2031
2032 *comp_len = uncomp_len+2;
2033
2034 return comp;
2035 }
2036
unlog2_data(char * x_comp,int comp_len,int * uncomp_len)2037 char *unlog2_data(char *x_comp,
2038 int comp_len,
2039 int *uncomp_len) {
2040 int i, u1, l1;
2041 char *uncomp;
2042 unsigned char *u_comp = (unsigned char *)x_comp;
2043
2044 uncomp = (char *)xmalloc(comp_len-2);
2045 if (!uncomp)
2046 return NULL;
2047
2048 u_comp+=2;
2049 comp_len-=2;
2050 *uncomp_len = comp_len;
2051
2052 for (i = 0; i < comp_len; i+=2) {
2053 l1 = ((u_comp[i ] << 8) |
2054 (u_comp[i+1] << 0));
2055 u1 = (int)pow(2.0, l1/10.0)-1;
2056 uncomp[i ] = (u1 >> 8) & 0xff;
2057 uncomp[i+1] = (u1 >> 0) & 0xff;
2058 }
2059
2060 return uncomp;
2061 }
2062
2063 /*
2064 * ---------------------------------------------------------------------------
2065 * ZTR_FORM_STHUFF
2066 * ---------------------------------------------------------------------------
2067 */
2068
2069 /*
2070 * Implements compression using a set of static huffman codes stored using
2071 * the Deflate algorithm (and so in this respect it's similar to zlib).
2072 *
2073 * The huffman codes though can be previously stored in the ztr object
2074 * using ztr_add_hcode(). "cset" indicates which numbered stored huffman
2075 * code set is to be used, or passing zero will use inline codes (ie they
2076 * are stored in the data stream itself, just as in standard deflate).
2077 *
2078 * Arguments:
2079 * ztr ztr_t pointer; used to find stored code-sets
2080 * uncomp The uncompressed input data
2081 * uncomp_len Length of uncomp
2082 * cset Stored code-set number, zero for inline
2083 * recsz Record size - only used when cset == 0.
2084 * comp_len Output: length of compressed data
2085 *
2086 * Returns:
2087 * Compressed data stream if successful + comp_len
2088 * NULL on failure
2089 */
sthuff(ztr_t * ztr,char * uncomp,int uncomp_len,int cset,int recsz,int * comp_len)2090 char *sthuff(ztr_t *ztr, char *uncomp, int uncomp_len,
2091 int cset, int recsz, int *comp_len) {
2092 block_t *blk = block_create(NULL, 2);
2093 unsigned char bytes[2];
2094 huffman_codeset_t *c = NULL;
2095 unsigned char *comp = NULL;
2096 ztr_hcode_t *hc = NULL;
2097
2098 if (cset >= CODE_USER) {
2099 if (NULL == (hc = ztr_find_hcode(ztr, cset)))
2100 return NULL;
2101 c = hc->codes;
2102 } else if (cset != CODE_INLINE) {
2103 c = generate_code_set(cset, 1, NULL, 0, 1, MAX_CODE_LEN, 0);
2104 }
2105
2106 if (!c) {
2107 /* No cached ones found, so inline some instead */
2108 cset = 0;
2109 c = generate_code_set(0, recsz, (unsigned char *)uncomp, uncomp_len, 1,
2110 MAX_CODE_LEN, 0);
2111 }
2112
2113 bytes[0] = ZTR_FORM_STHUFF;
2114 bytes[1] = cset;
2115 store_bytes(blk, bytes, 2);
2116
2117 if (hc) {
2118 if (!c->blk) {
2119 c->blk = block_create(NULL, 2);
2120 store_codes(c->blk, c, 1);
2121 }
2122 blk->bit = c->blk->bit;
2123 } else {
2124 store_codes(blk, c, 1);
2125 }
2126
2127 /*
2128 {int i;
2129 for (i = 0; i < c->ncodes; i++) {
2130 output_code_set(stderr, c->codes[i]);
2131 }}
2132 */
2133
2134
2135 /*
2136 * Unless CODE_INLINE, all we wanted to know is what bit number
2137 * to start on. The above is therefore somewhat inefficient.
2138 */
2139 if (cset != 0) {
2140 blk->byte = 2;
2141 memset(&blk->data[2], 0, blk->alloc - 2);
2142 }
2143
2144 if (0 == huffman_multi_encode(blk, c, cset,
2145 (unsigned char *)uncomp, uncomp_len)) {
2146 comp = blk->data;
2147 *comp_len = blk->byte + (blk->bit != 0);
2148 block_destroy(blk, 1);
2149 }
2150
2151 if (cset == 0)
2152 huffman_codeset_destroy(c);
2153
2154 return (char *)comp;
2155 }
2156
unsthuff(ztr_t * ztr,char * comp,int comp_len,int * uncomp_len)2157 char *unsthuff(ztr_t *ztr, char *comp, int comp_len, int *uncomp_len) {
2158 int cset = (unsigned char)(comp[1]);
2159 huffman_codeset_t *cs = NULL, *cs_free = NULL;
2160 block_t *blk_in = block_create(NULL, comp_len),
2161 *blk_out = block_create(NULL, 1000);
2162 int bfinal = 1, bit_num = 0;
2163 char *uncomp;
2164
2165 if (cset >= CODE_USER) {
2166 /* Scans through HUFF chunks */
2167 ztr_hcode_t *hc = ztr_find_hcode(ztr, cset);
2168 if (!hc)
2169 return NULL;
2170
2171 cs = hc->codes;
2172 bit_num = cs->bit_num;
2173 blk_in->bit = 0;
2174 } else if (cset > 0) {
2175 /* Create some temporary huffman_codes to stringify */
2176 cs_free = cs = generate_code_set(cset, 1, NULL, 0, 1, MAX_CODE_LEN, 0);
2177 if (!cs)
2178 return NULL;
2179
2180 bit_num = cs->bit_num;
2181 blk_in->bit = 0;
2182 } /* else inline codes */
2183
2184 /*
2185 * We need to know at what bit the huffman codes would have ended on
2186 * so we can store our huffman encoded symbols immediately following it.
2187 * For speed though this bit-number is cached.
2188 */
2189 blk_in->data[blk_in->byte++] |= *(comp+2);
2190 store_bytes(blk_in, (unsigned char *)comp+3, comp_len-3);
2191
2192
2193 /* Rewind */
2194 blk_in->byte = 0;
2195 blk_in->bit = bit_num;
2196
2197 do {
2198 block_t *out;
2199
2200 /*
2201 * We're either at the start of a block with codes to restore
2202 * (cset == INLINE or the 2nd onwards block) or we've already
2203 * got some codes in cs and we're at the position where huffman
2204 * encoded symbols are stored.
2205 */
2206 if (!cs)
2207 if (NULL == (cs = cs_free = restore_codes(blk_in, &bfinal)))
2208 return NULL;
2209
2210 /*
2211 {int i;
2212 for (i = 0; i < cs->ncodes; i++) {
2213 output_code_set(stderr, cs->codes[i]);
2214 }}
2215 */
2216
2217 if (NULL == (out = huffman_multi_decode(blk_in, cs))) {
2218 huffman_codeset_destroy(cs);
2219 return NULL;
2220 }
2221
2222 /* Could optimise this for the common case of only 1 block */
2223 store_bytes(blk_out, out->data, out->byte);
2224 block_destroy(out, 0);
2225 if (cs_free)
2226 huffman_codeset_destroy(cs_free);
2227 cs = cs_free = NULL;
2228 } while (!bfinal);
2229
2230 *uncomp_len = blk_out->byte;
2231 uncomp = (char *)blk_out->data;
2232
2233 block_destroy(blk_in, 0);
2234 block_destroy(blk_out, 1);
2235
2236 return uncomp;
2237 }
2238
2239
2240 #ifndef NDEBUG
2241 #define SYM_EOF 256
output_code_set(FILE * fp,huffman_codes_t * cds)2242 static void output_code_set(FILE *fp, huffman_codes_t *cds) {
2243 int i, j;
2244 int nbits_in = 0, nbits_out = 0;
2245 huffman_code_t *codes = cds->codes;
2246 int ncodes = cds->ncodes;
2247
2248 fprintf(fp, "static huffman_code_t codes_FIXME[] = {\n");
2249 for (i = j = 0; i < ncodes; i++) {
2250 nbits_out += codes[i].nbits * codes[i].freq;
2251 nbits_in += 8*codes[i].freq;
2252 if (j == 0)
2253 fprintf(fp, " ");
2254 if (codes[i].symbol == SYM_EOF) {
2255 fprintf(fp, "{SYM_EOF,%3d}, ", codes[i].nbits);
2256 j = 10;
2257 } else {
2258 if (isalnum(codes[i].symbol)) {
2259 fprintf(fp, "{'%c',%3d}, ", codes[i].symbol, codes[i].nbits);
2260 } else {
2261 fprintf(fp, "{%3d,%3d}, ", codes[i].symbol, codes[i].nbits);
2262 }
2263 }
2264 j++;
2265
2266 if (j >= 6) {
2267 fputc('\n', fp);
2268 j = 0;
2269 }
2270 }
2271 if (j)
2272 fputc('\n', fp);
2273 fprintf(fp, "};\n");
2274 fprintf(fp, "/* Expected compression to %f of input */\n",
2275 (double)nbits_out/nbits_in);
2276 }
2277 #endif
2278
2279 /*
2280 * Reorders quality data from its RAW format to an interleaved 4-byte
2281 * aligned format.
2282 *
2283 * Starting with sequence A1 C2 G3 the raw format is quality of called
2284 * bases followed by quality of remaining bases in triplets per base:
2285 * 0 (RAW format)
2286 * Q(A1) Q(C2) Q(G3)
2287 * Q(C1) Q(G1) Q(T1)
2288 * Q(A2) Q(G2) Q(T2)
2289 * Q(A3) Q(C3) Q(T3)
2290 *
2291 * We reorder it to:
2292 * ZTR_FORM_QSHIFT <any> <any> 0(raw)
2293 * Q(A1) Q(C1) Q(G1) Q(T1)
2294 * Q(C2) Q(A2) Q(G2) Q(T2)
2295 * Q(G3) Q(A3) Q(C3) Q(T3)
2296 *
2297 * Returns shifted data on success
2298 * NULL on failure
2299 */
qshift(char * qold,int qlen,int * new_len)2300 char *qshift(char *qold, int qlen, int *new_len) {
2301 int i, j, k;
2302 char *qnew;
2303 int nbases;
2304
2305 /*
2306 * Correct input is raw encoding + 4x nbases bytes
2307 */
2308 if ((qlen-1)%4 != 0 || *qold != 0)
2309 return NULL;
2310
2311 nbases = (qlen-1)/4;
2312 qnew = (char *)malloc((nbases+1)*4);
2313 qnew[0] = ZTR_FORM_QSHIFT; /* reorder code */
2314 qnew[1] = -40; /* pad */
2315 qnew[2] = -40; /* pad */
2316 qnew[3] = qold[0];
2317
2318 for (i = 0, j = 4, k = nbases; i < nbases; i++, j+=4, k+=3) {
2319 qnew[j ] = qold[1+i ];
2320 qnew[j+1] = qold[1+k ];
2321 qnew[j+2] = qold[1+k+1];
2322 qnew[j+3] = qold[1+k+2];
2323 }
2324
2325 *new_len = (nbases+1)*4;
2326 return qnew;
2327 }
2328
2329 /*
2330 * The opposite transform from qshift() above.
2331 *
2332 * Returns unshifted data on success
2333 * NULL on failure.
2334 */
unqshift(char * qold,int qlen,int * new_len)2335 char *unqshift(char *qold, int qlen, int *new_len) {
2336 int i, j, k;
2337 char *qnew;
2338 int nbases;
2339
2340 /*
2341 * Correct input is 4x (nbases+1) bytes
2342 */
2343 if (qlen%4 != 0 || *qold != ZTR_FORM_QSHIFT)
2344 return NULL;
2345
2346 nbases = qlen/4-1;
2347 qnew = (char *)malloc(nbases*4+1);
2348 qnew[0] = 0; /* raw byte */
2349
2350 for (i = 0, j = 4, k = nbases; i < nbases; i++, j+=4, k+=3) {
2351 qnew[1+i ] = qold[j ];
2352 qnew[1+k ] = qold[j+1];
2353 qnew[1+k+1] = qold[j+2];
2354 qnew[1+k+2] = qold[j+3];
2355 }
2356
2357 *new_len = nbases*4+1;
2358 return qnew;
2359 }
2360
2361 /*
2362 * Given a sequence ACTG this shifts trace data from the order:
2363 *
2364 * A1A2A3A4 C1C2C3C4 G1G2G3G4 T1T2T3T4
2365 *
2366 * to
2367 *
2368 * A1C1G1T1 C2A2G2T2 T3A3C3G3 G4C4C4T4
2369 *
2370 * Ie for each base it ouputs the signal for the called base first
2371 * followed by the remaining 3 signals in A,C,G,T order (minus the
2372 * called signal already output).
2373 */
2374
2375 /*
2376 * NCBI uses a rotation mechanism. Thus instead of converting acGt (G
2377 * called) to Gact they produce Gtac.
2378 *
2379 * Given 4 columns of data the two schemes produce value orderings as:
2380 *
2381 * Rotate: Acgt Shift: Acgt
2382 * Cgta Cagt
2383 * Gtac Gact
2384 * Tacg Tacg
2385 *
2386 * As a consequence any channel bias for a/c/g/t is better preserved
2387 * by shifting than it is by rotating as each column (except #1) are
2388 * only populated by two base types and 2 of those columns are on a
2389 * 3:1 ratio.
2390 */
2391 /* #define ROTATE_INSTEAD */
tshift(ztr_t * ztr,char * told_c,int tlen,int * new_len)2392 char *tshift(ztr_t *ztr, char *told_c, int tlen, int *new_len) {
2393 int nc, i;
2394 ztr_chunk_t **base = ztr_find_chunks(ztr, ZTR_TYPE_BASE, &nc);
2395 char *bases;
2396 int nbases;
2397 unsigned short *tnew, *told;
2398
2399 if (nc == 0)
2400 return NULL;
2401
2402 if (*told_c != 0)
2403 return NULL; /* assume RAW format trace input */
2404
2405 /* Use last BASE chunk if multiple are present */
2406 /* FIXME: ensure uncompressed first */
2407 bases = base[nc-1]->data+1;
2408 nbases = base[nc-1]->dlength-1;
2409
2410 if (nbases != (tlen-2)/8) {
2411 fprintf(stderr, "Mismatch in number of base calls to samples\n");
2412 return NULL;
2413 }
2414
2415 /* Allocate and initialise header */
2416 told = ((unsigned short *)told_c) + 1;
2417 *new_len = (nbases*4+4) * sizeof(*tnew);
2418 tnew = (unsigned short *)malloc(*new_len);
2419 for (i = 0; i < 4; i++) {
2420 tnew[i] = 0;
2421 }
2422 ((char *)tnew)[0] = ZTR_FORM_TSHIFT;
2423
2424 #ifdef ROTATE_INSTEAD
2425 /* Reorder */
2426 for (i = 0; i < nbases; i++) {
2427 switch(bases[i]) {
2428 case 'T':
2429 tnew[4+i*4+0] = told[3*nbases+i]; /* TACG */
2430 tnew[4+i*4+1] = told[0*nbases+i];
2431 tnew[4+i*4+2] = told[1*nbases+i];
2432 tnew[4+i*4+3] = told[2*nbases+i];
2433 break;
2434 case 'G':
2435 tnew[4+i*4+0] = told[2*nbases+i]; /* GTAC */
2436 tnew[4+i*4+1] = told[3*nbases+i];
2437 tnew[4+i*4+2] = told[0*nbases+i];
2438 tnew[4+i*4+3] = told[1*nbases+i];
2439 break;
2440 case 'C':
2441 tnew[4+i*4+0] = told[1*nbases+i]; /* CGTA */
2442 tnew[4+i*4+1] = told[2*nbases+i];
2443 tnew[4+i*4+2] = told[3*nbases+i];
2444 tnew[4+i*4+3] = told[0*nbases+i];
2445 break;
2446 default:
2447 tnew[4+i*4+0] = told[0*nbases+i]; /* ACGT */
2448 tnew[4+i*4+1] = told[1*nbases+i];
2449 tnew[4+i*4+2] = told[2*nbases+i];
2450 tnew[4+i*4+3] = told[3*nbases+i];
2451 break;
2452 }
2453 }
2454 #else
2455 /* Reorder */
2456 for (i = 0; i < nbases; i++) {
2457 switch(bases[i]) {
2458 case 'T':
2459 tnew[4+i*4+0] = told[3*nbases+i]; /* TACG */
2460 tnew[4+i*4+1] = told[0*nbases+i];
2461 tnew[4+i*4+2] = told[1*nbases+i];
2462 tnew[4+i*4+3] = told[2*nbases+i];
2463 break;
2464 case 'G':
2465 tnew[4+i*4+0] = told[2*nbases+i]; /* GACT */
2466 tnew[4+i*4+1] = told[0*nbases+i];
2467 tnew[4+i*4+2] = told[1*nbases+i];
2468 tnew[4+i*4+3] = told[3*nbases+i];
2469 break;
2470 case 'C':
2471 tnew[4+i*4+0] = told[1*nbases+i]; /* CAGT */
2472 tnew[4+i*4+1] = told[0*nbases+i];
2473 tnew[4+i*4+2] = told[2*nbases+i];
2474 tnew[4+i*4+3] = told[3*nbases+i];
2475 break;
2476 default:
2477 tnew[4+i*4+0] = told[0*nbases+i]; /* ACGT */
2478 tnew[4+i*4+1] = told[1*nbases+i];
2479 tnew[4+i*4+2] = told[2*nbases+i];
2480 tnew[4+i*4+3] = told[3*nbases+i];
2481 break;
2482 }
2483 }
2484 #endif
2485
2486 xfree(base);
2487
2488 return (char *)tnew;
2489 }
2490
untshift(ztr_t * ztr,char * told_c,int tlen,int * new_len)2491 char *untshift(ztr_t *ztr, char *told_c, int tlen, int *new_len) {
2492 unsigned short *tnew, *told = (unsigned short *)told_c;
2493 int nc, nbases, i;
2494 char *bases;
2495 ztr_chunk_t **base = ztr_find_chunks(ztr, ZTR_TYPE_BASE, &nc);
2496
2497 if (nc == 0)
2498 return NULL;
2499
2500 /* Use last BASE chunk if multiple are present */
2501 uncompress_chunk(ztr, base[nc-1]);
2502 bases = base[nc-1]->data+1;
2503 nbases = base[nc-1]->dlength-1;
2504
2505 *new_len = 2 + nbases*4 * sizeof(*tnew);
2506 tnew = (unsigned short *)malloc(*new_len);
2507 tnew[0] = 0;
2508 told += 4;
2509
2510 #ifdef ROTATE_INSTEAD
2511 /* Reorder */
2512 for (i = 0; i < nbases; i++) {
2513 switch(bases[i]) {
2514 case 'T':
2515 tnew[1+3*nbases+i] = *told++; /* TACG */
2516 tnew[1+0*nbases+i] = *told++;
2517 tnew[1+1*nbases+i] = *told++;
2518 tnew[1+2*nbases+i] = *told++;
2519 break;
2520 case 'G':
2521 tnew[1+2*nbases+i] = *told++; /* GTAC */
2522 tnew[1+3*nbases+i] = *told++;
2523 tnew[1+0*nbases+i] = *told++;
2524 tnew[1+1*nbases+i] = *told++;
2525 break;
2526 case 'C':
2527 tnew[1+1*nbases+i] = *told++; /* CGTA */
2528 tnew[1+2*nbases+i] = *told++;
2529 tnew[1+3*nbases+i] = *told++;
2530 tnew[1+0*nbases+i] = *told++;
2531 break;
2532 default:
2533 tnew[1+0*nbases+i] = *told++; /* ACGT */
2534 tnew[1+1*nbases+i] = *told++;
2535 tnew[1+2*nbases+i] = *told++;
2536 tnew[1+3*nbases+i] = *told++;
2537 break;
2538 }
2539 }
2540 #else
2541 /* Reorder */
2542 for (i = 0; i < nbases; i++) {
2543 switch(bases[i]) {
2544 case 'T':
2545 tnew[1+3*nbases+i] = *told++; /* TACG */
2546 tnew[1+0*nbases+i] = *told++;
2547 tnew[1+1*nbases+i] = *told++;
2548 tnew[1+2*nbases+i] = *told++;
2549 break;
2550 case 'G':
2551 tnew[1+2*nbases+i] = *told++; /* GACT */
2552 tnew[1+0*nbases+i] = *told++;
2553 tnew[1+1*nbases+i] = *told++;
2554 tnew[1+3*nbases+i] = *told++;
2555 break;
2556 case 'C':
2557 tnew[1+1*nbases+i] = *told++; /* CAGT */
2558 tnew[1+0*nbases+i] = *told++;
2559 tnew[1+2*nbases+i] = *told++;
2560 tnew[1+3*nbases+i] = *told++;
2561 break;
2562 default:
2563 tnew[1+0*nbases+i] = *told++; /* ACGT */
2564 tnew[1+1*nbases+i] = *told++;
2565 tnew[1+2*nbases+i] = *told++;
2566 tnew[1+3*nbases+i] = *told++;
2567 break;
2568 }
2569 }
2570 #endif
2571
2572 xfree(base);
2573
2574 return (char *)tnew;
2575 }
2576