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

Added fitting, fixed roi etc (#129)

Browse Source 
Co-authored-by: Patrick <patrick.sieberer@psi.ch>
Co-authored-by: JulianHeymes <julian.heymes@psi.ch>
This commit is contained in:
Erik Fröjdh
2025-02-12 16:50:31 +01:00
committed by GitHub
parent 7d6223d52d
commit dadf5f4869
55 changed files with 2931 additions and 693 deletions
Download Patch File Download Diff File Expand all files Collapse all files

337
src/ClusterFile.cpp
View File

@ -20,6 +20,12 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
throw std::runtime_error("Could not open file for writing: " +
fname.string());
}
} else if (mode == "a") {
fp = fopen(fname.c_str(), "ab");
if (!fp) {
throw std::runtime_error("Could not open file for appending: " +
fname.string());
}
} else {
throw std::runtime_error("Unsupported mode: " + mode);
}
@ -34,35 +40,35 @@ void ClusterFile::close() {
}
}
void ClusterFile::write_frame(int32_t frame_number,
const ClusterVector<int32_t> &clusters) {
if (m_mode != "w") {
void ClusterFile::write_frame(const ClusterVector<int32_t> &clusters) {
if (m_mode != "w" && m_mode != "a") {
throw std::runtime_error("File not opened for writing");
}
if (!(clusters.cluster_size_x() == 3) &&
!(clusters.cluster_size_y() == 3)) {
throw std::runtime_error("Only 3x3 clusters are supported");
}
int32_t frame_number = clusters.frame_number();
fwrite(&frame_number, sizeof(frame_number), 1, fp);
uint32_t n_clusters = clusters.size();
fwrite(&n_clusters, sizeof(n_clusters), 1, fp);
fwrite(clusters.data(), clusters.element_offset(), clusters.size(), fp);
// write clusters
// fwrite(clusters.data(), sizeof(Cluster), clusters.size(), fp);
fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp);
}
std::vector<Cluster3x3> 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");
}
std::vector<Cluster3x3> clusters(n_clusters);
ClusterVector<int32_t> clusters(3,3, n_clusters);
int32_t iframe = 0; // frame number needs to be 4 bytes!
size_t nph_read = 0;
uint32_t nn = m_num_left;
uint32_t nph = m_num_left; // number of clusters in frame needs to be 4
auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
// auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
auto buf = clusters.data();
// if there are photons left from previous frame read them first
if (nph) {
if (nph > n_clusters) {
@ -72,8 +78,8 @@ std::vector<Cluster3x3> ClusterFile::read_clusters(size_t n_clusters) {
} else {
nn = nph;
}
nph_read += fread(reinterpret_cast<void *>(buf + nph_read),
sizeof(Cluster3x3), nn, fp);
nph_read += fread((buf + nph_read*clusters.item_size()),
clusters.item_size(), nn, fp);
m_num_left = nph - nn; // write back the number of photons left
}
@ -87,8 +93,8 @@ std::vector<Cluster3x3> ClusterFile::read_clusters(size_t n_clusters) {
else
nn = nph;
nph_read += fread(reinterpret_cast<void *>(buf + nph_read),
sizeof(Cluster3x3), nn, fp);
nph_read += fread((buf + nph_read*clusters.item_size()),
clusters.item_size(), nn, fp);
m_num_left = nph - nn;
}
if (nph_read >= n_clusters)
@ -102,7 +108,7 @@ std::vector<Cluster3x3> ClusterFile::read_clusters(size_t n_clusters) {
return clusters;
}
std::vector<Cluster3x3> ClusterFile::read_frame(int32_t &out_fnum) {
ClusterVector<int32_t> ClusterFile::read_frame() {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
@ -110,8 +116,8 @@ std::vector<Cluster3x3> ClusterFile::read_frame(int32_t &out_fnum) {
throw std::runtime_error(
"There are still photons left in the last frame");
}
if (fread(&out_fnum, sizeof(out_fnum), 1, fp) != 1) {
int32_t frame_number;
if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
throw std::runtime_error("Could not read frame number");
}
@ -119,158 +125,163 @@ std::vector<Cluster3x3> ClusterFile::read_frame(int32_t &out_fnum) {
if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
throw std::runtime_error("Could not read number of clusters");
}
std::vector<Cluster3x3> clusters(n_clusters);
// std::vector<Cluster3x3> clusters(n_clusters);
ClusterVector<int32_t> clusters(3, 3, n_clusters);
clusters.set_frame_number(frame_number);
if (fread(clusters.data(), sizeof(Cluster3x3), n_clusters, fp) !=
if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
static_cast<size_t>(n_clusters)) {
throw std::runtime_error("Could not read clusters");
}
clusters.resize(n_clusters);
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");
// 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;
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);
// 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");
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;
}
// 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;
// }
NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
NDArray<double, 2> eta2({clusters.size(), 2});
//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++) {
// int32_t t2;
// auto* ptr = reinterpret_cast<int32_t*> (clusters.element_ptr(i) + 2 *
// sizeof(int16_t)); analyze_cluster(clusters.at<Cluster3x3>(i), &t2,
// nullptr, nullptr, &eta2(i,0), &eta2(i,1) , nullptr, nullptr);
auto [x, y] = calculate_eta2(clusters.at<Cluster3x3>(i));
eta2(i, 0) = x;
eta2(i, 1) = y;
auto e = calculate_eta2(clusters.at<Cluster3x3>(i));
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
}
return eta2;
}
std::array<double, 2> calculate_eta2(Cluster3x3 &cl) {
std::array<double, 2> eta2{};
/**
* @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2 struct
* containing etay, etax and the corner of the cluster.
*/
Eta2 calculate_eta2(Cluster3x3 &cl) {
Eta2 eta{};
std::array<int32_t, 4> tot2;
tot2[0] = cl.data[0] + cl.data[1] + cl.data[3] + cl.data[4];
@ -283,39 +294,43 @@ std::array<double, 2> calculate_eta2(Cluster3x3 &cl) {
switch (c) {
case cBottomLeft:
if ((cl.data[3] + cl.data[4]) != 0)
eta2[0] =
eta.x =
static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[1] + cl.data[4]) != 0)
eta2[1] =
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)
eta2[0] =
eta.x =
static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
if ((cl.data[1] + cl.data[4]) != 0)
eta2[1] =
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)
eta2[0] =
eta.x =
static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta2[1] =
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)
eta2[0] =
eta.x =
static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta2[1] =
eta.y =
static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.c = cTopRight;
break;
// default:;
// no default to allow compiler to warn about missing cases
}
return eta2;
return eta;
}
int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad,

