detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2025-06-13 15:57:14 +02:00
Code Issues Actions Packages Releases Activity

formatted main branch

Browse Source
This commit is contained in:
mazzol_a
2025-06-10 16:09:06 +02:00
parent efd2338f54
commit f9751902a2
87 changed files with 1710 additions and 1639 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
benchmarks/calculateeta_benchmark.cpp
View File

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

6
benchmarks/ndarray_benchmark.cpp
View File

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

15
include/aare/ArrayExpr.hpp
View File

@ -1,10 +1,9 @@
#pragma once
#include <cstdint>
#include <cstddef>
#include "aare/defs.hpp"
#include <array>
#include <cassert>
#include "aare/defs.hpp"
#include <cstddef>
#include <cstdint>
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 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>
@ -74,8 +75,6 @@ class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
};
template <typename A, typename B, ssize_t Ndim>
auto operator+(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &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);
}
} // namespace aare

6
include/aare/CircularFifo.hpp
View File

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

4
include/aare/ClusterFile.hpp
View File

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

3
include/aare/ClusterFinder.hpp
View File

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

6
include/aare/ClusterFinderMT.hpp
View File

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

8
include/aare/CtbRawFile.hpp
View File

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

51
include/aare/Dtype.hpp
View File

@ -6,31 +6,37 @@
namespace aare {
// The format descriptor is a single character that specifies the type of the data
// The format descriptor is a single character that specifies the type of the
// data
// - python documentation: https://docs.python.org/3/c-api/arg.html#numbers
// - py::format_descriptor<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.
// - numpy also doesn't use the same format. and also numpy associates the format
// with variable bitdepth types. (e.g. long is int64 on linux64 and int32 on win64)
// https://numpy.org/doc/stable/reference/arrays.scalars.html
// - numpy also doesn't use the same format. and also numpy associates the
// format
// with variable bitdepth types. (e.g. long is int64 on linux64 and int32 on
// win64) https://numpy.org/doc/stable/reference/arrays.scalars.html
//
// github issue discussing this:
// https://github.com/pybind/pybind11/issues/1908#issuecomment-658358767
//
// [IN LINUX] the difference is for int64 (long) and uint64 (unsigned long). The format
// descriptor is 'q' and 'Q' respectively and in the documentation it is 'l' and 'k'.
// [IN LINUX] the difference is for int64 (long) and uint64 (unsigned long). The
// format descriptor is 'q' and 'Q' respectively and in the documentation it is
// 'l' and 'k'.
// in practice numpy doesn't seem to care when reading buffer info: the library
// interprets 'q' or 'l' as int64 and 'Q' or 'L' as uint64.
// for this reason we decided to use the same format descriptor as pybind to avoid
// any further discrepancies.
// for this reason we decided to use the same format descriptor as pybind to
// avoid any further discrepancies.
// in the following order:
// int8, uint8, int16, uint16, int32, uint32, int64, uint64, float, double
const char DTYPE_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', 'I', 'q', 'Q', 'f', 'd'};
const char DTYPE_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i',
'I', 'q', 'Q', 'f', 'd'};
// on linux64 & apple
const char NUMPY_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', 'I', 'l', 'L', 'f', 'd'};
const char NUMPY_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i',
'I', 'l', 'L', 'f', 'd'};
/**
* @brief enum class to define the endianess of the system
*/
@ -52,12 +58,29 @@ enum class endian {
*/
class Dtype {
public:
enum TypeIndex { INT8, UINT8, INT16, UINT16, INT32, UINT32, INT64, UINT64, FLOAT, DOUBLE, ERROR, NONE };
enum TypeIndex {
INT8,
UINT8,
INT16,
UINT16,
INT32,
UINT32,
INT64,
UINT64,
FLOAT,
DOUBLE,
ERROR,
NONE
};
uint8_t bitdepth() const;
size_t bytes() const;
std::string format_descr() const { return std::string(1, DTYPE_FORMAT_DSC[static_cast<int>(m_type)]); }
std::string numpy_descr() const { return std::string(1, NUMPY_FORMAT_DSC[static_cast<int>(m_type)]); }
std::string format_descr() const {
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(std::string_view sv);

30
include/aare/File.hpp
View File

@ -6,11 +6,11 @@ namespace aare {
/**
* @brief RAII File class for reading, and in the future potentially writing
* image files in various formats. Minimal generic interface. For specail fuctions
* plase use the RawFile or NumpyFile classes directly.
* Wraps FileInterface to abstract the underlying file format
* @note **frame_number** refers the the frame number sent by the detector while **frame_index**
* is the position of the frame in the file
* image files in various formats. Minimal generic interface. For specail
* fuctions plase use the RawFile or NumpyFile classes directly. Wraps
* FileInterface to abstract the underlying file format
* @note **frame_number** refers the the frame number sent by the detector while
* **frame_index** is the position of the frame in the file
*/
class File {
std::unique_ptr<FileInterface> file_impl;
@ -25,9 +25,11 @@ class File {
* @throws std::invalid_argument if the file mode is not supported
*
*/
File(const std::filesystem::path &fname, const std::string &mode="r", const FileConfig &cfg = {});
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;
/**The same goes for copy assignment */
@ -39,15 +41,19 @@ class File {
// void close(); //!< close the file
Frame read_frame(); //!< read one frame from the file at the current position
Frame read_frame(size_t frame_index); //!< read one frame at the position given by frame number
std::vector<Frame> read_n(size_t n_frames); //!< read n_frames from the file at the current position
Frame
read_frame(); //!< read one frame from the file at the current position
Frame read_frame(size_t frame_index); //!< read one frame at the position
//!< given by frame number
std::vector<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, size_t n_frames);
size_t frame_number(); //!< get the frame number at the current position
size_t frame_number(size_t frame_index); //!< get the frame number at the given frame index
size_t frame_number(
size_t frame_index); //!< get the frame number at the given frame index
size_t bytes_per_frame() const;
size_t pixels_per_frame() const;
size_t bytes_per_pixel() const;
@ -59,8 +65,6 @@ class File {
size_t cols() const;
DetectorType detector_type() const;
};
} // namespace aare

26
include/aare/FileInterface.hpp
View File

@ -20,8 +20,10 @@ struct FileConfig {
uint64_t rows{};
uint64_t cols{};
bool operator==(const FileConfig &other) const {
return dtype == other.dtype && rows == other.rows && cols == other.cols && geometry == other.geometry &&
detector_type == other.detector_type && max_frames_per_file == other.max_frames_per_file;
return dtype == other.dtype && rows == other.rows &&
cols == other.cols && geometry == other.geometry &&
detector_type == other.detector_type &&
max_frames_per_file == other.max_frames_per_file;
}
bool operator!=(const FileConfig &other) const { return !(*this == other); }
@ -32,8 +34,11 @@ struct FileConfig {
int max_frames_per_file{};
size_t total_frames{};
std::string to_string() const {
return "{ dtype: " + dtype.to_string() + ", rows: " + std::to_string(rows) + ", cols: " + std::to_string(cols) +
", geometry: " + geometry.to_string() + ", detector_type: " + ToString(detector_type) +
return "{ dtype: " + dtype.to_string() +
", rows: " + std::to_string(rows) +
", cols: " + std::to_string(cols) +
", geometry: " + geometry.to_string() +
", detector_type: " + ToString(detector_type) +
", max_frames_per_file: " + std::to_string(max_frames_per_file) +
", total_frames: " + std::to_string(total_frames) + " }";
}
@ -42,7 +47,8 @@ struct FileConfig {
/**
* @brief FileInterface class to define the interface for file operations
* @note parent class for NumpyFile and RawFile
* @note all functions are pure virtual and must be implemented by the derived classes
* @note all functions are pure virtual and must be implemented by the derived
* classes
*/
class FileInterface {
public:
@ -64,17 +70,20 @@ class FileInterface {
* @param n_frames number of frames to read
* @return vector of frames
*/
virtual std::vector<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
* @return void
*/
virtual void read_into(std::byte *image_buf) = 0;
/**
* @brief read n_frames from the file at the current position and store them in the provided buffer
* @brief read n_frames from the file at the current position and store them
* in the provided buffer
* @param image_buf buffer to store the frames
* @param n_frames number of frames to read
* @return void
@ -134,7 +143,6 @@ class FileInterface {
*/
virtual size_t bitdepth() const = 0;
virtual DetectorType detector_type() const = 0;
// function to query the data type of the file

55
include/aare/Fit.hpp
View File

@ -23,16 +23,19 @@ NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par);
} // namespace func
/**
* @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> scurve2_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> scurve2_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
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);
/**
* @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, 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,
int n_threads = DEFAULT_NUM_THREADS);
/**
* @brief Fit a 1D Gaussian with error estimates
* @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
*/
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,
int n_threads = DEFAULT_NUM_THREADS
);
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
NDView<double, 2> chi2_out, int n_threads = DEFAULT_NUM_THREADS);
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);
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++
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,
int n_threads = DEFAULT_NUM_THREADS);
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
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, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_threads);
void fit_scurve(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);
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,
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
int n_threads);
void fit_scurve(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);
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);
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);
void fit_scurve2(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);
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,
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
int n_threads);
void fit_scurve2(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);
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);
} // namespace aare

