Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: NITF Read/Write Library
4 : * Purpose: ARIDPCM reading code.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : **********************************************************************
8 : * Copyright (c) 2007, Frank Warmerdam
9 : * Copyright (c) 2009, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "nitflib.h"
16 :
17 : #include <cstring>
18 :
19 : #include "gdal.h"
20 : #include "cpl_conv.h"
21 : #include "cpl_error.h"
22 :
23 : constexpr int neighbourhood_size_75[4] = {23, 47, 74, 173};
24 : constexpr int bits_per_level_by_busycode_75[4 /*busy code*/][4 /*level*/] = {
25 : {8, 5, 0, 0}, // BC = 00
26 : {8, 5, 2, 0}, // BC = 01
27 : {8, 6, 4, 0}, // BC = 10
28 : {8, 7, 4, 2}}; // BC = 11
29 :
30 : constexpr int CR075 = 1;
31 :
32 : // Level for each index value.
33 : constexpr int level_index_table[64] = {
34 : 0, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
35 : 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
36 : 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3};
37 :
38 : // List of i,j to linear index macros mappings.
39 : // Note that i is vertical and j is horizontal and progression is
40 : // right to left, bottom to top.
41 :
42 : #define IND(i, j) (ij_index[i + j * 8] - 1)
43 :
44 : constexpr int ij_index[64] = {
45 :
46 : 1, // 0, 0
47 : 18, // 1, 0
48 : 6, // 2, 0
49 : 30, // 3, 0
50 : 3, // 4, 0
51 : 42, // 5, 0
52 : 12, // 6, 0
53 : 54, // 7, 0
54 :
55 : 17, // 0, 1
56 : 19, // 1, 1
57 : 29, // 2, 1
58 : 31, // 3, 1
59 : 41, // 4, 1
60 : 43, // 5, 1
61 : 53, // 6, 1
62 : 55, // 7, 1
63 :
64 : 5, // 0, 2
65 : 21, // 1, 2
66 : 7, // 2, 2
67 : 33, // 3, 2
68 : 11, // 4, 2
69 : 45, // 5, 2
70 : 13, // 6, 2
71 : 57, // 7, 2
72 :
73 : 20, // 0, 3
74 : 22, // 1, 3
75 : 32, // 2, 3
76 : 34, // 3, 3
77 : 44, // 4, 3
78 : 46, // 5, 3
79 : 56, // 6, 3
80 : 58, // 7, 3
81 :
82 : 2, // 0, 4
83 : 24, // 1, 4
84 : 9, // 2, 4
85 : 36, // 3, 4
86 : 4, // 4, 4
87 : 48, // 5, 4
88 : 15, // 6, 4
89 : 60, // 7, 4
90 :
91 : 23, // 0, 5
92 : 25, // 1, 5
93 : 35, // 2, 5
94 : 37, // 3, 5
95 : 47, // 4, 5
96 : 49, // 5, 5
97 : 59, // 6, 5
98 : 61, // 7, 5
99 :
100 : 8, // 0, 6
101 : 27, // 1, 6
102 : 10, // 2, 6
103 : 39, // 3, 6
104 : 14, // 4, 6
105 : 51, // 5, 6
106 : 16, // 6, 6
107 : 63, // 7, 6
108 :
109 : 26, // 0, 7
110 : 28, // 1, 7
111 : 38, // 2, 7
112 : 40, // 3, 7
113 : 50, // 4, 7
114 : 52, // 5, 7
115 : 62, // 6, 7
116 : 64}; // 7, 7
117 :
118 : constexpr int delta_075_level_2_bc_0[32] = {
119 : -71, -49, -38, -32, -27, -23, -20, -17, -14, -12, -10, -8, -6, -4, -3, -1,
120 : 1, 2, 4, 6, 8, 12, 14, 16, 19, 22, 26, 31, 37, 46, 72};
121 : constexpr int delta_075_level_2_bc_1[32] = {
122 : -71, -49, -38, -32, -27, -23, -20, -17, -14, -12, -10, -8, -6, -4, -3, -1,
123 : 1, 2, 4, 6, 8, 12, 14, 16, 19, 22, 26, 31, 37, 46, 72};
124 : constexpr int delta_075_level_2_bc_2[64] = {
125 : -109, -82, -68, -59, -52, -46, -41, -37, -33, -30, -27, -25, -22,
126 : -20, -18, -16, -15, -13, -11, -10, -9, -8, -7, -6, -5, -4,
127 : -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
128 : 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 24,
129 : 26, 28, 31, 35, 38, 42, 47, 52, 60, 69, 85, 118};
130 : constexpr int delta_075_level_2_bc_3[128] = {
131 : -159, -134, -122, -113, -106, -100, -94, -88, -83, -79, -76, -72, -69,
132 : -66, -63, -61, -58, -56, -54, -52, -50, -48, -47, -45, -43, -42,
133 : -40, -39, -37, -36, -35, -33, -32, -31, -30, -29, -28, -27, -25,
134 : -24, -23, -22, -21, -20, -19, -18, -17, -16, -15, -14, -13, -12,
135 : -11, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1,
136 : 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
137 : 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
138 : 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
139 : 41, 42, 43, 45, 48, 52, 56, 60, 64, 68, 73, 79, 85,
140 : 92, 100, 109, 118, 130, 144, 159, 177, 196, 217, 236};
141 : static const int *const delta_075_level_2[4] = {
142 : delta_075_level_2_bc_0, delta_075_level_2_bc_1, delta_075_level_2_bc_2,
143 : delta_075_level_2_bc_3};
144 :
145 : constexpr int delta_075_level_3_bc_1[4] = {-24, -6, 6, 24};
146 : constexpr int delta_075_level_3_bc_2[16] = {-68, -37, -23, -15, -9, -6, -3, -1,
147 : 1, 4, 7, 10, 16, 24, 37, 70};
148 : constexpr int delta_075_level_3_bc_3[16] = {
149 : -117, -72, -50, -36, -25, -17, -10, -5, -1, 3, 7, 14, 25, 45, 82, 166};
150 : static const int *const delta_075_level_3[4] = {nullptr, delta_075_level_3_bc_1,
151 : delta_075_level_3_bc_2,
152 : delta_075_level_3_bc_3};
153 :
154 : constexpr int delta_075_level_4_bc_3[4] = {-47, -8, 4, 43};
155 : static const int *const delta_075_level_4[4] = {nullptr, nullptr, nullptr,
156 : delta_075_level_4_bc_3};
157 :
158 : static const int *const *const delta_075_by_level_by_bc[4] = {
159 : nullptr, delta_075_level_2, delta_075_level_3, delta_075_level_4};
160 :
161 : /************************************************************************/
162 : /* get_bits() */
163 : /************************************************************************/
164 :
165 0 : static int get_bits(unsigned char *buffer, int first_bit, int num_bits)
166 :
167 : {
168 0 : int total = 0;
169 :
170 0 : for (int i = first_bit; i < first_bit + num_bits; i++)
171 : {
172 0 : total *= 2;
173 0 : if (buffer[i >> 3] & (0x80 >> (i & 7)))
174 0 : total++;
175 : }
176 :
177 0 : return total;
178 : }
179 :
180 : /************************************************************************/
181 : /* get_delta() */
182 : /* */
183 : /* Compute the delta value for a particular (i,j) location. */
184 : /************************************************************************/
185 0 : static int get_delta(unsigned char *srcdata, int nInputBytes, int busy_code,
186 : CPL_UNUSED int comrat, int block_offset,
187 : CPL_UNUSED int block_size, int i, int j, int *pbError)
188 :
189 : {
190 0 : CPLAssert(comrat == CR075);
191 0 : const int pixel_index = IND(i, j);
192 0 : const int level_index = level_index_table[pixel_index];
193 0 : const int *bits_per_level = bits_per_level_by_busycode_75[busy_code];
194 0 : const int delta_bits = bits_per_level[level_index];
195 0 : int delta_offset = 0;
196 :
197 0 : *pbError = FALSE;
198 :
199 0 : if (delta_bits == 0)
200 0 : return 0;
201 :
202 0 : if (level_index == 3)
203 0 : delta_offset = bits_per_level[0] + bits_per_level[1] * 3 +
204 0 : bits_per_level[2] * 12 +
205 0 : (pixel_index - 16) * bits_per_level[3];
206 0 : else if (level_index == 2)
207 0 : delta_offset = bits_per_level[0] + bits_per_level[1] * 3 +
208 0 : (pixel_index - 4) * bits_per_level[2];
209 0 : else if (level_index == 1)
210 0 : delta_offset =
211 0 : bits_per_level[0] + (pixel_index - 1) * bits_per_level[1];
212 :
213 0 : if (nInputBytes * 8 < block_offset + delta_offset + delta_bits)
214 : {
215 0 : CPLError(CE_Failure, CPLE_AppDefined, "Input buffer too small");
216 0 : *pbError = TRUE;
217 0 : return 0;
218 : }
219 :
220 : const int delta_raw =
221 0 : get_bits(srcdata, block_offset + delta_offset, delta_bits);
222 :
223 : /* Should not happen as delta_075_by_level_by_bc[level_index] == NULL if and
224 : only if level_index == 0, which means that pixel_index == 0, which means
225 : (i, j) = (0, 0). That cannot happen as we are never called with those
226 : values
227 : */
228 0 : CPLAssert(delta_075_by_level_by_bc[level_index] != nullptr);
229 0 : const int *lookup_table = delta_075_by_level_by_bc[level_index][busy_code];
230 :
231 0 : CPLAssert(lookup_table != nullptr);
232 0 : int delta = lookup_table[delta_raw];
233 :
234 0 : return delta;
235 : }
236 :
237 : /************************************************************************/
238 : /* decode_block() */
239 : /* */
240 : /* Decode one 8x8 block. The 9x9 L buffer is pre-loaded with */
241 : /* the left and top values from previous blocks. */
242 : /************************************************************************/
243 0 : static int decode_block(unsigned char *srcdata, int nInputBytes, int busy_code,
244 : int comrat, int block_offset, int block_size,
245 : int left_side, int top_side, int L[9][9])
246 :
247 : {
248 : int bError;
249 :
250 : // Level 2
251 0 : L[0][4] = (L[0][0] + L[0][8]) / 2 +
252 0 : get_delta(srcdata, nInputBytes, busy_code, comrat, block_offset,
253 : block_size, 0, 4, &bError);
254 0 : if (bError)
255 0 : return FALSE;
256 0 : L[4][0] = (L[0][0] + L[8][0]) / 2 +
257 0 : get_delta(srcdata, nInputBytes, busy_code, comrat, block_offset,
258 : block_size, 4, 0, &bError);
259 0 : if (bError)
260 0 : return FALSE;
261 0 : L[4][4] = (L[0][0] + L[8][0] + L[0][8] + L[8][8]) / 4 +
262 0 : get_delta(srcdata, nInputBytes, busy_code, comrat, block_offset,
263 : block_size, 4, 4, &bError);
264 0 : if (bError)
265 0 : return FALSE;
266 :
267 0 : if (left_side)
268 0 : L[4][8] = L[4][0];
269 0 : if (top_side)
270 0 : L[8][4] = L[0][4];
271 :
272 : // Level 3
273 0 : for (int i = 0; i < 8; i += 4)
274 : {
275 0 : for (int j = 0; j < 8; j += 4)
276 : {
277 : // above
278 0 : L[i + 2][j] =
279 0 : (L[i][j] + L[i + 4][j]) / 2 +
280 0 : get_delta(srcdata, nInputBytes, busy_code, comrat, block_offset,
281 : block_size, i + 2, j, &bError);
282 0 : if (bError)
283 0 : return FALSE;
284 : // left
285 0 : L[i][j + 2] =
286 0 : (L[i][j] + L[i][j + 4]) / 2 +
287 0 : get_delta(srcdata, nInputBytes, busy_code, comrat, block_offset,
288 : block_size, i, j + 2, &bError);
289 0 : if (bError)
290 0 : return FALSE;
291 : // up-left
292 0 : L[i + 2][j + 2] =
293 0 : (L[i][j] + L[i][j + 4] + L[i + 4][j] + L[i + 4][j + 4]) / 4 +
294 0 : get_delta(srcdata, nInputBytes, busy_code, comrat, block_offset,
295 : block_size, i + 2, j + 2, &bError);
296 0 : if (bError)
297 0 : return FALSE;
298 : }
299 : }
300 :
301 0 : if (left_side)
302 : {
303 0 : L[2][8] = L[2][0];
304 0 : L[6][8] = L[6][0];
305 : }
306 0 : if (top_side)
307 : {
308 0 : L[8][2] = L[0][2];
309 0 : L[8][6] = L[0][6];
310 : }
311 :
312 : // Level 4
313 0 : for (int i = 0; i < 8; i += 2)
314 : {
315 0 : for (int j = 0; j < 8; j += 2)
316 : {
317 : // above
318 0 : L[i + 1][j] =
319 0 : (L[i][j] + L[i + 2][j]) / 2 +
320 0 : get_delta(srcdata, nInputBytes, busy_code, comrat, block_offset,
321 : block_size, i + 1, j, &bError);
322 0 : if (bError)
323 0 : return FALSE;
324 : // left
325 0 : L[i][j + 1] =
326 0 : (L[i][j] + L[i][j + 2]) / 2 +
327 0 : get_delta(srcdata, nInputBytes, busy_code, comrat, block_offset,
328 : block_size, i, j + 1, &bError);
329 0 : if (bError)
330 0 : return FALSE;
331 : // up-left
332 0 : L[i + 1][j + 1] =
333 0 : (L[i][j] + L[i][j + 2] + L[i + 2][j] + L[i + 2][j + 2]) / 4 +
334 0 : get_delta(srcdata, nInputBytes, busy_code, comrat, block_offset,
335 : block_size, i + 1, j + 1, &bError);
336 0 : if (bError)
337 0 : return FALSE;
338 : }
339 : }
340 :
341 0 : return TRUE;
342 : }
343 :
344 : /************************************************************************/
345 : /* NITFUncompressARIDPCM() */
346 : /************************************************************************/
347 :
348 0 : int NITFUncompressARIDPCM(NITFImage *psImage, GByte *pabyInputData,
349 : int nInputBytes, GByte *pabyOutputImage)
350 :
351 : {
352 : /* -------------------------------------------------------------------- */
353 : /* First, verify that we are a COMRAT 0.75 image, which is all */
354 : /* we currently support. */
355 : /* -------------------------------------------------------------------- */
356 0 : if (!EQUAL(psImage->szCOMRAT, "0.75"))
357 : {
358 0 : CPLError(CE_Failure, CPLE_AppDefined,
359 : "COMRAT=%s ARIDPCM is not supported.\n"
360 : "Currently only 0.75 is supported.",
361 0 : psImage->szCOMRAT);
362 0 : return FALSE;
363 : }
364 :
365 : /* -------------------------------------------------------------------- */
366 : /* Setup up the various info we need for each 8x8 neighbourhood */
367 : /* (which we call blocks in this context). */
368 : /* -------------------------------------------------------------------- */
369 0 : const int blocks_x = (psImage->nBlockWidth + 7) / 8;
370 0 : const int blocks_y = (psImage->nBlockHeight + 7) / 8;
371 0 : const int block_count = blocks_x * blocks_y;
372 0 : const int rowlen = blocks_x * 8;
373 :
374 0 : if (psImage->nBlockWidth > 1000 || /* to detect int overflow above */
375 0 : psImage->nBlockHeight > 1000 || block_count > 1000)
376 : {
377 0 : CPLError(CE_Failure, CPLE_AppDefined, "Block too large to be decoded");
378 0 : return FALSE;
379 : }
380 :
381 : int block_offset[1000];
382 : int block_size[1000];
383 : int busy_code[1000];
384 0 : const int busy_code_table_size = blocks_x * blocks_y * 2;
385 : unsigned char L00[1000];
386 :
387 : /* to make clang static analyzer happy */
388 0 : block_offset[0] = 0;
389 0 : block_size[0] = 0;
390 0 : busy_code[0] = 0;
391 0 : L00[0] = 0;
392 0 : CPL_IGNORE_RET_VAL(busy_code[0]);
393 0 : CPL_IGNORE_RET_VAL(block_size[0]);
394 0 : CPL_IGNORE_RET_VAL(busy_code[0]);
395 0 : CPL_IGNORE_RET_VAL(L00[0]);
396 :
397 : /* -------------------------------------------------------------------- */
398 : /* We allocate a working copy of the full image that may be a */
399 : /* bit larger than the output buffer if the width or height is */
400 : /* not divisible by 8. */
401 : /* -------------------------------------------------------------------- */
402 : GByte *full_image = reinterpret_cast<GByte *>(
403 0 : CPLMalloc(static_cast<size_t>(blocks_x) * blocks_y * 8 * 8));
404 :
405 : /* -------------------------------------------------------------------- */
406 : /* Scan through all the neighbourhoods determining the busyness */
407 : /* code, and the offset to each's data as well as the L00 value. */
408 : /* -------------------------------------------------------------------- */
409 0 : int total = busy_code_table_size;
410 :
411 0 : for (int i = 0; i < blocks_x * blocks_y; i++)
412 : {
413 0 : if (nInputBytes * 8 < i * 2 + 2)
414 : {
415 0 : CPLError(CE_Failure, CPLE_AppDefined, "Input buffer too small");
416 0 : CPLFree(full_image);
417 0 : return FALSE;
418 : }
419 0 : busy_code[i] = get_bits(pabyInputData, i * 2, 2);
420 :
421 0 : block_offset[i] = total;
422 0 : block_size[i] = neighbourhood_size_75[busy_code[i]];
423 :
424 0 : if (nInputBytes * 8 < block_offset[i] + 8)
425 : {
426 0 : CPLError(CE_Failure, CPLE_AppDefined, "Input buffer too small");
427 0 : CPLFree(full_image);
428 0 : return FALSE;
429 : }
430 0 : L00[i] = (unsigned char)get_bits(pabyInputData, block_offset[i], 8);
431 :
432 0 : total += block_size[i];
433 : }
434 :
435 : /* -------------------------------------------------------------------- */
436 : /* Process all the blocks, forming into a final image. */
437 : /* -------------------------------------------------------------------- */
438 0 : for (int iY = 0; iY < blocks_y; iY++)
439 : {
440 0 : for (int iX = 0; iX < blocks_x; iX++)
441 : {
442 0 : int iBlock = iX + iY * blocks_x;
443 : int L[9][9];
444 0 : unsigned char *full_tl = full_image + iX * 8 + iY * 8 * rowlen;
445 :
446 0 : L[0][0] = L00[iBlock];
447 0 : if (iX > 0)
448 : {
449 0 : L[0][8] = full_tl[rowlen * 7 - 1];
450 0 : L[2][8] = full_tl[rowlen * 5 - 1];
451 0 : L[4][8] = full_tl[rowlen * 3 - 1];
452 0 : L[6][8] = full_tl[rowlen * 1 - 1];
453 : }
454 : else
455 : {
456 0 : L[0][8] = L[0][0];
457 0 : L[2][8] = L[0][8]; // need to reconstruct the rest!
458 0 : L[4][8] = L[0][8];
459 0 : L[6][8] = L[0][8];
460 : }
461 :
462 0 : if (iY > 0)
463 : {
464 0 : L[8][0] = full_tl[7 - rowlen];
465 0 : L[8][2] = full_tl[5 - rowlen];
466 0 : L[8][4] = full_tl[3 - rowlen];
467 0 : L[8][6] = full_tl[1 - rowlen];
468 : }
469 : else
470 : {
471 0 : L[8][0] = L[0][0];
472 0 : L[8][2] = L[0][0]; // Need to reconstruct the rest!
473 0 : L[8][4] = L[0][0];
474 0 : L[8][5] = L[0][0];
475 : }
476 :
477 0 : if (iX == 0 || iY == 0)
478 0 : L[8][8] = L[0][0];
479 : else
480 0 : L[8][8] = full_tl[-1 - rowlen];
481 :
482 0 : if (!(decode_block(pabyInputData, nInputBytes, busy_code[iBlock],
483 : CR075, block_offset[iBlock], block_size[iBlock],
484 : iX == 0, iY == 0, L)))
485 : {
486 0 : CPLFree(full_image);
487 0 : return FALSE;
488 : }
489 :
490 : // Assign to output matrix.
491 0 : for (int i = 0; i < 8; i++)
492 : {
493 0 : for (int j = 0; j < 8; j++)
494 : {
495 0 : int value = L[i][j];
496 0 : if (value < 0)
497 0 : value = 0;
498 0 : if (value > 255)
499 0 : value = 255;
500 :
501 0 : full_tl[8 - j - 1 + (8 - i - 1) * rowlen] =
502 : (unsigned char)value;
503 : }
504 : }
505 : }
506 : }
507 :
508 : /* -------------------------------------------------------------------- */
509 : /* Copy full image back into target buffer, and free. */
510 : /* -------------------------------------------------------------------- */
511 0 : for (int iY = 0; iY < psImage->nBlockHeight; iY++)
512 : {
513 0 : memcpy(pabyOutputImage + iY * psImage->nBlockWidth,
514 0 : full_image + iY * rowlen, psImage->nBlockWidth);
515 : }
516 :
517 0 : CPLFree(full_image);
518 :
519 0 : return TRUE;
520 : }
|
