detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2026-02-03 01:48:38 +01:00
Code Issues Actions Packages Releases Activity

merged developer

Browse Source
This commit is contained in:
froejdh_e
2025-04-25 11:05:09 +02:00
parent febefd92ad 2e0424254c
commit 78406dc881
67 changed files with 2763 additions and 517 deletions
Download Patch File Download Diff File Expand all files Collapse all files

450
src/ClusterFile.cpp
View File

@@ -31,6 +31,24 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
}
}
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);
}
void ClusterFile::set_gain_map(const NDView<double, 2> gain_map){
m_gain_map = NDArray<double, 2>(gain_map);
// Gain map is passed as ADU/keV to avoid dividing in when applying the gain
// map we invert it here
for (auto &item : m_gain_map->view()) {
item = 1.0 / item;
}
}
ClusterFile::~ClusterFile() { close(); }
void ClusterFile::close() {
@@ -48,14 +66,37 @@ void ClusterFile::write_frame(const ClusterVector<int32_t> &clusters) {
!(clusters.cluster_size_y() == 3)) {
throw std::runtime_error("Only 3x3 clusters are supported");
}
//First write the frame number - 4 bytes
int32_t frame_number = clusters.frame_number();
fwrite(&frame_number, sizeof(frame_number), 1, fp);
if(fwrite(&frame_number, sizeof(frame_number), 1, fp)!=1){
throw std::runtime_error(LOCATION + "Could not write frame number");
}
//Then write the number of clusters - 4 bytes
uint32_t n_clusters = clusters.size();
fwrite(&n_clusters, sizeof(n_clusters), 1, fp);
fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp);
if(fwrite(&n_clusters, sizeof(n_clusters), 1, fp)!=1){
throw std::runtime_error(LOCATION + "Could not write number of clusters");
}
//Now write the clusters in the frame
if(fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp)!=clusters.size()){
throw std::runtime_error(LOCATION + "Could not write clusters");
}
}
ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters){
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
if (m_noise_map || m_roi){
return read_clusters_with_cut(n_clusters);
}else{
return read_clusters_without_cut(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");
}
@@ -86,6 +127,7 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
if (nph_read < n_clusters) {
// keep on reading frames and photons until reaching n_clusters
while (fread(&iframe, sizeof(iframe), 1, fp)) {
clusters.set_frame_number(iframe);
// read number of clusters in frame
if (fread(&nph, sizeof(nph), 1, fp)) {
if (nph > (n_clusters - nph_read))
@@ -105,10 +147,111 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
// Resize the vector to the number of clusters.
// No new allocation, only change bounds.
clusters.resize(nph_read);
if(m_gain_map)
clusters.apply_gain_map(m_gain_map->view());
return clusters;
}
ClusterVector<int32_t> ClusterFile::read_frame() {
ClusterVector<int32_t> ClusterFile::read_clusters_with_cut(size_t n_clusters) {
ClusterVector<int32_t> clusters(3,3);
clusters.reserve(n_clusters);
// if there are photons left from previous frame read them first
if (m_num_left) {
while(m_num_left && clusters.size() < n_clusters){
Cluster3x3 c = read_one_cluster();
if(is_selected(c)){
clusters.push_back(c.x, c.y, reinterpret_cast<std::byte*>(c.data));
}
}
}
// we did not have enough clusters left in the previous frame
// keep on reading frames until reaching 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");
}
int32_t frame_number = 0; // frame number needs to be 4 bytes!
while (fread(&frame_number, sizeof(frame_number), 1, fp)) {
if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) {
clusters.set_frame_number(frame_number); //cluster vector will hold the last frame number
while(m_num_left && clusters.size() < n_clusters){
Cluster3x3 c = read_one_cluster();
if(is_selected(c)){
clusters.push_back(c.x, c.y, reinterpret_cast<std::byte*>(c.data));
}
}
}
// we have enough clusters, break out of the outer while loop
if (clusters.size() >= n_clusters)
break;
}
}
if(m_gain_map)
clusters.apply_gain_map(m_gain_map->view());
return clusters;
}
Cluster3x3 ClusterFile::read_one_cluster(){
Cluster3x3 c;
auto rc = fread(&c, sizeof(c), 1, fp);
if (rc != 1) {
throw std::runtime_error(LOCATION + "Could not read cluster");
}
--m_num_left;
return c;
}
ClusterVector<int32_t> ClusterFile::read_frame(){
if (m_mode != "r") {
throw std::runtime_error(LOCATION + "File not opened for reading");
}
if (m_noise_map || m_roi){
return read_frame_with_cut();
}else{
return read_frame_without_cut();
}
}
ClusterVector<int32_t> ClusterFile::read_frame_without_cut() {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
if (m_num_left) {
throw std::runtime_error(
"There are still photons left in the last frame");
}
int32_t frame_number;
if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
throw std::runtime_error(LOCATION + "Could not read frame number");
}
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");
}
ClusterVector<int32_t> clusters(3, 3, n_clusters);
clusters.set_frame_number(frame_number);
if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
static_cast<size_t>(n_clusters)) {
throw std::runtime_error(LOCATION + "Could not read clusters");
}
clusters.resize(n_clusters);
if (m_gain_map)
clusters.apply_gain_map(m_gain_map->view());
return clusters;
}
ClusterVector<int32_t> ClusterFile::read_frame_with_cut() {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
@@ -121,158 +264,68 @@ ClusterVector<int32_t> ClusterFile::read_frame() {
throw std::runtime_error("Could not read frame number");
}
int32_t n_clusters; // Saved as 32bit integer in the cluster file
if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) {
throw std::runtime_error("Could not read number of clusters");
}
// std::vector<Cluster3x3> clusters(n_clusters);
ClusterVector<int32_t> clusters(3, 3, n_clusters);
ClusterVector<int32_t> clusters(3, 3);
clusters.reserve(m_num_left);
clusters.set_frame_number(frame_number);
if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
static_cast<size_t>(n_clusters)) {
throw std::runtime_error("Could not read clusters");
while(m_num_left){
Cluster3x3 c = read_one_cluster();
if(is_selected(c)){
clusters.push_back(c.x, c.y, reinterpret_cast<std::byte*>(c.data));
}
}
clusters.resize(n_clusters);
if (m_gain_map)
clusters.apply_gain_map(m_gain_map->view());
return clusters;
}
// std::vector<Cluster3x3> ClusterFile::read_cluster_with_cut(size_t n_clusters,
// double *noise_map,
// int nx, int ny) {
// if (m_mode != "r") {
// throw std::runtime_error("File not opened for reading");
// }
// std::vector<Cluster3x3> clusters(n_clusters);
// // size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
// // uint32_t *n_left, double *noise_map, int
// // nx, int ny) {
// int iframe = 0;
// // uint32_t nph = *n_left;
// uint32_t nph = m_num_left;
// // uint32_t nn = *n_left;
// uint32_t nn = m_num_left;
// size_t nph_read = 0;
// int32_t t2max, tot1;
// int32_t tot3;
// // Cluster *ptr = buf;
// Cluster3x3 *ptr = clusters.data();
// int good = 1;
// double noise;
// // read photons left from previous frame
// if (noise_map)
// printf("Using noise map\n");
bool ClusterFile::is_selected(Cluster3x3 &cl) {
//Should fail fast
if (m_roi) {
if (!(m_roi->contains(cl.x, cl.y))) {
return false;
}
}
if (m_noise_map){
int32_t sum_1x1 = cl.data[4]; // central pixel
int32_t sum_2x2 = cl.sum_2x2(); // highest sum of 2x2 subclusters
int32_t sum_3x3 = cl.sum(); // sum of all pixels
// if (nph) {
// if (nph > n_clusters) {
// // if we have more photons left in the frame then photons to
// // read we read directly the requested number
// nn = n_clusters;
// } else {
// nn = nph;
// }
// for (size_t iph = 0; iph < nn; iph++) {
// // read photons 1 by 1
// size_t n_read =
// fread(reinterpret_cast<void *>(ptr), sizeof(Cluster3x3), 1, fp);
// if (n_read != 1) {
// clusters.resize(nph_read);
// return clusters;
// }
// // TODO! error handling on read
// good = 1;
// if (noise_map) {
// if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) {
// tot1 = ptr->data[4];
// analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL,
// NULL);
// noise = noise_map[ptr->y * nx + ptr->x];
// if (tot1 > noise || t2max > 2 * noise || tot3 > 3 * noise) {
// ;
// } else {
// good = 0;
// printf("%d %d %f %d %d %d\n", ptr->x, ptr->y, noise,
// tot1, t2max, tot3);
// }
// } else {
// printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
// good = 0;
// }
// }
// if (good) {
// ptr++;
// nph_read++;
// }
// (m_num_left)--;
// if (nph_read >= n_clusters)
// break;
// }
// }
// if (nph_read < n_clusters) {
// // // keep on reading frames and photons until reaching
// // n_clusters
// while (fread(&iframe, sizeof(iframe), 1, fp)) {
// // // printf("%d\n",nph_read);
// if (fread(&nph, sizeof(nph), 1, fp)) {
// // // printf("** %d\n",nph);
// m_num_left = nph;
// for (size_t iph = 0; iph < nph; iph++) {
// // // read photons 1 by 1
// size_t n_read = fread(reinterpret_cast<void *>(ptr),
// sizeof(Cluster3x3), 1, fp);
// if (n_read != 1) {
// clusters.resize(nph_read);
// return clusters;
// // return nph_read;
// }
// good = 1;
// if (noise_map) {
// if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 &&
// ptr->y < ny) {
// tot1 = ptr->data[4];
// analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL,
// NULL, NULL, NULL);
// // noise = noise_map[ptr->y * nx + ptr->x];
// noise = noise_map[ptr->y + ny * ptr->x];
// if (tot1 > noise || t2max > 2 * noise ||
// tot3 > 3 * noise) {
// ;
// } else
// good = 0;
// } else {
// printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
// good = 0;
// }
// }
// if (good) {
// ptr++;
// nph_read++;
// }
// (m_num_left)--;
// if (nph_read >= n_clusters)
// break;
// }
// }
// if (nph_read >= n_clusters)
// break;
// }
// }
// // printf("%d\n",nph_read);
// clusters.resize(nph_read);
// return clusters;
// }
auto noise = (*m_noise_map)(cl.y, cl.x); //TODO! check if this is correct
if (sum_1x1 <= noise || sum_2x2 <= 2 * noise || sum_3x3 <= 3 * noise) {
return false;
}
}
//we passed all checks
return true;
}
NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
//TOTO! make work with 2x2 clusters
NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2(clusters.at<Cluster3x3>(i));
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2(clusters.at<Cluster3x3>(i));
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
}
}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;
eta2(i, 1) = e.y;
}
}else{
throw std::runtime_error("Only 3x3 and 2x2 clusters are supported");
}
return eta2;
}
@@ -290,7 +343,7 @@ Eta2 calculate_eta2(Cluster3x3 &cl) {
tot2[3] = cl.data[4] + cl.data[5] + cl.data[7] + cl.data[8];
auto c = std::max_element(tot2.begin(), tot2.end()) - tot2.begin();
eta.sum = tot2[c];
switch (c) {
case cBottomLeft:
if ((cl.data[3] + cl.data[4]) != 0)
@@ -333,110 +386,17 @@ Eta2 calculate_eta2(Cluster3x3 &cl) {
return eta;
}
int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x,
double *eta3y) {
return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y);
Eta2 calculate_eta2(Cluster2x2 &cl) {
Eta2 eta{};
if ((cl.data[0] + cl.data[1]) != 0)
eta.x = static_cast<double>(cl.data[1]) / (cl.data[0] + cl.data[1]);
if ((cl.data[0] + cl.data[2]) != 0)
eta.y = static_cast<double>(cl.data[2]) / (cl.data[0] + cl.data[2]);
eta.sum = cl.data[0] + cl.data[1] + cl.data[2]+ cl.data[3];
eta.c = cBottomLeft; //TODO! This is not correct, but need to put something
return eta;
}
int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x, double *eta3y) {
int ok = 1;
int32_t tot2[4];
int32_t t2max = 0;
char c = 0;
int32_t val, tot3;
tot3 = 0;
for (int i = 0; i < 4; i++)
tot2[i] = 0;
for (int ix = 0; ix < 3; ix++) {
for (int iy = 0; iy < 3; iy++) {
val = data[iy * 3 + ix];
// printf ("%d ",data[iy * 3 + ix]);
tot3 += val;
if (ix <= 1 && iy <= 1)
tot2[cBottomLeft] += val;
if (ix >= 1 && iy <= 1)
tot2[cBottomRight] += val;
if (ix <= 1 && iy >= 1)
tot2[cTopLeft] += val;
if (ix >= 1 && iy >= 1)
tot2[cTopRight] += val;
}
// printf ("\n");
}
// printf ("\n");
if (t2 || quad) {
t2max = tot2[0];
c = cBottomLeft;
for (int i = 1; i < 4; i++) {
if (tot2[i] > t2max) {
t2max = tot2[i];
c = i;
}
}
// printf("*** %d %d %d %d --
// %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max);
if (quad)
*quad = c;
if (t2)
*t2 = t2max;
}
if (t3)
*t3 = tot3;
if (eta2x || eta2y) {
if (eta2x)
*eta2x = 0;
if (eta2y)
*eta2y = 0;
switch (c) {
case cBottomLeft:
if (eta2x && (data[3] + data[4]) != 0)
*eta2x = static_cast<double>(data[4]) / (data[3] + data[4]);
if (eta2y && (data[1] + data[4]) != 0)
*eta2y = static_cast<double>(data[4]) / (data[1] + data[4]);
break;
case cBottomRight:
if (eta2x && (data[2] + data[5]) != 0)
*eta2x = static_cast<double>(data[5]) / (data[4] + data[5]);
if (eta2y && (data[1] + data[4]) != 0)
*eta2y = static_cast<double>(data[4]) / (data[1] + data[4]);
break;
case cTopLeft:
if (eta2x && (data[7] + data[4]) != 0)
*eta2x = static_cast<double>(data[4]) / (data[3] + data[4]);
if (eta2y && (data[7] + data[4]) != 0)
*eta2y = static_cast<double>(data[7]) / (data[7] + data[4]);
break;
case cTopRight:
if (eta2x && t2max != 0)
*eta2x = static_cast<double>(data[5]) / (data[5] + data[4]);
if (eta2y && t2max != 0)
*eta2y = static_cast<double>(data[7]) / (data[7] + data[4]);
break;
default:;
}
}
if (eta3x || eta3y) {
if (eta3x && (data[3] + data[4] + data[5]) != 0)
*eta3x = static_cast<double>(-data[3] + data[3 + 2]) /
(data[3] + data[4] + data[5]);
if (eta3y && (data[1] + data[4] + data[7]) != 0)
*eta3y = static_cast<double>(-data[1] + data[2 * 3 + 1]) /
(data[1] + data[4] + data[7]);
}
return ok;
}
} // namespace aare

