common/clitkXdrImageIOReader.cxx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   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.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18 /**
  19  * @file   clitkXdrImageIO.cxx
  20  * @author Simon Rit <simon.rit@gmail.com>
  21  * @date   Sun Jun  1 22:12:20 2008
  22  *
  23  * @brief
  24  *
  25  *
  26  */
  27
  28 #include "clitkXdrImageIO.h"
  29
  30 // std include
  31 #include <iostream>
  32 #include <fstream>
  33 #include <sstream>
  34
  35 //defines
  36 #define MAXDIM 5
  37 #ifdef _WIN32
  38 #ifdef memicmp
  39 #undef memicmp
  40 #endif
  41 #define memicmp _memicmp
  42 #else
  43 #define memicmp strncasecmp
  44 #endif
  45
  46
  47 // Based on a true story by the Nederlands Kanker Instituut (AVS_RXDR.CPP from the 20080414)
  48 #include <stdio.h>
  49 #include <string.h>
  50 #include <math.h>
  51 #include <stdlib.h>
  52 #include <ctype.h>
  53
  54 /************************************************************************/
  55 /*                    DEFINES, ENUMERATED TYPES AND CONSTANTS           */
  56 /************************************************************************/
  57 enum {
  58   OK,
  59   ER_ILLCOMMFUNCT,
  60   ER_INARGUMENTS,
  61   ER_XDR_NDIM,
  62   ER_XDR_DIM,
  63   ER_XDR_NSPACE,
  64   ER_XDR_VECLEN,
  65   ER_XDR_DATA,
  66   ER_XDR_FIELD,
  67   ER_XDR_OPEN,
  68   ER_XDR_NOCTRLL,
  69   ER_XDR_READ,
  70   ER_OUTOFMEMORY,
  71   ER_DECOMPRESSION,
  72   ER_NOT_HANDLED
  73 };
  74
  75 typedef struct {
  76   unsigned int iOrgSize;
  77   unsigned int iMode;
  78   unsigned int iCompressedSize;
  79   unsigned int iOrgCRC;
  80   unsigned int iCompressedCRC;  /* Excluding this header */
  81 } NKI_MODE2;
  82
  83
  84 /************************************************************************/
  85 /*                             GLOBAL VARIABLES                         */
  86 /************************************************************************/
  87 const char* gl_ErrorMsg[] = {
  88   "",
  89   "Command or function not supported in this way.",
  90   "Error in arguments",
  91   "XDR file header NDIM error",
  92   "XDR file header DIMn error",
  93   "XDR file header NSPACE error",
  94   "XDR file header VECLEN error",
  95   "XDR file header DATA(type) error",
  96   "XDR file header FIELD(coordinate type) error",
  97   "XDR file could not be opened",
  98   "XDR file header contains no ^L",
  99   "XDR file reading error",
 100   "Out of memory",
 101   "Decompression failed",
 102   "Format not handled by clitkXdrImageIO (RECTILINEAR or IRREGULAR field)"
 103 };
 104
 105
 106 static const unsigned long CRC32_table[256] = {
 107   0x00000000, 0x77073096, 0xee0e612c, 0x990951ba, 0x076dc419, 0x706af48f,
 108   0xe963a535, 0x9e6495a3, 0x0edb8832, 0x79dcb8a4, 0xe0d5e91e, 0x97d2d988,
 109   0x09b64c2b, 0x7eb17cbd, 0xe7b82d07, 0x90bf1d91, 0x1db71064, 0x6ab020f2,
 110   0xf3b97148, 0x84be41de, 0x1adad47d, 0x6ddde4eb, 0xf4d4b551, 0x83d385c7,
 111   0x136c9856, 0x646ba8c0, 0xfd62f97a, 0x8a65c9ec, 0x14015c4f, 0x63066cd9,
 112   0xfa0f3d63, 0x8d080df5, 0x3b6e20c8, 0x4c69105e, 0xd56041e4, 0xa2677172,
 113   0x3c03e4d1, 0x4b04d447, 0xd20d85fd, 0xa50ab56b, 0x35b5a8fa, 0x42b2986c,
 114   0xdbbbc9d6, 0xacbcf940, 0x32d86ce3, 0x45df5c75, 0xdcd60dcf, 0xabd13d59,
 115   0x26d930ac, 0x51de003a, 0xc8d75180, 0xbfd06116, 0x21b4f4b5, 0x56b3c423,
 116   0xcfba9599, 0xb8bda50f, 0x2802b89e, 0x5f058808, 0xc60cd9b2, 0xb10be924,
 117   0x2f6f7c87, 0x58684c11, 0xc1611dab, 0xb6662d3d, 0x76dc4190, 0x01db7106,
 118   0x98d220bc, 0xefd5102a, 0x71b18589, 0x06b6b51f, 0x9fbfe4a5, 0xe8b8d433,
 119   0x7807c9a2, 0x0f00f934, 0x9609a88e, 0xe10e9818, 0x7f6a0dbb, 0x086d3d2d,
 120   0x91646c97, 0xe6635c01, 0x6b6b51f4, 0x1c6c6162, 0x856530d8, 0xf262004e,
 121   0x6c0695ed, 0x1b01a57b, 0x8208f4c1, 0xf50fc457, 0x65b0d9c6, 0x12b7e950,
 122   0x8bbeb8ea, 0xfcb9887c, 0x62dd1ddf, 0x15da2d49, 0x8cd37cf3, 0xfbd44c65,
 123   0x4db26158, 0x3ab551ce, 0xa3bc0074, 0xd4bb30e2, 0x4adfa541, 0x3dd895d7,
 124   0xa4d1c46d, 0xd3d6f4fb, 0x4369e96a, 0x346ed9fc, 0xad678846, 0xda60b8d0,
 125   0x44042d73, 0x33031de5, 0xaa0a4c5f, 0xdd0d7cc9, 0x5005713c, 0x270241aa,
 126   0xbe0b1010, 0xc90c2086, 0x5768b525, 0x206f85b3, 0xb966d409, 0xce61e49f,
 127   0x5edef90e, 0x29d9c998, 0xb0d09822, 0xc7d7a8b4, 0x59b33d17, 0x2eb40d81,
 128   0xb7bd5c3b, 0xc0ba6cad, 0xedb88320, 0x9abfb3b6, 0x03b6e20c, 0x74b1d29a,
 129   0xead54739, 0x9dd277af, 0x04db2615, 0x73dc1683, 0xe3630b12, 0x94643b84,
 130   0x0d6d6a3e, 0x7a6a5aa8, 0xe40ecf0b, 0x9309ff9d, 0x0a00ae27, 0x7d079eb1,
 131   0xf00f9344, 0x8708a3d2, 0x1e01f268, 0x6906c2fe, 0xf762575d, 0x806567cb,
 132   0x196c3671, 0x6e6b06e7, 0xfed41b76, 0x89d32be0, 0x10da7a5a, 0x67dd4acc,
 133   0xf9b9df6f, 0x8ebeeff9, 0x17b7be43, 0x60b08ed5, 0xd6d6a3e8, 0xa1d1937e,
 134   0x38d8c2c4, 0x4fdff252, 0xd1bb67f1, 0xa6bc5767, 0x3fb506dd, 0x48b2364b,
 135   0xd80d2bda, 0xaf0a1b4c, 0x36034af6, 0x41047a60, 0xdf60efc3, 0xa867df55,
 136   0x316e8eef, 0x4669be79, 0xcb61b38c, 0xbc66831a, 0x256fd2a0, 0x5268e236,
 137   0xcc0c7795, 0xbb0b4703, 0x220216b9, 0x5505262f, 0xc5ba3bbe, 0xb2bd0b28,
 138   0x2bb45a92, 0x5cb36a04, 0xc2d7ffa7, 0xb5d0cf31, 0x2cd99e8b, 0x5bdeae1d,
 139   0x9b64c2b0, 0xec63f226, 0x756aa39c, 0x026d930a, 0x9c0906a9, 0xeb0e363f,
 140   0x72076785, 0x05005713, 0x95bf4a82, 0xe2b87a14, 0x7bb12bae, 0x0cb61b38,
 141   0x92d28e9b, 0xe5d5be0d, 0x7cdcefb7, 0x0bdbdf21, 0x86d3d2d4, 0xf1d4e242,
 142   0x68ddb3f8, 0x1fda836e, 0x81be16cd, 0xf6b9265b, 0x6fb077e1, 0x18b74777,
 143   0x88085ae6, 0xff0f6a70, 0x66063bca, 0x11010b5c, 0x8f659eff, 0xf862ae69,
 144   0x616bffd3, 0x166ccf45, 0xa00ae278, 0xd70dd2ee, 0x4e048354, 0x3903b3c2,
 145   0xa7672661, 0xd06016f7, 0x4969474d, 0x3e6e77db, 0xaed16a4a, 0xd9d65adc,
 146   0x40df0b66, 0x37d83bf0, 0xa9bcae53, 0xdebb9ec5, 0x47b2cf7f, 0x30b5ffe9,
 147   0xbdbdf21c, 0xcabac28a, 0x53b39330, 0x24b4a3a6, 0xbad03605, 0xcdd70693,
 148   0x54de5729, 0x23d967bf, 0xb3667a2e, 0xc4614ab8, 0x5d681b02, 0x2a6f2b94,
 149   0xb40bbe37, 0xc30c8ea1, 0x5a05df1b, 0x2d02ef8d
 150 };
 151
 152
 153 /************************************************************************/
 154 /*                             MODULE FUNCTIONS                         */
 155 /************************************************************************/
 156
 157 /* simple XDR file reader */
 158
 159 /* help routine, scans XDR file header for keyword entry, returns NULL
 160    if not found, else returns pointer to rest of line
 161 */
 162
 163 static char *scan_header(const char *file, const char *name, int offset, int removespaces)
 164 {
 165   int i, j, iStringLength;
 166   static char temp[512];
 167   FILE *f;
 168   char *p, *q;
 169
 170   if ((f = fopen(file, "rt")) == NULL) return NULL;
 171   if (offset) fseek(f, offset, SEEK_SET);
 172
 173   for (i=0; i<200; ) {
 174     if (fgets(temp, 500, f) == NULL      ) break;       /* end of file */
 175
 176     if (removespaces) {
 177       temp[500] = 0;
 178       p = q = temp;                     /* remove spaces */
 179       iStringLength = strlen(temp);
 180       for (j=0; j<iStringLength; j++)
 181         if (*q!=' ' && *q!=8) *p++ = *q++;
 182         else q++;
 183       *p++ = 0;
 184     }
 185
 186     if (temp[0] == 12                    ) break;       /* ^L end of header */
 187     if (temp[0] != '#') i++;                            /* The first 200 non comment lines must be read before data is opened. */
 188     if ((p = strchr(temp+1, '=')) == NULL) continue;    /* no '=' */
 189     if (memicmp(temp, name, p-temp)      ) continue;    /* no match */
 190
 191     p++;                                                /* match, skip = */
 192     if (p[strlen(p)-1] == '\n')                         /* remove \n */
 193       p[strlen(p)-1] = 0;
 194
 195     fclose (f);
 196     return p;
 197   }
 198
 199   fclose(f);
 200   return NULL;
 201 }
 202
 203 static int get_nki_compressed_size(FILE *f)
 204 {
 205   NKI_MODE2             Header;
 206   int                   iMode;
 207
 208   fread((void *)&Header, sizeof(Header), 1 , f);
 209
 210   iMode = Header.iMode;
 211
 212   switch (iMode) {
 213   case  1:
 214   case  3:
 215     return 0;
 216   case  2:
 217   case  4:
 218     return Header.iCompressedSize + sizeof(Header);
 219   default:
 220     return 0;
 221   }
 222 }
 223
 224 /* decoder for NKI private compressed pixel data
 225    arguments: dest    = (in) points to area where destination data is written (short)
 226               src     = (in) points compressed source data (byte stream)
 227               size    = (in) number of bytes source data (safety)
 228
 229    The return value is the number of pixels that have been processed.
 230
 231    The compressed data looks like:
 232    (number of pixels)-1 times:
 233      OR 1 byte   = LL     7  bits signed (difference pixel[1] - pixel[0]);
 234      OR 2 bytes  = HHLL   15 bits signed (difference pixel[1] - pixel[0]) xored with 0x4000;
 235      OR 3 bytes  = 7FHHLL 16 bits absolute pixel data if 15 bits difference is exceeded
 236      OR 2 bytes  = 80NN   run length encode NN zero differences (max 255)
 237 for mode 3 and 4 added:
 238      OR 2 bytes  = CONN   encode NN 4 bit differences (max 255)
 239
 240    Performance on typical CT or MRI is >2x compression and a very good speed
 241
 242    This code is not valid on HIGHENDIAN (high byte first) machines
 243 */
 244
 245 static int global_len;
 246
 247 static int nki_private_decompress(short int *dest, signed char *src, int size)
 248 {
 249   int                   npixels, retvalue, mode, iMode, val, j;
 250   NKI_MODE2*            pHeader = (NKI_MODE2*)src;
 251   unsigned long         iCRC=0, iCRC2=0;
 252   //unsigned char*      pDestStart = (unsigned char*)dest;
 253   signed char           *save, *end;
 254
 255   retvalue = npixels = pHeader->iOrgSize;
 256   iMode = pHeader->iMode;               // safety: this value is checked in case statement
 257
 258   if (npixels<1) return 0;              // safety: check for invalid npixels value
 259
 260   /* Up till now only Mode=1, 2, 3, and 4 are supported */
 261
 262   switch (iMode) {
 263   case 1:
 264     save = src;
 265
 266     src += 8;                           // mode 1 only has 8 bytes header: iOrgSize and iMode
 267     end  = src + size - 3;              // for overflow check if we are close to end of input buffer
 268
 269     *dest = *(short int *)src;
 270     src += 2;
 271     npixels--;
 272
 273     do {
 274       if (src > end) {          // check whether the last few messages fit in input buffer
 275         if (src<end+3) val = *src;
 276         else           val = 0;
 277
 278         if (val >= -64 && val <= 63)      mode = 1;     // 7 bit difference
 279         else if (val==0x7f)                 mode = 3;   // 16 bit value
 280         else if ((val&0xff)==0x80)          mode = 2;   // run length encoding
 281         else                                mode = 2;
 282
 283         if (src+mode > end+3)
 284           return 0;                     // safety: overflow input data
 285       }
 286
 287       val = *src;
 288
 289       if (val >= -64 && val <= 63) {    // 7 bit difference
 290         dest[1] = dest[0] + val;
 291         dest++;
 292         src++;
 293       } else if (val==0x7f) {   // 16 bit value
 294         dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
 295         dest++;
 296         src+=3;
 297       } else if ((val&0xff)==0x80) {    // run length encoding
 298         mode = ((unsigned char *)src)[1];
 299         npixels -= mode-1;
 300         if (npixels<=0) return 0;       // safety: overflow output data
 301         do {
 302           dest[1] = dest[0];
 303           dest++;
 304         } while (--mode);
 305         src+=2;
 306       } else {
 307         signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
 308         dest[1] = dest[0] + diff;       // 15 bit difference
 309         dest++;
 310         src+=2;
 311       }
 312     } while (--npixels);
 313
 314     global_len = src-save;
 315
 316     break;
 317
 318   case 2:
 319     src += sizeof(NKI_MODE2);
 320     save = src;
 321     end  = src + pHeader->iCompressedSize - 3;
 322
 323     if (end > src + size - 3)
 324       end = src + size - 3;             // may occur if pHeader is corrupted
 325
 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));
 329     src+=2;
 330
 331     npixels--;
 332
 333     do {
 334       if (src > end) {          // check whether the last few messages fit in input buffer
 335         if (src<end+3) val = *src;
 336         else           val = 0;
 337
 338         if (val >= -64 && val <= 63)      mode = 1;     // 7 bit difference
 339         else if (val==0x7f)                 mode = 3;   // 16 bit value
 340         else if ((val&0xff)==0x80)          mode = 2;   // run length encoding
 341         else                                mode = 2;
 342
 343         if (src+mode > end+3)
 344           break;                        // safety: overflow input data
 345       }
 346
 347       val = *src;
 348
 349       if (val >= -64 && val <= 63) {    // 7 bits difference
 350         dest[1] = val = dest[0] + val;
 351         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 352         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 353         dest++;
 354         src++;
 355       } else if (val==0x7f) {   // 16 bit value
 356         dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
 357
 358         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 359         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 360         dest++;
 361         src+=3;
 362       } else if ((val&0xff)==0x80) {    // run length encoding
 363         mode = ((unsigned char *)src)[1];
 364         npixels -= mode-1;
 365         if (npixels<=0) break;  // safety: overflow output data
 366         do {
 367           dest[1] = val = dest[0];
 368           iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 369           iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 370           dest++;
 371         } while (--mode);
 372         src+=2;
 373       } else {
 374         signed short diff = ((val^0x40)<<8) + ((unsigned char *)src)[1];
 375         dest[1] = val = dest[0] + diff; // 15 bit difference
 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));
 378         dest++;
 379         src+=2;
 380       }
 381     } while (--npixels);
 382
 383     if (iCRC2 != pHeader->iOrgCRC) {    // if error in output CRC:
 384       src = save;                       // check input CRC
 385       while (src < end) {
 386         iCRC = CRC32_table[(unsigned char)iCRC ^ (unsigned char)src[0]] ^ ((iCRC >> 8));
 387         src++;
 388       }
 389
 390       if (iCRC != pHeader->iCompressedCRC) {
 391         AVSerror("XDR decompression: the file is corrupted");
 392         retvalue=0;
 393       } else {
 394         AVSerror("XDR decompression: internal error");
 395         retvalue=0;
 396       }
 397     }
 398
 399     global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
 400
 401     break;
 402
 403   case 3:
 404     save = src;
 405
 406     src += 8;                           // mode 3 only has 8 bytes header: iOrgSize and iMode
 407     end  = src + size - 3;              // for overflow check if we are close to end of input buffer
 408
 409     *dest = *(short int *)src;
 410     src += 2;
 411     npixels--;
 412
 413     do {
 414       if (src > end) {          // check whether the last few messages fit in input buffer
 415         if (src<end+3) val = *src;
 416         else           val = 0;
 417
 418         if (val >= -63 && val <= 63)      mode = 1;     // 7 bit difference
 419         else if (val==0x7f)                 mode = 3;   // 16 bit value
 420         else if ((val&0xff)==0x80)          mode = 2;   // run length encoding
 421         else if ((val&0xff)==0xC0)          mode = 2;   // 4 bit encoding
 422         else                                mode = 2;
 423
 424         if (src+mode > end+3)
 425           return 0;                     // safety: overflow input data
 426       }
 427
 428       val = *src;
 429
 430       if (val >= -63 && val <= 63) {    // 7 bit difference
 431         dest[1] = dest[0] + val;
 432         dest++;
 433         src++;
 434       } else if (val==0x7f) {   // 16 bit value
 435         dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
 436         dest++;
 437         src+=3;
 438       } else if ((val&0xff)==0x80) {    // run length encoding
 439         mode = ((unsigned char *)src)[1];
 440         npixels -= mode-1;
 441         if (npixels<=0) return 0;       // safety: overflow output data
 442         do {
 443           dest[1] = dest[0];
 444           dest++;
 445         } while (--mode);
 446         src+=2;
 447       } else if ((val&0xff)==0xC0) {    // 4 bit run
 448         mode = ((unsigned char *)src)[1];
 449         npixels -= mode-1;
 450         mode/=2;
 451         src+=2;
 452         if (npixels<=0) return 0;       // safety: overflow output data
 453         do {
 454           val = *src++;
 455           dest[1] = dest[0] + (val>>4);
 456           dest++;
 457           if (val&8) val |= 0xfffffff0;
 458           else val &= 0x0f;
 459           dest[1] = dest[0] + val;
 460           dest++;
 461         } while (--mode);
 462       } else {
 463         signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
 464         dest[1] = dest[0] + diff;       // 15 bit difference
 465         dest++;
 466         src+=2;
 467       }
 468     } while (--npixels);
 469
 470     global_len = src-save;
 471
 472     break;
 473
 474   case 4:
 475     src += sizeof(NKI_MODE2);
 476     save = src;
 477     end  = src + pHeader->iCompressedSize - 3;
 478
 479     if (end > src + size - 3)
 480       end = src + size - 3;             // may occur if pHeader is corrupted
 481
 482     *dest = val = *(short int *)src;
 483     iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 484     iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 485     src += 2;
 486     npixels--;
 487
 488     do {
 489       if (src > end) {          // check whether the last few messages fit in input buffer
 490         if (src<end+3) val = *src;
 491         else           val = 0;
 492
 493         if (val >= -63 && val <= 63)      mode = 1;     // 7 bit difference
 494         else if (val==0x7f)                 mode = 3;   // 16 bit value
 495         else if ((val&0xff)==0x80)          mode = 2;   // run length encoding
 496         else if ((val&0xff)==0xC0)          mode = 2;   // 4 bit encoding
 497         else                                mode = 2;
 498
 499         if (src+mode > end+3)
 500           return 0;                     // safety: overflow input data
 501       }
 502
 503       val = *src;
 504
 505       if (val >= -63 && val <= 63) {    // 7 bit difference
 506         dest[1] = val = dest[0] + val;
 507         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 508         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 509         dest++;
 510         src++;
 511       } else if (val==0x7f) {   // 16 bit value
 512         dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
 513         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 514         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 515         dest++;
 516         src+=3;
 517       } else if ((val&0xff)==0x80) {    // run length encoding
 518         mode = ((unsigned char *)src)[1];
 519         npixels -= mode-1;
 520         if (npixels<=0) return 0;       // safety: overflow output data
 521         do {
 522           dest[1] = val = dest[0];
 523           iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 524           iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 525           dest++;
 526         } while (--mode);
 527         src+=2;
 528       } else if ((val&0xff)==0xC0) {    // 4 bit run
 529         mode = ((unsigned char *)src)[1];
 530         npixels -= mode-1;
 531         mode/=2;
 532         src+=2;
 533         if (npixels<=0) return 0;       // safety: overflow output data
 534         do {
 535           val = *src++;
 536           dest[1] = j = dest[0] + (val>>4);
 537           iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
 538           iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
 539           dest++;
 540           if (val&8) val |= 0xfffffff0;
 541           else val &= 0x0f;
 542           dest[1] = j = dest[0] + val;
 543           iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
 544           iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
 545           dest++;
 546         } while (--mode);
 547       } else {
 548         signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
 549         dest[1] = val = dest[0] + diff; // 15 bit difference
 550         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 551         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 552         dest++;
 553         src+=2;
 554       }
 555     } while (--npixels);
 556
 557     if (iCRC2 != pHeader->iOrgCRC)      // if error in output CRC:
 558       retvalue=0;
 559
 560     global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
 561
 562     break;
 563
 564
 565   default:
 566     AVSerror("XDR decompression: unsupported mode");
 567     return 0;
 568   }
 569
 570   return retvalue;
 571 }
 572
 573
 574 //====================================================================
 575 // Read image information (copied from XDRreader)
 576 int clitk::XdrImageIO::ReadImageInformationWithError()
 577 {
 578   int      offset=0;
 579   itk::Vector<int,MAXDIM> dim;
 580   int      veclen=1;
 581   int      total=1;
 582   unsigned int coords=0,i,j,ndim,nspace;
 583   char     temp[512];
 584   FILE     *fstream;
 585   char     *c;
 586
 587   long     swap_test = 0x1000000;      /* first byte is 1 when low-endian */
 588   forcenoswap=0;
 589   char     *file = const_cast<char *>(m_FileName.c_str());
 590   AVSType  field=UNIFORM;
 591
 592
 593   fstream = fopen(file, "rt");
 594   if (fstream == NULL) return ER_XDR_OPEN;
 595
 596   fgets(temp, 500, fstream);
 597   fclose(fstream);
 598
 599   if (memcmp(temp, "# AVS field file (produced by avs_nfwrite.c)", 44)==0) forcenoswap=1;
 600
 601   c      = scan_header(file, "ndim", offset, 1);
 602   if (!c) return ER_XDR_NDIM;
 603
 604   ndim   = atoi(c);
 605   if (ndim<1 || ndim>MAXDIM) return ER_XDR_NDIM;
 606   SetNumberOfDimensions(ndim);
 607
 608   nspace = ndim;
 609
 610   for (i=0; i<ndim; i++) {
 611     sprintf(temp, "dim%d", i+1);
 612     c = scan_header(file, temp, offset, 1);
 613     if (!c) return ER_XDR_DIM;
 614     dim[i]=atoi(c);
 615     if (dim[i]<1) return ER_XDR_DIM;
 616
 617     total  *= dim[i];
 618     coords += dim[i];
 619   }
 620   for (i=0; i<ndim; i++) {
 621     SetDimensions(i,dim[i]);
 622     SetSpacing(i,1.);
 623     SetOrigin(i,0.);
 624   }
 625
 626   c = scan_header(file, "nspace", offset, 1);
 627   if (c) nspace = atoi(c);
 628   if (nspace<1 || ndim > MAXDIM) return ER_XDR_NSPACE;
 629   if (nspace != ndim) return ER_NOT_HANDLED;
 630
 631   c = scan_header(file, "veclen", offset, 1);
 632   if (c) veclen = atoi(c);
 633   if (veclen<0 /*|| veclen>1000*/) return ER_XDR_VECLEN;
 634   SetNumberOfComponents(veclen);
 635   if (veclen==1) SetPixelType(itk::ImageIOBase::SCALAR);
 636   else          SetPixelType(itk::ImageIOBase::VECTOR);
 637
 638   c = scan_header(file, "data", offset, 1);
 639   if (c) {
 640     if (memicmp(c, "byte",  4) == 0 || memicmp(c, "xdr_byte",  8) == 0) SetComponentType(itk::ImageIOBase::CHAR);
 641     else if (memicmp(c, "short", 5) == 0 || memicmp(c, "xdr_short", 9) == 0) SetComponentType(itk::ImageIOBase::SHORT);
 642     else if (memicmp(c, "int" ,  3) == 0 || memicmp(c, "xdr_int" ,  7) == 0) SetComponentType(itk::ImageIOBase::INT);
 643     else if (memicmp(c, "real",  4) == 0 || memicmp(c, "xdr_real",  8) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
 644     else if (memicmp(c, "float", 5) == 0 || memicmp(c, "xdr_float", 9) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
 645     else if (memicmp(c, "double",6) == 0 || memicmp(c, "xdr_double",10)== 0) SetComponentType(itk::ImageIOBase::DOUBLE);
 646     else return ER_XDR_DATA;
 647
 648     if (memicmp(c, "xdr_",  4) == 0) forcenoswap=0;
 649   }
 650
 651   //Read coords here
 652   c = scan_header(file, "field", offset, 1);
 653   if (c) {
 654     if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace *2;
 655     else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
 656     else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
 657     else return ER_XDR_FIELD;
 658   } else
 659     coords=0;
 660
 661   if (coords) {                      /* expect AVS coordinates ? */
 662     coords *= sizeof(float);
 663     fstream = fopen(m_FileName.c_str(), "rb");
 664     if (fstream == NULL) return ER_XDR_OPEN;
 665
 666     float *points = (float *)malloc(coords);
 667     if (points == NULL) return ER_OUTOFMEMORY;
 668
 669     //Seek to coordinates position in file
 670     if (fseek(fstream,-static_cast<int>(coords),SEEK_END)) return ER_XDR_READ;
 671     if (fread( /*(*output)->*/points, 1, coords, fstream ) == coords) {
 672       /* swap data if read-ok and required (xdr is low-endian) */
 673       if (!(*(char *)(&swap_test)) && !forcenoswap) {
 674         c = (char *)/*(*output)->*/points;
 675         for (i=0; i<coords; i+=4) {
 676           j = c[i];
 677           c[i]   = c[i+3];
 678           c[i+3] = j;
 679           j = c[i+1];
 680           c[i+1] = c[i+2];
 681           c[i+2] = j;
 682         }
 683       }
 684     }
 685
 686     float *p = points;
 687     switch (field) {
 688     case UNIFORM:
 689       for (i=0; i<GetNumberOfDimensions(); i++) {
 690         SetSpacing(i,10.*(points[i*2+1]-points[i*2])/(GetDimensions(i)-1));
 691         SetOrigin(i,10.*points[i*2]);
 692       }
 693       break;
 694     case RECTILINEAR:
 695       //Rectilinear is reinterpreted as uniform because ITK does not know rectilinear
 696       //Error if fails
 697       for (i=0; i<GetNumberOfDimensions(); i++) {
 698         //Compute mean spacing
 699         SetSpacing(i,10*(p[GetDimensions(i)-1]-p[0])/(GetDimensions(i)-1));
 700         SetOrigin(i,10*p[0]);
 701
 702         //Test if rectilinear image is actually uniform (tolerance 0.1 mm)
 703         if(i<3) { // Only for first 3 dimensions because spacing is barely used in other dims
 704           for (j=0; j<GetDimensions(i)-1; j++) {
 705             if (fabs((p[j+1]-p[j])*10-GetSpacing(i))>0.1) {
 706               free(points);
 707               fclose(fstream);
 708               return ER_NOT_HANDLED;
 709             }
 710           }
 711         }
 712         p += (int)GetDimensions(i);
 713       }
 714       break;
 715     case IRREGULAR:
 716       free(points);
 717       fclose(fstream);
 718       return ER_NOT_HANDLED;
 719     }
 720     free(points);
 721     fclose(fstream);
 722   }
 723   return OK;
 724 }
 725
 726 //====================================================================
 727 // Read image information (copied from XDRreader)
 728 void clitk::XdrImageIO::ReadImageInformation()
 729 {
 730   int result = ReadImageInformationWithError();
 731   if (result) ITKError("clitk::XdrImageIO::ReadImageInformation",result);
 732 }
 733
 734
 735 //====================================================================
 736 // Read Image Content (copied from Xdr reader)
 737 int clitk::XdrImageIO::ReadWithError(void * buffer)
 738 {
 739   //AVSINT   dim[5];
 740   int      /*ndim,*/ nspace/*, veclen=1, data=AVS_TYPE_BYTE, field=UNIFORM*/;
 741   int      iNkiCompression = 0;
 742   int      j, coords=0,  datasize=0, HeaderSize;
 743   unsigned int i,iNumRead,total=1;
 744   char     temp[512];
 745   FILE     *fstream;
 746   char     *c;
 747   char     *buff;
 748   //AVSfield FieldTemplate;
 749   long     swap_test = 0x1000000;      /* first byte is 1 when low-endian */
 750   //int      forcenoswap=0;
 751   char     *file = const_cast<char *>(m_FileName.c_str());
 752   int      offset=0;
 753   // AVSType  field=UNIFORM;
 754
 755   for (i=0; i<GetNumberOfDimensions(); i++) coords += GetDimensions(i);
 756
 757   total = GetImageSizeInPixels();
 758   nspace = GetNumberOfDimensions();
 759
 760   c = scan_header(file, "field", offset, 1);
 761   if (c) {
 762     if (memicmp(c, "unifo", 5) == 0) /*field=UNIFORM,*/ coords=nspace*2;
 763     else if (memicmp(c, "recti", 5) == 0) /*field=RECTILINEAR*/;
 764     else if (memicmp(c, "irreg", 5) == 0) /*field=IRREGULAR,*/ coords=total*nspace;
 765     else return ER_XDR_FIELD;
 766   } else
 767     coords=0;
 768
 769   c = scan_header(file, "nki_compression", offset, 1);
 770   if (c) iNkiCompression = atoi(c);
 771
 772   c = scan_header(file, "coord1[0]", offset, 1);
 773   if (c) HeaderSize = 32768;
 774   else HeaderSize = 2048;
 775
 776   fstream = fopen(file, "rb");
 777   if (fstream == NULL)
 778     return ER_XDR_OPEN;
 779
 780   if (offset) fseek(fstream, offset, SEEK_SET);
 781
 782   while (1) {
 783     if (fgets(temp, 500, fstream) == NULL )
 784       return ER_XDR_NOCTRLL; /* end of file */
 785
 786     if (temp[0] == 10) continue;
 787
 788     if (temp[0] == 12) {
 789       fseek(fstream, -2, SEEK_CUR);
 790       break;
 791     } /* ^L end of header */
 792
 793     if (temp[0] != '#') break;
 794   }
 795
 796   buff = (char*)malloc(HeaderSize);
 797   if (buff == NULL) {
 798     return ER_OUTOFMEMORY;
 799   }
 800   memset(buff, 0, HeaderSize);
 801   iNumRead = fread(buff, 1, HeaderSize, fstream);
 802   if (iNumRead < 1) {
 803     free(buff);
 804     fclose(fstream);
 805     return ER_XDR_READ;
 806   }
 807
 808   for (i=0; i<iNumRead; i++) {
 809     if (buff[i] == 12) break;
 810   }
 811
 812   free(buff);
 813
 814   if (i==iNumRead) return ER_XDR_NOCTRLL;
 815
 816   total = GetImageSizeInBytes();
 817
 818   //We add casts because the resulting quantity can be negative.
 819   //There is no risk of looping because i and iNumRead are about the size of the header
 820   fseek(fstream, static_cast<int>(i)+2-static_cast<int>(iNumRead), SEEK_CUR);
 821
 822   if (total && iNkiCompression) {
 823     long                iCurPos;
 824     unsigned long iSize;
 825     signed char*        pCompressed;
 826
 827     /* Read or guess the size of the compressed data */
 828     iCurPos = ftell(fstream);
 829     iSize = get_nki_compressed_size(fstream);
 830
 831     if (iSize==0) {
 832       fseek(fstream, 0, SEEK_END);
 833       iSize = ftell(fstream);
 834       iSize = iSize - iCurPos - coords;
 835
 836       // Get compressed size from header if possible; else use uncompressed size as safe estimate
 837       if (iSize>total && offset) iSize=total+8;
 838     }
 839
 840     fseek(fstream, iCurPos, SEEK_SET);
 841
 842     /* Allocate space for the compressed pixels */
 843     pCompressed = (signed char*)malloc(iSize);
 844     if (!pCompressed) {
 845       fclose(fstream);
 846       return ER_OUTOFMEMORY;
 847     }
 848
 849     /* Read the compressed pixels */
 850     if (fread( (void *)pCompressed, 1, iSize, fstream ) != iSize) {
 851       fclose(fstream);
 852       return ER_XDR_READ;
 853     }
 854
 855     if (!nki_private_decompress((short*)buffer, pCompressed, iSize)) {
 856       fclose(fstream);
 857       return ER_DECOMPRESSION;
 858     }
 859
 860     // if (offset)
 861     fseek(fstream, iCurPos + global_len, SEEK_SET);
 862
 863     free(pCompressed);
 864     goto READ_COORDS;
 865   }
 866
 867
 868   if (total) {
 869     if (fread( (void *)buffer, 1, total, fstream ) != total) {
 870       fclose(fstream);
 871       return ER_XDR_READ;
 872     }
 873   }
 874
 875   /* swap data if required (xdr is low-endian) */
 876
 877   datasize = GetComponentSize();
 878   if (!(*(char *)(&swap_test)) && !forcenoswap) {
 879     if (datasize==2) {
 880       c = (char *)buffer;
 881       for (i=0; i<total; i+=2) {
 882         j = c[i];
 883         c[i]   = c[i+1];
 884         c[i+1] = j;
 885       }
 886     } else if (datasize==4) {
 887       c = (char *)buffer;
 888       for (i=0; i<total; i+=4) {
 889         j = c[i];
 890         c[i]   = c[i+3];
 891         c[i+3] = j;
 892         j = c[i+1];
 893         c[i+1] = c[i+2];
 894         c[i+2] = j;
 895       }
 896     } else if (datasize==8) {
 897       c = (char *)buffer;
 898       for (i=0; i<total; i+=8) {
 899         j = c[i];
 900         c[i]   = c[i+7];
 901         c[i+7] = j;
 902         j = c[i+1];
 903         c[i+1] = c[i+6];
 904         c[i+6] = j;
 905         j = c[i+2];
 906         c[i+2] = c[i+5];
 907         c[i+5] = j;
 908         j = c[i+3];
 909         c[i+3] = c[i+4];
 910         c[i+4] = j;
 911       }
 912     }
 913   }
 914
 915 READ_COORDS:
 916   fclose(fstream);
 917   return OK;
 918 }
 919
 920 //====================================================================
 921 // Read image information (copied from Xdr reader)
 922 void clitk::XdrImageIO::Read(void * buffer)
 923 {
 924   int result = ReadWithError(buffer);
 925   if (result) ITKError("clitk::XdrImageIO::Read",result);
 926 }
 927
 928 //====================================================================
 929 // Read Image Information
 930 bool clitk::XdrImageIO::CanReadFile(const char* FileNameToRead)
 931 {
 932   char     temp[512];
 933   FILE     *fstream;
 934
 935   fstream = fopen(FileNameToRead, "rt");
 936   if (fstream == NULL)
 937 //    {
 938 //        AVSerror("Couldn't open file " << FileNameToRead);
 939     return false;
 940 //    }
 941   fgets(temp, 500, fstream);
 942   fclose(fstream);
 943
 944   if (memcmp(temp, "# AVS", 5)==0)
 945     return true;
 946   else
 947     return false;
 948 } ////
 949
 950 void clitk::XdrImageIO::ITKError(std::string funcName, int msgID)
 951 {
 952   itkExceptionMacro(<< "Error in " << funcName << ". Message: " << gl_ErrorMsg[msgID]);
 953 }
 954