1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
19 * @file clitkXdrImageIO.cxx
20 * @author Simon Rit <simon.rit@gmail.com>
21 * @date Sun Jun 1 22:12:20 2008
28 #include "clitkXdrImageIO.h"
41 #define memicmp _memicmp
43 #define memicmp strncasecmp
47 // Based on a true story by the Nederlands Kanker Instituut (AVS_RXDR.CPP from the 20080414)
54 /************************************************************************/
55 /* DEFINES, ENUMERATED TYPES AND CONSTANTS */
56 /************************************************************************/
77 unsigned int iOrgSize;
79 unsigned int iCompressedSize;
81 unsigned int iCompressedCRC; /* Excluding this header */
85 /************************************************************************/
86 /* GLOBAL VARIABLES */
87 /************************************************************************/
88 const char* gl_ErrorMsg[] = {
90 "Command or function not supported in this way.",
92 "XDR file header NDIM error",
93 "XDR file header DIMn error",
94 "XDR file header NSPACE error",
95 "XDR file header VECLEN error",
96 "XDR file header DATA(type) error",
97 "XDR file header FIELD(coordinate type) error",
98 "XDR file could not be opened",
99 "XDR file header contains no ^L",
100 "XDR file reading error",
102 "Decompression failed",
103 "Format not handled by clitkXdrImageIO (RECTILINEAR or IRREGULAR field)"
107 static const unsigned long CRC32_table[256] = {
108 0x00000000, 0x77073096, 0xee0e612c, 0x990951ba, 0x076dc419, 0x706af48f,
109 0xe963a535, 0x9e6495a3, 0x0edb8832, 0x79dcb8a4, 0xe0d5e91e, 0x97d2d988,
110 0x09b64c2b, 0x7eb17cbd, 0xe7b82d07, 0x90bf1d91, 0x1db71064, 0x6ab020f2,
111 0xf3b97148, 0x84be41de, 0x1adad47d, 0x6ddde4eb, 0xf4d4b551, 0x83d385c7,
112 0x136c9856, 0x646ba8c0, 0xfd62f97a, 0x8a65c9ec, 0x14015c4f, 0x63066cd9,
113 0xfa0f3d63, 0x8d080df5, 0x3b6e20c8, 0x4c69105e, 0xd56041e4, 0xa2677172,
114 0x3c03e4d1, 0x4b04d447, 0xd20d85fd, 0xa50ab56b, 0x35b5a8fa, 0x42b2986c,
115 0xdbbbc9d6, 0xacbcf940, 0x32d86ce3, 0x45df5c75, 0xdcd60dcf, 0xabd13d59,
116 0x26d930ac, 0x51de003a, 0xc8d75180, 0xbfd06116, 0x21b4f4b5, 0x56b3c423,
117 0xcfba9599, 0xb8bda50f, 0x2802b89e, 0x5f058808, 0xc60cd9b2, 0xb10be924,
118 0x2f6f7c87, 0x58684c11, 0xc1611dab, 0xb6662d3d, 0x76dc4190, 0x01db7106,
119 0x98d220bc, 0xefd5102a, 0x71b18589, 0x06b6b51f, 0x9fbfe4a5, 0xe8b8d433,
120 0x7807c9a2, 0x0f00f934, 0x9609a88e, 0xe10e9818, 0x7f6a0dbb, 0x086d3d2d,
121 0x91646c97, 0xe6635c01, 0x6b6b51f4, 0x1c6c6162, 0x856530d8, 0xf262004e,
122 0x6c0695ed, 0x1b01a57b, 0x8208f4c1, 0xf50fc457, 0x65b0d9c6, 0x12b7e950,
123 0x8bbeb8ea, 0xfcb9887c, 0x62dd1ddf, 0x15da2d49, 0x8cd37cf3, 0xfbd44c65,
124 0x4db26158, 0x3ab551ce, 0xa3bc0074, 0xd4bb30e2, 0x4adfa541, 0x3dd895d7,
125 0xa4d1c46d, 0xd3d6f4fb, 0x4369e96a, 0x346ed9fc, 0xad678846, 0xda60b8d0,
126 0x44042d73, 0x33031de5, 0xaa0a4c5f, 0xdd0d7cc9, 0x5005713c, 0x270241aa,
127 0xbe0b1010, 0xc90c2086, 0x5768b525, 0x206f85b3, 0xb966d409, 0xce61e49f,
128 0x5edef90e, 0x29d9c998, 0xb0d09822, 0xc7d7a8b4, 0x59b33d17, 0x2eb40d81,
129 0xb7bd5c3b, 0xc0ba6cad, 0xedb88320, 0x9abfb3b6, 0x03b6e20c, 0x74b1d29a,
130 0xead54739, 0x9dd277af, 0x04db2615, 0x73dc1683, 0xe3630b12, 0x94643b84,
131 0x0d6d6a3e, 0x7a6a5aa8, 0xe40ecf0b, 0x9309ff9d, 0x0a00ae27, 0x7d079eb1,
132 0xf00f9344, 0x8708a3d2, 0x1e01f268, 0x6906c2fe, 0xf762575d, 0x806567cb,
133 0x196c3671, 0x6e6b06e7, 0xfed41b76, 0x89d32be0, 0x10da7a5a, 0x67dd4acc,
134 0xf9b9df6f, 0x8ebeeff9, 0x17b7be43, 0x60b08ed5, 0xd6d6a3e8, 0xa1d1937e,
135 0x38d8c2c4, 0x4fdff252, 0xd1bb67f1, 0xa6bc5767, 0x3fb506dd, 0x48b2364b,
136 0xd80d2bda, 0xaf0a1b4c, 0x36034af6, 0x41047a60, 0xdf60efc3, 0xa867df55,
137 0x316e8eef, 0x4669be79, 0xcb61b38c, 0xbc66831a, 0x256fd2a0, 0x5268e236,
138 0xcc0c7795, 0xbb0b4703, 0x220216b9, 0x5505262f, 0xc5ba3bbe, 0xb2bd0b28,
139 0x2bb45a92, 0x5cb36a04, 0xc2d7ffa7, 0xb5d0cf31, 0x2cd99e8b, 0x5bdeae1d,
140 0x9b64c2b0, 0xec63f226, 0x756aa39c, 0x026d930a, 0x9c0906a9, 0xeb0e363f,
141 0x72076785, 0x05005713, 0x95bf4a82, 0xe2b87a14, 0x7bb12bae, 0x0cb61b38,
142 0x92d28e9b, 0xe5d5be0d, 0x7cdcefb7, 0x0bdbdf21, 0x86d3d2d4, 0xf1d4e242,
143 0x68ddb3f8, 0x1fda836e, 0x81be16cd, 0xf6b9265b, 0x6fb077e1, 0x18b74777,
144 0x88085ae6, 0xff0f6a70, 0x66063bca, 0x11010b5c, 0x8f659eff, 0xf862ae69,
145 0x616bffd3, 0x166ccf45, 0xa00ae278, 0xd70dd2ee, 0x4e048354, 0x3903b3c2,
146 0xa7672661, 0xd06016f7, 0x4969474d, 0x3e6e77db, 0xaed16a4a, 0xd9d65adc,
147 0x40df0b66, 0x37d83bf0, 0xa9bcae53, 0xdebb9ec5, 0x47b2cf7f, 0x30b5ffe9,
148 0xbdbdf21c, 0xcabac28a, 0x53b39330, 0x24b4a3a6, 0xbad03605, 0xcdd70693,
149 0x54de5729, 0x23d967bf, 0xb3667a2e, 0xc4614ab8, 0x5d681b02, 0x2a6f2b94,
150 0xb40bbe37, 0xc30c8ea1, 0x5a05df1b, 0x2d02ef8d
154 /************************************************************************/
155 /* MODULE FUNCTIONS */
156 /************************************************************************/
158 /* simple XDR file reader */
160 /* help routine, scans XDR file header for keyword entry, returns NULL
161 if not found, else returns pointer to rest of line
164 static char *scan_header(const char *file, const char *name, int offset, int removespaces)
166 int i, j, iStringLength;
167 static char temp[512];
171 if ((f = fopen(file, "rt")) == NULL) return NULL;
172 if (offset) fseek(f, offset, SEEK_SET);
176 if (fgets(temp, 500, f) == NULL ) break; /* end of file */
181 p = q = temp; /* remove spaces */
182 iStringLength = strlen(temp);
183 for (j=0; j<iStringLength; j++)
184 if (*q!=' ' && *q!=8) *p++ = *q++;
189 if (temp[0] == 12 ) break; /* ^L end of header */
190 if (temp[0] != '#') i++; /* The first 200 non comment lines must be read before data is opened. */
191 if ((p = strchr(temp+1, '=')) == NULL) continue; /* no '=' */
192 if (memicmp(temp, name, p-temp) ) continue; /* no match */
194 p++; /* match, skip = */
195 if (p[strlen(p)-1] == '\n') /* remove \n */
206 static int get_nki_compressed_size(FILE *f)
211 fread((void *)&Header, sizeof(Header), 1 , f);
213 iMode = Header.iMode;
222 return Header.iCompressedSize + sizeof(Header);
228 /* decoder for NKI private compressed pixel data
229 arguments: dest = (in) points to area where destination data is written (short)
230 src = (in) points compressed source data (byte stream)
231 size = (in) number of bytes source data (safety)
233 The return value is the number of pixels that have been processed.
235 The compressed data looks like:
236 (number of pixels)-1 times:
237 OR 1 byte = LL 7 bits signed (difference pixel[1] - pixel[0]);
238 OR 2 bytes = HHLL 15 bits signed (difference pixel[1] - pixel[0]) xored with 0x4000;
239 OR 3 bytes = 7FHHLL 16 bits absolute pixel data if 15 bits difference is exceeded
240 OR 2 bytes = 80NN run length encode NN zero differences (max 255)
241 for mode 3 and 4 added:
242 OR 2 bytes = CONN encode NN 4 bit differences (max 255)
244 Performance on typical CT or MRI is >2x compression and a very good speed
246 This code is not valid on HIGHENDIAN (high byte first) machines
249 static int global_len;
251 static int nki_private_decompress(short int *dest, signed char *src, int size)
253 int npixels, retvalue, mode, iMode, val, j;
254 NKI_MODE2* pHeader = (NKI_MODE2*)src;
255 unsigned long iCRC=0, iCRC2=0;
256 //unsigned char* pDestStart = (unsigned char*)dest;
257 signed char *save, *end;
259 retvalue = npixels = pHeader->iOrgSize;
260 iMode = pHeader->iMode; // safety: this value is checked in case statement
262 if (npixels<1) return 0; // safety: check for invalid npixels value
264 /* Up till now only Mode=1, 2, 3, and 4 are supported */
271 src += 8; // mode 1 only has 8 bytes header: iOrgSize and iMode
272 end = src + size - 3; // for overflow check if we are close to end of input buffer
274 *dest = *(short int *)src;
280 if (src > end) // check whether the last few messages fit in input buffer
282 if (src<end+3) val = *src;
285 if (val >= -64 && val <= 63) mode = 1; // 7 bit difference
286 else if (val==0x7f) mode = 3; // 16 bit value
287 else if ((val&0xff)==0x80) mode = 2; // run length encoding
290 if (src+mode > end+3)
291 return 0; // safety: overflow input data
296 if (val >= -64 && val <= 63) // 7 bit difference
298 dest[1] = dest[0] + val;
302 else if (val==0x7f) // 16 bit value
304 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
308 else if ((val&0xff)==0x80) // run length encoding
310 mode = ((unsigned char *)src)[1];
312 if (npixels<=0) return 0; // safety: overflow output data
323 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
324 dest[1] = dest[0] + diff; // 15 bit difference
331 global_len = src-save;
336 src += sizeof(NKI_MODE2);
338 end = src + pHeader->iCompressedSize - 3;
340 if (end > src + size - 3)
341 end = src + size - 3; // may occur if pHeader is corrupted
343 *dest = val = *(short int *)src;
344 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
345 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
352 if (src > end) // check whether the last few messages fit in input buffer
354 if (src<end+3) val = *src;
357 if (val >= -64 && val <= 63) mode = 1; // 7 bit difference
358 else if (val==0x7f) mode = 3; // 16 bit value
359 else if ((val&0xff)==0x80) mode = 2; // run length encoding
362 if (src+mode > end+3)
363 break; // safety: overflow input data
368 if (val >= -64 && val <= 63) // 7 bits difference
370 dest[1] = val = dest[0] + val;
371 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
372 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
376 else if (val==0x7f) // 16 bit value
378 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
380 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
381 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
385 else if ((val&0xff)==0x80) // run length encoding
387 mode = ((unsigned char *)src)[1];
389 if (npixels<=0) break; // safety: overflow output data
392 dest[1] = val = dest[0];
393 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
394 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
402 signed short diff = ((val^0x40)<<8) + ((unsigned char *)src)[1];
403 dest[1] = val = dest[0] + diff; // 15 bit difference
404 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
405 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
412 if (iCRC2 != pHeader->iOrgCRC) // if error in output CRC:
414 src = save; // check input CRC
417 iCRC = CRC32_table[(unsigned char)iCRC ^ (unsigned char)src[0]] ^ ((iCRC >> 8));
421 if (iCRC != pHeader->iCompressedCRC)
423 AVSerror("XDR decompression: the file is corrupted");
428 AVSerror("XDR decompression: internal error");
433 global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
440 src += 8; // mode 3 only has 8 bytes header: iOrgSize and iMode
441 end = src + size - 3; // for overflow check if we are close to end of input buffer
443 *dest = *(short int *)src;
449 if (src > end) // check whether the last few messages fit in input buffer
451 if (src<end+3) val = *src;
454 if (val >= -63 && val <= 63) mode = 1; // 7 bit difference
455 else if (val==0x7f) mode = 3; // 16 bit value
456 else if ((val&0xff)==0x80) mode = 2; // run length encoding
457 else if ((val&0xff)==0xC0) mode = 2; // 4 bit encoding
460 if (src+mode > end+3)
461 return 0; // safety: overflow input data
466 if (val >= -63 && val <= 63) // 7 bit difference
468 dest[1] = dest[0] + val;
472 else if (val==0x7f) // 16 bit value
474 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
478 else if ((val&0xff)==0x80) // run length encoding
480 mode = ((unsigned char *)src)[1];
482 if (npixels<=0) return 0; // safety: overflow output data
491 else if ((val&0xff)==0xC0) // 4 bit run
493 mode = ((unsigned char *)src)[1];
497 if (npixels<=0) return 0; // safety: overflow output data
501 dest[1] = dest[0] + (val>>4);
503 if (val&8) val |= 0xfffffff0;
505 dest[1] = dest[0] + val;
512 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
513 dest[1] = dest[0] + diff; // 15 bit difference
520 global_len = src-save;
525 src += sizeof(NKI_MODE2);
527 end = src + pHeader->iCompressedSize - 3;
529 if (end > src + size - 3)
530 end = src + size - 3; // may occur if pHeader is corrupted
532 *dest = val = *(short int *)src;
533 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
534 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
540 if (src > end) // check whether the last few messages fit in input buffer
542 if (src<end+3) val = *src;
545 if (val >= -63 && val <= 63) mode = 1; // 7 bit difference
546 else if (val==0x7f) mode = 3; // 16 bit value
547 else if ((val&0xff)==0x80) mode = 2; // run length encoding
548 else if ((val&0xff)==0xC0) mode = 2; // 4 bit encoding
551 if (src+mode > end+3)
552 return 0; // safety: overflow input data
557 if (val >= -63 && val <= 63) // 7 bit difference
559 dest[1] = val = dest[0] + val;
560 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
561 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
565 else if (val==0x7f) // 16 bit value
567 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
568 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
569 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
573 else if ((val&0xff)==0x80) // run length encoding
575 mode = ((unsigned char *)src)[1];
577 if (npixels<=0) return 0; // safety: overflow output data
580 dest[1] = val = dest[0];
581 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
582 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
588 else if ((val&0xff)==0xC0) // 4 bit run
590 mode = ((unsigned char *)src)[1];
594 if (npixels<=0) return 0; // safety: overflow output data
598 dest[1] = j = dest[0] + (val>>4);
599 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
600 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
602 if (val&8) val |= 0xfffffff0;
604 dest[1] = j = dest[0] + val;
605 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
606 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
613 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
614 dest[1] = val = dest[0] + diff; // 15 bit difference
615 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
616 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
623 if (iCRC2 != pHeader->iOrgCRC) // if error in output CRC:
626 global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
632 AVSerror("XDR decompression: unsupported mode");
640 //====================================================================
641 // Read image information (copied from XDRreader)
642 int clitk::XdrImageIO::ReadImageInformationWithError()
645 itk::Vector<int,MAXDIM> dim;
648 unsigned int coords=0,i,j,ndim,nspace;
653 long swap_test = 0x1000000; /* first byte is 1 when low-endian */
655 char *file = const_cast<char *>(m_FileName.c_str());
656 AVSType field=UNIFORM;
659 fstream = fopen(file, "rt");
660 if (fstream == NULL) return ER_XDR_OPEN;
662 fgets(temp, 500, fstream);
665 if (memcmp(temp, "# AVS field file (produced by avs_nfwrite.c)", 44)==0) forcenoswap=1;
667 c = scan_header(file, "ndim", offset, 1);
668 if (!c) return ER_XDR_NDIM;
671 if (ndim<1 || ndim>MAXDIM) return ER_XDR_NDIM;
672 SetNumberOfDimensions(ndim);
676 for (i=0; i<ndim; i++)
678 sprintf(temp, "dim%d", i+1);
679 c = scan_header(file, temp, offset, 1);
680 if (!c) return ER_XDR_DIM;
682 if (dim[i]<1) return ER_XDR_DIM;
687 for (i=0; i<ndim; i++) {
688 SetDimensions(i,dim[i]);
693 c = scan_header(file, "nspace", offset, 1);
694 if (c) nspace = atoi(c);
695 if (nspace<1 || ndim > MAXDIM) return ER_XDR_NSPACE;
696 if (nspace != ndim) return ER_NOT_HANDLED;
698 c = scan_header(file, "veclen", offset, 1);
699 if (c) veclen = atoi(c);
700 if (veclen<0 /*|| veclen>1000*/) return ER_XDR_VECLEN;
701 SetNumberOfComponents(veclen);
702 if (veclen==1) SetPixelType(itk::ImageIOBase::SCALAR);
703 else SetPixelType(itk::ImageIOBase::VECTOR);
705 c = scan_header(file, "data", offset, 1);
708 if (memicmp(c, "byte", 4) == 0 || memicmp(c, "xdr_byte", 8) == 0) SetComponentType(itk::ImageIOBase::CHAR);
709 else if (memicmp(c, "short", 5) == 0 || memicmp(c, "xdr_short", 9) == 0) SetComponentType(itk::ImageIOBase::SHORT);
710 else if (memicmp(c, "int" , 3) == 0 || memicmp(c, "xdr_int" , 7) == 0) SetComponentType(itk::ImageIOBase::INT);
711 else if (memicmp(c, "real", 4) == 0 || memicmp(c, "xdr_real", 8) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
712 else if (memicmp(c, "float", 5) == 0 || memicmp(c, "xdr_float", 9) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
713 else if (memicmp(c, "double",6) == 0 || memicmp(c, "xdr_double",10)== 0) SetComponentType(itk::ImageIOBase::DOUBLE);
714 else return ER_XDR_DATA;
716 if (memicmp(c, "xdr_", 4) == 0) forcenoswap=0;
720 c = scan_header(file, "field", offset, 1);
723 if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace *2;
724 else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
725 else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
726 else return ER_XDR_FIELD;
731 if (coords) /* expect AVS coordinates ? */
733 coords *= sizeof(float);
734 fstream = fopen(m_FileName.c_str(), "rb");
735 if (fstream == NULL) return ER_XDR_OPEN;
737 float *points = (float *)malloc(coords);
738 if (points == NULL) return ER_OUTOFMEMORY;
740 //Seek to coordinates position in file
741 if (fseek(fstream,-static_cast<int>(coords),SEEK_END)) return ER_XDR_READ;
742 if (fread( /*(*output)->*/points, 1, coords, fstream ) == coords)
743 { /* swap data if read-ok and required (xdr is low-endian) */
744 if (!(*(char *)(&swap_test)) && !forcenoswap)
746 c = (char *)/*(*output)->*/points;
747 for (i=0; i<coords; i+=4)
761 for (i=0; i<GetNumberOfDimensions(); i++) {
762 SetSpacing(i,10.*(points[i*2+1]-points[i*2])/(GetDimensions(i)-1));
763 SetOrigin(i,10.*points[i*2]);
767 //Rectilinear is reinterpreted as uniform because ITK does not know rectilinear
769 for (i=0; i<GetNumberOfDimensions(); i++) {
770 //Compute mean spacing
771 SetSpacing(i,10*(points[GetDimensions(i)-1]-points[0])/(GetDimensions(i)-1));
772 SetOrigin(i,10*points[0]);
774 //Test if rectilinear image is actually uniform (tolerance 0.1 mm)
775 for (j=0; j<GetDimensions(i)-1; j++) {
776 if (fabs((points[j+1]-points[j])*10-GetSpacing(i))>0.1) {
779 return ER_NOT_HANDLED;
782 points += (int)GetDimensions(i);
784 for (i=0; i<GetNumberOfDimensions(); i++)
785 points -= GetDimensions(i);
790 return ER_NOT_HANDLED;
798 //====================================================================
799 // Read image information (copied from XDRreader)
800 void clitk::XdrImageIO::ReadImageInformation() {
801 int result = ReadImageInformationWithError();
802 if (result) ITKError("clitk::XdrImageIO::ReadImageInformation",result);
806 //====================================================================
807 // Read Image Content (copied from Xdr reader)
808 int clitk::XdrImageIO::ReadWithError(void * buffer)
810 int /*ndim,*/ nspace/*, veclen=1, data=AVS_TYPE_BYTE, field=UNIFORM*/;
811 int iNkiCompression = 0;
812 int j, coords=0, datasize=0, HeaderSize;
813 unsigned int i,iNumRead,total=1;
818 //AVSfield FieldTemplate;
819 long swap_test = 0x1000000; /* first byte is 1 when low-endian */
821 char *file = const_cast<char *>(m_FileName.c_str());
823 AVSType field=UNIFORM;
825 for (i=0; i<GetNumberOfDimensions(); i++) coords += GetDimensions(i);
827 total = GetImageSizeInPixels();
828 nspace = GetNumberOfDimensions();
830 c = scan_header(file, "field", offset, 1);
833 if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace*2;
834 else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
835 else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
836 else return ER_XDR_FIELD;
841 c = scan_header(file, "nki_compression", offset, 1);
842 if (c) iNkiCompression = atoi(c);
844 c = scan_header(file, "coord1[0]", offset, 1);
845 if (c) HeaderSize = 32768;
846 else HeaderSize = 2048;
848 fstream = fopen(file, "rb");
852 if (offset) fseek(fstream, offset, SEEK_SET);
856 if (fgets(temp, 500, fstream) == NULL )
857 return ER_XDR_NOCTRLL; /* end of file */
859 if (temp[0] == 10) continue;
863 fseek(fstream, -2, SEEK_CUR);
865 } /* ^L end of header */
867 if (temp[0] != '#') break;
870 buff = (char*)malloc(HeaderSize);
873 return ER_OUTOFMEMORY;
875 memset(buff, 0, HeaderSize);
876 iNumRead = fread(buff, 1, HeaderSize, fstream);
884 for (i=0; i<iNumRead; i++) {
885 if (buff[i] == 12) break;
890 if (i==iNumRead) return ER_XDR_NOCTRLL;
892 total = GetImageSizeInBytes();
894 //We add casts because the resulting quantity can be negative.
895 //There is no risk of looping because i and iNumRead are about the size of the header
896 fseek(fstream, static_cast<int>(i)+2-static_cast<int>(iNumRead), SEEK_CUR);
898 if (total && iNkiCompression)
902 signed char* pCompressed;
904 /* Read or guess the size of the compressed data */
905 iCurPos = ftell(fstream);
906 iSize = get_nki_compressed_size(fstream);
910 fseek(fstream, 0, SEEK_END);
911 iSize = ftell(fstream);
912 iSize = iSize - iCurPos - coords;
914 // Get compressed size from header if possible; else use uncompressed size as safe estimate
915 if (iSize>total && offset) iSize=total+8;
918 fseek(fstream, iCurPos, SEEK_SET);
920 /* Allocate space for the compressed pixels */
921 pCompressed = (signed char*)malloc(iSize);
925 return ER_OUTOFMEMORY;
928 /* Read the compressed pixels */
929 if (fread( (void *)pCompressed, 1, iSize, fstream ) != iSize)
935 if (!nki_private_decompress((short*)buffer, pCompressed, iSize))
938 return ER_DECOMPRESSION;
942 fseek(fstream, iCurPos + global_len, SEEK_SET);
951 if (fread( (void *)buffer, 1, total, fstream ) != total)
958 /* swap data if required (xdr is low-endian) */
960 datasize = GetComponentSize();
961 if (!(*(char *)(&swap_test)) && !forcenoswap)
966 for (i=0; i<total; i+=2)
973 else if (datasize==4)
976 for (i=0; i<total; i+=4)
986 else if (datasize==8)
989 for (i=0; i<total; i+=8)
1012 //====================================================================
1013 // Read image information (copied from Xdr reader)
1014 void clitk::XdrImageIO::Read(void * buffer) {
1015 int result = ReadWithError(buffer);
1016 if (result) ITKError("clitk::XdrImageIO::Read",result);
1019 //====================================================================
1020 // Read Image Information
1021 bool clitk::XdrImageIO::CanReadFile(const char* FileNameToRead)
1026 fstream = fopen(FileNameToRead, "rt");
1027 if (fstream == NULL)
1029 // AVSerror("Couldn't open file " << FileNameToRead);
1032 fgets(temp, 500, fstream);
1035 if (memcmp(temp, "# AVS", 5)==0)
1041 void clitk::XdrImageIO::ITKError(std::string funcName, int msgID) {
1042 itkExceptionMacro(<< "Error in " << funcName << ". Message: " << gl_ErrorMsg[msgID]);