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

WIP

Browse Source
This commit is contained in:
froejdh_e 2024-12-11 16:27:36 +01:00
parent 60534add92
commit b3a9e9576b
8 changed files with 301 additions and 183 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -344,6 +344,7 @@ if(AARE_TESTS)
${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterVector.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp

1
docs/conf.py.in
View File

@ -29,7 +29,6 @@ version = '@PROJECT_VERSION@'
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
# ones. # ones.
extensions = ['breathe', extensions = ['breathe',
'sphinx_rtd_theme',
'sphinx.ext.autodoc', 'sphinx.ext.autodoc',
'sphinx.ext.napoleon', 'sphinx.ext.napoleon',
] ]

6
docs/src/ClusterVector.rst Normal file
View File

@ -0,0 +1,6 @@
ClusterVector
=============
.. doxygenclass:: aare::ClusterVector
:members:
:undoc-members:

1
docs/src/index.rst
View File

@ -46,6 +46,7 @@ AARE
Dtype Dtype
ClusterFinder ClusterFinder
ClusterFile ClusterFile
ClusterVector
Pedestal Pedestal
RawFile RawFile
RawSubFile RawSubFile

305
include/aare/ClusterFinder.hpp
View File

@ -23,62 +23,83 @@ enum class eventType {
UNDEFINED_EVENT = -1 /** undefined */ UNDEFINED_EVENT = -1 /** undefined */
}; };
template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double, typename CT = int32_t> template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
typename CT = int32_t>
class ClusterFinder { class ClusterFinder {
Shape<2> m_image_size; Shape<2> m_image_size;
const int m_cluster_sizeX; const int m_cluster_sizeX;
const int m_cluster_sizeY; const int m_cluster_sizeY;
const double m_threshold; // const PEDESTAL_TYPE m_threshold;
const double m_nSigma; const PEDESTAL_TYPE m_nSigma;
const double c2; const PEDESTAL_TYPE c2;
const double c3; const PEDESTAL_TYPE c3;
Pedestal<PEDESTAL_TYPE> m_pedestal; Pedestal<PEDESTAL_TYPE> m_pedestal;
ClusterVector<CT> m_clusters; ClusterVector<CT> m_clusters;
public: public:
/**
* @brief Construct a new ClusterFinder object
* @param image_size size of the image
* @param cluster_size size of the cluster (x, y)
* @param nSigma number of sigma above the pedestal to consider a photon
* @param capacity initial capacity of the cluster vector
*
*/
ClusterFinder(Shape<2> image_size, Shape<2> cluster_size, ClusterFinder(Shape<2> image_size, Shape<2> cluster_size,
double nSigma = 5.0, double threshold = 0.0) PEDESTAL_TYPE nSigma = 5.0, size_t capacity = 1000000)
: m_image_size(image_size), m_cluster_sizeX(cluster_size[0]), : m_image_size(image_size), m_cluster_sizeX(cluster_size[0]),
m_cluster_sizeY(cluster_size[1]), m_threshold(threshold), m_cluster_sizeY(cluster_size[1]),
m_nSigma(nSigma), m_nSigma(nSigma),
c2(sqrt((m_cluster_sizeY + 1) / 2 * (m_cluster_sizeX + 1) / 2)), c2(sqrt((m_cluster_sizeY + 1) / 2 * (m_cluster_sizeX + 1) / 2)),
c3(sqrt(m_cluster_sizeX * m_cluster_sizeY)), c3(sqrt(m_cluster_sizeX * m_cluster_sizeY)),
m_pedestal(image_size[0], image_size[1]), m_pedestal(image_size[0], image_size[1]),
m_clusters(m_cluster_sizeX, m_cluster_sizeY, 1'000'000) { m_clusters(m_cluster_sizeX, m_cluster_sizeY, capacity) {};
// clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY, 2000);
};
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) { void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
m_pedestal.push(frame); m_pedestal.push(frame);
} }
NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); } NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); }
NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); } NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
ClusterVector<CT> steal_clusters() { /**
* @brief Move the clusters from the ClusterVector in the ClusterFinder to a
* new ClusterVector and return it.
* @param realloc_same_capacity if true the new ClusterVector will have the
* same capacity as the old one
*
*/
ClusterVector<CT> steal_clusters(bool realloc_same_capacity = false) {
ClusterVector<CT> tmp = std::move(m_clusters); ClusterVector<CT> tmp = std::move(m_clusters);
m_clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY, 2000); if (realloc_same_capacity)
m_clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY,
tmp.capacity());
else
m_clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY);
return tmp; return tmp;
} }
void void find_clusters(NDView<FRAME_TYPE, 2> frame) {
find_clusters(NDView<FRAME_TYPE, 2> frame) {
// // size_t capacity = 2000;
// // ClusterVector<CT> clusters(m_cluster_sizeX, m_cluster_sizeY, capacity);
eventType event_type = eventType::PEDESTAL;
// // TODO! deal with even size clusters // // TODO! deal with even size clusters
// // currently 3,3 -> +/- 1 // // currently 3,3 -> +/- 1
// // 4,4 -> +/- 2 // // 4,4 -> +/- 2
short dy = m_cluster_sizeY / 2; short dy = m_cluster_sizeY / 2;
short dx = m_cluster_sizeX / 2; short dx = m_cluster_sizeX / 2;
std::vector<CT> cluster_data(m_cluster_sizeX * m_cluster_sizeY);
for (int iy = 0; iy < frame.shape(0); iy++) { for (int iy = 0; iy < frame.shape(0); iy++) {
for (int ix = 0; ix < frame.shape(1); ix++) { for (int ix = 0; ix < frame.shape(1); ix++) {
PEDESTAL_TYPE max = std::numeric_limits<FRAME_TYPE>::min(); PEDESTAL_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
PEDESTAL_TYPE total = 0; PEDESTAL_TYPE total = 0;
// What can we short circuit here?
PEDESTAL_TYPE rms = m_pedestal.std(iy, ix);
PEDESTAL_TYPE value = (frame(iy, ix) - m_pedestal.mean(iy, ix));
if (value < -m_nSigma * rms)
continue; // NEGATIVE_PEDESTAL go to next pixel
// TODO! No pedestal update???
for (short ir = -dy; ir < dy + 1; ir++) { for (short ir = -dy; ir < dy + 1; ir++) {
for (short ic = -dx; ic < dx + 1; ic++) { for (short ic = -dx; ic < dx + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
@ -92,159 +113,157 @@ class ClusterFinder {
} }
} }
} }
PEDESTAL_TYPE rms = m_pedestal.std(iy, ix);
PEDESTAL_TYPE value = (frame(iy, ix) - m_pedestal.mean(iy, ix));
if (value < -m_nSigma * rms) { if ((max > m_nSigma * rms)) {
continue; // NEGATIVE_PEDESTAL go to next pixel
// TODO! No pedestal update???
} else if (max > m_nSigma * rms) {
event_type = eventType::PHOTON;
if (value < max) if (value < max)
continue; // Not max go to the next pixel continue; // Not max go to the next pixel
// but also no pedestal update
} else if (total > c3 * m_nSigma * rms) { } else if (total > c3 * m_nSigma * rms) {
event_type = eventType::PHOTON; // pass
} else { } else {
m_pedestal.push(iy, ix, frame(iy, ix)); m_pedestal.push(iy, ix, frame(iy, ix));
continue; // It was a pedestal value nothing to store continue; // It was a pedestal value nothing to store
} }
// Store cluster // Store cluster
if (event_type == eventType::PHOTON && value >= max) { if (value == max) {
event_type = eventType::PHOTON_MAX; // Zero out the cluster data
std::fill(cluster_data.begin(), cluster_data.end(), 0);
short i = 0; // Fill the cluster data since we have a photon to store
std::vector<CT> cluster_data(m_cluster_sizeX * // It's worth redoing the look since most of the time we
m_cluster_sizeY); // don't have a photon
int i = 0;
for (short ir = -dy; ir < dy + 1; ir++) { for (int ir = -dy; ir < dy + 1; ir++) {
for (short ic = -dx; ic < dx + 1; ic++) { for (int ic = -dx; ic < dx + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) { iy + ir >= 0 && iy + ir < frame.shape(0)) {
CT tmp = CT tmp =
static_cast<CT>( static_cast<CT>(frame(iy + ir, ix + ic)) -
frame(iy + ir, ix + ic)) -
m_pedestal.mean(iy + ir, ix + ic); m_pedestal.mean(iy + ir, ix + ic);
cluster_data[i] = tmp; //Watch for out of bounds access cluster_data[i] =
tmp; // Watch for out of bounds access
i++; i++;
} }
} }
} }
// Add the cluster to the output ClusterVector
m_clusters.push_back( m_clusters.push_back(
ix, iy, ix, iy,
reinterpret_cast<std::byte *>(cluster_data.data())); reinterpret_cast<std::byte *>(cluster_data.data()));
} }
} }
} }
// return clusters;
} }
// template <typename FRAME_TYPE, typename PEDESTAL_TYPE> // // template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
std::vector<DynamicCluster> // std::vector<DynamicCluster>
find_clusters_with_threshold(NDView<FRAME_TYPE, 2> frame, // find_clusters_with_threshold(NDView<FRAME_TYPE, 2> frame,
Pedestal<PEDESTAL_TYPE> &pedestal) { // Pedestal<PEDESTAL_TYPE> &pedestal) {
assert(m_threshold > 0); // assert(m_threshold > 0);
std::vector<DynamicCluster> clusters; // std::vector<DynamicCluster> clusters;
std::vector<std::vector<eventType>> eventMask; // std::vector<std::vector<eventType>> eventMask;
for (int i = 0; i < frame.shape(0); i++) { // for (int i = 0; i < frame.shape(0); i++) {
eventMask.push_back(std::vector<eventType>(frame.shape(1))); // eventMask.push_back(std::vector<eventType>(frame.shape(1)));
} // }
double tthr, tthr1, tthr2; // double tthr, tthr1, tthr2;
NDArray<FRAME_TYPE, 2> rest({frame.shape(0), frame.shape(1)}); // NDArray<FRAME_TYPE, 2> rest({frame.shape(0), frame.shape(1)});
NDArray<int, 2> nph({frame.shape(0), frame.shape(1)}); // NDArray<int, 2> nph({frame.shape(0), frame.shape(1)});
// convert to n photons // // convert to n photons
// nph = (frame-pedestal.mean()+0.5*m_threshold)/m_threshold; // can be // // nph = (frame-pedestal.mean()+0.5*m_threshold)/m_threshold; // can
// optimized with expression templates? // be
for (int iy = 0; iy < frame.shape(0); iy++) { // // optimized with expression templates?
for (int ix = 0; ix < frame.shape(1); ix++) { // for (int iy = 0; iy < frame.shape(0); iy++) {
auto val = frame(iy, ix) - pedestal.mean(iy, ix); // for (int ix = 0; ix < frame.shape(1); ix++) {
nph(iy, ix) = (val + 0.5 * m_threshold) / m_threshold; // auto val = frame(iy, ix) - pedestal.mean(iy, ix);
nph(iy, ix) = nph(iy, ix) < 0 ? 0 : nph(iy, ix); // nph(iy, ix) = (val + 0.5 * m_threshold) / m_threshold;
rest(iy, ix) = val - nph(iy, ix) * m_threshold; // nph(iy, ix) = nph(iy, ix) < 0 ? 0 : nph(iy, ix);
} // rest(iy, ix) = val - nph(iy, ix) * m_threshold;
} // }
// iterate over frame pixels // }
for (int iy = 0; iy < frame.shape(0); iy++) { // // iterate over frame pixels
for (int ix = 0; ix < frame.shape(1); ix++) { // for (int iy = 0; iy < frame.shape(0); iy++) {
eventMask[iy][ix] = eventType::PEDESTAL; // for (int ix = 0; ix < frame.shape(1); ix++) {
// initialize max and total // eventMask[iy][ix] = eventType::PEDESTAL;
FRAME_TYPE max = std::numeric_limits<FRAME_TYPE>::min(); // // initialize max and total
long double total = 0; // FRAME_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
if (rest(iy, ix) <= 0.25 * m_threshold) { // long double total = 0;
pedestal.push(iy, ix, frame(iy, ix)); // if (rest(iy, ix) <= 0.25 * m_threshold) {
continue; // pedestal.push(iy, ix, frame(iy, ix));
} // continue;
eventMask[iy][ix] = eventType::NEIGHBOUR; // }
// iterate over cluster pixels around the current pixel (ix,iy) // eventMask[iy][ix] = eventType::NEIGHBOUR;
for (short ir = -(m_cluster_sizeY / 2); // // iterate over cluster pixels around the current pixel
ir < (m_cluster_sizeY / 2) + 1; ir++) { // (ix,iy) for (short ir = -(m_cluster_sizeY / 2);
for (short ic = -(m_cluster_sizeX / 2); // ir < (m_cluster_sizeY / 2) + 1; ir++) {
ic < (m_cluster_sizeX / 2) + 1; ic++) { // for (short ic = -(m_cluster_sizeX / 2);
if (ix + ic >= 0 && ix + ic < frame.shape(1) && // ic < (m_cluster_sizeX / 2) + 1; ic++) {
iy + ir >= 0 && iy + ir < frame.shape(0)) { // if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
auto val = frame(iy + ir, ix + ic) - // iy + ir >= 0 && iy + ir < frame.shape(0)) {
pedestal.mean(iy + ir, ix + ic); // auto val = frame(iy + ir, ix + ic) -
total += val; // pedestal.mean(iy + ir, ix + ic);
if (val > max) { // total += val;
max = val; // if (val > max) {
} // max = val;
} // }
} // }
} // }
// }
auto rms = pedestal.std(iy, ix); // auto rms = pedestal.std(iy, ix);
if (m_nSigma == 0) { // if (m_nSigma == 0) {
tthr = m_threshold; // tthr = m_threshold;
tthr1 = m_threshold; // tthr1 = m_threshold;
tthr2 = m_threshold; // tthr2 = m_threshold;
} else { // } else {
tthr = m_nSigma * rms; // tthr = m_nSigma * rms;
tthr1 = m_nSigma * rms * c3; // tthr1 = m_nSigma * rms * c3;
tthr2 = m_nSigma * rms * c2; // tthr2 = m_nSigma * rms * c2;
if (m_threshold > 2 * tthr) // if (m_threshold > 2 * tthr)
tthr = m_threshold - tthr; // tthr = m_threshold - tthr;
if (m_threshold > 2 * tthr1) // if (m_threshold > 2 * tthr1)
tthr1 = tthr - tthr1; // tthr1 = tthr - tthr1;
if (m_threshold > 2 * tthr2) // if (m_threshold > 2 * tthr2)
tthr2 = tthr - tthr2; // tthr2 = tthr - tthr2;
} // }
if (total > tthr1 || max > tthr) { // if (total > tthr1 || max > tthr) {
eventMask[iy][ix] = eventType::PHOTON; // eventMask[iy][ix] = eventType::PHOTON;
nph(iy, ix) += 1; // nph(iy, ix) += 1;
rest(iy, ix) -= m_threshold; // rest(iy, ix) -= m_threshold;
} else { // } else {
pedestal.push(iy, ix, frame(iy, ix)); // pedestal.push(iy, ix, frame(iy, ix));
continue; // continue;
} // }
if (eventMask[iy][ix] == eventType::PHOTON && // if (eventMask[iy][ix] == eventType::PHOTON &&
frame(iy, ix) - pedestal.mean(iy, ix) >= max) { // frame(iy, ix) - pedestal.mean(iy, ix) >= max) {
eventMask[iy][ix] = eventType::PHOTON_MAX; // eventMask[iy][ix] = eventType::PHOTON_MAX;
DynamicCluster cluster(m_cluster_sizeX, m_cluster_sizeY, // DynamicCluster cluster(m_cluster_sizeX, m_cluster_sizeY,
Dtype(typeid(FRAME_TYPE))); // Dtype(typeid(FRAME_TYPE)));
cluster.x = ix; // cluster.x = ix;
cluster.y = iy; // cluster.y = iy;
short i = 0; // short i = 0;
for (short ir = -(m_cluster_sizeY / 2); // for (short ir = -(m_cluster_sizeY / 2);
ir < (m_cluster_sizeY / 2) + 1; ir++) { // ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2); // for (short ic = -(m_cluster_sizeX / 2);
ic < (m_cluster_sizeX / 2) + 1; ic++) { // ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && // if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) { // iy + ir >= 0 && iy + ir < frame.shape(0)) {
auto tmp = frame(iy + ir, ix + ic) - // auto tmp = frame(iy + ir, ix + ic) -
pedestal.mean(iy + ir, ix + ic); // pedestal.mean(iy + ir, ix + ic);
cluster.set<FRAME_TYPE>(i, tmp); // cluster.set<FRAME_TYPE>(i, tmp);
i++; // i++;
} // }
} // }
} // }
clusters.push_back(cluster); // clusters.push_back(cluster);
} // }
} // }
} // }
return clusters; // return clusters;
} // }
}; };
} // namespace aare } // namespace aare

