From bccf18fcf793fa5c9b0714c65df33e8990dd6d1b Mon Sep 17 00:00:00 2001 From: spencermanwell Date: Tue, 9 Apr 2024 11:36:52 -0400 Subject: [PATCH 1/6] Updates to the SPECT_dicom_to_interfile utility: - 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. --- src/utilities/SPECT_dicom_to_interfile.cxx | 311 +++++++++++++++------ 1 file changed, 225 insertions(+), 86 deletions(-) diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx index 3a65bbc7cc..e7353a02f1 100755 --- a/src/utilities/SPECT_dicom_to_interfile.cxx +++ b/src/utilities/SPECT_dicom_to_interfile.cxx @@ -45,6 +45,10 @@ enum class RadionuclideInfo CodingSchemDesignator, CodeMeaning }; +enum class DetectorInfo +{ + RadialPosition +}; enum class CalibrationInfo { }; @@ -53,6 +57,8 @@ stir::Succeeded GetEnergyWindowInfo(const gdcm::File& file, const EnergyWindowInfo request, std::string& dst, const int sequence_idx); stir::Succeeded GetRadionuclideInfo(const gdcm::File& file, const RadionuclideInfo request, std::string& dst, const int sequence_idx = 1); +stir::Succeeded +GetDetectorInfo(const gdcm::File& file, const DetectorInfo request, std::string& dst, const int sequence_idx = -1); class SPECTDICOMData { @@ -274,7 +280,74 @@ GetRadionuclideInfo(const gdcm::File& file, const RadionuclideInfo request, std: } stir::Succeeded -SPECTDICOMData::open_dicom_file() +GetDetectorInfo( const gdcm::File& file, const DetectorInfo request, std::string& dst, const int sequence_idx ) +{ + // Check for an unrecognized request. + if (request != DetectorInfo::RadialPosition) + { + dst = ""; + return stir::Succeeded::no; + } + + try + { + const gdcm::DataElement& de = file.GetDataSet().GetDataElement(gdcm::Tag(0x0054, 0x0022)); + const gdcm::SequenceOfItems* sqi = de.GetValueAsSQ(); + const int num_detector_information_squence_items = sqi->GetNumberOfItems(); + + for (int _sequence_id = 1; _sequence_id <= num_detector_information_squence_items; _sequence_id++) + { + // If a specific sequence item was requested, ignore all others. + if (sequence_idx != -1 && _sequence_id != sequence_idx) + continue; + + if (request == DetectorInfo::RadialPosition) + { + std::string radius_as_string; + if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1142), radius_as_string, _sequence_id) == stir::Succeeded::yes) + { + // Replace the \\ characters with commas. + std::replace(radius_as_string.begin(), radius_as_string.end(), char('\\'), char(',')); + // Remove empty spaces. + radius_as_string.erase(std::remove_if(radius_as_string.begin(), radius_as_string.end(), isspace), + radius_as_string.end()); + + + // For all but the first string copy we need to prepend the radius_as_string with "," + if (dst.empty()) + { + dst = (radius_as_string); + } + else + { + dst += ","; + dst += (radius_as_string); + } + std::cout << "Radii (Detector " << _sequence_id << "): " << radius_as_string << std::endl; + } + else + { + stir::warning("GetDetectorInfo: cannot read detector info"); + dst = ""; + return stir::Succeeded::no; + } + } + } + return stir::Succeeded::yes; + } + catch (...) + { + stir::warning("GetDetectorInfo: cannot read detector info"); + dst = ""; + return stir::Succeeded::no; + } + + dst = ""; + return stir::Succeeded::no; +} + +stir::Succeeded +SPECTDICOMData::open_dicom_file(bool is_planar) { stir::info(boost::format("SPECTDICOMData: opening file %1%") % dicom_filename); @@ -306,6 +379,7 @@ SPECTDICOMData::open_dicom_file() std::cout << "Patient name: " << patient_name << std::endl; } std::string no_of_proj_as_str; + std::string no_of_det_as_str; std::string start_angle_as_string; std::string angular_step_as_string; std::string extent_of_rotation_as_string; @@ -328,107 +402,172 @@ SPECTDICOMData::open_dicom_file() } this->num_of_projections = 1; - this->calibration_factor = -1; - if (GetRadionuclideInfo(file, RadionuclideInfo::CodeMeaning, isotope_name) == stir::Succeeded::yes) + if (!is_planar) { - std::cout << "Isotope name: " << isotope_name << std::endl; - } + // Determine the number of projections + // Given as the product of the number of frames in rotation and the number of detectors. + int num_of_frames_per_detector = 0; + if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x0053), no_of_proj_as_str) == stir::Succeeded::yes) + { + num_of_frames_per_detector = std::stoi(no_of_proj_as_str); + std::cout << "Number of frames per detector: " << num_of_frames_per_detector << 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); - } + int num_of_detectors = 0; + if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x0021), no_of_det_as_str) == stir::Succeeded::yes) + { + num_of_detectors = std::stoi(no_of_det_as_str); + } - { - 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 ((num_of_frames_per_detector * num_of_detectors) > 0) { - 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); + num_of_projections = num_of_frames_per_detector * num_of_detectors; + std::cout << "Number of projections: " << num_of_projections << std::endl; + } - 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 (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1140), direction_of_rotation) == stir::Succeeded::yes) + { - if (GetEnergyWindowInfo(file, EnergyWindowInfo::UpperThreshold, upper_window_as_string, w) == stir::Succeeded::yes) + if (direction_of_rotation == "CC") + direction_of_rotation = "CCW"; + + 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, 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; + } + 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) + { + 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; + } + + // Get the radial positions of the detector head(s) for all heads. + if (GetDetectorInfo(file, DetectorInfo::RadialPosition, radius_as_string) == stir::Succeeded::yes) + { + rotation_radius = radius_as_string; + std::cout << "Radial Positions: " << rotation_radius << 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) { - 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 (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 (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; - } + 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; + } - 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; - } - else - { - extent_of_rotation = 360; - stir::warning("Rotation was not present in DICOM file. Setting it to 360"); - } + 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; @@ -496,7 +635,7 @@ SPECTDICOMData::get_interfile_header(std::string& output_header, const std::stri ss << "!type of data := Tomographic" << std::endl; ss << "imagedata byte order := LITTLEENDIAN" << std::endl; ss << "!number format := float" << std::endl; - ss << "!number of bytes per pixel := 4" << std::endl; + ss << "!number of bytes per pixel := " << sizeof(float) << std::endl; ss << "calibration factor:= " << this->calibration_factor << std::endl; if (!this->isotope_name.empty()) ss << "isotope name:= " << this->isotope_name << std::endl; From 9ccd0e3248abffc78ce5b0b52b05041860a2db42 Mon Sep 17 00:00:00 2001 From: spencermanwell Date: Thu, 11 Apr 2024 10:32:49 -0400 Subject: [PATCH 2/6] When converting SPECT projection data from DICOM to interfile format, 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 63064ff608e5c3a8cc4e52ceb5355634b871f219) --- src/utilities/SPECT_dicom_to_interfile.cxx | 29 +++++++++++++++++----- 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx index e7353a02f1..36ec49cc77 100755 --- a/src/utilities/SPECT_dicom_to_interfile.cxx +++ b/src/utilities/SPECT_dicom_to_interfile.cxx @@ -32,6 +32,7 @@ #include "stir/error.h" #include "stir/warning.h" #include "stir/Succeeded.h" +#include "stir/IO/Interfileheader.h" enum class EnergyWindowInfo { @@ -81,7 +82,7 @@ class SPECTDICOMData int actual_frame_duration = 0; // frame duration in msec int num_of_rotations = 0; std::string direction_of_rotation, isotope_name; - int extent_of_rotation; + float extent_of_rotation; float calibration_factor; std::string rotation_radius; @@ -382,7 +383,7 @@ SPECTDICOMData::open_dicom_file(bool is_planar) std::string no_of_det_as_str; 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; std::string radius_as_string; std::string actual_frame_duration_as_string; std::string calib_factor_as_string; @@ -459,11 +460,27 @@ SPECTDICOMData::open_dicom_file(bool is_planar) 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) + // 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 + // do NOT overlap. + float scan_arc = 0; + if (GetDICOMTagInfo(file, gdcm::Tag(0x0018, 0x1143), scan_arc_as_string) == stir::Succeeded::yes) { - extent_of_rotation = std::stoi(extent_of_rotation_as_string); - std::cout << "Rotation extent: " << extent_of_rotation << std::endl; + scan_arc = std::stof(scan_arc_as_string); + std::cout << "Scan Arc: " << std::fixed << std::setprecision(6) << scan_arc << std::endl; + } + + // 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; + } + else + { + stir::warning(boost::format("SPECTDICOMData: cannot determine extent of rotation.")); + extent_of_rotation = stir::MinimalInterfileHeader::double_value_not_set; } + std::cout << "Rotation extent: " << std::fixed << std::setprecision(6) << extent_of_rotation << std::endl; // Get the radial positions of the detector head(s) for all heads. if (GetDetectorInfo(file, DetectorInfo::RadialPosition, radius_as_string) == stir::Succeeded::yes) @@ -626,7 +643,7 @@ SPECTDICOMData::get_interfile_header(std::string& output_header, const std::stri ss << "!version of keys := 3.3" << std::endl; ss << "name of data file := " << data_filename_only << std::endl; ss << "data offset in bytes := " - << this->matrix_size.at(0) * this->matrix_size.at(1) * this->num_of_projections * dataset_num * 4 // float hard-wired + << this->matrix_size.at(0) * this->matrix_size.at(1) * this->num_of_projections * dataset_num * sizeof(float) << std::endl; ss << std::endl; From 862adaf419755c100569ae8dd88cd24b0b0d52ec Mon Sep 17 00:00:00 2001 From: spencermanwell Date: Tue, 4 Jun 2024 10:06:12 -0400 Subject: [PATCH 3/6] SPECT_dicom_to_interfile utility failed in the presence of compressed 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. --- src/utilities/SPECT_dicom_to_interfile.cxx | 70 +++++++++++++++------- 1 file changed, 50 insertions(+), 20 deletions(-) diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx index 36ec49cc77..e11ce52ee8 100755 --- a/src/utilities/SPECT_dicom_to_interfile.cxx +++ b/src/utilities/SPECT_dicom_to_interfile.cxx @@ -22,11 +22,18 @@ */ #include #include +#include +#include #include #include #include +#include +#include +#include +#include +#include #include "stir/info.h" #include "stir/error.h" @@ -66,8 +73,10 @@ class SPECTDICOMData public: 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(); + stir::Succeeded get_proj_data(const std::string& output_file, const std::string& input_file = "" ) const; + stir::Succeeded open_dicom_file(bool is_planar); + void set_dicom_filename(const std::string& DICOM_filename) { dicom_filename = DICOM_filename; }; + bool is_planar; int num_energy_windows = 1; int num_frames = 0; int num_of_projections = 0; @@ -705,12 +714,19 @@ SPECTDICOMData::get_interfile_header(std::string& output_header, const std::stri } stir::Succeeded -SPECTDICOMData::get_proj_data(const std::string& output_file) const +SPECTDICOMData::get_proj_data(const std::string& output_file, const std::string& input_file ) const { - std::unique_ptr DICOM_reader(new gdcm::Reader); - DICOM_reader->SetFileName(dicom_filename.c_str()); - + std::unique_ptr DICOM_reader(new gdcm::ImageReader); + if (input_file.empty()) + { + DICOM_reader->SetFileName(dicom_filename.c_str()); + } + else + { + DICOM_reader->SetFileName(input_file.c_str()); + } + try { if (!DICOM_reader->Read()) @@ -726,24 +742,38 @@ SPECTDICOMData::get_proj_data(const std::string& output_file) const } const gdcm::File& file = DICOM_reader->GetFile(); - const gdcm::DataElement& de = file.GetDataSet().GetDataElement(gdcm::Tag(0x7fe0, 0x0010)); - const gdcm::ByteValue* bv = de.GetByteValue(); + const gdcm::ByteValue* bv = de.GetByteValue(); // Only works when data are not compressed. + + // Check if the pixel data pointer is NULL, if so the data may be compressed. 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); + stir::warning(boost::format("SPECTDICOMData: decompressing pixel data from: %1%") % dicom_filename); + gdcm::ImageChangeTransferSyntax change; + change.SetTransferSyntax(gdcm::TransferSyntax::ImplicitVRLittleEndian); + change.SetInput(DICOM_reader->GetImage()); - if (!outfile.is_open()) { - std::cerr << "Unable to write proj data to " << tmpFile; - return stir::Succeeded::no; - } - bv->WriteBuffer(outfile); - outfile.close();*/ + // Assert that the change in transfer syntax was successful. Return otherwise. + if (!change.Change()) + { + stir::error(boost::format("SPECTDICOMData: decompression failed.") % dicom_filename); + return stir::Succeeded::no; + } + + // Write a new file using the uncompressed transfer syntax. + const std::string uncompressed_dicom_filename = dicom_filename + std::string(".uncmp"); + gdcm::ImageWriter writer; + writer.SetImage(change.GetOutput()); + writer.SetFile(DICOM_reader->GetFile()); + writer.SetFileName(uncompressed_dicom_filename.c_str()); + if (!writer.Write()) + { + return stir::Succeeded::no; + } + + return get_proj_data(output_file, uncompressed_dicom_filename); + + } uint64_t len0 = (uint64_t)bv->GetLength() / 2; std::cout << "Number of data points = " << len0 << std::endl; From 6c449fe7866a70207ad701e6c7dd429d451ecb30 Mon Sep 17 00:00:00 2001 From: spencermanwell Date: Mon, 24 Mar 2025 11:53:06 -0400 Subject: [PATCH 4/6] For SPECT_dicom_to_interfile utility, replace const pointers with cosnt smart pointers for calls to gdcm::DataElement::GertValueAsSQ() --- src/utilities/SPECT_dicom_to_interfile.cxx | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx index e11ce52ee8..199b974379 100755 --- a/src/utilities/SPECT_dicom_to_interfile.cxx +++ b/src/utilities/SPECT_dicom_to_interfile.cxx @@ -143,7 +143,7 @@ GetDICOMTagInfo(const gdcm::File& file, const gdcm::Tag tag, std::string& dst, c try { const gdcm::DataElement& de = file.GetDataSet().GetDataElement(t); - const gdcm::SequenceOfItems* sqi = de.GetValueAsSQ(); + const gdcm::SmartPointer sqi = de.GetValueAsSQ(); if (!sqi) { // try next sequence @@ -156,6 +156,7 @@ GetDICOMTagInfo(const gdcm::File& file, const gdcm::Tag tag, std::string& dst, c // try next sequence continue; } + const gdcm::Item& item = sqi->GetItem(sequence_idx); element = item.GetDataElement(tag); @@ -200,12 +201,12 @@ GetEnergyWindowInfo(const gdcm::File& file, const EnergyWindowInfo request, std: // 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 sqi = de.GetValueAsSQ(); const gdcm::Item& item = sqi->GetItem(sequence_idx); // Get Energy Window Range Sequence const gdcm::DataElement& element = item.GetDataElement(energy_window_range_seq); - const gdcm::SequenceOfItems* sqi2 = element.GetValueAsSQ(); + const gdcm::SmartPointer sqi2 = element.GetValueAsSQ(); if (sqi2->GetNumberOfItems() > 1) stir::warning("Energy window sequence contains more than 1 window. Ignoring all later ones"); const gdcm::Item& item2 = sqi2->GetItem(1); @@ -254,7 +255,7 @@ GetRadionuclideInfo(const gdcm::File& file, const RadionuclideInfo request, std: // Get Radiopharmaceutical Info Sequence const gdcm::DataElement& de = file.GetDataSet().GetDataElement(radiopharm_info_seq_tag); - const gdcm::SequenceOfItems* sqi = de.GetValueAsSQ(); + const gdcm::SmartPointer sqi = de.GetValueAsSQ(); const gdcm::Item& item = sqi->GetItem(sequence_idx); // Get Radiopnuclide Code Sequence @@ -262,7 +263,7 @@ GetRadionuclideInfo(const gdcm::File& file, const RadionuclideInfo request, std: const gdcm::DataElement& element = item.GetDataElement(radionuclide_code_seq_tag); if (element.GetVL() > 0) { - const gdcm::SequenceOfItems* sqi2 = element.GetValueAsSQ(); + const gdcm::SmartPointer sqi2 = element.GetValueAsSQ(); const gdcm::Item& item2 = sqi2->GetItem(sequence_idx); // std::cout<< "num items"<< sqi2->GetNumberOfItems(); // std::cout << item2 << std::endl; @@ -302,7 +303,7 @@ GetDetectorInfo( const gdcm::File& file, const DetectorInfo request, std::string try { const gdcm::DataElement& de = file.GetDataSet().GetDataElement(gdcm::Tag(0x0054, 0x0022)); - const gdcm::SequenceOfItems* sqi = de.GetValueAsSQ(); + const gdcm::SmartPointer sqi = de.GetValueAsSQ(); const int num_detector_information_squence_items = sqi->GetNumberOfItems(); for (int _sequence_id = 1; _sequence_id <= num_detector_information_squence_items; _sequence_id++) From 9c3536b83d03e3b98eec50013af2751456a1b6f0 Mon Sep 17 00:00:00 2001 From: spencermanwell Date: Mon, 24 Mar 2025 14:23:02 -0400 Subject: [PATCH 5/6] Clean up SPECT_dicom_to_interfile.cxx following rebase of this branch 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. --- src/utilities/SPECT_dicom_to_interfile.cxx | 169 ++++++++------------- 1 file changed, 62 insertions(+), 107 deletions(-) diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx index 199b974379..eb15fc89ff 100755 --- a/src/utilities/SPECT_dicom_to_interfile.cxx +++ b/src/utilities/SPECT_dicom_to_interfile.cxx @@ -92,7 +92,7 @@ class SPECTDICOMData int num_of_rotations = 0; std::string direction_of_rotation, isotope_name; float extent_of_rotation; - float calibration_factor; + float calibration_factor = -1; std::string rotation_radius; std::vector lower_en_window_thres; @@ -403,14 +403,64 @@ SPECTDICOMData::open_dicom_file(bool is_planar) std::string lower_window_as_string; std::string upper_window_as_string; - { - std::string str; - if (GetDICOMTagInfo(file, gdcm::Tag(0x0028, 0x0008), str) == stir::Succeeded::yes) - { - num_frames = std::stoi(str); - std::cout << "Number of frames: " << num_frames << std::endl; - } - } + + std::string str; + if (GetDICOMTagInfo(file, gdcm::Tag(0x0028, 0x0008), str) == stir::Succeeded::yes) + { + num_frames = std::stoi(str); + std::cout << "Number of frames: " << num_frames << std::endl; + } + + 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; + } + } + } this->num_of_projections = 1; if (!is_planar) @@ -438,7 +488,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"; @@ -451,19 +500,6 @@ SPECTDICOMData::open_dicom_file(bool is_planar) 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; - } - 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) { angular_step = std::stof(angular_step_as_string); @@ -499,83 +535,6 @@ SPECTDICOMData::open_dicom_file(bool is_planar) std::cout << "Radial Positions: " << rotation_radius << 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 (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; - } - - 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; - } - - 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(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); @@ -587,13 +546,10 @@ SPECTDICOMData::open_dicom_file(bool is_planar) 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) + if (GetDICOMTagInfo(file, gdcm::Tag(0x0054, 0x1322), calib_factor_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); + calibration_factor = std::stof(calib_factor_as_string); + std::cout << "calibration factor: " << std::fixed << std::setprecision(6) << calibration_factor << std::endl; } } @@ -623,7 +579,6 @@ SPECTDICOMData::open_dicom_file(bool is_planar) { curr = pixel_size_as_string.find('\\', prev); std::string found = pixel_size_as_string.substr(prev, curr - prev); - // std::cout << "found = " << found << std::endl; pixel_sizes.push_back(std::stof(found)); prev = curr + 1; } From e971fbae2e276483020e8e5d001e87cd6b4765d2 Mon Sep 17 00:00:00 2001 From: spencermanwell Date: Wed, 26 Mar 2025 09:45:20 -0400 Subject: [PATCH 6/6] Additional cleanup activities for SPECT_dicom_to_interfile utility: - 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. --- src/utilities/SPECT_dicom_to_interfile.cxx | 110 ++++++++++++--------- 1 file changed, 61 insertions(+), 49 deletions(-) diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx index eb15fc89ff..c0eb3a9d6e 100755 --- a/src/utilities/SPECT_dicom_to_interfile.cxx +++ b/src/utilities/SPECT_dicom_to_interfile.cxx @@ -73,17 +73,16 @@ class SPECTDICOMData public: 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 std::string& input_file = "" ) const; + stir::Succeeded get_proj_data(const std::string& output_file ) const; stir::Succeeded open_dicom_file(bool is_planar); - void set_dicom_filename(const std::string& DICOM_filename) { dicom_filename = DICOM_filename; }; - bool is_planar; + std::vector store_pixel_buffer_as_float_vector(const gdcm::ByteValue*& bv) const; + bool is_planar = false; int num_energy_windows = 1; int num_frames = 0; int num_of_projections = 0; std::vector energy_window_name; private: - // shared_ptr proj_data_info_sptr; std::string dicom_filename; float start_angle = 0.0f; @@ -91,14 +90,13 @@ class SPECTDICOMData int actual_frame_duration = 0; // frame duration in msec int num_of_rotations = 0; std::string direction_of_rotation, isotope_name; - float extent_of_rotation; + float extent_of_rotation = 0.0f; float calibration_factor = -1; std::string rotation_radius; std::vector lower_en_window_thres; std::vector upper_en_window_thres; - int num_dimensions; std::vector matrix_labels; std::vector matrix_size; std::vector pixel_sizes; @@ -553,8 +551,6 @@ SPECTDICOMData::open_dicom_file(bool is_planar) } } - num_dimensions = 2; - matrix_labels.push_back("axial coordinate"); matrix_labels.push_back("bin coordinate"); @@ -669,26 +665,36 @@ SPECTDICOMData::get_interfile_header(std::string& output_header, const std::stri return stir::Succeeded::yes; } -stir::Succeeded -SPECTDICOMData::get_proj_data(const std::string& output_file, const std::string& input_file ) const +std::vector +SPECTDICOMData::store_pixel_buffer_as_float_vector(const gdcm::ByteValue*& bv) const { + std::vector pixel_data_as_float; + uint64_t n_elements = (uint64_t)bv->GetLength() / sizeof(uint16_t); + uint16_t* ptr = (uint16_t*)bv->GetPointer(); - std::unique_ptr DICOM_reader(new gdcm::ImageReader); - if (input_file.empty()) - { - DICOM_reader->SetFileName(dicom_filename.c_str()); - } - else + uint64_t ct = 0; + while (ct < n_elements) { - DICOM_reader->SetFileName(input_file.c_str()); + uint16_t val = *ptr; + pixel_data_as_float.push_back((float)val); + ptr++; + ct++; } - + + return pixel_data_as_float; +} + +stir::Succeeded +SPECTDICOMData::get_proj_data(const std::string& output_file) const +{ + + std::unique_ptr DICOM_reader(new gdcm::ImageReader); + DICOM_reader->SetFileName(dicom_filename.c_str()); try { if (!DICOM_reader->Read()) { - 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); } } catch (const std::string& e) @@ -700,51 +706,57 @@ SPECTDICOMData::get_proj_data(const std::string& output_file, const std::string& const gdcm::File& file = DICOM_reader->GetFile(); const gdcm::DataElement& de = file.GetDataSet().GetDataElement(gdcm::Tag(0x7fe0, 0x0010)); const gdcm::ByteValue* bv = de.GetByteValue(); // Only works when data are not compressed. - - // Check if the pixel data pointer is NULL, if so the data may be compressed. + + // Check if the pixel data pointer is NULL, if so the data were likely compressed. + std::vector pixel_data_as_float; if (!bv) { - stir::warning(boost::format("SPECTDICOMData: decompressing pixel data from: %1%") % dicom_filename); + stir::warning(boost::format("in SPECTDICOMData::get_proj_data() - attempting to decompress pixel data from: %1%") + % dicom_filename); gdcm::ImageChangeTransferSyntax change; change.SetTransferSyntax(gdcm::TransferSyntax::ImplicitVRLittleEndian); change.SetInput(DICOM_reader->GetImage()); // Assert that the change in transfer syntax was successful. Return otherwise. - if (!change.Change()) + try { - stir::error(boost::format("SPECTDICOMData: decompression failed.") % dicom_filename); + if (!change.Change()) + { + stir::error(boost::format("in SPECTDICOMData::get_proj_data() - decompression failed.")); + } + } + catch (const std::string& e) + { + std::cerr << e << std::endl; return stir::Succeeded::no; } - // Write a new file using the uncompressed transfer syntax. - const std::string uncompressed_dicom_filename = dicom_filename + std::string(".uncmp"); - gdcm::ImageWriter writer; - writer.SetImage(change.GetOutput()); - writer.SetFile(DICOM_reader->GetFile()); - writer.SetFileName(uncompressed_dicom_filename.c_str()); - if (!writer.Write()) + const gdcm::DataElement& de_uncmp = change.GetOutput().GetDataElement(); + const gdcm::ByteValue* bv_uncmp = de_uncmp.GetByteValue(); + // Assert that pixel buffer is now accessible. Return otherwise. + try { + if (!bv_uncmp) + { + stir::error(boost::format("in SPECTDICOMData::get_proj_data() - pixel data are not accessible.")); + } + else + { + std::cout << "In SPECTDICOMData::get_proj_data() - decompression of pixel data successful." << std::endl; + pixel_data_as_float = store_pixel_buffer_as_float_vector(bv_uncmp); + } + } + catch (const std::string& e) + { + std::cerr << e << std::endl; return stir::Succeeded::no; } - - return get_proj_data(output_file, uncompressed_dicom_filename); - } - - uint64_t len0 = (uint64_t)bv->GetLength() / 2; - std::cout << "Number of data points = " << len0 << std::endl; - - std::vector pixel_data_as_float; - - uint16_t* ptr = (uint16_t*)bv->GetPointer(); - - uint64_t ct = 0; - while (ct < len0) + else { - uint16_t val = *ptr; - pixel_data_as_float.push_back((float)val); - ptr++; - ct++; + // No uncompression required. + std::cout << "In SPECTDICOMData::get_proj_data() - pixel data are accessible." << std::endl; + pixel_data_as_float = store_pixel_buffer_as_float_vector(bv); } std::ofstream final_out(output_file.c_str(), std::ios::out | std::ofstream::binary);