Adding support for FileInterface

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

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:

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<size_t> m_frame_index{};
+    std::vector<size_t> 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<uint16_t> &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<std::byte *>(image.data()), &header);
-    }
-    NDArray<uint16_t> read_frame(JungfrauDataHeader& header) {
-        Shape<2> shape{rows(), cols()};
-        NDArray<uint16_t> image(shape);
+    // Extra functions needed for FileInterface
+    Frame read_frame() override;
+    Frame read_frame(size_t frame_number) override;
+    std::vector<Frame> 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<std::byte *>(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<uint16_t>* 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<uint16_t> 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);

python/src/jungfrau_data_file.hpp

@@ -60,8 +60,14 @@ void define_jungfrau_data_file_io_bindings(py::module &m) {
 
     py::class_<JungfrauDataFile>(m, "JungfrauDataFile")
         .def(py::init<const std::filesystem::path &>())
-        .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.
                )")

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<NumpyFile>(fname, mode, cfg);
+    }else if(fname.extension() == ".dat"){
+        file_impl = std::make_unique<JungfrauDataFile>(fname);
     } else {
         throw std::runtime_error("Unsupported file type");
     }

src/JungfrauDataFile.cpp

@@ -1,9 +1,9 @@
 #include "aare/JungfrauDataFile.hpp"
-#include "aare/defs.hpp"
 #include "aare/algorithm.hpp"
+#include "aare/defs.hpp"
 
-#include <fmt/format.h>
 #include <cerrno>
+#include <fmt/format.h>
 
 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<std::byte *>(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<std::byte *>(f.data()), nullptr);
+    return f;
+}
+
+std::vector<Frame> JungfrauDataFile::read_n(size_t n_frames) {
+    std::vector<Frame> 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)
         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 +
@@ -86,51 +133,54 @@ 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 +
+    // 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());
         }
-
-
-    //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());
+    } 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());
+    }
 
     // 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);
     }
@@ -138,22 +188,52 @@ void JungfrauDataFile::read_into(std::byte *image_buf,
 
 void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames,
                                  JungfrauDataHeader *header) {
-
-    for (size_t i = 0; i < n_frames; ++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<uint16_t>* 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<std::byte *>(image->data()), header);
 }
 
+NDArray<uint16_t> JungfrauDataFile::read_frame(JungfrauDataHeader* header) {
+    Shape<2> shape{rows(), cols()};
+    NDArray<uint16_t> image(shape);
+
+    read_into(reinterpret_cast<std::byte *>(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());
+    // 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;

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
+}