mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 67 Wiki Activity

ImagePreprocessor: Dedicated class to handle preprocessing of diffraction images in Receiver Lite workflow + remove ROI calculation in Lite workflow

Browse Source
This commit is contained in:
leonarski_f
2026-04-22 10:37:41 +02:00
parent 1f65050618
commit cf420d5622
8 changed files with 223 additions and 153 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/CMakeLists.txt
+2 -1
View File
@@ -41,5 +41,6 @@ ADD_SUBDIRECTORY(indexing)
ADD_SUBDIRECTORY(geom_refinement)
ADD_SUBDIRECTORY(lattice_search)
ADD_SUBDIRECTORY(scale_merge)
ADD_SUBDIRECTORY(image_preprocessing)
TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge gemmi)
TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge gemmi)
image_analysis/MXAnalysisWithoutFPGA.cpp
+35 -139
View File
@@ -9,6 +9,7 @@
#include "spot_finding/SpotUtils.h"
#include "spot_finding/ImageSpotFinderFactory.h"
#include "bragg_prediction/BraggPredictionFactory.h"
#include "image_preprocessing/ImagePreprocessorCPU.h"
MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_experiment,
const AzimuthalIntegration &in_integration,
@@ -16,24 +17,16 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp
IndexAndRefine &in_indexer)
: experiment(in_experiment),
integration(in_integration),
roi_map(experiment.ExportROIMap()),
roi_names(experiment.ROI().GetROINameMap()),
roi_count(experiment.ROI().size()),
npixels(experiment.GetPixelsNum()),
xpixels(experiment.GetXPixelsNum()),
mask_1bit(npixels, false),
spotFinder(CreateImageSpotFinder(experiment.GetXPixelsNum(), experiment.GetYPixelsNum())),
indexer(in_indexer),
prediction(CreateBraggPrediction(experiment.IsRotationIndexing())),
updated_image(spotFinder->GetInputBuffer()),
azint_bins(in_integration.GetBinNumber()),
saturation_limit(experiment.GetSaturationLimit()),
mask(in_mask),
mask_resolution(experiment.GetPixelsNum(), false),
mask_high_res(-1),
mask_low_res(-1) {
for (int i = 0; i < npixels; i++)
mask_1bit[i] = (in_mask.GetMask().at(i) != 0);
preprocessor = std::make_unique<ImagePreprocessorCPU>(in_experiment, in_integration, in_mask);
}
void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
@@ -50,28 +43,40 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
if (output.image.GetCompressionAlgorithm() != CompressionAlgorithm::NO_COMPRESSION)
output.compression_time_s = std::chrono::duration<float>(compression_end_time - compression_start_time).count();
switch (output.image.GetMode()) {
case CompressedImageMode::Int8:
Analyze<int8_t>(output, image_ptr, INT8_MIN, INT8_MAX, profile, spot_finding_settings);
break;
case CompressedImageMode::Int16:
Analyze<int16_t>(output, image_ptr, INT16_MIN, INT16_MAX, profile, spot_finding_settings);
break;
case CompressedImageMode::Int32:
Analyze<int32_t>(output, image_ptr, INT32_MIN, INT32_MAX, profile, spot_finding_settings);
break;
case CompressedImageMode::Uint8:
Analyze<uint8_t>(output, image_ptr, UINT8_MAX, UINT8_MAX, profile, spot_finding_settings);
break;
case CompressedImageMode::Uint16:
Analyze<uint16_t>(output, image_ptr, UINT16_MAX, UINT16_MAX, profile, spot_finding_settings);
break;
case CompressedImageMode::Uint32:
Analyze<uint32_t>(output, image_ptr, UINT32_MAX, UINT32_MAX, profile, spot_finding_settings);
break;
default:
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "RGB/float mode not supported");
const auto preprocessing_start_time = std::chrono::steady_clock::now();
auto ret = preprocessor->Analyze(image_ptr, output.image.GetMode());
const auto preprocessing_end_time = std::chrono::steady_clock::now();
output.preprocessing_time_s = std::chrono::duration<float>(preprocessing_end_time - preprocessing_start_time).count();
if (spot_finding_settings.enable) {
// Update resolution mask
if (mask_high_res != spot_finding_settings.high_resolution_limit
|| mask_low_res != spot_finding_settings.low_resolution_limit)
UpdateMaskResolution(spot_finding_settings);
const auto spot_finding_start_time = std::chrono::steady_clock::now();
memcpy(spotFinder->GetInputBuffer().data(), preprocessor->GetProcessedImage().data(), npixels * sizeof(int32_t));
const std::vector<DiffractionSpot> spots = spotFinder->Run(spot_finding_settings, mask_resolution);
SpotAnalyze(experiment, spot_finding_settings, spots, output);
const auto spot_finding_end_time = std::chrono::steady_clock::now();
output.spot_finding_time_s = std::chrono::duration<float>(spot_finding_end_time - spot_finding_start_time).count();
if (spot_finding_settings.indexing)
indexer.ProcessImage(output, spot_finding_settings,
CompressedImage(preprocessor->GetProcessedImage(), experiment.GetXPixelsNum(), experiment.GetYPixelsNum()),
*prediction);
}
preprocessor->Update(profile);
output.max_viable_pixel_value = ret.max_value;
output.min_viable_pixel_value = ret.min_value;
output.error_pixel_count = ret.error_pixel_count;
output.saturated_pixel_count = ret.saturated_pixel_count;
output.az_int_profile = profile.GetResult();
output.bkg_estimate = profile.GetBkgEstimate(integration.Settings());
}
void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &settings) {
@@ -81,112 +86,3 @@ void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &sett
for (int i = 0; i < mask_resolution.size(); i++)
mask_resolution[i] = (resolution_map[i] > mask_low_res) || (resolution_map[i] < mask_high_res);
}
template<class T>
void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
const uint8_t *in_image,
T err_pixel_val,
T sat_pixel_val,
AzimuthalIntegrationProfile &profile,
const SpotFindingSettings &settings) {
const auto preprocessing_start_time = std::chrono::steady_clock::now();
auto image = reinterpret_cast<const T *>(in_image);
std::vector<ROIMessage> roi(roi_count);
std::vector<float> azim_sum(azint_bins, 0.0f);
//std::vector<float> azim_sum2(integration.GetBinNumber(), 0.0f);
std::vector<uint32_t> azim_count(azint_bins, 0);
size_t err_pixels = 0;
size_t masked_pixels = 0;
size_t sat_pixels = 0;
int64_t min_value = INT64_MAX;
int64_t max_value = INT64_MIN;
if (sat_pixel_val > saturation_limit)
sat_pixel_val = static_cast<T>(saturation_limit);
auto &pixel_to_bin = integration.GetPixelToBin();
auto &corrections = integration.Corrections();
profile.Clear(integration);
for (int i = 0; i < npixels; i++) {
if (mask_1bit[i] != 0) {
updated_image[i] = INT32_MIN;
++masked_pixels;
} else if (image[i] >= sat_pixel_val) {
updated_image[i] = INT32_MIN;
++sat_pixels;
} else if (std::is_signed<T>::value && (image[i] == err_pixel_val)) {
// Error pixels are possible only for signed types
updated_image[i] = INT32_MIN;
++err_pixels;
} else {
updated_image[i] = static_cast<int32_t>(image[i]);
if (image[i] > max_value)
max_value = image[i];
if (image[i] < min_value)
min_value = image[i];
if (roi_count > 0 && (roi_map[i] != 0)) {
int64_t x = i % xpixels;
int64_t y = i / xpixels;
for (int8_t r = 0; r < roi_count; r++) {
if ((roi_map[i] & (1 << r)) != 0) {
roi[r].sum += image[i];
roi[r].sum_square += image[i] * image[i];
roi[r].pixels += 1;
if (image[i] > roi[r].max_count)
roi[r].max_count = image[i];
roi[r].x_weighted += x * image[i];
roi[r].y_weighted += y * image[i];
}
}
}
const uint16_t bin = pixel_to_bin[i];
if (bin < azint_bins) {
float val = image[i] * corrections[i];
azim_sum[bin] += val;
//azim_sum2[bin] += val * val;
++azim_count[bin];
}
}
}
const auto preprocessing_end_time = std::chrono::steady_clock::now();
output.preprocessing_time_s = std::chrono::duration<float>(preprocessing_end_time - preprocessing_start_time).count();
if (settings.enable) {
// Update resolution mask
if (mask_high_res != settings.high_resolution_limit
|| mask_low_res != settings.low_resolution_limit)
UpdateMaskResolution(settings);
const auto spot_finding_start_time = std::chrono::steady_clock::now();
const std::vector<DiffractionSpot> spots = spotFinder->Run(settings, mask_resolution);
SpotAnalyze(experiment, settings, spots, output);
const auto spot_finding_end_time = std::chrono::steady_clock::now();
output.spot_finding_time_s = std::chrono::duration<float>(spot_finding_end_time - spot_finding_start_time).count();
if (settings.indexing)
indexer.ProcessImage(output, settings,
CompressedImage(updated_image, experiment.GetXPixelsNum(), experiment.GetYPixelsNum()),
*prediction);
}
profile.Add(azim_sum, azim_count);
output.max_viable_pixel_value = max_value;
output.min_viable_pixel_value = min_value;
output.error_pixel_count = err_pixels;
output.saturated_pixel_count = sat_pixels;
output.az_int_profile = profile.GetResult();
output.bkg_estimate = profile.GetBkgEstimate(integration.Settings());
for (const auto &[key, val]: roi_names)
output.roi[key] = roi[val];
}
image_analysis/MXAnalysisWithoutFPGA.h
+2 -13
View File
@@ -15,6 +15,7 @@
#include "spot_finding/ImageSpotFinder.h"
#include "indexing/IndexerThreadPool.h"
#include "IndexAndRefine.h"
#include "image_preprocessing/ImagePreprocessor.h"
// MXAnalysisWithoutFPGA is not thread safe - it has to owned by a single thread
class MXAnalysisWithoutFPGA {
@@ -23,32 +24,20 @@ class MXAnalysisWithoutFPGA {
std::vector<uint8_t> decompression_buffer;
std::vector<uint16_t> roi_map;
std::map<std::string, uint16_t> roi_names;
size_t roi_count;
std::unique_ptr<ImagePreprocessor> preprocessor;
size_t npixels;
size_t xpixels;
std::vector<bool> mask_1bit;
std::unique_ptr<ImageSpotFinder> spotFinder;
IndexAndRefine &indexer;
std::unique_ptr<BraggPrediction> prediction;
std::vector<int32_t> &updated_image;
uint16_t azint_bins;
const int64_t saturation_limit;
const PixelMask &mask;
std::vector<bool> mask_resolution;
float mask_high_res;
float mask_low_res;
void UpdateMaskResolution(const SpotFindingSettings& settings);
template <class T>
void Analyze(DataMessage &output, const uint8_t *image, T err_pixel_val, T sat_pixel_val, AzimuthalIntegrationProfile &profile, const SpotFindingSettings &settings);
public:
MXAnalysisWithoutFPGA(const DiffractionExperiment &experiment, const AzimuthalIntegration &integration,
const PixelMask &mask, IndexAndRefine &indexer);
image_analysis/image_preprocessing/CMakeLists.txt
+5
View File
@@ -0,0 +1,5 @@
ADD_LIBRARY(JFJochImagePreprocessing
STATIC ImagePreprocessorCPU.cpp ImagePreprocessorCPU.h
ImagePreprocessor.cpp ImagePreprocessor.h)
TARGET_LINK_LIBRARIES(JFJochImagePreprocessing JFJochCommon)
image_analysis/image_preprocessing/ImagePreprocessor.cpp
+31
View File
@@ -0,0 +1,31 @@
//
// Created by jungfrau on 4/22/26.
//
#include "ImagePreprocessor.h"
ImagePreprocessor::ImagePreprocessor(const DiffractionExperiment &experiment,
const AzimuthalIntegration &integration,
const PixelMask &mask)
: npixels(experiment.GetPixelsNum()),
experiment(experiment),
integration(integration),
azint_sum(integration.GetBinNumber(), 0.0),
azint_sum2(integration.GetBinNumber(), 0.0),
azint_count(integration.GetBinNumber(), 0),
processed_image(npixels, INT32_MIN),
mask_1bit(npixels, false),
azint_bins(integration.GetBinNumber()),
saturation_limit(experiment.GetSaturationLimit()) {
for (int i = 0; i < npixels; i++)
mask_1bit[i] = (mask.GetMask().at(i) != 0);
}
const std::vector<int32_t> &ImagePreprocessor::GetProcessedImage() const {
return processed_image;
}
void ImagePreprocessor::Update(AzimuthalIntegrationProfile &profile) const {
profile.Clear(integration);
profile.Add(azint_sum, azint_count);
}
image_analysis/image_preprocessing/ImagePreprocessor.h
+47
View File
@@ -0,0 +1,47 @@
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include <cstdint>
#include <cstddef>
#include "../common/CompressedImage.h"
#include "../common/DiffractionExperiment.h"
#include "../common/AzimuthalIntegration.h"
#include "../common/PixelMask.h"
#include "../common/AzimuthalIntegrationProfile.h"
struct ImageStatistics {
size_t error_pixel_count = 0;
size_t saturated_pixel_count = 0;
size_t masked_pixel_count = 0;
int64_t max_value = INT64_MIN;
int64_t min_value = INT64_MAX;
};
class ImagePreprocessor {
protected:
const size_t npixels;
const DiffractionExperiment &experiment;
const AzimuthalIntegration &integration;
std::vector<float> azint_sum;
std::vector<float> azint_sum2;
std::vector<uint32_t> azint_count;
std::vector<int32_t> processed_image;
std::vector<bool> mask_1bit;
uint16_t azint_bins;
const int64_t saturation_limit;
public:
ImagePreprocessor(const DiffractionExperiment &experiment, const AzimuthalIntegration &integration, const PixelMask &mask);
virtual ~ImagePreprocessor() = default;
virtual ImageStatistics Analyze(const uint8_t *decompressed_image, CompressedImageMode image_mode) = 0;
[[nodiscard]] const std::vector<int32_t> &GetProcessedImage() const;
void Update(AzimuthalIntegrationProfile &profile) const;
};
image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp
+77
View File
@@ -0,0 +1,77 @@
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "ImagePreprocessorCPU.h"
ImagePreprocessorCPU::ImagePreprocessorCPU(const DiffractionExperiment &in_experiment,
const AzimuthalIntegration &in_integration,
const PixelMask &in_mask)
: ImagePreprocessor(in_experiment, in_integration, in_mask) {}
ImageStatistics ImagePreprocessorCPU::Analyze(const uint8_t *image_ptr, CompressedImageMode image_mode) {
switch (image_mode) {
case CompressedImageMode::Int8:
return Analyze<int8_t>(image_ptr, INT8_MIN, INT8_MAX);
case CompressedImageMode::Int16:
return Analyze<int16_t>(image_ptr, INT16_MIN, INT16_MAX);
case CompressedImageMode::Int32:
return Analyze<int32_t>(image_ptr, INT32_MIN, INT32_MAX);
case CompressedImageMode::Uint8:
return Analyze<uint8_t>(image_ptr, UINT8_MAX, UINT8_MAX);
case CompressedImageMode::Uint16:
return Analyze<uint16_t>(image_ptr, UINT16_MAX, UINT16_MAX);
case CompressedImageMode::Uint32:
return Analyze<uint32_t>(image_ptr, UINT32_MAX, UINT32_MAX);
default:
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "RGB/float mode not supported");
}
}
template<class T>
ImageStatistics ImagePreprocessorCPU::Analyze(const uint8_t *input, T err_pixel_val, T sat_pixel_val) {
auto image = reinterpret_cast<const T *>(input);
for (int i = 0; i < azint_count.size(); i++) {
azint_sum[i] = 0.0f;
azint_sum2[i] = 0.0f;
azint_count[i] = 0;
}
ImageStatistics ret{};
if (sat_pixel_val > saturation_limit)
sat_pixel_val = static_cast<T>(saturation_limit);
auto &pixel_to_bin = integration.GetPixelToBin();
auto &corrections = integration.Corrections();
for (int i = 0; i < npixels; i++) {
if (mask_1bit[i] != 0) {
processed_image[i] = INT32_MIN;
++ret.masked_pixel_count;
} else if (image[i] >= sat_pixel_val) {
processed_image[i] = INT32_MAX;
++ret.saturated_pixel_count;
} else if (std::is_signed<T>::value && (image[i] == err_pixel_val)) {
// Error pixels are possible only for signed types
processed_image[i] = INT32_MIN;
++ret.error_pixel_count;
} else {
processed_image[i] = static_cast<int32_t>(image[i]);
if (image[i] > ret.max_value)
ret.max_value = image[i];
if (image[i] < ret.min_value)
ret.min_value = image[i];
const uint16_t bin = pixel_to_bin[i];
if (bin < azint_bins) {
float val = image[i] * corrections[i];
azint_sum[bin] += val;
azint_sum2[bin] += val * val;
++azint_count[bin];
}
}
}
return ret;
}
image_analysis/image_preprocessing/ImagePreprocessorCPU.h
+24
View File
@@ -0,0 +1,24 @@
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include <vector>
#include "ImagePreprocessor.h"
#include "../common/JFJochMessages.h"
#include "../common/DiffractionExperiment.h"
#include "../common/AzimuthalIntegration.h"
#include "../common/PixelMask.h"
#include "../common/AzimuthalIntegrationProfile.h"
class ImagePreprocessorCPU : public ImagePreprocessor {
template <class T>
ImageStatistics Analyze(const uint8_t *input, T err_value, T sat_value);
public:
ImagePreprocessorCPU(const DiffractionExperiment &in_experiment,
const AzimuthalIntegration &in_integration,
const PixelMask &in_mask);
ImageStatistics Analyze(const uint8_t *decompressed_image, CompressedImageMode image_mode) override;
};