Trisurf Monte Carlo simulator
Samo Penic
2020-04-09 4f5ffcfd5ba38b6b7dd6d3b15de8bb8677537b9f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
#include<stdio.h>
#include<stdlib.h>
#include<general.h>
#include<string.h>
#include<zlib.h>
#include<inttypes.h>
#include<b64zlib_compression.h>
 
 
/* zlib compression base64 encoded */
/* compressed must not be pre-malloced */
/* taken from https://gist.github.com/arq5x/5315739 */
/* This function is not to be used in compressing paraview data. See below. */
ts_uint ts_compress_data(char *data, ts_uint data_len, char **compressed){
    z_stream defstream;
    defstream.zalloc = Z_NULL;
    defstream.zfree = Z_NULL;
    defstream.opaque = Z_NULL;
    defstream.avail_in = data_len+1;
    defstream.next_in = (unsigned char *)data;    
    char *compr=(char *)malloc(data_len*sizeof(char));
    defstream.avail_out = data_len+1;
    defstream.next_out = (unsigned char *)compr;
    deflateInit(&defstream, 6);
//    deflateInit(&defstream, Z_BEST_COMPRESSION);
        deflate(&defstream, Z_FINISH);
        deflateEnd(&defstream);
    *compressed=compr;
    return defstream.total_out;
}
 
/* This function may be used to compress lange strings, however it is not compatible with paraview */
ts_uint ts_compress_string64(char *data, ts_uint data_len, char **encoded_compressed){
    size_t nbase;
    char *compr;
    size_t number_of_compressed_bytes=ts_compress_data(data, data_len, &compr);
    *encoded_compressed=base64_encode((unsigned char *)compr,number_of_compressed_bytes,&nbase);
    free(compr);
    return nbase;
}
 
 
/** @brief Compresses and base64 encodes a specific section in VTK file. Function corresponds to rules prescribed by VTK standards.
 *
 *  The list of Int8, Int64 and Double precision values (or vectors) are compressed and base64 encoded according to the follow sheme:
 *
 * Binary data in the file are compressed when the VTKFile element is of the form
 * 
 * <VTKFile ... compressor="vtkZLibDataCompressor">
 *
 * The data corresponding to a data array are stored in a set of blocks which are each compressed using the zlib library. The block structure allows semi-random access without decompressing all data. In uncompressed form all the blocks have the same size except possibly the last block which might be smaller. The data for one array begin with a header of the form
 *
 * [#blocks][#u-size][#p-size][#c-size-1][#c-size-2]...[#c-size-#blocks][DATA]
 *
 * Each token is an integer value whose type is specified by "header_type" at the top of the file (UInt32 if no type specified). The token meanings are:
 *
 * [#blocks] = Number of blocks
 * [#u-size] = Block size before compression
 * [#p-size] = Size of last partial block (zero if it not needed)
 * [#c-size-i] = Size in bytes of block i after compression
 *
 * The [DATA] portion stores contiguously every block appended together. The offset from the beginning of the data section to the beginning of a block is computed by summing the compressed block sizes from preceding blocks according to the header.
 *
 *
 *  @param *data is a pointer of bytes to original vector containting double, long or char values.
 *  @param *data_len is a number of bytes in *data pointer, NOT the number of containing values!
 *  @returns a pointer to string containing base64 encoded compressed string. Should be freed somewhere afterwards.
*/
char *ts_compress(char *data, ts_uint data_len){
    size_t nbase1, nbase2;
    unsigned char *compr=(unsigned char *)malloc(data_len);
    size_t number_of_compressed_bytes=data_len;
 
    compress(compr,&number_of_compressed_bytes, (unsigned char *)data, data_len);
 
    char *encoded_compressed=base64_encode((unsigned char *)compr,number_of_compressed_bytes,&nbase1);
 
    free(compr);
 
    ts_uint header[4]={1, data_len, data_len, number_of_compressed_bytes};
    char *encoded_header=(char *)base64_encode((unsigned char *)header, 4*sizeof(ts_uint), &nbase2);
    char *return_value=malloc((nbase1+nbase2+1)*sizeof(char));
    strncpy(return_value,encoded_header,nbase2);
    strncpy(return_value+nbase2,encoded_compressed,nbase1);
    *(return_value+nbase1+nbase2)=0;
 
    free(encoded_compressed);
    free(encoded_header);
 
    //test decoding and decompression
    //char *test_data = ts_decompress((unsigned char *)return_value, nbase1+nbase2, data_len);
    //for(int i=0;i<data_len;i++){
    //    if(data[i]!=test_data[i]){
    //        fprintf(stderr,"(DE)COMPRESSION ERROR!\n");
    //    }
    //}
    //free(test_data);
    return return_value;
}
 
