diff --git a/benchmarks/calculateeta_benchmark.cpp b/benchmarks/calculateeta_benchmark.cpp index a320188..6c40c5c 100644 --- a/benchmarks/calculateeta_benchmark.cpp +++ b/benchmarks/calculateeta_benchmark.cpp @@ -41,8 +41,8 @@ BENCHMARK_F(ClusterFixture, Calculate2x2Eta)(benchmark::State &st) { } // almost takes double the time -BENCHMARK_F(ClusterFixture, - CalculateGeneralEtaFor2x2Cluster)(benchmark::State &st) { +BENCHMARK_F(ClusterFixture, CalculateGeneralEtaFor2x2Cluster) +(benchmark::State &st) { for (auto _ : st) { // This code gets timed Eta2 eta = calculate_eta2(cluster_2x2); @@ -59,8 +59,8 @@ BENCHMARK_F(ClusterFixture, Calculate3x3Eta)(benchmark::State &st) { } // almost takes double the time -BENCHMARK_F(ClusterFixture, - CalculateGeneralEtaFor3x3Cluster)(benchmark::State &st) { +BENCHMARK_F(ClusterFixture, CalculateGeneralEtaFor3x3Cluster) +(benchmark::State &st) { for (auto _ : st) { // This code gets timed Eta2 eta = calculate_eta2(cluster_3x3); diff --git a/benchmarks/ndarray_benchmark.cpp b/benchmarks/ndarray_benchmark.cpp index 55fa263..91a2d9b 100644 --- a/benchmarks/ndarray_benchmark.cpp +++ b/benchmarks/ndarray_benchmark.cpp @@ -1,136 +1,132 @@ -#include #include "aare/NDArray.hpp" - +#include using aare::NDArray; constexpr ssize_t size = 1024; class TwoArrays : public benchmark::Fixture { -public: - NDArray a{{size,size},0}; - NDArray b{{size,size},0}; - void SetUp(::benchmark::State& state) { - for(uint32_t i = 0; i < size; i++){ - for(uint32_t j = 0; j < size; j++){ - a(i, j)= i*j+1; - b(i, j)= i*j+1; - } + public: + NDArray a{{size, size}, 0}; + NDArray b{{size, size}, 0}; + void SetUp(::benchmark::State &state) { + for (uint32_t i = 0; i < size; i++) { + for (uint32_t j = 0; j < size; j++) { + a(i, j) = i * j + 1; + b(i, j) = i * j + 1; + } + } } - } - // void TearDown(::benchmark::State& state) { - // } + // void TearDown(::benchmark::State& state) { + // } }; - - - -BENCHMARK_F(TwoArrays, AddWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a+b; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, AddWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) + b(i); +BENCHMARK_F(TwoArrays, AddWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a + b; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, AddWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) + b(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } -BENCHMARK_F(TwoArrays, SubtractWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a-b; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, SubtractWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) - b(i); +BENCHMARK_F(TwoArrays, SubtractWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a - b; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, SubtractWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) - b(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } -BENCHMARK_F(TwoArrays, MultiplyWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a*b; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, MultiplyWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) * b(i); +BENCHMARK_F(TwoArrays, MultiplyWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a * b; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, MultiplyWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) * b(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } -BENCHMARK_F(TwoArrays, DivideWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a/b; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, DivideWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) / b(i); +BENCHMARK_F(TwoArrays, DivideWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a / b; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, DivideWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) / b(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } -BENCHMARK_F(TwoArrays, FourAddWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a+b+a+b; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, FourAddWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) + b(i) + a(i) + b(i); +BENCHMARK_F(TwoArrays, FourAddWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a + b + a + b; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, FourAddWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) + b(i) + a(i) + b(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } -BENCHMARK_F(TwoArrays, MultiplyAddDivideWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a*a+b/a; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, MultiplyAddDivideWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) * a(i) + b(i) / a(i); +BENCHMARK_F(TwoArrays, MultiplyAddDivideWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a * a + b / a; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, MultiplyAddDivideWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) * a(i) + b(i) / a(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } BENCHMARK_MAIN(); \ No newline at end of file diff --git a/include/aare/ArrayExpr.hpp b/include/aare/ArrayExpr.hpp index d326601..e5fb5d7 100644 --- a/include/aare/ArrayExpr.hpp +++ b/include/aare/ArrayExpr.hpp @@ -1,10 +1,9 @@ #pragma once -#include -#include -#include -#include #include "aare/defs.hpp" - +#include +#include +#include +#include namespace aare { @@ -15,7 +14,9 @@ template class ArrayExpr { auto operator[](size_t i) const { return static_cast(*this)[i]; } auto operator()(size_t i) const { return static_cast(*this)[i]; } auto size() const { return static_cast(*this).size(); } - std::array shape() const { return static_cast(*this).shape(); } + std::array shape() const { + return static_cast(*this).shape(); + } }; template @@ -47,7 +48,7 @@ class ArraySub : public ArrayExpr, Ndim> { }; template -class ArrayMul : public ArrayExpr,Ndim> { +class ArrayMul : public ArrayExpr, Ndim> { const A &arr1_; const B &arr2_; @@ -74,15 +75,13 @@ class ArrayDiv : public ArrayExpr, Ndim> { std::array shape() const { return arr1_.shape(); } }; - - template auto operator+(const ArrayExpr &arr1, const ArrayExpr &arr2) { return ArrayAdd, ArrayExpr, Ndim>(arr1, arr2); } template -auto operator-(const ArrayExpr &arr1, const ArrayExpr &arr2) { +auto operator-(const ArrayExpr &arr1, const ArrayExpr &arr2) { return ArraySub, ArrayExpr, Ndim>(arr1, arr2); } @@ -96,6 +95,4 @@ auto operator/(const ArrayExpr &arr1, const ArrayExpr &arr2) { return ArrayDiv, ArrayExpr, Ndim>(arr1, arr2); } - - } // namespace aare \ No newline at end of file diff --git a/include/aare/CircularFifo.hpp b/include/aare/CircularFifo.hpp index 8098082..853a89f 100644 --- a/include/aare/CircularFifo.hpp +++ b/include/aare/CircularFifo.hpp @@ -17,7 +17,8 @@ template class CircularFifo { public: CircularFifo() : CircularFifo(100){}; - CircularFifo(uint32_t size) : fifo_size(size), free_slots(size + 1), filled_slots(size + 1) { + CircularFifo(uint32_t size) + : fifo_size(size), free_slots(size + 1), filled_slots(size + 1) { // TODO! how do we deal with alignment for writing? alignas??? // Do we give the user a chance to provide memory locations? @@ -55,7 +56,8 @@ template class CircularFifo { bool try_pop_free(ItemType &v) { return free_slots.read(v); } - ItemType pop_value(std::chrono::nanoseconds wait, std::atomic &stopped) { + ItemType pop_value(std::chrono::nanoseconds wait, + std::atomic &stopped) { ItemType v; while (!filled_slots.read(v) && !stopped) { std::this_thread::sleep_for(wait); diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index e26e765..c0eca33 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -371,8 +371,8 @@ ClusterFile::read_frame_without_cut() { "Could not read number of clusters"); } - LOG(logDEBUG1) << "Reading " << n_clusters - << " clusters from frame " << frame_number; + LOG(logDEBUG1) << "Reading " << n_clusters << " clusters from frame " + << frame_number; ClusterVector clusters(n_clusters); clusters.set_frame_number(frame_number); diff --git a/include/aare/ClusterFileSink.hpp b/include/aare/ClusterFileSink.hpp index 09190fe..1900774 100644 --- a/include/aare/ClusterFileSink.hpp +++ b/include/aare/ClusterFileSink.hpp @@ -49,8 +49,8 @@ class ClusterFileSink { ClusterFileSink(ClusterFinderMT *source, const std::filesystem::path &fname) { LOG(logDEBUG) << "ClusterFileSink: " - << "source: " << source->sink() - << ", file: " << fname.string(); + << "source: " << source->sink() + << ", file: " << fname.string(); m_source = source->sink(); m_thread = std::thread(&ClusterFileSink::process, this); m_file.open(fname, std::ios::binary); diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index 7a34722..069d887 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -39,10 +39,9 @@ class ClusterFinder { c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)), c3(sqrt(ClusterSizeX * ClusterSizeY)), m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) { - LOG(logDEBUG ) << "ClusterFinder: " - << "image_size: " << image_size[0] << "x" << image_size[1] - << ", nSigma: " << nSigma - << ", capacity: " << capacity; + LOG(logDEBUG) << "ClusterFinder: " + << "image_size: " << image_size[0] << "x" << image_size[1] + << ", nSigma: " << nSigma << ", capacity: " << capacity; } void push_pedestal_frame(NDView frame) { diff --git a/include/aare/ClusterFinderMT.hpp b/include/aare/ClusterFinderMT.hpp index efc22d4..0340973 100644 --- a/include/aare/ClusterFinderMT.hpp +++ b/include/aare/ClusterFinderMT.hpp @@ -7,8 +7,8 @@ #include "aare/ClusterFinder.hpp" #include "aare/NDArray.hpp" -#include "aare/logger.hpp" #include "aare/ProducerConsumerQueue.hpp" +#include "aare/logger.hpp" namespace aare { @@ -125,10 +125,10 @@ class ClusterFinderMT { : m_n_threads(n_threads) { LOG(logDEBUG1) << "ClusterFinderMT: " - << "image_size: " << image_size[0] << "x" << image_size[1] - << ", nSigma: " << nSigma - << ", capacity: " << capacity - << ", n_threads: " << n_threads; + << "image_size: " << image_size[0] << "x" + << image_size[1] << ", nSigma: " << nSigma + << ", capacity: " << capacity + << ", n_threads: " << n_threads; for (size_t i = 0; i < n_threads; i++) { m_cluster_finders.push_back( diff --git a/include/aare/CtbRawFile.hpp b/include/aare/CtbRawFile.hpp index 68dab23..afae0a2 100644 --- a/include/aare/CtbRawFile.hpp +++ b/include/aare/CtbRawFile.hpp @@ -1,27 +1,27 @@ #pragma once #include "aare/FileInterface.hpp" -#include "aare/RawMasterFile.hpp" #include "aare/Frame.hpp" +#include "aare/RawMasterFile.hpp" #include #include -namespace aare{ +namespace aare { - -class CtbRawFile{ +class CtbRawFile { RawMasterFile m_master; std::ifstream m_file; size_t m_current_frame{0}; size_t m_current_subfile{0}; size_t m_num_subfiles{0}; -public: + + public: CtbRawFile(const std::filesystem::path &fname); - void read_into(std::byte *image_buf, DetectorHeader* header = nullptr); - void seek(size_t frame_index); //!< seek to the given frame index - size_t tell() const; //!< get the frame index of the file pointer + void read_into(std::byte *image_buf, DetectorHeader *header = nullptr); + void seek(size_t frame_index); //!< seek to the given frame index + size_t tell() const; //!< get the frame index of the file pointer // in the specific class we can expose more functionality @@ -29,13 +29,13 @@ public: size_t frames_in_file() const; RawMasterFile master() const; -private: + + private: void find_subfiles(); size_t sub_file_index(size_t frame_index) const { return frame_index / m_master.max_frames_per_file(); } void open_data_file(size_t subfile_index); - }; -} \ No newline at end of file +} // namespace aare \ No newline at end of file diff --git a/include/aare/Dtype.hpp b/include/aare/Dtype.hpp index 7e1e62a..7047264 100644 --- a/include/aare/Dtype.hpp +++ b/include/aare/Dtype.hpp @@ -6,31 +6,37 @@ namespace aare { -// The format descriptor is a single character that specifies the type of the data +// The format descriptor is a single character that specifies the type of the +// data // - python documentation: https://docs.python.org/3/c-api/arg.html#numbers -// - py::format_descriptor::format() (in pybind11) does not return the same format as +// - py::format_descriptor::format() (in pybind11) does not return the same +// format as // written in python.org documentation. -// - numpy also doesn't use the same format. and also numpy associates the format -// with variable bitdepth types. (e.g. long is int64 on linux64 and int32 on win64) -// https://numpy.org/doc/stable/reference/arrays.scalars.html +// - numpy also doesn't use the same format. and also numpy associates the +// format +// with variable bitdepth types. (e.g. long is int64 on linux64 and int32 on +// win64) https://numpy.org/doc/stable/reference/arrays.scalars.html // // github issue discussing this: // https://github.com/pybind/pybind11/issues/1908#issuecomment-658358767 // -// [IN LINUX] the difference is for int64 (long) and uint64 (unsigned long). The format -// descriptor is 'q' and 'Q' respectively and in the documentation it is 'l' and 'k'. +// [IN LINUX] the difference is for int64 (long) and uint64 (unsigned long). The +// format descriptor is 'q' and 'Q' respectively and in the documentation it is +// 'l' and 'k'. // in practice numpy doesn't seem to care when reading buffer info: the library // interprets 'q' or 'l' as int64 and 'Q' or 'L' as uint64. -// for this reason we decided to use the same format descriptor as pybind to avoid -// any further discrepancies. +// for this reason we decided to use the same format descriptor as pybind to +// avoid any further discrepancies. // in the following order: // int8, uint8, int16, uint16, int32, uint32, int64, uint64, float, double -const char DTYPE_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', 'I', 'q', 'Q', 'f', 'd'}; +const char DTYPE_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', + 'I', 'q', 'Q', 'f', 'd'}; // on linux64 & apple -const char NUMPY_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', 'I', 'l', 'L', 'f', 'd'}; +const char NUMPY_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', + 'I', 'l', 'L', 'f', 'd'}; /** * @brief enum class to define the endianess of the system */ @@ -52,12 +58,29 @@ enum class endian { */ class Dtype { public: - enum TypeIndex { INT8, UINT8, INT16, UINT16, INT32, UINT32, INT64, UINT64, FLOAT, DOUBLE, ERROR, NONE }; + enum TypeIndex { + INT8, + UINT8, + INT16, + UINT16, + INT32, + UINT32, + INT64, + UINT64, + FLOAT, + DOUBLE, + ERROR, + NONE + }; uint8_t bitdepth() const; size_t bytes() const; - std::string format_descr() const { return std::string(1, DTYPE_FORMAT_DSC[static_cast(m_type)]); } - std::string numpy_descr() const { return std::string(1, NUMPY_FORMAT_DSC[static_cast(m_type)]); } + std::string format_descr() const { + return std::string(1, DTYPE_FORMAT_DSC[static_cast(m_type)]); + } + std::string numpy_descr() const { + return std::string(1, NUMPY_FORMAT_DSC[static_cast(m_type)]); + } explicit Dtype(const std::type_info &t); explicit Dtype(std::string_view sv); diff --git a/include/aare/File.hpp b/include/aare/File.hpp index 1cef898..e8f1589 100644 --- a/include/aare/File.hpp +++ b/include/aare/File.hpp @@ -5,12 +5,12 @@ namespace aare { /** - * @brief RAII File class for reading, and in the future potentially writing - * image files in various formats. Minimal generic interface. For specail fuctions - * plase use the RawFile or NumpyFile classes directly. - * Wraps FileInterface to abstract the underlying file format - * @note **frame_number** refers the the frame number sent by the detector while **frame_index** - * is the position of the frame in the file + * @brief RAII File class for reading, and in the future potentially writing + * image files in various formats. Minimal generic interface. For specail + * fuctions plase use the RawFile or NumpyFile classes directly. Wraps + * FileInterface to abstract the underlying file format + * @note **frame_number** refers the the frame number sent by the detector while + * **frame_index** is the position of the frame in the file */ class File { std::unique_ptr file_impl; @@ -25,42 +25,46 @@ class File { * @throws std::invalid_argument if the file mode is not supported * */ - File(const std::filesystem::path &fname, const std::string &mode="r", const FileConfig &cfg = {}); - - /**Since the object is responsible for managing the file we disable copy construction */ - File(File const &other) = delete; + File(const std::filesystem::path &fname, const std::string &mode = "r", + const FileConfig &cfg = {}); + + /**Since the object is responsible for managing the file we disable copy + * construction */ + File(File const &other) = delete; /**The same goes for copy assignment */ - File& operator=(File const &other) = delete; + File &operator=(File const &other) = delete; File(File &&other) noexcept; - File& operator=(File &&other) noexcept; + File &operator=(File &&other) noexcept; ~File() = default; // void close(); //!< close the file - - Frame read_frame(); //!< read one frame from the file at the current position - Frame read_frame(size_t frame_index); //!< read one frame at the position given by frame number - std::vector read_n(size_t n_frames); //!< read n_frames from the file at the current position + + Frame + read_frame(); //!< read one frame from the file at the current position + Frame read_frame(size_t frame_index); //!< read one frame at the position + //!< given by frame number + std::vector read_n(size_t n_frames); //!< read n_frames from the file + //!< at the current position void read_into(std::byte *image_buf); void read_into(std::byte *image_buf, size_t n_frames); - - size_t frame_number(); //!< get the frame number at the current position - size_t frame_number(size_t frame_index); //!< get the frame number at the given frame index - size_t bytes_per_frame() const; - size_t pixels_per_frame() const; - size_t bytes_per_pixel() const; + + size_t frame_number(); //!< get the frame number at the current position + size_t frame_number( + size_t frame_index); //!< get the frame number at the given frame index + size_t bytes_per_frame() const; + size_t pixels_per_frame() const; + 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 + 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; DetectorType detector_type() const; - - }; } // namespace aare \ No newline at end of file diff --git a/include/aare/FileInterface.hpp b/include/aare/FileInterface.hpp index 3736c46..6ca4755 100644 --- a/include/aare/FileInterface.hpp +++ b/include/aare/FileInterface.hpp @@ -20,8 +20,10 @@ struct FileConfig { uint64_t rows{}; uint64_t cols{}; bool operator==(const FileConfig &other) const { - return dtype == other.dtype && rows == other.rows && cols == other.cols && geometry == other.geometry && - detector_type == other.detector_type && max_frames_per_file == other.max_frames_per_file; + return dtype == other.dtype && rows == other.rows && + cols == other.cols && geometry == other.geometry && + detector_type == other.detector_type && + max_frames_per_file == other.max_frames_per_file; } bool operator!=(const FileConfig &other) const { return !(*this == other); } @@ -32,8 +34,11 @@ struct FileConfig { int max_frames_per_file{}; size_t total_frames{}; std::string to_string() const { - return "{ dtype: " + dtype.to_string() + ", rows: " + std::to_string(rows) + ", cols: " + std::to_string(cols) + - ", geometry: " + geometry.to_string() + ", detector_type: " + ToString(detector_type) + + return "{ dtype: " + dtype.to_string() + + ", rows: " + std::to_string(rows) + + ", cols: " + std::to_string(cols) + + ", geometry: " + geometry.to_string() + + ", detector_type: " + ToString(detector_type) + ", max_frames_per_file: " + std::to_string(max_frames_per_file) + ", total_frames: " + std::to_string(total_frames) + " }"; } @@ -42,7 +47,8 @@ struct FileConfig { /** * @brief FileInterface class to define the interface for file operations * @note parent class for NumpyFile and RawFile - * @note all functions are pure virtual and must be implemented by the derived classes + * @note all functions are pure virtual and must be implemented by the derived + * classes */ class FileInterface { public: @@ -64,17 +70,20 @@ class FileInterface { * @param n_frames number of frames to read * @return vector of frames */ - virtual std::vector read_n(size_t n_frames) = 0; // Is this the right interface? + virtual std::vector + read_n(size_t n_frames) = 0; // Is this the right interface? /** - * @brief read one frame from the file at the current position and store it in the provided buffer + * @brief read one frame from the file at the current position and store it + * in the provided buffer * @param image_buf buffer to store the frame * @return void */ virtual void read_into(std::byte *image_buf) = 0; /** - * @brief read n_frames from the file at the current position and store them in the provided buffer + * @brief read n_frames from the file at the current position and store them + * in the provided buffer * @param image_buf buffer to store the frames * @param n_frames number of frames to read * @return void @@ -134,7 +143,6 @@ class FileInterface { */ virtual size_t bitdepth() const = 0; - virtual DetectorType detector_type() const = 0; // function to query the data type of the file diff --git a/include/aare/FilePtr.hpp b/include/aare/FilePtr.hpp index 4ddc76e..2ffc293 100644 --- a/include/aare/FilePtr.hpp +++ b/include/aare/FilePtr.hpp @@ -12,7 +12,7 @@ class FilePtr { public: FilePtr() = default; - FilePtr(const std::filesystem::path& fname, const std::string& mode); + FilePtr(const std::filesystem::path &fname, const std::string &mode); FilePtr(const FilePtr &) = delete; // we don't want a copy FilePtr &operator=(const FilePtr &) = delete; // since we handle a resource FilePtr(FilePtr &&other); diff --git a/include/aare/Fit.hpp b/include/aare/Fit.hpp index eb9ac22..1beec0a 100644 --- a/include/aare/Fit.hpp +++ b/include/aare/Fit.hpp @@ -23,16 +23,19 @@ NDArray scurve2(NDView x, NDView par); } // namespace func - /** * @brief Estimate the initial parameters for a Gaussian fit */ -std::array gaus_init_par(const NDView x, const NDView y); +std::array gaus_init_par(const NDView x, + const NDView y); -std::array pol1_init_par(const NDView x, const NDView y); +std::array pol1_init_par(const NDView x, + const NDView y); -std::array scurve_init_par(const NDView x, const NDView y); -std::array scurve2_init_par(const NDView x, const NDView y); +std::array scurve_init_par(const NDView x, + const NDView y); +std::array scurve2_init_par(const NDView x, + const NDView y); static constexpr int DEFAULT_NUM_THREADS = 4; @@ -43,7 +46,6 @@ static constexpr int DEFAULT_NUM_THREADS = 4; */ NDArray fit_gaus(NDView x, NDView y); - /** * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values] * @param x x values @@ -54,9 +56,6 @@ NDArray fit_gaus(NDView x, NDView y); NDArray fit_gaus(NDView x, NDView y, int n_threads = DEFAULT_NUM_THREADS); - - - /** * @brief Fit a 1D Gaussian with error estimates * @param x x values @@ -67,7 +66,7 @@ NDArray fit_gaus(NDView x, NDView y, */ void fit_gaus(NDView x, NDView y, NDView y_err, NDView par_out, NDView par_err_out, - double& chi2); + double &chi2); /** * @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout @@ -80,9 +79,8 @@ void fit_gaus(NDView x, NDView y, NDView y_err, * @param n_threads number of threads to use */ void fit_gaus(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads = DEFAULT_NUM_THREADS - ); + NDView par_out, NDView par_err_out, + NDView chi2_out, int n_threads = DEFAULT_NUM_THREADS); NDArray fit_pol1(NDView x, NDView y); @@ -90,26 +88,33 @@ NDArray fit_pol1(NDView x, NDView y, int n_threads = DEFAULT_NUM_THREADS); void fit_pol1(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2); + NDView par_out, NDView par_err_out, + double &chi2); // TODO! not sure we need to offer the different version in C++ void fit_pol1(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out,NDView chi2_out, - int n_threads = DEFAULT_NUM_THREADS); + NDView par_out, NDView par_err_out, + NDView chi2_out, int n_threads = DEFAULT_NUM_THREADS); -NDArray fit_scurve(NDView x, NDView y); -NDArray fit_scurve(NDView x, NDView y, int n_threads); -void fit_scurve(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2); -void fit_scurve(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads); +NDArray fit_scurve(NDView x, NDView y); +NDArray fit_scurve(NDView x, NDView y, + int n_threads); +void fit_scurve(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, double &chi2); +void fit_scurve(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, NDView chi2_out, + int n_threads); -NDArray fit_scurve2(NDView x, NDView y); -NDArray fit_scurve2(NDView x, NDView y, int n_threads); -void fit_scurve2(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2); -void fit_scurve2(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads); +NDArray fit_scurve2(NDView x, NDView y); +NDArray fit_scurve2(NDView x, NDView y, + int n_threads); +void fit_scurve2(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, double &chi2); +void fit_scurve2(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, NDView chi2_out, + int n_threads); } // namespace aare \ No newline at end of file diff --git a/include/aare/Frame.hpp b/include/aare/Frame.hpp index 02ea82f..27a2a4a 100644 --- a/include/aare/Frame.hpp +++ b/include/aare/Frame.hpp @@ -19,7 +19,7 @@ class Frame { uint32_t m_cols; Dtype m_dtype; std::byte *m_data; - //TODO! Add frame number? + // TODO! Add frame number? public: /** @@ -39,7 +39,7 @@ class Frame { * @param dtype data type of the pixels */ Frame(const std::byte *bytes, uint32_t rows, uint32_t cols, Dtype dtype); - ~Frame(){ delete[] m_data; }; + ~Frame() { delete[] m_data; }; /** @warning Copy is disabled to ensure performance when passing * frames around. Can discuss enabling it. @@ -52,7 +52,6 @@ class Frame { Frame &operator=(Frame &&other) noexcept; Frame(Frame &&other) noexcept; - Frame clone() const; //<- Explicit copy uint32_t rows() const; @@ -93,7 +92,7 @@ class Frame { if (row >= m_rows || col >= m_cols) { throw std::out_of_range("Invalid row or column index"); } - //TODO! add tests then reimplement using pixel_ptr + // TODO! add tests then reimplement using pixel_ptr T data; std::memcpy(&data, m_data + (row * m_cols + col) * m_dtype.bytes(), m_dtype.bytes()); @@ -102,9 +101,9 @@ class Frame { /** * @brief Return an NDView of the frame. This is the preferred way to access * data in the frame. - * + * * @tparam T type of the pixels - * @return NDView + * @return NDView */ template NDView view() { std::array shape = {static_cast(m_rows), @@ -113,7 +112,7 @@ class Frame { return NDView(data, shape); } - /** + /** * @brief Copy the frame data into a new NDArray. This is a deep copy. */ template NDArray image() { diff --git a/include/aare/JungfrauDataFile.hpp b/include/aare/JungfrauDataFile.hpp index 9b1bc48..f871b86 100644 --- a/include/aare/JungfrauDataFile.hpp +++ b/include/aare/JungfrauDataFile.hpp @@ -3,104 +3,113 @@ #include #include -#include "aare/FilePtr.hpp" -#include "aare/defs.hpp" -#include "aare/NDArray.hpp" #include "aare/FileInterface.hpp" +#include "aare/FilePtr.hpp" +#include "aare/NDArray.hpp" +#include "aare/defs.hpp" namespace aare { - -struct JungfrauDataHeader{ +struct JungfrauDataHeader { uint64_t framenum; uint64_t bunchid; }; class JungfrauDataFile : public FileInterface { - 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_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) + 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_last_frame_in_file{}; //!< Used for seeking to the correct file + 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 - FilePtr m_fp; //!< RAII wrapper for a FILE* + FilePtr m_fp; //!< RAII wrapper for a FILE* - using pixel_type = uint16_t; - static constexpr size_t header_size = sizeof(JungfrauDataHeader); - static constexpr size_t n_digits_in_file_index = 6; //!< to format file names + static constexpr size_t header_size = sizeof(JungfrauDataHeader); + 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() override; - size_t pixels_per_frame() override; - size_t bytes_per_pixel() const; + std::string base_name() + const; //!< get the base name of the file (without path and extension) + size_t bytes_per_frame() override; + size_t pixels_per_frame() override; + size_t bytes_per_pixel() const; 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 + 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; - std::array shape() const; - size_t n_files() const; //!< get the number of files in the series. + std::array shape() const; + size_t n_files() const; //!< get the number of files in the series. // 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; + 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; /** - * @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) + * @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); /** - * @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) + * @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); - - /** + 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); + void read_into(NDArray *image, + JungfrauDataHeader *header = nullptr); JungfrauDataHeader read_header(); - std::filesystem::path current_file() const { return fpath(m_current_file_index+m_offset); } + std::filesystem::path current_file() const { + return fpath(m_current_file_index + m_offset); + } - - private: + private: /** - * @brief Find the size of the frame in the file. (256x256, 256x1024, 512x1024) + * @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 + * @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 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); - std::filesystem::path fpath(size_t frame_index) const; - - - }; + void parse_fname(const std::filesystem::path &fname); + void scan_files(); + void open_file(size_t file_index); + std::filesystem::path fpath(size_t frame_index) const; +}; } // namespace aare \ No newline at end of file diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp index 3c08a3c..1a501eb 100644 --- a/include/aare/NDArray.hpp +++ b/include/aare/NDArray.hpp @@ -21,7 +21,6 @@ TODO! Add expression templates for operators namespace aare { - template class NDArray : public ArrayExpr, Ndim> { std::array shape_; @@ -34,7 +33,7 @@ class NDArray : public ArrayExpr, Ndim> { * @brief Default constructor. Will construct an empty NDArray. * */ - NDArray() : shape_(), strides_(c_strides(shape_)), data_(nullptr) {}; + NDArray() : shape_(), strides_(c_strides(shape_)), data_(nullptr){}; /** * @brief Construct a new NDArray object with a given shape. @@ -48,7 +47,6 @@ class NDArray : public ArrayExpr, Ndim> { std::multiplies<>())), data_(new T[size_]) {} - /** * @brief Construct a new NDArray object with a shape and value. * @@ -69,8 +67,8 @@ class NDArray : public ArrayExpr, Ndim> { std::copy(v.begin(), v.end(), begin()); } - template - NDArray(const std::array& arr) : NDArray({Size}) { + template + NDArray(const std::array &arr) : NDArray({Size}) { std::copy(arr.begin(), arr.end(), begin()); } @@ -79,7 +77,6 @@ class NDArray : public ArrayExpr, Ndim> { : shape_(other.shape_), strides_(c_strides(shape_)), size_(other.size_), data_(other.data_) { other.reset(); // TODO! is this necessary? - } // Copy constructor @@ -113,10 +110,10 @@ class NDArray : public ArrayExpr, Ndim> { NDArray &operator-=(const NDArray &other); NDArray &operator*=(const NDArray &other); - //Write directly to the data array, or create a new one - template - NDArray& operator=(const std::array &other){ - if(Size != size_){ + // Write directly to the data array, or create a new one + template + NDArray &operator=(const std::array &other) { + if (Size != size_) { delete[] data_; size_ = Size; data_ = new T[size_]; @@ -157,11 +154,6 @@ class NDArray : public ArrayExpr, Ndim> { NDArray &operator&=(const T & /*mask*/); - - - - - void sqrt() { for (int i = 0; i < size_; ++i) { data_[i] = std::sqrt(data_[i]); @@ -345,9 +337,6 @@ NDArray &NDArray::operator+=(const T &value) { return *this; } - - - template NDArray NDArray::operator+(const T &value) { NDArray result = *this; @@ -448,6 +437,4 @@ NDArray load(const std::string &pathname, return img; } - - } // namespace aare \ No newline at end of file diff --git a/include/aare/NDView.hpp b/include/aare/NDView.hpp index 56054e2..e7ad002 100644 --- a/include/aare/NDView.hpp +++ b/include/aare/NDView.hpp @@ -1,6 +1,6 @@ #pragma once -#include "aare/defs.hpp" #include "aare/ArrayExpr.hpp" +#include "aare/defs.hpp" #include #include @@ -17,7 +17,8 @@ namespace aare { template using Shape = std::array; // TODO! fix mismatch between signed and unsigned -template Shape make_shape(const std::vector &shape) { +template +Shape make_shape(const std::vector &shape) { if (shape.size() != Ndim) throw std::runtime_error("Shape size mismatch"); Shape arr; @@ -25,14 +26,18 @@ template Shape make_shape(const std::vector &shape) return arr; } -template ssize_t element_offset(const Strides & /*unused*/) { return 0; } +template +ssize_t element_offset(const Strides & /*unused*/) { + return 0; +} template ssize_t element_offset(const Strides &strides, ssize_t i, Ix... index) { return i * strides[Dim] + element_offset(strides, index...); } -template std::array c_strides(const std::array &shape) { +template +std::array c_strides(const std::array &shape) { std::array strides{}; std::fill(strides.begin(), strides.end(), 1); for (ssize_t i = Ndim - 1; i > 0; --i) { @@ -41,14 +46,16 @@ template std::array c_strides(const std::array std::array make_array(const std::vector &vec) { +template +std::array make_array(const std::vector &vec) { assert(vec.size() == Ndim); std::array arr{}; std::copy_n(vec.begin(), Ndim, arr.begin()); return arr; } -template class NDView : public ArrayExpr, Ndim> { +template +class NDView : public ArrayExpr, Ndim> { public: NDView() = default; ~NDView() = default; @@ -57,17 +64,23 @@ template class NDView : public ArrayExpr shape) : buffer_(buffer), strides_(c_strides(shape)), shape_(shape), - size_(std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>())) {} + size_(std::accumulate(std::begin(shape), std::end(shape), 1, + std::multiplies<>())) {} // NDView(T *buffer, const std::vector &shape) - // : buffer_(buffer), strides_(c_strides(make_array(shape))), shape_(make_array(shape)), - // size_(std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>())) {} + // : buffer_(buffer), + // strides_(c_strides(make_array(shape))), + // shape_(make_array(shape)), + // size_(std::accumulate(std::begin(shape), std::end(shape), 1, + // std::multiplies<>())) {} - template std::enable_if_t operator()(Ix... index) { + template + std::enable_if_t operator()(Ix... index) { return buffer_[element_offset(strides_, index...)]; } - template std::enable_if_t operator()(Ix... index) const { + template + std::enable_if_t operator()(Ix... index) const { return buffer_[element_offset(strides_, index...)]; } @@ -94,16 +107,21 @@ template class NDView : public ArrayExpr()); } NDView &operator-=(const T val) { return elemenwise(val, std::minus()); } - NDView &operator*=(const T val) { return elemenwise(val, std::multiplies()); } - NDView &operator/=(const T val) { return elemenwise(val, std::divides()); } + NDView &operator*=(const T val) { + return elemenwise(val, std::multiplies()); + } + NDView &operator/=(const T val) { + return elemenwise(val, std::divides()); + } - NDView &operator/=(const NDView &other) { return elemenwise(other, std::divides()); } + NDView &operator/=(const NDView &other) { + return elemenwise(other, std::divides()); + } - - template - NDView& operator=(const std::array &arr) { - if(size() != static_cast(arr.size())) - throw std::runtime_error(LOCATION + "Array and NDView size mismatch"); + template NDView &operator=(const std::array &arr) { + if (size() != static_cast(arr.size())) + throw std::runtime_error(LOCATION + + "Array and NDView size mismatch"); std::copy(arr.begin(), arr.end(), begin()); return *this; } @@ -147,13 +165,15 @@ template class NDView : public ArrayExpr shape_{}; uint64_t size_{}; - template NDView &elemenwise(T val, BinaryOperation op) { + template + NDView &elemenwise(T val, BinaryOperation op) { for (uint64_t i = 0; i != size_; ++i) { buffer_[i] = op(buffer_[i], val); } return *this; } - template NDView &elemenwise(const NDView &other, BinaryOperation op) { + template + NDView &elemenwise(const NDView &other, BinaryOperation op) { for (uint64_t i = 0; i != size_; ++i) { buffer_[i] = op(buffer_[i], other.buffer_[i]); } @@ -170,9 +190,8 @@ template void NDView::print_all() const { } } - template -std::ostream& operator <<(std::ostream& os, const NDView& arr){ +std::ostream &operator<<(std::ostream &os, const NDView &arr) { for (auto row = 0; row < arr.shape(0); ++row) { for (auto col = 0; col < arr.shape(1); ++col) { os << std::setw(3); @@ -183,10 +202,8 @@ std::ostream& operator <<(std::ostream& os, const NDView& arr){ return os; } - -template -NDView make_view(std::vector& vec){ - return NDView(vec.data(), {static_cast(vec.size())}); +template NDView make_view(std::vector &vec) { + return NDView(vec.data(), {static_cast(vec.size())}); } } // namespace aare \ No newline at end of file diff --git a/include/aare/NumpyFile.hpp b/include/aare/NumpyFile.hpp index 7381a76..481a1a0 100644 --- a/include/aare/NumpyFile.hpp +++ b/include/aare/NumpyFile.hpp @@ -1,9 +1,8 @@ #pragma once #include "aare/Dtype.hpp" -#include "aare/defs.hpp" #include "aare/FileInterface.hpp" #include "aare/NumpyHelpers.hpp" - +#include "aare/defs.hpp" #include #include @@ -11,13 +10,12 @@ namespace aare { - - /** * @brief NumpyFile class to read and write numpy files * @note derived from FileInterface * @note implements all the pure virtual functions from FileInterface - * @note documentation for the functions can also be found in the FileInterface class + * @note documentation for the functions can also be found in the FileInterface + * class */ class NumpyFile : public FileInterface { @@ -28,26 +26,35 @@ class NumpyFile : public FileInterface { * @param mode file mode (r, w) * @param cfg file configuration */ - explicit NumpyFile(const std::filesystem::path &fname, const std::string &mode = "r", FileConfig cfg = {}); + explicit NumpyFile(const std::filesystem::path &fname, + const std::string &mode = "r", FileConfig cfg = {}); void write(Frame &frame); Frame read_frame() override { return get_frame(this->current_frame++); } - Frame read_frame(size_t frame_number) override { return get_frame(frame_number); } + Frame read_frame(size_t frame_number) override { + return get_frame(frame_number); + } std::vector read_n(size_t n_frames) override; - void read_into(std::byte *image_buf) override { return get_frame_into(this->current_frame++, image_buf); } + void read_into(std::byte *image_buf) override { + return get_frame_into(this->current_frame++, image_buf); + } void read_into(std::byte *image_buf, size_t n_frames) override; size_t frame_number(size_t frame_index) override { return frame_index; }; size_t bytes_per_frame() override; size_t pixels_per_frame() override; - void seek(size_t frame_number) override { this->current_frame = frame_number; } + void seek(size_t frame_number) override { + this->current_frame = frame_number; + } size_t tell() override { return this->current_frame; } size_t total_frames() const override { return m_header.shape[0]; } size_t rows() const override { return m_header.shape[1]; } size_t cols() const override { return m_header.shape[2]; } size_t bitdepth() const override { return m_header.dtype.bitdepth(); } - DetectorType detector_type() const override { return DetectorType::Unknown; } + DetectorType detector_type() const override { + return DetectorType::Unknown; + } /** * @brief get the data type of the numpy file @@ -70,7 +77,8 @@ class NumpyFile : public FileInterface { template NDArray load() { NDArray arr(make_shape(m_header.shape)); if (fseek(fp, static_cast(header_size), SEEK_SET)) { - throw std::runtime_error(LOCATION + "Error seeking to the start of the data"); + throw std::runtime_error(LOCATION + + "Error seeking to the start of the data"); } size_t rc = fread(arr.data(), sizeof(T), arr.size(), fp); if (rc != static_cast(arr.size())) { @@ -78,16 +86,20 @@ class NumpyFile : public FileInterface { } return arr; } - template void write(NDView &frame) { + template + void write(NDView &frame) { write_impl(frame.data(), frame.total_bytes()); } - template void write(NDArray &frame) { + template + void write(NDArray &frame) { write_impl(frame.data(), frame.total_bytes()); } - template void write(NDView &&frame) { + template + void write(NDView &&frame) { write_impl(frame.data(), frame.total_bytes()); } - template void write(NDArray &&frame) { + template + void write(NDArray &&frame) { write_impl(frame.data(), frame.total_bytes()); } diff --git a/include/aare/NumpyHelpers.hpp b/include/aare/NumpyHelpers.hpp index 8ed0ec7..2facc4c 100644 --- a/include/aare/NumpyHelpers.hpp +++ b/include/aare/NumpyHelpers.hpp @@ -40,15 +40,18 @@ bool parse_bool(const std::string &in); std::string get_value_from_map(const std::string &mapstr); -std::unordered_map parse_dict(std::string in, const std::vector &keys); +std::unordered_map +parse_dict(std::string in, const std::vector &keys); -template bool in_array(T val, const std::array &arr) { +template +bool in_array(T val, const std::array &arr) { return std::find(std::begin(arr), std::end(arr), val) != std::end(arr); } bool is_digits(const std::string &str); aare::Dtype parse_descr(std::string typestring); -size_t write_header(const std::filesystem::path &fname, const NumpyHeader &header); +size_t write_header(const std::filesystem::path &fname, + const NumpyHeader &header); size_t write_header(std::ostream &out, const NumpyHeader &header); } // namespace NumpyHelpers diff --git a/include/aare/Pedestal.hpp b/include/aare/Pedestal.hpp index d6223c1..0efc4b7 100644 --- a/include/aare/Pedestal.hpp +++ b/include/aare/Pedestal.hpp @@ -18,15 +18,15 @@ template class Pedestal { uint32_t m_samples; NDArray m_cur_samples; - - //TODO! in case of int needs to be changed to uint64_t + + // TODO! in case of int needs to be changed to uint64_t NDArray m_sum; NDArray m_sum2; - //Cache mean since it is used over and over in the ClusterFinder - //This optimization is related to the access pattern of the ClusterFinder - //Relies on having more reads than pushes to the pedestal - NDArray m_mean; + // Cache mean since it is used over and over in the ClusterFinder + // This optimization is related to the access pattern of the ClusterFinder + // Relies on having more reads than pushes to the pedestal + NDArray m_mean; public: Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000) @@ -42,9 +42,7 @@ template class Pedestal { } ~Pedestal() = default; - NDArray mean() { - return m_mean; - } + NDArray mean() { return m_mean; } SUM_TYPE mean(const uint32_t row, const uint32_t col) const { return m_mean(row, col); @@ -71,8 +69,6 @@ template class Pedestal { return variance_array; } - - NDArray std() { NDArray standard_deviation_array({m_rows, m_cols}); for (uint32_t i = 0; i < m_rows * m_cols; i++) { @@ -83,8 +79,6 @@ template class Pedestal { return standard_deviation_array; } - - void clear() { m_sum = 0; m_sum2 = 0; @@ -92,16 +86,12 @@ template class Pedestal { m_mean = 0; } - - void clear(const uint32_t row, const uint32_t col) { m_sum(row, col) = 0; m_sum2(row, col) = 0; m_cur_samples(row, col) = 0; m_mean(row, col) = 0; } - - template void push(NDView frame) { assert(frame.size() == m_rows * m_cols); @@ -122,7 +112,7 @@ template class Pedestal { /** * Push but don't update the cached mean. Speeds up the process * when initializing the pedestal. - * + * */ template void push_no_update(NDView frame) { assert(frame.size() == m_rows * m_cols); @@ -140,9 +130,6 @@ template class Pedestal { } } - - - template void push(Frame &frame) { assert(frame.rows() == static_cast(m_rows) && frame.cols() == static_cast(m_cols)); @@ -170,7 +157,8 @@ template class Pedestal { m_sum(row, col) += val - m_sum(row, col) / m_samples; m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples; } - //Since we just did a push we know that m_cur_samples(row, col) is at least 1 + // Since we just did a push we know that m_cur_samples(row, col) is at + // least 1 m_mean(row, col) = m_sum(row, col) / m_cur_samples(row, col); } @@ -183,7 +171,8 @@ template class Pedestal { m_cur_samples(row, col)++; } else { m_sum(row, col) += val - m_sum(row, col) / m_cur_samples(row, col); - m_sum2(row, col) += val * val - m_sum2(row, col) / m_cur_samples(row, col); + m_sum2(row, col) += + val * val - m_sum2(row, col) / m_cur_samples(row, col); } } @@ -191,19 +180,16 @@ template class Pedestal { * @brief Update the mean of the pedestal. This is used after having done * push_no_update. It is not necessary to call this function after push. */ - void update_mean(){ - m_mean = m_sum / m_cur_samples; - } + void update_mean() { m_mean = m_sum / m_cur_samples; } - template - void push_fast(const uint32_t row, const uint32_t col, const T val_){ - //Assume we reached the steady state where all pixels have - //m_samples samples + template + void push_fast(const uint32_t row, const uint32_t col, const T val_) { + // Assume we reached the steady state where all pixels have + // m_samples samples SUM_TYPE val = static_cast(val_); m_sum(row, col) += val - m_sum(row, col) / m_samples; m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples; m_mean(row, col) = m_sum(row, col) / m_samples; } - }; } // namespace aare \ No newline at end of file diff --git a/include/aare/PixelMap.hpp b/include/aare/PixelMap.hpp index 1b7a890..9c30680 100644 --- a/include/aare/PixelMap.hpp +++ b/include/aare/PixelMap.hpp @@ -1,7 +1,7 @@ #pragma once -#include "aare/defs.hpp" #include "aare/NDArray.hpp" +#include "aare/defs.hpp" namespace aare { @@ -10,11 +10,11 @@ NDArray GenerateMoench05PixelMap(); NDArray GenerateMoench05PixelMap1g(); NDArray GenerateMoench05PixelMapOld(); -//Matterhorn02 -NDArrayGenerateMH02SingleCounterPixelMap(); +// Matterhorn02 +NDArray GenerateMH02SingleCounterPixelMap(); NDArray GenerateMH02FourCounterPixelMap(); -//Eiger -NDArrayGenerateEigerFlipRowsPixelMap(); +// Eiger +NDArray GenerateEigerFlipRowsPixelMap(); } // namespace aare \ No newline at end of file diff --git a/include/aare/ProducerConsumerQueue.hpp b/include/aare/ProducerConsumerQueue.hpp index 426b9e2..f189cec 100644 --- a/include/aare/ProducerConsumerQueue.hpp +++ b/include/aare/ProducerConsumerQueue.hpp @@ -18,9 +18,9 @@ // @author Jordan DeLong (delong.j@fb.com) // Changes made by PSD Detector Group: -// Copied: Line 34 constexpr std::size_t hardware_destructive_interference_size = 128; from folly/lang/Align.h -// Changed extension to .hpp -// Changed namespace to aare +// Copied: Line 34 constexpr std::size_t hardware_destructive_interference_size +// = 128; from folly/lang/Align.h Changed extension to .hpp Changed namespace to +// aare #pragma once @@ -45,15 +45,14 @@ template struct ProducerConsumerQueue { ProducerConsumerQueue(const ProducerConsumerQueue &) = delete; ProducerConsumerQueue &operator=(const ProducerConsumerQueue &) = delete; - - ProducerConsumerQueue(ProducerConsumerQueue &&other){ + ProducerConsumerQueue(ProducerConsumerQueue &&other) { size_ = other.size_; records_ = other.records_; other.records_ = nullptr; readIndex_ = other.readIndex_.load(std::memory_order_acquire); writeIndex_ = other.writeIndex_.load(std::memory_order_acquire); } - ProducerConsumerQueue &operator=(ProducerConsumerQueue &&other){ + ProducerConsumerQueue &operator=(ProducerConsumerQueue &&other) { size_ = other.size_; records_ = other.records_; other.records_ = nullptr; @@ -61,16 +60,17 @@ template struct ProducerConsumerQueue { writeIndex_ = other.writeIndex_.load(std::memory_order_acquire); return *this; } - - - ProducerConsumerQueue():ProducerConsumerQueue(2){}; + + ProducerConsumerQueue() : ProducerConsumerQueue(2){}; // size must be >= 2. // // Also, note that the number of usable slots in the queue at any // given time is actually (size-1), so if you start with an empty queue, // isFull() will return true after size-1 insertions. explicit ProducerConsumerQueue(uint32_t size) - : size_(size), records_(static_cast(std::malloc(sizeof(T) * size))), readIndex_(0), writeIndex_(0) { + : size_(size), + records_(static_cast(std::malloc(sizeof(T) * size))), + readIndex_(0), writeIndex_(0) { assert(size >= 2); if (!records_) { throw std::bad_alloc(); @@ -154,7 +154,8 @@ template struct ProducerConsumerQueue { } bool isEmpty() const { - return readIndex_.load(std::memory_order_acquire) == writeIndex_.load(std::memory_order_acquire); + return readIndex_.load(std::memory_order_acquire) == + writeIndex_.load(std::memory_order_acquire); } bool isFull() const { @@ -175,7 +176,8 @@ template struct ProducerConsumerQueue { // be removing items concurrently). // * It is undefined to call this from any other thread. size_t sizeGuess() const { - int ret = writeIndex_.load(std::memory_order_acquire) - readIndex_.load(std::memory_order_acquire); + int ret = writeIndex_.load(std::memory_order_acquire) - + readIndex_.load(std::memory_order_acquire); if (ret < 0) { ret += size_; } @@ -192,7 +194,7 @@ template struct ProducerConsumerQueue { // const uint32_t size_; uint32_t size_; // T *const records_; - T* records_; + T *records_; alignas(hardware_destructive_interference_size) AtomicIndex readIndex_; alignas(hardware_destructive_interference_size) AtomicIndex writeIndex_; diff --git a/include/aare/RawFile.hpp b/include/aare/RawFile.hpp index 1cca1fd..9ffdb7c 100644 --- a/include/aare/RawFile.hpp +++ b/include/aare/RawFile.hpp @@ -1,11 +1,10 @@ #pragma once #include "aare/FileInterface.hpp" -#include "aare/RawMasterFile.hpp" #include "aare/Frame.hpp" #include "aare/NDArray.hpp" //for pixel map +#include "aare/RawMasterFile.hpp" #include "aare/RawSubFile.hpp" - #include namespace aare { @@ -53,10 +52,10 @@ class RawFile : public FileInterface { void read_into(std::byte *image_buf) override; void read_into(std::byte *image_buf, size_t n_frames) override; - //TODO! do we need to adapt the API? + // TODO! do we need to adapt the API? void read_into(std::byte *image_buf, DetectorHeader *header); - void read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header); - + void read_into(std::byte *image_buf, size_t n_frames, + DetectorHeader *header); size_t frame_number(size_t frame_index) override; size_t bytes_per_frame() override; @@ -70,23 +69,20 @@ class RawFile : public FileInterface { size_t bitdepth() const override; xy geometry(); size_t n_modules() const; - + RawMasterFile master() const; - - DetectorType detector_type() const override; private: - /** * @brief read the frame at the given frame index into the image buffer * @param frame_number frame number to read * @param image_buf buffer to store the frame */ - void get_frame_into(size_t frame_index, std::byte *frame_buffer, DetectorHeader *header = nullptr); - + void get_frame_into(size_t frame_index, std::byte *frame_buffer, + DetectorHeader *header = nullptr); /** * @brief get the frame at the given frame index @@ -95,8 +91,6 @@ class RawFile : public FileInterface { */ Frame get_frame(size_t frame_index); - - /** * @brief read the header of the file * @param fname path to the data subfile @@ -108,5 +102,4 @@ class RawFile : public FileInterface { void find_geometry(); }; - } // namespace aare \ No newline at end of file diff --git a/include/aare/RawMasterFile.hpp b/include/aare/RawMasterFile.hpp index 4d143a6..2c64a90 100644 --- a/include/aare/RawMasterFile.hpp +++ b/include/aare/RawMasterFile.hpp @@ -45,7 +45,7 @@ class ScanParameters { int m_start = 0; int m_stop = 0; int m_step = 0; - //TODO! add settleTime, requires string to time conversion + // TODO! add settleTime, requires string to time conversion public: ScanParameters(const std::string &par); @@ -61,7 +61,6 @@ class ScanParameters { void increment_stop(); }; - /** * @brief Class for parsing a master file either in our .json format or the old * .raw format @@ -101,7 +100,6 @@ class RawMasterFile { std::optional m_roi; - public: RawMasterFile(const std::filesystem::path &fpath); @@ -129,10 +127,8 @@ class RawMasterFile { std::optional number_of_rows() const; std::optional quad() const; - std::optional roi() const; - ScanParameters scan_parameters() const; private: diff --git a/include/aare/RawSubFile.hpp b/include/aare/RawSubFile.hpp index 1059843..c38f540 100644 --- a/include/aare/RawSubFile.hpp +++ b/include/aare/RawSubFile.hpp @@ -10,32 +10,34 @@ namespace aare { /** - * @brief Class to read a singe subfile written in .raw format. Used from RawFile to read - * the entire detector. Can be used directly to read part of the image. + * @brief Class to read a singe subfile written in .raw format. Used from + * RawFile to read the entire detector. Can be used directly to read part of the + * image. */ class RawSubFile { protected: std::ifstream m_file; DetectorType m_detector_type; size_t m_bitdepth; - std::filesystem::path m_path; //!< path to the subfile - std::string m_base_name; //!< base name used for formatting file names - size_t m_offset{}; //!< file index of the first file, allow starting at non zero file - size_t m_total_frames{}; //!< total number of frames in the series of files + std::filesystem::path m_path; //!< path to the subfile + std::string m_base_name; //!< base name used for formatting file names + size_t m_offset{}; //!< file index of the first file, allow starting at non + //!< zero file + size_t m_total_frames{}; //!< total number of frames in the series of files size_t m_rows{}; size_t m_cols{}; size_t m_bytes_per_frame{}; - int m_module_index{}; - 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_last_frame_in_file{}; //!< Used for seeking to the correct 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_last_frame_in_file{}; //!< Used for seeking to the correct file uint32_t m_pos_row{}; uint32_t m_pos_col{}; - - + std::optional> m_pixel_map; public: @@ -49,12 +51,14 @@ class RawSubFile { * @throws std::invalid_argument if the detector,type pair is not supported */ RawSubFile(const std::filesystem::path &fname, DetectorType detector, - size_t rows, size_t cols, size_t bitdepth, uint32_t pos_row = 0, uint32_t pos_col = 0); + size_t rows, size_t cols, size_t bitdepth, uint32_t pos_row = 0, + uint32_t pos_col = 0); ~RawSubFile() = default; /** * @brief Seek to the given frame number - * @note Puts the file pointer at the start of the header, not the start of the data + * @note Puts the file pointer at the start of the header, not the start of + * the data * @param frame_index frame position in file to seek to * @throws std::runtime_error if the frame number is out of range */ @@ -62,14 +66,15 @@ class RawSubFile { size_t tell(); void read_into(std::byte *image_buf, DetectorHeader *header = nullptr); - void read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header= nullptr); + void read_into(std::byte *image_buf, size_t n_frames, + DetectorHeader *header = nullptr); void get_part(std::byte *buffer, size_t frame_index); - + void read_header(DetectorHeader *header); - + size_t rows() const; size_t cols() const; - + size_t frame_number(size_t frame_index); size_t bytes_per_frame() const { return m_bytes_per_frame; } @@ -78,15 +83,13 @@ class RawSubFile { size_t frames_in_file() const { return m_total_frames; } -private: - template - void read_with_map(std::byte *image_buf); - - void parse_fname(const std::filesystem::path &fname); - void scan_files(); - void open_file(size_t file_index); - std::filesystem::path fpath(size_t file_index) const; + private: + template void read_with_map(std::byte *image_buf); + void parse_fname(const std::filesystem::path &fname); + void scan_files(); + void open_file(size_t file_index); + std::filesystem::path fpath(size_t file_index) const; }; } // namespace aare \ No newline at end of file diff --git a/include/aare/VarClusterFinder.hpp b/include/aare/VarClusterFinder.hpp index 596bf06..3679d39 100644 --- a/include/aare/VarClusterFinder.hpp +++ b/include/aare/VarClusterFinder.hpp @@ -38,11 +38,13 @@ template class VarClusterFinder { bool use_noise_map = false; int peripheralThresholdFactor_ = 5; int current_label; - const std::array di{{0, -1, -1, -1}}; // row ### 8-neighbour by scaning from left to right - const std::array dj{{-1, -1, 0, 1}}; // col ### 8-neighbour by scaning from top to bottom + const std::array di{ + {0, -1, -1, -1}}; // row ### 8-neighbour by scaning from left to right + const std::array dj{ + {-1, -1, 0, 1}}; // col ### 8-neighbour by scaning from top to bottom const std::array di_{{0, 0, -1, 1, -1, 1, -1, 1}}; // row const std::array dj_{{-1, 1, 0, 0, 1, -1, -1, 1}}; // col - std::map child; // heirachy: key: child; val: parent + std::map child; // heirachy: key: child; val: parent std::unordered_map h_size; std::vector hits; // std::vector> row @@ -50,7 +52,8 @@ template class VarClusterFinder { public: VarClusterFinder(Shape<2> shape, T threshold) - : shape_(shape), labeled_(shape, 0), peripheral_labeled_(shape, 0), binary_(shape), threshold_(threshold) { + : shape_(shape), labeled_(shape, 0), peripheral_labeled_(shape, 0), + binary_(shape), threshold_(threshold) { hits.reserve(2000); } @@ -60,7 +63,9 @@ template class VarClusterFinder { noiseMap = noise_map; use_noise_map = true; } - void set_peripheralThresholdFactor(int factor) { peripheralThresholdFactor_ = factor; } + void set_peripheralThresholdFactor(int factor) { + peripheralThresholdFactor_ = factor; + } void find_clusters(NDView img); void find_clusters_X(NDView img); void rec_FillHit(int clusterIndex, int i, int j); @@ -144,7 +149,8 @@ template int VarClusterFinder::check_neighbours(int i, int j) { } } -template void VarClusterFinder::find_clusters(NDView img) { +template +void VarClusterFinder::find_clusters(NDView img) { original_ = img; labeled_ = 0; peripheral_labeled_ = 0; @@ -156,7 +162,8 @@ template void VarClusterFinder::find_clusters(NDView img) store_clusters(); } -template void VarClusterFinder::find_clusters_X(NDView img) { +template +void VarClusterFinder::find_clusters_X(NDView img) { original_ = img; int clusterIndex = 0; for (int i = 0; i < shape_[0]; ++i) { @@ -175,7 +182,8 @@ template void VarClusterFinder::find_clusters_X(NDView img h_size.clear(); } -template void VarClusterFinder::rec_FillHit(int clusterIndex, int i, int j) { +template +void VarClusterFinder::rec_FillHit(int clusterIndex, int i, int j) { // printf("original_(%d, %d)=%f\n", i, j, original_(i,j)); // printf("h_size[%d].size=%d\n", clusterIndex, h_size[clusterIndex].size); if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE) { @@ -203,11 +211,15 @@ template void VarClusterFinder::rec_FillHit(int clusterIndex, in } else { // if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE){ // h_size[clusterIndex].size += 1; - // h_size[clusterIndex].rows[h_size[clusterIndex].size] = row; - // h_size[clusterIndex].cols[h_size[clusterIndex].size] = col; - // h_size[clusterIndex].enes[h_size[clusterIndex].size] = original_(row, col); + // h_size[clusterIndex].rows[h_size[clusterIndex].size] = + // row; h_size[clusterIndex].cols[h_size[clusterIndex].size] + // = col; + // h_size[clusterIndex].enes[h_size[clusterIndex].size] = + // original_(row, col); // }// ? weather to include peripheral pixels - original_(row, col) = 0; // remove peripheral pixels, to avoid potential influence for pedestal updating + original_(row, col) = + 0; // remove peripheral pixels, to avoid potential influence + // for pedestal updating } } } @@ -275,8 +287,8 @@ template void VarClusterFinder::store_clusters() { for (int i = 0; i < shape_[0]; ++i) { for (int j = 0; j < shape_[1]; ++j) { if (labeled_(i, j) != 0 || false - // (i-1 >= 0 and labeled_(i-1, j) != 0) or // another circle of peripheral pixels - // (j-1 >= 0 and labeled_(i, j-1) != 0) or + // (i-1 >= 0 and labeled_(i-1, j) != 0) or // another circle of + // peripheral pixels (j-1 >= 0 and labeled_(i, j-1) != 0) or // (i+1 < shape_[0] and labeled_(i+1, j) != 0) or // (j+1 < shape_[1] and labeled_(i, j+1) != 0) ) { diff --git a/include/aare/algorithm.hpp b/include/aare/algorithm.hpp index be2018f..c590e91 100644 --- a/include/aare/algorithm.hpp +++ b/include/aare/algorithm.hpp @@ -1,9 +1,9 @@ #pragma once +#include #include #include #include -#include namespace aare { /** @@ -15,26 +15,24 @@ namespace aare { * @param last iterator to the last element * @param val value to compare * @return index of the last element that is smaller than val - * + * */ template -size_t last_smaller(const T* first, const T* last, T val) { - for (auto iter = first+1; iter != last; ++iter) { +size_t last_smaller(const T *first, const T *last, T val) { + for (auto iter = first + 1; iter != last; ++iter) { if (*iter >= val) { - return std::distance(first, iter-1); + return std::distance(first, iter - 1); } } - return std::distance(first, last-1); + return std::distance(first, last - 1); } -template -size_t last_smaller(const NDArray& arr, T val) { +template size_t last_smaller(const NDArray &arr, T val) { return last_smaller(arr.begin(), arr.end(), val); } -template -size_t last_smaller(const std::vector& vec, T val) { - return last_smaller(vec.data(), vec.data()+vec.size(), val); +template size_t last_smaller(const std::vector &vec, T val) { + return last_smaller(vec.data(), vec.data() + vec.size(), val); } /** @@ -48,65 +46,59 @@ size_t last_smaller(const std::vector& vec, T val) { * @return index of the first element that is larger than val */ template -size_t first_larger(const T* first, const T* last, T val) { +size_t first_larger(const T *first, const T *last, T val) { for (auto iter = first; iter != last; ++iter) { if (*iter > val) { return std::distance(first, iter); } } - return std::distance(first, last-1); + return std::distance(first, last - 1); } -template -size_t first_larger(const NDArray& arr, T val) { +template size_t first_larger(const NDArray &arr, T val) { return first_larger(arr.begin(), arr.end(), val); } -template -size_t first_larger(const std::vector& vec, T val) { - return first_larger(vec.data(), vec.data()+vec.size(), val); +template size_t first_larger(const std::vector &vec, T val) { + return first_larger(vec.data(), vec.data() + vec.size(), val); } /** * @brief Index of the nearest element to val. - * Requires a sorted array. If there is no difference it takes the first element. + * Requires a sorted array. If there is no difference it takes the first + * element. * @param first iterator to the first element * @param last iterator to the last element * @param val value to compare * @return index of the nearest element */ template -size_t nearest_index(const T* first, const T* last, T val) { - auto iter = std::min_element(first, last, - [val](T a, T b) { +size_t nearest_index(const T *first, const T *last, T val) { + auto iter = std::min_element(first, last, [val](T a, T b) { return std::abs(a - val) < std::abs(b - val); }); return std::distance(first, iter); } -template -size_t nearest_index(const NDArray& arr, T val) { +template size_t nearest_index(const NDArray &arr, T val) { return nearest_index(arr.begin(), arr.end(), val); } -template -size_t nearest_index(const std::vector& vec, T val) { - return nearest_index(vec.data(), vec.data()+vec.size(), val); +template size_t nearest_index(const std::vector &vec, T val) { + return nearest_index(vec.data(), vec.data() + vec.size(), val); } template -size_t nearest_index(const std::array& arr, T val) { - return nearest_index(arr.data(), arr.data()+arr.size(), val); +size_t nearest_index(const std::array &arr, T val) { + return nearest_index(arr.data(), arr.data() + arr.size(), val); } -template -std::vector cumsum(const std::vector& vec) { +template std::vector cumsum(const std::vector &vec) { std::vector result(vec.size()); std::partial_sum(vec.begin(), vec.end(), result.begin()); return result; } - template bool all_equal(const Container &c) { if (!c.empty() && std::all_of(begin(c), end(c), @@ -117,6 +109,4 @@ template bool all_equal(const Container &c) { return false; } - - } // namespace aare \ No newline at end of file diff --git a/include/aare/decode.hpp b/include/aare/decode.hpp index e784c4a..ec24447 100644 --- a/include/aare/decode.hpp +++ b/include/aare/decode.hpp @@ -1,26 +1,27 @@ #pragma once +#include #include #include -#include namespace aare { - uint16_t adc_sar_05_decode64to16(uint64_t input); uint16_t adc_sar_04_decode64to16(uint64_t input); -void adc_sar_05_decode64to16(NDView input, NDView output); -void adc_sar_04_decode64to16(NDView input, NDView output); - +void adc_sar_05_decode64to16(NDView input, + NDView output); +void adc_sar_04_decode64to16(NDView input, + NDView output); /** - * @brief Apply custom weights to a 16-bit input value. Will sum up weights[i]**i - * for each bit i that is set in the input value. + * @brief Apply custom weights to a 16-bit input value. Will sum up + * weights[i]**i for each bit i that is set in the input value. * @throws std::out_of_range if weights.size() < 16 * @param input 16-bit input value * @param weights vector of weights, size must be less than or equal to 16 */ double apply_custom_weights(uint16_t input, const NDView weights); -void apply_custom_weights(NDView input, NDView output, const NDView weights); +void apply_custom_weights(NDView input, NDView output, + const NDView weights); } // namespace aare diff --git a/include/aare/defs.hpp b/include/aare/defs.hpp index ccf07a5..71d8c49 100644 --- a/include/aare/defs.hpp +++ b/include/aare/defs.hpp @@ -3,16 +3,15 @@ #include "aare/Dtype.hpp" #include -#include #include #include #include +#include #include #include #include #include - /** * @brief LOCATION macro to get the current location in the code */ @@ -20,28 +19,24 @@ std::string(__FILE__) + std::string(":") + std::to_string(__LINE__) + \ ":" + std::string(__func__) + ":" - - #ifdef AARE_CUSTOM_ASSERT -#define AARE_ASSERT(expr)\ - if (expr)\ - {}\ - else\ +#define AARE_ASSERT(expr) \ + if (expr) { \ + } else \ aare::assert_failed(LOCATION + " Assertion failed: " + #expr + "\n"); #else -#define AARE_ASSERT(cond)\ - do { (void)sizeof(cond); } while(0) +#define AARE_ASSERT(cond) \ + do { \ + (void)sizeof(cond); \ + } while (0) #endif - namespace aare { inline constexpr size_t bits_per_byte = 8; void assert_failed(const std::string &msg); - - class DynamicCluster { public: int cluster_sizeX; @@ -55,7 +50,7 @@ class DynamicCluster { public: DynamicCluster(int cluster_sizeX_, int cluster_sizeY_, - Dtype dt_ = Dtype(typeid(int32_t))) + Dtype dt_ = Dtype(typeid(int32_t))) : cluster_sizeX(cluster_sizeX_), cluster_sizeY(cluster_sizeY_), dt(dt_) { m_data = new std::byte[cluster_sizeX * cluster_sizeY * dt.bytes()]{}; @@ -179,24 +174,24 @@ template struct t_xy { }; using xy = t_xy; - /** - * @brief Class to hold the geometry of a module. Where pixel 0 is located and the size of the module + * @brief Class to hold the geometry of a module. Where pixel 0 is located and + * the size of the module */ -struct ModuleGeometry{ +struct ModuleGeometry { int origin_x{}; int origin_y{}; int height{}; int width{}; int row_index{}; - int col_index{}; + int col_index{}; }; /** - * @brief Class to hold the geometry of a detector. Number of modules, their size and where pixel 0 - * for each module is located + * @brief Class to hold the geometry of a detector. Number of modules, their + * size and where pixel 0 for each module is located */ -struct DetectorGeometry{ +struct DetectorGeometry { int modules_x{}; int modules_y{}; int pixels_x{}; @@ -204,35 +199,34 @@ struct DetectorGeometry{ int module_gap_row{}; int module_gap_col{}; std::vector module_pixel_0; - + auto size() const { return module_pixel_0.size(); } }; -struct ROI{ +struct ROI { ssize_t xmin{}; ssize_t xmax{}; ssize_t ymin{}; ssize_t ymax{}; - + ssize_t height() const { return ymax - ymin; } ssize_t width() const { return xmax - xmin; } bool contains(ssize_t x, ssize_t y) const { return x >= xmin && x < xmax && y >= ymin && y < ymax; } - }; - +}; using dynamic_shape = std::vector; -//TODO! Can we uniform enums between the libraries? +// TODO! Can we uniform enums between the libraries? /** - * @brief Enum class to identify different detectors. + * @brief Enum class to identify different detectors. * The values are the same as in slsDetectorPackage * Different spelling to avoid confusion with the slsDetectorPackage */ enum class DetectorType { - //Standard detectors match the enum values from slsDetectorPackage + // Standard detectors match the enum values from slsDetectorPackage Generic, Eiger, Gotthard, @@ -243,8 +237,9 @@ enum class DetectorType { Gotthard2, Xilinx_ChipTestBoard, - //Additional detectors used for defining processing. Variants of the standard ones. - Moench03=100, + // Additional detectors used for defining processing. Variants of the + // standard ones. + Moench03 = 100, Moench03_old, Unknown }; diff --git a/include/aare/geo_helpers.hpp b/include/aare/geo_helpers.hpp index d0d5d1a..c6454a5 100644 --- a/include/aare/geo_helpers.hpp +++ b/include/aare/geo_helpers.hpp @@ -1,16 +1,15 @@ #pragma once -#include "aare/defs.hpp" #include "aare/RawMasterFile.hpp" //ROI refactor away -namespace aare{ +#include "aare/defs.hpp" +namespace aare { /** * @brief Update the detector geometry given a region of interest - * - * @param geo - * @param roi - * @return DetectorGeometry + * + * @param geo + * @param roi + * @return DetectorGeometry */ DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, ROI roi); - } // namespace aare \ No newline at end of file diff --git a/include/aare/logger.hpp b/include/aare/logger.hpp index b93c091..0bedd7a 100644 --- a/include/aare/logger.hpp +++ b/include/aare/logger.hpp @@ -1,7 +1,6 @@ #pragma once /*Utility to log to console*/ - #include #include #include @@ -27,7 +26,6 @@ namespace aare { #define RESET "\x1b[0m" #define BOLD "\x1b[1m" - enum TLogLevel { logERROR, logWARNING, @@ -37,7 +35,8 @@ enum TLogLevel { logINFOCYAN, logINFOMAGENTA, logINFO, - logDEBUG, // constructors, destructors etc. should still give too much output + logDEBUG, // constructors, destructors etc. should still give too much + // output logDEBUG1, logDEBUG2, logDEBUG3, @@ -47,7 +46,9 @@ enum TLogLevel { // Compiler should optimize away anything below this value #ifndef AARE_LOG_LEVEL -#define AARE_LOG_LEVEL "LOG LEVEL NOT SET IN CMAKE" //This is configured in the main CMakeLists.txt +#define AARE_LOG_LEVEL \ + "LOG LEVEL NOT SET IN CMAKE" // This is configured in the main + // CMakeLists.txt #endif #define __AT__ \ @@ -72,7 +73,8 @@ class Logger { std::clog << os.str() << std::flush; // Single write } - static TLogLevel &ReportingLevel() { // singelton eeh TODO! Do we need a runtime option? + static TLogLevel & + ReportingLevel() { // singelton eeh TODO! Do we need a runtime option? static TLogLevel reportingLevel = logDEBUG5; return reportingLevel; } diff --git a/include/aare/utils/ifstream_helpers.hpp b/include/aare/utils/ifstream_helpers.hpp index 0a842ed..a8d0d21 100644 --- a/include/aare/utils/ifstream_helpers.hpp +++ b/include/aare/utils/ifstream_helpers.hpp @@ -4,9 +4,9 @@ #include namespace aare { -/** +/** * @brief Get the error message from an ifstream object -*/ + */ std::string ifstream_error_msg(std::ifstream &ifs); } // namespace aare \ No newline at end of file diff --git a/include/aare/utils/par.hpp b/include/aare/utils/par.hpp index efb1c77..e52c897 100644 --- a/include/aare/utils/par.hpp +++ b/include/aare/utils/par.hpp @@ -1,18 +1,18 @@ #include -#include #include +#include namespace aare { - template - void RunInParallel(F func, const std::vector>& tasks) { - // auto tasks = split_task(0, y.shape(0), n_threads); - std::vector threads; - for (auto &task : tasks) { - threads.push_back(std::thread(func, task.first, task.second)); - } - for (auto &thread : threads) { - thread.join(); - } +template +void RunInParallel(F func, const std::vector> &tasks) { + // auto tasks = split_task(0, y.shape(0), n_threads); + std::vector threads; + for (auto &task : tasks) { + threads.push_back(std::thread(func, task.first, task.second)); } + for (auto &thread : threads) { + thread.join(); + } +} } // namespace aare \ No newline at end of file diff --git a/python/src/bind_Cluster.hpp b/python/src/bind_Cluster.hpp index daf0946..690d0e8 100644 --- a/python/src/bind_Cluster.hpp +++ b/python/src/bind_Cluster.hpp @@ -1,10 +1,10 @@ #include "aare/Cluster.hpp" #include -#include #include -#include +#include #include +#include #include #include diff --git a/python/src/bind_ClusterCollector.hpp b/python/src/bind_ClusterCollector.hpp index 4836e6e..84172cb 100644 --- a/python/src/bind_ClusterCollector.hpp +++ b/python/src/bind_ClusterCollector.hpp @@ -21,11 +21,9 @@ using namespace aare; #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused-parameter" - template -void define_ClusterCollector(py::module &m, - const std::string &typestr) { +void define_ClusterCollector(py::module &m, const std::string &typestr) { auto class_name = fmt::format("ClusterCollector_{}", typestr); using ClusterType = Cluster; diff --git a/python/src/bind_ClusterFile.hpp b/python/src/bind_ClusterFile.hpp index 8ce5360..c2c801d 100644 --- a/python/src/bind_ClusterFile.hpp +++ b/python/src/bind_ClusterFile.hpp @@ -21,8 +21,7 @@ using namespace ::aare; template -void define_ClusterFile(py::module &m, - const std::string &typestr) { +void define_ClusterFile(py::module &m, const std::string &typestr) { using ClusterType = Cluster; diff --git a/python/src/bind_ClusterFileSink.hpp b/python/src/bind_ClusterFileSink.hpp index 9b3a74d..f717de6 100644 --- a/python/src/bind_ClusterFileSink.hpp +++ b/python/src/bind_ClusterFileSink.hpp @@ -21,15 +21,9 @@ using namespace aare; #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused-parameter" - - - - - template -void define_ClusterFileSink(py::module &m, - const std::string &typestr) { +void define_ClusterFileSink(py::module &m, const std::string &typestr) { auto class_name = fmt::format("ClusterFileSink_{}", typestr); using ClusterType = Cluster; @@ -40,5 +34,4 @@ void define_ClusterFileSink(py::module &m, .def("stop", &ClusterFileSink::stop); } - #pragma GCC diagnostic pop diff --git a/python/src/bind_ClusterFinderMT.hpp b/python/src/bind_ClusterFinderMT.hpp index d1769db..0ecbbd1 100644 --- a/python/src/bind_ClusterFinderMT.hpp +++ b/python/src/bind_ClusterFinderMT.hpp @@ -23,8 +23,7 @@ using namespace aare; template -void define_ClusterFinderMT(py::module &m, - const std::string &typestr) { +void define_ClusterFinderMT(py::module &m, const std::string &typestr) { auto class_name = fmt::format("ClusterFinderMT_{}", typestr); using ClusterType = Cluster; @@ -49,9 +48,11 @@ void define_ClusterFinderMT(py::module &m, return; }, py::arg(), py::arg("frame_number") = 0) - .def_property_readonly("cluster_size", [](ClusterFinderMT &self){ - return py::make_tuple(ClusterSizeX, ClusterSizeY); - }) + .def_property_readonly( + "cluster_size", + [](ClusterFinderMT &self) { + return py::make_tuple(ClusterSizeX, ClusterSizeY); + }) .def("clear_pedestal", &ClusterFinderMT::clear_pedestal) .def("sync", &ClusterFinderMT::sync) @@ -77,5 +78,4 @@ void define_ClusterFinderMT(py::module &m, py::arg("thread_index") = 0); } - #pragma GCC diagnostic pop \ No newline at end of file diff --git a/python/src/bind_ClusterVector.hpp b/python/src/bind_ClusterVector.hpp index 550db9a..9e9c4ab 100644 --- a/python/src/bind_ClusterVector.hpp +++ b/python/src/bind_ClusterVector.hpp @@ -44,10 +44,11 @@ void define_ClusterVector(py::module &m, const std::string &typestr) { auto *vec = new std::vector(self.sum()); return return_vector(vec); }) - .def("sum_2x2", [](ClusterVector &self){ - auto *vec = new std::vector(self.sum_2x2()); - return return_vector(vec); - }) + .def("sum_2x2", + [](ClusterVector &self) { + auto *vec = new std::vector(self.sum_2x2()); + return return_vector(vec); + }) .def_property_readonly("size", &ClusterVector::size) .def("item_size", &ClusterVector::item_size) .def_property_readonly("fmt", diff --git a/python/src/ctb_raw_file.hpp b/python/src/ctb_raw_file.hpp index c9b5310..5eb9652 100644 --- a/python/src/ctb_raw_file.hpp +++ b/python/src/ctb_raw_file.hpp @@ -6,8 +6,8 @@ #include "aare/RawMasterFile.hpp" #include "aare/RawSubFile.hpp" -#include "aare/defs.hpp" #include "aare/decode.hpp" +#include "aare/defs.hpp" // #include "aare/fClusterFileV2.hpp" #include "np_helper.hpp" @@ -26,95 +26,103 @@ using namespace ::aare; void define_ctb_raw_file_io_bindings(py::module &m) { -m.def("adc_sar_05_decode64to16", [](py::array_t input) { + m.def("adc_sar_05_decode64to16", [](py::array_t input) { + if (input.ndim() != 2) { + throw std::runtime_error( + "Only 2D arrays are supported at this moment"); + } - - if(input.ndim() != 2){ - throw std::runtime_error("Only 2D arrays are supported at this moment"); - } + // Create a 2D output array with the same shape as the input + std::vector shape{input.shape(0), + input.shape(1) / + static_cast(bits_per_byte)}; + py::array_t output(shape); - //Create a 2D output array with the same shape as the input - std::vector shape{input.shape(0), input.shape(1)/static_cast(bits_per_byte)}; - py::array_t output(shape); + // Create a view of the input and output arrays + NDView input_view( + reinterpret_cast(input.mutable_data()), + {output.shape(0), output.shape(1)}); + NDView output_view(output.mutable_data(), + {output.shape(0), output.shape(1)}); - //Create a view of the input and output arrays - NDView input_view(reinterpret_cast(input.mutable_data()), {output.shape(0), output.shape(1)}); - NDView output_view(output.mutable_data(), {output.shape(0), output.shape(1)}); + adc_sar_05_decode64to16(input_view, output_view); - adc_sar_05_decode64to16(input_view, output_view); - - return output; -}); - - -m.def("adc_sar_04_decode64to16", [](py::array_t input) { - - - if(input.ndim() != 2){ - throw std::runtime_error("Only 2D arrays are supported at this moment"); - } - - //Create a 2D output array with the same shape as the input - std::vector shape{input.shape(0), input.shape(1)/static_cast(bits_per_byte)}; - py::array_t output(shape); - - //Create a view of the input and output arrays - NDView input_view(reinterpret_cast(input.mutable_data()), {output.shape(0), output.shape(1)}); - NDView output_view(output.mutable_data(), {output.shape(0), output.shape(1)}); - - adc_sar_04_decode64to16(input_view, output_view); - - return output; -}); - -m.def( - "apply_custom_weights", - [](py::array_t &input, - py::array_t - &weights) { - - - // Create new array with same shape as the input array (uninitialized values) - py::buffer_info buf = input.request(); - py::array_t output(buf.shape); - - // Use NDViews to call into the C++ library - auto weights_view = make_view_1d(weights); - NDView input_view(input.mutable_data(), {input.size()}); - NDView output_view(output.mutable_data(), {output.size()}); - - apply_custom_weights(input_view, output_view, weights_view); return output; }); -py::class_(m, "CtbRawFile") - .def(py::init()) - .def("read_frame", - [](CtbRawFile &self) { - size_t image_size = self.image_size_in_bytes(); - py::array image; - std::vector shape; - shape.reserve(2); - shape.push_back(1); - shape.push_back(image_size); + m.def("adc_sar_04_decode64to16", [](py::array_t input) { + if (input.ndim() != 2) { + throw std::runtime_error( + "Only 2D arrays are supported at this moment"); + } - py::array_t header(1); + // Create a 2D output array with the same shape as the input + std::vector shape{input.shape(0), + input.shape(1) / + static_cast(bits_per_byte)}; + py::array_t output(shape); - // always read bytes - image = py::array_t(shape); + // Create a view of the input and output arrays + NDView input_view( + reinterpret_cast(input.mutable_data()), + {output.shape(0), output.shape(1)}); + NDView output_view(output.mutable_data(), + {output.shape(0), output.shape(1)}); - self.read_into(reinterpret_cast(image.mutable_data()), - header.mutable_data()); + adc_sar_04_decode64to16(input_view, output_view); - return py::make_tuple(header, image); - }) - .def("seek", &CtbRawFile::seek) - .def("tell", &CtbRawFile::tell) - .def("master", &CtbRawFile::master) + return output; + }); - .def_property_readonly("image_size_in_bytes", - &CtbRawFile::image_size_in_bytes) + m.def("apply_custom_weights", + [](py::array_t + &input, + py::array_t + &weights) { + // Create new array with same shape as the input array + // (uninitialized values) + py::buffer_info buf = input.request(); + py::array_t output(buf.shape); - .def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file); + // Use NDViews to call into the C++ library + auto weights_view = make_view_1d(weights); + NDView input_view(input.mutable_data(), + {input.size()}); + NDView output_view(output.mutable_data(), + {output.size()}); + apply_custom_weights(input_view, output_view, weights_view); + return output; + }); + + py::class_(m, "CtbRawFile") + .def(py::init()) + .def("read_frame", + [](CtbRawFile &self) { + size_t image_size = self.image_size_in_bytes(); + py::array image; + std::vector shape; + shape.reserve(2); + shape.push_back(1); + shape.push_back(image_size); + + py::array_t header(1); + + // always read bytes + image = py::array_t(shape); + + self.read_into( + reinterpret_cast(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); } diff --git a/python/src/file.hpp b/python/src/file.hpp index f97db96..262b4f8 100644 --- a/python/src/file.hpp +++ b/python/src/file.hpp @@ -20,17 +20,13 @@ namespace py = pybind11; using namespace ::aare; - - - -//Disable warnings for unused parameters, as we ignore some -//in the __exit__ method +// Disable warnings for unused parameters, as we ignore some +// in the __exit__ method #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused-parameter" void define_file_io_bindings(py::module &m) { - py::enum_(m, "DetectorType") .value("Jungfrau", DetectorType::Jungfrau) .value("Eiger", DetectorType::Eiger) @@ -41,13 +37,10 @@ void define_file_io_bindings(py::module &m) { .value("ChipTestBoard", DetectorType::ChipTestBoard) .value("Unknown", DetectorType::Unknown); - PYBIND11_NUMPY_DTYPE(DetectorHeader, frameNumber, expLength, packetNumber, bunchId, timestamp, modId, row, column, reserved, debug, roundRNumber, detType, version, packetMask); - - py::class_(m, "File") .def(py::init([](const std::filesystem::path &fname) { return File(fname, "r", {}); @@ -112,45 +105,18 @@ void define_file_io_bindings(py::module &m) { reinterpret_cast(image.mutable_data())); return image; }) - .def("read_n", [](File &self, size_t n_frames) { - //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 shape{n_frames, self.rows(), self.cols()}; - - py::array image; - const uint8_t item_size = self.bytes_per_pixel(); - if (item_size == 1) { - image = py::array_t(shape); - } else if (item_size == 2) { - image = py::array_t(shape); - } else if (item_size == 4) { - image = py::array_t(shape); - } - self.read_into(reinterpret_cast(image.mutable_data()), - n_frames); - return image; - }) - .def("__enter__", [](File &self) { return &self; }) - .def("__exit__", - [](File &self, - const std::optional &exc_type, - const std::optional &exc_value, - const std::optional &traceback) { - // self.close(); - }) - .def("__iter__", [](File &self) { return &self; }) - .def("__next__", [](File &self) { + .def("read_n", + [](File &self, size_t n_frames) { + // 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 shape{n_frames, self.rows(), self.cols()}; - try{ - const uint8_t item_size = self.bytes_per_pixel(); py::array image; - std::vector shape; - shape.reserve(2); - shape.push_back(self.rows()); - shape.push_back(self.cols()); + const uint8_t item_size = self.bytes_per_pixel(); if (item_size == 1) { image = py::array_t(shape); } else if (item_size == 2) { @@ -159,14 +125,41 @@ void define_file_io_bindings(py::module &m) { image = py::array_t(shape); } self.read_into( - reinterpret_cast(image.mutable_data())); + reinterpret_cast(image.mutable_data()), + n_frames); return image; - }catch(std::runtime_error &e){ + }) + .def("__enter__", [](File &self) { return &self; }) + .def("__exit__", + [](File &self, const std::optional &exc_type, + const std::optional &exc_value, + const std::optional &traceback) { + // self.close(); + }) + .def("__iter__", [](File &self) { return &self; }) + .def("__next__", [](File &self) { + try { + const uint8_t item_size = self.bytes_per_pixel(); + py::array image; + std::vector shape; + shape.reserve(2); + shape.push_back(self.rows()); + shape.push_back(self.cols()); + if (item_size == 1) { + image = py::array_t(shape); + } else if (item_size == 2) { + image = py::array_t(shape); + } else if (item_size == 4) { + image = py::array_t(shape); + } + self.read_into( + reinterpret_cast(image.mutable_data())); + return image; + } catch (std::runtime_error &e) { throw py::stop_iteration(); } }); - py::class_(m, "FileConfig") .def(py::init<>()) .def_readwrite("rows", &FileConfig::rows) @@ -183,8 +176,6 @@ void define_file_io_bindings(py::module &m) { return ""; }); - - py::class_(m, "ScanParameters") .def(py::init()) .def(py::init()) @@ -195,7 +186,6 @@ void define_file_io_bindings(py::module &m) { .def_property_readonly("stop", &ScanParameters::stop) .def_property_readonly("step", &ScanParameters::step); - py::class_(m, "ROI") .def(py::init<>()) .def(py::init(), py::arg("xmin"), @@ -204,23 +194,21 @@ void define_file_io_bindings(py::module &m) { .def_readwrite("xmax", &ROI::xmax) .def_readwrite("ymin", &ROI::ymin) .def_readwrite("ymax", &ROI::ymax) - .def("__str__", [](const ROI& self){ - return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}", self.xmin, self.xmax, self.ymin, self.ymax); - }) - .def("__repr__", [](const ROI& self){ - return fmt::format("", self.xmin, self.xmax, self.ymin, self.ymax); - }) + .def("__str__", + [](const ROI &self) { + return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}", + self.xmin, self.xmax, self.ymin, self.ymax); + }) + .def("__repr__", + [](const ROI &self) { + return fmt::format( + "", self.xmin, + self.xmax, self.ymin, self.ymax); + }) .def("__iter__", [](const ROI &self) { - return py::make_iterator(&self.xmin, &self.ymax+1); //NOLINT + return py::make_iterator(&self.xmin, &self.ymax + 1); // NOLINT }); - - - - - - - #pragma GCC diagnostic pop // py::class_(m, "ClusterHeader") // .def(py::init<>()) diff --git a/python/src/fit.hpp b/python/src/fit.hpp index 97dafb5..47568d6 100644 --- a/python/src/fit.hpp +++ b/python/src/fit.hpp @@ -9,7 +9,6 @@ namespace py = pybind11; using namespace pybind11::literals; - void define_fit_bindings(py::module &m) { // TODO! Evaluate without converting to double @@ -61,7 +60,8 @@ void define_fit_bindings(py::module &m) { py::array_t par) { auto x_view = make_view_1d(x); auto par_view = make_view_1d(par); - auto y = new NDArray{aare::func::scurve(x_view, par_view)}; + auto y = + new NDArray{aare::func::scurve(x_view, par_view)}; return return_image_data(y); }, R"( @@ -82,7 +82,8 @@ void define_fit_bindings(py::module &m) { py::array_t par) { auto x_view = make_view_1d(x); auto par_view = make_view_1d(par); - auto y = new NDArray{aare::func::scurve2(x_view, par_view)}; + auto y = + new NDArray{aare::func::scurve2(x_view, par_view)}; return return_image_data(y); }, R"( @@ -139,7 +140,6 @@ n_threads : int, optional py::array_t y, py::array_t y_err, int n_threads) { - if (y.ndim() == 3) { // Allocate memory for the output // Need to have pointers to allow python to manage @@ -173,7 +173,6 @@ n_threads : int, optional auto y_view_err = make_view_1d(y_err); auto x_view = make_view_1d(x); - double chi2 = 0; aare::fit_gaus(x_view, y_view, y_view_err, par->view(), par_err->view(), chi2); @@ -248,11 +247,10 @@ n_threads : int, optional aare::fit_pol1(x_view, y_view, y_view_err, par->view(), par_err->view(), chi2->view(), n_threads); return py::dict("par"_a = return_image_data(par), - "par_err"_a = return_image_data(par_err), + "par_err"_a = return_image_data(par_err), "chi2"_a = return_image_data(chi2), "Ndf"_a = y.shape(2) - 2); - } else if (y.ndim() == 1) { auto par = new NDArray({2}); auto par_err = new NDArray({2}); @@ -289,7 +287,7 @@ n_threads : int, optional )", py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4); -//========= + //========= m.def( "fit_scurve", [](py::array_t x, @@ -333,13 +331,12 @@ n_threads : int, optional auto chi2 = new NDArray({y.shape(0), y.shape(1)}); aare::fit_scurve(x_view, y_view, y_view_err, par->view(), - par_err->view(), chi2->view(), n_threads); + par_err->view(), chi2->view(), n_threads); return py::dict("par"_a = return_image_data(par), - "par_err"_a = return_image_data(par_err), + "par_err"_a = return_image_data(par_err), "chi2"_a = return_image_data(chi2), "Ndf"_a = y.shape(2) - 2); - } else if (y.ndim() == 1) { auto par = new NDArray({2}); auto par_err = new NDArray({2}); @@ -351,7 +348,7 @@ n_threads : int, optional double chi2 = 0; aare::fit_scurve(x_view, y_view, y_view_err, par->view(), - par_err->view(), chi2); + par_err->view(), chi2); return py::dict("par"_a = return_image_data(par), "par_err"_a = return_image_data(par_err), "chi2"_a = chi2, "Ndf"_a = y.size() - 2); @@ -375,7 +372,6 @@ n_threads : int, optional The number of threads to use. Default is 4. )", py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4); - m.def( "fit_scurve2", @@ -420,13 +416,12 @@ n_threads : int, optional auto chi2 = new NDArray({y.shape(0), y.shape(1)}); aare::fit_scurve2(x_view, y_view, y_view_err, par->view(), - par_err->view(), chi2->view(), n_threads); + par_err->view(), chi2->view(), n_threads); return py::dict("par"_a = return_image_data(par), - "par_err"_a = return_image_data(par_err), + "par_err"_a = return_image_data(par_err), "chi2"_a = return_image_data(chi2), "Ndf"_a = y.shape(2) - 2); - } else if (y.ndim() == 1) { auto par = new NDArray({6}); auto par_err = new NDArray({6}); @@ -438,7 +433,7 @@ n_threads : int, optional double chi2 = 0; aare::fit_scurve2(x_view, y_view, y_view_err, par->view(), - par_err->view(), chi2); + par_err->view(), chi2); return py::dict("par"_a = return_image_data(par), "par_err"_a = return_image_data(par_err), "chi2"_a = chi2, "Ndf"_a = y.size() - 2); diff --git a/python/src/jungfrau_data_file.hpp b/python/src/jungfrau_data_file.hpp index 942f6a6..62a95c9 100644 --- a/python/src/jungfrau_data_file.hpp +++ b/python/src/jungfrau_data_file.hpp @@ -21,10 +21,7 @@ using namespace ::aare; auto read_dat_frame(JungfrauDataFile &self) { py::array_t header(1); - py::array_t image({ - self.rows(), - self.cols() - }); + py::array_t image({self.rows(), self.cols()}); self.read_into(reinterpret_cast(image.mutable_data()), header.mutable_data()); @@ -40,9 +37,7 @@ auto read_n_dat_frames(JungfrauDataFile &self, size_t n_frames) { } py::array_t header(n_frames); - py::array_t image({ - n_frames, self.rows(), - self.cols()}); + py::array_t image({n_frames, self.rows(), self.cols()}); self.read_into(reinterpret_cast(image.mutable_data()), n_frames, header.mutable_data()); diff --git a/python/src/module.cpp b/python/src/module.cpp index 681dd4b..fc04a9f 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -1,26 +1,26 @@ // Files with bindings to the different classes -//New style file naming +// New style file naming #include "bind_Cluster.hpp" #include "bind_ClusterCollector.hpp" -#include "bind_ClusterFinder.hpp" -#include "bind_ClusterFinderMT.hpp" #include "bind_ClusterFile.hpp" #include "bind_ClusterFileSink.hpp" +#include "bind_ClusterFinder.hpp" +#include "bind_ClusterFinderMT.hpp" #include "bind_ClusterVector.hpp" -//TODO! migrate the other names +// TODO! migrate the other names #include "ctb_raw_file.hpp" #include "file.hpp" #include "fit.hpp" #include "interpolation.hpp" -#include "raw_sub_file.hpp" -#include "raw_master_file.hpp" -#include "raw_file.hpp" -#include "pixel_map.hpp" -#include "var_cluster.hpp" -#include "pedestal.hpp" #include "jungfrau_data_file.hpp" +#include "pedestal.hpp" +#include "pixel_map.hpp" +#include "raw_file.hpp" +#include "raw_master_file.hpp" +#include "raw_sub_file.hpp" +#include "var_cluster.hpp" // Pybind stuff #include @@ -34,17 +34,18 @@ T - Storage type of the cluster data (int, float, double) N - Number of rows in the cluster M - Number of columns in the cluster U - Type of the pixel data (e.g., uint16_t) -TYPE_CODE - A character representing the type code (e.g., 'i' for int, 'd' for double, 'f' for float) +TYPE_CODE - A character representing the type code (e.g., 'i' for int, 'd' for +double, 'f' for float) */ -#define DEFINE_CLUSTER_BINDINGS(T, N, M, U, TYPE_CODE) \ - define_ClusterFile(m, "Cluster" #N "x" #M #TYPE_CODE); \ - define_ClusterVector(m, "Cluster" #N "x" #M #TYPE_CODE); \ - define_ClusterFinder(m, "Cluster" #N "x" #M #TYPE_CODE); \ - define_ClusterFinderMT(m, "Cluster" #N "x" #M #TYPE_CODE); \ - define_ClusterFileSink(m, "Cluster" #N "x" #M #TYPE_CODE); \ - define_ClusterCollector(m, "Cluster" #N "x" #M #TYPE_CODE); \ - define_Cluster(m, #N "x" #M #TYPE_CODE); \ +#define DEFINE_CLUSTER_BINDINGS(T, N, M, U, TYPE_CODE) \ + define_ClusterFile(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_ClusterVector(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_ClusterFinder(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_ClusterFinderMT(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_ClusterFileSink(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_ClusterCollector(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_Cluster(m, #N "x" #M #TYPE_CODE); \ register_calculate_eta(m); PYBIND11_MODULE(_aare, m) { diff --git a/python/src/pedestal.hpp b/python/src/pedestal.hpp index 77148dc..23d8247 100644 --- a/python/src/pedestal.hpp +++ b/python/src/pedestal.hpp @@ -9,7 +9,8 @@ namespace py = pybind11; -template void define_pedestal_bindings(py::module &m, const std::string &name) { +template +void define_pedestal_bindings(py::module &m, const std::string &name) { py::class_>(m, name.c_str()) .def(py::init()) .def(py::init()) @@ -19,16 +20,18 @@ template void define_pedestal_bindings(py::module &m, const *mea = self.mean(); return return_image_data(mea); }) - .def("variance", [](Pedestal &self) { - auto var = new NDArray{}; - *var = self.variance(); - return return_image_data(var); - }) - .def("std", [](Pedestal &self) { - auto std = new NDArray{}; - *std = self.std(); - return return_image_data(std); - }) + .def("variance", + [](Pedestal &self) { + auto var = new NDArray{}; + *var = self.variance(); + return return_image_data(var); + }) + .def("std", + [](Pedestal &self) { + auto std = new NDArray{}; + *std = self.std(); + return return_image_data(std); + }) .def("clear", py::overload_cast<>(&Pedestal::clear)) .def_property_readonly("rows", &Pedestal::rows) .def_property_readonly("cols", &Pedestal::cols) @@ -39,14 +42,19 @@ template void define_pedestal_bindings(py::module &m, const [&](Pedestal &pedestal) { return Pedestal(pedestal); }) - //TODO! add push for other data types - .def("push", [](Pedestal &pedestal, py::array_t &f) { - auto v = make_view_2d(f); - pedestal.push(v); - }) - .def("push_no_update", [](Pedestal &pedestal, py::array_t &f) { - auto v = make_view_2d(f); - pedestal.push_no_update(v); - }, py::arg().noconvert()) + // TODO! add push for other data types + .def("push", + [](Pedestal &pedestal, py::array_t &f) { + auto v = make_view_2d(f); + pedestal.push(v); + }) + .def( + "push_no_update", + [](Pedestal &pedestal, + py::array_t &f) { + auto v = make_view_2d(f); + pedestal.push_no_update(v); + }, + py::arg().noconvert()) .def("update_mean", &Pedestal::update_mean); } \ No newline at end of file diff --git a/python/src/pixel_map.hpp b/python/src/pixel_map.hpp index 46b1bc4..986728b 100644 --- a/python/src/pixel_map.hpp +++ b/python/src/pixel_map.hpp @@ -1,41 +1,46 @@ #include "aare/PixelMap.hpp" #include "np_helper.hpp" - #include #include #include #include - namespace py = pybind11; -using namespace::aare; - +using namespace ::aare; void define_pixel_map_bindings(py::module &m) { - m.def("GenerateMoench03PixelMap", []() { - auto ptr = new NDArray(GenerateMoench03PixelMap()); - return return_image_data(ptr); - }) - .def("GenerateMoench05PixelMap", []() { - auto ptr = new NDArray(GenerateMoench05PixelMap()); - return return_image_data(ptr); - }) - .def("GenerateMoench05PixelMap1g", []() { - auto ptr = new NDArray(GenerateMoench05PixelMap1g()); - return return_image_data(ptr); - }) - .def("GenerateMoench05PixelMapOld", []() { - auto ptr = new NDArray(GenerateMoench05PixelMapOld()); - return return_image_data(ptr); - }) - .def("GenerateMH02SingleCounterPixelMap", []() { - auto ptr = new NDArray(GenerateMH02SingleCounterPixelMap()); - return return_image_data(ptr); - }) - .def("GenerateMH02FourCounterPixelMap", []() { - auto ptr = new NDArray(GenerateMH02FourCounterPixelMap()); - return return_image_data(ptr); - }); - + m.def("GenerateMoench03PixelMap", + []() { + auto ptr = new NDArray(GenerateMoench03PixelMap()); + return return_image_data(ptr); + }) + .def("GenerateMoench05PixelMap", + []() { + auto ptr = new NDArray(GenerateMoench05PixelMap()); + return return_image_data(ptr); + }) + .def("GenerateMoench05PixelMap1g", + []() { + auto ptr = + new NDArray(GenerateMoench05PixelMap1g()); + return return_image_data(ptr); + }) + .def("GenerateMoench05PixelMapOld", + []() { + auto ptr = + new NDArray(GenerateMoench05PixelMapOld()); + return return_image_data(ptr); + }) + .def("GenerateMH02SingleCounterPixelMap", + []() { + auto ptr = new NDArray( + GenerateMH02SingleCounterPixelMap()); + return return_image_data(ptr); + }) + .def("GenerateMH02FourCounterPixelMap", []() { + auto ptr = + new NDArray(GenerateMH02FourCounterPixelMap()); + return return_image_data(ptr); + }); } \ No newline at end of file diff --git a/python/src/raw_file.hpp b/python/src/raw_file.hpp index 8d72220..689b84e 100644 --- a/python/src/raw_file.hpp +++ b/python/src/raw_file.hpp @@ -58,13 +58,14 @@ void define_raw_file_io_bindings(py::module &m) { throw std::runtime_error("No frames left in file"); } std::vector shape{n_frames, self.rows(), self.cols()}; - + // return headers from all subfiles py::array_t header; if (self.n_modules() == 1) { header = py::array_t(n_frames); } else { - header = py::array_t({self.n_modules(), n_frames}); + header = py::array_t( + {self.n_modules(), n_frames}); } // py::array_t header({self.n_mod(), n_frames}); diff --git a/python/src/raw_master_file.hpp b/python/src/raw_master_file.hpp index 943437f..9c2bd17 100644 --- a/python/src/raw_master_file.hpp +++ b/python/src/raw_master_file.hpp @@ -57,7 +57,8 @@ void define_raw_master_file_bindings(py::module &m) { .def_property_readonly("total_frames_expected", &RawMasterFile::total_frames_expected) .def_property_readonly("geometry", &RawMasterFile::geometry) - .def_property_readonly("analog_samples", &RawMasterFile::analog_samples, R"( + .def_property_readonly("analog_samples", &RawMasterFile::analog_samples, + R"( Number of analog samples Returns @@ -66,7 +67,7 @@ void define_raw_master_file_bindings(py::module &m) { The number of analog samples in the file (or None if not enabled) )") .def_property_readonly("digital_samples", - &RawMasterFile::digital_samples, R"( + &RawMasterFile::digital_samples, R"( Number of digital samples Returns diff --git a/python/src/raw_sub_file.hpp b/python/src/raw_sub_file.hpp index 2cb83fc..cff511b 100644 --- a/python/src/raw_sub_file.hpp +++ b/python/src/raw_sub_file.hpp @@ -24,8 +24,8 @@ auto read_frame_from_RawSubFile(RawSubFile &self) { py::array_t header(1); const uint8_t item_size = self.bytes_per_pixel(); std::vector shape{static_cast(self.rows()), - static_cast(self.cols())}; - + static_cast(self.cols())}; + py::array image; if (item_size == 1) { image = py::array_t(shape); @@ -43,12 +43,10 @@ auto read_frame_from_RawSubFile(RawSubFile &self) { auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) { py::array_t header(n_frames); const uint8_t item_size = self.bytes_per_pixel(); - std::vector shape{ - static_cast(n_frames), - static_cast(self.rows()), - static_cast(self.cols()) - }; - + std::vector shape{static_cast(n_frames), + static_cast(self.rows()), + static_cast(self.cols())}; + py::array image; if (item_size == 1) { image = py::array_t(shape); @@ -57,15 +55,14 @@ auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) { } else if (item_size == 4) { image = py::array_t(shape); } - self.read_into(reinterpret_cast(image.mutable_data()), n_frames, - header.mutable_data()); + self.read_into(reinterpret_cast(image.mutable_data()), + n_frames, header.mutable_data()); return py::make_tuple(header, image); } - -//Disable warnings for unused parameters, as we ignore some -//in the __exit__ method +// Disable warnings for unused parameters, as we ignore some +// in the __exit__ method #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused-parameter" @@ -76,7 +73,7 @@ void define_raw_sub_file_io_bindings(py::module &m) { .def_property_readonly("bytes_per_frame", &RawSubFile::bytes_per_frame) .def_property_readonly("pixels_per_frame", &RawSubFile::pixels_per_frame) - .def_property_readonly("bytes_per_pixel", &RawSubFile::bytes_per_pixel) + .def_property_readonly("bytes_per_pixel", &RawSubFile::bytes_per_pixel) .def("seek", &RawSubFile::seek) .def("tell", &RawSubFile::tell) .def_property_readonly("rows", &RawSubFile::rows) @@ -84,18 +81,17 @@ void define_raw_sub_file_io_bindings(py::module &m) { .def_property_readonly("frames_in_file", &RawSubFile::frames_in_file) .def("read_frame", &read_frame_from_RawSubFile) .def("read_n", &read_n_frames_from_RawSubFile) - .def("read", [](RawSubFile &self){ - self.seek(0); - auto n_frames = self.frames_in_file(); - return read_n_frames_from_RawSubFile(self, n_frames); - }) + .def("read", + [](RawSubFile &self) { + self.seek(0); + auto n_frames = self.frames_in_file(); + return read_n_frames_from_RawSubFile(self, n_frames); + }) .def("__enter__", [](RawSubFile &self) { return &self; }) .def("__exit__", - [](RawSubFile &self, - const std::optional &exc_type, + [](RawSubFile &self, const std::optional &exc_type, const std::optional &exc_value, - const std::optional &traceback) { - }) + const std::optional &traceback) {}) .def("__iter__", [](RawSubFile &self) { return &self; }) .def("__next__", [](RawSubFile &self) { try { @@ -104,7 +100,6 @@ void define_raw_sub_file_io_bindings(py::module &m) { throw py::stop_iteration(); } }); - } #pragma GCC diagnostic pop \ No newline at end of file diff --git a/python/src/var_cluster.hpp b/python/src/var_cluster.hpp index f7b373f..4e7302d 100644 --- a/python/src/var_cluster.hpp +++ b/python/src/var_cluster.hpp @@ -12,10 +12,8 @@ // #include // #include - namespace py = pybind11; -using namespace::aare; - +using namespace ::aare; void define_var_cluster_finder_bindings(py::module &m) { PYBIND11_NUMPY_DTYPE(VarClusterFinder::Hit, size, row, col, @@ -29,12 +27,12 @@ void define_var_cluster_finder_bindings(py::module &m) { return return_image_data(ptr); }) .def("set_noiseMap", - [](VarClusterFinder &self, + [](VarClusterFinder &self, py::array_t noise_map) { - auto noise_map_span = make_view_2d(noise_map); - self.set_noiseMap(noise_map_span); - }) + auto noise_map_span = make_view_2d(noise_map); + self.set_noiseMap(noise_map_span); + }) .def("set_peripheralThresholdFactor", &VarClusterFinder::set_peripheralThresholdFactor) .def("find_clusters", @@ -65,9 +63,7 @@ void define_var_cluster_finder_bindings(py::module &m) { return return_vector(ptr); }) .def("clear_hits", - [](VarClusterFinder &self) { - self.clear_hits(); - }) + [](VarClusterFinder &self) { self.clear_hits(); }) .def("steal_hits", [](VarClusterFinder &self) { auto ptr = new std::vector::Hit>( @@ -75,5 +71,4 @@ void define_var_cluster_finder_bindings(py::module &m) { return return_vector(ptr); }) .def("total_clusters", &VarClusterFinder::total_clusters); - } \ No newline at end of file diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp index d24e803..13f7364 100644 --- a/src/ClusterFile.cpp +++ b/src/ClusterFile.cpp @@ -31,17 +31,15 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size, } } -void ClusterFile::set_roi(ROI roi){ - m_roi = roi; -} +void ClusterFile::set_roi(ROI roi) { m_roi = roi; } -void ClusterFile::set_noise_map(const NDView noise_map){ +void ClusterFile::set_noise_map(const NDView noise_map) { m_noise_map = NDArray(noise_map); } -void ClusterFile::set_gain_map(const NDView gain_map){ +void ClusterFile::set_gain_map(const NDView gain_map) { m_gain_map = NDArray(gain_map); - + // Gain map is passed as ADU/keV to avoid dividing in when applying the gain // map we invert it here for (auto &item : m_gain_map->view()) { @@ -66,42 +64,44 @@ void ClusterFile::write_frame(const ClusterVector &clusters) { !(clusters.cluster_size_y() == 3)) { throw std::runtime_error("Only 3x3 clusters are supported"); } - //First write the frame number - 4 bytes + // First write the frame number - 4 bytes int32_t frame_number = clusters.frame_number(); - if(fwrite(&frame_number, sizeof(frame_number), 1, fp)!=1){ + if (fwrite(&frame_number, sizeof(frame_number), 1, fp) != 1) { throw std::runtime_error(LOCATION + "Could not write frame number"); } - //Then write the number of clusters - 4 bytes + // Then write the number of clusters - 4 bytes uint32_t n_clusters = clusters.size(); - if(fwrite(&n_clusters, sizeof(n_clusters), 1, fp)!=1){ - throw std::runtime_error(LOCATION + "Could not write number of clusters"); + if (fwrite(&n_clusters, sizeof(n_clusters), 1, fp) != 1) { + throw std::runtime_error(LOCATION + + "Could not write number of clusters"); } - //Now write the clusters in the frame - if(fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp)!=clusters.size()){ + // Now write the clusters in the frame + if (fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp) != + clusters.size()) { throw std::runtime_error(LOCATION + "Could not write clusters"); } } - -ClusterVector ClusterFile::read_clusters(size_t n_clusters){ +ClusterVector ClusterFile::read_clusters(size_t n_clusters) { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); } - if (m_noise_map || m_roi){ + if (m_noise_map || m_roi) { return read_clusters_with_cut(n_clusters); - }else{ + } else { return read_clusters_without_cut(n_clusters); } } -ClusterVector ClusterFile::read_clusters_without_cut(size_t n_clusters) { +ClusterVector +ClusterFile::read_clusters_without_cut(size_t n_clusters) { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); } - - ClusterVector clusters(3,3, n_clusters); + + ClusterVector clusters(3, 3, n_clusters); int32_t iframe = 0; // frame number needs to be 4 bytes! size_t nph_read = 0; @@ -119,7 +119,7 @@ ClusterVector ClusterFile::read_clusters_without_cut(size_t n_clusters) } else { nn = nph; } - nph_read += fread((buf + nph_read*clusters.item_size()), + nph_read += fread((buf + nph_read * clusters.item_size()), clusters.item_size(), nn, fp); m_num_left = nph - nn; // write back the number of photons left } @@ -135,7 +135,7 @@ ClusterVector ClusterFile::read_clusters_without_cut(size_t n_clusters) else nn = nph; - nph_read += fread((buf + nph_read*clusters.item_size()), + nph_read += fread((buf + nph_read * clusters.item_size()), clusters.item_size(), nn, fp); m_num_left = nph - nn; } @@ -147,22 +147,22 @@ ClusterVector ClusterFile::read_clusters_without_cut(size_t n_clusters) // Resize the vector to the number of clusters. // No new allocation, only change bounds. clusters.resize(nph_read); - if(m_gain_map) + if (m_gain_map) clusters.apply_gain_map(m_gain_map->view()); return clusters; } - ClusterVector ClusterFile::read_clusters_with_cut(size_t n_clusters) { - ClusterVector clusters(3,3); + ClusterVector clusters(3, 3); clusters.reserve(n_clusters); // if there are photons left from previous frame read them first if (m_num_left) { - while(m_num_left && clusters.size() < n_clusters){ + while (m_num_left && clusters.size() < n_clusters) { Cluster3x3 c = read_one_cluster(); - if(is_selected(c)){ - clusters.push_back(c.x, c.y, reinterpret_cast(c.data)); + if (is_selected(c)) { + clusters.push_back(c.x, c.y, + reinterpret_cast(c.data)); } } } @@ -172,17 +172,21 @@ ClusterVector ClusterFile::read_clusters_with_cut(size_t n_clusters) { if (clusters.size() < n_clusters) { // sanity check if (m_num_left) { - throw std::runtime_error(LOCATION + "Entered second loop with clusters left\n"); + throw std::runtime_error( + LOCATION + "Entered second loop with clusters left\n"); } - + int32_t frame_number = 0; // frame number needs to be 4 bytes! while (fread(&frame_number, sizeof(frame_number), 1, fp)) { if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) { - clusters.set_frame_number(frame_number); //cluster vector will hold the last frame number - while(m_num_left && clusters.size() < n_clusters){ + clusters.set_frame_number( + frame_number); // cluster vector will hold the last frame + // number + while (m_num_left && clusters.size() < n_clusters) { Cluster3x3 c = read_one_cluster(); - if(is_selected(c)){ - clusters.push_back(c.x, c.y, reinterpret_cast(c.data)); + if (is_selected(c)) { + clusters.push_back( + c.x, c.y, reinterpret_cast(c.data)); } } } @@ -191,15 +195,14 @@ ClusterVector ClusterFile::read_clusters_with_cut(size_t n_clusters) { if (clusters.size() >= n_clusters) break; } - } - if(m_gain_map) + if (m_gain_map) clusters.apply_gain_map(m_gain_map->view()); return clusters; } -Cluster3x3 ClusterFile::read_one_cluster(){ +Cluster3x3 ClusterFile::read_one_cluster() { Cluster3x3 c; auto rc = fread(&c, sizeof(c), 1, fp); if (rc != 1) { @@ -209,13 +212,13 @@ Cluster3x3 ClusterFile::read_one_cluster(){ return c; } -ClusterVector ClusterFile::read_frame(){ +ClusterVector ClusterFile::read_frame() { if (m_mode != "r") { throw std::runtime_error(LOCATION + "File not opened for reading"); } - if (m_noise_map || m_roi){ + if (m_noise_map || m_roi) { return read_frame_with_cut(); - }else{ + } else { return read_frame_without_cut(); } } @@ -235,7 +238,8 @@ ClusterVector ClusterFile::read_frame_without_cut() { int32_t n_clusters; // Saved as 32bit integer in the cluster file if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) { - throw std::runtime_error(LOCATION + "Could not read number of clusters"); + throw std::runtime_error(LOCATION + + "Could not read number of clusters"); } ClusterVector clusters(3, 3, n_clusters); @@ -264,18 +268,17 @@ ClusterVector ClusterFile::read_frame_with_cut() { throw std::runtime_error("Could not read frame number"); } - if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) { throw std::runtime_error("Could not read number of clusters"); } - + ClusterVector clusters(3, 3); clusters.reserve(m_num_left); clusters.set_frame_number(frame_number); - while(m_num_left){ + while (m_num_left) { Cluster3x3 c = read_one_cluster(); - if(is_selected(c)){ - clusters.push_back(c.x, c.y, reinterpret_cast(c.data)); + if (is_selected(c)) { + clusters.push_back(c.x, c.y, reinterpret_cast(c.data)); } } if (m_gain_map) @@ -283,56 +286,56 @@ ClusterVector ClusterFile::read_frame_with_cut() { return clusters; } - - bool ClusterFile::is_selected(Cluster3x3 &cl) { - //Should fail fast + // Should fail fast if (m_roi) { if (!(m_roi->contains(cl.x, cl.y))) { return false; } } - if (m_noise_map){ - int32_t sum_1x1 = cl.data[4]; // central pixel + if (m_noise_map) { + int32_t sum_1x1 = cl.data[4]; // central pixel int32_t sum_2x2 = cl.sum_2x2(); // highest sum of 2x2 subclusters - int32_t sum_3x3 = cl.sum(); // sum of all pixels + int32_t sum_3x3 = cl.sum(); // sum of all pixels - auto noise = (*m_noise_map)(cl.y, cl.x); //TODO! check if this is correct + auto noise = + (*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct if (sum_1x1 <= noise || sum_2x2 <= 2 * noise || sum_3x3 <= 3 * noise) { return false; } } - //we passed all checks + // we passed all checks return true; } NDArray calculate_eta2(ClusterVector &clusters) { - //TOTO! make work with 2x2 clusters + // TOTO! make work with 2x2 clusters NDArray eta2({static_cast(clusters.size()), 2}); - + if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) { for (size_t i = 0; i < clusters.size(); i++) { auto e = calculate_eta2(clusters.at(i)); eta2(i, 0) = e.x; eta2(i, 1) = e.y; } - }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){ + } else if (clusters.cluster_size_x() == 2 || + clusters.cluster_size_y() == 2) { for (size_t i = 0; i < clusters.size(); i++) { auto e = calculate_eta2(clusters.at(i)); eta2(i, 0) = e.x; eta2(i, 1) = e.y; } - }else{ + } else { throw std::runtime_error("Only 3x3 and 2x2 clusters are supported"); } - + return eta2; } -/** - * @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2 struct - * containing etay, etax and the corner of the cluster. -*/ +/** + * @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2 + * struct containing etay, etax and the corner of the cluster. + */ Eta2 calculate_eta2(Cluster3x3 &cl) { Eta2 eta{}; @@ -347,56 +350,46 @@ Eta2 calculate_eta2(Cluster3x3 &cl) { switch (c) { case cBottomLeft: if ((cl.data[3] + cl.data[4]) != 0) - eta.x = - static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); + eta.x = static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); if ((cl.data[1] + cl.data[4]) != 0) - eta.y = - static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); + eta.y = static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); eta.c = cBottomLeft; break; case cBottomRight: if ((cl.data[2] + cl.data[5]) != 0) - eta.x = - static_cast(cl.data[5]) / (cl.data[4] + cl.data[5]); + eta.x = static_cast(cl.data[5]) / (cl.data[4] + cl.data[5]); if ((cl.data[1] + cl.data[4]) != 0) - eta.y = - static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); + eta.y = static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); eta.c = cBottomRight; break; case cTopLeft: if ((cl.data[7] + cl.data[4]) != 0) - eta.x = - static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); + eta.x = static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); if ((cl.data[7] + cl.data[4]) != 0) - eta.y = - static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); + eta.y = static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); eta.c = cTopLeft; break; case cTopRight: if ((cl.data[5] + cl.data[4]) != 0) - eta.x = - static_cast(cl.data[5]) / (cl.data[5] + cl.data[4]); + eta.x = static_cast(cl.data[5]) / (cl.data[5] + cl.data[4]); if ((cl.data[7] + cl.data[4]) != 0) - eta.y = - static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); + eta.y = static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); eta.c = cTopRight; break; - // no default to allow compiler to warn about missing cases + // no default to allow compiler to warn about missing cases } return eta; } - Eta2 calculate_eta2(Cluster2x2 &cl) { Eta2 eta{}; if ((cl.data[0] + cl.data[1]) != 0) eta.x = static_cast(cl.data[1]) / (cl.data[0] + cl.data[1]); if ((cl.data[0] + cl.data[2]) != 0) eta.y = static_cast(cl.data[2]) / (cl.data[0] + cl.data[2]); - eta.sum = cl.data[0] + cl.data[1] + cl.data[2]+ cl.data[3]; - eta.c = cBottomLeft; //TODO! This is not correct, but need to put something + eta.sum = cl.data[0] + cl.data[1] + cl.data[2] + cl.data[3]; + eta.c = cBottomLeft; // TODO! This is not correct, but need to put something return eta; } - } // namespace aare \ No newline at end of file diff --git a/src/ClusterFile.test.cpp b/src/ClusterFile.test.cpp index 6254b5d..68e734d 100644 --- a/src/ClusterFile.test.cpp +++ b/src/ClusterFile.test.cpp @@ -10,9 +10,8 @@ using aare::Cluster; using aare::ClusterFile; using aare::ClusterVector; - TEST_CASE("Read one frame from a cluster file", "[.files]") { - //We know that the frame has 97 clusters + // We know that the frame has 97 clusters auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust"; REQUIRE(std::filesystem::exists(fpath)); @@ -27,7 +26,6 @@ TEST_CASE("Read one frame from a cluster file", "[.files]") { std::begin(expected_cluster_data))); } - TEST_CASE("Read one frame using ROI", "[.files]") { // We know that the frame has 97 clusters auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust"; @@ -60,8 +58,6 @@ TEST_CASE("Read one frame using ROI", "[.files]") { std::begin(expected_cluster_data))); } - - TEST_CASE("Read clusters from single frame file", "[.files]") { // frame_number, num_clusters [135] 97 diff --git a/src/CtbRawFile.cpp b/src/CtbRawFile.cpp index 4d9d895..a6a1d92 100644 --- a/src/CtbRawFile.cpp +++ b/src/CtbRawFile.cpp @@ -14,22 +14,24 @@ CtbRawFile::CtbRawFile(const std::filesystem::path &fname) : m_master(fname) { m_file.open(m_master.data_fname(0, 0), std::ios::binary); } -void CtbRawFile::read_into(std::byte *image_buf, DetectorHeader* header) { - if(m_current_frame >= m_master.frames_in_file()){ +void CtbRawFile::read_into(std::byte *image_buf, DetectorHeader *header) { + if (m_current_frame >= m_master.frames_in_file()) { throw std::runtime_error(LOCATION + " End of file reached"); } - if(m_current_frame != 0 && m_current_frame % m_master.max_frames_per_file() == 0){ - open_data_file(m_current_subfile+1); + if (m_current_frame != 0 && + m_current_frame % m_master.max_frames_per_file() == 0) { + open_data_file(m_current_subfile + 1); } - - if(header){ + + if (header) { m_file.read(reinterpret_cast(header), sizeof(DetectorHeader)); - }else{ + } else { m_file.seekg(sizeof(DetectorHeader), std::ios::cur); } - m_file.read(reinterpret_cast(image_buf), m_master.image_size_in_bytes()); + m_file.read(reinterpret_cast(image_buf), + m_master.image_size_in_bytes()); m_current_frame++; } @@ -38,13 +40,16 @@ void CtbRawFile::seek(size_t frame_number) { open_data_file(index); } size_t frame_number_in_file = frame_number % m_master.max_frames_per_file(); - m_file.seekg((sizeof(DetectorHeader)+m_master.image_size_in_bytes()) * frame_number_in_file); + m_file.seekg((sizeof(DetectorHeader) + m_master.image_size_in_bytes()) * + frame_number_in_file); m_current_frame = frame_number; } size_t CtbRawFile::tell() const { return m_current_frame; } -size_t CtbRawFile::image_size_in_bytes() const { return m_master.image_size_in_bytes(); } +size_t CtbRawFile::image_size_in_bytes() const { + return m_master.image_size_in_bytes(); +} size_t CtbRawFile::frames_in_file() const { return m_master.frames_in_file(); } @@ -63,12 +68,11 @@ void CtbRawFile::open_data_file(size_t subfile_index) { throw std::runtime_error(LOCATION + "Subfile index out of range"); } m_current_subfile = subfile_index; - m_file = std::ifstream(m_master.data_fname(0, subfile_index), std::ios::binary); // only one module for CTB + m_file = std::ifstream(m_master.data_fname(0, subfile_index), + std::ios::binary); // only one module for CTB if (!m_file.is_open()) { throw std::runtime_error(LOCATION + "Could not open data file"); } } - - } // namespace aare \ No newline at end of file diff --git a/src/Dtype.cpp b/src/Dtype.cpp index b818ea3..0fdffec 100644 --- a/src/Dtype.cpp +++ b/src/Dtype.cpp @@ -10,7 +10,8 @@ namespace aare { * @brief Construct a DType object from a type_info object * @param t type_info object * @throw runtime_error if the type is not supported - * @note supported types are: int8_t, uint8_t, int16_t, uint16_t, int32_t, uint32_t, int64_t, uint64_t, float, double + * @note supported types are: int8_t, uint8_t, int16_t, uint16_t, int32_t, + * uint32_t, int64_t, uint64_t, float, double * @note the type_info object is obtained using typeid (e.g. typeid(int)) */ Dtype::Dtype(const std::type_info &t) { @@ -35,7 +36,8 @@ Dtype::Dtype(const std::type_info &t) { else if (t == typeid(double)) m_type = TypeIndex::DOUBLE; else - throw std::runtime_error("Could not construct data type. Type not supported."); + throw std::runtime_error( + "Could not construct data type. Type not supported."); } /** @@ -63,7 +65,8 @@ uint8_t Dtype::bitdepth() const { case TypeIndex::NONE: return 0; default: - throw std::runtime_error(LOCATION + "Could not get bitdepth. Type not supported."); + throw std::runtime_error(LOCATION + + "Could not get bitdepth. Type not supported."); } } @@ -138,7 +141,8 @@ Dtype Dtype::from_bitdepth(uint8_t bitdepth) { case 64: return Dtype(TypeIndex::UINT64); default: - throw std::runtime_error("Could not construct data type from bitdepth."); + throw std::runtime_error( + "Could not construct data type from bitdepth."); } } /** @@ -175,17 +179,27 @@ std::string Dtype::to_string() const { case TypeIndex::DOUBLE: return "f8"; case TypeIndex::ERROR: - throw std::runtime_error("Could not get string representation. Type not supported."); + throw std::runtime_error( + "Could not get string representation. Type not supported."); case TypeIndex::NONE: - throw std::runtime_error("Could not get string representation. Type not supported."); + throw std::runtime_error( + "Could not get string representation. Type not supported."); } return {}; } -bool Dtype::operator==(const Dtype &other) const noexcept { return m_type == other.m_type; } -bool Dtype::operator!=(const Dtype &other) const noexcept { return !(*this == other); } +bool Dtype::operator==(const Dtype &other) const noexcept { + return m_type == other.m_type; +} +bool Dtype::operator!=(const Dtype &other) const noexcept { + return !(*this == other); +} -bool Dtype::operator==(const std::type_info &t) const { return Dtype(t) == *this; } -bool Dtype::operator!=(const std::type_info &t) const { return Dtype(t) != *this; } +bool Dtype::operator==(const std::type_info &t) const { + return Dtype(t) == *this; +} +bool Dtype::operator!=(const std::type_info &t) const { + return Dtype(t) != *this; +} } // namespace aare diff --git a/src/Dtype.test.cpp b/src/Dtype.test.cpp index b252267..256be64 100644 --- a/src/Dtype.test.cpp +++ b/src/Dtype.test.cpp @@ -51,4 +51,6 @@ TEST_CASE("Construct from string with endianess") { REQUIRE_THROWS(Dtype(">i4") == typeid(int32_t)); } -TEST_CASE("Convert to string") { REQUIRE(Dtype(typeid(int)).to_string() == "(fname, mode); - } - else if (fname.extension() == ".npy") { + } 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"){ + } else if (fname.extension() == ".dat") { file_impl = std::make_unique(fname); } else { throw std::runtime_error("Unsupported file type"); } } +File::File(File &&other) noexcept { std::swap(file_impl, other.file_impl); } -File::File(File &&other) noexcept{ - std::swap(file_impl, other.file_impl); -} - -File& File::operator=(File &&other) noexcept { +File &File::operator=(File &&other) noexcept { if (this != &other) { File tmp(std::move(other)); std::swap(file_impl, tmp.file_impl); @@ -70,15 +66,16 @@ size_t File::frame_number(size_t frame_index) { } size_t File::bytes_per_frame() const { return file_impl->bytes_per_frame(); } -size_t File::pixels_per_frame() const{ return file_impl->pixels_per_frame(); } +size_t File::pixels_per_frame() const { return file_impl->pixels_per_frame(); } void File::seek(size_t frame_index) { file_impl->seek(frame_index); } size_t File::tell() const { return file_impl->tell(); } size_t File::rows() const { return file_impl->rows(); } size_t File::cols() const { return file_impl->cols(); } size_t File::bitdepth() const { return file_impl->bitdepth(); } -size_t File::bytes_per_pixel() const { return file_impl->bitdepth() / bits_per_byte; } +size_t File::bytes_per_pixel() const { + return file_impl->bitdepth() / bits_per_byte; +} DetectorType File::detector_type() const { return file_impl->detector_type(); } - } // namespace aare \ No newline at end of file diff --git a/src/FilePtr.cpp b/src/FilePtr.cpp index e3cdb4b..f850080 100644 --- a/src/FilePtr.cpp +++ b/src/FilePtr.cpp @@ -6,10 +6,12 @@ namespace aare { -FilePtr::FilePtr(const std::filesystem::path& fname, const std::string& mode = "rb") { +FilePtr::FilePtr(const std::filesystem::path &fname, + const std::string &mode = "rb") { fp_ = fopen(fname.c_str(), mode.c_str()); if (!fp_) - throw std::runtime_error(fmt::format("Could not open: {}", fname.c_str())); + throw std::runtime_error( + fmt::format("Could not open: {}", fname.c_str())); } FilePtr::FilePtr(FilePtr &&other) { std::swap(fp_, other.fp_); } @@ -24,15 +26,16 @@ FILE *FilePtr::get() { return fp_; } ssize_t FilePtr::tell() { auto pos = ftell(fp_); if (pos == -1) - throw std::runtime_error(fmt::format("Error getting file position: {}", error_msg())); + throw std::runtime_error( + fmt::format("Error getting file position: {}", error_msg())); return pos; -} +} FilePtr::~FilePtr() { if (fp_) fclose(fp_); // check? } -std::string FilePtr::error_msg(){ +std::string FilePtr::error_msg() { if (feof(fp_)) { return "End of file reached"; } diff --git a/src/Fit.cpp b/src/Fit.cpp index 25000de..9d4b70b 100644 --- a/src/Fit.cpp +++ b/src/Fit.cpp @@ -1,13 +1,12 @@ #include "aare/Fit.hpp" -#include "aare/utils/task.hpp" #include "aare/utils/par.hpp" +#include "aare/utils/task.hpp" #include #include #include #include - namespace aare { namespace func { @@ -34,8 +33,10 @@ NDArray pol1(NDView x, NDView par) { return y; } -double scurve(const double x, const double * par) { - return (par[0] + par[1] * x) + 0.5 * (1 + erf((x - par[2]) / (sqrt(2) * par[3]))) * (par[4] + par[5] * (x - par[2])); +double scurve(const double x, const double *par) { + return (par[0] + par[1] * x) + + 0.5 * (1 + erf((x - par[2]) / (sqrt(2) * par[3]))) * + (par[4] + par[5] * (x - par[2])); } NDArray scurve(NDView x, NDView par) { @@ -46,8 +47,10 @@ NDArray scurve(NDView x, NDView par) { return y; } -double scurve2(const double x, const double * par) { - return (par[0] + par[1] * x) + 0.5 * (1 - erf((x - par[2]) / (sqrt(2) * par[3]))) * (par[4] + par[5] * (x - par[2])); +double scurve2(const double x, const double *par) { + return (par[0] + par[1] * x) + + 0.5 * (1 - erf((x - par[2]) / (sqrt(2) * par[3]))) * + (par[4] + par[5] * (x - par[2])); } NDArray scurve2(NDView x, NDView par) { @@ -91,7 +94,8 @@ NDArray fit_gaus(NDView x, NDView y, return result; } -std::array gaus_init_par(const NDView x, const NDView y) { +std::array gaus_init_par(const NDView x, + const NDView y) { std::array start_par{0, 0, 0}; auto e = std::max_element(y.begin(), y.end()); auto idx = std::distance(y.begin(), e); @@ -103,31 +107,29 @@ std::array gaus_init_par(const NDView x, const NDView *e / 2; }) * - delta / 2.35; + start_par[2] = std::count_if(y.begin(), y.end(), + [e](double val) { return val > *e / 2; }) * + delta / 2.35; return start_par; } +std::array pol1_init_par(const NDView x, + const NDView y) { + // Estimate the initial parameters for the fit + std::array start_par{0, 0}; -std::array pol1_init_par(const NDView x, const NDView y){ - // Estimate the initial parameters for the fit - std::array start_par{0, 0}; + auto y2 = std::max_element(y.begin(), y.end()); + auto x2 = x[std::distance(y.begin(), y2)]; + auto y1 = std::min_element(y.begin(), y.end()); + auto x1 = x[std::distance(y.begin(), y1)]; - - auto y2 = std::max_element(y.begin(), y.end()); - auto x2 = x[std::distance(y.begin(), y2)]; - auto y1 = std::min_element(y.begin(), y.end()); - auto x1 = x[std::distance(y.begin(), y1)]; - - start_par[0] = - (*y2 - *y1) / (x2 - x1); // For amplitude we use the maximum value - start_par[1] = - *y1 - ((*y2 - *y1) / (x2 - x1)) * - x1; // For the mean we use the x value of the maximum value - return start_par; + start_par[0] = + (*y2 - *y1) / (x2 - x1); // For amplitude we use the maximum value + start_par[1] = + *y1 - ((*y2 - *y1) / (x2 - x1)) * + x1; // For the mean we use the x value of the maximum value + return start_par; } void fit_gaus(NDView x, NDView y, NDView y_err, @@ -141,7 +143,6 @@ void fit_gaus(NDView x, NDView y, NDView y_err, "and par_out, par_err_out must have size 3"); } - // /* Collection of output parameters for status info. */ // typedef struct { // double fnorm; /* norm of the residue vector fvec. */ @@ -153,23 +154,32 @@ void fit_gaus(NDView x, NDView y, NDView y_err, // */ // } lm_status_struct; - lm_status_struct status; par_out = gaus_init_par(x, y); - std::array cov{0, 0, 0, 0, 0, 0, 0 , 0 , 0}; + std::array cov{0, 0, 0, 0, 0, 0, 0, 0, 0}; - // void lmcurve2( const int n_par, double *par, double *parerr, double *covar, const int m_dat, const double *t, const double *y, const double *dy, double (*f)( const double ti, const double *par ), const lm_control_struct *control, lm_status_struct *status); - // n_par - Number of free variables. Length of parameter vector par. - // par - Parameter vector. On input, it must contain a reasonable guess. On output, it contains the solution found to minimize ||r||. - // parerr - Parameter uncertainties vector. Array of length n_par or NULL. On output, unless it or covar is NULL, it contains the weighted parameter uncertainties for the found parameters. - // covar - Covariance matrix. Array of length n_par * n_par or NULL. On output, unless it is NULL, it contains the covariance matrix. - // m_dat - Number of data points. Length of vectors t, y, dy. Must statisfy n_par <= m_dat. - // t - Array of length m_dat. Contains the abcissae (time, or "x") for which function f will be evaluated. - // y - Array of length m_dat. Contains the ordinate values that shall be fitted. - // dy - Array of length m_dat. Contains the standard deviations of the values y. - // f - A user-supplied parametric function f(ti;par). - // control - Parameter collection for tuning the fit procedure. In most cases, the default &lm_control_double is adequate. If f is only computed with single-precision accuracy, &lm_control_float should be used. Parameters are explained in lmmin2(3). - // status - A record used to return information about the minimization process: For details, see lmmin2(3). + // void lmcurve2( const int n_par, double *par, double *parerr, double + // *covar, const int m_dat, const double *t, const double *y, const double + // *dy, double (*f)( const double ti, const double *par ), const + // lm_control_struct *control, lm_status_struct *status); n_par - Number of + // free variables. Length of parameter vector par. par - Parameter vector. + // On input, it must contain a reasonable guess. On output, it contains the + // solution found to minimize ||r||. parerr - Parameter uncertainties + // vector. Array of length n_par or NULL. On output, unless it or covar is + // NULL, it contains the weighted parameter uncertainties for the found + // parameters. covar - Covariance matrix. Array of length n_par * n_par or + // NULL. On output, unless it is NULL, it contains the covariance matrix. + // m_dat - Number of data points. Length of vectors t, y, dy. Must statisfy + // n_par <= m_dat. t - Array of length m_dat. Contains the abcissae (time, + // or "x") for which function f will be evaluated. y - Array of length + // m_dat. Contains the ordinate values that shall be fitted. dy - Array of + // length m_dat. Contains the standard deviations of the values y. f - A + // user-supplied parametric function f(ti;par). control - Parameter + // collection for tuning the fit procedure. In most cases, the default + // &lm_control_double is adequate. If f is only computed with + // single-precision accuracy, &lm_control_float should be used. Parameters + // are explained in lmmin2(3). status - A record used to return information + // about the minimization process: For details, see lmmin2(3). lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(), x.size(), x.data(), y.data(), y_err.data(), aare::func::gaus, @@ -178,12 +188,14 @@ void fit_gaus(NDView x, NDView y, NDView y_err, // Calculate chi2 chi2 = 0; for (ssize_t i = 0; i < y.size(); i++) { - chi2 += std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2); + chi2 += + std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2); } } void fit_gaus(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, + NDView par_out, NDView par_err_out, + NDView chi2_out, int n_threads) { @@ -197,10 +209,9 @@ void fit_gaus(NDView x, NDView y, NDView y_err, {par_out.shape(2)}); NDView par_err_out_view(&par_err_out(row, col, 0), {par_err_out.shape(2)}); - + fit_gaus(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col)); - } } }; @@ -210,7 +221,8 @@ void fit_gaus(NDView x, NDView y, NDView y_err, } void fit_pol1(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2) { + NDView par_out, NDView par_err_out, + double &chi2) { // Check that we have the correct sizes if (y.size() != x.size() || y.size() != y_err.size() || @@ -230,13 +242,14 @@ void fit_pol1(NDView x, NDView y, NDView y_err, // Calculate chi2 chi2 = 0; for (ssize_t i = 0; i < y.size(); i++) { - chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2); + chi2 += + std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2); } } void fit_pol1(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads) { + NDView par_out, NDView par_err_out, + NDView chi2_out, int n_threads) { auto process = [&](ssize_t first_row, ssize_t last_row) { for (ssize_t row = first_row; row < last_row; row++) { @@ -249,15 +262,14 @@ void fit_pol1(NDView x, NDView y, NDView y_err, NDView par_err_out_view(&par_err_out(row, col, 0), {par_err_out.shape(2)}); - fit_pol1(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col)); - + fit_pol1(x, y_view, y_err_view, par_out_view, par_err_out_view, + chi2_out(row, col)); } } }; auto tasks = split_task(0, y.shape(0), n_threads); RunInParallel(process, tasks); - } NDArray fit_pol1(NDView x, NDView y) { @@ -300,27 +312,29 @@ NDArray fit_pol1(NDView x, NDView y, // ~~ S-CURVES ~~ // SCURVE -- -std::array scurve_init_par(const NDView x, const NDView y){ - // Estimate the initial parameters for the fit - std::array start_par{0, 0, 0, 0, 0, 0}; +std::array scurve_init_par(const NDView x, + const NDView y) { + // Estimate the initial parameters for the fit + std::array start_par{0, 0, 0, 0, 0, 0}; - auto ymax = std::max_element(y.begin(), y.end()); - auto ymin = std::min_element(y.begin(), y.end()); - start_par[4] = *ymin + (*ymax - *ymin) / 2; - - // Find the first x where the corresponding y value is above the threshold (start_par[4]) - for (ssize_t i = 0; i < y.size(); ++i) { - if (y[i] >= start_par[4]) { - start_par[2] = x[i]; - break; // Exit the loop after finding the first valid x - } + auto ymax = std::max_element(y.begin(), y.end()); + auto ymin = std::min_element(y.begin(), y.end()); + start_par[4] = *ymin + (*ymax - *ymin) / 2; + + // Find the first x where the corresponding y value is above the threshold + // (start_par[4]) + for (ssize_t i = 0; i < y.size(); ++i) { + if (y[i] >= start_par[4]) { + start_par[2] = x[i]; + break; // Exit the loop after finding the first valid x } + } - start_par[3] = 2 * sqrt(start_par[2]); - start_par[0] = 100; - start_par[1] = 0.25; - start_par[5] = 1; - return start_par; + start_par[3] = 2 * sqrt(start_par[2]); + start_par[0] = 100; + start_par[1] = 0.25; + start_par[5] = 1; + return start_par; } // - No error @@ -334,7 +348,8 @@ NDArray fit_scurve(NDView x, NDView y) { return result; } -NDArray fit_scurve(NDView x, NDView y, int n_threads) { +NDArray fit_scurve(NDView x, NDView y, + int n_threads) { NDArray result({y.shape(0), y.shape(1), 6}, 0); auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) { @@ -358,8 +373,9 @@ NDArray fit_scurve(NDView x, NDView y, int n_th } // - Error -void fit_scurve(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2) { +void fit_scurve(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, double &chi2) { // Check that we have the correct sizes if (y.size() != x.size() || y.size() != y_err.size() || @@ -380,13 +396,15 @@ void fit_scurve(NDView x, NDView y, NDView y_er // Calculate chi2 chi2 = 0; for (ssize_t i = 0; i < y.size(); i++) { - chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2); + chi2 += + std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2); } } -void fit_scurve(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads) { +void fit_scurve(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, NDView chi2_out, + int n_threads) { auto process = [&](ssize_t first_row, ssize_t last_row) { for (ssize_t row = first_row; row < last_row; row++) { @@ -399,40 +417,41 @@ void fit_scurve(NDView x, NDView y, NDView y_er NDView par_err_out_view(&par_err_out(row, col, 0), {par_err_out.shape(2)}); - fit_scurve(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col)); - + fit_scurve(x, y_view, y_err_view, par_out_view, + par_err_out_view, chi2_out(row, col)); } } }; auto tasks = split_task(0, y.shape(0), n_threads); RunInParallel(process, tasks); - } // SCURVE2 --- -std::array scurve2_init_par(const NDView x, const NDView y){ - // Estimate the initial parameters for the fit - std::array start_par{0, 0, 0, 0, 0, 0}; +std::array scurve2_init_par(const NDView x, + const NDView y) { + // Estimate the initial parameters for the fit + std::array start_par{0, 0, 0, 0, 0, 0}; - auto ymax = std::max_element(y.begin(), y.end()); - auto ymin = std::min_element(y.begin(), y.end()); - start_par[4] = *ymin + (*ymax - *ymin) / 2; - - // Find the first x where the corresponding y value is above the threshold (start_par[4]) - for (ssize_t i = 0; i < y.size(); ++i) { - if (y[i] <= start_par[4]) { - start_par[2] = x[i]; - break; // Exit the loop after finding the first valid x - } + auto ymax = std::max_element(y.begin(), y.end()); + auto ymin = std::min_element(y.begin(), y.end()); + start_par[4] = *ymin + (*ymax - *ymin) / 2; + + // Find the first x where the corresponding y value is above the threshold + // (start_par[4]) + for (ssize_t i = 0; i < y.size(); ++i) { + if (y[i] <= start_par[4]) { + start_par[2] = x[i]; + break; // Exit the loop after finding the first valid x } + } - start_par[3] = 2 * sqrt(start_par[2]); - start_par[0] = 100; - start_par[1] = 0.25; - start_par[5] = -1; - return start_par; + start_par[3] = 2 * sqrt(start_par[2]); + start_par[0] = 100; + start_par[1] = 0.25; + start_par[5] = -1; + return start_par; } // - No error @@ -446,7 +465,8 @@ NDArray fit_scurve2(NDView x, NDView y) { return result; } -NDArray fit_scurve2(NDView x, NDView y, int n_threads) { +NDArray fit_scurve2(NDView x, NDView y, + int n_threads) { NDArray result({y.shape(0), y.shape(1), 6}, 0); auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) { @@ -470,8 +490,9 @@ NDArray fit_scurve2(NDView x, NDView y, int n_t } // - Error -void fit_scurve2(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2) { +void fit_scurve2(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, double &chi2) { // Check that we have the correct sizes if (y.size() != x.size() || y.size() != y_err.size() || @@ -492,13 +513,15 @@ void fit_scurve2(NDView x, NDView y, NDView y_e // Calculate chi2 chi2 = 0; for (ssize_t i = 0; i < y.size(); i++) { - chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2); + chi2 += + std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2); } } -void fit_scurve2(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads) { +void fit_scurve2(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, NDView chi2_out, + int n_threads) { auto process = [&](ssize_t first_row, ssize_t last_row) { for (ssize_t row = first_row; row < last_row; row++) { @@ -511,15 +534,14 @@ void fit_scurve2(NDView x, NDView y, NDView y_e NDView par_err_out_view(&par_err_out(row, col, 0), {par_err_out.shape(2)}); - fit_scurve2(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col)); - + fit_scurve2(x, y_view, y_err_view, par_out_view, + par_err_out_view, chi2_out(row, col)); } } }; auto tasks = split_task(0, y.shape(0), n_threads); RunInParallel(process, tasks); - } } // namespace aare \ No newline at end of file diff --git a/src/Frame.cpp b/src/Frame.cpp index d44bed5..ef7675f 100644 --- a/src/Frame.cpp +++ b/src/Frame.cpp @@ -29,8 +29,7 @@ uint64_t Frame::size() const { return m_rows * m_cols; } size_t Frame::bytes() const { return m_rows * m_cols * m_dtype.bytes(); } std::byte *Frame::data() const { return m_data; } - -std::byte *Frame::pixel_ptr(uint32_t row, uint32_t col) const{ +std::byte *Frame::pixel_ptr(uint32_t row, uint32_t col) const { if ((row >= m_rows) || (col >= m_cols)) { std::cerr << "Invalid row or column index" << '\n'; return nullptr; @@ -38,7 +37,6 @@ std::byte *Frame::pixel_ptr(uint32_t row, uint32_t col) const{ return m_data + (row * m_cols + col) * (m_dtype.bytes()); } - Frame &Frame::operator=(Frame &&other) noexcept { if (this == &other) { return *this; @@ -70,5 +68,4 @@ Frame Frame::clone() const { return frame; } - } // namespace aare diff --git a/src/Frame.test.cpp b/src/Frame.test.cpp index 4063701..bafb39c 100644 --- a/src/Frame.test.cpp +++ b/src/Frame.test.cpp @@ -65,7 +65,8 @@ TEST_CASE("Set a value in a 64 bit frame") { // only the value we did set should be non-zero for (size_t i = 0; i < rows; i++) { for (size_t j = 0; j < cols; j++) { - uint64_t *data = reinterpret_cast(frame.pixel_ptr(i, j)); + uint64_t *data = + reinterpret_cast(frame.pixel_ptr(i, j)); REQUIRE(data != nullptr); if (i == 5 && j == 7) { REQUIRE(*data == value); @@ -150,4 +151,3 @@ TEST_CASE("test explicit copy constructor") { REQUIRE(frame2.bytes() == rows * cols * bitdepth / 8); REQUIRE(frame2.data() != data); } - diff --git a/src/JungfrauDataFile.cpp b/src/JungfrauDataFile.cpp index 59a1a0a..5fb99a6 100644 --- a/src/JungfrauDataFile.cpp +++ b/src/JungfrauDataFile.cpp @@ -19,16 +19,15 @@ JungfrauDataFile::JungfrauDataFile(const std::filesystem::path &fname) { open_file(m_current_file_index); } - // FileInterface -Frame JungfrauDataFile::read_frame(){ +Frame JungfrauDataFile::read_frame() { Frame f(rows(), cols(), Dtype::UINT16); read_into(reinterpret_cast(f.data()), nullptr); return f; } -Frame JungfrauDataFile::read_frame(size_t frame_number){ +Frame JungfrauDataFile::read_frame(size_t frame_number) { seek(frame_number); Frame f(rows(), cols(), Dtype::UINT16); read_into(reinterpret_cast(f.data()), nullptr); @@ -37,7 +36,7 @@ Frame JungfrauDataFile::read_frame(size_t frame_number){ std::vector JungfrauDataFile::read_n(size_t n_frames) { std::vector frames; - for(size_t i = 0; i < n_frames; ++i){ + for (size_t i = 0; i < n_frames; ++i) { frames.push_back(read_frame()); } return frames; @@ -48,7 +47,7 @@ void JungfrauDataFile::read_into(std::byte *image_buf) { } 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); @@ -59,7 +58,9 @@ std::array JungfrauDataFile::shape() const { return {static_cast(rows()), static_cast(cols())}; } -DetectorType JungfrauDataFile::detector_type() const { return DetectorType::Jungfrau; } +DetectorType JungfrauDataFile::detector_type() const { + return DetectorType::Jungfrau; +} std::string JungfrauDataFile::base_name() const { return m_base_name; } @@ -195,22 +196,23 @@ void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames, JungfrauDataHeader *header) { if (header) { for (size_t i = 0; i < n_frames; ++i) - read_into(image_buf + i * m_bytes_per_frame, header + i); - }else{ + 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(image->shape()!=shape()){ - throw std::runtime_error(LOCATION + - "Image shape does not match file size: " + std::to_string(rows()) + "x" + std::to_string(cols())); +void JungfrauDataFile::read_into(NDArray *image, + JungfrauDataHeader *header) { + if (image->shape() != shape()) { + 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); } - JungfrauDataHeader JungfrauDataFile::read_header() { JungfrauDataHeader header; if (auto rc = fread(&header, 1, sizeof(header), m_fp.get()); diff --git a/src/JungfrauDataFile.test.cpp b/src/JungfrauDataFile.test.cpp index ce51168..21a4e32 100644 --- a/src/JungfrauDataFile.test.cpp +++ b/src/JungfrauDataFile.test.cpp @@ -1,21 +1,21 @@ #include "aare/JungfrauDataFile.hpp" -#include #include "test_config.hpp" +#include using aare::JungfrauDataFile; using aare::JungfrauDataHeader; TEST_CASE("Open a Jungfrau data file", "[.files]") { - //we know we have 4 files with 7, 7, 7, and 3 frames - //firs frame number if 1 and the bunch id is frame_number**2 - //so we can check the header + // we know we have 4 files with 7, 7, 7, and 3 frames + // firs frame number if 1 and the bunch id is frame_number**2 + // so we can check the header auto fpath = test_data_path() / "dat" / "AldoJF500k_000000.dat"; REQUIRE(std::filesystem::exists(fpath)); JungfrauDataFile f(fpath); REQUIRE(f.rows() == 512); REQUIRE(f.cols() == 1024); - REQUIRE(f.bytes_per_frame() == 1048576); + REQUIRE(f.bytes_per_frame() == 1048576); REQUIRE(f.pixels_per_frame() == 524288); REQUIRE(f.bytes_per_pixel() == 2); REQUIRE(f.bitdepth() == 16); @@ -25,7 +25,7 @@ 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 + // Check that the frame number and buch id is read correctly for (size_t i = 0; i < 24; ++i) { JungfrauDataHeader header; aare::NDArray image(f.shape()); @@ -37,65 +37,64 @@ TEST_CASE("Open a Jungfrau data file", "[.files]") { } } -TEST_CASE("Seek in a JungfrauDataFile", "[.files]"){ +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 + // The file should have 113 frames f.seek(19); REQUIRE(f.tell() == 19); auto h = f.read_header(); - REQUIRE(h.framenum == 19+1); + REQUIRE(h.framenum == 19 + 1); - //Reading again does not change the file pointer + // Reading again does not change the file pointer auto h2 = f.read_header(); - REQUIRE(h2.framenum == 19+1); + REQUIRE(h2.framenum == 19 + 1); f.seek(59); REQUIRE(f.tell() == 59); auto h3 = f.read_header(); - REQUIRE(h3.framenum == 59+1); + REQUIRE(h3.framenum == 59 + 1); JungfrauDataHeader h4; aare::NDArray image(f.shape()); f.read_into(&image, &h4); - REQUIRE(h4.framenum == 59+1); + REQUIRE(h4.framenum == 59 + 1); - //now we should be on the next frame + // now we should be on the next frame REQUIRE(f.tell() == 60); - REQUIRE(f.read_header().framenum == 60+1); + REQUIRE(f.read_header().framenum == 60 + 1); - REQUIRE_THROWS(f.seek(86356)); //out of range + REQUIRE_THROWS(f.seek(86356)); // out of range } -TEST_CASE("Open a Jungfrau data file with non zero file index", "[.files]"){ +TEST_CASE("Open a Jungfrau data file with non zero file index", "[.files]") { auto fpath = test_data_path() / "dat" / "AldoJF65k_000003.dat"; REQUIRE(std::filesystem::exists(fpath)); JungfrauDataFile f(fpath); - //18 files per data file, opening the 3rd file we ignore the first 3 - REQUIRE(f.total_frames() == 113-18*3); + // 18 files per data file, opening the 3rd file we ignore the first 3 + REQUIRE(f.total_frames() == 113 - 18 * 3); REQUIRE(f.tell() == 0); - //Frame numbers start at 1 in the first file - REQUIRE(f.read_header().framenum == 18*3+1); + // Frame numbers start at 1 in the first file + REQUIRE(f.read_header().framenum == 18 * 3 + 1); // moving relative to the third file f.seek(5); - REQUIRE(f.read_header().framenum == 18*3+1+5); + REQUIRE(f.read_header().framenum == 18 * 3 + 1 + 5); // ignoring the first 3 files REQUIRE(f.n_files() == 4); REQUIRE(f.current_file().stem() == "AldoJF65k_000003"); - } -TEST_CASE("Read into throws if size doesn't match", "[.files]"){ +TEST_CASE("Read into throws if size doesn't match", "[.files]") { auto fpath = test_data_path() / "dat" / "AldoJF65k_000000.dat"; REQUIRE(std::filesystem::exists(fpath)); @@ -109,6 +108,4 @@ TEST_CASE("Read into throws if size doesn't match", "[.files]"){ REQUIRE_THROWS(f.read_into(&image)); REQUIRE(f.tell() == 0); - - } \ No newline at end of file diff --git a/src/NDArray.test.cpp b/src/NDArray.test.cpp index 819a1a9..91b5933 100644 --- a/src/NDArray.test.cpp +++ b/src/NDArray.test.cpp @@ -35,7 +35,7 @@ TEST_CASE("Construct from an NDView") { } } -TEST_CASE("3D NDArray from NDView"){ +TEST_CASE("3D NDArray from NDView") { std::vector data(27); std::iota(data.begin(), data.end(), 0); NDView view(data.data(), Shape<3>{3, 3, 3}); @@ -44,9 +44,9 @@ TEST_CASE("3D NDArray from NDView"){ REQUIRE(image.size() == view.size()); REQUIRE(image.data() != view.data()); - for(ssize_t i=0; i MultiplyNDArrayUsingOperator(NDArray &a, NDArray &b) { // return a * a * b * b; - NDArrayc = a*b; + NDArray c = a * b; return c; } @@ -162,7 +162,6 @@ NDArray AddNDArrayUsingIndex(NDArray &a, NDArray &b) { return res; } - TEST_CASE("Compare two images") { NDArray a; NDArray b; @@ -222,7 +221,6 @@ TEST_CASE("Bitwise and on data") { REQUIRE(a(2) == 384); } - TEST_CASE("Elementwise operations on images") { std::array shape{5, 5}; double a_val = 3.0; @@ -258,7 +256,8 @@ TEST_CASE("Elementwise operations on images") { NDArray A(shape, a_val); NDArray B(shape, b_val); NDArray C = A - B; - // auto C = A - B; // This works but the result is a lazy ArraySub object + // auto C = A - B; // This works but the result is a lazy ArraySub + // object // Value of C matches for (uint32_t i = 0; i < C.size(); ++i) { @@ -282,7 +281,8 @@ TEST_CASE("Elementwise operations on images") { SECTION("Multiply two images") { NDArray A(shape, a_val); NDArray B(shape, b_val); - // auto C = A * B; // This works but the result is a lazy ArrayMul object + // auto C = A * B; // This works but the result is a lazy ArrayMul + // object NDArray C = A * B; // Value of C matches @@ -307,7 +307,8 @@ TEST_CASE("Elementwise operations on images") { SECTION("Divide two images") { NDArray A(shape, a_val); NDArray B(shape, b_val); - // auto C = A / B; // This works but the result is a lazy ArrayDiv object + // auto C = A / B; // This works but the result is a lazy ArrayDiv + // object NDArray C = A / B; // Value of C matches diff --git a/src/NDView.test.cpp b/src/NDView.test.cpp index 89e76e9..a65758c 100644 --- a/src/NDView.test.cpp +++ b/src/NDView.test.cpp @@ -2,8 +2,8 @@ #include #include -#include #include +#include using aare::NDView; using aare::Shape; @@ -151,8 +151,10 @@ TEST_CASE("divide with another span") { std::vector vec1{3, 2, 1}; std::vector result{3, 6, 3}; - NDView data0(vec0.data(), Shape<1>{static_cast(vec0.size())}); - NDView data1(vec1.data(), Shape<1>{static_cast(vec1.size())}); + NDView data0(vec0.data(), + Shape<1>{static_cast(vec0.size())}); + NDView data1(vec1.data(), + Shape<1>{static_cast(vec1.size())}); data0 /= data1; @@ -181,8 +183,7 @@ TEST_CASE("compare two views") { REQUIRE((view1 == view2)); } - -TEST_CASE("Create a view over a vector"){ +TEST_CASE("Create a view over a vector") { std::vector vec(12); std::iota(vec.begin(), vec.end(), 0); auto v = aare::make_view(vec); diff --git a/src/NumpyFile.cpp b/src/NumpyFile.cpp index e375ce3..4e0c215 100644 --- a/src/NumpyFile.cpp +++ b/src/NumpyFile.cpp @@ -4,16 +4,16 @@ namespace aare { - - -NumpyFile::NumpyFile(const std::filesystem::path &fname, const std::string &mode, FileConfig cfg) { +NumpyFile::NumpyFile(const std::filesystem::path &fname, + const std::string &mode, FileConfig cfg) { // TODO! add opts to constructor m_mode = mode; if (mode == "r") { fp = fopen(fname.string().c_str(), "rb"); if (!fp) { - throw std::runtime_error(fmt::format("Could not open: {} for reading", fname.string())); + throw std::runtime_error( + fmt::format("Could not open: {} for reading", fname.string())); } load_metadata(); } else if (mode == "w") { @@ -24,11 +24,15 @@ NumpyFile::NumpyFile(const std::filesystem::path &fname, const std::string &mode m_header.shape = {0, cfg.rows, cfg.cols}; fp = fopen(fname.string().c_str(), "wb"); if (!fp) { - throw std::runtime_error(fmt::format("Could not open: {} for reading", fname.string())); + throw std::runtime_error( + fmt::format("Could not open: {} for reading", fname.string())); } - initial_header_len = aare::NumpyHelpers::write_header(std::filesystem::path(fname.c_str()), m_header); + initial_header_len = aare::NumpyHelpers::write_header( + std::filesystem::path(fname.c_str()), m_header); } - m_pixels_per_frame = std::accumulate(m_header.shape.begin() + 1, m_header.shape.end(), 1, std::multiplies<>()); + m_pixels_per_frame = + std::accumulate(m_header.shape.begin() + 1, m_header.shape.end(), 1, + std::multiplies<>()); m_bytes_per_frame = m_header.dtype.bitdepth() / 8 * m_pixels_per_frame; } @@ -63,7 +67,8 @@ void NumpyFile::get_frame_into(size_t frame_number, std::byte *image_buf) { if (frame_number > m_header.shape[0]) { throw std::invalid_argument("Frame number out of range"); } - if (fseek(fp, header_size + frame_number * m_bytes_per_frame, SEEK_SET)) // NOLINT + if (fseek(fp, header_size + frame_number * m_bytes_per_frame, + SEEK_SET)) // NOLINT throw std::runtime_error("Could not seek to frame"); size_t const rc = fread(image_buf, m_bytes_per_frame, 1, fp); @@ -113,7 +118,8 @@ NumpyFile::~NumpyFile() noexcept { // write header size_t const rc = fwrite(header_str.c_str(), header_str.size(), 1, fp); if (rc != 1) { - std::cout << "Error writing header to numpy file in destructor" << std::endl; + std::cout << "Error writing header to numpy file in destructor" + << std::endl; } } @@ -140,8 +146,10 @@ void NumpyFile::load_metadata() { } // read version - rc = fread(reinterpret_cast(&major_ver_), sizeof(major_ver_), 1, fp); - rc += fread(reinterpret_cast(&minor_ver_), sizeof(minor_ver_), 1, fp); + rc = + fread(reinterpret_cast(&major_ver_), sizeof(major_ver_), 1, fp); + rc += + fread(reinterpret_cast(&minor_ver_), sizeof(minor_ver_), 1, fp); if (rc != 2) { throw std::runtime_error("Error reading numpy version"); } @@ -159,7 +167,8 @@ void NumpyFile::load_metadata() { if (rc != 1) { throw std::runtime_error("Error reading header length"); } - header_size = aare::NumpyHelpers::magic_string_length + 2 + header_len_size + header_len; + header_size = aare::NumpyHelpers::magic_string_length + 2 + + header_len_size + header_len; if (header_size % 16 != 0) { fmt::print("Warning: header length is not a multiple of 16\n"); } diff --git a/src/NumpyFile.test.cpp b/src/NumpyFile.test.cpp index bdd6451..e687bea 100644 --- a/src/NumpyFile.test.cpp +++ b/src/NumpyFile.test.cpp @@ -1,8 +1,8 @@ #include "aare/NumpyFile.hpp" #include "aare/NDArray.hpp" -#include #include "test_config.hpp" +#include using aare::Dtype; using aare::NumpyFile; @@ -23,7 +23,7 @@ TEST_CASE("Read a 1D numpy file with int32 data type", "[.integration]") { REQUIRE(data(i) == i); } } - + TEST_CASE("Read a 3D numpy file with np.double data type", "[.integration]") { auto fpath = test_data_path() / "numpy" / "test_3d_double.npy"; diff --git a/src/NumpyHelpers.cpp b/src/NumpyHelpers.cpp index b8414d9..106327e 100644 --- a/src/NumpyHelpers.cpp +++ b/src/NumpyHelpers.cpp @@ -29,7 +29,8 @@ namespace aare { std::string NumpyHeader::to_string() const { std::stringstream sstm; - sstm << "dtype: " << dtype.to_string() << ", fortran_order: " << fortran_order << ' '; + sstm << "dtype: " << dtype.to_string() + << ", fortran_order: " << fortran_order << ' '; sstm << "shape: ("; for (auto item : shape) sstm << item << ','; @@ -37,10 +38,10 @@ std::string NumpyHeader::to_string() const { return sstm.str(); } - namespace NumpyHelpers { -std::unordered_map parse_dict(std::string in, const std::vector &keys) { +std::unordered_map +parse_dict(std::string in, const std::vector &keys) { std::unordered_map map; if (keys.empty()) return map; @@ -100,7 +101,8 @@ aare::Dtype parse_descr(std::string typestring) { constexpr char little_endian_char = '<'; constexpr char big_endian_char = '>'; constexpr char no_endian_char = '|'; - constexpr std::array endian_chars = {little_endian_char, big_endian_char, no_endian_char}; + constexpr std::array endian_chars = { + little_endian_char, big_endian_char, no_endian_char}; constexpr std::array numtype_chars = {'f', 'i', 'u', 'c'}; const char byteorder_c = typestring[0]; @@ -139,7 +141,9 @@ std::string get_value_from_map(const std::string &mapstr) { return trim(tmp); } -bool is_digits(const std::string &str) { return std::all_of(str.begin(), str.end(), ::isdigit); } +bool is_digits(const std::string &str) { + return std::all_of(str.begin(), str.end(), ::isdigit); +} std::vector parse_tuple(std::string in) { std::vector v; @@ -215,20 +219,25 @@ inline std::string write_boolean(bool b) { return "False"; } -inline std::string write_header_dict(const std::string &descr, bool fortran_order, const std::vector &shape) { +inline std::string write_header_dict(const std::string &descr, + bool fortran_order, + const std::vector &shape) { std::string const s_fortran_order = write_boolean(fortran_order); std::string const shape_s = write_tuple(shape); - return "{'descr': '" + descr + "', 'fortran_order': " + s_fortran_order + ", 'shape': " + shape_s + ", }"; + return "{'descr': '" + descr + "', 'fortran_order': " + s_fortran_order + + ", 'shape': " + shape_s + ", }"; } -size_t write_header(const std::filesystem::path &fname, const NumpyHeader &header) { +size_t write_header(const std::filesystem::path &fname, + const NumpyHeader &header) { std::ofstream out(fname, std::ios::binary | std::ios::out); return write_header(out, header); } size_t write_header(std::ostream &out, const NumpyHeader &header) { - std::string const header_dict = write_header_dict(header.dtype.to_string(), header.fortran_order, header.shape); + std::string const header_dict = write_header_dict( + header.dtype.to_string(), header.fortran_order, header.shape); size_t length = magic_string_length + 2 + 2 + header_dict.length() + 1; @@ -247,17 +256,22 @@ size_t write_header(std::ostream &out, const NumpyHeader &header) { // write header length if (version_major == 1 && version_minor == 0) { - auto header_len = static_cast(header_dict.length() + padding.length() + 1); + auto header_len = + static_cast(header_dict.length() + padding.length() + 1); - std::array header_len_le16{static_cast((header_len >> 0) & 0xff), - static_cast((header_len >> 8) & 0xff)}; + std::array header_len_le16{ + static_cast((header_len >> 0) & 0xff), + static_cast((header_len >> 8) & 0xff)}; out.write(reinterpret_cast(header_len_le16.data()), 2); } else { - auto header_len = static_cast(header_dict.length() + padding.length() + 1); + auto header_len = + static_cast(header_dict.length() + padding.length() + 1); std::array header_len_le32{ - static_cast((header_len >> 0) & 0xff), static_cast((header_len >> 8) & 0xff), - static_cast((header_len >> 16) & 0xff), static_cast((header_len >> 24) & 0xff)}; + static_cast((header_len >> 0) & 0xff), + static_cast((header_len >> 8) & 0xff), + static_cast((header_len >> 16) & 0xff), + static_cast((header_len >> 24) & 0xff)}; out.write(reinterpret_cast(header_len_le32.data()), 4); } diff --git a/src/NumpyHelpers.test.cpp b/src/NumpyHelpers.test.cpp index 36fcfe6..ad55b17 100644 --- a/src/NumpyHelpers.test.cpp +++ b/src/NumpyHelpers.test.cpp @@ -19,7 +19,9 @@ TEST_CASE("Check for quotes and return stripped string") { REQUIRE(parse_str("''") == ""); } -TEST_CASE("parsing a string without quotes throws") { REQUIRE_THROWS(parse_str("hej")); } +TEST_CASE("parsing a string without quotes throws") { + REQUIRE_THROWS(parse_str("hej")); +} TEST_CASE("trim whitespace") { REQUIRE(trim(" hej ") == "hej"); @@ -53,7 +55,8 @@ TEST_CASE("is element in array") { } TEST_CASE("Parse numpy dict") { - std::string in = "{'descr': ' keys{"descr", "fortran_order", "shape"}; auto map = parse_dict(in, keys); REQUIRE(map["descr"] == "' #include +#include #include #include @@ -58,7 +57,8 @@ TEST_CASE("test pedestal push") { if (k < 5) { REQUIRE(pedestal.cur_samples()(i, j) == k + 1); REQUIRE(pedestal.get_sum()(i, j) == (k + 1) * (i + j)); - REQUIRE(pedestal.get_sum2()(i, j) == (k + 1) * (i + j) * (i + j)); + REQUIRE(pedestal.get_sum2()(i, j) == + (k + 1) * (i + j) * (i + j)); } else { REQUIRE(pedestal.cur_samples()(i, j) == 5); REQUIRE(pedestal.get_sum()(i, j) == 5 * (i + j)); @@ -95,9 +95,12 @@ TEST_CASE("test pedestal with normal distribution") { for (int i = 0; i < 3; i++) { for (int j = 0; j < 5; j++) { - REQUIRE_THAT(mean(i, j), Catch::Matchers::WithinAbs(MEAN, MEAN * TOLERANCE)); - REQUIRE_THAT(variance(i, j), Catch::Matchers::WithinAbs(VAR, VAR * TOLERANCE)); - REQUIRE_THAT(standard_deviation(i, j), Catch::Matchers::WithinAbs(STD, STD * TOLERANCE)); + REQUIRE_THAT(mean(i, j), + Catch::Matchers::WithinAbs(MEAN, MEAN * TOLERANCE)); + REQUIRE_THAT(variance(i, j), + Catch::Matchers::WithinAbs(VAR, VAR * TOLERANCE)); + REQUIRE_THAT(standard_deviation(i, j), + Catch::Matchers::WithinAbs(STD, STD * TOLERANCE)); } } } \ No newline at end of file diff --git a/src/PixelMap.cpp b/src/PixelMap.cpp index d62759a..5be5ed4 100644 --- a/src/PixelMap.cpp +++ b/src/PixelMap.cpp @@ -31,7 +31,7 @@ NDArray GenerateMoench03PixelMap() { } NDArray GenerateMoench05PixelMap() { - std::array adc_numbers = {5, 9, 1}; + std::array adc_numbers = {5, 9, 1}; NDArray order_map({160, 150}); int n_pixel = 0; for (int row = 0; row < 160; row++) { @@ -40,11 +40,11 @@ NDArray GenerateMoench05PixelMap() { for (int i_sc = 0; i_sc < 3; i_sc++) { int col = 50 * i_sc + i_col; int adc_nr = adc_numbers[i_sc]; - int i_analog = n_pixel * 12 + adc_nr; + int i_analog = n_pixel * 12 + adc_nr; - // analog_frame[row * 150 + col] = analog_data[i_analog] & 0x3FFF; + // analog_frame[row * 150 + col] = analog_data[i_analog] & + // 0x3FFF; order_map(row, col) = i_analog; - } } } @@ -52,7 +52,7 @@ NDArray GenerateMoench05PixelMap() { } NDArray GenerateMoench05PixelMap1g() { - std::array adc_numbers = {1, 2, 0}; + std::array adc_numbers = {1, 2, 0}; NDArray order_map({160, 150}); int n_pixel = 0; for (int row = 0; row < 160; row++) { @@ -61,12 +61,11 @@ NDArray GenerateMoench05PixelMap1g() { for (int i_sc = 0; i_sc < 3; i_sc++) { int col = 50 * i_sc + i_col; int adc_nr = adc_numbers[i_sc]; - int i_analog = n_pixel * 3 + adc_nr; + int i_analog = n_pixel * 3 + adc_nr; - - // analog_frame[row * 150 + col] = analog_data[i_analog] & 0x3FFF; + // analog_frame[row * 150 + col] = analog_data[i_analog] & + // 0x3FFF; order_map(row, col) = i_analog; - } } } @@ -85,42 +84,42 @@ NDArray GenerateMoench05PixelMapOld() { int adc_nr = adc_numbers[i_sc]; int i_analog = n_pixel * 32 + adc_nr; - - // analog_frame[row * 150 + col] = analog_data[i_analog] & 0x3FFF; + // analog_frame[row * 150 + col] = analog_data[i_analog] & + // 0x3FFF; order_map(row, col) = i_analog; - } } } return order_map; } -NDArrayGenerateEigerFlipRowsPixelMap(){ +NDArray GenerateEigerFlipRowsPixelMap() { NDArray order_map({256, 512}); - for(int row = 0; row < 256; row++){ - for(int col = 0; col < 512; col++){ - order_map(row, col) = 255*512-row*512 + col; + for (int row = 0; row < 256; row++) { + for (int col = 0; col < 512; col++) { + order_map(row, col) = 255 * 512 - row * 512 + col; } } return order_map; } -NDArrayGenerateMH02SingleCounterPixelMap(){ +NDArray GenerateMH02SingleCounterPixelMap() { NDArray order_map({48, 48}); - for(int row = 0; row < 48; row++){ - for(int col = 0; col < 48; col++){ - order_map(row, col) = row*48 + col; + for (int row = 0; row < 48; row++) { + for (int col = 0; col < 48; col++) { + order_map(row, col) = row * 48 + col; } } return order_map; } -NDArray GenerateMH02FourCounterPixelMap(){ +NDArray GenerateMH02FourCounterPixelMap() { NDArray order_map({4, 48, 48}); - for (int counter=0; counter<4; counter++){ - for(int row = 0; row < 48; row++){ - for(int col = 0; col < 48; col++){ - order_map(counter, row, col) = counter*48*48 + row*48 + col; + for (int counter = 0; counter < 4; counter++) { + for (int row = 0; row < 48; row++) { + for (int col = 0; col < 48; col++) { + order_map(counter, row, col) = + counter * 48 * 48 + row * 48 + col; } } } diff --git a/src/RawFile.cpp b/src/RawFile.cpp index 122cf96..22bf8dd 100644 --- a/src/RawFile.cpp +++ b/src/RawFile.cpp @@ -1,9 +1,9 @@ #include "aare/RawFile.hpp" -#include "aare/algorithm.hpp" #include "aare/PixelMap.hpp" +#include "aare/algorithm.hpp" #include "aare/defs.hpp" -#include "aare/logger.hpp" #include "aare/geo_helpers.hpp" +#include "aare/logger.hpp" #include #include @@ -17,8 +17,9 @@ RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode) m_mode = mode; if (mode == "r") { find_geometry(); - if (m_master.roi()){ - m_geometry = update_geometry_with_roi(m_geometry, m_master.roi().value()); + if (m_master.roi()) { + m_geometry = + update_geometry_with_roi(m_geometry, m_master.roi().value()); } open_subfiles(); } else { @@ -47,32 +48,31 @@ void RawFile::read_into(std::byte *image_buf) { return get_frame_into(m_current_frame++, image_buf); } - void RawFile::read_into(std::byte *image_buf, DetectorHeader *header) { return get_frame_into(m_current_frame++, image_buf, header); } -void RawFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header) { +void RawFile::read_into(std::byte *image_buf, size_t n_frames, + DetectorHeader *header) { // return get_frame_into(m_current_frame++, image_buf, header); for (size_t i = 0; i < n_frames; i++) { this->get_frame_into(m_current_frame++, image_buf, header); image_buf += bytes_per_frame(); - if(header) - header+=n_modules(); + if (header) + header += n_modules(); } - } size_t RawFile::n_modules() const { return m_master.n_modules(); } - size_t RawFile::bytes_per_frame() { - return m_geometry.pixels_x * m_geometry.pixels_y * m_master.bitdepth() / bits_per_byte; + return m_geometry.pixels_x * m_geometry.pixels_y * m_master.bitdepth() / + bits_per_byte; } -size_t RawFile::pixels_per_frame() { - // return m_rows * m_cols; +size_t RawFile::pixels_per_frame() { + // return m_rows * m_cols; return m_geometry.pixels_x * m_geometry.pixels_y; } @@ -128,27 +128,25 @@ DetectorHeader RawFile::read_header(const std::filesystem::path &fname) { return h; } - RawMasterFile RawFile::master() const { return m_master; } /** * @brief Find the geometry of the detector by opening all the subfiles and - * reading the headers. + * reading the headers. */ void RawFile::find_geometry() { - - //Hold the maximal row and column number found - //Later used for calculating the total number of rows and columns + + // Hold the maximal row and column number found + // Later used for calculating the total number of rows and columns uint16_t r{}; uint16_t c{}; - for (size_t i = 0; i < n_modules(); i++) { auto h = read_header(m_master.data_fname(i, 0)); r = std::max(r, h.row); c = std::max(c, h.column); // positions.push_back({h.row, h.column}); - + ModuleGeometry g; g.origin_x = h.column * m_master.pixels_x(); g.origin_y = h.row * m_master.pixels_y(); @@ -157,34 +155,30 @@ void RawFile::find_geometry() { g.width = m_master.pixels_x(); g.height = m_master.pixels_y(); m_geometry.module_pixel_0.push_back(g); - } r++; c++; m_geometry.pixels_y = (r * m_master.pixels_y()); - m_geometry.pixels_x = (c * m_master.pixels_x()); + m_geometry.pixels_x = (c * m_master.pixels_x()); m_geometry.modules_x = c; m_geometry.modules_y = r; m_geometry.pixels_y += static_cast((r - 1) * cfg.module_gap_row); - } - Frame RawFile::get_frame(size_t frame_index) { - auto f = Frame(m_geometry.pixels_y, m_geometry.pixels_x, Dtype::from_bitdepth(m_master.bitdepth())); + auto f = Frame(m_geometry.pixels_y, m_geometry.pixels_x, + Dtype::from_bitdepth(m_master.bitdepth())); std::byte *frame_buffer = f.data(); get_frame_into(frame_index, frame_buffer); return f; } +size_t RawFile::bytes_per_pixel() const { return m_master.bitdepth() / 8; } -size_t RawFile::bytes_per_pixel() const { - return m_master.bitdepth() / 8; -} - -void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, DetectorHeader *header) { +void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, + DetectorHeader *header) { LOG(logDEBUG) << "RawFile::get_frame_into(" << frame_index << ")"; if (frame_index >= total_frames()) { throw std::runtime_error(LOCATION + "Frame number out of range"); @@ -192,12 +186,12 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect std::vector frame_numbers(n_modules()); std::vector frame_indices(n_modules(), frame_index); - // sync the frame numbers - if (n_modules() != 1) { //if we have more than one module + if (n_modules() != 1) { // if we have more than one module for (size_t part_idx = 0; part_idx != n_modules(); ++part_idx) { - frame_numbers[part_idx] = m_subfiles[part_idx]->frame_number(frame_index); + frame_numbers[part_idx] = + m_subfiles[part_idx]->frame_number(frame_index); } // 1. if frame number vector is the same break @@ -218,7 +212,8 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect } frame_numbers[min_frame_idx] = - m_subfiles[min_frame_idx]->frame_number(frame_indices[min_frame_idx]); + m_subfiles[min_frame_idx]->frame_number( + frame_indices[min_frame_idx]); } } @@ -226,15 +221,18 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect // get the part from each subfile and copy it to the frame for (size_t part_idx = 0; part_idx != n_modules(); ++part_idx) { auto corrected_idx = frame_indices[part_idx]; - - // This is where we start writing - auto offset = (m_geometry.module_pixel_0[part_idx].origin_y * m_geometry.pixels_x + - m_geometry.module_pixel_0[part_idx].origin_x)*m_master.bitdepth()/8; - if (m_geometry.module_pixel_0[part_idx].origin_x!=0) - throw std::runtime_error(LOCATION + " Implementation error. x pos not 0."); - - //TODO! What if the files don't match? + // This is where we start writing + auto offset = (m_geometry.module_pixel_0[part_idx].origin_y * + m_geometry.pixels_x + + m_geometry.module_pixel_0[part_idx].origin_x) * + m_master.bitdepth() / 8; + + if (m_geometry.module_pixel_0[part_idx].origin_x != 0) + throw std::runtime_error(LOCATION + + " Implementation error. x pos not 0."); + + // TODO! What if the files don't match? m_subfiles[part_idx]->seek(corrected_idx); m_subfiles[part_idx]->read_into(frame_buffer + offset, header); if (header) @@ -242,7 +240,7 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect } } else { - //TODO! should we read row by row? + // TODO! should we read row by row? // create a buffer large enough to hold a full module auto bytes_per_part = m_master.pixels_y() * m_master.pixels_x() * @@ -260,7 +258,7 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect m_subfiles[part_idx]->seek(corrected_idx); m_subfiles[part_idx]->read_into(part_buffer, header); - if(header) + if (header) ++header; for (size_t cur_row = 0; cur_row < static_cast(pos.height); @@ -271,15 +269,13 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect auto dest = (irow * this->m_geometry.pixels_x + icol); dest = dest * m_master.bitdepth() / 8; memcpy(frame_buffer + dest, - part_buffer + cur_row * pos.width * - m_master.bitdepth() / 8, + part_buffer + + cur_row * pos.width * m_master.bitdepth() / 8, pos.width * m_master.bitdepth() / 8); - } } delete[] part_buffer; } - } std::vector RawFile::read_n(size_t n_frames) { @@ -299,5 +295,4 @@ size_t RawFile::frame_number(size_t frame_index) { return m_subfiles[0]->frame_number(frame_index); } - } // namespace aare diff --git a/src/RawFile.test.cpp b/src/RawFile.test.cpp index 9109985..1fb441f 100644 --- a/src/RawFile.test.cpp +++ b/src/RawFile.test.cpp @@ -1,18 +1,18 @@ +#include "aare/RawFile.hpp" #include "aare/File.hpp" #include "aare/RawMasterFile.hpp" //needed for ROI -#include "aare/RawFile.hpp" #include #include #include "test_config.hpp" - using aare::File; TEST_CASE("Read number of frames from a jungfrau raw file", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); @@ -20,7 +20,8 @@ TEST_CASE("Read number of frames from a jungfrau raw file", "[.integration]") { } TEST_CASE("Read frame numbers from a jungfrau raw file", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); @@ -36,7 +37,8 @@ TEST_CASE("Read frame numbers from a jungfrau raw file", "[.integration]") { } TEST_CASE("Read a frame number too high throws", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); @@ -49,8 +51,10 @@ TEST_CASE("Read a frame number too high throws", "[.integration]") { REQUIRE_THROWS(f.frame_number(10)); } -TEST_CASE("Read a frame numbers where the subfile is missing throws", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_missing_subfile_master_0.json"; +TEST_CASE("Read a frame numbers where the subfile is missing throws", + "[.integration]") { + auto fpath = test_data_path() / "jungfrau" / + "jungfrau_missing_subfile_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); @@ -58,7 +62,7 @@ TEST_CASE("Read a frame numbers where the subfile is missing throws", "[.integra // we know this file has 10 frames with frame numbers 1 to 10 // f0 1,2,3 // f1 4,5,6 - but files f1-f3 are missing - // f2 7,8,9 - gone + // f2 7,8,9 - gone // f3 10 - gone REQUIRE(f.frame_number(0) == 1); REQUIRE(f.frame_number(1) == 2); @@ -69,15 +73,18 @@ TEST_CASE("Read a frame numbers where the subfile is missing throws", "[.integra REQUIRE_THROWS(f.frame_number(10)); } - -TEST_CASE("Read data from a jungfrau 500k single port raw file", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; +TEST_CASE("Read data from a jungfrau 500k single port raw file", + "[.integration]") { + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); - // we know this file has 10 frames with pixel 0,0 being: 2123, 2051, 2109, 2117, 2089, 2095, 2072, 2126, 2097, 2102 - std::vector pixel_0_0 = {2123, 2051, 2109, 2117, 2089, 2095, 2072, 2126, 2097, 2102}; + // we know this file has 10 frames with pixel 0,0 being: 2123, 2051, 2109, + // 2117, 2089, 2095, 2072, 2126, 2097, 2102 + std::vector pixel_0_0 = {2123, 2051, 2109, 2117, 2089, + 2095, 2072, 2126, 2097, 2102}; for (size_t i = 0; i < 10; i++) { auto frame = f.read_frame(); CHECK(frame.rows() == 512); @@ -100,10 +107,12 @@ TEST_CASE("Read frame numbers from a raw file", "[.integration]") { } TEST_CASE("Compare reading from a numpy file with a raw file", "[.files]") { - auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_master_0.json"; + auto fpath_raw = + test_data_path() / "raw/jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath_raw)); - auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy"; + auto fpath_npy = + test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy"; REQUIRE(std::filesystem::exists(fpath_npy)); File raw(fpath_raw, "r"); @@ -121,17 +130,23 @@ TEST_CASE("Compare reading from a numpy file with a raw file", "[.files]") { } TEST_CASE("Read multipart files", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_double_master_0.json"; + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_double_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); // we know this file has 10 frames check read_multiport.py for the values - std::vector pixel_0_0 = {2099, 2121, 2108, 2084, 2084, 2118, 2066, 2108, 2112, 2116}; - std::vector pixel_0_1 = {2842, 2796, 2865, 2798, 2805, 2817, 2852, 2789, 2792, 2833}; - std::vector pixel_255_1023 = {2149, 2037, 2115, 2102, 2118, 2090, 2036, 2071, 2073, 2142}; - std::vector pixel_511_1023 = {3231, 3169, 3167, 3162, 3168, 3160, 3171, 3171, 3169, 3171}; - std::vector pixel_1_0 = {2748, 2614, 2665, 2629, 2618, 2630, 2631, 2634, 2577, 2598}; + std::vector pixel_0_0 = {2099, 2121, 2108, 2084, 2084, + 2118, 2066, 2108, 2112, 2116}; + std::vector pixel_0_1 = {2842, 2796, 2865, 2798, 2805, + 2817, 2852, 2789, 2792, 2833}; + std::vector pixel_255_1023 = {2149, 2037, 2115, 2102, 2118, + 2090, 2036, 2071, 2073, 2142}; + std::vector pixel_511_1023 = {3231, 3169, 3167, 3162, 3168, + 3160, 3171, 3171, 3169, 3171}; + std::vector pixel_1_0 = {2748, 2614, 2665, 2629, 2618, + 2630, 2631, 2634, 2577, 2598}; for (size_t i = 0; i < 10; i++) { auto frame = f.read_frame(); @@ -146,11 +161,9 @@ TEST_CASE("Read multipart files", "[.integration]") { } TEST_CASE("Read file with unordered frames", "[.integration]") { - //TODO! Better explanation and error message + // TODO! Better explanation and error message auto fpath = test_data_path() / "mythen" / "scan242_master_3.raw"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath); REQUIRE_THROWS((f.read_frame())); } - - diff --git a/src/RawMasterFile.cpp b/src/RawMasterFile.cpp index 8a2db87..508b396 100644 --- a/src/RawMasterFile.cpp +++ b/src/RawMasterFile.cpp @@ -1,5 +1,5 @@ #include "aare/RawMasterFile.hpp" -#include +#include namespace aare { RawFileNameComponents::RawFileNameComponents( @@ -37,18 +37,15 @@ 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) const { + std::string fmt = "{}_d{}_f{}_{}.raw"; - //Before version X we used to name the data files f000000000000 + // 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); + 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) { @@ -65,19 +62,19 @@ const std::string &RawFileNameComponents::ext() const { return m_ext; } int RawFileNameComponents::file_index() const { return m_file_index; } // "[enabled\ndac dac 4\nstart 500\nstop 2200\nstep 5\nsettleTime 100us\n]" -ScanParameters::ScanParameters(const std::string& par){ - std::istringstream iss(par.substr(1, par.size()-2)); +ScanParameters::ScanParameters(const std::string &par) { + std::istringstream iss(par.substr(1, par.size() - 2)); std::string line; - while(std::getline(iss, line)){ - if(line == "enabled"){ + while (std::getline(iss, line)) { + if (line == "enabled") { m_enabled = true; - }else if(line.find("dac") != std::string::npos){ + } else if (line.find("dac") != std::string::npos) { m_dac = line.substr(4); - }else if(line.find("start") != std::string::npos){ + } else if (line.find("start") != std::string::npos) { m_start = std::stoi(line.substr(6)); - }else if(line.find("stop") != std::string::npos){ + } else if (line.find("stop") != std::string::npos) { m_stop = std::stoi(line.substr(5)); - }else if(line.find("step") != std::string::npos){ + } else if (line.find("step") != std::string::npos) { m_step = std::stoi(line.substr(5)); } } @@ -85,14 +82,11 @@ ScanParameters::ScanParameters(const std::string& par){ int ScanParameters::start() const { return m_start; } int ScanParameters::stop() const { return m_stop; } -void ScanParameters::increment_stop(){ - m_stop += 1; -} +void ScanParameters::increment_stop() { m_stop += 1; } int ScanParameters::step() const { return m_step; } const std::string &ScanParameters::dac() const { return m_dac; } bool ScanParameters::enabled() const { return m_enabled; } - RawMasterFile::RawMasterFile(const std::filesystem::path &fpath) : m_fnc(fpath) { if (!std::filesystem::exists(fpath)) { @@ -163,10 +157,8 @@ ScanParameters RawMasterFile::scan_parameters() const { return m_scan_parameters; } - std::optional RawMasterFile::roi() const { return m_roi; } - void RawMasterFile::parse_json(const std::filesystem::path &fpath) { std::ifstream ifs(fpath); json j; @@ -205,17 +197,16 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) { // keep the optional empty } - // ---------------------------------------------------------------- // Special treatment of analog flag because of Moench03 - try{ + try { m_analog_flag = j.at("Analog Flag"); - }catch (const json::out_of_range &e) { + } catch (const json::out_of_range &e) { // if it doesn't work still set it to one // to try to decode analog samples (Old Moench03) m_analog_flag = 1; } - try { + try { if (m_analog_flag) { m_analog_samples = j.at("Analog Samples"); } @@ -248,27 +239,27 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) { // keep the optional empty } - try{ + try { m_transceiver_flag = j.at("Transceiver Flag"); - if(m_transceiver_flag){ + if (m_transceiver_flag) { m_transceiver_samples = j.at("Transceiver Samples"); } - }catch (const json::out_of_range &e) { + } catch (const json::out_of_range &e) { // keep the optional empty } - try{ + try { std::string scan_parameters = j.at("Scan Parameters"); m_scan_parameters = ScanParameters(scan_parameters); - if(v<7.21){ - m_scan_parameters.increment_stop(); //adjust for endpoint being included - } - }catch (const json::out_of_range &e) { + if (v < 7.21) { + m_scan_parameters + .increment_stop(); // adjust for endpoint being included + } + } catch (const json::out_of_range &e) { // not a scan } - - try{ + try { ROI tmp_roi; auto obj = j.at("Receiver Roi"); tmp_roi.xmin = obj.at("xmin"); @@ -276,37 +267,32 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) { tmp_roi.ymin = obj.at("ymin"); tmp_roi.ymax = obj.at("ymax"); - //if any of the values are set update the roi + // if any of the values are set update the roi if (tmp_roi.xmin != 4294967295 || tmp_roi.xmax != 4294967295 || tmp_roi.ymin != 4294967295 || tmp_roi.ymax != 4294967295) { - - if(v<7.21){ + + if (v < 7.21) { tmp_roi.xmax++; tmp_roi.ymax++; } - + m_roi = tmp_roi; } - - }catch (const json::out_of_range &e) { + } catch (const json::out_of_range &e) { // leave the optional empty } - //if we have an roi we need to update the geometry for the subfiles - if (m_roi){ - + // if we have an roi we need to update the geometry for the subfiles + if (m_roi) { } - - - - // Update detector type for Moench - // TODO! How does this work with old .raw master files? - #ifdef AARE_VERBOSE +// Update detector type for Moench +// TODO! How does this work with old .raw master files? +#ifdef AARE_VERBOSE fmt::print("Detecting Moench03: m_pixels_y: {}, m_analog_samples: {}\n", m_pixels_y, m_analog_samples.value_or(0)); - #endif +#endif if (m_type == DetectorType::Moench && !m_analog_samples && m_pixels_y == 400) { m_type = DetectorType::Moench03; @@ -332,19 +318,19 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) { if (key == "Version") { m_version = value; - //TODO!: How old versions can we handle? + // 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) + // 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(value); - if(m_type==DetectorType::Moench){ + if (m_type == DetectorType::Moench) { m_type = DetectorType::Moench03_old; } } else if (key == "Timing Mode") { @@ -381,10 +367,10 @@ 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"){ + } else if (key == "row") { pos = value.find('p'); m_pixels_y = std::stoi(value.substr(0, pos)); - }else if(key == "col"){ + } else if (key == "col") { pos = value.find('p'); m_pixels_x = std::stoi(value.substr(0, pos)); } else if (key == "Total Frames") { @@ -395,8 +381,8 @@ 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 + } 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(','); @@ -410,15 +396,14 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) { 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}; + // 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) + // TODO! Read files and find actual frames + if (m_frames_in_file == 0) m_frames_in_file = m_total_frames_expected; } } // namespace aare diff --git a/src/RawMasterFile.test.cpp b/src/RawMasterFile.test.cpp index 560217f..0e75ec4 100644 --- a/src/RawMasterFile.test.cpp +++ b/src/RawMasterFile.test.cpp @@ -1,12 +1,11 @@ #include "aare/RawMasterFile.hpp" -#include #include "test_config.hpp" +#include using namespace aare; - -TEST_CASE("Parse a master file fname"){ +TEST_CASE("Parse a master file fname") { RawFileNameComponents m("test_master_1.json"); REQUIRE(m.base_name() == "test"); REQUIRE(m.ext() == ".json"); @@ -14,7 +13,7 @@ TEST_CASE("Parse a master file fname"){ REQUIRE(m.base_path() == ""); } -TEST_CASE("Extraction of base path works"){ +TEST_CASE("Extraction of base path works") { RawFileNameComponents m("some/path/test_master_73.json"); REQUIRE(m.base_name() == "test"); REQUIRE(m.ext() == ".json"); @@ -22,7 +21,7 @@ TEST_CASE("Extraction of base path works"){ REQUIRE(m.base_path() == "some/path"); } -TEST_CASE("Construction of master file name and data files"){ +TEST_CASE("Construction of master file name and data files") { RawFileNameComponents m("test_master_1.json"); REQUIRE(m.master_fname() == "test_master_1.json"); REQUIRE(m.data_fname(0, 0) == "test_d0_f0_1.raw"); @@ -31,7 +30,7 @@ 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"){ +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"); @@ -41,16 +40,15 @@ TEST_CASE("Construction of master file name and data files using old scheme"){ REQUIRE(m.data_fname(1, 1) == "test_d1_f000000000001_1.raw"); } -TEST_CASE("Master file name does not fit pattern"){ +TEST_CASE("Master file name does not fit pattern") { REQUIRE_THROWS(RawFileNameComponents("somefile.json")); REQUIRE_THROWS(RawFileNameComponents("another_test_d0_f0_1.raw")); REQUIRE_THROWS(RawFileNameComponents("test_master_1.txt")); } - - -TEST_CASE("Parse scan parameters"){ - ScanParameters s("[enabled\ndac dac 4\nstart 500\nstop 2200\nstep 5\nsettleTime 100us\n]"); +TEST_CASE("Parse scan parameters") { + ScanParameters s("[enabled\ndac dac 4\nstart 500\nstop 2200\nstep " + "5\nsettleTime 100us\n]"); REQUIRE(s.enabled()); REQUIRE(s.dac() == "dac 4"); REQUIRE(s.start() == 500); @@ -58,7 +56,7 @@ TEST_CASE("Parse scan parameters"){ REQUIRE(s.step() == 5); } -TEST_CASE("A disabled scan"){ +TEST_CASE("A disabled scan") { ScanParameters s("[disabled]"); REQUIRE_FALSE(s.enabled()); REQUIRE(s.dac() == ""); @@ -67,9 +65,9 @@ TEST_CASE("A disabled scan"){ REQUIRE(s.step() == 0); } - -TEST_CASE("Parse a master file in .json format", "[.integration]"){ - auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; +TEST_CASE("Parse a master file in .json format", "[.integration]") { + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); RawMasterFile f(fpath); @@ -80,7 +78,7 @@ TEST_CASE("Parse a master file in .json format", "[.integration]"){ REQUIRE(f.detector_type() == DetectorType::Jungfrau); // "Timing Mode": "auto", REQUIRE(f.timing_mode() == TimingMode::Auto); - + // "Geometry": { // "x": 1, // "y": 1 @@ -100,10 +98,9 @@ TEST_CASE("Parse a master file in .json format", "[.integration]"){ // "Max Frames Per File": 3, REQUIRE(f.max_frames_per_file() == 3); - //Jungfrau doesn't write but it is 16 + // Jungfrau doesn't write but it is 16 REQUIRE(f.bitdepth() == 16); - // "Frame Discard Policy": "nodiscard", // "Frame Padding": 1, @@ -125,33 +122,35 @@ TEST_CASE("Parse a master file in .json format", "[.integration]"){ // "Frames in File": 10, REQUIRE(f.frames_in_file() == 10); - //TODO! Should we parse this? - // "Frame Header Format": { - // "Frame Number": "8 bytes", - // "SubFrame Number/ExpLength": "4 bytes", - // "Packet Number": "4 bytes", - // "Bunch ID": "8 bytes", - // "Timestamp": "8 bytes", - // "Module Id": "2 bytes", - // "Row": "2 bytes", - // "Column": "2 bytes", - // "Reserved": "2 bytes", - // "Debug": "4 bytes", - // "Round Robin Number": "2 bytes", - // "Detector Type": "1 byte", - // "Header Version": "1 byte", - // "Packets Caught Mask": "64 bytes" - // } - // } + // TODO! Should we parse this? + // "Frame Header Format": { + // "Frame Number": "8 bytes", + // "SubFrame Number/ExpLength": "4 bytes", + // "Packet Number": "4 bytes", + // "Bunch ID": "8 bytes", + // "Timestamp": "8 bytes", + // "Module Id": "2 bytes", + // "Row": "2 bytes", + // "Column": "2 bytes", + // "Reserved": "2 bytes", + // "Debug": "4 bytes", + // "Round Robin Number": "2 bytes", + // "Detector Type": "1 byte", + // "Header Version": "1 byte", + // "Packets Caught Mask": "64 bytes" + // } + // } REQUIRE_FALSE(f.analog_samples()); - REQUIRE_FALSE(f.digital_samples()); - + REQUIRE_FALSE(f.digital_samples()); } -TEST_CASE("Parse a master file in .raw format", "[.integration]"){ - - auto fpath = test_data_path() / "moench/moench04_noise_200V_sto_both_100us_no_light_thresh_900_master_0.raw"; +TEST_CASE("Parse a master file in .raw format", "[.integration]") { + + auto fpath = + test_data_path() / + "moench/" + "moench04_noise_200V_sto_both_100us_no_light_thresh_900_master_0.raw"; REQUIRE(std::filesystem::exists(fpath)); RawMasterFile f(fpath); @@ -209,80 +208,74 @@ TEST_CASE("Parse a master file in .raw format", "[.integration]"){ // Detector Type : 1 byte // Header Version : 1 byte // Packets Caught Mask : 64 bytes - - } - -TEST_CASE("Read eiger master file", "[.integration]"){ -auto fpath = test_data_path() / "eiger" / "eiger_500k_32bit_master_0.json"; +TEST_CASE("Read eiger master file", "[.integration]") { + auto fpath = test_data_path() / "eiger" / "eiger_500k_32bit_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); RawMasterFile f(fpath); - -// { -// "Version": 7.2, -REQUIRE(f.version() == "7.2"); -// "Timestamp": "Tue Mar 26 17:24:34 2024", -// "Detector Type": "Eiger", -REQUIRE(f.detector_type() == DetectorType::Eiger); -// "Timing Mode": "auto", -REQUIRE(f.timing_mode() == TimingMode::Auto); -// "Geometry": { -// "x": 2, -// "y": 2 -// }, -// "Image Size in bytes": 524288, -REQUIRE(f.image_size_in_bytes() == 524288); -// "Pixels": { -// "x": 512, -REQUIRE(f.pixels_x() == 512); -// "y": 256 -REQUIRE(f.pixels_y() == 256); -// }, -// "Max Frames Per File": 10000, -REQUIRE(f.max_frames_per_file() == 10000); -// "Frame Discard Policy": "nodiscard", -REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard); -// "Frame Padding": 1, -REQUIRE(f.frame_padding() == 1); - -// "Scan Parameters": "[disabled]", -// "Total Frames": 3, -// "Receiver Roi": { -// "xmin": 4294967295, -// "xmax": 4294967295, -// "ymin": 4294967295, -// "ymax": 4294967295 -// }, -// "Dynamic Range": 32, -// "Ten Giga": 0, -// "Exptime": "5s", -// "Period": "1s", -// "Threshold Energy": -1, -// "Sub Exptime": "2.62144ms", -// "Sub Period": "2.62144ms", -// "Quad": 0, -// "Number of rows": 256, -// "Rate Corrections": "[0, 0]", -// "Frames in File": 3, -// "Frame Header Format": { -// "Frame Number": "8 bytes", -// "SubFrame Number/ExpLength": "4 bytes", -// "Packet Number": "4 bytes", -// "Bunch ID": "8 bytes", -// "Timestamp": "8 bytes", -// "Module Id": "2 bytes", -// "Row": "2 bytes", -// "Column": "2 bytes", -// "Reserved": "2 bytes", -// "Debug": "4 bytes", -// "Round Robin Number": "2 bytes", -// "Detector Type": "1 byte", -// "Header Version": "1 byte", -// "Packets Caught Mask": "64 bytes" -// } -// } - + // { + // "Version": 7.2, + REQUIRE(f.version() == "7.2"); + // "Timestamp": "Tue Mar 26 17:24:34 2024", + // "Detector Type": "Eiger", + REQUIRE(f.detector_type() == DetectorType::Eiger); + // "Timing Mode": "auto", + REQUIRE(f.timing_mode() == TimingMode::Auto); + // "Geometry": { + // "x": 2, + // "y": 2 + // }, + // "Image Size in bytes": 524288, + REQUIRE(f.image_size_in_bytes() == 524288); + // "Pixels": { + // "x": 512, + REQUIRE(f.pixels_x() == 512); + // "y": 256 + REQUIRE(f.pixels_y() == 256); + // }, + // "Max Frames Per File": 10000, + REQUIRE(f.max_frames_per_file() == 10000); + // "Frame Discard Policy": "nodiscard", + REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard); + // "Frame Padding": 1, + REQUIRE(f.frame_padding() == 1); + // "Scan Parameters": "[disabled]", + // "Total Frames": 3, + // "Receiver Roi": { + // "xmin": 4294967295, + // "xmax": 4294967295, + // "ymin": 4294967295, + // "ymax": 4294967295 + // }, + // "Dynamic Range": 32, + // "Ten Giga": 0, + // "Exptime": "5s", + // "Period": "1s", + // "Threshold Energy": -1, + // "Sub Exptime": "2.62144ms", + // "Sub Period": "2.62144ms", + // "Quad": 0, + // "Number of rows": 256, + // "Rate Corrections": "[0, 0]", + // "Frames in File": 3, + // "Frame Header Format": { + // "Frame Number": "8 bytes", + // "SubFrame Number/ExpLength": "4 bytes", + // "Packet Number": "4 bytes", + // "Bunch ID": "8 bytes", + // "Timestamp": "8 bytes", + // "Module Id": "2 bytes", + // "Row": "2 bytes", + // "Column": "2 bytes", + // "Reserved": "2 bytes", + // "Debug": "4 bytes", + // "Round Robin Number": "2 bytes", + // "Detector Type": "1 byte", + // "Header Version": "1 byte", + // "Packets Caught Mask": "64 bytes" + // } + // } } \ No newline at end of file diff --git a/src/RawSubFile.cpp b/src/RawSubFile.cpp index a8d29ce..3ed2c6f 100644 --- a/src/RawSubFile.cpp +++ b/src/RawSubFile.cpp @@ -1,36 +1,30 @@ #include "aare/RawSubFile.hpp" #include "aare/PixelMap.hpp" #include "aare/algorithm.hpp" -#include "aare/utils/ifstream_helpers.hpp" #include "aare/logger.hpp" - +#include "aare/utils/ifstream_helpers.hpp" #include // memcpy #include #include #include - - - namespace aare { RawSubFile::RawSubFile(const std::filesystem::path &fname, DetectorType detector, size_t rows, size_t cols, size_t bitdepth, uint32_t pos_row, uint32_t pos_col) - : m_detector_type(detector), m_bitdepth(bitdepth), - m_rows(rows), m_cols(cols), - m_bytes_per_frame((m_bitdepth / 8) * m_rows * m_cols), m_pos_row(pos_row), - m_pos_col(pos_col) { + : m_detector_type(detector), m_bitdepth(bitdepth), m_rows(rows), + m_cols(cols), m_bytes_per_frame((m_bitdepth / 8) * m_rows * m_cols), + m_pos_row(pos_row), m_pos_col(pos_col) { - LOG(logDEBUG) << "RawSubFile::RawSubFile()"; + LOG(logDEBUG) << "RawSubFile::RawSubFile()"; if (m_detector_type == DetectorType::Moench03_old) { m_pixel_map = GenerateMoench03PixelMap(); } else if (m_detector_type == DetectorType::Eiger && m_pos_row % 2 == 0) { m_pixel_map = GenerateEigerFlipRowsPixelMap(); } - parse_fname(fname); scan_files(); open_file(m_current_file_index); // open the first file @@ -51,7 +45,8 @@ void RawSubFile::seek(size_t 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 + sizeof(DetectorHeader)); + auto byte_offset = + frame_offset * (m_bytes_per_frame + sizeof(DetectorHeader)); m_file.seekg(byte_offset); } @@ -69,7 +64,7 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) { m_file.seekg(sizeof(DetectorHeader), std::ios::cur); } - if (m_file.fail()){ + if (m_file.fail()) { throw std::runtime_error(LOCATION + ifstream_error_msg(m_file)); } @@ -78,14 +73,15 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) { // read into a temporary buffer and then copy the data to the buffer // in the correct order // TODO! add 4 bit support - if(m_bitdepth == 8){ + if (m_bitdepth == 8) { read_with_map(image_buf); - }else if (m_bitdepth == 16) { + } else if (m_bitdepth == 16) { read_with_map(image_buf); } else if (m_bitdepth == 32) { read_with_map(image_buf); - }else{ - throw std::runtime_error("Unsupported bitdepth for read with pixel map"); + } else { + throw std::runtime_error( + "Unsupported bitdepth for read with pixel map"); } } else { @@ -93,11 +89,11 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) { m_file.read(reinterpret_cast(image_buf), bytes_per_frame()); } - if (m_file.fail()){ + if (m_file.fail()) { throw std::runtime_error(LOCATION + ifstream_error_msg(m_file)); } - ++ m_current_frame_index; + ++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; @@ -105,7 +101,8 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) { } } -void RawSubFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header) { +void RawSubFile::read_into(std::byte *image_buf, size_t n_frames, + DetectorHeader *header) { for (size_t i = 0; i < n_frames; i++) { read_into(image_buf, header); image_buf += bytes_per_frame(); @@ -115,10 +112,7 @@ void RawSubFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader } } - - -template -void RawSubFile::read_with_map(std::byte *image_buf) { +template void RawSubFile::read_with_map(std::byte *image_buf) { auto part_buffer = new std::byte[bytes_per_frame()]; m_file.read(reinterpret_cast(part_buffer), bytes_per_frame()); auto *data = reinterpret_cast(image_buf); @@ -157,14 +151,17 @@ void RawSubFile::parse_fname(const std::filesystem::path &fname) { std::smatch match; if (std::regex_match(m_base_name, match, pattern)) { - m_offset = std::stoi(match[4].str()); // find the first file index in case of a truncated series - m_base_name = match[1].str() + match[2].str() + match[3].str() + "{}" + match[5].str(); + m_offset = std::stoi(match[4].str()); // find the first file index in + // case of a truncated series + m_base_name = match[1].str() + match[2].str() + match[3].str() + "{}" + + match[5].str(); LOG(logDEBUG) << "Base name: " << m_base_name; LOG(logDEBUG) << "Offset: " << m_offset; LOG(logDEBUG) << "Path: " << m_path.string(); } else { throw std::runtime_error( - LOCATION + fmt::format("Could not parse file name {}", fname.string())); + LOCATION + + fmt::format("Could not parse file name {}", fname.string())); } } @@ -175,12 +172,13 @@ std::filesystem::path RawSubFile::fpath(size_t file_index) const { void RawSubFile::open_file(size_t file_index) { m_file.close(); - auto fname = fpath(file_index+m_offset); + auto fname = fpath(file_index + m_offset); LOG(logDEBUG) << "RawSubFile::open_file(): " << fname.string(); m_file.open(fname, std::ios::binary); if (!m_file.is_open()) { throw std::runtime_error( - LOCATION + fmt::format("Could not open file {}", fpath(file_index).string())); + LOCATION + + fmt::format("Could not open file {}", fpath(file_index).string())); } m_current_file_index = file_index; } @@ -190,20 +188,21 @@ void RawSubFile::scan_files() { // find how many files we have and the number of frames in each file 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 + sizeof(DetectorHeader)); m_last_frame_in_file.push_back(n_frames); - LOG(logDEBUG) << "Found: " << n_frames << " frames in file: " << fpath(file_index).string(); + LOG(logDEBUG) << "Found: " << n_frames + << " frames in file: " << fpath(file_index).string(); ++file_index; } // find where we need to open the next file and total number of frames m_last_frame_in_file = cumsum(m_last_frame_in_file); - if(m_last_frame_in_file.empty()){ + if (m_last_frame_in_file.empty()) { m_total_frames = 0; - }else{ + } else { m_total_frames = m_last_frame_in_file.back(); } } diff --git a/src/RawSubFile.test.cpp b/src/RawSubFile.test.cpp index 89cf858..37d071b 100644 --- a/src/RawSubFile.test.cpp +++ b/src/RawSubFile.test.cpp @@ -1,13 +1,14 @@ #include "aare/RawSubFile.hpp" #include "aare/File.hpp" #include "aare/NDArray.hpp" -#include #include "test_config.hpp" +#include using namespace aare; -TEST_CASE("Read frames directly from a RawSubFile", "[.files]"){ - auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f0_0.raw"; +TEST_CASE("Read frames directly from a RawSubFile", "[.files]") { + auto fpath_raw = + test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f0_0.raw"; REQUIRE(std::filesystem::exists(fpath_raw)); RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16); @@ -17,19 +18,19 @@ TEST_CASE("Read frames directly from a RawSubFile", "[.files]"){ REQUIRE(f.bytes_per_frame() == 512 * 1024 * 2); REQUIRE(f.bytes_per_pixel() == 2); - - auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy"; + auto fpath_npy = + test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy"; REQUIRE(std::filesystem::exists(fpath_npy)); - //Numpy file with the same data to use as reference + // Numpy file with the same data to use as reference File npy(fpath_npy, "r"); CHECK(f.frames_in_file() == 10); CHECK(npy.total_frames() == 10); - DetectorHeader header{}; - NDArray image({static_cast(f.rows()), static_cast(f.cols())}); + NDArray image( + {static_cast(f.rows()), static_cast(f.cols())}); for (size_t i = 0; i < 10; ++i) { CHECK(f.tell() == i); f.read_into(image.buffer(), &header); @@ -38,38 +39,40 @@ TEST_CASE("Read frames directly from a RawSubFile", "[.files]"){ } } -TEST_CASE("Read frames directly from a RawSubFile starting at the second file", "[.files]"){ +TEST_CASE("Read frames directly from a RawSubFile starting at the second file", + "[.files]") { // we know this file has 10 frames with frame numbers 1 to 10 // f0 1,2,3 // f1 4,5,6 <-- starting here // f2 7,8,9 // f3 10 - - auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f1_0.raw"; + + auto fpath_raw = + test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f1_0.raw"; REQUIRE(std::filesystem::exists(fpath_raw)); RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16); - - auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy"; + auto fpath_npy = + test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy"; REQUIRE(std::filesystem::exists(fpath_npy)); - //Numpy file with the same data to use as reference + // Numpy file with the same data to use as reference File npy(fpath_npy, "r"); npy.seek(3); CHECK(f.frames_in_file() == 7); CHECK(npy.total_frames() == 10); - DetectorHeader header{}; - NDArray image({static_cast(f.rows()), static_cast(f.cols())}); + NDArray image( + {static_cast(f.rows()), static_cast(f.cols())}); for (size_t i = 0; i < 7; ++i) { CHECK(f.tell() == i); f.read_into(image.buffer(), &header); - // frame numbers start at 1 frame index at 0 + // frame numbers start at 1 frame index at 0 // adding 3 + 1 to verify the frame number - CHECK(header.frameNumber == i + 4); + CHECK(header.frameNumber == i + 4); auto npy_frame = npy.read_frame(); CHECK((image.view() == npy_frame.view())); } diff --git a/src/algorithm.test.cpp b/src/algorithm.test.cpp index bf49c52..3706aba 100644 --- a/src/algorithm.test.cpp +++ b/src/algorithm.test.cpp @@ -161,12 +161,10 @@ TEST_CASE("cumsum works with negative numbers", "[algorithm]") { REQUIRE(result[4] == -10); } - TEST_CASE("cumsum on an empty vector", "[algorithm]") { std::vector vec = {}; auto result = aare::cumsum(vec); REQUIRE(result.size() == 0); - } TEST_CASE("All equal on an empty vector is false", "[algorithm]") { @@ -184,7 +182,8 @@ TEST_CASE("All equal on a vector with 2 elements is true", "[algorithm]") { REQUIRE(aare::all_equal(vec) == true); } -TEST_CASE("All equal on a vector with two different elements is false", "[algorithm]") { +TEST_CASE("All equal on a vector with two different elements is false", + "[algorithm]") { std::vector vec = {1, 2}; REQUIRE(aare::all_equal(vec) == false); } diff --git a/src/decode.cpp b/src/decode.cpp index 436ad7b..0c78b70 100644 --- a/src/decode.cpp +++ b/src/decode.cpp @@ -2,9 +2,9 @@ #include namespace aare { -uint16_t adc_sar_05_decode64to16(uint64_t input){ +uint16_t adc_sar_05_decode64to16(uint64_t input) { - //we want bits 29,19,28,18,31,21,27,20,24,23,25,22 and then pad to 16 + // we want bits 29,19,28,18,31,21,27,20,24,23,25,22 and then pad to 16 uint16_t output = 0; output |= ((input >> 22) & 1) << 11; output |= ((input >> 25) & 1) << 10; @@ -21,19 +21,21 @@ uint16_t adc_sar_05_decode64to16(uint64_t input){ return output; } -void adc_sar_05_decode64to16(NDView input, NDView output){ - if(input.shape() != output.shape()){ - throw std::invalid_argument(LOCATION + " input and output shapes must match"); +void adc_sar_05_decode64to16(NDView input, + NDView output) { + if (input.shape() != output.shape()) { + throw std::invalid_argument(LOCATION + + " input and output shapes must match"); } - for(ssize_t i = 0; i < input.shape(0); i++){ - for(ssize_t j = 0; j < input.shape(1); j++){ - output(i,j) = adc_sar_05_decode64to16(input(i,j)); + for (ssize_t i = 0; i < input.shape(0); i++) { + for (ssize_t j = 0; j < input.shape(1); j++) { + output(i, j) = adc_sar_05_decode64to16(input(i, j)); } } } -uint16_t adc_sar_04_decode64to16(uint64_t input){ +uint16_t adc_sar_04_decode64to16(uint64_t input) { // bit_map = array([15,17,19,21,23,4,6,8,10,12,14,16] LSB->MSB uint16_t output = 0; @@ -52,20 +54,23 @@ uint16_t adc_sar_04_decode64to16(uint64_t input){ return output; } -void adc_sar_04_decode64to16(NDView input, NDView output){ - if(input.shape() != output.shape()){ - throw std::invalid_argument(LOCATION + " input and output shapes must match"); +void adc_sar_04_decode64to16(NDView input, + NDView output) { + if (input.shape() != output.shape()) { + throw std::invalid_argument(LOCATION + + " input and output shapes must match"); } - for(ssize_t i = 0; i < input.shape(0); i++){ - for(ssize_t j = 0; j < input.shape(1); j++){ - output(i,j) = adc_sar_04_decode64to16(input(i,j)); + for (ssize_t i = 0; i < input.shape(0); i++) { + for (ssize_t j = 0; j < input.shape(1); j++) { + output(i, j) = adc_sar_04_decode64to16(input(i, j)); } } } double apply_custom_weights(uint16_t input, const NDView weights) { - if(weights.size() > 16){ - throw std::invalid_argument("weights size must be less than or equal to 16"); + if (weights.size() > 16) { + throw std::invalid_argument( + "weights size must be less than or equal to 16"); } double result = 0.0; @@ -73,30 +78,30 @@ double apply_custom_weights(uint16_t input, const NDView weights) { result += ((input >> i) & 1) * std::pow(weights[i], i); } return result; - } -void apply_custom_weights(NDView input, NDView output, const NDView weights) { - if(input.shape() != output.shape()){ - throw std::invalid_argument(LOCATION + " input and output shapes must match"); +void apply_custom_weights(NDView input, NDView output, + const NDView weights) { + if (input.shape() != output.shape()) { + throw std::invalid_argument(LOCATION + + " input and output shapes must match"); } - //Calculate weights to avoid repeatedly calling std::pow + // Calculate weights to avoid repeatedly calling std::pow std::vector weights_powers(weights.size()); for (ssize_t i = 0; i < weights.size(); ++i) { weights_powers[i] = std::pow(weights[i], i); } // Apply custom weights to each element in the input array - for (ssize_t i = 0; i < input.shape(0); i++) { + for (ssize_t i = 0; i < input.shape(0); i++) { double result = 0.0; - for (size_t bit_index = 0; bit_index < weights_powers.size(); ++bit_index) { + for (size_t bit_index = 0; bit_index < weights_powers.size(); + ++bit_index) { result += ((input(i) >> bit_index) & 1) * weights_powers[bit_index]; } output(i) = result; } } - - } // namespace aare diff --git a/src/decode.test.cpp b/src/decode.test.cpp index 1e4b2fc..b4310ca 100644 --- a/src/decode.test.cpp +++ b/src/decode.test.cpp @@ -1,17 +1,16 @@ #include "aare/decode.hpp" -#include -#include #include "aare/NDArray.hpp" +#include +#include using Catch::Matchers::WithinAbs; #include -TEST_CASE("test_adc_sar_05_decode64to16"){ +TEST_CASE("test_adc_sar_05_decode64to16") { uint64_t input = 0; uint16_t output = aare::adc_sar_05_decode64to16(input); CHECK(output == 0); - // bit 29 on th input is bit 0 on the output input = 1UL << 29; output = aare::adc_sar_05_decode64to16(input); @@ -25,7 +24,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ CHECK(output == (1 << i)); } - // test a few "random" values input = 0; input |= (1UL << 29); @@ -34,7 +32,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ output = aare::adc_sar_05_decode64to16(input); CHECK(output == 7UL); - input = 0; input |= (1UL << 18); input |= (1UL << 27); @@ -47,10 +44,9 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ input |= (1UL << 22); output = aare::adc_sar_05_decode64to16(input); CHECK(output == 3072UL); - } +} - - TEST_CASE("test_apply_custom_weights") { +TEST_CASE("test_apply_custom_weights") { uint16_t input = 1; aare::NDArray weights_data({3}, 0.0); @@ -60,7 +56,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ auto weights = weights_data.view(); - double output = aare::apply_custom_weights(input, weights); CHECK_THAT(output, WithinAbs(1.0, 0.001)); @@ -68,7 +63,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ output = aare::apply_custom_weights(input, weights); CHECK_THAT(output, WithinAbs(2.1, 0.001)); - input = 1 << 2; output = aare::apply_custom_weights(input, weights); CHECK_THAT(output, WithinAbs(3.24, 0.001)); @@ -76,5 +70,4 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ input = 0b111; output = aare::apply_custom_weights(input, weights); CHECK_THAT(output, WithinAbs(6.34, 0.001)); - - } \ No newline at end of file +} \ No newline at end of file diff --git a/src/defs.cpp b/src/defs.cpp index 7c7cc45..b93ff06 100644 --- a/src/defs.cpp +++ b/src/defs.cpp @@ -5,15 +5,11 @@ #include namespace aare { - -void assert_failed(const std::string &msg) - { +void assert_failed(const std::string &msg) { fmt::print(msg); exit(1); } - - /** * @brief Convert a DetectorType to a string * @param type DetectorType @@ -40,7 +36,7 @@ template <> std::string ToString(DetectorType arg) { case DetectorType::Xilinx_ChipTestBoard: return "Xilinx_ChipTestBoard"; - //Custom ones + // Custom ones case DetectorType::Moench03: return "Moench03"; case DetectorType::Moench03_old: @@ -48,8 +44,8 @@ template <> std::string ToString(DetectorType arg) { case DetectorType::Unknown: return "Unknown"; - //no default case to trigger compiler warning if not all - //enum values are handled + // no default case to trigger compiler warning if not all + // enum values are handled } throw std::runtime_error("Could not decode detector to string"); } @@ -80,14 +76,14 @@ template <> DetectorType StringTo(const std::string &arg) { if (arg == "Xilinx_ChipTestBoard") return DetectorType::Xilinx_ChipTestBoard; - //Custom ones + // Custom ones if (arg == "Moench03") return DetectorType::Moench03; if (arg == "Moench03_old") return DetectorType::Moench03_old; if (arg == "Unknown") return DetectorType::Unknown; - + throw std::runtime_error("Could not decode detector from: \"" + arg + "\""); } @@ -102,7 +98,8 @@ template <> TimingMode StringTo(const std::string &arg) { return TimingMode::Auto; if (arg == "trigger") return TimingMode::Trigger; - throw std::runtime_error("Could not decode timing mode from: \"" + arg + "\""); + throw std::runtime_error("Could not decode timing mode from: \"" + arg + + "\""); } template <> FrameDiscardPolicy StringTo(const std::string &arg) { @@ -112,7 +109,8 @@ template <> FrameDiscardPolicy StringTo(const std::string &arg) { return FrameDiscardPolicy::Discard; if (arg == "discardpartial") return FrameDiscardPolicy::DiscardPartial; - throw std::runtime_error("Could not decode frame discard policy from: \"" + arg + "\""); + throw std::runtime_error("Could not decode frame discard policy from: \"" + + arg + "\""); } // template <> TimingMode StringTo(std::string mode); diff --git a/src/defs.test.cpp b/src/defs.test.cpp index 6ab9394..2106d86 100644 --- a/src/defs.test.cpp +++ b/src/defs.test.cpp @@ -3,12 +3,12 @@ #include #include -using aare::ToString; using aare::StringTo; +using aare::ToString; TEST_CASE("Enum to string conversion") { - // TODO! By the way I don't think the enum string conversions should be in the defs.hpp file - // but let's use this to show a test + // TODO! By the way I don't think the enum string conversions should be in + // the defs.hpp file but let's use this to show a test REQUIRE(ToString(aare::DetectorType::Generic) == "Generic"); REQUIRE(ToString(aare::DetectorType::Eiger) == "Eiger"); REQUIRE(ToString(aare::DetectorType::Gotthard) == "Gotthard"); @@ -17,30 +17,42 @@ TEST_CASE("Enum to string conversion") { REQUIRE(ToString(aare::DetectorType::Moench) == "Moench"); REQUIRE(ToString(aare::DetectorType::Mythen3) == "Mythen3"); REQUIRE(ToString(aare::DetectorType::Gotthard2) == "Gotthard2"); - REQUIRE(ToString(aare::DetectorType::Xilinx_ChipTestBoard) == "Xilinx_ChipTestBoard"); + REQUIRE(ToString(aare::DetectorType::Xilinx_ChipTestBoard) == + "Xilinx_ChipTestBoard"); REQUIRE(ToString(aare::DetectorType::Moench03) == "Moench03"); REQUIRE(ToString(aare::DetectorType::Moench03_old) == "Moench03_old"); REQUIRE(ToString(aare::DetectorType::Unknown) == "Unknown"); } -TEST_CASE("String to enum"){ - REQUIRE(StringTo("Generic") == aare::DetectorType::Generic); +TEST_CASE("String to enum") { + REQUIRE(StringTo("Generic") == + aare::DetectorType::Generic); REQUIRE(StringTo("Eiger") == aare::DetectorType::Eiger); - REQUIRE(StringTo("Gotthard") == aare::DetectorType::Gotthard); - REQUIRE(StringTo("Jungfrau") == aare::DetectorType::Jungfrau); - REQUIRE(StringTo("ChipTestBoard") == aare::DetectorType::ChipTestBoard); - REQUIRE(StringTo("Moench") == aare::DetectorType::Moench); - REQUIRE(StringTo("Mythen3") == aare::DetectorType::Mythen3); - REQUIRE(StringTo("Gotthard2") == aare::DetectorType::Gotthard2); - REQUIRE(StringTo("Xilinx_ChipTestBoard") == aare::DetectorType::Xilinx_ChipTestBoard); - REQUIRE(StringTo("Moench03") == aare::DetectorType::Moench03); - REQUIRE(StringTo("Moench03_old") == aare::DetectorType::Moench03_old); - REQUIRE(StringTo("Unknown") == aare::DetectorType::Unknown); + REQUIRE(StringTo("Gotthard") == + aare::DetectorType::Gotthard); + REQUIRE(StringTo("Jungfrau") == + aare::DetectorType::Jungfrau); + REQUIRE(StringTo("ChipTestBoard") == + aare::DetectorType::ChipTestBoard); + REQUIRE(StringTo("Moench") == + aare::DetectorType::Moench); + REQUIRE(StringTo("Mythen3") == + aare::DetectorType::Mythen3); + REQUIRE(StringTo("Gotthard2") == + aare::DetectorType::Gotthard2); + REQUIRE(StringTo("Xilinx_ChipTestBoard") == + aare::DetectorType::Xilinx_ChipTestBoard); + REQUIRE(StringTo("Moench03") == + aare::DetectorType::Moench03); + REQUIRE(StringTo("Moench03_old") == + aare::DetectorType::Moench03_old); + REQUIRE(StringTo("Unknown") == + aare::DetectorType::Unknown); } -TEST_CASE("Enum values"){ - //Since some of the enums are written to file we need to make sure - //they match the value in the slsDetectorPackage +TEST_CASE("Enum values") { + // Since some of the enums are written to file we need to make sure + // they match the value in the slsDetectorPackage REQUIRE(static_cast(aare::DetectorType::Generic) == 0); REQUIRE(static_cast(aare::DetectorType::Eiger) == 1); @@ -52,7 +64,7 @@ TEST_CASE("Enum values"){ REQUIRE(static_cast(aare::DetectorType::Gotthard2) == 7); REQUIRE(static_cast(aare::DetectorType::Xilinx_ChipTestBoard) == 8); - //Not included + // Not included REQUIRE(static_cast(aare::DetectorType::Moench03) == 100); } @@ -85,5 +97,6 @@ TEST_CASE("DynamicCluster creation") { // double v3 = c2.get(33 * 44 - 1); // REQUIRE(aare::compare_floats(123.11, v3)); -// REQUIRE_THROWS_AS(c2.set(0, 1), std::invalid_argument); // set int to double +// REQUIRE_THROWS_AS(c2.set(0, 1), std::invalid_argument); // set int to +// double // } \ No newline at end of file diff --git a/src/geo_helpers.cpp b/src/geo_helpers.cpp index 39086ec..96a9056 100644 --- a/src/geo_helpers.cpp +++ b/src/geo_helpers.cpp @@ -2,15 +2,16 @@ #include "aare/geo_helpers.hpp" #include "fmt/core.h" -namespace aare{ +namespace aare { DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) { - #ifdef AARE_VERBOSE +#ifdef AARE_VERBOSE fmt::println("update_geometry_with_roi() called with ROI: {} {} {} {}", roi.xmin, roi.xmax, roi.ymin, roi.ymax); - fmt::println("Geometry: {} {} {} {} {} {}", - geo.modules_x, geo.modules_y, geo.pixels_x, geo.pixels_y, geo.module_gap_row, geo.module_gap_col); - #endif + fmt::println("Geometry: {} {} {} {} {} {}", geo.modules_x, geo.modules_y, + geo.pixels_x, geo.pixels_y, geo.module_gap_row, + geo.module_gap_col); +#endif int pos_y = 0; int pos_y_increment = 0; for (int row = 0; row < geo.modules_y; row++) { @@ -41,9 +42,9 @@ DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) { if (m.origin_y + m.height < roi.ymin) { m.height = 0; } else { - if ((roi.ymin > m.origin_y) && (roi.ymin < m.origin_y + m.height)) { + if ((roi.ymin > m.origin_y) && + (roi.ymin < m.origin_y + m.height)) { m.height -= roi.ymin - m.origin_y; - } if (roi.ymax < m.origin_y + original_height) { m.height -= m.origin_y + original_height - roi.ymax; @@ -51,9 +52,10 @@ DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) { m.origin_y = pos_y; pos_y_increment = m.height; } - #ifdef AARE_VERBOSE - fmt::println("Module {} {} {} {}", m.origin_x, m.origin_y, m.width, m.height); - #endif +#ifdef AARE_VERBOSE + fmt::println("Module {} {} {} {}", m.origin_x, m.origin_y, m.width, + m.height); +#endif } // increment pos_y pos_y += pos_y_increment; @@ -65,7 +67,6 @@ DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) { geo.pixels_y = roi.height(); return geo; - } } // namespace aare \ No newline at end of file diff --git a/src/geo_helpers.test.cpp b/src/geo_helpers.test.cpp index 08ee96c..48ae9cf 100644 --- a/src/geo_helpers.test.cpp +++ b/src/geo_helpers.test.cpp @@ -1,6 +1,6 @@ #include "aare/File.hpp" -#include "aare/RawMasterFile.hpp" //needed for ROI #include "aare/RawFile.hpp" +#include "aare/RawMasterFile.hpp" //needed for ROI #include #include @@ -8,26 +8,24 @@ #include "aare/geo_helpers.hpp" #include "test_config.hpp" -TEST_CASE("Simple ROIs on one module"){ - // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) +TEST_CASE("Simple ROIs on one module") { + // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI + // roi) aare::DetectorGeometry geo; - - aare::ModuleGeometry mod; + aare::ModuleGeometry mod; mod.origin_x = 0; mod.origin_y = 0; mod.width = 1024; mod.height = 512; - - geo.pixels_x = 1024; geo.pixels_y = 512; geo.modules_x = 1; geo.modules_y = 1; geo.module_pixel_0.push_back(mod); - SECTION("ROI is the whole module"){ + SECTION("ROI is the whole module") { aare::ROI roi; roi.xmin = 0; roi.xmax = 1024; @@ -42,7 +40,7 @@ TEST_CASE("Simple ROIs on one module"){ REQUIRE(updated_geo.module_pixel_0[0].height == 512); REQUIRE(updated_geo.module_pixel_0[0].width == 1024); } - SECTION("ROI is the top left corner of the module"){ + SECTION("ROI is the top left corner of the module") { aare::ROI roi; roi.xmin = 100; roi.xmax = 200; @@ -58,7 +56,7 @@ TEST_CASE("Simple ROIs on one module"){ REQUIRE(updated_geo.module_pixel_0[0].width == 100); } - SECTION("ROI is a small square"){ + SECTION("ROI is a small square") { aare::ROI roi; roi.xmin = 1000; roi.xmax = 1010; @@ -73,7 +71,7 @@ TEST_CASE("Simple ROIs on one module"){ REQUIRE(updated_geo.module_pixel_0[0].height == 10); REQUIRE(updated_geo.module_pixel_0[0].width == 10); } - SECTION("ROI is a few columns"){ + SECTION("ROI is a few columns") { aare::ROI roi; roi.xmin = 750; roi.xmax = 800; @@ -90,14 +88,12 @@ TEST_CASE("Simple ROIs on one module"){ } } - - -TEST_CASE("Two modules side by side"){ - // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) +TEST_CASE("Two modules side by side") { + // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI + // roi) aare::DetectorGeometry geo; - - aare::ModuleGeometry mod; + aare::ModuleGeometry mod; mod.origin_x = 0; mod.origin_y = 0; mod.width = 1024; @@ -112,7 +108,7 @@ TEST_CASE("Two modules side by side"){ mod.origin_x = 1024; geo.module_pixel_0.push_back(mod); - SECTION("ROI is the whole image"){ + SECTION("ROI is the whole image") { aare::ROI roi; roi.xmin = 0; roi.xmax = 2048; @@ -125,7 +121,7 @@ TEST_CASE("Two modules side by side"){ REQUIRE(updated_geo.modules_x == 2); REQUIRE(updated_geo.modules_y == 1); } - SECTION("rectangle on both modules"){ + SECTION("rectangle on both modules") { aare::ROI roi; roi.xmin = 800; roi.xmax = 1300; @@ -141,11 +137,12 @@ TEST_CASE("Two modules side by side"){ REQUIRE(updated_geo.module_pixel_0[0].width == 224); REQUIRE(updated_geo.module_pixel_0[1].height == 299); REQUIRE(updated_geo.module_pixel_0[1].width == 276); - } + } } -TEST_CASE("Three modules side by side"){ - // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) +TEST_CASE("Three modules side by side") { + // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI + // roi) aare::DetectorGeometry geo; aare::ROI roi; roi.xmin = 700; @@ -153,7 +150,7 @@ TEST_CASE("Three modules side by side"){ roi.ymin = 0; roi.ymax = 123; - aare::ModuleGeometry mod; + aare::ModuleGeometry mod; mod.origin_x = 0; mod.origin_y = 0; mod.width = 1024; @@ -184,8 +181,9 @@ TEST_CASE("Three modules side by side"){ REQUIRE(updated_geo.module_pixel_0[2].width == 452); } -TEST_CASE("Four modules as a square"){ - // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) +TEST_CASE("Four modules as a square") { + // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI + // roi) aare::DetectorGeometry geo; aare::ROI roi; roi.xmin = 500; @@ -193,7 +191,7 @@ TEST_CASE("Four modules as a square"){ roi.ymin = 500; roi.ymax = 600; - aare::ModuleGeometry mod; + aare::ModuleGeometry mod; mod.origin_x = 0; mod.origin_y = 0; mod.width = 1024; diff --git a/src/utils/ifstream_helpers.cpp b/src/utils/ifstream_helpers.cpp index 74c56f3..79f006d 100644 --- a/src/utils/ifstream_helpers.cpp +++ b/src/utils/ifstream_helpers.cpp @@ -10,7 +10,7 @@ std::string ifstream_error_msg(std::ifstream &ifs) { return " Bad file stream"; } else if (state & std::ios_base::failbit) { return " File read failed"; - }else{ + } else { return " Unknown/no error"; } } diff --git a/src/utils/task.test.cpp b/src/utils/task.test.cpp index e19994a..3ca71c7 100644 --- a/src/utils/task.test.cpp +++ b/src/utils/task.test.cpp @@ -1,10 +1,9 @@ #include "aare/utils/task.hpp" -#include #include +#include - -TEST_CASE("Split a range into multiple tasks"){ +TEST_CASE("Split a range into multiple tasks") { auto tasks = aare::split_task(0, 10, 3); REQUIRE(tasks.size() == 3); @@ -22,11 +21,8 @@ TEST_CASE("Split a range into multiple tasks"){ tasks = aare::split_task(0, 10, 10); REQUIRE(tasks.size() == 10); - for (int i = 0; i < 10; i++){ + for (int i = 0; i < 10; i++) { REQUIRE(tasks[i].first == i); - REQUIRE(tasks[i].second == i+1); + REQUIRE(tasks[i].second == i + 1); } - - - } \ No newline at end of file diff --git a/tests/test.cpp b/tests/test.cpp index 513f690..1d4456f 100644 --- a/tests/test.cpp +++ b/tests/test.cpp @@ -2,8 +2,8 @@ #include #include #include -#include #include +#include TEST_CASE("Test suite can find data assets", "[.integration]") { auto fpath = test_data_path() / "numpy" / "test_numpy_file.npy";