Commit a3cb8a15 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Implemented parsing of binary data in VTKReader

fixes #70
parent e3a4ae49
Loading
Loading
Loading
Loading
+166 −101
Original line number Diff line number Diff line
@@ -17,7 +17,7 @@
#include <vector>

#include <TNL/Meshes/Readers/MeshReader.h>
#include <TNL/Exceptions/NotImplementedError.h>
#include <TNL/Endianness.h>

namespace TNL {
namespace Meshes {
@@ -39,33 +39,34 @@ public:

      std::ifstream inputFile( fileName );
      if( ! inputFile )
         throw MeshReaderError( "VTKReader", "failed to open the file '" + fileName + "'." );
         throw MeshReaderError( "VTKReader", "failed to open the file '" + fileName + "'" );

      if( ! parseHeader( inputFile ) )
         throw MeshReaderError( "VTKReader", "failed to parse the VTK file header." );
      const auto positionAfterHeading = inputFile.tellg();
      parseHeader( inputFile );

      if( dataset != "UNSTRUCTURED_GRID" )
         throw MeshReaderError( "VTKReader", "the dataset '" + dataset + "' is not supported." );
         throw MeshReaderError( "VTKReader", "the dataset '" + dataset + "' is not supported" );

      // TODO: implement binary parsing
      if( dataType == "BINARY" )
         throw Exceptions::NotImplementedError("VTKReader: parsing of BINARY data is not implemented yet.");
      // parse the file, find the starting positions of all relevant sections
      findSections( inputFile );

      std::string line, aux;
      std::istringstream iss;

      // parse points section
      if( ! findSection( inputFile, "POINTS" ) )
         throw MeshReaderError( "VTKReader", "unable to find the POINTS section, the file may be invalid or corrupted." );
      if( ! sectionPositions.count( "POINTS" ) )
         throw MeshReaderError( "VTKReader", "unable to find the POINTS section, the file may be invalid or corrupted" );
      inputFile.seekg( sectionPositions["POINTS"] );
      getline( inputFile, line );
      iss.clear();
      iss.str( line );
      iss >> aux;
      iss >> NumberOfPoints;
      iss >> pointsType;
      if( pointsType != "float" && pointsType != "double" )
         throw MeshReaderError( "VTKReader", "unsupported data type for POINTS: " + pointsType );

