detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2025-06-06 12:50:40 +02:00
Code Issues Actions Packages Releases Activity

cluster finder (#72)

Browse Source 
implement fixed size cluster finder algorithm 
clusterFile functions are commented and put on hold.

---------

Co-authored-by: Bechir <bechir.brahem420@gmail.com>
Co-authored-by: Erik Fröjdh <erik.frojdh@gmail.com>
This commit is contained in:
Bechir Braham 2024-06-17 10:57:23 +02:00 committed by GitHub
parent 1dd361a343
commit 5b7ab5a810
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
33 changed files with 845 additions and 111 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
.CMakeLists.txt.swp Normal file
View File

Binary file not shown.

2
CMakeLists.txt
View File

@ -1,4 +1,3 @@
cmake_minimum_required(VERSION 3.12)
set(CMAKE_CXX_STANDARD 17) #TODO! Global or per target?
@ -144,6 +143,7 @@ else()
-fno-sanitize-recover
# -D_FORTIFY_SOURCE=2
-fno-omit-frame-pointer
-lstdc++fs
)
endif()

2
examples/CMakeLists.txt
View File

@ -2,7 +2,7 @@
set(EXAMPLE_LIST "json_example;logger_example;numpy_read_example;multiport_example;raw_example;zmq_restream_example")
set(EXAMPLE_LIST "${EXAMPLE_LIST};mythen_example;numpy_write_example;zmq_receiver_example;zmq_sender_example;")
set(EXAMPLE_LIST "${EXAMPLE_LIST};cluster_example;zmq_multi_receiver;zmq_task_ventilator;zmq_worker;zmq_sink")
set(EXAMPLE_LIST "${EXAMPLE_LIST};testing_pedestal;")
set(EXAMPLE_LIST "${EXAMPLE_LIST};testing_pedestal;cluster_finder_example;")
foreach(example ${EXAMPLE_LIST})
add_executable(${example} ${example}.cpp)

57
examples/cluster_example.cpp
View File

