mirror of
https://github.com/slsdetectorgroup/aare.git
synced 2025-12-30 08:51:27 +01:00
217 lines
9.7 KiB
C++
217 lines
9.7 KiB
C++
#pragma once
|
|
#include "aare/ClusterFile.hpp"
|
|
#include "aare/ClusterVector.hpp"
|
|
#include "aare/Dtype.hpp"
|
|
#include "aare/NDArray.hpp"
|
|
#include "aare/NDView.hpp"
|
|
// #include "aare/Pedestal.hpp"
|
|
#include "aare/ChunkedPedestal.hpp"
|
|
#include "aare/defs.hpp"
|
|
#include <cstddef>
|
|
#include <iostream>
|
|
|
|
namespace aare {
|
|
|
|
template <typename ClusterType = Cluster<int32_t, 3, 3>,
|
|
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
|
|
class ClusterFinder {
|
|
Shape<2> m_image_size;
|
|
const PEDESTAL_TYPE m_nSigma;
|
|
const PEDESTAL_TYPE c2;
|
|
const PEDESTAL_TYPE c3;
|
|
ChunkedPedestal<PEDESTAL_TYPE> m_pedestal;
|
|
ClusterVector<ClusterType> m_clusters;
|
|
const uint32_t ClusterSizeX;
|
|
const uint32_t ClusterSizeY;
|
|
|
|
static const uint8_t SavedClusterSizeX = ClusterType::cluster_size_x;
|
|
static const uint8_t SavedClusterSizeY = ClusterType::cluster_size_y;
|
|
using CT = typename ClusterType::value_type;
|
|
|
|
public:
|
|
/**
|
|
* @brief Construct a new ClusterFinder object
|
|
* @param image_size size of the image
|
|
* @param cluster_size size of the cluster (x, y)
|
|
* @param nSigma number of sigma above the pedestal to consider a photon
|
|
* @param capacity initial capacity of the cluster vector
|
|
*
|
|
*/
|
|
ClusterFinder(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
|
|
size_t capacity = 1000000,
|
|
uint32_t chunk_size = 50000, uint32_t n_chunks = 10,
|
|
uint32_t cluster_size_x = 3, uint32_t cluster_size_y = 3)
|
|
: m_image_size(image_size), m_nSigma(nSigma),
|
|
c2(sqrt((cluster_size_y + 1) / 2 * (cluster_size_x + 1) / 2)),
|
|
c3(sqrt(cluster_size_x * cluster_size_y)),
|
|
ClusterSizeX(cluster_size_x), ClusterSizeY(cluster_size_y),
|
|
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_TYPE, 2> frame) {
|
|
// m_pedestal.push(frame);
|
|
// }
|
|
void push_pedestal_mean(NDView<PEDESTAL_TYPE, 2> frame, uint32_t chunk_number) {
|
|
m_pedestal.push_mean(frame, chunk_number);
|
|
}
|
|
void push_pedestal_std(NDView<PEDESTAL_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(); }
|
|
void clear_pedestal() { m_pedestal.clear(); }
|
|
|
|
/**
|
|
* @brief Move the clusters from the ClusterVector in the ClusterFinder to a
|
|
* new ClusterVector and return it.
|
|
* @param realloc_same_capacity if true the new ClusterVector will have the
|
|
* same capacity as the old one
|
|
*
|
|
*/
|
|
ClusterVector<ClusterType>
|
|
steal_clusters(bool realloc_same_capacity = false) {
|
|
ClusterVector<ClusterType> tmp = std::move(m_clusters);
|
|
if (realloc_same_capacity)
|
|
m_clusters = ClusterVector<ClusterType>(tmp.capacity());
|
|
else
|
|
m_clusters = ClusterVector<ClusterType>{};
|
|
return tmp;
|
|
}
|
|
void find_clusters(NDView<FRAME_TYPE, 2> frame, uint64_t frame_number = 0) {
|
|
// // currently 3,3 -> +/- 1
|
|
// // 4,4 -> +/- 2
|
|
int dy = ClusterSizeY / 2;
|
|
int dx = ClusterSizeX / 2;
|
|
int dy2 = SavedClusterSizeY / 2;
|
|
int dx2 = SavedClusterSizeX / 2;
|
|
|
|
int has_center_pixel_x =
|
|
ClusterSizeX %
|
|
2; // for even sized clusters there is no proper cluster center and
|
|
// even amount of pixels around the center
|
|
int has_center_pixel_y = ClusterSizeY % 2;
|
|
|
|
m_clusters.set_frame_number(frame_number);
|
|
m_pedestal.set_frame_number(frame_number);
|
|
auto mean_ptr = m_pedestal.get_mean_chunk_ptr();
|
|
auto std_ptr = m_pedestal.get_std_chunk_ptr();
|
|
|
|
for (int iy = 0; iy < frame.shape(0); iy++) {
|
|
size_t row_offset = iy * frame.shape(1);
|
|
for (int ix = 0; ix < frame.shape(1); ix++) {
|
|
|
|
// PEDESTAL_TYPE rms = m_pedestal.std(iy, ix);
|
|
PEDESTAL_TYPE rms = std_ptr[row_offset + ix];
|
|
if (rms == 0) continue;
|
|
|
|
PEDESTAL_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
|
|
PEDESTAL_TYPE total = 0;
|
|
|
|
// What can we short circuit here?
|
|
// PEDESTAL_TYPE value = (frame(iy, ix) - m_pedestal.mean(iy, ix));
|
|
PEDESTAL_TYPE value = (frame(iy, ix) - mean_ptr[row_offset + ix]);
|
|
|
|
if (value < -m_nSigma * rms)
|
|
continue; // NEGATIVE_PEDESTAL go to next pixel
|
|
// TODO! No pedestal update???
|
|
|
|
for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
|
|
size_t inner_row_offset = row_offset + (ir * frame.shape(1));
|
|
for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
|
|
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
|
|
iy + ir >= 0 && iy + ir < frame.shape(0)) {
|
|
// if (m_pedestal.std(iy + ir, ix + ic) == 0) continue;
|
|
if (std_ptr[inner_row_offset + ix + ic] == 0) continue;
|
|
|
|
// PEDESTAL_TYPE val = frame(iy + ir, ix + ic) - m_pedestal.mean(iy + ir, ix + ic);
|
|
PEDESTAL_TYPE val = frame(iy + ir, ix + ic) - mean_ptr[inner_row_offset + ix + ic];
|
|
|
|
total += val;
|
|
max = std::max(max, val);
|
|
}
|
|
}
|
|
}
|
|
|
|
// if ((max > m_nSigma * rms)) {
|
|
// if (value < max)
|
|
// continue; // Not max go to the next pixel
|
|
// // but also no pedestal update
|
|
// } else
|
|
if (total > c3 * m_nSigma * rms) {
|
|
// pass
|
|
} else {
|
|
// m_pedestal.push(iy, ix, frame(iy, ix)); // Safe option
|
|
|
|
//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
|
|
|
|
}
|
|
|
|
// Store cluster
|
|
if (value == max) {
|
|
|
|
// if (total < 0)
|
|
// {
|
|
// std::cout << "" << std::endl;
|
|
// std::cout << "frame_number: " << frame_number << std::endl;
|
|
// 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;
|
|
|
|
// Fill the cluster data since we have a photon to store
|
|
// It's worth redoing the look since most of the time we
|
|
// don't have a photon
|
|
int i = 0;
|
|
for (int ir = -dy2; ir < dy2 + has_center_pixel_y; ir++) {
|
|
size_t inner_row_offset = row_offset + (ir * frame.shape(1));
|
|
for (int ic = -dx2; ic < dx2 + has_center_pixel_y; ic++) {
|
|
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
|
|
iy + ir >= 0 && iy + ir < frame.shape(0)) {
|
|
// if (m_pedestal.std(iy + ir, ix + ic) == 0) continue;
|
|
// if (std_ptr[inner_row_offset + ix + ic] == 0) continue;
|
|
|
|
// CT tmp = static_cast<CT>(frame(iy + ir, ix + ic)) - static_cast<CT>(m_pedestal.mean(iy + ir, ix + ic));
|
|
|
|
// CT tmp = 0;
|
|
if (std_ptr[inner_row_offset + ix + ic] != 0)
|
|
{
|
|
CT tmp = static_cast<CT>(frame(iy + ir, ix + ic)) - static_cast<CT>(mean_ptr[inner_row_offset + ix + ic]);
|
|
cluster.data[i] = tmp; // Watch for out of bounds access
|
|
}
|
|
}
|
|
i++;
|
|
}
|
|
}
|
|
|
|
// Add the cluster to the output ClusterVector
|
|
m_clusters.push_back(cluster);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
};
|
|
|
|
} // namespace aare
|