112
src/ClusterVector.test.cpp
View File

@ -6,12 +6,14 @@
using aare::ClusterVector;
struct Cluster_i2x2 {
int16_t x;
int16_t y;
int32_t data[4];
};
TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read") {
struct Cluster_i2x2 {
int16_t x;
int16_t y;
int32_t data[4];
};
ClusterVector<int32_t> cv(2, 2, 4);
REQUIRE(cv.capacity() == 4);
@ -19,7 +21,7 @@ TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read") {
REQUIRE(cv.cluster_size_x() == 2);
REQUIRE(cv.cluster_size_y() == 2);
// int16_t, int16_t, 2x2 int32_t = 20 bytes
REQUIRE(cv.element_offset() == 20);
REQUIRE(cv.item_size() == 20);
//Create a cluster and push back into the vector
Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
@ -30,7 +32,7 @@ TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read") {
//Read the cluster back out using copy. TODO! Can we improve the API?
Cluster_i2x2 c2;
std::byte *ptr = cv.element_ptr(0);
std::copy(ptr, ptr + cv.element_offset(), reinterpret_cast<std::byte*>(&c2));
std::copy(ptr, ptr + cv.item_size(), reinterpret_cast<std::byte*>(&c2));
//Check that the data is the same
REQUIRE(c1.x == c2.x);
@ -83,8 +85,8 @@ TEST_CASE("Storing floats"){
float data[8];
};
ClusterVector<float> cv(2, 4, 2);
REQUIRE(cv.capacity() == 2);
ClusterVector<float> cv(2, 4, 10);
REQUIRE(cv.capacity() == 10);
REQUIRE(cv.size() == 0);
REQUIRE(cv.cluster_size_x() == 2);
REQUIRE(cv.cluster_size_y() == 4);
@ -92,17 +94,105 @@ TEST_CASE("Storing floats"){
//Create a cluster and push back into the vector
Cluster_f4x2 c1 = {1, 2, {3.0, 4.0, 5.0, 6.0,3.0, 4.0, 5.0, 6.0}};
cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.capacity() == 10);
REQUIRE(cv.size() == 1);
Cluster_f4x2 c2 = {6, 7, {8.0, 9.0, 10.0, 11.0,8.0, 9.0, 10.0, 11.0}};
cv.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.capacity() == 10);
REQUIRE(cv.size() == 2);
auto sums = cv.sum();
REQUIRE(sums.size() == 2);
REQUIRE_THAT(sums[0], Catch::Matchers::WithinAbs(36.0, 1e-6));
REQUIRE_THAT(sums[1], Catch::Matchers::WithinAbs(76.0, 1e-6));
}
TEST_CASE("Push back more than initial capacity"){
ClusterVector<int32_t> cv(2, 2, 2);
auto initial_data = cv.data();
Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
REQUIRE(cv.size() == 1);
REQUIRE(cv.capacity() == 2);
Cluster_i2x2 c2 = {6, 7, {8, 9, 10, 11}};
cv.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
REQUIRE(cv.size() == 2);
REQUIRE(cv.capacity() == 2);
Cluster_i2x2 c3 = {11, 12, {13, 14, 15, 16}};
cv.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
REQUIRE(cv.size() == 3);
REQUIRE(cv.capacity() == 4);
Cluster_i2x2* ptr = reinterpret_cast<Cluster_i2x2*>(cv.data());
REQUIRE(ptr[0].x == 1);
REQUIRE(ptr[0].y == 2);
REQUIRE(ptr[1].x == 6);
REQUIRE(ptr[1].y == 7);
REQUIRE(ptr[2].x == 11);
REQUIRE(ptr[2].y == 12);
//We should have allocated a new buffer, since we outgrew the initial capacity
REQUIRE(initial_data != cv.data());
}
TEST_CASE("Concatenate two cluster vectors where the first has enough capacity"){
ClusterVector<int32_t> cv1(2, 2, 12);
Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
cv1.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
Cluster_i2x2 c2 = {6, 7, {8, 9, 10, 11}};
cv1.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
ClusterVector<int32_t> cv2(2, 2, 2);
Cluster_i2x2 c3 = {11, 12, {13, 14, 15, 16}};
cv2.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
Cluster_i2x2 c4 = {16, 17, {18, 19, 20, 21}};
cv2.push_back(c4.x, c4.y, reinterpret_cast<std::byte*>(&c4.data[0]));
cv1 += cv2;
REQUIRE(cv1.size() == 4);
REQUIRE(cv1.capacity() == 12);
Cluster_i2x2* ptr = reinterpret_cast<Cluster_i2x2*>(cv1.data());
REQUIRE(ptr[0].x == 1);
REQUIRE(ptr[0].y == 2);
REQUIRE(ptr[1].x == 6);
REQUIRE(ptr[1].y == 7);
REQUIRE(ptr[2].x == 11);
REQUIRE(ptr[2].y == 12);
REQUIRE(ptr[3].x == 16);
REQUIRE(ptr[3].y == 17);
}
TEST_CASE("Concatenate two cluster vectors where we need to allocate"){
ClusterVector<int32_t> cv1(2, 2, 2);
Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
cv1.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
Cluster_i2x2 c2 = {6, 7, {8, 9, 10, 11}};
cv1.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
ClusterVector<int32_t> cv2(2, 2, 2);
Cluster_i2x2 c3 = {11, 12, {13, 14, 15, 16}};
cv2.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
Cluster_i2x2 c4 = {16, 17, {18, 19, 20, 21}};
cv2.push_back(c4.x, c4.y, reinterpret_cast<std::byte*>(&c4.data[0]));
cv1 += cv2;
REQUIRE(cv1.size() == 4);
REQUIRE(cv1.capacity() == 4);
Cluster_i2x2* ptr = reinterpret_cast<Cluster_i2x2*>(cv1.data());
REQUIRE(ptr[0].x == 1);
REQUIRE(ptr[0].y == 2);
REQUIRE(ptr[1].x == 6);
REQUIRE(ptr[1].y == 7);
REQUIRE(ptr[2].x == 11);
REQUIRE(ptr[2].y == 12);
REQUIRE(ptr[3].x == 16);
REQUIRE(ptr[3].y == 17);
}

