diff --git a/CMakeLists.txt b/CMakeLists.txt index a38e8fd..24b7b30 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -40,7 +40,7 @@ option(AARE_DOCS "Build documentation" OFF) option(AARE_VERBOSE "Verbose output" OFF) option(AARE_CUSTOM_ASSERT "Use custom assert" OFF) option(AARE_INSTALL_PYTHONEXT "Install the python extension in the install tree under CMAKE_INSTALL_PREFIX/aare/" OFF) - +option(AARE_ASAN "Enable AddressSanitizer" OFF) # Configure which of the dependencies to use FetchContent for option(AARE_FETCH_FMT "Use FetchContent to download fmt" ON) @@ -225,13 +225,6 @@ if(CMAKE_BUILD_TYPE STREQUAL "Release") target_compile_options(aare_compiler_flags INTERFACE -O3) else() message(STATUS "Debug build") - target_compile_options( - aare_compiler_flags - INTERFACE - -Og - -ggdb3 - ) - endif() # Common flags for GCC and Clang @@ -256,7 +249,21 @@ target_compile_options( endif() #GCC/Clang specific - +if(AARE_ASAN) + message(STATUS "AddressSanitizer enabled") + target_compile_options( + aare_compiler_flags + INTERFACE + -fsanitize=address,undefined,pointer-compare + -fno-omit-frame-pointer + ) + target_link_libraries( + aare_compiler_flags + INTERFACE + -fsanitize=address,undefined,pointer-compare + -fno-omit-frame-pointer + ) +endif() @@ -316,6 +323,8 @@ target_include_directories(aare_core PUBLIC "$" ) + + target_link_libraries( aare_core PUBLIC diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index 7111cf9..a98114d 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -82,8 +82,8 @@ class ClusterFinder { // // TODO! deal with even size clusters // // currently 3,3 -> +/- 1 // // 4,4 -> +/- 2 - short dy = m_cluster_sizeY / 2; - short dx = m_cluster_sizeX / 2; + int dy = m_cluster_sizeY / 2; + int dx = m_cluster_sizeX / 2; std::vector cluster_data(m_cluster_sizeX * m_cluster_sizeY); for (int iy = 0; iy < frame.shape(0); iy++) { @@ -100,8 +100,8 @@ class ClusterFinder { continue; // NEGATIVE_PEDESTAL go to next pixel // TODO! No pedestal update??? - for (short ir = -dy; ir < dy + 1; ir++) { - for (short ic = -dx; ic < dx + 1; ic++) { + for (int ir = -dy; ir < dy + 1; ir++) { + for (int ic = -dx; ic < dx + 1; ic++) { if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) { PEDESTAL_TYPE val = diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index 53aeed7..ce8d935 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -13,12 +13,12 @@ namespace aare { * contiguous memory buffer to store the clusters. * @note push_back can invalidate pointers to elements in the container * @tparam T data type of the pixels in the cluster + * @tparam CoordType data type of the x and y coordinates of the cluster (normally int16_t) */ -template class ClusterVector { +template class ClusterVector { using value_type = T; - using coord_t = int16_t; - coord_t m_cluster_size_x; - coord_t m_cluster_size_y; + size_t m_cluster_size_x; + size_t m_cluster_size_y; std::byte *m_data{}; size_t m_size{0}; size_t m_capacity; @@ -39,7 +39,7 @@ template class ClusterVector { * @param cluster_size_y size of the cluster in y direction * @param capacity initial capacity of the buffer in number of clusters */ - ClusterVector(coord_t cluster_size_x, coord_t cluster_size_y, + ClusterVector(size_t cluster_size_x, size_t cluster_size_y, size_t capacity = 1024) : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y), m_capacity(capacity) { @@ -94,21 +94,22 @@ template class ClusterVector { * @param data pointer to the data of the cluster * @warning The data pointer must point to a buffer of size cluster_size_x * cluster_size_y * sizeof(T) */ - void push_back(coord_t x, coord_t y, const std::byte *data) { + void push_back(CoordType x, CoordType y, const std::byte *data) { if (m_size == m_capacity) { allocate_buffer(m_capacity * 2); } std::byte *ptr = element_ptr(m_size); - *reinterpret_cast(ptr) = x; - ptr += sizeof(coord_t); - *reinterpret_cast(ptr) = y; - ptr += sizeof(coord_t); + *reinterpret_cast(ptr) = x; + ptr += sizeof(CoordType); + *reinterpret_cast(ptr) = y; + ptr += sizeof(CoordType); std::copy(data, data + m_cluster_size_x * m_cluster_size_y * sizeof(T), ptr); m_size++; } + /** * @brief Sum the pixels in each cluster * @return std::vector vector of sums for each cluster @@ -117,7 +118,7 @@ template class ClusterVector { std::vector sums(m_size); const size_t stride = element_offset(); const size_t n_pixels = m_cluster_size_x * m_cluster_size_y; - std::byte *ptr = m_data + 2 * sizeof(coord_t); // skip x and y + std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y for (size_t i = 0; i < m_size; i++) { sums[i] = @@ -135,7 +136,7 @@ template class ClusterVector { * @brief Return the offset in bytes for a single cluster */ size_t element_offset() const { - return sizeof(m_cluster_size_x) + sizeof(m_cluster_size_y) + + return 2*sizeof(CoordType) + m_cluster_size_x * m_cluster_size_y * sizeof(T); } /** @@ -148,8 +149,8 @@ template class ClusterVector { */ std::byte *element_ptr(size_t i) { return m_data + element_offset(i); } - int16_t cluster_size_x() const { return m_cluster_size_x; } - int16_t cluster_size_y() const { return m_cluster_size_y; } + size_t cluster_size_x() const { return m_cluster_size_x; } + size_t cluster_size_y() const { return m_cluster_size_y; } std::byte *data() { return m_data; } const std::byte *data() const { return m_data; } diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp index 346646c..15beb02 100644 --- a/include/aare/NDArray.hpp +++ b/include/aare/NDArray.hpp @@ -87,7 +87,7 @@ class NDArray : public ArrayExpr, Ndim> { // Conversion operator from array expression to array template NDArray(ArrayExpr &&expr) : NDArray(expr.shape()) { - for (int i = 0; i < size_; ++i) { + for (size_t i = 0; i < size_; ++i) { data_[i] = expr[i]; } } @@ -159,11 +159,11 @@ class NDArray : public ArrayExpr, Ndim> { } // TODO! is int the right type for index? - T &operator()(int i) { return data_[i]; } - const T &operator()(int i) const { return data_[i]; } + T &operator()(int64_t i) { return data_[i]; } + const T &operator()(int64_t i) const { return data_[i]; } - T &operator[](int i) { return data_[i]; } - const T &operator[](int i) const { return data_[i]; } + T &operator[](int64_t i) { return data_[i]; } + const T &operator[](int64_t i) const { return data_[i]; } T *data() { return data_; } std::byte *buffer() { return reinterpret_cast(data_); } diff --git a/include/aare/Pedestal.hpp b/include/aare/Pedestal.hpp index 216c204..bb2ea2c 100644 --- a/include/aare/Pedestal.hpp +++ b/include/aare/Pedestal.hpp @@ -23,31 +23,43 @@ template class Pedestal { NDArray m_sum; NDArray m_sum2; + //Cache mean since it is used over and over in the ClusterFinder + //This optimization is related to the access pattern of the ClusterFinder + //Relies on having more reads than pushes to the pedestal + NDArray m_mean; + public: Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000) : m_rows(rows), m_cols(cols), m_samples(n_samples), m_cur_samples(NDArray({rows, cols}, 0)), m_sum(NDArray({rows, cols})), - m_sum2(NDArray({rows, cols})) { + m_sum2(NDArray({rows, cols})), + m_mean(NDArray({rows, cols})) { assert(rows > 0 && cols > 0 && n_samples > 0); m_sum = 0; m_sum2 = 0; + m_mean = 0; } ~Pedestal() = default; NDArray mean() { - NDArray mean_array({m_rows, m_cols}); - for (uint32_t i = 0; i < m_rows * m_cols; i++) { - mean_array(i / m_cols, i % m_cols) = mean(i / m_cols, i % m_cols); - } - return mean_array; + return m_mean; } SUM_TYPE mean(const uint32_t row, const uint32_t col) const { + return m_mean(row, col); + } + + SUM_TYPE std(const uint32_t row, const uint32_t col) const { + return std::sqrt(variance(row, col)); + } + + SUM_TYPE variance(const uint32_t row, const uint32_t col) const { if (m_cur_samples(row, col) == 0) { return 0.0; } - return m_sum(row, col) / m_cur_samples(row, col); + return m_sum2(row, col) / m_cur_samples(row, col) - + mean(row, col) * mean(row, col); } NDArray variance() { @@ -59,13 +71,7 @@ template class Pedestal { return variance_array; } - SUM_TYPE variance(const uint32_t row, const uint32_t col) const { - if (m_cur_samples(row, col) == 0) { - return 0.0; - } - return m_sum2(row, col) / m_cur_samples(row, col) - - mean(row, col) * mean(row, col); - } + NDArray std() { NDArray standard_deviation_array({m_rows, m_cols}); @@ -77,14 +83,12 @@ template class Pedestal { return standard_deviation_array; } - SUM_TYPE std(const uint32_t row, const uint32_t col) const { - return std::sqrt(variance(row, col)); - } + void clear() { - for (uint32_t i = 0; i < m_rows * m_cols; i++) { - clear(i / m_cols, i % m_cols); - } + m_sum = 0; + m_sum2 = 0; + m_cur_samples = 0; } @@ -104,8 +108,8 @@ template class Pedestal { "Frame shape does not match pedestal shape"); } - for (uint32_t row = 0; row < m_rows; row++) { - for (uint32_t col = 0; col < m_cols; col++) { + for (size_t row = 0; row < m_rows; row++) { + for (size_t col = 0; col < m_cols; col++) { push(row, col, frame(row, col)); } } @@ -134,18 +138,17 @@ template class Pedestal { template void push(const uint32_t row, const uint32_t col, const T val_) { SUM_TYPE val = static_cast(val_); - const uint32_t idx = index(row, col); - if (m_cur_samples(idx) < m_samples) { - m_sum(idx) += val; - m_sum2(idx) += val * val; - m_cur_samples(idx)++; + if (m_cur_samples(row, col) < m_samples) { + m_sum(row, col) += val; + m_sum2(row, col) += val * val; + m_cur_samples(row, col)++; } else { - m_sum(idx) += val - m_sum(idx) / m_cur_samples(idx); - m_sum2(idx) += val * val - m_sum2(idx) / m_cur_samples(idx); + m_sum(row, col) += val - m_sum(row, col) / m_cur_samples(row, col); + m_sum2(row, col) += val * val - m_sum2(row, col) / m_cur_samples(row, col); } + //Since we just did a push we know that m_cur_samples(row, col) is at least 1 + m_mean(row, col) = m_sum(row, col) / m_cur_samples(row, col); } - uint32_t index(const uint32_t row, const uint32_t col) const { - return row * m_cols + col; - }; + }; } // namespace aare \ No newline at end of file diff --git a/python/aare/__init__.py b/python/aare/__init__.py index 5641d85..fb34c7a 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -3,9 +3,10 @@ from . import _aare from ._aare import File, RawMasterFile, RawSubFile -from ._aare import Pedestal, ClusterFinder, VarClusterFinder +from ._aare import Pedestal_d, Pedestal_f, ClusterFinder, VarClusterFinder from ._aare import DetectorType from ._aare import ClusterFile +from ._aare import hitmap from .CtbRawFile import CtbRawFile from .RawFile import RawFile diff --git a/python/src/cluster.hpp b/python/src/cluster.hpp index 0467c98..d11c706 100644 --- a/python/src/cluster.hpp +++ b/python/src/cluster.hpp @@ -80,7 +80,26 @@ void define_cluster_finder_bindings(py::module &m) { return; }); - + m.def("hitmap", [](std::array image_size, ClusterVector& cv){ + + py::array_t hitmap(image_size); + auto r = hitmap.mutable_unchecked<2>(); + + // Initialize hitmap to 0 + for (py::ssize_t i = 0; i < r.shape(0); i++) + for (py::ssize_t j = 0; j < r.shape(1); j++) + r(i, j) = 0; + + size_t stride = cv.element_offset(); + auto ptr = cv.data(); + for(size_t i=0; i(ptr); + auto y = *reinterpret_cast(ptr+sizeof(int16_t)); + r(y, x) += 1; + ptr += stride; + } + return hitmap; + }); define_cluster_vector(m, "i"); define_cluster_vector(m, "d"); define_cluster_vector(m, "f"); diff --git a/python/src/cluster_file.hpp b/python/src/cluster_file.hpp index 543073f..aa7fd23 100644 --- a/python/src/cluster_file.hpp +++ b/python/src/cluster_file.hpp @@ -10,6 +10,12 @@ #include #include +//Disable warnings for unused parameters, as we ignore some +//in the __exit__ method +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-parameter" + + namespace py = pybind11; using namespace ::aare; @@ -60,4 +66,6 @@ void define_cluster_file_io_bindings(py::module &m) { } return return_vector(vec); }); -} \ No newline at end of file +} + +#pragma GCC diagnostic pop \ No newline at end of file diff --git a/python/src/module.cpp b/python/src/module.cpp index 7963ac4..14a686a 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -22,8 +22,8 @@ PYBIND11_MODULE(_aare, m) { define_raw_master_file_bindings(m); define_var_cluster_finder_bindings(m); define_pixel_map_bindings(m); - define_pedestal_bindings(m, "Pedestal"); - define_pedestal_bindings(m, "Pedestal_float32"); + define_pedestal_bindings(m, "Pedestal_d"); + define_pedestal_bindings(m, "Pedestal_f"); define_cluster_finder_bindings(m); define_cluster_file_io_bindings(m); } \ No newline at end of file diff --git a/src/ClusterVector.test.cpp b/src/ClusterVector.test.cpp index ef4e0ee..24a482b 100644 --- a/src/ClusterVector.test.cpp +++ b/src/ClusterVector.test.cpp @@ -1,7 +1,7 @@ #include #include "aare/ClusterVector.hpp" -// #include +#include #include using aare::ClusterVector; @@ -44,10 +44,10 @@ TEST_CASE("Summing 3x1 clusters of int64"){ struct Cluster_l3x1{ int16_t x; int16_t y; - int64_t data[3]; + int32_t data[3]; }; - ClusterVector cv(3, 1, 2); + ClusterVector cv(3, 1, 2); REQUIRE(cv.capacity() == 2); REQUIRE(cv.size() == 0); REQUIRE(cv.cluster_size_x() == 3); @@ -74,4 +74,35 @@ TEST_CASE("Summing 3x1 clusters of int64"){ REQUIRE(sums[0] == 12); REQUIRE(sums[1] == 27); REQUIRE(sums[2] == 42); +} + +TEST_CASE("Storing floats"){ + struct Cluster_f4x2{ + int16_t x; + int16_t y; + float data[8]; + }; + + ClusterVector cv(2, 4, 2); + REQUIRE(cv.capacity() == 2); + REQUIRE(cv.size() == 0); + REQUIRE(cv.cluster_size_x() == 2); + REQUIRE(cv.cluster_size_y() == 4); + + //Create a cluster and push back into the vector + Cluster_f4x2 c1 = {1, 2, {3.0, 4.0, 5.0, 6.0,3.0, 4.0, 5.0, 6.0}}; + cv.push_back(c1.x, c1.y, reinterpret_cast(&c1.data[0])); + REQUIRE(cv.capacity() == 2); + REQUIRE(cv.size() == 1); + + + Cluster_f4x2 c2 = {6, 7, {8.0, 9.0, 10.0, 11.0,8.0, 9.0, 10.0, 11.0}}; + cv.push_back(c2.x, c2.y, reinterpret_cast(&c2.data[0])); + REQUIRE(cv.capacity() == 2); + REQUIRE(cv.size() == 2); + + auto sums = cv.sum(); + REQUIRE(sums.size() == 2); + REQUIRE_THAT(sums[0], Catch::Matchers::WithinAbs(36.0, 1e-6)); + REQUIRE_THAT(sums[1], Catch::Matchers::WithinAbs(76.0, 1e-6)); } \ No newline at end of file diff --git a/src/Frame.test.cpp b/src/Frame.test.cpp index 33bbbb6..4063701 100644 --- a/src/Frame.test.cpp +++ b/src/Frame.test.cpp @@ -19,7 +19,7 @@ TEST_CASE("Construct a frame") { // data should be initialized to 0 for (size_t i = 0; i < rows; i++) { for (size_t j = 0; j < cols; j++) { - uint8_t *data = (uint8_t *)frame.pixel_ptr(i, j); + uint8_t *data = reinterpret_cast(frame.pixel_ptr(i, j)); REQUIRE(data != nullptr); REQUIRE(*data == 0); } @@ -40,7 +40,7 @@ TEST_CASE("Set a value in a 8 bit frame") { // only the value we did set should be non-zero for (size_t i = 0; i < rows; i++) { for (size_t j = 0; j < cols; j++) { - uint8_t *data = (uint8_t *)frame.pixel_ptr(i, j); + uint8_t *data = reinterpret_cast(frame.pixel_ptr(i, j)); REQUIRE(data != nullptr); if (i == 5 && j == 7) { REQUIRE(*data == value); @@ -65,7 +65,7 @@ TEST_CASE("Set a value in a 64 bit frame") { // only the value we did set should be non-zero for (size_t i = 0; i < rows; i++) { for (size_t j = 0; j < cols; j++) { - uint64_t *data = (uint64_t *)frame.pixel_ptr(i, j); + uint64_t *data = reinterpret_cast(frame.pixel_ptr(i, j)); REQUIRE(data != nullptr); if (i == 5 && j == 7) { REQUIRE(*data == value); @@ -149,4 +149,5 @@ TEST_CASE("test explicit copy constructor") { REQUIRE(frame2.bitdepth() == bitdepth); REQUIRE(frame2.bytes() == rows * cols * bitdepth / 8); REQUIRE(frame2.data() != data); -} \ No newline at end of file +} + diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3170f7c..1906508 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -17,8 +17,8 @@ endif() list(APPEND CMAKE_MODULE_PATH ${Catch2_SOURCE_DIR}/extras) add_executable(tests test.cpp) -target_link_libraries(tests PRIVATE Catch2::Catch2WithMain) - +target_link_libraries(tests PRIVATE Catch2::Catch2WithMain aare_core aare_compiler_flags) +# target_compile_options(tests PRIVATE -fno-omit-frame-pointer -fsanitize=address) set_target_properties(tests PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR} OUTPUT_NAME run_tests @@ -34,7 +34,7 @@ set(TestSources target_sources(tests PRIVATE ${TestSources} ) #Work around to remove, this is not the way to do it =) -target_link_libraries(tests PRIVATE aare_core aare_compiler_flags) +# target_link_libraries(tests PRIVATE aare_core aare_compiler_flags) #configure a header to pass test file paths diff --git a/tests/test.cpp b/tests/test.cpp index 7c638e4..513f690 100644 --- a/tests/test.cpp +++ b/tests/test.cpp @@ -3,6 +3,7 @@ #include #include #include +#include TEST_CASE("Test suite can find data assets", "[.integration]") { auto fpath = test_data_path() / "numpy" / "test_numpy_file.npy"; @@ -18,4 +19,20 @@ TEST_CASE("Test suite can open data assets", "[.integration]") { TEST_CASE("Test float32 and char8") { REQUIRE(sizeof(float) == 4); REQUIRE(CHAR_BIT == 8); -} \ No newline at end of file +} + +/** + * Uncomment the following tests to verify that asan is working + */ + +// TEST_CASE("trigger asan stack"){ +// int arr[5] = {1,2,3,4,5}; +// int val = arr[7]; +// fmt::print("val: {}\n", val); +// } + +// TEST_CASE("trigger asan heap"){ +// auto *ptr = new int[5]; +// ptr[70] = 5; +// fmt::print("ptr: {}\n", ptr[70]); +// } \ No newline at end of file