52
include/aare/ClusterVector.hpp
View File

@ -2,9 +2,18 @@
#include <cstddef> #include <cstddef>
#include <cstdint> #include <cstdint>
#include <numeric> #include <numeric>
#include <vector>
#include <fmt/core.h> #include <fmt/core.h>
namespace aare {
/**
* @brief ClusterVector is a container for clusters of various sizes. It uses a
* contiguous memory buffer to store the clusters.
* @note push_back can invalidate pointers to elements in the container
* @tparam T data type of the pixels in the cluster
*/
template <typename T> class ClusterVector { template <typename T> class ClusterVector {
using value_type = T; using value_type = T;
using coord_t = int16_t; using coord_t = int16_t;
@ -24,6 +33,12 @@ template <typename T> class ClusterVector {
constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:" ; constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:" ;
public: public:
/**
* @brief Construct a new ClusterVector object
* @param cluster_size_x size of the cluster in x direction
* @param cluster_size_y size of the cluster in y direction
* @param capacity initial capacity of the buffer in number of clusters
*/
ClusterVector(coord_t cluster_size_x, coord_t cluster_size_y, ClusterVector(coord_t cluster_size_x, coord_t cluster_size_y,
size_t capacity = 1024) size_t capacity = 1024)
: m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y), : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y),
@ -31,8 +46,6 @@ template <typename T> class ClusterVector {
allocate_buffer(capacity); allocate_buffer(capacity);
} }
~ClusterVector() { ~ClusterVector() {
fmt::print("~ClusterVector - size: {}, capacity: {}\n", m_size,
m_capacity);
delete[] m_data; delete[] m_data;
} }
@ -63,13 +76,24 @@ template <typename T> class ClusterVector {
return *this; return *this;
} }
/**
* @brief Reserve space for at least capacity clusters
* @param capacity number of clusters to reserve space for
* @note If capacity is less than the current capacity, the function does nothing.
*/
void reserve(size_t capacity) { void reserve(size_t capacity) {
if (capacity > m_capacity) { if (capacity > m_capacity) {
allocate_buffer(capacity); allocate_buffer(capacity);
} }
} }
// data better hold data of the right size! /**
* @brief Add a cluster to the vector
* @param x x-coordinate of the cluster
* @param y y-coordinate of the cluster
* @param data pointer to the data of the cluster
* @warning The data pointer must point to a buffer of size cluster_size_x * cluster_size_y * sizeof(T)
*/
void push_back(coord_t x, coord_t y, const std::byte *data) { void push_back(coord_t x, coord_t y, const std::byte *data) {
if (m_size == m_capacity) { if (m_size == m_capacity) {
allocate_buffer(m_capacity * 2); allocate_buffer(m_capacity * 2);
@ -85,6 +109,10 @@ template <typename T> class ClusterVector {
m_size++; m_size++;
} }
/**
* @brief Sum the pixels in each cluster
* @return std::vector<T> vector of sums for each cluster
*/
std::vector<T> sum() { std::vector<T> sum() {
std::vector<T> sums(m_size); std::vector<T> sums(m_size);
const size_t stride = element_offset(); const size_t stride = element_offset();
@ -101,12 +129,23 @@ template <typename T> class ClusterVector {
} }
size_t size() const { return m_size; } size_t size() const { return m_size; }
size_t capacity() const { return m_capacity; }
/**
* @brief Return the offset in bytes for a single cluster
*/
size_t element_offset() const { size_t element_offset() const {
return sizeof(m_cluster_size_x) + sizeof(m_cluster_size_y) + return sizeof(m_cluster_size_x) + sizeof(m_cluster_size_y) +
m_cluster_size_x * m_cluster_size_y * sizeof(T); m_cluster_size_x * m_cluster_size_y * sizeof(T);
} }
/**
* @brief Return the offset in bytes for the i-th cluster
*/
size_t element_offset(size_t i) const { return element_offset() * i; } size_t element_offset(size_t i) const { return element_offset() * i; }
/**
* @brief Return a pointer to the i-th cluster
*/
std::byte *element_ptr(size_t i) { return m_data + element_offset(i); } std::byte *element_ptr(size_t i) { return m_data + element_offset(i); }
int16_t cluster_size_x() const { return m_cluster_size_x; } int16_t cluster_size_x() const { return m_cluster_size_x; }
@ -123,13 +162,12 @@ template <typename T> class ClusterVector {
private: private:
void allocate_buffer(size_t new_capacity) { void allocate_buffer(size_t new_capacity) {
size_t num_bytes = element_offset() * new_capacity; size_t num_bytes = element_offset() * new_capacity;
fmt::print(
"ClusterVector allocating {} elements for a total of {} bytes\n",
new_capacity, num_bytes);
std::byte *new_data = new std::byte[num_bytes]{}; std::byte *new_data = new std::byte[num_bytes]{};
std::copy(m_data, m_data + element_offset() * m_size, new_data); std::copy(m_data, m_data + element_offset() * m_size, new_data);
delete[] m_data; delete[] m_data;
m_data = new_data; m_data = new_data;
m_capacity = new_capacity; m_capacity = new_capacity;
} }
}; };
} // namespace aare

41
python/src/cluster.hpp
View File

@ -10,7 +10,7 @@
#include <pybind11/stl.h> #include <pybind11/stl.h>
namespace py = pybind11; namespace py = pybind11;
using pd_type = double; using pd_type = float;
template <typename T> template <typename T>
void define_cluster_vector(py::module &m, const std::string &typestr) { void define_cluster_vector(py::module &m, const std::string &typestr) {
@ -46,7 +46,9 @@ void define_cluster_vector(py::module &m, const std::string &typestr) {
void define_cluster_finder_bindings(py::module &m) { void define_cluster_finder_bindings(py::module &m) {
py::class_<ClusterFinder<uint16_t, pd_type>>(m, "ClusterFinder") py::class_<ClusterFinder<uint16_t, pd_type>>(m, "ClusterFinder")
.def(py::init<Shape<2>, Shape<2>>()) .def(py::init<Shape<2>, Shape<2>, pd_type, size_t>(), py::arg("image_size"),
py::arg("cluster_size"), py::arg("n_sigma") = 5.0,
py::arg("capacity") = 1'000'000)
.def("push_pedestal_frame", .def("push_pedestal_frame",
[](ClusterFinder<uint16_t, pd_type> &self, [](ClusterFinder<uint16_t, pd_type> &self,
py::array_t<uint16_t> frame) { py::array_t<uint16_t> frame) {
@ -66,10 +68,10 @@ void define_cluster_finder_bindings(py::module &m) {
return return_image_data(arr); return return_image_data(arr);
}) })
.def("steal_clusters", .def("steal_clusters",
[](ClusterFinder<uint16_t, pd_type> &self) { [](ClusterFinder<uint16_t, pd_type> &self, bool realloc_same_capacity) {
auto v = new ClusterVector<int>(self.steal_clusters()); auto v = new ClusterVector<int>(self.steal_clusters(realloc_same_capacity));
return v; return v;
}) }, py::arg("realloc_same_capacity") = false)
.def("find_clusters", .def("find_clusters",
[](ClusterFinder<uint16_t, pd_type> &self, [](ClusterFinder<uint16_t, pd_type> &self,
py::array_t<uint16_t> frame) { py::array_t<uint16_t> frame) {
@ -78,36 +80,11 @@ void define_cluster_finder_bindings(py::module &m) {
return; return;
}); });
m.def("hello", []() {
fmt::print("Hello from C++\n");
auto v = new ClusterVector<int>(3, 3);
int data[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
v->push_back(5, 30, reinterpret_cast<std::byte *>(data));
v->push_back(5, 55, reinterpret_cast<std::byte *>(data));
v->push_back(5, 20, reinterpret_cast<std::byte *>(data));
v->push_back(5, 30, reinterpret_cast<std::byte *>(data));
return v;
});
define_cluster_vector<int>(m, "i"); define_cluster_vector<int>(m, "i");
define_cluster_vector<double>(m, "d"); define_cluster_vector<double>(m, "d");
define_cluster_vector<float>(m, "f");
// py::class_<ClusterVector<int>>(m, "ClusterVector", py::buffer_protocol())
// .def(py::init<int, int>())
// .def("size", &ClusterVector<int>::size)
// .def("element_offset",
// py::overload_cast<>(&ClusterVector<int>::element_offset, py::const_))
// .def_buffer([](ClusterVector<int> &v) -> py::buffer_info {
// return py::buffer_info(
// v.data(), /* Pointer to buffer */
// v.element_offset(), /* Size of one scalar */
// fmt::format("h:x:\nh:y:\n({},{})i:data:", v.cluster_size_x(),
// v.cluster_size_y()), /* Format descriptor */ 1, /* Number of
// dimensions */ {v.size()}, /* Buffer dimensions */
// {v.element_offset()} /* Strides (in bytes) for each index */
// );
// });
py::class_<DynamicCluster>(m, "DynamicCluster", py::buffer_protocol()) py::class_<DynamicCluster>(m, "DynamicCluster", py::buffer_protocol())
.def(py::init<int, int, Dtype>()) .def(py::init<int, int, Dtype>())

77
src/ClusterVector.test.cpp Normal file
View File

@ -0,0 +1,77 @@
#include <cstdint>
#include "aare/ClusterVector.hpp"
// #include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/catch_test_macros.hpp>
using aare::ClusterVector;
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);
REQUIRE(cv.size() == 0);
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);
//Create a cluster and push back into the vector
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() == 4);
//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));
//Check that the data is the same
REQUIRE(c1.x == c2.x);
REQUIRE(c1.y == c2.y);
for(size_t i = 0; i < 4; i++) {
REQUIRE(c1.data[i] == c2.data[i]);
}
}
TEST_CASE("Summing 3x1 clusters of int64"){
struct Cluster_l3x1{
int16_t x;
int16_t y;
int64_t data[3];
};
ClusterVector<int64_t> cv(3, 1, 2);
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 0);
REQUIRE(cv.cluster_size_x() == 3);
REQUIRE(cv.cluster_size_y() == 1);
//Create a cluster and push back into the vector
Cluster_l3x1 c1 = {1, 2, {3, 4, 5}};
cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 1);
Cluster_l3x1 c2 = {6, 7, {8, 9, 10}};
cv.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 2);
Cluster_l3x1 c3 = {11, 12, {13, 14, 15}};
cv.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
REQUIRE(cv.capacity() == 4);
REQUIRE(cv.size() == 3);
auto sums = cv.sum();
REQUIRE(sums.size() == 3);
REQUIRE(sums[0] == 12);
REQUIRE(sums[1] == 27);
REQUIRE(sums[2] == 42);
}