      // global index type is not stored in legacy VTK files
      // (binary VTK files don't support int64)
      connectivityType = offsetsType = "std::int32_t";
      // only std::uint8_t makes sense for entity types
      typesType = "std::uint8_t";
@@ -78,22 +79,18 @@ public:
      // read points
      worldDimension = 0;
      for( std::size_t pointIndex = 0; pointIndex < NumberOfPoints; pointIndex++ ) {
         if( ! inputFile ) {
            reset();
            throw MeshReaderError( "VTKReader", "unable to read enough vertices, the file may be invalid or corrupted." );
         }
         getline( inputFile, line );
         if( ! inputFile )
            throw MeshReaderError( "VTKReader", "unable to read enough vertices, the file may be invalid or corrupted" );

         // read the coordinates and compute the world dimension
         iss.clear();
         iss.str( line );
         for( int i = 0; i < 3; i++ ) {
            double aux;
            iss >> aux;
            if( ! iss ) {
               reset();
               throw MeshReaderError( "VTKReader", "unable to read " + std::to_string(i) + "th component of the vertex number " + std::to_string(pointIndex) + "." );
            }
            double aux = 0;
            if( pointsType == "float" )
               aux = readValue< float >( dataFormat, inputFile );
            else
               aux = readValue< double >( dataFormat, inputFile );
            if( ! inputFile )
               throw MeshReaderError( "VTKReader", "unable to read " + std::to_string(i) + "th component of the vertex number " + std::to_string(pointIndex) );
            if( aux != 0.0 )
               worldDimension = std::max( worldDimension, i + 1 );
            pointsArray.push_back( aux );
@@ -101,30 +98,22 @@ public:
      }

      // skip to the CELL_TYPES section
      if( ! findSection( inputFile, "CELL_TYPES" ) ) {
         reset();
         throw MeshReaderError( "VTKReader", "unable to find the CELL_TYPES section, the file may be invalid or corrupted." );
      }
      if( ! sectionPositions.count( "CELL_TYPES" ) )
         throw MeshReaderError( "VTKReader", "unable to find the CELL_TYPES section, the file may be invalid or corrupted" );
      inputFile.seekg( sectionPositions["CELL_TYPES"] );
      getline( inputFile, line );
      iss.clear();
      iss.str( line );
      iss >> aux;
      std::size_t NumberOfEntities;
      std::size_t NumberOfEntities = 0;
      iss >> NumberOfEntities;

      // read entity types
      for( std::size_t entityIndex = 0; entityIndex < NumberOfEntities; entityIndex++ ) {
         if( ! inputFile ) {
            reset();
            throw MeshReaderError( "VTKReader", "unable to read enough cell types, the file may be invalid or corrupted." );
         }
         getline( inputFile, line );

         // get entity type
         int typeId;
         iss.clear();
         iss.str( line );
         iss >> typeId;
         if( ! inputFile )
            throw MeshReaderError( "VTKReader", "unable to read enough cell types, the file may be invalid or corrupted" );
         // cell types are stored with great redundancy as int32 in the VTK file
         const std::uint8_t typeId = readValue< std::int32_t >( dataFormat, inputFile );
         typesArray.push_back( typeId );
      }

@@ -146,10 +135,8 @@ public:
         }
      }

      if( meshDimension == 0 ) {
         reset();
      if( meshDimension == 0 )
         throw MeshReaderError( "VTKReader", "Mesh dimension cannot be 0. Are there any entities at all?" );
      }

      // filter out cell shapes
      std::vector< std::uint8_t > cellTypes;
@@ -165,7 +152,6 @@ public:
         const std::string msg = "invalid number of cells (" + std::to_string(NumberOfCells) + "). Counts of entities for each dimension (0,1,2,3) are: "
                               + std::to_string(entitiesCounts[0]) + ", " + std::to_string(entitiesCounts[1]) + ", "
                               + std::to_string(entitiesCounts[2]) + ", " + std::to_string(entitiesCounts[3]);
         reset();
         throw MeshReaderError( "VTKReader", msg );
      }

