merge from formatted main

This commit is contained in:
2025-06-11 15:11:47 +02:00
93 changed files with 1861 additions and 1618 deletions

View File

@ -2,7 +2,10 @@ name: Build the package using cmake then documentation
on: on:
workflow_dispatch: workflow_dispatch:
push: pull_request:
release:
types:
- published
permissions: permissions:
@ -55,7 +58,7 @@ jobs:
url: ${{ steps.deployment.outputs.page_url }} url: ${{ steps.deployment.outputs.page_url }}
runs-on: ubuntu-latest runs-on: ubuntu-latest
needs: build needs: build
if: github.event_name == 'release' && github.event.action == 'published' if: (github.event_name == 'release' && github.event.action == 'published') || (github.event_name == 'workflow_dispatch' )
steps: steps:
- name: Deploy to GitHub Pages - name: Deploy to GitHub Pages
id: deployment id: deployment

22
RELEASE.md Normal file
View File

@ -0,0 +1,22 @@
# Release notes
### head
Features:
- Cluster finder now works with 5x5, 7x7 and 9x9 clusters
### 2025.05.22
Features:
- Added scurve fitting
Bugfixes:
- Fixed crash when opening raw files with large number of data files

View File

@ -41,8 +41,8 @@ BENCHMARK_F(ClusterFixture, Calculate2x2Eta)(benchmark::State &st) {
} }
// almost takes double the time // almost takes double the time
BENCHMARK_F(ClusterFixture, BENCHMARK_F(ClusterFixture, CalculateGeneralEtaFor2x2Cluster)
CalculateGeneralEtaFor2x2Cluster)(benchmark::State &st) { (benchmark::State &st) {
for (auto _ : st) { for (auto _ : st) {
// This code gets timed // This code gets timed
Eta2 eta = calculate_eta2<int, 2, 2>(cluster_2x2); Eta2 eta = calculate_eta2<int, 2, 2>(cluster_2x2);
@ -59,8 +59,8 @@ BENCHMARK_F(ClusterFixture, Calculate3x3Eta)(benchmark::State &st) {
} }
// almost takes double the time // almost takes double the time
BENCHMARK_F(ClusterFixture, BENCHMARK_F(ClusterFixture, CalculateGeneralEtaFor3x3Cluster)
CalculateGeneralEtaFor3x3Cluster)(benchmark::State &st) { (benchmark::State &st) {
for (auto _ : st) { for (auto _ : st) {
// This code gets timed // This code gets timed
Eta2 eta = calculate_eta2<int, 3, 3>(cluster_3x3); Eta2 eta = calculate_eta2<int, 3, 3>(cluster_3x3);

View File

@ -1,6 +1,5 @@
#include <benchmark/benchmark.h>
#include "aare/NDArray.hpp" #include "aare/NDArray.hpp"
#include <benchmark/benchmark.h>
using aare::NDArray; using aare::NDArray;
@ -22,9 +21,6 @@ public:
// } // }
}; };
BENCHMARK_F(TwoArrays, AddWithOperator)(benchmark::State &st) { BENCHMARK_F(TwoArrays, AddWithOperator)(benchmark::State &st) {
for (auto _ : st) { for (auto _ : st) {
// This code gets timed // This code gets timed

47
docs/src/Philosophy.rst Normal file
View File

@ -0,0 +1,47 @@
****************
Philosophy
****************
Fast code with a simple interface
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Aare should be fast and efficient, but also easy to use. We strive to keep a simple interface that feels intuitive.
Internally we use C++ for performance and the ability to integrate the library in other programs, but we see most
users using the Python interface.
Live at head
~~~~~~~~~~~~~~~~~~
As a user of the library you should be able to, and is expected to, use the latest version. Bug fixes will rarely be backported
to older releases. By upgrading frequently you will benefit from the latest features and minimize the effort to maintain your scripts/code
by doing several small upgrades instead of one big upgrade.
API
~~~~~~~~~~~~~~~~~~
We aim to keep the API stable and only break it for good reasons. But specially now in the early stages of development
the API will change. On those occasions it will be clearly stated in the release notes. However, the norm should be a
backward compatible API.
Documentation
~~~~~~~~~~~~~~~~~~
Being a library it is important to have a well documented API. We use Doxygen to generate the C++ documentation
and Sphinx for the Python part. Breathe is used to integrate the two into one Sphinx html site. The documentation is built
automatically on release by the CI and published to GitHub pages. In addition to the generated API documentation,
certain classes might need more descriptions of the usage. This is then placed in the .rst files in the docs/src directory.
.. attention::
The code should be well documented, but using descriptive names is more important. In the same spirit
if a function is called `getNumberOfFrames()` you don't need to write a comment saying that it gets the
number of frames.
Dependencies
~~~~~~~~~~~~~~~~~~
Deployment in the scientific community is often tricky. Either due to old OS versions or the lack of package managers.
We strive to keep the dependencies to a minimum and will vendor some libraries to simplify deployment even though it comes
at a cost of build time.

View File

@ -2,18 +2,21 @@ Requirements
============================================== ==============================================
- C++17 compiler (gcc 8/clang 7) - C++17 compiler (gcc 8/clang 7)
- CMake 3.14+ - CMake 3.15+
**Internally used libraries** **Internally used libraries**
.. note :: .. note ::
These can also be picked up from the system/conda environment by specifying: To save compile time some of the dependencies can also be picked up from the system/conda environment by specifying:
-DAARE_SYSTEM_LIBRARIES=ON during the cmake configuration. -DAARE_SYSTEM_LIBRARIES=ON during the cmake configuration.
- pybind11 To simplify deployment we build and statically link a few libraries.
- fmt - fmt
- lmfit - https://jugit.fz-juelich.de/mlz/lmfit
- nlohmann_json - nlohmann_json
- pybind11
- ZeroMQ - ZeroMQ
**Extra dependencies for building documentation** **Extra dependencies for building documentation**

86
docs/src/Workflow.rst Normal file
View File

@ -0,0 +1,86 @@
****************
Workflow
****************
This page describes how we develop aare.
GitHub centric
~~~~~~~~~~~~~~~~~~
We use GitHub for all development. Issues and pull requests provide a platform for collaboration as well
as a record of the development process. Even if we discuss things in person, we record the outcome in an issue.
If a particular implementation is chosen over another, the reason should be recorded in the pull request.
Branches
~~~~~~~~~~~~~~~~~~
We aim for an as lightweight branching strategy as possible. Short-lived feature branches are merged back into main.
The main branch is expected to always be in a releasable state. A release is simply a tag on main which provides a
reference and triggers the CI to build the release artifacts (conda, pypi etc.). For large features consider merging
smaller chunks into main as they are completed, rather than waiting for the entire feature to be finished. Worst case
make sure your feature branch merges with main regularly to avoid large merge conflicts later on.
.. note::
The main branch is expected to always work. Feel free to pull from main instead of sticking to a
release
Releases
~~~~~~~~~~~~~~~~~~
Release early, release often. As soon as "enough" new features have been implemented, a release is created.
A release should not be a big thing, rather a routine part of development that does not require any special person or
unfamiliar steps.
Checklists for deployment
~~~~~~~~~~~~~~~~~~
**Feature:**
#. Create a new issue for the feature (label feature)
#. Create a new branch from main.
#. Implement the feature including test and documentation
#. Add the feature to RELEASE.md under head
#. Create a pull request linked to the issue
#. Code is reviewed by at least one other person
#. Once approved, the branch is merged into main
**BugFix:**
Essentially the same as for a feature, if possible start with
a failing test that demonstrates the bug.
#. Create a new issue for the bug (label bug)
#. Create a new branch from main.
#. **Write a test that fails for the bug**
#. Implement the fix
#. **Run the test to ensure it passes**
#. Add the bugfix to RELEASE.md under head
#. Create a pull request linked to the issue.
#. Code is reviewed by at least one other person
#. Once approved, the branch is merged into main
**Release:**
#. Once "enough" new features have been implemented, a release is created
#. Update RELEASE.md with the tag of the release and verify that it is complete
#. Create the release in GitHub describing the new features and bug fixes
#. CI makes magic
**Update documentation only:**
.. attention::
It's possible to update the documentation without changing the code, but take
care since the docs will reflect the code in main and not the latest release.
#. Create a PR to main with the documentation changes
#. Create a pull request linked to the issue.
#. Code is reviewed by at least one other person
#. Once merged you can manually trigger the CI workflow for documentation

View File

@ -67,4 +67,6 @@ AARE
:caption: Developer :caption: Developer
:maxdepth: 3 :maxdepth: 3
Philosophy
Workflow
Tests Tests

View File

@ -1,10 +1,9 @@
#pragma once #pragma once
#include <cstdint> #include "aare/defs.hpp"
#include <cstddef>
#include <array> #include <array>
#include <cassert> #include <cassert>
#include "aare/defs.hpp" #include <cstddef>
#include <cstdint>
namespace aare { namespace aare {
@ -15,7 +14,9 @@ template <typename E, ssize_t Ndim> class ArrayExpr {
auto operator[](size_t i) const { return static_cast<E const &>(*this)[i]; } auto operator[](size_t i) const { return static_cast<E const &>(*this)[i]; }
auto operator()(size_t i) const { return static_cast<E const &>(*this)[i]; } auto operator()(size_t i) const { return static_cast<E const &>(*this)[i]; }
auto size() const { return static_cast<E const &>(*this).size(); } auto size() const { return static_cast<E const &>(*this).size(); }
std::array<ssize_t, Ndim> shape() const { return static_cast<E const &>(*this).shape(); } std::array<ssize_t, Ndim> shape() const {
return static_cast<E const &>(*this).shape();
}
}; };
template <typename A, typename B, ssize_t Ndim> template <typename A, typename B, ssize_t Ndim>
@ -74,8 +75,6 @@ class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); } std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
}; };
template <typename A, typename B, ssize_t Ndim> template <typename A, typename B, ssize_t Ndim>
auto operator+(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) { auto operator+(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArrayAdd<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2); return ArrayAdd<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
@ -96,6 +95,4 @@ auto operator/(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArrayDiv<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2); return ArrayDiv<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
} }
} // namespace aare } // namespace aare

View File

@ -17,7 +17,8 @@ template <class ItemType> class CircularFifo {
public: public:
CircularFifo() : CircularFifo(100){}; 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??? // TODO! how do we deal with alignment for writing? alignas???
// Do we give the user a chance to provide memory locations? // Do we give the user a chance to provide memory locations?
@ -55,7 +56,8 @@ template <class ItemType> class CircularFifo {
bool try_pop_free(ItemType &v) { return free_slots.read(v); } bool try_pop_free(ItemType &v) { return free_slots.read(v); }
ItemType pop_value(std::chrono::nanoseconds wait, std::atomic<bool> &stopped) { ItemType pop_value(std::chrono::nanoseconds wait,
std::atomic<bool> &stopped) {
ItemType v; ItemType v;
while (!filled_slots.read(v) && !stopped) { while (!filled_slots.read(v) && !stopped) {
std::this_thread::sleep_for(wait); std::this_thread::sleep_for(wait);

View File

@ -371,8 +371,8 @@ ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
"Could not read number of clusters"); "Could not read number of clusters");
} }
LOG(logDEBUG1) << "Reading " << n_clusters LOG(logDEBUG1) << "Reading " << n_clusters << " clusters from frame "
<< " clusters from frame " << frame_number; << frame_number;
ClusterVector<ClusterType> clusters(n_clusters); ClusterVector<ClusterType> clusters(n_clusters);
clusters.set_frame_number(frame_number); clusters.set_frame_number(frame_number);

View File

@ -41,8 +41,7 @@ class ClusterFinder {
m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) { m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {
LOG(logDEBUG) << "ClusterFinder: " LOG(logDEBUG) << "ClusterFinder: "
<< "image_size: " << image_size[0] << "x" << image_size[1] << "image_size: " << image_size[0] << "x" << image_size[1]
<< ", nSigma: " << nSigma << ", nSigma: " << nSigma << ", capacity: " << capacity;
<< ", capacity: " << capacity;
} }
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) { void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {

View File

@ -7,8 +7,8 @@
#include "aare/ClusterFinder.hpp" #include "aare/ClusterFinder.hpp"
#include "aare/NDArray.hpp" #include "aare/NDArray.hpp"
#include "aare/logger.hpp"
#include "aare/ProducerConsumerQueue.hpp" #include "aare/ProducerConsumerQueue.hpp"
#include "aare/logger.hpp"
namespace aare { namespace aare {
@ -125,8 +125,8 @@ class ClusterFinderMT {
: m_n_threads(n_threads) { : m_n_threads(n_threads) {
LOG(logDEBUG1) << "ClusterFinderMT: " LOG(logDEBUG1) << "ClusterFinderMT: "
<< "image_size: " << image_size[0] << "x" << image_size[1] << "image_size: " << image_size[0] << "x"
<< ", nSigma: " << nSigma << image_size[1] << ", nSigma: " << nSigma
<< ", capacity: " << capacity << ", capacity: " << capacity
<< ", n_threads: " << n_threads; << ", n_threads: " << n_threads;

View File

@ -1,21 +1,21 @@
#pragma once #pragma once
#include "aare/FileInterface.hpp" #include "aare/FileInterface.hpp"
#include "aare/RawMasterFile.hpp"
#include "aare/Frame.hpp" #include "aare/Frame.hpp"
#include "aare/RawMasterFile.hpp"
#include <filesystem> #include <filesystem>
#include <fstream> #include <fstream>
namespace aare { namespace aare {
class CtbRawFile { class CtbRawFile {
RawMasterFile m_master; RawMasterFile m_master;
std::ifstream m_file; std::ifstream m_file;
size_t m_current_frame{0}; size_t m_current_frame{0};
size_t m_current_subfile{0}; size_t m_current_subfile{0};
size_t m_num_subfiles{0}; size_t m_num_subfiles{0};
public: public:
CtbRawFile(const std::filesystem::path &fname); CtbRawFile(const std::filesystem::path &fname);
@ -29,13 +29,13 @@ public:
size_t frames_in_file() const; size_t frames_in_file() const;
RawMasterFile master() const; RawMasterFile master() const;
private: private:
void find_subfiles(); void find_subfiles();
size_t sub_file_index(size_t frame_index) const { size_t sub_file_index(size_t frame_index) const {
return frame_index / m_master.max_frames_per_file(); return frame_index / m_master.max_frames_per_file();
} }
void open_data_file(size_t subfile_index); void open_data_file(size_t subfile_index);
}; };
} } // namespace aare

View File

@ -6,31 +6,37 @@
namespace aare { 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 // - python documentation: https://docs.python.org/3/c-api/arg.html#numbers
// - py::format_descriptor<T>::format() (in pybind11) does not return the same format as // - py::format_descriptor<T>::format() (in pybind11) does not return the same
// format as
// written in python.org documentation. // written in python.org documentation.
// - numpy also doesn't use the same format. and also numpy associates the format // - numpy also doesn't use the same format. and also numpy associates the
// with variable bitdepth types. (e.g. long is int64 on linux64 and int32 on win64) // format
// https://numpy.org/doc/stable/reference/arrays.scalars.html // 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: // github issue discussing this:
// https://github.com/pybind/pybind11/issues/1908#issuecomment-658358767 // https://github.com/pybind/pybind11/issues/1908#issuecomment-658358767
// //
// [IN LINUX] the difference is for int64 (long) and uint64 (unsigned long). The format // [IN LINUX] the difference is for int64 (long) and uint64 (unsigned long). The
// descriptor is 'q' and 'Q' respectively and in the documentation it is 'l' and 'k'. // 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 // 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. // 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 // for this reason we decided to use the same format descriptor as pybind to
// any further discrepancies. // avoid any further discrepancies.
// in the following order: // in the following order:
// int8, uint8, int16, uint16, int32, uint32, int64, uint64, float, double // 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 // 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 * @brief enum class to define the endianess of the system
*/ */
@ -52,12 +58,29 @@ enum class endian {
*/ */
class Dtype { class Dtype {
public: 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; uint8_t bitdepth() const;
size_t bytes() const; size_t bytes() const;
std::string format_descr() const { return std::string(1, DTYPE_FORMAT_DSC[static_cast<int>(m_type)]); } std::string format_descr() const {
std::string numpy_descr() const { return std::string(1, NUMPY_FORMAT_DSC[static_cast<int>(m_type)]); } return std::string(1, DTYPE_FORMAT_DSC[static_cast<int>(m_type)]);
}
std::string numpy_descr() const {
return std::string(1, NUMPY_FORMAT_DSC[static_cast<int>(m_type)]);
}
explicit Dtype(const std::type_info &t); explicit Dtype(const std::type_info &t);
explicit Dtype(std::string_view sv); explicit Dtype(std::string_view sv);

View File

@ -7,7 +7,11 @@ namespace aare {
/** /**
* @brief RAII File class for reading, and in the future potentially writing * @brief RAII File class for reading, and in the future potentially writing
* image files in various formats. Minimal generic interface. For specail * image files in various formats. Minimal generic interface. For specail
<<<<<<< HEAD
* fuctions plase use the RawFile, NumpyFile or Hdf5File classes directly. Wraps * fuctions plase use the RawFile, NumpyFile or Hdf5File classes directly. Wraps
=======
* fuctions plase use the RawFile or NumpyFile classes directly. Wraps
>>>>>>> main
* FileInterface to abstract the underlying file format * FileInterface to abstract the underlying file format
* @note **frame_number** refers the the frame number sent by the detector while * @note **frame_number** refers the the frame number sent by the detector while
* **frame_index** is the position of the frame in the file * **frame_index** is the position of the frame in the file
@ -25,9 +29,11 @@ class File {
* @throws std::invalid_argument if the file mode is not supported * @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 = {}); 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 */ /**Since the object is responsible for managing the file we disable copy
* construction */
File(File const &other) = delete; File(File const &other) = delete;
/**The same goes for copy assignment */ /**The same goes for copy assignment */
@ -39,15 +45,19 @@ class File {
// void close(); //!< close the file // void close(); //!< close the file
Frame read_frame(); //!< read one frame from the file at the current position Frame
Frame read_frame(size_t frame_index); //!< read one frame at the position given by frame number read_frame(); //!< read one frame from the file at the current position
std::vector<Frame> read_n(size_t n_frames); //!< read n_frames 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<Frame> 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);
void read_into(std::byte *image_buf, size_t n_frames); 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(); //!< 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 frame_number(
size_t frame_index); //!< get the frame number at the given frame index
size_t bytes_per_frame() const; size_t bytes_per_frame() const;
size_t pixels_per_frame() const; size_t pixels_per_frame() const;
size_t bytes_per_pixel() const; size_t bytes_per_pixel() const;
@ -59,8 +69,6 @@ class File {
size_t cols() const; size_t cols() const;
DetectorType detector_type() const; DetectorType detector_type() const;
}; };
} // namespace aare } // namespace aare

View File

@ -20,8 +20,10 @@ struct FileConfig {
uint64_t rows{}; uint64_t rows{};
uint64_t cols{}; uint64_t cols{};
bool operator==(const FileConfig &other) const { bool operator==(const FileConfig &other) const {
return dtype == other.dtype && rows == other.rows && cols == other.cols && geometry == other.geometry && return dtype == other.dtype && rows == other.rows &&
detector_type == other.detector_type && max_frames_per_file == other.max_frames_per_file; 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); } bool operator!=(const FileConfig &other) const { return !(*this == other); }
@ -32,8 +34,11 @@ struct FileConfig {
int max_frames_per_file{}; int max_frames_per_file{};
size_t total_frames{}; size_t total_frames{};
std::string to_string() const { std::string to_string() const {
return "{ dtype: " + dtype.to_string() + ", rows: " + std::to_string(rows) + ", cols: " + std::to_string(cols) + return "{ dtype: " + dtype.to_string() +
", geometry: " + geometry.to_string() + ", detector_type: " + ToString(detector_type) + ", 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) + ", max_frames_per_file: " + std::to_string(max_frames_per_file) +
", total_frames: " + std::to_string(total_frames) + " }"; ", total_frames: " + std::to_string(total_frames) + " }";
} }
@ -65,17 +70,20 @@ class FileInterface {
* @param n_frames number of frames to read * @param n_frames number of frames to read
* @return vector of frames * @return vector of frames
*/ */
virtual std::vector<Frame> read_n(size_t n_frames) = 0; // Is this the right interface? virtual std::vector<Frame>
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 * @param image_buf buffer to store the frame
* @return void * @return void
*/ */
virtual void read_into(std::byte *image_buf) = 0; 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 image_buf buffer to store the frames
* @param n_frames number of frames to read * @param n_frames number of frames to read
* @return void * @return void
@ -135,7 +143,6 @@ class FileInterface {
*/ */
virtual size_t bitdepth() const = 0; virtual size_t bitdepth() const = 0;
virtual DetectorType detector_type() const = 0; virtual DetectorType detector_type() const = 0;
// function to query the data type of the file // function to query the data type of the file

View File

@ -23,16 +23,19 @@ NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par);
} // namespace func } // namespace func
/** /**
* @brief Estimate the initial parameters for a Gaussian fit * @brief Estimate the initial parameters for a Gaussian fit
*/ */
std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<double, 1> y); std::array<double, 3> gaus_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
std::array<double, 2> pol1_init_par(const NDView<double, 1> x, const NDView<double, 1> y); std::array<double, 2> pol1_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
std::array<double, 6> scurve_init_par(const NDView<double, 1> x, const NDView<double, 1> y); std::array<double, 6> scurve_init_par(const NDView<double, 1> x,
std::array<double, 6> scurve2_init_par(const NDView<double, 1> x, const NDView<double, 1> y); const NDView<double, 1> y);
std::array<double, 6> scurve2_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
static constexpr int DEFAULT_NUM_THREADS = 4; static constexpr int DEFAULT_NUM_THREADS = 4;
@ -43,7 +46,6 @@ static constexpr int DEFAULT_NUM_THREADS = 4;
*/ */
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y); NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
/** /**
* @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values] * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
* @param x x values * @param x x values
@ -54,9 +56,6 @@ NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
int n_threads = DEFAULT_NUM_THREADS); int n_threads = DEFAULT_NUM_THREADS);
/** /**
* @brief Fit a 1D Gaussian with error estimates * @brief Fit a 1D Gaussian with error estimates
* @param x x values * @param x x values
@ -80,9 +79,8 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
* @param n_threads number of threads to use * @param n_threads number of threads to use
*/ */
void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err, void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out, NDView<double, 3> par_out, NDView<double, 3> par_err_out,
int n_threads = DEFAULT_NUM_THREADS NDView<double, 2> chi2_out, int n_threads = DEFAULT_NUM_THREADS);
);
NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y); NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y);
@ -90,26 +88,33 @@ NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
int n_threads = DEFAULT_NUM_THREADS); int n_threads = DEFAULT_NUM_THREADS);
void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err, void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2); NDView<double, 1> par_out, NDView<double, 1> par_err_out,
double &chi2);
// TODO! not sure we need to offer the different version in C++ // TODO! not sure we need to offer the different version in C++
void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err, void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out,NDView<double, 2> chi2_out, NDView<double, 3> par_out, NDView<double, 3> par_err_out,
int n_threads = DEFAULT_NUM_THREADS); NDView<double, 2> chi2_out, int n_threads = DEFAULT_NUM_THREADS);
NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y); NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_threads); NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err, int n_threads);
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2); void fit_scurve(NDView<double, 1> x, NDView<double, 1> y,
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err, NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out, NDView<double, 1> par_err_out, double &chi2);
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads); int n_threads);
NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y); NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, int n_threads); NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err, int n_threads);
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2); void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y,
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err, NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out, NDView<double, 1> par_err_out, double &chi2);
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads); int n_threads);
} // namespace aare } // namespace aare

View File

