Bug Fixed for Chunked Pedestal

This commit is contained in:
2025-08-11 15:24:42 +02:00
parent 836dddbc26
commit 7aa3fcfcd0
4 changed files with 67 additions and 29 deletions

View File

@@ -2,7 +2,6 @@
#include "aare/Frame.hpp"
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
#include <fmath>
#include <cstddef>
//JMulvey
@@ -26,34 +25,34 @@ namespace aare {
template <typename SUM_TYPE = double> class ChunkedPedestal {
uint32_t m_rows;
uint32_t m_cols;
uint32_t m_n_chunks
uint64_t m_current_frame_number
uint64_t m_current_chunk_number
uint32_t m_n_chunks;
uint64_t m_current_frame_number;
uint64_t m_current_chunk_number;
NDArray<SUM_TYPE, 3> m_mean;
NDArray<SUM_TYPE, 3> m_std;
uint32_t m_chunk_size
uint32_t m_chunk_size;
public:
Pedestal(uint32_t rows, uint32_t cols, uint32_t chunk_size = 50000, uint32_t n_chunks = 10)
: m_rows(rows), m_cols(cols), m_chunk_size(chunk_size), m_n_chunks(n_chunks)
m_mean(NDArray<SUM_TYPE, 2>({rows, cols, n_chunks})) {
ChunkedPedestal(uint32_t rows, uint32_t cols, uint32_t chunk_size = 50000, uint32_t n_chunks = 10)
: m_rows(rows), m_cols(cols), m_chunk_size(chunk_size), m_n_chunks(n_chunks),
m_mean(NDArray<SUM_TYPE, 3>({rows, cols, n_chunks})), m_std(NDArray<SUM_TYPE, 3>({rows, cols, n_chunks})) {
assert(rows > 0 && cols > 0 && chunk_size > 0);
m_mean = 0;
m_std = 0;
}
~Pedestal() = default;
~ChunkedPedestal() = default;
NDArray<SUM_TYPE, 3> mean() { return m_mean; }
NDArray<SUM_TYPE, 3> std() { return m_std; }
void set_frame_number (uint64_t frame_number) {
m_current_frame_number = frame_number
uint32_t chunk_number = floor(frame_number / m_chunk_size)
m_current_frame_number = frame_number;
m_current_chunk_number = std::floor(frame_number / m_chunk_size);
if (chunk_number >= m_n_chunks)
if (m_current_chunk_number >= m_n_chunks)
{
chunk_number = 0;
m_current_chunk_number = 0;
throw std::runtime_error(
"Chunk number exceeds the number of chunks");
}
@@ -64,11 +63,7 @@ template <typename SUM_TYPE = double> class ChunkedPedestal {
}
SUM_TYPE std(const uint32_t row, const uint32_t col) const {
uint32_t chunk_number = floor(frame_number / m_chunk_size)
if (chunk_number >= m_n_chunks) return 0;
return m_std(row, col, chunk_number);
return m_std(row, col, m_current_chunk_number);
}
void clear() {
@@ -81,6 +76,10 @@ template <typename SUM_TYPE = double> class ChunkedPedestal {
template <typename T> void push_mean(NDView<T, 2> frame, uint32_t chunk_number) {
assert(frame.size() == m_rows * m_cols);
if (chunk_number >= m_n_chunks)
throw std::runtime_error(
"Chunk number is larger than the number of chunks");
// TODO! move away from m_rows, m_cols
if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
throw std::runtime_error(
@@ -92,11 +91,16 @@ template <typename SUM_TYPE = double> class ChunkedPedestal {
push_mean<T>(row, col, chunk_number, frame(row, col));
}
}
}
template <typename T> void push_std(NDView<T, 2> frame, uint32_t chunk_number) {
assert(frame.size() == m_rows * m_cols);
if (chunk_number >= m_n_chunks)
throw std::runtime_error(
"Chunk number is larger than the number of chunks");
// TODO! move away from m_rows, m_cols
if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
throw std::runtime_error(
@@ -108,18 +112,19 @@ template <typename SUM_TYPE = double> class ChunkedPedestal {
push_std<T>(row, col, chunk_number, frame(row, col));
}
}
}
// pixel level operations (should be refactored to allow users to implement
// their own pixel level operations)
template <typename T>
void push_mean(const uint32_t row, const uint32_t col, const uint32_t chunk_number, const T val_) {
m_mean(row, col, chunk_number) = val_
m_mean(row, col, chunk_number) = val_;
}
template <typename T>
void push_std(const uint32_t row, const uint32_t col, const uint32_t chunk_number, const T val_) {
m_std(row, col, chunk_number) = val_
m_std(row, col, chunk_number) = val_;
}
// getter functions

View File

@@ -8,6 +8,7 @@
#include "aare/ChunkedPedestal.hpp"
#include "aare/defs.hpp"
#include <cstddef>
#include <iostream>
namespace aare {
@@ -55,6 +56,8 @@ class ClusterFinder {
void push_pedestal_std(NDView<FRAME_TYPE, 2> frame, uint32_t chunk_number) {
m_pedestal.push_std(frame, chunk_number);
}
//This is here purely to keep the compiler happy for now
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {}
NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); }
NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
@@ -110,7 +113,7 @@ class ClusterFinder {
iy + ir >= 0 && iy + ir < frame.shape(0)) {
PEDESTAL_TYPE val =
frame(iy + ir, ix + ic) -
m_pedestal.mean(iy + ir, ix + ic, frame_number);
m_pedestal.mean(iy + ir, ix + ic);
total += val;
max = std::max(max, val);
@@ -131,13 +134,28 @@ class ClusterFinder {
// m_pedestal.push_fast(
// iy, ix,
// frame(iy,
// ix)); // Assume we have reached n_samples in the
// // pedestal, slight performance improvement
continue; // It was a pedestal value nothing to store
}
// ix)); /std::cout << "max: " << max << std::endl;
// Store cluster
if (value == max) {
/*
if (total < 0)
{
// std::cout << "";
std::cout << "ix: " << ix << std::endl;
std::cout << "iy: " << iy << std::endl;
std::cout << "frame(iy, ix): " << frame(iy, ix) << std::endl;
std::cout << "m_pedestal.mean(iy, ix): " << m_pedestal.mean(iy, ix) << std::endl;
std::cout << "m_pedestal.std(iy, ix): " << m_pedestal.std(iy, ix) << std::endl;
std::cout << "max: " << max << std::endl;
std::cout << "value: " << value << std::endl;
std::cout << "m_nSigma * rms: " << (m_nSigma * rms) << std::endl;
std::cout << "total: " << total << std::endl;
std::cout << "c3 * m_nSigma * rms: " << (c3 * m_nSigma * rms) << std::endl;
}
*/
ClusterType cluster{};
cluster.x = ix;
cluster.y = iy;

View File

@@ -26,13 +26,13 @@ def _get_class(name, cluster_size, dtype):
def ClusterFinder(image_size, cluster_size, n_sigma=5, dtype = np.int32, capacity = 1024):
def ClusterFinder(image_size, cluster_size, n_sigma=5, dtype = np.int32, capacity = 1024, chunk_size=50000, n_chunks = 10):
"""
Factory function to create a ClusterFinder object. Provides a cleaner syntax for
the templated ClusterFinder in C++.
"""
cls = _get_class("ClusterFinder", cluster_size, dtype)
return cls(image_size, n_sigma=n_sigma, capacity=capacity)
return cls(image_size, n_sigma=n_sigma, capacity=capacity, chunk_size=chunk_size, n_chunks=n_chunks)

View File

@@ -30,14 +30,29 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) {
py::class_<ClusterFinder<ClusterType, uint16_t, pd_type>>(
m, class_name.c_str())
.def(py::init<Shape<2>, pd_type, size_t>(), py::arg("image_size"),
py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000)
.def(py::init<Shape<2>, pd_type, size_t, uint32_t, uint32_t>(), py::arg("image_size"),
py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000,
py::arg("chunk_size") = 50'000, py::arg("n_chunks") = 10)
.def("push_pedestal_frame",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame) {
auto view = make_view_2d(frame);
self.push_pedestal_frame(view);
})
.def("push_pedestal_mean",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame, uint32_t chunk_number) {
auto view = make_view_2d(frame);
self.push_pedestal_mean(view, chunk_number);
})
.def("push_pedestal_std",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame, uint32_t chunk_number) {
auto view = make_view_2d(frame);
self.push_pedestal_std(view, chunk_number);
})
.def("clear_pedestal",
&ClusterFinder<ClusterType, uint16_t, pd_type>::clear_pedestal)
.def_property_readonly(