84
src/ClusterFile.test.cpp Normal file
View File

@@ -0,0 +1,84 @@
#include "aare/ClusterFile.hpp"
#include "test_config.hpp"
#include "aare/defs.hpp"
#include <catch2/catch_test_macros.hpp>
#include <filesystem>
using aare::ClusterFile;
TEST_CASE("Read one frame from a a cluster file", "[.files]") {
//We know that the frame has 97 clusters
auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
REQUIRE(std::filesystem::exists(fpath));
ClusterFile f(fpath);
auto clusters = f.read_frame();
REQUIRE(clusters.size() == 97);
REQUIRE(clusters.frame_number() == 135);
}
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";
REQUIRE(std::filesystem::exists(fpath));
ClusterFile f(fpath);
aare::ROI roi;
roi.xmin = 0;
roi.xmax = 50;
roi.ymin = 200;
roi.ymax = 249;
f.set_roi(roi);
auto clusters = f.read_frame();
REQUIRE(clusters.size() == 49);
REQUIRE(clusters.frame_number() == 135);
//Check that all clusters are within the ROI
for (size_t i = 0; i < clusters.size(); i++) {
auto c = clusters.at<aare::Cluster3x3>(i);
REQUIRE(c.x >= roi.xmin);
REQUIRE(c.x <= roi.xmax);
REQUIRE(c.y >= roi.ymin);
REQUIRE(c.y <= roi.ymax);
}
}
TEST_CASE("Read clusters from single frame file", "[.files]") {
auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
REQUIRE(std::filesystem::exists(fpath));
SECTION("Read fewer clusters than available") {
ClusterFile f(fpath);
auto clusters = f.read_clusters(50);
REQUIRE(clusters.size() == 50);
REQUIRE(clusters.frame_number() == 135);
}
SECTION("Read more clusters than available") {
ClusterFile f(fpath);
// 100 is the maximum number of clusters read
auto clusters = f.read_clusters(100);
REQUIRE(clusters.size() == 97);
REQUIRE(clusters.frame_number() == 135);
}
SECTION("Read all clusters") {
ClusterFile f(fpath);
auto clusters = f.read_clusters(97);
REQUIRE(clusters.size() == 97);
REQUIRE(clusters.frame_number() == 135);
}
}

