Parallelize phantom creation for improved performance

CMakeLists.txt

@@ -28,6 +28,8 @@ else()
     message("CUDA not found. Compiling without GPU support.")
 endif()
 
+find_package(OpenMP)
+
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 add_compile_definitions(MINI_CASE_SENSITIVE)
@@ -46,13 +48,9 @@ add_subdirectory(./src/phantom)
 add_subdirectory(./src/config)
 add_subdirectory(./src/sim)
 
-# Add the executable
+ # Add the sources files
 set(SOURCES ./src/spinwalk.cu ${SUBCOMMAND_PHANTOM} ${SUBCOMMAND_SIM} ${SUBCOMMAND_DWI} ${SUBCOMMAND_CONFIG})
-
-if(CMAKE_CUDA_COMPILER)
-    add_executable(${project} ${SOURCES})
-    target_link_libraries(${project} ${CUDA_LIBRARIES} ${Boost_LIBRARIES} ${HDF5_CXX_LIBRARIES} TBB::tbb)
-else()    
+if(NOT  CMAKE_CUDA_COMPILER) 
     set(RENAMED_SOURCES)
     foreach(OLD_FILE ${SOURCES})
         # Get the file name without extension
@@ -64,12 +62,28 @@ else()
         file(COPY_FILE ${OLD_FILE} ${NEW_FILE})
         list(APPEND RENAMED_SOURCES ${NEW_FILE})
     endforeach()
-    add_executable(${project} ${RENAMED_SOURCES})
-    target_link_libraries(${project} ${Boost_LIBRARIES} ${HDF5_CXX_LIBRARIES} TBB::tbb)
+    set(SOURCES ${RENAMED_SOURCES})
+endif()
+
+# Add the executable
+add_executable(${project} ${SOURCES})
+
+# Add the libraries
+if(CMAKE_CUDA_COMPILER)
+    target_link_libraries(${project} PRIVATE ${CUDA_LIBRARIES})
+endif()
+if(OpenMP_CXX_FOUND)
+    message(STATUS "Found OpenMP, adding to target link libraries.")
+    target_link_libraries(${project} PRIVATE OpenMP::OpenMP_CXX)
+else()
+    message(STATUS "OpenMP not found, skipping.")
 endif()
+target_link_libraries(${project} PRIVATE ${Boost_LIBRARIES} ${HDF5_CXX_LIBRARIES} TBB::tbb)
 
+# Set the output name
 set_target_properties(${project} PROPERTIES OUTPUT_NAME "spinwalk")
 
+# Install the executable
 if (UNIX)
     install(TARGETS ${project} DESTINATION bin) 
 endif()

src/definitions.h

@@ -2,8 +2,8 @@
 #define DEFINITIONS_H
 
 #define VERSION_MAJOR 1
-#define VERSION_MINOR 17
-#define VERSION_PATCH 2
+#define VERSION_MINOR 18
+#define VERSION_PATCH 0
 
 // Helper macros to stringify values
 #define STRINGIFY_HELPER(x) #x

src/phantom/phantom_base.cpp

@@ -23,15 +23,15 @@ phantom_base::phantom_base()
     m_fov = 0;
     m_resolution = 0;
     m_seed = std::random_device{}();
-    set_blood_parameters(0.273e-6 * 0.4, 0, 10.0);
+    set_parameters(0.273e-6 * 0.4, 0, 10.0);
     set_filename(); 
 }
 
 phantom_base::phantom_base(float fov_um, size_t resolution, float dChi, float Y, float BVF, int32_t seed, std::string filename)
 :phantom_base()
 {
     set_space(fov_um, resolution);
-    set_blood_parameters(dChi, Y, BVF);
+    set_parameters(dChi, Y, BVF);
     set_filename(filename);
     if(seed >= 0)
         m_seed = seed;
@@ -47,11 +47,11 @@ void phantom_base::set_space(float fov_um, size_t resolution)
     this->m_resolution = resolution;
 }
 
-void phantom_base::set_blood_parameters(float dChi, float Y, float BVF)
+void phantom_base::set_parameters(float dChi, float Y, float volume_fraction)
 {
     this->m_dChi = dChi;
     this->m_Y = Y;
-    this->m_BVF = BVF;
+    this->m_volume_fraction = volume_fraction;
     m_calc_fieldmap = Y >= 0;
 }
 
@@ -92,7 +92,7 @@ bool phantom_base::save() const
     // blood volume fraction
     std::vector<size_t> dims_1(1, 1);
     HighFive::DataSet dataset_BVF = file.createDataSet<float>("bvf", HighFive::DataSpace(dims_1));
-    dataset_BVF.write_raw(&m_BVF);
+    dataset_BVF.write_raw(&m_volume_fraction);
 
     return true;
 }

src/phantom/phantom_base.h

