Line data Source code
1 : /* libblx - Magellan BLX topo reader/writer library
2 : *
3 : * Copyright (c) 2008, Henrik Johansson <henrik@johome.net>
4 : * Copyright (c) 2008-2009, Even Rouault <even dot rouault at spatialys.com>
5 : *
6 : * SPDX-License-Identifier: MIT
7 : */
8 :
9 : #include "blx.h"
10 : #include <stdio.h>
11 : #include <string.h>
12 : #include <stdlib.h>
13 :
14 : #include "cpl_port.h"
15 :
16 : /* Constants */
17 : #define MAXLEVELS 5
18 : #define MAXCOMPONENTS 4
19 :
20 : static const int table1[] = {
21 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
22 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
23 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
24 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
25 : 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
26 : 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
27 : 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2,
28 : 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
29 : 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3,
30 : 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4,
31 : 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5,
32 : 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7,
33 : 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10,
34 : 10, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14, 14,
35 : 15, 15, 16, 16, 17, 17, 18, 18, 19, 19, 20, 21, 22, 23, 24,
36 : 25, 26, 27, 28, 29, 30, 31, 255, 255, 255, 255, 255, 255, 255, 255,
37 : 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
38 : 255};
39 :
40 : /* { byte, n of bits when compressed, bit pattern << (13-n of bits) } */
41 : static const int table2[][3] = {
42 : {0, 2, 0}, {255, 3, 2048}, {1, 3, 3072}, {2, 4, 4096},
43 : {3, 4, 4608}, {254, 5, 5120}, {4, 5, 5376}, {5, 5, 5632},
44 : {253, 6, 5888}, {6, 6, 6016}, {252, 6, 6144}, {7, 6, 6272},
45 : {251, 6, 6400}, {8, 6, 6528}, {9, 7, 6656}, {250, 7, 6720},
46 : {10, 7, 6784}, {249, 7, 6848}, {11, 7, 6912}, {248, 7, 6976},
47 : {12, 8, 7040}, {247, 8, 7072}, {16, 8, 7104}, {246, 8, 7136},
48 : {13, 8, 7168}, {245, 8, 7200}, {14, 8, 7232}, {244, 8, 7264},
49 : {15, 8, 7296}, {243, 8, 7328}, {242, 8, 7360}, {241, 8, 7392},
50 : {17, 9, 7424}, {18, 9, 7440}, {240, 9, 7456}, {239, 9, 7472},
51 : {19, 9, 7488}, {238, 9, 7504}, {20, 9, 7520}, {237, 9, 7536},
52 : {21, 9, 7552}, {236, 9, 7568}, {22, 9, 7584}, {235, 9, 7600},
53 : {234, 9, 7616}, {23, 9, 7632}, {233, 9, 7648}, {24, 10, 7664},
54 : {232, 10, 7672}, {231, 10, 7680}, {25, 10, 7688}, {230, 10, 7696},
55 : {229, 10, 7704}, {26, 10, 7712}, {228, 10, 7720}, {27, 10, 7728},
56 : {227, 10, 7736}, {225, 10, 7744}, {226, 10, 7752}, {28, 10, 7760},
57 : {29, 10, 7768}, {224, 10, 7776}, {30, 10, 7784}, {31, 10, 7792},
58 : {223, 10, 7800}, {32, 10, 7808}, {222, 10, 7816}, {33, 10, 7824},
59 : {221, 11, 7832}, {220, 11, 7836}, {34, 11, 7840}, {219, 11, 7844},
60 : {35, 11, 7848}, {218, 11, 7852}, {256, 11, 7856}, {36, 11, 7860},
61 : {217, 11, 7864}, {216, 11, 7868}, {37, 11, 7872}, {215, 11, 7876},
62 : {38, 11, 7880}, {214, 11, 7884}, {193, 11, 7888}, {213, 11, 7892},
63 : {39, 11, 7896}, {128, 11, 7900}, {212, 11, 7904}, {40, 11, 7908},
64 : {194, 11, 7912}, {211, 11, 7916}, {210, 11, 7920}, {41, 11, 7924},
65 : {209, 11, 7928}, {208, 11, 7932}, {42, 11, 7936}, {207, 11, 7940},
66 : {43, 11, 7944}, {195, 11, 7948}, {206, 11, 7952}, {205, 11, 7956},
67 : {204, 11, 7960}, {44, 11, 7964}, {203, 11, 7968}, {192, 11, 7972},
68 : {196, 11, 7976}, {45, 11, 7980}, {201, 11, 7984}, {200, 11, 7988},
69 : {197, 11, 7992}, {202, 11, 7996}, {127, 11, 8000}, {199, 11, 8004},
70 : {198, 11, 8008}, {46, 12, 8012}, {47, 12, 8014}, {48, 12, 8016},
71 : {49, 12, 8018}, {50, 12, 8020}, {51, 12, 8022}, {191, 12, 8024},
72 : {52, 12, 8026}, {183, 12, 8028}, {53, 12, 8030}, {54, 12, 8032},
73 : {55, 12, 8034}, {190, 12, 8036}, {56, 12, 8038}, {57, 12, 8040},
74 : {189, 12, 8042}, {58, 12, 8044}, {176, 12, 8046}, {59, 12, 8048},
75 : {126, 12, 8050}, {60, 12, 8052}, {188, 12, 8054}, {61, 12, 8056},
76 : {63, 12, 8058}, {62, 12, 8060}, {64, 12, 8062}, {129, 12, 8064},
77 : {187, 12, 8066}, {186, 12, 8068}, {65, 12, 8070}, {66, 12, 8072},
78 : {185, 12, 8074}, {184, 12, 8076}, {68, 12, 8078}, {174, 12, 8080},
79 : {67, 12, 8082}, {182, 13, 8084}, {69, 13, 8085}, {180, 13, 8086},
80 : {181, 13, 8087}, {71, 13, 8088}, {70, 13, 8089}, {179, 13, 8090},
81 : {125, 13, 8091}, {72, 13, 8092}, {130, 13, 8093}, {178, 13, 8094},
82 : {177, 13, 8095}, {73, 13, 8096}, {74, 13, 8097}, {124, 13, 8098},
83 : {76, 13, 8099}, {175, 13, 8100}, {75, 13, 8101}, {131, 13, 8102},
84 : {132, 13, 8103}, {79, 13, 8104}, {77, 13, 8105}, {123, 13, 8106},
85 : {80, 13, 8107}, {172, 13, 8108}, {171, 13, 8109}, {78, 13, 8110},
86 : {173, 13, 8111}, {81, 13, 8112}, {169, 13, 8113}, {122, 13, 8114},
87 : {82, 13, 8115}, {133, 13, 8116}, {168, 13, 8117}, {84, 13, 8118},
88 : {164, 13, 8119}, {167, 13, 8120}, {85, 13, 8121}, {170, 13, 8122},
89 : {166, 13, 8123}, {165, 13, 8124}, {121, 13, 8125}, {160, 13, 8126},
90 : {134, 13, 8127}, {136, 13, 8128}, {161, 13, 8129}, {120, 13, 8130},
91 : {88, 13, 8131}, {83, 13, 8132}, {119, 13, 8133}, {163, 13, 8134},
92 : {162, 13, 8135}, {159, 13, 8136}, {91, 13, 8137}, {135, 13, 8138},
93 : {90, 13, 8139}, {86, 13, 8140}, {137, 13, 8141}, {87, 13, 8142},
94 : {89, 13, 8143}, {158, 13, 8144}, {152, 13, 8145}, {138, 13, 8146},
95 : {139, 13, 8147}, {116, 13, 8148}, {140, 13, 8149}, {92, 13, 8150},
96 : {96, 13, 8151}, {157, 13, 8152}, {153, 13, 8153}, {97, 13, 8154},
97 : {94, 13, 8155}, {93, 13, 8156}, {117, 13, 8157}, {156, 13, 8158},
98 : {155, 13, 8159}, {95, 13, 8160}, {118, 13, 8161}, {143, 13, 8162},
99 : {151, 13, 8163}, {142, 13, 8164}, {104, 13, 8165}, {100, 13, 8166},
100 : {148, 13, 8167}, {144, 13, 8168}, {154, 13, 8169}, {115, 13, 8170},
101 : {113, 13, 8171}, {98, 13, 8172}, {146, 13, 8173}, {112, 13, 8174},
102 : {145, 13, 8175}, {149, 13, 8176}, {141, 13, 8177}, {150, 13, 8178},
103 : {103, 13, 8179}, {147, 13, 8180}, {99, 13, 8181}, {108, 13, 8182},
104 : {101, 13, 8183}, {114, 13, 8184}, {105, 13, 8185}, {102, 13, 8186},
105 : {107, 13, 8187}, {109, 13, 8188}, {110, 13, 8189}, {111, 13, 8190},
106 : {106, 13, 8191}, {0, 0, 8192}};
107 :
108 : static const int table3[] = {0x20, 0x2f, 0x44, 0x71, 0x95, 0x101};
109 :
110 32 : STATIC int compress_chunk(unsigned char *inbuf, int inlen,
111 : unsigned char *outbuf, int outbuflen)
112 : {
113 32 : int next, m = 0, j, outlen = 0;
114 32 : unsigned reg = 0;
115 :
116 32 : next = *inbuf++;
117 32 : inlen--;
118 541230 : while (next >= 0)
119 : {
120 : /* Find index of input byte in table2 and put it in j */
121 541198 : j = 0;
122 7158480 : while (next != table2[j][0])
123 6617280 : j++;
124 :
125 541198 : if (inlen)
126 : {
127 541134 : next = *inbuf++;
128 541134 : inlen--;
129 : }
130 : else
131 : {
132 64 : if (next == 0x100)
133 32 : next = -1;
134 : else
135 32 : next = 0x100;
136 : }
137 541198 : reg = (reg << table2[j][1]) | (table2[j][2] >> (13 - table2[j][1]));
138 541198 : m += table2[j][1];
139 :
140 893230 : while (m >= 8)
141 : {
142 352032 : if (outlen >= outbuflen)
143 0 : return -1;
144 352032 : *outbuf++ = (unsigned char)((reg >> (m - 8)) & 0xff);
145 352032 : outlen++;
146 352032 : m -= 8;
147 : }
148 : }
149 32 : if (outlen >= outbuflen)
150 0 : return -1;
151 32 : *outbuf++ = (unsigned char)((reg << (8 - m)) & 0xff);
152 :
153 32 : return outlen + 1;
154 : }
155 :
156 192 : STATIC int uncompress_chunk(unsigned char *inbuf, int inlen,
157 : unsigned char *outbuf, int outbuflen)
158 : {
159 192 : int i, j, k, m = 0, outlen = 0;
160 : unsigned reg, newdata;
161 :
162 192 : if (inlen < 4)
163 0 : return -1;
164 :
165 192 : reg = *(inbuf + 3) | (*(inbuf + 2) << 8) | (*(inbuf + 1) << 16) |
166 192 : ((unsigned)*(inbuf + 0) << 24);
167 192 : inbuf += 4;
168 192 : inlen -= 4;
169 :
170 192 : newdata = (reg >> 19) & 0x1fff;
171 :
172 : while (1)
173 : {
174 2579890 : j = newdata >> 5;
175 :
176 2579890 : if (table1[j] == 0xff)
177 : {
178 114604 : i = 1;
179 203816 : while ((int)newdata >= table2[table3[i]][2])
180 89212 : i++;
181 :
182 114604 : j = table3[i - 1];
183 :
184 114604 : k = j + ((newdata - table2[j][2]) >> (13 - table2[j][1]));
185 :
186 114604 : if (table2[k][0] == 0x100)
187 192 : break;
188 : else
189 : {
190 114412 : if (outlen >= outbuflen)
191 0 : return -1;
192 114412 : *outbuf++ = (unsigned char)table2[k][0];
193 114412 : outlen++;
194 : }
195 : }
196 : else
197 : {
198 2465290 : j = table1[j];
199 2465290 : if (outlen >= outbuflen)
200 0 : return -1;
201 2465290 : *outbuf++ = (unsigned char)table2[j][0];
202 2465290 : outlen++;
203 : }
204 :
205 2579700 : m += table2[j][1];
206 :
207 2579700 : if (m >= 19)
208 : {
209 616582 : if (m >= 8)
210 : {
211 1929530 : for (i = m >> 3; i; i--)
212 : {
213 1312940 : if (inlen)
214 : {
215 1312750 : reg = (reg << 8) | *inbuf++;
216 1312750 : inlen--;
217 : }
218 : else
219 192 : reg = reg << 8;
220 : }
221 : }
222 616582 : m = m & 7;
223 : }
224 2579700 : newdata = (reg >> (19 - m)) & 0x1fff;
225 : }
226 192 : return outlen;
227 : }
228 :
229 : /*
230 : Reconstruct a new detail level with double resolution in the horizontal
231 : direction from data from the previous detail level and plus new differential
232 : data.
233 : */
234 1600 : STATIC blxdata *reconstruct_horiz(blxdata *base, blxdata *diff, unsigned rows,
235 : unsigned cols, blxdata *out)
236 : {
237 : unsigned int i, j;
238 : blxdata tmp;
239 :
240 : /* Last column */
241 36672 : for (i = 0; i < rows; i++)
242 35072 : out[2 * (cols * i + cols - 1)] =
243 35072 : diff[cols * i + cols - 1] +
244 35072 : (((short)(base[cols * i + cols - 2] - base[cols * i + cols - 1] -
245 : 1)) >>
246 : 2);
247 :
248 : /* Intermediate columns */
249 36672 : for (i = 0; i < rows; i++)
250 1419010 : for (j = cols - 2; j > 0; j--)
251 1383940 : out[2 * (cols * i + j)] =
252 1383940 : diff[cols * i + j] +
253 1383940 : (((short)(base[cols * i + j] +
254 1383940 : 2 * (base[cols * i + j - 1] -
255 1383940 : out[2 * (cols * i + j + 1)]) -
256 1383940 : 3 * base[cols * i + j + 1] + 1)) >>
257 : 3);
258 :
259 : /* First column */
260 36672 : for (i = 0; i < rows; i++)
261 35072 : out[2 * cols * i] =
262 35072 : diff[cols * i] +
263 35072 : (((short)(base[cols * i] - base[cols * i + 1] + 1)) >> 2);
264 :
265 36672 : for (i = 0; i < rows; i++)
266 1489150 : for (j = 0; j < cols; j++)
267 : {
268 1454080 : tmp = base[cols * i + j] +
269 1454080 : (((short)(out[2 * (cols * i + j)] + 1)) >> 1);
270 1454080 : out[2 * cols * i + 2 * j + 1] = tmp - out[2 * (cols * i + j)];
271 1454080 : out[2 * cols * i + 2 * j] = tmp;
272 : }
273 :
274 1600 : return out;
275 : }
276 :
277 : /*
278 : Reconstruct a new detail level with double resolution in the vertical
279 : direction from data from the previous detail level and plus new differential
280 : data.
281 : */
282 800 : STATIC blxdata *reconstruct_vert(blxdata *base, blxdata *diff, unsigned rows,
283 : unsigned cols, blxdata *out)
284 : {
285 : unsigned int i, j;
286 : blxdata tmp;
287 :
288 : /* Last row */
289 35872 : for (i = 0; i < cols; i++)
290 35072 : out[2 * cols * (rows - 1) + i] =
291 35072 : diff[cols * (rows - 1) + i] +
292 35072 : (((short)(base[cols * (rows - 2) + i] -
293 35072 : base[cols * (rows - 1) + i] - 1)) >>
294 : 2);
295 :
296 : /* Intermediate rows */
297 35872 : for (i = 0; i < cols; i++)
298 1419010 : for (j = rows - 2; j > 0; j--)
299 1383940 : out[2 * cols * j + i] =
300 1383940 : diff[cols * j + i] +
301 1383940 : ((short)((base[cols * j + i] +
302 1383940 : 2 * (base[cols * (j - 1) + i] -
303 1383940 : out[2 * cols * (j + 1) + i]) -
304 1383940 : 3 * base[cols * (j + 1) + i] + 1)) >>
305 : 3);
306 :
307 : /* First row */
308 35872 : for (i = 0; i < cols; i++)
309 35072 : out[i] = diff[i] + (((short)(base[i] - base[i + cols] + 1)) >> 2);
310 :
311 35872 : for (i = 0; i < cols; i++)
312 1489150 : for (j = 0; j < rows; j++)
313 : {
314 1454080 : tmp = base[cols * j + i] +
315 1454080 : (((short)(out[2 * cols * j + i] + 1)) >> 1);
316 1454080 : out[cols * (2 * j + 1) + i] = tmp - out[2 * cols * j + i];
317 1454080 : out[cols * 2 * j + i] = tmp;
318 : }
319 800 : return out;
320 : }
321 :
322 : /*
323 : Inverse of reconstruct_horiz
324 : */
325 320 : STATIC void decimate_horiz(blxdata *in, unsigned int rows, unsigned int cols,
326 : blxdata *base, blxdata *diff)
327 : {
328 : unsigned int i, j;
329 : blxdata tmp;
330 :
331 8256 : for (i = 0; i < rows; i++)
332 : {
333 357120 : for (j = 0; j < cols; j += 2)
334 : {
335 349184 : tmp = in[i * cols + j] - in[i * cols + j + 1];
336 349184 : diff[i * cols / 2 + j / 2] = tmp;
337 349184 : base[i * cols / 2 + j / 2] =
338 349184 : in[i * cols + j] - (((short)(tmp + 1)) >> 1);
339 : }
340 : }
341 :
342 : /* First column */
343 8256 : for (i = 0; i < rows; i++)
344 : {
345 7936 : diff[cols / 2 * i] -=
346 7936 : ((short)(base[i * cols / 2] - base[i * cols / 2 + 1] + 1)) >> 2;
347 : }
348 :
349 : /* Intermediate columns */
350 8256 : for (i = 0; i < rows; i++)
351 341248 : for (j = 1; j < cols / 2 - 1; j++)
352 333312 : diff[cols / 2 * i + j] -=
353 333312 : ((short)(base[cols / 2 * i + j] +
354 333312 : 2 * (base[cols / 2 * i + j - 1] -
355 333312 : diff[cols / 2 * i + j + 1]) -
356 333312 : 3 * base[cols / 2 * i + j + 1] + 1)) >>
357 : 3;
358 :
359 : /* Last column */
360 8256 : for (i = 0; i < rows; i++)
361 7936 : diff[cols / 2 * i + cols / 2 - 1] -=
362 7936 : ((short)(base[i * cols / 2 + cols / 2 - 2] -
363 7936 : base[i * cols / 2 + cols / 2 - 1] - 1)) >>
364 : 2;
365 320 : }
366 :
367 : /*
368 : Inverse of reconstruct_vert
369 : */
370 160 : STATIC void decimate_vert(blxdata *in, unsigned int rows, unsigned int cols,
371 : blxdata *base, blxdata *diff)
372 : {
373 : unsigned int i, j;
374 : blxdata tmp;
375 :
376 4128 : for (i = 0; i < rows; i += 2)
377 353152 : for (j = 0; j < cols; j++)
378 : {
379 349184 : tmp = in[i * cols + j] - in[(i + 1) * cols + j];
380 349184 : diff[i / 2 * cols + j] = tmp;
381 349184 : base[i / 2 * cols + j] =
382 349184 : in[i * cols + j] - (((short)(tmp + 1)) >> 1);
383 : }
384 :
385 : /* First row */
386 8096 : for (j = 0; j < cols; j++)
387 7936 : diff[j] -= ((short)(base[j] - base[cols + j] + 1)) >> 2;
388 :
389 : /* Intermediate rows */
390 3808 : for (i = 1; i < rows / 2 - 1; i++)
391 336960 : for (j = 0; j < cols; j++)
392 333312 : diff[cols * i + j] -= ((short)(base[cols * i + j] +
393 333312 : 2 * (base[cols * (i - 1) + j] -
394 333312 : diff[cols * (i + 1) + j]) -
395 333312 : 3 * base[cols * (i + 1) + j] + 1)) >>
396 : 3;
397 :
398 : /* Last row */
399 8096 : for (j = 0; j < cols; j++)
400 7936 : diff[cols * (rows / 2 - 1) + j] -=
401 7936 : ((short)(base[cols * (rows / 2 - 2) + j] -
402 7936 : base[cols * (rows / 2 - 1) + j] - 1)) >>
403 : 2;
404 160 : }
405 :
406 : typedef union
407 : {
408 : short s;
409 : unsigned short u;
410 : } unionshort;
411 :
412 57615 : static int get_short_le(const unsigned char **data)
413 : {
414 : /* We assume two's complement representation for this to work */
415 57615 : unionshort result = {0};
416 57615 : result.u = (unsigned short)(*(*data) | (*(*data + 1) << 8));
417 57615 : *data += 2;
418 57615 : return result.s;
419 : }
420 :
421 9270 : static int get_short_be(const unsigned char **data)
422 : {
423 : /* We assume two's complement representation for this to work */
424 9270 : unionshort result = {0};
425 9270 : result.u = (unsigned short)(*(*data + 1) | (*(*data) << 8));
426 9270 : *data += 2;
427 9270 : return result.s;
428 : }
429 :
430 17644 : static void put_short_le(short data, unsigned char **bufptr)
431 : {
432 : /* We assume two's complement representation for this to work */
433 17644 : unionshort us = {0};
434 17644 : us.s = data;
435 17644 : *(*bufptr)++ = (unsigned char)(us.u & 0xff);
436 17644 : *(*bufptr)++ = (unsigned char)((us.u >> 8) & 0xff);
437 17644 : }
438 :
439 3602 : static void put_short_be(short data, unsigned char **bufptr)
440 : {
441 : /* We assume two's complement representation for this to work */
442 3602 : unionshort us = {0};
443 3602 : us.s = data;
444 3602 : *(*bufptr)++ = (unsigned char)((us.u >> 8) & 0xff);
445 3602 : *(*bufptr)++ = (unsigned char)(us.u & 0xff);
446 3602 : }
447 :
448 224 : static int get_unsigned_short_le(const unsigned char **data)
449 : {
450 : int result;
451 :
452 224 : result = *(*data) | (*(*data + 1) << 8);
453 224 : *data += 2;
454 224 : return result;
455 : }
456 :
457 192 : static int get_unsigned_short_be(const unsigned char **data)
458 : {
459 : int result;
460 :
461 192 : result = *(*data + 1) | (*(*data) << 8);
462 192 : *data += 2;
463 192 : return result;
464 : }
465 :
466 64 : static void put_unsigned_short_le(unsigned short data, unsigned char **bufptr)
467 : {
468 64 : *(*bufptr)++ = (unsigned char)(data & 0xff);
469 64 : *(*bufptr)++ = (unsigned char)((data >> 8) & 0xff);
470 64 : }
471 :
472 64 : static void put_unsigned_short_be(unsigned short data, unsigned char **bufptr)
473 : {
474 64 : *(*bufptr)++ = (unsigned char)((data >> 8) & 0xff);
475 64 : *(*bufptr)++ = (unsigned char)(data & 0xff);
476 64 : }
477 :
478 22619 : static int get_short(blxcontext_t *ctx, const unsigned char **data)
479 : {
480 :
481 22619 : if (ctx->endian == LITTLEENDIAN)
482 13361 : return get_short_le(data);
483 : else
484 9258 : return get_short_be(data);
485 : }
486 :
487 416 : static int get_unsigned_short(blxcontext_t *ctx, const unsigned char **data)
488 : {
489 :
490 416 : if (ctx->endian == LITTLEENDIAN)
491 224 : return get_unsigned_short_le(data);
492 : else
493 192 : return get_unsigned_short_be(data);
494 : }
495 :
496 7204 : static void put_short(blxcontext_t *ctx, short data, unsigned char **bufptr)
497 : {
498 7204 : if (ctx->endian == LITTLEENDIAN)
499 3602 : put_short_le(data, bufptr);
500 : else
501 3602 : put_short_be(data, bufptr);
502 7204 : }
503 :
504 128 : static void put_unsigned_short(blxcontext_t *ctx, unsigned short data,
505 : unsigned char **bufptr)
506 : {
507 128 : if (ctx->endian == LITTLEENDIAN)
508 64 : put_unsigned_short_le(data, bufptr);
509 : else
510 64 : put_unsigned_short_be(data, bufptr);
511 128 : }
512 :
513 : typedef union
514 : {
515 : int i;
516 : unsigned int u;
517 : } unionint;
518 :
519 39 : static int get_int32(blxcontext_t *ctx, const unsigned char **data)
520 : {
521 : /* We assume two's complement representation for this to work */
522 39 : unionint result = {0};
523 :
524 39 : if (ctx->endian == LITTLEENDIAN)
525 21 : result.u = *(*data) | (*(*data + 1) << 8) | (*(*data + 2) << 16) |
526 21 : ((unsigned)*(*data + 3) << 24);
527 : else
528 18 : result.u = *(*data + 3) | (*(*data + 2) << 8) | (*(*data + 1) << 16) |
529 18 : ((unsigned)*(*data) << 24);
530 39 : *data += 4;
531 39 : return result.i;
532 : }
533 :
534 76 : static void put_int32(blxcontext_t *ctx, int data, unsigned char **bufptr)
535 : {
536 : /* We assume two's complement representation for this to work */
537 76 : unionint ui = {0};
538 76 : ui.i = data;
539 76 : if (ctx->endian == LITTLEENDIAN)
540 : {
541 38 : *(*bufptr)++ = (unsigned char)(ui.u & 0xff);
542 38 : *(*bufptr)++ = (unsigned char)((ui.u >> 8) & 0xff);
543 38 : *(*bufptr)++ = (unsigned char)((ui.u >> 16) & 0xff);
544 38 : *(*bufptr)++ = (unsigned char)((ui.u >> 24) & 0xff);
545 : }
546 : else
547 : {
548 38 : *(*bufptr)++ = (unsigned char)((ui.u >> 24) & 0xff);
549 38 : *(*bufptr)++ = (unsigned char)((ui.u >> 16) & 0xff);
550 38 : *(*bufptr)++ = (unsigned char)((ui.u >> 8) & 0xff);
551 38 : *(*bufptr)++ = (unsigned char)(ui.u & 0xff);
552 : }
553 76 : }
554 :
555 208 : static int get_unsigned32(blxcontext_t *ctx, const unsigned char **data)
556 : {
557 : int result;
558 :
559 208 : if (ctx->endian == LITTLEENDIAN)
560 112 : result = *(*data) | (*(*data + 1) << 8) | (*(*data + 2) << 16) |
561 112 : ((unsigned)*(*data + 3) << 24);
562 : else
563 96 : result = *(*data + 3) | (*(*data + 2) << 8) | (*(*data + 1) << 16) |
564 96 : ((unsigned)*(*data) << 24);
565 208 : *data += 4;
566 208 : return result;
567 : }
568 :
569 : /* Check native endian order */
570 136 : static int is_big_endian(void)
571 : {
572 136 : short int word = 0x0001;
573 136 : char *byte = (char *)&word;
574 136 : return (byte[0] ? 0 : 1);
575 : }
576 :
577 32 : static double doubleSWAP(double df)
578 : {
579 32 : CPL_SWAP64PTR(&df);
580 32 : return df;
581 : }
582 :
583 52 : static double get_double(blxcontext_t *ctx, const unsigned char **data)
584 : {
585 : double result;
586 52 : memcpy(&result, *data, sizeof(double));
587 104 : if ((is_big_endian() && ctx->endian == LITTLEENDIAN) ||
588 104 : (!is_big_endian() && ctx->endian == BIGENDIAN))
589 24 : result = doubleSWAP(result);
590 :
591 52 : *data += sizeof(double);
592 :
593 52 : return result;
594 : }
595 :
596 16 : static void put_double(blxcontext_t *ctx, double data, unsigned char **bufptr)
597 : {
598 32 : if ((is_big_endian() && ctx->endian == LITTLEENDIAN) ||
599 32 : (!is_big_endian() && ctx->endian == BIGENDIAN))
600 8 : data = doubleSWAP(data);
601 16 : memcpy(*bufptr, &data, sizeof(double));
602 16 : *bufptr += sizeof(double);
603 16 : }
604 :
605 64 : static void put_cellindex_entry(blxcontext_t *ctx, struct cellindex_s *ci,
606 : unsigned char **buffer)
607 : {
608 64 : put_int32(ctx, (int)ci->offset, buffer);
609 64 : put_unsigned_short(ctx, (unsigned short)ci->datasize, buffer);
610 64 : put_unsigned_short(ctx, (unsigned short)ci->compdatasize, buffer);
611 64 : }
612 :
613 : /* Transpose matrix in-place */
614 728 : static void transpose(blxdata *data, int rows, int cols)
615 : {
616 : int i, j;
617 : blxdata tmp;
618 :
619 23448 : for (i = 0; i < rows; i++)
620 532320 : for (j = i + 1; j < cols; j++)
621 : {
622 509600 : tmp = data[i * cols + j];
623 509600 : data[i * cols + j] = data[j * cols + i];
624 509600 : data[j * cols + i] = tmp;
625 : }
626 728 : }
627 :
628 : struct lutentry_s
629 : {
630 : blxdata value;
631 : int frequency;
632 : };
633 :
634 66234 : static int lutcmp(const void *aa, const void *bb)
635 : {
636 66234 : const struct lutentry_s *a = aa, *b = bb;
637 :
638 66234 : return b->frequency - a->frequency;
639 : }
640 :
641 32 : int blx_encode_celldata(blxcontext_t *ctx, blxdata *indata, int side,
642 : unsigned char *outbuf, CPL_UNUSED int outbufsize)
643 : {
644 32 : unsigned char *p = outbuf, *tmpdata, *coutstart, *cout = NULL;
645 : int level, cn, coutsize, zeros;
646 32 : blxdata *vdec = NULL, *vdiff = NULL, *c[4] = {NULL}, *tc1, *clut,
647 : *indata_scaled;
648 :
649 : struct lutentry_s lut[256];
650 32 : int lutsize = 0;
651 :
652 : int i, j;
653 :
654 32 : memset(&lut, 0, sizeof(lut));
655 32 : lut[0].value = 0;
656 :
657 32 : *p++ = (unsigned char)(side / 32 - 4); /* Resolution */
658 :
659 : /* Allocated memory for buffers */
660 32 : indata_scaled = BLXmalloc(sizeof(blxdata) * side * side);
661 32 : vdec = BLXmalloc(sizeof(blxdata) * side * side / 2);
662 32 : vdiff = BLXmalloc(sizeof(blxdata) * side * side / 2);
663 160 : for (cn = 0; cn < 4; cn++)
664 128 : c[cn] = BLXmalloc(sizeof(blxdata) * side * side / 4);
665 32 : tc1 = BLXmalloc(sizeof(blxdata) * side * side / 4);
666 32 : tmpdata = BLXmalloc(5 * 4 * side * side / 4);
667 :
668 : /* Scale indata and process undefined values*/
669 524320 : for (i = 0; i < side * side; i++)
670 : {
671 524288 : if ((indata[i] == BLX_UNDEF) && ctx->fillundef)
672 0 : indata[i] = (blxdata)ctx->fillundefval;
673 : /* cppcheck-suppress uninitdata */
674 524288 : indata_scaled[i] = (blxdata)(indata[i] / ctx->zscale);
675 : }
676 :
677 32 : indata = indata_scaled;
678 :
679 32 : cout = tmpdata;
680 :
681 192 : for (level = 0; level < 5; level++)
682 : {
683 160 : if (ctx->debug)
684 : {
685 0 : BLXdebug1("\nlevel=%d\n", level);
686 : }
687 160 : decimate_vert(indata, side, side, vdec, vdiff);
688 160 : decimate_horiz(vdec, side / 2, side, c[0], c[1]);
689 160 : decimate_horiz(vdiff, side / 2, side, c[2], c[3]);
690 :
691 : /* For some reason the matrix is transposed if the lut is used for
692 : * vdec_hdiff */
693 4128 : for (i = 0; i < side / 2; i++)
694 178560 : for (j = 0; j < side / 2; j++)
695 : {
696 174592 : tc1[j * side / 2 + i] = c[1][i * side / 2 + j];
697 174592 : tc1[i * side / 2 + j] = c[1][j * side / 2 + i];
698 : }
699 :
700 640 : for (cn = 1; cn < 4; cn++)
701 : {
702 : /* Use the possibly transposed version of c when building lut */
703 480 : if (cn == 1)
704 160 : clut = tc1;
705 : else
706 320 : clut = c[cn];
707 :
708 480 : lutsize = 0;
709 480 : coutstart = cout;
710 524256 : for (i = 0; i < side * side / 4; i++)
711 : {
712 : /* Find element in lookup table */
713 7448770 : for (j = 0; (j < lutsize) && (lut[j].value != clut[i]); j++)
714 : ;
715 :
716 523776 : if (clut[i] != 0)
717 : {
718 456638 : if (j == lutsize)
719 : {
720 18636 : lut[lutsize].value = clut[i];
721 18636 : lut[lutsize].frequency = 1;
722 18636 : lutsize++;
723 18636 : if (lutsize >= 255)
724 0 : break;
725 : }
726 : else
727 438002 : lut[j].frequency++;
728 : }
729 : }
730 480 : if (lutsize < 255)
731 : {
732 : /* Since the Huffman table is arranged to let smaller number
733 : occupy less bits after the compression, the lookup table is sorted on
734 : frequency */
735 480 : qsort(lut, lutsize, sizeof(struct lutentry_s), lutcmp);
736 :
737 480 : zeros = 0;
738 524256 : for (i = 0; i < side * side / 4; i++)
739 : {
740 523776 : if (clut[i] == 0)
741 67138 : zeros++;
742 523776 : if (((zeros > 0) && (clut[i] != 0)) ||
743 475588 : (zeros >= 0x100 - lutsize))
744 : {
745 48188 : *cout++ = (unsigned char)(0x100 - zeros);
746 48188 : zeros = 0;
747 : }
748 523776 : if (clut[i] != 0)
749 : {
750 3788590 : for (j = 0; (j < lutsize) && (lut[j].value != clut[i]);
751 3331950 : j++)
752 : ;
753 456638 : *cout++ = (unsigned char)j;
754 : }
755 : }
756 480 : if (zeros > 0)
757 24 : *cout++ = (unsigned char)(0x100 - zeros);
758 : }
759 : /* Use the lookuptable only when it pays off to do do.
760 : For some reason there cannot be lookup tables in level 4.
761 : Otherwise Mapsend crashes. */
762 480 : coutsize = (int)(cout - coutstart);
763 480 : if ((lutsize < 255) &&
764 480 : (coutsize + 2 * lutsize + 1 < 2 * side * side / 4) &&
765 : (level < 4))
766 : {
767 304 : *p++ = (unsigned char)(lutsize + 1);
768 14042 : for (j = 0; j < lutsize; j++)
769 13738 : put_short_le(lut[j].value, &p);
770 304 : put_short_le((short)coutsize, &p);
771 :
772 304 : if (ctx->debug)
773 : {
774 0 : BLXdebug2("n=%d dlen=%d\n", lutsize + 1, coutsize);
775 0 : BLXdebug0("lut={");
776 0 : for (i = 0; i < lutsize; i++)
777 0 : BLXdebug1("%d, ", lut[i].value);
778 0 : BLXdebug0("}\n");
779 : }
780 : }
781 : else
782 : {
783 176 : *p++ = 0;
784 176 : cout = coutstart;
785 6832 : for (i = 0; i < side * side / 4; i++)
786 6656 : put_short(ctx, c[cn][i], &cout);
787 : }
788 : }
789 :
790 160 : side >>= 1;
791 160 : indata = c[0];
792 : }
793 :
794 32 : memcpy(p, tmpdata, cout - tmpdata);
795 32 : p += cout - tmpdata;
796 :
797 544 : for (i = 0; i < side * side; i++)
798 512 : put_short(ctx, indata[i], &p);
799 :
800 32 : *p++ = 0;
801 :
802 32 : BLXfree(indata_scaled);
803 32 : BLXfree(vdec);
804 32 : BLXfree(vdiff);
805 160 : for (cn = 0; cn < 4; cn++)
806 128 : BLXfree(c[cn]);
807 32 : BLXfree(tc1);
808 32 : BLXfree(tmpdata);
809 :
810 32 : return (int)(p - outbuf);
811 : }
812 :
813 192 : STATIC blxdata *decode_celldata(blxcontext_t *ctx, const unsigned char *inbuf,
814 : int len, int *side, blxdata *outbuf,
815 : int outbufsize, int overviewlevel)
816 : {
817 192 : const unsigned char *inptr = inbuf;
818 : int resolution, l_div, level, c, n, i, j, dpos, tmp, a, value, l_index,
819 : step, cellsize;
820 192 : int baseside[12] = {0};
821 : blxdata *base, *diff;
822 : struct component_s linfo[MAXLEVELS][MAXCOMPONENTS];
823 :
824 192 : if (len < 1)
825 : {
826 0 : BLXerror0("Cell corrupt");
827 0 : return NULL;
828 : }
829 192 : resolution = *inptr++;
830 192 : len--;
831 :
832 192 : tmp = (resolution + 4) * 32;
833 2304 : for (l_div = 1; l_div < 12; l_div++)
834 2112 : baseside[l_div - 1] = tmp >> l_div;
835 :
836 192 : if (side != NULL)
837 0 : *side = tmp >> overviewlevel;
838 :
839 192 : cellsize = tmp * tmp;
840 192 : if (outbufsize < cellsize * (int)sizeof(blxdata))
841 : {
842 0 : BLXerror0("Cell will not fit in output buffer\n");
843 0 : return NULL;
844 : }
845 :
846 192 : if (outbuf == NULL)
847 : {
848 0 : BLXerror0("outbuf is NULL");
849 0 : return NULL;
850 : }
851 :
852 192 : if (ctx->debug)
853 : {
854 0 : BLXdebug0("==============================\n");
855 : }
856 :
857 : /* Clear level info structure */
858 192 : memset(linfo, 0, sizeof(linfo));
859 :
860 192 : base = BLXmalloc(sizeof(blxdata) * 2 * baseside[0] * baseside[0]);
861 192 : diff = BLXmalloc(sizeof(blxdata) * 2 * baseside[0] * baseside[0]);
862 192 : if (base == NULL || diff == NULL)
863 : {
864 0 : BLXerror0("Not enough memory\n");
865 0 : outbuf = NULL;
866 0 : goto error;
867 : }
868 :
869 1152 : for (level = 0; level < 5; level++)
870 : {
871 3840 : for (c = 1; c < 4; c++)
872 : {
873 2880 : if (len < 1)
874 : {
875 0 : BLXerror0("Cell corrupt");
876 0 : outbuf = NULL;
877 0 : goto error;
878 : }
879 2880 : n = *inptr++;
880 2880 : len--;
881 2880 : linfo[level][c].n = n;
882 2880 : if (n > 0)
883 : {
884 2144 : linfo[level][c].lut = BLXmalloc(sizeof(blxdata) * (n - 1));
885 2144 : if (len < (int)sizeof(short) * n)
886 : {
887 0 : BLXerror0("Cell corrupt");
888 0 : outbuf = NULL;
889 0 : goto error;
890 : }
891 44228 : for (i = 0; i < n - 1; i++)
892 42084 : linfo[level][c].lut[i] = (blxdata)get_short_le(&inptr);
893 2144 : linfo[level][c].dlen = get_short_le(&inptr);
894 2144 : if (linfo[level][c].dlen < 0)
895 : {
896 0 : BLXerror0("Cell corrupt");
897 0 : outbuf = NULL;
898 0 : goto error;
899 : }
900 2144 : len -= sizeof(short) * n;
901 : }
902 : else
903 : {
904 736 : linfo[level][c].dlen = 0;
905 : }
906 : }
907 : }
908 :
909 1152 : for (level = 0; level < 5; level++)
910 : {
911 960 : if (ctx->debug)
912 : {
913 0 : BLXdebug1("\nlevel=%d\n", level);
914 : }
915 :
916 960 : linfo[level][0].data =
917 960 : BLXmalloc(sizeof(blxdata) * baseside[level] * baseside[level]);
918 960 : if (linfo[level][0].data == NULL)
919 : {
920 0 : BLXerror0("Not enough memory\n");
921 0 : outbuf = NULL;
922 0 : goto error;
923 : }
924 :
925 3840 : for (c = 1; c < 4; c++)
926 : {
927 2880 : if (ctx->debug)
928 : {
929 0 : BLXdebug2("n=%d dlen=%d\n", linfo[level][c].n,
930 : linfo[level][c].dlen);
931 0 : BLXdebug0("lut={");
932 0 : for (i = 0; i < linfo[level][c].n - 1; i++)
933 0 : BLXdebug1("%d, ", linfo[level][c].lut[i]);
934 0 : BLXdebug0("}\n");
935 : }
936 :
937 2880 : linfo[level][c].data =
938 2880 : BLXmalloc(sizeof(blxdata) * baseside[level] * baseside[level]);
939 2880 : if (linfo[level][c].data == NULL)
940 : {
941 0 : BLXerror0("Not enough memory\n");
942 0 : outbuf = NULL;
943 0 : goto error;
944 : }
945 :
946 2880 : if (linfo[level][c].n == 0)
947 : {
948 736 : if (len <
949 736 : (int)sizeof(short) * baseside[level] * baseside[level])
950 : {
951 0 : BLXerror0("Cell corrupt");
952 0 : outbuf = NULL;
953 0 : goto error;
954 : }
955 20192 : for (i = 0; i < baseside[level] * baseside[level]; i++)
956 19456 : linfo[level][c].data[i] = (blxdata)get_short(ctx, &inptr);
957 736 : len -= sizeof(short) * baseside[level] * baseside[level];
958 : }
959 : else
960 : {
961 2144 : dpos = 0;
962 2144 : if (len < linfo[level][c].dlen)
963 : {
964 0 : BLXerror0("Cell corrupt");
965 0 : outbuf = NULL;
966 0 : goto error;
967 : }
968 2445070 : for (i = 0; i < linfo[level][c].dlen; i++)
969 : {
970 2442920 : unsigned char v = *inptr++;
971 2442920 : if (v >= linfo[level][c].n - 1)
972 : {
973 549452 : if (dpos + 256 - v > baseside[level] * baseside[level])
974 : {
975 0 : BLXerror0("Cell corrupt\n");
976 0 : outbuf = NULL;
977 0 : goto error;
978 : }
979 : /* coverity[tainted_data] */
980 1779180 : for (j = 0; j < 256 - v; j++)
981 1229730 : linfo[level][c].data[dpos++] = 0;
982 : }
983 : else
984 : {
985 1893470 : if (dpos + 1 > baseside[level] * baseside[level])
986 : {
987 0 : BLXerror0("Cell corrupt\n");
988 0 : outbuf = NULL;
989 0 : goto error;
990 : }
991 1893470 : linfo[level][c].data[dpos++] = linfo[level][c].lut[v];
992 : }
993 : }
994 2144 : len -= linfo[level][c].dlen;
995 2144 : if (c == 1)
996 728 : transpose(linfo[level][c].data, baseside[level],
997 : baseside[level]);
998 : }
999 : if (0 && ctx->debug)
1000 : {
1001 : BLXdebug1("baseside:%d\n", baseside[level]);
1002 : BLXdebug0("data={");
1003 : for (i = 0; i < baseside[level] * baseside[level]; i++)
1004 : BLXdebug1("%d, ", linfo[level][c].data[i]);
1005 : BLXdebug0("}\n");
1006 : }
1007 : }
1008 : }
1009 :
1010 192 : if (len < (int)sizeof(short) * baseside[4] * baseside[4])
1011 : {
1012 0 : BLXerror0("Cell corrupt");
1013 0 : outbuf = NULL;
1014 0 : goto error;
1015 : }
1016 3264 : for (i = 0; i < baseside[4] * baseside[4]; i++)
1017 3072 : linfo[4][0].data[i] = (blxdata)get_short(ctx, &inptr);
1018 192 : len -= sizeof(short) * baseside[4] * baseside[4];
1019 :
1020 992 : for (level = 4; level >= overviewlevel; level--)
1021 : {
1022 800 : if (ctx->debug)
1023 : {
1024 0 : BLXdebug1("baseside:%d\n", baseside[level]);
1025 0 : BLXdebug0("inbase={");
1026 0 : for (i = 0; i < baseside[level] * baseside[level]; i++)
1027 0 : BLXdebug1("%d, ", linfo[level][0].data[i]);
1028 0 : BLXdebug0("}\n");
1029 0 : BLXdebug0("indiff={");
1030 0 : for (i = 0; i < baseside[level] * baseside[level]; i++)
1031 0 : BLXdebug1("%d, ", linfo[level][1].data[i]);
1032 0 : BLXdebug0("}\n");
1033 : }
1034 :
1035 800 : reconstruct_horiz(linfo[level][0].data, linfo[level][1].data,
1036 800 : baseside[level], baseside[level], base);
1037 800 : if (ctx->debug)
1038 : {
1039 0 : BLXdebug0("base={");
1040 0 : for (i = 0; i < baseside[level] * baseside[level]; i++)
1041 0 : BLXdebug1("%d, ", base[i]);
1042 0 : BLXdebug0("}\n");
1043 : }
1044 :
1045 800 : reconstruct_horiz(linfo[level][2].data, linfo[level][3].data,
1046 800 : baseside[level], baseside[level], diff);
1047 800 : if (ctx->debug)
1048 : {
1049 0 : BLXdebug0("diff={");
1050 0 : for (i = 0; i < baseside[level] * baseside[level]; i++)
1051 0 : BLXdebug1("%d, ", diff[i]);
1052 0 : BLXdebug0("}\n");
1053 : }
1054 800 : if (level > overviewlevel)
1055 608 : reconstruct_vert(base, diff, baseside[level], 2 * baseside[level],
1056 608 : linfo[level - 1][0].data);
1057 : else
1058 192 : reconstruct_vert(base, diff, baseside[level], 2 * baseside[level],
1059 : outbuf);
1060 : }
1061 :
1062 192 : if (overviewlevel == 0)
1063 : {
1064 128 : if (len < 1)
1065 : {
1066 0 : BLXerror0("Cell corrupt");
1067 0 : outbuf = NULL;
1068 0 : goto error;
1069 : }
1070 128 : a = *((char *)inptr++);
1071 128 : len--;
1072 128 : l_index = 0;
1073 128 : while (len >= 3)
1074 : {
1075 0 : step = inptr[0] | (inptr[1] << 8);
1076 0 : inptr += 2;
1077 0 : value = *((char *)inptr++);
1078 0 : len -= 3;
1079 :
1080 0 : l_index += step;
1081 :
1082 0 : if (value & 1)
1083 0 : value = (value - 1) / 2 - a;
1084 : else
1085 0 : value = value / 2 + a;
1086 :
1087 0 : if (l_index >= cellsize)
1088 : {
1089 0 : BLXerror0("Cell data corrupt\n");
1090 0 : outbuf = NULL;
1091 0 : goto error;
1092 : }
1093 :
1094 0 : outbuf[l_index] = outbuf[l_index] + (blxdata)value;
1095 : }
1096 :
1097 128 : if (len != 0)
1098 0 : BLXdebug1("remaining len=%d", len);
1099 : }
1100 : else
1101 : {
1102 64 : if (len != 1)
1103 0 : BLXdebug1("remaining len=%d", len);
1104 : }
1105 :
1106 : /* Scale data */
1107 3145920 : for (i = 0; i < cellsize; i++)
1108 : {
1109 3145730 : int val = outbuf[i] * (blxdata)ctx->zscale;
1110 3145730 : if (val < SHRT_MIN)
1111 66429 : outbuf[i] = SHRT_MIN;
1112 3079300 : else if (val > SHRT_MAX)
1113 673282 : outbuf[i] = SHRT_MAX;
1114 : else
1115 2406020 : outbuf[i] = (blxdata)val;
1116 : }
1117 :
1118 192 : error:
1119 192 : if (base != NULL)
1120 192 : BLXfree(base);
1121 192 : if (diff != NULL)
1122 192 : BLXfree(diff);
1123 :
1124 : /* Free allocated memory */
1125 1152 : for (level = 4; level >= 0; level--)
1126 4800 : for (c = 0; c < 4; c++)
1127 : {
1128 3840 : if (linfo[level][c].lut)
1129 2144 : BLXfree(linfo[level][c].lut);
1130 3840 : if (linfo[level][c].data)
1131 3840 : BLXfree(linfo[level][c].data);
1132 : }
1133 :
1134 192 : return outbuf;
1135 : }
1136 :
1137 17 : blxcontext_t *blx_create_context()
1138 : {
1139 : blxcontext_t *c;
1140 :
1141 17 : c = BLXmalloc(sizeof(blxcontext_t));
1142 :
1143 17 : memset(c, 0, sizeof(blxcontext_t));
1144 :
1145 17 : c->cell_ysize = c->cell_xsize = 128;
1146 :
1147 17 : c->minval = 32767;
1148 17 : c->maxval = -32768;
1149 :
1150 17 : c->debug = 0;
1151 :
1152 17 : c->zscale = 1;
1153 :
1154 17 : c->fillundef = 1;
1155 17 : c->fillundefval = 0;
1156 :
1157 17 : return c;
1158 : }
1159 :
1160 17 : void blx_free_context(blxcontext_t *ctx)
1161 : {
1162 17 : if (ctx->cellindex)
1163 15 : BLXfree(ctx->cellindex);
1164 :
1165 17 : BLXfree(ctx);
1166 17 : }
1167 :
1168 0 : void blxprintinfo(blxcontext_t *ctx)
1169 : {
1170 0 : BLXnotice2("Lat: %f Lon: %f\n", ctx->lat, ctx->lon);
1171 0 : BLXnotice2("Pixelsize: Lat: %f Lon: %f\n", 3600 * ctx->pixelsize_lat,
1172 : 3600 * ctx->pixelsize_lon);
1173 0 : BLXnotice2("Size %dx%d\n", ctx->xsize, ctx->ysize);
1174 0 : BLXnotice2("Cell size %dx%d\n", ctx->cell_xsize, ctx->cell_ysize);
1175 0 : BLXnotice2("Cell grid %dx%d\n", ctx->cell_cols, ctx->cell_rows);
1176 0 : BLXnotice1("Ysize scale factor: %d\n", ctx->zscale);
1177 0 : BLXnotice1("Max Ysize: %d\n", ctx->zscale * ctx->maxval);
1178 0 : BLXnotice1("Min Ysize: %d\n", ctx->zscale * ctx->minval);
1179 0 : BLXnotice1("Max chunksize: %d\n", ctx->maxchunksize);
1180 0 : }
1181 :
1182 2410 : int blx_checkheader(const char *header)
1183 : {
1184 2410 : const short *signature = (const short *)header;
1185 :
1186 4813 : return ((signature[0] == 0x4) && (signature[1] == 0x66)) ||
1187 2403 : ((signature[0] == 0x400) && (signature[1] == 0x6600));
1188 : }
1189 :
1190 4 : static void blx_generate_header(blxcontext_t *ctx, unsigned char *header)
1191 : {
1192 4 : unsigned char *hptr = header;
1193 :
1194 4 : memset(header, 0, 102);
1195 :
1196 : /* Write signature */
1197 4 : put_short(ctx, 0x4, &hptr); // 0
1198 4 : put_short(ctx, 0x66, &hptr); // 2
1199 :
1200 4 : put_int32(ctx, ctx->cell_xsize * ctx->cell_cols, &hptr); // 4
1201 4 : put_int32(ctx, ctx->cell_ysize * ctx->cell_rows, &hptr); // 8
1202 :
1203 4 : put_short(ctx, (short)ctx->cell_xsize, &hptr); // 12
1204 4 : put_short(ctx, (short)ctx->cell_ysize, &hptr); // 14
1205 :
1206 4 : put_short(ctx, (short)ctx->cell_cols, &hptr); // 16
1207 4 : put_short(ctx, (short)ctx->cell_rows, &hptr); // 18
1208 :
1209 4 : put_double(ctx, ctx->lon, &hptr); // 20
1210 4 : put_double(ctx, -ctx->lat, &hptr); // 28
1211 :
1212 4 : put_double(ctx, ctx->pixelsize_lon, &hptr); // 36
1213 4 : put_double(ctx, -ctx->pixelsize_lat, &hptr); // 44
1214 :
1215 4 : put_short(ctx, (short)ctx->minval, &hptr); // 52
1216 4 : put_short(ctx, (short)ctx->maxval, &hptr); // 54
1217 4 : put_short(ctx, (short)ctx->zscale, &hptr); // 56
1218 4 : put_int32(ctx, ctx->maxchunksize, &hptr); // 58
1219 : // 62
1220 4 : }
1221 :
1222 32 : int blx_writecell(blxcontext_t *ctx, blxdata *cell, int cellrow, int cellcol)
1223 : {
1224 32 : unsigned char *uncompbuf = NULL, *outbuf = NULL;
1225 : int bufsize, uncompsize, compsize;
1226 32 : int status = 0;
1227 : int i, allundef;
1228 :
1229 : /* Calculate statistics and find out if all elements have undefined values
1230 : */
1231 32 : allundef = 1;
1232 524320 : for (i = 0; i < ctx->cell_xsize * ctx->cell_ysize; i++)
1233 : {
1234 524288 : if (cell[i] > ctx->maxval)
1235 224 : ctx->maxval = cell[i];
1236 524288 : if (cell[i] < ctx->minval)
1237 14 : ctx->minval = cell[i];
1238 524288 : if (cell[i] != BLX_UNDEF)
1239 524288 : allundef = 0;
1240 : }
1241 :
1242 32 : if (allundef)
1243 0 : return status;
1244 :
1245 32 : if (ctx->debug)
1246 0 : BLXdebug2("Writing cell (%d,%d)\n", cellrow, cellcol);
1247 :
1248 32 : if (!ctx->open)
1249 : {
1250 0 : status = -3;
1251 0 : goto error;
1252 : }
1253 :
1254 32 : if ((cellrow >= ctx->cell_rows) || (cellcol >= ctx->cell_cols))
1255 : {
1256 0 : status = -2;
1257 0 : goto error;
1258 : }
1259 :
1260 32 : bufsize = sizeof(blxdata) * ctx->cell_xsize * ctx->cell_ysize + 1024;
1261 32 : uncompbuf = BLXmalloc(bufsize);
1262 32 : outbuf = BLXmalloc(bufsize);
1263 :
1264 : uncompsize =
1265 32 : blx_encode_celldata(ctx, cell, ctx->cell_xsize, uncompbuf, bufsize);
1266 32 : compsize = compress_chunk(uncompbuf, uncompsize, outbuf, bufsize);
1267 32 : if (compsize < 0)
1268 : {
1269 0 : BLXerror0("Couldn't compress chunk");
1270 0 : status = -1;
1271 0 : goto error;
1272 : }
1273 :
1274 32 : if (uncompsize > ctx->maxchunksize)
1275 12 : ctx->maxchunksize = uncompsize;
1276 :
1277 32 : ctx->cellindex[cellrow * ctx->cell_cols + cellcol].offset =
1278 32 : (int)BLXftell(ctx->fh);
1279 32 : ctx->cellindex[cellrow * ctx->cell_cols + cellcol].datasize = uncompsize;
1280 32 : ctx->cellindex[cellrow * ctx->cell_cols + cellcol].compdatasize = compsize;
1281 :
1282 32 : if ((int)BLXfwrite(outbuf, 1, compsize, ctx->fh) != compsize)
1283 : {
1284 0 : status = -1;
1285 0 : goto error;
1286 : }
1287 :
1288 32 : error:
1289 32 : if (uncompbuf)
1290 32 : BLXfree(uncompbuf);
1291 32 : if (outbuf)
1292 32 : BLXfree(outbuf);
1293 32 : return status;
1294 : }
1295 :
1296 17 : int blxopen(blxcontext_t *ctx, const char *filename, const char *rw)
1297 : {
1298 : unsigned char header[102];
1299 17 : int signature[2] = {0};
1300 : int i, j;
1301 : struct cellindex_s *ci;
1302 :
1303 17 : if (!strcmp(rw, "r") || !strcmp(rw, "rb"))
1304 13 : ctx->write = 0;
1305 4 : else if (!strcmp(rw, "w") || !strcmp(rw, "wb"))
1306 4 : ctx->write = 1;
1307 : else
1308 0 : goto error;
1309 :
1310 17 : ctx->fh = BLXfopen(filename, rw);
1311 :
1312 17 : if (ctx->fh == NULL)
1313 2 : goto error;
1314 :
1315 15 : if (ctx->write)
1316 : {
1317 2 : unsigned char *hptr = header;
1318 2 : blx_generate_header(ctx, header);
1319 :
1320 2 : if (BLXfwrite(header, 1, 102, ctx->fh) != 102)
1321 0 : goto error;
1322 :
1323 4 : ctx->cellindex = BLXmalloc(sizeof(struct cellindex_s) * ctx->cell_rows *
1324 2 : ctx->cell_cols);
1325 2 : if (ctx->cellindex == NULL)
1326 : {
1327 0 : goto error;
1328 : }
1329 2 : memset(ctx->cellindex, 0,
1330 2 : sizeof(struct cellindex_s) * ctx->cell_rows * ctx->cell_cols);
1331 :
1332 : /* Write the empty cell index (this will be overwritten when we have
1333 : * actual data)*/
1334 10 : for (i = 0; i < ctx->cell_rows; i++)
1335 40 : for (j = 0; j < ctx->cell_cols; j++)
1336 : {
1337 32 : hptr = header;
1338 32 : put_cellindex_entry(
1339 32 : ctx, ctx->cellindex + i * ctx->cell_cols + j, &hptr);
1340 32 : if ((int)BLXfwrite(header, 1, hptr - header, ctx->fh) !=
1341 32 : (int)(hptr - header))
1342 0 : goto error;
1343 : }
1344 : }
1345 : else
1346 : {
1347 13 : const unsigned char *hptr = header;
1348 :
1349 : /* Read header */
1350 13 : if (BLXfread(header, 1, 102, ctx->fh) != 102)
1351 0 : goto error;
1352 :
1353 13 : signature[0] = get_short_le(&hptr);
1354 13 : signature[1] = get_short_le(&hptr);
1355 :
1356 : /* Determine if the endianness of the BLX file */
1357 13 : if ((signature[0] == 0x4) && (signature[1] == 0x66))
1358 7 : ctx->endian = LITTLEENDIAN;
1359 : else
1360 : {
1361 6 : hptr = header;
1362 6 : signature[0] = get_short_be(&hptr);
1363 6 : signature[1] = get_short_be(&hptr);
1364 6 : if ((signature[0] == 0x4) && (signature[1] == 0x66))
1365 6 : ctx->endian = BIGENDIAN;
1366 : else
1367 0 : goto error;
1368 : }
1369 :
1370 13 : ctx->xsize = get_int32(ctx, &hptr);
1371 13 : ctx->ysize = get_int32(ctx, &hptr);
1372 13 : if (ctx->xsize <= 0 || ctx->ysize <= 0)
1373 : {
1374 0 : BLXerror0("Invalid raster size");
1375 0 : goto error;
1376 : }
1377 :
1378 13 : ctx->cell_xsize = get_short(ctx, &hptr);
1379 13 : ctx->cell_ysize = get_short(ctx, &hptr);
1380 13 : if (ctx->cell_xsize <= 0 || ctx->cell_ysize <= 0)
1381 : {
1382 0 : BLXerror0("Invalid cell size");
1383 0 : goto error;
1384 : }
1385 :
1386 13 : ctx->cell_cols = get_short(ctx, &hptr);
1387 13 : ctx->cell_rows = get_short(ctx, &hptr);
1388 13 : if (ctx->cell_cols <= 0 || ctx->cell_cols > 10000 ||
1389 13 : ctx->cell_rows <= 0 || ctx->cell_rows > 10000)
1390 : {
1391 0 : BLXerror0("Invalid cell number");
1392 0 : goto error;
1393 : }
1394 :
1395 13 : ctx->lon = get_double(ctx, &hptr);
1396 13 : ctx->lat = -get_double(ctx, &hptr);
1397 :
1398 13 : ctx->pixelsize_lon = get_double(ctx, &hptr);
1399 13 : ctx->pixelsize_lat = -get_double(ctx, &hptr);
1400 :
1401 13 : ctx->minval = get_short(ctx, &hptr);
1402 13 : ctx->maxval = get_short(ctx, &hptr);
1403 13 : ctx->zscale = get_short(ctx, &hptr);
1404 13 : ctx->maxchunksize = get_int32(ctx, &hptr);
1405 :
1406 26 : ctx->cellindex = BLXmalloc(sizeof(struct cellindex_s) * ctx->cell_rows *
1407 13 : ctx->cell_cols);
1408 13 : if (ctx->cellindex == NULL)
1409 : {
1410 0 : goto error;
1411 : }
1412 :
1413 65 : for (i = 0; i < ctx->cell_rows; i++)
1414 260 : for (j = 0; j < ctx->cell_cols; j++)
1415 : {
1416 : /* Read cellindex entry */
1417 208 : if (BLXfread(header, 1, 8, ctx->fh) != 8)
1418 0 : goto error;
1419 208 : hptr = header;
1420 :
1421 208 : ci = &ctx->cellindex[i * ctx->cell_cols + j];
1422 208 : ci->offset = get_unsigned32(ctx, &hptr);
1423 208 : ci->datasize = get_unsigned_short(ctx, &hptr);
1424 208 : ci->compdatasize = get_unsigned_short(ctx, &hptr);
1425 : }
1426 : }
1427 15 : ctx->open = 1;
1428 :
1429 15 : return 0;
1430 :
1431 2 : error:
1432 2 : return -1;
1433 : }
1434 :
1435 15 : int blxclose(blxcontext_t *ctx)
1436 : {
1437 : unsigned char header[102], *hptr;
1438 15 : int i, j, status = 0;
1439 :
1440 15 : if (ctx->write)
1441 : {
1442 : /* Write updated header and cellindex */
1443 2 : if (BLXfseek(ctx->fh, 0, SEEK_SET) != 0)
1444 : {
1445 0 : status = -1;
1446 0 : goto error;
1447 : }
1448 :
1449 2 : blx_generate_header(ctx, header);
1450 :
1451 2 : if (BLXfwrite(header, 1, 102, ctx->fh) != 102)
1452 : {
1453 0 : status = -1;
1454 0 : goto error;
1455 : }
1456 10 : for (i = 0; i < ctx->cell_rows; i++)
1457 40 : for (j = 0; j < ctx->cell_cols; j++)
1458 : {
1459 32 : hptr = header;
1460 32 : put_cellindex_entry(
1461 32 : ctx, ctx->cellindex + i * ctx->cell_cols + j, &hptr);
1462 32 : if ((int)BLXfwrite(header, 1, hptr - header, ctx->fh) !=
1463 32 : (int)(hptr - header))
1464 : {
1465 0 : status = -1;
1466 0 : break;
1467 : }
1468 : }
1469 : }
1470 15 : ctx->open = 1;
1471 :
1472 15 : error:
1473 15 : if (ctx->fh)
1474 15 : BLXfclose(ctx->fh);
1475 :
1476 15 : return status;
1477 : }
1478 :
1479 192 : short *blx_readcell(blxcontext_t *ctx, int row, int col, short *buffer,
1480 : int bufsize, int overviewlevel)
1481 : {
1482 : struct cellindex_s *ci;
1483 192 : unsigned char *chunk = NULL, *cchunk = NULL;
1484 192 : blxdata *tmpbuf = NULL;
1485 : int tmpbufsize, i;
1486 192 : short *result = NULL;
1487 : int npoints;
1488 :
1489 192 : if ((ctx == NULL) || (row >= ctx->cell_rows) || (col >= ctx->cell_cols))
1490 0 : return NULL;
1491 :
1492 192 : ci = &ctx->cellindex[row * ctx->cell_cols + col];
1493 :
1494 192 : npoints = (ctx->cell_xsize * ctx->cell_ysize) >> (2 * overviewlevel);
1495 192 : if (bufsize < npoints * (int)sizeof(short))
1496 0 : return NULL;
1497 :
1498 192 : if (ci->datasize == 0)
1499 : {
1500 0 : for (i = 0; i < npoints; i++)
1501 0 : buffer[i] = BLX_UNDEF;
1502 : }
1503 : else
1504 : {
1505 192 : if (BLXfseek(ctx->fh, ci->offset, SEEK_SET) != 0)
1506 0 : goto error;
1507 :
1508 192 : chunk = BLXmalloc(ci->datasize);
1509 192 : cchunk = BLXmalloc(ci->compdatasize);
1510 :
1511 192 : if ((chunk == NULL) || (cchunk == NULL))
1512 0 : goto error;
1513 :
1514 192 : if (BLXfread(cchunk, 1, ci->compdatasize, ctx->fh) != ci->compdatasize)
1515 0 : goto error;
1516 :
1517 384 : if ((unsigned int)uncompress_chunk(cchunk, ci->compdatasize, chunk,
1518 192 : ci->datasize) != ci->datasize)
1519 0 : goto error;
1520 :
1521 192 : tmpbufsize = sizeof(blxdata) * ctx->cell_xsize * ctx->cell_ysize;
1522 192 : tmpbuf = BLXmalloc(tmpbufsize);
1523 192 : if (tmpbuf == NULL)
1524 0 : goto error;
1525 :
1526 192 : if (decode_celldata(ctx, chunk, ci->datasize, NULL, tmpbuf, tmpbufsize,
1527 : overviewlevel) == NULL)
1528 0 : goto error;
1529 :
1530 2184380 : for (i = 0; i < npoints; i++)
1531 2184190 : buffer[i] = tmpbuf[i];
1532 : }
1533 :
1534 192 : result = buffer;
1535 :
1536 192 : error:
1537 192 : if (chunk)
1538 192 : BLXfree(chunk);
1539 192 : if (cchunk)
1540 192 : BLXfree(cchunk);
1541 192 : if (tmpbuf)
1542 192 : BLXfree(tmpbuf);
1543 :
1544 192 : return result;
1545 : }
|
