Enhanced support multi head spect projection data #1579

smanwell · 2025-03-26T14:04:20Z

Changes in this pull request

This changes is related to #1184: SPECT projection DICOM: support multiple heads.
Changes have been introduced in the SPECT_dicom_to_interfile utility to accommodate such datasets.
Relevant details have been described here.

Additional changes include:

support for decompressing pixel data, if required (using gdcm API, no new dependencies introduced).
removal of unused variable and commented code

Testing performed

The code was executed on several multi-head SPECT projection datasets, including those with and without compressed pixel data.
White-box testing was performed to confirm that exception handling was working as expected in the function: SPECTDICOMData::get_proj_data(const std::string& output_file).

Related issues

This change fixes #1184.

Checklist before requesting a review

I have performed a self-review of my code
The code builds and runs on my machine (built on Visual Studio 2022, Windows 11 OS)

- Trying to enhance support for the multi-head SPECT systems. - Modifications to the open_read_file function: - set the number of projections as the number of detector heads times the number of frames per detector. - Retrieve the axial positions array from every item (i.e detector head) in the Detector Information Sequence. -- A new function GetDetectorInfo() was created in the fashion of existing methods, e.g. GetRadionuclideInfo(), to handle this.

… the SPECT_dicom_to_interfile utilitly was failing to set the full extent of rotation in the presence of multiple detector heads. The Scan Arc DICOM Tag (0018,1143), represents the "The effective angular range of the scan data" and is limited to the range of motion of only a single head. Since STIR determine the azimuthal position of the each projection view using the extent of rotation and the number of views, the total angular range among all detector heads is required. This commit introduces a change to the SPECT_dicom_to_interfile utilitly such that the extent of rotation written to the output interfile is determined as the product of the number of detectors and the scan arc. (cherry picked from commit 69178f599610bbb14db5471f636ecf204aa61c54) (cherry picked from commit 63064ff)

… DICOM data. - It now detects whether the data are compressed, and if so it rewrites the file using a non-compressed transfer syntax UID and re-reads the file. - The get_proj_data function is now recursive as it calls itself to re-read the non-compressed DICOM projection data.

…nt smart pointers for calls to gdcm::DataElement::GertValueAsSQ()

… on to the latest STIR master. - remove duplicate code in SPECTDICOMData::open_dicom_file() -- Gave precedence to my changes related to setting the num_of_projections variable which accounts for multiple detector heads -- Gave precedence to my changes related to setting the rotation_radius variable, which accounts for multiple detector heads - remove unecessary { } blocks in SPECTDICOMData::open_dicom_file() - moved the DICOM reads associated with isotope name, actual frame duration (string), and energy windows outside of the if (!is_planar) block. Consistent with STIR master.

- modified SPECTDICOMData::open_proj_data() such that it no longer needs to write the uncompressed pixel data to a file. The change of DICOM transfer syntax is managed in memory. As a result this function no longer calls itself recursively when decompression is required. - removed the optional input argument from the SPECTDICOMData::open_proj_data() method, no longer needed due to previous change. - refactorerd the code block in SPECTDICOMData::open_proj_data() that was responsible for populate the float vector - pixel_data_as_float, into a separate function as it is now invoked in separate if-blocks depending on whether decompression was required. - included initialization values for several more member variavles of the SPECTDICOMData class. - removed the unused member variable SPECTDICOMData::num_dimensions.

KrisThielemans

some code clean-up. I haven't checked actual logic...

KrisThielemans · 2025-05-12T15:08:03Z

src/utilities/SPECT_dicom_to_interfile.cxx

-  stir::Succeeded open_dicom_file();
+  stir::Succeeded get_proj_data(const std::string& output_file ) const;
+  stir::Succeeded open_dicom_file(bool is_planar);
+  std::vector<float> store_pixel_buffer_as_float_vector(const gdcm::ByteValue*& bv) const;


why a reference? Can this not simply be

Suggested change

std::vector<float> store_pixel_buffer_as_float_vector(const gdcm::ByteValue*& bv) const;

std::vector<float> store_pixel_buffer_as_float_vector(const gdcm::ByteValue* bv) const;

KrisThielemans · 2025-05-12T15:10:50Z