1
include/aare/Frame.hpp
View File

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

55
include/aare/JungfrauDataFile.hpp
View File

@ -3,13 +3,12 @@
#include <filesystem>
#include <vector>
#include "aare/FilePtr.hpp"
#include "aare/defs.hpp"
#include "aare/NDArray.hpp"
#include "aare/FileInterface.hpp"
#include "aare/FilePtr.hpp"
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
namespace aare {
struct JungfrauDataHeader {
uint64_t framenum;
uint64_t bunchid;
@ -18,33 +17,39 @@ struct JungfrauDataHeader{
class JungfrauDataFile : public FileInterface {
size_t m_rows{}; //!< number of rows in the image, from find_frame_size();
size_t m_cols{}; //!< number of columns in the image, from find_frame_size();
size_t
m_cols{}; //!< number of columns in the image, from find_frame_size();
size_t m_bytes_per_frame{}; //!< number of bytes per frame excluding header
size_t m_total_frames{}; //!< total number of frames in the series of files
size_t m_offset{}; //!< file index of the first file, allow starting at non zero file
size_t m_offset{}; //!< file index of the first file, allow starting at non
//!< zero file
size_t m_current_file_index{}; //!< The index of the open file
size_t m_current_frame_index{}; //!< The index of the current frame (with reference to all files)
size_t m_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::string m_base_name; //!< base name used for formatting file names
FilePtr m_fp; //!< RAII wrapper for a FILE*
using pixel_type = uint16_t;
static constexpr size_t header_size = sizeof(JungfrauDataHeader);
static constexpr size_t n_digits_in_file_index = 6; //!< to format file names
static constexpr size_t n_digits_in_file_index =
6; //!< to format file names
public:
JungfrauDataFile(const std::filesystem::path &fname);
std::string base_name() const; //!< get the base name of the file (without path and extension)
std::string base_name()
const; //!< get the base name of the file (without path and extension)
size_t bytes_per_frame() override;
size_t pixels_per_frame() override;
size_t bytes_per_pixel() const;
size_t bitdepth() const override;
void seek(size_t frame_index) override; //!< seek to the given frame index (note not byte offset)
void seek(size_t frame_index)
override; //!< seek to the given frame index (note not byte offset)
size_t tell() override; //!< get the frame index of the file pointer
size_t total_frames() const override;
size_t rows() const override;
@ -63,44 +68,48 @@ class JungfrauDataFile : public FileInterface {
/**
* @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)
*/
void read_into(std::byte *image_buf, JungfrauDataHeader *header = nullptr);
/**
* @brief Read a multiple frames from the file into the given buffer.
* @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer)
* @param image_buf buffer to read the frame into. (Note the caller is
* responsible for allocating the buffer)
* @param n_frames number of frames to read
* @param header pointer to a JungfrauDataHeader or nullptr to skip header)
*/
void read_into(std::byte *image_buf, size_t n_frames, JungfrauDataHeader *header = nullptr);
void read_into(std::byte *image_buf, size_t n_frames,
JungfrauDataHeader *header = nullptr);
/**
* @brief Read a single frame from the file into the given NDArray
* @param image NDArray to read the frame into.
*/
void read_into(NDArray<uint16_t>* image, JungfrauDataHeader* header = nullptr);
void read_into(NDArray<uint16_t> *image,
JungfrauDataHeader *header = nullptr);
JungfrauDataHeader read_header();
std::filesystem::path current_file() const { return fpath(m_current_file_index+m_offset); }
std::filesystem::path current_file() const {
return fpath(m_current_file_index + m_offset);
}
private:
/**
* @brief Find the size of the frame in the file. (256x256, 256x1024, 512x1024)
* @brief Find the size of the frame in the file. (256x256, 256x1024,
* 512x1024)
* @param fname path to the file
* @throws std::runtime_error if the file is empty or the size cannot be determined
* @throws std::runtime_error if the file is empty or the size cannot be
* determined
*/
void find_frame_size(const std::filesystem::path &fname);
void parse_fname(const std::filesystem::path &fname);
void scan_files();
void open_file(size_t file_index);
std::filesystem::path fpath(size_t frame_index) const;
};
} // namespace aare

13
include/aare/NDArray.hpp
View File

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

65
include/aare/NDView.hpp
View File

@ -1,6 +1,6 @@
#pragma once
#include "aare/defs.hpp"
#include "aare/ArrayExpr.hpp"
#include "aare/defs.hpp"
#include <algorithm>
#include <array>
@ -17,7 +17,8 @@ namespace aare {
template <ssize_t Ndim> using Shape = std::array<ssize_t, Ndim>;
// 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)
throw std::runtime_error("Shape size mismatch");
Shape<Ndim> arr;
@ -25,14 +26,18 @@ template <ssize_t Ndim> Shape<Ndim> make_shape(const std::vector<size_t> &shape)
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>
ssize_t element_offset(const Strides &strides, ssize_t i, Ix... 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::fill(strides.begin(), strides.end(), 1);
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;
}
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);
std::array<ssize_t, Ndim> arr{};
std::copy_n(vec.begin(), Ndim, arr.begin());
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:
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)
: 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)
// : buffer_(buffer), 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<>())) {}
// : buffer_(buffer),
// 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...)];
}
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...)];
}
@ -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::minus<T>()); }
NDView &operator*=(const T val) { return elemenwise(val, std::multiplies<T>()); }
NDView &operator/=(const T val) { return elemenwise(val, std::divides<T>()); }
NDView &operator*=(const T val) {
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()))
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());
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_{};
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) {
buffer_[i] = op(buffer_[i], val);
}
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) {
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>
std::ostream &operator<<(std::ostream &os, const NDView<T, Ndim> &arr) {
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;
}
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())});
}

42
include/aare/NumpyFile.hpp
View File

@ -1,9 +1,8 @@
#pragma once
#include "aare/Dtype.hpp"
#include "aare/defs.hpp"
#include "aare/FileInterface.hpp"
#include "aare/NumpyHelpers.hpp"
#include "aare/defs.hpp"
#include <filesystem>
#include <iostream>
@ -11,13 +10,12 @@
namespace aare {
/**
* @brief NumpyFile class to read and write numpy files
* @note derived from FileInterface
* @note implements all the pure virtual functions from FileInterface
* @note documentation for the functions can also be found in the FileInterface class
* @note documentation for the functions can also be found in the FileInterface
* class
*/
class NumpyFile : public FileInterface {
@ -28,26 +26,35 @@ class NumpyFile : public FileInterface {
* @param mode file mode (r, w)
* @param cfg file configuration
*/
explicit NumpyFile(const std::filesystem::path &fname, const std::string &mode = "r", FileConfig cfg = {});
explicit NumpyFile(const std::filesystem::path &fname,
const std::string &mode = "r", FileConfig cfg = {});
void write(Frame &frame);
Frame read_frame() override { return get_frame(this->current_frame++); }
Frame read_frame(size_t frame_number) override { return get_frame(frame_number); }
Frame read_frame(size_t frame_number) override {
return get_frame(frame_number);
}
std::vector<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;
size_t frame_number(size_t frame_index) override { return frame_index; };
size_t bytes_per_frame() override;
size_t pixels_per_frame() override;
void seek(size_t frame_number) override { this->current_frame = frame_number; }
void seek(size_t frame_number) override {
this->current_frame = frame_number;
}
size_t tell() override { return this->current_frame; }
size_t total_frames() const override { return m_header.shape[0]; }
size_t rows() const override { return m_header.shape[1]; }
size_t cols() const override { return m_header.shape[2]; }
size_t bitdepth() const override { return m_header.dtype.bitdepth(); }
DetectorType detector_type() const override { return DetectorType::Unknown; }
DetectorType detector_type() const override {
return DetectorType::Unknown;
}
/**
* @brief get the data type of the numpy file
@ -70,7 +77,8 @@ class NumpyFile : public FileInterface {
template <typename T, size_t NDim> NDArray<T, NDim> load() {
NDArray<T, NDim> arr(make_shape<NDim>(m_header.shape));
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);
if (rc != static_cast<size_t>(arr.size())) {
@ -78,16 +86,20 @@ class NumpyFile : public FileInterface {
}
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());
}
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());
}
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());
}
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());
}

9
include/aare/NumpyHelpers.hpp
View File

@ -40,15 +40,18 @@ bool parse_bool(const std::string &in);
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);
}
bool is_digits(const std::string &str);
aare::Dtype parse_descr(std::string typestring);
size_t write_header(const std::filesystem::path &fname, const NumpyHeader &header);
size_t write_header(const std::filesystem::path &fname,
const NumpyHeader &header);
size_t write_header(std::ostream &out, const NumpyHeader &header);
} // namespace NumpyHelpers

26
include/aare/Pedestal.hpp
View File

@ -42,9 +42,7 @@ template <typename SUM_TYPE = double> class Pedestal {
}
~Pedestal() = default;
NDArray<SUM_TYPE, 2> mean() {
return m_mean;
}
NDArray<SUM_TYPE, 2> mean() { return m_mean; }
SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
return m_mean(row, col);
@ -71,8 +69,6 @@ template <typename SUM_TYPE = double> class Pedestal {
return variance_array;
}
NDArray<SUM_TYPE, 2> std() {
NDArray<SUM_TYPE, 2> standard_deviation_array({m_rows, m_cols});
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;
}
void clear() {
m_sum = 0;
m_sum2 = 0;
@ -92,8 +86,6 @@ template <typename SUM_TYPE = double> class Pedestal {
m_mean = 0;
}
void clear(const uint32_t row, const uint32_t col) {
m_sum(row, col) = 0;
m_sum2(row, col) = 0;
@ -101,8 +93,6 @@ template <typename SUM_TYPE = double> class Pedestal {
m_mean(row, col) = 0;
}
template <typename T> void push(NDView<T, 2> frame) {
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) {
assert(frame.rows() == static_cast<size_t>(m_rows) &&
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_sum2(row, col) += val * val - m_sum2(row, col) / m_samples;
}
//Since we just did a push we know that m_cur_samples(row, col) is at least 1
// Since we just did a push we know that m_cur_samples(row, col) is at
// least 1
m_mean(row, col) = m_sum(row, col) / m_cur_samples(row, col);
}
@ -183,7 +171,8 @@ template <typename SUM_TYPE = double> class Pedestal {
m_cur_samples(row, col)++;
} else {
m_sum(row, col) += val - m_sum(row, col) / m_cur_samples(row, col);
m_sum2(row, col) += val * val - m_sum2(row, col) / m_cur_samples(row, col);
m_sum2(row, col) +=
val * val - m_sum2(row, col) / m_cur_samples(row, col);
}
}
@ -191,9 +180,7 @@ template <typename SUM_TYPE = double> class Pedestal {
* @brief Update the mean of the pedestal. This is used after having done
* push_no_update. It is not necessary to call this function after push.
*/
void update_mean(){
m_mean = m_sum / m_cur_samples;
}
void update_mean() { m_mean = m_sum / m_cur_samples; }
template <typename T>
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_mean(row, col) = m_sum(row, col) / m_samples;
}
};
} // namespace aare