@@ -25,16 +25,17 @@ class phantom_base
 {
     public:
         phantom_base();
-        phantom_base(float fov_um, size_t resolution, float dChi, float Y, float BVF, int32_t seed, std::string filename);
+        phantom_base(float fov_um, size_t resolution, float dChi, float Y, float volume_fraction, int32_t seed, std::string filename);
         virtual ~phantom_base();
         void set_space(float fov_um, size_t resolution);
-        void set_blood_parameters(float dChi, float Y, float BVF = 10.0); 
+        void set_parameters(float dChi, float Y, float volume_fraction = 10.0); 
         void set_filename(std::string filename = "shape.h5");      
         virtual bool run(){return true;}; 
         virtual bool save() const;
         virtual bool create_grid();
         virtual bool generate_shapes() = 0;
         virtual bool generate_mask_fieldmap() = 0;
+        float get_actual_volume_fraction() const {return m_volume_fraction;}
 
         friend std::ostream& operator<<(std::ostream& os, const phantom_base& obj);
 
@@ -44,7 +45,7 @@ class phantom_base
         std::vector<int8_t> m_mask;
         size_t m_resolution;
         float m_fov, m_dChi, m_Y;
-        float m_BVF; // blood volume fraction
+        float m_volume_fraction; // volume fraction
         std::string m_filename;
         float B0[3] = {0.f, 0.f, 1.f};
         bool m_calc_fieldmap;

src/phantom/phantom_cylinder.cpp

@@ -20,14 +20,51 @@
 
 namespace phantom
 {
+
+bool check_cylinder_overlap(const std::vector<std::vector<float>>& m_cylinder_points, const std::vector<float>& m_cylinder_radii, const float* cyl_pnt, float& radius, bool is_random_radius) {
+    int c = 0;
+    bool stop = false;
+
+    #pragma omp parallel for shared(stop, radius)
+    for (c = 0; c < m_cylinder_radii.size(); c++) {
+        if (stop) continue; // Skip computation if condition already met.
+
+        float p2p1[3];
+        subtract(cyl_pnt, m_cylinder_points[c].data(), p2p1);
+        float distance = sqrtf(p2p1[0]*p2p1[0] + p2p1[1]*p2p1[1]);;
+        // Check the conditions
+        if (distance <= m_cylinder_radii[c] || distance <= radius) {
+            #pragma omp critical
+            {
+                stop = true; // Signal all threads to stop further checks.
+            }
+        } else if (distance < m_cylinder_radii[c] + radius) {
+            if (!is_random_radius) {
+                #pragma omp critical
+                {
+                    stop = true; // Signal all threads to stop further checks.
+                }
+            } else {
+                #pragma omp critical
+                {
+                    if (!stop) 
+                        radius = distance - m_cylinder_radii[c]; // Adjust radius.
+                }
+            }
+        }
+    }
+    return stop; // Return whether an overlap condition was found.
+}
+
+
 cylinder::cylinder()
 {
     m_radius = 0;
     m_orientation = 0;
 }
 
-cylinder::cylinder(float fov_um, size_t resolution, float dChi, float Y, float radius_um, float BVF, float orientation, int32_t seed, std::string filename)
-: phantom_base(fov_um, resolution, dChi, Y, BVF, seed, filename)
+cylinder::cylinder(float fov_um, size_t resolution, float dChi, float Y, float radius_um, float volume_fraction, float orientation, int32_t seed, std::string filename)
+: phantom_base(fov_um, resolution, dChi, Y, volume_fraction, seed, filename)
 {
     set_cylinder_parameters(radius_um, orientation);
 }
@@ -52,7 +89,7 @@ bool cylinder::generate_shapes()
         BOOST_LOG_TRIVIAL(error) << "Error: The radius of the cylinder is too large for the given FOV!\n";
         return false;
     }
-    BOOST_LOG_TRIVIAL(info) << "Generating coordinates...for target BVF = " << m_BVF << "% ...\n"; 
+    BOOST_LOG_TRIVIAL(info) << "Generating coordinates...for target BVF = " << m_volume_fraction << "% ...\n"; 
     bool is_random_radius = m_radius < 0;
     float max_radius    = m_radius>0 ? m_radius:-m_radius;
     m_cylinder_points.clear();
@@ -69,37 +106,20 @@ bool cylinder::generate_shapes()
     while(progress < 100)
     {
         cyl_rad = is_random_radius ? dist(gen) * max_radius : max_radius;
-        for (size_t i = 0; i < 3; i++) // generate a random point for a sphere which fit in the FOV
+        for (size_t i = 0; i < 3; i++) // generate a random point for a cylinder which fit in the FOV
             cyl_pnt[i] = dist(gen) * (m_fov+2*cyl_rad) - cyl_rad;   
 
-        // check if sphere coordinate is ok
-        size_t c;
-        for (c=0; c<m_cylinder_radii.size(); c++)
-        {   
-            float p2p1[3];
-            subtract(cyl_pnt, m_cylinder_points[c].data(), p2p1);
-            distance = sqrtf(p2p1[0]*p2p1[0] + p2p1[1]*p2p1[1]);
-            // if the sphere is inside another sphere, generate a new sphere
-            if(distance <= m_cylinder_radii[c] ||  distance <= cyl_rad)
-                break;
-            // adjust the radius of the sphere to avoid overlap
-            if (distance < m_cylinder_radii[c] + cyl_rad)
-            {
-                if (!is_random_radius)
-                    break;            
-                cyl_rad = distance - m_cylinder_radii[c];
-            }
-        }
-        if (c < m_cylinder_radii.size())
-            continue;
+        // check if cylinder coordinate is ok
+        if (check_cylinder_overlap(m_cylinder_points, m_cylinder_radii, cyl_pnt, cyl_rad, is_random_radius)) 
+            continue;   
 
         vol_cyl = calculate_volume(cyl_pnt, cyl_rad);
         // if the total volume of the cylinders is more than the target BVF or the cylinder is outside of volume, skip this cylinder
-        if (100*(vol_cyl + vol_cyl_total) / vol_tol > 1.02*m_BVF || vol_cyl < 0)
+        if (100*(vol_cyl + vol_cyl_total) / vol_tol > 1.02*m_volume_fraction || vol_cyl < 0)
             continue;
 
         vol_cyl_total += vol_cyl;
-        progress = 100*(100.*vol_cyl_total/vol_tol/m_BVF);
+        progress = 100*(100.*vol_cyl_total/vol_tol/m_volume_fraction);
         m_cylinder_points.push_back({cyl_pnt[0], cyl_pnt[1], cyl_pnt[2]});
         m_cylinder_radii.push_back(cyl_rad);     
     }
@@ -125,7 +145,7 @@ float cylinder::calculate_volume(float *cyl_pnt, float cyl_rad)
     if (intersect == false)
         return M_PI * cyl_rad*cyl_rad * m_fov;    
 
-    // find the bounding box of the sphere
+    // find the bounding box of the cylinder
     float v_size = m_fov / m_resolution;
     float cyl_rad2  = cyl_rad*cyl_rad;
     size_t res1 = m_resolution;
@@ -198,7 +218,7 @@ bool cylinder::generate_mask_fieldmap()
         cyl_pnt_vox[1] = int32_t(cyl_pnt[1]/v_size);
         cyl_pnt_vox[2] = int32_t(cyl_pnt[2]/v_size);
 
-        // find the bounding box of the sphere
+        // find the bounding box of the cylinder
         z_min = 0;
         z_max = m_resolution;
         if (m_calc_fieldmap)
@@ -247,8 +267,8 @@ bool cylinder::generate_mask_fieldmap()
         }        
     } 
     bar->done();    
