2 * @file clitkXdrImageIO.cxx
3 * @author Simon Rit <simon.rit@gmail.com>
4 * @date Sun Jun 1 22:12:20 2008
11 #include "clitkXdrImageIO.h"
24 #define memicmp _memicmp
26 #define memicmp strncasecmp
29 /************************************************************************/
31 /* file : AVS_RXDR.CPP */
33 /* purpose : AVS module for reading XDRs */
35 /* authors : Lambert Zijp (Conquest) */
36 /* Based on a true story by Marcel van Herk */
40 /* portability: AVS requires sizeof(void *)==sizeof(int) */
41 /* This module assumes sizeof(int)>=4 */
46 /************************************************************************/
49 19980528 ljz Copied from mbfield1
50 19981209 mvh Taken into action by removing _; added READ_XDR_HEADER
51 19981210 mvh+lsp Now proper freeing of buff in case of error
52 19990130 mvh Added ENUM_XDR_HEADER
53 19990206 mvh Added READ_XDR_PREVIEW
54 19990208 mvh Fixed this, the while(1) fgets loop is rquired to skip long headers
55 19990406 ljz Removed sizelimit of dimensions
56 19991003 mvh Loosened sizelimit for veclen to 1000
57 19991005 mvh Remove spaces from AVS own header lines (ndim etc) in scan_header
58 20000808 lsp No longer mix file streams and handles
59 20000821 ljz Reader can now handle NkiCompression
60 20000822 lsp+nd Initialize data and datasize for dots, added "const" to avoid
61 writing into a literal string
62 20000823 mvh Optimized nki_private_decompression mode 2 (from 6.4 to 4.9 s)
63 reading mode 1: compression makes read faster
64 reading mode 2: read compressed at same speed as uncompressed
65 20000824 mvh Speed up mode 2 to 4.2 s: only check compressedCRC if other failed
66 Added some safeties to avoid mode 1 or 2 crash for corrupted data
67 20000825 mvh Pass buffer size to nki_private_decompress as extra safety
68 Fix mode 1, added full safety against input buffer overflow
69 Added some comments and notes
70 20010507 mvh Support larger headers when coordinate information has been written
71 20010726 mvh Fix for decompression when information is offset in file
72 20011207 ljz Removed check on veclen>1000
73 20020124 mvh Added test for fields by kg: non-portable types as "integer" then not swapped
74 20020903 ljz Made scan_header much faster
75 20030122 mvh Fix read of coords in compression mode 2
76 20030430 mvh Fix of read coords in compression mode 2 when offset is 0
77 20030717 ljz Added support for NkiCompressionModes 3 and 4
78 20040426 mvh ELEKTA NKI-XVI0.1 RELEASE
79 20050302 mvh Estimate size of compressed data to accelerate speed of reading embedded fields
80 20050308 ljz+mvh ELEKTA NKI-XVI0.1j RELEASE
81 20071024 mvh Adapted for 64 bits
82 20080110 mvh ELEKTA NKI-XVI3.09 RELEASE
83 20080414 lsp+mgw __sun__ doesn't know <io.h>
92 /************************************************************************/
93 /* DEFINES, ENUMERATED TYPES AND CONSTANTS */
94 /************************************************************************/
115 unsigned int iOrgSize;
117 unsigned int iCompressedSize;
118 unsigned int iOrgCRC;
119 unsigned int iCompressedCRC; /* Excluding this header */
123 /************************************************************************/
124 /* GLOBAL VARIABLES */
125 /************************************************************************/
126 const char* gl_ErrorMsg[] = {
128 "Command or function not supported in this way.",
129 "Error in arguments",
130 "XDR file header NDIM error",
131 "XDR file header DIMn error",
132 "XDR file header NSPACE error",
133 "XDR file header VECLEN error",
134 "XDR file header DATA(type) error",
135 "XDR file header FIELD(coordinate type) error",
136 "XDR file could not be opened",
137 "XDR file header contains no ^L",
138 "XDR file reading error",
140 "Decompression failed",
141 "Format not handled by clitkXdrImageIO (RECTILINEAR or IRREGULAR field)"
145 static const unsigned long CRC32_table[256] = {
146 0x00000000, 0x77073096, 0xee0e612c, 0x990951ba, 0x076dc419, 0x706af48f,
147 0xe963a535, 0x9e6495a3, 0x0edb8832, 0x79dcb8a4, 0xe0d5e91e, 0x97d2d988,
148 0x09b64c2b, 0x7eb17cbd, 0xe7b82d07, 0x90bf1d91, 0x1db71064, 0x6ab020f2,
149 0xf3b97148, 0x84be41de, 0x1adad47d, 0x6ddde4eb, 0xf4d4b551, 0x83d385c7,
150 0x136c9856, 0x646ba8c0, 0xfd62f97a, 0x8a65c9ec, 0x14015c4f, 0x63066cd9,
151 0xfa0f3d63, 0x8d080df5, 0x3b6e20c8, 0x4c69105e, 0xd56041e4, 0xa2677172,
152 0x3c03e4d1, 0x4b04d447, 0xd20d85fd, 0xa50ab56b, 0x35b5a8fa, 0x42b2986c,
153 0xdbbbc9d6, 0xacbcf940, 0x32d86ce3, 0x45df5c75, 0xdcd60dcf, 0xabd13d59,
154 0x26d930ac, 0x51de003a, 0xc8d75180, 0xbfd06116, 0x21b4f4b5, 0x56b3c423,
155 0xcfba9599, 0xb8bda50f, 0x2802b89e, 0x5f058808, 0xc60cd9b2, 0xb10be924,
156 0x2f6f7c87, 0x58684c11, 0xc1611dab, 0xb6662d3d, 0x76dc4190, 0x01db7106,
157 0x98d220bc, 0xefd5102a, 0x71b18589, 0x06b6b51f, 0x9fbfe4a5, 0xe8b8d433,
158 0x7807c9a2, 0x0f00f934, 0x9609a88e, 0xe10e9818, 0x7f6a0dbb, 0x086d3d2d,
159 0x91646c97, 0xe6635c01, 0x6b6b51f4, 0x1c6c6162, 0x856530d8, 0xf262004e,
160 0x6c0695ed, 0x1b01a57b, 0x8208f4c1, 0xf50fc457, 0x65b0d9c6, 0x12b7e950,
161 0x8bbeb8ea, 0xfcb9887c, 0x62dd1ddf, 0x15da2d49, 0x8cd37cf3, 0xfbd44c65,
162 0x4db26158, 0x3ab551ce, 0xa3bc0074, 0xd4bb30e2, 0x4adfa541, 0x3dd895d7,
163 0xa4d1c46d, 0xd3d6f4fb, 0x4369e96a, 0x346ed9fc, 0xad678846, 0xda60b8d0,
164 0x44042d73, 0x33031de5, 0xaa0a4c5f, 0xdd0d7cc9, 0x5005713c, 0x270241aa,
165 0xbe0b1010, 0xc90c2086, 0x5768b525, 0x206f85b3, 0xb966d409, 0xce61e49f,
166 0x5edef90e, 0x29d9c998, 0xb0d09822, 0xc7d7a8b4, 0x59b33d17, 0x2eb40d81,
167 0xb7bd5c3b, 0xc0ba6cad, 0xedb88320, 0x9abfb3b6, 0x03b6e20c, 0x74b1d29a,
168 0xead54739, 0x9dd277af, 0x04db2615, 0x73dc1683, 0xe3630b12, 0x94643b84,
169 0x0d6d6a3e, 0x7a6a5aa8, 0xe40ecf0b, 0x9309ff9d, 0x0a00ae27, 0x7d079eb1,
170 0xf00f9344, 0x8708a3d2, 0x1e01f268, 0x6906c2fe, 0xf762575d, 0x806567cb,
171 0x196c3671, 0x6e6b06e7, 0xfed41b76, 0x89d32be0, 0x10da7a5a, 0x67dd4acc,
172 0xf9b9df6f, 0x8ebeeff9, 0x17b7be43, 0x60b08ed5, 0xd6d6a3e8, 0xa1d1937e,
173 0x38d8c2c4, 0x4fdff252, 0xd1bb67f1, 0xa6bc5767, 0x3fb506dd, 0x48b2364b,
174 0xd80d2bda, 0xaf0a1b4c, 0x36034af6, 0x41047a60, 0xdf60efc3, 0xa867df55,
175 0x316e8eef, 0x4669be79, 0xcb61b38c, 0xbc66831a, 0x256fd2a0, 0x5268e236,
176 0xcc0c7795, 0xbb0b4703, 0x220216b9, 0x5505262f, 0xc5ba3bbe, 0xb2bd0b28,
177 0x2bb45a92, 0x5cb36a04, 0xc2d7ffa7, 0xb5d0cf31, 0x2cd99e8b, 0x5bdeae1d,
178 0x9b64c2b0, 0xec63f226, 0x756aa39c, 0x026d930a, 0x9c0906a9, 0xeb0e363f,
179 0x72076785, 0x05005713, 0x95bf4a82, 0xe2b87a14, 0x7bb12bae, 0x0cb61b38,
180 0x92d28e9b, 0xe5d5be0d, 0x7cdcefb7, 0x0bdbdf21, 0x86d3d2d4, 0xf1d4e242,
181 0x68ddb3f8, 0x1fda836e, 0x81be16cd, 0xf6b9265b, 0x6fb077e1, 0x18b74777,
182 0x88085ae6, 0xff0f6a70, 0x66063bca, 0x11010b5c, 0x8f659eff, 0xf862ae69,
183 0x616bffd3, 0x166ccf45, 0xa00ae278, 0xd70dd2ee, 0x4e048354, 0x3903b3c2,
184 0xa7672661, 0xd06016f7, 0x4969474d, 0x3e6e77db, 0xaed16a4a, 0xd9d65adc,
185 0x40df0b66, 0x37d83bf0, 0xa9bcae53, 0xdebb9ec5, 0x47b2cf7f, 0x30b5ffe9,
186 0xbdbdf21c, 0xcabac28a, 0x53b39330, 0x24b4a3a6, 0xbad03605, 0xcdd70693,
187 0x54de5729, 0x23d967bf, 0xb3667a2e, 0xc4614ab8, 0x5d681b02, 0x2a6f2b94,
188 0xb40bbe37, 0xc30c8ea1, 0x5a05df1b, 0x2d02ef8d
192 /************************************************************************/
193 /* MODULE FUNCTIONS */
194 /************************************************************************/
196 /* simple XDR file reader */
198 /* help routine, scans XDR file header for keyword entry, returns NULL
199 if not found, else returns pointer to rest of line
202 static char *scan_header(const char *file, const char *name, int offset, int removespaces)
204 int i, j, iStringLength;
205 static char temp[512];
209 if ((f = fopen(file, "rt")) == NULL) return NULL;
210 if (offset) fseek(f, offset, SEEK_SET);
214 if (fgets(temp, 500, f) == NULL ) break; /* end of file */
219 p = q = temp; /* remove spaces */
220 iStringLength = strlen(temp);
221 for (j=0; j<iStringLength; j++)
222 if (*q!=' ' && *q!=8) *p++ = *q++;
227 if (temp[0] == 12 ) break; /* ^L end of header */
228 if (temp[0] != '#') i++; /* The first 200 non comment lines must be read before data is opened. */
229 if ((p = strchr(temp+1, '=')) == NULL) continue; /* no '=' */
230 if (memicmp(temp, name, p-temp) ) continue; /* no match */
232 p++; /* match, skip = */
233 if (p[strlen(p)-1] == '\n') /* remove \n */
244 static int get_nki_compressed_size(FILE *f)
249 fread((void *)&Header, sizeof(Header), 1 , f);
251 iMode = Header.iMode;
260 return Header.iCompressedSize + sizeof(Header);
266 /* decoder for NKI private compressed pixel data
267 arguments: dest = (in) points to area where destination data is written (short)
268 src = (in) points compressed source data (byte stream)
269 size = (in) number of bytes source data (safety)
271 The return value is the number of pixels that have been processed.
273 The compressed data looks like:
274 (number of pixels)-1 times:
275 OR 1 byte = LL 7 bits signed (difference pixel[1] - pixel[0]);
276 OR 2 bytes = HHLL 15 bits signed (difference pixel[1] - pixel[0]) xored with 0x4000;
277 OR 3 bytes = 7FHHLL 16 bits absolute pixel data if 15 bits difference is exceeded
278 OR 2 bytes = 80NN run length encode NN zero differences (max 255)
279 for mode 3 and 4 added:
280 OR 2 bytes = CONN encode NN 4 bit differences (max 255)
282 Performance on typical CT or MRI is >2x compression and a very good speed
284 This code is not valid on HIGHENDIAN (high byte first) machines
287 static int global_len;
289 static int nki_private_decompress(short int *dest, signed char *src, int size)
291 int npixels, retvalue, mode, iMode, val, j;
292 NKI_MODE2* pHeader = (NKI_MODE2*)src;
293 unsigned long iCRC=0, iCRC2=0;
294 //unsigned char* pDestStart = (unsigned char*)dest;
295 signed char *save, *end;
297 retvalue = npixels = pHeader->iOrgSize;
298 iMode = pHeader->iMode; // safety: this value is checked in case statement
300 if (npixels<1) return 0; // safety: check for invalid npixels value
302 /* Up till now only Mode=1, 2, 3, and 4 are supported */
309 src += 8; // mode 1 only has 8 bytes header: iOrgSize and iMode
310 end = src + size - 3; // for overflow check if we are close to end of input buffer
312 *dest = *(short int *)src;
318 if (src > end) // check whether the last few messages fit in input buffer
320 if (src<end+3) val = *src;
323 if (val >= -64 && val <= 63) mode = 1; // 7 bit difference
324 else if (val==0x7f) mode = 3; // 16 bit value
325 else if ((val&0xff)==0x80) mode = 2; // run length encoding
328 if (src+mode > end+3)
329 return 0; // safety: overflow input data
334 if (val >= -64 && val <= 63) // 7 bit difference
336 dest[1] = dest[0] + val;
340 else if (val==0x7f) // 16 bit value
342 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
346 else if ((val&0xff)==0x80) // run length encoding
348 mode = ((unsigned char *)src)[1];
350 if (npixels<=0) return 0; // safety: overflow output data
361 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
362 dest[1] = dest[0] + diff; // 15 bit difference
369 global_len = src-save;
374 src += sizeof(NKI_MODE2);
376 end = src + pHeader->iCompressedSize - 3;
378 if (end > src + size - 3)
379 end = src + size - 3; // may occur if pHeader is corrupted
381 *dest = val = *(short int *)src;
382 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
383 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
390 if (src > end) // check whether the last few messages fit in input buffer
392 if (src<end+3) val = *src;
395 if (val >= -64 && val <= 63) mode = 1; // 7 bit difference
396 else if (val==0x7f) mode = 3; // 16 bit value
397 else if ((val&0xff)==0x80) mode = 2; // run length encoding
400 if (src+mode > end+3)
401 break; // safety: overflow input data
406 if (val >= -64 && val <= 63) // 7 bits difference
408 dest[1] = val = dest[0] + val;
409 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
410 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
414 else if (val==0x7f) // 16 bit value
416 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
418 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
419 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
423 else if ((val&0xff)==0x80) // run length encoding
425 mode = ((unsigned char *)src)[1];
427 if (npixels<=0) break; // safety: overflow output data
430 dest[1] = val = dest[0];
431 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
432 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
440 signed short diff = ((val^0x40)<<8) + ((unsigned char *)src)[1];
441 dest[1] = val = dest[0] + diff; // 15 bit difference
442 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
443 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
450 if (iCRC2 != pHeader->iOrgCRC) // if error in output CRC:
452 src = save; // check input CRC
455 iCRC = CRC32_table[(unsigned char)iCRC ^ (unsigned char)src[0]] ^ ((iCRC >> 8));
459 if (iCRC != pHeader->iCompressedCRC)
461 AVSerror("XDR decompression: the file is corrupted");
466 AVSerror("XDR decompression: internal error");
471 global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
478 src += 8; // mode 3 only has 8 bytes header: iOrgSize and iMode
479 end = src + size - 3; // for overflow check if we are close to end of input buffer
481 *dest = *(short int *)src;
487 if (src > end) // check whether the last few messages fit in input buffer
489 if (src<end+3) val = *src;
492 if (val >= -63 && val <= 63) mode = 1; // 7 bit difference
493 else if (val==0x7f) mode = 3; // 16 bit value
494 else if ((val&0xff)==0x80) mode = 2; // run length encoding
495 else if ((val&0xff)==0xC0) mode = 2; // 4 bit encoding
498 if (src+mode > end+3)
499 return 0; // safety: overflow input data
504 if (val >= -63 && val <= 63) // 7 bit difference
506 dest[1] = dest[0] + val;
510 else if (val==0x7f) // 16 bit value
512 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
516 else if ((val&0xff)==0x80) // run length encoding
518 mode = ((unsigned char *)src)[1];
520 if (npixels<=0) return 0; // safety: overflow output data
529 else if ((val&0xff)==0xC0) // 4 bit run
531 mode = ((unsigned char *)src)[1];
535 if (npixels<=0) return 0; // safety: overflow output data
539 dest[1] = dest[0] + (val>>4);
541 if (val&8) val |= 0xfffffff0;
543 dest[1] = dest[0] + val;
550 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
551 dest[1] = dest[0] + diff; // 15 bit difference
558 global_len = src-save;
563 src += sizeof(NKI_MODE2);
565 end = src + pHeader->iCompressedSize - 3;
567 if (end > src + size - 3)
568 end = src + size - 3; // may occur if pHeader is corrupted
570 *dest = val = *(short int *)src;
571 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
572 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
578 if (src > end) // check whether the last few messages fit in input buffer
580 if (src<end+3) val = *src;
583 if (val >= -63 && val <= 63) mode = 1; // 7 bit difference
584 else if (val==0x7f) mode = 3; // 16 bit value
585 else if ((val&0xff)==0x80) mode = 2; // run length encoding
586 else if ((val&0xff)==0xC0) mode = 2; // 4 bit encoding
589 if (src+mode > end+3)
590 return 0; // safety: overflow input data
595 if (val >= -63 && val <= 63) // 7 bit difference
597 dest[1] = val = dest[0] + val;
598 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
599 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
603 else if (val==0x7f) // 16 bit value
605 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
606 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
607 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
611 else if ((val&0xff)==0x80) // run length encoding
613 mode = ((unsigned char *)src)[1];
615 if (npixels<=0) return 0; // safety: overflow output data
618 dest[1] = val = dest[0];
619 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
620 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
626 else if ((val&0xff)==0xC0) // 4 bit run
628 mode = ((unsigned char *)src)[1];
632 if (npixels<=0) return 0; // safety: overflow output data
636 dest[1] = j = dest[0] + (val>>4);
637 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
638 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
640 if (val&8) val |= 0xfffffff0;
642 dest[1] = j = dest[0] + val;
643 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
644 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
651 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
652 dest[1] = val = dest[0] + diff; // 15 bit difference
653 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
654 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
661 if (iCRC2 != pHeader->iOrgCRC) // if error in output CRC:
664 global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
670 AVSerror("XDR decompression: unsupported mode");
678 //====================================================================
679 // Read image information (copied from XDRreader)
680 int clitk::XdrImageIO::ReadImageInformationWithError()
683 itk::Vector<int,MAXDIM> dim;
686 unsigned int coords=0,i,j,ndim,nspace;
691 long swap_test = 0x1000000; /* first byte is 1 when low-endian */
693 char *file = const_cast<char *>(m_FileName.c_str());
694 AVSType field=UNIFORM;
697 fstream = fopen(file, "rt");
698 if (fstream == NULL) return ER_XDR_OPEN;
700 fgets(temp, 500, fstream);
703 if (memcmp(temp, "# AVS field file (produced by avs_nfwrite.c)", 44)==0) forcenoswap=1;
705 c = scan_header(file, "ndim", offset, 1);
706 if (!c) return ER_XDR_NDIM;
709 if (ndim<1 || ndim>MAXDIM) return ER_XDR_NDIM;
710 SetNumberOfDimensions(ndim);
714 for (i=0; i<ndim; i++)
716 sprintf(temp, "dim%d", i+1);
717 c = scan_header(file, temp, offset, 1);
718 if (!c) return ER_XDR_DIM;
720 if (dim[i]<1) return ER_XDR_DIM;
725 for (i=0; i<ndim; i++) {
726 SetDimensions(i,dim[i]);
731 c = scan_header(file, "nspace", offset, 1);
732 if (c) nspace = atoi(c);
733 if (nspace<1 || ndim > MAXDIM) return ER_XDR_NSPACE;
734 if (nspace != ndim) return ER_NOT_HANDLED;
736 c = scan_header(file, "veclen", offset, 1);
737 if (c) veclen = atoi(c);
738 if (veclen<0 /*|| veclen>1000*/) return ER_XDR_VECLEN;
739 SetNumberOfComponents(veclen);
740 if (veclen==1) SetPixelType(itk::ImageIOBase::SCALAR);
741 else SetPixelType(itk::ImageIOBase::VECTOR);
743 c = scan_header(file, "data", offset, 1);
746 if (memicmp(c, "byte", 4) == 0 || memicmp(c, "xdr_byte", 8) == 0) SetComponentType(itk::ImageIOBase::CHAR);
747 else if (memicmp(c, "short", 5) == 0 || memicmp(c, "xdr_short", 9) == 0) SetComponentType(itk::ImageIOBase::SHORT);
748 else if (memicmp(c, "int" , 3) == 0 || memicmp(c, "xdr_int" , 7) == 0) SetComponentType(itk::ImageIOBase::INT);
749 else if (memicmp(c, "real", 4) == 0 || memicmp(c, "xdr_real", 8) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
750 else if (memicmp(c, "float", 5) == 0 || memicmp(c, "xdr_float", 9) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
751 else if (memicmp(c, "double",6) == 0 || memicmp(c, "xdr_double",10)== 0) SetComponentType(itk::ImageIOBase::DOUBLE);
752 else return ER_XDR_DATA;
754 if (memicmp(c, "xdr_", 4) == 0) forcenoswap=0;
758 c = scan_header(file, "field", offset, 1);
761 if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace *2;
762 else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
763 else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
764 else return ER_XDR_FIELD;
769 if (coords) /* expect AVS coordinates ? */
771 coords *= sizeof(float);
772 fstream = fopen(m_FileName.c_str(), "rb");
773 if (fstream == NULL) return ER_XDR_OPEN;
775 float *points = (float *)malloc(coords);
776 if (points == NULL) return ER_OUTOFMEMORY;
778 //Seek to coordinates position in file
779 if (fseek(fstream,-static_cast<int>(coords),SEEK_END)) return ER_XDR_READ;
780 if (fread( /*(*output)->*/points, 1, coords, fstream ) == coords)
781 { /* swap data if read-ok and required (xdr is low-endian) */
782 if (!(*(char *)(&swap_test)) && !forcenoswap)
784 c = (char *)/*(*output)->*/points;
785 for (i=0; i<coords; i+=4)
799 for (i=0; i<GetNumberOfDimensions(); i++) {
800 SetSpacing(i,10.*(points[i*2+1]-points[i*2])/(GetDimensions(i)-1));
801 SetOrigin(i,10.*points[i*2]);
805 //Rectilinear is reinterpreted as uniform because ITK does not know rectilinear
807 for (i=0; i<GetNumberOfDimensions(); i++) {
808 //Compute mean spacing
809 SetSpacing(i,10*(points[GetDimensions(i)-1]-points[0])/(GetDimensions(i)-1));
810 SetOrigin(i,10*points[0]);
812 //Test if rectilinear image is actually uniform (tolerance 0.1 mm)
813 for (j=0; j<GetDimensions(i)-1; j++) {
814 if (fabs((points[j+1]-points[j])*10-GetSpacing(i))>0.1) {
817 return ER_NOT_HANDLED;
820 points += (int)GetDimensions(i);
822 for (i=0; i<GetNumberOfDimensions(); i++)
823 points -= GetDimensions(i);
828 return ER_NOT_HANDLED;
836 //====================================================================
837 // Read image information (copied from XDRreader)
838 void clitk::XdrImageIO::ReadImageInformation() {
839 int result = ReadImageInformationWithError();
840 if (result) ITKError("clitk::XdrImageIO::ReadImageInformation",result);
844 //====================================================================
845 // Read Image Content (copied from Xdr reader)
846 int clitk::XdrImageIO::ReadWithError(void * buffer)
848 int /*ndim,*/ nspace/*, veclen=1, data=AVS_TYPE_BYTE, field=UNIFORM*/;
849 int iNkiCompression = 0;
850 int j, coords=0, datasize=0, HeaderSize;
851 unsigned int i,iNumRead,total=1;
856 //AVSfield FieldTemplate;
857 long swap_test = 0x1000000; /* first byte is 1 when low-endian */
859 char *file = const_cast<char *>(m_FileName.c_str());
861 AVSType field=UNIFORM;
863 for (i=0; i<GetNumberOfDimensions(); i++) coords += GetDimensions(i);
865 total = GetImageSizeInPixels();
866 nspace = GetNumberOfDimensions();
868 c = scan_header(file, "field", offset, 1);
871 if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace*2;
872 else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
873 else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
874 else return ER_XDR_FIELD;
879 c = scan_header(file, "nki_compression", offset, 1);
880 if (c) iNkiCompression = atoi(c);
882 c = scan_header(file, "coord1[0]", offset, 1);
883 if (c) HeaderSize = 32768;
884 else HeaderSize = 2048;
886 fstream = fopen(file, "rb");
890 if (offset) fseek(fstream, offset, SEEK_SET);
894 if (fgets(temp, 500, fstream) == NULL )
895 return ER_XDR_NOCTRLL; /* end of file */
897 if (temp[0] == 10) continue;
901 fseek(fstream, -2, SEEK_CUR);
903 } /* ^L end of header */
905 if (temp[0] != '#') break;
908 buff = (char*)malloc(HeaderSize);
911 return ER_OUTOFMEMORY;
913 memset(buff, 0, HeaderSize);
914 iNumRead = fread(buff, 1, HeaderSize, fstream);
922 for (i=0; i<iNumRead; i++) {
923 if (buff[i] == 12) break;
928 if (i==iNumRead) return ER_XDR_NOCTRLL;
930 total = GetImageSizeInBytes();
932 //We add casts because the resulting quantity can be negative.
933 //There is no risk of looping because i and iNumRead are about the size of the header
934 fseek(fstream, static_cast<int>(i)+2-static_cast<int>(iNumRead), SEEK_CUR);
936 if (total && iNkiCompression)
940 signed char* pCompressed;
942 /* Read or guess the size of the compressed data */
943 iCurPos = ftell(fstream);
944 iSize = get_nki_compressed_size(fstream);
948 fseek(fstream, 0, SEEK_END);
949 iSize = ftell(fstream);
950 iSize = iSize - iCurPos - coords;
952 // Get compressed size from header if possible; else use uncompressed size as safe estimate
953 if (iSize>total && offset) iSize=total+8;
956 fseek(fstream, iCurPos, SEEK_SET);
958 /* Allocate space for the compressed pixels */
959 pCompressed = (signed char*)malloc(iSize);
963 return ER_OUTOFMEMORY;
966 /* Read the compressed pixels */
967 if (fread( (void *)pCompressed, 1, iSize, fstream ) != iSize)
973 if (!nki_private_decompress((short*)buffer, pCompressed, iSize))
976 return ER_DECOMPRESSION;
980 fseek(fstream, iCurPos + global_len, SEEK_SET);
989 if (fread( (void *)buffer, 1, total, fstream ) != total)
996 /* swap data if required (xdr is low-endian) */
998 datasize = GetComponentSize();
999 if (!(*(char *)(&swap_test)) && !forcenoswap)
1004 for (i=0; i<total; i+=2)
1011 else if (datasize==4)
1014 for (i=0; i<total; i+=4)
1024 else if (datasize==8)
1027 for (i=0; i<total; i+=8)
1050 //====================================================================
1051 // Read image information (copied from Xdr reader)
1052 void clitk::XdrImageIO::Read(void * buffer) {
1053 int result = ReadWithError(buffer);
1054 if (result) ITKError("clitk::XdrImageIO::Read",result);
1057 //====================================================================
1058 // Read Image Information
1059 bool clitk::XdrImageIO::CanReadFile(const char* FileNameToRead)
1064 fstream = fopen(FileNameToRead, "rt");
1065 if (fstream == NULL)
1067 AVSerror("Couldn't open file " << FileNameToRead);
1070 fgets(temp, 500, fstream);
1073 if (memcmp(temp, "# AVS", 5)==0)
1079 void clitk::XdrImageIO::ITKError(std::string funcName, int msgID) {
1080 itkExceptionMacro(<< "Error in " << funcName << ". Message: " << gl_ErrorMsg[msgID]);