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