2
src/File.cpp
View File

@ -45,6 +45,8 @@ File& File::operator=(File &&other) noexcept {
return *this;
}
// void File::close() { file_impl->close(); }
Frame File::read_frame() { return file_impl->read_frame(); }
Frame File::read_frame(size_t frame_index) {
return file_impl->read_frame(frame_index);

300
src/Fit.cpp Normal file
View File

@ -0,0 +1,300 @@
#include "aare/Fit.hpp"
#include "aare/utils/task.hpp"
#include <lmcurve2.h>
#include <lmfit.hpp>
#include <thread>
namespace aare {
namespace func {
double gaus(const double x, const double *par) {
return par[0] * exp(-pow(x - par[1], 2) / (2 * pow(par[2], 2)));
}
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++) {
y(i) = gaus(x(i), par.data());
}
return y;
}
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++) {
y(i) = pol1(x(i), par.data());
}
return y;
}
} // namespace func
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y) {
NDArray<double, 1> result({3}, 0);
lm_control_struct control = lm_control_double;
// Estimate the initial parameters for the fit
std::vector<double> start_par{0, 0, 0};
auto e = std::max_element(y.begin(), y.end());
auto idx = std::distance(y.begin(), e);
start_par[0] = *e; // For amplitude we use the maximum value
start_par[1] =
x[idx]; // For the mean we use the x value of the maximum value
// 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(),
[e, delta](double val) { return val > *e / 2; }) *
delta / 2.35;
lmfit::result_t res(start_par);
lmcurve(res.par.size(), res.par.data(), x.size(), x.data(), y.data(),
aare::func::gaus, &control, &res.status);
result(0) = res.par[0];
result(1) = res.par[1];
result(2) = res.par[2];
return result;
}
NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 3}, 0);
auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
auto res = fit_gaus(x, values);
result(row, col, 0) = res(0);
result(row, col, 1) = res(1);
result(row, col, 2) = res(2);
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads;
for (auto &task : tasks) {
threads.push_back(std::thread(process, task.first, task.second));
}
for (auto &thread : threads) {
thread.join();
}
return result;
}
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,
int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
NDView<double, 1> y_err_view(&y_err(row, col, 0),
{y_err.shape(2)});
NDView<double, 1> par_out_view(&par_out(row, col, 0),
{par_out.shape(2)});
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)});
fit_gaus(x, y_view, y_err_view, par_out_view, par_err_out_view);
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads;
for (auto &task : tasks) {
threads.push_back(std::thread(process, task.first, task.second));
}
for (auto &thread : threads) {
thread.join();
}
}
void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out) {
// Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() ||
par_out.size() != 3 || par_err_out.size() != 3) {
throw std::runtime_error("Data, x, data_err must have the same size "
"and par_out, par_err_out must have size 3");
}
lm_control_struct control = lm_control_double;
// Estimate the initial parameters for the fit
std::vector<double> start_par{0, 0, 0};
std::vector<double> start_par_err{0, 0, 0};
std::vector<double> start_cov{0, 0, 0, 0, 0, 0, 0, 0, 0};
auto e = std::max_element(y.begin(), y.end());
auto idx = std::distance(y.begin(), e);
start_par[0] = *e; // For amplitude we use the maximum value
start_par[1] =
x[idx]; // For the mean we use the x value of the maximum value
// 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(),
[e, delta](double val) { return val > *e / 2; }) *
delta / 2.35;
lmfit::result_t res(start_par);
lmfit::result_t res_err(start_par_err);
lmfit::result_t cov(start_cov);
// TODO can we make lmcurve write the result directly where is should be?
lmcurve2(res.par.size(), res.par.data(), res_err.par.data(), cov.par.data(),
x.size(), x.data(), y.data(), y_err.data(), aare::func::gaus,
&control, &res.status);
par_out(0) = res.par[0];
par_out(1) = res.par[1];
par_out(2) = res.par[2];
par_err_out(0) = res_err.par[0];
par_err_out(1) = res_err.par[1];
par_err_out(2) = res_err.par[2];
}
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) {
// Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() ||
par_out.size() != 2 || par_err_out.size() != 2) {
throw std::runtime_error("Data, x, data_err must have the same size "
"and par_out, par_err_out must have size 2");
}
lm_control_struct control = lm_control_double;
// Estimate the initial parameters for the fit
std::vector<double> start_par{0, 0};
std::vector<double> start_par_err{0, 0};
std::vector<double> start_cov{0, 0, 0, 0};
auto y2 = std::max_element(y.begin(), y.end());
auto x2 = x[std::distance(y.begin(), y2)];
auto y1 = std::min_element(y.begin(), y.end());
auto x1 = x[std::distance(y.begin(), y1)];
start_par[0] =
(*y2 - *y1) / (x2 - x1); // For amplitude we use the maximum value
start_par[1] =
*y1 - ((*y2 - *y1) / (x2 - x1)) *
x1; // For the mean we use the x value of the maximum value
lmfit::result_t res(start_par);
lmfit::result_t res_err(start_par_err);
lmfit::result_t cov(start_cov);
lmcurve2(res.par.size(), res.par.data(), res_err.par.data(), cov.par.data(),
x.size(), x.data(), y.data(), y_err.data(), aare::func::pol1,
&control, &res.status);
par_out(0) = res.par[0];
par_out(1) = res.par[1];
par_err_out(0) = res_err.par[0];
par_err_out(1) = res_err.par[1];
}
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,
int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
NDView<double, 1> y_err_view(&y_err(row, col, 0),
{y_err.shape(2)});
NDView<double, 1> par_out_view(&par_out(row, col, 0),
{par_out.shape(2)});
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);
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads;
for (auto &task : tasks) {
threads.push_back(std::thread(process, task.first, task.second));
}
for (auto &thread : threads) {
thread.join();
}
}
NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y) {
// // Check that we have the correct sizes
// if (y.size() != x.size() || y.size() != y_err.size() ||
// par_out.size() != 2 || par_err_out.size() != 2) {
// throw std::runtime_error("Data, x, data_err must have the same size "
// "and par_out, par_err_out must have size 2");
// }
NDArray<double, 1> par({2}, 0);
lm_control_struct control = lm_control_double;
// Estimate the initial parameters for the fit
std::vector<double> start_par{0, 0};
auto y2 = std::max_element(y.begin(), y.end());
auto x2 = x[std::distance(y.begin(), y2)];
auto y1 = std::min_element(y.begin(), y.end());
auto x1 = x[std::distance(y.begin(), y1)];
start_par[0] = (*y2 - *y1) / (x2 - x1);
start_par[1] = *y1 - ((*y2 - *y1) / (x2 - x1)) * x1;
lmfit::result_t res(start_par);
lmcurve(res.par.size(), res.par.data(), x.size(), x.data(), y.data(),
aare::func::pol1, &control, &res.status);
par(0) = res.par[0];
par(1) = res.par[1];
return par;
}
NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 2}, 0);
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
auto res = fit_pol1(x, values);
result(row, col, 0) = res(0);
result(row, col, 1) = res(1);
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads;
for (auto &task : tasks) {
threads.push_back(std::thread(process, task.first, task.second));
}
for (auto &thread : threads) {
thread.join();
}
return result;
}
} // namespace aare