@ -52,7 +52,6 @@ class Frame {
Frame &operator=(Frame &&other) noexcept; Frame &operator=(Frame &&other) noexcept;
Frame(Frame &&other) noexcept; Frame(Frame &&other) noexcept;
Frame clone() const; //<- Explicit copy Frame clone() const; //<- Explicit copy
uint32_t rows() const; uint32_t rows() const;

View File

@ -3,13 +3,12 @@
#include <filesystem> #include <filesystem>
#include <vector> #include <vector>
#include "aare/FilePtr.hpp"
#include "aare/defs.hpp"
#include "aare/NDArray.hpp"
#include "aare/FileInterface.hpp" #include "aare/FileInterface.hpp"
#include "aare/FilePtr.hpp"
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
namespace aare { namespace aare {
struct JungfrauDataHeader { struct JungfrauDataHeader {
uint64_t framenum; uint64_t framenum;
uint64_t bunchid; uint64_t bunchid;
@ -18,33 +17,39 @@ struct JungfrauDataHeader{
class JungfrauDataFile : public FileInterface { class JungfrauDataFile : public FileInterface {
size_t m_rows{}; //!< number of rows 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_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_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_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_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_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_current_frame_index{}; //!< The index of the current frame (with
//!< reference to all files)
std::vector<size_t> m_last_frame_in_file{}; //!< Used for seeking to the correct file std::vector<size_t>
m_last_frame_in_file{}; //!< Used for seeking to the correct file
std::filesystem::path m_path; //!< path to the files std::filesystem::path m_path; //!< path to the files
std::string m_base_name; //!< base name used for formatting file names 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; using pixel_type = uint16_t;
static constexpr size_t header_size = sizeof(JungfrauDataHeader); 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 n_digits_in_file_index =
6; //!< to format file names
public: public:
JungfrauDataFile(const std::filesystem::path &fname); JungfrauDataFile(const std::filesystem::path &fname);
std::string base_name() const; //!< get the base name of the file (without path and extension) std::string base_name()
const; //!< get the base name of the file (without path and extension)
size_t bytes_per_frame() override; size_t bytes_per_frame() override;
size_t pixels_per_frame() override; size_t pixels_per_frame() override;
size_t bytes_per_pixel() const; size_t bytes_per_pixel() const;
size_t bitdepth() const override; size_t bitdepth() const override;
void seek(size_t frame_index) override; //!< seek to the given frame index (note not byte offset) 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 tell() override; //!< get the frame index of the file pointer
size_t total_frames() const override; size_t total_frames() const override;
size_t rows() const override; size_t rows() const override;
@ -63,44 +68,48 @@ class JungfrauDataFile : public FileInterface {
/** /**
* @brief Read a single frame from the file into the given 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 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) * @param header pointer to a JungfrauDataHeader or nullptr to skip header)
*/ */
void read_into(std::byte *image_buf, JungfrauDataHeader *header = nullptr); void read_into(std::byte *image_buf, JungfrauDataHeader *header = nullptr);
/** /**
* @brief Read a multiple frames from the file into the given 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 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 n_frames number of frames to read
* @param header pointer to a JungfrauDataHeader or nullptr to skip header) * @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 * @brief Read a single frame from the file into the given NDArray
* @param image NDArray to read the frame into. * @param image NDArray to read the frame into.
*/ */
void read_into(NDArray<uint16_t>* image, JungfrauDataHeader* header = nullptr); void read_into(NDArray<uint16_t> *image,
JungfrauDataHeader *header = nullptr);
JungfrauDataHeader read_header(); 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 * @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 parse_fname(const std::filesystem::path &fname);
void scan_files(); void scan_files();
void open_file(size_t file_index); void open_file(size_t file_index);
std::filesystem::path fpath(size_t frame_index) const; std::filesystem::path fpath(size_t frame_index) const;
}; };
} // namespace aare } // namespace aare

View File

@ -21,7 +21,6 @@ TODO! Add expression templates for operators
namespace aare { namespace aare {
template <typename T, ssize_t Ndim = 2> template <typename T, ssize_t Ndim = 2>
class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> { class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
std::array<ssize_t, Ndim> shape_; std::array<ssize_t, Ndim> shape_;
@ -48,7 +47,6 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
std::multiplies<>())), std::multiplies<>())),
data_(new T[size_]) {} data_(new T[size_]) {}
/** /**
* @brief Construct a new NDArray object with a shape and value. * @brief Construct a new NDArray object with a shape and value.
* *
@ -79,7 +77,6 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
: shape_(other.shape_), strides_(c_strides<Ndim>(shape_)), : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
size_(other.size_), data_(other.data_) { size_(other.size_), data_(other.data_) {
other.reset(); // TODO! is this necessary? other.reset(); // TODO! is this necessary?
} }
// Copy constructor // Copy constructor
@ -157,11 +154,6 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
NDArray &operator&=(const T & /*mask*/); NDArray &operator&=(const T & /*mask*/);
void sqrt() { void sqrt() {
for (int i = 0; i < size_; ++i) { for (int i = 0; i < size_; ++i) {
data_[i] = std::sqrt(data_[i]); data_[i] = std::sqrt(data_[i]);
@ -345,9 +337,6 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
return *this; return *this;
} }
template <typename T, ssize_t Ndim> template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) { NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) {
NDArray result = *this; NDArray result = *this;
@ -448,6 +437,4 @@ NDArray<T, Ndim> load(const std::string &pathname,
return img; return img;
} }
} // namespace aare } // namespace aare

View File

@ -1,6 +1,6 @@
#pragma once #pragma once
#include "aare/defs.hpp"
#include "aare/ArrayExpr.hpp" #include "aare/ArrayExpr.hpp"
#include "aare/defs.hpp"
#include <algorithm> #include <algorithm>
#include <array> #include <array>
@ -17,7 +17,8 @@ namespace aare {
template <ssize_t Ndim> using Shape = std::array<ssize_t, Ndim>; template <ssize_t Ndim> using Shape = std::array<ssize_t, Ndim>;
// TODO! fix mismatch between signed and unsigned // TODO! fix mismatch between signed and unsigned
template <ssize_t Ndim> Shape<Ndim> make_shape(const std::vector<size_t> &shape) { template <ssize_t Ndim>
Shape<Ndim> make_shape(const std::vector<size_t> &shape) {
if (shape.size() != Ndim) if (shape.size() != Ndim)
throw std::runtime_error("Shape size mismatch"); throw std::runtime_error("Shape size mismatch");
Shape<Ndim> arr; Shape<Ndim> arr;
@ -25,14 +26,18 @@ template <ssize_t Ndim> Shape<Ndim> make_shape(const std::vector<size_t> &shape)
return arr; return arr;
} }
template <ssize_t Dim = 0, typename Strides> ssize_t element_offset(const Strides & /*unused*/) { return 0; } template <ssize_t Dim = 0, typename Strides>
ssize_t element_offset(const Strides & /*unused*/) {
return 0;
}
template <ssize_t Dim = 0, typename Strides, typename... Ix> template <ssize_t Dim = 0, typename Strides, typename... Ix>
ssize_t element_offset(const Strides &strides, ssize_t i, Ix... index) { ssize_t element_offset(const Strides &strides, ssize_t i, Ix... index) {
return i * strides[Dim] + element_offset<Dim + 1>(strides, index...); return i * strides[Dim] + element_offset<Dim + 1>(strides, index...);
} }
template <ssize_t Ndim> std::array<ssize_t, Ndim> c_strides(const std::array<ssize_t, Ndim> &shape) { template <ssize_t Ndim>
std::array<ssize_t, Ndim> c_strides(const std::array<ssize_t, Ndim> &shape) {
std::array<ssize_t, Ndim> strides{}; std::array<ssize_t, Ndim> strides{};
std::fill(strides.begin(), strides.end(), 1); std::fill(strides.begin(), strides.end(), 1);
for (ssize_t i = Ndim - 1; i > 0; --i) { for (ssize_t i = Ndim - 1; i > 0; --i) {
@ -41,14 +46,16 @@ template <ssize_t Ndim> std::array<ssize_t, Ndim> c_strides(const std::array<ssi
return strides; return strides;
} }
template <ssize_t Ndim> std::array<ssize_t, Ndim> make_array(const std::vector<ssize_t> &vec) { template <ssize_t Ndim>
std::array<ssize_t, Ndim> make_array(const std::vector<ssize_t> &vec) {
assert(vec.size() == Ndim); assert(vec.size() == Ndim);
std::array<ssize_t, Ndim> arr{}; std::array<ssize_t, Ndim> arr{};
std::copy_n(vec.begin(), Ndim, arr.begin()); std::copy_n(vec.begin(), Ndim, arr.begin());
return arr; return arr;
} }
template <typename T, ssize_t Ndim = 2> class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> { template <typename T, ssize_t Ndim = 2>
class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
public: public:
NDView() = default; NDView() = default;
~NDView() = default; ~NDView() = default;
@ -57,17 +64,23 @@ template <typename T, ssize_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
NDView(T *buffer, std::array<ssize_t, Ndim> shape) NDView(T *buffer, std::array<ssize_t, Ndim> shape)
: buffer_(buffer), strides_(c_strides<Ndim>(shape)), shape_(shape), : buffer_(buffer), strides_(c_strides<Ndim>(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<ssize_t> &shape) // NDView(T *buffer, const std::vector<ssize_t> &shape)
// : buffer_(buffer), strides_(c_strides<Ndim>(make_array<Ndim>(shape))), shape_(make_array<Ndim>(shape)), // : buffer_(buffer),
// size_(std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>())) {} // strides_(c_strides<Ndim>(make_array<Ndim>(shape))),
// shape_(make_array<Ndim>(shape)),
// size_(std::accumulate(std::begin(shape), std::end(shape), 1,
// std::multiplies<>())) {}
template <typename... Ix> std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) { template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
return buffer_[element_offset(strides_, index...)]; return buffer_[element_offset(strides_, index...)];
} }
template <typename... Ix> std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const { template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
return buffer_[element_offset(strides_, index...)]; return buffer_[element_offset(strides_, index...)];
} }
@ -94,16 +107,21 @@ template <typename T, ssize_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
NDView &operator+=(const T val) { return elemenwise(val, std::plus<T>()); } NDView &operator+=(const T val) { return elemenwise(val, std::plus<T>()); }
NDView &operator-=(const T val) { return elemenwise(val, std::minus<T>()); } NDView &operator-=(const T val) { return elemenwise(val, std::minus<T>()); }
NDView &operator*=(const T val) { return elemenwise(val, std::multiplies<T>()); } NDView &operator*=(const T val) {
NDView &operator/=(const T val) { return elemenwise(val, std::divides<T>()); } return elemenwise(val, std::multiplies<T>());
}
NDView &operator/=(const T val) {
return elemenwise(val, std::divides<T>());
}
NDView &operator/=(const NDView &other) { return elemenwise(other, std::divides<T>()); } NDView &operator/=(const NDView &other) {
return elemenwise(other, std::divides<T>());
}
template <size_t Size> NDView &operator=(const std::array<T, Size> &arr) {
template<size_t Size>
NDView& operator=(const std::array<T, Size> &arr) {
if (size() != static_cast<ssize_t>(arr.size())) if (size() != static_cast<ssize_t>(arr.size()))
throw std::runtime_error(LOCATION + "Array and NDView size mismatch"); throw std::runtime_error(LOCATION +
"Array and NDView size mismatch");
std::copy(arr.begin(), arr.end(), begin()); std::copy(arr.begin(), arr.end(), begin());
return *this; return *this;
} }
@ -147,13 +165,15 @@ template <typename T, ssize_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
std::array<ssize_t, Ndim> shape_{}; std::array<ssize_t, Ndim> shape_{};
uint64_t size_{}; uint64_t size_{};
template <class BinaryOperation> NDView &elemenwise(T val, BinaryOperation op) { template <class BinaryOperation>
NDView &elemenwise(T val, BinaryOperation op) {
for (uint64_t i = 0; i != size_; ++i) { for (uint64_t i = 0; i != size_; ++i) {
buffer_[i] = op(buffer_[i], val); buffer_[i] = op(buffer_[i], val);
} }
return *this; return *this;
} }
template <class BinaryOperation> NDView &elemenwise(const NDView &other, BinaryOperation op) { template <class BinaryOperation>
NDView &elemenwise(const NDView &other, BinaryOperation op) {
for (uint64_t i = 0; i != size_; ++i) { for (uint64_t i = 0; i != size_; ++i) {
buffer_[i] = op(buffer_[i], other.buffer_[i]); buffer_[i] = op(buffer_[i], other.buffer_[i]);
} }
@ -170,7 +190,6 @@ template <typename T, ssize_t Ndim> void NDView<T, Ndim>::print_all() const {
} }
} }
template <typename T, ssize_t Ndim> template <typename T, ssize_t Ndim>
std::ostream &operator<<(std::ostream &os, const NDView<T, Ndim> &arr) { std::ostream &operator<<(std::ostream &os, const NDView<T, Ndim> &arr) {
for (auto row = 0; row < arr.shape(0); ++row) { for (auto row = 0; row < arr.shape(0); ++row) {
@ -183,9 +202,7 @@ std::ostream& operator <<(std::ostream& os, const NDView<T, Ndim>& arr){
return os; return os;
} }
template <typename T> NDView<T, 1> make_view(std::vector<T> &vec) {
template <typename T>
NDView<T,1> make_view(std::vector<T>& vec){
return NDView<T, 1>(vec.data(), {static_cast<ssize_t>(vec.size())}); return NDView<T, 1>(vec.data(), {static_cast<ssize_t>(vec.size())});
} }

View File

@ -1,9 +1,8 @@
#pragma once #pragma once
#include "aare/Dtype.hpp" #include "aare/Dtype.hpp"
#include "aare/defs.hpp"
#include "aare/FileInterface.hpp" #include "aare/FileInterface.hpp"
#include "aare/NumpyHelpers.hpp" #include "aare/NumpyHelpers.hpp"
#include "aare/defs.hpp"
#include <filesystem> #include <filesystem>
#include <iostream> #include <iostream>
@ -11,13 +10,12 @@
namespace aare { namespace aare {
/** /**
* @brief NumpyFile class to read and write numpy files * @brief NumpyFile class to read and write numpy files
* @note derived from FileInterface * @note derived from FileInterface
* @note implements all the pure virtual functions 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 { class NumpyFile : public FileInterface {
@ -28,26 +26,35 @@ class NumpyFile : public FileInterface {
* @param mode file mode (r, w) * @param mode file mode (r, w)
* @param cfg file configuration * @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); void write(Frame &frame);
Frame read_frame() override { return get_frame(this->current_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<Frame> read_n(size_t n_frames) override; std::vector<Frame> 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; 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 frame_number(size_t frame_index) override { return frame_index; };
size_t bytes_per_frame() override; size_t bytes_per_frame() override;
size_t pixels_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 tell() override { return this->current_frame; }
size_t total_frames() const override { return m_header.shape[0]; } size_t total_frames() const override { return m_header.shape[0]; }
size_t rows() const override { return m_header.shape[1]; } size_t rows() const override { return m_header.shape[1]; }
size_t cols() const override { return m_header.shape[2]; } size_t cols() const override { return m_header.shape[2]; }
size_t bitdepth() const override { return m_header.dtype.bitdepth(); } 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 * @brief get the data type of the numpy file
@ -70,7 +77,8 @@ class NumpyFile : public FileInterface {
template <typename T, size_t NDim> NDArray<T, NDim> load() { template <typename T, size_t NDim> NDArray<T, NDim> load() {
NDArray<T, NDim> arr(make_shape<NDim>(m_header.shape)); NDArray<T, NDim> arr(make_shape<NDim>(m_header.shape));
if (fseek(fp, static_cast<long>(header_size), SEEK_SET)) { if (fseek(fp, static_cast<long>(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); size_t rc = fread(arr.data(), sizeof(T), arr.size(), fp);
if (rc != static_cast<size_t>(arr.size())) { if (rc != static_cast<size_t>(arr.size())) {
@ -78,16 +86,20 @@ class NumpyFile : public FileInterface {
} }
return arr; return arr;
} }
template <typename A, typename TYPENAME, A Ndim> void write(NDView<TYPENAME, Ndim> &frame) { template <typename A, typename TYPENAME, A Ndim>
void write(NDView<TYPENAME, Ndim> &frame) {
write_impl(frame.data(), frame.total_bytes()); write_impl(frame.data(), frame.total_bytes());
} }
template <typename A, typename TYPENAME, A Ndim> void write(NDArray<TYPENAME, Ndim> &frame) { template <typename A, typename TYPENAME, A Ndim>
void write(NDArray<TYPENAME, Ndim> &frame) {
write_impl(frame.data(), frame.total_bytes()); write_impl(frame.data(), frame.total_bytes());
} }
template <typename A, typename TYPENAME, A Ndim> void write(NDView<TYPENAME, Ndim> &&frame) { template <typename A, typename TYPENAME, A Ndim>
void write(NDView<TYPENAME, Ndim> &&frame) {
write_impl(frame.data(), frame.total_bytes()); write_impl(frame.data(), frame.total_bytes());
} }
template <typename A, typename TYPENAME, A Ndim> void write(NDArray<TYPENAME, Ndim> &&frame) { template <typename A, typename TYPENAME, A Ndim>
void write(NDArray<TYPENAME, Ndim> &&frame) {
write_impl(frame.data(), frame.total_bytes()); write_impl(frame.data(), frame.total_bytes());
} }

View File

@ -40,15 +40,18 @@ bool parse_bool(const std::string &in);
std::string get_value_from_map(const std::string &mapstr); std::string get_value_from_map(const std::string &mapstr);
std::unordered_map<std::string, std::string> parse_dict(std::string in, const std::vector<std::string> &keys); std::unordered_map<std::string, std::string>
parse_dict(std::string in, const std::vector<std::string> &keys);
template <typename T, size_t N> bool in_array(T val, const std::array<T, N> &arr) { template <typename T, size_t N>
bool in_array(T val, const std::array<T, N> &arr) {
return std::find(std::begin(arr), std::end(arr), val) != std::end(arr); return std::find(std::begin(arr), std::end(arr), val) != std::end(arr);
} }
bool is_digits(const std::string &str); bool is_digits(const std::string &str);
aare::Dtype parse_descr(std::string typestring); 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); size_t write_header(std::ostream &out, const NumpyHeader &header);
} // namespace NumpyHelpers } // namespace NumpyHelpers

View File

@ -42,9 +42,7 @@ template <typename SUM_TYPE = double> class Pedestal {
} }
~Pedestal() = default; ~Pedestal() = default;
NDArray<SUM_TYPE, 2> mean() { NDArray<SUM_TYPE, 2> mean() { return m_mean; }
return m_mean;
}
SUM_TYPE mean(const uint32_t row, const uint32_t col) const { SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
return m_mean(row, col); return m_mean(row, col);
@ -71,8 +69,6 @@ template <typename SUM_TYPE = double> class Pedestal {
return variance_array; return variance_array;
} }
NDArray<SUM_TYPE, 2> std() { NDArray<SUM_TYPE, 2> std() {
NDArray<SUM_TYPE, 2> standard_deviation_array({m_rows, m_cols}); NDArray<SUM_TYPE, 2> standard_deviation_array({m_rows, m_cols});
for (uint32_t i = 0; i < m_rows * m_cols; i++) { for (uint32_t i = 0; i < m_rows * m_cols; i++) {
@ -83,8 +79,6 @@ template <typename SUM_TYPE = double> class Pedestal {
return standard_deviation_array; return standard_deviation_array;
} }
void clear() { void clear() {
m_sum = 0; m_sum = 0;
m_sum2 = 0; m_sum2 = 0;
@ -92,8 +86,6 @@ template <typename SUM_TYPE = double> class Pedestal {
m_mean = 0; m_mean = 0;
} }
void clear(const uint32_t row, const uint32_t col) { void clear(const uint32_t row, const uint32_t col) {
m_sum(row, col) = 0; m_sum(row, col) = 0;
m_sum2(row, col) = 0; m_sum2(row, col) = 0;
@ -101,8 +93,6 @@ template <typename SUM_TYPE = double> class Pedestal {
m_mean(row, col) = 0; m_mean(row, col) = 0;
} }
template <typename T> void push(NDView<T, 2> frame) { template <typename T> void push(NDView<T, 2> frame) {
assert(frame.size() == m_rows * m_cols); assert(frame.size() == m_rows * m_cols);
@ -140,9 +130,6 @@ template <typename SUM_TYPE = double> class Pedestal {
} }
} }
template <typename T> void push(Frame &frame) { template <typename T> void push(Frame &frame) {
assert(frame.rows() == static_cast<size_t>(m_rows) && assert(frame.rows() == static_cast<size_t>(m_rows) &&
frame.cols() == static_cast<size_t>(m_cols)); frame.cols() == static_cast<size_t>(m_cols));
@ -170,7 +157,8 @@ template <typename SUM_TYPE = double> class Pedestal {
m_sum(row, col) += val - m_sum(row, col) / m_samples; m_sum(row, col) += val - m_sum(row, col) / m_samples;
m_sum2(row, col) += val * val - m_sum2(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); m_mean(row, col) = m_sum(row, col) / m_cur_samples(row, col);
} }
@ -183,7 +171,8 @@ template <typename SUM_TYPE = double> class Pedestal {
m_cur_samples(row, col)++; m_cur_samples(row, col)++;
} else { } else {
m_sum(row, col) += val - m_sum(row, col) / m_cur_samples(row, col); 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,9 +180,7 @@ template <typename SUM_TYPE = double> class Pedestal {
* @brief Update the mean of the pedestal. This is used after having done * @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. * push_no_update. It is not necessary to call this function after push.
*/ */
void update_mean(){ void update_mean() { m_mean = m_sum / m_cur_samples; }
m_mean = m_sum / m_cur_samples;
}
template <typename T> template <typename T>
void push_fast(const uint32_t row, const uint32_t col, const T val_) { void push_fast(const uint32_t row, const uint32_t col, const T val_) {
@ -204,6 +191,5 @@ template <typename SUM_TYPE = double> class Pedestal {
m_sum2(row, col) += val * val - m_sum2(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; m_mean(row, col) = m_sum(row, col) / m_samples;
} }
}; };
} // namespace aare } // namespace aare

View File

@ -1,7 +1,7 @@
#pragma once #pragma once
#include "aare/defs.hpp"
#include "aare/NDArray.hpp" #include "aare/NDArray.hpp"
#include "aare/defs.hpp"
namespace aare { namespace aare {

View File

@ -18,9 +18,9 @@
// @author Jordan DeLong (delong.j@fb.com) // @author Jordan DeLong (delong.j@fb.com)
// Changes made by PSD Detector Group: // Changes made by PSD Detector Group:
// Copied: Line 34 constexpr std::size_t hardware_destructive_interference_size = 128; from folly/lang/Align.h // Copied: Line 34 constexpr std::size_t hardware_destructive_interference_size
// Changed extension to .hpp // = 128; from folly/lang/Align.h Changed extension to .hpp Changed namespace to
// Changed namespace to aare // aare
#pragma once #pragma once
@ -45,7 +45,6 @@ template <class T> struct ProducerConsumerQueue {
ProducerConsumerQueue(const ProducerConsumerQueue &) = delete; ProducerConsumerQueue(const ProducerConsumerQueue &) = delete;
ProducerConsumerQueue &operator=(const ProducerConsumerQueue &) = delete; ProducerConsumerQueue &operator=(const ProducerConsumerQueue &) = delete;
ProducerConsumerQueue(ProducerConsumerQueue &&other) { ProducerConsumerQueue(ProducerConsumerQueue &&other) {
size_ = other.size_; size_ = other.size_;
records_ = other.records_; records_ = other.records_;
@ -62,7 +61,6 @@ template <class T> struct ProducerConsumerQueue {
return *this; return *this;
} }
ProducerConsumerQueue() : ProducerConsumerQueue(2){}; ProducerConsumerQueue() : ProducerConsumerQueue(2){};
// size must be >= 2. // size must be >= 2.
// //
@ -70,7 +68,9 @@ template <class T> struct ProducerConsumerQueue {
// given time is actually (size-1), so if you start with an empty queue, // given time is actually (size-1), so if you start with an empty queue,
// isFull() will return true after size-1 insertions. // isFull() will return true after size-1 insertions.
explicit ProducerConsumerQueue(uint32_t size) explicit ProducerConsumerQueue(uint32_t size)
: size_(size), records_(static_cast<T *>(std::malloc(sizeof(T) * size))), readIndex_(0), writeIndex_(0) { : size_(size),
records_(static_cast<T *>(std::malloc(sizeof(T) * size))),
readIndex_(0), writeIndex_(0) {
assert(size >= 2); assert(size >= 2);
if (!records_) { if (!records_) {
throw std::bad_alloc(); throw std::bad_alloc();
@ -154,7 +154,8 @@ template <class T> struct ProducerConsumerQueue {
} }
bool isEmpty() const { 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 { bool isFull() const {
@ -175,7 +176,8 @@ template <class T> struct ProducerConsumerQueue {
// be removing items concurrently). // be removing items concurrently).
// * It is undefined to call this from any other thread. // * It is undefined to call this from any other thread.
size_t sizeGuess() const { 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) { if (ret < 0) {
ret += size_; ret += size_;
} }

View File

@ -1,11 +1,10 @@
#pragma once #pragma once
#include "aare/FileInterface.hpp" #include "aare/FileInterface.hpp"
#include "aare/RawMasterFile.hpp"
#include "aare/Frame.hpp" #include "aare/Frame.hpp"
#include "aare/NDArray.hpp" //for pixel map #include "aare/NDArray.hpp" //for pixel map
#include "aare/RawMasterFile.hpp"
#include "aare/RawSubFile.hpp" #include "aare/RawSubFile.hpp"
#include <optional> #include <optional>
namespace aare { namespace aare {
@ -55,8 +54,8 @@ class RawFile : public FileInterface {
// 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, 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 frame_number(size_t frame_index) override;
size_t bytes_per_frame() override; size_t bytes_per_frame() override;
@ -73,20 +72,17 @@ class RawFile : public FileInterface {
RawMasterFile master() const; RawMasterFile master() const;
DetectorType detector_type() const override; DetectorType detector_type() const override;
private: private:
/** /**
* @brief read the frame at the given frame index into the image buffer * @brief read the frame at the given frame index into the image buffer
* @param frame_number frame number to read * @param frame_number frame number to read
* @param image_buf buffer to store the frame * @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 * @brief get the frame at the given frame index
@ -95,8 +91,6 @@ class RawFile : public FileInterface {
*/ */
Frame get_frame(size_t frame_index); Frame get_frame(size_t frame_index);
/** /**
* @brief read the header of the file * @brief read the header of the file
* @param fname path to the data subfile * @param fname path to the data subfile
@ -108,5 +102,4 @@ class RawFile : public FileInterface {
void find_geometry(); void find_geometry();
}; };
} // namespace aare } // namespace aare

View File

@ -79,7 +79,6 @@ class RawMasterFile {
std::optional<ROI> m_roi; std::optional<ROI> m_roi;
public: public:
RawMasterFile(const std::filesystem::path &fpath); RawMasterFile(const std::filesystem::path &fpath);
@ -107,7 +106,6 @@ class RawMasterFile {
std::optional<size_t> number_of_rows() const; std::optional<size_t> number_of_rows() const;
std::optional<uint8_t> quad() const; std::optional<uint8_t> quad() const;
std::optional<ROI> roi() const; std::optional<ROI> roi() const;
ScanParameters scan_parameters() const; ScanParameters scan_parameters() const;

View File

@ -10,8 +10,9 @@
namespace aare { namespace aare {
/** /**
* @brief Class to read a singe subfile written in .raw format. Used from RawFile to read * @brief Class to read a singe subfile written in .raw format. Used from
* the entire detector. Can be used directly to read part of the image. * RawFile to read the entire detector. Can be used directly to read part of the
* image.
*/ */
class RawSubFile { class RawSubFile {
protected: protected:
@ -20,22 +21,23 @@ class RawSubFile {
size_t m_bitdepth; size_t m_bitdepth;
std::filesystem::path m_path; //!< path to the subfile std::filesystem::path m_path; //!< path to the subfile
std::string m_base_name; //!< base name used for formatting file names 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_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_total_frames{}; //!< total number of frames in the series of files
size_t m_rows{}; size_t m_rows{};
size_t m_cols{}; size_t m_cols{};
size_t m_bytes_per_frame{}; size_t m_bytes_per_frame{};
int m_module_index{}; int m_module_index{};
size_t m_current_file_index{}; //!< The index of the open 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_current_frame_index{}; //!< The index of the current frame (with
std::vector<size_t> m_last_frame_in_file{}; //!< Used for seeking to the correct file //!< reference to all files)
std::vector<size_t>
m_last_frame_in_file{}; //!< Used for seeking to the correct file
uint32_t m_pos_row{}; uint32_t m_pos_row{};
uint32_t m_pos_col{}; uint32_t m_pos_col{};
std::optional<NDArray<ssize_t, 2>> m_pixel_map; std::optional<NDArray<ssize_t, 2>> m_pixel_map;
public: public:
@ -49,12 +51,14 @@ class RawSubFile {
* @throws std::invalid_argument if the detector,type pair is not supported * @throws std::invalid_argument if the detector,type pair is not supported
*/ */
RawSubFile(const std::filesystem::path &fname, DetectorType detector, 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; ~RawSubFile() = default;
/** /**
* @brief Seek to the given frame number * @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 * @param frame_index frame position in file to seek to
* @throws std::runtime_error if the frame number is out of range * @throws std::runtime_error if the frame number is out of range
*/ */
@ -62,7 +66,8 @@ class RawSubFile {
size_t tell(); size_t tell();
void read_into(std::byte *image_buf, DetectorHeader *header = nullptr); 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 get_part(std::byte *buffer, size_t frame_index);
void read_header(DetectorHeader *header); void read_header(DetectorHeader *header);
@ -79,14 +84,12 @@ class RawSubFile {
size_t frames_in_file() const { return m_total_frames; } size_t frames_in_file() const { return m_total_frames; }
private: private:
template <typename T> template <typename T> void read_with_map(std::byte *image_buf);
void read_with_map(std::byte *image_buf);
void parse_fname(const std::filesystem::path &fname); void parse_fname(const std::filesystem::path &fname);
void scan_files(); void scan_files();
void open_file(size_t file_index); void open_file(size_t file_index);
std::filesystem::path fpath(size_t file_index) const; std::filesystem::path fpath(size_t file_index) const;
}; };
} // namespace aare } // namespace aare

View File

@ -38,8 +38,10 @@ template <typename T> class VarClusterFinder {
bool use_noise_map = false; bool use_noise_map = false;
int peripheralThresholdFactor_ = 5; int peripheralThresholdFactor_ = 5;
int current_label; int current_label;
const std::array<int, 4> di{{0, -1, -1, -1}}; // row ### 8-neighbour by scaning from left to right const std::array<int, 4> di{
const std::array<int, 4> dj{{-1, -1, 0, 1}}; // col ### 8-neighbour by scaning from top to bottom {0, -1, -1, -1}}; // row ### 8-neighbour by scaning from left to right
const std::array<int, 4> dj{
{-1, -1, 0, 1}}; // col ### 8-neighbour by scaning from top to bottom
const std::array<int, 8> di_{{0, 0, -1, 1, -1, 1, -1, 1}}; // row const std::array<int, 8> di_{{0, 0, -1, 1, -1, 1, -1, 1}}; // row
const std::array<int, 8> dj_{{-1, 1, 0, 0, 1, -1, -1, 1}}; // col const std::array<int, 8> dj_{{-1, 1, 0, 0, 1, -1, -1, 1}}; // col
std::map<int, int> child; // heirachy: key: child; val: parent std::map<int, int> child; // heirachy: key: child; val: parent
@ -50,7 +52,8 @@ template <typename T> class VarClusterFinder {
public: public:
VarClusterFinder(Shape<2> shape, T threshold) 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); hits.reserve(2000);
} }
@ -60,7 +63,9 @@ template <typename T> class VarClusterFinder {
noiseMap = noise_map; noiseMap = noise_map;
use_noise_map = true; use_noise_map = true;
} }
void set_peripheralThresholdFactor(int factor) { peripheralThresholdFactor_ = factor; } void set_peripheralThresholdFactor(int factor) {
peripheralThresholdFactor_ = factor;
}
void find_clusters(NDView<T, 2> img); void find_clusters(NDView<T, 2> img);
void find_clusters_X(NDView<T, 2> img); void find_clusters_X(NDView<T, 2> img);
void rec_FillHit(int clusterIndex, int i, int j); void rec_FillHit(int clusterIndex, int i, int j);
@ -144,7 +149,8 @@ template <typename T> int VarClusterFinder<T>::check_neighbours(int i, int j) {
} }
} }
template <typename T> void VarClusterFinder<T>::find_clusters(NDView<T, 2> img) { template <typename T>
void VarClusterFinder<T>::find_clusters(NDView<T, 2> img) {
original_ = img; original_ = img;
labeled_ = 0; labeled_ = 0;
peripheral_labeled_ = 0; peripheral_labeled_ = 0;
@ -156,7 +162,8 @@ template <typename T> void VarClusterFinder<T>::find_clusters(NDView<T, 2> img)
store_clusters(); store_clusters();
} }
template <typename T> void VarClusterFinder<T>::find_clusters_X(NDView<T, 2> img) { template <typename T>
void VarClusterFinder<T>::find_clusters_X(NDView<T, 2> img) {
original_ = img; original_ = img;
int clusterIndex = 0; int clusterIndex = 0;
for (int i = 0; i < shape_[0]; ++i) { for (int i = 0; i < shape_[0]; ++i) {
@ -175,7 +182,8 @@ template <typename T> void VarClusterFinder<T>::find_clusters_X(NDView<T, 2> img
h_size.clear(); h_size.clear();
} }
template <typename T> void VarClusterFinder<T>::rec_FillHit(int clusterIndex, int i, int j) { template <typename T>
void VarClusterFinder<T>::rec_FillHit(int clusterIndex, int i, int j) {
// printf("original_(%d, %d)=%f\n", i, j, original_(i,j)); // printf("original_(%d, %d)=%f\n", i, j, original_(i,j));
// printf("h_size[%d].size=%d\n", clusterIndex, h_size[clusterIndex].size); // printf("h_size[%d].size=%d\n", clusterIndex, h_size[clusterIndex].size);
if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE) { if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE) {
@ -203,11 +211,15 @@ template <typename T> void VarClusterFinder<T>::rec_FillHit(int clusterIndex, in
} else { } else {
// if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE){ // if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE){
// h_size[clusterIndex].size += 1; // h_size[clusterIndex].size += 1;
// h_size[clusterIndex].rows[h_size[clusterIndex].size] = row; // h_size[clusterIndex].rows[h_size[clusterIndex].size] =
// h_size[clusterIndex].cols[h_size[clusterIndex].size] = col; // row; h_size[clusterIndex].cols[h_size[clusterIndex].size]
// h_size[clusterIndex].enes[h_size[clusterIndex].size] = original_(row, col); // = col;
// h_size[clusterIndex].enes[h_size[clusterIndex].size] =
// original_(row, col);
// }// ? weather to include peripheral pixels // }// ? 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 <typename T> void VarClusterFinder<T>::store_clusters() {
for (int i = 0; i < shape_[0]; ++i) { for (int i = 0; i < shape_[0]; ++i) {
for (int j = 0; j < shape_[1]; ++j) { for (int j = 0; j < shape_[1]; ++j) {
if (labeled_(i, j) != 0 || false if (labeled_(i, j) != 0 || false
// (i-1 >= 0 and labeled_(i-1, j) != 0) or // another circle of peripheral pixels // (i-1 >= 0 and labeled_(i-1, j) != 0) or // another circle of
// (j-1 >= 0 and labeled_(i, j-1) != 0) or // peripheral pixels (j-1 >= 0 and labeled_(i, j-1) != 0) or
// (i+1 < shape_[0] and labeled_(i+1, j) != 0) or // (i+1 < shape_[0] and labeled_(i+1, j) != 0) or
// (j+1 < shape_[1] and labeled_(i, j+1) != 0) // (j+1 < shape_[1] and labeled_(i, j+1) != 0)
) { ) {

View File

@ -1,9 +1,9 @@
#pragma once #pragma once
#include <aare/NDArray.hpp>
#include <algorithm> #include <algorithm>
#include <array> #include <array>
#include <vector> #include <vector>
#include <aare/NDArray.hpp>
namespace aare { namespace aare {
/** /**
@ -27,13 +27,11 @@ size_t last_smaller(const T* first, const T* last, T val) {
return std::distance(first, last - 1); return std::distance(first, last - 1);
} }
template <typename T> template <typename T> size_t last_smaller(const NDArray<T, 1> &arr, T val) {
size_t last_smaller(const NDArray<T, 1>& arr, T val) {
return last_smaller(arr.begin(), arr.end(), val); return last_smaller(arr.begin(), arr.end(), val);
} }
template <typename T> template <typename T> size_t last_smaller(const std::vector<T> &vec, T val) {
size_t last_smaller(const std::vector<T>& vec, T val) {
return last_smaller(vec.data(), vec.data() + vec.size(), val); return last_smaller(vec.data(), vec.data() + vec.size(), val);
} }
@ -57,19 +55,18 @@ size_t first_larger(const T* first, const T* last, T val) {
return std::distance(first, last - 1); return std::distance(first, last - 1);
} }
template <typename T> template <typename T> size_t first_larger(const NDArray<T, 1> &arr, T val) {
size_t first_larger(const NDArray<T, 1>& arr, T val) {
return first_larger(arr.begin(), arr.end(), val); return first_larger(arr.begin(), arr.end(), val);
} }
template <typename T> template <typename T> size_t first_larger(const std::vector<T> &vec, T val) {
size_t first_larger(const std::vector<T>& vec, T val) {
return first_larger(vec.data(), vec.data() + vec.size(), val); return first_larger(vec.data(), vec.data() + vec.size(), val);
} }
/** /**
* @brief Index of the nearest element to 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 first iterator to the first element
* @param last iterator to the last element * @param last iterator to the last element
* @param val value to compare * @param val value to compare
@ -77,20 +74,17 @@ size_t first_larger(const std::vector<T>& vec, T val) {
*/ */
template <typename T> template <typename T>
size_t nearest_index(const T *first, const T *last, T val) { size_t nearest_index(const T *first, const T *last, T val) {
auto iter = std::min_element(first, last, auto iter = std::min_element(first, last, [val](T a, T b) {
[val](T a, T b) {
return std::abs(a - val) < std::abs(b - val); return std::abs(a - val) < std::abs(b - val);
}); });
return std::distance(first, iter); return std::distance(first, iter);
} }
template <typename T> template <typename T> size_t nearest_index(const NDArray<T, 1> &arr, T val) {
size_t nearest_index(const NDArray<T, 1>& arr, T val) {
return nearest_index(arr.begin(), arr.end(), val); return nearest_index(arr.begin(), arr.end(), val);
} }
template <typename T> template <typename T> size_t nearest_index(const std::vector<T> &vec, T val) {
size_t nearest_index(const std::vector<T>& vec, T val) {
return nearest_index(vec.data(), vec.data() + vec.size(), val); return nearest_index(vec.data(), vec.data() + vec.size(), val);
} }
@ -99,14 +93,12 @@ size_t nearest_index(const std::array<T,N>& arr, T val) {
return nearest_index(arr.data(), arr.data() + arr.size(), val); return nearest_index(arr.data(), arr.data() + arr.size(), val);
} }
template <typename T> template <typename T> std::vector<T> cumsum(const std::vector<T> &vec) {
std::vector<T> cumsum(const std::vector<T>& vec) {
std::vector<T> result(vec.size()); std::vector<T> result(vec.size());
std::partial_sum(vec.begin(), vec.end(), result.begin()); std::partial_sum(vec.begin(), vec.end(), result.begin());
return result; return result;
} }
template <typename Container> bool all_equal(const Container &c) { template <typename Container> bool all_equal(const Container &c) {
if (!c.empty() && if (!c.empty() &&
std::all_of(begin(c), end(c), std::all_of(begin(c), end(c),
@ -117,6 +109,4 @@ template <typename Container> bool all_equal(const Container &c) {
return false; return false;
} }
} // namespace aare } // namespace aare

View File

@ -1,26 +1,27 @@
#pragma once #pragma once
#include <aare/NDView.hpp>
#include <cstdint> #include <cstdint>
#include <vector> #include <vector>
#include <aare/NDView.hpp>
namespace aare { namespace aare {
uint16_t adc_sar_05_decode64to16(uint64_t input); uint16_t adc_sar_05_decode64to16(uint64_t input);
uint16_t adc_sar_04_decode64to16(uint64_t input); uint16_t adc_sar_04_decode64to16(uint64_t input);
void adc_sar_05_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output); void adc_sar_05_decode64to16(NDView<uint64_t, 2> input,
void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output); NDView<uint16_t, 2> output);
void adc_sar_04_decode64to16(NDView<uint64_t, 2> input,
NDView<uint16_t, 2> output);
/** /**
* @brief Apply custom weights to a 16-bit input value. Will sum up weights[i]**i * @brief Apply custom weights to a 16-bit input value. Will sum up
* for each bit i that is set in the input value. * weights[i]**i for each bit i that is set in the input value.
* @throws std::out_of_range if weights.size() < 16 * @throws std::out_of_range if weights.size() < 16
* @param input 16-bit input value * @param input 16-bit input value
* @param weights vector of weights, size must be less than or equal to 16 * @param weights vector of weights, size must be less than or equal to 16
*/ */
double apply_custom_weights(uint16_t input, const NDView<double, 1> weights); double apply_custom_weights(uint16_t input, const NDView<double, 1> weights);
void apply_custom_weights(NDView<uint16_t, 1> input, NDView<double, 1> output, const NDView<double, 1> weights); void apply_custom_weights(NDView<uint16_t, 1> input, NDView<double, 1> output,
const NDView<double, 1> weights);
} // namespace aare } // namespace aare

View File

@ -25,28 +25,24 @@
std::string(__FILE__) + std::string(":") + std::to_string(__LINE__) + \ std::string(__FILE__) + std::string(":") + std::to_string(__LINE__) + \
":" + std::string(__func__) + ":" ":" + std::string(__func__) + ":"
#ifdef AARE_CUSTOM_ASSERT #ifdef AARE_CUSTOM_ASSERT
#define AARE_ASSERT(expr) \ #define AARE_ASSERT(expr) \
if (expr)\ if (expr) { \
{}\ } else \
else\
aare::assert_failed(LOCATION + " Assertion failed: " + #expr + "\n"); aare::assert_failed(LOCATION + " Assertion failed: " + #expr + "\n");
#else #else
#define AARE_ASSERT(cond) \ #define AARE_ASSERT(cond) \
do { (void)sizeof(cond); } while(0) do { \
(void)sizeof(cond); \
} while (0)
#endif #endif
namespace aare { namespace aare {
inline constexpr size_t bits_per_byte = 8; inline constexpr size_t bits_per_byte = 8;
void assert_failed(const std::string &msg); void assert_failed(const std::string &msg);
class DynamicCluster { class DynamicCluster {
public: public:
int cluster_sizeX; int cluster_sizeX;
@ -184,9 +180,9 @@ template <typename T> struct t_xy {
}; };
using xy = t_xy<uint32_t>; using xy = t_xy<uint32_t>;
/** /**
* @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_x{};
@ -198,8 +194,8 @@ struct ModuleGeometry{
}; };
/** /**
* @brief Class to hold the geometry of a detector. Number of modules, their size and where pixel 0 * @brief Class to hold the geometry of a detector. Number of modules, their
* for each module is located * size and where pixel 0 for each module is located
*/ */
struct DetectorGeometry { struct DetectorGeometry {
int modules_x{}; int modules_x{};
@ -226,7 +222,6 @@ struct ROI{
} }
}; };
using dynamic_shape = std::vector<ssize_t>; using dynamic_shape = std::vector<ssize_t>;
@ -288,7 +283,8 @@ enum class DetectorType {
Gotthard2, Gotthard2,
Xilinx_ChipTestBoard, Xilinx_ChipTestBoard,
//Additional detectors used for defining processing. Variants of the standard ones. // Additional detectors used for defining processing. Variants of the
// standard ones.
Moench03 = 100, Moench03 = 100,
Moench03_old, Moench03_old,
Unknown Unknown

View File

@ -1,6 +1,6 @@
#pragma once #pragma once
#include "aare/defs.hpp"
#include "aare/RawMasterFile.hpp" //ROI refactor away #include "aare/RawMasterFile.hpp" //ROI refactor away
#include "aare/defs.hpp"
namespace aare { namespace aare {
/** /**
@ -12,5 +12,4 @@ namespace aare{
*/ */
DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, ROI roi); DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, ROI roi);
} // namespace aare } // namespace aare

View File

@ -1,7 +1,6 @@
#pragma once #pragma once
/*Utility to log to console*/ /*Utility to log to console*/
#include <iostream> #include <iostream>
#include <sstream> #include <sstream>
#include <sys/time.h> #include <sys/time.h>
@ -27,7 +26,6 @@ namespace aare {
#define RESET "\x1b[0m" #define RESET "\x1b[0m"
#define BOLD "\x1b[1m" #define BOLD "\x1b[1m"
enum TLogLevel { enum TLogLevel {
logERROR, logERROR,
logWARNING, logWARNING,
@ -37,7 +35,8 @@ enum TLogLevel {
logINFOCYAN, logINFOCYAN,
logINFOMAGENTA, logINFOMAGENTA,
logINFO, logINFO,
logDEBUG, // constructors, destructors etc. should still give too much output logDEBUG, // constructors, destructors etc. should still give too much
// output
logDEBUG1, logDEBUG1,
logDEBUG2, logDEBUG2,
logDEBUG3, logDEBUG3,
@ -47,7 +46,9 @@ enum TLogLevel {
// Compiler should optimize away anything below this value // Compiler should optimize away anything below this value
#ifndef AARE_LOG_LEVEL #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 #endif
#define __AT__ \ #define __AT__ \
@ -72,7 +73,8 @@ class Logger {
std::clog << os.str() << std::flush; // Single write 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; static TLogLevel reportingLevel = logDEBUG5;
return reportingLevel; return reportingLevel;
} }

View File

@ -1,6 +1,6 @@
#include <thread> #include <thread>
#include <vector>
#include <utility> #include <utility>
#include <vector>
namespace aare { namespace aare {

View File

@ -1,10 +1,10 @@
#include "aare/Cluster.hpp" #include "aare/Cluster.hpp"
#include <cstdint> #include <cstdint>
#include <fmt/format.h>
#include <filesystem> #include <filesystem>
#include <pybind11/pybind11.h> #include <fmt/format.h>
#include <pybind11/numpy.h> #include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h> #include <pybind11/stl.h>
#include <pybind11/stl_bind.h> #include <pybind11/stl_bind.h>

View File

@ -21,11 +21,9 @@ using namespace aare;
#pragma GCC diagnostic push #pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter" #pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY, template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t> typename CoordType = uint16_t>
void define_ClusterCollector(py::module &m, void define_ClusterCollector(py::module &m, const std::string &typestr) {
const std::string &typestr) {
auto class_name = fmt::format("ClusterCollector_{}", typestr); auto class_name = fmt::format("ClusterCollector_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>; using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;

View File

@ -21,8 +21,7 @@ using namespace ::aare;
template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY, template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
typename CoordType = uint16_t> typename CoordType = uint16_t>
void define_ClusterFile(py::module &m, void define_ClusterFile(py::module &m, const std::string &typestr) {
const std::string &typestr) {
using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>; using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;

View File

@ -21,15 +21,9 @@ using namespace aare;
#pragma GCC diagnostic push #pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter" #pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY, template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t> typename CoordType = uint16_t>
void define_ClusterFileSink(py::module &m, void define_ClusterFileSink(py::module &m, const std::string &typestr) {
const std::string &typestr) {
auto class_name = fmt::format("ClusterFileSink_{}", typestr); auto class_name = fmt::format("ClusterFileSink_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>; using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
@ -40,5 +34,4 @@ void define_ClusterFileSink(py::module &m,
.def("stop", &ClusterFileSink<ClusterType>::stop); .def("stop", &ClusterFileSink<ClusterType>::stop);
} }
#pragma GCC diagnostic pop #pragma GCC diagnostic pop

View File

@ -23,8 +23,7 @@ using namespace aare;
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY, template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t> typename CoordType = uint16_t>
void define_ClusterFinderMT(py::module &m, void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
const std::string &typestr) {
auto class_name = fmt::format("ClusterFinderMT_{}", typestr); auto class_name = fmt::format("ClusterFinderMT_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>; using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
@ -49,7 +48,9 @@ void define_ClusterFinderMT(py::module &m,
return; return;
}, },
py::arg(), py::arg("frame_number") = 0) py::arg(), py::arg("frame_number") = 0)
.def_property_readonly("cluster_size", [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self){ .def_property_readonly(
"cluster_size",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self) {
return py::make_tuple(ClusterSizeX, ClusterSizeY); return py::make_tuple(ClusterSizeX, ClusterSizeY);
}) })
.def("clear_pedestal", .def("clear_pedestal",
@ -77,5 +78,4 @@ void define_ClusterFinderMT(py::module &m,
py::arg("thread_index") = 0); py::arg("thread_index") = 0);
} }
#pragma GCC diagnostic pop #pragma GCC diagnostic pop

View File

@ -44,7 +44,8 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
auto *vec = new std::vector<Type>(self.sum()); auto *vec = new std::vector<Type>(self.sum());
return return_vector(vec); return return_vector(vec);
}) })
.def("sum_2x2", [](ClusterVector<ClusterType> &self){ .def("sum_2x2",
[](ClusterVector<ClusterType> &self) {
auto *vec = new std::vector<Type>(self.sum_2x2()); auto *vec = new std::vector<Type>(self.sum_2x2());
return return_vector(vec); return return_vector(vec);
}) })

View File

@ -6,8 +6,8 @@
#include "aare/RawMasterFile.hpp" #include "aare/RawMasterFile.hpp"
#include "aare/RawSubFile.hpp" #include "aare/RawSubFile.hpp"
#include "aare/defs.hpp"
#include "aare/decode.hpp" #include "aare/decode.hpp"
#include "aare/defs.hpp"
// #include "aare/fClusterFileV2.hpp" // #include "aare/fClusterFileV2.hpp"
#include "np_helper.hpp" #include "np_helper.hpp"
@ -27,61 +27,69 @@ using namespace ::aare;
void define_ctb_raw_file_io_bindings(py::module &m) { void define_ctb_raw_file_io_bindings(py::module &m) {
m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) { m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) {
if (input.ndim() != 2) { if (input.ndim() != 2) {
throw std::runtime_error("Only 2D arrays are supported at this moment"); throw std::runtime_error(
"Only 2D arrays are supported at this moment");
} }
// Create a 2D output array with the same shape as the input // Create a 2D output array with the same shape as the input
std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<ssize_t>(bits_per_byte)}; std::vector<ssize_t> shape{input.shape(0),
input.shape(1) /
static_cast<ssize_t>(bits_per_byte)};
py::array_t<uint16_t> output(shape); py::array_t<uint16_t> output(shape);
// Create a view of the input and output arrays // Create a view of the input and output arrays
NDView<uint64_t, 2> input_view(reinterpret_cast<uint64_t*>(input.mutable_data()), {output.shape(0), output.shape(1)}); NDView<uint64_t, 2> input_view(
NDView<uint16_t, 2> output_view(output.mutable_data(), {output.shape(0), output.shape(1)}); reinterpret_cast<uint64_t *>(input.mutable_data()),
{output.shape(0), output.shape(1)});
NDView<uint16_t, 2> 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; return output;
}); });
m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) { m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
if (input.ndim() != 2) { if (input.ndim() != 2) {
throw std::runtime_error("Only 2D arrays are supported at this moment"); throw std::runtime_error(
"Only 2D arrays are supported at this moment");
} }
// Create a 2D output array with the same shape as the input // Create a 2D output array with the same shape as the input
std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<ssize_t>(bits_per_byte)}; std::vector<ssize_t> shape{input.shape(0),
input.shape(1) /
static_cast<ssize_t>(bits_per_byte)};
py::array_t<uint16_t> output(shape); py::array_t<uint16_t> output(shape);
// Create a view of the input and output arrays // Create a view of the input and output arrays
NDView<uint64_t, 2> input_view(reinterpret_cast<uint64_t*>(input.mutable_data()), {output.shape(0), output.shape(1)}); NDView<uint64_t, 2> input_view(
NDView<uint16_t, 2> output_view(output.mutable_data(), {output.shape(0), output.shape(1)}); reinterpret_cast<uint64_t *>(input.mutable_data()),
{output.shape(0), output.shape(1)});
NDView<uint16_t, 2> output_view(output.mutable_data(),
{output.shape(0), output.shape(1)});
adc_sar_04_decode64to16(input_view, output_view); adc_sar_04_decode64to16(input_view, output_view);
return output; return output;
}); });
m.def( m.def("apply_custom_weights",
"apply_custom_weights", [](py::array_t<uint16_t, py::array::c_style | py::array::forcecast>
[](py::array_t<uint16_t, py::array::c_style | py::array::forcecast> &input, &input,
py::array_t<double, py::array::c_style | py::array::forcecast> py::array_t<double, py::array::c_style | py::array::forcecast>
&weights) { &weights) {
// Create new array with same shape as the input array
// (uninitialized values)
// Create new array with same shape as the input array (uninitialized values)
py::buffer_info buf = input.request(); py::buffer_info buf = input.request();
py::array_t<double> output(buf.shape); py::array_t<double> output(buf.shape);
// Use NDViews to call into the C++ library // Use NDViews to call into the C++ library
auto weights_view = make_view_1d(weights); auto weights_view = make_view_1d(weights);
NDView<uint16_t, 1> input_view(input.mutable_data(), {input.size()}); NDView<uint16_t, 1> input_view(input.mutable_data(),
NDView<double, 1> output_view(output.mutable_data(), {output.size()}); {input.size()});
NDView<double, 1> output_view(output.mutable_data(),
{output.size()});
apply_custom_weights(input_view, output_view, weights_view); apply_custom_weights(input_view, output_view, weights_view);
return output; return output;
@ -103,7 +111,8 @@ py::class_<CtbRawFile>(m, "CtbRawFile")
// always read bytes // always read bytes
image = py::array_t<uint8_t>(shape); image = py::array_t<uint8_t>(shape);
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()), self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()),
header.mutable_data()); header.mutable_data());
return py::make_tuple(header, image); return py::make_tuple(header, image);
@ -116,5 +125,4 @@ py::class_<CtbRawFile>(m, "CtbRawFile")
&CtbRawFile::image_size_in_bytes) &CtbRawFile::image_size_in_bytes)
.def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file); .def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
} }

View File

@ -25,9 +25,6 @@
namespace py = pybind11; namespace py = pybind11;
using namespace ::aare; using namespace ::aare;
// Disable warnings for unused parameters, as we ignore some // Disable warnings for unused parameters, as we ignore some
// in the __exit__ method // in the __exit__ method
#pragma GCC diagnostic push #pragma GCC diagnostic push
@ -35,7 +32,6 @@ using namespace ::aare;
void define_file_io_bindings(py::module &m) { void define_file_io_bindings(py::module &m) {
py::enum_<DetectorType>(m, "DetectorType") py::enum_<DetectorType>(m, "DetectorType")
.value("Jungfrau", DetectorType::Jungfrau) .value("Jungfrau", DetectorType::Jungfrau)
.value("Eiger", DetectorType::Eiger) .value("Eiger", DetectorType::Eiger)
@ -46,13 +42,10 @@ void define_file_io_bindings(py::module &m) {
.value("ChipTestBoard", DetectorType::ChipTestBoard) .value("ChipTestBoard", DetectorType::ChipTestBoard)
.value("Unknown", DetectorType::Unknown); .value("Unknown", DetectorType::Unknown);
PYBIND11_NUMPY_DTYPE(DetectorHeader, frameNumber, expLength, packetNumber, PYBIND11_NUMPY_DTYPE(DetectorHeader, frameNumber, expLength, packetNumber,
bunchId, timestamp, modId, row, column, reserved, bunchId, timestamp, modId, row, column, reserved,
debug, roundRNumber, detType, version, packetMask); debug, roundRNumber, detType, version, packetMask);
py::class_<File>(m, "File") py::class_<File>(m, "File")
.def(py::init([](const std::filesystem::path &fname) { .def(py::init([](const std::filesystem::path &fname) {
return File(fname, "r", {}); return File(fname, "r", {});
@ -117,9 +110,11 @@ void define_file_io_bindings(py::module &m) {
reinterpret_cast<std::byte *>(image.mutable_data())); reinterpret_cast<std::byte *>(image.mutable_data()));
return image; return image;
}) })
.def("read_n", [](File &self, size_t n_frames) { .def("read_n",
[](File &self, size_t n_frames) {
// adjust for actual frames left in the file // adjust for actual frames left in the file
n_frames = std::min(n_frames, self.total_frames()-self.tell()); n_frames =
std::min(n_frames, self.total_frames() - self.tell());
if (n_frames == 0) { if (n_frames == 0) {
throw std::runtime_error("No frames left in file"); throw std::runtime_error("No frames left in file");
} }
@ -134,21 +129,20 @@ void define_file_io_bindings(py::module &m) {
} else if (item_size == 4) { } else if (item_size == 4) {
image = py::array_t<uint32_t>(shape); image = py::array_t<uint32_t>(shape);
} }
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()), self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()),
n_frames); n_frames);
return image; return image;
}) })
.def("__enter__", [](File &self) { return &self; }) .def("__enter__", [](File &self) { return &self; })
.def("__exit__", .def("__exit__",
[](File &self, [](File &self, const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::object> &exc_value, const std::optional<pybind11::object> &exc_value,
const std::optional<pybind11::object> &traceback) { const std::optional<pybind11::object> &traceback) {
// self.close(); // self.close();
}) })
.def("__iter__", [](File &self) { return &self; }) .def("__iter__", [](File &self) { return &self; })
.def("__next__", [](File &self) { .def("__next__", [](File &self) {
try { try {
const uint8_t item_size = self.bytes_per_pixel(); const uint8_t item_size = self.bytes_per_pixel();
py::array image; py::array image;
@ -171,7 +165,6 @@ void define_file_io_bindings(py::module &m) {
} }
}); });
py::class_<FileConfig>(m, "FileConfig") py::class_<FileConfig>(m, "FileConfig")
.def(py::init<>()) .def(py::init<>())
.def_readwrite("rows", &FileConfig::rows) .def_readwrite("rows", &FileConfig::rows)
@ -188,8 +181,6 @@ void define_file_io_bindings(py::module &m) {
return "<FileConfig: " + a.to_string() + ">"; return "<FileConfig: " + a.to_string() + ">";
}); });
py::class_<ScanParameters>(m, "ScanParameters") py::class_<ScanParameters>(m, "ScanParameters")
.def(py::init<const std::string &>()) .def(py::init<const std::string &>())
.def(py::init<const ScanParameters &>()) .def(py::init<const ScanParameters &>())
@ -200,7 +191,6 @@ void define_file_io_bindings(py::module &m) {
.def_property_readonly("stop", &ScanParameters::stop) .def_property_readonly("stop", &ScanParameters::stop)
.def_property_readonly("step", &ScanParameters::step); .def_property_readonly("step", &ScanParameters::step);
py::class_<ROI>(m, "ROI") py::class_<ROI>(m, "ROI")
.def(py::init<>()) .def(py::init<>())
.def(py::init<ssize_t, ssize_t, ssize_t, ssize_t>(), py::arg("xmin"), .def(py::init<ssize_t, ssize_t, ssize_t, ssize_t>(), py::arg("xmin"),
@ -209,23 +199,21 @@ void define_file_io_bindings(py::module &m) {
.def_readwrite("xmax", &ROI::xmax) .def_readwrite("xmax", &ROI::xmax)
.def_readwrite("ymin", &ROI::ymin) .def_readwrite("ymin", &ROI::ymin)
.def_readwrite("ymax", &ROI::ymax) .def_readwrite("ymax", &ROI::ymax)
.def("__str__", [](const ROI& self){ .def("__str__",
return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}", self.xmin, self.xmax, self.ymin, self.ymax); [](const ROI &self) {
return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}",
self.xmin, self.xmax, self.ymin, self.ymax);
}) })
.def("__repr__", [](const ROI& self){ .def("__repr__",
return fmt::format("<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin, self.xmax, self.ymin, self.ymax); [](const ROI &self) {
return fmt::format(
"<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin,
self.xmax, self.ymin, self.ymax);
}) })
.def("__iter__", [](const ROI &self) { .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 #pragma GCC diagnostic pop
// py::class_<ClusterHeader>(m, "ClusterHeader") // py::class_<ClusterHeader>(m, "ClusterHeader")
// .def(py::init<>()) // .def(py::init<>())

View File

@ -9,7 +9,6 @@
namespace py = pybind11; namespace py = pybind11;
using namespace pybind11::literals; using namespace pybind11::literals;
void define_fit_bindings(py::module &m) { void define_fit_bindings(py::module &m) {
// TODO! Evaluate without converting to double // TODO! Evaluate without converting to double
@ -61,7 +60,8 @@ void define_fit_bindings(py::module &m) {
py::array_t<double, py::array::c_style | py::array::forcecast> par) { py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x); auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par); auto par_view = make_view_1d(par);
auto y = new NDArray<double, 1>{aare::func::scurve(x_view, par_view)}; auto y =
new NDArray<double, 1>{aare::func::scurve(x_view, par_view)};
return return_image_data(y); return return_image_data(y);
}, },
R"( R"(
@ -82,7 +82,8 @@ void define_fit_bindings(py::module &m) {
py::array_t<double, py::array::c_style | py::array::forcecast> par) { py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x); auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par); auto par_view = make_view_1d(par);
auto y = new NDArray<double, 1>{aare::func::scurve2(x_view, par_view)}; auto y =
new NDArray<double, 1>{aare::func::scurve2(x_view, par_view)};
return return_image_data(y); return return_image_data(y);
}, },
R"( R"(
@ -139,7 +140,6 @@ n_threads : int, optional
py::array_t<double, py::array::c_style | py::array::forcecast> y, py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::array_t<double, py::array::c_style | py::array::forcecast> y_err, py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
int n_threads) { int n_threads) {
if (y.ndim() == 3) { if (y.ndim() == 3) {
// Allocate memory for the output // Allocate memory for the output
// Need to have pointers to allow python to manage // 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 y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x); auto x_view = make_view_1d(x);
double chi2 = 0; double chi2 = 0;
aare::fit_gaus(x_view, y_view, y_view_err, par->view(), aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2); par_err->view(), chi2);
@ -252,7 +251,6 @@ n_threads : int, optional
"chi2"_a = return_image_data(chi2), "chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2); "Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) { } else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({2}); auto par = new NDArray<double, 1>({2});
auto par_err = new NDArray<double, 1>({2}); auto par_err = new NDArray<double, 1>({2});
@ -339,7 +337,6 @@ n_threads : int, optional
"chi2"_a = return_image_data(chi2), "chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2); "Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) { } else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({2}); auto par = new NDArray<double, 1>({2});
auto par_err = new NDArray<double, 1>({2}); auto par_err = new NDArray<double, 1>({2});
@ -376,7 +373,6 @@ n_threads : int, optional
)", )",
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4); py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
m.def( m.def(
"fit_scurve2", "fit_scurve2",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x, [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
@ -426,7 +422,6 @@ n_threads : int, optional
"chi2"_a = return_image_data(chi2), "chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2); "Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) { } else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({6}); auto par = new NDArray<double, 1>({6});
auto par_err = new NDArray<double, 1>({6}); auto par_err = new NDArray<double, 1>({6});

View File

@ -21,10 +21,7 @@ using namespace ::aare;
auto read_dat_frame(JungfrauDataFile &self) { auto read_dat_frame(JungfrauDataFile &self) {
py::array_t<JungfrauDataHeader> header(1); py::array_t<JungfrauDataHeader> header(1);
py::array_t<uint16_t> image({ py::array_t<uint16_t> image({self.rows(), self.cols()});
self.rows(),
self.cols()
});
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()), self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
header.mutable_data()); header.mutable_data());
@ -40,9 +37,7 @@ auto read_n_dat_frames(JungfrauDataFile &self, size_t n_frames) {
} }
py::array_t<JungfrauDataHeader> header(n_frames); py::array_t<JungfrauDataHeader> header(n_frames);
py::array_t<uint16_t> image({ py::array_t<uint16_t> image({n_frames, self.rows(), self.cols()});
n_frames, self.rows(),
self.cols()});
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()), self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
n_frames, header.mutable_data()); n_frames, header.mutable_data());

View File

@ -3,10 +3,10 @@
// New style file naming // New style file naming
#include "bind_Cluster.hpp" #include "bind_Cluster.hpp"
#include "bind_ClusterCollector.hpp" #include "bind_ClusterCollector.hpp"
#include "bind_ClusterFinder.hpp"
#include "bind_ClusterFinderMT.hpp"
#include "bind_ClusterFile.hpp" #include "bind_ClusterFile.hpp"
#include "bind_ClusterFileSink.hpp" #include "bind_ClusterFileSink.hpp"
#include "bind_ClusterFinder.hpp"
#include "bind_ClusterFinderMT.hpp"
#include "bind_ClusterVector.hpp" #include "bind_ClusterVector.hpp"
// TODO! migrate the other names // TODO! migrate the other names
@ -14,17 +14,17 @@
#include "file.hpp" #include "file.hpp"
#include "fit.hpp" #include "fit.hpp"
#include "interpolation.hpp" #include "interpolation.hpp"
#include "raw_sub_file.hpp"
#include "raw_master_file.hpp"
#include "raw_file.hpp"
#ifdef HDF5_FOUND #ifdef HDF5_FOUND
#include "hdf5_file.hpp" #include "hdf5_file.hpp"
#include "hdf5_master_file.hpp" #include "hdf5_master_file.hpp"
#endif #endif
#include "pixel_map.hpp"
#include "var_cluster.hpp"
#include "pedestal.hpp"
#include "jungfrau_data_file.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 // Pybind stuff
#include <pybind11/pybind11.h> #include <pybind11/pybind11.h>
@ -38,7 +38,8 @@ T - Storage type of the cluster data (int, float, double)
N - Number of rows in the cluster N - Number of rows in the cluster
M - Number of columns in the cluster M - Number of columns in the cluster
U - Type of the pixel data (e.g., uint16_t) 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 DEFINE_CLUSTER_BINDINGS(T, N, M, U, TYPE_CODE) \

View File

@ -9,7 +9,8 @@
namespace py = pybind11; namespace py = pybind11;
template <typename SUM_TYPE> void define_pedestal_bindings(py::module &m, const std::string &name) { template <typename SUM_TYPE>
void define_pedestal_bindings(py::module &m, const std::string &name) {
py::class_<Pedestal<SUM_TYPE>>(m, name.c_str()) py::class_<Pedestal<SUM_TYPE>>(m, name.c_str())
.def(py::init<int, int, int>()) .def(py::init<int, int, int>())
.def(py::init<int, int>()) .def(py::init<int, int>())
@ -19,12 +20,14 @@ template <typename SUM_TYPE> void define_pedestal_bindings(py::module &m, const
*mea = self.mean(); *mea = self.mean();
return return_image_data(mea); return return_image_data(mea);
}) })
.def("variance", [](Pedestal<SUM_TYPE> &self) { .def("variance",
[](Pedestal<SUM_TYPE> &self) {
auto var = new NDArray<SUM_TYPE, 2>{}; auto var = new NDArray<SUM_TYPE, 2>{};
*var = self.variance(); *var = self.variance();
return return_image_data(var); return return_image_data(var);
}) })
.def("std", [](Pedestal<SUM_TYPE> &self) { .def("std",
[](Pedestal<SUM_TYPE> &self) {
auto std = new NDArray<SUM_TYPE, 2>{}; auto std = new NDArray<SUM_TYPE, 2>{};
*std = self.std(); *std = self.std();
return return_image_data(std); return return_image_data(std);
@ -40,13 +43,18 @@ template <typename SUM_TYPE> void define_pedestal_bindings(py::module &m, const
return Pedestal<SUM_TYPE>(pedestal); return Pedestal<SUM_TYPE>(pedestal);
}) })
// TODO! add push for other data types // TODO! add push for other data types
.def("push", [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> &f) { .def("push",
[](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> &f) {
auto v = make_view_2d(f); auto v = make_view_2d(f);
pedestal.push(v); pedestal.push(v);
}) })
.def("push_no_update", [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t, py::array::c_style> &f) { .def(
"push_no_update",
[](Pedestal<SUM_TYPE> &pedestal,
py::array_t<uint16_t, py::array::c_style> &f) {
auto v = make_view_2d(f); auto v = make_view_2d(f);
pedestal.push_no_update(v); pedestal.push_no_update(v);
}, py::arg().noconvert()) },
py::arg().noconvert())
.def("update_mean", &Pedestal<SUM_TYPE>::update_mean); .def("update_mean", &Pedestal<SUM_TYPE>::update_mean);
} }

View File

@ -1,41 +1,46 @@
#include "aare/PixelMap.hpp" #include "aare/PixelMap.hpp"
#include "np_helper.hpp" #include "np_helper.hpp"
#include <cstdint> #include <cstdint>
#include <pybind11/numpy.h> #include <pybind11/numpy.h>
#include <pybind11/pybind11.h> #include <pybind11/pybind11.h>
#include <pybind11/stl.h> #include <pybind11/stl.h>
namespace py = pybind11; namespace py = pybind11;
using namespace ::aare; using namespace ::aare;
void define_pixel_map_bindings(py::module &m) { void define_pixel_map_bindings(py::module &m) {
m.def("GenerateMoench03PixelMap", []() { m.def("GenerateMoench03PixelMap",
[]() {
auto ptr = new NDArray<ssize_t, 2>(GenerateMoench03PixelMap()); auto ptr = new NDArray<ssize_t, 2>(GenerateMoench03PixelMap());
return return_image_data(ptr); return return_image_data(ptr);
}) })
.def("GenerateMoench05PixelMap", []() { .def("GenerateMoench05PixelMap",
[]() {
auto ptr = new NDArray<ssize_t, 2>(GenerateMoench05PixelMap()); auto ptr = new NDArray<ssize_t, 2>(GenerateMoench05PixelMap());
return return_image_data(ptr); return return_image_data(ptr);
}) })
.def("GenerateMoench05PixelMap1g", []() { .def("GenerateMoench05PixelMap1g",
auto ptr = new NDArray<ssize_t,2>(GenerateMoench05PixelMap1g()); []() {
auto ptr =
new NDArray<ssize_t, 2>(GenerateMoench05PixelMap1g());
return return_image_data(ptr); return return_image_data(ptr);
}) })
.def("GenerateMoench05PixelMapOld", []() { .def("GenerateMoench05PixelMapOld",
auto ptr = new NDArray<ssize_t,2>(GenerateMoench05PixelMapOld()); []() {
auto ptr =
new NDArray<ssize_t, 2>(GenerateMoench05PixelMapOld());
return return_image_data(ptr); return return_image_data(ptr);
}) })
.def("GenerateMH02SingleCounterPixelMap", []() { .def("GenerateMH02SingleCounterPixelMap",
auto ptr = new NDArray<ssize_t,2>(GenerateMH02SingleCounterPixelMap()); []() {
auto ptr = new NDArray<ssize_t, 2>(
GenerateMH02SingleCounterPixelMap());
return return_image_data(ptr); return return_image_data(ptr);
}) })
.def("GenerateMH02FourCounterPixelMap", []() { .def("GenerateMH02FourCounterPixelMap", []() {
auto ptr = new NDArray<ssize_t,3>(GenerateMH02FourCounterPixelMap()); auto ptr =
new NDArray<ssize_t, 3>(GenerateMH02FourCounterPixelMap());
return return_image_data(ptr); return return_image_data(ptr);
}); });
} }

View File

@ -64,7 +64,8 @@ void define_raw_file_io_bindings(py::module &m) {
if (self.n_modules() == 1) { if (self.n_modules() == 1) {
header = py::array_t<DetectorHeader>(n_frames); header = py::array_t<DetectorHeader>(n_frames);
} else { } else {
header = py::array_t<DetectorHeader>({self.n_modules(), n_frames}); header = py::array_t<DetectorHeader>(
{self.n_modules(), n_frames});
} }
// py::array_t<DetectorHeader> header({self.n_mod(), n_frames}); // py::array_t<DetectorHeader> header({self.n_mod(), n_frames});

View File

@ -57,7 +57,8 @@ void define_raw_master_file_bindings(py::module &m) {
.def_property_readonly("total_frames_expected", .def_property_readonly("total_frames_expected",
&RawMasterFile::total_frames_expected) &RawMasterFile::total_frames_expected)
.def_property_readonly("geometry", &RawMasterFile::geometry) .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 Number of analog samples
Returns Returns

View File

@ -43,11 +43,9 @@ auto read_frame_from_RawSubFile(RawSubFile &self) {
auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) { auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) {
py::array_t<DetectorHeader> header(n_frames); py::array_t<DetectorHeader> header(n_frames);
const uint8_t item_size = self.bytes_per_pixel(); const uint8_t item_size = self.bytes_per_pixel();
std::vector<ssize_t> shape{ std::vector<ssize_t> shape{static_cast<ssize_t>(n_frames),
static_cast<ssize_t>(n_frames),
static_cast<ssize_t>(self.rows()), static_cast<ssize_t>(self.rows()),
static_cast<ssize_t>(self.cols()) static_cast<ssize_t>(self.cols())};
};
py::array image; py::array image;
if (item_size == 1) { if (item_size == 1) {
@ -57,13 +55,12 @@ auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) {
} else if (item_size == 4) { } else if (item_size == 4) {
image = py::array_t<uint32_t>(shape); image = py::array_t<uint32_t>(shape);
} }
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()), n_frames, self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
header.mutable_data()); n_frames, header.mutable_data());
return py::make_tuple(header, image); return py::make_tuple(header, image);
} }
// Disable warnings for unused parameters, as we ignore some // Disable warnings for unused parameters, as we ignore some
// in the __exit__ method // in the __exit__ method
#pragma GCC diagnostic push #pragma GCC diagnostic push
@ -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_property_readonly("frames_in_file", &RawSubFile::frames_in_file)
.def("read_frame", &read_frame_from_RawSubFile) .def("read_frame", &read_frame_from_RawSubFile)
.def("read_n", &read_n_frames_from_RawSubFile) .def("read_n", &read_n_frames_from_RawSubFile)
.def("read", [](RawSubFile &self){ .def("read",
[](RawSubFile &self) {
self.seek(0); self.seek(0);
auto n_frames = self.frames_in_file(); auto n_frames = self.frames_in_file();
return read_n_frames_from_RawSubFile(self, n_frames); return read_n_frames_from_RawSubFile(self, n_frames);
}) })
.def("__enter__", [](RawSubFile &self) { return &self; }) .def("__enter__", [](RawSubFile &self) { return &self; })
.def("__exit__", .def("__exit__",
[](RawSubFile &self, [](RawSubFile &self, const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::object> &exc_value, const std::optional<pybind11::object> &exc_value,
const std::optional<pybind11::object> &traceback) { const std::optional<pybind11::object> &traceback) {})
})
.def("__iter__", [](RawSubFile &self) { return &self; }) .def("__iter__", [](RawSubFile &self) { return &self; })
.def("__next__", [](RawSubFile &self) { .def("__next__", [](RawSubFile &self) {
try { try {
@ -104,7 +100,6 @@ void define_raw_sub_file_io_bindings(py::module &m) {
throw py::stop_iteration(); throw py::stop_iteration();
} }
}); });
} }
#pragma GCC diagnostic pop #pragma GCC diagnostic pop

View File

@ -12,11 +12,9 @@
// #include <pybind11/stl/filesystem.h> // #include <pybind11/stl/filesystem.h>
// #include <string> // #include <string>
namespace py = pybind11; namespace py = pybind11;
using namespace ::aare; using namespace ::aare;
void define_var_cluster_finder_bindings(py::module &m) { void define_var_cluster_finder_bindings(py::module &m) {
PYBIND11_NUMPY_DTYPE(VarClusterFinder<double>::Hit, size, row, col, PYBIND11_NUMPY_DTYPE(VarClusterFinder<double>::Hit, size, row, col,
reserved, energy, max, rows, cols, enes); reserved, energy, max, rows, cols, enes);
@ -65,9 +63,7 @@ void define_var_cluster_finder_bindings(py::module &m) {
return return_vector(ptr); return return_vector(ptr);
}) })
.def("clear_hits", .def("clear_hits",
[](VarClusterFinder<double> &self) { [](VarClusterFinder<double> &self) { self.clear_hits(); })
self.clear_hits();
})
.def("steal_hits", .def("steal_hits",
[](VarClusterFinder<double> &self) { [](VarClusterFinder<double> &self) {
auto ptr = new std::vector<VarClusterFinder<double>::Hit>( auto ptr = new std::vector<VarClusterFinder<double>::Hit>(
@ -75,5 +71,4 @@ void define_var_cluster_finder_bindings(py::module &m) {
return return_vector(ptr); return return_vector(ptr);
}) })
.def("total_clusters", &VarClusterFinder<double>::total_clusters); .def("total_clusters", &VarClusterFinder<double>::total_clusters);
} }

View File

@ -31,9 +31,7 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
} }
} }
void ClusterFile::set_roi(ROI roi){ void ClusterFile::set_roi(ROI roi) { m_roi = roi; }
m_roi = roi;
}
void ClusterFile::set_noise_map(const NDView<int32_t, 2> noise_map) { void ClusterFile::set_noise_map(const NDView<int32_t, 2> noise_map) {
m_noise_map = NDArray<int32_t, 2>(noise_map); m_noise_map = NDArray<int32_t, 2>(noise_map);
@ -75,16 +73,17 @@ void ClusterFile::write_frame(const ClusterVector<int32_t> &clusters) {
// Then write the number of clusters - 4 bytes // Then write the number of clusters - 4 bytes
uint32_t n_clusters = clusters.size(); uint32_t n_clusters = clusters.size();
if (fwrite(&n_clusters, sizeof(n_clusters), 1, fp) != 1) { if (fwrite(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
throw std::runtime_error(LOCATION + "Could not write number of clusters"); throw std::runtime_error(LOCATION +
"Could not write number of clusters");
} }
// Now write the clusters in the frame // Now write the clusters in the frame
if(fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp)!=clusters.size()){ if (fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp) !=
clusters.size()) {
throw std::runtime_error(LOCATION + "Could not write clusters"); throw std::runtime_error(LOCATION + "Could not write clusters");
} }
} }
ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) { ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
if (m_mode != "r") { if (m_mode != "r") {
throw std::runtime_error("File not opened for reading"); throw std::runtime_error("File not opened for reading");
@ -96,7 +95,8 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters){
} }
} }
ClusterVector<int32_t> ClusterFile::read_clusters_without_cut(size_t n_clusters) { ClusterVector<int32_t>
ClusterFile::read_clusters_without_cut(size_t n_clusters) {
if (m_mode != "r") { if (m_mode != "r") {
throw std::runtime_error("File not opened for reading"); throw std::runtime_error("File not opened for reading");
} }
@ -152,7 +152,6 @@ ClusterVector<int32_t> ClusterFile::read_clusters_without_cut(size_t n_clusters)
return clusters; return clusters;
} }
ClusterVector<int32_t> ClusterFile::read_clusters_with_cut(size_t n_clusters) { ClusterVector<int32_t> ClusterFile::read_clusters_with_cut(size_t n_clusters) {
ClusterVector<int32_t> clusters(3, 3); ClusterVector<int32_t> clusters(3, 3);
clusters.reserve(n_clusters); clusters.reserve(n_clusters);
@ -162,7 +161,8 @@ ClusterVector<int32_t> ClusterFile::read_clusters_with_cut(size_t n_clusters) {
while (m_num_left && clusters.size() < n_clusters) { while (m_num_left && clusters.size() < n_clusters) {
Cluster3x3 c = read_one_cluster(); Cluster3x3 c = read_one_cluster();
if (is_selected(c)) { if (is_selected(c)) {
clusters.push_back(c.x, c.y, reinterpret_cast<std::byte*>(c.data)); clusters.push_back(c.x, c.y,
reinterpret_cast<std::byte *>(c.data));
} }
} }
} }
@ -172,17 +172,21 @@ ClusterVector<int32_t> ClusterFile::read_clusters_with_cut(size_t n_clusters) {
if (clusters.size() < n_clusters) { if (clusters.size() < n_clusters) {
// sanity check // sanity check
if (m_num_left) { 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! int32_t frame_number = 0; // frame number needs to be 4 bytes!
while (fread(&frame_number, sizeof(frame_number), 1, fp)) { while (fread(&frame_number, sizeof(frame_number), 1, fp)) {
if (fread(&m_num_left, sizeof(m_num_left), 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 clusters.set_frame_number(
frame_number); // cluster vector will hold the last frame
// number
while (m_num_left && clusters.size() < n_clusters) { while (m_num_left && clusters.size() < n_clusters) {
Cluster3x3 c = read_one_cluster(); Cluster3x3 c = read_one_cluster();
if (is_selected(c)) { if (is_selected(c)) {
clusters.push_back(c.x, c.y, reinterpret_cast<std::byte*>(c.data)); clusters.push_back(
c.x, c.y, reinterpret_cast<std::byte *>(c.data));
} }
} }
} }
@ -191,7 +195,6 @@ ClusterVector<int32_t> ClusterFile::read_clusters_with_cut(size_t n_clusters) {
if (clusters.size() >= n_clusters) if (clusters.size() >= n_clusters)
break; break;
} }
} }
if (m_gain_map) if (m_gain_map)
clusters.apply_gain_map(m_gain_map->view()); clusters.apply_gain_map(m_gain_map->view());
@ -235,7 +238,8 @@ ClusterVector<int32_t> ClusterFile::read_frame_without_cut() {
int32_t n_clusters; // Saved as 32bit integer in the cluster file int32_t n_clusters; // Saved as 32bit integer in the cluster file
if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) { 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<int32_t> clusters(3, 3, n_clusters); ClusterVector<int32_t> clusters(3, 3, n_clusters);
@ -264,7 +268,6 @@ ClusterVector<int32_t> ClusterFile::read_frame_with_cut() {
throw std::runtime_error("Could not read frame number"); throw std::runtime_error("Could not read frame number");
} }
if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) { if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) {
throw std::runtime_error("Could not read number of clusters"); throw std::runtime_error("Could not read number of clusters");
} }
@ -283,8 +286,6 @@ ClusterVector<int32_t> ClusterFile::read_frame_with_cut() {
return clusters; return clusters;
} }
bool ClusterFile::is_selected(Cluster3x3 &cl) { bool ClusterFile::is_selected(Cluster3x3 &cl) {
// Should fail fast // Should fail fast
if (m_roi) { if (m_roi) {
@ -297,7 +298,8 @@ bool ClusterFile::is_selected(Cluster3x3 &cl) {
int32_t sum_2x2 = cl.sum_2x2(); // highest sum of 2x2 subclusters 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) { if (sum_1x1 <= noise || sum_2x2 <= 2 * noise || sum_3x3 <= 3 * noise) {
return false; return false;
} }
@ -316,7 +318,8 @@ NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
eta2(i, 0) = e.x; eta2(i, 0) = e.x;
eta2(i, 1) = e.y; 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++) { for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2(clusters.at<Cluster2x2>(i)); auto e = calculate_eta2(clusters.at<Cluster2x2>(i));
eta2(i, 0) = e.x; eta2(i, 0) = e.x;
@ -330,8 +333,8 @@ NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
} }
/** /**
* @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2 struct * @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2
* containing etay, etax and the corner of the cluster. * struct containing etay, etax and the corner of the cluster.
*/ */
Eta2 calculate_eta2(Cluster3x3 &cl) { Eta2 calculate_eta2(Cluster3x3 &cl) {
Eta2 eta{}; Eta2 eta{};
@ -347,38 +350,30 @@ Eta2 calculate_eta2(Cluster3x3 &cl) {
switch (c) { switch (c) {
case cBottomLeft: case cBottomLeft:
if ((cl.data[3] + cl.data[4]) != 0) if ((cl.data[3] + cl.data[4]) != 0)
eta.x = eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[1] + cl.data[4]) != 0) if ((cl.data[1] + cl.data[4]) != 0)
eta.y = eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.c = cBottomLeft; eta.c = cBottomLeft;
break; break;
case cBottomRight: case cBottomRight:
if ((cl.data[2] + cl.data[5]) != 0) if ((cl.data[2] + cl.data[5]) != 0)
eta.x = eta.x = static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
if ((cl.data[1] + cl.data[4]) != 0) if ((cl.data[1] + cl.data[4]) != 0)
eta.y = eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.c = cBottomRight; eta.c = cBottomRight;
break; break;
case cTopLeft: case cTopLeft:
if ((cl.data[7] + cl.data[4]) != 0) if ((cl.data[7] + cl.data[4]) != 0)
eta.x = eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0) if ((cl.data[7] + cl.data[4]) != 0)
eta.y = eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.c = cTopLeft; eta.c = cTopLeft;
break; break;
case cTopRight: case cTopRight:
if ((cl.data[5] + cl.data[4]) != 0) if ((cl.data[5] + cl.data[4]) != 0)
eta.x = eta.x = static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0) if ((cl.data[7] + cl.data[4]) != 0)
eta.y = eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.c = cTopRight; eta.c = cTopRight;
break; break;
// no default to allow compiler to warn about missing cases // no default to allow compiler to warn about missing cases
@ -386,7 +381,6 @@ Eta2 calculate_eta2(Cluster3x3 &cl) {
return eta; return eta;
} }
Eta2 calculate_eta2(Cluster2x2 &cl) { Eta2 calculate_eta2(Cluster2x2 &cl) {
Eta2 eta{}; Eta2 eta{};
if ((cl.data[0] + cl.data[1]) != 0) if ((cl.data[0] + cl.data[1]) != 0)
@ -398,5 +392,4 @@ Eta2 calculate_eta2(Cluster2x2 &cl) {
return eta; return eta;
} }
} // namespace aare } // namespace aare

View File

@ -10,7 +10,6 @@ using aare::Cluster;
using aare::ClusterFile; using aare::ClusterFile;
using aare::ClusterVector; using aare::ClusterVector;
TEST_CASE("Read one frame from a cluster file", "[.files]") { 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"; auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
@ -27,7 +26,6 @@ TEST_CASE("Read one frame from a cluster file", "[.files]") {
std::begin(expected_cluster_data))); std::begin(expected_cluster_data)));
} }
TEST_CASE("Read one frame using ROI", "[.files]") { TEST_CASE("Read one frame using ROI", "[.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"; 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))); std::begin(expected_cluster_data)));
} }
TEST_CASE("Read clusters from single frame file", "[.files]") { TEST_CASE("Read clusters from single frame file", "[.files]") {
// frame_number, num_clusters [135] 97 // frame_number, num_clusters [135] 97

View File

@ -19,7 +19,8 @@ void CtbRawFile::read_into(std::byte *image_buf, DetectorHeader* header) {
throw std::runtime_error(LOCATION + " End of file reached"); throw std::runtime_error(LOCATION + " End of file reached");
} }
if(m_current_frame != 0 && m_current_frame % m_master.max_frames_per_file() == 0){ if (m_current_frame != 0 &&
m_current_frame % m_master.max_frames_per_file() == 0) {
open_data_file(m_current_subfile + 1); open_data_file(m_current_subfile + 1);
} }
@ -29,7 +30,8 @@ void CtbRawFile::read_into(std::byte *image_buf, DetectorHeader* header) {
m_file.seekg(sizeof(DetectorHeader), std::ios::cur); m_file.seekg(sizeof(DetectorHeader), std::ios::cur);
} }
m_file.read(reinterpret_cast<char *>(image_buf), m_master.image_size_in_bytes()); m_file.read(reinterpret_cast<char *>(image_buf),
m_master.image_size_in_bytes());
m_current_frame++; m_current_frame++;
} }
@ -38,13 +40,16 @@ void CtbRawFile::seek(size_t frame_number) {
open_data_file(index); open_data_file(index);
} }
size_t frame_number_in_file = frame_number % m_master.max_frames_per_file(); 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; m_current_frame = frame_number;
} }
size_t CtbRawFile::tell() const { return m_current_frame; } 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(); } 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"); throw std::runtime_error(LOCATION + "Subfile index out of range");
} }
m_current_subfile = subfile_index; 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()) { if (!m_file.is_open()) {
throw std::runtime_error(LOCATION + "Could not open data file"); throw std::runtime_error(LOCATION + "Could not open data file");
} }
} }
} // namespace aare } // namespace aare