@@ -174,41 +160,31 @@ public:
      for( auto c : cellTypes )
         if( (VTK::EntityShape) c != cellShape ) {
            const std::string msg = "Mixed unstructured meshes are not supported. There are cells with type "
                                  + VTK::getShapeName(cellShape) + " and " + VTK::getShapeName((VTK::EntityShape) c) + ".";
            reset();
                                  + VTK::getShapeName(cellShape) + " and " + VTK::getShapeName((VTK::EntityShape) c);
            throw MeshReaderError( "VTKReader", msg );
         }

      // find to the CELLS section
      if( ! findSection( inputFile, "CELLS", positionAfterHeading ) ) {
         reset();
         throw MeshReaderError( "VTKReader", "unable to find the CELLS section, the file may be invalid or corrupted." );
      }
      if( ! sectionPositions.count( "CELLS" ) )
         throw MeshReaderError( "VTKReader", "unable to find the CELLS section, the file may be invalid or corrupted" );
      inputFile.seekg( sectionPositions["CELLS"] );
      getline( inputFile, line );

      // read entities
      for( std::size_t entityIndex = 0; entityIndex < NumberOfEntities; entityIndex++ ) {
         if( ! inputFile ) {
            reset();
            throw MeshReaderError( "VTKReader", "unable to read enough cells, the file may be invalid or corrupted."
         if( ! inputFile )
            throw MeshReaderError( "VTKReader", "unable to read enough cells, the file may be invalid or corrupted"
                                                " (entityIndex = " + std::to_string(entityIndex) + ")" );
         }
         getline( inputFile, line );

         if( (VTK::EntityShape) typesArray[ entityIndex ] == cellShape ) {
            iss.clear();
            iss.str( line );
            // read number of subvertices
            int subvertices = 0;
            iss >> subvertices;
            const std::int32_t subvertices = readValue< std::int32_t >( dataFormat, inputFile );
            for( int v = 0; v < subvertices; v++ ) {
               std::size_t vid;
               iss >> vid;
               if( ! iss ) {
                  reset();
                  throw MeshReaderError( "VTKReader", "unable to read enough cells, the file may be invalid or corrupted."
               // legacy VTK files do not support 64-bit integers, even in the BINARY format
               const std::int32_t vid = readValue< std::int32_t >( dataFormat, inputFile );
               if( ! inputFile )
                  throw MeshReaderError( "VTKReader", "unable to read enough cells, the file may be invalid or corrupted"
                                                      " (entityIndex = " + std::to_string(entityIndex) + ", subvertex = " + std::to_string(v) + ")" );
               }
               connectivityArray.push_back( vid );
            }
            offsetsArray.push_back( connectivityArray.size() );
@@ -231,79 +207,168 @@ public:
   virtual void reset() override
   {
      resetBase();
      dataType = "";
      dataFormat = VTK::FileFormat::ascii;
      dataset = "";
      sectionPositions.clear();
   }

protected:
   // output of parseHeader
   std::string dataType;
   VTK::FileFormat dataFormat = VTK::FileFormat::ascii;
   std::string dataset;

   bool parseHeader( std::istream& str )
   // output of findSections
   std::map< std::string, std::ios::pos_type > sectionPositions;

   void parseHeader( std::istream& str )
   {
      std::string line;
      std::istringstream iss;

      // check header
      getline( str, line );
      if( line != "# vtk DataFile Version 2.0" ) {
         std::cerr << "VTKReader: unsupported VTK header: '" << line << "'." << std::endl;
         return false;
      }
      if( line != "# vtk DataFile Version 2.0" )
         throw MeshReaderError( "VTKReader", "failed to parse the VTK file header: unsupported VTK header '" + line + "'" );

      // skip title
      if( ! str )
         return false;
         throw MeshReaderError( "VTKReader", "failed to parse the VTK file header" );
      getline( str, line );

      // parse data type
      if( ! str )
         return false;
      getline( str, dataType );
      if( dataType != "ASCII" && dataType != "BINARY" ) {
         std::cerr << "VTKReader: unknown data type: '" << dataType << "'." << std::endl;
         return false;
      }
         throw MeshReaderError( "VTKReader", "failed to parse the VTK file header" );
      std::string format;
      getline( str, format );
      if( format == "ASCII" )
         dataFormat = VTK::FileFormat::ascii;
      else if( format == "BINARY" )
         dataFormat = VTK::FileFormat::binary;
      else
         throw MeshReaderError( "VTKReader", "unknown data format: '" + format + "'" );

      // parse dataset
      if( ! str )
         return false;
         throw MeshReaderError( "VTKReader", "failed to parse the VTK file header" );
      getline( str, line );
      iss.clear();
      iss.str( line );
      std::string tmp;
      iss >> tmp;
      if( tmp != "DATASET" ) {
         std::cerr << "VTKReader: wrong dataset specification: '" << line << "'." << std::endl;
         return false;
      }
      if( tmp != "DATASET" )
         throw MeshReaderError( "VTKReader", "wrong dataset specification: '" + line + "'" );
      iss >> dataset;

      return true;
   }

   bool findSection( std::istream& str, const std::string& section, std::ios::pos_type begin = -1 )
   void findSections( std::istream& str )
   {
      std::string line, aux;
      std::istringstream iss;

      if( begin >= 0 )
         str.seekg( begin );

      while( str ) {
         std::ios::pos_type currentPosition = str.tellg();
         // drop all whitespace (empty lines etc) before saving a position and reading a line
         str >> std::ws;
         if( str.eof() )
            break;

         // read a line which should contain the following section header
         const std::ios::pos_type currentPosition = str.tellg();
         std::string line;
         getline( str, line );
         if( ! str )
            throw MeshReaderError( "VTKReader", "failed to parse sections of the VTK file" );

         // parse the section name
         std::istringstream iss( line );
         std::string name;
         iss >> name;

         if( name == "FIELD" ) {
            sectionPositions.insert( {"FIELD", currentPosition} );
            // parse the rest of the line: FIELD FieldData <count>
            std::string aux;
            int count = 0;
            iss >> aux >> count;
            // skip the FieldData arrays
            for( int i = 0; i < count; i++ ) {
               getline( str, line );
               iss.clear();
               iss.str( line );
         iss >> aux;
         if( aux == section ) {
            str.seekg( currentPosition );
            return true;
               // <name> <components> <tuples> <datatype>
               std::int32_t components = 0;
               std::int32_t tuples = 0;
               std::string datatype;
               iss >> aux >> components >> tuples >> datatype;
               if( ! iss )
                  throw MeshReaderError( "VTKReader", "failed to extract FieldData information from line '" + line + "'" );
               // skip the points coordinates
               for( std::int32_t j = 0; j < components * tuples; j++ )
                  skipValue( dataFormat, str, datatype );
            }
         }
         else if( name == "POINTS" ) {
            sectionPositions.insert( {"POINTS", currentPosition} );
            // parse the rest of the line: POINTS <count> <datatype>
            std::int32_t count = 0;
            std::string datatype;
            iss >> count >> datatype;
            // skip the values
            for( std::int32_t j = 0; j < 3 * count; j++ )
               skipValue( dataFormat, str, datatype );
         }
         else if( name == "CELLS" ) {
            sectionPositions.insert( {"CELLS", currentPosition} );
            // parse the rest of the line: CELLS <cells_count> <values_count>
            std::int32_t cells_count = 0;
            std::int32_t values_count = 0;
            iss >> cells_count >> values_count;
            // skip the values
            for( std::int32_t j = 0; j < values_count; j++ )
               skipValue( dataFormat, str, "int" );
         }
         else if( name == "CELL_TYPES" ) {
            sectionPositions.insert( {"CELL_TYPES", currentPosition} );
            // parse the rest of the line: CELL_TYPES <count>
            std::int32_t count = 0;
            iss >> count;
            // skip the values
            for( std::int32_t j = 0; j < count; j++ )
               // cell types are stored with great redundancy as int32 in the VTK file
               skipValue( dataFormat, str, "int" );
         }
         else if( name == "CELL_DATA" || name == "POINT_DATA" ) {
            // TODO: implement reading values from these sections
            break;
         }
         else
            throw MeshReaderError( "VTKReader", "parsing error: unexpected section start at byte " + std::to_string(currentPosition)
                                    + " (section name is '" + name + "')" );
      }
   }

   static void skipValue( VTK::FileFormat format, std::istream& str, std::string datatype )
   {
      if( datatype == "int" )
         readValue< std::int32_t >( format, str );
      else if( datatype == "float" )
         readValue< float >( format, str );
      else if( datatype == "double" )
         readValue< double >( format, str );
      else
         throw MeshReaderError( "VTKReader", "found data type which is not implemented in the reader: " + datatype );
   }

      return false;
   template< typename T >
   static T readValue( VTK::FileFormat format, std::istream& str )
   {
      T value;
      if( format == VTK::FileFormat::binary ) {
         str.read( reinterpret_cast<char*>(&value), sizeof(T) );
         // forceBigEndian = swapIfLittleEndian, i.e. here it forces a big-endian
         // value to the correct system endianness
         value = forceBigEndian( value );
      }
      else {
         str >> value;
      }
      return value;
   }
};