common/clitkXdrImageIOReader.cxx

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