src/utilities/SPECT_dicom_to_interfile.cxx

            {
              const gdcm::DataElement& de = file.GetDataSet().GetDataElement(t);
-              const gdcm::SequenceOfItems* sqi = de.GetValueAsSQ();
+              const gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = de.GetValueAsSQ();


the following should work and be less dependent on GDCM version

Suggested change

const gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = de.GetValueAsSQ();

const auto sqi = de.GetValueAsSQ();

this would now work for de as well (as we require C++-17 by now)

const auto& de = file.GetDataSet().GetDataElement(t);

KrisThielemans · 2025-05-12T15:13:39Z

src/utilities/SPECT_dicom_to_interfile.cxx

      // Get Energy Window Info Sequence
      const gdcm::DataElement& de = file.GetDataSet().GetDataElement(energy_window_info_seq);
-      const gdcm::SequenceOfItems* sqi = de.GetValueAsSQ();
+      const gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = de.GetValueAsSQ();


KrisThielemans · 2025-05-12T15:13:45Z

src/utilities/SPECT_dicom_to_interfile.cxx

      // Get Energy Window Range Sequence
      const gdcm::DataElement& element = item.GetDataElement(energy_window_range_seq);
-      const gdcm::SequenceOfItems* sqi2 = element.GetValueAsSQ();
+      const gdcm::SmartPointer<gdcm::SequenceOfItems> sqi2 = element.GetValueAsSQ();


as above, I'll stop flagging remaining ones

KrisThielemans · 2025-05-12T15:20:46Z

src/utilities/SPECT_dicom_to_interfile.cxx

      actual_frame_duration = std::stoi(actual_frame_duration_as_string);
    }

-  {


Please keep this {} level if you can. It's there to "isolate" bits of code. E.g. the str is a local variable, but also make it clearer to the reader (hopefully) that this is a bit of code initialising the num_energy_windows code (only)

KrisThielemans · 2025-05-12T15:23:04Z

src/utilities/SPECT_dicom_to_interfile.cxx

-      std::replace(rotation_radius.begin(), rotation_radius.end(), slash, comma);
-    }
+      // Set the extent of rotation as the product of the scan arc and the number of detectors.
+      // Note that in the presence of multiple detector heads, this assumes that the angular extent of


truncated sentence

KrisThielemans · 2025-05-12T15:36:37Z

src/utilities/SPECT_dicom_to_interfile.cxx

+  std::vector<float> pixel_data_as_float;
+  uint64_t n_elements = (uint64_t)bv->GetLength() / sizeof(uint16_t);
+  uint16_t* ptr = (uint16_t*)bv->GetPointer();


this would be rather slow. Best to pre-allocate. Also, best to use C++ casts

Suggested change

std::vector<float> pixel_data_as_float;

uint64_t n_elements = (uint64_t)bv->GetLength() / sizeof(uint16_t);

uint16_t* ptr = (uint16_t*)bv->GetPointer();

const auto n_elements = static_cast<uint64_t>(bv->GetLength()) / sizeof(uint16_t);

std::vector<float> pixel_data_as_float(n_elements);

auto ptr = reinterpret_cast<uint16_t const*>(bv->GetPointer());

auto to_float = [](const int i) { return float{i}; };

std::transform(ptr, ptr+n_elements, pixel_data_as_float.begin(), to_float);

would do the whole thing I believe

KrisThielemans · 2025-05-12T15:37:48Z

