diff options
Diffstat (limited to 'src/compression/coder.c')
| -rw-r--r-- | src/compression/coder.c | 481 |
1 files changed, 481 insertions, 0 deletions
diff --git a/src/compression/coder.c b/src/compression/coder.c new file mode 100644 index 0000000..bed73ca --- /dev/null +++ b/src/compression/coder.c @@ -0,0 +1,481 @@ +/* This code is part of the tng compression routines. + * + * Written by Daniel Spangberg + * Copyright (c) 2010, 2013, The GROMACS development team. + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + */ + + +#if HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdio.h> +#include <stdlib.h> +#include "tng_compress.h" +#include "bwlzh.h" +#include "coder.h" +#include "warnmalloc.h" + +struct coder *Ptngc_coder_init(void) +{ + struct coder *coder=warnmalloc(sizeof *coder); + coder->pack_temporary_bits=0; + return coder; +} + +void Ptngc_coder_deinit(struct coder *coder) +{ + free(coder); +} + +void Ptngc_out8bits(struct coder *coder, unsigned char **output) +{ + while (coder->pack_temporary_bits>=8) + { + unsigned int mask=~(0xFFU<<(coder->pack_temporary_bits-8)); + unsigned char out=(unsigned char)(coder->pack_temporary>>(coder->pack_temporary_bits-8)); + **output=out; + (*output)++; + coder->pack_temporary_bits-=8; + coder->pack_temporary&=mask; + } +} + +void Ptngc_write_pattern(struct coder *coder,unsigned int pattern, int nbits, unsigned char **output) +{ + unsigned int mask1,mask2; + mask1=1; + mask2=1<<(nbits-1); + coder->pack_temporary<<=nbits; /* Make room for new data. */ + coder->pack_temporary_bits+=nbits; + while (nbits) + { + if (pattern & mask1) + coder->pack_temporary|=mask2; + nbits--; + mask1<<=1; + mask2>>=1; + } + Ptngc_out8bits(coder,output); +} + +/* Write up to 24 bits */ +void Ptngc_writebits(struct coder *coder,unsigned int value,int nbits, unsigned char **output_ptr) +{ + /* Make room for the bits. */ + coder->pack_temporary<<=nbits; + coder->pack_temporary_bits+=nbits; + coder->pack_temporary|=value; + Ptngc_out8bits(coder,output_ptr); +} + +/* Write up to 32 bits */ +void Ptngc_write32bits(struct coder *coder,unsigned int value,int nbits, unsigned char **output_ptr) +{ + unsigned int mask; + if (nbits>=8) + mask=0xFFU<<(nbits-8); + else + mask=0xFFU>>(8-nbits); + while (nbits>8) + { + /* Make room for the bits. */ + coder->pack_temporary<<=8; + coder->pack_temporary_bits+=8; + coder->pack_temporary|=(value&mask)>>(nbits-8); + Ptngc_out8bits(coder,output_ptr); + nbits-=8; + mask>>=8; + } + if (nbits) + Ptngc_writebits(coder,value&mask,nbits,output_ptr); +} + +/* Write "arbitrary" number of bits */ +void Ptngc_writemanybits(struct coder *coder,unsigned char *value,int nbits, unsigned char **output_ptr) +{ + int vptr=0; + while (nbits>=24) + { + unsigned int v=((((unsigned int)value[vptr])<<16)| + (((unsigned int)value[vptr+1])<<8)| + (((unsigned int)value[vptr+2]))); + Ptngc_writebits(coder,v,24,output_ptr); + vptr+=3; + nbits-=24; + } + while (nbits>=8) + { + Ptngc_writebits(coder,(unsigned int)value[vptr],8,output_ptr); + vptr++; + nbits-=8; + } + if (nbits) + { + Ptngc_writebits(coder,(unsigned int)value[vptr],nbits,output_ptr); + } +} + +static int write_stop_bit_code(struct coder *coder, unsigned int s,unsigned int coding_parameter, unsigned char **output) +{ + do { + unsigned int extract=~(0xffffffffU<<coding_parameter); + unsigned int this=(s&extract)<<1; + s>>=coding_parameter; + if (s) + { + this|=1U; + coder->stat_overflow++; + } + coder->pack_temporary<<=(coding_parameter+1); + coder->pack_temporary_bits+=coding_parameter+1; + coder->pack_temporary|=this; + Ptngc_out8bits(coder,output); + if (s) + { + coding_parameter>>=1; + if (coding_parameter<1) + coding_parameter=1; + } + } while (s); + coder->stat_numval++; + return 0; +} + +static int pack_stopbits_item(struct coder *coder,int item, unsigned char **output, int coding_parameter) +{ + /* Find this symbol in table. */ + int s=0; + if (item>0) + s=1+(item-1)*2; + else if (item<0) + s=2+(-item-1)*2; + return write_stop_bit_code(coder,s,coding_parameter,output); + return 0; +} + +static int pack_triplet(struct coder *coder,unsigned int *s, unsigned char **output, int coding_parameter, + unsigned int max_base, int maxbits) +{ + /* Determine base for this triplet. */ + unsigned int min_base=1U<<coding_parameter; + unsigned int this_base=min_base; + int i; + unsigned int jbase=0; + unsigned int bits_per_value; + for (i=0; i<3; i++) + while (s[i]>=this_base) + { + this_base*=2; + jbase++; + } + bits_per_value=coding_parameter+jbase; + if (jbase>=3) + { + if (this_base>max_base) + return 1; + this_base=max_base; + bits_per_value=maxbits; + jbase=3; + } + /* 2 bits selects the base */ + coder->pack_temporary<<=2; + coder->pack_temporary_bits+=2; + coder->pack_temporary|=jbase; + Ptngc_out8bits(coder,output); + for (i=0; i<3; i++) + Ptngc_write32bits(coder,s[i],bits_per_value,output); + return 0; +} + +void Ptngc_pack_flush(struct coder *coder,unsigned char **output) +{ + /* Zero-fill just enough. */ + if (coder->pack_temporary_bits>0) + Ptngc_write_pattern(coder,0,8-coder->pack_temporary_bits,output); +} + +unsigned char *Ptngc_pack_array(struct coder *coder,int *input, int *length, int coding, int coding_parameter,int natoms, int speed) +{ + if ((coding==TNG_COMPRESS_ALGO_BWLZH1) || (coding==TNG_COMPRESS_ALGO_BWLZH2)) + { + unsigned char *output=warnmalloc(4+bwlzh_get_buflen(*length)); + int i,j,k,n=*length; + unsigned int *pval=warnmalloc(n*sizeof *pval); + int nframes=n/natoms/3; + int cnt=0; + int most_negative=2147483647; + for (i=0; i<n; i++) + if (input[i]<most_negative) + most_negative=input[i]; + most_negative=-most_negative; + output[0]=((unsigned int)most_negative)&0xFFU; + output[1]=(((unsigned int)most_negative)>>8)&0xFFU; + output[2]=(((unsigned int)most_negative)>>16)&0xFFU; + output[3]=(((unsigned int)most_negative)>>24)&0xFFU; + for (i=0; i<natoms; i++) + for (j=0; j<3; j++) + for (k=0; k<nframes; k++) + { + int item=input[k*3*natoms+i*3+j]; + pval[cnt++]=(unsigned int)(item+most_negative); + + } + if (speed>=5) + bwlzh_compress(pval,n,output+4,length); + else + bwlzh_compress_no_lz77(pval,n,output+4,length); + (*length)+=4; + free(pval); + return output; + } + else if (coding==TNG_COMPRESS_ALGO_POS_XTC3) + return Ptngc_pack_array_xtc3(input,length,natoms,speed); + else if (coding==TNG_COMPRESS_ALGO_POS_XTC2) + return Ptngc_pack_array_xtc2(coder,input,length); + else + { + unsigned char *output=NULL; + unsigned char *output_ptr=NULL; + int i; + int output_length=0; + + coder->stat_numval=0; + coder->stat_overflow=0; + /* Allocate enough memory for output */ + output=warnmalloc(8* *length*sizeof *output); + output_ptr=output; + if ((coding==TNG_COMPRESS_ALGO_TRIPLET) || + (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) || + (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE)) + { + /* Pack triplets. */ + int ntriplets=*length/3; + /* Determine max base and maxbits */ + unsigned int max_base=1U<<coding_parameter; + unsigned int maxbits=coding_parameter; + unsigned int intmax=0; + for (i=0; i<*length; i++) + { + int item=input[i]; + unsigned int s=0; + if (item>0) + s=1+(item-1)*2; + else if (item<0) + s=2+(-item-1)*2; + if (s>intmax) + intmax=s; + } + /* Store intmax */ + coder->pack_temporary_bits=32; + coder->pack_temporary=intmax; + Ptngc_out8bits(coder,&output_ptr); + while (intmax>=max_base) + { + max_base*=2; + maxbits++; + } + for (i=0; i<ntriplets; i++) + { + int j; + unsigned int s[3]; + for (j=0; j<3; j++) + { + int item=input[i*3+j]; + /* Find this symbol in table. */ + s[j]=0; + if (item>0) + s[j]=1+(item-1)*2; + else if (item<0) + s[j]=2+(-item-1)*2; + } + if (pack_triplet(coder,s,&output_ptr,coding_parameter,max_base,maxbits)) + { + free(output); + return NULL; + } + } + } + else + for (i=0; i<*length; i++) + if (pack_stopbits_item(coder,input[i],&output_ptr,coding_parameter)) + { + free(output); + return NULL; + } + Ptngc_pack_flush(coder,&output_ptr); + output_length=(int)(output_ptr-output); + *length=output_length; + return output; + } +} + +static int unpack_array_stop_bits(struct coder *coder,unsigned char *packed,int *output, int length, int coding_parameter) +{ + int i,j; + unsigned int extract_mask=0x80; + unsigned char *ptr=packed; + for (i=0; i<length; i++) + { + unsigned int pattern=0; + int numbits=coding_parameter; + unsigned int bit; + int s; + unsigned int insert_mask=1U<<(numbits-1); + int inserted_bits=numbits; + do { + for (j=0; j<numbits; j++) + { + bit=*ptr & extract_mask; + if (bit) + pattern|=insert_mask; + insert_mask>>=1; + extract_mask>>=1; + if (!extract_mask) + { + extract_mask=0x80; + ptr++; + } + } + /* Check stop bit */ + bit=*ptr & extract_mask; + extract_mask>>=1; + if (!extract_mask) + { + extract_mask=0x80; + ptr++; + } + if (bit) + { + numbits>>=1; + if (numbits<1) + numbits=1; + inserted_bits+=numbits; + insert_mask=1U<<(inserted_bits-1); + } + } while (bit); + s=(pattern+1)/2; + if ((pattern%2)==0) + s=-s; + output[i]=s; + } + return 0; +} + +static int unpack_array_triplet(struct coder *coder,unsigned char *packed,int *output, int length, int coding_parameter) +{ + int i,j; + unsigned int extract_mask=0x80; + unsigned char *ptr=packed; + /* Determine max base and maxbits */ + unsigned int max_base=1U<<coding_parameter; + unsigned int maxbits=coding_parameter; + unsigned int intmax; + /* Get intmax */ + intmax=((unsigned int)ptr[0])<<24| + ((unsigned int)ptr[1])<<16| + ((unsigned int)ptr[2])<<8| + ((unsigned int)ptr[3]); + ptr+=4; + while (intmax>=max_base) + { + max_base*=2; + maxbits++; + } + length/=3; + for (i=0; i<length; i++) + { + /* Find base */ + unsigned int jbase=0; + unsigned int numbits; + unsigned int bit; + for (j=0; j<2; j++) + { + bit=*ptr & extract_mask; + jbase<<=1; + if (bit) + jbase|=1U; + extract_mask>>=1; + if (!extract_mask) + { + extract_mask=0x80; + ptr++; + } + } + if (jbase==3) + numbits=maxbits; + else + numbits=coding_parameter+jbase; + for (j=0; j<3; j++) + { + int s; + unsigned int jbit; + unsigned int pattern=0; + for (jbit=0; jbit<numbits; jbit++) + { + bit=*ptr & extract_mask; + pattern<<=1; + if (bit) + pattern|=1U; + extract_mask>>=1; + if (!extract_mask) + { + extract_mask=0x80; + ptr++; + } + } + s=(pattern+1)/2; + if ((pattern%2)==0) + s=-s; + output[i*3+j]=s; + } + } + return 0; +} + +static int unpack_array_bwlzh(struct coder *coder,unsigned char *packed,int *output, int length, int natoms) +{ + int i,j,k,n=length; + unsigned int *pval=warnmalloc(n*sizeof *pval); + int nframes=n/natoms/3; + int cnt=0; + int most_negative=(int)(((unsigned int)packed[0]) | + (((unsigned int)packed[1])<<8) | + (((unsigned int)packed[2])<<16) | + (((unsigned int)packed[3])<<24)); + bwlzh_decompress(packed+4,length,pval); + for (i=0; i<natoms; i++) + for (j=0; j<3; j++) + for (k=0; k<nframes; k++) + { + unsigned int s=pval[cnt++]; + output[k*3*natoms+i*3+j]=(int)s-most_negative; + } + free(pval); + return 0; +} + +int Ptngc_unpack_array(struct coder *coder,unsigned char *packed,int *output, int length, int coding, int coding_parameter, int natoms) +{ + if ((coding==TNG_COMPRESS_ALGO_STOPBIT) || + (coding==TNG_COMPRESS_ALGO_VEL_STOPBIT_INTER)) + return unpack_array_stop_bits(coder, packed, output, length, coding_parameter); + else if ((coding==TNG_COMPRESS_ALGO_TRIPLET) || + (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_INTRA) || + (coding==TNG_COMPRESS_ALGO_POS_TRIPLET_ONETOONE)) + return unpack_array_triplet(coder, packed, output, length, coding_parameter); + else if (coding==TNG_COMPRESS_ALGO_POS_XTC2) + return Ptngc_unpack_array_xtc2(coder, packed, output, length); + else if ((coding==TNG_COMPRESS_ALGO_BWLZH1) || (coding==TNG_COMPRESS_ALGO_BWLZH2)) + return unpack_array_bwlzh(coder, packed, output, length,natoms); + else if (coding==TNG_COMPRESS_ALGO_POS_XTC3) + return Ptngc_unpack_array_xtc3(packed, output, length,natoms); + return 1; +} + |