Merge pull request #1503 from KrisThielemans/GERDF9_DMI

GE Discovery MI updates

documentation/release_6.3.htm

@@ -25,6 +25,11 @@ <h3>New functionality</h3>
     is disabled for backwards compatibility.<br>
     <a href=https://github.com/UCL/STIR/pull/1291>PR #1291</a>
   </li>
+  <li>
+    Data from GE Discovery MI systems in RDF9 should now be readable. TOF information on these scanners has also been added.
+    However, projection data is currently still always returned as non-TOF (but list-mode data is read as TOF).<br>
+    <a href=https://github.com/UCL/STIR/pull/1503>PR #1503</a>
+  </li>
 </ul>
 
 

src/IO/GEHDF5Wrapper.cxx

@@ -65,6 +65,17 @@ read_float(const H5::H5File& file, const std::string& dataset)
   return tmp;
 }
 
+static std::string
+read_string(const H5::H5File& file, const std::string& dataset)
+{
+  H5::DataSet ds = file.openDataSet(dataset.c_str());
+  H5::StrType datatype = ds.getStrType();
+
+  std::string value;
+  ds.read(value, datatype);
+  return value;
+}
+
 bool
 GEHDF5Wrapper::check_GE_signature(const std::string& filename)
 {
@@ -90,13 +101,7 @@ GEHDF5Wrapper::check_GE_signature(H5::H5File& file)
   if (file.getId() == -1)
     error("File is not open. Aborting");
 
-  H5::StrType vlst(
-      0,
-      37); // 37 here is the length of the string (PW got it from the text file generated by list2txt with the LIST000_decomp.BLF
-  std::string read_str_manufacturer;
-
-  H5::DataSet dataset2 = file.openDataSet("/HeaderData/ExamData/manufacturer");
-  dataset2.read(read_str_manufacturer, vlst);
+  std::string read_str_manufacturer = read_string(file, "/HeaderData/ExamData/manufacturer");
 
   if (read_str_manufacturer == "GE MEDICAL SYSTEMS")
     {
@@ -277,14 +282,7 @@ GEHDF5Wrapper::open(const std::string& filename)
 shared_ptr<Scanner>
 GEHDF5Wrapper::get_scanner_from_HDF5()
 {
-  std::string read_str_scanner;
-  H5::StrType vlst(
-      0,
-      37); // 37 here is the length of the string (PW got it from the text file generated by list2txt with the LIST000_decomp.BLF
-
-  H5::DataSet dataset = file.openDataSet("/HeaderData/ExamData/scannerDesc");
-  dataset.read(read_str_scanner, vlst);
-
+  std::string read_str_scanner = read_string(this->file, "/HeaderData/ExamData/scannerDesc");
   float effective_ring_diameter;
   int num_transaxial_blocks_per_bucket = 0;
   int num_axial_blocks_per_bucket = 0;
@@ -425,6 +423,11 @@ GEHDF5Wrapper::initialise_proj_data_info_from_HDF5()
 {
   shared_ptr<Scanner> scanner_sptr = get_scanner_from_HDF5();
 
+  // TODO get TOF mashing when reading sinos as TOF
+  const auto num_tof_bins = read_dataset_uint32("/HeaderData/Sorter/numTOF_bins");
+  if (num_tof_bins > 1)
+    warning("GE RDF data currently still read as non-TOF");
+
   this->proj_data_info_sptr
       = ProjDataInfo::construct_proj_data_info(scanner_sptr,
                                                /*span*/ 2,
@@ -654,11 +657,20 @@ GEHDF5Wrapper::initialise_proj_data(const unsigned int view_num)
   if (view_num == 0 || view_num > static_cast<unsigned>(this->get_scanner_sptr()->get_num_detectors_per_ring() / 2))
     error("internal error in GE HDF5 code: view number " + std::to_string(view_num) + " is incorrect");
 
+  const auto num_tof_bins = read_dataset_uint32("/HeaderData/Sorter/numTOF_bins");
+
   if (rdf_ver == 9)
     {
-      m_address = "/SegmentData/Segment2/3D_TOF_Sinogram/view" + std::to_string(view_num);
-
-      m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address)));
+      if (num_tof_bins > 1)
+        {
+          m_address = "/SegmentData/Segment2/3D_TOF_Sinogram/view" + std::to_string(view_num);
+          m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address)));
+        }
+      else
+        {
+          m_address = "/SegmentData/Segment2/3D_Sinogram/view" + std::to_string(view_num);
+          m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address)));
+        }
       m_dataspace = m_dataset_sptr->getSpace();
       // Create an array to host the size of the dimensions
       const int rank = m_dataspace.getSimpleExtentNdims();
