Developer (#108)

- Support for very old moench
- read_n in RawFile

conda-recipe/meta.yaml

@@ -1,6 +1,6 @@
 package:
   name: aare
-  version: 2024.11.26.dev0 #TODO! how to not duplicate this?
+  version: 2024.11.27.dev0 #TODO! how to not duplicate this?
 
 
 source:

include/aare/RawMasterFile.hpp

@@ -14,6 +14,7 @@ namespace aare {
  * @brief Implementation used in RawMasterFile to parse the file name
  */
 class RawFileNameComponents {
+    bool m_old_scheme{false};
     std::filesystem::path m_base_path{};
     std::string m_base_name{};
     std::string m_ext{};
@@ -35,6 +36,7 @@ class RawFileNameComponents {
     const std::string &base_name() const;
     const std::string &ext() const;
     int file_index() const;
+    void set_old_scheme(bool old_scheme);
 };
 
 class ScanParameters {
@@ -88,7 +90,7 @@ class RawMasterFile {
     size_t m_pixels_x{};
     size_t m_bitdepth{};
 
-    xy m_geometry;
+    xy m_geometry{};
 
     size_t m_max_frames_per_file{};
     // uint32_t m_adc_mask{}; // TODO! implement reading

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "aare"
-version = "2024.11.26.dev0"
+version = "2024.11.27.dev0"
 
 
 [tool.scikit-build]

python/CMakeLists.txt

@@ -28,8 +28,10 @@ target_link_libraries(_aare PRIVATE aare_core aare_compiler_flags)
 set( PYTHON_FILES
     aare/__init__.py
     aare/CtbRawFile.py
+    aare/RawFile.py
     aare/transform.py
     aare/ScanParameters.py
+    aare/utils.py
 )
 
 # Copy the python files to the build directory

python/aare/CtbRawFile.py

@@ -8,14 +8,15 @@ class CtbRawFile(_aare.CtbRawFile):
     
     Args:  
             fname (pathlib.Path | str): Path to the file to be read.
+            chunk_size (int): Number of frames to read at a time. Default is 1.
             transform (function): Function to apply to the data after reading it. 
                 The function should take a numpy array of type uint8 and return one
                 or several numpy arrays.
     """
-    def __init__(self, fname, transform = None, chunk_size = 1):
+    def __init__(self, fname, chunk_size = 1, transform = None):
         super().__init__(fname)
-        self.transform = transform
         self._chunk_size = chunk_size
+        self._transform = transform
 
 
     def read_frame(self, frame_index: int | None = None ) -> tuple:
@@ -44,8 +45,8 @@ class CtbRawFile(_aare.CtbRawFile):
             header = header[0]
 
 
-        if self.transform:
+        if self._transform:
-            res = self.transform(data)
+            res = self._transform(data)
             if isinstance(res, tuple):
                 return header, *res
             else:

python/aare/RawFile.py

@@ -0,0 +1,66 @@
+from . import _aare
+import numpy as np
+from .ScanParameters import ScanParameters
+
+class RawFile(_aare.RawFile):
+    def __init__(self, fname, chunk_size = 1):
+        super().__init__(fname)
+        self._chunk_size = chunk_size
+
+
+    def read(self) -> tuple:
+        """Read the entire file.
+        Seeks to the beginning of the file before reading.
+
+        Returns:
+            tuple: header, data
+        """
+        self.seek(0)
+        return self.read_n(self.total_frames)
+
+    @property
+    def scan_parameters(self):
+        """Return the scan parameters.
+
+        Returns:
+            ScanParameters: Scan parameters.
+        """
+        return ScanParameters(self.master.scan_parameters)
+    
+    @property
+    def master(self):
+        """Return the master file.
+
+        Returns:
+            RawMasterFile: Master file.
+        """
+        return super().master()
+    
+    def __len__(self) -> int:
+        """Return the number of frames in the file.
+
+        Returns:
+            int: Number of frames in file.
+        """
+        return super().frames_in_file
+
+    def __enter__(self):
+        return self
+    
+    def __exit__(self, exc_type, exc_value, traceback):
+        pass
+
+    def __iter__(self):
+        return self
+    
+    def __next__(self):
+        try:
+            if self._chunk_size == 1:
+                return self.read_frame()
+            else:
+                return self.read_n(self._chunk_size)
+            
+                
+        except RuntimeError:
+            # TODO! find a good way to check that we actually have the right exception
+            raise StopIteration

python/aare/__init__.py

@@ -2,10 +2,13 @@
 from . import _aare
 
 
-from ._aare import File, RawFile, RawMasterFile, RawSubFile
+from ._aare import File, RawMasterFile, RawSubFile
 from ._aare import Pedestal, ClusterFinder, VarClusterFinder
 from ._aare import DetectorType
 from ._aare import ClusterFile
 
 from .CtbRawFile import CtbRawFile
+from .RawFile import RawFile
 from .ScanParameters import ScanParameters
+
+from .utils import random_pixels, random_pixel

python/aare/utils.py

@@ -0,0 +1,23 @@
+import numpy as np
+
+def random_pixels(n_pixels, xmin=0, xmax=512, ymin=0, ymax=1024):
+    """Return a list of random pixels.
+    
+    Args:
+        n_pixels (int): Number of pixels to return.
+        rows (int): Number of rows in the image.
+        cols (int): Number of columns in the image.
+        
+    Returns:
+        list: List of (row, col) tuples.
+    """
+    return [(np.random.randint(ymin, ymax), np.random.randint(xmin, xmax)) for _ in range(n_pixels)]
+
+
+def random_pixel(xmin=0, xmax=512, ymin=0, ymax=1024):
+    """Return a random pixel.
+    
+    Returns:
+        tuple: (row, col)
+    """
+    return random_pixels(1, xmin, xmax, ymin, ymax)[0]

python/src/cluster.hpp

@@ -21,9 +21,9 @@ void define_cluster_finder_bindings(py::module &m) {
              })
         .def("pedestal",
              [](ClusterFinder<uint16_t, double> &self) {
-                 auto m = new NDArray<double, 2>{};
+                 auto pd = new NDArray<double, 2>{};
-                 *m = self.pedestal();
+                 *pd = self.pedestal();
-                 return return_image_data(m);
+                 return return_image_data(pd);
              })
         .def("find_clusters_without_threshold",
              [](ClusterFinder<uint16_t, double> &self,

python/src/ctb_raw_file.hpp

@@ -0,0 +1,57 @@
+
+#include "aare/CtbRawFile.hpp"
+#include "aare/File.hpp"
+#include "aare/Frame.hpp"
+#include "aare/RawFile.hpp"
+#include "aare/RawMasterFile.hpp"
+#include "aare/RawSubFile.hpp"
+
+#include "aare/defs.hpp"
+// #include "aare/fClusterFileV2.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/iostream.h>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl/filesystem.h>
+#include <string>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+void define_ctb_raw_file_io_bindings(py::module &m) {
+
+ py::class_<CtbRawFile>(m, "CtbRawFile")
+        .def(py::init<const std::filesystem::path &>())
+        .def("read_frame",
+             [](CtbRawFile &self) {
+                 size_t image_size = self.image_size_in_bytes();
+                 py::array image;
+                 std::vector<ssize_t> shape;
+                 shape.reserve(2);
+                 shape.push_back(1);
+                 shape.push_back(image_size);
+
+                 py::array_t<DetectorHeader> header(1);
+
+                 // always read bytes
+                 image = py::array_t<uint8_t>(shape);
+
+                 self.read_into(
+                     reinterpret_cast<std::byte *>(image.mutable_data()),
+                     header.mutable_data());
+
+                 return py::make_tuple(header, image);
+             })
+        .def("seek", &CtbRawFile::seek)
+        .def("tell", &CtbRawFile::tell)
+        .def("master", &CtbRawFile::master)
+
+        .def_property_readonly("image_size_in_bytes",
+                               &CtbRawFile::image_size_in_bytes)
+
+        .def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
+
+}

python/src/file.hpp

@@ -38,36 +38,7 @@ void define_file_io_bindings(py::module &m) {
                          bunchId, timestamp, modId, row, column, reserved,
                          debug, roundRNumber, detType, version, packetMask);
 
-    py::class_<CtbRawFile>(m, "CtbRawFile")
-        .def(py::init<const std::filesystem::path &>())
-        .def("read_frame",
-             [](CtbRawFile &self) {
-                 size_t image_size = self.image_size_in_bytes();
-                 py::array image;
-                 std::vector<ssize_t> shape;
-                 shape.reserve(2);
-                 shape.push_back(1);
-                 shape.push_back(image_size);
    
-                 py::array_t<DetectorHeader> header(1);
-
-                 // always read bytes
-                 image = py::array_t<uint8_t>(shape);
-
-                 self.read_into(
-                     reinterpret_cast<std::byte *>(image.mutable_data()),
-                     header.mutable_data());
-
-                 return py::make_tuple(header, image);
-             })
-        .def("seek", &CtbRawFile::seek)
-        .def("tell", &CtbRawFile::tell)
-        .def("master", &CtbRawFile::master)
-
-        .def_property_readonly("image_size_in_bytes",
-                               &CtbRawFile::image_size_in_bytes)
-
-        .def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
 
     py::class_<File>(m, "File")
         .def(py::init([](const std::filesystem::path &fname) {
@@ -133,13 +104,15 @@ void define_file_io_bindings(py::module &m) {
                  return image;
              })
         .def("read_n", [](File &self, size_t n_frames) {
-            const uint8_t item_size = self.bytes_per_pixel();
+            //adjust for actual frames left in the file
+            n_frames = std::min(n_frames, self.total_frames()-self.tell());
+            if(n_frames == 0){
+                throw std::runtime_error("No frames left in file");
+            }
+            std::vector<size_t> shape{n_frames, self.rows(), self.cols()};
+            
             py::array image;
-            std::vector<ssize_t> shape;
+            const uint8_t item_size = self.bytes_per_pixel();
-            shape.reserve(3);
-            shape.push_back(n_frames);
-            shape.push_back(self.rows());
-            shape.push_back(self.cols());
             if (item_size == 1) {
                 image = py::array_t<uint8_t>(shape);
             } else if (item_size == 2) {

python/src/module.cpp

@@ -1,6 +1,7 @@
 //Files with bindings to the different classes
 #include "file.hpp"
 #include "raw_file.hpp"
+#include "ctb_raw_file.hpp"
 #include "raw_master_file.hpp"
 #include "var_cluster.hpp"
 #include "pixel_map.hpp"
@@ -17,6 +18,7 @@ namespace py = pybind11;
 PYBIND11_MODULE(_aare, m) {
     define_file_io_bindings(m);
     define_raw_file_io_bindings(m);
+    define_ctb_raw_file_io_bindings(m);
     define_raw_master_file_bindings(m);
     define_var_cluster_finder_bindings(m);
     define_pixel_map_bindings(m);

python/src/pedestal.hpp

@@ -15,19 +15,19 @@ template <typename SUM_TYPE> void define_pedestal_bindings(py::module &m, const
         .def(py::init<int, int>())
         .def("mean",
              [](Pedestal<SUM_TYPE> &self) {
-                 auto m = new NDArray<SUM_TYPE, 2>{};
+                 auto mea = new NDArray<SUM_TYPE, 2>{};
-                 *m = self.mean();
+                 *mea = self.mean();
-                 return return_image_data(m);
+                 return return_image_data(mea);
              })
         .def("variance", [](Pedestal<SUM_TYPE> &self) {
-            auto m = new NDArray<SUM_TYPE, 2>{};
+            auto var = new NDArray<SUM_TYPE, 2>{};
-            *m = self.variance();
+            *var = self.variance();
-            return return_image_data(m);
+            return return_image_data(var);
         })
         .def("std", [](Pedestal<SUM_TYPE> &self) {
-            auto m = new NDArray<SUM_TYPE, 2>{};
+            auto std = new NDArray<SUM_TYPE, 2>{};
-            *m = self.std();
+            *std = self.std();
-            return return_image_data(m);
+            return return_image_data(std);
         })
         .def("clear", py::overload_cast<>(&Pedestal<SUM_TYPE>::clear))
         .def_property_readonly("rows", &Pedestal<SUM_TYPE>::rows)

python/src/raw_file.hpp

@@ -48,32 +48,44 @@ void define_raw_file_io_bindings(py::module &m) {
 
                  return py::make_tuple(header, image);
              })
-        .def("read_n",
+        .def(
-             [](RawFile &self, size_t n_frames) {
+            "read_n",
-                 py::array image;
+            [](RawFile &self, size_t n_frames) {
-                 std::vector<ssize_t> shape;
+                // adjust for actual frames left in the file
-                 shape.reserve(3);
+                n_frames =
-                 shape.push_back(n_frames);
+                    std::min(n_frames, self.total_frames() - self.tell());
-                 shape.push_back(self.rows());
+                if (n_frames == 0) {
-                 shape.push_back(self.cols());
+                    throw std::runtime_error("No frames left in file");
+                }
+                std::vector<size_t> shape{n_frames, self.rows(), self.cols()};
             
-                 // return headers from all subfiles
+                // return headers from all subfiles
-                 py::array_t<DetectorHeader> header({self.n_mod(), n_frames});
+                py::array_t<DetectorHeader> header;
+                if (self.n_mod() == 1) {
+                    header = py::array_t<DetectorHeader>(n_frames);
+                } else {
+                    header = py::array_t<DetectorHeader>({self.n_mod(), n_frames});
+                }
+                // py::array_t<DetectorHeader> header({self.n_mod(), n_frames});
 
-                 const uint8_t item_size = self.bytes_per_pixel();
+                py::array image;
-                 if (item_size == 1) {
+                const uint8_t item_size = self.bytes_per_pixel();
-                     image = py::array_t<uint8_t>(shape);
+                if (item_size == 1) {
-                 } else if (item_size == 2) {
+                    image = py::array_t<uint8_t>(shape);
-                     image = py::array_t<uint16_t>(shape);
+                } else if (item_size == 2) {
-                 } else if (item_size == 4) {
+                    image = py::array_t<uint16_t>(shape);
-                     image = py::array_t<uint32_t>(shape);
+                } else if (item_size == 4) {
-                 }
+                    image = py::array_t<uint32_t>(shape);
-                 self.read_into(
+                }
-                     reinterpret_cast<std::byte *>(image.mutable_data()),
+                self.read_into(
-                     n_frames, header.mutable_data());
+                    reinterpret_cast<std::byte *>(image.mutable_data()),
+                    n_frames, header.mutable_data());
 
-                 return py::make_tuple(header, image);
+                return py::make_tuple(header, image);
-             })
+            },
+            R"(
+             Read n frames from the file.
+             )")
         .def("frame_number", &RawFile::frame_number)
         .def_property_readonly("bytes_per_frame", &RawFile::bytes_per_frame)
         .def_property_readonly("pixels_per_frame", &RawFile::pixels_per_frame)

src/RawFile.cpp

@@ -17,6 +17,9 @@ RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode)
         n_subfiles = find_number_of_subfiles(); // f0,f1...fn
         n_subfile_parts =
             m_master.geometry().col * m_master.geometry().row; // d0,d1...dn
+
+        
+
         find_geometry();
         update_geometry_with_roi();
 

src/RawMasterFile.cpp

@@ -37,11 +37,24 @@ std::filesystem::path RawFileNameComponents::master_fname() const {
 }
 
 std::filesystem::path RawFileNameComponents::data_fname(size_t mod_id,
-                                                        size_t file_id) const {
+                                                        size_t file_id
-    return m_base_path / fmt::format("{}_d{}_f{}_{}.raw", m_base_name, mod_id,
+                                                        ) const {
+    
+    
+    
+    std::string fmt = "{}_d{}_f{}_{}.raw";
+    //Before version X we used to name the data files f000000000000
+    if (m_old_scheme) {
+        fmt = "{}_d{}_f{:012}_{}.raw";
+    }
+    return m_base_path / fmt::format(fmt, m_base_name, mod_id,
                                      file_id, m_file_index);
 }
 
+void RawFileNameComponents::set_old_scheme(bool old_scheme) {
+    m_old_scheme = old_scheme;
+}
+
 const std::filesystem::path &RawFileNameComponents::base_path() const {
     return m_base_path;
 }
@@ -314,10 +327,22 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) {
             // do the actual parsing
             if (key == "Version") {
                 m_version = value;
+
+                //TODO!: How old versions can we handle?
+                auto v = std::stod(value);
+
+                //TODO! figure out exactly when we did the change
+                //This enables padding of f to 12 digits
+                if (v<4.0)
+                    m_fnc.set_old_scheme(true);
+
             } else if (key == "TimeStamp") {
 
             } else if (key == "Detector Type") {
                 m_type = StringTo<DetectorType>(value);
+                if(m_type==DetectorType::Moench){
+                    m_type = DetectorType::Moench03_old;
+                }
             } else if (key == "Timing Mode") {
                 m_timing_mode = StringTo<TimingMode>(value);
             } else if (key == "Image Size") {
@@ -352,6 +377,12 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) {
                 pos = value.find(',');
                 m_pixels_x = std::stoi(value.substr(1, pos));
                 m_pixels_y = std::stoi(value.substr(pos + 1));
+            }else if(key == "row"){
+                pos = value.find('p');
+                m_pixels_y = std::stoi(value.substr(0, pos));
+            }else if(key == "col"){
+                pos = value.find('p');
+                m_pixels_x = std::stoi(value.substr(0, pos));
             } else if (key == "Total Frames") {
                 m_total_frames_expected = std::stoi(value);
             } else if (key == "Dynamic Range") {
@@ -360,6 +391,9 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) {
                 m_quad = std::stoi(value);
             } else if (key == "Max Frames Per File") {
                 m_max_frames_per_file = std::stoi(value);
+            }else if(key == "Max. Frames Per File"){
+                //Version 3.0 way of writing it
+                m_max_frames_per_file = std::stoi(value);
             } else if (key == "Geometry") {
                 pos = value.find(',');
                 m_geometry = {
@@ -368,5 +402,19 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) {
             }
         }
     }
+    if (m_pixels_x == 400 && m_pixels_y == 400) {
+        m_type = DetectorType::Moench03_old;
+    }
+
+
+    //TODO! Look for d0, d1...dn and update geometry
+    if(m_geometry.col == 0 && m_geometry.row == 0){
+        m_geometry = {1,1};
+        fmt::print("Warning: No geometry found in master file. Assuming 1x1\n");
+    }
+
+    //TODO! Read files and find actual frames
+    if(m_frames_in_file==0)
+        m_frames_in_file = m_total_frames_expected;
 }
 } // namespace aare

src/RawMasterFile.test.cpp

@@ -31,6 +31,16 @@ TEST_CASE("Construction of master file name and data files"){
     REQUIRE(m.data_fname(1, 1) == "test_d1_f1_1.raw");
 }
 
+TEST_CASE("Construction of master file name and data files using old scheme"){
+    RawFileNameComponents m("test_master_1.raw");
+    m.set_old_scheme(true);
+    REQUIRE(m.master_fname() == "test_master_1.raw");
+    REQUIRE(m.data_fname(0, 0) == "test_d0_f000000000000_1.raw");
+    REQUIRE(m.data_fname(1, 0) == "test_d1_f000000000000_1.raw");
+    REQUIRE(m.data_fname(0, 1) == "test_d0_f000000000001_1.raw");
+    REQUIRE(m.data_fname(1, 1) == "test_d1_f000000000001_1.raw");
+}
+
 TEST_CASE("Master file name does not fit pattern"){
     REQUIRE_THROWS(RawFileNameComponents("somefile.json"));
     REQUIRE_THROWS(RawFileNameComponents("another_test_d0_f0_1.raw"));