View File

@ -10,7 +10,8 @@ namespace aare {
* @brief Construct a DType object from a type_info object * @brief Construct a DType object from a type_info object
* @param t type_info object * @param t type_info object
* @throw runtime_error if the type is not supported * @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)) * @note the type_info object is obtained using typeid (e.g. typeid(int))
*/ */
Dtype::Dtype(const std::type_info &t) { Dtype::Dtype(const std::type_info &t) {
@ -35,7 +36,8 @@ Dtype::Dtype(const std::type_info &t) {
else if (t == typeid(double)) else if (t == typeid(double))
m_type = TypeIndex::DOUBLE; m_type = TypeIndex::DOUBLE;
else 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: case TypeIndex::NONE:
return 0; return 0;
default: 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: case 64:
return Dtype(TypeIndex::UINT64); return Dtype(TypeIndex::UINT64);
default: 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: case TypeIndex::DOUBLE:
return "f8"; return "f8";
case TypeIndex::ERROR: 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: 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 {}; return {};
} }
bool Dtype::operator==(const Dtype &other) const noexcept { return m_type == other.m_type; } bool Dtype::operator==(const Dtype &other) const noexcept {
bool Dtype::operator!=(const Dtype &other) const noexcept { return !(*this == other); } 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 {
bool Dtype::operator!=(const std::type_info &t) const { return Dtype(t) != *this; } return Dtype(t) == *this;
}
bool Dtype::operator!=(const std::type_info &t) const {
return Dtype(t) != *this;
}
} // namespace aare } // namespace aare

View File

@ -51,4 +51,6 @@ TEST_CASE("Construct from string with endianess") {
REQUIRE_THROWS(Dtype(">i4") == typeid(int32_t)); REQUIRE_THROWS(Dtype(">i4") == typeid(int32_t));
} }
TEST_CASE("Convert to string") { REQUIRE(Dtype(typeid(int)).to_string() == "<i4"); } TEST_CASE("Convert to string") {
REQUIRE(Dtype(typeid(int)).to_string() == "<i4");
}

