v1.0.0-rc.137 (#46)

This is an UNSTABLE release. The release has significant modifications and bug fixes, if things go wrong, it is better to revert to 1.0.0-rc.132. * jfjoch_broker: Better track time for each operation in the processing stack * jfjoch_broker: Rewrite preprocessing of diffraction images in the non-FPGA workflow to better use GPUs (work in progress) * jfjoch_broker: Remove ROI calculation in the non-FPGA workflow (work in progress) * jfjoch_viewer: Toolbar displays image number starting from 1 (instead of 0) Reviewed-on: #46
2026-04-25 19:59:21 +02:00
parent c1c170112c
commit c981e1b91c
184 changed files with 1190 additions and 486 deletions
@@ -7,8 +7,19 @@
 #include "../compression/JFJochDecompress.h"

 #include "spot_finding/SpotUtils.h"
-#include "spot_finding/ImageSpotFinderFactory.h"
 #include "bragg_prediction/BraggPredictionFactory.h"
+#include "image_preprocessing/ImagePreprocessorCPU.h"
+
+#include "azint/AzIntEngineCPU.h"
+#include "spot_finding/ImageSpotFinderCPU.h"
+#ifdef JFJOCH_USE_CUDA
+#include "azint/AzIntEngineGPU.h"
+#include "spot_finding/ImageSpotFinderGPU.h"
+#include "image_preprocessing/ImagePreprocessorGPU.h"
+#include "image_preprocessing/ImagePreprocessorBufferGPU.h"
+#include "../common/CUDAWrapper.h"
+#endif
+

 MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_experiment,
                                             const AzimuthalIntegration &in_integration,
@@ -16,27 +27,33 @@ 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);
+#ifdef JFJOCH_USE_CUDA
+    if (get_gpu_count() == 0) {
+#endif
+        preprocessor_buffer = std::make_unique<ImagePreprocessorBuffer>(experiment.GetPixelsNum());
+        spotFinder = std::make_unique<ImageSpotFinderCPU>(experiment.GetXPixelsNum(), experiment.GetYPixelsNum());
+        azint = std::make_unique<AzIntEngineCPU>(integration);
+        preprocessor = std::make_unique<ImagePreprocessorCPU>(in_experiment, in_mask);
+#ifdef JFJOCH_USE_CUDA
+    } else {
+        auto stream = std::make_shared<CudaStream>();
+        preprocessor_buffer = std::make_unique<ImagePreprocessorBufferGPU>(experiment.GetPixelsNum());
+        preprocessor = std::make_unique<ImagePreprocessorGPU>(in_experiment, in_mask, stream);
+        spotFinder = std::make_unique<ImageSpotFinderGPU>(experiment.GetXPixelsNum(), experiment.GetYPixelsNum(), stream);
+        azint = std::make_unique<AzIntEngineGPU>(integration, stream);
+    }
+#endif
 }

-void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, std::vector<uint8_t> &image,
+void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
                                    AzimuthalIntegrationProfile &profile,
                                    const SpotFindingSettings &spot_finding_settings) {
    if ((output.image.GetWidth() != xpixels)
@@ -44,30 +61,46 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, std::vector<uint8_t> &i
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Mismatch in pixel size");

-    const uint8_t *image_ptr = output.image.GetUncompressedPtr(image);
+    const auto compression_start_time = std::chrono::steady_clock::now();
+    const uint8_t *image_ptr = output.image.GetUncompressedPtr(decompression_buffer);
+    const auto compression_end_time = std::chrono::steady_clock::now();
+    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(*preprocessor_buffer, 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();
+
+    const auto azint_start_time = std::chrono::steady_clock::now();
+    azint->Run(*preprocessor_buffer, profile);
+    const auto azint_end_time = std::chrono::steady_clock::now();
+    output.azint_time_s = std::chrono::duration<float>(azint_end_time - azint_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();
+        const std::vector<DiffractionSpot> spots = spotFinder->Run(*preprocessor_buffer, 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_buffer->getBuffer(), experiment.GetXPixelsNum(), experiment.GetYPixelsNum()),
+                *prediction);
    }
+
+    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) {
@@ -77,108 +110,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) {
-    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];
-            }
-        }
-    }
-
-    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];
-}