common/rtkEdfImageIO.cxx

   1 /*=========================================================================
   2  *
   3  *  Copyright RTK Consortium
   4  *
   5  *  Licensed under the Apache License, Version 2.0 (the "License");
   6  *  you may not use this file except in compliance with the License.
   7  *  You may obtain a copy of the License at
   8  *
   9  *         http://www.apache.org/licenses/LICENSE-2.0.txt
  10  *
  11  *  Unless required by applicable law or agreed to in writing, software
  12  *  distributed under the License is distributed on an "AS IS" BASIS,
  13  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14  *  See the License for the specific language governing permissions and
  15  *  limitations under the License.
  16  *
  17  *=========================================================================*/
  18
  19 #include "rtkEdfImageIO.h"
  20
  21 // itk include (for itkReadRawBytesAfterSwappingMacro)
  22 #include <itkRawImageIO.h>
  23 #include <itk_zlib.h>
  24
  25 //--------------------------------------------------------------------
  26 /* Find value_ptr as pointer to the parameter of the given key in the header.
  27  * Returns NULL on success.
  28  */
  29 char*
  30 rtk::EdfImageIO::edf_findInHeader( char* header, const char* key )
  31 {
  32   char *value_ptr = strstr( header, key );
  33
  34   if (!value_ptr) return NULL;
  35   /* an edf line is "key     = value ;" */
  36   value_ptr = 1 + strchr( value_ptr + strlen(key), '=' );
  37   while (isspace(*value_ptr) ) value_ptr++;
  38   return value_ptr;
  39 }
  40
  41 //--------------------------------------------------------------------
  42 // Read Image Information
  43 void rtk::EdfImageIO::ReadImageInformation()
  44 {
  45   int    k;
  46   char * header = NULL;
  47   int    header_size = 0;
  48   char * p;
  49   gzFile inp;
  50
  51   inp = gzopen(m_FileName.c_str(), "rb");
  52   if (!inp)
  53     itkGenericExceptionMacro(<< "Cannot open input file " << m_FileName);
  54
  55   // read header: it is a multiple of 512 B ending by "}\n"
  56   while (header_size == 0 || strncmp(&header[header_size-2],"}\n",2) ) {
  57     int header_size_prev = header_size;
  58     header_size += 512;
  59     if (!header)
  60       header = (char*)malloc(header_size+1);
  61     else
  62       header = (char*)realloc(header, header_size+1);
  63     header[header_size_prev] = 0; /* protection against empty file */
  64     // fread(header+header_size_prev, 512, 1, fp);
  65     k = gzread(inp, header+header_size_prev, 512);
  66     if (k < 512) { /* protection against infinite loop */
  67       gzclose(inp);
  68       free(header);
  69       itkGenericExceptionMacro(<< "Damaged EDF header of "
  70                                << m_FileName
  71                                << ": not multiple of 512 B.");
  72       }
  73     header[header_size] = 0; /* end of string: protection against strstr later
  74                                on */
  75     }
  76
  77   // parse the header
  78   int   dim1 = -1, dim2 = -1, datalen = -1;
  79   char *otherfile_name = 0; // this file, or another file with the data (EDF vs
  80                             // EHF formats)
  81   int   otherfile_skip = 0;
  82
  83   if ( (p = edf_findInHeader(header, "EDF_BinaryFileName") ) ) {
  84     int plen = strcspn(p, " ;\n");
  85     otherfile_name = (char*)realloc(otherfile_name, plen+1);
  86     strncpy(otherfile_name, p, plen);
  87     otherfile_name[plen] = '\0';
  88     if ( (p = edf_findInHeader(header, "EDF_BinaryFilePosition") ) )
  89       otherfile_skip = atoi(p);
  90     }
  91
  92   if ( (p = edf_findInHeader(header, "Dim_1") ) )
  93     dim1 = atoi(p);
  94   if ( (p = edf_findInHeader(header, "Dim_2") ) )
  95     dim2 = atoi(p);
  96
  97 //  int orig1 = -1, orig2 = -1;
  98 //  if ((p = edf_findInHeader(header, "row_beg")))
  99 //    orig1 = atoi(p);
 100 //  if ((p = edf_findInHeader(header, "col_beg")))
 101 //    orig2 = atoi(p);
 102
 103   static const struct table3 edf_datatype_table[] =
 104     {
 105           { "UnsignedByte",    U_CHAR_DATATYPE,  1 },
 106           { "SignedByte",      CHAR_DATATYPE,    1 },
 107           { "UnsignedShort",   U_SHORT_DATATYPE, 2 },
 108           { "SignedShort",     SHORT_DATATYPE,   2 },
 109           { "UnsignedInteger", U_INT_DATATYPE,   4 },
 110           { "SignedInteger",   INT_DATATYPE,     4 },
 111           { "UnsignedLong",    U_L_INT_DATATYPE, 4 },
 112           { "SignedLong",      L_INT_DATATYPE,   4 },
 113           { "FloatValue",      FLOAT_DATATYPE,   4 },
 114           { "DoubleValue",     DOUBLE_DATATYPE,  8 },
 115           { "Float",           FLOAT_DATATYPE,   4 }, // Float and FloatValue
 116                                                       // are synonyms
 117           { "Double",          DOUBLE_DATATYPE,  8 }, // Double and DoubleValue
 118                                                       // are synonyms
 119           { NULL, -1, -1 }
 120     };
 121   if ( (p = edf_findInHeader(header, "DataType") ) ) {
 122     k = lookup_table3_nth(edf_datatype_table, p);
 123     if (k < 0) { // unknown EDF DataType
 124       gzclose(inp);
 125       free(header);
 126       itkGenericExceptionMacro( <<"Unknown EDF datatype \""
 127                                 << p
 128                                 << "\"");
 129       }
 130     datalen = edf_datatype_table[k].sajzof;
 131     switch(edf_datatype_table[k].value) {
 132       case U_CHAR_DATATYPE:
 133         SetComponentType(itk::ImageIOBase::UCHAR);
 134         break;
 135       case CHAR_DATATYPE:
 136         SetComponentType(itk::ImageIOBase::CHAR);
 137         break;
 138       case U_SHORT_DATATYPE:
 139         SetComponentType(itk::ImageIOBase::USHORT);
 140         break;
 141       case SHORT_DATATYPE:
 142         SetComponentType(itk::ImageIOBase::SHORT);
 143         break;
 144       case U_INT_DATATYPE:
 145         SetComponentType(itk::ImageIOBase::UINT);
 146         break;
 147       case INT_DATATYPE:
 148         SetComponentType(itk::ImageIOBase::INT);
 149         break;
 150       case U_L_INT_DATATYPE:
 151         SetComponentType(itk::ImageIOBase::UINT);
 152         break;
 153       case L_INT_DATATYPE:
 154         SetComponentType(itk::ImageIOBase::INT);
 155         break;
 156       case FLOAT_DATATYPE:
 157         SetComponentType(itk::ImageIOBase::FLOAT);
 158         break;
 159       case DOUBLE_DATATYPE:
 160         SetComponentType(itk::ImageIOBase::DOUBLE);
 161         break;
 162       }
 163     }
 164
 165   static const struct table edf_byteorder_table[] =
 166     {
 167           { "LowByteFirst",  LittleEndian }, /* little endian */
 168           { "HighByteFirst", BigEndian },    /* big endian */
 169           { NULL, -1 }
 170     };
 171
 172   int byteorder = LittleEndian;
 173   if ( (p = edf_findInHeader(header, "ByteOrder") ) ) {
 174     k = lookup_table_nth(edf_byteorder_table, p);
 175     if (k >= 0) {
 176
 177       byteorder = edf_byteorder_table[k].value;
 178       if(byteorder==LittleEndian)
 179         this->SetByteOrder(LittleEndian);
 180       else
 181         this->SetByteOrder(BigEndian);
 182       }
 183     } else
 184     itkWarningMacro(<<"ByteOrder not specified in the header! Not swapping bytes (figure may not be correct).");
 185   // Get and verify size of the data:
 186   int datasize = dim1 * dim2 * datalen;
 187   if ( (p = edf_findInHeader(header, "Size") ) ) {
 188     int d = atoi(p);
 189     if (d != datasize) {
 190       itkWarningMacro(<< "Size " << datasize << " is not "
 191                       << dim1 << 'x' << dim2 << "x" << datalen
 192                       << " = " << d << ". Supposing the latter.");
 193       }
 194     }
 195
 196   // EHF files: binary data are in another file than the header file
 197   m_BinaryFileName = m_FileName;
 198   m_BinaryFileSkip = header_size;
 199   if (otherfile_name) {
 200     m_BinaryFileName = std::string(otherfile_name);
 201     m_BinaryFileSkip = otherfile_skip;
 202     }
 203
 204   double spacing = 1.;
 205   if ( (p = edf_findInHeader(header, "optic_used") ) )
 206     {
 207     spacing = atof(p);
 208     if(spacing == 0.)
 209       spacing = 1.;
 210     }
 211
 212
 213   free(header);
 214   gzclose(inp);
 215
 216   SetNumberOfDimensions(2);
 217   SetDimensions(0, dim1);
 218   SetDimensions(1, dim2);
 219   SetSpacing(0, spacing);
 220   SetSpacing(1, spacing);
 221   SetOrigin(0, 0.);
 222   SetOrigin(1, 0.);
 223 } ////
 224
 225 //--------------------------------------------------------------------
 226 // Read Image Information
 227 bool rtk::EdfImageIO::CanReadFile(const char* FileNameToRead)
 228 {
 229   std::string                  filename(FileNameToRead);
 230   const std::string::size_type it = filename.find_last_of( "." );
 231   std::string                  fileExt( filename, it+1, filename.length() );
 232
 233   if (fileExt != std::string("edf") ) return false;
 234   return true;
 235 } ////
 236
 237 //--------------------------------------------------------------------
 238 // Read Image Content
 239 void rtk::EdfImageIO::Read(void * buffer)
 240 {
 241   gzFile inp;
 242
 243   inp = gzopen(m_BinaryFileName.c_str(), "rb");
 244   if (!inp)
 245     itkGenericExceptionMacro(<< "Cannot open file \"" << m_FileName << "\"");
 246   gzseek(inp, m_BinaryFileSkip, SEEK_SET);
 247
 248   // read the data (image)
 249   long numberOfBytesToBeRead = GetComponentSize();
 250   for(unsigned int i=0; i<GetNumberOfDimensions(); i++) numberOfBytesToBeRead *= GetDimensions(i);
 251
 252   if (numberOfBytesToBeRead != gzread(inp, buffer, numberOfBytesToBeRead) )
 253     itkGenericExceptionMacro(<< "The image " << m_BinaryFileName << " cannot be read completely.");
 254
 255   gzclose(inp);
 256
 257   // Adapted from itkRawImageIO
 258     {
 259     using namespace itk;
 260     // Swap bytes if necessary
 261     if itkReadRawBytesAfterSwappingMacro( unsigned short, USHORT )
 262     else if itkReadRawBytesAfterSwappingMacro( short, SHORT )
 263     else if itkReadRawBytesAfterSwappingMacro( char, CHAR )
 264     else if itkReadRawBytesAfterSwappingMacro( unsigned char, UCHAR )
 265     else if itkReadRawBytesAfterSwappingMacro( unsigned int, UINT )
 266     else if itkReadRawBytesAfterSwappingMacro( int, INT )
 267     else if itkReadRawBytesAfterSwappingMacro( unsigned int, UINT )
 268     else if itkReadRawBytesAfterSwappingMacro( int, INT )
 269     else if itkReadRawBytesAfterSwappingMacro( float, FLOAT )
 270     else if itkReadRawBytesAfterSwappingMacro( double, DOUBLE );
 271     }
 272 }
 273
 274 //--------------------------------------------------------------------
 275 // Write Image Information
 276 void rtk::EdfImageIO::WriteImageInformation( bool itkNotUsed(keepOfStream) )
 277 {
 278 }
 279
 280 //--------------------------------------------------------------------
 281 // Write Image Information
 282 bool rtk::EdfImageIO::CanWriteFile( const char* itkNotUsed(FileNameToWrite) )
 283 {
 284   return false;
 285 }
 286
 287 //--------------------------------------------------------------------
 288 // Write Image
 289 void rtk::EdfImageIO::Write( const void * itkNotUsed(buffer) )
 290 {
 291 } ////