diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt
index 08054bb4..ed9dca95 100644
--- a/image_analysis/CMakeLists.txt
+++ b/image_analysis/CMakeLists.txt
@@ -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)
diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp
index def9ab24..094f1762 100644
--- a/image_analysis/MXAnalysisWithoutFPGA.cpp
+++ b/image_analysis/MXAnalysisWithoutFPGA.cpp
@@ -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];
-}
diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h
index 454a4209..21d4d06a 100644
--- a/image_analysis/MXAnalysisWithoutFPGA.h
+++ b/image_analysis/MXAnalysisWithoutFPGA.h
@@ -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);
diff --git a/image_analysis/image_preprocessing/CMakeLists.txt b/image_analysis/image_preprocessing/CMakeLists.txt
new file mode 100644
index 00000000..fd796a7e
--- /dev/null
+++ b/image_analysis/image_preprocessing/CMakeLists.txt
@@ -0,0 +1,5 @@
+ADD_LIBRARY(JFJochImagePreprocessing
+        STATIC ImagePreprocessorCPU.cpp ImagePreprocessorCPU.h
+        ImagePreprocessor.cpp ImagePreprocessor.h)
+
+TARGET_LINK_LIBRARIES(JFJochImagePreprocessing JFJochCommon)
\ No newline at end of file
diff --git a/image_analysis/image_preprocessing/ImagePreprocessor.cpp b/image_analysis/image_preprocessing/ImagePreprocessor.cpp
new file mode 100644
index 00000000..dd4413de
--- /dev/null
+++ b/image_analysis/image_preprocessing/ImagePreprocessor.cpp
@@ -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);
+}
diff --git a/image_analysis/image_preprocessing/ImagePreprocessor.h b/image_analysis/image_preprocessing/ImagePreprocessor.h
new file mode 100644
index 00000000..61891a66
--- /dev/null
+++ b/image_analysis/image_preprocessing/ImagePreprocessor.h
@@ -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;
+};
diff --git a/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp b/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp
new file mode 100644
index 00000000..af513d0c
--- /dev/null
+++ b/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp
@@ -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;
+}
diff --git a/image_analysis/image_preprocessing/ImagePreprocessorCPU.h b/image_analysis/image_preprocessing/ImagePreprocessorCPU.h
new file mode 100644
index 00000000..4f8d5b7c
--- /dev/null
+++ b/image_analysis/image_preprocessing/ImagePreprocessorCPU.h
@@ -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;
+};