View File

@ -27,8 +27,7 @@ File::File(const std::filesystem::path &fname, const std::string &mode,
if (fname.extension() == ".raw" || fname.extension() == ".json") { if (fname.extension() == ".raw" || fname.extension() == ".json") {
// file_impl = new RawFile(fname, mode, cfg); // file_impl = new RawFile(fname, mode, cfg);
file_impl = std::make_unique<RawFile>(fname, mode); file_impl = std::make_unique<RawFile>(fname, mode);
} } else if (fname.extension() == ".npy") {
else if (fname.extension() == ".npy") {
// file_impl = new NumpyFile(fname, mode, cfg); // file_impl = new NumpyFile(fname, mode, cfg);
file_impl = std::make_unique<NumpyFile>(fname, mode, cfg); file_impl = std::make_unique<NumpyFile>(fname, mode, cfg);
} }
@ -48,10 +47,7 @@ File::File(const std::filesystem::path &fname, const std::string &mode,
} }
} }
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) { if (this != &other) {
@ -89,9 +85,10 @@ size_t File::tell() const { return file_impl->tell(); }
size_t File::rows() const { return file_impl->rows(); } size_t File::rows() const { return file_impl->rows(); }
size_t File::cols() const { return file_impl->cols(); } size_t File::cols() const { return file_impl->cols(); }
size_t File::bitdepth() const { return file_impl->bitdepth(); } 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(); } DetectorType File::detector_type() const { return file_impl->detector_type(); }
} // namespace aare } // namespace aare

