diff --git a/CMakeLists.txt b/CMakeLists.txt index bff2afe..4772f0b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -386,6 +386,7 @@ endif() if(AARE_TESTS) set(TestSources + ${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp index 310d070..45d3a83 100644 --- a/include/aare/NDArray.hpp +++ b/include/aare/NDArray.hpp @@ -102,6 +102,9 @@ class NDArray : public ArrayExpr, Ndim> { auto begin() { return data_; } auto end() { return data_ + size_; } + auto begin() const { return data_; } + auto end() const { return data_ + size_; } + using value_type = T; NDArray &operator=(NDArray &&other) noexcept; // Move assign diff --git a/include/aare/algorithm.hpp b/include/aare/algorithm.hpp new file mode 100644 index 0000000..5d6dc57 --- /dev/null +++ b/include/aare/algorithm.hpp @@ -0,0 +1,55 @@ + +#pragma once +#include +#include +#include +#include + +namespace aare { +/** + * @brief Find the index of the last element smaller than val + * assume a sorted array + */ +template +size_t last_smaller(const T* first, const T* last, T val) { + for (auto iter = first+1; iter != last; ++iter) { + if (*iter > val) { + return std::distance(first, iter-1); + } + } + return std::distance(first, last-1); +} + +template +size_t last_smaller(const NDArray& arr, T val) { + return last_smaller(arr.begin(), arr.end(), val); +} + + +template +size_t nearest_index(const T* first, const T* last, T val) { + auto iter = std::min_element(first, last, + [val](T a, T b) { + return std::abs(a - val) < std::abs(b - val); + }); + return std::distance(first, iter); +} + +template +size_t nearest_index(const NDArray& arr, T val) { + return nearest_index(arr.begin(), arr.end(), val); +} + +template +size_t nearest_index(const std::vector& vec, T val) { + return nearest_index(vec.data(), vec.data()+vec.size(), val); +} + +template +size_t nearest_index(const std::array& arr, T val) { + return nearest_index(arr.data(), arr.data()+arr.size(), val); +} + + + +} // namespace aare \ No newline at end of file diff --git a/python/examples/play.py b/python/examples/play.py index 05f3c82..b2c368b 100644 --- a/python/examples/play.py +++ b/python/examples/play.py @@ -1,40 +1,77 @@ import sys sys.path.append('/home/l_msdetect/erik/aare/build') +from aare._aare import ClusterVector_i, Interpolator -#Our normal python imports -from pathlib import Path -import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1 import make_axes_locatable +import pickle import numpy as np +import matplotlib.pyplot as plt import boost_histogram as bh +import torch +import math import time -import tifffile -#Directly import what we need from aare -from aare import File, ClusterFile, hitmap -from aare._aare import calculate_eta2, ClusterFinderMT, ClusterCollector +def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2): + """ + Generate a 2D gaussian as position mx, my, with sigma=sigma. + The gaussian is placed on a 2x2 pixel matrix with resolution + res in one dimesion. + """ + x = torch.linspace(0, pixel_size*grid_size, res) + x,y = torch.meshgrid(x,x, indexing="ij") + return 1 / (2*math.pi*sigma**2) * \ + torch.exp(-((x - my)**2 / (2*sigma**2) + (y - mx)**2 / (2*sigma**2))) + +scale = 1000 #Scale factor when converting to integer +pixel_size = 25 #um +grid = 2 +resolution = 100 +sigma_um = 10 +xa = np.linspace(0,grid*pixel_size,resolution) +ticks = [0, 25, 50] + +hit = np.array((20,20)) +etahist_fname = "/home/l_msdetect/erik/tmp/test_hist.pkl" + +local_resolution = 99 +grid_size = 3 +xaxis = np.linspace(0,grid_size*pixel_size, local_resolution) +t = gaussian_2d(hit[0],hit[1], grid_size = grid_size, sigma = 10, res = local_resolution) +pixels = t.reshape(grid_size, t.shape[0] // grid_size, grid_size, t.shape[1] // grid_size).sum(axis = 3).sum(axis = 1) +pixels = pixels.numpy() +pixels = (pixels*scale).astype(np.int32) +v = ClusterVector_i(3,3) +v.push_back(1,1, pixels) + +with open(etahist_fname, "rb") as f: + hist = pickle.load(f) +eta = hist.view().copy() +etabinsx = np.array(hist.axes.edges.T[0].flat) +etabinsy = np.array(hist.axes.edges.T[1].flat) +ebins = np.array(hist.axes.edges.T[2].flat) +p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1]) -base = Path('/mnt/sls_det_storage/moench_data/tomcat_nanoscope_21042020/09_Moench_650um/') -# for f in base.glob('*'): -# print(f.name) +#Generate the hit -cluster_fname = base/'acq_interp_center_3.8Mfr_200V.clust' -flatfield_fname = base/'flatfield_center_200_d0_f000000000000_0.clust' -cluster_fname.stat().st_size/1e6/4 -image = np.zeros((400,400)) -with ClusterFile(cluster_fname, chunk_size = 1000000) as f: - for clusters in f: - test = hitmap(image.shape, clusters) - break - # image += hitmap(image.shape, clusters) - # break -print('We are back in python') -# fig, ax = plt.subplots(figsize = (7,7)) -# im = ax.imshow(image) -# im.set_clim(0,1) \ No newline at end of file + +tmp = p.interpolate(v) +print(f'tmp:{tmp}') +pos = np.array((tmp['x'], tmp['y']))*25 + + +print(pixels) +fig, ax = plt.subplots(figsize = (7,7)) +ax.pcolormesh(xaxis, xaxis, t) +ax.plot(*pos, 'o') +ax.set_xticks([0,25,50,75]) +ax.set_yticks([0,25,50,75]) +ax.set_xlim(0,75) +ax.set_ylim(0,75) +ax.grid() +print(f'{hit=}') +print(f'{pos=}') \ No newline at end of file diff --git a/src/Interpolator.cpp b/src/Interpolator.cpp index 0e72849..85e0b5d 100644 --- a/src/Interpolator.cpp +++ b/src/Interpolator.cpp @@ -1,4 +1,5 @@ #include "aare/Interpolator.hpp" +#include "aare/algorithm.hpp" namespace aare { @@ -68,16 +69,17 @@ std::vector Interpolator::interpolate(const ClusterVector& clus photon.y = cluster.y; photon.energy = eta.sum; - //Now do some actual interpolation. - //Find which energy bin the cluster is in - //TODO! Could we use boost-histogram Axis.index here? - ssize_t idx = std::lower_bound(m_energy_bins.begin(), m_energy_bins.end(), photon.energy)-m_energy_bins.begin(); - auto ix = std::lower_bound(m_etabinsx.begin(), m_etabinsx.end(), eta.x)- m_etabinsx.begin(); - auto iy = std::lower_bound(m_etabinsy.begin(), m_etabinsy.end(), eta.y)- m_etabinsy.begin(); - + // auto ie = nearest_index(m_energy_bins, photon.energy)-1; + // auto ix = nearest_index(m_etabinsx, eta.x)-1; + // auto iy = nearest_index(m_etabinsy, eta.y)-1; + //Finding the index of the last element that is smaller + //should work fine as long as we have many bins + auto ie = last_smaller(m_energy_bins, photon.energy); + auto ix = last_smaller(m_etabinsx, eta.x); + auto iy = last_smaller(m_etabinsy, eta.y); + + // fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy); - // ibx=(np.abs(etabinsx - ex)).argmin() #Find out which bin the eta should land in - // iby=(np.abs(etabinsy - ey)).argmin() double dX, dY; int ex, ey; // cBottomLeft = 0, @@ -98,21 +100,16 @@ std::vector Interpolator::interpolate(const ClusterVector& clus dY = -1.; break; case cBottomRight: - dX = 0; + dX = 0.; dY = -1.; break; } - photon.x += m_ietax(ix, iy, idx) + dX + 0.5; - photon.y += m_ietay(ix, iy, idx) + dY + 0.5; - - - // fmt::print("x: {}, y: {}, energy: {}\n", photon.x, photon.y, photon.energy); + photon.x += m_ietax(ix, iy, 0)*2 + dX; + photon.y += m_ietay(ix, iy, 0)*2 + dY; photons.push_back(photon); } }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){ - //TODO! Implement 2x2 interpolation for (size_t i = 0; i(i); Eta2 eta= calculate_eta2(cluster); @@ -123,30 +120,17 @@ std::vector Interpolator::interpolate(const ClusterVector& clus //Now do some actual interpolation. //Find which energy bin the cluster is in - //TODO! Could we use boost-histogram Axis.index here? - ssize_t idx = std::lower_bound(m_energy_bins.begin(), m_energy_bins.end(), photon.energy)-m_energy_bins.begin(); - // auto ix = std::lower_bound(m_etabinsx.begin(), m_etabinsx.end(), eta.x)- m_etabinsx.begin(); - // auto iy = std::lower_bound(m_etabinsy.begin(), m_etabinsy.end(), eta.y)- m_etabinsy.begin(); - // if(ix<0) ix=0; - // if(iy<0) iy=0; - - auto find_index = [](NDArray& etabins, double val){ - auto iter = std::min_element(etabins.begin(), etabins.end(), - [val,etabins](double a, double b) { - return std::abs(a - val) < std::abs(b - val); - }); - return std::distance(etabins.begin(), iter); - }; - auto ix = find_index(m_etabinsx, eta.x)-1; - auto iy = find_index(m_etabinsy, eta.y)-1; - if(ix<0) ix=0; - if(iy<0) iy=0; + // auto ie = nearest_index(m_energy_bins, photon.energy)-1; + // auto ix = nearest_index(m_etabinsx, eta.x)-1; + // auto iy = nearest_index(m_etabinsy, eta.y)-1; + //Finding the index of the last element that is smaller + //should work fine as long as we have many bins + auto ie = last_smaller(m_energy_bins, photon.energy); + auto ix = last_smaller(m_etabinsx, eta.x); + auto iy = last_smaller(m_etabinsy, eta.y); photon.x += m_ietax(ix, iy, 0)*2; //eta goes between 0 and 1 but we could move the hit anywhere in the 2x2 photon.y += m_ietay(ix, iy, 0)*2; - - // photon.x = ix; - // photon.y = idx; photons.push_back(photon); } diff --git a/src/algorithm.test.cpp b/src/algorithm.test.cpp new file mode 100644 index 0000000..9e75eb9 --- /dev/null +++ b/src/algorithm.test.cpp @@ -0,0 +1,63 @@ + + +#include +#include + + +TEST_CASE("Find the closed index in a 1D array", "[algorithm]") { + aare::NDArray arr({5}); + for (size_t i = 0; i < arr.size(); i++) { + arr[i] = i; + } + // arr 0, 1, 2, 3, 4 + REQUIRE(aare::nearest_index(arr, 2.3) == 2); + REQUIRE(aare::nearest_index(arr, 2.6) == 3); + REQUIRE(aare::nearest_index(arr, 45.0) == 4); + REQUIRE(aare::nearest_index(arr, 0.0) == 0); + REQUIRE(aare::nearest_index(arr, -1.0) == 0); +} + +TEST_CASE("Passing integers to nearest_index works"){ + aare::NDArray arr({5}); + for (size_t i = 0; i < arr.size(); i++) { + arr[i] = i; + } + // arr 0, 1, 2, 3, 4 + REQUIRE(aare::nearest_index(arr, 2) == 2); + REQUIRE(aare::nearest_index(arr, 3) == 3); + REQUIRE(aare::nearest_index(arr, 45) == 4); + REQUIRE(aare::nearest_index(arr, 0) == 0); + REQUIRE(aare::nearest_index(arr, -1) == 0); +} + + +TEST_CASE("nearest_index works with std::vector"){ + std::vector vec = {0, 1, 2, 3, 4}; + REQUIRE(aare::nearest_index(vec, 2.123) == 2); + REQUIRE(aare::nearest_index(vec, 2.66) == 3); + REQUIRE(aare::nearest_index(vec, 4555555.0) == 4); + REQUIRE(aare::nearest_index(vec, 0.0) == 0); + REQUIRE(aare::nearest_index(vec, -10.0) == 0); +} + +TEST_CASE("nearest index works with std::array"){ + std::array arr = {0, 1, 2, 3, 4}; + REQUIRE(aare::nearest_index(arr, 2.123) == 2); + REQUIRE(aare::nearest_index(arr, 2.501) == 3); + REQUIRE(aare::nearest_index(arr, 4555555.0) == 4); + REQUIRE(aare::nearest_index(arr, 0.0) == 0); + REQUIRE(aare::nearest_index(arr, -10.0) == 0); +} + + +TEST_CASE("last smaller"){ + aare::NDArray arr({5}); + for (size_t i = 0; i < arr.size(); i++) { + arr[i] = i; + } + // arr 0, 1, 2, 3, 4 + REQUIRE(aare::last_smaller(arr, -10.0) == 0); + REQUIRE(aare::last_smaller(arr, 0.0) == 0); + REQUIRE(aare::last_smaller(arr, 2.3) == 2); + REQUIRE(aare::last_smaller(arr, 253.) == 4); +} \ No newline at end of file