@ -7,48 +7,33 @@
using namespace aare;
int main() {
auto PROJECT_ROOT_DIR = std::filesystem::path(getenv("AARE_ROOT_DIR"));
std::filesystem::path const fpath(PROJECT_ROOT_DIR / "data" / "clusters" / "single_frame_97_clustrers.clust");
std::filesystem::path const fpath("/mnt/sls_det_storage/moench_data/testNewFW20230714/cu_half_speed_master_4.json");
// reading a file
aare::ClusterFile cf(fpath, "r");
std::cout << "file opened " << '\n';
std::cout << "n_clusters " << cf.count() << '\n';
std::cout << "frame_number " << cf.frame() << '\n';
std::cout << "file size: " << cf.count() << '\n';
cf.seek(0); // seek to the beginning of the file (this is the default behavior of the constructor)
NDArray<double, 2> frame({10, 10});
frame = 0;
auto cluster = cf.read(97);
std::cout << "read 10 clusters" << '\n';
int offset = 0;
int data_offset = 0;
for (auto c : cluster) {
assert(c.y == offset + 200);
for (int i = 0; i < 9; i++) {
assert(c.data[i] == data_offset + i);
for (int i = 5; i < 8; i++) {
for (int j = 5; j < 8; j++) {
frame(i, j) = (i + j) % 10;
}
offset++;
data_offset += 9;
}
// writing a file
std::filesystem::path const fpath_out("/tmp/cluster_example_file.clust");
aare::ClusterFile cf_out(fpath_out, "w", ClusterFileConfig(1, 44));
std::cout << "file opened for writing" << '\n';
std::vector<Cluster> clusters;
for (int i = 0; i < 1084; i++) {
Cluster c;
c.x = i;
c.y = i + 200;
for (int j = 0; j < 9; j++) {
c.data[j] = j;
for (int i = 0; i < frame.shape(0); i++) {
for (int j = 0; j < frame.shape(1); j++) {
std::cout << frame(i, j) << " ";
}
clusters.push_back(c);
std::cout << std::endl;
}
cf_out.write(clusters);
std::cout << "wrote 10 clusters" << '\n';
cf_out.update_header();
return 0;
ClusterFinder clusterFinder(3, 3, 1, 1); // 3x3 cluster, 1 nSigma, 1 threshold
Pedestal p(10, 10);
auto clusters = clusterFinder.find_clusters(frame.span(), p);
aare::logger::info("nclusters:", clusters.size());
for (auto &cluster : clusters) {
aare::logger::info("cluster center:", cluster.to_string<double>());
}
}

63
examples/cluster_finder_example.cpp Normal file
View File

@ -0,0 +1,63 @@
#include "aare.hpp"
#include "aare/examples/defs.hpp"
#include <cassert>
#include <iostream>
using namespace aare;
int main() {
auto PROJECT_ROOT_DIR = std::filesystem::path(getenv("AARE_ROOT_DIR"));
std::filesystem::path const fpath_frame("/home/l_bechir/tmp/testNewFW20230714/cu_half_speed_master_4.json");
std::filesystem::path const fpath_cluster("/home/l_bechir/tmp/testNewFW20230714/clust/cu_half_speed_d0_f0_4.clust");
auto file = RawFile(fpath_frame, "r");
logger::info("RAW file");
logger::info("file rows:", file.rows(), "cols:", file.cols());
logger::info(file.total_frames());
Frame frame(0, 0, 0);
for (auto i = 0; i < 10000; i++) {
if (file.frame_number(i) == 23389) {
logger::info("frame_number:", file.frame_number(i));
frame = file.read();
}
}
logger::info("frame", frame.rows(), frame.cols(), frame.bitdepth());
ClusterFileV2 cf(fpath_cluster, "r");
auto anna_clusters = cf.read();
logger::info("Cluster file");
logger::info("nclusters:", anna_clusters.size());
logger::info("frame_number", anna_clusters[0].frame_number);
ClusterFinder clusterFinder(3, 3, 5, 0);
Pedestal p(file.rows(), file.cols());
file.seek(0);
logger::info("Starting Pedestal calculation");
for (auto i = 0; i < 1000; i++) {
p.push(file.read().view<uint16_t>());
}
logger::info("Pedestal calculation done");
logger::info("Pedestal mean:", p.mean(0, 0), "std:", p.standard_deviation(0, 0));
logger::info("Pedestal mean:", p.mean(200, 200), "std:", p.standard_deviation(200, 200));
FileConfig cfg;
cfg.dtype = DType(typeid(double));
cfg.rows = p.rows();
cfg.cols = p.cols();
NumpyFile np_pedestal("/home/l_bechir/tmp/testNewFW20230714/pedestal.npy", "w", cfg);
cfg.dtype = DType(typeid(uint16_t));
NumpyFile np_frame("/home/l_bechir/tmp/testNewFW20230714/frame.npy", "w", cfg);
np_pedestal.write(p.mean());
np_frame.write(frame.view<uint16_t>());
auto clusters = clusterFinder.find_clusters(frame.view<uint16_t>(), p);
logger::info("nclusters:", clusters.size());
// aare::logger::info("nclusters:", clusters.size());
// for (auto &cluster : clusters) {
// aare::logger::info("cluster center:", cluster.to_string<double>());
// }
}

7
extra/.gitignore vendored
View File

@ -1 +1,6 @@
sqlite_clusters.db
fixed_size_clusters.bin
numpy_clusters.npy
sqlite_clusters.zip
fixed_size_clusters.zip
numpy_clusters.zip

179
extra/compare_file_formats.py Normal file
View File

@ -0,0 +1,179 @@
import numpy as np
import zipfile
from zipfile import ZipFile as zf
import time
from pathlib import Path
import sqlite3
def timing_val(f):
def wrapper(*arg, **kw):
t1 = time.time()
res = f(*arg, **kw)
t2 = time.time()
return (t2 - t1), res, f.__name__
return wrapper
N_CLUSTERS = 1_000_000
"""
fixed size clusters:
header:
- magic_string: 4 bytes
- version: 1 byte
- n_records: 4 bytes
- indexed: 1 byte
- metadata_length: 1 byte (number of chars)
- metadata: metadata_length bytes (json string)
- field_count: 1 byte
- fields:
- field_label_length: 1 byte
- field_label: field_label_length bytes (string)
- dtype: 3 bytes (string)
- is_array: 1 byte (0: not array, 1:fixed_length_array, 2:variable_length_array)
- array_length: 4 bytes (used if is_array == 1)
data:
- field: (field_1_dtype_length bytes) or
"""
header_length = 4 + 1 + 4 + 1 + 1 + 2 + 1 + (1 + 5 + 3 + 1 + 4) * 3
cluster_dt = np.dtype([("x", "int16"), ("y", "int16"), ("data", "int32", (3, 3))])
def write_binary_clusters():
with open("fixed_size_clusters.bin", "wb") as f:
f.write(b"H" * header_length)
arr = np.zeros(N_CLUSTERS, dtype=cluster_dt)
f.write(arr.tobytes())
def write_numpy_clusters():
np.save("numpy_clusters.npy", np.zeros(N_CLUSTERS, dtype=cluster_dt))
def write_sqlite_clusters():
data = np.zeros(9, dtype=np.int32).tobytes()
c = conn.cursor()
c.execute("CREATE TABLE clusters (x int, y int, data blob)")
c.executemany("INSERT INTO clusters VALUES (?, ?, ?)", [(0, 0, data)] * N_CLUSTERS)
conn.commit()
READ_N_CLUSTERS = N_CLUSTERS
def read_binary_clusters():
with open("fixed_size_clusters.bin", "rb") as f:
f.read(header_length)
f.read(READ_N_CLUSTERS * cluster_dt.itemsize)
def read_numpy_clusters():
arr = np.load("numpy_clusters.npy")
def read_sqlite_clusters():
c = conn.cursor()
c.execute("SELECT * FROM clusters LIMIT ?", (READ_N_CLUSTERS,))
arr = c.fetchall()
N_APPEND_CLUSTERS = 100_000
def append_binary_clusters():
with open("fixed_size_clusters.bin", "ab") as f:
arr = np.zeros(N_APPEND_CLUSTERS, dtype=cluster_dt)
f.write(arr.tobytes())
def append_sqlite_clusters():
data = np.zeros(9, dtype=np.int32).tobytes()
c = conn.cursor()
c.executemany(
"INSERT INTO clusters VALUES (?, ?, ?)", [(0, 0, data)] * N_APPEND_CLUSTERS
)
conn.commit()
def p(write_time, file_size):
file_size = file_size / 1024 / 1024
print("%.3fs" % write_time, "%.3fMB" % file_size)
if __name__ == "__main__":
# setup
Path("fixed_size_clusters.bin").unlink(missing_ok=True)
Path("numpy_clusters.npy").unlink(missing_ok=True)
Path("sqlite_clusters.db").unlink(missing_ok=True)
Path("fixed_size_clusters.zip").unlink(missing_ok=True)
Path("numpy_clusters.zip").unlink(missing_ok=True)
Path("sqlite_clusters.zip").unlink(missing_ok=True)
conn = sqlite3.connect("sqlite_clusters.db")
# run
print("Testing file creation", f"(N_CLUSTERS={N_CLUSTERS}):")
print("Binary clusters:", end=" ")
bin_time, _, _ = timing_val(write_binary_clusters)()
bin_size = Path("fixed_size_clusters.bin").stat().st_size
p(bin_time, bin_size)
print("Numpy clusters:", end=" ")
np_time, _, _ = timing_val(write_numpy_clusters)()
np_size = Path("numpy_clusters.npy").stat().st_size
p(np_time, np_size)
print("SQLite clusters:", end=" ")
sql_time, _, _ = timing_val(write_sqlite_clusters)()
sql_size = Path("sqlite_clusters.db").stat().st_size
p(sql_time, sql_size)
print("\nTesting file reading", f"(READ_N_CLUSTERS={READ_N_CLUSTERS}):")
print("Binary clusters:", end=" ")
print("%.5fs" % timing_val(read_binary_clusters)()[0])
print("Numpy clusters:", end=" ")
print("%.5fs" % timing_val(read_numpy_clusters)()[0])
print("SQLite clusters:", end=" ")
print("%.5fs" % timing_val(read_sqlite_clusters)()[0])
print("\nTesting appending to file:")
print("Binary clusters:", end=" ")
print("%.5fs" % timing_val(append_binary_clusters)()[0])
print("SQLite clusters:", end=" ")
print("%.5fs" % timing_val(append_sqlite_clusters)()[0])
print("\nTesting zip compression:")
print("Binary clusters compressed:", end=" ")
with zf("fixed_size_clusters.zip", "w", zipfile.ZIP_DEFLATED) as z:
z.write("fixed_size_clusters.bin")
print(
"%.3fMB" % (Path("fixed_size_clusters.zip").stat().st_size / 1024 / 1024),
end=" ",
)
rate = (1 - Path("fixed_size_clusters.zip").stat().st_size / bin_size) * 100
print("rate:", "%.2f" % rate + "%")
print("Numpy clusters compressed:", end=" ")
with zf("numpy_clusters.zip", "w", zipfile.ZIP_DEFLATED) as z:
z.write("numpy_clusters.npy")
print("%.3fMB" % (Path("numpy_clusters.zip").stat().st_size / 1024 / 1024), end=" ")
rate = (1 - Path("numpy_clusters.zip").stat().st_size / bin_size) * 100
print("rate:", "%.2f" % rate + "%")
print("SQLite clusters compressed:", end=" ")
with zf("sqlite_clusters.zip", "w", zipfile.ZIP_DEFLATED) as z:
z.write("sqlite_clusters.db")
print(
"%.3fMB" % (Path("sqlite_clusters.zip").stat().st_size / 1024 / 1024), end=" "
)
rate = (1 - Path("sqlite_clusters.zip").stat().st_size / bin_size) * 100
print("rate:", "%.2f" % rate + "%")
# clean
conn.close()
# Path("fixed_size_clusters.bin").unlink(missing_ok=True)
# Path("numpy_clusters.npy").unlink(missing_ok=True)
# Path("sqlite_clusters.db").unlink(missing_ok=True)

2
include/aare.hpp
View File

@ -1,6 +1,8 @@
// This is the top level header to include and what most users will use
// include all header files
// #define private public
// #define protected public
#include "aare/core.hpp"
#include "aare/file_io.hpp"
#include "aare/network_io.hpp"

1
include/aare/file_io.hpp
View File

@ -1,4 +1,5 @@
#include "aare/file_io/ClusterFile.hpp"
#include "aare/file_io/ClusterFileV2.hpp"
#include "aare/file_io/File.hpp"
#include "aare/file_io/FileInterface.hpp"
#include "aare/file_io/NumpyFile.hpp"

1
include/aare/processing.hpp
View File

@ -1 +1,2 @@
#include "aare/processing/ClusterFinder.hpp"
#include "aare/processing/Pedestal.hpp"

3
src/core/include/aare/core/DType.hpp
View File

@ -32,6 +32,7 @@ class DType {
enum TypeIndex { INT8, UINT8, INT16, UINT16, INT32, UINT32, INT64, UINT64, FLOAT, DOUBLE, ERROR };
uint8_t bitdepth() const;
uint8_t bytes() const;
explicit DType(const std::type_info &t);
explicit DType(std::string_view sv);
@ -46,7 +47,7 @@ class DType {
// bool operator==(DType::TypeIndex ti) const;
// bool operator!=(DType::TypeIndex ti) const;
std::string str() const;
std::string to_string() const;
private:
TypeIndex m_type{TypeIndex::ERROR};

3
src/core/include/aare/core/Frame.hpp
View File

@ -58,6 +58,9 @@ class Frame {
m_rows = other.rows();
m_cols = other.cols();
m_bitdepth = other.bitdepth();
if (m_data != nullptr) {
delete[] m_data;
}
m_data = other.m_data;
other.m_data = nullptr;
other.m_rows = other.m_cols = other.m_bitdepth = 0;

1
src/core/include/aare/core/NDView.hpp
View File

@ -67,6 +67,7 @@ template <typename T, ssize_t Ndim = 2> class NDView {
}
ssize_t size() const { return size_; }
size_t total_bytes() const { return size_ * sizeof(T); }
T *begin() { return buffer_; }
T *end() { return buffer_ + size_; }

72
src/core/include/aare/core/defs.hpp
View File

@ -1,9 +1,14 @@
#pragma once
#include "aare/core/DType.hpp"
#include "aare/utils/logger.hpp"
#include <array>
#include <stdexcept>
#include <cassert>
#include <cstdint>
#include <cstring>
#include <string>
#include <string_view>
#include <variant>
@ -11,18 +16,73 @@
namespace aare {
struct Cluster {
class Cluster {
public:
int cluster_sizeX;
int cluster_sizeY;
int16_t x;
int16_t y;
std::array<int32_t, 9> data;
std::string to_string() const {
std::string s = "x: " + std::to_string(x) + " y: " + std::to_string(y) + "\ndata: [";
for (auto d : data) {
s += std::to_string(d) + " ";
DType dt;
private:
std::byte *m_data;
public:
Cluster(int cluster_sizeX_, int cluster_sizeY_, DType dt_ = DType(typeid(int32_t)))
: cluster_sizeX(cluster_sizeX_), cluster_sizeY(cluster_sizeY_), dt(dt_) {
m_data = new std::byte[cluster_sizeX * cluster_sizeY * dt.bytes()]{};
}
Cluster() : Cluster(3, 3) {}
Cluster(const Cluster &other) : Cluster(other.cluster_sizeX, other.cluster_sizeY, other.dt) {
if (this == &other)
return;
x = other.x;
y = other.y;
memcpy(m_data, other.m_data, other.size());
}
Cluster &operator=(const Cluster &other) {
if (this == &other)
return *this;
this->~Cluster();
new (this) Cluster(other);
return *this;
}
Cluster(Cluster &&other) noexcept
: cluster_sizeX(other.cluster_sizeX), cluster_sizeY(other.cluster_sizeY), x(other.x), y(other.y), dt(other.dt),
m_data(other.m_data) {
other.m_data = nullptr;
other.dt = DType(DType::TypeIndex::ERROR);
}
~Cluster() { delete[] m_data; }
template <typename T> T get(int idx) {
(sizeof(T) == dt.bytes()) ? 0 : throw std::invalid_argument("[ERROR] Type size mismatch");
return *reinterpret_cast<T *>(m_data + idx * dt.bytes());
}
template <typename T> auto set(int idx, T val) {
(sizeof(T) == dt.bytes()) ? 0 : throw std::invalid_argument("[ERROR] Type size mismatch");
return memcpy(m_data + idx * dt.bytes(), &val, (size_t)dt.bytes());
}
// auto x() const { return x; }
// auto y() const { return y; }
// auto x(int16_t x_) { return x = x_; }
// auto y(int16_t y_) { return y = y_; }
template <typename T> std::string to_string() const {
(sizeof(T) == dt.bytes()) ? 0 : throw std::invalid_argument("[ERROR] Type size mismatch");
std::string s = "x: " + std::to_string(x) + " y: " + std::to_string(y) + "\nm_data: [";
for (int i = 0; i < cluster_sizeX * cluster_sizeY; i++) {
s += std::to_string(*reinterpret_cast<T *>(m_data + i * dt.bytes())) + " ";
}
s += "]";
return s;
}
/**
* @brief size of the cluster in bytes when saved to a file
*/
size_t size() const { return cluster_sizeX * cluster_sizeY * dt.bytes(); }
auto begin() const { return m_data; }
auto end() const { return m_data + cluster_sizeX * cluster_sizeY * dt.bytes(); }
std::byte *data() { return m_data; }
};
/**

7
src/core/src/DType.cpp
View File

@ -66,6 +66,11 @@ uint8_t DType::bitdepth() const {
}
}
/**
* @brief Get the number of bytes of the data type
*/
uint8_t DType::bytes() const { return bitdepth() / 8; }
/**
* @brief Construct a DType object from a TypeIndex
* @param ti TypeIndex
@ -125,7 +130,7 @@ DType::DType(std::string_view sv) {
* @brief Get the string representation of the data type
* @return string representation
*/
std::string DType::str() const {
std::string DType::to_string() const {
char ec{};
if (endian::native == endian::little)

2
src/core/test/DType.test.cpp
View File

@ -51,4 +51,4 @@ TEST_CASE("Construct from string with endianess") {
REQUIRE_THROWS(DType(">i4") == typeid(int32_t));
}
TEST_CASE("Convert to string") { REQUIRE(DType(typeid(int)).str() == "<i4"); }
TEST_CASE("Convert to string") { REQUIRE(DType(typeid(int)).to_string() == "<i4"); }

33
src/core/test/defs.test.cpp
View File

@ -1,8 +1,41 @@
#include "aare/core/defs.hpp"
#include "aare/utils/floats.hpp"
#include <catch2/catch_test_macros.hpp>
#include <string>
TEST_CASE("Enum to string conversion") {
// By the way I don't think the enum string conversions should be in the defs.hpp file
// but let's use this to show a test
REQUIRE(toString(aare::DetectorType::Jungfrau) == "Jungfrau");
}
TEST_CASE("Cluster creation") {
aare::Cluster c(13, 15);
REQUIRE(c.cluster_sizeX == 13);
REQUIRE(c.cluster_sizeY == 15);
REQUIRE(c.dt == aare::DType(typeid(int32_t)));
REQUIRE(c.data() != nullptr);
aare::Cluster c2(c);
REQUIRE(c2.cluster_sizeX == 13);
REQUIRE(c2.cluster_sizeY == 15);
REQUIRE(c2.dt == aare::DType(typeid(int32_t)));
REQUIRE(c2.data() != nullptr);
}
TEST_CASE("cluster set and get data") {
aare::Cluster c2(33, 44, aare::DType(typeid(double)));
REQUIRE(c2.cluster_sizeX == 33);
REQUIRE(c2.cluster_sizeY == 44);
REQUIRE(c2.dt == aare::DType::DOUBLE);
double v = 3.14;
c2.set<double>(0, v);
double v2 = c2.get<double>(0);
REQUIRE(aare::compare_floats<double>(v, v2));
c2.set<double>(33 * 44 - 1, 123.11);
double v3 = c2.get<double>(33 * 44 - 1);
REQUIRE(aare::compare_floats<double>(123.11, v3));
REQUIRE_THROWS_AS(c2.set(0, 1), std::invalid_argument); // set int to double
}

27
src/file_io/include/aare/file_io/ClusterFile.hpp
View File

@ -12,7 +12,7 @@
* typedef struct {
* int16_t x;
* int16_t y;
* int32_t data[9];
* int32_t data[n];
*} Cluster ;
*
*/
@ -26,18 +26,31 @@ namespace aare {
struct ClusterFileConfig {
int32_t frame_number;
int32_t n_clusters;
ClusterFileConfig(int32_t frame_number_, int32_t n_clusters_)
: frame_number(frame_number_), n_clusters(n_clusters_) {}
ClusterFileConfig() : frame_number(0), n_clusters(0) {}
int cluster_sizeX;
int cluster_sizeY;
DType dt;
ClusterFileConfig(int32_t frame_number_, int32_t n_clusters_, int cluster_sizeX_ = 3, int cluster_sizeY_ = 3,
DType dt_ = DType::INT32)
: frame_number{frame_number_}, n_clusters{n_clusters_}, cluster_sizeX{cluster_sizeX_},
cluster_sizeY{cluster_sizeY_}, dt{dt_} {}
ClusterFileConfig() : ClusterFileConfig(0, 0) {}
bool operator==(const ClusterFileConfig &other) const {
return frame_number == other.frame_number && n_clusters == other.n_clusters;
return frame_number == other.frame_number && n_clusters == other.n_clusters && dt == other.dt &&
cluster_sizeX == other.cluster_sizeX && cluster_sizeY == other.cluster_sizeY;
}
bool operator!=(const ClusterFileConfig &other) const { return !(*this == other); }
std::string to_string() const {
return "frame_number: " + std::to_string(frame_number) + " n_clusters: " + std::to_string(n_clusters) + "\n";
return "frame_number: " + std::to_string(frame_number) + ", n_clusters: " + std::to_string(n_clusters) +
", dt: " + dt.to_string() + "\n cluster_sizeX: " + std::to_string(cluster_sizeX) +
", cluster_sizeY: " + std::to_string(cluster_sizeY);
}
};
struct ClusterFileHeader {
int32_t frame_number;
int32_t n_clusters;
};
/**
* @brief Class to read and write clusters to a file
*/
@ -65,6 +78,8 @@ class ClusterFile {
int32_t frame_number{};
int32_t n_clusters{};
static const int HEADER_BYTES = 8;
DType dt;
size_t m_cluster_size{};
};
} // namespace aare

68
src/file_io/include/aare/file_io/ClusterFileV2.hpp Normal file
View File

@ -0,0 +1,68 @@
#pragma once
#include "aare/core/defs.hpp"
#include <filesystem>
#include <string>
namespace aare {
struct ClusterHeader {
int32_t frame_number;
int32_t n_clusters;
};
struct ClusterV2_ {
int16_t x;
int16_t y;
std::array<int32_t, 9> data;
};
struct ClusterV2 {
ClusterV2_ cluster;
int16_t frame_number;
};
class ClusterFileV2 {
private:
bool m_closed = true;
std::filesystem::path m_fpath;
std::string m_mode;
FILE *fp;
public:
ClusterFileV2(std::filesystem::path const &fpath, std::string const &mode) {
m_fpath = fpath;
m_mode = mode;
fp = fopen(fpath.c_str(), "rb");
m_closed = false;
}
~ClusterFileV2() { close(); }
std::vector<ClusterV2> read() {
ClusterHeader header;
fread(&header, sizeof(ClusterHeader), 1, fp);
std::vector<ClusterV2_> clusters_(header.n_clusters);
fread(clusters_.data(), sizeof(ClusterV2_), header.n_clusters, fp);
std::vector<ClusterV2> clusters;
for (auto &c : clusters_) {
ClusterV2 cluster;
cluster.cluster = std::move(c);
cluster.frame_number = header.frame_number;
clusters.push_back(cluster);
}
return clusters;
}
std::vector<std::vector<ClusterV2>> read(int n_frames) {
std::vector<std::vector<ClusterV2>> clusters;
for (int i = 0; i < n_frames; i++) {
clusters.push_back(read());
}
return clusters;
}
void close() {
if (!m_closed) {
fclose(fp);
m_closed = true;
}
}
};
} // namespace aare

13
src/file_io/include/aare/file_io/NumpyFile.hpp
View File

@ -72,6 +72,18 @@ class NumpyFile : public FileInterface {
}
return arr;
}
template <typename A, typename TYPENAME, A Ndim> void write(NDView<TYPENAME, Ndim> &frame) {
write_impl(frame.data(), frame.total_bytes());
}
template <typename A, typename TYPENAME, A Ndim> void write(NDArray<TYPENAME, Ndim> &frame) {
write_impl(frame.data(), frame.total_bytes());
}
template <typename A, typename TYPENAME, A Ndim> void write(NDView<TYPENAME, Ndim> &&frame) {
write_impl(frame.data(), frame.total_bytes());
}
template <typename A, typename TYPENAME, A Ndim> void write(NDArray<TYPENAME, Ndim> &&frame) {
write_impl(frame.data(), frame.total_bytes());
}
~NumpyFile() noexcept override;
@ -91,6 +103,7 @@ class NumpyFile : public FileInterface {
void load_metadata();
void get_frame_into(size_t /*frame_number*/, std::byte * /*image_buf*/);
Frame get_frame(size_t frame_number);
void write_impl(void *data, uint64_t size);
};
} // namespace aare

10
src/file_io/include/aare/file_io/RawFile.hpp
View File

@ -57,15 +57,15 @@ class RawFile : public FileInterface {
*/
size_t pixels_per_frame() override { return m_rows * m_cols; }
// goto frame number
void seek(size_t frame_number) override {
// goto frame index
void seek(size_t frame_index) override {
// check if the frame number is greater than the total frames
// if frame_number == total_frames, then the next read will throw an error
if (frame_number > this->total_frames()) {
if (frame_index > this->total_frames()) {
throw std::runtime_error(
fmt::format("frame number {} is greater than total frames {}", frame_number, m_total_frames));
fmt::format("frame number {} is greater than total frames {}", frame_index, m_total_frames));
}
this->current_frame = frame_number;
this->current_frame = frame_index;
};
// return the position of the file pointer (in number of frames)

93
src/file_io/src/ClusterFile.cpp
View File

@ -15,7 +15,7 @@ namespace aare {
* @param config Configuration for the file header.
*/
ClusterFile::ClusterFile(const std::filesystem::path &fname_, const std::string &mode_, ClusterFileConfig config)
: fname{fname_}, mode{mode_}, frame_number{config.frame_number}, n_clusters{config.n_clusters} {
: fname{fname_}, mode{mode_}, frame_number{config.frame_number}, n_clusters{config.n_clusters}, dt{config.dt} {
// check if the file has the .clust extension
if (fname.extension() != ".clust") {
@ -30,12 +30,8 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname_, const std::string
if (not std::filesystem::is_regular_file(fname)) {
throw std::invalid_argument(fmt::format("file {} is not a regular file", fname.c_str()));
}
// check if the file size is a multiple of the cluster size
if ((std::filesystem::file_size(fname) - HEADER_BYTES) % sizeof(Cluster) != 0) {
aare::logger::warn("file", fname, "size is not a multiple of cluster size");
}
if (config != ClusterFileConfig()) {
aare::logger::warn("ignored ClusterFileConfig for read mode");
if (dt == DType(DType::TypeIndex::ERROR)) {
throw std::invalid_argument("data type not set in ClusterFileConfig");
}
// open file
fp = fopen(fname.c_str(), "rb");
@ -43,12 +39,13 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname_, const std::string
throw std::runtime_error(fmt::format("could not open file {}", fname.c_str()));
}
// read header
const size_t rc = fread(&config, sizeof(config), 1, fp);
ClusterFileHeader header{};
const size_t rc = fread(&header, sizeof(header), 1, fp);
if (rc != 1) {
throw std::runtime_error(fmt::format("could not read header from file {}", fname.c_str()));
}
frame_number = config.frame_number;
n_clusters = config.n_clusters;
frame_number = header.frame_number;
n_clusters = header.n_clusters;
} else if (mode == "w") {
// open file
fp = fopen(fname.c_str(), "wb");
@ -57,36 +54,56 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname_, const std::string
}
// write header
if (fwrite(&config, sizeof(config), 1, fp) != 1) {
ClusterFileHeader header{config.frame_number, config.n_clusters};
if (fwrite(&header, sizeof(header), 1, fp) != 1) {
throw std::runtime_error(fmt::format("could not write header to file {}", fname.c_str()));
}
} else {
throw std::invalid_argument("mode must be 'r' or 'w'");
}
m_cluster_size =
2 * sizeof(int16_t) + static_cast<size_t>(config.cluster_sizeX * config.cluster_sizeY * dt.bytes());
}
/**
* @brief Writes a vector of clusters to the file.
*
* Each cluster is written as a binary block of size sizeof(Cluster).
*
* @param clusters The vector of clusters to write to the file.
*/
void ClusterFile::write(std::vector<Cluster> &clusters) {
fwrite(clusters.data(), sizeof(Cluster), clusters.size(), fp);
if (clusters.empty()) {
return;
}
assert(clusters[0].dt == dt && "cluster data type mismatch");
// prepare buffer to write to file
auto *buffer = new std::byte[m_cluster_size * clusters.size()];
for (size_t i = 0; i < clusters.size(); i++) {
auto &cluster = clusters[i];
memcpy(buffer + i * m_cluster_size, &cluster.x, sizeof(int16_t));
memcpy(buffer + i * m_cluster_size + sizeof(int16_t), &cluster.y, sizeof(int16_t));
memcpy(buffer + i * m_cluster_size + 2 * sizeof(int16_t), cluster.data(), cluster.size());
}
fwrite(buffer, m_cluster_size * clusters.size(), 1, fp);
delete[] buffer;
}
/**
* @brief Writes a single cluster to the file.
*
* The cluster is written as a binary block of size sizeof(Cluster).
*
* @param cluster The cluster to write to the file.
*/
void ClusterFile::write(Cluster &cluster) { fwrite(&cluster, sizeof(Cluster), 1, fp); }
void ClusterFile::write(Cluster &cluster) {
// prepare buffer to write to file
auto *buffer = new std::byte[m_cluster_size];
memcpy(buffer, &cluster.x, sizeof(int16_t));
memcpy(buffer + sizeof(int16_t), &cluster.y, sizeof(int16_t));
memcpy(buffer + 2 * sizeof(int16_t), cluster.data(), cluster.size());
fwrite(buffer, m_cluster_size, 1, fp);
delete[] buffer;
}
/**
* @brief Reads a single cluster from the file.
*
* The cluster is read as a binary block of size sizeof(Cluster).
*
* @return The cluster read from the file.
*/
Cluster ClusterFile::read() {
@ -94,15 +111,20 @@ Cluster ClusterFile::read() {
throw std::runtime_error("cluster number out of range");
}
Cluster cluster{};
fread(&cluster, sizeof(Cluster), 1, fp);
Cluster cluster(3, 3, DType::INT32);
auto *tmp = new std::byte[cluster.size() + 2 * sizeof(int16_t)];
fread(tmp, cluster.size() + 2 * sizeof(int16_t), 1, fp);
memcpy(&cluster.x, tmp, sizeof(int16_t));
memcpy(&cluster.y, tmp + sizeof(int16_t), sizeof(int16_t));
memcpy(cluster.data(), tmp + 2 * sizeof(int16_t), cluster.size());
delete[] tmp;
return cluster;
}
/**
* @brief Reads a specific cluster from the file.
*
* The file pointer is moved to the specific cluster, and the cluster is read as a binary block of size sizeof(Cluster).
* The file pointer is moved to the specific cluster
*
* @param cluster_number The number of the cluster to read from the file.
* @return The cluster read from the file.
@ -114,8 +136,7 @@ Cluster ClusterFile::iread(size_t cluster_number) {
auto old_pos = ftell(fp);
this->seek(cluster_number);
Cluster cluster{};
fread(&cluster, sizeof(Cluster), 1, fp);
Cluster cluster = read();
fseek(fp, old_pos, SEEK_SET); // restore the file position
return cluster;
}
@ -132,8 +153,12 @@ std::vector<Cluster> ClusterFile::read(size_t n_clusters_) {
if (n_clusters_ + tell() > count()) {
throw std::runtime_error("cluster number out of range");
}
std::vector<Cluster> clusters(n_clusters_);
fread(clusters.data(), sizeof(Cluster), n_clusters, fp);
std::vector<Cluster> clusters(n_clusters_, Cluster(3, 3, DType::INT32));
// TODO: read all clusters at once if possible
for (size_t i = 0; i < n_clusters_; i++) {
clusters[i] = read();
}
return clusters;
}
@ -148,21 +173,19 @@ void ClusterFile::seek(size_t cluster_number) {
if (cluster_number > count()) {
throw std::runtime_error("cluster number out of range");
}
const auto offset = static_cast<int64_t>(sizeof(ClusterFileConfig) + cluster_number * sizeof(Cluster));
const auto offset = static_cast<int64_t>(sizeof(ClusterFileHeader) + cluster_number * m_cluster_size);
fseek(fp, offset, SEEK_SET);
}
/**
* @brief Gets the current position of the file pointer in terms of clusters.
*
* The position is calculated as the number of clusters from the beginning of the file to the current position of the
* file pointer.
* The position is calculated as the number of clusters from the beginning of the file to the current position of
* the file pointer.
*
* @return The current position of the file pointer in terms of clusters.
*/
size_t ClusterFile::tell() const { return ftell(fp) / sizeof(Cluster); }
size_t ClusterFile::tell() const { return (ftell(fp) - sizeof(ClusterFileHeader)) / m_cluster_size; }
/**
* @brief Counts the number of clusters in the file.
@ -178,7 +201,7 @@ size_t ClusterFile::count() noexcept {
// save the current position
auto old_pos = ftell(fp);
fseek(fp, 0, SEEK_END);
const size_t n_clusters_ = ftell(fp) / sizeof(Cluster);
const size_t n_clusters_ = (ftell(fp) - sizeof(ClusterFileHeader)) / m_cluster_size;
// restore the file position
fseek(fp, old_pos, SEEK_SET);
return n_clusters_;
@ -193,8 +216,8 @@ void ClusterFile::update_header() {
aare::logger::debug("updating header with correct number of clusters", count());
auto tmp_n_clusters = count();
fseek(fp, 0, SEEK_SET);
ClusterFileConfig config(frame_number, static_cast<int32_t>(tmp_n_clusters));
if (fwrite(&config, sizeof(config), 1, fp) != 1) {
ClusterFileHeader header{frame_number, static_cast<int32_t>(tmp_n_clusters)};
if (fwrite(&header, sizeof(ClusterFileHeader), 1, fp) != 1) {
throw std::runtime_error("could not write header to file");
}
if (fflush(fp) != 0) {

18
src/file_io/src/ClusterFileV2.cpp Normal file
View File

@ -0,0 +1,18 @@
#include "aare/file_io/ClusterFileV2.hpp"
namespace aare {
// ClusterFileV2::ClusterFileV2(std::filesystem::path const &fpath, std::string const &mode,
// ClusterFileConfig const &config)
// : m_fpath(fpath), m_mode(mode) {
// if (mode == "w") {
// fp = fopen(fpath.c_str(), "wb");
// write_header(config);
// } else if (mode == "r") {
// fp = fopen(fpath.c_str(), "rb");
// read_header();
// } else {
// throw std::runtime_error("invalid mode");
// }
// }
} // namespace aare

7
src/file_io/src/NumpyFile.cpp
View File

@ -30,8 +30,9 @@ NumpyFile::NumpyFile(const std::filesystem::path &fname, const std::string &mode
m_bytes_per_frame = m_header.dtype.bitdepth() / 8 * m_pixels_per_frame;
}
void NumpyFile::write(Frame &frame) { write_impl(frame.data(), frame.size()); }
void NumpyFile::write_impl(void *data, uint64_t size) {
void NumpyFile::write(Frame &frame) {
if (fp == nullptr) {
throw std::runtime_error("File not open");
}
@ -40,10 +41,12 @@ void NumpyFile::write(Frame &frame) {
}
if (fseek(fp, 0, SEEK_END))
throw std::runtime_error("Could not seek to end of file");
size_t const rc = fwrite(frame.data(), frame.size(), 1, fp);
size_t const rc = fwrite(data, size, 1, fp);
if (rc != 1) {
throw std::runtime_error("Error writing frame to file");
}
m_header.shape[0]++;
}
Frame NumpyFile::get_frame(size_t frame_number) {

4
src/file_io/src/NumpyHelpers.cpp
View File

@ -29,7 +29,7 @@ namespace aare {
std::string NumpyHeader::to_string() const {
std::stringstream sstm;
sstm << "dtype: " << dtype.str() << ", fortran_order: " << fortran_order << ' ';
sstm << "dtype: " << dtype.to_string() << ", fortran_order: " << fortran_order << ' ';
sstm << "shape: (";
for (auto item : shape)
sstm << item << ',';
@ -227,7 +227,7 @@ size_t write_header(const std::filesystem::path &fname, const NumpyHeader &heade
}
size_t write_header(std::ostream &out, const NumpyHeader &header) {
std::string const header_dict = write_header_dict(header.dtype.str(), header.fortran_order, header.shape);
std::string const header_dict = write_header_dict(header.dtype.to_string(), header.fortran_order, header.shape);
size_t length = magic_string_length + 2 + 2 + header_dict.length() + 1;

18
src/file_io/test/ClusterFile.test.cpp
View File

@ -22,7 +22,7 @@ TEST_CASE("Read a cluster file") {
REQUIRE(c.x == 1);
REQUIRE(c.y == 200);
for (int i = 0; i < 9; i++) {
REQUIRE(c.data[i] == i);
REQUIRE(c.get<int32_t>(i) == i);
}
}
SECTION("Read a single cluster using iread") {
@ -30,7 +30,7 @@ TEST_CASE("Read a cluster file") {
REQUIRE(c.x == 1);
REQUIRE(c.y == 200);
for (int i = 0; i < 9; i++) {
REQUIRE(c.data[i] == i);
REQUIRE(c.get<int32_t>(i) == i);
}
}
SECTION("Read a cluster using seek") {
@ -39,13 +39,13 @@ TEST_CASE("Read a cluster file") {
REQUIRE(c.x == 1);
REQUIRE(c.y == 200);
for (int i = 0; i < 9; i++) {
REQUIRE(c.data[i] == i);
REQUIRE(c.get<int32_t>(i) == i);
}
c = cf.read();
REQUIRE(c.x == 2);
REQUIRE(c.y == 201);
for (int i = 0; i < 9; i++) {
REQUIRE(c.data[i] == i + 9);
REQUIRE(c.get<int32_t>(i) == i + 9);
}
}
SECTION("check out of bound reading") {
@ -66,7 +66,7 @@ TEST_CASE("Read a cluster file") {
REQUIRE(c.x == offset + 1);
REQUIRE(c.y == offset + 200);
for (int i = 0; i < 9; i++) {
REQUIRE(c.data[i] == data_offset + i);
REQUIRE(c.get<int32_t>(i) == data_offset + i);
}
offset++;
@ -90,13 +90,13 @@ TEST_CASE("write a cluster file") {
std::vector<Cluster> clusters(TOTAL_CLUSTERS);
for (int32_t i = 0; i < TOTAL_CLUSTERS; i++) {
Cluster c;
c.x = INT16_MAX - offset;
c.y = INT16_MAX - (offset + 200);
c.x = (INT16_MAX - offset);
c.y = (INT16_MAX - (offset + 200));
for (int32_t j = 0; j < 9; j++) {
if (j % 2 == 0)
c.data[j] = -(offset * 2);
c.set<int32_t>(j, -(offset * 2));
else
c.data[j] = (offset * 2);
c.set<int32_t>(j, (offset * 2));
}
clusters[i] = c;
offset++;

2
src/network_io/include/aare/network_io/ZmqHeader.hpp
View File

@ -36,7 +36,7 @@ namespace simdjson {
* adds a check for 32bit overflow
*/
template <> simdjson_inline simdjson::simdjson_result<uint32_t> simdjson::ondemand::value::get() noexcept {
size_t val = 0;
uint64_t val = 0;
auto error = get_uint64().get(val);
if (error) {
return error;

2
src/network_io/include/aare/network_io/ZmqSocket.hpp
View File

@ -6,7 +6,7 @@
// needs to be in sync with the main library (or maybe better use the versioning in the header)
// forward declare zmq_msg_t to avoid including zmq.h in the header
class zmq_msg_t;
struct zmq_msg_t;
namespace aare {

2
src/network_io/include/aare/network_io/ZmqSocketReceiver.hpp
View File

@ -9,7 +9,7 @@
#include <string>
// forward declare zmq_msg_t to avoid including zmq.h in the header
class zmq_msg_t;
struct zmq_msg_t;
namespace aare {

3
src/processing/CMakeLists.txt
View File

@ -10,7 +10,8 @@ endif()
if(AARE_TESTS)
set(TestSources
${CMAKE_CURRENT_SOURCE_DIR}/test/pedestal.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/test/Pedestal.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/test/ClusterFinder.test.cpp
)
target_sources(tests PRIVATE ${TestSources} )
target_link_libraries(tests PRIVATE processing)

197
src/processing/include/aare/processing/ClusterFinder.hpp Normal file
View File

@ -0,0 +1,197 @@
#pragma once
#include "aare/core/DType.hpp"
#include "aare/core/NDArray.hpp"
#include "aare/core/NDView.hpp"
#include "aare/processing/Pedestal.hpp"
#include <cstddef>
namespace aare {
/** enum to define the event types */
enum eventType {
PEDESTAL, /** pedestal */
NEIGHBOUR, /** neighbour i.e. below threshold, but in the cluster of a photon */
PHOTON, /** photon i.e. above threshold */
PHOTON_MAX, /** maximum of a cluster satisfying the photon conditions */
NEGATIVE_PEDESTAL, /** negative value, will not be accounted for as pedestal in order to
avoid drift of the pedestal towards negative values */
UNDEFINED_EVENT = -1 /** undefined */
};
class ClusterFinder {
public:
ClusterFinder(int cluster_sizeX, int cluster_sizeY, double nSigma = 5.0, double threshold = 0.0)
: m_cluster_sizeX(cluster_sizeX), m_cluster_sizeY(cluster_sizeY), m_threshold(threshold), m_nSigma(nSigma) {
c2 = sqrt((cluster_sizeY + 1) / 2 * (cluster_sizeX + 1) / 2);
c3 = sqrt(cluster_sizeX * cluster_sizeY);
};
template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
std::vector<Cluster> find_clusters(NDView<FRAME_TYPE, 2> frame, Pedestal<PEDESTAL_TYPE> &pedestal) {
std::vector<Cluster> clusters;
std::vector<std::vector<eventType>> eventMask;
for (int i = 0; i < frame.shape(0); i++) {
eventMask.push_back(std::vector<eventType>(frame.shape(1)));
}
long double val;
long double max;
for (int iy = 0; iy < frame.shape(0); iy++) {
for (int ix = 0; ix < frame.shape(1); ix++) {
// initialize max and total
max = std::numeric_limits<FRAME_TYPE>::min();
long double total = 0;
eventMask[iy][ix] = PEDESTAL;
for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
val = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
total += val;
if (val > max) {
max = val;
}
}
}
}
auto rms = pedestal.standard_deviation(iy, ix);
if (frame(iy, ix) < -m_nSigma * rms) {
eventMask[iy][ix] = NEGATIVE_PEDESTAL;
continue;
} else if (max > m_nSigma * rms) {
eventMask[iy][ix] = PHOTON;
} else if (total > c3 * m_nSigma * rms) {
eventMask[iy][ix] = PHOTON;
} else {
pedestal.push(iy, ix, frame(iy, ix));
continue;
}
if (eventMask[iy][ix] == PHOTON and frame(iy, ix) - pedestal.mean(iy, ix) >= max) {
eventMask[iy][ix] = PHOTON_MAX;
Cluster cluster(m_cluster_sizeX, m_cluster_sizeY, DType(typeid(FRAME_TYPE)));
cluster.x = ix;
cluster.y = iy;
short i = 0;
for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
FRAME_TYPE tmp = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
cluster.set<FRAME_TYPE>(i, tmp);
i++;
}
}
}
clusters.push_back(cluster);
}
}
}
return clusters;
}
template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
std::vector<Cluster> find_clusters_with_threshold(NDView<FRAME_TYPE, 2> frame, Pedestal<PEDESTAL_TYPE> &pedestal) {
assert(m_threshold > 0);
std::vector<Cluster> clusters;
std::vector<std::vector<eventType>> eventMask;
for (int i = 0; i < frame.shape(0); i++) {
eventMask.push_back(std::vector<eventType>(frame.shape(1)));
}
double tthr, tthr1, tthr2;
NDArray<FRAME_TYPE, 2> rest({frame.shape(0), frame.shape(1)});
NDArray<int, 2> nph({frame.shape(0), frame.shape(1)});
// convert to n photons
// nph = (frame-pedestal.mean()+0.5*m_threshold)/m_threshold; // can be optimized with expression templates?
for (int iy = 0; iy < frame.shape(0); iy++) {
for (int ix = 0; ix < frame.shape(1); ix++) {
auto val = frame(iy, ix) - pedestal.mean(iy, ix);
nph(iy, ix) = (val + 0.5 * m_threshold) / m_threshold;
nph(iy, ix) = nph(iy, ix) < 0 ? 0 : nph(iy, ix);
rest(iy, ix) = val - nph(iy, ix) * m_threshold;
}
}
// iterate over frame pixels
for (int iy = 0; iy < frame.shape(0); iy++) {
for (int ix = 0; ix < frame.shape(1); ix++) {
eventMask[iy][ix] = PEDESTAL;
// initialize max and total
FRAME_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
long double total = 0;
if (rest(iy, ix) <= 0.25 * m_threshold) {
pedestal.push(iy, ix, frame(iy, ix));
continue;
}
eventMask[iy][ix] = NEIGHBOUR;
// iterate over cluster pixels around the current pixel (ix,iy)
for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
auto val = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
total += val;
if (val > max) {
max = val;
}
}
}
}
auto rms = pedestal.standard_deviation(iy, ix);
if (m_nSigma == 0) {
tthr = m_threshold;
tthr1 = m_threshold;
tthr2 = m_threshold;
} else {
tthr = m_nSigma * rms;
tthr1 = m_nSigma * rms * c3;
tthr2 = m_nSigma * rms * c2;
if (m_threshold > 2 * tthr)
tthr = m_threshold - tthr;
if (m_threshold > 2 * tthr1)
tthr1 = tthr - tthr1;
if (m_threshold > 2 * tthr2)
tthr2 = tthr - tthr2;
}
if (total > tthr1 || max > tthr) {
eventMask[iy][ix] = PHOTON;
nph(iy, ix) += 1;
rest(iy, ix) -= m_threshold;
} else {
pedestal.push(iy, ix, frame(iy, ix));
continue;
}
if (eventMask[iy][ix] == PHOTON and frame(iy, ix) - pedestal.mean(iy, ix) == max) {
eventMask[iy][ix] = PHOTON_MAX;
Cluster cluster(m_cluster_sizeX, m_cluster_sizeY, DType(typeid(FRAME_TYPE)));
cluster.x = ix;
cluster.y = iy;
short i = 0;
for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
auto tmp = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
cluster.set<FRAME_TYPE>(i, tmp);
i++;
}
}
}
clusters.push_back(cluster);
}
}
}
return clusters;
}
protected:
int m_cluster_sizeX;
int m_cluster_sizeY;
double m_threshold;
double m_nSigma;
double c2;
double c3;
};
} // namespace aare

57
src/processing/test/ClusterFinder.test.cpp Normal file
View File

@ -0,0 +1,57 @@
#include "aare/processing/ClusterFinder.hpp"
#include "aare/processing/Pedestal.hpp"
#include "aare/utils/floats.hpp"
#include <catch2/catch_test_macros.hpp>
#include <chrono>
#include <random>
using namespace aare;
class ClusterFinderUnitTest : public ClusterFinder {
public:
ClusterFinderUnitTest(int cluster_sizeX, int cluster_sizeY, double nSigma = 5.0, double threshold = 0.0)
: ClusterFinder(cluster_sizeX, cluster_sizeY, nSigma, threshold) {}
double get_c2() { return c2; }
double get_c3() { return c3; }
auto get_threshold() { return m_threshold; }
auto get_nSigma() { return m_nSigma; }
auto get_cluster_sizeX() { return m_cluster_sizeX; }
auto get_cluster_sizeY() { return m_cluster_sizeY; }
};
TEST_CASE("test ClusterFinder constructor") {
ClusterFinderUnitTest cf(55, 100);
REQUIRE(cf.get_cluster_sizeX() == 55);
REQUIRE(cf.get_cluster_sizeY() == 100);
REQUIRE(cf.get_threshold() == 0.0);
REQUIRE(cf.get_nSigma() == 5.0);
double c2 = sqrt((100 + 1) / 2 * (55 + 1) / 2);
double c3 = sqrt(55 * 100);
REQUIRE(compare_floats<double>(cf.get_c2(), c2));
REQUIRE(compare_floats<double>(cf.get_c3(), c3));
}
TEST_CASE("test cluster finder") {
aare::Pedestal pedestal(10, 10, 5);
NDArray<double, 2> frame({10, 10});
frame = 0;
ClusterFinder clusterFinder(3, 3, 1, 1); // 3x3 cluster, 1 nSigma, 1 threshold
auto clusters = clusterFinder.find_clusters(frame.span(), pedestal);
REQUIRE(clusters.size() == 0);
frame(5, 5) = 10;
clusters = clusterFinder.find_clusters(frame.span(), pedestal);
REQUIRE(clusters.size() == 1);
REQUIRE(clusters[0].x == 5);
REQUIRE(clusters[0].y == 5);
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
if (i == 1 && j == 1)
REQUIRE(clusters[0].get<double>(i * 3 + j) == 10);
else
REQUIRE(clusters[0].get<double>(i * 3 + j) == 0);
}
}
}

0
src/processing/test/pedestal.test.cpp → src/processing/test/Pedestal.test.cpp
View File