From 4bd5124fc4d84744d0d67283cc30aa3faa48fb91 Mon Sep 17 00:00:00 2001 From: nmdicom-recon Date: Wed, 4 Jan 2023 11:19:05 +0000 Subject: [PATCH 1/5] adding utility that extracts interfile+binary from dicom data --- src/utilities/CMakeLists.txt | 1 + src/utilities/SPECT_dicom_to_interfile.cxx | 505 +++++++++++++++++++++ 2 files changed, 506 insertions(+) create mode 100644 src/utilities/SPECT_dicom_to_interfile.cxx diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index 65fe3c290..4f4c1702c 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -73,6 +73,7 @@ set(${dir_EXE_SOURCES} write_sinogram_to_txt.cxx find_sum_projection_of_viewgram_and_sinogram.cxx separate_true_from_random_scatter_for_necr.cxx + SPECT_dicom_to_interfile.cxx ) if (AVW_FOUND) diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx new file mode 100644 index 000000000..73bb92111 --- /dev/null +++ b/src/utilities/SPECT_dicom_to_interfile.cxx @@ -0,0 +1,505 @@ +#include +#include + +#include +#include + +#include +#include + +#include "stir/info.h" +#include "stir/error.h" +#include "stir/Succeeded.h" + +enum class EnergyWindowInfo { LowerThreshold, UpperThreshold, WindowName }; +enum class RadionuclideInfo { CodeValue, CodingSchemDesignator, CodeMeaning }; +enum class CalibrationInfo {}; +stir::Succeeded GetDICOMTagInfo(const gdcm::File &file, const gdcm::Tag tag, std::string &dst); +stir::Succeeded GetEnergyWindowInfo(const gdcm::File &file, const EnergyWindowInfo request, std::string &dst); +stir::Succeeded GetRadionuclideInfo(const gdcm::File &file, const RadionuclideInfo request, std::string &dst); + +class SPECTDICOMData +{ +public: + SPECTDICOMData(const std::string& DICOM_filename){ dicom_filename = DICOM_filename; }; + stir::Succeeded get_interfile_header(std::string &output_header) const; + stir::Succeeded get_proj_data(const std::string &output_file) const; + stir::Succeeded open_dicom_file(bool is_planar); + void set_data_filename(const std::string &data_file){ data_filename = data_file; }; + bool is_planar; +private: + //shared_ptr proj_data_info_sptr; + std::string dicom_filename; + std::string data_filename; + + float start_angle = 0.0f; + int num_of_projections = 0; + float angular_step = 0.0f; + 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 calibration_factor; + std::string rotation_radius ; + + float lower_en_window_thres = 0.0f; + float upper_en_window_thres = 0.0f; + std::string energy_window_name; + + int num_dimensions; + std::vector matrix_labels; + std::vector matrix_size; + std::vector pixel_sizes; +}; + +stir::Succeeded GetDICOMTagInfo(const gdcm::File &file, const gdcm::Tag tag, std::string &dst){ + + //Extracts information for a given DICOM tag from a gdcm dataset. + //Tag contents are returned as a string in dst variable. + + //Tries to read the element associated with the tag. If the read fails, the + //DataElement should have a ByteValue of NULL. + + try { + const gdcm::DataSet &ds = file.GetDataSet(); + + gdcm::StringFilter sf; + sf.SetFile(file); + + std::stringstream inStream; + inStream.exceptions(std::ios::badbit); + + // First try and see if this is a standard tag. + gdcm::DataElement element = ds.GetDataElement(tag); + + if (element.GetByteValue() != NULL) { + dst = sf.ToString(tag); + return stir::Succeeded::yes; + } + + // Try: RotationInformationSequence (0054,0052) + // DetectorInformationSequence (0054,0022) + // EnergyWindowInformationSequence (0054,0012) + std::vector seqs = { gdcm::Tag(0x0054,0x0052), gdcm::Tag(0x0054,0x0022), gdcm::Tag(0x0054,0x0012)}; + + for (const auto& t : seqs) { + const gdcm::DataElement &de = file.GetDataSet().GetDataElement(t); + const gdcm::SequenceOfItems *sqi = de.GetValueAsSQ(); + const gdcm::Item &item = sqi->GetItem(1); + + element = item.GetDataElement(tag); + + if (element.GetByteValue() != NULL) { + dst = sf.ToString(element); + return stir::Succeeded::yes; + } + } + + } catch (std::bad_alloc){ + stir::error(boost::format("GetDICOMTagInfo: cannot read tag %1%") % tag); + return stir::Succeeded::no; + } + + return stir::Succeeded::no; +} + +stir::Succeeded GetEnergyWindowInfo(const gdcm::File &file, const EnergyWindowInfo request, std::string &dst){ + + if (request == EnergyWindowInfo::WindowName){ + return GetDICOMTagInfo(file, gdcm::Tag(0x0054,0x0018), dst); + } + + try { + const gdcm::Tag energy_window_info_seq = gdcm::Tag(0x0054,0x0012); + const gdcm::Tag energy_window_range_seq = gdcm::Tag(0x0054,0x0013); + + const gdcm::Tag lower_energy_window_tag = gdcm::Tag(0x0054,0x0014); + const gdcm::Tag upper_energy_window_tag = gdcm::Tag(0x0054,0x0015); + + //Get Energy Window Info Sequence + const gdcm::DataElement &de = file.GetDataSet().GetDataElement(energy_window_info_seq); + const gdcm::SequenceOfItems *sqi = de.GetValueAsSQ(); + const gdcm::Item &item = sqi->GetItem(1); + + //Get Energy Window Range Sequence + const gdcm::DataElement &element = item.GetDataElement(energy_window_range_seq); + const gdcm::SequenceOfItems *sqi2 = element.GetValueAsSQ(); + const gdcm::Item &item2 = sqi2->GetItem(1); + + //std::cout << item2 << std::endl; + + gdcm::DataElement window_element; + + if (request == EnergyWindowInfo::LowerThreshold) + window_element = item2.GetDataElement(lower_energy_window_tag); + else + window_element = item2.GetDataElement(upper_energy_window_tag); + + if (window_element.GetByteValue() != NULL) { + gdcm::StringFilter sf; + sf.SetFile(file); + dst = sf.ToString(window_element); + return stir::Succeeded::yes; + } + + } catch (std::bad_alloc){ + stir::error(boost::format("GetEnergyWindowInfo: cannot read energy info")); + return stir::Succeeded::no; + } + + return stir::Succeeded::no; +} + +stir::Succeeded GetRadionuclideInfo(const gdcm::File &file, const RadionuclideInfo request, std::string &dst){ + +// if (request == RadionuclideInfo::CodeMeaning){ +// return GetDICOMTagInfo(file, gdcm::Tag(0x0008,0x0018), dst); +// } + + try { + const gdcm::Tag radiopharm_info_seq_tag = gdcm::Tag(0x0054,0x0016); + const gdcm::Tag radionuclide_code_seq_tag = gdcm::Tag(0x0054,0x0300); + + const gdcm::Tag radionuclide_codemeaning_tag = gdcm::Tag(0x0008,0x0104); +// const gdcm::Tag upper_radionuclide_tag = gdcm::Tag(0x0054,0x0015); + + //Get Radiopharmaceutical Info Sequence + const gdcm::DataElement &de = file.GetDataSet().GetDataElement(radiopharm_info_seq_tag); + const gdcm::SequenceOfItems *sqi = de.GetValueAsSQ(); + const gdcm::Item &item = sqi->GetItem(1); + + + //Get Radiopnuclide Code Sequence + gdcm::DataElement nuclide_element; + const gdcm::DataElement &element = item.GetDataElement(radionuclide_code_seq_tag); + if(element.GetVL()>0){ + const gdcm::SequenceOfItems *sqi2 = element.GetValueAsSQ(); + const gdcm::Item &item2 = sqi2->GetItem(1); + //std::cout<< "num items"<< sqi2->GetNumberOfItems(); + //std::cout << item2 << std::endl; + +// gdcm::DataElement nuclide_element; + nuclide_element = item2.GetDataElement(radionuclide_codemeaning_tag); + } + + if (nuclide_element.GetByteValue() != NULL) { + gdcm::StringFilter sf; + sf.SetFile(file); + dst = sf.ToString(nuclide_element); + return stir::Succeeded::yes; + } + } catch (std::bad_alloc){ + stir::error(boost::format("GetRadionuclideInfo: cannot read radiopharaceutical info")); + return stir::Succeeded::no; + } + + return stir::Succeeded::no; +} + +stir::Succeeded SPECTDICOMData::open_dicom_file(bool is_planar) +{ + + stir::info(boost::format("SPECTDICOMData: opening file %1%") % dicom_filename); + + std::unique_ptr DICOM_reader(new gdcm::Reader); + 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; + } + } catch (const std::string &e){ + std::cerr << e << std::endl; + return stir::Succeeded::no; + } + + const gdcm::File &file = DICOM_reader->GetFile(); + + std::cout << std::endl; + + std::string patient_name; + if (GetDICOMTagInfo(file, gdcm::Tag(0x0010,0x0010), patient_name) == stir::Succeeded::yes){ + std::cout << "Patient name: " << patient_name << std::endl; + } + std::string no_of_proj_as_str; + std::string start_angle_as_string; + std::string angular_step_as_string; + std::string extent_of_rotation_as_string; + std::string radius_as_string; + std::string actual_frame_duration_as_string; + std::string calib_factor_as_string; + + std::string matrix_size_as_string; + std::string pixel_size_as_string; + std::string lower_window_as_string; + std::string upper_window_as_string; + + this->num_of_projections=1; + if(!is_planar){ + 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; + } + 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; + } + + 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 <<" " <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 :="; + + output_header = ss.str(); + + return stir::Succeeded::yes; +} + +stir::Succeeded SPECTDICOMData::get_proj_data(const std::string &output_file) const{ + + std::unique_ptr DICOM_reader(new gdcm::Reader); + 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; + } + } catch (const std::string &e){ + std::cerr << e << std::endl; + return stir::Succeeded::no; + } + + const gdcm::File &file = DICOM_reader->GetFile(); + + const gdcm::DataElement &de = file.GetDataSet().GetDataElement(gdcm::Tag(0x7fe0,0x0010)); + const gdcm::ByteValue *bv = de.GetByteValue(); + + /* + std::string tmpFile = "tmp.s"; + std::ofstream outfile(tmpFile.c_str(), std::ios::out | std::ios::binary); + + if (!outfile.is_open()) { + std::cerr << "Unable to write proj data to " << tmpFile; + return stir::Succeeded::no; + } + bv->WriteBuffer(outfile); + outfile.close();*/ + + uint64_t len0 = (uint64_t)bv->GetLength()/2; + std::cout << "Length = " << len0 << std::endl; + + std::vector pixel_data_as_float; + + uint16_t *ptr = (uint16_t*)bv->GetPointer(); + + uint64_t ct = 0; + while (ct < len0){ + uint16_t val = *ptr; + pixel_data_as_float.push_back((float)val); + ptr++; + ct++; + } + + std::cout << "pixel_data_as_float length = " << pixel_data_as_float.size() << std::endl; + + std::ofstream final_out(output_file.c_str(), std::ios::out | std::ofstream::binary); + final_out.write(reinterpret_cast(&pixel_data_as_float[0]), pixel_data_as_float.size() * sizeof(float)); + final_out.close(); + + return stir::Succeeded::yes; +} + +int main(int argc, char * argv[]) +{ + + if ( argc!=5) { + std::cerr << "Usage: " << argv[0] << " "<< " "<<""<<"is_planar?\n"; + exit(EXIT_FAILURE); + } + + std::unique_ptr + spect(new SPECTDICOMData(argv[1])); + + spect->is_planar=atoi(argv[4]); + + try{ + if (spect->open_dicom_file(spect->is_planar) == stir::Succeeded::no){ + std::cerr << "Failed to read!" << std::endl; + return EXIT_FAILURE; + } + } catch(const std::string& e){ + std::cerr << e << std::endl; + return EXIT_FAILURE; + } + + std::cout << "Read successfully!" << std::endl; + + spect->set_data_filename(argv[3]); + stir::Succeeded s = spect->get_proj_data(argv[3]); + + std::string header; + s = spect->get_interfile_header(header); + + std::cout << header << std::endl; + + std::filebuf fb; + fb.open (argv[2],std::ios::out); + std::ostream os(&fb); + os << header; + fb.close(); + + return EXIT_SUCCESS; +} From 49556608da933cca57b6fa5f3f0b0e0e69125d64 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Fri, 31 Mar 2023 14:10:51 +0100 Subject: [PATCH 2/5] remove some compiler warnings --- src/utilities/SPECT_dicom_to_interfile.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx index 73bb92111..ea93d47ae 100644 --- a/src/utilities/SPECT_dicom_to_interfile.cxx +++ b/src/utilities/SPECT_dicom_to_interfile.cxx @@ -142,7 +142,7 @@ stir::Succeeded GetEnergyWindowInfo(const gdcm::File &file, const EnergyWindowIn return stir::Succeeded::yes; } - } catch (std::bad_alloc){ + } catch (...){ stir::error(boost::format("GetEnergyWindowInfo: cannot read energy info")); return stir::Succeeded::no; } @@ -188,7 +188,7 @@ stir::Succeeded GetRadionuclideInfo(const gdcm::File &file, const RadionuclideIn dst = sf.ToString(nuclide_element); return stir::Succeeded::yes; } - } catch (std::bad_alloc){ + } catch (...){ stir::error(boost::format("GetRadionuclideInfo: cannot read radiopharaceutical info")); return stir::Succeeded::no; } @@ -279,7 +279,7 @@ stir::Succeeded SPECTDICOMData::open_dicom_file(bool is_planar) if (GetDICOMTagInfo(file, gdcm::Tag(0x0018,0x1142), radius_as_string) == stir::Succeeded::yes){ rotation_radius =(radius_as_string); char slash='\\'; - char comma='\,'; + char comma=','; std::cout << "Radius: " << radius_as_string <<" " < Date: Fri, 31 Mar 2023 19:15:41 +0100 Subject: [PATCH 3/5] handle multiple energy windows in SPECT DICOM --- src/utilities/SPECT_dicom_to_interfile.cxx | 153 ++++++++++++--------- 1 file changed, 91 insertions(+), 62 deletions(-) diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx index ea93d47ae..700694c25 100644 --- a/src/utilities/SPECT_dicom_to_interfile.cxx +++ b/src/utilities/SPECT_dicom_to_interfile.cxx @@ -9,28 +9,29 @@ #include "stir/info.h" #include "stir/error.h" +#include "stir/warning.h" #include "stir/Succeeded.h" enum class EnergyWindowInfo { LowerThreshold, UpperThreshold, WindowName }; enum class RadionuclideInfo { CodeValue, CodingSchemDesignator, CodeMeaning }; enum class CalibrationInfo {}; -stir::Succeeded GetDICOMTagInfo(const gdcm::File &file, const gdcm::Tag tag, std::string &dst); -stir::Succeeded GetEnergyWindowInfo(const gdcm::File &file, const EnergyWindowInfo request, std::string &dst); -stir::Succeeded GetRadionuclideInfo(const gdcm::File &file, const RadionuclideInfo request, std::string &dst); +stir::Succeeded GetDICOMTagInfo(const gdcm::File &file, const gdcm::Tag tag, std::string &dst, const int sequence_idx = 1); +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); class SPECTDICOMData { public: SPECTDICOMData(const std::string& DICOM_filename){ dicom_filename = DICOM_filename; }; - stir::Succeeded get_interfile_header(std::string &output_header) const; + 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); - void set_data_filename(const std::string &data_file){ data_filename = data_file; }; bool is_planar; + int num_energy_windows = 1; + std::vector energy_window_name; private: //shared_ptr proj_data_info_sptr; std::string dicom_filename; - std::string data_filename; float start_angle = 0.0f; int num_of_projections = 0; @@ -42,9 +43,8 @@ class SPECTDICOMData float calibration_factor; std::string rotation_radius ; - float lower_en_window_thres = 0.0f; - float upper_en_window_thres = 0.0f; - std::string energy_window_name; + std::vector lower_en_window_thres; + std::vector upper_en_window_thres; int num_dimensions; std::vector matrix_labels; @@ -52,7 +52,7 @@ class SPECTDICOMData std::vector pixel_sizes; }; -stir::Succeeded GetDICOMTagInfo(const gdcm::File &file, const gdcm::Tag tag, std::string &dst){ +stir::Succeeded GetDICOMTagInfo(const gdcm::File &file, const gdcm::Tag tag, std::string &dst, const int sequence_idx){ //Extracts information for a given DICOM tag from a gdcm dataset. //Tag contents are returned as a string in dst variable. @@ -83,30 +83,37 @@ stir::Succeeded GetDICOMTagInfo(const gdcm::File &file, const gdcm::Tag tag, std std::vector seqs = { gdcm::Tag(0x0054,0x0052), gdcm::Tag(0x0054,0x0022), gdcm::Tag(0x0054,0x0012)}; for (const auto& t : seqs) { - const gdcm::DataElement &de = file.GetDataSet().GetDataElement(t); - const gdcm::SequenceOfItems *sqi = de.GetValueAsSQ(); - const gdcm::Item &item = sqi->GetItem(1); - - element = item.GetDataElement(tag); - - if (element.GetByteValue() != NULL) { - dst = sf.ToString(element); - return stir::Succeeded::yes; - } + try + { + const gdcm::DataElement &de = file.GetDataSet().GetDataElement(t); + const gdcm::SequenceOfItems *sqi = de.GetValueAsSQ(); + const gdcm::Item &item = sqi->GetItem(sequence_idx); + + element = item.GetDataElement(tag); + + if (element.GetByteValue() != NULL) { + dst = sf.ToString(element); + return stir::Succeeded::yes; + } + } + catch (...) + { + // ignore error and try next sequence + } } - } catch (std::bad_alloc){ - stir::error(boost::format("GetDICOMTagInfo: cannot read tag %1%") % tag); + } catch (...){ + stir::error(boost::format("GetDICOMTagInfo: cannot read tag %1% sequence index %2%") % tag % sequence_idx); return stir::Succeeded::no; } return stir::Succeeded::no; } -stir::Succeeded GetEnergyWindowInfo(const gdcm::File &file, const EnergyWindowInfo request, std::string &dst){ +stir::Succeeded GetEnergyWindowInfo(const gdcm::File &file, const EnergyWindowInfo request, std::string &dst, const int sequence_idx){ if (request == EnergyWindowInfo::WindowName){ - return GetDICOMTagInfo(file, gdcm::Tag(0x0054,0x0018), dst); + return GetDICOMTagInfo(file, gdcm::Tag(0x0054,0x0018), dst, sequence_idx); } try { @@ -119,11 +126,13 @@ stir::Succeeded GetEnergyWindowInfo(const gdcm::File &file, const EnergyWindowIn //Get Energy Window Info Sequence const gdcm::DataElement &de = file.GetDataSet().GetDataElement(energy_window_info_seq); const gdcm::SequenceOfItems *sqi = de.GetValueAsSQ(); - const gdcm::Item &item = sqi->GetItem(1); + 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(); + 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); //std::cout << item2 << std::endl; @@ -150,7 +159,7 @@ stir::Succeeded GetEnergyWindowInfo(const gdcm::File &file, const EnergyWindowIn return stir::Succeeded::no; } -stir::Succeeded GetRadionuclideInfo(const gdcm::File &file, const RadionuclideInfo request, std::string &dst){ +stir::Succeeded GetRadionuclideInfo(const gdcm::File &file, const RadionuclideInfo request, std::string &dst, const int sequence_idx){ // if (request == RadionuclideInfo::CodeMeaning){ // return GetDICOMTagInfo(file, gdcm::Tag(0x0008,0x0018), dst); @@ -166,7 +175,7 @@ stir::Succeeded GetRadionuclideInfo(const gdcm::File &file, const RadionuclideIn //Get Radiopharmaceutical Info Sequence const gdcm::DataElement &de = file.GetDataSet().GetDataElement(radiopharm_info_seq_tag); const gdcm::SequenceOfItems *sqi = de.GetValueAsSQ(); - const gdcm::Item &item = sqi->GetItem(1); + const gdcm::Item &item = sqi->GetItem(sequence_idx); //Get Radiopnuclide Code Sequence @@ -174,7 +183,7 @@ stir::Succeeded GetRadionuclideInfo(const gdcm::File &file, const RadionuclideIn const gdcm::DataElement &element = item.GetDataElement(radionuclide_code_seq_tag); if(element.GetVL()>0){ const gdcm::SequenceOfItems *sqi2 = element.GetValueAsSQ(); - const gdcm::Item &item2 = sqi2->GetItem(1); + const gdcm::Item &item2 = sqi2->GetItem(sequence_idx); //std::cout<< "num items"<< sqi2->GetNumberOfItems(); //std::cout << item2 << std::endl; @@ -288,20 +297,30 @@ stir::Succeeded SPECTDICOMData::open_dicom_file(bool is_planar) actual_frame_duration = std::stoi(actual_frame_duration_as_string); } - if (GetEnergyWindowInfo(file, EnergyWindowInfo::WindowName , energy_window_name) == stir::Succeeded::yes){ - // actual_frame_duration = std::stoi(matrix_size_as_string); - std::cout << "Energy window: " << energy_window_name << std::endl; - } - - if (GetEnergyWindowInfo(file, EnergyWindowInfo::LowerThreshold , lower_window_as_string) == stir::Succeeded::yes){ - lower_en_window_thres = std::stof(lower_window_as_string); - std::cout << "Lower energy window limit: " << std::fixed << std::setprecision(6) << lower_en_window_thres << std::endl; - } - - if (GetEnergyWindowInfo(file, EnergyWindowInfo::UpperThreshold , upper_window_as_string) == stir::Succeeded::yes){ - upper_en_window_thres = std::stof(upper_window_as_string); - std::cout << "Upper energy window limit: " << std::fixed << std::setprecision(6) << upper_en_window_thres << std::endl; + { + 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; + } + + 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; @@ -338,7 +357,7 @@ stir::Succeeded SPECTDICOMData::open_dicom_file(bool is_planar) } -stir::Succeeded SPECTDICOMData::get_interfile_header(std::string &output_header) const{ +stir::Succeeded SPECTDICOMData::get_interfile_header(std::string &output_header, const std::string& data_filename, const int dataset_num) const{ std::string data_filename_only = data_filename; @@ -354,7 +373,9 @@ stir::Succeeded SPECTDICOMData::get_interfile_header(std::string &output_header) ss << "!imaging modality := nucmed" << std::endl; ss << "!version of keys := 3.3" << std::endl; ss << "name of data file := " << data_filename_only << std::endl; - ss << "data offset in bytes := 0" << 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 + << std::endl; ss << std::endl; ss << "!GENERAL IMAGE DATA :=" << std::endl; @@ -366,6 +387,10 @@ stir::Succeeded SPECTDICOMData::get_interfile_header(std::string &output_header) ss << "isotope name:= " << this->isotope_name<< std::endl; ss << std::endl; + // one energy window per header + ss << "number of energy windows :=1\n"; + ss <<"energy window lower level[1] := " << this->lower_en_window_thres[dataset_num] << "\n"; + ss <<"energy window upper level[1] := " << this->upper_en_window_thres[dataset_num] << "\n"; ss << "!SPECT STUDY (General) :=" << std::endl; ss << "number of dimensions := 2" << std::endl; ss << "matrix axis label [2] := axial coordinate" << std::endl; @@ -465,41 +490,45 @@ stir::Succeeded SPECTDICOMData::get_proj_data(const std::string &output_file) co int main(int argc, char * argv[]) { - if ( argc!=5) { - std::cerr << "Usage: " << argv[0] << " "<< " "<<""<<"is_planar?\n"; + if ( argc!=4) { + std::cerr << "Usage: " << argv[0] << " sinogram(dcm)> is_planar?\n"; exit(EXIT_FAILURE); } std::unique_ptr - spect(new SPECTDICOMData(argv[1])); + spect(new SPECTDICOMData(argv[2])); + const std::string output_prefix(argv[1]); - spect->is_planar=atoi(argv[4]); + spect->is_planar=atoi(argv[3]); try{ if (spect->open_dicom_file(spect->is_planar) == stir::Succeeded::no){ std::cerr << "Failed to read!" << std::endl; return EXIT_FAILURE; } - } catch(const std::string& e){ - std::cerr << e << std::endl; + } catch(const std::exception& e){ + std::cerr << e.what() << std::endl; return EXIT_FAILURE; } std::cout << "Read successfully!" << std::endl; - spect->set_data_filename(argv[3]); - stir::Succeeded s = spect->get_proj_data(argv[3]); - - std::string header; - s = spect->get_interfile_header(header); - - std::cout << header << std::endl; - - std::filebuf fb; - fb.open (argv[2],std::ios::out); - std::ostream os(&fb); - os << header; - fb.close(); + const std::string data_filename(output_prefix + ".s"); + stir::Succeeded s = spect->get_proj_data(data_filename); + + for (int w=0; wnum_energy_windows; ++w) + { + std::string header; + s = spect->get_interfile_header(header, data_filename, w); + + std::cout << header << std::endl; + const std::string header_filename = output_prefix + spect->energy_window_name[w] + ".hdr"; + std::filebuf fb; + fb.open (header_filename,std::ios::out); + std::ostream os(&fb); + os << header; + fb.close(); + } return EXIT_SUCCESS; } From 88315bd6bddd9a8279ad567d8152b7b102bc7a6b Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Fri, 31 Mar 2023 20:13:18 +0100 Subject: [PATCH 4/5] check SPECT number of frames add warning if values don't fit --- src/utilities/SPECT_dicom_to_interfile.cxx | 56 ++++++++++++++++------ 1 file changed, 42 insertions(+), 14 deletions(-) diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx index 700694c25..c778a1836 100644 --- a/src/utilities/SPECT_dicom_to_interfile.cxx +++ b/src/utilities/SPECT_dicom_to_interfile.cxx @@ -1,8 +1,26 @@ +/* + Copyright (C) 2018, 2023, University College London + Copyright (C) 2019-2023, National Physical Laboratory + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +/*! + \file + \ingroup utilities + \brief Convert SPECT DICOM projection data to Interfile + + \author Benjamin Thomas + \author Daniel Deidda + \author Kris Thielemans +*/ #include #include #include -#include #include #include @@ -28,13 +46,14 @@ class SPECTDICOMData stir::Succeeded open_dicom_file(bool is_planar); bool is_planar; 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; - int num_of_projections = 0; float angular_step = 0.0f; int actual_frame_duration = 0; //frame duration in msec int num_of_rotations = 0; @@ -244,6 +263,15 @@ stir::Succeeded 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; + } + } + this->num_of_projections=1; if(!is_planar){ if (GetDICOMTagInfo(file, gdcm::Tag(0x0054,0x0053), no_of_proj_as_str) == stir::Succeeded::yes){ @@ -464,7 +492,7 @@ stir::Succeeded SPECTDICOMData::get_proj_data(const std::string &output_file) co outfile.close();*/ uint64_t len0 = (uint64_t)bv->GetLength()/2; - std::cout << "Length = " << len0 << std::endl; + std::cout << "Number of data points = " << len0 << std::endl; std::vector pixel_data_as_float; @@ -478,8 +506,6 @@ stir::Succeeded SPECTDICOMData::get_proj_data(const std::string &output_file) co ct++; } - std::cout << "pixel_data_as_float length = " << pixel_data_as_float.size() << std::endl; - std::ofstream final_out(output_file.c_str(), std::ios::out | std::ofstream::binary); final_out.write(reinterpret_cast(&pixel_data_as_float[0]), pixel_data_as_float.size() * sizeof(float)); final_out.close(); @@ -495,14 +521,13 @@ int main(int argc, char * argv[]) exit(EXIT_FAILURE); } - std::unique_ptr - spect(new SPECTDICOMData(argv[2])); + SPECTDICOMData spect(argv[2]); const std::string output_prefix(argv[1]); - spect->is_planar=atoi(argv[3]); + spect.is_planar=atoi(argv[3]); try{ - if (spect->open_dicom_file(spect->is_planar) == stir::Succeeded::no){ + if (spect.open_dicom_file(spect.is_planar) == stir::Succeeded::no){ std::cerr << "Failed to read!" << std::endl; return EXIT_FAILURE; } @@ -513,16 +538,19 @@ int main(int argc, char * argv[]) std::cout << "Read successfully!" << std::endl; + if (spect.num_of_projections * spect.num_energy_windows != spect.num_frames) + stir::warning(boost::format("num_projections (%1%) * num_energy_windows (%2%) != num_frames (%3%)\n" + "This could be gated or dynamic data. The interfile header(s) will be incorrect!") + % spect.num_of_projections % spect.num_energy_windows % spect.num_frames); const std::string data_filename(output_prefix + ".s"); - stir::Succeeded s = spect->get_proj_data(data_filename); + stir::Succeeded s = spect.get_proj_data(data_filename); - for (int w=0; wnum_energy_windows; ++w) + for (int w=0; wget_interfile_header(header, data_filename, w); + s = spect.get_interfile_header(header, data_filename, w); - std::cout << header << std::endl; - const std::string header_filename = output_prefix + spect->energy_window_name[w] + ".hdr"; + const std::string header_filename = output_prefix + spect.energy_window_name[w] + ".hdr"; std::filebuf fb; fb.open (header_filename,std::ios::out); std::ostream os(&fb); From cfc52afc70e9f6a41d4e2dc6deced0e77190ad8b Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Fri, 31 Mar 2023 20:26:59 +0100 Subject: [PATCH 5/5] only build SPECT DICOM file if ITK is present --- src/utilities/CMakeLists.txt | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index 4f4c1702c..3d69eb71e 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -73,9 +73,12 @@ set(${dir_EXE_SOURCES} write_sinogram_to_txt.cxx find_sum_projection_of_viewgram_and_sinogram.cxx separate_true_from_random_scatter_for_necr.cxx - SPECT_dicom_to_interfile.cxx ) +if (ITK_FOUND) + list(APPEND ${dir_EXE_SOURCES} SPECT_dicom_to_interfile.cxx) +endif() + if (AVW_FOUND) list(APPEND ${dir_EXE_SOURCES} conv_AVW.cxx) endif()