Merge pull request #481 from OpenGATE/consolidate_dose_actor

Tests pass, both locally and remotely. I merge this so we can move on with the uncertainty threshold implementation in the DoseActor.

core/opengate_core/opengate_lib/GateDoseActor.h

@@ -49,90 +49,83 @@ class GateDoseActor : public GateVActor {
 
   inline void SetToWaterFlag(const bool b) { fToWaterFlag = b; }
 
-  inline bool GetSquareFlag() const { return fSquareFlag; }
+  inline bool GetEdepSquaredFlag() const { return fEdepSquaredFlag; }
 
-  inline void SetSquareFlag(const bool b) { fSquareFlag = b; }
+  inline void SetEdepSquaredFlag(const bool b) { fEdepSquaredFlag = b; }
 
-  inline void SetDensityFlag(const bool b) { fDensityFlag = b; }
+  inline void SetDoseFlag(const bool b) { fDoseFlag = b; }
 
-  inline bool GetDensityFlag() const { return fDensityFlag; }
+  inline bool GetDoseFlag() const { return fDoseFlag; }
+
+  inline void SetDoseSquaredFlag(const bool b) { fDoseSquaredFlag = b; }
+
+  inline bool GetDoseSquaredFlag() const { return fDoseSquaredFlag; }
+
+  inline void SetCountsFlag(const bool b) { fCountsFlag = b; }
+
+  inline bool GetCountsFlag() const { return fCountsFlag; }
 
   inline std::string GetPhysicalVolumeName() const {
     return fPhysicalVolumeName;
   }
 
   inline void SetPhysicalVolumeName(std::string s) { fPhysicalVolumeName = s; }
 
-  // virtual void EndSimulationAction();
-
   // Image type needs to be 3D double by default
   typedef itk::Image<double, 3> Image3DType;
 
   int sub2ind(Image3DType::IndexType index3D);
   void ind2sub(int index, Image3DType::IndexType &index3D);
-  double ComputeMeanUncertainty();
   double GetMaxValueOfImage(Image3DType::Pointer imageP);
 
   // The image is accessible on py side (shared by all threads)
   Image3DType::Pointer cpp_edep_image;
-  //  Image3DType::Pointer cpp_dose_image;
-  Image3DType::Pointer cpp_square_image;
-  Image3DType::SizeType size_edep{};
+  Image3DType::Pointer cpp_edep_squared_image;
+  Image3DType::Pointer cpp_dose_image;
+  Image3DType::Pointer cpp_dose_squared_image;
   Image3DType::Pointer cpp_density_image;
+  Image3DType::Pointer cpp_counts_image;
+  Image3DType::SizeType size_edep{};
 
-  //  // Option: indicate if we must compute uncertainty
-  //  bool fUncertaintyFlag;
+  struct threadLocalT {
+    G4EmCalculator emcalc;
+    std::vector<double> squared_worker_flatimg;
+    std::vector<int> lastid_worker_flatimg;
+  };
 
-  // Option: indicate if we must compute square
-  bool fSquareFlag{};
+  void ScoreSquaredValue(threadLocalT &data, Image3DType::Pointer cpp_image,
+                         double value, int event_id,
+                         Image3DType::IndexType index);
 
-  //  // Option: indicate if we must compute dose in Gray also
-  //  bool fDoseFlag;
+  void FlushSquaredValue(threadLocalT &data, Image3DType::Pointer cpp_image);
 
-  std::string fScoreIn;
+  void PrepareLocalDataForRun(threadLocalT &data, int numberOfVoxels);
 
   // Option: indicate we must convert to dose to water
   bool fToWaterFlag{};
 
-  // Option: indicate the density is needed
-  bool fDensityFlag{};
+  // Option: indicate if we must compute edep squared
+  bool fEdepSquaredFlag{};
 
-  //  // Option: calculate dose in stepping action. If False, calc only edep and
-  //  // divide by mass at the end of the simulation, on py side
-  //  bool fOnFlyCalcFlag;
+  // Option: Is dose to be scored?
+  bool fDoseFlag{};
+  bool fDoseSquaredFlag{};
 
-  //  // Option: cp image for each thread
-  //  bool fcpImageForThreadsFlag;
-  //
-  //  // Option: calculate the standard error of the mean
-  //  bool fSTEofMeanFlag;
-
-  // For uncertainty computation, we need temporary images
+  // Option: Are counts to be scored
+  bool fCountsFlag{};
 
   double fVoxelVolume{};
   int NbOfEvent = 0;
   int NbOfThreads = 0;
 
-  //  double goalUncertainty;
-  double threshEdepPerc{};
-  // struct timeval mTimeOfLastSaveEvent;
-
   std::string fPhysicalVolumeName;
 
   G4ThreeVector fTranslation;
   std::string fHitType;
 
 protected:
-  struct threadLocalT {
-    G4EmCalculator emcalc;
-    std::vector<double> edep_worker_flatimg;
-    std::vector<double> edepSquared_worker_flatimg;
-    std::vector<int> lastid_worker_flatimg;
-    //    int NbOfEvent_worker = 0;
-    // Image3DType::IndexType index3D;
-    // int index_flat;
-  };
-  G4Cache<threadLocalT> fThreadLocalData;
+  G4Cache<threadLocalT> fThreadLocalDataEdep;
+  G4Cache<threadLocalT> fThreadLocalDataDose;
 };
 
 #endif // GateDoseActor_h