-    m_BVF = std::accumulate(m_mask.begin(), m_mask.end(), 0) * 100.0 / m_mask.size();
-    BOOST_LOG_TRIVIAL(info) << "Actual Volume Fraction = " << m_BVF << "% ...\n";   
+    m_volume_fraction = std::accumulate(m_mask.begin(), m_mask.end(), 0) * 100.0 / m_mask.size();
+    BOOST_LOG_TRIVIAL(info) << "Actual Volume Fraction = " << m_volume_fraction << "% ...\n";   
     auto end = std::chrono::high_resolution_clock::now();
     BOOST_LOG_TRIVIAL(info) << "Cylinders generated successfully! " << "Elapsed Time: " << std::chrono::duration_cast<std::chrono::seconds>(end - start).count() << " s\n";
     return true;
@@ -258,7 +278,7 @@ std::ostream& operator<<(std::ostream& os, const cylinder& obj)
 {
     os << static_cast<const phantom_base&>(obj) 
        << "  Radius: " << obj.m_radius << " \xC2\xB5m\n"
-       << "  Volume Fraction: " << obj.m_BVF << "\n"
+       << "  Volume Fraction: " << obj.m_volume_fraction << "\n"
        << "  Orientation: " << obj.m_orientation << " rad\n";
     return os;
 }

src/phantom/phantom_cylinder.h

@@ -19,7 +19,7 @@ class cylinder : public phantom_base
 {
     public:
     cylinder();
-    cylinder(float fov_um, size_t resolution, float dChi, float Y, float radius_um = 50, float BVF = 10., float orientation = -1.0f, int32_t seed=-1, std::string filename = "shape.h5");
+    cylinder(float fov_um, size_t resolution, float dChi, float Y, float radius_um = 50, float volume_fraction = 10., float orientation = -1.0f, int32_t seed=-1, std::string filename = "shape.h5");
     ~cylinder();
 
     virtual bool run() override;