diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 862094d5..c51b8a44 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -2,8 +2,6 @@ ADD_LIBRARY(ImageAnalysis STATIC CrystalLattice.cpp CrystalLattice.h IndexerWrapper.cpp IndexerWrapper.h GPUImageAnalysis.h - RadialIntegration.cpp RadialIntegration.h - RadialIntegrationMappingFPGA.cpp RadialIntegrationMappingFPGA.h RadialIntegrationMapping.cpp RadialIntegrationMapping.h StrongPixelSet.cpp StrongPixelSet.h GPUImageAnalysis.cpp RadialIntegrationProfile.cpp RadialIntegrationProfile.h PredictSpotsOnDetector.h) diff --git a/image_analysis/GPUImageAnalysis.cpp b/image_analysis/GPUImageAnalysis.cpp index e336b79f..ae4487a9 100644 --- a/image_analysis/GPUImageAnalysis.cpp +++ b/image_analysis/GPUImageAnalysis.cpp @@ -12,15 +12,6 @@ GPUImageAnalysis::GPUImageAnalysis(int32_t in_xpixels, int32_t in_ypixels, const xpixels(in_xpixels), ypixels(in_ypixels), gpu_out(nullptr), rad_integration_nbins(0) { } -GPUImageAnalysis::GPUImageAnalysis(int32_t xpixels, int32_t ypixels, const std::vector &mask, - const std::vector &rad_int_mapping, uint16_t rad_int_nbins) - : GPUImageAnalysis(xpixels, ypixels, mask) {} - -GPUImageAnalysis::GPUImageAnalysis(int32_t xpixels, int32_t ypixels, const std::vector &mask, - const RadialIntegrationMapping& mapping) - : GPUImageAnalysis(xpixels, ypixels, mask, mapping.GetPixelToBinMapping(), - mapping.GetBinNumber()) {} - GPUImageAnalysis::~GPUImageAnalysis() {} void GPUImageAnalysis::SetInputBuffer(void *ptr) {} @@ -43,18 +34,4 @@ void GPUImageAnalysis::UnregisterBuffer() {} void GPUImageAnalysis::LoadDataToGPU(bool apply_pixel_mask_on_gpu) {} -void GPUImageAnalysis::RunRadialIntegration() {} - -void GPUImageAnalysis::GetRadialIntegrationProfile(std::vector &result) {} - -std::vector GPUImageAnalysis::GetRadialIntegrationSum() const { - return {}; -} - -std::vector GPUImageAnalysis::GetRadialIntegrationCount() const { - return {}; -} - -void GPUImageAnalysis::LoadRadialIntegrationCorr(const std::vector& v) {} - #endif diff --git a/image_analysis/GPUImageAnalysis.cu b/image_analysis/GPUImageAnalysis.cu index d2dabedf..0157b90e 100644 --- a/image_analysis/GPUImageAnalysis.cu +++ b/image_analysis/GPUImageAnalysis.cu @@ -178,36 +178,6 @@ __global__ void analyze_pixel(const int16_t *in, uint8_t *out, const spot_parame } while (back < rmax); } -__global__ void gpu_radial_integration(const int16_t *image, const uint16_t *rad_integration_mapping, - float *corr, float *bin_sum, float *bin_count, uint32_t npixel, - uint16_t nbins) { - - extern __shared__ float shared_mem_fp[]; - float* shared_sum = shared_mem_fp; // shared buffer [nbins] - float* shared_count = &shared_sum[nbins]; // shared buffer [nbins] - - uint32_t idx = blockDim.x*blockIdx.x + threadIdx.x; - - for (uint32_t i = threadIdx.x; i < 2 * nbins; i += blockDim.x) - shared_mem_fp[i] = 0; - __syncthreads(); - - for (uint32_t i = idx; i < npixel; i += blockDim.x * gridDim.x) { - uint16_t bin = rad_integration_mapping[i]; - float value = static_cast(image[i]) * corr[i]; - if ((image[i] > INT16_MIN + 4) && (image[i] < INT16_MAX - 4) && (bin < nbins)) { - atomicAdd(&shared_sum[bin], value); - atomicAdd(&shared_count[bin], 1.0f); - } - } - __syncthreads(); - - for (uint32_t i = threadIdx.x; i < nbins; i += blockDim.x) { - atomicAdd(&bin_sum[i], shared_sum[i]); - atomicAdd(&bin_count[i], shared_count[i]); - } -} - __global__ void spot_finder_reduce(uint32_t *output, uint32_t* counter, const uint8_t* input, uint32_t npixel, uint32_t output_size) { uint32_t idx = blockDim.x*blockIdx.x + threadIdx.x; @@ -229,7 +199,7 @@ __global__ void apply_pixel_mask(int16_t *image, const uint8_t *mask, uint32_t n } GPUImageAnalysis::GPUImageAnalysis(int32_t in_xpixels, int32_t in_ypixels, const std::vector &mask) : - xpixels(in_xpixels), ypixels(in_ypixels), gpu_out(nullptr), rad_integration_nbins(0), numberOfSMs(1) { + xpixels(in_xpixels), ypixels(in_ypixels), gpu_out(nullptr), numberOfSMs(1) { if (get_gpu_count() == 0) throw JFJochException(JFJochExceptionCategory::GPUCUDAError, "No CUDA devices found"); @@ -260,49 +230,10 @@ GPUImageAnalysis::GPUImageAnalysis(int32_t in_xpixels, int32_t in_ypixels, const cudaMemcpy(gpu_mask, mask.data(), xpixels*ypixels, cudaMemcpyHostToDevice); } -GPUImageAnalysis::GPUImageAnalysis(int32_t xpixels, int32_t ypixels, const std::vector &mask, - const std::vector &rad_int_mapping, uint16_t rad_int_nbins) - : GPUImageAnalysis(xpixels, ypixels, mask) { - rad_integration_nbins = rad_int_nbins; - - if (rad_int_mapping.size() != xpixels * ypixels) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Mismatch in rad. int. mapping size"); - - if (rad_integration_nbins > 0) { - cuda_err(cudaMalloc(&gpu_rad_integration_bin_map, xpixels * ypixels * sizeof(int16_t))); - cuda_err(cudaMalloc(&gpu_rad_integration_corr, xpixels * ypixels * sizeof(float))); - - cuda_err(cudaMalloc(&gpu_rad_integration_count, rad_integration_nbins * sizeof(float))); - cuda_err(cudaMalloc(&gpu_rad_integration_sum, rad_integration_nbins * sizeof(float))); - cuda_err(cudaHostAlloc(&host_rad_integration_count, rad_integration_nbins * sizeof(float), cudaHostAllocPortable)); - cuda_err(cudaHostAlloc(&host_rad_integration_sum, rad_integration_nbins * sizeof(float), cudaHostAllocPortable)); - - cudaMemcpy(gpu_rad_integration_bin_map, rad_int_mapping.data(), - xpixels*ypixels * sizeof(uint16_t), cudaMemcpyHostToDevice); - - std::vector corr(xpixels * ypixels, 1.0f); - cudaMemcpy(gpu_rad_integration_corr, corr.data(), - xpixels * ypixels * sizeof(float), cudaMemcpyHostToDevice); - } -} - -GPUImageAnalysis::GPUImageAnalysis(int32_t xpixels, int32_t ypixels, const std::vector &mask, - const RadialIntegrationMapping& mapping) - : GPUImageAnalysis(xpixels, ypixels, mask, mapping.GetPixelToBinMapping(),mapping.GetBinNumber()) {} - GPUImageAnalysis::~GPUImageAnalysis() { cudaStreamDestroy(cudastream->v); delete(cudastream); - if (rad_integration_nbins > 0) { - cudaFree(gpu_rad_integration_bin_map); - cudaFree(gpu_rad_integration_count); - cudaFree(gpu_rad_integration_sum); - cudaFree(gpu_rad_integration_corr); - cudaFreeHost(host_rad_integration_count); - cudaFreeHost(host_rad_integration_sum); - } - cudaFreeHost(host_out_reduced); cudaFreeHost(host_out_reduced_counter); cudaFree(gpu_mask); @@ -419,70 +350,3 @@ void GPUImageAnalysis::LoadDataToGPU(bool apply_pixel_mask_on_gpu) { } } -void GPUImageAnalysis::LoadRadialIntegrationCorr(const std::vector& v) { - if (rad_integration_nbins == 0) - throw JFJochException(JFJochExceptionCategory::SpotFinderError, "Radial integration not initialized"); - if (v.size() != xpixels * ypixels) - throw JFJochException(JFJochExceptionCategory::SpotFinderError, "Mismatch in correction input size"); - - cudaMemcpy(gpu_rad_integration_corr, v.data(), xpixels * ypixels * sizeof(float), cudaMemcpyHostToDevice); -} - -void GPUImageAnalysis::RunRadialIntegration() { - if (rad_integration_nbins == 0) - throw JFJochException(JFJochExceptionCategory::SpotFinderError, "Radial integration not initialized"); - - cuda_err(cudaMemsetAsync(gpu_rad_integration_sum, 0, rad_integration_nbins * sizeof(int32_t), cudastream->v)); - cuda_err(cudaMemsetAsync(gpu_rad_integration_count, 0, rad_integration_nbins * sizeof(int32_t), cudastream->v)); - - gpu_radial_integration<<v>>>( - gpu_in, gpu_rad_integration_bin_map, gpu_rad_integration_corr, - gpu_rad_integration_sum, gpu_rad_integration_count, xpixels*ypixels, - rad_integration_nbins); - - cuda_err(cudaMemcpyAsync(host_rad_integration_count, gpu_rad_integration_count, - rad_integration_nbins * sizeof(int32_t), - cudaMemcpyDeviceToHost, cudastream->v)); - cuda_err(cudaMemcpyAsync(host_rad_integration_sum, gpu_rad_integration_sum, - rad_integration_nbins * sizeof(int32_t), - cudaMemcpyDeviceToHost, cudastream->v)); -} - -void GPUImageAnalysis::GetRadialIntegrationProfile(std::vector &result) { - if (rad_integration_nbins == 0) - throw JFJochException(JFJochExceptionCategory::SpotFinderError, "Radial integration not initialized"); - - cuda_err(cudaStreamSynchronize(cudastream->v)); - - result.resize(rad_integration_nbins); - - for (int i = 0; i < rad_integration_nbins; i++) { - if (host_rad_integration_count[i] > 0) - result[i] = static_cast(host_rad_integration_sum[i]) - / static_cast(host_rad_integration_count[i]); - else - result[i] = 0; - } -} - -std::vector GPUImageAnalysis::GetRadialIntegrationSum() const { - if (rad_integration_nbins == 0) - throw JFJochException(JFJochExceptionCategory::SpotFinderError, "Radial integration not initialized"); - - cuda_err(cudaStreamSynchronize(cudastream->v)); - - std::vector out(rad_integration_nbins); - memcpy(out.data(), host_rad_integration_sum, rad_integration_nbins * sizeof(int32_t)); - return out; -} - -std::vector GPUImageAnalysis::GetRadialIntegrationCount() const { - if (rad_integration_nbins == 0) - throw JFJochException(JFJochExceptionCategory::SpotFinderError, "Radial integration not initialized"); - - cuda_err(cudaStreamSynchronize(cudastream->v)); - - std::vector out(rad_integration_nbins); - memcpy(out.data(), host_rad_integration_count, rad_integration_nbins * sizeof(int32_t)); - return out; -} diff --git a/image_analysis/GPUImageAnalysis.h b/image_analysis/GPUImageAnalysis.h index aea1e19d..5b4d0bb1 100644 --- a/image_analysis/GPUImageAnalysis.h +++ b/image_analysis/GPUImageAnalysis.h @@ -6,7 +6,6 @@ #include #include "StrongPixelSet.h" -#include "RadialIntegrationMapping.h" struct spot_parameters { float strong_pixel_threshold2; @@ -36,14 +35,6 @@ class GPUImageAnalysis { uint32_t *host_out_reduced = nullptr; uint32_t *host_out_reduced_counter = nullptr; - uint16_t rad_integration_nbins; - uint16_t *gpu_rad_integration_bin_map = nullptr; - float *gpu_rad_integration_sum = nullptr; - float *gpu_rad_integration_count = nullptr; - float *gpu_rad_integration_corr = nullptr; - - float *host_rad_integration_sum = nullptr, *host_rad_integration_count = nullptr; - int numberOfSMs; const int numberOfCudaThreads = 128; // #threads per block that works well for Nvidia T4 const int numberOfWaves = 40; // #waves that works well for Nvidia T4 @@ -52,11 +43,6 @@ class GPUImageAnalysis { public: GPUImageAnalysis(int32_t xpixels, int32_t ypixels, const std::vector &mask); - GPUImageAnalysis(int32_t xpixels, int32_t ypixels, const std::vector &mask, - const std::vector &rad_int_mapping, uint16_t rad_int_nbins); - GPUImageAnalysis(int32_t xpixels, int32_t ypixels, const std::vector &mask, - const RadialIntegrationMapping& mapping); - ~GPUImageAnalysis(); void SetInputBuffer(void *host_in); @@ -67,12 +53,6 @@ public: void GetSpotFinderResults(const DiffractionExperiment &experiment, const JFJochProtoBuf::DataProcessingSettings &settings, std::vector &vec); - void LoadRadialIntegrationCorr(const std::vector& v); - void RunRadialIntegration(); - void GetRadialIntegrationProfile(std::vector &result); - [[nodiscard]] std::vector GetRadialIntegrationSum() const; - [[nodiscard]] std::vector GetRadialIntegrationCount() const; - static bool GPUPresent(); void RegisterBuffer(); void UnregisterBuffer(); diff --git a/image_analysis/RadialIntegration.cpp b/image_analysis/RadialIntegration.cpp deleted file mode 100644 index 27b56b0a..00000000 --- a/image_analysis/RadialIntegration.cpp +++ /dev/null @@ -1,98 +0,0 @@ -// Copyright (2019-2023) Paul Scherrer Institute - -#include "RadialIntegration.h" -#include "../common/JFJochException.h" - -RadialIntegration::RadialIntegration(const std::vector& in_mapping, uint32_t in_nbins) : - pixel_to_bin(in_mapping), nbins(in_nbins), sum(in_nbins, 0), count(in_nbins, 0) -{ - - coeff = (float *) std::aligned_alloc(64, in_mapping.size() * sizeof(float)); - if (coeff == nullptr) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Memory allocation error"); - - for (int i = 0; i < in_mapping.size(); i++) - coeff[i] = 1.0f; -} - -RadialIntegration::RadialIntegration(const RadialIntegrationMapping &mapping) : - RadialIntegration(mapping.GetPixelToBinMapping(), mapping.GetBinNumber()) -{} - -RadialIntegration::~RadialIntegration() { - std::free(coeff); -} - -void RadialIntegration::LoadRadialIntegrationCorr(const std::vector &v) { - if (v.size() != pixel_to_bin.size()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Mismatch in size of pixel-to-bin mapping and correction"); - memcpy(coeff, v.data(), pixel_to_bin.size() * sizeof(float)); -} - -void RadialIntegration::Clear() { - for (auto &i : sum) - i = 0; - - for (auto &i : count) - i = 0; -} - -void RadialIntegration::Process(const int16_t *__restrict data, size_t npixel) { - if (npixel != pixel_to_bin.size()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Mismatch in size of pixel-to-bin mapping and image"); - - const auto coeff_aligned = std::assume_aligned<64>(coeff); - - for (int i = 0; i < npixel; i++) { - auto value = data[i]; - if ((value > INT16_MIN + 4) && (value < INT16_MAX - 4)) { - auto bin = pixel_to_bin[i]; - if (bin < nbins) { - sum[bin] += coeff_aligned[i] * value; - count[bin] += 1.0f; - } - } - } - -} - -void RadialIntegration::ProcessOneImage(const int16_t *data, size_t npixel) { - Clear(); - Process(data, npixel); -} - -void RadialIntegration::GetResult(std::vector &result) const { - result.resize(nbins); - - for (int i = 0; i < nbins; i++) { - if (count[i] > 0) - result[i] = static_cast(sum[i]) / static_cast(count[i]); - else - result[i] = 0; - } -} - -const std::vector &RadialIntegration::GetSum() const { - return sum; -} - -const std::vector &RadialIntegration::GetCount() const { - return count; -} - -float RadialIntegration::GetRangeValue(uint32_t min_bin, uint32_t max_bin) { - float ret_sum = 0; - float ret_count = 0; - - for (int i = std::min(nbins,min_bin); i <= std::min(nbins-1,max_bin); i++) { - ret_sum += sum[i]; - ret_count += count[i]; - } - if (ret_count == 0) - return 0; - else - return static_cast(ret_sum) / static_cast(ret_count); -} \ No newline at end of file diff --git a/image_analysis/RadialIntegration.h b/image_analysis/RadialIntegration.h deleted file mode 100644 index 91d15eaf..00000000 --- a/image_analysis/RadialIntegration.h +++ /dev/null @@ -1,34 +0,0 @@ -// Copyright (2019-2023) Paul Scherrer Institute - -#ifndef JUNGFRAUJOCH_RADIALINTEGRATION_H -#define JUNGFRAUJOCH_RADIALINTEGRATION_H - -#include -#include -#include - -#include "../common/DiffractionExperiment.h" -#include "RadialIntegrationMapping.h" - -class RadialIntegration { - const std::vector pixel_to_bin; - float *coeff; - const uint32_t nbins; - std::vector sum; - std::vector count; -public: - RadialIntegration(const RadialIntegrationMapping& mapping); - RadialIntegration(const std::vector& mapping, uint32_t nbins); - ~RadialIntegration(); - void Clear(); - void Process(const int16_t *data, size_t npixel); - void ProcessOneImage(const int16_t *data, size_t npixel); // Process + Clear - void GetResult(std::vector &result) const; - void LoadRadialIntegrationCorr(const std::vector &v); - [[nodiscard]] float GetRangeValue(uint32_t min_bin, uint32_t max_bin); - [[nodiscard]] const std::vector& GetSum() const; - [[nodiscard]] const std::vector& GetCount() const; -}; - - -#endif //JUNGFRAUJOCH_RADIALINTEGRATION_H diff --git a/image_analysis/RadialIntegrationMapping.cpp b/image_analysis/RadialIntegrationMapping.cpp index 6da3f21a..2ea4daff 100644 --- a/image_analysis/RadialIntegrationMapping.cpp +++ b/image_analysis/RadialIntegrationMapping.cpp @@ -10,7 +10,7 @@ RadialIntegrationMapping::RadialIntegrationMapping(const DiffractionExperiment& low_q(experiment.GetLowQForRadialInt_recipA()), high_q(experiment.GetHighQForRadialInt_recipA()), q_spacing(experiment.GetQSpacingForRadialInt_recipA()), - pixel_to_bin(experiment.GetPixelsNum()), + pixel_to_bin_raw(experiment.GetModulesNum() * RAW_MODULE_SIZE), max_bin_number(0) { @@ -30,21 +30,20 @@ RadialIntegrationMapping::RadialIntegrationMapping(const DiffractionExperiment& throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot accommodate more than 65536 rings"); - for (int y = 0; y < experiment.GetYPixelsNum(); y++) { - for (int x = 0; x < experiment.GetXPixelsNum(); x++) { - int64_t pixel_number = y * experiment.GetXPixelsNum() + x; - + for (int m = 0; m < experiment.GetModulesNum(); m++) { + for (int pxl = 0; pxl < RAW_MODULE_SIZE; pxl++) { + auto [x,y] = RawToConvertedCoordinate(experiment, m, pxl); double pixel_q = 2 * M_PI / PxlToRes(experiment, x, y); if ((pixel_q < low_q) || (pixel_q >= high_q)) - pixel_to_bin[pixel_number] = UINT16_MAX; + pixel_to_bin_raw[m * RAW_MODULE_SIZE + pxl] = UINT16_MAX; else - pixel_to_bin[pixel_number] = std::floor((pixel_q - low_q) / q_spacing); + pixel_to_bin_raw[m * RAW_MODULE_SIZE + pxl] = std::floor((pixel_q - low_q) / q_spacing); } } // In principle max bin number is equal to std::floor((high_q - low_q) / q_spacing), but it might be lower than this // depending on detector distance and beam center position - for (const auto &i: pixel_to_bin) { + for (const auto &i: pixel_to_bin_raw) { if ((i != UINT16_MAX) && (i > max_bin_number)) max_bin_number = i; } @@ -62,8 +61,8 @@ uint16_t RadialIntegrationMapping::GetBinNumber() const { return max_bin_number + 1; } -const std::vector &RadialIntegrationMapping::GetPixelToBinMapping() const { - return pixel_to_bin; +const std::vector &RadialIntegrationMapping::GetPixelToBinMappingRaw() const { + return pixel_to_bin_raw; } const std::vector &RadialIntegrationMapping::GetBinToQ() const { diff --git a/image_analysis/RadialIntegrationMapping.h b/image_analysis/RadialIntegrationMapping.h index a2d1d7c3..df1a11e1 100644 --- a/image_analysis/RadialIntegrationMapping.h +++ b/image_analysis/RadialIntegrationMapping.h @@ -10,12 +10,12 @@ class RadialIntegrationMapping { const double low_q, high_q, q_spacing; std::vector bin_to_q; std::vector solid_angle_corr; - std::vector pixel_to_bin; + std::vector pixel_to_bin_raw; uint16_t max_bin_number; public: RadialIntegrationMapping(const DiffractionExperiment& experiment); [[nodiscard]] uint16_t GetBinNumber() const; - [[nodiscard]] const std::vector &GetPixelToBinMapping() const; + [[nodiscard]] const std::vector& GetPixelToBinMappingRaw() const; [[nodiscard]] const std::vector &GetBinToQ() const; [[nodiscard]] const std::vector &GetSolidAngleCorr() const; [[nodiscard]] double QToBin(double q) const; diff --git a/image_analysis/RadialIntegrationMappingFPGA.cpp b/image_analysis/RadialIntegrationMappingFPGA.cpp deleted file mode 100644 index 3a341677..00000000 --- a/image_analysis/RadialIntegrationMappingFPGA.cpp +++ /dev/null @@ -1,78 +0,0 @@ -// Copyright (2019-2023) Paul Scherrer Institute - -#include -#include "RadialIntegrationMappingFPGA.h" -#include "../common/RawToConvertedGeometry.h" -#include "../common/DiffractionGeometry.h" - -RadialIntegrationMappingFPGA::RadialIntegrationMappingFPGA(const DiffractionExperiment &experiment, uint16_t data_stream) : - low_q(experiment.GetLowQForRadialInt_recipA()), - high_q(experiment.GetHighQForRadialInt_recipA()), - q_spacing(experiment.GetQSpacingForRadialInt_recipA()), - pixel_to_bin(experiment.GetModulesNum(data_stream) * RAW_MODULE_SIZE), - max_bin_number(0) -{ - - if (q_spacing < 0.0) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Q-spacing has to be >= 0.0"); - - if (low_q >= high_q) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Low-Q must be smaller than high-Q"); - - if (low_q <= 0) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Low-Q must be > 0"); - - if (std::ceil((high_q-low_q) / q_spacing ) >= FPGA_INTEGRATION_BIN_COUNT) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Cannot accommodate more than " + std::to_string(FPGA_INTEGRATION_BIN_COUNT) + " rings"); - - for (int m = 0; m < experiment.GetModulesNum(data_stream); m++) { - for (int pxl = 0; pxl < RAW_MODULE_SIZE; pxl++) { - auto [x,y] = RawToConvertedCoordinate(experiment, m + experiment.GetFirstModuleOfDataStream(data_stream), pxl); - double pixel_q = 2 * M_PI / PxlToRes(experiment, x, y); - if ((pixel_q < low_q) || (pixel_q >= high_q)) - pixel_to_bin[m * RAW_MODULE_SIZE + pxl] = UINT16_MAX; - else - pixel_to_bin[m * RAW_MODULE_SIZE + pxl] = std::floor((pixel_q - low_q) / q_spacing); - } - } - - // In principle max bin number is equal to std::floor((high_q - low_q) / q_spacing), but it might be lower than this - // depending on detector distance and beam center position - for (const auto &i: pixel_to_bin) { - if ((i != UINT16_MAX) && (i > max_bin_number)) - max_bin_number = i; - } - - bin_to_q.resize(max_bin_number + 1); - for (int i = 0; i < max_bin_number + 1; i++) - bin_to_q[i] = static_cast(q_spacing * (i + 0.5) + low_q); - - solid_angle_corr.resize(max_bin_number + 1); - for (int i = 0; i < max_bin_number + 1; i++) - solid_angle_corr[i] = 1.0f / CalcRadIntSolidAngleCorr(experiment, bin_to_q[i]); - -} - -uint16_t RadialIntegrationMappingFPGA::GetBinNumber() const { - return max_bin_number + 1; -} - -const std::vector &RadialIntegrationMappingFPGA::GetPixelToBinMapping() const { - return pixel_to_bin; -} - -const std::vector &RadialIntegrationMappingFPGA::GetBinToQ() const { - return bin_to_q; -} - -const std::vector &RadialIntegrationMappingFPGA::GetSolidAngleCorr() const { - return solid_angle_corr; -} - -double RadialIntegrationMappingFPGA::QToBin(double q) const { - return std::min(static_cast(max_bin_number), std::max(0.0, (q - low_q) / q_spacing)); -} diff --git a/image_analysis/RadialIntegrationMappingFPGA.h b/image_analysis/RadialIntegrationMappingFPGA.h deleted file mode 100644 index 4957db8e..00000000 --- a/image_analysis/RadialIntegrationMappingFPGA.h +++ /dev/null @@ -1,25 +0,0 @@ -// Copyright (2019-2023) Paul Scherrer Institute - -#ifndef JUNGFRAUJOCH_RADIALINTEGRATIONMAPPINGFPGA_H -#define JUNGFRAUJOCH_RADIALINTEGRATIONMAPPINGFPGA_H - -#include -#include -#include "../common/DiffractionExperiment.h" - -class RadialIntegrationMappingFPGA { - const double low_q, high_q, q_spacing; - std::vector bin_to_q; - std::vector solid_angle_corr; - std::vector pixel_to_bin; - uint16_t max_bin_number; -public: - explicit RadialIntegrationMappingFPGA(const DiffractionExperiment& experiment, uint16_t data_stream); - [[nodiscard]] uint16_t GetBinNumber() const; - [[nodiscard]] const std::vector &GetPixelToBinMapping() const; - [[nodiscard]] const std::vector &GetBinToQ() const; - [[nodiscard]] const std::vector &GetSolidAngleCorr() const; - [[nodiscard]] double QToBin(double q) const; -}; - -#endif //JUNGFRAUJOCH_RADIALINTEGRATIONMAPPINGFPGA_H diff --git a/image_analysis/RadialIntegrationProfile.cpp b/image_analysis/RadialIntegrationProfile.cpp index 282a931a..b7c28604 100644 --- a/image_analysis/RadialIntegrationProfile.cpp +++ b/image_analysis/RadialIntegrationProfile.cpp @@ -3,6 +3,13 @@ #include "RadialIntegrationProfile.h" #include "../common/JFJochException.h" +inline float sum_to_count(int64_t sum, uint64_t count) { + if (count == 0) + return 0.0f; + else + return static_cast(static_cast(sum) / (static_cast(count * (1LU<<24)))); +} + RadialIntegrationProfile::RadialIntegrationProfile(RadialIntegrationMapping &mapping, const DiffractionExperiment& experiment) : bin_to_q(mapping.GetBinToQ()), @@ -10,7 +17,7 @@ RadialIntegrationProfile::RadialIntegrationProfile(RadialIntegrationMapping &map count(mapping.GetBinNumber(), 0){ } -void RadialIntegrationProfile::Add(const std::vector &in_sum, const std::vector &in_count) { +void RadialIntegrationProfile::Add(const std::vector &in_sum, const std::vector &in_count) { std::unique_lock ul(m); if ((in_sum.size() == sum.size()) && (in_count.size() == count.size())) { @@ -24,10 +31,8 @@ void RadialIntegrationProfile::Add(const std::vector &in_sum, const std:: std::vector RadialIntegrationProfile::GetResult() const { std::vector rad_int_profile(sum.size(), 0); - for (int i = 0; i < sum.size(); i++) { - if (count[i] > 0) - rad_int_profile[i] = static_cast(sum[i]) / static_cast(count[i]); - } + for (int i = 0; i < sum.size(); i++) + rad_int_profile[i] = sum_to_count(sum[i], count[i]); return rad_int_profile; } @@ -43,8 +48,8 @@ void RadialIntegrationProfile::GetPlot(JFJochProtoBuf::Plot &plot) const { float RadialIntegrationProfile::GetMeanValueOfBins(uint16_t min_bin, uint16_t max_bin) const { std::unique_lock ul(m); - float ret_sum = 0; - float ret_count = 0; + int64_t ret_sum = 0; + uint64_t ret_count = 0; for (int i = std::min((uint16_t)sum.size(),min_bin); i <= std::min((uint16_t)(sum.size()-1),max_bin); @@ -53,10 +58,7 @@ float RadialIntegrationProfile::GetMeanValueOfBins(uint16_t min_bin, uint16_t ma ret_count += count[i]; } - if (ret_count == 0) - return 0; - else - return ret_sum / ret_count; + return sum_to_count(ret_sum, ret_count); } RadialIntegrationProfile &RadialIntegrationProfile::operator+=(const RadialIntegrationProfile &other) { @@ -84,7 +86,7 @@ void RadialIntegrationProfile::Add(const DeviceOutput &result) { "Error in getting result from FPGA"); for (int i = 0; i < sum.size(); i++) { - sum[i] += static_cast(static_cast(result.integration_result[i].sum) / (1LU << 24)); + sum[i] += result.integration_result[i].sum; count[i] += result.integration_result[i].count; } } \ No newline at end of file diff --git a/image_analysis/RadialIntegrationProfile.h b/image_analysis/RadialIntegrationProfile.h index fc9d4326..5959322e 100644 --- a/image_analysis/RadialIntegrationProfile.h +++ b/image_analysis/RadialIntegrationProfile.h @@ -12,13 +12,13 @@ class RadialIntegrationProfile { mutable std::mutex m; - std::vector sum; - std::vector count; + std::vector sum; + std::vector count; std::vector bin_to_q; public: explicit RadialIntegrationProfile(RadialIntegrationMapping &mapping, const DiffractionExperiment &experiment); void Add(const DeviceOutput &result); - void Add(const std::vector &sum, const std::vector &count); + void Add(const std::vector &sum, const std::vector &count); std::vector GetResult() const; float GetMeanValueOfBins(uint16_t min_bin, uint16_t max_bin) const; void GetPlot(JFJochProtoBuf::Plot &plot) const; diff --git a/receiver/AcquisitionDevice.cpp b/receiver/AcquisitionDevice.cpp index dcc4a284..1f677e2a 100644 --- a/receiver/AcquisitionDevice.cpp +++ b/receiver/AcquisitionDevice.cpp @@ -13,7 +13,6 @@ #include "../common/JFJochException.h" #include "AcquisitionDevice.h" #include "../common/NetworkAddressConvert.h" -#include "../image_analysis/RadialIntegrationMappingFPGA.h" void *mmap_acquisition_buffer(size_t size, int16_t numa_node) { void *ret = mmap(nullptr, size, PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_ANONYMOUS, -1, 0); diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index cd43de36..ce05c621 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -93,6 +93,7 @@ JFJochReceiver::JFJochReceiver(const JFJochProtoBuf::ReceiverInput &settings, rad_int_mapping = std::make_unique(experiment); rad_int_profile = std::make_unique(*rad_int_mapping, experiment); rad_int_corr = CalcRadIntCorr(experiment); + rad_int_corr_raw = CalcRadIntCorrRawCoord(experiment); for (int i = 0; i < experiment.GetDataFileCount(); i++) rad_int_profile_per_file.emplace_back( @@ -102,8 +103,6 @@ JFJochReceiver::JFJochReceiver(const JFJochProtoBuf::ReceiverInput &settings, } for (int d = 0; d < ndatastreams; d++) { - if (calib) - acquisition_device[d]->InitializeCalibration(experiment, calib.value()); acquisition_device[d]->PrepareAction(experiment); logger.Debug("Acquisition device {} prepared", d); @@ -253,6 +252,13 @@ void JFJochReceiver::AcquireThread(uint16_t data_stream) { } try { + if (calib) + acquisition_device[data_stream]->InitializeCalibration(experiment, calib.value()); + + if (rad_int_mapping) + acquisition_device[data_stream]->InitializeIntegrationMap(experiment, + rad_int_mapping->GetPixelToBinMappingRaw(), + rad_int_corr_raw); frame_transformation_ready.wait(); logger.Debug("Device thread {} start FPGA action", data_stream); acquisition_device[data_stream]->StartAction(experiment); @@ -337,7 +343,8 @@ void JFJochReceiver::MeasurePedestalThread(uint16_t data_stream, uint16_t module void JFJochReceiver::MiniSummationThread(int d, int m, size_t image_number, bool &send_image, FrameTransformation &transformation, DataMessage &message, - std::vector &conversion) { + std::vector &conversion, + RadialIntegrationProfile *profile) { for (int j = 0; j < experiment.GetSummation(); j++) { size_t frame_number = image_number * experiment.GetSummation() + j; acquisition_device[d]->Counters().WaitForFrame(frame_number + 2); @@ -360,6 +367,9 @@ void JFJochReceiver::MiniSummationThread(int d, int m, size_t image_number, bool } else src = acquisition_device[d]->GetErrorFrameBuffer(); + if (profile) + profile->Add(*acquisition_device[d]->GetDeviceOutput(frame_number, m)); + if (experiment.GetConversionOnCPU()) { auto &conv = conversion.at(experiment.GetFirstModuleOfDataStream(d) + m); transformation.ProcessModule(conv, src, m, d); @@ -393,15 +403,9 @@ void JFJochReceiver::FrameTransformationThread() { std::unique_ptr spot_finder; try { - if (rad_int_mapping) { - spot_finder = std::make_unique(experiment.GetXPixelsNum(), - experiment.GetYPixelsNum(), - one_byte_mask, *rad_int_mapping); - spot_finder->LoadRadialIntegrationCorr(rad_int_corr); - } else - spot_finder = std::make_unique(experiment.GetXPixelsNum(), - experiment.GetYPixelsNum(), - one_byte_mask); + spot_finder = std::make_unique(experiment.GetXPixelsNum(), + experiment.GetYPixelsNum(), + one_byte_mask); spot_finder->SetInputBuffer(transformation.GetPreview16BitImage()); spot_finder->RegisterBuffer(); } catch (const JFJochException& e) { @@ -440,11 +444,14 @@ void JFJochReceiver::FrameTransformationThread() { bool send_image = false; // We send image if at least one module was collected in full + std::unique_ptr rad_int_profile_image; + + if (rad_int_mapping) + rad_int_profile_image = std::make_unique(*rad_int_mapping, experiment); + if (GPUImageAnalysis::GPUPresent()) { if (((spotfinder_stride > 0) && (image_number % spotfinder_stride == 0)) || send_preview) { calculate_spots = true; - if (rad_int_mapping) - send_bkg_estimate = true; } } @@ -455,7 +462,8 @@ void JFJochReceiver::FrameTransformationThread() { mini_summation_threads.emplace_back(std::async(std::launch::async, &JFJochReceiver::MiniSummationThread, this, d, m, image_number, std::ref(send_image), std::ref(transformation), - std::ref(message), std::ref(*conversion))); + std::ref(message), std::ref(*conversion), + rad_int_profile_image.get())); message.receiver_aq_dev_delay = max_delay; } else { for (int j = 0; j < experiment.GetSummation(); j++) { @@ -489,6 +497,9 @@ void JFJochReceiver::FrameTransformationThread() { } else transformation.ProcessModule(src, m, d); + if (rad_int_profile_image) + rad_int_profile_image->Add(*acquisition_device[d]->GetDeviceOutput(frame_number, m)); + acquisition_device[d]->FrameBufferRelease(frame_number, m); } auto delay = acquisition_device[d]->Counters().CalculateDelay(frame_number); @@ -512,10 +523,7 @@ void JFJochReceiver::FrameTransformationThread() { if (calculate_spots) spot_finder->RunSpotFinder(local_data_processing_settings); - if (send_bkg_estimate) - spot_finder->RunRadialIntegration(); - - if (calculate_spots) { + if (calculate_spots) { spot_finder->GetSpotFinderResults(experiment, GetDataProcessingSettings(), spots); for (const auto &spot: spots) message.spots.push_back(spot); @@ -547,25 +555,21 @@ void JFJochReceiver::FrameTransformationThread() { } } - if (send_bkg_estimate) { - - RadialIntegrationProfile rad_int_profile_image(*rad_int_mapping, experiment); - rad_int_profile_image.Add(spot_finder->GetRadialIntegrationSum(), - spot_finder->GetRadialIntegrationCount()); + if (rad_int_profile_image) { uint16_t rad_int_min_bin = std::floor( rad_int_mapping->QToBin(local_data_processing_settings.bkg_estimate_low_q())); uint16_t rad_int_max_bin = std::ceil( rad_int_mapping->QToBin(local_data_processing_settings.bkg_estimate_high_q())); - float bkg_estimate_val = rad_int_profile_image.GetMeanValueOfBins(rad_int_min_bin, rad_int_max_bin); + float bkg_estimate_val = rad_int_profile_image->GetMeanValueOfBins(rad_int_min_bin, rad_int_max_bin); bkg_estimate.AddElement(image_number, bkg_estimate_val); - message.rad_int_profile = rad_int_profile_image.GetResult(); + message.rad_int_profile = rad_int_profile_image->GetResult(); - *rad_int_profile += rad_int_profile_image; + *rad_int_profile += *rad_int_profile_image; if (image_number % experiment.GetDataFileCount() < rad_int_profile_per_file.size()) - *rad_int_profile_per_file[image_number % experiment.GetDataFileCount()] += rad_int_profile_image; + *rad_int_profile_per_file[image_number % experiment.GetDataFileCount()] += *rad_int_profile_image; } if (send_preview) diff --git a/receiver/JFJochReceiver.h b/receiver/JFJochReceiver.h index 680e25cf..9fed7c17 100644 --- a/receiver/JFJochReceiver.h +++ b/receiver/JFJochReceiver.h @@ -46,6 +46,7 @@ class JFJochReceiver { std::unique_ptr rad_int_profile; std::vector> rad_int_profile_per_file; std::vector rad_int_corr; + std::vector rad_int_corr_raw; std::vector spot_finder_mask; @@ -109,7 +110,8 @@ class JFJochReceiver { void MeasurePedestalThread(uint16_t data_stream, uint16_t module_number, uint16_t storage_cell); void MiniSummationThread(int d, int m, size_t image_number, bool &send_image, FrameTransformation &transformation, DataMessage &message, - std::vector &conversion); + std::vector &conversion, + RadialIntegrationProfile *profile); void Cancel(const JFJochException &e); void FinalizeMeasurement(); JFJochProtoBuf::DataProcessingSettings GetDataProcessingSettings(); diff --git a/tests/JFJochFullIntegrationTest.cpp b/tests/JFJochFullIntegrationTest.cpp index 676b2265..669012a2 100644 --- a/tests/JFJochFullIntegrationTest.cpp +++ b/tests/JFJochFullIntegrationTest.cpp @@ -9,6 +9,7 @@ #include "../receiver/JFJochReceiverService.h" #include "FPGAUnitTest.h" #include "../receiver/MockAcquisitionDevice.h" +#include "../receiver/HLSSimulatedDevice.h" #include "../common/ZMQImagePusher.h" #include "../common/jsonToGrpc.h" @@ -1135,128 +1136,6 @@ TEST_CASE("JFJochIntegrationTest_ZMQ_with_preview_no_writer_binning2x2", "[JFJoc writer_server->Shutdown(); } -TEST_CASE("JFJochIntegrationTest_ZMQ_background_estimation", "[JFJochReceiver]") { - Logger logger("JFJochIntegrationTest_ZMQ_background_estimation"); - - RegisterHDF5Filter(); - - size_t nimages = 10; - int64_t ndatastream = 1; - int64_t nmodules = 8; - - JFJochServices services(logger); - JFJochStateMachine state_machine(services, logger); - - state_machine.AddDetectorSetup(DetectorGeometry(ndatastream * nmodules, 2, 8, 36)); - state_machine.NotThreadSafe_Experiment().DataStreams(ndatastream); - state_machine.NotThreadSafe_Experiment().PedestalG0Frames(0).PedestalG1Frames(0).PedestalG2Frames(0) - .MaskChipEdges(true).MaskModuleEdges(true).SpotFindingPeriod(1ms).FrameTime(1ms); - services.Writer("unix:writer_test", "inproc://#1").Receiver("unix:fpga_receiver_test"); - - logger.Verbose(true); - - std::vector image_raw_geom(RAW_MODULE_SIZE); - - std::vector> aq_devices; - - auto *test = new MockAcquisitionDevice(0, 256); - - for (int i = 0; i < nimages; i++) { - for (auto &j: image_raw_geom) - j = 2 * i + 5; - for (int m = 0; m < state_machine.NotThreadSafe_Experiment().GetModulesNum(); m++) - test->AddModule(i + 1, m, image_raw_geom.data()); - } - - test->Terminate(); - - aq_devices.emplace_back(test); - - ZMQContext zmq_context; - - std::vector tmp_devices; - for (const auto &i: aq_devices) - tmp_devices.emplace_back(i.get()); - - ZMQSocket rcv_spots_socket(zmq_context, ZMQSocketType::Sub); - REQUIRE_NOTHROW(rcv_spots_socket.Connect("inproc://#2")); - rcv_spots_socket.SubscribeAll(); - - ZMQImagePusher pusher(zmq_context, {"inproc://#1"}); - JFJochReceiverService fpga_receiver(tmp_devices, logger, pusher); - JFJochWriterService writer(zmq_context, logger); - - auto fpga_receiver_server = gRPCServer("unix:fpga_receiver_test", fpga_receiver); - auto writer_server = gRPCServer("unix:writer_test", writer); - - REQUIRE_NOTHROW(state_machine.Initialize()); - logger.Info("Initialized"); - - JFJochProtoBuf::DatasetSettings setup; - setup.set_ntrigger(1); - setup.set_detector_distance_mm(75); - setup.set_summation(1); - setup.set_file_prefix("bkg_estimate_test"); - setup.set_images_per_trigger(nimages); - setup.set_photon_energy_kev(12.4); - setup.set_beam_x_pxl(1090); - setup.set_beam_y_pxl(1136); - setup.mutable_unit_cell()->set_a(36.9); - setup.mutable_unit_cell()->set_b(78.95); - setup.mutable_unit_cell()->set_c(78.95); - setup.mutable_unit_cell()->set_alpha(90.0); - setup.mutable_unit_cell()->set_beta(90.0); - setup.mutable_unit_cell()->set_gamma(90.0); - setup.set_data_file_count(1); - - JFJochProtoBuf::DataProcessingSettings settings; - settings.set_signal_to_noise_threshold(4); - settings.set_photon_count_threshold(5); - settings.set_min_pix_per_spot(3); - settings.set_max_pix_per_spot(200); - settings.set_low_resolution_limit(80.0); - settings.set_high_resolution_limit(2.5); - settings.set_local_bkg_size(5); - - settings.set_bkg_estimate_low_q(2 * M_PI / 4.0); - settings.set_bkg_estimate_high_q(2 * M_PI / 3.0); - - REQUIRE_NOTHROW(state_machine.Start(setup)); - logger.Info("Started measurement"); - - REQUIRE_NOTHROW(state_machine.Stop()); - logger.Info("Stopped measurement"); - - JFJochProtoBuf::PlotRequest req; - req.set_type(JFJochProtoBuf::BKG_ESTIMATE); - req.set_binning(10); - - JFJochProtoBuf::Plot plot; - REQUIRE(fpga_receiver.GetDataProcessingPlots(nullptr, &req, &plot).ok()); - - REQUIRE(plot.x_size() == 1); - REQUIRE(plot.y_size() == 1); - REQUIRE(plot.x(0) == Approx(4.5)); - REQUIRE(plot.y(0) == Approx((5 + 23) / 2.0)); - - req.set_type(JFJochProtoBuf::BKG_ESTIMATE); - req.set_binning(1); - REQUIRE(fpga_receiver.GetDataProcessingPlots(nullptr, &req, &plot).ok()); - - REQUIRE(plot.x_size() == 10); - REQUIRE(plot.y_size() == 10); - - auto tmp = state_machine.GetMeasurementStatistics(); - REQUIRE(tmp.has_value()); - auto statistics = tmp.value(); - - REQUIRE(statistics.collection_efficiency() == 1.0); - REQUIRE(statistics.images_collected() == 10); - - fpga_receiver_server->Shutdown(); - writer_server->Shutdown(); -} - TEST_CASE("JFJochIntegrationTest_ZMQ_lysozyme_spot", "[JFJochReceiver]") { Logger logger("JFJochIntegrationTest_ZMQ_lysozyme_spot_and_index"); diff --git a/tests/RadialIntegrationTest.cpp b/tests/RadialIntegrationTest.cpp index dc9c3c4c..19ec84b6 100644 --- a/tests/RadialIntegrationTest.cpp +++ b/tests/RadialIntegrationTest.cpp @@ -3,8 +3,6 @@ #include #include "../image_analysis/RadialIntegrationProfile.h" -#include "../image_analysis/RadialIntegrationMappingFPGA.h" -#include "../image_analysis/RadialIntegration.h" TEST_CASE("RadialIntegrationMapping_Constructor","[RadialIntegration]") { DiffractionExperiment x; @@ -15,15 +13,6 @@ TEST_CASE("RadialIntegrationMapping_Constructor","[RadialIntegration]") { REQUIRE_NOTHROW(radial = std::make_unique(x)); } -TEST_CASE("RadialIntegrationMappingFPGA_Constructor","[RadialIntegration]") { - DiffractionExperiment x; - - std::unique_ptr radial; - - x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(5); - REQUIRE_NOTHROW(radial = std::make_unique(x, 0)); -} - TEST_CASE("RadialIntegrationMapping_GetBinNumber","[RadialIntegration]") { DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); @@ -32,14 +21,6 @@ TEST_CASE("RadialIntegrationMapping_GetBinNumber","[RadialIntegration]") { REQUIRE(mapping.GetBinNumber() == 39); } -TEST_CASE("RadialIntegrationMappingFPGA_GetBinNumber","[RadialIntegration]") { - DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); - x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); - x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(4); - RadialIntegrationMappingFPGA mapping(x, 0); - REQUIRE(mapping.GetBinNumber() == 39); -} - TEST_CASE("RadialIntegrationMapping_GetBinNumber_DetectorLimit","[RadialIntegration]") { DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); @@ -49,15 +30,6 @@ TEST_CASE("RadialIntegrationMapping_GetBinNumber_DetectorLimit","[RadialIntegrat REQUIRE(mapping.GetBinNumber() < 80); } -TEST_CASE("RadialIntegrationMappingFPGA_GetBinNumber_DetectorLimit","[RadialIntegration]") { - DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); - x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); - - x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(9.9); - RadialIntegrationMappingFPGA mapping(x, 0); - REQUIRE(mapping.GetBinNumber() < 80); -} - TEST_CASE("RadialIntegrationMapping_GetBinToQ","[RadialIntegration]") { DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); @@ -72,20 +44,6 @@ TEST_CASE("RadialIntegrationMapping_GetBinToQ","[RadialIntegration]") { REQUIRE(bin_to_q[38] == Approx(3.95)); } -TEST_CASE("RadialIntegrationMappingFPGA_GetBinToQ","[RadialIntegration]") { - DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); - x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); - x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(4); - - RadialIntegrationMappingFPGA mapping(x, 0); - auto bin_to_q = mapping.GetBinToQ(); - - REQUIRE(bin_to_q[0] == Approx(0.15)); - REQUIRE(bin_to_q[1] == Approx(0.25)); - REQUIRE(bin_to_q[15] == Approx(1.65)); - REQUIRE(bin_to_q[38] == Approx(3.95)); -} - TEST_CASE("RadialIntegrationMapping_QToBin","[RadialIntegration]") { DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); @@ -99,19 +57,6 @@ TEST_CASE("RadialIntegrationMapping_QToBin","[RadialIntegration]") { REQUIRE(mapping.QToBin(50.0) == Approx(38)); } -TEST_CASE("RadialIntegrationMappingFPGA_QToBin","[RadialIntegration]") { - DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); - x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); - x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(4); - - RadialIntegrationMappingFPGA mapping(x, 0); - - REQUIRE(mapping.QToBin(0.0) == 0); - REQUIRE(std::floor(mapping.QToBin(0.200001)) == 1); - REQUIRE(mapping.QToBin(0.6) == Approx(5)); - REQUIRE(mapping.QToBin(50.0) == Approx(38)); -} - TEST_CASE("RadialIntegrationProfile","[RadialIntegration]") { DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); @@ -121,17 +66,17 @@ TEST_CASE("RadialIntegrationProfile","[RadialIntegration]") { RadialIntegrationProfile profile(mapping, x); - std::vector sum(mapping.GetBinNumber()); - std::vector count(mapping.GetBinNumber()); + std::vector sum(mapping.GetBinNumber()); + std::vector count(mapping.GetBinNumber()); for (int i = 0; i < mapping.GetBinNumber(); i++) { - sum[i] = i * i * 4; + sum[i] = (i * i * 4) * (1LU<<24); count[i] = i; } REQUIRE_NOTHROW(profile.Add(sum, count)); REQUIRE_NOTHROW(profile.Add(sum, count)); - std::vector sum_wr(mapping.GetBinNumber() - 1); + std::vector sum_wr(mapping.GetBinNumber() - 1); REQUIRE_THROWS(profile.Add(sum_wr, count)); JFJochProtoBuf::Plot plot; @@ -154,11 +99,11 @@ TEST_CASE("RadialIntegrationProfile_operatorAdd","[RadialIntegration]") { RadialIntegrationProfile profile0(mapping, x), profile1(mapping, x); - std::vector sum(mapping.GetBinNumber()); - std::vector count(mapping.GetBinNumber()); + std::vector sum(mapping.GetBinNumber()); + std::vector count(mapping.GetBinNumber()); for (int i = 0; i < mapping.GetBinNumber(); i++) { - sum[i] = i * i * 4; + sum[i] = (i * i * 4) * (1LU<<24); count[i] = i; } REQUIRE_NOTHROW(profile0.Add(sum, count)); @@ -184,166 +129,15 @@ TEST_CASE("RadialIntegrationProfile_GetMeanValueOfBins","[RadialIntegration]") { RadialIntegrationProfile profile(mapping, x); - std::vector sum(mapping.GetBinNumber()); - std::vector count(mapping.GetBinNumber()); + std::vector sum(mapping.GetBinNumber()); + std::vector count(mapping.GetBinNumber()); for (int i = 0; i < mapping.GetBinNumber(); i++) { - sum[i] = i * i * 4; + sum[i] = (i * i * 4) * (1LU<<24); count[i] = i; } REQUIRE_NOTHROW(profile.Add(sum, count)); - REQUIRE(profile.GetMeanValueOfBins(0,2) == Approx ((sum[0] + sum[1] + sum[2]) / (count[0] + count[1] + count[2]))); - REQUIRE(profile.GetMeanValueOfBins(5,7) == Approx ((sum[5] + sum[6] + sum[7]) / (count[5] + count[6] + count[7]))); + REQUIRE(profile.GetMeanValueOfBins(0,2) == Approx ((sum[0] + sum[1] + sum[2]) /(1LU<<24) / double(count[0] + count[1] + count[2]))); + REQUIRE(profile.GetMeanValueOfBins(5,7) == Approx ((sum[5] + sum[6] + sum[7]) /(1LU<<24) / double (count[5] + count[6] + count[7]))); } - -TEST_CASE("RadialIntegration_Process","[RadialIntegration]") { - std::vector pixel_to_bin = {0,1,2,4,3,1,2,3}; - std::vector test_image = {7,6,5,4,3,2,1,0}; - std::vector result; - - RadialIntegration radial(pixel_to_bin, 4); - radial.ProcessOneImage(test_image.data(), 8); - - REQUIRE(radial.GetCount().size() == 4); - REQUIRE(radial.GetSum().size() == 4); - - size_t count = 0; - for (const auto &i: radial.GetCount()) - count += i; - REQUIRE(count == 7); // bin number 4 is excluded - - REQUIRE(radial.GetCount()[0] == 1); - REQUIRE(radial.GetCount()[1] == 2); - - REQUIRE(radial.GetSum()[0] == 7); - REQUIRE(radial.GetSum()[1] == 6+2); -} - -TEST_CASE("RadialIntegration_GetResult","[RadialIntegration]") { - std::vector pixel_to_bin = {0,1,2,4,3,1,2,3}; - std::vector test_image = {7,6,5,4,3,2,1,0}; - std::vector result; - - RadialIntegration radial(pixel_to_bin, 4); - radial.ProcessOneImage(test_image.data(), 8); - radial.GetResult(result); - - REQUIRE(result.size() == 4); - - REQUIRE(result[0] == 7); - REQUIRE(result[1] == (6 + 2) / 2.0); - REQUIRE(result[2] == (5 + 1) / 2.0); - REQUIRE(result[3] == (3 + 0) / 2.0); -} - -TEST_CASE("RadialIntegration_GetRangeValue","[RadialIntegration]") { - std::vector pixel_to_bin = {0,1,2,4,3,1,2,3}; - std::vector test_image = {7,6,5,4,3,2,1,0}; - std::vector result; - - RadialIntegration radial(pixel_to_bin, 4); - radial.ProcessOneImage(test_image.data(), 8); - REQUIRE(radial.GetRangeValue(0, 7) == Approx((0+1+2+3+5+6+7) / 7.0)); - REQUIRE(radial.GetRangeValue(2, 2) == Approx((5+1) / 2.0)); - REQUIRE(radial.GetRangeValue(2, 3) == Approx((5+3+1+0) / 4.0)); - REQUIRE(radial.GetRangeValue(15, 15) == 0); // Empty set -} - -#include "../image_analysis/GPUImageAnalysis.h" - -TEST_CASE("RadialIntegrationGPU_Process","[RadialIntegration]") { - std::vector pixel_to_bin = {0,1,2,4,3,1,2,3}; - std::vector test_image = {7,6,5,4,3,2,1,0}; - std::vector one_byte_mask = {1,1,1,1,1,1,1,1}; - - GPUImageAnalysis image_analysis(8, 1, one_byte_mask, pixel_to_bin, 4); - image_analysis.SetInputBuffer(test_image.data()); - image_analysis.LoadDataToGPU(); - image_analysis.RunRadialIntegration(); - - REQUIRE(image_analysis.GetRadialIntegrationCount().size() == 4); - REQUIRE(image_analysis.GetRadialIntegrationSum().size() == 4); - - size_t count = 0; - for (const auto &i: image_analysis.GetRadialIntegrationCount()) - count += i; - REQUIRE(count == 7); // bin number 4 is excluded - - REQUIRE(image_analysis.GetRadialIntegrationCount()[0] == 1); - REQUIRE(image_analysis.GetRadialIntegrationCount()[1] == 2); - - REQUIRE(image_analysis.GetRadialIntegrationSum()[0] == 7); - REQUIRE(image_analysis.GetRadialIntegrationSum()[1] == 6+2); -} - -TEST_CASE("RadialIntegrationGPU_Process_Corr","[RadialIntegration]") { - std::vector pixel_to_bin = {0,1,2,4,3,1,2,3}; - std::vector test_image = {7,6,5,4,3,2,1,0}; - std::vector one_byte_mask = {1,1,1,1,1,1,1,1}; - - std::vector corr = {0.5f, 0.3f, 0.5f, 0.2f, 0.8f, 0.9f, 1.0f, 0.5f}; - GPUImageAnalysis image_analysis(8, 1, one_byte_mask, pixel_to_bin, 4); - image_analysis.SetInputBuffer(test_image.data()); - image_analysis.LoadRadialIntegrationCorr(corr); - image_analysis.LoadDataToGPU(); - image_analysis.RunRadialIntegration(); - - REQUIRE(image_analysis.GetRadialIntegrationCount().size() == 4); - REQUIRE(image_analysis.GetRadialIntegrationSum().size() == 4); - - REQUIRE(image_analysis.GetRadialIntegrationCount()[0] == 1); - REQUIRE(image_analysis.GetRadialIntegrationCount()[1] == 2); - - REQUIRE(image_analysis.GetRadialIntegrationSum()[0] == Approx(7 * 0.5f)); - REQUIRE(image_analysis.GetRadialIntegrationSum()[1] == Approx(6*0.3f + 2*0.9f)); - REQUIRE(image_analysis.GetRadialIntegrationSum()[2] == Approx(5*0.5f + 1*1.0f)); - REQUIRE(image_analysis.GetRadialIntegrationSum()[3] == Approx(3*0.8f)); -} - -TEST_CASE("RadialIntegrationGPU_Process_Mask","[RadialIntegration]") { - std::vector pixel_to_bin = {0,1,2,4,3,1,2,3}; - std::vector test_image = {7,6,5,4,3,2,1,0}; - std::vector one_byte_mask = {0,0,1,1,1,1,1,1}; - - GPUImageAnalysis image_analysis(8, 1, one_byte_mask, pixel_to_bin, 4); - image_analysis.SetInputBuffer(test_image.data()); - image_analysis.LoadDataToGPU(); - image_analysis.RunRadialIntegration(); - - REQUIRE(image_analysis.GetRadialIntegrationCount().size() == 4); - REQUIRE(image_analysis.GetRadialIntegrationSum().size() == 4); - - size_t count = 0; - for (const auto &i: image_analysis.GetRadialIntegrationCount()) - count += i; - REQUIRE(count == 8 - 3); // bin number 0,1,4 is excluded - - REQUIRE(image_analysis.GetRadialIntegrationCount()[0] == 0); - REQUIRE(image_analysis.GetRadialIntegrationCount()[1] == 1); - - REQUIRE(image_analysis.GetRadialIntegrationSum()[0] == 0); - REQUIRE(image_analysis.GetRadialIntegrationSum()[1] == 2); -} - -TEST_CASE("RadialIntegrationGPU_GetResult","[RadialIntegration]") { - std::vector pixel_to_bin = {0,1,2,4,3,1,2,3}; - std::vector test_image = {7,6,5,4,3,2,1,0}; - std::vector one_byte_mask = {1,1,1,1,1,1,1,1}; - std::vector result; - - GPUImageAnalysis image_analysis(8, 1, one_byte_mask, pixel_to_bin, 4); - image_analysis.SetInputBuffer(test_image.data()); - image_analysis.LoadDataToGPU(); - image_analysis.RunRadialIntegration(); - image_analysis.GetRadialIntegrationProfile(result); - - REQUIRE(result.size() == 4); - - REQUIRE(result[0] == 7); - REQUIRE(result[1] == (6 + 2) / 2.0); - REQUIRE(result[2] == (5 + 1) / 2.0); - REQUIRE(result[3] == (3 + 0) / 2.0); -} - - diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 6c3ad3c0..3c6c95db 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -1,5 +1,3 @@ -ADD_EXECUTABLE(RadialIntDataset RadialIntDataset.cpp) -TARGET_LINK_LIBRARIES(RadialIntDataset ImageAnalysis CommonFunctions HDF5Wrappers) ADD_EXECUTABLE(jfjoch_udp_simulator jfjoch_udp_simulator.cpp UDPSimulator.cpp UDPSimulator.h) TARGET_LINK_LIBRARIES(jfjoch_udp_simulator CommonFunctions) @@ -13,14 +11,11 @@ target_link_libraries(HDF5DatasetWriteTest JFJochWriter CommonFunctions) ADD_EXECUTABLE(DataAnalysisPerfTest DataAnalysisPerfTest.cpp) TARGET_LINK_LIBRARIES(DataAnalysisPerfTest ImageAnalysis JFJochWriter CommonFunctions) -ADD_EXECUTABLE(RadialIntegrationCPUTest RadialIntegrationCPUTest.cpp) -TARGET_LINK_LIBRARIES(RadialIntegrationCPUTest ImageAnalysis JFJochWriter CommonFunctions) - ADD_EXECUTABLE(PreviewTest PreviewTest.cpp) TARGET_LINK_LIBRARIES(PreviewTest JFJochWriter CommonFunctions) ADD_EXECUTABLE(JFCalibrationPerfTest JFCalibrationPerfTest.cpp) TARGET_LINK_LIBRARIES(JFCalibrationPerfTest CommonFunctions) -INSTALL(TARGETS jfjoch_udp_simulator RadialIntDataset CompressionBenchmark HDF5DatasetWriteTest DataAnalysisPerfTest +INSTALL(TARGETS jfjoch_udp_simulator CompressionBenchmark HDF5DatasetWriteTest DataAnalysisPerfTest PreviewTest JFCalibrationPerfTest RUNTIME) \ No newline at end of file diff --git a/tools/DataAnalysisPerfTest.cpp b/tools/DataAnalysisPerfTest.cpp index 6b6d50c3..6c1ea348 100644 --- a/tools/DataAnalysisPerfTest.cpp +++ b/tools/DataAnalysisPerfTest.cpp @@ -10,7 +10,6 @@ #include "../image_analysis/IndexerWrapper.h" -#include "../image_analysis/RadialIntegration.h" #include "../common/Logger.h" #include "../image_analysis/PredictSpotsOnDetector.h" @@ -32,12 +31,10 @@ auto TestAll(const DiffractionExperiment &experiment, const JFJochProtoBuf::Data spot_finder.LoadDataToGPU(); spot_finder.RunSpotFinder(settings); spot_finder.GetSpotFinderResults(experiment, settings, spots); - spot_finder.RunRadialIntegration(); std::vector recip; for (const auto& s: spots) recip.emplace_back(s.ReciprocalCoord(experiment)); indexer.Run(recip); - spot_finder.GetRadialIntegrationProfile(result); } auto end_time = std::chrono::system_clock::now(); @@ -68,12 +65,10 @@ auto TestAllWithROI(const DiffractionExperiment &experiment, const JFJochProtoBu spot_finder.LoadDataToGPU(); spot_finder.RunSpotFinder(settings); spot_finder.GetSpotFinderResults(experiment, settings, spots); - spot_finder.RunRadialIntegration(); std::vector recip; for (const auto& s: spots) recip.emplace_back(s.ReciprocalCoord(experiment)); auto indexer_ret = indexer.Run(recip); - spot_finder.GetRadialIntegrationProfile(result); if (!indexer_ret.empty()) { PredictSpotsOnDetector(filter, experiment, indexer_ret[0].l); filter.Apply(roi_image, INT16_MIN); @@ -178,89 +173,6 @@ auto TestSpotFinderWithoutCopyToDevice(const DiffractionExperiment &experiment, return strstream.str(); } -auto TestRadialIntegrationGPU(const DiffractionExperiment &x, - const JFJochProtoBuf::DataProcessingSettings &settings, - int16_t* image, size_t nimages) { - - uint32_t nredo = 20; - RadialIntegrationMapping mapping(x); - - std::vector one_byte_mask(x.GetPixelsNum(), 1); - GPUImageAnalysis integration(x.GetXPixelsNum(), x.GetYPixelsNum(), one_byte_mask, mapping); - - integration.SetInputBuffer(image); - - std::vector result; - - auto start_time = std::chrono::system_clock::now(); - for (int redo = 0; redo < nredo; redo++) { - for (int i = 0; i < nimages; i++) { - integration.LoadDataToGPU(); - integration.RunRadialIntegration(); - integration.GetRadialIntegrationProfile(result); - } - } - auto end_time = std::chrono::system_clock::now(); - auto elapsed = std::chrono::duration_cast(end_time - start_time); - - return elapsed.count() / (1000.0 * (double) nimages * nredo); -} - -auto TestRadialIntegrationGPUWithoutCopyToDevice(const DiffractionExperiment &x, - const JFJochProtoBuf::DataProcessingSettings &settings, - int16_t* image, size_t nimages) { - - uint32_t nredo = 20; - RadialIntegrationMapping mapping(x); - - std::vector one_byte_mask(x.GetPixelsNum(), 1); - GPUImageAnalysis integration(x.GetXPixelsNum(), x.GetYPixelsNum(), one_byte_mask, mapping); - - integration.SetInputBuffer(image); - - std::vector result; - - integration.LoadDataToGPU(); - - auto start_time = std::chrono::system_clock::now(); - for (int redo = 0; redo < nredo; redo++) { - for (int i = 0; i < nimages; i++) { - integration.RunRadialIntegration(); - integration.GetRadialIntegrationProfile(result); - } - } - auto end_time = std::chrono::system_clock::now(); - auto elapsed = std::chrono::duration_cast(end_time - start_time); - - return elapsed.count() / (1000.0 * (double) nimages * nredo); -} - - -auto TestRadialIntegration(const DiffractionExperiment &experiment, - const JFJochProtoBuf::DataProcessingSettings &settings, - int16_t* image, size_t nimages) { - - uint32_t nredo = 20; - RadialIntegrationMapping mapping(experiment); - std::unique_ptr integration; - - integration = std::make_unique(mapping); - - std::vector result; - - auto start_time = std::chrono::system_clock::now(); - for (int redo = 0; redo < nredo; redo++) { - for (int i = 0; i < nimages; i++) { - integration->ProcessOneImage(image + i * experiment.GetPixelsNum(), experiment.GetPixelsNum()); - //integration.GetResult(result); - } - } - auto end_time = std::chrono::system_clock::now(); - auto elapsed = std::chrono::duration_cast(end_time - start_time); - - return elapsed.count() / (1000.0 * (double) nimages * nredo); -} - int main(int argc, char **argv) { RegisterHDF5Filter(); @@ -315,7 +227,6 @@ int main(int argc, char **argv) { settings.set_local_bkg_size(5); std::vector one_byte_mask(x.GetPixelsNum(), 1); - RadialIntegrationMapping mapping(x); logger.Info("COLSPOT NBX=NBY=5 (GPU)"); if (GPUImageAnalysis::GPUPresent()) { @@ -335,19 +246,11 @@ int main(int argc, char **argv) { TestSpotFinder(x, settings, local_peakfinder_gpu,image_conv.data(), nimages); } - logger.Info("{:30s} {:8.1f} ms/image", "Radial int. (GPU)", TestRadialIntegrationGPU(x, settings, - image_conv.data(), nimages)); - logger.Info("{:30s} {:8.1f} ms/image", "Radial int. (GPU/nocpy)", TestRadialIntegrationGPUWithoutCopyToDevice(x, settings, - image_conv.data(), nimages)); - - logger.Info("{:30s} {:8.1f} ms/image", "Radial int. (CPU)", TestRadialIntegration(x, settings, - image_conv.data(), nimages)); - TestIndexing(); logger.Info("Full package"); if (GPUImageAnalysis::GPUPresent()) { - GPUImageAnalysis local_peakfinder_gpu(x.GetXPixelsNum(), x.GetYPixelsNum(), one_byte_mask, mapping); + GPUImageAnalysis local_peakfinder_gpu(x.GetXPixelsNum(), x.GetYPixelsNum(), one_byte_mask); TestAll(x, settings, local_peakfinder_gpu,image_conv.data(), nimages); TestAllWithROI(x, settings, local_peakfinder_gpu,image_conv.data(), nimages); } diff --git a/tools/RadialIntDataset.cpp b/tools/RadialIntDataset.cpp deleted file mode 100644 index d8b13f18..00000000 --- a/tools/RadialIntDataset.cpp +++ /dev/null @@ -1,97 +0,0 @@ -// Copyright (2019-2023) Paul Scherrer Institute - -#include - -#include -#include -#include "../writer/HDF5Objects.h" -#include "../image_analysis/RadialIntegration.h" - -std::vector GetOneByteMask(DiffractionExperiment &x, HDF5Object &master_file) { - std::vector ret(x.GetPixelsNum(), 1); - - HDF5DataSet pixelmask_dataset(master_file, "/entry/instrument/detector/pixel_mask"); - - std::vector pm_start = {0, 0}; - std::vector pm_dim = {(hsize_t) x.GetYPixelsNum(), (hsize_t) x.GetXPixelsNum()}; - - std::vector pixel_mask(x.GetPixelsNum()); - - pixelmask_dataset.ReadVector(pixel_mask, pm_start, pm_dim); - - for (int i = 0; i < x.GetPixelsNum(); i++) { - if (pixel_mask[i] != 0) - ret[i] = 0; - } - return ret; -} - -void CalcRadialIntegration(DiffractionExperiment &x, RadialIntegration &rad, const char *file_name, - hsize_t first, hsize_t last) { - HDF5ReadOnlyFile data_file(file_name); - HDF5DataSet dataset(data_file, "/entry/data/data"); - HDF5DataSpace file_space(dataset); - - if (file_space.GetNumOfDimensions() != 3) { - std::cout << "/entry/data/data must be 3D" << std::endl; - exit(EXIT_FAILURE); - } - - std::vector source(x.GetPixelsNum()); - - first = std::min(first, file_space.GetDimensions()[0] - 1); - last = std::min(last, file_space.GetDimensions()[0] - 1); - std::cout << "# " << first << " " << last << std::endl; - for (hsize_t i = first; i < last; i++) { - std::vector start = {i, 0, 0}; - std::vector dim = {1, (hsize_t) x.GetYPixelsNum(), (hsize_t) x.GetXPixelsNum()}; - dataset.ReadVector(source, start, dim); - rad.Process(source.data(), x.GetPixelsNum()); - } -} - -void GetGeometry(DiffractionExperiment &x, HDF5Object &master_file) { - x.BeamX_pxl( - HDF5DataSet(master_file, "/entry/instrument/detector/beam_center_x").ReadScalar()); - x.BeamY_pxl( - HDF5DataSet(master_file, "/entry/instrument/detector/beam_center_y").ReadScalar()); - x.DetectorDistance_mm( - HDF5DataSet(master_file, "/entry/instrument/detector/detector_distance").ReadScalar() *1e3); -} - -int main(int argc, char **argv) { - - hsize_t first, last; - DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); - - RegisterHDF5Filter(); - - if (argc == 3) { - first = 0; - last = INT32_MAX; - } else if (argc == 5) { - first = std::strtol(argv[3], nullptr, 10); - last = std::strtol(argv[4], nullptr, 10); - } else { - std::cout << "Usage ./ImageAverage { }" << std::endl; - exit(EXIT_FAILURE); - } - - HDF5ReadOnlyFile master_file(argv[1]); - auto one_byte_mask = GetOneByteMask(x, master_file); - - GetGeometry(x, master_file); - - //TODO: Setup Mask - RadialIntegrationMapping rad_int_map(x); - RadialIntegration rad_int(rad_int_map); - - CalcRadialIntegration(x, rad_int, argv[2], first, last); - - auto result_x = rad_int_map.GetBinToQ(); - std::vector result_y; - rad_int.GetResult(result_y); - - for (int i = 0; i < result_y.size(); i++) - std::cout << result_x[i] << " " << result_y[i] << std::endl; -} diff --git a/tools/RadialIntegrationCPUTest.cpp b/tools/RadialIntegrationCPUTest.cpp deleted file mode 100644 index 15408fca..00000000 --- a/tools/RadialIntegrationCPUTest.cpp +++ /dev/null @@ -1,88 +0,0 @@ -// Copyright (2019-2023) Paul Scherrer Institute - -#include - -#include "../image_analysis/RadialIntegration.h" -#include "../common/Logger.h" -#include "../writer/HDF5Objects.h" - -Logger logger{"RadialIntegrationCPUTest"}; - -void RunRadialIntegrationThread(RadialIntegration *integration, - int16_t* image, size_t nimages, - size_t image0, size_t stride, - size_t npixel) { - for (size_t i = image0; i < nimages; i += stride) - integration->ProcessOneImage(image + i * npixel, npixel); -} - -auto TestRadialIntegration(const DiffractionExperiment &experiment, - int16_t* image, size_t nimages, - uint32_t nthreads, - uint32_t pixel_split) { - uint32_t nredo = 5; - RadialIntegrationMapping mapping(experiment); - - std::vector result; - std::vector> integration; - for (int i = 0; i < nthreads; i++) { - if (pixel_split == 1) - integration.emplace_back(std::make_unique(mapping)); - else - integration.emplace_back(std::make_unique(mapping.GetPixelToBinMappingSplitTo4(), - mapping.GetBinNumber(), 4)); - } - - auto start_time = std::chrono::system_clock::now(); - for (int redo = 0; redo < nredo; redo++) { - std::vector> futures; - for (int i = 0; i < nthreads; i++) { - futures.emplace_back(std::async(std::launch::async, &RunRadialIntegrationThread, - integration[i].get(), image, nimages, i, nthreads, - experiment.GetPixelsNum())); - } - for (auto &f: futures) - f.get(); - } - auto end_time = std::chrono::system_clock::now(); - auto elapsed = std::chrono::duration_cast(end_time - start_time); - - return elapsed.count() / (1000.0 * (double) nimages * nredo); -} - -int main(int argc, char **argv) { - if (argc > 3) { - logger.Info("Usage: ./DataAnalysisPerfTest {} {}"); - exit(EXIT_FAILURE); - } - - uint32_t nthreads = 1; - if (argc >= 2) - nthreads = atol(argv[1]); - - uint64_t nimages = 100; - if (argc == 3) - nimages = atol(argv[2]); - - if ((nthreads <= 0) || (nimages <= 0)) { - std::cerr << "Error in input parameters" << std::endl; - exit(EXIT_FAILURE); - } - - DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); - - logger.Info("Number of threads: {}", nthreads); - logger.Info("Number of images in the dataset: {}", nimages); - - x.Mode(DetectorMode::Conversion); - - std::vector image_conv ( nimages * x.GetPixelsNum(), 30); - - x.BeamX_pxl(1090).BeamY_pxl(1136).DetectorDistance_mm(75).PhotonEnergy_keV(WVL_1A_IN_KEV); - std::vector one_byte_mask(x.GetPixelsNum(), 1); - - logger.Info("{:30s} {:8.2f} ms/image", "Radial int. pxlspl 1 (CPU)", - TestRadialIntegration(x, image_conv.data(), nimages, nthreads, 1)); - logger.Info("{:30s} {:8.2f} ms/image", "Radial int. pxlspl 2 (CPU)", - TestRadialIntegration(x, image_conv.data(), nimages, nthreads, 4)); -}