core/opengate_core/opengate_lib/pyGateDoseActor.cpp

@@ -37,18 +37,27 @@ void init_GateDoseActor(py::module &m) {
            &GateDoseActor::BeginOfRunActionMasterThread)
       .def("EndOfRunActionMasterThread",
            &GateDoseActor::EndOfRunActionMasterThread)
-      .def("GetSquareFlag", &GateDoseActor::GetSquareFlag)
-      .def("SetSquareFlag", &GateDoseActor::SetSquareFlag)
+      .def("GetEdepSquaredFlag", &GateDoseActor::GetEdepSquaredFlag)
+      .def("SetEdepSquaredFlag", &GateDoseActor::SetEdepSquaredFlag)
+      .def("GetDoseFlag", &GateDoseActor::GetDoseFlag)
+      .def("SetDoseFlag", &GateDoseActor::SetDoseFlag)
+      .def("GetDoseSquaredFlag", &GateDoseActor::GetDoseSquaredFlag)
+      .def("SetDoseSquaredFlag", &GateDoseActor::SetDoseSquaredFlag)
       .def("GetToWaterFlag", &GateDoseActor::GetToWaterFlag)
       .def("SetToWaterFlag", &GateDoseActor::SetToWaterFlag)
-      .def("GetDensityFlag", &GateDoseActor::GetDensityFlag)
-      .def("SetDensityFlag", &GateDoseActor::SetDensityFlag)
+      .def("GetCountsFlag", &GateDoseActor::GetCountsFlag)
+      .def("SetCountsFlag", &GateDoseActor::SetCountsFlag)
       .def("GetPhysicalVolumeName", &GateDoseActor::GetPhysicalVolumeName)
       .def("SetPhysicalVolumeName", &GateDoseActor::SetPhysicalVolumeName)
       .def_readwrite("NbOfEvent", &GateDoseActor::NbOfEvent)
       .def_readwrite("cpp_edep_image", &GateDoseActor::cpp_edep_image)
+      .def_readwrite("cpp_edep_squared_image",
+                     &GateDoseActor::cpp_edep_squared_image)
+      .def_readwrite("cpp_dose_image", &GateDoseActor::cpp_dose_image)
+      .def_readwrite("cpp_dose_squared_image",
+                     &GateDoseActor::cpp_dose_squared_image)
       .def_readwrite("cpp_density_image", &GateDoseActor::cpp_density_image)
-      .def_readwrite("cpp_square_image", &GateDoseActor::cpp_square_image)
+      .def_readwrite("cpp_counts_image", &GateDoseActor::cpp_counts_image)
       .def_readwrite("fPhysicalVolumeName",
                      &GateDoseActor::fPhysicalVolumeName);
 }

opengate/actors/actoroutput.py

@@ -504,6 +504,8 @@ def set_active(self, value, item=0):
             self.data_item_config[i]["active"] = bool(value)
 
     def get_active(self, item=0):
+        if item == "any":
+            item = "all"
         items = self._collect_item_identifiers(item)
         return any([self.data_item_config[k]["active"] is True for k in items])
 
@@ -535,6 +537,7 @@ def get_item_suffix(self, item=0, **kwargs):
             )
         else:
             try:
+                # FIXME: the .get() method implicitly defines a default value, but it should not. Is this a workaround?
                 return self.data_item_config[item].get("suffix", str(item))
             except KeyError:
                 self._fatal_unknown_item(item)

opengate/actors/base.py

@@ -157,6 +157,12 @@ def __setstate__(self, state):
     #     # output_filename is a property
     #     self.known_attributes.add("output_filename")
 
+    def configure_like(self, other_obj):
+        super().configure_like(other_obj)
+        # also pick up the configuration of the user output
+        for k, v in self.user_output.items():
+            v.configure_like(other_obj.user_output[k])
+
     def to_dictionary(self):
         d = super().to_dictionary()
         d["user_output"] = dict(

opengate/actors/dataitems.py

@@ -675,9 +675,9 @@ def get_variance_or_uncertainty(self, which_quantity):
             elif self.data[1] is None or self.data[1].data is None:
                 warning(
                     "This data item does not contain squared values so no variance can be calculated. "
-                    "The variance will be set to 1 everywhere. "
+                    "The variance will be set to 0 everywhere. "
                 )
-                output_arr = np.ones_like(value_array)
+                output_arr = np.zeros_like(value_array)
             else:
                 squared_value_array = np.asarray(self.data[1].data)
                 output_arr = calculate_variance(
@@ -692,7 +692,7 @@ def get_variance_or_uncertainty(self, which_quantity):
                     output_arr = np.divide(
                         output_arr,
                         value_array / number_of_samples,
-                        out=np.ones_like(output_arr),
+                        out=np.zeros_like(output_arr),
                         where=value_array != 0,
                     )
             output_image = itk.image_view_from_array(output_arr)