Switch image module to SimpleITK - part 1

core/opengate_core/opengate_lib/helpers_itk_image_py.h

@@ -13,9 +13,60 @@
 #include <pybind11/stl_bind.h>
 
 #include "itkImage.h"
+#include "itkImageRegionConstIterator.h"
 #include "itkImportImageFilter.h"
 #include "itkSmartPointer.h"
 
+void print_sum_of_array(py::array array) {
+  // (function for debug only)
+  // Get the size of the array (number of elements)
+  ssize_t num_elements = array.size();
+
+  // Determine the data type of the array
+  py::buffer_info buf_info = array.request();
+  const void *data_ptr = buf_info.ptr;
+
+  // Initialize the sum variable
+  double sum = 0.0;
+  std::cout << "num_elements = " << num_elements << std::endl;
+  // Calculate the sum based on the data type
+  if (buf_info.format == py::format_descriptor<unsigned char>::format()) {
+    const auto *data = static_cast<const unsigned char *>(data_ptr);
+    for (ssize_t i = 0; i < num_elements; ++i) {
+      if (i == 0)
+        std::cout << "uc data = " << data[i] << std::endl;
+      sum += data[i];
+    }
+  } else if (buf_info.format ==
+             py::format_descriptor<unsigned short>::format()) {
+    const auto *data = static_cast<const unsigned short *>(data_ptr);
+    for (ssize_t i = 0; i < num_elements; ++i) {
+      if (i == 0)
+        std::cout << "us data = " << data[i] << std::endl;
+      sum += data[i];
+    }
+  } else if (buf_info.format == py::format_descriptor<float>::format()) {
+    const auto *data = static_cast<const float *>(data_ptr);
+    for (ssize_t i = 0; i < num_elements; ++i) {
+      if (i == 0)
+        std::cout << "f data = " << data[i] << std::endl;
+      sum += data[i];
+    }
+  } else if (buf_info.format == py::format_descriptor<double>::format()) {
+    const auto *data = static_cast<const double *>(data_ptr);
+    for (ssize_t i = 0; i < num_elements; ++i) {
+      if (i == 0)
+        std::cout << "d data = " << data[i] << std::endl;
+      sum += data[i];
+    }
+  } else {
+    throw std::runtime_error("Unsupported data type");
+  }
+
+  // Print the sum of the array elements
+  std::cout << "Sum of array elements: " << sum << std::endl;
+}
+
 /** For now make copies of data to numpy arrays.
  * Some examples (pybind11 source code docs are non-existant)
  * From: https://github.com/pybind/pybind11/issues/323
@@ -177,9 +228,37 @@ void declare_itk_image_ptr(pybind11::module &m, const std::string &typestr) {
                 (contiguous == "F")
                     ? std::vector<size_t>{size[2], size[1], size[0]}
                     : std::vector<size_t>{size[0], size[1], size[2]};
-            return py::array(
+
+            // Calculate strides
+            using PixelType = typename TImagePointer::ObjectType::PixelType;
+            std::vector<unsigned long> strides;
+            if (contiguous == "F") {
+              strides = {
+                  sizeof(PixelType), // Stride for the first dimension (z-axis)
+                  sizeof(PixelType) *
+                      size[2], // Stride for the second dimension (y-axis)
+                  sizeof(PixelType) * size[2] *
+                      size[1] // Stride for the third dimension (x-axis)
+              };
+            } else {
+              strides = {
+                  sizeof(PixelType) * size[1] *
+                      size[2], // Stride for the first dimension (x-axis)
+                  sizeof(PixelType) *
+                      size[2],      // Stride for the second dimension (y-axis)
+                  sizeof(PixelType) // Stride for the third dimension (z-axis)
+              };
+            }
+
+            auto a = py::array(
                 py::dtype::of<typename TImagePointer::ObjectType::PixelType>(),
-                shape, img->GetBufferPointer());
+                shape, strides, img->GetBufferPointer());
+
+            // print_sum_of_array(a);
+            // a = a.attr("copy")();
+            // print_sum_of_array(a);
+
+            return a;
           },
           py::arg("contig") = "F")
       // TODO: Create a view (non-copy) of the data

opengate/actors/digitizers.py

@@ -12,7 +12,7 @@
 from ..image import (
     align_image_with_physical_volume,
     update_image_py_to_cpp,
-    get_cpp_image,
+    update_image_cpp_to_py,
     get_info_from_image,
     create_3d_image,
     write_itk_image,
@@ -637,7 +637,7 @@ def StartSimulationAction(self):
     def EndSimulationAction(self):
         g4.GateDigitizerProjectionActor.EndSimulationAction(self)
         # retrieve the image
-        self.output_image = get_cpp_image(self.fImage)
+        self.output_image = update_image_cpp_to_py(self.fImage)
         # put back the origin
         self.output_image.SetOrigin(self.start_output_origin)
         info = get_info_from_image(self.output_image)

opengate/actors/doseactors.py

@@ -1,27 +1,12 @@
-import itk
-import numpy as np
 import opengate_core as g4
 from .base import ActorBase
-from ..exception import fatal, warning
+from ..exception import warning
 from ..utility import (
     g4_units,
-    ensure_filename_is_str,
     standard_error_c4_correction,
 )
-from ..image import (
-    create_3d_image,
-    align_image_with_physical_volume,
-    update_image_py_to_cpp,
-    create_image_like,
-    get_info_from_image,
-    get_origin_wrt_images_g4_position,
-    get_cpp_image,
-    itk_image_view_from_array,
-    divide_itk_images,
-    scale_itk_image,
-    write_itk_image,
-)
-from ..geometry.materials import create_mass_img, create_density_img
+from ..image import *
+from ..geometry.materials import create_density_img
 
 
 class DoseActor(g4.GateDoseActor, ActorBase):
@@ -160,7 +145,9 @@ def initialize(self, volume_engine=None):
         # create itk image (py side)
         size = np.array(self.user_info.size)
         spacing = np.array(self.user_info.spacing)
-        self.py_edep_image = create_3d_image(size, spacing, pixel_type="double")
+        self.py_edep_image = create_3d_image(
+            size, spacing, pixel_type="double", fill_value=0
+        )
         # compute the center, using translation and half pixel spacing
         self.img_origin_during_run = (
             -size * spacing / 2.0 + spacing / 2.0 + self.user_info.translation
@@ -208,7 +195,7 @@ def StartSimulationAction(self):
             or self.user_info.ste_of_mean
         ):
             self.py_square_image = create_image_like(
-                self.py_edep_image, pixel_type="double"
+                self.py_edep_image, pixel_type="double", fill_value=0
             )
             update_image_py_to_cpp(
                 self.py_square_image, self.cpp_square_image, self.first_run
@@ -259,7 +246,7 @@ def EndSimulationAction(self):
 
         # Get the itk image from the cpp side
         # Currently a copy. Maybe later as_pyarray ?
-        self.py_edep_image = get_cpp_image(self.cpp_edep_image)
+        self.py_edep_image = update_image_cpp_to_py(self.cpp_edep_image)
 
         # set the property of the output image:
         # in the coordinate system of the attached volume
@@ -355,8 +342,8 @@ def compute_dose_from_edep_img(self):
             )
 
     def fetch_square_image_from_cpp(self):
-        if self.py_square_image == None:
-            self.py_square_image = get_cpp_image(self.cpp_square_image)
+        if self.py_square_image is None:
+            self.py_square_image = update_image_cpp_to_py(self.cpp_square_image)
             self.py_square_image.SetOrigin(self.output_origin)
             self.py_square_image.CopyInformation(self.py_edep_image)
 
@@ -391,17 +378,17 @@ def create_uncertainty_img(self):
 
         self.fetch_square_image_from_cpp()
 
-        edep = itk.array_view_from_image(self.py_edep_image)
-        square = itk.array_view_from_image(self.py_square_image)
+        edep = itk_array_view_from_image(self.py_edep_image)
+        square = itk_array_view_from_image(self.py_square_image)
 
-        self.py_edep_image_tmp = itk_image_view_from_array(edep)
-        self.py_edep_image_tmp.CopyInformation(self.py_edep_image)
-        self.py_edep_image = self.py_edep_image_tmp
-        del self.py_edep_image_tmp
+        # self.py_edep_image_tmp = itk_image_view_from_array(edep)
+        # self.py_edep_image_tmp.CopyInformation(self.py_edep_image)
+        # self.py_edep_image = self.py_edep_image_tmp
+        # del self.py_edep_image_tmp
 
         # uncertainty image
         self.uncertainty_image = create_image_like(
-            self.py_edep_image, pixel_type="double"
+            self.py_edep_image, pixel_type="double", fill_value=0
         )
         # unc = itk.array_view_from_image(self.uncertainty_image)
 
@@ -616,8 +603,8 @@ def EndSimulationAction(self):
 
         # Get the itk image from the cpp side
         # Currently a copy. Maybe later as_pyarray ?
-        self.py_numerator_image = get_cpp_image(self.cpp_numerator_image)
-        self.py_denominator_image = get_cpp_image(self.cpp_denominator_image)
+        self.py_numerator_image = update_image_cpp_to_py(self.cpp_numerator_image)
+        self.py_denominator_image = update_image_cpp_to_py(self.cpp_denominator_image)
 
         # set the property of the output image:
         # in the coordinate system of the attached volume
@@ -757,7 +744,7 @@ def EndSimulationAction(self):
 
         # Get the itk image from the cpp side
         # Currently a copy. Maybe later as_pyarray ?
-        self.py_fluence_image = get_cpp_image(self.cpp_fluence_image)
+        self.py_fluence_image = update_image_cpp_to_py(self.cpp_fluence_image)
 
         # set the property of the output image:
         origin = self.img_origin_during_run