char  *ts_decompress(unsigned char *compressed_data, unsigned int data_len, size_t *decompressed_length){
    const size_t encoded_header_len=24;
    ts_uint *header;
    size_t nbase1=0,nbase2=0;
    header=(ts_uint *)base64_decode((const char *)compressed_data, encoded_header_len, &nbase2);
    if(header==NULL) fprintf(stderr, "Error decoding header\n");
    //fprintf(stderr,"Header=%d, %d, %d, %d\n", header[0],header[1], header[2], header[3]);
    unsigned long result_len=header[2];
    unsigned char *decoded_data=base64_decode((const char *)&compressed_data[encoded_header_len], data_len-encoded_header_len, &nbase1);
    if(decoded_data==NULL) fprintf(stderr, "Error decoding compressed data\n");
    unsigned char *return_value=(unsigned char *)malloc(result_len*sizeof(char));
    uncompress(return_value, &result_len, decoded_data, header[3]);
    *decompressed_length=result_len;
    free(decoded_data);
    free(header);
    return (char *)return_value;
}
 
ts_uint ts_decompress_string64(char *b64, ts_uint data_len, char **decompressed){
return TS_SUCCESS;
}
 
/* base64 encoding, taken from http://stackoverflow.com/questions/342409/how-do-i-base64-encode-decode-in-c */
static char encoding_table[] = {'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H',
                                'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P',
                                'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X',
                                'Y', 'Z', 'a', 'b', 'c', 'd', 'e', 'f',
                                'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n',
                                'o', 'p', 'q', 'r', 's', 't', 'u', 'v',
                                'w', 'x', 'y', 'z', '0', '1', '2', '3',
                                '4', '5', '6', '7', '8', '9', '+', '/'};
static int mod_table[] = {0, 2, 1};
 
 
char *base64_encode(const unsigned char *data, size_t input_length, size_t *output_length) {
    *output_length = 4 * ((input_length + 2) / 3);
    int i,j;
    char *encoded_data = malloc(*output_length);
    if (encoded_data == NULL) return NULL;
    
    for (i = 0, j = 0; i < input_length;) {
 
           uint32_t octet_a = i < input_length ? (unsigned char)data[i++] : 0;
        uint32_t octet_b = i < input_length ? (unsigned char)data[i++] : 0;
           uint32_t octet_c = i < input_length ? (unsigned char)data[i++] : 0;
        uint32_t triple = (octet_a << 0x10) + (octet_b << 0x08) + octet_c;
 
           encoded_data[j++] = encoding_table[(triple >> 3 * 6) & 0x3F];
           encoded_data[j++] = encoding_table[(triple >> 2 * 6) & 0x3F];
           encoded_data[j++] = encoding_table[(triple >> 1 * 6) & 0x3F];
           encoded_data[j++] = encoding_table[(triple >> 0 * 6) & 0x3F];
    }
 
    for (i = 0; i < mod_table[input_length % 3]; i++)
           encoded_data[*output_length - 1 - i] = '=';
 
    return encoded_data;
}
 
 
unsigned char *base64_decode(const char *data, size_t input_length, size_t *output_length) {
    int i,j;
    char *decoding_table = build_decoding_table();
 
    if (input_length % 4 != 0) return NULL;
 
    *output_length = input_length / 4 * 3;
    if (data[input_length - 1] == '=') (*output_length)--;
    if (data[input_length - 2] == '=') (*output_length)--;
 
    unsigned char *decoded_data = malloc(*output_length);
    if (decoded_data == NULL) return NULL;
 
    for (i = 0, j = 0; i < input_length;) {
           uint32_t sextet_a = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
           uint32_t sextet_b = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
           uint32_t sextet_c = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
           uint32_t sextet_d = data[i] == '=' ? 0 & i++ : decoding_table[(int)data[i++]];
 
           uint32_t triple = (sextet_a << 3 * 6)
            + (sextet_b << 2 * 6)
            + (sextet_c << 1 * 6)
            + (sextet_d << 0 * 6);
 
           if (j < *output_length) decoded_data[j++] = (triple >> 2 * 8) & 0xFF;
           if (j < *output_length) decoded_data[j++] = (triple >> 1 * 8) & 0xFF;
           if (j < *output_length) decoded_data[j++] = (triple >> 0 * 8) & 0xFF;
    }
    free(decoding_table);
    return decoded_data;
}
 
 
char *build_decoding_table() {
    char *decoding_table = malloc(256);
    int i;
    for (i = 0; i < 64; i++)
        decoding_table[(unsigned char) encoding_table[i]] = i;
    return decoding_table;
}