common/clitkXdrImageIOReader.cxx

   1 /**
   2  * @file   clitkXdrImageIO.cxx
   3  * @author Simon Rit <simon.rit@gmail.com>
   4  * @date   Sun Jun  1 22:12:20 2008
   5  *
   6  * @brief
   7  *
   8  *
   9  */
  10
  11 #include "clitkXdrImageIO.h"
  12
  13 // std include
  14 #include <iostream>
  15 #include <fstream>
  16 #include <sstream>
  17
  18 //defines
  19 #define MAXDIM 5
  20 #ifdef _WIN32
  21 #ifdef memicmp
  22 #undef memicmp
  23 #endif
  24 #define memicmp _memicmp
  25 #else
  26 #define memicmp strncasecmp
  27 #endif
  28
  29
  30 // Based on a true story by the Nederlands Kanker Instituut (AVS_RXDR.CPP from the 20080414)
  31 #include <stdio.h>
  32 #include <string.h>
  33 #include <math.h>
  34 #include <stdlib.h>
  35 #include <ctype.h>
  36
  37 /************************************************************************/
  38 /*                    DEFINES, ENUMERATED TYPES AND CONSTANTS           */
  39 /************************************************************************/
  40 enum {
  41     OK,
  42     ER_ILLCOMMFUNCT,
  43     ER_INARGUMENTS,
  44     ER_XDR_NDIM,
  45     ER_XDR_DIM,
  46     ER_XDR_NSPACE,
  47     ER_XDR_VECLEN,
  48     ER_XDR_DATA,
  49     ER_XDR_FIELD,
  50     ER_XDR_OPEN,
  51     ER_XDR_NOCTRLL,
  52     ER_XDR_READ,
  53     ER_OUTOFMEMORY,
  54     ER_DECOMPRESSION,
  55     ER_NOT_HANDLED
  56 };
  57
  58 typedef struct
  59 {
  60     unsigned int iOrgSize;
  61     unsigned int iMode;
  62     unsigned int iCompressedSize;
  63     unsigned int iOrgCRC;
  64     unsigned int iCompressedCRC;        /* Excluding this header */
  65 } NKI_MODE2;
  66
  67
  68 /************************************************************************/
  69 /*                             GLOBAL VARIABLES                         */
  70 /************************************************************************/
  71 const char* gl_ErrorMsg[] = {
  72     "",
  73     "Command or function not supported in this way.",
  74     "Error in arguments",
  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",
  84     "Out of memory",
  85     "Decompression failed",
  86     "Format not handled by clitkXdrImageIO (RECTILINEAR or IRREGULAR field)"
  87 };
  88
  89
  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
 134 };
 135
 136
 137 /************************************************************************/
 138 /*                             MODULE FUNCTIONS                         */
 139 /************************************************************************/
 140
 141 /* simple XDR file reader */
 142
 143 /* help routine, scans XDR file header for keyword entry, returns NULL
 144    if not found, else returns pointer to rest of line
 145 */
 146
 147 static char *scan_header(const char *file, const char *name, int offset, int removespaces)
 148 {
 149     int i, j, iStringLength;
 150     static char temp[512];
 151     FILE *f;
 152     char *p, *q;
 153
 154     if ((f = fopen(file, "rt")) == NULL) return NULL;
 155     if (offset) fseek(f, offset, SEEK_SET);
 156
 157     for (i=0; i<200; )
 158     {
 159         if (fgets(temp, 500, f) == NULL      ) break;       /* end of file */
 160
 161         if (removespaces)
 162         {
 163             temp[500] = 0;
 164             p = q = temp;                       /* remove spaces */
 165             iStringLength = strlen(temp);
 166             for (j=0; j<iStringLength; j++)
 167                 if (*q!=' ' && *q!=8) *p++ = *q++;
 168                 else q++;
 169             *p++ = 0;
 170         }
 171
 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 */
 176
 177         p++;                                                /* match, skip = */
 178         if (p[strlen(p)-1] == '\n')                         /* remove \n */
 179             p[strlen(p)-1] = 0;
 180
 181         fclose (f);
 182         return p;
 183     }
 184
 185     fclose(f);
 186     return NULL;
 187 }
 188
 189 static int get_nki_compressed_size(FILE *f)
 190 {
 191     NKI_MODE2           Header;
 192     int                 iMode;
 193
 194     fread((void *)&Header, sizeof(Header), 1 , f);
 195
 196     iMode = Header.iMode;
 197
 198     switch (iMode)
 199     {
 200     case  1:
 201     case  3:
 202         return 0;
 203     case  2:
 204     case  4:
 205         return Header.iCompressedSize + sizeof(Header);
 206     default:
 207         return 0;
 208     }
 209 }
 210
 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)
 215
 216    The return value is the number of pixels that have been processed.
 217
 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)
 226
 227    Performance on typical CT or MRI is >2x compression and a very good speed
 228
 229    This code is not valid on HIGHENDIAN (high byte first) machines
 230 */
 231
 232 static int global_len;
 233
 234 static int nki_private_decompress(short int *dest, signed char *src, int size)
 235 {
 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;
 241
 242     retvalue = npixels = pHeader->iOrgSize;
 243     iMode = pHeader->iMode;             // safety: this value is checked in case statement
 244
 245     if (npixels<1) return 0;            // safety: check for invalid npixels value
 246
 247     /* Up till now only Mode=1, 2, 3, and 4 are supported */
 248
 249     switch (iMode)
 250     {
 251     case 1:
 252         save = src;
 253
 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
 256
 257         *dest = *(short int *)src;
 258         src += 2;
 259         npixels--;
 260
 261         do
 262         {
 263             if (src > end)                      // check whether the last few messages fit in input buffer
 264             {
 265                 if (src<end+3) val = *src;
 266                 else           val = 0;
 267
 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
 271                 else                                mode = 2;
 272
 273                 if (src+mode > end+3)
 274                     return 0;                   // safety: overflow input data
 275             }
 276
 277             val = *src;
 278
 279             if (val >= -64 && val <= 63)        // 7 bit difference
 280             {
 281                 dest[1] = dest[0] + val;
 282                 dest++;
 283                 src++;
 284             }
 285             else if (val==0x7f)         // 16 bit value
 286             {
 287                 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
 288                 dest++;
 289                 src+=3;
 290             }
 291             else if ((val&0xff)==0x80)  // run length encoding
 292             {
 293                 mode = ((unsigned char *)src)[1];
 294                 npixels -= mode-1;
 295                 if (npixels<=0) return 0;       // safety: overflow output data
 296                 do
 297                 {
 298                     dest[1] = dest[0];
 299                     dest++;
 300                 }
 301                 while (--mode);
 302                 src+=2;
 303             }
 304             else
 305             {
 306                 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
 307                 dest[1] = dest[0] + diff;       // 15 bit difference
 308                 dest++;
 309                 src+=2;
 310             }
 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         {
 335             if (src > end)                      // check whether the last few messages fit in input buffer
 336             {
 337                 if (src<end+3) val = *src;
 338                 else           val = 0;
 339
 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
 343                 else                                mode = 2;
 344
 345                 if (src+mode > end+3)
 346                     break;                      // safety: overflow input data
 347             }
 348
 349             val = *src;
 350
 351             if (val >= -64 && val <= 63)        // 7 bits difference
 352             {
 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));
 356                 dest++;
 357                 src++;
 358             }
 359             else if (val==0x7f)         // 16 bit value
 360             {
 361                 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
 362
 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));
 365                 dest++;
 366                 src+=3;
 367             }
 368             else if ((val&0xff)==0x80)  // run length encoding
 369             {
 370                 mode = ((unsigned char *)src)[1];
 371                 npixels -= mode-1;
 372                 if (npixels<=0) break;  // safety: overflow output data
 373                 do
 374                 {
 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));
 378                     dest++;
 379                 }
 380                 while (--mode);
 381                 src+=2;
 382             }
 383             else
 384             {
 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));
 389                 dest++;
 390                 src+=2;
 391             }
 392         }
 393         while (--npixels);
 394
 395         if (iCRC2 != pHeader->iOrgCRC)  // if error in output CRC:
 396         {
 397             src = save;                 // check input CRC
 398             while (src < end)
 399             {
 400                 iCRC = CRC32_table[(unsigned char)iCRC ^ (unsigned char)src[0]] ^ ((iCRC >> 8));
 401                 src++;
 402             }
 403
 404             if (iCRC != pHeader->iCompressedCRC)
 405             {
 406                 AVSerror("XDR decompression: the file is corrupted");
 407                 retvalue=0;
 408             }
 409             else
 410             {
 411                 AVSerror("XDR decompression: internal error");
 412                 retvalue=0;
 413             }
 414         }
 415
 416         global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
 417
 418         break;
 419
 420     case 3:
 421         save = src;
 422
 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
 425
 426         *dest = *(short int *)src;
 427         src += 2;
 428         npixels--;
 429
 430         do
 431         {
 432             if (src > end)                      // check whether the last few messages fit in input buffer
 433             {
 434                 if (src<end+3) val = *src;
 435                 else           val = 0;
 436
 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
 441                 else                                mode = 2;
 442
 443                 if (src+mode > end+3)
 444                     return 0;                   // safety: overflow input data
 445             }
 446
 447             val = *src;
 448
 449             if (val >= -63 && val <= 63)        // 7 bit difference
 450             {
 451                 dest[1] = dest[0] + val;
 452                 dest++;
 453                 src++;
 454             }
 455             else if (val==0x7f)         // 16 bit value
 456             {
 457                 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
 458                 dest++;
 459                 src+=3;
 460             }
 461             else if ((val&0xff)==0x80)  // run length encoding
 462             {
 463                 mode = ((unsigned char *)src)[1];
 464                 npixels -= mode-1;
 465                 if (npixels<=0) return 0;       // safety: overflow output data
 466                 do
 467                 {
 468                     dest[1] = dest[0];
 469                     dest++;
 470                 }
 471                 while (--mode);
 472                 src+=2;
 473             }
 474             else if ((val&0xff)==0xC0)  // 4 bit run
 475             {
 476                 mode = ((unsigned char *)src)[1];
 477                 npixels -= mode-1;
 478                 mode/=2;
 479                 src+=2;
 480                 if (npixels<=0) return 0;       // safety: overflow output data
 481                 do
 482                 {
 483                     val = *src++;
 484                     dest[1] = dest[0] + (val>>4);
 485                     dest++;
 486                     if (val&8) val |= 0xfffffff0;
 487                     else val &= 0x0f;
 488                     dest[1] = dest[0] + val;
 489                     dest++;
 490                 }
 491                 while (--mode);
 492             }
 493             else
 494             {
 495                 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
 496                 dest[1] = dest[0] + diff;       // 15 bit difference
 497                 dest++;
 498                 src+=2;
 499             }
 500         }
 501         while (--npixels);
 502
 503         global_len = src-save;
 504
 505         break;
 506
 507     case 4:
 508         src += sizeof(NKI_MODE2);
 509         save = src;
 510         end  = src + pHeader->iCompressedSize - 3;
 511
 512         if (end > src + size - 3)
 513             end = src + size - 3;               // may occur if pHeader is corrupted
 514
 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));
 518         src += 2;
 519         npixels--;
 520
 521         do
 522         {
 523             if (src > end)                      // check whether the last few messages fit in input buffer
 524             {
 525                 if (src<end+3) val = *src;
 526                 else           val = 0;
 527
 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
 532                 else                                mode = 2;
 533
 534                 if (src+mode > end+3)
 535                     return 0;                   // safety: overflow input data
 536             }
 537
 538             val = *src;
 539
 540             if (val >= -63 && val <= 63)        // 7 bit difference
 541             {
 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));
 545                 dest++;
 546                 src++;
 547             }
 548             else if (val==0x7f)         // 16 bit value
 549             {
 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));
 553                 dest++;
 554                 src+=3;
 555             }
 556             else if ((val&0xff)==0x80)  // run length encoding
 557             {
 558                 mode = ((unsigned char *)src)[1];
 559                 npixels -= mode-1;
 560                 if (npixels<=0) return 0;       // safety: overflow output data
 561                 do
 562                 {
 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));
 566                     dest++;
 567                 }
 568                 while (--mode);
 569                 src+=2;
 570             }
 571             else if ((val&0xff)==0xC0)  // 4 bit run
 572             {
 573                 mode = ((unsigned char *)src)[1];
 574                 npixels -= mode-1;
 575                 mode/=2;
 576                 src+=2;
 577                 if (npixels<=0) return 0;       // safety: overflow output data
 578                 do
 579                 {
 580                     val = *src++;
 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));
 584                     dest++;
 585                     if (val&8) val |= 0xfffffff0;
 586                     else val &= 0x0f;
 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));
 590                     dest++;
 591                 }
 592                 while (--mode);
 593             }
 594             else
 595             {
 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));
 600                 dest++;
 601                 src+=2;
 602             }
 603         }
 604         while (--npixels);
 605
 606         if (iCRC2 != pHeader->iOrgCRC)  // if error in output CRC:
 607             retvalue=0;
 608
 609         global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
 610
 611         break;
 612
 613
 614     default:
 615         AVSerror("XDR decompression: unsupported mode");
 616         return 0;
 617     }
 618
 619     return retvalue;
 620 }
 621
 622
 623 //====================================================================
 624 // Read image information (copied from XDRreader)
 625 int clitk::XdrImageIO::ReadImageInformationWithError()
 626 {
 627     int      offset=0;
 628     itk::Vector<int,MAXDIM> dim;
 629     int      veclen=1;
 630     int      total=1;
 631     unsigned int coords=0,i,j,ndim,nspace;
 632     char     temp[512];
 633     FILE     *fstream;
 634     char     *c;
 635
 636     long     swap_test = 0x1000000;      /* first byte is 1 when low-endian */
 637     forcenoswap=0;
 638     char     *file = const_cast<char *>(m_FileName.c_str());
 639     AVSType  field=UNIFORM;
 640
 641
 642     fstream = fopen(file, "rt");
 643     if (fstream == NULL) return ER_XDR_OPEN;
 644
 645     fgets(temp, 500, fstream);
 646     fclose(fstream);
 647
 648     if (memcmp(temp, "# AVS field file (produced by avs_nfwrite.c)", 44)==0) forcenoswap=1;
 649
 650     c      = scan_header(file, "ndim", offset, 1);
 651     if (!c) return ER_XDR_NDIM;
 652
 653     ndim   = atoi(c);
 654     if (ndim<1 || ndim>MAXDIM) return ER_XDR_NDIM;
 655     SetNumberOfDimensions(ndim);
 656
 657     nspace = ndim;
 658
 659     for (i=0; i<ndim; i++)
 660     {
 661         sprintf(temp, "dim%d", i+1);
 662         c = scan_header(file, temp, offset, 1);
 663         if (!c) return ER_XDR_DIM;
 664         dim[i]=atoi(c);
 665         if (dim[i]<1) return ER_XDR_DIM;
 666
 667         total  *= dim[i];
 668         coords += dim[i];
 669     }
 670     for (i=0; i<ndim; i++) {
 671         SetDimensions(i,dim[i]);
 672         SetSpacing(i,1.);
 673         SetOrigin(i,0.);
 674     }
 675
 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;
 680
 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);
 687
 688     c = scan_header(file, "data", offset, 1);
 689     if (c)
 690     {
 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;
 698
 699         if (memicmp(c, "xdr_",  4) == 0) forcenoswap=0;
 700     }
 701
 702     //Read coords here
 703     c = scan_header(file, "field", offset, 1);
 704     if (c)
 705     {
 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;
 710     }
 711     else
 712         coords=0;
 713
 714     if (coords)                        /* expect AVS coordinates ? */
 715     {
 716         coords *= sizeof(float);
 717         fstream = fopen(m_FileName.c_str(), "rb");
 718         if (fstream == NULL) return ER_XDR_OPEN;
 719
 720         float *points = (float *)malloc(coords);
 721         if (points == NULL) return ER_OUTOFMEMORY;
 722
 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)
 728             {
 729                 c = (char *)/*(*output)->*/points;
 730                 for (i=0; i<coords; i+=4)
 731                 {
 732                     j = c[i];
 733                     c[i]   = c[i+3];
 734                     c[i+3] = j;
 735                     j = c[i+1];
 736                     c[i+1] = c[i+2];
 737                     c[i+2] = j;
 738                 }
 739             }
 740         }
 741
 742         switch (field) {
 743         case UNIFORM:
 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]);
 747             }
 748             break;
 749         case RECTILINEAR:
 750             //Rectilinear is reinterpreted as uniform because ITK does not know rectilinear
 751             //Error if fails
 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]);
 756
 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) {
 760                         free(points);
 761                         fclose(fstream);
 762                         return ER_NOT_HANDLED;
 763                     }
 764                 }
 765                 points += (int)GetDimensions(i);
 766             }
 767             for (i=0; i<GetNumberOfDimensions(); i++)
 768                 points -= GetDimensions(i);
 769             break;
 770         case IRREGULAR:
 771             free(points);
 772             fclose(fstream);
 773             return ER_NOT_HANDLED;
 774         }
 775         free(points);
 776         fclose(fstream);
 777     }
 778     return OK;
 779 }
 780
 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);
 786 }
 787
 788
 789 //====================================================================
 790 // Read Image Content (copied from Xdr reader)
 791 int clitk::XdrImageIO::ReadWithError(void * buffer)
 792 { //AVSINT   dim[5];
 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;
 797     char     temp[512];
 798     FILE     *fstream;
 799     char     *c;
 800     char     *buff;
 801     //AVSfield FieldTemplate;
 802     long     swap_test = 0x1000000;      /* first byte is 1 when low-endian */
 803     //int      forcenoswap=0;
 804     char     *file = const_cast<char *>(m_FileName.c_str());
 805     int      offset=0;
 806     AVSType  field=UNIFORM;
 807
 808     for (i=0; i<GetNumberOfDimensions(); i++) coords += GetDimensions(i);
 809
 810     total = GetImageSizeInPixels();
 811     nspace = GetNumberOfDimensions();
 812
 813     c = scan_header(file, "field", offset, 1);
 814     if (c)
 815     {
 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;
 820     }
 821     else
 822         coords=0;
 823
 824     c = scan_header(file, "nki_compression", offset, 1);
 825     if (c) iNkiCompression = atoi(c);
 826
 827     c = scan_header(file, "coord1[0]", offset, 1);
 828     if (c) HeaderSize = 32768;
 829     else HeaderSize = 2048;
 830
 831     fstream = fopen(file, "rb");
 832     if (fstream == NULL)
 833         return ER_XDR_OPEN;
 834
 835     if (offset) fseek(fstream, offset, SEEK_SET);
 836
 837     while (1)
 838     {
 839         if (fgets(temp, 500, fstream) == NULL )
 840             return ER_XDR_NOCTRLL; /* end of file */
 841
 842         if (temp[0] == 10) continue;
 843
 844         if (temp[0] == 12)
 845         {
 846             fseek(fstream, -2, SEEK_CUR);
 847             break;
 848         } /* ^L end of header */
 849
 850         if (temp[0] != '#') break;
 851     }
 852
 853     buff = (char*)malloc(HeaderSize);
 854     if (buff == NULL)
 855     {
 856         return ER_OUTOFMEMORY;
 857     }
 858     memset(buff, 0, HeaderSize);
 859     iNumRead = fread(buff, 1, HeaderSize, fstream);
 860     if (iNumRead < 1)
 861     {
 862         free(buff);
 863         fclose(fstream);
 864         return ER_XDR_READ;
 865     }
 866
 867     for (i=0; i<iNumRead; i++) {
 868         if (buff[i] == 12) break;
 869     }
 870
 871     free(buff);
 872
 873     if (i==iNumRead) return ER_XDR_NOCTRLL;
 874
 875     total = GetImageSizeInBytes();
 876
 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);
 880
 881     if (total && iNkiCompression)
 882     {
 883         long            iCurPos;
 884         unsigned long iSize;
 885         signed char*    pCompressed;
 886
 887         /* Read or guess the size of the compressed data */
 888         iCurPos = ftell(fstream);
 889         iSize = get_nki_compressed_size(fstream);
 890
 891         if (iSize==0)
 892         {
 893             fseek(fstream, 0, SEEK_END);
 894             iSize = ftell(fstream);
 895             iSize = iSize - iCurPos - coords;
 896
 897             // Get compressed size from header if possible; else use uncompressed size as safe estimate
 898             if (iSize>total && offset) iSize=total+8;
 899         }
 900
 901         fseek(fstream, iCurPos, SEEK_SET);
 902
 903         /* Allocate space for the compressed pixels */
 904         pCompressed = (signed char*)malloc(iSize);
 905         if (!pCompressed)
 906         {
 907             fclose(fstream);
 908             return ER_OUTOFMEMORY;
 909         }
 910
 911         /* Read the compressed pixels */
 912         if (fread( (void *)pCompressed, 1, iSize, fstream ) != iSize)
 913         {
 914             fclose(fstream);
 915             return ER_XDR_READ;
 916         }
 917
 918         if (!nki_private_decompress((short*)buffer, pCompressed, iSize))
 919         {
 920             fclose(fstream);
 921             return ER_DECOMPRESSION;
 922         }
 923
 924         // if (offset)
 925         fseek(fstream, iCurPos + global_len, SEEK_SET);
 926
 927         free(pCompressed);
 928         goto READ_COORDS;
 929     }
 930
 931
 932     if (total)
 933     {
 934         if (fread( (void *)buffer, 1, total, fstream ) != total)
 935         {
 936             fclose(fstream);
 937             return ER_XDR_READ;
 938         }
 939     }
 940
 941     /* swap data if required (xdr is low-endian) */
 942
 943     datasize = GetComponentSize();
 944     if (!(*(char *)(&swap_test)) && !forcenoswap)
 945     {
 946         if (datasize==2)
 947         {
 948             c = (char *)buffer;
 949             for (i=0; i<total; i+=2)
 950             {
 951                 j = c[i];
 952                 c[i]   = c[i+1];
 953                 c[i+1] = j;
 954             }
 955         }
 956         else if (datasize==4)
 957         {
 958             c = (char *)buffer;
 959             for (i=0; i<total; i+=4)
 960             {
 961                 j = c[i];
 962                 c[i]   = c[i+3];
 963                 c[i+3] = j;
 964                 j = c[i+1];
 965                 c[i+1] = c[i+2];
 966                 c[i+2] = j;
 967             }
 968         }
 969         else if (datasize==8)
 970         {
 971             c = (char *)buffer;
 972             for (i=0; i<total; i+=8)
 973             {
 974                 j = c[i];
 975                 c[i]   = c[i+7];
 976                 c[i+7] = j;
 977                 j = c[i+1];
 978                 c[i+1] = c[i+6];
 979                 c[i+6] = j;
 980                 j = c[i+2];
 981                 c[i+2] = c[i+5];
 982                 c[i+5] = j;
 983                 j = c[i+3];
 984                 c[i+3] = c[i+4];
 985                 c[i+4] = j;
 986             }
 987         }
 988     }
 989
 990 READ_COORDS:
 991     fclose(fstream);
 992     return OK;
 993 }
 994
 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);
1000 }
1001
1002 //====================================================================
1003 // Read Image Information
1004 bool clitk::XdrImageIO::CanReadFile(const char* FileNameToRead)
1005 {
1006     char     temp[512];
1007     FILE     *fstream;
1008
1009     fstream = fopen(FileNameToRead, "rt");
1010     if (fstream == NULL)
1011     {
1012         AVSerror("Couldn't open file " << FileNameToRead);
1013         return false;
1014     }
1015     fgets(temp, 500, fstream);
1016     fclose(fstream);
1017
1018     if (memcmp(temp, "# AVS", 5)==0)
1019         return true;
1020     else
1021         return false;
1022 } ////
1023
1024 void clitk::XdrImageIO::ITKError(std::string funcName, int msgID) {
1025     itkExceptionMacro(<< "Error in " << funcName << ". Message: " << gl_ErrorMsg[msgID]);
1026 }
1027