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