159
src/RawFile.cpp
View File

@ -1,6 +1,7 @@
#include "aare/RawFile.hpp"
#include "aare/PixelMap.hpp"
#include "aare/defs.hpp"
#include "aare/geo_helpers.hpp"
#include <fmt/format.h>
#include <nlohmann/json.hpp>
@ -21,8 +22,11 @@ RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode)
find_geometry();
update_geometry_with_roi();
if (m_master.roi()){
m_geometry = update_geometry_with_roi(m_geometry, m_master.roi().value());
}
open_subfiles();
} else {
throw std::runtime_error(LOCATION +
@ -72,9 +76,13 @@ 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_rows * m_cols * m_master.bitdepth() / 8;
return m_geometry.pixels_x * m_geometry.pixels_y * m_master.bitdepth() / 8;
}
size_t RawFile::pixels_per_frame() {
// return m_rows * m_cols;
return m_geometry.pixels_x * m_geometry.pixels_y;
}
size_t RawFile::pixels_per_frame() { return m_rows * m_cols; }
DetectorType RawFile::detector_type() const { return m_master.detector_type(); }
@ -92,8 +100,8 @@ void RawFile::seek(size_t frame_index) {
size_t RawFile::tell() { return m_current_frame; };
size_t RawFile::total_frames() const { return m_master.frames_in_file(); }
size_t RawFile::rows() const { return m_rows; }
size_t RawFile::cols() const { return m_cols; }
size_t RawFile::rows() const { return m_geometry.pixels_y; }
size_t RawFile::cols() const { return m_geometry.pixels_x; }
size_t RawFile::bitdepth() const { return m_master.bitdepth(); }
xy RawFile::geometry() { return m_master.geometry(); }
@ -102,11 +110,11 @@ void RawFile::open_subfiles() {
for (size_t i = 0; i != n_subfiles; ++i) {
auto v = std::vector<RawSubFile *>(n_subfile_parts);
for (size_t j = 0; j != n_subfile_parts; ++j) {
auto pos = m_module_pixel_0[j];
auto pos = m_geometry.module_pixel_0[j];
v[j] = new RawSubFile(m_master.data_fname(j, i),
m_master.detector_type(), pos.height,
pos.width, m_master.bitdepth(),
positions[j].row, positions[j].col);
pos.row_index, pos.col_index);
}
subfiles.push_back(v);
@ -149,112 +157,49 @@ int RawFile::find_number_of_subfiles() {
RawMasterFile RawFile::master() const { return m_master; }
/**
* @brief Find the geometry of the detector by opening all the subfiles and
* reading the headers.
*/
void RawFile::find_geometry() {
//Hold the maximal row and column number found
//Later used for calculating the total number of rows and columns
uint16_t r{};
uint16_t c{};
for (size_t i = 0; i < n_subfile_parts; i++) {
auto h = this->read_header(m_master.data_fname(i, 0));
auto h = read_header(m_master.data_fname(i, 0));
r = std::max(r, h.row);
c = std::max(c, h.column);
positions.push_back({h.row, h.column});
// positions.push_back({h.row, h.column});
ModuleGeometry g;
g.x = h.column * m_master.pixels_x();
g.y = h.row * m_master.pixels_y();
g.origin_x = h.column * m_master.pixels_x();
g.origin_y = h.row * m_master.pixels_y();
g.row_index = h.row;
g.col_index = h.column;
g.width = m_master.pixels_x();
g.height = m_master.pixels_y();
m_module_pixel_0.push_back(g);
m_geometry.module_pixel_0.push_back(g);
}
r++;
c++;
m_rows = (r * m_master.pixels_y());
m_cols = (c * m_master.pixels_x());
m_rows += static_cast<size_t>((r - 1) * cfg.module_gap_row);
#ifdef AARE_VERBOSE
fmt::print("\nRawFile::find_geometry()\n");
for (size_t i = 0; i < m_module_pixel_0.size(); i++) {
fmt::print("Module {} at position: (r:{},c:{})\n", i,
m_module_pixel_0[i].y, m_module_pixel_0[i].x);
}
fmt::print("Image size: {}x{}\n\n", m_rows, m_cols);
#endif
}
void RawFile::update_geometry_with_roi() {
// TODO! implement this
if (m_master.roi()) {
auto roi = m_master.roi().value();
// TODO! can we do this cleaner?
int pos_y = 0;
int pos_y_increment = 0;
for (size_t row = 0; row < m_master.geometry().row; row++) {
int pos_x = 0;
for (size_t col = 0; col < m_master.geometry().col; col++) {
auto &m = m_module_pixel_0[row * m_master.geometry().col + col];
auto original_height = m.height;
auto original_width = m.width;
// module is to the left of the roi
if (m.x + m.width < roi.xmin) {
m.width = 0;
// roi is in module
} else {
// here we only arrive when the roi is in or to the left of
// the module
if (roi.xmin > m.x) {
m.width -= roi.xmin - m.x;
}
if (roi.xmax < m.x + m.width) {
m.width -= m.x + original_width - roi.xmax;
}
m.x = pos_x;
pos_x += m.width;
}
if (m.y + m.height < roi.ymin) {
m.height = 0;
} else {
if ((roi.ymin > m.y) && (roi.ymin < m.y + m.height)) {
m.height -= roi.ymin - m.y;
}
if (roi.ymax < m.y + m.height) {
m.height -= m.y + original_height - roi.ymax;
}
m.y = pos_y;
pos_y_increment = m.height;
}
}
// increment pos_y
pos_y += pos_y_increment;
}
m_rows = roi.height();
m_cols = roi.width();
}
#ifdef AARE_VERBOSE
fmt::print("RawFile::update_geometry_with_roi()\n");
for (const auto &m : m_module_pixel_0) {
fmt::print("Module at position: (r:{}, c:{}, h:{}, w:{})\n", m.y, m.x,
m.height, m.width);
}
fmt::print("Updated image size: {}x{}\n\n", m_rows, m_cols);
fmt::print("\n");
#endif
m_geometry.pixels_y = (r * m_master.pixels_y());
m_geometry.pixels_x = (c * m_master.pixels_x());
m_geometry.modules_x = c;
m_geometry.modules_y = r;
m_geometry.pixels_y += static_cast<size_t>((r - 1) * cfg.module_gap_row);
}
Frame RawFile::get_frame(size_t frame_index) {
auto f = Frame(m_rows, m_cols, 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;
@ -278,6 +223,10 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect
if (n_subfile_parts != 1) {
for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) {
auto subfile_id = frame_index / m_master.max_frames_per_file();
if (subfile_id >= subfiles.size()) {
throw std::runtime_error(LOCATION +
" Subfile out of range. Possible missing data.");
}
frame_numbers[part_idx] =
subfiles[subfile_id][part_idx]->frame_number(
frame_index % m_master.max_frames_per_file());
@ -311,12 +260,16 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect
for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) {
auto corrected_idx = frame_indices[part_idx];
auto subfile_id = corrected_idx / m_master.max_frames_per_file();
if (subfile_id >= subfiles.size()) {
throw std::runtime_error(LOCATION +
" Subfile out of range. Possible missing data.");
}
// This is where we start writing
auto offset = (m_module_pixel_0[part_idx].y * m_cols +
m_module_pixel_0[part_idx].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_module_pixel_0[part_idx].x!=0)
if (m_geometry.module_pixel_0[part_idx].origin_x!=0)
throw std::runtime_error(LOCATION + "Implementation error. x pos not 0.");
//TODO! Risk for out of range access
@ -340,9 +293,13 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect
// level
for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) {
auto pos = m_module_pixel_0[part_idx];
auto pos = m_geometry.module_pixel_0[part_idx];
auto corrected_idx = frame_indices[part_idx];
auto subfile_id = corrected_idx / m_master.max_frames_per_file();
if (subfile_id >= subfiles.size()) {
throw std::runtime_error(LOCATION +
" Subfile out of range. Possible missing data.");
}
subfiles[subfile_id][part_idx]->seek(corrected_idx % m_master.max_frames_per_file());
subfiles[subfile_id][part_idx]->read_into(part_buffer, header);
@ -352,9 +309,9 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect
for (size_t cur_row = 0; cur_row < static_cast<size_t>(pos.height);
cur_row++) {
auto irow = (pos.y + cur_row);
auto icol = pos.x;
auto dest = (irow * this->m_cols + icol);
auto irow = (pos.origin_y + cur_row);
auto icol = pos.origin_x;
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 *
@ -400,4 +357,8 @@ RawFile::~RawFile() {
}
}
} // namespace aare

5
src/RawFile.test.cpp
View File

@ -1,10 +1,13 @@
#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]") {
@ -148,3 +151,5 @@ TEST_CASE("Read file with unordered frames", "[.integration]") {
File f(fpath);
REQUIRE_THROWS((f.read_frame()));
}

50
src/RawSubFile.cpp
View File

@ -9,11 +9,13 @@ 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_fname(fname), 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_fname(fname),
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) {
if (m_detector_type == DetectorType::Moench03_old) {
m_pixel_map = GenerateMoench03PixelMap();
}else if(m_detector_type == DetectorType::Eiger && m_pos_row % 2 == 0){
} else if (m_detector_type == DetectorType::Eiger && m_pos_row % 2 == 0) {
m_pixel_map = GenerateEigerFlipRowsPixelMap();
}
@ -42,7 +44,7 @@ RawSubFile::RawSubFile(const std::filesystem::path &fname,
void RawSubFile::seek(size_t frame_index) {
if (frame_index >= n_frames) {
throw std::runtime_error("Frame number out of range");
throw std::runtime_error(LOCATION + fmt::format("Frame index {} out of range in a file with {} frames", frame_index, n_frames));
}
m_file.seekg((sizeof(DetectorHeader) + bytes_per_frame()) * frame_index);
}
@ -51,37 +53,48 @@ size_t RawSubFile::tell() {
return m_file.tellg() / (sizeof(DetectorHeader) + bytes_per_frame());
}
void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) {
if(header){
m_file.read(reinterpret_cast<char *>(header), sizeof(DetectorHeader));
if (header) {
m_file.read(reinterpret_cast<char *>(header), sizeof(DetectorHeader));
} else {
m_file.seekg(sizeof(DetectorHeader), std::ios::cur);
}
//TODO! expand support for different bitdepths
if(m_pixel_map){
// TODO! expand support for different bitdepths
if (m_pixel_map) {
// read into a temporary buffer and then copy the data to the buffer
// in the correct order
// currently this only supports 16 bit data!
auto part_buffer = new std::byte[bytes_per_frame()];
m_file.read(reinterpret_cast<char *>(part_buffer), bytes_per_frame());
auto *data = reinterpret_cast<uint16_t *>(image_buf);
auto *part_data = reinterpret_cast<uint16_t *>(part_buffer);
for (size_t i = 0; i < pixels_per_frame(); i++) {
data[i] = part_data[(*m_pixel_map)(i)];
// TODO! add 4 bit support
if(m_bitdepth == 8){
read_with_map<uint8_t>(image_buf);
}else if (m_bitdepth == 16) {
read_with_map<uint16_t>(image_buf);
} else if (m_bitdepth == 32) {
read_with_map<uint32_t>(image_buf);
}else{
throw std::runtime_error("Unsupported bitdepth for read with pixel map");
}
delete[] part_buffer;
} else {
// read directly into the buffer
m_file.read(reinterpret_cast<char *>(image_buf), bytes_per_frame());
}
}
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);
auto *part_data = reinterpret_cast<T *>(part_buffer);
for (size_t i = 0; i < pixels_per_frame(); i++) {
data[i] = part_data[(*m_pixel_map)(i)];
}
delete[] part_buffer;
}
size_t RawSubFile::rows() const { return m_rows; }
size_t RawSubFile::cols() const { return m_cols; }
void RawSubFile::get_part(std::byte *buffer, size_t frame_index) {
seek(frame_index);
read_into(buffer, nullptr);
@ -94,5 +107,4 @@ size_t RawSubFile::frame_number(size_t frame_index) {
return h.frameNumber;
}
} // namespace aare

61
src/decode.cpp Normal file
View File

@ -0,0 +1,61 @@
#include "aare/decode.hpp"
namespace aare {
uint16_t adc_sar_05_decode64to16(uint64_t input){
//we want bits 29,19,28,18,31,21,27,20,24,23,25,22 and then pad to 16
uint16_t output = 0;
output |= ((input >> 22) & 1) << 11;
output |= ((input >> 25) & 1) << 10;
output |= ((input >> 23) & 1) << 9;
output |= ((input >> 24) & 1) << 8;
output |= ((input >> 20) & 1) << 7;
output |= ((input >> 27) & 1) << 6;
output |= ((input >> 21) & 1) << 5;
output |= ((input >> 31) & 1) << 4;
output |= ((input >> 18) & 1) << 3;
output |= ((input >> 28) & 1) << 2;
output |= ((input >> 19) & 1) << 1;
output |= ((input >> 29) & 1) << 0;
return output;
}
void adc_sar_05_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output){
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));
}
}
}
uint16_t adc_sar_04_decode64to16(uint64_t input){
// bit_map = array([15,17,19,21,23,4,6,8,10,12,14,16] LSB->MSB
uint16_t output = 0;
output |= ((input >> 16) & 1) << 11;
output |= ((input >> 14) & 1) << 10;
output |= ((input >> 12) & 1) << 9;
output |= ((input >> 10) & 1) << 8;
output |= ((input >> 8) & 1) << 7;
output |= ((input >> 6) & 1) << 6;
output |= ((input >> 4) & 1) << 5;
output |= ((input >> 23) & 1) << 4;
output |= ((input >> 21) & 1) << 3;
output |= ((input >> 19) & 1) << 2;
output |= ((input >> 17) & 1) << 1;
output |= ((input >> 15) & 1) << 0;
return output;
}
void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output){
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));
}
}
}
} // namespace aare

