Fixed Chunked Pedestal. Now should work as intended, giving sensible results compared to the previous version

This commit is contained in:
2025-08-11 16:44:21 +02:00
parent 7aa3fcfcd0
commit efb16ea8c1
3 changed files with 35 additions and 24 deletions

View File

@@ -40,6 +40,8 @@ template <typename SUM_TYPE = double> class ChunkedPedestal {
assert(rows > 0 && cols > 0 && chunk_size > 0); assert(rows > 0 && cols > 0 && chunk_size > 0);
m_mean = 0; m_mean = 0;
m_std = 0; m_std = 0;
m_current_frame_number = 0;
m_current_chunk_number = 0;
} }
~ChunkedPedestal() = default; ~ChunkedPedestal() = default;

View File

@@ -50,10 +50,10 @@ class ClusterFinder {
// 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);
// } // }
void push_pedestal_mean(NDView<FRAME_TYPE, 2> frame, uint32_t chunk_number) { void push_pedestal_mean(NDView<PEDESTAL_TYPE, 2> frame, uint32_t chunk_number) {
m_pedestal.push_mean(frame, chunk_number); m_pedestal.push_mean(frame, chunk_number);
} }
void push_pedestal_std(NDView<FRAME_TYPE, 2> frame, uint32_t chunk_number) { void push_pedestal_std(NDView<PEDESTAL_TYPE, 2> frame, uint32_t chunk_number) {
m_pedestal.push_std(frame, chunk_number); m_pedestal.push_std(frame, chunk_number);
} }
//This is here purely to keep the compiler happy for now //This is here purely to keep the compiler happy for now
@@ -96,11 +96,13 @@ class ClusterFinder {
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 rms = m_pedestal.std(iy, ix);
if (rms == 0) continue;
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? // 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)); PEDESTAL_TYPE value = (frame(iy, ix) - m_pedestal.mean(iy, ix));
if (value < -m_nSigma * rms) if (value < -m_nSigma * rms)
@@ -111,6 +113,8 @@ class ClusterFinder {
for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) { for (int ic = -dx; ic < dx + has_center_pixel_x; 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)) {
if (m_pedestal.std(iy + ir, ix + ic) == 0) continue;
PEDESTAL_TYPE val = PEDESTAL_TYPE val =
frame(iy + ir, ix + ic) - frame(iy + ir, ix + ic) -
m_pedestal.mean(iy + ir, ix + ic); m_pedestal.mean(iy + ir, ix + ic);
@@ -134,27 +138,30 @@ class ClusterFinder {
// m_pedestal.push_fast( // m_pedestal.push_fast(
// iy, ix, // iy, ix,
// frame(iy, // frame(iy,
// ix)); /std::cout << "max: " << max << std::endl; // ix)); // Assume we have reached n_samples in the
// // pedestal, slight performance improvement
continue; // It was a pedestal value nothing to store
}
// Store cluster // Store cluster
if (value == max) { if (value == max) {
/* // if (total < 0)
if (total < 0) // {
{ // std::cout << "" << std::endl;
// std::cout << ""; // std::cout << "frame_number: " << frame_number << std::endl;
std::cout << "ix: " << ix << std::endl; // std::cout << "ix: " << ix << std::endl;
std::cout << "iy: " << iy << std::endl; // std::cout << "iy: " << iy << std::endl;
std::cout << "frame(iy, ix): " << frame(iy, ix) << 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.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 << "m_pedestal.std(iy, ix): " << m_pedestal.std(iy, ix) << std::endl;
std::cout << "max: " << max << std::endl; // std::cout << "max: " << max << std::endl;
std::cout << "value: " << value << std::endl; // std::cout << "value: " << value << std::endl;
std::cout << "m_nSigma * rms: " << (m_nSigma * rms) << std::endl; // std::cout << "m_nSigma * rms: " << (m_nSigma * rms) << std::endl;
std::cout << "total: " << total << std::endl; // std::cout << "total: " << total << std::endl;
std::cout << "c3 * m_nSigma * rms: " << (c3 * m_nSigma * rms) << std::endl; // std::cout << "c3 * m_nSigma * rms: " << (c3 * m_nSigma * rms) << std::endl;
} // }
*/
ClusterType cluster{}; ClusterType cluster{};
cluster.x = ix; cluster.x = ix;
@@ -168,6 +175,8 @@ class ClusterFinder {
for (int ic = -dx; ic < dx + has_center_pixel_y; ic++) { for (int ic = -dx; ic < dx + has_center_pixel_y; 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)) {
if (m_pedestal.std(iy + ir, ix + ic) == 0) continue;
CT tmp = CT tmp =
static_cast<CT>(frame(iy + ir, ix + ic)) - static_cast<CT>(frame(iy + ir, ix + ic)) -
static_cast<CT>( static_cast<CT>(

View File

@@ -42,13 +42,13 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) {
.def("push_pedestal_mean", .def("push_pedestal_mean",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self, [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame, uint32_t chunk_number) { py::array_t<double> frame, uint32_t chunk_number) {
auto view = make_view_2d(frame); auto view = make_view_2d(frame);
self.push_pedestal_mean(view, chunk_number); self.push_pedestal_mean(view, chunk_number);
}) })
.def("push_pedestal_std", .def("push_pedestal_std",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self, [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame, uint32_t chunk_number) { py::array_t<double> frame, uint32_t chunk_number) {
auto view = make_view_2d(frame); auto view = make_view_2d(frame);
self.push_pedestal_std(view, chunk_number); self.push_pedestal_std(view, chunk_number);
}) })