vv/nkitkXDRImageIOReader.cxx

   1 #ifndef NKITKXDRIMAGEIO_CXX
   2 #define NKITKXDRIMAGEIO_CXX
   3
   4 /**
   5  * @file   nkitkXDRImageIO.cxx
   6  * @author Simon Rit <simon.rit@gmail.com>
   7  * @date   Sun Jun  1 22:12:20 2008
   8  *
   9  * @brief
  10  *
  11  *
  12  */
  13
  14 #include "nkitkXDRImageIO.h"
  15
  16 // std include
  17 #include <iostream>
  18 #include <fstream>
  19 #include <sstream>
  20
  21 //defines
  22 #define MAXDIM 5
  23 #define AVSerror(v) std::cerr << "Error in nkitk::XDRImageIO. Message:" << v << std::endl;
  24 #ifdef _WIN32
  25 #ifdef memicmp
  26 #undef memicmp
  27 #endif
  28 #define memicmp _memicmp
  29 #else
  30 #define memicmp strncasecmp
  31 #endif
  32
  33 /************************************************************************/
  34 /*                                                                      */
  35 /*      file       : AVS_RXDR.CPP                                       */
  36 /*                                                                      */
  37 /*      purpose    : AVS module for reading XDRs                        */
  38 /*                                                                      */
  39 /*      authors    : Lambert Zijp (Conquest)                            */
  40 /*                   Based on a true story by Marcel van Herk           */
  41 /*                                                                      */
  42 /*      date       : 19980528                                           */
  43 /*                                                                      */
  44 /*      portability: AVS requires sizeof(void *)==sizeof(int)           */
  45 /*                   This module assumes sizeof(int)>=4                 */
  46 /*                                                                      */
  47 /*      notes      :                                                    */
  48 /*                                                                      */
  49 /*                                                                      */
  50 /************************************************************************/
  51 /* Updates:
  52 When       Who   What
  53 19980528   ljz   Copied from mbfield1
  54 19981209   mvh   Taken into action by removing _; added READ_XDR_HEADER
  55 19981210 mvh+lsp Now proper freeing of buff in case of error
  56 19990130   mvh   Added ENUM_XDR_HEADER
  57 19990206   mvh   Added READ_XDR_PREVIEW
  58 19990208   mvh   Fixed this, the while(1) fgets loop is rquired to skip long headers
  59 19990406   ljz   Removed sizelimit of dimensions
  60 19991003   mvh   Loosened sizelimit for veclen to 1000
  61 19991005   mvh   Remove spaces from AVS own header lines (ndim etc) in scan_header
  62 20000808   lsp   No longer mix file streams and handles
  63 20000821   ljz   Reader can now handle NkiCompression
  64 20000822 lsp+nd  Initialize data and datasize for dots, added "const" to avoid
  65                  writing into a literal string
  66 20000823   mvh   Optimized nki_private_decompression mode 2 (from 6.4 to 4.9 s)
  67                  reading mode 1: compression makes read faster
  68                  reading mode 2: read compressed at same speed as uncompressed
  69 20000824   mvh   Speed up mode 2 to 4.2 s: only check compressedCRC if other failed
  70                  Added some safeties to avoid mode 1 or 2 crash for corrupted data
  71 20000825   mvh   Pass buffer size to nki_private_decompress as extra safety
  72                  Fix mode 1, added full safety against input buffer overflow
  73                  Added some comments and notes
  74 20010507   mvh   Support larger headers when coordinate information has been written
  75 20010726   mvh   Fix for decompression when information is offset in file
  76 20011207   ljz   Removed check on veclen>1000
  77 20020124   mvh   Added test for fields by kg: non-portable types as "integer" then not swapped
  78 20020903   ljz   Made scan_header much faster
  79 20030122   mvh   Fix read of coords in compression mode 2
  80 20030430   mvh   Fix of read coords in compression mode 2 when offset is 0
  81 20030717   ljz   Added support for NkiCompressionModes 3 and 4
  82 20040426   mvh   ELEKTA NKI-XVI0.1 RELEASE
  83 20050302   mvh   Estimate size of compressed data to accelerate speed of reading embedded fields
  84 20050308 ljz+mvh ELEKTA NKI-XVI0.1j RELEASE
  85 20071024        mvh     Adapted for 64 bits
  86 20080110        mvh     ELEKTA NKI-XVI3.09 RELEASE
  87 20080414 lsp+mgw __sun__ doesn't know <io.h>
  88 */
  89
  90 /************************************************************************/
  91 /*                         MODULE DOCUMENTATION                         */
  92 /************************************************************************/
  93 /*
  94    READ_XDR                     Read XDR file (may be compressed) into field
  95         name                    output field handle
  96         file_expression         name of file
  97         numerical_expression    start XDR data in file offset (default 0)
  98         NOTE: compressed XDR data may not be part of a larger file
  99
 100    READ_XDR_HEADER              Get entry from xdr header
 101         %name                   Output = string block
 102         file_expression         file name default extension ''
 103         %string_expression      name of element to load
 104
 105    READ_XDR_PREVIEW             Read and downsize XDR file (for bitmap)
 106         name                    output field handle
 107         file_expression         name of file
 108         numerical_expression    start XDR data in file offset (default 0)
 109         NOTE: this command is not supported for compressed XDR files
 110
 111    ENUM_XDR_HEADER
 112         %name                   Output = string block
 113         file_expression         file name default extension ''
 114         numerical_expression    number of element to find name of
 115 */
 116 /************************************************************************/
 117 /*                             INCLUDE FILES                            */
 118 /************************************************************************/
 119
 120 #include <stdio.h>
 121 #include <string.h>
 122 #include <math.h>
 123 #include <stdlib.h>
 124 #include <ctype.h>
 125 #if 0
 126 #ifndef __sun__
 127 #include <io.h>
 128 #endif
 129 #include "mbavs2q.h"
 130
 131 #ifndef QUIRT
 132 #include <memory.h>
 133 #include <avs/avs.h>
 134 #include <avs/field.h>
 135 #else
 136 #include "mbfield.h"
 137 #endif
 138 #endif
 139
 140 /************************************************************************/
 141 /*                    DEFINES, ENUMERATED TYPES AND CONSTANTS           */
 142 /************************************************************************/
 143 enum {
 144     OK,
 145     ER_ILLCOMMFUNCT,
 146     ER_INARGUMENTS,
 147     ER_XDR_NDIM,
 148     ER_XDR_DIM,
 149     ER_XDR_NSPACE,
 150     ER_XDR_VECLEN,
 151     ER_XDR_DATA,
 152     ER_XDR_FIELD,
 153     ER_XDR_OPEN,
 154     ER_XDR_NOCTRLL,
 155     ER_XDR_READ,
 156     ER_OUTOFMEMORY,
 157     ER_DECOMPRESSION,
 158     ER_NOT_HANDLED
 159 };
 160
 161 typedef struct
 162 {
 163     unsigned int iOrgSize;
 164     unsigned int iMode;
 165     unsigned int iCompressedSize;
 166     unsigned int iOrgCRC;
 167     unsigned int iCompressedCRC;        /* Excluding this header */
 168 } NKI_MODE2;
 169
 170
 171 /************************************************************************/
 172 /*                             GLOBAL VARIABLES                         */
 173 /************************************************************************/
 174 const char* gl_ErrorMsg[] = {
 175     "",
 176     "Command or function not supported in this way.",
 177     "Error in arguments",
 178     "XDR file header NDIM error",
 179     "XDR file header DIMn error",
 180     "XDR file header NSPACE error",
 181     "XDR file header VECLEN error",
 182     "XDR file header DATA(type) error",
 183     "XDR file header FIELD(coordinate type) error",
 184     "XDR file could not be opened",
 185     "XDR file header contains no ^L",
 186     "XDR file reading error",
 187     "Out of memory",
 188     "Decompression failed",
 189     "Format not handled by nkitkXDRImageIO (RECTILINEAR or IRREGULAR field)"
 190 };
 191
 192
 193 static const unsigned long CRC32_table[256] = {
 194     0x00000000, 0x77073096, 0xee0e612c, 0x990951ba, 0x076dc419, 0x706af48f,
 195     0xe963a535, 0x9e6495a3, 0x0edb8832, 0x79dcb8a4, 0xe0d5e91e, 0x97d2d988,
 196     0x09b64c2b, 0x7eb17cbd, 0xe7b82d07, 0x90bf1d91, 0x1db71064, 0x6ab020f2,
 197     0xf3b97148, 0x84be41de, 0x1adad47d, 0x6ddde4eb, 0xf4d4b551, 0x83d385c7,
 198     0x136c9856, 0x646ba8c0, 0xfd62f97a, 0x8a65c9ec, 0x14015c4f, 0x63066cd9,
 199     0xfa0f3d63, 0x8d080df5, 0x3b6e20c8, 0x4c69105e, 0xd56041e4, 0xa2677172,
 200     0x3c03e4d1, 0x4b04d447, 0xd20d85fd, 0xa50ab56b, 0x35b5a8fa, 0x42b2986c,
 201     0xdbbbc9d6, 0xacbcf940, 0x32d86ce3, 0x45df5c75, 0xdcd60dcf, 0xabd13d59,
 202     0x26d930ac, 0x51de003a, 0xc8d75180, 0xbfd06116, 0x21b4f4b5, 0x56b3c423,
 203     0xcfba9599, 0xb8bda50f, 0x2802b89e, 0x5f058808, 0xc60cd9b2, 0xb10be924,
 204     0x2f6f7c87, 0x58684c11, 0xc1611dab, 0xb6662d3d, 0x76dc4190, 0x01db7106,
 205     0x98d220bc, 0xefd5102a, 0x71b18589, 0x06b6b51f, 0x9fbfe4a5, 0xe8b8d433,
 206     0x7807c9a2, 0x0f00f934, 0x9609a88e, 0xe10e9818, 0x7f6a0dbb, 0x086d3d2d,
 207     0x91646c97, 0xe6635c01, 0x6b6b51f4, 0x1c6c6162, 0x856530d8, 0xf262004e,
 208     0x6c0695ed, 0x1b01a57b, 0x8208f4c1, 0xf50fc457, 0x65b0d9c6, 0x12b7e950,
 209     0x8bbeb8ea, 0xfcb9887c, 0x62dd1ddf, 0x15da2d49, 0x8cd37cf3, 0xfbd44c65,
 210     0x4db26158, 0x3ab551ce, 0xa3bc0074, 0xd4bb30e2, 0x4adfa541, 0x3dd895d7,
 211     0xa4d1c46d, 0xd3d6f4fb, 0x4369e96a, 0x346ed9fc, 0xad678846, 0xda60b8d0,
 212     0x44042d73, 0x33031de5, 0xaa0a4c5f, 0xdd0d7cc9, 0x5005713c, 0x270241aa,
 213     0xbe0b1010, 0xc90c2086, 0x5768b525, 0x206f85b3, 0xb966d409, 0xce61e49f,
 214     0x5edef90e, 0x29d9c998, 0xb0d09822, 0xc7d7a8b4, 0x59b33d17, 0x2eb40d81,
 215     0xb7bd5c3b, 0xc0ba6cad, 0xedb88320, 0x9abfb3b6, 0x03b6e20c, 0x74b1d29a,
 216     0xead54739, 0x9dd277af, 0x04db2615, 0x73dc1683, 0xe3630b12, 0x94643b84,
 217     0x0d6d6a3e, 0x7a6a5aa8, 0xe40ecf0b, 0x9309ff9d, 0x0a00ae27, 0x7d079eb1,
 218     0xf00f9344, 0x8708a3d2, 0x1e01f268, 0x6906c2fe, 0xf762575d, 0x806567cb,
 219     0x196c3671, 0x6e6b06e7, 0xfed41b76, 0x89d32be0, 0x10da7a5a, 0x67dd4acc,
 220     0xf9b9df6f, 0x8ebeeff9, 0x17b7be43, 0x60b08ed5, 0xd6d6a3e8, 0xa1d1937e,
 221     0x38d8c2c4, 0x4fdff252, 0xd1bb67f1, 0xa6bc5767, 0x3fb506dd, 0x48b2364b,
 222     0xd80d2bda, 0xaf0a1b4c, 0x36034af6, 0x41047a60, 0xdf60efc3, 0xa867df55,
 223     0x316e8eef, 0x4669be79, 0xcb61b38c, 0xbc66831a, 0x256fd2a0, 0x5268e236,
 224     0xcc0c7795, 0xbb0b4703, 0x220216b9, 0x5505262f, 0xc5ba3bbe, 0xb2bd0b28,
 225     0x2bb45a92, 0x5cb36a04, 0xc2d7ffa7, 0xb5d0cf31, 0x2cd99e8b, 0x5bdeae1d,
 226     0x9b64c2b0, 0xec63f226, 0x756aa39c, 0x026d930a, 0x9c0906a9, 0xeb0e363f,
 227     0x72076785, 0x05005713, 0x95bf4a82, 0xe2b87a14, 0x7bb12bae, 0x0cb61b38,
 228     0x92d28e9b, 0xe5d5be0d, 0x7cdcefb7, 0x0bdbdf21, 0x86d3d2d4, 0xf1d4e242,
 229     0x68ddb3f8, 0x1fda836e, 0x81be16cd, 0xf6b9265b, 0x6fb077e1, 0x18b74777,
 230     0x88085ae6, 0xff0f6a70, 0x66063bca, 0x11010b5c, 0x8f659eff, 0xf862ae69,
 231     0x616bffd3, 0x166ccf45, 0xa00ae278, 0xd70dd2ee, 0x4e048354, 0x3903b3c2,
 232     0xa7672661, 0xd06016f7, 0x4969474d, 0x3e6e77db, 0xaed16a4a, 0xd9d65adc,
 233     0x40df0b66, 0x37d83bf0, 0xa9bcae53, 0xdebb9ec5, 0x47b2cf7f, 0x30b5ffe9,
 234     0xbdbdf21c, 0xcabac28a, 0x53b39330, 0x24b4a3a6, 0xbad03605, 0xcdd70693,
 235     0x54de5729, 0x23d967bf, 0xb3667a2e, 0xc4614ab8, 0x5d681b02, 0x2a6f2b94,
 236     0xb40bbe37, 0xc30c8ea1, 0x5a05df1b, 0x2d02ef8d
 237 };
 238
 239
 240 /************************************************************************/
 241 /*                              PROTOTYPES                              */
 242 /************************************************************************/
 243 #if 0
 244 int Rxdr_compute(AVSfield** ppOut, char* pszFileName, int iOffset);
 245
 246 int RxdrPreview_compute(AVSfield** ppOut, char* pszFileName, int iOffset);
 247
 248 int RxdrHeader_compute(char *pszOut, char *pszFileName, char *pszEntry);
 249
 250 int RxdrEnum_compute(char *pszOut, char *pszFileName, int iEntry);
 251
 252 /******************************************************************************/
 253 /*                          AVS/QUIRT INTERFACE                               */
 254 /******************************************************************************/
 255
 256 void Rxdr_desc(void)
 257 {
 258     int  param;
 259
 260     /* Set the module name and type */
 261     AVSset_module_name("XDR reader", MODULE_DATA);
 262
 263     /* Create output ports for the resulting fields */
 264     AVScreate_output_port("Output", "field");
 265
 266     /* declare widgets */
 267     QUIRT_NEXT_PARAMETER_FILE("");
 268     param = AVSadd_parameter("FileName", "string", "/data/AVS/", NULL, "");
 269     AVSconnect_widget(param, "browser");
 270     AVSadd_parameter_prop(param, "height", "integer", 8);
 271
 272     param = AVSadd_parameter("Offset", "integer", 0, 0, INT_UNBOUND);
 273     AVSconnect_widget(param, "typein_integer");
 274
 275     AVSset_compute_proc(CF Rxdr_compute);
 276 }
 277 AVS_TO_QUIRT(READ_XDR, Rxdr_desc);
 278
 279 void RxdrPreview_desc(void)
 280 {
 281     int  param;
 282
 283     /* Set the module name and type */
 284     AVSset_module_name("XDR preview reader", MODULE_DATA);
 285
 286     /* Create output ports for the resulting fields */
 287     AVScreate_output_port("Output", "field");
 288
 289     /* declare widgets */
 290     QUIRT_NEXT_PARAMETER_FILE("");
 291     param = AVSadd_parameter("FileName", "string", "/data/AVS/", NULL, "");
 292     AVSconnect_widget(param, "browser");
 293     AVSadd_parameter_prop(param, "height", "integer", 8);
 294
 295     param = AVSadd_parameter("Offset", "integer", 0, 0, INT_UNBOUND);
 296     AVSconnect_widget(param, "typein_integer");
 297
 298     AVSset_compute_proc(CF RxdrPreview_compute);
 299 }
 300 AVS_TO_QUIRT(READ_XDR_PREVIEW, RxdrPreview_desc);
 301
 302 void RxdrHeader_desc(void)
 303 { /* Set the module name and type */
 304     AVSset_module_name("Read XDR header", MODULE_DATA);
 305
 306     /* Create output ports for the resulting fields */
 307     AVSadd_parameter("Output", "string_block", "", NULL, "");
 308
 309     QUIRT_NEXT_PARAMETER_FILE("");
 310     AVSadd_parameter("FileName", "string", "", NULL, "");
 311
 312     AVSadd_parameter("Entry", "string", "", NULL, "");
 313
 314     AVSset_compute_proc(CF RxdrHeader_compute);
 315 }
 316 AVS_TO_QUIRT(READ_XDR_HEADER, RxdrHeader_desc);
 317
 318 void RxdrEnum_desc(void)
 319 { /* Set the module name and type */
 320     AVSset_module_name("Enumerate XDR header", MODULE_DATA);
 321
 322     /* Create output ports for the resulting fields */
 323     AVSadd_parameter("Output", "string_block", "", NULL, "");
 324
 325     QUIRT_NEXT_PARAMETER_FILE("");
 326     AVSadd_parameter("FileName", "string", "", NULL, "");
 327
 328     AVSadd_parameter("Entry", "integer", 0, INT_UNBOUND, INT_UNBOUND);
 329
 330     AVSset_compute_proc(CF RxdrEnum_compute);
 331 }
 332 AVS_TO_QUIRT(ENUM_XDR_HEADER, RxdrEnum_desc);
 333
 334 #ifndef QUIRT
 335 AVSinit_modules(void)
 336 {
 337     AVSmodule_from_desc( (int_desc_func)Rxdr_desc );
 338     AVSmodule_from_desc( (int_desc_func)RxdrHeader_desc );
 339     AVSmodule_from_desc( (int_desc_func)RxdrEnum_desc );
 340 }
 341 #endif
 342 #endif
 343 /************************************************************************/
 344 /*                             MODULE FUNCTIONS                         */
 345 /************************************************************************/
 346
 347 /* simple XDR file reader */
 348
 349 /* help routine, scans XDR file header for keyword entry, returns NULL
 350    if not found, else returns pointer to rest of line
 351 */
 352
 353 static char *scan_header(const char *file, const char *name, int offset, int removespaces)
 354 {
 355     int i, j, iStringLength;
 356     static char temp[512];
 357     FILE *f;
 358     char *p, *q;
 359
 360     if ((f = fopen(file, "rt")) == NULL) return NULL;
 361     if (offset) fseek(f, offset, SEEK_SET);
 362
 363     for (i=0; i<200; )
 364     {
 365         if (fgets(temp, 500, f) == NULL      ) break;       /* end of file */
 366
 367         if (removespaces)
 368         {
 369             temp[500] = 0;
 370             p = q = temp;                       /* remove spaces */
 371             iStringLength = strlen(temp);
 372             for (j=0; j<iStringLength; j++)
 373                 if (*q!=' ' && *q!=8) *p++ = *q++;
 374                 else q++;
 375             *p++ = 0;
 376         }
 377
 378         if (temp[0] == 12                    ) break;       /* ^L end of header */
 379         if (temp[0] != '#') i++;                            /* The first 200 non comment lines must be read before data is opened. */
 380         if ((p = strchr(temp+1, '=')) == NULL) continue;    /* no '=' */
 381         if (memicmp(temp, name, p-temp)      ) continue;    /* no match */
 382
 383         p++;                                                /* match, skip = */
 384         if (p[strlen(p)-1] == '\n')                         /* remove \n */
 385             p[strlen(p)-1] = 0;
 386
 387         fclose (f);
 388         return p;
 389     }
 390
 391     fclose(f);
 392     return NULL;
 393 }
 394
 395 /* help routine, enumerates XDR file for names of keyword entrys, returns NULL
 396    if not found, else returns pointer to start of line upto '='
 397 */
 398 // //Commented out because it does not seem to be needed for vv
 399 //static char *enum_header(char *file, int iEntry, int offset)
 400 //{
 401 //    int i, count;
 402 //    static char temp[512];
 403 //    FILE *f;
 404 //    char *p;
 405 //
 406 //    if ((f = fopen(file, "rt")) == NULL) return NULL;
 407 //    if (offset) fseek(f, offset, SEEK_SET);
 408 //
 409 //    count = 0;
 410 //
 411 //    for (i=0; i<200; )
 412 //    {
 413 //        if (fgets(temp, 500, f) == NULL      ) break;       /* end of file */
 414 //
 415 //        if (temp[0] == 12                    ) break;       /* ^L end of header */
 416 //        if (temp[0] != '#') i++;                            /* The first 200 non comment lines must be read before data is opened. */
 417 //        if ((p = strchr(temp+1, '=')) == NULL) continue;    /* no '=' */
 418 //        if (count++ != iEntry)                 continue;    /* no match */
 419 //
 420 //        *p=0;                                               /* match, end at  = */
 421 //
 422 //        fclose (f);
 423 //        return temp;
 424 //    }
 425 //
 426 //    fclose(f);
 427 //    return NULL;
 428 //}
 429
 430
 431 static int get_nki_compressed_size(FILE *f)
 432 {
 433     NKI_MODE2           Header;
 434     int                 iMode;
 435
 436     fread((void *)&Header, sizeof(Header), 1 , f);
 437
 438     iMode = Header.iMode;
 439
 440     switch (iMode)
 441     {
 442     case  1:
 443     case  3:
 444         return 0;
 445     case  2:
 446     case  4:
 447         return Header.iCompressedSize + sizeof(Header);
 448     default:
 449         return 0;
 450     }
 451 }
 452
 453 /* decoder for NKI private compressed pixel data
 454    arguments: dest    = (in) points to area where destination data is written (short)
 455               src     = (in) points compressed source data (byte stream)
 456               size    = (in) number of bytes source data (safety)
 457
 458    The return value is the number of pixels that have been processed.
 459
 460    The compressed data looks like:
 461    (number of pixels)-1 times:
 462      OR 1 byte   = LL     7  bits signed (difference pixel[1] - pixel[0]);
 463      OR 2 bytes  = HHLL   15 bits signed (difference pixel[1] - pixel[0]) xored with 0x4000;
 464      OR 3 bytes  = 7FHHLL 16 bits absolute pixel data if 15 bits difference is exceeded
 465      OR 2 bytes  = 80NN   run length encode NN zero differences (max 255)
 466 for mode 3 and 4 added:
 467      OR 2 bytes  = CONN   encode NN 4 bit differences (max 255)
 468
 469    Performance on typical CT or MRI is >2x compression and a very good speed
 470
 471    This code is not valid on HIGHENDIAN (high byte first) machines
 472 */
 473
 474 static int global_len;
 475
 476 static int nki_private_decompress(short int *dest, signed char *src, int size)
 477 {
 478     int                 npixels, retvalue, mode, iMode, val, j;
 479     NKI_MODE2*          pHeader = (NKI_MODE2*)src;
 480     unsigned long               iCRC=0, iCRC2=0;
 481     //unsigned char*    pDestStart = (unsigned char*)dest;
 482     signed char           *save, *end;
 483
 484     retvalue = npixels = pHeader->iOrgSize;
 485     iMode = pHeader->iMode;             // safety: this value is checked in case statement
 486
 487     if (npixels<1) return 0;            // safety: check for invalid npixels value
 488
 489     /* Up till now only Mode=1, 2, 3, and 4 are supported */
 490
 491     switch (iMode)
 492     {
 493     case 1:
 494         save = src;
 495
 496         src += 8;                               // mode 1 only has 8 bytes header: iOrgSize and iMode
 497         end  = src + size - 3;          // for overflow check if we are close to end of input buffer
 498
 499         *dest = *(short int *)src;
 500         src += 2;
 501         npixels--;
 502
 503         do
 504         {
 505             if (src > end)                      // check whether the last few messages fit in input buffer
 506             {
 507                 if (src<end+3) val = *src;
 508                 else           val = 0;
 509
 510                 if (val >= -64 && val <= 63)      mode = 1;     // 7 bit difference
 511                 else if (val==0x7f)                 mode = 3;   // 16 bit value
 512                 else if ((val&0xff)==0x80)          mode = 2;   // run length encoding
 513                 else                                mode = 2;
 514
 515                 if (src+mode > end+3)
 516                     return 0;                   // safety: overflow input data
 517             }
 518
 519             val = *src;
 520
 521             if (val >= -64 && val <= 63)        // 7 bit difference
 522             {
 523                 dest[1] = dest[0] + val;
 524                 dest++;
 525                 src++;
 526             }
 527             else if (val==0x7f)         // 16 bit value
 528             {
 529                 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
 530                 dest++;
 531                 src+=3;
 532             }
 533             else if ((val&0xff)==0x80)  // run length encoding
 534             {
 535                 mode = ((unsigned char *)src)[1];
 536                 npixels -= mode-1;
 537                 if (npixels<=0) return 0;       // safety: overflow output data
 538                 do
 539                 {
 540                     dest[1] = dest[0];
 541                     dest++;
 542                 }
 543                 while (--mode);
 544                 src+=2;
 545             }
 546             else
 547             {
 548                 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
 549                 dest[1] = dest[0] + diff;       // 15 bit difference
 550                 dest++;
 551                 src+=2;
 552             }
 553         }
 554         while (--npixels);
 555
 556         global_len = src-save;
 557
 558         break;
 559
 560     case 2:
 561         src += sizeof(NKI_MODE2);
 562         save = src;
 563         end  = src + pHeader->iCompressedSize - 3;
 564
 565         if (end > src + size - 3)
 566             end = src + size - 3;               // may occur if pHeader is corrupted
 567
 568         *dest = val = *(short int *)src;
 569         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 570         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 571         src+=2;
 572
 573         npixels--;
 574
 575         do
 576         {
 577             if (src > end)                      // check whether the last few messages fit in input buffer
 578             {
 579                 if (src<end+3) val = *src;
 580                 else           val = 0;
 581
 582                 if (val >= -64 && val <= 63)      mode = 1;     // 7 bit difference
 583                 else if (val==0x7f)                 mode = 3;   // 16 bit value
 584                 else if ((val&0xff)==0x80)          mode = 2;   // run length encoding
 585                 else                                mode = 2;
 586
 587                 if (src+mode > end+3)
 588                     break;                      // safety: overflow input data
 589             }
 590
 591             val = *src;
 592
 593             if (val >= -64 && val <= 63)        // 7 bits difference
 594             {
 595                 dest[1] = val = dest[0] + val;
 596                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 597                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 598                 dest++;
 599                 src++;
 600             }
 601             else if (val==0x7f)         // 16 bit value
 602             {
 603                 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
 604
 605                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 606                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 607                 dest++;
 608                 src+=3;
 609             }
 610             else if ((val&0xff)==0x80)  // run length encoding
 611             {
 612                 mode = ((unsigned char *)src)[1];
 613                 npixels -= mode-1;
 614                 if (npixels<=0) break;  // safety: overflow output data
 615                 do
 616                 {
 617                     dest[1] = val = dest[0];
 618                     iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 619                     iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 620                     dest++;
 621                 }
 622                 while (--mode);
 623                 src+=2;
 624             }
 625             else
 626             {
 627                 signed short diff = ((val^0x40)<<8) + ((unsigned char *)src)[1];
 628                 dest[1] = val = dest[0] + diff; // 15 bit difference
 629                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 630                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 631                 dest++;
 632                 src+=2;
 633             }
 634         }
 635         while (--npixels);
 636
 637         if (iCRC2 != pHeader->iOrgCRC)  // if error in output CRC:
 638         {
 639             src = save;                 // check input CRC
 640             while (src < end)
 641             {
 642                 iCRC = CRC32_table[(unsigned char)iCRC ^ (unsigned char)src[0]] ^ ((iCRC >> 8));
 643                 src++;
 644             }
 645
 646             if (iCRC != pHeader->iCompressedCRC)
 647             {
 648                 AVSerror("XDR decompression: the file is corrupted");
 649                 retvalue=0;
 650             }
 651             else
 652             {
 653                 AVSerror("XDR decompression: internal error");
 654                 retvalue=0;
 655             }
 656         }
 657
 658         global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
 659
 660         break;
 661
 662     case 3:
 663         save = src;
 664
 665         src += 8;                               // mode 3 only has 8 bytes header: iOrgSize and iMode
 666         end  = src + size - 3;          // for overflow check if we are close to end of input buffer
 667
 668         *dest = *(short int *)src;
 669         src += 2;
 670         npixels--;
 671
 672         do
 673         {
 674             if (src > end)                      // check whether the last few messages fit in input buffer
 675             {
 676                 if (src<end+3) val = *src;
 677                 else           val = 0;
 678
 679                 if (val >= -63 && val <= 63)      mode = 1;     // 7 bit difference
 680                 else if (val==0x7f)                 mode = 3;   // 16 bit value
 681                 else if ((val&0xff)==0x80)          mode = 2;   // run length encoding
 682                 else if ((val&0xff)==0xC0)          mode = 2;   // 4 bit encoding
 683                 else                                mode = 2;
 684
 685                 if (src+mode > end+3)
 686                     return 0;                   // safety: overflow input data
 687             }
 688
 689             val = *src;
 690
 691             if (val >= -63 && val <= 63)        // 7 bit difference
 692             {
 693                 dest[1] = dest[0] + val;
 694                 dest++;
 695                 src++;
 696             }
 697             else if (val==0x7f)         // 16 bit value
 698             {
 699                 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
 700                 dest++;
 701                 src+=3;
 702             }
 703             else if ((val&0xff)==0x80)  // run length encoding
 704             {
 705                 mode = ((unsigned char *)src)[1];
 706                 npixels -= mode-1;
 707                 if (npixels<=0) return 0;       // safety: overflow output data
 708                 do
 709                 {
 710                     dest[1] = dest[0];
 711                     dest++;
 712                 }
 713                 while (--mode);
 714                 src+=2;
 715             }
 716             else if ((val&0xff)==0xC0)  // 4 bit run
 717             {
 718                 mode = ((unsigned char *)src)[1];
 719                 npixels -= mode-1;
 720                 mode/=2;
 721                 src+=2;
 722                 if (npixels<=0) return 0;       // safety: overflow output data
 723                 do
 724                 {
 725                     val = *src++;
 726                     dest[1] = dest[0] + (val>>4);
 727                     dest++;
 728                     if (val&8) val |= 0xfffffff0;
 729                     else val &= 0x0f;
 730                     dest[1] = dest[0] + val;
 731                     dest++;
 732                 }
 733                 while (--mode);
 734             }
 735             else
 736             {
 737                 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
 738                 dest[1] = dest[0] + diff;       // 15 bit difference
 739                 dest++;
 740                 src+=2;
 741             }
 742         }
 743         while (--npixels);
 744
 745         global_len = src-save;
 746
 747         break;
 748
 749     case 4:
 750         src += sizeof(NKI_MODE2);
 751         save = src;
 752         end  = src + pHeader->iCompressedSize - 3;
 753
 754         if (end > src + size - 3)
 755             end = src + size - 3;               // may occur if pHeader is corrupted
 756
 757         *dest = val = *(short int *)src;
 758         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 759         iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 760         src += 2;
 761         npixels--;
 762
 763         do
 764         {
 765             if (src > end)                      // check whether the last few messages fit in input buffer
 766             {
 767                 if (src<end+3) val = *src;
 768                 else           val = 0;
 769
 770                 if (val >= -63 && val <= 63)      mode = 1;     // 7 bit difference
 771                 else if (val==0x7f)                 mode = 3;   // 16 bit value
 772                 else if ((val&0xff)==0x80)          mode = 2;   // run length encoding
 773                 else if ((val&0xff)==0xC0)          mode = 2;   // 4 bit encoding
 774                 else                                mode = 2;
 775
 776                 if (src+mode > end+3)
 777                     return 0;                   // safety: overflow input data
 778             }
 779
 780             val = *src;
 781
 782             if (val >= -63 && val <= 63)        // 7 bit difference
 783             {
 784                 dest[1] = val = dest[0] + val;
 785                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 786                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 787                 dest++;
 788                 src++;
 789             }
 790             else if (val==0x7f)         // 16 bit value
 791             {
 792                 dest[1] = val = ((int)(((unsigned char *)src)[1])<<8) + ((unsigned char*)src)[2];
 793                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 794                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 795                 dest++;
 796                 src+=3;
 797             }
 798             else if ((val&0xff)==0x80)  // run length encoding
 799             {
 800                 mode = ((unsigned char *)src)[1];
 801                 npixels -= mode-1;
 802                 if (npixels<=0) return 0;       // safety: overflow output data
 803                 do
 804                 {
 805                     dest[1] = val = dest[0];
 806                     iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 807                     iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 808                     dest++;
 809                 }
 810                 while (--mode);
 811                 src+=2;
 812             }
 813             else if ((val&0xff)==0xC0)  // 4 bit run
 814             {
 815                 mode = ((unsigned char *)src)[1];
 816                 npixels -= mode-1;
 817                 mode/=2;
 818                 src+=2;
 819                 if (npixels<=0) return 0;       // safety: overflow output data
 820                 do
 821                 {
 822                     val = *src++;
 823                     dest[1] = j = dest[0] + (val>>4);
 824                     iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
 825                     iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
 826                     dest++;
 827                     if (val&8) val |= 0xfffffff0;
 828                     else val &= 0x0f;
 829                     dest[1] = j = dest[0] + val;
 830                     iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)j] ^ ((iCRC2 >> 8));
 831                     iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(j>>8)] ^ ((iCRC2 >> 8));
 832                     dest++;
 833                 }
 834                 while (--mode);
 835             }
 836             else
 837             {
 838                 signed short diff = ((val^0x40)<<8) + (unsigned char)(src[1]);
 839                 dest[1] = val = dest[0] + diff; // 15 bit difference
 840                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)val] ^ ((iCRC2 >> 8));
 841                 iCRC2 = CRC32_table[(unsigned char)iCRC2 ^ (unsigned char)(val>>8)] ^ ((iCRC2 >> 8));
 842                 dest++;
 843                 src+=2;
 844             }
 845         }
 846         while (--npixels);
 847
 848         if (iCRC2 != pHeader->iOrgCRC)  // if error in output CRC:
 849             retvalue=0;
 850
 851         global_len = sizeof(NKI_MODE2) + pHeader->iCompressedSize;
 852
 853         break;
 854
 855
 856     default:
 857         AVSerror("XDR decompression: unsupported mode");
 858         return 0;
 859     }
 860
 861     return retvalue;
 862 }
 863
 864
 865 //====================================================================
 866 // Read image information (copied from XDRreader)
 867 int nkitk::XDRImageIO::ReadImageInformationWithError()
 868 {
 869     int      offset=0;
 870     itk::Vector<int,MAXDIM> dim;
 871     int      veclen=1; //, data=AVS_TYPE_BYTE, field=UNIFORM;
 872 #if 0
 873     int      iNkiCompression = 0;
 874 #endif
 875     int      total=1/*, datasize=0, iNumRead, HeaderSize*/;
 876     unsigned int coords=0,i,j,ndim,nspace;
 877     char     temp[512];
 878     FILE     *fstream;
 879     char     *c;
 880 #if 0
 881     char     *buff;
 882     AVSfield FieldTemplate;
 883 #endif
 884     long     swap_test = 0x1000000;      /* first byte is 1 when low-endian */
 885     forcenoswap=0;
 886     char     *file = const_cast<char *>(m_FileName.c_str());
 887     AVSType  field=UNIFORM;
 888
 889
 890     fstream = fopen(file, "rt");
 891     if (fstream == NULL) return ER_XDR_OPEN;
 892
 893     fgets(temp, 500, fstream);
 894     fclose(fstream);
 895
 896     if (memcmp(temp, "# AVS field file (produced by avs_nfwrite.c)", 44)==0) forcenoswap=1;
 897
 898     c      = scan_header(file, "ndim", offset, 1);
 899     if (!c) return ER_XDR_NDIM;
 900
 901     ndim   = atoi(c);
 902     if (ndim<1 || ndim>MAXDIM) return ER_XDR_NDIM;
 903     SetNumberOfDimensions(ndim);
 904
 905     nspace = ndim;
 906
 907     for (i=0; i<ndim; i++)
 908     {
 909         sprintf(temp, "dim%d", i+1);
 910         c = scan_header(file, temp, offset, 1);
 911         if (!c) return ER_XDR_DIM;
 912         dim[i]=atoi(c);
 913         if (dim[i]<1) return ER_XDR_DIM;
 914
 915         total  *= dim[i];
 916         coords += dim[i];
 917     }
 918     for (i=0; i<ndim; i++) {
 919         SetDimensions(i,dim[i]);
 920         SetSpacing(i,1.);
 921         SetOrigin(i,0.);
 922     }
 923
 924     c = scan_header(file, "nspace", offset, 1);
 925     if (c) nspace = atoi(c);
 926     if (nspace<1 || ndim > MAXDIM) return ER_XDR_NSPACE;
 927     if (nspace != ndim) return ER_NOT_HANDLED;
 928
 929     c = scan_header(file, "veclen", offset, 1);
 930     if (c) veclen = atoi(c);
 931     if (veclen<0 /*|| veclen>1000*/) return ER_XDR_VECLEN;
 932     SetNumberOfComponents(veclen);
 933     if (veclen==1) SetPixelType(itk::ImageIOBase::SCALAR);
 934     else          SetPixelType(itk::ImageIOBase::VECTOR);
 935
 936     c = scan_header(file, "data", offset, 1);
 937     if (c)
 938     {
 939         if (memicmp(c, "byte",  4) == 0 || memicmp(c, "xdr_byte",  8) == 0) SetComponentType(itk::ImageIOBase::CHAR);
 940         else if (memicmp(c, "short", 5) == 0 || memicmp(c, "xdr_short", 9) == 0) SetComponentType(itk::ImageIOBase::SHORT);
 941         else if (memicmp(c, "int" ,  3) == 0 || memicmp(c, "xdr_int" ,  7) == 0) SetComponentType(itk::ImageIOBase::INT);
 942         else if (memicmp(c, "real",  4) == 0 || memicmp(c, "xdr_real",  8) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
 943         else if (memicmp(c, "float", 5) == 0 || memicmp(c, "xdr_float", 9) == 0) SetComponentType(itk::ImageIOBase::FLOAT);
 944         else if (memicmp(c, "double",6) == 0 || memicmp(c, "xdr_double",10)== 0) SetComponentType(itk::ImageIOBase::DOUBLE);
 945         else return ER_XDR_DATA;
 946
 947         if (memicmp(c, "xdr_",  4) == 0) forcenoswap=0;
 948     }
 949
 950     //Read coords here
 951     c = scan_header(file, "field", offset, 1);
 952     if (c)
 953     {
 954         if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace *2;
 955         else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
 956         else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
 957         else return ER_XDR_FIELD;
 958     }
 959     else
 960         coords=0;
 961
 962     if (coords)                        /* expect AVS coordinates ? */
 963     {
 964         coords *= sizeof(float);
 965         fstream = fopen(m_FileName.c_str(), "rb");
 966         if (fstream == NULL) return ER_XDR_OPEN;
 967
 968         float *points = (float *)malloc(coords);
 969         if (points == NULL) return ER_OUTOFMEMORY;
 970
 971         //Seek to coordinates position in file
 972         if (fseek(fstream,-static_cast<int>(coords),SEEK_END)) return ER_XDR_READ;
 973         if (fread( /*(*output)->*/points, 1, coords, fstream ) == coords)
 974         { /* swap data if read-ok and required (xdr is low-endian) */
 975             if (!(*(char *)(&swap_test)) && !forcenoswap)
 976             {
 977                 c = (char *)/*(*output)->*/points;
 978                 for (i=0; i<coords; i+=4)
 979                 {
 980                     j = c[i];
 981                     c[i]   = c[i+3];
 982                     c[i+3] = j;
 983                     j = c[i+1];
 984                     c[i+1] = c[i+2];
 985                     c[i+2] = j;
 986                 }
 987             }
 988         }
 989
 990         switch (field) {
 991         case UNIFORM:
 992             for (i=0; i<GetNumberOfDimensions(); i++) {
 993                 SetSpacing(i,10.*(points[i*2+1]-points[i*2])/(GetDimensions(i)-1));
 994                 SetOrigin(i,10.*points[i*2]);
 995             }
 996             break;
 997         case RECTILINEAR:
 998             //Rectilinear is reinterpreted as uniform because ITK does not know rectilinear
 999             //Error if fails
1000             for (i=0; i<GetNumberOfDimensions(); i++) {
1001                 //Compute mean spacing
1002                 SetSpacing(i,10*(points[GetDimensions(i)-1]-points[0])/(GetDimensions(i)-1));
1003                 SetOrigin(i,10*points[0]);
1004
1005                 //Test if rectilinear image is actually uniform (tolerance 0.1 mm)
1006                 for (j=0; j<GetDimensions(i)-1; j++) {
1007                     if (fabs((points[j+1]-points[j])*10-GetSpacing(i))>0.1) {
1008                         free(points);
1009                         fclose(fstream);
1010                         return ER_NOT_HANDLED;
1011                     }
1012                 }
1013                 points += (int)GetDimensions(i);
1014             }
1015             for (i=0; i<GetNumberOfDimensions(); i++)
1016                 points -= GetDimensions(i);
1017             break;
1018         case IRREGULAR:
1019             free(points);
1020             fclose(fstream);
1021             return ER_NOT_HANDLED;
1022         }
1023         free(points);
1024         fclose(fstream);
1025     }
1026     return OK;
1027 }
1028
1029 //====================================================================
1030 // Read image information (copied from XDRreader)
1031 void nkitk::XDRImageIO::ReadImageInformation() {
1032     int result = ReadImageInformationWithError();
1033     if (result) ITKError("nkitk::XDRImageIO::ReadImageInformation",result);
1034 }
1035
1036
1037 //====================================================================
1038 // Read Image Content (copied from XDRreader)
1039 int nkitk::XDRImageIO::ReadWithError(void * buffer)
1040 { //AVSINT   dim[5];
1041     int      /*ndim,*/ nspace/*, veclen=1, data=AVS_TYPE_BYTE, field=UNIFORM*/;
1042     int      iNkiCompression = 0;
1043     int      j, coords=0,  datasize=0, HeaderSize;
1044     unsigned int i,iNumRead,total=1;
1045     char     temp[512];
1046     FILE     *fstream;
1047     char     *c;
1048     char     *buff;
1049     //AVSfield FieldTemplate;
1050     long     swap_test = 0x1000000;      /* first byte is 1 when low-endian */
1051     //int      forcenoswap=0;
1052     char     *file = const_cast<char *>(m_FileName.c_str());
1053     int      offset=0;
1054     AVSType  field=UNIFORM;
1055
1056     for (i=0; i<GetNumberOfDimensions(); i++) coords += GetDimensions(i);
1057
1058     total = GetImageSizeInPixels();
1059     nspace = GetNumberOfDimensions();
1060
1061     c = scan_header(file, "field", offset, 1);
1062     if (c)
1063     {
1064         if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace*2;
1065         else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
1066         else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
1067         else return ER_XDR_FIELD;
1068     }
1069     else
1070         coords=0;
1071
1072     c = scan_header(file, "nki_compression", offset, 1);
1073     if (c) iNkiCompression = atoi(c);
1074
1075     c = scan_header(file, "coord1[0]", offset, 1);
1076     if (c) HeaderSize = 32768;
1077     else HeaderSize = 2048;
1078
1079 #if 0
1080     FIELDdefault(&FieldTemplate);
1081     FieldTemplate.ndim    = ndim;
1082     FieldTemplate.nspace  = nspace;
1083     FieldTemplate.veclen  = veclen;
1084     FieldTemplate.type    = data;
1085     FieldTemplate.uniform = field;
1086     FieldTemplate.size    = datasize;
1087 #endif
1088
1089     fstream = fopen(file, "rb");
1090     if (fstream == NULL)
1091         return ER_XDR_OPEN;
1092
1093     if (offset) fseek(fstream, offset, SEEK_SET);
1094
1095     while (1)
1096     {
1097         if (fgets(temp, 500, fstream) == NULL )
1098             return ER_XDR_NOCTRLL; /* end of file */
1099
1100         if (temp[0] == 10) continue;
1101
1102         if (temp[0] == 12)
1103         {
1104             fseek(fstream, -2, SEEK_CUR);
1105             break;
1106         } /* ^L end of header */
1107
1108         if (temp[0] != '#') break;
1109     }
1110
1111     buff = (char*)malloc(HeaderSize);
1112     if (buff == NULL)
1113     {
1114         return ER_OUTOFMEMORY;
1115     }
1116     memset(buff, 0, HeaderSize);
1117     iNumRead = fread(buff, 1, HeaderSize, fstream);
1118     if (iNumRead < 1)
1119     {
1120         free(buff);
1121         fclose(fstream);
1122         return ER_XDR_READ;
1123     }
1124
1125     for (i=0; i<iNumRead; i++) {
1126         if (buff[i] == 12) break;
1127     }
1128
1129     free(buff);
1130
1131     if (i==iNumRead) return ER_XDR_NOCTRLL;
1132
1133 #if 0
1134     if (*output) AVSfield_free(*output);
1135     *output = (AVSfield *) AVSfield_alloc(&FieldTemplate, dim);
1136     if (*output == NULL)
1137     {
1138         fclose(fstream);
1139         return ER_OUTOFMEMORY;
1140     }
1141
1142     total  *= datasize * veclen;
1143     coords *= sizeof(float);
1144 #endif
1145     total = GetImageSizeInBytes();
1146
1147     //We add casts because the resulting quantity can be negative.
1148     //There is no risk of looping because i and iNumRead are about the size of the header
1149     fseek(fstream, static_cast<int>(i)+2-static_cast<int>(iNumRead), SEEK_CUR);
1150
1151     if (total && iNkiCompression)
1152     {
1153         long            iCurPos;
1154         unsigned long iSize;
1155         signed char*    pCompressed;
1156
1157         /* Read or guess the size of the compressed data */
1158         iCurPos = ftell(fstream);
1159         iSize = get_nki_compressed_size(fstream);
1160
1161         if (iSize==0)
1162         {
1163             fseek(fstream, 0, SEEK_END);
1164             iSize = ftell(fstream);
1165             iSize = iSize - iCurPos - coords;
1166
1167             // Get compressed size from header if possible; else use uncompressed size as safe estimate
1168             if (iSize>total && offset) iSize=total+8;
1169         }
1170
1171         fseek(fstream, iCurPos, SEEK_SET);
1172
1173         /* Allocate space for the compressed pixels */
1174         pCompressed = (signed char*)malloc(iSize);
1175         if (!pCompressed)
1176         {
1177             fclose(fstream);
1178 #if 0
1179             if (*output) AVSfield_free(*output);
1180             *output = NULL;
1181 #endif
1182             return ER_OUTOFMEMORY;
1183         }
1184
1185         /* Read the compressed pixels */
1186         if (fread( (void *)pCompressed, 1, iSize, fstream ) != iSize)
1187         {
1188             fclose(fstream);
1189 #if 0
1190             if (*output) AVSfield_free(*output);
1191             *output = NULL;
1192 #endif
1193             return ER_XDR_READ;
1194         }
1195
1196         if (!nki_private_decompress((short*)buffer, pCompressed, iSize))
1197         {
1198             fclose(fstream);
1199 #if 0
1200             if (*output) AVSfield_free(*output);
1201             *output = NULL;
1202 #endif
1203             return ER_DECOMPRESSION;
1204         }
1205
1206         // if (offset)
1207         fseek(fstream, iCurPos + global_len, SEEK_SET);
1208
1209         free(pCompressed);
1210         goto READ_COORDS;
1211     }
1212
1213
1214     if (total)
1215     {
1216         if (fread( (void *)buffer, 1, total, fstream ) != total)
1217         {
1218             fclose(fstream);
1219 #if 0
1220             if (*output) AVSfield_free(*output);
1221             *output = NULL;
1222 #endif
1223             return ER_XDR_READ;
1224         }
1225     }
1226
1227     /* swap data if required (xdr is low-endian) */
1228
1229     datasize = GetComponentSize();
1230     if (!(*(char *)(&swap_test)) && !forcenoswap)
1231     {
1232         if (datasize==2)
1233         {
1234             c = (char *)buffer;
1235             for (i=0; i<total; i+=2)
1236             {
1237                 j = c[i];
1238                 c[i]   = c[i+1];
1239                 c[i+1] = j;
1240             }
1241         }
1242         else if (datasize==4)
1243         {
1244             c = (char *)buffer;
1245             for (i=0; i<total; i+=4)
1246             {
1247                 j = c[i];
1248                 c[i]   = c[i+3];
1249                 c[i+3] = j;
1250                 j = c[i+1];
1251                 c[i+1] = c[i+2];
1252                 c[i+2] = j;
1253             }
1254         }
1255         else if (datasize==8)
1256         {
1257             c = (char *)buffer;
1258             for (i=0; i<total; i+=8)
1259             {
1260                 j = c[i];
1261                 c[i]   = c[i+7];
1262                 c[i+7] = j;
1263                 j = c[i+1];
1264                 c[i+1] = c[i+6];
1265                 c[i+6] = j;
1266                 j = c[i+2];
1267                 c[i+2] = c[i+5];
1268                 c[i+5] = j;
1269                 j = c[i+3];
1270                 c[i+3] = c[i+4];
1271                 c[i+4] = j;
1272             }
1273         }
1274     }
1275
1276 READ_COORDS:
1277 #if 0
1278     if (coords)                        /* expect AVS coordinates ? */
1279     {
1280         if (fread( (*output)->points, 1, coords, fstream ) == coords)
1281         { /* swap data if read-ok and required (xdr is low-endian) */
1282
1283             if (!(*(char *)(&swap_test)) && !forcenoswap)
1284             {
1285                 c = (char *)(*output)->points;
1286                 for (i=0; i<coords; i+=4)
1287                 {
1288                     j = c[i];
1289                     c[i]   = c[i+3];
1290                     c[i+3] = j;
1291                     j = c[i+1];
1292                     c[i+1] = c[i+2];
1293                     c[i+2] = j;
1294                 }
1295             }
1296         }
1297     }
1298 #endif
1299     fclose(fstream);
1300     return OK;
1301 }
1302
1303 //====================================================================
1304 // Read image information (copied from XDRreader)
1305 void nkitk::XDRImageIO::Read(void * buffer) {
1306     int result = ReadWithError(buffer);
1307     if (result) ITKError("nkitk::XDRImageIO::Read",result);
1308 }
1309
1310 #if 0
1311 static int XDRreader_preview(AVSfield **output, char *file, int offset)
1312 {
1313     AVSINT   dim[MAXDIM], dim2[MAXDIM], ds[MAXDIM];
1314     int      ndim, nspace, veclen=1, data=AVS_TYPE_BYTE, field=UNIFORM;
1315     int      iNkiCompression = 0;
1316     int      total=1, i, j, k, l, m, coords=0,  datasize=0, iNumRead, len, start, sliceax;
1317     char     temp[512];
1318     FILE     *fstream;
1319     char     *c;
1320     char     *buff;
1321     AVSfield FieldTemplate;
1322     long     swap_test = 0x1000000;      /* first byte is 1 when low-endian */
1323     int      forcenoswap=0;
1324
1325     fstream = fopen(file, "rt");
1326     if (fstream == NULL) return ER_XDR_OPEN;
1327     fgets(temp, 500, fstream);
1328     fclose(fstream);
1329
1330     if (memcmp(temp, "# AVS field file (produced by avs_nfwrite.c)", 44)==0) forcenoswap=1;
1331
1332     c      = scan_header(file, "ndim", offset, 1);
1333     if (!c) return ER_XDR_NDIM;
1334     ndim   = atoi(c);
1335     if (ndim<1 || ndim>MAXDIM) return ER_XDR_NDIM;
1336     nspace = ndim;
1337
1338     /* defaults for dimensions and downsize */
1339
1340     for (i=0; i<MAXDIM; i++)
1341         dim[i] = dim2[i] = ds[i] = 1;
1342
1343     for (i=0; i<ndim; i++)
1344     {
1345         sprintf(temp, "dim%d", i+1);
1346         c = scan_header(file, temp, offset, 1);
1347         if (!c) return ER_XDR_DIM;
1348         dim[i]=atoi(c);
1349         if (dim[i]<1 || dim[i]>200000L) return ER_XDR_DIM;
1350
1351         total  *= dim[i];
1352         coords += dim[i];
1353     }
1354
1355     c = scan_header(file, "nspace", offset, 1);
1356     if (c) nspace = atoi(c);
1357     if (nspace<1 || ndim > MAXDIM) return ER_XDR_NSPACE;
1358
1359     c = scan_header(file, "veclen", offset, 1);
1360     if (c) veclen = atoi(c);
1361     if (veclen<0 || veclen>100) return ER_XDR_VECLEN;
1362
1363     c = scan_header(file, "data", offset, 1);
1364
1365     if (c)
1366     {
1367         if (memicmp(c, "byte",  4) == 0)     data=AVS_TYPE_BYTE,   datasize=1;
1368         else if (memicmp(c, "short", 5) == 0)     data=AVS_TYPE_SHORT,  datasize=2;
1369         else if (memicmp(c, "int" ,  3) == 0)     data=AVS_TYPE_INTEGER,datasize=4;
1370         else if (memicmp(c, "real",  4) == 0)     data=AVS_TYPE_REAL,   datasize=4;
1371         else if (memicmp(c, "float", 5) == 0)     data=AVS_TYPE_REAL,   datasize=4;
1372         else if (memicmp(c, "double",6) == 0)     data=AVS_TYPE_DOUBLE, datasize=8;
1373
1374         else if (memicmp(c, "xdr_byte",  8) == 0) data=AVS_TYPE_BYTE,   datasize=1, forcenoswap=0;
1375         else if (memicmp(c, "xdr_short", 9) == 0) data=AVS_TYPE_SHORT,  datasize=2, forcenoswap=0;
1376         else if (memicmp(c, "xdr_int" ,  7) == 0) data=AVS_TYPE_INTEGER,datasize=4, forcenoswap=0;
1377         else if (memicmp(c, "xdr_real",  8) == 0) data=AVS_TYPE_REAL,   datasize=4, forcenoswap=0;
1378         else if (memicmp(c, "xdr_float", 9) == 0) data=AVS_TYPE_REAL,   datasize=4, forcenoswap=0;
1379         else if (memicmp(c, "xdr_double",10)== 0) data=AVS_TYPE_DOUBLE, datasize=8, forcenoswap=0;
1380         else                                      return ER_XDR_DATA;
1381     }
1382
1383     c = scan_header(file, "field", offset, 1);
1384     if (c)
1385     {
1386         if (memicmp(c, "unifo", 5) == 0) field=UNIFORM, coords=nspace*2;
1387         else if (memicmp(c, "recti", 5) == 0) field=RECTILINEAR;
1388         else if (memicmp(c, "irreg", 5) == 0) field=IRREGULAR, coords=total*nspace;
1389         else return ER_XDR_FIELD;
1390     }
1391     else
1392         coords=0;
1393
1394     c = scan_header(file, "nki_compression", offset, 1);
1395     if (c) iNkiCompression = atoi(c);
1396     if (iNkiCompression)
1397     {
1398         fclose(fstream);
1399         return ER_ILLCOMMFUNCT;
1400     }
1401
1402     FIELDdefault(&FieldTemplate);
1403     FieldTemplate.ndim    = ndim;
1404     FieldTemplate.nspace  = nspace;
1405     FieldTemplate.veclen  = veclen;
1406     FieldTemplate.type    = data;
1407     FieldTemplate.uniform = field;
1408     FieldTemplate.size    = datasize;
1409
1410     fstream = fopen(file, "rb");
1411
1412     if (fstream == NULL)
1413         return ER_XDR_OPEN;
1414
1415     buff = (char *)malloc(8192);
1416     if (buff == NULL)
1417     {
1418         return ER_OUTOFMEMORY;
1419     }
1420     memset(buff, 0, 8192);
1421     fseek(fstream, offset, SEEK_SET);
1422
1423     while (1)
1424     {
1425         if (fgets(temp, 500, fstream) == NULL )
1426             return ER_XDR_NOCTRLL; /* end of file */
1427
1428         if (temp[0] == 10) continue;
1429
1430         if (temp[0] == 12)
1431         {
1432             fseek(fstream, -2, SEEK_CUR);
1433             break;
1434         } /* ^L end of header */
1435
1436         if (temp[0] != '#') break;
1437     }
1438     start = ftell(fstream);
1439
1440     iNumRead = fread(buff, 1, 8192, fstream);
1441     if (iNumRead < 1)
1442     {
1443         free(buff);
1444         fclose(fstream);
1445         return ER_XDR_READ;
1446     }
1447
1448     for (i=0; i<iNumRead; i++) {
1449         if (buff[i] == 12) break;
1450     }
1451
1452     free(buff);
1453
1454     if (i==iNumRead) return ER_XDR_NOCTRLL;
1455
1456     start += i+2;
1457
1458     /* determine slice axis and downsize slice axis to 3, others to 32 */
1459
1460     sliceax = 2;
1461
1462     if (ndim>=3)
1463     {
1464         if (dim[0]==dim[1]) sliceax = 2;
1465         else                sliceax = 1;
1466     }
1467
1468     total = 1;
1469
1470     for (i=0; i<ndim; i++)
1471     {
1472         if (i==sliceax) ds[i] = dim[i] / 3;
1473         else            ds[i] = dim[i] / 32;
1474         if (ds[i]==0) ds[i]=1;
1475
1476         dim2[i] = dim[i]/ds[i];
1477         if (dim2[i]==0) dim2[i] = 1;
1478
1479         total *= dim2[i];
1480     }
1481
1482     if (*output) AVSfield_free(*output);
1483     *output = (AVSfield *) AVSfield_alloc(&FieldTemplate, dim2);
1484     if (*output == NULL)
1485     {
1486         fclose(fstream);
1487         return ER_OUTOFMEMORY;
1488     }
1489
1490     total  *= datasize * veclen;
1491     coords *= sizeof(float);
1492
1493     /* Read and downsize the data */
1494
1495     if (total)
1496     {
1497         c = (char *)((*output)->field_data);
1498
1499         buff = (char *)malloc(dim[0] * datasize * veclen);
1500
1501         for (i=0; i<dim2[4]; i++)
1502             for (j=0; j<dim2[3]; j++)
1503                 for (k=0; k<dim2[2]; k++)
1504                     for (l=0; l<dim2[1]; l++)
1505                     {
1506                         len = i*ds[4];
1507                         len = len * dim[3] + j*ds[3];
1508                         len = len * dim[2] + k*ds[2];
1509                         len = len * dim[1] + l*ds[1];
1510                         len = len * dim[0] * datasize * veclen;
1511
1512                         // This is the time-critical statement
1513                         fseek(fstream, start + len, SEEK_SET);
1514
1515                         fread(buff, 1, dim[0]*datasize*veclen, fstream);
1516                         for (m=0; m<dim2[0]; m++)
1517                         {
1518                             memcpy(c, buff+m*datasize*veclen*ds[0], datasize*veclen);
1519                             c += datasize*veclen;
1520                         };
1521                     };
1522
1523         free (buff);
1524     }
1525
1526     /* swap data if required (xdr is low-endian) */
1527
1528     if (!(*(char *)(&swap_test))  && !forcenoswap)
1529     {
1530         if (datasize==2)
1531         {
1532             c = (char *)(*output)->field_data;
1533             for (i=0; i<total; i+=2)
1534             {
1535                 j = c[i];
1536                 c[i]   = c[i+1];
1537                 c[i+1] = j;
1538             }
1539         }
1540         else if (datasize==4)
1541         {
1542             c = (char *)(*output)->field_data;
1543             for (i=0; i<total; i+=4)
1544             {
1545                 j = c[i];
1546                 c[i]   = c[i+3];
1547                 c[i+3] = j;
1548                 j = c[i+1];
1549                 c[i+1] = c[i+2];
1550                 c[i+2] = j;
1551             }
1552         }
1553         else if (datasize==8)
1554         {
1555             c = (char *)(*output)->field_data;
1556             for (i=0; i<total; i+=8)
1557             {
1558                 j = c[i];
1559                 c[i]   = c[i+7];
1560                 c[i+7] = j;
1561                 j = c[i+1];
1562                 c[i+1] = c[i+6];
1563                 c[i+6] = j;
1564                 j = c[i+2];
1565                 c[i+2] = c[i+5];
1566                 c[i+5] = j;
1567                 j = c[i+3];
1568                 c[i+3] = c[i+4];
1569                 c[i+4] = j;
1570             }
1571         }
1572     }
1573
1574
1575     /*
1576       if (coords)
1577       { if (read( fHandle,(*output)->points, coords ) == coords)
1578         { if (!(*(char *)(&swap_test))  && !forcenoswap)
1579           { c = (char *)(*output)->points;
1580         for (i=0; i<coords; i+=4)
1581         { j = c[i];   c[i]   = c[i+3]; c[i+3] = j;
1582           j = c[i+1]; c[i+1] = c[i+2]; c[i+2] = j;
1583         }
1584           }
1585         }
1586       }
1587     */
1588
1589     fclose(fstream);
1590     return OK;
1591 }
1592
1593 int Rxdr_compute(AVSfield** ppOut, char* pszFileName, int iOffset)
1594 {
1595     int   rc;
1596     char  szMsg[64];
1597
1598     rc = XDRreader(ppOut, pszFileName, iOffset);
1599     if (rc == OK)
1600         rc = AVS_OK;
1601     else
1602     {
1603         strcpy(szMsg, "Avs_rxdr: ");
1604         strcat(szMsg, gl_ErrorMsg[rc]);
1605         AVSerror(szMsg);
1606         rc = AVS_ERROR;
1607     }
1608     return rc;
1609 }
1610
1611 int RxdrPreview_compute(AVSfield** ppOut, char* pszFileName, int iOffset)
1612 {
1613     int   rc;
1614     char  szMsg[64];
1615
1616     rc = XDRreader_preview(ppOut, pszFileName, iOffset);
1617     if (rc == OK)
1618         rc = AVS_OK;
1619     else
1620     {
1621         strcpy(szMsg, "Avs_rxdr: ");
1622         strcat(szMsg, gl_ErrorMsg[rc]);
1623         AVSerror(szMsg);
1624         rc = AVS_ERROR;
1625     }
1626     return rc;
1627 }
1628
1629 int RxdrHeader_compute(char *pszOut, char *pszFileName, char *pszEntry)
1630 {
1631     char *cc;
1632
1633     cc = scan_header(pszFileName, pszEntry, 0, 0);
1634     if (cc) strcpy(pszOut, cc);
1635     else    *pszOut = 0;
1636
1637     return AVS_OK;
1638 }
1639
1640 int RxdrEnum_compute(char *pszOut, char *pszFileName, int iEntry)
1641 {
1642     char *cc;
1643
1644     cc = enum_header(pszFileName, iEntry, 0);
1645     if (cc) strcpy(pszOut, cc);
1646     else    *pszOut = 0;
1647
1648     return AVS_OK;
1649 }
1650 #endif
1651
1652 //====================================================================
1653 // Read Image Information
1654 bool nkitk::XDRImageIO::CanReadFile(const char* FileNameToRead)
1655 {
1656     char     temp[512];
1657     FILE     *fstream;
1658
1659     fstream = fopen(FileNameToRead, "rt");
1660     if (fstream == NULL)
1661         return false;
1662     fgets(temp, 500, fstream);
1663     fclose(fstream);
1664
1665     if (memcmp(temp, "# AVS", 5)==0)
1666         return true;
1667     else
1668         return false;
1669 } ////
1670
1671 void nkitk::XDRImageIO::ITKError(std::string funcName, int msgID) {
1672     itkExceptionMacro(<< "Error in " << funcName << ". Message: " << gl_ErrorMsg[msgID]);
1673 }
1674
1675 #endif /* end #define NKITKXDRIMAGEIO_CXX */
1676