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
30 // Based on a true story by the Nederlands Kanker Instituut (AVS_RXDR.CPP from the 20080414)
37 /************************************************************************/
38 /* DEFINES, ENUMERATED TYPES AND CONSTANTS */
39 /************************************************************************/
60 unsigned int iOrgSize;
62 unsigned int iCompressedSize;
64 unsigned int iCompressedCRC; /* Excluding this header */
68 /************************************************************************/
69 /* GLOBAL VARIABLES */
70 /************************************************************************/
71 const char* gl_ErrorMsg[] = {
73 "Command or function not supported in this way.",
75 "XDR file header NDIM error",
76 "XDR file header DIMn error",
77 "XDR file header NSPACE error",
78 "XDR file header VECLEN error",
79 "XDR file header DATA(type) error",
80 "XDR file header FIELD(coordinate type) error",
81 "XDR file could not be opened",
82 "XDR file header contains no ^L",
83 "XDR file reading error",
85 "Decompression failed",
86 "Format not handled by clitkXdrImageIO (RECTILINEAR or IRREGULAR field)"
90 static const unsigned long CRC32_table[256] = {
91 0x00000000, 0x77073096, 0xee0e612c, 0x990951ba, 0x076dc419, 0x706af48f,
92 0xe963a535, 0x9e6495a3, 0x0edb8832, 0x79dcb8a4, 0xe0d5e91e, 0x97d2d988,
93 0x09b64c2b, 0x7eb17cbd, 0xe7b82d07, 0x90bf1d91, 0x1db71064, 0x6ab020f2,
94 0xf3b97148, 0x84be41de, 0x1adad47d, 0x6ddde4eb, 0xf4d4b551, 0x83d385c7,
95 0x136c9856, 0x646ba8c0, 0xfd62f97a, 0x8a65c9ec, 0x14015c4f, 0x63066cd9,
96 0xfa0f3d63, 0x8d080df5, 0x3b6e20c8, 0x4c69105e, 0xd56041e4, 0xa2677172,
97 0x3c03e4d1, 0x4b04d447, 0xd20d85fd, 0xa50ab56b, 0x35b5a8fa, 0x42b2986c,
98 0xdbbbc9d6, 0xacbcf940, 0x32d86ce3, 0x45df5c75, 0xdcd60dcf, 0xabd13d59,
99 0x26d930ac, 0x51de003a, 0xc8d75180, 0xbfd06116, 0x21b4f4b5, 0x56b3c423,
100 0xcfba9599, 0xb8bda50f, 0x2802b89e, 0x5f058808, 0xc60cd9b2, 0xb10be924,
101 0x2f6f7c87, 0x58684c11, 0xc1611dab, 0xb6662d3d, 0x76dc4190, 0x01db7106,
102 0x98d220bc, 0xefd5102a, 0x71b18589, 0x06b6b51f, 0x9fbfe4a5, 0xe8b8d433,
103 0x7807c9a2, 0x0f00f934, 0x9609a88e, 0xe10e9818, 0x7f6a0dbb, 0x086d3d2d,
104 0x91646c97, 0xe6635c01, 0x6b6b51f4, 0x1c6c6162, 0x856530d8, 0xf262004e,
105 0x6c0695ed, 0x1b01a57b, 0x8208f4c1, 0xf50fc457, 0x65b0d9c6, 0x12b7e950,
106 0x8bbeb8ea, 0xfcb9887c, 0x62dd1ddf, 0x15da2d49, 0x8cd37cf3, 0xfbd44c65,
107 0x4db26158, 0x3ab551ce, 0xa3bc0074, 0xd4bb30e2, 0x4adfa541, 0x3dd895d7,
108 0xa4d1c46d, 0xd3d6f4fb, 0x4369e96a, 0x346ed9fc, 0xad678846, 0xda60b8d0,
109 0x44042d73, 0x33031de5, 0xaa0a4c5f, 0xdd0d7cc9, 0x5005713c, 0x270241aa,
110 0xbe0b1010, 0xc90c2086, 0x5768b525, 0x206f85b3, 0xb966d409, 0xce61e49f,
111 0x5edef90e, 0x29d9c998, 0xb0d09822, 0xc7d7a8b4, 0x59b33d17, 0x2eb40d81,
112 0xb7bd5c3b, 0xc0ba6cad, 0xedb88320, 0x9abfb3b6, 0x03b6e20c, 0x74b1d29a,
113 0xead54739, 0x9dd277af, 0x04db2615, 0x73dc1683, 0xe3630b12, 0x94643b84,
114 0x0d6d6a3e, 0x7a6a5aa8, 0xe40ecf0b, 0x9309ff9d, 0x0a00ae27, 0x7d079eb1,
115 0xf00f9344, 0x8708a3d2, 0x1e01f268, 0x6906c2fe, 0xf762575d, 0x806567cb,
116 0x196c3671, 0x6e6b06e7, 0xfed41b76, 0x89d32be0, 0x10da7a5a, 0x67dd4acc,
117 0xf9b9df6f, 0x8ebeeff9, 0x17b7be43, 0x60b08ed5, 0xd6d6a3e8, 0xa1d1937e,
118 0x38d8c2c4, 0x4fdff252, 0xd1bb67f1, 0xa6bc5767, 0x3fb506dd, 0x48b2364b,
119 0xd80d2bda, 0xaf0a1b4c, 0x36034af6, 0x41047a60, 0xdf60efc3, 0xa867df55,
120 0x316e8eef, 0x4669be79, 0xcb61b38c, 0xbc66831a, 0x256fd2a0, 0x5268e236,
121 0xcc0c7795, 0xbb0b4703, 0x220216b9, 0x5505262f, 0xc5ba3bbe, 0xb2bd0b28,
122 0x2bb45a92, 0x5cb36a04, 0xc2d7ffa7, 0xb5d0cf31, 0x2cd99e8b, 0x5bdeae1d,
123 0x9b64c2b0, 0xec63f226, 0x756aa39c, 0x026d930a, 0x9c0906a9, 0xeb0e363f,
124 0x72076785, 0x05005713, 0x95bf4a82, 0xe2b87a14, 0x7bb12bae, 0x0cb61b38,
125 0x92d28e9b, 0xe5d5be0d, 0x7cdcefb7, 0x0bdbdf21, 0x86d3d2d4, 0xf1d4e242,
126 0x68ddb3f8, 0x1fda836e, 0x81be16cd, 0xf6b9265b, 0x6fb077e1, 0x18b74777,
127 0x88085ae6, 0xff0f6a70, 0x66063bca, 0x11010b5c, 0x8f659eff, 0xf862ae69,
128 0x616bffd3, 0x166ccf45, 0xa00ae278, 0xd70dd2ee, 0x4e048354, 0x3903b3c2,
129 0xa7672661, 0xd06016f7, 0x4969474d, 0x3e6e77db, 0xaed16a4a, 0xd9d65adc,
130 0x40df0b66, 0x37d83bf0, 0xa9bcae53, 0xdebb9ec5, 0x47b2cf7f, 0x30b5ffe9,
131 0xbdbdf21c, 0xcabac28a, 0x53b39330, 0x24b4a3a6, 0xbad03605, 0xcdd70693,
132 0x54de5729, 0x23d967bf, 0xb3667a2e, 0xc4614ab8, 0x5d681b02, 0x2a6f2b94,
133 0xb40bbe37, 0xc30c8ea1, 0x5a05df1b, 0x2d02ef8d
137 /************************************************************************/
138 /* MODULE FUNCTIONS */
139 /************************************************************************/
141 /* simple XDR file reader */
143 /* help routine, scans XDR file header for keyword entry, returns NULL
144 if not found, else returns pointer to rest of line
147 static char *scan_header(const char *file, const char *name, int offset, int removespaces)
149 int i, j, iStringLength;
150 static char temp[512];
154 if ((f = fopen(file, "rt")) == NULL) return NULL;
155 if (offset) fseek(f, offset, SEEK_SET);
159 if (fgets(temp, 500, f) == NULL ) break; /* end of file */
164 p = q = temp; /* remove spaces */
165 iStringLength = strlen(temp);
166 for (j=0; j<iStringLength; j++)
167 if (*q!=' ' && *q!=8) *p++ = *q++;
172 if (temp[0] == 12 ) break; /* ^L end of header */
173 if (temp[0] != '#') i++; /* The first 200 non comment lines must be read before data is opened. */
174 if ((p = strchr(temp+1, '=')) == NULL) continue; /* no '=' */
175 if (memicmp(temp, name, p-temp) ) continue; /* no match */
177 p++; /* match, skip = */
178 if (p[strlen(p)-1] == '\n') /* remove \n */
189 static int get_nki_compressed_size(FILE *f)
194 fread((void *)&Header, sizeof(Header), 1 , f);
196 iMode = Header.iMode;
205 return Header.iCompressedSize + sizeof(Header);
211 /* decoder for NKI private compressed pixel data
212 arguments: dest = (in) points to area where destination data is written (short)
213 src = (in) points compressed source data (byte stream)
214 size = (in) number of bytes source data (safety)
216 The return value is the number of pixels that have been processed.
218 The compressed data looks like:
219 (number of pixels)-1 times:
220 OR 1 byte = LL 7 bits signed (difference pixel[1] - pixel[0]);
221 OR 2 bytes = HHLL 15 bits signed (difference pixel[1] - pixel[0]) xored with 0x4000;
222 OR 3 bytes = 7FHHLL 16 bits absolute pixel data if 15 bits difference is exceeded
223 OR 2 bytes = 80NN run length encode NN zero differences (max 255)
224 for mode 3 and 4 added:
225 OR 2 bytes = CONN encode NN 4 bit differences (max 255)
227 Performance on typical CT or MRI is >2x compression and a very good speed
229 This code is not valid on HIGHENDIAN (high byte first) machines
232 static int global_len;
234 static int nki_private_decompress(short int *dest, signed char *src, int size)
236 int npixels, retvalue, mode, iMode, val, j;
237 NKI_MODE2* pHeader = (NKI_MODE2*)src;
238 unsigned long iCRC=0, iCRC2=0;
239 //unsigned char* pDestStart = (unsigned char*)dest;
240 signed char *save, *end;
242 retvalue = npixels = pHeader->iOrgSize;
243 iMode = pHeader->iMode; // safety: this value is checked in case statement
245 if (npixels<1) return 0; // safety: check for invalid npixels value
247 /* Up till now only Mode=1, 2, 3, and 4 are supported */
254 src += 8; // mode 1 only has 8 bytes header: iOrgSize and iMode
255 end = src + size - 3; // for overflow check if we are close to end of input buffer
257 *dest = *(short int *)src;
263 if (src > end) // check whether the last few messages fit in input buffer
265 if (src<end+3) val = *src;
268 if (val >= -64 && val <= 63) mode = 1; // 7 bit difference
269 else if (val==0x7f) mode = 3; // 16 bit value
270 else if ((val&0xff)==0x80) mode = 2; // run length encoding
273 if (src+mode > end+3)
274 return 0; // safety: overflow input data
279 if (val >= -64 && val <= 63) // 7 bit difference
281 dest[1] = dest[0] + val;
285 else if (val==0x7f) // 16 bit value
287 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
291 else if ((val&0xff)==0x80) // run length encoding
293 mode = ((unsigned char *)src)[1];
295 if (npixels<=0) return 0; // safety: overflow output data
306 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
307 dest[1] = dest[0] + diff; // 15 bit difference
314 global_len = src-save;
319 src += sizeof(NKI_MODE2);
321 end = src + pHeader->iCompressedSize - 3;
323 if (end > src + size - 3)
324 end = src + size - 3; // may occur if pHeader is corrupted
326 *dest = val = *(short int *)src;
327 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
328 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
335 if (src > end) // check whether the last few messages fit in input buffer
337 if (src<end+3) val = *src;
340 if (val >= -64 && val <= 63) mode = 1; // 7 bit difference
341 else if (val==0x7f) mode = 3; // 16 bit value
342 else if ((val&0xff)==0x80) mode = 2; // run length encoding
345 if (src+mode > end+3)
346 break; // safety: overflow input data
351 if (val >= -64 && val <= 63) // 7 bits difference
353 dest[1] = val = dest[0] + val;
354 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
355 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
359 else if (val==0x7f) // 16 bit value
361 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
363 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
364 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
368 else if ((val&0xff)==0x80) // run length encoding
370 mode = ((unsigned char *)src)[1];
372 if (npixels<=0) break; // safety: overflow output data
375 dest[1] = val = dest[0];
376 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
377 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
385 signed short diff = ((val^0x40)<<8) + ((unsigned char *)src)[1];
386 dest[1] = val = dest[0] + diff; // 15 bit difference
387 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
388 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
395 if (iCRC2 != pHeader->iOrgCRC) // if error in output CRC:
397 src = save; // check input CRC
400 iCRC = CRC32_table[(unsigned char)iCRC ^ (unsigned char)src[0]] ^ ((iCRC >> 8));
404 if (iCRC != pHeader->iCompressedCRC)
406 AVSerror("XDR decompression: the file is corrupted");
411 AVSerror("XDR decompression: internal error");
416 global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
423 src += 8; // mode 3 only has 8 bytes header: iOrgSize and iMode
424 end = src + size - 3; // for overflow check if we are close to end of input buffer
426 *dest = *(short int *)src;
432 if (src > end) // check whether the last few messages fit in input buffer
434 if (src<end+3) val = *src;
437 if (val >= -63 && val <= 63) mode = 1; // 7 bit difference
438 else if (val==0x7f) mode = 3; // 16 bit value
439 else if ((val&0xff)==0x80) mode = 2; // run length encoding
440 else if ((val&0xff)==0xC0) mode = 2; // 4 bit encoding
443 if (src+mode > end+3)
444 return 0; // safety: overflow input data
449 if (val >= -63 && val <= 63) // 7 bit difference
451 dest[1] = dest[0] + val;
455 else if (val==0x7f) // 16 bit value
457 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
461 else if ((val&0xff)==0x80) // run length encoding
463 mode = ((unsigned char *)src)[1];
465 if (npixels<=0) return 0; // safety: overflow output data
474 else if ((val&0xff)==0xC0) // 4 bit run
476 mode = ((unsigned char *)src)[1];
480 if (npixels<=0) return 0; // safety: overflow output data
484 dest[1] = dest[0] + (val>>4);
486 if (val&8) val |= 0xfffffff0;
488 dest[1] = dest[0] + val;
495 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
496 dest[1] = dest[0] + diff; // 15 bit difference
503 global_len = src-save;
508 src += sizeof(NKI_MODE2);
510 end = src + pHeader->iCompressedSize - 3;
512 if (end > src + size - 3)
513 end = src + size - 3; // may occur if pHeader is corrupted
515 *dest = val = *(short int *)src;
516 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
517 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
523 if (src > end) // check whether the last few messages fit in input buffer
525 if (src<end+3) val = *src;
528 if (val >= -63 && val <= 63) mode = 1; // 7 bit difference
529 else if (val==0x7f) mode = 3; // 16 bit value
530 else if ((val&0xff)==0x80) mode = 2; // run length encoding
531 else if ((val&0xff)==0xC0) mode = 2; // 4 bit encoding
534 if (src+mode > end+3)
535 return 0; // safety: overflow input data
540 if (val >= -63 && val <= 63) // 7 bit difference
542 dest[1] = val = dest[0] + val;
543 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
544 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
548 else if (val==0x7f) // 16 bit value
550 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
551 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
552 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
556 else if ((val&0xff)==0x80) // run length encoding
558 mode = ((unsigned char *)src)[1];
560 if (npixels<=0) return 0; // safety: overflow output data
563 dest[1] = val = dest[0];
564 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
565 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
571 else if ((val&0xff)==0xC0) // 4 bit run
573 mode = ((unsigned char *)src)[1];
577 if (npixels<=0) return 0; // safety: overflow output data
581 dest[1] = j = dest[0] + (val>>4);
582 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
583 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
585 if (val&8) val |= 0xfffffff0;
587 dest[1] = j = dest[0] + val;
588 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
589 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
596 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
597 dest[1] = val = dest[0] + diff; // 15 bit difference
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));
606 if (iCRC2 != pHeader->iOrgCRC) // if error in output CRC:
609 global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
615 AVSerror("XDR decompression: unsupported mode");
623 //====================================================================
624 // Read image information (copied from XDRreader)
625 int clitk::XdrImageIO::ReadImageInformationWithError()
628 itk::Vector<int,MAXDIM> dim;
631 unsigned int coords=0,i,j,ndim,nspace;
636 long swap_test = 0x1000000; /* first byte is 1 when low-endian */
638 char *file = const_cast<char *>(m_FileName.c_str());
639 AVSType field=UNIFORM;
642 fstream = fopen(file, "rt");
643 if (fstream == NULL) return ER_XDR_OPEN;
645 fgets(temp, 500, fstream);
648 if (memcmp(temp, "# AVS field file (produced by avs_nfwrite.c)", 44)==0) forcenoswap=1;
650 c = scan_header(file, "ndim", offset, 1);
651 if (!c) return ER_XDR_NDIM;
654 if (ndim<1 || ndim>MAXDIM) return ER_XDR_NDIM;
655 SetNumberOfDimensions(ndim);
659 for (i=0; i<ndim; i++)
661 sprintf(temp, "dim%d", i+1);
662 c = scan_header(file, temp, offset, 1);
663 if (!c) return ER_XDR_DIM;
665 if (dim[i]<1) return ER_XDR_DIM;
670 for (i=0; i<ndim; i++) {
671 SetDimensions(i,dim[i]);
676 c = scan_header(file, "nspace", offset, 1);
677 if (c) nspace = atoi(c);
678 if (nspace<1 || ndim > MAXDIM) return ER_XDR_NSPACE;
679 if (nspace != ndim) return ER_NOT_HANDLED;
681 c = scan_header(file, "veclen", offset, 1);
682 if (c) veclen = atoi(c);
683 if (veclen<0 /*|| veclen>1000*/) return ER_XDR_VECLEN;
684 SetNumberOfComponents(veclen);
685 if (veclen==1) SetPixelType(itk::ImageIOBase::SCALAR);
686 else SetPixelType(itk::ImageIOBase::VECTOR);
688 c = scan_header(file, "data", offset, 1);
691 if (memicmp(c, "byte", 4) == 0 || memicmp(c, "xdr_byte", 8) == 0) SetComponentType(itk::ImageIOBase::CHAR);
692 else if (memicmp(c, "short", 5) == 0 || memicmp(c, "xdr_short", 9) == 0) SetComponentType(itk::ImageIOBase::SHORT);
693 else if (memicmp(c, "int" , 3) == 0 || memicmp(c, "xdr_int" , 7) == 0) SetComponentType(itk::ImageIOBase::INT);
694 else if (memicmp(c, "real", 4) == 0 || memicmp(c, "xdr_real", 8) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
695 else if (memicmp(c, "float", 5) == 0 || memicmp(c, "xdr_float", 9) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
696 else if (memicmp(c, "double",6) == 0 || memicmp(c, "xdr_double",10)== 0) SetComponentType(itk::ImageIOBase::DOUBLE);
697 else return ER_XDR_DATA;
699 if (memicmp(c, "xdr_", 4) == 0) forcenoswap=0;
703 c = scan_header(file, "field", offset, 1);
706 if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace *2;
707 else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
708 else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
709 else return ER_XDR_FIELD;
714 if (coords) /* expect AVS coordinates ? */
716 coords *= sizeof(float);
717 fstream = fopen(m_FileName.c_str(), "rb");
718 if (fstream == NULL) return ER_XDR_OPEN;
720 float *points = (float *)malloc(coords);
721 if (points == NULL) return ER_OUTOFMEMORY;
723 //Seek to coordinates position in file
724 if (fseek(fstream,-static_cast<int>(coords),SEEK_END)) return ER_XDR_READ;
725 if (fread( /*(*output)->*/points, 1, coords, fstream ) == coords)
726 { /* swap data if read-ok and required (xdr is low-endian) */
727 if (!(*(char *)(&swap_test)) && !forcenoswap)
729 c = (char *)/*(*output)->*/points;
730 for (i=0; i<coords; i+=4)
744 for (i=0; i<GetNumberOfDimensions(); i++) {
745 SetSpacing(i,10.*(points[i*2+1]-points[i*2])/(GetDimensions(i)-1));
746 SetOrigin(i,10.*points[i*2]);
750 //Rectilinear is reinterpreted as uniform because ITK does not know rectilinear
752 for (i=0; i<GetNumberOfDimensions(); i++) {
753 //Compute mean spacing
754 SetSpacing(i,10*(points[GetDimensions(i)-1]-points[0])/(GetDimensions(i)-1));
755 SetOrigin(i,10*points[0]);
757 //Test if rectilinear image is actually uniform (tolerance 0.1 mm)
758 for (j=0; j<GetDimensions(i)-1; j++) {
759 if (fabs((points[j+1]-points[j])*10-GetSpacing(i))>0.1) {
762 return ER_NOT_HANDLED;
765 points += (int)GetDimensions(i);
767 for (i=0; i<GetNumberOfDimensions(); i++)
768 points -= GetDimensions(i);
773 return ER_NOT_HANDLED;
781 //====================================================================
782 // Read image information (copied from XDRreader)
783 void clitk::XdrImageIO::ReadImageInformation() {
784 int result = ReadImageInformationWithError();
785 if (result) ITKError("clitk::XdrImageIO::ReadImageInformation",result);
789 //====================================================================
790 // Read Image Content (copied from Xdr reader)
791 int clitk::XdrImageIO::ReadWithError(void * buffer)
793 int /*ndim,*/ nspace/*, veclen=1, data=AVS_TYPE_BYTE, field=UNIFORM*/;
794 int iNkiCompression = 0;
795 int j, coords=0, datasize=0, HeaderSize;
796 unsigned int i,iNumRead,total=1;
801 //AVSfield FieldTemplate;
802 long swap_test = 0x1000000; /* first byte is 1 when low-endian */
804 char *file = const_cast<char *>(m_FileName.c_str());
806 AVSType field=UNIFORM;
808 for (i=0; i<GetNumberOfDimensions(); i++) coords += GetDimensions(i);
810 total = GetImageSizeInPixels();
811 nspace = GetNumberOfDimensions();
813 c = scan_header(file, "field", offset, 1);
816 if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace*2;
817 else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
818 else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
819 else return ER_XDR_FIELD;
824 c = scan_header(file, "nki_compression", offset, 1);
825 if (c) iNkiCompression = atoi(c);
827 c = scan_header(file, "coord1[0]", offset, 1);
828 if (c) HeaderSize = 32768;
829 else HeaderSize = 2048;
831 fstream = fopen(file, "rb");
835 if (offset) fseek(fstream, offset, SEEK_SET);
839 if (fgets(temp, 500, fstream) == NULL )
840 return ER_XDR_NOCTRLL; /* end of file */
842 if (temp[0] == 10) continue;
846 fseek(fstream, -2, SEEK_CUR);
848 } /* ^L end of header */
850 if (temp[0] != '#') break;
853 buff = (char*)malloc(HeaderSize);
856 return ER_OUTOFMEMORY;
858 memset(buff, 0, HeaderSize);
859 iNumRead = fread(buff, 1, HeaderSize, fstream);
867 for (i=0; i<iNumRead; i++) {
868 if (buff[i] == 12) break;
873 if (i==iNumRead) return ER_XDR_NOCTRLL;
875 total = GetImageSizeInBytes();
877 //We add casts because the resulting quantity can be negative.
878 //There is no risk of looping because i and iNumRead are about the size of the header
879 fseek(fstream, static_cast<int>(i)+2-static_cast<int>(iNumRead), SEEK_CUR);
881 if (total && iNkiCompression)
885 signed char* pCompressed;
887 /* Read or guess the size of the compressed data */
888 iCurPos = ftell(fstream);
889 iSize = get_nki_compressed_size(fstream);
893 fseek(fstream, 0, SEEK_END);
894 iSize = ftell(fstream);
895 iSize = iSize - iCurPos - coords;
897 // Get compressed size from header if possible; else use uncompressed size as safe estimate
898 if (iSize>total && offset) iSize=total+8;
901 fseek(fstream, iCurPos, SEEK_SET);
903 /* Allocate space for the compressed pixels */
904 pCompressed = (signed char*)malloc(iSize);
908 return ER_OUTOFMEMORY;
911 /* Read the compressed pixels */
912 if (fread( (void *)pCompressed, 1, iSize, fstream ) != iSize)
918 if (!nki_private_decompress((short*)buffer, pCompressed, iSize))
921 return ER_DECOMPRESSION;
925 fseek(fstream, iCurPos + global_len, SEEK_SET);
934 if (fread( (void *)buffer, 1, total, fstream ) != total)
941 /* swap data if required (xdr is low-endian) */
943 datasize = GetComponentSize();
944 if (!(*(char *)(&swap_test)) && !forcenoswap)
949 for (i=0; i<total; i+=2)
956 else if (datasize==4)
959 for (i=0; i<total; i+=4)
969 else if (datasize==8)
972 for (i=0; i<total; i+=8)
995 //====================================================================
996 // Read image information (copied from Xdr reader)
997 void clitk::XdrImageIO::Read(void * buffer) {
998 int result = ReadWithError(buffer);
999 if (result) ITKError("clitk::XdrImageIO::Read",result);
1002 //====================================================================
1003 // Read Image Information
1004 bool clitk::XdrImageIO::CanReadFile(const char* FileNameToRead)
1009 fstream = fopen(FileNameToRead, "rt");
1010 if (fstream == NULL)
1012 AVSerror("Couldn't open file " << FileNameToRead);
1015 fgets(temp, 500, fstream);
1018 if (memcmp(temp, "# AVS", 5)==0)
1024 void clitk::XdrImageIO::ITKError(std::string funcName, int msgID) {
1025 itkExceptionMacro(<< "Error in " << funcName << ". Message: " << gl_ErrorMsg[msgID]);