detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2025-12-23 21:41:24 +01:00
Code Issues Actions Packages Releases Activity

added g0 calibration, pedestal and pixel counting

Browse Source
This commit is contained in:
froejdh_e
2025-07-22 16:42:09 +02:00
parent b898e1c8d0
commit 9a7713e98a
8 changed files with 379 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@@ -368,6 +368,7 @@ set(PUBLICHEADERS
set(SourceFiles
${CMAKE_CURRENT_SOURCE_DIR}/src/calibration.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/decode.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp

27
include/aare/NDView.hpp
View File

@@ -73,10 +73,21 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
}
template <typename... Ix>
const std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
std::enable_if_t<sizeof...(Ix) == 1 && (Ndim > 1), NDView<T, Ndim - 1>> operator()(Ix... index) {
// return a view of the next dimension
std::array<ssize_t, Ndim - 1> new_shape{};
std::copy_n(shape_.begin() + 1, Ndim - 1, new_shape.begin());
return NDView<T, Ndim - 1>(&buffer_[element_offset(strides_, index...)],
new_shape);
}
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, const T &> operator()(Ix... index) const {
return buffer_[element_offset(strides_, index...)];
}
ssize_t size() const { return static_cast<ssize_t>(size_); }
size_t total_bytes() const { return size_ * sizeof(T); }
std::array<ssize_t, Ndim> strides() const noexcept { return strides_; }
@@ -85,9 +96,19 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
T *end() { return buffer_ + size_; }
T const *begin() const { return buffer_; }
T const *end() const { return buffer_ + size_; }
T &operator()(ssize_t i) { return buffer_[i]; }
/**
* @brief Access element at index i.
*/
T &operator[](ssize_t i) { return buffer_[i]; }
const T &operator()(ssize_t i) const { return buffer_[i]; }
/**
* @brief Access element at index i.
*/
const T &operator[](ssize_t i) const { return buffer_[i]; }
bool operator==(const NDView &other) const {

6
include/aare/VarClusterFinder.hpp
View File

@@ -240,14 +240,14 @@ template <typename T> void VarClusterFinder<T>::first_pass() {
for (ssize_t i = 0; i < original_.size(); ++i) {
if (use_noise_map)
threshold_ = 5 * noiseMap(i);
binary_(i) = (original_(i) > threshold_);
threshold_ = 5 * noiseMap[i];
binary_[i] = (original_[i] > threshold_);
}
for (int i = 0; i < shape_[0]; ++i) {
for (int j = 0; j < shape_[1]; ++j) {
// do we have someting to process?
// do we have something to process?
if (binary_(i, j)) {
auto tmp = check_neighbours(i, j);
if (tmp != 0) {

207
include/aare/calibration.hpp
View File

@@ -1,5 +1,7 @@
#pragma once
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
#include "aare/defs.hpp"
#include "aare/utils/task.hpp"
#include <cstdint>
@@ -61,26 +63,221 @@ void apply_calibration_impl(NDView<T, 3> res, NDView<uint16_t, 3> raw_data,
for (int frame_nr = start; frame_nr != stop; ++frame_nr) {
for (int row = 0; row != raw_data.shape(1); ++row) {
for (int col = 0; col != raw_data.shape(2); ++col) {
auto [value, gain] = get_value_and_gain(raw_data(frame_nr, row, col));
auto [value, gain] =
get_value_and_gain(raw_data(frame_nr, row, col));
// Using multiplication does not seem to speed up the code here
// ADU/keV is the standard unit for the calibration which
// means rewriting the formula is not worth it.
res(frame_nr, row, col) =
(value - ped(gain, row, col)) / cal(gain, row, col); //TODO! use multiplication
(value - ped(gain, row, col)) / cal(gain, row, col);
}
}
}
}
template <class T>
void apply_calibration_impl(NDView<T, 3> res, NDView<uint16_t, 3> raw_data,
NDView<T, 2> ped, NDView<T, 2> cal, int start,
int stop) {
for (int frame_nr = start; frame_nr != stop; ++frame_nr) {
for (int row = 0; row != raw_data.shape(1); ++row) {
for (int col = 0; col != raw_data.shape(2); ++col) {
auto [value, gain] =
get_value_and_gain(raw_data(frame_nr, row, col));
// Using multiplication does not seem to speed up the code here
// ADU/keV is the standard unit for the calibration which
// means rewriting the formula is not worth it.
// ignore anything that is not gain 0
// TODO! deal with fixed gain?
if (gain == 0)
res(frame_nr, row, col) =
(value - ped(row, col)) / cal(row, col);
}
}
}
}
template <class T, ssize_t Ndim=3>
void apply_calibration(NDView<T, 3> res, NDView<uint16_t, 3> raw_data,
NDView<T, 3> ped, NDView<T, 3> cal,
NDView<T, Ndim> ped, NDView<T, Ndim> cal,
ssize_t n_threads = 4) {
std::vector<std::future<void>> futures;
futures.reserve(n_threads);
auto limits = split_task(0, raw_data.shape(0), n_threads);
for (const auto &lim : limits)
futures.push_back(std::async(&apply_calibration_impl<T>, res, raw_data, ped, cal,
lim.first, lim.second));
futures.push_back(std::async(
static_cast<void (*)(NDView<T, 3>, NDView<uint16_t, 3>,
NDView<T, Ndim>, NDView<T, Ndim>, int, int)>(
apply_calibration_impl),
res, raw_data, ped, cal, lim.first, lim.second));
for (auto &f : futures)
f.get();
}
std::pair<NDArray<size_t, 3>, NDArray<size_t, 3>>
sum_and_count_per_gain(NDView<uint16_t, 3> raw_data);
std::pair<NDArray<size_t, 2>, NDArray<size_t, 2>>
sum_and_count_g0(NDView<uint16_t, 3> raw_data);
template <typename T>
NDArray<T, 3> calculate_pedestal(NDView<uint16_t, 3> raw_data) {
auto [accumulator, count] = sum_and_count_per_gain(raw_data);
NDArray<T, 3> pedestal(
std::array<ssize_t, 3>{3, raw_data.shape(1), raw_data.shape(2)}, 0);
for (int gain = 0; gain < 3; ++gain) {
for (int row = 0; row < raw_data.shape(1); ++row) {
for (int col = 0; col < raw_data.shape(2); ++col) {
if (count(gain, row, col) != 0) {
pedestal(gain, row, col) =
static_cast<T>(accumulator(gain, row, col)) /
static_cast<T>(count(gain, row, col));
}
}
}
}
return pedestal;
}
template <typename T>
NDArray<T, 2> calculate_pedestal_g0(NDView<uint16_t, 3> raw_data) {
auto [accumulator, count] = sum_and_count_g0(raw_data);
NDArray<T, 2> pedestal(
std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
for (int row = 0; row < raw_data.shape(1); ++row) {
for (int col = 0; col < raw_data.shape(2); ++col) {
if (count(row, col) != 0) {
pedestal(row, col) =
static_cast<T>(accumulator(row, col)) /
static_cast<T>(count(row, col));
}
}
}
return pedestal;
}
template <typename T>
NDArray<T, 2> calculate_pedestal_g0(NDView<uint16_t, 3> raw_data,
ssize_t n_threads) {
std::vector<std::future<std::pair<NDArray<size_t, 2>, NDArray<size_t, 2>>>>
futures;
futures.reserve(n_threads);
auto limits = split_task(0, raw_data.shape(0), n_threads);
// make subviews for each thread
std::vector<NDView<uint16_t, 3>> subviews;
for (const auto &lim : limits) {
subviews.emplace_back(
raw_data.data() + lim.first * raw_data.strides()[0],
std::array<ssize_t, 3>{lim.second - lim.first, raw_data.shape(1),
raw_data.shape(2)});
}
for (auto view : subviews) {
futures.push_back(std::async(
static_cast<std::pair<NDArray<size_t, 2>, NDArray<size_t, 2>> (*)(
NDView<uint16_t, 3>)>(&sum_and_count_g0),
view));
}
NDArray<size_t, 2> accumulator(
std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
NDArray<size_t, 2> count(
std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
for (auto &f : futures) {
auto [acc, cnt] = f.get();
accumulator += acc;
count += cnt;
}
NDArray<T, 2> pedestal(
std::array<ssize_t,2>{raw_data.shape(1), raw_data.shape(2)}, 0);
for (int gain = 0; gain < 3; ++gain) {
for (int row = 0; row < raw_data.shape(1); ++row) {
for (int col = 0; col < raw_data.shape(2); ++col) {
if (count(row, col) != 0) {
pedestal(row, col) =
static_cast<T>(accumulator(row, col)) /
static_cast<T>(count(row, col));
}
}
}
}
return pedestal;
}
template <typename T>
NDArray<T, 3> calculate_pedestal(NDView<uint16_t, 3> raw_data,
ssize_t n_threads) {
NDArray<int, 2> switched(
std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
std::vector<std::future<std::pair<NDArray<size_t, 3>, NDArray<size_t, 3>>>>
futures;
futures.reserve(n_threads);
auto limits = split_task(0, raw_data.shape(0), n_threads);
// make subviews for each thread
std::vector<NDView<uint16_t, 3>> subviews;
for (const auto &lim : limits) {
subviews.emplace_back(
raw_data.data() + lim.first * raw_data.strides()[0],
std::array<ssize_t, 3>{lim.second - lim.first, raw_data.shape(1),
raw_data.shape(2)});
}
for (auto view : subviews) {
futures.push_back(std::async(
static_cast<std::pair<NDArray<size_t, 3>, NDArray<size_t, 3>> (*)(
NDView<uint16_t, 3>)>(&sum_and_count_per_gain),
view));
}
NDArray<size_t, 3> accumulator(
std::array<ssize_t, 3>{3, raw_data.shape(1), raw_data.shape(2)}, 0);
NDArray<size_t, 3> count(
std::array<ssize_t, 3>{3, raw_data.shape(1), raw_data.shape(2)}, 0);
for (auto &f : futures) {
auto [acc, cnt] = f.get();
accumulator += acc;
count += cnt;
}
NDArray<T, 3> pedestal(
std::array<ssize_t, 3>{3, raw_data.shape(1), raw_data.shape(2)}, 0);
for (int gain = 0; gain < 3; ++gain) {
for (int row = 0; row < raw_data.shape(1); ++row) {
for (int col = 0; col < raw_data.shape(2); ++col) {
if (count(gain, row, col) != 0) {
pedestal(gain, row, col) =
static_cast<T>(accumulator(gain, row, col)) /
static_cast<T>(count(gain, row, col));
}
}
}
}
return pedestal;
}
/**
* @brief Count the number of switching pixels in the raw data.
* This function counts the number of pixels that switch between G1 and G2 gain.
* It returns an NDArray with the number of switching pixels per pixel.
* @param raw_data The NDView containing the raw data
* @return An NDArray with the number of switching pixels per pixel
*/
NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data);
/**
* @brief Count the number of switching pixels in the raw data.
* This function counts the number of pixels that switch between G1 and G2 gain.
* It returns an NDArray with the number of switching pixels per pixel.
* @param raw_data The NDView containing the raw data
* @param n_threads The number of threads to use for parallel processing
* @return An NDArray with the number of switching pixels per pixel
*/
NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data,
ssize_t n_threads);
} // namespace aare

3
python/aare/__init__.py
View File

@@ -32,6 +32,7 @@ from .utils import random_pixels, random_pixel, flat_list, add_colorbar
from .func import *
from .calibration import *
from ._aare import apply_calibration
from ._aare import apply_calibration, count_switching_pixels
from ._aare import calculate_pedestal, calculate_pedestal_float, calculate_pedestal_g0, calculate_pedestal_g0_float
from ._aare import VarClusterFinder

71
python/src/bind_calibration.hpp
View File

@@ -17,19 +17,60 @@ py::array_t<DataType> pybind_apply_calibration(
calibration,
int n_threads = 4) {
auto data_span = make_view_3d(data);
auto ped = make_view_3d(pedestal);
auto cal = make_view_3d(calibration);
auto data_span = make_view_3d(data); // data is always 3D
/* No pointer is passed, so NumPy will allocate the buffer */
auto result = py::array_t<DataType>(data_span.shape());
auto res = make_view_3d(result);
aare::apply_calibration<DataType>(res, data_span, ped, cal, n_threads);
if (data.ndim() == 3 && pedestal.ndim() == 3 && calibration.ndim() == 3) {
auto ped = make_view_3d(pedestal);
auto cal = make_view_3d(calibration);
aare::apply_calibration<DataType, 3>(res, data_span, ped, cal,
n_threads);
} else if (data.ndim() == 3 && pedestal.ndim() == 2 &&
calibration.ndim() == 2) {
auto ped = make_view_2d(pedestal);
auto cal = make_view_2d(calibration);
aare::apply_calibration<DataType, 2>(res, data_span, ped, cal,
n_threads);
} else {
throw std::runtime_error(
"Invalid number of dimensions for data, pedestal or calibration");
}
return result;
}
py::array_t<int> pybind_count_switching_pixels(
py::array_t<uint16_t, py::array::c_style | py::array::forcecast> data,
ssize_t n_threads = 4) {
auto data_span = make_view_3d(data);
auto arr = new NDArray<int, 2>{};
*arr = aare::count_switching_pixels(data_span, n_threads);
return return_image_data(arr);
}
template <typename T>
py::array_t<T> pybind_calculate_pedestal(
py::array_t<uint16_t, py::array::c_style | py::array::forcecast> data,
ssize_t n_threads) {
auto data_span = make_view_3d(data);
auto arr = new NDArray<T, 3>{};
*arr = aare::calculate_pedestal<T>(data_span, n_threads);
return return_image_data(arr);
}
template <typename T>
py::array_t<T> pybind_calculate_pedestal_g0(
py::array_t<uint16_t, py::array::c_style | py::array::forcecast> data,
ssize_t n_threads) {
auto data_span = make_view_3d(data);
auto arr = new NDArray<T, 2>{};
*arr = aare::calculate_pedestal_g0<T>(data_span, n_threads);
return return_image_data(arr);
}
void bind_calibration(py::module &m) {
m.def("apply_calibration", &pybind_apply_calibration<float>,
py::arg("raw_data").noconvert(), py::kw_only(),
@@ -40,4 +81,20 @@ void bind_calibration(py::module &m) {
py::arg("raw_data").noconvert(), py::kw_only(),
py::arg("pd").noconvert(), py::arg("cal").noconvert(),
py::arg("n_threads") = 4);
m.def("count_switching_pixels", &pybind_count_switching_pixels,
py::arg("raw_data").noconvert(), py::kw_only(),
py::arg("n_threads") = 4);
m.def("calculate_pedestal_float", &pybind_calculate_pedestal<float>,
py::arg("raw_data").noconvert(), py::arg("n_threads") = 4);
m.def("calculate_pedestal", &pybind_calculate_pedestal<double>,
py::arg("raw_data").noconvert(), py::arg("n_threads") = 4);
m.def("calculate_pedestal_g0", &pybind_calculate_pedestal_g0<double>,
py::arg("raw_data").noconvert(), py::arg("n_threads") = 4);
m.def("calculate_pedestal_g0_float", &pybind_calculate_pedestal_g0<float>,
py::arg("raw_data").noconvert(), py::arg("n_threads") = 4);
}

4
src/NDArray.test.cpp
View File

@@ -25,13 +25,13 @@ TEST_CASE("Construct from an NDView") {
REQUIRE(image.data() != view.data());
for (uint32_t i = 0; i < image.size(); ++i) {
REQUIRE(image(i) == view(i));
REQUIRE(image[i] == view[i]);
}
// Changing the image doesn't change the view
image = 43;
for (uint32_t i = 0; i < image.size(); ++i) {
REQUIRE(image(i) != view(i));
REQUIRE(image[i] != view[i]);
}
}

78
src/calibration.cpp Normal file
View File

@@ -0,0 +1,78 @@
#include "aare/calibration.hpp"
namespace aare {
std::pair<NDArray<size_t, 3>, NDArray<size_t,3>> sum_and_count_per_gain(NDView<uint16_t,3> raw_data){
NDArray<size_t, 3> accumulator(std::array<ssize_t, 3>{3, raw_data.shape(1), raw_data.shape(2)}, 0);
NDArray<size_t, 3> count(std::array<ssize_t, 3>{3, raw_data.shape(1), raw_data.shape(2)}, 0);
for (int frame_nr = 0; frame_nr != raw_data.shape(0); ++frame_nr) {
for (int row = 0; row != raw_data.shape(1); ++row) {
for (int col = 0; col != raw_data.shape(2); ++col) {
auto [value, gain] = get_value_and_gain(raw_data(frame_nr, row, col));
accumulator(gain, row, col) += value;
count(gain, row, col) += 1;
}
}
}
return {std::move(accumulator), std::move(count)};
}
std::pair<NDArray<size_t, 2>, NDArray<size_t,2>> sum_and_count_g0(NDView<uint16_t, 3> raw_data){
NDArray<size_t, 2> accumulator(std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
NDArray<size_t, 2> count(std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
for (int frame_nr = 0; frame_nr != raw_data.shape(0); ++frame_nr) {
for (int row = 0; row != raw_data.shape(1); ++row) {
for (int col = 0; col != raw_data.shape(2); ++col) {
auto [value, gain] = get_value_and_gain(raw_data(frame_nr, row, col));
if (gain != 0)
continue; // we only care about gain 0
accumulator(row, col) += value;
count(row, col) += 1;
}
}
}
return {std::move(accumulator), std::move(count)};
}
NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data) {
NDArray<int, 2> switched(std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)},0);
for (int frame_nr = 0; frame_nr != raw_data.shape(0); ++frame_nr) {
for (int row = 0; row != raw_data.shape(1); ++row) {
for (int col = 0; col != raw_data.shape(2); ++col) {
auto [value, gain] = get_value_and_gain(raw_data(frame_nr, row, col));
if (gain != 0) {
switched(row, col) += 1;
}
}
}
}
return switched;
}
NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data, ssize_t n_threads){
NDArray<int, 2> switched(std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)},0);
std::vector<std::future<NDArray<int, 2>>> futures;
futures.reserve(n_threads);
auto limits = split_task(0, raw_data.shape(0), n_threads);
// make subviews for each thread
std::vector<NDView<uint16_t, 3>> subviews;
for (const auto &lim : limits) {
subviews.emplace_back(raw_data.data() + lim.first * raw_data.strides()[0],
std::array<ssize_t, 3>{lim.second-lim.first, raw_data.shape(1), raw_data.shape(2)});
}
for (auto view : subviews) {
futures.push_back(std::async(static_cast<NDArray<int, 2>(*)(NDView<uint16_t, 3>)>(&count_switching_pixels), view));
}
for (auto &f : futures) {
switched += f.get();
}
return switched;
}
}// namespace aare