diff --git a/.CMakeLists.txt.swp b/.CMakeLists.txt.swp new file mode 100644 index 0000000..e6e9f10 Binary files /dev/null and b/.CMakeLists.txt.swp differ diff --git a/CMakeLists.txt b/CMakeLists.txt index ca110b0..a3cc152 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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() diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 78b7143..5e9e46c 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -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) diff --git a/examples/cluster_example.cpp b/examples/cluster_example.cpp index f25dc06..447f2cb 100644 --- a/examples/cluster_example.cpp +++ b/examples/cluster_example.cpp @@ -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 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 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()); + } } \ No newline at end of file diff --git a/examples/cluster_finder_example.cpp b/examples/cluster_finder_example.cpp new file mode 100644 index 0000000..d9de551 --- /dev/null +++ b/examples/cluster_finder_example.cpp @@ -0,0 +1,63 @@ + +#include "aare.hpp" +#include "aare/examples/defs.hpp" + +#include +#include + +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()); + } + 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()); + + auto clusters = clusterFinder.find_clusters(frame.view(), p); + logger::info("nclusters:", clusters.size()); + + // aare::logger::info("nclusters:", clusters.size()); + + // for (auto &cluster : clusters) { + // aare::logger::info("cluster center:", cluster.to_string()); + // } +} \ No newline at end of file diff --git a/extra/.gitignore b/extra/.gitignore index 8b13789..6aac717 100644 --- a/extra/.gitignore +++ b/extra/.gitignore @@ -1 +1,6 @@ - +sqlite_clusters.db +fixed_size_clusters.bin +numpy_clusters.npy +sqlite_clusters.zip +fixed_size_clusters.zip +numpy_clusters.zip diff --git a/extra/compare_file_formats.py b/extra/compare_file_formats.py new file mode 100644 index 0000000..4d24ca6 --- /dev/null +++ b/extra/compare_file_formats.py @@ -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) diff --git a/include/aare.hpp b/include/aare.hpp index 5017d34..c0d8e80 100644 --- a/include/aare.hpp +++ b/include/aare.hpp @@ -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" diff --git a/include/aare/file_io.hpp b/include/aare/file_io.hpp index 6c132f6..350a504 100644 --- a/include/aare/file_io.hpp +++ b/include/aare/file_io.hpp @@ -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" diff --git a/include/aare/processing.hpp b/include/aare/processing.hpp index f1d1592..8b85e1a 100644 --- a/include/aare/processing.hpp +++ b/include/aare/processing.hpp @@ -1 +1,2 @@ +#include "aare/processing/ClusterFinder.hpp" #include "aare/processing/Pedestal.hpp" \ No newline at end of file diff --git a/src/core/include/aare/core/DType.hpp b/src/core/include/aare/core/DType.hpp index d5aa42c..122cd1f 100644 --- a/src/core/include/aare/core/DType.hpp +++ b/src/core/include/aare/core/DType.hpp @@ -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}; diff --git a/src/core/include/aare/core/Frame.hpp b/src/core/include/aare/core/Frame.hpp index 1f2c2a8..85c0b4c 100644 --- a/src/core/include/aare/core/Frame.hpp +++ b/src/core/include/aare/core/Frame.hpp @@ -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; diff --git a/src/core/include/aare/core/NDView.hpp b/src/core/include/aare/core/NDView.hpp index a39831f..2a54e67 100644 --- a/src/core/include/aare/core/NDView.hpp +++ b/src/core/include/aare/core/NDView.hpp @@ -67,6 +67,7 @@ template 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_; } diff --git a/src/core/include/aare/core/defs.hpp b/src/core/include/aare/core/defs.hpp index 05fa509..3681258 100644 --- a/src/core/include/aare/core/defs.hpp +++ b/src/core/include/aare/core/defs.hpp @@ -1,9 +1,14 @@ #pragma once +#include "aare/core/DType.hpp" +#include "aare/utils/logger.hpp" + #include #include +#include #include +#include #include #include #include @@ -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 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 T get(int idx) { + (sizeof(T) == dt.bytes()) ? 0 : throw std::invalid_argument("[ERROR] Type size mismatch"); + return *reinterpret_cast(m_data + idx * dt.bytes()); + } + template 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 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(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; } }; /** diff --git a/src/core/src/DType.cpp b/src/core/src/DType.cpp index 544f213..6411cbc 100644 --- a/src/core/src/DType.cpp +++ b/src/core/src/DType.cpp @@ -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) diff --git a/src/core/test/DType.test.cpp b/src/core/test/DType.test.cpp index 98ec528..0c1e961 100644 --- a/src/core/test/DType.test.cpp +++ b/src/core/test/DType.test.cpp @@ -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() == " #include 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(0, v); + double v2 = c2.get(0); + REQUIRE(aare::compare_floats(v, v2)); + + c2.set(33 * 44 - 1, 123.11); + double v3 = c2.get(33 * 44 - 1); + REQUIRE(aare::compare_floats(123.11, v3)); + + REQUIRE_THROWS_AS(c2.set(0, 1), std::invalid_argument); // set int to double } \ No newline at end of file diff --git a/src/file_io/include/aare/file_io/ClusterFile.hpp b/src/file_io/include/aare/file_io/ClusterFile.hpp index 71d4970..4296d47 100644 --- a/src/file_io/include/aare/file_io/ClusterFile.hpp +++ b/src/file_io/include/aare/file_io/ClusterFile.hpp @@ -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 \ No newline at end of file diff --git a/src/file_io/include/aare/file_io/ClusterFileV2.hpp b/src/file_io/include/aare/file_io/ClusterFileV2.hpp new file mode 100644 index 0000000..d701ad8 --- /dev/null +++ b/src/file_io/include/aare/file_io/ClusterFileV2.hpp @@ -0,0 +1,68 @@ +#pragma once +#include "aare/core/defs.hpp" +#include +#include + +namespace aare { +struct ClusterHeader { + int32_t frame_number; + int32_t n_clusters; +}; + +struct ClusterV2_ { + int16_t x; + int16_t y; + std::array 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 read() { + ClusterHeader header; + fread(&header, sizeof(ClusterHeader), 1, fp); + std::vector clusters_(header.n_clusters); + fread(clusters_.data(), sizeof(ClusterV2_), header.n_clusters, fp); + std::vector 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> read(int n_frames) { + std::vector> 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 \ No newline at end of file diff --git a/src/file_io/include/aare/file_io/NumpyFile.hpp b/src/file_io/include/aare/file_io/NumpyFile.hpp index 3288626..e11f7d8 100644 --- a/src/file_io/include/aare/file_io/NumpyFile.hpp +++ b/src/file_io/include/aare/file_io/NumpyFile.hpp @@ -72,6 +72,18 @@ class NumpyFile : public FileInterface { } return arr; } + template void write(NDView &frame) { + write_impl(frame.data(), frame.total_bytes()); + } + template void write(NDArray &frame) { + write_impl(frame.data(), frame.total_bytes()); + } + template void write(NDView &&frame) { + write_impl(frame.data(), frame.total_bytes()); + } + template void write(NDArray &&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 \ No newline at end of file diff --git a/src/file_io/include/aare/file_io/RawFile.hpp b/src/file_io/include/aare/file_io/RawFile.hpp index 1047d1a..feb65db 100644 --- a/src/file_io/include/aare/file_io/RawFile.hpp +++ b/src/file_io/include/aare/file_io/RawFile.hpp @@ -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) diff --git a/src/file_io/src/ClusterFile.cpp b/src/file_io/src/ClusterFile.cpp index 877e974..3937efe 100644 --- a/src/file_io/src/ClusterFile.cpp +++ b/src/file_io/src/ClusterFile.cpp @@ -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(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 &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 ClusterFile::read(size_t n_clusters_) { if (n_clusters_ + tell() > count()) { throw std::runtime_error("cluster number out of range"); } - std::vector clusters(n_clusters_); - fread(clusters.data(), sizeof(Cluster), n_clusters, fp); + std::vector 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(sizeof(ClusterFileConfig) + cluster_number * sizeof(Cluster)); - + const auto offset = static_cast(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(tmp_n_clusters)); - if (fwrite(&config, sizeof(config), 1, fp) != 1) { + ClusterFileHeader header{frame_number, static_cast(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) { diff --git a/src/file_io/src/ClusterFileV2.cpp b/src/file_io/src/ClusterFileV2.cpp new file mode 100644 index 0000000..83c8238 --- /dev/null +++ b/src/file_io/src/ClusterFileV2.cpp @@ -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 \ No newline at end of file diff --git a/src/file_io/src/NumpyFile.cpp b/src/file_io/src/NumpyFile.cpp index aab1ebe..9c2e716 100644 --- a/src/file_io/src/NumpyFile.cpp +++ b/src/file_io/src/NumpyFile.cpp @@ -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) { diff --git a/src/file_io/src/NumpyHelpers.cpp b/src/file_io/src/NumpyHelpers.cpp index cc154de..2a0ebb1 100644 --- a/src/file_io/src/NumpyHelpers.cpp +++ b/src/file_io/src/NumpyHelpers.cpp @@ -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; diff --git a/src/file_io/test/ClusterFile.test.cpp b/src/file_io/test/ClusterFile.test.cpp index 581b054..22eab6c 100644 --- a/src/file_io/test/ClusterFile.test.cpp +++ b/src/file_io/test/ClusterFile.test.cpp @@ -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(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(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(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(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(i) == data_offset + i); } offset++; @@ -90,13 +90,13 @@ TEST_CASE("write a cluster file") { std::vector 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(j, -(offset * 2)); else - c.data[j] = (offset * 2); + c.set(j, (offset * 2)); } clusters[i] = c; offset++; diff --git a/src/network_io/include/aare/network_io/ZmqHeader.hpp b/src/network_io/include/aare/network_io/ZmqHeader.hpp index 6562fec..7f3c17e 100644 --- a/src/network_io/include/aare/network_io/ZmqHeader.hpp +++ b/src/network_io/include/aare/network_io/ZmqHeader.hpp @@ -36,7 +36,7 @@ namespace simdjson { * adds a check for 32bit overflow */ template <> simdjson_inline simdjson::simdjson_result simdjson::ondemand::value::get() noexcept { - size_t val = 0; + uint64_t val = 0; auto error = get_uint64().get(val); if (error) { return error; diff --git a/src/network_io/include/aare/network_io/ZmqSocket.hpp b/src/network_io/include/aare/network_io/ZmqSocket.hpp index f10bfaf..e082492 100644 --- a/src/network_io/include/aare/network_io/ZmqSocket.hpp +++ b/src/network_io/include/aare/network_io/ZmqSocket.hpp @@ -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 { diff --git a/src/network_io/include/aare/network_io/ZmqSocketReceiver.hpp b/src/network_io/include/aare/network_io/ZmqSocketReceiver.hpp index e708acb..7a3a388 100644 --- a/src/network_io/include/aare/network_io/ZmqSocketReceiver.hpp +++ b/src/network_io/include/aare/network_io/ZmqSocketReceiver.hpp @@ -9,7 +9,7 @@ #include // forward declare zmq_msg_t to avoid including zmq.h in the header -class zmq_msg_t; +struct zmq_msg_t; namespace aare { diff --git a/src/processing/CMakeLists.txt b/src/processing/CMakeLists.txt index 0800ca4..689355a 100644 --- a/src/processing/CMakeLists.txt +++ b/src/processing/CMakeLists.txt @@ -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) diff --git a/src/processing/include/aare/processing/ClusterFinder.hpp b/src/processing/include/aare/processing/ClusterFinder.hpp new file mode 100644 index 0000000..63772d8 --- /dev/null +++ b/src/processing/include/aare/processing/ClusterFinder.hpp @@ -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 + +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 + std::vector find_clusters(NDView frame, Pedestal &pedestal) { + std::vector clusters; + std::vector> eventMask; + for (int i = 0; i < frame.shape(0); i++) { + eventMask.push_back(std::vector(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::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(i, tmp); + i++; + } + } + } + clusters.push_back(cluster); + } + } + } + return clusters; + } + + template + std::vector find_clusters_with_threshold(NDView frame, Pedestal &pedestal) { + assert(m_threshold > 0); + std::vector clusters; + std::vector> eventMask; + for (int i = 0; i < frame.shape(0); i++) { + eventMask.push_back(std::vector(frame.shape(1))); + } + double tthr, tthr1, tthr2; + + NDArray rest({frame.shape(0), frame.shape(1)}); + NDArray 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::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(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 \ No newline at end of file diff --git a/src/processing/test/ClusterFinder.test.cpp b/src/processing/test/ClusterFinder.test.cpp new file mode 100644 index 0000000..f6a9425 --- /dev/null +++ b/src/processing/test/ClusterFinder.test.cpp @@ -0,0 +1,57 @@ +#include "aare/processing/ClusterFinder.hpp" +#include "aare/processing/Pedestal.hpp" +#include "aare/utils/floats.hpp" +#include +#include +#include + +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(cf.get_c2(), c2)); + REQUIRE(compare_floats(cf.get_c3(), c3)); +} + +TEST_CASE("test cluster finder") { + aare::Pedestal pedestal(10, 10, 5); + NDArray 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(i * 3 + j) == 10); + else + REQUIRE(clusters[0].get(i * 3 + j) == 0); + } + } +} \ No newline at end of file diff --git a/src/processing/test/pedestal.test.cpp b/src/processing/test/Pedestal.test.cpp similarity index 100% rename from src/processing/test/pedestal.test.cpp rename to src/processing/test/Pedestal.test.cpp