View File

@ -6,10 +6,12 @@
namespace aare { 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()); fp_ = fopen(fname.c_str(), mode.c_str());
if (!fp_) 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_); } FilePtr::FilePtr(FilePtr &&other) { std::swap(fp_, other.fp_); }
@ -24,7 +26,8 @@ FILE *FilePtr::get() { return fp_; }
ssize_t FilePtr::tell() { ssize_t FilePtr::tell() {
auto pos = ftell(fp_); auto pos = ftell(fp_);
if (pos == -1) 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; return pos;
} }
FilePtr::~FilePtr() { FilePtr::~FilePtr() {

View File

@ -1,13 +1,12 @@
#include "aare/Fit.hpp" #include "aare/Fit.hpp"
#include "aare/utils/task.hpp"
#include "aare/utils/par.hpp" #include "aare/utils/par.hpp"
#include "aare/utils/task.hpp"
#include <lmcurve2.h> #include <lmcurve2.h>
#include <lmfit.hpp> #include <lmfit.hpp>
#include <thread> #include <thread>
#include <array> #include <array>
namespace aare { namespace aare {
namespace func { namespace func {
@ -35,7 +34,9 @@ NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
} }
double scurve(const double x, const double *par) { 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])); 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<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par) { NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par) {
@ -47,7 +48,9 @@ NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par) {
} }
double scurve2(const double x, const double *par) { 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])); 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<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par) { NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par) {
@ -91,7 +94,8 @@ NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
return result; return result;
} }
std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<double, 1> y) { std::array<double, 3> gaus_init_par(const NDView<double, 1> x,
const NDView<double, 1> y) {
std::array<double, 3> start_par{0, 0, 0}; std::array<double, 3> start_par{0, 0, 0};
auto e = std::max_element(y.begin(), y.end()); auto e = std::max_element(y.begin(), y.end());
auto idx = std::distance(y.begin(), e); auto idx = std::distance(y.begin(), e);
@ -103,20 +107,18 @@ std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<doub
// For sigma we estimate the fwhm and divide by 2.35 // For sigma we estimate the fwhm and divide by 2.35
// assuming equally spaced x values // assuming equally spaced x values
auto delta = x[1] - x[0]; auto delta = x[1] - x[0];
start_par[2] = start_par[2] = std::count_if(y.begin(), y.end(),
std::count_if(y.begin(), y.end(),
[e](double val) { return val > *e / 2; }) * [e](double val) { return val > *e / 2; }) *
delta / 2.35; delta / 2.35;
return start_par; return start_par;
} }
std::array<double, 2> pol1_init_par(const NDView<double, 1> x,
std::array<double, 2> pol1_init_par(const NDView<double, 1> x, const NDView<double, 1> y){ const NDView<double, 1> y) {
// Estimate the initial parameters for the fit // Estimate the initial parameters for the fit
std::array<double, 2> start_par{0, 0}; std::array<double, 2> start_par{0, 0};
auto y2 = std::max_element(y.begin(), y.end()); auto y2 = std::max_element(y.begin(), y.end());
auto x2 = x[std::distance(y.begin(), y2)]; auto x2 = x[std::distance(y.begin(), y2)];
auto y1 = std::min_element(y.begin(), y.end()); auto y1 = std::min_element(y.begin(), y.end());
@ -141,7 +143,6 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
"and par_out, par_err_out must have size 3"); "and par_out, par_err_out must have size 3");
} }
// /* Collection of output parameters for status info. */ // /* Collection of output parameters for status info. */
// typedef struct { // typedef struct {
// double fnorm; /* norm of the residue vector fvec. */ // double fnorm; /* norm of the residue vector fvec. */
@ -153,23 +154,32 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
// */ // */
// } lm_status_struct; // } lm_status_struct;
lm_status_struct status; lm_status_struct status;
par_out = gaus_init_par(x, y); par_out = gaus_init_par(x, y);
std::array<double, 9> cov{0, 0, 0, 0, 0, 0, 0, 0, 0}; std::array<double, 9> 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); // void lmcurve2( const int n_par, double *par, double *parerr, double
// n_par - Number of free variables. Length of parameter vector par. // *covar, const int m_dat, const double *t, const double *y, const double
// par - Parameter vector. On input, it must contain a reasonable guess. On output, it contains the solution found to minimize ||r||. // *dy, double (*f)( const double ti, const double *par ), const
// 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. // lm_control_struct *control, lm_status_struct *status); n_par - Number of
// covar - Covariance matrix. Array of length n_par * n_par or NULL. On output, unless it is NULL, it contains the covariance matrix. // free variables. Length of parameter vector par. par - Parameter vector.
// m_dat - Number of data points. Length of vectors t, y, dy. Must statisfy n_par <= m_dat. // On input, it must contain a reasonable guess. On output, it contains the
// t - Array of length m_dat. Contains the abcissae (time, or "x") for which function f will be evaluated. // solution found to minimize ||r||. parerr - Parameter uncertainties
// y - Array of length m_dat. Contains the ordinate values that shall be fitted. // vector. Array of length n_par or NULL. On output, unless it or covar is
// dy - Array of length m_dat. Contains the standard deviations of the values y. // NULL, it contains the weighted parameter uncertainties for the found
// f - A user-supplied parametric function f(ti;par). // parameters. covar - Covariance matrix. Array of length n_par * n_par or
// 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). // NULL. On output, unless it is NULL, it contains the covariance matrix.
// status - A record used to return information about the minimization process: For details, see lmmin2(3). // 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(), 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, x.size(), x.data(), y.data(), y_err.data(), aare::func::gaus,
@ -178,12 +188,14 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
// Calculate chi2 // Calculate chi2
chi2 = 0; chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) { 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<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err, void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out, NDView<double, 3> par_out, NDView<double, 3> par_err_out,
NDView<double, 2> chi2_out,
int n_threads) { int n_threads) {
@ -200,7 +212,6 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
fit_gaus(x, y_view, y_err_view, par_out_view, par_err_out_view, fit_gaus(x, y_view, y_err_view, par_out_view, par_err_out_view,
chi2_out(row, col)); chi2_out(row, col));
} }
} }
}; };
@ -210,7 +221,8 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
} }
void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err, void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) { NDView<double, 1> par_out, NDView<double, 1> par_err_out,
double &chi2) {
// Check that we have the correct sizes // Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() || if (y.size() != x.size() || y.size() != y_err.size() ||
@ -230,13 +242,14 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
// Calculate chi2 // Calculate chi2
chi2 = 0; chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) { 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<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err, void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out, NDView<double, 3> par_out, NDView<double, 3> par_err_out,
int n_threads) { NDView<double, 2> chi2_out, int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) { auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) { for (ssize_t row = first_row; row < last_row; row++) {
@ -249,15 +262,14 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0), NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)}); {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); auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks); RunInParallel(process, tasks);
} }
NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y) { NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y) {
@ -300,7 +312,8 @@ NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
// ~~ S-CURVES ~~ // ~~ S-CURVES ~~
// SCURVE -- // SCURVE --
std::array<double, 6> scurve_init_par(const NDView<double, 1> x, const NDView<double, 1> y){ std::array<double, 6> scurve_init_par(const NDView<double, 1> x,
const NDView<double, 1> y) {
// Estimate the initial parameters for the fit // Estimate the initial parameters for the fit
std::array<double, 6> start_par{0, 0, 0, 0, 0, 0}; std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
@ -308,7 +321,8 @@ std::array<double, 6> scurve_init_par(const NDView<double, 1> x, const NDView<do
auto ymin = std::min_element(y.begin(), y.end()); auto ymin = std::min_element(y.begin(), y.end());
start_par[4] = *ymin + (*ymax - *ymin) / 2; start_par[4] = *ymin + (*ymax - *ymin) / 2;
// Find the first x where the corresponding y value is above the threshold (start_par[4]) // 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) { for (ssize_t i = 0; i < y.size(); ++i) {
if (y[i] >= start_par[4]) { if (y[i] >= start_par[4]) {
start_par[2] = x[i]; start_par[2] = x[i];
@ -334,7 +348,8 @@ NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y) {
return result; return result;
} }
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_threads) { NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0); NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) { auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
@ -358,8 +373,9 @@ NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_th
} }
// - Error // - Error
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err, void fit_scurve(NDView<double, 1> x, NDView<double, 1> y,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) { NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 1> par_err_out, double &chi2) {
// Check that we have the correct sizes // Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() || if (y.size() != x.size() || y.size() != y_err.size() ||
@ -380,12 +396,14 @@ void fit_scurve(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_er
// Calculate chi2 // Calculate chi2
chi2 = 0; chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) { 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<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err, void fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out, NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) { int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) { auto process = [&](ssize_t first_row, ssize_t last_row) {
@ -399,20 +417,20 @@ void fit_scurve(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_er
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0), NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)}); {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); auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks); RunInParallel(process, tasks);
} }
// SCURVE2 --- // SCURVE2 ---
std::array<double, 6> scurve2_init_par(const NDView<double, 1> x, const NDView<double, 1> y){ std::array<double, 6> scurve2_init_par(const NDView<double, 1> x,
const NDView<double, 1> y) {
// Estimate the initial parameters for the fit // Estimate the initial parameters for the fit
std::array<double, 6> start_par{0, 0, 0, 0, 0, 0}; std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
@ -420,7 +438,8 @@ std::array<double, 6> scurve2_init_par(const NDView<double, 1> x, const NDView<d
auto ymin = std::min_element(y.begin(), y.end()); auto ymin = std::min_element(y.begin(), y.end());
start_par[4] = *ymin + (*ymax - *ymin) / 2; start_par[4] = *ymin + (*ymax - *ymin) / 2;
// Find the first x where the corresponding y value is above the threshold (start_par[4]) // 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) { for (ssize_t i = 0; i < y.size(); ++i) {
if (y[i] <= start_par[4]) { if (y[i] <= start_par[4]) {
start_par[2] = x[i]; start_par[2] = x[i];
@ -446,7 +465,8 @@ NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y) {
return result; return result;
} }
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, int n_threads) { NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0); NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) { auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
@ -470,8 +490,9 @@ NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, int n_t
} }
// - Error // - Error
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err, void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) { NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 1> par_err_out, double &chi2) {
// Check that we have the correct sizes // Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() || if (y.size() != x.size() || y.size() != y_err.size() ||
@ -492,12 +513,14 @@ void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_e
// Calculate chi2 // Calculate chi2
chi2 = 0; chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) { 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<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err, void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out, NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) { int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) { auto process = [&](ssize_t first_row, ssize_t last_row) {
@ -511,15 +534,14 @@ void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_e
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0), NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)}); {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); auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks); RunInParallel(process, tasks);
} }
} // namespace aare } // namespace aare

