mirror of
https://github.com/slsdetectorgroup/aare.git
synced 2025-06-24 12:21:07 +02:00
moved anglecalibration stuff to own repository still keep some changes
This commit is contained in:
@ -62,7 +62,6 @@ option(AARE_FETCH_CATCH "Use FetchContent to download catch2" ON)
|
||||
option(AARE_FETCH_JSON "Use FetchContent to download nlohmann::json" ON)
|
||||
option(AARE_FETCH_ZMQ "Use FetchContent to download libzmq" ON)
|
||||
option(AARE_FETCH_LMFIT "Use FetchContent to download lmfit" ON)
|
||||
option(AARE_FETCH_HDF5 "Use FetchContent to download hdf5-devel" OFF)
|
||||
|
||||
|
||||
#Convenience option to use system libraries only (no FetchContent)
|
||||
@ -74,7 +73,6 @@ if(AARE_SYSTEM_LIBRARIES)
|
||||
set(AARE_FETCH_CATCH OFF CACHE BOOL "Disabled FetchContent for catch2" FORCE)
|
||||
set(AARE_FETCH_JSON OFF CACHE BOOL "Disabled FetchContent for nlohmann::json" FORCE)
|
||||
set(AARE_FETCH_ZMQ OFF CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
|
||||
set(AARE_FETCH_HDF5 OFF CACHE BOOL "Disabled FetchContent for hdf5" FORCE)
|
||||
# Still fetch lmfit when setting AARE_SYSTEM_LIBRARIES since this is not available
|
||||
# on conda-forge
|
||||
endif()
|
||||
@ -221,23 +219,6 @@ else()
|
||||
find_package(nlohmann_json 3.11.3 REQUIRED)
|
||||
endif()
|
||||
|
||||
if(AARE_FETCH_HDF5)
|
||||
message(FATAL_ERROR "Fetching HDF5 via FetchContent is not supported here. Please install it via your system.
|
||||
For Ubuntu: sudo apt install libhdf5-dev
|
||||
For Red Hat: sudo dnf install hdf5-devel
|
||||
For MacOS: brew install hdf5")
|
||||
else()
|
||||
find_package(HDF5 QUIET COMPONENTS CXX)
|
||||
if (HDF5_FOUND)
|
||||
message(STATUS "Found HDF5: ${HDF5_INCLUDE_DIRS}")
|
||||
else()
|
||||
message(FATAL_ERROR "HDF5 was NOT found! Please install it via your system
|
||||
For Ubuntu: sudo apt install libhdf5-dev
|
||||
For Red Hat: sudo dnf install hdf5-devel
|
||||
For MacOS: brew install hdf5")
|
||||
endif()
|
||||
endif()
|
||||
|
||||
include(GNUInstallDirs)
|
||||
|
||||
# If conda build, always set lib dir to 'lib'
|
||||
@ -359,7 +340,6 @@ endif()
|
||||
###------------------------------------------------------------------------------------------
|
||||
|
||||
set(PUBLICHEADERS
|
||||
include/aare/AngleCalibration.hpp
|
||||
include/aare/ArrayExpr.hpp
|
||||
include/aare/CalculateEta.hpp
|
||||
include/aare/Cluster.hpp
|
||||
@ -367,7 +347,6 @@ set(PUBLICHEADERS
|
||||
include/aare/ClusterFile.hpp
|
||||
include/aare/CtbRawFile.hpp
|
||||
include/aare/ClusterVector.hpp
|
||||
include/aare/CustomFiles.hpp
|
||||
include/aare/decode.hpp
|
||||
include/aare/defs.hpp
|
||||
include/aare/Dtype.hpp
|
||||
@ -375,14 +354,10 @@ set(PUBLICHEADERS
|
||||
include/aare/Fit.hpp
|
||||
include/aare/FileInterface.hpp
|
||||
include/aare/FilePtr.hpp
|
||||
include/aare/FlatField.hpp
|
||||
include/aare/Frame.hpp
|
||||
include/aare/GainMap.hpp
|
||||
include/aare/geo_helpers.hpp
|
||||
include/aare/Hdf5FileReader.hpp
|
||||
include/aare/JungfrauDataFile.hpp
|
||||
include/aare/MythenDetectorSpecifications.hpp
|
||||
include/aare/MythenFileReader.hpp
|
||||
include/aare/NDArray.hpp
|
||||
include/aare/NDView.hpp
|
||||
include/aare/NumpyFile.hpp
|
||||
@ -398,9 +373,7 @@ set(PUBLICHEADERS
|
||||
|
||||
|
||||
set(SourceFiles
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/AngleCalibration.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/CustomFiles.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/decode.cpp
|
||||
@ -417,7 +390,6 @@ set(SourceFiles
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.cpp
|
||||
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/ifstream_helpers.cpp
|
||||
)
|
||||
@ -436,13 +408,11 @@ target_link_libraries(
|
||||
PUBLIC
|
||||
fmt::fmt
|
||||
nlohmann_json::nlohmann_json
|
||||
HDF5::HDF5
|
||||
${STD_FS_LIB} # from helpers.cmake
|
||||
PRIVATE
|
||||
aare_compiler_flags
|
||||
Threads::Threads
|
||||
$<BUILD_INTERFACE:lmfit>
|
||||
|
||||
)
|
||||
|
||||
set_target_properties(aare_core PROPERTIES
|
||||
@ -457,14 +427,12 @@ endif()
|
||||
if(AARE_TESTS)
|
||||
set(TestSources
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/AngleCalibration.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/decode.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/geo_helpers.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/Hdf5FileReader.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp
|
||||
@ -475,7 +443,6 @@ if(AARE_TESTS)
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinderMT.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/MythenFileReader.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.test.cpp
|
||||
|
||||
@ -1,64 +0,0 @@
|
||||
from slsdet import Mythen3,timingMode,detectorSettings,runStatus,dacIndex,scanParameters
|
||||
#from patterntools.zmqreceiver import ZmqReceiver
|
||||
#from detConf_module import *
|
||||
import numpy as np
|
||||
from slsdet.lookup import view, find
|
||||
from epics import caput, caget
|
||||
import time
|
||||
|
||||
|
||||
d = Mythen3()
|
||||
#d.detsize = [10, 1]
|
||||
|
||||
d.hostname = 'localhost:2000+localhost:2002'
|
||||
|
||||
d.rx_hostname = 'localhost'
|
||||
|
||||
#d.rx_tcpport = 2012
|
||||
|
||||
d.udp_dstport = 5006
|
||||
|
||||
d.udp_dstip = 'auto'
|
||||
|
||||
d.udp_srcip = 'auto'
|
||||
|
||||
#rx=makeReceiver(d)
|
||||
|
||||
nmod=len(d.hostname)
|
||||
print("num modules: ", nmod)
|
||||
|
||||
d.frames = 3
|
||||
d.exptime = 5
|
||||
|
||||
d.fwrite=1
|
||||
#d.counters = [0, 1]
|
||||
print(d.counters)
|
||||
|
||||
d.fpath = "~/Documents/tmp/Flatfieldacquisition"
|
||||
d.fwrite = 1
|
||||
|
||||
#where do i set the strips?
|
||||
#d.detsize = 2 do I have to set this
|
||||
|
||||
d.startReceiver()
|
||||
d.startDetector()
|
||||
for angle in [87,2]:
|
||||
d.fname = 'run_'+str(angle)
|
||||
caput('BL11I-MO-DIFF-01:DELTA.VAL',angle,wait=False)
|
||||
print("moving detector to ",angle)
|
||||
d.acquire()
|
||||
time.sleep(d.exptime)
|
||||
while d.status != runStatus.IDLE:
|
||||
time.sleep(0.01)
|
||||
|
||||
nf=np.min(d.rx_framescaught)
|
||||
ang=caget('BL11I-MO-DIFF-01:DELTA.RBV')
|
||||
print("angle is: ", ang)
|
||||
print("caught frames: ", nf)
|
||||
|
||||
d.stopReceiver()
|
||||
|
||||
#dd, hh = rx[imod].receive_one_frame() what is this - does this exist in cpp?
|
||||
|
||||
|
||||
|
||||
@ -1,164 +0,0 @@
|
||||
#pragma once
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <cstdint>
|
||||
|
||||
#include <fstream>
|
||||
#include <iomanip>
|
||||
#include <iostream>
|
||||
#include <memory>
|
||||
#include <sstream>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
#include "FlatField.hpp"
|
||||
#include "MythenDetectorSpecifications.hpp"
|
||||
#include "MythenFileReader.hpp"
|
||||
#include "NDArray.hpp"
|
||||
|
||||
namespace aare {
|
||||
|
||||
using parameters =
|
||||
std::tuple<std::vector<double>, std::vector<double>, std::vector<double>>;
|
||||
|
||||
class AngleCalibration {
|
||||
|
||||
public:
|
||||
AngleCalibration(
|
||||
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_,
|
||||
std::shared_ptr<FlatField> flat_field_,
|
||||
std::shared_ptr<MythenFileReader> mythen_file_reader_);
|
||||
|
||||
/** set the histogram bin width [degrees] */
|
||||
void set_histogram_bin_width(double bin_width);
|
||||
|
||||
double get_histogram_bin_width() const;
|
||||
|
||||
ssize_t get_new_num_bins() const;
|
||||
|
||||
/** reads the historical Detector Group (DG) parameters from file **/
|
||||
void read_initial_calibration_from_file(const std::string &filename);
|
||||
|
||||
std::vector<double> get_centers() const;
|
||||
std::vector<double> get_conversions() const;
|
||||
|
||||
std::vector<double> get_offsets() const;
|
||||
|
||||
NDView<double, 1> get_new_photon_counts() const;
|
||||
|
||||
NDView<double, 1> get_new_statistical_errors() const;
|
||||
|
||||
/** converts DG parameters to easy EE parameters e.g.geometric
|
||||
* parameters */
|
||||
parameters convert_to_EE_parameters() const;
|
||||
|
||||
std::tuple<double, double, double>
|
||||
convert_to_EE_parameters(const size_t module_index) const;
|
||||
|
||||
std::tuple<double, double, double>
|
||||
convert_to_EE_parameters(const double center, const double conversion,
|
||||
const double offset) const;
|
||||
|
||||
/** converts DG parameters to BC parameters e.g. best computing
|
||||
* parameters */
|
||||
parameters convert_to_BC_parameters() const;
|
||||
|
||||
/**
|
||||
* calculates new histogram with fixed sized angle bins
|
||||
* for several acquisitions at different detector angles for given frame
|
||||
* indices
|
||||
* @param start_frame_index, end_frame_index gives range of frames
|
||||
*/
|
||||
void
|
||||
calculate_fixed_bin_angle_width_histogram(const size_t start_frame_index,
|
||||
const size_t end_frame_index);
|
||||
|
||||
void write_to_file(const std::string &filename,
|
||||
const bool store_nonzero_bins = false,
|
||||
const std::filesystem::path &filepath =
|
||||
std::filesystem::current_path()) const;
|
||||
|
||||
/** calculates diffraction angle from EE module parameters (used in
|
||||
* Beer's Law)
|
||||
* @param strip_index local strip index of module
|
||||
*/
|
||||
double diffraction_angle_from_EE_parameters(
|
||||
const double module_center_distance, const double normal_distance,
|
||||
const double angle, const size_t strip_index,
|
||||
const double distance_to_strip = 0) const;
|
||||
|
||||
/** calculates diffraction angle from EE module parameters (used in
|
||||
* Beer's Law)
|
||||
* @param center module center
|
||||
* @param conversion module conversion
|
||||
* @param offset module offset
|
||||
* @param strip_index local strip index of module
|
||||
* @param distance_to_strip distance to strip given by strip_index and
|
||||
* module -> note needs to be small enough to be in the respective module
|
||||
*/
|
||||
double diffraction_angle_from_DG_parameters(
|
||||
const double center, const double conversion, const double offset,
|
||||
const size_t strip_index, const double distance_to_strip = 0) const;
|
||||
|
||||
/** calculated the strip width expressed as angle [degrees]
|
||||
* @param strip_index local strip index of module
|
||||
*/
|
||||
double angular_strip_width_from_DG_parameters(
|
||||
const double center, const double conversion, const double offset,
|
||||
const size_t local_strip_index) const;
|
||||
|
||||
double angular_strip_width_from_EE_parameters(
|
||||
const double module_center_distance, const double normal_distance,
|
||||
const double angle, const size_t local_strip_index) const;
|
||||
|
||||
protected:
|
||||
/** converts global strip index to local strip index of that module */
|
||||
size_t global_to_local_strip_index_conversion(
|
||||
const size_t global_strip_index) const;
|
||||
|
||||
/**
|
||||
* redistributes photon counts with of histogram using one bin per strip
|
||||
* to histogram with fixed size angle bins
|
||||
* @param frame MythenFrame storing data from image
|
||||
* @param bin_counts accumulate new photon counts
|
||||
* @param new_statistical_weights accumulate new statistical weights
|
||||
* @param new_errors accumulate new_errors
|
||||
*/
|
||||
void redistribute_photon_counts_to_fixed_angle_bins(
|
||||
const MythenFrame &frame, NDView<double, 1> bin_counts,
|
||||
NDView<double, 1> new_statistical_weights, NDView<double, 1> new_errors,
|
||||
NDView<double, 1> inverse_nromalized_flatfield) const;
|
||||
|
||||
private:
|
||||
// TODO: Design maybe have a struct with three vectors, store all three
|
||||
// sets of parameters as member variables
|
||||
|
||||
// TODO: check if interpretation and units are correct
|
||||
// historical DG parameters
|
||||
// TODO change to NDArray
|
||||
std::vector<double> centers; // orthogonal projection of sample onto
|
||||
// detector (given in strip number) [mm]
|
||||
// D/pitch
|
||||
std::vector<double> conversions; // pitch/(normal distance from sample
|
||||
// to detector (R)) [mm]
|
||||
// //used for easy conversion
|
||||
std::vector<double>
|
||||
offsets; // position of strip zero relative to sample [degrees] phi
|
||||
// - 180/pi*D/R TODO: expected an arcsin(D/R)?
|
||||
|
||||
std::shared_ptr<MythenDetectorSpecifications> mythen_detector;
|
||||
|
||||
std::shared_ptr<FlatField> flat_field;
|
||||
|
||||
NDArray<double, 1> new_photon_counts;
|
||||
NDArray<double, 1> new_photon_count_errors;
|
||||
|
||||
double histogram_bin_width = 0.0036; // [degrees]
|
||||
|
||||
ssize_t num_bins{};
|
||||
|
||||
std::shared_ptr<MythenFileReader>
|
||||
mythen_file_reader; // TODO replace by FileInterface ptr
|
||||
};
|
||||
|
||||
} // namespace aare
|
||||
@ -1,59 +0,0 @@
|
||||
#include "Dtype.hpp"
|
||||
#include "FileInterface.hpp"
|
||||
#include <cstdint>
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
namespace aare {
|
||||
|
||||
// TODO a lot to overload
|
||||
class CustomMythenFile : public FileInterface {
|
||||
|
||||
public:
|
||||
CustomMythenFile(const std::string &filename, const ssize_t rows,
|
||||
const ssize_t cols = 1, const std::string &mode = "r");
|
||||
|
||||
~CustomMythenFile();
|
||||
|
||||
Frame read_frame() override;
|
||||
|
||||
Frame read_frame(size_t frame_number) override;
|
||||
|
||||
std::vector<Frame> read_n(size_t n_frames) override;
|
||||
|
||||
void read_into(std::byte *image_buf) override;
|
||||
|
||||
void read_into(std::byte *image_buf, size_t n_frames) override;
|
||||
|
||||
size_t frame_number(size_t frame_index) override;
|
||||
|
||||
size_t bytes_per_frame() override;
|
||||
|
||||
size_t pixels_per_frame() override;
|
||||
|
||||
void seek(size_t frame_number) override;
|
||||
|
||||
size_t tell() override;
|
||||
|
||||
size_t total_frames() const override;
|
||||
|
||||
size_t rows() const override;
|
||||
size_t cols() const override;
|
||||
|
||||
size_t bitdepth() const override;
|
||||
|
||||
DetectorType detector_type() const override;
|
||||
|
||||
private:
|
||||
std::string m_filename{};
|
||||
std::ifstream m_file{};
|
||||
ssize_t m_num_strips{};
|
||||
// uint8_t m_num_counts{}; TODO extend
|
||||
ssize_t m_rows{};
|
||||
ssize_t m_cols{};
|
||||
static const Dtype m_dtype;
|
||||
static const DetectorType det_type = DetectorType::Mythen3;
|
||||
};
|
||||
} // namespace aare
|
||||
@ -1,215 +0,0 @@
|
||||
|
||||
/**
|
||||
* stores flatfield for angle calibration
|
||||
*/
|
||||
|
||||
#pragma once
|
||||
#include <cmath>
|
||||
#include <cstdint>
|
||||
|
||||
#include <filesystem>
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
#include <memory>
|
||||
#include <sstream>
|
||||
#include <string>
|
||||
|
||||
#include "File.hpp"
|
||||
#include "MythenDetectorSpecifications.hpp"
|
||||
#include "NDArray.hpp"
|
||||
|
||||
namespace aare {
|
||||
|
||||
template <class CustomFile> struct custom_file_compatibility {
|
||||
static constexpr bool value =
|
||||
std::is_base_of<FileInterface, CustomFile>::value &&
|
||||
std::is_constructible<CustomFile, std::string, ssize_t, ssize_t,
|
||||
std::string>::value;
|
||||
};
|
||||
|
||||
template <class CustomFile>
|
||||
constexpr bool custom_file_compatibility_v =
|
||||
custom_file_compatibility<CustomFile>::value;
|
||||
|
||||
class FlatField {
|
||||
|
||||
public:
|
||||
FlatField(std::shared_ptr<MythenDetectorSpecifications> mythen_detector_)
|
||||
: mythen_detector(mythen_detector_) {
|
||||
|
||||
flat_field = NDArray<uint32_t, 1>(
|
||||
std::array<ssize_t, 1>{mythen_detector->num_strips()}, 0);
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief sums up the photon counts for multiple acquisitions
|
||||
* @param file_path: path to filesystem - the filesystem should contain
|
||||
* multiple acqisitions for different detector angles acquired using
|
||||
* slsReceiver and mythen3Detector //TODO: constructor needs to be the same
|
||||
* - ugly
|
||||
*/
|
||||
// TODO unsure about design - maybe scientist has to give file with paths
|
||||
// one wants to accumulate - what to do with strange file format?
|
||||
// TODO: adjust for different counters!!!!
|
||||
void create_flatfield_from_rawfilesystem(
|
||||
const std::filesystem::path &file_path) {
|
||||
|
||||
try {
|
||||
for (const auto &file_in_path :
|
||||
std::filesystem::recursive_directory_iterator(file_path)) {
|
||||
if (file_in_path.is_regular_file()) {
|
||||
std::string filename =
|
||||
file_in_path.path().filename().string();
|
||||
if (filename.find("master") != std::string::npos) {
|
||||
File f(file_in_path);
|
||||
auto frames = f.read_n(f.total_frames());
|
||||
for (const auto &frame : frames) {
|
||||
if (frame.rows() * frame.cols() !=
|
||||
mythen_detector->num_strips() *
|
||||
mythen_detector->num_counters()) {
|
||||
throw std::runtime_error(fmt::format(
|
||||
"sizes mismatch. Expect a size of "
|
||||
"{} - frame has a size of {}",
|
||||
mythen_detector->num_strips() *
|
||||
mythen_detector->num_counters(),
|
||||
frame.rows() * frame.cols()));
|
||||
}
|
||||
for (ssize_t row = 0; row < frame.rows(); ++row)
|
||||
for (ssize_t col = 0; col < frame.cols();
|
||||
++col) {
|
||||
flat_field(row * frame.cols() + col) +=
|
||||
*reinterpret_cast<uint32_t *>(
|
||||
frame.pixel_ptr(
|
||||
row,
|
||||
col)); // TODO inefficient as
|
||||
// one has to copy twice
|
||||
// into frame and into
|
||||
// flat_field
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
} catch (const std::filesystem::filesystem_error &e) {
|
||||
std::cerr << "Filesystem error: " << e.what()
|
||||
<< '\n'; // TODO replace with log
|
||||
} catch (const std::exception &e) {
|
||||
std::cerr << "Runtime error: " << e.what() << '\n';
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief sums up the photon counts for multiple acquisitions
|
||||
* @param filelist: path to file that stores the file paths to the aquired
|
||||
* data the list should contain multiple acquisitions for different detector
|
||||
* angles
|
||||
* @tparam CustomFile: Fileclass that inherits from aare::FileInterface
|
||||
* class needs to overload read_frame and constructor needs to have the
|
||||
* signature CustomFile(std::string filename, ssize_t rows, ssize_t cols,
|
||||
* std::string reading_mode) rows, and cols define the intended image size
|
||||
* to read
|
||||
*/
|
||||
template <class CustomFile,
|
||||
typename = std::enable_if_t<
|
||||
custom_file_compatibility_v<CustomFile>, void>>
|
||||
void create_flatfield_from_filelist(const std::filesystem::path &filelist) {
|
||||
std::ifstream file_filelist(filelist);
|
||||
|
||||
try {
|
||||
std::string filename;
|
||||
while (std::getline(file_filelist, filename)) {
|
||||
CustomFile file(filename, mythen_detector->num_strips(), 1);
|
||||
Frame frame = file.read_frame();
|
||||
if (frame.rows() * frame.cols() !=
|
||||
mythen_detector->num_strips()) {
|
||||
throw std::runtime_error(
|
||||
fmt::format("sizes mismatch. Expect a size of "
|
||||
"{} - frame has a size of {}",
|
||||
mythen_detector->num_strips(),
|
||||
frame.rows() * frame.cols()));
|
||||
}
|
||||
for (ssize_t row = 0; row < frame.rows(); ++row)
|
||||
for (ssize_t col = 0; col < frame.cols(); ++col) {
|
||||
flat_field(row * frame.cols() + col) +=
|
||||
*reinterpret_cast<uint32_t *>(frame.pixel_ptr(
|
||||
row,
|
||||
col)); // TODO inefficient as one has to copy
|
||||
// twice into frame and into flat_field
|
||||
}
|
||||
}
|
||||
file_filelist.close();
|
||||
} catch (const std::exception &e) {
|
||||
std::cerr << "Error: " << e.what() << '\n';
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* @tparam CustomFile: Fileclass that inherits from aare::FileInterface
|
||||
* class needs to overload read_frame and constructor needs to have the
|
||||
* signature CustomFile(std::string filename, ssize_t rows, ssize_t cols,
|
||||
* std::string reading_mode) rows, and cols define the intended image size
|
||||
* to read
|
||||
*/
|
||||
|
||||
template <class CustomFile,
|
||||
typename = std::enable_if_t<
|
||||
custom_file_compatibility_v<CustomFile>, void>>
|
||||
void read_flatfield_from_file(const std::string &filename) {
|
||||
CustomFile file(filename, mythen_detector->num_strips(), 1);
|
||||
file.read_into(reinterpret_cast<std::byte *>(flat_field.data()));
|
||||
}
|
||||
|
||||
NDView<uint32_t, 1> get_flatfield() const { return flat_field.view(); }
|
||||
|
||||
// TODO: remove tolerance
|
||||
double calculate_mean(double tolerance) const {
|
||||
auto [sum, count] = std::accumulate(
|
||||
flat_field.begin(), flat_field.end(),
|
||||
std::make_pair<double, ssize_t>(0.0, 0),
|
||||
[&tolerance](std::pair<double, ssize_t> acc, const auto &element) {
|
||||
return element == 0 ? acc
|
||||
: std::make_pair(acc.first + element,
|
||||
acc.second + 1);
|
||||
});
|
||||
|
||||
return sum / count;
|
||||
}
|
||||
|
||||
NDArray<double, 1>
|
||||
inverse_normalized_flatfield(double tolerance = 0.0001) const {
|
||||
double mean = calculate_mean(tolerance);
|
||||
|
||||
NDArray<double, 1> inverse_normalized_flatfield(flat_field.shape());
|
||||
|
||||
for (ssize_t i = 0; i < flat_field.size(); ++i) {
|
||||
inverse_normalized_flatfield[i] =
|
||||
(flat_field[i] <= tolerance ? 0.0 : mean / flat_field[i]);
|
||||
if (inverse_normalized_flatfield[i] < tolerance)
|
||||
mythen_detector->get_bad_channels()[i] = true;
|
||||
}
|
||||
|
||||
return inverse_normalized_flatfield; // TODO: better to have a copy in
|
||||
// this context but unneccessary
|
||||
// for angle calibration code
|
||||
// maybe provide inplace and copy option
|
||||
// maybe store as member variable access with view
|
||||
}
|
||||
|
||||
NDArray<double, 1> normalized_flatfield(double tolerance = 0.0001) const {
|
||||
double mean = calculate_mean(tolerance);
|
||||
|
||||
NDArray<double, 1> normalized_flatfield(flat_field.shape());
|
||||
|
||||
for (ssize_t i = 0; i < flat_field.size(); ++i) {
|
||||
normalized_flatfield[i] = (flat_field[i] == flat_field[i] / mean);
|
||||
if (normalized_flatfield[i] < tolerance)
|
||||
mythen_detector->get_bad_channels()[i] = true;
|
||||
}
|
||||
return normalized_flatfield;
|
||||
}
|
||||
|
||||
private:
|
||||
NDArray<uint32_t, 1> flat_field; // TODO: should be 2d
|
||||
std::shared_ptr<MythenDetectorSpecifications> mythen_detector;
|
||||
};
|
||||
} // namespace aare
|
||||
@ -1,212 +0,0 @@
|
||||
/************************************************
|
||||
* @file HD5FFileReader.hpp
|
||||
* @short HDF5FileReader based on H5File object
|
||||
***********************************************/
|
||||
|
||||
#pragma once
|
||||
|
||||
#include "Frame.hpp"
|
||||
#include "NDArray.hpp"
|
||||
#include <H5Cpp.h>
|
||||
#include <array>
|
||||
#include <cxxabi.h>
|
||||
#include <optional>
|
||||
|
||||
namespace aare {
|
||||
|
||||
// return std::type_info
|
||||
inline const std::type_info &deduce_cpp_type(const H5::DataType datatype) {
|
||||
if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_UINT8.getId())) {
|
||||
return typeid(uint8_t);
|
||||
} else if (H5Tequal(datatype.getId(),
|
||||
H5::PredType::NATIVE_UINT16.getId())) {
|
||||
return typeid(uint16_t);
|
||||
} else if (H5Tequal(datatype.getId(),
|
||||
H5::PredType::NATIVE_UINT32.getId())) {
|
||||
return typeid(uint32_t);
|
||||
} else if (H5Tequal(datatype.getId(),
|
||||
H5::PredType::NATIVE_UINT64.getId())) {
|
||||
return typeid(uint64_t);
|
||||
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT8.getId())) {
|
||||
return typeid(int8_t);
|
||||
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT16.getId())) {
|
||||
return typeid(int16_t);
|
||||
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT32.getId())) {
|
||||
return typeid(int32_t);
|
||||
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT64.getId())) {
|
||||
return typeid(int64_t);
|
||||
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT.getId())) {
|
||||
return typeid(int);
|
||||
} else if (H5Tequal(datatype.getId(), H5::PredType::IEEE_F64LE.getId())) {
|
||||
return typeid(double);
|
||||
} else if (H5Tequal(datatype.getId(), H5::PredType::IEEE_F32LE.getId())) {
|
||||
return typeid(float);
|
||||
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_FLOAT.getId())) {
|
||||
return typeid(float);
|
||||
} else if (H5Tequal(datatype.getId(),
|
||||
H5::PredType::NATIVE_DOUBLE.getId())) {
|
||||
return typeid(float);
|
||||
} else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_CHAR.getId()) &&
|
||||
datatype.getId() == H5::PredType::NATIVE_CHAR.getId()) {
|
||||
return typeid(char);
|
||||
} else {
|
||||
throw std::runtime_error("c++ type cannot be deduced");
|
||||
}
|
||||
}
|
||||
|
||||
struct Subset {
|
||||
std::vector<hsize_t> shape;
|
||||
std::vector<hsize_t> offset; // index where data subset should start
|
||||
};
|
||||
|
||||
class HDF5Dataset {
|
||||
|
||||
public:
|
||||
HDF5Dataset(const std::string &datasetname_, const H5::DataSet dataset_)
|
||||
: datasetname(datasetname_), dataset(dataset_) {
|
||||
datatype = dataset.getDataType();
|
||||
|
||||
cpp_type = &deduce_cpp_type(datatype);
|
||||
|
||||
dataspace = dataset.getSpace();
|
||||
rank = dataspace.getSimpleExtentNdims(); // number of dimensions
|
||||
|
||||
shape.resize(rank);
|
||||
dataspace.getSimpleExtentDims(shape.data(), nullptr);
|
||||
}
|
||||
|
||||
hsize_t get_shape(ssize_t index) const { return shape[index]; }
|
||||
|
||||
std::vector<hsize_t> get_shape() const { return shape; }
|
||||
|
||||
H5::DataType get_datatype() const { return datatype; }
|
||||
|
||||
const std::type_info *get_cpp_type() const { return cpp_type; }
|
||||
|
||||
/**
|
||||
* Reads subset of dataset into the buffer
|
||||
* e.g. to read one 2d frame pass Subset({shape[1], shape[2]}, {frame_index,
|
||||
* 0,0})
|
||||
*/
|
||||
void
|
||||
read_into_buffer(std::byte *buffer,
|
||||
std::optional<const Subset> subset = std::nullopt) const {
|
||||
|
||||
if (subset) {
|
||||
// TODO treat scalar cases
|
||||
if (static_cast<ssize_t>(subset->offset.size()) != rank) {
|
||||
throw std::runtime_error("provide an offset for" +
|
||||
std::to_string(rank) + "dimensions");
|
||||
}
|
||||
for (ssize_t i = 0; i < rank; ++i) {
|
||||
hsize_t size =
|
||||
i < rank - static_cast<ssize_t>(subset->shape.size())
|
||||
? 0
|
||||
: subset->shape[i - (rank - subset->shape.size())];
|
||||
if ((size + subset->offset[i]) > shape[i]) {
|
||||
throw std::runtime_error(
|
||||
"subset is too large or offset is out of bounds");
|
||||
}
|
||||
}
|
||||
|
||||
H5::DataSpace memspace(static_cast<int>(subset->shape.size()),
|
||||
subset->shape.data());
|
||||
|
||||
dataspace.selectHyperslab(H5S_SELECT_SET, subset->shape.data(),
|
||||
subset->offset.data());
|
||||
dataset.read(buffer, datatype, memspace, dataspace);
|
||||
} else {
|
||||
dataset.read(buffer, datatype);
|
||||
}
|
||||
}
|
||||
|
||||
Frame store_as_frame() const {
|
||||
uint32_t rows{}, cols{};
|
||||
if (rank == 1) {
|
||||
rows = 1;
|
||||
// TODO overflow
|
||||
cols = static_cast<uint32_t>(shape[0]);
|
||||
} else if (rank == 2) {
|
||||
rows = static_cast<uint32_t>(shape[0]);
|
||||
cols = static_cast<uint32_t>(shape[1]);
|
||||
} else {
|
||||
throw std::runtime_error("Frame only supports 2d images");
|
||||
}
|
||||
|
||||
Frame frame(rows, cols, Dtype(*cpp_type));
|
||||
|
||||
read_into_buffer(frame.data());
|
||||
|
||||
return frame;
|
||||
}
|
||||
|
||||
template <typename T, ssize_t NDim>
|
||||
NDArray<T, NDim> store_as_ndarray() const {
|
||||
if (NDim != rank) {
|
||||
std::cout
|
||||
<< "Warning: dataset dimension and array dimension mismatch"
|
||||
<< std::endl; // TODO: replace with log - still valid if we
|
||||
// want subset
|
||||
}
|
||||
if (typeid(T) != *cpp_type) {
|
||||
std::cout << "Warning: dataset and array type mismatch"
|
||||
<< std::endl;
|
||||
}
|
||||
std::array<ssize_t, NDim> array_shape{};
|
||||
std::transform(
|
||||
shape.begin(), shape.end(), array_shape.begin(),
|
||||
[](const auto dim) { return static_cast<ssize_t>(dim); });
|
||||
|
||||
aare::NDArray<T, NDim> dataset_array(array_shape);
|
||||
|
||||
read_into_buffer(reinterpret_cast<std::byte *>(dataset_array.data()));
|
||||
|
||||
return dataset_array;
|
||||
}
|
||||
|
||||
// getMemDataSize()
|
||||
|
||||
private:
|
||||
std::string datasetname{};
|
||||
H5::DataSet dataset;
|
||||
H5::DataSpace dataspace;
|
||||
H5::DataType datatype;
|
||||
const std::type_info *cpp_type;
|
||||
ssize_t rank{};
|
||||
std::vector<hsize_t> shape{};
|
||||
};
|
||||
|
||||
class HDF5FileReader {
|
||||
|
||||
public:
|
||||
HDF5FileReader() = default;
|
||||
|
||||
void open_file(const std::string &filename_) {
|
||||
filename = filename_;
|
||||
try {
|
||||
file = H5::H5File(filename, H5F_ACC_RDONLY);
|
||||
} catch (H5::Exception &e) {
|
||||
std::cerr << "Error: " << e.getDetailMsg() << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
void close_file() { file.close(); }
|
||||
|
||||
HDF5Dataset get_dataset(const std::string &dataset_name) const {
|
||||
H5::DataSet dataset;
|
||||
try {
|
||||
dataset = file.openDataSet(dataset_name);
|
||||
} catch (H5::Exception &e) {
|
||||
std::cerr << "Error: " << e.getDetailMsg() << std::endl;
|
||||
}
|
||||
|
||||
// TODO use optional to handle error
|
||||
return HDF5Dataset(dataset_name, dataset);
|
||||
}
|
||||
|
||||
private:
|
||||
std::string filename{};
|
||||
H5::H5File file;
|
||||
};
|
||||
|
||||
} // namespace aare
|
||||
@ -1,157 +0,0 @@
|
||||
#pragma once
|
||||
#include <cmath>
|
||||
#include <cstdint>
|
||||
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
#include <sstream>
|
||||
#include <string>
|
||||
|
||||
#include "NDArray.hpp"
|
||||
|
||||
namespace aare {
|
||||
|
||||
class MythenDetectorSpecifications {
|
||||
|
||||
public:
|
||||
// TODO: constructor that reads from a config file
|
||||
|
||||
MythenDetectorSpecifications() {
|
||||
num_strips_ = max_modules_ * strips_per_module_;
|
||||
|
||||
num_connected_modules_ = max_modules_;
|
||||
|
||||
bad_channels =
|
||||
NDArray<bool, 1>(std::array<ssize_t, 1>{num_strips_}, false);
|
||||
|
||||
connected_modules = NDArray<bool, 1>(
|
||||
std::array<ssize_t, 1>{static_cast<ssize_t>(max_modules_)}, true);
|
||||
}
|
||||
|
||||
MythenDetectorSpecifications(const size_t max_modules,
|
||||
const double exposure_time,
|
||||
const double num_counters = 1,
|
||||
double bloffset = 1.532)
|
||||
: max_modules_(max_modules), num_counters_(num_counters),
|
||||
exposure_time_(exposure_time), bloffset_(bloffset) {
|
||||
num_strips_ = max_modules_ * strips_per_module_;
|
||||
|
||||
num_connected_modules_ = max_modules_;
|
||||
|
||||
bad_channels =
|
||||
NDArray<bool, 1>(std::array<ssize_t, 1>{num_strips_}, false);
|
||||
|
||||
connected_modules = NDArray<bool, 1>(
|
||||
std::array<ssize_t, 1>{static_cast<ssize_t>(max_modules_)}, true);
|
||||
}
|
||||
|
||||
// TODO templated on filereader
|
||||
void read_bad_channels_from_file(const std::string &filename) {
|
||||
std::string line;
|
||||
|
||||
try {
|
||||
std::ifstream file(filename, std::ios_base::in);
|
||||
if (!file.good()) {
|
||||
throw std::logic_error("file does not exist");
|
||||
}
|
||||
|
||||
while (std::getline(file, line)) {
|
||||
std::size_t pos = line.find("-");
|
||||
|
||||
if (pos == std::string::npos) {
|
||||
bad_channels(std::stoi(line)) = true;
|
||||
} else {
|
||||
size_t line_size = line.size();
|
||||
for (int i = std::stoi(line.substr(0, pos));
|
||||
i <= std::stoi(line.substr(pos + 1, line_size - pos));
|
||||
++i)
|
||||
bad_channels(i) = true;
|
||||
}
|
||||
}
|
||||
|
||||
file.close();
|
||||
} catch (const std::exception &e) {
|
||||
std::cerr << "Error: " << e.what() << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
// TODO template on filereader
|
||||
void read_unconnected_modules_from_file(const std::string &filename) {
|
||||
std::string line;
|
||||
|
||||
try {
|
||||
std::ifstream file(filename, std::ios_base::in);
|
||||
if (!file.good()) {
|
||||
throw std::logic_error("file does not exist");
|
||||
}
|
||||
|
||||
std::stringstream file_buffer;
|
||||
file_buffer << file.rdbuf();
|
||||
|
||||
file_buffer >> line;
|
||||
num_connected_modules_ -= std::stoi(line);
|
||||
|
||||
while (file_buffer >> line) {
|
||||
size_t module = std::stoi(line);
|
||||
connected_modules[module] = false;
|
||||
for (size_t i = module * strips_per_module_;
|
||||
i < (module + 1) * strips_per_module_; ++i)
|
||||
bad_channels[i] = true;
|
||||
}
|
||||
} catch (const std::exception &e) {
|
||||
std::cerr << "Error: " << e.what() << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
NDView<bool, 1> get_bad_channels() const { return bad_channels.view(); }
|
||||
|
||||
NDView<bool, 1> get_connected_modules() const {
|
||||
return connected_modules.view();
|
||||
}
|
||||
|
||||
static constexpr double pitch() { return pitch_; }
|
||||
|
||||
static constexpr size_t strips_per_module() { return strips_per_module_; }
|
||||
|
||||
size_t max_modules() const { return max_modules_; }
|
||||
|
||||
size_t num_counters() const { return num_counters_; }
|
||||
|
||||
double exposure_time() const { return exposure_time_; }
|
||||
|
||||
double bloffset() const { return bloffset_; }
|
||||
|
||||
double dtt0() const { return dtt0_; }
|
||||
|
||||
static constexpr double min_angle() { return min_angle_; }
|
||||
|
||||
static constexpr double max_angle() { return max_angle_; }
|
||||
|
||||
ssize_t num_strips() const { return num_strips_; }
|
||||
|
||||
private:
|
||||
static constexpr size_t strips_per_module_ = 1280;
|
||||
static constexpr double pitch_ = 0.05; // strip width [mm]
|
||||
static constexpr double min_angle_ =
|
||||
-180.0; // maybe shoudnt be static but configurable
|
||||
static constexpr double max_angle_ = 180.0;
|
||||
static constexpr double dtt0_ =
|
||||
0.0; // No idea what this is - probably configurable
|
||||
|
||||
size_t max_modules_ = 48;
|
||||
|
||||
size_t num_counters_ = 1;
|
||||
|
||||
double exposure_time_ = 5.0; // TODO: could read from acquired file but
|
||||
// maybe should be configurable
|
||||
double bloffset_ = 1.532; // what is this? detector offset relative to what?
|
||||
|
||||
size_t num_connected_modules_{};
|
||||
|
||||
ssize_t num_strips_{};
|
||||
|
||||
NDArray<bool, 1> bad_channels;
|
||||
NDArray<bool, 1> connected_modules; // connected modules
|
||||
};
|
||||
|
||||
} // namespace aare
|
||||
@ -1,82 +0,0 @@
|
||||
/************************************************
|
||||
* @file MythenFileReader.hpp
|
||||
* @short minimal file reader to read mythen files
|
||||
***********************************************/
|
||||
|
||||
#include <bitset>
|
||||
#include <filesystem>
|
||||
#include <string>
|
||||
|
||||
#include "Hdf5FileReader.hpp"
|
||||
#include "NDArray.hpp"
|
||||
|
||||
namespace aare {
|
||||
|
||||
struct MythenFrame {
|
||||
NDArray<uint32_t, 1> photon_counts;
|
||||
double detector_angle{};
|
||||
// double reference_intensity{}; not needed
|
||||
std::array<uint8_t, 3> channel_mask{};
|
||||
};
|
||||
|
||||
/** minimal version for a mythen file reader */
|
||||
class MythenFileReader : public HDF5FileReader {
|
||||
|
||||
public:
|
||||
MythenFileReader(const std::filesystem::path &file_path_,
|
||||
const std::string &file_prefix_)
|
||||
: m_base_path(file_path_), file_prefix(file_prefix_) {};
|
||||
|
||||
MythenFrame read_frame(ssize_t frame_index) {
|
||||
// TODO not a good design fixed number of digits in file name for frame
|
||||
// number -> pad with zeros
|
||||
// not even sure if files have the same name
|
||||
std::string current_file_name =
|
||||
m_base_path / (file_prefix + std::to_string(frame_index) + ".h5");
|
||||
|
||||
MythenFrame frame;
|
||||
open_file(current_file_name);
|
||||
|
||||
auto dataset_photon_count =
|
||||
get_dataset("/entry/instrument/detector/data");
|
||||
|
||||
frame.photon_counts =
|
||||
dataset_photon_count.store_as_ndarray<uint32_t, 1>();
|
||||
|
||||
++frame.photon_counts; // Why though?
|
||||
|
||||
auto dataset_detector_angle =
|
||||
get_dataset("/entry/instrument/NDAttributes/DetectorAngle");
|
||||
|
||||
dataset_detector_angle.read_into_buffer(
|
||||
reinterpret_cast<std::byte *>(&frame.detector_angle));
|
||||
|
||||
auto dataset_channel_number =
|
||||
get_dataset("/entry/instrument/NDAttributes/CounterMask");
|
||||
|
||||
uint8_t channel_number;
|
||||
|
||||
dataset_channel_number.read_into_buffer(
|
||||
reinterpret_cast<std::byte *>(&channel_number));
|
||||
|
||||
std::bitset<3> binary_channel_numbers(channel_number); // 1 0 0
|
||||
|
||||
// binary_channel_numbers.flip(); // TODO not sure where most
|
||||
// significant
|
||||
// bit is ask Anna again
|
||||
|
||||
frame.channel_mask = std::array<uint8_t, 3>{binary_channel_numbers[0],
|
||||
binary_channel_numbers[1],
|
||||
binary_channel_numbers[2]};
|
||||
|
||||
close_file();
|
||||
|
||||
return frame;
|
||||
}
|
||||
|
||||
private:
|
||||
std::filesystem::path m_base_path{};
|
||||
std::string file_prefix{};
|
||||
};
|
||||
|
||||
} // namespace aare
|
||||
@ -1,377 +0,0 @@
|
||||
#include "aare/AngleCalibration.hpp"
|
||||
|
||||
namespace aare {
|
||||
|
||||
AngleCalibration::AngleCalibration(
|
||||
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_,
|
||||
std::shared_ptr<FlatField> flat_field_,
|
||||
std::shared_ptr<MythenFileReader> mythen_file_reader_)
|
||||
: mythen_detector(mythen_detector_), flat_field(flat_field_),
|
||||
mythen_file_reader(mythen_file_reader_) {
|
||||
centers.reserve(mythen_detector->max_modules());
|
||||
conversions.reserve(mythen_detector->max_modules());
|
||||
offsets.reserve(mythen_detector->max_modules());
|
||||
|
||||
num_bins = std::floor(mythen_detector->max_angle() / histogram_bin_width) -
|
||||
std::floor(mythen_detector->min_angle() / histogram_bin_width) +
|
||||
1; // TODO only works if negative
|
||||
// and positive angle
|
||||
}
|
||||
|
||||
void AngleCalibration::set_histogram_bin_width(double bin_width) {
|
||||
histogram_bin_width = bin_width;
|
||||
|
||||
num_bins = std::floor(mythen_detector->max_angle() / histogram_bin_width) -
|
||||
std::floor(mythen_detector->min_angle() / histogram_bin_width) +
|
||||
1; // TODO only works if negative
|
||||
// and positive angle
|
||||
}
|
||||
|
||||
double AngleCalibration::get_histogram_bin_width() const {
|
||||
return histogram_bin_width;
|
||||
}
|
||||
|
||||
ssize_t AngleCalibration::get_new_num_bins() const { return num_bins; }
|
||||
|
||||
std::vector<double> AngleCalibration::get_centers() const { return centers; }
|
||||
|
||||
std::vector<double> AngleCalibration::get_conversions() const {
|
||||
return conversions;
|
||||
}
|
||||
|
||||
std::vector<double> AngleCalibration::get_offsets() const { return offsets; }
|
||||
|
||||
NDView<double, 1> AngleCalibration::get_new_photon_counts() const {
|
||||
return new_photon_counts.view();
|
||||
}
|
||||
|
||||
NDView<double, 1> AngleCalibration::get_new_statistical_errors() const {
|
||||
return new_photon_count_errors.view();
|
||||
}
|
||||
|
||||
void AngleCalibration::read_initial_calibration_from_file(
|
||||
const std::string &filename) {
|
||||
|
||||
std::string line;
|
||||
uint32_t module_number{};
|
||||
|
||||
try {
|
||||
std::ifstream file(filename, std::ios_base::in);
|
||||
if (!file.good()) {
|
||||
throw std::logic_error("file does not exist");
|
||||
}
|
||||
|
||||
std::stringstream file_buffer;
|
||||
file_buffer << file.rdbuf();
|
||||
|
||||
while (file_buffer >> line) {
|
||||
if (line == "module") {
|
||||
file_buffer >> line;
|
||||
module_number = std::stoi(line);
|
||||
}
|
||||
if (line == "center") {
|
||||
file_buffer >> line;
|
||||
centers.insert(centers.begin() + module_number,
|
||||
std::stod(line));
|
||||
}
|
||||
if (line == "conversion") {
|
||||
file_buffer >> line;
|
||||
conversions.insert(conversions.begin() + module_number,
|
||||
std::stod(line));
|
||||
}
|
||||
if (line == "offset") {
|
||||
file_buffer >> line;
|
||||
offsets.insert(offsets.begin() + module_number,
|
||||
std::stod(line));
|
||||
}
|
||||
}
|
||||
|
||||
file.close();
|
||||
} catch (const std::exception &e) {
|
||||
std::cerr << "Error: " << e.what() << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
parameters AngleCalibration::convert_to_EE_parameters() const {
|
||||
|
||||
// normal distance between sample and detector (R)
|
||||
std::vector<double> normal_distances(centers.size());
|
||||
// distances between intersection point of sample normal and module origin
|
||||
// (D)
|
||||
std::vector<double> module_center_distances(centers.size());
|
||||
// angles between undiffracted beam and orthogonal sample projection on
|
||||
// detector (phi)
|
||||
std::vector<double> angles(centers.size());
|
||||
|
||||
for (size_t i = 0; i < centers.size(); ++i) {
|
||||
auto [module_center_distance, normal_distance, angle] =
|
||||
convert_to_EE_parameters(i);
|
||||
normal_distances[i] = normal_distance;
|
||||
module_center_distances[i] = module_center_distance;
|
||||
angles[i] = angle;
|
||||
}
|
||||
|
||||
return std::make_tuple(module_center_distances, normal_distances, angles);
|
||||
}
|
||||
|
||||
std::tuple<double, double, double>
|
||||
AngleCalibration::convert_to_EE_parameters(const size_t module_index) const {
|
||||
return convert_to_EE_parameters(centers[module_index],
|
||||
conversions[module_index],
|
||||
offsets[module_index]);
|
||||
}
|
||||
|
||||
std::tuple<double, double, double> AngleCalibration::convert_to_EE_parameters(
|
||||
const double center, const double conversion, const double offset) const {
|
||||
const double module_center_distance =
|
||||
center * MythenDetectorSpecifications::pitch();
|
||||
const double normal_distance =
|
||||
MythenDetectorSpecifications::pitch() / std::abs(conversion);
|
||||
const double angle = offset + 180.0 / M_PI * center * std::abs(conversion);
|
||||
|
||||
return std::make_tuple(module_center_distance, normal_distance, angle);
|
||||
}
|
||||
|
||||
size_t AngleCalibration::global_to_local_strip_index_conversion(
|
||||
const size_t global_strip_index) const {
|
||||
const size_t module_index =
|
||||
global_strip_index / MythenDetectorSpecifications::strips_per_module();
|
||||
// local strip index in module
|
||||
size_t local_strip_index =
|
||||
global_strip_index -
|
||||
module_index * MythenDetectorSpecifications::strips_per_module();
|
||||
// switch if indexing is in clock-wise direction
|
||||
local_strip_index =
|
||||
std::signbit(conversions[module_index])
|
||||
? MythenDetectorSpecifications::strips_per_module() - 1 -
|
||||
local_strip_index
|
||||
: local_strip_index;
|
||||
|
||||
return local_strip_index;
|
||||
}
|
||||
|
||||
/*
|
||||
parameters
|
||||
AngleCalibration::convert_to_BC_parameters() {}
|
||||
*/
|
||||
|
||||
double AngleCalibration::diffraction_angle_from_DG_parameters(
|
||||
const double center, const double conversion, const double offset,
|
||||
const size_t strip_index, const double distance_to_strip) const {
|
||||
|
||||
return offset + 180.0 / M_PI *
|
||||
(center * std::abs(conversion) -
|
||||
atan((center - (strip_index + distance_to_strip)) *
|
||||
std::abs(conversion)));
|
||||
}
|
||||
|
||||
double AngleCalibration::diffraction_angle_from_EE_parameters(
|
||||
const double module_center_distance, const double normal_distance,
|
||||
const double angle, const size_t strip_index,
|
||||
const double distance_to_strip) const {
|
||||
|
||||
return angle - 180.0 / M_PI *
|
||||
atan((module_center_distance -
|
||||
MythenDetectorSpecifications::pitch() *
|
||||
(strip_index + distance_to_strip)) /
|
||||
normal_distance); // TODO: why is it minus
|
||||
// is it defined counter
|
||||
// clockwise? thought
|
||||
// should have a flipped
|
||||
// sign
|
||||
}
|
||||
|
||||
double AngleCalibration::angular_strip_width_from_DG_parameters(
|
||||
const double center, const double conversion, const double offset,
|
||||
const size_t local_strip_index) const {
|
||||
|
||||
return std::abs(diffraction_angle_from_DG_parameters(
|
||||
center, conversion, offset, local_strip_index, -0.5) -
|
||||
diffraction_angle_from_DG_parameters(
|
||||
center, conversion, offset, local_strip_index, 0.5));
|
||||
}
|
||||
|
||||
double AngleCalibration::angular_strip_width_from_EE_parameters(
|
||||
const double module_center_distance, const double normal_distance,
|
||||
const double angle, const size_t local_strip_index) const {
|
||||
|
||||
return std::abs(diffraction_angle_from_EE_parameters(
|
||||
module_center_distance, normal_distance, angle,
|
||||
local_strip_index, -0.5) -
|
||||
diffraction_angle_from_EE_parameters(
|
||||
module_center_distance, normal_distance, angle,
|
||||
local_strip_index, 0.5));
|
||||
|
||||
// TODO: again not sure about division order - taking abs anyway
|
||||
}
|
||||
|
||||
void AngleCalibration::calculate_fixed_bin_angle_width_histogram(
|
||||
const size_t start_frame_index, const size_t end_frame_index) {
|
||||
|
||||
new_photon_counts = NDArray<double, 1>(std::array<ssize_t, 1>{num_bins});
|
||||
|
||||
new_photon_count_errors =
|
||||
NDArray<double, 1>(std::array<ssize_t, 1>{num_bins});
|
||||
|
||||
// TODO: maybe group these into a 2d array - better cache usage
|
||||
NDArray<double, 1> bin_counts(std::array<ssize_t, 1>{num_bins}, 0.0);
|
||||
NDArray<double, 1> new_statistical_weights(std::array<ssize_t, 1>{num_bins},
|
||||
0.0);
|
||||
NDArray<double, 1> new_errors(std::array<ssize_t, 1>{num_bins}, 0.0);
|
||||
|
||||
NDArray<double, 1> inverse_normalized_flatfield =
|
||||
flat_field->inverse_normalized_flatfield();
|
||||
|
||||
for (size_t frame_index = start_frame_index; frame_index < end_frame_index;
|
||||
++frame_index) {
|
||||
MythenFrame frame = mythen_file_reader->read_frame(frame_index);
|
||||
redistribute_photon_counts_to_fixed_angle_bins(
|
||||
frame, bin_counts.view(), new_statistical_weights.view(),
|
||||
new_errors.view(), inverse_normalized_flatfield.view());
|
||||
}
|
||||
|
||||
for (ssize_t i = 0; i < new_photon_counts.size(); ++i) {
|
||||
new_photon_counts[i] = (new_statistical_weights[i] <=
|
||||
std::numeric_limits<double>::epsilon())
|
||||
? 0.
|
||||
: bin_counts[i] / new_statistical_weights[i];
|
||||
new_photon_count_errors[i] =
|
||||
(bin_counts[i] <= std::numeric_limits<double>::epsilon())
|
||||
? 0.
|
||||
: 1.0 / std::sqrt(bin_counts[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
|
||||
const MythenFrame &frame, NDView<double, 1> bin_counts,
|
||||
NDView<double, 1> new_statistical_weights, NDView<double, 1> new_errors,
|
||||
NDView<double, 1> inverse_normalized_flatfield) const {
|
||||
|
||||
ssize_t channel = 0; // TODO handle mask - FlatField still 1d
|
||||
|
||||
if (frame.photon_counts.shape()[0] != mythen_detector->num_strips()) {
|
||||
throw std::runtime_error("wrong number of strips read");
|
||||
}
|
||||
|
||||
ssize_t num_bins1 = mythen_detector->min_angle() / histogram_bin_width;
|
||||
ssize_t num_bins2 = mythen_detector->max_angle() / histogram_bin_width;
|
||||
|
||||
std::cout << "position: " << frame.detector_angle
|
||||
<< std::endl; // replace with log
|
||||
|
||||
double exposure_rate = 1. / mythen_detector->exposure_time();
|
||||
|
||||
for (ssize_t strip_index = 0; strip_index < mythen_detector->num_strips();
|
||||
++strip_index) {
|
||||
|
||||
size_t module_index =
|
||||
strip_index / MythenDetectorSpecifications::strips_per_module();
|
||||
|
||||
if (mythen_detector->get_bad_channels()[strip_index] ||
|
||||
!mythen_detector->get_connected_modules()[module_index])
|
||||
continue;
|
||||
|
||||
double poisson_error = std::sqrt(frame.photon_counts(strip_index)) *
|
||||
inverse_normalized_flatfield(strip_index) *
|
||||
exposure_rate;
|
||||
|
||||
double corrected_photon_count =
|
||||
frame.photon_counts(strip_index) *
|
||||
inverse_normalized_flatfield(strip_index) * exposure_rate;
|
||||
|
||||
size_t local_strip_index =
|
||||
global_to_local_strip_index_conversion(strip_index);
|
||||
|
||||
double diffraction_angle = diffraction_angle_from_DG_parameters(
|
||||
centers[module_index], conversions[module_index],
|
||||
offsets[module_index], local_strip_index);
|
||||
|
||||
diffraction_angle += (frame.detector_angle + mythen_detector->dtt0() +
|
||||
mythen_detector->bloffset());
|
||||
|
||||
if (diffraction_angle < mythen_detector->min_angle() ||
|
||||
diffraction_angle > mythen_detector->max_angle())
|
||||
continue;
|
||||
|
||||
double angle_covered_by_strip = angular_strip_width_from_DG_parameters(
|
||||
centers[module_index], conversions[module_index],
|
||||
offsets[module_index], local_strip_index);
|
||||
|
||||
double photon_count_per_bin = histogram_bin_width *
|
||||
corrected_photon_count /
|
||||
angle_covered_by_strip;
|
||||
double error_photon_count_per_bin =
|
||||
histogram_bin_width * poisson_error / angle_covered_by_strip;
|
||||
|
||||
double statistical_weights =
|
||||
1.0 / std::pow(error_photon_count_per_bin, 2); // 1./sigma²
|
||||
|
||||
double strip_boundary_left =
|
||||
diffraction_angle - 0.5 * angle_covered_by_strip;
|
||||
double strip_boundary_right =
|
||||
diffraction_angle + 0.5 * angle_covered_by_strip;
|
||||
|
||||
ssize_t left_bin_index = std::max(
|
||||
num_bins1,
|
||||
static_cast<ssize_t>(
|
||||
std::floor(strip_boundary_left / histogram_bin_width) - 1));
|
||||
ssize_t right_bin_index = std::min(
|
||||
num_bins2,
|
||||
static_cast<ssize_t>(
|
||||
std::ceil(strip_boundary_right / histogram_bin_width) + 1));
|
||||
|
||||
// TODO should it be < or <=
|
||||
for (ssize_t bin = left_bin_index; bin <= right_bin_index; ++bin) {
|
||||
double bin_coverage = std::min(strip_boundary_right,
|
||||
(bin + 0.5) * histogram_bin_width) -
|
||||
std::max(strip_boundary_left,
|
||||
(bin - 0.5) * histogram_bin_width);
|
||||
|
||||
double bin_coverage_factor = bin_coverage / histogram_bin_width;
|
||||
|
||||
ssize_t bin_index = bin - num_bins1;
|
||||
// TODO: maybe have this threshold configurable
|
||||
if (bin_coverage >= 0.0001) {
|
||||
new_statistical_weights(bin_index) +=
|
||||
statistical_weights * bin_coverage_factor;
|
||||
bin_counts(bin_index) += statistical_weights *
|
||||
bin_coverage_factor *
|
||||
photon_count_per_bin;
|
||||
new_errors(bin_index) += statistical_weights *
|
||||
bin_coverage_factor *
|
||||
std::pow(photon_count_per_bin, 2);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void AngleCalibration::write_to_file(
|
||||
const std::string &filename, const bool store_nonzero_bins,
|
||||
const std::filesystem::path &filepath) const {
|
||||
|
||||
std::ofstream output_file(filepath / filename);
|
||||
|
||||
if (!output_file) {
|
||||
std::cerr << "Error opening file!"
|
||||
<< std::endl; // TODO: replace with log
|
||||
}
|
||||
|
||||
output_file << std::fixed << std::setprecision(15);
|
||||
|
||||
for (ssize_t i = 0; i < num_bins; ++i) {
|
||||
if (new_photon_counts[i] <= std::numeric_limits<double>::epsilon() &&
|
||||
store_nonzero_bins) {
|
||||
continue;
|
||||
}
|
||||
|
||||
output_file << std::floor(mythen_detector->min_angle() /
|
||||
histogram_bin_width) *
|
||||
histogram_bin_width +
|
||||
i * histogram_bin_width
|
||||
<< " " << new_photon_counts[i] << " "
|
||||
<< new_photon_count_errors[i] << std::endl;
|
||||
}
|
||||
output_file.close();
|
||||
}
|
||||
|
||||
} // namespace aare
|
||||
@ -1,265 +0,0 @@
|
||||
/************************************************
|
||||
* @file AngleCalibration.test.cpp
|
||||
* @short test case for angle calibration class
|
||||
***********************************************/
|
||||
|
||||
#include "aare/AngleCalibration.hpp"
|
||||
#include "aare/CustomFiles.hpp"
|
||||
#include "aare/FlatField.hpp"
|
||||
|
||||
#include <filesystem>
|
||||
|
||||
#include "test_config.hpp"
|
||||
|
||||
#include <iomanip>
|
||||
#include <type_traits>
|
||||
|
||||
#include <catch2/catch_all.hpp>
|
||||
#include <catch2/catch_test_macros.hpp>
|
||||
#include <catch2/matchers/catch_matchers_floating_point.hpp>
|
||||
|
||||
using namespace aare;
|
||||
|
||||
TEST_CASE("read initial angle calibration file",
|
||||
"[anglecalibration] [.files]") {
|
||||
|
||||
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
|
||||
std::make_shared<MythenDetectorSpecifications>();
|
||||
|
||||
AngleCalibration anglecalibration(mythen_detector_ptr,
|
||||
std::shared_ptr<FlatField>{},
|
||||
std::shared_ptr<MythenFileReader>{});
|
||||
|
||||
std::string filename = test_data_path() / "AngleCalibration_Test_Data" /
|
||||
"Angcal_2E_Feb2023_P29.off";
|
||||
|
||||
REQUIRE(std::filesystem::exists(filename));
|
||||
|
||||
anglecalibration.read_initial_calibration_from_file(filename);
|
||||
|
||||
auto centers = anglecalibration.get_centers();
|
||||
auto conversions = anglecalibration.get_conversions();
|
||||
auto offsets = anglecalibration.get_offsets();
|
||||
|
||||
std::cout.precision(17);
|
||||
|
||||
CHECK(centers.size() == 48);
|
||||
CHECK(conversions.size() == 48);
|
||||
CHECK(offsets.size() == 48);
|
||||
|
||||
CHECK(centers[9] == Catch::Approx(660.342326));
|
||||
CHECK(offsets[47] == Catch::Approx(5.8053312));
|
||||
CHECK(conversions[27] == Catch::Approx(-0.6581179125e-4));
|
||||
}
|
||||
|
||||
TEST_CASE("read bad channels",
|
||||
"[anglecalibration][mythenspecifications][.files]") {
|
||||
|
||||
MythenDetectorSpecifications mythen_detector;
|
||||
|
||||
std::string bad_channels_filename = test_data_path() /
|
||||
"AngleCalibration_Test_Data" /
|
||||
"bc2023_003_RING.chans";
|
||||
|
||||
REQUIRE(std::filesystem::exists(bad_channels_filename));
|
||||
|
||||
mythen_detector.read_bad_channels_from_file(bad_channels_filename);
|
||||
|
||||
CHECK(mythen_detector.get_bad_channels().size() == 61440);
|
||||
|
||||
CHECK(mythen_detector.get_bad_channels()[61437] == true);
|
||||
CHECK(std::all_of(mythen_detector.get_bad_channels().begin() + 30720,
|
||||
mythen_detector.get_bad_channels().begin() + 61439,
|
||||
[](const bool element) { return element; }));
|
||||
}
|
||||
|
||||
TEST_CASE("read unconnected modules",
|
||||
"[anglecalibration][mythenspecifications][.files]") {
|
||||
|
||||
MythenDetectorSpecifications mythen_detector;
|
||||
|
||||
std::string unconnected_modules_filename =
|
||||
test_data_path() / "AngleCalibration_Test_Data" / "ModOut.txt";
|
||||
|
||||
REQUIRE(std::filesystem::exists(unconnected_modules_filename));
|
||||
|
||||
mythen_detector.read_unconnected_modules_from_file(
|
||||
unconnected_modules_filename);
|
||||
|
||||
CHECK(mythen_detector.get_connected_modules().size() == 48);
|
||||
|
||||
CHECK(std::all_of(mythen_detector.get_connected_modules().begin(),
|
||||
mythen_detector.get_connected_modules().end(),
|
||||
[](const bool element) { return element; }));
|
||||
}
|
||||
|
||||
TEST_CASE("read flatfield", "[anglecalibration][flatfield][.files]") {
|
||||
|
||||
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
|
||||
std::make_shared<MythenDetectorSpecifications>();
|
||||
|
||||
FlatField flatfield(mythen_detector_ptr);
|
||||
|
||||
std::string flatfield_filename =
|
||||
test_data_path() / "AngleCalibration_Test_Data" /
|
||||
"Flatfield_E22p0keV_T11000eV_up_48M_a_LONG_Feb2023_open_WS_SUMC.raw";
|
||||
|
||||
REQUIRE(std::filesystem::exists(flatfield_filename));
|
||||
|
||||
flatfield.read_flatfield_from_file<CustomMythenFile>(flatfield_filename);
|
||||
|
||||
auto flatfield_data = flatfield.get_flatfield();
|
||||
|
||||
CHECK(flatfield_data.size() == 61440);
|
||||
|
||||
CHECK(flatfield_data[0] == 0);
|
||||
CHECK(flatfield_data[21] == 4234186);
|
||||
}
|
||||
|
||||
TEST_CASE("create flatfield", "[anglecalibration], [flatfield], [.files]") {
|
||||
|
||||
ssize_t n_modules = 2;
|
||||
double exposure_time = 5.0;
|
||||
size_t n_counters = 1;
|
||||
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
|
||||
std::make_shared<MythenDetectorSpecifications>(n_modules, exposure_time,
|
||||
n_counters);
|
||||
|
||||
FlatField flatfield(mythen_detector_ptr);
|
||||
|
||||
auto data_path = test_data_path() / "AngleCalibration_Test_Data" /
|
||||
"Flatfieldacquisition";
|
||||
|
||||
REQUIRE(std::filesystem::exists(data_path));
|
||||
|
||||
CHECK_NOTHROW(flatfield.create_flatfield_from_rawfilesystem(data_path));
|
||||
|
||||
auto flatfield_data = flatfield.get_flatfield();
|
||||
|
||||
CHECK(flatfield_data.size() == n_modules * 1280);
|
||||
|
||||
CHECK(flatfield_data[0] == 0);
|
||||
CHECK(flatfield_data[21] ==
|
||||
2 * 3 *
|
||||
21); // virtual data 2 angles, 3 frames - 21 as increasing numbers
|
||||
}
|
||||
|
||||
TEST_CASE("compare result with python code", "[anglecalibration] [.files]") {
|
||||
|
||||
auto fpath = test_data_path() / "AngleCalibration_Test_Data";
|
||||
|
||||
REQUIRE(std::filesystem::exists(fpath));
|
||||
|
||||
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
|
||||
std::make_shared<MythenDetectorSpecifications>();
|
||||
|
||||
std::string bad_channels_filename = fpath / "bc2023_003_RING.chans";
|
||||
|
||||
REQUIRE(std::filesystem::exists(bad_channels_filename));
|
||||
|
||||
mythen_detector_ptr->read_bad_channels_from_file(bad_channels_filename);
|
||||
|
||||
std::string unconnected_modules_filename = fpath / "ModOut.txt";
|
||||
|
||||
REQUIRE(std::filesystem::exists(unconnected_modules_filename));
|
||||
|
||||
mythen_detector_ptr->read_unconnected_modules_from_file(
|
||||
unconnected_modules_filename);
|
||||
|
||||
std::shared_ptr<FlatField> flat_field_ptr =
|
||||
std::make_shared<FlatField>(mythen_detector_ptr);
|
||||
|
||||
std::string flatfield_filename =
|
||||
fpath /
|
||||
"Flatfield_E22p0keV_T11000eV_up_48M_a_LONG_Feb2023_open_WS_SUMC.raw";
|
||||
|
||||
REQUIRE(std::filesystem::exists(flatfield_filename));
|
||||
|
||||
flat_field_ptr->read_flatfield_from_file<CustomMythenFile>(
|
||||
flatfield_filename);
|
||||
|
||||
std::shared_ptr<MythenFileReader> mythen_file_reader_ptr =
|
||||
std::make_shared<MythenFileReader>(fpath,
|
||||
"ang1up_22keV_LaB60p3mm_48M_a_0");
|
||||
|
||||
AngleCalibration anglecalibration(mythen_detector_ptr, flat_field_ptr,
|
||||
mythen_file_reader_ptr);
|
||||
|
||||
std::string initial_angles_filename = fpath / "Angcal_2E_Feb2023_P29.off";
|
||||
|
||||
REQUIRE(std::filesystem::exists(initial_angles_filename));
|
||||
|
||||
anglecalibration.read_initial_calibration_from_file(
|
||||
initial_angles_filename);
|
||||
|
||||
anglecalibration.calculate_fixed_bin_angle_width_histogram(320, 340);
|
||||
|
||||
// anglecalibration.write_to_file("cpp_new_photon_counts.xye");
|
||||
|
||||
auto expected_filename_photons =
|
||||
test_data_path() / "AngleCalibration_Test_Data" / "new_photons.bin";
|
||||
|
||||
REQUIRE(std::filesystem::exists(expected_filename_photons));
|
||||
|
||||
auto expected_filename_errors =
|
||||
test_data_path() / "AngleCalibration_Test_Data" / "new_errors.bin";
|
||||
|
||||
REQUIRE(std::filesystem::exists(expected_filename_errors));
|
||||
|
||||
ssize_t new_num_bins = anglecalibration.get_new_num_bins();
|
||||
|
||||
auto python_output_errors = load<double, 1>(
|
||||
expected_filename_errors, std::array<ssize_t, 1>{new_num_bins});
|
||||
|
||||
auto python_output_photons = load<double, 1>(
|
||||
expected_filename_photons, std::array<ssize_t, 1>{new_num_bins});
|
||||
|
||||
CHECK(anglecalibration.get_new_photon_counts().equals(
|
||||
python_output_photons.view(),
|
||||
1e-8)); // not sure about precision does not exactly match to all
|
||||
// decimal digits
|
||||
|
||||
CHECK(anglecalibration.get_new_statistical_errors().equals(
|
||||
python_output_errors.view(),
|
||||
1e-8)); //
|
||||
}
|
||||
|
||||
TEST_CASE("check conversion from DG to EE parameters", "[anglecalibration]") {
|
||||
|
||||
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
|
||||
std::make_shared<MythenDetectorSpecifications>();
|
||||
|
||||
AngleCalibration anglecalibration(mythen_detector_ptr,
|
||||
std::shared_ptr<FlatField>{},
|
||||
std::shared_ptr<MythenFileReader>{});
|
||||
|
||||
// DG test parameters
|
||||
const double center = 642.197591224993;
|
||||
const double conversion = 0.657694036246975e-4;
|
||||
const double offset = 5.004892881251670;
|
||||
const ssize_t local_strip_index = 1;
|
||||
|
||||
double diffraction_angle_DG_param =
|
||||
anglecalibration.diffraction_angle_from_DG_parameters(
|
||||
center, conversion, offset, local_strip_index);
|
||||
|
||||
auto [distance_center, normal_distance, angle] =
|
||||
anglecalibration.convert_to_EE_parameters(center, conversion, offset);
|
||||
|
||||
double diffraction_angle_EE_param =
|
||||
anglecalibration.diffraction_angle_from_EE_parameters(
|
||||
distance_center, normal_distance, angle, local_strip_index);
|
||||
|
||||
CHECK(diffraction_angle_EE_param ==
|
||||
Catch::Approx(diffraction_angle_DG_param));
|
||||
|
||||
double strip_width_DG_param =
|
||||
anglecalibration.angular_strip_width_from_DG_parameters(
|
||||
center, conversion, offset, local_strip_index);
|
||||
|
||||
double strip_width_EE_param =
|
||||
anglecalibration.angular_strip_width_from_EE_parameters(
|
||||
distance_center, normal_distance, angle, local_strip_index);
|
||||
|
||||
CHECK(strip_width_DG_param == Catch::Approx(strip_width_EE_param));
|
||||
}
|
||||
@ -1,89 +0,0 @@
|
||||
#include "aare/CustomFiles.hpp"
|
||||
|
||||
namespace aare {
|
||||
|
||||
CustomMythenFile::CustomMythenFile(const std::string &filename,
|
||||
const ssize_t rows, const ssize_t cols,
|
||||
const std::string &mode)
|
||||
: m_filename(filename), m_rows(rows), m_cols(cols) {
|
||||
|
||||
m_mode = mode;
|
||||
if (m_mode == "r") {
|
||||
try {
|
||||
m_file.open(m_filename, std::ios_base::in);
|
||||
if (!m_file.good()) {
|
||||
throw std::logic_error("file does not exist");
|
||||
}
|
||||
} catch (std::exception &e) {
|
||||
|
||||
std::cerr << "Error: " << e.what()
|
||||
<< std::endl; // TODO replace with log
|
||||
}
|
||||
|
||||
} else {
|
||||
throw std::runtime_error(LOCATION +
|
||||
"Unsupported mode. Can only read RawFiles.");
|
||||
}
|
||||
}
|
||||
|
||||
CustomMythenFile::~CustomMythenFile() { m_file.close(); }
|
||||
|
||||
Frame CustomMythenFile::read_frame() {
|
||||
auto f = Frame(m_rows, m_cols, m_dtype);
|
||||
uint32_t *frame_buffer = reinterpret_cast<uint32_t *>(f.data());
|
||||
uint32_t strip_index, photon_count;
|
||||
while (m_file >> strip_index >> photon_count) {
|
||||
*frame_buffer = photon_count;
|
||||
++frame_buffer;
|
||||
}
|
||||
return f;
|
||||
}
|
||||
|
||||
void CustomMythenFile::read_into(std::byte *image_buf) {
|
||||
uint32_t strip_index, photon_count;
|
||||
while (m_file >> strip_index >> photon_count) {
|
||||
std::memcpy(image_buf, &photon_count, sizeof(photon_count));
|
||||
image_buf += sizeof(photon_count);
|
||||
}
|
||||
}
|
||||
|
||||
size_t CustomMythenFile::bytes_per_frame() {
|
||||
return m_num_strips * m_dtype.bytes(); // TODO do i want m_counts?
|
||||
}
|
||||
|
||||
Frame CustomMythenFile::read_frame(size_t frame_number) {
|
||||
return read_frame(); // maybe give count as frame_number
|
||||
}
|
||||
|
||||
std::vector<Frame> CustomMythenFile::read_n(size_t n_frames) {
|
||||
std::vector<Frame> vec;
|
||||
vec.reserve(1);
|
||||
vec.push_back(read_frame());
|
||||
return vec; // std::vector<Frame>{read_frame()};
|
||||
}
|
||||
|
||||
void CustomMythenFile::read_into(std::byte *image_buf, size_t n_frames) {
|
||||
read_into(image_buf);
|
||||
}
|
||||
|
||||
size_t CustomMythenFile::frame_number(size_t frame_index) { return 1; }
|
||||
|
||||
size_t CustomMythenFile::pixels_per_frame() { return m_rows * m_cols; }
|
||||
|
||||
void CustomMythenFile::seek(size_t frame_number) {}
|
||||
|
||||
size_t CustomMythenFile::tell() { return 1; }
|
||||
|
||||
size_t CustomMythenFile::total_frames() const { return 1; }
|
||||
|
||||
size_t CustomMythenFile::rows() const { return m_rows; }
|
||||
|
||||
size_t CustomMythenFile::cols() const { return m_cols; }
|
||||
|
||||
size_t CustomMythenFile::bitdepth() const { return m_dtype.bitdepth(); }
|
||||
|
||||
DetectorType CustomMythenFile::detector_type() const { return det_type; }
|
||||
|
||||
const Dtype CustomMythenFile::m_dtype(Dtype::TypeIndex::UINT32);
|
||||
|
||||
} // namespace aare
|
||||
@ -1,110 +0,0 @@
|
||||
/************************************************
|
||||
* @file Hdf5FileReader.test.cpp
|
||||
* @short test case for reading hdf5 files
|
||||
***********************************************/
|
||||
|
||||
#include <filesystem>
|
||||
|
||||
#include "test_config.hpp"
|
||||
|
||||
#include <H5Cpp.h>
|
||||
|
||||
#include "aare/Hdf5FileReader.hpp"
|
||||
#include "aare/NDArray.hpp"
|
||||
|
||||
#include <catch2/catch_all.hpp>
|
||||
#include <catch2/catch_test_macros.hpp>
|
||||
#include <catch2/matchers/catch_matchers_floating_point.hpp>
|
||||
|
||||
using namespace aare;
|
||||
|
||||
TEST_CASE("read hdf5 file", "[.hdf5file][.files]") {
|
||||
|
||||
// TODO generalize datasetpath
|
||||
std::string filename = test_data_path() / "AngleCalibration_Test_Data" /
|
||||
"ang1up_22keV_LaB60p3mm_48M_a_0320.h5";
|
||||
|
||||
REQUIRE(std::filesystem::exists(filename));
|
||||
|
||||
HDF5FileReader file_reader;
|
||||
|
||||
file_reader.open_file(filename);
|
||||
|
||||
auto dataset = file_reader.get_dataset("/entry/data/data");
|
||||
|
||||
auto shape = dataset.get_shape();
|
||||
|
||||
CHECK(shape[0] == 61440);
|
||||
|
||||
auto type = dataset.get_datatype();
|
||||
|
||||
const std::type_info *type_info = dataset.get_cpp_type();
|
||||
|
||||
CHECK(*type_info == typeid(uint32_t));
|
||||
|
||||
SECTION("read dataset into NDArray") {
|
||||
|
||||
NDArray<uint32_t, 1> dataset_array =
|
||||
dataset.store_as_ndarray<uint32_t, 1>();
|
||||
|
||||
CHECK(dataset_array(0) == 866);
|
||||
CHECK(dataset_array(61439) == 1436);
|
||||
}
|
||||
|
||||
SECTION("read dataset into Frame") {
|
||||
Frame frame = dataset.store_as_frame();
|
||||
CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 0))) == 866);
|
||||
CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 61439))) ==
|
||||
1436);
|
||||
}
|
||||
SECTION("read subset of dataset") {
|
||||
Frame frame(1, 10, Dtype(typeid(uint32_t)));
|
||||
|
||||
Subset subset{std::vector<hsize_t>{10}, std::vector<hsize_t>{10}};
|
||||
|
||||
dataset.read_into_buffer(frame.data(), subset);
|
||||
|
||||
CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 0))) == 664);
|
||||
CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 9))) == 654);
|
||||
}
|
||||
/*
|
||||
SECTION("read scalar") {
|
||||
}
|
||||
*/
|
||||
}
|
||||
|
||||
TEST_CASE("test datatypes", "[.hdf5file]") {
|
||||
|
||||
auto [dtype, expected_type_info] = GENERATE(
|
||||
std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT), &typeid(int)),
|
||||
std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT8),
|
||||
&typeid(int8_t)),
|
||||
std::make_tuple(H5::DataType(H5::PredType::NATIVE_UINT16),
|
||||
&typeid(uint16_t)),
|
||||
std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT16),
|
||||
&typeid(int16_t)),
|
||||
std::make_tuple(H5::DataType(H5::PredType::STD_U32LE),
|
||||
&typeid(uint32_t)),
|
||||
std::make_tuple(H5::DataType(H5::PredType::STD_I32LE),
|
||||
&typeid(int32_t)),
|
||||
std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT32),
|
||||
&typeid(int32_t)),
|
||||
std::make_tuple(H5::DataType(H5::PredType::IEEE_F64LE),
|
||||
&typeid(double)),
|
||||
std::make_tuple(H5::DataType(H5::PredType::IEEE_F32LE), &typeid(float)),
|
||||
std::make_tuple(H5::DataType(H5::PredType::NATIVE_FLOAT),
|
||||
&typeid(float)),
|
||||
std::make_tuple(H5::DataType(H5::PredType::NATIVE_DOUBLE),
|
||||
&typeid(double)),
|
||||
std::make_tuple(H5::DataType(H5::PredType::NATIVE_CHAR),
|
||||
&typeid(int8_t)));
|
||||
|
||||
const std::type_info &type_info = deduce_cpp_type(dtype);
|
||||
|
||||
CHECK(type_info == *expected_type_info);
|
||||
|
||||
// TODO: handle bit swapping
|
||||
REQUIRE_THROWS(deduce_cpp_type(
|
||||
H5::DataType(H5::PredType::IEEE_F32BE))); // does not convert from big
|
||||
// to little endian
|
||||
}
|
||||
@ -1,33 +0,0 @@
|
||||
/************************************************
|
||||
* @file MythenFileReader.test.cpp
|
||||
* @short test case for angle calibration class
|
||||
***********************************************/
|
||||
|
||||
#include "aare/MythenFileReader.hpp"
|
||||
|
||||
#include <filesystem>
|
||||
|
||||
#include "test_config.hpp"
|
||||
|
||||
#include <catch2/catch_all.hpp>
|
||||
#include <catch2/catch_test_macros.hpp>
|
||||
#include <catch2/matchers/catch_matchers_floating_point.hpp>
|
||||
|
||||
using namespace aare;
|
||||
|
||||
TEST_CASE("test mythenfile_reader", "[.mythenfilereader][.files]") {
|
||||
|
||||
auto fpath = test_data_path() / "AngleCalibration_Test_Data";
|
||||
|
||||
REQUIRE(std::filesystem::exists(fpath));
|
||||
|
||||
MythenFileReader file_reader(fpath, "ang1up_22keV_LaB60p3mm_48M_a_0");
|
||||
|
||||
auto frame = file_reader.read_frame(320);
|
||||
|
||||
CHECK(frame.detector_angle == 0.99955);
|
||||
|
||||
CHECK(frame.channel_mask == std::array<uint8_t, 3>{0, 0, 1});
|
||||
|
||||
CHECK(frame.photon_counts.size() == 61440);
|
||||
}
|
||||
Reference in New Issue
Block a user