71
src/geo_helpers.cpp Normal file
View File

@ -0,0 +1,71 @@
#include "aare/geo_helpers.hpp"
#include "fmt/core.h"
namespace aare{
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);
#endif
int pos_y = 0;
int pos_y_increment = 0;
for (int row = 0; row < geo.modules_y; row++) {
int pos_x = 0;
for (int col = 0; col < geo.modules_x; col++) {
auto &m = geo.module_pixel_0[row * geo.modules_x + col];
auto original_height = m.height;
auto original_width = m.width;
// module is to the left of the roi
if (m.origin_x + m.width < roi.xmin) {
m.width = 0;
// roi is in module
} else {
// here we only arrive when the roi is in or to the left of
// the module
if (roi.xmin > m.origin_x) {
m.width -= roi.xmin - m.origin_x;
}
if (roi.xmax < m.origin_x + original_width) {
m.width -= m.origin_x + original_width - roi.xmax;
}
m.origin_x = pos_x;
pos_x += m.width;
}
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)) {
m.height -= roi.ymin - m.origin_y;
}
if (roi.ymax < m.origin_y + original_height) {
m.height -= m.origin_y + original_height - roi.ymax;
}
m.origin_y = pos_y;
pos_y_increment = m.height;
}
#ifdef AARE_VERBOSE
fmt::println("Module {} {} {} {}", m.origin_x, m.origin_y, m.width, m.height);
#endif
}
// increment pos_y
pos_y += pos_y_increment;
}
// m_rows = roi.height();
// m_cols = roi.width();
geo.pixels_x = roi.width();
geo.pixels_y = roi.height();
return geo;
}
} // namespace aare

