diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index a98114d..8bd77cc 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -121,7 +121,8 @@ class ClusterFinder { } else if (total > c3 * m_nSigma * rms) { // pass } else { - m_pedestal.push(iy, ix, frame(iy, ix)); + // m_pedestal.push(iy, ix, frame(iy, ix)); + m_pedestal.push_fast(iy, ix, frame(iy, ix)); continue; // It was a pedestal value nothing to store } diff --git a/include/aare/Pedestal.hpp b/include/aare/Pedestal.hpp index bb2ea2c..bda94f2 100644 --- a/include/aare/Pedestal.hpp +++ b/include/aare/Pedestal.hpp @@ -98,7 +98,9 @@ template class Pedestal { m_sum2(row, col) = 0; m_cur_samples(row, col) = 0; } - // frame level operations + + + template void push(NDView frame) { assert(frame.size() == m_rows * m_cols); @@ -113,12 +115,32 @@ template class Pedestal { push(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(index / m_cols, index % m_cols, frame(index)); - // } } + + /** + * Push but don't update the cached mean. Speeds up the process + * when intitializing the pedestal. + * + */ + template void push_no_update(NDView frame) { + 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_no_update(row, col, frame(row, col)); + } + } + } + + + + template void push(Frame &frame) { assert(frame.rows() == static_cast(m_rows) && frame.cols() == static_cast(m_cols)); @@ -150,5 +172,36 @@ template class Pedestal { m_mean(row, col) = m_sum(row, col) / m_cur_samples(row, col); } + template + void push_no_update(const uint32_t row, const uint32_t col, const T val_) { + SUM_TYPE val = static_cast(val_); + if (m_cur_samples(row, col) < m_samples) { + m_sum(row, col) += val; + m_sum2(row, col) += val * val; + m_cur_samples(row, col)++; + } else { + m_sum(row, col) += val - m_sum(row, col) / m_cur_samples(row, col); + m_sum2(row, col) += val * val - m_sum2(row, col) / m_cur_samples(row, col); + } + } + + /** + * @brief Update the mean of the pedestal. This is used after having done + * push_no_update. It is not necessary to call this function after push. + */ + void update_mean(){ + m_mean = m_sum / m_cur_samples; + } + + template + void push_fast(const uint32_t row, const uint32_t col, const T val_){ + //Assume we reached the steady state where all pixels have + //m_samples samples + SUM_TYPE val = static_cast(val_); + m_sum(row, col) += val - m_sum(row, col) / m_samples; + m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples; + m_mean(row, col) = m_sum(row, col) / m_samples; + } + }; } // namespace aare \ No newline at end of file diff --git a/python/src/pedestal.hpp b/python/src/pedestal.hpp index 4d5d043..77148dc 100644 --- a/python/src/pedestal.hpp +++ b/python/src/pedestal.hpp @@ -43,5 +43,10 @@ template void define_pedestal_bindings(py::module &m, const .def("push", [](Pedestal &pedestal, py::array_t &f) { auto v = make_view_2d(f); pedestal.push(v); - }); + }) + .def("push_no_update", [](Pedestal &pedestal, py::array_t &f) { + auto v = make_view_2d(f); + pedestal.push_no_update(v); + }, py::arg().noconvert()) + .def("update_mean", &Pedestal::update_mean); } \ No newline at end of file