diff --git a/CMakeLists.txt b/CMakeLists.txt index b90c4ac..6db9314 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -342,7 +342,7 @@ set(PUBLICHEADERS include/aare/File.hpp include/aare/Fit.hpp include/aare/FileInterface.hpp - include/aare/FiltPtr.hpp + include/aare/FilePtr.hpp include/aare/Frame.hpp include/aare/geo_helpers.hpp include/aare/JungfrauDataFile.hpp diff --git a/docs/src/JungfrauDataFile.rst b/docs/src/JungfrauDataFile.rst index 6fb7dab..78d473f 100644 --- a/docs/src/JungfrauDataFile.rst +++ b/docs/src/JungfrauDataFile.rst @@ -1,6 +1,23 @@ JungfrauDataFile ================== +JungfrauDataFile is a class to read the .dat files that are produced by Aldo's receiver. +It is mostly used for calibration. + +The structure of the file is: + +* JungfrauDataHeader +* Binary data (256x256, 256x1024 or 512x1024) +* JungfrauDataHeader +* ... + +There is no metadata indicating number of frames or the size of the image, but this +will be infered by this reader. + +.. doxygenstruct:: aare::JungfrauDataHeader + :members: + :undoc-members: + :private-members: .. doxygenclass:: aare::JungfrauDataFile :members: diff --git a/include/aare/JungfrauDataFile.hpp b/include/aare/JungfrauDataFile.hpp index 301cb9f..bba5403 100644 --- a/include/aare/JungfrauDataFile.hpp +++ b/include/aare/JungfrauDataFile.hpp @@ -6,77 +6,101 @@ #include "aare/FilePtr.hpp" #include "aare/defs.hpp" #include "aare/NDArray.hpp" +#include "aare/FileInterface.hpp" namespace aare { + struct JungfrauDataHeader{ uint64_t framenum; uint64_t bunchid; }; -class JungfrauDataFile { +class JungfrauDataFile : public FileInterface { - size_t m_rows{}; - size_t m_cols{}; - size_t m_bytes_per_frame{}; - size_t m_total_frames{}; - size_t m_offset{}; // file index of the first file, allow starting at non zero file - size_t m_current_file_index{}; - size_t m_current_frame{}; + size_t m_rows{}; //!< number of rows in the image, from find_frame_size(); + size_t m_cols{}; //!< number of columns in the image, from find_frame_size(); + size_t m_bytes_per_frame{}; //!< number of bytes per frame excluding header + size_t m_total_frames{}; //!< total number of frames in the series of files + size_t m_offset{}; //!< file index of the first file, allow starting at non zero file + size_t m_current_file_index{}; //!< The index of the open file + size_t m_current_frame_index{}; //!< The index of the current frame (with reference to all files) - std::vector m_frame_index{}; + std::vector m_last_frame_in_file{}; //!< Used for seeking to the correct file + std::filesystem::path m_path; //!< path to the files + std::string m_base_name; //!< base name used for formatting file names - std::filesystem::path m_path; - std::string m_base_name; + FilePtr m_fp; //!< RAII wrapper for a FILE* - FilePtr m_fp; - - //Data sizes to guess the frame size + using pixel_type = uint16_t; static constexpr size_t header_size = sizeof(JungfrauDataHeader); - static constexpr size_t module_data_size = header_size + sizeof(pixel_type) * 512 * 1024; - static constexpr size_t half_data_size = header_size + sizeof(pixel_type) * 256 * 1024; - static constexpr size_t chip_data_size = header_size + sizeof(pixel_type) * 256 * 256; - static constexpr size_t n_digits_in_file_index = 6; // number of digits in the file index + static constexpr size_t n_digits_in_file_index = 6; //!< to format file names public: JungfrauDataFile(const std::filesystem::path &fname); std::string base_name() const; //!< get the base name of the file (without path and extension) - size_t bytes_per_frame() const; - size_t pixels_per_frame() const; + size_t bytes_per_frame() override; + size_t pixels_per_frame() override; size_t bytes_per_pixel() const; - size_t bitdepth() const; - void seek(size_t frame_index); //!< seek to the given frame index - size_t tell() const; //!< get the frame index of the file pointer - size_t total_frames() const; - size_t rows() const; - size_t cols() const; - size_t n_files() const { return m_frame_index.size(); } //!< get the number of files + size_t bitdepth() const override; + void seek(size_t frame_index) override; //!< seek to the given frame index (note not byte offset) + size_t tell() override; //!< get the frame index of the file pointer + size_t total_frames() const override; + size_t rows() const override; + size_t cols() const override; + size_t n_files() const; //!< get the number of files in the series. - void read_into(std::byte *image_buf, JungfrauDataHeader *header); - void read_into(std::byte *image_buf, size_t n_frames, JungfrauDataHeader *header); - void read_into(NDArray &image, JungfrauDataHeader &header) { - if(!(rows() == image.shape(0) && cols() == image.shape(1))){ - throw std::runtime_error(LOCATION + - "Image shape does not match file size: " + std::to_string(rows()) + "x" + std::to_string(cols())); - } - read_into(reinterpret_cast(image.data()), &header); - } - NDArray read_frame(JungfrauDataHeader& header) { - Shape<2> shape{rows(), cols()}; - NDArray image(shape); + // Extra functions needed for FileInterface + Frame read_frame() override; + Frame read_frame(size_t frame_number) override; + std::vector read_n(size_t n_frames=0) override; + void read_into(std::byte *image_buf) override; + void read_into(std::byte *image_buf, size_t n_frames) override; + size_t frame_number(size_t frame_index) override; + DetectorType detector_type() const override; - read_into(reinterpret_cast(image.data()), - &header); + /** + * @brief Read a single frame from the file into the given buffer. + * @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer) + * @param header pointer to a JungfrauDataHeader or nullptr to skip header) + */ + void read_into(std::byte *image_buf, JungfrauDataHeader *header = nullptr); - return image; - } + /** + * @brief Read a multiple frames from the file into the given buffer. + * @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer) + * @param n_frames number of frames to read + * @param header pointer to a JungfrauDataHeader or nullptr to skip header) + */ + void read_into(std::byte *image_buf, size_t n_frames, JungfrauDataHeader *header = nullptr); + + /** + * @brief Read a single frame from the file into the given NDArray + * @param image NDArray to read the frame into. + */ + void read_into(NDArray* image, JungfrauDataHeader* header = nullptr); + /** + * @brief Read a single frame from the file. Allocated a new NDArray for the output data + * @param header pointer to a JungfrauDataHeader or nullptr to skip header) + * @return NDArray with the image data + */ + NDArray read_frame(JungfrauDataHeader* header = nullptr); + + JungfrauDataHeader read_header(); std::filesystem::path current_file() const { return fpath(m_current_file_index+m_offset); } private: + /** + * @brief Find the size of the frame in the file. (256x256, 256x1024, 512x1024) + * @param fname path to the file + * @throws std::runtime_error if the file is empty or the size cannot be determined + */ void find_frame_size(const std::filesystem::path &fname); + + void parse_fname(const std::filesystem::path &fname); void scan_files(); void open_file(size_t file_index); diff --git a/python/src/jungfrau_data_file.hpp b/python/src/jungfrau_data_file.hpp index 47e64f1..13a3202 100644 --- a/python/src/jungfrau_data_file.hpp +++ b/python/src/jungfrau_data_file.hpp @@ -60,8 +60,14 @@ void define_jungfrau_data_file_io_bindings(py::module &m) { py::class_(m, "JungfrauDataFile") .def(py::init()) - .def("seek", &JungfrauDataFile::seek) - .def("tell", &JungfrauDataFile::tell) + .def("seek", &JungfrauDataFile::seek, + R"( + Seek to the given frame index. + )") + .def("tell", &JungfrauDataFile::tell, + R"( + Get the current frame index. + )") .def_property_readonly("rows", &JungfrauDataFile::rows) .def_property_readonly("cols", &JungfrauDataFile::cols) .def_property_readonly("base_name", &JungfrauDataFile::base_name) @@ -75,12 +81,11 @@ void define_jungfrau_data_file_io_bindings(py::module &m) { .def_property_readonly("current_file", &JungfrauDataFile::current_file) .def_property_readonly("total_frames", &JungfrauDataFile::total_frames) .def_property_readonly("n_files", &JungfrauDataFile::n_files) - .def("read_frame", - [](JungfrauDataFile &self) { return read_dat_frame(self); }) - .def("read_n", - [](JungfrauDataFile &self, size_t n_frames) { - return read_n_dat_frames(self, n_frames); - }, + .def("read_frame", &read_dat_frame, + R"( + Read a single frame from the file. + )") + .def("read_n", &read_n_dat_frames, R"( Read maximum n_frames frames from the file. )") diff --git a/src/File.cpp b/src/File.cpp index 3c68eff..eb04893 100644 --- a/src/File.cpp +++ b/src/File.cpp @@ -1,4 +1,5 @@ #include "aare/File.hpp" +#include "aare/JungfrauDataFile.hpp" #include "aare/NumpyFile.hpp" #include "aare/RawFile.hpp" @@ -27,6 +28,8 @@ File::File(const std::filesystem::path &fname, const std::string &mode, else if (fname.extension() == ".npy") { // file_impl = new NumpyFile(fname, mode, cfg); file_impl = std::make_unique(fname, mode, cfg); + }else if(fname.extension() == ".dat"){ + file_impl = std::make_unique(fname); } else { throw std::runtime_error("Unsupported file type"); } diff --git a/src/JungfrauDataFile.cpp b/src/JungfrauDataFile.cpp index 20ce969..6e1ccd6 100644 --- a/src/JungfrauDataFile.cpp +++ b/src/JungfrauDataFile.cpp @@ -1,9 +1,9 @@ #include "aare/JungfrauDataFile.hpp" -#include "aare/defs.hpp" #include "aare/algorithm.hpp" +#include "aare/defs.hpp" -#include #include +#include namespace aare { @@ -19,42 +19,89 @@ JungfrauDataFile::JungfrauDataFile(const std::filesystem::path &fname) { open_file(m_current_file_index); } -std::string JungfrauDataFile::base_name() const { return m_base_name; } -size_t JungfrauDataFile::bytes_per_frame() const { - return m_bytes_per_frame; +// FileInterface + +Frame JungfrauDataFile::read_frame(){ + Frame f(rows(), cols(), Dtype::UINT16); + read_into(reinterpret_cast(f.data()), nullptr); + return f; } -size_t JungfrauDataFile::pixels_per_frame() const { return m_rows * m_cols; } +Frame JungfrauDataFile::read_frame(size_t frame_number){ + seek(frame_number); + Frame f(rows(), cols(), Dtype::UINT16); + read_into(reinterpret_cast(f.data()), nullptr); + return f; +} + +std::vector JungfrauDataFile::read_n(size_t n_frames) { + std::vector frames; + throw std::runtime_error(LOCATION + + "Not implemented yet"); + return frames; +} + +void JungfrauDataFile::read_into(std::byte *image_buf) { + read_into(image_buf, nullptr); +} +void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames) { + read_into(image_buf, n_frames, nullptr); +} + +size_t JungfrauDataFile::frame_number(size_t frame_index) { + seek(frame_index); + return read_header().framenum; +} + +DetectorType JungfrauDataFile::detector_type() const { return DetectorType::Jungfrau; } + +std::string JungfrauDataFile::base_name() const { return m_base_name; } + +size_t JungfrauDataFile::bytes_per_frame() { return m_bytes_per_frame; } + +size_t JungfrauDataFile::pixels_per_frame() { return m_rows * m_cols; } size_t JungfrauDataFile::bytes_per_pixel() const { return sizeof(pixel_type); } -size_t JungfrauDataFile::bitdepth() const { return bytes_per_pixel()*bits_per_byte; } +size_t JungfrauDataFile::bitdepth() const { + return bytes_per_pixel() * bits_per_byte; +} void JungfrauDataFile::seek(size_t frame_index) { if (frame_index >= m_total_frames) { - throw std::runtime_error(LOCATION + - "Frame index out of range: " + std::to_string(frame_index)); + throw std::runtime_error(LOCATION + "Frame index out of range: " + + std::to_string(frame_index)); } - m_current_frame = frame_index; - auto file_index = first_larger(m_frame_index, frame_index); + m_current_frame_index = frame_index; + auto file_index = first_larger(m_last_frame_in_file, frame_index); - if(file_index != m_current_file_index) + if (file_index != m_current_file_index) open_file(file_index); - auto frame_offset = (file_index) ? frame_index - m_frame_index[file_index-1] : frame_index; + auto frame_offset = (file_index) + ? frame_index - m_last_frame_in_file[file_index - 1] + : frame_index; auto byte_offset = frame_offset * (m_bytes_per_frame + header_size); m_fp.seek(byte_offset); -}; +}; -size_t JungfrauDataFile::tell() const { - return m_current_frame; -} +size_t JungfrauDataFile::tell() { return m_current_frame_index; } size_t JungfrauDataFile::total_frames() const { return m_total_frames; } size_t JungfrauDataFile::rows() const { return m_rows; } size_t JungfrauDataFile::cols() const { return m_cols; } +size_t JungfrauDataFile::n_files() const { return m_last_frame_in_file.size(); } + void JungfrauDataFile::find_frame_size(const std::filesystem::path &fname) { + + static constexpr size_t module_data_size = + header_size + sizeof(pixel_type) * 512 * 1024; + static constexpr size_t half_data_size = + header_size + sizeof(pixel_type) * 256 * 1024; + static constexpr size_t chip_data_size = + header_size + sizeof(pixel_type) * 256 * 256; + auto file_size = std::filesystem::file_size(fname); if (file_size == 0) { throw std::runtime_error(LOCATION + @@ -64,15 +111,15 @@ void JungfrauDataFile::find_frame_size(const std::filesystem::path &fname) { if (file_size % module_data_size == 0) { m_rows = 512; m_cols = 1024; - m_bytes_per_frame = module_data_size-header_size; + m_bytes_per_frame = module_data_size - header_size; } else if (file_size % half_data_size == 0) { m_rows = 256; m_cols = 1024; - m_bytes_per_frame = half_data_size-header_size; + m_bytes_per_frame = half_data_size - header_size; } else if (file_size % chip_data_size == 0) { m_rows = 256; m_cols = 256; - m_bytes_per_frame = chip_data_size-header_size; + m_bytes_per_frame = chip_data_size - header_size; } else { throw std::runtime_error(LOCATION + "Cannot find frame size: file size is not a " @@ -86,74 +133,107 @@ void JungfrauDataFile::parse_fname(const std::filesystem::path &fname) { // find file index, then remove if from the base name if (auto pos = m_base_name.find_last_of('_'); pos != std::string::npos) { - m_offset = - std::stoul(m_base_name.substr(pos+1)); + m_offset = std::stoul(m_base_name.substr(pos + 1)); m_base_name.erase(pos); } } void JungfrauDataFile::scan_files() { // find how many files we have and the number of frames in each file - m_frame_index.clear(); + m_last_frame_in_file.clear(); size_t file_index = m_offset; while (std::filesystem::exists(fpath(file_index))) { - auto n_frames = - std::filesystem::file_size(fpath(file_index)) / (m_bytes_per_frame + - header_size); - m_frame_index.push_back(n_frames); + auto n_frames = std::filesystem::file_size(fpath(file_index)) / + (m_bytes_per_frame + header_size); + m_last_frame_in_file.push_back(n_frames); ++file_index; } // find where we need to open the next file and total number of frames - m_frame_index = cumsum(m_frame_index); - m_total_frames = m_frame_index.back(); - + m_last_frame_in_file = cumsum(m_last_frame_in_file); + m_total_frames = m_last_frame_in_file.back(); } void JungfrauDataFile::read_into(std::byte *image_buf, JungfrauDataHeader *header) { - //read header - if(auto rc = fread(header, sizeof(JungfrauDataHeader), 1, m_fp.get()); rc!= 1){ - throw std::runtime_error(LOCATION + - "Could not read header from file:" + m_fp.error_msg()); + // read header if not passed nullptr + if (header) { + if (auto rc = fread(header, sizeof(JungfrauDataHeader), 1, m_fp.get()); + rc != 1) { + throw std::runtime_error( + LOCATION + + "Could not read header from file:" + m_fp.error_msg()); + } + } else { + m_fp.seek(header_size, SEEK_CUR); } - - //read data - if(auto rc = fread(image_buf, 1, m_bytes_per_frame, m_fp.get()); rc != m_bytes_per_frame){ - throw std::runtime_error(LOCATION + - "Could not read image from file" + m_fp.error_msg()); + // read data + if (auto rc = fread(image_buf, 1, m_bytes_per_frame, m_fp.get()); + rc != m_bytes_per_frame) { + throw std::runtime_error(LOCATION + "Could not read image from file" + + m_fp.error_msg()); } - // prepare for next read // if we are at the end of the file, open the next file - ++ m_current_frame; - if(m_current_frame >= m_frame_index[m_current_file_index] && (m_current_frame < m_total_frames)){ + ++m_current_frame_index; + if (m_current_frame_index >= m_last_frame_in_file[m_current_file_index] && + (m_current_frame_index < m_total_frames)) { ++m_current_file_index; open_file(m_current_file_index); - } + } } void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames, JungfrauDataHeader *header) { - - for (size_t i = 0; i < n_frames; ++i) { - read_into(image_buf + i * m_bytes_per_frame, header + i); + if (header) { + for (size_t i = 0; i < n_frames; ++i) + read_into(image_buf + i * m_bytes_per_frame, header + i); + }else{ + for (size_t i = 0; i < n_frames; ++i) + read_into(image_buf + i * m_bytes_per_frame, nullptr); } - } +void JungfrauDataFile::read_into(NDArray* image, JungfrauDataHeader* header) { + if(!(rows() == image->shape(0) && cols() == image->shape(1))){ + throw std::runtime_error(LOCATION + + "Image shape does not match file size: " + std::to_string(rows()) + "x" + std::to_string(cols())); + } + read_into(reinterpret_cast(image->data()), header); +} + +NDArray JungfrauDataFile::read_frame(JungfrauDataHeader* header) { + Shape<2> shape{rows(), cols()}; + NDArray image(shape); + + read_into(reinterpret_cast(image.data()), + header); + + return image; +} + +JungfrauDataHeader JungfrauDataFile::read_header() { + JungfrauDataHeader header; + if (auto rc = fread(&header, 1, sizeof(header), m_fp.get()); + rc != sizeof(header)) { + throw std::runtime_error(LOCATION + "Could not read header from file" + + m_fp.error_msg()); + } + m_fp.seek(-header_size, SEEK_CUR); + return header; +} void JungfrauDataFile::open_file(size_t file_index) { - // fmt::print(stderr, "Opening file: {}\n", fpath(file_index+m_offset).string()); - m_fp = FilePtr(fpath(file_index+m_offset), "rb"); + // fmt::print(stderr, "Opening file: {}\n", + // fpath(file_index+m_offset).string()); + m_fp = FilePtr(fpath(file_index + m_offset), "rb"); m_current_file_index = file_index; } -std::filesystem::path -JungfrauDataFile::fpath(size_t file_index) const{ +std::filesystem::path JungfrauDataFile::fpath(size_t file_index) const { auto fname = fmt::format("{}_{:0{}}.dat", m_base_name, file_index, n_digits_in_file_index); return m_path / fname; diff --git a/src/JungfrauDataFile.test.cpp b/src/JungfrauDataFile.test.cpp index 69f3b72..db5464b 100644 --- a/src/JungfrauDataFile.test.cpp +++ b/src/JungfrauDataFile.test.cpp @@ -25,12 +25,45 @@ TEST_CASE("Open a Jungfrau data file", "[.files]") { REQUIRE(f.total_frames() == 24); REQUIRE(f.current_file() == fpath); + //Check that the frame number and buch id is read correctly for (size_t i = 0; i < 24; ++i) { JungfrauDataHeader header; - auto image = f.read_frame(header); + auto image = f.read_frame(&header); REQUIRE(header.framenum == i + 1); REQUIRE(header.bunchid == (i + 1) * (i + 1)); REQUIRE(image.shape(0) == 512); REQUIRE(image.shape(1) == 1024); } +} + +TEST_CASE("Seek in a JungfrauDataFile", "[.files]"){ + auto fpath = test_data_path() / "dat" / "AldoJF65k_000000.dat"; + REQUIRE(std::filesystem::exists(fpath)); + + JungfrauDataFile f(fpath); + + //The file should have 113 frames + f.seek(19); + REQUIRE(f.tell() == 19); + auto h = f.read_header(); + REQUIRE(h.framenum == 19+1); + + //Reading again does not change the file pointer + auto h2 = f.read_header(); + REQUIRE(h2.framenum == 19+1); + + f.seek(59); + REQUIRE(f.tell() == 59); + auto h3 = f.read_header(); + REQUIRE(h3.framenum == 59+1); + + JungfrauDataHeader h4; + auto image = f.read_frame(&h4); + REQUIRE(h4.framenum == 59+1); + + //now we should be on the next frame + REQUIRE(f.tell() == 60); + REQUIRE(f.read_header().framenum == 60+1); + + REQUIRE_THROWS(f.seek(86356)); //out of range } \ No newline at end of file