230
src/geo_helpers.test.cpp Normal file
View File

@ -0,0 +1,230 @@
#include "aare/File.hpp"
#include "aare/RawMasterFile.hpp" //needed for ROI
#include "aare/RawFile.hpp"
#include <catch2/catch_test_macros.hpp>
#include <filesystem>
#include "aare/geo_helpers.hpp"
#include "test_config.hpp"
TEST_CASE("Simple ROIs on one module"){
// DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi)
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;
geo.modules_y = 1;
geo.module_pixel_0.push_back(mod);
SECTION("ROI is the whole module"){
aare::ROI roi;
roi.xmin = 0;
roi.xmax = 1024;
roi.ymin = 0;
roi.ymax = 512;
auto updated_geo = aare::update_geometry_with_roi(geo, roi);
REQUIRE(updated_geo.pixels_x == 1024);
REQUIRE(updated_geo.pixels_y == 512);
REQUIRE(updated_geo.modules_x == 1);
REQUIRE(updated_geo.modules_y == 1);
REQUIRE(updated_geo.module_pixel_0[0].height == 512);
REQUIRE(updated_geo.module_pixel_0[0].width == 1024);
}
SECTION("ROI is the top left corner of the module"){
aare::ROI roi;
roi.xmin = 100;
roi.xmax = 200;
roi.ymin = 150;
roi.ymax = 200;
auto updated_geo = aare::update_geometry_with_roi(geo, roi);
REQUIRE(updated_geo.pixels_x == 100);
REQUIRE(updated_geo.pixels_y == 50);
REQUIRE(updated_geo.modules_x == 1);
REQUIRE(updated_geo.modules_y == 1);
REQUIRE(updated_geo.module_pixel_0[0].height == 50);
REQUIRE(updated_geo.module_pixel_0[0].width == 100);
}
SECTION("ROI is a small square"){
aare::ROI roi;
roi.xmin = 1000;
roi.xmax = 1010;
roi.ymin = 500;
roi.ymax = 510;
auto updated_geo = aare::update_geometry_with_roi(geo, roi);
REQUIRE(updated_geo.pixels_x == 10);
REQUIRE(updated_geo.pixels_y == 10);
REQUIRE(updated_geo.modules_x == 1);
REQUIRE(updated_geo.modules_y == 1);
REQUIRE(updated_geo.module_pixel_0[0].height == 10);
REQUIRE(updated_geo.module_pixel_0[0].width == 10);
}
SECTION("ROI is a few columns"){
aare::ROI roi;
roi.xmin = 750;
roi.xmax = 800;
roi.ymin = 0;
roi.ymax = 512;
auto updated_geo = aare::update_geometry_with_roi(geo, roi);
REQUIRE(updated_geo.pixels_x == 50);
REQUIRE(updated_geo.pixels_y == 512);
REQUIRE(updated_geo.modules_x == 1);
REQUIRE(updated_geo.modules_y == 1);
REQUIRE(updated_geo.module_pixel_0[0].height == 512);
REQUIRE(updated_geo.module_pixel_0[0].width == 50);
}
}
TEST_CASE("Two modules side by side"){
// 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 = 2048;
geo.pixels_y = 512;
geo.modules_x = 2;
geo.modules_y = 1;
geo.module_pixel_0.push_back(mod);
mod.origin_x = 1024;
geo.module_pixel_0.push_back(mod);
SECTION("ROI is the whole image"){
aare::ROI roi;
roi.xmin = 0;
roi.xmax = 2048;
roi.ymin = 0;
roi.ymax = 512;
auto updated_geo = aare::update_geometry_with_roi(geo, roi);
REQUIRE(updated_geo.pixels_x == 2048);
REQUIRE(updated_geo.pixels_y == 512);
REQUIRE(updated_geo.modules_x == 2);
REQUIRE(updated_geo.modules_y == 1);
}
SECTION("rectangle on both modules"){
aare::ROI roi;
roi.xmin = 800;
roi.xmax = 1300;
roi.ymin = 200;
roi.ymax = 499;
auto updated_geo = aare::update_geometry_with_roi(geo, roi);
REQUIRE(updated_geo.pixels_x == 500);
REQUIRE(updated_geo.pixels_y == 299);
REQUIRE(updated_geo.modules_x == 2);
REQUIRE(updated_geo.modules_y == 1);
REQUIRE(updated_geo.module_pixel_0[0].height == 299);
REQUIRE(updated_geo.module_pixel_0[0].width == 224);
REQUIRE(updated_geo.module_pixel_0[1].height == 299);
REQUIRE(updated_geo.module_pixel_0[1].width == 276);
}
}
TEST_CASE("Three modules side by side"){
// DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi)
aare::DetectorGeometry geo;
aare::ROI roi;
roi.xmin = 700;
roi.xmax = 2500;
roi.ymin = 0;
roi.ymax = 123;
aare::ModuleGeometry mod;
mod.origin_x = 0;
mod.origin_y = 0;
mod.width = 1024;
mod.height = 512;
geo.pixels_x = 3072;
geo.pixels_y = 512;
geo.modules_x = 3;
geo.modules_y = 1;
geo.module_pixel_0.push_back(mod);
mod.origin_x = 1024;
geo.module_pixel_0.push_back(mod);
mod.origin_x = 2048;
geo.module_pixel_0.push_back(mod);
auto updated_geo = aare::update_geometry_with_roi(geo, roi);
REQUIRE(updated_geo.pixels_x == 1800);
REQUIRE(updated_geo.pixels_y == 123);
REQUIRE(updated_geo.modules_x == 3);
REQUIRE(updated_geo.modules_y == 1);
REQUIRE(updated_geo.module_pixel_0[0].height == 123);
REQUIRE(updated_geo.module_pixel_0[0].width == 324);
REQUIRE(updated_geo.module_pixel_0[1].height == 123);
REQUIRE(updated_geo.module_pixel_0[1].width == 1024);
REQUIRE(updated_geo.module_pixel_0[2].height == 123);
REQUIRE(updated_geo.module_pixel_0[2].width == 452);
}
TEST_CASE("Four modules as a square"){
// DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi)
aare::DetectorGeometry geo;
aare::ROI roi;
roi.xmin = 500;
roi.xmax = 2000;
roi.ymin = 500;
roi.ymax = 600;
aare::ModuleGeometry mod;
mod.origin_x = 0;
mod.origin_y = 0;
mod.width = 1024;
mod.height = 512;
geo.pixels_x = 2048;
geo.pixels_y = 1024;
geo.modules_x = 2;
geo.modules_y = 2;
geo.module_pixel_0.push_back(mod);
mod.origin_x = 1024;
geo.module_pixel_0.push_back(mod);
mod.origin_x = 0;
mod.origin_y = 512;
geo.module_pixel_0.push_back(mod);
mod.origin_x = 1024;
geo.module_pixel_0.push_back(mod);
auto updated_geo = aare::update_geometry_with_roi(geo, roi);
REQUIRE(updated_geo.pixels_x == 1500);
REQUIRE(updated_geo.pixels_y == 100);
REQUIRE(updated_geo.modules_x == 2);
REQUIRE(updated_geo.modules_y == 2);
REQUIRE(updated_geo.module_pixel_0[0].height == 12);
REQUIRE(updated_geo.module_pixel_0[0].width == 524);
REQUIRE(updated_geo.module_pixel_0[1].height == 12);
REQUIRE(updated_geo.module_pixel_0[1].width == 976);
REQUIRE(updated_geo.module_pixel_0[2].height == 88);
REQUIRE(updated_geo.module_pixel_0[2].width == 524);
REQUIRE(updated_geo.module_pixel_0[3].height == 88);
REQUIRE(updated_geo.module_pixel_0[3].width == 976);
}