View File

@ -29,7 +29,6 @@ uint64_t Frame::size() const { return m_rows * m_cols; }
size_t Frame::bytes() const { return m_rows * m_cols * m_dtype.bytes(); } size_t Frame::bytes() const { return m_rows * m_cols * m_dtype.bytes(); }
std::byte *Frame::data() const { return m_data; } 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)) { if ((row >= m_rows) || (col >= m_cols)) {
std::cerr << "Invalid row or column index" << '\n'; std::cerr << "Invalid row or column index" << '\n';
@ -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()); return m_data + (row * m_cols + col) * (m_dtype.bytes());
} }
Frame &Frame::operator=(Frame &&other) noexcept { Frame &Frame::operator=(Frame &&other) noexcept {
if (this == &other) { if (this == &other) {
return *this; return *this;
@ -70,5 +68,4 @@ Frame Frame::clone() const {
return frame; return frame;
} }
} // namespace aare } // namespace aare

View File

@ -65,7 +65,8 @@ TEST_CASE("Set a value in a 64 bit frame") {
// only the value we did set should be non-zero // only the value we did set should be non-zero
for (size_t i = 0; i < rows; i++) { for (size_t i = 0; i < rows; i++) {
for (size_t j = 0; j < cols; j++) { for (size_t j = 0; j < cols; j++) {
uint64_t *data = reinterpret_cast<uint64_t *>(frame.pixel_ptr(i, j)); uint64_t *data =
reinterpret_cast<uint64_t *>(frame.pixel_ptr(i, j));
REQUIRE(data != nullptr); REQUIRE(data != nullptr);
if (i == 5 && j == 7) { if (i == 5 && j == 7) {
REQUIRE(*data == value); REQUIRE(*data == value);
@ -150,4 +151,3 @@ TEST_CASE("test explicit copy constructor") {
REQUIRE(frame2.bytes() == rows * cols * bitdepth / 8); REQUIRE(frame2.bytes() == rows * cols * bitdepth / 8);
REQUIRE(frame2.data() != data); REQUIRE(frame2.data() != data);
} }

View File

@ -19,7 +19,6 @@ JungfrauDataFile::JungfrauDataFile(const std::filesystem::path &fname) {
open_file(m_current_file_index); open_file(m_current_file_index);
} }
// FileInterface // FileInterface
Frame JungfrauDataFile::read_frame() { Frame JungfrauDataFile::read_frame() {
@ -59,7 +58,9 @@ std::array<ssize_t, 2> JungfrauDataFile::shape() const {
return {static_cast<ssize_t>(rows()), static_cast<ssize_t>(cols())}; return {static_cast<ssize_t>(rows()), static_cast<ssize_t>(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; } std::string JungfrauDataFile::base_name() const { return m_base_name; }
@ -202,15 +203,16 @@ void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames,
} }
} }
void JungfrauDataFile::read_into(NDArray<uint16_t>* image, JungfrauDataHeader* header) { void JungfrauDataFile::read_into(NDArray<uint16_t> *image,
JungfrauDataHeader *header) {
if (image->shape() != shape()) { if (image->shape() != shape()) {
throw std::runtime_error(LOCATION + throw std::runtime_error(
"Image shape does not match file size: " + std::to_string(rows()) + "x" + std::to_string(cols())); LOCATION + "Image shape does not match file size: " +
std::to_string(rows()) + "x" + std::to_string(cols()));
} }
read_into(reinterpret_cast<std::byte *>(image->data()), header); read_into(reinterpret_cast<std::byte *>(image->data()), header);
} }
JungfrauDataHeader JungfrauDataFile::read_header() { JungfrauDataHeader JungfrauDataFile::read_header() {
JungfrauDataHeader header; JungfrauDataHeader header;
if (auto rc = fread(&header, 1, sizeof(header), m_fp.get()); if (auto rc = fread(&header, 1, sizeof(header), m_fp.get());

View File

@ -1,7 +1,7 @@
#include "aare/JungfrauDataFile.hpp" #include "aare/JungfrauDataFile.hpp"
#include <catch2/catch_test_macros.hpp>
#include "test_config.hpp" #include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
using aare::JungfrauDataFile; using aare::JungfrauDataFile;
using aare::JungfrauDataHeader; using aare::JungfrauDataHeader;
@ -92,7 +92,6 @@ TEST_CASE("Open a Jungfrau data file with non zero file index", "[.files]"){
REQUIRE(f.n_files() == 4); REQUIRE(f.n_files() == 4);
REQUIRE(f.current_file().stem() == "AldoJF65k_000003"); 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]") {
@ -109,6 +108,4 @@ TEST_CASE("Read into throws if size doesn't match", "[.files]"){
REQUIRE_THROWS(f.read_into(&image)); REQUIRE_THROWS(f.read_into(&image));
REQUIRE(f.tell() == 0); REQUIRE(f.tell() == 0);
} }

View File

@ -162,7 +162,6 @@ NDArray<int> AddNDArrayUsingIndex(NDArray<int> &a, NDArray<int> &b) {
return res; return res;
} }
TEST_CASE("Compare two images") { TEST_CASE("Compare two images") {
NDArray<int> a; NDArray<int> a;
NDArray<int> b; NDArray<int> b;
@ -222,7 +221,6 @@ TEST_CASE("Bitwise and on data") {
REQUIRE(a(2) == 384); REQUIRE(a(2) == 384);
} }
TEST_CASE("Elementwise operations on images") { TEST_CASE("Elementwise operations on images") {
std::array<ssize_t, 2> shape{5, 5}; std::array<ssize_t, 2> shape{5, 5};
double a_val = 3.0; double a_val = 3.0;
@ -258,7 +256,8 @@ TEST_CASE("Elementwise operations on images") {
NDArray<double> A(shape, a_val); NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val); NDArray<double> B(shape, b_val);
NDArray<double> C = A - B; NDArray<double> 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 // Value of C matches
for (uint32_t i = 0; i < C.size(); ++i) { for (uint32_t i = 0; i < C.size(); ++i) {
@ -282,7 +281,8 @@ TEST_CASE("Elementwise operations on images") {
SECTION("Multiply two images") { SECTION("Multiply two images") {
NDArray<double> A(shape, a_val); NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val); NDArray<double> 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<double> C = A * B; NDArray<double> C = A * B;
// Value of C matches // Value of C matches
@ -307,7 +307,8 @@ TEST_CASE("Elementwise operations on images") {
SECTION("Divide two images") { SECTION("Divide two images") {
NDArray<double> A(shape, a_val); NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val); NDArray<double> 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<double> C = A / B; NDArray<double> C = A / B;
// Value of C matches // Value of C matches

View File

@ -2,8 +2,8 @@
#include <catch2/catch_test_macros.hpp> #include <catch2/catch_test_macros.hpp>
#include <iostream> #include <iostream>
#include <vector>
#include <numeric> #include <numeric>
#include <vector>
using aare::NDView; using aare::NDView;
using aare::Shape; using aare::Shape;
@ -151,8 +151,10 @@ TEST_CASE("divide with another span") {
std::vector<int> vec1{3, 2, 1}; std::vector<int> vec1{3, 2, 1};
std::vector<int> result{3, 6, 3}; std::vector<int> result{3, 6, 3};
NDView<int, 1> data0(vec0.data(), Shape<1>{static_cast<ssize_t>(vec0.size())}); NDView<int, 1> data0(vec0.data(),
NDView<int, 1> data1(vec1.data(), Shape<1>{static_cast<ssize_t>(vec1.size())}); Shape<1>{static_cast<ssize_t>(vec0.size())});
NDView<int, 1> data1(vec1.data(),
Shape<1>{static_cast<ssize_t>(vec1.size())});
data0 /= data1; data0 /= data1;
@ -181,7 +183,6 @@ TEST_CASE("compare two views") {
REQUIRE((view1 == view2)); REQUIRE((view1 == view2));
} }
TEST_CASE("Create a view over a vector") { TEST_CASE("Create a view over a vector") {
std::vector<int> vec(12); std::vector<int> vec(12);
std::iota(vec.begin(), vec.end(), 0); std::iota(vec.begin(), vec.end(), 0);

View File

@ -4,16 +4,16 @@
namespace aare { 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 // TODO! add opts to constructor
m_mode = mode; m_mode = mode;
if (mode == "r") { if (mode == "r") {
fp = fopen(fname.string().c_str(), "rb"); fp = fopen(fname.string().c_str(), "rb");
if (!fp) { 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(); load_metadata();
} else if (mode == "w") { } 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}; m_header.shape = {0, cfg.rows, cfg.cols};
fp = fopen(fname.string().c_str(), "wb"); fp = fopen(fname.string().c_str(), "wb");
if (!fp) { 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; 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]) { if (frame_number > m_header.shape[0]) {
throw std::invalid_argument("Frame number out of range"); 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"); throw std::runtime_error("Could not seek to frame");
size_t const rc = fread(image_buf, m_bytes_per_frame, 1, fp); size_t const rc = fread(image_buf, m_bytes_per_frame, 1, fp);
@ -113,7 +118,8 @@ NumpyFile::~NumpyFile() noexcept {
// write header // write header
size_t const rc = fwrite(header_str.c_str(), header_str.size(), 1, fp); size_t const rc = fwrite(header_str.c_str(), header_str.size(), 1, fp);
if (rc != 1) { 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 // read version
rc = fread(reinterpret_cast<char *>(&major_ver_), sizeof(major_ver_), 1, fp); rc =
rc += fread(reinterpret_cast<char *>(&minor_ver_), sizeof(minor_ver_), 1, fp); fread(reinterpret_cast<char *>(&major_ver_), sizeof(major_ver_), 1, fp);
rc +=
fread(reinterpret_cast<char *>(&minor_ver_), sizeof(minor_ver_), 1, fp);
if (rc != 2) { if (rc != 2) {
throw std::runtime_error("Error reading numpy version"); throw std::runtime_error("Error reading numpy version");
} }
@ -159,7 +167,8 @@ void NumpyFile::load_metadata() {
if (rc != 1) { if (rc != 1) {
throw std::runtime_error("Error reading header length"); 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) { if (header_size % 16 != 0) {
fmt::print("Warning: header length is not a multiple of 16\n"); fmt::print("Warning: header length is not a multiple of 16\n");
} }

View File

@ -1,8 +1,8 @@
#include "aare/NumpyFile.hpp" #include "aare/NumpyFile.hpp"
#include "aare/NDArray.hpp" #include "aare/NDArray.hpp"
#include <catch2/catch_test_macros.hpp>
#include "test_config.hpp" #include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
using aare::Dtype; using aare::Dtype;
using aare::NumpyFile; using aare::NumpyFile;

View File

@ -29,7 +29,8 @@ namespace aare {
std::string NumpyHeader::to_string() const { std::string NumpyHeader::to_string() const {
std::stringstream sstm; std::stringstream sstm;
sstm << "dtype: " << dtype.to_string() << ", fortran_order: " << fortran_order << ' '; sstm << "dtype: " << dtype.to_string()
<< ", fortran_order: " << fortran_order << ' ';
sstm << "shape: ("; sstm << "shape: (";
for (auto item : shape) for (auto item : shape)
sstm << item << ','; sstm << item << ',';
@ -37,10 +38,10 @@ std::string NumpyHeader::to_string() const {
return sstm.str(); return sstm.str();
} }
namespace NumpyHelpers { namespace NumpyHelpers {
std::unordered_map<std::string, std::string> parse_dict(std::string in, const std::vector<std::string> &keys) { std::unordered_map<std::string, std::string>
parse_dict(std::string in, const std::vector<std::string> &keys) {
std::unordered_map<std::string, std::string> map; std::unordered_map<std::string, std::string> map;
if (keys.empty()) if (keys.empty())
return map; return map;
@ -100,7 +101,8 @@ aare::Dtype parse_descr(std::string typestring) {
constexpr char little_endian_char = '<'; constexpr char little_endian_char = '<';
constexpr char big_endian_char = '>'; constexpr char big_endian_char = '>';
constexpr char no_endian_char = '|'; constexpr char no_endian_char = '|';
constexpr std::array<char, 3> endian_chars = {little_endian_char, big_endian_char, no_endian_char}; constexpr std::array<char, 3> endian_chars = {
little_endian_char, big_endian_char, no_endian_char};
constexpr std::array<char, 4> numtype_chars = {'f', 'i', 'u', 'c'}; constexpr std::array<char, 4> numtype_chars = {'f', 'i', 'u', 'c'};
const char byteorder_c = typestring[0]; const char byteorder_c = typestring[0];
@ -139,7 +141,9 @@ std::string get_value_from_map(const std::string &mapstr) {
return trim(tmp); 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<std::string> parse_tuple(std::string in) { std::vector<std::string> parse_tuple(std::string in) {
std::vector<std::string> v; std::vector<std::string> v;
@ -215,20 +219,25 @@ inline std::string write_boolean(bool b) {
return "False"; return "False";
} }
inline std::string write_header_dict(const std::string &descr, bool fortran_order, const std::vector<size_t> &shape) { inline std::string write_header_dict(const std::string &descr,
bool fortran_order,
const std::vector<size_t> &shape) {
std::string const s_fortran_order = write_boolean(fortran_order); std::string const s_fortran_order = write_boolean(fortran_order);
std::string const shape_s = write_tuple(shape); 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); std::ofstream out(fname, std::ios::binary | std::ios::out);
return write_header(out, header); return write_header(out, header);
} }
size_t write_header(std::ostream &out, const NumpyHeader &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; 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 // write header length
if (version_major == 1 && version_minor == 0) { if (version_major == 1 && version_minor == 0) {
auto header_len = static_cast<uint16_t>(header_dict.length() + padding.length() + 1); auto header_len =
static_cast<uint16_t>(header_dict.length() + padding.length() + 1);
std::array<uint8_t, 2> header_len_le16{static_cast<uint8_t>((header_len >> 0) & 0xff), std::array<uint8_t, 2> header_len_le16{
static_cast<uint8_t>((header_len >> 0) & 0xff),
static_cast<uint8_t>((header_len >> 8) & 0xff)}; static_cast<uint8_t>((header_len >> 8) & 0xff)};
out.write(reinterpret_cast<char *>(header_len_le16.data()), 2); out.write(reinterpret_cast<char *>(header_len_le16.data()), 2);
} else { } else {
auto header_len = static_cast<uint32_t>(header_dict.length() + padding.length() + 1); auto header_len =
static_cast<uint32_t>(header_dict.length() + padding.length() + 1);
std::array<uint8_t, 4> header_len_le32{ std::array<uint8_t, 4> header_len_le32{
static_cast<uint8_t>((header_len >> 0) & 0xff), static_cast<uint8_t>((header_len >> 8) & 0xff), static_cast<uint8_t>((header_len >> 0) & 0xff),
static_cast<uint8_t>((header_len >> 16) & 0xff), static_cast<uint8_t>((header_len >> 24) & 0xff)}; static_cast<uint8_t>((header_len >> 8) & 0xff),
static_cast<uint8_t>((header_len >> 16) & 0xff),
static_cast<uint8_t>((header_len >> 24) & 0xff)};
out.write(reinterpret_cast<char *>(header_len_le32.data()), 4); out.write(reinterpret_cast<char *>(header_len_le32.data()), 4);
} }

View File

@ -19,7 +19,9 @@ TEST_CASE("Check for quotes and return stripped string") {
REQUIRE(parse_str("''") == ""); 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") { TEST_CASE("trim whitespace") {
REQUIRE(trim(" hej ") == "hej"); REQUIRE(trim(" hej ") == "hej");
@ -53,7 +55,8 @@ TEST_CASE("is element in array") {
} }
TEST_CASE("Parse numpy dict") { TEST_CASE("Parse numpy dict") {
std::string in = "{'descr': '<f4', 'fortran_order': False, 'shape': (3, 4)}"; std::string in =
"{'descr': '<f4', 'fortran_order': False, 'shape': (3, 4)}";
std::vector<std::string> keys{"descr", "fortran_order", "shape"}; std::vector<std::string> keys{"descr", "fortran_order", "shape"};
auto map = parse_dict(in, keys); auto map = parse_dict(in, keys);
REQUIRE(map["descr"] == "'<f4'"); REQUIRE(map["descr"] == "'<f4'");

View File

@ -1,8 +1,7 @@
#include "aare/Pedestal.hpp" #include "aare/Pedestal.hpp"
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/catch_test_macros.hpp> #include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <chrono> #include <chrono>
#include <random> #include <random>
@ -58,7 +57,8 @@ TEST_CASE("test pedestal push") {
if (k < 5) { if (k < 5) {
REQUIRE(pedestal.cur_samples()(i, j) == k + 1); REQUIRE(pedestal.cur_samples()(i, j) == k + 1);
REQUIRE(pedestal.get_sum()(i, j) == (k + 1) * (i + j)); 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 { } else {
REQUIRE(pedestal.cur_samples()(i, j) == 5); REQUIRE(pedestal.cur_samples()(i, j) == 5);
REQUIRE(pedestal.get_sum()(i, j) == 5 * (i + j)); 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 i = 0; i < 3; i++) {
for (int j = 0; j < 5; j++) { for (int j = 0; j < 5; j++) {
REQUIRE_THAT(mean(i, j), Catch::Matchers::WithinAbs(MEAN, MEAN * TOLERANCE)); REQUIRE_THAT(mean(i, j),
REQUIRE_THAT(variance(i, j), Catch::Matchers::WithinAbs(VAR, VAR * TOLERANCE)); Catch::Matchers::WithinAbs(MEAN, MEAN * TOLERANCE));
REQUIRE_THAT(standard_deviation(i, j), Catch::Matchers::WithinAbs(STD, STD * TOLERANCE)); REQUIRE_THAT(variance(i, j),
Catch::Matchers::WithinAbs(VAR, VAR * TOLERANCE));
REQUIRE_THAT(standard_deviation(i, j),
Catch::Matchers::WithinAbs(STD, STD * TOLERANCE));
} }
} }
} }

View File

@ -42,9 +42,9 @@ NDArray<ssize_t, 2> GenerateMoench05PixelMap() {
int adc_nr = adc_numbers[i_sc]; 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; order_map(row, col) = i_analog;
} }
} }
} }
@ -63,10 +63,9 @@ NDArray<ssize_t, 2> GenerateMoench05PixelMap1g() {
int adc_nr = adc_numbers[i_sc]; 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] &
// analog_frame[row * 150 + col] = analog_data[i_analog] & 0x3FFF; // 0x3FFF;
order_map(row, col) = i_analog; order_map(row, col) = i_analog;
} }
} }
} }
@ -85,10 +84,9 @@ NDArray<ssize_t, 2> GenerateMoench05PixelMapOld() {
int adc_nr = adc_numbers[i_sc]; int adc_nr = adc_numbers[i_sc];
int i_analog = n_pixel * 32 + adc_nr; int i_analog = n_pixel * 32 + adc_nr;
// analog_frame[row * 150 + col] = analog_data[i_analog] &
// analog_frame[row * 150 + col] = analog_data[i_analog] & 0x3FFF; // 0x3FFF;
order_map(row, col) = i_analog; order_map(row, col) = i_analog;
} }
} }
} }
@ -120,7 +118,8 @@ NDArray<ssize_t, 3> GenerateMH02FourCounterPixelMap(){
for (int counter = 0; counter < 4; counter++) { for (int counter = 0; counter < 4; counter++) {
for (int row = 0; row < 48; row++) { for (int row = 0; row < 48; row++) {
for (int col = 0; col < 48; col++) { for (int col = 0; col < 48; col++) {
order_map(counter, row, col) = counter*48*48 + row*48 + col; order_map(counter, row, col) =
counter * 48 * 48 + row * 48 + col;
} }
} }
} }

View File

@ -1,9 +1,9 @@
#include "aare/RawFile.hpp" #include "aare/RawFile.hpp"
#include "aare/algorithm.hpp"
#include "aare/PixelMap.hpp" #include "aare/PixelMap.hpp"
#include "aare/algorithm.hpp"
#include "aare/defs.hpp" #include "aare/defs.hpp"
#include "aare/logger.hpp"
#include "aare/geo_helpers.hpp" #include "aare/geo_helpers.hpp"
#include "aare/logger.hpp"
#include <fmt/format.h> #include <fmt/format.h>
#include <nlohmann/json.hpp> #include <nlohmann/json.hpp>
@ -18,7 +18,8 @@ RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode)
if (mode == "r") { if (mode == "r") {
find_geometry(); find_geometry();
if (m_master.roi()) { if (m_master.roi()) {
m_geometry = update_geometry_with_roi(m_geometry, m_master.roi().value()); m_geometry =
update_geometry_with_roi(m_geometry, m_master.roi().value());
} }
open_subfiles(); open_subfiles();
} else { } else {
@ -47,13 +48,13 @@ void RawFile::read_into(std::byte *image_buf) {
return get_frame_into(m_current_frame++, image_buf); return get_frame_into(m_current_frame++, image_buf);
} }
void RawFile::read_into(std::byte *image_buf, DetectorHeader *header) { void RawFile::read_into(std::byte *image_buf, DetectorHeader *header) {
return get_frame_into(m_current_frame++, image_buf, 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); // return get_frame_into(m_current_frame++, image_buf, header);
for (size_t i = 0; i < n_frames; i++) { for (size_t i = 0; i < n_frames; i++) {
@ -62,14 +63,13 @@ void RawFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *h
if (header) if (header)
header += n_modules(); header += n_modules();
} }
} }
size_t RawFile::n_modules() const { return m_master.n_modules(); } size_t RawFile::n_modules() const { return m_master.n_modules(); }
size_t RawFile::bytes_per_frame() { 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() { size_t RawFile::pixels_per_frame() {
// return m_rows * m_cols; // return m_rows * m_cols;
@ -128,7 +128,6 @@ DetectorHeader RawFile::read_header(const std::filesystem::path &fname) {
return h; return h;
} }
RawMasterFile RawFile::master() const { return m_master; } RawMasterFile RawFile::master() const { return m_master; }
/** /**
@ -142,7 +141,6 @@ void RawFile::find_geometry() {
uint16_t r{}; uint16_t r{};
uint16_t c{}; uint16_t c{};
for (size_t i = 0; i < n_modules(); i++) { for (size_t i = 0; i < n_modules(); i++) {
auto h = read_header(m_master.data_fname(i, 0)); auto h = read_header(m_master.data_fname(i, 0));
r = std::max(r, h.row); r = std::max(r, h.row);
@ -157,7 +155,6 @@ void RawFile::find_geometry() {
g.width = m_master.pixels_x(); g.width = m_master.pixels_x();
g.height = m_master.pixels_y(); g.height = m_master.pixels_y();
m_geometry.module_pixel_0.push_back(g); m_geometry.module_pixel_0.push_back(g);
} }
r++; r++;
@ -168,23 +165,20 @@ void RawFile::find_geometry() {
m_geometry.modules_x = c; m_geometry.modules_x = c;
m_geometry.modules_y = r; m_geometry.modules_y = r;
m_geometry.pixels_y += static_cast<size_t>((r - 1) * cfg.module_gap_row); m_geometry.pixels_y += static_cast<size_t>((r - 1) * cfg.module_gap_row);
} }
Frame RawFile::get_frame(size_t frame_index) { 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(); std::byte *frame_buffer = f.data();
get_frame_into(frame_index, frame_buffer); get_frame_into(frame_index, frame_buffer);
return f; return f;
} }
size_t RawFile::bytes_per_pixel() const { return m_master.bitdepth() / 8; }
size_t RawFile::bytes_per_pixel() const { void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer,
return m_master.bitdepth() / 8; 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 << ")"; LOG(logDEBUG) << "RawFile::get_frame_into(" << frame_index << ")";
if (frame_index >= total_frames()) { if (frame_index >= total_frames()) {
throw std::runtime_error(LOCATION + "Frame number out of range"); 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<size_t> frame_numbers(n_modules()); std::vector<size_t> frame_numbers(n_modules());
std::vector<size_t> frame_indices(n_modules(), frame_index); std::vector<size_t> frame_indices(n_modules(), frame_index);
// sync the frame numbers // 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) { 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 // 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] = 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]);
} }
} }
@ -228,11 +223,14 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect
auto corrected_idx = frame_indices[part_idx]; auto corrected_idx = frame_indices[part_idx];
// This is where we start writing // This is where we start writing
auto offset = (m_geometry.module_pixel_0[part_idx].origin_y * m_geometry.pixels_x + auto offset = (m_geometry.module_pixel_0[part_idx].origin_y *
m_geometry.module_pixel_0[part_idx].origin_x)*m_master.bitdepth()/8; 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) if (m_geometry.module_pixel_0[part_idx].origin_x != 0)
throw std::runtime_error(LOCATION + " Implementation error. x pos not 0."); throw std::runtime_error(LOCATION +
" Implementation error. x pos not 0.");
// TODO! What if the files don't match? // TODO! What if the files don't match?
m_subfiles[part_idx]->seek(corrected_idx); m_subfiles[part_idx]->seek(corrected_idx);
@ -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); auto dest = (irow * this->m_geometry.pixels_x + icol);
dest = dest * m_master.bitdepth() / 8; dest = dest * m_master.bitdepth() / 8;
memcpy(frame_buffer + dest, memcpy(frame_buffer + dest,
part_buffer + cur_row * pos.width * part_buffer +
m_master.bitdepth() / 8, cur_row * pos.width * m_master.bitdepth() / 8,
pos.width * m_master.bitdepth() / 8); pos.width * m_master.bitdepth() / 8);
} }
} }
delete[] part_buffer; delete[] part_buffer;
} }
} }
std::vector<Frame> RawFile::read_n(size_t n_frames) { std::vector<Frame> 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); return m_subfiles[0]->frame_number(frame_index);
} }
} // namespace aare } // namespace aare

View File

@ -1,18 +1,18 @@
#include "aare/RawFile.hpp"
#include "aare/File.hpp" #include "aare/File.hpp"
#include "aare/RawMasterFile.hpp" //needed for ROI #include "aare/RawMasterFile.hpp" //needed for ROI
#include "aare/RawFile.hpp"
#include <catch2/catch_test_macros.hpp> #include <catch2/catch_test_macros.hpp>
#include <filesystem> #include <filesystem>
#include "test_config.hpp" #include "test_config.hpp"
using aare::File; using aare::File;
TEST_CASE("Read number of frames from a jungfrau raw file", "[.integration]") { 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)); REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r"); 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]") { 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)); REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r"); 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]") { 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)); REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r"); File f(fpath, "r");
@ -49,8 +51,10 @@ TEST_CASE("Read a frame number too high throws", "[.integration]") {
REQUIRE_THROWS(f.frame_number(10)); REQUIRE_THROWS(f.frame_number(10));
} }
TEST_CASE("Read a frame numbers where the subfile is missing throws", "[.integration]") { TEST_CASE("Read a frame numbers where the subfile is missing throws",
auto fpath = test_data_path() / "jungfrau" / "jungfrau_missing_subfile_master_0.json"; "[.integration]") {
auto fpath = test_data_path() / "jungfrau" /
"jungfrau_missing_subfile_master_0.json";
REQUIRE(std::filesystem::exists(fpath)); REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r"); File f(fpath, "r");
@ -69,15 +73,18 @@ TEST_CASE("Read a frame numbers where the subfile is missing throws", "[.integra
REQUIRE_THROWS(f.frame_number(10)); REQUIRE_THROWS(f.frame_number(10));
} }
TEST_CASE("Read data from a jungfrau 500k single port raw file",
TEST_CASE("Read data from a jungfrau 500k single port raw file", "[.integration]") { "[.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)); REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r"); 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 // we know this file has 10 frames with pixel 0,0 being: 2123, 2051, 2109,
std::vector<uint16_t> pixel_0_0 = {2123, 2051, 2109, 2117, 2089, 2095, 2072, 2126, 2097, 2102}; // 2117, 2089, 2095, 2072, 2126, 2097, 2102
std::vector<uint16_t> pixel_0_0 = {2123, 2051, 2109, 2117, 2089,
2095, 2072, 2126, 2097, 2102};
for (size_t i = 0; i < 10; i++) { for (size_t i = 0; i < 10; i++) {
auto frame = f.read_frame(); auto frame = f.read_frame();
CHECK(frame.rows() == 512); 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]") { 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)); 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)); REQUIRE(std::filesystem::exists(fpath_npy));
File raw(fpath_raw, "r"); 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]") { 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)); REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r"); File f(fpath, "r");
// we know this file has 10 frames check read_multiport.py for the values // we know this file has 10 frames check read_multiport.py for the values
std::vector<uint16_t> pixel_0_0 = {2099, 2121, 2108, 2084, 2084, 2118, 2066, 2108, 2112, 2116}; std::vector<uint16_t> pixel_0_0 = {2099, 2121, 2108, 2084, 2084,
std::vector<uint16_t> pixel_0_1 = {2842, 2796, 2865, 2798, 2805, 2817, 2852, 2789, 2792, 2833}; 2118, 2066, 2108, 2112, 2116};
std::vector<uint16_t> pixel_255_1023 = {2149, 2037, 2115, 2102, 2118, 2090, 2036, 2071, 2073, 2142}; std::vector<uint16_t> pixel_0_1 = {2842, 2796, 2865, 2798, 2805,
std::vector<uint16_t> pixel_511_1023 = {3231, 3169, 3167, 3162, 3168, 3160, 3171, 3171, 3169, 3171}; 2817, 2852, 2789, 2792, 2833};
std::vector<uint16_t> pixel_1_0 = {2748, 2614, 2665, 2629, 2618, 2630, 2631, 2634, 2577, 2598}; std::vector<uint16_t> pixel_255_1023 = {2149, 2037, 2115, 2102, 2118,
2090, 2036, 2071, 2073, 2142};
std::vector<uint16_t> pixel_511_1023 = {3231, 3169, 3167, 3162, 3168,
3160, 3171, 3171, 3169, 3171};
std::vector<uint16_t> pixel_1_0 = {2748, 2614, 2665, 2629, 2618,
2630, 2631, 2634, 2577, 2598};
for (size_t i = 0; i < 10; i++) { for (size_t i = 0; i < 10; i++) {
auto frame = f.read_frame(); auto frame = f.read_frame();
@ -152,5 +167,3 @@ TEST_CASE("Read file with unordered frames", "[.integration]") {
File f(fpath); File f(fpath);
REQUIRE_THROWS((f.read_frame())); REQUIRE_THROWS((f.read_frame()));
} }

View File

@ -37,18 +37,15 @@ std::filesystem::path RawFileNameComponents::master_fname() const {
} }
std::filesystem::path RawFileNameComponents::data_fname(size_t mod_id, std::filesystem::path RawFileNameComponents::data_fname(size_t mod_id,
size_t file_id size_t file_id) const {
) const {
std::string fmt = "{}_d{}_f{}_{}.raw"; 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) { if (m_old_scheme) {
fmt = "{}_d{}_f{:012}_{}.raw"; fmt = "{}_d{}_f{:012}_{}.raw";
} }
return m_base_path / fmt::format(fmt, m_base_name, mod_id, return m_base_path /
file_id, m_file_index); fmt::format(fmt, m_base_name, mod_id, file_id, m_file_index);
} }
void RawFileNameComponents::set_old_scheme(bool old_scheme) { void RawFileNameComponents::set_old_scheme(bool old_scheme) {
@ -135,10 +132,8 @@ ScanParameters RawMasterFile::scan_parameters() const {
return m_scan_parameters; return m_scan_parameters;
} }
std::optional<ROI> RawMasterFile::roi() const { return m_roi; } std::optional<ROI> RawMasterFile::roi() const { return m_roi; }
void RawMasterFile::parse_json(const std::filesystem::path &fpath) { void RawMasterFile::parse_json(const std::filesystem::path &fpath) {
std::ifstream ifs(fpath); std::ifstream ifs(fpath);
json j; json j;
@ -177,7 +172,6 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) {
// keep the optional empty // keep the optional empty
} }
// ---------------------------------------------------------------- // ----------------------------------------------------------------
// Special treatment of analog flag because of Moench03 // Special treatment of analog flag because of Moench03
try { try {
@ -233,13 +227,13 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) {
std::string scan_parameters = j.at("Scan Parameters"); std::string scan_parameters = j.at("Scan Parameters");
m_scan_parameters = ScanParameters(scan_parameters); m_scan_parameters = ScanParameters(scan_parameters);
if (v < 7.21) { if (v < 7.21) {
m_scan_parameters.increment_stop(); //adjust for endpoint being included m_scan_parameters
.increment_stop(); // adjust for endpoint being included
} }
} catch (const json::out_of_range &e) { } catch (const json::out_of_range &e) {
// not a scan // not a scan
} }
try { try {
ROI tmp_roi; ROI tmp_roi;
auto obj = j.at("Receiver Roi"); auto obj = j.at("Receiver Roi");
@ -260,19 +254,14 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) {
m_roi = tmp_roi; m_roi = tmp_roi;
} }
} catch (const json::out_of_range &e) { } catch (const json::out_of_range &e) {
// leave the optional empty // leave the optional empty
} }
// if we have an roi we need to update the geometry for the subfiles // if we have an roi we need to update the geometry for the subfiles
if (m_roi) { if (m_roi) {
} }
// Update detector type for Moench // Update detector type for Moench
// TODO! How does this work with old .raw master files? // TODO! How does this work with old .raw master files?
#ifdef AARE_VERBOSE #ifdef AARE_VERBOSE
@ -382,7 +371,6 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) {
m_type = DetectorType::Moench03_old; m_type = DetectorType::Moench03_old;
} }
// TODO! Look for d0, d1...dn and update geometry // TODO! Look for d0, d1...dn and update geometry
if (m_geometry.col == 0 && m_geometry.row == 0) { if (m_geometry.col == 0 && m_geometry.row == 0) {
m_geometry = {1, 1}; m_geometry = {1, 1};

View File

@ -1,11 +1,10 @@
#include "aare/RawMasterFile.hpp" #include "aare/RawMasterFile.hpp"
#include <catch2/catch_test_macros.hpp>
#include "test_config.hpp" #include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
using namespace aare; using namespace aare;
TEST_CASE("Parse a master file fname") { TEST_CASE("Parse a master file fname") {
RawFileNameComponents m("test_master_1.json"); RawFileNameComponents m("test_master_1.json");
REQUIRE(m.base_name() == "test"); REQUIRE(m.base_name() == "test");
@ -47,10 +46,9 @@ TEST_CASE("Master file name does not fit pattern"){
REQUIRE_THROWS(RawFileNameComponents("test_master_1.txt")); REQUIRE_THROWS(RawFileNameComponents("test_master_1.txt"));
} }
TEST_CASE("Parse scan parameters") { TEST_CASE("Parse scan parameters") {
ScanParameters s("[enabled\ndac dac 4\nstart 500\nstop 2200\nstep 5\nsettleTime 100us\n]"); ScanParameters s("[enabled\ndac dac 4\nstart 500\nstop 2200\nstep "
"5\nsettleTime 100us\n]");
REQUIRE(s.enabled()); REQUIRE(s.enabled());
REQUIRE(s.dac() == "dac 4"); REQUIRE(s.dac() == "dac 4");
REQUIRE(s.start() == 500); REQUIRE(s.start() == 500);
@ -67,9 +65,9 @@ TEST_CASE("A disabled scan"){
REQUIRE(s.step() == 0); REQUIRE(s.step() == 0);
} }
TEST_CASE("Parse a master file in .json format", "[.integration]") { TEST_CASE("Parse a master file in .json format", "[.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)); REQUIRE(std::filesystem::exists(fpath));
RawMasterFile f(fpath); RawMasterFile f(fpath);
@ -103,7 +101,6 @@ TEST_CASE("Parse a master file in .json format", "[.integration]"){
// Jungfrau doesn't write but it is 16 // Jungfrau doesn't write but it is 16
REQUIRE(f.bitdepth() == 16); REQUIRE(f.bitdepth() == 16);
// "Frame Discard Policy": "nodiscard", // "Frame Discard Policy": "nodiscard",
// "Frame Padding": 1, // "Frame Padding": 1,
@ -146,12 +143,14 @@ TEST_CASE("Parse a master file in .json format", "[.integration]"){
REQUIRE_FALSE(f.analog_samples()); 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]") { 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"; auto fpath =
test_data_path() /
"moench/"
"moench04_noise_200V_sto_both_100us_no_light_thresh_900_master_0.raw";
REQUIRE(std::filesystem::exists(fpath)); REQUIRE(std::filesystem::exists(fpath));
RawMasterFile f(fpath); RawMasterFile f(fpath);
@ -209,11 +208,8 @@ TEST_CASE("Parse a master file in .raw format", "[.integration]"){
// Detector Type : 1 byte // Detector Type : 1 byte
// Header Version : 1 byte // Header Version : 1 byte
// Packets Caught Mask : 64 bytes // Packets Caught Mask : 64 bytes
} }
TEST_CASE("Read eiger master file", "[.integration]") { TEST_CASE("Read eiger master file", "[.integration]") {
auto fpath = test_data_path() / "eiger" / "eiger_500k_32bit_master_0.json"; auto fpath = test_data_path() / "eiger" / "eiger_500k_32bit_master_0.json";
REQUIRE(std::filesystem::exists(fpath)); REQUIRE(std::filesystem::exists(fpath));
@ -282,7 +278,4 @@ REQUIRE(f.frame_padding() == 1);
// "Packets Caught Mask": "64 bytes" // "Packets Caught Mask": "64 bytes"
// } // }
// } // }
} }

View File

@ -1,27 +1,22 @@
#include "aare/RawSubFile.hpp" #include "aare/RawSubFile.hpp"
#include "aare/PixelMap.hpp" #include "aare/PixelMap.hpp"
#include "aare/algorithm.hpp" #include "aare/algorithm.hpp"
#include "aare/utils/ifstream_helpers.hpp"
#include "aare/logger.hpp" #include "aare/logger.hpp"
#include "aare/utils/ifstream_helpers.hpp"
#include <cstring> // memcpy #include <cstring> // memcpy
#include <fmt/core.h> #include <fmt/core.h>
#include <iostream> #include <iostream>
#include <regex> #include <regex>
namespace aare { namespace aare {
RawSubFile::RawSubFile(const std::filesystem::path &fname, RawSubFile::RawSubFile(const std::filesystem::path &fname,
DetectorType detector, size_t rows, size_t cols, DetectorType detector, size_t rows, size_t cols,
size_t bitdepth, uint32_t pos_row, uint32_t pos_col) size_t bitdepth, uint32_t pos_row, uint32_t pos_col)
: m_detector_type(detector), m_bitdepth(bitdepth), : m_detector_type(detector), m_bitdepth(bitdepth), m_rows(rows),
m_rows(rows), m_cols(cols), m_cols(cols), m_bytes_per_frame((m_bitdepth / 8) * m_rows * m_cols),
m_bytes_per_frame((m_bitdepth / 8) * m_rows * m_cols), m_pos_row(pos_row), m_pos_row(pos_row), m_pos_col(pos_col) {
m_pos_col(pos_col) {
LOG(logDEBUG) << "RawSubFile::RawSubFile()"; LOG(logDEBUG) << "RawSubFile::RawSubFile()";
if (m_detector_type == DetectorType::Moench03_old) { if (m_detector_type == DetectorType::Moench03_old) {
@ -30,7 +25,6 @@ RawSubFile::RawSubFile(const std::filesystem::path &fname,
m_pixel_map = GenerateEigerFlipRowsPixelMap(); m_pixel_map = GenerateEigerFlipRowsPixelMap();
} }
parse_fname(fname); parse_fname(fname);
scan_files(); scan_files();
open_file(m_current_file_index); // open the first file 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) auto frame_offset = (file_index)
? frame_index - m_last_frame_in_file[file_index - 1] ? frame_index - m_last_frame_in_file[file_index - 1]
: frame_index; : 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); m_file.seekg(byte_offset);
} }
@ -85,7 +80,8 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) {
} else if (m_bitdepth == 32) { } else if (m_bitdepth == 32) {
read_with_map<uint32_t>(image_buf); read_with_map<uint32_t>(image_buf);
} else { } else {
throw std::runtime_error("Unsupported bitdepth for read with pixel map"); throw std::runtime_error(
"Unsupported bitdepth for read with pixel map");
} }
} else { } else {
@ -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++) { for (size_t i = 0; i < n_frames; i++) {
read_into(image_buf, header); read_into(image_buf, header);
image_buf += bytes_per_frame(); image_buf += bytes_per_frame();
@ -115,10 +112,7 @@ void RawSubFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader
} }
} }
template <typename T> void RawSubFile::read_with_map(std::byte *image_buf) {
template <typename T>
void RawSubFile::read_with_map(std::byte *image_buf) {
auto part_buffer = new std::byte[bytes_per_frame()]; auto part_buffer = new std::byte[bytes_per_frame()];
m_file.read(reinterpret_cast<char *>(part_buffer), bytes_per_frame()); m_file.read(reinterpret_cast<char *>(part_buffer), bytes_per_frame());
auto *data = reinterpret_cast<T *>(image_buf); auto *data = reinterpret_cast<T *>(image_buf);
@ -157,14 +151,17 @@ void RawSubFile::parse_fname(const std::filesystem::path &fname) {
std::smatch match; std::smatch match;
if (std::regex_match(m_base_name, match, pattern)) { 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_offset = std::stoi(match[4].str()); // find the first file index in
m_base_name = match[1].str() + match[2].str() + match[3].str() + "{}" + match[5].str(); // 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) << "Base name: " << m_base_name;
LOG(logDEBUG) << "Offset: " << m_offset; LOG(logDEBUG) << "Offset: " << m_offset;
LOG(logDEBUG) << "Path: " << m_path.string(); LOG(logDEBUG) << "Path: " << m_path.string();
} else { } else {
throw std::runtime_error( throw std::runtime_error(
LOCATION + fmt::format("Could not parse file name {}", fname.string())); LOCATION +
fmt::format("Could not parse file name {}", fname.string()));
} }
} }
@ -180,7 +177,8 @@ void RawSubFile::open_file(size_t file_index) {
m_file.open(fname, std::ios::binary); m_file.open(fname, std::ios::binary);
if (!m_file.is_open()) { if (!m_file.is_open()) {
throw std::runtime_error( 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; m_current_file_index = file_index;
} }
@ -195,7 +193,8 @@ void RawSubFile::scan_files() {
auto n_frames = std::filesystem::file_size(fpath(file_index)) / auto n_frames = std::filesystem::file_size(fpath(file_index)) /
(m_bytes_per_frame + sizeof(DetectorHeader)); (m_bytes_per_frame + sizeof(DetectorHeader));
m_last_frame_in_file.push_back(n_frames); 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; ++file_index;
} }

View File

@ -1,13 +1,14 @@
#include "aare/RawSubFile.hpp" #include "aare/RawSubFile.hpp"
#include "aare/File.hpp" #include "aare/File.hpp"
#include "aare/NDArray.hpp" #include "aare/NDArray.hpp"
#include <catch2/catch_test_macros.hpp>
#include "test_config.hpp" #include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
using namespace aare; using namespace aare;
TEST_CASE("Read frames directly from a RawSubFile", "[.files]") { TEST_CASE("Read frames directly from a RawSubFile", "[.files]") {
auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f0_0.raw"; auto fpath_raw =
test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f0_0.raw";
REQUIRE(std::filesystem::exists(fpath_raw)); REQUIRE(std::filesystem::exists(fpath_raw));
RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16); RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16);
@ -17,8 +18,8 @@ TEST_CASE("Read frames directly from a RawSubFile", "[.files]"){
REQUIRE(f.bytes_per_frame() == 512 * 1024 * 2); REQUIRE(f.bytes_per_frame() == 512 * 1024 * 2);
REQUIRE(f.bytes_per_pixel() == 2); REQUIRE(f.bytes_per_pixel() == 2);
auto fpath_npy =
auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy"; test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
REQUIRE(std::filesystem::exists(fpath_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
@ -27,9 +28,9 @@ TEST_CASE("Read frames directly from a RawSubFile", "[.files]"){
CHECK(f.frames_in_file() == 10); CHECK(f.frames_in_file() == 10);
CHECK(npy.total_frames() == 10); CHECK(npy.total_frames() == 10);
DetectorHeader header{}; DetectorHeader header{};
NDArray<uint16_t, 2> image({static_cast<ssize_t>(f.rows()), static_cast<ssize_t>(f.cols())}); NDArray<uint16_t, 2> image(
{static_cast<ssize_t>(f.rows()), static_cast<ssize_t>(f.cols())});
for (size_t i = 0; i < 10; ++i) { for (size_t i = 0; i < 10; ++i) {
CHECK(f.tell() == i); CHECK(f.tell() == i);
f.read_into(image.buffer(), &header); f.read_into(image.buffer(), &header);
@ -38,20 +39,22 @@ 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 // we know this file has 10 frames with frame numbers 1 to 10
// f0 1,2,3 // f0 1,2,3
// f1 4,5,6 <-- starting here // f1 4,5,6 <-- starting here
// f2 7,8,9 // f2 7,8,9
// f3 10 // 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)); REQUIRE(std::filesystem::exists(fpath_raw));
RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16); RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16);
auto fpath_npy =
auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy"; test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
REQUIRE(std::filesystem::exists(fpath_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
@ -61,9 +64,9 @@ TEST_CASE("Read frames directly from a RawSubFile starting at the second file",
CHECK(f.frames_in_file() == 7); CHECK(f.frames_in_file() == 7);
CHECK(npy.total_frames() == 10); CHECK(npy.total_frames() == 10);
DetectorHeader header{}; DetectorHeader header{};
NDArray<uint16_t, 2> image({static_cast<ssize_t>(f.rows()), static_cast<ssize_t>(f.cols())}); NDArray<uint16_t, 2> image(
{static_cast<ssize_t>(f.rows()), static_cast<ssize_t>(f.cols())});
for (size_t i = 0; i < 7; ++i) { for (size_t i = 0; i < 7; ++i) {
CHECK(f.tell() == i); CHECK(f.tell() == i);
f.read_into(image.buffer(), &header); f.read_into(image.buffer(), &header);

View File

@ -161,12 +161,10 @@ TEST_CASE("cumsum works with negative numbers", "[algorithm]") {
REQUIRE(result[4] == -10); REQUIRE(result[4] == -10);
} }
TEST_CASE("cumsum on an empty vector", "[algorithm]") { TEST_CASE("cumsum on an empty vector", "[algorithm]") {
std::vector<double> vec = {}; std::vector<double> vec = {};
auto result = aare::cumsum(vec); auto result = aare::cumsum(vec);
REQUIRE(result.size() == 0); REQUIRE(result.size() == 0);
} }
TEST_CASE("All equal on an empty vector is false", "[algorithm]") { 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); 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<int> vec = {1, 2}; std::vector<int> vec = {1, 2};
REQUIRE(aare::all_equal(vec) == false); REQUIRE(aare::all_equal(vec) == false);
} }

View File

@ -21,9 +21,11 @@ uint16_t adc_sar_05_decode64to16(uint64_t input){
return output; return output;
} }
void adc_sar_05_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output){ void adc_sar_05_decode64to16(NDView<uint64_t, 2> input,
NDView<uint16_t, 2> output) {
if (input.shape() != output.shape()) { if (input.shape() != output.shape()) {
throw std::invalid_argument(LOCATION + " input and output shapes must match"); throw std::invalid_argument(LOCATION +
" input and output shapes must match");
} }
for (ssize_t i = 0; i < input.shape(0); i++) { for (ssize_t i = 0; i < input.shape(0); i++) {
@ -52,9 +54,11 @@ uint16_t adc_sar_04_decode64to16(uint64_t input){
return output; return output;
} }
void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output){ void adc_sar_04_decode64to16(NDView<uint64_t, 2> input,
NDView<uint16_t, 2> output) {
if (input.shape() != output.shape()) { if (input.shape() != output.shape()) {
throw std::invalid_argument(LOCATION + " input and output shapes must match"); throw std::invalid_argument(LOCATION +
" input and output shapes must match");
} }
for (ssize_t i = 0; i < input.shape(0); i++) { for (ssize_t i = 0; i < input.shape(0); i++) {
for (ssize_t j = 0; j < input.shape(1); j++) { for (ssize_t j = 0; j < input.shape(1); j++) {
@ -65,7 +69,8 @@ void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> outpu
double apply_custom_weights(uint16_t input, const NDView<double, 1> weights) { double apply_custom_weights(uint16_t input, const NDView<double, 1> weights) {
if (weights.size() > 16) { if (weights.size() > 16) {
throw std::invalid_argument("weights size must be less than or equal to 16"); throw std::invalid_argument(
"weights size must be less than or equal to 16");
} }
double result = 0.0; double result = 0.0;
@ -73,12 +78,13 @@ double apply_custom_weights(uint16_t input, const NDView<double, 1> weights) {
result += ((input >> i) & 1) * std::pow(weights[i], i); result += ((input >> i) & 1) * std::pow(weights[i], i);
} }
return result; return result;
} }
void apply_custom_weights(NDView<uint16_t, 1> input, NDView<double, 1> output, const NDView<double,1> weights) { void apply_custom_weights(NDView<uint16_t, 1> input, NDView<double, 1> output,
const NDView<double, 1> weights) {
if (input.shape() != output.shape()) { if (input.shape() != output.shape()) {
throw std::invalid_argument(LOCATION + " input and output shapes must match"); 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
@ -90,13 +96,12 @@ void apply_custom_weights(NDView<uint16_t, 1> input, NDView<double, 1> output, c
// Apply custom weights to each element in the input array // 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; 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]; result += ((input(i) >> bit_index) & 1) * weights_powers[bit_index];
} }
output(i) = result; output(i) = result;
} }
} }
} // namespace aare } // namespace aare

View File

@ -1,8 +1,8 @@
#include "aare/decode.hpp" #include "aare/decode.hpp"
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/catch_test_macros.hpp>
#include "aare/NDArray.hpp" #include "aare/NDArray.hpp"
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
using Catch::Matchers::WithinAbs; using Catch::Matchers::WithinAbs;
#include <vector> #include <vector>
@ -11,7 +11,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){
uint16_t output = aare::adc_sar_05_decode64to16(input); uint16_t output = aare::adc_sar_05_decode64to16(input);
CHECK(output == 0); CHECK(output == 0);
// bit 29 on th input is bit 0 on the output // bit 29 on th input is bit 0 on the output
input = 1UL << 29; input = 1UL << 29;
output = aare::adc_sar_05_decode64to16(input); output = aare::adc_sar_05_decode64to16(input);
@ -25,7 +24,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){
CHECK(output == (1 << i)); CHECK(output == (1 << i));
} }
// test a few "random" values // test a few "random" values
input = 0; input = 0;
input |= (1UL << 29); input |= (1UL << 29);
@ -34,7 +32,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){
output = aare::adc_sar_05_decode64to16(input); output = aare::adc_sar_05_decode64to16(input);
CHECK(output == 7UL); CHECK(output == 7UL);
input = 0; input = 0;
input |= (1UL << 18); input |= (1UL << 18);
input |= (1UL << 27); input |= (1UL << 27);
@ -49,7 +46,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){
CHECK(output == 3072UL); CHECK(output == 3072UL);
} }
TEST_CASE("test_apply_custom_weights") { TEST_CASE("test_apply_custom_weights") {
uint16_t input = 1; uint16_t input = 1;
@ -60,7 +56,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){
auto weights = weights_data.view(); auto weights = weights_data.view();
double output = aare::apply_custom_weights(input, weights); double output = aare::apply_custom_weights(input, weights);
CHECK_THAT(output, WithinAbs(1.0, 0.001)); 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); output = aare::apply_custom_weights(input, weights);
CHECK_THAT(output, WithinAbs(2.1, 0.001)); CHECK_THAT(output, WithinAbs(2.1, 0.001));
input = 1 << 2; input = 1 << 2;
output = aare::apply_custom_weights(input, weights); output = aare::apply_custom_weights(input, weights);
CHECK_THAT(output, WithinAbs(3.24, 0.001)); CHECK_THAT(output, WithinAbs(3.24, 0.001));
@ -76,5 +70,4 @@ TEST_CASE("test_adc_sar_05_decode64to16"){
input = 0b111; input = 0b111;
output = aare::apply_custom_weights(input, weights); output = aare::apply_custom_weights(input, weights);
CHECK_THAT(output, WithinAbs(6.34, 0.001)); CHECK_THAT(output, WithinAbs(6.34, 0.001));
} }

View File

@ -6,15 +6,11 @@
#include <fmt/core.h> #include <fmt/core.h>
namespace aare { namespace aare {
void assert_failed(const std::string &msg) {
void assert_failed(const std::string &msg)
{
fmt::print(msg); fmt::print(msg);
exit(1); exit(1);
} }
/** /**
* @brief Convert a DetectorType to a string * @brief Convert a DetectorType to a string
* @param type DetectorType * @param type DetectorType
@ -121,7 +117,8 @@ template <> TimingMode StringTo(const std::string &arg) {
return TimingMode::Auto; return TimingMode::Auto;
if (arg == "trigger") if (arg == "trigger")
return TimingMode::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 +
"\"");
} }
/** /**
@ -151,7 +148,8 @@ template <> FrameDiscardPolicy StringTo(const std::string &arg) {
return FrameDiscardPolicy::Discard; return FrameDiscardPolicy::Discard;
if (arg == "discardpartial") if (arg == "discardpartial")
return FrameDiscardPolicy::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 + "\"");
} }
/** /**

View File

@ -3,12 +3,12 @@
#include <catch2/catch_test_macros.hpp> #include <catch2/catch_test_macros.hpp>
#include <string> #include <string>
using aare::ToString;
using aare::StringTo; using aare::StringTo;
using aare::ToString;
TEST_CASE("Enum to string conversion") { 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 // TODO! By the way I don't think the enum string conversions should be in
// but let's use this to show a test // the defs.hpp file but let's use this to show a test
REQUIRE(ToString(aare::DetectorType::Generic) == "Generic"); REQUIRE(ToString(aare::DetectorType::Generic) == "Generic");
REQUIRE(ToString(aare::DetectorType::Eiger) == "Eiger"); REQUIRE(ToString(aare::DetectorType::Eiger) == "Eiger");
REQUIRE(ToString(aare::DetectorType::Gotthard) == "Gotthard"); REQUIRE(ToString(aare::DetectorType::Gotthard) == "Gotthard");
@ -17,25 +17,37 @@ TEST_CASE("Enum to string conversion") {
REQUIRE(ToString(aare::DetectorType::Moench) == "Moench"); REQUIRE(ToString(aare::DetectorType::Moench) == "Moench");
REQUIRE(ToString(aare::DetectorType::Mythen3) == "Mythen3"); REQUIRE(ToString(aare::DetectorType::Mythen3) == "Mythen3");
REQUIRE(ToString(aare::DetectorType::Gotthard2) == "Gotthard2"); 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) == "Moench03");
REQUIRE(ToString(aare::DetectorType::Moench03_old) == "Moench03_old"); REQUIRE(ToString(aare::DetectorType::Moench03_old) == "Moench03_old");
REQUIRE(ToString(aare::DetectorType::Unknown) == "Unknown"); REQUIRE(ToString(aare::DetectorType::Unknown) == "Unknown");
} }
TEST_CASE("String to enum") { TEST_CASE("String to enum") {
REQUIRE(StringTo<aare::DetectorType>("Generic") == aare::DetectorType::Generic); REQUIRE(StringTo<aare::DetectorType>("Generic") ==
aare::DetectorType::Generic);
REQUIRE(StringTo<aare::DetectorType>("Eiger") == aare::DetectorType::Eiger); REQUIRE(StringTo<aare::DetectorType>("Eiger") == aare::DetectorType::Eiger);
REQUIRE(StringTo<aare::DetectorType>("Gotthard") == aare::DetectorType::Gotthard); REQUIRE(StringTo<aare::DetectorType>("Gotthard") ==
REQUIRE(StringTo<aare::DetectorType>("Jungfrau") == aare::DetectorType::Jungfrau); aare::DetectorType::Gotthard);
REQUIRE(StringTo<aare::DetectorType>("ChipTestBoard") == aare::DetectorType::ChipTestBoard); REQUIRE(StringTo<aare::DetectorType>("Jungfrau") ==
REQUIRE(StringTo<aare::DetectorType>("Moench") == aare::DetectorType::Moench); aare::DetectorType::Jungfrau);
REQUIRE(StringTo<aare::DetectorType>("Mythen3") == aare::DetectorType::Mythen3); REQUIRE(StringTo<aare::DetectorType>("ChipTestBoard") ==
REQUIRE(StringTo<aare::DetectorType>("Gotthard2") == aare::DetectorType::Gotthard2); aare::DetectorType::ChipTestBoard);
REQUIRE(StringTo<aare::DetectorType>("Xilinx_ChipTestBoard") == aare::DetectorType::Xilinx_ChipTestBoard); REQUIRE(StringTo<aare::DetectorType>("Moench") ==
REQUIRE(StringTo<aare::DetectorType>("Moench03") == aare::DetectorType::Moench03); aare::DetectorType::Moench);
REQUIRE(StringTo<aare::DetectorType>("Moench03_old") == aare::DetectorType::Moench03_old); REQUIRE(StringTo<aare::DetectorType>("Mythen3") ==
REQUIRE(StringTo<aare::DetectorType>("Unknown") == aare::DetectorType::Unknown); aare::DetectorType::Mythen3);
REQUIRE(StringTo<aare::DetectorType>("Gotthard2") ==
aare::DetectorType::Gotthard2);
REQUIRE(StringTo<aare::DetectorType>("Xilinx_ChipTestBoard") ==
aare::DetectorType::Xilinx_ChipTestBoard);
REQUIRE(StringTo<aare::DetectorType>("Moench03") ==
aare::DetectorType::Moench03);
REQUIRE(StringTo<aare::DetectorType>("Moench03_old") ==
aare::DetectorType::Moench03_old);
REQUIRE(StringTo<aare::DetectorType>("Unknown") ==
aare::DetectorType::Unknown);
} }
TEST_CASE("Enum values") { TEST_CASE("Enum values") {
@ -85,5 +97,6 @@ TEST_CASE("DynamicCluster creation") {
// double v3 = c2.get<double>(33 * 44 - 1); // double v3 = c2.get<double>(33 * 44 - 1);
// REQUIRE(aare::compare_floats<double>(123.11, v3)); // REQUIRE(aare::compare_floats<double>(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
// } // }

View File

@ -8,8 +8,9 @@ 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: {} {} {} {}", fmt::println("update_geometry_with_roi() called with ROI: {} {} {} {}",
roi.xmin, roi.xmax, roi.ymin, roi.ymax); roi.xmin, roi.xmax, roi.ymin, roi.ymax);
fmt::println("Geometry: {} {} {} {} {} {}", fmt::println("Geometry: {} {} {} {} {} {}", geo.modules_x, geo.modules_y,
geo.modules_x, geo.modules_y, geo.pixels_x, geo.pixels_y, geo.module_gap_row, geo.module_gap_col); geo.pixels_x, geo.pixels_y, geo.module_gap_row,
geo.module_gap_col);
#endif #endif
int pos_y = 0; int pos_y = 0;
int pos_y_increment = 0; int pos_y_increment = 0;
@ -41,9 +42,9 @@ DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) {
if (m.origin_y + m.height < roi.ymin) { if (m.origin_y + m.height < roi.ymin) {
m.height = 0; m.height = 0;
} else { } 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; m.height -= roi.ymin - m.origin_y;
} }
if (roi.ymax < m.origin_y + original_height) { if (roi.ymax < m.origin_y + original_height) {
m.height -= m.origin_y + original_height - roi.ymax; m.height -= m.origin_y + original_height - roi.ymax;
@ -52,7 +53,8 @@ DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) {
pos_y_increment = m.height; pos_y_increment = m.height;
} }
#ifdef AARE_VERBOSE #ifdef AARE_VERBOSE
fmt::println("Module {} {} {} {}", m.origin_x, m.origin_y, m.width, m.height); fmt::println("Module {} {} {} {}", m.origin_x, m.origin_y, m.width,
m.height);
#endif #endif
} }
// increment pos_y // increment pos_y
@ -65,7 +67,6 @@ DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) {
geo.pixels_y = roi.height(); geo.pixels_y = roi.height();
return geo; return geo;
} }
} // namespace aare } // namespace aare

View File

@ -1,6 +1,6 @@
#include "aare/File.hpp" #include "aare/File.hpp"
#include "aare/RawMasterFile.hpp" //needed for ROI
#include "aare/RawFile.hpp" #include "aare/RawFile.hpp"
#include "aare/RawMasterFile.hpp" //needed for ROI
#include <catch2/catch_test_macros.hpp> #include <catch2/catch_test_macros.hpp>
#include <filesystem> #include <filesystem>
@ -9,18 +9,16 @@
#include "test_config.hpp" #include "test_config.hpp"
TEST_CASE("Simple ROIs on one module") { TEST_CASE("Simple ROIs on one module") {
// DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI
// roi)
aare::DetectorGeometry geo; aare::DetectorGeometry geo;
aare::ModuleGeometry mod; aare::ModuleGeometry mod;
mod.origin_x = 0; mod.origin_x = 0;
mod.origin_y = 0; mod.origin_y = 0;
mod.width = 1024; mod.width = 1024;
mod.height = 512; mod.height = 512;
geo.pixels_x = 1024; geo.pixels_x = 1024;
geo.pixels_y = 512; geo.pixels_y = 512;
geo.modules_x = 1; geo.modules_x = 1;
@ -90,13 +88,11 @@ TEST_CASE("Simple ROIs on one module"){
} }
} }
TEST_CASE("Two modules side by side") { TEST_CASE("Two modules side by side") {
// DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI
// roi)
aare::DetectorGeometry geo; aare::DetectorGeometry geo;
aare::ModuleGeometry mod; aare::ModuleGeometry mod;
mod.origin_x = 0; mod.origin_x = 0;
mod.origin_y = 0; mod.origin_y = 0;
@ -145,7 +141,8 @@ TEST_CASE("Two modules side by side"){
} }
TEST_CASE("Three modules side by side") { TEST_CASE("Three modules side by side") {
// DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI
// roi)
aare::DetectorGeometry geo; aare::DetectorGeometry geo;
aare::ROI roi; aare::ROI roi;
roi.xmin = 700; roi.xmin = 700;
@ -185,7 +182,8 @@ TEST_CASE("Three modules side by side"){
} }
TEST_CASE("Four modules as a square") { TEST_CASE("Four modules as a square") {
// DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI
// roi)
aare::DetectorGeometry geo; aare::DetectorGeometry geo;
aare::ROI roi; aare::ROI roi;
roi.xmin = 500; roi.xmin = 500;

View File

@ -1,8 +1,7 @@
#include "aare/utils/task.hpp" #include "aare/utils/task.hpp"
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/catch_test_macros.hpp> #include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
TEST_CASE("Split a range into multiple tasks") { TEST_CASE("Split a range into multiple tasks") {
@ -26,7 +25,4 @@ TEST_CASE("Split a range into multiple tasks"){
REQUIRE(tasks[i].first == i); REQUIRE(tasks[i].first == i);
REQUIRE(tasks[i].second == i + 1); REQUIRE(tasks[i].second == i + 1);
} }
} }

View File

@ -2,8 +2,8 @@
#include <catch2/catch_test_macros.hpp> #include <catch2/catch_test_macros.hpp>
#include <climits> #include <climits>
#include <filesystem> #include <filesystem>
#include <fstream>
#include <fmt/core.h> #include <fmt/core.h>
#include <fstream>
TEST_CASE("Test suite can find data assets", "[.integration]") { TEST_CASE("Test suite can find data assets", "[.integration]") {
auto fpath = test_data_path() / "numpy" / "test_numpy_file.npy"; auto fpath = test_data_path() / "numpy" / "test_numpy_file.npy";