From dd81ff17c41317f1db6720d8b3c6e4326f004163 Mon Sep 17 00:00:00 2001
From: Kris Thielemans <k.thielemans@ucl.ac.uk>
Date: Thu, 18 Jul 2024 17:14:07 +0100
Subject: [PATCH 1/2] SPECT DICOM: safer reading
- do checks if sequences are present to avoid seg-faults
- read more data if planar, but avoid reading data that doesn't make sense
---
src/utilities/SPECT_dicom_to_interfile.cxx | 122 ++++++++++++---------
1 file changed, 73 insertions(+), 49 deletions(-)
diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx
index 356550ecf6..ad08f258c2 100644
--- a/src/utilities/SPECT_dicom_to_interfile.cxx
+++ b/src/utilities/SPECT_dicom_to_interfile.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2018, 2023, University College London
+ Copyright (C) 2018, 2023-2024, University College London
Copyright (C) 2019-2023, National Physical Laboratory
This file is part of STIR.
@@ -129,6 +129,18 @@ GetDICOMTagInfo(const gdcm::File& file, const gdcm::Tag tag, std::string& dst, c
{
const gdcm::DataElement& de = file.GetDataSet().GetDataElement(t);
const gdcm::SequenceOfItems* sqi = de.GetValueAsSQ();
+ if (!sqi)
+ {
+ // try next sequence
+ continue;
+ }
+ const auto sqi_length = sqi->GetNumberOfItems();
+ if (static_cast<unsigned>(sequence_idx) > sqi_length)
+ {
+ // stir::warning("sequence_id larger than length");
+ // try next sequence
+ continue;
+ }
const gdcm::Item& item = sqi->GetItem(sequence_idx);
element = item.GetDataElement(tag);
@@ -317,6 +329,57 @@ SPECTDICOMData::open_dicom_file(bool is_planar)
}
this->num_of_projections = 1;
+ this->calibration_factor = -1;
+ if (GetRadionuclideInfo(file, RadionuclideInfo::CodeMeaning, isotope_name) == stir::Succeeded::yes)
+ {
+ std::cout << "Isotope name: " << isotope_name << std::endl;
+ }
+
+ if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1242), actual_frame_duration_as_string) == stir::Succeeded::yes)
+ {
+ actual_frame_duration = std::stoi(actual_frame_duration_as_string);
+ }
+
+ {
+ this->num_energy_windows = 0;
+ std::string str;
+ if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x0011), str) == stir::Succeeded::yes)
+ this->num_energy_windows = std::stoi(str);
+ }
+ if (this->num_energy_windows == 0)
+ {
+ std::cout << "No energy window information found\n";
+ }
+ else
+ {
+ energy_window_name.resize(num_energy_windows);
+ lower_en_window_thres.resize(num_energy_windows);
+ upper_en_window_thres.resize(num_energy_windows);
+ for (int w = 1; w <= num_energy_windows; ++w)
+ {
+ if (GetEnergyWindowInfo(file, EnergyWindowInfo::WindowName, energy_window_name[w - 1], w) == stir::Succeeded::yes)
+ {
+ std::cout << "Energy window: " << energy_window_name[w - 1] << std::endl;
+ }
+ else
+ energy_window_name[w - 1] = "en" + std::to_string(w);
+
+ if (GetEnergyWindowInfo(file, EnergyWindowInfo::LowerThreshold, lower_window_as_string, w) == stir::Succeeded::yes)
+ {
+ lower_en_window_thres[w - 1] = std::stof(lower_window_as_string);
+ std::cout << "Lower energy window limit: " << std::fixed << std::setprecision(6) << lower_en_window_thres[w - 1]
+ << std::endl;
+ }
+
+ if (GetEnergyWindowInfo(file, EnergyWindowInfo::UpperThreshold, upper_window_as_string, w) == stir::Succeeded::yes)
+ {
+ upper_en_window_thres[w - 1] = std::stof(upper_window_as_string);
+ std::cout << "Upper energy window limit: " << std::fixed << std::setprecision(6) << upper_en_window_thres[w - 1]
+ << std::endl;
+ }
+ }
+ }
+
if (!is_planar)
{
if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x0053), no_of_proj_as_str) == stir::Succeeded::yes)
@@ -327,7 +390,6 @@ SPECTDICOMData::open_dicom_file(bool is_planar)
if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1140), direction_of_rotation) == stir::Succeeded::yes)
{
-
if (direction_of_rotation == "CC")
direction_of_rotation = "CCW";
@@ -345,13 +407,6 @@ SPECTDICOMData::open_dicom_file(bool is_planar)
calibration_factor = std::stof(calib_factor_as_string);
std::cout << "calibration factor: " << std::fixed << std::setprecision(6) << calibration_factor << std::endl;
}
- else
- calibration_factor = -1;
-
- if (GetRadionuclideInfo(file, RadionuclideInfo::CodeMeaning, isotope_name) == stir::Succeeded::yes)
- {
- std::cout << "Isotope name: " << isotope_name << std::endl;
- }
if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1144), angular_step_as_string) == stir::Succeeded::yes)
{
@@ -373,43 +428,6 @@ SPECTDICOMData::open_dicom_file(bool is_planar)
std::cout << "Radius: " << radius_as_string << " " << slash << std::endl;
std::replace(rotation_radius.begin(), rotation_radius.end(), slash, comma);
}
-
- if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1242), actual_frame_duration_as_string) == stir::Succeeded::yes)
- {
- actual_frame_duration = std::stoi(actual_frame_duration_as_string);
- }
-
- {
- std::string str;
- if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x0011), str) == stir::Succeeded::yes)
- num_energy_windows = std::stoi(str);
- }
- energy_window_name.resize(num_energy_windows);
- lower_en_window_thres.resize(num_energy_windows);
- upper_en_window_thres.resize(num_energy_windows);
- for (int w = 1; w <= num_energy_windows; ++w)
- {
- if (GetEnergyWindowInfo(file, EnergyWindowInfo::WindowName, energy_window_name[w - 1], w) == stir::Succeeded::yes)
- {
- std::cout << "Energy window: " << energy_window_name[w - 1] << std::endl;
- }
- else
- energy_window_name[w - 1] = "en" + std::to_string(w);
-
- if (GetEnergyWindowInfo(file, EnergyWindowInfo::LowerThreshold, lower_window_as_string, w) == stir::Succeeded::yes)
- {
- lower_en_window_thres[w - 1] = std::stof(lower_window_as_string);
- std::cout << "Lower energy window limit: " << std::fixed << std::setprecision(6) << lower_en_window_thres[w - 1]
- << std::endl;
- }
-
- if (GetEnergyWindowInfo(file, EnergyWindowInfo::UpperThreshold, upper_window_as_string, w) == stir::Succeeded::yes)
- {
- upper_en_window_thres[w - 1] = std::stof(upper_window_as_string);
- std::cout << "Upper energy window limit: " << std::fixed << std::setprecision(6) << upper_en_window_thres[w - 1]
- << std::endl;
- }
- }
}
num_dimensions = 2;
@@ -497,12 +515,14 @@ SPECTDICOMData::get_interfile_header(std::string& output_header, const std::stri
ss << "!number of projections := " << this->num_of_projections << std::endl;
ss << "number of time frames := 1" << std::endl;
ss << "!image duration (sec)[1] := " << this->num_of_projections * this->actual_frame_duration / 1000 << std::endl;
- ss << "!extent of rotation := " << this->extent_of_rotation << std::endl;
+ if (!is_planar)
+ ss << "!extent of rotation := " << this->extent_of_rotation << std::endl;
ss << "!process status := acquired" << std::endl;
ss << std::endl;
ss << "!SPECT STUDY (acquired data) :=" << std::endl;
- ss << "!direction of rotation := " << this->direction_of_rotation << std::endl;
+ if (!is_planar)
+ ss << "!direction of rotation := " << this->direction_of_rotation << std::endl;
ss << "start angle := " << this->start_angle << std::endl;
if (is_planar)
@@ -558,7 +578,11 @@ SPECTDICOMData::get_proj_data(const std::string& output_file) const
const gdcm::DataElement& de = file.GetDataSet().GetDataElement(gdcm::Tag(0x7fe0, 0x0010));
const gdcm::ByteValue* bv = de.GetByteValue();
-
+ if (!bv)
+ {
+ stir::warning("GDCM GetByteValue failed on x7fe0, 0x0010");
+ return stir::Succeeded::no;
+ }
/*
std::string tmpFile = "tmp.s";
std::ofstream outfile(tmpFile.c_str(), std::ios::out | std::ios::binary);
From 844a02ad6a5bf9798509afa475bcbabadaeefe43 Mon Sep 17 00:00:00 2001
From: Kris Thielemans <k.thielemans@ucl.ac.uk>
Date: Fri, 19 Jul 2024 15:23:21 +0100
Subject: [PATCH 2/2] SPECT DICOM: remove requirement for is_planar
As STIR can only read SPECT sinograms, we now read/set all fields
from a planar scan as well. There is therefore no need anymore for
the boolean.
Output of a conversion of planar data is now directly readable into STIR.
---
documentation/release_6.2.htm | 15 +++
src/utilities/SPECT_dicom_to_interfile.cxx | 134 ++++++++++-----------
2 files changed, 79 insertions(+), 70 deletions(-)
mode change 100644 => 100755 src/utilities/SPECT_dicom_to_interfile.cxx
diff --git a/documentation/release_6.2.htm b/documentation/release_6.2.htm
index 6bc5032970..5fdf4cd75f 100644
--- a/documentation/release_6.2.htm
+++ b/documentation/release_6.2.htm
@@ -91,6 +91,21 @@ <h3>Changed functionality</h3>
<br>
<a href=https://github.com/UCL/STIR/pull/1464>PR #1464</a>
</li>
+ <li>
+ <tt>SPECT_dicom_to_interfile</tt> improvements:
+ <ul>
+ <li>remove requirement for the <tt>is_planar</tt> parameters.
+ As STIR can only read SPECT sinograms, we now read/set all fields
+ from a planar scan as well. There is therefore no need anymore for
+ the boolean, and it is just ignored.
+ Output of a conversion of planar data is now directly readable into STIR.
+ </li>
+ <li>
+ do checks if sequences are present to avoid seg-faults
+ </li>
+ </ul>
+ See <a href=https://github.com/UCL/STIR/pull/1473>PR #1473</a>
+ </li>
</ul>
diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx
old mode 100644
new mode 100755
index ad08f258c2..3a65bbc7cc
--- a/src/utilities/SPECT_dicom_to_interfile.cxx
+++ b/src/utilities/SPECT_dicom_to_interfile.cxx
@@ -60,8 +60,7 @@ class SPECTDICOMData
SPECTDICOMData(const std::string& DICOM_filename) { dicom_filename = DICOM_filename; };
stir::Succeeded get_interfile_header(std::string& output_header, const std::string& data_filename, const int dataset_num) const;
stir::Succeeded get_proj_data(const std::string& output_file) const;
- stir::Succeeded open_dicom_file(bool is_planar);
- bool is_planar;
+ stir::Succeeded open_dicom_file();
int num_energy_windows = 1;
int num_frames = 0;
int num_of_projections = 0;
@@ -275,7 +274,7 @@ GetRadionuclideInfo(const gdcm::File& file, const RadionuclideInfo request, std:
}
stir::Succeeded
-SPECTDICOMData::open_dicom_file(bool is_planar)
+SPECTDICOMData::open_dicom_file()
{
stir::info(boost::format("SPECTDICOMData: opening file %1%") % dicom_filename);
@@ -380,54 +379,56 @@ SPECTDICOMData::open_dicom_file(bool is_planar)
}
}
- if (!is_planar)
+ if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x0053), no_of_proj_as_str) == stir::Succeeded::yes)
{
- if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x0053), no_of_proj_as_str) == stir::Succeeded::yes)
- {
- num_of_projections = std::stoi(no_of_proj_as_str);
- std::cout << "Number of projections: " << num_of_projections << std::endl;
- }
+ num_of_projections = std::stoi(no_of_proj_as_str);
+ std::cout << "Number of projections: " << num_of_projections << std::endl;
+ }
- if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1140), direction_of_rotation) == stir::Succeeded::yes)
- {
- if (direction_of_rotation == "CC")
- direction_of_rotation = "CCW";
+ if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1140), direction_of_rotation) == stir::Succeeded::yes)
+ {
+ if (direction_of_rotation == "CC")
+ direction_of_rotation = "CCW";
- std::cout << "Direction of rotation: " << direction_of_rotation << std::endl;
- }
+ std::cout << "Direction of rotation: " << direction_of_rotation << std::endl;
+ }
- if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x0200), start_angle_as_string) == stir::Succeeded::yes)
- {
- start_angle = std::stof(start_angle_as_string);
- std::cout << "Starting angle: " << std::fixed << std::setprecision(6) << start_angle << std::endl;
- }
+ if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x0200), start_angle_as_string) == stir::Succeeded::yes)
+ {
+ start_angle = std::stof(start_angle_as_string);
+ std::cout << "Starting angle: " << std::fixed << std::setprecision(6) << start_angle << std::endl;
+ }
- if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x1322), calib_factor_as_string) == stir::Succeeded::yes)
- {
- calibration_factor = std::stof(calib_factor_as_string);
- std::cout << "calibration factor: " << std::fixed << std::setprecision(6) << calibration_factor << std::endl;
- }
+ if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x1322), calib_factor_as_string) == stir::Succeeded::yes)
+ {
+ calibration_factor = std::stof(calib_factor_as_string);
+ std::cout << "calibration factor: " << std::fixed << std::setprecision(6) << calibration_factor << std::endl;
+ }
- if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1144), angular_step_as_string) == stir::Succeeded::yes)
- {
- angular_step = std::stof(angular_step_as_string);
- std::cout << "Angular step: " << std::fixed << std::setprecision(6) << angular_step << std::endl;
- }
+ if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1144), angular_step_as_string) == stir::Succeeded::yes)
+ {
+ angular_step = std::stof(angular_step_as_string);
+ std::cout << "Angular step: " << std::fixed << std::setprecision(6) << angular_step << std::endl;
+ }
- if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1143), extent_of_rotation_as_string) == stir::Succeeded::yes)
- {
- extent_of_rotation = std::stoi(extent_of_rotation_as_string);
- std::cout << "Rotation extent: " << extent_of_rotation << std::endl;
- }
+ if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1143), extent_of_rotation_as_string) == stir::Succeeded::yes)
+ {
+ extent_of_rotation = std::stoi(extent_of_rotation_as_string);
+ std::cout << "Rotation extent: " << extent_of_rotation << std::endl;
+ }
+ else
+ {
+ extent_of_rotation = 360;
+ stir::warning("Rotation was not present in DICOM file. Setting it to 360");
+ }
- if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1142), radius_as_string) == stir::Succeeded::yes)
- {
- rotation_radius = (radius_as_string);
- char slash = '\\';
- char comma = ',';
- std::cout << "Radius: " << radius_as_string << " " << slash << std::endl;
- std::replace(rotation_radius.begin(), rotation_radius.end(), slash, comma);
- }
+ if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1142), radius_as_string) == stir::Succeeded::yes)
+ {
+ rotation_radius = (radius_as_string);
+ char slash = '\\';
+ char comma = ',';
+ std::cout << "Radius: " << radius_as_string << " " << slash << std::endl;
+ std::replace(rotation_radius.begin(), rotation_radius.end(), slash, comma);
}
num_dimensions = 2;
@@ -491,6 +492,7 @@ SPECTDICOMData::get_interfile_header(std::string& output_header, const std::stri
ss << std::endl;
ss << "!GENERAL IMAGE DATA :=" << std::endl;
+ // We will write planar data as SPECT (i.e. "tomographic") anyway
ss << "!type of data := Tomographic" << std::endl;
ss << "imagedata byte order := LITTLEENDIAN" << std::endl;
ss << "!number format := float" << std::endl;
@@ -515,35 +517,28 @@ SPECTDICOMData::get_interfile_header(std::string& output_header, const std::stri
ss << "!number of projections := " << this->num_of_projections << std::endl;
ss << "number of time frames := 1" << std::endl;
ss << "!image duration (sec)[1] := " << this->num_of_projections * this->actual_frame_duration / 1000 << std::endl;
- if (!is_planar)
- ss << "!extent of rotation := " << this->extent_of_rotation << std::endl;
+ ss << "!extent of rotation := " << this->extent_of_rotation << std::endl;
ss << "!process status := acquired" << std::endl;
ss << std::endl;
ss << "!SPECT STUDY (acquired data) :=" << std::endl;
- if (!is_planar)
+ if (!this->direction_of_rotation.empty())
ss << "!direction of rotation := " << this->direction_of_rotation << std::endl;
ss << "start angle := " << this->start_angle << std::endl;
- if (is_planar)
- {
- ss << "orbit := circular" << std::endl;
- ss << "radius := " << 0 << std::endl;
- }
- else
- {
- if (this->rotation_radius.find(",") != std::string::npos)
- {
- ss << "orbit := non-circular" << std::endl;
- ss << "radii := {" << this->rotation_radius << "}" << std::endl;
- std::cout << "orbit := non-circular" << std::endl;
- }
- else
- {
- ss << "orbit := circular" << std::endl;
- ss << "radius := " << this->rotation_radius << "" << std::endl;
- }
- }
+ {
+ if (this->rotation_radius.find(",") != std::string::npos)
+ {
+ ss << "orbit := non-circular" << std::endl;
+ ss << "radii := {" << this->rotation_radius << "}" << std::endl;
+ std::cout << "orbit := non-circular" << std::endl;
+ }
+ else
+ {
+ ss << "orbit := circular" << std::endl;
+ ss << "radius := " << this->rotation_radius << "" << std::endl;
+ }
+ }
ss << std::endl;
ss << "!END OF INTERFILE :=";
@@ -638,20 +633,19 @@ int
main(int argc, char* argv[])
{
- if (argc != 4)
+ if (argc != 3 && argc != 4)
{
- std::cerr << "Usage: " << argv[0] << " <output_interfile_prefix> sinogram(dcm)> is_planar?\n";
+ std::cerr << "Usage: " << argv[0] << " <output_interfile_prefix> sinogram(dcm)> [is_planar?]\n"
+ << "Note: the is_planar value is ignored. All data will be written as SPECT.\n";
exit(EXIT_FAILURE);
}
SPECTDICOMData spect(argv[2]);
const std::string output_prefix(argv[1]);
- spect.is_planar = atoi(argv[3]);
-
try
{
- if (spect.open_dicom_file(spect.is_planar) == stir::Succeeded::no)
+ if (spect.open_dicom_file() == stir::Succeeded::no)
{
std::cerr << "Failed to read!" << std::endl;
return EXIT_FAILURE;