diff --git a/recon_test_pack/root_header.hroot b/recon_test_pack/root_header.hroot index 354c5df6ea..d21a1265d6 100644 --- a/recon_test_pack/root_header.hroot +++ b/recon_test_pack/root_header.hroot @@ -34,6 +34,14 @@ number of crystals_Z := 4 Singles readout depth := 1 exclude scattered events := ${EXCLUDE_SCATTERED} + +; +; If set only a certain order of scatter will be kept +exclude scattered events := ${MAXIMUM_ORDER_OF_SCATTER} + +; +; If set the random events will be skipped + exclude random events := ${EXCLUDE_RANDOM} check energy window information := 1 ; Default 1 offset (num of detectors) := 0 diff --git a/src/IO/InputStreamFromROOTFile.cxx b/src/IO/InputStreamFromROOTFile.cxx index 749eca9789..d612aa013a 100644 --- a/src/IO/InputStreamFromROOTFile.cxx +++ b/src/IO/InputStreamFromROOTFile.cxx @@ -42,12 +42,13 @@ InputStreamFromROOTFile() InputStreamFromROOTFile:: InputStreamFromROOTFile(std::string filename, std::string chain_name, - bool exclude_scattered, bool exclude_randoms, - float low_energy_window, float up_energy_window, + int maximum_order_of_scatter, bool exclude_randoms, + float low_energy_window_1, float up_energy_window_1, + float low_energy_window_2, float up_energy_window_2, int offset_dets) : filename(filename), chain_name(chain_name), - exclude_scattered(exclude_scattered), exclude_randoms(exclude_randoms), - low_energy_window(low_energy_window), up_energy_window(up_energy_window), offset_dets(offset_dets) + maximum_order_of_scatter(maximum_order_of_scatter), exclude_randoms(exclude_randoms), + low_energy_window_1(low_energy_window_1), up_energy_window_1(up_energy_window_1),low_energy_window_2(low_energy_window_2), up_energy_window_2(up_energy_window_2), offset_dets(offset_dets) { set_defaults(); error("This constructor is incorrect"); //TODO set_defaults() will override the above @@ -60,13 +61,16 @@ InputStreamFromROOTFile::set_defaults() { starting_stream_position = 0; singles_readout_depth = -1; + maximum_order_of_scatter = -1; exclude_nonrandom = false; exclude_scattered = false; exclude_unscattered = false; exclude_randoms = false; check_energy_window_information = true; - low_energy_window = 0.f; - up_energy_window = 1000.f; + low_energy_window_1 = 0.f; + up_energy_window_1 = 1000.f; + low_energy_window_2 = 0.f; + up_energy_window_2 = 1000.f; read_optional_root_fields=false; crystal_repeater_x = -1; crystal_repeater_y = -1; @@ -81,14 +85,17 @@ InputStreamFromROOTFile::initialise_keymap() this->parser.add_key("name of data file", &this->filename); this->parser.add_key("Singles readout depth", &this->singles_readout_depth); this->parser.add_key("name of input TChain", &this->chain_name); + this->parser.add_key("maximum order of scatter", &this->maximum_order_of_scatter); this->parser.add_key("exclude non-random events", &this->exclude_nonrandom); this->parser.add_key("exclude scattered events", &this->exclude_scattered); this->parser.add_key("exclude unscattered events", &this->exclude_unscattered); this->parser.add_key("exclude random events", &this->exclude_randoms); this->parser.add_key("check energy window information", &this->check_energy_window_information); this->parser.add_key("offset (num of detectors)", &this->offset_dets); - this->parser.add_key("low energy window (keV)", &this->low_energy_window); - this->parser.add_key("upper energy window (keV)", &this->up_energy_window); + this->parser.add_key("low energy window 1 (MeV)", &this->low_energy_window_1); + this->parser.add_key("upper energy window 1 (MeV)", &this->up_energy_window_1); + this->parser.add_key("low energy window 2 (MeV)", &this->low_energy_window_2); + this->parser.add_key("upper energy window 2 (MeV)", &this->up_energy_window_2); this->parser.add_key("read optional ROOT fields", &this->read_optional_root_fields); this->parser.add_key("number of crystals X", &this->crystal_repeater_x); @@ -203,7 +210,7 @@ InputStreamFromROOTFile::check_brentry_randoms_scatter_energy_conditions(Long64_ return false; // Scattered/unscattered event exclusion condition. - if (this->exclude_scattered || this->exclude_unscattered) { + if (this->exclude_scattered || this->exclude_unscattered || this->maximum_order_of_scatter >= 0) { GetEntryCheck(br_comptonPhantom1->GetEntry(brentry)); GetEntryCheck(br_comptonPhantom2->GetEntry(brentry)); @@ -212,6 +219,8 @@ InputStreamFromROOTFile::check_brentry_randoms_scatter_energy_conditions(Long64_ return false; if (this->exclude_unscattered && (this->comptonphantom1 == 0 && this->comptonphantom2 == 0)) return false; + if ((this->maximum_order_of_scatter >= 0) && (this->comptonphantom1 + this->comptonphantom2 > this->maximum_order_of_scatter)) + return false; } // Trues/Random event exclusion condition. @@ -234,10 +243,10 @@ InputStreamFromROOTFile::check_brentry_randoms_scatter_energy_conditions(Long64_ GetEntryCheck(br_energy2->GetEntry(brentry)); // Check both energy values are within window - if (this->get_energy1_in_keV() < this->low_energy_window || - this->get_energy1_in_keV() > this->up_energy_window || - this->get_energy2_in_keV() < this->low_energy_window || - this->get_energy2_in_keV() > this->up_energy_window) + if (this->get_energy1_in_keV() < this->low_energy_window_1 || + this->get_energy1_in_keV() > this->up_energy_window_1 || + this->get_energy2_in_keV() < this->low_energy_window_2 || + this->get_energy2_in_keV() > this->up_energy_window_2) { return false; } diff --git a/src/IO/InputStreamFromROOTFileForCylindricalPET.cxx b/src/IO/InputStreamFromROOTFileForCylindricalPET.cxx index 053aa62c2c..f5e9536248 100644 --- a/src/IO/InputStreamFromROOTFileForCylindricalPET.cxx +++ b/src/IO/InputStreamFromROOTFileForCylindricalPET.cxx @@ -30,8 +30,9 @@ InputStreamFromROOTFileForCylindricalPET(std::string _filename, int submodule_repeater_x, int submodule_repeater_y, int submodule_repeater_z, int module_repeater_x, int module_repeater_y, int module_repeater_z, int rsector_repeater, - bool _exclude_scattered, bool _exclude_randoms, - float _low_energy_window, float _up_energy_window, + int _maximum_order_of_scatter, bool _exclude_randoms, + float _low_energy_window_1, float _up_energy_window_1, + float _low_energy_window_2, float _up_energy_window_2, int _offset_dets): base_type(), crystal_repeater_x(crystal_repeater_x), crystal_repeater_y(crystal_repeater_y), crystal_repeater_z(crystal_repeater_z), @@ -44,10 +45,12 @@ InputStreamFromROOTFileForCylindricalPET(std::string _filename, filename = _filename; chain_name = _chain_name; - exclude_scattered = _exclude_scattered; + maximum_order_of_scatter = _maximum_order_of_scatter; exclude_randoms = _exclude_randoms; - low_energy_window = _low_energy_window; - up_energy_window = _up_energy_window; + low_energy_window_1 = _low_energy_window_1; + up_energy_window_1 = _up_energy_window_1; + low_energy_window_2 = _low_energy_window_2; + up_energy_window_2 = _up_energy_window_2; offset_dets = _offset_dets; half_block = module_repeater_y * submodule_repeater_y * crystal_repeater_y / 2 - 1; @@ -123,7 +126,8 @@ get_next_record(CListRecordROOT& record) record.init_from_data(ring1, ring2, crystal1, crystal2, time1, time2, - eventID1, eventID2); + eventID1, eventID2, + energy1,energy2); } std::string diff --git a/src/IO/InputStreamFromROOTFileForECATPET.cxx b/src/IO/InputStreamFromROOTFileForECATPET.cxx index 2bded5af8c..6adb417647 100644 --- a/src/IO/InputStreamFromROOTFileForECATPET.cxx +++ b/src/IO/InputStreamFromROOTFileForECATPET.cxx @@ -30,8 +30,9 @@ InputStreamFromROOTFileForECATPET(std::string _filename, std::string _chain_name, int crystal_repeater_x, int crystal_repeater_y, int crystal_repeater_z, int block_repeater_y, int block_repeater_z, - bool _exclude_scattered, bool _exclude_randoms, - float _low_energy_window, float _up_energy_window, + int _maximum_order_of_scater, bool _exclude_randoms, + float _low_energy_window_1, float _up_energy_window_1, + float _low_energy_window_2, float _up_energy_window_2, int _offset_dets): base_type(), crystal_repeater_x(crystal_repeater_x), crystal_repeater_y(crystal_repeater_y), crystal_repeater_z(crystal_repeater_z), @@ -42,10 +43,12 @@ InputStreamFromROOTFileForECATPET(std::string _filename, filename = _filename; chain_name = _chain_name; - exclude_scattered = _exclude_scattered; + maximum_order_of_scatter = _maximum_order_of_scater; exclude_randoms = _exclude_randoms; - low_energy_window = _low_energy_window; - up_energy_window = _up_energy_window; + low_energy_window_1 = _low_energy_window_1; + up_energy_window_1 = _up_energy_window_1; + low_energy_window_2 = _low_energy_window_2; + up_energy_window_2 = _up_energy_window_2; offset_dets = _offset_dets; half_block = crystal_repeater_y / 2 - 1; @@ -107,7 +110,8 @@ get_next_record(CListRecordROOT& record) record.init_from_data(ring1, ring2, crystal1, crystal2, time1, time2, - eventID1, eventID2); + eventID1, eventID2, + energy1,energy2); } std::string diff --git a/src/IO/InterfileHeader.cxx b/src/IO/InterfileHeader.cxx index 94925dd08b..67d0177498 100644 --- a/src/IO/InterfileHeader.cxx +++ b/src/IO/InterfileHeader.cxx @@ -99,6 +99,7 @@ void MinimalInterfileHeader::set_version_specific_keys() InterfileHeader::InterfileHeader() : MinimalInterfileHeader() { + number_format_values.push_back("bit"); number_format_values.push_back("ascii"); number_format_values.push_back("signed integer"); @@ -169,6 +170,7 @@ InterfileHeader::InterfileHeader() + add_key("name of data file", &data_file_name); add_key("originating system", &exam_info_sptr->originating_system); ignore_key("GENERAL DATA"); @@ -222,8 +224,8 @@ InterfileHeader::InterfileHeader() // support for Louvain la Neuve's extension of 3.3 add_key("quantification units", &lln_quantification_units); - add_key("number of energy windows", - KeyArgument::INT, (KeywordProcessor)&InterfileHeader::read_num_energy_windows,&num_energy_windows); + add_key("energy window pair", KeyArgument::LIST_OF_INTS, (KeywordProcessor)&InterfileHeader::en_window_pair_set, &energy_window_pair); + add_key("number of energy windows", KeyArgument::INT, (KeywordProcessor)&InterfileHeader::read_num_energy_windows,&num_energy_windows); add_vectorised_key("energy window lower level", &lower_en_window_thresholds); add_vectorised_key("energy window upper level", &upper_en_window_thresholds); @@ -231,6 +233,7 @@ InterfileHeader::InterfileHeader() add_key("start horizontal bed position (mm)", &bed_position_horizontal); bed_position_vertical = 0.F; add_key("start vertical bed position (mm)", &bed_position_vertical); + } void InterfileHeader::set_version_specific_keys() @@ -240,16 +243,17 @@ void InterfileHeader::set_version_specific_keys() { info("Setting energy window keys as in STIR3.0"); // only a single energy window, and non-vectorised - remove_key("energy window lower level"); + /*remove_key("energy window lower level"); remove_key("energy window upper level"); add_key("energy window lower level", &lower_en_window_thresholds[0]); - add_key("energy window upper level", &upper_en_window_thresholds[0]); + add_key("energy window upper level", &upper_en_window_thresholds[0]);*/ } } // MJ 17/05/2000 made bool bool InterfileHeader::post_processing() { + if(type_of_data_index<0) { warning("Interfile Warning: 'type_of_data' keyword required"); @@ -373,19 +377,45 @@ bool InterfileHeader::post_processing() lln_quantification_units); } } // lln_quantification_units - if (num_energy_windows>0) - { - if (num_energy_windows>1) - warning("Currently only reading the first energy window."); - if (upper_en_window_thresholds[0] > 0 && lower_en_window_thresholds[0] > 0 ) - { - exam_info_sptr->set_high_energy_thres(upper_en_window_thresholds[0]); - exam_info_sptr->set_low_energy_thres(lower_en_window_thresholds[0]); - } - } + + + //set the lower and the higher energy thresholds for all the energy windows available. Default: 1. + + //set the number of energy windows + + exam_info_sptr->set_num_energy_windows(num_energy_windows); + + //set the number of energy window pair + if (energy_window_pair.size() > 0) { + + if (energy_window_pair.size() != 2) + error("should have two."); + if (energy_window_pair[0] < 0) + error("first window should be >= 0."); + if (energy_window_pair[1] < 0) + error("second window should be >= 0."); + if (energy_window_pair[0] > num_energy_windows) + error("The selected window %d exceeds the number of energy windows %d.\n",energy_window_pair[0],num_energy_windows); + if (energy_window_pair[1] > num_energy_windows) + error("The selected window %d exceeds the number of energy windows %d.\n",energy_window_pair[1],num_energy_windows); + + exam_info_sptr->set_energy_window_pair(energy_window_pair); + } + + //set the high and low energy window threshold + + for (int i = 0; i < num_energy_windows; ++i) + { + + if (upper_en_window_thresholds[i] >0 && lower_en_window_thresholds[i] >0) + { + exam_info_sptr->set_high_energy_thres(upper_en_window_thresholds[i],i); + exam_info_sptr->set_low_energy_thres(lower_en_window_thresholds[i],i); + } + } exam_info_sptr->time_frame_definitions = - TimeFrameDefinitions(image_relative_start_times, image_durations); + TimeFrameDefinitions(image_relative_start_times, image_durations); return false; @@ -401,6 +431,14 @@ void InterfileHeader::read_matrix_info() } +void InterfileHeader::en_window_pair_set() +{ + + energy_window_pair.resize(2); + set_variable(); + + +} void InterfileHeader::read_num_energy_windows() { set_variable(); @@ -428,7 +466,7 @@ void InterfileHeader::set_type_of_data() &PET_data_type_values); ignore_key("process status"); ignore_key("IMAGE DATA DESCRIPTION"); - // TODO rename keyword + // TODO rename keyword add_vectorised_key("data offset in bytes", &data_offset_each_dataset); } @@ -671,6 +709,9 @@ void InterfilePDFSHeader::resize_segments_and_set() } + + + int InterfilePDFSHeader::find_storage_order() { @@ -1409,7 +1450,7 @@ bool InterfilePDFSHeader::post_processing() scanner_ptr_from_file->parameter_info().c_str()); } - + // float azimuthal_angle_sampling =_PI/num_views; diff --git a/src/IO/interfile.cxx b/src/IO/interfile.cxx index 2dbe67b2ff..8e6e385efd 100644 --- a/src/IO/interfile.cxx +++ b/src/IO/interfile.cxx @@ -469,6 +469,11 @@ static void write_interfile_time_frame_definitions(std::ostream& output_header, // Write energy window lower and upper thresholds, if they are not -1 static void write_interfile_energy_windows(std::ostream& output_header, const ExamInfo& exam_info) { + //Write the number of energy windows + if (exam_info.get_num_energy_windows() > 0) + { + output_header <<"number of energy windows := " << exam_info.get_num_energy_windows() << '\n'; +/* if (exam_info.get_high_energy_thres() > 0 && exam_info.get_low_energy_thres() >= 0) { @@ -477,8 +482,35 @@ static void write_interfile_energy_windows(std::ostream& output_header, const Ex exam_info.get_low_energy_thres() << '\n'; output_header << "energy window upper level[1] := " << exam_info.get_high_energy_thres() << '\n'; +*/ } -} + + else + { // need to write this anyway to allow vectored keys below + //output_header <<"number of energy windows := 1"; + } + + //Write energy window thresholds + for (unsigned int num_windows = 0; num_windows < exam_info.get_num_energy_windows(); ++num_windows) + { + if (exam_info.get_high_energy_thres(num_windows) > 0 && + exam_info.get_low_energy_thres(num_windows) >= 0) + { + output_header << "energy window lower level [" << num_windows +1 << "] := " << + exam_info.get_low_energy_thres(num_windows) << '\n'; + output_header << "energy window upper level [" << num_windows +1 << "] := " << + exam_info.get_high_energy_thres(num_windows) << '\n'; + } + } + + if (exam_info.get_energy_window_pair().first> 0 && + exam_info.get_energy_window_pair().second> 0) + { + output_header << "energy window pair :="<<" {"<< exam_info.get_energy_window_pair().first << + ',' << exam_info.get_energy_window_pair().second <<"}\n"; + } + + } // Write data type descriptions (if there are any) static void write_interfile_image_data_descriptions(std::ostream& output_header, const std::vector& data_type_descriptions) @@ -1152,15 +1184,19 @@ read_interfile_PDFS(const string& filename, const ios::openmode open_mode) { ifstream image_stream(filename.c_str()); + + if (!image_stream) { error("read_interfile_PDFS: couldn't open file %s\n", filename.c_str()); } char directory_name[max_filename_length]; + get_directory_name(directory_name, filename.c_str()); - + return read_interfile_PDFS(image_stream, directory_name, open_mode); + } Succeeded @@ -1509,6 +1545,7 @@ write_basic_interfile_PDFS_header(const string& header_file_name, write_interfile_time_frame_definitions(output_header, pdfs.get_exam_info()); write_interfile_energy_windows(output_header, pdfs.get_exam_info()); + if (pdfs.get_scale_factor()!=1.F) output_header <<"image scaling factor[1] := " <time_frame_definitions.get_duration(1) << ")\n"; } + // TODOENERGY s << "number of energy windows:=1\n" << "energy window lower level[1] := " << this->get_low_energy_thres() << '\n' @@ -69,8 +70,8 @@ ExamInfo::parameter_info() const bool ExamInfo::operator == (const ExamInfo &p1) const { - return abs(this->up_energy_thres - p1.up_energy_thres )<=1 && // keV - abs(this->low_energy_thres - p1.low_energy_thres) <=1 &&// keV + const bool equal_without_en = + this->num_windows==p1.num_windows && this->radionuclide==p1.radionuclide && this->time_frame_definitions==p1.time_frame_definitions && // this->branching_ratio==p1.branching_ratio && @@ -78,7 +79,18 @@ ExamInfo::operator == (const ExamInfo &p1) const { abs(this->calibration_factor/p1.calibration_factor -1.)<=1E-3) && this->imaging_modality==p1.imaging_modality && this->patient_position==p1.patient_position && - abs(this->start_time_in_secs_since_1970 - p1.start_time_in_secs_since_1970)<=.5;// sec + abs(this->start_time_in_secs_since_1970 - p1.start_time_in_secs_since_1970)<=.5;// sec + if (!equal_without_en) + return false; + for (int en=0; en < this->num_windows; ++en) + { + if (abs(this->up_energy_thres[en] - p1.up_energy_thres[en] ) > 1) // keV + return false; + if (abs(this->low_energy_thres[en] - p1.low_energy_thres[en]) > 1) // keV + return false; + } + return this->en_win_pair[0] == p1.en_win_pair[0] && + this->en_win_pair[1] == p1.en_win_pair[1]; } END_NAMESPACE_STIR diff --git a/src/buildblock/KeyParser.cxx b/src/buildblock/KeyParser.cxx index 3cb75ec335..c84fd8cd05 100644 --- a/src/buildblock/KeyParser.cxx +++ b/src/buildblock/KeyParser.cxx @@ -528,6 +528,7 @@ void KeyParser::add_key(const string& keyword, void* variable, const ASCIIlist_type * const list_of_values) { + add_in_keymap(keyword, map_element(t, function, variable, 0, list_of_values)); } @@ -549,6 +550,16 @@ void KeyParser::add_key(const string& keyword, add_in_keymap(keyword, map_element(t, &KeyParser::set_variable, variable, 0, list_of_values)); } +#if 0 +void KeyParser::add_vectored_key(const string& keyword, + KeyArgument::type t, + KeywordProcessor function, + void* variable, + const ASCIIlist_type * const list_of_values) +{ + add_in_keymap(keyword, map_element(t, function, variable, list_of_values, true)); +} +#endif void KeyParser::add_key(const string& keyword, KeyArgument::type t, void* variable, @@ -858,6 +869,12 @@ Succeeded KeyParser::parse_value_in_line(const string& line, const bool write_wa // maps keyword to appropriate map_element (sets current) if(map_keyword(keyword)==Succeeded::yes) { +#if 0 + if (current_index != 0 && current->is_vector.is_false()) + error(boost::format("Error parsing keyword %1%: should not be vectored"), keyword); + if (current_index == 0 && current->is_vector.is_true()) + error(boost::format("Error parsing keyword %1%: should be vectored"), keyword); +#endif switch(current->type) // depending on the par_type, gets the correct value from the line { // and sets the right temporary variable case KeyArgument::NONE : diff --git a/src/buildblock/ProjData.cxx b/src/buildblock/ProjData.cxx index 270f866802..923a77fb1a 100644 --- a/src/buildblock/ProjData.cxx +++ b/src/buildblock/ProjData.cxx @@ -102,6 +102,8 @@ ProjData:: read_from_file(const string& filename, const std::ios::openmode openmode) { + + std::string actual_filename = filename; // parse filename to see if it's like filename,options { @@ -115,7 +117,6 @@ read_from_file(const string& filename, fstream * input = new fstream(actual_filename.c_str(), openmode | ios::binary); if (! *input) error("ProjData::read_from_file: error opening file %s", actual_filename.c_str()); - const FileSignature file_signature(actual_filename); const char * signature = file_signature.get_signature(); @@ -190,7 +191,9 @@ read_from_file(const string& filename, #ifndef NDEBUG warning("ProjData::read_from_file trying to read %s as Interfile", filename.c_str()); #endif + shared_ptr ptr(read_interfile_PDFS(filename, openmode)); + if (!is_null_ptr(ptr)) return ptr; } diff --git a/src/buildblock/ProjDataFromStream.cxx b/src/buildblock/ProjDataFromStream.cxx index aaf5894935..27f8d0e282 100644 --- a/src/buildblock/ProjDataFromStream.cxx +++ b/src/buildblock/ProjDataFromStream.cxx @@ -848,6 +848,8 @@ ProjDataFromStream::get_offset_segment(const int segment_num) const SegmentBySinogram ProjDataFromStream::get_segment_by_sinogram(const int segment_num) const { + + if(is_null_ptr(sino_stream)) { error("ProjDataFromStream::get_segment_by_sinogram: stream ptr is 0\n"); @@ -878,6 +880,7 @@ ProjDataFromStream::get_segment_by_sinogram(const int segment_num) const succeeded = read_data(*sino_stream, segment, on_disk_data_type, scale, on_disk_byte_order); } } // end of critical section + if (succeeded == Succeeded::no) error("ProjDataFromStream: error reading data\n"); if(scale != 1) @@ -914,6 +917,7 @@ ProjDataFromStream::get_segment_by_view(const int segment_num) const #ifdef STIR_OPENMP #pragma omp critical(PROJDATAFROMSTREAMIO) #endif + { sino_stream->seekg(segment_offset, ios::beg); if (! *sino_stream) diff --git a/src/buildblock/ProjDataInfoCylindricalArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalArcCorr.cxx index 27b1396e74..2aad69527e 100644 --- a/src/buildblock/ProjDataInfoCylindricalArcCorr.cxx +++ b/src/buildblock/ProjDataInfoCylindricalArcCorr.cxx @@ -115,9 +115,13 @@ ProjDataInfoCylindricalArcCorr::parameter_info() const Bin ProjDataInfoCylindricalArcCorr:: -get_bin(const LOR& lor) const +get_bin(const LOR& lor, const std::pair &energy_window_pair) const { + if ((energy_window_pair.first!= 1)||(energy_window_pair.second!= 1)) + { + error("TODO NOT IMPLEMENTED YET FOR MULTIPLE ENERGY WINDOWS"); + } Bin bin; LORInAxialAndSinogramCoordinates lor_coords; if (lor.change_representation(lor_coords, get_ring_radius()) == Succeeded::no) diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx index f9089707b4..3c52edbac4 100644 --- a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx @@ -540,7 +540,7 @@ find_bin_given_cartesian_coordinates_of_detection(Bin& bin, Bin ProjDataInfoCylindricalNoArcCorr:: -get_bin(const LOR& lor) const +get_bin(const LOR& lor,const std::pair &energy_window_pair) const { Bin bin; #ifndef STIR_DEVEL @@ -568,7 +568,7 @@ get_bin(const LOR& lor) const if (ring1 >=0 && ring1=0 && ring2= get_min_tangential_pos_num() && bin.tangential_pos_num() <= get_max_tangential_pos_num()) { diff --git a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx index bff19f716b..589625faf2 100644 --- a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx @@ -394,9 +394,13 @@ find_cartesian_coordinates_given_scanner_coordinates(CartesianCoordinate3D& lor) const +get_bin(const LOR& lor, const std::pair &energy_window_pair) const { - Bin bin; + if ((energy_window_pair.first!= 1)||(energy_window_pair.second!= 1)) + { + error("TODO NOT IMPLEMENTED YET FOR MULTIPLE ENERGY WINDOWS"); + } + Bin bin; const LORAs2Points & lor_as_2points = dynamic_cast &>(lor); diff --git a/src/buildblock/ProjDataInfoSubsetByView.cxx b/src/buildblock/ProjDataInfoSubsetByView.cxx index 7d5889d14b..65d7e480ad 100644 --- a/src/buildblock/ProjDataInfoSubsetByView.cxx +++ b/src/buildblock/ProjDataInfoSubsetByView.cxx @@ -221,9 +221,9 @@ ProjDataInfoSubsetByView::get_sampling_in_s(const Bin& bin) const } Bin -ProjDataInfoSubsetByView::get_bin(const LOR& lor) const +ProjDataInfoSubsetByView::get_bin(const LOR& lor, const std::pair &energy_window_pair) const { - return get_bin_from_original(this->org_proj_data_info_sptr->get_bin(lor)); + return get_bin_from_original(this->org_proj_data_info_sptr->get_bin(lor, energy_window_pair)); } bool diff --git a/src/buildblock/extend_projdata.cxx b/src/buildblock/extend_projdata.cxx index ae9496f01e..0b1f6084cb 100644 --- a/src/buildblock/extend_projdata.cxx +++ b/src/buildblock/extend_projdata.cxx @@ -33,6 +33,7 @@ namespace detail to find data in the negative segment if necessary. However, it needs testing if it would work for non-direct sinograms. */ + inline static Array<2,float> extend_sinogram_in_views(const Array<2,float>& sino_positive_segment, @@ -61,12 +62,12 @@ namespace detail if (fabs(min_phi)< .01) { - min_in[1]-=min_view_extension; + min_in[1]-=min_view_extension; //here the new min is lower than the original min_is_extended=true; } if (fabs(max_phi-(_PI-sampling_phi))<.01) { - max_in[1]+=max_view_extension; + max_in[1]+=max_view_extension; //here the new max is bigger than the original max_is_extended=true; } @@ -83,13 +84,15 @@ namespace detail { bool use_extension=false; int symmetric_view_num=0; - if (view_numorg_max_view_num && max_is_extended==true) + else if (view_num>org_max_view_num && max_is_extended==true) //if view number is bigger than the old maximum (right hand side?) it's been extended { + // std::cout << "view_num>org_max_view_num " << '\n'; use_extension=true; symmetric_view_num = view_num - num_views_for_180; } @@ -113,10 +116,50 @@ namespace detail input_extended_view[view_num][symmetric_max]; } } // loop over views + return input_extended_view; } -} // end of namespace detail + + + inline static + Array<2,float> + transpose_extend_sinogram_in_views(const Array<2,float>& sino_segment, + const ProjDataInfo& proj_data_info, + const float min_view_compression, const float max_view_compression) + { + //* Check if projdata are from 0 to pi-phi + + BasicCoordinate<2,int> min_in, max_in; + if (!sino_segment.get_regular_range(min_in, max_in)) + { + warning("input segment 0 should have a regular range"); + } + + + min_in[1]+=min_view_compression; //here the new min is bigger than the original + max_in[1]-=max_view_compression; //here the new max is smaller than the original + + const float min_phi = proj_data_info.get_phi(Bin(0,0,0,0)); + + const float sampling_phi = + proj_data_info.get_phi(Bin(0,1,0,0)) - min_phi; + const int num_views_for_180 = round(_PI/sampling_phi); + + IndexRange<2> compressed_range(min_in, max_in); + Array<2,float> input_compressed_view(compressed_range); + + for (int view_num=min_in[1]; view_num<=max_in[1]; ++view_num) + { + input_compressed_view[view_num]=sino_segment[view_num]; //here we cut everything bigger than max_in and smaller than min_in + + } // loop over views + + input_compressed_view[min_in[1]]+=sino_segment[min_in[1]+num_views_for_180];//boundary condition + return input_compressed_view; + } + +} Array<3,float> extend_segment_in_views(const SegmentBySinogram& sino, const int min_view_extension, const int max_view_extension) @@ -159,4 +202,46 @@ extend_sinogram_in_views(const Sinogram& sino, min_view_extension, max_view_extension); } +Array<3,float> +transpose_extend_segment_in_views(const SegmentBySinogram& sino, + const float min_view_compression, const float max_view_compression) +{ + if (sino.get_segment_num()!=0) + error("extend_segment with single segment works only for segment 0"); + + BasicCoordinate<3,int> min, max; + + min[1]=sino.get_min_axial_pos_num(); + max[1]=sino.get_max_axial_pos_num(); + min[2]=sino.get_min_view_num(); + max[2]=sino.get_max_view_num(); + min[3]=sino.get_min_tangential_pos_num(); + max[3]=sino.get_max_tangential_pos_num(); + const IndexRange<3> out_range(min,max); + Array<3,float> out(out_range); + for (int ax_pos_num=min[1]; ax_pos_num <=max[1] ; ++ax_pos_num) + { + out[ax_pos_num] = + detail:: + transpose_extend_sinogram_in_views(sino[ax_pos_num], + *(sino.get_proj_data_info_sptr()), + min_view_compression, max_view_compression); + } + return out; +} + +Array<2,float> +transpose_extend_sinogram_in_views(const Sinogram& sino, + const float min_view_compression, const float max_view_compression) +{ + if (sino.get_segment_num()!=0) + error("extend_segment with single segment works only for segment 0"); + + return + detail:: + transpose_extend_sinogram_in_views(sino, + *(sino.get_proj_data_info_sptr()), + min_view_compression, max_view_compression); +} + END_NAMESPACE_STIR diff --git a/src/buildblock/interpolate_projdata.cxx b/src/buildblock/interpolate_projdata.cxx index 66c36297a8..a0b7c0f3c9 100644 --- a/src/buildblock/interpolate_projdata.cxx +++ b/src/buildblock/interpolate_projdata.cxx @@ -32,8 +32,13 @@ #include "stir/interpolate_projdata.h" #include "stir/extend_projdata.h" #include "stir/numerics/sampling_functions.h" +#include "stir_experimental/motion/Transform3DObjectImageProcessor.h" +#include "stir_experimental/motion/transform_3d_object.h" +#include "stir_experimental/numerics/more_interpolators.h" #include +#include "stir/IO/write_data.h" + START_NAMESPACE_STIR namespace detail_interpolate_projdata @@ -61,6 +66,34 @@ namespace detail_interpolate_projdata return new_proj_data_info_sptr; } + static shared_ptr + make_extended_proj_data_info(const ProjDataInfo& proj_data_info) + { + + if (dynamic_cast(&proj_data_info) == NULL) + error("make_extended_proj_data is only appropriate for non-arccorrected data"); + + shared_ptr new_proj_data_info_sptr( + proj_data_info.clone()); + new_proj_data_info_sptr-> + set_num_views(proj_data_info.get_num_views()*2); + return new_proj_data_info_sptr; + } + + static shared_ptr + transpose_make_non_interleaved_proj_data_info(const ProjDataInfo& proj_data_info) + { + + // TODO: arc-correct? + if (dynamic_cast(&proj_data_info) == NULL) + error("make_non_interleaved_proj_data is only appropriate for non-arccorrected data"); + + shared_ptr new_proj_data_info_sptr( + proj_data_info.clone()); + new_proj_data_info_sptr-> + set_num_views(proj_data_info.get_num_views()/2); + return new_proj_data_info_sptr; + } // access Sinogram element with wrap-around and boundary conditions static float sino_element(const Sinogram& sinogram, const int view_num, const int tangential_pos_num) { @@ -83,7 +116,7 @@ namespace detail_interpolate_projdata assert(out_sinogram.get_min_view_num() == 0); assert(in_sinogram.get_min_view_num() == 0); - assert(out_sinogram.get_num_views() == in_sinogram.get_num_views()*2); + assert(out_sinogram.get_num_views() == in_sinogram.get_num_views()*2); //check that the views are dowbled assert(in_sinogram.get_segment_num() == 0); assert(out_sinogram.get_segment_num() == 0); @@ -95,6 +128,7 @@ namespace detail_interpolate_projdata tangential_pos_num <= out_sinogram.get_max_tangential_pos_num()-1; ++tangential_pos_num) { + // here we're filling a bigger grid. if ((view_num+tangential_pos_num)%2 == 0) { const int in_view_num = @@ -133,6 +167,76 @@ namespace detail_interpolate_projdata } #endif + static void + transpose_make_non_interleaved_sinogram(Sinogram& out_sinogram, + const Sinogram& in_sinogram) + { + if (is_null_ptr(dynamic_pointer_cast(in_sinogram.get_proj_data_info_sptr()))) + error("make_non_interleaved_proj_data is only appropriate for non-arccorrected data"); + + // Sinogram& out_sinogram2 = out_sinogram.get_empty_copy(); + assert(out_sinogram.get_min_view_num() == 0); + assert(in_sinogram.get_min_view_num() == 0); + assert(out_sinogram.get_num_views() == in_sinogram.get_num_views()/2); + assert(in_sinogram.get_segment_num() == 0); + assert(out_sinogram.get_segment_num() == 0); + + if (out_sinogram.get_num_views() != in_sinogram.get_num_views()/2) + error("views need to be reduced of a factor 2"); + + //const int in_num_views = in_sinogram.get_num_views(); + const int out_num_views = out_sinogram.get_num_views(); + + for (int view_num = in_sinogram.get_min_view_num()+1; + view_num <= in_sinogram.get_max_view_num()-1; + ++view_num) + { + // TODO don't put in outer tangential poss for now to avoid boundary stuff + for (int tangential_pos_num = in_sinogram.get_min_tangential_pos_num()+1; + tangential_pos_num <= in_sinogram.get_max_tangential_pos_num()-1; + ++tangential_pos_num) + { + const int out_view_num = + view_num%2==0 ? view_num/2 : (view_num+1)/2; + int first_view=in_sinogram.get_min_view_num(); + int last_view=in_sinogram.get_max_view_num(); + + // here i need to take the bigger grid an + if ((view_num+tangential_pos_num)%2 == 0) //if it's even + { + + out_sinogram[out_view_num%out_num_views][tangential_pos_num] += in_sinogram[view_num][tangential_pos_num] + + (in_sinogram[view_num-1][tangential_pos_num] + + in_sinogram[view_num+1][tangential_pos_num] + + in_sinogram[view_num][(tangential_pos_num-1)] + + in_sinogram[view_num][(tangential_pos_num+1)])/4; + + //BOUNDARY CONDITIONS FOR MIN AND MAX TAN POS + out_sinogram[out_view_num%out_num_views][in_sinogram.get_min_tangential_pos_num()] = + (in_sinogram[view_num][(in_sinogram.get_min_tangential_pos_num()+1)])/4; + out_sinogram[out_view_num%out_num_views][in_sinogram.get_max_tangential_pos_num()] = + (in_sinogram[view_num][(in_sinogram.get_max_tangential_pos_num()-1)])/4; + } + + + //BOUNDARY FOR MIN AND MAX VIEWS + out_sinogram[first_view][tangential_pos_num] = in_sinogram[first_view][tangential_pos_num]+ + (in_sinogram[first_view+1][tangential_pos_num] + + in_sinogram[last_view][tangential_pos_num]+ + in_sinogram[first_view][(tangential_pos_num-1)] + + in_sinogram[first_view][(tangential_pos_num+1)])/4; + + //CORNERS + out_sinogram[first_view][in_sinogram.get_min_tangential_pos_num()] = + (in_sinogram[first_view][(in_sinogram.get_max_tangential_pos_num())])/4; + out_sinogram[first_view][in_sinogram.get_max_tangential_pos_num()] = + (in_sinogram[first_view][(in_sinogram.get_min_tangential_pos_num())])/4; + + + } + } + } + static void make_non_interleaved_segment(SegmentBySinogram& out_segment, const SegmentBySinogram& in_segment) @@ -151,35 +255,75 @@ namespace detail_interpolate_projdata } } + + + + static void + transpose_make_non_interleaved_segment(SegmentBySinogram& out_segment, + const SegmentBySinogram& in_segment) + { + if (is_null_ptr(dynamic_pointer_cast(in_segment.get_proj_data_info_sptr()))) + error("make_non_interleaved_proj_data is only appropriate for non-arccorrected data"); + + for (int axial_pos_num = out_segment.get_min_axial_pos_num(); + axial_pos_num <= out_segment.get_max_axial_pos_num(); + ++axial_pos_num) + { + Sinogram out_sinogram = out_segment.get_sinogram(axial_pos_num); + transpose_make_non_interleaved_sinogram(out_sinogram, + in_segment.get_sinogram(axial_pos_num)); + out_segment.set_sinogram(out_sinogram); + } + } + static SegmentBySinogram make_non_interleaved_segment(const ProjDataInfo& non_interleaved_proj_data_info, const SegmentBySinogram& in_segment) { + SegmentBySinogram out_segment = non_interleaved_proj_data_info.get_empty_segment_by_sinogram(in_segment.get_segment_num()); make_non_interleaved_segment(out_segment, in_segment); + return out_segment; } + static SegmentBySinogram + transpose_make_non_interleaved_segment(const ProjDataInfo& compressed_proj_data_info, + const SegmentBySinogram& in_segment) + { + + SegmentBySinogram out_segment = + compressed_proj_data_info.get_empty_segment_by_sinogram(in_segment.get_segment_num()); + + transpose_make_non_interleaved_segment(out_segment, in_segment); + + return out_segment; + } } // end namespace detail_interpolate_projdata using namespace detail_interpolate_projdata; -Succeeded +Succeeded interpolate_projdata(ProjData& proj_data_out, const ProjData& proj_data_in, const BSpline::BSplineType these_types, const bool remove_interleaving, const bool use_view_offset) { + BasicCoordinate<3, BSpline::BSplineType> these_types_3; + these_types_3[1]=these_types_3[2]=these_types_3[3]=these_types; + interpolate_projdata(proj_data_out,proj_data_in,these_types_3, remove_interleaving, use_view_offset); + + return Succeeded::yes; } -Succeeded +Succeeded interpolate_projdata(ProjData& proj_data_out, const ProjData& proj_data_in, const BasicCoordinate<3, BSpline::BSplineType> & these_types, @@ -187,6 +331,8 @@ interpolate_projdata(ProjData& proj_data_out, const bool use_view_offset) { + + if (use_view_offset) warning("interpolate_projdata with use_view_offset is EXPERIMENTAL and NOT TESTED."); @@ -212,8 +358,9 @@ interpolate_projdata(ProjData& proj_data_out, BSpline::BSplinesRegularGrid<3, float, float> proj_data_interpolator(these_types); + BasicCoordinate<3, double> offset, step ; - + // find relation between out_index and in_index such that they correspond to the same physical position // out_index * m_zoom + m_offset = in_index const float in_sampling_m = proj_data_in_info.get_sampling_in_m(Bin(0,0,0,0)); @@ -255,7 +402,10 @@ interpolate_projdata(ProjData& proj_data_out, // initialise interpolator if (remove_interleaving) + { + + shared_ptr non_interleaved_proj_data_info_sptr = make_non_interleaved_proj_data_info(proj_data_in_info); @@ -263,6 +413,7 @@ interpolate_projdata(ProjData& proj_data_out, make_non_interleaved_segment(*non_interleaved_proj_data_info_sptr, proj_data_in.get_segment_by_sinogram(0)); // display(non_interleaved_segment, non_interleaved_segment.find_max(),"non-inter"); + Array<3,float> extended = extend_segment_in_views(non_interleaved_segment, 2, 2); for (int z=extended.get_min_index(); z<= extended.get_max_index(); ++z) @@ -276,6 +427,8 @@ interpolate_projdata(ProjData& proj_data_out, extended[z][y][old_max+1] = extended[z][y][old_max]; } } + std::cerr << "MAX UP "<< extended.find_max() << '\n'; + std::cerr << "MIN UP "<< extended.find_min() << '\n'; proj_data_interpolator.set_coef(extended); } else @@ -293,17 +446,271 @@ interpolate_projdata(ProjData& proj_data_out, extended[z][y][old_max+1] = extended[z][y][old_max]; } } + std::cerr << "MAX UP "<< extended.find_max() << '\n'; + std::cerr << "MIN UP "<< extended.find_min() << '\n'; proj_data_interpolator.set_coef(extended); } - // now do interpolation + // now do interpolation + SegmentBySinogram sino_3D_out = proj_data_out.get_empty_segment_by_sinogram(0) ; sample_function_on_regular_grid(sino_3D_out, proj_data_interpolator, offset, step); proj_data_out.set_segment(sino_3D_out); + if (proj_data_out.set_segment(sino_3D_out) == Succeeded::no) return Succeeded::no; return Succeeded::yes; } + + + +Succeeded +interpolate_projdata_pull(ProjData& proj_data_out, + const ProjData& proj_data_in, + const bool remove_interleaving, + const bool use_view_offset) +{ + + + SegmentBySinogram sino_3D_out = proj_data_out.get_empty_segment_by_sinogram(0) ; + if (use_view_offset) + warning("interpolate_projdata with use_view_offset is EXPERIMENTAL and NOT TESTED."); + + const ProjDataInfo & proj_data_in_info = + *proj_data_in.get_proj_data_info_sptr(); + const ProjDataInfo & proj_data_out_info = + *proj_data_out.get_proj_data_info_sptr(); + + if (typeid(proj_data_in_info) != typeid(proj_data_out_info)) + { + error("interpolate_projdata needs both projection data to be of the same type\n" + "(e.g. both arc-corrected or both not arc-corrected)"); + } + // check for the same ring radius + // This is strictly speaking only necessary for non-arccorrected data, but + // we leave it in for all cases. + if (fabs(proj_data_in_info.get_scanner_ptr()->get_inner_ring_radius() - + proj_data_out_info.get_scanner_ptr()->get_inner_ring_radius()) > 1) + { + error("interpolate_projdata needs both projection to be of a scanner with the same ring radius"); + } + + + + + + BasicCoordinate<3, double> offset, step ; + + // find relation between out_index and in_index such that they correspond to the same physical position + // out_index * m_zoom + m_offset = in_index + const float in_sampling_m = proj_data_in_info.get_sampling_in_m(Bin(0,0,0,0)); + const float out_sampling_m = proj_data_out_info.get_sampling_in_m(Bin(0,0,0,0)); + // offset in 'in' index units + offset[1] = + (proj_data_in_info.get_m(Bin(0,0,0,0)) - + proj_data_out_info.get_m(Bin(0,0,0,0))) / in_sampling_m; + step[1]= + out_sampling_m/in_sampling_m; + + const float in_sampling_phi = + (proj_data_in_info.get_phi(Bin(0,1,0,0)) - proj_data_in_info.get_phi(Bin(0,0,0,0))) / + (remove_interleaving ? 2 : 1); + + const float out_sampling_phi = + proj_data_out_info.get_phi(Bin(0,1,0,0)) - proj_data_out_info.get_phi(Bin(0,0,0,0)); + + const float out_view_offset = + use_view_offset + ? proj_data_out_info.get_scanner_ptr()->get_intrinsic_azimuthal_tilt() + : 0.F; + const float in_view_offset = + use_view_offset + ? proj_data_in_info.get_scanner_ptr()->get_intrinsic_azimuthal_tilt() + : 0.F; + offset[2] = + (proj_data_in_info.get_phi(Bin(0,0,0,0)) + in_view_offset - proj_data_out_info.get_phi(Bin(0,0,0,0)) - out_view_offset) / in_sampling_phi; + step[2] = + out_sampling_phi/in_sampling_phi; + + const float in_sampling_s = proj_data_in_info.get_sampling_in_s(Bin(0,0,0,0)); + const float out_sampling_s = proj_data_out_info.get_sampling_in_s(Bin(0,0,0,0)); + offset[3] = + (proj_data_out_info.get_s(Bin(0,0,0,0)) - + proj_data_in_info.get_s(Bin(0,0,0,0))) / in_sampling_s; + step[3]= + out_sampling_s/in_sampling_s; + + + if (remove_interleaving) + + { + shared_ptr non_interleaved_proj_data_info_sptr = make_non_interleaved_proj_data_info(proj_data_in_info); //create non interleaved projdata info + const SegmentBySinogram non_interleaved_segment = make_non_interleaved_segment(*non_interleaved_proj_data_info_sptr, proj_data_in.get_segment_by_sinogram(0)); //create non interleaved segment + + Array<3,float> extended = extend_segment_in_views(non_interleaved_segment, 2, 2); //here we are doubling the number of views + extend_tangential_position(extended); // here we are adding one tangential position before and after max and min index (for the interpolation) + extend_axial_position(extended); + sample_function_on_regular_grid_pull(sino_3D_out,extended, offset, step); //pull interpolation + proj_data_out.set_segment(sino_3D_out); //fill the output + if (proj_data_out.set_segment(sino_3D_out) == Succeeded::no) + return Succeeded::no; + } + else + { + Array<3,float> extended = extend_segment_in_views(proj_data_in.get_segment_by_sinogram(0), 2, 2); //here we are extending the number of views + extend_tangential_position(extended); // here we are adding one tangential position before and after max and min index (for the interpolation) + extend_axial_position(extended); + sample_function_on_regular_grid_pull(sino_3D_out,extended, offset, step); //pull interpolation + proj_data_out.set_segment(sino_3D_out); //fill the output + if (proj_data_out.set_segment(sino_3D_out) == Succeeded::no) + return Succeeded::no; + + } + + return Succeeded::yes; +} + + +Succeeded +interpolate_projdata_push(ProjData& proj_data_out, + const ProjData& proj_data_in, + const bool remove_interleaving, + const bool use_view_offset) +{ + + + //SegmentBySinogram sino_3D_out = proj_data_out.get_empty_segment_by_sinogram(0) ; + + if (use_view_offset) + warning("interpolate_projdata with use_view_offset is EXPERIMENTAL and NOT TESTED."); + + const ProjDataInfo & proj_data_in_info = + *proj_data_in.get_proj_data_info_sptr(); + const ProjDataInfo & proj_data_out_info = + *proj_data_out.get_proj_data_info_sptr(); + + if (typeid(proj_data_in_info) != typeid(proj_data_out_info)) + { + error("interpolate_projdata needs both projection data to be of the same type\n" + "(e.g. both arc-corrected or both not arc-corrected)"); + } + + if (fabs(proj_data_in_info.get_scanner_ptr()->get_inner_ring_radius() - + proj_data_out_info.get_scanner_ptr()->get_inner_ring_radius()) > 1) + { + error("interpolate_projdata needs both projection to be of a scanner with the same ring radius"); + } + + + + BasicCoordinate<3, double> offset, step ; + + // find relation between out_index and in_index such that they correspond to the same physical position + // out_index * m_zoom + m_offset = in_index + const float in_sampling_m = proj_data_in_info.get_sampling_in_m(Bin(0,0,0,0)); + const float out_sampling_m = proj_data_out_info.get_sampling_in_m(Bin(0,0,0,0)); + // offset in 'in' index units + offset[1] = + (proj_data_out_info.get_m(Bin(0,0,0,0)) - + proj_data_in_info.get_m(Bin(0,0,0,0))) / out_sampling_m; //divide by sampling: conversion from mm to voxel units + step[1]= + in_sampling_m/out_sampling_m; + + const float out_sampling_phi = + (proj_data_out_info.get_phi(Bin(0,1,0,0)) - proj_data_out_info.get_phi(Bin(0,0,0,0))) / + (remove_interleaving ? 2 : 1); + + const float in_sampling_phi = + proj_data_in_info.get_phi(Bin(0,1,0,0)) - proj_data_in_info.get_phi(Bin(0,0,0,0)); + + const float out_view_offset = + use_view_offset + ? proj_data_out_info.get_scanner_ptr()->get_intrinsic_azimuthal_tilt() + : 0.F; + const float in_view_offset = + use_view_offset + ? proj_data_in_info.get_scanner_ptr()->get_intrinsic_azimuthal_tilt() + : 0.F; + offset[2] = + (proj_data_out_info.get_phi(Bin(0,0,0,0)) + out_view_offset - proj_data_in_info.get_phi(Bin(0,0,0,0)) - in_view_offset) / out_sampling_phi; + + step[2] = + in_sampling_phi/out_sampling_phi; + + const float out_sampling_s = proj_data_out_info.get_sampling_in_s(Bin(0,0,0,0)); + const float in_sampling_s = proj_data_in_info.get_sampling_in_s(Bin(0,0,0,0)); + offset[3] = + (proj_data_in_info.get_s(Bin(0,0,0,0)) - + proj_data_out_info.get_s(Bin(0,0,0,0))) / out_sampling_s; + step[3]= + in_sampling_s/out_sampling_s; + + if (remove_interleaving) + + { + // =================== RESIZE 'EXTENDED' OUTPUT ========================= + shared_ptr non_interleaved_proj_data_info_sptr = make_non_interleaved_proj_data_info(proj_data_out_info); + const SegmentBySinogram non_interleaved_segment = make_non_interleaved_segment(*non_interleaved_proj_data_info_sptr, proj_data_out.get_segment_by_sinogram(0)); + Array<3,float> extended = extend_segment_in_views(non_interleaved_segment,2,2); //extend the views + extend_tangential_position(extended); //extend the tangential positions + extend_axial_position(extended); + // ===========================PUSH ================================== + + SegmentBySinogram sino_3D_in = proj_data_in.get_segment_by_sinogram(0); //TODO: check if we need a for loop over segments + sample_function_on_regular_grid_push(extended,sino_3D_in, offset, step); //here the output of the push is 'extended' + + // =================== TRANSPOSE EXTENDED =========================== + + transpose_extend_axial_position(extended); + transpose_extend_tangential_position(extended); //reduce tangential positions + SegmentBySinogram extended_segment_sino(extended, non_interleaved_proj_data_info_sptr, 0); //create SegmentBySinogram with extended + Array<3,float> out = transpose_extend_segment_in_views(extended_segment_sino,2,2); // here we do the tranpose : extended -> sino_out + // =================TRANSPOSE REMOVE INTERLEAVING ==================== + + SegmentBySinogram compressed_output(out, non_interleaved_proj_data_info_sptr, 0); + + shared_ptr new_proj_data_info_sptr = proj_data_out_info.create_shared_clone(); + SegmentBySinogram transpose_interleaved_segment = transpose_make_non_interleaved_segment(*new_proj_data_info_sptr, compressed_output); + + + // ======================== CREATE OUTPUT ============================= + proj_data_out.set_segment(transpose_interleaved_segment); + if (proj_data_out.set_segment(transpose_interleaved_segment) == Succeeded::no) + return Succeeded::no; + } + else + { + // =================== RESIZE 'EXTENDED' OUTPUT ========================= + + Array<3,float> extended = extend_segment_in_views(proj_data_out.get_segment_by_sinogram(0), 2, 2); //extend the views + extend_tangential_position(extended); //extend the tangential positions + extend_axial_position(extended); + // ===========================PUSH ================================== + + SegmentBySinogram sino_3D_in = proj_data_in.get_segment_by_sinogram(0); + sample_function_on_regular_grid_push(extended,sino_3D_in, offset, step); //here the output of the push is 'extended' + // =================== TRANSPOSE EXTENDED =========================== + + transpose_extend_axial_position(extended); + transpose_extend_tangential_position(extended); + shared_ptr extended_proj_data_info_sptr(proj_data_out_info.clone()); //create extended projdata inf + SegmentBySinogram extended_segment_sino(extended, extended_proj_data_info_sptr, 0); //create SegmentBySinogram with extended + Array<3,float> out = transpose_extend_segment_in_views(extended_segment_sino,2,2); // here we do the tranpose : extended -> sino_out + + // ======================== CREATE OUTPUT ============================= + + SegmentBySinogram compressed_output(out, extended_proj_data_info_sptr, 0); + + proj_data_out.set_segment(compressed_output); + if (proj_data_out.set_segment(compressed_output) == Succeeded::no) + + return Succeeded::no; + + } + + return Succeeded::yes; +} + END_NAMESPACE_STIR diff --git a/src/buildblock/inverse_SSRB.cxx b/src/buildblock/inverse_SSRB.cxx index 943c0441ce..9daeec470c 100644 --- a/src/buildblock/inverse_SSRB.cxx +++ b/src/buildblock/inverse_SSRB.cxx @@ -58,10 +58,11 @@ inverse_SSRB(ProjData& proj_data_4D, Sinogram sino_4D = proj_data_4D. get_empty_sinogram(proj_data_4D.get_min_axial_pos_num(0) , 0); + sino_4D.fill(0); Sinogram sino_3D = proj_data_3D. get_empty_sinogram(proj_data_3D.get_min_axial_pos_num(0) , 0); - + sino_3D.fill(0); for (int out_segment_num = proj_data_4D.get_min_segment_num(); out_segment_num <= proj_data_4D.get_max_segment_num(); ++out_segment_num) @@ -70,6 +71,10 @@ inverse_SSRB(ProjData& proj_data_4D, out_ax_pos_num <= proj_data_4D.get_max_axial_pos_num(out_segment_num); ++out_ax_pos_num ) { + + // if((out_segment_num+out_ax_pos_num )>proj_data_4D.get_max_axial_pos_num(out_segment_num)||(out_segment_num+out_ax_pos_num) @@ -111,6 +116,84 @@ inverse_SSRB(ProjData& proj_data_4D, out_segment_num, out_ax_pos_num); } } + // Sinogram test = proj_data_4D.get_sinogram(1,0); + //std::cout << "expeted:1\n" << proj_data_4D.get_sinogram(1,0)[0][0]; + //std::cout << "expected: ?\n" << proj_data_4D.get_sinogram(1,-2)[0][0]; + // std::cout << "expected: /=1\n" << proj_data_4D.get_sinogram(1,2)[0][0]; return Succeeded::yes; } + + +Succeeded +transpose_inverse_SSRB(ProjData& proj_data_3D, + const ProjData& proj_data_4D) +{ + //construct projdata info + const shared_ptr proj_data_3D_info_sptr = + dynamic_pointer_cast + (proj_data_3D.get_proj_data_info_sptr()); + const shared_ptr proj_data_4D_info_sptr = + dynamic_pointer_cast + (proj_data_4D.get_proj_data_info_sptr()); + + // keep sinograms out of the loop to avoid reallocations. initialise to something because there's no default constructor + Sinogram sino_3D = proj_data_3D.get_empty_sinogram(proj_data_3D.get_min_axial_pos_num(0) , 0); + Sinogram sino_3D2 = proj_data_3D.get_empty_sinogram(proj_data_3D.get_min_axial_pos_num(0) , 0); + Sinogram sino_4D = proj_data_4D.get_empty_sinogram(proj_data_4D.get_min_axial_pos_num(0) , 0); + + sino_3D.fill(0); + sino_3D2.fill(0); + sino_4D.fill(0); + for (int out_ax_pos_num = proj_data_3D.get_min_axial_pos_num(0); out_ax_pos_num <= proj_data_3D.get_max_axial_pos_num(0); ++out_ax_pos_num ) + { + sino_3D = proj_data_3D.get_empty_sinogram(out_ax_pos_num, 0); + const float out_m = proj_data_3D_info_sptr->get_m(Bin(0, 0, out_ax_pos_num, 0)); + const float out_m_next = out_ax_pos_num == proj_data_3D.get_max_axial_pos_num(0) ? + -1000000.F : proj_data_3D_info_sptr->get_m(Bin(0, 0, out_ax_pos_num+1, 0)); + const float out_m_prev = out_ax_pos_num == proj_data_3D.get_min_axial_pos_num(0) ? + -1000000.F : proj_data_3D_info_sptr->get_m(Bin(0, 0, out_ax_pos_num-1, 0)); + + for (int in_segment_num = proj_data_4D.get_min_segment_num(); in_segment_num <= proj_data_4D.get_max_segment_num(); ++in_segment_num) + { + for (int in_ax_pos_num = proj_data_4D.get_min_axial_pos_num(in_segment_num);in_ax_pos_num <= proj_data_4D.get_max_axial_pos_num(in_segment_num); ++in_ax_pos_num ) + { + + //if((in_segment_num+out_ax_pos_num )>proj_data_4D.get_max_axial_pos_num(in_segment_num)||(in_segment_num+out_ax_pos_num)get_m(Bin(in_segment_num, 0, in_ax_pos_num, 0)); + + if (fabs(out_m - in_m) < 1E-2) + { + sino_4D = proj_data_4D.get_sinogram(in_ax_pos_num,in_segment_num); + sino_3D += sino_4D; + } + + + if ((fabs(in_m - .5F*(out_m + out_m_next)) < 1E-2) || (fabs(in_m - .5F*(out_m + out_m_prev)) < 1E-2)) + { + // for (int i = out_ax_pos_num + 1; i <= proj_data_3D.get_max_axial_pos_num(0); ++i) + //{ + // sino_4D *= .5F; + sino_4D = proj_data_4D.get_sinogram(in_ax_pos_num,in_segment_num); + sino_3D += sino_4D/2; + // sino_3D2 = proj_data_3D.get_empty_sinogram(i, 0); + // sino_3D2 += sino_4D; + + } + + + } + + } + + proj_data_3D.set_sinogram(sino_3D); + //proj_data_3D.set_sinogram(sino_3D2); + } + + + + return Succeeded::yes; +} + END_NAMESPACE_STIR diff --git a/src/buildblock/zoom.cxx b/src/buildblock/zoom.cxx index 4eaee9ce94..dab975382d 100644 --- a/src/buildblock/zoom.cxx +++ b/src/buildblock/zoom.cxx @@ -3,23 +3,18 @@ Copyright (C) 2000- 2011, Hammersmith Imanet Ltd Copyright (C) 2018-2019, University College London This file is part of STIR. - SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license See STIR/LICENSE.txt for details */ /*! - \file + \file \ingroup buildblock - \brief Implementations of the stir::zoom functions - \author Kris Thielemans \author Claire Labbe \author PARAPET project \author Ludovica Brusaferri - - */ /* Modification history: - First versions by CL and KT (sinogram version based on C code by Matthias Egger (using linear interpolation). @@ -28,13 +23,13 @@ - KT introduced 3D zooming for images. - KT converted to new design */ - + #include "stir/interpolate.h" #include "stir/zoom.h" #include "stir/DataProcessor.h" #include "stir/DiscretisedDensity.h" -#include "stir/VoxelsOnCartesianGrid.h" -#include "stir/PixelsOnCartesianGrid.h" +#include "stir/VoxelsOnCartesianGrid.h" +#include "stir/PixelsOnCartesianGrid.h" #include "stir/Viewgram.h" #include "stir/RelatedViewgrams.h" #include "stir/ProjDataInfoCylindricalArcCorr.h" @@ -49,42 +44,42 @@ START_NAMESPACE_STIR // also they need to be converted to the new design // as we don't need them at the moment, I can't be bothered... #if 0 -void zoom_segment (SegmentByView& segment, - const float zoom, const float azimuthal_angle_sampling, const float y_offset_in_mm, +void zoom_segment (SegmentByView& segment, + const float zoom, const float azimuthal_angle_sampling, const float y_offset_in_mm, const int new_size, const float azimuthal_angle_sampling) { - // First check if there is anything to do at all, + // First check if there is anything to do at all, // if not, return original segment. if (new_size == segment.get_num_bins() && - zoom == 1.0 && x_offset_in_mm == 0.0 && y_offset_in_mm == 0.0) + zoom == 1.0 && x_offset_in_mm == 0.0 && y_offset_in_mm == 0.0) return; // KT 17/01/2000 use local copy of scan_info instead of segment.scan_info ScanInfo scan_info = segment.scan_info; scan_info.set_num_bins(new_size); scan_info.set_bin_size(segment.scan_info.get_bin_size() / zoom); - + const int minsize = -new_size/2; const int maxsize = minsize+new_size-1; - + // TODO replace by get_empty_segment or so - SegmentByView + SegmentByView out_segment(Tensor3D(segment.get_min_view(), segment.get_max_view(), - segment.get_min_axial_pos_num(), segment.get_max_axial_pos_num(), - minsize, maxsize), - scan_info, - segment.get_segment_num(), - segment.get_min_axial_pos_num_difference(), - segment.get_max_ring_difference()); - - for (int view = segment.get_min_view(); view <= segment.get_max_view(); view++) + segment.get_min_axial_pos_num(), segment.get_max_axial_pos_num(), + minsize, maxsize), + scan_info, + segment.get_segment_num(), + segment.get_min_axial_pos_num_difference(), + segment.get_max_ring_difference()); + + for (int view = segment.get_min_view(); view <= segment.get_max_view(); view++) { Viewgram viewgram = segment.get_viewgram(view); zoom_viewgram(viewgram, - zoom, x_offset_in_mm, y_offset_in_mm, - new_size, azimuthal_angle_sampling); + zoom, x_offset_in_mm, y_offset_in_mm, + new_size, azimuthal_angle_sampling); out_segment.set_viewgram(viewgram); } @@ -95,16 +90,16 @@ void zoom_segment (SegmentByView& segment, #endif void -zoom_viewgrams (RelatedViewgrams& in_viewgrams, - const float zoom, - const int min_tang_pos_num, const int max_tang_pos_num, - const float x_offset_in_mm, const float y_offset_in_mm) +zoom_viewgrams (RelatedViewgrams& in_viewgrams, + const float zoom, + const int min_tang_pos_num, const int max_tang_pos_num, + const float x_offset_in_mm, const float y_offset_in_mm) { if (min_tang_pos_num == in_viewgrams.get_min_tangential_pos_num() && max_tang_pos_num == in_viewgrams.get_max_tangential_pos_num() && - zoom == 1.0 && x_offset_in_mm == 0.0 && y_offset_in_mm == 0.0) + zoom == 1.0 && x_offset_in_mm == 0.0 && y_offset_in_mm == 0.0) return; - + shared_ptr new_proj_data_info_sptr(in_viewgrams.get_proj_data_info_sptr()->clone()); ProjDataInfoCylindricalArcCorr* new_proj_data_info_arccorr_sptr = @@ -120,37 +115,37 @@ zoom_viewgrams (RelatedViewgrams& in_viewgrams, set_tangential_sampling(new_proj_data_info_arccorr_sptr-> get_tangential_sampling() / zoom); - shared_ptr + shared_ptr symmetries_sptr(in_viewgrams.get_symmetries_ptr()->clone()); RelatedViewgrams out_viewgrams = new_proj_data_info_arccorr_sptr-> get_empty_related_viewgrams(in_viewgrams.get_basic_view_segment_num(), - symmetries_sptr); + symmetries_sptr); { RelatedViewgrams::iterator out_iter = out_viewgrams.begin(); RelatedViewgrams::const_iterator in_iter = in_viewgrams.begin(); for (; out_iter != out_viewgrams.end(); ++out_iter, ++in_iter) zoom_viewgram(*out_iter, *in_iter, - x_offset_in_mm, y_offset_in_mm); + x_offset_in_mm, y_offset_in_mm); } in_viewgrams = out_viewgrams; } void -zoom_viewgram (Viewgram& in_view, - const float zoom, - const int min_tang_pos_num, const int max_tang_pos_num, - const float x_offset_in_mm, const float y_offset_in_mm) -{ +zoom_viewgram (Viewgram& in_view, + const float zoom, + const int min_tang_pos_num, const int max_tang_pos_num, + const float x_offset_in_mm, const float y_offset_in_mm) +{ if (min_tang_pos_num == in_view.get_min_tangential_pos_num() && max_tang_pos_num == in_view.get_max_tangential_pos_num() && - zoom == 1.0 && x_offset_in_mm == 0.0 && y_offset_in_mm == 0.0) + zoom == 1.0 && x_offset_in_mm == 0.0 && y_offset_in_mm == 0.0) return; - + shared_ptr new_proj_data_info_sptr(in_view.get_proj_data_info_sptr()->clone()); ProjDataInfoCylindricalArcCorr* new_proj_data_info_arccorr_sptr = @@ -169,20 +164,20 @@ zoom_viewgram (Viewgram& in_view, Viewgram out_view = new_proj_data_info_arccorr_sptr-> get_empty_viewgram( - in_view.get_view_num(), - in_view.get_segment_num()); + in_view.get_view_num(), + in_view.get_segment_num()); - zoom_viewgram(out_view, in_view, - x_offset_in_mm, y_offset_in_mm); + zoom_viewgram(out_view, in_view, + x_offset_in_mm, y_offset_in_mm); in_view = out_view; } void -zoom_viewgram (Viewgram& out_view, - const Viewgram& in_view, - const float x_offset_in_mm, const float y_offset_in_mm) -{ +zoom_viewgram (Viewgram& out_view, + const Viewgram& in_view, + const float x_offset_in_mm, const float y_offset_in_mm) +{ // minimal checks on compatibility assert(in_view.get_proj_data_info_sptr()->get_num_views() == out_view.get_proj_data_info_sptr()->get_num_views()); @@ -213,21 +208,21 @@ zoom_viewgram (Viewgram& out_view, if (out_view.get_min_tangential_pos_num() == in_view.get_min_tangential_pos_num() && out_view.get_max_tangential_pos_num() == in_view.get_max_tangential_pos_num() && - zoom == 1.0F && x_offset_in_mm == 0.0F && y_offset_in_mm == 0.0F) + zoom == 1.0F && x_offset_in_mm == 0.0F && y_offset_in_mm == 0.0F) return; - + const float phi = in_proj_data_info_arccorr_sptr-> get_phi(Bin(in_view.get_segment_num(), in_view.get_view_num(), 0,0)); // compute offset in tangential_sampling_in units - const float offset = + const float offset = (x_offset_in_mm*cos(phi) +y_offset_in_mm*sin(phi))/ in_bin_size; for (int axial_pos_num= out_view.get_min_axial_pos_num(); axial_pos_num <= out_view.get_max_axial_pos_num(); ++axial_pos_num) { overlap_interpolate(out_view[axial_pos_num], in_view[axial_pos_num], zoom, offset); - } + } } @@ -236,9 +231,9 @@ zoom_viewgram (Viewgram& out_view, static VoxelsOnCartesianGrid construct_new_image_from_zoom_parameters(const VoxelsOnCartesianGrid &image, - const CartesianCoordinate3D& zooms, - const CartesianCoordinate3D& offsets_in_mm, - const BasicCoordinate<3,int>& new_sizes_arg) + const CartesianCoordinate3D& zooms, + const CartesianCoordinate3D& offsets_in_mm, + const BasicCoordinate<3,int>& new_sizes_arg) { CartesianCoordinate3D new_sizes = new_sizes_arg; assert(new_sizes.x()>=0); @@ -250,14 +245,14 @@ construct_new_image_from_zoom_parameters(const VoxelsOnCartesianGrid &ima // first set origin to 0 CartesianCoordinate3D origin(0.F,0.F,0.F); - - VoxelsOnCartesianGrid + + VoxelsOnCartesianGrid new_image(image.get_exam_info_sptr(), IndexRange3D(0, new_sizes.z()-1, - -new_sizes.y()/2, -new_sizes.y()/2+new_sizes.y()-1, - -new_sizes.x()/2, -new_sizes.x()/2+new_sizes.x()-1), - origin, - voxel_size); + -new_sizes.y()/2, -new_sizes.y()/2+new_sizes.y()-1, + -new_sizes.x()/2, -new_sizes.x()/2+new_sizes.x()-1), + origin, + voxel_size); // find coordinates of middle of images BasicCoordinate<3,int> min_indices, max_indices; @@ -288,9 +283,9 @@ construct_new_image_from_zoom_parameters(const VoxelsOnCartesianGrid &ima } void zoom_image_in_place(VoxelsOnCartesianGrid &image, - const float zoom, - const float x_offset_in_mm, const float y_offset_in_mm, - const int new_size, + const float zoom, + const float x_offset_in_mm, const float y_offset_in_mm, + const int new_size, const ZoomOptions zoom_options) { VoxelsOnCartesianGrid new_image = @@ -301,14 +296,14 @@ zoom_image_in_place(VoxelsOnCartesianGrid &image, VoxelsOnCartesianGrid zoom_image(const VoxelsOnCartesianGrid &image, const float zoom, - const float x_offset_in_mm, const float y_offset_in_mm, - const int new_size, + const float x_offset_in_mm, const float y_offset_in_mm, + const int new_size, const ZoomOptions zoom_options) { assert(new_size>=0); - if(zoom==1 && x_offset_in_mm==0 && y_offset_in_mm==0 && new_size== image.get_x_size()) + if(zoom==1 && x_offset_in_mm==0 && y_offset_in_mm==0 && new_size== image.get_x_size()) return image; - + const CartesianCoordinate3D zooms(1,zoom,zoom); const CartesianCoordinate3D offsets_in_mm(0.F, y_offset_in_mm, x_offset_in_mm); const BasicCoordinate<3,int> new_sizes = @@ -316,18 +311,18 @@ zoom_image(const VoxelsOnCartesianGrid &image, VoxelsOnCartesianGrid new_image = construct_new_image_from_zoom_parameters(image, - zooms, - offsets_in_mm, - new_sizes); + zooms, + offsets_in_mm, + new_sizes); - PixelsOnCartesianGrid + PixelsOnCartesianGrid new_image2D = new_image.get_plane(new_image.get_min_z()); for (int plane = image.get_min_z(); plane <= image.get_max_z(); plane++) { zoom_image(new_image2D, image.get_plane(plane), zoom_options); new_image.set_plane(new_image2D, plane); } - + assert(norm(new_image.get_voxel_size() - image.get_voxel_size()/zooms)<1); return new_image; @@ -335,9 +330,9 @@ zoom_image(const VoxelsOnCartesianGrid &image, void zoom_image_in_place(VoxelsOnCartesianGrid &image, - const CartesianCoordinate3D& zooms, - const CartesianCoordinate3D& offsets_in_mm, - const BasicCoordinate<3,int>& new_sizes, + const CartesianCoordinate3D& zooms, + const CartesianCoordinate3D& offsets_in_mm, + const BasicCoordinate<3,int>& new_sizes, const ZoomOptions zoom_options) { const VoxelsOnCartesianGrid new_image = @@ -347,23 +342,23 @@ zoom_image_in_place(VoxelsOnCartesianGrid &image, VoxelsOnCartesianGrid zoom_image(const VoxelsOnCartesianGrid &image, - const CartesianCoordinate3D& zooms, - const CartesianCoordinate3D& offsets_in_mm, - const BasicCoordinate<3,int>& new_sizes, + const CartesianCoordinate3D& zooms, + const CartesianCoordinate3D& offsets_in_mm, + const BasicCoordinate<3,int>& new_sizes, const ZoomOptions zoom_options) { VoxelsOnCartesianGrid new_image = construct_new_image_from_zoom_parameters(image, - zooms, - offsets_in_mm, - new_sizes); + zooms, + offsets_in_mm, + new_sizes); zoom_image(new_image, image, zoom_options); return new_image; } -void -zoom_image(VoxelsOnCartesianGrid &image_out, +void +zoom_image(VoxelsOnCartesianGrid &image_out, const VoxelsOnCartesianGrid &image_in, const ZoomOptions zoom_options) { @@ -371,47 +366,44 @@ zoom_image(VoxelsOnCartesianGrid &image_out, /* interpolation routine uses the following relation: x_in_index = x_out_index/zoom + offset - compare to 'physical' coordinates x_phys = (x_index) * voxel_size.x + origin.x - + as x_in_phys == x_out_phys, we find (x_in_index)* voxel_size_in.x + origin_in.x == (x_out_index )* voxel_size_out.x + origin_out.x <=> - x_in_index = (x_out_index * voxel_size_out.x - + origin_out.x - origin_in.x) - / voxel_size_in.x - + x_in_index = (x_out_index * voxel_size_out.x + + origin_out.x - origin_in.x) + / voxel_size_in.x so, zoom= voxel_size_in.x/ voxel_size_out.x offset = (origin_out.x - origin_in.x)/ voxel_size_in.x - */ // check relation between indices and physical coordinates { const BasicCoordinate<3,int> indices = make_coordinate(1,2,3); if (norm(image_in.get_physical_coordinates_for_indices(indices) - - (image_in.get_voxel_size() * BasicCoordinate<3,float>(indices) + image_in.get_origin()) - ) > 2.F) + - (image_in.get_voxel_size() * BasicCoordinate<3,float>(indices) + image_in.get_origin()) + ) > 2.F) error("zoom_image is confused about the relation between indices and physical coordinates"); } - const float zoom_x = + const float zoom_x = image_in.get_voxel_size().x() / image_out.get_voxel_size().x(); - const float zoom_y = + const float zoom_y = image_in.get_voxel_size().y() / image_out.get_voxel_size().y(); - const float zoom_z = + const float zoom_z = image_in.get_voxel_size().z() / image_out.get_voxel_size().z(); - const float x_offset = - (image_out.get_origin().x() - image_in.get_origin().x()) + const float x_offset = + (image_out.get_origin().x() - image_in.get_origin().x()) / image_in.get_voxel_size().x(); - const float y_offset = - (image_out.get_origin().y() - image_in.get_origin().y()) + const float y_offset = + (image_out.get_origin().y() - image_in.get_origin().y()) / image_in.get_voxel_size().y(); - const float z_offset = - (image_out.get_origin().z() - image_in.get_origin().z()) + const float z_offset = + (image_out.get_origin().z() - image_in.get_origin().z()) / image_in.get_voxel_size().z(); - + if(zoom_x==1.0F && zoom_y==1.0F && zoom_z==1.0F && x_offset == 0.F && y_offset == 0.F && z_offset == 0.F && image_in.get_index_range() == image_out.get_index_range() @@ -422,18 +414,18 @@ zoom_image(VoxelsOnCartesianGrid &image_out, } // TODO creating a lot of new images here... - Array<3,float> + Array<3,float> temp(IndexRange3D(image_in.get_min_z(), image_in.get_max_z(), - image_in.get_min_y(), image_in.get_max_y(), - image_out.get_min_x(), image_out.get_max_x())); + image_in.get_min_y(), image_in.get_max_y(), + image_out.get_min_x(), image_out.get_max_x())); for (int z=image_in.get_min_z(); z<=image_in.get_max_z(); z++) for (int y=image_in.get_min_y(); y<=image_in.get_max_y(); y++) overlap_interpolate(temp[z][y], image_in[z][y], zoom_x, x_offset); Array<3,float> temp2(IndexRange3D(image_in.get_min_z(), image_in.get_max_z(), - image_out.get_min_y(), image_out.get_max_y(), - image_out.get_min_x(), image_out.get_max_x())); + image_out.get_min_y(), image_out.get_max_y(), + image_out.get_min_x(), image_out.get_max_x())); for (int z=image_in.get_min_z(); z<=image_in.get_max_z(); z++) overlap_interpolate(temp2[z], temp[z], zoom_y, y_offset); @@ -472,7 +464,7 @@ zoom_image(VoxelsOnCartesianGrid &image_out, } void -zoom_image(PixelsOnCartesianGrid &image2D_out, +zoom_image(PixelsOnCartesianGrid &image2D_out, const PixelsOnCartesianGrid &image2D_in, const ZoomOptions zoom_options) { @@ -481,19 +473,19 @@ zoom_image(PixelsOnCartesianGrid &image2D_out, see above for how to find zoom and offsets */ - const float zoom_x = + const float zoom_x = image2D_in.get_pixel_size().x() / image2D_out.get_pixel_size().x(); - const float zoom_y = + const float zoom_y = image2D_in.get_pixel_size().y() / image2D_out.get_pixel_size().y(); - const float x_offset = - ((image2D_out.get_origin().x() - image2D_in.get_origin().x()) + const float x_offset = + ((image2D_out.get_origin().x() - image2D_in.get_origin().x()) / image2D_in.get_pixel_size().x() ); - const float y_offset = - ((image2D_out.get_origin().y() - image2D_in.get_origin().y()) + const float y_offset = + ((image2D_out.get_origin().y() - image2D_in.get_origin().y()) / image2D_in.get_pixel_size().y() ); - + if(zoom_x==1.0F && zoom_y==1.0F && x_offset == 0.F && y_offset == 0.F && image2D_in.get_index_range() == image2D_out.get_index_range() @@ -503,14 +495,14 @@ zoom_image(PixelsOnCartesianGrid &image2D_out, return; } - Array<2,float> + Array<2,float> temp(IndexRange2D(image2D_in.get_min_y(), image2D_in.get_max_y(), - image2D_out.get_min_x(), image2D_out.get_max_x())); + image2D_out.get_min_x(), image2D_out.get_max_x())); for (int y=image2D_in.get_min_y(); y<=image2D_in.get_max_y(); y++) overlap_interpolate(temp[y], image2D_in[y], zoom_x, x_offset); - overlap_interpolate(image2D_out, temp, zoom_y, y_offset); + overlap_interpolate(image2D_out, temp, zoom_y, y_offset); float scale_image = 1.F; @@ -541,4 +533,125 @@ zoom_image(PixelsOnCartesianGrid &image2D_out, } +void +transpose_zoom_image(VoxelsOnCartesianGrid &image_out, + const VoxelsOnCartesianGrid &image_in, + const ZoomOptions zoom_options) +{ + image_out.set_exam_info(image_in.get_exam_info()); + /* + interpolation routine uses the following relation: + x_in_index = x_out_index/zoom + offset + compare to 'physical' coordinates + x_phys = (x_index) * voxel_size.x + origin.x + as x_in_phys == x_out_phys, we find + (x_in_index)* voxel_size_in.x + origin_in.x == + (x_out_index )* voxel_size_out.x + origin_out.x + <=> + x_in_index = (x_out_index * voxel_size_out.x + + origin_out.x - origin_in.x) + / voxel_size_in.x + so, zoom= voxel_size_in.x/ voxel_size_out.x + offset = (origin_out.x - origin_in.x)/ voxel_size_in.x + */ + // check relation between indices and physical coordinates + { + const BasicCoordinate<3,int> indices = make_coordinate(1,2,3); + if (norm(image_in.get_physical_coordinates_for_indices(indices) + - (image_in.get_voxel_size() * BasicCoordinate<3,float>(indices) + image_in.get_origin()) + ) > 2.F) + error("zoom_image is confused about the relation between indices and physical coordinates"); + } + const float zoom_x = + image_in.get_voxel_size().x() / image_out.get_voxel_size().x(); + const float zoom_y = + image_in.get_voxel_size().y() / image_out.get_voxel_size().y(); + const float zoom_z = + image_in.get_voxel_size().z() / image_out.get_voxel_size().z(); + const float x_offset = + (image_out.get_origin().x() - image_in.get_origin().x()) + / image_in.get_voxel_size().x(); + const float y_offset = + (image_out.get_origin().y() - image_in.get_origin().y()) + / image_in.get_voxel_size().y(); + const float z_offset = + (image_out.get_origin().z() - image_in.get_origin().z()) + / image_in.get_voxel_size().z(); + + + if(zoom_x==1.0F && zoom_y==1.0F && zoom_z==1.0F && + x_offset == 0.F && y_offset == 0.F && z_offset == 0.F && + image_in.get_index_range() == image_out.get_index_range() + ) + { + image_out = image_in; + return; + } + + // TODO creating a lot of new images here... + Array<3,float> + temp(IndexRange3D(image_in.get_min_z(), image_in.get_max_z(), + image_in.get_min_y(), image_in.get_max_y(), + image_out.get_min_x(), image_out.get_max_x())); + + for (int z=image_in.get_min_z(); z<=image_in.get_max_z(); z++) + for (int y=image_in.get_min_y(); y<=image_in.get_max_y(); y++) + overlap_interpolate(temp[z][y], image_in[z][y], zoom_x, x_offset); + + Array<3,float> temp2(IndexRange3D(image_in.get_min_z(), image_in.get_max_z(), + image_out.get_min_y(), image_out.get_max_y(), + image_out.get_min_x(), image_out.get_max_x())); + + for (int z=image_in.get_min_z(); z<=image_in.get_max_z(); z++) + overlap_interpolate(temp2[z], temp[z], zoom_y, y_offset); + + temp.recycle(); + + overlap_interpolate(image_out, temp2, zoom_z, z_offset); + + float scale_image = 1.F; + + switch (zoom_options.get_scaling_option()) + { + case ZoomOptions::preserve_values: + { + return; + } + + case ZoomOptions::preserve_projections: + + { + scale_image = zoom_x; + break; + } + + case ZoomOptions::preserve_sum: + { + scale_image = zoom_x*zoom_y*zoom_z; + break; + } + + } + + if (scale_image != 1.F) + image_out*= scale_image; + +} + +void +zoom_image_swig(VoxelsOnCartesianGrid &image_out, + const VoxelsOnCartesianGrid &image_in, const int zoom_option) +{ + ZoomOptions::Scaling zo = ZoomOptions::Scaling(zoom_option); + zoom_image(image_out, image_in,zo); +} + +void +transpose_zoom_image_swig(VoxelsOnCartesianGrid &image_out, + const VoxelsOnCartesianGrid &image_in, const int zoom_option) +{ + ZoomOptions::Scaling zo = ZoomOptions::Scaling(zoom_option); + transpose_zoom_image(image_out, image_in,zo); +} + END_NAMESPACE_STIR diff --git a/src/experimental/buildblock/local_buildblock_registries.cxx b/src/experimental/buildblock/local_buildblock_registries.cxx index bd02a654b7..20725951ad 100644 --- a/src/experimental/buildblock/local_buildblock_registries.cxx +++ b/src/experimental/buildblock/local_buildblock_registries.cxx @@ -3,6 +3,7 @@ See STIR/LICENSE.txt for detail $, IRSL See STIR/LICENSE.txt for details */ + #include "stir_experimental/SeparableGaussianImageFilter.h" #ifdef SANIDA #include "stir_experimental/DAVImageFilter3D.h" diff --git a/src/include/stir/Bin.h b/src/include/stir/Bin.h index 989e3d7a33..e0746be5fd 100644 --- a/src/include/stir/Bin.h +++ b/src/include/stir/Bin.h @@ -46,8 +46,11 @@ class Bin inline Bin(); //! A constructor : constructs a bin with value (defaulting to 0) - inline Bin(int segment_num,int view_num, int axial_pos_num, - int tangential_pos_num,float bin_value=0); + //inline Bin(int segment_num,int view_num, int axial_pos_num, + // int tangential_pos_num,float bin_value=0); +// + inline Bin(int segment_num,int view_num, int axial_pos_num, + int tangential_pos_num, float bin_value = 0); //!get axial position number inline int axial_pos_num()const; @@ -59,14 +62,19 @@ class Bin inline int view_num() const; //! get time-frame number (1-based) inline int time_frame_num() const; - + //! get first_energy_window_number + inline int first_energy_window_num() const; + //! get second_energy_window_number + inline int second_energy_window_num() const; + inline int& axial_pos_num(); inline int& segment_num(); inline int& tangential_pos_num(); inline int& view_num(); inline int& time_frame_num(); - - + inline int& first_energy_window_num(); + inline int& second_energy_window_num(); + //! get an empty copy inline Bin get_empty_copy() const; @@ -74,6 +82,16 @@ class Bin inline float get_bin_value()const; //! set the value to be back projected inline void set_bin_value( float v ); + + //! get the first energy window number + inline float get_first_energy_window_num()const; + //! set the first energy window number + inline void set_first_energy_window_num( int v ); + //! get the second energy window number + inline float get_second_energy_window_num()const; + //! set the second energy window number + inline void set_second_energy_window_num( int v ); + //! accumulate voxel's contribution during forward projection inline Bin& operator+=(const float dx); @@ -96,6 +114,8 @@ private : int axial_pos; int tangential_pos; float bin_value; + int first_energy_window; + int second_energy_window; int time_frame; diff --git a/src/include/stir/Bin.inl b/src/include/stir/Bin.inl index 663f9853f6..7cf716b187 100644 --- a/src/include/stir/Bin.inl +++ b/src/include/stir/Bin.inl @@ -30,9 +30,9 @@ Bin::Bin() Bin::Bin(int segment_num,int view_num, int axial_pos_num,int tangential_pos_num,float bin_value) :segment(segment_num),view(view_num), - axial_pos(axial_pos_num),tangential_pos(tangential_pos_num),bin_value(bin_value),time_frame(1) + axial_pos(axial_pos_num),tangential_pos(tangential_pos_num),bin_value(bin_value), + first_energy_window(0), second_energy_window(0), time_frame(1) {} - int Bin:: axial_pos_num()const @@ -50,6 +50,14 @@ int Bin::view_num() const { return view;} +int + Bin::first_energy_window_num() const +{ return first_energy_window;} + +int + Bin::second_energy_window_num() const +{ return second_energy_window;} + int Bin:: time_frame_num() const {return time_frame;} @@ -70,6 +78,14 @@ int& Bin:: view_num() { return view;} +int& + Bin::first_energy_window_num() +{ return first_energy_window;} + +int& + Bin::second_energy_window_num() +{ return second_energy_window;} + int& Bin:: time_frame_num() {return time_frame;} @@ -99,6 +115,24 @@ void Bin::set_bin_value( float v ) { bin_value = v ;} + +float +Bin::get_first_energy_window_num()const +{ return first_energy_window;} + +void +Bin::set_first_energy_window_num( int v ) +{ first_energy_window = v ;} + +float +Bin::get_second_energy_window_num()const +{ return second_energy_window;} + +void +Bin::set_second_energy_window_num( int v ) +{ second_energy_window = v ;} + + Bin& Bin::operator+=(const float dx) { bin_value+=dx; diff --git a/src/include/stir/ExamInfo.h b/src/include/stir/ExamInfo.h index 69e1a7f594..d0258fbbe5 100644 --- a/src/include/stir/ExamInfo.h +++ b/src/include/stir/ExamInfo.h @@ -26,8 +26,6 @@ #include "stir/Radionuclide.h" #include "stir/shared_ptr.h" -#include "stir/shared_ptr.h" - START_NAMESPACE_STIR @@ -53,13 +51,20 @@ public : ExamInfo() : start_time_in_secs_since_1970(0.), + calibration_factor(-1.F) + { - calibration_factor(-1.F), - low_energy_thres(-1.F), - up_energy_thres(-1.F) + num_windows = 1; + low_energy_thres.resize(1); + up_energy_thres.resize(1); + en_win_pair.resize(2); + low_energy_thres[0]=-1.F; + up_energy_thres[0]=-1.F; + en_win_pair[0]=1; + en_win_pair[1]=1; + + } - { - } std::string originating_system; @@ -83,21 +88,34 @@ public : //! \name Functions that return info related on the acquisition settings //@{ //! Get the low energy boundary - inline float get_low_energy_thres() const; + inline float get_low_energy_thres(int en_window = 0) const; //! Get the high energy boundary - inline float get_high_energy_thres() const; + inline float get_high_energy_thres(int en_window = 0) const; + + //! Get the number of energy windows + inline int get_num_energy_windows() const; + inline std::pair get_energy_window_pair() const; //! Get the calibration factor inline float get_calibration_factor() const; //! Get the radionuclide name inline Radionuclide get_radionuclide() const; //@} + //! \name Functions that set values related on the acquisition settings //@{ //! Set the low energy boundary - inline void set_low_energy_thres(float new_val); + inline void set_low_energy_thres(float new_val,int en_window = 0); //! Set the high energy boundary - inline void set_high_energy_thres(float new_val); + inline void set_high_energy_thres(float new_val,int en_window = 0); + //! Set the low energy boundary as a vector (currently only used in listmode code) + inline void set_low_energy_thres_vect(std::vector new_val); + //! Set the high energy boundary as a vector (currently only used in listmode code) + inline void set_high_energy_thres_vect(std::vector new_val); + //! Set the number of energy windows + inline void set_num_energy_windows(int n_win); + //! Set the energy window pair + inline void set_energy_window_pair(std::vector val); //! Set the Calibration factor inline void set_calibration_factor(const float cal_val); @@ -109,7 +127,8 @@ public : inline bool has_energy_information() const { - return (low_energy_thres > 0.f)&&(up_energy_thres > 0.f); + //TODO ADD LOOP: for (int i=0; i< num_windows; i++) + return (low_energy_thres[0] > 0.f)&&(up_energy_thres[0] > 0.f); } //! Standard trick for a 'virtual copy-constructor' @@ -139,15 +158,23 @@ public : float calibration_factor; private: - //! + //! + //! \brief num_windows + //! \author Ludovica Brusaferri + //! \details This is the value of the number of enegry windows + //! This parameter was initially introduced for scatter simulation. + //! If scatter simulation is not needed, can default to 1 + int num_windows; + + //! //! \brief low_energy_thres //! \author Nikos Efthimiou //! \details This is the value of low energy threshold of the energy window. //! The units are keV //! This parameter was initially introduced for scatter simulation. //! If scatter simulation is not needed, can default to -1 - float low_energy_thres; - + std::vector low_energy_thres; + //! //! \brief up_energy_thres //! \author Nikos Efthimiou @@ -155,7 +182,15 @@ public : //! The units are keV //! This parameter was initially introduced for scatter simulation //! If scatter simulation is not needed, can default to -1 - float up_energy_thres; + std::vector up_energy_thres; + + //! + //! \brief en_win_pair + //! \author Ludovica Brusaferri + //! \details This is the value of the energy window pair for a certain sinogram + //! Can default to {1,1} + std::vector en_win_pair; + }; END_NAMESPACE_STIR diff --git a/src/include/stir/ExamInfo.inl b/src/include/stir/ExamInfo.inl index df55c7320b..7a17504208 100644 --- a/src/include/stir/ExamInfo.inl +++ b/src/include/stir/ExamInfo.inl @@ -31,17 +31,46 @@ create_shared_clone() const } void -ExamInfo::set_low_energy_thres(float new_val) +ExamInfo::set_low_energy_thres(float new_val,int en_window) +{ + low_energy_thres[en_window] = new_val; + +} + +void +ExamInfo::set_high_energy_thres(float new_val,int en_window) +{ + up_energy_thres[en_window] = new_val; + +} + +void +ExamInfo::set_low_energy_thres_vect(std::vector new_val) { low_energy_thres = new_val; } void -ExamInfo::set_high_energy_thres(float new_val) +ExamInfo::set_high_energy_thres_vect(std::vector new_val) { up_energy_thres = new_val; } +void +ExamInfo::set_num_energy_windows(int n_win) +{ + num_windows = n_win; + +} + + +void +ExamInfo::set_energy_window_pair(std::vector val) +{ + en_win_pair=val; +} + + void ExamInfo::set_calibration_factor( const float cal_val) { @@ -54,18 +83,54 @@ ExamInfo::set_radionuclide(const Radionuclide& arg) radionuclide = arg; } +//Get the lower energy boundary for all the energy windows. en_window is set to 0 by default +//So that it will work also in the case of 1 energy window float -ExamInfo::get_low_energy_thres() const +ExamInfo::get_low_energy_thres(int en_window) const { - return low_energy_thres; + return low_energy_thres[en_window]; + } +//Get the high energy boundary for all the energy windows. en_window is set to 0 by default +//So that it will work also in the case of 1 energy window + float -ExamInfo::get_high_energy_thres() const +ExamInfo::get_high_energy_thres(int en_window) const { - return up_energy_thres; + + return up_energy_thres[en_window]; + +} + + +//Get the number of energy windows +int +ExamInfo::get_num_energy_windows() const +{ + + return num_windows; + + } +std::pair +ExamInfo::get_energy_window_pair() const +{ + +std::pair pair; +pair.first = 0; +pair.second = 0; + + +pair.first=en_win_pair[0]; +pair.second=en_win_pair[1]; + + return pair; + +} + + float ExamInfo::get_calibration_factor() const { @@ -85,4 +150,5 @@ ExamInfo::set_energy_information_from(const ExamInfo& other) this->low_energy_thres = other.low_energy_thres; } + END_NAMESPACE_STIR diff --git a/src/include/stir/IO/InputStreamFromROOTFile.h b/src/include/stir/IO/InputStreamFromROOTFile.h index 2b62cee969..fa2e5c9a47 100644 --- a/src/include/stir/IO/InputStreamFromROOTFile.h +++ b/src/include/stir/IO/InputStreamFromROOTFile.h @@ -76,8 +76,9 @@ class InputStreamFromROOTFile : public RegisteredObject< InputStreamFromROOTFile //! constructor InputStreamFromROOTFile(std::string filename, std::string chain_name, - bool exclude_scattered, bool exclude_randoms, - float low_energy_window, float up_energy_window, + int maximum_order_of_scatter, bool exclude_randoms, + float low_energy_window_1, float up_energy_window_1, + float low_energy_window_2, float up_energy_window_2, int offset_dets); #endif @@ -138,9 +139,11 @@ class InputStreamFromROOTFile : public RegisteredObject< InputStreamFromROOTFile //@} //! Lower energy threshold - inline float get_low_energy_thres() const; + inline std::vector get_low_energy_thres_in_keV() const; //! Upper energy threshold - inline float get_up_energy_thres() const; + inline std::vector get_up_energy_thres_in_keV() const; + //! Number of energy windows + inline int get_number_of_energy_windows() const; //! Set singles_readout_depth inline void set_singles_readout_depth(int); @@ -149,6 +152,7 @@ class InputStreamFromROOTFile : public RegisteredObject< InputStreamFromROOTFile inline void set_chain_name(const std::string&); + inline void set_maximum_order_of_scatter(int); inline void set_exclude_true_events(bool); inline void set_exclude_scattered_events(bool); @@ -159,9 +163,9 @@ class InputStreamFromROOTFile : public RegisteredObject< InputStreamFromROOTFile inline void set_detectors_offset(int); - inline void set_low_energy_window(float); + inline void set_low_energy_window(const std::vector &); - inline void set_upper_energy_window(float); + inline void set_upper_energy_window(const std::vector &); //! Set the read_optional_root_fields flag inline void set_optional_ROOT_fields(bool); @@ -254,20 +258,29 @@ class InputStreamFromROOTFile : public RegisteredObject< InputStreamFromROOTFile int num_virtual_axial_crystals_per_block; int num_virtual_transaxial_crystals_per_block; //@} + //! Skip True events (eventID1 == eventID2). Default is false bool exclude_nonrandom; //! Skip scattered events (comptonphantom1 > 0 && comptonphantom2 > 0). Default is false bool exclude_scattered; //! Skip unscattered events (comptonphantom1 == 0 && comptonphantom2 == 0)). Default is false bool exclude_unscattered; + //! Skip scattered events (comptonphantom1 + comptonphantom2 > maximum_order_of_scatter); + int maximum_order_of_scatter; //! Skip random events (eventID1 != eventID2). Default is false bool exclude_randoms; //! Check energy window information (low_energy_window < energy < up_energy_window). Default is true bool check_energy_window_information; - //! Lower energy threshold. Default is 1000 (keV) - float low_energy_window; - //! Upper energy threshold. Default is 0 (keV) - float up_energy_window; + //! Lower energy threshold photon 1. Default is 1000 (keV) + float low_energy_window_1; + //! Upper energy threshold photon 1. Default is 0 (keV) + float up_energy_window_1; + //! Lower energy threshold photon 2. Default is 1000 (keV) + float low_energy_window_2; + //! Upper energy threshold photon 2. Default is 0 (keV) + float up_energy_window_2; + //! Number of energy windows + int num_en_windows; //! This value will apply a rotation on the detectors' id in the same ring. int offset_dets; //!For the singles_readout_depth from GATE's online documentation: diff --git a/src/include/stir/IO/InputStreamFromROOTFile.inl b/src/include/stir/IO/InputStreamFromROOTFile.inl index 99fe7840c1..b1cb5bb78a 100644 --- a/src/include/stir/IO/InputStreamFromROOTFile.inl +++ b/src/include/stir/IO/InputStreamFromROOTFile.inl @@ -65,6 +65,13 @@ get_saved_get_positions() const return saved_get_positions; } +int +InputStreamFromROOTFile:: +get_number_of_energy_windows() const +{ + return num_en_windows; +} + void InputStreamFromROOTFile:: set_saved_get_positions(const std::vector &poss) @@ -72,17 +79,25 @@ set_saved_get_positions(const std::vector &poss) saved_get_positions = poss; } -float +std::vector InputStreamFromROOTFile:: -get_low_energy_thres() const +get_low_energy_thres_in_keV() const { + std::vector low_energy_window; + low_energy_window.resize(2); + low_energy_window[0]=1e3*low_energy_window_1; + low_energy_window[1]=1e3*low_energy_window_2; return low_energy_window; } -float +std::vector InputStreamFromROOTFile:: -get_up_energy_thres() const +get_up_energy_thres_in_keV() const { + std::vector up_energy_window; + up_energy_window.resize(2); + up_energy_window[0]=1e3*up_energy_window_1; + up_energy_window[1]=1e3*up_energy_window_2; return up_energy_window; } @@ -111,6 +126,12 @@ InputStreamFromROOTFile::set_chain_name(const std::string& val) chain_name = val; } +void +InputStreamFromROOTFile::set_maximum_order_of_scatter(int val) +{ + maximum_order_of_scatter = val; +} + void InputStreamFromROOTFile::set_exclude_true_events(bool val) { @@ -120,7 +141,7 @@ InputStreamFromROOTFile::set_exclude_true_events(bool val) void InputStreamFromROOTFile::set_exclude_scattered_events(bool val) { - exclude_scattered = val; + exclude_scatter = val; } void @@ -142,15 +163,17 @@ InputStreamFromROOTFile::set_detectors_offset(int val) } void -InputStreamFromROOTFile::set_low_energy_window(float val) +InputStreamFromROOTFile::set_low_energy_window(const std::vector &val) { - low_energy_window = val; + low_energy_window_1 = val[0]; + low_energy_window_2 = val[1]; } void -InputStreamFromROOTFile::set_upper_energy_window(float val) +InputStreamFromROOTFile::set_upper_energy_window(const std::vector &val) { - up_energy_window = val; + up_energy_window_1 = val[0]; + up_energy_window_1 = val[1]; } void diff --git a/src/include/stir/IO/InputStreamFromROOTFileForCylindricalPET.h b/src/include/stir/IO/InputStreamFromROOTFileForCylindricalPET.h index 6f239b75b2..2d5e8f5dbb 100644 --- a/src/include/stir/IO/InputStreamFromROOTFileForCylindricalPET.h +++ b/src/include/stir/IO/InputStreamFromROOTFileForCylindricalPET.h @@ -94,8 +94,9 @@ class InputStreamFromROOTFileForCylindricalPET : public int submodule_repeater_x, int submodule_repeater_y, int submodule_repeater_z, int module_repeater_x, int module_repeater_y, int module_repeater_z, int rsector_repeater, - bool exclude_scattered, bool exclude_randoms, - float low_energy_window, float up_energy_window, + int maximum_order_of_scatter, bool exclude_randoms, + float low_energy_window_1, float up_energy_window_1, + float low_energy_window_2, float up_energy_window_2, int offset_dets); #endif diff --git a/src/include/stir/IO/InputStreamFromROOTFileForECATPET.h b/src/include/stir/IO/InputStreamFromROOTFileForECATPET.h index b4a3342c8f..6132beacd9 100644 --- a/src/include/stir/IO/InputStreamFromROOTFileForECATPET.h +++ b/src/include/stir/IO/InputStreamFromROOTFileForECATPET.h @@ -88,8 +88,9 @@ class InputStreamFromROOTFileForECATPET : public std::string chain_name, int crystal_repeater_x, int crystal_repeater_y, int crystal_repeater_z, int blocks_repeater_y, int blocks_repeater_z, - bool exclude_scattered, bool exclude_randoms, - float low_energy_window, float up_energy_window, + int maximum_order_of_scatter, bool exclude_randoms, + float low_energy_window_1, float up_energy_window_1, + float low_energy_window_2, float up_energy_window_2, int offset_dets); #endif diff --git a/src/include/stir/IO/InterfileHeader.h b/src/include/stir/IO/InterfileHeader.h index f15e3c1a2d..0847690468 100644 --- a/src/include/stir/IO/InterfileHeader.h +++ b/src/include/stir/IO/InterfileHeader.h @@ -135,6 +135,7 @@ class InterfileHeader : public MinimalInterfileHeader virtual void set_version_specific_keys(); virtual void read_matrix_info(); virtual void read_num_energy_windows(); + virtual void en_window_pair_set(); void read_frames_info(); //! \brief Get the number of datasets /*! To be overloaded by derived classes if multiple "dimensions" are supported. @@ -182,6 +183,8 @@ public : //! \brief upper_en_window_thresholds //! \details High energy window limit std::vector upper_en_window_thresholds; + + std::vector energy_window_pair; // end acquisition parameters protected: @@ -243,7 +246,7 @@ class InterfilePDFSHeader : public InterfileHeader std::vector segment_sequence; std::vector min_ring_difference; - std::vector max_ring_difference; + std::vector max_ring_difference; std::vector num_rings_per_segment; std::vector applied_corrections; diff --git a/src/include/stir/IO/InterfileHeaderSiemens.h b/src/include/stir/IO/InterfileHeaderSiemens.h index c492da24a3..d17a96cbda 100644 --- a/src/include/stir/IO/InterfileHeaderSiemens.h +++ b/src/include/stir/IO/InterfileHeaderSiemens.h @@ -92,6 +92,7 @@ class InterfileImageHeader : public InterfileHeaderSiemens protected: virtual void read_matrix_info(); + virtual void read_energy_window_info(); //! Returns false if OK, true if not. virtual bool post_processing(); diff --git a/src/include/stir/KeyParser.h b/src/include/stir/KeyParser.h index 96db2e2dcf..9706ad8dfa 100644 --- a/src/include/stir/KeyParser.h +++ b/src/include/stir/KeyParser.h @@ -77,6 +77,10 @@ class map_element { public : KeyArgument::type type; +#if 0 + trilian is_vectored; +#endif + void (KeyParser::*p_object_member)(); // pointer to a member function //TODO void (*p_object_member)(); void *p_object_variable; // pointer to a variable diff --git a/src/include/stir/ProjDataInfo.h b/src/include/stir/ProjDataInfo.h index 09faf76d56..d2f91d6853 100644 --- a/src/include/stir/ProjDataInfo.h +++ b/src/include/stir/ProjDataInfo.h @@ -320,7 +320,7 @@ class ProjDataInfo */ virtual Bin - get_bin(const LOR&) const = 0; + get_bin(const LOR&, const std::pair &energy_window_pair = std::pair(1,1)) const = 0; //! \name Equality of ProjDataInfo objects //@{ diff --git a/src/include/stir/ProjDataInfoCylindrical.inl b/src/include/stir/ProjDataInfoCylindrical.inl index 5297e87316..12486cfa54 100644 --- a/src/include/stir/ProjDataInfoCylindrical.inl +++ b/src/include/stir/ProjDataInfoCylindrical.inl @@ -246,7 +246,11 @@ get_segment_axial_pos_num_for_ring_pair(int& segment_num, const int ring2) const { assert(0<=ring1); + //std::cout << "RING1" << ring1 << std::endl; + //std::cout << "NUM RINGS" << get_scanner_ptr()->get_num_rings()<< std::endl; + assert(ring1get_num_rings()); + assert(0<=ring2); assert(ring2get_num_rings()); diff --git a/src/include/stir/ProjDataInfoCylindricalArcCorr.h b/src/include/stir/ProjDataInfoCylindricalArcCorr.h index 6dd02b4176..d69673b222 100644 --- a/src/include/stir/ProjDataInfoCylindricalArcCorr.h +++ b/src/include/stir/ProjDataInfoCylindricalArcCorr.h @@ -70,7 +70,7 @@ class ProjDataInfoCylindricalArcCorr : public ProjDataInfoCylindrical virtual Bin - get_bin(const LOR&) const; + get_bin(const LOR&, const std::pair&) const; virtual std::string parameter_info() const; private: diff --git a/src/include/stir/ProjDataInfoCylindricalNoArcCorr.h b/src/include/stir/ProjDataInfoCylindricalNoArcCorr.h index 477d8536e3..2f13a3f57e 100644 --- a/src/include/stir/ProjDataInfoCylindricalNoArcCorr.h +++ b/src/include/stir/ProjDataInfoCylindricalNoArcCorr.h @@ -229,7 +229,7 @@ class ProjDataInfoCylindricalNoArcCorr : public ProjDataInfoCylindrical inline Succeeded get_bin_for_det_pair(Bin&, const int det1_num, const int ring1_num, - const int det2_num, const int ring2_num) const; + const int det2_num, const int ring2_num, const std::pair &energy_window_pair = std::pair(1,1)) const; //! This routine gets the detector pair corresponding to a bin. @@ -249,7 +249,7 @@ class ProjDataInfoCylindricalNoArcCorr : public ProjDataInfoCylindrical virtual Bin - get_bin(const LOR&) const; + get_bin(const LOR&, const std::pair&) const; //! \name set of obsolete functions to go between bins<->LORs (will disappear!) diff --git a/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl b/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl index 596b1269d6..eac49ec32a 100644 --- a/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl +++ b/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl @@ -127,12 +127,20 @@ Succeeded ProjDataInfoCylindricalNoArcCorr:: get_bin_for_det_pair(Bin& bin, const int det_num1, const int ring_num1, - const int det_num2, const int ring_num2) const + const int det_num2, const int ring_num2, const std::pair& energy_window_pair) const { if (get_view_tangential_pos_num_for_det_num_pair(bin.view_num(), bin.tangential_pos_num(), det_num1, det_num2)) + { + bin.first_energy_window_num() = energy_window_pair.first;//energy_window_pair.first; + bin.second_energy_window_num()=energy_window_pair.second; return get_segment_axial_pos_num_for_ring_pair(bin.segment_num(), bin.axial_pos_num(), ring_num1, ring_num2); + } else + { + bin.first_energy_window_num()=energy_window_pair.second; + bin.second_energy_window_num()=energy_window_pair.first; return get_segment_axial_pos_num_for_ring_pair(bin.segment_num(), bin.axial_pos_num(), ring_num2, ring_num1); + } } Succeeded diff --git a/src/include/stir/ProjDataInfoGenericNoArcCorr.h b/src/include/stir/ProjDataInfoGenericNoArcCorr.h index 621413f01c..399a2ae1d9 100644 --- a/src/include/stir/ProjDataInfoGenericNoArcCorr.h +++ b/src/include/stir/ProjDataInfoGenericNoArcCorr.h @@ -226,7 +226,7 @@ class ProjDataInfoGenericNoArcCorr : public ProjDataInfoGeneric virtual Bin - get_bin(const LOR&) const; + get_bin(const LOR&, const std::pair&) const; //! \name set of obsolete functions to go between bins<->LORs (will disappear!) diff --git a/src/include/stir/ProjDataInfoSubsetByView.h b/src/include/stir/ProjDataInfoSubsetByView.h index a497d07c1e..2a5a789aab 100644 --- a/src/include/stir/ProjDataInfoSubsetByView.h +++ b/src/include/stir/ProjDataInfoSubsetByView.h @@ -159,7 +159,7 @@ class ProjDataInfoSubsetByView : public ProjDataInfo //! Find the bin in the projection data that 'contains' an LOR /*! Forwards ProjDataInfo::get_bin */ - Bin get_bin(const LOR&) const override; + Bin get_bin(const LOR&, const std::pair&) const override; //! Check if \c *this contains \c proj /*! diff --git a/src/include/stir/extend_projdata.h b/src/include/stir/extend_projdata.h index 53ef731886..c827f9f300 100644 --- a/src/include/stir/extend_projdata.h +++ b/src/include/stir/extend_projdata.h @@ -35,6 +35,15 @@ extend_segment_in_views(const SegmentBySinogram& sino, Array<2,float> extend_sinogram_in_views(const Sinogram& sino, const int min_view_extension, const int max_view_extension); + +Array<3,float> +transpose_extend_segment_in_views(const SegmentBySinogram& sino, + const float min_view_compression, const float max_view_compression); +Array<2,float> +transpose_extend_sinogram_in_views(const Sinogram& sino, + const float min_view_compression, const float max_view_compression); + + //@} END_NAMESPACE_STIR diff --git a/src/include/stir/interpolate_projdata.h b/src/include/stir/interpolate_projdata.h index 114f533a2c..a10c004bcb 100644 --- a/src/include/stir/interpolate_projdata.h +++ b/src/include/stir/interpolate_projdata.h @@ -49,18 +49,31 @@ template class SegmentBySinogram; the input sinogram. */ //@{ -Succeeded +Succeeded interpolate_projdata(ProjData& proj_data_out, const ProjData& proj_data_in, const BSpline::BSplineType spline_type, const bool remove_interleaving = false, - const bool use_view_offset = false); -Succeeded + const bool use_view_offset = false); +Succeeded interpolate_projdata(ProjData& proj_data_out, const ProjData& proj_data_in, const BasicCoordinate<3, BSpline::BSplineType> & spline_type, const bool remove_interleaving = false, - const bool use_view_offset = false); + const bool use_view_offset = false); + +Succeeded +interpolate_projdata_pull(ProjData& proj_data_out, + const ProjData& proj_data_in, + const bool remove_interleaving = false, + const bool use_view_offset = false); + +Succeeded +interpolate_projdata_push(ProjData& proj_data_out, + const ProjData& proj_data_in, + const bool remove_interleaving = false, + const bool use_view_offset = false); + //@} END_NAMESPACE_STIR diff --git a/src/include/stir/inverse_SSRB.h b/src/include/stir/inverse_SSRB.h index b6b832cefb..f55c93913e 100644 --- a/src/include/stir/inverse_SSRB.h +++ b/src/include/stir/inverse_SSRB.h @@ -48,5 +48,8 @@ Succeeded inverse_SSRB(ProjData& proj_data_4D, const ProjData& proj_data_3D); +Succeeded +transpose_inverse_SSRB(ProjData& proj_data_4D, + const ProjData& proj_data_3D); END_NAMESPACE_STIR diff --git a/src/include/stir/listmode/CListModeDataROOT.h b/src/include/stir/listmode/CListModeDataROOT.h index 64b73d7409..a3dc4a498b 100644 --- a/src/include/stir/listmode/CListModeDataROOT.h +++ b/src/include/stir/listmode/CListModeDataROOT.h @@ -21,6 +21,7 @@ #define __stir_listmode_CListModeDataROOT_H__ #include "stir/listmode/CListModeData.h" +#include "stir/listmode/ListEnergy.h" #include "stir/listmode/CListRecordROOT.h" #include "stir/IO/InputStreamFromROOTFile.h" #include "stir/KeyParser.h" diff --git a/src/include/stir/listmode/CListRecordECAT8_32bit.h b/src/include/stir/listmode/CListRecordECAT8_32bit.h index ec834510c6..cd962745e3 100644 --- a/src/include/stir/listmode/CListRecordECAT8_32bit.h +++ b/src/include/stir/listmode/CListRecordECAT8_32bit.h @@ -28,6 +28,7 @@ #include "boost/static_assert.hpp" #include "boost/cstdint.hpp" #include "stir/DetectionPositionPair.h" +#include "stir/listmode/ListEnergy.h" START_NAMESPACE_STIR namespace ecat { @@ -186,6 +187,20 @@ class CListTimeECAT8_32bit : public ListTime }; }; +//! A class for storing and using an energy 'event' from a listmode file from the ECAT 8_32bit scanner +/*! \ingroup listmode + */ +class CListEnergyECAT8_32bit : public ListEnergy +{ + public: + bool is_energy() const + { return true; } + inline double get_energyA_in_keV() const + { return 0.F; } + inline double get_energyB_in_keV() const + { return 0.F; } +}; + //! A class for a general element of a listmode file for a Siemens scanner using the ECAT8 32bit format. /*! \ingroup listmode We currently only support coincidence events and a timing flag. @@ -202,6 +217,8 @@ class CListTimeECAT8_32bit : public ListTime bool is_time() const { return this->any_data.is_time(); } + bool is_energy() const + { return true; } /* bool is_gating_input() const { return this->is_time(); } @@ -216,7 +233,10 @@ class CListTimeECAT8_32bit : public ListTime { return this->time_data; } virtual const CListTimeECAT8_32bit& time() const { return this->time_data; } - + virtual CListEnergyECAT8_32bit& energy() + { return this->energy_data; } + virtual const CListEnergyECAT8_32bit& energy() const + { return this->energy_data; } bool operator==(const CListRecord& e2) const { return dynamic_cast(&e2) != 0 && @@ -256,6 +276,7 @@ class CListTimeECAT8_32bit : public ListTime private: CListEventECAT8_32bit event_data; CListTimeECAT8_32bit time_data; + CListEnergyECAT8_32bit energy_data; CListDataAnyECAT8_32bit any_data; boost::int32_t raw; // this raw field isn't strictly necessary, get rid of it? diff --git a/src/include/stir/listmode/CListRecordECAT962.h b/src/include/stir/listmode/CListRecordECAT962.h index 2c06b48493..67b5cc7642 100644 --- a/src/include/stir/listmode/CListRecordECAT962.h +++ b/src/include/stir/listmode/CListRecordECAT962.h @@ -23,6 +23,7 @@ #include "stir/listmode/CListRecord.h" #include "stir/listmode/ListGatingInput.h" #include "stir/listmode/ListTime.h" +#include "stir/listmode/ListEnergy.h" #include "stir/listmode/CListEventCylindricalScannerWithViewTangRingRingEncoding.h" #include "stir/ProjDataInfoCylindrical.h" #include "stir/ProjDataInfoCylindricalNoArcCorr.h" @@ -175,10 +176,18 @@ class CListTimeDataECAT962 #endif }; +//! A class for storing and using an energy 'event' from a listmode file +/*! \ingroup listmode + */ +class CListEnergyDataECAT962 +{ + public: +}; + //! A class for a general element of a listmode file /*! \ingroup listmode For the 962 it's either a coincidence event, or a timing flag.*/ -class CListRecordECAT962 : public CListRecordWithGatingInput, public ListTime, public ListGatingInput, +class CListRecordECAT962 : public CListRecordWithGatingInput, public ListTime, public ListEnergy, public ListGatingInput, public CListEventCylindricalScannerWithViewTangRingRingEncoding { public: @@ -192,6 +201,8 @@ class CListRecordECAT962 : public CListRecordWithGatingInput, public ListTime, p bool is_time() const { return time_data.type == 1U; } + bool is_energy() const + { return true; } bool is_gating_input() const { return this->is_time(); } bool is_event() const @@ -204,6 +215,10 @@ class CListRecordECAT962 : public CListRecordWithGatingInput, public ListTime, p { return *this; } virtual const ListTime& time() const { return *this; } + virtual ListEnergy& energy() + { return *this; } + virtual const ListEnergy& energy() const + { return *this; } virtual ListGatingInput& gating_input() { return *this; } virtual const ListGatingInput& gating_input() const @@ -225,6 +240,12 @@ class CListRecordECAT962 : public CListRecordWithGatingInput, public ListTime, p inline Succeeded set_gating(unsigned int g) { return time_data.set_gating(g); } + // energy + inline double get_energyA_in_keV() const + { return 0.F; } + inline double get_energyB_in_keV() const + { return 0.F; } + // event inline bool is_prompt() const { return event_data.is_prompt(); } inline Succeeded set_prompt(const bool prompt = true) @@ -252,6 +273,7 @@ class CListRecordECAT962 : public CListRecordWithGatingInput, public ListTime, p union { CListEventDataECAT962 event_data; CListTimeDataECAT962 time_data; + CListEnergyDataECAT962 energy_data; boost::int32_t raw; }; BOOST_STATIC_ASSERT(sizeof(boost::int32_t)==4); diff --git a/src/include/stir/listmode/CListRecordECAT966.h b/src/include/stir/listmode/CListRecordECAT966.h index 093c31823b..462af4a58c 100644 --- a/src/include/stir/listmode/CListRecordECAT966.h +++ b/src/include/stir/listmode/CListRecordECAT966.h @@ -23,6 +23,7 @@ #include "stir/listmode/CListRecord.h" #include "stir/listmode/ListGatingInput.h" #include "stir/listmode/ListTime.h" +#include "stir/listmode/ListEnergy.h" #include "stir/listmode/CListEventCylindricalScannerWithViewTangRingRingEncoding.h" #include "stir/ProjDataInfoCylindrical.h" #include "stir/IO/stir_ecat_common.h" // for namespace macros @@ -192,6 +193,20 @@ class CListTimeECAT966 : public ListTime, public ListGatingInput }; }; +//! A class for storing and using a energy'event' from a listmode file from the ECAT 966 scanner +/*! \ingroup listmode + */ +class CListEnergyECAT966 : public ListEnergy +{ + public: + bool is_energy() const + {return true; } + inline double get_energyA_in_keV() const + { return 0.F; } + inline double get_energyB_in_keV() const + { return 0.F; } +}; + //! A class for a general element of a listmode file /*! \ingroup listmode For the 966 it's either a coincidence event, or a timing flag.*/ @@ -202,6 +217,8 @@ class CListRecordECAT966 : public CListRecordWithGatingInput bool is_time() const { return this->time_data.is_time(); } + bool is_energy() const + { return true; } bool is_gating_input() const { return this->is_time(); } bool is_event() const @@ -214,6 +231,10 @@ class CListRecordECAT966 : public CListRecordWithGatingInput { return this->time_data; } virtual const CListTimeECAT966& time() const { return this->time_data; } + virtual CListEnergyECAT966& energy() + { return this->energy_data; } + virtual const CListEnergyECAT966& energy() const + { return this->energy_data; } virtual CListTimeECAT966& gating_input() { return this->time_data; } virtual const CListTimeECAT966& gating_input() const @@ -252,6 +273,7 @@ class CListRecordECAT966 : public CListRecordWithGatingInput private: CListEventECAT966 event_data; CListTimeECAT966 time_data; + CListEnergyECAT966 energy_data; boost::int32_t raw; }; diff --git a/src/include/stir/listmode/CListRecordGEHDF5.h b/src/include/stir/listmode/CListRecordGEHDF5.h index 86a0bec20b..16b198deaa 100644 --- a/src/include/stir/listmode/CListRecordGEHDF5.h +++ b/src/include/stir/listmode/CListRecordGEHDF5.h @@ -19,6 +19,7 @@ #define __stir_listmode_CListRecordGEHDF5_H__ #include "stir/listmode/CListRecord.h" +#include "stir/listmode/ListEnergy.h" #include "stir/listmode/CListEventCylindricalScannerWithDiscreteDetectors.h" #include "stir/Succeeded.h" #include "stir/ByteOrder.h" @@ -181,7 +182,13 @@ namespace RDF_HDF5 { All types of records are stored in a (private) union with the "basic" classes such as CListEventDataGEHDF5. This class essentially just forwards the work to the "basic" classes. */ -class CListRecordGEHDF5 : public CListRecord, public ListTime, // public CListGatingInput, + + class CListEnergyDataGEHDF5 + { + public: + }; + +class CListRecordGEHDF5 : public CListRecord, public ListTime, public ListEnergy, // public CListGatingInput, public CListEventCylindricalScannerWithDiscreteDetectors { typedef detail::CListEventDataGEHDF5 DataType; @@ -213,6 +220,8 @@ class CListRecordGEHDF5 : public CListRecord, public ListTime, // public CListGa bool is_event() const { return this->event_data.is_event(); } + bool is_energy() const + { return true; } virtual CListEvent& event() { return *this; } virtual const CListEvent& event() const @@ -221,6 +230,10 @@ class CListRecordGEHDF5 : public CListRecord, public ListTime, // public CListGa { return *this; } virtual const ListTime& time() const { return *this; } + virtual ListEnergy& energy() + { return *this; } + virtual const ListEnergy& energy() const + { return *this; } #if 0 virtual CListGatingInput& gating_input() { return *this; } @@ -238,6 +251,12 @@ dynamic_cast(&e2) != 0 && #endif } + // energy + inline double get_energyA_in_keV() const + { return 0.F; } + inline double get_energyB_in_keV() const + { return 0.F; } + // time inline unsigned long get_time_in_millisecs() const { return time_data.get_time_in_millisecs() - first_time_stamp; } diff --git a/src/include/stir/listmode/CListRecordROOT.h b/src/include/stir/listmode/CListRecordROOT.h index 1de386da54..29607b6f8f 100644 --- a/src/include/stir/listmode/CListRecordROOT.h +++ b/src/include/stir/listmode/CListRecordROOT.h @@ -119,6 +119,45 @@ class CListTimeROOT : public ListTime double timeB; }; +//! A class for storing and using an energy 'event' from a listmode file +/*! \ingroup listmode + */ +class CListEnergyROOT : public ListEnergy +{ +public: + + void init_energy_from_data(double energy1, double energy2) + { + energyA = energy1; + energyB = energy2; + } + + //! Returns always true + bool is_energy() const + { return true; } + //! Get the detection energy of the first photon + //! in keV + inline double get_energyA_in_keV() const + { return 1e3*energyA; } + //! Get the detection energy of the second photon + //! in keV + inline double get_energyB_in_keV() const + { return 1e3*energyB; } + +private: + + //! + //! \brief energyA + //! \details The detected energy of the first of the two photons, in MeV + double energyA; + + //! + //! \brief energyB + //! \details The detected energy of the second of the two photons, in MeV + double energyB; + + +}; //! A class for a general element of a listmode file for a Siemens scanner using the ROOT files class CListRecordROOT : public CListRecord // currently no gating yet { @@ -128,6 +167,8 @@ class CListRecordROOT : public CListRecord // currently no gating yet //! Returns always true bool inline is_event() const; //! Returns always true + bool inline is_energy() const; + //! Returns always true bool inline is_full_event() const; virtual CListEventROOT& event() @@ -150,6 +191,16 @@ class CListRecordROOT : public CListRecord // currently no gating yet return this->time_data; } + virtual CListEnergyROOT& energy() + { + return this->energy_data; + } + + virtual const CListEnergyROOT& energy() const + { + return this->energy_data; + } + bool operator==(const CListRecord& e2) const { return dynamic_cast(&e2) != 0 && @@ -166,7 +217,8 @@ class CListRecordROOT : public CListRecord // currently no gating yet const int& crystal1, const int& crystal2, double time1, double time2, - const int& event1, const int& event2) + const int& event1, const int& event2, + double energy1 = 0.511F, double energy2 = 0.511F) { /// \warning ROOT data are time and event at the same time. @@ -176,6 +228,7 @@ class CListRecordROOT : public CListRecord // currently no gating yet this->time_data.init_from_data( time1,time2); + this->energy_data.init_energy_from_data(energy1,energy2); // We can make a singature raw based on the two events IDs. // It is pretty unique. raw[0] = event1; @@ -187,6 +240,7 @@ class CListRecordROOT : public CListRecord // currently no gating yet private: CListEventROOT event_data; CListTimeROOT time_data; + CListEnergyROOT energy_data; boost::int32_t raw[2]; // this raw field isn't strictly necessary, get rid of it? }; diff --git a/src/include/stir/listmode/CListRecordROOT.inl b/src/include/stir/listmode/CListRecordROOT.inl index 970c0fbc5c..1bbe58d7e8 100644 --- a/src/include/stir/listmode/CListRecordROOT.inl +++ b/src/include/stir/listmode/CListRecordROOT.inl @@ -25,6 +25,9 @@ bool CListRecordROOT::is_time() const bool CListRecordROOT::is_event() const { return true; } +bool CListRecordROOT::is_energy() const +{ return true; } + bool CListRecordROOT::is_full_event() const { return true; } diff --git a/src/include/stir/listmode/CListRecordSAFIR.h b/src/include/stir/listmode/CListRecordSAFIR.h index cc0e4bc5ec..0482113b15 100644 --- a/src/include/stir/listmode/CListRecordSAFIR.h +++ b/src/include/stir/listmode/CListRecordSAFIR.h @@ -177,8 +177,16 @@ class CListTimeDataSAFIR #endif }; +//! Class for record with time data +class CListEnergyDataSAFIR +{ + +public: + +}; + //! Class for general record, containing a union of data, time and raw record and providing access to certain elements. -class CListRecordSAFIR : public CListRecord, public ListTime, public CListEventSAFIR +class CListRecordSAFIR : public CListRecord, public ListTime, public ListEnergy, public CListEventSAFIR { public: typedef CListEventDataSAFIR DataType; @@ -197,6 +205,9 @@ class CListRecordSAFIR : public CListRecord, public ListTime, public CListEventS virtual bool is_event() const { return time_data.type == 0; } + virtual bool is_energy() const + { return true;} + virtual CListEvent& event() { return *this; } @@ -215,6 +226,12 @@ class CListRecordSAFIR : public CListRecord, public ListTime, public CListEventS virtual const ListTime& time() const { return *this; } + virtual ListEnergy& energy() + { return *this; } + + virtual const ListEnergy& energy() const + { return *this; } + virtual bool operator==(const CListRecord& e2) const { return dynamic_cast(&e2) != 0 && @@ -227,6 +244,11 @@ class CListRecordSAFIR : public CListRecord, public ListTime, public CListEventS inline Succeeded set_time_in_millisecs(const unsigned long time_in_millisecs) { return time_data.set_time_in_millisecs(time_in_millisecs); } + inline double get_energyA_in_keV() const + { return 0.F; } + inline double get_energyB_in_keV() const + { return 0.F; } + inline bool is_prompt() const { return !(event_data.isRandom); } Succeeded @@ -254,6 +276,7 @@ class CListRecordSAFIR : public CListRecord, public ListTime, public CListEventS union { CListEventDataSAFIR event_data; CListTimeDataSAFIR time_data; + CListEnergyDataSAFIR energy_data; boost::int64_t raw; }; BOOST_STATIC_ASSERT(sizeof(boost::uint64_t)==8); diff --git a/src/include/stir/listmode/ListEnergy.h b/src/include/stir/listmode/ListEnergy.h new file mode 100644 index 0000000000..d7128387a4 --- /dev/null +++ b/src/include/stir/listmode/ListEnergy.h @@ -0,0 +1,56 @@ +// +/*! + \file + \ingroup listmode + \brief Declarations of class stir::ListEnergy which + is used for list mode data. + + \author Ludovica Brusaferri + +*/ +/* + Copyright (C) 2019, University College of London + This file is part of STIR. + + This file is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + This file is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + See STIR/LICENSE.txt for details +*/ + +#ifndef __stir_listmode_ListEnergy_H__ +#define __stir_listmode_ListEnergy_H__ + +#include "stir/Succeeded.h" +#include "stir/round.h" + +START_NAMESPACE_STIR +//class Succeeded; + +//! A class for storing and using an energy record from a listmode file +/*! \ingroup listmode + ListEnergy is used to provide an interface to the 'energy' events + in the list mode stream. + + \todo this is still under development. +*/ +class ListEnergy +{ +public: + virtual ~ListEnergy() {} + + virtual double get_energyA_in_keV() const = 0; + virtual double get_energyB_in_keV() const = 0; + +}; + +END_NAMESPACE_STIR + +#endif diff --git a/src/include/stir/listmode/ListEvent.h b/src/include/stir/listmode/ListEvent.h index 6bff69fd2c..4c206756ac 100644 --- a/src/include/stir/listmode/ListEvent.h +++ b/src/include/stir/listmode/ListEvent.h @@ -84,7 +84,7 @@ class ListEvent */ virtual void - get_bin(Bin& bin, const ProjDataInfo&) const; + get_bin(Bin& bin, const ProjDataInfo&, const std::pair &energy_window_pair = std::pair(1,1)) const; //! This method checks if the template is valid for LmToProjData /*! Used before the actual processing of the data (see issue #61), before calling get_bin() diff --git a/src/include/stir/listmode/ListRecord.h b/src/include/stir/listmode/ListRecord.h index 3e82a47b2b..8e680ac236 100644 --- a/src/include/stir/listmode/ListRecord.h +++ b/src/include/stir/listmode/ListRecord.h @@ -26,6 +26,7 @@ #include "ListEvent.h" #include "ListTime.h" +#include "ListEnergy.h" #include "stir/Succeeded.h" START_NAMESPACE_STIR @@ -52,10 +53,14 @@ class ListRecord virtual bool is_event() const = 0; + virtual bool is_energy() const = 0; + virtual ListEvent& event() = 0; virtual const ListEvent& event() const = 0; virtual ListTime& time() = 0; - virtual const ListTime& time() const = 0; + virtual const ListTime& time() const = 0; + virtual ListEnergy& energy() = 0; + virtual const ListEnergy& energy() const = 0; }; END_NAMESPACE_STIR diff --git a/src/include/stir/listmode/LmToProjData.h b/src/include/stir/listmode/LmToProjData.h index 58171f03fa..4bc1f592d2 100644 --- a/src/include/stir/listmode/LmToProjData.h +++ b/src/include/stir/listmode/LmToProjData.h @@ -37,6 +37,7 @@ START_NAMESPACE_STIR class ListEvent; class ListTime; +class ListEnergy; /*! \ingroup listmode @@ -227,7 +228,7 @@ class LmToProjData : public LmToProjDataAbstract (on top of anything done by normalisation_ptr). \todo Would need timing info or so for e.g. time dependent normalisation or angle info for a rotating scanner.*/ - virtual void get_bin_from_event(Bin& bin, const ListEvent&) const; + virtual void get_bin_from_event(Bin& bin, const ListEvent&, const std::pair &energy_window_pair = std::pair(1,1)) const; //! A function that should return the number of uncompressed bins in the current bin /*! \todo it is not compatiable with e.g. HiDAC doesn't belong here anyway @@ -268,6 +269,8 @@ class LmToProjData : public LmToProjDataAbstract bool interactive; shared_ptr template_proj_data_info_ptr; + + //shared_ptr template_exam_info_sptr; //! This will be used for pre-normalisation shared_ptr normalisation_ptr; //! This will be used for post-normalisation diff --git a/src/include/stir/listmode/SPECTListRecord.h b/src/include/stir/listmode/SPECTListRecord.h index d77766c762..f54d4c558c 100644 --- a/src/include/stir/listmode/SPECTListRecord.h +++ b/src/include/stir/listmode/SPECTListRecord.h @@ -25,6 +25,7 @@ #include "stir/listmode/ListTime.h" +#include "stir/listmode/ListEnergy.h" #include "ListRecord.h" #include "stir/listmode/SPECTListEvent.h" @@ -51,11 +52,14 @@ class SPECTListRecord : public ListRecord virtual bool is_event() const = 0; + //virtual bool is_energy() const = 0; + virtual SPECTListEvent& event() = 0; virtual const SPECTListEvent& event() const = 0; virtual ListTime& time() = 0; virtual const ListTime& time() const = 0; - + virtual ListEnergy & energy() = 0; + virtual const ListEnergy& energy() const = 0; virtual bool operator==(const SPECTListRecord& e2) const = 0; bool operator!=(const SPECTListRecord& e2) const { return !(*this == e2); } diff --git a/src/include/stir/numerics/sampling_functions.h b/src/include/stir/numerics/sampling_functions.h index ddc2f75668..046b465ab9 100644 --- a/src/include/stir/numerics/sampling_functions.h +++ b/src/include/stir/numerics/sampling_functions.h @@ -49,6 +49,44 @@ void sample_function_on_regular_grid(Array<3,elemT>& out, const BasicCoordinate<3, positionT>& offset, const BasicCoordinate<3, positionT>& step); +template +inline +void sample_function_on_regular_grid_pull(Array<3,elemT>& out, + const Array<3,elemT>& in, + const BasicCoordinate<3, positionT>& offset, + const BasicCoordinate<3, positionT>& step); + +template +inline +void sample_function_on_regular_grid_push(Array<3,elemT>& out, + const Array<3,elemT>& in, + const BasicCoordinate<3, positionT>& offset, + const BasicCoordinate<3, positionT>& step); +/*template +inline +void axial_position_boundary_conditions(Array<3,elemT>& out); + +template +inline +void views_position_boundary_conditions(Array<3,elemT>& out);*/ + +template +inline +void extend_tangential_position(Array<3,elemT>& out); + +template +inline +void transpose_extend_tangential_position(Array<3,elemT>& out); + + +template +inline +void extend_axial_position(Array<3,elemT>& out); + +template +inline +void transpose_extend_axial_position(Array<3,elemT>& out); + END_NAMESPACE_STIR #include "stir/numerics/sampling_functions.inl" diff --git a/src/include/stir/numerics/sampling_functions.inl b/src/include/stir/numerics/sampling_functions.inl index b504727f97..63606a68bc 100644 --- a/src/include/stir/numerics/sampling_functions.inl +++ b/src/include/stir/numerics/sampling_functions.inl @@ -15,6 +15,7 @@ This file is part of STIR. \brief implementation of stir::sample_function_on_regular_grid */ +#include "stir_experimental/numerics/more_interpolators.h" START_NAMESPACE_STIR template @@ -54,4 +55,179 @@ void sample_function_on_regular_grid(Array<3,elemT>& out, } } +template +void sample_function_on_regular_grid_pull(Array<3,elemT>& out, + const Array<3,elemT>& in, + const BasicCoordinate<3, positionT>& offset, + const BasicCoordinate<3, positionT>& step) +{ + BasicCoordinate<3,int> min_out, max_out; + IndexRange<3> out_range = out.get_index_range(); + if (!out_range.get_regular_range(min_out,max_out)) + warning("Output must be regular range!"); + + BasicCoordinate<3, int> index_out; + BasicCoordinate<3, positionT> relative_positions; + index_out[1]=min_out[1]; + relative_positions[1]= index_out[1] * step[1] - offset[1] ; + const BasicCoordinate<3, positionT> max_relative_positions= + (BasicCoordinate<3,positionT>(max_out)+static_cast(.001)) * step + offset; + for (; + index_out[1]<=max_out[1] && relative_positions[1]<=max_relative_positions[1]; + ++index_out[1], relative_positions[1]+= step[1]) + { + index_out[2]=min_out[2]; + relative_positions[2]= index_out[2] * step[2] + offset[2] ; + for (; + index_out[2]<=max_out[2] && relative_positions[2]<=max_relative_positions[2]; + ++index_out[2], relative_positions[2]+= step[2]) + { + index_out[3]=min_out[3]; + relative_positions[3]= index_out[3] * step[3] + offset[3] ; + for (; + index_out[3]<=max_out[3] && relative_positions[3]<=max_relative_positions[3]; + ++index_out[3], relative_positions[3]+= step[3]) + out[index_out] = pull_linear_interpolate(in, + relative_positions); + } + } +} + + +template +void sample_function_on_regular_grid_push(Array<3,elemT>& out, + const Array<3,elemT>& in, + const BasicCoordinate<3, positionT>& offset, + const BasicCoordinate<3, positionT>& step) +{ + +BasicCoordinate<3,int> min_in, max_in; + + IndexRange<3> in_range = in.get_index_range(); + if (!in_range.get_regular_range(min_in,max_in)) + warning("Output must be regular range!"); + + BasicCoordinate<3, int> index_in; + BasicCoordinate<3, positionT> relative_positions; + index_in[1]=min_in[1]; + relative_positions[1]= index_in[1] *step[1] - offset[1] ; + const BasicCoordinate<3, positionT> max_relative_positions= + (BasicCoordinate<3,positionT>(max_in)+static_cast(.001)) * step + offset; + for (; + index_in[1]<=max_in[1] && relative_positions[1]<=max_relative_positions[1]; + ++index_in[1], relative_positions[1]+= step[1]) + { + index_in[2]=min_in[2]; + relative_positions[2]= index_in[2] * step[2] + offset[2] ; + for (; + index_in[2]<=max_in[2] && relative_positions[2]<=max_relative_positions[2]; + ++index_in[2], relative_positions[2]+= step[2]) + { + index_in[3]=min_in[3]; + relative_positions[3]= index_in[3] * step[3] + offset[3] ; + for (; + index_in[3]<=max_in[3] && relative_positions[3]<=max_relative_positions[3]; + ++index_in[3], relative_positions[3]+= step[3]) + push_transpose_linear_interpolate(out, + relative_positions, + in[index_in]); + } + } + //out*=(step[1]*step[2]*step[3]); //very important +} + +template +void extend_tangential_position(Array<3,elemT>& array) +{ + for (int z=array.get_min_index(); z<= array.get_max_index(); ++z) + { + for (int y=array[z].get_min_index(); y<= array[z].get_max_index(); ++y) + { + const int old_min = array[z][y].get_min_index(); + const int old_max = array[z][y].get_max_index(); + array[z][y].grow(old_min-1, old_max+1); + + array[z][y][old_min-1] = array[z][y][old_min]; + array[z][y][old_max+1] = array[z][y][old_max]; + + } + } +} + + +template +void transpose_extend_tangential_position(Array<3,elemT>& array) +{ + for (int z=array.get_min_index(); z<= array.get_max_index(); ++z) + { + for (int y=array[z].get_min_index(); y<= array[z].get_max_index(); ++y) + { + const int old_min = array[z][y].get_min_index(); + const int old_max = array[z][y].get_max_index(); + array[z][y][old_min+1] += array[z][y][old_min]; + array[z][y][old_max-1] += array[z][y][old_max]; + array[z][y].grow(old_min+1, old_max-1); //resize + } + } +} + + + +template +void extend_axial_position(Array<3,elemT>& array) +{ + + BasicCoordinate<3,int> min, max; + + min[1]=array.get_min_index()-1; + max[1]=array.get_max_index()+1; + min[2]=array[0].get_min_index();; + max[2]=array[0].get_max_index();; + min[3]=array[0][0].get_min_index(); + max[3]=array[0][0].get_max_index(); + + array.grow(IndexRange<3>(min, max)); + + for (int i=array[0].get_min_index(); i<= array[0].get_max_index(); ++i) + { + for (int j=array[0][0].get_min_index(); j<= array[0][0].get_max_index(); ++j) + { + const int min = array.get_min_index(); + const int max = array.get_max_index(); + array[min][i][j]=array[min+1][i][j]; + array[max][i][j]=array[max-1][i][j]; + } + } + +} + + + + +template +void transpose_extend_axial_position(Array<3,elemT>& array) +{ + BasicCoordinate<3,int> min, max; + + min[1]=array.get_min_index()+1; + max[1]=array.get_max_index()-1; + min[2]=array[0].get_min_index();; + max[2]=array[0].get_max_index();; + min[3]=array[0][0].get_min_index(); + max[3]=array[0][0].get_max_index(); + + + for (int i=array[0].get_min_index(); i<= array[0].get_max_index(); ++i) + { + for (int j=array[0][0].get_min_index(); j<= array[0][0].get_max_index(); ++j) + { + + array[min[1]][i][j]+=array[min[1]-1][i][j]; + array[min[1]][i][j]+=array[min[1]+1][i][j]-array[min[1]][i][j]; + array[max[1]][i][j]+=array[max[1]+1][i][j]; + } + } + + array.grow(IndexRange<3>(min, max)); +} END_NAMESPACE_STIR diff --git a/src/include/stir/scatter/ScatterEstimation.h b/src/include/stir/scatter/ScatterEstimation.h index 4fcb8cdf49..b4aa71e2f7 100644 --- a/src/include/stir/scatter/ScatterEstimation.h +++ b/src/include/stir/scatter/ScatterEstimation.h @@ -104,6 +104,41 @@ class ScatterEstimation: public ParsingObject const bool remove_interleaving = true); + static void + upsample_scatter_estimate(ProjData& scaled_scatter_proj_data, + const ProjData& emission_proj_data, + const ProjData& scatter_proj_data, + const bool remove_interleaving = true); + + static void + pull_scatter_estimate(ProjData& scaled_scatter_proj_data, + const ProjData& emission_proj_data, + const ProjData& scatter_proj_data, + const bool remove_interleaving = true); + + static void + push_scatter_estimate(ProjData& scaled_scatter_proj_data, + const ProjData& emission_proj_data, + const ProjData& scatter_proj_data, + const bool remove_interleaving = true); + + static void + pull_scatter_estimate(ProjData& scaled_scatter_proj_data, + const ProjData& emission_proj_data, + const ProjData& scatter_proj_data, + const ProjData& norm, + const bool remove_interleaving=true); + + static void + push_scatter_estimate(ProjData& scaled_scatter_proj_data, + const ProjData& emission_proj_data, + const ProjData& scatter_proj_data, + const ProjData& norm, + const bool remove_interleaving=true); + + static void + apply_norm(ProjData& projdata,const ProjData& norm); + //! Default constructor (calls set_defaults()) ScatterEstimation(); //! Overloaded constructor with parameter file and initialisation diff --git a/src/include/stir/scatter/ScatterSimulation.h b/src/include/stir/scatter/ScatterSimulation.h index 8a2e5209eb..909d57156d 100644 --- a/src/include/stir/scatter/ScatterSimulation.h +++ b/src/include/stir/scatter/ScatterSimulation.h @@ -23,6 +23,7 @@ */ #include "stir/ParsingObject.h" +#include "stir/DataProcessor.h" #include "stir/RegisteredObject.h" #include "stir/ProjData.h" #include "stir/VoxelsOnCartesianGrid.h" @@ -217,6 +218,18 @@ class ScatterSimulation : public RegisteredObject total_Compton_cross_section_relative_to_511keV(const float energy); //@} + float detection_efficiency(const float incoming_photon_energy, const int en_window = 0) const; + float detection_efficiency_numerical_formulation(const float energy, const int en_window = 0) const; + std::vector detection_spectrum_numerical_formulation(const float LLD, const float HLD, const float size, const float incoming_photon_energy) const; + float detection_model_with_fitted_parameters(const float x, const float energy) const; + float photoelectric(const float K, const float std_peak, const float x, const float energy) const; + float compton_plateau(const float K, const float std_peak, const float x, const float energy, const float scaling_std_compton,const float shift_compton) const; + float flat_continuum(const float K, const float std_peak, const float x, const float energy) const; + float exponential_tail(const float K, const float std_peak, const float x, const float energy, const float beta) const; + //! find scatter points + /*! This function sets scatt_points_vector and scatter_volume. It will also + remove any cached integrals as they would be incorrect otherwise. + */ virtual Succeeded set_up(); //! Output the log of the process. @@ -239,7 +252,6 @@ class ScatterSimulation : public RegisteredObject const CartesianCoordinate3D& detector_coord); - virtual void set_defaults(); virtual void initialise_keymap(); //! \warning post_processing will set everything that has a file name in @@ -280,7 +292,6 @@ class ScatterSimulation : public RegisteredObject * @{ */ - float detection_efficiency(const float energy) const; //! maximum angle to consider above which detection after Compton scatter is considered too small @@ -377,6 +388,7 @@ class ScatterSimulation : public RegisteredObject std::string output_proj_data_filename; //! Shared ptr to hold the simulated data. shared_ptr output_proj_data_sptr; + shared_ptr HR_output_proj_data_sptr; //! threshold below which a voxel in the attenuation image will not be considered as a candidate scatter point float attenuation_threshold; diff --git a/src/include/stir/scatter/SingleScatterLikelihoodAndGradient.h b/src/include/stir/scatter/SingleScatterLikelihoodAndGradient.h new file mode 100644 index 0000000000..a1cc2ffcd0 --- /dev/null +++ b/src/include/stir/scatter/SingleScatterLikelihoodAndGradient.h @@ -0,0 +1,106 @@ +/* + Copyright (C) 2016 University College London + This file is part of STIR. + + This file is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + This file is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup scatter + \brief Definition of class stir::SingleScatterLikelihoodAndGradient + + \author Ludovica Brusaferri +*/ + +#ifndef __stir_scatter_SingleScatterLikelihoodAndGradient_H__ +#define __stir_scatter_SingleScatterLikelihoodAndGradient_H__ + +#include "stir/Succeeded.h" +#include "stir/scatter/ScatterSimulation.h" +#include "stir/scatter/SingleScatterSimulation.h" +#include "stir/RegisteredParsingObject.h" +#include "stir/zoom.h" + +START_NAMESPACE_STIR + +class SingleScatterLikelihoodAndGradient : public + RegisteredParsingObject< + SingleScatterLikelihoodAndGradient, + SingleScatterSimulation, + SingleScatterSimulation> +{ +private: + typedef RegisteredParsingObject< + SingleScatterLikelihoodAndGradient, + SingleScatterSimulation, + SingleScatterSimulation> base_type; +public: + + + + //! Name which will be used when parsing a ScatterSimulation object + static const char * const registered_name; + + //! Default constructor + SingleScatterLikelihoodAndGradient(); + + //! Constructor with initialisation from parameter file + + explicit + SingleScatterLikelihoodAndGradient(const std::string& parameter_filename); + + virtual ~SingleScatterLikelihoodAndGradient(); + + + double L_G_function(const ProjData& data,VoxelsOnCartesianGrid& gradient_image, const bool compute_gradient = true ,const bool isgradient_mu = true,const float rescale = 1.F); + double L_G_function(const ProjData& data,const ProjData &add_sino,VoxelsOnCartesianGrid& gradient_image,const bool compute_gradient = true ,const bool isgradient_mu = true,const float rescale = 1.F); + ProjDataInMemory likelihood_and_gradient_scatter(const ProjData &projdata,const ProjData& norm, const ProjData &add_projdata, VoxelsOnCartesianGrid& gradient_image_HR, VoxelsOnCartesianGrid& gradient_image_LR,const bool compute_gradient, const bool isgradient_mu); + ProjDataInMemory get_jacobian(std::vector > &gradient_image_array,const bool compute_gradient, const bool isgradient_mu); + ProjDataInMemory get_ratio(const ProjData& projdata,const ProjData& norm,const ProjData &add_projdata, const ProjData &est_projdata, std::vector &ratio_vector); + protected: + + void + line_contribution(VoxelsOnCartesianGrid& gradient_image,const float scale, + const CartesianCoordinate3D& scatter_point, + const CartesianCoordinate3D& detector_coord, + const float C); + + void + line_contribution_act(VoxelsOnCartesianGrid& gradient_image, + const CartesianCoordinate3D& scatter_point, + const CartesianCoordinate3D& detector_coord, + const float C); + + void + s_contribution(VoxelsOnCartesianGrid& gradient_image, + const CartesianCoordinate3D& scatter_point, + const float D); + float + L_G_for_one_scatter_point(VoxelsOnCartesianGrid& gradient, + const std::size_t scatter_point_num, + const unsigned det_num_A, + const unsigned det_num_B, + const bool compute_gradient, + const bool isgradient_mu); + + double L_G_estimate(VoxelsOnCartesianGrid& gradient_image_bin,const Bin bin, const bool compute_gradient, const bool isgradient_mu); + double L_G_for_view_segment_number(const ProjData&data,const ProjData&add_sino,VoxelsOnCartesianGrid& gradient_image,const ViewSegmentNumbers& vs_num, const float rescale, const bool compute_gradient,const bool isgradient_mu); + double L_G_for_viewgram(const Viewgram& viewgram,const Viewgram& v_add, Viewgram& v_est,VoxelsOnCartesianGrid& gradient_image,const float rescale, const bool compute_gradient,const bool isgradient_mu); + void get_jacobian_for_view_segment_number(std::vector > &gradient_image_array,ProjData &est_data,const ViewSegmentNumbers& vs_num, const bool compute_gradient,const bool isgradient_mu); + void get_jacobian_for_viewgram(Viewgram& v_est, std::vector > &gradient_image_array,const bool compute_gradient, const bool isgradient_mu); + +}; + +END_NAMESPACE_STIR + +#endif diff --git a/src/include/stir/scatter/SingleScatterSimulation.h b/src/include/stir/scatter/SingleScatterSimulation.h index c330910d2c..74bbee7ac1 100644 --- a/src/include/stir/scatter/SingleScatterSimulation.h +++ b/src/include/stir/scatter/SingleScatterSimulation.h @@ -62,7 +62,6 @@ class SingleScatterSimulation : public virtual void ask_parameters(); //! Perform checks and intialisations virtual Succeeded set_up(); -protected: void initialise(const std::string& parameter_filename); @@ -73,6 +72,8 @@ class SingleScatterSimulation : public virtual bool post_processing(); +protected: + //! //! \brief simulate single scatter for one scatter point double diff --git a/src/include/stir/zoom.h b/src/include/stir/zoom.h index cbc78f9dd1..41a4711cc4 100644 --- a/src/include/stir/zoom.h +++ b/src/include/stir/zoom.h @@ -232,6 +232,17 @@ zoom_image_in_place(VoxelsOnCartesianGrid &image, //@} +void +transpose_zoom_image(VoxelsOnCartesianGrid &image_out, + const VoxelsOnCartesianGrid &image_in, + const ZoomOptions = ZoomOptions::preserve_sum); + +void +zoom_image_swig(VoxelsOnCartesianGrid &image_out, + const VoxelsOnCartesianGrid &image_in, const int zoom_option); +void +transpose_zoom_image_swig(VoxelsOnCartesianGrid &image_out, + const VoxelsOnCartesianGrid &image_in, const int zoom_option); END_NAMESPACE_STIR #endif diff --git a/src/include/stir_experimental/listmode/CListRecordLMF.h b/src/include/stir_experimental/listmode/CListRecordLMF.h index ced054ad59..f49ca9f988 100644 --- a/src/include/stir_experimental/listmode/CListRecordLMF.h +++ b/src/include/stir_experimental/listmode/CListRecordLMF.h @@ -71,6 +71,14 @@ class CListTimeDataLMF unsigned long time; // in millisecs TODO }; +//! A class for storing and using an energy 'event' from a listmode file +/*! \ingroup ClearPET_utilities + */ +class CListEnergyDataLMF +{ + public: +}; + //! A class for a general element of a listmode file /*! \ingroup ClearPET_utilities @@ -80,7 +88,7 @@ I store both time and CListEvent info in one CListRecordLMF. That's obviously not necessary nor desirable. */ -class CListRecordLMF : public CListRecord, public ListTime, public CListEvent +class CListRecordLMF : public CListRecord, public ListTime, public ListEnergy, public CListEvent { public: @@ -93,6 +101,8 @@ class CListRecordLMF : public CListRecord, public ListTime, public CListEvent { return is_time_flag == true; } bool is_event() const { return is_time_flag == false; } + bool is_energy() const + { return true; } virtual CListEvent& event() { return *this; } virtual const CListEvent& event() const @@ -101,7 +111,10 @@ class CListRecordLMF : public CListRecord, public ListTime, public CListEvent { return *this; } virtual const ListTime& time() const { return *this; } - + virtual ListEnergy& energy() + { return *this; } + virtual const ListEnergy& energy() const + { return *this; } bool operator==(const CListRecord& e2) const; // time @@ -125,6 +138,7 @@ class CListRecordLMF : public CListRecord, public ListTime, public CListEvent CListEventDataLMF event_data; CListTimeDataLMF time_data; + CListEnergyDataLMF energy_data; bool is_time_flag; }; diff --git a/src/include/stir_experimental/numerics/more_interpolators.h b/src/include/stir_experimental/numerics/more_interpolators.h index 9be43b2001..bdfb7f3fc7 100644 --- a/src/include/stir_experimental/numerics/more_interpolators.h +++ b/src/include/stir_experimental/numerics/more_interpolators.h @@ -17,6 +17,7 @@ */ #include "stir/BasicCoordinate.h" +#include "stir/ProjData.h" #include "stir/Array.h" @@ -80,6 +81,7 @@ PullLinearInterpolator this->_input_ptr = &input; } + template elemT operator()(const BasicCoordinate<3, positionT>& point_in_input_coords) const { diff --git a/src/include/stir_experimental/numerics/more_interpolators.inl b/src/include/stir_experimental/numerics/more_interpolators.inl index c04de3c1d6..d26d20f96d 100644 --- a/src/include/stir_experimental/numerics/more_interpolators.inl +++ b/src/include/stir_experimental/numerics/more_interpolators.inl @@ -13,6 +13,7 @@ */ #include "stir/Coordinate3D.h" +#include "stir/ProjData.h" #include "stir/Array.h" #include "stir/round.h" #include @@ -65,16 +66,17 @@ push_nearest_neighbour_interpolate(Array& out, + template elemT -pull_linear_interpolate(const Array<3, elemT>& in, - const BasicCoordinate<3, positionT>& point_in_input_coords) +pull_linear_interpolate(const Array<3, elemT>& in, + const BasicCoordinate<3, positionT>& point_in_input_coords) { - // find left neighbour - const Coordinate3D + // find left neighbour + const Coordinate3D left_neighbour(round(std::floor(point_in_input_coords[1])), - round(std::floor(point_in_input_coords[2])), - round(std::floor(point_in_input_coords[3]))); + round(std::floor(point_in_input_coords[2])), + round(std::floor(point_in_input_coords[3]))); // TODO handle boundary conditions if (left_neighbour[1] < in.get_max_index() && @@ -97,17 +99,17 @@ pull_linear_interpolate(const Array<3, elemT>& in, const positionT iyc = 1 - iy; const positionT izc = 1 - iz; return - static_cast - ( - ixc * (iyc * (izc * in[z1][y1][x1] - + iz * in[z2][y1][x1]) - + iy * (izc * in[z1][y2][x1] - + iz * in[z2][y2][x1])) - + ix * (iyc * (izc * in[z1][y1][x2] - + iz * in[z2][y1][x2]) - + iy * (izc * in[z1][y2][x2] - + iz * in[z2][y2][x2])) - ); + static_cast + ( + ixc * (iyc * (izc * in[z1][y1][x1] + + iz * in[z2][y1][x1]) + + iy * (izc * in[z1][y2][x1] + + iz * in[z2][y2][x1])) + + ix * (iyc * (izc * in[z1][y1][x2] + + iz * in[z2][y1][x2]) + + iy * (izc * in[z1][y2][x2] + + iz * in[z2][y2][x2])) + ); } else return 0; @@ -117,7 +119,7 @@ template void push_transpose_linear_interpolate(Array<3, elemT>& out, const BasicCoordinate<3, positionT>& point_in_output_coords, - valueT value) + valueT value) { if (value==0) return; @@ -127,6 +129,7 @@ push_transpose_linear_interpolate(Array<3, elemT>& out, round(std::floor(point_in_output_coords[2])), round(std::floor(point_in_output_coords[3]))); + // TODO handle boundary conditions if (left_neighbour[1] < out.get_max_index() && left_neighbour[1] >= out.get_min_index() && diff --git a/src/listmode_buildblock/CListModeDataROOT.cxx b/src/listmode_buildblock/CListModeDataROOT.cxx index b14bd91c69..c278eaebb3 100644 --- a/src/listmode_buildblock/CListModeDataROOT.cxx +++ b/src/listmode_buildblock/CListModeDataROOT.cxx @@ -73,8 +73,8 @@ CListModeDataROOT(const std::string& hroot_filename) // Only PET scanners supported _exam_info_sptr->imaging_modality = ImagingModality::PT; _exam_info_sptr->originating_system = this->originating_system; - _exam_info_sptr->set_low_energy_thres(this->root_file_sptr->get_low_energy_thres()); - _exam_info_sptr->set_high_energy_thres(this->root_file_sptr->get_up_energy_thres()); + _exam_info_sptr->set_high_energy_thres_vect(this->root_file_sptr->get_up_energy_thres_in_keV()); + _exam_info_sptr->set_low_energy_thres_vect(this->root_file_sptr->get_low_energy_thres_in_keV()); this->exam_info_sptr = _exam_info_sptr; diff --git a/src/listmode_buildblock/ListEvent.cxx b/src/listmode_buildblock/ListEvent.cxx index ff9315bdcf..725f1b3c20 100644 --- a/src/listmode_buildblock/ListEvent.cxx +++ b/src/listmode_buildblock/ListEvent.cxx @@ -30,9 +30,9 @@ START_NAMESPACE_STIR void ListEvent:: -get_bin(Bin& bin, const ProjDataInfo& proj_data_info) const +get_bin(Bin& bin, const ProjDataInfo& proj_data_info, const std::pair &energy_window_pair) const { - bin = proj_data_info.get_bin(get_LOR()); + bin = proj_data_info.get_bin(get_LOR(),energy_window_pair); } END_NAMESPACE_STIR diff --git a/src/listmode_buildblock/LmToProjData.cxx b/src/listmode_buildblock/LmToProjData.cxx index 1524ee6563..ce0cbd596f 100644 --- a/src/listmode_buildblock/LmToProjData.cxx +++ b/src/listmode_buildblock/LmToProjData.cxx @@ -460,12 +460,12 @@ LmToProjData(const char * const par_filename) ***************************************************************/ void LmToProjData:: -get_bin_from_event(Bin& bin, const ListEvent& event) const +get_bin_from_event(Bin& bin, const ListEvent& event, const std::pair &energy_window_pair) const { if (do_pre_normalisation) { Bin uncompressed_bin; - event.get_bin(uncompressed_bin, *proj_data_info_cyl_uncompressed_ptr); + event.get_bin(uncompressed_bin, *proj_data_info_cyl_uncompressed_ptr, energy_window_pair); if (uncompressed_bin.get_bin_value()<=0) return; // rejected for some strange reason @@ -499,7 +499,7 @@ get_bin_from_event(Bin& bin, const ListEvent& event) const const float bin_value = 1/bin_efficiency; // TODO wasteful: we decode the event twice. replace by something like // template_proj_data_info_ptr->get_bin_from_uncompressed(bin, uncompressed_bin); - event.get_bin(bin, *template_proj_data_info_ptr); + event.get_bin(bin, *template_proj_data_info_ptr, energy_window_pair); if (bin.get_bin_value()>0) { @@ -509,7 +509,7 @@ get_bin_from_event(Bin& bin, const ListEvent& event) const } else { - event.get_bin(bin, *template_proj_data_info_ptr); + event.get_bin(bin, *template_proj_data_info_ptr, energy_window_pair); } } @@ -632,6 +632,8 @@ process_data() // we have to do this because the first time tag might occur only after a // few coincidence events (as happens with ECAT scanners) current_time = 0; + shared_ptr template_proj_data_ptr = + ProjData::read_from_file(template_proj_data_name); double time_of_last_stored_event = 0; long num_stored_events = 0; @@ -670,7 +672,16 @@ process_data() char rest[50]; sprintf(rest, "_f%dg1d0b0", current_frame_num); const string output_filename = output_filename_prefix + rest; - + this_frame_exam_info.set_num_energy_windows(template_proj_data_ptr->get_exam_info_sptr()->get_num_energy_windows()); + std::vector energy_window_pair(2); + energy_window_pair.at(0) = template_proj_data_ptr->get_exam_info_sptr()->get_energy_window_pair().first; + energy_window_pair.at(1) = template_proj_data_ptr->get_exam_info_sptr()->get_energy_window_pair().second; + this_frame_exam_info.set_energy_window_pair(energy_window_pair); + for (int i = 0; i < template_proj_data_ptr->get_exam_info_sptr()->get_num_energy_windows(); ++i ) + { + this_frame_exam_info.set_high_energy_thres(template_proj_data_ptr->get_exam_info_sptr()->get_low_energy_thres(i),i); + this_frame_exam_info.set_low_energy_thres(template_proj_data_ptr->get_exam_info_sptr()->get_high_energy_thres(i),i); + } proj_data_sptr = construct_proj_data(output, output_filename, this_frame_exam_info, template_proj_data_info_ptr); } @@ -758,17 +769,23 @@ process_data() // set value in case the event decoder doesn't touch it // otherwise it would be 0 and all events will be ignored bin.set_bin_value(1); - bin.time_frame_num() = current_frame_num; - get_bin_from_event(bin, record.event()); + bin.time_frame_num() = current_frame_num; + get_bin_from_event(bin, record.event(), template_proj_data_ptr->get_exam_info_sptr()->get_energy_window_pair()); // check if it's inside the range we want to store if (bin.get_bin_value()>0 && bin.tangential_pos_num()>= proj_data_sptr->get_min_tangential_pos_num() && bin.tangential_pos_num()<= proj_data_sptr->get_max_tangential_pos_num() && bin.axial_pos_num()>=proj_data_sptr->get_min_axial_pos_num(bin.segment_num()) - && bin.axial_pos_num()<=proj_data_sptr->get_max_axial_pos_num(bin.segment_num()) - ) + && bin.axial_pos_num()<=proj_data_sptr->get_max_axial_pos_num(bin.segment_num()) + && record.energy().get_energyA_in_keV() >= template_proj_data_ptr->get_exam_info_sptr()->get_low_energy_thres(bin.first_energy_window_num()-1) + && record.energy().get_energyA_in_keV() <= template_proj_data_ptr->get_exam_info_sptr()->get_high_energy_thres(bin.first_energy_window_num()-1) + && record.energy().get_energyB_in_keV() >= template_proj_data_ptr->get_exam_info_sptr()->get_low_energy_thres(bin.second_energy_window_num()-1) + && record.energy().get_energyB_in_keV() <= template_proj_data_ptr->get_exam_info_sptr()->get_high_energy_thres(bin.second_energy_window_num()-1)) { + + std::cout<< "energy A new: " << record.energy().get_energyA_in_keV()<< '\n'; + std::cout<< "energy B new: " << record.energy().get_energyB_in_keV() << '\n'; assert(bin.view_num()>=proj_data_sptr->get_min_view_num()); assert(bin.view_num()<=proj_data_sptr->get_max_view_num()); diff --git a/src/scatter_buildblock/CMakeLists.txt b/src/scatter_buildblock/CMakeLists.txt index abc5527cf4..e0841b9d81 100644 --- a/src/scatter_buildblock/CMakeLists.txt +++ b/src/scatter_buildblock/CMakeLists.txt @@ -17,6 +17,7 @@ set(${dir_LIB_SOURCES} CreateTailMaskFromACFs.cxx ScatterSimulation.cxx SingleScatterSimulation.cxx + SingleScatterLikelihoodAndGradient.cxx ) #$(dir)_REGISTRY_SOURCES:= scatter_buildblock_registries diff --git a/src/scatter_buildblock/ScatterEstimation.cxx b/src/scatter_buildblock/ScatterEstimation.cxx index 9a1112de1b..9ac010ddb2 100644 --- a/src/scatter_buildblock/ScatterEstimation.cxx +++ b/src/scatter_buildblock/ScatterEstimation.cxx @@ -130,6 +130,10 @@ initialise_keymap() &this->recon_template_par_filename); this->parser.add_parsing_key("reconstruction type", &this->reconstruction_template_sptr); + + this->parser.add_key("reconstruction parameter template file", + &this->recon_template_par_filename); + // END RECONSTRUCTION RELATED this->parser.add_key("number of scatter iterations", @@ -209,6 +213,8 @@ post_processing() %this->recon_template_par_filename); return true; } + //std::cerr << "print"<< local_parser.parameter_info()<< '\n'; + } if (!this->atten_image_filename.empty()) @@ -917,6 +923,7 @@ process_data() local_max_scale_value, this->half_filter_width, spline_type, true); + if(this->run_debug_mode) { std::stringstream convert; // stream used for the conversion @@ -1047,7 +1054,6 @@ process_data() error("ScatterEstimation: You should not be here. This is not 2D."); } current_activity_image_sptr->fill(1.f); - iterative_method ? reconstruct_iterative(i_scat_iter, this->current_activity_image_sptr): reconstruct_analytic(i_scat_iter, this->current_activity_image_sptr); @@ -1056,6 +1062,7 @@ process_data() } + info("ScatterEstimation: Scatter Estimation finished !!!"); return Succeeded::yes; diff --git a/src/scatter_buildblock/ScatterSimulation.cxx b/src/scatter_buildblock/ScatterSimulation.cxx index a3bb4adac2..286e989693 100644 --- a/src/scatter_buildblock/ScatterSimulation.cxx +++ b/src/scatter_buildblock/ScatterSimulation.cxx @@ -36,6 +36,7 @@ #include #include #include +#include #include "stir/zoom.h" #include "stir/SSRB.h" @@ -93,6 +94,32 @@ process_data() error("ScatterSimulation: output projection data incompatible with what was used for set_up()"); // TODO enable check on exam_info but this has no operator== yet } + + //show energy window information + std::cerr << "number of energy windows:= "<< this->template_exam_info_sptr->get_num_energy_windows() << '\n'; + + if(this->template_exam_info_sptr->get_energy_window_pair().first!= -1 && + this->template_exam_info_sptr->get_energy_window_pair().second!= -1 ) + { + std::cerr << "energy window pair :="<<" {"<< this->template_exam_info_sptr->get_energy_window_pair().first << ',' << this->template_exam_info_sptr->get_energy_window_pair().second <<"}\n"; + + } + + + for (int i = 1; i <= this->template_exam_info_sptr->get_num_energy_windows(); ++i) + { + std::cerr << "energy window lower level"<<"["<attenuation_threshold = 0.01f ; - this->randomly_place_scatter_points = true; + this->randomly_place_scatter_points = false; this->use_cache = true; this->zoom_xy = -1.f; this->zoom_z = -1.f; @@ -247,7 +274,7 @@ ScatterSimulation:: ask_parameters() { this->attenuation_threshold = ask_num("attenuation threshold(cm^-1)",0.0f, 5.0f, 0.01f); - this->randomly_place_scatter_points = ask_num("random place scatter points?",0, 1, 1); + this->randomly_place_scatter_points = ask_num("random place scatter points?",0, 1, 0); this->use_cache = ask_num(" Use cache?",0, 1, 1); this->density_image_filename = ask_string("density image filename", ""); this->activity_image_filename = ask_string("activity image filename", ""); @@ -331,7 +358,7 @@ set_up() error("ScatterSimulation: projection data info not set. Aborting."); if(!template_exam_info_sptr->has_energy_information()) - error("ScatterSimulation: template energy window information not set. Aborting."); + error("ScatterSimulation: template energy window information not set. Aborting."); if(is_null_ptr(activity_image_sptr)) error("ScatterSimulation: activity image not set. Aborting."); @@ -697,12 +724,11 @@ void ScatterSimulation:: set_output_proj_data(const std::string& filename) { - + if(is_null_ptr(this->proj_data_info_cyl_noarc_cor_sptr)) { error("ScatterSimulation: Template ProjData has not been set. Aborting."); } - this->output_proj_data_filename = filename; shared_ptr tmp_sptr; @@ -753,7 +779,6 @@ set_template_proj_data_info(const std::string& filename) shared_ptr template_proj_data_sptr(ProjData::read_from_file(this->template_proj_data_filename)); this->set_exam_info(template_proj_data_sptr->get_exam_info()); - this->set_template_proj_data_info(*template_proj_data_sptr->get_proj_data_info_sptr()); } @@ -788,12 +813,20 @@ ScatterSimulation::set_template_proj_data_info(const ProjDataInfo& arg) this->remove_cache_for_integrals_over_activity(); } + + void ScatterSimulation:: set_exam_info(const ExamInfo& arg) { this->_already_set_up = false; this->template_exam_info_sptr = arg.create_shared_clone(); + + for (int i = 0; i < arg.get_num_energy_windows(); ++i) + { + this->template_exam_info_sptr->set_high_energy_thres(arg.get_high_energy_thres(i),i); + this->template_exam_info_sptr->set_low_energy_thres(arg.get_low_energy_thres(i),i); + } } void diff --git a/src/scatter_buildblock/SingleScatterLikelihoodAndGradient.cxx b/src/scatter_buildblock/SingleScatterLikelihoodAndGradient.cxx new file mode 100644 index 0000000000..29cd9a9cca --- /dev/null +++ b/src/scatter_buildblock/SingleScatterLikelihoodAndGradient.cxx @@ -0,0 +1,780 @@ +/* + Copyright (C) 2016 University College London + This file is part of STIR. + + This file is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + This file is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + See STIR/LICENSE.txt for details +*/ +#include "stir/scatter/SingleScatterLikelihoodAndGradient.h" +#include "stir/scatter/ScatterEstimation.h" +#include "stir/scatter/ScatterSimulation.h" +#include "stir/ViewSegmentNumbers.h" +#include "stir/Bin.h" + +#include "stir/Viewgram.h" +#include "stir/is_null_ptr.h" +#include "stir/IO/read_from_file.h" +#include "stir/IO/write_to_file.h" +#include "stir/info.h" +#include "stir/error.h" +#include +#include + +#include "stir/zoom.h" +#include "stir/SSRB.h" + +#include "stir/stir_math.h" +#include "stir/NumericInfo.h" + +START_NAMESPACE_STIR + +const char * const +SingleScatterLikelihoodAndGradient::registered_name = + "Single Scatter Likelihood And Gradient"; + + +SingleScatterLikelihoodAndGradient:: +SingleScatterLikelihoodAndGradient() : + base_type() +{ + this->set_defaults(); +} + + + +SingleScatterLikelihoodAndGradient:: +SingleScatterLikelihoodAndGradient(const std::string& parameter_filename) +{ + this->initialise(parameter_filename); +} + +SingleScatterLikelihoodAndGradient:: +~SingleScatterLikelihoodAndGradient() +{} + +static const float total_Compton_cross_section_511keV = +ScatterSimulation:: +total_Compton_cross_section(511.F); + +double +SingleScatterLikelihoodAndGradient:: +L_G_function(const ProjData& data,VoxelsOnCartesianGrid& gradient_image, const bool compute_gradient, const bool isgradient_mu, const float rescale) +{ + + shared_ptr add_sino(new ProjDataInMemory(this->output_proj_data_sptr->get_exam_info_sptr(), + this->output_proj_data_sptr->get_proj_data_info_sptr()->create_shared_clone())); + add_sino->fill(0.000000000000000000001); //automatically fills an additive sinogram to avoid division by 0 + L_G_function(data,*add_sino,gradient_image,compute_gradient,isgradient_mu,rescale); +} + +double +SingleScatterLikelihoodAndGradient:: +L_G_function(const ProjData& data,const ProjData &add_sino,VoxelsOnCartesianGrid& gradient_image,const bool compute_gradient, const bool isgradient_mu, const float rescale) +{ + + this->output_proj_data_sptr->fill(0.f); + + std::cout << "number of energy windows:= "<< this->template_exam_info_sptr->get_num_energy_windows() << '\n'; + + if(this->template_exam_info_sptr->get_energy_window_pair().first!= -1 && + this->template_exam_info_sptr->get_energy_window_pair().second!= -1 ) + { + std::cout << "energy window pair :="<<" {"<< this->template_exam_info_sptr->get_energy_window_pair().first << ',' << this->template_exam_info_sptr->get_energy_window_pair().second <<"}\n"; + } + + + for (int i = 0; i < this->template_exam_info_sptr->get_num_energy_windows(); ++i) + { + std::cout << "energy window lower level"<<"["<remove_cache_for_integrals_over_activity(); + this->remove_cache_for_integrals_over_attenuation(); + this->sample_scatter_points(); + this->initialise_cache_for_scattpoint_det_integrals_over_attenuation(); + this->initialise_cache_for_scattpoint_det_integrals_over_activity(); + + ViewSegmentNumbers vs_num; + + int bin_counter = 0; + int axial_bins = 0 ; + double sum = 0; + + for (vs_num.segment_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_segment_num(); + vs_num.segment_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_segment_num(); + ++vs_num.segment_num()) + axial_bins += this->proj_data_info_cyl_noarc_cor_sptr->get_num_axial_poss(vs_num.segment_num()); + + const int total_bins = + this->proj_data_info_cyl_noarc_cor_sptr->get_num_views() * axial_bins * + this->proj_data_info_cyl_noarc_cor_sptr->get_num_tangential_poss(); + /* Currently, proj_data_info.find_cartesian_coordinates_of_detection() returns + coordinate in a coordinate system where z=0 in the first ring of the scanner. + We want to shift this to a coordinate system where z=0 in the middle + of the scanner. + We can use get_m() as that uses the 'middle of the scanner' system. + (sorry) + */ +#ifndef NDEBUG + { + CartesianCoordinate3D detector_coord_A, detector_coord_B; + // check above statement + this->proj_data_info_cyl_noarc_cor_sptr->find_cartesian_coordinates_of_detection( + detector_coord_A, detector_coord_B, Bin(0, 0, 0, 0)); + assert(detector_coord_A.z() == 0); + assert(detector_coord_B.z() == 0); + // check that get_m refers to the middle of the scanner + const float m_first = + this->proj_data_info_cyl_noarc_cor_sptr->get_m(Bin(0, 0, this->proj_data_info_cyl_noarc_cor_sptr->get_min_axial_pos_num(0), 0)); + const float m_last = + this->proj_data_info_cyl_noarc_cor_sptr->get_m(Bin(0, 0, this->proj_data_info_cyl_noarc_cor_sptr->get_max_axial_pos_num(0), 0)); + assert(fabs(m_last + m_first) < m_last * 10E-4); + } +#endif + this->shift_detector_coordinates_to_origin = + CartesianCoordinate3D(this->proj_data_info_cyl_noarc_cor_sptr->get_m(Bin(0, 0, 0, 0)), 0, 0); + + info("ScatterSimulator: Initialization finished ..."); + + for (vs_num.segment_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_segment_num(); + vs_num.segment_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_segment_num(); + ++vs_num.segment_num()) + { + for (vs_num.view_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_view_num(); + vs_num.view_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_view_num(); + ++vs_num.view_num()) + { + //info(boost::format("ScatterSimulator: %d / %d") % bin_counter% total_bins); + + + sum+=this->L_G_for_view_segment_number(data, add_sino,gradient_image,vs_num,rescale,compute_gradient,isgradient_mu); + + bin_counter += + this->proj_data_info_cyl_noarc_cor_sptr->get_num_axial_poss(vs_num.segment_num()) * + this->proj_data_info_cyl_noarc_cor_sptr->get_num_tangential_poss(); + //info(boost::format("ScatterSimulator: %d / %d") % bin_counter% total_bins); + + std::cout<< bin_counter << " / "<< total_bins < 1.E20) + warning("KL large at a=%g b=%g, threshold %g\n",a,b,threshold_a); +#ifdef ICHANGEDIT +#undef NDEBUG +#endif +#endif + assert(res>=-1.e-4); + return res; +} + + +double +SingleScatterLikelihoodAndGradient:: +L_G_for_viewgram(const Viewgram& viewgram, const Viewgram& v_add,Viewgram& v_est,VoxelsOnCartesianGrid& gradient_image,const float rescale,const bool compute_gradient, const bool isgradient_mu) +{ + + const ViewSegmentNumbers vs_num(viewgram.get_view_num(),viewgram.get_segment_num()); + + // First construct a vector of all bins that we'll process. + // The reason for making this list before the actual calculation is that we can then parallelise over all bins + // without having to think about double loops. + std::vector all_bins; + { + Bin bin(vs_num.segment_num(), vs_num.view_num(), 0, 0); + + for (bin.axial_pos_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_axial_pos_num(bin.segment_num()); + bin.axial_pos_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_axial_pos_num(bin.segment_num()); + ++bin.axial_pos_num()) + { + for (bin.tangential_pos_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_tangential_pos_num(); + bin.tangential_pos_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_tangential_pos_num(); + ++bin.tangential_pos_num()) + { + all_bins.push_back(bin); + } + } + } + // now compute scatter for all bins + + double sum = 0; + + + VoxelsOnCartesianGrid tmp_gradient_image(gradient_image); + + for (int i = 0; i < static_cast(all_bins.size()); ++i) + { + //creates a template image to fill + tmp_gradient_image.fill(0); + + const Bin bin = all_bins[i]; + + //forward model + const double y = L_G_estimate(tmp_gradient_image,bin,compute_gradient,isgradient_mu); + + v_est[bin.axial_pos_num()][bin.tangential_pos_num()] = static_cast(rescale*y); + //in case a scaling factor for the data is needed,i.e. for adding different level of noise. By default is set to 1. + + float eps = v_add[bin.axial_pos_num()][bin.tangential_pos_num()]; + sum+=viewgram[bin.axial_pos_num()][bin.tangential_pos_num()]*log(v_est[bin.axial_pos_num()][bin.tangential_pos_num()]+eps)- v_est[bin.axial_pos_num()][bin.tangential_pos_num()]-eps; + gradient_image += tmp_gradient_image*(viewgram[bin.axial_pos_num()][bin.tangential_pos_num()]/(v_est[bin.axial_pos_num()][bin.tangential_pos_num()]+eps)-1); + } + + return sum; +} + + + + +double +SingleScatterLikelihoodAndGradient:: +L_G_estimate(VoxelsOnCartesianGrid& gradient_image_bin,const Bin bin, const bool compute_gradient,const bool isgradient_mu) +{ + double scatter_ratio_singles = 0; + unsigned det_num_B=0; + unsigned det_num_A=0; + + this->find_detectors(det_num_A, det_num_B,bin); + + for(std::size_t scatter_point_num =0; + scatter_point_num < this->scatt_points_vector.size(); + ++scatter_point_num) + { + scatter_ratio_singles += + L_G_for_one_scatter_point(gradient_image_bin, + scatter_point_num, + det_num_A, det_num_B,compute_gradient, isgradient_mu); + } + + return scatter_ratio_singles; +} + + + +float +SingleScatterLikelihoodAndGradient:: +L_G_for_one_scatter_point(VoxelsOnCartesianGrid& gradient, + const std::size_t scatter_point_num, + const unsigned det_num_A, + const unsigned det_num_B, const bool compute_gradient,const bool isgradient_mu) +{ + + // The code now supports more than one energy window: the low energy threshold has to correspond to lowest window. + + int low = 0; + + if (this->template_exam_info_sptr->get_num_energy_windows()>1) + + { + + int first_window=this->template_exam_info_sptr->get_energy_window_pair().first-1; + int second_window=this->template_exam_info_sptr->get_energy_window_pair().second-1; + + if(this->template_exam_info_sptr->get_low_energy_thres(first_window) <= this->template_exam_info_sptr->get_low_energy_thres(second_window) ) + + { + low = first_window; + } + + else if(this->template_exam_info_sptr->get_low_energy_thres(first_window) >= this->template_exam_info_sptr->get_low_energy_thres(second_window) ) + + { + low = second_window; + } + + } + + static const float max_single_scatter_cos_angle=max_cos_angle(this->template_exam_info_sptr->get_low_energy_thres(low), + 2.f, + this->proj_data_info_cyl_noarc_cor_sptr->get_scanner_ptr()->get_energy_resolution()); + + //static const float min_energy=energy_lower_limit(lower_energy_threshold,2.,energy_resolution); + + const CartesianCoordinate3D& scatter_point = + this->scatt_points_vector[scatter_point_num].coord; + const CartesianCoordinate3D& detector_coord_A = + this->detection_points_vector[det_num_A]; + const CartesianCoordinate3D& detector_coord_B = + this->detection_points_vector[det_num_B]; + // note: costheta is -cos_angle such that it is 1 for zero scatter angle + const float costheta = static_cast( + -cos_angle(detector_coord_A - scatter_point, + detector_coord_B - scatter_point)); + // note: costheta is identical for scatter to A or scatter to B + // Hence, the Compton_cross_section and energy are identical for both cases as well. + if(max_single_scatter_cos_angle>costheta) + return 0; + const float new_energy = + photon_energy_after_Compton_scatter_511keV(costheta); + + // The detection efficiency varies with respect to the energy window. + //The code can now compute the scatter for a combination of two windows X and Y + //Default: one window -> The code will combine the window with itself + + + //compute the probability of detection for two given energy windows X and Y + + + std::vectordetection_efficiency_scattered; + std::vectordetection_efficiency_unscattered; + + + detection_efficiency_scattered.push_back(0); + detection_efficiency_unscattered.push_back(0); + + + + //detection efficiency of each window for that energy + for (int i = 0; i < this->template_exam_info_sptr->get_num_energy_windows(); ++i) + { + detection_efficiency_scattered[i] = detection_efficiency(new_energy,i); + detection_efficiency_unscattered[i] = detection_efficiency(511.F,i); + } + + + int index0 = 0; + int index1 = 0; + + if (this->template_exam_info_sptr->get_num_energy_windows()>1) + { + index0 = this->template_exam_info_sptr->get_energy_window_pair().first-1; + index1 = this->template_exam_info_sptr->get_energy_window_pair().second-1; + + } + + + float detection_probability_XY=detection_efficiency_scattered[index0]*detection_efficiency_unscattered[index1]; + float detection_probability_YX=detection_efficiency_scattered[index1]*detection_efficiency_unscattered[index0]; + + + if ((detection_probability_XY==0)&&(detection_probability_YX==0)) + return 0; + + const float emiss_to_detA = + cached_integral_over_activity_image_between_scattpoint_det + (static_cast (scatter_point_num), + det_num_A); + const float emiss_to_detB = + cached_integral_over_activity_image_between_scattpoint_det + (static_cast (scatter_point_num), + det_num_B); + if (emiss_to_detA==0 && emiss_to_detB==0) + return 0; + const float atten_to_detA = + cached_exp_integral_over_attenuation_image_between_scattpoint_det + (scatter_point_num, + det_num_A); + const float atten_to_detB = + cached_exp_integral_over_attenuation_image_between_scattpoint_det + (scatter_point_num, + det_num_B); + + const float dif_Compton_cross_section_value = + dif_Compton_cross_section(costheta, 511.F); + + const float rA_squared=static_cast(norm_squared(scatter_point-detector_coord_A)); + const float rB_squared=static_cast(norm_squared(scatter_point-detector_coord_B)); + + const float scatter_point_mu= + scatt_points_vector[scatter_point_num].mu_value; + + const CartesianCoordinate3D + detA_to_ring_center(0,-detector_coord_A[2],-detector_coord_A[3]); + const CartesianCoordinate3D + detB_to_ring_center(0,-detector_coord_B[2],-detector_coord_B[3]); + const float cos_incident_angle_AS = static_cast( + cos_angle(scatter_point - detector_coord_A, + detA_to_ring_center)); + const float cos_incident_angle_BS = static_cast( + cos_angle(scatter_point - detector_coord_B, + detB_to_ring_center)); + + if (cos_incident_angle_AS*cos_incident_angle_BS<0) + return 0; + +#ifndef NEWSCALE + /* projectors work in pixel units, so convert attenuation data + from cm^-1 to pixel_units^-1 */ + const float rescale = + dynamic_cast &>(*density_image_sptr). + get_grid_spacing()[3]/10; +#else + const float rescale = + 0.1F; +#endif + //normalisation + // we will divide by the solid angle factors for unscattered photons + // (computed with the same detection model as used in the scatter code) + // the energy dependency is left out + + + const double common_factor = + 1/detection_efficiency_no_scatter(det_num_A, det_num_B) * + scatter_volume/total_Compton_cross_section_511keV; + + + // Single ScatterForward Model + + const double line_integral1 = detection_probability_XY*(1.F/rB_squared)*pow(atten_to_detB,total_Compton_cross_section_relative_to_511keV(new_energy)-1); + const double line_integral2 = detection_probability_YX*(1.F/rA_squared)*pow(atten_to_detA,total_Compton_cross_section_relative_to_511keV(new_energy)-1); + const double line_integral1_times_activityA = line_integral1*emiss_to_detA; + const double line_integral2_times_activityB = line_integral2*emiss_to_detB; + const double global_factor = atten_to_detB*atten_to_detA*cos_incident_angle_AS*cos_incident_angle_BS*dif_Compton_cross_section_value*common_factor; + const double global_factor_times_mu = global_factor*scatter_point_mu; + + float scatter_ratio=0; + + scatter_ratio = (line_integral1_times_activityA+line_integral2_times_activityB)*global_factor_times_mu; + + + /*Single Scatter Forward model Jacobian w.r.t. attenuation: + * The derivative is given by three terms, respectively in [A,S], [B,S] and [S] */ + + float contribution_AS_mu = (line_integral1_times_activityA+line_integral2_times_activityB*total_Compton_cross_section_relative_to_511keV(new_energy)) + *global_factor_times_mu; + + float contribution_BS_mu = (line_integral1_times_activityA*total_Compton_cross_section_relative_to_511keV(new_energy)+line_integral2_times_activityB*pow(atten_to_detA,total_Compton_cross_section_relative_to_511keV(new_energy)-1)) + *global_factor_times_mu; + + float contribution_S = (line_integral1_times_activityA+line_integral2_times_activityB) + *global_factor; + + + /*Single Scatter Forward model Jacobian w.r.t. activity: + * The derivative is given by two terms, respectively in [A,S] and [B,S] */ + + + float contribution_AS_act = line_integral1*global_factor_times_mu; + float contribution_BS_act = line_integral2*global_factor_times_mu; + + //Fill gradient image along [A,S], [B,S] and in [S] + +if(compute_gradient) +{ + if (isgradient_mu) + + { + line_contribution(gradient,rescale,scatter_point,detector_coord_B,contribution_BS_mu); + line_contribution(gradient,rescale,scatter_point, detector_coord_A,contribution_AS_mu); + s_contribution(gradient,scatter_point,contribution_S); + } + else + + { + line_contribution_act(gradient,scatter_point,detector_coord_B,contribution_BS_act); + line_contribution_act(gradient,scatter_point, detector_coord_A,contribution_AS_act); + } + +} + return scatter_ratio; +} + +ProjDataInMemory +SingleScatterLikelihoodAndGradient:: +likelihood_and_gradient_scatter(const ProjData &projdata, const ProjData& norm , const ProjData &add_projdata, VoxelsOnCartesianGrid& gradient_image_HR, VoxelsOnCartesianGrid& gradient_image_LR,const bool compute_gradient, const bool isgradient_mu) +{ + gradient_image_LR.fill(0); + gradient_image_HR.fill(0); + int length = this->output_proj_data_sptr->get_num_views()*this->output_proj_data_sptr->get_num_axial_poss(0)*this->output_proj_data_sptr->get_num_tangential_poss(); //TODO: the code is for segment zero only + std::vector > jacobian_array; + std::vector ratio; + jacobian_array.reserve(length); + ratio.reserve(length); + for (int i = 0 ; i < length ; ++i) + { + jacobian_array.push_back(gradient_image_LR); + ratio.push_back(0); + } + + const ProjDataInMemory est_data_LR = get_jacobian(jacobian_array, compute_gradient, isgradient_mu); + const ProjDataInMemory est_data_HR = get_ratio(projdata,norm,add_projdata,est_data_LR,ratio); + + for (int i = 0 ; i < length ; ++i) + { + gradient_image_LR += jacobian_array[i]*ratio[i]; + } + + if((gradient_image_HR.get_x_size()!=gradient_image_LR.get_x_size())||(gradient_image_HR.get_y_size()!=gradient_image_LR.get_y_size())||(gradient_image_HR.get_z_size()!=gradient_image_LR.get_z_size())) + { + if(isgradient_mu==true) + transpose_zoom_image(gradient_image_HR,gradient_image_LR,ZoomOptions::preserve_values); + else + transpose_zoom_image(gradient_image_HR,gradient_image_LR,ZoomOptions::preserve_projections); + + } + + return est_data_HR; + +} + +ProjDataInMemory +SingleScatterLikelihoodAndGradient:: +get_jacobian(std::vector > &gradient_image_array,const bool compute_gradient, const bool isgradient_mu) +{ + + ProjDataInMemory est_data(this->get_output_proj_data_sptr()->get_exam_info_sptr(),this->get_output_proj_data_sptr()->get_proj_data_info_sptr()); + est_data.fill(0); + + this->remove_cache_for_integrals_over_activity(); + this->remove_cache_for_integrals_over_attenuation(); + this->sample_scatter_points(); + this->initialise_cache_for_scattpoint_det_integrals_over_attenuation(); + this->initialise_cache_for_scattpoint_det_integrals_over_activity(); + + int bin_counter = 0; + int axial_bins = 0 ; + //double sum = 0; TODO : return the likelihood rather than estimated projdata + +// #ifdef STIR_OPENMP +// #pragma omp parallel for reduction(+:axial_bins) schedule(dynamic) +// #endif + + ViewSegmentNumbers vs_num; + for (vs_num.segment_num()= this->proj_data_info_cyl_noarc_cor_sptr->get_min_segment_num(); + vs_num.segment_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_segment_num(); + ++vs_num.segment_num()) + + axial_bins += this->proj_data_info_cyl_noarc_cor_sptr->get_num_axial_poss(vs_num.segment_num()); + + const int total_bins = + this->proj_data_info_cyl_noarc_cor_sptr->get_num_views() * axial_bins * + this->proj_data_info_cyl_noarc_cor_sptr->get_num_tangential_poss(); + /* Currently, proj_data_info.find_cartesian_coordinates_of_detection() returns + coordinate in a coordinate system where z=0 in the first ring of the scanner. + We want to shift this to a coordinate system where z=0 in the middle + of the scanner. + We can use get_m() as that uses the 'middle of the scanner' system. + (sorry) + */ +#ifndef NDEBUG + { + CartesianCoordinate3D detector_coord_A, detector_coord_B; + // check above statement + this->proj_data_info_cyl_noarc_cor_sptr->find_cartesian_coordinates_of_detection( + detector_coord_A, detector_coord_B, Bin(0, 0, 0, 0)); + assert(detector_coord_A.z() == 0); + assert(detector_coord_B.z() == 0); + // check that get_m refers to the middle of the scanner + const float m_first = + this->proj_data_info_cyl_noarc_cor_sptr->get_m(Bin(0, 0, this->proj_data_info_cyl_noarc_cor_sptr->get_min_axial_pos_num(0), 0)); + const float m_last = + this->proj_data_info_cyl_noarc_cor_sptr->get_m(Bin(0, 0, this->proj_data_info_cyl_noarc_cor_sptr->get_max_axial_pos_num(0), 0)); + assert(fabs(m_last + m_first) < m_last * 10E-4); + } +#endif + this->shift_detector_coordinates_to_origin = + CartesianCoordinate3D(this->proj_data_info_cyl_noarc_cor_sptr->get_m(Bin(0, 0, 0, 0)), 0, 0); + + info("ScatterSimulator: Initialization finished ..."); + + for (vs_num.segment_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_segment_num(); + vs_num.segment_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_segment_num(); + ++vs_num.segment_num()) + { + for (vs_num.view_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_view_num(); + vs_num.view_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_view_num(); + ++vs_num.view_num()) + { + //info(boost::format("ScatterSimulator: %d / %d") % bin_counter% total_bins); + + + this->get_jacobian_for_view_segment_number(gradient_image_array,est_data,vs_num,compute_gradient,isgradient_mu); + + bin_counter += + this->proj_data_info_cyl_noarc_cor_sptr->get_num_axial_poss(vs_num.segment_num()) * + this->proj_data_info_cyl_noarc_cor_sptr->get_num_tangential_poss(); + //info(boost::format("ScatterSimulator: %d / %d") % bin_counter% total_bins); + + std::cout<< bin_counter << " / "<< total_bins < > &gradient_image_array, ProjData &est_data, const ViewSegmentNumbers& vs_num, const bool compute_gradient,const bool isgradient_mu) +{ + + Viewgram v_est = est_data.get_empty_viewgram(vs_num.view_num(), vs_num.segment_num()); + get_jacobian_for_viewgram(v_est,gradient_image_array, compute_gradient,isgradient_mu); + est_data.set_viewgram(v_est); + +} + +void +SingleScatterLikelihoodAndGradient:: +get_jacobian_for_viewgram(Viewgram& v_est,std::vector > &gradient_image_array,const bool compute_gradient, const bool isgradient_mu) +{ + + const ViewSegmentNumbers vs_num(v_est.get_view_num(),v_est.get_segment_num()); + + // First construct a vector of all bins that we'll process. + // The reason for making this list before the actual calculation is that we can then parallelise over all bins + // without having to think about double loops. + std::vector all_bins; + { + Bin bin(vs_num.segment_num(), vs_num.view_num(), 0, 0); + + for (bin.axial_pos_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_axial_pos_num(bin.segment_num()); + bin.axial_pos_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_axial_pos_num(bin.segment_num()); + ++bin.axial_pos_num()) + { + for (bin.tangential_pos_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_tangential_pos_num(); + bin.tangential_pos_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_tangential_pos_num(); + ++bin.tangential_pos_num()) + { + all_bins.push_back(bin); + } + } + } + // now compute scatter for all bins + if(gradient_image_array.size()!=static_cast(all_bins.size())*this->output_proj_data_sptr->get_num_views()) + error("SIZE is %d , but it should be %d",gradient_image_array.size(),static_cast(all_bins.size())); + + + #ifdef STIR_OPENMP + #pragma omp parallel for schedule(dynamic) + #endif + for (int i = 0; i < static_cast(all_bins.size()); ++i) + { // this needs to be defined inside here to be thread-safe + VoxelsOnCartesianGrid tmp_gradient_image(gradient_image_array[0]); + tmp_gradient_image.fill(0); + + const Bin bin = all_bins[i]; + const double y = L_G_estimate(tmp_gradient_image,bin,compute_gradient,isgradient_mu); + + v_est[bin.axial_pos_num()][bin.tangential_pos_num()] = static_cast(y); //this is passed as reference and filled in the loop + gradient_image_array[i] += tmp_gradient_image; + } + +} + +ProjDataInMemory +SingleScatterLikelihoodAndGradient:: +get_ratio(const ProjData& projdata,const ProjData& norm,const ProjData &add_projdata, const ProjData &est_projdata, std::vector &ratio_vector) +{ + + ProjDataInMemory ratio_HR(projdata); ratio_HR.fill(0); + ProjDataInMemory est_projdata_HR(projdata); est_projdata_HR.fill(0); + ProjDataInMemory ratio_LR(est_projdata); ratio_LR.fill(0); + + if((est_projdata.get_num_views()!=projdata.get_num_views())||(est_projdata.get_num_tangential_poss()!=projdata.get_num_tangential_poss())) + ScatterEstimation::pull_scatter_estimate(est_projdata_HR,projdata,est_projdata,norm,true); + else + est_projdata_HR.fill(est_projdata); + + Bin bin; + { + for (bin.segment_num()=projdata.get_min_segment_num(); bin.segment_num()<=projdata.get_max_segment_num(); ++bin.segment_num()) + for (bin.axial_pos_num() = projdata.get_min_axial_pos_num(bin.segment_num()); bin.axial_pos_num()<=projdata.get_max_axial_pos_num(bin.segment_num()); ++bin.axial_pos_num()) + { + Sinogram sino = projdata.get_sinogram(bin.axial_pos_num(),bin.segment_num()); + Sinogram add_sino = add_projdata.get_sinogram(bin.axial_pos_num(),bin.segment_num()); + Sinogram est_sino = est_projdata_HR.get_sinogram(bin.axial_pos_num(),bin.segment_num()); + Sinogram ratio_sino = ratio_HR.get_empty_sinogram(bin.axial_pos_num(),bin.segment_num()); + + for (bin.view_num()=sino.get_min_view_num(); + bin.view_num()<=sino.get_max_view_num(); + ++bin.view_num()) + { + for (bin.tangential_pos_num()= sino.get_min_tangential_pos_num(); bin.tangential_pos_num()<= sino.get_max_tangential_pos_num(); ++bin.tangential_pos_num()) + if(est_sino[bin.axial_pos_num()][bin.tangential_pos_num()]==0) + ratio_sino[bin.view_num()][bin.tangential_pos_num()] = 0; + else + ratio_sino[bin.view_num()][bin.tangential_pos_num()] = sino[bin.axial_pos_num()][bin.tangential_pos_num()]/(est_sino[bin.axial_pos_num()][bin.tangential_pos_num()]+add_sino[bin.axial_pos_num()][bin.tangential_pos_num()])-1; + ratio_HR.set_sinogram(ratio_sino); + + } + + } + } + + if((est_projdata.get_num_views()!=projdata.get_num_views())||(est_projdata.get_num_tangential_poss()!=projdata.get_num_tangential_poss())) + ScatterEstimation::push_scatter_estimate(ratio_LR,est_projdata,ratio_HR,norm,true); + else + ratio_LR.fill(ratio_HR); + + ViewSegmentNumbers vs_num; + int counter = 0; + + for (vs_num.segment_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_segment_num(); vs_num.segment_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_segment_num(); ++vs_num.segment_num()) + { + for (vs_num.view_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_view_num();vs_num.view_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_view_num(); ++vs_num.view_num()) + { + + Viewgram viewgram = ratio_LR.get_viewgram(vs_num.view_num(), vs_num.segment_num(),false); + const ViewSegmentNumbers vs_num(viewgram.get_view_num(),viewgram.get_segment_num()); + std::vector all_bins; + { + Bin bin(vs_num.segment_num(), vs_num.view_num(), 0, 0); + for (bin.axial_pos_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_axial_pos_num(bin.segment_num()); bin.axial_pos_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_axial_pos_num(bin.segment_num()); ++bin.axial_pos_num()) + { + for (bin.tangential_pos_num() = this->proj_data_info_cyl_noarc_cor_sptr->get_min_tangential_pos_num(); bin.tangential_pos_num() <= this->proj_data_info_cyl_noarc_cor_sptr->get_max_tangential_pos_num(); ++bin.tangential_pos_num()) + { + all_bins.push_back(bin); + } + } + } + + for (int i = 0; i < static_cast(all_bins.size()); ++i) + { + ++ counter; + + const Bin bin = all_bins[i]; + + ratio_vector[i] = viewgram[bin.axial_pos_num()][bin.tangential_pos_num()]; + } + + } + } + + if(ratio_vector.size()!=counter) + error("SIZE is %d , but it should be %d",ratio_vector.size(),counter); +return est_projdata_HR;// +} +END_NAMESPACE_STIR + diff --git a/src/scatter_buildblock/scatter_detection_modelling.cxx b/src/scatter_buildblock/scatter_detection_modelling.cxx index 25edc73f08..05ff966ad1 100644 --- a/src/scatter_buildblock/scatter_detection_modelling.cxx +++ b/src/scatter_buildblock/scatter_detection_modelling.cxx @@ -25,7 +25,11 @@ #include "stir/ProjDataInfoBlocksOnCylindricalNoArcCorr.h" #include "stir/numerics/erf.h" #include "stir/info.h" +#include "stir/CPUTimer.h" #include +#ifdef STIR_OPENMP +#include +#endif START_NAMESPACE_STIR unsigned @@ -110,7 +114,7 @@ compute_emis_to_det_points_solid_angle_factor( float ScatterSimulation:: -detection_efficiency(const float energy) const +detection_efficiency(const float energy, const int en_window) const { #ifndef NDEBUG if (!this->_already_set_up) @@ -125,12 +129,132 @@ detection_efficiency(const float energy) const // sigma_times_sqrt2= sqrt(2) * sigma // resolution proportional to FWHM const float efficiency = - 0.5f*( erf((this->template_exam_info_sptr->get_high_energy_thres()-energy)/sigma_times_sqrt2) - - erf((this->template_exam_info_sptr->get_low_energy_thres()-energy)/sigma_times_sqrt2 )); + 0.5f*( erf((this->template_exam_info_sptr->get_high_energy_thres(en_window)-energy)/sigma_times_sqrt2) + - erf((this->template_exam_info_sptr->get_low_energy_thres(en_window)-energy)/sigma_times_sqrt2 )); /* Maximum efficiency is 1.*/ + return efficiency; } +float +ScatterSimulation::detection_efficiency_numerical_formulation(const float incoming_photon_energy, const int en_window) const +{ + + const float HLD = this->template_exam_info_sptr->get_high_energy_thres(en_window); + const float LLD = this->template_exam_info_sptr->get_low_energy_thres(en_window); + float sum = 0; + const int size = 30; + double increment_x = (HLD - LLD) / (size - 1); + + #ifdef STIR_OPENMP + #pragma omp parallel for reduction(+:sum) schedule(dynamic) + #endif + for(int i = 0 ; i< size; ++i) + { + const float energy_range = LLD+i*increment_x; + sum+= detection_model_with_fitted_parameters(energy_range, incoming_photon_energy); + } + sum*=increment_x; + return sum; +} + +std::vector +ScatterSimulation::detection_spectrum_numerical_formulation(const float LLD, const float HLD, const float size, const float incoming_photon_energy) const +{ + + std::vector output(size); + double increment_x = (HLD - LLD) / (size - 1); + + for(int i = 0; i < size; ++i) + { + output[i] = LLD + (i * increment_x); + std::cout << output[i] << std::endl; + } + + for(int i = 0; i < size; ++i) + { + output[i] = detection_model_with_fitted_parameters(output[i], incoming_photon_energy); + } + + return output; +} + + +float +ScatterSimulation:: +detection_model_with_fitted_parameters(const float x, const float energy) const +{ + //! Brief + //! All the parameters are obtained by fitting the model to the energy spectrum obtained with GATE. + //! The crystal used here is LSO, the one for the Siemens mMR (atomic number Z = 66). + //! We consider to have four terms: (i) gaussian model for the photopeak, (ii) compton plateau, (iii) flat continuum, (iv) exponential tale + //! The model has be trained with 511 keV and tested with 370 keV. + + //const int Z = 66; // atomic number of LSO + //const float H_1 = pow(Z,5)/energy; //the height of the photopeak is prop. to the photoelectric cross section + //const float H_2 = 9.33*pow(10,25)*total_Compton_cross_section(energy)*Z; // the eight of the compton plateau is proportional to the compton cross section + //const float H_3 = 7; //fitting parameter + //const float H_4 = 29.4; //fitting parameter + //const float beta = -0.8401; //fitting parameter + //const float global_scale = 0.29246*0.8*1e-06;//2.33*1e-07; //fitting parameter + //const float fwhm = this->proj_data_info_cyl_noarc_cor_sptr->get_scanner_ptr()->get_energy_resolution(); + //const float fwhm = 0.14; + //const float std_peak = energy*fwhm/2.35482; + //const float scaling_std_compton = 28.3; //fitting parameter + //const float shift_compton = 0.597; //fitting parameter + const float f1 = photoelectric((66*66*66*66*66)/energy, (energy*0.14f)/2.35482f, x, energy); + const float f2 = compton_plateau(9.33f*1e+25*total_Compton_cross_section(energy)*66, (energy*0.14f)/2.35482f, x, energy, 28.3f,0.597f); + const float f3 = flat_continuum(7,(energy*0.14f)/2.35482f, x, energy); + const float f4 = exponential_tail(29.4f,(energy*0.14f)/2.35482f, x, energy,-0.8401f); + + return 0.29246f*0.8f*1e-06*(f1+f2+f3+f4); +} + +float +ScatterSimulation:: +photoelectric(const float K, const float std_peak, const float x, const float energy) const +{ + const float diff = x - energy; + const float pow_diff = diff*diff; + const float pow_std_peak = std_peak*std_peak; + return K/(std_peak*2.5066f)*exp(-pow_diff/(2*pow_std_peak)); +} + +float +ScatterSimulation:: +compton_plateau(const float K, const float std_peak, const float x, const float energy, const float scaling_std_compton,const float shift_compton) const +{ + const float m_0_c_2 = 511.0f; + const float alpha = energy/m_0_c_2; + const float E_1 = energy/(1+alpha*(2)); + const float mean = energy*shift_compton; + const float x_minus_mean = x - mean; + return ((energy/E_1)+(E_1/energy)-2)*(K*exp(-(x_minus_mean*x_minus_mean)/(4*scaling_std_compton*std_peak))); +} +float +ScatterSimulation:: +flat_continuum(const float K, const float std_peak, const float x, const float energy) const +{ + const float den = 1.4142f*std_peak; + if (x<=energy) + return K* erfc((x-energy)/den); + else + return 0; +} + +float +ScatterSimulation:: +exponential_tail(const float K, const float std_peak, const float x, const float energy, const float beta) const +{ + const float den1 = 1.4142f*_PI*std_peak*beta; + const float den2 = 1.4142f*std_peak; + const float den3 = 2*beta; + if (x > 210) //i am not sure of the behaviour of the function at too low energies + return K * exp((x-energy)/den1)*erfc((x-energy)/den2+1/den3); + else + return 0; +} + float ScatterSimulation:: max_cos_angle(const float low, const float approx, const float resolution_at_511keV) @@ -154,18 +278,7 @@ ScatterSimulation:: detection_efficiency_no_scatter(const unsigned det_num_A, const unsigned det_num_B) const { -#ifndef NDEBUG - if (!this->_already_set_up) - error("ScatterSimulation::find_detectors: need to call set_up() first"); -#endif - if (detector_efficiency_no_scatter <= 0.F) // set to negative value by set_up(), so recompute - { - detector_efficiency_no_scatter = - detection_efficiency(511.F) > 0 - ? detection_efficiency(511.F) - : (info("Zero detection efficiency for 511. Will normalise to 1"), 1.F); - } const CartesianCoordinate3D& detector_coord_A = detection_points_vector[det_num_A]; const CartesianCoordinate3D& detector_coord_B = @@ -186,7 +299,7 @@ detection_efficiency_no_scatter(const unsigned det_num_A, return 1./( 0.75/2./_PI * rAB_squared - /detector_efficiency_no_scatter/ + /pow(1,2.0)/ (cos_incident_angle_A* cos_incident_angle_B)); } diff --git a/src/scatter_buildblock/scatter_estimate_for_one_scatter_point.cxx b/src/scatter_buildblock/scatter_estimate_for_one_scatter_point.cxx index 3345e4c63b..67c588af35 100644 --- a/src/scatter_buildblock/scatter_estimate_for_one_scatter_point.cxx +++ b/src/scatter_buildblock/scatter_estimate_for_one_scatter_point.cxx @@ -42,14 +42,44 @@ SingleScatterSimulation:: const std::size_t scatter_point_num, const unsigned det_num_A, const unsigned det_num_B) -{ - if (this->max_single_scatter_cos_angle <= 0.F) // set to negative value by set_up(), so recompute +{ + + + + + // The code now supports more than one energy window: the low energy threshold has to correspond to lowest window. + + int low = 0; + + // TODO: check that the selected windows don't overcome the max num windows + + if (this->template_exam_info_sptr->get_num_energy_windows()>1) + { - this->max_single_scatter_cos_angle=max_cos_angle(this->template_exam_info_sptr->get_low_energy_thres(), - 2.f, - this->proj_data_info_cyl_noarc_cor_sptr->get_scanner_ptr()->get_energy_resolution()); + + int first_window=this->template_exam_info_sptr->get_energy_window_pair().first-1; + int second_window=this->template_exam_info_sptr->get_energy_window_pair().second-1; + + if(this->template_exam_info_sptr->get_low_energy_thres(first_window) <= this->template_exam_info_sptr->get_low_energy_thres(second_window) ) + + { + low = first_window; + } + else if(this->template_exam_info_sptr->get_low_energy_thres(first_window) >= this->template_exam_info_sptr->get_low_energy_thres(second_window) ) + + { + low = second_window; + } + //std::cerr << "low:= "<< low<< '\n'; + } + + + static const float max_single_scatter_cos_angle=max_cos_angle(this->template_exam_info_sptr->get_low_energy_thres(low), + 2.f, + this->proj_data_info_cyl_noarc_cor_sptr->get_scanner_ptr()->get_energy_resolution()); + //static const float min_energy=energy_lower_limit(lower_energy_threshold,2.,energy_resolution); const CartesianCoordinate3D& scatter_point = @@ -64,15 +94,51 @@ SingleScatterSimulation:: detector_coord_B - scatter_point)); // note: costheta is identical for scatter to A or scatter to B // Hence, the Compton_cross_section and energy are identical for both cases as well. - if(this->max_single_scatter_cos_angle>costheta) - return 0; + if(max_single_scatter_cos_angle>costheta) + return 0; const float new_energy = photon_energy_after_Compton_scatter_511keV(costheta); - const float detection_efficiency_scatter = - detection_efficiency(new_energy); - if (detection_efficiency_scatter==0) - return 0; + + + +// The detection efficiency varies with respect to the energy window. +//The code can now compute the scatter for a combination of two windows X and Y +//Default: one window -> The code will combine the window with itself + +std::vectordetection_efficiency_scattered; +std::vectordetection_efficiency_unscattered; + + +detection_efficiency_scattered.push_back(0); +detection_efficiency_unscattered.push_back(0); + + + +//detection efficiency of each window for that energy + for (int i = 0; i < this->template_exam_info_sptr->get_num_energy_windows(); ++i) + { + detection_efficiency_scattered[i] = detection_efficiency(new_energy,i); + detection_efficiency_unscattered[i] = detection_efficiency(511.F,i); + } + + + int index0 = 0; + int index1 = 0; + + if (this->template_exam_info_sptr->get_num_energy_windows()>1) + { + index0 = this->template_exam_info_sptr->get_energy_window_pair().first-1; + index1 = this->template_exam_info_sptr->get_energy_window_pair().second-1; + + } + + float detection_probability_XY=detection_efficiency_scattered[index0]*detection_efficiency_unscattered[index1]; + float detection_probability_YX=detection_efficiency_scattered[index1]*detection_efficiency_unscattered[index0]; + + if ((detection_probability_XY==0)&&(detection_probability_YX==0)) + return 0; + const float emiss_to_detA = cached_integral_over_activity_image_between_scattpoint_det @@ -102,6 +168,19 @@ SingleScatterSimulation:: const float scatter_point_mu= scatt_points_vector[scatter_point_num].mu_value; + const CartesianCoordinate3D + detA_to_ring_center(0,-detector_coord_A[2],-detector_coord_A[3]); + const CartesianCoordinate3D + detB_to_ring_center(0,-detector_coord_B[2],-detector_coord_B[3]); + const float cos_incident_angle_AS = static_cast( + cos_angle(scatter_point - detector_coord_A, + detA_to_ring_center)) ; + + const float cos_incident_angle_BS = static_cast( + cos_angle(scatter_point - detector_coord_B, + detB_to_ring_center)) ; + if (cos_incident_angle_AS*cos_incident_angle_BS<0) + return 0; #ifndef NDEBUG { // check if mu-value ok @@ -118,29 +197,33 @@ SingleScatterSimulation:: } #endif + //normalisation + // we will divide by the solid angle factors for unscattered photons + // (computed with the same detection model as used in the scatter code) + // the energy dependency is left out + + const double common_factor = + 1/detection_efficiency_no_scatter(det_num_A, det_num_B) * + scatter_volume/total_Compton_cross_section_511keV; + + double scatter_ratio=0 ; - scatter_ratio= - (emiss_to_detA*(1./rB_squared)*pow(atten_to_detB,total_Compton_cross_section_relative_to_511keV(new_energy)-1) - +emiss_to_detB*(1./rA_squared)*pow(atten_to_detA,total_Compton_cross_section_relative_to_511keV(new_energy)-1)) - *atten_to_detB - *atten_to_detA - *scatter_point_mu - *detection_efficiency_scatter; - + scatter_ratio= + (detection_probability_XY*emiss_to_detA*(1.F/rB_squared)*pow(atten_to_detB,total_Compton_cross_section_relative_to_511keV(new_energy)-1) + +detection_probability_YX*emiss_to_detB*(1.F/rA_squared)*pow(atten_to_detA,total_Compton_cross_section_relative_to_511keV(new_energy)-1)) + *atten_to_detB + *atten_to_detA + *scatter_point_mu + *cos_incident_angle_AS + *cos_incident_angle_BS + *dif_Compton_cross_section_value + *common_factor; + + + return scatter_ratio; + - const CartesianCoordinate3D - detA_to_ring_center(0,-detector_coord_A[2],-detector_coord_A[3]); - const CartesianCoordinate3D - detB_to_ring_center(0,-detector_coord_B[2],-detector_coord_B[3]); - const float cos_incident_angle_AS = static_cast( - cos_angle(scatter_point - detector_coord_A, - detA_to_ring_center)) ; - const float cos_incident_angle_BS = static_cast( - cos_angle(scatter_point - detector_coord_B, - detB_to_ring_center)) ; - - return scatter_ratio*cos_incident_angle_AS*cos_incident_angle_BS*dif_Compton_cross_section_value; } diff --git a/src/scatter_buildblock/single_scatter_estimate.cxx b/src/scatter_buildblock/single_scatter_estimate.cxx index 41e3ba5719..6be449cc47 100644 --- a/src/scatter_buildblock/single_scatter_estimate.cxx +++ b/src/scatter_buildblock/single_scatter_estimate.cxx @@ -19,9 +19,7 @@ */ #include "stir/scatter/SingleScatterSimulation.h" START_NAMESPACE_STIR -static const float total_Compton_cross_section_511keV = -ScatterSimulation:: - total_Compton_cross_section(511.F); + double SingleScatterSimulation:: @@ -60,20 +58,6 @@ actual_scatter_estimate(double& scatter_ratio_singles, det_num_A, det_num_B); } - - // we will divide by the effiency of the detector pair for unscattered photons - // (computed with the same detection model as used in the scatter code) - // This way, the scatter estimate will correspond to a 'normalised' scatter estimate. - - // there is a scatter_volume factor for every scatter point, as the sum over scatter points - // is an approximation for the integral over the scatter point. - - // the factors total_Compton_cross_section_511keV should probably be moved to the scatter_computation code - const double common_factor = - 1/detection_efficiency_no_scatter(det_num_A, det_num_B) * - scatter_volume/total_Compton_cross_section_511keV; - - scatter_ratio_singles *= common_factor; } END_NAMESPACE_STIR diff --git a/src/scatter_buildblock/single_scatter_integrals.cxx b/src/scatter_buildblock/single_scatter_integrals.cxx index d270b2f287..25574925e3 100644 --- a/src/scatter_buildblock/single_scatter_integrals.cxx +++ b/src/scatter_buildblock/single_scatter_integrals.cxx @@ -23,6 +23,7 @@ */ #include "stir/scatter/ScatterSimulation.h" +#include "stir/scatter/SingleScatterLikelihoodAndGradient.h" #include "stir/VoxelsOnCartesianGrid.h" #include "stir/recon_buildblock/ProjMatrixElemsForOneBin.h" #include "stir/recon_buildblock/RayTraceVoxelsOnCartesianGrid.h" @@ -133,6 +134,137 @@ integral_between_2_points(const DiscretisedDensity<3,float>& density, } } return sum; -} +} + + + +void +SingleScatterLikelihoodAndGradient:: +line_contribution(VoxelsOnCartesianGrid& gradient_image,const float rescale, + const CartesianCoordinate3D& scatter_point, + const CartesianCoordinate3D& detector_coord, + const float C) +{ + + + const DiscretisedDensity<3,float>& density=*density_image_sptr; + const VoxelsOnCartesianGrid& image = + dynamic_cast& > + (density); + + const CartesianCoordinate3D voxel_size = image.get_grid_spacing(); + + CartesianCoordinate3D origin = + image.get_origin(); + const float z_to_middle = + (image.get_max_index() + image.get_min_index())*voxel_size.z()/2.F; + origin.z() -= z_to_middle; + /* TODO replace with image.get_index_coordinates_for_physical_coordinates */ + ProjMatrixElemsForOneBin lor; + RayTraceVoxelsOnCartesianGrid(lor, + (scatter_point-origin)/voxel_size, // should be in voxel units + (detector_coord-origin)/voxel_size, // should be in voxel units + voxel_size, //should be in mm +#ifdef NEWSCALE + 1.F // normalise to mm +#else + 1/voxel_size.x() // normalise to some kind of 'pixel units' +#endif + ); + lor.sort(); + + + //VoxelsOnCartesianGrid gradient_image(image); + //gradient_image.fill(0.F); + { + ProjMatrixElemsForOneBin::iterator element_ptr =lor.begin() ; + bool we_have_been_within_the_image = false; + while (element_ptr != lor.end()) + { + const BasicCoordinate<3,int> coords = element_ptr->get_coords(); + if (coords[1] >= image.get_min_index() && + coords[1] <= image.get_max_index() && + coords[2] >= image[coords[1]].get_min_index() && + coords[2] <= image[coords[1]].get_max_index() && + coords[3] >= image[coords[1]][coords[2]].get_min_index() && + coords[3] <= image[coords[1]][coords[2]].get_max_index()) + { + we_have_been_within_the_image = true; + //gradient_image[coords] += -C*element_ptr->get_value() + gradient_image[coords] += -C*rescale*element_ptr->get_value(); + } + else if (we_have_been_within_the_image) + { + // we jump out of the loop as we are now at the other side of + // the image + // break; + } + ++element_ptr; + } + } + +} + + + + + +void +SingleScatterLikelihoodAndGradient:: +s_contribution(VoxelsOnCartesianGrid& gradient_image, +const CartesianCoordinate3D& scatter_point, + const float D) + +{ + +const DiscretisedDensity<3,float>& density=*density_image_sptr; +const VoxelsOnCartesianGrid& image = + dynamic_cast& > + (density); + +const CartesianCoordinate3D voxel_size = image.get_grid_spacing(); + +CartesianCoordinate3D origin = + image.get_origin(); + +const float z_to_middle = + (image.get_max_index() + image.get_min_index())*voxel_size.z()/2.F; + origin.z() -= z_to_middle; + +CartesianCoordinate3D coord = + (scatter_point-origin)/voxel_size; // voxel units + + +const BasicCoordinate<3,int> coords=round(coord); + + +gradient_image[coords] += D; + + + +} + + +void +SingleScatterLikelihoodAndGradient:: +line_contribution_act(VoxelsOnCartesianGrid& gradient_image, + const CartesianCoordinate3D& scatter_point, + const CartesianCoordinate3D& detector_coord, + const float C) +{ + + const CartesianCoordinate3D dist_vector = scatter_point - detector_coord ; + const float dist_sp1_det_squared = norm_squared(dist_vector); + const float solid_angle_factor = + std::min(static_cast(_PI/2), 1.F / dist_sp1_det_squared) ; + + line_contribution(gradient_image,solid_angle_factor, + scatter_point, + detector_coord, + -C); + +} + + END_NAMESPACE_STIR diff --git a/src/scatter_buildblock/upsample_and_fit_scatter_estimate.cxx b/src/scatter_buildblock/upsample_and_fit_scatter_estimate.cxx index aae648189f..4780c8f869 100644 --- a/src/scatter_buildblock/upsample_and_fit_scatter_estimate.cxx +++ b/src/scatter_buildblock/upsample_and_fit_scatter_estimate.cxx @@ -19,6 +19,8 @@ #include "stir/ProjDataInfo.h" #include "stir/ExamInfo.h" #include "stir/ProjDataInMemory.h" +#include "stir/ViewSegmentNumbers.h" +#include "stir/Viewgram.h" #include "stir/inverse_SSRB.h" #include "stir/scale_sinograms.h" #include "stir/scatter/ScatterEstimation.h" @@ -55,15 +57,17 @@ upsample_and_fit_scatter_estimate(ProjData& scaled_scatter_proj_data, interpolated_direct_scatter_proj_data_info_sptr(emission_proj_data.get_proj_data_info_sptr()->clone()); interpolated_direct_scatter_proj_data_info_sptr->reduce_segment_range(0,0); - info("upsample_and_fit_scatter_estimate: Interpolating scatter estimate to size of emission data"); - ProjDataInMemory interpolated_direct_scatter(emission_proj_data.get_exam_info_sptr(), + + info("upsample_and_fit_scatter_estimate: Interpolating scatter estimate to size of emission data"); + ProjDataInMemory interpolated_direct_scatter(emission_proj_data.get_exam_info_sptr(), interpolated_direct_scatter_proj_data_info_sptr); - interpolate_projdata(interpolated_direct_scatter, scatter_proj_data, spline_type, remove_interleaving); - const TimeFrameDefinitions& time_frame_defs = - emission_proj_data.get_exam_info_sptr()->time_frame_definitions; + interpolate_projdata(interpolated_direct_scatter, scatter_proj_data, spline_type, remove_interleaving); + + + const TimeFrameDefinitions& time_frame_defs = emission_proj_data.get_exam_info_sptr()->time_frame_definitions; - if (min_scale_factor != 1 || max_scale_factor != 1 || !scatter_normalisation.is_trivial()) + if (min_scale_factor != 1 || max_scale_factor != 1 || !scatter_normalisation.is_trivial()) { ProjDataInMemory interpolated_scatter(emission_proj_data.get_exam_info_sptr(), emission_proj_data.get_proj_data_info_sptr()->create_shared_clone()); @@ -139,4 +143,155 @@ upsample_and_fit_scatter_estimate(ProjData& scaled_scatter_proj_data, } } + +void +ScatterEstimation:: +upsample_scatter_estimate(ProjData& scaled_scatter_proj_data, + const ProjData& emission_proj_data, + const ProjData& scatter_proj_data, + const bool remove_interleaving) +{ + stir::BSpline::BSplineType spline_type = stir::BSpline::linear; + shared_ptr interpolated_direct_scatter_proj_data_info_sptr(emission_proj_data.get_proj_data_info_sptr()->clone()); + interpolated_direct_scatter_proj_data_info_sptr->reduce_segment_range(0,0); + + + info("upsample_and_fit_scatter_estimate: Interpolating scatter estimate to size of emission data"); + ProjDataInMemory interpolated_direct_scatter(emission_proj_data.get_exam_info_sptr(), + interpolated_direct_scatter_proj_data_info_sptr); + + // interpolate projdata + interpolate_projdata(interpolated_direct_scatter, scatter_proj_data, spline_type, remove_interleaving); + + // Perform Inverse Single Slice Rebinning + inverse_SSRB(scaled_scatter_proj_data, interpolated_direct_scatter); + +} + + + +void +ScatterEstimation:: +pull_scatter_estimate(ProjData& scaled_scatter_proj_data, + const ProjData& emission_proj_data, + const ProjData& scatter_proj_data, + const bool remove_interleaving) +{ + + shared_ptr interpolated_direct_scatter_proj_data_info_sptr(emission_proj_data.get_proj_data_info_sptr()->clone()); + interpolated_direct_scatter_proj_data_info_sptr->reduce_segment_range(0,0); //create the output template + + + info("upsample_and_fit_scatter_estimate: Interpolating scatter estimate to size of emission data"); + ProjDataInMemory interpolated_direct_scatter(emission_proj_data.get_exam_info_sptr(), + interpolated_direct_scatter_proj_data_info_sptr); + + // interpolate projdata + interpolate_projdata_pull(interpolated_direct_scatter, scatter_proj_data, remove_interleaving); + + // Perform Inverse Single Slice Rebinning + inverse_SSRB(scaled_scatter_proj_data, interpolated_direct_scatter); + +} + +void +ScatterEstimation:: +push_scatter_estimate(ProjData& scaled_scatter_proj_data, + const ProjData& emission_proj_data, + const ProjData& scatter_proj_data, + const bool remove_interleaving) +{ + + shared_ptr new_input_proj_data_info_sptr(scatter_proj_data.get_proj_data_info_sptr()->clone()); + new_input_proj_data_info_sptr->reduce_segment_range(0,0); //create input template + + ProjDataInMemory new_input(scatter_proj_data.get_exam_info_sptr(),new_input_proj_data_info_sptr); + transpose_inverse_SSRB(new_input, scatter_proj_data); + + + interpolate_projdata_push(scaled_scatter_proj_data, new_input, remove_interleaving); + + +} + + +void +ScatterEstimation:: +pull_scatter_estimate(ProjData& scaled_scatter_proj_data, + const ProjData& emission_proj_data, + const ProjData& scatter_proj_data, + const ProjData& norm, + const bool remove_interleaving) +{ + + shared_ptr interpolated_direct_scatter_proj_data_info_sptr(emission_proj_data.get_proj_data_info_sptr()->clone()); + interpolated_direct_scatter_proj_data_info_sptr->reduce_segment_range(0,0); //create the output template + + + info("upsample_and_fit_scatter_estimate: Interpolating scatter estimate to size of emission data"); + ProjDataInMemory interpolated_direct_scatter(emission_proj_data.get_exam_info_sptr(), + interpolated_direct_scatter_proj_data_info_sptr); + + // interpolate projdata + interpolate_projdata_pull(interpolated_direct_scatter, scatter_proj_data, remove_interleaving); + + // Perform Inverse Single Slice Rebinning + inverse_SSRB(scaled_scatter_proj_data, interpolated_direct_scatter); + + apply_norm(scaled_scatter_proj_data,norm); + +} + +void +ScatterEstimation:: +push_scatter_estimate(ProjData& scaled_scatter_proj_data, + const ProjData& emission_proj_data, + const ProjData& scatter_proj_data, + const ProjData& norm, + const bool remove_interleaving) +{ + + ProjDataInMemory scatter_proj_data_in_memory(scatter_proj_data); + apply_norm(scatter_proj_data_in_memory,norm); + + shared_ptr new_input_proj_data_info_sptr(scatter_proj_data_in_memory.get_proj_data_info_sptr()->clone()); + new_input_proj_data_info_sptr->reduce_segment_range(0,0); //create input template + + ProjDataInMemory new_input(scatter_proj_data_in_memory.get_exam_info_sptr(),new_input_proj_data_info_sptr); + + transpose_inverse_SSRB(new_input, scatter_proj_data_in_memory); + + + interpolate_projdata_push(scaled_scatter_proj_data, new_input, remove_interleaving); + + +} + +void +ScatterEstimation:: +apply_norm(ProjData& projdata,const ProjData& norm) +{ + +if((projdata.get_num_views()!=norm.get_num_views())||(projdata.get_num_tangential_poss()!=norm.get_num_tangential_poss())) + error("sinograms have to have the same dimensions"); +ProjDataInMemory projdata_out(projdata); +projdata_out.fill(0); + +ViewSegmentNumbers vs_num; + +for (vs_num.segment_num() = norm.get_min_segment_num(); vs_num.segment_num() <= norm.get_max_segment_num(); ++vs_num.segment_num()) +{ + for (vs_num.view_num() = norm.get_min_view_num();vs_num.view_num() <= norm.get_max_view_num(); ++vs_num.view_num()) + { + + Viewgram viewgram_n = norm.get_viewgram(vs_num.view_num(), vs_num.segment_num()); + Viewgram viewgram_in = projdata.get_viewgram(vs_num.view_num(), vs_num.segment_num()); + viewgram_in *= viewgram_n; + projdata_out.set_viewgram(viewgram_in); + } +} + +projdata.fill(projdata_out); + +} END_NAMESPACE_STIR diff --git a/src/scatter_utilities/CMakeLists.txt b/src/scatter_utilities/CMakeLists.txt index 617721bc10..6793d1f6e4 100644 --- a/src/scatter_utilities/CMakeLists.txt +++ b/src/scatter_utilities/CMakeLists.txt @@ -18,6 +18,7 @@ set(${dir_EXE_SOURCES} create_tail_mask_from_ACFs.cxx upsample_and_fit_single_scatter.cxx simulate_scatter.cxx + scatter_likelihood_and_gradient.cxx ) include(stir_exe_targets) diff --git a/src/scatter_utilities/scatter_likelihood_and_gradient.cxx b/src/scatter_utilities/scatter_likelihood_and_gradient.cxx new file mode 100644 index 0000000000..fb375a1dbe --- /dev/null +++ b/src/scatter_utilities/scatter_likelihood_and_gradient.cxx @@ -0,0 +1,103 @@ +// +// +/* + Copyright (C) 2016, UCL + This file is part of STIR. + + This file is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + This file is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup utilities + \ingroup scatter + \brief Simulates a coarse scatter sinogram + + \author Nikos Efthimiou + + \par Usage: + \code + simulate_scatter parfile + \endcode + See stir::ScatterSimulation documentation for the format + of the parameter file. +*/ + +#include "stir/scatter/ScatterSimulation.h" +#include "stir/scatter/SingleScatterLikelihoodAndGradient.h".h" +#include "stir/Succeeded.h" +#include "stir/CPUTimer.h" +#include "stir/HighResWallClockTimer.h" +#include "stir/IO/read_from_file.h" + +using std::cerr; +using std::cout; +using std::endl; + +static void print_usage_and_exit() +{ + std::cerr<<"This executable runs a Scatter simulation method based on the options " + "in a parameter file"; + std::cerr<<"\nUsage:\n simulate_scatter scatter_simulation.par\n"; + std::cerr<<"Example parameter file:\n\n" + "Scatter Simulation :=\n" + "Simulation method := Single Scatter Simulation\n" + "Scatter Simulation Parameters :=\n" + " template projdata filename :=\n" + "attenuation image filename := \n" + "attenuation image for scatter points filename := \n" + "activity image filename :=\n" + "output filename prefix := ${OUTPUT_PREFIX}\n" + " scatter level := 1\n" + "attenuation threshold := 0.01\n" + "random := 1\n" + "use cache := 1\n" + "End Scatter Simulation Parameters :=\n" + "End Scatter Simulation:="<< std::endl; + exit(EXIT_FAILURE); +} +/***********************************************************/ + +int main(int argc, const char *argv[]) +{ + USING_NAMESPACE_STIR + + HighResWallClockTimer t; + t.reset(); + t.start(); + + if (argc!=2) + print_usage_and_exit(); + + shared_ptr < SingleScatterLikelihoodAndGradient > + simulation_method_sptr; + + KeyParser parser; + parser.add_start_key("Scatter Simulation"); + parser.add_stop_key("End Scatter Simulation"); + parser.add_parsing_key("Simulation method", &simulation_method_sptr); + parser.parse(argv[1]); + + shared_ptr data = ProjData::read_from_file("simulated_scatter_sino.hs"); + //const float rescale = 1; + shared_ptr > a(read_from_file >("true_atn_image.hv")); + VoxelsOnCartesianGrid& gradient_image = dynamic_cast< VoxelsOnCartesianGrid& > (*a); + gradient_image.fill(0); + + + double g= simulation_method_sptr->L_G_function(*data,gradient_image, false); + + t.stop(); + cout << "Total Wall clock timetime: " << t.value() << " seconds" << endl; + return EXIT_SUCCESS; +} + diff --git a/src/scatter_utilities/simulate_scatter.cxx b/src/scatter_utilities/simulate_scatter.cxx index 781a257b01..98d0d59551 100644 --- a/src/scatter_utilities/simulate_scatter.cxx +++ b/src/scatter_utilities/simulate_scatter.cxx @@ -62,10 +62,8 @@ int main(int argc, const char *argv[]) if (argc!=2) print_usage_and_exit(); - shared_ptr < ScatterSimulation > simulation_method_sptr; - KeyParser parser; parser.add_start_key("Scatter Simulation Parameters"); parser.add_stop_key("End Scatter Simulation Parameters"); diff --git a/src/swig/CMakeLists.txt b/src/swig/CMakeLists.txt index 1ede5d7ccb..cdf64042df 100644 --- a/src/swig/CMakeLists.txt +++ b/src/swig/CMakeLists.txt @@ -218,7 +218,7 @@ if (BUILD_SWIG_MATLAB) set(MATLAB_DEST ${CMAKE_INSTALL_PREFIX}/matlab CACHE PATH "Destination for Matlab module (relative to CMAKE_INSTALL_PREFIX)") INSTALL(TARGETS ${SWIG_MODULE_stirMATLAB_REAL_NAME} DESTINATION ${MATLAB_DEST}) INSTALL(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/+${module_name} DESTINATION ${MATLAB_DEST}) - file(GLOB SwigMatlabFiles "${CMAKE_CURRENT_BINARY_DIR}/Swig*.m") + file(GLOB SwigMatlabFiles "${CMAKE_CURRENT_BINARY_DIR}/*.m") INSTALL(FILES ${SwigMatlabFiles} DESTINATION ${MATLAB_DEST}) -endif (BUILD_SWIG_MATLAB) +endif (BUILD_SWIG_MATLAB) \ No newline at end of file diff --git a/src/swig/Makefile b/src/swig/Makefile index 2b54b1bdd6..2797d94206 100644 --- a/src/swig/Makefile +++ b/src/swig/Makefile @@ -60,4 +60,4 @@ _stir.so: stir_PYTHONwrap.o $(LIBS) Makefile stir.oct: stir_OCTAVEwrap.cpp $(LIBS) Makefile - CXXFLAGS=-g mkoctfile -v -I../include/ -g -DSWIG -o $@ $< $(LIBSa) + CXXFLAGS=-g mkoctfile -v -I../include/ -g -DSWIG -o $@ $< $(LIBSa) \ No newline at end of file diff --git a/src/swig/pyfragments.swg b/src/swig/pyfragments.swg index 33e78385fa..e933497bb6 100644 --- a/src/swig/pyfragments.swg +++ b/src/swig/pyfragments.swg @@ -16,53 +16,53 @@ */ %fragment(SWIG_AsVal_frag(long), "header", - fragment="SWIG_CanCastAsInteger", + fragment="SWIG_CanCastAsInteger", fragment="NumPy_Backward_Compatibility") { - SWIGINTERN int - SWIG_AsVal_dec(long)(PyObject * obj, long * val) - { - if (PyInt_Check(obj)) { - if (val) *val = PyInt_AsLong(obj); - return SWIG_OK; - } else if (PyLong_Check(obj)) { - long v = PyLong_AsLong(obj); - if (!PyErr_Occurred()) { - if (val) *val = v; - return SWIG_OK; - } else { - PyErr_Clear(); - } - } -%#ifdef SWIG_PYTHON_CAST_MODE - { - int dispatch = 0; - long v = PyInt_AsLong(obj); - if (!PyErr_Occurred()) { - if (val) *val = v; - return SWIG_AddCast(SWIG_OK); - } else { - PyErr_Clear(); - } - if (!dispatch) { - double d; - int res = SWIG_AddCast(SWIG_AsVal(double)(obj,&d)); - if (SWIG_IsOK(res) && SWIG_CanCastAsInteger(&d, LONG_MIN, LONG_MAX)) { - if (val) *val = (long)(d); - return res; - } - } - } -%#endif - if (PyArray_IsScalar(obj,Integer)) + SWIGINTERN int + SWIG_AsVal_dec(long)(PyObject * obj, long * val) { - PyArray_Descr * descr = PyArray_DescrNewFromType(NPY_LONG); - PyArray_CastScalarToCtype(obj, (void*)val, descr); - Py_DECREF(descr); - return SWIG_AddCast(SWIG_OK); + if (PyInt_Check(obj)) { + if (val) *val = PyInt_AsLong(obj); + return SWIG_OK; + } else if (PyLong_Check(obj)) { + long v = PyLong_AsLong(obj); + if (!PyErr_Occurred()) { + if (val) *val = v; + return SWIG_OK; + } else { + PyErr_Clear(); + } + } + %#ifdef SWIG_PYTHON_CAST_MODE + { + int dispatch = 0; + long v = PyInt_AsLong(obj); + if (!PyErr_Occurred()) { + if (val) *val = v; + return SWIG_AddCast(SWIG_OK); + } else { + PyErr_Clear(); + } + if (!dispatch) { + double d; + int res = SWIG_AddCast(SWIG_AsVal(double)(obj,&d)); + if (SWIG_IsOK(res) && SWIG_CanCastAsInteger(&d, LONG_MIN, LONG_MAX)) { + if (val) *val = (long)(d); + return res; + } + } + } + %#endif + if (PyArray_IsScalar(obj,Integer)) + { + PyArray_Descr * descr = PyArray_DescrNewFromType(NPY_LONG); + PyArray_CastScalarToCtype(obj, (void*)val, descr); + Py_DECREF(descr); + return SWIG_AddCast(SWIG_OK); + } + return SWIG_TypeError; } - return SWIG_TypeError; - } } @@ -72,119 +72,119 @@ */ %fragment(SWIG_AsVal_frag(unsigned long),"header", - fragment="SWIG_CanCastAsInteger", + fragment="SWIG_CanCastAsInteger", fragment="NumPy_Backward_Compatibility") { - SWIGINTERN int - SWIG_AsVal_dec(unsigned long)(PyObject *obj, unsigned long *val) - { - if (PyInt_Check(obj)) { - long v = PyInt_AsLong(obj); - if (v >= 0) { - if (val) *val = v; - return SWIG_OK; - } else { - return SWIG_OverflowError; - } - } else if (PyLong_Check(obj)) { - unsigned long v = PyLong_AsUnsignedLong(obj); - if (!PyErr_Occurred()) { - if (val) *val = v; - return SWIG_OK; - } else { - PyErr_Clear(); - } - } -%#ifdef SWIG_PYTHON_CAST_MODE + SWIGINTERN int + SWIG_AsVal_dec(unsigned long)(PyObject *obj, unsigned long *val) { - int dispatch = 0; - unsigned long v = PyLong_AsUnsignedLong(obj); - if (!PyErr_Occurred()) { - if (val) *val = v; - return SWIG_AddCast(SWIG_OK); - } else { - PyErr_Clear(); - } - if (!dispatch) { - double d; - int res = SWIG_AddCast(SWIG_AsVal(double)(obj,&d)); - if (SWIG_IsOK(res) && SWIG_CanCastAsInteger(&d, 0, ULONG_MAX)) { - if (val) *val = (unsigned long)(d); - return res; - } - } + if (PyInt_Check(obj)) { + long v = PyInt_AsLong(obj); + if (v >= 0) { + if (val) *val = v; + return SWIG_OK; + } else { + return SWIG_OverflowError; + } + } else if (PyLong_Check(obj)) { + unsigned long v = PyLong_AsUnsignedLong(obj); + if (!PyErr_Occurred()) { + if (val) *val = v; + return SWIG_OK; + } else { + PyErr_Clear(); + } + } + %#ifdef SWIG_PYTHON_CAST_MODE + { + int dispatch = 0; + unsigned long v = PyLong_AsUnsignedLong(obj); + if (!PyErr_Occurred()) { + if (val) *val = v; + return SWIG_AddCast(SWIG_OK); + } else { + PyErr_Clear(); + } + if (!dispatch) { + double d; + int res = SWIG_AddCast(SWIG_AsVal(double)(obj,&d)); + if (SWIG_IsOK(res) && SWIG_CanCastAsInteger(&d, 0, ULONG_MAX)) { + if (val) *val = (unsigned long)(d); + return res; + } + } + } + %#endif + if (PyArray_IsScalar(obj,Integer)) + { + PyArray_Descr * descr = PyArray_DescrNewFromType(NPY_ULONG); + PyArray_CastScalarToCtype(obj, (void*)val, descr); + Py_DECREF(descr); + return SWIG_AddCast(SWIG_OK); + } + return SWIG_TypeError; } -%#endif - if (PyArray_IsScalar(obj,Integer)) - { - PyArray_Descr * descr = PyArray_DescrNewFromType(NPY_ULONG); - PyArray_CastScalarToCtype(obj, (void*)val, descr); - Py_DECREF(descr); - return SWIG_AddCast(SWIG_OK); - } - return SWIG_TypeError; - } } %fragment(SWIG_AsVal_frag(double),"header", - fragment="NumPy_Backward_Compatibility") -{ -SWIGINTERN int -SWIG_AsVal_dec(double)(PyObject *obj, double *val) + fragment="NumPy_Backward_Compatibility") { - int res = SWIG_TypeError; - if (PyFloat_Check(obj)) { - if (val) *val = PyFloat_AsDouble(obj); - return SWIG_OK; - } else if (PyInt_Check(obj)) { - if (val) *val = PyInt_AsLong(obj); - return SWIG_OK; - } else if (PyLong_Check(obj)) { - double v = PyLong_AsDouble(obj); - if (!PyErr_Occurred()) { - if (val) *val = v; - return SWIG_OK; - } else { - PyErr_Clear(); - } - } -%#ifdef SWIG_PYTHON_CAST_MODE - { - int dispatch = 0; - double d = PyFloat_AsDouble(obj); - if (!PyErr_Occurred()) { - if (val) *val = d; - return SWIG_AddCast(SWIG_OK); - } else { - PyErr_Clear(); - } - if (!dispatch) { - long v = PyLong_AsLong(obj); - if (!PyErr_Occurred()) { - if (val) *val = v; - return SWIG_AddCast(SWIG_AddCast(SWIG_OK)); - } else { - PyErr_Clear(); - } - } - } -%#endif - if (PyArray_IsScalar(obj,Double)) + SWIGINTERN int + SWIG_AsVal_dec(double)(PyObject *obj, double *val) { - PyArray_Descr * descr = PyArray_DescrNewFromType(NPY_DOUBLE); - PyArray_CastScalarToCtype(obj, (void*)val, descr); - Py_DECREF(descr); - return SWIG_AddCast(SWIG_OK); + int res = SWIG_TypeError; + if (PyFloat_Check(obj)) { + if (val) *val = PyFloat_AsDouble(obj); + return SWIG_OK; + } else if (PyInt_Check(obj)) { + if (val) *val = PyInt_AsLong(obj); + return SWIG_OK; + } else if (PyLong_Check(obj)) { + double v = PyLong_AsDouble(obj); + if (!PyErr_Occurred()) { + if (val) *val = v; + return SWIG_OK; + } else { + PyErr_Clear(); + } + } + %#ifdef SWIG_PYTHON_CAST_MODE + { + int dispatch = 0; + double d = PyFloat_AsDouble(obj); + if (!PyErr_Occurred()) { + if (val) *val = d; + return SWIG_AddCast(SWIG_OK); + } else { + PyErr_Clear(); + } + if (!dispatch) { + long v = PyLong_AsLong(obj); + if (!PyErr_Occurred()) { + if (val) *val = v; + return SWIG_AddCast(SWIG_AddCast(SWIG_OK)); + } else { + PyErr_Clear(); + } + } + } + %#endif + if (PyArray_IsScalar(obj,Double)) + { + PyArray_Descr * descr = PyArray_DescrNewFromType(NPY_DOUBLE); + PyArray_CastScalarToCtype(obj, (void*)val, descr); + Py_DECREF(descr); + return SWIG_AddCast(SWIG_OK); + } + if (PyArray_IsScalar(obj,Float)) + { + float fval; + PyArray_Descr * descr = PyArray_DescrNewFromType(NPY_FLOAT32); + PyArray_CastScalarToCtype(obj, (void*)&fval, descr); + *val=(double)fval; + Py_DECREF(descr); + return SWIG_AddCast(SWIG_OK); + } + return res; } - if (PyArray_IsScalar(obj,Float)) - { - float fval; - PyArray_Descr * descr = PyArray_DescrNewFromType(NPY_FLOAT32); - PyArray_CastScalarToCtype(obj, (void*)&fval, descr); - *val=(double)fval; - Py_DECREF(descr); - return SWIG_AddCast(SWIG_OK); - } - return res; -} } diff --git a/src/swig/stir.i b/src/swig/stir.i index f5a3a73541..38bbdccf52 100644 --- a/src/swig/stir.i +++ b/src/swig/stir.i @@ -9,13 +9,13 @@ See STIR/LICENSE.txt for details */ /*! - \file + \file \brief Interface file for SWIG \author Kris Thielemans \author Markus Jehl */ - +// %module stir %{ @@ -78,7 +78,7 @@ #include "stir/zoom.h" #include "stir/GeneralisedPoissonNoiseGenerator.h" - + #include "stir/IO/read_from_file.h" #include "stir/IO/write_to_file.h" #include "stir/IO/InterfileOutputFileFormat.h" @@ -104,7 +104,7 @@ #include "stir/HUToMuImageProcessor.h" #endif -#include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h" +#include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h" #include "stir/OSMAPOSL/OSMAPOSLReconstruction.h" #include "stir/OSSPS/OSSPSReconstruction.h" #include "stir/recon_buildblock/ForwardProjectorByBinUsingProjMatrixByBin.h" @@ -142,6 +142,13 @@ #include #include +#include "stir/IO/write_to_file.h" + +#include "stir/scatter/ScatterSimulation.h" +#include "stir/scatter/SingleScatterSimulation.h" +#include "stir/scatter/SingleScatterLikelihoodAndGradient.h" +#include "stir/scatter/ScatterEstimation.h" + // TODO need this (bug in swig) // this bug occurs (only?) when using "%template(name) someclass;" inside the namespace // as opposed to "%template(name) stir::someclass" outside the namespace @@ -155,7 +162,7 @@ SWIG_AsVal_double (PyObject * obj, double *val); #endif - // local helper functions for conversions etc. These are not "exposed" to the target language + // local helper functions for conversions etc. These are not "exposed" to the target language // (but only enter in the wrapper) namespace swigstir { #if defined(SWIGPYTHON) @@ -216,15 +223,15 @@ { PyObject *iterator = PyObject_GetIter(arg); - + PyObject *item; typename stir::Array::full_iterator array_iter = array_ptr->begin_all(); - while ((item = PyIter_Next(iterator)) && array_iter != array_ptr->end_all()) + while ((item = PyIter_Next(iterator)) && array_iter != array_ptr->end_all()) { double val; // TODO currently hard-wired as double which might imply extra conversions int ecode = SWIG_AsVal_double(item, &val); - if (SWIG_IsOK(ecode)) + if (SWIG_IsOK(ecode)) { *array_iter++ = static_cast(val); } @@ -246,7 +253,7 @@ } Py_DECREF(iterator); - if (PyErr_Occurred()) + if (PyErr_Occurred()) { throw std::runtime_error("Error during fill()"); } @@ -255,7 +262,7 @@ } #if 0 - + // TODO does not work yet. // it doesn't compile as includes are in init section, which follows after this in the wrapper // Even if it did compile, it might not work anyway as I haven't tested it. @@ -364,7 +371,7 @@ float const* data_ptr = (float *)mxGetData(pm); a.fill(static_cast(*data_ptr)); } else - { + { throw std::runtime_error("currently only supporting double or single arrays for filling a stir array"); } } @@ -387,8 +394,8 @@ } if (matlab_num_dims > static_cast(num_dimensions)) { - throw std::runtime_error(boost::str(boost::format("number of dimensions in matlab array is incorrect for constructing a stir array of dimension %d") % - num_dimensions)); + throw std::runtime_error(boost::str(boost::format("number of dimensions in matlab array is incorrect for constructing a stir array of dimension %d") % + num_dimensions)); } if (do_resize) { @@ -399,7 +406,7 @@ a.resize(sizes); } else - { + { // check sizes BasicCoordinate minind,maxind; a.get_regular_range(minind,maxind); @@ -415,7 +422,7 @@ throw std::runtime_error("sizes of first dimensions of the stir array have to be 1 if initialising from a lower dimensional matlab array"); } } - } + } if (mxIsDouble(pm)) { double * data_ptr = mxGetPr(pm); @@ -425,10 +432,10 @@ float * data_ptr = (float *)mxGetData(pm); std::copy(data_ptr, data_ptr+a.size_all(), a.begin_all()); } else - { + { throw std::runtime_error("currently only supporting double or single arrays for constructing a stir array"); } - } + } //////////// same for Coordinate @@ -457,7 +464,7 @@ float const* data_ptr = (float *)mxGetData(pm); a.fill(static_cast(*data_ptr)); } else - { + { throw std::runtime_error("currently only supporting double or single arrays for filling a stir coordinate"); } } @@ -480,8 +487,8 @@ } if (matlab_num_dims != static_cast(1)) { - throw std::runtime_error(boost::str(boost::format("number of dimensions %d of matlab array is incorrect for constructing a stir coordinate of dimension %d (expecting a column vector)") % - matlab_num_dims % num_dimensions)); + throw std::runtime_error(boost::str(boost::format("number of dimensions %d of matlab array is incorrect for constructing a stir coordinate of dimension %d (expecting a column vector)") % + matlab_num_dims % num_dimensions)); } if (m_sizes[0]!=static_cast(num_dimensions)) { @@ -496,10 +503,10 @@ float * data_ptr = (float *)mxGetData(pm); std::copy(data_ptr, data_ptr+a.size(), a.begin()); } else - { + { throw std::runtime_error("currently only supporting double or single arrays for constructing a stir array"); } - } + } #endif } // end namespace swigstir %} @@ -545,7 +552,7 @@ } catch (const std::string& e) { SWIG_exception(SWIG_RuntimeError, e.c_str()); } -} +} // declare some functions that return a new pointer such that SWIG can release memory properly %newobject *::clone; @@ -557,7 +564,7 @@ %ignore *::read_from_stream; #if defined(SWIGPYTHON) -%rename(__assign__) *::operator=; +%rename(__assign__) *::operator=; #endif // include standard swig support for some bits of the STL (i.e. standard C++ lib) @@ -609,7 +616,7 @@ namespace std { %template(StringList) list; } -// section for helper classes for creating new iterators. +// section for helper classes for creating new iterators. // The code here is nearly a copy of what's in PyIterators.swg, // except that the decr() function isn't defined. This is because we need it for some STIR iterators // which are forward_iterators. @@ -620,7 +627,7 @@ namespace std { %{ namespace swigstir { #ifdef SWIGPYTHON - template::value_type, typename FromOper = swig::from_oper > class SwigPyForwardIteratorClosed_T : public swig::SwigPyIterator_T @@ -629,14 +636,14 @@ namespace std { FromOper from; typedef OutIterator out_iterator; typedef ValueType value_type; - typedef swig::SwigPyIterator_T base; + typedef swig::SwigPyIterator_T base; typedef SwigPyForwardIteratorClosed_T self_type; - + SwigPyForwardIteratorClosed_T(out_iterator curr, out_iterator first, out_iterator last, PyObject *seq) : swig::SwigPyIterator_T(curr, seq), begin(first), end(last) { } - + PyObject *value() const { if (base::current == end) { throw swig::stop_iteration(); @@ -644,7 +651,7 @@ namespace std { return swig::from(static_cast(*(base::current))); } } - + swig::SwigPyIterator *copy() const { return new self_type(*this); @@ -697,13 +704,13 @@ namespace std { copy_to(proj_data, array_iter); return array; } - + } // end of namespace %} // end of initial code specification for inclusino in the SWIG wrapper // doesn't work (yet?) because of bug in int template arguments -// %rename(__getitem__) *::at; +// %rename(__getitem__) *::at; // MACROS to define index access (Work-in-progress) %define %ADD_indexaccess(INDEXTYPE,RETTYPE,TYPE...) @@ -883,7 +890,7 @@ T * operator-> () const; // use simpler version for SWIG to make the hierarchy a bit easier namespace stir { template - class RegisteredObject + class RegisteredObject { public: //! List all possible registered names to the stream @@ -1010,14 +1017,14 @@ T * operator-> () const; static stir::VoxelsOnCartesianGrid * read_from_file(const std::string& filename) { using namespace stir; - unique_ptr > + unique_ptr > ret(read_from_file >(filename)); return dynamic_cast *>(ret.release()); } } //%ADD_indexaccess(int,stir::BasicCoordinate::value_type,stir::BasicCoordinate); -namespace stir { +namespace stir { #ifdef SWIGPYTHON // add extra features to the coordinates to make them a bit more Python friendly %extend BasicCoordinate { @@ -1034,7 +1041,7 @@ namespace stir { // print as (1,2,3) as opposed to non-informative default provided by SWIG std::string __str__() - { + { std::ostringstream s; s<<'('; for (int d=1; d<=num_dimensions-1; ++d) @@ -1045,7 +1052,7 @@ namespace stir { // print as classname((1,2,3)) as opposed to non-informative default provided by SWIG std::string __repr__() - { + { #if SWIG_VERSION < 0x020009 // don't know how to get the Python typename std::string repr = "stir.Coordinate"; @@ -1076,8 +1083,8 @@ namespace stir { return $self->size(); } #if defined(SWIGPYTHON_BUILTIN) - %feature("python:slot", "tp_str", functype="reprfunc") __str__; - %feature("python:slot", "tp_repr", functype="reprfunc") __repr__; + %feature("python:slot", "tp_str", functype="reprfunc") __str__; + %feature("python:slot", "tp_repr", functype="reprfunc") __repr__; %feature("python:slot", "nb_nonzero", functype="inquiry") __nonzero__; %feature("python:slot", "sq_length", functype="lenfunc") __len__; #endif // SWIGPYTHON_BUILTIN @@ -1087,14 +1094,14 @@ namespace stir { %extend BasicCoordinate { // print as [1;2;3] as opposed to non-informative default provided by SWIG void disp() - { + { std::ostringstream s; s<<'['; for (int d=1; d<=num_dimensions-1; ++d) s << (*$self)[d] << "; "; s << (*$self)[num_dimensions] << "]\n"; mexPrintf(s.str().c_str()); - + } //%feature("autodoc", "construct from vector, e.g. [2;3;4] for a 3d coordinate") BasicCoordinate(const mxArray *pm) @@ -1120,7 +1127,7 @@ namespace stir { %template(Float3Coordinate) Coordinate3D< float >; %template(FloatCartesianCoordinate3D) CartesianCoordinate3D; %template(IntCartesianCoordinate3D) CartesianCoordinate3D; - + %template(Int2BasicCoordinate) BasicCoordinate<2,int>; %template(Size2BasicCoordinate) BasicCoordinate<2,std::size_t>; %template(Float2BasicCoordinate) BasicCoordinate<2,float>; @@ -1206,7 +1213,7 @@ namespace stir { snprintf(str, 1000, "Wrong argument-type used for fill(): should be a scalar or an iterator or so, but is of type %s", arg->ob_type->tp_name); throw std::invalid_argument(str); - } + } } } #endif @@ -1252,7 +1259,7 @@ namespace stir { { swigstir::fill_Array_from_matlab(*$self, pm, false /*do not resize */); } } #endif - // TODO next line doesn't give anything useful as SWIG doesn't recognise that + // TODO next line doesn't give anything useful as SWIG doesn't recognise that // the return value is an array. So, we get a wrapped object that we cannot handle //%ADD_indexaccess(int,Array::value_type, Array); @@ -1303,7 +1310,7 @@ namespace stir { %rename (get_back_projector) *::get_back_projector_sptr; %rename (get_kappa) *::get_kappa_sptr; %rename (get_attenuation_image) *::get_attenuation_image_sptr; -/* would be nice, but needs swig to be compiled with PCRE support +/* would be nice, but needs swig to be compiled with PCRE support %rename("rstrip:[_ptr]") %rename("rstrip:[_sptr]") */ @@ -1331,6 +1338,11 @@ namespace stir { %include "stir/IO/write_to_file.h" %template(write_image_to_file) stir::write_to_file >; +%include "stir/IO/write_to_file.h" +%template(write_image_to_file) stir::write_to_file >; + +%include "stir/zoom.h" +%include "stir/ZoomOptions.h" #ifdef STIRSWIG_SHARED_PTR #define DataT stir::DiscretisedDensity<3,float> %shared_ptr(stir::OutputFileFormat >); @@ -1361,6 +1373,8 @@ namespace stir { /* Now do ProjDataInfo, Sinogram et al */ %include "stir/TimeFrameDefinitions.h" +%include "stir/zoom.h" +%include "stir/ZoomOptions.h" %include "stir/ExamInfo.h" %include "stir/ExamData.h" @@ -1458,7 +1472,7 @@ stir::CartesianCoordinate3D const int num_views, const int num_tangential_poss, const bool arc_corrected = true) { - return + return construct_proj_data_info(scanner_sptr, span, max_delta, num_views, num_tangential_poss, arc_corrected).get(); @@ -1485,7 +1499,7 @@ namespace stir { %feature("autodoc", "create a stir 3D Array from the projection data (internal)") to_array; %newobject to_array; Array<3,float> to_array() - { + { Array<3,float> array = swigstir::projdata_to_3D(*$self); return array; } @@ -1506,19 +1520,19 @@ namespace stir { snprintf(str, 1000, "Wrong argument-type used for fill(): should be a scalar or an iterator or so, but is of type %s", arg->ob_type->tp_name); throw std::invalid_argument(str); - } + } } #elif defined(SWIGMATLAB) %newobject to_matlab; mxArray * to_matlab() - { + { Array<3,float> array = swigstir::projdata_to_3D(*$self); - return swigstir::Array_to_matlab(array); + return swigstir::Array_to_matlab(array); } void fill(const mxArray *pm) - { + { Array<3,float> array; swigstir::fill_Array_from_matlab(array, pm, true); fill_from(*$self, array.begin_all(), array.end_all()); @@ -1565,7 +1579,7 @@ namespace stir { %include "stir/ProjDataInterfile.h" %include "stir/ProjDataInMemory.h" -namespace stir { +namespace stir { %template(FloatViewgram) Viewgram; %template(FloatSinogram) Sinogram; // TODO don't want to give a name @@ -1636,6 +1650,12 @@ namespace stir { #undef elemT #endif +/*#define elemT float +extend stir::DataProcessor { + stir::DataProcessor >* ptr() const { + return this; }} +#undef elemT*/ + %include "stir/DataProcessor.h" %include "stir/ChainedDataProcessor.h" %include "stir/SeparableCartesianMetzImageFilter.h" @@ -1646,6 +1666,8 @@ namespace stir { #endif #define elemT float + + %template(DataProcessor3DFloat) stir::DataProcessor >; %template(RPChainedDataProcessor3DFloat) stir::RegisteredParsingObject< stir::ChainedDataProcessor >, @@ -1903,7 +1925,7 @@ stir::RegisteredParsingObject< stir::LogcoshPrior, %include "stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h" -%template (internalRPForwardProjectorByBinUsingProjMatrixByBin) +%template (internalRPForwardProjectorByBinUsingProjMatrixByBin) stir::RegisteredParsingObject; %include "stir/recon_buildblock/ForwardProjectorByBinUsingProjMatrixByBin.h" @@ -1969,6 +1991,15 @@ void multiply_crystal_factors(stir::ProjData& proj_data, const stir::Array<2,flo stir::ScatterSimulation, stir::ScatterSimulation>; %include "stir/scatter/SingleScatterSimulation.h" +%shared_ptr(stir::RegisteredParsingObject< + stir::SingleScatterLikelihoodAndGradient, + stir::SingleScatterSimulation, + stir::SingleScatterSimulation + >); +%shared_ptr(stir::SingleScatterLikelihoodAndGradient); + +%include "stir/scatter/SingleScatterLikelihoodAndGradient.h" + %shared_ptr(stir::ScatterEstimation); %include "stir/scatter/ScatterEstimation.h" @@ -2009,3 +2040,4 @@ void multiply_crystal_factors(stir::ProjData& proj_data, const stir::Array<2,flo stir::BackProjectorByBin>; %include "stir/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.h" #endif + diff --git a/src/swig/stirextra.py b/src/swig/stirextra.py index ce475ec018..b7faa641c3 100644 --- a/src/swig/stirextra.py +++ b/src/swig/stirextra.py @@ -16,8 +16,8 @@ def get_bounding_box_indices(image): """ - return [min,max] tuple with indices of first and last voxel in a STIR image - """ + return [min,max] tuple with indices of first and last voxel in a STIR image + """ num_dims = image.get_num_dimensions(); if num_dims == 2: minind=stir.Int2BasicCoordinate() @@ -27,14 +27,13 @@ def get_bounding_box_indices(image): maxind=stir.Int3BasicCoordinate() else: raise exceptions.NotImplementedError('need to handle dimensions different from 2 and 3') - image.get_regular_range(minind, maxind); return [minind, maxind] def get_physical_coordinates_for_bounding_box(image): """ - return physical coordinates (in mm) for the first and last voxel in a STIR image - """ + return physical coordinates (in mm) for the first and last voxel in a STIR image + """ [minind,maxind]=get_bounding_box_indices(image); minphys=image.get_physical_coordinates_for_indices(minind); maxphys=image.get_physical_coordinates_for_indices(maxind); @@ -42,8 +41,8 @@ def get_physical_coordinates_for_bounding_box(image): def get_physical_coordinates_for_grid(image): """ - return [z,y,x] tuple of grid coordinates for a STIR image - """ + return [z,y,x] tuple of grid coordinates for a STIR image + """ [minind,maxind]=get_bounding_box_indices(image); sizes=maxind-minind+1; minphys=image.get_physical_coordinates_for_indices(minind); @@ -60,8 +59,8 @@ def get_physical_coordinates_for_grid(image): def to_numpy(stirdata): """ - return the data in a STIR image or other Array as a numpy array - """ + return the data in a STIR image or other Array as a numpy array + """ # construct a numpy array using the "flat" STIR iterator try: npstirdata=numpy.fromiter(stirdata.flat(), dtype=numpy.float32); diff --git a/src/test/CMakeLists.txt b/src/test/CMakeLists.txt index c92350c3e7..adec80d33b 100644 --- a/src/test/CMakeLists.txt +++ b/src/test/CMakeLists.txt @@ -37,6 +37,9 @@ Set(${dir_INVOLVED_TEST_EXE_SOURCES} # the next 2 are interactive, so we don't add a test for it, but only compile them test_display.cxx test_interpolate.cxx + test_upsample.cxx + test_SSRB.cxx + test_ScatterGradient.cxx ) set(buildblock_simple_tests @@ -65,6 +68,7 @@ set(buildblock_simple_tests test_export_array.cxx test_GeneralisedPoissonNoiseGenerator.cxx test_multiple_proj_data.cxx + test_zoom_image_adjoint.cxx ) include(stir_test_exe_targets) diff --git a/src/test/test_SSRB.cxx b/src/test/test_SSRB.cxx new file mode 100644 index 0000000000..86c869193a --- /dev/null +++ b/src/test/test_SSRB.cxx @@ -0,0 +1,203 @@ +// +// +/*! + + \file + \ingroup test + + \author Ludovica Brusaferri + +*/ +/* + Copyright (C) 2015, University College London + This file is part of STIR. + + This file is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + This file is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + See STIR/LICENSE.txt for details +*/ + +#include "stir/ProjDataInMemory.h" +#include "stir/inverse_SSRB.h" +#include "stir/ExamInfo.h" +#include "stir/ProjDataInfo.h" +#include "stir/Sinogram.h" +#include "stir/Viewgram.h" +#include "stir/Succeeded.h" +#include "stir/RunTests.h" +#include "stir/Scanner.h" +#include "stir/ProjDataInfoCylindricalNoArcCorr.h" +#include "stir/scatter/SingleScatterSimulation.h" +#include "stir/scatter/ScatterEstimation.h" +#include "stir/Sinogram.h" +#include "stir/Array.h" + +#include "stir/IO/write_to_file.h" + +#include + + +START_NAMESPACE_STIR + + +class SSRBTests: public RunTests +{ +public: + void run_tests(); + void fill_projdata_with_random(ProjData & projdata); +}; + +void +SSRBTests:: +fill_projdata_with_random(ProjData & projdata) +{ + + projdata.fill(0); + std::random_device random_device; + std::mt19937 random_number_generator(random_device()); + std::uniform_real_distribution number_distribution(0,1); + Bin bin; + { + for (bin.segment_num()=projdata.get_min_segment_num(); + bin.segment_num()<=projdata.get_max_segment_num(); + ++bin.segment_num()) + for (bin.axial_pos_num()= + projdata.get_min_axial_pos_num(bin.segment_num()); + bin.axial_pos_num()<=projdata.get_max_axial_pos_num(bin.segment_num()); + ++bin.axial_pos_num()) + { + + if((bin.segment_num()+bin.axial_pos_num())>projdata.get_max_axial_pos_num(bin.segment_num())||(bin.segment_num()+bin.axial_pos_num()) sino = projdata.get_empty_sinogram(bin.axial_pos_num(),bin.segment_num()); + + for (bin.view_num()=sino.get_min_view_num(); + bin.view_num()<=sino.get_max_view_num(); + ++bin.view_num()) + { + for (bin.tangential_pos_num()= + sino.get_min_tangential_pos_num(); + bin.tangential_pos_num()<= + sino.get_max_tangential_pos_num(); + ++bin.tangential_pos_num()) + sino[bin.view_num()][bin.tangential_pos_num()]= number_distribution(random_number_generator); + + projdata.set_sinogram(sino); + + } + + } + + } +} +void +SSRBTests:: +run_tests() +{ + + + std::cout << "========== TEST SSRB =========== \n"; + + shared_ptr scanner_sptr(new Scanner(Scanner::Siemens_mMR)); + + //creating proj data info + shared_ptr proj_data_info_sptr_4D(ProjDataInfo::ProjDataInfoCTI(scanner_sptr,/*span*/1, 2,/*views*/ 252, /*tang_pos*/344, /*arc_corrected*/ false)); + shared_ptr exam_info_sptr_4D(new ExamInfo); + ProjDataInMemory projdata_4D(exam_info_sptr_4D, proj_data_info_sptr_4D); + //shared_ptr proj_data_info_sptr_3D(ProjDataInfo::ProjDataInfoCTI(scanner_sptr,/*span*/1, 0,/*views*/ 252, /*tang_pos*/344, /*arc_corrected*/ false)); + + shared_ptr proj_data_info_sptr_3D(projdata_4D.get_proj_data_info_sptr()->clone()); + proj_data_info_sptr_3D->reduce_segment_range(0,0); //create input template + + ProjDataInMemory projdata_3D(projdata_4D.get_exam_info_sptr(),proj_data_info_sptr_3D); + + ProjDataInMemory A_4D(exam_info_sptr_4D, proj_data_info_sptr_4D); + ProjDataInMemory A_3D(projdata_3D.get_exam_info_sptr(),proj_data_info_sptr_3D); + // ProjDataInterfile A_4D(exam_info_sptr_4D, proj_data_info_sptr_4D, "4D.hs",std::ios::out|std::ios::in|std::ios::app); + //ProjDataInterfile A_3D(projdata_4D.get_exam_info_sptr(),proj_data_info_sptr_3D, "3D.hs",std::ios::out|std::ios::in|std::ios::app); + + fill_projdata_with_random(projdata_4D); + //projdata_4D.fill(1); + fill_projdata_with_random(projdata_3D); + //projdata_3D.fill(1); + A_4D.fill(0); //initialise output + A_3D.fill(0); //initialise output + + + std::cout << "-------- Testing Inverse SSRB --------\n"; + inverse_SSRB(A_4D,projdata_3D); + + std::cout << "-------- Testing Transpose Inverse SSRB --------\n"; + + transpose_inverse_SSRB(A_3D,projdata_4D); + + std::cout << "-------- = --------\n"; + + float cdot1 = 0; + float cdot2 = 0; + + for ( int segment_num = projdata_4D.get_min_segment_num(); segment_num <= projdata_4D.get_max_segment_num(); ++segment_num) + { + for (int axial_pos = projdata_4D.get_min_axial_pos_num(segment_num); axial_pos <= projdata_4D.get_max_axial_pos_num(segment_num); ++axial_pos) + { + + const Sinogram yy_sinogram = projdata_4D.get_sinogram(axial_pos, segment_num); + const Sinogram At_sinogram = A_4D.get_sinogram(axial_pos, segment_num); + + for ( int view_num = projdata_4D.get_min_view_num(); view_num <= projdata_4D.get_max_view_num();view_num++) + { + for ( int tang_pos = projdata_4D.get_min_tangential_pos_num(); tang_pos <= projdata_4D.get_max_tangential_pos_num(); tang_pos++) + { + cdot1 += At_sinogram[view_num][tang_pos]*yy_sinogram[view_num][tang_pos]; + } + } + } + } + + for ( int segment_num = projdata_3D.get_min_segment_num(); segment_num <= projdata_3D.get_max_segment_num(); ++segment_num) + { + for (int axial_pos = projdata_3D.get_min_axial_pos_num(segment_num); axial_pos <= projdata_3D.get_max_axial_pos_num(segment_num); ++axial_pos) + { + + const Sinogram xx_sinogram = projdata_3D.get_sinogram(axial_pos, segment_num); + const Sinogram A_sinogram = A_3D.get_sinogram(axial_pos, segment_num); + + for ( int view_num = projdata_3D.get_min_view_num(); view_num <= projdata_3D.get_max_view_num();view_num++) + { + for ( int tang_pos = projdata_3D.get_min_tangential_pos_num(); tang_pos <= projdata_3D.get_max_tangential_pos_num(); tang_pos++) + { + cdot2 += xx_sinogram[view_num][tang_pos]*A_sinogram[view_num][tang_pos]; + } + } + } + } + + std::cout << cdot1 << "=" << cdot2 << '\n'; + set_tolerance(0.008); + check_if_equal(cdot1, cdot2, "test adjoint"); + + + + +} + +END_NAMESPACE_STIR + + +USING_NAMESPACE_STIR + +int main() +{ + SSRBTests tests; + tests.run_tests(); + return tests.main_return_value(); +} diff --git a/src/test/test_ScatterGradient.cxx b/src/test/test_ScatterGradient.cxx new file mode 100644 index 0000000000..64d59929a4 --- /dev/null +++ b/src/test/test_ScatterGradient.cxx @@ -0,0 +1,95 @@ +// +// +/* + Copyright (C) 2019, University of Hull + This file is part of STIR. + + This file is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + This file is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + See STIR/LICENSE.txt for details +*/ +/*! + + \file + \ingroup test + + \brief Test program for ScatterGradient + + \author Ludovica Brusaferri + +*/ + +#include "stir/RunTests.h" +#include "stir/Scanner.h" +#include "stir/Viewgram.h" +#include "stir/Succeeded.h" +#include "stir/ProjDataInfoCylindricalNoArcCorr.h" +#include "stir/VoxelsOnCartesianGrid.h" +#include "stir/scatter/SingleScatterSimulation.h" +#include "stir/scatter/SingleScatterLikelihoodAndGradient.h" +#include "stir/recon_buildblock/ForwardProjectorByBin.h" +#include "stir/recon_buildblock/ForwardProjectorByBinUsingRayTracing.h" +#include "stir/recon_buildblock/ForwardProjectorByBinUsingProjMatrixByBin.h" +#include "stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h" +#include "stir/recon_buildblock/BinNormalisationFromAttenuationImage.h" +#include "stir/Shape/EllipsoidalCylinder.h" +#include "stir/Shape/Box3D.h" +#include "stir/IO/write_to_file.h" +#include +#include +#include "stir/centre_of_gravity.h" + +using std::cerr; +using std::endl; +using std::string; + +START_NAMESPACE_STIR + +/*! + \ingroup test + \brief Test class for Scanner +*/ +class ScatterGradientTests: public RunTests +{ +public: + void run_tests(); +private: + void test_scatter_gradient(); + +}; + + + +void +ScatterGradientTests::test_scatter_gradient() +{ + //TODO +} + +void +ScatterGradientTests:: +run_tests() +{ + + test_scatter_gradient(); +} + + +END_NAMESPACE_STIR + + +int main() +{ + USING_NAMESPACE_STIR + + ScatterGradientTests tests; + tests.run_tests(); + return tests.main_return_value(); +} diff --git a/src/test/test_upsample.cxx b/src/test/test_upsample.cxx new file mode 100644 index 0000000000..16f4cabb6d --- /dev/null +++ b/src/test/test_upsample.cxx @@ -0,0 +1,279 @@ +// +// +/*! + + \file + \ingroup test + + \author Ludovica Brusaferri + +*/ +/* + Copyright (C) 2015, University College London + This file is part of STIR. + + This file is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + This file is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + See STIR/LICENSE.txt for details +*/ + +#include "stir/ProjDataInMemory.h" +#include "stir/inverse_SSRB.h" +#include "stir/ExamInfo.h" +#include "stir/ProjDataInfo.h" +#include "stir/Sinogram.h" +#include "stir/Viewgram.h" +#include "stir/Succeeded.h" +#include "stir/RunTests.h" +#include "stir/Scanner.h" +#include "stir/ProjDataInfoCylindricalNoArcCorr.h" +#include "stir/scatter/SingleScatterSimulation.h" +#include "stir/scatter/ScatterEstimation.h" +#include "stir/Sinogram.h" +#include "stir/Array.h" + +#include "stir/IO/write_to_file.h" + +#include + + + +START_NAMESPACE_STIR + + +class UpsampleDownsampleTests: public RunTests +{ +public: + void run_tests(); + void fill_projdata_with_random(ProjData & projdata); +}; + +void +UpsampleDownsampleTests:: +fill_projdata_with_random(ProjData & projdata) +{ + + std::random_device random_device; + std::mt19937 random_number_generator(random_device()); + std::uniform_real_distribution number_distribution(0,10); + Bin bin; + { + for (bin.segment_num()=projdata.get_min_segment_num(); + bin.segment_num()<=projdata.get_max_segment_num(); + ++bin.segment_num()) + for (bin.axial_pos_num()= + projdata.get_min_axial_pos_num(bin.segment_num()); + bin.axial_pos_num()<=projdata.get_max_axial_pos_num(bin.segment_num()); + ++bin.axial_pos_num()) + { + + Sinogram sino = projdata.get_empty_sinogram(bin.axial_pos_num(),bin.segment_num()); + + for (bin.view_num()=sino.get_min_view_num(); + bin.view_num()<=sino.get_max_view_num(); + ++bin.view_num()) + { + for (bin.tangential_pos_num()= + sino.get_min_tangential_pos_num(); + bin.tangential_pos_num()<= + sino.get_max_tangential_pos_num(); + ++bin.tangential_pos_num()) + sino[bin.view_num()][bin.tangential_pos_num()]= number_distribution(random_number_generator); + + projdata.set_sinogram(sino); + + } + + } + + } +} +void +UpsampleDownsampleTests:: +run_tests() +{ +// std::cout << "-------- Testing Upsampling and Downsampling ---------\n"; + + shared_ptr scanner_sptr(new Scanner(Scanner::Siemens_mMR)); + shared_ptr exam_info_sptr(new ExamInfo); + shared_ptr LR_exam_info_sptr(new ExamInfo); + + //creating proj data info + shared_ptr proj_data_info_sptr(ProjDataInfo::ProjDataInfoCTI(scanner_sptr,/*span*/1, 0,/*views*/ 252, /*tang_pos*/344, /*arc_corrected*/ false)); + + // shared_ptr LR = ProjData::read_from_file("simulated_scatter_sino_UU2.hs"); + + + // construct y + ProjDataInMemory y(exam_info_sptr, proj_data_info_sptr); + + + unique_ptr sss(new SingleScatterSimulation()); + sss->set_template_proj_data_info(*y.get_proj_data_info_sptr()); + int down_rings = static_cast(scanner_sptr->get_num_rings()/6); + int down_dets = static_cast(scanner_sptr->get_max_num_views()/6); + + sss->downsample_scanner(down_rings, down_dets); + + + //shared_ptr sss_projdata(sss->get_template_proj_data_info_sptr()); + + + // construct x + ProjDataInMemory x(sss->get_exam_info_sptr(), sss->get_template_proj_data_info_sptr()); + + // construct Ax + + ProjDataInMemory Ax(exam_info_sptr, proj_data_info_sptr); + // construct x + ProjDataInMemory Aty(sss->get_exam_info_sptr(), sss->get_template_proj_data_info_sptr()); + + ProjDataInMemory N(exam_info_sptr, proj_data_info_sptr); + + fill_projdata_with_random(x); + fill_projdata_with_random(y); + + Ax.fill(0); //initialise output + Aty.fill(0); //initialise output + + //N.fill(10); + fill_projdata_with_random(N); + + bool remove_interleaving = false; + + std::cout << "========== REMOVE INTERLEAVING = FALSE =========== \n"; + + std::cout << "-------- Testing Pull --------\n"; + ScatterEstimation::pull_scatter_estimate(Ax,y,x,N,remove_interleaving); + + std::cout << "-------- Testing Push --------\n"; + + ScatterEstimation::push_scatter_estimate(Aty,x,y,N,remove_interleaving); + + std::cout << "-------- = --------\n"; + + float cdot1 = 0; + float cdot2 = 0; + + for ( int segment_num = y.get_min_segment_num(); segment_num <= y.get_max_segment_num(); ++segment_num) + { + for (int axial_pos = y.get_min_axial_pos_num(segment_num); axial_pos <= y.get_max_axial_pos_num(segment_num); ++axial_pos) + { + + const Sinogram y_sinogram = y.get_sinogram(axial_pos, segment_num); + const Sinogram Ax_sinogram = Ax.get_sinogram(axial_pos, segment_num); + + for ( int view_num = y.get_min_view_num(); view_num <= y.get_max_view_num();view_num++) + { + for ( int tang_pos = y.get_min_tangential_pos_num(); tang_pos <= y.get_max_tangential_pos_num(); tang_pos++) + { + cdot1 += Ax_sinogram[view_num][tang_pos]*y_sinogram[view_num][tang_pos]; + } + } + } + } + + for ( int segment_num = x.get_min_segment_num(); segment_num <= x.get_max_segment_num(); ++segment_num) + { + for (int axial_pos = x.get_min_axial_pos_num(segment_num); axial_pos <= x.get_max_axial_pos_num(segment_num); ++axial_pos) + { + + const Sinogram x_sinogram = x.get_sinogram(axial_pos, segment_num); + const Sinogram Aty_sinogram = Aty.get_sinogram(axial_pos, segment_num); + + for ( int view_num = x.get_min_view_num(); view_num <= x.get_max_view_num();view_num++) + { + for ( int tang_pos = x.get_min_tangential_pos_num(); tang_pos <= x.get_max_tangential_pos_num(); tang_pos++) + { + cdot2 += x_sinogram[view_num][tang_pos]*Aty_sinogram[view_num][tang_pos]; + } + } + } + } + + std::cout << cdot1 << "=" << cdot2 << '\n'; + set_tolerance(0.006); + check_if_equal(cdot1, cdot2, "test adjoint"); + + + std::cout << "========== REMOVE INTERLEAVING = TRUE =========== \n"; + + Ax.fill(0); //initialise output + Aty.fill(0); //initialise output + + remove_interleaving = true; + + std::cout << "-------- Testing Pull --------\n"; + ScatterEstimation::pull_scatter_estimate(Ax,y,x,N,remove_interleaving); + + std::cout << "-------- Testing Push --------\n"; + + ScatterEstimation::push_scatter_estimate(Aty,x,y,N,remove_interleaving); + + std::cout << "-------- = --------\n"; + + cdot1 = 0; + cdot2 = 0; + + for ( int segment_num = y.get_min_segment_num(); segment_num <= y.get_max_segment_num(); ++segment_num) + { + for (int axial_pos = y.get_min_axial_pos_num(segment_num); axial_pos <= y.get_max_axial_pos_num(segment_num); ++axial_pos) + { + + const Sinogram y_sinogram = y.get_sinogram(axial_pos, segment_num); + const Sinogram Ax_sinogram = Ax.get_sinogram(axial_pos, segment_num); + + for ( int view_num = y.get_min_view_num(); view_num <= y.get_max_view_num();view_num++) + { + for ( int tang_pos = y.get_min_tangential_pos_num(); tang_pos <= y.get_max_tangential_pos_num(); tang_pos++) + { + cdot1 += Ax_sinogram[view_num][tang_pos]*y_sinogram[view_num][tang_pos]; + } + } + } + } + + for ( int segment_num = x.get_min_segment_num(); segment_num <= x.get_max_segment_num(); ++segment_num) + { + for (int axial_pos = x.get_min_axial_pos_num(segment_num); axial_pos <= x.get_max_axial_pos_num(segment_num); ++axial_pos) + { + + const Sinogram x_sinogram = x.get_sinogram(axial_pos, segment_num); + const Sinogram Aty_sinogram = Aty.get_sinogram(axial_pos, segment_num); + + for ( int view_num = x.get_min_view_num(); view_num <= x.get_max_view_num();view_num++) + { + for ( int tang_pos = x.get_min_tangential_pos_num(); tang_pos <= x.get_max_tangential_pos_num(); tang_pos++) + { + cdot2 += x_sinogram[view_num][tang_pos]*Aty_sinogram[view_num][tang_pos]; + } + } + } + } + + std::cout << cdot1 << "=" << cdot2 << '\n'; + set_tolerance(0.006); + check_if_equal(cdot1, cdot2, "test adjoint"); + + +} + +END_NAMESPACE_STIR + + +USING_NAMESPACE_STIR + +int main() +{ + UpsampleDownsampleTests tests; + tests.run_tests(); + return tests.main_return_value(); +} diff --git a/src/test/test_zoom_image_adjoint.cxx b/src/test/test_zoom_image_adjoint.cxx new file mode 100644 index 0000000000..9f8911ddde --- /dev/null +++ b/src/test/test_zoom_image_adjoint.cxx @@ -0,0 +1,206 @@ +// +// +/* + Copyright (C) 2006- 2007, Hammersmith Imanet Ltd + This file is part of STIR. + This file is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + This file is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup test + \brief Test program for stir::zoom_image (and stir::centre_of_gravity) + \author Kris Thielemans +*/ + +#include "stir/VoxelsOnCartesianGrid.h" +#include "stir/IndexRange.h" +#include "stir/zoom.h" +#include "stir/centre_of_gravity.h" + +#include +#include +#include "stir/RunTests.h" + +#include + +START_NAMESPACE_STIR + +/*! + \ingroup test + \brief Test class for zoom_image (and centre_of_gravity) + The tests check if a point source remains in the same physical location + after zooming. This is done by checking the centre of gravity of the + zoomed image. +*/ +class zoom_imageAdjointTests : public RunTests +{ +public: + void run_tests(); + void fill_image_with_random(VoxelsOnCartesianGrid & image); + float dot_product(VoxelsOnCartesianGrid & image1, VoxelsOnCartesianGrid & image2); +}; + +void +zoom_imageAdjointTests::fill_image_with_random(VoxelsOnCartesianGrid & image) +{ + std::random_device random_device; + std::mt19937 random_number_generator(random_device()); + std::uniform_real_distribution number_distribution(0,10); + + for(int i=0 ; i & image1,VoxelsOnCartesianGrid & image2) +{ + float cdot = 0; + for(int i=0 ; i origin (0,1,2); + CartesianCoordinate3D grid_spacing (3,4,5); + + IndexRange<3> + range(CartesianCoordinate3D(0,-15,-14), + CartesianCoordinate3D(4,14,15)); + + VoxelsOnCartesianGrid image(range,origin, grid_spacing); + image.fill(0.F); + + const BasicCoordinate<3,int> indices = make_coordinate(1,2,3); + image[indices] = 1.F; + const CartesianCoordinate3D coord = + image.get_physical_coordinates_for_indices(indices); + + { + // check if centre_of_gravity_in_mm returns same point + check_if_equal(coord, + find_centre_of_gravity_in_mm(image), + "test on get_physical_coordinates_for_indices and find_centre_of_gravity_in_mm"); + } + + // we cannot have very good accuracy in the centre of gravity + // calculation for zooming a single pixel + // the following threshold seems very reasonable (.3mm distance) and works. + const double tolerance_for_distance = .3; + const double old_tolerance = this->get_tolerance(); + + // test 2 arg zoom_image + { + CartesianCoordinate3D new_origin (4.F,5.F,6.F); + CartesianCoordinate3D new_grid_spacing (2.2F,3.1F,4.3F); + + IndexRange<3> + new_range(CartesianCoordinate3D(-1,-16,-17), + CartesianCoordinate3D(5,15,20)); + + VoxelsOnCartesianGrid new_image(new_range,new_origin, new_grid_spacing); + + VoxelsOnCartesianGrid A(new_range,new_origin, new_grid_spacing); + VoxelsOnCartesianGrid At(range,origin, grid_spacing); + + zoom_image(new_image, image); + { + // check if centre_of_gravity_in_mm returns same point + this->set_tolerance(tolerance_for_distance); + check_if_equal(coord, + find_centre_of_gravity_in_mm(new_image), + "test on 2-argument zoom_image"); + this->set_tolerance(old_tolerance); + check_if_equal(new_range, new_image.get_index_range(), + "test on 2-argument argument zoom_image: index range"); + check_if_equal(new_grid_spacing, new_image.get_voxel_size(), + "test on 2-argument argument zoom_image: voxel size"); + check_if_equal(new_origin, new_image.get_origin(), + "test on 2-argument argument zoom_image: origin"); + + } + + + fill_image_with_random(image); + fill_image_with_random(new_image); + + //test preserve projections + zoom_image(A, image,ZoomOptions::preserve_projections); + transpose_zoom_image(At, new_image,ZoomOptions::preserve_projections); + float cdot1 = dot_product(A,new_image); + float cdot2 = dot_product(At,image); + + set_tolerance(0.004); + check_if_equal(cdot1,cdot2,"test on zoom option : preserve_projections"); + + //test preserve values + zoom_image(A, image,ZoomOptions::preserve_values); + transpose_zoom_image(At, new_image,ZoomOptions::preserve_values); + cdot1 = dot_product(A,new_image); + cdot2 = dot_product(At,image); + + set_tolerance(0.004); + check_if_equal(cdot1,cdot2,"test on zoom option : preserve_values"); + + //test preserve sum + zoom_image(A, image,ZoomOptions::preserve_sum); + transpose_zoom_image(At, new_image,ZoomOptions::preserve_sum); + cdot1 = dot_product(A,new_image); + cdot2 = dot_product(At,image); + + set_tolerance(0.004); + check_if_equal(cdot1,cdot2,"test on zoom option : preserve_sum"); + + + + + } + + + + +} + +END_NAMESPACE_STIR + + +USING_NAMESPACE_STIR + + +int main() +{ + zoom_imageAdjointTests tests; + tests.run_tests(); + return tests.main_return_value(); +}