diff --git a/include/ChunkedPedestal.hpp b/include/ChunkedPedestal.hpp new file mode 100644 index 0000000..1100eb0 --- /dev/null +++ b/include/ChunkedPedestal.hpp @@ -0,0 +1,131 @@ +#pragma once +#include "aare/Frame.hpp" +#include "aare/NDArray.hpp" +#include "aare/NDView.hpp" +#include +#include + +//JMulvey +//This is a new way to do pedestals (inspired by Dominic's cluster finder) +//Instead of pedestal tracking, we split the data (photon data) up into chunks (say 50K frames) +//For each chunk, we look at the spectra and fit to the noise peak. When we run the cluster finder, we then use this chunked pedestal data +//The smaller the chunk size, the more accurate, but also the longer it takes to process. +//It is essentially a pre-processing step. +//Ideally this new class will do that processing. +//But for now we will just implement a method to pass in the chunked pedestal values directly (I have my own script which does it for now) +//I've cut this down a lot, knowing full well it'll need changing if we want to merge it with main (happy to do that once I get it work for what I need) + +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 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 + + NDArray m_mean; + NDArray m_std; + 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({rows, cols, n_chunks})) { + assert(rows > 0 && cols > 0 && chunk_size > 0); + m_mean = 0; + m_std = 0; + } + ~Pedestal() = default; + + NDArray mean() { return m_mean; } + NDArray 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) + + if (chunk_number >= m_n_chunks) + { + chunk_number = 0; + throw std::runtime_error( + "Chunk number exceeds the number of chunks"); + } + } + + SUM_TYPE mean(const uint32_t row, const uint32_t col) const { + return m_mean(row, col, m_current_chunk_number); + } + + 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); + } + + void clear() { + m_mean = 0; + m_std = 0; + m_n_chunks = 0; + } + + //Probably don't need to do this one at a time, but let's keep it simple for now + template void push_mean(NDView frame, uint32_t chunk_number) { + assert(frame.size() == m_rows * m_cols); + + // TODO! move away from m_rows, m_cols + if (frame.shape() != std::array{m_rows, m_cols}) { + throw std::runtime_error( + "Frame shape does not match pedestal shape"); + } + + for (size_t row = 0; row < m_rows; row++) { + for (size_t col = 0; col < m_cols; col++) { + push_mean(row, col, chunk_number, frame(row, col)); + } + } + } + + template void push_std(NDView frame, uint32_t chunk_number) { + assert(frame.size() == m_rows * m_cols); + + // TODO! move away from m_rows, m_cols + if (frame.shape() != std::array{m_rows, m_cols}) { + throw std::runtime_error( + "Frame shape does not match pedestal shape"); + } + + for (size_t row = 0; row < m_rows; row++) { + for (size_t col = 0; col < m_cols; col++) { + push_std(row, col, chunk_number, frame(row, col)); + } + } + } + + // pixel level operations (should be refactored to allow users to implement + // their own pixel level operations) + template + 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_ + } + + template + 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_ + } + + // getter functions + uint32_t rows() const { return m_rows; } + uint32_t cols() const { return m_cols; } + +}; + +} // namespace aare \ No newline at end of file diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index 069d887..ae35740 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -4,7 +4,8 @@ #include "aare/Dtype.hpp" #include "aare/NDArray.hpp" #include "aare/NDView.hpp" -#include "aare/Pedestal.hpp" +// #include "aare/Pedestal.hpp" +#include "aare/ChunkedPedestal.hpp" #include "aare/defs.hpp" #include @@ -17,7 +18,7 @@ class ClusterFinder { const PEDESTAL_TYPE m_nSigma; const PEDESTAL_TYPE c2; const PEDESTAL_TYPE c3; - Pedestal m_pedestal; + ChunkedPedestal m_pedestal; ClusterVector m_clusters; static const uint8_t ClusterSizeX = ClusterType::cluster_size_x; @@ -34,18 +35,25 @@ class ClusterFinder { * */ ClusterFinder(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0, - size_t capacity = 1000000) + size_t capacity = 1000000, + uint32_t chunk_size = 50000, uint32_t n_chunks = 10) : m_image_size(image_size), m_nSigma(nSigma), c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)), c3(sqrt(ClusterSizeX * ClusterSizeY)), - m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) { + m_pedestal(image_size[0], image_size[1], chunk_size, n_chunks), m_clusters(capacity) { LOG(logDEBUG) << "ClusterFinder: " << "image_size: " << image_size[0] << "x" << image_size[1] << ", nSigma: " << nSigma << ", capacity: " << capacity; } - void push_pedestal_frame(NDView frame) { - m_pedestal.push(frame); + // void push_pedestal_frame(NDView frame) { + // m_pedestal.push(frame); + // } + void push_pedestal_mean(NDView frame, uint32_t chunk_number) { + m_pedestal.push_mean(frame, chunk_number); + } + void push_pedestal_std(NDView frame, uint32_t chunk_number) { + m_pedestal.push_std(frame, chunk_number); } NDArray pedestal() { return m_pedestal.mean(); } @@ -81,6 +89,7 @@ class ClusterFinder { int has_center_pixel_y = ClusterSizeY % 2; m_clusters.set_frame_number(frame_number); + m_pedestal.set_frame_number(frame_number); for (int iy = 0; iy < frame.shape(0); iy++) { for (int ix = 0; ix < frame.shape(1); ix++) { @@ -101,7 +110,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); + m_pedestal.mean(iy + ir, ix + ic, frame_number); total += val; max = std::max(max, val); @@ -117,11 +126,13 @@ class ClusterFinder { // pass } else { // m_pedestal.push(iy, ix, frame(iy, ix)); // Safe option - m_pedestal.push_fast( - iy, ix, - frame(iy, - ix)); // Assume we have reached n_samples in the - // pedestal, slight performance improvement + + //Not needed for chunked pedestal + // 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 }