30
src/utils/task.cpp Normal file
View File

@ -0,0 +1,30 @@
#include "aare/utils/task.hpp"
namespace aare {
std::vector<std::pair<int, int>> split_task(int first, int last,
int n_threads) {
std::vector<std::pair<int, int>> vec;
vec.reserve(n_threads);
int n_frames = last - first;
if (n_threads >= n_frames) {
for (int i = 0; i != n_frames; ++i) {
vec.push_back({i, i + 1});
}
return vec;
}
int step = (n_frames) / n_threads;
for (int i = 0; i != n_threads; ++i) {
int start = step * i;
int stop = step * (i + 1);
if (i == n_threads - 1)
stop = last;
vec.push_back({start, stop});
}
return vec;
}
} // namespace aare

32
src/utils/task.test.cpp Normal file
View File

@ -0,0 +1,32 @@
#include "aare/utils/task.hpp"
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/catch_test_macros.hpp>
TEST_CASE("Split a range into multiple tasks"){
auto tasks = aare::split_task(0, 10, 3);
REQUIRE(tasks.size() == 3);
REQUIRE(tasks[0].first == 0);
REQUIRE(tasks[0].second == 3);
REQUIRE(tasks[1].first == 3);
REQUIRE(tasks[1].second == 6);
REQUIRE(tasks[2].first == 6);
REQUIRE(tasks[2].second == 10);
tasks = aare::split_task(0, 10, 1);
REQUIRE(tasks.size() == 1);
REQUIRE(tasks[0].first == 0);
REQUIRE(tasks[0].second == 10);
tasks = aare::split_task(0, 10, 10);
REQUIRE(tasks.size() == 10);
for (int i = 0; i < 10; i++){
REQUIRE(tasks[i].first == i);
REQUIRE(tasks[i].second == i+1);
}
}