2
include/aare/PixelMap.hpp
View File

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

18
include/aare/ProducerConsumerQueue.hpp
View File

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

17
include/aare/RawFile.hpp
View File

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

4
include/aare/RawMasterFile.hpp
View File

@ -61,7 +61,6 @@ class ScanParameters {
void increment_stop();
};
/**
* @brief Class for parsing a master file either in our .json format or the old
* .raw format
@ -101,7 +100,6 @@ class RawMasterFile {
std::optional<ROI> m_roi;
public:
RawMasterFile(const std::filesystem::path &fpath);
@ -129,10 +127,8 @@ class RawMasterFile {
std::optional<size_t> number_of_rows() const;
std::optional<uint8_t> quad() const;
std::optional<ROI> roi() const;
ScanParameters scan_parameters() const;
private:

29
include/aare/RawSubFile.hpp
View File

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

38
include/aare/VarClusterFinder.hpp
View File

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

32
include/aare/algorithm.hpp
View File

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

17
include/aare/decode.hpp
View File

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

29
include/aare/defs.hpp
View File

@ -3,16 +3,15 @@
#include "aare/Dtype.hpp"
#include <array>
#include <stdexcept>
#include <cassert>
#include <cstdint>
#include <cstring>
#include <stdexcept>
#include <string>
#include <string_view>
#include <variant>
#include <vector>
/**
* @brief LOCATION macro to get the current location in the code
*/
@ -20,28 +19,24 @@
std::string(__FILE__) + std::string(":") + std::to_string(__LINE__) + \
":" + std::string(__func__) + ":"
#ifdef AARE_CUSTOM_ASSERT
#define AARE_ASSERT(expr) \
if (expr)\
{}\
else\
if (expr) { \
} else \
aare::assert_failed(LOCATION + " Assertion failed: " + #expr + "\n");
#else
#define AARE_ASSERT(cond) \
do { (void)sizeof(cond); } while(0)
do { \
(void)sizeof(cond); \
} while (0)
#endif
namespace aare {
inline constexpr size_t bits_per_byte = 8;
void assert_failed(const std::string &msg);
class DynamicCluster {
public:
int cluster_sizeX;
@ -179,9 +174,9 @@ template <typename T> struct t_xy {
};
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 {
int origin_x{};
@ -193,8 +188,8 @@ struct ModuleGeometry{
};
/**
* @brief Class to hold the geometry of a detector. Number of modules, their size and where pixel 0
* for each module is located
* @brief Class to hold the geometry of a detector. Number of modules, their
* size and where pixel 0 for each module is located
*/
struct DetectorGeometry {
int modules_x{};
@ -221,7 +216,6 @@ struct ROI{
}
};
using dynamic_shape = std::vector<ssize_t>;
// TODO! Can we uniform enums between the libraries?
@ -243,7 +237,8 @@ enum class DetectorType {
Gotthard2,
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_old,
Unknown

3
include/aare/geo_helpers.hpp
View File

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

12
include/aare/logger.hpp
View File

@ -1,7 +1,6 @@
#pragma once
/*Utility to log to console*/
#include <iostream>
#include <sstream>
#include <sys/time.h>
@ -27,7 +26,6 @@ namespace aare {
#define RESET "\x1b[0m"
#define BOLD "\x1b[1m"
enum TLogLevel {
logERROR,
logWARNING,
@ -37,7 +35,8 @@ enum TLogLevel {
logINFOCYAN,
logINFOMAGENTA,
logINFO,
logDEBUG, // constructors, destructors etc. should still give too much output
logDEBUG, // constructors, destructors etc. should still give too much
// output
logDEBUG1,
logDEBUG2,
logDEBUG3,
@ -47,7 +46,9 @@ enum TLogLevel {
// Compiler should optimize away anything below this value
#ifndef AARE_LOG_LEVEL
#define AARE_LOG_LEVEL "LOG LEVEL NOT SET IN CMAKE" //This is configured in the main CMakeLists.txt
#define AARE_LOG_LEVEL \
"LOG LEVEL NOT SET IN CMAKE" // This is configured in the main
// CMakeLists.txt
#endif
#define __AT__ \
@ -72,7 +73,8 @@ class Logger {
std::clog << os.str() << std::flush; // Single write
}
static TLogLevel &ReportingLevel() { // singelton eeh TODO! Do we need a runtime option?
static TLogLevel &
ReportingLevel() { // singelton eeh TODO! Do we need a runtime option?
static TLogLevel reportingLevel = logDEBUG5;
return reportingLevel;
}

2
include/aare/utils/par.hpp
View File

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

4
python/src/bind_Cluster.hpp
View File

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

4
python/src/bind_ClusterCollector.hpp
View File

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

3
python/src/bind_ClusterFile.hpp
View File

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

9
python/src/bind_ClusterFileSink.hpp
View File

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

8
python/src/bind_ClusterFinderMT.hpp
View File

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

3
python/src/bind_ClusterVector.hpp
View File

@ -44,7 +44,8 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
auto *vec = new std::vector<Type>(self.sum());
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());
return return_vector(vec);
})

56
python/src/ctb_raw_file.hpp
View File

@ -6,8 +6,8 @@
#include "aare/RawMasterFile.hpp"
#include "aare/RawSubFile.hpp"
#include "aare/defs.hpp"
#include "aare/decode.hpp"
#include "aare/defs.hpp"
// #include "aare/fClusterFileV2.hpp"
#include "np_helper.hpp"
@ -27,61 +27,69 @@ using namespace ::aare;
void define_ctb_raw_file_io_bindings(py::module &m) {
m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) {
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
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);
// 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<uint16_t, 2> output_view(output.mutable_data(), {output.shape(0), output.shape(1)});
NDView<uint64_t, 2> input_view(
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);
return output;
});
m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
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
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);
// 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<uint16_t, 2> output_view(output.mutable_data(), {output.shape(0), output.shape(1)});
NDView<uint64_t, 2> input_view(
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);
return output;
});
m.def(
"apply_custom_weights",
[](py::array_t<uint16_t, py::array::c_style | py::array::forcecast> &input,
m.def("apply_custom_weights",
[](py::array_t<uint16_t, py::array::c_style | py::array::forcecast>
&input,
py::array_t<double, py::array::c_style | py::array::forcecast>
&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::array_t<double> output(buf.shape);
// Use NDViews to call into the C++ library
auto weights_view = make_view_1d(weights);
NDView<uint16_t, 1> input_view(input.mutable_data(), {input.size()});
NDView<double, 1> output_view(output.mutable_data(), {output.size()});
NDView<uint16_t, 1> input_view(input.mutable_data(),
{input.size()});
NDView<double, 1> output_view(output.mutable_data(),
{output.size()});
apply_custom_weights(input_view, output_view, weights_view);
return output;
@ -103,7 +111,8 @@ py::class_<CtbRawFile>(m, "CtbRawFile")
// always read bytes
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());
return py::make_tuple(header, image);
@ -116,5 +125,4 @@ py::class_<CtbRawFile>(m, "CtbRawFile")
&CtbRawFile::image_size_in_bytes)
.def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
}

44
python/src/file.hpp
View File

@ -20,9 +20,6 @@
namespace py = pybind11;
using namespace ::aare;
// Disable warnings for unused parameters, as we ignore some
// in the __exit__ method
#pragma GCC diagnostic push
@ -30,7 +27,6 @@ using namespace ::aare;
void define_file_io_bindings(py::module &m) {
py::enum_<DetectorType>(m, "DetectorType")
.value("Jungfrau", DetectorType::Jungfrau)
.value("Eiger", DetectorType::Eiger)
@ -41,13 +37,10 @@ void define_file_io_bindings(py::module &m) {
.value("ChipTestBoard", DetectorType::ChipTestBoard)
.value("Unknown", DetectorType::Unknown);
PYBIND11_NUMPY_DTYPE(DetectorHeader, frameNumber, expLength, packetNumber,
bunchId, timestamp, modId, row, column, reserved,
debug, roundRNumber, detType, version, packetMask);
py::class_<File>(m, "File")
.def(py::init([](const std::filesystem::path &fname) {
return File(fname, "r", {});
@ -112,9 +105,11 @@ void define_file_io_bindings(py::module &m) {
reinterpret_cast<std::byte *>(image.mutable_data()));
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
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) {
throw std::runtime_error("No frames left in file");
}
@ -129,21 +124,20 @@ void define_file_io_bindings(py::module &m) {
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()),
n_frames);
return image;
})
.def("__enter__", [](File &self) { return &self; })
.def("__exit__",
[](File &self,
const std::optional<pybind11::type> &exc_type,
[](File &self, const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::object> &exc_value,
const std::optional<pybind11::object> &traceback) {
// self.close();
})
.def("__iter__", [](File &self) { return &self; })
.def("__next__", [](File &self) {
try {
const uint8_t item_size = self.bytes_per_pixel();
py::array image;
@ -166,7 +160,6 @@ void define_file_io_bindings(py::module &m) {
}
});
py::class_<FileConfig>(m, "FileConfig")
.def(py::init<>())
.def_readwrite("rows", &FileConfig::rows)
@ -183,8 +176,6 @@ void define_file_io_bindings(py::module &m) {
return "<FileConfig: " + a.to_string() + ">";
});
py::class_<ScanParameters>(m, "ScanParameters")
.def(py::init<const std::string &>())
.def(py::init<const ScanParameters &>())
@ -195,7 +186,6 @@ void define_file_io_bindings(py::module &m) {
.def_property_readonly("stop", &ScanParameters::stop)
.def_property_readonly("step", &ScanParameters::step);
py::class_<ROI>(m, "ROI")
.def(py::init<>())
.def(py::init<ssize_t, ssize_t, ssize_t, ssize_t>(), py::arg("xmin"),
@ -204,23 +194,21 @@ void define_file_io_bindings(py::module &m) {
.def_readwrite("xmax", &ROI::xmax)
.def_readwrite("ymin", &ROI::ymin)
.def_readwrite("ymax", &ROI::ymax)
.def("__str__", [](const ROI& self){
return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}", self.xmin, self.xmax, self.ymin, self.ymax);
.def("__str__",
[](const ROI &self) {
return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}",
self.xmin, self.xmax, self.ymin, self.ymax);
})
.def("__repr__", [](const ROI& self){
return fmt::format("<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin, self.xmax, self.ymin, self.ymax);
.def("__repr__",
[](const ROI &self) {
return fmt::format(
"<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin,
self.xmax, self.ymin, self.ymax);
})
.def("__iter__", [](const ROI &self) {
return py::make_iterator(&self.xmin, &self.ymax + 1); // NOLINT
});
#pragma GCC diagnostic pop
// py::class_<ClusterHeader>(m, "ClusterHeader")
// .def(py::init<>())

13
python/src/fit.hpp
View File

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

9
python/src/jungfrau_data_file.hpp
View File

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

19
python/src/module.cpp
View File

@ -3,10 +3,10 @@
// New style file naming
#include "bind_Cluster.hpp"
#include "bind_ClusterCollector.hpp"
#include "bind_ClusterFinder.hpp"
#include "bind_ClusterFinderMT.hpp"
#include "bind_ClusterFile.hpp"
#include "bind_ClusterFileSink.hpp"
#include "bind_ClusterFinder.hpp"
#include "bind_ClusterFinderMT.hpp"
#include "bind_ClusterVector.hpp"
// TODO! migrate the other names
@ -14,13 +14,13 @@
#include "file.hpp"
#include "fit.hpp"
#include "interpolation.hpp"
#include "raw_sub_file.hpp"
#include "raw_master_file.hpp"
#include "raw_file.hpp"
#include "pixel_map.hpp"
#include "var_cluster.hpp"
#include "pedestal.hpp"
#include "jungfrau_data_file.hpp"
#include "pedestal.hpp"
#include "pixel_map.hpp"
#include "raw_file.hpp"
#include "raw_master_file.hpp"
#include "raw_sub_file.hpp"
#include "var_cluster.hpp"
// Pybind stuff
#include <pybind11/pybind11.h>
@ -34,7 +34,8 @@ T - Storage type of the cluster data (int, float, double)
N - Number of rows in the cluster
M - Number of columns in the cluster
U - Type of the pixel data (e.g., uint16_t)
TYPE_CODE - A character representing the type code (e.g., 'i' for int, 'd' for double, 'f' for float)
TYPE_CODE - A character representing the type code (e.g., 'i' for int, 'd' for
double, 'f' for float)
*/
#define DEFINE_CLUSTER_BINDINGS(T, N, M, U, TYPE_CODE) \

20
python/src/pedestal.hpp
View File

@ -9,7 +9,8 @@
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())
.def(py::init<int, 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();
return return_image_data(mea);
})
.def("variance", [](Pedestal<SUM_TYPE> &self) {
.def("variance",
[](Pedestal<SUM_TYPE> &self) {
auto var = new NDArray<SUM_TYPE, 2>{};
*var = self.variance();
return return_image_data(var);
})
.def("std", [](Pedestal<SUM_TYPE> &self) {
.def("std",
[](Pedestal<SUM_TYPE> &self) {
auto std = new NDArray<SUM_TYPE, 2>{};
*std = self.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);
})
// 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);
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);
pedestal.push_no_update(v);
}, py::arg().noconvert())
},
py::arg().noconvert())
.def("update_mean", &Pedestal<SUM_TYPE>::update_mean);
}

31
python/src/pixel_map.hpp
View File

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

3
python/src/raw_file.hpp
View File

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

3
python/src/raw_master_file.hpp
View File

@ -57,7 +57,8 @@ void define_raw_master_file_bindings(py::module &m) {
.def_property_readonly("total_frames_expected",
&RawMasterFile::total_frames_expected)
.def_property_readonly("geometry", &RawMasterFile::geometry)
.def_property_readonly("analog_samples", &RawMasterFile::analog_samples, R"(
.def_property_readonly("analog_samples", &RawMasterFile::analog_samples,
R"(
Number of analog samples
Returns

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

7
python/src/var_cluster.hpp
View File

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

71
src/ClusterFile.cpp
View File

@ -31,9 +31,7 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
}
}
void ClusterFile::set_roi(ROI roi){
m_roi = roi;
}
void ClusterFile::set_roi(ROI roi) { m_roi = roi; }
void ClusterFile::set_noise_map(const NDView<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
uint32_t n_clusters = clusters.size();
if (fwrite(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
throw std::runtime_error(LOCATION + "Could not write number of clusters");
throw std::runtime_error(LOCATION +
"Could not write number of clusters");
}
// Now write the clusters in the frame
if(fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp)!=clusters.size()){
if (fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp) !=
clusters.size()) {
throw std::runtime_error(LOCATION + "Could not write clusters");
}
}
ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
if (m_mode != "r") {
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") {
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;
}
ClusterVector<int32_t> ClusterFile::read_clusters_with_cut(size_t n_clusters) {
ClusterVector<int32_t> clusters(3, 3);
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) {
Cluster3x3 c = read_one_cluster();
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) {
// sanity check
if (m_num_left) {
throw std::runtime_error(LOCATION + "Entered second loop with clusters left\n");
throw std::runtime_error(
LOCATION + "Entered second loop with clusters left\n");
}
int32_t frame_number = 0; // frame number needs to be 4 bytes!
while (fread(&frame_number, sizeof(frame_number), 1, fp)) {
if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) {
clusters.set_frame_number(frame_number); //cluster vector will hold the last frame number
clusters.set_frame_number(
frame_number); // cluster vector will hold the last frame
// number
while (m_num_left && clusters.size() < n_clusters) {
Cluster3x3 c = read_one_cluster();
if (is_selected(c)) {
clusters.push_back(c.x, c.y, reinterpret_cast<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)
break;
}
}
if (m_gain_map)
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
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);
@ -264,7 +268,6 @@ ClusterVector<int32_t> ClusterFile::read_frame_with_cut() {
throw std::runtime_error("Could not read frame number");
}
if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) {
throw std::runtime_error("Could not read number of clusters");
}
@ -283,8 +286,6 @@ ClusterVector<int32_t> ClusterFile::read_frame_with_cut() {
return clusters;
}
bool ClusterFile::is_selected(Cluster3x3 &cl) {
// Should fail fast
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_3x3 = cl.sum(); // sum of all pixels
auto noise = (*m_noise_map)(cl.y, cl.x); //TODO! check if this is correct
auto noise =
(*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct
if (sum_1x1 <= noise || sum_2x2 <= 2 * noise || sum_3x3 <= 3 * noise) {
return false;
}
@ -316,7 +318,8 @@ NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
}
}else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){
} else if (clusters.cluster_size_x() == 2 ||
clusters.cluster_size_y() == 2) {
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2(clusters.at<Cluster2x2>(i));
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
* containing etay, etax and the corner of the cluster.
* @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2
* struct containing etay, etax and the corner of the cluster.
*/
Eta2 calculate_eta2(Cluster3x3 &cl) {
Eta2 eta{};
@ -347,38 +350,30 @@ Eta2 calculate_eta2(Cluster3x3 &cl) {
switch (c) {
case cBottomLeft:
if ((cl.data[3] + cl.data[4]) != 0)
eta.x =
static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[1] + cl.data[4]) != 0)
eta.y =
static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.c = cBottomLeft;
break;
case cBottomRight:
if ((cl.data[2] + cl.data[5]) != 0)
eta.x =
static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
eta.x = static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
if ((cl.data[1] + cl.data[4]) != 0)
eta.y =
static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.c = cBottomRight;
break;
case cTopLeft:
if ((cl.data[7] + cl.data[4]) != 0)
eta.x =
static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta.y =
static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.c = cTopLeft;
break;
case cTopRight:
if ((cl.data[5] + cl.data[4]) != 0)
eta.x =
static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
eta.x = static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta.y =
static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.c = cTopRight;
break;
// no default to allow compiler to warn about missing cases
@ -386,7 +381,6 @@ Eta2 calculate_eta2(Cluster3x3 &cl) {
return eta;
}
Eta2 calculate_eta2(Cluster2x2 &cl) {
Eta2 eta{};
if ((cl.data[0] + cl.data[1]) != 0)
@ -398,5 +392,4 @@ Eta2 calculate_eta2(Cluster2x2 &cl) {
return eta;
}
} // namespace aare

4
src/ClusterFile.test.cpp
View File

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

18
src/CtbRawFile.cpp
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");
}
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);
}
@ -29,7 +30,8 @@ void CtbRawFile::read_into(std::byte *image_buf, DetectorHeader* header) {
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++;
}
@ -38,13 +40,16 @@ void CtbRawFile::seek(size_t frame_number) {
open_data_file(index);
}
size_t frame_number_in_file = frame_number % m_master.max_frames_per_file();
m_file.seekg((sizeof(DetectorHeader)+m_master.image_size_in_bytes()) * frame_number_in_file);
m_file.seekg((sizeof(DetectorHeader) + m_master.image_size_in_bytes()) *
frame_number_in_file);
m_current_frame = frame_number;
}
size_t CtbRawFile::tell() const { return m_current_frame; }
size_t CtbRawFile::image_size_in_bytes() const { return m_master.image_size_in_bytes(); }
size_t CtbRawFile::image_size_in_bytes() const {
return m_master.image_size_in_bytes();
}
size_t CtbRawFile::frames_in_file() const { return m_master.frames_in_file(); }
@ -63,12 +68,11 @@ void CtbRawFile::open_data_file(size_t subfile_index) {
throw std::runtime_error(LOCATION + "Subfile index out of range");
}
m_current_subfile = subfile_index;
m_file = std::ifstream(m_master.data_fname(0, subfile_index), std::ios::binary); // only one module for CTB
m_file = std::ifstream(m_master.data_fname(0, subfile_index),
std::ios::binary); // only one module for CTB
if (!m_file.is_open()) {
throw std::runtime_error(LOCATION + "Could not open data file");
}
}
} // namespace aare

34
src/Dtype.cpp
View File

@ -10,7 +10,8 @@ namespace aare {
* @brief Construct a DType object from a type_info object
* @param t type_info object
* @throw runtime_error if the type is not supported
* @note supported types are: int8_t, uint8_t, int16_t, uint16_t, int32_t, uint32_t, int64_t, uint64_t, float, double
* @note supported types are: int8_t, uint8_t, int16_t, uint16_t, int32_t,
* uint32_t, int64_t, uint64_t, float, double
* @note the type_info object is obtained using typeid (e.g. typeid(int))
*/
Dtype::Dtype(const std::type_info &t) {
@ -35,7 +36,8 @@ Dtype::Dtype(const std::type_info &t) {
else if (t == typeid(double))
m_type = TypeIndex::DOUBLE;
else
throw std::runtime_error("Could not construct data type. Type not supported.");
throw std::runtime_error(
"Could not construct data type. Type not supported.");
}
/**
@ -63,7 +65,8 @@ uint8_t Dtype::bitdepth() const {
case TypeIndex::NONE:
return 0;
default:
throw std::runtime_error(LOCATION + "Could not get bitdepth. Type not supported.");
throw std::runtime_error(LOCATION +
"Could not get bitdepth. Type not supported.");
}
}
@ -138,7 +141,8 @@ Dtype Dtype::from_bitdepth(uint8_t bitdepth) {
case 64:
return Dtype(TypeIndex::UINT64);
default:
throw std::runtime_error("Could not construct data type from bitdepth.");
throw std::runtime_error(
"Could not construct data type from bitdepth.");
}
}
/**
@ -175,17 +179,27 @@ std::string Dtype::to_string() const {
case TypeIndex::DOUBLE:
return "f8";
case TypeIndex::ERROR:
throw std::runtime_error("Could not get string representation. Type not supported.");
throw std::runtime_error(
"Could not get string representation. Type not supported.");
case TypeIndex::NONE:
throw std::runtime_error("Could not get string representation. Type not supported.");
throw std::runtime_error(
"Could not get string representation. Type not supported.");
}
return {};
}
bool Dtype::operator==(const Dtype &other) const noexcept { return m_type == other.m_type; }
bool Dtype::operator!=(const Dtype &other) const noexcept { return !(*this == other); }
bool Dtype::operator==(const Dtype &other) const noexcept {
return m_type == other.m_type;
}
bool Dtype::operator!=(const Dtype &other) const noexcept {
return !(*this == other);
}
bool Dtype::operator==(const std::type_info &t) const { return Dtype(t) == *this; }
bool Dtype::operator!=(const std::type_info &t) const { return Dtype(t) != *this; }
bool Dtype::operator==(const std::type_info &t) const {
return Dtype(t) == *this;
}
bool Dtype::operator!=(const std::type_info &t) const {
return Dtype(t) != *this;
}
} // namespace aare

4
src/Dtype.test.cpp
View File

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

13
src/File.cpp
View File

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

9
src/FilePtr.cpp
View File

@ -6,10 +6,12 @@
namespace aare {
FilePtr::FilePtr(const std::filesystem::path& fname, const std::string& mode = "rb") {
FilePtr::FilePtr(const std::filesystem::path &fname,
const std::string &mode = "rb") {
fp_ = fopen(fname.c_str(), mode.c_str());
if (!fp_)
throw std::runtime_error(fmt::format("Could not open: {}", fname.c_str()));
throw std::runtime_error(
fmt::format("Could not open: {}", fname.c_str()));
}
FilePtr::FilePtr(FilePtr &&other) { std::swap(fp_, other.fp_); }
@ -24,7 +26,8 @@ FILE *FilePtr::get() { return fp_; }
ssize_t FilePtr::tell() {
auto pos = ftell(fp_);
if (pos == -1)
throw std::runtime_error(fmt::format("Error getting file position: {}", error_msg()));
throw std::runtime_error(
fmt::format("Error getting file position: {}", error_msg()));
return pos;
}
FilePtr::~FilePtr() {

134
src/Fit.cpp
View File

@ -1,13 +1,12 @@
#include "aare/Fit.hpp"
#include "aare/utils/task.hpp"
#include "aare/utils/par.hpp"
#include "aare/utils/task.hpp"
#include <lmcurve2.h>
#include <lmfit.hpp>
#include <thread>
#include <array>
namespace aare {
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) {
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) {
@ -47,7 +48,9 @@ NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> 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) {
@ -91,7 +94,8 @@ NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
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};
auto e = std::max_element(y.begin(), y.end());
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
// assuming equally spaced x values
auto delta = x[1] - x[0];
start_par[2] =
std::count_if(y.begin(), y.end(),
start_par[2] = std::count_if(y.begin(), y.end(),
[e](double val) { return val > *e / 2; }) *
delta / 2.35;
return start_par;
}
std::array<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) {
// Estimate the initial parameters for the fit
std::array<double, 2> start_par{0, 0};
auto y2 = std::max_element(y.begin(), y.end());
auto x2 = x[std::distance(y.begin(), y2)];
auto y1 = std::min_element(y.begin(), y.end());
@ -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");
}
// /* Collection of output parameters for status info. */
// typedef struct {
// double fnorm; /* norm of the residue vector fvec. */
@ -153,23 +154,32 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
// */
// } lm_status_struct;
lm_status_struct status;
par_out = gaus_init_par(x, y);
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);
// n_par - Number of free variables. Length of parameter vector par.
// par - Parameter vector. On input, it must contain a reasonable guess. On output, it contains the solution found to minimize ||r||.
// parerr - Parameter uncertainties vector. Array of length n_par or NULL. On output, unless it or covar is NULL, it contains the weighted parameter uncertainties for the found parameters.
// covar - Covariance matrix. Array of length n_par * n_par or NULL. On output, unless it is NULL, it contains the covariance matrix.
// m_dat - Number of data points. Length of vectors t, y, dy. Must statisfy n_par <= m_dat.
// t - Array of length m_dat. Contains the abcissae (time, or "x") for which function f will be evaluated.
// y - Array of length m_dat. Contains the ordinate values that shall be fitted.
// dy - Array of length m_dat. Contains the standard deviations of the values y.
// f - A user-supplied parametric function f(ti;par).
// control - Parameter collection for tuning the fit procedure. In most cases, the default &lm_control_double is adequate. If f is only computed with single-precision accuracy, &lm_control_float should be used. Parameters are explained in lmmin2(3).
// status - A record used to return information about the minimization process: For details, see lmmin2(3).
// void lmcurve2( const int n_par, double *par, double *parerr, double
// *covar, const int m_dat, const double *t, const double *y, const double
// *dy, double (*f)( const double ti, const double *par ), const
// lm_control_struct *control, lm_status_struct *status); n_par - Number of
// free variables. Length of parameter vector par. par - Parameter vector.
// On input, it must contain a reasonable guess. On output, it contains the
// solution found to minimize ||r||. parerr - Parameter uncertainties
// vector. Array of length n_par or NULL. On output, unless it or covar is
// NULL, it contains the weighted parameter uncertainties for the found
// parameters. covar - Covariance matrix. Array of length n_par * n_par or
// NULL. On output, unless it is NULL, it contains the covariance matrix.
// m_dat - Number of data points. Length of vectors t, y, dy. Must statisfy
// n_par <= m_dat. t - Array of length m_dat. Contains the abcissae (time,
// or "x") for which function f will be evaluated. y - Array of length
// m_dat. Contains the ordinate values that shall be fitted. dy - Array of
// length m_dat. Contains the standard deviations of the values y. f - A
// user-supplied parametric function f(ti;par). control - Parameter
// collection for tuning the fit procedure. In most cases, the default
// &lm_control_double is adequate. If f is only computed with
// single-precision accuracy, &lm_control_float should be used. Parameters
// are explained in lmmin2(3). status - A record used to return information
// about the minimization process: For details, see lmmin2(3).
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
x.size(), x.data(), y.data(), y_err.data(), aare::func::gaus,
@ -178,12 +188,14 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
// Calculate chi2
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 += std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2);
chi2 +=
std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_gaus(NDView<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) {
@ -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,
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,
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
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
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
chi2 +=
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_pol1(NDView<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) {
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
NDView<double, 2> chi2_out, int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
@ -249,15 +262,14 @@ void fit_pol1(NDView<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),
{par_err_out.shape(2)});
fit_pol1(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
fit_pol1(x, y_view, y_err_view, par_out_view, par_err_out_view,
chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
NDArray<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 ~~
// 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
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());
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) {
if (y[i] >= start_par[4]) {
start_par[2] = x[i];
@ -334,7 +348,8 @@ NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y) {
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);
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
void fit_scurve(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) {
void fit_scurve(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) {
// Check that we have the correct sizes
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
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
chi2 +=
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_scurve(NDView<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,
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) {
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),
{par_err_out.shape(2)});
fit_scurve(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
fit_scurve(x, y_view, y_err_view, par_out_view,
par_err_out_view, chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
// SCURVE2 ---
std::array<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
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());
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) {
if (y[i] <= start_par[4]) {
start_par[2] = x[i];
@ -446,7 +465,8 @@ NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y) {
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);
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
void fit_scurve2(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) {
void fit_scurve2(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) {
// Check that we have the correct sizes
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
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
chi2 +=
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_scurve2(NDView<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,
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) {
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),
{par_err_out.shape(2)});
fit_scurve2(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
fit_scurve2(x, y_view, y_err_view, par_out_view,
par_err_out_view, chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
} // namespace aare

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

4
src/Frame.test.cpp
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
for (size_t i = 0; i < rows; i++) {
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);
if (i == 5 && j == 7) {
REQUIRE(*data == value);
@ -150,4 +151,3 @@ TEST_CASE("test explicit copy constructor") {
REQUIRE(frame2.bytes() == rows * cols * bitdepth / 8);
REQUIRE(frame2.data() != data);
}

14
src/JungfrauDataFile.cpp
View File

@ -19,7 +19,6 @@ JungfrauDataFile::JungfrauDataFile(const std::filesystem::path &fname) {
open_file(m_current_file_index);
}
// FileInterface
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())};
}
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; }
@ -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()) {
throw std::runtime_error(LOCATION +
"Image shape does not match file size: " + std::to_string(rows()) + "x" + std::to_string(cols()));
throw std::runtime_error(
LOCATION + "Image shape does not match file size: " +
std::to_string(rows()) + "x" + std::to_string(cols()));
}
read_into(reinterpret_cast<std::byte *>(image->data()), header);
}
JungfrauDataHeader JungfrauDataFile::read_header() {
JungfrauDataHeader header;
if (auto rc = fread(&header, 1, sizeof(header), m_fp.get());

5
src/JungfrauDataFile.test.cpp
View File

@ -1,7 +1,7 @@
#include "aare/JungfrauDataFile.hpp"
#include <catch2/catch_test_macros.hpp>
#include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
using aare::JungfrauDataFile;
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.current_file().stem() == "AldoJF65k_000003");
}
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(f.tell() == 0);
}

11
src/NDArray.test.cpp
View File

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

9
src/NDView.test.cpp
View File

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

33
src/NumpyFile.cpp
View File

@ -4,16 +4,16 @@
namespace aare {
NumpyFile::NumpyFile(const std::filesystem::path &fname, const std::string &mode, FileConfig cfg) {
NumpyFile::NumpyFile(const std::filesystem::path &fname,
const std::string &mode, FileConfig cfg) {
// TODO! add opts to constructor
m_mode = mode;
if (mode == "r") {
fp = fopen(fname.string().c_str(), "rb");
if (!fp) {
throw std::runtime_error(fmt::format("Could not open: {} for reading", fname.string()));
throw std::runtime_error(
fmt::format("Could not open: {} for reading", fname.string()));
}
load_metadata();
} else if (mode == "w") {
@ -24,11 +24,15 @@ NumpyFile::NumpyFile(const std::filesystem::path &fname, const std::string &mode
m_header.shape = {0, cfg.rows, cfg.cols};
fp = fopen(fname.string().c_str(), "wb");
if (!fp) {
throw std::runtime_error(fmt::format("Could not open: {} for reading", fname.string()));
throw std::runtime_error(
fmt::format("Could not open: {} for reading", fname.string()));
}
initial_header_len = aare::NumpyHelpers::write_header(std::filesystem::path(fname.c_str()), m_header);
initial_header_len = aare::NumpyHelpers::write_header(
std::filesystem::path(fname.c_str()), m_header);
}
m_pixels_per_frame = std::accumulate(m_header.shape.begin() + 1, m_header.shape.end(), 1, std::multiplies<>());
m_pixels_per_frame =
std::accumulate(m_header.shape.begin() + 1, m_header.shape.end(), 1,
std::multiplies<>());
m_bytes_per_frame = m_header.dtype.bitdepth() / 8 * m_pixels_per_frame;
}
@ -63,7 +67,8 @@ void NumpyFile::get_frame_into(size_t frame_number, std::byte *image_buf) {
if (frame_number > m_header.shape[0]) {
throw std::invalid_argument("Frame number out of range");
}
if (fseek(fp, header_size + frame_number * m_bytes_per_frame, SEEK_SET)) // NOLINT
if (fseek(fp, header_size + frame_number * m_bytes_per_frame,
SEEK_SET)) // NOLINT
throw std::runtime_error("Could not seek to frame");
size_t const rc = fread(image_buf, m_bytes_per_frame, 1, fp);
@ -113,7 +118,8 @@ NumpyFile::~NumpyFile() noexcept {
// write header
size_t const rc = fwrite(header_str.c_str(), header_str.size(), 1, fp);
if (rc != 1) {
std::cout << "Error writing header to numpy file in destructor" << std::endl;
std::cout << "Error writing header to numpy file in destructor"
<< std::endl;
}
}
@ -140,8 +146,10 @@ void NumpyFile::load_metadata() {
}
// read version
rc = fread(reinterpret_cast<char *>(&major_ver_), sizeof(major_ver_), 1, fp);
rc += fread(reinterpret_cast<char *>(&minor_ver_), sizeof(minor_ver_), 1, fp);
rc =
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) {
throw std::runtime_error("Error reading numpy version");
}
@ -159,7 +167,8 @@ void NumpyFile::load_metadata() {
if (rc != 1) {
throw std::runtime_error("Error reading header length");
}
header_size = aare::NumpyHelpers::magic_string_length + 2 + header_len_size + header_len;
header_size = aare::NumpyHelpers::magic_string_length + 2 +
header_len_size + header_len;
if (header_size % 16 != 0) {
fmt::print("Warning: header length is not a multiple of 16\n");
}

2
src/NumpyFile.test.cpp
View File

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

42
src/NumpyHelpers.cpp
View File

@ -29,7 +29,8 @@ namespace aare {
std::string NumpyHeader::to_string() const {
std::stringstream sstm;
sstm << "dtype: " << dtype.to_string() << ", fortran_order: " << fortran_order << ' ';
sstm << "dtype: " << dtype.to_string()
<< ", fortran_order: " << fortran_order << ' ';
sstm << "shape: (";
for (auto item : shape)
sstm << item << ',';
@ -37,10 +38,10 @@ std::string NumpyHeader::to_string() const {
return sstm.str();
}
namespace NumpyHelpers {
std::unordered_map<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;
if (keys.empty())
return map;
@ -100,7 +101,8 @@ aare::Dtype parse_descr(std::string typestring) {
constexpr char little_endian_char = '<';
constexpr char big_endian_char = '>';
constexpr char no_endian_char = '|';
constexpr std::array<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'};
const char byteorder_c = typestring[0];
@ -139,7 +141,9 @@ std::string get_value_from_map(const std::string &mapstr) {
return trim(tmp);
}
bool is_digits(const std::string &str) { return std::all_of(str.begin(), str.end(), ::isdigit); }
bool is_digits(const std::string &str) {
return std::all_of(str.begin(), str.end(), ::isdigit);
}
std::vector<std::string> parse_tuple(std::string in) {
std::vector<std::string> v;
@ -215,20 +219,25 @@ inline std::string write_boolean(bool b) {
return "False";
}
inline std::string write_header_dict(const std::string &descr, bool fortran_order, const std::vector<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 shape_s = write_tuple(shape);
return "{'descr': '" + descr + "', 'fortran_order': " + s_fortran_order + ", 'shape': " + shape_s + ", }";
return "{'descr': '" + descr + "', 'fortran_order': " + s_fortran_order +
", 'shape': " + shape_s + ", }";
}
size_t write_header(const std::filesystem::path &fname, const NumpyHeader &header) {
size_t write_header(const std::filesystem::path &fname,
const NumpyHeader &header) {
std::ofstream out(fname, std::ios::binary | std::ios::out);
return write_header(out, header);
}
size_t write_header(std::ostream &out, const NumpyHeader &header) {
std::string const header_dict = write_header_dict(header.dtype.to_string(), header.fortran_order, header.shape);
std::string const header_dict = write_header_dict(
header.dtype.to_string(), header.fortran_order, header.shape);
size_t length = magic_string_length + 2 + 2 + header_dict.length() + 1;
@ -247,17 +256,22 @@ size_t write_header(std::ostream &out, const NumpyHeader &header) {
// write header length
if (version_major == 1 && version_minor == 0) {
auto header_len = static_cast<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)};
out.write(reinterpret_cast<char *>(header_len_le16.data()), 2);
} 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{
static_cast<uint8_t>((header_len >> 0) & 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)};
static_cast<uint8_t>((header_len >> 0) & 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);
}

7
src/NumpyHelpers.test.cpp
View File

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

15
src/Pedestal.test.cpp
View File

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

17
src/PixelMap.cpp
View File

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

55
src/RawFile.cpp
View File

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

57
src/RawFile.test.cpp
View File

@ -1,18 +1,18 @@
#include "aare/RawFile.hpp"
#include "aare/File.hpp"
#include "aare/RawMasterFile.hpp" //needed for ROI
#include "aare/RawFile.hpp"
#include <catch2/catch_test_macros.hpp>
#include <filesystem>
#include "test_config.hpp"
using aare::File;
TEST_CASE("Read number of frames from a jungfrau raw file", "[.integration]") {
auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
auto fpath =
test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r");
@ -20,7 +20,8 @@ TEST_CASE("Read number of frames from a jungfrau raw file", "[.integration]") {
}
TEST_CASE("Read frame numbers from a jungfrau raw file", "[.integration]") {
auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
auto fpath =
test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r");
@ -36,7 +37,8 @@ TEST_CASE("Read frame numbers from a jungfrau raw file", "[.integration]") {
}
TEST_CASE("Read a frame number too high throws", "[.integration]") {
auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
auto fpath =
test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r");
@ -49,8 +51,10 @@ TEST_CASE("Read a frame number too high throws", "[.integration]") {
REQUIRE_THROWS(f.frame_number(10));
}
TEST_CASE("Read a frame numbers where the subfile is missing throws", "[.integration]") {
auto fpath = test_data_path() / "jungfrau" / "jungfrau_missing_subfile_master_0.json";
TEST_CASE("Read a frame numbers where the subfile is missing throws",
"[.integration]") {
auto fpath = test_data_path() / "jungfrau" /
"jungfrau_missing_subfile_master_0.json";
REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r");
@ -69,15 +73,18 @@ TEST_CASE("Read a frame numbers where the subfile is missing throws", "[.integra
REQUIRE_THROWS(f.frame_number(10));
}
TEST_CASE("Read data from a jungfrau 500k single port raw file", "[.integration]") {
auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
TEST_CASE("Read data from a jungfrau 500k single port raw file",
"[.integration]") {
auto fpath =
test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r");
// we know this file has 10 frames with pixel 0,0 being: 2123, 2051, 2109, 2117, 2089, 2095, 2072, 2126, 2097, 2102
std::vector<uint16_t> pixel_0_0 = {2123, 2051, 2109, 2117, 2089, 2095, 2072, 2126, 2097, 2102};
// we know this file has 10 frames with pixel 0,0 being: 2123, 2051, 2109,
// 2117, 2089, 2095, 2072, 2126, 2097, 2102
std::vector<uint16_t> pixel_0_0 = {2123, 2051, 2109, 2117, 2089,
2095, 2072, 2126, 2097, 2102};
for (size_t i = 0; i < 10; i++) {
auto frame = f.read_frame();
CHECK(frame.rows() == 512);
@ -100,10 +107,12 @@ TEST_CASE("Read frame numbers from a raw file", "[.integration]") {
}
TEST_CASE("Compare reading from a numpy file with a raw file", "[.files]") {
auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_master_0.json";
auto fpath_raw =
test_data_path() / "raw/jungfrau" / "jungfrau_single_master_0.json";
REQUIRE(std::filesystem::exists(fpath_raw));
auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
auto fpath_npy =
test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
REQUIRE(std::filesystem::exists(fpath_npy));
File raw(fpath_raw, "r");
@ -121,17 +130,23 @@ TEST_CASE("Compare reading from a numpy file with a raw file", "[.files]") {
}
TEST_CASE("Read multipart files", "[.integration]") {
auto fpath = test_data_path() / "jungfrau" / "jungfrau_double_master_0.json";
auto fpath =
test_data_path() / "jungfrau" / "jungfrau_double_master_0.json";
REQUIRE(std::filesystem::exists(fpath));
File f(fpath, "r");
// we know this file has 10 frames check read_multiport.py for the values
std::vector<uint16_t> pixel_0_0 = {2099, 2121, 2108, 2084, 2084, 2118, 2066, 2108, 2112, 2116};
std::vector<uint16_t> pixel_0_1 = {2842, 2796, 2865, 2798, 2805, 2817, 2852, 2789, 2792, 2833};
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};
std::vector<uint16_t> pixel_0_0 = {2099, 2121, 2108, 2084, 2084,
2118, 2066, 2108, 2112, 2116};
std::vector<uint16_t> pixel_0_1 = {2842, 2796, 2865, 2798, 2805,
2817, 2852, 2789, 2792, 2833};
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++) {
auto frame = f.read_frame();
@ -152,5 +167,3 @@ TEST_CASE("Read file with unordered frames", "[.integration]") {
File f(fpath);
REQUIRE_THROWS((f.read_frame()));
}

27
src/RawMasterFile.cpp
View File

@ -37,18 +37,15 @@ std::filesystem::path RawFileNameComponents::master_fname() const {
}
std::filesystem::path RawFileNameComponents::data_fname(size_t mod_id,
size_t file_id
) const {
size_t file_id) const {
std::string fmt = "{}_d{}_f{}_{}.raw";
// Before version X we used to name the data files f000000000000
if (m_old_scheme) {
fmt = "{}_d{}_f{:012}_{}.raw";
}
return m_base_path / fmt::format(fmt, m_base_name, mod_id,
file_id, m_file_index);
return m_base_path /
fmt::format(fmt, m_base_name, mod_id, file_id, m_file_index);
}
void RawFileNameComponents::set_old_scheme(bool old_scheme) {
@ -85,14 +82,11 @@ ScanParameters::ScanParameters(const std::string& par){
int ScanParameters::start() const { return m_start; }
int ScanParameters::stop() const { return m_stop; }
void ScanParameters::increment_stop(){
m_stop += 1;
}
void ScanParameters::increment_stop() { m_stop += 1; }
int ScanParameters::step() const { return m_step; }
const std::string &ScanParameters::dac() const { return m_dac; }
bool ScanParameters::enabled() const { return m_enabled; }
RawMasterFile::RawMasterFile(const std::filesystem::path &fpath)
: m_fnc(fpath) {
if (!std::filesystem::exists(fpath)) {
@ -163,10 +157,8 @@ ScanParameters RawMasterFile::scan_parameters() const {
return m_scan_parameters;
}
std::optional<ROI> RawMasterFile::roi() const { return m_roi; }
void RawMasterFile::parse_json(const std::filesystem::path &fpath) {
std::ifstream ifs(fpath);
json j;
@ -205,7 +197,6 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) {
// keep the optional empty
}
// ----------------------------------------------------------------
// Special treatment of analog flag because of Moench03
try {
@ -261,13 +252,13 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) {
std::string scan_parameters = j.at("Scan Parameters");
m_scan_parameters = ScanParameters(scan_parameters);
if (v < 7.21) {
m_scan_parameters.increment_stop(); //adjust for endpoint being included
m_scan_parameters
.increment_stop(); // adjust for endpoint being included
}
} catch (const json::out_of_range &e) {
// not a scan
}
try {
ROI tmp_roi;
auto obj = j.at("Receiver Roi");
@ -288,19 +279,14 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) {
m_roi = tmp_roi;
}
} catch (const json::out_of_range &e) {
// leave the optional empty
}
// if we have an roi we need to update the geometry for the subfiles
if (m_roi) {
}
// Update detector type for Moench
// TODO! How does this work with old .raw master files?
#ifdef AARE_VERBOSE
@ -410,7 +396,6 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) {
m_type = DetectorType::Moench03_old;
}
// TODO! Look for d0, d1...dn and update geometry
if (m_geometry.col == 0 && m_geometry.row == 0) {
m_geometry = {1, 1};

25
src/RawMasterFile.test.cpp
View File

@ -1,11 +1,10 @@
#include "aare/RawMasterFile.hpp"
#include <catch2/catch_test_macros.hpp>
#include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
using namespace aare;
TEST_CASE("Parse a master file fname") {
RawFileNameComponents m("test_master_1.json");
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"));
}
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.dac() == "dac 4");
REQUIRE(s.start() == 500);
@ -67,9 +65,9 @@ TEST_CASE("A disabled scan"){
REQUIRE(s.step() == 0);
}
TEST_CASE("Parse a master file in .json format", "[.integration]") {
auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
auto fpath =
test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
REQUIRE(std::filesystem::exists(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
REQUIRE(f.bitdepth() == 16);
// "Frame Discard Policy": "nodiscard",
// "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.digital_samples());
}
TEST_CASE("Parse a master file in .raw format", "[.integration]") {
auto fpath = test_data_path() / "moench/moench04_noise_200V_sto_both_100us_no_light_thresh_900_master_0.raw";
auto fpath =
test_data_path() /
"moench/"
"moench04_noise_200V_sto_both_100us_no_light_thresh_900_master_0.raw";
REQUIRE(std::filesystem::exists(fpath));
RawMasterFile f(fpath);
@ -209,11 +208,8 @@ TEST_CASE("Parse a master file in .raw format", "[.integration]"){
// Detector Type : 1 byte
// Header Version : 1 byte
// Packets Caught Mask : 64 bytes
}
TEST_CASE("Read eiger master file", "[.integration]") {
auto fpath = test_data_path() / "eiger" / "eiger_500k_32bit_master_0.json";
REQUIRE(std::filesystem::exists(fpath));
@ -282,7 +278,4 @@ REQUIRE(f.frame_padding() == 1);
// "Packets Caught Mask": "64 bytes"
// }
// }
}

43
src/RawSubFile.cpp
View File

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

27
src/RawSubFile.test.cpp
View File

@ -1,13 +1,14 @@
#include "aare/RawSubFile.hpp"
#include "aare/File.hpp"
#include "aare/NDArray.hpp"
#include <catch2/catch_test_macros.hpp>
#include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
using namespace aare;
TEST_CASE("Read frames directly from a RawSubFile", "[.files]") {
auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f0_0.raw";
auto fpath_raw =
test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f0_0.raw";
REQUIRE(std::filesystem::exists(fpath_raw));
RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16);
@ -17,8 +18,8 @@ TEST_CASE("Read frames directly from a RawSubFile", "[.files]"){
REQUIRE(f.bytes_per_frame() == 512 * 1024 * 2);
REQUIRE(f.bytes_per_pixel() == 2);
auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
auto fpath_npy =
test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
REQUIRE(std::filesystem::exists(fpath_npy));
// Numpy file with the same data to use as reference
@ -27,9 +28,9 @@ TEST_CASE("Read frames directly from a RawSubFile", "[.files]"){
CHECK(f.frames_in_file() == 10);
CHECK(npy.total_frames() == 10);
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) {
CHECK(f.tell() == i);
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
// f0 1,2,3
// f1 4,5,6 <-- starting here
// f2 7,8,9
// f3 10
auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f1_0.raw";
auto fpath_raw =
test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f1_0.raw";
REQUIRE(std::filesystem::exists(fpath_raw));
RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16);
auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
auto fpath_npy =
test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
REQUIRE(std::filesystem::exists(fpath_npy));
// Numpy file with the same data to use as reference
@ -61,9 +64,9 @@ TEST_CASE("Read frames directly from a RawSubFile starting at the second file",
CHECK(f.frames_in_file() == 7);
CHECK(npy.total_frames() == 10);
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) {
CHECK(f.tell() == i);
f.read_into(image.buffer(), &header);

5
src/algorithm.test.cpp
View File

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

27
src/decode.cpp
View File

@ -21,9 +21,11 @@ uint16_t adc_sar_05_decode64to16(uint64_t input){
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()) {
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++) {
@ -52,9 +54,11 @@ uint16_t adc_sar_04_decode64to16(uint64_t input){
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()) {
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 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) {
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;
@ -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);
}
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()) {
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
@ -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
for (ssize_t i = 0; i < input.shape(0); i++) {
double result = 0.0;
for (size_t bit_index = 0; bit_index < weights_powers.size(); ++bit_index) {
for (size_t bit_index = 0; bit_index < weights_powers.size();
++bit_index) {
result += ((input(i) >> bit_index) & 1) * weights_powers[bit_index];
}
output(i) = result;
}
}
} // namespace aare

11
src/decode.test.cpp
View File

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

12
src/defs.cpp
View File

@ -5,15 +5,11 @@
#include <fmt/core.h>
namespace aare {
void assert_failed(const std::string &msg)
{
void assert_failed(const std::string &msg) {
fmt::print(msg);
exit(1);
}
/**
* @brief Convert a DetectorType to a string
* @param type DetectorType
@ -102,7 +98,8 @@ template <> TimingMode StringTo(const std::string &arg) {
return TimingMode::Auto;
if (arg == "trigger")
return TimingMode::Trigger;
throw std::runtime_error("Could not decode timing mode from: \"" + arg + "\"");
throw std::runtime_error("Could not decode timing mode from: \"" + arg +
"\"");
}
template <> FrameDiscardPolicy StringTo(const std::string &arg) {
@ -112,7 +109,8 @@ template <> FrameDiscardPolicy StringTo(const std::string &arg) {
return FrameDiscardPolicy::Discard;
if (arg == "discardpartial")
return FrameDiscardPolicy::DiscardPartial;
throw std::runtime_error("Could not decode frame discard policy from: \"" + arg + "\"");
throw std::runtime_error("Could not decode frame discard policy from: \"" +
arg + "\"");
}
// template <> TimingMode StringTo<TimingMode>(std::string mode);

45
src/defs.test.cpp
View File

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

13
src/geo_helpers.cpp
View File

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

20
src/geo_helpers.test.cpp
View File

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

6
src/utils/task.test.cpp
View File

@ -1,8 +1,7 @@
#include "aare/utils/task.hpp"
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
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].second == i + 1);
}
}

2
tests/test.cpp
View File

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