Line data    Source code

       1             : #include "EDFHandler.h"
       2             : #include "Logger.h"
       3             : #include "VolumeDescriptor.h"
       4             : 
       5             : #include <stdexcept>
       6             : 
       7             : namespace elsa
       8             : {
       9             :     template <typename data_t>
      10           0 :     DataContainer<data_t> EDF::read(std::string filename)
      11             :     {
      12           0 :         Logger::get("EDF")->info("Reading data from {}", filename);
      13             : 
      14             :         // open the file
      15           0 :         std::ifstream file(filename, std::ios::binary);
      16           0 :         if (!file.good())
      17           0 :             throw Error("EDF::read: cannot read from '" + filename + "'");
      18             : 
      19           0 :         return EDF::read<data_t>(file);
      20           0 :     }
      21             : 
      22             :     template <typename data_t>
      23           0 :     DataContainer<data_t> EDF::read(std::istream& input)
      24             :     {
      25             :         // get the meta data from the header
      26           0 :         auto properties = readHeader(input);
      27           0 :         auto [descriptor, dataType] = parseHeader(properties);
      28             : 
      29             :         // read in the data
      30           0 :         DataContainer<data_t> dataContainer(*descriptor);
      31             : 
      32           0 :         if (dataType == DataUtils::DataType::UINT16)
      33           0 :             DataUtils::parseRawData<uint16_t, data_t>(input, dataContainer);
      34           0 :         else if (dataType == DataUtils::DataType::FLOAT32)
      35           0 :             DataUtils::parseRawData<float, data_t>(input, dataContainer);
      36           0 :         else if (dataType == DataUtils::DataType::FLOAT64)
      37           0 :             DataUtils::parseRawData<double, data_t>(input, dataContainer);
      38             :         else
      39           0 :             throw Error("EDF::read: invalid/unsupported data type");
      40             : 
      41           0 :         return dataContainer;
      42           0 :     }
      43             : 
      44             :     template <typename data_t>
      45           0 :     void EDF::write(const DataContainer<data_t>& data, std::string filename)
      46             :     {
      47           0 :         Logger::get("EDF")->info("Writing data to {}", filename);
      48             : 
      49             :         // open the file
      50           0 :         std::ofstream file(filename, std::ios::binary);
      51           0 :         if (!file.good())
      52           0 :             throw Error("EDF::write: cannot write to '" + filename + "'");
      53             : 
      54           0 :         EDF::write<data_t>(data, file);
      55           0 :     }
      56             : 
      57             :     template <typename data_t>
      58           0 :     void EDF::write(const DataContainer<data_t>& data, std::ostream& output)
      59             :     {
      60             :         // output the header
      61           0 :         writeHeader(output, data);
      62             : 
      63             :         // output the raw data
      64             :         // TODO: this would be more efficient if we had a data pointer...
      65           0 :         for (index_t i = 0; i < data.getSize(); ++i)
      66           0 :             output.write(reinterpret_cast<const char*>(&data[i]), sizeof(data_t));
      67           0 :     }
      68             : 
      69           0 :     std::map<std::string, std::string> EDF::readHeader(std::istream& file)
      70             :     {
      71           0 :         std::map<std::string, std::string> properties;
      72             : 
      73             :         // read a single character and make sure that a header is opened
      74           0 :         if (file.eof() || file.get() != '{')
      75           0 :             throw InvalidArgumentError("EDF::readHeader: no header opening marker");
      76             : 
      77             :         // read header data
      78           0 :         while (!file.eof()) {
      79             :             // skip whitespace
      80           0 :             while (!file.eof()) {
      81           0 :                 const int chr = file.peek();
      82           0 :                 if (chr != '\r' && chr != '\n' && chr != ' ' && chr != '\t')
      83           0 :                     break;
      84           0 :                 file.ignore(1);
      85             :             }
      86             : 
      87             :             // abort if the header is closed
      88           0 :             if (file.eof() || file.peek() == '}')
      89           0 :                 break;
      90             : 
      91             :             // extract the property assignment
      92           0 :             bool quotesSingle = false, quotesDouble = false;
      93           0 :             std::string assignment;
      94           0 :             while (!file.eof()) {
      95           0 :                 auto chr = static_cast<char>(file.get());
      96             : 
      97             :                 // abort on end-of-assignment
      98           0 :                 if (chr == ';' && !(quotesSingle || quotesDouble))
      99           0 :                     break;
     100             : 
     101             :                 // check for quote characters
     102           0 :                 if (chr == '\'')
     103           0 :                     quotesSingle = !quotesSingle;
     104           0 :                 if (chr == '\"')
     105           0 :                     quotesDouble = !quotesDouble;
     106             : 
     107           0 :                 assignment += chr;
     108             :             }
     109             : 
     110             :             // split the assignment
     111           0 :             auto delim = assignment.find('=');
     112           0 :             if (delim == std::string::npos)
     113           0 :                 throw InvalidArgumentError("failed reading name/value delimiter");
     114             : 
     115           0 :             std::string name = assignment.substr(0, delim);
     116           0 :             StringUtils::trim(name);
     117             : 
     118           0 :             std::string value = assignment.substr(delim + 1);
     119           0 :             StringUtils::trim(value);
     120             : 
     121             :             // remove quotes (if they exist)
     122           0 :             if (value[0] == value[value.size() - 1] && (value[0] == '\'' || value[0] == '\"'))
     123           0 :                 value = value.substr(1, value.size() - 2);
     124             : 
     125           0 :             StringUtils::toLower(name);
     126           0 :             properties[name] = value;
     127           0 :         }
     128           0 :         file.ignore(2); // end of header marker
     129             : 
     130           0 :         return properties;
     131           0 :     }
     132             : 
     133             :     std::pair<std::unique_ptr<DataDescriptor>, DataUtils::DataType>
     134           0 :         EDF::parseHeader(const std::map<std::string, std::string>& properties)
     135             :     {
     136             :         // read the dimensions
     137           0 :         std::vector<index_t> dim;
     138           0 :         for (index_t i = 1;; i++) {
     139             :             // assemble the property name
     140           0 :             std::stringstream aux;
     141           0 :             aux << "dim_" << i;
     142             : 
     143             :             // try to find the property
     144           0 :             auto dimIt = properties.find(aux.str());
     145           0 :             if (dimIt == properties.end())
     146           0 :                 break;
     147             : 
     148           0 :             dim.push_back(DataUtils::parse<index_t>(dimIt->second));
     149           0 :         }
     150           0 :         const auto nDims = static_cast<index_t>(dim.size());
     151           0 :         if (nDims == 0u)
     152           0 :             throw Error("EDF::parseHeader: dimension information not found");
     153             : 
     154             :         // parse the (non-standard) spacing tag
     155           0 :         std::vector<real_t> spacing;
     156           0 :         auto spacingIt = properties.find("spacing");
     157           0 :         if (spacingIt != properties.end())
     158           0 :             spacing = DataUtils::parseVector<real_t>(spacingIt->second);
     159             : 
     160             :         // check for a byte order tag, but fall back to the default value
     161           0 :         auto byteorderIt = properties.find("byteorder");
     162           0 :         if (byteorderIt != properties.end()) {
     163           0 :             std::string byteorderValue = byteorderIt->second;
     164           0 :             StringUtils::toLower(byteorderValue);
     165             : 
     166           0 :             if (byteorderValue != "lowbytefirst")
     167           0 :                 throw Error("EDF::parseHeader: unsupported byte order value");
     168           0 :         }
     169             : 
     170             :         // check for the 'element type' value
     171             :         DataUtils::DataType dataType;
     172           0 :         auto datatypeIt = properties.find("datatype");
     173           0 :         if (datatypeIt != properties.end()) {
     174           0 :             std::string datatypeValue = datatypeIt->second;
     175           0 :             StringUtils::toLower(datatypeValue);
     176             : 
     177           0 :             if (datatypeValue == "signedbyte")
     178           0 :                 dataType = DataUtils::DataType::INT8;
     179           0 :             else if (datatypeValue == "unsignedbyte")
     180           0 :                 dataType = DataUtils::DataType::UINT8;
     181           0 :             else if (datatypeValue == "signedshort")
     182           0 :                 dataType = DataUtils::DataType::INT16;
     183           0 :             else if (datatypeValue == "unsignedshort")
     184           0 :                 dataType = DataUtils::DataType::UINT16;
     185           0 :             else if (datatypeValue == "float" || datatypeValue == "floatvalue"
     186           0 :                      || datatypeValue == "real")
     187           0 :                 dataType = DataUtils::DataType::FLOAT32;
     188           0 :             else if (datatypeValue == "double" || datatypeValue == "doublevalue")
     189           0 :                 dataType = DataUtils::DataType::FLOAT64;
     190             :             else
     191           0 :                 throw Error("EDF::parseHeader: invalid/unsupported data type");
     192           0 :         } else
     193           0 :             throw Error("EDF::parseHeader: data type not found");
     194             : 
     195           0 :         auto compressionIt = properties.find("compression");
     196           0 :         if (compressionIt != properties.end())
     197           0 :             throw Error("EDF::parseHeader: compression not supported");
     198             : 
     199           0 :         index_t size = 0;
     200           0 :         auto sizeIt = properties.find("size");
     201           0 :         if (sizeIt != properties.end())
     202           0 :             size = DataUtils::parse<index_t>(sizeIt->second);
     203             : 
     204           0 :         auto imageIt = properties.find("image");
     205           0 :         if (imageIt != properties.end() && DataUtils::parse<index_t>(imageIt->second) != 1)
     206           0 :             throw Error("EDF::parseHeader: image not set to 1");
     207             : 
     208             :         // convert size
     209           0 :         IndexVector_t dimSizeVec(nDims);
     210           0 :         for (index_t i = 0; i < nDims; ++i)
     211           0 :             dimSizeVec[i] = dim[static_cast<std::size_t>(i)];
     212           0 :         if (dimSizeVec.prod() * DataUtils::getSizeOfDataType(dataType) != size)
     213           0 :             throw Error("EDF::parseHeader: size inconsistency");
     214             : 
     215             :         // convert spacing
     216           0 :         RealVector_t dimSpacingVec(RealVector_t::Ones(nDims));
     217           0 :         if (!spacing.empty()) {
     218           0 :             if (nDims != static_cast<index_t>(spacing.size()))
     219           0 :                 throw Error("EDF::parseHeader: spacing inconsistency");
     220           0 :             for (index_t i = 0; i < nDims; ++i)
     221           0 :                 dimSpacingVec[i] = spacing[static_cast<std::size_t>(i)];
     222             :         }
     223             : 
     224           0 :         return std::make_pair(std::make_unique<VolumeDescriptor>(dimSizeVec, dimSpacingVec),
     225           0 :                               dataType);
     226           0 :     }
     227             : 
     228             :     template <typename data_t>
     229           0 :     void EDF::writeHeader(std::ostream& file, const DataContainer<data_t>& data)
     230             :     {
     231             :         // open the header
     232           0 :         file << "{\n";
     233             : 
     234           0 :         file << "HeaderID = EH:000001:000000:000000;\n";
     235           0 :         file << "Image = " << 1 << ";\n";
     236           0 :         file << "ByteOrder = LowByteFirst;\n";
     237           0 :         file << "DataType = " << getDataTypeName(data) << ";\n";
     238             : 
     239           0 :         auto& descriptor = data.getDataDescriptor();
     240             : 
     241             :         // write dimension and size
     242           0 :         for (index_t i = 0; i < descriptor.getNumberOfDimensions(); ++i)
     243           0 :             file << "Dim_" << (i + 1) << " = "
     244           0 :                  << descriptor.getNumberOfCoefficientsPerDimension()[i] << ";\n";
     245           0 :         file << "Size = "
     246           0 :              << descriptor.getNumberOfCoefficients() * static_cast<index_t>(sizeof(data_t))
     247           0 :              << ";\n";
     248             : 
     249             :         // write spacing
     250           0 :         file << "Spacing =";
     251           0 :         for (index_t i = 0; i < descriptor.getNumberOfDimensions(); ++i)
     252           0 :             file << ' ' << descriptor.getSpacingPerDimension()[i];
     253           0 :         file << ";\n";
     254             : 
     255             :         // pad the header by adding spaces such that the header ends on a kilobyte boundary
     256           0 :         index_t n = 1024;
     257           0 :         while (n < (static_cast<index_t>(file.tellp()) + 3))
     258           0 :             n += 1024;
     259           0 :         n -= static_cast<index_t>(file.tellp()) + 3;
     260           0 :         while (n > 0) {
     261           0 :             file.put(' ');
     262           0 :             n--;
     263             :         }
     264             : 
     265             :         // close the header
     266           0 :         file << "\n}\n";
     267           0 :     }
     268             : 
     269             :     template <typename data_t>
     270           0 :     std::string EDF::getDataTypeName([[maybe_unused]] const DataContainer<data_t>& data)
     271             :     {
     272           0 :         throw InvalidArgumentError("EDF::getDataTypeName: invalid/unsupported data type");
     273             :     }
     274             : 
     275             :     template <>
     276           0 :     std::string EDF::getDataTypeName([[maybe_unused]] const DataContainer<int8_t>& data)
     277             :     {
     278           0 :         return "SignedByte";
     279             :     }
     280             :     template <>
     281           0 :     std::string EDF::getDataTypeName([[maybe_unused]] const DataContainer<uint8_t>& data)
     282             :     {
     283           0 :         return "UnsignedByte";
     284             :     }
     285             :     template <>
     286           0 :     std::string EDF::getDataTypeName([[maybe_unused]] const DataContainer<int16_t>& data)
     287             :     {
     288           0 :         return "SignedShort";
     289             :     }
     290             :     template <>
     291           0 :     std::string EDF::getDataTypeName([[maybe_unused]] const DataContainer<uint16_t>& data)
     292             :     {
     293           0 :         return "UnsignedShort";
     294             :     }
     295             :     template <>
     296           0 :     std::string EDF::getDataTypeName([[maybe_unused]] const DataContainer<float>& data)
     297             :     {
     298           0 :         return "FloatValue";
     299             :     }
     300             :     template <>
     301           0 :     std::string EDF::getDataTypeName([[maybe_unused]] const DataContainer<double>& data)
     302             :     {
     303           0 :         return "DoubleValue";
     304             :     }
     305             : 
     306             :     // ------------------------------------------
     307             :     // explicit template instantiation
     308             :     template DataContainer<float> EDF::read(std::string);
     309             :     template DataContainer<double> EDF::read(std::string);
     310             :     template DataContainer<index_t> EDF::read(std::string);
     311             : 
     312             :     template DataContainer<float> EDF::read(std::istream&);
     313             :     template DataContainer<double> EDF::read(std::istream&);
     314             :     template DataContainer<index_t> EDF::read(std::istream&);
     315             : 
     316             :     template void EDF::write(const DataContainer<float>&, std::string);
     317             :     template void EDF::write(const DataContainer<double>&, std::string);
     318             :     template void EDF::write(const DataContainer<index_t>&, std::string);
     319             : 
     320             :     template void EDF::write(const DataContainer<float>&, std::ostream&);
     321             :     template void EDF::write(const DataContainer<double>&, std::ostream&);
     322             :     template void EDF::write(const DataContainer<index_t>&, std::ostream&);
     323             : 
     324             :     template void EDF::writeHeader(std::ostream&, const DataContainer<float>&);
     325             :     template void EDF::writeHeader(std::ostream&, const DataContainer<double>&);
     326             :     template void EDF::writeHeader(std::ostream&, const DataContainer<index_t>&);
     327             : 
     328             : } // namespace elsa