Merge branch 'dev240207-resnet' into 'main'

Add deep learning resolution estimation model from Stanford

See merge request jungfraujoch/nextgendcu!26

CMakeLists.txt

@@ -27,6 +27,13 @@ SET(JFJOCH_WRITER_ONLY    OFF CACHE BOOL "Compile HDF5 writer only")
 INCLUDE_DIRECTORIES(include)
 INCLUDE(CheckIncludeFile)
 
+#This is to supress error in TORCH
+
+IF ((NOT EXISTS /usr/bin/python) AND (EXISTS /usr/bin/python3))
+    SET(PYTHON_EXECUTABLE /usr/bin/python3)
+ENDIF()
+
+FIND_PACKAGE(Torch HINTS /opt/libtorch/share/cmake/Torch/)
 FIND_LIBRARY(NUMA_LIBRARY NAMES numa DOC "NUMA Library")
 CHECK_INCLUDE_FILE(numaif.h HAS_NUMAIF)
 CHECK_INCLUDE_FILE(numa.h HAS_NUMA_H)
@@ -54,6 +61,7 @@ ELSE()
     ADD_SUBDIRECTORY(tests)
     ADD_SUBDIRECTORY(tools)
     ADD_SUBDIRECTORY(export_images)
+    ADD_SUBDIRECTORY(resonet)
     SET(jfjoch_executables jfjoch_broker jfjoch_writer CatchTest CompressionBenchmark HDF5DatasetWriteTest jfjoch_udp_simulator sls_detector_put sls_detector_get)
 ENDIF()
 

README.md

@@ -38,6 +38,7 @@ Optional:
 * CUDA compiler version 11 or newer - image analysis features won't work without it
 * NUMA library
 * Node.js - to make frontend
+* libtorch - for resolution estimation using model from Stanford - see below
 
 provided as GIT submodules: SLS Detector Package, tinycbor (Intel), Zstandard (Facebook), Pistache webserver and fast replacement for bitshuffle filter (DECTRIS). 
 
@@ -90,5 +91,16 @@ Automated test routine is then accessible as `tests/CatchTest`. There are also b
 
 In addition, tests are executed to verify that datasets written by Jungfraujoch are readable with XDS Durin plugin and CrystFEL. Input files for these programs are placed in `xds_durin` and `crystfel` folders. See `.gitlab-ci.yml` for details.
 
+## Resolution estimation
 
