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