diff --git a/CMakeLists.txt b/CMakeLists.txt index 2406a81b..747826ac 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -273,6 +273,7 @@ set(XRF_IO_HEADERS src/io/file/netcdf_io.h src/io/file/csv_io.h src/io/file/aps/aps_fit_params_import.h + src/io/file/aps/aps_roi.h src/io/file/file_scan.h src/io/file/hl_file_io.h src/io/net/basic_serializer.h @@ -299,6 +300,7 @@ set(XRF_IO_SOURCE src/io/file/netcdf_io.cpp src/io/file/file_scan.cpp src/io/file/hl_file_io.cpp + src/io/file/aps/aps_roi.cpp src/io/net/basic_serializer.cpp src/workflow/xrf/spectra_file_source.cpp src/workflow/xrf/spectra_net_source.cpp diff --git a/docs/README.md b/docs/README.md index 2d670e3a..d9536a13 100644 --- a/docs/README.md +++ b/docs/README.md @@ -10,11 +10,11 @@ #### Description ROI fitting -### NNLS +### Non-Negative Least Squares (NNLS) -### SVD +### Singualr Value Decomposition (SVD) -### MATRIX +### Per Pixel (MATRIX) ## Quantification ### Option diff --git a/src/core/main.cpp b/src/core/main.cpp index 73158e41..e3de6839 100644 --- a/src/core/main.cpp +++ b/src/core/main.cpp @@ -74,6 +74,7 @@ void help() " 1 = matrix batch fit\n 2 = batch fit without tails\n 3 = batch fit with tails\n 4 = batch fit with free E, everything else fixed \n"; logit_s<<"--optimize-fit-routine : General (default): passes elements amplitudes as fit parameters. Hybrid only passes fit parameters and fits element amplitudes using NNLS\n"; logit_s<<"--optimizer : Choose which optimizer to use for --optimize-fit-override-params or matrix fit routine \n"; + logit_s<<"--optimize-rois : Looks in 'rois' directory and performs --optimize-fit-override-params on each roi separately. \n"; logit_s<<"Fitting Routines: \n"; logit_s<< "--fit comma seperated \n"; logit_s<<" roi : element energy region of interest \n"; @@ -108,6 +109,15 @@ void set_optimizer(Command_Line_Parser& clp, data_struct::Analysis_Job& { analysis_job.set_optimizer(clp.get_option("--optimizer")); } + + if (clp.option_exists("--optimize-fit-routine")) + { + std::string opt = clp.get_option("--optimize-fit-routine"); + if (opt == "hybrid") + { + analysis_job.optimize_fit_routine = OPTIMIZE_FIT_ROUTINE::HYBRID; + } + } } // ---------------------------------------------------------------------------- @@ -212,9 +222,13 @@ template void set_fit_routines(Command_Line_Parser& clp, data_struct::Analysis_Job& analysis_job) { //Look for which analysis types we want to run - if (clp.option_exists("--fit")) + if (clp.option_exists("--fit") || clp.option_exists("--quantify-fit")) { std::string fitting = clp.get_option("--fit"); + if (fitting.length() < 1) + { + fitting = clp.get_option("--quantify-fit"); + } if (fitting.find(',') != std::string::npos) { // if we found a comma, split the string to get list of dataset files @@ -438,15 +452,6 @@ int run_optimization(Command_Line_Parser& clp) } } - if (clp.option_exists("--optimize-fit-routine")) - { - std::string opt = clp.get_option("--optimize-fit-routine"); - if (opt == "hybrid") - { - analysis_job.optimize_fit_routine = OPTIMIZE_FIT_ROUTINE::HYBRID; - } - } - if (io::file::init_analysis_job_detectors(&analysis_job)) { io::file::File_Scan::inst()->populate_netcdf_hdf5_files(analysis_job.dataset_directory); @@ -477,16 +482,16 @@ int run_quantification(Command_Line_Parser& clp) return -1; } set_fit_routines(clp, analysis_job); - + if (analysis_job.fitting_routines.size() == 0) + { + logE << "Please specify fit routines with --quantify-fit roi,nnls,matrix \n"; + return -1; + } //Check if we want to quantifiy with a standard if (clp.option_exists("--quantify-with")) { analysis_job.quantification_standard_filename = clp.get_option("--quantify-with"); } - if (clp.option_exists("--quantify-fit")) - { - //TODO: parse save as --fit so we can generate calibration curve without refitting - } if (io::file::init_analysis_job_detectors(&analysis_job)) { @@ -515,6 +520,22 @@ int run_quantification(Command_Line_Parser& clp) // ---------------------------------------------------------------------------- +int run_optimize_rois(Command_Line_Parser& clp) +{ + data_struct::Analysis_Job analysis_job; + + if (set_general_options(clp, analysis_job) == -1) + { + return -1; + } + + optimize_rois(analysis_job); + + return 0; +} + +// ---------------------------------------------------------------------------- + int run_streaming(Command_Line_Parser& clp) { data_struct::Analysis_Job analysis_job; @@ -713,7 +734,12 @@ int run_h5_file_updates(Command_Line_Parser& clp) int main(int argc, char* argv[]) { - //std::string whole_command_line = ""; + std::string whole_command_line = ""; + for (int i = 0; i < argc; i++) + { + whole_command_line += std::string(argv[i]) + " "; + } + logI << whole_command_line << "\n"; //Performance measure std::chrono::time_point start, end; @@ -770,6 +796,11 @@ int main(int argc, char* argv[]) run_quantification(clp); } + if (clp.option_exists("--optimize-rois")) + { + run_optimize_rois(clp); + } + run_h5_file_updates(clp); end = std::chrono::system_clock::now(); diff --git a/src/core/process_whole.cpp b/src/core/process_whole.cpp index ee07ce17..e39e8705 100644 --- a/src/core/process_whole.cpp +++ b/src/core/process_whole.cpp @@ -53,23 +53,17 @@ using namespace std::placeholders; //for _1, _2, // ---------------------------------------------------------------------------- bool optimize_integrated_fit_params(data_struct::Analysis_Job * analysis_job, - std::string dataset_filename, + data_struct::Spectra& int_spectra, size_t detector_num, data_struct::Params_Override* params_override, + std::string save_filename, data_struct::Fit_Parameters& out_fitp) { fitting::models::Gaussian_Model model; bool ret_val = false; - data_struct::Spectra int_spectra; if (params_override != nullptr) { - //load the quantification standard dataset - if (false == io::file::load_and_integrate_spectra_volume(analysis_job->dataset_directory, dataset_filename, detector_num, &int_spectra, params_override)) - { - logE << "In optimize_integrated_dataset loading dataset" << dataset_filename << " for detector" << detector_num << "\n"; - return false; - } //Range of energy in spectra to fit fitting::models::Range energy_range = data_struct::get_energy_range(int_spectra.size(), &(params_override->fit_params)); @@ -124,7 +118,7 @@ bool optimize_integrated_fit_params(data_struct::Analysis_Job * analysis ret_val = false; break; } - io::file::save_optimized_fit_params(analysis_job->dataset_directory, dataset_filename, detector_num, result, &out_fitp, &int_spectra, &(params_override->elements_to_fit)); + io::file::save_optimized_fit_params(analysis_job->dataset_directory, save_filename, detector_num, result, &out_fitp, &int_spectra, &(params_override->elements_to_fit)); delete fit_routine; } @@ -141,6 +135,7 @@ void generate_optimal_params(data_struct::Analysis_Job* analysis_job) std::unordered_map*> params; std::unordered_map detector_file_cnt; data_struct::Params_Override* params_override = nullptr; + data_struct::Spectra int_spectra; std::string full_path = analysis_job->dataset_directory + DIR_END_CHAR + "maps_fit_parameters_override.txt"; @@ -179,19 +174,27 @@ void generate_optimal_params(data_struct::Analysis_Job* analysis_job) } - data_struct::Fit_Parameters out_fitp; - if (optimize_integrated_fit_params(analysis_job, itr, detector_num, params_override, out_fitp)) + //load the int spectra from the dataset. + if (io::file::load_and_integrate_spectra_volume(analysis_job->dataset_directory, itr, detector_num, &int_spectra, params_override)) { - detector_file_cnt[detector_num] += 1.0; - if (fit_params_avgs.count(detector_num) > 0) + data_struct::Fit_Parameters out_fitp; + if (optimize_integrated_fit_params(analysis_job, int_spectra, detector_num, params_override, itr, out_fitp)) { - fit_params_avgs[detector_num].sum_values(out_fitp); - } - else - { - fit_params_avgs[detector_num] = out_fitp; + detector_file_cnt[detector_num] += 1.0; + if (fit_params_avgs.count(detector_num) > 0) + { + fit_params_avgs[detector_num].sum_values(out_fitp); + } + else + { + fit_params_avgs[detector_num] = out_fitp; + } } } + else + { + logE << "In optimize_integrated_dataset loading dataset" << itr << " for detector" << detector_num << "\n"; + } } } @@ -515,6 +518,100 @@ bool perform_quantification(data_struct::Analysis_Job* analysis_job) } +// ---------------------------------------------------------------------------- + +void find_and_optimize_roi(data_struct::Analysis_Job& analysis_job, int detector_num, std::map>>& rois, std::string search_filename) +{ + std::vector files = io::file::File_Scan::inst()->find_all_dataset_files(analysis_job.dataset_directory + "img.dat", search_filename); + if (files.size() > 0) + { + for (auto& roi_itr : rois) + { + Spectra int_spectra; + std::string file_path = analysis_job.dataset_directory + "img.dat" + DIR_END_CHAR + files[0]; + if (io::file::HDF5_IO::inst()->load_integrated_spectra_analyzed_h5_roi(file_path, &int_spectra, roi_itr.second)) + { + data_struct::Params_Override params_override; + //load override parameters + if (false == io::file::load_override_params(analysis_job.dataset_directory, detector_num, ¶ms_override)) + { + if (false == io::file::load_override_params(analysis_job.dataset_directory, -1, ¶ms_override)) + { + logE << "Loading maps_fit_parameters_override.txt\n"; + return; + } + } + data_struct::Fit_Parameters out_fitp; + std::string roi_name = std::to_string(roi_itr.first); + if (false == optimize_integrated_fit_params(&analysis_job, int_spectra, detector_num, ¶ms_override, files[0]+ "_roi_" + roi_name, out_fitp)) + { + logE << "Failed to optimize ROI "<< file_path<<" : "<< roi_name<<".\n"; + } + } + } + } + else + { + logE << "Could not find dataset file for " << search_filename << "\n"; + } +} + +// ---------------------------------------------------------------------------- + +void optimize_single_roi(data_struct::Analysis_Job& analysis_job, std::string roi_file_name) +{ + //std::map>> rois; + std::map>> rois; + std::string search_filename; + //data_struct::Detector* detector; + if (io::file::aps::load_v9_rois(analysis_job.dataset_directory + "rois" + DIR_END_CHAR + roi_file_name, rois)) + { + logI << "Loaded roi file "< 0) + { + // iterate through all + for (size_t detector_num : analysis_job.detector_num_arr) + { + //detector = analysis_job.get_detector(detector_num); + if (detector_num != -1) + { + std::string str_detector_num = std::to_string(detector_num); + search_filename = dataset_num + ".h5" + str_detector_num; + find_and_optimize_roi(analysis_job, detector_num, rois, search_filename); + } + } + //now for the avg h5 + //detector = analysis_job.get_detector(-1); + search_filename = dataset_num + ".h5"; + dataset_num + ".h5"; + find_and_optimize_roi(analysis_job, -1, rois, search_filename); + } + else + { + logE<<"Error parsing dataset number for "<< roi_file_name << "\n"; + } + + } + else + { + logE << "Error loading roi file "<< roi_file_name << "\n"; + } +} + +// ---------------------------------------------------------------------------- + +void optimize_rois(data_struct::Analysis_Job& analysis_job) +{ + for (auto& itr : io::file::File_Scan::inst()->find_all_dataset_files(analysis_job.dataset_directory + "rois", ".roi")) + { + optimize_single_roi(analysis_job, itr); + } +} + // ---------------------------------------------------------------------------- // ---------------------------------------------------------------------------- // ---------------------------------------------------------------------------- diff --git a/src/core/process_whole.h b/src/core/process_whole.h index 66175603..3b6d7fd8 100644 --- a/src/core/process_whole.h +++ b/src/core/process_whole.h @@ -99,6 +99,8 @@ POSSIBILITY OF SUCH DAMAGE. #include "data_struct/quantification_standard.h" +#include "io/file/aps/aps_roi.h" + using namespace std::placeholders; //for _1, _2, @@ -117,6 +119,10 @@ DLL_EXPORT void generate_optimal_params(data_struct::Analysis_Job* analy void load_and_fit_quatification_datasets(data_struct::Analysis_Job* analysis_job, size_t detector_num, vector>& standard_element_weights, unordered_map& quant_map); +void optimize_single_roi(data_struct::Analysis_Job& analysis_job, std::string roi_file_name); + +DLL_EXPORT void optimize_rois(data_struct::Analysis_Job& analysis_job); + // ---------------------------------------------------------------------------- template @@ -639,6 +645,7 @@ DLL_EXPORT void iterate_datasets_and_update(data_struct::Analysis_Job& a } } } + // ---------------------------------------------------------------------------- #endif diff --git a/src/io/file/aps/aps_fit_params_import.h b/src/io/file/aps/aps_fit_params_import.h index c984d523..b1ad6bfb 100644 --- a/src/io/file/aps/aps_fit_params_import.h +++ b/src/io/file/aps/aps_fit_params_import.h @@ -202,8 +202,8 @@ DLL_EXPORT bool load_parameters_override(std::string path, Params_Override sorted_params; + std::map sorted_elements; + + for (auto itr = fit_params.begin(); itr != fit_params.end(); itr++) + { + std::string::size_type pos = itr->second.name.find('_'); + if (pos == std::string::npos) + { + sorted_elements[itr->second.name] = itr->second.value; + } + else if (pos < 3) + { + if (itr->second.name.find("F_") == 0 || itr->second.name.find("KB_F_") == 0) + { + sorted_params[itr->second.name] = itr->second.value; + } + else + { + sorted_elements[itr->second.name] = itr->second.value; + } + } + else + { + sorted_params[itr->second.name] = itr->second.value; + } + } + if (out_stream.is_open()) { out_stream << "Fitting_Result," << result << "\n"; - for (auto itr = fit_params.begin(); itr != fit_params.end(); itr++) + for (auto itr : sorted_params) { - out_stream << itr->second.name << "," << itr->second.value << "\n"; + out_stream << itr.first << "," << itr.second << "\n"; + } + + //out_stream << "\n"; + + for (auto itr : sorted_elements) + { + out_stream << itr.first << "," << pow((T_real)10.0, itr.second) << "\n"; } out_stream.close(); return true; - } logE << "Couldn't opening file " << path << "\n"; diff --git a/src/io/file/aps/aps_roi.cpp b/src/io/file/aps/aps_roi.cpp new file mode 100644 index 00000000..46ad7b57 --- /dev/null +++ b/src/io/file/aps/aps_roi.cpp @@ -0,0 +1,143 @@ +/*** +Copyright (c) 2016, UChicago Argonne, LLC. All rights reserved. + +Copyright 2016. UChicago Argonne, LLC. This software was produced +under U.S. Government contract DE-AC02-06CH11357 for Argonne National +Laboratory (ANL), which is operated by UChicago Argonne, LLC for the +U.S. Department of Energy. The U.S. Government has rights to use, +reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR +UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR +ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is +modified to produce derivative works, such modified software should +be clearly marked, so as not to confuse it with the version available +from ANL. + +Additionally, redistribution and use in source and binary forms, with +or without modification, are permitted provided that the following +conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the + distribution. + + * Neither the name of UChicago Argonne, LLC, Argonne National + Laboratory, ANL, the U.S. Government, nor the names of its + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago +Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. +***/ + +/// Initial Author <2022>: Arthur Glowacki + + + +#include "aps_roi.h" +#include + +namespace io +{ +namespace file +{ +namespace aps +{ + +union Num +{ + char buffer[4]; + unsigned int num; +} num; + +void swapChars(char* pChar1, char* pChar2) +{ + char temp = *pChar1; + *pChar1 = *pChar2; + *pChar2 = temp; +} + +int swapOrder(Num num) +{ + swapChars(&num.buffer[0], &num.buffer[3]); + swapChars(&num.buffer[1], &num.buffer[2]); + + return num.num; +} + +bool load_v9_rois(std::string path, std::map>& rois) +{ + std::ifstream fileStream(path, std::ios::in | std::ios::binary); + + logI << "Loading: " << path << "\n"; + Num val; + unsigned int width; + unsigned int height; + unsigned int mask; + std::vector myData; + if (fileStream.is_open()) + { + fileStream.read((char*)&val, sizeof(val)); + width = swapOrder(val); + fileStream.read((char*)&val, sizeof(val)); + height = swapOrder(val); + + myData.resize(width * height); + fileStream.read((char*)myData.data(), (width * height) * sizeof(unsigned int)); + fileStream.close(); + + unsigned long i = 0; + for (unsigned int y = 0; y < height; y++) + { + for (unsigned int x = 0; x < width; x++) + { + if (myData[i]) + { + val.num = myData[i]; + mask = swapOrder(val); + unsigned int m = 1; + for (int idx = 1; idx < 12; idx++) + { + if ((m & mask) == m) + { + rois[idx-1].push_back(int_point(x, y)); + } + m = m << 1; + } + + } + i++; + } + } + fileStream.close(); + } + else + { + return false; + } + + + return true; + +} + + + + + +} //end namespace aps +} //end namespace file +}// end namespace io diff --git a/src/io/file/aps/aps_roi.h b/src/io/file/aps/aps_roi.h new file mode 100644 index 00000000..d886390f --- /dev/null +++ b/src/io/file/aps/aps_roi.h @@ -0,0 +1,73 @@ +/*** +Copyright (c) 2016, UChicago Argonne, LLC. All rights reserved. + +Copyright 2016. UChicago Argonne, LLC. This software was produced +under U.S. Government contract DE-AC02-06CH11357 for Argonne National +Laboratory (ANL), which is operated by UChicago Argonne, LLC for the +U.S. Department of Energy. The U.S. Government has rights to use, +reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR +UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR +ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is +modified to produce derivative works, such modified software should +be clearly marked, so as not to confuse it with the version available +from ANL. + +Additionally, redistribution and use in source and binary forms, with +or without modification, are permitted provided that the following +conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the + distribution. + + * Neither the name of UChicago Argonne, LLC, Argonne National + Laboratory, ANL, the U.S. Government, nor the names of its + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago +Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. +***/ + +/// Initial Author <2022>: Arthur Glowacki + + + +#ifndef _APS_ROI_ +#define _APS_ROI_ + +#include "core/defines.h" +#include +#include +#include + +namespace io +{ +namespace file +{ +namespace aps +{ + +typedef std::pair int_point; + +DLL_EXPORT bool load_v9_rois(std::string path, std::map> &rois); + +}// end namespace aps +}// end namespace file +}// end namespace io + +#endif // APS_ROI diff --git a/src/io/file/csv_io.h b/src/io/file/csv_io.h index 074ade3f..0b24bca8 100644 --- a/src/io/file/csv_io.h +++ b/src/io/file/csv_io.h @@ -356,9 +356,9 @@ namespace csv // ---------------------------------------------------------------------------- template - DLL_EXPORT bool save_fit_and_int_spectra(const std::string fullpath, const data_struct::ArrayTr* energy, const data_struct::ArrayTr* spectra, const data_struct::ArrayTr* spectra_model, const data_struct::ArrayTr* background) + DLL_EXPORT bool save_fit_and_int_spectra(const std::string fullpath, const data_struct::ArrayTr* energy, const data_struct::ArrayTr* spectra, const data_struct::ArrayTr* spectra_model, const data_struct::ArrayTr* background, unordered_map>* labeled_spectras = nullptr) { - return save_fit_and_int_spectra_labeled(fullpath, energy, spectra, spectra_model, background, nullptr); + return save_fit_and_int_spectra_labeled(fullpath, energy, spectra, spectra_model, background, labeled_spectras); } // ---------------------------------------------------------------------------- diff --git a/src/io/file/hdf5_io.h b/src/io/file/hdf5_io.h index bd7441b4..e4c5fac9 100644 --- a/src/io/file/hdf5_io.h +++ b/src/io/file/hdf5_io.h @@ -54,6 +54,7 @@ POSSIBILITY OF SUCH DAMAGE. #include #include #include +#include #include #include #include "hdf5.h" @@ -85,6 +86,8 @@ enum H5_SPECTRA_LAYOUTS {MAPS_RAW, MAPS_V9, MAPS_V10, XSPRESS, APS_SEC20}; enum GSE_CARS_SAVE_VER {UNKNOWN, XRFMAP, XRMMAP}; +using ROI_Vec = std::vector>; + template int parse_str_val_to_int(std::string start_delim, std::string end_delim, std::string lookup_str) { @@ -2802,20 +2805,9 @@ class DLL_EXPORT HDF5_IO //----------------------------------------------------------------------------- template - bool load_integrated_spectra_analyzed_h5(std::string path, data_struct::Spectra* spectra, bool log_error=true) + bool load_integrated_spectra_analyzed_h5(std::string path, data_struct::Spectra* spectra, ROI_Vec* roi = nullptr, bool log_error=true) { std::lock_guard lock(_mutex); - - //_is_loaded = ERROR_LOADING; - //std::chrono::time_point start, end; - //start = std::chrono::system_clock::now(); - /* - - if (detector_num > -1) - { - path += ".h5" + std::to_string(detector_num); - } - */ hid_t file_id; file_id = H5Fopen(path.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); @@ -2833,6 +2825,259 @@ class DLL_EXPORT HDF5_IO //----------------------------------------------------------------------------- + template + bool load_integrated_spectra_analyzed_h5_roi(std::string path, data_struct::Spectra* int_spectra, ROI_Vec& roi) + { + std::lock_guard lock(_mutex); + + bool is_v9 = false; + //_is_loaded = ERROR_LOADING; + std::chrono::time_point start, end; + start = std::chrono::system_clock::now(); + + std::stack > close_map; + + logI << path << "\n"; + + hid_t file_id, dset_id, dataspace_id, spec_grp_id, memoryspace_id, memoryspace_meta_id, dset_incnt_id, dset_outcnt_id; + hid_t memoryspace_1; + hid_t dset_rt_id, dset_lt_id, dset_scalers, dset_scaler_names; + hid_t dataspace_lt_id, dataspace_rt_id, dataspace_inct_id, dataspace_outct_id, dataspace_scalers, dataspace_scaler_names; + herr_t error; + hsize_t dims_in[3] = { 0,0,0 }; + hsize_t offset[3] = { 0,0,0 }; + hsize_t count[3] = { 1,1,1 }; + hsize_t offset_time[2] = { 0,0 }; + hsize_t count_time[2] = { 1,1 }; + hsize_t offset_1[1] = { 0 }; + hsize_t count_1[1] = { 1 }; + + hsize_t elt_off = -1, ert_off = -1, in_off = -1, out_off = -1; + + if (false == _open_h5_object(file_id, H5O_FILE, close_map, path, -1, false)) + return false; + + if (false == _open_h5_object(spec_grp_id, H5O_GROUP, close_map, "/MAPS/Spectra", file_id, false, false)) + { + if (false == _open_h5_object(spec_grp_id, H5O_GROUP, close_map, "/MAPS", file_id)) + { + return false; + } + } + + if (false == _open_h5_object(dset_id, H5O_DATASET, close_map, "mca_arr", spec_grp_id)) + { + return false; + } + + dataspace_id = H5Dget_space(dset_id); + close_map.push({ dataspace_id, H5O_DATASPACE }); + + if (false == _open_h5_object(dset_lt_id, H5O_DATASET, close_map, "Elapsed_Livetime", spec_grp_id, false, false)) + { + if (false == _open_h5_object(dset_scalers, H5O_DATASET, close_map, "scalers", spec_grp_id)) + { + return false; + } + if (false == _open_h5_object(dset_scaler_names, H5O_DATASET, close_map, "scaler_names", spec_grp_id)) + { + return false; + } + dataspace_scalers = H5Dget_space(dset_scalers); + close_map.push({ dataspace_scalers, H5O_DATASPACE }); + + dataspace_scaler_names = H5Dget_space(dset_scaler_names); + close_map.push({ dataspace_scaler_names, H5O_DATASPACE }); + hid_t memtype = H5Dget_type(dset_scaler_names); + + // read scaler names and search for elt1, ert1, incnt1, outcnt1 + hsize_t dims_out[1]; + unsigned int status_n = H5Sget_simple_extent_dims(dataspace_scaler_names, &dims_out[0], nullptr); + char tmp_name[256] = { 0 }; + memoryspace_1 = H5Screate_simple(1, count, nullptr); + for (hsize_t idx = 0; idx < dims_out[0]; idx++) + { + offset_1[0] = idx; + memset(&tmp_name[0], 0, 254); + H5Sselect_hyperslab(dataspace_scaler_names, H5S_SELECT_SET, offset_1, nullptr, count_1, nullptr); + error = H5Dread(dset_scaler_names, memtype, memoryspace_1, dataspace_scaler_names, H5P_DEFAULT, (void*)&tmp_name[0]); + if (error == 0) + { + std::string name = std::string(tmp_name); + + if (name == STR_ELT + "1") + { + elt_off = idx; + } + else if (name == STR_ERT + "1") + { + ert_off = idx; + } + else if (name == STR_ICR + "1") + { + in_off = idx; + } + else if (name == STR_OCR + "1") + { + out_off = idx; + } + } + } + is_v9 = true; + H5Tclose(memtype); + } + else + { + dataspace_lt_id = H5Dget_space(dset_lt_id); + close_map.push({ dataspace_lt_id, H5O_DATASPACE }); + + if (false == _open_h5_object(dset_rt_id, H5O_DATASET, close_map, "Elapsed_Realtime", spec_grp_id)) + return false; + dataspace_rt_id = H5Dget_space(dset_rt_id); + close_map.push({ dataspace_rt_id, H5O_DATASPACE }); + + if (false == _open_h5_object(dset_incnt_id, H5O_DATASET, close_map, "Input_Counts", spec_grp_id)) + return false; + dataspace_inct_id = H5Dget_space(dset_incnt_id); + close_map.push({ dataspace_inct_id, H5O_DATASPACE }); + + if (false == _open_h5_object(dset_outcnt_id, H5O_DATASET, close_map, "Output_Counts", spec_grp_id)) + return false; + dataspace_outct_id = H5Dget_space(dset_outcnt_id); + close_map.push({ dataspace_outct_id, H5O_DATASPACE }); + } + + int rank = H5Sget_simple_extent_ndims(dataspace_id); + if (rank != 3) + { + _close_h5_objects(close_map); + logE << "Dataset /MAPS/Spectra/mca_arr rank != 3. rank = " << rank << ". Can't load dataset. returning" << "\n"; + return false; + //throw exception ("Dataset is not a volume"); + } + + int status_n = H5Sget_simple_extent_dims(dataspace_id, &dims_in[0], nullptr); + if (status_n < 0) + { + _close_h5_objects(close_map); + logE << "getting dataset dims for /MAPS/Spectra/mca_arr" << "\n"; + return false; + } + + for (int i = 0; i < rank; i++) + { + + offset[i] = 0; + count[i] = dims_in[i]; + } + + + int_spectra->resize(dims_in[0]); + int_spectra->setZero(dims_in[0]); + + count[0] = dims_in[0]; + count[1] = 1; + count[2] = 1; + + memoryspace_id = H5Screate_simple(3, count, nullptr); + memoryspace_meta_id = H5Screate_simple(2, count_time, nullptr); + H5Sselect_hyperslab(memoryspace_id, H5S_SELECT_SET, offset, nullptr, count, nullptr); + H5Sselect_hyperslab(memoryspace_meta_id, H5S_SELECT_SET, offset_time, nullptr, count_time, nullptr); + + T_real live_time = 1.0; + T_real real_time = 1.0; + T_real in_cnt = 1.0; + T_real out_cnt = 1.0; + + + for (auto& itr : roi) + { + hsize_t xoffset = itr.first; + hsize_t yoffset = itr.second; + + offset[0] = 0; + offset[1] = yoffset; + offset_time[0] = yoffset; + offset[2] = xoffset; + offset_time[1] = xoffset; + + count[0] = dims_in[0]; + + H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, offset, nullptr, count, nullptr); + + data_struct::Spectra spectra; + spectra.resize(dims_in[0]); + + error = _read_h5d(dset_id, memoryspace_id, dataspace_id, H5P_DEFAULT, (void*)spectra.data()); + if (error > 0) + { + logW << "Counld not read row " << yoffset << " col " << xoffset << "\n"; + } + + if (is_v9) + { + count[0] = 1; + + offset[0] = elt_off; + H5Sselect_hyperslab(dataspace_scalers, H5S_SELECT_SET, offset, nullptr, count, nullptr); + error = _read_h5d(dset_scalers, memoryspace_meta_id, dataspace_scalers, H5P_DEFAULT, (void*)&live_time); + if (error > 0) + live_time = 0; + + offset[0] = ert_off; + H5Sselect_hyperslab(dataspace_scalers, H5S_SELECT_SET, offset, nullptr, count, nullptr); + error = _read_h5d(dset_scalers, memoryspace_meta_id, dataspace_scalers, H5P_DEFAULT, (void*)&real_time); + if (error > 0) + real_time = 0; + + offset[0] = in_off; + H5Sselect_hyperslab(dataspace_scalers, H5S_SELECT_SET, offset, nullptr, count, nullptr); + error = _read_h5d(dset_scalers, memoryspace_meta_id, dataspace_scalers, H5P_DEFAULT, (void*)&in_cnt); + if (error > 0) + in_cnt = 0; + + offset[0] = out_off; + H5Sselect_hyperslab(dataspace_scalers, H5S_SELECT_SET, offset, nullptr, count, nullptr); + error = _read_h5d(dset_scalers, memoryspace_meta_id, dataspace_scalers, H5P_DEFAULT, (void*)&out_cnt); + if (error > 0) + out_cnt = 0; + } + else + { + H5Sselect_hyperslab(dataspace_lt_id, H5S_SELECT_SET, offset_time, nullptr, count_time, nullptr); + H5Sselect_hyperslab(dataspace_rt_id, H5S_SELECT_SET, offset_time, nullptr, count_time, nullptr); + H5Sselect_hyperslab(dataspace_inct_id, H5S_SELECT_SET, offset_time, nullptr, count_time, nullptr); + H5Sselect_hyperslab(dataspace_outct_id, H5S_SELECT_SET, offset_time, nullptr, count_time, nullptr); + + error = _read_h5d(dset_rt_id, memoryspace_meta_id, dataspace_rt_id, H5P_DEFAULT, (void*)&real_time); + error = _read_h5d(dset_lt_id, memoryspace_meta_id, dataspace_lt_id, H5P_DEFAULT, (void*)&live_time); + error = _read_h5d(dset_incnt_id, memoryspace_meta_id, dataspace_inct_id, H5P_DEFAULT, (void*)&in_cnt); + error = _read_h5d(dset_outcnt_id, memoryspace_meta_id, dataspace_outct_id, H5P_DEFAULT, (void*)&out_cnt); + } + + spectra.elapsed_livetime(live_time); + spectra.elapsed_realtime(real_time); + spectra.input_counts(in_cnt); + spectra.output_counts(out_cnt); + + int_spectra->add(spectra); + } + + int_spectra->recalc_elapsed_livetime(); + + _close_h5_objects(close_map); + + end = std::chrono::system_clock::now(); + std::chrono::duration elapsed_seconds = end - start; + //std::time_t end_time = std::chrono::system_clock::to_time_t(end); + + logI << "elapsed time: " << elapsed_seconds.count() << "s\n"; + + return true; + } + + //----------------------------------------------------------------------------- + template bool load_quantification_scalers_analyzed_h5(std::string path, data_struct::Params_Override *override_values) { @@ -2856,15 +3101,30 @@ class DLL_EXPORT HDF5_IO if (false == _open_h5_object(file_id, H5O_FILE, close_map, path, -1)) return false; - if (false == _open_h5_object(ds_ic_id, H5O_DATASET, close_map, "/MAPS/Quantification/Scalers/DS_IC", file_id)) - return false; - - if (false == _open_h5_object(us_ic_id, H5O_DATASET, close_map, "/MAPS/Quantification/Scalers/US_IC", file_id)) - return false; + if (false == _open_h5_object(ds_ic_id, H5O_DATASET, close_map, "/MAPS/Quantification/Standard0/Scalers/DS_IC", file_id, false, false)) + { + if (false == _open_h5_object(ds_ic_id, H5O_DATASET, close_map, "/MAPS/Quantification/Standard1/Scalers/DS_IC", file_id)) + { + return false; + } + } + - if (false == _open_h5_object(srcurrent_id, H5O_DATASET, close_map, "/MAPS/Quantification/Scalers/SR_Current", file_id)) - return false; + if (false == _open_h5_object(us_ic_id, H5O_DATASET, close_map, "/MAPS/Quantification/Standard0/Scalers/US_IC", file_id, false, false)) + { + if (false == _open_h5_object(us_ic_id, H5O_DATASET, close_map, "/MAPS/Quantification/Standard1/Scalers/US_IC", file_id)) + { + return false; + } + } + if (false == _open_h5_object(srcurrent_id, H5O_DATASET, close_map, "/MAPS/Quantification/Standard0/Scalers/SR_Current", file_id, false, false)) + { + if (false == _open_h5_object(srcurrent_id, H5O_DATASET, close_map, "/MAPS/Quantification/Standard1/Scalers/SR_Current", file_id)) + { + return false; + } + } //read in scaler hid_t d_space = H5Dget_space(srcurrent_id); diff --git a/src/io/file/hl_file_io.cpp b/src/io/file/hl_file_io.cpp index 48ce3bb5..c0d036ad 100644 --- a/src/io/file/hl_file_io.cpp +++ b/src/io/file/hl_file_io.cpp @@ -255,7 +255,9 @@ void save_optimized_fit_params(std::string dataset_dir, std::string dataset_file fitting::models::Range energy_range = data_struct::get_energy_range(spectra->size(), fit_params); data_struct::Spectra snip_spectra = spectra->sub_spectra(energy_range.min, energy_range.count()); - data_struct::Spectra model_spectra = model.model_spectrum_mp(fit_params, elements_to_fit, energy_range); + unordered_map> labeled_spectras; + data_struct::Spectra model_spectra = model.model_spectrum(fit_params, elements_to_fit, &labeled_spectras, energy_range); + data_struct::ArrayTr background; double energy_offset = fit_params->value(STR_ENERGY_OFFSET); @@ -304,7 +306,7 @@ void save_optimized_fit_params(std::string dataset_dir, std::string dataset_file visual::SavePlotSpectrasFromConsole(str_path, &ev, &snip_spectra, &model_spectra, &background, true); #endif - io::file::csv::save_fit_and_int_spectra(full_path, &ev, &snip_spectra, &model_spectra, &background); + io::file::csv::save_fit_and_int_spectra(full_path, &ev, &snip_spectra, &model_spectra, &background, &labeled_spectras); io::file::aps::save_fit_parameters_override(fp_full_path, *fit_params, result); std::unordered_map scaler_map; scaler_map[STR_ENERGY_OFFSET] = energy_offset; diff --git a/src/io/file/hl_file_io.h b/src/io/file/hl_file_io.h index 68ce33a7..c1746163 100644 --- a/src/io/file/hl_file_io.h +++ b/src/io/file/hl_file_io.h @@ -442,7 +442,7 @@ DLL_EXPORT bool load_and_integrate_spectra_volume(std::string dataset_directory, { fullpath += std::to_string(detector_num); } - if (true == io::file::HDF5_IO::inst()->load_integrated_spectra_analyzed_h5(fullpath, integrated_spectra, false)) + if (true == io::file::HDF5_IO::inst()->load_integrated_spectra_analyzed_h5(fullpath, integrated_spectra, nullptr, false)) { logI << "Loaded integradted spectra from h5.\n"; if (params_override != nullptr) diff --git a/src/io/file/mca_io.h b/src/io/file/mca_io.h index 9b84ab2c..b393325c 100644 --- a/src/io/file/mca_io.h +++ b/src/io/file/mca_io.h @@ -144,7 +144,7 @@ DLL_EXPORT bool load_integrated_spectra(std::string path, data_struct::Spectrasize() << "\n"; - paramFileStream << "T_realIME: " << spectra->elapsed_realtime() << "\n"; + paramFileStream << "REAL_TIME: " << spectra->elapsed_realtime() << "\n"; paramFileStream << "LIVE_TIME: " << spectra->elapsed_livetime() << "\n"; paramFileStream << "CAL_OFFSET: " << offset << "\n"; paramFileStream << "CAL_SLOPE: " << slope << "\n"; diff --git a/vcpkg b/vcpkg index a106de33..203f5f3b 160000 --- a/vcpkg +++ b/vcpkg @@ -1 +1 @@ -Subproject commit a106de33bbee694e3be6243718aa2a549a692832 +Subproject commit 203f5f3bef28ffccb667fdce97bcde7a88aba4e6