+Resolution estimation can be done with a recent deep learning model by D. Mendez et al. 
+(see [Acta Cryst D, 80, 26-43](https://doi.org/10.1107/S2059798323010586)), adapted to Jungfraujoch. 
+Model used in the original paper is located in the [resonet/](resonet) directory, after converting to TorchScript format.
+
+To use the feature it is necessary to install [libtorch](https://pytorch.org/) library, preferably in `/opt/libtorch` location. The C++11 ABI version needs to be chosen.
+For RHEL 8 systems, please download older version 2.1.0, as version 2.2.0 requires newer `glibc` library than available with the operating system. 
+Version 2.1.0 can be downloaded with the following command:
+```
+wget https://download.pytorch.org/libtorch/cu121/libtorch-cxx11-abi-shared-with-deps-2.1.0%2Bcu121.zip
+```
 

acquisition_device/CMakeLists.txt

@@ -1,4 +1,4 @@
-ADD_LIBRARY(JungfraujochAcqusitionDevice STATIC
+ADD_LIBRARY(JFJochAcquisitionDevice STATIC
         AcquisitionDevice.cpp AcquisitionDevice.h
         AcquisitionCounters.cpp AcquisitionCounters.h
         HLSSimulatedDevice.cpp HLSSimulatedDevice.h
@@ -8,4 +8,4 @@ ADD_LIBRARY(JungfraujochAcqusitionDevice STATIC
         AcquisitionDeviceGroup.cpp
         AcquisitionDeviceGroup.h)
 
-TARGET_LINK_LIBRARIES(JungfraujochAcqusitionDevice JungfraujochDevice CommonFunctions HLSSimulation JFCalibration)
+TARGET_LINK_LIBRARIES(JFJochAcquisitionDevice JFJochDevice JFJochCommon JFJochHLSSimulation JFCalibration)

broker/CMakeLists.txt

@@ -9,7 +9,7 @@ ADD_LIBRARY(JFJochBroker STATIC
         JFJochServices.cpp JFJochServices.h
         JFJochBrokerHttp.cpp JFJochBrokerHttp.h JFJochBrokerParser.cpp JFJochBrokerParser.h)
 
-TARGET_LINK_LIBRARIES(JFJochBroker JFJochReceiver JFJochDetector CommonFunctions JFJochAPI)
+TARGET_LINK_LIBRARIES(JFJochBroker JFJochReceiver JFJochDetector JFJochCommon JFJochAPI)
 
 ADD_EXECUTABLE(jfjoch_broker jfjoch_broker.cpp)
 TARGET_LINK_LIBRARIES(jfjoch_broker JFJochBroker)

broker/JFJochBrokerHttp.cpp

@@ -703,3 +703,9 @@ void JFJochBrokerHttp::preview_jpeg_settings_put(const org::openapitools::server
     state_machine.SetPreviewJPEGSettings(Convert(previewSettings));
     response.send(Pistache::Http::Code::Ok);
 }
+
+void JFJochBrokerHttp::plot_resolution_estimate_histogram_get(Pistache::Http::ResponseWriter &response) {
+    PlotRequest req{.type = PlotType::ResEstimation, .binning = 0};
+    auto plot = state_machine.GetPlots(req);
+    ProcessOutput(Convert(plot), response);
+}

broker/JFJochBrokerHttp.h

@@ -71,6 +71,8 @@ class JFJochBrokerHttp : public org::openapitools::server::api::DefaultApi {
     void plot_spot_count_post(const org::openapitools::server::model::Plot_request &plotRequest,
                               Pistache::Http::ResponseWriter &response) override;
 
+    void plot_resolution_estimate_histogram_get(Pistache::Http::ResponseWriter &response) override;
+
     void statistics_calibration_get(Pistache::Http::ResponseWriter &response) override;
 
     void statistics_data_collection_get(Pistache::Http::ResponseWriter &response) override;
@@ -102,6 +104,7 @@ class JFJochBrokerHttp : public org::openapitools::server::api::DefaultApi {
     void preview_jpeg_settings_put(const org::openapitools::server::model::Preview_settings &previewSettings,
                                    Pistache::Http::ResponseWriter &response) override;
 
+
     void GetStaticFile(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
     std::pair<Pistache::Http::Code, std::string> handleOperationException(const std::exception &ex) const noexcept override;
 

broker/JFJochBrokerParser.cpp

@@ -273,6 +273,7 @@ void ParseFacilityConfiguration(const nlohmann::json &input, const std::string&
                                       "omega_axis must be float array of 3");
         }
 
+        experiment.NeuralNetModelPath(GET_STR(j, "neural_net_model", ""));
         experiment.PedestalG0Frames(GET_I64(j, "pedestal_g0_frames"));
         experiment.PedestalG1Frames(GET_I64(j, "pedestal_g1_frames"));
         experiment.PedestalG2Frames(GET_I64(j, "pedestal_g2_frames"));

broker/gen/api/DefaultApi.cpp

@@ -54,6 +54,7 @@ void DefaultApi::setupRoutes() {
     Routes::Get(*router, base + "/plot/rad_int", Routes::bind(&DefaultApi::plot_rad_int_get_handler, this));
     Routes::Get(*router, base + "/plot/rad_int_per_file", Routes::bind(&DefaultApi::plot_rad_int_per_file_get_handler, this));
     Routes::Post(*router, base + "/plot/receiver_delay", Routes::bind(&DefaultApi::plot_receiver_delay_post_handler, this));
+    Routes::Get(*router, base + "/plot/resolution_estimate_histogram", Routes::bind(&DefaultApi::plot_resolution_estimate_histogram_get_handler, this));
     Routes::Post(*router, base + "/plot/roi_sum", Routes::bind(&DefaultApi::plot_roi_sum_post_handler, this));
     Routes::Post(*router, base + "/plot/saturated_pixel", Routes::bind(&DefaultApi::plot_saturated_pixel_post_handler, this));
     Routes::Post(*router, base + "/plot/spot_count", Routes::bind(&DefaultApi::plot_spot_count_post_handler, this));
@@ -628,6 +629,26 @@ void DefaultApi::plot_receiver_delay_post_handler(const Pistache::Rest::Request
         response.send(Pistache::Http::Code::Internal_Server_Error, e.what());
     }
 
+}
+void DefaultApi::plot_resolution_estimate_histogram_get_handler(const Pistache::Rest::Request &, Pistache::Http::ResponseWriter response) {
+    try {
+
+
+    try {
+        this->plot_resolution_estimate_histogram_get(response);
+    } catch (Pistache::Http::HttpError &e) {
+        response.send(static_cast<Pistache::Http::Code>(e.code()), e.what());
+        return;
+    } catch (std::exception &e) {
+        const std::pair<Pistache::Http::Code, std::string> errorInfo = this->handleOperationException(e);
+        response.send(errorInfo.first, errorInfo.second);
+        return;
+    }
+
+    } catch (std::exception &e) {
+        response.send(Pistache::Http::Code::Internal_Server_Error, e.what());
+    }
+
 }
 void DefaultApi::plot_roi_sum_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response) {
     try {

broker/gen/api/DefaultApi.h

@@ -78,6 +78,7 @@ private:
     void plot_rad_int_get_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
     void plot_rad_int_per_file_get_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
     void plot_receiver_delay_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
+    void plot_resolution_estimate_histogram_get_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
     void plot_roi_sum_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
     void plot_saturated_pixel_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
     void plot_spot_count_post_handler(const Pistache::Rest::Request &request, Pistache::Http::ResponseWriter response);
@@ -268,6 +269,13 @@ private:
     /// <param name="plotRequest"> (optional)</param>
     virtual void plot_receiver_delay_post(const org::openapitools::server::model::Plot_request &plotRequest, Pistache::Http::ResponseWriter &response) = 0;
     /// <summary>
+    /// Generate resolution estimate histogram
+    /// </summary>
+    /// <remarks>
+    /// Generate histogram of crystal resolutions from 1.0 to 5.0 A based on ML model
+    /// </remarks>
+    virtual void plot_resolution_estimate_histogram_get(Pistache::Http::ResponseWriter &response) = 0;
+    /// <summary>
     /// Generate ROI sum plot
     /// </summary>
     /// <remarks>

broker/gen/model/Detector_status.cpp

@@ -22,14 +22,9 @@ namespace org::openapitools::server::model
 Detector_status::Detector_status()
 {
     m_State = "";
-    m_StateIsSet = false;
     m_Powerchip = "";
-    m_PowerchipIsSet = false;
     m_Server_version = "";
-    m_Server_versionIsSet = false;
     m_Number_of_triggers_left = 0L;
-    m_Number_of_triggers_leftIsSet = false;
-    m_Fpga_temp_degCIsSet = false;
     
 }
 
@@ -53,8 +48,8 @@ bool Detector_status::validate(std::stringstream& msg, const std::string& pathPr
     const std::string _pathPrefix = pathPrefix.empty() ? "Detector_status" : pathPrefix;
 
                          
-    if (fpgaTempDegCIsSet())
-    {
+    
+    /* Fpga_temp_degC */ {
         const std::vector<int64_t>& value = m_Fpga_temp_degC;
         const std::string currentValuePath = _pathPrefix + ".fpgaTempDegC";
                 
@@ -82,20 +77,20 @@ bool Detector_status::operator==(const Detector_status& rhs) const
     return
     
     
+    (getState() == rhs.getState())
+     &&
     
-    ((!stateIsSet() && !rhs.stateIsSet()) || (stateIsSet() && rhs.stateIsSet() && getState() == rhs.getState())) &&
+    (getPowerchip() == rhs.getPowerchip())
+     &&
     
+    (getServerVersion() == rhs.getServerVersion())
+     &&
     
-    ((!powerchipIsSet() && !rhs.powerchipIsSet()) || (powerchipIsSet() && rhs.powerchipIsSet() && getPowerchip() == rhs.getPowerchip())) &&
+    (getNumberOfTriggersLeft() == rhs.getNumberOfTriggersLeft())
+     &&
     
+    (getFpgaTempDegC() == rhs.getFpgaTempDegC())
     
-    ((!serverVersionIsSet() && !rhs.serverVersionIsSet()) || (serverVersionIsSet() && rhs.serverVersionIsSet() && getServerVersion() == rhs.getServerVersion())) &&
-    
-    
-    ((!numberOfTriggersLeftIsSet() && !rhs.numberOfTriggersLeftIsSet()) || (numberOfTriggersLeftIsSet() && rhs.numberOfTriggersLeftIsSet() && getNumberOfTriggersLeft() == rhs.getNumberOfTriggersLeft())) &&
-    
-    
-    ((!fpgaTempDegCIsSet() && !rhs.fpgaTempDegCIsSet()) || (fpgaTempDegCIsSet() && rhs.fpgaTempDegCIsSet() && getFpgaTempDegC() == rhs.getFpgaTempDegC()))
     
     ;
 }
@@ -108,46 +103,21 @@ bool Detector_status::operator!=(const Detector_status& rhs) const
 void to_json(nlohmann::json& j, const Detector_status& o)
 {
     j = nlohmann::json();
-    if(o.stateIsSet())
-        j["state"] = o.m_State;
-    if(o.powerchipIsSet())
-        j["powerchip"] = o.m_Powerchip;
-    if(o.serverVersionIsSet())
-        j["server_version"] = o.m_Server_version;
-    if(o.numberOfTriggersLeftIsSet())
-        j["number_of_triggers_left"] = o.m_Number_of_triggers_left;
-    if(o.fpgaTempDegCIsSet() || !o.m_Fpga_temp_degC.empty())
-        j["fpga_temp_degC"] = o.m_Fpga_temp_degC;
+    j["state"] = o.m_State;
+    j["powerchip"] = o.m_Powerchip;
+    j["server_version"] = o.m_Server_version;
+    j["number_of_triggers_left"] = o.m_Number_of_triggers_left;
+    j["fpga_temp_degC"] = o.m_Fpga_temp_degC;
     
 }
 
 void from_json(const nlohmann::json& j, Detector_status& o)
 {
-    if(j.find("state") != j.end())
-    {
-        j.at("state").get_to(o.m_State);
-        o.m_StateIsSet = true;
-    } 
-    if(j.find("powerchip") != j.end())
-    {
-        j.at("powerchip").get_to(o.m_Powerchip);
-        o.m_PowerchipIsSet = true;
-    } 
-    if(j.find("server_version") != j.end())
-    {
-        j.at("server_version").get_to(o.m_Server_version);
-        o.m_Server_versionIsSet = true;
-    } 
-    if(j.find("number_of_triggers_left") != j.end())
-    {
-        j.at("number_of_triggers_left").get_to(o.m_Number_of_triggers_left);
-        o.m_Number_of_triggers_leftIsSet = true;
-    } 
-    if(j.find("fpga_temp_degC") != j.end())
-    {
-        j.at("fpga_temp_degC").get_to(o.m_Fpga_temp_degC);
-        o.m_Fpga_temp_degCIsSet = true;
-    } 
+    j.at("state").get_to(o.m_State);
+    j.at("powerchip").get_to(o.m_Powerchip);
+    j.at("server_version").get_to(o.m_Server_version);
+    j.at("number_of_triggers_left").get_to(o.m_Number_of_triggers_left);
+    j.at("fpga_temp_degC").get_to(o.m_Fpga_temp_degC);
     
 }
 
@@ -158,15 +128,6 @@ std::string Detector_status::getState() const
 void Detector_status::setState(std::string const& value)
 {
     m_State = value;
-    m_StateIsSet = true;
-}
-bool Detector_status::stateIsSet() const
-{
-    return m_StateIsSet;
-}
-void Detector_status::unsetState()
-{
-    m_StateIsSet = false;
 }
 std::string Detector_status::getPowerchip() const
 {
@@ -175,15 +136,6 @@ std::string Detector_status::getPowerchip() const
 void Detector_status::setPowerchip(std::string const& value)
 {
     m_Powerchip = value;
-    m_PowerchipIsSet = true;
-}
-bool Detector_status::powerchipIsSet() const
-{
-    return m_PowerchipIsSet;
-}
-void Detector_status::unsetPowerchip()
-{
-    m_PowerchipIsSet = false;
 }
 std::string Detector_status::getServerVersion() const
 {
@@ -192,15 +144,6 @@ std::string Detector_status::getServerVersion() const
 void Detector_status::setServerVersion(std::string const& value)
 {
     m_Server_version = value;
-    m_Server_versionIsSet = true;
-}
-bool Detector_status::serverVersionIsSet() const
-{
-    return m_Server_versionIsSet;
-}
-void Detector_status::unsetServer_version()
-{
-    m_Server_versionIsSet = false;
 }
 int64_t Detector_status::getNumberOfTriggersLeft() const
 {
@@ -209,15 +152,6 @@ int64_t Detector_status::getNumberOfTriggersLeft() const
 void Detector_status::setNumberOfTriggersLeft(int64_t const value)
 {
     m_Number_of_triggers_left = value;
-    m_Number_of_triggers_leftIsSet = true;
-}
-bool Detector_status::numberOfTriggersLeftIsSet() const
-{
-    return m_Number_of_triggers_leftIsSet;
-}
-void Detector_status::unsetNumber_of_triggers_left()
-{
-    m_Number_of_triggers_leftIsSet = false;
 }
 std::vector<int64_t> Detector_status::getFpgaTempDegC() const
 {
@@ -226,15 +160,6 @@ std::vector<int64_t> Detector_status::getFpgaTempDegC() const
 void Detector_status::setFpgaTempDegC(std::vector<int64_t> const value)
 {
     m_Fpga_temp_degC = value;
-    m_Fpga_temp_degCIsSet = true;
-}
-bool Detector_status::fpgaTempDegCIsSet() const
-{
-    return m_Fpga_temp_degCIsSet;
-}
-void Detector_status::unsetFpga_temp_degC()
-{
-    m_Fpga_temp_degCIsSet = false;
 }
 
 

broker/gen/model/Detector_status.h

@@ -64,50 +64,40 @@ public:
     /// </summary>
     std::string getState() const;
     void setState(std::string const& value);
-    bool stateIsSet() const;
-    void unsetState();
     /// <summary>
     /// Power on of ASICs
     /// </summary>
     std::string getPowerchip() const;
     void setPowerchip(std::string const& value);
-    bool powerchipIsSet() const;
-    void unsetPowerchip();
     /// <summary>
     /// Detector server (on read-out boards) version
     /// </summary>
     std::string getServerVersion() const;
     void setServerVersion(std::string const& value);
-    bool serverVersionIsSet() const;
-    void unsetServer_version();
     /// <summary>
     /// Remaining triggers to the detector (max of all modules)
     /// </summary>
     int64_t getNumberOfTriggersLeft() const;
     void setNumberOfTriggersLeft(int64_t const value);
-    bool numberOfTriggersLeftIsSet() const;
-    void unsetNumber_of_triggers_left();
     /// <summary>
     /// Temperature of detector FPGAs
     /// </summary>
     std::vector<int64_t> getFpgaTempDegC() const;
     void setFpgaTempDegC(std::vector<int64_t> const value);
-    bool fpgaTempDegCIsSet() const;
-    void unsetFpga_temp_degC();
 
     friend void to_json(nlohmann::json& j, const Detector_status& o);
     friend void from_json(const nlohmann::json& j, Detector_status& o);
 protected:
     std::string m_State;
-    bool m_StateIsSet;
+
     std::string m_Powerchip;
-    bool m_PowerchipIsSet;
+
     std::string m_Server_version;
-    bool m_Server_versionIsSet;
+
     int64_t m_Number_of_triggers_left;
-    bool m_Number_of_triggers_leftIsSet;
+
     std::vector<int64_t> m_Fpga_temp_degC;
-    bool m_Fpga_temp_degCIsSet;
+
     
 };
 

broker/jfjoch_api.yaml

@@ -1131,6 +1131,17 @@ paths:
             application/json:
               schema:
                 $ref: '#/components/schemas/plot'
+  /plot/resolution_estimate_histogram:
+    get:
+      summary: Generate resolution estimate histogram
+      description: Generate histogram of crystal resolutions from 1.0 to 5.0 A based on ML model
+      responses:
+        "200":
+          description: Everything OK
+          content:
+            application/json:
+              schema:
+                $ref: '#/components/schemas/plot'
   /plot/rad_int:
     get:
       summary: Generate radial integration profile

common/CMakeLists.txt

@@ -21,7 +21,7 @@ MESSAGE(STATUS "Git date: ${GIT_DATE}")
 
 CONFIGURE_FILE("${CMAKE_CURRENT_SOURCE_DIR}/GitInfo.cpp.in" "${CMAKE_CURRENT_BINARY_DIR}/GitInfo.cpp" @ONLY)
 
-ADD_LIBRARY( CommonFunctions STATIC
+ADD_LIBRARY( JFJochCommon STATIC
         Logger.cpp Logger.h
         Coord.cpp Coord.h
         DiffractionExperiment.cpp DiffractionExperiment.h
@@ -45,16 +45,16 @@ ADD_LIBRARY( CommonFunctions STATIC
         ../fpga/include/jfjoch_fpga.h
         ZMQWrappers.cpp ZMQWrappers.h)
 
-TARGET_LINK_LIBRARIES(CommonFunctions Compression JFCalibration "$<BUILD_INTERFACE:libzmq-static>" -lrt)
+TARGET_LINK_LIBRARIES(JFJochCommon Compression JFCalibration "$<BUILD_INTERFACE:libzmq-static>" -lrt)
 
 IF (CMAKE_CUDA_COMPILER)
-    TARGET_SOURCES(CommonFunctions PRIVATE CUDAWrapper.cu )
-    TARGET_LINK_LIBRARIES(CommonFunctions ${CUDART_LIBRARY} ${CMAKE_DL_LIBS} rt)
+    TARGET_SOURCES(JFJochCommon PRIVATE CUDAWrapper.cu )
+    TARGET_LINK_LIBRARIES(JFJochCommon ${CUDART_LIBRARY} ${CMAKE_DL_LIBS} rt)
 ENDIF()
 
 IF(HAS_NUMAIF AND HAS_NUMA_H AND NUMA_LIBRARY)
-    TARGET_COMPILE_DEFINITIONS(CommonFunctions PUBLIC JFJOCH_USE_NUMA)
-    TARGET_LINK_LIBRARIES(CommonFunctions ${NUMA_LIBRARY})
+    TARGET_COMPILE_DEFINITIONS(JFJochCommon PUBLIC JFJOCH_USE_NUMA)
+    TARGET_LINK_LIBRARIES(JFJochCommon ${NUMA_LIBRARY})
     MESSAGE(STATUS "NUMA memory/CPU pinning enabled")
 ELSE()
     MESSAGE(WARNING "NUMA memory/CPU pinning disabled")

common/DiffractionExperiment.cpp

@@ -1147,3 +1147,12 @@ std::chrono::nanoseconds DiffractionExperiment::GetDetectorDelay() const {
 const DetectorSetup &DiffractionExperiment::GetDetectorSetup() const {
     return detector;
 }
+
+DiffractionExperiment &DiffractionExperiment::NeuralNetModelPath(const std::string &input) {
+    neural_net_model_path = input;
+    return *this;
+}
+
+std::string DiffractionExperiment::GetNeuralNetModelPath() const {
+    return neural_net_model_path;
+}

common/DiffractionExperiment.h

@@ -91,6 +91,8 @@ class DiffractionExperiment {
 
     uint64_t series_id;
 
+    std::string neural_net_model_path;
+
     // Dataset settings
     struct {
         int64_t images_per_trigger;
@@ -197,6 +199,7 @@ public:
     DiffractionExperiment& IncrementSeriesID();
     DiffractionExperiment& ROISummation(const std::optional<ROIRectangle>& input);
     DiffractionExperiment& ConversionOnFPGA(bool input);
+    DiffractionExperiment& NeuralNetModelPath(const std::string& input);
 
     void FillMessage(StartMessage &message) const;
 
@@ -343,6 +346,7 @@ public:
     bool IsConversionOnFPGA() const;
 
     const DetectorSetup& GetDetectorSetup() const;
+    std::string GetNeuralNetModelPath() const;
 };
 
 inline int64_t CalculateStride(const std::chrono::microseconds &frame_time, const std::chrono::microseconds &preview_time) {

common/Histogram.h

@@ -7,15 +7,17 @@
 #include <cstddef>
 #include <vector>
 #include <mutex>
+
 #include "Plot.h"
+#include "../common/JFJochException.h"
 
 template <class T>
-class Histogram {
+class SetAverage {
     std::vector<T> sum;
     std::vector<uint64_t> count;
     std::mutex m;
 public:
-    explicit Histogram(size_t bins) : sum(bins), count(bins) {}
+    explicit SetAverage(size_t bins) : sum(bins), count(bins) {}
     void Add(size_t bin, T val) {
         std::unique_lock<std::mutex> ul(m);
 
@@ -42,5 +44,39 @@ public:
     }
 };
 
+class FloatHistogram {
+    std::vector<uint64_t> count;
+    std::mutex m;
+    float min;
+    float max;
+    float spacing;
+public:
+    explicit FloatHistogram(size_t bins, float in_min, float in_max) : count(bins), min(in_min), max(in_max) {
+        if (bins == 0)
+            throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, "FloatHistogram neds to have non-zero bins");
+        spacing = (max - min) / bins;
+    }
+    void Add(float val) {
+        std::unique_lock<std::mutex> ul(m);
+
+        size_t bin = (val - min) / spacing;
+
+        if (bin < count.size())
+            count[bin] += 1;
+    }
+
+    Plot GetPlot() {
+        std::unique_lock<std::mutex> ul(m);
+
+        Plot ret;
+        ret.x.resize(count.size());
+        ret.y.resize(count.size());
+        for (int i = 0; i < count.size(); i++) {
+            ret.x[i] = min + (i + 0.5f) * spacing;
+            ret.y[i] = static_cast<float>(count[i]);
+        }
+        return ret;
+    }
+};
 
 #endif //JUNGFRAUJOCH_HISTOGRAM_H

common/Plot.h

@@ -8,7 +8,8 @@
 #include <string>
 
 enum class PlotType {BkgEstimate, RadInt, SpotCount, IndexingRate, IndexingRatePerFile, SaturatedPixels,
-                     ErrorPixels, ImageCollectionEfficiency, ReceiverDelay, StrongPixels, ROISum};
+                     ErrorPixels, ImageCollectionEfficiency, ReceiverDelay, StrongPixels, ROISum,
+                     ResEstimation};
 
 struct PlotRequest {
     PlotType type;

detector_control/CMakeLists.txt

@@ -2,6 +2,6 @@ ADD_SUBDIRECTORY(slsDetectorPackage)
 INSTALL(TARGETS sls_detector_put sls_detector_get RUNTIME)
 
 ADD_LIBRARY(JFJochDetector STATIC DetectorWrapper.cpp DetectorWrapper.h)
-TARGET_LINK_LIBRARIES(JFJochDetector CommonFunctions slsSupportShared slsDetectorShared)
+TARGET_LINK_LIBRARIES(JFJochDetector JFJochCommon slsSupportShared slsDetectorShared)
 
 

export_images/CMakeLists.txt

@@ -1,24 +1,24 @@
 FIND_PACKAGE(TIFF COMPONENTS CXX REQUIRED)
 FIND_PACKAGE(JPEG REQUIRED)
 
-ADD_LIBRARY(ExportImage STATIC
+ADD_LIBRARY(JFJochImageExport STATIC
         WriteTIFF.cpp WriteTIFF.h
         WriteJPEG.cpp WriteJPEG.h
         PreviewImage.cpp PreviewImage.h)
 
-TARGET_LINK_LIBRARIES(ExportImage PUBLIC CommonFunctions)
+TARGET_LINK_LIBRARIES(JFJochImageExport PUBLIC JFJochCommon)
 
 IF((EXISTS ${TIFF_INCLUDE_DIR}/tiffio.hxx) AND (EXISTS ${TIFF_INCLUDE_DIR}/tiffio.h))
-    TARGET_INCLUDE_DIRECTORIES(ExportImage PRIVATE ${TIFF_INCLUDE_DIR})
-    TARGET_LINK_LIBRARIES(ExportImage PUBLIC ${TIFF_LIBRARIES})
+    TARGET_INCLUDE_DIRECTORIES(JFJochImageExport PRIVATE ${TIFF_INCLUDE_DIR})
+    TARGET_LINK_LIBRARIES(JFJochImageExport PUBLIC ${TIFF_LIBRARIES})
     MESSAGE(STATUS "TIFF headers present and library included")
 ELSE()
     MESSAGE(FATAL_ERROR "TIFF headers tiffio.h and tiffio.hxx not present")
 ENDIF()
 
 IF (EXISTS ${JPEG_INCLUDE_DIR}/jpeglib.h)
-    TARGET_INCLUDE_DIRECTORIES(ExportImage PRIVATE ${JPEG_INCLUDE_DIR})
-    TARGET_LINK_LIBRARIES(ExportImage PUBLIC ${JPEG_LIBRARIES})
+    TARGET_INCLUDE_DIRECTORIES(JFJochImageExport PRIVATE ${JPEG_INCLUDE_DIR})
+    TARGET_LINK_LIBRARIES(JFJochImageExport PUBLIC ${JPEG_LIBRARIES})
     MESSAGE(STATUS "JPEG headers present and library included")
 ELSE()
     MESSAGE(FATAL_ERROR "JPEG header jpeglib.h not found")

fpga/hls/CMakeLists.txt

@@ -1,4 +1,4 @@
-ADD_LIBRARY( HLSSimulation STATIC
+ADD_LIBRARY( JFJochHLSSimulation STATIC
         jf_conversion.cpp
         data_collection_fsm.cpp
         timer.cpp
@@ -30,9 +30,9 @@ ADD_LIBRARY( HLSSimulation STATIC
         spot_finder_connectivity.cpp
         spot_finder_mask.cpp)
 
-TARGET_INCLUDE_DIRECTORIES(HLSSimulation PUBLIC ../include)
-TARGET_LINK_LIBRARIES(HLSSimulation CommonFunctions)
-TARGET_COMPILE_DEFINITIONS(HLSSimulation PUBLIC -DJFJOCH_HLS_NOSYNTH)
+TARGET_INCLUDE_DIRECTORIES(JFJochHLSSimulation PUBLIC ../include)
+TARGET_LINK_LIBRARIES(JFJochHLSSimulation JFJochCommon)
+TARGET_COMPILE_DEFINITIONS(JFJochHLSSimulation PUBLIC -DJFJOCH_HLS_NOSYNTH)
 
 ADD_EXECUTABLE(frame_summation_tb frame_summation_tb.cpp)
 TARGET_LINK_LIBRARIES(frame_summation_tb HLSSimulation)

fpga/host_library/CMakeLists.txt

@@ -1,6 +1,6 @@
-ADD_LIBRARY(JungfraujochDevice STATIC JungfraujochDevice.cpp JungfraujochDevice.h ../include/jfjoch_fpga.h)
+ADD_LIBRARY(JFJochDevice STATIC JungfraujochDevice.cpp JungfraujochDevice.h ../include/jfjoch_fpga.h)
 
-TARGET_LINK_LIBRARIES(JungfraujochDevice CommonFunctions)
+TARGET_LINK_LIBRARIES(JFJochDevice JFJochCommon)
 
 ADD_EXECUTABLE(jfjoch_pcie_status jfjoch_pcie_status.cpp)
 TARGET_LINK_LIBRARIES(jfjoch_pcie_status JungfraujochDevice )

frame_serialize/CMakeLists.txt

@@ -24,4 +24,4 @@ ADD_LIBRARY(ImagePusher STATIC
 
 FIND_PACKAGE(OpenSSL REQUIRED)
 
-TARGET_LINK_LIBRARIES(ImagePusher CBORStream2FrameSerialize CommonFunctions Compression OpenSSL::Crypto)
+TARGET_LINK_LIBRARIES(ImagePusher CBORStream2FrameSerialize JFJochCommon Compression OpenSSL::Crypto)

frame_serialize/JFJochMessages.h

@@ -89,6 +89,8 @@ struct DataMessage {
     float beam_center_y;
     float incident_energy;
 
+    std::optional<float> resolution_estimation;
+
     std::optional<int64_t> roi_sum;
 };
 

frontend_ui/src/components/DataProcessingPlot.tsx

@@ -15,7 +15,8 @@ export enum PlotType {
     ERROR_PIXELS,
     RECEIVER_DELAY,
     STRONG_PIXELS,
-    ROI_SUM
+    ROI_SUM,
+    RES_ESTIMATION
 }
 
 type MyProps = {
@@ -120,6 +121,13 @@ class DataProcessingPlot extends Component<MyProps, MyState> {
                         this.setState({connection_error: true});
                     });
                 break;
+            case PlotType.RES_ESTIMATION:
+                DefaultService.getPlotResolutionEstimateHistogram()
+                    .then(data => this.setState({plot: data, connection_error: false}))
+                    .catch(error => {
+                        this.setState({connection_error: true});
+                    });
+                break;
         }
     }
     componentDidMount() {

frontend_ui/src/components/DataProcessingPlots.tsx

@@ -44,13 +44,13 @@ class DataProcessingPlots extends Component<MyProps, MyState> {
                 this.setState({type: PlotType.SPOT_COUNT, xlabel: "Image number", ylabel: "Spot count" });
                 break;
             case "3":
-                this.setState({type: PlotType.RAD_INT, xlabel: "Q [A<sup>-1</sup>]", ylabel: "Photon count"});
+                this.setState({type: PlotType.RAD_INT, xlabel: "Q [&#8491;<sup>-1</sup>]", ylabel: "Photon count"});
                 break;
             case "4":
                 this.setState({type: PlotType.INDEXING_RATE_PER_FILE, xlabel: "File number", ylabel: "Indexing rate"});
                 break;
             case "5":
-                this.setState({type: PlotType.RAD_INT_PER_FILE, xlabel: "Q [A^-1]", ylabel: "Photon count"});
+                this.setState({type: PlotType.RAD_INT_PER_FILE, xlabel: "Q [&#8491;<sup>-1</sup>]", ylabel: "Photon count"});
                 break;
             case "6":
                 this.setState({type: PlotType.SATURATED_PIXELS, xlabel: "Image number", ylabel: "Count"});
@@ -70,6 +70,9 @@ class DataProcessingPlots extends Component<MyProps, MyState> {
             case "11":
                 this.setState({type: PlotType.ROI_SUM, xlabel: "Image number", ylabel: "Photon count"});
                 break;
+            case "12":
+                this.setState({type: PlotType.RES_ESTIMATION, xlabel: "Resolution [&#8491;]", ylabel: "Number of crystals"});
+                break;
         }
     };
 
@@ -94,6 +97,7 @@ class DataProcessingPlots extends Component<MyProps, MyState> {
                             <MenuItem value={2}>Spot count</MenuItem>
                             <MenuItem value={3}>Azimuthal integration profile</MenuItem>
                             <MenuItem value={11}>ROI area sum</MenuItem>
+                            <MenuItem value={12}>Crystal resolution histogram</MenuItem>
                             <MenuItem value={4}>Indexing rate (per file)</MenuItem>
                             <MenuItem value={5}>Azimuthal integration profile (per file)</MenuItem>
                             <MenuItem value={6}>Saturated pixels</MenuItem>

frontend_ui/src/openapi/services/DefaultService.ts

@@ -544,6 +544,19 @@ export class DefaultService {
         });
     }
 
+    /**
+     * Generate resolution estimate histogram
+     * Generate histogram of crystal resolutions from 1.0 to 5.0 A based on ML model
+     * @returns plot Everything OK
+     * @throws ApiError
+     */
+    public static getPlotResolutionEstimateHistogram(): CancelablePromise<plot> {
+        return __request(OpenAPI, {
+            method: 'GET',
+            url: '/plot/resolution_estimate_histogram',
+        });
+    }
+
     /**
      * Generate radial integration profile
      * Generate average radial integration profile

image_analysis/CMakeLists.txt

@@ -1,30 +1,30 @@
-ADD_LIBRARY(ImageAnalysis STATIC
+ADD_LIBRARY(JFJochImageAnalysis STATIC
         CrystalLattice.cpp CrystalLattice.h
         IndexerWrapper.cpp IndexerWrapper.h
         AzimuthalIntegrationMapping.cpp AzimuthalIntegrationMapping.h
         StrongPixelSet.cpp StrongPixelSet.h AzimuthalIntegrationProfile.cpp AzimuthalIntegrationProfile.h
         SpotFindingSettings.h)
 
-TARGET_LINK_LIBRARIES(ImageAnalysis CommonFunctions)
+TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochCommon)
 
-TARGET_INCLUDE_DIRECTORIES(ImageAnalysis PUBLIC fast-feedback-indexer/eigen)
+TARGET_INCLUDE_DIRECTORIES(JFJochImageAnalysis PUBLIC fast-feedback-indexer/eigen)
 
 IF (CMAKE_CUDA_COMPILER)
-    TARGET_SOURCES(ImageAnalysis PRIVATE
+    TARGET_SOURCES(JFJochImageAnalysis PRIVATE
             fast-feedback-indexer/indexer/src/indexer.cpp
             fast-feedback-indexer/indexer/src/indexer_gpu.cu
             fast-feedback-indexer/indexer/src/log.cpp)
 
-    TARGET_SOURCES(ImageAnalysis PUBLIC
+    TARGET_SOURCES(JFJochImageAnalysis PUBLIC
             fast-feedback-indexer/indexer/src/ffbidx/indexer.h
             fast-feedback-indexer/indexer/src/ffbidx/indexer_gpu.h
             fast-feedback-indexer/indexer/src/ffbidx/refine.h
             fast-feedback-indexer/indexer/src/ffbidx/log.h
             fast-feedback-indexer/indexer/src/ffbidx/exception.h)
 
-    TARGET_INCLUDE_DIRECTORIES(ImageAnalysis PUBLIC fast-feedback-indexer/indexer/src/)
+    TARGET_INCLUDE_DIRECTORIES(JFJochImageAnalysis PUBLIC fast-feedback-indexer/indexer/src/)
 
-    TARGET_LINK_LIBRARIES(ImageAnalysis ${CUDART_LIBRARY} ${CMAKE_DL_LIBS} rt)
+    TARGET_LINK_LIBRARIES(JFJochImageAnalysis ${CUDART_LIBRARY} ${CMAKE_DL_LIBS} rt)
 ELSE()
     MESSAGE(WARNING "CUDA is strongly recommended for image analysis." )
 ENDIF()

image_analysis/IndexerWrapper.h

@@ -23,9 +23,7 @@ struct IndexingResult {
 
 class IndexerWrapper {
 #ifdef JFJOCH_USE_CUDA
-    fast_feedback::config_runtime<float> crt{
-            .num_sample_points = 32768
-    };
+    fast_feedback::config_runtime<float> crt{};
     fast_feedback::config_persistent<float> cpers{
             .max_output_cells = 4,
             .max_input_cells = 1,

receiver/CMakeLists.txt

@@ -11,7 +11,7 @@ ADD_LIBRARY(JFJochReceiver STATIC
         PreviewCounter.cpp
         PreviewCounter.h)
 
-TARGET_LINK_LIBRARIES(JFJochReceiver ImagePusher ImageAnalysis JungfraujochAcqusitionDevice CommonFunctions HLSSimulation ExportImage)
+TARGET_LINK_LIBRARIES(JFJochReceiver ImagePusher JFJochImageAnalysis JFJochAcquisitionDevice JFJochCommon JFJochHLSSimulation JFJochImageExport JFJochResonet)
 
 ADD_EXECUTABLE(jfjoch_action_test jfjoch_action_test.cpp)
 TARGET_LINK_LIBRARIES(jfjoch_action_test JFJochReceiver)

receiver/JFJochReceiver.cpp

@@ -7,6 +7,7 @@
 #include "../image_analysis/IndexerWrapper.h"
 #include "../common/DiffractionGeometry.h"
 #include "ImageMetadata.h"
+#include "../resonet/NeuralNetResPredictor.h"
 
 inline std::string time_UTC(const std::chrono::time_point<std::chrono::system_clock> &input) {
     auto time_ms = std::chrono::duration_cast<std::chrono::milliseconds>(input.time_since_epoch()).count();
@@ -269,6 +270,19 @@ void JFJochReceiver::FrameTransformationThread() {
         Cancel(e);
     }
 
+    std::unique_ptr<NeuralNetResPredictor> neural_net_predictor;
+
+#ifdef JFJOCH_USE_TORCH
+    if (!experiment.GetNeuralNetModelPath().empty()) {
+        try {
+            neural_net_predictor = std::make_unique<NeuralNetResPredictor>(experiment);
+        } catch (const JFJochException &e) {
+            logger.Warning("Torch error {}", e.what());
+            neural_net_predictor = nullptr;
+        }
+    }
+#endif
+
     FrameTransformation transformation(experiment);
     std::vector<char> writer_buffer(experiment.GetMaxCompressedSize());
 
@@ -382,6 +396,14 @@ void JFJochReceiver::FrameTransformationThread() {
             if (experiment.GetROISummation())
                 message.roi_sum = transformation.CalculateROISum(experiment.GetROISummation().value());
 
+            if (neural_net_predictor) {
+                try {
+                    message.resolution_estimation = neural_net_predictor->Inference(transformation.GetImage());
+                } catch (const std::exception &e) {
+                    logger.ErrorException(e);
+                    message.resolution_estimation.reset();
+                }
+            }
             plots.Add(message, image_number % experiment.GetDataFileCount(), az_int_profile_image);
 
             message.image = transformation.GetCompressedImage();
@@ -603,8 +625,8 @@ JFJochReceiverStatus JFJochReceiver::GetStatus() const {
     };
 }
 
-std::string JFJochReceiver::GetTIFF(bool calibration) const {
-    if (calibration)
+std::string JFJochReceiver::GetTIFF(bool calibration_run) const {
+    if (calibration_run)
         return preview_image.GenerateTIFFDioptas();
     else
         return preview_image.GenerateTIFF();

receiver/JFJochReceiver.h

@@ -20,8 +20,6 @@
 #include "../common/ThreadSafeFIFO.h"
 #include "../frame_serialize/ZMQStream2PreviewPublisher.h"
 #include "../common/NUMAHWPolicy.h"
-#include "../common/StatusVector.h"
-#include "../common/Histogram.h"
 
 #include "../image_analysis/StrongPixelSet.h"
 #include "../image_analysis/AzimuthalIntegrationMapping.h"

receiver/JFJochReceiverPlots.cpp

@@ -6,7 +6,7 @@ JFJochReceiverPlots::JFJochReceiverPlots(const DiffractionExperiment &experiment
                                          const AzimuthalIntegrationMapping &mapping)
         : indexing_solution_per_file(experiment.GetDataFileCount()),
           default_binning(experiment.GetSpotFindingBin()),
-          rad_int_profile(mapping) {
+          rad_int_profile(mapping), resolution_estimation(50, 1.0, 5.0) {
     for (int i = 0; i < experiment.GetDataFileCount(); i++)
         rad_int_profile_per_file.emplace_back(mapping);
 }
@@ -26,6 +26,8 @@ void JFJochReceiverPlots::Add(const DataMessage &msg, uint64_t file_number, cons
 
     if (msg.roi_sum)
         roi_sum.AddElement(file_number, msg.roi_sum.value());
+    if (msg.resolution_estimation)
+        resolution_estimation.Add(msg.resolution_estimation.value());
 
     rad_int_profile += profile;
     rad_int_profile_per_file.at(file_number) += profile;
@@ -64,6 +66,8 @@ Plot JFJochReceiverPlots::GetPlots(const PlotRequest &request) {
             return  strong_pixels.GetPlot(nbins);
         case PlotType::ROISum:
             return roi_sum.GetPlot(nbins);
+        case PlotType::ResEstimation:
+            return resolution_estimation.GetPlot();
         default:
             // Do nothing
             break;

receiver/JFJochReceiverPlots.h

@@ -22,8 +22,9 @@ class JFJochReceiverPlots {
     StatusVector<uint64_t> receiver_delay;
     StatusVector<float> image_collection_efficiency;
     StatusVector<float> unit_cell[6];
-    Histogram<uint64_t> indexing_solution_per_file;
+    SetAverage<uint64_t> indexing_solution_per_file;
     StatusVector<int64_t> roi_sum;
+    FloatHistogram resolution_estimation;
 
     AzimuthalIntegrationProfile rad_int_profile;
     std::vector<AzimuthalIntegrationProfile> rad_int_profile_per_file;

receiver/jfjoch_action_test.cpp

@@ -19,6 +19,7 @@ void print_usage(Logger &logger) {
     logger.Info("   -p<num> data processing period");
     logger.Info("   -N<num> number of image processing threads");
     logger.Info("   -P<txt> NUMA Policy (none|n2g2|n8g4|n8g4_hbm), none is default");
+    logger.Info("   -D<path> use resonet deep learning model for resolution estimation - path to TorchScript");
 }
 
 int main(int argc, char **argv) {
@@ -35,6 +36,7 @@ int main(int argc, char **argv) {
     std::string numa_policy_name;
     bool raw_data = false;
     bool force_32bit = false;
+    std::string resonet_path;
 
     if (argc == 1) {
         print_usage(logger);
@@ -71,6 +73,9 @@ int main(int argc, char **argv) {
             case 'I':
                 force_32bit = true;
                 break;
+            case 'D':
+                resonet_path = std::string(optarg);
+                break;
             default: /* '?' */
                 print_usage(logger);
                 exit(EXIT_FAILURE);
@@ -94,6 +99,8 @@ int main(int argc, char **argv) {
     x.Compression(CompressionAlgorithm::BSHUF_LZ4).DataStreams(nstreams);
     if (force_32bit)
         x.FPGAOutputMode(FPGAPixelOutput::Int32);
+    if (!resonet_path.empty())
+        x.NeuralNetModelPath(resonet_path);
 
     logger.Info("Data streams {} Total modules {} Total images {} Threads {}", nstreams, nmodules, nimages, nthreads);
 

resonet/CMakeLists.txt

@@ -0,0 +1,10 @@
+ADD_LIBRARY(JFJochResonet STATIC NeuralNetResPredictor.cpp NeuralNetResPredictor.h)
+TARGET_LINK_LIBRARIES(JFJochResonet JFJochCommon)
+
+IF (${TORCH_FOUND})
+    TARGET_COMPILE_DEFINITIONS(JFJochResonet PUBLIC -DJFJOCH_USE_TORCH)
+    TARGET_LINK_LIBRARIES(JFJochResonet ${TORCH_LIBRARIES})
+    TARGET_INCLUDE_DIRECTORIES(JFJochResonet PUBLIC ${TORCH_INCLUDE_DIRS})
+    ADD_EXECUTABLE(resonet_test resonet_test.cpp)
+    TARGET_LINK_LIBRARIES(resonet_test JFJochResonet JFJochWriter)
+ENDIF()

resonet/NeuralNetResPredictor.cpp

@@ -0,0 +1,93 @@
+// Copyright (2019-2024) Paul Scherrer Institute
+
+#include "NeuralNetResPredictor.h"
+
+#include <cmath>
+#include "../common/JFJochException.h"
+
+NeuralNetResPredictor::NeuralNetResPredictor(const DiffractionExperiment &in_experiment)
+        : experiment(in_experiment), model_input(512*512)
+#ifdef JFJOCH_USE_TORCH
+          , device(torch::kCUDA)
+#endif
+        {
+    float max_direction = std::max(in_experiment.GetXPixelsNum(), in_experiment.GetYPixelsNum()) / 2.0;
+    pool_factor = std::lround(max_direction / 512.0f);
+    if (pool_factor <= 0)
+        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Detector size is too small");
+
+    if (pool_factor > 8)
+        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Detector size is too large");
+#ifdef JFJOCH_USE_TORCH
+    module = torch::jit::load(experiment.GetNeuralNetModelPath());
+    module.to(device);
+#endif
+}
+
+template<class T>
+void NeuralNetResPredictor::PrepareInternal(const T* image) {
+    size_t xpixel = experiment.GetXPixelsNum();
+    size_t ypixel = experiment.GetYPixelsNum();
+    size_t start_x = std::lround(experiment.GetBeamX_pxl());
+    size_t start_y = std::lround(experiment.GetBeamY_pxl());
+
+    for (size_t y = 0; y < 512; y++) {
+        size_t y0 = y * pool_factor + start_y;
+        size_t min_yp = std::min(y0 + pool_factor, ypixel);
+        for (size_t x = 0; x < 512; x++) {
+            float val = 0.0;
+            size_t x0 = x * pool_factor + start_x;
+            size_t min_xp = std::min(x0 + pool_factor, ypixel);
+
+            for (size_t yp = y0; yp < min_yp; yp++) {
+                for (size_t xp = x0; xp < min_xp; xp++) {
+                    int16_t pxl = image[yp * xpixel + xp];
+                    if (val < pxl)
+                        val = pxl;
+                }
+            }
+            float max_pool = floorf(sqrtf(val));
+            model_input[512 * y + x] = max_pool;
+        }
+    }
+}
+
+void NeuralNetResPredictor::Prepare(const int16_t *image) {
+    PrepareInternal(image);
+}
+
+void NeuralNetResPredictor::Prepare(const int32_t *image) {
+    PrepareInternal(image);
+}
+
+size_t NeuralNetResPredictor::GetMaxPoolFactor() const {
+    return pool_factor;
+}
+
+float NeuralNetResPredictor::Inference(const void *image) {
+    if (experiment.GetPixelDepth() == 2)
+        Prepare((int16_t *) image);
+    else
+        Prepare((int32_t *) image);
+
+#ifdef JFJOCH_USE_TORCH
+    auto options = torch::TensorOptions().dtype(at::kFloat);
+    auto model_input_tensor = torch::from_blob(model_input.data(), {1,1,512,512}, options).to(device);
+
+    std::vector<torch::jit::IValue> pixels;
+    pixels.emplace_back(std::move(model_input_tensor));
+
+    auto output = module.forward(pixels).toTensor();
+    auto tensor_output = output[0].item<float>();
+    float two_theta = atanf(((2.0f * experiment.GetPixelSize_mm() / experiment.GetDetectorDistance_mm()) * tensor_output));
+    float stheta = sinf(two_theta * 0.5f);
+    float resolution = experiment.GetWavelength_A() / (2.0f * stheta);
+#else
+    float resolution = 50.0;
+#endif
+    return resolution;
+}
+
+const std::vector<float> &NeuralNetResPredictor::GetModelInput() const {
+    return model_input;
+}

resonet/NeuralNetResPredictor.h

@@ -0,0 +1,37 @@
+// Copyright (2019-2024) Paul Scherrer Institute
+
+#ifndef JUNGFRAUJOCH_NEURALNETRESPREDICTOR_H
+#define JUNGFRAUJOCH_NEURALNETRESPREDICTOR_H
+
+#ifdef JFJOCH_USE_TORCH
+#include <torch/script.h>
+#endif
+
+#include "../common/DiffractionExperiment.h"
+
+// Based on model described in:
+// Mendez, D., Holton, J. M., Lyubimov, A. Y., Hollatz, S., Mathews, I. I., Cichosz, A., Martirosyan, V.,
+// Zeng, T., Stofer, R., Liu, R., Song, J., McPhillips, S., Soltis, M. & Cohen, A. E. (2024).
+// Acta Cryst. D80, 26-43.
+
+class NeuralNetResPredictor {
+    const DiffractionExperiment& experiment;
+    size_t pool_factor;
+    std::vector<float> model_input;
+#ifdef JFJOCH_USE_TORCH
+    torch::Device device;
+    torch::jit::script::Module module;
+#endif
+    template<class T>
+    void PrepareInternal(const T* image);
+public:
+    explicit NeuralNetResPredictor(const DiffractionExperiment& experiment);
+    void Prepare(const int16_t* image);
+    void Prepare(const int32_t* image);
+    float Inference(const void* image);
+    size_t GetMaxPoolFactor() const;
+    const std::vector<float> &GetModelInput() const;
+};
+
+
+#endif //JUNGFRAUJOCH_NEURALNETRESPREDICTOR_H

resonet/resonet_test.cpp

@@ -0,0 +1,42 @@
+// Copyright (2019-2024) Paul Scherrer Institute
+
+#include <iostream>
+
+#include "NeuralNetResPredictor.h"
+#include "../writer/HDF5Objects.h"
+
+int main(int argc, char **argv) {
+    DiffractionExperiment experiment(DetectorGeometry(8, 2, 8, 36));
+    experiment.DetectorDistance_mm(75).PhotonEnergy_keV(12.4).BeamY_pxl(1136).BeamX_pxl(1090);
+    experiment.NeuralNetModelPath(argv[1]);
+    if (argc == 0)
+        exit(EXIT_FAILURE);
+
+    NeuralNetResPredictor predictor(experiment);
+
+    HDF5ReadOnlyFile data("../../tests/test_data/compression_benchmark.h5");
+    HDF5DataSet dataset(data, "/entry/data/data");
+    HDF5DataSpace file_space(dataset);
+
+    std::vector<int16_t> image_conv (file_space.GetDimensions()[1] * file_space.GetDimensions()[2]);
+
+    std::vector<hsize_t> start = {4,0,0};
+    std::vector<hsize_t> file_size = {1, file_space.GetDimensions()[1], file_space.GetDimensions()[2]};
+    dataset.ReadVector(image_conv, start, file_size);
+
+    std::cout << "Max pooling " << predictor.GetMaxPoolFactor() << std::endl;
+
+    float x = 0;
+    size_t niterations = 25;
+
+    auto start_time = std::chrono::system_clock::now();
+    for (int i = 0; i < niterations; i++)
+        x += predictor.Inference(image_conv.data());
+    auto end_time = std::chrono::system_clock::now();
+    auto elapsed = std::chrono::duration_cast<std::chrono::microseconds>(end_time - start_time);
+    int64_t frequency_Hz = (niterations * 1e6) / (double) (elapsed.count());
+
+    std::cout << ((double) elapsed.count()) / (1e3 * niterations) << " ms" << std::endl;
+    std::cout << frequency_Hz << " Hz" << std::endl;
+    std::cout << "Resolution " << predictor.Inference(image_conv.data()) << std::endl;
+}

tests/CMakeLists.txt

@@ -30,9 +30,11 @@ ADD_EXECUTABLE(CatchTest
         TIFFTest.cpp
         JFJochReceiverProcessingTest.cpp
         JPEGTest.cpp
+        NeuralNetResPredictorTest.cpp
+        HistogramTest.cpp
 )
 
-target_link_libraries(CatchTest JFJochBroker JFJochReceiver JFJochWriter ImageAnalysis CommonFunctions HLSSimulation ExportImage)
+target_link_libraries(CatchTest JFJochBroker JFJochReceiver JFJochWriter JFJochImageAnalysis JFJochCommon JFJochHLSSimulation JFJochImageExport JFJochResonet)
 target_include_directories(CatchTest PRIVATE .)
 
 INSTALL(TARGETS CatchTest RUNTIME)

tests/HistogramTest.cpp

@@ -0,0 +1,53 @@
+// Copyright (2019-2024) Paul Scherrer Institute
+
+#include <catch2/catch.hpp>
+
+#include "../common/Histogram.h"
+
+TEST_CASE("SetAverage") {
+    SetAverage<float> h(100);
+    h.Add(80, 30.0);
+    h.Add(50, 20.0);
+    h.Add(50, 30.0);
+    h.Add(0, -2.0);
+    h.Add(100, 50.0);
+
+    auto p = h.GetPlot();
+    REQUIRE(p.x.size() == 100);
+    REQUIRE(p.y.size() == 100);
+
+    CHECK(p.x[0] == 0.0);
+    CHECK(p.x[80] == 80.0);
+    CHECK(p.x[30] == 30.0);
+    CHECK(p.x[99] == 99.0);
+
+    CHECK(p.y[0] == -2.0);
+    CHECK(p.y[50] == 25.0);
+    CHECK(p.y[80] == 30.0);
+    CHECK(p.y[99] == 0.0);
+
+}
+
+TEST_CASE("FloatHistogram") {
+    FloatHistogram h(100, 100.0, 200.0);
+    h.Add(150.5);
+    h.Add(100.2);
+    h.Add(100.9999);
+    h.Add(100.0);
+    h.Add(100.00001);
+    h.Add(101.0);
+    h.Add(200.0);
+
+    auto p = h.GetPlot();
+    REQUIRE(p.x.size() == 100);
+    REQUIRE(p.y.size() == 100);
+
+    CHECK(p.x[0] == 100.5);
+    CHECK(p.x[34] == 134.5);
+    CHECK(p.x[99] == 199.5);
+
+    CHECK(p.y[0] == 4.0);
+    CHECK(p.y[1] == 1.0);
+    CHECK(p.y[50] == 1.0);
+    CHECK(p.y[99] == 0.0);
+}

tests/JFJochReceiverProcessingTest.cpp

@@ -173,3 +173,65 @@ TEST_CASE("GenerateResolutionMap") {
     WriteTIFFToFile("ResolutionMap.tiff", spot_finder_resolution_map_int_conv.data(), experiment.GetXPixelsNum(),
                     experiment.GetYPixelsNum(), 4);
 }
+
+TEST_CASE("JFJochIntegrationTest_ZMQ_lysozyme_resolution", "[JFJochReceiver]") {
+    Logger logger("JFJochIntegrationTest_ZMQ_lysozyme_resolution");
+
+    RegisterHDF5Filter();
+
+    const uint16_t nthreads = 4;
+
+    DiffractionExperiment experiment(DetectorGeometry(8,2,8,36));
+    experiment.ImagesPerTrigger(5).NumTriggers(1).DataFileCount(1).UseInternalPacketGenerator(true)
+            .FilePrefix("lyso_test_resolution").ConversionOnFPGA(false)
+            .DetectorDistance_mm(75).BeamY_pxl(1136).BeamX_pxl(1090).PhotonEnergy_keV(12.4)
+            .SetUnitCell(UnitCell{.a = 36.9, .b = 78.95, .c = 78.95, .alpha =90, .beta = 90, .gamma = 90})
+            .NeuralNetModelPath("../../resonet/traced_resnet_model.pt");
+
+    // Load example image
+    HDF5ReadOnlyFile data("../../tests/test_data/compression_benchmark.h5");
+    HDF5DataSet dataset(data, "/entry/data/data");
+    HDF5DataSpace file_space(dataset);
+
+    REQUIRE(file_space.GetDimensions()[2] == experiment.GetXPixelsNum());
+    REQUIRE(file_space.GetDimensions()[1] == experiment.GetYPixelsNum());
+    std::vector<int16_t> image_conv (file_space.GetDimensions()[1] * file_space.GetDimensions()[2]);
+
+    std::vector<hsize_t> start = {4,0,0};
+    std::vector<hsize_t> file_size = {1, file_space.GetDimensions()[1], file_space.GetDimensions()[2]};
+    dataset.ReadVector(image_conv, start, file_size);
+
+    std::vector<int16_t> image_raw_geom(experiment.GetModulesNum() * RAW_MODULE_SIZE);
+    ConvertedToRawGeometry(experiment, image_raw_geom.data(), image_conv.data());
+    logger.Info("Loaded image");
+
+    // Setup acquisition device
+    AcquisitionDeviceGroup aq_devices;
+    std::unique_ptr<HLSSimulatedDevice> test = std::make_unique<HLSSimulatedDevice>(0, 64);
+    test->SetInternalGeneratorFrame(image_raw_geom);
+    aq_devices.Add(std::move(test));
+
+    JFCalibration calib(experiment);
+    HDF5ImagePusher pusher;
+
+    ZMQContext context;
+    context.NumThreads(4);
+
+    JFJochReceiverService service(aq_devices, logger, pusher);
+    service.NumThreads(nthreads);
+
+    SpotFindingSettings settings = DiffractionExperiment::DefaultDataProcessingSettings();
+    settings.signal_to_noise_threshold = 2.5;
+    settings.photon_count_threshold = 5;
+    settings.min_pix_per_spot = 1;
+    settings.max_pix_per_spot = 200;
+    service.SetSpotFindingSettings(settings);
+
+    service.Start(experiment, nullptr);
+    auto receiver_out = service.Stop();
+
+    CHECK(receiver_out.efficiency == 1.0);
+    CHECK(receiver_out.status.indexing_rate == 1.0);
+    CHECK(receiver_out.status.images_sent == experiment.GetImageNum());
+    CHECK(!receiver_out.status.cancelled);
+}

tests/NeuralNetResPredictorTest.cpp

@@ -0,0 +1,53 @@
+// Copyright (2019-2024) Paul Scherrer Institute
+
+#include <catch2/catch.hpp>
+
+#include "../writer/HDF5Objects.h"
+#include "../resonet/NeuralNetResPredictor.h"
+
+TEST_CASE("NeuralNetResPredictor_Prepare", "[LinearAlgebra][Coord]") {
+    DiffractionExperiment experiment(DetectorGeometry(8, 2, 8, 36));
+    experiment.DetectorDistance_mm(75).PhotonEnergy_keV(12.4).BeamX_pxl(1000).BeamY_pxl(1000);
+    experiment.NeuralNetModelPath("../../resonet/traced_resnet_model.pt");
+
+    std::vector<int16_t> v(experiment.GetPixelsNum(),0);
+    v[1000 * experiment.GetXPixelsNum() + 1000] = 100;
+    v[1000 * experiment.GetXPixelsNum() + 1001] = 20;
+    v[1001 * experiment.GetXPixelsNum() + 1000] = 30;
+    v[1001 * experiment.GetXPixelsNum() + 1001] = INT16_MIN;
+
+    v[1050 * experiment.GetXPixelsNum() + 1050] = 52;
+    v[2000 * experiment.GetXPixelsNum() + 1500] = 160;
+    NeuralNetResPredictor predictor(experiment);
+
+    REQUIRE(predictor.GetMaxPoolFactor() == 2);
+
+    predictor.Prepare(v.data());
+    auto nn_input = predictor.GetModelInput();
+    REQUIRE(nn_input[0] == 10);
+    REQUIRE(nn_input[25 * 512 + 25] == 7);
+    REQUIRE(nn_input[500 * 512 + 250] == 12);
+}
+
+TEST_CASE("NeuralNetResPredictor_Inference", "[LinearAlgebra][Coord]") {
+    DiffractionExperiment experiment(DetectorGeometry(8, 2, 8, 36));
+    experiment.DetectorDistance_mm(75).PhotonEnergy_keV(12.4).BeamY_pxl(1136).BeamX_pxl(1090);
+    experiment.NeuralNetModelPath("../../resonet/traced_resnet_model.pt");
+
+    NeuralNetResPredictor predictor(experiment);
+
+    HDF5ReadOnlyFile data("../../tests/test_data/compression_benchmark.h5");
+    HDF5DataSet dataset(data, "/entry/data/data");
+    HDF5DataSpace file_space(dataset);
+
+    std::vector<int16_t> image_conv (file_space.GetDimensions()[1] * file_space.GetDimensions()[2]);
+
+    std::vector<hsize_t> start = {4,0,0};
+    std::vector<hsize_t> file_size = {1, file_space.GetDimensions()[1], file_space.GetDimensions()[2]};
+    dataset.ReadVector(image_conv, start, file_size);
+
+    auto res = predictor.Inference(image_conv.data());
+    std::cout << res << std::endl;
+    REQUIRE(res < 1.5);
+    REQUIRE(res > 1.4);
+}

tools/CMakeLists.txt

@@ -1,23 +1,23 @@
 
 ADD_EXECUTABLE(jfjoch_udp_simulator jfjoch_udp_simulator.cpp UDPSimulator.cpp UDPSimulator.h)
-TARGET_LINK_LIBRARIES(jfjoch_udp_simulator CommonFunctions)
+TARGET_LINK_LIBRARIES(jfjoch_udp_simulator JFJochCommon)
 
 add_executable(CompressionBenchmark CompressionBenchmark.cpp)
-target_link_libraries(CompressionBenchmark JFJochWriter JFJochReceiver CommonFunctions)
+target_link_libraries(CompressionBenchmark JFJochWriter JFJochReceiver JFJochCommon)
 
 add_executable(HDF5DatasetWriteTest HDF5DatasetWriteTest.cpp)
-target_link_libraries(HDF5DatasetWriteTest JFJochWriter JFJochReceiver CommonFunctions)
+target_link_libraries(HDF5DatasetWriteTest JFJochWriter JFJochReceiver JFJochCommon)
 
 ADD_EXECUTABLE(jfjoch_preview_test jfjoch_preview_test.cpp)
-TARGET_LINK_LIBRARIES(jfjoch_preview_test JFJochReceiver ImagePusher CommonFunctions)
+TARGET_LINK_LIBRARIES(jfjoch_preview_test JFJochReceiver ImagePusher JFJochCommon)
 
 ADD_EXECUTABLE(jfjoch_writer_test jfjoch_writer_test.cpp)
-TARGET_LINK_LIBRARIES(jfjoch_writer_test JFJochReceiver JFJochWriter  ImagePusher CommonFunctions)
+TARGET_LINK_LIBRARIES(jfjoch_writer_test JFJochReceiver JFJochWriter  ImagePusher JFJochCommon)
 
 ADD_EXECUTABLE(gain_file_statistics gain_file_statistics.cpp)
-TARGET_LINK_LIBRARIES(gain_file_statistics CommonFunctions)
+TARGET_LINK_LIBRARIES(gain_file_statistics JFJochCommon)
 
 ADD_EXECUTABLE(SpotFindingPerformanceTest SpotFindingPerformanceTest.cpp)
-TARGET_LINK_LIBRARIES(SpotFindingPerformanceTest ImageAnalysis CommonFunctions)
+TARGET_LINK_LIBRARIES(SpotFindingPerformanceTest JFJochImageAnalysis JFJochCommon)
 
 INSTALL(TARGETS jfjoch_udp_simulator CompressionBenchmark HDF5DatasetWriteTest jfjoch_writer_test RUNTIME)

writer/CMakeLists.txt

@@ -1,8 +1,8 @@
 FIND_PACKAGE(HDF5 1.10 REQUIRED)
 
-ADD_LIBRARY(HDF5Wrappers STATIC HDF5Objects.cpp HDF5Objects.h HDF5DataFile.h  HDF5DataFile.cpp  ../compression/bitshuffle/bshuf_h5filter.c)
-TARGET_LINK_LIBRARIES(HDF5Wrappers Compression ${HDF5_LIBRARIES})
-TARGET_INCLUDE_DIRECTORIES(HDF5Wrappers PUBLIC ${HDF5_INCLUDE_DIRS})
+ADD_LIBRARY(JFJochHDF5Wrappers STATIC HDF5Objects.cpp HDF5Objects.h HDF5DataFile.h  HDF5DataFile.cpp  ../compression/bitshuffle/bshuf_h5filter.c)
+TARGET_LINK_LIBRARIES(JFJochHDF5Wrappers Compression ${HDF5_LIBRARIES})
+TARGET_INCLUDE_DIRECTORIES(JFJochHDF5Wrappers PUBLIC ${HDF5_INCLUDE_DIRS})
 
 ADD_LIBRARY(JFJochWriter STATIC
         HDF5DataFile.h  HDF5DataFile.cpp
@@ -13,10 +13,10 @@ ADD_LIBRARY(JFJochWriter STATIC
         HDF5ImagePusher.cpp
         HDF5ImagePusher.h)
 
-TARGET_LINK_LIBRARIES(JFJochWriter HDF5Wrappers CommonFunctions CBORStream2FrameSerialize)
+TARGET_LINK_LIBRARIES(JFJochWriter JFJochHDF5Wrappers JFJochCommon CBORStream2FrameSerialize)
 
 ADD_EXECUTABLE(HDF5Sum HDF5Sum.cpp)
-TARGET_LINK_LIBRARIES(HDF5Sum HDF5Wrappers)
+TARGET_LINK_LIBRARIES(HDF5Sum JFJochHDF5Wrappers)
 
 AUX_SOURCE_DIRECTORY(gen/model MODEL_SOURCES)
 ADD_LIBRARY(WriterAPI STATIC ${MODEL_SOURCES} gen/api/DefaultApi.cpp gen/api/DefaultApi.h JFJochWriterHttp.h JFJochWriterHttp.cpp)

writer/HDF5DataFile.cpp

@@ -33,12 +33,16 @@ HDF5DataFile::~HDF5DataFile() {
     if (data_file) {
         HDF5Group group_exp(*data_file, "/entry/jungfrau");
         group_exp.NXClass("NXcollection");
+
         group_exp.SaveVector("info", jf_info);
         group_exp.SaveVector("timestamp", timestamp);
         group_exp.SaveVector("storage_cell", storage_cell);
         group_exp.SaveVector("exptime", exptime);
-        group_exp.SaveVector("receiver_aq_dev_delay", receiver_aq_dev_delay);
-        data_file->SaveVector("/entry/result/strong_pixel_count", strong_pixel_count);
+        group_exp.SaveVector("receiverAqDevDelay", receiver_aq_dev_delay);
+
+        data_file->SaveVector("/entry/result/strongPixelCount", strong_pixel_count);
+        if (!resolution_estimation.empty())
+            data_file->SaveVector("/entry/result/resolutionEstimation", resolution_estimation);
     }
 }
 
@@ -58,21 +62,18 @@ void HDF5DataFile::CreateFile() {
     data_set_indexed_lattice = std::make_unique<HDF5DataSet>(*data_file, "/entry/result/latticeIndexed");
     if (!rad_int_bin_to_q.empty())
         data_set_az_int = std::make_unique<HDF5DataSet>(*data_file, "/entry/result/azimIntegration");
-
-    jf_info.resize(1);
-    receiver_aq_dev_delay.resize(1);
-    timestamp.resize(1);
-    storage_cell.resize(1);
-    exptime.resize(1);
-    strong_pixel_count.resize(1);
 }
 
 void HDF5DataFile::Write(const DataMessage &msg, uint64_t image_number) {
     std::lock_guard<std::mutex> lock(hdf5_mutex);
-    if (!data_file)
-        CreateFile();
+    bool new_file = false;
 
-    if (image_number > max_image_number) {
+    if (!data_file) {
+        CreateFile();
+        new_file = true;
+    }
+
+    if (new_file || (static_cast<int64_t>(image_number) > max_image_number)) {
         max_image_number = image_number;
         data_set->SetExtent({max_image_number+1, ypixel, xpixel});
         data_set_spot_x->SetExtent({max_image_number+1, max_spots});
@@ -93,6 +94,9 @@ void HDF5DataFile::Write(const DataMessage &msg, uint64_t image_number) {
         timestamp.resize(max_image_number + 1);
         exptime.resize(max_image_number + 1);
         storage_cell.resize(max_image_number + 1);
+
+        if (msg.resolution_estimation)
+            resolution_estimation.resize(max_image_number + 1);
     }
 
     nimages++;
@@ -128,6 +132,9 @@ void HDF5DataFile::Write(const DataMessage &msg, uint64_t image_number) {
     storage_cell[image_number] = msg.storage_cell;
     exptime[image_number] = msg.exptime;
 
+    if (msg.resolution_estimation)
+        resolution_estimation[image_number] = msg.resolution_estimation.value();
+
     if (!msg.az_int_profile.empty() && (msg.az_int_profile.size() == rad_int_bin_to_q.size()))
         data_set_az_int->WriteVec(msg.az_int_profile, {image_number, 0}, {1, rad_int_bin_to_q.size()});
 

writer/HDF5DataFile.h

@@ -13,7 +13,7 @@
 
 struct HDF5DataFileStatistics {
     std::string filename;
-    uint64_t max_image_number;
+    size_t max_image_number;
     uint64_t total_images;
 };
 
@@ -49,6 +49,7 @@ class HDF5DataFile {
     std::vector<uint32_t> strong_pixel_count;
     std::vector<int64_t> receiver_aq_dev_delay;
     std::vector<float> rad_int_bin_to_q;
+    std::vector<float> resolution_estimation;
 
     const size_t max_spots;