2
src/Dtype.cpp
View File

@@ -70,7 +70,7 @@ uint8_t Dtype::bitdepth() const {
/**
* @brief Get the number of bytes of the data type
*/
size_t Dtype::bytes() const { return bitdepth() / 8; }
size_t Dtype::bytes() const { return bitdepth() / bits_per_byte; }
/**
* @brief Construct a DType object from a TypeIndex

5
src/File.cpp
View File

@@ -1,4 +1,5 @@
#include "aare/File.hpp"
#include "aare/JungfrauDataFile.hpp"
#include "aare/NumpyFile.hpp"
#include "aare/RawFile.hpp"
@@ -27,6 +28,8 @@ File::File(const std::filesystem::path &fname, const std::string &mode,
else if (fname.extension() == ".npy") {
// file_impl = new NumpyFile(fname, mode, cfg);
file_impl = std::make_unique<NumpyFile>(fname, mode, cfg);
}else if(fname.extension() == ".dat"){
file_impl = std::make_unique<JungfrauDataFile>(fname);
} else {
throw std::runtime_error("Unsupported file type");
}
@@ -73,7 +76,7 @@ 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() / 8; }
size_t File::bytes_per_pixel() const { return file_impl->bitdepth() / bits_per_byte; }
DetectorType File::detector_type() const { return file_impl->detector_type(); }

44
src/FilePtr.cpp Normal file
View File

@@ -0,0 +1,44 @@
#include "aare/FilePtr.hpp"
#include <fmt/format.h>
#include <stdexcept>
#include <utility>
namespace aare {
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()));
}
FilePtr::FilePtr(FilePtr &&other) { std::swap(fp_, other.fp_); }
FilePtr &FilePtr::operator=(FilePtr &&other) {
std::swap(fp_, other.fp_);
return *this;
}
FILE *FilePtr::get() { return fp_; }
int64_t FilePtr::tell() {
auto pos = ftell(fp_);
if (pos == -1)
throw std::runtime_error(fmt::format("Error getting file position: {}", error_msg()));
return pos;
}
FilePtr::~FilePtr() {
if (fp_)
fclose(fp_); // check?
}
std::string FilePtr::error_msg(){
if (feof(fp_)) {
return "End of file reached";
}
if (ferror(fp_)) {
return fmt::format("Error reading file: {}", std::strerror(errno));
}
return "";
}
} // namespace aare

8
src/Fit.cpp
View File

@@ -18,7 +18,7 @@ double gaus(const double x, const double *par) {
NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par) {
NDArray<double, 1> y({x.shape(0)}, 0);
for (size_t i = 0; i < x.size(); i++) {
for (ssize_t i = 0; i < x.size(); i++) {
y(i) = gaus(x(i), par.data());
}
return y;
@@ -28,7 +28,7 @@ double pol1(const double x, const double *par) { return par[0] * x + par[1]; }
NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
NDArray<double, 1> y({x.shape()}, 0);
for (size_t i = 0; i < x.size(); i++) {
for (ssize_t i = 0; i < x.size(); i++) {
y(i) = pol1(x(i), par.data());
}
return y;
@@ -177,7 +177,7 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
// Calculate chi2
chi2 = 0;
for (size_t i = 0; i < y.size(); i++) {
for (ssize_t i = 0; i < y.size(); i++) {
chi2 += std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2);
}
}
@@ -229,7 +229,7 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
// Calculate chi2
chi2 = 0;
for (size_t i = 0; i < y.size(); i++) {
for (ssize_t i = 0; i < y.size(); i++) {
chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}

139
src/Interpolator.cpp Normal file
View File

@@ -0,0 +1,139 @@
#include "aare/Interpolator.hpp"
#include "aare/algorithm.hpp"
namespace aare {
Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
NDView<double, 1> ybins, NDView<double, 1> ebins)
: m_ietax(etacube), m_ietay(etacube), m_etabinsx(xbins), m_etabinsy(ybins), m_energy_bins(ebins) {
if (etacube.shape(0) != xbins.size() || etacube.shape(1) != ybins.size() ||
etacube.shape(2) != ebins.size()) {
throw std::invalid_argument(
"The shape of the etacube does not match the shape of the bins");
}
// Cumulative sum in the x direction
for (ssize_t i = 1; i < m_ietax.shape(0); i++) {
for (ssize_t j = 0; j < m_ietax.shape(1); j++) {
for (ssize_t k = 0; k < m_ietax.shape(2); k++) {
m_ietax(i, j, k) += m_ietax(i - 1, j, k);
}
}
}
// Normalize by the highest row, if norm less than 1 don't do anything
for (ssize_t i = 0; i < m_ietax.shape(0); i++) {
for (ssize_t j = 0; j < m_ietax.shape(1); j++) {
for (ssize_t k = 0; k < m_ietax.shape(2); k++) {
auto val = m_ietax(m_ietax.shape(0) - 1, j, k);
double norm = val < 1 ? 1 : val;
m_ietax(i, j, k) /= norm;
}
}
}
// Cumulative sum in the y direction
for (ssize_t i = 0; i < m_ietay.shape(0); i++) {
for (ssize_t j = 1; j < m_ietay.shape(1); j++) {
for (ssize_t k = 0; k < m_ietay.shape(2); k++) {
m_ietay(i, j, k) += m_ietay(i, j - 1, k);
}
}
}
// Normalize by the highest column, if norm less than 1 don't do anything
for (ssize_t i = 0; i < m_ietay.shape(0); i++) {
for (ssize_t j = 0; j < m_ietay.shape(1); j++) {
for (ssize_t k = 0; k < m_ietay.shape(2); k++) {
auto val = m_ietay(i, m_ietay.shape(1) - 1, k);
double norm = val < 1 ? 1 : val;
m_ietay(i, j, k) /= norm;
}
}
}
}
std::vector<Photon> Interpolator::interpolate(const ClusterVector<int32_t>& clusters) {
std::vector<Photon> photons;
photons.reserve(clusters.size());
if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
for (size_t i = 0; i<clusters.size(); i++){
auto cluster = clusters.at<Cluster3x3>(i);
Eta2 eta= calculate_eta2(cluster);
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = eta.sum;
//Finding the index of the last element that is smaller
//should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, photon.energy);
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
double dX{}, dY{};
// cBottomLeft = 0,
// cBottomRight = 1,
// cTopLeft = 2,
// cTopRight = 3
switch (eta.c) {
case cTopLeft:
dX = -1.;
dY = 0.;
break;
case cTopRight:;
dX = 0.;
dY = 0.;
break;
case cBottomLeft:
dX = -1.;
dY = -1.;
break;
case cBottomRight:
dX = 0.;
dY = -1.;
break;
}
photon.x += m_ietax(ix, iy, ie)*2 + dX;
photon.y += m_ietay(ix, iy, ie)*2 + dY;
photons.push_back(photon);
}
}else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){
for (size_t i = 0; i<clusters.size(); i++){
auto cluster = clusters.at<Cluster2x2>(i);
Eta2 eta= calculate_eta2(cluster);
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = eta.sum;
//Now do some actual interpolation.
//Find which energy bin the cluster is in
// auto ie = nearest_index(m_energy_bins, photon.energy)-1;
// auto ix = nearest_index(m_etabinsx, eta.x)-1;
// auto iy = nearest_index(m_etabinsy, eta.y)-1;
//Finding the index of the last element that is smaller
//should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, photon.energy);
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
photon.x += m_ietax(ix, iy, ie)*2; //eta goes between 0 and 1 but we could move the hit anywhere in the 2x2
photon.y += m_ietay(ix, iy, ie)*2;
photons.push_back(photon);
}
}else{
throw std::runtime_error("Only 3x3 and 2x2 clusters are supported for interpolation");
}
return photons;
}
} // namespace aare

238
src/JungfrauDataFile.cpp Normal file
View File

@@ -0,0 +1,238 @@
#include "aare/JungfrauDataFile.hpp"
#include "aare/algorithm.hpp"
#include "aare/defs.hpp"
#include <cerrno>
#include <fmt/format.h>
namespace aare {
JungfrauDataFile::JungfrauDataFile(const std::filesystem::path &fname) {
if (!std::filesystem::exists(fname)) {
throw std::runtime_error(LOCATION +
"File does not exist: " + fname.string());
}
find_frame_size(fname);
parse_fname(fname);
scan_files();
open_file(m_current_file_index);
}
// FileInterface
Frame JungfrauDataFile::read_frame(){
Frame f(rows(), cols(), Dtype::UINT16);
read_into(reinterpret_cast<std::byte *>(f.data()), nullptr);
return f;
}
Frame JungfrauDataFile::read_frame(size_t frame_number){
seek(frame_number);
Frame f(rows(), cols(), Dtype::UINT16);
read_into(reinterpret_cast<std::byte *>(f.data()), nullptr);
return f;
}
std::vector<Frame> JungfrauDataFile::read_n(size_t n_frames) {
std::vector<Frame> frames;
for(size_t i = 0; i < n_frames; ++i){
frames.push_back(read_frame());
}
return frames;
}
void JungfrauDataFile::read_into(std::byte *image_buf) {
read_into(image_buf, nullptr);
}
void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames) {
read_into(image_buf, n_frames, nullptr);
}
size_t JungfrauDataFile::frame_number(size_t frame_index) {
seek(frame_index);
return read_header().framenum;
}
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; }
std::string JungfrauDataFile::base_name() const { return m_base_name; }
size_t JungfrauDataFile::bytes_per_frame() { return m_bytes_per_frame; }
size_t JungfrauDataFile::pixels_per_frame() { return m_rows * m_cols; }
size_t JungfrauDataFile::bytes_per_pixel() const { return sizeof(pixel_type); }
size_t JungfrauDataFile::bitdepth() const {
return bytes_per_pixel() * bits_per_byte;
}
void JungfrauDataFile::seek(size_t frame_index) {
if (frame_index >= m_total_frames) {
throw std::runtime_error(LOCATION + "Frame index out of range: " +
std::to_string(frame_index));
}
m_current_frame_index = frame_index;
auto file_index = first_larger(m_last_frame_in_file, frame_index);
if (file_index != m_current_file_index)
open_file(file_index);
auto frame_offset = (file_index)
? frame_index - m_last_frame_in_file[file_index - 1]
: frame_index;
auto byte_offset = frame_offset * (m_bytes_per_frame + header_size);
m_fp.seek(byte_offset);
};
size_t JungfrauDataFile::tell() { return m_current_frame_index; }
size_t JungfrauDataFile::total_frames() const { return m_total_frames; }
size_t JungfrauDataFile::rows() const { return m_rows; }
size_t JungfrauDataFile::cols() const { return m_cols; }
size_t JungfrauDataFile::n_files() const { return m_last_frame_in_file.size(); }
void JungfrauDataFile::find_frame_size(const std::filesystem::path &fname) {
static constexpr size_t module_data_size =
header_size + sizeof(pixel_type) * 512 * 1024;
static constexpr size_t half_data_size =
header_size + sizeof(pixel_type) * 256 * 1024;
static constexpr size_t chip_data_size =
header_size + sizeof(pixel_type) * 256 * 256;
auto file_size = std::filesystem::file_size(fname);
if (file_size == 0) {
throw std::runtime_error(LOCATION +
"Cannot guess frame size: file is empty");
}
if (file_size % module_data_size == 0) {
m_rows = 512;
m_cols = 1024;
m_bytes_per_frame = module_data_size - header_size;
} else if (file_size % half_data_size == 0) {
m_rows = 256;
m_cols = 1024;
m_bytes_per_frame = half_data_size - header_size;
} else if (file_size % chip_data_size == 0) {
m_rows = 256;
m_cols = 256;
m_bytes_per_frame = chip_data_size - header_size;
} else {
throw std::runtime_error(LOCATION +
"Cannot find frame size: file size is not a "
"multiple of any known frame size");
}
}
void JungfrauDataFile::parse_fname(const std::filesystem::path &fname) {
m_path = fname.parent_path();
m_base_name = fname.stem();
// find file index, then remove if from the base name
if (auto pos = m_base_name.find_last_of('_'); pos != std::string::npos) {
m_offset = std::stoul(m_base_name.substr(pos + 1));
m_base_name.erase(pos);
}
}
void JungfrauDataFile::scan_files() {
// find how many files we have and the number of frames in each file
m_last_frame_in_file.clear();
size_t file_index = m_offset;
while (std::filesystem::exists(fpath(file_index))) {
auto n_frames = std::filesystem::file_size(fpath(file_index)) /
(m_bytes_per_frame + header_size);
m_last_frame_in_file.push_back(n_frames);
++file_index;
}
// find where we need to open the next file and total number of frames
m_last_frame_in_file = cumsum(m_last_frame_in_file);
m_total_frames = m_last_frame_in_file.back();
}
void JungfrauDataFile::read_into(std::byte *image_buf,
JungfrauDataHeader *header) {
// read header if not passed nullptr
if (header) {
if (auto rc = fread(header, sizeof(JungfrauDataHeader), 1, m_fp.get());
rc != 1) {
throw std::runtime_error(
LOCATION +
"Could not read header from file:" + m_fp.error_msg());
}
} else {
m_fp.seek(header_size, SEEK_CUR);
}
// read data
if (auto rc = fread(image_buf, 1, m_bytes_per_frame, m_fp.get());
rc != m_bytes_per_frame) {
throw std::runtime_error(LOCATION + "Could not read image from file" +
m_fp.error_msg());
}
// prepare for next read
// if we are at the end of the file, open the next file
++m_current_frame_index;
if (m_current_frame_index >= m_last_frame_in_file[m_current_file_index] &&
(m_current_frame_index < m_total_frames)) {
++m_current_file_index;
open_file(m_current_file_index);
}
}
void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames,
JungfrauDataHeader *header) {
if (header) {
for (size_t i = 0; i < n_frames; ++i)
read_into(image_buf + i * m_bytes_per_frame, header + i);
}else{
for (size_t i = 0; i < n_frames; ++i)
read_into(image_buf + i * m_bytes_per_frame, nullptr);
}
}
void JungfrauDataFile::read_into(NDArray<uint16_t>* image, JungfrauDataHeader* header) {
if(image->shape()!=shape()){
throw std::runtime_error(LOCATION +
"Image shape does not match file size: " + std::to_string(rows()) + "x" + std::to_string(cols()));
}
read_into(reinterpret_cast<std::byte *>(image->data()), header);
}
JungfrauDataHeader JungfrauDataFile::read_header() {
JungfrauDataHeader header;
if (auto rc = fread(&header, 1, sizeof(header), m_fp.get());
rc != sizeof(header)) {
throw std::runtime_error(LOCATION + "Could not read header from file" +
m_fp.error_msg());
}
m_fp.seek(-header_size, SEEK_CUR);
return header;
}
void JungfrauDataFile::open_file(size_t file_index) {
// fmt::print(stderr, "Opening file: {}\n",
// fpath(file_index+m_offset).string());
m_fp = FilePtr(fpath(file_index + m_offset), "rb");
m_current_file_index = file_index;
}
std::filesystem::path JungfrauDataFile::fpath(size_t file_index) const {
auto fname = fmt::format("{}_{:0{}}.dat", m_base_name, file_index,
n_digits_in_file_index);
return m_path / fname;
}
} // namespace aare

114
src/JungfrauDataFile.test.cpp Normal file
View File

@@ -0,0 +1,114 @@
#include "aare/JungfrauDataFile.hpp"
#include <catch2/catch_test_macros.hpp>
#include "test_config.hpp"
using aare::JungfrauDataFile;
using aare::JungfrauDataHeader;
TEST_CASE("Open a Jungfrau data file", "[.files]") {
//we know we have 4 files with 7, 7, 7, and 3 frames
//firs frame number if 1 and the bunch id is frame_number**2
//so we can check the header
auto fpath = test_data_path() / "dat" / "AldoJF500k_000000.dat";
REQUIRE(std::filesystem::exists(fpath));
JungfrauDataFile f(fpath);
REQUIRE(f.rows() == 512);
REQUIRE(f.cols() == 1024);
REQUIRE(f.bytes_per_frame() == 1048576);
REQUIRE(f.pixels_per_frame() == 524288);
REQUIRE(f.bytes_per_pixel() == 2);
REQUIRE(f.bitdepth() == 16);
REQUIRE(f.base_name() == "AldoJF500k");
REQUIRE(f.n_files() == 4);
REQUIRE(f.tell() == 0);
REQUIRE(f.total_frames() == 24);
REQUIRE(f.current_file() == fpath);
//Check that the frame number and buch id is read correctly
for (size_t i = 0; i < 24; ++i) {
JungfrauDataHeader header;
aare::NDArray<uint16_t> image(f.shape());
f.read_into(&image, &header);
REQUIRE(header.framenum == i + 1);
REQUIRE(header.bunchid == (i + 1) * (i + 1));
REQUIRE(image.shape(0) == 512);
REQUIRE(image.shape(1) == 1024);
}
}
TEST_CASE("Seek in a JungfrauDataFile", "[.files]"){
auto fpath = test_data_path() / "dat" / "AldoJF65k_000000.dat";
REQUIRE(std::filesystem::exists(fpath));
JungfrauDataFile f(fpath);
//The file should have 113 frames
f.seek(19);
REQUIRE(f.tell() == 19);
auto h = f.read_header();
REQUIRE(h.framenum == 19+1);
//Reading again does not change the file pointer
auto h2 = f.read_header();
REQUIRE(h2.framenum == 19+1);
f.seek(59);
REQUIRE(f.tell() == 59);
auto h3 = f.read_header();
REQUIRE(h3.framenum == 59+1);
JungfrauDataHeader h4;
aare::NDArray<uint16_t> image(f.shape());
f.read_into(&image, &h4);
REQUIRE(h4.framenum == 59+1);
//now we should be on the next frame
REQUIRE(f.tell() == 60);
REQUIRE(f.read_header().framenum == 60+1);
REQUIRE_THROWS(f.seek(86356)); //out of range
}
TEST_CASE("Open a Jungfrau data file with non zero file index", "[.files]"){
auto fpath = test_data_path() / "dat" / "AldoJF65k_000003.dat";
REQUIRE(std::filesystem::exists(fpath));
JungfrauDataFile f(fpath);
//18 files per data file, opening the 3rd file we ignore the first 3
REQUIRE(f.total_frames() == 113-18*3);
REQUIRE(f.tell() == 0);
//Frame numbers start at 1 in the first file
REQUIRE(f.read_header().framenum == 18*3+1);
// moving relative to the third file
f.seek(5);
REQUIRE(f.read_header().framenum == 18*3+1+5);
// ignoring the first 3 files
REQUIRE(f.n_files() == 4);
REQUIRE(f.current_file().stem() == "AldoJF65k_000003");
}
TEST_CASE("Read into throws if size doesn't match", "[.files]"){
auto fpath = test_data_path() / "dat" / "AldoJF65k_000000.dat";
REQUIRE(std::filesystem::exists(fpath));
JungfrauDataFile f(fpath);
aare::NDArray<uint16_t> image({39, 85});
JungfrauDataHeader header;
REQUIRE_THROWS(f.read_into(&image, &header));
REQUIRE_THROWS(f.read_into(&image, nullptr));
REQUIRE_THROWS(f.read_into(&image));
REQUIRE(f.tell() == 0);
}

23
src/NDArray.test.cpp
View File

@@ -2,6 +2,7 @@
#include <array>
#include <catch2/benchmark/catch_benchmark.hpp>
#include <catch2/catch_test_macros.hpp>
#include <numeric>
using aare::NDArray;
using aare::NDView;
@@ -34,6 +35,24 @@ TEST_CASE("Construct from an NDView") {
}
}
TEST_CASE("3D NDArray from NDView"){
std::vector<int> data(27);
std::iota(data.begin(), data.end(), 0);
NDView<int, 3> view(data.data(), Shape<3>{3, 3, 3});
NDArray<int, 3> image(view);
REQUIRE(image.shape() == view.shape());
REQUIRE(image.size() == view.size());
REQUIRE(image.data() != view.data());
for(int64_t i=0; i<image.shape(0); i++){
for(int64_t j=0; j<image.shape(1); j++){
for(int64_t k=0; k<image.shape(2); k++){
REQUIRE(image(i, j, k) == view(i, j, k));
}
}
}
}
TEST_CASE("1D image") {
std::array<int64_t, 1> shape{{20}};
NDArray<short, 1> img(shape, 3);
@@ -164,14 +183,14 @@ TEST_CASE("Size and shape matches") {
int64_t h = 75;
std::array<int64_t, 2> shape{w, h};
NDArray<double> a{shape};
REQUIRE(a.size() == static_cast<uint64_t>(w * h));
REQUIRE(a.size() == w * h);
REQUIRE(a.shape() == shape);
}
TEST_CASE("Initial value matches for all elements") {
double v = 4.35;
NDArray<double> a{{5, 5}, v};
for (uint32_t i = 0; i < a.size(); ++i) {
for (int i = 0; i < a.size(); ++i) {
REQUIRE(a(i) == v);
}
}

47
src/NDView.test.cpp
View File

@@ -3,6 +3,7 @@
#include <iostream>
#include <vector>
#include <numeric>
using aare::NDView;
using aare::Shape;
@@ -21,10 +22,8 @@ TEST_CASE("Element reference 1D") {
}
TEST_CASE("Element reference 2D") {
std::vector<int> vec;
for (int i = 0; i != 12; ++i) {
vec.push_back(i);
}
std::vector<int> vec(12);
std::iota(vec.begin(), vec.end(), 0);
NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
REQUIRE(vec.size() == static_cast<size_t>(data.size()));
@@ -58,10 +57,8 @@ TEST_CASE("Element reference 3D") {
}
TEST_CASE("Plus and miuns with single value") {
std::vector<int> vec;
for (int i = 0; i != 12; ++i) {
vec.push_back(i);
}
std::vector<int> vec(12);
std::iota(vec.begin(), vec.end(), 0);
NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
data += 5;
int i = 0;
@@ -116,10 +113,8 @@ TEST_CASE("elementwise assign") {
}
TEST_CASE("iterators") {
std::vector<int> vec;
for (int i = 0; i != 12; ++i) {
vec.push_back(i);
}
std::vector<int> vec(12);
std::iota(vec.begin(), vec.end(), 0);
NDView<int, 1> data(vec.data(), Shape<1>{12});
int i = 0;
for (const auto item : data) {
@@ -167,27 +162,31 @@ TEST_CASE("divide with another span") {
}
TEST_CASE("Retrieve shape") {
std::vector<int> vec;
for (int i = 0; i != 12; ++i) {
vec.push_back(i);
}
std::vector<int> vec(12);
std::iota(vec.begin(), vec.end(), 0);
NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
REQUIRE(data.shape()[0] == 3);
REQUIRE(data.shape()[1] == 4);
}
TEST_CASE("compare two views") {
std::vector<int> vec1;
for (int i = 0; i != 12; ++i) {
vec1.push_back(i);
}
std::vector<int> vec1(12);
std::iota(vec1.begin(), vec1.end(), 0);
NDView<int, 2> view1(vec1.data(), Shape<2>{3, 4});
std::vector<int> vec2;
for (int i = 0; i != 12; ++i) {
vec2.push_back(i);
}
std::vector<int> vec2(12);
std::iota(vec2.begin(), vec2.end(), 0);
NDView<int, 2> view2(vec2.data(), Shape<2>{3, 4});
REQUIRE((view1 == view2));
}
TEST_CASE("Create a view over a vector"){
std::vector<int> vec(12);
std::iota(vec.begin(), vec.end(), 0);
auto v = aare::make_view(vec);
REQUIRE(v.shape()[0] == 12);
REQUIRE(v[0] == 0);
REQUIRE(v[11] == 11);
}

3
src/RawFile.cpp
View File

@@ -76,8 +76,7 @@ size_t RawFile::n_mod() const { return n_subfile_parts; }
size_t RawFile::bytes_per_frame() {
// return m_rows * m_cols * m_master.bitdepth() / 8;
return m_geometry.pixels_x * m_geometry.pixels_y * m_master.bitdepth() / 8;
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;

31
src/RawSubFile.cpp
View File

@@ -1,9 +1,12 @@
#include "aare/RawSubFile.hpp"
#include "aare/PixelMap.hpp"
#include "aare/utils/ifstream_helpers.hpp"
#include <cstring> // memcpy
#include <fmt/core.h>
#include <iostream>
namespace aare {
RawSubFile::RawSubFile(const std::filesystem::path &fname,
@@ -20,7 +23,7 @@ RawSubFile::RawSubFile(const std::filesystem::path &fname,
}
if (std::filesystem::exists(fname)) {
n_frames = std::filesystem::file_size(fname) /
m_num_frames = std::filesystem::file_size(fname) /
(sizeof(DetectorHeader) + rows * cols * bitdepth / 8);
} else {
throw std::runtime_error(
@@ -35,7 +38,7 @@ RawSubFile::RawSubFile(const std::filesystem::path &fname,
}
#ifdef AARE_VERBOSE
fmt::print("Opened file: {} with {} frames\n", m_fname.string(), n_frames);
fmt::print("Opened file: {} with {} frames\n", m_fname.string(), m_num_frames);
fmt::print("m_rows: {}, m_cols: {}, m_bitdepth: {}\n", m_rows, m_cols,
m_bitdepth);
fmt::print("file size: {}\n", std::filesystem::file_size(fname));
@@ -43,8 +46,8 @@ RawSubFile::RawSubFile(const std::filesystem::path &fname,
}
void RawSubFile::seek(size_t frame_index) {
if (frame_index >= n_frames) {
throw std::runtime_error(LOCATION + fmt::format("Frame index {} out of range in a file with {} frames", frame_index, n_frames));
if (frame_index >= m_num_frames) {
throw std::runtime_error(LOCATION + fmt::format("Frame index {} out of range in a file with {} frames", frame_index, m_num_frames));
}
m_file.seekg((sizeof(DetectorHeader) + bytes_per_frame()) * frame_index);
}
@@ -60,6 +63,10 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) {
m_file.seekg(sizeof(DetectorHeader), std::ios::cur);
}
if (m_file.fail()){
throw std::runtime_error(LOCATION + ifstream_error_msg(m_file));
}
// TODO! expand support for different bitdepths
if (m_pixel_map) {
// read into a temporary buffer and then copy the data to the buffer
@@ -79,8 +86,24 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) {
// read directly into the buffer
m_file.read(reinterpret_cast<char *>(image_buf), bytes_per_frame());
}
if (m_file.fail()){
throw std::runtime_error(LOCATION + ifstream_error_msg(m_file));
}
}
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();
if (header) {
++header;
}
}
}
template <typename T>
void RawSubFile::read_with_map(std::byte *image_buf) {
auto part_buffer = new std::byte[bytes_per_frame()];

159
src/algorithm.test.cpp Normal file
View File

@@ -0,0 +1,159 @@
#include <catch2/catch_test_macros.hpp>
#include <aare/algorithm.hpp>
TEST_CASE("Find the closed index in a 1D array", "[algorithm]") {
aare::NDArray<double, 1> arr({5});
for (ssize_t i = 0; i < arr.size(); i++) {
arr[i] = i;
}
// arr 0, 1, 2, 3, 4
REQUIRE(aare::nearest_index(arr, 2.3) == 2);
REQUIRE(aare::nearest_index(arr, 2.6) == 3);
REQUIRE(aare::nearest_index(arr, 45.0) == 4);
REQUIRE(aare::nearest_index(arr, 0.0) == 0);
REQUIRE(aare::nearest_index(arr, -1.0) == 0);
}
TEST_CASE("Passing integers to nearest_index works", "[algorithm]"){
aare::NDArray<int, 1> arr({5});
for (ssize_t i = 0; i < arr.size(); i++) {
arr[i] = i;
}
// arr 0, 1, 2, 3, 4
REQUIRE(aare::nearest_index(arr, 2) == 2);
REQUIRE(aare::nearest_index(arr, 3) == 3);
REQUIRE(aare::nearest_index(arr, 45) == 4);
REQUIRE(aare::nearest_index(arr, 0) == 0);
REQUIRE(aare::nearest_index(arr, -1) == 0);
}
TEST_CASE("nearest_index works with std::vector", "[algorithm]"){
std::vector<double> vec = {0, 1, 2, 3, 4};
REQUIRE(aare::nearest_index(vec, 2.123) == 2);
REQUIRE(aare::nearest_index(vec, 2.66) == 3);
REQUIRE(aare::nearest_index(vec, 4555555.0) == 4);
REQUIRE(aare::nearest_index(vec, 0.0) == 0);
REQUIRE(aare::nearest_index(vec, -10.0) == 0);
}
TEST_CASE("nearest index works with std::array", "[algorithm]"){
std::array<double, 5> arr = {0, 1, 2, 3, 4};
REQUIRE(aare::nearest_index(arr, 2.123) == 2);
REQUIRE(aare::nearest_index(arr, 2.501) == 3);
REQUIRE(aare::nearest_index(arr, 4555555.0) == 4);
REQUIRE(aare::nearest_index(arr, 0.0) == 0);
REQUIRE(aare::nearest_index(arr, -10.0) == 0);
}
TEST_CASE("nearest index when there is no different uses the first element", "[algorithm]"){
std::vector<int> vec = {5, 5, 5, 5, 5};
REQUIRE(aare::nearest_index(vec, 5) == 0);
}
TEST_CASE("nearest index when there is no different uses the first element also when all smaller", "[algorithm]"){
std::vector<int> vec = {5, 5, 5, 5, 5};
REQUIRE(aare::nearest_index(vec, 10) == 0);
}
TEST_CASE("last smaller", "[algorithm]"){
aare::NDArray<double, 1> arr({5});
for (ssize_t i = 0; i < arr.size(); i++) {
arr[i] = i;
}
// arr 0, 1, 2, 3, 4
REQUIRE(aare::last_smaller(arr, -10.0) == 0);
REQUIRE(aare::last_smaller(arr, 0.0) == 0);
REQUIRE(aare::last_smaller(arr, 2.3) == 2);
REQUIRE(aare::last_smaller(arr, 253.) == 4);
}
TEST_CASE("returns last bin strictly smaller", "[algorithm]"){
aare::NDArray<double, 1> arr({5});
for (ssize_t i = 0; i < arr.size(); i++) {
arr[i] = i;
}
// arr 0, 1, 2, 3, 4
REQUIRE(aare::last_smaller(arr, 2.0) == 1);
}
TEST_CASE("last_smaller with all elements smaller returns last element", "[algorithm]"){
aare::NDArray<double, 1> arr({5});
for (ssize_t i = 0; i < arr.size(); i++) {
arr[i] = i;
}
// arr 0, 1, 2, 3, 4
REQUIRE(aare::last_smaller(arr, 50.) == 4);
}
TEST_CASE("last_smaller with all elements bigger returns first element", "[algorithm]"){
aare::NDArray<double, 1> arr({5});
for (ssize_t i = 0; i < arr.size(); i++) {
arr[i] = i;
}
// arr 0, 1, 2, 3, 4
REQUIRE(aare::last_smaller(arr, -50.) == 0);
}
TEST_CASE("last smaller with all elements equal returns the first element", "[algorithm]"){
std::vector<int> vec = {5,5,5,5,5,5,5};
REQUIRE(aare::last_smaller(vec, 5) == 0);
}
TEST_CASE("first_lager with vector", "[algorithm]"){
std::vector<double> vec = {0, 1, 2, 3, 4};
REQUIRE(aare::first_larger(vec, 2.5) == 3);
}
TEST_CASE("first_lager with all elements smaller returns last element", "[algorithm]"){
std::vector<double> vec = {0, 1, 2, 3, 4};
REQUIRE(aare::first_larger(vec, 50.) == 4);
}
TEST_CASE("first_lager with all elements bigger returns first element", "[algorithm]"){
std::vector<double> vec = {0, 1, 2, 3, 4};
REQUIRE(aare::first_larger(vec, -50.) == 0);
}
TEST_CASE("first_lager with all elements the same as the check returns last", "[algorithm]"){
std::vector<int> vec = {14, 14, 14, 14, 14};
REQUIRE(aare::first_larger(vec, 14) == 4);
}
TEST_CASE("first larger with the same element", "[algorithm]"){
std::vector<int> vec = {7,8,9,10,11};
REQUIRE(aare::first_larger(vec, 9) == 3);
}
TEST_CASE("cumsum works", "[algorithm]"){
std::vector<double> vec = {0, 1, 2, 3, 4};
auto result = aare::cumsum(vec);
REQUIRE(result.size() == vec.size());
REQUIRE(result[0] == 0);
REQUIRE(result[1] == 1);
REQUIRE(result[2] == 3);
REQUIRE(result[3] == 6);
REQUIRE(result[4] == 10);
}
TEST_CASE("cumsum works with empty vector", "[algorithm]"){
std::vector<double> vec = {};
auto result = aare::cumsum(vec);
REQUIRE(result.size() == 0);
}
TEST_CASE("cumsum works with negative numbers", "[algorithm]"){
std::vector<double> vec = {0, -1, -2, -3, -4};
auto result = aare::cumsum(vec);
REQUIRE(result.size() == vec.size());
REQUIRE(result[0] == 0);
REQUIRE(result[1] == -1);
REQUIRE(result[2] == -3);
REQUIRE(result[3] == -6);
REQUIRE(result[4] == -10);
}

43
src/decode.cpp
View File

@@ -1,5 +1,5 @@
#include "aare/decode.hpp"
#include <cmath>
namespace aare {
uint16_t adc_sar_05_decode64to16(uint64_t input){
@@ -22,6 +22,10 @@ uint16_t adc_sar_05_decode64to16(uint64_t input){
}
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");
}
for(int64_t i = 0; i < input.shape(0); i++){
for(int64_t j = 0; j < input.shape(1); j++){
output(i,j) = adc_sar_05_decode64to16(input(i,j));
@@ -49,6 +53,9 @@ uint16_t adc_sar_04_decode64to16(uint64_t input){
}
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");
}
for(int64_t i = 0; i < input.shape(0); i++){
for(int64_t j = 0; j < input.shape(1); j++){
output(i,j) = adc_sar_04_decode64to16(input(i,j));
@@ -56,6 +63,40 @@ 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");
}
double result = 0.0;
for (ssize_t i = 0; i < weights.size(); ++i) {
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) {
if(input.shape() != output.shape()){
throw std::invalid_argument(LOCATION + " input and output shapes must match");
}
//Calculate weights to avoid repeatedly calling std::pow
std::vector<double> weights_powers(weights.size());
for (ssize_t i = 0; i < weights.size(); ++i) {
weights_powers[i] = std::pow(weights[i], i);
}
// Apply custom weights to each element in the input array
for (ssize_t i = 0; i < input.shape(0); i++) {
double result = 0.0;
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

80
src/decode.test.cpp Normal file
View File

@@ -0,0 +1,80 @@
#include "aare/decode.hpp"
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/catch_test_macros.hpp>
#include "aare/NDArray.hpp"
using Catch::Matchers::WithinAbs;
#include <vector>
TEST_CASE("test_adc_sar_05_decode64to16"){
uint64_t input = 0;
uint16_t output = aare::adc_sar_05_decode64to16(input);
CHECK(output == 0);
// bit 29 on th input is bit 0 on the output
input = 1UL << 29;
output = aare::adc_sar_05_decode64to16(input);
CHECK(output == 1);
// test all bits by iteratting through the bitlist
std::vector<int> bitlist = {29, 19, 28, 18, 31, 21, 27, 20, 24, 23, 25, 22};
for (size_t i = 0; i < bitlist.size(); i++) {
input = 1UL << bitlist[i];
output = aare::adc_sar_05_decode64to16(input);
CHECK(output == (1 << i));
}
// test a few "random" values
input = 0;
input |= (1UL << 29);
input |= (1UL << 19);
input |= (1UL << 28);
output = aare::adc_sar_05_decode64to16(input);
CHECK(output == 7UL);
input = 0;
input |= (1UL << 18);
input |= (1UL << 27);
input |= (1UL << 25);
output = aare::adc_sar_05_decode64to16(input);
CHECK(output == 1096UL);
input = 0;
input |= (1UL << 25);
input |= (1UL << 22);
output = aare::adc_sar_05_decode64to16(input);
CHECK(output == 3072UL);
}
TEST_CASE("test_apply_custom_weights") {
uint16_t input = 1;
aare::NDArray<double, 1> weights_data({3}, 0.0);
weights_data(0) = 1.7;
weights_data(1) = 2.1;
weights_data(2) = 1.8;
auto weights = weights_data.view();
double output = aare::apply_custom_weights(input, weights);
CHECK_THAT(output, WithinAbs(1.0, 0.001));
input = 1 << 1;
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));
input = 0b111;
output = aare::apply_custom_weights(input, weights);
CHECK_THAT(output, WithinAbs(6.34, 0.001));
}

18
src/utils/ifstream_helpers.cpp Normal file
View File

@@ -0,0 +1,18 @@
#include "aare/utils/ifstream_helpers.hpp"
namespace aare {
std::string ifstream_error_msg(std::ifstream &ifs) {
std::ios_base::iostate state = ifs.rdstate();
if (state & std::ios_base::eofbit) {
return " End of file reached";
} else if (state & std::ios_base::badbit) {
return " Bad file stream";
} else if (state & std::ios_base::failbit) {
return " File read failed";
}else{
return " Unknown/no error";
}
}
} // namespace aare