1 #ifndef NKITKXDRIMAGEIO_CXX
2 #define NKITKXDRIMAGEIO_CXX
5 * @file nkitkXDRImageIO.cxx
6 * @author Simon Rit <simon.rit@gmail.com>
7 * @date Sun Jun 1 22:12:20 2008
14 #include "nkitkXDRImageIO.h"
23 #define AVSerror(v) std::cerr << "Error in nkitk::XDRImageIO. Message:" << v << std::endl;
28 #define memicmp _memicmp
30 #define memicmp strncasecmp
33 /************************************************************************/
35 /* file : AVS_RXDR.CPP */
37 /* purpose : AVS module for reading XDRs */
39 /* authors : Lambert Zijp (Conquest) */
40 /* Based on a true story by Marcel van Herk */
44 /* portability: AVS requires sizeof(void *)==sizeof(int) */
45 /* This module assumes sizeof(int)>=4 */
50 /************************************************************************/
53 19980528 ljz Copied from mbfield1
54 19981209 mvh Taken into action by removing _; added READ_XDR_HEADER
55 19981210 mvh+lsp Now proper freeing of buff in case of error
56 19990130 mvh Added ENUM_XDR_HEADER
57 19990206 mvh Added READ_XDR_PREVIEW
58 19990208 mvh Fixed this, the while(1) fgets loop is rquired to skip long headers
59 19990406 ljz Removed sizelimit of dimensions
60 19991003 mvh Loosened sizelimit for veclen to 1000
61 19991005 mvh Remove spaces from AVS own header lines (ndim etc) in scan_header
62 20000808 lsp No longer mix file streams and handles
63 20000821 ljz Reader can now handle NkiCompression
64 20000822 lsp+nd Initialize data and datasize for dots, added "const" to avoid
65 writing into a literal string
66 20000823 mvh Optimized nki_private_decompression mode 2 (from 6.4 to 4.9 s)
67 reading mode 1: compression makes read faster
68 reading mode 2: read compressed at same speed as uncompressed
69 20000824 mvh Speed up mode 2 to 4.2 s: only check compressedCRC if other failed
70 Added some safeties to avoid mode 1 or 2 crash for corrupted data
71 20000825 mvh Pass buffer size to nki_private_decompress as extra safety
72 Fix mode 1, added full safety against input buffer overflow
73 Added some comments and notes
74 20010507 mvh Support larger headers when coordinate information has been written
75 20010726 mvh Fix for decompression when information is offset in file
76 20011207 ljz Removed check on veclen>1000
77 20020124 mvh Added test for fields by kg: non-portable types as "integer" then not swapped
78 20020903 ljz Made scan_header much faster
79 20030122 mvh Fix read of coords in compression mode 2
80 20030430 mvh Fix of read coords in compression mode 2 when offset is 0
81 20030717 ljz Added support for NkiCompressionModes 3 and 4
82 20040426 mvh ELEKTA NKI-XVI0.1 RELEASE
83 20050302 mvh Estimate size of compressed data to accelerate speed of reading embedded fields
84 20050308 ljz+mvh ELEKTA NKI-XVI0.1j RELEASE
85 20071024 mvh Adapted for 64 bits
86 20080110 mvh ELEKTA NKI-XVI3.09 RELEASE
87 20080414 lsp+mgw __sun__ doesn't know <io.h>
90 /************************************************************************/
91 /* MODULE DOCUMENTATION */
92 /************************************************************************/
94 READ_XDR Read XDR file (may be compressed) into field
95 name output field handle
96 file_expression name of file
97 numerical_expression start XDR data in file offset (default 0)
98 NOTE: compressed XDR data may not be part of a larger file
100 READ_XDR_HEADER Get entry from xdr header
101 %name Output = string block
102 file_expression file name default extension ''
103 %string_expression name of element to load
105 READ_XDR_PREVIEW Read and downsize XDR file (for bitmap)
106 name output field handle
107 file_expression name of file
108 numerical_expression start XDR data in file offset (default 0)
109 NOTE: this command is not supported for compressed XDR files
112 %name Output = string block
113 file_expression file name default extension ''
114 numerical_expression number of element to find name of
116 /************************************************************************/
118 /************************************************************************/
134 #include <avs/field.h>
140 /************************************************************************/
141 /* DEFINES, ENUMERATED TYPES AND CONSTANTS */
142 /************************************************************************/
163 unsigned int iOrgSize;
165 unsigned int iCompressedSize;
166 unsigned int iOrgCRC;
167 unsigned int iCompressedCRC; /* Excluding this header */
171 /************************************************************************/
172 /* GLOBAL VARIABLES */
173 /************************************************************************/
174 const char* gl_ErrorMsg[] = {
176 "Command or function not supported in this way.",
177 "Error in arguments",
178 "XDR file header NDIM error",
179 "XDR file header DIMn error",
180 "XDR file header NSPACE error",
181 "XDR file header VECLEN error",
182 "XDR file header DATA(type) error",
183 "XDR file header FIELD(coordinate type) error",
184 "XDR file could not be opened",
185 "XDR file header contains no ^L",
186 "XDR file reading error",
188 "Decompression failed",
189 "Format not handled by nkitkXDRImageIO (RECTILINEAR or IRREGULAR field)"
193 static const unsigned long CRC32_table[256] = {
194 0x00000000, 0x77073096, 0xee0e612c, 0x990951ba, 0x076dc419, 0x706af48f,
195 0xe963a535, 0x9e6495a3, 0x0edb8832, 0x79dcb8a4, 0xe0d5e91e, 0x97d2d988,
196 0x09b64c2b, 0x7eb17cbd, 0xe7b82d07, 0x90bf1d91, 0x1db71064, 0x6ab020f2,
197 0xf3b97148, 0x84be41de, 0x1adad47d, 0x6ddde4eb, 0xf4d4b551, 0x83d385c7,
198 0x136c9856, 0x646ba8c0, 0xfd62f97a, 0x8a65c9ec, 0x14015c4f, 0x63066cd9,
199 0xfa0f3d63, 0x8d080df5, 0x3b6e20c8, 0x4c69105e, 0xd56041e4, 0xa2677172,
200 0x3c03e4d1, 0x4b04d447, 0xd20d85fd, 0xa50ab56b, 0x35b5a8fa, 0x42b2986c,
201 0xdbbbc9d6, 0xacbcf940, 0x32d86ce3, 0x45df5c75, 0xdcd60dcf, 0xabd13d59,
202 0x26d930ac, 0x51de003a, 0xc8d75180, 0xbfd06116, 0x21b4f4b5, 0x56b3c423,
203 0xcfba9599, 0xb8bda50f, 0x2802b89e, 0x5f058808, 0xc60cd9b2, 0xb10be924,
204 0x2f6f7c87, 0x58684c11, 0xc1611dab, 0xb6662d3d, 0x76dc4190, 0x01db7106,
205 0x98d220bc, 0xefd5102a, 0x71b18589, 0x06b6b51f, 0x9fbfe4a5, 0xe8b8d433,
206 0x7807c9a2, 0x0f00f934, 0x9609a88e, 0xe10e9818, 0x7f6a0dbb, 0x086d3d2d,
207 0x91646c97, 0xe6635c01, 0x6b6b51f4, 0x1c6c6162, 0x856530d8, 0xf262004e,
208 0x6c0695ed, 0x1b01a57b, 0x8208f4c1, 0xf50fc457, 0x65b0d9c6, 0x12b7e950,
209 0x8bbeb8ea, 0xfcb9887c, 0x62dd1ddf, 0x15da2d49, 0x8cd37cf3, 0xfbd44c65,
210 0x4db26158, 0x3ab551ce, 0xa3bc0074, 0xd4bb30e2, 0x4adfa541, 0x3dd895d7,
211 0xa4d1c46d, 0xd3d6f4fb, 0x4369e96a, 0x346ed9fc, 0xad678846, 0xda60b8d0,
212 0x44042d73, 0x33031de5, 0xaa0a4c5f, 0xdd0d7cc9, 0x5005713c, 0x270241aa,
213 0xbe0b1010, 0xc90c2086, 0x5768b525, 0x206f85b3, 0xb966d409, 0xce61e49f,
214 0x5edef90e, 0x29d9c998, 0xb0d09822, 0xc7d7a8b4, 0x59b33d17, 0x2eb40d81,
215 0xb7bd5c3b, 0xc0ba6cad, 0xedb88320, 0x9abfb3b6, 0x03b6e20c, 0x74b1d29a,
216 0xead54739, 0x9dd277af, 0x04db2615, 0x73dc1683, 0xe3630b12, 0x94643b84,
217 0x0d6d6a3e, 0x7a6a5aa8, 0xe40ecf0b, 0x9309ff9d, 0x0a00ae27, 0x7d079eb1,
218 0xf00f9344, 0x8708a3d2, 0x1e01f268, 0x6906c2fe, 0xf762575d, 0x806567cb,
219 0x196c3671, 0x6e6b06e7, 0xfed41b76, 0x89d32be0, 0x10da7a5a, 0x67dd4acc,
220 0xf9b9df6f, 0x8ebeeff9, 0x17b7be43, 0x60b08ed5, 0xd6d6a3e8, 0xa1d1937e,
221 0x38d8c2c4, 0x4fdff252, 0xd1bb67f1, 0xa6bc5767, 0x3fb506dd, 0x48b2364b,
222 0xd80d2bda, 0xaf0a1b4c, 0x36034af6, 0x41047a60, 0xdf60efc3, 0xa867df55,
223 0x316e8eef, 0x4669be79, 0xcb61b38c, 0xbc66831a, 0x256fd2a0, 0x5268e236,
224 0xcc0c7795, 0xbb0b4703, 0x220216b9, 0x5505262f, 0xc5ba3bbe, 0xb2bd0b28,
225 0x2bb45a92, 0x5cb36a04, 0xc2d7ffa7, 0xb5d0cf31, 0x2cd99e8b, 0x5bdeae1d,
226 0x9b64c2b0, 0xec63f226, 0x756aa39c, 0x026d930a, 0x9c0906a9, 0xeb0e363f,
227 0x72076785, 0x05005713, 0x95bf4a82, 0xe2b87a14, 0x7bb12bae, 0x0cb61b38,
228 0x92d28e9b, 0xe5d5be0d, 0x7cdcefb7, 0x0bdbdf21, 0x86d3d2d4, 0xf1d4e242,
229 0x68ddb3f8, 0x1fda836e, 0x81be16cd, 0xf6b9265b, 0x6fb077e1, 0x18b74777,
230 0x88085ae6, 0xff0f6a70, 0x66063bca, 0x11010b5c, 0x8f659eff, 0xf862ae69,
231 0x616bffd3, 0x166ccf45, 0xa00ae278, 0xd70dd2ee, 0x4e048354, 0x3903b3c2,
232 0xa7672661, 0xd06016f7, 0x4969474d, 0x3e6e77db, 0xaed16a4a, 0xd9d65adc,
233 0x40df0b66, 0x37d83bf0, 0xa9bcae53, 0xdebb9ec5, 0x47b2cf7f, 0x30b5ffe9,
234 0xbdbdf21c, 0xcabac28a, 0x53b39330, 0x24b4a3a6, 0xbad03605, 0xcdd70693,
235 0x54de5729, 0x23d967bf, 0xb3667a2e, 0xc4614ab8, 0x5d681b02, 0x2a6f2b94,
236 0xb40bbe37, 0xc30c8ea1, 0x5a05df1b, 0x2d02ef8d
240 /************************************************************************/
242 /************************************************************************/
244 int Rxdr_compute(AVSfield** ppOut, char* pszFileName, int iOffset);
246 int RxdrPreview_compute(AVSfield** ppOut, char* pszFileName, int iOffset);
248 int RxdrHeader_compute(char *pszOut, char *pszFileName, char *pszEntry);
250 int RxdrEnum_compute(char *pszOut, char *pszFileName, int iEntry);
252 /******************************************************************************/
253 /* AVS/QUIRT INTERFACE */
254 /******************************************************************************/
260 /* Set the module name and type */
261 AVSset_module_name("XDR reader", MODULE_DATA);
263 /* Create output ports for the resulting fields */
264 AVScreate_output_port("Output", "field");
266 /* declare widgets */
267 QUIRT_NEXT_PARAMETER_FILE("");
268 param = AVSadd_parameter("FileName", "string", "/data/AVS/", NULL, "");
269 AVSconnect_widget(param, "browser");
270 AVSadd_parameter_prop(param, "height", "integer", 8);
272 param = AVSadd_parameter("Offset", "integer", 0, 0, INT_UNBOUND);
273 AVSconnect_widget(param, "typein_integer");
275 AVSset_compute_proc(CF Rxdr_compute);
277 AVS_TO_QUIRT(READ_XDR, Rxdr_desc);
279 void RxdrPreview_desc(void)
283 /* Set the module name and type */
284 AVSset_module_name("XDR preview reader", MODULE_DATA);
286 /* Create output ports for the resulting fields */
287 AVScreate_output_port("Output", "field");
289 /* declare widgets */
290 QUIRT_NEXT_PARAMETER_FILE("");
291 param = AVSadd_parameter("FileName", "string", "/data/AVS/", NULL, "");
292 AVSconnect_widget(param, "browser");
293 AVSadd_parameter_prop(param, "height", "integer", 8);
295 param = AVSadd_parameter("Offset", "integer", 0, 0, INT_UNBOUND);
296 AVSconnect_widget(param, "typein_integer");
298 AVSset_compute_proc(CF RxdrPreview_compute);
300 AVS_TO_QUIRT(READ_XDR_PREVIEW, RxdrPreview_desc);
302 void RxdrHeader_desc(void)
303 { /* Set the module name and type */
304 AVSset_module_name("Read XDR header", MODULE_DATA);
306 /* Create output ports for the resulting fields */
307 AVSadd_parameter("Output", "string_block", "", NULL, "");
309 QUIRT_NEXT_PARAMETER_FILE("");
310 AVSadd_parameter("FileName", "string", "", NULL, "");
312 AVSadd_parameter("Entry", "string", "", NULL, "");
314 AVSset_compute_proc(CF RxdrHeader_compute);
316 AVS_TO_QUIRT(READ_XDR_HEADER, RxdrHeader_desc);
318 void RxdrEnum_desc(void)
319 { /* Set the module name and type */
320 AVSset_module_name("Enumerate XDR header", MODULE_DATA);
322 /* Create output ports for the resulting fields */
323 AVSadd_parameter("Output", "string_block", "", NULL, "");
325 QUIRT_NEXT_PARAMETER_FILE("");
326 AVSadd_parameter("FileName", "string", "", NULL, "");
328 AVSadd_parameter("Entry", "integer", 0, INT_UNBOUND, INT_UNBOUND);
330 AVSset_compute_proc(CF RxdrEnum_compute);
332 AVS_TO_QUIRT(ENUM_XDR_HEADER, RxdrEnum_desc);
335 AVSinit_modules(void)
337 AVSmodule_from_desc( (int_desc_func)Rxdr_desc );
338 AVSmodule_from_desc( (int_desc_func)RxdrHeader_desc );
339 AVSmodule_from_desc( (int_desc_func)RxdrEnum_desc );
343 /************************************************************************/
344 /* MODULE FUNCTIONS */
345 /************************************************************************/
347 /* simple XDR file reader */
349 /* help routine, scans XDR file header for keyword entry, returns NULL
350 if not found, else returns pointer to rest of line
353 static char *scan_header(const char *file, const char *name, int offset, int removespaces)
355 int i, j, iStringLength;
356 static char temp[512];
360 if ((f = fopen(file, "rt")) == NULL) return NULL;
361 if (offset) fseek(f, offset, SEEK_SET);
365 if (fgets(temp, 500, f) == NULL ) break; /* end of file */
370 p = q = temp; /* remove spaces */
371 iStringLength = strlen(temp);
372 for (j=0; j<iStringLength; j++)
373 if (*q!=' ' && *q!=8) *p++ = *q++;
378 if (temp[0] == 12 ) break; /* ^L end of header */
379 if (temp[0] != '#') i++; /* The first 200 non comment lines must be read before data is opened. */
380 if ((p = strchr(temp+1, '=')) == NULL) continue; /* no '=' */
381 if (memicmp(temp, name, p-temp) ) continue; /* no match */
383 p++; /* match, skip = */
384 if (p[strlen(p)-1] == '\n') /* remove \n */
395 /* help routine, enumerates XDR file for names of keyword entrys, returns NULL
396 if not found, else returns pointer to start of line upto '='
398 // //Commented out because it does not seem to be needed for vv
399 //static char *enum_header(char *file, int iEntry, int offset)
402 // static char temp[512];
406 // if ((f = fopen(file, "rt")) == NULL) return NULL;
407 // if (offset) fseek(f, offset, SEEK_SET);
411 // for (i=0; i<200; )
413 // if (fgets(temp, 500, f) == NULL ) break; /* end of file */
415 // if (temp[0] == 12 ) break; /* ^L end of header */
416 // if (temp[0] != '#') i++; /* The first 200 non comment lines must be read before data is opened. */
417 // if ((p = strchr(temp+1, '=')) == NULL) continue; /* no '=' */
418 // if (count++ != iEntry) continue; /* no match */
420 // *p=0; /* match, end at = */
431 static int get_nki_compressed_size(FILE *f)
436 fread((void *)&Header, sizeof(Header), 1 , f);
438 iMode = Header.iMode;
447 return Header.iCompressedSize + sizeof(Header);
453 /* decoder for NKI private compressed pixel data
454 arguments: dest = (in) points to area where destination data is written (short)
455 src = (in) points compressed source data (byte stream)
456 size = (in) number of bytes source data (safety)
458 The return value is the number of pixels that have been processed.
460 The compressed data looks like:
461 (number of pixels)-1 times:
462 OR 1 byte = LL 7 bits signed (difference pixel[1] - pixel[0]);
463 OR 2 bytes = HHLL 15 bits signed (difference pixel[1] - pixel[0]) xored with 0x4000;
464 OR 3 bytes = 7FHHLL 16 bits absolute pixel data if 15 bits difference is exceeded
465 OR 2 bytes = 80NN run length encode NN zero differences (max 255)
466 for mode 3 and 4 added:
467 OR 2 bytes = CONN encode NN 4 bit differences (max 255)
469 Performance on typical CT or MRI is >2x compression and a very good speed
471 This code is not valid on HIGHENDIAN (high byte first) machines
474 static int global_len;
476 static int nki_private_decompress(short int *dest, signed char *src, int size)
478 int npixels, retvalue, mode, iMode, val, j;
479 NKI_MODE2* pHeader = (NKI_MODE2*)src;
480 unsigned long iCRC=0, iCRC2=0;
481 //unsigned char* pDestStart = (unsigned char*)dest;
482 signed char *save, *end;
484 retvalue = npixels = pHeader->iOrgSize;
485 iMode = pHeader->iMode; // safety: this value is checked in case statement
487 if (npixels<1) return 0; // safety: check for invalid npixels value
489 /* Up till now only Mode=1, 2, 3, and 4 are supported */
496 src += 8; // mode 1 only has 8 bytes header: iOrgSize and iMode
497 end = src + size - 3; // for overflow check if we are close to end of input buffer
499 *dest = *(short int *)src;
505 if (src > end) // check whether the last few messages fit in input buffer
507 if (src<end+3) val = *src;
510 if (val >= -64 && val <= 63) mode = 1; // 7 bit difference
511 else if (val==0x7f) mode = 3; // 16 bit value
512 else if ((val&0xff)==0x80) mode = 2; // run length encoding
515 if (src+mode > end+3)
516 return 0; // safety: overflow input data
521 if (val >= -64 && val <= 63) // 7 bit difference
523 dest[1] = dest[0] + val;
527 else if (val==0x7f) // 16 bit value
529 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
533 else if ((val&0xff)==0x80) // run length encoding
535 mode = ((unsigned char *)src)[1];
537 if (npixels<=0) return 0; // safety: overflow output data
548 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
549 dest[1] = dest[0] + diff; // 15 bit difference
556 global_len = src-save;
561 src += sizeof(NKI_MODE2);
563 end = src + pHeader->iCompressedSize - 3;
565 if (end > src + size - 3)
566 end = src + size - 3; // may occur if pHeader is corrupted
568 *dest = val = *(short int *)src;
569 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
570 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
577 if (src > end) // check whether the last few messages fit in input buffer
579 if (src<end+3) val = *src;
582 if (val >= -64 && val <= 63) mode = 1; // 7 bit difference
583 else if (val==0x7f) mode = 3; // 16 bit value
584 else if ((val&0xff)==0x80) mode = 2; // run length encoding
587 if (src+mode > end+3)
588 break; // safety: overflow input data
593 if (val >= -64 && val <= 63) // 7 bits difference
595 dest[1] = val = dest[0] + val;
596 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
597 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
601 else if (val==0x7f) // 16 bit value
603 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
605 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
606 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
610 else if ((val&0xff)==0x80) // run length encoding
612 mode = ((unsigned char *)src)[1];
614 if (npixels<=0) break; // safety: overflow output data
617 dest[1] = val = dest[0];
618 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
619 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
627 signed short diff = ((val^0x40)<<8) + ((unsigned char *)src)[1];
628 dest[1] = val = dest[0] + diff; // 15 bit difference
629 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
630 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
637 if (iCRC2 != pHeader->iOrgCRC) // if error in output CRC:
639 src = save; // check input CRC
642 iCRC = CRC32_table[(unsigned char)iCRC ^ (unsigned char)src[0]] ^ ((iCRC >> 8));
646 if (iCRC != pHeader->iCompressedCRC)
648 AVSerror("XDR decompression: the file is corrupted");
653 AVSerror("XDR decompression: internal error");
658 global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
665 src += 8; // mode 3 only has 8 bytes header: iOrgSize and iMode
666 end = src + size - 3; // for overflow check if we are close to end of input buffer
668 *dest = *(short int *)src;
674 if (src > end) // check whether the last few messages fit in input buffer
676 if (src<end+3) val = *src;
679 if (val >= -63 && val <= 63) mode = 1; // 7 bit difference
680 else if (val==0x7f) mode = 3; // 16 bit value
681 else if ((val&0xff)==0x80) mode = 2; // run length encoding
682 else if ((val&0xff)==0xC0) mode = 2; // 4 bit encoding
685 if (src+mode > end+3)
686 return 0; // safety: overflow input data
691 if (val >= -63 && val <= 63) // 7 bit difference
693 dest[1] = dest[0] + val;
697 else if (val==0x7f) // 16 bit value
699 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
703 else if ((val&0xff)==0x80) // run length encoding
705 mode = ((unsigned char *)src)[1];
707 if (npixels<=0) return 0; // safety: overflow output data
716 else if ((val&0xff)==0xC0) // 4 bit run
718 mode = ((unsigned char *)src)[1];
722 if (npixels<=0) return 0; // safety: overflow output data
726 dest[1] = dest[0] + (val>>4);
728 if (val&8) val |= 0xfffffff0;
730 dest[1] = dest[0] + val;
737 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
738 dest[1] = dest[0] + diff; // 15 bit difference
745 global_len = src-save;
750 src += sizeof(NKI_MODE2);
752 end = src + pHeader->iCompressedSize - 3;
754 if (end > src + size - 3)
755 end = src + size - 3; // may occur if pHeader is corrupted
757 *dest = val = *(short int *)src;
758 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
759 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
765 if (src > end) // check whether the last few messages fit in input buffer
767 if (src<end+3) val = *src;
770 if (val >= -63 && val <= 63) mode = 1; // 7 bit difference
771 else if (val==0x7f) mode = 3; // 16 bit value
772 else if ((val&0xff)==0x80) mode = 2; // run length encoding
773 else if ((val&0xff)==0xC0) mode = 2; // 4 bit encoding
776 if (src+mode > end+3)
777 return 0; // safety: overflow input data
782 if (val >= -63 && val <= 63) // 7 bit difference
784 dest[1] = val = dest[0] + val;
785 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
786 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
790 else if (val==0x7f) // 16 bit value
792 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
793 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
794 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
798 else if ((val&0xff)==0x80) // run length encoding
800 mode = ((unsigned char *)src)[1];
802 if (npixels<=0) return 0; // safety: overflow output data
805 dest[1] = val = dest[0];
806 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
807 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
813 else if ((val&0xff)==0xC0) // 4 bit run
815 mode = ((unsigned char *)src)[1];
819 if (npixels<=0) return 0; // safety: overflow output data
823 dest[1] = j = dest[0] + (val>>4);
824 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
825 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
827 if (val&8) val |= 0xfffffff0;
829 dest[1] = j = dest[0] + val;
830 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
831 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
838 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
839 dest[1] = val = dest[0] + diff; // 15 bit difference
840 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
841 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
848 if (iCRC2 != pHeader->iOrgCRC) // if error in output CRC:
851 global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
857 AVSerror("XDR decompression: unsupported mode");
865 //====================================================================
866 // Read image information (copied from XDRreader)
867 int nkitk::XDRImageIO::ReadImageInformationWithError()
870 itk::Vector<int,MAXDIM> dim;
871 int veclen=1; //, data=AVS_TYPE_BYTE, field=UNIFORM;
873 int iNkiCompression = 0;
875 int total=1/*, datasize=0, iNumRead, HeaderSize*/;
876 unsigned int coords=0,i,j,ndim,nspace;
882 AVSfield FieldTemplate;
884 long swap_test = 0x1000000; /* first byte is 1 when low-endian */
886 char *file = const_cast<char *>(m_FileName.c_str());
887 AVSType field=UNIFORM;
890 fstream = fopen(file, "rt");
891 if (fstream == NULL) return ER_XDR_OPEN;
893 fgets(temp, 500, fstream);
896 if (memcmp(temp, "# AVS field file (produced by avs_nfwrite.c)", 44)==0) forcenoswap=1;
898 c = scan_header(file, "ndim", offset, 1);
899 if (!c) return ER_XDR_NDIM;
902 if (ndim<1 || ndim>MAXDIM) return ER_XDR_NDIM;
903 SetNumberOfDimensions(ndim);
907 for (i=0; i<ndim; i++)
909 sprintf(temp, "dim%d", i+1);
910 c = scan_header(file, temp, offset, 1);
911 if (!c) return ER_XDR_DIM;
913 if (dim[i]<1) return ER_XDR_DIM;
918 for (i=0; i<ndim; i++) {
919 SetDimensions(i,dim[i]);
924 c = scan_header(file, "nspace", offset, 1);
925 if (c) nspace = atoi(c);
926 if (nspace<1 || ndim > MAXDIM) return ER_XDR_NSPACE;
927 if (nspace != ndim) return ER_NOT_HANDLED;
929 c = scan_header(file, "veclen", offset, 1);
930 if (c) veclen = atoi(c);
931 if (veclen<0 /*|| veclen>1000*/) return ER_XDR_VECLEN;
932 SetNumberOfComponents(veclen);
933 if (veclen==1) SetPixelType(itk::ImageIOBase::SCALAR);
934 else SetPixelType(itk::ImageIOBase::VECTOR);
936 c = scan_header(file, "data", offset, 1);
939 if (memicmp(c, "byte", 4) == 0 || memicmp(c, "xdr_byte", 8) == 0) SetComponentType(itk::ImageIOBase::CHAR);
940 else if (memicmp(c, "short", 5) == 0 || memicmp(c, "xdr_short", 9) == 0) SetComponentType(itk::ImageIOBase::SHORT);
941 else if (memicmp(c, "int" , 3) == 0 || memicmp(c, "xdr_int" , 7) == 0) SetComponentType(itk::ImageIOBase::INT);
942 else if (memicmp(c, "real", 4) == 0 || memicmp(c, "xdr_real", 8) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
943 else if (memicmp(c, "float", 5) == 0 || memicmp(c, "xdr_float", 9) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
944 else if (memicmp(c, "double",6) == 0 || memicmp(c, "xdr_double",10)== 0) SetComponentType(itk::ImageIOBase::DOUBLE);
945 else return ER_XDR_DATA;
947 if (memicmp(c, "xdr_", 4) == 0) forcenoswap=0;
951 c = scan_header(file, "field", offset, 1);
954 if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace *2;
955 else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
956 else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
957 else return ER_XDR_FIELD;
962 if (coords) /* expect AVS coordinates ? */
964 coords *= sizeof(float);
965 fstream = fopen(m_FileName.c_str(), "rb");
966 if (fstream == NULL) return ER_XDR_OPEN;
968 float *points = (float *)malloc(coords);
969 if (points == NULL) return ER_OUTOFMEMORY;
971 //Seek to coordinates position in file
972 if (fseek(fstream,-static_cast<int>(coords),SEEK_END)) return ER_XDR_READ;
973 if (fread( /*(*output)->*/points, 1, coords, fstream ) == coords)
974 { /* swap data if read-ok and required (xdr is low-endian) */
975 if (!(*(char *)(&swap_test)) && !forcenoswap)
977 c = (char *)/*(*output)->*/points;
978 for (i=0; i<coords; i+=4)
992 for (i=0; i<GetNumberOfDimensions(); i++) {
993 SetSpacing(i,10.*(points[i*2+1]-points[i*2])/(GetDimensions(i)-1));
994 SetOrigin(i,10.*points[i*2]);
998 //Rectilinear is reinterpreted as uniform because ITK does not know rectilinear
1000 for (i=0; i<GetNumberOfDimensions(); i++) {
1001 //Compute mean spacing
1002 SetSpacing(i,10*(points[GetDimensions(i)-1]-points[0])/(GetDimensions(i)-1));
1003 SetOrigin(i,10*points[0]);
1005 //Test if rectilinear image is actually uniform (tolerance 0.1 mm)
1006 for (j=0; j<GetDimensions(i)-1; j++) {
1007 if (fabs((points[j+1]-points[j])*10-GetSpacing(i))>0.1) {
1010 return ER_NOT_HANDLED;
1013 points += (int)GetDimensions(i);
1015 for (i=0; i<GetNumberOfDimensions(); i++)
1016 points -= GetDimensions(i);
1021 return ER_NOT_HANDLED;
1029 //====================================================================
1030 // Read image information (copied from XDRreader)
1031 void nkitk::XDRImageIO::ReadImageInformation() {
1032 int result = ReadImageInformationWithError();
1033 if (result) ITKError("nkitk::XDRImageIO::ReadImageInformation",result);
1037 //====================================================================
1038 // Read Image Content (copied from XDRreader)
1039 int nkitk::XDRImageIO::ReadWithError(void * buffer)
1041 int /*ndim,*/ nspace/*, veclen=1, data=AVS_TYPE_BYTE, field=UNIFORM*/;
1042 int iNkiCompression = 0;
1043 int j, coords=0, datasize=0, HeaderSize;
1044 unsigned int i,iNumRead,total=1;
1049 //AVSfield FieldTemplate;
1050 long swap_test = 0x1000000; /* first byte is 1 when low-endian */
1051 //int forcenoswap=0;
1052 char *file = const_cast<char *>(m_FileName.c_str());
1054 AVSType field=UNIFORM;
1056 for (i=0; i<GetNumberOfDimensions(); i++) coords += GetDimensions(i);
1058 total = GetImageSizeInPixels();
1059 nspace = GetNumberOfDimensions();
1061 c = scan_header(file, "field", offset, 1);
1064 if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace*2;
1065 else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
1066 else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
1067 else return ER_XDR_FIELD;
1072 c = scan_header(file, "nki_compression", offset, 1);
1073 if (c) iNkiCompression = atoi(c);
1075 c = scan_header(file, "coord1[0]", offset, 1);
1076 if (c) HeaderSize = 32768;
1077 else HeaderSize = 2048;
1080 FIELDdefault(&FieldTemplate);
1081 FieldTemplate.ndim = ndim;
1082 FieldTemplate.nspace = nspace;
1083 FieldTemplate.veclen = veclen;
1084 FieldTemplate.type = data;
1085 FieldTemplate.uniform = field;
1086 FieldTemplate.size = datasize;
1089 fstream = fopen(file, "rb");
1090 if (fstream == NULL)
1093 if (offset) fseek(fstream, offset, SEEK_SET);
1097 if (fgets(temp, 500, fstream) == NULL )
1098 return ER_XDR_NOCTRLL; /* end of file */
1100 if (temp[0] == 10) continue;
1104 fseek(fstream, -2, SEEK_CUR);
1106 } /* ^L end of header */
1108 if (temp[0] != '#') break;
1111 buff = (char*)malloc(HeaderSize);
1114 return ER_OUTOFMEMORY;
1116 memset(buff, 0, HeaderSize);
1117 iNumRead = fread(buff, 1, HeaderSize, fstream);
1125 for (i=0; i<iNumRead; i++) {
1126 if (buff[i] == 12) break;
1131 if (i==iNumRead) return ER_XDR_NOCTRLL;
1134 if (*output) AVSfield_free(*output);
1135 *output = (AVSfield *) AVSfield_alloc(&FieldTemplate, dim);
1136 if (*output == NULL)
1139 return ER_OUTOFMEMORY;
1142 total *= datasize * veclen;
1143 coords *= sizeof(float);
1145 total = GetImageSizeInBytes();
1147 //We add casts because the resulting quantity can be negative.
1148 //There is no risk of looping because i and iNumRead are about the size of the header
1149 fseek(fstream, static_cast<int>(i)+2-static_cast<int>(iNumRead), SEEK_CUR);
1151 if (total && iNkiCompression)
1154 unsigned long iSize;
1155 signed char* pCompressed;
1157 /* Read or guess the size of the compressed data */
1158 iCurPos = ftell(fstream);
1159 iSize = get_nki_compressed_size(fstream);
1163 fseek(fstream, 0, SEEK_END);
1164 iSize = ftell(fstream);
1165 iSize = iSize - iCurPos - coords;
1167 // Get compressed size from header if possible; else use uncompressed size as safe estimate
1168 if (iSize>total && offset) iSize=total+8;
1171 fseek(fstream, iCurPos, SEEK_SET);
1173 /* Allocate space for the compressed pixels */
1174 pCompressed = (signed char*)malloc(iSize);
1179 if (*output) AVSfield_free(*output);
1182 return ER_OUTOFMEMORY;
1185 /* Read the compressed pixels */
1186 if (fread( (void *)pCompressed, 1, iSize, fstream ) != iSize)
1190 if (*output) AVSfield_free(*output);
1196 if (!nki_private_decompress((short*)buffer, pCompressed, iSize))
1200 if (*output) AVSfield_free(*output);
1203 return ER_DECOMPRESSION;
1207 fseek(fstream, iCurPos + global_len, SEEK_SET);
1216 if (fread( (void *)buffer, 1, total, fstream ) != total)
1220 if (*output) AVSfield_free(*output);
1227 /* swap data if required (xdr is low-endian) */
1229 datasize = GetComponentSize();
1230 if (!(*(char *)(&swap_test)) && !forcenoswap)
1235 for (i=0; i<total; i+=2)
1242 else if (datasize==4)
1245 for (i=0; i<total; i+=4)
1255 else if (datasize==8)
1258 for (i=0; i<total; i+=8)
1278 if (coords) /* expect AVS coordinates ? */
1280 if (fread( (*output)->points, 1, coords, fstream ) == coords)
1281 { /* swap data if read-ok and required (xdr is low-endian) */
1283 if (!(*(char *)(&swap_test)) && !forcenoswap)
1285 c = (char *)(*output)->points;
1286 for (i=0; i<coords; i+=4)
1303 //====================================================================
1304 // Read image information (copied from XDRreader)
1305 void nkitk::XDRImageIO::Read(void * buffer) {
1306 int result = ReadWithError(buffer);
1307 if (result) ITKError("nkitk::XDRImageIO::Read",result);
1311 static int XDRreader_preview(AVSfield **output, char *file, int offset)
1313 AVSINT dim[MAXDIM], dim2[MAXDIM], ds[MAXDIM];
1314 int ndim, nspace, veclen=1, data=AVS_TYPE_BYTE, field=UNIFORM;
1315 int iNkiCompression = 0;
1316 int total=1, i, j, k, l, m, coords=0, datasize=0, iNumRead, len, start, sliceax;
1321 AVSfield FieldTemplate;
1322 long swap_test = 0x1000000; /* first byte is 1 when low-endian */
1325 fstream = fopen(file, "rt");
1326 if (fstream == NULL) return ER_XDR_OPEN;
1327 fgets(temp, 500, fstream);
1330 if (memcmp(temp, "# AVS field file (produced by avs_nfwrite.c)", 44)==0) forcenoswap=1;
1332 c = scan_header(file, "ndim", offset, 1);
1333 if (!c) return ER_XDR_NDIM;
1335 if (ndim<1 || ndim>MAXDIM) return ER_XDR_NDIM;
1338 /* defaults for dimensions and downsize */
1340 for (i=0; i<MAXDIM; i++)
1341 dim[i] = dim2[i] = ds[i] = 1;
1343 for (i=0; i<ndim; i++)
1345 sprintf(temp, "dim%d", i+1);
1346 c = scan_header(file, temp, offset, 1);
1347 if (!c) return ER_XDR_DIM;
1349 if (dim[i]<1 || dim[i]>200000L) return ER_XDR_DIM;
1355 c = scan_header(file, "nspace", offset, 1);
1356 if (c) nspace = atoi(c);
1357 if (nspace<1 || ndim > MAXDIM) return ER_XDR_NSPACE;
1359 c = scan_header(file, "veclen", offset, 1);
1360 if (c) veclen = atoi(c);
1361 if (veclen<0 || veclen>100) return ER_XDR_VECLEN;
1363 c = scan_header(file, "data", offset, 1);
1367 if (memicmp(c, "byte", 4) == 0) data=AVS_TYPE_BYTE, datasize=1;
1368 else if (memicmp(c, "short", 5) == 0) data=AVS_TYPE_SHORT, datasize=2;
1369 else if (memicmp(c, "int" , 3) == 0) data=AVS_TYPE_INTEGER,datasize=4;
1370 else if (memicmp(c, "real", 4) == 0) data=AVS_TYPE_REAL, datasize=4;
1371 else if (memicmp(c, "float", 5) == 0) data=AVS_TYPE_REAL, datasize=4;
1372 else if (memicmp(c, "double",6) == 0) data=AVS_TYPE_DOUBLE, datasize=8;
1374 else if (memicmp(c, "xdr_byte", 8) == 0) data=AVS_TYPE_BYTE, datasize=1, forcenoswap=0;
1375 else if (memicmp(c, "xdr_short", 9) == 0) data=AVS_TYPE_SHORT, datasize=2, forcenoswap=0;
1376 else if (memicmp(c, "xdr_int" , 7) == 0) data=AVS_TYPE_INTEGER,datasize=4, forcenoswap=0;
1377 else if (memicmp(c, "xdr_real", 8) == 0) data=AVS_TYPE_REAL, datasize=4, forcenoswap=0;
1378 else if (memicmp(c, "xdr_float", 9) == 0) data=AVS_TYPE_REAL, datasize=4, forcenoswap=0;
1379 else if (memicmp(c, "xdr_double",10)== 0) data=AVS_TYPE_DOUBLE, datasize=8, forcenoswap=0;
1380 else return ER_XDR_DATA;
1383 c = scan_header(file, "field", offset, 1);
1386 if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace*2;
1387 else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
1388 else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
1389 else return ER_XDR_FIELD;
1394 c = scan_header(file, "nki_compression", offset, 1);
1395 if (c) iNkiCompression = atoi(c);
1396 if (iNkiCompression)
1399 return ER_ILLCOMMFUNCT;
1402 FIELDdefault(&FieldTemplate);
1403 FieldTemplate.ndim = ndim;
1404 FieldTemplate.nspace = nspace;
1405 FieldTemplate.veclen = veclen;
1406 FieldTemplate.type = data;
1407 FieldTemplate.uniform = field;
1408 FieldTemplate.size = datasize;
1410 fstream = fopen(file, "rb");
1412 if (fstream == NULL)
1415 buff = (char *)malloc(8192);
1418 return ER_OUTOFMEMORY;
1420 memset(buff, 0, 8192);
1421 fseek(fstream, offset, SEEK_SET);
1425 if (fgets(temp, 500, fstream) == NULL )
1426 return ER_XDR_NOCTRLL; /* end of file */
1428 if (temp[0] == 10) continue;
1432 fseek(fstream, -2, SEEK_CUR);
1434 } /* ^L end of header */
1436 if (temp[0] != '#') break;
1438 start = ftell(fstream);
1440 iNumRead = fread(buff, 1, 8192, fstream);
1448 for (i=0; i<iNumRead; i++) {
1449 if (buff[i] == 12) break;
1454 if (i==iNumRead) return ER_XDR_NOCTRLL;
1458 /* determine slice axis and downsize slice axis to 3, others to 32 */
1464 if (dim[0]==dim[1]) sliceax = 2;
1470 for (i=0; i<ndim; i++)
1472 if (i==sliceax) ds[i] = dim[i] / 3;
1473 else ds[i] = dim[i] / 32;
1474 if (ds[i]==0) ds[i]=1;
1476 dim2[i] = dim[i]/ds[i];
1477 if (dim2[i]==0) dim2[i] = 1;
1482 if (*output) AVSfield_free(*output);
1483 *output = (AVSfield *) AVSfield_alloc(&FieldTemplate, dim2);
1484 if (*output == NULL)
1487 return ER_OUTOFMEMORY;
1490 total *= datasize * veclen;
1491 coords *= sizeof(float);
1493 /* Read and downsize the data */
1497 c = (char *)((*output)->field_data);
1499 buff = (char *)malloc(dim[0] * datasize * veclen);
1501 for (i=0; i<dim2[4]; i++)
1502 for (j=0; j<dim2[3]; j++)
1503 for (k=0; k<dim2[2]; k++)
1504 for (l=0; l<dim2[1]; l++)
1507 len = len * dim[3] + j*ds[3];
1508 len = len * dim[2] + k*ds[2];
1509 len = len * dim[1] + l*ds[1];
1510 len = len * dim[0] * datasize * veclen;
1512 // This is the time-critical statement
1513 fseek(fstream, start + len, SEEK_SET);
1515 fread(buff, 1, dim[0]*datasize*veclen, fstream);
1516 for (m=0; m<dim2[0]; m++)
1518 memcpy(c, buff+m*datasize*veclen*ds[0], datasize*veclen);
1519 c += datasize*veclen;
1526 /* swap data if required (xdr is low-endian) */
1528 if (!(*(char *)(&swap_test)) && !forcenoswap)
1532 c = (char *)(*output)->field_data;
1533 for (i=0; i<total; i+=2)
1540 else if (datasize==4)
1542 c = (char *)(*output)->field_data;
1543 for (i=0; i<total; i+=4)
1553 else if (datasize==8)
1555 c = (char *)(*output)->field_data;
1556 for (i=0; i<total; i+=8)
1577 { if (read( fHandle,(*output)->points, coords ) == coords)
1578 { if (!(*(char *)(&swap_test)) && !forcenoswap)
1579 { c = (char *)(*output)->points;
1580 for (i=0; i<coords; i+=4)
1581 { j = c[i]; c[i] = c[i+3]; c[i+3] = j;
1582 j = c[i+1]; c[i+1] = c[i+2]; c[i+2] = j;
1593 int Rxdr_compute(AVSfield** ppOut, char* pszFileName, int iOffset)
1598 rc = XDRreader(ppOut, pszFileName, iOffset);
1603 strcpy(szMsg, "Avs_rxdr: ");
1604 strcat(szMsg, gl_ErrorMsg[rc]);
1611 int RxdrPreview_compute(AVSfield** ppOut, char* pszFileName, int iOffset)
1616 rc = XDRreader_preview(ppOut, pszFileName, iOffset);
1621 strcpy(szMsg, "Avs_rxdr: ");
1622 strcat(szMsg, gl_ErrorMsg[rc]);
1629 int RxdrHeader_compute(char *pszOut, char *pszFileName, char *pszEntry)
1633 cc = scan_header(pszFileName, pszEntry, 0, 0);
1634 if (cc) strcpy(pszOut, cc);
1640 int RxdrEnum_compute(char *pszOut, char *pszFileName, int iEntry)
1644 cc = enum_header(pszFileName, iEntry, 0);
1645 if (cc) strcpy(pszOut, cc);
1652 //====================================================================
1653 // Read Image Information
1654 bool nkitk::XDRImageIO::CanReadFile(const char* FileNameToRead)
1659 fstream = fopen(FileNameToRead, "rt");
1660 if (fstream == NULL)
1662 fgets(temp, 500, fstream);
1665 if (memcmp(temp, "# AVS", 5)==0)
1671 void nkitk::XDRImageIO::ITKError(std::string funcName, int msgID) {
1672 itkExceptionMacro(<< "Error in " << funcName << ". Message: " << gl_ErrorMsg[msgID]);
1675 #endif /* end #define NKITKXDRIMAGEIO_CXX */