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