diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index 65fe3c290..3d69eb71e 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -75,6 +75,10 @@ set(${dir_EXE_SOURCES} separate_true_from_random_scatter_for_necr.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() diff --git a/src/utilities/SPECT_dicom_to_interfile.cxx b/src/utilities/SPECT_dicom_to_interfile.cxx new file mode 100644 index 000000000..c778a1836 --- /dev/null +++ b/src/utilities/SPECT_dicom_to_interfile.cxx @@ -0,0 +1,562 @@ +/* + 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 "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, 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 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; + 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; + 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 ; + + 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; +}; + +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. + + //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) { + 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 (...){ + 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, const int sequence_idx){ + + if (request == EnergyWindowInfo::WindowName){ + return GetDICOMTagInfo(file, gdcm::Tag(0x0054,0x0018), dst, sequence_idx); + } + + 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(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; + + 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 (...){ + 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, const int sequence_idx){ + +// 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(sequence_idx); + + + //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(sequence_idx); + //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 (...){ + 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; + + + { + 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){ + 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 << "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){ + uint16_t val = *ptr; + pixel_data_as_float.push_back((float)val); + ptr++; + ct++; + } + + 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!=4) { + std::cerr << "Usage: " << argv[0] << " sinogram(dcm)> is_planar?\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){ + std::cerr << "Failed to read!" << std::endl; + return EXIT_FAILURE; + } + } catch(const std::exception& e){ + std::cerr << e.what() << std::endl; + return EXIT_FAILURE; + } + + 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); + + for (int w=0; w