@@ -673,13 +685,7 @@ GEHDF5Wrapper::initialise_proj_data(const unsigned int view_num)
       // AB for signa, these where [1981,27,357] and [45,448,357]
       m_NX_SUB = dims[0]; // hyperslab dimensions
       m_NY_SUB = dims[1];
-      m_NZ_SUB = dims[2];
-#if 0
-        // AB todo: ??? why are these different?
-        m_NX = 45;        // output buffer dimensions
-        m_NY = 448;
-        m_NZ = 357;
-#endif
+      m_NZ_SUB = rank > 2 ? dims[2] : 1;
     }
   else
     return Succeeded::no;
@@ -831,7 +837,6 @@ GEHDF5Wrapper::read_sinogram(Array<3, unsigned char>& output,
                 }
             }
         }
-      output.release_data_ptr();
     }
 
   return Succeeded::yes;

src/buildblock/Scanner.cxx

@@ -864,9 +864,10 @@ Scanner::Scanner(Type scanner_type)
                  1,
                  0.0944F, // energy resolution from Hsu et al. 2017
                  511.F,
-                 (short int)(0),
-                 (float)(0), // TODO
-                 (float)(0));
+                 (short int)(2 * 188 + 1), // from RDF file
+                 (float)(13),              // from RDF file
+                 (float)(375.4)            // Hsu et al.
+      );
       break;
 
     case DiscoveryMI4ring: // This is the 4-ring DMI
@@ -895,9 +896,10 @@ Scanner::Scanner(Type scanner_type)
                  1,
                  0.0944F, // energy resolution from Hsu et al. 2017
                  511.F,
-                 (short int)(0),
-                 (float)(0), // TODO
-                 (float)(0));
+                 (short int)(2 * 188 + 1), // from RDF file
+                 (float)(13),              // from RDF file
+                 (float)(375.4)            // Hsu et al.
+      );
       break;
 
     case DiscoveryMI5ring: // This is the 5-ring DMI
@@ -926,9 +928,42 @@ Scanner::Scanner(Type scanner_type)
                  1,
                  0.0944F, // energy resolution from Hsu et al. 2017
                  511.F,
-                 (short int)(0),
-                 (float)(0), // TODO
-                 (float)(0));
+                 (short int)(2 * 188 + 1), // from RDF file
+                 (float)(13),              // from RDF file
+                 (float)(375.4)            // Hsu et al.
+      );
+      break;
+
+    case DiscoveryMI6ring: // This is the 6-ring DMI
+      // as above, but one extra block
+      // Hsu et al. 2017 JNM
+      // crystal size 3.95 x 5.3 x 25
+      set_params(DiscoveryMI6ring,
+                 string_list("GE Discovery MI 6 rings",
+                             "Discovery MI6",
+                             "Discovery MI"), // needs to include last value as used by GE in RDF files
+                 54,
+                 415,
+                 401, // TODO should compute num_arccorrected_bins from effective_FOV/default_bin_size
+                 2 * 272,
+                 380.5F - 9.4F, // TODO inner_ring_radius and DOI, currently set such that effective ring-radius is correct
+                 9.4F,          // TODO DOI
+                 5.52296F,      // ring-spacing
+                 2.206F,        // TODO currently using the central bin size default bin size. GE might be using something else
+                 static_cast<float>(-4.399 * _PI / 180), // TODO check sign
+                 6,
+                 4,
+                 9,
+                 4,
+                 1,
+                 1,
+                 1,
+                 0.0944F, // energy resolution from Hsu et al. 2017
+                 511.F,
+                 (short int)(2 * 188 + 1), // from RDF file
+                 (float)(13),              // from RDF file
+                 (float)(375.4)            // Hsu et al.
+      );
       break;
     case HZLR:
 

src/include/stir/Scanner.h

@@ -165,6 +165,7 @@ class Scanner
     DiscoveryMI3ring,
     DiscoveryMI4ring,
     DiscoveryMI5ring,
+    DiscoveryMI6ring,
     HZLR,
     RATPET,
     PANDA,