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

ClusterFinder

Browse Source
This commit is contained in:
Erik Fröjdh
2024-11-06 12:41:41 +01:00
parent 5b2809d6b0
commit cbfd1f0b6c
10 changed files with 358 additions and 156 deletions
Download Patch File Download Diff File Expand all files Collapse all files

147
include/aare/ClusterFinder.hpp
View File

@ -3,33 +3,60 @@
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "aare/defs.hpp"
#include <cstddef>
namespace aare {
/** enum to define the event types */
enum eventType {
PEDESTAL, /** pedestal */
NEIGHBOUR, /** neighbour i.e. below threshold, but in the cluster of a photon */
PHOTON, /** photon i.e. above threshold */
PHOTON_MAX, /** maximum of a cluster satisfying the photon conditions */
NEGATIVE_PEDESTAL, /** negative value, will not be accounted for as pedestal in order to
avoid drift of the pedestal towards negative values */
PEDESTAL, /** pedestal */
NEIGHBOUR, /** neighbour i.e. below threshold, but in the cluster of a
photon */
PHOTON, /** photon i.e. above threshold */
PHOTON_MAX, /** maximum of a cluster satisfying the photon conditions */
NEGATIVE_PEDESTAL, /** negative value, will not be accounted for as pedestal
in order to avoid drift of the pedestal towards
negative values */
UNDEFINED_EVENT = -1 /** undefined */
};
template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
class ClusterFinder {
Shape<2> m_image_size;
const int m_cluster_sizeX;
const int m_cluster_sizeY;
const double m_threshold;
const double m_nSigma;
const double c2;
const double c3;
Pedestal<PEDESTAL_TYPE> m_pedestal;
public:
ClusterFinder(int cluster_sizeX, int cluster_sizeY, double nSigma = 5.0, double threshold = 0.0)
: m_cluster_sizeX(cluster_sizeX), m_cluster_sizeY(cluster_sizeY), m_threshold(threshold), m_nSigma(nSigma) {
ClusterFinder(Shape<2> image_size, Shape<2>cluster_size, double nSigma = 5.0,
double threshold = 0.0)
: m_image_size(image_size), m_cluster_sizeX(cluster_size[0]), m_cluster_sizeY(cluster_size[1]),
m_threshold(threshold), m_nSigma(nSigma),
c2(sqrt((m_cluster_sizeY + 1) / 2 * (m_cluster_sizeX + 1) / 2)),
c3(sqrt(m_cluster_sizeX * m_cluster_sizeY)),
m_pedestal(image_size[0], image_size[1]) {
c2 = sqrt((cluster_sizeY + 1) / 2 * (cluster_sizeX + 1) / 2);
c3 = sqrt(cluster_sizeX * cluster_sizeY);
};
// c2 = sqrt((cluster_sizeY + 1) / 2 * (cluster_sizeX + 1) / 2);
// c3 = sqrt(cluster_sizeX * cluster_sizeY);
};
template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
std::vector<Cluster> find_clusters_without_threshold(NDView<FRAME_TYPE, 2> frame, Pedestal<PEDESTAL_TYPE> &pedestal,
bool late_update = false) {
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
m_pedestal.push(frame);
}
NDArray<PEDESTAL_TYPE, 2> pedestal() {
return m_pedestal.mean();
}
std::vector<Cluster>
find_clusters_without_threshold(NDView<FRAME_TYPE, 2> frame,
// Pedestal<PEDESTAL_TYPE> &pedestal,
bool late_update = false) {
struct pedestal_update {
int x;
int y;
@ -52,10 +79,14 @@ class ClusterFinder {
long double total = 0;
eventMask[iy][ix] = PEDESTAL;
for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
val = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
for (short ir = -(m_cluster_sizeY / 2);
ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2);
ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
val = frame(iy + ir, ix + ic) -
m_pedestal.mean(iy + ir, ix + ic);
total += val;
if (val > max) {
max = val;
@ -63,9 +94,9 @@ class ClusterFinder {
}
}
}
auto rms = pedestal.standard_deviation(iy, ix);
auto rms = m_pedestal.std(iy, ix);
if (frame(iy, ix) - pedestal.mean(iy, ix) < -m_nSigma * rms) {
if (frame(iy, ix) - m_pedestal.mean(iy, ix) < -m_nSigma * rms) {
eventMask[iy][ix] = NEGATIVE_PEDESTAL;
continue;
} else if (max > m_nSigma * rms) {
@ -73,26 +104,33 @@ class ClusterFinder {
} else if (total > c3 * m_nSigma * rms) {
eventMask[iy][ix] = PHOTON;
} else{
} else {
if (late_update) {
pedestal_updates.push_back({ix, iy, frame(iy, ix)});
} else {
pedestal.push(iy, ix, frame(iy, ix));
m_pedestal.push(iy, ix, frame(iy, ix));
}
continue;
}
if (eventMask[iy][ix] == PHOTON && (frame(iy, ix) - pedestal.mean(iy, ix)) >= max) {
if (eventMask[iy][ix] == PHOTON &&
(frame(iy, ix) - m_pedestal.mean(iy, ix)) >= max) {
eventMask[iy][ix] = PHOTON_MAX;
Cluster cluster(m_cluster_sizeX, m_cluster_sizeY, Dtype(typeid(PEDESTAL_TYPE)));
Cluster cluster(m_cluster_sizeX, m_cluster_sizeY,
Dtype(typeid(PEDESTAL_TYPE)));
cluster.x = ix;
cluster.y = iy;
short i = 0;
for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
PEDESTAL_TYPE tmp = static_cast<PEDESTAL_TYPE>(frame(iy + ir, ix + ic)) -
pedestal.mean(iy + ir, ix + ic);
for (short ir = -(m_cluster_sizeY / 2);
ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2);
ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
PEDESTAL_TYPE tmp =
static_cast<PEDESTAL_TYPE>(
frame(iy + ir, ix + ic)) -
m_pedestal.mean(iy + ir, ix + ic);
cluster.set<PEDESTAL_TYPE>(i, tmp);
i++;
}
@ -104,14 +142,16 @@ class ClusterFinder {
}
if (late_update) {
for (auto &update : pedestal_updates) {
pedestal.push(update.y, update.x, update.value);
m_pedestal.push(update.y, update.x, update.value);
}
}
return clusters;
}
template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
std::vector<Cluster> find_clusters_with_threshold(NDView<FRAME_TYPE, 2> frame, Pedestal<PEDESTAL_TYPE> &pedestal) {
// template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
std::vector<Cluster>
find_clusters_with_threshold(NDView<FRAME_TYPE, 2> frame,
Pedestal<PEDESTAL_TYPE> &pedestal) {
assert(m_threshold > 0);
std::vector<Cluster> clusters;
std::vector<std::vector<eventType>> eventMask;
@ -123,7 +163,8 @@ class ClusterFinder {
NDArray<FRAME_TYPE, 2> rest({frame.shape(0), frame.shape(1)});
NDArray<int, 2> nph({frame.shape(0), frame.shape(1)});
// convert to n photons
// nph = (frame-pedestal.mean()+0.5*m_threshold)/m_threshold; // can be optimized with expression templates?
// nph = (frame-pedestal.mean()+0.5*m_threshold)/m_threshold; // can be
// optimized with expression templates?
for (int iy = 0; iy < frame.shape(0); iy++) {
for (int ix = 0; ix < frame.shape(1); ix++) {
auto val = frame(iy, ix) - pedestal.mean(iy, ix);
@ -145,10 +186,14 @@ class ClusterFinder {
}
eventMask[iy][ix] = NEIGHBOUR;
// iterate over cluster pixels around the current pixel (ix,iy)
for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
auto val = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
for (short ir = -(m_cluster_sizeY / 2);
ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2);
ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
auto val = frame(iy + ir, ix + ic) -
pedestal.mean(iy + ir, ix + ic);
total += val;
if (val > max) {
max = val;
@ -157,7 +202,7 @@ class ClusterFinder {
}
}
auto rms = pedestal.standard_deviation(iy, ix);
auto rms = pedestal.std(iy, ix);
if (m_nSigma == 0) {
tthr = m_threshold;
tthr1 = m_threshold;
@ -182,16 +227,22 @@ class ClusterFinder {
pedestal.push(iy, ix, frame(iy, ix));
continue;
}
if (eventMask[iy][ix] == PHOTON && frame(iy, ix) - pedestal.mean(iy, ix) >= max) {
if (eventMask[iy][ix] == PHOTON &&
frame(iy, ix) - pedestal.mean(iy, ix) >= max) {
eventMask[iy][ix] = PHOTON_MAX;
Cluster cluster(m_cluster_sizeX, m_cluster_sizeY, Dtype(typeid(FRAME_TYPE)));
Cluster cluster(m_cluster_sizeX, m_cluster_sizeY,
Dtype(typeid(FRAME_TYPE)));
cluster.x = ix;
cluster.y = iy;
short i = 0;
for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
auto tmp = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
for (short ir = -(m_cluster_sizeY / 2);
ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2);
ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
auto tmp = frame(iy + ir, ix + ic) -
pedestal.mean(iy + ir, ix + ic);
cluster.set<FRAME_TYPE>(i, tmp);
i++;
}
@ -203,14 +254,6 @@ class ClusterFinder {
}
return clusters;
}
protected:
int m_cluster_sizeX;
int m_cluster_sizeY;
double m_threshold;
double m_nSigma;
double c2;
double c3;
};
} // namespace aare

138
include/aare/Pedestal.hpp
View File

@ -6,11 +6,27 @@
namespace aare {
/**
* @brief Calculate the pedestal of a series of frames. Can be used as
* standalone but mostly used in the ClusterFinder.
*
* @tparam SUM_TYPE type of the sum
*/
template <typename SUM_TYPE = double> class Pedestal {
uint32_t m_rows;
uint32_t m_cols;
uint32_t m_samples;
NDArray<uint32_t, 2> m_cur_samples;
NDArray<SUM_TYPE, 2> m_sum;
NDArray<SUM_TYPE, 2> m_sum2;
public:
Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000)
: m_rows(rows), m_cols(cols), m_freeze(false), m_samples(n_samples), m_cur_samples(NDArray<uint32_t, 2>({rows, cols}, 0)),m_sum(NDArray<SUM_TYPE, 2>({rows, cols})),
m_sum2(NDArray<SUM_TYPE, 2>({rows, cols})) {
: m_rows(rows), m_cols(cols), m_samples(n_samples),
m_cur_samples(NDArray<uint32_t, 2>({rows, cols}, 0)),
m_sum(NDArray<SUM_TYPE, 2>({rows, cols})),
m_sum2(NDArray<SUM_TYPE, 2>({rows, cols})) {
assert(rows > 0 && cols > 0 && n_samples > 0);
m_sum = 0;
m_sum2 = 0;
@ -25,44 +41,51 @@ template <typename SUM_TYPE = double> class Pedestal {
return mean_array;
}
NDArray<SUM_TYPE, 2> variance() {
NDArray<SUM_TYPE, 2> variance_array({m_rows, m_cols});
for (uint32_t i = 0; i < m_rows * m_cols; i++) {
variance_array(i / m_cols, i % m_cols) = variance(i / m_cols, i % m_cols);
}
return variance_array;
}
NDArray<SUM_TYPE, 2> standard_deviation() {
NDArray<SUM_TYPE, 2> standard_deviation_array({m_rows, m_cols});
for (uint32_t i = 0; i < m_rows * m_cols; i++) {
standard_deviation_array(i / m_cols, i % m_cols) = standard_deviation(i / m_cols, i % m_cols);
}
return standard_deviation_array;
}
void clear() {
for (uint32_t i = 0; i < m_rows * m_cols; i++) {
clear(i / m_cols, i % m_cols);
}
}
/*
* index level operations
*/
SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
if (m_cur_samples(row, col) == 0) {
return 0.0;
}
return m_sum(row, col) / m_cur_samples(row, col);
}
NDArray<SUM_TYPE, 2> variance() {
NDArray<SUM_TYPE, 2> variance_array({m_rows, m_cols});
for (uint32_t i = 0; i < m_rows * m_cols; i++) {
variance_array(i / m_cols, i % m_cols) =
variance(i / m_cols, i % m_cols);
}
return variance_array;
}
SUM_TYPE variance(const uint32_t row, const uint32_t col) const {
if (m_cur_samples(row, col) == 0) {
return 0.0;
}
return m_sum2(row, col) / m_cur_samples(row, col) - mean(row, col) * mean(row, col);
return m_sum2(row, col) / m_cur_samples(row, col) -
mean(row, col) * mean(row, col);
}
SUM_TYPE standard_deviation(const uint32_t row, const uint32_t col) const { return std::sqrt(variance(row, col)); }
NDArray<SUM_TYPE, 2> std() {
NDArray<SUM_TYPE, 2> standard_deviation_array({m_rows, m_cols});
for (uint32_t i = 0; i < m_rows * m_cols; i++) {
standard_deviation_array(i / m_cols, i % m_cols) =
std(i / m_cols, i % m_cols);
}
return standard_deviation_array;
}
SUM_TYPE std(const uint32_t row, const uint32_t col) const {
return std::sqrt(variance(row, col));
}
void clear() {
for (uint32_t i = 0; i < m_rows * m_cols; i++) {
clear(i / m_cols, i % m_cols);
}
}
void clear(const uint32_t row, const uint32_t col) {
m_sum(row, col) = 0;
@ -72,29 +95,42 @@ template <typename SUM_TYPE = double> class Pedestal {
// frame level operations
template <typename T> void push(NDView<T, 2> frame) {
assert(frame.size() == m_rows * m_cols);
// TODO: test the effect of #pragma omp parallel for
for (uint32_t index = 0; index < m_rows * m_cols; index++) {
push<T>(index / m_cols, index % m_cols, frame(index));
// TODO! move away from m_rows, m_cols
if (frame.shape() != std::array<int64_t, 2>{m_rows, m_cols}) {
throw std::runtime_error(
"Frame shape does not match pedestal shape");
}
for (uint32_t row = 0; row < m_rows; row++) {
for (uint32_t col = 0; col < m_cols; col++) {
push<T>(row, col, frame(row, col));
}
}
// // TODO: test the effect of #pragma omp parallel for
// for (uint32_t index = 0; index < m_rows * m_cols; index++) {
// push<T>(index / m_cols, index % m_cols, frame(index));
// }
}
template <typename T> void push(Frame &frame) {
assert(frame.rows() == static_cast<size_t>(m_rows) && frame.cols() == static_cast<size_t>(m_cols));
assert(frame.rows() == static_cast<size_t>(m_rows) &&
frame.cols() == static_cast<size_t>(m_cols));
push<T>(frame.view<T>());
}
// getter functions
inline uint32_t rows() const { return m_rows; }
inline uint32_t cols() const { return m_cols; }
inline uint32_t n_samples() const { return m_samples; }
inline NDArray<uint32_t, 2> cur_samples() const { return m_cur_samples; }
inline NDArray<SUM_TYPE, 2> get_sum() const { return m_sum; }
inline NDArray<SUM_TYPE, 2> get_sum2() const { return m_sum2; }
uint32_t rows() const { return m_rows; }
uint32_t cols() const { return m_cols; }
uint32_t n_samples() const { return m_samples; }
NDArray<uint32_t, 2> cur_samples() const { return m_cur_samples; }
NDArray<SUM_TYPE, 2> get_sum() const { return m_sum; }
NDArray<SUM_TYPE, 2> get_sum2() const { return m_sum2; }
// pixel level operations (should be refactored to allow users to implement their own pixel level operations)
template <typename T> inline void push(const uint32_t row, const uint32_t col, const T val_) {
if (m_freeze) {
return;
}
// pixel level operations (should be refactored to allow users to implement
// their own pixel level operations)
template <typename T>
void push(const uint32_t row, const uint32_t col, const T val_) {
SUM_TYPE val = static_cast<SUM_TYPE>(val_);
const uint32_t idx = index(row, col);
if (m_cur_samples(idx) < m_samples) {
@ -106,16 +142,8 @@ template <typename SUM_TYPE = double> class Pedestal {
m_sum2(idx) += val * val - m_sum2(idx) / m_cur_samples(idx);
}
}
inline uint32_t index(const uint32_t row, const uint32_t col) const { return row * m_cols + col; };
void set_freeze(bool freeze) { m_freeze = freeze; }
private:
uint32_t m_rows;
uint32_t m_cols;
bool m_freeze;
uint32_t m_samples;
NDArray<uint32_t, 2> m_cur_samples;
NDArray<SUM_TYPE, 2> m_sum;
NDArray<SUM_TYPE, 2> m_sum2;
uint32_t index(const uint32_t row, const uint32_t col) const {
return row * m_cols + col;
};
};
} // namespace aare