src/utilities/SPECT_dicom_to_interfile.cxx

        {
-          stir::error(boost::format("SPECTDICOMData: cannot read file %1%") % dicom_filename);
-          // return stir::Succeeded::no;
+          stir::error(boost::format("in SPECTDICOMData::get_proj_data() - cannot read file %1%") % dicom_filename);


Suggested change

stir::error(boost::format("in SPECTDICOMData::get_proj_data() - cannot read file %1%") % dicom_filename);

stir::error("in SPECTDICOMData::get_proj_data() - cannot read file " + dicom_filename);

avoiding boost::format as much as possible (see the fmt PR which gets rid of it)

danieldeidda · 2025-05-22T13:14:18Z

BTW @smanwell I notice recently with some siemens data that was acquired with energy windows for scatter estimation, that every other sinogram is saved as if there was a change of rotation orientation. Did you have the same issue?

smanwell · 2025-05-22T13:25:37Z

@KrisThielemans Thanks for your feedback, I'll incorporate your suggestions as soon as I can return my focus to this project.

smanwell · 2025-05-22T13:28:52Z

@danieldeidda I'm not quite sure what you mean. Is it possible to provide a screen capture of some of the DICOM fields that are involved?

danieldeidda · 2025-05-22T13:22:24Z

src/utilities/SPECT_dicom_to_interfile.cxx

  std::string start_angle_as_string;
  std::string angular_step_as_string;
-  std::string extent_of_rotation_as_string;
+  std::string scan_arc_as_string;


why the change of name? isn't the old one the same in the DICOM tag?

actually nevermind this is not a dicom tag

danieldeidda · 2025-05-22T13:44:59Z

src/utilities/SPECT_dicom_to_interfile.cxx

+        }
+
+      if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1143), extent_of_rotation_as_string) == stir::Succeeded::yes)
+        {


is this not a repetition then, you have it already in scan arc

danieldeidda · 2025-05-22T13:48:06Z

src/utilities/SPECT_dicom_to_interfile.cxx

+      // Assert that the factors were valid before calculating the extent of rotation.
+      if (num_of_detectors > 0 && scan_arc > 0)
+        {
+          extent_of_rotation = num_of_detectors * scan_arc;


i am not sure this is correct. The definition of scan arc in the dicom standard is the following: The effective angular range of the scan data in degrees. The value shall be positive

this makes me think that scan arc=num detectors* number of projection frames. In all the data I have this is 360 degrees. So scan arc =extent of rotation

Before I made changes to this code, the extent of rotation was equated to the DICOM scan arc. For the dual-head datasets that I had used for testing, the scan arc (and thus the extent of rotation) was simply 180 degrees. As a result, all of the projection images were interpreted as belonging to erroneous angular positions and the reconstructed had severe artifacts.

Take a look at a dataset with multiple (dual) heads. If, for example, the detectors are separated by 180 deg wrt each other, each will have a scan arc of 180 deg. But the full extent of the rotation is given as the product of the number of detectors with the scan arc. For a single headed system, then yes I'd expect the scan arc to be 360 deg and equal to the extent of rotation.

If there were three heads and each was positioned at equal angles wrt each other - 120 deg. the full extent of rotation, once again, would be given as scan arc * number of detectors.

I've not come across cases with more complex detector arrangements, so I haven't seen if my change holds for all cases.

smanwell added 6 commits March 24, 2025 10:14

For SPECT_dicom_to_interfile utility, replace const pointers with cos…

6c449fe

…nt smart pointers for calls to gdcm::DataElement::GertValueAsSQ()

KrisThielemans requested changes May 12, 2025

View reviewed changes

danieldeidda reviewed May 22, 2025

View reviewed changes

	std::vector<float> store_pixel_buffer_as_float_vector(const gdcm::ByteValue*& bv) const;
	std::vector<float> store_pixel_buffer_as_float_vector(const gdcm::ByteValue* bv) const;

	const gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = de.GetValueAsSQ();
	const auto sqi = de.GetValueAsSQ();

-  std::vector<float> pixel_data_as_float;
-  uint64_t n_elements = (uint64_t)bv->GetLength() / sizeof(uint16_t);
-  uint16_t* ptr = (uint16_t*)bv->GetPointer();
+  const auto n_elements = static_cast<uint64_t>(bv->GetLength()) / sizeof(uint16_t);
+  std::vector<float> pixel_data_as_float(n_elements);
+  auto ptr = reinterpret_cast<uint16_t const*>(bv->GetPointer());
+  auto to_float = [](const int i) { return float{i}; };
+  std::transform(ptr, ptr+n_elements, pixel_data_as_float.begin(), to_float);

	stir::error(boost::format("in SPECTDICOMData::get_proj_data() - cannot read file %1%") % dicom_filename);
	stir::error("in SPECTDICOMData::get_proj_data() - cannot read file " + dicom_filename);

Enhanced support multi head spect projection data #1579

Are you sure you want to change the base?

Enhanced support multi head spect projection data #1579

Uh oh!

Conversation

smanwell commented Mar 26, 2025

Changes in this pull request

Testing performed

Related issues

Checklist before requesting a review

Uh oh!

KrisThielemans left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KrisThielemans May 12, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

danieldeidda commented May 22, 2025

Uh oh!

smanwell commented May 22, 2025

Uh oh!

smanwell commented May 22, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

KrisThielemans May 12, 2025 •

edited

Loading