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

modified algo

Browse Source
This commit is contained in:
froejdh_e 2025-03-14 11:07:09 +01:00
parent 3a987319d4
commit 332bdeb02b
6 changed files with 206 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -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

3
include/aare/NDArray.hpp
View File

@ -102,6 +102,9 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, 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

55
include/aare/algorithm.hpp Normal file
View File

@ -0,0 +1,55 @@
#pragma once
#include <algorithm>
#include <array>
#include <vector>
#include <aare/NDArray.hpp>
namespace aare {
/**
* @brief Find the index of the last element smaller than val
* assume a sorted array
*/
template <typename T>
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 <typename T>
size_t last_smaller(const NDArray<T, 1>& arr, T val) {
return last_smaller(arr.begin(), arr.end(), val);
}
template <typename T>
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 <typename T>
size_t nearest_index(const NDArray<T, 1>& arr, T val) {
return nearest_index(arr.begin(), arr.end(), val);
}
template <typename T>
size_t nearest_index(const std::vector<T>& vec, T val) {
return nearest_index(vec.data(), vec.data()+vec.size(), val);
}
template <typename T, size_t N>
size_t nearest_index(const std::array<T,N>& arr, T val) {
return nearest_index(arr.data(), arr.data()+arr.size(), val);
}
} // namespace aare

87
python/examples/play.py
View File

@ -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)
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=}')

58
src/Interpolator.cpp
View File

@ -1,4 +1,5 @@
#include "aare/Interpolator.hpp"
#include "aare/algorithm.hpp"
namespace aare {
@ -68,16 +69,17 @@ std::vector<Photon> Interpolator::interpolate(const ClusterVector<int32_t>& 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<Photon> Interpolator::interpolate(const ClusterVector<int32_t>& 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<clusters.size(); i++){
auto cluster = clusters.at<Cluster2x2>(i);
Eta2 eta= calculate_eta2(cluster);
@ -123,30 +120,17 @@ std::vector<Photon> Interpolator::interpolate(const ClusterVector<int32_t>& 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<double, 1>& 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);
}

63
src/algorithm.test.cpp Normal file
View File

@ -0,0 +1,63 @@
#include <catch2/catch_test_macros.hpp>
#include <aare/algorithm.hpp>
TEST_CASE("Find the closed index in a 1D array", "[algorithm]") {
aare::NDArray<double, 1> 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<int, 1> 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<double> 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<double, 5> 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<double, 1> 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);
}